text:
% dtiComputePlotDiffusionPropertiesAlongFiberGroup

% OVERVIEW:
% This script takes a group of subjects with a given fiber group(s) and set
% of ROIs and:
%   1. Loads the fibers and rois, as well as the dt6.mat file created
%      during preprocessing.
%   2. Computes diffusion properties along the fiber group using
%      dtiComputeDiffusionPropertiesAlongFG. (see that function for more
%      info.
%   3. Saves out a structure (dtval) with the tensor values for 'fa','md',
%      'rd', and 'ad'.
%   4. Plots the data
%   5. Saves the images to "saveDir"

% VARIABLES:
%   saveDir = the directory wherein the data and figures will be saved.
%   subs = the indetifying codes for each subject in your directory.
%         (e.g., {subs = 'sub1','sub2','sub3'...})
%   baseDir = the level of your directory structure that has your data at
%             the first level.
%   dirs = the name of the directory that contains each subject's dt6.mat
%          file - usually named for the number of diffusion directions.
%   fgName = cell array containg the names of each fiber group you wish to
%            compute properties for.
%   rois = cell array containing the names of the first of 2 rois used to
%          generate the fibers in fgName. NO FILE EXTENSIONS
%   roi2 = cell array containing the names of the second of 2 rois used to
%          generate the fibers in fgName. NO FILE EXTENTIONS
%   plotNames = a cell array contining one name for each fiber group, which
%               will be the name for each plot associated with that fiber
%               group.
%   numberOfNodes = the number of steps along the fg at which properties
%                   will be computed.
%   propertyofinterest = the properties that will be computed for each
%                        group.

% NOTES:
%  You will have to set variables in the first two sections. If  you have a
%  probelm with data not being found check the directory structure within the
%  loop.

%  You must provide the ROIs in rois and rois2 that were used to track the
%  fibers in fgName.
%   For example, if you used roiA and roiB to track fibergroup AB in
%   contrack, then the cell arrays should look like:
%       fgName  = {'fibergroupAB.pdb'};
%       rois    = {'roiA'};
%       rois2   = {'roiB'};

%  The name in the rois cell array will be used for naming the figures.
%  This can be changed in sections IV - VI if you don't wish to use this
%  default behavior.

% HISTORY:
% 2011.1.10 - LMP wrote the code.

%% I. Set directories, rois and plot names

saveDir = '/path/to/your/data/directory';
           if ~exist(saveDir,'file'), mkdir(saveDir); end

subs    = {'sub1','sub2','sub3'};
baseDir = '/path/to/your/subjects"/data/';
dirs    = 'dti06trilinrt';

fgName  = {'fibergroupAB.pdb','fibergroupCD.pdb'};
rois    = {'roiA','roiC'};
rois2   = {'roiB','roiD'};

plotName = {'Fiber Group AB Plot Name','Fiber Group CD Plot Name'};

%% II. Set up parameters

numberOfNodes      = 50;
propertyofinterest = {'fa','md', 'rd', 'ad'};
roi1name           = rois;
roi2name           = rois2;

%% III Loop over subjects and compute fiber properties
wb = mrvWaitbar(0,'Overall Script Progress');
for zz = 1:numel(propertyofinterest)
    for kk = 1:numel(fgName)
        for ii=1:numel(subs)
            sub      = dir(fullfile(baseDir,[subs{ii} '*']));
            subDir   = fullfile(baseDir,sub.name);
            dt6Dir   = fullfile(subDir,dirs);
            fiberDir = fullfile(dt6Dir,'fibers','conTrack');
            roiDir   = fullfile(dt6Dir,'ROIs');
            try
                % III. 1 LOAD THE DATA
                fibersFile = fullfile(fiberDir,fgName{kk});
                fg         = mtrImportFibers(fibersFile);
                roi1File   = fullfile(roiDir, [roi1name{kk} '.mat']);
                roi2File   = fullfile(roiDir, [roi2name{kk} '.mat']);
                roi1       = dtiReadRoi(roi1File);
                roi2       = dtiReadRoi(roi2File);
                dt         = dtiLoadDt6(fullfile(dt6Dir,'dt6.mat'));

                % III. 2 Compute
                [fa(:, ii),md(:, ii),rd(:, ii),ad(:, ii)]=...
                    dtiComputeDiffusionPropertiesAlongFG(fg, dt, roi1, roi2, numberOfNodes);
            catch ME
                disp(ME.message);
                disp(['PROBLEM loading data for subject: ' subs{ii} '. This subject will NOT be included in the graph!']);
            end
        end
        dtval.(propertyofinterest{zz}).(rois{kk}) = eval(propertyofinterest{zz});

    end
mrvWaitbar(zz/(numel(propertyofinterest)),wb);

end
save((fullfile(saveDir,'dtval.mat')),'dtval');

%% PLOT DATA

%% IV Plot results: Average for each roi.

load(fullfile(saveDir,'dtval.mat'));

    for zz = 1:numel(propertyofinterest) % fa,md,rd,ad
        for kk = 1:numel(rois)

            AVG = mean(dtval.(propertyofinterest{zz}).(rois{kk}),2);

            figure;
            hold on
            plot(AVG,'b','LineWidth',5);
            set(gca,'PlotBoxAspectRatio',[1,.7,1]);

           switch propertyofinterest{zz}
                case 'fa'
                    yText = 'FA (weighted)';
                    titleText = ['Fractional Anisotropy Along ' plotName{kk}];
                case 'md'
                    yText = 'MD \mum^2/msec (weighted)';
                    titleText = ['Mean Diffusivity Along ' plotName{kk}];
                case 'rd'
                    yText = 'RD \mum^2/msec (weighted)';
                    titleText = ['Radial Diffusivity Along ' plotName{kk}];
                case 'ad'
                    yText = 'AD \mum^2/msec (weighted)';
                    titleText = ['Axial Diffusivity Along ' plotName{kk}];
            end

            title(titleText);
            ylabel(yText);
            xlabel('Fiber Group Trajectory');
            set(gca,'xtick',[0 numberOfNodes]);
            set(gca,'xticklabel',{rois{kk},rois2{kk}})
            saveName = [plotName{kk} '_' propertyofinterest{zz}];
            saveas(gcf,(fullfile(saveDir,saveName)),'epsc2');
            hold off
        end
    end

%% V Plot results: Each fiber group on the same grapgh.

load(fullfile(saveDir,'dtval.mat'));
col = jet(numel(rois)); % set colors for graph.

    for zz = 1:numel(propertyofinterest) % fa,md,rd,ad
        figure;
        for kk = 1:numel(rois)

            AVG = mean(dtval.(propertyofinterest{zz}).(rois{kk}),2);

            hold on
            plot(AVG,'LineWidth',5,'color',col(kk,:));
            set(gca,'PlotBoxAspectRatio',[1,.7,1]);

            switch propertyofinterest{zz}
                case 'fa'
                    yText = 'FA (weighted)';
                    titleText = ['Fractional Anisotropy'];
                case 'md'
                    yText = 'MD \mum^2/msec (weighted)';
                    titleText = ['Mean Diffusivity'];
                case 'rd'
                    yText = 'RD \mum^2/msec (weighted)';
                    titleText = ['Radial Diffusivity'];
                case 'ad'
                    yText = 'AD \mum^2/msec (weighted)';
                    titleText = ['Axial Diffusivity '];
            end

            title(titleText);
            ylabel(yText);
            xlabel('Fiber Group Trajectory');

        end
            set(gca,'xtick',[0 numberOfNodes]);
            set(gca,'xticklabel',{rois{kk},rois2{kk}})
            ld = legend(plotName{:});
            set(ld,'Interpreter','tex','Location','NorthEast');
            saveName = [propertyofinterest{zz} '_allGroups'];
            saveas(gcf,(fullfile(saveDir,saveName)),'epsc2');
    end

%% VI Plot results for each roi sorted by subject

load(fullfile(saveDir,'dtval.mat'));

col = jet(numel(subs)); % set colors for graph.

for zz = 1:numel(propertyofinterest) % fa,md,rd,ad
    for kk = 1:numel(rois)

        AVG = mean(dtval.(propertyofinterest{zz}).(rois{kk}),2);

        figure;
        hold on
        % Loop over each subject and plot their data using a different color line (col)
        for ss = 1:numel(subs)
            plot(dtval.(propertyofinterest{zz}).(rois{kk})(:,ss),'color',col(ss,:));
        end
        ld = legend(subs{:});
        set(ld,'Interpreter','tex','Location','NorthEastOutSide');
        plot(AVG,'g','LineWidth',5);
        set(gca,'PlotBoxAspectRatio',[1,.7,1]);

        switch propertyofinterest{zz}
            case 'fa'
                yText = 'FA (weighted)';
                titleText = ['Fractional Anisotropy Along ' plotName{kk}];
            case 'md'
                yText = 'MD \mum^2/msec (weighted)';
                titleText = ['Mean Diffusivity Along ' plotName{kk}];
            case 'rd'
                yText = 'RD \mum^2/msec (weighted)';
                titleText = ['Radial Diffusivity Along ' plotName{kk}];
            case 'ad'
                yText = 'AD \mum^2/msec (weighted)';
                titleText = ['Axial Diffusivity Along ' plotName{kk}];
        end

        title(titleText);
        ylabel(yText);
        xlabel('Fiber Group Trajectory');
        set(gca,'xtick',[0 numberOfNodes]);
        set(gca,'xticklabel',{rois{kk},rois2{kk}});
        whitebg;
        saveName = [plotName{kk} '_' propertyofinterest{zz} '_SUBJECT_avgLine'];
        saveas(gcf,(fullfile(saveDir,saveName)),'epsc2');
    end
end
close(wb)
close all

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [n,f,a,b]=lpcar2fm(ar,t)
%LPCAR2RF Convert autoregressive coefficients to formant freq+amp+bw [N,F,A,B]=(AR,T)
%
% Input:   ar(:,p+1)  Autoregressive coefficients
%          t          Threshold (see below)
% Output:  n          Number of formants found
%          f          Formant frequencies in normalized Hz (in increasing order)
%          a          Formant amplitudes
%          b          Formant bandwidths in normalized Hz
%
% The number of columns in the output arrays f, a and b is max(n); surplus positions
% in any given row have f=b=0.
%
% In determining formants, poles are ignored if any of the following hold:
%        (a) they are on the real axis
%        (b) they have bandwidth > t*frequency (if t>0)
%        (c) they have bandwidth > -t (if t<=0)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcar2fm.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
p=p1-1;
d=(1:nf)';
zz=lpcar2zz(ar);
ig=imag(zz)<=0;
n=p1-1-sum(ig,2);
mn=max(n);

% remove redundant columns

if mn<p
   [ig,ix]=sort(ig,2);
   zz=reshape(zz(d(:,ones(1,mn))+nf*(ix(:,1:mn)-1)),nf,mn);
   ig(:,mn+1:end)=[];
end

zz(ig)=1;      % to prevent infinities
f=angle(zz)*0.5/pi;
b=-log(abs(zz))/pi;
if nargin > 1
   if t>0
      ig=ig | b>t*f;
   else
      ig=ig | b+t>0;
   end
end
f(ig)=0;
b(ig)=0;
n=mn-sum(ig,2);
m=max(n);

% remove redundant columns

[igf,ix]=sort(ig+f,2);
dd=d(:,ones(1,m))+nf*(ix(:,1:m)-1);
zz=reshape(zz(dd),nf,m);
f=reshape(f(dd),nf,m);
b=reshape(b(dd),nf,m);
ig=reshape(ig(dd),nf,m);

% now calculate gain
ap=permute(ar,[1 3 2]);
pw=permute(-2*pi*1i*(0:p),[1 3 2]);
a=abs(sum(ap(:,ones(1,m),:).*exp(pw(ones(1,nf),ones(1,m),:).*f(:,:,ones(1,p1))),3)).^(-1);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
%   atlasNew = updateTinC(atlas,X,Y)
%
%Author: B. Fischer
%Purpose:
%   Deform an image, say atlas, so that values in the implicit function,
%   atlas(1:M,1:N) are mapped to the positions in atlas(X,Y).
%   This function is used in the elastic matching algorithms in the
%   Analysis/Atlas directory.
%
%   It was written in C- by Bernd.  It appears that the values in X,Y
%   must match the coordinates (size) of atlas.  No linear interpolation is
%   performed, just a deformation of the cartesian coordinates.
%
%  These notes by BW.
%


text:
function [x,y,z] = getLMpolygon3D(polygon)
%
% Utility function that gives the coordinates of the polygon's vertices

x = str2num(char({polygon.pt.x})); % get X polygon coordinates
y = str2num(char({polygon.pt.y})); % get Y polygon coordinates
z = str2num(char({polygon.pt.z})); % get Y polygon coordinates

author: CSAILVision

repo: LabelMeToolbox


text:
function [dist]=find3DNeighbourDists(mesh,scaling);
% function [dist]=find3DNeighbourDists(mesh,scaling);
% AUTHOR: WADE
% DATE : 062700
% PURPOSE : Return a list of the distances between each node and its neighbours.
% If there are n connections in the mesh, there should be n entries in distSQ
% See also find2DNeighboutDists - does the same thing with the 2d mesh.
% ARW 031201 - Now takes a scaling argument (equiv of dimdist for mrManDist)
% This scales distances along the appropriate matrix dimensions
% ARW 082702 - Now returns the abs distance not the squared one.
if (~exist('scaling'))
	scaling=[1 1 1];
end

spX=mesh.connectionMatrix;
spY=mesh.connectionMatrix;
spZ=mesh.connectionMatrix;
nVerts=length(mesh.connectionMatrix);

% What am I doing here?
 %   S = SPARSE(i,j,s,m,n,nzmax) uses the rows of [i,j,s] to generate an
 %   m-by-n sparse matrix with space allocated for nzmax nonzeros.  The
 %   two integer index vectors, i and j, and the real or complex entries
 %   vector, s, all have the same length, nnz, which is the number of
 %   nonzeros in the resulting sparse matrix S .  Any elements of s
 %   which have duplicate values of i and j are added together.

xI=sparse((1:nVerts),(1:nVerts),mesh.uniqueVertices(:,1),nVerts,nVerts)*scaling(1);
yI=sparse((1:nVerts),(1:nVerts),mesh.uniqueVertices(:,2),nVerts,nVerts)*scaling(2);
zI=sparse((1:nVerts),(1:nVerts),mesh.uniqueVertices(:,3),nVerts,nVerts)*scaling(3);
% So we're creating 3 diagonal sparse matrices.
% They contain the x,y,z coordinates of the uniqueVertices on their diagonals.

spX=spX*xI;
spY=spY*yI;
spZ=spZ*zI;

% What does this mean? If there's a connection at between nodes i and j, there will be a '1' in
% the conmat at i,j (and j,i...).
% The above multiplication means that all the '1's (connections) in column j will be
% replaced by the X, Y, or Z ordinates of the j'th node.

% Cols of sp(X,Y,Z) are now ordinates of corresponding vertices
spX=spX-spX'; % Smart eh? So the i,jth entry is the X ordinate of the j'th node. And the j,i'th entry is the X ordinate ofthe
              % ith node. So this subtraction calculates the X distance between them.
spY=spY-spY';
spZ=spZ-spZ';
dist=sqrt(spX.^2+spY.^2+spZ.^2); % And this just sums the squares of the X,Y and Z distances. Take the sqrt of this matrix to get the
                               % edge distances.


text:
classdef matRad_SquaredDeviation < DoseObjectives.matRad_DoseObjective
% matRad_SquaredDeviation Implements a penalized least squares objective
%   See matRad_DoseObjective for interface description
%
% References
%     -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties (Constant)
        name = 'Squared Deviation';
        parameterNames = {'d^{ref}'};
        parameterTypes = {'dose'};
    end

    properties
        parameters = {60};
        penalty = 1;
    end

    methods
        function obj = matRad_SquaredDeviation(penalty,dRef)

            %If we have a struct in first argument
            if nargin == 1 && isstruct(penalty)
                inputStruct = penalty;
                initFromStruct = true;
            else
                initFromStruct = false;
                inputStruct = [];
            end

            %Call Superclass Constructor (for struct initialization)
            obj@DoseObjectives.matRad_DoseObjective(inputStruct);

            %now handle initialization from other parameters
            if ~initFromStruct
                if nargin == 2 && isscalar(dRef)
                    obj.parameters{1} = dRef;
                end

                if nargin >= 1 && isscalar(penalty)
                    obj.penalty = penalty;
                end
            end
        end

        %% Calculates the Objective Function value
        function fDose = computeDoseObjectiveFunction(obj,dose)
            % deviation : dose minus prefered dose
            deviation = dose - obj.parameters{1};
            % claculate objective function
            fDose = obj.penalty/numel(dose) * (deviation'*deviation);
        end

        %% Calculates the Objective Function gradient
        function fDoseGrad   = computeDoseObjectiveGradient(obj,dose)
            % deviation : Dose minus prefered dose
            deviation = dose - obj.parameters{1};

            % calculate delta
            fDoseGrad = 2 * obj.penalty/numel(dose) * deviation;
        end
    end

end

repo: matRad


text:
function x=randiscr(p,n,a)
%RANDISCR Generate discrete random numbers with specified probabiities [X]=(P,N,A)
%
% Inputs: P  vector of probabilities (not necessarily normalized)
%         N  number of random values to generate [default = 1]
%         A  output alphabet [default = 1:length(p)]
%
% Outputs: X  vector of values taken from alphabet A
%
% The vector P is internally normalized by dividing by its sum.
% If P is an M-dimensional matrix (and A is unspecified), then X will
% have dimensions (N,M) with the corresponding indices for each dimension.

% Somewhat similar in function to RANDSRC in the comms toolbox

%   Copyright (c) 2005 Mike Brookes,  mike.brookes@ic.ac.uk
%      Version: $Id: randiscr.m,v 1.3 2007/05/04 07:01:39 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    n=1;
end
d=length(p(:)); % size of output alphabet
z=zeros(d+n-1,1); % array to hold random numbers
z(1:d)=cumsum(p(:)/sum(p(:))); % last value is actually overwritten in the next line
z(d:d+n-1)=rand(n,1);
[y,iy]=sort(z);
y(iy)=(1:d+n-1)';
m=zeros(d+n-1,1);
m(y(1:d-1))=1;
m(1)=m(1)+1;
mc=cumsum(m);
x=mc(y(d:d+n-1));
% lo=[0;z(1:d-1)];hi=[z(1:d-1);1]; % print boundaries
% [lo(x) z(d:d+n-1) hi(x) x] % print out random values and results
if nargin>2
    x=a(x);
elseif length(p(:))>length(p) % need multiple dimensions
    v=x-1;
    s=cumprod(size(p));
    m=length(s);
    s(2:end)=s(1:end-1);
    s(1)=1;
    x=zeros(n,m);
    for i=m:-1:1
        x(:,i)=1+floor(v/s(i));
        v=rem(v,s(i));
    end
end

author: decouples


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

predictions =  sigmoid(X*theta);

leftPart = -y' * log(predictions);

rightPart = (1 - y') * log(1 - predictions);

J = (1 / m) * (leftPart - rightPart);

grad = (1 / m) * ((predictions - y)' * X);

% =============================================================

end

author: vkosuri

repo: CourseraMachineLearning


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: demo for 2D data setup and visualization
%
% load images                        (hands-?.jpg)
% visualize ij and omega based data  (viewImage)
% creates multi-lebel representation (getMultilevel)
%
%==============================================================================

clear, close all, help(mfilename);

% load data
Tij = imread('hands-T.jpg'); Tdata = double(flipud(Tij))';
Rij = imread('hands-R.jpg'); Rdata = double(flipud(Rij))';

% specify domain Omega = [omega(1),omega(2)]x[omega(3),omega(3)];
omega = [0,20,0,25];
m     = size(Tdata);

% setup image viewer
viewImage('reset','viewImage','viewImage2D','colormap','gray(256)');

% visualize
FAIRfigure(1);  colormap(gray(256));
subplot(2,2,1); imagesc(Tij);              title('original T data, uint8, ij');
subplot(2,2,2); viewImage(Tdata,omega,m);  title('FAIR T data on \Omega, xy');
subplot(2,2,3); imagesc(Rij);              title('original R data, uint8, ij');
subplot(2,2,4); viewImage(Rdata,omega,m);  title('FAIR R data on \Omega, xy');

% create multi-level representaion
ML = getMultilevel({Tdata,Rdata},omega,m,'fig',2);

%==============================================================================


text:
function [y, sigma, p] = linRegPred(model, X, t)
% Compute linear regression model reponse y = w'*X+w0 and likelihood
% Input:
%   model: trained model structure
%   X: d x n testing data
%   t (optional): 1 x n testing response
% Output:
%   y: 1 x n prediction
%   sigma: variance
%   p: 1 x n likelihood of t
% Written by Mo Chen (sth4nth@gmail.com).
w = model.w;
w0 = model.w0;
y = w'*X+w0;
%% probability prediction
if nargout > 1
    beta = model.beta;
    if isfield(model,'U')
        U = model.U;        % 3.54
        Xo = bsxfun(@minus,X,model.xbar);
        XU = U'\Xo;
        sigma = sqrt((1+dot(XU,XU,1))/beta);   % 3.59
    else
        sigma = sqrt(1/beta)*ones(1,size(X,2));
    end
end

if nargin == 3 && nargout == 3
    p = exp(-0.5*(((t-y)./sigma).^2+log(2*pi))-log(sigma));
end

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [x, infos] = mu_conv_nmf(V, rank, t, in_options)
% Multiplicative update (MU) based convolutive non-negative matrix factorization (MU-Conv-NMF).
%
% The problem of interest is defined as
%
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       matrix      V
%       rank        rank
%
% Output:
%       w           solution of w
%       infos       information
%
% References:
%
%
% This file is part of NMFLibrary
%
% This file has been ported by H.Kasai from
% convNMF_MM1.m and convNMF_MM2.m at https://github.com/lyn202206/ADMM-Convolutive-NMF.
%
% Ported by H.Kasai on June 29, 2022
%
% Change log:
%
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.metric_type = 'beta-div';
    local_options.d_beta = 2;
    local_options.alg = 'type1';

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    if ~isfield(options, 'x_init')
        W = zeros(m, rank, t);
        for i = 1 : t
            [init_factors, ~] = generate_init_factors(V, rank, init_options);
            W(:, :, i) = init_factors.W;
        end
        H = init_factors.H;
    else
        W = init_options.x_init.W;
        H = init_options.x_init.H;
    end

    % initialize
    epoch = 0;
    grad_calc_count = 0;

    options = check_divergence(options);
    sub_mode = sprintf('beta=%.1f', options.d_beta);
    if ~strcmp(options.metric_type, 'beta-div')
        sub_mode = options.metric_type;
    end
    method_name = sprintf('MU-Conv (%s)', sub_mode);

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    V_hat = zeros(m, n);
    for i = 0 : t-1
        V_hat = V_hat + W(:, :, i+1) * shift_t(H, i);
    end

    % store initial info
    clear infos;

    [Wcon, Hcon] = reconstruct_wh(W, H, t);
    [infos, f_val, optgap] = store_nmf_info(V, Wcon, Hcon, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('MU-Conv (%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H sequentially
        for i = 1 : n

            num = zeros(rank, 1);
            denom = num;
            if i <= n-t+1
                for n_prime = i : i+t-1
                    num = num + W(:, :, n_prime-i+1)' * ((V(:, n_prime) + eps) .* (V_hat(:, n_prime) + eps) .^ (options.d_beta - 2));
                    denom = denom + W(:, :, n_prime-i+1)'*(V_hat(:, n_prime)+eps).^(options.d_beta-1);
                end
            else
                for n_prime = i : n
                    num = num + W(:, :, n_prime-i+1)' * ((V(:, n_prime) + eps) .* (V_hat(:, n_prime) + eps) .^ (options.d_beta - 2));
                    denom = denom + W(:, :, n_prime-i+1)'*(V_hat(:, n_prime) + eps) .^ (options.d_beta-1);
                end
            end

            if strcmp(options.alg, 'type1')
                h_n_old = H(:, i);
                H(:, i) = H(:, i) .* (num./denom) .^ gamma_beta(options.d_beta);

                if i <= n-t+1
                    for n_prime = i:i+t-1
                        V_hat(:, n_prime) = max(V_hat(:, n_prime) + W(:, :, n_prime - i + 1) * (H(:, i) - h_n_old), 0);
                        % max(.,0) ensures the nonnegativity
                    end
                else
                    for n_prime = i:n
                        V_hat(:, n_prime) = max(V_hat(:, n_prime) + W(:, :, n_prime - i + 1) * (H(:, i) - h_n_old), 0);
                    end
                end
            else
                H(:, i) = H(:, i) .* (num./denom) .^ gamma_beta(options.d_beta);
            end

        end

        if strcmp(options.alg, 'type2')

            V_hat = zeros(m, n);
            for i = 0 : t-1
                V_hat = V_hat + W(:, :, i+1) * shift_t(H, i);
            end
        end

        % update W
        for i = 0 : t - 1
            W_t_old = W(:, :, i+1);
            H_shift_t = shift_t(H, i);
            W(:, :, i+1) = W(:, :, i+1) .* (((((V + eps) .* (V_hat + eps) .^ (options.d_beta - 2)) * H_shift_t') ...
                ./ ((V_hat + eps).^(options.d_beta - 1) * H_shift_t'))) .^ gamma_beta(options.d_beta);
            V_hat = max(V_hat + (W(:, :, i+1) - W_t_old) * H_shift_t, 0);
            % max(.,0) ensures nonnegativity
        end

        [W, H] = renormalize_convNMF(W,H);

        if strcmp(options.alg, 'type1')
            % recalate the V_hat
            X_hat = zeros(size(V));
            for i=0:t-1
                tW = W(:, :, i+1);
                tH = shift_t(H, i);
                X_hat = X_hat + tW * tH;
            end
            V_hat = max(X_hat,0);
        end

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        [Wcon, Hcon] = reconstruct_wh(W, H, t);
        infos = store_nmf_info(V, Wcon, Hcon, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

function [W_concat, H_concat] = reconstruct_wh(W, H, t)

    W_concat = [];
    H_concat = [];
    for j = 1 : t
        W_concat = [W_concat W(:, :,j)];
        H_concat = [H_concat; shift_t(H, j-1)];
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function y = relu(x,dzdy)

  if nargin <= 1 || isempty(dzdy)
    y = max(x, single(0)) ;
  else
    y = dzdy .* (x > single(0)) ;
  end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function cheby_v_poly_values_test ( )

%*****************************************************************************80
%
%% CHEBY_V_POLY_VALUES_TEST demonstrates the use of CHEBY_V_POLY_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_V_POLY_VALUES_TEST:\n' );
  fprintf ( 1, '  CHEBYS_V_POLY_VALUES returns values of\n' );
  fprintf ( 1, '  the Chebyshev V polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N       X      V(N)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = cheby_v_poly_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %8f  %24.16f\n', n, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%run this test case with the command
%results = runtests('importExportTests.m')
function tests = importExportTests
tests = functiontests(localfunctions);
end

function testExcelImport(testCase)
sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
excelFile=fullfile(sourceDir,'tutorial','empty.xlsx');
%%Prepare test results, uncomment and run when new reference is needed
% modelExcel=importExcelModel(excelFile);
% sbmlFile=fullfile(sourceDir,'tutorial','empty.xml');
% modelSBML=importModel(sbmlFile);
% save(fullfile(sourceDir,'testing','unit_tests','test_data','importExportResults.mat'),'modelExcel','modelSBML');
evalc('model=importExcelModel(excelFile)'); % Repress warnings
load(fullfile(sourceDir,'testing','unit_tests','test_data','importExportResults.mat'), 'modelExcel');
verifyEqual(testCase,model,modelExcel)
end

function testSBMLImport(testCase)
sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
sbmlFile=fullfile(sourceDir,'tutorial','empty.xml');
evalc('model=importModel(sbmlFile)'); % Repress warnings
load(fullfile(sourceDir,'testing','unit_tests','test_data','importExportResults.mat'), 'modelSBML');
verifyEqual(testCase,model,modelSBML)
end

function testExcelExport(testCase)
sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model');
exportToExcelFormat(model,fullfile(sourceDir,'testing','unit_tests','test_data','_test.xlsx'));
%File will not be exactly equal as it contains the current date and time,
%so md5 or similar would not work. Just check whether file is reasonably
%sized.
s = dir(fullfile(sourceDir,'testing','unit_tests','test_data','_test.xlsx'));
filesize = s.bytes;
verifyTrue(testCase,filesize>17000);
delete(fullfile(sourceDir,'testing','unit_tests','test_data','_test.xlsx'));
end

function testSBMLExport(testCase)
sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model');
evalc('exportModel(model,fullfile(sourceDir,''testing'',''unit_tests'',''test_data'',''_test.xml''))');
%File will not be exactly equal as it contains the current date and time,
%so md5 or similar would not work. Just check whether file is reasonably
%sized.
s = dir(fullfile(sourceDir,'testing','unit_tests','test_data','_test.xml'));
filesize = s.bytes;
verifyTrue(testCase,filesize>18500);
delete(fullfile(sourceDir,'testing','unit_tests','test_data','_test.xml'));
end

author: SysBioChalmers

repo: RAVEN


text:
function dat = probability_maps_to_region_index(dat)
% Use dat.probability_maps to rebuild integer vector of index labels (dat.dat)
%
% dat = probability_maps_to_region_index(dat)
%

% Start: dat has one image per region, with probability values
% convert to integer vector

[maxval, condf] = max(double(full(dat.probability_maps)),[], 2);   % double is safer

allempty = all(dat.probability_maps == 0, 2) | isnan(maxval);  % some out-of-mask values may get NaNs

condf(allempty) = 0;

dat.dat = int32(condf);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function center = curveCentroid(varargin)
%CURVECENTROID compute centroid of a curve defined by a series of points
%
%   PT = curveCentroid(POINTS);
%   Computes center of mass of a curve defined by POINTS. POINTS is a [NxD]
%   array of double, N D-dimensional points.
%
%   PT = curveCentroid(PTX, PTY);
%   PT = curveCentroid(PTX, PTY, PTZ);
%   Specifies points as separate column vectors
%
%   PT = curveCentroid(..., TYPE);
%   Specifies if the last point is connected to the first one. TYPE can be
%   either 'closed' or 'open'.
%
%
%   See also :
%   polygons2d, centroid, polygonCentroid, curveLength
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 22/05/2006.
%

%   HISTORY
%   23/07/2009 deprecate and replace by 'reverseLine'.

% deprecation warning
warning('geom2d:deprecated', ...
    '''curveCentroid'' is deprecated, use ''polylineCentroid'' instead');

%% process input arguments

% check whether the curve is closed
closed = false;
var = varargin{end};
if ischar(var)
    if strcmpi(var, 'closed')
        closed = true;
    end
    varargin = varargin(1:end-1);
end

% extract point coordinates
if length(varargin)==1
    points = varargin{1};
elseif length(varargin)==2
    points = [varargin{1} varargin{2}];
end

% compute centers and lengths composing the curve
if closed
    centers = (points + points([2:end 1],:))/2;
    lengths = sqrt(sum(diff(points([1:end 1],:)).^2, 2));
else
    centers = (points(1:end-1,:) + points(2:end,:))/2;
    lengths = sqrt(sum(diff(points).^2, 2));
end

% centroid of edge centers weighted by edge length
%weigths = repmat(lengths/sum(lengths), [1 size(points, 2)]);
center = sum(centers.*repmat(lengths, [1 size(points, 2)]), 1)/sum(lengths);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function mtrClusterResampleSISPathways(fgFilenames,samplerOptsFilename,numSamplesPerCluster,clusterROIFilename,dt6Filename,fgResampledFilename, hintNumTotalFibers, bCalcStats)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% mtrClusterResampleSISPathways(fgFilenames,samplerOptsFilename,
%%                                numSamplesPerCluster,clusterROI,
%%                                fgResampledFilename)
%%
%% fgFilenames - filenames for metrotrac pathways
%% samplerOptsFilename - parameters file for collecting paths and stats
%% numSamplesPerCluster - total weight of all resamples
%% clusterROIFilename - ROI where we will resample pathways around each coordinate
%%              of the ROI.
%% fgResampledFilename - filename for resampled pathways
%% dt6Filename - Filename for dt6.
%% hintNumTotalFibers - Number hopefully greater than the total number of
%% fibers that we will load.  This will speed up processing significantly
%% as it will allow us to preallocate space for validity matrix.
%% bCalcStats - True then calculate the statistics False then they already
%% exist.
%%
%%
%% Author: Anthony Sherbondy
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Constants
DISTANCE_THRESH = 2;

%% MUST MAKE SURE THIS PROGRAM IS COMPILED AND IN RIGHT LOCATION
if(ispc)
    executable = which('updatestats.exe');
else
    error('Not compiled for Linux.');
end

% Get all the needed files
if(ieNotDefined('dt6Filename'))
    [f, p] = uigetfile({'*.mat';'*.*'}, 'Load DT6 file...');
    if(isnumeric(f)), disp('Resampling canceled.'); return; end
    dt6Filename = fullfile(p,f);
end

dt6 = load(dt6Filename);

% Load ROI
roi = load(clusterROIFilename);
roi = roi.roi;

% Get importance weight file that we need to resample with also load the
% pathfile into matlab
totalFiberCount = 1;
fg_input = dtiNewFiberGroup;
fg_temp = dtiNewFiberGroup;

% Allocate validity map with hint of number of total fibers
if(ieNotDefined('hintNumTotalFibers'))
    hinNumTotalFibers = 50000;
end
validFiberToROIMap = logical(zeros(size(roi.coords,1),hintNumTotalFibers));

for ii = 1:length(fgFilenames)
     tic;
    statvec_names{ii} = sprintf('temp%d_statvec_iw.dat',ii);
    if(bCalcStats)
        disp('Calculating Statistics.');
        args = sprintf(' -v %s -i %s -o temp.dat -s 4 -sf %s', samplerOptsFilename, fgFilenames{ii}, statvec_names{ii});
        cmd = [executable args];
        disp(cmd); disp('...')
        [s, ret_info] = system(cmd,'-echo');
        disp('Done.')
    end

    disp('Loading pathway file into matlab...')
    fg_temp = mtrImportFibers(fgFilenames{ii}, dt6.xformToAcPc);
    disp('Done')

    % For all of the ROI points see which pathways are close enough
    disp('Finding closest ROI points for each path...')

    for ff = 1:length(fg_temp.fibers)
        [indices, bestSqDist] = nearpoints(roi.coords', fg_temp.fibers{ff});
        % Coordinates are in mm
        validFiberToROIMap(:,totalFiberCount) = bestSqDist < DISTANCE_THRESH;
        totalFiberCount = totalFiberCount+1;
    end
    disp('Done')

    % Keep aggregating all of the fibers together into one group
    if(isempty(fg_input.fibers))
        fg_input.fibers = fg_temp.fibers;
    else
        fg_input.fibers = {fg_input.fibers{:}, fg_temp.fibers{:}};
    end
     toc;
end
clear fg_temp;

% Send the valid matrix in with the resample computations
[fg_output, sample_w] = resampleClusterSISPathways(numSamplesPerCluster, statvec_names, fg_input, roi, validFiberToROIMap);

% Output the pathways in dtifiberUI space
dtiWriteFiberGroup(fg_output,fgResampledFilename);

% Output the pathways in metrotrac space
%mtrExportFiberGroupToMetrotrac(fgResampledFilename,fg_output,dt6.t1NormParams,dt6.mmPerVox,dt6.xformToAcPc,size(dt6.b0));

% Cleanup temp files
delete('temp.dat');

% Why should I delete these stats?? they are yummy!!
% for ii = 1:length(statvec_names)
%     delete(statvec_names{ii});
% end

return;

function [fg_output, sample_w_out] = resampleClusterSISPathways(numSamplesPerCluster, weightFilenames, fg_input, roi, validFiberToROIMap)

numFiles = length(weightFilenames);
sizeFiles = [];

% Setup the output file
fg_output = dtiNewFiberGroup;

% Load multiple importance weight files
iw = [];
for f = 1:numFiles
    iwtemp = mtrLoadStatvec(weightFilenames{f});
    iw = [iw; iwtemp];
    sizeFiles(f) = length(iwtemp);
end
% Give everything a small chance of sampling
iw = exp(iw)+realmin;

% For each ROI point lets resample from the database
for rr = 1:size(roi.coords,1)
    iwROI = validFiberToROIMap(rr,:) .* iw(:)';
    % If no fibers are near enough to this position forget about it
    if(sum(iwROI)>0)
        % Resample pathways according to normalized weight and whether they are
        % near this ROI point
        sample_ind = randsample([1:length(iwROI)],numSamplesPerCluster,true,iwROI);
        unique_sample_ind = unique(sample_ind);
        unique_sample_w = 0;
        for i = 1:length(unique_sample_ind)
            unique_sample_w(i) =  sum(sample_ind == unique_sample_ind(i));
        end
        % Add fibers sampled from this ROI to the output fiber group
        if(isempty(fg_output.fibers))
            fg_output.fibers = {fg_input.fibers{unique_sample_ind}};
            sample_w_out = unique_sample_w;
        else
            fg_output.fibers = {fg_output.fibers{:}, fg_input.fibers{unique_sample_ind}};
            sample_w_out = [sample_w_out, unique_sample_w];
        end
    else
    end
end

return;


text:
function [q] = rotMatToQuat(R)

[r,c] = size( R );
R=R';

if( r ~= 3 | c ~= 3 )
    fprintf( 'R must be a 3x3 matrix\n\r' );
    return;
end

Rxx = R(1,1); Rxy = R(1,2); Rxz = R(1,3);
Ryx = R(2,1); Ryy = R(2,2); Ryz = R(2,3);
Rzx = R(3,1); Rzy = R(3,2); Rzz = R(3,3);

w = sqrt( trace( R ) + 1 ) / 2;

% check if w is real. Otherwise, zero it.
if( imag( w ) > 0 )
     w = 0;
end

x = sqrt( 1 + Rxx - Ryy - Rzz ) / 2;
y = sqrt( 1 + Ryy - Rxx - Rzz ) / 2;
z = sqrt( 1 + Rzz - Ryy - Rxx ) / 2;

[~, i ] = max( [w,x,y,z] );

if( i == 1 )
    x = ( Rzy - Ryz ) / (4*w);
    y = ( Rxz - Rzx ) / (4*w);
    z = ( Ryx - Rxy ) / (4*w);
end

if( i == 2 )
    w = ( Rzy - Ryz ) / (4*x);
    y = ( Rxy + Ryx ) / (4*x);
    z = ( Rzx + Rxz ) / (4*x);
end

if( i == 3 )
    w = ( Rxz - Rzx ) / (4*y);
    x = ( Rxy + Ryx ) / (4*y);
    z = ( Ryz + Rzy ) / (4*y);
end

if( i == 4 )
    w = ( Ryx - Rxy ) / (4*z);
    x = ( Rzx + Rxz ) / (4*z);
    y = ( Ryz + Rzy ) / (4*z);
end

q = [x; y; z; w];

end

repo: MSCKF


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
    %% EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
    %where you select the optimal (C, sigma) learning parameters to use for SVM
    %with RBF kernel
    %   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
    %   sigma. You should complete this function to return the optimal C and
    %   sigma based on a cross-validation set.
    %

    % ====================== YOUR CODE HERE ======================
    % Instructions: Fill in this function to return the optimal C and sigma
    %               learning parameters found using the cross validation set.
    %               You can use svmPredict to predict the labels on the cross
    %               validation set. For example,
    %                   predictions = svmPredict(model, Xval);
    %               will return the predictions on the cross validation set.
    %
    %  Note: You can compute the prediction error using
    %        mean(double(predictions ~= yval))
    %
    c_values = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30];
    sigma_values = c_values(:);
    best_accuracy = 0.0;

    for i = 1 : length(c_values)

        for j = 1 : length(sigma_values)

            % test this grid entry on our validation set
            test_C = c_values(i);
            test_sigma = sigma_values(j);
            model = svmTrain(X, y, test_C, @(x1, x2) gaussianKernel(x1, x2, test_sigma));
            predictions = svmPredict(model, Xval);
            accuracy = sum(predictions == yval) / length(yval);

            % is this the best performing so far?
            if accuracy > best_accuracy
                best_accuracy = accuracy;
                C = test_C;
                sigma = test_sigma;
            end
        end
    end

    fprintf('Best hyperparameters were C = %f, sigma = %f => best cv acc = %f\n', ...
        C, sigma, best_accuracy);

end

author: worldveil

repo: coursera-ml


text:
function value = func ( x )

%*****************************************************************************80
%
%% FUNC evaluates a function of X, as chosen by the user.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the function.
%
%    Output, real FUNC, the value of the function.
%
  ifunc = func_set ( 'GET', 'DUMMY' );

  if ( ifunc == 1 )
    value = 1.0;
  elseif ( ifunc == 2 )
    value = x;
  elseif ( ifunc == 3 )
    value = x.^2;
  elseif ( ifunc == 4 )
    value = x.^3;
  elseif ( ifunc == 5 )
    value = x.^4;
  elseif ( ifunc == 6 )
    value = x.^5;
  elseif ( ifunc == 7 )
    value = x.^6;
  elseif ( ifunc == 8 )
    value = x.^7;
  elseif ( ifunc == 9 )
    value = sin ( x );
  elseif ( ifunc == 10 )
    value = exp ( x );
  elseif ( ifunc == 11 )
    value = sqrt ( abs ( x ) );
  else
    value = 0.0;
  end

  return
end


text:
function [fx,tt]=fxrapt(s,fs,mode,q)
%FXRAPT RAPT pitch tracker [FX,VUV]=(S,FS,M,Q)
%
% Input:   s(ns)      Speech signal
%          fs         Sample frequency (Hz)
%          mode       'g' will plot a graph [default if no output arguments]
%                     'u' will include unvoiced fames (with fx=NaN)
%          q          stucture with parameter values (e.g. q.f0min=40); see below for a list
%
% Outputs: fx(nframe)     Larynx frequency for each fram,e (or NaN for silent/unvoiced)
%          tt(nframe,3)  Start and end samples of each frame. tt(*,3)=1 at the start of each talk spurt
%
% Plots a graph if no outputs are specified showing lag candidates and selected path
%

% Bugs/Suggestions:
%   (1) Include backward DP pass and output the true cost for each candidate.
%   (2) Add an extra state to distinguish between voiceless and silent
%   (3) N-best DP to allow longer term penalties (e.g. for frequent pitch doubling/halving)

% The algorithm is taken from [1] with the following differences:
%
%      (a)  the factor AFACT which in the Talkin algorithm corresponds roughly
%           to the absolute level of harmonic noise in the correlation window. This value
%           is here calculated as the maximum of three figures:
%                   (i) an absolute floor set by p.absnoise
%                  (ii) a multiple of the peak signal set by p.signoise
%                 (iii) a multiple of the noise floor set by p.relnoise
%      (b) The LPC used in calculating the Itakura distance uses a Hamming window rather than
%          a Hanning window.
%
% A C implementation of this algorithm by Derek Lin and David Talkin is included as  "get_f0.c"
% in the esps.zip package available from http://www.speech.kth.se/esps/esps.zip under the BSD
% license.
%
% Refs:
%      [1]   D. Talkin, "A Robust Algorithm for Pitch Tracking (RAPT)"
%            in "Speech Coding & Synthesis", W B Kleijn, K K Paliwal eds,
%            Elsevier ISBN 0444821694, 1995

%      Copyright (C) Mike Brookes 2006-2013
%      Version: $Id: fxrapt.m 9312 2017-01-19 13:19:13Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=s(:); % force s to be a column
if nargin<4
    q=[];
    if nargin<3
        mode=' ';
    end
end
doback=0;   % don't do backwards DP for now

% set default parameters

p0.f0min=50;           % Min F0 (Hz)
p0.f0max=500;          % Max F0 (Hz)
p0.tframe=0.01;        % frame size (s)
p0.tlpw=0.005;         % low pass filter window size (s)
p0.tcorw=0.0075;       % correlation window size (s)
p0.candtr=0.3;         % minimum peak in NCCF
p0.lagwt=0.3;          % linear lag taper factor
p0.freqwt=0.02;        % cost factor for F0 change
p0.vtranc=0.005;       % fixed voice-state transition cost
p0.vtrac=0.5;          % delta amplitude modulated transition cost
p0.vtrsc=0.5;          % delta spectrum modulated transition cost
p0.vobias=0.0;         % bias to encourage voiced hypotheses
p0.doublec=0.35;       % cost of exact doubling or halving
p0.absnoise=0;         % absolute rms noise level
p0.relnoise=2;         % rms noise level relative to noise floor
p0.signoise=0.001;     % ratio of peak signal rms to noise floor (0.001 = 60dB)
p0.ncands=20;          % max hypotheses at each frame
p0.trms=0.03;                      % window length for rms measurement
p0.dtrms=0.02;                     % window spacing for rms measurement
p0.preemph=-7000;                  % s-plane position of preemphasis zero
p0.nfullag=7;                      % number of full lags to try (must be odd)
p=paramsetch(p0,q);

% redefine commonly used parameters

candtr=p.candtr;          % minimum peak in NCCF                      [0.3]
vtranc=p.vtranc;          % fixed voice-state transition cost         [0.005]
vtrac=p.vtrac;            % delta amplitude modulated transition cost [0.5]
vtrsc=p.vtrsc;            % delta spectrum modulated transition cost  [0.5]
vobias=p.vobias;          % bias to encourage voiced hypotheses       [0.0]
doublec=p.doublec;        % cost of exact doubling or halving         [0.35]
ncands=p.ncands;          % max hypotheses at each frame              [20]
nfullag=p.nfullag;        % number of full lags to try (must be odd)  [7]

% derived parameters (mostly dependent on sample rate fs)

krms=round(p.trms*fs);            % window length for rms measurement
kdrms=round(p.dtrms*fs);          % window spacing for rms measurement
rmswin=hanning(krms).^2;
kdsmp=round(0.25*fs/p.f0max);
hlpw=round(p.tlpw*fs/2);          % force window to be an odd length
blp=sinc((-hlpw:hlpw)/kdsmp).*hamming(2*hlpw+1).';
fsd=fs/kdsmp;
kframed=round(fsd*p.tframe);      % downsampled frame length
kframe=kframed*kdsmp;           % frame increment at full rate
rmsix=(1:krms)+floor((kdrms-kframe)/2); % rms index according to Talkin; better=(1:krms)+floor((kdrms-krms+1)/2)
minlag=ceil(fsd/p.f0max);
maxlag=round(fsd/p.f0min);        % use round() only because that is what Talkin does
kcorwd=round(fsd*p.tcorw);        % downsampled correlation window
kcorw=kcorwd*kdsmp;             % full rate correlation window
spoff=max(hlpw-floor(kdsmp/2),1+kdrms-rmsix(1)-kframe);  % offset for first speech frame at full rate
sfoff=spoff-hlpw+floor(kdsmp/2); % offset for downsampling filter
sfi=1:kcorwd;                   % initial decimated correlation window index array
sfhi=1:kcorw;                   % initial correlation window index array
sfj=1:kcorwd+maxlag;
sfmi=repmat((minlag:maxlag)',1,kcorwd)+repmat(sfi,maxlag-minlag+1,1);
lagoff=(minlag-1)*kdsmp;        % lag offset when converting to high sample rate
beta=p.lagwt*p.f0min/fs;            % bias towards low lags
log2=log(2);
lpcord=2+round(fs/1000);        % lpc order for itakura distance
hnfullag=floor(nfullag/2);
jumprat=exp((doublec+log2)/2);  % lag ratio at which octave jump cost is lowest
ssq=s.^2;
csssq=cumsum(ssq);
sqrt(min(csssq(kcorw+1:end)-csssq(1:end-kcorw))/kcorw);
afact=max([p.absnoise^2,max(ssq)*p.signoise^2,min(csssq(kcorw+1:end)-csssq(1:end-kcorw))*(p.relnoise/kcorw)^2])^2*kcorw^2;

% downsample signal to approx 2 kHz to speed up autocorrelation calculation
% kdsmp is the downsample factor

sf=filter(blp/sum(blp),1,s(sfoff+1:end));
sp=filter([1 exp(p.preemph/fs)],1,s); % preemphasised speech for LPC calculation
sf(1:length(blp)-1)=[];         % remove startup transient
sf=sf(1:kdsmp:end);             % downsample to =~2kHz
nsf=length(sf);                 % length of downsampled speech
ns=length(s);                   % length of full rate speech

% Calculate the frame limit to ensure we don't run off the end of the speech or decimated speech:
%   (a) For decimated autocorrelation when calculating sff():  (nframe-1)*kframed+kcorwd+maxlag <= nsf
%   (b) For full rate autocorrelation when calculating sfh():  max(fho)+kcorw+maxlag*kdsamp+hnfllag <= ns
%   (c) For rms ratio window when calculating rr            :  max(fho)+rmsix(end) <= ns
% where max(fho) = (nframe-1)*kframe + spoff

nframe=floor(1+min((nsf-kcorwd-maxlag)/kframed,(ns-spoff-max(kcorw-maxlag*kdsmp-hnfullag,rmsix(end)))/kframe));

% now search for autocorrelation peaks in the downsampled signal

cost=zeros(nframe,ncands);      % cumulative cost
prev=zeros(nframe,ncands);      % traceback pointer
mcands=zeros(nframe,1);         % number of actual candidates excluding voiceless
lagval=repmat(NaN,nframe,ncands-1);    % lag of each voiced candidate
tv=zeros(nframe,3);             % diagnostics: 1=voiceless cost, 2=min voiced cost, 3:cumulative voiceless-min voiced
if doback
    costms=cell(nframe,1);
end

% Main processing loop for each 10 ms frame

for iframe=1:nframe       % loop for each frame (~10 ms)

    % Find peaks in the normalized autocorrelation of subsampled (2Khz) speech
    % only keep peaks that are > 30% of highest peak

    fho=(iframe-1)*kframe+spoff;
    sff=sf((iframe-1)*kframed+sfj);
    sffdc=mean(sff(sfi));       % mean of initial correlation window length
    sff=sff-sffdc;              % subtract off the mean
    nccfd=normxcor(sff(1:kcorwd),sff(minlag+1:end));
    [ipkd,vpkd]=v_findpeaks(nccfd,'q');

    % Debugging: execute the line below to plot the autocorrelation peaks.
    % v_findpeaks(nccfd,'q'); xlabel(sprintf('Lag = (x+%d)*%g ms',minlag-1,1000*kdsmp/fs)); ylabel('Normalized Cross Correlation'); title (sprintf('Frame %d/%d',iframe,nframe));

    vipkd=[vpkd ipkd];
    vipkd(vpkd<max(vpkd)*candtr,:)=[];          % eliminate peaks that are small
    if size(vipkd,1)
        if size(vipkd,1)>ncands-1
            vipkd=sortrows(vipkd);
            vipkd(1:size(vipkd,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end
        lagcan=round(vipkd(:,2)*kdsmp+lagoff);        % convert the lag candidate values to the full sample rate
        nlcan=length(lagcan);
    else
        nlcan=0;
    end

    % If there are any candidate lag values (nlcan>0) then refine their accuracy at the full sample rate

    if nlcan
        laglist=reshape(repmat(lagcan(:)',nfullag,1)+repmat((-hnfullag:hnfullag)',1,nlcan),nfullag*nlcan,1);
        sfh=s(fho+(1:kcorw+max(lagcan)+hnfullag));
        sfhdc=mean(sfh(sfhi));
        sfh=sfh-sfhdc;
        e0=sum(sfh(sfhi).^2);                     % energy of initial correlation window (only needed to store in tv(:,6)
        lagl2=repmat(lagcan(:)',nfullag+kcorw-1,1)+repmat((1-hnfullag:hnfullag+kcorw)',1,nlcan);
        nccf=normxcor(sfh(1:kcorw),sfh(lagl2),afact);

        [maxcc,maxcci]=max(nccf,[],1);
        vipk=[maxcc(:) lagcan(:)+maxcci(:)-hnfullag-1];
        vipk=vipk(:,[1 2 2]);
        maxccj=maxcci(:)'+nfullag*(0:nlcan-1);    % vector index into nccf array
        msk=mod(maxcci,nfullag-1)~=1 & 2*nccf(maxccj)-nccf(mod(maxccj-2,nfullag*nlcan)+1)-nccf(mod(maxccj,nfullag*nlcan)+1)>0;  % don't do quadratic interpolation for the end ones
        if any(msk)
            maxccj=maxccj(msk);
            vipk(msk,3)=vipk(msk,3)+(nccf(maxccj+1)-nccf(maxccj-1))'./(2*(2*nccf(maxccj)-nccf(maxccj-1)-nccf(maxccj+1)))';
        end
        vipk(maxcc<max(maxcc)*candtr,:)=[];          % eliminate peaks that are small
        if size(vipk,1)>ncands-1
            vipk=sortrows(vipk);
            vipk(1:size(vipk,1)-ncands+1,:)=[];   % eliminate lowest to leave only ncands-1
        end

        % vipk(:,1) has NCCF value, vipk(:,2) has integer peak position, vipk(:,3) has refined peak position

        mc=size(vipk,1);
    else
        mc=0;
    end

    % We now have mc lag candidates at the full sample rate

    mc1=mc+1;               % total number of candidates including "unvoiced" possibility
    mcands(iframe)=mc;      % save number of lag candidates (needed for pitch consistency cost calculation)
    if mc
        lagval(iframe,1:mc)=vipk(:,3)';
        cost(iframe,1)=vobias+max(vipk(:,1));   % voiceless cost
        cost(iframe,2:mc1)=1-vipk(:,1)'.*(1-beta*vipk(:,3)');   % local voiced costs
        tv(iframe,2)=min(cost(iframe,2:mc1));
    else
        cost(iframe,1)=vobias;          % if no lag candidates (mc=0), then the voiceless case is the only possibility
    end
    tv(iframe,1)=cost(iframe,1);
    if iframe>1                         % if it is not the first frame, then calculate pitch consistency and v/uv transition costs
        mcp=mcands(iframe-1);
        costm=zeros(mcp+1,mc1);         % cost matrix: rows and cols correspond to candidates in previous and current frames (incl voiceless)

        % if both frames have at least one lag candidate, then calculate a pitch consistency cost

        if mc*mcp
            lrat=abs(log(repmat(lagval(iframe,1:mc),mcp,1)./repmat(lagval(iframe-1,1:mcp)',1,mc)));
            costm(2:end,2:end)=p.freqwt*min(lrat,doublec+abs(lrat-log2));  % allow pitch doubling/halving
        end

        % if either frame has a lag candidate, then calculate the cost of voiced/voiceless transition and vice versa

        if mc+mcp
            rr=sqrt((rmswin'*s(fho+rmsix).^2)/(rmswin'*s(fho+rmsix-kdrms).^2)); % amplitude "gradient"
            ss=0.2/(distitar(lpcauto(sp(fho+rmsix),lpcord),lpcauto(sp(fho+rmsix-kdrms),lpcord),'e')-0.8);   % Spectral stationarity: note: Talkin uses Hanning instead of Hamming windows for LPC
            costm(1,2:end)= vtranc+vtrsc*ss+vtrac/rr;   % voiceless -> voiced cost
            costm(2:end,1)= vtranc+vtrsc*ss+vtrac*rr;
            tv(iframe,4:5)=[costm(1,mc1) costm(mcp+1,1)];
        end
        costm=costm+repmat(cost(iframe-1,1:mcp+1)',1,mc1);  % add in cumulative costs
        [costi,previ]=min(costm,[],1);
        cost(iframe,1:mc1)=cost(iframe,1:mc1)+costi;
        prev(iframe,1:mc1)=previ;
    else                            % first ever frame
        costm=zeros(1,mc1); % create a cost matrix in case doing a backward recursion
    end
    if mc
        tv(iframe,3)=cost(iframe,1)-min(cost(iframe,2:mc1));
        tv(iframe,6)=5*log10(e0*e0/afact);
    end
    if doback
        costms{iframe}=costm; % need to add repmatted cost into this
    end
end

% now do traceback

best=zeros(nframe,1);
[cbest,best(nframe)]=min(cost(nframe,1:mcands(nframe)+1));
for i=nframe:-1:2
    best(i-1)=prev(i,best(i));
end
vix=find(best>1);
fx=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
fx(vix)=fs*lagval(vix+nframe*(best(vix)-2)).^(-1); % leave as NaN if unvoiced
tt=zeros(nframe,3);
tt(:,1)=(1:nframe)'*kframe+spoff;       % find frame times
tt(:,2)=tt(:,1)+kframe-1;
jratm=(jumprat+1/jumprat)/2;
tt(2:end,3)=abs(fx(2:end)./fx(1:end-1)-jratm)>jumprat-jratm;    % new spurt if frequency ratio is outside (1/jumprat,jumprat)
tt(1,3)=1;           % first frame always starts a spurt
tt(1+find(isnan(fx(1:end-1))),3)=1; % NaN always forces a new spurt

% plot results if there are no output arguments of if the 'g' mode option is specified

if ~nargout || any(mode=='g')
    tf=spoff+(0:nframe-1)'*kframe;      % one sample before start of each frame
    blag=repmat(NaN,nframe,1);                        % unvoiced frames will be NaN
    blag(vix)=lagval(vix+nframe*(best(vix)-2)); % leave as NaN if unvoiced
    ts=(1:ns)/fs;                       % time scale for speech samples
    tsa=[1:tf(1) tf(end)+kframe+1:ns];  % indexes for unprocessed speech [-1 term is an error methinks]
    sup=repmat(NaN,ns,1);               % unprocessed speech - plot in black
    sup(tsa)=s(tsa);
    sv=reshape(s(tf(1)+1:tf(end)+kframe),kframe,nframe);               % processed speech
    su=sv;
    su(:,best>1)=NaN;                   % delete all voiced samples
    sv(:,best==1)=NaN;                  % delete all unvoiced samples
    tsuv=(tf(1)+1:tf(end)+kframe)/fs;
    su=su(:);
    sv=sv(:);
    ax=zeros(2,1);
    ax(1)=subplot(211);
    plot(ts,sup,'-k',tsuv,su,'r-',tsuv,sv,'b-');
    title('Speech');
    ax(2)=subplot(212);
    plot((tf+(kframe+1)/2)/fs,lagval*1000/fs,'xr',(tf+(kframe+1)/2)/fs,blag*1000/fs,'-b')
    xlabel('Time (s)');
    ylabel('Period (ms)');
    title('Lag Candidates');
    linkaxes(ax,'x');
end
if ~any(mode=='u')
    tt(isnan(fx),:)=[];    % remove NaN spurts
    fx(isnan(fx),:)=[];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function v=normxcor(x,y,d)
% Calculate the normalized cross correlation of column vectors x and y
% we can calculate this in two ways but fft is much faster even for nx small
% We must have nx<=ny and the output length is ny-nx+1
% note that this routine does not do mean subtraction even though this is normally a good idea
% if y is a matrix, we correlate with each column
% d is a constant added onto the normalization factor
% v(j)=x'*yj/sqrt(d + x'*x * yj'*yj) where yj=y(j:j+nx-1) for j=1:ny-nx+1

if nargin<3
    d=0;
end
nx=length(x);
[ny,my]=size(y);
nv=1+ny-nx;
if nx>ny
    error('second argument is shorter than the first');
end

nf=pow2(nextpow2(ny));
w=irfft(repmat(conj(rfft(x,nf,1)),1,my).*rfft(y,nf,1));
s=zeros(ny+1,my);
s(2:end,:)=cumsum(y.^2,1);
v=w(1:nv,:)./sqrt(d+(x'*x).*(s(nx+1:end,:)-s(1:end-nx,:)));


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Copyright 2015-2018 Jonas Koennemanm, Giovanni Licitra
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function p_vec = pendulum_sim

  conf = struct;
  conf.l = 1;
  conf.m = 1;

  simulator = ocl.Simulator(@ocl.examples.pendulum.varsfun, ...
    @ocl.examples.pendulum.daefun, ...
    @ocl.examples.pendulum.icfun, ...
    'userdata', conf);

  x0 = simulator.getStates();
  x0.p.set([0;conf.l]);
  x0.v.set([-0.5;-1]);

  times = 0:0.1:4;

  figure
  simulator.reset(x0);

  p_vec = zeros(2,length(times));
  p_vec(:,1) = x0.p.value;

  for k=1:length(times)-1

    dt = times(k+1)-times(k);
    [x,~] = simulator.step(10, dt);

    p_vec(:,k+1) = x.p.value;
  end

  ocl.examples.pendulum.animate(conf.l, p_vec, times);
  snapnow;
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function bool= str_matchesHead(head, str)
%STR_MATCHESHEAD - Is one string the head of another string?
%
%Synopsis:
% TF= str_matchesHead(HEAD, STR)
%
%Arguments:
%  HEAD:  CHAR - Short string
%  STR :  CHAR - String that is presumably shorter
%
%Description:
%  This function decides whether the string HEAD is the head, i.e.,
%  the first part of the second string STR.
%
%Returns:
%  TF : BOOL - Truth value of HEAD being the head of STR

% 06-2012 Some Genius

bool= strncmp(head, str, length(head));

author: PatternRecognition

repo: OpenBMI


text:
function [inplane_pix_size] = mrGetIPPixSize()
%function [inplane_pix_size] = mrGetIPPixSize()
%
% PURPOSE: Prompt the user for the inplane anatomy pixel size
% AUTHOR:  Poirson
% DATE:    07.16.97
% HISTORY: Based on routine by Geoff Boynton
% NOTES:   The default values [256/260,256/260,1/4] mean that
%          26 cm of brain image are interpolated on the 256 pixels
%          in the x and y dimensions (units are pixels/mm).
%          And that 1 pixel in the inplane direction equals 4mm thickness.
%

% 04/13/00 huk and nestares -- gets defPixSize from mrSESSION now
% 07/27/00 wandell, brewer  -- added extra check for reconParams field

global mrSESSION

qt=''''; %single quote character

if    isfield(mrSESSION,'fullInplaneSize') & ...
      isfield(mrSESSION,'reconParams') & ...
      isfield(mrSESSION.reconParams(1),'FOV') & ...
      isfield(mrSESSION.reconParams(1),'sliceThickness')
  defPixSize =  [mrSESSION.fullInplaneSize(1)/ mrSESSION.reconParams(1).FOV, ...
                 mrSESSION.fullInplaneSize(2)/mrSESSION.reconParams(1).FOV, ...
                 1/mrSESSION.reconParams(1).sliceThickness];
  disp('Calculating inplane anatomy voxel size from fullInplaneSize, FOV, and sliceThickness.');
else
  defPixSize = [256/260,256/260,1/4];
  disp('mrSESSION inplane voxel size info does not exist.  Default values set; check your protocol parameters.');
end

disp('Enter size of inplane anatomy pixels/mm in x,y and z directions');
%inplane_pix_size=input('Default is [256/260,256/260,1/4]: ');
defStr = mat2str(defPixSize);
inplane_pix_size=input(['Default is: ' defStr ': ']);
%Or go the default
if isempty(inplane_pix_size)
	inplane_pix_size= defPixSize;
end

return


text:
function out = proj_l1ball_box(x,w,r,u)
%PROJ_L1BALL_BOX computes the orthogonal projection onto the intersection of an l1 ball and a box
%                                                  {x: norm(w(:).*x(:),1)<=r, -u<=x<=u}
%
%  Usage:
%  out = PROJ_L1BALL_BOX(x,[w],[r],[u])
%  =============================================================
%  INPUT:
%  x - point to be projected (vector/matrix)
%  w - weights vecotr - [default: 1]
%  r - positive scalar - [default: 1]
%  u - box paramters (scalar/vector/matrix)- [default: Inf]
%  ==============================================================
%  Assumptions:
%  r,u >= 0
%  w >= 0
%  ==============================================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%reading the user x and setting defalut values when required.
if (nargin < 1)
    error ('usage: proj_l1_inter_box(x,[w],[r],[u])') ;
end

if (nargin < 4)
    %box parameters are not given, setting to defalut value :inf
    u = inf*(ones(size(x))) ;
end

if ((nargin < 3) || (isempty(r)))
    %ball radius is not given, setting to defalut value :1
    r = 1;
end

if ((nargin < 2) || (isempty(w)))
    %weight vector is not given, setting to default value: ones
    w=ones(size(x)) ;
end

if  ((r < 0) || (min(min(u)) < 0))
      error('Set is infeasible') ;
end

if (min(min(w)) < 0)
    error('usage: proj_l1_inter_box(x,[w],[r],[u]) - w should be a non-negative vector or matrix');
end

%checking the simple projection
out = min(max(x,-u),u) ;
if ((sum(sum(w .* abs (out)))) <= r)
    return;
end

%defining f on lambda - to be used by the bisetion

f= @(lam)   trace(w' * abs(min(max(abs(x)-lam*w,0),u) .* sign(x)))- r;
lambda_min = 0;

lambda_max = 1;
while(f(lambda_max)>0)
    lambda_max = lambda_max *2  ;
end

eps = 1e-10 ;
final_lam = bisection(f,lambda_min,lambda_max,eps) ;
out=  min(max(abs(x)-final_lam*w,0),u) .* sign(x) ;

repo: SGDLibrary


text:
function [L1t,L2t,L3t]=triplet_find(DT1,DT2,DT3,Vn1,Vn2,Vn3,T_dist,opt)

switch opt
    case 1 %Full Delaunay input based
        [~,Dc12] = nearestNeighbor(DT2,Vn1);
        [~,Dc13] = nearestNeighbor(DT3,Vn1);
        L1t=Dc12<=T_dist & Dc13<=T_dist;

        [~,Dc21] = nearestNeighbor(DT1,Vn2);
        [~,Dc23] = nearestNeighbor(DT3,Vn2);
        L2t=Dc21<=T_dist & Dc23<=T_dist;

        [~,Dc31] = nearestNeighbor(DT1,Vn3);
        [~,Dc32] = nearestNeighbor(DT2,Vn3);
        L3t=Dc31<=T_dist & Dc32<=T_dist;

    case 2 %Accelerating renew with cropped Delaunay if possible otherwise full

        [~,Dc12] = nearestNeighbor(DT2,Vn1);
        L12t=Dc12<=T_dist;
        if any(L12t)
            [~,Dc13] = nearestNeighbor(DT3,Vn1);
            L1t=L12t & Dc13<=T_dist;
        else
            L1t=false(size(Vn1,1),1);
        end

        if any(L1t)
            if nnz(L1t)>=4
                DT=delaunayTriangulation(Vn1(L1t,:)); %Tesselation of selected coordinates
                if ~isempty(DT.ConnectivityList)
                    DT1=DT;
                end
            end

            [~,Dc21] = nearestNeighbor(DT1,Vn2);
            L21t=Dc21<=T_dist;
            if any(L21t)
                [~,Dc23] = nearestNeighbor(DT3,Vn2);
                L2t=L21t & Dc23<=T_dist;
            else
                L2t=false(size(Vn2,1),1);
            end
        else
            L2t=false(size(Vn2,1),1);
        end

        if any(L2t)
            if nnz(L2t)>=4
                DT=delaunayTriangulation(Vn2(L2t,:)); %Tesselation of selected coordinates
                if ~isempty(DT.ConnectivityList)
                    DT2=DT;
                end
            end

            [~,Dc31] = nearestNeighbor(DT1,Vn3);
            L31t=Dc31<=T_dist;
            if any(L31t)
                [~,Dc32] = nearestNeighbor(DT2,Vn3);
                L3t=L31t & Dc32<=T_dist;
            else
                L3t=false(size(Vn3,1),1);
            end
        else
            L3t=false(size(Vn3,1),1);
        end

    case 3 %Accelerating renew with cropped Delaunay/non-Delaunay
        [~,Dc12] = nearestNeighbor(DT2,Vn1);
        [~,Dc13] = nearestNeighbor(DT3,Vn1);
        L1t=Dc12<=T_dist & Dc13<=T_dist;
        if any(L1t)
            if nnz(L1t)>=4
                DT=delaunayTriangulation(Vn1(L1t,:)); %Tesselation of selected coordinates
                if ~isempty(DT.ConnectivityList)
                    DT1=DT;
                end
                [~,Dc21] = nearestNeighbor(DT1,Vn2);
                [~,Dc31] = nearestNeighbor(DT1,Vn3);
            else
                [INDc,Dc12] = dsearchn(Vn2,Vn1(L1t,:));
                Dc21=NaN(size(Vn2,1),1); Dc21(INDc)=Dc12;

                [INDc,Dc13] = dsearchn(Vn3,Vn1(L1t,:));
                Dc31=NaN(size(Vn3,1),1); Dc31(INDc)=Dc13;
            end
            [~,Dc23] = nearestNeighbor(DT3,Vn2);
            L2t=Dc21<=T_dist & Dc23<=T_dist;
        else
            L2t=false(size(Vn2,1),1);
        end
        if any(L2t)
            if nnz(L2t)>=4 && nnz(L1t)>=4
                 DT=delaunayTriangulation(Vn2(L2t,:)); %Tesselation of selected coordinates
                if ~isempty(DT.ConnectivityList)
                    DT2=DT;
                end
                [~,Dc32] = nearestNeighbor(DT2,Vn3);
            else
                [INDc,Dc23] = dsearchn(Vn3,Vn2(L2t,:));
                Dc32=NaN(size(Vn3,1),1); Dc32(INDc)=Dc23;
            end
            L3t=Dc31<=T_dist & Dc32<=T_dist;
        else
            L3t=false(size(Vn3,1),1);
        end
end
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% FUNCTION Least_SparseTrace
%   Incoherent Sparse and Low-Rank Learning with Least Squares Loss.
%
%% OBJECTIVE
%   argmin_W { sum_i^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%            + gamma * \|P\|_1 }
%   subject to: W = P + Q, \|Q\|_* <= tau
%   where \|Q\|_*      = sum( svd(Q, 0) )
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   gamma: sparse component P L1,2-norm sprasity controlling parameter
%   tau: low rank component Q trace-norm regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   funcVal: function value vector.
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Jianhui and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Chen, J., Liu, J. and Ye, J. Learning Incoherent Sparse and Low-Rank
%   Patterns from Multiple Tasks, KDD 2010
%
%% RELATED FUNCTIONS
%   init_opts

function [W funcVal Tp_val Tq_val] = Least_SparseTrace(X, Y, gamma, tau, opts )

if nargin <4
    error('\n Inputs: X, Y, rho1, and rho2 should be specified!\n');
end
X = multi_transpose(X);

if nargin <5
    opts = [];
end

% initialize options.
opts=init_opts(opts);

zero_threshold = 0;

d = size(X{1}, 1); m = length(Y);

% Initialize L and S
if opts.init==2
    P0 = zeros(d, m);
    Q0 = zeros(d, m);
elseif opts.init == 0
    P0 = randn(d, m);
    Q0 = P0;
else
    if isfield(opts,'P0')
        P0=opts.P0;
        if (nnz(size(P0)-[d, m]))
            error('\n Check the input .L0');
        end
    else
        P0=randn(d, m);
    end

    if isfield(opts,'Q0')
        Q0=opts.Q0;
        if (nnz(size(Q0)-[d, m]))
            error('\n Check the input .S0');
        end
    else
        Q0=P0;
    end

end

% Initialize
Tp_i = P0;
Tq_i = Q0;
Tp_i_min_1 = P0;
Tq_i_min_1 = Q0;

% initial parameters setting
%  ++++++++++++++++++++++++++++++++++++++++++

% Initialize L
Li = 1;

% Set an array to save the objective value
funcVal = [ ];

% Set auxiliary parameters

t_i_min_2 = 0;
t_i_min_1 = 1;
XY = cell(m, 1);
for jj = 1:m
    XY{jj} = X{jj} * Y{jj};
end

% main loop
% ++++++++++++++++++++++++++++++++++++++++++
for iter = 1 : opts.maxIter

    alpha_i = ( t_i_min_2 - 1 ) / t_i_min_1;

    Sp = ( 1 + alpha_i ) * Tp_i - alpha_i * Tp_i_min_1;
    Sq = ( 1 + alpha_i ) * Tq_i - alpha_i * Tq_i_min_1;

    derivative = zeros(d, m);
    if opts.pFlag
        parfor kk = 1:m
            Xi = X{kk};
            col = Sp(:, kk) + Sq(:, kk);
            derivative(:, kk) = 2 * Xi * (Xi' * col)  - 2 * XY{kk};
        end
    else
        for kk = 1:m
            Xi = X{kk};
            col = Sp(:, kk) + Sq(:, kk);
            derivative(:, kk) = 2 * Xi * (Xi' * col)  - 2 * XY{kk};
        end
    end

    lambda_old = 1;
    sign = 1;

    while (sign)

        beta = Li / 2;

        hat_Sp = Sp - derivative / Li;
        hat_Sq = Sq - derivative / Li;

        % Compute Tp
        Tp = Solve_OneNorm(hat_Sp, beta, gamma);

        % Compute Tq
        [U S V] = svd(hat_Sq, 'econ');
        s_bar = diag(S);
        s_bar = s_bar( s_bar > zero_threshold );
        U = U( :, 1:length(s_bar));
        V = V( :, 1:length(s_bar));

        [s_hat, lambda_new ]=eplb(s_bar, size(s_bar, 1), tau, lambda_old);
        Tq = U * diag(s_hat) * V';

        [sign qp_term]= line_search_cond_sparse_lowrank(Tp, Tq, Sp, Sq, X, Y, derivative, Li, opts);

        if (sign)
            Li = Li * 2;
            lambda_old = lambda_new;
        end

    end

    funcVal = cat(1, funcVal,  qp_term +  gamma * sum( abs(Tp(:)) ) );

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    t_i_min_2 = t_i_min_1;
    t_i_min_1 = ( 1 + ( 1+4*t_i_min_1^2 )^0.5 ) / 2;

    Tp_i_min_1 = Tp_i;
    Tq_i_min_1 = Tq_i;

    Tp_i = Tp;
    Tq_i = Tq;

end % for loop

Tp_val = Tp;
Tq_val = Tq;
W = Tp_val + Tq_val;

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
%Works reasonably well since double represents integers <= about 2^50 exactly
function I = colon(a,b,c)
	if nargin == 2
		I  = double(a):double(b);
	else
		I = double(a):double(b):double(c);
	end
	I = int64(I);
end

author: Sable

repo: mcbench-benchmarks


text:
classdef AcceleratedForwardBackward < SplittingAlgorithm

    properties (Access = private)
        gradientMethod
        proximal
        momentum
    end

    methods (Access = public)

        function obj = AcceleratedForwardBackward(cParams)
            obj.proximal = cParams.proximal;
            obj.createGradientMethods(cParams);
            obj.createMomentum(cParams);
        end

        function update(obj)
            obj.momentum.apply();
            obj.gradientMethod.compute();
            obj.proximal.solve();
        end

    end

    methods (Access = private)

        function createGradientMethods(obj,cParams)
            s = cParams.gradientMethodParams;
            obj.gradientMethod = GradientMethod(s);
        end

        function createMomentum(obj,cParams)
            s = cParams.momentumParams;
            obj.momentum = Momentum(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = subsref(a,L)
% function B=subsref(A,L);
%
% DESCRIPTION
%   Subsreference for polynomial objects.
%
% INPUTS
%   A: polynomial
%   L: a structure array with the fields:
%    type -- string containing '()'  or '.' specifying the
%              subscript type.
%    subs -- Cell array or string containing the actual subscripts.
%
% OUTPUTS
%   B: object after referencing
%
% SYNTAX
%   B=subsref(A,L)

% 6/7/2002: PJS  Initial Coding

switch L(1).type
    case '.'
        b = get(a,L(1).subs);
    case '()'
        sza = a.matdim;
        nra = sza(1);
        nca = sza(2);
        acoef = a.coefficient;
        b = a;

        if length(L(1).subs)==1
            % Check Indices
            subsidx = L(1).subs{1};

            if strcmp(subsidx,':')
                subsidx = (1:nra*nca)';
            end
            if islogical(subsidx)
                subsidx = find(subsidx);
            end

            if ~isempty(subsidx) && ...
                    (max(subsidx)>nra*nca || min(subsidx)<1)
                error('Index exceeds matrix dimensions.');
            end

            % Do subsref
            b.coefficient = acoef(:,subsidx);
            if isempty(subsidx)
                b = polynomial(zeros(size(subsidx)));
            elseif nca==1
                % Column vector should remain column vector
                b.matdim = [length(subsidx) 1];
            else
                % All other single index references convert to row vector
                b.matdim = size(subsidx);
            end

        elseif length(L(1).subs)==2
            ridx = L(1).subs{1};
            if strcmp(ridx,':')
                ridx = 1:nra;
            end
            cidx = L(1).subs{2};
            if strcmp(cidx,':');
                cidx = 1:nca;
            end

            % Check Indices
            if ~isempty(ridx) && (max(ridx)>nra || min(ridx)<1)
                error('Index exceeds matrix dimensions.');
            elseif ~isempty(cidx) && (max(cidx)>nca || min(cidx)<1)
                error('Index exceeds matrix dimensions.');
            end

            % Do subsref
            bcoef = acoef;
            idx = [];
            for i1 = 1:length(cidx)
                idx = [idx ridx+(cidx(i1)-1)*nra];
            end
            b.coefficient = bcoef(:,idx);
            b.matdim = [length(ridx) length(cidx)];

        else
            error('Invalid subsref for polynomials');
        end

%         % Remove any terms with zero coeffs and then combine down
%         % (This is faster than directly calling combine without
%         %  pruning out terms with zero coefs)
%         bcoef = b.coefficient;
%         [ridx,cidx] = find( bcoef );
%         ridx = unique(ridx);
%         b.coefficient = bcoef(ridx,:);
%         b.degmat = b.degmat(ridx,:);

        if isempty(b.coefficient)
            b = polynomial(zeros(b.matdim));
        else
            b = combine(b);
        end

    case '{}'
        error('{}- like subsreference is not supported for polynomial objects.');
end

if length(L)>1
    b = subsref(b,L(2:end));
end
varargout{1}=b;

author: yu-jiang

repo: radpbook

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


text:
function save_as_pdf(plot_obj,outfilename)
% Saves the plot to a PDF file.
% Inputs:
%   outfilename: The name of the PDF file in which the plot will be
%     saved.

assert(ischar(outfilename))
saveas(plot_obj.fh,outfilename,'pdf');
end

author: nesl


text:
function writeKaldiNetworks_FeedForward(layer, kaldiFileName)

fprintf('Save DNN into Kaldi network format: %s!\n', kaldiFileName);

FID = fopen(kaldiFileName, 'w');
fprintf(FID, '<Nnet>\n');

for i=1:length(layer)
    fprintf('Write layer %d - %s - %s\n', i, layer{i}.name, datestr(now));
    switch lower(layer{i}.name)
        case 'affine'
            [outputSize, inputSize] = size(layer{i}.W);
            fprintf(FID, '<AffineTransform> %d %d\n', outputSize, inputSize);
            fprintf(FID, '<LearnRateCoef> 1 <BiasLearnRateCoef> 1 <MaxNorm> 0  [\n');

            for j=1:outputSize-1
                fprintf(FID, '%f ', layer{i}.W(j,:));
                fprintf(FID,'\n');
            end
            fprintf(FID, '%f ', layer{i}.W(end,:));
            fprintf(FID,']\n');

            fprintf(FID, '[ ');
            fprintf(FID, '%f ', layer{i}.b);
            fprintf(FID, ']\n');
        case 'sigmoid'
            fprintf(FID, '<Sigmoid> %d %d\n', outputSize, outputSize);
        case 'softmax'
            fprintf(FID, '<Softmax> %d %d\n', outputSize, outputSize);
    end
end

fprintf(FID, '</Nnet>\n');
fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function [sys,x0,str,ts] = ArticulatedVehicleSFunction(t,x,u,flag)
% This file is a s-function template for simulating the articulated vehicle
% model in Simulink.

% Choosing tire model
TireModel = VehicleDynamicsLateral.TirePacejka();
% Defining tire parameters
TireModel.a0        = 1;
TireModel.a1        = 2;
TireModel.a2        = 700;
TireModel.a3        = 5000;
TireModel.a4        = 80;
TireModel.a5        = 0;
TireModel.a6        = 0;
TireModel.a7        = 0.6;
TireModel.a8        = 0;
TireModel.a9        = 0;
TireModel.a10       = 0;
TireModel.a11       = 0;
TireModel.a12       = 0;
TireModel.a13       = 0;

% Choosing vehicle model
VehicleModel = VehicleDynamicsLateral.VehicleArticulatedNonlinear();
% Defining vehicle parameters
VehicleModel.mF0    = 5200;
VehicleModel.mR0    = 2400;
VehicleModel.mF     = 6000;
VehicleModel.mR     = 10000;
VehicleModel.mM     = 17000;
VehicleModel.IT     = 46000;
VehicleModel.IS     = 450000;
VehicleModel.lT     = 3.5;
VehicleModel.lS     = 7.7;
VehicleModel.c      = -0.3;
VehicleModel.nF     = 2;
VehicleModel.nR     = 4;
VehicleModel.nM     = 8;
VehicleModel.wT     = 2.6;
VehicleModel.wS     = 2.4;
VehicleModel.muy    = 0.8;
VehicleModel.deltaf = 0;
VehicleModel.Fxf    = 0;
VehicleModel.Fxr    = 0;
VehicleModel.Fxm    = 0;
VehicleModel.tire   = TireModel;

switch flag

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0
    [sys,x0,str,ts]=mdlInitializeSizes();

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u,VehicleModel);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3
    sys=mdlOutputs(t,x,u,VehicleModel);

  %%%%%%%%%%%%%%%%%%%
  % Unhandled flags %
  %%%%%%%%%%%%%%%%%%%
  case { 2, 4, 9 }
    sys = [];

  %%%%%%%%%%%%%%%%%
  % Vehicle model %
  %%%%%%%%%%%%%%%%%
  % Case 5 returns the vehicle model.
  case 5
    sys = VehicleModel;
    x0  =1; % Dummy
    str =1; % Dummy
    ts  =1; % Dummy

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end
% end csfunc

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%

function [sys,x0,str,ts]=mdlInitializeSizes()

% Definitions
sizes = simsizes;
sizes.NumContStates  = 8;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 8;
sizes.NumInputs      = 6;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

% Setting initial conditions
VEL0 = 20; % Initial velocity [m/s]

x0  = [0 0 0 0 VEL0 0 0 0];
str = [];
ts  = [0 0];

% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%

function sys = mdlDerivatives(t,x,u,vehicle)

% Defining input
vehicle.deltaf      = u(1);
vehicle.deltar      = u(2);
vehicle.deltam      = u(3);
vehicle.Fxf         = u(4);
vehicle.Fxr         = u(5);
vehicle.Fxm         = u(6);

% Getting the vehicle model function (state equations)
ModelFunction       = @vehicle.Model;
% Getting the mass matrix function of the model
MassMatrixFunction  = @vehicle.MassMatrix;

% The derivatives are calculated with the inverse of the MassMatrix
sys = MassMatrixFunction(t,x)\ModelFunction(t,x,0);

% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%

function sys = mdlOutputs(~,x,~,~)

% Output all state variables
sys = x;

% end mdlOutputs

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function solution = CS_GWO(UAV, SearchAgents, Max_iter)
%CS_GWO 灰狼-布谷鸟优化算法
%Gray Wolf Cuckoo Optimization

% 超参数
pa = 0.25;   % 布谷鸟搜索参数

% 算法初始化
[WolfPops, Tracks] = PopsInit(UAV, SearchAgents, false);   % 随机生成 初始狼群 和 轨迹们
dim = WolfPops.PosDim;                                                         % 状态变量维度

% 初始化解
Alpha_pos = zeros(1, dim);   % α解
Alpha_score = inf;                  % α解适应度
Alpha_no = 1;                         % α解编号

Beta_pos = zeros(1, dim);      % β解
Beta_score = inf;                     % β解适应度
Beta_no = 1;                             % β解编号

Delta_pos = zeros(1, dim);     % δ解
Delta_score = inf;                    % δ解适应度
Delta_no = 1;                           % δ解编号

Fitness_list = zeros(1, Max_iter);

% 迭代求解
tic
fprintf('>>CS_GWO 优化中    00.00%%')
for iter = 1 : Max_iter

    % ①  计算每只狼的适应度，更改其种群等级
    ProbPoints = cell(SearchAgents, 1);
    for i = 1 : SearchAgents
        % 计算目标函数
        [fitness, ~, Data] = ObjFun(Tracks{i}, UAV);    % 一个智能体的目标函数
        ProbPoints{i} = Data.ProbPoint;                       % 所有智能体不符合条件的状态

        % 更新 Alpha、Beta 和 Delta 解
        if fitness <= Alpha_score  % 适应能力最强（因为性能指标越小越好，因此为小于号）
            Alpha_score = fitness;
            Alpha_pos = WolfPops.Pos(i, :);
            Alpha_no = i;
        end
        if fitness > Alpha_score && fitness <= Beta_score
            Beta_score = fitness;
            Beta_pos = WolfPops.Pos(i, :);
            Beta_no = i;
        end
        if fitness > Alpha_score && fitness > Beta_score && fitness <= Delta_score
            Delta_score = fitness;
            Delta_pos = WolfPops.Pos(i, :);
            Delta_no = i;
        end
    end

    % ②  更新参数a
    a = 2 - iter * 2 / Max_iter;                  % 线性递减
    %a = 2 * cos((iter / Max_iter) * pi/2);   % 非线性递减

    % ③  更新位置（朝着前三只狼位置前进）
    for i = 1 : SearchAgents
        for j = 1 : dim

            r1 = rand();
            r2 = rand();
            A1 = 2*a*r1 - a;
            C1 = 2*r2;
            D_alpha = abs(C1*Alpha_pos(j) - WolfPops.Pos(i, j));
            X1(i, j) = Alpha_pos(j) - A1*D_alpha;

            r1 = rand();
            r2 = rand();
            A2 = 2*a*r1 - a;
            C2 = 2*r2;
            D_beta = abs(C2*Beta_pos(j) - WolfPops.Pos(i, j));
            X2(i, j) = Beta_pos(j) - A2*D_beta;

            r1 = rand();
            r2 = rand();
            A3 = 2*a*r1 - a;
            C3 = 2*r2;
            D_delta = abs(C3*Delta_pos(j) - WolfPops.Pos(i, j));
            X3(i, j) = Delta_pos(j) - A3*D_delta;

            %更新
            %WolfPops.Pos(i, j) = (X1(i, j) + X2(i, j) + X3(i, j)) / 3;

        end

    end

    % ④  Cuckoo 搜索
    fitness = nan(SearchAgents, 1);
    for i = 1 : SearchAgents
        [fitness(i), ~, ~] = ObjFun(Tracks{i}, UAV);
    end
    [~, index] = min(fitness);
    best = WolfPops.Pos(index, :);
    X1 = get_cuckoos(X1, best, WolfPops.lb, WolfPops.ub);
    X2 = get_cuckoos(X2, best, WolfPops.lb, WolfPops.ub);
    X3 = get_cuckoos(X3, best, WolfPops.lb, WolfPops.ub);
    X1 = empty_nests(X1, WolfPops.lb, WolfPops.ub, pa);
    X2 = empty_nests(X2, WolfPops.lb, WolfPops.ub, pa);
    X3 = empty_nests(X3, WolfPops.lb, WolfPops.ub, pa);
    WolfPops.Pos = (X1 + X2 + X3) / 3;

    % ⑤  调整不符合要求的状态变量
    [WolfPops, Tracks] = BoundAdjust(WolfPops, ProbPoints, UAV);

    % ⑥  存储适应度
    Fitness_list(iter) = Alpha_score;

    if iter/Max_iter*100 < 10
        fprintf('\b\b\b\b\b%.2f%%', iter/Max_iter*100)
    else
        fprintf('\b\b\b\b\b\b%.2f%%', iter/Max_iter*100)
    end
end
fprintf('\n\n>>计算完成！\n\n')
toc

% 输出值
solution.method = 'GWO';                % 算法
solution.WolfPops = WolfPops;       % 所有解种群信息
solution.Tracks = Tracks;                  % 所有解航迹信息
solution.Fitness_list = Fitness_list;   % α解适应度曲线
solution.Alpha_Data = Data;            % α解的威胁信息
solution.Alpha_no = Alpha_no;        % α解的位置
solution.Beta_no = Beta_no;             % β解的位置
solution.Delta_no = Delta_no;          % δ解的位置

end

%%%%% 布谷鸟搜索 %%%%%
function nest = get_cuckoos(nest, best, Lb, Ub)
    n = size(nest, 1);
    beta = 3/2;
    sigma = (gamma(1 + beta) * sin(pi*beta/2)/(gamma((1 + beta)/2)*beta*2^((beta - 1)/2)))^(1/beta);
    for j = 1:n
        s = nest(j, :);
        u = randn(size(s))*sigma;
        v = randn(size(s));
        step = u./abs(v).^(1/beta);
        stepsize = 0.01*step.*(s - best);
        s = s + stepsize.*randn(size(s));
        nest(j, :) = BoundClamp(s, Lb, Ub);
    end
end
function new_nest = empty_nests(nest, Lb, Ub, pa)
    n = size(nest, 1);
    K = rand(size(nest)) > pa;
    stepsize = rand*(nest(randperm(n), :) - nest(randperm(n), :));
    new_nest = nest + stepsize.*K;
    for j = 1:size(new_nest, 1)
        s = new_nest(j, :);
        new_nest(j, :) = BoundClamp(s, Lb, Ub);
    end
end
function x = BoundClamp(x, lb, ub)
    Flag4ub = x > ub;
    Flag4lb = x < lb;
    x = x .* ( ~(Flag4ub + Flag4lb) ) + ub .* Flag4ub + lb .* Flag4lb;
end

repo: Grey-Wolf-Optimizer-for-Path-Planning


text:
function seed = get_seed ( )

%*****************************************************************************80
%
%% GET_SEED returns a random seed for the random number generator.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer SEED, a random seed value.
%
  time_array = clock;

  hour = time_array(4);
  minute = time_array(5);
  second = time_array(6);

  seed = second + 60 * ( minute + 60 * hour );
%
%  We want values in [1,43200], not [0,43199].
%
  seed = seed + 1;
%
%  Remap SEED from [1,43200] to [1,I_MAX].
%
  i4_huge = 2147483647;
  seed = i4_huge * ( seed  / ( 60.0 * 60.0 * 24.0 ) );

  seed = floor ( seed );
%
%  Never use a seed of 0.
%
  if ( seed == 0 )
    seed = 1;
  end

  return
end

repo: jburkardt-m


text:
function ebsd = fill(ebsd,varargin)
% fill EBSD data by nearest neighbour
%
% Syntax
%   ebsd_filled = fill(ebsd)
%   ebsd_filled = fill(ebsd,grains)
%
% Input
%  ebsd - @EBSD
%  grains - @grain2d
%
% Options
%  extrapolate - extrapolate up the the outer boundaries
%

if ~(isa(ebsd,'EBSDsquare') || isa(ebsd,'EBSDhex')), ebsd = ebsd.gridify; end

% the values to be filled
nanId = isnan(ebsd.phaseId);

if check_option(varargin,'extrapolate')
  opt = 'nearest';
else
  opt = 'none';
end

F = scatteredInterpolant([ebsd.prop.x(~nanId),ebsd.prop.y(~nanId)],...
  find(~nanId),'nearest',opt);

newId = F(ebsd.prop.x(nanId),ebsd.prop.y(nanId));

nanId(nanId) = ~isnan(newId);
newId(isnan(newId)) = [];

% interpolate phaseId
ebsd.phaseId(nanId) = ebsd.phaseId(newId);
ebsd.rotations(nanId) = ebsd.rotations(newId);

% interpolate grainId
try
  ebsd.prop.grainId(nanId) = ebsd.prop.grainId(newId);
end

grains = getClass(varargin,'grain2d',[]);
if isempty(grains), return; end

grains = grains(ismember(grains.id,unique(ebsd.grainId)));

nanId = find(nanId);

% check for whether the pixels are within certain grains
isInside = checkInside(grains,ebsd.subSet(nanId));

% set phase to not indexed if not inside any grain
ebsd.phaseId(nanId(~any(isInside,2))) = 1;
ebsd.grainId(nanId(~any(isInside,2))) = 0;

% the values to be filled
[ebsdId,hostId] = find(isInside);

wrongGrainId = ebsd.grainId(nanId(ebsdId)) ~= grains.id(hostId);

ebsd.phaseId(nanId(ebsdId)) = grains.phaseId(hostId);
ebsd.grainId(nanId(ebsdId)) = grains.id(hostId);
ebsd.rotations(nanId(ebsdId(wrongGrainId))) = grains.meanRotation(hostId(wrongGrainId));

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% Cobb-Douglas production function
% Y=F(K,N)=K^{alpha}*N^{1-alpha}
% In per capita form y=k^{alpha}

function y=cobb(k)
global alpha
if alpha>0 & alpha<1
   y=k.^alpha;
else
   y=-inf*abs(k);
   disp('invalid capital share')
end


text:
%CANNY2  Finds edges in an image using the Canny algorithm with custom image gradient
%
%     edges = cv.Canny2(dx, dy, thresh)
%     edges = cv.Canny2(dx, dy, thresh, 'OptionName', optionValue, ...)
%
% ## Input
% * __dx__ 16-bit x derivative of input image (1 or 3 channels of type `int16`).
% * __dy__ 16-bit y derivative of input image (same size and type as `dx`).
% * __thresh__ Threshold for the hysteresis procedure. Scalar or 2-element
%   vector `[low_thresh,high_thresh]`.
%
% ## Output
% * __edges__ Output edge map; single channels 8-bit image, which has the same
%   size as the input image.
%
% ## Options
% * __L2Gradient__ Flag indicating whether a more accurate L2 norm
%   `sqrt((dI/dx)^2 + (dI/dy)^2)` should be used to compute the image gradient
%   magnitude (`L2gradient=true`), or whether the default L1 norm
%   `abs(dI/dx) + abs(dI/dy)` is enough (`L2gradient=false`). Default false
%
% See also: cv.Canny, cv.Sobel, cv.Scharr
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sz = sizeJava(jv)
% sizeJava(jv) returns the size of a java array, in a manner corresponding
% to 'size' on matlab arrays.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

sz(1) = size(jv,1);
sz(2) = size(jv(1,:),1);
if sz(2) == 0
  % workaround matlab's bug
  sz(2) = 1;
end

author: andrewssobral

repo: lrslibrary