text:
% logimagesc() - make an imagesc(0) plot with log y-axis values (ala semilogy())
%
% Usage:  >> [logfreqs,dataout] = logimagesc(times,freqs,data);
%
% Input:
%   times = vector of x-axis values
%   freqs = vector of y-axis values
%   data  = matrix of size (freqs,times)
%
% Optional Input:
%   plot = ['on'|'off'] plot image or return output (default 'on').
%
% Note: Entering text() onto the image requires specifying (x,log(y)).

% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 4/2000

% Copyright (C) 4/2000 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% 08-07-00 made ydir normal -sm
% 01-25-02 reformated help & license -ad

function [lgfreqs,datout, h, yt, yl] = logimagesc(times,freqs,data,varargin)

  if nargin < 1
      help logimagesc;
      return
  end;
  if size(data,1) ~= length(freqs)
      fprintf('logfreq(): data matrix must have %d rows!\n',length(freqs));
      datout = data;
      return
  end
  if size(data,2) ~= length(times)
      fprintf('logfreq(): data matrix must have %d columns!\n',length(times));
      datout = data;
      return
  end
  if min(freqs)<= 0
      fprintf('logfreq(): frequencies must be > 0!\n');
      datout = data;
      return
  end

  try, icadefs; catch, warning('Using MATLAB default colormap'); end

  lfreqs = log(freqs);
  lgfreqs = linspace(lfreqs(1),lfreqs(end),length(lfreqs));
  lgfreqs = lgfreqs(:);
  lfreqs = lfreqs(:);
  [mesht meshf] = meshgrid(times,lfreqs);
  try
      datout = griddata(mesht,meshf,double(data),times,lgfreqs);
  catch
      fprintf('error in logimagesc.m calling griddata.m, trying v4 method.');
      datout = griddata(mesht,meshf,data,times,lgfreqs,'v4');
  end
  datout(find(isnan(datout(:)))) = 0;

  if ~isempty(varargin)
      plot = varargin{2};
  else
      plot = 'on';
  end

  if strcmp(plot, 'on')
      imagesc(times,freqs,data);
      try colormap(DEFAULT_COLORMAP); catch, end;
      nt = ceil(min(freqs)); % new tick - round up min y to int
      ht = floor(max(freqs)); % high freq - round down

      yt=get(gca,'ytick');
      yl=get(gca,'yticklabel');

      h=imagesc(times,lgfreqs,datout); % plot the image
      set(gca,'ydir','normal')

      i = 0; yt = [];
      yl = cell(1,100);

      tickscale = 1.618; % log scaling power for frequency ticks
      while (nt*tickscale^i < ht )
        yt = [yt log(round(nt*tickscale^i))];
        yl{i+1}=int2str(round(nt*tickscale^i));
        i=i+1;
      end

      if ht/(nt*tickscale^(i-1)) > 1.35
         yt = [yt log(ht)];
         yl{i+1} = ht;
      else
         i=i-1;
      end
     yl = {yl{1:i+1}};
     set(gca,'ytick',yt);
     set(gca,'yticklabel',yl);

%     if nt > min(yt),
%         set(gca,'ytick',log([nt yt]));
%         set(gca,'yticklabel',{int2str(nt) yl});
%      else
%         set(gca,'ytick',log([yt]));
%         set(gca,'yticklabel',{yl});
%      end

  end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function net = vgg_s(labelNum,lossWeight,opts)
channel_num=512;
mag=0.1;
netname='vgg-s';

partRate=1;
textureRate=0;
output_num=labelNum;
net=load('../nets/imagenet-vgg-s.mat');

net.layers=net.layers(1:end-7);
net.layers{end-1}.learningRate=[1,5];
net=add_block_mask(net,lossWeight, opts,'6',3,3,512,channel_num,1,1,partRate,textureRate,mag); %%%%%%%%%%%%%%%%%
net.layers{end+1} = struct('type', 'pool', 'name', 'pool6', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 3, ...
                           'pad', [0 1 0 1]) ;
net=add_block_mask(net,lossWeight,opts,'7',6,6,channel_num,4096,1,0,partRate,textureRate,mag); %%%%%%%%%%%%%%%%%

net = add_dropout(net, opts, '7') ;
net = add_block(net, opts, '8', 1, 1, 4096, 4096, 1, 0) ;
net = add_dropout(net, opts, '8') ;

net = add_block(net, opts, '9', 1, 1, 4096, output_num, 1, 0) ;
net.layers(end) = [] ;
net.meta.classes.name={'pos'};
net.meta.classes.description={'pos'};
net.meta.normalization.imageSize=net.meta.normalization.imageSize(1:3);
net.meta.netname=netname;
end

repo: interpretableCNN


text:
function genes = drosFindTargets(droschip),

% DROSFINDTARGETS returns a sorted list of targets
%
% Usage:
%   genes = drosFindTargets(droschip)
% COPYRIGHT : Antti Honkela, 2008

% SHEFFIELDML

[foo, I] = sort(-sum(droschip.data ~= 0, 2));
genes = droschip.genes(I);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [] = shapeSIRFS(class,jobID)
%GETDEPTHIMAGE Summary of this function goes here
%   Detailed explanation goes here

statesDir = jobDirs(class,jobID,'state');
sirfsstatedir = jobDirs(class,jobID,'sirfs');
dmapDir = jobDirs(class,jobID,'dmap');
mkdirOptional(sirfsstatedir);
fnames = getFileNamesFromDirectory(dmapDir,'types',{'.mat'});

fprintf('\n%%%%%%%%%%%% SIRFS %%%%%%%%%%%%\n');
%% Loading precomputed shape
for i=1:length(fnames)
%for i=1:length(fnames)
    sirfsdmapFile = fullfile(sirfsstatedir,fnames{i});
    if(exist(sirfsdmapFile,'file'))
        continue;
    end
    stateFile = fullfile(statesDir,fnames{i});
    state = load(stateFile);state=state.state;
    dmapFile = fullfile(dmapDir,fnames{i});
    dmap = load(dmapFile); dmap = dmap.dmap;
    sirfsstate = getsirfsstate(state,dmap);
    savefunc(sirfsdmapFile,sirfsstate);
    %vis_PSIRFS(sirfsstate); pause
end
end

function savefunc(file,state)
    save(file,'state');
end

function sirfsstate = getsirfsstate(state,depthIm)
    sigma = 3; % Controls the "bandwidth" of the input depth (the standard deviation of a Gaussian at which point the signal becomes reliable)
    mult = 5; % Controls the importance of the input depth (the multiplier on the loss)
    niters = 200;
    depthIm(isinf(depthIm)) = nan;
    depthIm = -depthIm;
    input_image = im2double(state.im);
    input_image(input_image<1/255) = 1/255;
    dmapMask = ~isinf(depthIm) & ~isnan(depthIm);
    %imshow(color_seg(dmapMask,state.im));
    sirfsstate = SIRFS(input_image, dmapMask, depthIm, ...
        ['params.DO_DISPLAY = 0; params.N_ITERS_OPTIMIZE = ' num2str(niters) ';params.USE_INIT_Z = true; params.INIT_Z_SIGMA = ',...
        num2str(sigma), ';params.multipliers.height.init = { ', num2str(mult), ' };']);
%     sirfsstate = SIRFS(input_image, (state.mask), [], ...
%         ['params.DO_DISPLAY = 0; params.N_ITERS_OPTIMIZE = ' num2str(niters) ';']);
    sirfsstate.mask = dmapMask;
    sirfsstate.im = input_image;
    sirfsstate.dmap = sirfsstate.height;
end

repo: CategoryShapes


text:
function [Ec,Rad,Diam,Cv,Pv]=grEccentricity(E)
% Function Ec=grEccentricity(E) find the (weighted)
% eccentricity of all vertexes of graph.
% Input parameter:
%   E(m,2) or (m,3) - the edges of graph and their weight;
%     1st and 2nd elements of each row is numbers of vertexes;
%     3rd elements of each row is weight of arrow;
%     m - number of arrows.
%     If we set the array E(m,2), then all weights is 1.
% Output parameter:
%   Ec(1,n) - the (weighted) eccentricity of all vertexes.
% [Ec,Rad,Diam]=grEccentricity(E) find also the radius Rad and
%   diameter Diam of the graph (the scalars).
% [Ec,Rad,Diam,Cv,Pv]=grEccentricity(E) find also
%   the center vertexes Cv and the periphery vertexes Pv
%   of the graph (the vector-rows with numbers of vertexes).
% Author: Sergiy Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru

if nargin<1,
  error('There are no input data!')
end
dSP=grDistances(E); % the matrix of distances
Ec=max(dSP); % the eccentricity of all vertexes
Rad=min(Ec); % the radius
Diam=max(Ec); % the diameter
Cv=find(Ec==Rad); % the center vertexes of the graph
Pv=find(Ec==Diam); % the periphery vertexes of the graph
return

author: HuangCongQing


text:
function sir = SIR(Ori,Est)

  s = Ori(:);

  t = Est(:);

  sir = 10*log10( (s'*s) / ((s-t)'*(s-t)) );

author: andrewssobral


text:
function [modelJoint] = createMultipleSpeciesModel(models, varargin)
% Based on the implementation from *Klitgord and Segre 2010, PMID 21124952*.
% The present implementation has been used in *PMID 23022739*, *PMID 25841013*,
% *PMID 25901891*, *PMID 27893703*.
%
% Joins one or more COBRA models with or without another COBRA model
% representing the host. The created setup when a host is entered is
% depicted schematically in Figures 1 and 2 in *PMID 27893703*.
%
% Creates a common space u (lumen) through which all cells can feed and exchange metabolites,
% and separate extracellular spaces for all joined models.
% If a host model is entered, a separate compartment `b` (body fluids) which has no
% connection to extracellular space is created for the host. Metabolites can be transported from
% the lumen to `e`, from `e` to the cytosol and from the cytosol to `b`,
% but not from body fluids to the lumen.
%
% USAGE:
%
%    [modelJoint] = createMultipleSpeciesModel(models, varargin)
%
% INPUTS:
%    models:            cell array of COBRA models(at least one).
%                       Format
%
%                         * models{1,1} = model 1
%                         * models{2,1} = model 2...
%
% OPTIONAL INPUTS:
%    nameTagsModels:    cell array of tags for reaction/metabolite abbreviation
%                       corresponding to each model.
%                       Format
%                         * nameTagsModels{1,1} = 'name tag 1'
%                         * nameTagsModels{2,1} = 'name tag 2'...
%    modelHost:         COBRA model for host
%    nameTagHost:       string of tag for reaction/metabolite abbreviation of host model
%    mergeGenesFlag:    If true, the gene associations in both models are
%                       included in the joined model. If false, empty fields are created
%                       instead (default:false). Note: merging genes is time-consuming
%                       and may crash certain models.
%    remCyclesFlag:     If true, the function will attempt to remove futile
%                       cycles that appear after merging the models (default: false).
%
% OUTPUT:
%    modelJoint:        model structure for joint model
%
% .. Authors:
%       - Ines Thiele and Almut Heinken, 2011-2018
%       - Almut Heinken, 07.02.2018-included option whether or not genes are
%         merged
%       - Almut Heinken, 21.02.2018-fixed compatibility issue with reconstructions
%         from BIGG Models database that have _e instead of [e] as compartment IDs
%       - Almut Heinken, 06.03.2018-changed to parameter-input pairs
%       - Laurent Heirendt, 16/3/2018 - backward compatibility
%       - Almut Heinken, 15.01.2019-fixed compatibility issue with reconstructions
%         from KBase database that have [e0] instead of [e] as compartment IDs
%       - Almut Heinken, 01/2020-remove futile cycles in multi-species
%       AGORA models
%       - Almut Heinken, 07/2022-made futile cycles removal optional
%
% NOTE:
%    This function assumes, that exchange reactions are identified by
%    containing 'EX' in the reaction name and that no other reactions do have this property!

oldOptionalOrder = {'nameTagsModels', 'modelHost', 'nameTagHost', 'mergeGenesFlag'};

% ensure backward compatibility
if numel(varargin) > 0 && ~ischar(varargin{1})
    tempargin = cell(0);
    for i = 1:numel(varargin)
        if ~isempty(varargin{i})
            tempargin{end + 1} = oldOptionalOrder{i};
            tempargin{end + 1} = varargin{i};
        end
    end
    varargin = tempargin;
