Code covered by the BSD License  

Highlights from
Optimal Inverse Function Creation

image thumbnail

Optimal Inverse Function Creation

by

 

Creates continuous, optimal inverse-functions: given a desired output, finds the optimal input

InvFun_Test

Contents

% Set up problem

% % % domain
x_min=[-1;-1.00]; %min values for each dimension
x_max=[ 1; 1.25]; %max values for each dimension

% % % functions, label one Y and Y_Grad, and another set J and J_Grad

% straight up parabola
J     =@(x) +((x(1)-0).^2+(x(2)-0.2).^2);
J_Grad=@(x) [2*(x(1)-0);2*(x(2)-0.2)];

% linear in x1 and quadratic in x2
% Y=@(x) 1.25*(x(1)-0.2).*(x(2).^2);
% Y_Grad=@(x) [1.25*(x(2).^2);1.25*(x(1)-0.2).*2*x(2)];

% linear in 1d and parabola in other dimension
% Y     =@(x) (x(1).^1).*(sin(2*pi/2*x(2)));
% Y_Grad=@(x) [(sin(2*pi/2*x(2)));2*pi/2*(x(1).^1).*(cos(2*pi/2*x(2)))];

% Multimodal gausian, peaks(a1*2,a2*2)/7, equation is used to show gradient
Y     =@(x) (3*(1-2*x(1))^2*exp(-4*x(1)^2-(2*x(2)+1)^2)+...
    -exp(-(2*x(1)+1)^2-4*x(2)^2)/3+...
    -10*(2*x(1)/5-8*x(1)^3-32*x(2)^5)*exp(-4*x(1)^2-4*x(2)^2))/7;
Y_Grad=@(x) [-(12+24*x(1)*(1-2*x(1)))*(1-2*x(1))*exp(-4*x(1)^2-(2*x(2)+1)^2)+...
    4*(2*x(1)+1)*exp(-(2*x(1)+1)^2-4*x(2)^2)/3+...
    10*(-2/5+24*x(1)^2+8*x(1)*(2*x(1)/5-8*x(1)^3-32*x(2)^5))*exp(-4*x(1)^2-4*x(2)^2);...
    3*(1-2*x(1))^2*(-4*(2*x(2)+1))*exp(-4*x(1)^2-(2*x(2)+1)^2)+...
    +8*x(2)*exp(-(2*x(1)+1)^2-4*x(2)^2)/3+...
    +10*(160*x(2)^4+8*x(2)*(2*x(1)/5-8*x(1)^3-32*x(2)^5))*exp(-4*x(1)^2-4*x(2)^2)]/7;

optimization problem statement and parameters

% ProbStrut.J=@(x) x(1).^2+x(2).^2; % cost function
% ProbStrut.J_Grad=@(x) [2*x(1);2*x(2)]; % gradient of the cost function
% ProbStrut.y=@(x) x(2); % output function
% ProbStrut.y_Grad=@(x) [0;1]; % gradient of the cost function
% ProbStrut.y=@(x) x(1)+x(2); % output function
% ProbStrut.y_Grad=@(x) [1;1]; % gradient of the cost function

ProbStrut.J     =J;
ProbStrut.y     =Y;
ProbStrut.J_Grad=J_Grad;
ProbStrut.y_Grad=Y_Grad;

ProbStrut.sat = @(x,beta) x/max(norm(x),beta); %saturation function used

% domain, used to size initial population
ProbStrut.x_min=x_min;
ProbStrut.x_max=x_max;

% step parameters
ProbStrut.beta_J=0.05; % saturation length for cost gradient
ProbStrut.beta_y=0.025; % saturation length for output gradient
ProbStrut.beta_Delta=0.5; % saturation length for step direction
ProbStrut.beta_delta=0.88; % step length reduction factor for Armijo condition
ProbStrut.sigma=0.2; % Armijo rule factor
ProbStrut.gamma_J=0.05; % max step length for cost reduction steps
ProbStrut.gamma_y=0.03; % max step length for output change steps

% cluster parameters
ProbStrut.knots_min=7.5e-6; %If new knot is within this, cluster is terminated
ProbStrut.knots_max=0.1;%If new knot is further than this, cluster is terminated
ProbStrut.knots_itter_max=150; %If new knot does not converge within this, cluster is terminated

% other parameters
ProbStrut.pop_init_size=60; % number of initial agents to search space
ProbStrut.Delta_min=1e-4; % if step length is less, agent is converged
ProbStrut.pop_dist_kill=0.15; % points within this distance are terminated
ProbStrut.pop_itter_max=400; % number of itterations for population steps

ProbStrut.BoundViolateOK=1e-10; %so there isn't numerical issues with bounds

set up function to watch particles move

plot function is optional argument first argument is a list of the population (columns are agents) second argument is a style for different phases: 'b.' for moving, 'r+' for cluster,'c*' for cluster start point and 'go' for agent removed due to proximity.

% give the figure a background
figure(61);clf;
[temp_xx,temp_yy]=meshgrid(linspace(ProbStrut.x_min(1),ProbStrut.x_max(1),40),...
    linspace(ProbStrut.x_min(2),ProbStrut.x_max(2),55));
YY=zeros(size(temp_xx));JJ=zeros(size(temp_xx));
for ella=1:numel(temp_xx)
    YY(ella)=ProbStrut.y([temp_xx(ella);temp_yy(ella)]);
    JJ(ella)=ProbStrut.J([temp_xx(ella);temp_yy(ella)]);
