text:
function [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, fullMipRes] = ftINIT(prepData, tissue, celltype, hpaData, transcrData, metabolomicsData, INITSteps, removeGenes, useScoresForTasks, paramsFT, verbose)
% ftINIT
%   Main function for generates a model using the ftINIT algorithm, based
%   on proteomics and/or transcriptomics and/or metabolomics and/or metabolic
%   tasks. The algorithm is not designed for running with metabolomics only.
%   The function prepINITModel needs to be run first for the template
%   model (such as the generic Human-GEM), but only need to be run once. This
%   function precalculates things independent of the omics to speed up the model
%   generation process, and outputs the prepData, which is input to this function.
%
%   prepData            The prepdata for the model.
%   tissue              tissue to score for. Should exist in either
%                       hpaData.tissues or transcrData.tissues
%   celltype            cell type to score for. Should exist in either
%                       hpaData.celltypes or transcrData.celltypes for this
%                       tissue (opt, default is to use the max values
%                       among all the cell types for the tissue.
%   hpaData             HPA data structure from parseHPA (opt if transcrData
%                       is supplied, default [])
%   transcrData         gene expression data structure (opt if hpaData is
%                       supplied, default []). Used to be called arrayData.
%       genes           cell array with the unique gene names
%       tissues         cell array with the tissue names. The list may not
%                       be unique, as there can be multiple cell types per
%                       tissue.
%       celltypes       cell array with the cell type names for each tissue
%       levels          GENESxTISSUES array with the expression level for
%                       each gene in each tissue/celltype. NaN should be
%                       used when no measurement was performed
%       threshold       a single value or a vector of gene expression
%                       thresholds, above which genes are considered to be
%                       "expressed". default = 1(opt, by default, the mean expression
%                       levels of each gene across all tissues in transcrData
%                       will be used as the threshold values)
%       singleCells     binary value selecting whether to use the
%                       single-cell algorithm to identify expressed genes.
%                       If used, specify cell subpopulations in CELLTYPES
%                       (opt, default [])
%       plotResults     true if single cell probability distributions
%                       should be plotted (opt, default = false)
%   metabolomicsData    cell array with metabolite names that the model
%                       should produce (opt, default [])
%   INITSteps           Specifies the steps in the algorithm. For more info,
%                       see INITStepDesc and getINITSteps.
%                       (opt, default getINITSteps(), which is the standard ftINIT).
%   removeGenes         if true, low-abundance genes will be removed from
%                       grRules, unless they are the only gene associated
%                       with a reaction, or a subunit of an enzyme complex
%                       (see "removeLowScoreGenes" function for details).
%                       If false, grRules will not be modified; however,
%                       genes that were associated only with removed
%                       reactions will not be present in the final model.
%                       (opt, default true).
%   useScoresForTasks   true if the calculated reaction scored should be
%                       used as weights when fitting to tasks (opt, default
%                       true)
%   paramsFT            parameter structure as used by getMILPParams. This
%                       is for the fitTasks step. For the INIT algorithm,
%                       see params (opt, default [])
%   verbose             if true, the MILP progression will be shown.
%                       (opt, default true)
%
%   model                   the resulting model structure
%   metProduction           array that indicates which of the
%                           metabolites in metabolomicsData that could be
%                           produced. Note that this is before the
%                           gap-filling process to enable defined tasks. To
%                           see which metabolites that can be produced in
%                           the final model, use canProduce.
%                           -2: metabolite name not found in model
%                           -1: metabolite found, but it could not be produced
%                           1: metabolite could be produced
%   addedRxnsForTasks       cell array of the reactions which were added in
%                           order to perform the tasks
%   deletedRxnsInINIT       cell array of reactions deleted because they
%                           could not carry flux (INIT requires a
%                           functional input model)
%   fullMipRes              The solver results from the last MILP step run
%
%   This is the main function for automatic reconstruction of models based
%   on the ftINIT algorithm ().
%
%   NOTE: Exchange metabolites should normally not be removed from the model
%   when using this approach, since checkTasks/fitTasks rely on putting specific
%   constraints for each task. The INIT algorithm will remove exchange metabolites
%   if any are present. Use importModel(file,false) to import a model with
%   exchange metabolites remaining.
%
%   Usage: [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, ...
%               fullMipRes] = ...
%               ftINIT(prepData, tissue, celltype, hpaData, transcrData, ...
%               metabolomicsData, INITSteps, removeGenes, useScoresForTasks, ...
%               paramsFT);
%

if nargin < 5
    transcrData = [];
end
if nargin < 6
    metabolomicsData = [];
end
if nargin < 7 || isempty(INITSteps)
    INITSteps = getINITSteps([],'1+1');
end
if nargin < 8 || isempty(removeGenes)
    removeGenes = true;
end
if nargin < 9 || isempty(useScoresForTasks)
    useScoresForTasks = true;
end
if nargin < 10
    paramsFT = [];
end

if nargin < 11
    verbose = true;
end
%Handle detected mets:
%Previously, this was handled by giving a bonus for secreting those metabolites,
%but that doesn't work since the metabolite secretion and uptake can be lost when
%we merge linearly dependent reactions.
%Instead, we need to figure out which reactions either produce or take up the mets.
%We then give a bonus if any of them carry flux.
%To simplify things, we focus on reactions that produce the metabolite (since there must be one such reaction).
%It is still a bit complicated though. In this step, we focus on identifying
%producer reactions. We further reason that the direction doesn't matter -
%we can force one of these reactions in any direction - if it becomes a consumer, it will
%automatically force another producer on as well (otherwise we'll have a net consumption).

if (~isempty(metabolomicsData))
    if length(unique(upper(metabolomicsData))) ~= length(metabolomicsData)
        dispEM('Metabolomics contains the same metabolite multiple times');
    end
    metData = false(numel(metabolomicsData), length(prepData.minModel.rxns)); %one row per metabolite that is a boolean vector
    for i=1:numel(metabolomicsData)
        %Get the matching mets
        metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i}));
        prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers
                     (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers
        %convert the production rxns from refModel to minModel
        prepData.groupIds;
        [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns);
        grpIdsMerged = nan(length(prepData.minModel.rxns),1);
        grpIdsMerged(ia) = prepData.groupIds(ib);

        groupIdsPos = unique(prepData.groupIds(prodRxnsSel));%gets all group ids which includes a production rxn
        groupIdsPos = groupIdsPos(groupIdsPos ~= 0);%remove the 0 id, it means there is no group
        %the other option is that there is a direct match between the rxn id in minModel and refModel:
        posRxns = prepData.refModel.rxns(prodRxnsSel);
        directMatch = ismember(prepData.minModel.rxns, posRxns).';

        metData(i,:) = ismember(grpIdsMerged, groupIdsPos).' | directMatch;
    end
    metData = sparse(metData);
else
    metData = [];
end

% Get rxn scores and adapt them to the minimized model
origRxnScores = scoreComplexModel(prepData.refModel,hpaData,transcrData,tissue,celltype);
origRxnScores(origRxnScores > -0.1 & origRxnScores <= 0) = -0.1;%we don't want reaction scores that are exactly 0 (or close), this causes problems in the milp
origRxnScores(origRxnScores < 0.1 & origRxnScores > 0) = 0.1;

rxnsTurnedOn = false(length(prepData.minModel.rxns),1);
fluxes = zeros(length(prepData.minModel.rxns),1);

rxnsToIgnoreLastStep = [1;1;1;1;1;1;1;1];

%We assume that all essential rxns are irrev - this is taken care of in
%prepINITModel. We then use an initial flux "from last run" of 0.1 for all
%reactions. This is used for knowing what flux should be forced through an
%essential rxn.
fluxes = ones(length(prepData.minModel.rxns), 1).*0.1;

for initStep = 1:length(INITSteps)
    disp(['ftINIT: Running step ' num2str(initStep)])
    stp = INITSteps{initStep};

    if any ((rxnsToIgnoreLastStep - stp.RxnsToIgnoreMask) < 0)
        dispEM('RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.');
    end
    rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask;

    mm = prepData.minModel;

    if (~isempty(stp.MetsToIgnore))
        if (~isempty(stp.MetsToIgnore.simpleMets))
            %Here, we remove simple metabolites that will not really affect the milp but
            %are very common in the S matrix. For example H2O, H+, etc.
            %It is also possible to leave compartments untouched, for example the i compartment in the mitochondria (for H+).
            metsToRem = ismember(mm.metNames,stp.MetsToIgnore.simpleMets.mets);
            compsToKeep = find(ismember(mm.comps, stp.MetsToIgnore.simpleMets.compsToKeep));
            metsToRem = metsToRem & ~ismember(mm.metComps, compsToKeep);
            mm.S(metsToRem,:) = 0;
        end
    end

    %Set up the reaction scores and essential rxns
    rxnsToIgnore = getRxnsFromPattern(stp.RxnsToIgnoreMask, prepData);
    rxnScores = groupRxnScores(prepData.minModel, origRxnScores, prepData.refModel.rxns, prepData.groupIds, rxnsToIgnore);

    essentialRxns = prepData.essentialRxns;
    toRev = false(numel(mm.rxns),1);
    %Handle the results from previous steps ('ignore', 'exclude', 'essential')
    if strcmp(stp.HowToUsePrevResults, 'exclude')
        rxnScores(rxnsTurnedOn) = 0; %This is not used anymore in any step setup.
    elseif strcmp(stp.HowToUsePrevResults, 'essential')
        %Make all reversible reactions turned on in previous steps reversible
        %in the direction that they were previously carrying flux

        %first reverse the reactions that need to be reversed
        rev = mm.rev == 1;
        toRev = rxnsTurnedOn & rev & fluxes < 0;
        mm = reverseRxns(mm, mm.rxns(toRev));

        %Then make them irreversible
        mm.rev(rxnsTurnedOn) = 0;
        mm.lb(rxnsTurnedOn) = 0;

        essentialRxns = unique([prepData.essentialRxns;mm.rxns(rxnsTurnedOn)]);
    end

    mipGap = 1;
    first = true;
    success = false;
    fullMipRes = [];
    for rn = 1:length(stp.MILPParams)
        params = stp.MILPParams{rn};
        if ~isfield(params, 'MIPGap')
            params.MIPGap = 0.0004;
        end

        if ~isfield(params, 'TimeLimit')
            params.TimeLimit = 5000;
        end

        if ~first
            %There is sometimes a problem with that the objective function becomes close to zero,
            %which leads to that a small percentage of that (which is the MIPGap sent in) is very small
            %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap,
            %converted to a percentage using the last value of the objective function.
            params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1);
            params.seed = 1234;%use another seed, may work better

            if mipGap <= params.MIPGap
                success = true;
                break; %we're done - this will not happen the first time
            else
                disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)])
            end
        end

        first = false;

        %now run the MILP
        try
            %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small?
            %There could be a need to specify this somewhere in the call at some point.
            %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case
            startVals = [];
            if ~isempty(fullMipRes)
                startVals = fullMipRes.full;
            end
            [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose);
            %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed
            fluxes1(toRev) = -fluxes1(toRev);

            mipGap = fullMipRes.mipgap;
            lastObjVal = fullMipRes.obj;
        catch e
            mipGap = Inf;
            lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play
        end

        success = mipGap <= params.MIPGap;
    end

    if ~success
        dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]);
    end

    %save the reactions turned on and their fluxes for the next step
    rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.';
    %The fluxes are a bit tricky - what if they change direction between the steps?
    %The fluxes are used to determine the direction in which reactions are forced on
    %(to simplify the problem it is good if they are unidirectional).
    %We use the following strategy:
    %1. Use the fluxes from the most recent step.
    %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps
    %This could in theory cause problems, but seems to work well practically
    fluxesOld = fluxes;
    fluxes = fluxes1;
    %make sure that all reactions that are on actually has a flux - otherwise
    %things could go bad, since the flux will be set to essential in a random direction
    %This sometimes happens for rxns with negative score - let's just accept that.
    %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn))
    %    dispEM('There are rxns turned on without flux - this might cause problems');
    %end
    %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9);
end

%get the essential rxns
essential = ismember(prepData.minModel.rxns,prepData.essentialRxns);
%So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions
% (this cannot happen by chance, taken care of in the function groupRxnScores)
rxnsToIgn = rxnScores == 0;
deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential);
deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel);

%Here we need to figure out which original reactions (before the linear merge)
%that were removed. These are all reactions with the same group ids as the removed reactions
groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above
groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group!
rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well

initModel = removeReactions(prepData.refModel,rxnsToRem,false,true);

% remove metabolites separately to avoid removing those needed for tasks
unusedMets = initModel.mets(all(initModel.S == 0,2));
initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks));

%if printReport == true
%    printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype);
%    printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype);
%end