end

% Define default input parameters if not specified
parser = inputParser();
parser.addRequired('models', @iscell);
parser.addParameter('nameTagsModels', {}, @iscell);
parser.addParameter('modelHost', {}, @isstruct);
parser.addParameter('nameTagHost', '', @(x) ischar(x) || iscell(x))
parser.addParameter('mergeGenesFlag', false, @(x) isnumeric(x) || islogical(x))
parser.addParameter('remCyclesFlag', false, @(x) isnumeric(x) || islogical(x))

parser.parse(models, varargin{:});

models = parser.Results.models;
nameTagsModels = parser.Results.nameTagsModels;
modelHost = parser.Results.modelHost;
nameTagHost = parser.Results.nameTagHost;
mergeGenesFlag = parser.Results.mergeGenesFlag;
remCyclesFlag = parser.Results.remCyclesFlag;

if isempty(models)
   error('Please enter at least one model!')
end
% prepare the model structures and assign name tags for each model structure if not provided
modelNumber = size(models, 1);

if isempty(nameTagsModels)
    % assign default name tags for microbes
    for i = 1:modelNumber
        nameTagsModels{i, 1} = strcat('model', num2str(i), '_');
    end
else
    if size(nameTagsModels, 1) ~= modelNumber
        error('Number of name tags and joint models needs to be identical!')
    end
end

if isempty(nameTagHost)
    % assign default name tag for host
    nameTagHost = 'Host_';
end

%% ensure compatibility with reconstructions from BIGG Models database
for i = 1:modelNumber
    model = models{i, 1};
    metIndices =~cellfun(@isempty, regexp(model.mets, '_e$'));
    model.mets(metIndices) = strrep(model.mets(metIndices), '_e', '[e]');
    models{i, 1} = model;
end
if ~isempty(modelHost)
metIndices =~cellfun(@isempty, regexp(modelHost.mets, '_e$'));
modelHost.mets(metIndices) = strrep(modelHost.mets(metIndices), '_e', '[e]');
end

%% Ensure compatibility with reconstructions from KBase database
for i = 1:modelNumber
    model = models{i, 1};
    metIndices =~cellfun(@isempty, regexp(model.mets, '\[e0\]$'));
    model.mets(metIndices) = strrep(model.mets(metIndices), '[e0]', '[e]');
    % need workaround for biomass metabolite, otherwise the resulting joint
    % model will be unable to carry biomass flux
    if ~isempty(find(ismember(model.mets, 'cpd11416[c0]')))
        model = addDemandReaction(model, 'cpd11416[c0]');
    end
    models{i, 1} = model;
end
if ~isempty(modelHost)
metIndices =~cellfun(@isempty, regexp(modelHost.mets, '\[e0\]$'));
modelHost.mets(metIndices) = strrep(modelHost.mets(metIndices), '[e0]', '[e]');
end
%% Remove futile cycles that may appear after joining (optional)

if remCyclesFlag
    database = loadVMHDatabase;

    unionRxns={};
    for i = 1:modelNumber
        model=models{i, 1};
        unionRxns=vertcat(unionRxns,model.rxns);
    end
    unionRxns=unique(unionRxns);
    for i = 1:modelNumber
        model=models{i, 1};
        biomassReaction=model.rxns(find(strncmp(model.rxns, 'bio', 3)));
        model = removeFutileCycles(model, biomassReaction, database,unionRxns);
        models{i, 1} = model;
    end
end
%% define some variables
eTag = 'u';
exTag = 'e';
% find exchange reactions for models, but leaves demand and sink reactions
% do this for all models to be added, remove exchange reactions from host while leaving demand and sink reactions
% First, find the minimal number of fields common to all models.
presentinallModels = fieldnames(models{1});
missingFields = {};
for i = 2:modelNumber
    cfields = fieldnames(models{i});
    missingFields = union(missingFields, setxor(cfields, presentinallModels));
    presentinallModels = intersect(presentinallModels, cfields);
end
fprintf('The following fields are missing in several models, they will not be merged:\n');
disp(missingFields);
models = restrictModelsToFields(models, presentinallModels);

modelStorage = cell(modelNumber, 1);
for i = 1:modelNumber
    % a new model each turn
    model = models{i, 1};
    % find exchange reactions and external metabolites
    exmod = model.rxns(strmatch('EX', model.rxns));
    % remove all previously defined exchange reactions
    model = removeRxns(model, exmod);
    % make sure the exchange reactions and changed model are saved under correct name
    modelStorage{i, 1} = model;
end

if ~isempty(modelHost)
    %% with a host
    exmod = modelHost.rxns(strmatch('EX', modelHost.rxns));

    % modelHost = removeRxns(modelHost,ExRH);
    % ExRH = modelHost.rxns(selExcH);
    % ExRH(strmatch('sink',modelHost.rxns(selExcH)))=[];
    % ExRH(strmatch('DM',modelHost.rxns(selExcH)))=[];

    % create a new extracellular space for host
    % find all metabolites in e
    relMetIndex = cellfun(@(x) ~isempty(strfind(x, '[e]')), modelHost.mets);
    relMets = modelHost.mets(relMetIndex);
    relRxns = findRxnsFromMets(modelHost, relMets);  % These are all reactions which are relevant
    rxnIndices = ismember(modelHost.rxns, relRxns);
    Stoich = modelHost.S(:, rxnIndices);
    changedMets = regexprep(modelHost.mets, '\[e\]', '\[b\]');
    modelHost = addMultipleMetabolites(modelHost, setdiff(changedMets, modelHost.mets));
    modelHost = addMultipleReactions(modelHost, strcat(modelHost.rxns(rxnIndices), 'b'), changedMets, Stoich, 'lb', modelHost.lb(rxnIndices), ...
                                    'ub', modelHost.ub(rxnIndices), 'c', modelHost.c(rxnIndices), 'subSystems', repmat({'Host Exchange'}, numel(relRxns), 1));

    % remove exchange reactions from host while leaving demand and sink reactions
    modelHost = removeRxns(modelHost, exmod);

    %% create intercellular space
    % will need to find all extracellular metabolites and duplicate reactions
    % if a host model was input, create the shared compartment for the microbes
    model = modelStorage{1, 1};
    nameTag = nameTagsModels{1, 1};
    [modelJoint, MexGJoint] = createInterSpace(model, nameTag, eTag, exTag);

    % if more than one microbe was input
    if modelNumber > 1
        for i = 1:modelNumber
            model = modelStorage{i, 1};
            nameTag = nameTagsModels{i, 1};
            [model, MexG] = createInterSpace(model, nameTag, eTag, exTag);
            % make sure the changed model is saved under correct name
            modelStorage{i, 1} = model;
            MexGJoint = union(MexG, MexGJoint);
        end
     end

    [modelHost,MexGHost] = createInterSpace(modelHost, nameTagHost, eTag, exTag);
    MexGJoint = union(MexGJoint, MexGHost);

    %% merge the models
    % if more than one microbe was input
    if modelNumber > 1
        for i = 2:modelNumber
            model = modelStorage{i, 1};
            [modelJoint] = mergeTwoModels(modelJoint, model, 1, mergeGenesFlag);
        end
    end
    [modelJoint] = mergeTwoModels(modelJoint,modelHost, 1, mergeGenesFlag);

    modelJoint = addExchangeRxn(modelJoint, unique(MexGJoint));

else
    %% without a host
     % create the shared compartment for the microbes
    model = modelStorage{1, 1};
    nameTag = nameTagsModels{1, 1};
    [modelJoint, MexGJoint] = createInterSpace(model, nameTag, eTag, exTag);
     if modelNumber > 1
        for i = 2:modelNumber
            model = modelStorage{i, 1};
            nameTag = nameTagsModels{i, 1};
            [model, MexG] = createInterSpace(model, nameTag, eTag, exTag);
            modelStorage{i, 1}=model;
            MexGJoint = union(MexG, MexGJoint);
        end
     end

     if modelNumber > 1
        for i = 2:modelNumber
            model = modelStorage{i, 1};
            [modelJoint] = mergeTwoModels(modelJoint, model, 1, mergeGenesFlag);
        end
     end

    modelJoint = addExchangeRxn(modelJoint, unique(MexGJoint));
end
end

%%
function [modelNew,MexG] = createInterSpace(model, nameTag, eTag, exTag)
% create intercellular space
% will need to find all extracellular metabolites and duplicate reactions using them
modelNew = model;
% add name tag to all metabolites and reactions in model
modelNew.mets = strcat(nameTag, model.mets);
modelNew.rxns = strcat(nameTag, model.rxns);
% Get the relevant metabolites
relMetsIndex = cellfun(@(x) ~isempty(strfind(x,'biomass[c]')) || ~isempty(strfind(x,['[', exTag, ']'])),modelNew.mets);
relMets = modelNew.mets(relMetsIndex);
% Define the names of the interspace metabolites
MexG = regexprep(strrep(relMets,nameTag,''),strcat('\[', exTag, '\]'), strcat('\[', eTag, '\]'));
varinput = {};
% Add metNames and metFormulas, if present in the original model
if isfield(modelNew,'metNames')
    varinput{end+1} = 'metNames';
    varinput{end+1} = modelNew.metNames(relMetsIndex);
end
if isfield(modelNew,'metFormulas')
    varinput{end+1} = 'metFormulas';
    varinput{end+1} = modelNew.metFormulas(relMetsIndex);
end
% Add all new Metabolites
modelNew = addMultipleMetabolites(modelNew, MexG, varinput{:});

nExchange = numel(relMets);
% Set the exchanger Stoichiometries (met[e] -> met[u])
stoich = [-speye(nExchange);speye(nExchange)];
% Set the names of the exchangers
rxnNames =  strcat(nameTag, 'IEX_', MexG, 'tr');
% Set the bounds
lbs = repmat(-1000,nExchange,1);
ubs = repmat(1000,nExchange,1);
% Set the subSystem
subSystems = repmat({'Transport, intercellular'},nExchange,1);
% Add all Reactions in one go.
modelNew = addMultipleReactions(modelNew,rxnNames,[relMets;MexG],stoich,'lb',lbs,'ub',ubs,'subSystems',subSystems);

end

author: opencobra

repo: cobratoolbox