end

contour(temp_xx,temp_yy,YY);hold on; %need hold to be on to see progression
contour(temp_xx,temp_yy,JJ);
PlotStep_axe=gca; % this ensures same figure is used each time
title({'Progress of optimization';'Legend: b.: path, r+: cluster, ';...
    'c*: Initial cluster point and go for killed path'});

ProbStrut.Plot_Step=@(Pops,Symb) plot(PlotStep_axe,Pops(1,:),Pops(2,:),Symb);

ProbStrut.Plot_final=@(Pops,Symb) pause(0.1); %just to slow it down to watch

Initialize population

Pop_I=bsxfun(@plus,bsxfun(@times,rand(length(ProbStrut.x_min),ProbStrut.pop_init_size),...
    ProbStrut.x_max-ProbStrut.x_min),ProbStrut.x_min);

Call optimization

InvFuntions=InvFun(Pop_I,ProbStrut);
% InvFuntions=InvFun([0;0],ProbStrut)

Show results

Comparison of all inverse functions

temp_color=jet(length(InvFuntions));
figure(41);clf;
for carly=1:length(InvFuntions)
for dani=1:length(ProbStrut.x_min)
subplot(length(ProbStrut.x_min)+1,1,dani)
hold on;grid on;
plot(InvFuntions{carly}.y,InvFuntions{carly}.x(dani,:),'-+','color',temp_color(carly,:));
ylabel(['Input x_' num2str(dani)]);
end
subplot(length(ProbStrut.x_min)+1,1,length(ProbStrut.x_min)+1)
hold on;grid on;
plot(InvFuntions{carly}.y,InvFuntions{carly}.J,'-+','color',temp_color(carly,:));
ylabel('J, cost');
xlabel('y, output');
end

Accuracy and optimality of all inverse functions

Pnts_on=50;
Pnts_off=200;
Pnts_off_baseNoise=0.01; %additional noise for near points
Pnts_mrk=5; %marker size for dots in plots

for carly=1:length(InvFuntions)
figure(100+carly);clf;

% points on the inverse function
Pnts_on_test_y=rand(1,Pnts_on)*(InvFuntions{carly}.y(end)-InvFuntions{carly}.y(1))+InvFuntions{carly}.y(1);
Pnts_on_test_x=pchip(InvFuntions{carly}.y,InvFuntions{carly}.x,Pnts_on_test_y);
Pnts_on_test_y_act=zeros(size(Pnts_on_test_y));
Pnts_on_test_J=zeros(size(Pnts_on_test_y));
for amy=1:size(Pnts_on_test_x,2)
    Pnts_on_test_y_act(amy)=ProbStrut.y(Pnts_on_test_x(:,amy));
    Pnts_on_test_J(amy)=ProbStrut.J(Pnts_on_test_x(:,amy));
end

% nearby points
% using testing points should be more uniform in output... maybe
temp_cov=cov(Pnts_on_test_x.').'+Pnts_off_baseNoise*eye(length(ProbStrut.x_min));
Pnts_off_test_x=bsxfun(@plus,sqrtm(temp_cov)*randn(length(ProbStrut.x_min),Pnts_off),...
    mean(Pnts_on_test_x,2));
Pnts_off_test_x=bsxfun(@max,bsxfun(@min,Pnts_off_test_x,ProbStrut.x_max),ProbStrut.x_min);
Pnts_off_test_y=zeros(1,size(Pnts_off_test_x,2));
Pnts_off_test_J=zeros(1,size(Pnts_off_test_x,2));
for amy=1:size(Pnts_off_test_x,2)
    Pnts_off_test_y(amy)=ProbStrut.y(Pnts_off_test_x(:,amy));
    Pnts_off_test_J(amy)=ProbStrut.J(Pnts_off_test_x(:,amy));
end


subplot(1,2,1)
plot(Pnts_off_test_x(1,:),Pnts_off_test_x(2,:),'g.','markersize',Pnts_mrk);
hold on;
plot(InvFuntions{carly}.x(1,:),InvFuntions{carly}.x(2,:),'r+');
plot(Pnts_on_test_x(1,:),Pnts_on_test_x(2,:),'m.','markersize',Pnts_mrk);
contour(temp_xx,temp_yy,YY);
contour(temp_xx,temp_yy,JJ);
xlabel('x_1');ylabel('x_2'); title('Location of test and near points');
legend({'Near points','Inverse function','Test points'},'location','best')
axis([ProbStrut.x_min(1),ProbStrut.x_max(1),ProbStrut.x_min(2),ProbStrut.x_max(2)]);

subplot(2,2,2)
plot([InvFuntions{carly}.y(end),InvFuntions{carly}.y(1)],[0,0],'r-');hold on;
plot(Pnts_on_test_y,Pnts_on_test_y-Pnts_on_test_y_act,'m.','markersize',Pnts_mrk);
xlabel('Output, y');ylabel('Error, predicted-actual'); title('Accuracy of test points');

subplot(2,2,4)
plot(Pnts_off_test_y,Pnts_off_test_J,'g.','markersize',Pnts_mrk);hold on;
plot(InvFuntions{carly}.y,InvFuntions{carly}.J,'r+');
plot(Pnts_on_test_y,Pnts_on_test_J,'m.','markersize',Pnts_mrk);
xlabel('Output, y');ylabel('Cost'); title('Optimality of test points');
legend({'Near points','Inverse function','Test points'})
end

Contact us