%The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls
%ftINITFillGaps, which automatically removes exchange metabolites (because it
%assumes that the reactions are constrained when appropriate). In this case the
%uptakes/outputs are retrieved from the task sheet instead. To prevent
%exchange reactions being used to fill gaps, they are deleted from the
%reference model here.
initModel.id = 'INITModel';

%If gaps in the model should be filled using a task list
if ~isempty(prepData.taskStruct)
    %Remove exchange reactions and reactions already included in the INIT
    %model
    %We changed strategy and instead include all rxns except the exchange rxns in the ref model
    %But we do keep the exchange rxns that are essential.
    %Let's test to remove all, that should work

    %At this stage the model is fully connected and most of the genes with
    %good scores should have been included. The final gap-filling should
    %take the scores of the genes into account, so that "rather bad"
    %reactions are preferred to "very bad" reactions. However, reactions
    %with positive scores will be included even if they are not connected
    %in the current formulation. Therefore, such reactions will have to be
    %assigned a small negative score instead.
    exchRxns = getExchangeRxns(prepData.refModel);
    refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,true,true);
    exchRxns = getExchangeRxns(initModel);
    initModelNoExc = removeReactions(closeModel(initModel),exchRxns,true,true);

    if useScoresForTasks == true
        %map the rxn scores to the model without exchange rxns
        [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns);
        rxnScores2nd = NaN(length(refModelNoExc.rxns),1);
        rxnScores2nd(ia) = origRxnScores(ib);
        %all(rxnScores2nd == refRxnScores);%should be the same, ok!
        [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose);
    else
        [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose);
    end
    %if printReport == true
    %    printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype);
    %    printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype);
    %end

    addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2));
else
    outModel = initModel;
    addedRxnMat = [];
    addedRxnsForTasks = {};
end

% The model can now perform all the tasks defined in the task list.
model = outModel;

% At this stage the model will contain some exchange reactions but probably
% not all (and maybe zero). This can be inconvenient, so all exchange
% reactions from the reference model are added, except for those which
% involve metabolites that are not in the model.

%Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out
%from the problem is not complete, it doesn't have GRPs etc.
%The logic below is a bit complicated. We identify the reactions that should be removed from the full model as
%reactions that have been removed in the init model except the ones that were added back. In addition, we make
%sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other
%reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can
%be performed also without manipulating the b vector in the model (which is what is done in the gap-filling).
exchRxns = getExchangeRxns(prepData.refModel);
deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns));
outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable

% If requested, attempt to remove negative-score genes from the model,
% depending on their role (isozyme or complex subunit) in each grRule.
% See the "removeLowScoreGenes" function more more details, and to adjust
% any default parameters therein.
if ( removeGenes )
    [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype);
    outModel = removeLowScoreGenes(outModel,geneScores);
end

model = outModel;

end

%This is for printing a summary of a model
function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype)
    [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype);
    rxnS = mean(a);
    geneS = mean(b,'omitnan');
    fprintf([name ':\n']);
    fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
    fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
    fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
    fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
end

function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData)
    rxnsToIgnore = false(length(prepData.toIgnoreExch),1);
    if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end;
    if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end;
    if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end;
    if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end;
    if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end;
    if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end;
    if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end;
    if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end;
end

author: SysBioChalmers

repo: RAVEN


text:
function ret = logistic(input)
    ret = 1 ./ (1 + exp(-input));
end

repo: neural-nets


text:
function [img, mmPerVox, hdr, filename, fullHeader] = loadAnalyze(filename, endianType, nativeTypeFlag, overrideScale)
% function [img, mmPerVox, hdr, filename, fullHeader] = loadAnalyze(analHeaderFileName, ...
%   [endianType], [nativeTypeFlag], [overrideScale]);
%
% Reads Analyze7.5 image data into a matlab array 'img'.
% Certain parts of the Analyze header info are stored in 'hdr'.
%
% If you pass in the (optional) endianType it will use this when
% reading 16bit data files. Defaults to 'ieee-le' (little-endian- intel).
% Another common option is 'ieee-be' (big-endian, ie. non-intel).
% NOTE: as of 2004.02.25, the endianType is automatiacally determined, so
% this parameter is not needed and should go away. It is currently silently
% ignored.
%
% If nativeTypeFlag is 1 (default = 0), then the data are returned in an
% array of the same type as the file data. Eg. int16, uint8, etc.
%
% If overrideScale is passed, the passed value is used to scale the data
% rather than the value in the Analyze file header.  This is useful for
% incorrectly formated Analyze files.
%
% After this routine scales the image data, it sets the scale factor in
% returned hdr.pinfo(1) = 1.  This allows backward compatibility with
% routines that performed the scaling outside of this routine, before
% we folded the scaling in.  We don't see any obvious problem with this
% action, but it does deviate from what you would get if you used the
% routine pair hdr = spm_vol(filename); img = spm_read_vols(hdr).  This
% would result in the returned data being scaled (as here), but the hdr
% would in that case still have the scale factor used.
%
% NOTE: as of 2004.02.25, the hdr.dim(4) value is no longer spm-compatible.
% Spm uses this 4th dimension to store the data type. We now return that
% separately. Thus, we are ready to allow 4-dimensional data files.
%
% REQUIRES:
%   No longer requires spm code. We now use a function adapted from
%   Worsley's fmristat (http://www.math.mcgill.ca/~keith/fmristat/)
%   This function (fmris_read_analyze) should be in our repository.
%
% SEE ALSO: analyze2mrGray, mrGray2Analyze, spm_read_vols
%  In particular, analyze2mrGray provides an example of how to orient the
%  data for mrGray.
% s
% HISTORY:
%   2001.12.13 RFD (bob@white.stanford.edu) wrote it.
%   2003.06.13 DHB (brainard@psych.upenn.edu)  Change comment above,
%       as default for nativeTypeFlag is in fact 0 as coded, not 1 as
%       the previous comment said.
%  2003.06.18  DHB (brainard@psych.upenn.edu)  Multiply read image data
%       by scale factor in header.  Set this scale factor to 1 in returned
%       data.  Also added overrideScale arg.   This change made after
%       consultation with RFD.
%  2003.07.03 RFD Fixed typo-looking bug that prevented the scale factor
%       from getting set to 1.
%  2004.01.06 RFD: Replaced "hdr = spm_vol(filename);" with the relevant
%  lines of code from that spm function. We do this because the mmPerVox
%  value in the Analyze header is more likely to be accurate than the one
%  in the spm space-transform 'mat' file. That scale factor is still
%  available in the transform matrix (hdr.mat). Note that this change will
%  only affect those of us who dare to use the spm-style transform files
%  with our analyze data. (This is a third file, in addition to the .img
%  and .hdr files, that is a matlab 'mat' file containing extra spatial
%  transform info that doesn't fit into the Analyze header, such as a full
%  affine transform.)
%  2004.02.25 RFD: Removed spm-dependent code. We now use a function
%  adapted from FMRISTAT (fmris_read_analyze).
%

disp('Obsolete.  Use analyzeRead');
evalin('caller','mfilename');

if ~exist('filename','var') | isempty(filename)
   [fname, fpath] = uigetfile({'*.img','Analyze 7.5 format (*.img)'}, 'Select analyze volume file...');
   filename = [fpath fname];
   if fname == 0  % user cancelled
      return;
   end
end
%[p,f,ext] = fileparts(filename);

% Fix filename to point to the .img file.
[p,f,e] = fileparts(filename);
if(strcmpi(e,'.hdr'))
    filename = fullfile(p,[f '.img']);
elseif(isempty(e) | ~strcmpi(e,'.img'))
    % This will allow dots in the filename
    filename = [filename '.img'];
end
baseFilename = filename(1:end-4);

