Code covered by the BSD License  

Highlights from
Optimizing breakpoints for Tables

image thumbnail

Optimizing breakpoints for Tables

by

 

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')

Contact us