gtFitJacobianPattern.m

function JacPatt = gtFitJacobianPattern(newvars, spot, gr, par)
% GTFITJACOBIANPATTERN Jacobian matrix structure for fitting.
%
% The JacPatt matrix is sparse and created by collecting all indices which 
% will be non-zero.
%
% Jacobian has the following structure (ng = no. of grains; ns = no. of 
% spots):
%   Each target variable (deviations in u,v,w) represents a row: 
%     [all du; all dv, all dw], size (3*ng,1)
%   Each fitted variable represents a column, ordered this way:
%     - global variables (geometry, constant energy, etc.)
%     - drift variables:
%         - X position
%         - Y position
%         - Z position
%         - X rotation
%         - Y rotation
%         - Z rotation
%         - energy
%     - grain centroid positions
%     - grain orientation
%     - grain strain components
%     Grain parameters are ordered by grain. 
%
% OUTPUT
%   JacPatt - Jacobian pattern (logical sparse matrix)
%
%
% Version 001 18-04-2014 by P.Reischig
%

%%%%%%%%%%%%%%%%%%%%%%
%% Global variables
%%%%%%%%%%%%%%%%%%%%%%
% Matrix size for global parameters: 3*numtarg x numglob.
% First build the matrix only for the spot U corrdinates, then repeat it 
% another 2 times.

% No. of target parameters 
numspot = sum([gr(:).numspots]);

% Indices of fitted global parameters (numglob: no. of fitted global parameters) 
indglob = 1:par.nvarglobal;

% Row indices
ir = (1:numspot)';
ir = ir(:, ones(1,length(indglob)));
ir = reshape(ir, [], 1);

% Column indices
ic = indglob(ones(numspot,1), :);
ic = reshape(ic, [], 1);

% No. of all non-zero components
nonz = length(ir);


%%%%%%%%%%%%%%%%%%%%%%
%% Drift variables
%%%%%%%%%%%%%%%%%%%%%%

switch par.drift.method

    case 'linear'
        error('Linear interpolation not yet supported.')
      
    case 'pchip'

        ird = cell(1,length(par.drift.omega));
        nom = length(par.drift.omega);
        
        % ROW INDICES
        
        % Omega range low and high boundaries where the drift variables 
        % have an effect (in degrees)
        omlo(2)     = par.drift.omega(1);
        omlo(nom-1) = par.drift.omega(nom-3);
        omlo(nom)   = par.drift.omega(nom-2);
        
        omhi(2)     = par.drift.omega(4);
        omhi(nom-1) = par.drift.omega(nom);
        omhi(nom)   = par.drift.omega(nom);
        
        for ii = 3 : (nom-2)           
            omlo(ii) = par.drift.omega(ii-2);
            omhi(ii) = par.drift.omega(ii+2);           
        end
        
        % Change from degrees to images
        omlo = omlo/par.omstep;
        omhi = omhi/par.omstep;
        
        % Find spot indices in the omega ranges (based on measured omega);
        % ird is indexed according to fitted variables, not to omega ranges
        for ii = 1:(nom-1)
            ird{ii} = find((omlo(ii+1) < spot.mesw) & (spot.mesw < omhi(ii+1)));     
        end
        
        
        % COLUMN INDICES
        
        % Active drift modes
        fitv = [par.fit_drift_posX, par.fit_drift_posY, par.fit_drift_posZ, ...
                par.fit_drift_rotX, par.fit_drift_rotY, par.fit_drift_rotZ,...
                par.fit_drift_energy];
           
        % Last column index    
        lastcolind  = par.nvarglobal;
        
        % Loop through drift modes
        for ii = 1:length(fitv) 
            if ~fitv(ii)
                continue
            end
            
            % Column indices of current drift variables 
            icd = (lastcolind + 1) : (lastcolind + nom - 1) ;
            
            % Current maximum column index
            lastcolind = icd(end);
            
            % Loop through time point (omega) variables within drift mode
            for jj = 1:length(icd)               
                % Add row indices to list
                ir = [ir; ird{jj}];
                
                % Add the same number of column indices to list
                icd_jj = icd(jj);
                ic     = [ic; icd_jj(ones(length(ird{jj}),1),1)]; 
            end
        end
          
    otherwise
        error('Interpolation method not recognised.')
        
end



% error
% add drift vars here
% par.drift

%%%%%%%%%%%%%%%%%%%%%%
%% Grain variables
%%%%%%%%%%%%%%%%%%%%%%
% Loop through grains, get corresponding row and column indices of 
% Jacobian for each grain:
%   - position
%   - orientation
%   - strain (if any)
% Assemble them in a vector.

if par.fit_grain_position || par.fit_grain_orientation || par.fit_grain_strain
    
    for ii = 1:length(gr)
        % Spot index tells which row indices (spots) correspond to grain
        irg = gr(ii).indspot';

        if isempty(gr(ii).indstrain)
            % No strain, only position and orientation
            irg = irg(:, ones(1,6));   % repeat 6 times
            
            icg = [gr(ii).indcenter, gr(ii).indRvec];
            icg = icg(ones(size(irg,1),1), :);
        else
            % Position, orientation, strain 
            irg = irg(:, ones(1,12));  % repeat 12 times
            
            icg = [gr(ii).indcenter, gr(ii).indRvec, gr(ii).indstrain];
            icg = icg(ones(size(irg,1),1), :);
        end
        
        ir = [ir; reshape(irg, [], 1)];
        ic = [ic; reshape(icg, [], 1)];
        
        nonz = nonz + numel(irg);
    end
    
end


% Create sparse matrix from indices and component values
JacPatt0 = sparse(ir, ic, true, numspot, length(newvars)); %, nonz);

% Repeat to extend to V and W:
JacPatt  = [JacPatt0; JacPatt0; JacPatt0];

    
end % function