text:
function [quantized] = quantizeContinuous(values, min_v, max_v, num_bins)

    step_size = (max_v - min_v) / num_bins;
    bin_centres = min_v + step_size/2 : step_size : max_v;

    bin_centres = repmat(bin_centres, numel(values),1);
    values = repmat(values, 1, num_bins);

    [~,quantized] = min(abs(values - bin_centres)');
    quantized = quantized';
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function[struct]=structindex(struct,index)
%STRUCTINDEX  Applies an index to all array-valued fields in a structure.
%
%   STRUCTOUT=STRUCTINDEX(STRUCT,INDEX), where STRUCT is a structure
%   containing a set of fields that are non-string arrays all having the
%   same length N>1, applies INDEX to all of them and returns the result in
%   STRUCTOUT.  All other fields of STRUCT are left unchanged.
%
%   STRUCTINDEX(STRUCT,BOOL) where BOOL is a boolean array of length N also
%   works.
%
%   STRUCTINDEX is useful for dealing with the ridge structures output by
%   EDDYRIDGES.
%
%   See also EDDYRIDGES.
%
%   Usage: struct=structindex(struct,index);
%          struct=structindex(struct,bool);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2019 J.M. Lilly --- type 'help jlab_license' for details

names=fieldnames(struct);
for i=1:length(names)
    temp=getfield(struct,names{i});
    if length(temp)~=1&&~ischar(temp)
        %names{i}
        temp=temp(index);
        struct=setfield(struct,names{i},temp);
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [x,y]=wfmread(filename,mode)
% wfmread   :   Reads binary *.WFM file from a Tektronix TDS digitizer
% Usage:
% [x,y]=wfmread(filename,mode);
% x and y are one dimensional arrays containing the horizontal and vertical
% data from the signal contained in filename.
%
% filename is a character string name of a *.WFM file.  If filename is
% omitted or empty, a graphical window is launched to help the user choose
% the desired file.
%
% mode allows the user to get some extra information about the digitizer
% settings.  This feature is disabled by default, but can be activated by
% setting mode to 'verbose'.
%
%
% Notes:
% I wrote this a number of years ago as a graduate student so I could store
% many digitizer signals on a 3.5" floppy disk.  I cannot guarantee that
% this will work on more modern digitizers as Tektronix might change their
% format.  As with any function, use at your own risk
%
% Created 6/10/2000 by Daniel Dolan
% Last updated 9/14/2004
%

%%%%%%%%%%%%%%%%%%%%%
% argument checking %
%%%%%%%%%%%%%%%%%%%%%
if nargout<2
   error('Two outputs required');
end

switch nargin
case 0
   filename='';
   mode='quiet';
case 1
   mode='quiet';
case 2
   if ~strcmp(mode,'verbose')& ~strcmp(mode,'quiet')
      error('Invalid mode');
   end
otherwise
   error('Too many arguments');
end

if isempty(filename)
    [fname,pname]=uigetfile('*.wfm;*.WFM','Choose Tektronix WFM file');
    filename=[pname fname];
end

if exist(filename)~=2
   error('Invalid file name');
end

%%%%%%%%%%%%%%%%
% main program %
%%%%%%%%%%%%%%%%
[fp,message]=fopen(filename,'rb','b'); % big-endian format
if fp==-1
    error(['Unable to open file: ' filename]);
end

% main header
a=fread(fp,1,'uchar');
if char(a)==':'
   start=fread(fp,7,'uchar');
else
   start=fread(fp,6,'uchar');
end
start=char([a start']);
count_bytes=str2num(char(fread(fp,1,'uchar'))'); % number of bytes in count
bytes=str2num(char(fread(fp,count_bytes,'uchar'))'); % bytes to end of file
magic_num=fread(fp,1,'int32'); % "magic number"
length=fread(fp,1,'int32'); % length of header+curve data

% reference header
vertPos=fread(fp,1,'double');
horzPos=fread(fp,1,'double');
vertZoom=fread(fp,1,'double');
horzZoom=fread(fp,1,'double');

% waveform header
acqmode=fread(fp,1,'int16');
minMaxFormat=fread(fp,1,'int16');
duration=fread(fp,1,'double');
vertCpl=fread(fp,1,'int16');
horzUnit=fread(fp,1,'int16');
horzScalePerPoint=fread(fp,1,'double');
vertUnit=fread(fp,1,'int16');
vertOffset=fread(fp,1,'double');
vertPos=fread(fp,1,'double');
vertGain=fread(fp,1,'double');
%recordLength=fread(fp,1,'uint64');
recordLength=fread(fp,1,'uint32');
trigPos=fread(fp,1,'int16');
wfmHeaderVersion=fread(fp,1,'int16');
sampleDensity=fread(fp,1,'int16');
burstSegmentLength=fread(fp,1,'int16');
sourceWfm=fread(fp,1,'int16');
video1=fread(fp,3,'int16');
video2=fread(fp,1,'double');
video3=fread(fp,1,'int16');

% check to see extended header is present
if (length-2*recordLength-2*32)==196
   % stuff
end

% print some important info for the user
if strcmp(mode,'verbose')
   switch acqmode
   case 285
      fprintf('Sample Mode (number of data points= %d)\n',recordLength);
   end

   switch horzUnit
   case 609
      fprintf('Horizontal units are volts \n');
   case 610
      fprintf('Horizontal units are seconds \n');
      fprintf('\t Sampling rate= %g Hz \n',horzScalePerPoint)
      fprintf('\t Duration of sweep= %g s\n',duration);
   end

   switch vertUnit
   case 609
      fprintf('Vertical units are volts \n');
   case 610
      fprintf('Vertical units are seconds \n');
   end

   switch vertCpl
   case 565
      fprintf('DC coupling \n');
   case 566
      fprintf('AC coupling \n');
   end

end

% extract the data
preamble=fread(fp,16,'int16');
data=fread(fp,recordLength,'int16');
postamble=fread(fp,16,'int16');
checksum=fread(fp,1,'int16');
fclose(fp);

% calculate the x and y data series
ii=0:(recordLength-1);
x=(ii'-recordLength*trigPos/100)*horzScalePerPoint;
y=data*vertGain/25/256+vertOffset-vertPos*vertGain;

author: Sable

repo: mcbench-benchmarks


text:
function [channels,overall_amplitude] = eeg_amplitudearea(EEG, channels, resrate, wstart, wend)
% EEG_AMPLITUDEAREA - Resamples an ERP average using spline interpolation
%                       at a new sample rate (resrate) in Hz to get the exact limits
%                       of the window of integration. Finely samples the window
%                       and adds together very narrow rectangles capped by
%                       right-angled triangles under the window. Output is in uV.
%                       Trade-off between speed and number of resamples and number of
%                        channels selected occurs.
% Usage:
%      >> [channels, amplitude] = eeg_amplitudearea(EEG,channels, resrate, wstart, wend);
% Inputs:
%   EEG         - EEGLAB data struct containing a (3-D) epoched data matrix
%   channels    - vector of channel indices
%   resrate     - resampling rate for window of integration in Hz
%   wstart      - start of window of integration in ms post stimulus-onset
%   wend        - end of window of integration in ms post stimulus-onset
%
% Outputs:
%   channels    - a vector of channel indices.
%   amplitude   - 1-dimensional array in uV for the channels
%
% Example
%    >> [channels, amplitude] = eeg_amplitudearea(EEG,[12 18 25 29], 2000, 90.52, 120.52);
%
% Author: Tom Campbell, Helsinki Collegium for Advanced Studies, Biomag Laboratory,
%         Engineering Centre, Helsinki University Central Hospital Helsinki Brain
%         Research Centre (tom.campbell@helsinki.fi) Spartam nanctus es: Hanc exorna.
%         Combined with AMPLITUDEAREA_MSUV by Darren Weber, UCSF 28/1/05
%         Retested and debugged Tom Campbell 2/2/05
%         Reconceived, factored somewhat, tested and debugged Tom Campbell 13:24 23.3.2005

if wstart > wend
    error ('ERROR: wstart must be greater than wend')
else
    [channels, overall_amplitude] = eeg_amplitudearea_msuV (EEG,channels, resrate, wstart, wend);
    overall_amplitude = overall_amplitude/(wend - wstart);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [channels, overall_area] = eeg_amplitudearea_msuV (EEG, channels, resrate, wstart, wend)

%if ndim(EEG.data) ~= 3
%  error('EEG.data must be 3-D data epochs');
%end
erp = mean(EEG.data,3);

[tmp ind1] =min( abs( EEG.times - wstart ) ); % closest time index to wstart
[tmp ind2] =min( abs( EEG.times - wend ) ); % closest time index to wend
restep = 1/resrate;

if EEG.times(ind1) > wstart
    ind1= ind1 -1;
end

if EEG.times(ind2) < wend
    ind2= ind2 +1;
end

for x= ind1:ind2
    t = (x -ind1)+1;
    tim(t) = EEG.times(x);
end

tr = 1;
timr(tr) = wstart;
while timr(tr) < wend
    tr = tr + 1;
    timr(tr) = timr(tr-1)+ restep;
end

for x = 1:size(channels,2)
    channel = channels(x);
    %resamples
    rerp(x, 1:tr) = spline(tim(:),erp(channel, ind1:ind2), timr(1:tr));
    pent = timr(tr - 1);
    overall_area(x) = 0;
    for y = 1:(tr -1)
        v1 =  rerp(x,(y));
        v2 =  rerp(x,(y+1));
        if ((v1 > 0) && (v2 < 0)) || ((v1 < 0) && (v2 > 0))
            if (y == (tr-1)) && (timr(y+1)> wend)
                area1 = zero_crossing_truncated(v1, v2, restep, wend, pent);
            else
                area1 = zero_crossing(v1, v2, restep);
            end
        else
            if( y == (tr-1)) && (timr(y+1)> wend)
                area1 = rect_tri_truncated(v1, v2, restep,wend,pent);
            else
                area1 = rect_tri(v1, v2, restep);
            end
        end
        overall_area(x) = overall_area(x) + area1;
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = zero_crossing(v1,v2,step)
    if (v1 > v2)
        T1 = v1;
        T2 = v2;
    else
        T1 = v2;
        T2 = v1;
    end
    tantheta = (abs(T1)+ abs(T2))/step;
    if (v1 > v2)
        %decline
        z = abs(T1)/tantheta;
        tr1= abs(T1)*(z/2);
        tr2= abs(T2)*((step-z)/2);
    else
        %incline
        z = abs(T2)/tantheta;
        tr2= abs(T2)*(z/2);
        tr1= abs(T1)*((step-z)/2);
    end
    [area] = (tr1 - tr2);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = zero_crossing_truncated(v1,v2,step,wend,pent)
    if (v1 > v2)
        T1 = v1;
        T2 = v2;
    else
        T1 = v2;
        T2 = v1;
    end
    tantheta = (abs(T1)+ abs(T2))/step;
    s = wend - pent;
    if (v1 > v2)
        z = abs(T1)/tantheta
        if s < z
            %decline,truncated before zerocrossing
            t1 = tantheta * s;
            r1 = abs(T1)-abs(t1);
            tr1= abs(t1)*(s/2);
            tr2= 0;
            rect1 = r1*s;
            rect2 = 0;
        else
            %decline,truncated after zerocrossing
            t2= tantheta*(s-z);
            tr1= abs(T1)*(z/2);
            tr2 = abs(t2)*((s-z)/2);
            rect1 = 0;
            rect2 = 0;
        end
    else
        z = abs(T2)/tantheta;
        if s < z
            %incline,truncated before zerocrossing
            t2 = tantheta * s;
            r2 = abs(T2)-abs(t2);
            tr1= 0;
            tr2= abs(t2)*(s/2);
            rect1 = 0;
            rect2 = r2*s;
        else
            %incline,truncated after zerocrossing
            t1= tantheta*(s-z);
            tr1 = abs(t1)*((s-z)/2);
            tr2 = abs(T2) * (z/2);
            rect1 = 0;
            rect2 = 0;
        end
    end

[area] = ((rect1 + tr1) - (rect2 + tr2));
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = rect_tri(v1,v2,step)
    if (abs(v1) > abs(v2))
        T = abs(v1)-abs(v2);
        R = abs(v2);
    else
        T = abs(v2)-abs(v1);
        R = abs(v1);
    end
    rect = R*step;
    tri = T*(step/2);
    if v1 > 0
        area = 1* (rect+tri);
    else
        area = -1 * (rect+tri);
    end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [area] = rect_tri_truncated(v1,v2,step,wend,pent)
    if (abs(v1) > abs(v2))
        T = abs(v1)-abs(v2);
        R = abs(v2);
    else
        T = abs(v2)-abs(v1);
        R = abs(v1);
    end

    tantheta = abs(T)/step;
    s = wend -pent;

    if (v1>0)
        if v1 >v2
        %positive decline
            t = tantheta*s;
            e = abs(T)-abs(t);
            rect = s*R;
            exrect = s*e;
            tri = (s/2)*R;
        else
        %positive incline
            t = tantheta*s;
            rect = s*R;
            exrect = 0;
            tri = (s/2)*R;
        end
    else
       if v1 >v2
        %negative decline
            t = tantheta*s;
            rect = s*R;
            exrect = 0;
            tri = (s/2)*R;
        else
        %negative incline
            t = tantheta*s;
            e = abs(T)-abs(t);
            rect = s*R;
            exrect = s*e;
            tri = (s/2)*R;
        end
    end
    tri = T*(step/2);
    if v1 > 0
        area = 1* (rect+exrect+tri);
    else
        area = -1 * (rect+exrect+tri);
    end
return

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function submit(partId)
%SUBMIT Submit your code and output to the ml-class servers
%   SUBMIT() will connect to the ml-class server and submit your solution

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

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

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

  fprintf('\n== Connecting to ml-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    % Get Challenge
    [login, ch, signature] = getChallenge(login);
    if isempty(login) || isempty(ch) || isempty(signature)
      % Some error occured, error string in first return element.
      fprintf('\n!! Error: %s\n\n', login);
      return
    end

    % Attempt Submission with Challenge
    ch_resp = challengeResponse(login, password, ch);
    [result, str] = submitSolution(login, ch_resp, partId, output(partId), ...
                                 source(partId), signature);

    fprintf('\n== [ml-class] Submitted Homework %s - Part %d - %s\n', ...
            homework_id(), partId, partNames{partId});
    fprintf('== %s\n', strtrim(str));
    if exist('OCTAVE_VERSION')
      fflush(stdout);
    end
  end

end

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

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

function [partNames] = validParts()
  partNames = { 'Regularized Linear Regression Cost Function', ...
                'Regularized Linear Regression Gradient', ...
                'Learning Curve', ...
                'Polynomial Feature Mapping' ...
                'Validation Curve' ...
                };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'linearRegCostFunction.m' }, ...
           { 'linearRegCostFunction.m' }, ...
           { 'learningCurve.m' }, ...
           { 'polyFeatures.m' }, ...
           { 'validationCurve.m' } };
end

function out = output(partId)
  % Random Test Cases
  X = [ones(10,1) sin(1:1.5:15)' cos(1:1.5:15)'];
  y = sin(1:3:30)';
  Xval = [ones(10,1) sin(0:1.5:14)' cos(0:1.5:14)'];
  yval = sin(1:10)';
  if partId == 1
    [J] = linearRegCostFunction(X, y, [0.1 0.2 0.3]', 0.5);
    out = sprintf('%0.5f ', J);
  elseif partId == 2
    [J, grad] = linearRegCostFunction(X, y, [0.1 0.2 0.3]', 0.5);
    out = sprintf('%0.5f ', grad);
  elseif partId == 3
    [error_train, error_val] = ...
        learningCurve(X, y, Xval, yval, 1);
    out = sprintf('%0.5f ', [error_train(:); error_val(:)]);
  elseif partId == 4
    [X_poly] = polyFeatures(X(2,:)', 8);
    out = sprintf('%0.5f ', X_poly);
  elseif partId == 5
    [lambda_vec, error_train, error_val] = ...
        validationCurve(X, y, Xval, yval);
    out = sprintf('%0.5f ', ...
        [lambda_vec(:); error_train(:); error_val(:)]);
  end
end

function url = challenge_url()
  url = 'http://www.ml-class.org/course/homework/challenge';
end

function url = submit_url()
  url = 'http://www.ml-class.org/course/homework/submit';
end

% ========================= CHALLENGE HELPERS =========================

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

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

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

function [email,ch,signature] = getChallenge(email)
  str = urlread(challenge_url(), 'post', {'email_address', email});

  str = strtrim(str);
  [email, str] = strtok (str, '|');
  [ch, str] = strtok (str, '|');
  [signature, str] = strtok (str, '|');
end

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

  params = {'homework', homework_id(), ...
            'part', num2str(part), ...
            'email', email, ...
            'output', output, ...
            'source', source, ...
            'challenge_response', ch_resp, ...
            'signature', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

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

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

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

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

function [str] = challengeResponse(email, passwd, challenge)
  salt = ')~/|]QMB3[!W`?OVt7qC"@+}';
  str = sha1([challenge sha1([salt email passwd])]);
  sel = randperm(numel(str));
  sel = sort(sel(1:16));
  str = str(sel);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function  [mls,row,col] = mls(n)

%[mls,row,col] = mls(n);
%
%Generates a Maximum Length Sequence of n bits by utilizing a
%linear feedback shift register with an XOR gate on the tap bits
%Also given are the permutation vector row and col, used by prior and after
%the autocorreltion calculation via FHT.
%
%Function can accept bit lengths of between 2 and 24
%
%y is a vector of 1's & -1's that is (2^n)-1 in length.
%
%reference:
%	Davies, W.D.T. (June, July, August, 1966). Generation and
%properties of maximum-length sequences. Control, 302-4, 364-5,431-3.
%
%Spring 2001, Christopher Brown, cbrown@phi.luc.edu

switch n							%assign taps which will yeild a maximum
case 2								%length sequence for a given bit length
   taps=2;							%I forget the reference I used, but theres
   tap1=1;							%a list of appropriate tap values in
   tap2=2;							%Vanderkooy, JAES, 42(4), 1994.
case 3
   taps=2;
   tap1=1;
   tap2=3;
case 4
   taps=2;
   tap1=1;
   tap2=4;
case 5
   taps=2;
   tap1=2;
   tap2=5;
case 6
   taps=2;
   tap1=1;
   tap2=6;
case 7
   taps=2;
   tap1=1;
   tap2=7;
case 8
   taps=4;
   tap1=1;
   tap2=5;
   tap3=6;
   tap4=8;
case 9
   taps=2;
   tap1=4;
   tap2=9;
case 10
   taps=2;
   tap1=3;
   tap2=10;
case 11
   taps=2;
   tap1=2;
   tap2=11;
case 12
   taps=4;
   tap1=3;
   tap2=4;
   tap3=7;
   tap4=12;
case 13
   taps=4;
   tap1=1;
   tap2=3;
   tap3=4;
   tap4=13;
case 14
   taps=4;
   tap1=1;
   tap2=11;
   tap3=12;
   tap4=14;
case 15
   taps=2;
   tap1=1;
   tap2=15;
case 16
   taps=4;
   tap1=2;
   tap2=3;
   tap3=5;
   tap4=16;
case 17
   taps=2;
   tap1=3;
   tap2=17;
case 18
   taps=2;
   tap1=7;
   tap2=18;
case 19
   taps=4;
   tap1=1;
   tap2=5;
   tap3=6;
   tap4=19;
case 20
   taps=2;
   tap1=3;
   tap2=20;
case 21
   taps=2;
   tap1=2;
   tap2=21;
case 22
   taps=2;
   tap1=1;
   tap2=22;
case 23
   taps=2;
   tap1=5;
   tap2=23;
case 24
   taps=4;
   tap1=1;
   tap2=3;
   tap3=4;
   tap4=24;
case 25
  taps=2;
  tap1=3;
  tap2=25;
case 26
  taps=4;
  tap1=1;
  tap2=7;
  tap3=8;
  tap4=26;
case 27
  taps=4;
  tap1=1;
  tap2=7;
  tap3=8;
  tap4=27;
case 28
  taps=2;
  tap1=3;
  tap2=28;
case 29
  taps=2;
  tap1=2;
  tap2=29;
case 30
  taps=4;
  tap1=1;
  tap2=15;
  tap3=16;
  tap4=30;
case 31
  taps=2;
  tap1=3;
  tap2=31;
case 32
  taps=4;
  tap1=1;
  tap2=27;
  tap3=28;
  tap4=32;
otherwise
   disp(' ');
   error('input bits must be between 2 and 32');
end

if taps == 2
    [mls,row,col]=mls2tap(tap2,tap1,tap2);
elseif taps == 4
    [mls,row,col]=mls4tap(tap4,tap1,tap2,tap3,tap4);
end


text:
function y = vl_sigmoid(x)
% VL_SIGMOID Sigmoid function
%   Y = VL_SIGMOID(X) returns
%
%    Y = 1 ./ (1 + EXP(X)) ;
%
%   Remark::
%     Useful properties of the sigmoid function are:
%
%     -  1 - VL_SIGMOID(X) = VL_SIGMOID(-X)
%     -  Centered sigmoid: 2 * VL_SIGMOID(X) - 1 ;
%     -  VL_SIGMOID(X) = (EXP(X/2) - EXP(X/2)) / (EXP(X/2) + EXP(X/2))
%
%   See also: VL_DSIGMOID(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% 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).

y = 1 ./ (1 + exp(-x)) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function scale = p16_scale ( neqn )

%*****************************************************************************80
%
%% P16_SCALE returns scale factors for problem p16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function determ = fibonacci3_determinant ( n )

%*****************************************************************************80
%
%% FIBONACCI3_DETERMINANT returns the determinant of the FIBONACCI3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  f1 = 0;
  f2 = 0;
  f3 = 1;

  for i = 1 : n
    f1 = f2;
    f2 = f3;
    f3 = f1 + f2;
  end

  determ = f3;

  return
end


text:
function varargout = spm_mip_ui(varargin)
% GUI for displaying MIPs with interactive pointers
% FORMAT hMIPax = spm_mip_ui(Z,XYZ,M,DIM,F,units)
% Z       - {1 x ?} vector point list of SPM values for MIP
% XYZ     - {3 x ?} matrix of coordinates of points (Talairach coordinates)
% M       - voxels - > mm matrix
% DIM     - image dimensions {voxels}
% F       - Figure (or axes) to work in [Defaults to gcf]
% hMIPax  - handle of MIP axes
% units   - units of space
%
% FORMAT xyz = spm_mip_ui('GetCoords',h)
% h       - Handle of MIP axes, or figure containing MIP axis [default gcf]
% xyz     - Current Talairach coordinates of cursor
%
% FORMAT [xyz,d] = spm_mip_ui('SetCoords',xyz,h,hC)
% xyz     - (Input) {3 x 1} vector of desired Talairach coordinates
% h       - Handle of MIP axes, or figure containing MIP axis [default gcf]
% hC      - Handle of calling object, if used as a callback. [Default 0]
% xyz     - (Output) {3 x 1} vector of voxel centre nearest desired xyz co-ords
% d       - Euclidian distance from desired co-ords & nearest voxel
%__________________________________________________________________________
%
% spm_mip_ui displays a maximum intensity projection (using spm_mip)
% with draggable cursors.
%
% See spm_mip.m for details of MIP construction, display, and the brain
% outlines used.
%                           ----------------
%
% The cursor can be dragged to new locations in three ways:
%
% (1) Point & drop: Using the primary "select" mouse button, click on a
%     cursor and drag the crosshair which appears to the desired location.
%     On dropping, the cursors jump to the voxel centre nearest the drop
%     site.
%
% (2) Dynamic drag & drop: Using the middle "extend" mouse button, click on
%     a cursor and drag it about. The cursor follows the mouse, jumping to
%     the voxel centre nearest the pointer. A dynamically updating
%     information line appears above the MIP and gives the current
%     co-ordinates. If the current voxel centre is in the XYZ pointlist,
%     then the corresponding image value is also printed.
%
% (3) Magnetic drag & drop: As with "Dynamic drag & drop", except the cursors
%     jump to the voxel centre in the pointlist nearest to the cursor. Use
%     the right "alt" mouse button for "magnetic drag & drop".
%
% In addition a ContextMenu is provided, giving the option to jump the
% cursors to the nearest suprathreshold voxel, the nearest local
% maximum, or to the global maximum. (Right click on the MIP to bring up
% the ContextMenu.) A message in the MATLAB command window describes the
% jump.
%
%                           ----------------
%
% The current cursor position (constrained to lie on a voxel) can be
% obtained by xyz=spm_mip_ui('GetCoords',hMIPax), and set with
% xyz=spm_mip_ui('SetCoords',xyz,hMIPax), where hMIPax is the handle of
% the MIP axes, or of the figure containing a single MIP [default gcf].
% The latter rounds xyz to the nearest voxel center, returning the
% result.
%
% spm_mip_ui handles all the callbacks required for moving the cursors, and
% is "registry" enabled (See spm_XYZreg.m). Programmers help is below in the
% main body of the function.
%
%__________________________________________________________________________
% Copyright (C) 1996-2018 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_mip_ui.m 7388 2018-08-06 12:04:26Z guillaume $

%==========================================================================
% - FORMAT specifications for embedded CallBack functions
%==========================================================================
%( This is a multi function function, the first argument is an action  )
%( string, specifying the particular action function to take.          )
%
% FORMAT hMIPax = spm_mip_ui(Z,XYZ,V,F,units)
% [ShortCut] Defaults to hMIPax=spm_mip_ui('Display',Z,XYZ,V,F,units)
%
% FORMAT hMIPax = spm_mip_ui('Display',Z,XYZ,M,DIM,F,units)
% Displays the MIP and sets up cursors
% Z       - {1 x ?} vector point list of SPM values for MIP
% XYZ     - {3 x ?} matrix of coordinates of points (Talairach coordinates)
% M       - voxels - > mm matrix
% DIM     - image dimensions {voxels}
% F       - Handle of figure (or axes) to work in [Defaults to gcf]
% hMIPax  - handle of MIP axes
% units   - units of space
%
% FORMAT xyz = spm_mip_ui('GetCoords',h)
% Returns coordinates of current cursor position
% h       - Handle of MIP axes [defaults to spm_mip_ui('FindMIPax')]
% xyz     - Current Talairach coordinates of cursor
%
% FORMAT [xyz,d] = spm_mip_ui('SetCoords',xyz,h,hC)
% Sets cursor position
% xyz     - (Input) {3 x 1} vector of desired Talairach coordinates
% h       - Handle of MIP axes [defaults to spm_mip_ui('FindMIPax')]
% hC      - Handle of calling object, if used as a callback. [Default 0]
% xyz     - (Output) {3 x 1} vector of voxel centre nearest desired xyz co-ords
% d       - Euclidian distance from desired co-ords & nearest voxel
%
% FORMAT spm_mip_ui('PosnMarkerPoints',xyz,h)
% Utility routine: Positions cursor markers
% xyz     - {3 x 1} vector of Talairach coordinates for cursor
% h       - Handle of MIP axes [defaults to spm_mip_ui('FindMIPax')]
%
% FORMAT [xyz,d] = spm_mip_ui('Jump',h,loc)
% Utility routine (CallBack of UIcontextMenu) to jump cursor
% h       - Handle of MIP axes [defaults to spm_mip_ui('FindMIPax')]
% loc     - String defining jump: 'dntmv' - don't move
%                                 'nrvox' - nearest suprathreshold voxel
%                                 'nrmax' - nearest local maximum
%                                 'glmax' - global maximum
% xyz     - co-ordinates of voxel centre jumped to {3 x 1} vector
% d       - (square) Euclidian distance jumped
%
% FORMAT hMIPax = spm_mip_ui('FindMIPax',h)
% Looks for / checks MIP axes 'Tag'ged 'hMIPax'... errors if no valid MIP axes
% h       - Handle of MIP axes, or figure containing MIP axis [default gcf]
% hMIPax  - Handle of valid MIP axis found (errors if multiple found)
%
% FORMAT spm_mip_ui('MoveStart')
% Utility routine: CallBack for starting cursor dragging
% This is the ButtonDownFcn for the cursor markers
%
% FORMAT spm_mip_ui('Move',DragType)
% Utility routine: Initiate cursor move
% DragType - 0 = Point & drop (no dragging)
%            1 = Drag'n'drop with dynamic coordinate/value updating
%            2 = Magnetic drag'n'drop with dynamic coordinate updating
%
% FORMAT spm_mip_ui('MoveEnd')
% Utility routine: End cursor move
%__________________________________________________________________________

%-Condition arguments
%==========================================================================
if nargin==0
    error('Insufficient arguments')
elseif ~ischar(varargin{1})
    varargout={spm_mip_ui('Display',varargin{1:end})}; return
end

%-Axis offsets for 3d MIPs:
%==========================================================================
%-MIP pane dimensions, origin offsets and #pixels per mm
%-See spm_project.c for derivation
mipmat = char(spm_get_defaults('stats.results.mipmat'));
load(mipmat, 'DXYZ', 'CXYZ', 'scale');
% DMIP = [DXYZ(2)+DXYZ(1), DXYZ(1)+DXYZ(3)];
%-Coordinates of Talairach origin in multipane MIP image (Axes are 'ij' + rot90)
% Transverse: [Po(1), Po(2)]
% Sagittal  : [Po(1), Po(3)]
% Coronal   : [Po(4), Po(3)]
% 4 voxel offsets in Y since using character '<' as a pointer.
Po(1)  =                  CXYZ(2) -2;
Po(2)  = DXYZ(3)+DXYZ(1) -CXYZ(1) +2;
Po(3)  = DXYZ(3)         -CXYZ(3) +2;
Po(4)  = DXYZ(2)         +CXYZ(1) -2;

%==========================================================================
switch lower(varargin{1}), case 'display'
%==========================================================================
    % hMIPax = spm_mip_ui('Display',Z,XYZ,M,DIM,F,units)
    if nargin<5
        F      = gcf;
        hMIPax = [];
    else
        F = varargin{6};
        if ischar(F), F=spm_figure('FindWin',F); end
        if ~ishandle(F), error('Invalid handle'), end
        switch get(F,'Type'), case 'figure'
            hMIPax = [];
            case 'axes'
                hMIPax = F;
                F      = get(hMIPax,'Parent');
            otherwise
                error('F not a figure or axis handle')
        end
    end
    if nargin<4, error('Insufficient arguments'), end
    Z       = varargin{2};
    XYZ     = varargin{3};
    M       = varargin{4};
    DIM     = varargin{5};
    try
        units = varargin{7};
    catch
        units = {'mm' 'mm' 'mm'};
    end

    xyz = spm_XYZreg('RoundCoords',[0;0;0],M,DIM);

    %-Display (MIP) transformation matrix
    %----------------------------------------------------------------------
    % Md maps various non-spatial data
    % Ms maps spatial (mm) to pixel coordinates, it must not change
    % mapping of non-spatial coordinates
    % Ms is left-multiplied to Md
    if isequal(units,{'mm' 'mm' 'mm'})
        Md      = eye(4);
        Ms      = diag([scale(1:3) 1]);
    elseif isequal(units,{'mm' 'mm' ''})
        Md      = eye(4);
        Md(3,4) = -100;  % out of field of view (MNI)
        Ms      = diag([scale(1:2) 1 1]);
    elseif isequal(units,{'mm' 'mm' 'ms'}) || isequal(units,{'mm' 'mm' 'Hz'})
        Md      = eye(4);
        Md(3,3) = 80 / (M(3,3)*DIM(3));
        Md(3,4) = -80 * M(3,4) / (M(3,3)*DIM(3));
        Ms      = diag([scale(1:2) 1 1]);
    elseif isequal(units,{'Hz' 'ms' ''})
        Md      = eye(4);
        Md(1,1) = -136 / (M(1,1)*DIM(1));
        Md(1,4) = M(1,4)*136 / (M(1,1)*DIM(1)) + 68;
        Md(2,2) = 172 / (M(2,2)*DIM(2));
        Md(2,4) = -M(2,4)*172 / (M(2,2)*DIM(2)) - 100;
        Md(3,4) = -100;
        Ms      = eye(4);
    elseif isequal(units,{'Hz' 'Hz' ''})
        Md      = eye(4);
        Md(1,1) = 136 / (M(1,1)*DIM(1));
        Md(1,4) = -M(1,4)*136 / (M(1,1)*DIM(1)) - 68;
        Md(2,2) = 172 / (M(2,2)*DIM(2));
        Md(2,4) = -M(2,4)*172 / (M(2,2)*DIM(2)) - 100;
        Md(3,4) = -100;
        Ms      = eye(4);
    elseif isequal(units,{'mm' 'mm' '%'})
        warning('Handling of data units changed: please re-estimate model.');
        units   = {'mm' 'mm' 'ms'};
        Md      = eye(4);
        Md(3,3) = 80 / (M(3,3)*DIM(3));
        Md(3,4) = -80 * M(3,4) / (M(3,3)*DIM(3));
        Ms      = diag([scale(1:2) 1 1]);
    else
        error('Unknown data units.');
    end

    %-Display MIP
    %----------------------------------------------------------------------
    Funits = get(F,'Units');
    set(F,'Units','normalized')
    if isempty(hMIPax)
        hMIPax = axes('Position',[0.24 0.54 0.62 0.42],'Parent',F);
    else
        axes(hMIPax), cla reset
    end

    %-NB: spm_mip's `image` uses a newplot, & screws stuff without the figure.
    figure(F)
    pXYZ = Ms*Md*[XYZ;ones(1,size(XYZ,2))];
    spm_mip(Z,pXYZ(1:3,:),Ms*Md*M,units);
    hMIPim = get(gca,'Children');

    %-Print coordinates
    %----------------------------------------------------------------------
    hMIPxyz = text(0,max(get(hMIPax,'YLim'))/2,...
        {'{\bfSPM}{\itmip}',sprintf('[%g, %g, %g]',xyz(1:3))},...
        'Interpreter','TeX','FontName',spm_platform('font','times'),...
        'Color',[1,1,1]*.7,...
        'HorizontalAlignment','Center',...
        'VerticalAlignment','Bottom',...
        'Rotation',90,...
        'Tag','hMIPxyz',...
        'UserData',xyz);

    %-Create point markers
    %----------------------------------------------------------------------
    xyz = Ms*Md*[xyz(:);1];
    xyz = xyz(1:3);
    hX1r  = text(Po(1)+xyz(2),Po(2)+xyz(1),'<',...
        'Color','r','Fontsize',16,...
        'Tag','hX1r',...
        'ButtonDownFcn','spm_mip_ui(''MoveStart'')');
    hX2r  = text(Po(1)+xyz(2),Po(3)-xyz(3),'<',...
        'Color','r','Fontsize',16,...
        'Tag','hX2r',...
        'ButtonDownFcn','spm_mip_ui(''MoveStart'')');
    hX3r  = text(Po(4)+xyz(1),Po(3)-xyz(3),'<',...
        'Color','r','Fontsize',16,...
        'Tag','hX3r',...
        'ButtonDownFcn','spm_mip_ui(''MoveStart'')');
    hXr   = [hX1r,hX2r,hX3r];

    if DIM(3) == 1
        %-2 dimensional data
        %------------------------------------------------------------------
        set(hXr(3),'Visible','off');
        set(hXr(2),'Visible','off');

    end

    %-Create UIContextMenu for marker jumping
    %-----------------------------------------------------------------------
    h = uicontextmenu('Tag','MIPconmen','UserData',hMIPax,'Parent',F);
    if isempty(XYZ), str='off'; else str='on'; end
    uimenu(h,'Separator','off','Label','goto nearest suprathreshold voxel',...
        'CallBack',['spm_mip_ui(''Jump'',',...
        'get(get(gcbo,''Parent''),''UserData''),''nrvox'');'],...
        'Interruptible','off','BusyAction','Cancel','Enable',str);
    uimenu(h,'Separator','off','Label','goto nearest local maximum',...
        'CallBack',['spm_mip_ui(''Jump'',',...
        'get(get(gcbo,''Parent''),''UserData''),''nrmax'');'],...
        'Interruptible','off','BusyAction','Cancel','Enable',str);
    uimenu(h,'Separator','off','Label','goto global maximum',...
        'CallBack',['spm_mip_ui(''Jump'',',...
        'get(get(gcbo,''Parent''),''UserData''),''glmax'');'],...
        'Interruptible','off','BusyAction','Cancel','Enable',str);
    uimenu(h,'Separator','on','Label','save MIP as...',...
        'CallBack',['spm_mip_ui(''Save'', ',...
        'get(get(gcbo,''Parent''),''UserData''));'],...
        'Interruptible','off','BusyAction','Cancel');
    h1 = uimenu(h,'Separator','on','Label','Extract data');
    h2 = uimenu(h1,'Label','raw y');
    uimenu(h2,'Label','This voxel',...
        'CallBack','y=spm_mip_ui(''Extract'', ''Y'', ''voxel'')',...
        'Interruptible','off','BusyAction','Cancel');
    uimenu(h2,'Label','This cluster',...
        'CallBack','y=spm_mip_ui(''Extract'', ''Y'', ''cluster'')',...
        'Interruptible','off','BusyAction','Cancel');
    h2 = uimenu(h1,'Label','whitened and filtered y');
    uimenu(h2,'Label','This voxel',...
        'CallBack','y=spm_mip_ui(''Extract'', ''KWY'', ''voxel'')',...
        'Interruptible','off','BusyAction','Cancel');
    uimenu(h2,'Label','This cluster',...
        'CallBack','y=spm_mip_ui(''Extract'', ''KWY'', ''cluster'')',...
        'Interruptible','off','BusyAction','Cancel');
    h2 = uimenu(h1,'Label','betas');
    uimenu(h2,'Label','This voxel',...
        'CallBack','beta=spm_mip_ui(''Extract'', ''beta'', ''voxel'')',...
        'Interruptible','off','BusyAction','Cancel');
    uimenu(h2,'Label','This cluster',...
        'CallBack','beta=spm_mip_ui(''Extract'', ''beta'', ''cluster'')',...
        'Interruptible','off','BusyAction','Cancel');

    % overlay sensor positions for M/EEG
    %----------------------------------------------------------------------
    if strcmp(spm('CheckModality'), 'EEG') && ...
            (isequal(units,{'mm' 'mm' 'ms'}) || isequal(units,{'mm' 'mm' 'Hz'}))
        uimenu(h,'Separator','on','Label','display/hide channels',...
            'CallBack',['spm_mip_ui(''Channels'', ',...
            'get(get(gcbo,''Parent''),''UserData''));'],...
            'Interruptible','off','BusyAction','Cancel','Enable',str);
        uimenu(h,'Separator','off','Label','goto nearest suprathreshold channel',...
            'CallBack',['spm_mip_ui(''Jump'',',...
        'get(get(gcbo,''Parent''),''UserData''),''nrchan'');'],...
            'Interruptible','off','BusyAction','Cancel','Enable',str);
    end

    set(hMIPim,'UIContextMenu',h)

    %-Save handles and data
    %----------------------------------------------------------------------
    set(hMIPax,'Tag','hMIPax','UserData',...
        struct(...
        'hReg',     [],...
        'hMIPxyz',  hMIPxyz,...
        'XYZ',      XYZ,...
        'Z',        Z,...
        'M',        M,...
        'Md',       Md,...
        'Ms',       Ms,...
        'DIM',      DIM,...
        'units',    {units},...
        'hXr',      hXr))

    varargout = {hMIPax};

    %======================================================================
    case 'getcoords'
    %======================================================================
        % xyz = spm_mip_ui('GetCoords',h)
        if nargin<2, h=spm_mip_ui('FindMIPax'); else h=varargin{2}; end
        varargout = {get(getfield(get(h,'UserData'),'hMIPxyz')  ,'UserData')};

    %======================================================================
    case 'setcoords'
    %======================================================================
        % [xyz,d] = spm_mip_ui('SetCoords',xyz,h,hC)
        if nargin<4, hC=NaN; else hC=varargin{4}; end
        if nargin<3, h=spm_mip_ui('FindMIPax'); else h=varargin{3}; end
        if nargin<2, error('Set coords to what?'), else xyz=varargin{2}; end

        MD  = get(h,'UserData');

        %-Check validity of coords only when called without a caller handle
        %------------------------------------------------------------------
        if ~ishandle(hC)
            [xyz,d] = spm_XYZreg('RoundCoords',xyz,MD.M,MD.DIM);
            if d>0 && nargout<2, warning(sprintf(...
                    '%s: Co-ords rounded to nearest voxel center: Discrepancy %.2f',...
                    mfilename,d));
            end
        else
            d = [];
        end

        %-Move marker points, update internal cache in hMIPxyz
        %------------------------------------------------------------------
        spm_mip_ui('PosnMarkerPoints',MD.Md(1:3,:)*[xyz;1],h);
        set(MD.hMIPxyz,'UserData',reshape(xyz(1:3),3,1))
        set(MD.hMIPxyz,'String',{'{\bfSPM}{\itmip}',sprintf('[%g, %g, %g]',xyz(1:3))})

        %-Tell the registry, if we've not been called by the registry...
        %------------------------------------------------------------------
        if ~isempty(MD.hReg) && ~isequal(MD.hReg,hC), spm_XYZreg('SetCoords',xyz,MD.hReg,h); end

        %-Return arguments
        %------------------------------------------------------------------
        varargout = {xyz,d};

    %======================================================================
    case 'posnmarkerpoints'
    %======================================================================
        % spm_mip_ui('PosnMarkerPoints',xyz,h)
        if nargin<3, h=spm_mip_ui('FindMIPax'); else h=varargin{3}; end
        if nargin<2, xyz = spm_mip_ui('GetCoords',h); else xyz = varargin{2}; end
        MD = get(h,'UserData');

        %-Get handles of marker points from UserData of hMIPax
        %------------------------------------------------------------------
        hX = MD.hXr;

        %-Set marker points
        %------------------------------------------------------------------
        set(hX,'Units','Data')
        if length(hX)==1
            vx  = sqrt(sum(MD.M(1:3,1:3).^2));
            tmp = MD.M\[xyz ; 1];
            tmp = tmp(1:2).*vx(1:2)';
            set(hX,'Position',[tmp(1), tmp(2), 1])
        else
            pxyz = MD.Ms*[xyz(1:3);1];
            set(hX(1),'Position',[ Po(1) + pxyz(2), Po(2) + pxyz(1), 0])
            set(hX(2),'Position',[ Po(1) + pxyz(2), Po(3) - pxyz(3), 0])
            set(hX(3),'Position',[ Po(4) + pxyz(1), Po(3) - pxyz(3), 0])
        end

    %======================================================================
    case 'jump'
    %======================================================================
        % [xyz,d] = spm_mip_ui('Jump',h,loc)
        if nargin<3, loc='nrvox'; else loc=varargin{3}; end
        if nargin<2, h=spm_mip_ui('FindMIPax'); else h=varargin{2}; end

        %-Get current location & MipData
        %------------------------------------------------------------------
        oxyz = spm_mip_ui('GetCoords',h);
        MD   = get(h,'UserData');

        %-Compute location to jump to
        %------------------------------------------------------------------
        if isempty(MD.XYZ), loc='dntmv'; end
        switch lower(loc), case 'dntmv'
            spm('alert!','No suprathreshold voxels to jump to!',mfilename,0);
            varargout = {oxyz, 0};
            return
            case 'nrvox'
                str       = 'nearest suprathreshold voxel';
                [xyz,i,d] = spm_XYZreg('NearestXYZ',oxyz,MD.XYZ);
            case 'nrmax'
                str       = 'nearest local maximum';
                iM        = inv(MD.M);
                XYZvox    = iM(1:3,:)*[MD.XYZ; ones(1,size(MD.XYZ,2))];
                [null,null,XYZvox] = spm_max(MD.Z,XYZvox);
                XYZ       = MD.M(1:3,:)*[XYZvox; ones(1,size(XYZvox,2))];
                [xyz,i,d] = spm_XYZreg('NearestXYZ',oxyz,XYZ);
            case 'glmax'
                str       = 'global maximum';
                [null, i] = max(MD.Z); i = i(1);
                xyz       = MD.XYZ(:,i);
                d         = sqrt(sum((oxyz-xyz).^2));
            case 'nrchan'
                str       = 'nearest suprathreshold channel';
                if ~isfield(MD, 'hChanPlot'), spm_mip_ui('Channels',h); end
                MD        = get(h,'UserData');
                if ~isfield(MD, 'hChanPlot'), return; end
                [xyz,i,d] = spm_XYZreg('NearestXYZ',[oxyz(1); oxyz(2); 0],MD.Channels.pos);
                xyz(3)    = oxyz(3);
                str       = [str sprintf(' (%s)',MD.Channels.name{i})];
            otherwise
                warning('Unknown jumpmode')
                varargout = {xyz,0};
                return
        end

        %-Write jump report, jump, and return arguments
        %------------------------------------------------------------------
        fprintf(['\n\t%s:\tJumped %0.2fmm from [%3.0f, %3.0f, %3.0f],\n\t\t\t',...
            'to %s at [%3.0f, %3.0f, %3.0f]\n'],...
            mfilename, d, oxyz, str, xyz)

        spm_mip_ui('SetCoords',xyz,h,h);
        varargout = {xyz, d};

    %======================================================================
    case 'findmipax'
    %======================================================================
        % hMIPax = spm_mip_ui('FindMIPax',h)
        % Checks / finds hMIPax handles
        %-**** h is handle of hMIPax, or figure containing MIP (default gcf)
        if nargin<2, h=get(0,'CurrentFigure'); else h=varargin{2}; end
        if ischar(h), h=spm_figure('FindWin',h); end
        if ~ishandle(h), error('invalid handle'), end
        if ~strcmp(get(h,'Tag'),'hMIPax'), h=findobj(h,'Tag','hMIPax'); end
        if isempty(h), error('MIP axes not found'), end
        if length(h)>1, error('Multiple MIPs in this figure'), end
        varargout = {h};

    %======================================================================
    case 'movestart'
    %======================================================================
        % spm_mip_ui('MoveStart')
        [cO,cF] = gcbo;
        hMIPax  = get(cO,'Parent');
        MD      = get(hMIPax,'UserData');

        %-Store useful quantities in UserData of gcbo, the object to be dragged
        %------------------------------------------------------------------
        set(hMIPax,'Units','Pixels')
        set(cO,'UserData',struct(...
            'hReg',     MD.hReg,...
            'xyz',      spm_mip_ui('GetCoords',hMIPax),...
            'MIPaxPos', get(hMIPax,'Position')*[1,0;0,1;0,0;0,0],...
            'hMIPxyz',  MD.hMIPxyz,...
            'M',        MD.M,...
            'Md',       MD.Md,...
            'Ms',       MD.Ms,...
            'DIM',      MD.DIM,...
            'hX',       MD.hXr))

        %-Initiate dragging
        %------------------------------------------------------------------
        if strcmp(get(cF,'SelectionType'),'normal') || isempty(MD.XYZ)
            %-Set Figure callbacks for drop but no drag (DragType 0)
            %--------------------------------------------------------------
            set(MD.hMIPxyz,'Visible','on','String',...
                {'{\bfSPM}{\itmip}','\itPoint & drop...'})
            set(cF,'WindowButtonUpFcn',    'spm_mip_ui(''Move'',0)',...
                'Interruptible','off')
            set(cF,'Pointer','CrossHair')

       %-Set Figure callbacks for drag'n'drop (DragType 1)
       %-------------------------------------------------------------------
        elseif strcmp(get(cF,'SelectionType'),'extend')
            set(MD.hMIPxyz,'Visible','on','String',...
                {'{\bfSPM}{\itmip}','\itDynamic drag & drop...'})
            set(cF,'WindowButtonMotionFcn','spm_mip_ui(''Move'',1)',...
                'Interruptible','off')
            set(cF,'WindowButtonUpFcn',    'spm_mip_ui(''MoveEnd'')',...
                'Interruptible','off')
            set(cF,'Pointer','Fleur')

        %-Set Figure callbacks for drag'n'drop with co-ord updating (DragType 2)
        %------------------------------------------------------------------
        elseif strcmp(get(cF,'SelectionType'),'alt')
            set(MD.hMIPxyz,'Visible','on','String',...
                {'{\bfSPM}{\itmip}','\itMagnetic drag & drop...'})
            set(cF,'WindowButtonMotionFcn','spm_mip_ui(''Move'',2)',...
                'Interruptible','off')
            set(cF,'WindowButtonUpFcn',    'spm_mip_ui(''MoveEnd'')',...
                'Interruptible','off')
            set(cF,'Pointer','Fleur')
        end

    %======================================================================
    case 'move'
    %======================================================================
        % spm_mip_ui('Move',DragType)
        if nargin<2, DragType = 2; else DragType = varargin{2}; end
        cF = gcbf;
        cO = gco(cF);

        %-Get useful data from UserData of gcbo, the object to be dragged
        %------------------------------------------------------------------
        MS  = get(cO,'UserData');

        %-Work out where we are moving to - Use HandleGraphics to give position
        %------------------------------------------------------------------
        set(cF,'Units','pixels')
        d = get(cF,'CurrentPoint') - MS.MIPaxPos;
        set(cO,'Units','pixels')
        set(cO,'Position',d)
        set(cO,'Units','data')
        d = get(cO,'Position');

        %-Work out xyz, depending on which view is being manipulated
        %------------------------------------------------------------------
        sMarker = get(cO,'Tag');
        pxyz = MS.Ms*[MS.xyz;1];
        if strcmp(sMarker,'hX1r')
            xyz = [d(2) - Po(2); d(1) - Po(1); pxyz(3)];
        elseif strcmp(sMarker,'hX2r')
            xyz = [pxyz(1); d(1) - Po(1); Po(3) - d(2)];
        elseif strcmp(sMarker,'hX3r')
            xyz = [d(1) - Po(4); pxyz(2); Po(3) - d(2)];
        else
            error('Can''t work out which marker point')
        end
        xyz = inv(MS.Ms*MS.Md) * [xyz(:);1]; xyz = xyz(1:3);

        %-Round coordinates according to DragType & set in hMIPxyz's UserData
        %------------------------------------------------------------------
        if DragType==0
            xyz    = spm_XYZreg('RoundCoords',xyz,MS.M,MS.DIM);
        elseif DragType==1
            xyz    = spm_XYZreg('RoundCoords',xyz,MS.M,MS.DIM);
            hMIPax = get(cO,'Parent');
            MD     = get(hMIPax,'UserData');
            i      = spm_XYZreg('FindXYZ',xyz,MD.XYZ);
        elseif DragType==2
            hMIPax = get(cO,'Parent');
            MD     = get(hMIPax,'UserData');
            [xyz,i,d] = spm_XYZreg('NearestXYZ',xyz,MD.XYZ);
        end
        set(MS.hMIPxyz,'UserData',xyz)

        %-Move marker points
        %------------------------------------------------------------------
        set(MS.hX,'Units','Data')
        xyz2 = MS.Ms * MS.Md * [xyz(:);1];
        xyz2 = xyz2(1:3);
        set(MS.hX(1),'Position',[ Po(1) + xyz2(2), Po(2) + xyz2(1), 0])
        set(MS.hX(2),'Position',[ Po(1) + xyz2(2), Po(3) - xyz2(3), 0])
        set(MS.hX(3),'Position',[ Po(4) + xyz2(1), Po(3) - xyz2(3), 0])

        %-Update dynamic co-ordinate strings (if appropriate DragType)
        %------------------------------------------------------------------
        if DragType==0
            spm_mip_ui('MoveEnd')
        elseif DragType==1
            if isempty(i)
                str = {'{\bfSPM}{\itmip}',sprintf('[%g, %g, %g]',xyz)};
            else
                str = {'{\bfSPM}{\itmip}: ',...
                    sprintf('[%g, %g, %g]: %.4f',xyz,MD.Z(i))};
            end
            set(MD.hMIPxyz,'String',str)
        elseif DragType==2
            set(MD.hMIPxyz,'String',...
                {'{\bfSPM}{\itmip}',sprintf('[%g, %g, %g]: %.4f',xyz,MD.Z(i))})
        else
            error('Illegal DragType')
        end

    %======================================================================
    case 'moveend'
    %======================================================================
        % spm_mip_ui('MoveEnd')
        cF = gcbf;
        cO = gco(cF);
        hMIPax  = get(cO,'Parent');
        MS      = get(cO,'UserData');

        %-Reset WindowButton functions, pointer & SPMmip info-string
        %------------------------------------------------------------------
        set(gcbf,'WindowButtonMotionFcn',' ')
        set(gcbf,'WindowButtonUpFcn',' ')
        set(gcbf,'Pointer','arrow')
        set(MS.hMIPxyz,'String',...
            {'{\bfSPM}{\itmip}',sprintf('[%g, %g, %g]',get(MS.hMIPxyz,'UserData'))})

        %-Set coordinates after drag'n'drop, tell registry
        %------------------------------------------------------------------
        % don't need to set internal coordinates 'cos 'move' does that
        if ~isempty(MS.hReg)
            spm_XYZreg('SetCoords',get(MS.hMIPxyz,'UserData'),MS.hReg,hMIPax);
        end

    %======================================================================
    case 'channels'
    %======================================================================
        % spm_mip_ui('Channels',h)
        % M/EEG specific to display channels on 1-slice MIP

        if nargin<2, h=spm_mip_ui('FindMIPax'); else h=varargin{2}; end
        MD  = get(h,'UserData');

        if ~isfield(MD, 'hChanPlot')
            % first time call
            D = spm_eeg_load;
            if isempty(D), return; end

            DIM = get(findobj('Tag','hFxyz'), 'UserData');

            [mod, Cind] = spm_eeg_modality_ui(D, 1, 1);
            otherind = setdiff(1:nchannels(D), Cind);
            if ~isempty(otherind)
                D = chantype(D, otherind, 'Other');
            end
            [Cel x, y] = spm_eeg_locate_channels(D, DIM.DIM(1), Cind);
            Cel = DIM.M * [Cel'; ones(2,size(Cel,1))];
            Cel = Cel(1:2,:)';
            pos = [Cel'; zeros(1,size(Cel,1))];
            Cel(:,1) = Cel(:,1) + Po(2);
            Cel(:,2) = Cel(:,2) + Po(1);
            hold(h,'on'), hChanPlot = plot(h,Cel(:, 2), Cel(:, 1), 'b*');

            hChanText = cell(1,size(Cel,1));
            name = D.chanlabels(Cind);
            figure(spm_figure('FindWin'));
            for i = 1:size(Cel, 1)
                hChanText{i} = text(Cel(i, 2)+0.5, Cel(i, 1), name{i}, 'Color', 'b');
            end

            MD.hChanPlot     = hChanPlot;
            MD.hChanText     = hChanText;
            MD.Channels.pos  = pos;
            MD.Channels.name = D.chanlabels(Cind);
            set(h, 'UserData', MD);
        else
            if strcmp(get(MD.hChanPlot, 'Visible'), 'on');
                % switch it off
                set(MD.hChanPlot, 'Visible', 'off');
                for i = 1:length(MD.hChanText)
                    set(MD.hChanText{i}, 'Visible', 'off');
                end
            else
                % switch it on
                set(MD.hChanPlot, 'Visible', 'on');
                for i = 1:length(MD.hChanText)
                    set(MD.hChanText{i}, 'Visible', 'on');
                end
            end

        end

    %======================================================================
    case 'save'
    %======================================================================
        % spm_mip_ui('Save',h)
        if nargin<2, h=spm_mip_ui('FindMIPax'); else h=varargin{2}; end
        %MD  = get(h,'UserData');
        %pXYZ = MD.Ms*MD.Md*[MD.XYZ;ones(1,size(MD.XYZ,2))];
        %mip  = spm_mip(MD.Z,pXYZ(1:3,:),MD.Ms*MD.Md*MD.M,MD.units);
        mip = get(findobj(h,'Type','image'),'CData');
        [f, p] = uiputfile('*.png', 'Save as', 'mip.png');
        if ~isequal(f,0) && ~isequal(p,0)
            if ndims(mip) == 3
                imwrite(mip,fullfile(p,f),'png');
            else
                imwrite(mip,gray(64),fullfile(p,f),'png');
            end
            fprintf('Saving SPM MIP in %s\n',fullfile(p,f));
        end

    %======================================================================
    case 'extract'
    %======================================================================
        % data = spm_mip_ui('Extract',what,where)
        if nargin<2, what='beta'; else what=varargin{2}; end
        if nargin<3, where='voxel'; else where=varargin{3}; end
        xSPM = evalin('base','xSPM');
        SPM  = evalin('base','SPM');

        XYZmm = spm_results_ui('GetCoords');
        [XYZmm,i] = spm_XYZreg('NearestXYZ',XYZmm,xSPM.XYZmm);
        if isempty(i), varargout = {[]}; return; end
        spm_results_ui('SetCoords',xSPM.XYZmm(:,i));

        switch lower(where)
            case 'voxel'
                % current voxel
                XYZ = SPM.xVol.iM(1:3,:)*[XYZmm;1];
            case 'cluster'
                % current cluster
                A   = spm_clusters(xSPM.XYZ);
                j   = find(A == A(i));
                XYZ = xSPM.XYZ(:,j);
            otherwise
                error('Unknown action.');
        end

        switch lower(what)
            case {'y','kwy'}
                data = spm_get_data(SPM.xY.VY,XYZ);
                if strcmpi(what,'kwy')
                    data = spm_filter(SPM.xX.K,SPM.xX.W*data);
                end
            case 'beta'
                data = spm_get_data(SPM.Vbeta,XYZ);
            otherwise
                error('Unknown action.');
        end

        varargout = {data};

    %======================================================================
    otherwise
    %======================================================================
        error('Unknown action string')

end

author: spm


text:
function [ fr1, fr2 ] = splitDouble( tt )
% breaks up a double number into a double integer part and a fractional part
%
% void split (double tt,
%
%             double *fr)
%
% ------------------------------------------------------------------------
%
%    PURPOSE:
%       This function breaks up a double number into a double integer
%       part and a fractional part.
%
%    REFERENCES:
%       Standish, E.M. and Newhall, X X (1988). "The JPL Export
%          Planetary Ephemeris"; JPL document dated 17 June 1988.
%
%    INPUT
%    ARGUMENTS:
%       tt (double)
%          Input number.
%
%    OUTPUT
%    ARGUMENTS:
%       *fr (double)
%          2-element output array;
%             fr[0] contains integer part,
%             fr[1] contains fractional part.
%          For negative input numbers,
%             fr[0] contains the next more negative integer;
%             fr[1] contains a positive fraction.
%
%    RETURNED
%    VALUE:
%       None.
%
%    GLOBALS
%    USED:
%       None.
%
%    FUNCTIONS
%    CALLED:
%       None.
%
%    VER./DATE/
%    PROGRAMMER:
%       V1.0/06-90/JAB (USNO/NA): CA coding standards
%       V1.1/03-93/WTH (USNO/AA): Convert to C.
%       V1.2/07-93/WTH (USNO/AA): Update to C standards.
%       V1.3/10-10/WKP (USNO/AA): Renamed function to lowercase to
%                                 comply with coding standards.
%
%    NOTES:
%       None.
%
% ------------------------------------------------------------------------

%   Get integer and fractional parts.

   ir = floor(tt);
   fr2 = tt - ir;
   fr1 = ir;

end

author: goGPS-Project


text:
function [COU, FA_STD, ADC_STD] = wild_bootstrapping_DTI( I, bval, Mask, Nreps, DT, VectorField )
% WILD_BOOTSTRAPPING_DTI Returns uncertainty measures from diffusion tensor
% images via wild bootstrapping. See Whitcher 2008 - Using the Wild
% Bootstrap to Quantify Uncertainty in Diffusion Tensor Imaging
%
%
% Inputs:
%
%   I is the input image, of any dimensionality, with the diffusion
%   scans in the last dimension
%
%   BVAL is the b-matrix, of size [3 3 N] (same as in fit_DT)
%
%   MASK is the mask of voxels for analysis, same size as the first N-1
%   dimensions of I.
%
%   NREPS is the number of Monte Carlo repetitions (default 1000)
%
%   DT and VectorField are optional - if you have a ground truth that you
%   want to compare against, instead of the input signal.
%
% Outputs:
%
%   COU is the cone of uncertainty (in degrees) (See Jones 2003 -
%   Determining and Visualizing Uncertainty in Estimates of Fiber
%   Orientation From Diffusion Tensor MRI). Here, the 95th percentile is
%   returned. The COU of the second and third eigenvectors is concatenated
%   in the last dimension.
%
%   FA_STD is the standard deviation of the fractional anisotropy over
%   the repetitions
%
%   ADC_STD is the standard deviation of the ADC over the repetitions

% Author: Darryl McClymont <darryl.mcclymont@gmail.com>
% Copyright � 2014 University of Oxford
% Version: 0.1.2
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2,6);
nargoutchk(0, 3);

sz = size(I);

if nargin < 3
    Mask = ones(sz(1:end-1));
end

if nargin < 4
    Nreps = 1000;
end

% save memory and time by fitting in vector form
I = reshape(I, [prod(sz(1:end-1)), sz(end)]);

if isscalar(Mask)
    Mask = I(:,1) >= Mask;
end

Ivector = I(Mask(:), :);

% if we haven't been given a ground truth, fit the tensor
if nargin < 5
    [DT, ~, ~, VectorField, ~] = fit_DT(Ivector, bval);
end

% Model fitted image
Ifit = dt2image(DT, bval);
% Residuals
Resids = Ivector - Ifit;

% DT_reps = zeros([size(DT), Nreps]);
% FA_reps = zeros([size(FA), Nreps]);
% ADC_reps = zeros([size(ADC), Nreps]);
% VectorField_reps = zeros([size(VectorField), Nreps]);
% EigVals_reps = zeros([size(EigVals), Nreps]);

% for N iterations, randomly multiply the residuals by 1 or -1 and compute
% parameters
for n = Nreps:-1:1
    if rem(n, 10) == 0
        fprintf('%d, ', n);
    end

    % either 1 or -1, with 50% probability each (Rademacher)
    F = rand(size(Resids)) > 0.5;
    F = F * 2 - 1;

    Resids_to_add = Resids .* F;

    Inew = Ifit + Resids_to_add;

    [~, FA2, ADC2, VectorField2, ~] = fit_DT(Inew, bval);

    VectorField2 = real(VectorField2);

    %DT_reps(:,:,n) = DT2;
    FA_reps(:,1,n) = FA2;
    ADC_reps(:,1,n) = ADC2;
    VectorField_reps(:,:,:,n) = VectorField2;
    %EigVals_reps(:,:,n) = EigVals2;

end

fprintf('done.\n')

FA_STD = zeros(size(Mask));
FA_STD(Mask) = std(FA_reps, [], 3);

ADC_STD = zeros(size(Mask));
ADC_STD(Mask) = std(ADC_reps, [], 3);

% angle between original data and bootstrapped data
Angle_deviation_primary = zeros(size(FA_reps));
Angle_deviation_secondary = zeros(size(FA_reps));
Angle_deviation_tertiary = zeros(size(FA_reps));

for n = 1:Nreps
    % primary eigenvectors
    v1 = real(squeeze(VectorField(:,:,1)));
    v2 = real(squeeze(VectorField_reps(:,:,1,n)));

    % ensure unit magnitude
    v1 = bsxfun(@rdivide, v1, sqrt(sum(v1.^2, 2))+eps);
    v2 = bsxfun(@rdivide, v2, sqrt(sum(v2.^2, 2))+eps);

    theta = dot(v1, v2, 2);
    Angle_deviation_primary(:,1,n) = acos(theta) / pi * 180;

    % secondary eigenvectors
    v1 = real(squeeze(VectorField(:,:,2)));
    v2 = real(squeeze(VectorField_reps(:,:,2,n)));

    % ensure unit magnitude
    v1 = bsxfun(@rdivide, v1, sqrt(sum(v1.^2, 2))+eps);
    v2 = bsxfun(@rdivide, v2, sqrt(sum(v2.^2, 2))+eps);

    theta = dot(v1, v2, 2);
    Angle_deviation_secondary(:,1,n) = acos(theta) / pi * 180;

    % tertiary eigenvectors
    v1 = real(squeeze(VectorField(:,:,3)));
    v2 = real(squeeze(VectorField_reps(:,:,3,n)));

    % ensure unit magnitude
    v1 = bsxfun(@rdivide, v1, sqrt(sum(v1.^2, 2))+eps);
    v2 = bsxfun(@rdivide, v2, sqrt(sum(v2.^2, 2))+eps);

    theta = dot(v1, v2, 2);
    Angle_deviation_tertiary(:,1,n) = acos(theta) / pi * 180;

end

% bigger than 90 degrees? flip it over
Angle_deviation_primary(Angle_deviation_primary > 90) = 180 - Angle_deviation_primary(Angle_deviation_primary > 90);
Angle_deviation_secondary(Angle_deviation_secondary > 90) = 180 - Angle_deviation_secondary(Angle_deviation_secondary > 90);
Angle_deviation_tertiary(Angle_deviation_tertiary > 90) = 180 - Angle_deviation_tertiary(Angle_deviation_tertiary > 90);

% 95th percentile
COU_primary = zeros(size(Mask));
COU_primary(Mask) = prctile(Angle_deviation_primary, 95, 3);
COU_secondary = zeros(size(Mask));
COU_secondary(Mask) = prctile(Angle_deviation_secondary, 95, 3);
COU_tertiary = zeros(size(Mask));
COU_tertiary(Mask) = prctile(Angle_deviation_tertiary, 95, 3);

n = ndims(COU_primary);
if (n == 2) && (size(COU_primary,2) == 1), n = 1; end % ndims gives a vector 2

COU = cat(n+1, COU_primary, COU_secondary, COU_tertiary);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function estimation_results = se3EKF_SLAM(data, NumberOfSteps)
% R-EKF SLAM

% load pre-given data: odometry and observations
if nargin < 1
    load('./data.mat');
end

data_matrix = data.state;

odom_sigma = data.odom_sigma;
obsv_sigma = data.obsv_sigma;

% odoCov = data.odom_cov;   % constant variable
% obsCov = data.obse_cov;   % constant variable

%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% In every step, all elements of Estimation_X will be changed %%%%%%%%%%%%
Estimation_X.orientation = data.poses.orientation(1:3,1:3);
Estimation_X.position    = data.poses.position(:,1);
Estimation_X.cov         = sparse(6,6);
Estimation_X.landmarks   = [];       % the landmarks observed until this step (included), 4*N format, the 4-th row is the index
Estimation_X0.IndexObservedNew=[];
Estimation_X0.IndexObservedAlreadyThis=[];
%%%%%%%%%%%%%%%%%%%% Estimation_X is used to save the state in each step %%%%%%%%%%%%%%%%%%%%

% Initialize
if nargin < 2
    NumberOfSteps = max(data_matrix(:,4));  % step instead of pose,  hence, it does not include pose 0
elseif NumberOfSteps > max(data_matrix(:,4))
    NumberOfSteps = max(data_matrix(:,4));
end
estimation_results = cell(1, NumberOfSteps+1);
estimation_results{1} = Estimation_X;
row_idx = (data_matrix(:, end) <= NumberOfSteps+1);
data_matrix = data_matrix(row_idx, :);

for i = 0:NumberOfSteps
    IndexOfCurrentStepInDataMatrix = find(data_matrix(:,4) == i);
    m = size(IndexOfCurrentStepInDataMatrix, 1);
    if ( mod(i, 50) == 0 )
        disp(['Processing pose ', int2str(i)]);
    end
    % det(Estimation_X.cov)
    if i==NumberOfSteps-1
    a=1;
    end

    if i ~= NumberOfSteps
        OdometryFromThis2Next = data_matrix(IndexOfCurrentStepInDataMatrix(m-5):IndexOfCurrentStepInDataMatrix(m),1);
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 1 ),...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(m-6) , 3 )];
           [Estimation_X] = se3EKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end

        estimation_results{i+1} = Estimation_X;

%        propagation using odometry info
        [Estimation_X] = se3EKF_propagate(Estimation_X, OdometryFromThis2Next, odom_sigma );

    else
        a=2;
        if m > 6
            CameraMeasurementThis = [ data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 1 ) , ...
                                      data_matrix( IndexOfCurrentStepInDataMatrix(1): IndexOfCurrentStepInDataMatrix(end) , 3 )];
            [Estimation_X] = se3EKF_update(Estimation_X, CameraMeasurementThis, obsv_sigma );
        end
        estimation_results{i+1} = Estimation_X;
    end
end
clearvars -except estimation_results

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function varargout = interp1_internal(varargin)
%INTERP1_INTERNAL (overloaded)

switch class(varargin{1})

    case 'double'
        if isequal(varargin{4},'graph') || isequal(varargin{4},'lp') || isequal(varargin{4},'milp') || isequal(varargin{4},'sos2')
            varargout{1} = interp1(varargin{2},varargin{3},varargin{1},'linear');
        else
            varargout{1} = interp1(varargin{2},varargin{3},varargin{1},varargin{4});
        end

    case 'char'

        t = varargin{2};
        x = varargin{3};
        xi = varargin{4};
        yi = varargin{5};
        method = varargin{6};
        dyi = varargin{7};

        if strcmpi(varargin{1},'graph')
            % Convexity propagation wants epi-graph models, so is that what
            % we have used for modelling?
            if strcmpi(method,'envelope')
                   [xi,yi] = convhullprune(xi,yi);
                   [Model,mono,def] = convexGraph(xi,yi,x,t,dyi);
                   operator = struct('convexity','convex','monotonicity',mono,'definiteness',def,'model','graph');
            elseif strcmpi(method,'graph') || strcmpi(method,'lp')
                if isconvexdata(xi,yi)
                    % Convex case, create epi-graph
                    [Model,mono,def] = convexGraph(xi,yi,x,t,dyi);
                    operator = struct('convexity','convex','monotonicity',mono,'definiteness',def,'model','graph');
                elseif isconvexdata(xi,-yi)
                    [Model,mono,def] = concaveGraph(xi,yi,x,t,dyi);
                    operator = struct('convexity','concave','monotonicity',mono,'definiteness',def,'model','graph');
                else
                    Model = [];
                    operator = [];
                end
            else
                Model = [];
                operator = [];
            end
        else
            if strcmpi(method,'sos2') || strcmpi(method,'milp') || strcmpi(method,'lp')
                lambda = sdpvar(length(xi),1);
                Model = [sos2(lambda), x == lambda'*xi(:), t == lambda'*yi(:),lambda>=0, sum(lambda)==1];
                [mono,def] = classifyData(yi);
                operator = struct('convexity','none','monotonicity',mono,'definiteness',def,'model','integer');
            elseif strcmpi(method,'graph')
                Model = [];
                operator = [];
            else
                Model = [min(xi) <= varargin{3} <= max(xi)];
                operator = struct('convexity','none','monotonicity','none','definiteness','none','model','callback');
                operator.bounds = @bounds;
            end
        end

        varargout{1} = Model;
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/INTERP1_INTERNAL called with strange argument');
end

function [L,U] = bounds(xL,xU,varargin)

xv = varargin{1};
yv = varargin{2};
if xL <= xv(1)
    index1 = 1;
else
    index1 = find(xL > xv);index1 = max(1,max(index1));
end
if xU >= xv(end)
    index2 = length(xv);
else
    index2 = find(xv < xU);index2 = min(length(xv),max(index2)+1);
end

if isequal(varargin{3},'linear')
    Z = yv(index1:index2);
    L = min(yv(index1:index2));
    U = max(yv(index1:index2));
else
    N = ceil((xv(index2)-xv(index1))/(mean(diff(xv))/100));
    z = linspace(xv(index1),xv(index2),N);
    yz = interp1(xv,yv,z,varargin{3});
    % To account for finite grid, we add a precision dependent margin
    dz = (z(2)-z(1));
    L = min(yz)-dz;
    U = max(yz)+dz;
end

function isconvex = isconvexdata(xi,yi)
finterp = (yi(1:end-2) + yi(3:end))/2;
if all(finterp >= yi(2:end-1))
    isconvex = 1;
else
    isconvex = 0;
end

function [Model,mono,def] = convexGraph(xi,yi,x,t,dyi)

Model = [min(xi) <= x <= max(xi)];
if isempty(dyi)
    grads = diff(yi)./diff(xi);
    Model = [Model, yi(1:end-1) + grads.*(x - xi(1:end-1)) <= t];
else
    grads = dyi;
    Model = [Model, yi(:) + grads(:).*(x - xi(:)) <= t];
end
[mono,def] = classifyData(yi);

function [Model,mono,def] = concaveGraph(xi,yi,x,t,dyi)

Model = [min(xi) <= x <= max(xi)];
if isempty(dyi)
    grads = diff(yi)./diff(xi);
    Model = [Model, yi(1:end-1) + grads.*(x - xi(1:end-1)) >= t];
else
    grads = dyi;
    Model = [Model, yi(:) + grads(:).*(x - xi(:)) >= t];
end

[mono,def] = classifyData(yi);

function [xout,yout] = convhullprune(xin,yin)
% Prunes data down to the lower convex envelope

% First prune by convex hull
cj = convhulln([xin(:) yin(:)]);
cj = sort(unique(cj));
x = xin(cj);
y = yin(cj);

% This can probably be done in linear complexity, I am just lazy here
keep = ones(length(x),1);
for i = 1:length(x)-1
    if keep(i)
        for j = length(x):-1:3
            if keep(j)
                grad = (y(j)-y(i))/(x(j)-x(i));
                k = i+1:j-1;
                fail =  y(k) > y(i) + grad*(x(k)-x(i));
                keep(k(fail)) = 0;
            end
        end
    end
end
keep = cj(find(keep));
xout = xin(keep);
yout = yin(keep);

function [mono,def] = classifyData(yi);
if all(diff(yi)) >= 0
    mono = 'increasing';
elseif all(diff(yi) <= 0)
    mono = 'decreasing';
else
    mono = 'none';
end
if all(yi>=0)
    def = 'positive';
elseif all(yi <= 0)
    def = 'negative';
else
    def = 0;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP