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

InvFuntions=InvFun(Pop_I,ProbStrut)
function InvFuntions=InvFun(Pop_I,ProbStrut)
%INVFUN Finds inverse optimal functions
% Pop_I are the initial population (every column is an agent vector) 
% ProbStrut contains all problem parameters 

% ProbStrut, provides all parameters of optimization
% 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
% 
% % domain, used to size initial population 
% ProbStrut.x_min=[-1;-1.50]; %min values for each dimension
% ProbStrut.x_max=[ 1; 1.25]; %max values for each dimension
% 
% % 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.02; % 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-6; % 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

Pops=Pop_I;
InvFuntions=cell(0);
KillList=zeros(length(ProbStrut.x_min),0);
if ~isfield(ProbStrut,'Plot_Step'); %assign blank function otherwise
    ProbStrut.Plot_Step=@(x,style) [];
end
if ~isfield(ProbStrut,'Plot_final'); %assign blank function otherwise
    ProbStrut.Plot_final=@(x,style) [];
end



%% do each itteration 
for abby=1:ProbStrut.pop_itter_max

%% do each agent individually 
for barbie=size(Pops,2):-1:1 %reverse ordering so agents can be removed

Pnt=Pops(:,barbie);
%% step to reduce cost 
[Pnt_new, Pnt_step] = InvFun_Armijo(Pnt, ProbStrut);

%% check for settling condition

if norm(Pnt_step)<=ProbStrut.Delta_min
InvFuntions{end+1}=InvFun_Cluster(Pnt_new, ProbStrut); %#ok<AGROW> %form cluster

ProbStrut.Plot_Step(InvFuntions{end}.x,'r+-');
ProbStrut.Plot_Step(Pops(:,barbie),'c*');

Pops(:,barbie)=[];% remove agent
KillList=[KillList,InvFuntions{end}.x]; %#ok<AGROW>
else 
Pops(:,barbie)=Pnt_new;
end

end

if size(Pops,2)<1
    break
end

%% test for removal by proximity
ListAll=[Pops, KillList];
tempA=repmat(reshape(ListAll.',1,size(ListAll,2),[]),size(Pops,2),1);
tempB=repmat(reshape(Pops.',size(Pops,2),1,[]),1,size(ListAll,2));
% The relative direction of any particle with respect to any other
Direction_all=tempA-tempB;
% Find the distance, the lower diag inf makes it so only the lower
% indexed particle is labeled as violating the distance, prevents
% double deletion, though allows for chain deletion
Dist=sum(Direction_all.^2,3).^.5+...
    [tril(inf+zeros(size(Pops,2))), ...
    zeros(size(Pops,2),size(ListAll,2)-size(Pops,2))];
% [size(Dist), size(min(Dist,[],2)), size(Pops)]

ProbStrut.Plot_Step(Pops(:,min(Dist,[],2)<ProbStrut.pop_dist_kill),'go');

Pops(:,min(Dist,[],2)<ProbStrut.pop_dist_kill)=[];

if size(Pops,2)<1
    break
end

ProbStrut.Plot_Step(Pops,'b.');
ProbStrut.Plot_final(Pops,'b.')
drawnow;

end
end
%% %% %% Subfunctions %% %% %% 

%% Cost reduction step
function [Pnt_new, Pnt_step] = InvFun_Armijo(Pnt, ProbStrut)
%INVFUN_OUTSTEP Produces a step in the output direction

%% run loop around checking boundary violations
looping=true;
Voilate_dims=false(size(ProbStrut.x_min));
Pnt_new=Pnt;
while looping
Pnt_old=Pnt_new;

%% Calculate and apply step 
J_hat=ProbStrut.sat(ProbStrut.J_Grad(Pnt_old).*~Voilate_dims,...
    ProbStrut.beta_J); % equation #1 
f_hat=ProbStrut.sat(ProbStrut.y_Grad(Pnt_old).*~Voilate_dims,...
    ProbStrut.beta_y); % equation #2 
Delta_J=-J_hat+(J_hat.'*f_hat)*f_hat; % eqn #3


%% Check Armijo condition
checking_A=true;
temp_m=1;
while checking_A
temp_step=ProbStrut.gamma_J*(ProbStrut.beta_delta)^temp_m*...
    ProbStrut.sat(Delta_J,ProbStrut.beta_Delta); 
Pnt_new=Pnt_old+temp_step; %equation #4
if (ProbStrut.J(Pnt_new)-ProbStrut.J(Pnt_old))<=...
        (ProbStrut.sigma*temp_step.'*ProbStrut.J_Grad(Pnt_old)) %eqn # 5
    checking_A=false;
else
    temp_m=temp_m+1;
end
end

%% Check for domain violation
Violate=max(max(ProbStrut.x_min-Pnt_new,Pnt_new-ProbStrut.x_max),0);
[Violate_reduce,Violate_dim]=max(Violate);
if Violate_reduce<=ProbStrut.BoundViolateOK
    looping=false; %no violations, so end
else
    Voilate_dims(Violate_dim)=true;
    temp_step=temp_step*(1-abs(Violate_reduce/temp_step(Violate_dim)));
    Pnt_new=Pnt_old+temp_step;
end

end

Pnt_new=max(min(Pnt_new,ProbStrut.x_max),ProbStrut.x_min);
Pnt_step=Pnt_new-Pnt;
end

%% Clustering step 
function [ InvFuntions ] = InvFun_Cluster(Pnt, ProbStrut)
%INVFUN_CLUSTER finds the optimal inverse function at the input point
% Pnt, is the current point
% ProbStrut, provides all parameters of optimization
% InvFuntions has fields .x, .y, .J which are the ordered sets

%% set up 
ProbStrut.pop_itter_max=ProbStrut.knots_itter_max; 
% clusters are required to converge earlier
% scope of the varible is such it doesn't return to the main function
InvFuntions.x=Pnt;
InvFuntions.y=ProbStrut.y(Pnt);
InvFuntions.J=ProbStrut.J(Pnt);

%% loop in each direction 
for Dir=[-1,1];
looping = true; 

if Dir>0
Pnt_new=InvFuntions.x(:,end);
else
Pnt_new=InvFuntions.x(:,1);
end 

while looping
%% step in output direction
Pnt_old = Pnt_new;
Pnt_new = InvFun_OutStep(Pnt_new, ProbStrut, Dir);
%% loop over itterations 
failed=true;
for abby=1:ProbStrut.pop_itter_max
%% step to reduce cost 
% Pnt_new
[Pnt_new, Pnt_step] = InvFun_Armijo(Pnt_new, ProbStrut);
%% check for settling condition
if norm(Pnt_step)<=ProbStrut.Delta_min
    failed=false;
    break;
end
end %ending optimization of new end point 
%% clasify new point as termination, top or bottom 
% norm(Pnt_new-Pnt_old)
if failed||(norm(Pnt_new-Pnt_old)<ProbStrut.knots_min)||...
        (norm(Pnt_new-Pnt_old)>ProbStrut.knots_max)||...
        ((ProbStrut.y(Pnt_new)-ProbStrut.y(Pnt_old))*Dir<0) % condition for increasing/decreasing in output
looping=false;
elseif Dir>0
InvFuntions.x(:,end+1)=Pnt_new;
InvFuntions.y(:,end+1)=ProbStrut.y(Pnt_new);
InvFuntions.J(:,end+1)=ProbStrut.J(Pnt_new);
else
InvFuntions.x=[Pnt_new,InvFuntions.x];
InvFuntions.y=[ProbStrut.y(Pnt_new),InvFuntions.y];
InvFuntions.J=[ProbStrut.J(Pnt_new),InvFuntions.J];
end 

end %finishing extension, going to next direction


end % finished doing positive and negative side
end


%% output changing step 
function [Pnt_new, Pnt_step] = InvFun_OutStep(Pnt, ProbStrut, Dir)
%INVFUN_OUTSTEP Produces a step in the output direction
% Pnt, is the current point
% Dir, indicates the positive or negative direction (cannot = 0) 

%% run loop around checking boundary violations
looping=true;
Voilate_dims=false(size(ProbStrut.x_min));
Pnt_new=Pnt;
while looping
Pnt_old=Pnt_new;
%% Calculate and apply step 
y_grad=sign(Dir)*ProbStrut.gamma_y*...
    ProbStrut.sat(ProbStrut.y_Grad(Pnt_old).*~Voilate_dims,...
    ProbStrut.beta_y); % equation #6 
Pnt_new=Pnt_old+y_grad; %equation #7 or 8
%% Check for domain violation
Violate=max(max(ProbStrut.x_min-Pnt_new,Pnt_new-ProbStrut.x_max),0);
[Violate_reduce,Violate_dim]=max(Violate);
if Violate_reduce<=ProbStrut.BoundViolateOK
    looping=false; %no violations, so end
else
    Voilate_dims(Violate_dim)=true;
    y_grad=y_grad*(1-abs(Violate_reduce/y_grad(Violate_dim)));
    Pnt_new=Pnt_old+y_grad;
end

end
Pnt_new=max(min(Pnt_new,ProbStrut.x_max),ProbStrut.x_min);
Pnt_step=Pnt_new-Pnt;
end

Contact us