Skip to content

Commit

Permalink
cleaned up parallel step tremendously #84
Browse files Browse the repository at this point in the history
  • Loading branch information
@timueh
timueh committed Apr 17, 2020
1 parent 06b0f2f commit 3464fe0
Showing 1 changed file with 95 additions and 67 deletions.
162 changes: 95 additions & 67 deletions src/core/parallelStepCentral.m
View file
Original file line number Diff line number Diff line change
@@ -1,79 +1,102 @@
function [ timers, opts, iter ] = parallelStepCentral( sProb, iter, timers, opts )
global use_fmincon
%PARALLELSTEP Summary of this function goes here
NsubSys = length(sProb.AA);

% solve local problems
for j=1:NsubSys % parfor???
nngi{j} = size(sProb.locFunsCas.ggi{j},1);
nnhi{j} = size(sProb.locFunsCas.hhi{j},1);
global use_fmincon
NsubSys = length(sProb.AA);

loc = [];
% solve local problems
for j=1:NsubSys % parfor???
[Neq, Nineq] = get_constraint_dimensions(sProb, j);
pNum = setup_parameters(sProb, iter, j);

tic
fprintf('\n\nSolving NLP in region %i\n', j);
[x0, z, rho, lambda, Sigma, problem] = unpack(sProb, iter, opts, j);
sol = problem.solve_nlp(x0, z, rho, lambda, Sigma, problem.pars);
loc = assign_solution(loc, sol, j);
timers.NLPtotTime = timers.NLPtotTime + toc;
loc = detect_active_set(sProb, loc, opts, j);
[loc, opts] = update_sigma(loc, iter, opts, j);

tic
fprintf('\nEvaluating sensitivities in region %i\n', j);
[loc, iter, Jac] = evaluate_sensitivities(sProb, loc, iter, opts, j);
timers.sensEvalT = timers.sensEvalT + toc;

fprintf('\nComputing reduced QP in region %i\n', j);
[loc, iter, timers] = compute_reduced_qp(sProb, loc, iter, opts, timers, Jac, j);
end

loc = save_for_bfgs(loc, opts);
iter.loc = loc;
end

function [Neq, Nineq] = get_constraint_dimensions(prob, j)
x = prob.zz0{j};
eq = prob.locFuns.ggi{j};
ineq = prob.locFuns.hhi{j};

Neq = length(eq(x));
Nineq = length(ineq(x));
end

function p = setup_parameters(prob, iter, j)
% set up parameter vector for local NLPs
p = [ iter.stepSizes.rho; iter.lam; iter.yy{j}];

% set up parameter vector for local NLP's
if ~isfield(sProb, 'p')
pNum = [ iter.stepSizes.rho;
iter.lam;
iter.yy{j}];
else
pNum = [ iter.stepSizes.rho;
iter.lam;
iter.yy{j};
sProb.p];
if isfield(prob, 'p')
p = [p; prob.p];
end
end


% solve local NLP's
tic
%sol = sProb.nnlp{j}('x0' , iter.yy{j},... this should be correct
%but is not working for BFGS example

function [x0, z, rho, lambda, Sigma, problem] = unpack(prob, iter, opts, j)
x0 = iter.yy{j};
z = iter.yy{j};
rho = iter.stepSizes.rho;
lambda = iter.lam;
Sigma = sparse(opts.SSig{j});

fprintf('Solving NLP in region %i\n', j);
problem = sProb.nnlp{j};
sol = problem.solve_nlp(x0, z, rho, lambda, Sigma, problem.pars);

% collect variables
[ loc.xx{j}, loc.KKapp{j}, loc.LLam_x{j} ] = deal(full(sol.x), ...
full(sol.lam_g), full(sol.lam_x));
timers.NLPtotTime = timers.NLPtotTime + toc;
problem = prob.nnlp{j};
end

% primal active set detection
loc.inact{j} = logical([false(nngi{j},1); ...
full(sProb.locFuns.hhi{j}(loc.xx{j}) < opts.actMargin)]);
KKapp{j}(loc.inact{j}) = 0;
function loc = assign_solution(loc, sol, j)
[ loc.xx{j}, loc.KKapp{j}, loc.LLam_x{j} ] = deal(full(sol.x), full(sol.lam_g), full(sol.lam_x));
end

function loc = detect_active_set(prob, loc, opts, j)
x = loc.xx{j};
ineq_eval = prob.locFuns.hhi{j}(x);
eq_eval = prob.locFuns.ggi{j}(x);

loc.inact{j} = logical([false(size(eq_eval)); full(ineq_eval < opts.actMargin)]);
end

function [loc, opts] = update_sigma(loc, iter, opts, j)
% dynamically changing \Sigma?
if strcmp(opts.Sig,'dyn')
if strcmp(opts.Sig, 'dyn')
% after second iteration
if iter.i > 1
[opts.SSig{j}, loc.locStep{j}] = computeDynSig(opts.SSig{j},...
iter.yy{j} - loc.xx{j},iter.loc.locStep{j}, 'Sig');
[opts.SSig{j}, loc.locStep{j}] = computeDynSig(opts.SSig{j}, iter.yy{j} - loc.xx{j},iter.loc.locStep{j}, 'Sig');
else
loc.locStep{j} = iter.yy{j} - loc.xx{j};
end
end
end

% evaluate sensitivities locally
tic
% gradient of the objective
loc.sensEval.ggiEval{j} = sProb.sens.gg{j}(loc.xx{j});
function loc = evaluate_gradient(prob, loc, j)
loc.sensEval.ggiEval{j} = prob.sens.gg{j}(loc.xx{j});
end

function [loc, iter] = evaluate_hessian(prob, loc, opts, iter, j)
% Hessian approximations
if strcmp(opts.Hess, 'BFGS') || strcmp(opts.Hess, 'DBFGS')
loc.sensEval.gLiEval{j} = sProb.sens.gL{j}(loc.xx{j},loc.KKapp{j});
loc.sensEval.gLiEval{j} = prob.sens.gL{j}(loc.xx{j},loc.KKapp{j});
if ~isfield(iter.loc, 'sensEval')
if strcmp(opts.BFGSinit, 'ident')
% initialize BFGS with identity matrix
loc.sensEval.HHiEval{j} = eye(length(sProb.zz0{j}));
loc.sensEval.HHiEval{j} = eye(length(prob.zz0{j}));
elseif strcmp(opts.BFGSinit, 'exact')
% initialize BFGS with exact Hessian
loc.sensEval.HHiEval{j} = ...
sProb.sens.HH{j}(loc.xx{j},loc.KKapp{j},iter.stepSizes.rho);
prob.sens.HH{j}(loc.xx{j},loc.KKapp{j},iter.stepSizes.rho);
end
else
loc.sensEval.HHiEval{j} = BFGS(iter.loc.sensEval.HHiEval{j},...
Expand All @@ -91,7 +114,7 @@
iter.comm.globF.primVal{j} = [ iter.comm.globF.primVal{j} length(iter.loc.xx{j}) ];
end
else
loc.sensEval.HHiEval{j} = sProb.sens.HH{j}(loc.xx{j},loc.KKapp{j},iter.stepSizes.rho);
loc.sensEval.HHiEval{j} = prob.sens.HH{j}(loc.xx{j},loc.KKapp{j},iter.stepSizes.rho);

if strcmp(opts.commCount, 'true') && ~strcmp(opts.slack,'redSpace') && strcmp(opts.innerAlg, 'none')
% communication for xx and the gradient of the Lagrangian and
Expand All @@ -101,26 +124,34 @@
iter.comm.globF.primVal{j} = [ iter.comm.globF.primVal{j} length(iter.loc.xx{j}) ];
end
end
end

function [loc, JacCon] = evaluate_jacobian(prob, loc, opts, j)
% Jacobians of active nonlinear constraints/bounds
JacCon = full(sProb.sens.JJac{j}(loc.xx{j}));
JacCon = full(prob.sens.JJac{j}(loc.xx{j}));
JacBounds = eye(size(loc.xx{j},1));

% eliminate inactive entries
JJacCon{j} = sparse(JacCon(~loc.inact{j},:));
JacBounds = JacBounds((sProb.llbx{j} - loc.xx{j}) ...
> opts.actMargin |(loc.xx{j}-sProb.uubx{j}) > opts.actMargin,:);
loc.sensEval.JJacCon{j} = [JJacCon{j}; JacBounds];
timers.sensEvalT = timers.sensEvalT + toc;

JacBounds = JacBounds((prob.llbx{j} - loc.xx{j}) > opts.actMargin |(loc.xx{j}-prob.uubx{j}) > opts.actMargin,:);
loc.sensEval.JJacCon{j} = [JJacCon{j}; JacBounds];
end

function [loc, iter, Jac] = evaluate_sensitivities(prob, loc, iter, opts, j)
loc = evaluate_gradient(prob, loc, j);
[loc, iter] = evaluate_hessian(prob, loc, opts, iter, j);
[loc, Jac] = evaluate_jacobian(prob, loc, opts, j);
end

function [loc, iter, timers] = compute_reduced_qp(prob, loc, iter, opts, timers, Jac, j)
% for reduced-space method, compute reduced QP
if strcmp(opts.slack,'redSpace') && strcmp(opts.innerAlg, 'none')

loc.sensEval.ZZ{j} = null(full(JJacCon{j}));
loc.sensEval.ZZ{j} = null(full(Jac));
loc.sensEval.HHred{j} = loc.sensEval.ZZ{j}'* ...
full(loc.sensEval.HHiEval{j})*loc.sensEval.ZZ{j};

loc.sensEval.AAred{j} = sProb.AA{j}*loc.sensEval.ZZ{j};
loc.sensEval.AAred{j} = prob.AA{j}*loc.sensEval.ZZ{j};
loc.sensEval.ggred{j} = loc.sensEval.ZZ{j}'*full(loc.sensEval.ggiEval{j});

% regularize reduced Hessian
Expand Down Expand Up @@ -156,15 +187,12 @@
iter.comm.globF.Jac{j} = [ iter.comm.globF.Jac{j} sJ(1)*sJ(2) ];
end
end
end

% save information for next BFGS iteration
if strcmp(opts.Hess, 'BFGS') || strcmp(opts.Hess, 'DBFGS')
loc.sensEval.gLiEvalOld = loc.sensEval.gLiEval;
loc.xxOld = loc.xx;
end

iter.loc = loc;

end

function loc = save_for_bfgs(loc, opts)
% save information for next BFGS iteration
if strcmp(opts.Hess, 'BFGS') || strcmp(opts.Hess, 'DBFGS')
loc.sensEval.gLiEvalOld = loc.sensEval.gLiEval;
loc.xxOld = loc.xx;
end
end

0 comments on commit 3464fe0

Please sign in to comment.