text:
function [InitialObservation,StockSaved] = myResetFunction(trainData)
% Copyright 2020 The MathWorks, Inc.

%% init
[total_step,n_step] = size(trainData);
init_invest = 20000;

%% reset the environment
cur_step = 1;
stock_price = trainData(cur_step, :);

%randomize stock holdings at the beginning
Ir = rand(1);

%% Init cash in and cash invested
cash_in_hand = round(init_invest*Ir);
invest_cash = init_invest-cash_in_hand;

stock_init = rand(1,3);
stock_init = stock_init/sum(stock_init);
stock_init_buy = invest_cash.*stock_init;

%% Save it in shared structure
StockSaved.stock_owned  = round(stock_init_buy./stock_price{:,:});
StockSaved.cash_in_hand = cash_in_hand;

%% Init Indicators
Ind1 = zeros(1,3);
Ind2 = Ind1;
Ind3 = Ind1;
Ind4 = Ind1;

%% Init price differentiation ratio
diffPricesBought = zeros(1,3);
%% Init first state and buffer
StockSaved.State = [StockSaved.stock_owned,diffPricesBought,StockSaved.cash_in_hand,Ind1,Ind2,Ind3,Ind4];
StockSaved.last7 = zeros(7,3);
StockSaved.last7(1,:) = stock_price{:,:};
StockSaved.prevBoughtPrices = Ind1;
StockSaved.profits = [];
%% Adding parameters to the structure saved
StockSaved.total_step = total_step;
StockSaved.n_step = n_step;
StockSaved.cur_step = cur_step;
StockSaved.cur_val = [];

%% Updating initial Observation
InitialObservation = StockSaved.State;

author: matlab-deep-learning


text:
function distance = getStructureDistance(structNum1,structNum2,planC)
%distance = getStructureDistance.m(structNum1,structNum1,planC)
%
%This function returns the distance(cm) between center of masses of
%structNum1 and structNum1
%
%APA, 10/08/2009
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Check if plan passed, if not use global.
if ~exist('planC')
    global planC;
end
indexS = planC{end};

%Compute centroid of structNum1
[x1,y1,z1] = calcIsocenter(structNum1, 'COM', planC);

%Compute centroid of structNum2
[x2,y2,z2] = calcIsocenter(structNum2, 'COM', planC);

