text:
function ROI = grayLineROI(nodes, edges, startPoint, endPoint, mmPerVox, mask);
% Create a line ROI representing the shortest path ("geodesic") between two
% gray nodes. Requires mrManDist.
%
%    ROI = grayLineROI([nodes, edges], startPoint, endPoint, [mmPerVox=1,1,1], [mask]);
%
% This function simply adapts some test code provided into the comments
% for mrManDist, to create a separate function to find the geodesic.
%
% Returns a mrVista-format ROI. (See roiCreate1 for format.)
%
% INPUTS:
% nodes, edges: gray nodes and edges from a gray graph. You can also
% provide a path to a NIFTI classification file or mrGray gray graph as the
% nodes argument, and it will load/grow them from the file. [default:
% prompt for a file.]
%
% startPoint: node index or 3x1 [x y z] volume coords of the start point for
% the line ROI.
%
% endPoint: node index or 3x1 [x y z] volume coords of the end point for
% the line ROI.
%
% mask: optional [1 x N] mask vector, where N is the number of nodes. Will
% only look for the geodesic along nodes where the mask is 1 (true).
%
% OUTPUTS:
%
% ROI: mrVista-format ROI containing the line coordinates.
%
% NOTE: There is an error on some versions of MATLAB (such as 2007a,b on
% Windows) where mrManDist returns an invalid set of distances -- it
% doesn't error, but also doesn't measure distances. If that's the case,
% this function will not work.
%
%
% ras, 07/2009. I implemented this separate from the mrVista views, so it
% will hopefully be more portable in the future.
if notDefined('nodes') & notDefined('edges')
	[pth ok] = mrvSelectFile('r', {'gray' 'Gray' 'nii' 'nii.gz'}, ...
					    'Select a gray graph or NIFTI file');
	if ~ok, disp('Aborted.'); return; end
	[p f ext] = fileparts(pth);
	if strncmp( lower(ext), '.gray', 5 )
		[nodes edges] = readGrayGraph(pth);
	elseif strncmp( lower(ext), 'nii', 4 );
		%% TODO: support NIFTI
	else
		error('Invalid file type %s.', pth);
	end
end

if notDefined('startPoint') | notDefined('endPoint')
	error('Need start and end points.')
end

if notDefined('mmPerVox'),	mmPerVox = [1 1 1];		end

%% TODO: parse 3x1 start/end point specifications into nodes
if length(startPoint)==3
	startPoint = nearestGrayNode(nodes, startPoint);
end

if length(endPoint)==3
	endPoint = nearestGrayNode(nodes, endPoint);
end

%% core part of code: get geodisc using mrManDist.
% find distances from start point, and connected list of points from start
% point.
[dist npts lastPoint] = mrManDist( double(nodes), double(edges), ...
								   startPoint, mmPerVox, -1, 0 );

% track back along the 'lastPoint' list, from the start point to the end
% point, to construct the geodesic:
geodesic = [];
nextPoint = endPoint;
while nextPoint ~= startPoint
	geodesic(end+1) = lastPoint(nextPoint);
	nextPoint = lastPoint(nextPoint);
end

%% create line ROI from geodesic.
ROI = roiCreate1;
ROI.name = sprintf('Line ROI %s %s', num2str(startPoint), num2str(endPoint));
ROI.color = [1 1 0];
ROI.coords = nodes([2 1 3],geodesic);
ROI.comments = [sprintf('Created by %s %s.\n', mfilename, datestr(now)) ...
				sprintf('Start point: %s', num2str(startPoint)) ...
				sprintf('End point: %s', num2str(endPoint))];

return
% /--------------------------------------------------------------/ %

% /--------------------------------------------------------------/ %
function node = nearestGrayNode(nodes, pt);
% find the index of the gray node nearest the point at 3x1 coordinate pt.
% we assume here the point is specified using mrVista conventions: [axi cor
% sag].

%% first, check if the point is in the nodes list.
% the mrVista coordinate specification is [axi cor sag], while the
% convention for the gray nodes is [cor axi sag]. Hence the flip.
node = find( nodes(2,:)==pt(1) & nodes(1,:)==pt(2) & nodes(3,:)==pt(3) );

%% if the node is found, we're done. Otherwise, find the nearest gray node.
if isempty(node)
	if prefsVerboseCheck >= 1
		fprintf('[%s]: looking for nearest gray node to %s...', ...
				mfilename, num2str(pt));
	end

% 	% this doesn't seem to work...
% 	[node bestSqDist] = nearpoints(double(pt), nodes([2 1 3],:));

	% use the algorithm I use in segGet to find the nearest node
	tolerance = 5; % mm
	dist = sqrt( [nodes(2,:) - pt(1)] .^ 2 + ...
		[nodes(1,:) - pt(2)] .^ 2 + ...
		[nodes(3,:) - pt(3)] .^ 2 );
	if min(dist) > tolerance
		warning(sprintf('[%s]: no node found within tolerance [%i mm]', ...
			mfilename, tolerance));
		varargout{1} = [];
	else
		I = find(dist==min(dist));
		node = I(1);
	end

	if prefsVerboseCheck >= 1
		fprintf('done. Best distance: %.1f.\n', min(dist));
	end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ p, Fstat, df1, df2 ] = ftest(n,np1,np2,chi1,chi2)
% [ p ] = ftest(n,np1,np2,chi1,chi2)
% function to test whether addition of model parameters to
%     fit data is warranted by level of misfit improvement
%
% Inputs
% n = # of data
% np1, np2 = number of free parameters for each fit.
% chi1 & chi2 = sum of squares of misfits for two models.
%     must be positive values
%     normalization of chis by n not required
%     (normalization divides out in Fstatistic)
%
% Outputs
% p = probability (between 0 - 1) that improvement to fit from
%     addition of parameters is due to chance. I.e.,
%     0 means certainty that extra parameters are warranted
%     1 means improvement undoubtedly attributable to chance
% If desired, will also return F-statistic (Fstat) and
%     degrees of freedom (df1, df2) for f-distribution
%
% Written by James Conder, Southern Illinois University, Oct. 2010
% Citation:
%     Anderson, K.B. & Conder, J.A., Discussion of Multicyclic Hubbert
% Modeling as a Method for Forecasting Future Petroleum Production,
% Energy & Fuels, dx.doi.org/10.1021/ef1012648, 2011
%
% Update June, 2011
%   Check added for df1=1 (giving -Inf). Use MatrixLabs approximation
%   for f when df1=1. Accuracy ok, but some small differences against
%   fcdf
%
% Update May 31, 2012
%   Allow Fstat, df1, & df2 as outputs.
%   Use fcdf from matlab statistics toolbox if available (fast)
%   Some cosmetic clean up
%
% Comments & questions should be directed to conder@siu.edu
%----------------------

%%% check ordering. np2 should be > np1 and chi2 should be < chi1
% i.e., second model should have more parameters and better fit
if np2 == np1
    disp('number of model parameters are the same in both cases!')
    p = 1;
    return
elseif np2 < np1        % np1 should be less than np2. If not, just swap.
    nptemp = np1; np1 = np2; np2 = nptemp;
    chitemp = chi1; chi1 = chi2; chi2 = chitemp;
end

if chi2 >= chi1
    disp('misfit higher for model with more parameters!')
    p = 1;
    return
end

%%% number of degrees of freedom for f-distribution
df1 = np2 - np1;		% number of degrees of freedom
df2 = n - np2 - 1;

%%% F-statistic
Fstat = df2*(chi1 - chi2)/(df1*chi2);

%%% find p by determination of cumulative f-distribution at Fstat
% first check if fcdf from matlab statistics toolbox is present (fast).
% if not, numerically integrate f-distribution.
% cdff is equivalent to result from fcdf in Matlab statistics package.

if exist('fcdf.m','file') == 2       % check for availability of fcdf from statistics toolbox
    p = 1 - fcdf(Fstat,df1,df2);
else
    if df1 ~= 1         % numerically integrate f-distribution
        ifpt = 1000001; % large number of slices for accurate numerical integration
        dx = Fstat/(ifpt-1);
        x = 0:dx:1.2*Fstat;

        fnumgam = gammaln((df1+df2)/2);		% gamma func factors can be very large, use ln
        fdengam = gammaln(df1/2) + gammaln(df2/2);
        fgam = exp(fnumgam - fdengam);

        fnum = fgam*((df1/df2)^(df1/2)).*(x.^(0.5*df1 -1));
        fden = ((1 + df1*x/df2).^(0.5*(df1+df2)));
        f = fnum./fden;	% f distribution for df1, df2
        %fF = f(ifpt);		% f at Fstat

        cdff = cumsum(f)*dx;	% numerical integration of f distribution
        p = 1 - cdff(ifpt);
    else    % when df1=1 use F-dist approximation from MatrixLabs to avoid -Inf
        if Fstat <= 0.5    % Compute using inverse for small F-values
            s = df2;
            t = df1;
            z = 1/Fstat;
        else
            s = df1;
            t = df2;
            z = Fstat;
        end
        j = 2/(9*s);
        k = 2/(9*t);

        % Use approximation formulas
        y = abs((1 - k)*z^(1/3) - 1 + j)/sqrt(k*z^(2/3) + j);
        if t < 4
            y = y*(1 + 0.08*y^4/t^3);
        end

        a1 = 0.196854;
        a2 = 0.115194;
        a3 = 0.000344;
        a4 = 0.019527;
        p = 0.5/(1 + y*(a1 + y*(a2 + y*(a3 + y*a4))))^4;

        % Adjust if inverse was computed
        if Fstat <= 0.5
            p = 1 - p;
        end
     end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [fSet] = featureSetReduction_LGG(pathMINE,outcome,setSize,nonTextStruct,textCells,textCellsName,paramAll,paramUsed,baseline,alpha,delta,nBoot,seed,batchNum)
% -------------------------------------------------------------------------
% function [fSet] = featureSetReduction_LGG(pathMINE,outcome,setSize,nonTextStruct,textCells,textCellsName,paramAll,paramUsed,baseline,alpha,delta,nBoot,seed,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set reduction for an experiment with a
% specific degree of freedom on texture extraction parameters as defined by
% 'paramUsed', according to the methodology described in ref. [1].
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathMINE: Full path to the MINE.jar executable. The executable can be
%              downloaded at: <http://www.exploredata.net/>.
% 2. outcome: Column vector of size [nInst X 1] specifying the outcome status
%             (1 or 0) for all instances.
% 3. setSize: Size of the output feature set (typically set to 25 in ref. [1]).
% 4. nonTextStruct: Structure data for non-texture features. This structure
%                   is of the same format as the one saved as output to
%                   computeAllNonTextureFeatures_LGG.m, for example.
% 5. textCells: Cell vector of organized texture cells for all the
%               different scans. Format: textCells = {cell1, cell2, etc.},
%               where cell1, cell2, etc. are the files saved as output to
%               organizeData_LGG.m
% 6. textCellsName: Cell of strings corresponding to the name of the
%                   corresponding cells in textCells.
% 7. paramAll: Cell vector incorporating all texture extraction parameters
%              tested in textCells. See EXAMPLE below for more details.
% 8. paramUsed: Vector of 1's and 0's to specify the degree of freedom on
%               texture extraction parameters for the current experiment.
%               For example, for an experiment where extraction parameters
%               1, 2 and 3 in paramAll are allowed to vary, use
%               paramUsed = [1,1,1].
% 9. baseline: Vector of numerical values specifying the baseline texture
%              extraction parameters for each entry in paramAll. See EXAMPLE
%              below for more details.
% 10. alpha: Numerical values specifying the coefficient of the first part of
%            the Gain equation, as defined in ref. [1].
% 11. delta: Numerical values specifying the coefficient of the second part
%            of the Gain equation, as defined in ref. [1] (third part is set
%            to 0 in this function).
% 12. nBoot: Number of bootstrap samples to use.
% 13. batchNum: (optional input). If present, integer that specifies the
%               batch number for parallelization purposes
%
% See masterScript_LGG.m for a complete example of how to utilize the
% current function.
% -------------------------------------------------------------------------
% OUTPUTS:
% 1. fSET: Structure specifying the resulting feature set for the current
%          experiment.
%          --> fSET.Data: Array of size [nPatient X setSize], specifying the
%                         numerical data of the chosen features, where 'nInst'
%                         refers to the number of instances.
%          --> fSET.Info: Vector of cells with strings specifying information
%                         about the chosen features.
% -------------------------------------------------------------------------
% EXAMPLE:
% scale_mat = [1,2,3,4,5];
% algo_cell = {'Equal','Uniform'};
% Ng_mat = [8,16,32,64];
%
% paramAll = {scale_mat,algo_cell,Ng_mat};
% paramUsed = [1 1 1]; (example for a given experiment)
% baseline = [1,2,3];
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% RANDOM NUMBER GENERATOR SEED
rng(seed);

% INITIALIZATION
nInst = numel(outcome); % Number of patients
nScan = numel(textCellsName); % Number of scans included in that feature set experiment
fSet = struct; fSet.Data=zeros(nInst,setSize); fSet.Info=cell(setSize,1); % Initializing final feature set structure
if nargin == 15
    micName = ['batch',num2str(batchNum)];
else
    micName = 'master';
end
if ~isempty(strfind(textCellsName{1},'CT'))
    name = 'CT';
else
    name = 'MR';
end

% Initializing names
nonTextName = fieldnames(nonTextStruct); nNonText = numel(nonTextName);
textTypeName = fieldnames(textCells{1}{1}); nTextType = numel(textTypeName);
nameF = cell(1,nTextType);
nText = 0;
for t = 1: nTextType
    nameF{t} = fieldnames(textCells{1}{1}.(textTypeName{t}));
    nText = nText + numel(nameF{t});
end
textName = cell(nText,1);
count = 1;
for t = 1:nTextType
    for f = 1:numel(nameF{t})
        textName{count} = [textTypeName{t},'-',nameF{t}{f}];
        count = count +1;
    end
end
nFeatures = nNonText + nScan*nText;
cellNames = cell(nFeatures,2);
count = 1;
for i = 1:nScan
    for j = 1:nText
        cellNames{count,1} = textCellsName{i};
        cellNames{count,2} = textName{j};
        count = count + 1;
    end
end
for i = 1:nNonText
    cellNames{count,1} = nonTextName{i};
    count = count + 1;
end

% Extraction parameters used in that experiment
nParamType = numel(paramUsed);
param = cell(1,nParamType);
paramSize = ones(1,nParamType);
nParam = 1;
for i = 1:nParamType
    if paramUsed(i)
        temp = length(paramAll{i});
        nParam = nParam * temp;
        paramSize(i) = temp;
        param{i} = paramAll{i};
    end
end

% FILLING UP WORKING DATA MATRIX
matData = zeros(nFeatures,nParam,nInst); % Matrix containing all the data vectors for each feature (used for part 2 of Gain equation)
count = 1;
for i = 1:nScan
    textCell = textCells{i};
    textCell = removeParam(textCell,paramUsed,baseline);
    for t = 1:nTextType
        for f = 1:numel(nameF{t})
            for p = 1:numel(textCell)
                matData(count,p,:) = textCell{p}.(textTypeName{t}).(nameF{t}{f}).Data;
            end
            count = count + 1;
        end
    end
end
for i = 1:nNonText
    matData(nText*nScan + i,1,:) = nonTextStruct.(nonTextName{i}).Data;
    for j = 2:nParam
        matData(nText*nScan + i,j,:) = matData(nText*nScan + i,1,:);
    end
end

% IMPORTANT: IF NAN, IT WOULD MEAN SOMETHING WENT WRONG IN SOME TEXTURE CALCULATIONS
matData(isnan(matData(:))) = 0;

% GETTING BOOTSTRAP SAMPLES TO BE USED FOR ALL EXPERIMENTS (using imbalance-adjusted resampling)
[bootSam,~] = buildBootSet(outcome,nBoot,'adjust');

% GETTING BOOTSTRAP RESULTS FOR THE FIRST PART OF THE GAIN EQUATION (ref. [1])
matCorr = zeros(nFeatures,nParam); % Matrix of Spearman correlation of each feature with the outcome (part 1 of Gain equation)
if sum(paramUsed(:)) % Need a transpose
    for n = 1:nBoot
        dataBoot = matData(:,:,bootSam(:,n));
        outcomeBoot = outcome(bootSam(:,n));
        for i = 1:size(dataBoot,1)
            matCorr(i,:) = matCorr(i,:) + corr(squeeze(dataBoot(i,:,:))',outcomeBoot,'type','Spearman','rows','pairwise')';
        end
    end
else % No transpose
    for n = 1:nBoot
        dataBoot = matData(:,:,bootSam(:,n));
        outcomeBoot = outcome(bootSam(:,n));
        for i = 1:size(dataBoot,1)
            matCorr(i,:) = matCorr(i,:) + corr(squeeze(dataBoot(i,:,:)),outcomeBoot,'type','Spearman','rows','pairwise')';
        end
    end
end
matCorr = abs(matCorr./nBoot);

% CHOOSING FIRST FEATURE(depends only on Spearman's correlation)
indChosen = zeros(setSize-1,1);

% Choosing feature
[~,index] = max(matCorr(:));
[row,col] = ind2sub([nFeatures,nParam],index);
indChosen(1) = row;
fSet.Data(:,1) = squeeze(matData(row,col,:));

% Obtaining feature name
if row <= nScan*nText
    [ind1,ind2,ind3] = ind2sub(paramSize,col);
    indBase = [ind1,ind2,ind3];
    indFinal = zeros(1,3);
    for i = 1:3
        if isempty(param{i})
            indFinal(i) = baseline(i);
        else
            indFinal(i) = indBase(i);
        end
    end
    string = [cellNames{row,1},'(R=',num2str(1,'%.2f'),',Scale=',num2str(paramAll{1}(indFinal(1))),',Quant.algo=',paramAll{2}{indFinal(2)},',Ng=',num2str(paramAll{3}(indFinal(3))),')--',cellNames{row,2}];
else % This is a nonTexture feature
    string = cellNames{row,1};
end
fSet.Info{1} = string;

% COMPUTING FOR OTHER FEATURES
PICtest = zeros(nFeatures,setSize-1);
for f = 1:setSize-1

    % Computing PIC for all bootstrap samples
    for n = 1:nBoot
        varBoot = fSet.Data(bootSam(:,n),f);
        dataBoot = matData(:,:,bootSam(:,n));
        dataBootAv = squeeze(mean(dataBoot,2))';
        try
            PICtest(:,f) = PICtest(:,f) + (1 - applyMIC(pathMINE,varBoot,dataBootAv,micName));
        catch % SOLVE THAT PROBLEM
            PICtest(:,f) = PICtest(:,f) + zeros(nFeatures,1);
        end
    end
    PICtest(:,f) = PICtest(:,f)./nBoot;
    PICtemp = zeros(nFeatures,1);
    for k = 1:f
        PICtemp = PICtemp + 2*(f-k+1)/(f*(f+1)).*PICtest(:,k);
    end
    PIC = repmat(PICtemp,1,nParam);

    % Choosing feature
    Gain = alpha.*matCorr + delta.*PIC;
    [~,index] = sort(Gain(:),'descend'); best = 1;
    [row,col] = ind2sub([nFeatures,nParam],index(best));
    while ~isempty(find(indChosen==row))
        best = best + 1;
        [row,col] = ind2sub([nFeatures,nParam],index(best));
    end
    indChosen(f+1) = row;
    fSet.Data(:,f+1) = squeeze(matData(row,col,:));

    % Obtaining feature name
    if row <= nScan*nText
        [ind1,ind2,ind3] = ind2sub(paramSize,col);
        indBase = [ind1,ind2,ind3];
        indFinal = zeros(1,3);
        for i = 1:3
            if isempty(param{i})
                indFinal(i) = baseline(i);
            else
                indFinal(i) = indBase(i);
            end
        end
        string = [cellNames{row,1},'(R=',num2str(1,'%.2f'),',Scale=',num2str(paramAll{1}(indFinal(1))),',Quant.algo=',paramAll{2}{indFinal(2)},',Ng=',num2str(paramAll{3}(indFinal(3))),')--',cellNames{row,2}];
    else % This is a nonTexture feature
        string = cellNames{row,1};
    end
    fSet.Info{f+1} = string;
end

end

author: mvallieres

repo: radiomics


text:
function b = padarray(varargin)
%PADARRAY Pad an array.
%   B = PADARRAY(A,PADSIZE) pads array A with PADSIZE(k) number of zeros
%   along the k-th dimension of A.  PADSIZE should be a vector of
%   positive integers.
%
%   B = PADARRAY(A,PADSIZE,PADVAL) pads array A with PADVAL (a scalar)
%   instead of with zeros.
%
%   B = PADARRAY(A,PADSIZE,PADVAL,DIRECTION) pads A in the direction
%   specified by the string DIRECTION.  DIRECTION can be one of the
%   following strings.
%
%       String values for DIRECTION
%       'pre'         Pads before the first array element along each
%                     dimension .
%       'post'        Pads after the last array element along each
%                     dimension.
%       'both'        Pads before the first array element and after the
%                     last array element along each dimension.
%
%   By default, DIRECTION is 'post'.
%
%   B = PADARRAY(A,PADSIZE,METHOD,DIRECTION) pads array A using the
%   specified METHOD.  METHOD can be one of these strings:
%
%       String values for METHOD
%       'circular'    Pads with circular repetion of elements.
%       'replicate'   Repeats border elements of A.
%       'symmetric'   Pads array with mirror reflections of itself.
%
%   Class Support
%   -------------
%   When padding with a constant value, A can be numeric or logical.
%   When padding using the 'circular', 'replicate', or 'symmetric'
%   methods, A can be of any class.  B is of the same class as A.
%
%   Example
%   -------
%   Add three elements of padding to the beginning of a vector.  The
%   padding elements contain mirror copies of the array.
%
%       b = padarray([1 2 3 4],3,'symmetric','pre')
%
%   Add three elements of padding to the end of the first dimension of
%   the array and two elements of padding to the end of the second
%   dimension.  Use the value of the last array element as the padding
%   value.
%
%       B = padarray([1 2; 3 4],[3 2],'replicate','post')
%
%   Add three elements of padding to each dimension of a
%   three-dimensional array.  Each pad element contains the value 0.
%
%       A = [1 2; 3 4];
%       B = [5 6; 7 8];
%       C = cat(3,A,B)
%       D = padarray(C,[3 3],0,'both')
%
%   See also CIRCSHIFT, IMFILTER.

%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.11 $  $Date: 2002/03/15 15:28:49 $
%
%   This function is part of the Image Processing Toolbox.
%   At the Householder Meeting in Peebles, Scotland, May 17--21, 2002,
%   Cleve Moler has given us permission to include this function,
%   and mkconstarray.m, in our toolbox.  We have modified the
%   subfunction ParseInputs.m, to remove some dependencies on
%   other

[a, method, padSize, padVal, direction] = ParseInputs(varargin{:});

if isempty(a),% treat empty matrix similar for any method

   if strcmp(direction,'both')
      sizeB = size(a) + 2*padSize;
   else
      sizeB = size(a) + padSize;
   end

   b = mkconstarray(class(a), padVal, sizeB);

else
  switch method
    case 'constant'
        b = ConstantPad(a, padSize, padVal, direction);

    case 'circular'
        b = CircularPad(a, padSize, direction);

    case 'symmetric'
        b = SymmetricPad(a, padSize, direction);

    case 'replicate'
        b = ReplicatePad(a, padSize, direction);
  end
end

if (islogical(a))
    b = logical(b);
end

%%%
%%% ConstantPad
%%%
function b = ConstantPad(a, padSize, padVal, direction)

numDims = prod(size(padSize));

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
sizeB = zeros(1,numDims);
for k = 1:numDims
    M = size(a,k);
    switch direction
        case 'pre'
            idx{k}   = (1:M) + padSize(k);
            sizeB(k) = M + padSize(k);

        case 'post'
            idx{k}   = 1:M;
            sizeB(k) = M + padSize(k);

        case 'both'
            idx{k}   = (1:M) + padSize(k);
            sizeB(k) = M + 2*padSize(k);
    end
end

% Initialize output array with the padding value.  Make sure the
% output array is the same type as the input.
b         = mkconstarray(class(a), padVal, sizeB);
b(idx{:}) = a;

%%%
%%% CircularPad
%%%
function b = CircularPad(a, padSize, direction)

numDims = prod(size(padSize));

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
  M = size(a,k);
  dimNums = [1:M];
  p = padSize(k);

  switch direction
    case 'pre'
       idx{k}   = dimNums(mod([-p:M-1], M) + 1);

    case 'post'
      idx{k}   = dimNums(mod([0:M+p-1], M) + 1);

    case 'both'
      idx{k}   = dimNums(mod([-p:M+p-1], M) + 1);

  end
end
b = a(idx{:});

%%%
%%% SymmetricPad
%%%
function b = SymmetricPad(a, padSize, direction)

numDims = prod(size(padSize));

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
  M = size(a,k);
  dimNums = [1:M M:-1:1];
  p = padSize(k);

  switch direction
    case 'pre'
      idx{k}   = dimNums(mod([-p:M-1], 2*M) + 1);

    case 'post'
      idx{k}   = dimNums(mod([0:M+p-1], 2*M) + 1);

    case 'both'
      idx{k}   = dimNums(mod([-p:M+p-1], 2*M) + 1);
  end
end
b = a(idx{:});

%%%
%%% ReplicatePad
%%%
function b = ReplicatePad(a, padSize, direction)

numDims = prod(size(padSize));

% Form index vectors to subsasgn input array into output array.
% Also compute the size of the output array.
idx   = cell(1,numDims);
for k = 1:numDims
  M = size(a,k);
  p = padSize(k);
  onesVector = ones(1,p);

  switch direction
    case 'pre'
      idx{k}   = [onesVector 1:M];

    case 'post'
      idx{k}   = [1:M M*onesVector];

    case 'both'
      idx{k}   = [onesVector 1:M M*onesVector];
  end
end
 b = a(idx{:});

%%%
%%% ParseInputs
%%%
function [a, method, padSize, padVal, direction] = ParseInputs(varargin)

% default values
a         = [];
method    = 'constant';
padSize   = [];
padVal    = 0;
direction = 'both';

%checknargin(2,4,nargin,mfilename);

a = varargin{1};

padSize = varargin{2};
%checkinput(padSize, {'double'}, {'real' 'vector' 'nonnan' 'nonnegative' ...
%                    'integer'}, mfilename, 'PADSIZE', 2);

% Preprocess the padding size
if (prod(size(padSize)) < ndims(a))
    padSize           = padSize(:);
    padSize(ndims(a)) = 0;
end

if nargin > 2

    firstStringToProcess = 3;

    if ~ischar(varargin{3})
        % Third input must be pad value.
        padVal = varargin{3};
%        checkinput(padVal, {'numeric' 'logical'}, {'scalar'}, ...
%                   mfilename, 'PADVAL', 3);

        firstStringToProcess = 4;

    end

    for k = firstStringToProcess:nargin
%        validStrings = {'circular' 'replicate' 'symmetric' 'pre' ...
%                        'post' 'both'};
%        string = checkstrs(varargin{k}, validStrings, mfilename, ...
%                           'METHOD or DIRECTION', k);
        string = varargin{k};
        switch string
         case {'circular' 'replicate' 'symmetric'}
          method = string;

         case {'pre' 'post' 'both'}
          direction = string;

         otherwise
          error('Images:padarray:unexpectedError', '%s', ...
                'Unexpected logic error.')
        end
    end
end

% Check the input array type
if strcmp(method,'constant') & ~(isnumeric(a) | islogical(a))
    id = sprintf('Images:%s:badTypeForConstantPadding', mfilename);
    msg1 = sprintf('Function %s expected A (argument 1)',mfilename);
    msg2 = 'to be numeric or logical for constant padding.';
    error(id,'%s\n%s',msg1,msg2);
end

repo: IRtools


text:
function varargout = createSoccerBall()
%CREATESOCCERBALL Create a 3D mesh representing a soccer ball.
%
%   It is basically a wrapper of the 'bucky' function in matlab.
%   [V, E, F] = createSoccerBall
%   return vertices, edges and faces that constitute a soccerball
%   V is a 60-by-3 array containing vertex coordinates
%   E is a 90-by-2 array containing indices of neighbor vertices
%   F is a 32-by-1 cell array containing vertex indices of each face
%   Example
%   [v, f] = createSoccerBall;
%   drawMesh(v, f);
%
%   See also
%   meshes, drawMesh, bucky

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2006-08-09
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   2007-01-04 remove unused variables, enhance output processing
%   2010-12-07 clean up edges, uses formatMeshOutput

% get vertices and adjacency matrix of the buckyball
[b, n] = bucky;

% compute edges
[i, j] = find(b);
e = [i j];
e = unique(sort(e, 2), 'rows');

% compute polygons that correspond to each 3D face
f = minConvexHull(n)';

% format output
varargout = formatMeshOutput(nargout, n, e, f);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function K = kernel(X1, X2, param)

n1 = size(X1,1);
n2 = size(X2,1);
d = size(X1,2);
switch param.type
    case 'linear'
        K = X1*X2';

    case 'rbf'
        norm1 = sum(X1.^2,2);
        norm2 = sum(X2.^2,2);
        dist = (repmat(norm1 ,1,size(X2,1)) + ...
            repmat(norm2',size(X1,1),1) - ...
            2*X1*X2');
        K = exp(-0.5/param.sig^2 * dist);

    case 'dist'
        norm1 = sum(X1.^2,2);
        norm2 = sum(X2.^2,2);
        dist = (repmat(norm1 ,1,size(X2,1)) + ...
            repmat(norm2',size(X1,1),1) - ...
            2*X1*X2');
        K = param.maxdist - dist;

    case 'Lp'
        K = zeros(n1,n2);
        for n = 1:n1
            K(n,:) = real(sum(bsxfun(@minus, X2, X1(n,:)).^param.p,2));
        end
        K = exp(-0.5/param.sig^2 * K);

    case 'L1'
        K = zeros(n1,n2);
        for n = 1:n1
            K(n,:) = real(sum(abs(bsxfun(@minus, X2, X1(n,:))),2));
        end
        K = exp(-0.5/param.sig * K);

    case 'histintersection'
		K = zeros(n1,n2);
        for n = 1:n1
            K(n,:) = sum(bsxfun(@min, X2, X1(n,:)),2);
        end

    case 'L1_james'
        K = zeros(n1,n2);
        for n = 1:n1
            K(n,:) = real(sum(abs(bsxfun(@minus, X2, X1(n,:))),2));
        end
        K = max(max(K)) - K;

    otherwise
        error('Unknown kernel');
end;

author: CSAILVision

repo: LabelMeToolbox


text:
function [ ptab, b, c, d, eps, ierror ] =  least_set_old ( ...
  ntab, xtab, ytab, ndeg )

%*****************************************************************************80
%
%% LEAST_SET_OLD constructs the least squares polynomial approximation to data.
%
%  Discussion:
%
%    The least squares polynomial is not returned directly as a simple
%    polynomial.  Instead, it is represented in terms of a set of
%    orthogonal polynomials appopriate for the given data.  This makes
%    the computation more accurate, but means that the user can not
%    easily evaluate the computed polynomial.  Instead, the routine
%    LEAST_EVAL should be used to evaluate the least squares polynomial
%    at any point.  (However, the value of the least squares polynomial
%    at each of the data points is returned as part of this computation.)
%
%
%    A discrete unweighted inner product is used, so that
%
%      ( F(X), G(X) ) = sum ( 1 <= I <= NTAB ) F(XTAB(I)) * G(XTAB(I)).
%
%    The least squares polynomial is determined using a set of
%    orthogonal polynomials PHI.  These polynomials can be defined
%    recursively by:
%
%      PHI(0)(X) = 1
%      PHI(1)(X) = X - B(1)
%      PHI(I)(X) = ( X - B(I) ) * PHI(I-1)(X) - D(I) * PHI(I-2)(X)
%
%    The array B(1:NDEG) contains the values
%
%      B(I) = ( X*PHI(I-1), PHI(I-1) ) / ( PHI(I-1), PHI(I-1) )
%
%    The array D(2:NDEG) contains the values
%
%      D(I) = ( PHI(I-1), PHI(I-1) ) / ( PHI(I-2), PHI(I-2) )
%
%    Using this basis, the least squares polynomial can be represented as
%
%      P(X)(I) = sum ( 0 <= I <= NDEG ) C(I) * PHI(I)(X)
%
%    The array C(0:NDEG) contains the values
%
%      C(I) = ( YTAB(I), PHI(I) ) / ( PHI(I), PHI(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Gisela Engeln-Muellges and Frank Uhlig,
%    Numerical Algorithms with C, pages 191-193.
%    Springer, 1996.
%
%  Parameters:
%
%    Input, integer NTAB, the number of data points.
%
%    Input, real XTAB(NTAB), the X data.  The values in XTAB
%    should be distinct, and in increasing order.
%
%    Input, real YTAB(NTAB), the Y data values corresponding
%    to the X data in XTAB.
%
%    Input, integer NDEG, the degree of the polynomial which the
%    program is to use.  NDEG must be at least 1, and less than or
%    equal to NTAB-1.
%
%    Output, real PTAB(NTAB), the value of the least squares polynomial
%    at the points XTAB(1:NTAB).
%
%    Output, real B(1:NDEG), C(1:NDEG+1), D(1:NDEG-1), arrays containing
%    data about the polynomial.
%
%    Output, real EPS, the root-mean-square discrepancy of the
%    polynomial fit.
%
%    Output, integer IERROR, error flag.
%    zero, no error occurred;
%    nonzero, an error occurred, and the polynomial could not be computed.
%
  ierror = 0;
  C_OFFSET = 1;
  D_OFFSET = -1;
%
%  Check NDEG.
%
  if ( ndeg < 1 )
    ierror = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEAST_SET_OLD - Fatal error!\n' );
    fprintf ( 1, '  NDEG < 1.\n' );
    error ( 'LEAST_SET_OLD - Fatal error!' );
  end

  if ( ntab <= ndeg )
    ierror = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEAST_SET_OLD - Fatal error!\n' );
    fprintf ( 1, '  NTAB <= NDEG.\n' );
    error ( 'LEAST_SET_OLD - Fatal error!' );
  end
%
%  Check that the abscissas are strictly increasing.
%
  for i = 1 : ntab-1
    if ( xtab(i+1) <= xtab(i) )
      ierror = 1;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LEAST_SET_OLD - Fatal error!\n' );
      fprintf ( 1, '  XTAB must be strictly increasing, but\n' );
      fprintf ( 1, '  XTAB(%d) = %f\n', i,   xtab(i)   );
      fprintf ( 1, '  XTAB(%d) = %f\n', i+1, xtab(i+1) );
      error ( 'LEAST_SET_OLD - Fatal error!' );
    end
  end

  i0l1 = 0;
  i1l1 = ntab;
%
%  The polynomial is of degree at least 0.
%
  y_sum = sum ( ytab(1:ntab) );
  rn0 = ntab;
  c(0+C_OFFSET) = y_sum / ntab;

  ptab(1:ntab) = y_sum / ntab;

  if ( ndeg == 0 )
    eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
    eps = sqrt ( eps / ntab );
    b = [];
    d = [];
    return;
  end
%
%  The polynomial is of degree at least 1.
%
  b(1) = sum ( xtab(1:ntab) ) / ntab;

  s = 0.0E+00;
  sum2 = 0.0E+00;
  for i = 1 : ntab
    ztab(i1l1+i) = xtab(i) - b(1);
    s = s + ztab(i1l1+i)^2;
    sum2 = sum2 + ztab(i1l1+i) * ( ytab(i) - ptab(i) );
  end

  rn1 = s;
  c(1+C_OFFSET) = sum2 / s;

  for i = 1 : ntab
    ptab(i) = ptab(i) + c(1+C_OFFSET) * ztab(i1l1+i);
  end

  if ( ndeg == 1 )
    eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
    eps = sqrt ( eps / ntab );
    d = [];
    return;
  end

  ztable(1:ntab) = 1.0;

  mdeg = 2;
  k = 2;

  while ( 1 )

    d(k+D_OFFSET) = rn1 / rn0;

    sum2 = 0.0;
    for i = 1 : ntab
      sum2 = sum2 + xtab(i) * ztab(i1l1+i) * ztab(i1l1+i);
    end

    b(k) = sum2 / rn1;

    s = 0.0;
    sum2 = 0.0;
    for i = 1 : ntab
      ztab(i0l1+i) = ( xtab(i) - b(k) ) * ztab(i1l1+i) ...
        - d(k+D_OFFSET) * ztab(i0l1+i);
      s = s + ztab(i0l1+i) * ztab(i0l1+i);
      sum2 = sum2 + ztab(i0l1+i) * ( ytab(i) - ptab(i) );
    end

    rn0 = rn1;
    rn1 = s;

    c(k+C_OFFSET) = sum2 / rn1;

    it = i0l1;
    i0l1 = i1l1;
    i1l1 = it;

    for i = 1 : ntab
      ptab(i) = ptab(i) + c(k+C_OFFSET) * ztab(i1l1+i);
    end

    if ( ndeg <= mdeg )
      break;
    end

    mdeg = mdeg + 1;
    k = k + 1;

  end
%
%  Compute the RMS error.
%
  eps = sum ( ( ptab(1:ntab) - ytab(1:ntab) ).^2 );
  eps = sqrt ( eps / ntab );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function printOptions(fid, expInfo)
%PRINTOPTIONS   Print problem options when they are exported.
%
% Calling sequence:
%   PRINTOPTIONS(FID, EXPINFO)
% where
%   FID:        ID of a file-writing stream.
%   EXPINFO:    Struct containing information for printing the problem.

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

% Extract info from the expInfo struct:
tol = expInfo.tol;
pdeflag = expInfo.pdeflag;
doplot = expInfo.doplot;
dohold = expInfo.dohold;
ylim1 = expInfo.ylim1;
ylim2 = expInfo.ylim2;
fixN = expInfo.fixN;

% Start with an empty struct for opts
opts = [];

% Option for tolerance
opts = [opts, '''Eps'', ', tol];

% Is PDEFLAG set?
if ( ~all(pdeflag) )
    opts = [opts, ', ''PDEflag'', ', 'pdeflag'];
end

% Options for plotting
if ( strcmpi(doplot, 'off') )
    opts = [opts, ', ''Plot'', ', '''off'''];
else
    if ( dohold )
        opts = [opts, ', ''HoldPlot'', ', '''on'''];
    end
    if ( ~isempty(ylim1) && ~isempty(ylim2) )
        opts = [opts, ', ''Ylim'', [', ylim1, ',', ylim2,']'];
    end
end

% Options for fixed N
if ( ~isempty(fixN) )
    N = str2double(fixN);
    opts = [opts, ',''N'',', N];
end

% Set up preferences

fprintf(fid, '\n%%%% Setup preferences for solving the problem.\n');
fprintf(fid, 'opts = pdeset');
if ( isempty(opts) )
    fprintf(fid, ';\n', opts);
else
    fprintf(fid, '(%s);\n', opts);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function X = pinv(A, tol)
%PINV   Pseudoinverse of a column CHEBFUN.
%   X = PINV(A) produces a row CHEBFUN X so that A*X*A = A and X*A*X = X.
%
%   X = PINV(A, TOL) uses the tolerance TOL. The computation uses SVD(A) and any
%   singular value less than the tolerance TOL is treated as zero.
%
% See also SVD, RANK.

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

if ( A(1).isTransposed )
    error('CHEBFUN:CHEBFUN:pinv:row', ...
        'PINV only defined for column CHEBFUN objects.')
end

% Compute the SVD:
[U, S, V] = svd(A, 0);
s = diag(S);

% Choose a tolerance if none is given:
if ( nargin == 1 )
	tol = max(length(A)*eps(max(s)), vscale(A)*eps);
end

% Compute the rank:
r = sum(s > tol);

if ( r == 0 )
	X = 0*A';
else
    U = extractColumns(U, 1:r);
    S = diag(ones(r,1)./s(1:r));
    V = V(:,1:r);
	X = V*S*U';
end


text:
function [cost, grad] = getCostGrad(problem, x, storedb, key)
% Computes the cost function and the gradient at x in one call if possible.
%
% function [cost, grad] = getCostGrad(problem, x)
% function [cost, grad] = getCostGrad(problem, x, storedb)
% function [cost, grad] = getCostGrad(problem, x, storedb, key)
%
% Returns the value at x of the cost function described in the problem
% structure, as well as the gradient at x.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: canGetCost canGetGradient getCost getGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    if isfield(problem, 'costgrad')
    %% Compute the cost/grad pair using costgrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.costgrad)
            case 1
                [cost, grad] = problem.costgrad(x);
            case 2
                % Obtain, pass along, and save the store for x.
                store = storedb.getWithShared(key);
                [cost, grad, store] = problem.costgrad(x, store);
                storedb.setWithShared(store, key);
            case 3
                % Pass along the whole storedb (by reference), with key.
                [cost, grad] = problem.costgrad(x, storedb, key);
            otherwise
                up = MException('manopt:getCostGrad:badcostgrad', ...
                    'costgrad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

    else
    %% Revert to calling getCost and getGradient separately

        cost = getCost(problem, x, storedb, key);
        grad = getGradient(problem, x, storedb, key);

    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function meshCanopy(I,Mdata,Mmap,height)
%function meshCanopy: display a mesh above a grayscale image
%SCd 11/18/2010
%
%Updates:
%   -11/22/2010: Added example (requires IPT)
%                Added height input argument
%
%Input Arguments:
%   -I: 2-dimensional grayscale image slice.  The values are expected to
%       range from 0 to 255.  If the maximum value is greater than 255 or
%       the minimum value is less than 0, it will be scaled to 0:255.  Else
%       it will remain unchanged.
%   -Mdata: 2-dimensional mesh data (Z for a meshplot).
%       NOTE: the mesh command will be called with:
%       >>[ii jj] = meshgrid(ceil(cumsum(diff([0 linspace(1,size(I,2),size(Mdata,2))]))),ceil(cumsum(diff([0 linspace(1,size(I,1),size(Mdata,1))]))));
%       >>mesh(ii,jj,Mdata);
%       and thus does not need to be the same size as the image!
%   -Mmap: string, function_handle or nx3, mesh color map.  See:
%       >>doc colormap
%       for valid options.  The argument is optional and defaults to 'jet'
%       Examples: 'jet', @jet, [0 0 1; 0.5 0;.1 .1 .1]
%   -height: scalar height of the mesh above the image so you can see both.
%       Optional defaults to 80.
%
%Output Arguments:
%   -None!
%
%Example: (Requires the Image Processing Toolbox)
%   %Display a Mesh Canopy of a standard deviation image, above the original image
%   I = imread('cameraman.tif');
%   M = stdfilt(I);
%   meshCanopy(I,M,@spring)
%
%See also: mesh colormap
%

%Error Checking:
assert(nargin==2||nargin==3||nargin==4,'The number of input arguments is expected to be 2, 3 or 4.\n  It was %i',nargin);
assert(ndims(I)==2,'The first input argument, I, is required to be 2-dimensional');
assert(ndims(Mdata)==2,'The second input argument, Mdata, is required to be 2-dimensional');

%Assigning and checking the mesh colormap/height
if nargin == 2
    Cmap = [gray(256); jet(256)]; %Default
elseif ischar(Mmap)
    %String is used, assert it's right.
    valid_maps = {'jet','hsv','hot','cool','spring','summer','autumn','winter','gray','bone','copper','pink','lines'};
    assert(any(ismember(valid_maps,Mmap)),'If a string is used as a colormap, it is expected to match one of these:\n%s',sprintf('%s\n',valid_maps{:}));
    Mmap = str2func(Mmap);
    Cmap = [gray(256); Mmap(256)];
elseif isa(Mmap,'function_handle')
    %Function handle is user, assert it's right
    valid_maps = {'jet','hsv','hot','cool','spring','summer','autumn','winter','gray','bone','copper','pink','lines'};
    assert(any(ismember(valid_maps,func2str(Mmap))),'If a function_handle is used as a colormap, it is expected to match one of these:\n%s',sprintf('%s\n',valid_maps{:}));
    Cmap = [gray(256); Mmap(256)];
else
    %Explicit color map is used, make sure it's ok
    assert(size(Mmap,2)==3,'If a matrix colormap is used the second dimension must be 3');
    assert(all(Mmap(:)<=1&Mmap(:)>=0),'If a matrix colormap is used the values must all: 0 <= Mmap <= 1');
    Cmap = [gray(256); Mmap];
end
if ~exist('height','var');
    height = 80;
else
    assert(isscalar(height),'The fourth argument, height, is expected to be a scalar');
end

%Making required pieces
I = double(I); %Needs to be double for slice() and all other calculations
Mdata = double(Mdata);
[ii jj] = meshgrid(ceil(cumsum(diff([0 linspace(1,size(I,2),size(Mdata,2))]))),ceil(cumsum(diff([0 linspace(1,size(I,1),size(Mdata,1))]))));

Mdata = (Mdata - min(Mdata(:)))+height; %Scale so minimum is height (so it doesn't conflict with image)
Midx = ceil(min((length(Cmap)-256),round((length(Cmap)-255)*(Mdata-min(Mdata(:)))/(max(Mdata(:))-min(Mdata(:))))+1))+256;

if any(I(:)<0)||any(I(:)>255)
    %Scale whole image to 1:256 for the map (only if it was out of bounds before!
    I = ceil(min(256,round((255)*(I-min(I(:)))/(max(I(:))-min(I(:))))+1));
else
    %Else adjust to 1:256 integer increment
    I = ceil(I+1);
end

%Plotting7
figure;
H(1) = slice(repmat(I,[1 1 2]),[],[],1); %slice() requires at least 2x2x2
set(H(1),'EdgeColor','none') %required so image isn't just an edge
hold on
H(2) = mesh(ii,jj,double(Mdata));
hold off

%Plot Properties
axis vis3d
axis ij
axis tight
colormap(Cmap)
set(H(1),'CData',I);
set(H(2),'CData',Midx);
caxis([1 length(Cmap)])


text:
% magnitude of velocity of center
function [data,units] = compute_velmag_ctr(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  % change in center position
  dx = diff(trx(fly).x_mm,1,2);
  dy = diff(trx(fly).y_mm,1,2);

  if trx(fly).nframes < 2,
    data{i} = [];
  else
    % magnitude of velocity vector
    data{i} = sqrt(dx.^2 + dy.^2)./trx(fly).dt;
  end
end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
function [model] = em_single_instance(model_3d,test, instance, max_em_iter)
%EM_SINGLE_INSTANCE Summary of this function goes here
%   Detailed explanation goes here
params = get_params();

%% Changing test data for just the particular instance
N = size(test.points,1)/2;
testInstance.points = test.points([instance N+instance],:);
testInstance.voc_image_id{1} = test.voc_image_id{instance};
testInstance.labels = test.labels;
testInstance.bbox = test.bbox(instance,:);
testInstance.imsize = test.imsize(instance,:);
try
    testInstance.poly_x = test.poly_x(instance);
    testInstance.poly_y = test.poly_y(instance);
catch
    testInstance.mask = test.mask{instance};
end
testInstance.bbox   = test.bbox(instance,:);
if(isfield(test,'voc_rec_id'))
    testInstance.voc_rec_id(1) = test.voc_rec_id(instance);
end
if(isfield(test,'flip'))
    testInstance.flip(1) = test.flip(instance);
end

%% Obtaining the model
model = test_nrsfm(testInstance,model_3d,max_em_iter);

end

repo: CategoryShapes


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Romano Cicchetti

%params
dataset = 'thun';
households = [1,2,3,4,5,6];
evalDays_type = 'completeSM-first90';
granularity = 1;

power_pct = {};
existing = {};
for h = 1:length(households)
    house = households(h);
    path_to_evalDays = strcat(pwd, '/input/evaluation_days/', dataset, '/', evalDays_type, '/', num2str(house, '%02d'), '.mat');
    load(path_to_evalDays); % evalDays

    % total consumption
    smartmeter_consumption = read_smartmeter_data(dataset, house, evalDays, granularity, 'powerallphases');
    missing_values_idx_sm = smartmeter_consumption == -1;
    total_consumption = sum(smartmeter_consumption(~missing_values_idx_sm));

    % consumption of each appliance
    appliances = findAppliances(house, 'eco');
    for appliance = appliances
       plug_consumption = read_plug_data(dataset, house, appliance, evalDays, granularity);
       missing_values_idx_plug = plug_consumption == -1;
       sum_plug_consumption = sum(plug_consumption(~missing_values_idx_plug));
       ratio = nnz(~missing_values_idx_sm) / nnz(~missing_values_idx_plug);
       extrapolated_plug_consumption = ratio*sum_plug_consumption;
       power_pct{appliance, h} = strcat(num2str(100*extrapolated_plug_consumption / total_consumption, '%.0f'), ' \%');
       existing{appliance, h} = strcat(num2str(100*nnz(~missing_values_idx_plug) / nnz(~missing_values_idx_sm), '%.0f'), ' \%');
    end

    %power_pct{18,h} = sum(power_pct(~isnan(power_pct(:,h)),h));
end

rowLabels = {'Fridge', 'Freezer', 'Microwave', 'Dishwasher', 'Entertainment',...
        'Water kettle', 'Stove', 'Coffee machine', 'Washing machine', 'Dryer', 'Lamp', 'PC', 'Laptop', 'TV', 'Stereo', 'Tablet', 'Router'};
columnLabels = {'House 1', 'House 2', 'House 3', 'House 4', 'House 5', 'House 6'};
matrix2latex(power_pct, 'dataset.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'format', '%.0f');
matrix2latex(existing, 'dataset2.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'format', '%0.2f');

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function [secUT] = convertYearDayHrMnSec2Sec(year, day, hour, min, sec)
%convertYearDayHrMnSec2Sec Summary of this function goes here
%   Detailed explanation goes here

    [secInMin, secInHr, secInDay, secInYear] = getSecondsInVariousTimeUnits();

    secUT = (year-1)*secInYear;
    secUT = secUT + (day-1)*secInDay;
    secUT = secUT + hour*secInHr;
    secUT = secUT + min*secInMin;
    secUT = secUT + sec;
end


text:
function test_num = r8mat_exp_test_num ( )

%*****************************************************************************80
%
%% R8MAT_EXP_TEST_NUM returns the number of matrix exponential tests.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer TEST_NUM, the number of tests.
%
  test_num = 14;

  return
end


text:
function title = p39_title ( )

%*****************************************************************************80
%
%% P39_TITLE returns a title for problem 39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Bohachevsky Function #2';

  return
end


text:
function varargout = chol( f, varargin )
%CHOL    Cholesky factorization of a SEPARABLEAPPROX.
%
% R = CHOL( F ), if F is a nonnegative definite SEPARABLEAPPROX then this
% returns an upper triangular quasimatrix so that R'*R is a
% decomposition of F. If F is not nonnegative definite then an error is thrown.
%
% L = CHOL(F, 'lower'), if F is a nonnegative definite SEPARABLEAPPROX then this
% produces a lower triangular quasimatrix so that L*L' is a decomposition of F.
% If F is not nonnegative definite then an error is thrown.
%
% [R, p] = CHOL( F ), with two outputs never throwns an error message. If F is
% nonnegative definite then p is 0 and R is the same as above. If F is
% symmetric but negative definite or semidefinite then p is a positive
% integer such that R has p columns and R'*R is a rank p nonnegative definite
% SEPARABLEAPPROX that approximates F.
% This is particular useful when F is nonnegative definite, but rounding errors
% have perturbed it to be semidefinite.
%
% [L, p] = CHOL(F, 'lower') same as above but the first argument is lower
% triangular.
%
% For more information about the factorization:
% A. Townsend and L. N. Trefethen, Continuous analogues of matrix
% factorizations, Proc. Royal Soc. A., 2015.
%
% See also LU, and QR.

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

% Empty check:
if ( isempty( f ) )
    varargout = cell(1, nargout);
    return
end

% Is the separableApprox on a square domain?:
if ( ~domainCheck(f.cols, f.rows) )
    error('CHEBFUN:SEPARABLEAPPROX:chol:domain', ...
        'SEPARABLEAPPROX is not on a square domain.');
end

% Get rank of f:
k = length( f );

% All the pivots should be on the y = x line.
PivLoc = f.pivotLocations;

if ( isempty( PivLoc ) )
    % For some reason (probably because f was made with sampling data, rather
    % than a handle) f didn't have any pivot information. Make it now:
    f = compose(f, @plus, 0 );
    varargout = {chol( f )};
    return
end

% Find the first pivot location is off-diagonal:
Diagk = find( PivLoc(:,1) ~= PivLoc(:,2), 1, 'first');
k = min( k, Diagk );

% Check that the pivots are positive:
Posk = find( pivots( f ) < 0, 1, 'first');
k = min(k, Posk);

% Were all the pivots nonnegative and the locations on the diagonal?
posdef = ( k == length( f ) );

% Return an error if the function is not nonnegative definite:
if ( nargout < 2 && ~isempty( posdef ) && ~posdef )
    error('CHEBFUN:SEPARABLEAPPROX:chol:definite', ...
        'SEPARABLEAPPROX is not nonnegative definite.');
end

% Get the CDR decomposition (already computed by constructor):
[C, D, R] = cdr( f );

% Return an error if the function is not symmetric:
dom = f.domain;
r = 0.0192475;   % arbitrary point in [-1,1]
s = -.34756987;  % arbitrary point in [-1,1]
r = diff(dom(1:2))/2*r + mean(dom(1:2));
s = diff(dom(3:4))/2*s + mean(dom(3:4));
symTest = abs(feval(f,r,s) - feval(f,s,r)) < 1e2*eps;
if ( ~symTest )
    error('CHEBFUN:SEPARABLEAPPROX:chol:symmetric', ...
        'The SEPARABLEAPPROX must be a symmetric function.');
end

% How many terms are posdef:
p = k;
if ( isempty( p ) )
    p = length( f );
end
% Extract out posdef part:
C = C( :, 1:p );
D = D( 1:p, 1:p );
R = R( :, 1:p );

% Balance scaling:
C = C * sqrt( D );
R = R * sqrt( D );

% Output to user:
if ( ( nargout < 2 ) && ( nargin == 1 ) )
    varargout = { R.' };
elseif ( ( nargout < 2 ) && strcmpi( varargin{ 1 }, 'lower') )
    varargout = { C };
elseif ( ( nargout == 2 ) && ( nargin == 1 ) )
    varargout = { R.', p };
elseif ( ( nargout == 2 ) && strcmpi( varargin{ 1 }, 'lower') )
    varargout = { C, p };
end

end


text:
function obj = updateBar3h(obj, surfaceIndex)

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-SURFACE DATA STRUCTURE- %
bar_data = get(obj.State.Plot(surfaceIndex).Handle);
figure_data = get(obj.State.Figure.Handle);

%-AXIS STRUCTURE-%
axis_data = get(ancestor(bar_data.Parent,'axes'));

%-GET SCENE-%
eval(['scene = obj.layout.scene' num2str(xsource) ';']);

%-------------------------------------------------------------------------%

%-associate scene-%
obj.data{surfaceIndex}.scene = sprintf('scene%d', xsource);

%-------------------------------------------------------------------------%

%-surface type-%
obj.data{surfaceIndex}.type = 'mesh3d';

%-------------------------------------------------------------------------%

%-FORMAT DATA-%
xdata = bar_data.XData;
ydata = bar_data.ZData;
zdata = bar_data.YData;
cdata = bar_data.CData;

%-parse xedges-%
xedges = xdata(2, 1:2:end);

%-parse yedges-%
yedges = ydata(2:6:end, 2);
yedges = [yedges', mean(diff(yedges(1:2)))];

%-parse values-%
values = [];
for n = 1:6:size(zdata, 1)
    values = [values, diff(zdata(n:n+1, 2))];
end

%-parse offsets-%
offsets = zdata(1:6:end, 2)';

%-------------------------------------------------------------------------%

%-get the values to use plotly's mesh3D-%
bargap = diff(yedges(1:2)) - diff(ydata(2:3));
[X, Y, Z, I, J, K] = get_plotly_mesh3d(xedges, yedges, values, bargap);

%---------------------------------------------------------------------%

%-reformat Z according to offsets-%
m = 1;
lz2 = 0.5*length(Z);

for n = 1:4:lz2
    Z(n:n+3) = Z(n:n+3)+offsets(m);
    Z(n+lz2:n+lz2+3) = Z(n+lz2:n+lz2+3)+offsets(m);
    m = m + 1;
end

%-------------------------------------------------------------------------%

%-set mesh3d data-%
obj.data{surfaceIndex}.x = X;
obj.data{surfaceIndex}.y = Z;
obj.data{surfaceIndex}.z = Y;
obj.data{surfaceIndex}.i = int16(I-1);
obj.data{surfaceIndex}.j = int16(J-1);
obj.data{surfaceIndex}.k = int16(K-1);

%-------------------------------------------------------------------------%

%-coloring-%
cmap = figure_data.Colormap;

if isnumeric(bar_data.FaceColor)

    %-paper_bgcolor-%
    col = 255*bar_data.FaceColor;
    col = sprintf('rgb(%f,%f,%f)', col);

else
    switch bar_data.FaceColor

        case 'none'
            col = 'rgba(0,0,0,0)';

        case {'flat','interp'}

            switch bar_data.CDataMapping

                case 'scaled'
                    capCD = max(min(cdata(1,1),axis_data.CLim(2)),axis_data.CLim(1));
                    scalefactor = (capCD - axis_data.CLim(1))/diff(axis_data.CLim);
                    col =  255*(cmap(1+ floor(scalefactor*(length(cmap)-1)),:));
                case 'direct'
                    col =  255*(cmap(cdata(1,1),:));

            end

            col = sprintf('rgb(%f,%f,%f)', col);

        case 'auto'
            col = 'rgb(0,113.985,188.955)';
    end
end

obj.data{surfaceIndex}.color = col;

%-------------------------------------------------------------------------%

%-some settings-%
obj.data{surfaceIndex}.contour.show = true;
obj.data{surfaceIndex}.contour.width = 6;
obj.data{surfaceIndex}.contour.color='rgb(0,0,0)';
obj.data{surfaceIndex}.flatshading = false;

%-------------------------------------------------------------------------%

%-lighting settings-%
obj.data{surfaceIndex}.lighting.diffuse = 0.8;
obj.data{surfaceIndex}.lighting.ambient = 0.65;
obj.data{surfaceIndex}.lighting.specular = 1.42;
obj.data{surfaceIndex}.lighting.roughness = 0.52;
obj.data{surfaceIndex}.lighting.fresnel = 0.2;
obj.data{surfaceIndex}.lighting.vertexnormalsepsilon = 1e-12;
obj.data{surfaceIndex}.lighting.facenormalsepsilon = 1e-6;

obj.data{surfaceIndex}.lightposition.x = 0;
obj.data{surfaceIndex}.lightposition.y = 0;
obj.data{surfaceIndex}.lightposition.z = 0;

%-------------------------------------------------------------------------%

%-surface name-%
obj.data{surfaceIndex}.name = bar_data.DisplayName;

%-------------------------------------------------------------------------%

%-surface visible-%
obj.data{surfaceIndex}.visible = strcmp(bar_data.Visible,'on');

%-------------------------------------------------------------------------%

leg = get(bar_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{surfaceIndex}.showlegend = showleg;

%-------------------------------------------------------------------------%

%-SETTING SCENE-%

%-------------------------------------------------------------------------%

%-aspect ratio-%
ar = obj.PlotOptions.AspectRatio;

if ~isempty(ar)
    if ischar(ar)
        scene.aspectmode = ar;
    elseif isvector(ar) && length(ar) == 3
        xar = ar(1);
        yar = ar(2);
        zar = ar(3);
    end
else

    %-define as default-%
    xar = max(xedges(:));
    zar = max(yedges(:));
    yar = 0.7*max([xar, zar]);
end

scene.aspectratio.x = xar;
scene.aspectratio.y = yar;
scene.aspectratio.z = zar;

%-------------------------------------------------------------------------%

%-camera eye-%
ey = obj.PlotOptions.CameraEye;

if ~isempty(ey)
    if isvector(ey) && length(ey) == 3
        scene.camera.eye.x = ey(1);
        scene.camera.eye.y = ey(2);
        scene.camera.eye.z = ey(3);
    end
else

    %-define as default-%
    scene.camera.eye.x = xar + 7;
    scene.camera.eye.y = yar + 0;
    scene.camera.eye.z = zar + 0.5;
end

%-------------------------------------------------------------------------%

%-axis configuration-%
scene.xaxis.range = axis_data.XLim(end:-1:1);
scene.yaxis.range = axis_data.YLim;
scene.zaxis.range = axis_data.ZLim;

scene.xaxis.tickvals = axis_data.XTick;
scene.xaxis.ticktext = axis_data.XTickLabel;

scene.yaxis.tickvals = axis_data.YTick;
scene.yaxis.ticktext = axis_data.YTickLabel;

scene.zaxis.tickvals = axis_data.ZTick;
scene.zaxis.ticktext = axis_data.ZTickLabel;

scene.xaxis.zeroline = false;
scene.yaxis.zeroline = false;
scene.zaxis.zeroline = false;

scene.xaxis.showline = true;
scene.yaxis.showline = true;
scene.zaxis.showline = true;

scene.xaxis.tickcolor = 'rgba(0,0,0,1)';
scene.yaxis.tickcolor = 'rgba(0,0,0,1)';
scene.zaxis.tickcolor = 'rgba(0,0,0,1)';

scene.xaxis.ticklabelposition = 'outside';
scene.yaxis.ticklabelposition = 'outside';
scene.zaxis.ticklabelposition = 'outside';

scene.xaxis.title = axis_data.XLabel.String;
scene.yaxis.title = axis_data.YLabel.String;
scene.zaxis.title = axis_data.ZLabel.String;

%-------------------------------------------------------------------------%

%-SET SCENE TO LAYOUT-%
obj.layout = setfield(obj.layout, sprintf('scene%d', xsource), scene);

%-------------------------------------------------------------------------%

end

function bar_ = bar_data(position3d, size_)
    % position3d - 3-list or array of shape (3,) that represents the point of coords (x, y, 0), where a bar is placed
    % size = a 3-tuple whose elements are used to scale a unit cube to get a paralelipipedic bar
    % returns - an array of shape(8,3) representing the 8 vertices of  a bar at position3d

    if nargin < 2
        size_ = [1, 1, 1];
    end

    bar_ = [...
            0, 0, 0; ...
            1, 0, 0; ...
            1, 1, 0; ...
            0, 1, 0; ...
            0, 0, 1; ...
            1, 0, 1; ...
            1, 1, 1; ...
            0, 1, 1 ...
            ]; % the vertices of the unit cube

    for n =1:size(bar_, 1)
        bar_(n,:) = bar_(n,:) .* size_; % scale the cube to get the vertices of a parallelipipedic bar_
    end

    bar_ = bar_ + position3d; %translate each  bar_ on the directio OP, with P=position3d
end

function [vertices, I, J, K] = triangulate_bar_faces(positions, sizes)
    % positions - array of shape (N, 3) that contains all positions in the plane z=0, where a histogram bar is placed
    % sizes -  array of shape (N,3); each row represents the sizes to scale a unit cube to get a bar
    % returns the array of unique vertices, and the lists i, j, k to be used in instantiating the go.Mesh3d class

    if nargin < 2
        sizes = ones(size(positions,1), 3); %[(1,1,1)]*len(positions)
    else
        sizes;
        % if isinstance(sizes, (list, np.ndarray)) and len(sizes) != len(positions):
        %     raise ValueError('Your positions and sizes lists/arrays do not have the same length')
    end

    c = 1;
    for n = 1:size(positions, 1)
        if sizes(n, 3) ~= 0
            all_bars(:,:,c) = bar_data(positions(n,:), sizes(n,:))';
            c = c+1;
        end
    end

    % all_bars = [bar_data(pos, size)  for pos, size in zip(positions, sizes) if size[2]!=0]
    [r, q, p] = size(all_bars);

    % extract unique vertices from the list of all bar vertices
    all_bars = reshape(all_bars, [r, p*q])';
    [vertices, ~, ixr] = unique(all_bars, 'rows');

    %for each bar, derive the sublists of indices i, j, k assocated to its chosen  triangulation
    I = [];
    J = [];
    K = [];

    for k = 0:p-1
        aux = ixr([1+8*k, 1+8*k+2,1+8*k, 1+8*k+5,1+8*k, 1+8*k+7, 1+8*k+5, 1+8*k+2, 1+8*k+3, 1+8*k+6, 1+8*k+7, 1+8*k+5]);
        I = [ I;  aux(:)];
        aux = ixr([1+8*k+1, 1+8*k+3, 1+8*k+4, 1+8*k+1, 1+8*k+3, 1+8*k+4, 1+8*k+1, 1+8*k+6, 1+8*k+7, 1+8*k+2, 1+8*k+4, 1+8*k+6]);
        J = [ J;  aux(:)];
        aux = ixr([1+8*k+2, 1+8*k, 1+8*k+5, 1+8*k, 1+8*k+7, 1+8*k, 1+8*k+2, 1+8*k+5, 1+8*k+6, 1+8*k+3, 1+8*k+5, 1+8*k+7]);
        K = [ K;  aux(:)];
    end

end

function [X, Y, Z, I, J, K] = get_plotly_mesh3d(xedges, yedges, values, bargap)
    % x, y- array-like of shape (n,), defining the x, and y-ccordinates of data set for which we plot a 3d hist

    xsize = xedges(2)-xedges(1)-bargap;
    ysize = yedges(2)-yedges(1)-bargap;
    [xe, ye]= meshgrid(xedges(1:end-1), yedges(1:end-1));
    ze = zeros(size(xe));

    positions = zeros([size(xe), 3]);
    positions(:,:,1) = xe;
    positions(:,:,2) = ye;
    positions(:,:,3) = ze;

    [m, n, p] = size(positions);
    positions = reshape(positions, [m*n, p]);

    h = values'; h = h(:);
    sizes = [];
    for n = 1:length(h)
        sizes = [sizes; ysize, ysize, h(n)];
    end

    [vertices, I, J, K]  = triangulate_bar_faces(positions, sizes);
    X = vertices(:,1);
    Y = vertices(:,2);
    Z = vertices(:,3);

end

author: plotly


text:
function [meg_loc, meg_dir] = msi_file_MEG_loc(set_fid)
%function [meg_loc, meg_dir] = msi_file_MEG_loc(set_fid, chans)
%	MSI_FILE_GET_LONG(set_fid, chans)	Skip to a keyword in a file FID
%	and return the single float value
%

key_count = msi_file_find_keyword(set_fid, 'MSI.MegChanCount:');
if (key_count == 1)
    MegChanCount = msi_file_get_long(set_fid, 'MSI.MegChanCount:');
end

frewind(set_fid);
keyword = 'MSI.Meg_Position_Information.Begin:';

while feof(set_fid) == 0
    	line = fgetl(set_fid);
	matches = findstr(line, keyword);
	num = length(matches);
	if num > 0
	    fprintf(1, 'I found a %s\n', line);
	    for (ii=1:MegChanCount)
		line = fgetl(set_fid);
		[name, count, errmsg, nextindex] = sscanf(line, '%s\t', 1);
		[token, nextline] = strtok(line, '	');
		[bigresult, count, errmsg, nextindex] = sscanf(nextline, '\t%f\t%f\t%f\t%f\t%f\t%f', 6);
		loc(1) = bigresult(1);
		loc(2) = bigresult(2);
		loc(3) = bigresult(3);
		dir(1) = bigresult(4);
		dir(2) = bigresult(5);
		dir(3) = bigresult(6);
		%fprintf(1, 'chan: %s\n',  name);
		%fprintf(1, '  %f %f %f   ',  loc(1), loc(2), loc(3));
		%fprintf(1, '  %f %f %f \n',  dir(1), dir(2), dir(3));
		%fprintf(1, 'chan: %s (%f,%f,%f) (%f,%f,%f)\n',  name, loc[1], loc[2], loc[3], dir[1], dir[2], dir[3]);
		outloc(ii,1) = loc(1);
		outloc(ii,2) = loc(2);
		outloc(ii,3) = loc(3);

		outdir(ii,1) = dir(1);
		outdir(ii,2) = dir(2);
		outdir(ii,3) = dir(3);
	    end
meg_loc = outloc;
meg_dir = outdir;
	end
end

author: PatternRecognition

repo: OpenBMI


text:
function exact = p00_exact ( problem )

%*****************************************************************************80
%
%% P00_EXACT returns the exact integral for any problem.
%
%  Discussion:
%
%    This routine provides a "generic" interface to the exact integral
%    routines for the various problems, and allows a problem to be called
%    by index (PROBLEM) rather than by name.
%
%    In most cases, the "exact" value of the integral is not given;
%    instead a "respectable" approximation is available.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the index of the problem.
%
%    Output, real EXACT, the exact value of the integral.
%
  if ( problem == 1 )
    exact = p01_exact ( );
  elseif ( problem == 2 )
    exact = p02_exact ( );
  elseif ( problem == 3 )
    exact = p03_exact ( );
  elseif ( problem == 4 )
    exact = p04_exact ( );
  elseif ( problem == 5 )
    exact = p05_exact ( );
  elseif ( problem == 6 )
    exact = p06_exact ( );
  elseif ( problem == 7 )
    exact = p07_exact ( );
  elseif ( problem == 8 )
    exact = p08_exact ( );
  elseif ( problem == 9 )
    exact = p09_exact ( );
  elseif ( problem == 10 )
    exact = p10_exact ( );
  elseif ( problem == 11 )
    exact = p11_exact ( );
  elseif ( problem == 12 )
    exact = p12_exact ( );
  elseif ( problem == 13 )
    exact = p13_exact ( );
  elseif ( problem == 14 )
    exact = p14_exact ( );
  elseif ( problem == 15 )
    exact = p15_exact ( );
  elseif ( problem == 16 )
    exact = p16_exact ( );
  elseif ( problem == 17 )
    exact = p17_exact ( );
  elseif ( problem == 18 )
    exact = p18_exact ( );
  elseif ( problem == 19 )
    exact = p19_exact ( );
  elseif ( problem == 20 )
    exact = p20_exact ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_EXACT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', problem );
    error ( 'P00_EXACT - Fatal error!' );
  end

  return
end


text:
% [INPUT]
% txt = A string representing the text to clear.
%
% [OUTPUT]
% out = A string representing the escaped text.
%
% [NOTES]
% The function strips HTML tags and escapes backslashes.

function out = clear_text(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('txt',@(x)validateattributes(x,{'char'},{}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    txt = ipr.txt;

    nargoutchk(1,1);

    out = clear_text_internal(txt);

end

function out = clear_text_internal(txt)

    if (isempty(txt))
        out = txt;
    else
        out = regexprep(txt,'<[^>]*>','');
        out = strrep(out,filesep(),[filesep() filesep()]);
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


function res=read_cell(filename,linesep,cellsep)
if nargin < 2, linesep='\n'; end
if nargin < 3, cellsep = '\t'; end
if exist(filename,'file')
    fid = fopen(filename);
else
    % Assume that filename is either a file ide or a string
    fid = filename;
end

fileLines = textscan(fid,'%s','delimiter',linesep,'BufSize',100000);

fileLines = fileLines{1};

if regexp(fileLines{1},cellsep,'once')
    fileLines = regexprep(fileLines,['^' cellsep '|' cellsep '$'],'');
    res = regexp(fileLines,cellsep,'split');
    res = cell2matcell(res);
else
    res = fileLines;
end


text:
function x = i4vec_inc ( n, x )

%*****************************************************************************80
%
%% I4VEC_INC increments an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input/output, integer X(N), the array to be incremented.
%
  x(1:n) = x(1:n) + 1;

  return
end


text:
function bezedit(alpha, haxes, color, varargin)
    % BEZEDIT Draw and edit Bezier polynomial by dragging control points
    %
    %   BEZEDIT(ALPHA) draws the graph of a 1D Bezier polynomial or
    %   arbitrary degree. Draggable control points are also shown. As the
    %   control points are moved, the graph is updated.
    %
    %   Written by Brian G. Buss

    if nargin<2
        figure;
        haxes = gca;
        set(gca, 'XLim', [-0.2 1.2]);
    end
    if nargin<3
        color = 'b';
    end

    s = 0:0.01:1;
    line(s, bezier_vec(alpha, s), 'LineStyle', ':', 'Color', color, 'Parent', haxes);
    hline=line(s, bezier_vec(alpha, s), 'LineStyle', '-', 'Color', color, 'Parent', haxes, varargin{:});

    for k=1:length(alpha)
        hpoints(k)=impoint(haxes, [(k-1)/5 alpha(k)]);
        setPositionConstraintFcn(hpoints(k), makeConstrainToRectFcn('impoint', [(k-1)/5 (k-1)/5], [-inf inf]));
        addNewPositionCallback(hpoints(k), @(pos)update(haxes));
        setColor(hpoints(k), color);

        hcmenu = get(hpoints(k), 'uicontextmenu');
        exportmenu = uimenu(hcmenu, 'Label', 'Export', 'Callback', @(a,b)export(haxes));
        set(hpoints(k), 'uicontextmenu', hcmenu);
    end

    setappdata(haxes, 'alpha', alpha);
    setappdata(haxes, 'impoints', hpoints);
    setappdata(haxes, 'hline', hline);
end

function update(haxes)
    hpoints = getappdata(haxes, 'impoints');
    hline = getappdata(haxes, 'hline');

    alpha = zeros(1, length(hpoints));
    for k=1:length(hpoints)
        pos = getPosition(hpoints(k));
        alpha(k) = pos(2);
    end

    s = 0:0.01:1;
    set(hline, 'YData', bezier_vec(alpha, s));
end

function export(haxes)
    hpoints = getappdata(haxes, 'impoints');
    hline = getappdata(haxes, 'hline');

    alpha = zeros(1, length(hpoints));
    for k=1:length(hpoints)
        pos = getPosition(hpoints(k));
        alpha(k) = pos(2);
    end

    name=inputdlg('Enter variable name:', 'Save Bezier parameters to base workspace', 1);
    while (~isvarname(name{1}) && ~isempty(name))
        name=inputdlg(sprintf('''%s'' is not a valid variable name.\nEnter variable name:', name{1}), 'Save Bezier parameters to workspace', 1);
    end
    if isvarname(name{1})
        assignin('base',name{1},alpha);
    end
end

function y = bezier_vec(alpha, s)
    y = zeros(size(alpha,1), size(s,2));
    for k=1:length(s)
        y(:,k) = bezier(alpha, s(k));
    end
end


text:
function [g] = computeGradAtPoint(fun, x0, fAtX0, h, diffType, numPts, sparsity, useParallel)
%computeGradAtPoint Summary of this function goes here
%   Detailed explanation goes here

    switch diffType
        case FiniteDiffTypeEnum.Central
            if(mod(numPts,2) == 0) %even
                numPtsPerSide = numPts/2;
            else %odd
                numPtsPerSide = (numPts-1)/2;
            end

            if(numPtsPerSide <= 0)
                numPtsPerSide = 1;
            end

            xPts = [-numPtsPerSide:1:0, 1:1:numPtsPerSide];

        case FiniteDiffTypeEnum.Forward
            xPts = 0:1:(numPts-1);

        case FiniteDiffTypeEnum.Backward
            xPts = 0:-1:(numPts-1);

        otherwise
            error('Invalid finite difference type!  Only forward, backward, and central allowed!');
    end
    xPts = xPts(:)';

    [diffCoeff,~,~] = TT(xPts,1);

    if(useParallel)
        p = gcp('nocreate');
        if(isempty(p))
            error('Cannot run gradient in parallel: no parallel pool exists!');
        else
%             M = p.NumWorkers;
            M = parforOptions(p, 'RangePartitionMethod','fixed', 'SubrangeSize',1);
            C = parallel.pool.Constant({fun});
        end
    else
        M = 0;
        C.Value{1} = fun;
    end

    numFunOutputs = length(fAtX0);
    x0 = x0(:);
    g = nan([numel(x0),numFunOutputs]);
    zeroArr = zeros(1,size(g,1));
	parfor(i=1:size(g,1), M)
% 	for(i=1:size(g,1))
        if(isempty(sparsity) || (not(isempty(sparsity)) && sparsity(i) ~= 0))
            varArr = zeroArr;
            varArr(i) = 1;

            xDeltas = h.*(varArr(:) .* xPts); %consider FMINCON style: delta = v.*sign?(x).*max(abs(x),TypicalX); or delta = v.*max(abs(x),TypicalX);

            xToEvalAt = bsxfun(@plus, x0, xDeltas);

            numPtsToEval = size(xToEvalAt, 2);

            numerator = zeros(1,numFunOutputs);
            for(j=1:numPtsToEval) %#ok<*NO4LP>
                if(diffCoeff(j) ~= 0) %#ok<PFBNS> %otherwise we're just adding zero regardless
                    if(not(isempty(fAtX0)) && all(x0 == xToEvalAt(:,j)))
                        fAtX = fAtX0;
                    else
%                         if(useParallel)
                            f = C.Value{1}; %#ok<PFBNS>
                            fAtX = f(xToEvalAt(:,j)); %#ok<PFBNS>
%                         else
%                             fAtX = fun(xToEvalAt(:,j)); %#ok<PFBNS>
%                         end
                    end
                    numerator = numerator + diffCoeff(j) .* fAtX(:)';
                end
            end

            g(i,:) = numerator/h;
        else
            g(i,:) = zeros(1,numFunOutputs);
        end
	end
end


text:
function value = epn_lag_monomial_integral ( n, expon )

%*****************************************************************************80
%
%% EPN_LAG_MONOMIAL_INTEGRAL: integral of monomial with Laguerre weight on EPN.
%
%  Discussion:
%
%    EPN is the N-dimensional positive space [0,+oo)^N with exponential
%    or Laguerre weight function:
%
%      w(x(1:n)) = exp ( - sum ( x(1:n) ) )
%
%    value = integral ( EPN )
%      product ( 1 <= i <= n ) x(I)^expon(i) exp ( -x(i) ) dx(i)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, integer EXPON(N), the exponents.
%
%    Output, real VALUE, the value of the integral.
%
  value = 1.0;
  for i = 1 : n
    value2 = ep1_lag_monomial_integral ( expon(i) );
    value = value * value2;
  end

  return
end


text:
function pink_noise_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests RANH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  RANH is a random hold function.\n' );
  fprintf ( 1, '  Given a value U and a delay D, it returns the value\n' );
  fprintf ( 1, '  U for D calls, then resets U.\n' );

  for d = 5 : -1 : 1
    fprintf ( 1, '\n' );
    fprintf ( 1, '   I   D   Q      U           Y\n' );
    fprintf ( 1, '\n' );
    u = 0.5;
    q = 3;
    for i = 1 : 20
      [ y, u, q ] = ranh ( d, u, q );
      fprintf ( 1, '  %2d  %2d  %2d  %10f  %10f\n', i, d, q, u, y );
    end
  end

  return
end


text:
function fileExists = matRad_checkMexFileExists(filename,linkOctave)
% Checks if a matching mex file exists, and can create a link
% if a matching custom, system specific precompiled octave mex file is found
%
% call
%   matRad_checkMexFileExists(filename)
%   matRad_checkMexFileExists(filename,linkOctave)
%
% input:
%   filename:   name of the mex file (without extension)
%   linkOctave: (optional: default true) If set to true, the function will check
%               for a custom build mex file for octave with our custom extension
%               mexoct<version><system>. If such a file exists, it will create
%               a link to the file as filename.mex since octave not uses system-
%               specific extensions by default. If false, the function will only
%               check for an existing .mex file for octave.
% output:
%   fileExists: true if the mex file exists (or can be linked), and false
%               otherwise
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2020 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matRad_cfg = MatRad_Config.instance();

if nargin < 2
    linkOctave = true;
end

%Explicitly check for matching mex file (id 3)
fileExists = (exist(filename,'file') == 3);

%For octave we have experimental precompiled files for Octave 5 64 bit
[env,ver] = matRad_getEnvironment();

if ~fileExists && matRad_cfg.isOctave && linkOctave

    %versionstring
    verStripped = erase(ver,'.');

    %Check Architecture
    [~,maxArraySize] = computer();

    if maxArraySize > 2^31
        bitExt = '64';
    else
        bitExt = '32';
    end

    systemext = ['mexoct' verStripped];
    if ispc
        systemext = [systemext 'w'];
    elseif ismac
        systemext = [systemext 'mac'];
    elseif isunix
        systemext = [systemext 'a'];
    else
        %No file for unknown operating system
        fileExists = false;
        return;
    end

    %Build our octfilename
    octfilename = [filename '.' systemext bitExt];

    %Check if matching precompiled octave file exists
    fileExists = (exist(octfilename,'file') == 2);

    %If it exists, we create a link with the ending "mex" which is used by octave
    if fileExists
        mexFileName = [filename '.' mexext];

        %Make the link in the right directory
        cFilePath = which(octfilename);
        [MexFolder,~,~] = fileparts(cFilePath);

        oldpath = fullfile(MexFolder,octfilename);
        linkpath = fullfile(MexFolder,mexFileName);

        %Create link

        [status,msg] = link(oldpath,linkpath);

        %We need this so the file directory cache is updated to recognize the linked file
        rehash;

        if status == 0 && exist([filename '.mex']) == 3
            matRad_cfg.dispWarning('Trying to use a precompiled mex for Octave %s. This is experimental!',ver);
            fileExists = true;
        else
            matRad_cfg.dispWarning('Could not link existing precompiled mex file %s to %s! Reason: %s',octfilename,mexFileName, msg);
            fileExists = false;
        end
    end
end

repo: matRad


text:
%%*****************************************************************
%% HSDsqlpCpert: perturb C.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

    function [At,Cpert] = HSDsqlpCpert(blk,At,par,C,X,Cpert,runhist);

    iter = length(runhist.pinfeas);
    prim_infeas = runhist.pinfeas(iter);
    dual_infeas = runhist.dinfeas(iter);
    relgap      = runhist.relgap(iter);
    infeas      = runhist.infeas(iter);
    theta       = runhist.theta(iter);
%%
    Cpertold = Cpert;
    err = max(relgap,infeas);
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       n = sum(pblk{2});
       tmp = max(1,norm(C{p},'fro'))/sqrt(n);
       if (err < 1e-6)
          if (norm(X{p},'fro') < 1e2); const=0.2; else; const=0.3; end
          Cpert(p) = max(const*Cpert(p),1e-10*tmp);
       elseif (err < 1e-2)
          if (norm(X{p},'fro') < 1e2); const=0.4; else; const=0.5; end
          Cpert(p) = max(const*Cpert(p),1e-8*tmp);
       else
	  Cpert(p) = max(0.9*Cpert(p),1e-6*tmp);
       end
       Cpert = min(Cpert,Cpertold);
       if (prim_infeas < min([0.1*dual_infeas, 1e-7*runhist.pinfeas(1)])) ...
      	  & (iter > 1 & dual_infeas > 0.8*runhist.dinfeas(iter-1) & relgap < 1e-4)
          Cpert(p) = 0.5*Cpert(p);
       elseif (dual_infeas < min([0.1*prim_infeas, 1e-7*runhist.dinfeas(1)])) ...
          & (iter > 1 & prim_infeas > 0.8*runhist.pinfeas(iter-1) & relgap < 1e-4)
          Cpert(p) = 0.5*Cpert(p);
       elseif (max(relgap,1e-2*infeas) < 1e-6 & relgap < 0.1*infeas)
          Cpert(p) = 0.5*Cpert(p);
       end
       if (prim_infeas < min([1e-4*dual_infeas,1e-7]) & theta < 1e-6) ...
	  | (prim_infeas < 1e-4 & theta < 1e-10)
          Cpert(p) = 0.1*Cpert(p);
       elseif (dual_infeas < min([1e-4*prim_infeas,1e-7]) & theta < 1e-6) ...
	  | (dual_infeas < 1e-4 & theta < 1e-10)
          Cpert(p) = 0.1*Cpert(p);
       elseif (iter > 1 & theta > 0.9*runhist.theta(iter-1) & infeas < 1e-3)
          Cpert(p) = 0.1*Cpert(p);
       end
       if strcmp(pblk{1},'s')
          Cnew = C{p} + Cpert(p)*speye(n);
          At{p}(:,par.invpermA(p,end-1)) = -svec(pblk,Cnew,1);
       else
          Cnew = C{p} + Cpert(p)*ones(n,1);
          At{p}(:,par.invpermA(p,end-1)) = -Cnew;
       end
    end
%%*****************************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%% Copyright (C) 2014-2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym symprod (@var{f}, @var{n}, @var{a}, @var{b})
%% @defmethodx @@sym symprod (@var{f}, @var{n}, [@var{a} @var{b}])
%% @defmethodx @@sym symprod (@var{f}, @var{a}, @var{b})
%% @defmethodx @@sym symprod (@var{f}, [@var{a} @var{b}])
%% @defmethodx @@sym symprod (@var{f}, @var{n})
%% @defmethodx @@sym symprod (@var{f})
%% Symbolic product.
%%
%% The product of the expression @var{f} as variable @var{n} changes
%% from @var{a} to @var{b}.  When @var{n} is omitted it is determined
%% using @code{symvar} and defaults to @code{x} if @var{f} is
%% constant. The limits @var{a} and @var{b} default to @code{1} and
%% @var{n} respectively.
%%
%% Examples:
%% @example
%% @group
%% syms n m x
%% symprod(sin(n*x), n, [1 3])
%%   @result{} (sym) sin(x)⋅sin(2⋅x)⋅sin(3⋅x)
%% symprod(n, n, 1, m)
%%   @result{} (sym) m!
%% @end group
%% @end example
%%
%% Unevaluated product:
%% @example
%% @group
%% syms x m
%% symprod(sin(x), x, [1 m])
%%   @result{} (sym)
%%         m
%%       ─┬─┬─
%%        │ │ sin(x)
%%        │ │
%%       x = 1
%% @end group
%% @end example
%%
%% @seealso{@@sym/symsum, @@sym/prod}
%% @end defmethod

function S = symprod(f, n, a, b)

  if (nargin > 4)
    print_usage ();
  end

  idx1.type = '()';
  idx1.subs = {1};
  idx2.type = '()';
  idx2.subs = {2};

  if (nargin == 1)
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
    a = sym(1);
    b = n;
  elseif (nargin == 2) && (length(n) == 2)
    f = sym(f);
    %a = n(1);  % issue #17
    %b = n(2);
    a = subsref(n, idx1);
    b = subsref(n, idx2);
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
  elseif (nargin == 2)
    f = sym(f);
    n = sym(n);
    a = sym(1);
    b = n;
  elseif (nargin == 3) && (length(a) == 2)
    f = sym(f);
    n = sym(n);
    %b = a(2);  % issue #17
    %a = a(1);
    b = subsref(a, idx2);
    a = subsref(a, idx1);
  elseif (nargin == 3)
    f = sym(f);
    b = a;
    a = n;
    n = symvar(f, 1);
    if (isempty(n))
      n = sym('x');
    end
  else
    f = sym(f);
    n = sym(n);
    a = sym(a);
    b = sym(b);
  end

  cmd = { '(f, n, a, b) = _ins'
          'S = sp.product(f, (n, a, b))'
          'return S,' };

  S = pycall_sympy__ (cmd, sym(f), sym(n), sym(a), sym(b));

end

%!error symprod (sym(1), 2, 3, 4, 5)

%!test
%! % simple
%! syms n
%! assert (isequal (symprod(n, n, 1, 10), factorial(sym(10))))
%! assert (isequal (symprod(n, n, sym(1), sym(10)), factorial(10)))

%!test
%! % one input
%! syms n
%! f = symprod (n);
%! g = factorial (n);
%! assert (isequal (f, g))
%! f = symprod (2*n);
%! g = 2^n * factorial (n);
%! assert (isequal (f, g))

%!test
%! % constant input
%! f = symprod (sym(2));
%! syms x
%! g = 2^x;
%! assert (isequal (f, g))

%!test
%! % two inputs
%! syms n
%! f = symprod (2*n, n);
%! g = 2^n * factorial (n);
%! assert (isequal (f, g))

%!test
%! % two inputs, second is range
%! syms n
%! f = symprod (n, [1 6]);
%! g = 720;
%! assert (isequal (f, g))
%! f = symprod (n, [sym(1) 6]);
%! g = 720;
%! assert (isequal (f, g))
%! f = symprod (2*n, [1 6]);
%! g = sym(2)^6*720;
%! assert (isequal (f, g))

%!test
%! % three inputs, last is range
%! syms n
%! f = symprod (2*n, n, [1 4]);
%! g = sym(384);
%! assert (isequal (f, g))
%! f = symprod (2*n, n, [sym(1) 4]);
%! g = sym(384);
%! assert (isequal (f, g))
%! f = symprod (2, n, [sym(1) 4]);
%! g = sym(16);
%! assert (isequal (f, g))

%!test
%! % three inputs, no range
%! syms n
%! f = symprod (2*n, 1, 4);
%! g = sym(384);
%! assert (isequal (f, g))
%! f = symprod (5, sym(1), 3);
%! g = sym(125);
%! assert (isequal (f, g))

%!test
%! % infinite product
%! syms a n oo
%! zoo = sym('zoo');
%! assert (isequal (symprod(a, n, 1, oo), a^oo))
%! assert (isequal (symprod(a, n, 1, inf), a^oo))

%%!test
%%! % FIXME: commented out test...
%%! % SymPy 0.7.6: nan
%%! % SymPy git: interesting that 1**oo is nan but this is still 1
%%! assert (isequal (symprod(1, n, 1, oo), sym(1)))

repo: octsympy


% tracker_ensemble: Correlation filter tracking with convolutional features
%
% Input:
%   - video_path:          path to the image sequence
%   - img_files:           list of image names
%   - pos:                 intialized center position of the target in (row, col)
%   - target_sz:           intialized target size in (Height, Width)
% 	- padding:             padding parameter for the search area
%   - lambda:              regularization term for ridge regression
%   - output_sigma_factor: spatial bandwidth for the Gaussian label
%   - interp_factor:       learning rate for model update
%   - cell_size:           spatial quantization level
%   - show_visualization:  set to True for showing intermediate results
% Output:
%   - positions:           predicted target position at each frame
%   - time:                time spent for tracking
%
%   It is provided for educational/researrch purpose only.
%   If you find the software useful, please consider cite our paper.
%
%   Hierarchical Convolutional Features for Visual Tracking
%   Chao Ma, Jia-Bin Huang, Xiaokang Yang, and Ming-Hsuan Yang
%   IEEE International Conference on Computer Vision, ICCV 2015
%
% Contact:
%   Chao Ma (chaoma99@gmail.com), or
%   Jia-Bin Huang (jbhuang1@illinois.edu).

function [positions, time] = tracker_ensemble(video_path, img_files, pos, target_sz, ...
    padding, lambda, output_sigma_factor, interp_factor, cell_size, show_visualization)

% ================================================================================
% Environment setting
% ================================================================================

indLayers = [37, 28, 19];   % The CNN layers Conv5-4, Conv4-4, and Conv3-4 in VGG Net
nweights  = [1, 0.5, 0.25]; % Weights for combining correlation filter responses
numLayers = length(indLayers);

% Get image size and search window size
im_sz     = size(imread([video_path img_files{1}]));
window_sz = get_search_window(target_sz, im_sz, padding);

% Compute the sigma for the Gaussian function label
output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;

%create regression labels, gaussian shaped, with a bandwidth
%proportional to target size    d=bsxfun(@times,c,[1 2]);

l1_patch_num = floor(window_sz/ cell_size);

% Pre-compute the Fourier Transform of the Gaussian function label
yf = fft2(gaussian_shaped_labels(output_sigma, l1_patch_num));

% Pre-compute and cache the cosine window (for avoiding boundary discontinuity)
cos_window = hann(size(yf,1)) * hann(size(yf,2))';

% Create video interface for visualization
if(show_visualization)
    update_visualization = show_video(img_files, video_path);
end

% Initialize variables for calculating FPS and distance precision
time      = 0;
rects = zeros(numel(img_files), 4);
nweights  = reshape(nweights,1,1,[]);

% Note: variables ending with 'f' are in the Fourier domain.
model_xf     = cell(1, numLayers);
model_alphaf = cell(1, numLayers);

current_scale_factor=1;

% ================================================================================
% Start tracking
% ================================================================================
for frame = 1:numel(img_files),
    im = imread([video_path img_files{frame}]); % Load the image at the current frame
    if ismatrix(im)
        im = cat(3, im, im, im);
    end

    tic();
    % ================================================================================
    % Predicting the object position from the learned object model
    % ================================================================================
    if frame > 1
        % Extracting hierarchical convolutional features
        feat = extractFeature(im, pos, window_sz, cos_window, indLayers);
        % Predict position
        pos  = predictPosition(feat, pos, indLayers, nweights, cell_size, l1_patch_num, ...
            model_xf, model_alphaf);

        % Scale estimation
        current_scale_factor = estimate_scale( rgb2gray(im), pos, current_scale_factor);
    else
        init_scale_para(rgb2gray(im), target_sz, pos);
    end

    % ================================================================================
    % Learning correlation filters over hierarchical convolutional features
    % ================================================================================
    % Extracting hierarchical convolutional features
    feat  = extractFeature(im, pos, window_sz, cos_window, indLayers);
    % Model update
    [model_xf, model_alphaf] = updateModel(feat, yf, interp_factor, lambda, frame, ...
        model_xf, model_alphaf);

    % ================================================================================
    % Save predicted position and timing
    % ================================================================================
    % positions(frame,:) = pos;

    target_sz_t=target_sz*current_scale_factor;
    box = [pos([2,1]) - target_sz_t([2,1])/2, target_sz_t([2,1])];
    rects(frame,:)=box;

    time = time + toc();

    % Visualization
    if show_visualization,
        %box = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])];
        stop = update_visualization(frame, box);
        if stop, break, end  %user pressed Esc, stop early
        drawnow
        % 			pause(0.05)  % uncomment to run slower
    end
end

end

function pos = predictPosition(feat, pos, indLayers, nweights, cell_size, l1_patch_num, ...
    model_xf, model_alphaf)

% ================================================================================
% Compute correlation filter responses at each layer
% ================================================================================
res_layer = zeros([l1_patch_num, length(indLayers)]);

for ii = 1 : length(indLayers)
    zf = fft2(feat{ii});
    kzf=sum(zf .* conj(model_xf{ii}), 3) / numel(zf);

    temp= real(fftshift(ifft2(model_alphaf{ii} .* kzf)));  %equation for fast detection
    res_layer(:,:,ii)=temp/max(temp(:));
end

% Combine responses from multiple layers (see Eqn. 5)
response = sum(bsxfun(@times, res_layer, nweights), 3);

% ================================================================================
% Find target location
% ================================================================================
% Target location is at the maximum response. we must take into
% account the fact that, if the target doesn't move, the peak
% will appear at the top-left corner, not at the center (this is
% discussed in the KCF paper). The responses wrap around cyclically.
[vert_delta, horiz_delta] = find(response == max(response(:)), 1);
vert_delta  = vert_delta  - floor(size(zf,1)/2);
horiz_delta = horiz_delta - floor(size(zf,2)/2);

% Map the position to the image space
pos = pos + cell_size * [vert_delta - 1, horiz_delta - 1];

end

function [model_xf, model_alphaf] = updateModel(feat, yf, interp_factor, lambda, frame, ...
    model_xf, model_alphaf)

numLayers = length(feat);

% ================================================================================
% Initialization
% ================================================================================
xf       = cell(1, numLayers);
alphaf   = cell(1, numLayers);

% ================================================================================
% Model update
% ================================================================================
for ii=1 : numLayers
    xf{ii} = fft2(feat{ii});
    kf = sum(xf{ii} .* conj(xf{ii}), 3) / numel(xf{ii});
    alphaf{ii} = yf./ (kf+ lambda);   % Fast training
end

% Model initialization or update
if frame == 1,  % First frame, train with a single image
    for ii=1:numLayers
        model_alphaf{ii} = alphaf{ii};
        model_xf{ii} = xf{ii};
    end
else
    % Online model update using learning rate interp_factor
    for ii=1:numLayers
        model_alphaf{ii} = (1 - interp_factor) * model_alphaf{ii} + interp_factor * alphaf{ii};
        model_xf{ii}     = (1 - interp_factor) * model_xf{ii}     + interp_factor * xf{ii};
    end
end

end

function feat  = extractFeature(im, pos, window_sz, cos_window, indLayers)

% Get the search window from previous detection
patch = get_subwindow(im, pos, window_sz);
% Extracting hierarchical convolutional features
feat  = get_features(patch, cos_window, indLayers);

end


text:
% Run experiments for simulated datasets.
% Usage: Q = metoffice_experiment_pca('hadsst2d1', 1, 40, 4)

function [ Q, perfmeas ] = metoffice_experiment_pca(dataset, anomalies, D, remove_val)

%
% Load data
%

if nargin < 4 || isempty(remove_val)
    remove_val = 1;
end
[data,folder,maskfile] = metoffice_getdata(dataset, anomalies, remove_val);

% Form the data matrix
Y = data.data;
[M,N] = size(Y);
Obs = ~isnan(Y);

%
% PCA inference
%

% Number of components
if nargin < 3 || isempty(D)
    D = 120;
end

% Filename for saving the results
%folder = [folder '/pca'];
mkdir(folder);
filename = sprintf('%s/%s/pca_D=%d_anomalies=%d_remval=%02d_%s', ...
                   folder, ...
                   dataset, ...
                   D, ...
                   anomalies, ...
                   remove_val, ...
                   datestr(now,'yyyymmdd'))

% PCA module for X (one constant component for modeling bias)
prior.mu = [1; zeros(D-1,1)];
prior.CovX = diag([1e-6; ones(D-1,1)]);
X_module = factor_module_iid('prior', prior);

% ARD module for W
W_module = factor_module_ard();

% Isotropic noise (precisions weighted proportionally to grid size)
[LON,LAT] = meshgrid(data.lon, data.lat);
weights = cosd(LAT(:));
weights = metoffice_remove_bins(weights,maskfile);
weights = repmat(weights, [1, N]);
noise_module = noise_module_isotropic('init', struct('tau', 10), ...
                                      'weights', weights);

% Run VB PCA
Q = vbfa(D, Y, W_module, X_module, noise_module, ...
         'maxiter', 2, ...
         'rotate', true, ...
         'autosavefile', filename, ...
         'autosave', [1 20:20:2000]);

if nargout == 2
    % Reconstruct
    Yrec = Q.W'*Q.X;

    % Remove the autosave file
    delete( [ filename '.mat' ] )

    addpath /home/alexilin/matlab/metoffice
    Yrec = add_land( Yrec, dataset );
    Yrec = add_climatology( Yrec, dataset, anomalies );
    perfmeas = compute_rmse( Yrec, dataset, remove_val );

else
    % Reconstruct
    Yrec = Q.W'*Q.X;

    % Some performance measures
    fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
            rmsew(Y(Obs)-Yrec(Obs),weights(Obs)));

end

% Save the results
%save(filename, '-struct', 'Q');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [var, cut, children, catsplit] = tree_getParameters(dt)

if iscell(dt.cut)
    var = dt.var;
    children = dt.children;

    cut = zeros(size(dt.cut));
    catind = var<0;
    ncat = sum(catind);
    catsplit = cell(ncat, 2);
    ncat = 0;
    for k = 1:numel(cut)
        if var(k)>=0
            cut(k) = dt.cut{k};
        else
            ncat = ncat+1;
            cut(k) = ncat;
            catsplit(ncat, :) = dt.cut{k};
        end
    end
else
    var = dt.var;
    cut = dt.cut;
    children = dt.children;
    catsplit = dt.catsplit;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm