function g = matRad_daoGradFunc(apertureInfoVec,dij,cst,options)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% matRad IPOPT callback: gradient function for direct aperture optimization
%
% call
%   g = matRad_daoGradFunc(apertureInfoVec,apertureInfo,dij,cst,type)
%
% input
%   apertureInfoVec:  aperture info in form of vector
%   dij:              matRad dij struct as generated by bixel-based dose calculation
%   cst:              matRad cst struct
%   options:          option struct defining the type of optimization
%
% output
%   g: gradient
%
% References
%   [1] http://dx.doi.org/10.1118/1.4914863
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% read in the global apertureInfo and apertureVector variables
global matRad_global_apertureInfo;
% update apertureInfo from the global variable
apertureInfo = matRad_global_apertureInfo;

% update apertureInfo, bixel weight vector an mapping of leafes to bixels
if ~isequal(apertureInfoVec,apertureInfo.apertureVector)
    apertureInfo = matRad_daoVec2ApertureInfo(apertureInfo,apertureInfoVec);
    matRad_global_apertureInfo = apertureInfo;
end

% bixel based gradient calculation
bixelG = matRad_gradFuncWrapper(apertureInfo.bixelWeights,dij,cst,options);

% allocate gradient vector for aperture weights and leaf positions
%changed from NaNs to zeros
g = zeros(size(apertureInfoVec,1),1);

if apertureInfo.runVMAT
    
    % use the Jacobian calculated in daoVec2ApertureInfo.
    % should also do this for non-VMAT
    g = g+apertureInfo.bixelJApVec * bixelG;
else
    %we're not doing VMAT
    
    % 1. calculate aperatureGrad
    % loop over all beams
    offset = 0;
    for i = 1:numel(apertureInfo.beam)
        
        % get used bixels in beam
        ix = ~isnan(apertureInfo.beam(i).bixelIndMap);
        
        % loop over all shapes and add up the gradients x openingFrac for this shape
        for j = 1:apertureInfo.beam(i).numOfShapes
            g(j+offset) = apertureInfo.beam(i).shape(j).shapeMap(ix)' ...
                * bixelG(apertureInfo.beam(i).bixelIndMap(ix))./apertureInfo.beam(i).shape(j).jacobiScale;
        end
        
        % increment offset
        offset = offset + apertureInfo.beam(i).numOfShapes;
        
    end
    
    % 2. find corresponding bixel to the leaf Positions and aperture
    % weights to calculate the gradient
    g(apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs*2) = ...
        apertureInfoVec(apertureInfo.mappingMx(apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs*2,2)) ...
        .* bixelG(apertureInfo.bixelIndices(1:apertureInfo.totalNumOfLeafPairs*2)) ./ ...
        (apertureInfo.bixelWidth.*apertureInfo.jacobiScale(apertureInfo.mappingMx(apertureInfo.totalNumOfShapes+1:apertureInfo.totalNumOfShapes+apertureInfo.totalNumOfLeafPairs*2,2)));
    
    % correct the sign for the left leaf positions
    g(apertureInfo.totalNumOfShapes+(1:(apertureInfo.totalNumOfLeafPairs))) = ...
        -g(apertureInfo.totalNumOfShapes+(1:(apertureInfo.totalNumOfLeafPairs)));
    
end


end