% Note that the following applies the scale factor for us.
d = fmris_read_analyze(filename);
img = d.data;
mmPerVox = d.vox(1:3);
filename = d.file_name;
fullHeader = d.hdr;
hdr.fname = d.file_name;
hdr.dim = d.dim;
% build spm-style mat field
d.vox = d.vox(1:3);
%hdr.mat = [[diag(d.vox), [d.origin.*-d.vox-d.vox./2]']; [0 0 0 1]];
hdr.mat = [[diag(d.vox), [d.origin.*-d.vox]']; [0 0 0 1]];
hdr.pinfo = [d.scale;0;d.vox_offset];
hdr.descrip = d.descrip;

% NOTE: no longer needed- fmris_read_analyze automatically determins this.
if (exist('endianType','var') & ~isempty(endianType))
    disp([mfilename ': endianType no longer needed- it will be silently ignored.']);
end
if ~exist('nativeTypeFlag','var') | isempty(nativeTypeFlag)
    nativeTypeFlag = 0;
end
if ~exist('overrideScale','var') | isempty(overrideScale)
    overrideScale = 0;
end

% load an spm-style transform file, if it exists
xformFile = [baseFilename '.mat'];
if(exist(xformFile, 'file'))
    xform = load(xformFile);
    hdr.mat = xform.M;
end

% We could check the file size (d.bytes) against the expected number of bytes.
% d = dir(hdr.fname);
% if(length(d)<1)
%     error (['Could not find image file (',hdr.fname,') for this header file (',analHeaderFileName,').']);
% end

% Override scale factor in header?
% Note that the scale factor has already been applied, so we only have to
% apply and 'ajustment' if a different scale is desired.
if (overrideScale)
    hdr.pinfo(1) = overrideScale;
    img = overrideScale./hdr.pinfo(1).*img;
end

% Since the scale factor has been applied, we set it to 1 so that old code
% will still work.
hdr.pinfo(1) = 1;

if (nativeTypeFlag)
    eval(['img=' d.precision '(img);'], ...
        ['disp(''could not convert to native type "' d.precision '"'');']);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [I] = spm_voice_onsets(Y,FS,C,U)
% Identify intervals containing acoustic energy and post onset minima
% FORMAT [I] = spm_voice_onsets(Y,FS,C,U)
%
% Y    - timeseries
% FS   - sampling frequency
% C    - Convolution kernel [Default: 1/16 sec]
% U    - crossing threshold [Default: 1/8  a.u]

%
% I{i} - cell array of intervals (time bins) containing spectral energy
%
% This routine identifies epochs constaining spectral energy of the power
% envelope, defined as the root mean square (RMS) power. The onset and
% offset of words is evaluated in terms of threshold crossings before and
% after the midpoint of a one second epoch. These are supplemented with
% internal minima (after the spectral peak).
%
% see also: spm_voice_onset.m
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_voice_onsets.m 7750 2019-12-05 17:54:29Z spm $

% find the interval that contains spectral energy
%==========================================================================
global VOX

if nargin < 2, FS = VOX.FS; end                  % onset  threshold
if nargin < 4, U  = 1/8;    end                  % height threshold
if nargin < 3, C  = 1/16;   end                  % smoothing

% identify threshold crossings in power
%--------------------------------------------------------------------------
[G,Y] = spm_voice_check(Y,FS,C);                 % smooth envelope
G     = G/mean(G);                               % normalise to mean
n     = length(G);                               % number of bins

% find zero crossings (of U), for isolated words
%--------------------------------------------------------------------------
j0    = find(G(1:end - 1) < U & G(2:end) > U);
jT    = find(G(1:end - 1) > U & G(2:end) < U);

% find minima
%--------------------------------------------------------------------------
j     = find(diff(G(1:end - 1)) < 0 & diff(G(2:end)) > 0);
if numel(j)
    j = j([FS; diff(j)] > FS/32);                % remove minima < 1/32
    j = j(G(j) < 2);                             % remove minima > 2*mean
end

% find onsets preceded by silence
%--------------------------------------------------------------------------
k0    = [];
for k = 1:numel(j0)
    i = (j0(k) - FS/8):(j0(k) - 1);
    i = fix(i(i < n & i > 1));
    if all(G(i) < U) && j0(k) < FS/2
        k0(end + 1) = j0(k);
    end
end

% find the last onset preceded by silence
%--------------------------------------------------------------------------
if isempty(k0), j0 = []; else, j0 = k0(end); end

% find offsets followed by silence
%--------------------------------------------------------------------------
kT    = [];
for k = 1:numel(jT)
    i = (jT(k) + 1):(jT(k) + FS/8);
    i = fix(i(i < n & i > 1));
    if all(G(i) < U) && jT(k) > FS/2
        kT(end + 1) = jT(k);
    end
end

% find the first offset followed by silence
%--------------------------------------------------------------------------
if isempty(kT), jT = []; else, jT = kT(1); end

% use the first minima in the absence of zero crossings
%--------------------------------------------------------------------------
if isempty(j0)
    try
        j0 = j(1);
    catch
        j0 = 1;
    end
end

% use the last minima in the absence of zero crossings
%--------------------------------------------------------------------------
if isempty(jT)
    try
        jT = j(end);
    catch
        jT = n;
    end
end

% use the last sample if offset precedes onset
%--------------------------------------------------------------------------
if (jT - j0) < 1, jT = n; end

% add internal minima
%--------------------------------------------------------------------------
i   = j(j < jT(end) & j > (j0 + FS/16));
jT  = sort(unique([jT; i]));

% indices of interval containing spectral energy
%--------------------------------------------------------------------------
I     = {};
for i = 1:numel(j0)
    for j = 1:numel(jT)

        % k-th interval
        %------------------------------------------------------------------
        k  = j0(i):jT(j);
        ni = numel(k);

        % retain intervals of plausible length
        %------------------------------------------------------------------
        if ni > FS/16 && ni < FS
            I{end + 1} = k;
        end
    end
end

% sort lengths (longest last), with 3 minima or less
%--------------------------------------------------------------------------
for i = 1:numel(I)
    ni(i) = numel(I{i});
end
[d,j] = sort(ni,'ascend');
I     = I(j(1:min(3,end)));

% graphics(if requested)
%==========================================================================
if ~VOX.onsets
    return
else
    spm_figure('GetWin','onsets'); clf;
end

% timeseries
%--------------------------------------------------------------------------
pst   = (1:n)/FS;
subplot(2,1,1)
plot(pst,Y/max(Y), 'b'),     hold on
plot(pst,G/max(G),':b'),     hold on
plot([1 1]/2,[-1 1],'b'),    hold on
for i = 1:numel(I)
    x = [I{i}(1),I{i}(end),I{i}(end),I{i}(1)]/FS;
    y = [-1,-1,1,1];
    c = spm_softmax(rand(3,1))';
    h = fill(x,y,c);
    set(h,'Facealpha',1/8,'EdgeAlpha',1/8);
end
title('Onsets and offsets','FontSize',16)
xlabel('peristimulus time (seconds)'), spm_axis tight, hold off

% envelope and threshold crossings
%--------------------------------------------------------------------------
subplot(2,1,2)
plot(pst,G,'r'), hold on
plot(pst,0*pst + U,'-.'),         hold on
plot([1 1]/2,[0 max(G)],'b'),     hold on
for i = 1:numel(j0), plot(pst(j0(i)),G(j0(i)),'og'), end
for i = 1:numel(jT), plot(pst(jT(i)),G(jT(i)),'or'), end
title('Spectral envelope','FontSize',16)
xlabel('peristimulus time (secs)'), spm_axis tight, hold off
drawnow, pause(1/4)

% uncomment to play interval
%--------------------------------------------------------------------------
% sound(Y(i),FS)

author: spm


text:
function [grad,err,finaldelta] = gradest(fun,x0)
% gradest: estimate of the gradient vector of an analytical function of n variables
% usage: [grad,err,finaldelta] = gradest(fun,x0)
%
% Uses derivest to provide both derivative estimates
% and error estimates. fun needs not be vectorized.
%
% arguments: (input)
%  fun - analytical function to differentiate. fun must
%        be a function of the vector or array x0.
%
%  x0  - vector location at which to differentiate fun
%        If x0 is an nxm array, then fun is assumed to be
%        a function of n*m variables.
%
% arguments: (output)
%  grad - vector of first partial derivatives of fun.
%        grad will be a row vector of length numel(x0).
%
%  err - vector of error estimates corresponding to
%        each partial derivative in grad.
%
%  finaldelta - vector of final step sizes chosen for
%        each partial derivative.
%
%
% Example:
%  [grad,err] = gradest(@(x) sum(x.^2),[1 2 3])
%  grad =
%      2     4     6
%  err =
%      5.8899e-15    1.178e-14            0
%
%
% Example:
%  At [x,y] = [1,1], compute the numerical gradient
%  of the function sin(x-y) + y*exp(x)
%
%  z = @(xy) sin(diff(xy)) + xy(2)*exp(xy(1))
%
%  [grad,err ] = gradest(z,[1 1])
%  grad =
%       1.7183       3.7183
%  err =
%    7.537e-14   1.1846e-13
%
%
% Example:
%  At the global minimizer (1,1) of the Rosenbrock function,
%  compute the gradient. It should be essentially zero.
%
%  rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%  [g,err] = gradest(rosen,[1 1])
%  g =
%    1.0843e-20            0
%  err =
%    1.9075e-18            0
%
%
% See also: derivest, gradient
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 2/9/2007

% get the size of x0 so we can reshape
% later.
sx = size(x0);

% total number of derivatives we will need to take
nx = numel(x0);

grad = zeros(1,nx);
err = grad;
finaldelta = grad;
for ind = 1:nx
  [grad(ind),err(ind),finaldelta(ind)] = derivest( ...
    @(xi) fun(swapelement(x0,ind,xi)), ...
    x0(ind),'deriv',1,'vectorized','no', ...
    'methodorder',2);
end

end % mainline function end

% =======================================
%      sub-functions
% =======================================
function vec = swapelement(vec,ind,val)
% swaps val as element ind, into the vector vec
vec(ind) = val;

end % sub-function end

author: yorgoon

repo: minimum-snap-geometric-control


text:
function wiggleinterferogram(c,scale,type,norm,range)

% Private method. See ../plot for details.

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% EXTRACT IMAGE
w = get(c,'WAVES');
if isempty(get(w(1),'interferogram_time'))
    error('PLOT INTERFER method can only be used following INTERFEROGRAM function');
end
I.time  = get(w(1),'INTERFEROGRAM_TIME');
I.index = get(w(1),'INTERFEROGRAM_INDEX');
I.CC    = get(w(1),'INTERFEROGRAM_MAXCORR');
I.LL    = get(w(1),'INTERFEROGRAM_LAG');

% SET COLOR SCALE FOR LAG PLOT
I.LL = I.LL ./ range;
f = find((I.LL)<-1); I.LL(f)=-1;      % eliminate outliers
f = find((I.LL)> 1); I.LL(f)= 1;
I.LL = 31.5 * I.LL;
I.LL = round(I.LL + 32.5);                       % shift to positive indices

I.TRANS = I.CC - 0.6;
f = find(I.TRANS<0);
I.TRANS(f) = 0;
I.TRANS = round(0.75*10*I.TRANS);
I.LL = I.LL + 64*I.TRANS;

% PREP PLOT
figure('Color','w','Position',[50 50 850 1100]);
box on; hold on;

% FIX ORDER OF TRACES
ord = 1:get(c,'TRACES');

% GET MEAN TRACE AMPLITUDE FOR SCALING BELOW (when norm = 0)
maxlist = [];
for i = ord
    maxlist(end+1) = max(abs( get(c.W(i),'DATA') ));
end;
normval = mean(maxlist);

% ADD IMAGE
if strncmpi(type,'C',1)
    imagesc(I.time,I.index,I.CC);
    colormap(c,'LTC');
    colorbar;
    hold on;
elseif strncmpi(type,'L',1)
    h = image(I.time,I.index,I.LL);
    colormap(jet);
    %caxis([-1 1]);
    cmap = load('colormap_lag.txt');
    invcmap = 1 - cmap;
    cmap = 1 - [0*invcmap ; 0.33*invcmap ; 0.66*invcmap ; 1*invcmap];
    colormap(cmap);
    hcb = colorbar;
    set(hcb,'YLim',[193 256]);
    set(hcb,'YTick',193+64*[0:.125:1]);
    set(hcb,'YTickLabel',range*[-1:.25:1]);

    hold on;
else
    error('Plot type not recognized.');
end

% SET COLOR

% LOOP THROUGH WAVEFORMS
% tmin =  999999;
% tmax = -999999;
% count = 0;
% for i = ord
% 	count = count + 1;
%     d = get(c.W(i),'DATA');            %%%d = c.w(:,i);
% 	if norm==0
%         d = scale * d/normval;			% do not normalize trace amplitudes
%     else
%         if max(abs(d))==0
%             d = scale * d;              	% ignore zero traces
%         else
%             d = scale * d/max(abs(d));		% normalize trace amplitudes
%         end
%     end
%     d = -1 * d; 				% because scale is reversed below
% 	wstartrel = 86400*(get(c.W(i),'START_MATLAB')-c.trig(i));	% relative start time (trigger is at zero)
% 	tr = wstartrel + [ 0:length(d)-1]'/get(c.W(i),'Fs');
% 	plot(tr,d+count,'k-','LineWidth',1.5);
%     % save min and max relative trace times
% 	if tr(1) < tmin
% 		tmin = tr(1);
% 	end;
% 	if tr(end) > tmax
% 		tmax = tr(end);
% 	end;
%
% end;

% --------------------------------
wstartrel = 86400 *( get(c.W(ord),'START_MATLAB') - c.trig(ord));% relative start time (trigger is at zero)
freq = get(c.W(ord),'Fs');
lengths = get(c.W(ord),'data_length');
tr = nan(max(lengths),numel(ord)); %pre allocate with nan to not plot
abs_max =  max(abs(c.W(ord)));
for count = 1:numel(ord)
    tr(1:lengths(count),ord(count)) = ...
        wstartrel(count) + [ 0:lengths(count)-1]'/freq(count);
end;

% scale is negative because it is reversed below
if norm==0
    % GET MEAN TRACE AMPLITUDE FOR SCALING BELOW (when norm = 0)
    maxlist = max(abs(c.W(ord)));
    normval = mean(maxlist);
    d =  double(c.W(ord) .*( -scale ./ normval)+ [1:numel(ord)]','nan'); % do not normalize trace amplitudes
else
    abs_max(abs_max==0) = 1; % ignore zero traces

    d = double(c.W(ord) .* (-scale ./ abs_max)+[1:numel(ord)]','nan'); % normalize trace amplitudes
end

plot(tr,d,'k-','LineWidth',1.5);
% ------------------------------------------------------
% adjust figure
axis([min(tr(:)) max(tr(:)) 0 length(ord)+1]);
%axis([tmin tmax 0 length(ord)+1]);
set(gca,'YDir','reverse');
set(gca,'YTick',1:length(ord));
set(gca,'YTickLabel',datestr(c.trig(ord)),'FontSize',6);
xlabel('Relative Time,(s)','FontSize',8);

% replace dates with station names if stations are different
if ~check(c,'STA')
    sta  = get(c,'STA');
    chan = get(c,'CHAN');

    for i=1:get(c,'TRACES')
       labels(i) = strcat( sta(i) , '_' , chan(i) );
    end
    set( gca , 'YTick' , [1:1:get(c,'TRACES')] );
    set( gca , 'YTickLabel' , labels );
end

%PRINT OUT FIGURE
set(gcf, 'paperorientation', 'landscape');
set(gcf, 'paperposition', [.25 .25 10.5 8] );
try
if strncmpi(type,'C',1)
    print(gcf, '-depsc2', 'FIG_interferogram_corr.ps')
else
    print(gcf, '-depsc2', 'FIG_interferogram_lag.ps')
end
catch
disp('Warning: Unable to save figure in current directory.');
end

author: geoscience-community-codes

repo: GISMO


text:
function Kss=cp_SSgain(A,Q,H,R,P, mtd)
num_st=size(A,1);

if (mtd==0)     %method to use when dare fails to reach a solution
    % [Beq, yy]=chol(Q,'lower');
    Beq=Q^0.5;

    %basis for observable space
    Ob=obsv(A,H);
    [u,s,v]=svd(Ob'*Ob);
    num_null_obs=length(find(diag(s)<1e-6));
    num_range_obs=num_st-num_null_obs;
    base_null_obs=v(:,(num_range_obs+1):num_st);
    base_range_obs=u(:,1:num_range_obs);
    Wo=[base_null_obs base_range_obs];   %x=Wz

    %basis for controllable and uncontrollable spaces
    Co=ctrb(A,Beq);
    [u,s,v]=svd(Co*Co');
    num_null_ctr=length(find(diag(s)<1e-12));
    num_range_ctr=size(A,1)-num_null_ctr;
    base_null_ctr=v(:,(num_range_ctr+1):size(A,1));
    base_range_ctr=u(:,1:num_range_ctr);
    Wc=[base_null_ctr base_range_ctr];   %x=Wz

    %%%compute the intersection of observable space with cont/uncont spaces
    %obs+ctr space (min real)
    if (num_null_ctr>0)
        mx_a=[base_range_obs base_range_ctr];
        [mx_b,pvt] = rref(mx_a);
        ind=setdiff(1:(num_range_obs+num_range_ctr),pvt);
        if (~isempty(ind))
            Woc=base_range_obs*mx_b(1:num_range_obs,ind);
            Woc*diag(diag(Woc'*Woc).^-0.5); %normalize
            num_base_oc=size(Woc,2);
        else
            Woc=[];
            num_base_oc=0;
        end
        %obs+unctr space
        mx_a=[base_range_obs base_null_ctr];
        [mx_b,pvt] = rref(mx_a);
        ind=setdiff(1:(num_range_obs+num_null_ctr),pvt);
        if (~isempty(ind))
            Wou=base_range_obs*mx_b(1:num_range_obs,ind);
            Wou=diag(diag(Wou'*Wou).^-0.5); %normalize
            num_base_ou=size(Wou,2);
        else
            Wou=[];
            num_base_ou=0;
        end

        %%%overall transformation matrix
        W=[base_null_obs Wou Woc];
    else    %no uncontrollable subspace
        W=[base_null_obs base_range_obs];
        num_base_ou=0;
        num_base_oc=base_range_obs;
    end

    %transform the system
    W_inv=inv(W);
    A_t=W_inv*A*W;
    Beq_t=W_inv*Beq;
    Qeq_t=Beq_t*Beq_t';
    H_t=H*W;

    %partitions
    in1=num_null_obs;
    in2=num_null_obs+num_base_ou;
    in3=num_null_obs+num_base_ou+num_base_oc;

    Au=A_t(1:in1,1:in1);   %unobservable part
    Au_cross=A_t(1:in1,in1+1:end);   %cross effects of observable part on unobservable part
    Aou=A_t(in1+1:in2,in1+1:in2);   %observable but unctrollable part
    Aoc=A_t(in2+1:end,in2+1:end);   %observable and controllable
    Aoc_cross=A_t(in2+1:end,in1+1:in2); %cross effects of ou part on oc part (not used)
    Ao=A_t(in1+1:end,in1+1:end);   %observable part

    Qu=Qeq_t(1:in1,1:in1);
    Qu_cross=Qeq_t(1:in1,in1+1:end);    %cross disturbance between o and uno parts
    Qoc=Qeq_t(in2+1:end,in2+1:end);
    Qo=Qeq_t(in1+1:end,in1+1:end);

    Hou=H_t(:,in1+1:in2);
    Hoc=H_t(:,in2+1:end);
    Ho=H_t(:,in1+1:end);    %observable

    %steady state solution for ubsorvable and controllable parts (including
    %uncont. parts except the unit circle)
    [Poc,clp,Koc]=dare(Aoc',Hoc', Qoc, R);
    Koc=Poc*Hoc'*inv(Hoc*Poc*Hoc'+R);

    %Add SS gain values for the unit circle uncontrollable parts but observable
    %parts
    Ko=[zeros(num_base_ou,size(Koc,2));Koc];
    Po=diag_mat(Poc,1,zeros(num_base_ou));

    % %%Steady state solution for observable part
    % [Po,clp,Ko]=dare(A3',H2', Q3, R);
    % Ko=Po*H2'*inv(H2*Po*H2'+R);

    %%%%Compute the cross gain
    mx_a=Au;
    mx_b=(eye(num_range_obs)-Ho'*inv(Ho*Po*Ho'+R)*Ho*Po)*Ao';
    mx_c=Au_cross*(Po-Ko*Ho*Po)'*Ao'+Qu_cross;
    %%check the multiplication of eigenvalues
    if (isempty(mx_a)) %no unobservable space
        eig_mul=0;
    else
        lef=eig(mx_a);
        rig=eig(mx_b);
        eig_mul=lef*rig';
    end

    if (isempty(find(eig_mul==1,1,'first')))
        Pu_cross = dlyap(mx_a,mx_b,mx_c);
        Ku_cross=Pu_cross*Ho'*inv(Ho*Po*Ho'+R);
        %%overall gain
        Kss=W*[Ku_cross;Ko];
    else %in this case stein equation has no unique solution (solution depends initial P value)
        %Perform standard KF iteration to find a K
        disp('comp_ss_gain: NO SS!!');
        K_pre=zeros(num_st,size(H,1));
        K=P*H'*inv(H*P*H'+R);
        while (max(abs(K(:)-K_pre(:))>1e-7))
            K_pre=K;
            P=(eye(size(A))-K*H)*P;
            P=A*P*A'+Beq*Beq';
            K=P*H'*inv(H*P*H'+R);
        end
        Kss=K;
    end
elseif (mtd==1) %iterations
    K_pre=zeros(num_st,size(H,1));
    K=P*H'*inv(H*P*H'+R);
    in=0;
    while (max(abs(K(:)-K_pre(:))>1e-7))
        K_pre=K;
        P=(eye(size(A))-K*H)*P;
        P=A*P*A'+Q;
        K=P*H'*inv(H*P*H'+R);
        in=in+1;
    end
    disp(['cpSSgain 145: SS iteration =' num2str(in)]);
    Kss=K;
end

repo: ekfmonoslam


text:
function x = triangle_area_sample ( n, v )

%*****************************************************************************80
%
%% TRIANGLE_AREA_SAMPLE samples points from the interior of a triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Greg Turk,
%    Generating Random Points in a Triangle,
%    in Graphics Gems,
%    edited by Andrew Glassner,
%    AP Professional, 1990, pages 24-28.
%
%  Parameters:
%
%    Input, integer N, the number of samples to compute.
%
%    Input, real V(2,3), the vertices.
%
%    Output, real X(2,N), the points.
%
  r = rand ( 2, n );

  c = zeros ( 3, n );

  c(1,1:n) = 1.0                 - sqrt ( r(2,1:n) );
  c(2,1:n) = ( 1.0 - r(1,1:n) ) .* sqrt ( r(2,1:n) );
  c(3,1:n) =         r(1,1:n)   .* sqrt ( r(2,1:n) );

  x = v * c;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [jrt, jst]=jacprojRTS(j, i, rt, xyz, r0, a)
% symbolic projection function Jacobian
% code automatically generated with maple

  qr0=r0(j*4+1:(j+1)*4);

  t1 = (rt(1) ^ 2);
  t2 = (rt(2) ^ 2);
  t3 = (rt(3) ^ 2);
  t5 = sqrt((1 - t1 - t2 - t3));
  t6 = 0.1e1 / t5;
  t7 = t6 * qr0(2);
  t9 = -t7 * rt(1) + qr0(1);
  t11 = t6 * qr0(3);
  t13 = -t11 * rt(1) - qr0(4);
  t15 = t6 * qr0(4);
  t17 = -t15 * rt(1) + qr0(3);
  t19 = -t9 * xyz(1) - t13 * xyz(2) - t17 * xyz(3);
  t24 = -t5 * qr0(2) - qr0(1) * rt(1) - rt(2) * qr0(4) + rt(3) * qr0(3);
  t31 = t5 * qr0(3) + qr0(1) * rt(2) + rt(3) * qr0(2) - rt(1) * qr0(4);
  t37 = t5 * qr0(4) + qr0(1) * rt(3) + rt(1) * qr0(3) - rt(2) * qr0(2);
  t39 = t24 * xyz(1) - t31 * xyz(2) - t37 * xyz(3);
  t41 = t6 * qr0(1);
  t43 = -t41 * rt(1) - qr0(2);
  t48 = t5 * qr0(1) - rt(1) * qr0(2) - rt(2) * qr0(3) - rt(3) * qr0(4);
  t52 = t48 * xyz(1) + t31 * xyz(3) - t37 * xyz(2);
  t57 = t43 * xyz(1) + t13 * xyz(3) - t17 * xyz(2);
  t62 = t43 * xyz(2) + t17 * xyz(1) - t9 * xyz(3);
  t67 = t48 * xyz(2) + t37 * xyz(1) + t24 * xyz(3);
  t72 = t43 * xyz(3) + t9 * xyz(2) - t13 * xyz(1);
  t77 = t48 * xyz(3) - t24 * xyz(2) - t31 * xyz(1);
  t89 = -t19 * t31 - t39 * t13 + t43 * t67 + t48 * t62 + t72 * t24 - t77 * t9 + t57 * t37 + t52 * t17;
  t99 = -t19 * t37 - t39 * t17 + t43 * t77 + t48 * t72 - t57 * t31 - t52 * t13 - t62 * t24 + t67 * t9;
  t106 = -t39 * t37 + t48 * t77 - t52 * t31 - t67 * t24 + rt(6);
  t107 = 0.1e1 / t106;
  t119 = -t39 * t31 + t48 * t67 + t77 * t24 + t52 * t37 + rt(5);
  t123 = t106 ^ 2;
  t124 = 0.1e1 / t123;
  t125 = (a(1) * (t39 * t24 + t48 * t52 - t67 * t37 + t77 * t31 + rt(4)) + a(2) * t119 + a(3) * t106) * t124;
  t129 = -t7 * rt(2) + qr0(4);
  t132 = -t11 * rt(2) + qr0(1);
  t135 = -t15 * rt(2) - qr0(2);
  t137 = -t129 * xyz(1) - t132 * xyz(2) - t135 * xyz(3);
  t141 = -t41 * rt(2) - qr0(3);
  t146 = t141 * xyz(1) + t132 * xyz(3) - t135 * xyz(2);
  t151 = t141 * xyz(2) + t135 * xyz(1) - t129 * xyz(3);
  t157 = t141 * xyz(3) + t129 * xyz(2) - t132 * xyz(1);
  t170 = -t137 * t31 - t39 * t132 + t141 * t67 + t48 * t151 + t157 * t24 - t77 * t129 + t146 * t37 + t52 * t135;
  t180 = -t137 * t37 - t39 * t135 + t141 * t77 + t48 * t157 - t146 * t31 - t52 * t132 - t151 * t24 + t67 * t129;
  t187 = -t7 * rt(3) - qr0(3);
  t190 = -t11 * rt(3) + qr0(2);
  t193 = -t15 * rt(3) + qr0(1);
  t195 = -t187 * xyz(1) - t190 * xyz(2) - t193 * xyz(3);
  t199 = -t41 * rt(3) - qr0(4);
  t204 = t199 * xyz(1) + t190 * xyz(3) - t193 * xyz(2);
  t209 = t199 * xyz(2) + t193 * xyz(1) - t187 * xyz(3);
  t215 = t199 * xyz(3) + t187 * xyz(2) - t190 * xyz(1);
  t228 = -t195 * t31 - t39 * t190 + t199 * t67 + t48 * t209 + t215 * t24 - t77 * t187 + t204 * t37 + t52 * t193;
  t238 = -t195 * t37 - t39 * t193 + t199 * t77 + t48 * t215 - t204 * t31 - t52 * t190 - t209 * t24 + t67 * t187;
  t255 = (a(4) * t119 + a(5) * t106) * t124;
  jrt(1) = (a(1) * (t19 * t24 - t39 * t9 + t43 * t52 + t48 * t57 - t62 * t37 - t67 * t17 + t72 * t31 + t77 * t13) + a(2) * t89 + a(3) * t99) * t107 - t125 * t99;
  jrt(2) = (a(1) * (t137 * t24 - t39 * t129 + t141 * t52 + t48 * t146 - t151 * t37 - t67 * t135 + t157 * t31 + t77 * t132) + a(2) * t170 + a(3) * t180) * t107 - t125 * t180;
  jrt(3) = (a(1) * (t195 * t24 - t39 * t187 + t199 * t52 + t48 * t204 - t209 * t37 - t67 * t193 + t215 * t31 + t77 * t190) + a(2) * t228 + a(3) * t238) * t107 - t125 * t238;
  jrt(4) = a(1) * t107;
  jrt(5) = a(2) * t107;
  jrt(6) = a(3) * t107 - t125;
  jrt(7) = (a(4) * t89 + a(5) * t99) * t107 - t255 * t99;
  jrt(8) = (a(4) * t170 + a(5) * t180) * t107 - t255 * t180;
  jrt(9) = (a(4) * t228 + a(5) * t238) * t107 - t255 * t238;
  jrt(10) = 0.0e0;
  jrt(11) = a(4) * t107;
  jrt(12) = a(5) * t107 - t255;

  t1 = (rt(1) ^ 2);
  t2 = (rt(2) ^ 2);
  t3 = (rt(3) ^ 2);
  t5 = sqrt((1 - t1 - t2 - t3));
  t10 = -t5 * qr0(2) - qr0(1) * rt(1) - rt(2) * qr0(4) + rt(3) * qr0(3);
  t11 = t10 ^ 2;
  t16 = t5 * qr0(1) - rt(1) * qr0(2) - rt(2) * qr0(3) - rt(3) * qr0(4);
  t17 = t16 ^ 2;
  t22 = t5 * qr0(4) + qr0(1) * rt(3) + rt(1) * qr0(3) - rt(2) * qr0(2);
  t28 = -t5 * qr0(3) - qr0(1) * rt(2) - rt(3) * qr0(2) + rt(1) * qr0(4);
  t29 = t28 ^ 2;
  t32 = t10 * t28;
  t35 = -t16 * t22;
  t36 = 0.2e1 * t32 + t16 * t22 - t35;
  t38 = -t10 * t22;
  t39 = t16 * t28;
  t42 = t38 + 0.2e1 * t39 - t10 * t22;
  t48 = t10 * xyz(1) + t28 * xyz(2) - t22 * xyz(3);
  t53 = t16 * xyz(3) - t10 * xyz(2) + t28 * xyz(1);
  t58 = t16 * xyz(1) - t28 * xyz(3) - t22 * xyz(2);
  t63 = t16 * xyz(2) + t22 * xyz(1) + t10 * xyz(3);
  t65 = -t48 * t22 + t16 * t53 + t58 * t28 - t63 * t10 + rt(6);
  t66 = 0.1e1 / t65;
  t78 = t48 * t28 + t16 * t63 + t53 * t10 + t58 * t22 + rt(5);
  t82 = t65 ^ 2;
  t83 = 0.1e1 / t82;
  t84 = (a(1) * (t48 * t10 + t16 * t58 - t63 * t22 - t53 * t28 + rt(4)) + a(2) * t78 + a(3) * t65) * t83;
  t92 = t22 ^ 2;
  t93 = t29 + t17 - t10 ^ 2 - t92;
  t95 = -t28 * t22;
  t98 = t16 * t10;
  t99 = 0.2e1 * t95 - t16 * t10 - t98;
  t111 = t95 + 0.2e1 * t98 - t28 * t22;
  t114 = t92 + t17 - t28 ^ 2 - t11;
  t127 = (a(4) * t78 + a(5) * t65) * t83;
  jst(1) = (a(1) * (t11 + t17 - t22 ^ 2 - t29) + a(2) * t36 + a(3) * t42) * t66 - t84 * t42;
  jst(2) = (a(1) * (t32 + 0.2e1 * t35 + t10 * t28) + a(2) * t93 + a(3) * t99) * t66 - t84 * t99;
  jst(3) = (a(1) * (0.2e1 * t38 - t16 * t28 - t39) + a(2) * t111 + a(3) * t114) * t66 - t84 * t114;
  jst(4) = (a(4) * t36 + a(5) * t42) * t66 - t127 * t42;
  jst(5) = (a(4) * t93 + a(5) * t99) * t66 - t127 * t99;
  jst(6) = (a(4) * t111 + a(5) * t114) * t66 - t127 * t114;

author: vrabaud


text:
%% housekeeping
clear
close all
clc
%% bring in the data
[num,txt] = xlsread('Smets_Wouters_data.xlsx');
vnames=strrep(txt,' ','');
date_start='1947q3';
data=ts(date_start,num(:,2:end),vnames);
data=pages2struct(data);
figure('name','observed data')
for ii=1:numel(vnames)
    subplot(3,3,ii)
    v=vnames{ii};
    plot(data.(v),'linewidth',2)
    title(v)
end
xrotate(45)

%% Set the RFVAR
endo_aliases=struct('robs','interest rates',...
    'dy','GDP growth',...
    'labobs','hours worked',...
    'pinfobs','inflation',...
    'dw','wages growth',...
    'dc','Consumption growth',...
    'dinve','Investment growth');

endog=fieldnames(endo_aliases);

nlags=3;

const=true;

exog={};

rv=rfvar(endog,exog,nlags,const); % formerly redfvar

%% Estimate the reduced-form VAR

rv1=estimate(rv,data);%,{db.LGDP.start,db.LGDP.finish}

%% identification setup

shock_aliases=struct('mp','monetary policy',...
    'ad','aggreg. demand',...
    'as','aggreg. supply',...
    'ls','labor supply');

structural_shocks=fieldnames(shock_aliases);

ident_restr={
    %     'dy{inf}@mp',0
    %     'dy{inf}@ad',0
    %     'dy{inf}@as',0
    %     'dy{inf}@ls',0
    'robs{0}@mp','+'
    'dy@mp','-'
    'pinfobs@mp','-'
    'robs{0}@ad','+'
    'pinfobs@ad','+'
    'dy@ad','+'
    %     'robs@as','-'
    'dy@as','+'
    'pinfobs@as','-'
    'robs@ls','-'
    'dy@ls','+'
    'labobs@ls','+'
    %     'pinfobs@ls','-'
    'dw@ls','-'
    };

agnostic=true;

max_trials=6000;

[Rfunc,ident]=identification(rv1,ident_restr,structural_shocks,...
        agnostic,max_trials);

%% Impulse responses

myirfs=irf(rv1,structural_shocks,40,[],Rfunc);

%% Plot the IRFs for one particular rotation

for ishock=1:numel(structural_shocks)

    shock=structural_shocks{ishock};

    figure('name',['Orthogonalized responses to a ',...
        shock_aliases.(shock),' shock']);

    for ivar=1:numel(endog)

        vname=endog{ivar};

        subplot(3,3,ivar)

        plot(myirfs.(shock).(vname),'linewidth',2)

        title(endo_aliases.(vname))

    end

end

author: jmaih


text:
function [cost, info, x, A] = test22(n)
% Test for the preconditioner support of RCG and RTR.
%
%

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    clc;
%     reset(RandStream.getDefaultStream);
%     randnfoo = randn(123456, 1); %#ok<NASGU>

    if ~exist('n', 'var') || isempty(n)
        n = 3000;
    end

    % Define the problem data : a random symmetric, positive definite
    % matrix with an ill-conditioned diagonal part.
    [Q, ~] = qr(randn(n));
    A = Q*diag(rand(n, 1))*Q';
    A = A + 150*diag(logspace(-3, 3, n));
%     A = -A;

    % Compute a preconditioner for A
%     P = diag(1./diag(A));
%     fprintf('Cond. of  A : %e\n', cond(A));
%     fprintf('Cond. of PA : %e\n', cond(A*P));

%     keyboard;

    % Create the problem structure
    M = spherefactory(n);
    problem.M = M;

    % The cost and gradient
    problem.costgrad = @costgrad;
    function [cost grad store] = costgrad(x, store)

        if ~isfield(store, 'cost')
            Ax = A*x;
            store.cost = -x'*Ax;
        end
        cost = store.cost;

        if ~isfield(store, 'grad')
            store.grad = -2*(Ax + cost*x);
        end
        grad = store.grad;

    end

    % The preconditioner, which is an approximation for the inverse of the
    % Hessian and should be a symmetric, positive definite linear operator.
    % See notes March 18, 2013
    problem.precon = @precon;
    function [Pu store] = precon(x, u, store)
        cost = store.cost;
%         approx_hess = 2*(diag(diag(A)) - cost*eye(n));
        approx_hess = 2*(-cost - diag(A));
        approx_hess(approx_hess < 1e-8) = 1;
%         if any(approx_hess < 0)
%             fprintf('oops\n');
%         end
%         Pu = M.proj(x, approx_hess\u);
        Pu = M.proj(x, u./approx_hess);
%         Pu = (-2*(A+cost*eye(n)))\u;
%         H = 2*(x*x'*A + (x'*A*x)*eye(n) - A)*(eye(n)-x*x');
%         keyboard;
%         Pu = H\u;
%         Pu = u;
%         Pu = M.proj(x, Pu);
%         hist(log10(approx_hess));
%         drawnow;
    end

    % If the optimization algorithms require Hessians, since we do not
    % provide it, it will go for a standard approximation of it. This line
    % tells Matlab not to issue a warning when this happens.
    warning('off', 'manopt:getHessian:approx');

    % Check gradient consistency.
%     checkgradient(problem);

    % Solve
    fprintf('\n ---- With preconditioning ---- \n');

    options = struct();
    options.maxinner = n;
    options.tolgradnorm = 1e-5;
    options.minstepsize = 1e-15;
    options.beta_type = 'H-Z';
    [x cost info] = trustregions(problem, [], options);
    fprintf('Inner work: %d\n', sum([info.numinner]));
    [x cost info] = conjugategradient(problem, [], options);

    fprintf('\n ---- Without preconditioning ---- \n');
    problem = rmfield(problem, 'precon');
    options.maxinner = n;
    [x cost info] = trustregions(problem, [], options);
    fprintf('Inner work: %d\n', sum([info.numinner]));
    [x cost info] = conjugategradient(problem, [], options);

% keyboard;

end

author: NicolasBoumal

repo: manopt


text:
clear all

Nsamples = 41500;
EulerSaved = zeros(Nsamples, 3);

dt = 0.01;

for k=1:Nsamples
  [p q r] = GetGyro();
  A = eye(4) + dt*1/2*[ 0  -p  -q  -r;
                        p   0   r  -q;
                        q  -r   0   p;
                        r   q  -p   0
                      ];

  [ax ay] = GetAccel();
  [phi theta] = EulerAccel(ax, ay);
  z = EulerToQuaternion(phi, theta, 0);

  [phi theta psi] = EulerKalman(A, z);

  EulerSaved(k, :) = [ phi theta psi ];
end

PhiSaved   = EulerSaved(:, 1) * 180/pi;
ThetaSaved = EulerSaved(:, 2) * 180/pi;
PsiSaved   = EulerSaved(:, 3) * 180/pi;

t = 0:dt:Nsamples*dt-dt;

figure
plot(t, PhiSaved)

figure
plot(t, ThetaSaved)

figure
plot(t, PsiSaved)

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function saveresults_data(res,peval, name_res_file)
fprintf('saving data \n');
dir_res = [peval.path_res];
if ~(strcmp(peval.path_res, peval.path_data)) %not identical
    mkdir (dir_res);
end
save ([dir_res '/' name_res_file],'res', 'peval', 'p')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


function info = vl_simplenn_display(net, res)
% VL_SIMPLENN_DISPLAY  Simple CNN statistics
%    VL_SIMPLENN_DISPLAY(NET) prints statistics about the network NET.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

fields={'layer', 'type', 'support', 'stride', 'pad', 'dim', 'fdim', 'field', 'mem'};
if nargin > 1
  fields = {fields{:}, 'xwhd', 'xmem', 'dxmem'} ;
end

do_print = (nargout == 0) ;

for w=fields
  switch char(w)
    case 'type', s = 'type' ;
    case 'stride', s = 'stride' ;
    case 'padding', s = 'pad' ;
    case 'field', s = 'rec. field' ;
    case 'dim', s = 'out dim' ;
    case 'fdim', s = 'filt dim' ;
    case 'mem', s = 'c/g net KB' ;
    case 'xwhd', s = 'x w/h/d' ;
    case 'xmem', s = 'c/g x MB' ;
    case 'dxmem', s = 'c/g dx MB' ;
    otherwise, s = char(w) ;
  end
  if do_print, fprintf('%10s',s) ; end
  for l=1:numel(net.layers)
    ly=net.layers{l} ;
    switch char(w)
      case 'layer', s=sprintf('%d', l) ;
      case 'type'
        switch ly.type
          case 'normalize', s='nrm';
          case 'pool', if strcmpi(ly.method,'avg'), s='apool'; else s='mpool'; end
          case 'conv', s='cnv' ;
          case 'softmax', s='sftm' ;
          case 'loss', s='lloss' ;
          case 'softmaxloss', 'sftml' ;
          otherwise s=ly.type ;
        end
      case 'support'
        switch ly.type
          case 'conv', support(1:2,l) = max([size(ly.filters,1) ; size(ly.filters,2)],1) ;
          case 'pool', support(1:2,l) = ly.pool(:) ;
          otherwise, support(1:2,l) = [1;1] ;
        end
        s=sprintf('%dx%d', support(1,l), support(2,l)) ;
      case 'fdim'
        switch ly.type
          case 'conv'
            filterDimension(l) = size(ly.filters,3) ;
            s=sprintf('%d', filterDimension(l)) ;
          otherwise
            filterDimension(l) = 0 ;
            s='n/a' ;
        end
      case 'stride'
        switch ly.type
          case {'conv', 'pool'}
            if numel(ly.stride) == 1
              stride(1:2,l) = ly.stride ;
            else
              stride(1:2,l) = ly.stride(:) ;
            end
          otherwise, stride(1:2,l)=1 ;
        end
        if all(stride(:,l)==stride(1,l))
          s=sprintf('%d', stride(1,l)) ;
        else
          s=sprintf('%dx%d', stride(1,l), stride(2,l)) ;
        end
      case 'pad'
        switch ly.type
          case {'conv', 'pool'}
            if numel(ly.pad) == 1
              pad(1:4,l) = ly.pad ;
            else
              pad(1:4,l) = ly.pad(:) ;
            end
          otherwise, pad(1:4,l)=0 ;
        end
        if all(pad(:,l)==pad(1,l))
          s=sprintf('%d', pad(1,l)) ;
        else
          s=sprintf('%d,%dx%d,%d', pad(1,l), pad(2,l), pad(3,l), pad(4,l)) ;
        end
      case 'field'
        for i=1:2
          field(i,l) = sum(cumprod([1 stride(i,1:l-1)]).*(support(i,1:l)-1))+1 ;
        end
        if all(field(:,l)==field(1,l))
          s=sprintf('%d', field(1,l)) ;
        else
          s=sprintf('%dx%d', field(1,l), field(2,l)) ;
        end
      case 'mem'
        [a,b] = xmem(ly) ;
        mem(1:2,l) = [a;b] ;
        s=sprintf('%.0f/%.0f', a/1024, b/1024) ;
      case 'dim'
        switch ly.type
          case 'conv', dimension(1,l) = size(ly.filters,4) ;
          otherwise
            if l > 1
              dimension(1,l) = dimension(1,l-1) ;
            end
        end
        s=sprintf('%d', dimension(1,l)) ;
      case 'xwhd'
        for d=1:4, sz(d) = size(res(l+1).x,1) ; end
        s=sprintf('%dx%dx%d%d', sz(1), sz(2), sz(3), sz(4)) ;
      case 'xmem'
        [a,b]=xmem(res(l+1).x) ;
        rmem(1:2,l) = [a;b] ;
        s=sprintf('%.0f/%.0f', a/1024^2, b/1024^2) ;
      case 'dxmem'
        [a,b]=xmem(res(l+1).dzdx) ;
        rmem(1:2,l) = [a;b] ;
        s=sprintf('%.0f/%.0f', a/1024^2, b/1024^2) ;
    end
    if do_print, fprintf('|%7s', s) ; end
  end
  if do_print, fprintf('|\n') ; end
end
if do_print
  [a,b] = xmem(net) ;
  fprintf('total network CPU/GPU memory: %.1f/%1.f MB\n', a/1024^2, b/1024^2) ;
  if nargin > 1
    [a,b] = xmem(res) ;
    fprintf('total result CPU/GPU memory: %.1f/%1.f MB\n', a/1024^2, b/1024^2) ;
  end
end

if nargout > 0
  info.support = support ;
  info.receptiveField = field ;
  info.stride = stride ;
  info.pad = pad ;
end

% -------------------------------------------------------------------------
function [cpuMem,gpuMem] = xmem(s, cpuMem, gpuMem)
% -------------------------------------------------------------------------
if nargin <= 1
  cpuMem = 0 ;
  gpuMem = 0 ;
end
if isstruct(s)
  for f=fieldnames(s)'
    f = char(f) ;
    for i=1:numel(s)
      [cpuMem,gpuMem] = xmem(s(i).(f), cpuMem, gpuMem) ;
    end
  end
elseif iscell(s)
  for i=1:numel(s)
    [cpuMem,gpuMem] = xmem(s{i}, cpuMem, gpuMem) ;
  end
elseif isnumeric(s)
  if isa(s, 'single')
    mult = 4 ;
  else
    mult = 8 ;
  end
  if isa(s,'gpuArray')
    gpuMem = gpuMem + mult * numel(s) ;
  else
    cpuMem = cpuMem + mult * numel(s) ;
  end
end


text:
function a = lt( x, y )

%   Disciplined convex programming information for LT (<):
%      The left-hand side of a less-than constraint must be convex. The
%      right-hand side must be concave. Of course, real constant and
%      affine expressions are both convex and concave and can be used on
%      either side as well.
%
%   Disciplined geometric programming information for LT (<):
%      The left-hand side of a less-than constraint must be log-convex,
%      including positive constants, monomials, posynomials, generalized
%      posynomials, and products thereof. The right-hand side must be
%      log-concave---including positive constants, monomials,
%      reciprocals of log-convex expressions, and products thereof.
%
%   Note that CVX does not distinguish between strict less-than (<) and
%   less-than-or-equal (<=) constraints; they are treated identically.
%   Feasible interior-point solvers tend to return points which satisfy
%   strict inequality, but not all solvers do.

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '<' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function test_min_test04 ( )

%*****************************************************************************80
%
%% TEST_MIN_TEST04 compares the eact and approximate second derivatives.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_MIN_TEST04\n' );
  fprintf ( 1, '  For each problem, compare the eact and\n' );
  fprintf ( 1, '  approximate second derivatives at the starting point.\n' );
%
%  Get the number of problems.
%
  problem_num = p00_problem_num ( );

  for problem = 1 : problem_num

    title = p00_title ( problem );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %d\n', problem );
    fprintf ( 1, '  %s\n', title );

    x = p00_start ( problem );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  X:\n' );
    fprintf ( 1, '  %e\n', x );

    f2 = p00_f2 ( problem, x );

    fprintf ( 1, '  F"(X) (exact):\n' );
    fprintf ( 1, '  %e\n', f2 );

    f2_dif = p00_f2_dif ( problem, x );

    fprintf ( 1, '  F"(X) (difference):\n' );
    fprintf ( 1, '  %e\n', f2_dif );

  end

  return
end


text:
function [out] = Combination(n, k)
% Produces the array of combinations possible picking k from n
% adapted from Combinadics
% http://msdn.microsoft.com/en-us/library/aa289166(VS.71).aspx
%
% USAGE:
%
%    [out] = Combination(n, k)
%
% INPUTS:
%    n:      number of elements in the pool
%    k:      number of elements to pick from `n`
%
% OUTPUT:
%    out:    array of combinations

    if (n < 0 || k < 0) % normally n >= k
      disp('Negative parameter in constructor');
      out = [];
      return
    end

    data(k) = 0;
    for i = 1:k;
      data(i) = i;
    end
   % Combination(n,k)
    out.choose = getNumberCombinations(n,k);

   % determine the combinations
   for i = 1:out.choose
       out.all(:,i) = getCombination(n,k,i-1);
   end
return;

function [m] = getNumberCombinations(n,k)
   % find number of combinations by choose k items from n
    if (n < k)
      m = 0;  % special case
    else
        if (n == k)
            m = 1;
        else
            m = factorial(n) / (factorial(k) * factorial(n-k));
        end
    end

return;

function [c] = getCombination(n,k,index)
    c(1) = 1;
    x = 1;
    for i = 1:n
        if (k == 0)
            return;
        end

        threshold = getNumberCombinations(n-i,k-1);
        %disp(sprintf('index = %d, threshold = %d',index,threshold));
        if (index < threshold)
            c(x) = i;
            x = x+1;
            k = k-1;
        else
            if (index >= threshold)
                index = index - threshold;
            end
        end
    end
return;

author: opencobra

repo: cobratoolbox


text:
function s = resampres(p,m,n);
%RESAMPRES Residual resampling
%
%   Description:
%   S = RESAMPRES(P) returns a new set of indices according to
%   the probabilities P. P is array of probabilities, which are
%   not necessarily normalized, though they must be non-negative,
%   and not all zero. The size of S is the size of P.
%
%   S = RESAMPRES(P,M,N) return M by N matrix.
%
%   S = RESAMPRES(P,M) returns M by M matrix.
%
%   Note that stratified and deterministic resampling have smaller
%   variance.
%
%   Residual resampling retains k_j=floor(n*p_j) copies of index j,
%   where n is length of P and p is normalized probabilities.
%   Additionally n_r=n-k_1-...-k_n indices are sampled with
%   probabilities proportional to n*p_j-k_j. See, e.g., Liu, J. S.,
%   Monte Carlo Strategies in Scientific Computing, Springer, 2001,
%   p. 72.
%
%   See also RESAMPSIM, RESAMPSTR, RESAMPDET
%
% Copyright (c) 2003-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin<2
    [m,n] = size(p);
elseif nargin==2
    n = m;
end
nin=m.*n;
p=p(:);
pn=p./sum(p).*nin;
fpn=floor(pn);
s=zeros(m,n);
k=0;
for i=1:numel(pn)
  if pn(i)>=1
    a=fpn(i);
    pn(i)=pn(i)-a;
    s(k+[1:a])=i;
    k=k+a;
  end
end
pc=cumsum(pn./(nin-k));
pc(end)=1;
s((k+1):nin)=binsgeq(pc,rand(1,nin-k));

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function imdb = Extract_PPG_feautures(PPG,HRVparams,SigName)

% function imdb = PPG_feat_ext(PPG,HRVparams)
% This function is a modified version of CRC_features_calculation used in
% CV_SleepStaging designed to dela with PPG signals insted ECG
%
% Inputs:
%        PPG : vector containing the PPG signals
%  HRVparams : structur with
%    SigName : name use to save derived features
%
% Output:
%       imdb : structur containing the extracted feature to be used for
%              classifictaion
%
% This function is a part of the PPG_sleepstagin project
%         https://github.com/cliffordlab/PPG-sleepstaging
%
% A list of dependencies is given in the readme file
%
% Author : Giulia Da Poian (giulia.dap@gmail.com), 07-Dec-2018
% Part of this code was originally included in the CRC_features_ectract
% in the CV_SleepStaging code, and has been modified to work with PPG
% signals
%

imdb=[];
imdb.images.data=[];
imdb.meta.features=[];

SavePath = HRVparams.writedata;
Fs_resamp = 4;
NFFT = 2000;

% 512 points for 2 Hz (Fs_resamp/2)
very_low_band = 1:3;  % <0.01
low_band = 4:25;      % 0.01-0.1
high_band = 26:103;   % 0.1-0.4

win_inc = HRVparams.increment; % sliding window length is 30s
win_len = HRVparams.windowlength; % segment length 300s (5 min)

% try

   % Preprocess the PPG signal and extract IBI
   if ~exist([SavePath filesep 'FiltSigs' filesep SigName '_filt.mat'],'file')
       ppg_filt = get_filtered_PPG(PPG,HRVparams,1,SigName);
   else
       load([SavePath filesep 'FiltSigs' filesep SigName '_filt.mat'],'ppg_filt');
   end

    % check if annotations are already presents
    if exist([SavePath filesep 'Annotation' filesep SigName '.sqippg'],'file')
        [onsets,~,sqinum] = read_ann([SavePath filesep 'Annotation' filesep SigName],'sqippg');
%%%%%        sqinum = sqinum(2:end); % convert sqi to range 0-100
    else
        % PPG Detection - qppg
        onsets = qppg_fast(ppg_filt,HRVparams.Fs);
        % PPG SQI
        [sqinum, ~, ppgsqi] = calculate_ppgsqi(onsets,ppg_filt,HRVparams.Fs);
        % Write PPG  annotations
        AnnotationFolder = [SavePath filesep 'Annotation'];
        if ~exist(AnnotationFolder,'dir'); mkdir(AnnotationFolder); end
        write_ann([AnnotationFolder filesep SigName],HRVparams,'sqippg',onsets(1:length(sqinum)),char(ppgsqi),sqinum);
    end

    RIAV = Extract_Respiratory_Components(ppg_filt,onsets,HRVparams.Fs); % respiratory induced  amplitude  variation
%%%%% adding sqi_resamp processing
    [NN_resamp, RIAV_resamp, sqi_resamp] = clean_resample_timeseries(onsets,RIAV,sqinum,HRVparams.Fs,Fs_resamp);

    % segmnets of 5-min-length, sliding forwrard every 30s
    N_win = floor(length(RIAV_resamp)/Fs_resamp-win_len)/win_inc+1;

    for idx = 1:N_win
        % image data number

        try
            rr_in_win = NN_resamp((idx-1)*win_inc*Fs_resamp+1:min(((idx-1)*win_inc+win_len)*Fs_resamp,length(NN_resamp)));
            RIAV_in_win = RIAV_resamp((idx-1)*win_inc*Fs_resamp+1:min(((idx-1)*win_inc+win_len)*Fs_resamp,length(RIAV_resamp)));

            rr = rr_in_win;
            rrs = (1/Fs_resamp):(1/Fs_resamp):length(rr_in_win)/Fs_resamp;

            %%% HRV parameters from NN intervals
            SDNN=std(rr);

            % Using HRV_toolbox functions
            [AC,DC] = prsa(rr, rrs, HRVparams, [], 1);
            % SampEn
            sampEn_10 = EvalEntropyMetrics(rr, cumsum(rr),2,0.10, HRVparams, 1, []);
            sampEn_15 = EvalEntropyMetrics(rr, cumsum(rr),2,0.15, HRVparams, 1, []);
            sampEn_20 = EvalEntropyMetrics(rr, cumsum(rr),2,0.20, HRVparams, 1, []);
            % LF/HF
            [PSD,F] = CalcLomb(rrs',detrend(rr'), 1/NFFT:1/NFFT:(Fs_resamp/2) ,[],0);
            LF_HF_ratio=sum(PSD(F>0.04&F<=0.15))/sum(PSD(F>0.15&F<=0.4));

            x = detrend(rr_in_win);
            y = detrend(RIAV_in_win);
            Pxy_amp = abs(csd_amp(x,y,2^10,Fs_resamp,hanning(512),256));
            Pxy_phase = abs(csd_phase(x,y,2^10,Fs_resamp,hanning(512),256));

            % get max amplitude and normalize
            max_amp = max(Pxy_amp);
            max_phase = max(Pxy_phase);
            z1 = (Pxy_amp./max_amp).^2;
            z2 = (Pxy_phase./max_phase).^2;
            Cxy = z1.*z2;

            % call crc % from Robert Joseph Thomas, An Electrocardiogram-Based
            % Technique to Assess Cardiopulmonary Coupling During Sleep

            cross_sp=cpsd(x,y,hanning(512),[],1024);
            co=coherence(x(1:1024),y(1:1024),hanning(512),1024);
            crc = abs(cross_sp).^2 .* co';
            crc=crc./max(crc);

            [low_peaks1, ~]=findpeaks(Cxy(low_band),'SORTSTR','descend');
            [high_peaks1, ~]=findpeaks(Cxy(high_band),'SORTSTR','descend');
            if ~isempty(low_peaks1) && ~isempty(low_peaks1)
                Lo_Hi_ratio1=sum(low_peaks1(1:min(2,length(low_peaks1))))/sum(high_peaks1(1:min(2,length(low_peaks1))));
            else
                if isempty(low_peaks1)
                    Lo_Hi_ratio1=0;
                end
                if isempty(high_peaks1)
                    Lo_Hi_ratio1=NaN;
                end
            end
            [low_peaks2, ~] = findpeaks(crc(low_band), 'SORTSTR' ,'descend');
            [high_peaks2, ~] = findpeaks(crc(high_band), 'SORTSTR' ,'descend');
            if ~isempty(low_peaks2) && ~isempty(low_peaks2)
                Lo_Hi_ratio2=sum(low_peaks2(1:min(2,length(low_peaks2))))/sum(high_peaks2(1:min(2,length(low_peaks2))));
            else
                if isempty(low_peaks2)
                    Lo_Hi_ratio2=0;
                end
                if isempty(high_peaks2)
                    Lo_Hi_ratio2=NaN;
                end
            end

            % area ratio
            sum_very_low=sum(crc(very_low_band));
            sum_low=sum(crc(low_band));
            sum_high=sum(crc(high_band));
            area_ratio_v_l=sum_very_low/sum_low;
            area_ratio_v_h=sum_very_low/sum_high;
            area_ratio_l_h=sum_low/sum_high;

            % set each image by a 512 window and slide forward every 10 seconds
            moving_window=512;
            sliding=10*Fs_resamp;
            for k=1:18
                x=rr_in_win((k-1)*sliding+1:(k-1)*sliding+moving_window);
                y=RIAV_in_win((k-1)*sliding+1:(k-1)*sliding+moving_window);
                x=detrend(x);
                y=detrend(y);

                % call crc % from Robert Joseph Thomas, An Electrocardiogram-Based
                % Technique to Assess Cardiopulmonary Coupling During Sleep

                cross_sp=cpsd(x,y,hanning(256),[],NFFT);
                co=coherence(x,y,hanning(256),NFFT);
                crc = abs(cross_sp).^2 .* co';
                crc = crc./max(crc);
                imagedata = crc(1:250);
                % resample crc data to 50 points
                tmpImg = [imagedata(1:25); resample(imagedata(26:end),25,225)];
                imdb.images.data(:,k,1,idx)= tmpImg;
            end

            % set class of image
            % SQI_i = sqinum((idx-1)*win_inc+1:min(((idx-1)*win_inc+win_len),length(sqinum)));
            SQI_i = nanmean(sqi_resamp((idx-1)*win_inc*Fs_resamp+1:min(((idx-1)*win_inc+win_len)*Fs_resamp,length(sqi_resamp))));

            imdb.meta.sqi(:,idx) = SQI_i;
            imdb.meta.features(idx,1) = AC;
            imdb.meta.features(idx,2) = DC;
            imdb.meta.features(idx,3) = sampEn_10;
            imdb.meta.features(idx,4) = sampEn_15;
            imdb.meta.features(idx,5) = sampEn_20;
            imdb.meta.features(idx,6) = NaN;
            imdb.meta.features(idx,7) = SDNN;
            imdb.meta.features(idx,8) = LF_HF_ratio;
            imdb.meta.features(idx,9) = Lo_Hi_ratio1;
            imdb.meta.features(idx,10) = Lo_Hi_ratio2;
            imdb.meta.features(idx,11) = area_ratio_v_l;
            imdb.meta.features(idx,12) = area_ratio_v_h;
            imdb.meta.features(idx,13) = area_ratio_l_h;
            imdb.meta.rr{idx} = rr;
        catch

            error=1
            imdb.meta.sqi(:,idx) = NaN;
            imdb.meta.features(idx,1) = NaN;
            imdb.meta.features(idx,2) = NaN;
            imdb.meta.features(idx,3) = NaN;
            imdb.meta.features(idx,4) = NaN;
            imdb.meta.features(idx,5) = NaN;
            imdb.meta.features(idx,6) = NaN;
            imdb.meta.features(idx,7) = NaN;
            imdb.meta.features(idx,8) = NaN;
            imdb.meta.features(idx,9) = NaN;
            imdb.meta.features(idx,10) = NaN;
            imdb.meta.features(idx,11) = NaN;
            imdb.meta.features(idx,12) = NaN;
            imdb.meta.features(idx,13) = NaN;
%             imagedata = zeros(1,250);
            imdb.meta.rr{idx} = zeros(1,1200)*NaN;%[imagedata(1:25); resample(imagedata(26:end),25,225)];
            continue;
        end
    end

% catch
% end
imdb.images.data=single(imdb.images.data);

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
if(exist('D:/Datasets/Bosphorus/', 'file'))
    Bosphorus_dir = 'D:\Datasets\Bosphorus/';
elseif(exist('E:/datasets/Bosphorus/', 'file'))
    Bosphorus_dir = 'E:\datasets\Bosphorus/';
else
    fprintf('Bosphorus dataset location not found (or not defined)\n');
end

hog_data_dir = ['E:\Datasets\face_datasets_processed\bosph'];

all_recs = dir([Bosphorus_dir, '/BosphorusDB/BosphorusDB/bs*']);
all_recs_mat = cat(1, all_recs.name);
all_recs = cell(numel(all_recs), 1);

for i=1:size(all_recs_mat,1)

    all_recs{i} = all_recs_mat(i,:);

end

devel_recs = all_recs(1:3:end);
train_recs = setdiff(all_recs, devel_recs);

author: TadasBaltrusaitis

repo: OpenFace


text:
function pass = test_biharm( )

tol = 1e4*chebfunpref().cheb2Prefs.chebfun2eps;

%  Test some spherical harmonics
k = 1;
for ell = [1 2 4 5]
    for m = 0:ell
        f = spherefun.sphharm(ell, m);
        lap2 = biharm(f);
        pass(k, 1) = numel(lap2.pivotValues) == numel(f.pivotValues);
        err(k) = SampleError((ell*(ell+1))^2*f, lap2)/(ell*(ell+1))^2;
        pass(k, 2) = SampleError((ell*(ell+1))^2*f, lap2) < (ell*(ell+1))^2*tol;
        k = k+1;
    end
end
pass = pass(:)';

end

function sample_error = SampleError(h, g)
m = 6;
n = m;
[x, y] = getPoints(m, n);
[L2, T2] = meshgrid(x, y);
F = feval(h, L2, T2);
approx = fevalm(g, x, y);
sample_error = norm(F(:) - approx(:), inf);
end

function [x, y] = getPoints(m, n)

x = trigpts(2*n, [-pi pi]);
y = linspace(0, pi, m).';

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function view=rmSearchFit_sigmaOnly(view,params,matFileName,sigmaExtend);
% rmSearchFit_sigmaOnly - find minimum for retinotopic model per voxel
%
% model=rmSearchFit(view,params,matFileName,sigmaExtend);
%
% Refine sigma estimate over range originalsigma*sigmaExtend
%
% sigmaExtend defaults to [1./3 3];
%
% 2006/04 SOD: wrote it.
% 2006/12 SOD: further optimizations fminsearch->fmincon

% Programming notes:
% This will probably only work for data in 'Gray'-view.

%-----------------------------------
%--- input handling
%-----------------------------------
if ieNotDefined('view'),   error('Need view struct'); end;
if ieNotDefined('params'),
    % See first if they are stored in the view struct
    params = viewGet(view,'rmParams');
    % if not loaded load them:
    if isempty(params),
        view = rmLoadParameters(view);
        params = viewGet(view,'rmParams');
    end;
end;
if ieNotDefined('sigmaExtend'),   sigmaExtend = [1./3 3]; end;

% get rmFile. This is the model definition that will start as a
% starting point for our search.
try,
    rmFile = viewGet(view,'rmFile');
catch,
    disp(sprintf('[%s]:No file selected',mfilename));
    view = rmSelect(view);
    rmFile = viewGet(view,'rmFile');
end;

% Load previous model, but not params since these are allowed to be
% redefined every time. It is the model that is transferable
% between every scan of the same subject but not the exact
% parameters.
% actually we are going to load it so we can use the "grid" to confine our
% nonlinear minimization.
tmp = load(rmFile);
model = tmp.model;

% lastly we need to define the output file name
if ieNotDefined('matFileName'),
    if ieNotDefined('params.matFileName'),
        params.matFileName = ['retModel-',datestr(now,'yyyymmdd'),'-sFit.mat'];
        % else, use existing params.matFileName
    end;
else,
    params.matFileName = matFileName;
end;

% roi check
switch lower(params.wData),
    case {'roi'},

        % if no roi is selected: select one
        if view.selectedROI == 0,
            switch lower(view.viewType),

                case 'inplane',
                    % for inplanes default to gray matter
                    filename = 'gray.mat';
                    try,
                        view   = loadROI(view,filename);
                    catch
                        view   = loadROI(view,filename);
                    end;

                otherwise,
                    % otherwise ask
                    filename = getROIfilename(view);
                    view     = loadROI(view,filename);

            end;
        end;
        ROIcoords = view.ROIs(view.selectedROI).coords;

    otherwise,
        ROIcoords = [];
        % do nothing
end;

%-----------------------------------
%--- make trends to fit with the model (dc, linear and sinewaves)
%-----------------------------------
[trends nTrends] = rmMakeTrends(params);

%-----------------------------------
%--- now loop over slices
%--- but initiate stuff first
%-----------------------------------
switch lower(params.wData),
    case {'fig','roi'},
        loopSlices = 1;
    otherwise,
        % ras 07/07: now uses a separate analysis param, instead of stim
        % param:
        loopSlices = 1:params.analysis.nSlices;

end;
numSlices    = length(loopSlices);

%--- parameters for search fit (fminsearch)
% t-threshold above which to do search. This limit's the search
% algorithms to voxels that will have 'good' data.
params.analysis.fmins.tthresh = 4.42; %p=1e-5; 3.3 = p = 1e-3 (uncorrected)

% display iterations?
params.analysis.fmins.Display = 'none';%'none','iter','final'

% maximum iterations:
params.analysis.fmins.MaxIter = 50;

% Precision of output (degrees). That is, stop if the estimate is
% within TolX degrees:
params.analysis.fmins.TolX    = 1e-2; % degrees

% Precision of evaluation function. We define RMS improvement
% relative to the initial raw 'no-fit' data RMS. So, 1 means
% stop if there is less than 1% improvement on the fit:
params.analysis.fmins.TolFun  = 1e-2; % percent

% fminsearch options
searchOptions.TolX    = params.analysis.fmins.TolX;
searchOptions.MaxIter = params.analysis.fmins.MaxIter;
searchOptions.Display = params.analysis.fmins.Display;
searchOptions.tolFun  = params.analysis.fmins.TolFun;
tthresh       = params.analysis.fmins.tthresh;

% for backward compatibility:
if ~exist('params.analysis.relativeGridStep'),
    params.analysis.relativeGridStep = 1;
end;

% give some feedback so we know we are going
if isempty(ROIcoords),
    fprintf(1,'[%s]:Processing voxels with |t| >= %.2f\n',...
        mfilename,tthresh);
else,
    fprintf(1,'[%s]:Processing voxels with |t| >= %.2f in ROI: %s\n',...
        mfilename,tthresh,view.ROIs(view.selectedROI).name);
end;

% go loop over slices
for slice=loopSlices,

    %-----------------------------------
    % Place datasets behind each other. This is a rather crude way of
    % stimultaneously fitting both. Due to this we cannot
    % prewhiten (we could zeropadd/let the trends deal with this/not care).
    %-----------------------------------
    % we get all the data -
    p2 = params;
    p2.wData = 'all';
    [data, p2] = rmLoadData(view,p2,slice);
    %  [data, params] = rmLoadData(view,params,slice);

    %-----------------------------------
    % now we extract only the data from that slice and put it in a
    % temporary structure that will be modified throughout.
    %-----------------------------------
    s = extractSlice(model,slice);
    % We take the predefined model and remake the params. They can
    % give problems if the 'betas' are different. This may be because
    % there is different amount of detrending and/or different amount
    % of data-sets were used. Anyway, we need to make sure that the
    % amount of 'betas' are set to the current conditions.
    % This means we may have to add more trends to the data:
    if nTrends+1 > size(s{1}.b,1),
        sz = size(s{1}.b,1);
        s{1}.b(sz+1:nTrends+1,:) = 0;
        % or we may have to remove some betas from the model data-set:
    elseif nTrends+1 < size(s{1}.b,1),
        s{1}.b = s{1}.b(1:nTrends+1,:);
    end;
    % sometimes the fit return negative values.

    % amount of models
    nModels  = numel(s);
    % and which ones we are processing (positive only)
    idModels = 1:2:nModels;
    % amount of negative fit
    nNegFit  = 0;

    % for double gaussian model if it exists
    if nModels>2,
        if nTrends+2 > size(s{3}.b,1),
            sz = size(s{3}.b,1);
            s{3}.b(sz+1:nTrends+2,:) = 0;
            % or we may have to remove some betas from the model data-set:
        elseif nTrends+2 < size(s{3}.b,1),
            s{3}.b = s{3}.b(1:nTrends+2,:);
        end;
    end;

    % premake contrasts:
    C{1} = zeros(1,nTrends+1);
    C{1}(1) = 1;
    C{2} = zeros(4,nTrends+2);
    C{2}(:,1:2) = [1 1;1 0;0 1;1 -1];

    %-----------------------------------
    % Now find voxels (|voxel|>tthresh AND in ROI) that will be processed
    %-----------------------------------
    if isempty(ROIcoords),
        wProcess = find(abs(s{1}.t)>=tthresh);
    else,
        %wProcess  = zeros(1,size(ROIcoords,2));
        allcoords = viewGet(view,'coords');
        [tmp, wProcess] = intersectCols(allcoords,ROIcoords);
        wProcess = wProcess(find(s{1}.t(wProcess)>=tthresh));
    end;

    %-----------------------------------
    % Go for each voxel
    %-----------------------------------
    progress = 0;tic;
    for ii = 1:numel(wProcess),

        % progress monitor (10 dots)
        if floor(ii./numel(wProcess)*10)>progress,
            % print out estimated time left
            if slice==1 && progress==0,
                esttime = toc.*10.*numSlices;
                if floor(esttime./3600)>0,
                    fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                        mfilename,numel(wProcess),ceil(esttime./3600));
                else,
                    fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                        mfilename,numel(wProcess),ceil(esttime./60));
                end;
                fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma):',mfilename);
            end;
            fprintf(1,'.');drawnow;
            progress = progress + 1;
        end;

        % volume index
        vi = wProcess(ii);

        % raw rss value
        rawrss     = norm(data(:,vi)).^2;

        % start point from original fit
        xy          = [s{1}.x(vi);...
            s{1}.y(vi)];
        startParams = [s{1}.s(vi)];

        % tight search region [lowerbound upperbound]
        % gridSigmas==startParams(3), somehow this fails sometimes so we'll
        % look for the closest one.
        bndParams      = startParams*sigmaExtend;

        % actual fitting routine
        outParams = ...
            fmincon(@(x) rmModelSearchFit_oneGaussianSigmaOnly(x,data(:,vi),trends,...
            params.analysis.X,...
            params.analysis.Y,...
            params.analysis.allstimimages,...
            rawrss,xy),...
            startParams,[],[],[],[],bndParams(1),bndParams(2),...
            [],searchOptions);

        % compute t
        pred = rfMakePrediction(params,[outParams(1) outParams(1) 0 ...
            xy(1) xy(2)]);
        [t, tmp, rss, b] = rmGLM(data(:,vi),[pred trends],C{1});

        if b(1)>0,
            % store results only if the first beta is positive, somehow fmincon
            % outputs negative fits. If the fit is negative keep old (grid) fit.
            s{1}.s(vi)    = outParams(1);
            s{1}.t(vi)    = t;
            s{1}.rss(vi)  = rss;
            s{1}.rawrss(vi) = rawrss;
            s{1}.b(:,vi)  = b;
        else,
            %disp('Negative beta!?! - ignoring fit');
            nNegFit = nNegFit + 1;
        end

        %--- now for double Gaussian model:
        if numel(s)>2,
            % start point from grid fit
            startParams = [s{3}.x(vi); ...
                s{3}.y(vi); ...
                s{3}.s(vi);...
                s{3}.s2(vi)];

            % tight search region [lowerbound upperbound]
            bndParams   = [-1 1; -1 1; ...
                0.01 params.analysis.sigmaRatioMaxVal;...
                0.1  params.analysis.sigmaRatioInfVal];
            bndParams(1:2,1:2) = startParams(1:2)*[1 1] + ...
                [-1 1;-1 1].*(params.analysis.relativeGridStep.*startParams(3));

            % actual fitting routine
            outParams = ...
                fmincon(@(x) rmModelSearchFit_twoGaussians(x,data(:,vi),trends,...
                params.analysis.X,...
                params.analysis.Y,...
                params.analysis.allstimimages,...
                rawrss),...
                startParams,[],[],[],[],bndParams(:,1),bndParams(:,2),...
                [],searchOptions);
            % compute t
            pred = rfMakePrediction(params,[outParams(3) outParams(3) 0 ...
                outParams(1) outParams(2); ...
                outParams(4) outParams(4) 0 ...
                outParams(1) outParams(2)]);
            [t, tmp, rss, b] = rmGLM(data(:,vi),[pred trends],C{2});

            if b(1)>0 && b(1)+b(2)>0,
                % store results only if the first beta is positive, somehow fmincon
                % outputs negative fits. If the fit is negative keep old (grid) fit.
                s{3}.x(vi)    = outParams(1);
                s{3}.y(vi)    = outParams(2);
                s{3}.s(vi)    = outParams(3);
                s{3}.s2(vi)   = outParams(4);
                s{3}.t(vi)    = t;
                s{3}.rss(vi)  = rss;
                s{3}.b(:,vi)  = b;
            end;
        end;
    end;

    %-----------------------------------
    % now we put back the temporary data from that slice
    %-----------------------------------
    model = putSlice(model,s,slice);

    % end time monitor
    et  = toc;
    if floor(et/3600)>0,
        fprintf(1,'Done [%d hours].\n',ceil(et/3600));
    else,
        fprintf(1,'Done [%d minutes].\n',ceil(et/60));
    end;
    fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
        mfilename,nNegFit,nNegFit./numel(wProcess).*100);
    drawnow;
