text:
function planC = createPCimage(scanNum,structNumV,pcaParamsFile,planC)
% function createPCimage(scanNum,structNumV,pcaParamsFile,planC)
%
% APA, 2/4/2017

if ~iscell(scanNum)

    rowMargin = 100; % extend rows by this amount
    colMargin = 512; % extend cols by this amount
    slcMargin = 7; % extend slcss by this amount  %Changed from 15

    if ~exist('planC','var')
        global planC
    end

    indexS = planC{end};

    % get ROI
    randomFlg = 1;
    [volToEval,maskBoundingBox3M,mask3M,minr,maxr,minc,maxc,mins,maxs,uniqueSlices] = ...
        getROI(structNumV,rowMargin,colMargin,slcMargin,planC,randomFlg);

    sliceThickNessV = ...
        [planC{indexS.scan}(scanNum).scanInfo(mins:maxs).sliceThickness];
    [xVals, yVals, zVals] = getScanXYZVals(planC{indexS.scan}(scanNum));

else

    volToEval = scanNum{1};
    maskBoundingBox3M = scanNum{2};
    mask3M = scanNum{3};
    minr = scanNum{4};
    maxr = scanNum{5};
    minc = scanNum{6};
    maxc = scanNum{7};
    mins = scanNum{8};
    maxs = scanNum{9};
    uniqueSlices = scanNum{10};
    sliceThickNessV = scanNum{11};
    xVals = scanNum{12};
    yVals = scanNum{13};
    zVals = scanNum{14};
    rowMargin = [];
    colMargin = [];
    slcMargin = [];
end

minIntensity = -200;   % Clipping min
maxIntensity = 400; % Clipping max
harOnlyFlg = 1;
featFlagsV = [1,1,1,1,1,1,1,1,1];
numLevsV = 64;
patchRadiusV = [1,2];
compNumV = 1:2;

load(pcaParamsFile)

% Clip intensities
% volToEval(volToEval < minIntensity) = minIntensity;
% volToEval(volToEval > maxIntensity) = maxIntensity;
nanIntenityV = volToEval < -400;

%%% AA commented for writing cropped CT scans

% Get Law's and Haralick features for this structure
newFeaturesM = getLawsAndHaralickFeatures({volToEval,maskBoundingBox3M},...
    rowMargin,colMargin,slcMargin,minIntensity,maxIntensity,planC);
% rowMargin = [];
% colMargin = [];
% slcMargin = [];
% newFeaturesM = gpuGetLawsAndHaralickFeatures({volToEval,maskBoundingBox3M},...
%         rowMargin,colMargin,slcMargin,minIntensity,maxIntensity,planC,...
%         harOnlyFlg,featFlagsV,numLevsV,patchRadiusV);
newFeaturesM = bsxfun(@minus,newFeaturesM, featureMeanV);
newFeaturesM = bsxfun(@rdivide,newFeaturesM, featureStdV);
%newFeaturesM = newFeaturesM(:,indSignifV);

siz = size(maskBoundingBox3M);

% Parameters for storing component image to planC
deltaXYZv = [abs(xVals(1)-xVals(2)) abs(yVals(1)-yVals(2)) abs(zVals(1)-zVals(2))];
zV = zVals(mins:maxs);
regParamsS.horizontalGridInterval = deltaXYZv(1);
regParamsS.verticalGridInterval   = deltaXYZv(2); %(-)ve for dose
regParamsS.coord1OFFirstPoint   = xVals(minc);
%regParamsS.coord2OFFirstPoint   = yVals(minr); % for dose
regParamsS.coord2OFFirstPoint   = yVals(maxr);
regParamsS.zValues  = zV;
regParamsS.sliceThickness = sliceThickNessV;
assocTextureUID = '';
for compNum = compNumV

    %% Plot components, raw image and segmentation
    compV = zeros(prod(siz),1);

    compV(~nanIntenityV) = newFeaturesM * coeff(:,compNum) * 1.0 + newFeaturesM * coeff(:,2) * 0.0 + ...
        newFeaturesM * coeff(:,3) * 0.0 + newFeaturesM * coeff(:,4) * 0.0;

    %compV = max(compV) - compV; % reverse intensities

    compV(nanIntenityV) = min(compV);
    comp3M = reshape(compV,siz);

    % figure, imagesc(comp3M(:,:,30))
    % figure, imagesc(volToEval(:,:,30))

    %%% AA commented for writing cropped CT scans ends

    %%
    % Save comp3M to planC
    %dose2CERR(entropy3M,[], 'entropy3voxls_Ins3_NI14','test','test','non CT',regParamsS,'no',assocScanUID)
    planC = scan2CERR(comp3M,['PC_',num2str(compNum)],'Passed',regParamsS,assocTextureUID,planC);
    % planC = scan2CERR(volToEval,'PC','Passed',regParamsS,assocTextureUID,planC);

end

% for iFeat = 1:size(newFeaturesM,2)
%     compV = zeros(prod(siz),1);
%     compV(nanIntenityV) = min(compV);
%     compV(~nanIntenityV) = newFeaturesM(:,iFeat);
%         planC = scan2CERR(reshape(compV,siz),['Feat_',num2str(iFeat)],'Passed',regParamsS,assocTextureUID,planC);
% end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Ex0301: Perform Surface Meshes (CALD)')

SpharmMatDir = '../';
CodeDir = [SpharmMatDir 'code'];
addpath(CodeDir);

DataFolder = 'data';

%% (1) Input data directory
InDataDir = [SpharmMatDir DataFolder '/Ex0301/mesh01_obj'];

