text:
function imageCha = recon_main(obj, iRep, iAvg)
% main reconstruction function
% prepare constraints and select algorithm depending on image dimensionality
%
% input:
% obj       CS reconstruction object (holding all parameters)
% iRep      current repetition
% iAvg      current average
%
% output:
% imageCha  individual reconstructed channel images
%
% (c) Thomas Kuestner
% -------------------------------------------------------------------------

nCha = obj.measPara.dim(5);
nSlices = size(obj.kSpace,1);
nTime = obj.measPara.dim(4);
imageCha = cell(nSlices,nCha);

if(strcmp(obj.measPara.dimension,'2D'))

elseif(strcmp(obj.measPara.dimension,'3D'))

    % prepare k-space
    % => k_y - k_z - k_x - cha
    kSpaceL = cell2mat(shiftdim(obj.kSpace(:,:,iRep,iAvg),-2));
    kSpaceL = permute(kSpaceL,[1 3 2 4]);
    if(strcmp(obj.measPara.precision,'double'))
        kSpaceL = double(kSpaceL);
    end
    % same undersampling for all coils and fully sampled in k_x direction
    obj.fullMask = abs(kSpaceL) > 0; % k_y - k_z - k_x
    if(obj.measPara.oversampling{2,1})
        kSpaceL = ifftnshift(kSpaceL, 3); % -> k_y - k_z - x - cha
        kSpaceL = kSpaceL(:, :, obj.measPara.oversampling{1,1}, :); % anti-aliasing
        obj.fullMask = obj.fullMask(:,:, obj.measPara.oversampling{1,1}, :);
        kSpaceL = fftnshift(kSpaceL,3); % -> k_y - k_z - k_x - cha
        kSpaceL = kSpaceL .* obj.fullMask;
        obj.measPara.dim(2) = size(kSpaceL,3);
    end
    [nPha, nZ, nFreq, nCha] = size(kSpaceL);

    if(~isempty(obj.calibSize))
        % adjust calibSize due to anti-aliasing
        obj.calibSize(3) = nFreq;
    end

    % create a 3D kernel fomr the calibration region
    fprintf('Performing calibration\n');
    if(isempty(obj.calibSize))
        % automatically determine calibration region
        obj.calibSize = obj.calibrationSize(obj.fullMask(:,:,:,1)); % k_y - k_z - k_x
        if(size(kSpaceL,3) == 1), obj.calibSize = [obj.calibSize, 1]; obj.kernelSize = obj.kernelSize([1 3 2]); end; % k_y - k_x sparse!
        if(~obj.flagToolbox)
            obj.kCalib = crop(kSpaceL, [obj.calibSize, nCha]); % k_y - k_z - k_x - cha
        end
    else
        if(~obj.flagToolbox)
            % fixed calibration region
            plonk = [nPha, nZ, nFreq];
            idx = cell(1,length(plonk));
            for iInner=1:length(plonk)
                if(mod(obj.calibSize(iInner),2) == 0)
                    idx{iInner} = floor(plonk(iInner)/2)+1+ceil(-obj.calibSize(iInner)/2):floor(plonk(iInner)/2)+ceil(obj.calibSize(iInner)/2);
                else
                    idx{iInner} = floor(plonk(iInner)/2)+ceil(-obj.calibSize(iInner)/2):floor(plonk(iInner)/2)+ceil(obj.calibSize(iInner)/2)-1;
                end

                tmp = [idx{iInner}(1) <= 0, idx{iInner}(end) > plonk(iInner)];
                if(any(tmp))
                    if(all(tmp)), error('crop(): Both index out of bounds'); end;
                    hShift = [abs(idx{iInner}(1)) + 1, idx{iInner}(end) - plonk(iInner)];
                    op = {idx{iInner}(1), idx{iInner}(end); '+', '-'; '<', '>'; plonk(iInner) + 1, 0};
                    eval(sprintf('if(op{1,~helper} %s hShift(helper) %s %d), idx{iInner} = idx{iInner} %s hShift(helper); else idx{iInner} = idx{iInner}(idx{iInner} %s %d);end;',op{2,helper}, op{3,helper}, op{4,helper}, op{2,helper}, op{3,~helper}, op{4,~helper}));
                end
            end
            obj.kCalib = kSpaceL(idx{1},idx{2},idx{3},:); % k_y - k_z - k_x - cha
        end
    end

    if(~obj.flagToolbox)
        % windowing for cut-out region
    %     [~,w1] = window2D(prop.window.type,size(obj.kCalib,1),prop.window.windowOpt{1},prop.window.windowOpt{2});
    %     [~,w2] = window2D(prop.window.type,size(obj.kCalib,2),prop.window.windowOpt{1},prop.window.windowOpt{2});
    %     window =  w1(:) * w2(:).';
        if(obj.window.windowingOn)
            window = windowND(obj.window.type,[size(obj.kCalib,1), size(obj.kCalib,2)],obj.window.windowOpt{1},obj.window.windowOpt{2});
            obj.kCalib = obj.kCalib .* repmat(window,[1 1 size(obj.kCalib,3) nCha]);
        end

        % compute eigen-value maps
        [hybridKernel,S] = dat3Kernel(obj.kCalib,obj.kernelSize);
        hybridKernel = ifftnshift(zpad(hybridKernel(:,:,end:-1:1,:,:),[obj.kernelSize(1),obj.kernelSize(2),nFreq + obj.kernelSize(3)-1,nCha,size(hybridKernel,5)]),3);
        idx = max(find(S >= S(1)*obj.eigThresh_k));

        % compute inverse fourier transform along readout direction
        kSpaceL = ifftnshift(kSpaceL,3); % k_y - k_z - x - cha

        % convolution with GRAPPA kernel in k-space -> multiplication in image
        % domain
        obj.method = 'fft';

        % prepare wavelet transformation
        if(strcmp(obj.solver,'L1'))
            if(strcmp(obj.trafo.trafoType,'wavelet_lab'))
                ssx = 2^ceil(log2(nPha));
                ssy = 2^ceil(log2(nZ));
                ss = max(ssx, ssy);
                XOP = Wavelet(obj.trafo.waveletFilter,obj.trafo.waveletFilterSize,log2(ss/2^(obj.trafo.waveletStages)));
            end
        end

        image = zeros(nPha,nFreq,nZ,obj.n_maps);
    %     AtA_sub = zeros(prod(obj.kernelSize)*nCha); % subtract this matrix from AtA_mask, due to accumulation of previous samples (A'*A multiplication)
    %     A_sub = zeros((size(obj.kCalib,1)-obj.kernelSize(1)+1) * (size(obj.kCalib,2)-obj.kernelSize(2)+1) * (size(obj.kCalib,3)-obj.kernelSize(3)+1),prod(obj.kernelSize)*nCha);

        fprintf('ESPIRiT reconstruction\n');
        dispProgress('Frequency',0,nFreq);
        for iFreq=1:nFreq

            % extract 2D kernel from complete 3D kernel
            obj.kernel = squeeze(hybridKernel(:,:,iFreq,:,:));

            % crop kernels and compute eigen-value decomposition in image space to get maps
            [M,W] = kernelEig(obj.kernel(:,:,:,1:idx),[nPha,nZ]);

            % Compute Soft-SENSE ESPIRiT Maps
            % crop sensitivity maps according to eigenvalues == 1
            % weigth with n maps of eigen-values

            maps = M(:,:,:,end-(obj.n_maps-1):end);

            % Weight the eigenvectors with soft-senses eigen-values
            weights = W(:,:,end-(obj.n_maps-1):end);
            weights = (weights - obj.eigThresh_im)./(1-obj.eigThresh_im).* (W(:,:,end-(obj.n_maps-1):end) > obj.eigThresh_im);
            weights = -cos(pi*weights)/2 + 1/2;

            % create an ESPIRiT operator
            ESP = ESPIRiT(maps,weights);

            if(strcmp(obj.solver,'L1'))
                FT = p2DFT(squeeze(obj.fullMask(:,:,iFreq,:)),[nPha,nZ,nCha]);
                if(strcmp(obj.trafo.trafoType,'fft'))
                    XOP = FT;
                end
                image(:,iFreq,:,:) = permute(obj.cgL1ESPIRiT(squeeze(kSpaceL(:,:,iFreq,:)), zeros([size(squeeze(kSpaceL(:,:,iFreq,1))),obj.n_maps]), FT, ESP, obj.iNINNER, XOP, obj.trafo.wavWeight, obj.splitWeight, obj.iNIterSplit),[1 4 2 3]); % => y-x-z-cha
            else
                [tmpK,tmp] = obj.cgESPIRiT(squeeze(kSpaceL(:,:,iFreq,:)), ESP, obj.iNINNER, obj.calibTyk, zeros(size(squeeze(kSpaceL(:,:,iFreq,:)))));
                image(:,iFreq,:,:) = permute(tmp,[1 4 2 3]);
            end
            dispProgress('Frequency', iFreq/nFreq);
        end
        dispProgress('Frequency','Close');
    %     image = ifftnshift(image,[1,3]);
        imageCha(:,1:obj.n_maps) = shiftdim(mat2cell(image,nPha,nFreq,nZ,ones(1,obj.n_maps)),2);
        [imageCha{:,obj.n_maps+1:end}] = deal(zeros(size(imageCha{1,1})));

    else
        fprintf('ESPIRiT reconstruction\n');
        if(ispc || ismac)
            eval(sprintf('ESP = bart(''ecalib -r %d:%d -k %d:%d -m %d'',kSpaceL);', obj.calibSize(1), obj.calibSize(2), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps ));
            if(strcmp(obj.solver,'L1'))
                eval(sprintf('image = permute(squeeze(bart(''pics -l1 -r %f'', kSpaceL, ESP)),[1 3 2 4]);', obj.trafo.wavWeight));
            else
                image = permute(squeeze(bart('pics -l2', kSpaceL, ESP)),[1 3 2 4]);
            end
        else
            % create temporary dir if not yet existing
            if(~exist([obj.path,filesep,'tmpESPIRiT'],'dir'))
                mkdir([obj.path,filesep,'tmpESPIRiT']);
            end
            fftdims = size(kSpaceL);
            fftdims = fftdims(1:3);
            if(any(mod(fftdims,2) ~= 0))
                writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], zpad(kSpaceL,size(kSpaceL) + [mod(fftdims,2),0])); % k_y - k_z - k_x - cha
            else
                writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL); % k_y - k_z - k_x - cha
            end
            obj.calibSize = obj.calibSize - mod(obj.calibSize,2);
            eval(sprintf('! ecalib -r %d:%d:%d -k %d:%d:%d -m %d %s %s', obj.calibSize(1), obj.calibSize(2), obj.calibSize(3), obj.kernelSize(1), obj.kernelSize(2), obj.kernelSize(3), obj.n_maps, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps']));
            ESP = [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps'];

            if(strcmp(obj.solver,'L1'))
                eval(sprintf('! sense -l1 -r %f %s %s %s', obj.trafo.wavWeight, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
            else
                eval(sprintf('! sense -l2 %s %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
            end
            image = permute(squeeze(readcfl([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut'])),[1 3 2 4]);
        end
%         % OR:
%         kSpaceL = permute(ifftnshift(kSpaceL,3), [3 1 2 4]); % x - k_y - k_z - cha
%         writecfl([obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], kSpaceL); % x - k_y - k_z - cha
%         eval(sprintf('! ecalib -r %d:%d:%d -k %d:%d:%d -m %d %s %s', obj.calibSize(3), obj.calibSize(1), obj.calibSize(2), obj.kernelSize(3), obj.kernelSize(1), obj.kernelSize(2), obj.n_maps, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps']));
%         ESP = [obj.path,filesep,'tmpESPIRiT',filesep,'espiritmaps'];
%          if(strcmp(obj.solver,'L1'))
%             eval(sprintf('! rsense -l1 -r %f %s %s %s', obj.trafo.wavWeight, [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
%         else
%             eval(sprintf('! rsense -l2 %s %s %s', [obj.path,filesep,'tmpESPIRiT',filesep,'kspace'], ESP, [obj.path,filesep,'tmpESPIRiT',filesep,'imgOut']));
%         end
%         image = permute(squeeze(readcfl([obj.path,filesep,'tmpESPIRiT',filesep,'imgOut'])), [2 1 3 4]);

        image = crop(image,[size(kSpaceL,1),size(kSpaceL,3),size(kSpaceL,2),obj.n_maps]);
        imageCha(:,1:obj.n_maps) = shiftdim(mat2cell(image,nPha,nFreq,nZ,ones(1,obj.n_maps)),2);
        [imageCha{:,obj.n_maps+1:end}] = deal(zeros(size(imageCha{1,1})));
    end

elseif(strcmp(obj.measPara.dimension,'4D'))

    % prepare k-space
    % => k_x - k_y - k_z - cha - ... - t
    kSpaceL = cell2mat(shiftdim(obj.kSpace(:,:,iRep,iAvg),-3));
    kSpaceL = permute(kSpaceL,[2 1 3 5 4]);
    if(strcmp(obj.measPara.precision,'double'))
        kSpaceL = double(kSpaceL);
    end
    % same undersampling for all coils and fully sampled in k_x direction
    obj.fullMask = abs(kSpaceL) > 0; % k_x - k_y - k_z - cha - t
    if(obj.measPara.oversampling{2,1})
        kSpaceL = ifftnshift(kSpaceL, 1);
        kSpaceL = kSpaceL(obj.measPara.oversampling{1,1}, :, :, :, :); % anti-aliasing
        obj.fullMask = obj.fullMask(obj.measPara.oversampling{1,1}, :, :, :, :);
        kSpaceL = fftnshift(kSpaceL,1);
        kSpaceL = kSpaceL .* obj.fullMask;
        obj.measPara.dim(2) = size(kSpaceL,1);
    end

    if(~isempty(obj.calibSize))
        % adjust calibSize due to anti-aliasing
        obj.calibSize(3) = size(kSpaceL,1);
    end

    lEnsureCenter = false;
    calibSize = zeros(size(obj.fullMask,5),3);
    for iTime = 1:size(obj.fullMask,5)
        calibSize(iTime,:) = obj.calibrationSize(obj.fullMask(:,:,:,1,iTime));
        if(any(calibSize(iTime,:) <= 4)) % or <= 2
            lEnsureCenter = true;
        end
    end

    if(lEnsureCenter) % ensure a center region is found!
        dPercentage = [0.065,0.2];
        m = [size(obj.fullMask,1), size(obj.fullMask,2), size(obj.fullMask,3)];
        s = [size(obj.fullMask,1)/2, round(dPercentage(1)*size(obj.fullMask,2)), round(dPercentage(2)*size(obj.fullMask,3))];
        idx = cell(1,length(s));
        for n=1:length(s)
            if(mod(s(n),2) == 0)
                idx{n} = floor(m(n)/2)+1+ceil(-s(n)/2) : floor(m(n)/2)+ceil(s(n)/2);
            else
                idx{n} = floor(m(n)/2)+ceil(-s(n)/2) : floor(m(n)/2)+ceil(s(n)/2)-1;
            end
        end
        lMaskNew = false(size(obj.fullMask));
        lMaskNew(idx{1},idx{2},idx{3},:,:) = true;
        lMaskNew = lMaskNew & ~obj.fullMask;
        kSpaceL(lMaskNew) = eps;
    end

    if(obj.mc.jointMCCS)
        % coarse motion estimation (in lowres image)
%         kSpaceL = squeeze(kSpaceL); % x-y-z-cha-t
        nTime = size(kSpaceL,5);
%         % get SENSE map
%         pos = 'idx{1},idx{2},idx{3},:,t';
%         m = size(kSpaceL);
%         kSpaceCenter = false(size(kSpaceL));
%         for t=1:nTime
%             s = calibSize(t,:);
%             idx = cell(1,length(s));
%             for n=1:length(s)
%                 if(mod(s(n),2) == 0)
%                     idx{n} = floor(m(n)/2)+1+ceil(-s(n)/2) : floor(m(n)/2)+ceil(s(n)/2);
%                 else
%                     idx{n} = floor(m(n)/2)+ceil(-s(n)/2) : floor(m(n)/2)+ceil(s(n)/2)-1;
%                 end
%
%                 tmp = [idx{n}(1) <= 0, idx{n}(end) > m(n)];
%                 if(any(tmp))
%                     if(all(tmp)), error('crop(): Both index out of bounds'); end;
%                     hShift = [abs(idx{n}(1)) + 1, idx{n}(end) - m(n)];
%                     op = {idx{n}(1), idx{n}(end); '+', '-'; '<', '>'; m(n) + 1, 0};
%                     eval(sprintf('if(op{1,~tmp} %s hShift(tmp) %s %d), idx{n} = idx{n} %s hShift(tmp); else idx{n} = idx{n}(idx{n} %s %d);end;',op{2,tmp}, op{3,tmp}, op{4,tmp}, op{2,tmp}, op{3,~tmp}, op{4,~tmp}));
%                 end
%             end
%             eval(sprintf('kSpaceCenter(%s) = true;',pos));
%         end
%
%         kSpaceLow = zeros(size(kSpaceL)); kSpaceLow(kSpaceCenter) = kSpaceL(kSpaceCenter);
%         dImgCenter = ifftnshift(kSpaceLow,1:3); % x-y-z-cha-t
%         clear 'kSpaceLow';
%         dImgCenterSOS = sqrt(sum(abs(dImgCenter).^2,4)); % x-y-z-1-t
%         dSensemap = dImgCenter./(repmat(dImgCenterSOS,[1 1 1 size(dImgCenter,4) 1]));
%         clear 'dImgCenterSOS' 'dImgCenter';
%         dWindow = windowND(@hamming, [size(dSensemap,1), size(dSensemap,2), size(dSensemap,3)]);
%         dSensemap = fftnshift(dSensemap,1:3);
%         dSensemap = dSensemap .* repmat(dWindow, [1 1 1 size(dSensemap,4) size(dSensemap,5)]);
%         dSensemap = ifftnshift(dSensemap,1:3);
%         dSensemapFctr = sqrt(sum(conj(dSensemap).*dSensemap,4)).^(-1);
%         dSensemapFctr(isinf(dSensemapFctr)) = 0;
%         dSensemap = conj(repmat(dSensemapFctr,[1 1 1 size(dSensemap,4) 1]) .* dSensemap);
%         clear 'dSensemapFctr';

%         dSensemap = complex(zeros(size(kSpaceL),'single'),zeros(size(kSpaceL),'single'));
%         for t=1:nTime
%             dSensemap(:,:,:,:,t) = bart('ecalib -c0. -I -S -m1', kSpaceL(:,:,:,:,t)); % y-x-z-cha-t
%         end

        lSumMask = logical(sum(obj.fullMask,5));
%         dImgIn = ifftnshift(kSpaceL,1:3) .* dSensemap;

        % registration pairs
%         nPairs = nchoosek(nTime,2);
        [p, q] = meshgrid(1:nTime,1:nTime);
        mask   = triu(ones(nTime), 1) > 0.5;
        pairs = [p(mask) q(mask)];

        % multiresolution approach
        nRes = 5;
        alpha = 0.5;
        beta = [linspace(0.1,0.65,nRes-1),1];
        smallestRes = 64; % in readout direction
        iSize = size(kSpaceL); % y-x-z-cha-t
        iResolutions = [linspace(smallestRes, iSize(1), nRes).', linspace(round(iSize(2)/iSize(1) * smallestRes), iSize(2), nRes).',  linspace(round(iSize(3)/iSize(1) * smallestRes), iSize(3), nRes).'];
        % init
        u = cell(size(pairs,1),2);
        for iI=1:size(u,1)
            u{iI,1} = cell(1,3);
            [u{iI,1}{:}] = deal(zeros(iResolutions(1,:),'single'));
            u{iI,2} = pairs(iI,:);
        end
        matlabpool('open', min([12, size(pairs,1)]));

        kSpaceRes = kSpaceL;
        for iRes=1:nRes
            % first iteration -> downsample from full
            % remaining -> upsample from low_res (with linear interpolation)
            % last iteration -> no downsampling
            kSpaceRes = updownsample(kSpaceRes, false, iResolutions(iRes,:));
            dSenseRes = ones(size(kSpaceRes));
%             dSenseRes = updownsample(dSensemap, false, iResolutions(iRes,:)); % always downsampling from full
            lSumMaskRes = updownsample(lSumMask, false, iResolutions(iRes,:), false); % always downsampling from full, no AA filter
            dImgRes = ifftnshift(kSpaceRes,1:3) .* conj(dSenseRes);
%             dImgDiff = ifftnshift(fftnshift(dImgRes .* conj(dSenseRes),1:3).*obj.fullMask - kSpaceL,1:3) .* dSensemap;
            dImgDiff = ifftnshift(updownsample(updownsample(kSpaceRes,false, size(kSpaceL)).*obj.fullMask - kSpaceL,false,iResolutions(iRes,:)),1:3) .* conj(dSenseRes);

            % LAP registration
            if(iRes > 1) % upsample motion field from last iteration
                [X,Y,Z] = ndgrid(linspace(1,iResolutions(iRes,1),iResolutions(iRes-1,1)), linspace(1,iResolutions(iRes,2),iResolutions(iRes-1,2)), linspace(1,iResolutions(iRes,3),iResolutions(iRes-1,3)));
                [Xq,Yq,Zq] = ndgrid(1:iResolutions(iRes,1),1:iResolutions(iRes,2),   1:iResolutions(iRes,3));
                for iI=1:size(u,1)
                    for iInner=1:3
                        u{iI,1}{iInner} = interpn(X,Y,Z,u{iI,1}{iInner},Xq,Yq,Zq);
                    end
                end
            end
            u = fJointMCCS(squeeze(abs(sum(dImgRes,4))), u, pairs);

            % image update
            for t=1:nTime
                fprintf('Image %d/%d update', t, nTime);
                dSumImg = complex(zeros(size(dImgRes,1),size(dImgRes,2),size(dImgRes,3),'single'),zeros(size(dImgRes,1),size(dImgRes,2),size(dImgRes,3),'single'));
                tmap = 1:nTime; tmap = tmap(~ismember(tmap,t));
                for tI=1:length(tmap)
                    fprintf('.');
                    tInner = tmap(tI);
                    lInd = [tInner, t];
                    iRow = find(ismember(pairs,lInd, 'rows'));
                    if(~isempty(iRow))
                        uDeform = u{iRow,1};
                    else
                        iRow = find(ismember(pairs,lInd(end:-1:1), 'rows'));
                        uDeform = u{iRow,1}; %#ok<FNDSB>
                        uDeform = cellfun(@(x) -x, uDeform, 'UniformOutput', false);
                    end
                    dSumImg = dSumImg + imshift_3D(sum(dImgRes(:,:,:,:,tInner),4),uDeform, 'cubicOMOMS');
                end
                dImgRes(:,:,:,:,t) = dImgRes(:,:,:,:,t) + alpha .* ((1-beta(iRes)) .* dImgDiff(:,:,:,:,t) + beta(iRes)./nTime .* repmat(dSumImg,[1 1 1 size(dSenseRes,4)]) .* dSenseRes(:,:,:,:,t));
                fprintf('\n');
            end

            if(iRes == nRes)
                % CS recon
                kSpaceRes = permute(fftnshift(dImgRes,1:3).* repmat(lSumMaskRes,[1 1 1 1 nTime]),[1 2 3 4 6 7 8 9 10 11 5]) ;
                maps = bart(sprintf('ecalib -c0. -m%d',obj.n_maps), kSpaceRes);
                bartcmd = sprintf('pics -d5 -S -R W:7:0:%f',obj.lambda);
                if(obj.lambdaMC > 0)
                    bartcmd = [bartcmd, sprintf(' -R T:1024:0:%f', obj.lambdaMC)];
                %         bartcmd = [bartcmd, sprintf(' -R L:1024:0:%f', obj.lambdaMC)];
                end
                if(obj.lambdaTV > 0)
                    bartcmd = [bartcmd, sprintf(' -R T:7:0:%f', obj.lambdaTV)];
                end
                imgOut = bart(bartcmd, kSpaceRes, maps);
%                 kSpaceRes = fftnshift(squeeze(repmat(imgOut,[1 1 1 size(dSenseRes,4) 1 1 1 1 1 1 1])) .* dSenseRes,1:3); % y-x-z-cha-t
            else
                kSpaceRes = fftnshift(dImgRes .* dSenseRes,1:3);
            end
        end
        matlabpool('close');
    else
        kSpaceL = permute(kSpaceL,[1 2 3 4 6 7 8 9 10 11 5]);

        maps = bart(sprintf('ecalib -c0. -m%d',obj.n_maps), kSpaceL);

        bartcmd = sprintf('pics -d 5 -R W:7:0:%f',obj.lambda);
        if(obj.lambdaMC > 0)
            bartcmd = [bartcmd, sprintf(' -R T:1024:1024:%f -u 0.1', obj.lambdaMC)];
    %         bartcmd = [bartcmd, sprintf(' -R L:1024:1024:%f', obj.lambdaMC)];
        end
        if(obj.lambdaTV > 0)
            bartcmd = [bartcmd, sprintf(' -R T:7:0:%f -u 0.25', obj.lambdaTV)];
        end
        imgOut = bart(bartcmd, kSpaceL, maps);
    end
    imgOut = permute(imgOut,[2 1 3 11 5 4 6 7 8 9 10]);

    imageCha{1} = imgOut;
%     imageCha = shiftdim(mat2cell(imgOut,size(imgOut,1),size(imgOut,2),size(imgOut,3),size(imgOut,4),ones(1,size(imgOut,5))),3);
end

end

author: thomaskuestner


text:
function ival = r8mat_is_adjacency ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_IS_ADJACENCY checks whether A is an adjacency matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the row and column dimensions of
%    the matrix.  M and N must be positive.
%
%    Input, real A(M,N), the matrix.
%
%    Output, integer IVAL:
%    -1, the matrix is not an adjacency matrix.
%    1, the matrix is an adjacency matrix.
%
  tol = 0.00001;

  if ( m ~= n )
    ival = -1;
    return
  end

  error_frobenius = r8mat_is_symmetric ( m, n, a );

  if ( tol < error_frobenius )
    ival = -1;
    return
  end

  jval = r8mat_is_zero_one ( m, n, a );

  if ( jval ~= 1 )
    ival = -1;
    return
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% .............................................................
    function [shape,naturalDerivatives]=shapeFunctionL2(xi)

    % shape function and derivatives for L2 elements
    % shape : Shape functions
    % naturalDerivatives: derivatives w.r.t. xi
    % xi: natural coordinates (-1 ... +1)

    shape=([1-xi,1+xi]/2);

    naturalDerivatives=[-1,1]/2;

    end % end function shapeFunctionL2

repo: Aero-matlab


text:
% this program validates NN_GradientFunction of a NN from the toolbox

%% Set up network
load house_dataset
net = newff(houseInputs,houseTargets,20);
net = train(net,houseInputs,houseTargets);
% and testing parameters.
TestIndex=110;
TestRange=linspace(-8,8,500); %should be scaled appropriate for data

% Create Gradient network.
[net_ActivateDeriv,OutputOffset] = NN_GradientFunction(net);

% functions used
Prime=@(X,Y) [(-Y(3)+4*Y(2)-3*Y(1))./(-X(3)+4*X(2)-3*X(1)),...
    (Y(3:end)-Y(1:end-2))./(X(3:end)-X(1:end-2)),...
	(3*Y(end)-4*Y(end-1)+Y(end-2))./(3*X(end)-4*X(end-1)+X(end-2))];
% Numerical center difference approximation, implies evenly spaced samples,
% but not structurally required

%% Set up truth comparison, single dimension

Base     =zeros(size(houseInputs,1),length(TestRange));
Base_Diff=zeros(size(houseInputs,1),length(TestRange));
ActivateDeriv=zeros(size(houseInputs,1),length(TestRange));
for alley=3:6%1:size(houseInputs,1); %abbreviated to not clutter document
% define the test points along the dimension specified by alley
TestPoints=repmat(houseInputs(:,TestIndex),1,length(TestRange));
TestPoints(alley,:)=TestPoints(alley,:)+TestRange(:).';

temp = sim(net_ActivateDeriv,TestPoints)-repmat(OutputOffset,1,size(TestPoints,2));
ActivateDeriv(alley,:)=temp(alley,:);

Base(alley,:)= sim(net,TestPoints);
Base_Diff(alley,:)=Prime(TestRange,Base(alley,:));

% Now Show that you did it right
figure(20+100*alley)
subplot(311)
plot(TestPoints(alley,:),Base(alley,:),'g.')
title(['Original Points, ',num2str(alley), 'th dim'] )
subplot(312)
plot(TestPoints(alley,:),Base_Diff(alley,:),'b.')
title('Numerically Differentiated Points')
subplot(313)
plot(TestPoints(alley,:),ActivateDeriv(alley,:),'r.')
title('Neural Network Differentiated Points')
figure(1+100*alley)
b=Base_Diff(alley,:)-ActivateDeriv(alley,:);
plot(b)
title(['Difference from two methods, ',num2str(alley), 'th dim'] )
end
%% Activation function derivative check, symbolic
% x_s=sym('x_s');
% y_s= 2/(1+exp(-2*x_s))-1; %Eqn of hyperbolic tangent, from apply_transfer
% dy_s=diff(y_s,x_s); % Put into apply_transfer of modified file
% ddy_s=diff(dy_s,x_s); % Put into derivative of modified file
% fprintf('Derivative expression:\n');
% pretty(dy_s);
% fprintf('Derivative expression, simplified:\n');
% pretty(simple(dy_s));
% fprintf('Second derivative expression:\n');
% pretty(ddy_s);
% fprintf('Second derivative expression, simplified:\n');
% pretty(simple(ddy_s));
%% Activation function derivative numerical check
% n1 = linspace(-5,5,1000);
% a1 = tansig(n1); % replace function with desired
% figure(165)
% subplot(311);
% plot(n1,a1)
% title('Original activation function');
% da1=Prime(n1,a1);
% da2=dtansig_0(n1);
% subplot(323);
% plot(n1,da1,'b.',n1,da2,'ro')
% title('First Derivative')
% legend('Numeric','By Function');
% subplot(324)
% plot(n1,da1-da2,'k.')
% title('Difference of 1st derivative')
% dda1=Prime(n1,da1);
% dda2=dtansig_0('dn',n1);
% subplot(325);
% plot(n1,dda1,'b.',n1,dda2,'ro')
% title('Second Derivative')
% legend('Numeric','By Function');
% subplot(326)
% plot(n1,dda1-dda2,'k.')
% title('Difference of 2nd derivative')
% % Note that magnitude of the difference is primarily due to numerical
% % approximation error. Decreasing the step size of n1 will decrease the
% % magnitude of the error.

author: Sable

repo: mcbench-benchmarks


text:
function y = cvx_constant( x )
error( nargchk( 1, 1, nargin ) );
y = x;

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function submitWeb(partId)
%SUBMITWEB Generates a base64 encoded string for web-based submissions
%   SUBMITWEB() will generate a base64 encoded string so that you can submit your
%   solutions via a web form

  fprintf('==\n== [ml-class] Submitting Solutions | Programming Exercise %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid homework part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames));
    fprintf('!! Submission Cancelled\n');
    return
  end

  [login] = loginPrompt();
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  [result] = submitSolution(login, partId, output(partId), ...
                            source(partId));
  result = base64encode(result);

  fprintf('\nSave as submission file [submit_ex%s_part%d.txt]: ', ...
          homework_id(), partId);
  saveAsFile = input('', 's');
  if (isempty(saveAsFile))
    saveAsFile = sprintf('submit_ex%s_part%d.txt', homework_id(), partId);
  end

  fid = fopen(saveAsFile, 'w');
  if (fid)
    fwrite(fid, result);
    fclose(fid);
    fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
    fprintf(['You can now submit your solutions through the web \n' ...
             'form in the programming exercises. Select the corresponding \n' ...
             'programming exercise to access the form.\n']);

  else
    fprintf('Unable to save to %s\n\n', saveAsFile);
    fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
    pause;
    fprintf(result);
  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '7';
end

function [partNames] = validParts()
  partNames = {
                'Find Closest Centroids (k-Means)', ...
                'Compute Centroid Means (k-Means)' ...
                'PCA', ...
                'Project Data (PCA)', ...
                'Recover Data (PCA)' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'findClosestCentroids.m' }, ...
           { 'computeCentroids.m' }, ...
           { 'pca.m' }, ...
           { 'projectData.m' }, ...
           { 'recoverData.m' } ...
           };
end

function out = output(partId)
  % Random Test Cases
  X = reshape(sin(1:165), 15, 11);
  Z = reshape(cos(1:121), 11, 11);
  C = Z(1:5, :);
  idx = (1 + mod(1:15, 3))';
  if partId == 1
    idx = findClosestCentroids(X, C);
    out = sprintf('%0.5f ', idx(:));
  elseif partId == 2
    centroids = computeCentroids(X, idx, 3);
    out = sprintf('%0.5f ', centroids(:));
  elseif partId == 3
    [U, S] = pca(X);
    out = sprintf('%0.5f ', abs([U(:); S(:)]));
  elseif partId == 4
    X_proj = projectData(X, Z, 5);
    out = sprintf('%0.5f ', X_proj(:));
  elseif partId == 5
    X_rec = recoverData(X(:,1:5), Z, 5);
    out = sprintf('%0.5f ', X_rec(:));
  end
end

% ========================= SUBMIT HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          while ~feof(fid)
            line = fgets(fid);
            src = [src line];
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames));
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('\nEnter your choice [1-%d]: ', ...
          numel(partNames));
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [result, str] = submitSolution(email, part, output, source)

  result = ['a:5:{' ...
            p_s('homework') p_s64(homework_id()) ...
            p_s('part') p_s64(part) ...
            p_s('email') p_s64(email) ...
            p_s('output') p_s64(output) ...
            p_s('source') p_s64(source) ...
            '}'];

end

function s = p_s(str)
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

function s = p_s64(str)
   str = base64encode(str, '');
   s = ['s:' num2str(numel(str)) ':"' str '";'];
end

% =========================== LOGIN HELPERS ===========================

function [login] = loginPrompt()
  % Prompt for password
  [login] = basicPrompt();
end

function [login] = basicPrompt()
  login = input('Login (Email address): ', 's');
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function [flowx,flowy] = specific_wind(x,y,nel)
%circular_wind   Reference problem 3.4 convective wind
%   [flowx,flowy] = specific_wind(x,y,nel);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          nel        number of elements
%
%   specifies circular wind /Morton pp.10/
%   IFISS function: DJS; 5 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      flowx =  2*y.*(1-x.*x);
      flowy = -2*x.*(1-y.*y);
      return


text:
function [fun,lb,ub,nlcon,cl,cu,xtype] = GetBARON(B)
%Build BARON arguments from SymBuilder Object
%
%   Called By SYMBUILDER Class

%   Copyright (C) 2012 Jonathan Currie (I2C2)

%Check object is built
if(~IsBuilt(B))
    error('Please build the object using Draft() or Build()');
end

%Check we have an objective (no good without!)
if(B.noObjs == 0)
    error('You can only get BARON arguments when you have specified an objective!');
elseif(B.noObjs > 1)
    error('You can only get BARON arguments with a single objective');
end

%Clear existing functions from memory
clear symb_grad symb_hess symb_nlcon symb_nljac symb_obj

%Objective
fun = GetAllObj(B);
%Bounds
lb = B.lb;
ub = B.ub;
%Nonlinear Constraints
opts.preallocate = false;
[nlcon,cl,cu] = GetAllCon(B,opts);
%Integer Constraints
xtype = B.xtype;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%%% show_2dvideo(2dmatrix,int,int)
%
function show_2dvideo(M,m,n)
  for i = 1 : size(M,2)
    I = reshape(M(:,i),m,n);
    imshow(I,[],'InitialMagnification','fit');
    %imshow(I,'InitialMagnification','fit');
    disp(i);
    pause(0.01);
  end
end

author: andrewssobral


text:
function [V,gradV,HessianV]=spherHarmonicSetEval(C,S,point,a,c,systemType,spherDerivs,scalFactor,algorithm)
%%SPHERHARMONICSETEVAL Evaluate a set of real or complex potentials,
%                   gradients and Hessians when the potentials are
%                   expressed in terms of a set of real or complex
%                   spherical harmonic coefficients. This function is very
%                   similar to the spherHarmonicEval function, but it
%                   allows C and S to be matrices, where each column
%                   contains a set of spherical harmonic coefficients. The
%                   evaluation of multiple sets at once is faster than
%                   making multiple calls to spherHarmonicEval. The
%                   ordering of the elements of the outputs of this
%                   function is different than that used in
%                   spherHarmonicEval.
%
%INPUTS: C An ((M+2)*(M+1)/2)XnumSets matrix holding numSets real or
%          complex coefficient terms that are multiplied by cosines in the
%          harmonic expansion. The coefficients must be fully normalized
%          using the type of full normalization that is used in the EGM2008
%          model. Their normalization type can be changed using the
%          changeSpherHarmonicNorm function If given to a
%          CountingClusterSet class, C(n+1,m+1) is the 1XnumSets vector of
%          the coefficient of degree n and order m for each set. When a
%          maximum degree of M is used, all C must have values for all n
%          from 0 to M and for all m from 0 to n for each n. If
%          coefficients are not present for certain degrees, then insert a
%          0. It is assumed that M>=3.
%        S An ((M+2)*(M+1)/2)XnumSets matrix holding numSets sets of
%          coefficient terms that are multiplied by sines in the harmonic
%          expansion. The requirements on S are the same as those on C.
%    point The 3XN set of N real points at which the potential and/or
%          gradient should be evaluated given in SPHERICAL coordinates
%          consisting of [r;azimuth;elevation]; When evaluating points on a
%          grid, the algorithm will be fastest if the points are provided
%          presorted by range and then by azimuth. This reduces the amount
%          of recomputation of certain values. Alternatively, if C and S
%          are for evaluating terrain heights, then points are 2XN having
%          the format [azimuth;elevation] and it is best if the points are
%          sorted by azimuth. Azimuth is measured counterclockwise from the
%          x-axis in the x-y plane. Elevation is measured up from the x-y
%          plane (towards the z-axis). Angles must be given in radians.
%        a The real or complex numerator in the (a/r)^n term in the
%          spherical harmonic sum. Normally, this is some type of a
%          reference radius. For example, when using most gravitational
%          models, a is the semi-major axis of the reference ellipsoid. If
%          this parameter is omitted, it is assumed that one is using the
%          spherical harmonics with something like the National Geospatial
%          Intelligence Agency's (NGA's) EGM96 or EGM2008 models, in which
%          case a=Constants.EGM2008SemiMajorAxis is used unless point is
%          2D, in which case c=1 is used.
%        c The real or complex constant value by which the spherical
%          harmonic series is multiplied. For example, for gravitational
%          potentials, the value is usually GM where G is the universal
%          gravitational constant and  M is the mass of the Earth. When
%          using the International Geomagnetic Reference Field (IGRF), the
%          constant is a^2, where a is the same as the numerator in the a/r
%          term. If this parameter is omitted, it is assumed that one is
%          using the spherical harmonics with something like the NGA's
%          EGM96 or EGM2008 models, in which case c=Constants.EGM2008GM is
%          used unless point is 2D, in which case c=1 is used.
% systemType An optional parameter specifying the axis from which the
%           angles are measured in radians. Possible values are
%           0 (The default if omitted) Azimuth is measured
%             counterclockwise from the x-axis in the x-y plane. Elevation
%             is measured up from the x-y plane (towards the z-axis). This
%             is consistent with common spherical coordinate systems for
%             specifying longitude (azimuth) and geocentric latitude
%             (elevation).
%           2 This is the same as 0 except instead of being given
%             elevation, one is given the angle away from the z-axis, which
%             is (pi/2-elevation).
% spherDerivs This parameter specifies whether the gradient and Hessian
%          terms should be returned in spherical coordinates rather than
%          Cartesian coordinates. Possible values are true and false. The
%          default if this parameter is omitted or an empty matrix is
%          passed is false. The spherical coordinate system used by the
%          gradient and Hessian will match that specified by the systemType
%          input.
% scalFactor An optional real scale factor used in computing the normalized
%          associated Legendre polynomials. Generally, the default value
%          (if the scalFactor parameter is omitted) of 10^(-280) that is
%          suggested in [1] is sufficient. When very high-order models are
%          used, this scale factor prevents overflows.
% algorithm An optional parameter that selects which algorithm to use.
%          Possible values are:
%          0 (The default if omitted or an empty matrix is passed) If only
%            the potential is desired, always use Legendre's algorithm.
%            Otherwise, for points above 88 degrees spherical latitude, use
%            Pines' algorithm and for all other points use Legendre's
%            algorithm.
%          1 Only use Legendre's algorithm regardless of the location of
%            the points. Note that the gradient is singular at the poles,
%            so numerical problems will arise.
%          2 Only use Pines' algorithm.
%
%OUTPUTS: V The numSetsXN scalar potentials as obtained from the spherical
%           harmonic series for each set and point. This is real is all of
%           the inputs are real.
%     gradV The 3XnumSetsXN set of gradients of the potential in Cartesian
%           coordinates for each set and point as obtained using the
%           spherical harmonic series. The derivatives are in the order
%           [dV/dx;dV/dy;dV/dz] for Cartesian values and in the order
%           [dV/dr;dV/dAz;dV/dEl] for spherical values.
%     HessV The 3X3XnumSetsXN collection of Hessian matrices of the
%           potential for each set and point.  If Cartesian derivatives are
%           used, then the ordering is
%           [d2/(dxdx),d2/(dxdy),d2/(dxdz);
%            d2/(dydx),d2/(dydy),d2/(dydz);
%            d2/(dzdx),d2/(dzdy),d2/(dxdx)]; If spherical derivatives are
%            used, then the ordering is the same with (x,y,z) replaced by
%            (r,Az,El).
%
%This function is useful for the evaluation of scalar spherical harmonic
%series with vector coefficients, as can be used for modelling 3D far-field
%antenna response patterns as in [1].
%
%This function is implemented in the same manner as spherHarmonicEval,
%except minor changes to allow C and S to be matrices rather than vector
%have been made. Additionally, for numSet=1, the shape of the matrices
%outputted by this function differs from that outputted by
%spherHarmonicEval.
%
%It is recommended that the C++ version of this function be compiled and
%used in place of this as it is thousands of times faster when C and S are
%very large.
%
%REFERENCES:
%[1] J. Rahola, F. Belloni, and A. Richter, "Modelling of radiation
%    patterns using scalar spherical harmonics with vector coefficients,"
%    in Proceedings of the 3rd European Conference on Antennas and
%    Propagation, Berlin, Germany, 23-27 Mar. 2009, pp. 3361-3365.
%
%December 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<9||isempty(algorithm))
    algorithm=0;
end

if(nargin<8||isempty(scalFactor))
    scalFactor=10^(-280);
end

if(nargin<7||isempty(spherDerivs))
    spherDerivs=false;
end

if(nargin<6||isempty(systemType))
    systemType=0;
end

if(nargin<5||isempty(c))
    if(size(point,1)==2)
        c=1;
    else
        c=Constants.EGM2008GM;
    end
end

if(nargin<4||isempty(a))
    if(size(point,1)==2)
        a=1;
    else
        a=Constants.EGM2008SemiMajorAxis;
    end
end

if(systemType~=0&&systemType~=2)
    error('An unsupported systemType was specified.')
end

M=(1/2)*(sqrt(1+8*length(C))-1)-1;

if(M<3)
    error('The coefficients must be provided to at least degree 3. To use a lower degree, one can insert zero coefficients.');
end

numPoints=size(point,2);
%If we are evaluating terrain heights.
switch(size(point,1))
    case 2
        point=[ones(1,numPoints);point(1,:);point(2,:)];
    case 3
    otherwise
        error('Invalid point length');
end

%Using a CountingClusterSet simplifies the indexation of the coefficients.
C=CountingClusterSet(C);
S=CountingClusterSet(S);
numSets=C.numSets();

%Preallocate space used by the modified forward row algorithm when
%evaluating over multiple values with the same range and latitude but
%different longitudes.
if(algorithm==0||algorithm==1)
    A=zeros(numSets,M+1);
    B=zeros(numSets,M+1);
    if(nargout>1)
        Ar=zeros(numSets,M+1);
        Br=zeros(numSets,M+1);
        ATheta=zeros(numSets,M+1);
        BTheta=zeros(numSets,M+1);

        if(nargout>2)
            AThetaTheta=zeros(numSets,M+1);
            BThetaTheta=zeros(numSets,M+1);

            Arr=zeros(numSets,M+1);
            Brr=zeros(numSets,M+1);

            AThetar=zeros(numSets,M+1);
            BThetar=zeros(numSets,M+1);
        end
    end
end

V=zeros(numSets,numPoints);
gradV=zeros(3,numSets,numPoints);
HessianV=zeros(3,3,numSets,numPoints);

rPrev=Inf;
thetaPrev=Inf;
for curPoint=1:numPoints
    pointCur=point(:,curPoint);

    if(systemType==2)
        pointCur(3)=pi/2-pointCur(3);
    end

    r=pointCur(1);
    lambda=pointCur(2);
    thetaCur=pointCur(3);

    rChanged=r~=rPrev;
    thetaChanged=thetaCur~=thetaPrev;
    rPrev=r;
    thetaPrev=thetaCur;

    if(rChanged)
        crScal=(c/r)/scalFactor;

        %This stores all of the powers of a/r needed for the sum,
        %regardless of which algorithm is used.
        nCoeff=zeros(M+1,1);
        nCoeff(1)=1;
        for n=1:M
            nCoeff(n+1)=nCoeff(n)*(a/r);
        end
    end

    switch(algorithm)
        case 0
        %At latitudes that are not near the poles, the Legendre method is
        %used. It cannot be used for the gradient or Hessian near the
        %poles, because of the singularity of the spherical coordinate
        %system.
            useLegendre=abs(thetaCur)<88*(pi/180)||nargout<2;
        case 1
            useLegendre=true;
        case 2
            useLegendre=false;
        otherwise
            error('Unknown Algorithm option specified.');
    end

    if(useLegendre)
        %Compute the sine and cosine terms.
        [SinVec,CosVec]=calcSinCosTerms(lambda,M);
        %The formulae for spherical harmonic synthesis with Legendre's
        %method uses clolatitude, pi/2-elevation
        theta=pi/2-thetaCur;
        u=sin(theta);
        if(thetaChanged)
            if(nargout==2)
                [PBarUVals,dPBarUValsdTheta]=NALegendreCosRat(theta,M,scalFactor);
            elseif(nargout>2)
                [PBarUVals,dPBarUValsdTheta,d2PBarUValsdTheta2]=NALegendreCosRat(theta,M,scalFactor);
            else
                PBarUVals=NALegendreCosRat(theta,M,scalFactor);
            end
        end

        %Evaluate Equation 7 from the Holmes and Featherstone paper.
        if(rChanged||thetaChanged)
            A(:)=0;
            B(:)=0;
            for m=0:M
                for n=m:M
                    CScal=nCoeff(n+1)*C(n+1,m+1).';
                    SScal=nCoeff(n+1)*S(n+1,m+1).';

                    %From Table 4
                    A(:,m+1)=A(:,m+1)+CScal*PBarUVals(n+1,m+1);
                    B(:,m+1)=B(:,m+1)+SScal*PBarUVals(n+1,m+1);
                end
            end

            %If additional terms should be computed so a gradient can
            %be determined.
            if(nargout>1)
                Ar(:)=0;
                Br(:)=0;
                ATheta(:)=0;
                BTheta(:)=0;
                for m=0:M
                    for n=m:M
                        CScal=nCoeff(n+1)*C(n+1,m+1).';
                        SScal=nCoeff(n+1)*S(n+1,m+1).';

                        %From Table 4
                        Ar(:,m+1)=Ar(:,m+1)+(n+1)*CScal*PBarUVals(n+1,m+1);
                        Br(:,m+1)=Br(:,m+1)+(n+1)*SScal*PBarUVals(n+1,m+1);

                        %From Table 4
                        ATheta(:,m+1)=ATheta(:,m+1)+CScal*dPBarUValsdTheta(n+1,m+1);
                        BTheta(:,m+1)=BTheta(:,m+1)+SScal*dPBarUValsdTheta(n+1,m+1);
                    end
                end

                %If additional terms should be computed so a Hessian can
                %be determined.
                if(nargout>2)
                    Arr(:)=0;
                    Brr(:)=0;
                    AThetar(:)=0;
                    BThetar(:)=0;
                    AThetaTheta(:)=0;
                    BThetaTheta(:)=0;
                    for m=0:M
                        for n=m:M
                            CScal=nCoeff(n+1)*C(n+1,m+1).';
                            SScal=nCoeff(n+1)*S(n+1,m+1).';

                            %From Table 5, with the correction from the
                            %erratum.
                            Arr(:,m+1)=Arr(:,m+1)+(n+1)*(n+2)*CScal*PBarUVals(n+1,m+1);
                            Brr(:,m+1)=Brr(:,m+1)+(n+1)*(n+2)*SScal*PBarUVals(n+1,m+1);

                            %From Table 5
                            AThetar(:,m+1)=AThetar(:,m+1)+(n+1)*CScal*dPBarUValsdTheta(n+1,m+1);
                            BThetar(:,m+1)=BThetar(:,m+1)+(n+1)*SScal*dPBarUValsdTheta(n+1,m+1);

                            %From Table 5
                            AThetaTheta(:,m+1)=AThetaTheta(:,m+1)+CScal*d2PBarUValsdTheta2(n+1,m+1);
                            BThetaTheta(:,m+1)=BThetaTheta(:,m+1)+SScal*d2PBarUValsdTheta2(n+1,m+1);
                        end
                    end
                end
            end
        end

        %Use Horner's method to compute V all values.
        V(:,curPoint)=0;
        for m=M:-1:0
            V(:,curPoint)=V(:,curPoint)*u+(A(:,m+1)*CosVec(m+1)+B(:,m+1)*SinVec(m+1));
        end

        V(:,curPoint)=crScal*V(:,curPoint);

        if(nargout>1)
            %Use Horner's method to compute all dV values.
            dVdr=zeros(numSets,1);
            dVdLambda=zeros(numSets,1);
            dVdTheta=zeros(numSets,1);

            %The following first-order derivative formulae are from
            %Table 1 (expressed using Horner's method).
            for m=M:-1:0
                dVdr=dVdr*u-(Ar(:,m+1)*CosVec(m+1)+Br(:,m+1)*SinVec(m+1));
                dVdLambda=dVdLambda*u-m*(A(:,m+1)*SinVec(m+1)-B(:,m+1)*CosVec(m+1));
                dVdTheta=dVdTheta*u+(ATheta(:,m+1)*CosVec(m+1)+BTheta(:,m+1)*SinVec(m+1));
            end

            dVdr=crScal*dVdr.'/r;
            dVdLambda=crScal*dVdLambda.';
            %The minus sign adjusts for the coordinate system change.
            dVdTheta=-crScal*dVdTheta.';

            if(spherDerivs)
                %The sign change on the theta terms deals with the
                %different input coordinate system used.
                gradV(:,:,curPoint)=[dVdr;dVdLambda;dVdTheta];
            else%Convert the derivatives to Cartesian coordinates.
                J=calcSpherConvJacob(pointCur);

                for curSet=1:numSets
                    gradV(:,curSet,curPoint)=J'*[dVdr(curSet);dVdLambda(curSet);dVdTheta(curSet)];
                end
            end
        end

        if(nargout>2)
            %Use Horner's method to compute d2V all values.
            d2VdLambdadLambda=zeros(numSets,1);
            d2VdLambdadTheta=zeros(numSets,1);
            d2VdrdLambda=zeros(numSets,1);
            d2VdThetadTheta=zeros(numSets,1);
            d2VdrdTheta=zeros(numSets,1);
            d2Vdrdr=zeros(numSets,1);

            %The following second-order derivative formulae are from
            %Table 2 (expressed using Horner's method).
            for m=M:-1:0
                d2VdLambdadLambda=d2VdLambdadLambda*u-m^2*(A(:,m+1)*CosVec(m+1)+B(:,m+1)*SinVec(m+1));
                d2VdLambdadTheta=d2VdLambdadTheta*u+m*(ATheta(:,m+1)*SinVec(m+1)-BTheta(:,m+1)*CosVec(m+1));
                d2VdrdLambda=d2VdrdLambda*u+m*(Ar(:,m+1)*SinVec(m+1)-Br(:,m+1)*CosVec(m+1));
                d2VdThetadTheta=d2VdThetadTheta*u+(AThetaTheta(:,m+1)*CosVec(m+1)+BThetaTheta(:,m+1)*SinVec(m+1));
                d2VdrdTheta=d2VdrdTheta*u-(AThetar(:,m+1)*CosVec(m+1)+BThetar(:,m+1)*SinVec(m+1));
                d2Vdrdr=d2Vdrdr*u+(Arr(:,m+1)*CosVec(m+1)+Brr(:,m+1)*SinVec(m+1));
            end

            %The minus signs in the following equations adjust for the
            %spherical coordinate system difference.
            d2Vdrdr=crScal*d2Vdrdr.'/r^2;
            d2VdLambdadLambda=crScal*d2VdLambdadLambda.';
            d2VdThetadTheta=crScal*d2VdThetadTheta.';
            d2VdrdLambda=crScal*d2VdrdLambda.'/r;
            d2VdrdTheta=-crScal*d2VdrdTheta.'/r;
            d2VdLambdadTheta=crScal*d2VdLambdadTheta.';

            if(spherDerivs)
                HessianV(1,1,:,curPoint)=d2Vdrdr;
                HessianV(2,2,:,curPoint)=d2VdLambdadLambda;
                HessianV(3,3,:,curPoint)=d2VdThetadTheta;
                HessianV(1,2,:,curPoint)=d2VdrdLambda;
                HessianV(2,1,:,curPoint)=HessianV(1,2,:,curPoint);
                HessianV(1,3,:,curPoint)=d2VdrdTheta;
                HessianV(3,1,:,curPoint)=HessianV(1,3,:,curPoint);
                HessianV(2,3,:,curPoint)=d2VdLambdadTheta;
                HessianV(3,2,:,curPoint)=HessianV(2,3,:,curPoint);
            else%Convert the Hessian to Cartesian coordinates.
                drdx=J(1,1);
                drdy=J(1,2);
                drdz=J(1,3);
                dLambdadx=J(2,1);
                dLambdady=J(2,2);
                dLambdadz=J(2,3);
                dPhidx=J(3,1);
                dPhidy=J(3,2);
                dPhidz=J(3,3);

                H=calcSpherConvHessian(pointCur);

                drdxdx=H(1,1,1);
                drdydy=H(2,2,1);
                drdzdz=H(3,3,1);
                drdxdy=H(1,2,1);
                drdxdz=H(1,3,1);
                drdydz=H(2,3,1);

                dLambdadxdx=H(1,1,2);
                dLambdadydy=H(2,2,2);
                dLambdadzdz=H(3,3,2);
                dLambdadxdy=H(1,2,2);
                dLambdadxdz=H(1,3,2);
                dLambdadydz=H(2,3,2);

                dPhidxdx=H(1,1,3);
                dPhidydy=H(2,2,3);
                dPhidzdz=H(3,3,3);
                dPhidxdy=H(1,2,3);
                dPhidxdz=H(1,3,3);
                dPhidydz=H(2,3,3);

                HessianV(1,1,:,curPoint)=d2VdLambdadLambda*dLambdadx^2+2*d2VdLambdadTheta*dLambdadx*dPhidx+d2VdThetadTheta*dPhidx^2+2*dPhidx*drdx*d2VdrdTheta+2*dLambdadx*drdx*d2VdrdLambda+dVdLambda*dLambdadxdx+dVdTheta*dPhidxdx+dVdr*drdxdx+drdx^2*d2Vdrdr;
                HessianV(2,2,:,curPoint)=d2VdThetadTheta*dPhidy^2+2*dLambdady*dPhidy*d2VdLambdadTheta+dVdLambda*dLambdadydy+dVdTheta*dPhidydy+dLambdady^2*d2VdLambdadLambda+drdydy*dVdr+2*dPhidy*drdy*d2VdrdTheta+2*dLambdady*drdy*d2VdrdLambda+drdy^2*d2Vdrdr;
                HessianV(3,3,:,curPoint)=dVdTheta*dPhidzdz+dPhidz^2*d2VdThetadTheta+dLambdadzdz*dVdLambda+2*dLambdadz*dPhidz*d2VdLambdadTheta+(dLambdadz)^2*d2VdLambdadLambda+drdzdz*dVdr+2*dPhidz*drdz*d2VdrdTheta+2*dLambdadz*drdz*d2VdrdLambda+drdz^2*d2Vdrdr;
                HessianV(1,2,:,curPoint)=dPhidy*d2VdLambdadTheta*dLambdadx+dLambdady*d2VdLambdadLambda*dLambdadx+d2VdThetadTheta*dPhidy*dPhidx+dLambdady*d2VdLambdadTheta*dPhidx+drdy*dPhidx*d2VdrdTheta+dPhidy*drdx*d2VdrdTheta+dVdLambda*dLambdadxdy+dVdTheta*dPhidxdy+dVdr*drdxdy+drdy*dLambdadx*d2VdrdLambda+dLambdady*drdx*d2VdrdLambda+drdy*drdx*d2Vdrdr;
                HessianV(2,1,:,curPoint)=HessianV(1,2,:,curPoint);
                HessianV(1,3,:,curPoint)=dPhidz*d2VdLambdadTheta*dLambdadx+dLambdadz*d2VdLambdadLambda*dLambdadx+dPhidz*d2VdThetadTheta*dPhidx+dLambdadz*d2VdLambdadTheta*dPhidx+dVdLambda*dLambdadxdz+dVdTheta*dPhidxdz+dVdr*drdxdz+drdz*dPhidx*d2VdrdTheta+dPhidz*drdx*d2VdrdTheta+drdz*dLambdadx*d2VdrdLambda+dLambdadz*drdx*d2VdrdLambda+drdz*drdx*d2Vdrdr;
                HessianV(3,1,:,curPoint)=HessianV(1,3,:,curPoint);
                HessianV(2,3,:,curPoint)=dPhidz*d2VdThetadTheta*dPhidy+dVdLambda*dLambdadydz+dVdTheta*dPhidydz+dPhidz*dLambdady*d2VdLambdadTheta+dLambdadz*dPhidy*d2VdLambdadTheta+dLambdadz*dLambdady*d2VdLambdadLambda+drdydz*dVdr+drdz*dPhidy*d2VdrdTheta+dPhidz*drdy*d2VdrdTheta+drdz*dLambdady*d2VdrdLambda+dLambdadz*drdy*d2VdrdLambda+drdz*drdy*d2Vdrdr;
                HessianV(3,2,:,curPoint)=HessianV(2,3,:,curPoint);
            end
        end
    else
        %At latitudes that are near the poles, the non-singular algorithm
        %of Pines using the fully normalized Helmholtz equations from
        %Fantino and Casotto is used. The algorithm has been slightly
        %modified so that the c/r term is out front and the fully
        %normalized Helmholtz polynomials can be scaled. Also, lumped
        %coefficients are not used. The Pines algorithm is generally slower
        %than the algorithm of Holmes and Featherstone and it suffers a
        %loss of precision near the equator. Thus, the Pines algorithm is
        %only used near the poles where the other algorithm has issues with
        %a singularity.

        CartPoint=spher2Cart(pointCur);

        x=CartPoint(1);
        y=CartPoint(2);
        z=CartPoint(3);

        %Get the direction cosines used by Pines' algorithm.
        s=x/r;
        t=y/r;
        u=z/r;

        %Compute the fully normalized Helmholtz polynomials.
        if(thetaChanged)
            if(nargout>1)
                [HBar,dHBardu,d2HBardu2]=normHelmholtz(u,M,scalFactor);
            else
                HBar=normHelmholtz(u,M,scalFactor);
            end
        end

        %Recursively compute the rm and im terms for the sums.
        rm=zeros(M+1,1);
        im=zeros(M+1,1);
        rm(0+1)=1;
        im(0+1)=0;
        for m=1:M
            %These are equation 49 in the Fantino and Casotto paper.
            rm(m+1)=s*rm(m-1+1)-t*im(m-1+1);
            im(m+1)=s*im(m-1+1)+t*rm(m-1+1);
        end

        %Perform the sum for the potential from Equation 44 in the Fantino
        %and Casotto paper.
        V(:,curPoint)=0;
        for n=0:M
            innerTerm=0;
            for m=0:n
                innerTerm=innerTerm+(C(n+1,m+1)*rm(m+1)+S(n+1,m+1)*im(m+1))*HBar(n+1,m+1);
            end
            V(:,curPoint)=V(:,curPoint)+nCoeff(n+1)*innerTerm.';
        end

        V(:,curPoint)=crScal*V(:,curPoint);

        %If only the gradient is desired as the next output
        if(nargout==2)
            %If the gradient and not the Hessian is desired.
            a1=zeros(1,numSets);
            a2=zeros(1,numSets);
            a3=zeros(1,numSets);
            a4=zeros(1,numSets);

            for m=0:M
                A1=zeros(1,numSets);
                A2=zeros(1,numSets);
                A3=zeros(1,numSets);

                B1=zeros(1,numSets);
                B2=zeros(1,numSets);
                B3=zeros(1,numSets);

                %Compute the lumped coefficients for Pine's method from
                %Table 13 for the current m.
                for n=m:M
                    HVal=HBar(n+1,m+1);
                    dHVal=dHBardu(n+1,m+1);

                    %The expressions for Lmn, is from Table 14
                    Lmn=(n+m+1)*HVal+u*dHVal;

                    rhoC=nCoeff(n+1)*C(n+1,m+1);
                    rhoS=nCoeff(n+1)*S(n+1,m+1);

                    A1=A1+rhoC*HVal;
                    A2=A2+rhoC*dHVal;
                    A3=A3+rhoC*Lmn;

                    B1=B1+rhoS*HVal;
                    B2=B2+rhoS*dHVal;
                    B3=B3+rhoS*Lmn;
                end
                if(m>=1)
                    a1=a1+m*(A1*rm(m-1+1)+B1*im(m-1+1));
                    a2=a2+m*(B1*rm(m-1+1)-A1*im(m-1+1));
                end
                a3=a3+(A2*rm(m+1)+B2*im(m+1));
                a4=a4-(A3*rm(m+1)+B3*im(m+1));
            end

            a1=a1/r;
            a2=a2/r;
            a3=a3/r;
            a4=a4/r;

            dVdx=crScal*(a1+s*a4);
            dVdy=crScal*(a2+t*a4);
            dVdz=crScal*(a3+u*a4);

            if(spherDerivs)
                %Convert the derivatives to spherical coordinates.
                J=calcSpherInvJacob(pointCur)';

                gradV(:,:,curPoint)=J*[dVdx;dVdy;dVdz];
            else%If a gradient in Cartesian coordinates is desired.
                gradV(1,:,curPoint)=dVdx;
                gradV(2,:,curPoint)=dVdy;
                gradV(3,:,curPoint)=dVdz;
            end
        else
            %If the gradient and the Hessian are desired.
            a1=zeros(1,numSets);
            a2=zeros(1,numSets);
            a3=zeros(1,numSets);
            a4=zeros(1,numSets);
            a11=zeros(1,numSets);
            a12=zeros(1,numSets);
            a13=zeros(1,numSets);
            a14=zeros(1,numSets);
            a23=zeros(1,numSets);
            a24=zeros(1,numSets);
            a33=zeros(1,numSets);
            a34=zeros(1,numSets);
            a44=zeros(1,numSets);

            for m=0:M
                A1=zeros(1,numSets);
                A2=zeros(1,numSets);
                A3=zeros(1,numSets);
                A4=zeros(1,numSets);
                A5=zeros(1,numSets);
                A6=zeros(1,numSets);

                B1=zeros(1,numSets);
                B2=zeros(1,numSets);
                B3=zeros(1,numSets);
                B4=zeros(1,numSets);
                B5=zeros(1,numSets);
                B6=zeros(1,numSets);

                %Compute the lumped coefficients for Pine's method from
                %Table 13 for the current m.
                for n=m:M
                    HVal=HBar(n+1,m+1);
                    dHVal=dHBardu(n+1,m+1);
                    d2HVal=d2HBardu2(n+1,m+1);

                    %The expressions for Lmn, dLmn, and Omn are from
                    %Table 14
                    Lmn=(n+m+1)*HVal+u*dHVal;
                    dLmn=(n+m+2)*dHVal+u*d2HVal;
                    Omn=(n+m+1)*(n+m+2)*HVal+2*u*(n+m+2)*dHVal+u^2*d2HVal;

                    rhoC=nCoeff(n+1)*C(n+1,m+1);
                    rhoS=nCoeff(n+1)*S(n+1,m+1);

                    A1=A1+rhoC*HVal;
                    A2=A2+rhoC*dHVal;
                    A3=A3+rhoC*Lmn;
                    A4=A4+rhoC*d2HVal;
                    A5=A5+rhoC*dLmn;
                    A6=A6+rhoC*Omn;

                    B1=B1+rhoS*HVal;
                    B2=B2+rhoS*dHVal;
                    B3=B3+rhoS*Lmn;
                    B4=B4+rhoS*d2HVal;
                    B5=B5+rhoS*dLmn;
                    B6=B6+rhoS*Omn;
                end
                if(m>=1)
                    a1=a1+m*(A1*rm(m-1+1)+B1*im(m-1+1));
                    a2=a2+m*(B1*rm(m-1+1)-A1*im(m-1+1));
                end
                a3=a3+(A2*rm(m+1)+B2*im(m+1));
                a4=a4-(A3*rm(m+1)+B3*im(m+1));

                if(m>=2)
                    a11=a11+m*(m-1)*(A1*rm(m-2+1)+B1*im(m-2+1));
                    a12=a12+m*(m-1)*(B1*rm(m-2+1)-A1*im(m-2+1));
                end
                if(m>=1)
                    a13=a13+m*(A2*rm(m-1+1)+B2*im(m-1+1));
                    a14=a14-m*(A3*rm(m-1+1)+B3*im(m-1+1));
                    a23=a23+m*(B2*rm(m-1+1)-A2*im(m-1+1));
                    a24=a24-m*(B3*rm(m-1+1)-A3*im(m-1+1));
                end
                a33=a33+(A4*rm(m+1)+B4*im(m+1));
                a34=a34-(A5*rm(m+1)+B5*im(m+1));
                a44=a44+(A6*rm(m+1)+B6*im(m+1));
            end

            a1=a1/r;
            a2=a2/r;
            a3=a3/r;
            a4=a4/r;
            a11=a11/r^2;
            a12=a12/r^2;
            a13=a13/r^2;
            a14=a14/r^2;
            a23=a23/r^2;
            a24=a24/r^2;
            a33=a33/r^2;
            a34=a34/r^2;
            a44=a44/r^2;
            a22=-a11;

            crScal=(c/r)/scalFactor;

            dVdx=crScal*(a1+s*a4);
            dVdy=crScal*(a2+t*a4);
            dVdz=crScal*(a3+u*a4);

            if(spherDerivs)
                %Convert the derivatives to spherical coordinates.
                J=calcSpherInvJacob(pointCur)';

                gradV(:,:,curPoint)=J*[dVdx;dVdy;dVdz];
            else%If a gradient in Cartesian coordinates is desired.
                gradV(1,:,curPoint)=dVdx;
                gradV(2,:,curPoint)=dVdy;
                gradV(3,:,curPoint)=dVdz;
            end

            d2Vdxdx=crScal*(a11+2*s*a14+a4/r+s^2*a44-s^2*a4/r);
            d2Vdydy=crScal*(a22+2*t*a24+a4/r+t^2*a44-t^2*a4/r);
            d2Vdzdz=crScal*(a33+2*u*a34+a4/r+u^2*a44-u^2*a4/r);

            d2Vdxdy=crScal*(a12+s*t*a44+s*a24+t*a14-s*t*a4/r);
            d2Vdxdz=crScal*(a13+s*u*a44+s*a34+u*a14-s*u*a4/r);
            d2Vdydz=crScal*(a23+t*u*a44+t*a34+u*a24-t*u*a4/r);

            if(spherDerivs)
                dxdr=J(1,1);
                dxdAz=J(2,1);
                dxdEl=J(3,1);

                dydr=J(1,2);
                dydAz=J(2,2);
                dydEl=J(3,2);

                dzdr=J(1,3);
                dzdAz=J(2,3);
                dzdEl=J(3,3);

                H=calcSpherInvHessian(pointCur);

                d2xdrdr=H(1,1,1);
                d2xdAzdAz=H(2,2,1);
                d2xdEldEl=H(3,3,1);
                d2xdrdAz=H(1,2,1);
                d2xdrdEl=H(1,3,1);
                d2xdAzdEl=H(2,3,1);

                d2ydrdr=H(1,1,2);
                d2ydAzdAz=H(2,2,2);
                d2ydEldEl=H(3,3,2);
                d2ydrdAz=H(1,2,2);
                d2ydrdEl=H(1,3,2);
                d2ydAzdEl=H(2,3,2);

                d2zdrdr=H(1,1,3);
                d2zdAzdAz=H(2,2,3);
                d2zdEldEl=H(3,3,3);
                d2zdrdAz=H(1,2,3);
                d2zdrdEl=H(1,3,3);
                d2zdAzdEl=H(2,3,3);

                %d2Vdrdr
                HessianV(1,1,:,curPoint)=d2Vdydy*dydr^2+2*d2Vdydz*dydr*dzdr+d2Vdzdz*dzdr^2+2*dxdr*dzdr*d2Vdxdz+2*dxdr*dydr*d2Vdxdy+dxdr^2*d2Vdxdx+dVdx*d2xdrdr+dVdy*d2ydrdr+dVdz*d2zdrdr;
                %d2VdAzdAz
                HessianV(2,2,:,curPoint)=d2Vdzdz*dzdAz^2+2*dydAz*dzdAz*d2Vdydz+dydAz^2*d2Vdydy+dVdy*d2ydAzdAz+dVdz*d2zdAzdAz+d2xdAzdAz*dVdx+2*dxdAz*dzdAz*d2Vdxdz+2*dxdAz*dydAz*d2Vdxdy+dxdAz^2*d2Vdxdx;
                %d2VdEldEl
                HessianV(3,3,:,curPoint)=dzdEl^2*d2Vdzdz+dVdz*d2zdEldEl+d2ydEldEl*dVdy+2*dydEl*dzdEl*d2Vdydz+dydEl^2*d2Vdydy+d2xdEldEl*dVdx+2*dxdEl*dzdEl*d2Vdxdz+2*dxdEl*dydEl*d2Vdxdy+dxdEl^2*d2Vdxdx;
                %d2VdrdAz
                HessianV(1,2,:,curPoint)=dzdAz*d2Vdydz*dydr+dydAz*d2Vdydy*dydr+d2Vdzdz*dzdAz*dzdr+dydAz*d2Vdydz*dzdr+dzdAz*dxdr*d2Vdxdz+dxdAz*dzdr*d2Vdxdz+dydAz*dxdr*d2Vdxdy+dxdAz*dydr*d2Vdxdy+dVdx*d2xdrdAz+dVdy*d2ydrdAz+dVdz*d2zdrdAz+dxdAz*dxdr*d2Vdxdx;
                HessianV(2,1,:,curPoint)=HessianV(1,2,:,curPoint);
                %d2VdrdEl
                HessianV(1,3,:,curPoint)=dzdEl*d2Vdydz*dydr+dydEl*d2Vdydy*dydr+dzdEl*d2Vdzdz*dzdr+dydEl*d2Vdydz*dzdr+dzdEl*dxdr*d2Vdxdz+dxdEl*dzdr*d2Vdxdz+dVdx*d2xdrdEl+dVdy*d2ydrdEl+dVdz*d2zdrdEl+dydEl*dxdr*d2Vdxdy+dxdEl*dydr*d2Vdxdy+dxdEl*dxdr*d2Vdxdx;
                HessianV(3,1,:,curPoint)=HessianV(1,3,:,curPoint);
                %d2VdAzdEl
                HessianV(2,3,:,curPoint)=dzdEl*d2Vdzdz*dzdAz+dzdEl*dydAz*d2Vdydz+dydEl*dzdAz*d2Vdydz+dVdy*d2ydAzdEl+dVdz*d2zdAzdEl+dydEl*dydAz*d2Vdydy+d2xdAzdEl*dVdx+dzdEl*dxdAz*d2Vdxdz+dxdEl*dzdAz*d2Vdxdz+dydEl*dxdAz*d2Vdxdy+dxdEl*dydAz*d2Vdxdy+dxdEl*dxdAz*d2Vdxdx;
                HessianV(3,2,:,curPoint)=HessianV(2,3,:,curPoint);
            else
                HessianV(1,1,:,curPoint)=d2Vdxdx;
                HessianV(2,2,:,curPoint)=d2Vdydy;
                HessianV(3,3,:,curPoint)=d2Vdzdz;

                HessianV(1,2,:,curPoint)=d2Vdxdy;
                HessianV(2,1,:,curPoint)=HessianV(1,2,:,curPoint);
                HessianV(1,3,:,curPoint)=d2Vdxdz;
                HessianV(3,1,:,curPoint)=HessianV(1,3,:,curPoint);
                HessianV(2,3,:,curPoint)=d2Vdydz;
                HessianV(3,2,:,curPoint)=HessianV(2,3,:,curPoint);
            end
        end
    end
end

if(spherDerivs&&systemType==2)
    %Flip signs of the elevation terms reflecting the
    %difference in definition between systems 0 and 2.
    gradV(3,:)=-gradV(3,:);
    HessianV(1,3,:)=-HessianV(1,3,:);
    HessianV(3,1,:)=-HessianV(3,1,:);
    HessianV(2,3,:)=-HessianV(2,3,:);
    HessianV(3,2,:)=-HessianV(3,2,:);
end
end

function [SinVec,CosVec]=calcSinCosTerms(lambda,M)
    %Compute sin(m*lambda) and cos(m*lambda) for m=0 to m=M.
    SinVec=zeros(M+1,1);
    CosVec=zeros(M+1,1);
    %Explicitly set the first two terms.
    SinVec(0+1)=0;
    CosVec(0+1)=1;
    SinVec(1+1)=sin(lambda);
    CosVec(1+1)=cos(lambda);
    %Use a double angle identity to get the second order term.
    SinVec(2+1)=2*SinVec(1+1)*CosVec(1+1);
    CosVec(2+1)=1-2*SinVec(1+1)^2;
    %Use a two-part recursion for the rest of the terms.
    for m=3:M
        SinVec(m+1)=2*CosVec(1+1)*SinVec(m-1+1)-SinVec(m-2+1);
        CosVec(m+1)=2*CosVec(1+1)*CosVec(m-1+1)-CosVec(m-2+1);
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [mu, varsigma] = vargplvmPosteriorMeanVarPar(model, X, varX);

% VARGPLVMPOSTERIORMEANVARPAR Mean and variances of the posterior at points given by X, parallel version.
% FORMAT
% DESC returns the posterior mean and variance for a given set of
% points.
% ARG model : the model for which the posterior will be computed.
% ARG X : variational mean in the latent space for which posterior is computed.
% ARG varX : variational variances in the latent space for which posterior is computed (assumed zero if not present).
% RETURN mu : the mean of the posterior distribution.
% RETURN sigma : the variances of the posterior distributions.
%
% SEEALSO : vargplvmPosteriorMeanVar, gpPosteriorMeanVar, vargplvmCreate
%
% COPYRIGHT : Michalis K. Titsias and Neil D. Lawrence, 2009, 2011

% VARGPLVM

% do prediction by replacing the variational distribution with a delta function
%model.K_uf = kernCompute(model.kern, model.X_u, model.vardist.means);
%model.A = (1/model.beta)*model.K_uu + model.K_uf*model.K_uf';
%[model.Ainv, U] = pdinv(model.A);
%[mu1, varsigma1] = gpPosteriorMeanVar(model, vardistX.means);

% Find exactly the mean and the variances of the predictive distribution
% (which is not Gaussian, however its moments can be computed in closed-form)

if nargin < 3
  vardistX.covars = repmat(0.0, size(X, 1), size(X, 2));%zeros(size(X, 1), size(X, 2));
else
  vardistX.covars = varX;
end
vardistX.latentDimension = size(X, 2);
vardistX.numData = size(X, 1);
%model.vardist.covars = 0*model.vardist.covars;
vardistX.means = X;
%model = vargplvmUpdateStats(model, model.X_u);

Ainv = model.P1' * model.P1; % size: NxN

if ~isfield(model,'alpha')
    model.alpha = Ainv*model.Psi1'*model.m; % size: 1xD
end
Psi1_star = kernVardistPsi1Compute(model.kern, vardistX, model.X_u);

% mean prediction
mu = Psi1_star*model.alpha; % size: 1xD

if nargout > 1
   %
   % precomputations
   vard = vardistCreate(zeros(1,model.q), model.q, 'gaussian');
   Kinvk = (model.invK_uu - (1/model.beta)*Ainv);
   %
   for i=1:size(vardistX.means,1)
      %
      vard.means = vardistX.means(i,:);
      vard.covars = vardistX.covars(i,:);
      % compute psi0 term
      Psi0_star = kernVardistPsi0Compute(model.kern, vard);
      % compute psi2 term
      Psi2_star = kernVardistPsi2Compute(model.kern, vard, model.X_u);

      vars = Psi0_star - sum(sum(Kinvk.*Psi2_star));

     % parfor j=1:model.d
     %    %[model.alpha(:,j)'*(Psi2_star*model.alpha(:,j)), mu(i,j)^2]
     %    varsigma(i,j) = model.alpha(:,j)'*(Psi2_star*model.alpha(:,j)) - mu(i,j)^2;
     % end
     % varsigma(i,:) = varsigma(i,:) + vars;
      %
   end
   %
   %if isfield(model, 'beta')
   %   varsigma = varsigma + (1/model.beta);
   %end
   %
end

% Rescale the mean
mu = mu.*repmat(model.scale, size(vardistX.means,1), 1);

% Add the bias back in
mu = mu + repmat(model.bias, size(vardistX.means,1), 1);

% rescale the variances
if nargout > 1
    varsigma = varsigma.*repmat(model.scale.*model.scale, size(vardistX.means,1), 1);
end

author: SheffieldML

repo: vargplvm


text:
function [C,dt] = VBA_getNtuples(k,n,verbose)
% gets all n-draws (with replacement) from a k-urn.
% [C] = VBA_getNtuples(k,n,verbose)
% This function construct all the n-tuples that can be obtained (with
% replacement) from an urn containing k marbles.
% A couple of examples is worth a long description:
%
% VBA_getNtuples(2,3)
%
% ans =
%
%      1     1     1     1     2     2     2     2
%      1     1     2     2     1     1     2     2
%      1     2     1     2     1     2     1     2
%
% VBA_getNtuples(3,2)
%
% ans =
%
%      1     1     1     2     2     2     3     3     3
%      1     2     3     1     2     3     1     2     3
%
% IN:
%   - k: the cardinality of the set, whose permutations are replicated over
%   n dimensions
%   - n: the dimension of the tuple, whose entries are elements of the set
%   {1,2,...,k}.
% OUT:
%   - C: a nXn^p matrix, whose columns are the n-draws
% Thanks to Antonin "magic" Ancelle!

try; verbose;catch;verbose=0;end
tic;
np = k^n;
if verbose
    fprintf(1,['Getting the ',num2str(np),' ',num2str(n),'-draws from ',num2str(k),'-urn...'])
end
C = zeros(n,np); % pre-allocate k-permutations
for ind = 1:n
    for j = 1:k^(ind-1)
        nbch = np/(k^(ind-1));
        C(ind,(1+(j-1)*nbch):j*nbch) = ceil((1:nbch)/(nbch/k));
    end
end
dt = toc;
if verbose
    fprintf(1,[' OK. (took ',num2str(dt),' seconds)'] )
    fprintf(1,'\n')
end

return

author: MBB-team

repo: VBA-toolbox


text:
function geometry_test024 ( )

%*****************************************************************************80
%
%% TEST024 tests DUAL_SHAPE_3D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST024\n' );
  fprintf ( 1, '  DUAL_SHAPE_3D finds the dual of a polyhedron.\n' );
%
%  Get the dodecahedron shape.
%
  [ point_num1, edge_num1, face_num1, face_order_max1 ] = dodec_size_3d ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num1 );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num1 );
  fprintf ( 1, '  Number of faces =    %d\n', face_num1 );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max1 );

  [ point_coord1, face_order1, face_point1 ] = dodec_shape_3d ( point_num1, ...
    face_num1, face_order_max1 );
%
%  Get the dual.
%
  [ point_num2, edge_num2, face_num2, face_order_max2 ] = dual_size_3d ( ...
    point_num1, edge_num1, face_num1, face_order_max1, point_coord1, ...
    face_order1, face_point1 );

  [ point_coord2, face_order2, face_point2 ] = dual_shape_3d ( point_num1, ...
    face_num1, face_order_max1, point_coord1, face_order1, face_point1, ...
    point_num2, face_num2, face_order_max2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num2 );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num2 );
  fprintf ( 1, '  Number of faces =    %d\n', face_num2 );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max2 );

  shape_print_3d ( point_num2, face_num2, face_order_max2, ...
    point_coord2, face_order2, face_point2 );

  return
end


text:
function A = meanGPexact(mean,cov,x,y, hypz,z,i)

% Mean function being the predictive mean of a GP model:
%
% m(z) = posterior mean of GP at location z as given by
% m(z) = gp(hyp,@infExact,mean,cov,@likGauss,x,y, z) where
% hyp.mean = hyp_mean; hyp.lik = log(sn); hyp.cov = hyp.cov;
%
% The hyperparameters are:
%
% hypz = [ hyp_cov
%          log(sn)
%          hyp_mean ]
%
% where hyp_cov are the covariance function hyperparameters, sn is the
% noise variance of the Gaussian likelihood and hyp_mean are the mean
% function hyperparameters.
%
% Copyright (c) by Hannes Nickisch, 2014-11-01.
%
% See also MEANFUNCTIONS.M and MEANGP.M.

if nargin<4, error('GP must be specified.'), end           % check for dimension
if isempty(mean), mean = @meanZero; end              % set default and make cell
if ~iscell(mean), mean = {mean};    end
if isempty(cov),  cov  = @covSEiso; end
if ~iscell(cov),  cov  = {cov};     end
nms = feval(mean{:}); ncs = feval(cov{:});     % number of hyperparameter string
if nargin<6, A = [ncs,'+1+',nms]; return, end % report number of hyperparameters

[nz,D] = size(z); n = size(x,1);
nc = eval(ncs); nm = eval(nms);
hyp = rewrap(struct('cov',zeros(nc,1),'lik',0,'mean',zeros(nm,1)),hypz);

m  = feval(mean{:},hyp.mean,x);
mz = feval(mean{:},hyp.mean,z);
K  = feval(cov{:}, hyp.cov, x);
kz = feval(cov{:}, hyp.cov, x,z);

sn2 = exp(2*hyp.lik);                               % noise variance of likGauss
if sn2<1e-6                        % very tiny sn2 can lead to numerical trouble
  L = chol(K+sn2*eye(n)); sl =   1;   % Cholesky factor of covariance with noise
else
  L = chol(K/sn2+eye(n)); sl = sn2;                       % Cholesky factor of B
end
iKs = @(t) solve_chol(L,t)/sl;                       % iKs(t) = (K+sn2*eye(n))\t
alpha = iKs(y-m);
if nargin==6                                               % eval posterior mean
  A = mz+kz'*alpha;
else
  if i<=nc                                          % covariance function hypers
    dK = feval(cov{:},hyp.cov,x,[],i);
    dkz = feval(cov{:},hyp.cov,x,z,i);
    A = dkz'*alpha-kz'*iKs(dK*alpha);
  elseif i==nc+1                              % likelihood function parameter sn
    A = -2*sn2*kz'*iKs(alpha);
  else                                                 % mean function parameter
    dm  = feval(mean{:},hyp.mean,x,i-nc-1);
    dmz = feval(mean{:},hyp.mean,z,i-nc-1);
    A = dmz-kz'*iKs(dm);
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Replication of the trivariate VAR in Stock and Watson (2001, JEP).
% Figure 1 and Table 1.B.
%==========================================================================
% The VAR Toolbox 3.0 is required to run this code. To get the
% latest version of the toolboxes visit:
% https://github.com/ambropo/VAR-Toolbox
%==========================================================================
% Ambrogio Cesa Bianchi, November 2020
% ambrogio.cesabianchi@gmail.com

%% PRELIMINARIES
%==========================================================================
clear all; clear session; close all; clc
warning off all

% Load data
[xlsdata, xlstext] = xlsread('SW2001_Data.xlsx','Sheet1');
X = xlsdata;
dates = xlstext(3:end,1);
vnames_long = xlstext(1,2:end);
vnames = xlstext(2,2:end);
nvar = length(vnames);
data   = Num2NaN(xlsdata);
% Store variables in the structure DATA
for ii=1:length(vnames)
    DATA.(vnames{ii}) = data(:,ii);
end
% Convert the first date to numeric
year = str2double(xlstext{3,1}(1:4));
quarter = str2double(xlstext{3,1}(6));
% Observations
nobs = size(data,1);

%% VAR ESTIMATION
%==========================================================================
% Set deterministics for the VAR
det = 1;
% Set number of nlags
nlags = 4;
% Estimate VAR
[VAR, VARopt] = VARmodel(X,nlags,det);
% Print estimation on screen
VARopt.vnames = vnames;
[TABLE, beta] = VARprint(VAR,VARopt,2);

%% COMPUTE IR AND VD
%==========================================================================
% Set options some options for IRF calculation
VARopt.nsteps = 24;
VARopt.ident = 'short';
VARopt.vnames = vnames_long;
VARopt.FigSize = [26,12];
% Compute IRF
[IRF, VAR] = VARir(VAR,VARopt);
% Compute error bands
[IRinf,IRsup,IRmed,IRbar] = VARirband(VAR,VARopt);
% Plot
VARirplot(IRbar,VARopt,IRinf,IRsup);

% Compute VD
[VD, VAR] = VARvd(VAR,VARopt);
% Compute VD error bands
[VDinf,VDsup,VDmed,VDbar] = VARvdband(VAR,VARopt);
% Plot VD
VARvdplot(VDbar,VARopt);

%% Print Table 1.B on screen
%==========================================================================
% Retrieve Forecast Error Variance Decomposition
FEVD_Table(1, :) = VD(1,:,1);
FEVD_Table(2, :) = VD(4,:,1);
FEVD_Table(3, :) = VD(8,:,1);
FEVD_Table(4, :) = VD(12,:,1);

FEVD_Table(5, :) = VD(1,:,2);
FEVD_Table(6, :) = VD(4,:,2);
FEVD_Table(7, :) = VD(8,:,2);
FEVD_Table(8, :) = VD(12,:,2);

FEVD_Table(9, :) = VD(1,:,3);
FEVD_Table(10,:) = VD(4,:,3);
FEVD_Table(11,:) = VD(8,:,3);
FEVD_Table(12,:) = VD(12,:,3);

% Print on screen
disp(' ')
disp('Variance Decomposition of Inflation (t=1,4,8,12)')
disp('---------------------------------------------------')
mprint(FEVD_Table(1:4,:));
disp('Variance Decomposition of Unemployment (t=1,4,8,12)')
disp('---------------------------------------------------')
mprint(FEVD_Table(5:8,:));
disp('Variance Decomposition of Fed Funds (t=1,4,8,12)')
disp('---------------------------------------------------')
mprint(FEVD_Table(9:12,:));

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function R = quatToMat(b, c, d, qx, qy, qz, qfac, pixdim)
%
% R = quatToMat(qb, qc, qd, qx, qy, qz, qfac, pixdim)
%
% Convert a quaternion to a 4x4 rotation/scaling/translation matrix. Based on the
% nifti1 spec implemented in nifti1_io.c.
%
% 2008.10.01 RFD: wrote it.
%

if(isstruct(b))
    c = b.quatern_c;
    d = b.quatern_d;
    qx = b.quatern_x;
    qy = b.quatern_y;
    qz = b.quatern_z;
    qfac = b.qfac;
    pixdim = [b.dx b.dy b.dz];
    b = b.quatern_b;
end

R = zeros(4);

R(4,:) = [0 0 0 1];

% compute a from b,c,d
a = 1.0 - (b*b + c*c + d*d);
if( a < 1.e-7 )                     % special case
     a = 1.0 / sqrt(b*b+c*c+d*d);
     b = b*a; c = c*a; d = d*a;     % normalize (b,c,d) vector
     a = 0.0;                       % a = 0 ==> 180 degree rotation
else
     a = sqrt(a);                   % angle = 2*arccos(a)
end

% build rotation matrix, including scaling factors for voxel sizes
if length(pixdim) == 3
    xd = pixdim(1); yd = pixdim(2); zd = pixdim(3);
else
    xd = pixdim(1); yd = pixdim(2); zd = 1;
end
if(xd<=0.0), xd = 1.0; end
if(yd<=0.0), yd = 1.0; end
if(zd<=0.0), zd = 1.0; end

if(qfac < 0.0), zd = -zd; end       % left handedness?

R(1,1) =       (a*a+b*b-c*c-d*d) * xd ;
R(1,2) = 2.0 * (b*c-a*d        ) * yd ;
R(1,3) = 2.0 * (b*d+a*c        ) * zd ;
R(2,1) = 2.0 * (b*c+a*d        ) * xd ;
R(2,2) =       (a*a+c*c-b*b-d*d) * yd ;
R(2,3) = 2.0 * (c*d-a*b        ) * zd ;
R(3,1) = 2.0 * (b*d-a*c        ) * xd ;
R(3,2) = 2.0 * (c*d+a*b        ) * yd ;
R(3,3) =       (a*a+d*d-c*c-b*b) * zd ;

% load offsets

R(1:3,4) = [qx;qy;qz];

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
obSize gives the size.
a = reshape(obSlice,42,51);
b = flipup(a);

function [sagSlice,obSlice,sagX,sagY] = ...
	mrUpSagShowObl(volume,sagSize,numSlices,curSag,obPts,obSize,obX,obY,volselpts)
%
%NAME: [sagSlice,obSlice,sagX,sagY] = ...
%	mrUpSagShowObl(volume,sagSize,numSlices,curSag,obPts,obSize,obX,obY,volselpts)
%AUTHOR:  Poirson
%DATE:    08.07.96
%HISTORY  Started with mrUpSagShowObl.m
%BUGS:

% relevant code to display the interpolated oblique plane
global obwin
if ~isempty(obPts)
	figure(obwin);
	[obSlice,sagX,sagY] = ...
		mrDispOblVol(volume,sagSize,numSlices,obPts,obSize,curSag,volselpts);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [obSlice, sagX, sagY] = mrDispOblVol(volume,sagSize,numSlices,obPts,obSize,curSag,volselpts)
%
global obwin volslimin2 volslimax2;

if ~isempty(obPts)
	sagX = [curSag,curSag];
	sagY = [0,obSize(1)];
	tmp = volume;
	if ~isempty(volselpts)
		tmp(volselpts) = -1*ones(1,length(volselpts));
	end
	obSlice = mrExtractImgVol(tmp, sagSize, numSlices, obPts);
	figure(obwin);
	myShowImageVol(obSlice,obSize,max(obSlice)*get(volslimin2,'value'),max(obSlice)*get(volslimax2,'value'),sagX,sagY);
end


text:
function [ output_meta ] = meta2sicd_s1manifest( domnode )
%META2SICD_S1MANIFEST Converts Sentinel-1 manifest.safe description into a SICD-style metadata structure
%
% In general, the Sentinel-1 manifest.safe file doesn't have much real
% metadata directly in it.  However, it is the document that points to all
% of the data and metadata files in the package.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Setup
% Overall wrapping node is in a namespace that makes things messier
% (than they already are) in a XPath call, so we go down one level.
domnode = domnode.getFirstChild;
xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();

%% CollectionInfo
platform_xp_str = ['metadataSection/' ...
    'metadataObject[@ID="platform"]/' ...
    'metadataWrap/' ...
    'xmlData/' ...
    '*[local-name()=''platform'']/'];
output_meta.CollectionInfo.CollectorName = [
    char(xp.evaluate(...
    [platform_xp_str '*[local-name()=''familyName'']'],domnode)) ...
    char(xp.evaluate(...
    [platform_xp_str '*[local-name()=''number'']'],domnode))];
% The following core name is unique within all Sentinel-1 collections.
% Perhaps we would rather use the NGA CoreName standard.
% We also may need to further add to this CoreName, if each collect has
% multiple swaths and bursts which must each go into their own SICD, each
% requiring their own unique CoreNames.
standalone_xp_str = ['metadataSection/' ...
    'metadataObject[@ID="generalProductInformation"]/' ...
    'metadataWrap/' ...
    'xmlData/' ...
    '*[local-name()=''standAloneProductInformation'']/'];
% Since each TOPSAR product has multiple CDPs (bursts), we really don't
% have enough information here to make a unique identifier for each CDP.
% We will defer generating the CoreName until meta2sicd_s1product.m
% output_meta.CollectionInfo.CoreName = char(xp.evaluate(...
%     [standalone_xp_str '*[local-name()=''missionDataTakeID'']'],domnode));
output_meta.CollectionInfo.CollectType='MONOSTATIC';
output_meta.CollectionInfo.RadarMode.ModeID = ...
    char(xp.evaluate(...
    [platform_xp_str ...
    '*[local-name()=''instrument'']/' ...
    '*[local-name()=''extension'']/' ...
    '*[local-name()=''instrumentMode'']/' ...
    '*[local-name()=''mode'']'],domnode));
if strcmp(output_meta.CollectionInfo.RadarMode.ModeID, 'SM')
    output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
else
    % Actually TOPSAR.  Not what we normally think of for DYNAMIC STRIPMAP,
    % but it is definitely not SPOTLIGHT, and doesn't seem to be regular
    % STRIPMAP either.
    output_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
end
output_meta.CollectionInfo.Classification='UNCLASSIFIED';

%% ImageCreation
% These fields describe the final SLC post processing step, not the SLC
% processing itself, although in datasets we have seen, this information
% has been the same.
processing_xp_str = ['metadataSection/' ...
    'metadataObject[@ID="processing"]/' ...
    'metadataWrap/' ...
    'xmlData/' ...
    '*[local-name()=''processing'']/'];
facility_xp_str = [processing_xp_str '*[local-name()=''facility'']/'];
software_xp_str = [facility_xp_str '*[local-name()=''software'']/'];
output_meta.ImageCreation.Application = [...
    char(xp.evaluate([software_xp_str '@name'],domnode)) ' '...
    char(xp.evaluate([software_xp_str '@version'],domnode))];
output_meta.ImageCreation.DateTime = datenum( ...
    char(xp.evaluate([processing_xp_str '@stop'],domnode)), ...
    'yyyy-mm-ddTHH:MM:SS.FFF');
output_meta.ImageCreation.Site = [ ...
    char(xp.evaluate([facility_xp_str '@name'],domnode)) ', ' ...
    char(xp.evaluate([facility_xp_str '@site'],domnode)) ', ' ...
    char(xp.evaluate([facility_xp_str '@country'],domnode))];
output_meta.ImageCreation.Profile='Prototype';

%% Timeline
% The start/stop times given in the manifest file appear to be the zero
% Doppler times for the first and the last lines, not the total collect
% time, which seems to be given by firstLineSensingTime/lastLineSensingTime
% in the annotation XML file. These start/stop times in the manifest file
% also seem to be replicated in the annotation XML file.
% acquistion_xp_str = ['metadataSection/' ...
%     'metadataObject[@ID="acquisitionPeriod"]/' ...
%     'metadataWrap/' ...
%     'xmlData/' ...
%     '*[local-name()=''acquisitionPeriod'']/'];
% char(xp.evaluate(...
%     [acquistion_xp_str '*[local-name()=''startTime'']'],domnode));
% char(xp.evaluate(...
%     [acquistion_xp_str '*[local-name()=''stopTime'']'],domnode));

%% RadarCollection
num_pols = str2double(xp.evaluate(...
    ['count(' standalone_xp_str '*[local-name()=''transmitterReceiverPolarisation''])'],domnode));
for i = 1:num_pols
    pol = char(xp.evaluate(...
        [standalone_xp_str '*[local-name()=''transmitterReceiverPolarisation''][' num2str(i) ']'],domnode));
    output_meta.RadarCollection.RcvChannels.ChanParameters(i).TxRcvPolarization=[pol(1) ':' pol(2)];
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function spp_cache_convX_features(imdb, varargin)
% spp_cache_convX_features(imdb, varargin)
%   Computes last conv features and saves them to disk.
%
%   Keys that can be passed in:
%
%   imdb_copy_from    An imdb with partial scales cache which can copy from
%   start             Index of the first image in imdb to process
%   end               Index of the last image in imdb to process
%   spm_im_size       Scales of image resize to
%   spp_params_def    spp pooling def file, only for adaptive scale
%   roidb             only for adaptive scale
%   net_def_file      Path to the Caffe CNN to use
%   net_file          Path to the Caffe CNN to use
%   cache             Path to the precomputed feature cache
%   cache_copy_from   Path to the precomputed feature cache to copy from
%
% Adapted from spp code written by Ross Girshick
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

ip = inputParser;
ip.addRequired('imdb', @isstruct);
ip.addOptional('imdb_copy_from', struct(), @isstruct);
ip.addOptional('start', 1, @isscalar);
ip.addOptional('end', 0, @isscalar);
ip.addOptional('spm_im_size', 688, @ismatrix);
ip.addOptional('spp_params_def', '', @isstr);
ip.addOptional('roidb', '', @isstruct);
ip.addOptional('net_def_file', '', @isstr);
ip.addOptional('net_file', '', @isstr);
ip.addOptional('cache', 'feat_cache', @isstr);
ip.addOptional('cache_copy_from', 'feat_cache', @isstr);

ip.parse(imdb, varargin{:});
opts = ip.Results;

% load caffe model
spp_load_model(spp_create_model(opts.net_def_file, opts.net_file, opts.spp_params_def));

% Where to save feature cache
opts.output_dir = ['./feat_cache/' opts.cache '/' imdb.name '/'];
copy_from_previous_cache = false;
if ~isempty(opts.imdb_copy_from) && isfield(opts.imdb_copy_from, 'name');
    copy_from_previous_cache = true;
    opts.output_dir_copy_from = ['./feat_cache/' opts.cache_copy_from '/' opts.imdb_copy_from.name '/'];
end
mkdir_if_missing(opts.output_dir);
conf = spp_config();

% Log feature extraction
timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
diary_file = [opts.output_dir 'spp_train_' opts.cache '_'  timestamp '.txt'];
diary(diary_file);
fprintf('Logging output in %s\n', diary_file);

fprintf('\n\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n');
fprintf('Feature caching options:\n');
disp(opts);
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\n');

image_ids = imdb.image_ids;
if opts.end == 0
  opts.end = length(image_ids);
end

total_time = 0;
count = 0;
for i = opts.start:opts.end
  fprintf('%s: cache features: %d/%d\n', procid(), i, opts.end);

  save_file = [opts.output_dir image_ids{i} '.mat'];
  if exist(save_file, 'file') ~= 0
    fprintf(' [already exists]\n');

    broken = false;
    if 1 % check all cache files
        try
            d = load(save_file);
            assert(isfield(d, 'feat'));
            assert(~isempty(d.feat));
            if ~isempty(opts.spm_im_size)
                assert(isequal(d.feat.scale, opts.spm_im_size));
                assert(length(d.feat.rsp) == length(opts.spm_im_size));
            end
        catch
            broken = true;
        end
    end

    if ~broken
        continue;
    end
  end
  mkdir_if_missing(fileparts(save_file));

  feat_cache = [];
  if copy_from_previous_cache
      load_file = [opts.output_dir_copy_from image_ids{i} '.mat'];
      if exist(load_file, 'file') == 0
         fprintf(' [missing]\n');
      else
         cache = load(load_file);
         feat_cache = cache.feat;
      end
  end

  tot_th = tic;
  count = count + 1;

  d.gt = []; d.overlap = []; d.boxes = []; d.class = []; d.feat = [];
  try
        im = imread(imdb.image_at(i));
        th = tic;

        d.feat = spp_features_convX(im, opts.spm_im_size, feat_cache, conf.use_gpu);

        fprintf(' [features: %.3fs]\n', toc(th));
  catch err
        fprintf(' extract feature error : %s \n', err.message);
  end

  th = tic;
  save(save_file, '-struct', 'd');
  fprintf(' [saving:   %.3fs]\n', toc(th));

  total_time = total_time + toc(tot_th);
  fprintf(' [avg time: %.3fs (total: %.3fs)]\n', ...
      total_time/count, total_time);
end

author: ShaoqingRen


text:
function pass = test_isNotMultOrDiff(~)
%TEST_ISNOTMULTORDIFF    Test that the property of operator not being
%diff/integration is dealt with correctly.

% Construct the primitive operatorBlocks
dom = [0, 2];
[Z, I, D, C, M] = linop.primitiveOperators(dom);
M = M(chebfun(@sin, dom));

% Check that they have the expected value of the property
pass(1) = ( Z.isNotDiffOrInt == 1 );
pass(2) = ( I.isNotDiffOrInt == 1 );
pass(3) = ( D.isNotDiffOrInt == 0 );
pass(4) = ( C.isNotDiffOrInt == 0 );
pass(5) = ( M.isNotDiffOrInt == 1 );

% Operations on operatorBlocks
A = Z + I;
pass(6) = ( A.isNotDiffOrInt == 1 );
A = Z + D;
pass(7) = ( A.isNotDiffOrInt == 0 );
A = M*Z;
pass(8) = ( A.isNotDiffOrInt == 1 );
A = C*Z;
pass(9) = ( A.isNotDiffOrInt == 1 );
A = C + M;
pass(10) = ( A.isNotDiffOrInt == 0 );
A = 2*M;
pass(11) = ( A.isNotDiffOrInt == 1 );
A = 2*D;
pass(12) = ( A.isNotDiffOrInt == 0 );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [scaled] = NewScale(Data, Lower, Upper)
% [scaled] = Scale(Data);
% [scaled] = Scale(Data, Lower, Upper);
%
% scale the elements of all the column vectors in
% the range of [Lower Upper]. (default [-1 1])

if (nargin<3)
   Lower = -1;
   Upper = 1;
elseif (Lower > Upper)
   disp (['Wrong Lower or Upper values!']);
end

[MaxV, I]=max(Data);
[MinV, I]=min(Data);

[R,C]= size(Data);

scaled=(Data-ones(R,1)*MinV).*(ones(R,1)*((Upper-Lower)*ones(1,C)./(MaxV-MinV)))+Lower;

scaled(isnan(scaled)) = 1/sqrt(size(scaled, 1));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function xyr = Xr(s)

x = 1 ;

y = s ;

xyr = [x ; y] ;