end;

%-----------------------------------
% save
%-----------------------------------
for n=1:length(model),
    model{n} = rmSet(model{n},'coords',[]);
end;
output = rmSave(view,model,params,1);
view   = viewSet(view,'rmFile',output);

% that's it
return;
%-----------------------------------

%-----------------------------------
function tmp = extractSlice(model,slice);
% some code to extract slice from model struct and place it
% into temporary struct
f = {'x','y','s','rss','t','b','rawrss'};

% loop over models
for n=1:length(model),

    if n==1 || n==2 || n==5, % retinotopy and full field
        f = {'x','y','s','rss','t','rawrss'};
    else,
        if n==3 || n==4,
            f = {'x','y','s','s2','rss','rawrss'};
        else,
            f = {'x','y','s','rss','rawrss'};
        end;
        % put all t values in one matrix
        ts         = {'tall','trm','tf','trmf'};
        tmp{n}.t   = zeros(length(ts),size(val,2));
        for fn = 1:length(ts),
            val            = rmGet(model{n},ts{fn});
            tmp{n}.t(fn,:) = val(slice,:);
        end;
    end;

    % for all models
    tmp{n}.desc = rmGet(model{n},'desc');
    tmp{n}.df   = rmGet(model{n},'dfglm');
    for fn = 1:length(f),
        val    = rmGet(model{n},f{fn});
        if ~isempty(val),
            tmp{n}.(f{fn}) = val(slice,:);
        end;
    end;

    % put all beta values in one matrix
    val      = rmGet(model{n},'b');
    tmp{n}.b = zeros(size(val,3),size(val,2));
    for fn = 1:size(val,3),
        tmp{n}.b(fn,:) = val(slice,:,fn);
    end;