%% (2) List input file names
inFiles = dir([InDataDir '/*_obj.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [InDataDir '/' inFiles(i).name];
end

%% (3) Display input objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'none';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

%% (4) Perform Surface Meshes (CALD)
confs.MeshGridSize = 50;
confs.MaxSPHARMDegree = 6;
confs.Tolerance = 2;
confs.Smoothing = 2;
confs.Iteration = 100;
confs.LocalIteration = 10;
    % Available values for t_major- 'x';'y'
confs.t_major = 'x';
    % Available values for SelectDiagonal- 'ShortDiag';'LongDiag'
confs.SelectDiagonal = 'ShortDiag';

confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0301/mesh02_smo'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

outNames = SpharmMatParameterization(confs, inNames, 'ParamCALD');

%% (5) Display output objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'both';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'orig';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'adc_paramap';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, outNames, CodeDir);

inFiles = dir([confs.OutDirectory '/initParamCALD/*.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [confs.OutDirectory '/initParamCALD/' inFiles(i).name];
end
SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

rmpath(CodeDir);

clear;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function linpack_s_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests SGBFA and SGBSL.
%
%  Discussion:
%
%    The problem solved here is a larger version of this one:
%
%    Matrix A is ( 2 -1  0  0  0)    right hand side B is  (1)
%                (-1  2 -1  0  0)                          (0)
%                ( 0 -1  2 -1  0)                          (0)
%                ( 0  0 -1  2 -1)                          (0)
%                ( 0  0  0 -1  2)                          (1)
%
%
%    Solution is   (1)
%                  (1)
%                  (1)
%                  (1)
%                  (1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    N is the number of equations.
%
%    ML is the number of subdiagonals,
%    MU the number of superdiagonals.
%
%    LDA is the leading dimension of the array used to store the
%    matrix, which must be at least 2*ML+MU+1.
%
  n = 10;
  ml = 1;
  mu = 1;
  lda = 2*ml+mu+1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  SGBFA factors the matrix,\n' );
  fprintf ( 1, '  SGBSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Set the right hand side B.
%
  b(1) = 1.0;
  b(2:n-1) = 0.0;
  b(n) = 1.0;
%
%  Force B to be a column vector.
%
  b = b';
%
%  Set the matrix A.
%
  m = ml + mu + 1;
  fprintf ( 1, '  The bandwidth of the matrix is %d\n', m );

  for j = 1 : n
    a(m-1,j) = -1.0;
    a(m,j) = 2.0;
    a(m+1,j) = -1.0;
  end
%
%  Call SGBFA to factor the matrix A.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipivot, info ] = sgbfa ( a, lda, n, ml, mu );

  if ( info ~= 0 )
    fprintf ( 1, '  Error!  SGBFA returns INFO = %d\n', info );
    return
  end
%
%  Call SGBSL to solve the linear system.  The solution
%  is returned in B, that is, it overwrites the right hand side.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  job = 0;
  b = sgbsl ( a, lda, n, ml, mu, ipivot, b, job );
%
%  Print the results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last 5 entries of the solution:\n' );
  fprintf ( 1, '  (All should be 1):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = r8_choose ( n, k )

%*****************************************************************************80
%
%% R8_CHOOSE computes the combinatorial coefficient C(N,K).
%
%  Discussion:
%
%    Real arithmetic is used, and C(N,K) is computed directly, via
%    Gamma functions, rather than recursively.
%
%    C(N,K) is the number of distinct combinations of K objects
%    chosen from a set of N distinct objects.  A combination is
%    like a set, in that order does not matter.
%
%  Example:
%
%    The number of combinations of 2 things chosen from 5 is 10.
%
%    C(5,2) = ( 5 * 4 * 3 * 2 * 1 ) / ( ( 3 * 2 * 1 ) * ( 2 * 1 ) ) = 10.
%
%    The actual combinations may be represented as:
%
%      (1,2), (1,3), (1,4), (1,5), (2,3),
%      (2,4), (2,5), (3,4), (3,5), (4,5).
%
%  Formula:
%
%    C(N,K) = N! / ( (N-K)! * K! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the value of N.
%
%    Input, integer K, the value of K.
%
%    Output, real VALUE, the value of C(N,K)
%
  n = floor ( n );
  k = floor ( k );

  if ( n < 0 )

    value = 0.0;

  elseif ( k == 0 )

    value = 1.0;

  elseif ( k == 1 )

    value = n;

  elseif ( 1 < k && k < n-1 )

    facn = gammaln ( n + 1 );
    fack = gammaln ( k + 1 );
    facnmk = gammaln ( n - k + 1 );

    value = round ( exp ( facn - fack - facnmk ) );

  elseif ( k == n-1 )

    value = n;

  elseif ( k == n )

    value = 1.0;

  else

    value = 0.0;

  end

  return
end


text:
classdef prtKernelStep < prtKernel

    properties (SetAccess = private)
        name = 'Kernel Step'; % RBF Kernel
        nameAbbreviation = 'Step'; % RBF
    end

    properties
        sigma = 1; % The inverse kernel width
    end

    methods (Access = protected, Hidden = true)
        function Obj = trainAction(Obj,ds)
            Obj.internalDataSet = ds;
            Obj.isTrained = true;
        end

        function dsOut = runAction(Obj,ds)
            if ~Obj.isTrained
                error('prtKernelRbf:run','Attempt to run an untrained kernel; use kernel.train(ds) to train');
            end
            if Obj.internalDataSet.nObservations == 0
                dsOut = prtDataSetClass;
            else
                gram = prtKernelStep.kernelFn(ds.getObservations,Obj.internalDataSet.getObservations);
                dsOut = ds.setObservations(gram);
            end
        end
    end

    methods
        function Obj = prtKernelStep(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

    end

    methods(Hidden = true)
        function varargout = plot(obj)
            x = obj.internalDataSet.getObservations;

            if size(x,2) <= 3
                if size(x,2) == 1 && obj.internalDataSet.isLabeled
                    xy = cat(2,x,obj.internalDataSet.getTargets);
                    h = prtPlotUtilScatter(xy, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                else
                    h = prtPlotUtilScatter(x, {}, obj.plotOptions.symbol, obj.plotOptions.markerFaceColor, obj.plotOptions.color, obj.plotOptions.symbolLineWidth, obj.plotOptions.symbolSize);
                end
            else
                h = nan;
            end

            varargout = {};
            if nargout
                varargout = {h};
            end
        end
    end

    methods (Static, Hidden = true)
        function gram = kernelFn(x,y)
            [n1, d] = size(x);
            [n2, nin] = size(y);
            if d ~= nin
                error('size(x,2) must equal size(y,2)');
            end

            gram = true(n1,n2);
            for dim = 1:d;
                gram = gram & bsxfun(@(x,y) x > y,x(:,d),y(:,d)');
            end
            gram = double(gram);
%             gram = imfilter(gram,fspecial('gaussian',[11,1],4));

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [CORR, TABLE] = CorrTable(x,vnames,pairwise)
% =======================================================================
% Computes the correlation matrix of a matrix of time series with titles
% =======================================================================
% out = CorrTable(x,vnames)
% -----------------------------------------------------------------------
% INPUT
%   - x: matrix (rows x cols) [double]
%   - vnames: array (1 x cols) with names of the (cols) series [cell]
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - pairwise: default 0, change to 1 to compute pairwise correlations
% -----------------------------------------------------------------------
% OUTPUT
%   - CORR: correlation matrix [double]
%   - TABLE: table of correlation matrix with titles [cell]
% -----------------------------------------------------------------------
% EXAMPLE
%   x = rand(50,2);
%   [CORR, TABLE] = CorrTable(x,{'Consumption','Investment'})
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

% If no dimension is specified, set it to 1
if ~exist('pairwise','var')
    pairwise = 0;
end

% vnames must be a row vector
if size(vnames,1)>1
    vnames = vnames';
end

% Compute correlation
if pairwise==0
    CORR = corr(x);
elseif pairwise==1
    CORR = corr(x,'rows','pairwise');
end
[r, ~] = size(CORR);

% Place NaNs on the upper triangular
CORR(find(triu(CORR))) = NaN;
jj = 1;
for ii=1:r
    CORR(ii,jj) = 1;
    jj = jj +1;
end

% Save position of NaNs
nans = isnan(CORR);

% Transform table into cell array
TAB = num2cell(CORR);

% Add vnames on top
TAB = [vnames ; TAB];

% Add vnames on the left
aux = [{''} vnames]';
TABLE = [aux TAB];

% Print using mprint
info.cvnames = char(vnames);
info.rvnames = char([{' '} vnames]);
mprint(CORR,info);

author: ambropo

repo: VAR-Toolbox

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


text:
function val = rmGetVoxelData(param, vw, roi, varargin)
% rmGetVoxelData - Get Retinotopy Model data for a selected set of voxels.
%
% val = rmGetVoxelData(param, [view or RM model], [voxels=cur ROI], [options]);
%
%
% INPUTS:
%   param: name of the parameter to get for each voxel. Possible values:
%       'prf', 'rf': constructed 2D population receptive fields for each voxel.
%       If several voxels are specified, returns a 3D matrix in which the
%       slices correspond to the pRF for each voxel.
%
%		'prfvector', 'rfvector': same as 'prf', but each voxel's RF is
%		provided in a single vector format. The resulting matrix is size
%		pixels (in visual space) x voxels.
%
%       'x', 'x0': X position center of the pRF for the voxel.
%
%       'y', 'y0': Y position center of the pRF for the voxel.
%
%		'ecc', eccentricity of pRF center.
%
%		'pol', polar angle of pRF center.
%
%       'sigma', 'sigmamajor': pRF size of the first/only Gaussian.
%
%       'sigmaminor': pRF size of the second Gaussian, if it exists.
%
%       'theta': relative rotation of the major/minor axes.
%
%       'sig', 'log10p': significance level (-log10(p)) for the model for
%        each voxel.
%
%       'beta': scaling coefficient ("beta value") for the Gaussian pRF for
%       each voxel.
%
%		'amp': amplitude of the response. This is the scale factor
%		necessary to take a predicted response for each voxel, and fit it
%		to the time series. (This is the "beta value" for a GLM, constructed
%		from a pRF-derived esign matrix).
%
%
%   model: can either be a retinotopy model struct, or a mrVista view.
%   If the latter, grabs the first model loaded into the view (prompting
%   the user to select and RM file if none are loaded). Can also be
%   a number, indexing into loaded model structs. [Default: cur view
%   model]
%
%   voxels: ROI specification. Can be a numeric index into the view's ROIs
%   field, the ROI struct itself, a 3xN set of coordinates relative to the
%   view, or the ROI name of a loaded ROI. (Uses tc_roiStruct to
%   disambiguate these possibilities.) [Default: cur ROI of cur view]
%
% options include:
%   'rot', [value]: rotate RF params clockwise by [value] degrees. Ths will
%                   happen if the stimulus specification doesn't quite
%                   match the actual presentation regime. (Shouldn't
%                   happen, I know, but it seems to for some data sets, to
%                   a small degree.)
%   'xRange', [value]: range of X values for estimating the pRF. Defaults
%                   to -14:.2:14.
%   'yRange', [value]: same as for xRange, but along the Y axis.
%
%
% ras, 10/2006.
if notDefined('param'),     error('Need to specify param name.');	end
if notDefined('vw'),        vw = getCurView;						end
if notDefined('roi'),       roi = viewGet(vw, 'curROI');            end

% params/defaults
rot = 0;

% parse options
varargin = unNestCell(varargin); % allow recursive passing of options
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch lower(varargin{i})
            case 'rot', rot = varargin{i+1};
            case 'x', X = varargin{i+1};
            case 'y', Y = varargin{i+1};
            case 'rmparams', rmParams = varargin{i+1};
        end
    end
end

% test for a view struct being passed in
if ~isfield(vw, 'x0') && isfield(vw, 'viewType')
    if ~checkfields(vw, 'rm', 'retinotopyModels')
        vw = rmSelect(vw, 1);
    end
    rmParams = vw.rm.retinotopyParams;
    model = vw.rm.retinotopyModels{1};
else
    % this is messy, need to reorganize...
    model = vw;
    vw = getCurView;
end

% disambiguate ROI sepecification
if ~isstruct(roi), roi = tc_roiStruct(vw, roi); end

% get indices of coords which correspond to selected voxels
if ismember(roi.viewType, {'Volume' 'Gray'})
    % get the indices corresponding to the ROI, preserving the voxel order
    % (this may take some time and memory for large ROIs)
    I = roiIndices(vw, roi.coords, 1);

    ok = find( ~isnan(I) );

else
    % not yet implemented for other view types
    error('Sorry, Not Yet Implemented.')
end

% plural/singular flexibility: ignore any 's' at the end of the param name
if lower(param(end))=='s'
    param = param(1:end-1);
end

% get the relevant value for the specified parameter
val = NaN(1, length(I));
switch lower(param)
    case {'amp' 'voxelamplitude' 'voxamp' 'glmbeta'}
        % compute the scale factor which maps from a normalized (-1 - 1
        % range) prediction to the fitted data. This runs a GLM for each
        % voxel in the ROI.
        %
        % I'm going to go ahead and assume vw is a view structure; it's too
        % hard to do this otherwise.
        [tSeries coords] = voxelTSeries(vw, roi.coords, [], 0, 0);

        pred = rmGetVoxelData('predictedtseries', vw, roi, varargin);

        % we called voxelTSeries with the preserveCoords flag set to zero.
        % This means the tSeries will have some voxels discarded (if
        % there's no data), and the coord order shuffled. Make sure the
        % order of predictions matches this order.
        [coords ok map2roi] = intersectCols(roi.coords, coords); %#ok<*ASGLU>
        pred = pred(:,ok);

        % take the prediction and make a design matrix (add trends terms)
        trends = rmMakeTrends(rmParams, prefsVerboseCheck);

        % initialize an empty output value
        val = NaN(1, length(I));

        % run a GLM for each non-NaN voxel, returning the scale factor
        for v = 1:length(ok)
            p = pred(:,v);
            if all(p==0) || all( isnan(p) )
                continue
            end
            p = p ./ abs(max(p));
            X = [p trends];
            [t df RSS B] = rmGLM(tSeries(:,v), X);

            % because of the sorting we introduced in voxelTSeries, the
            % order of tSeries may not match the order of our output amps
            % value. Map this back to the original order of the ROI coords.
            whichVoxel = map2roi(v);
            val(whichVoxel) = B(1);
        end

    case {'pred' 'prediction' 'predictedtserie'}
        % predicted time series of response based on the pRF for each voxel
        % and the stimulus.

        % we need pRF params for each voxel
        X = rmParams.analysis.X;
        Y = rmParams.analysis.Y;
        x0 = zeros(1, length(I));
        y0 = x0;
        sigma = x0;
        beta = x0;
        x0(ok) = model.x0(I(ok));
        y0(ok) = model.y0(I(ok));
        sigma(ok) = model.sigma.major(I(ok));
        beta(ok) = model.beta(1,I(ok),1);

        % rotation compensation if selected
        if rot ~= 0
            R = sqrt(x0.^2 + y0.^2);
            theta = atan2(y0, x0);
            theta = theta - deg2rad(rot);
            theta = mod(theta, 2*pi);
            x0 = R .* cos(theta);
            y0 = R .* sin(theta);
        end

        % initalize an empty output matrix
        nFrames = size(rmParams.analysis.allstimimages, 1);
        val = NaN(nFrames, length(I));

        for v = ok
            pRF = rfGaussian2D(X, Y, sigma(v), sigma(v), 0, x0(v), y0(v));
            val(:,v) = rmParams.analysis.allstimimages * pRF;
        end

    case {'prf' 'rf' 'prfs' 'rfs'}
        % we need to find a sampling grid for the pRF.
        if notDefined('X') || notDefined('Y')
            [X Y] = prfSamplingGrid(rmParams);
        end

        x0 = zeros(1, length(I));
        y0 = x0;
        sigma = x0;
        beta = x0;
        x0(ok) = model.x0(I(ok));
        y0(ok) = model.y0(I(ok));
        sigma(ok) = model.sigma.major(I(ok));
        beta(ok) = model.beta(1,I(ok),1);

        % rotation compensation if selected
        if rot ~= 0
            R = sqrt(x0.^2 + y0.^2);
            theta = atan2(y0, x0);
            theta = theta - deg2rad(rot);
            theta = mod(theta, 2*pi);
            x0 = R .* cos(theta);
            y0 = R .* sin(theta);
        end

        % re-initialize the return value as a matrix of NaNs
        val = NaN(size(X, 1), size(X, 2), length(I));

        for v = ok
            pRF = rfGaussian2D(X, Y, sigma(v), sigma(v), 0, x0(v), y0(v));
            val(:,:,v) = pRF;
        end

    case {'prfvector' 'rfvector' 'prfsvector' 'rfsvector' 'prfvectors' 'rfvectors'}
        % we need to find a sampling grid for the pRF.
        if notDefined('X') || notDefined('Y')
            [X Y] = prfSamplingGrid(rmParams);
        end

        x0 = zeros(1, length(I));
        y0 = x0;
        sigma = x0;
        beta = x0;
        x0(ok) = model.x0(I(ok));
        y0(ok) = model.y0(I(ok));
        sigma(ok) = model.sigma.major(I(ok));
        beta(ok) = model.beta(1,I(ok),1);

        % rotation compensation if selected
        if rot ~= 0
            R = sqrt(x0.^2 + y0.^2);
            theta = atan2(y0, x0);
            theta = theta - deg2rad(rot);
            theta = mod(theta, 2*pi);
            x0 = R .* cos(theta);
            y0 = R .* sin(theta);
        end

        % re-initialize the return value as a matrix of NaNs
        val = NaN(numel(X), length(I));

        for v = ok
            pRF = rfGaussian2D(X(:), Y(:), sigma(v), sigma(v), 0, x0(v), y0(v));
            val(:,v) = pRF;
        end

    case {'x' 'x0'}, val(ok) = model.x0(I(ok));

    case {'y' 'y0'}, val(ok) = model.y0(I(ok));

    case {'sigma' 'sigmamajor'}, val(ok) = model.sigma.major(I(ok));

    case {'sigmaminor'}, val(ok) = model.sigma.minor(I(ok));

    case {'theta'}, val(ok) = model.sigma.theta(I(ok));

    case {'beta'}, val(ok) = model.beta(1,I(ok),1);

    case {'sig' 'log10p' 'logp'}
        allsig = rmGet(model, 'log10p');
        val(ok) = allsig(I(ok));

    otherwise,
        try
            allvals = rmGet(model, param);
            val(ok) = allvals(I(ok));
        catch %#ok<*CTCH>
            fprintf('[%s] Unknown parameter name %s \n', mfilename, param);
        end

end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%MODESEEK Clustering by mode-seeking
%
% 	[LAB,J] = MODESEEK(D,K)
%
% INPUT
%   D       Distance matrix or distance dataset (square)
%   K       Number of neighbours to search for local mode (default: 10)
%
% OUTPUT
%   LAB     Cluster assignments, 1..K
%   J       Indices of modal samples
%
% DESCRIPTION
% A K-NN modeseeking method is used to assign each object to its nearest mode.
%
% REFERENCES
% 1. Cheng, Y. "Mean shift, mode Seeking, and clustering", IEEE Transactions
% on Pattern Analysis and Machine Intelligence, vol. 17, no. 8, pp. 790-799,
% 1995.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, PRKMEANS, HCLUST, KCENTRES, PROXM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: modeseek.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function [assign,J] = modeseek (d,k)

		if (nargin < 2)
		prwarning(1,'No k supplied, assuming k = 10');
		k = 10;
  end

	[m,n] = size(d);

  if numel(k) > 1
    assign = zeros(m,numel(k));
    for j = 1:numel(k)
      assign(:,j) = feval(mfilename,d,k(j));
    end
    return
  end

  d = d'; % correction to analyse asymmetric matrices horizontally

	if (m ~= n), error('distance matrix should be square'); end
	if (k < 2),  error('neighborhood size should be at least 2'); end
	if (k > n),  error('k too large for this dataset'); end

	[d,J] = sort(+d,1);		   % Find neighbours.
	f = 1./(d(k,:)+realmin); % Calculate densities.
	J(k+1:end,:) = [];   	   % Just retain indices of neighbours.

	% Find indices of local modes in neighbourhood.
	[dummy,I] = max(reshape(f(J),size(J)));

	% Translate back to indices in all the data. N now contains the
	% index of the nearest neighbour in the K-neighbourhood.
	N = J(I+[0:k:k*(m-1)]);

	% Re-assign samples to the sample their nearest neighbour is assigned to.
	% Iterate until assignments don't change anymore. Samples that then point
	% to themselves are modes; all other samples point to the closest mode.

	M = N(N);
	while (any(M~=N))
		N = M; M = N(N);
	end

	% Use renumlab to obtain assignments 1, 2, ... and the list of unique
	% assignments (the modes).

	[assign,J] = renumlab(M');

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ state ] =movePTP_ConditionalTorque_Arc_AC(t_Kuka,theta,c,k,VEL,joints_indices,max_torque,min_torque)
%% This function is used for moving the end-effector on an arc, for the KUKA iiwa 7 R 800.
% this motion is interruptible when one or more of the joints-torques exceed the predefined limits

%% Syntax
% [ state ] =movePTP_ConditionalTorque_Arc_AC(t_Kuka,theta,c,k,relVel,joints_indices,max_torque,min_torque)

%% About:
% This function is used to move the end-effector on an arc,

%% Arreguments:
% t_Kuka: is the TCP/IP connection
% theta: is the arc angle, in radians.
% c: the x,y,z coordinates of the center of the circle, it is 1x3 vector.
% k: is the normal vector on the plane of the arc, it is 1x3
% vector.
% VEL : is a double, defines the motion velocity mm/sec.
% joints_indices: is a vector of the indices of the joints where the
% torques limits are to be imposed, the joints are indexed starting from one.
% max_torque: is a vector of the maximum torque limits for the joints
% specified in the (joints_indices) vector.
% min_torque: is a vecotr of the minimum torque limits for the joints
% specified in the (joints_indices) vector.

%% Return value:
% state: a number equals to one if the motion is completed sccessfully or a
% zero if the motion was interrupted due to external contact with the
% robot.
% if there is an error, the return value is minus one

% Copy right, Mohammad SAFEEA, 9th of April 2018

if((size(VEL,1)==1)&&(size(VEL,2)==1))
else
    disp('Error, velocity is a double and shall not be an array');
    state=-1;
    return;
end
% Check if the inputs "joints_indices,max_torque,min_torque" are vectors
if(isVec(joints_indices)==0)
    disp('Error, joints_indices shall be a vector')
    state=-1;
    return;
end
if(isVec(max_torque)==0)
    disp('Error, max_torque shall be a vector')
    state=-1;
    return;
end
if(isVec(min_torque)==0)
    disp('Error, min_torque shall be a vector')
    state=-1;
    return;
end
joints_indices=colVec(joints_indices);
max_torque=colVec(max_torque);
min_torque=colVec(min_torque);

% Check if size of vectrs are equal
j_num=size(joints_indices,1);
n=size(min_torque,1);
if(j_num==n)
else
    disp('Error, sizes of vectors "joints_indices" and "min_torque" shall be equal');
    state=-1;
    return;
end
n=size(max_torque,1);
if(j_num==n)
else
    disp('Error, sizes of vectors "joints_indices" and "max_torque" shall be equal');
    state=-1;
    return;
end

c=colVec(c);
k=colVec(k);

pos=getEEFPos( t_Kuka );
p1=[pos{1};pos{2};pos{3}];
p1=colVec(p1);

r=norm(p1-c);

if(or(r==0,theta==0))
 return;
end
normK=norm(k);
if(normK==0)
 fprintf('Norm of direction vector k shall not be zero \n');
 return;
end
k=k/normK;

s=(p1-c)/r;
n=cross(k,s);

c1=r*cos(theta/2)*s+r*sin(theta/2)*n+c;
c2=r*cos(theta)*s+r*sin(theta)*n+c;

f1=pos;
i=1;
f1{i}=c1(i);
i=i+1;
f1{i}=c1(i);
i=i+1;
f1{i}=c1(i);

f2=pos;
i=1;
f2{i}=c2(i);
i=i+1;
f2{i}=c2(i);
i=i+1;
f2{i}=c2(i);

state=movePTP_ConditionalTorque_Circ1OrintationInter( t_Kuka , f1,f2,...
    VEL,joints_indices,max_torque,min_torque);

end

function y=colVec(x)
if(size(x,2)==1)
    y=x;
else
    y=x';
end
end

function y=isVec(x)
y=0;
if(size(x,2)==1)
    y=1;
end
if(size(x,1)==1)
    y=1;
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% SB2_POSTERIORMODE  Posterior mode-finder for the SPARSEBAYES algorithm
%
% [MU, U, BETA, LIKELIHOODMODE, BADHESS] = ...
%    SB2_POSTERIORMODE(LIKELIHOOD,BASIS,TARGETS,ALPHA,MU,ITSMAX,OPTIONS)
%
% OUTPUT ARGUMENTS:
%
%	MU				Parameter values at the mode
%	U				Cholesky factor of the covariance at the mode
%	BETA			Vector of pseudo-noise variances at the mode
%	LIKELIHOODMODE	Data likelihood at the mode
%	BADHESS			Returns true if Hessian is "bad" (becomes
%					non-positive-definite during maximisation)
%
% INPUT ARGUMENTS:
%
%	LIKELIHOOD	LIKELIHOOD structure
%	BASIS		Current relevant basis matrix
%	TARGETS		N-vector with target output values
%	ALPHA		Current hyperparameters
%	MU			Current weights
%	ITSMAX		Maximum number of iterations to run
%	OPTIONS		Standard OPTIONS structure (only affects diagnostics)
%
% NOTES:
%
% SB2_POSTERIORMODE finds the posterior mode (with respect to the
% weights) of the likelihood function in the non-Gaussian case to
% facilitate subsequent Laplace approximation.
%
% This function is intended for internal use by SPARSEBAYES only (within
% SB2_FullStatistics).
%

%
% Copyright 2009, Vector Anomaly Ltd
%
% This file is part of the SPARSEBAYES library for Matlab (V2.0).
%
% SPARSEBAYES 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.
%
% SPARSEBAYES 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 SPARSEBAYES in the accompanying file "licence.txt"; if not, write to
% the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
% MA 02110-1301 USA
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function [Mu, U, Beta, likelihoodMode, badHess] = ...
    SB2_PosteriorMode(LIKELIHOOD,BASIS,Targets,Alpha,Mu,itsMax,OPTIONS)

% TOLERANCES
%
% Termination criterion for each gradient dimension
%
GRADIENT_MIN	= 1e-6;
%
% Minimum fraction of the full Newton step considered
%
STEP_MIN	= 1/(2^8);
%

[N, M]	= size(BASIS);
A		= diag(Alpha);

% NB: for historical reasons, we work in term of error here (negative
% log-liklihood) and minimise

%
% Get current model output and data error
%
BASIS_Mu		= BASIS*Mu; % Linear output
[dataError, y]	= SB2_DataError(LIKELIHOOD, BASIS_Mu, Targets);
%
% Add on the weight penalty
%
regulariser		= (Alpha'*(Mu.^2))/2;
newTotalError	= dataError + regulariser;
%
badHess		= false;
errorLog	= zeros(itsMax,1);
%
for iteration=1:itsMax
  %
  % Log the error value each iteration
  %
  errorLog(iteration)	= newTotalError;
  SB2_Diagnostic(OPTIONS,4,'PM cycle: %2d\t error: %.6f\n', ...
				 iteration, errorLog(iteration));

  % Construct the gradient
  %
  e	= (Targets-y);
  g	= BASIS'*e - Alpha.*Mu;
  %
  % Compute the likelihood-dependent analogue of the noise precision.
  % NB: Beta now a vector.
  %
  switch LIKELIHOOD.InUse
   case LIKELIHOOD.Bernoulli
    Beta	= y.*(1-y);
   case LIKELIHOOD.Poisson
    Beta	= y;
  end
  % Compute the Hessian
  BASIS_B	= BASIS .* (Beta * ones(1,M));
  H			= (BASIS_B'*BASIS + A);
  % Invert Hessian via Cholesky, watching out for ill-conditioning
  [U, pdErr]	= chol(H);
  % Make sure its positive definite
  if pdErr
	% If you see this, it's *probably* the result of a bad choice of
    % basis. e.g. a kernel with too large a "width"
	%
    SB2_Diagnostic(OPTIONS, 1,...
				   '** Warning ** Ill-conditioned Hessian (%g)\n', ...
				   rcond(H));
    badHess			= true;
    U				= [];
    Beta			= [];
    likelihoodMode	= [];
    return
  end
  %
  % Before progressing, check for termination based on the gradient norm
  %
  if all(abs(g)<GRADIENT_MIN)
    errorLog	= errorLog(1:iteration);
    SB2_Diagnostic(OPTIONS,4,['PM convergence (<%g) after ' ...
					'%d iterations, |g| = %g\n'], ...
				   GRADIENT_MIN,iteration,max(abs(g)));
    break
  end

  % If all OK, compute full Newton step: H^{-1} * g
  %
  DeltaMu	= U \ (U' \ g);
  step		= 1;
  %
  while step>STEP_MIN
	% Follow gradient to get new value of parameters
	%
    Mu_new		= Mu + step*DeltaMu;
    BASIS_Mu	= BASIS*Mu_new;
	%
    % Compute outputs and error at new point
	%
    [dataError,y]	= SB2_DataError(LIKELIHOOD,BASIS_Mu,Targets);
    regulariser		= (Alpha'*(Mu_new.^2))/2;
    newTotalError	= dataError + regulariser;
	%
    % Test that we haven't made things worse
	%
    if newTotalError>=errorLog(iteration)
      % If so, back off!
	  %
      step	= step/2;
      SB2_Diagnostic(OPTIONS, 4,['PM error increase! Backing off to l=' ...
					'%.3f\n'], step);
    else
      Mu	= Mu_new;
      step	= 0;			% this will force exit from the "while" loop
    end
  end
  %
  % If we get here with non-zero "step", it means that the smallest
  % offset from the current point along the "downhill" direction did not
  % lead to a decrease in error. In other words, we must be
  % infinitesimally close to a minimum (which is OK).
  %
  if step
    SB2_Diagnostic(OPTIONS, 4, ...
				   'PM stopping due to back-off limit (|g| = %g)\n', ...
				   max(abs(g)));
    break;
  end
end
%
% Simple computation of return value of log likelihood at mode
%
likelihoodMode	= -dataError;

%%%%%%%%%%%%%%%%%%%%%%
%
% Support function
%
%%%%%%%%%%%%%%%%%%%%%%

function [e,y] = SB2_DataError(LIKELIHOOD,BASIS_Mu,Targets)

switch LIKELIHOOD.InUse
 case LIKELIHOOD.Bernoulli
  y	= SB2_Sigmoid(BASIS_Mu);
  % Handle probability zero cases
  y0	= (y==0);
  y1	= (y==1);
  if any(y0(Targets>0)) || any(y1(Targets<1))
    % Error infinite when model gives zero probability in
    % contradiction to data
    e	= inf;
  else
    % Any y=0 or y=1 cases must now be accompanied by appropriate
    % output=0 or output=1 values, so can be excluded.
    e	= -(Targets(~y0)'*log(y(~y0)) + (1-Targets(~y1))'*log(1-y(~y1)));
  end
  %
 case LIKELIHOOD.Poisson
  y	= exp(BASIS_Mu);
  e	= -sum(Targets.*BASIS_Mu - y);
  %
end

repo: REST


%==============================================================================
% This code is part of the Finite Element Method app for the Matlab-based toolbox
%  FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR/FAIRFEM
%==============================================================================
%
% classdef TriMesh3 < handle
%
% Finite Element Mesh based on triangular subdivision of rectangular mesh.
%
% Each Cell is divided into four triangles:
%
%  o---------o---------o---------o---------o
%  | \     / | \     / | \     / | \     / |
%  |   \ /   |   \ /   |   \ /   |   \ /   |
%  |    X    |    X    |    X    |    X    |
%  |  /  \   |  /  \   |  /  \   |  /  \   |
%  |/      \ |/      \ |/      \ |/      \ |
%  o---------o---------o---------o---------o
%  | \     / | \     / | \     / | \     / |
%  |   \ /   |   \ /   |   \ /   |   \ /   |
%  |    X    |    X    |    X    |    X    |
%  |  /  \   |  /  \   |  /  \   |  /  \   |
%  |/      \ |/      \ |/      \ |/      \ |
%  o---------o---------o---------o---------o
%  | \     / | \     / | \     / | \     / |
%  |   \ /   |   \ /   |   \ /   |   \ /   |
%  |    X    |    X    |    X    |    X    |
%  |  /  \   |  /  \   |  /  \   |  /  \   |
%  |/      \ |/      \ |/      \ |/      \ |
%  o---------o---------o---------o---------o
%
%  To construct an instance of this class type:
%
%  >> Mesh = TriMesh3(omega,m)
%
% Input:
% 	omega - description of spatial domain
%   m     - number of cells
%
% Properties:
%   xn     - node list
%   tri    - triangle list
%   dim    - space dimension
%   omega  - description of spatial domain
%   m      - number of cells
%   type   - type of partition
%   vol    - volume of triangles
%   nnodes - number of nodes
%   ntri   - number of triangles
%   dx1    - partial derivative operator
%   dx2    - partial derivative operator
%   GRAD   - gradient operator
%   P1     - projection operator for node 1
%   P2     - projection operator for node 2
%   P3     - projection operator for node 3
%   PC     - projection operator for Barycentrum
%   P      - prolongation operator
%   Pt     - restriction operator
%
%  Methods:
%   mfPu   - matrix free prolongation/restriction
%   mfPi   - matrix free edge projector
%   getP   - builds prolongation operator
%   tri2cc - averaging
%
%
% see also
% =========================================================================
classdef TriMesh3 < handle

    properties
        % ===================================
        % node list
        % ===================================
        xn
        % ===================================
        % triangle list
        % ===================================
        tri
        % ===================================
        % space dimension
        % ===================================
        dim  = 2;
        % ===================================
        % description of computational domain
        % ===================================
        omega
        % ===================================
        % number of cells
        % ===================================
        m
        % ===================================
        % type of partition
        % ===================================
        type = 3;
        % ===================================
        % function handle to myself
        % ===================================
        me = @TriMesh3
        % ===================================
        % volume of triangles
        % ===================================
        vol
        % ===================================
        % number of nodes
        % ===================================
        nnodes
        % ===================================
        % number of triangles in mesh
        % ===================================
        ntri
    end

    properties (Access = public, Dependent) % These will be created when first callend and stores persistently
        % ===================================
        % dx1  - partial derivative operator
        % ===================================
        dx1
        % ===================================
        % dx2  - Partial derivative operator
        % ===================================
        dx2
        % ===================================
        % GRAD - Gradient operator
        %
        %         | dx1 |
        %  GRAD = |     |
        %         | dx2 |
        %
        % ===================================
        GRAD
        % ===================================
        % B - Vector gradient operator
        %
        %         | GRAD  0    |
        %  B =    |            |
        %         | 0     GRAD |
        %
        % ===================================
        B
        % ===================================
        % P1 - Projection operator on Node 1
        % ===================================
        P1
        % ===================================
        % P2 - Projection operator on Node 2
        % ===================================
        P2
        % ===================================
        % P3 - Projection operator on Node 3
        % ===================================
        P3
        % ===================================
        % PC - Projection operator on Barycenter
        % ===================================
        PC
        % ===================================
        % P  - Prolongation operator
        % ===================================
        P
        % ===================================
        % Pt - Restriction operator
        % ===================================
        Pt
        % ===================================
        % Boundary indices
        % ===================================
        boundaryIdx
        % ===================================
        % Boundary projector
        % ===================================
        boundaryProj
        AvN
        mfdx1
        mfdx2
        mfGRAD
    end

    properties (Access = private)
        % These are where the dependent data is actually stored
        dx1_
        dx2_
        GRAD_
        B_
        P1_
        P2_
        P3_
        PC_
        P_
        Pt_
        mfdx1_
        mfdx2_
        mfGRAD_
        boundaryIdx_
        boundaryProj_
        AvN_
    end

    methods
        function this = TriMesh3(omega,m)
            if nargin==0,
                help(mfilename);
                this.runMinimalExample;
                return;
            end
            this.omega = omega;
            this.m     = m;

            this.xn = [ reshape(getNodalGrid(omega,m),[],2);...
               reshape(getCellCenteredGrid(omega,m),[],2)];

            % get indices of bottom left vertices
            nodes = reshape(1:prod(m+1),m+1);
            nodes = reshape(nodes(1:end-1,1:end-1),1,[]);
            % get indices of cell-centered points
            cc    = prod(m+1)+(1:prod(m));
            % specify triangles
            this.tri   = [...
                nodes;          nodes+1;        cc; ...
                nodes+1;        nodes+m(1)+2;   cc; ...
                nodes+m(1)+2;   nodes+m(1)+1;   cc; ...
                nodes+m(1)+1;   nodes;          cc; ...
                ];
            this.tri   = reshape(this.tri,3,[])';

            this.nnodes = prod(m+1)+prod(m);
            this.ntri  = size(this.tri,1);
            this.vol = prod((omega(2:2:end)-omega(1:2:end))./m)/4*ones(this.ntri,1);

        end

        function runMinimalExample(~)
            omega = [0 4 2 6]; m = [3 6];
            Mesh  = feval(mfilename,omega,m);
            figure(3); clf;
            set(gcf,'Name',sprintf('%s, m=[%d,%d]',mfilename,m));
            plotTriMesh(Mesh,Mesh.xn);
        end

        function x = mfPi(this,x,i)
            % =============================================================
            % function x = mfPi(this,x,i)
            %
            % matrix free edge projector
            % =============================================================
            switch i
                case 1
                    P = this.P1;
                case 2
                    P = this.P2;
                case 3
                    P = this.P3;
                case 'C'
                    P = this.PC;
            end
            if size(x,1) == this.ntri,
                % ajoint
                x = P'*x;
            else
                x = P * x;
            end
        end

        function x = tri2cc(this,x)
            % =============================================================
            % function x = tri2cc(x)
            %
            % averaging or adjoint
            % =============================================================
            if numel(x)==this.ntri,
                x = mean(reshape(x,4,[]),1);
            else
                x = reshape(x,1,[]);
                x = .25*repmat(x,[1 4]);
                x = x(:);
            end

        end

        % ===================================
        % x = getNodalGridData(this,x)
        %
        % gets values at triangle nodes that belong the underlying nodal
        % grid
        % ===================================
        function x = getNodalGridData(this,x)
            x = reshape(x,prod(this.nnodes),[]);
            x = x(1:prod(this.m+1),:);
            x = x(:);
        end

        % ===================================
        % x = getCCGridData(this,x)
        %
        % gets values at cell-centers that belong the underlying nodal
        % grid
        % ===================================
        function x = getCCGridData(this,x)
            x = reshape(x,prod(this.nnodes),[]);
            x = x(prod(this.m+1)+1:end,:);
            x = x(:);
        end

        function Pu = mfPuNodal(~,yn,m,flag)
            % =============================================================
            % function Pu = mfPuNodal(~,yn,m,flag)
            %
            % matrix free prolongation/restriction for nodal quantities
            % =============================================================
            switch flag
                case 'Pu' % coarse --> fine
                    yn = reshape(yn,[],2);
                    yc = reshape(yn(prod(m+1)+1:end,:),[m 2]);
                    yn = reshape(yn(1:prod(m+1),:),[m+1 2]);

                    Pyn = zeros([2*m+1 2]);
                    Pyc = zeros([2*m 2]);

                    % include existing nodes
                    Pyn(1:2:end,1:2:end,:) = yn;
                    % prolongate in x direction
                    Pyn(2:2:end-1,:,:) =  .5*(Pyn(1:2:end-2,:,:) + Pyn(3:2:end,:,:));
                    % prolongate in y direction
                    Pyn(1:2:end,2:2:end-1,:) =  .5*(Pyn(1:2:end,1:2:end-2,:) + Pyn(1:2:end,3:2:end,:));
                    % midpoints
                    Pyn(2:2:end,2:2:end,:)   = yc;

                    % compute cell-centers
                    Pyc(1:2:end-1,1:2:end-1,:) = .5*(yn(1:end-1,1:end-1,:) + yc);
                    Pyc(2:2:end  ,1:2:end-1,:) = .5*(yn(2:end  ,1:end-1,:) + yc);
                    Pyc(1:2:end-1,2:2:end  ,:) = .5*(yn(1:end-1,2:end  ,:) + yc);
                    Pyc(2:2:end  ,2:2:end  ,:) = .5*(yn(2:end  ,2:end  ,:) + yc);

                    Pu = [reshape(Pyn,[],2); reshape(Pyc,[],2)];
                    Pu = Pu(:);
                case 'PTu'
                    yn = reshape(yn,[],2);
                    yc = reshape(yn(prod(m+1)+1:end,:),[m 2]);
                    yn = reshape(yn(1:prod(m+1),:),[m+1 2]);

                    % include parent nodes
                    Pyn = yn(1:2:end,1:2:end,:);
                    % distribute in x direction
                    t = .5* yn(2:2:end-1,1:2:end,:);
                    Pyn(1:end-1,:,:)   = Pyn(1:end-1,:,:) + t;
                    Pyn(2:end  ,:,:)   = Pyn(2:end  ,:,:) + t;

                    % distribute in y direction
                    t = .5* yn(1:2:end,2:2:end-1,:);
                    Pyn(:,1:end-1,:) = Pyn(:,1:end-1,:) + t;
                    Pyn(:,2:end  ,:) = Pyn(:,2:end  ,:) + t;

                    % distribute cell-centers
                    t1 = .5*yc(1:2:end-1,1:2:end-1,:);
                    t2 = .5*yc(2:2:end  ,1:2:end-1,:);
                    t3 = .5*yc(1:2:end-1,2:2:end  ,:);
                    t4 = .5*yc(2:2:end  ,2:2:end  ,:);
                    Pyc = yn(2:2:end,2:2:end,:) + (t1+t2+t3+t4);
                    Pyn(1:end-1,1:end-1,:) = Pyn(1:end-1,1:end-1,:) + t1;
                    Pyn(2:end  ,1:end-1,:) = Pyn(2:end  ,1:end-1,:) + t2;
                    Pyn(1:end-1,2:end  ,:) = Pyn(1:end-1,2:end  ,:) + t3;
                    Pyn(2:end  ,2:end  ,:) = Pyn(2:end  ,2:end  ,:) + t4;

                    Pu = [reshape(Pyn,[],2); reshape(Pyc,[],2)];
                    Pu = Pu(:);
            end
        end

        function P = getPuNodal(~,m)
            % =============================================================
            % function P = getPuNodal(~,m)
            %
            % returns prolongation operator for input of cell-width m
            % =============================================================

            mf = 2*m;
            % indices of fine and coarse grid nodes
            indfn = reshape(1:prod(mf+1),mf+1);
            indfc = reshape(prod(mf+1)+(1:prod(mf)),mf);
            indcn = reshape(1:prod(m+1),m+1);
            indcc = reshape(prod(m+1)+(1:prod(m)),m);

            % allocate space
            I = []; J = []; W = [];

            % include existing nodes
            ii = indfn(1:2:end,1:2:end);
            jj = indcn(:);
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % include midpoint as node
            ii = indfn(2:2:end,2:2:end);
            jj = indcc(:);
            ww = ones(size(jj));
            I = [I; ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate in x direction
            ii   = indfn(2:2:end-1,1:2:end);
            jj   = [reshape(indcn(1:end-1,:),[],1); reshape(indcn(2:end,:),[],1)];
            ww = .5*ones(size(jj));
            I = [I; ii(:);ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate in y direction
            ii   = indfn(1:2:end,2:2:end-1);
            jj   = [reshape(indcn(:,1:end-1),[],1); reshape(indcn(:,2:end),[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate from bottom-left to top-right
            ii   = indfc(1:2:end,1:2:end);
            jj   = [reshape(indcn(1:end-1,1:end-1),[],1); reshape(indcc,[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            ii   = indfc(2:2:end,2:2:end);
            jj   = [reshape(indcn(2:end,2:end),[],1); reshape(indcc,[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            % prolongate from top-left to bottom-right
            ii   = indfc(2:2:end,1:2:end);
            jj   = [reshape(indcn(2:end,1:end-1),[],1); reshape(indcc,[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];
            ii   = indfc(1:2:end,2:2:end);
            jj   = [reshape(indcn(1:end-1,2:end),[],1); reshape(indcc,[],1)];
            ww   = .5*ones(size(jj));
            I = [I; ii(:); ii(:)]; J = [J; jj(:)]; W = [W; ww(:)];

            P = sparse(I,J,W,prod(mf+1)+prod(mf),prod(m+1)+prod(m));
        end

        % ========== get methods ========================================
        function dx1 = get.dx1(this)
            if isempty(this.dx1_),
                [this.dx1_, this.dx2_] = getGradientMatrixFEM(this,0);
            end
            dx1 = this.dx1_;
        end

        function dx2 = get.dx2(this)
            if isempty(this.dx2_),
                [this.dx1_, this.dx2_] = getGradientMatrixFEM(this,0);
            end
            dx2 = this.dx2_;
        end

        function GRAD = get.GRAD(this)
            if isempty(this.GRAD_),
                this.GRAD_ = [this.dx1;this.dx2];
            end
            GRAD = this.GRAD_;
        end
        function B = get.B(this)
            if isempty(this.B_),
                this.B_ = blkdiag(this.GRAD,this.GRAD);
            end
            B = this.B_;
        end

        function P1 = get.P1(this)
            if isempty(this.P1_),
                A = speye(this.nnodes);
                this.P1_ = A(this.tri(:,1),:);
            end
            P1 = this.P1_;
        end
        function P2 = get.P2(this)
            if isempty(this.P2_),
                A = speye(this.nnodes);
                this.P2_ = A(this.tri(:,2),:);
            end
            P2 = this.P2_;
        end

        function P3 = get.P3(this)
            if isempty(this.P3_),
                A = speye(this.nnodes);
                this.P3_ = A(this.tri(:,3),:);
            end
            P3 = this.P3_;
        end
        function PC = get.PC(this)
            if isempty(this.PC_),
                this.PC_ = (this.P1+this.P2+this.P3)/3;
            end
            PC = this.PC_;
        end
        function P = get.P(this)
            if isempty(this.P_),
                this.P_ = this.getPuNodal(this.m);
            end
            P = this.P_;
        end
        function Pt = get.Pt(this)
            if isempty(this.Pt_),
                this.Pt_ = this.getPuNodal(this.m/2);
            end
            Pt = this.Pt_;
        end

        function mfdx1 = get.mfdx1(this)
            if isempty(this.mfdx1_),
                [this.mfdx1_, this.mfdx2_] = getGradientMatrixFEM(this,1);
            end
            mfdx1 = this.mfdx1_;
        end

        function mfdx2 = get.mfdx2(this)
            if isempty(this.mfdx2_),
                [this.mfdx1_, this.mfdx2_] = getGradientMatrixFEM(this,1);
            end
            mfdx2 = this.mfdx2_;
        end

        function AvN = get.AvN(this)
            if isempty(this.AvN_),
                av  = @(n) spdiags(ones(n+1,1)*[0.5 0.5],[0,1],n,n+1);
                this.AvN_ = kron(speye(2), [speye(prod(this.m+1)); kron(av(this.m(2)),av(this.m(1)))]);
            end
            AvN = this.AvN_;
        end

        function mfGRAD = get.mfGRAD(this)
            if isempty(this.mfGRAD_),
                this.mfGRAD_ = getGradientMatrixFEM(this,1);
            end
            mfGRAD = this.mfGRAD_;
        end
        function idx = get.boundaryIdx(this)
            if isempty(this.boundaryIdx_),
                id = reshape(1:prod(this.m+1),this.m+1);
                idx = [reshape(id([1,end],:),[],1);reshape(id(2:end-1,[1,end]),[],1)] ;
                this.boundaryIdx_ = idx;
            end
            idx = this.boundaryIdx_;
        end

        function idx = get.boundaryProj(this)
            if isempty(this.boundaryProj_),
                idx = this.boundaryIdx;

                P = speye(prod(this.m+1));
                P = P(idx,:);
                P = [P,sparse(size(P,1),prod(this.m))];
                P = kron(speye(2),P);

                this.boundaryProj_ = P;
            end
            idx = this.boundaryProj_;
        end

        % ========== set methods ========================================
        function set.xn(this,xn)
            if numel(xn)~=(2*(prod(this.m+1)+prod(this.m))),
                error('Invalid number of nodes');
            end
            if isempty(this.xn),
                this.xn = xn;
            else

                this.xn = xn;
                % delete operators that are sensitive to nodes
                this.dx1_  = [];
                this.dx2_  = [];
                this.GRAD_ = [];
                this.B_    = [];

                this.vol = volTetraGrid(this,this.xn,'matrixFree',1);
            end
        end
    end
end


text:
function [X,T] = read_spm_file(spm_file,mat_file)
% Load an (continuous data) SPM file with filename specified in the spm_file
% input parameter, and put the data into X in matrix format ready to be used by
% the HMM. If mat_file is specified, X is saved in the file mat_file
% Importantly, if there are bad time points in the time series, it removes
% them and returns a vector of segment lengths (T) for the HMM to take
% into account that the time series is not seamlessly continuous anymore.
% In this case, although the hmmmar function can receive SPM files
% directly, it is strongly recommended to run this function for each
% subject before running the HMM, concatenate T across subjects,
% and then pass the mat_file's to the hmmmar function.
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

D = spm_eeg_load(spm_file);
X = D(:,:,:);
X = reshape(X,[D.nchannels,D.nsamples*D.ntrials]);

% select only good data
goodsamples = good_samples(D,[],[],[],0);
goodsamples = reshape(goodsamples,1,D.nsamples*D.ntrials);
X = X(:,goodsamples)';
options = struct();
T = cell2mat(getStateLifeTimes (goodsamples',length(goodsamples),options,1));

if min(T)<20
    warning('There are time segments with less than 20 time points')
end

if nargin > 1
   save(mat_file,'X','T')
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function calpak_test41 ( )

%*****************************************************************************80
%
%% CALPAK_TEST41 tests MONTH_TO_MONTH_NAME_HINDU_LUNAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST41\n' );
  fprintf ( 1, '  For the Hindu lunar calendar,\n' );
  fprintf ( 1, '  MONTH_TO_MONTH_NAME_HINDU_LUNAR names the months.\n' );
  fprintf ( 1, '\n' );

  y = 1;
  months = year_length_months_hindu_lunar ( y );

  for m = 1 : months
    month_name = month_to_month_name_hindu_lunar ( m );
    fprintf ( 1, '  %2d  %s\n', m, month_name );
  end

  return
end


text:
% load camera results (for setting active images, n_ima,...)
load camera_results;

% Projection settings:

wintx = 12; %8;
winty = 12; %8;
nx = 1024;
ny = 768;
dX = 32;
dY = 32;
dXoff=511.5;
dYoff=383.5;

proj_name = 'proj'; %input('Basename projector calibration images (without number nor suffix): ','s');
camera_name = 'cam';

xr_list = NaN*ones(1,n_ima);
yr_list = NaN*ones(1,n_ima);

ind_proc = ind_active;

DEBUG = 1;
%ind_ima_proj = [18];

for i = ind_proc,

   projector_marker;

end;

fprintf(1,'\nExtraction of the grid corners on the image\n');

recompute_corner = 0;

if recompute_corner & ~exist(['xproj_' num2str(ind_proc(1))]),
   if exist('projector_data'),
      load projector_data;
   else
      recompute_corner = 0;
      disp('WARNING: Cannot recompute corners. Data need to be extracted at least once');
   end;
end;

if ~recompute_corner,
   disp('Manual extraction mode');
else
   disp('Automatic recomputation of the corners');
end;

% extract the projector corners:

for kk = ind_proc,

   projector_ima_corners

end;

string_save = 'save projector_data ';

for kk = ind_active,
   string_save = [string_save ' xproj_' num2str(kk) ' x_proj_' num2str(kk)];
end;
eval(string_save);

return;

i = 18;

figure(4)
image(I_29);
colormap(gray(256));
hold on;
plot(x_29(1,:)+1,x_29(2,:)+1,'r+','markersize',13,'linewidth',3)
hold off;
axis image
axis off;

repo: ekfmonoslam


text:
function XvalResult = ncrossdecomp(Method,X,FacMin,FacMax,Segments,Cent,Show);

%NCROSSDECOMP crossvalidation of PARAFAC/Tucker/PCA
%
% See also:
% 'ncrossreg'
%
% This file performs cross-validation of decomposition models
% PARAFAC, PCA, and Tucker. The cross-validation is performed
% such that part of the data are set to missing, the model is
% fitted to the remaining data, and the residuals between fitted
% and true left-out elements is calculated. This is performed
% 'Segments' times such that all elements are left out once.
% The segments are chosen by taking every 'Segments' element of
% X(:), i.e. from the vectorized array. If X is of size 5 x 7,
% and three segemnts are chosen ('Segments' = 3), then in the
% first of three models, the model is fitted to the matrix
%
% |x 0 0 x 0 0 x|
% |0 x 0 0 x 0 0|
% |0 0 x 0 0 x 0|
% |x 0 0 x 0 0 x|
% |0 x 0 0 x 0 0|
%
% where x's indicate missing elements. After fitting the residuals
% in the locations of missing values are calculated. After fitting
% all three models, all residuals have been calculated.
%
% Note that the number of segments must be chosen such that no columns
% or rows contain only missing elements (the algorithm will check this).
% Using 'Segments' = 7, 9, or 13 will usually achieve that.
%
% I/O
% XvalResult = ncrossdecomp(Method,X,FacMin,FacMax,Segments,Cent,Show);
%
% INPUT
% Method   : 'parafac', 'tucker', 'pca', or 'nipals'
%            For Tucker only Tucker3 models with equal
%            number of components is currently available.
%            For PCA the least squares model is calculated.
%            Thus, offsets and parameters are calculated in
%            a least squares sense unlike the method NIPALS,
%            which calculates the PCA model using an ad hoc
%            approach for handling missing data (as in
%            standard chemometric software).
% X        : Multi-way array of data
% FacMin   : Lowest number of factors to use
% FacMax   : Highest number of factors (note that for Tucker only models
%            with the same number of components in each mode are
%            calculated currently
% Segments : The number of segments to use. Try many!
% Cent     : If set of one, the data are centered across samples,
%            i.e. ordinary centering. Note, however, that the centering
%            is not performed in a least squares sense but as preprocessing.
%            This is not optimal because the data have missing data because
%            of the way the elements are left out. This can give
%            significantly lower fit than reasonable if you have few samples
%            or use few segments. Alternatively, you can center the data
%            beforehand and perform cross-validation on the centered data
% Show     : If set to 0, no plot is given
%
% OUTPUT
% Structure XvalResult holding:
%           Fit: The fitted percentage of variation explaind (as a
%                function of component number)
%          Xval: The cross-validated percentage of variation explaind
%                (as a function of component number)
%   FittedModel: The fitted model (as a function of component number)
%     XvalModel: The cross-validated model (as a function of component number)

% $ Version 1.0301 $ Date 28. June 1999 $ Not compiled $
% $ Version 2.00 $ May 2001 $ Changed to array notation $ RB $ Not compiled $

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% 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., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

% uses NANSUM,

%GT
Perm  = [1:ndims(X)];
DimX  = size(X);
ord   = ndims(X);
delta = zeros(1,ndims(X));
for i = 1:ndims(X)
    c(i)     = isempty(intersect(factor(DimX(i)),factor(Segments)));
    index{i} = 1:DimX(i);
end
delta = zeros(1,ndims(X));
for i = 1:ndims(X)
    c(i)     = isempty(intersect(factor(DimX(i)),factor(Segments)));
    index{i} = 1:DimX(i);
end
P = find(~c);
if sum(c) < length(DimX) - 1
    for i = 1:sum(~c)
        while ~c(P(i)) & delta(P(i)) < (abs(DimX(P(i))- Segments) - 1)
            if ~rem(DimX(P(i)),2) & ~rem(Segments,2) & delta(P(i)) > 0
                Off = 2;
            else
                Off = 1;
            end
            delta(P(i)) = delta(P(i)) + Off;
            c(P(i)) = isempty(intersect(factor(DimX(P(i)) + delta(P(i))),factor(Segments)));
        end
        if sum(c) == length(DimX) - 1
            return
        end
    end
    if sum(~c) > 1
        error('The chosen segmentation leads to tubes of only missing values')
    end
end
if sum(c) == length(DimX) - 1
    [c,Perm] = sort(~c);
end
[nil,PermI] = sort(Perm);
X           = reshape(X,DimX(1),prod(DimX(2:end)));
%GT end

[I,J] = size(X);

if exist('Show')~=1
   Show = 1;
end

if lower(Method(1:3)=='tuc')
   if length(FacMin)==1
      FacMin = ones(1,ord)*FacMin;
   elseif length(FacMin)~=ord
      error('Error in FacMin: When fitting Tucker models, the number of factors should be given for each mode')
   end
   if length(FacMax)==1
      FacMax = ones(1,ord)*FacMax;
   elseif length(FacMax)~=ord
      error('Error in FacMax: When fitting Tucker models, the number of factors should be given for each mode')
   end
end

%RB
% Check if the selected segmentation works (does not produce rows/columns of only missing)
%out = ones(I,J);
%out(1:Segments:end)=NaN;
%out(find(isnan(X))) = NaN;
%if any(sum(isnan(out))==I)
%   error(' The chosen segmentation leads to columns of only missing elements')
%elseif any(sum(isnan(out'))==J)
%   error(' The chosen segmentation leads to rows of only missing elements')
%end
%RB end

if lower(Method(1:3)~='tuc')

   XvalResult.Fit = zeros(FacMax,2)*NaN;
   XvalResult.Xval = zeros(FacMax,2)*NaN;
   for f = 1:FacMin-1
      XvalResult.XvalModel{f} =  'Not fitted';
      XvalResult.FittedModel{f} = 'Not fitted';
   end

   for f = FacMin:FacMax

      % Fitted model
      disp([' Total model - Comp. ',num2str(f),'/',num2str(FacMax)])
      [M,Mean,Param] = decomp(Method,X,DimX,f,1,Segments,Cent,I,J);
      Model = M + ones(I,1)*Mean';
      id    = find(~isnan(X));
      OffsetCorrectedData = X - ones(I,1)*Mean';
      XvalResult.Fit(f,:) = [100*(1 - sum( (X(id) - Model(id)).^2)/sum(OffsetCorrectedData(id).^2)) f];
      XvalResult.FittedModel{f} = Model;

      % Xvalidated Model of data
      ModelXval = zeros(I,J)*NaN;
      for s = 1:Segments
         disp([' Segment ',num2str(s),'/',num2str(Segments),' - Comp. ',num2str(f),'/',num2str(FacMax)])
         %GT
         Xnow    = permute(zeros(DimX),Perm);
         dimsadd = size(Xnow);
         for j = 1:length(DimX)
            dimsadd(j) = delta(Perm(j));
            Xnow       = cat(j,Xnow,zeros(dimsadd));
            index2{j}  = 1:DimX(Perm(j));
            dimsadd    = size(Xnow);
         end
         Xnow(s:Segments:end)     = NaN;
         Xnow                     = reshape(permute(Xnow(index2{:}),PermI),DimX(1),prod(DimX(2:end)));
         Pos                      = isnan(Xnow);
         [M,Mean]                 = decomp(Method,Xnow + X,DimX,f,s,Segments,Cent,I,J,Param);
         model                    = M + ones(I,1)*Mean';
         ModelXval(Pos)           = model(Pos);
         %GT end
      end

      XvalResult.Xval(f,:) = [100*(1 - sum( (X(id) - ModelXval(id)).^2)/sum(OffsetCorrectedData(id).^2)) f];
      XvalResult.XvalModel{f} =  ModelXval;
   end

else % Do Tucker model

   %GT
   if length(FacMin) ~= ord
      FacMin = ones(1,ord)*min(FacMin);
   end
   if length(FacMax) ~= ord
      FacMax = ones(1,ord)*min(FacMax);
   end
   for i=1:length(FacMin)
      ind{i} = [FacMin(i):FacMax(i)]';
      ind2{i} = ones(length(ind{i}),1);
   end
   NCombs = cellfun('length',ind);
   for i=1:length(FacMin)
      o = [1:i-1,i+1:length(FacMin)];
      t = ipermute(ind{i}(:,ind2{o}),[i,o]);
      possibleCombs(1:prod(NCombs),i) = t(:);
   end
   FeasCombs = sort(possibleCombs');
   f2        = prod(FeasCombs(1:size(FeasCombs,1)-1,:))<FeasCombs(end,:);
   possibleCombs(f2,:) = [];
   possibleCombs(:,end+1) = find(~f2(:));
   %GTend
   %RB
   % Find all
   %PossibleNumber = [min(FacMin):max(FacMax)]'*ones(1,ord);
   %possibleCombs = unique(nchoosek(PossibleNumber(:),ord),'rows');
   %remove useless
   %f2 = [];
   %for f1 = 1:size(possibleCombs,1)
   %   if (prod(possibleCombs(f1,:))/max(possibleCombs(f1,:)))<max(possibleCombs(f1,:)) % Check that the largest mode is larger than the product of the other
   %      f2 = [f2;f1];
   %   elseif any(possibleCombs(f1,:)>FacMax)  % Chk the model is desired,
   %      f2 = [f2;f1];
   %   end
   %end
   %possibleCombs(f2,:)=[];
   %[f1,f2]=sort(sum(possibleCombs'));
   %possibleCombs = [possibleCombs(f2,:) f1'];
   %RBend

   XvalResult.Fit = zeros(size(possibleCombs,1),ord+1)*NaN;
   XvalResult.Xval = zeros(size(possibleCombs,1),ord+1)*NaN;
   for f = 1:size(possibleCombs,1)
      XvalResult.XvalModel{f} =  'Not fitted';
      XvalResult.FittedModel{f} = 'Not fitted';
   end

   for f1 = 1:size(possibleCombs,1)

      % Fitted model
      %GT
      disp([' Total model - Comp. ',num2str(possibleCombs(f1,:)),'/',num2str(FacMax)])
      [M,Mean,Param] = decomp(Method,X,DimX,possibleCombs(f1,:),1,Segments,Cent,I,J);
      %GTend
      %RB
      %disp([' Total model - Comp. ',num2str(possibleCombs(f1,1:end-1)),'/',num2str(FacMax)])
      %[M,Mean,Param] = decomp(Method,X,DimX,possibleCombs(f1,1:end-1),1,Segments,Cent,I,J);
      %RBend
      Model = M + ones(I,1)*Mean';
      id    = find(~isnan(X));
      OffsetCorrectedData = X - ones(I,1)*Mean';
      XvalResult.Fit(f1,:) = [100*(1 - sum( (X(id) - Model(id)).^2)/sum(OffsetCorrectedData(id).^2)) possibleCombs(f1,1:end-1)];
      XvalResult.FittedModel{f1} = Model;

      % Xvalidated Model of data
      ModelXval = zeros(I,J)*NaN;
      for s = 1:Segments
         disp([' Segment ',num2str(s),'/',num2str(Segments),' - Comp. ',num2str(possibleCombs(f1,1:end-1)),'/',num2str(FacMax)])
         %GT
         Xnow    = permute(zeros(DimX),Perm);
         dimsadd = DimX;
         for j = 1:length(DimX)
             dimsadd(j) = delta(j);
             Xnow       = cat(j,Xnow,zeros(dimsadd));
             index2{j}  = 1:DimX(j);
             dimsadd(j) = dimsadd(j) + DimX(j);
         end
         Xnow(s:Segments:end)     = NaN;
         Xnow                     = reshape(permute(Xnow(index2{:}),PermI),DimX(1),prod(DimX(2:end)));
         Pos                      = isnan(Xnow);
         [M,Mean]                 = decomp(Method,Xnow + X,DimX,possibleCombs(f1,1:end-1),s,Segments,Cent,I,J,Param);
         model                    = M + ones(I,1)*Mean';
         ModelXval(Pos)           = model(Pos);
         %GT end
      end

      XvalResult.Xval(f1,:) = [100*(1 - sum( (X(id) - ModelXval(id)).^2)/sum(OffsetCorrectedData(id).^2)) possibleCombs(f1,1:end-1)];
      XvalResult.XvalModel{f1} =  ModelXval;
   end

end

if Show&FacMin-FacMax~=0
   if Method(1:3) == 'pca'
      Nam = 'PCA';
   elseif Method(1:3) == 'tuc'
      Nam = 'Tucker';
   elseif Method(1:3) == 'par'
      Nam = 'PARAFAC';
   elseif Method(1:3) == 'nip'
      Nam = 'NIPALS';
   end

   figure
   save jjj
   if lower(Method(1:3))~='tuc'
      bar(FacMin:FacMax,[XvalResult.Fit(FacMin:FacMax,1) XvalResult.Xval(FacMin:FacMax,1)],.76,'grouped')
   else
      % extract the ones with lowest Xval fit (for each # total comp) for plotting
      fx = [];
      f5 =[];
      for f1 = 1:max(possibleCombs(:,end))
         f2 = find(possibleCombs(:,end)==f1);
         if length(f2)
            [f3,f4] = max(XvalResult.Xval(f2));
            f5 = [f5,f2(f4)];
         end
      end
      fx = [possibleCombs(f5,end) XvalResult.Fit(f5,1) XvalResult.Xval(f5,1)];
      bar(fx(:,1),fx(:,2:3),.76,'grouped')
      for f1 = 1:size(fx,1)
         f6=text(fx(f1,1),95,['[',num2str(possibleCombs(f5(f1),1:end-1)),']']);
         set(f6,'Rotation',270)
      end
   end

   g=get(gca,'YLim');
   set(gca,'YLim',[max(-20,g(1)) 100])
   legend('Fitted','Xvalidated',0)
   titl = ['Xvalidation results (',Nam,')'];
   if Cent
      titl = [titl ,' - centering'];
   else
      titl = [titl ,' - no centering'];
   end
   title(titl,'FontWeight','Bold')
   xlabel('Total number of components')
   ylabel('Percent variance explained')
end

function [M,Mean,parameters] = decomp(Method,X,DimX,f,s,Segments,Cent,I,J,parameters);

Conv = 0;
it = 0;
maxit = 500;
% Initialize
if Cent
   Mean = nanmean(X)';
else
   Mean = zeros(J,1);
end

if lower(Method(1:3)) == 'par'
   Xc = reshape(X- ones(I,1)*Mean',DimX);
   if exist('parameters')==1
      fact = parafac(Xc,f,[1e-5 10 0 0 NaN maxit],[],parameters.fact);
   else
      fact = parafac(Xc,f,[1e-5 10 0 0 NaN maxit]);
   end
   M = reshape(nmodel(fact),DimX(1),prod(DimX(2:end)));
   parameters.fact=fact;
elseif lower(Method(1:3)) == 'tuc'
   Xc = reshape(X- ones(I,1)*Mean',DimX);
   if exist('parameters')==1
      [fact,G] = tucker(Xc,f,[1e-2 0 0 0 NaN maxit],[],[],parameters.fact,parameters.G);
   else
      [fact,G] = tucker(Xc,f,[1e-2 0 0 0 NaN maxit]);
   end
   parameters.fact=fact;
   parameters.G=G;
   M = reshape(nmodel(fact,G),DimX(1),prod(DimX(2:end)))      ;
elseif lower(Method) == 'pca'|lower(Method) == 'nip'
   Xc = reshape(X- ones(I,1)*Mean',DimX(1),prod(DimX(2:end)));
   [t,p] = pcanipals(X- ones(I,1)*Mean',f,0);
   parameters.t=t;
   parameters.p=p;
   M = t*p';
else
   error(' Name of method not recognized')
end
Fit = X - M - ones(I,1)*Mean';
Fit = sum(Fit(find(~isnan(X))).^2);

% Iterate
while ~Conv
   it     = it+1;
   FitOld = Fit;

   % Fit multilinear part
   Xcent = X - ones(I,1)*Mean';

   if Method(1:3) == 'par'
      fact = parafac(reshape(Xcent,DimX),f,[1e-2 0 0 0 NaN maxit],[],fact);
      M = reshape(nmodel(fact),DimX(1),prod(DimX(2:end)));
   elseif Method(1:3) == 'tuc'
      [fact,G] = tucker(reshape(Xcent,DimX),f,[1e-2 0 0 0 NaN maxit],[0 0 0],zeros(size(G)),fact,G);
      M = reshape(nmodel(fact,G),DimX(1),prod(DimX(2:end)));
   elseif Method == 'pca'
      [t,p] = pcals(Xcent,f,0,t,p,0);
      M = t*p';
   elseif Method == 'nip'
      [t,p] = pcanipals(Xcent,f,0);
      M = t*p';
   end

   % Find offsets
   if Cent
      x = X;
      mm=M+ones(I,1)*Mean';
      x(find(isnan(X)))=mm(find(isnan(X)));
      Mean = mean(x)';
   end

   %Find fit
   Fit = X - M - ones(I,1)*Mean';
   Fit = sum(Fit(find(~isnan(X))).^2);
   if abs(Fit-FitOld)/FitOld<1e-8 | it > 1500
      Conv = 1;
   end

end
disp([' Fit ',num2str(Fit),' using ',num2str(it),' it.'])

function [t,p] = pcals(X,F,cent,t,p,show);

%  LEAST SQUARES PCA WITH MISSING ELEMENTS
%  20-6-1999
%
%  Calculates a least squares PCA model. Missing elements
%  are denoted NaN. The solution is NOT nested, so one has
%  to calculate a new model for each number of components.

ShowMeFitEvery = 20;
MaxIterations  = 5;
[I,J]=size(X);
Xorig      = X;
Miss       = find(isnan(X));
NotMiss    = find(~isnan(X));
m          = t*p';
X(Miss)    = m(Miss);
ssX    = sum(X(NotMiss).^2);

Fit    = 3;
OldFit = 6;
it     = 0;

while abs(Fit-OldFit)/OldFit>1e-3 & it < MaxIterations;
   it      = it +1;
   OldFit  = Fit;

   [t,s,p] = svds(X,F);
   t       = t*s;

   Model   = t*p';
   X(Miss) = Model(Miss);
   Fit     = sum(sum( (Xorig(NotMiss) - Model(NotMiss)).^2));

   if ~rem(it,ShowMeFitEvery)&show
      disp(['    Fit after ',num2str(it),' it. :',num2str(RelFit),'%'])
   end

end

function [t,p,Mean] = pcanipals(X,F,cent);

% NIPALS-PCA WITH MISSING ELEMENTS
% cent: One if centering is to be included, else zero

[I,J]=size(X);
rand('state',sum(100*clock))

Xorig      = X;
Miss       = isnan(X);
NotMiss    = ~isnan(X);
ssX    = sum(X(find(NotMiss)).^2);
Mean   = zeros(1,J);
if cent
   Mean    = nanmean(X);
end
X      = X - ones(I,1)*Mean;

t=[];
p=[];

for f=1:F
   Fit    = 3;
   OldFit = 6;
   it     = 0;
   T      = rand(I,1);
   P      = rand(J,1);
   Fit    = 2;
   FitOld = 3;

   while abs(Fit-FitOld)/FitOld>1e-7 & it < 100;
      FitOld  = Fit;
      it      = it +1;

      for j = 1:J
         id=find(NotMiss(:,j));
         if length(id)==0
            id,end
         P(j) = T(id)'*X(id,j)/(T(id)'*T(id));
      end
      P = P/norm(P);

      for i = 1:I
         id=find(NotMiss(i,:));
         T(i) = P(id)'*X(i,id)'/(P(id)'*P(id));
      end

      Fit = X-T*P';
      Fit = sum(Fit(find(NotMiss)).^2);
   end
   t = [t T];
   p = [p P];
   X = X - T*P';

end

function Xc = nanmean(X)

if isempty(X)
   Xc = NaN;
   return
end

i = isnan(X);
j = find(i);
i = sum(i);
X(j) = 0;
Num = size(X,1)-i;
Xc = sum(X);
i = find(Num);
Xc(i) = Xc(i)./Num(i);
Xc(find(~Num))=NaN;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
clear;  clc;
% P2.5
n = 0:10;
x = 10*exp(-j*0.4*pi*n);
[xe, xo, m] = evenoddcomplex(x,n);
subplot(3,2,1); stem(n,real(x),'ko'); title('Real Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(3,2,3); stem(m,real(xe),'ko'); title('Real Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(3,2,5); stem(m,real(xo),'ko'); title('Real Odd Part')
xlabel('n'); ylabel('xo(n)');
subplot(3,2,2); stem(n,imag(x),'ko'); title('Imaginary Original Sequence')
xlabel('n'); ylabel('x(n)');
subplot(3,2,4); stem(m,imag(xe),'ko'); title('Imaginary Even Part')
xlabel('n'); ylabel('xe(n)');
subplot(3,2,6); stem(m,imag(xo),'ko'); title('Imaginary Odd Part')
xlabel('n'); ylabel('xo(n)');

author: Sable

repo: mcbench-benchmarks


text:
%SCATTERD Display scatterplot
%
%   H = SCATTERD(A)
%   H = SCATTERD(A,DIM,S,CMAP,FONTSIZE,'label','both','legend','gridded')
%
% INPUT
%   A     Dataset or matrix
%   DIM   Number of dimensions: 1,2 or 3 (optional; default: 2)
%   S     String specifying the colors and markers (optional)
%   CMAP  Matrix with a color map (optional)
%
% OUTPUT
%   H     Vector of handles
%
% DESCRIPTION
% SCATTERD(A) displays a 2D scatterplot of the first two features of the
% dataset A. If the number of dimensions DIM is specified (1..3), it plots
% the first D features in a D-dimensional plot (D<4). If the plot string S
% is provided, e.g. S = 'w+', all points are plotted accordingly. If given,
% different plot strings are used for different classes. See PLOT for the
% specification of plot strings.
%
% If CMAP is specified, the color of the object symbols is determined by
% CMAP indexed by the object labels. A colormap has the size [C x 3], where
% C is the number of classes. The three components of CMAP(I,:) determine
% the red, green and blue components of the color. For instance:
%
%   MAP = HSV; [M,K] = SIZE(A); LABELS = CEIL(64*[1:M]'/M);
%   A = DATASET(A,LABELS); SCATTERD(A,'.','COLORMAP',MAP);
%
% This may be used for tracking ordered objects.
%
% FONTSIZE may be a vector with three elements: fontsize, markersize and
% size of the label font in case of a label plot.
%
% Various other options are:
%   'LABEL'  : plot labels instead of symbols
%   'BOTH'   : plot labels next to each sample
%   'LEGEND' : place a legend in the figure
%   'GRIDDED': make a grid of 2D scatterplots of each pair of features
%
% All the parameters, except for the dataset A can be specified in any
% order or can be left out.
%
% Classifiers can be plot in the scatterplot by PLOTC.
% Note that PLOTC does not work for 1D and 3D scatterplots.
%
% EXAMPLES
% See PREX_CONFMAT, PREX_DENSITY, PREX_PLOTC, PREX_MCPLOT.
%
% SEE ALSO
% DATASETS, COLORMAP, PLOT, PLOTC

% Copyright: D. de Ridder, R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: scatterd.m,v 1.5 2009/01/04 21:43:17 duin Exp $

% REVISIONS
% DR1 - Dick, 05-10-2004
% Added plotting of unlabeled data as 'k.'.

function handle = scatterd(a,p1,p2,p3,p4,p5,p6,p7,p8)

	prtrace(mfilename);

	% Defaults
	d = min(size(a,2),2);		% Dimensionality of plot
	s = []; 					% Plot symbol(s)
	map = []; 					% Color map
	plotlab = 0;				% Put text labels instead of or next to samples
	plotsym = 1;				% Plot symbols for samples
	plotlegend = 0; 		    % Plot legend
	gridded = 0;				% Make a gridded plot
	gridrun = 0;				% Inner loop in a gridded plot?
	font_size = [];
	mark_size = [];
	lab_size  = [];
	hold_axis = ishold; % A flag to check if 'hold on' is set for the current axis

	a = prdataset(a); 		% Allow for a non-dataset data

	if (nargin < 9), par{8} = []; else, par{8} = p8; end
	if (nargin < 8), par{7} = []; else, par{7} = p7; end
	if (nargin < 7), par{6} = []; else, par{6} = p6; end
	if (nargin < 6), par{5} = []; else, par{5} = p5; end
	if (nargin < 5), par{4} = []; else, par{4} = p4; end
	if (nargin < 4), par{3} = []; else, par{3} = p3; end
	if (nargin < 3), par{2} = []; else, par{2} = p2; end
	if (nargin < 2), par{1} = []; else, par{1} = p1; end

	% Set up default values.
	for i = 1:5
		if (~isempty(par{i}))
			if (length(par{i}) == 1) & par{i} < 5 & (~ischar(par{i})) 	% Dimensionality
				d   = par{i}; par{i} = 2; 		% Necessary for gridded: D needs to be 2.
			elseif ((size(par{i},1) > 1) & (size(par{i},2)==3) & (~ischar(par{i}))) 	% Color map
				map = par{i};
			elseif (strcmp(par{i},'label'))
				plotlab = 1; plotsym = 0;
			elseif (strcmp(par{i},'both'))
				plotlab = 1; plotsym = 1;
			elseif (strcmp(par{i},'legend'))
				plotlegend = 1;
			elseif (strcmp(par{i},'gridded'))
				gridded = 1;
				par{i} = 'gridrun'; 					% Necessary for gridded: otherwise an infinite recursion :)
			elseif (strcmp(par{i},'gridrun'))
				gridrun = 1;
			elseif ~isstr(par{i}) & length(par{i}) <= 3 & par{i}(1) >= 5
				font_size = par{i}(1);
				if length(par{i}) >= 2
					mark_size = par{i}(2);
				end
				if length(par{i}) == 3
					lab_size = par{i}(3);
				end
			else
				s	= par{i};
			end
		end
	end

	if (gridrun)
		if isempty(font_size), font_size = 10; end
		if isempty(mark_size), mark_size =  5; end
		if isempty(lab_size),   lab_size =  8; end
	elseif (~hold_axis)
		%clf;
		cla;
		if isempty(font_size), font_size = 16; end
		if isempty(mark_size), mark_size =  7; end
		if isempty(lab_size),   lab_size = 14; end
	else
		if isempty(font_size), font_size = get(gca,'fontsize'); end
		if isempty(mark_size), mark_size = font_size/2; end
		if isempty(lab_size),   lab_size = font_size-2; end
	end

	feats = getfeatlab(a,'string');

	if isempty(feats)
		for i=1:size(a,2)
			feats = strvcat(feats,sprintf('Feature %d',i));
		end
	else
		if size(feats,2) == 1
			feats = [repmat('Feature ',size(feats,1),1) feats];
		end
	end

	if (gridded)
		clf;
		gs = size(a,2);
		for i = 1:gs
			for j = 1:gs
				subplot(gs,gs,(i-1)*gs+j);
				gridrun = 1;
				h = feval(mfilename,a(:,[j i]),par{1},par{2},par{3},par{4},par{5},par{6},par{7});
				gridrun = 0;
				if (i~=gs), xlabel(''); end
				if (j~=1),  ylabel(''); end
			end
			subplot(111); 		% Takes care of clf for the next scatterplot.
		end
		return;
	end

	if (isa(a,'prdataset')) & (~isempty(getlablist(a)))
        [m,k,c] = getsize(a);
        cs_a = classsizes(a);
        classes_present = cs_a > 0;
		lab = getnlab(a);
		lablist = getlablist(a,'string');
		ident = getident(a); %DXD: needed for prcursor
		dataset_name = getname(a);
		a = double(a);
	else
		[m,k] = size(a);
		lab = ones(m,1);
		ident = (1:m)'; %DXD: needed for prcursor
		dataset_name = [];
		c = 1;
		a = double(a);
	end

	% Character string defining the plotting setup in terms of symbols and colors.
	if (isempty(s))
		vers = version;
		if (str2num(vers(1)) < 5)
			col = 'brmw';
			sym = ['+*xo.']';
			i   = [1:20];
			ss  = [col(i-floor((i-1)/4)*4)' sym(i-floor((i-1)/5)*5)];
		else
			col = 'brmk';
			sym = ['+*oxsdv^<>p']';
			i   = [1:44];
			ss  = [col(i-floor((i-1)/4)*4)' sym(i-floor((i-1)/11)*11)];
		end
    ss  = ['k.'; ss];  								% DR1 - Add symbol for "unlabeled" data.
		[ms,ns] = size(ss);
		if ms == 1, ss = setstr(ones(m,1)*ss); end
	else
		if size(s,1) == 1
			ss = repmat(s,c,1); s = [];
		else
			ss = s; s = [];
		end
    ss  = char('k.', ss);  								% DR1 - Add symbol for "unlabeled" data.
	                                                %DXD - changed [.] into char(.)
	end

	% Define some 'space' OY to be added around the data plotted in symbols.
	oy = zeros(1,3);
	if (plotsym)
		oy = 0.02*(max(a)-min(a));
	else
		s = 'w.'; 				% Plot white spot instead of symbols.
	end
	oy(2) = 0;

	% Make a plot.
	lhandle = []; thandle = [];

    classes_present_idx = find(classes_present);

    [dummy cs_a_idx] = sort(cs_a(classes_present_idx),'descend');

  % Also plot label "0" (unlabeled).
	for i = [0 classes_present_idx(cs_a_idx) ]
		J = find(lab==i);
		if (isempty(s)), symbol = ss(i+1,:); else, symbol = s; end
		if ((d == 1) & ~isempty(J))
			h = plot(a(J,1),zeros(length(J),1),symbol);
			hold on;
			set(h,'markersize',mark_size);
			lhandle = [lhandle h];
			if (plotlab)
				for j = 1:length(J)
					h = text(a(J(j),1)+oy(1),oy(2),lablist(lab(J(j)),:));
					set(h,'fontsize',lab_size);
					thandle = [thandle h]; if (~isempty(map)), set (h, 'color', map(i+1,:)); end
				end
			end
		elseif ((d == 2) & ~isempty(J))
			h = plot(a(J,1),a(J,2),symbol);
			hold on;
			set(h,'markersize',mark_size);
			lhandle = [lhandle h];
			if (plotlab)
				for j = 1:length(J)
					h = text(a(J(j),1)+oy(1),a(J(j),2)+oy(2),lablist(lab(J(j)),:));
					set(h,'fontsize',lab_size);
					thandle = [thandle h]; if (~isempty(map)), set (h, 'color', map(i+1,:)); end
				end
			end
		elseif (~isempty(J))
			h = plot3(a(J,1),a(J,2),a(J,3),symbol);
			hold on;
			set(h,'markersize',mark_size);
			lhandle = [lhandle h];
			if (plotlab)
				for j = 1:length(J)
					h = text(a(J(j),1)+oy(1),a(J(j),2)+oy(2),a(J(j),3)+oy(3),lablist(lab(J(j)),:));
					set(h,'fontsize',lab_size);
					thandle = [thandle h]; if (~isempty(map)), set (h, 'color', map(i+1,:)); end
				end
			end
		end
		%DXD: store the object identifiers in the userdata such that you
		%can retrieve them by prcursor:
		if ~isempty(J)
			ud = get(h,'UserData');
			ud.ident = ident(J);
			set(h,'UserData',ud);
		end
    end

	if (plotsym)
		if (~isempty(map))
			for i = 0:c, set (lhandle(i+1), 'color', map(i+1,:)); end
		end
		if (plotlegend),
			[ht, hl] = legend(lhandle(:), [ num2str(cs_a(classes_present_idx(cs_a_idx))') repmat(' ', sum(classes_present),1) lablist(classes_present_idx(cs_a_idx),:) ]);
			hl = hl(:)';
			%set(hl(1:2*c),'markersize',mark_size);
			lhandle = [lhandle hl];
			thandle = [thandle ht(:)'];
		end
	end

	% !%_%*!_% Matlab

	set(gca,'fontsize',font_size);

	if (~hold_axis)
		dd = (max(a) - min(a))*0.05;		% offset, avoiding points on plot box.
		J = find(dd==0);
		dd(J) = dd(J)+1;

%feats
		if (d == 1),
			axis ([min(a(:,1))-dd(1) max(a(:,1))+dd(1) -0.5 0.5]);
			hx = xlabel(feats(1,:));
			thandle = [thandle hx];
		elseif (d == 2),
			axis ([min(a(:,1))-dd(1) max(a(:,1))+dd(1) min(a(:,2))-dd(2) max(a(:,2))+dd(2)]);
			hx = xlabel(feats(1,:)); hy = ylabel(feats(2,:));
			thandle = [thandle(:)' hx hy];
			view(2);
		else		% D = 3
			axis ([min(a(:,1))-dd(1) max(a(:,1))+dd(1) min(a(:,2))-dd(2) max(a(:,2))+dd(2) min(a(:,3))-dd(3) max(a(:,3))+dd(3)]);
			hx = xlabel(feats(1,:)); hy = ylabel(feats(2,:)); hz = zlabel(feats(3,:));
			thandle = [thandle hx hy hz];
			view(3);
		end
	end

	if (~gridrun) & (length(get(gcf,'children')) < 3 & any(get(gca,'position')>0.80))
		set(gca,'position',[0.13 0.13 0.79 0.78]); % make axis labels visible
	end

	if (~gridrun) & (~isempty(dataset_name))
		title(dataset_name);
	end

	hold off;

	if (nargout > 0)
		handle = [lhandle thandle];
	end

return;


text:
function displayIVPinfo(u, isIVP, varargin)
%DISPLAYINFOIVP    Utility routine for displaying IVP solving
%
% Calling sequence:
%   DISPLAYINFO(VARARGIN)
% where
%   U:          The solution computed.
%   ISIVP:      Equal to 1 if we've just solved an IVP, 0 otherwise (FVP).
%   VARARGIN:   Various useful information passed from the Newton solver
%               to this method.
%
% See also: chebop/displayInfo.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developers note: This method was introduced to allow the IVP solver to use
% similar syntax as the BVP solver for displaying useful information about the
% solution process.
%
% At the moment, we only display information at the end of solving IVPs, rather
% than during the Newton iteration which is the case for BVPs. Thus, this method
% only gets called at the end of the solveivp method, which is why we don't need
% the 'mode' argument like we do in the BVP case.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Are we dealing with a CHEBFUN or a CHEBMATRIX?
if (isa(u, 'chebfun'))
    len = length(u);
else
    len = max(cellfun(@length, u.blocks));
end

% Did we just solve an IVP or an FVP?
if ( isIVP == 1 )
    fprintf('Initial value problem detected.\n');
else
    fprintf('Final value problem detected.\n');
end

% Show information about the number of pieces and the total length of the
% solution:
fprintf('Number of pieces of the solution: %i.\n', length(u.domain) - 1);
fprintf('Total length of solution: %i.\n', len);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [flag] = IsRow(v,n)
% Default setting
flag = false;

% Get vector dimensions
[a,b] = size(v);

% Check it is a column
if a == 1 && nargin == 1
    flag = true;
    return;
end

% Otherwise check its length against n
if b == n           % Is not same length
    flag = true;
    return
end
end

repo: OpenMAS


text:
function kern = tensorKernSlash(kern, index)

% TENSORKERNSLASH Tensor kernel created by removing ith component.

% KERN

kern.nParams = kern.nParams - kern.comp{index}.nParams;
kern.comp = kern.comp([1:index-1 index+1:end]);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function dz = rollingCartPoleDynamics(~,z,P)
%DZ = ROLLINGCARTPOLEDYNAMICS(T,Z,P)
%
%FUNCTION:  This function computes the dynamics of a rolling
%    cart-pole, and is designed to be called from ode45.
%
%INPUTS:
%    t = time. Dummy input for ode45. Not used.
%    z = [4xn] matrix of states.
%    P = struct of parameters
%OUTPUTS:
%    dz = [4xn] matrix of state derivatives
%
%NOTES:
%    This file was automatically generated by writeRollingCartPoleDynamics

m1 = P.m1; %disk mass
m2 = P.m2; %bob mass
I = P.I; %disk moment of inertia
g  = P.g ; %gravity
l = P.l; %pendulum length
r = P.r; %disk radius

phi = z(1,:); %link one absolute angle
th = z(2,:); %link one angular rate
dphi = z(3,:); %link two absolute angle
dth = z(4,:); %link two angular rate

f1 = - m2.*r.*sin(th).*(dth.^2.*l + g.*cos(th)) - dth.^2.*l.*m1.*r.*sin(th);
f2 = -g.*m2.*sin(th);

M11 = - I - m1.*r.^2 - m2.*r.^2.*sin(th).^2;
M12 = l.*m1.*r.*cos(th);
M21 = m2.*r.*cos(th);
M22 = -l.*m2;

D = M11.*M22 - M12.*M21;

ddphi = (f2.*M12 - f1.*M22)./D;
ddth = -(f2.*M11 - f1.*M21)./D;

dz = [...
    dphi; %derivative of link one absolute angle
    dth; %derivative of link one angular rate
    ddphi; %derivative of link two absolute angle
    ddth; %derivative of link two angular rate
];
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function nfgCheckBMVolume(phantomDir)
%Verify that the BlueMatter volume calculation is correct
%
%   nfgCheckBMVolume(phantomDir)
%
%
% AUTHORS:
%   2009.08.05 : AJS wrote it
%
% NOTES:

% Directories
volcheckDir = nfgGetName('volcheckDir',phantomDir);
% Input Files
goldPDBFile = nfgGetName('goldPDBFile',phantomDir);
goldInfoFile = nfgGetName('goldInfoFile',phantomDir);
b0File = nfgGetName('b0File',phantomDir);
ctrparamsFile = nfgGetName('ctrparamsFile',phantomDir);
noisyImg = nfgGetName('noisyImg',phantomDir);
bvalsFile = nfgGetName('bvalsFile',phantomDir);
bvecsFile = nfgGetName('bvecsFile',phantomDir);
wmROIFile = nfgGetName('wmROIFile',phantomDir);

mkdir(volcheckDir);

% Get bundle ID and radius of each gold pathway from info file
gi = load(goldInfoFile);
g_bundleID = zeros(1,length(gi.strand_info));
g_radius = zeros(1,length(gi.strand_info));
for ii=1:length(gi.strand_info)
    g_bundleID(ii) = gi.strand_info(ii).bundleID+1;
    g_radius(ii) = gi.strand_info(ii).radius;
end
% Get gold fibers
fgG = mtrImportFibers(goldPDBFile);

% Radii to search through
vR = 0.01:0.01:0.04;
volG = zeros(1,length(vR));
volBM = zeros(1,length(vR));
% Save a .pdb file for each bundle
for bb=1:20
    disp(['Examining bundle ' num2str(bb) ' ...']);
    vcPDB = fullfile(volcheckDir,['b' num2str(bb) '.pdb']);
    vcSBfloat = fullfile(volcheckDir,['b' num2str(bb) '_0.SBfloat']);
    fFile = fullfile(volcheckDir,['f_b' num2str(bb) '.nii.gz']);
    eFile = fullfile(volcheckDir,['e_b' num2str(bb) '.nii.gz']);
    % Write the pdb
    fg = dtiNewFiberGroup();
    fg.fibers = fgG.fibers(g_bundleID==bb);
    if ~isempty(fg.fibers)
        mtrExportFibers(fg,vcPDB);
        % Convert it to the SBfloat format
        pParamFile = [' -i ' ctrparamsFile];
        pOutFile = [' -p ' vcSBfloat];
        pInFile = [' ' vcPDB];
        pThresh = [' --thresh ' num2str(length(fg.fibers))];
        cmd = ['contrack_score' pParamFile pOutFile pThresh ' --seq' pInFile];
        disp(cmd);
        [s,r] = system(cmd);
        for rr=1:length(vR)
            % Run error calculation for each radii
            disp(' '); disp(['Calculating for raidus ' num2str(vR(rr)) ' ...']);
            b0 = niftiRead(b0File);
            fLambda = 0;
            fDiameter = vR(rr)*2;
            argSubSize = [' -s 0,' num2str(b0.dim(1)-1) ',0,' num2str(b0.dim(2)-1) ',0,' num2str(b0.dim(3)-1)];
            argDatabase = [' -d ' vcSBfloat];
            argRaw = [' -r ' noisyImg];
            argBvals = [' --val ' bvalsFile];
            argBvecs = [' --vec ' bvecsFile];
            argB0 = [' -0 ' b0File];
            argMatter = [' -m ' wmROIFile];
            argLambda = [' -w ' num2str(fLambda)];
            argDiameter = [' --diameter ' num2str(fDiameter)];
            argEFile = [' -e ' eFile];
            argFFile = [' --fraction ' fFile];
            argError = [argDatabase argRaw argBvals argBvecs argB0 argMatter argSubSize argLambda argDiameter argEFile argFFile];
            cmdError = ['trueError ' argError];
            disp(cmdError);
            [s,r] = system(cmdError);
            % Get BlueMatter volume calculation
            f = niftiRead(fFile);
            volBM(rr) = volBM(rr)+sum(f.data(:));
            % Get gold volume calculation
            volG(rr) = volG(rr)+volBundle(fg.fibers,vR(rr));
        end
    end
end

% Plot the volume functions
figure; plot(vR,volBM,'b'); hold on; plot(vR,volG,'g');

return;

function [arcL] = arclength(fc)
arcL = sum(sqrt(sum((fc(:,2:end) - fc(:,1:end-1)).^2,1)));
return;

function vol = volBundle(fibers,radius)
lenFibers=0;
for ff=1:length(fibers)
    lenFibers = lenFibers + arclength(fibers{ff});
end
vol = lenFibers*pi*radius^2;
return;