Code covered by the BSD License

# Optimizing breakpoints for Tables

### Richard Willey (view profile)

MATLAB code to support the "Generating Optimal Tables using MATLAB" webinar.

TableDemo.m
```% Copyright 2009 The MathWorks, Inc.

%% Introduction

% TableDemo.m contains the code used during the "Generating Optimal Tables
% using MATLAB Products" webinar.

%% Generate a reference model using the peaks function

% In the live webinar I did this using the Surface Fitting Tool
% To make life easy, here's some MATLAB code

% Use a halton set to generate some random x, y data
P = haltonset( 2 );
P = scramble( P, 'RR2' );
X = net( P, 30);

x = ((X(:,1)) * 6) - 3;
y = (X(:,2) * 6) - 3;

% Use the peaks function to generate a z vector

z = peaks(x,y);

% Use the fit command to create a fit object

ft = fittype( 'a*(1-x).^2.*exp(-(x.^2) - (y+1).^2) +  b*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) + c*exp(-(x+1).^2 - y.^2)', 'indep', {'x', 'y'}, 'depend', 'z' );
opts = fitoptions( ft );
opts.Display = 'Off';
opts.Lower = [-Inf -Inf -Inf];
opts.StartPoint = [0.698493238453756 0.502175213611344 0.687201121067495];
opts.Upper = [Inf Inf Inf];
opts.Weights = zeros(1,0);
[fittedmodel, gof] = fit( [x, y], z, ft, opts );

%% Demonstrate some capabilities of a fit object

plot(fittedmodel)
plot(fittedmodel, 'style','contour')
fittedmodel(0,0)
methods(fittedmodel)

%% Problem Definition

% A 2D lookup table is an approximation of a response with points defined
% by a row vector and column vector of containing M and N elements,
% respectively.  For example, we can approximate our model using a uniform
% 13 x 13 grid

a = -3; b = 3;
x = linspace(a,b,13);
y = linspace(a,b,13);
[X,Y] = meshgrid(x,y);

plot(fittedmodel, 'style','contour')
hold on
scatter(X(:),Y(:),'k','filled');

%% Generate a Fit Object to describe our 13 x 13 uniform table

Z = fittedmodel(X,Y);
fitObj = fit([X(:), Y(:)], Z(:), 'linearinterp')

hold off
plot(fitObj)

%% Visually compare the two surfaces

% Generate a reference grid to evaluate the two surfaces
xr = linspace(a,b,100);
yr = linspace(a,b,100);
[Xr, Yr] = meshgrid(xr,yr);

% Calculate the difference between the two surfaces
resid = fittedmodel(Xr(:), Yr(:)) - fitObj(Xr(:), Yr(:));

% Create a Fit Object describing the residuals
Diff_Contour = fit([Xr(:), Yr(:)], resid, 'linearinterp');

% Generate a contour plot of the residuals
figure('Numbertitle', 'off', 'name', 'Contour Map of Residuals:  Uniform Grid')

xlim = [-3, 3];
ylim = [-3, 3];
obj = Diff_Contour;

[xi, yi] = meshgrid( ...
linspace( xlim(1), xlim(2), 49 ), ...
linspace( ylim(1), ylim(2), 51 ) );
zi = feval( obj, xi, yi );

[~, h] = contourf( xi, yi, zi, 21 );

grid on
colorbar
caxis ([-1.3 1.3])

% Fill in the break points
hold on
scatter(X(:), Y(:), 'k', 'filled')

%% Objective Function Definition
% We want to minimize the residuals, which can be done using a least
% squares approach, to identify the best x and y locations.
type objFcn
%%
% Evaluate it
M = length(x)-2; % remove endpoints
N = length(y)-2;
xrange = [min(x) max(x)];
yrange = [min(y) max(y)];
deltax = diff(x);
deltay = diff(y);
delta  = [deltax(1:end-1) deltay(1:end-1)];
sqError = @(delta) objFcn(delta,M,xrange,yrange,fittedmodel,Xr,Yr);
sqError(delta)

% Alternatively, we can use the "norm" command
% Use the "norm" command to create a summary statistics that describes the
% goodness of fit

% norm calculates the sum of the squared errors
% norm allows us to collapse the residual matrix into a scalar

%% Define Bounds
% values of XY must be >= 0
lb = zeros(size(delta));

%% Define Constraint Definition
% equality A*x < b
A = zeros(2,M+N);
A(1,1:M) = 1;
A(2,M+1:end) = 1;
b = [x(end) - x(1);
y(end) - y(1)];
%% Set Options
plotFcn = @(x,o,s) plotCurrentFit(x,o,s,fittedmodel,xrange,yrange,M);
options = optimset('Disp','Iter','Algorithm','interior-point',...
'UseParallel','Always','OutputFcn',plotFcn,'TolX',1e-3);
%% Solve the problem

% matlabpool open speedy-optim 16

delta0 = delta;
[delta lsqerror] = fmincon(sqError,delta0,A,b,[],[],lb,[],[],options);

% matlabpool close

% New x/y
[newx,newy] = delta2xy(delta,M,xrange,yrange);

%% Visually compare the two surfaces

% Create a new z variable
[NewX,NewY] = meshgrid(newx,newy);
NewZ = fittedmodel(NewX,NewY);

% Create a new fit object

fitObj2 = fit([NewX(:), NewY(:)], NewZ(:), 'linearinterp')

% Calculate the difference between the two surfaces
resid = fittedmodel(Xr(:), Yr(:)) - fitObj2(Xr(:), Yr(:));

% Create a Fit Object describing the residuals
Diff_Contour2 = fit([Xr(:), Yr(:)], resid, 'linearinterp');

% Generate a contour plot of the residuals
figure('Numbertitle', 'off', 'name', 'Contour Map of Residuals:  Optimized Grid')

xlim = [-3, 3];
ylim = [-3, 3];
obj = Diff_Contour2;

[xi, yi] = meshgrid( ...
linspace( xlim(1), xlim(2), 49 ), ...
linspace( ylim(1), ylim(2), 51 ) );
zi = feval( obj, xi, yi );

[~, h] = contourf( xi, yi, zi, 21 );

grid on
colorbar
caxis ([-1.3 1.3])

% Fill in the break points
hold on
scatter(NewX(:), NewY(:), 'k', 'filled')

%% Generate a set of Pareto frontiers comparing the various techniques
% open Pareto.mat
%
% figure('Numbertitle', 'off', 'name', 'Pareto Chart')
% Plot1 = plot(ans.Size, ans.SplineOptimized, 'color', 'b', 'linewidth', 2);
%
% xlabel('Size')
% ylabel('Norm(Residual Matrix)')
%
% hold on
%
% Plot2 = plot(ans.Size, ans.SplineUniform, 'color', 'b', 'linestyle', '--', 'linewidth', 2);
% Plot3 = plot(ans.Size, ans.LinearOptimized, 'color', 'k', 'linewidth', 2);
%
% legend('SplineOptimized','SplineUniform','LinearOptimized','LinearUniform')
```