end;

return;
%-----------------------------------

%-----------------------------------
function model = putSlice(model,tmp,slice);
% some code to put slice info into model struct

% loop over models
for n=1:length(model),
    if n<=2, % retinotopy only
        ftmp = {'x','y','s','rss','t','rawrss'};
        fput = {'x','y','s','rss','trm','rawrss'};
        model{n}       = rmSet(model{n},'dfcorr',tmp{n}.df-3); % for x,y,s
    elseif n==5, % full field only
        ftmp = {'x','y','s','rss','t','rawrss'};
        fput = {'x','y','s','rss','tf','rawrss'};
        model{n}       = rmSet(model{n},'dfcorr',tmp{n}.df);
    else, % both
        if n==3 || n==4,
            ftmp = {'x','y','s','s2','rss','rawrss'};
            model{n}       = rmSet(model{n},'dfcorr',tmp{n}.df-4); % for x,y,s,s2
        else,
            ftmp = {'x','y','s','rss','rawrss'};
            model{n}       = rmSet(model{n},'dfcorr',tmp{n}.df-3); % for x,y,s
        end;
        fput = ftmp;

        % distribute t values
        ts  = {'tall','trm','tf','trmf'};
        for fn = 1:4,
            val           = rmGet(model{n},ts{fn});
            val(slice,:)  = tmp{n}.t(fn,:);
            model{n}      = rmSet(model{n},ts{fn},val);
        end;
    end;

    % now get values from model and put in new slice values
    for fn = 1:length(ftmp),
        val             = rmGet(model{n},fput{fn});
        if ~isempty(val),
            val(slice,:)  = tmp{n}.(ftmp{fn});
            model{n}      = rmSet(model{n},fput{fn},val);
        end;
    end;
    model{n}       = rmSet(model{n},'dfglm',tmp{n}.df);

    % distribute beta values
    val      = rmGet(model{n},'b');
    for fn = 1:size(tmp{n}.b,1),
        val(slice,:,fn) = tmp{n}.b(fn,:);
    end;
    model{n} = rmSet(model{n},'b',val);

