Home > source > @penlab > eval_aldx.m

# eval_aldx

## PURPOSE

As eval_alx(), evaluate gradient of the Augmented Lagrangian and Jacobian

## SYNOPSIS

function [status] = eval_aldx(obj)

## DESCRIPTION

     As eval_alx(), evaluate gradient of the Augmented Lagrangian and Jacobian
     of the equality constraints. The results are stored in obj.ALdx, obj.eqdx
     and tickers updated.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SOURCE CODE

    0001 % As eval_alx(), evaluate gradient of the Augmented Lagrangian and Jacobian
    0002 % of the equality constraints. The results are stored in obj.ALdx, obj.eqdx
    0003 % and tickers updated.
    0004 function [status] = eval_aldx(obj)
    0005
    0006   % status TODO
    0007   status = 0;
    0008
    0009   if (obj.ALdxtck < obj.ticker)
    0010     starttime = cputime;
    0011
    0012     % create local copies of obj.x,obj.Y to avoid checking repetitively if they
    0013     % are 'up-to-date' with obj.xall, they get used many times in calls to
    0014     % user's functions
    0015     x=obj.x;
    0016     Y=obj.Y;
    0017     % same for obj.userdata, it might be quite expensive to get&store it in
    0018     % the object
    0019     userdata=obj.userdata;
    0020
    0021     % TODO reuse it + 'pointchanged' flag
    0022     [fx, userdata] = obj.objfun(x, Y, userdata);
    0023     [fdx, userdata] = obj.objgrad(x, Y, userdata);
    0024
    0025     % function inequal & equal
    0026     if (obj.NgNLN + obj.NgLIN>0)
    0027       [gx, userdata] = obj.confun(x, Y, userdata);
    0028       ineqx = obj.ineqshift + obj.ineqmlt .* gx(obj.ineqmap);
    0029       obj.eqx = obj.eqshift + gx(obj.eqmap);
    0030
    0031       [gdx, userdata] = obj.congrad(x, Y, userdata);
    0032       % because of stupid spdiags() which doesn't accept ([],0,0,0)
    0033       if (obj.Nineq>0)
    0034         ineqdx = gdx(:,obj.ineqmap) * spdiags(obj.ineqmlt,0,obj.Nineq,obj.Nineq);
    0035       else
    0036         ineqdx = [];
    0037       end
    0038       % TODO update ticker eqdx
    0039       obj.eqdx = gdx(:,obj.eqmap);
    0040     end
    0041
    0042     % box constraints
    0043     if (obj.Nxbox>0)
    0044       xboxx = obj.xboxshift + obj.xboxmlt .* obj.xall(obj.xboxmap);
    0045       xboxdx = sparse(obj.xboxmap,[1:obj.Nxbox],obj.xboxmlt,obj.Nx,obj.Nxbox);
    0046     end
    0047
    0048     ALdx=fdx;
    0049
    0050     % box constraints
    0051     ind=obj.xboxindbar;
    0052     if (~isempty(ind))
    0053       ALdx = ALdx + xboxdx(:,ind)*(obj.uxbox(ind).*obj.pxbox(ind).*obj.phibar_D(xboxx(ind)));
    0054     end
    0055
    0056     ind=obj.xboxindphi;
    0057     if (~isempty(ind))
    0058       ALdx = ALdx + xboxdx(:,ind)*(obj.uxbox(ind).*obj.phi2_D(xboxx(ind)./obj.pxbox(ind)));
    0059     end
    0060
    0061     % function inequalitites
    0062     %ind=obj.ineqindbar;
    0063     %... TODO
    0064
    0065     ind=obj.ineqindphi;
    0066     if (~isempty(ind))
    0067       ALdx = ALdx + ineqdx(:,ind)*(obj.uineq(ind).*obj.phi2_D(ineqx(ind)./obj.pineq(ind)));
    0068     end
    0069
    0070     % (function) equalities
    0071     if (obj.Neq>0)
    0072       ALdx = ALdx + obj.eqdx*obj.ueq;
    0073     end
    0074
    0075     % matrix variable - log barrier (strict feasibility)
    0076     if ~isempty(obj.Yboxindbar)
    0077     for k=obj.Yboxindbar
    0078       % convert the matrix box constraint to the form:   +/-Y +/-bound >=0
    0079       pkx=obj.PYbox(k);
    0080       Ykx = Y{obj.Yboxmap(k)};
    0081       Akx=-obj.Yboxshift(k)*speye(size(Ykx)) - obj.Yboxmlt(k)*Ykx;
    0082
    0083       % can assume that ALx was computed before ALdx --> Akx must be pos. def.
    0084       chol(Akx);
    0085       invAkx = full(inv(Akx));
    0086
    0087       % compute partial derivatives of the penalty term:
    0088       %    d/dxi... = -p*trace( (A^-1) * (d/dxi A) )
    0089
    0090       % there should be just one dependant with obvious value (+-1) in each
    0091       % triangle of every Y{k} derivative
    0092       mlt=-obj.Yboxmlt(k);  % +/-1
    0093       mapper=obj.vec2Ymap{obj.Yboxmap(k)};
    0094       dim=mapper.dim;
    0095       offset=obj.Nx + mapper.xmap(1) - 1;
    0096       irow=mapper.irow;
    0097       icol=mapper.icol;
    0098       for idx=1:mapper.nelem
    0099         if (irow(idx)==icol(idx))
    0100           % diagonal element
    0101           Akdx = sparse(irow(idx),icol(idx),mlt,dim,dim);
    0102         else
    0103           % nondiag element --> add two
    0104           Akdx = sparse([irow(idx),icol(idx)],[icol(idx),irow(idx)],[mlt,mlt],dim,dim);
    0105         end
    0106         % but this can be done directly...!
    0107         ALdx(offset+idx) = ALdx(offset+idx) - pkx*trace(invAkx*Akdx);
    0108       end
    0109
    0110     end
    0111     end
    0112
    0113     % matrix variable - pen/bar
    0114     if ~isempty(obj.Yboxindphi)
    0115     for k=obj.Yboxindphi
    0116       pkx=obj.PYbox(k);            %2*p(sdpdata.Ng+k);
    0117       Ykx = obj.Y{obj.Yboxmap(k)};
    0118       Akx=obj.Yboxshift(k)*speye(size(Ykx)) + obj.Yboxmlt(k)*Ykx;
    0119       umatk=obj.UYbox{k};
    0120
    0121       %Z=(pkx*speye(size(Akx))-Akx);
    0122       Z=(pkx*speye(size(Akx))-Akx);
    0123       %invZ=inv(Z);
    0124       invZ=full(inv(Z));
    0125       pZUZ=pkx^2*invZ*umatk*invZ;
    0126
    0127       % there should be just one dependant with obvious value (+-1) in each
    0128       % triangle of every Y{k} derivative
    0129       mlt=obj.Yboxmlt(k);  % +/-1
    0130       mapper=obj.vec2Ymap{obj.Yboxmap(k)};
    0131       dim=mapper.dim;
    0132       offset=obj.Nx + mapper.xmap(1) - 1;
    0133       irow=mapper.irow;
    0134       icol=mapper.icol;
    0135       for idx=1:mapper.nelem
    0136         if (irow(idx)==icol(idx))
    0137           % diagonal element
    0138           Akdx = sparse(irow(idx),icol(idx),mlt,dim,dim);
    0139         else
    0140           % nondiag element --> add two
    0141           Akdx = sparse([irow(idx),icol(idx)],[icol(idx),irow(idx)],[mlt,mlt],dim,dim);
    0142         end
    0143         % but this can be done directly...!
    0144         ALdx(offset+idx) = ALdx(offset+idx) + trace(pZUZ*Akdx);
    0145       end
    0146     end
    0147     end
    0148
    0149     % matrix constraints - pen/bar
    0150     if ~isempty(obj.Aindphi)
    0151     for k=obj.Aindphi
    0152       pkx=obj.PA(k);  % I used to use 2*         !!!!!!!!
    0153       % TODO need to map the matrix first! - is it correct???
    0154       kuser=obj.Amap(k);
    0155       [Akuserx, userdata] = obj.mconfun(x, Y, kuser, userdata);
    0156       Akx = obj.Ashift(k)*speye(size(Akuserx)) + obj.Amlt(k) .* Akuserx;
    0157       umatk=obj.UA{k};
    0158
    0159       %Z=(pkx*speye(size(Akx))-Akx);
    0160       Z=(pkx*speye(size(Akx))-Akx);
    0161       %invZ=inv(Z);
    0162       invZ=full(inv(Z));
    0163       pZUZ=pkx^2*invZ*umatk*invZ;
    0164
    0165       % MAPPING & transformation !!!
    0166       for i=obj.Adep{kuser}
    0167         [Akdx, userdata] = obj.mcongrad(x,Y,kuser,i,userdata);
    0168         Akdx = obj.Amlt(k)*Akdx;
    0169         %ALdx(i) = ALdx(i) + trace(pZUZ*Akdx);
    0170         ALdx(i) = ALdx(i) + pZUZ(:)'*Akdx(:);
    0171       end
    0172     end
    0173     end
    0174
    0175     % store user's data back in the object if it got changed
    0176     obj.userdata=userdata;
    0177
    0178     % update ticker
    0179     obj.ALdx = ALdx;
    0180     obj.ALdxtck = obj.ticker;
    0181
    0182     % update stats
    0183     obj.stats_ncall_aldx = obj.stats_ncall_aldx + 1;
    0184     obj.stats_time_aldx = obj.stats_time_aldx + cputime - starttime;
    0185
    0186   end
    0187
    0188