%Compute distance between centroid1 and centroid2
distance = sqrt((x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Convert feature prediction results fot deep image reconstruction
%
% Output files:
%
% decodedfeatures -+- imagery --- matconvnet -+- cnn1 -+- Subject1 -+- FFA --- Mat files
%                  |                          |        |            |
%                  |                          |        |            +- LOC --- Mat files
%                  |                          |        |            |
%                  |                          |        |           ...
%                  |                          |        |
%                  |                          |        +- Subject2 --- ...
%                  |                          |        |
%                  |                          |       ...
%                  |                          |
%                  |                          +- cnn2 --- ...
%                  |                          |
%                  |                         ...
%                  |
%                  +- perception --- ...
%

clear all;

%% Initial settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Data settings
% subjectList  : List of subject IDs [cell array]
% roiList      : List of ROIs [cell array]
% featureList  : List of image features [cell array]

subjectList  = {'Subject1', 'Subject2', 'Subject3', 'Subject4', 'Subject5'};
roiList      = {'V1', 'V2', 'V3', 'V4', 'FFA', 'LOC', 'PPA', 'LVC', 'HVC',  'VC'};
featureList  = {'cnn1', 'cnn2', 'cnn3', 'cnn4', 'cnn5', 'cnn6', 'cnn7', 'cnn8'};

network = 'matconvnet';

%% Directory settings
workDir = pwd;
resultsDir = fullfile(workDir, 'results'); % Directory to save analysis results
outputDir = fullfile(workDir, 'decodedfeatures');  % Directory to save converted decoded features

%% File name settings
resultFileNameFormat = @(s, r, f) fullfile(resultsDir, sprintf('%s/%s/%s.mat', s, r, f));

%% Main %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for isub = 1:length(subjectList)
for iroi = 1:length(roiList)
for ifeat = 1:length(featureList)
    subject = subjectList{isub};
    roi = roiList{iroi};
    feature = featureList{ifeat};

    resultfile = resultFileNameFormat(subject, roi, feature);

    dirDecodedFeatruePercept = fullfile(outputDir, 'perception', network, ...
                                        featureList{ifeat}, subjectList{isub}, roiList{iroi});
    dirDecodedFeatrueImagery = fullfile(outputDir, 'imagery', network, ...
                                        featureList{ifeat}, subjectList{isub}, roiList{iroi});

    if ~exist(resultfile)
        error(sprintf('Result file not found: %s\n', resultfile));
    end

    res = load(resultfile);

    categoryTestPercept = res.categoryTestPercept;
    categoryTestImagery = res.categoryTestImagery;
    predfeatPerceptAve = res.predictPerceptAveraged;
    predfeatImageryAve = res.predictImageryAveraged;

    %% Save decoded features for image reconstruction
    setupdir(dirDecodedFeatruePercept);
    setupdir(dirDecodedFeatrueImagery);

    subsPercept.type = '()';
    subsPercept.subs = repmat({':'}, 1, ndims(predfeatPerceptAve));

    subsImagery.type = '()';
    subsImagery.subs = repmat({':'}, 1, ndims(predfeatImageryAve));

    for p = 1:length(categoryTestPercept)
        catid = categoryTestPercept(p);
        subsPercept.subs{1} = p;

        feat = squeeze(shiftdim(subsref(predfeatPerceptAve, subsPercept), 1));
        resfilePercept = fullfile(dirDecodedFeatruePercept, ...
                                  sprintf('%s-%s-%s-%s-%s-%s-%d.mat', ...
                                          'decodedfeatures', 'perception', network, ...
                                          featureList{ifeat}, subjectList{isub}, roiList{iroi}, catid));
        save(resfilePercept, 'feat', '-v7.3');
    end
    for p = 1:length(categoryTestImagery)
        catid = categoryTestImagery(p);
        subsImagery.subs{1} = p;

        feat = squeeze(shiftdim(subsref(predfeatImageryAve, subsImagery), 1));
        resfileImagery = fullfile(dirDecodedFeatrueImagery, ...
                                  sprintf('%s-%s-%s-%s-%s-%s-%d.mat', ...
                                          'decodedfeatures', 'imagery', network, ...
                                          featureList{ifeat}, subjectList{isub}, roiList{iroi}, catid));
        save(resfileImagery, 'feat', '-v7.3');
    end

end
end
end

author: KamitaniLab

repo: GenericObjectDecoding


text:
function [x,fval,exitflag,info,lambda,Opt] = opti_quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,opts)
%OPTI_QUADPROG Solve a QP using an OPTI QP Solver (Matlab Overload)
%
%   [x,fval,exitflag,info] = opti_quadprog(H,f,A,b,Aeq,beq,lb,ub,x0) solves
%   the quadratic program min 1/2x'Hx + f'x where A,b are the inequality
%   constraints, Aeq,beq are the equality constraints and lb,ub are the
%   bounds.
%
%   [x,fval,exitflag,info] = opti_quadprog(H,...,ub,opts) allows the user
%   to specify optiset options. This includes specifying a solver via the
%   'solver' field of optiset.
%
%   [x,...,info,lambda] = opti_quadprog(H,...) returns a structure of the
%   Lagrange multipliers (dual solution).
%
%   [x,...,lambda,Opt] = opti_quadprog(H,...) returns the internally built
%   OPTI object.

%   Copyright (C) 2011 Jonathan Currie (IPL)

% Handle missing arguments
if nargin < 10, opts = optiset; end
if nargin < 9, x0 = []; end
if nargin < 8, ub = []; end
if nargin < 7, lb = []; end
if nargin < 6, beq = []; end
if nargin < 5, Aeq = []; end
if nargin < 4, error('You must supply at least 4 arguments to opti_quadprog'); end

%Build OPTI Object
Opt = opti('qp',H,f,'ineq',A,b,'eq',Aeq,beq,'bounds',lb,ub,'x0',x0,'options',opts);

%Solve
[x,fval,exitflag,info] = solve(Opt);

%Extra Dual Solution
if(isfield(info,'Lambda'))
    lambda = info.Lambda;
else
    lambda = [];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% //This software is licensed under the BSD 3 Clause license: http://opensource.org/licenses/BSD-3-Clause
%
%
% //Copyright (c) 2013, University of Oxford
% //All rights reserved.
%
% //Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
%
% //Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
% //Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% //Neither the name of the University of Oxford nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
% //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%   The method implemented in this file has been patented by their original
%   authors. Commercial use of this code is thus strongly not
%   recocomended.
%
% //Authors: 	Gari D Clifford -
% //            Roberta Colloca -
% //			Julien Oster	-

function AFEv = comput_AFEv2(segment)

%comput_AFEv takes the RR intervals segment (column vector)
%and computes AFEv

%Inputs
%segment:   RR intervals segment (column vector)

%Output
%AFEv:   feature for detection of Atrial Fibrillation

%Compute dRR intervals series
dRR=comp_dRR(segment);

%Compute metrics
[OriginCount,IrrEv,PACEv]=metrics(dRR);

%Compute AFEvidence
AFEv=IrrEv-OriginCount-2*PACEv;

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function cvx_optval = sum_largest( x, k, dim )

%SUM_LARGEST Sum of the largest k values of a vector.
%   For a real vector X and an integer k between 1 and length(X) inclusive,
%   y = SUM_LARGEST(X,k) is the sum of the k largest elements of X; e.g.,
%       temp = sort( x )
%       y = sum( temp( 1 : k ) )
%   If k=1, then SUM_LARGEST(X,k) is equivalent to MAX(X); if k=length(X),
%   then SUM_LARGEST(X,k) is equivalent to SUM(X).
%
%   Both X and k must be real, and k must be a scalar. But k is not, in
%   fact, constrained to be an integer between 1 and length(X); the
%   function is extended continuously and logically to all real k. For
%   example, if k <= 0, then SUM_LARGEST(X,k)=0. If k > length(X), then
%   SUM_LARGEST(X,k)=SUM(X). Non-integer values of k interpolate linearly
%   between their integral neighbors.
%
%   For matrices, SUM_LARGEST(X,k) is a row vector containing the
%   application of SUM_LARGEST to each column. For N-D arrays, the
%   SUM_LARGEST operation is applied to the first non-singleton dimension
%   of X.
%
%   SUM_LARGEST(X,k,DIM) performs the operation along dimension DIM of X.
%
%   Disciplined convex programming information:
%       SUM_LARGEST(X,...) is convex and nondecreasing in X. Thus, when
%       used in CVX expressions, X must be convex (or affine). k and DIM
%       must both be constant.

%
% Check arguments
%

narginchk(2,3);
if ~isreal( x ),
    error( 'First argument must be real.' );
elseif ~isnumeric( k ) || ~isreal( k ) || length( k ) ~= 1,
    error( 'Second argument must be a real scalar.' );
elseif nargin < 3 || isempty( dim ),
    dim = cvx_default_dimension( size( x ) );
elseif ~cvx_check_dimension( dim, false ),
    error( 'Third argument, if supplied, must be a positive integer.' );
end

%
% Determine output size
%

sx = size( x );
nd = max( dim, length( sx ) );
sx = [ sx, ones( 1, dim - nd ) ];
sy = sx;
sy( dim ) = 1;

%
% Compute results
%

if k <= 0,

    cvx_optval = zeros( sy );

elseif k <= 1,

    cvx_optval = k * max( x, [], dim );

elseif k >= sx( dim ),

    cvx_optval = sum( x, dim );

else

    ck = ceil( k );
    x = sort( x, dim );
    ndxs = cell( 1, nd );
    [ ndxs{:} ] = deal( ':' );
    ndxs{ dim } = size( x, dim ) - ( 0 : ck - 1 );
    x = x( ndxs{ : } );
    if k ~= ck,
        ndxs{ dim } = ck;
        x( ndxs{ : } ) = ( k - floor( k ) ) * x( ndxs{ : } );
    end
    cvx_optval = sum( x, dim );

end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% CRTP.m - Create an initial population
%
% This function creates a binary population of given size and structure.
%
% Syntax: [Chrom Lind BaseV] = crtbp(Nind, Lind, Base)
%
% Input Parameters:
%
%		Nind	- Either a scalar containing the number of individuals
%			  in the new population or a row vector of length two
%			  containing the number of individuals and their length.
%
%		Lind	- A scalar containing the length of the individual
%			  chromosomes.
%
%		Base	- A scalar containing the base of the chromosome
%			  elements or a row vector containing the base(s)
%			  of the loci of the chromosomes.
%
% Output Parameters:
%
%		Chrom	- A matrix containing the random valued chromosomes
%			  row wise.
%
%		Lind	- A scalar containing the length of the chromosome.
%
%		BaseV	- A row vector containing the base of the
%			  chromosome loci.

% Author: Andrew Chipperfield
% Date:	19-Jan-94

function [Chrom, Lind, BaseV] = crtbp(Nind, Lind, Base)
nargs = nargin ;

% Check parameter consistency

if nargs >= 1, [mN, nN] = size(Nind) ; end
if nargs >= 2, [mL, nL] = size(Lind) ; end
if nargs == 3, [mB, nB] = size(Base) ; end

if nN == 2
   if (nargs == 1)
      Lind = Nind(2) ; Nind = Nind(1) ; BaseV = crtbase(Lind) ;
   elseif (nargs == 2 & nL == 1)
      BaseV = crtbase(Nind(2),Lind) ; Lind = Nind(2) ; Nind = Nind(1) ;
   elseif (nargs == 2 & nL > 1)
      if Lind ~= length(Lind), error('Lind and Base disagree'); end
      BaseV = Lind ; Lind = Nind(2) ; Nind = Nind(1) ;
   end
elseif nN == 1
   if nargs == 2
      if nL == 1, BaseV = crtbase(Lind) ;
      else, BaseV = Lind ; Lind = nL ; end
   elseif nargs == 3
      if nB == 1, BaseV = crtbase(Lind,Base) ;
      elseif nB ~= Lind, error('Lind and Base disagree') ;
      else BaseV = Base ; end
   end
else
   error('Input parameters inconsistent') ;
end

% Create a structure of random chromosomes in row wise order, dimensions
% Nind by Lind. The base of each chromosomes loci is given by the value
% of the corresponding element of the row vector base.

Chrom = floor(rand(Nind,Lind).*BaseV(ones(Nind,1),:)) ;

% End of file

author: vonsylvia


text:
function DEM = reclassify(DEM,varargin)

%RECLASSIFY generate univariate class intervals for an instance of GRIDobj
%
% Syntax
%
%     C = reclassify(DEM);
%     C = reclassify(DEM,'method',value)
%
% Description
%
%     reclassify bins continous values of an instance of GRIDobj by setting
%     class intervals based on different classification methods. The
%     default method is equal interval classification with 10 classes.
%
% Input Arguments
%
%     DEM       Grid (class = GRIDobj)
%
%     'Methods' and [values]
%
%     'equalintervals'   [number of classes]
%     'definedintervals' [vector of class breaks]
%     'equalquantiles'   [number of classes]
%     'definedquantiles' [vector of quantiles]
%                        e.g. [0.1 0.9] results in three classes
%     'kmeans'           [number of classes]
%                        uses the k_means algorithm of Yi Cao
%                        FEX submission 19344 included in this function
%                        ! may take a while for large grids !
%     'std'              [s = scale of standard deviation, e.g. 2]
%                        equal intervals with a width of std/s. Bins are
%                        centered around sample mean
%     'otsu'             [number of classes, must be (2 4 8 16 etc)]
%                        recursive Otsu thresholding (see function
%                        graythresh)
%
% Output argument
%
%     C        Classified grid (class = GRIDobj)
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     C   = reclassify(DEM,'equalquantiles',10);
%     imageschs(DEM,C)
%
% See also: k_means, graythresh
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 18. January, 2013

narginchk(1,3)

allowedmethods = {'equalintervals',...
                  'definedintervals',...
                  'definedquantiles',...
                  'equalquantiles',...
                  'kmeans',...
                  'standarddeviation',...
                  'std',...
                  'otsu'...
                 };
% nr of classes
if nargin == 1;
    method = 'equalinterval';
    num = 10;
else
    method = validatestring(varargin{1},allowedmethods);

    if nargin==3
        num = varargin{2};
    end
end

INAN = isnan(DEM.Z) | isinf(DEM.Z);
DEM.Z(INAN) = nan;
inan = any(INAN(:));

switch method
    case 'equalintervals'
        validateattributes(num,{'numeric'},{'scalar'});

        if inan
            DEM.Z(~INAN) = mat2gray(DEM.Z(~INAN));
            DEM.Z = grayslice(DEM.Z,num);
            if num == 256;
                DEM.Z = double(DEM.Z) + 1;
                DEM.Z(INAN) = nan;
            elseif num < 256;
                DEM.Z = DEM.Z + 1;
                DEM.Z(INAN) = 0;
            else
                DEM.Z(INAN) = nan;
            end

        else
            DEM.Z = mat2gray(DEM.Z);
            DEM.Z = grayslice(DEM.Z,num);

            if num <= 255
                DEM.Z = DEM.Z + 1;
            elseif num == 256
                DEM.Z = double(DEM.Z) + 1;
            end
        end
    case 'definedintervals'
        validateattributes(num,{'numeric'},{'vector'});
        num(end+1) = inf;
        num = [-inf; num(:)];
        siz = size(DEM.Z);
        [~,DEM.Z] = histc(DEM.Z(:),num);
        DEM.Z = reshape(DEM.Z,siz);
    case 'equalquantiles'
        validateattributes(num,{'numeric'},{'scalar','integer','>',1});
        z = DEM.Z(~INAN);
        z = z(:);
        z = sort(z,'ascend');
        q = (1:num)/num;
        q = ceil(q*numel(z));
        edges = z(q);
        edges(end) = inf;
        [~,DEM.Z] = histc(DEM.Z(:),[-inf; edges]);
        DEM.Z = reshape(DEM.Z,DEM.size);
        DEM.Z(INAN) = nan;
    case 'definedquantiles'
        validateattributes(num,{'numeric'},{'vector','>',0,'<=',1});
        p = num(:);
        if p(end) < 1;
            p(end+1) = 1;
        end

        z = DEM.Z(~INAN);
        z = z(:);
        n = numel(z);
        edges = interp1((0:n-1).'/(n-1), sort(z), p);

        edges(end) = inf;
        [~,DEM.Z] = histc(DEM.Z(:),[-inf; edges]);
        DEM.Z = reshape(DEM.Z,DEM.size);
        DEM.Z(INAN) = nan;

    case {'standarddeviation','std'}
        validateattributes(num,{'numeric'},{'vector','>',0,'<=',1});
        z = DEM.Z(~INAN);
        z = z(:);
        s = std(z(:));
        s = s*num;
        mz = mean(z);
        minz = min(z);
        maxz = max(z);

        edges1 = [mz-s/2 :-s: minz];
        if edges1(end) > minz
            edges1(end+1) = -inf;
        else
            edges1(end) = -inf;
        end
        edges2 = mz+s/2 :s: maxz;

        edges2(end) = inf;

        [~,DEM.Z] = histc(DEM.Z(:),[edges1(end:-1:1) edges2]);
        DEM.Z = reshape(DEM.Z,DEM.size);
        DEM.Z(INAN) = nan;

    case 'kmeans';

        validateattributes(num,{'numeric'},{'scalar','integer','>',1});
        z = DEM.Z(~INAN);
        z = z(:);

        IX = k_means(z,num);

        DEM.Z(:,:) = nan;
        DEM.Z(~INAN) = IX;
    case 'otsu';
        validateattributes(num,{'numeric'},{'scalar','integer','>',1});
        if ceil(log2(num)) ~= log2(num);
            error('TopToolbox:GRIDobj','log2(value) must be an integer')
        end

        INAN = ~INAN;
        z  = mat2gray(DEM.Z(INAN));
        DEM.Z(INAN) = cast(otsu(z,num),class(DEM.Z));
end

end

% subfunctions

function [gIdx,c]=k_means(X,k)
% K_MEANS    k-means clustring
%   IDX = k_means(X, K) partititions the N x P data matrix X into K
%   clusters through a fully vectorized algorithm, where N is the number of
%   data points and P is the number of dimensions (variables). The
%   partition minimizes the sum of point-to-cluster-centroid Euclidean
%   distances of all clusters. The returned N x 1 vector IDX contains the
%   cluster indices of each point.
%
%   IDX = k_means(X, C) works with the initial centroids, C, (K x P).
%
%   [IDX, C] = k_means(X, K) also returns the K cluster centroid locations
%   in the K x P matrix, C.
%
% See also kmeans

% Version 2.0, by Yi Cao at Cranfield University on 27 March 2008.

% Example 1: small data set
%{
N=200;
X = [randn(N,2)+ones(N,2); randn(N,2)-ones(N,2)];
[cidx, ctrs] = k_means(X, 2);
plot(X(cidx==1,1),X(cidx==1,2),'r.',X(cidx==2,1),X(cidx==2,2),'b.', ctrs(:,1),ctrs(:,2),'kx');
%}

% Example 2: large data set
%{
N=20000;
X = [randn(N,2)+ones(N,2); randn(N,2)-ones(N,2)];
tic
[cidx, ctrs] = k_means(X, 2);
toc
plot(X(cidx==1,1),X(cidx==1,2),'r.',X(cidx==2,1),X(cidx==2,2),'b.', ctrs(:,1),ctrs(:,2),'kx');
%}

% Example 3: large data set with 5 centroids
%{
N=20000;
X = [randn(N,2)+ones(N,2); randn(N,2)-ones(N,2)];
tic
[cidx, ctrs] = k_means(X, 5);
toc
plot(X(cidx==1,1),X(cidx==1,2),'.',...
X(cidx==2,1),X(cidx==2,2),'.',...
X(cidx==3,1),X(cidx==3,2),'.',...
X(cidx==4,1),X(cidx==4,2),'.',...
X(cidx==5,1),X(cidx==5,2),'.',...
ctrs(:,1),ctrs(:,2),'+','linewidth',2)
%}

% Example 4: Comparison with kmeans in Statistics Toolbox
%{
N=20000;
X = [randn(N,2)+ones(N,2); randn(N,2)-ones(N,2)];
rand('state',0);
tic
cidx = k_means(X, 20);
toc
% Compare with kmeans in Statistis Toolbox
rand('state',0);
tic,
cidx1 = kmeans(X, 20, 'Option', statset('MaxIter',200));
toc
%}

% Check input and output
narginchk(2,2);

[n,m]=size(X);

% Check if second input is centroids
if ~isscalar(k)
    c=k;
    k=size(c,1);
else
    c=X(ceil(rand(k,1)*n),:);
end

% allocating variables
g0=ones(n,1);
gIdx=zeros(n,1);
D=zeros(n,k);

% Main loop converge if previous partition is the same as current
while any(g0~=gIdx)
%     disp(sum(g0~=gIdx))
    g0=gIdx;
    % Loop for each centroid
    for t=1:k
        d=zeros(n,1);
        % Loop for each dimension
        for s=1:m
            d=d+(X(:,s)-c(t,s)).^2;
        end
        D(:,t)=d;
    end
    % Partition data to closest centroids
    [~,gIdx]=min(D,[],2);
    % Update centroids using means of partitions
    for t=1:k
        c(t,:)=mean(X(gIdx==t,:));
    end
%     for t=1:m
%         c(:,t)=accumarray(gIdx,X(:,t),[],@mean);
%     end
end

end

function ix = otsu(x,n)

if mod(n,2)~=0
    error('number of classes is odd');
end

I = x >= graythresh(x);
ix = I*n/2;% + 1;
% do until all classes are
if n > 2
ix(I) = otsu(x(I),n/2) + ix(I);
ix(~I) = otsu(x(~I),n/2) + ix(~I);
else
    ix = ix+1;

end

end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function demoMultipleModelFiltering()
%%DEMOMULTIPLEMODELFILTERING This demonstrates how multiple model filters
%           can be used. The use of multiple dynamic models can improve the
%           accuracy of tracking targets than can maneuver. The scenario
%           considered is the simple 2D tracking of an aircraft. The
%           separate problems of data association, scheduling, track
%           initiation and termination are not addressed here.
%
%The example scenario is the air traffic control (ATC) scenario discussed
%in Chapter 11.7.4 of [1]. However, rather than using a fixed mode
%transition probability matrix, the use of the function
%getMarkovPTransProbMat to get the mode transition probability matrix given
%mean sojorn times and the time between samples is demonstrated. This
%allows the algorithms to be used with variable sampling rates, though in
%this instance, the sampling rate is held constant.
%
%Six scenarios are considered. One is a baseline Kalman filter with a
%polynomial dynamic model and a large process noise. The second is a Kalman
%filter with a first-order Gauss-Markov model (the integrated Ornstein-
%Uhlenbeck model), which, unlike the standard polynomial model, includes a
%correlation constant to try to better cover turns. The third is the
%reduced state filter. The fourth is the separated covariance filter. The
%other two scenarios utilize the interacting multiple model (IMM) filter.
%The first IMM scenario is as in [1], where two nearly constant velocity
%dynamic models with differing process noises are used in a standard Kalman
%filter (with the IMM). The second IMM scenario is where a nearly constant
%velocity model in a Kalman filter is used as well as a coordinated turn
%model in an extended Kalman filter (EKF). The functions multipleModelPred
%and multipleModelUpdate are used to handle the prediction and mixing in
%the IMM, though a change of the algorithm selection would allow one to
%just as easily use a generalized pseudo-Bayesian 2 filter. The IMM
%estimators can be more difficult to use as there are more parameters to
%tune.
%
%The measurements are generated in Cartesian space as measurement
%conversion and tracking using non-Cartesian measurements is not the focus
%of this example.
%
%REFERENCES:
%[1] Y. Bar-Shalom, X. R. Li, and T. Kiruabarajan, Estimation with
%    Applications to Tracking and Navigation. New York: Wiley Interscience,
%    2001.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%This can be changed to try other multiple model algorithms. If GPB1 is
%selected, the multiple model filter using the coordinate dturn model is
%not used, because the GPB1 requires that all models have the same
%dimensionality.
AlgSel='GPB2';

%Interval between samples in seconds.
T=5;

%%STEP 1: Generate the simulation scenario.
disp('Generating the "true" track for the simulation scenario.')
TTotal=125+90+125+30+125;
%The extra sample is the very first step in the trajectory.
numSamples=TTotal/T+1;
xTrue=zeros(4,numSamples);

%An indicator of the maneuvering mode 1 for not maneuvering, and 2 for
%maneuvering.
mode=zeros(numSamples,1);

%The plane first flies westward at 120 m/s for 125 seconds.
xInit=[25e3;10e3; -120; 0];
xTrue(:,1)=xInit;
mode(1)=1;
numStepSeg=125/T;%The number of steps in this segment.

%The state transition matrix for constant velocity motion.
F=FPolyKal(T,4,1);

%Propagate the constant velocity model for 125 seconds.
baseStep=1;
for curStep=(baseStep+1):(baseStep+numStepSeg)
    xTrue(:,curStep)=F*xTrue(:,curStep-1);
    mode(curStep)=1;
end
baseStep=baseStep+numStepSeg;

%Next, the plane performs a 1 degree per second coordinated turn for 90
%seconds.
turnRate=1*(pi/180);%Convert from degrees per second to radians per second.

numStepSeg=90/T;
for curStep=(baseStep+1):(baseStep+numStepSeg)
    %Get the discrete-time propagation matrix for the state. The
    %FCoordTurn2D function assumes that the turn rate is part of the state.
    %In this instance, we just add it to the state, and then remove the
    %extra elements from the returned matrix.
    F=FCoordTurn2D(T,[xTrue(:,curStep-1);turnRate]);
    F=F(1:4,1:4);

    xTrue(:,curStep)=F*xTrue(:,curStep-1);
    mode(curStep)=2;
end
baseStep=baseStep+numStepSeg;

%Next, go straight (South after the previous turn) for another 125 seconds.
numStepSeg=125/T;%The number of steps in this segment.

%The state transition matrix for constant velocity motion.
F=FPolyKal(T,4,1);

for curStep=(baseStep+1):(baseStep+numStepSeg)
    xTrue(:,curStep)=F*xTrue(:,curStep-1);
    mode(curStep)=1;
end
baseStep=baseStep+numStepSeg;

%Make another turn. This time, it is at a rate of -three degrees per second
%for 30 seconds.
turnRate=-3*(pi/180);%Degrees per second to radians per second.

numStepSeg=30/T;%The number of steps in this segment.
for curStep=(baseStep+1):(baseStep+numStepSeg)
    F=FCoordTurn2D(T,[xTrue(:,curStep-1);turnRate]);
    F=F(1:4,1:4);

    xTrue(:,curStep)=F*xTrue(:,curStep-1);
    mode(curStep)=2;
end
baseStep=baseStep+numStepSeg;

%Finally, go straight again for another 125 seconds.
numStepSeg=125/T;%The number of steps in this segment.

%The state transition matrix for constant velocity motion.
F=FPolyKal(T,4,1);

for curStep=(baseStep+1):(baseStep+numStepSeg)
    xTrue(:,curStep)=F*xTrue(:,curStep-1);
    mode(curStep)=1;
end

%%STEP 2: Generate measurements.
disp('Generating the measurements for the simulation; A simple Cartesian model is used.')
%The measurements are just Cartesian with 100 meters of noise added per
%dimension, with no correlation between them.
H=[1,0,0,0;
   0,1,0,0];%Measurement matrix

SR=[100, 0;
    0,  100];%Square root covariance matrix.
R=SR*SR';%The covariance matrix.

z=zeros(2,numSamples);
for curSamp=1:numSamples
    z(:,curSamp)=H*xTrue(:,curSamp)+SR*randn(2,1);
end

%Plot the trajectory and the measurements
disp('Plotting the trajectory the measurements, and the true maneuver mode.')
figure(1)
clf
hold on
plot(xTrue(1,:),xTrue(2,:),'-b','linewidth',2)
scatter(z(1,:),z(2,:),'.k')

h1=xlabel('Meters West->East');
h2=ylabel('Meters South->North');
h3=title('The True Trajectory and the Detections');
set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')

%Plot when it is in Mode 2, 1 for being in it and 2 for not.
figure(2)
clf
hold on
plot(mode-1,'-k')
h1=xlabel('Discrete Step');
h2=ylabel('Maneuver mode');
h3=title('Indicating when the target is performing maneuvers.');
set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')

%STEP 3: Get the baseline scenario: Run a Kalman filter with a standard
%first-order white-noise process model, a Kalman filter with a first-order
%Gauss-Markov model, a reduced state estimator, and the
%separated covariance filter.
disp('Computing the baseline scenarios with just a Kalman filter and with the')
disp('reduced state estimator.')

%A direct-discrete model with a 1m/2s^2 process noise standard deviation is
%used (normally, a discretized model is preferred, because it is
%more consistent when using a variable sampling rate, and the covariance
%matrix is not singular).
Q=QPolyKalDirectDisc(T,4,1,1^2);%1m/s^2 process noise

%The state transition matrix for constant velocity motion.
F=FPolyKal(T,4,1);

%Parameters for the Gauss-Markov dynamic model. We are using a first-order
%Gauss-Markov model.
tau=20;%20 seconds; the assumed maneuver decorrelation time.
maxAccel=9.8*3;%Assume max 3G turn
%Rule of thumb- process noise suggestion.
q=processNoiseSuggest('PolyKal-ROT',maxAccel,T);
QGM=QGaussMarkov(T,4,q,tau,1);%Process noise matrix.
FGM=FGaussMarkov(T,4,tau,1);%State Transition matrix

%Allocate space for the polynomial Kalman filter estimates
xKalman=zeros(4,numSamples);%The states
PKalman=zeros(4,4,numSamples);%The covariance matrices

%Allocate space for the Gauss-Markov Kalman filter estimates
xKalmanGM=zeros(4,numSamples);%The states
PKalmanGM=zeros(4,4,numSamples);%The covariance matrices

%Allocate space for the reduced state estimator.
xRedState=zeros(4,numSamples);%The states

%Allocate space for the separated covariance filter.
xSCFState=zeros(4,numSamples);%The states

%We are using the version of the reduced state estimator that needs a
%maximum linear acceleration (directed along the direction of motion of the
%target) for the target along with a maximum turn rate acceleration.
%A 0.5m/s^2 maximum linear acceleration with 1 deg/s^2 maximum turn rate
%acceleration is used.
ARedState=0.5;
OmegaRedState=1*(pi/180);

%The separated covariance filter needs a matrix introducing the effects of
%the maximum acceleration onto the state (for covariance computation
%results. For a state consisting of position and velocity in 2D, this can
%be
aMax=9.8/2;
Ba=[T^2/2,   0;
    0,       T^2/2;
	T,       0;
	0,       T]*aMax;
%where we chose the maximum acceleration in each dimension to be 0.5G
%(1G=9.8m/s^2).

%Track initiation is by one-point differencing using just assumed maximum
%values as the standard deviations for the velocity covariance matrix. The
%assumed maximum value is 300m/s.
[xKalman(:,1),PKalman(:,:,1)]=onePointCartInit(z(:,1),SR,300);
xKalmanGM(:,1)=xKalman(:,1);
PKalmanGM(:,1)=PKalman(:,1);

%The reduced state estimator keeps the covariance matrix broken into parts
%M and D, which will be propagated and updated.
xRedState(1:2)=z(:,1);
M=PKalman(:,:,1);
D=zeros(4,2);

%The separated covariance filter. It has been chosen to start the filter
%with zero lag and all uncertainty in the covariance matrix.
xSCFState(1:2)=z(:,1);
PSCF=PKalman(:,:,1);
LSCF=zeros(4,2);

absErrKalman=zeros(numSamples,1);%Allocate space
absErrKalman(1)=norm(xKalman(1:2,1)-xTrue(1:2,1));

absErrKalmanGM=zeros(numSamples,1);%Allocate space
absErrKalmanGM(1)=norm(xKalman(1:2,1)-xTrue(1:2,1));

absErrRedState=zeros(numSamples,1);%Allocate space
absErrRedState(1)=norm(xRedState(1:2,1)-xTrue(1:2,1));

absErrSCF=zeros(numSamples,1);%Allocate space
absErrSCF(1)=norm(xSCFState(1:2,1)-xTrue(1:2,1));

%Run the filter
for curSamp=2:numSamples
    %The Kalman filter with the polynomial model
    %Predict the state forward
    [xPred, PPred]=discKalPred(xKalman(:,curSamp-1),PKalman(:,:,curSamp-1),F,Q);
    %Update the state with a measurement
    [xKalman(:,curSamp),PKalman(:,:,curSamp)]=KalmanUpdate(xPred,PPred,z(:,curSamp),R,H);

    absErrKalman(curSamp)=norm(xKalman(1:2,curSamp)-xTrue(1:2,curSamp));

    %The Kalman filter with the Gauss-Markov model
    [xPred, PPred]=discKalPred(xKalmanGM(:,curSamp-1),PKalmanGM(:,:,curSamp-1),FGM,QGM);
    %Update the state with a measurement
    [xKalmanGM(:,curSamp),PKalmanGM(:,:,curSamp)]=KalmanUpdate(xPred,PPred,z(:,curSamp),R,H);

    absErrKalmanGM(curSamp)=norm(xKalmanGM(1:2,curSamp)-xTrue(1:2,curSamp));

    %The reduced state estimator.
    %Predict the state forward
    modParams=[];
    modParams.T=T;
    modParams.A=ARedState;
    modParams.Omega=OmegaRedState;
    [xPred,M,D,PPred]=reducedStateDiscPred(xRedState(:,curSamp-1),M,D,[],[],'GenTurn',modParams);
    %Update the state with a measurement
    [xRedState(:,curSamp),M,D]=reducedStateUpdate(xPred,PPred,M,D,z(:,curSamp),R);

    absErrRedState(curSamp)=norm(xRedState(1:2,curSamp)-xTrue(1:2,curSamp));

    %The separated covariance filter.
    %Predict the state forward
    [xPred,LSCF,TSCF]=separatedCovDiscPred(xSCFState(:,curSamp-1),PSCF,LSCF,F,Ba);
    %Update the state with a measurement
    [xSCFState(:,curSamp),LSCF,PSCF]=separatedCovUpdate(xPred,LSCF,TSCF,z(:,curSamp),R,H);
    absErrSCF(curSamp)=norm(xSCFState(1:2,curSamp)-xTrue(1:2,curSamp));
end

figure(1)
hold on
plot(xKalman(1,:),xKalman(2,:),'-g','linewidth',2)
plot(xKalmanGM(1,:),xKalmanGM(2,:),'-c','linewidth',2)
plot(xRedState(1,:),xRedState(2,:),'--b','linewidth',2)
plot(xSCFState(1,:),xSCFState(2,:),'-.m','linewidth',2)

figure(3)
clf
hold on
plot(absErrKalman,'-g','linewidth',2)
plot(absErrKalmanGM,'-c','linewidth',2)
plot(absErrRedState,'--b','linewidth',2)
plot(absErrSCF,'-.m','linewidth',2)

h1=xlabel('Discrete Step');
h2=ylabel('Absolute distance error');
set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')

%%STEP 4: Run the algorithms with the two linear models.
%The linear model pair consists of trackers using direct discrete linear
%dynamic models.
disp('Computing the IMM scenario with two linear Kalman filters with different process noises.')

%The covariance matrix for the low-noise model.
QLow=QPolyKalDirectDisc(T,4,1,0.1^2);%0.1m/s^2 process noise
%The covariance matrix for the high-noice model.
QHigh=QPolyKalDirectDisc(T,4,1,2^2);%2m/s^2 process noise.

%The propagation routines for each model
transFuns{1}=@(xPrev,PPrev)discKalPred(xPrev,PPrev,F,QLow);
transFuns{2}=@(xPrev,PPrev)discKalPred(xPrev,PPrev,F,QHigh);

%The measurement update function is the same for both models.
measUpdateFuns=@(x,P,z,R)KalmanUpdate(x,P,z,R,H);

%Allocate space for the IMM model states. This is the size of the
%largest state by the number of models, except when using the GPB1
%estimator, in which case it is just the size of the state (all states have
%to be the same size).
if(~strcmp(AlgSel,'GPB1'))
    xIMML=zeros(4,2);
    PIMML=zeros(4,4,2);

    %Initialize the states using one-point differencing. The covariance matrix
    %is set using the assumed velocity for each of the models, which is
    %300m/s^2.
    %The first model.
    [xIMML(:,1),PIMML(:,:,1)]=onePointCartInit(z(:,1),SR,300);
    %The second model has the same initialization.
    xIMML(:,2)=xIMML(:,1);
    PIMML(:,:,2)=PIMML(:,:,1);
else
    [xIMML,PIMML]=onePointCartInit(z(:,1),SR,300);
end

%Allocate space for the states to display. These will be the merged values
%from the IMM model states.
xIMMLDisp=zeros(4,numSamples);
PIMMLDisp=zeros(4,4,numSamples);

%Initially, both models have the same initialization, which will thus be
%the one to display at the first step.
xIMMLDisp(:,1)=xIMML(:,1);
PIMMLDisp(:,:,1)=PIMML(:,:,1);

%Allocate space for the mode probabilities. These are saved over all scans
%so that they can be plotted.
muMode=zeros(2,numSamples);%There is one for each model at each step.
%Initially, the mode probabilities are uniform: we do not know which mode
%it is in.
muMode(:,1)=1/2;

%The mode transition probability matrix for the two linear models must be
%determined. This is based on a matrix relating the mean sojourn times of
%the models and the time between samples.
%The mean sojourn time matrix is a design parameter. The negative inverse
%of the diagonal terms is the mean time spent in that state. Each row must
%sum to one. If the matrix were larger, the relation between the
%non-diagonal values in the rows would affect to which state it
%transitions. This is a design parameter.
A=[-0.01/2, 0.01/2;
    0.02*2, -0.02*2];

%Compute the mode transition probability matrix for the given time between
%samples. Because the sample interval is constant, this does not change
%over time.
LambdaL=getMarkovPTransProbMat(A,T);

absErrIMML=zeros(numSamples,1);%Allocate space
absErrIMML(1)=norm(xIMMLDisp(1:2,1)-xTrue(1:2,1));
%Run the filter
for curSamp=2:numSamples
   %Predict the states forward
   [xIMML,PIMML]=multipleModelPred(AlgSel,xIMML,PIMML,transFuns);

   %Update the state with a measurement
   [xIMML,PIMML,muMode(:,curSamp),xIMMLDisp(:,curSamp),PIMMLDisp(:,:,curSamp)]=multipleModelUpdate(AlgSel,xIMML,PIMML,z(:,curSamp),R,measUpdateFuns,muMode(:,curSamp-1),LambdaL);
   absErrIMML(curSamp)=norm(xIMMLDisp(1:2,curSamp)-xTrue(1:2,curSamp));
end

figure(1)
plot(xIMMLDisp(1,:),xIMMLDisp(2,:),'-r','linewidth',2)

figure(2)
plot(muMode(2,:),'-r','linewidth',2)

figure(3)
plot(absErrIMML,'-r','linewidth',2)

if(strcmp(AlgSel,'GPB1'))
    return
end

%STEP 5: Run the algorithms with the linear and the coordinated turn model.
disp('Computing the IMM scenario with a linear Kalman filter and an EKF.')

%The covariance matrix for the linear model.
QLinear=QPolyKalDirectDisc(T,4,1,0.01^2);%0.01m/s^2 process noise.

%0.025m/s^2 linear acceleration (in 3D, not just directed along the
%direction of motion of the target) with 0.5 deg/s^2 turning acceleration
%for the coordinated turn model as the standard deviations going into the
%covariance matrix. Note that these assumptions differ from those used in
%the reduced state filter. The filters should be tuned separately for a
%fair comparison. Often, improving the IMM's ability to recognize a turn
%worsens its overall track performance. Thus, a tradeoff must be present.
QCT=QCoordTurn(T,zeros(5,1),0.025^2,(0.5*(pi/180))^2);

%The propagation routines for each model
transFuns=[];
transFuns{1}=@(xPrev,PPrev)discKalPred(xPrev,PPrev,F,QLinear);
f=@(x)(FCoordTurn2D(T,x)*x);
transFuns{2}=@(xPrev,PPrev)discEKFPred(xPrev,PPrev,f,@(x)JacobCoordTurn2D(T,x),QCT);

%The measurement update function differs for the different models, because
%H is different for the coordinated turn model, due to the extra element
%in the state.
HCT=[1,0,0,0,0;
     0,1,0,0,0];
measUpdateFuns=[];
measUpdateFuns{1}=@(x,P,z,R)KalmanUpdate(x,P,z,R,H);
measUpdateFuns{2}=@(x,P,z,R)KalmanUpdate(x,P,z,R,HCT);

%The number of elements in each state and the number of dimensions of each
%state that play a role in mixing the states.
numStateDims=[4;5];
numMixDims=[4;4];

%Allocate space for the IMM model states. The number of rows is the size of
%the longest state.
xLCT=zeros(5,2);
PLCT=zeros(5,2);

%Initialize the states using one-point differencing. The covariance matrix
%is set using the assumed velocity for each of the models, which is
%300m/s^2 and the maximum turn rate, which is assumed to be 20 degrees per
%second.
%The first model.
xLCT(1:2,1)=z(:,1);
PLCT(1:2,1:2,1)=R;
PLCT(3,3,1)=300^2;
PLCT(4,4,1)=300^2;
%The second model has the same initialization.
xLCT(:,2)=xLCT(:,1);
PLCT(:,:,2)=PLCT(:,:,1);
PLCT(5,5,2)=(20*(pi/180))^2;

%Allocate space for the states to display. These will be the merged values
%from the IMM model states. The size is the size of the state with the
%least number of mixing components.
xLCTDisp=zeros(4,numSamples);
PLCTDisp=zeros(4,4,numSamples);

%Initially, both models have the same initialization, which will thus be
%the one to display at the first step.
xLCTDisp(:,1)=xLCT(1:4,1);
PLCTDisp(:,:,1)=PLCT(1:4,1:4,1);

%Allocate space for the mode probabilities. These are saved over all scans
%so that they can be plotted.
muModeLCT=zeros(2,numSamples);%There is one for each model at each step.
%Initially, the mode probabilities are uniform: we do not know which mode
%it is in.
muModeLCT(:,1)=1/2;

%Compute the mode transition probability matrix.
A=[-0.01/2, 0.01/2;
    0.02*2, -0.02*2];
LambdaCT=getMarkovPTransProbMat(A,T);

absErrLCT=zeros(numSamples,1);%Allocate space
absErrLCT(1)=norm(xLCTDisp(1:2,1)-xTrue(1:2,1));
%Run the filter
for curSamp=2:numSamples
   %Predict the states forward.
   [xLCT,PLCT]=multipleModelPred(AlgSel,xLCT,PLCT,transFuns,numStateDims,numMixDims);

   %Update the state with a measurement.
   [xLCT,PLCT,muModeLCT(:,curSamp),xLCTDisp(:,curSamp),PLCTDisp(:,:,curSamp)]=multipleModelUpdate(AlgSel,xLCT,PLCT,z(:,curSamp),R,measUpdateFuns,muModeLCT(:,curSamp-1),LambdaCT,numStateDims,numMixDims);
   absErrLCT(curSamp)=norm(xLCTDisp(1:2,curSamp)-xTrue(1:2,curSamp));
end

figure(1)
plot(xLCTDisp(1,:),xLCTDisp(2,:),'--c','linewidth',2)
legend('True Trajectory','Raw Detections','Basic Kalman Filter', 'Gauss-Markov Kalman Filter','Reduced State Filter', 'Separated Covariance Filter', 'Two Linear Models','Linear and Maneuvering Models','Location','SouthEast')

figure(2)
plot(muModeLCT(2,:),'--c','linewidth',2)
legend('True Mode', 'Two Linear Models','Linear and Maneuvering Models','Location','NorthWest')

figure(3)
plot(absErrLCT,'--c','linewidth',2)
legend('Polynomial Kalman Filter','Gauss-Markov Kalman Filter','Reduced State Filter', 'Separated Covariance Filter', 'Two Linear Models','Linear and Maneuvering Models','Location','NorthWest')
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% File conversion
% Reading images from CASIA database and converting to 'jpg' format
% Author: Qingbao Guo
% Gj?vik University College, GUC, Norway

clear all;

path='CASIA Iris Image Database (version 1.0)\00'; %please update to your directory
subpath='00'; %please update to your directory
for i=1:108

            if i>=10
                path='CASIA Iris Image Database (version 1.0)\0';
                subpath='0';
            end

            if i>=100
                 path='CASIA Iris Image Database (version 1.0)\';
                 subpath='';
            end

    for j=1:2

        if j==1
        for k=1:3

        filesrcpath1 = strcat(path,num2str(i),'\',num2str(j),'\',subpath,num2str(i),'_',num2str(j),'_',num2str(k),'.bmp');
        filedespath1 = strcat(path,num2str(i),'_',num2str(j),'_',num2str(k),'.jpg');
        im = imread(filesrcpath1,'bmp');
        data(:,7*(i-1)+k) = reshape(im,1,280*320);
        imwrite(im,filedespath1,'jpg');
        end
        end

        if j==2
        for h=1:4
        filesrcpath2 = strcat(path,num2str(i),'\',num2str(j),'\',subpath,num2str(i),'_',num2str(j),'_',num2str(h),'.bmp');
        filedespath2 = strcat(path,num2str(i),'_',num2str(j),'_',num2str(h),'.jpg');
        im = imread(filesrcpath2,'bmp');
        data(:,7*i-4+h) = reshape(im,1,280*320);
        imwrite(im,filedespath2,'jpg');
        end
        end

    end
end
save data.mat data;

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris


text:
classdef TestLineSegmentDetector
    %TestLineSegmentDetector

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','blox.jpg');
    end

    methods (Static)
        function test_1
            img = cv.imread(TestLineSegmentDetector.im, ...
                'Grayscale',true, 'ReduceScale',2);
            lsd = cv.LineSegmentDetector('Refine','Advanced');

            [lines, width, prec, nfa] = lsd.detect(img);
            validateattributes(lines, {'cell'}, {'vector'});
            cellfun(@(v) validateattributes(v, {'numeric'}, ...
                {'vector', 'numel',4, 'real'}), lines);
            validateattributes(width, {'numeric'}, {'vector', 'numel',numel(lines)});
            validateattributes(prec, {'numeric'}, {'vector', 'numel',numel(lines)});
            validateattributes(nfa, {'numeric'}, {'vector', 'numel',numel(lines)});

            out = lsd.drawSegments(img, lines);
            validateattributes(out, {class(img)}, {'3d', 'size',[size(img) 3]});
        end

        function test_2
            img = cv.imread(TestLineSegmentDetector.im, ...
                'Grayscale',true, 'ReduceScale',2);
            lsd1 = cv.LineSegmentDetector('Refine','Standard');
            lsd2 = cv.LineSegmentDetector('Refine','None');
            lines1 = lsd1.detect(img);
            lines2 = lsd2.detect(img);
            [comp,c] = lsd1.compareSegments([size(img,2) size(img,1)], lines1, lines2);
            validateattributes(comp, {class(img)}, {'3d', 'size',[size(img) 3]});
            validateattributes(c, {'numeric'}, {'scalar', 'integer', 'nonnegative'});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function p_sensor = directionalResponse(kgrid, sensor, sensor_mask_index, p_k)
%DIRECTIONALRESPONSE    Apply directivity to the sensor measurements.
%
% DESCRIPTION:
%       directionalResponse takes a pressure field in k-space, p_k,
%       multiplies by a directivity function, and outputs the grid points
%       indicated in sensor.mask.
%
% USAGE:
%       p_sensor = directionalResponse(kgrid, sensor, p_k)
%
% INPUTS:
%       kgrid       - k-Wave grid structure returned by makeGrid
%
%       sensor      - k-Wave sensor structure containing the following
%                     fields:
%
%           .directivity_angle
%                   - a matrix with the same structure as
%                     sensor.mask that allocates a directivity angle to
%                     each sensor element as defined in sensor.mask. The
%                     angles are in radians:
%                     0 = max sensitivity in y direction (up/down)
%                     pi/2 or -pi/2 = max sensitivity in x direction (left/right)
%
%       	.directivity_unique_angles
%                   - list of the unique directivity angles returned by
%                     sensor.directivity_unique_angles =
%                     unique(sensor.directivity_angle(sensor.mask == 1));
%
%       	.directivity_pattern
%                   - a text string with currently only one
%                     option. 'pressure' indicates that the directional
%                     response should be of the kind due to spatial
%                     averaging over a sensor surface, so a sinc function
%                     in 2D.
%
%       	.directivity_size
%                   - the directivity pattern used is what
%                     would be the directivity if the sensor were this
%                     length (width). The larger this is the more
%                     directional the response.
%
%       	.directivity_wavenumbers
%                   - this is set to [kgrid.ky(:)'; kgrid.kx(:)']. It is
%                     precomputed to allow data casting, as kgrid.kx (etc)
%                     are computed on the fly.
%
%       sensor_mask_index
%                   - indices of the active sensor elements in sensor.mask
%
%       p_k         - the acoustic pressure field in the k-space domain.
%
% Currently works for binary sensor_mask, but not when sensor_mask is given
% as Cartesian coordinates. Also, currently works only in 2D.
%
% ABOUT:
%       author: Ben Cox and Bradley Treeby
%       date: 21st January 2010
%       last update: 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% sensor mask indices
Ns = length(sensor_mask_index);

% pre-allocate p_sensor vector
p_sensor = zeros(Ns, 1, class(sensor.directivity_angle));

% loop over number of unique angles
for loop = 1:length(sensor.directivity_unique_angles);

   % get current angle
   theta = sensor.directivity_unique_angles(loop);

   % find which of the sensors have this directivity
   indices = find(sensor.directivity_angle(sensor_mask_index) == theta);

   switch sensor.directivity_pattern
      case 'pressure'
         % calculate magnitude of component of wavenumber along sensor face
         k_tangent = reshape([cos(theta) -sin(theta)]*sensor.directivity_wavenumbers, kgrid.Nx, kgrid.Ny);
         directionality = fftshift(sinc(k_tangent*sensor.directivity_size/2));
      case 'gradient'
         % calculate magnitude of component of wavenumber normal to the sensor face
         k_normal = reshape([sin(theta), cos(theta)]*[kgrid.ky(:)'; kgrid.kx(:)'], kgrid.Nx, kgrid.Ny);
         temp = k_normal./kgrid.k;
         temp(kgrid.k==0) = 0;
         directionality = fftshift(temp);
      otherwise
         error('Unsupported directivity pattern')
   end

   % apply the directivity response to the pressure field (in k-space)
   p_directivity = real(ifft2(p_k.*directionality));

   % pick out the response at the sensor points
   p_sensor(indices) = p_directivity(sensor_mask_index(indices));

end

repo: Aero-matlab


text:
function [Rob,Sen,Lmk,Obs,Frm,Fac] = solveGraphQR(Rob,Sen,Lmk,Obs,Frm,Fac,options)

% SOLVEGRAPHQR Solves the SLAM graph using QR decomposition.
%   [Rob,Sen,Lmk,Obs,Frm,Fac] = SOLVEGRAPHQR(Rob,Sen,Lmk,Obs,Frm,Fac)
%   solves the graph-SLAM problem using QR decomposition of the
%   Hessian matrix.
%
%   IMPORTANT NOTE: This method is illustrative and constitutes the
%   motivation for this toolbox. One can achieve better performances, both
%   in computing time and possibly in robustness and accuracy, by using
%   Matlab's built-in nonlinear optimization tools, such as LSQNONLIN.
%
%   See courseSLAM.pdf in the documentation for details about the QR
%   decomposition for solving the graph-SLAM problem.
%
%   See also SOLVEGRAPHCHOLESKY, ERRORSTATEJACOBIANS, UPDATESTATES,
%   COMPUTEERROR, COMPUTERESIDUAL, COLAMD, '\', MLDIVIDE, LSQNONLIN.

% Copyright 2015-    Joan Sola @ IRI-UPC-CSIC.

global Map

% Control of iterations and exit conditions
n_iter      = options.niterations; % exit criterion of number of iterations
target_dres = options.target_dres; % exit criterion for error variation
target_res  = options.target_res;  % exit criterion for current residual
res_old     = 1e10;                   % last iteration's error

% Map states range
Map.mr = find(Map.used);
% Map factors range
errs = [Fac.err];
Map.fr = 1:sum([errs.size]);

for it = 1:n_iter

    % Compute Jacobians for projection onto the manifold
    [Frm,Lmk] = errorStateJacobians(Frm,Lmk);

    % Build Hessian A and rhs vector b, in global Map
    Fac       = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac);

    if it == 1 % do this only once:

        % Column permutation
        p  = colamd(Map.A(Map.fr,Map.mr))';

        % Permutated map range
        pr = Map.mr(p);
        Map.pr = pr;

    end

    % Decomposition
    [Map.d, Map.R] = qr(Map.A(Map.fr,pr), Map.b(Map.fr), 0);

    % Solve for dx and reorder:
    %   - dx is Map.x(mr)
    %   - reordered dx is Map.x(pr)
    Map.x(pr) = -Map.R\Map.d; % solve for dx;

    % NOTE: Matlab is able to do all the reordering and QR factorization
    % for you. If you just use the operator '\', as in 'dx = -A\b', Matlab
    % will reorder A, then factor it to get R and d, then solve the
    % factored problem, then reorder back the result into dx. Use the
    % following line to accomplish this, and comment out the code from line
    % 'if it == 1' until here:
    %
    %     Map.x(Map.mr) = -Map.A(Map.fr,Map.mr)\Map.b(Map.fr);

    % Update nominal states
    [Rob,Lmk,Frm] = updateStates(Rob,Lmk,Frm);

    % Check resulting errors
    [res, err_max] = computeResidual(Rob,Sen,Lmk,Obs,Frm,Fac);
    dres           = res - res_old;
    res_old        = res;

    % Test and exit
    if ( ( -dres <= target_dres ) || (err_max <= target_res) ) %&& ( abs(derr) < target_derr) )
        break;
    end

end

% Compute full covariance matrix.
Map.P(pr,pr) = inv(full(Map.R))*inv(full(Map.R))';

end

function Fac = buildProblem(Rob,Sen,Lmk,Obs,Frm,Fac)

% BUILDPROBLEM Build least squares problem's matrix A and vector b
%   Fac = BUILDPROBLEM(Rob,Sen,Lmk,Obs,Frm,Fac) Builds the least squares
%   problem's matrix A and vector b for a solution using sparse QR
%   factorization of A.

global Map

% Reset Hessian and rhs vector
Map.A(Map.fr,Map.mr) = 0;
Map.b(Map.fr)    = 0;

% Iterate all factors
facCount = 1;
for fac = find([Fac.used])

    % Extract some pointers
    rob    = Fac(fac).rob;
    sen    = Fac(fac).sen;
    lmk    = Fac(fac).lmk;
    frames = Fac(fac).frames;

    % Compute factor error, info mat, and Jacobians
    [Fac(fac), e, ~, Wsqrt, J1, J2, r1, r2] = computeError(...
        Rob(rob),       ...
        Sen(sen),       ...
        Lmk(lmk),       ...
        Obs(sen,lmk),   ...
        Frm(frames),    ...
        Fac(fac));

    % row band matrix size
    m  = numel(e);
    mr = (facCount : facCount + m - 1);

    % Update A and b
    Map.A(mr,r1) = Wsqrt * J1;
    Map.A(mr,r2) = Wsqrt * J2;

    Map.b(mr,1)  = Wsqrt * e;

    % Advance to next row band
    facCount = facCount + m;

end

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function momftfrt
%MOMFTFRT Unit test for the function MOMFTFR.

%	O. Lemoine - March 1996.

N=128;

% For a perfect line
[sig,ifl]=fmlin(N);
tfr=tfrideal(ifl);
[tm,D2]=momftfr(tfr);
tmth=(1:N);
if any(abs(tm-tmth')>sqrt(eps)),
  error('momftfr test 1 failed');
elseif any(abs(D2)>sqrt(eps)),
  error('momftfr test 2 failed');
end

% For a sinusoid
[sig,ifl]=fmsin(N,0,.5,N*2,1,0);
tfr=tfrideal(ifl);
[tm,D2]=momftfr(tfr');
tmth=round(ifl*2*(N-1))+1;
if any(abs(tm-tmth)>sqrt(1/N)),
  error('momftfr test 3 failed');
elseif any(abs(D2)>sqrt(eps)),
  error('momftfr test 4 failed');
end

N=117;

% For a perfect line
[sig,ifl]=fmlin(N);
tfr=tfrideal(ifl);
[tm,D2]=momftfr(tfr);
tmth=(1:N);
if any(abs(tm-tmth')>sqrt(eps)),
  error('momftfr test 5 failed');
elseif any(abs(D2)>sqrt(eps)),
  error('momftfr test 6 failed');
end

% For a sinusoid
[sig,ifl]=fmsin(N,0,.5,N*2,1,0);
tfr=tfrideal(ifl);
[tm,D2]=momftfr(tfr');
tmth=round(ifl*2*(N-1))+1;
if any(abs(tm-tmth)>sqrt(1/N)),
  error('momftfr test 7 failed');
elseif any(abs(D2)>sqrt(eps)),
  error('momftfr test 8 failed');
end

author: HeLiangHIT


text:
function [ gx ] = ObsRecGen(x,P,u,in)
% observation function for k-ToM's bet about her opponent's next move
% [ gx ] = ObsRecGen(x,P,u,in)
% Marie Devaine wrote this function in November 2015 (comments: JD).
% A k-ToM learner bases her decision (a=1 or a=0) upon her prediction of
% her opponent's next move, given the game payoff table. When k>1, k-ToM
% maintains more than one such prediction, which depends upon the possible
% level of her opponent. Let P(o=1) be the probability that k-ToM's
% opponent will pick the first alternative option. Then:
% P(o=1) = sum_k P(o=1|k)*P(k)
% where P(o=1|k) is the probability that k-ToM's opponent will pick
% the first alternative option if he was a k-ToM, and P(k) is the
% probability that k-ToM's opponent is a k-ToM.
% IN:
%   - x: hidden states (see indexing in inG.indlev)
%   - P: observation param:
%       P(1) = (log-) temperature
%       P(2) = bias [optional]
%   - u: [useless]
%   - inG: input structure (see prepare_kToM.m)
% OUT:
%   - gx: proba that the agent will pick the first option, i.e. gx=P(y=1).

player = in.player; % 1 or 2: role of the player
ntotPar = in.npara; % only for k-ToM with k>0
level = in.lev; % depth of k-ToM's recursive beliefs
game = in.game; % payoff table
a = 0.36; % for E[s(x)] when x~n(mu,Sig)
indlev = in.indlev; % hidden-states indexing [see defIndlev.m]

% Get the agent's prediction about her opponent's next move, ie P(o=1).
if level==0 % 0-ToM

    mx = x(1); % E[log-odds of P(o=1)]
    Vx = exp(x(2)); % V[log-odds of P(o=1)]
    Po = VBA_sigmoid(mx/(sqrt(1+a*Vx))); % P(o=1)

else

    % Get P(k'). Note: if the agent is k-ToM, then, by definition, she
    % considers that her opponent's sophistication is k' < k. In addition,
    % there is a constraint of normalization, ie sum_k' P(k') = 1. Thus,
    % one only needs to keep track of k'-1 probabilities (the last one is,
    % by construction, 1-sum_k' P(k')).
    Pk = VBA_sigmoid(x(1:(level-1))); % P(k'), with k'=0,...,k-1
    Pk = [Pk;max(0,1-sum(Pk))]; % insert last P(k'=k-1)

    % Get P(o=1|k'). Note: the agent's prediction P(o=1|k') depends upon
    % her estimate of her opponent's parameters (learning rate, tmperature,
    % bias...). Uncertainty Re: these parameters eventually results in
    % blurring her prediction.
    % Note: hidden states encode x(theta), the log-odds of P(o=1|k',theta)
    % evaluated at the agent's estimate of theta (mu). In addition, they
    % encode the gradient of x wrt to theta (dx/dtheta), and V[theta]. This
    % then serves to derive P(o=1|k') as follows:
    % P(o=1|k') = E[sigm(x(theta))]
    %           = sigm(E[x(theta)]/sqrt(1+a*V[x(theta)])
    %           = sigm(E[x(theta)]/sqrt(1+a*V[theta]*(dx/dtheta)^2)
    f = zeros(level,1); % E[x(theta)]
    Vx = zeros(level,1); % V[x(theta)]
    for j=1:level % loop over possible opponent's levels  (k'=j-1)
        f(j) = x(indlev(j).f); % E[x(theta)|k'=j-1]
        df = x(indlev(j).df); % d[x(theta)]/dtheta for k'=j-1
        Sig = exp(x(indlev(j).Par(2:2:2*ntotPar))); % V[theta|k'=j-1]
        Vx(j) = sum(Sig.*df.^2); % V[x(theta)|k'=j-1]
    end
    Es = VBA_sigmoid(f./sqrt(1+a*Vx)); % E[sigm(x(theta))]

    % Get P(o=1) = sum_k P(o=1|k')*P(k')
    Po = Pk'*Es; % k-ToM's belief about her opponent's next move

end

% Make decision based upon P(o=1)
DV = fplayer(Po,exp(P(1)),player,game); % incentive for a=1
if length(P)==1
    gx = VBA_sigmoid(DV); % P(a=1)
else % P(2) = bias
    gx = VBA_sigmoid(DV+P(2)); % P(a=1) with bias
end

author: MBB-team

repo: VBA-toolbox


text:
function d = month_to_nones_roman ( m )

%*****************************************************************************80
%
%% MONTH_TO_NONES_ROMAN returns the day of the nones of a Roman month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the month index.
%
%    Output, integer D, the day of the nones of the month.
%
  nones = [ ...
    5, 5, 7, 5, 7, 5, 7, 5, 5, 7, 5, 5 ];

  if ( m < 1 || 12 < m )
    d = -1;
  else
    d = nones(m);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Test file for chebtech/trigcoeffs.m

function pass = test_trigcoeffs(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Check a few simple examples.

    f = testclass.make(@(x) zeros(size(x)), [], pref);
    p = trigcoeffs(f);
    pass(n,1) = (norm(p, inf) <= 10*vscale(f)*eps);

    f = testclass.make(@(x) 3*ones(size(x)), [], pref);
    p = trigcoeffs(f);
    pass(n,2) = (norm(p - 3, inf) < 10*vscale(f)*eps);

    f = testclass.make(@(x) 1+cos(pi*x), [], pref);
    p = trigcoeffs(f,3);
    pass(n,3) = (norm(p - [0.5 1 0.5]', inf) < 10*vscale(f)*eps);
    p = trigcoeffs(f,5);
    pass(n,4) = (norm(p - [0 0.5 1 0.5 0]', inf) < 10*vscale(f)*eps);
    p = trigcoeffs(f,1);
    pass(n,5) = (norm(p - 1, inf) < 10*vscale(f)*eps);

    f = testclass.make(@(x) 1 + exp(2*1i*pi*x) + exp(-1i*pi*x), [], pref);
    p = trigcoeffs(f,5);
    pass(n,6) = (norm(p - [0 1 1 0 1 ]', inf) ...
        < 10*vscale(f)*eps);
    p = trigcoeffs(f,9);
    pass(n,7) = (norm(p - [0 0 0 1 1 0 1 0 0 ]', inf) ...
        < 10*vscale(f)*eps);
    p = trigcoeffs(f,3);
    pass(n,8) = (norm(p - [1 1 0]', inf) ...
        < 10*vscale(f)*eps);

    %%
    % Verify operation for array-valued chebtech objects.

    f = testclass.make(@(x) [3*ones(size(x)), 1+cos(pi*x), ...
        1 + exp(2*1i*pi*x) + exp(-1i*pi*x)], [], pref);
    p = trigcoeffs(f,5);
    p_exact = [0 0   0;...
               0 0.5 1;...
               3 1   1;...
               0 0.5 0;...
               0 0   1];
    pass(n,9) = (norm(p(:) - p_exact(:), inf) < 10*max(vscale(f)*eps));

    p = trigcoeffs(f,7);
    p_exact = [0 0   0;...
               0 0   0;...
               0 0.5 1;...
               3 1   1;...
               0 0.5 0;...
               0 0   1;...
               0 0   0];
    pass(n,10) = (norm(p(:) - p_exact(:), inf) < 10*max(vscale(f)*eps));

    p = trigcoeffs(f,3);
    p_exact = [0 0.5 1;...
               3 1   1;...
               0 0.5 0];
    pass(n,11) = (norm(p(:) - p_exact(:), inf) < 10*max(vscale(f)*eps));

    p = trigcoeffs(f,0);
    pass(n,12) = isempty(p);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This folder contains scripts/utilities to support the functionality of classifiers
%
%
%CLSUTIL_SHRINKAGE - Determine shrinkage parameter of a classifier.
%                    Example: the degree of regularization of the estimated covariance matrix
%                    towards a unit sphere (used in LDA classifiers)

author: bbci


text:
function pde = Kelloggdata
%%  KELLOGGDATA data of Kellogg Problem
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',0,'exactu',@exactu,'g_D',@exactu,'Du',@Du,'d',@d);

    function K =  d(p)  % Diffusion constant
    idx = (p(:,1).*p(:,2) >0);
    K  = ones(size(p,1),1);
    K(idx) = 161.4476387975881;
    end

    function u =  exactu(p) % exact solution
    gamma = 0.1;
    sigma = -14.9225565104455152;
    rho = pi/4;
    r = sqrt(sum(p.^2,2));
    theta = atan2(p(:,2),p(:,1));
    theta = (theta>= 0).*theta + (theta<0).*(theta+2*pi);
    mu = (theta>= 0 & theta<pi/2).*cos((pi/2-sigma)*gamma).*cos((theta-pi/2+rho)*gamma)...
        +(theta>= pi/2 & theta<pi).*cos(rho*gamma).*cos((theta-pi+sigma)*gamma)...
        +(theta>= pi & theta<1.5*pi).*cos(sigma*gamma).*cos((theta-pi-rho)*gamma)...
        +(theta>= 1.5*pi & theta<2*pi).*cos((pi/2-rho)*gamma).*cos((theta-1.5*pi-sigma)*gamma);
    u =  r.^gamma.*mu;
    end

    function uprime =  Du(p)
    % the gradient of exact solution
    gamma = 0.1;
    sigma = -14.92256510455152;
    rho = pi/4;
    theta = atan2(p(:,2),p(:,1));  % jiaodu
    theta = (theta>= 0).*theta +(theta<0).*(theta+2*pi);
    t = 1+(p(:,2).^2)./(p(:,1).^2);
    r = sqrt(sum(p.^2,2));
    rg = r.^gamma;

    ux1 = (p(:,1)>= 0.0 & p(:,2)>= 0.0).*(rg.*gamma./r.*cos((pi/2-sigma)*gamma)./r.*p(:,1)...
        .*cos((theta-pi/2+rho)*gamma)...
        +rg.*cos((pi/2-sigma)*gamma).*sin((theta-pi/2+rho)*gamma)...
        *gamma.*p(:,2)./(p(:,1).^2)./t);

    uy1 = (p(:,1)>= 0.0 & p(:,2)>= 0.0).*(rg*gamma./r.*cos((pi/2-sigma).*gamma)...
        .*cos((theta-pi/2+rho).*gamma)./r.*p(:,2)...
        -rg.*cos((pi/2-sigma).*gamma).*sin((theta-pi/2+rho)*gamma)...
        *gamma./p(:,1)./t);

    ux2 = (p(:,1)<= 0.0 & p(:,2)>= 0.0).*(r.^(-1.9).*p(:,1)*gamma.*cos(rho.*gamma)...
        .*cos((theta-pi+sigma).*gamma)...
        +rg.*cos(rho*gamma).*sin((theta-pi+sigma).*gamma)*gamma...
        .*p(:,2)./(p(:,1).^2)./t);
    uy2 = (p(:,1)<= 0.0 & p(:,2)>= 0.0).*( r.^(-1.9).*p(:,2)*gamma.*cos(rho.*gamma)...
        .*cos((theta-pi+sigma)*gamma)...
        -rg.*cos(rho*gamma).*sin((theta-pi+sigma).*gamma)*gamma./p(:,1)./t);

    ux3 = (p(:,1)<= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,1).*gamma.*cos(sigma.*gamma)...
        .*cos((theta-pi-rho).*gamma) ...
        +rg.*cos(sigma.*gamma).*sin((theta-pi-rho).*gamma).*gamma...
        .*p(:,2)./(p(:,1).^2)./t);
    uy3 = (p(:,1)<= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,2)*gamma.*cos(sigma*gamma)...
        .*cos((theta-pi-rho).*gamma) ...
        -rg.*cos(sigma*gamma).*sin((theta-pi-rho).*gamma)*gamma./p(:,1)./t);

    ux4 = (p(:,1)>= 0.0& p(:,2)<= 0.0).*(r.^(-1.9).*p(:,1).*gamma.*cos((pi/2-rho)*gamma)...
        .*cos((theta-3*pi/2-sigma)*gamma) ...
         +rg.*cos((pi/2-rho).*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
         *gamma.*p(:,2)./(p(:,1).^2)./t);

    uy4 = (p(:,1)>= 0.0 & p(:,2)<= 0.0).*(r.^(-1.9).*p(:,2)*gamma.*cos((pi/2-rho)*gamma)...
        .*cos((theta-3*pi/2-sigma)*gamma)...
        -rg.*cos((pi/2-rho)*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
        *gamma./p(:,1)./t);

    uprime(:,1) =  ux1+ux2+ux3+ux4;
    uprime(:,2) =  uy1+uy2+uy3+uy4;
    end
end


text:
function egrad = getEuclideanGradient(problem, x, storedb, key)
% Computes the Euclidean gradient of the cost function at x.
%
% function egrad = getEuclideanGradient(problem, x)
% function egrad = getEuclideanGradient(problem, x, storedb)
% function egrad = getEuclideanGradient(problem, x, storedb, key)
%
% Returns the Euclidean gradient at x of the cost function described in the
% problem structure.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% Because computing the Hessian based on the Euclidean Hessian will require
% the Euclidean gradient every time, to avoid overly redundant
% computations, if the egrad function does not use the store caching
% capabilites, this implements an automatic caching functionality. Writing
% egrad to accept the optional store or storedb parameter will disable
% automatic caching, but allow user controlled caching.
%
% See also: getGradient canGetGradient canGetEuclideanGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 9, 2013.
% 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, 'egrad')
    %% Compute the Euclidean gradient using egrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.egrad)
            case 1
                % If it does not want to deal with the store structure,
                % then we do some caching of our own. There is a small
                % performance hit for this is some cases, but we expect
                % that this is most often the preferred choice.
                store = storedb.get(key);
                if ~isfield(store, 'egrad__')
                    store.egrad__ = problem.egrad(x);
                    storedb.set(store, key);
                end
                egrad = store.egrad__;
            case 2
                % Obtain, pass along, and save the store for x.
                % If the user deals with the store structure, then we don't
                % do any automatic caching: the user is in control.
                store = storedb.getWithShared(key);
                [egrad, store] = problem.egrad(x, store);
                storedb.setWithShared(store, key);
            case 3
                % Pass along the whole storedb (by reference), with key.
                % Same here: no automatic caching.
                egrad = problem.egrad(x, storedb, key);
            otherwise
                up = MException('manopt:getEuclideanGradient:badegrad', ...
                    'egrad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

    else
    %% Abandon computing the Euclidean gradient

        up = MException('manopt:getEuclideanGradient:fail', ...
            ['The problem description is not explicit enough to ' ...
             'compute the Euclidean gradient of the cost.']);
        throw(up);

    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   FANUC MATE M710iC/50, FANUC Robotics Europe.
%
%   Authors: Carlos Carrazoni Perez, Juan Jose Perez Hernandez & David
%   Martinez Pascual
%   Universidad Miguel Hernandez de Elche.
%   date:   7/1/2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2019, by  Carlos Carrazoni Perez, Juan Jose Perez Hernandez
% & David Martinez Pascual
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function robot = parameters()

robot.name= 'Fanuc_M_710iC_50';

%Path where everything is stored for this robot
robot.path = 'robots/fanuc/M_710iC_50';
%Tabla de D-H
robot.DH.theta= '[q(1)+pi/2  q(2)+pi/2  q(3)  q(4)  q(5)  q(6)]';
robot.DH.d='[0.565 0 0 1.016 0 0.175]';
robot.DH.a='[0.15 0.87 0.17 0 0 0]';
robot.DH.alpha= '[pi/2  0  pi/2  pi/2  -pi/2  0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_fanuc_m710(robot, T)';

%number of degrees of freedom
robot.DOF = 6;

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-135) deg2rad(90); %Axis 2, minimum, maximum
                deg2rad(-160) deg2rad(280); %Axis 3
                deg2rad(-360) deg2rad(360); %Axis 4
                deg2rad(-125) deg2rad(125); %Axis 5
                deg2rad(-360) deg2rad(360)]; %Axis 6

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(175); %Axis 1, rad/s
                deg2rad(175); %Axis 2, rad/s
                deg2rad(175); %Axis 3, rad/s
                deg2rad(250); %Axis 4, rad/s
                deg2rad(250); %Axis 5, rad/s
                deg2rad(355)];%Axis 6, rad/s
% end effectors maximum velocity
robot.linear_velmax = 1.0; %m/s, not specified
robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.tool_activated=1;

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [0.9 0.9 0]%./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-1.5 1.5 -1.5 1.5 0 2];
%read graphics files
robot = read_graphics(robot);

%robot.tool=load_robot('equipment/end_tools','water_cutter');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[191.865 71.802 111.5 27.55 5.705 0.674];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.109  -0.099  0.053; %(rx, ry, rz) link 1
                      -0.466  -0.001 -0.287; %(rx, ry, rz) link 2
                      -0.086   0.032  0.007;  %(rx, ry, rz) link 3
                       0       0.342  0.011;  %(rx, ry, rz) link 4
                       0      -0.015  0.076 ; %(rx, ry, rz) link 5
                       0       0      0.010];%(rx, ry, rz) link 6

%Inertia matrices of each link with respect to its D-H reference system.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, for each row
robot.dynamics.Inertia=[8.411      9.434	9.499   -2.810	0.691	1.396;
                        6.297     28.162	22.396  -0.049	-0.012	-9.398;
                        2.313      3.367	2.986	0.467	0.232	-0.285;
                        4.863	0.072	4.846	-0.002	0	0;
                        0.065	0.058	0.022	0	0.004   0;
                        0.000768	0.000767	0.001368	0	0	0];

%Speed reductor at each joint
robot.motors.G=[10 20 120 200 300 1];

repo: ARTE


text:
function varargout = spm_preproc_run(job,action)
% Segment a bunch of images
% FORMAT spm_preproc_run(job)
% job.channel(n).vols{m}
% job.channel(n).biasreg
% job.channel(n).biasfwhm
% job.channel(n).write
% job.tissue(k).tpm
% job.tissue(k).ngaus
% job.tissue(k).native
% job.tissue(k).warped
% job.warp.mrf
% job.warp.cleanup
% job.warp.affreg
% job.warp.reg
% job.warp.fwhm
% job.warp.samp
% job.warp.write
% job.warp.bb
% job.warp.vox
% job.iterations
% job.alpha
%
% See the batch interface for a description of the fields.
%
% See also spm_preproc8.m amd spm_preproc_write8.m
%__________________________________________________________________________
% Copyright (C) 2008-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_preproc_run.m 7670 2019-10-01 16:55:13Z john $

SVNid = '$Rev: 7670 $';

if nargin == 1, action = 'run'; end

switch lower(action)
    case 'run'
        spm('FnBanner',mfilename,SVNid);
        varargout{1} = run_job(job);
        fprintf('%-40s: %30s\n','Completed',spm('time'))                %-#
    case 'check'
        varargout{1} = check_job(job);
    case 'vfiles'
        varargout{1} = vfiles_job(job);
    case 'vout'
        varargout{1} = vout_job(job);
    otherwise
        error('Unknown argument ("%s").', action);
end

%==========================================================================
% Run
%==========================================================================
function vout = run_job(job)

vout = vout_job(job);
tpm  = strvcat(cat(1,job.tissue(:).tpm));
tpm  = spm_load_priors8(tpm);

if ~isfield(job,'iterations'),   nit              =  1; else
                                 nit              = job.iterations; end
if ~isfield(job,'alpha'),        alpha            = 12; else
                                 alpha            = job.alpha;      end
if ~isfield(job.warp,'fwhm'),    job.warp.fwhm    =  1; end
if ~isfield(job.warp,'bb'),      job.warp.bb      =  NaN(2,3); end
if ~isfield(job.warp,'vox'),     job.warp.vox     =  NaN; end
if ~isfield(job.warp,'cleanup'), job.warp.cleanup =  0; end
if ~isfield(job.warp,'mrf'),     job.warp.mrf     =  0; end

if nit > 1
    orig_priors = tpm;
end

for iter=1:nit
    for subj=1:numel(job.channel(1).vols)
        fprintf('Segment %s\n',spm_file(job.channel(1).vols{subj},...
            'link','spm_image(''display'',''%s'')'));

        images = cell(numel(job.channel),1);
        for n=1:numel(job.channel)
            images{n} = job.channel(n).vols{subj};
        end
        obj.image    = spm_vol(char(images));
        spm_check_orientations(obj.image);

        obj.fwhm     = job.warp.fwhm;
        obj.biasreg  = cat(1,job.channel(:).biasreg);
        obj.biasfwhm = cat(1,job.channel(:).biasfwhm);
        obj.tpm      = tpm;
        obj.lkp      = [];
        if all(isfinite(cat(1,job.tissue.ngaus)))
            for k=1:numel(job.tissue)
                obj.lkp = [obj.lkp ones(1,job.tissue(k).ngaus)*k];
            end
        end
        obj.reg      = job.warp.reg;
        obj.samp     = job.warp.samp;

        if iter==1
            % Initial affine registration.
            Affine  = eye(4);
            if ~isempty(job.warp.affreg)
                if isfield(job.warp,'Affine')
                    Affine = job.warp.Affine;
                else
                    % Sometimes the image origins are poorly specified, in which case it might be worth trying
                    % the centre of the field of view instead. The idea here is to run a coarse registration
                    % using two sets of starting estimates, and pick the one producing the better objective function.

                    % Run using origin at centre of the field of view
                    im1            = obj.image(1);
                    M              = im1.mat;
                    c              = (im1.dim+1)/2;
                    im1.mat(1:3,4) = -M(1:3,1:3)*c(:);
                    [Affine1,ll1]  = spm_maff8(im1,8,(obj.fwhm+1)*16,tpm,[],job.warp.affreg); % Closer to rigid
                    Affine1        = Affine1*(im1.mat/M);

                    % Run using the origin from the header
                    im1            = obj.image(1);
                    [Affine2,ll2]  = spm_maff8(im1,8,(obj.fwhm+1)*16,tpm,[],job.warp.affreg); % Closer to rigid

                    % Pick the result with the best fit and use as starting estimate
                    if ll1>ll2
                        Affine  = Affine1;
                    else
                        Affine  = Affine2;
                    end
                end
                Affine = spm_maff8(obj.image(1),job.warp.samp,(obj.fwhm+1)*16,tpm,Affine,job.warp.affreg); % Closer to rigid
                Affine = spm_maff8(obj.image(1),job.warp.samp, obj.fwhm,      tpm,Affine,job.warp.affreg);
            end
            obj.Affine = Affine;
        else
            % Load results from previous iteration for use with next round of
            % iterations, with the new group-specific tissue probability map.
            [pth,nam] = fileparts(job.channel(1).vols{subj});
            res       = load(fullfile(pth,[nam '_seg8.mat']));
            obj.Affine = res.Affine;
            obj.Twarp  = res.Twarp;
            obj.Tbias  = res.Tbias;
            if ~isempty(obj.lkp)
                obj.mg     = res.mg;
                obj.mn     = res.mn;
                obj.vr     = res.vr;
            end
        end

        % in case masking is needed (e.g. CFM for lesions)
        if isfield(job,'msk')
            obj.msk = job.msk ;
        end

        res = spm_preproc8(obj);

        if ~isfield(job,'savemat') || job.savemat==1
            try
                [pth,nam] = fileparts(job.channel(1).vols{subj});
                save(fullfile(pth,[nam '_seg8.mat']),'-struct','res', spm_get_defaults('mat.format'));
            catch
            end
        end

        if iter==nit
            % Final iteration, so write out the required data.
            tmp1 = [cat(1,job.tissue(:).native) cat(1,job.tissue(:).warped)];
            tmp2 =  cat(1,job.channel(:).write);
            tmp3 = job.warp.write;
            spm_preproc_write8(res,tmp1,tmp2,tmp3,job.warp.mrf,job.warp.cleanup,job.warp.bb,job.warp.vox);
        else
            % Not the final iteration, so compute sufficient statistics for
            % re-estimating the template data.
            N    = numel(job.channel);
            K    = numel(job.tissue);
            [cls,M1] = spm_preproc_write8(res,zeros(K,4),zeros(N,2),[0 0],job.warp.mrf,...
                                          job.warp.cleanup,job.warp.bb,job.warp.vox);

            if subj==1
                % Sufficient statistics for possible generation of group-specific
                % template data.
                SS = zeros([size(cls{1}),numel(cls)],'single');
            end

            for k=1:K
                SS(:,:,:,k) = SS(:,:,:,k) + cls{k};
            end
        end

    end
    if iter<nit && nit>1
         % Treat the tissue probability maps as Dirichlet priors, and compute the
         % MAP estimate of group tissue probability map using the sufficient
         % statistics.
         [x1,x2] = ndgrid(1:size(SS,1),1:size(SS,2));
         for i=1:size(SS,3)
             M  = orig_priors.M\M1;
             y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*i + M(1,4);
             y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*i + M(2,4);
             y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*i + M(3,4);
             b  = spm_sample_priors8(orig_priors,y1,y2,y3);
             msk = (y1<1) | (y1>orig_priors.V(1).dim(1)) | ...
                   (y2<1) | (y2>orig_priors.V(1).dim(2)) | ...
                   (y3<1) | (y3>orig_priors.V(1).dim(3));
             for k=1:K
                 bk      = b{k}*alpha;
                 bk(msk) = bk(msk)*0.01;
                 SS(:,:,i,k) = SS(:,:,i,k) + bk;
             end
         end
         save SS.mat SS M1
         tpm.M = M1;
         s     = sum(SS,4);
         for k=1:K
             tmp        = SS(:,:,:,k)./s;
             tpm.bg1(k) = mean(mean(tmp(:,:,1)));
             tpm.bg2(k) = mean(mean(tmp(:,:,end)));
             tpm.dat{k} = spm_bsplinc(log(tmp+tpm.tiny),[ones(1,3)*(tpm.deg-1)  0 0 0]);
         end
    end
end

%==========================================================================
% Check
%==========================================================================
function msg = check_job(job)
msg = {};
if numel(job.channel) >1
    k = numel(job.channel(1).vols);
    for i=2:numel(job.channel)
        if numel(job.channel(i).vols)~=k
            msg = {['Incompatible number of images in channel ' num2str(i)]};
            break
        end
    end
elseif numel(job.channel)==0
    msg = {'No data'};
end

%==========================================================================
% Vout
%==========================================================================
function vout = vout_job(job)

n     = numel(job.channel(1).vols);
parts = cell(n,4);

channel = struct('biasfield',{},'biascorr',{});
for i=1:numel(job.channel)
    for j=1:n
        [parts{j,:}] = spm_fileparts(job.channel(i).vols{j});
    end
    if job.channel(i).write(1)
        channel(i).biasfield = cell(n,1);
        for j=1:n
            channel(i).biasfield{j} = fullfile(parts{j,1},['BiasField_',parts{j,2},'.nii']);
        end
    end
    if job.channel(i).write(2)
        channel(i).biascorr = cell(n,1);
        for j=1:n
            channel(i).biascorr{j} = fullfile(parts{j,1},['m',parts{j,2},'.nii']);
        end
    end
end

for j=1:n
    [parts{j,:}] = spm_fileparts(job.channel(1).vols{j});
end
param = cell(n,1);
for j=1:n
    param{j} = fullfile(parts{j,1},[parts{j,2},'_seg8.mat']);
end

tiss = struct('c',{},'rc',{},'wc',{},'mwc',{});
for i=1:numel(job.tissue)
    if job.tissue(i).native(1)
        tiss(i).c = cell(n,1);
        for j=1:n
            tiss(i).c{j} = fullfile(parts{j,1},['c',num2str(i),parts{j,2},'.nii']);
        end
    end
    if job.tissue(i).native(2)
        tiss(i).rc = cell(n,1);
        for j=1:n
            tiss(i).rc{j} = fullfile(parts{j,1},['rc',num2str(i),parts{j,2},'.nii']);
        end
    end
    if job.tissue(i).warped(1)
        tiss(i).wc = cell(n,1);
        for j=1:n
            tiss(i).wc{j} = fullfile(parts{j,1},['wc',num2str(i),parts{j,2},'.nii']);
        end
    end
    if job.tissue(i).warped(2)
        tiss(i).mwc = cell(n,1);
        for j=1:n
            tiss(i).mwc{j} = fullfile(parts{j,1},['mwc',num2str(i),parts{j,2},'.nii']);
        end
    end
end

if job.warp.write(1)
    invdef = cell(n,1);
    for j=1:n
        invdef{j} = fullfile(parts{j,1},['iy_',parts{j,2},'.nii']);
    end
else
    invdef = {};
end

if job.warp.write(2)
    fordef = cell(n,1);
    for j=1:n
        fordef{j} = fullfile(parts{j,1},['y_',parts{j,2},'.nii']);
    end
else
    fordef = {};
end

vout  = struct('channel',channel,'tiss',tiss,'param',{param},'invdef',{invdef},'fordef',{fordef});

%==========================================================================
% Vfiles
%==========================================================================
function vf = vfiles_job(job)
vout = vout_job(job);
vf   = vout.param;
if ~isempty(vout.invdef), vf = {vf{:}, vout.invdef{:}}; end
if ~isempty(vout.fordef), vf = {vf{:}, vout.fordef{:}}; end
for i=1:numel(vout.channel)
    if ~isempty(vout.channel(i).biasfield), vf = {vf{:}, vout.channel(i).biasfield{:}}; end
    if ~isempty(vout.channel(i).biascorr),  vf = {vf{:}, vout.channel(i).biascorr{:}};  end
end

for i=1:numel(vout.tiss)
    if ~isempty(vout.tiss(i).c),   vf = {vf{:}, vout.tiss(i).c{:}};   end
    if ~isempty(vout.tiss(i).rc),  vf = {vf{:}, vout.tiss(i).rc{:}};  end
    if ~isempty(vout.tiss(i).wc),  vf = {vf{:}, vout.tiss(i).wc{:}};  end
    if ~isempty(vout.tiss(i).mwc), vf = {vf{:}, vout.tiss(i).mwc{:}}; end
end
vf = reshape(vf,numel(vf),1);

author: spm


text:
function  [X] =  WNNM( Y, C, NSig, m, Iter )
    [U,SigmaY,V] =   svd(full(Y),'econ');
    PatNum       = size(Y,2);
    TempC  = C*sqrt(PatNum)*2*NSig^2;
    [SigmaX,svp] = ClosedWNNM(SigmaY,TempC,eps);
    X =  U(:,1:svp)*diag(SigmaX)*V(:,1:svp)' + m;
return;

author: lbasek

repo: image-denoising-benchmark


text:
y       = load('data/motorcycle.dat')';
time    = y(1, :);
delta   = y(2, [2:end 1]);
y       = y(3, :);

%% Spline smoothing model %%
spline              = ssm_spline(delta);
spline              = estimate(y, spline, [1 0.1]);
[alphahat V]        = statesmo(y, spline);
conf                = squeeze(1.96*realsqrt(V(1, 1, :)))';
[eps eta epsvar]    = disturbsmo(y, spline);

figure('Name', 'Motorcycle acceleration data analyzed by a cubic spline');
subplot(2, 1, 1), plot(time, alphahat(1, :), 'b'), hold all, plot(time, [alphahat(1, :)+conf; alphahat(1, :)-conf], 'b:'), scatter(time, y, 10, 'r', 's', 'filled'), hold off, title('Spline and 95% confidence intervals'), ylim([-140 80]), set(gca,'YGrid','on');
subplot(2, 1, 2), scatter(time, eps./realsqrt(epsvar), 10, 'r', 's', 'filled'), title('Standardized irregular'), set(gca,'YGrid','on');
if ispc, set(gcf, 'WindowStyle', 'docked'); end

%% Continuous local level model %%
abseps          = abs(eps);
contllm         = [ssm_gaussian ssmodel('continuous local level', 0, 1, 1, 1, ssmat(0, [], true, zeros(size(delta)), true), 'Qd', {@(X) exp(2*X)*delta}, {[]}, ssparam({'zeta var'}, '1/2 log'))];
[contllm logL]  = estimate(abseps, contllm, [1 0.1]);
alphahat        = statesmo(abseps, contllm);

figure('Name', 'Correction for heteroscedasticity');
subplot(3, 1, 1), plot(time, alphahat, 'b'), hold all, scatter(time, abseps, 10, 'r', 's', 'filled'), hold off, title('Absolute smoothed irregular and h^\ast_t'), ylim([0 87.5]);

%% Correction for heteroscedasticity %%
h2                  = (alphahat/alphahat(1)).^2;
splineh             = [ssmodel('Heteroscedastic noise', ssmat(0, [], true, zeros(size(h2)), true), zeros(1, 0), [], [], [], 'Hd', {@(X) exp(2*X)*h2}, {[]}, ssparam({'epsilon var'}, '1/2 log')) spline];
splineh             = estimate(y, splineh, [1 0.1]);
[alphahath Vh]      = statesmo(y, splineh);
confh               = squeeze(1.96*realsqrt(Vh(1, 1, :)))';
[epsh eta epsvarh]  = disturbsmo(y, splineh);

subplot(3, 1, 2), plot(time, alphahath(1, :), 'b'), hold all, plot(time, [alphahath(1, :)+confh; alphahath(1, :)-confh], 'b:'), scatter(time, y, 10, 'r', 's', 'filled'), hold off, title('Spline and 95% confidence intervals'), ylim([-140 80]), set(gca,'YGrid','on');
subplot(3, 1, 3), scatter(time, epsh./realsqrt(epsvarh), 10, 'r', 's', 'filled'), title('Standardized irregular'), set(gca,'YGrid','on');
if ispc, set(gcf, 'WindowStyle', 'docked'); end

repo: REST


text:
function stroud_test27 ( )

%*****************************************************************************80
%
%% TEST27 tests SIMPLEX_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_ND_INDEX;

  num = function_nd_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST27\n' );
  fprintf ( 1, '  SIMPLEX_ND approximates integrals inside an\n' );
  fprintf ( 1, '    arbitrary simplex in ND.\n' );
  fprintf ( 1, '\n' );

  for n = 2 : 4

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Spatial dimension N = %d\n', n );
%
%  Restore values of simplex.
%
    v = setsim ( n );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Simplex vertices:\n' );
    fprintf ( 1, '\n' );

    for i = 1 : n+1
      for j = 1 : n
        fprintf ( 1, '  %4f', v(i,j) );
      end
      fprintf ( 1, '\n' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  F(X)    SIMPLEX_ND\n' );
    fprintf ( 1, '\n' );

    for i = 1 : num

      FUNC_ND_INDEX = i;

      result = simplex_nd ( 'function_nd', n, v );

      fname = function_nd_name ( i );

      fprintf ( 1, '  %s  %14f\n', fname, result );

      v = setsim ( n );

    end

  end

  return
end


text:
function  generateProposalsForImages()
%boxes have this format: [xmin ymin xmax ymax]
addpath(genpath('../'));
addpath(genpath('../../../utils/'));
config=createConfig();

imageLoc=config.path.input;
outputLoc=config.path.output;
%extension inclues .
imageExt=config.opts.imageExt;
images=dir([imageLoc '*' imageExt]);

for i=1:length(images)
	imName=images(i).name;
	image=imread([imageLoc imName]);
	[ranked_regions superpixels image_data]=generate_proposals(image);
	boxes=zeros(length(ranked_regions),4);
	for j=1:length(ranked_regions)
		boxes(j,:)=mask2box(ismember(superpixels, ranked_regions{j}));
	end

	if(isfield(config.opts,'numProposals'))
        numProposals=config.opts.numProposals;

        if(size(boxes,1)>=numProposals)
                boxes=boxes(1:numProposals,:);
                ranked_regions=ranked_regions(1:numProposals);
        else
                fprintf('Only %d proposals were generated for image:%s\n',size(boxes,1),imName);
        end
end
	%mask2box outputs  [ymin xmin ymax xmax]
	boxes=[boxes(:,2) boxes(:,1) boxes(:,4) boxes(:,3)];
	proposals.boxes=boxes;
	proposals.scores = [size(boxes,1):-1:1]`;
	proposals.regions.ranked_regions=ranked_regions;
	proposals.regions.superpixels=superpixels;
	proposals.regions.image_data=image_data;

	proposalFileName=strrep(imName,imageExt,'.mat');
	save([outputLoc proposalFileName],'proposals');
end

author: Cloud-CV

repo: object-proposals


text:
function [meanError, all_rot_preds, all_rot_gts, meanErrors, all_errors, rels_all, seq_ids] = calcIctError(resDir, gtDir)
%CALCICTERROR Summary of this function goes here
%   Detailed explanation goes here

    polhemus = 'polhemusNorm.csv';

    sequences = dir([resDir '*.csv']);

    rotMeanErr = zeros(numel(sequences),3);
    rotRMS = zeros(numel(sequences),3);
    rot = cell(1,numel(sequences));
    rotg = cell(1,numel(sequences));

    rels_all = [];

    seq_ids = {};

    for i = 1:numel(sequences)

        [~, name,~] = fileparts(sequences(i).name);

        fname = [resDir '/' sequences(i).name];
        if(i == 1)
            % First read in the column names
            tab = readtable(fname);
            column_names = tab.Properties.VariableNames;

            confidence_id = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'confidence'));
            rot_ids = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'pose_R'));
        end

        all_params  = dlmread(fname, ',', 1, 0);

        rot{i} = all_params(:, rot_ids);
        rels = all_params(:, confidence_id);

        % the reliabilities of head pose
        rels_all = cat(1, rels_all, rels);
        [txg tyg tzg rxg ryg rzg] =  textread([gtDir name '/'  polhemus], '%f,%f,%f,%f,%f,%f');

        rotg{i} = [rxg ryg rzg];

        % Correct the first frame so it corresponds to (0,0,0), as slightly
        % different pose might be assumed frontal and this corrects for
        % that

        % Work out the correction matrix for ground truth
        rot_corr_gt = Euler2Rot(rotg{i}(1,:));
        for r_e = 1:size(rotg{i},1)
            rot_curr_gt = Euler2Rot(rotg{i}(r_e,:));
            rot_new_gt = rot_corr_gt' * rot_curr_gt;
            rotg{i}(r_e,:) = Rot2Euler(rot_new_gt);
        end

        % Work out the correction matrix for estimates
        rot_corr_est = Euler2Rot(rot{i}(1,:));
        for r_e = 1:size(rot{i},1)
            rot_curr_est = Euler2Rot(rot{i}(r_e,:));
            rot_new_est = rot_corr_est' * rot_curr_est;
            rot{i}(r_e,:) = Rot2Euler(rot_new_est);
        end

        % Convert the ground truth and estimates to degrees
        rot{i} = rot{i} * (180/ pi);
        rotg{i} = rotg{i} * (180/ pi);

        % Now compute the errors
        rotMeanErr(i,:) = mean(abs((rot{i}(:,:)-rotg{i}(:,:))));
        rotRMS(i,:) = sqrt(mean(((rot{i}(:,:)-rotg{i}(:,:))).^2));

        seq_ids = cat(1, seq_ids, repmat({[name 'ict']}, size(rot{i},1), 1));

    end
    allRot = cell2mat(rot');
    allRotg = cell2mat(rotg');
    meanErrors = rotMeanErr;
    meanError = mean(abs((allRot(:,:)-allRotg(:,:))));
    all_errors = abs(allRot-allRotg);
    rmsError = sqrt(mean(((allRot(:,:)-allRotg(:,:))).^2));
    errorVariance = var(abs((allRot(:,:)-allRotg(:,:))));

    all_rot_preds = allRot;
    all_rot_gts = allRotg;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [xyVec] = createXYCoordVec(start,pix,dim,dir)
% createXYCoordVec
%
% Create one dimensional vector

if strcmpi(dir,'pos')% positive direction
    for i=1:dim
        xyVec(i)=start+(pix.*i)-pix;
    end

elseif strcmpi(dir,'neg')% negative direction
    for i=1:dim
        xyVec(i)=start-(pix.*i)+pix;
    end
end

author: cerr

repo: CERR