end;
return;
%-----------------------------------


text:
function [km3] = cl2km3(cl)
% Convert volume from centiliters to cubic kilometers.
% Chad Greene 2012
km3 = cl*1e-14;

author: Sable

repo: mcbench-benchmarks


text:
classdef prtDecisionBinary < prtDecision
    % prtDecisionBinary Base class for all prtDecisionBinary objects
    %
    % A prtDecisionBinary object is an abstract class and cannot be
    % instantiated.
    %
    % prtBinaryDecsion objects find a threshold value that is used to make
    % decisions based on certain criteria.
    %
    % prtDecisionBinary objects all have the following function:
    %
    % getThreshold - return the prtDecisionBinary objects decision
    %                threshold
    %
    % See also: prtDecisionBinaryMinPe, prtDecisionBinarySpecifiedPd,
    % prtDecisionBinarySpecifiedPf, prtDecisionMap

    methods (Abstract)
        threshold = getThreshold(Obj)
        % THRESH = getThreshold returns the objects threshold
    end
    methods
        function obj = prtDecisionBinary()
            obj.classTrain = 'prtDataSetClass';
            obj.classRun = 'prtDataSetStandard';
            obj.classRunRetained = true;

            obj.isSupervised = true;
        end
    end

    methods (Access=protected,Hidden=true)
        function DS = runAction(Obj,DS)
            theClasses = Obj.classList;
            DS = DS.setObservations(theClasses((DS.getObservations >= Obj.getThreshold) + 1));
        end
         function xOut = runActionFast(Obj,xIn,ds) %#ok<INUSD>
            theClasses = Obj.classList;
            xOut = theClasses((xIn >= Obj.getThreshold) + 1);
        end
    end

    methods (Access = protected, Hidden = true)
        function ClassObj = preTrainProcessing(ClassObj, DataSet)
            % Overload preTrainProcessing() so that we can determine mary
            % output status
            assert(DataSet.isLabeled & DataSet.nClasses > 1,'The prtDataSetClass input to the train() method of a prtDecisionBinary must have non-empty targets and have more than one class.');

            ClassObj = preTrainProcessing@prtAction(ClassObj,DataSet);
        end
    end

    methods (Hidden)
        function str = exportSimpleText(self) %#ok<MANU>
            titleText = sprintf('%% prtDecisionBinaryMinPe\n');
            decisionThresholdText = prtUtilMatrixToText(self.threshold,'varName','minPeDecisionThreshold');
            str = sprintf('%s%s%s',titleText,decisionThresholdText);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT