text:
%% PASCAL3D evaluation
% This script runs evaluation on PASCAL3D.
% We assume that the heatmaps (keypoint localizations)
% are already computed and located in the folder 'predpath'.

clear
startup

% path for annotations
datapath = 'pose-hg/pose-hg-demo/data/pascal3d/annot/';
% path for network output
predpath = 'pose-hg/pose-hg-demo/exp/pascal3d/';
% path where the results are stored
savepath = 'result/cad/';
mkdir(savepath);
annotfile = sprintf('%s/valid.mat',datapath);
load(annotfile);

% determines shape model used for pose optimization
% if set to 1 then we use the particular cad model instance
% if set to 0 then we use the deformable shape model
cadSpecific = 1;

% visualization flag
vis = 0;

% evaluate only on the clean set
testlist = find(~annot.occluded & ~annot.truncated);

for ID = testlist'

    % input
    imgname = annot.imgname{ID};
    center = annot.center(ID,:);
    scale = annot.scale(ID);
    class = annot.class{ID};
    indices = annot.indices{ID};
    cadID = annot.cad_index(ID);

    cad = load(sprintf('cad/%s.mat',class));
    cad = cad.(class);
    cad = cad(cadID);

    savefile = sprintf('%s/valid_%d.mat',savepath,ID);

    if cadSpecific
        dict = getPascalTemplate(cad);
    else
        dict = load(sprintf('dict/pca-%s.mat',class));
    end
    % read heatmaps and detect maximum responses
    heatmap = h5read(sprintf('%s/valid_%d.h5',predpath,ID),'/heatmaps');
    heatmap = permute(heatmap(:,:,indices(dict.kpt_id)),[2,1,3]);
    [W_hp,score] = findWmax(heatmap);

    % pose optimization - weak perspective
    output_wp = PoseFromKpts_WP(W_hp,dict,'weight',score,'verb',false,'lam',1);

    % visualization
    if vis
        img = imread(sprintf('%s/../images/%s.jpg',datapath,imgname));
        vis_wp(img,output_wp,heatmap,center,scale,cad,dict);
        pause
        close all
    end

    % save output
    save(savefile,'output_wp');
end

% results
pascal3d_res(datapath, savepath);

author: geopavlakos


text:
function f = p17_f ( n, x )

%*****************************************************************************80
%
%% P17_F evaluates the objective function for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Brent,
%    Algorithms for Minimization with Derivatives,
%    Dover, 2002,
%    ISBN: 0-486-41998-3,
%    LC: QA402.5.B74.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  f1 = x(2) - x(1) * x(1);
  f2 = 1.0 - x(1);
  f3 = x(4) - x(3) * x(3);
  f4 = 1.0 - x(3);
  f5 = x(2) + x(4) - 2.0;
  f6 = x(2) - x(4);

  f = 100.0 * f1 * f1 ...
    +         f2 * f2 ...
    +  90.0 * f3 * f3 ...
    +         f4 * f4 ...
    +  10.0 * f5 * f5 ...
    +   0.1 * f6 * f6;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [b,a,fx,bx,gd]=gammabank(n,fs,w,fc,bw,ph,k)
%GAMMABANK gammatone filter bank [b,a,fx,bx,gd]=(n,fs,w,fc,bw,ph,k)
%
% Usage:
%          (1) [b,a,fx,bx,gd]=gammabank(0.35,fs,'',[100 6000]);
%              y = filterbank(b,a,s,gd);
%
%              Will create an erb-spaced filterbank between 100 Hz and 6kHz
%              with a filter spacing of 0.35 erb and a default bandwidth
%              of 1.019 erb. Omitting the "y =" from the second line will plot
%              a spectrogram.
% Inputs:
%       n   number of filters in filterbank or the filter spacing in
%           bark/erb/k-mel. Set n=0 if fc lists centre frequencies explicitly.
%		fs  sample rate in Hz
%		fc  centre frequencies [default = [100 6000] ]
%		bw  bandwidths [default = 1.019*erb(fc) ]
%       ph  phase of gammatone impulse repsonse [default = 0]
%       k   filter ordere [default = 4]
%		w   any sensible combination of the following:
%             'e' = erb scale for filter spacing, frequencies ('F' option)
%             and bandwidths ('W' option but can be overridden by 'EMBLH')
%                   'm','b','l','h' = mel, bark, log10 and Hz scale
%             'E' = erb scale for bandwidths ('W' option)
%                   'M','B','L','H' = mel, bark, log10 and Hz scale
%
%             'n' = n input gives number of filters [default if n>=1]
%             'N' = n input gives filter spacing  [default if n<1]
%
%             'f' = fc is in Hz [default]
%             'F' = fc is in mel/erb-rate/bark/log10
%             'w' = bw inputs are in Hz [default]
%             'W' = bw inputs are in multiples of df/dx where x=bark/erb/mel etc
%
%             'k' = force a filter at 1kHz
%             ['d' = choose ph() so that all filters have zero DC gain]
%             ['a' = use all-pole gammtone funtion: see [1]]
%             ['s' = use Slaney gammatone approximation: see [2]]
%             ['z' = use one-zero gammatone function: see [1]]
%
%             'g' = plot filter responses [default if no output arguments present]
%             'G' = plot frequency responses on a log axis
%
% Outputs:
%       b/a    filter coefficients: one filter per row
%       fx,bx  centre frequencies and bandwidths in Hz
%       gd     group delay at the centre frequencies (in samples)
%

% The impulse response of filter i is proportional to:
%       h(n)=((n/fs).^(k-1))*cos(2*pi*fx(i)*n/fs+ph(i))*exp(-2*pi*bx(i)*n/fs)
% where n=0,1,2,...
% Note that the DC gain is only equal to zero for one particular value of ph(i)
% The filters are normalized to have unity gain at the centre frequencies
%
% References
%  [1]	R. F. Lyon, A. G. Katsiamis, and E. M. Drakakis.
%       History and future of auditory filter models.
%       In Proc Intl Symp Circuits and Systems, pages 38093812, 2010.
%       doi: 10.1109/ISCAS.2010.5537724.
%  [2]	M. Slaney.
%       An efficient implementation of the patterson-holdsworth auditory filter bank.
%       Technical report, Apple Computer, Perception Group, Tech. Rep, 1993.

%      Copyright (C) Mike Brookes 2009-2010
%      Version: $Id: gammabank.m 2190 2012-07-20 13:47:42Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<7
    k=[];
    if nargin<6
        ph=[];
        if nargin<5
            bw=[];
            if nargin<4
                fc=[];
                if nargin<3
                    w='';
                end
            end
        end
    end
end
fx=fc(:);
bx=bw(:);
if ~numel(k)
    k=4;
end
if ~numel(fx)
    fx=[100; 6000]; % default
end
wr='e';   % default frequency warping is erb
for i=1:length(w)
    if any(w(i)=='bmlef');
        wr=w(i);
    end
end
if any(w=='k')
    fk=1000;
    switch wr              % convert 1kHz to spacing units
        case 'b'
            fk=frq2bark(fk);
        case 'l'
            fk=log10(fk);
        case 'e'
            fk=frq2erb(fk);
    end
else
    fk=0;
end
if any(w)=='W'
    wb=wr;
else
    wb='h';     % default bandwidth units are Hz
end
for i=1:length(w)
    if any(w(i)=='BMLEF');
        wb=w(i)+'a'-'A';        % convert to lower case
    end
end
if ~numel(bx)
    bx=1.019;
    wb='e';
end
if any(w)=='F'          % convert centre frequencies to Hz
    switch wr
        case 'b'
            fx=bark2frq(fx);
        case 'm'
            fx=mel2frq(fx);
        case 'l'
            fx=10.^(fx);
        case 'e'
            fx=erb2frq(fx);
    end
end

% now sort out the centre frequencies

if n>0                      % n>0: filter end points specified
    bx=bx(1);               % only use the first bx element
    if n==1             % only one filter requested
        fx=fx(1);           % just use the first frequency
    else
        switch wr               % convert end frequencies to spacing units
            case 'b'
                fx=frq2bark(fx);
            case 'm'
                fx=frq2mel(fx);
            case 'l'
                fx=log10(fx);
            case 'e'
                fx=frq2erb(fx);
        end
        if n<1 || any(w=='N')       % n = filter spacing
            if fk               % force filter to 1 kHz
                f0=fk-n*floor((fk-fx(1))/n);
            else                % centre filters in range
                f0=(fx(2)+fx(1)-n*floor((fx(2)-fx(1))/n))/2;
            end
            fx=(f0:n:fx(2))';

        else                        % n = number of filters specified
            % Multiple filters - evenly spaced
            fx=linspace(fx(1),fx(2),n)';     % centre frequencies in spacing units
            if fk              % force a filter at 1kHz
                ik=1+ceil((fk-fx(1))*(n-1)/(fx(n)-fx(1))); % index of centre freq immediately above 1 kHz
                if ik>n || ik>1 && ((fk-fx(1))*(fx(n)-fx(ik-1))>(fx(n)-fk)*(fx(ik)-fx(1)))
                    fx=fx(1)+(fx-fx(1))*(fk-fx(1))/(fx(ik)-fx(1));
                else
                    fx=fx(n)+(fx-fx(n))*(fx(n)-fk)/(fx(n)-fx(ik-1));
                end
            end
        end
        switch wr % convert back to Hz
            case 'b'
                fx=bark2frq(fx);
            case 'm'
                fx=mel2frq(fx);
            case 'l'
                fx=10.^(fx);
            case 'e'
                fx=erb2frq(fx);
        end
    end

end
% now sort out the bandwidths
nf=numel(fx);
if numel(bx)==1
    bx=bx(ones(nf,1));      % replicate if necessary
end
switch wb               % convert bandwidth to Hz
    case 'b'
        [dum,bwf]=frq2bark(fx);
    case 'm'
        [dum,bwf]=frq2mel(fx);
    case 'l'
        bwf=fx*log(10);
    case 'e'
        [dum,bwf]=frq2erb(fx);
    case 'h'
        bwf=ones(nf,1);
end
bx=bx.*bwf;
if ~numel(ph)
    ph=0;
end
if numel(ph)==1
    ph=ph(ones(nf,1));      % replicate if necessary
else
    ph=ph(:);
end
%
% t=(0:ceil(10*fs/(2*pi*bnd)))/fs;  % five time constants
% gt=t.^(n-1).*cos(2*pi*cfr*t+phi).*exp(-2*pi*bnd*t);
% gt=gt/sqrt(mean(gt.^2)); % normalize
% figure(1);
% plot(t,gt);
% title('Desired Impulse response');
% xlabel(['Time (' xticksi 's)']);
%
ww=exp((1i*fx-bx)*2*pi/fs);
a=round([1 cumprod((-k:-1)./(1:k))]);   % create binomial coefficients
b=conv(a,(0:k-1).^(k-1));
b=exp(1i*ph)*b(1:k);
wwp=repmat(ww,1,k+1).^repmat(0:k,nf,1);
denc=repmat(a,nf,1).*wwp;
numc=b.*wwp(:,1:k);
b=zeros(nf,2*k);
a=zeros(nf,2*k+1);
gd=zeros(nf,1);
ww=exp(2i*fx*pi/fs);
for i=1:nf
    b(i,:)=real(conv(numc(i,:),conj(denc(i,:))));
    a(i,:)=real(conv(denc(i,:),conj(denc(i,:))));
    u=polyval(b(i,:),ww(i));
    v=polyval(a(i,:),ww(i));
    ud=polyval(b(i,:).*(0:2*k-1),ww(i));
    vd=polyval(a(i,:).*(0:2*k),ww(i));
    b(i,:)=b(i,:)*abs(v/u);
    gd(i)=real((v*ud-u*vd)/(u*v));     % group delay at centre freq in samples
end

% now plot graph

if ~nargout || any(w=='g') || any(w=='G')
    ng=200;      %number of points to plot
    if any(w=='G')
        fax=logspace(log10(fx(1)/4),log10(fs/2),ng);
    else
        fax=linspace(0,fs/2,ng);
    end
    ww=exp(2i*pi*fax/fs);
    gg=zeros(nf,ng);
    for i=1:nf
        gg(i,:)=10*log10(abs(polyval(b(i,:),ww)./polyval(a(i,:),ww)));
    end
    if any(w=='G')
        semilogx(fax,gg','-b');
        set(gca,'xlim',[fax(1) fax(end)]);
    else
        plot(fax,gg','-b');
    end
    xlabel(['Frequency (' xticksi 'Hz)']);
    set(gca,'ylim',[-50 1]);
    title(sprintf('%d Gammatone Filters (Opt=%s)',nf,w));
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
classdef H0Kernel < Kernel
    properties
        k;
        C; % Such that G(r) = C*besselh(0,1,k*r)
        % Hankel function of the first kind.
    end
    methods
        function[this] = H0Kernel(kk,CC)
            if nargin == 0
                kk = 1;
            end
            if nargin <= 1
                CC = 1;
            end
            this.k = kk;
            this.C = CC;
            Jk = J0Kernel(kk);
            Yk = Y0Kernel(kk);
            modelKern = this.C*(Jk + 1i*Yk);
            this.func = modelKern.func;
            this.der = modelKern.der;
            this.scalFunc = modelKern.scalFunc;
            this.lim0 = this.C*(Jk.lim0 + 1i*Yk.lim0);
        end
        function[out] = dilatation(this,lambda)
            lim0keep = this.lim0;
            out = H0Kernel(this.k*lambda,this.C);
            out.lim0 = lim0keep;
        end
        function[out] = mtimes(this,mu)
            if isa(this,'Kernel')
                assert(and(isa(mu,'double'),isscalar(mu)));
                out = H0Kernel(this.k,this.C*mu);
            else
                out = mtimes(mu,this);
            end
        end
        function[rq] = radialQuadKernel(this,a,tol,varargin)
            modelKern = J0Kernel(this.k) + 1i*Y0Kernel(this.k);
            rqTemp = modelKern.radialQuadKernel(a,tol/abs(this.C),varargin{:});
            rq = this.C*rqTemp;
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function s = logsumexp(a, dim)
% Returns log(sum(exp(a),dim)) while avoiding numerical underflow.
% Default is dim = 1 (columns).
% logsumexp(a, 2) will sum across rows instead of columns.
% Unlike matlab's "sum", it will not switch the summing direction
% if you provide a row vector.

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

if nargin < 2
  dim = 1;
end

% subtract the largest in each column
[y, i] = max(a,[],dim);
dims = ones(1,ndims(a));
dims(dim) = size(a,dim);
a = a - repmat(y, dims);
s = y + log(sum(exp(a),dim));
i = find(~finite(y));
if ~isempty(i)
  s(i) = y(i);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%CODEGENERATOR.GENMFUNFKINE Generate M-function for forward kinematics
%
% cGen.genmfunfkine() generates a robot-specific M-function to compute
% forward kinematics.
%
% Notes::
% - Is called by CodeGenerator.genfkine if cGen has active flag genmfun
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.genjacobian.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [] = genmfunfkine(CGen)

%% Does robot class exist?
if ~exist(fullfile(CGen.robjpath,[CGen.getrobfname,'.m']),'file')
    CGen.logmsg([datestr(now),'\tCreating ',CGen.getrobfname,' m-constructor ']);
    CGen.createmconstructor;
    CGen.logmsg('\t%s\n',' done!');
end

%% Forward kinematics up to tool center point
CGen.logmsg([datestr(now),'\tGenerating forward kinematics m-function up to the end-effector frame: ']);
symname = 'fkine';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genmfunfkine:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

funfilename = fullfile(CGen.robjpath,[symname,'.m']);
q = CGen.rob.gencoords;

matlabFunction(tmpStruct.(symname).T,'file',funfilename,...              % generate function m-file
    'outputs', {symname},...
    'vars', {'rob',[q]});
hStruct = createHeaderStructFkine(CGen.rob,symname);                 % replace autogenerated function header
replaceheader(CGen,hStruct,funfilename);
CGen.logmsg('\t%s\n',' done!');

%% Individual joint forward kinematics
CGen.logmsg([datestr(now),'\tGenerating forward kinematics m-function up to joint: ']);
for iJoints=1:CGen.rob.n

    CGen.logmsg(' %i ',iJoints);
    symname = ['T0_',num2str(iJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    tmpStruct = struct;
    tmpStruct = load(fname);

    funfilename = fullfile(CGen.robjpath,[symname,'.m']);
    q = CGen.rob.gencoords;

    matlabFunction(tmpStruct.(symname).T,'file',funfilename,...              % generate function m-file
        'outputs', {symname},...
        'vars', {'rob',[q]});
    hStruct = createHeaderStruct(CGen.rob,iJoints,symname);                 % replace autogenerated function header
    CGen.replaceheader(hStruct,funfilename);
end
CGen.logmsg('\t%s\n',' done!');

end

%% Definition of the header contents for each generated file
function hStruct = createHeaderStruct(rob,curBody,fname)
[~,hStruct.funName] = fileparts(fname);
hStruct.shortDescription = ['Forward kinematics for the ',rob.name,' arm up to frame ',int2str(curBody),' of ',int2str(rob.n),'.'];
hStruct.calls = {['T = ',hStruct.funName,'(rob,q)'],...
    ['T = rob.',hStruct.funName,'(q)']};
hStruct.detailedDescription = {['Given a set of joint variables up to joint number ',int2str(curBody),' the function'],...
    'computes the pose belonging to that joint with respect to the base frame.'};
hStruct.inputs = {['q:  ',int2str(curBody),'-element vector of generalized coordinates.'],...
    'Angles have to be given in radians!'};
hStruct.outputs = {['T:  [4x4] Homogenous transformation matrix relating the pose of joint ',int2str(curBody),' of ',int2str(rob.n)],...
    '          for the given joint values to the base frame.'};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
     'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {rob.name};
end

%% Definition of the header contents for each generated file
function hStruct = createHeaderStructFkine(rob,fname)
[~,hStruct.funName] = fileparts(fname);
hStruct.shortDescription = ['Forward kinematics solution including tool transformation for the ',rob.name,' arm.'];
hStruct.calls = {['T = ',hStruct.funName,'(rob,q)'],...
    ['T = rob.',hStruct.funName,'(q)']};
hStruct.detailedDescription = {['Given a full set of joint variables the function'],...
    'computes the pose belonging to that joint with respect to the base frame.'};
hStruct.inputs = { ['rob: robot object of ', rob.name, ' specific class'],...
                   ['q:  ',int2str(rob.n),'-element vector of generalized'],...
                   '     coordinates',...
                   'Angles have to be given in radians!'};
hStruct.outputs = {['T:  [4x4] Homogenous transformation matrix relating the pose of the tool'],...
    '          for the given joint values to the base frame.'};
hStruct.references = {'1) Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    '2) Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    '3) Introduction to Robotics, Mechanics and Control - Craig',...
    '4) Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn',...
    '2012 RST, Technische Universitaet Dortmund, Germany',...
     'http://www.rst.e-technik.tu-dortmund.de'};
hStruct.seeAlso = {'jacob0'};
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [hdr, mrGrayFileName] = analyze2mrGray(analHeaderFileName,endianType);
% function [V,mrGrayFileName]=analyze2mrGray(analHeaderFileName,endianType);
% Converts Analyze7.5 images to Stanford VISTA lab / mrGray vAnatomy.dat format.
% Analyze headers hold far more information that mrGray volumes. We just toss this
% information out when we write the mrGray file, but this routine returns the
% Analyze header info in 'V'.
% Uses the SPM99 spm_vol command to parse the header. Uses misc. VISTA lab routines
% to clip and write the volume.
% See also mrGray2Analyze
% If you pass in the (optional) endianType it will use this when reading 16bit data
% files.
%
% HISTORY:
% ARW 081701
% ARW Added >16 bit data handling 082101
% 02.07.16 RFD (bob@white.stanford.edu) Rewrote code in a more modular format.
% Also added proper image reorientation to make SPM-style analyze format
% files into mrGray orientation.

if (~exist('endianType','var'))
    endianType='';
end
if (~exist('analHeaderFileName','var'))
    analHeaderFileName = '';
end
[img, mmPerPix, hdr] = loadAnalyze(analHeaderFileName, endianType);

% to make analyze data the same orientation as vAnatomy, we reorient the data.
% Note that we assume Analyze orientation code '0'- transverse, unflipped.
% If that is the case, then the following will make the data look correct in mrGray.
% It will even get left/right correct (left will be left and right right).
img = permute(img,[3,2,1]);

mmPerPix = [mmPerPix(3),mmPerPix(2),mmPerPix(1)];
% flip each slice ud (ie. flip along matlab's first dimension, which is our x-axis)
for(jj=1:size(img,3))
    img(:,:,jj) = flipud(squeeze(img(:,:,jj)));
end
% flip each slice lr(ie. flip along matlab's second dimension, which is our y-axis)
for(jj=1:size(img,3))
    img(:,:,jj) = fliplr(squeeze(img(:,:,jj)));
end
% flip along matlab's third dimension, which is our z-axis
% NO NEED TO DO THIS IF DATA ARE IN STANDARD ANALYZE ORIENTATION
% for(jj=1:size(img,1))
%     img(jj,:,:) = fliplr(squeeze(img(jj,:,:)));
% end

% launch a GUI to find the optimal clip values
img = makeClippedVolume(img);

% Now write out a mrGray file
mrGrayFileName = writeVolAnat(img, mmPerPix);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef ParameterAnalysis
    %PARAMETERANALYSIS Summary of this class goes here
    %   Detailed explanation goes here

    properties
        % 1 for multiplication, 2 for addition
        StepSizeMode
    end

    properties
        % 1 for 1-order.
        % 2 for 2-order (cartesian product).
        % A vector to indicate size for each parameter, e.g. (5,5) means
        % 5 neighbors (including current value) for each of the 2
        % parameters.
        NeighborhoodMode
    end

    properties
        PrintListSearch
    end

    properties
        % Use 'Greedy' or 'BrutalForce' or 'BlockCoordinateDescent'
        Algorithm
    end

    %% parameters used by brutal force search
    properties
        BrutalForceDepth
    end

    %% parameters used by BCD
    properties
        BlockCoordinateDescentDepth
    end

    properties
        BlockCoordinateDescentGroup
    end

    properties
        MaxNumberOfIterations
    end

    %% parameters used by all
    properties
        CheckValueInListSearchOptions
    end

    methods
        function obj = ParameterAnalysis()
            obj.StepSizeMode = 1;
            obj.NeighborhoodMode = 1;
            obj.PrintListSearch = 0;
            obj.Algorithm = 'Greedy';
            obj.BrutalForceDepth = 3;
            obj.CheckValueInListSearchOptions = 1;

            obj.MaxNumberOfIterations = 50;
            obj.BlockCoordinateDescentDepth = 5;
        end

        function options = autoParamSelection(obj, fcn_run, options, ...
                list_params, step_size, range)
            disp('Start searching best parameters. This could take a while.');

            switch obj.Algorithm
                case 'Greedy'
                    options = autoParamSelectionByGreedyAlgo(obj, ...
                        fcn_run, options, list_params, step_size, range);
                case 'BrutalForce'
                    options = autoParamSelectionByBrutalForce(obj, ...
                        fcn_run, options, list_params, step_size, range);
                case 'BlockCoordinateDescent'
                    options = autoParamSelectionByBCD(obj, ...
                        fcn_run, options, list_params, step_size, range);
                otherwise
            end
        end

        function options = autoParamSelectionByGreedyAlgo(obj, fcn_run, ...
                options, list_params, step_size, range)
%             if 0
            curr_val = fcn_run(options);
            list_search_options = options;
            list_search_value = curr_val;
            while 1
                nbhds = obj.findNeighboringOptions(options, list_params, ...
                    step_size, range, obj.NeighborhoodMode);
                nbhds_vals = zeros(1,length(nbhds));
                for i = 1:length(nbhds)
                    [val, list_search_options, list_search_value] = ...
                        obj.run(fcn_run, list_params, nbhds(i), ...
                        list_search_options, list_search_value);
                    nbhds_vals(i) = val;
                end

                if obj.testLocalMinimum(curr_val, nbhds_vals)
                    break;
                else
                    [min_val,ind] = min(nbhds_vals);
                    options = nbhds(ind);
                    curr_val = min_val;
                end

                if obj.PrintListSearch
                    obj.printListSearch(list_search_options, list_search_value, ...
                        list_params);
                end
            end
            save('param_analysis.mat','list_search_options',...
                'list_search_value','options');
%             end

            disp('The best searched parameters are');
            options
        end

        function options = autoParamSelectionByBrutalForce(obj, fcn_run, ...
                options, list_params, step_size, range)
            list_search_options = [];
            list_search_value = [];
            for k = 1:obj.BrutalForceDepth
                nbhds = obj.findNeighboringOptions(options, list_params, ...
                    step_size, range, obj.NeighborhoodMode);
                nbhds_vals = zeros(1,length(nbhds));
                for i = 1:length(nbhds)
                    [val, list_search_options, list_search_value] = ...
                        obj.run(fcn_run, list_params, nbhds(i), ...
                        list_search_options, list_search_value);
                    nbhds_vals(i) = val;
                end

                [min_val,ind] = min(nbhds_vals);
                options = nbhds(ind);
                step_size = step_size / 2;
            end
            save('param_analysis.mat','list_search_options',...
                'list_search_value','options');

            if obj.PrintListSearch
                obj.printListSearch(list_search_options, list_search_value, ...
                    list_params);
            end
            disp('The best searched parameters are');
            options

        end

        function options = autoParamSelectionByBCD(obj, ...
                fcn_run, options, list_params, step_sizes, range)
            last_options = options;
            for iter = 1:obj.MaxNumberOfIterations
                for j = 1:length(obj.BlockCoordinateDescentGroup)
                    inds = obj.BlockCoordinateDescentGroup{j};
                    list_search_options = [];
                    list_search_value = [];
                    step_size = step_sizes(inds);
                    for k = 1:obj.BlockCoordinateDescentDepth
                        nbhds = obj.findNeighboringOptions(options, list_params(inds), ...
                            step_size, range(:,inds), obj.NeighborhoodMode(inds));
                        nbhds_vals = zeros(1,length(nbhds));
                        for i = 1:length(nbhds)
                            [val, list_search_options, list_search_value] = ...
                                obj.run(fcn_run, list_params, nbhds(i), ...
                                list_search_options, list_search_value);
                            nbhds_vals(i) = val;
                        end

                        [min_val,ind] = min(nbhds_vals);
                        options = nbhds(ind);
                        step_size = step_size / 2;
                    end
                end
                if ParameterAnalysis.isOptionsEqual(options, last_options, list_params)
                    disp(['BCD stopped at iteration ',num2str(iter)]);
                    break;
                else
                    last_options = options;
                end
            end

            disp('The best searched parameters are');
            options
        end

        function printListSearch(obj, list_search_options, list_search_value,...
                list_params)
            disp('--------------------- Search Results -----------------');
            list_params = [list_params,{'value'}];
            disp(list_params);

            N = length(list_search_options);
            M = length(list_params);
            X = zeros(N, M);
            for i = 1:N
                options = list_search_options(i);
                value = list_search_value(i);
                for j = 1:M-1
                    name = list_params{j};
                    X(i,j) = options.(name);
                end
                X(i,M) = value;
            end
            disp(X);
        end

        function [val, list_search_options, list_search_value] = run( ...
                obj, fcn_run, list_params, options, list_search_options,...
                list_search_value)
            find = 0;
            if obj.CheckValueInListSearchOptions
                for i = 1:length(list_search_options)
                    if ParameterAnalysis.isOptionsEqual(options, ...
                            list_search_options(i), list_params)
                        find = 1;
                        val = list_search_value(i);
                        break;
                    end
                end
            end
            if ~find
                val = fcn_run(options);
                list_search_options = [list_search_options, options];
                list_search_value = [list_search_value, val];
            end
        end

        function T = testLocalMinimum(obj, curr_val, nbhds_vals)
            T = all(curr_val <= nbhds_vals);
        end

        function nbhds = findNeighboringOptions(obj, options, list_params, ...
                step_size, range, nbhd_mode)
            if length(obj.StepSizeMode) == 1
                step_mode = repmat(obj.StepSizeMode, [length(step_size),1]);
            else
                step_mode = obj.StepSizeMode;
            end

            if isscalar(nbhd_mode) && nbhd_mode == 1
                % 1 order neighborhood will not consider the current point
                nbhds = find1OrderNeighborhood(obj, options, ...
                    list_params, step_size, range, step_mode);
            elseif isscalar(nbhd_mode) && nbhd_mode == 2
                nbhds = find2OrderNeighborhood(obj, options, ...
                    list_params, step_size, range, step_mode);
            else % it's a vector specifying the size of nbhd
                nbhds = findBigNeighborhood(obj, options, ...
                    list_params, step_size, range, step_mode, ...
                    nbhd_mode);
            end
        end

        function nbhds = find1OrderNeighborhood(obj, options, ...
                list_params, step_size, range, modes)
            nbhds = [];
            options_start = options;
            for i = 1:length(list_params)
                options = options_start;
                name = list_params{i};
                forw_val = obj.findNeighboringValue(options_start.(name),...
                    step_size, i, 1, modes, range);
                back_val = obj.findNeighboringValue(options_start.(name),...
                    step_size, i, 0, modes, range);

                options.(name) = forw_val;
                nbhds = [nbhds, options];
                options.(name) = back_val;
                nbhds = [nbhds, options];
            end
        end

        function nbhds = find2OrderNeighborhood(obj, options, ...
                list_params, step_size, range, modes)
            % modes are the step mode (multiplication or addition)
            nbhd_size = 3 * ones(size(list_params));
            nbhds = findBigNeighborhood(obj, options, ...
                list_params, step_size, range, modes, nbhd_size);
        end

        function nbhds = findBigNeighborhood(obj, options, ...
                list_params, step_size, range, modes, nbhd_size)
            % num_steps takes the form like (3,3) for 2 parameters where
            % 3 values including the current value are considered as nbhd
            num_steps = floor(nbhd_size/2);
            options_start = options;
            params = cell(1,length(list_params));
            for i = 1:length(list_params)
                name = list_params{i};
                curr_val = options_start.(name);
                vals = curr_val;
                for j = 1:num_steps(i)
                    back_val = obj.findNeighboringValue(vals(1), ...
                        step_size, i, 0, modes, range);
                    vals = cat(1,back_val,vals);
                end
                for j = 1:num_steps(i)
                    forw_val = obj.findNeighboringValue(vals(end), ...
                        step_size, i, 1, modes, range);
                    vals = cat(1,vals,forw_val);
                end
                params{i} = vals;
            end
            % cartprod can't handle duplicate entries
%             C = cartprod(params{:});
            C = ParameterAnalysis.cartesian(params{:});

            nbhds = repmat(options_start, [size(C,1),1]);
            for i = 1:size(C,1)
                for j = 1:length(list_params)
                    nbhds(i).(list_params{j}) = C(i,j);
                end
            end
        end

        function new_val = findNeighboringValue(obj, curr_val, ...
                step_size, i, direction, modes, range)
            order_of_magnitude = @(x) floor(log10(x));

            if isnumeric(step_size)
                size = step_size(i);
                mode = modes(i);
                new_val = obj.calcNewValue(direction, curr_val, size, mode);
            elseif isequal(step_size,'medium')
                possible_vals = [1,2,5];
                size = 2;
                mode = 1;
                new_val = obj.calcNewValue(direction, curr_val, size, mode);
                order = order_of_magnitude(new_val);
                new_val = new_val/(10^order);
                [~,ind] = min(abs(new_val - possible_vals));
                new_val = possible_vals(ind);
                new_val = new_val * (10^order);
            end

            if direction == 1
                upper_bd = range(2,i);
                new_val = min(new_val, upper_bd);
            else
                lower_bd = range(1,i);
                new_val = max(new_val, lower_bd);
            end
        end

        function new_val = calcNewValue(obj, direction, curr_val, size, mode)
            if mode == 1 % multiplication
                if direction == 1
                    new_val = curr_val * size;
                else
                    new_val = curr_val / size;
                end
            elseif mode == 2 % addition
                if direction == 1
                    new_val = curr_val + size;
                else
                    new_val = curr_val - size;
                end
            end
        end
    end
    methods (Static)
        function C = cartesian(varargin)
            args = varargin;
            n = nargin;

            [F{1:n}] = ndgrid(args{:});

            for i=n:-1:1
                G(:,i) = F{i}(:);
            end

            C = unique(G , 'rows');
        end

        function t = isOptionsEqual(options1, options2, list_params)
            t = 1;
            for i = 1:length(list_params)
                if options1.(list_params{i}) ~= options2.(list_params{i})
                    t = 0;
                    break;
                end
            end
        end
    end
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [y,s]=gammalns(x)
% GAMMALNS Log of Gamma(x) for positive or negative real x [y,s]=(x)
%
%  Inputs: x      real matrix of input values
%
% Outputs: y      log(gamma(x)) or, if s is present, log(abs(gamma(x)))
%          s      sign(gamma(x))
%
% If s is not specified then y will b complex if gamma(x) is negative (i.e. if floor(x) is negative and odd).
% Uses gamma(x) = pi/gamma(1-z)/sin(pi*z) which is 5.5.3 from [1]
%
% [1]	F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.
%       NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.

%      Copyright (C) Mike Brookes 2017
%      Version: $Id: gammalns.m 9975 2017-07-02 09:57:41Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent l
m=x<=0;
y=zeros(size(x));
s=ones(size(x));
if ~all(m(:)) % non-negative x values
    y(~m)=gammaln(x(~m));
end
f=m & (x==fix(x)); % find and non-positive integers
if any(f(:))
    y(f)=Inf; % set output to infinity
    m=m & ~f; % do not consider these values furhter
end
if any(m(:)) % negative x values
    if isempty(l)
        l=log(pi);
    end
    t=sin(pi*x(m));
    if nargout>1
        p=t<0;
        s(m)=1-2*(p);
        y(m)=l-gammaln(1-x(m))-log(abs(t));
    else
        y(m)=l-gammaln(1-x(m))-log(t);
    end
end

author: jtkim-kaist

repo: Speech-enhancement


text:
% design and visualize arrays
clear(); close all;

wavelength = 1; % normalized
d = wavelength / 2;
design_ula = design_array_1d('ula', 12, d);
design_cp = design_array_1d('coprime', [4 5], d);
design_nested = design_array_1d('nested', [5 7], d);
design_mra = design_array_1d('mra', 12, d);

visualize_array(design_ula, 'VisualizeCoarray', true);
visualize_array(design_cp, 'VisualizeCoarray', true);
visualize_array(design_nested, 'VisualizeCoarray', true);
visualize_array(design_mra, 'VisualizeCoarray', true);

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
classdef PTKCoords < handle
    % PTKCoords. Class for storing a set of coordinates.
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        Coords
    end

    methods
        function obj = PTKCoords(coords)
            if nargin > 0
                obj.Coords = coords;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef JointPD < Controller
    % This class defines a PD feedback control on the joint space
    %
    % @author ayonga @date 2016-10-14
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    methods

        function obj = JointPD(name)
            % The controller class constructor function
            %
            % Parameters:
            % name: the controller name @type char

            % call superclass constructor
            obj = obj@Controller(name);

            % initialize default control parameters
            ep = 10;
            obj.Param.kp = ep^2;
            obj.Param.kd = 2*ep;
        end

        function u = calcControl(obj, t, x, vfc, gfc, plant, params, logger)
            % Computes the PD feedback control law on joint space
            %
            % Parameters:
            % t: the time instant @type double
            % x: the states @type colvec
            % vfc: the vector field f(x) @type colvec
            % gfc: the vector field g(x) @type colvec
            % plant: the continuous domain @type DynamicalSystem
            % params: the control parameters @type struct
            % logger: the data logger object @type SimLogger
            %
            % Return values:
            % u: the computed torque @type colvec

            [qd, dqd] = calcDesiredStates(plant, t, x, params);

            % compute error terms
            qerr = qa - qd;
            dqerr = dqa - dqd;

            % feedback controller
            u = - obj.Param.kp*qerr - obj.Param.kd*dqerr;

            if ~isempty(logger)
                calc = logger.calc;

                calc.qerr = qerr;
                calc.dqerr = dqerr;
                calc.qd = qd;
                calc.dqd = dqd;
                calc.u = u;

                logger.calc = calc;
            end

        end
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [fhat, fola] = blocksyn(F, c , Lb, fola)
%BLOCKSYN Blockwise synthesis interface
%   Usage: blocksyn(F, c, Lb)
%
%   Input parameters:
%      F    : Synthesis frame object.
%      c    : Coefficients of a block.
%      Lb   : Length of the block.
%      fola : Explicitly defined overlap.
%   Output parameters:
%      fhat : Reconstructed block of signal.
%      fola : Stored overlap.
%
%   `fhat=blocksyn(F,c,Lb)` reconstructs the signal block *fhat* from the
%   coefficients *c* using the frame defined by *F*. If some overlap is used,
%   it is stored internally using block_interface. Note that the function is
%   capable of handling overlaps internally only for a single call
%   to the function in the blockproc. loop.%
%
%   `[fhat,fola]=blocksyn(F,c,Lb,fola)` does the same, but the block algorithm
%   uses `fola` to read and store overlap explicitly. `fola` can be empty.
%
%   **Note:** To get perfect reconstruction, the synthesis frame *F* must
%   be a dual frame of the analysis frame used in |blockana|.
%
%   See also: block, blockana, framedual
%
%   References: dogrhove12 ltfatnote026

complainif_notenoughargs(nargin,3,'BLOCKSYN');
complainif_notvalidframeobj(F,'BLOCKSYN');

    if ~isfield(F,'blockalg')
        F.blockalg = 'naive';
    end

    if ~isfield(F,'L')
         error(['%s: The frame object was not accelerated. See ',...
                'BLOCKFRAMEACCEL.'],upper(mfilename));
    end

    if nargin<4
       fola = [];
    end

    % Next block index start (from a global point of view, starting with zero)
    nextSb = block_interface('getPos');
    % Block index start (from a global point of view, starting with zero)
    Sb = nextSb-Lb;

    if strcmp(F.blockalg,'naive')
       if F.L < Lb
          error(['%s: The frame object was accelerated with incompatible ',...
                 'length. The block length is %i but the accelerated ',...
                 'length is %i.'],upper(mfilename),Lb,F.L);
       end
       % Most general. Should work for anything.
       % Produces awful block artifacts when coefficients are altered.
       fhat = F.frsyn(c);
       fhat = fhat(1:Lb,:);
    elseif strcmp(F.blockalg,'sliced')
        if F.L < 2*Lb
          error(['%s: The frame object was accelerated with incompatible ',...
                 'length. The block length is %i but the accelerated ',...
                 'length is %i.'],upper(mfilename),Lb,F.L);
        end
        % General processing
        % Equal block length assumption
        % Reconstruct
        f = F.frsyn(c);
        % If the transform length differs from the 2*Lb,
        % Pick the correct part from the result
        pad = size(f,1) - 2*Lb;
        sIdx = floor(pad/2);
        % Result should not be longer than 2*Lb
        f = f(1+sIdx:sIdx+2*Lb,:);
        % Multiply by a slicing window
        f = bsxfun(@times,F.sliwin,f);
        % Load and add overlap (first half)
        ol = loadOverlap(Lb,size(f,2),fola);
        olLen = size(ol,1);
        f(1:olLen,:) = f(1:olLen,:) + ol;
        % Store overlap (second half)
        if nargout>1
           fola=storeOverlap(f,Lb);
        else
           storeOverlap(f,Lb);
        end
        % Return first half
        fhat = f(1:Lb,:);
    elseif strcmp(F.blockalg,'segola')
       if ~isfield(F,'winLen')
          error('%s: Frame does not have FIR windows.',upper(mfilename));
       end
       Lw = F.winLen;
       switch(F.type)
         case 'fwt'
           % The SegDWT algorithm
           J = F.J;
           w = F.g;
           m = numel(w.g{1}.h);
           a = w.a(1);
           blocksize = a^J;
           r = (a^J-1)/(a-1)*(m-1);
           Lbrec = (floor(nextSb/blocksize) - floor(Sb/blocksize))*blocksize;
           rSb = (a^J-1)/(a-1)*(m-a) + mod(Sb,a^J);
           over = r - rSb;
           f = block_ifwt(c,w,J,Lbrec);
           ol = loadOverlap(r-mod(Sb, a^J),size(c,2),fola);
           olLen = size(ol,1);
           f(1:olLen-over,:) = f(1:olLen-over,:) + ol(1+over:end,:);
           f = [ol(1:over,:);f];
           if nargout>1
              fola=storeOverlap(f,r-mod(nextSb, a^J));
           else
              storeOverlap(f,r-mod(nextSb, a^J));
           end
           fhat = f(1:Lb,:);
          case {'dgt','dgtreal','dwilt','wmdct'}
           % Time step
           a = F.a;
           % Length of the left half of the window
           Lwl = floor(Lw/2);

           Sbonelmax =  ceil((Lw-1)/a)*a + a-1;
           Sbolen = ceil((Lw-1)/a)*a + mod(Sb,a);
           % Next block overlap length
           nextSbolen = ceil((Lw-1)/a)*a + mod(nextSb,a);
           Lext = Sbolen + Lb - mod(nextSb,a);
           Lextc = Sbolen + Lb - nextSbolen + Lwl;

           startc = ceil(Lwl/a)+1;
           endc = ceil((Lextc)/a);

           cc = F.coef2native(c,size(c));
           chat = zeros(size(cc,1),ceil(F.L/a),size(cc,3),class(cc));
           chat(:,startc:endc,:) = cc;
           chat = F.native2coef(chat);
           f = F.frsyn(chat);
           f = f(1:Lext,:);
           over = Sbonelmax - Sbolen;

           ol = loadOverlap(Sbonelmax-mod(Sb,a),size(c,2),fola);
           olLen = size(ol,1);
           f(1:olLen-over,:) = f(1:olLen-over,:) + ol(1+over:end,:);
           f = [ol(1:over,:);f];
           if nargout>1
              fola=storeOverlap(f,Sbonelmax-mod(nextSb,a));
           else
              storeOverlap(f,Sbonelmax-mod(nextSb,a));
           end
           fhat = f(1:Lb,:);
          case {'filterbank','filterbankreal'}
           lcma = F.lcma;
           % Time step
           a = F.a(:,1);
           % Length of the left half of the window
           Lwl = max(-F.g_info.offset);

                if Lw-1 < a
           Sbonelmax =   lcma-1;
           Sbolen =  mod(Sb,lcma);
           nextSbolen = mod(nextSb,lcma);
                else
           Sbonelmax =  ceil((Lw-1)/lcma)*lcma + lcma-1;
           Sbolen = ceil((Lw-1)/lcma)*lcma + mod(Sb,lcma);
           nextSbolen = ceil((Lw-1)/lcma)*lcma + mod(nextSb,lcma);
                end

           Lext = Sbolen + Lb - mod(nextSb,lcma);
           Lextc = Sbolen + Lb - nextSbolen + Lwl;

           startc = ceil(Lwl./a)+1;
           endc = ceil((Lextc)./a);

           cc = F.coef2native(c,size(c));

           chat = cell(numel(cc),1);
           for ii=1:numel(cc)
              chat{ii} = zeros(ceil(F.L./a(ii)),size(cc{ii},2));
              chat{ii}(startc(ii):endc(ii),:) = cc{ii};
           end

           chat = F.native2coef(chat);
           f = F.frsyn(chat);
           f = f(1:Lext,:);
           over = Sbonelmax - Sbolen;

           ol = loadOverlap(Sbonelmax-mod(Sb,lcma),size(c,2),fola);
           olLen = size(ol,1);
           f(1:olLen-over,:) = f(1:olLen-over,:) + ol(1+over:end,:);
           f = [ol(1:over,:);f];
           if nargout>1
              fola=storeOverlap(f,Sbonelmax-mod(nextSb,lcma));
           else
              storeOverlap(f,Sbonelmax-mod(nextSb,lcma));
           end
           fhat = f(1:Lb,:);
          otherwise
           error('%s: Unsupported frame.',upper(mfilename));
        end

    else
       error('%s: Frame was not created with blockaccel.',upper(mfilename));
    end

end

function overlap = loadOverlap(L,chan,overlap)
%LOADOVERLAP Loads overlap
%
%
   if isempty(overlap)
      overlap = block_interface('getSynOverlap');
   end

   if isempty(overlap)
     overlap = zeros(L,chan,block_interface('getClassId'));
   end
   Lo = size(overlap,1);
   if nargin<1
      L = Lo;
   end
   if L>Lo
      % Required more samples than stored
      if 0
         error('%s: Required more samples than stored.',upper(mfilename));
      else
         % pad with zeros
         overlapTmp = zeros(L,chan,block_interface('getClassId'));
         overlapTmp(end-Lo+1:end,:) = overlap;
         overlap = overlapTmp;
      end
   end
   overlap = overlap(end-L+1:end,:);
end

function overlap = storeOverlap(fext,L)
%STOREOVERLAP Stores overlap
%
%
    if L>size(fext,1)
        error('%s: Storing more samples than passed.',upper(mfilename));
    end
    overlap = fext(end-L+1:end,:);

    if nargout<1
       block_interface('setSynOverlap',overlap);
    end
end % STOREOVERLAP

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function demoImOrientedBox(varargin)
%DEMOIMORIENTEDBOX Demo file for using function imOrientedBox
%
%   Syntax
%   demoImOrientedBox
%   The demo runs automatically. A basic segmentation is performed on a
%   demo image, then oriented boxes are computed for each particle and
%   displayed on the original image.
%
%   Example
%   demoImOrientedBox
%
%   See also
%   demoImMaxFeretDiameter, imOrientedBox

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-09,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Image segmentation

% read image
img = imread('rice.png');
figure(1); clf;
imshow(img);

% compute background
bg = imopen(img, ones(30, 30));
img2 = img - bg;
figure(2); clf;
imshow(img2);

% display histogram, to identify the threshold value
figure(3); clf;
imHistogram(img2);

% image binarisation, and remove particles touching border
bin = img2 > 50;
bin = imclearborder(bin, 4);
imshow(bin);

% compute image labels, using minimal connectivity
lbl = bwlabel(bin, 4);
nLabels = max(lbl(:));

% display label image
rgb = label2rgb(lbl, jet(nLabels), 'w', 'shuffle');
figure(4); clf;
imshow(rgb);

%% Compute enclosing oriented boxes

% call the function
boxes = imOrientedBox(lbl);

% display result
hold on;
drawOrientedBox(boxes, 'linewidth', 2);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function results = vl_test_imwbackward(varargin)
% VL_TEST_IMWBACKWARD
vl_test_init ;

function s = setup()
s.I = im2double(imread(fullfile(vl_root,'data','spots.jpg'))) ;

function test_identity(s)
xr = 1:size(s.I,2) ;
yr = 1:size(s.I,1) ;
[x,y] = meshgrid(xr,yr) ;
vl_assert_almost_equal(s.I, vl_imwbackward(xr,yr,s.I,x,y)) ;

function test_invalid_args(s)
xr = 1:size(s.I,2) ;
yr = 1:size(s.I,1) ;
[x,y] = meshgrid(xr,yr) ;
vl_assert_exception(@() vl_imwbackward(xr,yr,single(s.I),x,y), 'vl:invalidArgument') ;

author: Cloud-CV

repo: object-proposals


text:
classdef CorrelateCross < dagnn.Filter
  properties
    opts = {'cuDNN'}
    ImageArea = [];
    BatchSize = [];
  end

  methods
    function outputs = forward(obj, inputs, params)
      % Assume the same size
      obj.BatchSize = size(inputs{1});
      obj.ImageArea =  obj.BatchSize(1) * obj.BatchSize(2);
      InputsAll = gpuArray.zeros(obj.ImageArea * obj.BatchSize(end), 1, obj.BatchSize(3), 1, 'single');
      for i = 1:obj.BatchSize(end)
          Start = (i - 1) * obj.ImageArea + 1;
          End =  i * obj.ImageArea;
          InputsAll(Start:End,:,:) = reshape(inputs{1}(:,:,:,i), [obj.ImageArea  1 obj.BatchSize(3)]);
      end
      outputs = obj.CorrelateFoward({InputsAll});
    end

    function outputs = CorrelateFoward(obj, inputs)
        [h,w,d,batchsize]= size(inputs{1});
        inputsfilter = reshape(permute(inputs{1},[3,1,2,4]),1,1,d,h*w,[]);
        outputs{1} = gpuArray.zeros(h,w,h*w,batchsize,'single');
        for i=1:batchsize
            outputs{1}(:,:,:,i) = vl_nnconv(inputs{1}(:,:,:,i), inputsfilter(:,:,:,:,i),[]) ;
        end
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      InputsAll = gpuArray.zeros(obj.ImageArea * obj.BatchSize(end), 1, obj.BatchSize(3), 1, 'single');
      for i = 1:obj.BatchSize(end)
          Start = (i - 1) * obj.ImageArea + 1;
          End =  i * obj.ImageArea;
          InputsAll(Start:End,:,:) = reshape(inputs{1}(:,:,:,i), [obj.ImageArea  1 obj.BatchSize(3)]);
      end
      derInputsAll = obj.CorrelateBackward({InputsAll}, derOutputs);
      derInputs{1} = gpuArray.zeros(obj.BatchSize, 'single');
      for i = 1:obj.BatchSize(end)
          Start = (i - 1) * obj.ImageArea + 1;
          End =  i * obj.ImageArea;
          derInputs{1}(:,:,:,i) = reshape(derInputsAll(Start:End,:,:), obj.BatchSize(1:3));
      end

      derParams = {} ;
    end

    function derInputs = CorrelateBackward(obj, inputs, derOutputs)
        [h,w,d,batchsize]= size(inputs{1});
        inputsfilter = reshape(permute(inputs{1},[3,1,2,4]),1,1,d,h*w,[]);
        derInputs1 = gpuArray.zeros(h,w,d,batchsize,'single');
        derInputs2 = gpuArray.zeros(h,w,d,batchsize,'single');

        for i=1:batchsize
            [dI1, dI2, ~] = vl_nnconv(...
                inputs{1}(:,:,:,i), inputsfilter(:,:,:,:,i), [], derOutputs{1}(:,:,:,i)) ;
            derInputs1(:,:,:,i) = dI1;
            derInputs2(:,:,:,i) = permute(reshape(dI2,d,h,w,[]),[2 3 1 4]);
        end
        derInputs = derInputs1 + derInputs2;
    end

    function obj = CorrelateCross(varargin)
      obj.load(varargin) ;
    end
  end
end

author: KuangJuiHsu


text:
m = find(tm(:,22)<0);
avr = sum(tm(:,22))/size(m,1);
var1 = (tm(:,22)-avr).^2;
var = sum(var1)/size(m,1);
stdev = var.^(1/2);

author: CelsoReyes


text:
function OUT = Date2Num(DATA)
% =======================================================================
% Substitute cell dates 1989Q4 with numeric dates 1989.75. Only for
% quarterly data.
% =======================================================================
% OUT = Date2Num(DATA)
% -----------------------------------------------------------------------
% INPUTS
%	- DATA: a (T x 1) vector with dates in cell format
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT:  a (T x 1) vector with dates in numeric format
% =========================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogiocesabianchi@gmail.com
%-------------------------------------------------------------------------

OUT = nan(size(DATA,1),size(DATA,2));
for i=1:size(DATA,1)
    for j=1:size(DATA,2)
        aux = DATA{i,j};
        year = str2double(aux(1:4));
        quarter = str2double(aux(6));
        OUT(i,j) = year + quarter/4 -0.25;
    end
end

author: ambropo

repo: VAR-Toolbox

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


text:
function Smooth = my_conv_local(S1, sig)
NN = size(S1,1);
NT = size(S1,2);
dt = -4*sig:1:4*sig;
gaus = exp( - dt.^2/(2*sig^2));
gaus = gaus'/sum(gaus);
Smooth = filter(gaus, 1, [S1' ones(NT,1); zeros(4*sig, NN+1)]);
Smooth = Smooth(1+4*sig:end, :);
Smooth = Smooth(:,1:NN) ./ (Smooth(:, NN+1) * ones(1,NN));
Smooth = Smooth';

author: cortex-lab


text:
%NEWJOR Create an Jordan backpropagation network.
%
%  Syntax
%
%    net = newjor(P,T,[S1...S(N-l)],{TF1...TFN},BTF,BLF,PF,IPF,OPF,DDF)
%
%  Description
%
%     NET = NEWJOR(P,T,[S1...SNl],{TF1...TFN},BTF,BLF,PF,IPF,OPF,DDF) takes,
%      P  - RxQ1 matrix of Q1 representative R-element input vectors.
%      T  - SNxQ2 matrix of Q2 representative SN-element target vectors.
%      Si  - Sizes of N-1 hidden layers, S1 to S(N-1), default = [].
%            (Output layer size SN is determined from T.)
%      TFi - Transfer function of ith layer. Default is 'tansig' for
%            hidden layers, and 'purelin' for output layer.
%      BTF - Backprop network training function, default = 'traingdx'.
%      BLF - Backprop weight/bias learning function, default = 'learngdm'.
%      PF  - Performance function, default = 'mse'.
%      IPF - Row cell array of input processing functions.
%            Default is {'fixunknowns','removeconstantrows','mapminmax'}.
%      OPF - Row cell array of output processing functions.
%            Default is {'removeconstantrows','mapminmax'}.
%      DDF - Data division function, default = 'dividerand';
%    and returns an N-layer Jordan network.
%
%    The training function BTF can be any of the backprop training
%    functions such as TRAINGD, TRAINGDM, TRAINGDA, TRAINGDX, etc.
%
%    *WARNING*: Algorithms which take large step sizes, such as TRAINLM,
%    and TRAINRP, etc., are not recommended for Jordan networks.  Because
%    of the delays in Jordan networks the gradient of performance used
%    by these algorithms is only approximated making learning difficult
%    for large step algorithms.
%
%    The learning function BLF can be either of the backpropagation
%    learning functions such as LEARNGD, or LEARNGDM.
%
%    The performance function can be any of the differentiable performance
%    functions such as MSE or MSEREG.
% 	Mark Beale, 11-31-97
% 	Copyright 1992-2008 The MathWorks, Inc.
% 	$Revision: 1.1.6.9 $  $Date: 2009/08/08 01:11:34 $

function net = newjor(varargin)

% Arnav Goel, 11-06-12
% Copyright 2012 BITS Pilani KK Birla Goa Campus
% $Revision 0.0.3 $Date 21-11-12

%This portion of the code refers to error declaration of the code inorder
%to make sure that the no of the arguments and the nature of the arguments
%is compatible to the NNT ver > 5.0 for the jordan_p1

if nargin < 1, error('NNET:Arguments','Not enough input arguments'),end

v1 = varargin{1};
if isa(v1,'cell'), v1 = cell2mat(v1); end
v2 = varargin{2};
if nargin > 2, v3 = varargin{3}; end

if (nargin<= 6) && (size(v1,2)==2) && (~iscell(v2)) && (size(v2,1)==1) && ((nargin<3)||iscell(v3))
  nntobsu(mfilename,['See help for ' upper(mfilename) ' to update calls to the new argument list.']);
  net = jordan_p0(varargin{:});
else
  net = jordan_p1(varargin{:});
end

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

%this portion of the code is the main part with designing the Architecture,
%Simulation, Adaptation, Training & Initialization

function net = jordan_p1(p,t,s,tf,btf,blf,pf,ipf,tpf,ddf)
if nargin < 1, error('NNET:Arguments','Not enough arguments.'),end

% Defaults
if (nargin < 3), s = []; end
if (nargin < 4), tf = {}; end
if (nargin < 5), btf = 'traingdx'; end
if (nargin < 6), blf = 'learngdm'; end
if (nargin < 7), pf = 'mse'; end
if (nargin < 8), ipf = {'fixunknowns','removeconstantrows','mapminmax'}; end
if (nargin < 9), tpf = {'removeconstantrows','mapminmax'}; end
if (nargin < 10), ddf = 'dividerand'; end

% Format
if isa(p,'cell'), p = cell2mat(p); end
if isa(t,'cell'), t = cell2mat(t); end

% Error checking
if (~isa(p,'double')) || ~isreal(p)
  error('NNET:Arguments','Inputs are not a matrix or cell array with a single matrix.')
end
if isa(s,'cell')
  if (size(s,1) ~= 1)
    error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
  end
  for i=1:length(s)
    si = s{i};
    if ~isa(si,'double') || ~isreal(si) || any(size(si) ~= 1) || any(si<1) || any(round(si) ~= si)
      error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
    end
  end
  s = cell2mat(s);
end
if (~isa(s,'double')) || ~isreal(s) || (size(s,1) ~= 1) || any(s<1) || any(round(s) ~= s)
  error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
end

% Architecture

Nl = length(s)+1; %N-1 hidden + 1 output to decide the no of layers for the network
net = network(1,Nl); % 1 input N-1 HL 1 output layer network
net.biasConnect = ones(Nl,1);
net.inputConnect(1,1) = 1;
[j,i] = meshgrid(1:Nl,1:Nl);
net.layerConnect = (j == (i-1)) | ((j-1) == i); % here LC(1,2) = 1 and LC(2,1) = 1
net.outputConnect(Nl) = 1;

% Simulation
net.inputs{1}.exampleInput = p;
net.inputs{1}.processFcns = ipf; %pre-processing conditions for less arguments in the newjor

for i=1:Nl
    if (i<Nl)
        net.layerWeights{i,i+1}.delays = 1;
    else
       net.layerWeights{Nl,Nl}.delays = 1 ;
    end

  if (i<Nl), net.layers{i}.size = s(i); end
  if (length(tf) < i) || all(isnan(tf{i}))
    if (i<Nl)
      net.layers{i}.transferFcn = 'tansig';
    else
      net.layers{i}.transferFcn = 'purelin';
    end
  else
    net.layers{i}.transferFcn = tf{i};
  end
end
net.outputs{Nl}.exampleOutput = t;
net.outputs{Nl}.processFcns = tpf;

% Adaption
net.adaptfcn = 'trains';
net.inputWeights{1,1}.learnFcn = blf;
for i=1:Nl
  net.biases{i}.learnFcn = blf;
  net.layerWeights{i,:}.learnFcn = blf;
end

% Training
net.trainfcn = btf;
net.dividefcn = ddf;
net.performFcn = pf;
net.gradientFcn = feval(btf,'gdefaults',0);

% Initialization
net.initFcn = 'initlay';
for i=1:Nl
  net.layers{i}.initFcn = 'initnw';
end
net = init(net);

% Plots
net.plotFcns = {'plotperform','plottrainstate'};

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

%this is similar to the above function written for the lower end MATLAB
%version of NNT tooblox

function net = jordan_p0(pr,s,tf,btf,blf,pf)
% Backward compatible to NNT 5.0

if nargin < 1, error('NNET:Arguments','Not enough arguments.'),end

% Defaults
Nl = length(s);
if nargin < 3, tf = {'tansig'}; tf = tf(ones(1,Nl)); end
if nargin < 4, btf = 'traingdx'; end
if nargin < 5, blf = 'learngdm'; end
if nargin < 6, pf = 'mse'; end

% Error checking
if (~isa(pr,'double')) || ~isreal(pr) || (size(pr,2) ~= 2)
  error('NNET:Arguments','Input ranges is not a two column matrix.')
end
if any(pr(:,1) > pr(:,2))
  error('NNET:Arguments','Input ranges has values in the second column larger in the values in the same row of the first column.')
end
if isa(s,'cell')
  if (size(s,1) ~= 1)
    error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
  end
  for i=1:length(s)
    si = s{i};
    if ~isa(si,'double') || ~isreal(si) || any(size(si) ~= 1) || any(si<1) || any(round(si) ~= si)
      error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
    end
  end
  s = cell2mat(s);
end
if (~isa(s,'double')) || ~isreal(s) || (size(s,1) ~= 1) || any(s<1) || any(round(s) ~= s)
  error('NNET:Arguments','Layer sizes is not a row vector of positive integers.')
end

% Architecture
net = network(1,Nl);
net.biasConnect = ones(Nl,1);
net.inputConnect(1,1) = 1;
[j,i] = meshgrid(1:Nl,1:Nl);
net.layerConnect = (j == (i-1)) | ((j-1) == i);
net.outputConnect(Nl) = 1;

% Simulation
net.inputs{1}.range = pr;
for i=1:Nl
  if (i<Nl)
        net.layerWeights{i,i+1}.delays = [1];
    else
       net.layerWeights{Nl,Nl}.delays = [1] ;
    end
  net.layers{i}.size = s(i);
  net.layers{i}.transferFcn = tf{i};
end

% Performance
net.performFcn = pf;

% Adaption
net.adaptfcn = 'trains';
net.inputWeights{1,1}.learnFcn = blf;
for i=1:Nl
  net.biases{i}.learnFcn = blf;
  net.layerWeights{i,:}.learnFcn = blf;
end

% Training
net.trainfcn = btf;
net.gradientFcn = feval(btf,'gdefaults',0);

% Initialization
net.initFcn = 'initlay';
for i=1:Nl
  net.layers{i}.initFcn = 'initnw';
end
net = init(net);

% Plots
net.plotFcns = {'plotperform','plottrainstate'};

author: Sable

repo: mcbench-benchmarks


text:
function result = make_epochs_per_sample(weights)
%MAKE_EPOCHS_PER_SAMPLE Given a set of weights generate the number of
% epochs per sample for each weight.
%
% result = MAKE_EPOCHS_PER_SAMPLE(weights)
%
% Parameters
% ----------
% weights: array of size (n_1_simplices, 1)
%     The weights of how much we wish to sample each 1-simplex.
%
% Note that the total number of epochs does not impact this result.
%
% Returns
% -------
% result: array of size (n_1_simplices, 1)
%     The number of epochs per sample, one for each 1-simplex.
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

result = -1*ones(size(weights, 1), 1);
L = weights > 0;
result(L) = max(weights)*ones(sum(L), 1)./weights(L);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function nvar = p12_nvar ( option )

%*****************************************************************************80
%
%% P12_NVAR sets the number of variables for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, integer NVAR, the number of variables.
%
  if ( option == 1 )
    npolys = 2;
    nderiv = 1;
  elseif ( option == 2 )
    npolys = 4;
    nderiv = 1;
  elseif ( option == 3 )
    npolys = 4;
    nderiv = 2;
  elseif ( option == 4 )
    npolys = 6;
    nderiv = 1;
  elseif ( option == 5 )
    npolys = 6;
    nderiv = 2;
  elseif ( option == 6 )
    npolys = 6;
    nderiv = 3;
  end

  nint = 8;
  nvary = nint * npolys;
  nbcz = 1;
  nbco = 1;
  nvarz = nbcz + ( nint - 1 ) * nderiv + nbco;
  nvar = nvary + nvarz + 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function structout = setstructfields(structin, newfields)
%SETSTRUCTFIELDS Set fields of a structure using another structure
%   SETSTRUCTFIELDS(STRUCTIN, NEWFIELDS) Set fields of STRUCTIN using
%   another structure NEWFIELDS fields.  If fields exist in STRUCTIN
%   but not in NEWFIELDS, they will not be changed.

%   Author(s): J. Schickler
%   Copyright 1988-2004 The MathWorks, Inc.
%   $Revision: 1.3.4.3 $  $Date: 2007/12/14 15:16:12 $
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%narginchk(2,2);

% Set the output to the input
structout = structin;

% If newfields is not a structure return
if ~isstruct(newfields), return; end

% Loop over all the fields of newfields and assign them to the output
fields = fieldnames(newfields);
for i = 1:length(fields),
    value = newfields.(fields{i});
    if isstruct(value) && isfield(structout, fields{i}) && isstruct(structout.(fields{i}))
        structout.(fields{i}) = setstructfields(structout.(fields{i}), value);
    else
        structout.(fields{i}) = newfields.(fields{i});
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [um] = A2um(A)
% Convert length from angstroms to micrometers (or microns).
% Chad A. Greene 2012
um = A*.0001;


text:
function [U,S,output] = lmlra_rnd(size_tens,size_core,options)
%LMLRA_RND Pseudorandom initialization for LMLRA.
%   [U,S] = lmlra_rnd(size_tens,size_core) generates pseudorandom
%   unitary factor matrices U{1}, ..., U{N} of dimensions size_tens(n)-by-
%   size_core(n) and a core tensor of dimensions size_core that can be used
%   to initialize algorithms that compute a low multilinear rank
%   approximation of an N-th order tensor.
%
%   lmlra_rnd(T,size_core) is shorthand for lmlra_rnd(size(T),size_core) if
%   T is real. If T is complex, then by default S and U{n} will be
%   generated as a complex tensor and matrices, respectively (cf. options).
%
%   lmlra_rnd(size_tens,size_core,options) and lmlra_rnd(T,size_core, ...
%   options) may be used to set the following options:
%
%      options.Real =         - The type of random number generator used to
%      [{@randn}|@rand|0]       generate the real part of the core tensor S
%                               and matrices U{n}. If 0, there is no real
%                               part.
%      options.Imag =         - The type of random number generator used to
%      [@randn|@rand|0|...      generate the imaginary part of the core
%       {'auto'}]               tensor S and matrices U{n}. If 0, there is
%                               no imaginary part. On 'auto', options.Imag
%                               is 0 unless the first argument is a complex
%                               tensor T, in which case it is equal to
%                               options.Real.
%      options.Unitary = true - Set equal to false to generate the factor
%                               matrices U{n} without converting them to
%                               unitary matrices afterwards.
%
%   See also btd_rnd, cpd_rnd, lmlragen.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

% Process input.
isSizeVector = any(size(size_tens) == numel(size_tens));
if ~isSizeVector, T = size_tens; size_tens = size(size_tens); end
N = length(size_tens);
if N ~= length(size_core)
    error('lmlra_rnd:size_core', ...
         ['length(size_core) should be equal to ndims(T) or ' ...
          'length(size_tens).']);
end

% Check the options structure.
isfunc = @(f)isa(f,'function_handle');
if nargin < 3, options = struct; end
if ~isfield(options,'Real'), options.Real = @randn; end
if ~isfunc(options.Real), options.Real = @zeros; end
if ~isfield(options,'Imag'), options.Imag = 'auto'; end
if ~isfield(options,'Unitary'), options.Unitary = true; end
if ischar(options.Imag) && strcmpi(options.Imag,'auto')
    if ~isSizeVector && ~isreal(T)
        options.Imag = options.Real;
    else
        options.Imag = 0;
    end
end

% Generate factor matrices and core tensor.
U = arrayfun(@(n)options.Real(size_tens(n),size_core(n)),1:N, ...
        'UniformOutput',0);
if isfunc(options.Imag)
    Ui = arrayfun(@(n)options.Imag(size_tens(n),size_core(n)),1:N, ...
            'UniformOutput',0);
    U = cellfun(@(ur,ui)ur+ui*1i,U,Ui,'UniformOutput',0);
end
if options.Unitary
    for n = 1:N
        if size(U{n},1) >= size(U{n},2), [U{n},~] = qr(U{n},0);
        else [Q,~] = qr(U{n}.',0); U{n} = Q.'; end
    end
end
S = options.Real(size_core(:).');
if isfunc(options.Imag), S = S+1i*options.Imag(size_core(:).'); end
output = struct;

author: andrewssobral

repo: lrslibrary


text:
classdef ZDT5 < PROBLEM
% <multi> <binary> <large/none> <expensive/none>
% Benchmark MOP proposed by Zitzler, Deb, and Thiele

%------------------------------- Reference --------------------------------
% E. Zitzler, K. Deb, and L. Thiele, Comparison of multiobjective
% evolutionary algorithms: Empirical results, Evolutionary computation,
% 2000, 8(2): 173-195.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 80; end
            obj.D        = ceil(max(obj.D-30,1)/5)*5 + 30;
            obj.encoding = 4 + zeros(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            u      = zeros(size(PopDec,1),1+(size(PopDec,2)-30)/5);
            u(:,1) = sum(PopDec(:,1:30),2);
            for i = 2 : size(u,2)
                u(:,i) = sum(PopDec(:,(i-2)*5+31:(i-2)*5+35),2);
            end
            v           = zeros(size(u));
            v(u<5)      = 2 + u(u<5);
            v(u==5)     = 1;
            PopObj(:,1) = 1 + u(:,1);
            g           = sum(v(:,2:end),2);
            h           = 1./PopObj(:,1);
            PopObj(:,2) = g.*h;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = 1 : 31;
            R(:,2) = (obj.D-30)./5./R(:,1);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function res = vg_method(args)

% Variational Garrote: performs linear regression with L0-norm penalty (Spike and Slab model)
%   See file test_vg.m for an example of use
%
% required parameters (n input dimension, p samples)
%   x     : n x p (training set, input)
%   y     : 1 x s (training set, output)
%   xv    : n x p (validation set, input)
%   yv    : 1 x s (validation set, output)
%
% optional parameters (default)
%   method      : method for optimization 'dual' or 'regression' for fixed gamma ('dual')
%   maxiter     : maximum number of iterations for optimization for fixed gamma (1e4)
%   max_sum_m   : increases gamma values until sum(m)=max_sum_m  (n/2)
%   beta_max    : increases gamma values until beta=beta_max (1e3)
%   n_gamma     : number of gamma values to scan (50)
%   dmmin       : convergence threshold for mean field error (1e-12)

%----------------
% REQUIRED PARAMS
ok = true;
if ~isfield(args, 'x') disp('train input x not provided'); ok = false; else x=args.x; end
if ~isfield(args, 'y') disp('train output y not provided'); ok = false; else y=args.y; end
if ~isfield(args, 'xv') disp('val set input xv not provided'); ok = false; else xv=args.xv; end
if ~isfield(args, 'yv') disp('val set output yv not provided'); ok = false; else yv=args.yv; end
if isfield(args, 'xt') && isfield(args, 'yt')
    xt=args.xt;
    yt=args.yt;
    pt = size(xt,2);
else
    pt = 0;
end

if ~ok
    return;
end

n = size(x,1);
p = size(x,2);
pv = size(xv,2);

%----------------
% OPTIONAL PARAMS

% method for optimization {dual or regression} for fixed gamma
if ~isfield(args, 'method') method='dual'; else method=args.method; end

% maximum number of iterations for optimization for fixed gamma
if ~isfield(args, 'maxiter') maxiter=1e4; else maxiter=args.maxiter; end

% increases gamma values until sum(m)=max_sum_m
if ~isfield(args, 'max_sum_m') max_sum_m=n/2; else max_sum_m=args.max_sum_m; end

% increases gamma values until beta=beta_max
if ~isfield(args, 'beta_max') beta_max=1e3; else beta_max=args.beta_max; end

% number of gamma values to scan
if ~isfield(args, 'n_gamma') n_gamma=50; else n_gamma=args.n_gamma; end

% convergence threshold for mean field error
if ~isfield(args, 'dmmin') dmmin=1e-12; else dmmin=args.dmmin; end

%----------------
% compute garrote solution for range of gammas
% first from gamma_min to gamma_max and then in
% a second pass from gamma_max to gamma_min.

% C is input data covariance matrix.
if strcmp(method, 'regression')
    if n<=1500,
        C=x*x'/p;
    end;
end

% b is input output covariance
b=x*y'/p;

% sigma is output variance
sigmay=y*y'/p;

% set gamma range (min, max and step size)
delta=1e-8;
[b2sort,isort]=sort(b.^2,'descend');
bsort=b(isort);
gamma_min=log(delta*sigmay/p/max(abs(b)));
eps_gamma=0.001;
gamma_max=eps_gamma*gamma_min;
gamma_all =linspace(gamma_min,gamma_max,n_gamma);

% initial step size of mean field update
eta0=1; %e-2;
% initial step size for change in w in dual.m
eta_w0=0.02;

% input data variance
chi_ii=1/p*sum(x.^2,2);
if sum(abs(chi_ii-1)>1e-10),
    fprintf('input design matrix is not normalized\n');
    pause
end;

lg=n_gamma;
kl_all=inf(lg,2);
v_all=inf(lg,2,n);
m_all=inf(lg,2,n);
beta_all=inf(lg,2);
v_mf_all=inf(lg,n);
m_mf_all=inf(lg,n);
iter_all=inf(lg,2);
beta_mf_all=inf(1,lg);
error_mf_all=inf(1,lg);
errorv_mf_all=inf(1,lg);
errort_mf_all=inf(1,lg);
m=zeros(1,n);

% the estimated inverse noise variance beta is initialized as the
% output variance
beta=1/sigmay;
i=0;

% for gamma is gamma_min to gamma_max, or when some criteria are
% satisfied
while (beta<beta_max)&&(i<n_gamma)&&(sum(m)<max_sum_m),
	i=i+1;
	gamma=gamma_all(i);
    eval(method);
	v_all(i,1,:)=v;
	m_all(i,1,:)=m;
	beta_all(i,1)=beta;
	iter_all(i,1)=iter;
	kl_all(i,1)=kl1;
	fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
end;

if beta>=beta_max
 fprintf('-----------------------------------------------------------\n');
 fprintf('beta > beta_max (%.3f > %.3f)\n', beta, beta_max);
 if i>1
    m = squeeze(m_all(i-1,1,:))';
 end
end
if sum(m)>=max_sum_m
 fprintf('-----------------------------------------------------------\n');
 fprintf('sum(m) > max_sum_m (%.3f > %.3f)\n', sum(m), max_sum_m);
end

% for gamma is current gamma decreasing to gamma_min
imax=i-1;
for i=imax:-1:1,
	gamma=gamma_all(i);
    eval(method);
	v_all(i,2,:)=v;
	m_all(i,2,:)=m;
	beta_all(i,2)=beta;
	iter_all(i,2)=iter;
	kl_all(i,2)=kl1;
	fprintf('gamma = %f beta = %f sum(m) = %f iter = %d kl = %f\n',gamma,beta,sum(m),iter,kl1);
end;

% select for each gamma from these two solutions the one with lowest KL
[klmin, imin]=min(kl_all,[],2);
for i=1:imax,
	v_mf_all(i,:)=squeeze(v_all(i,imin(i),:))';
	m_mf_all(i,:)=squeeze(m_all(i,imin(i),:))';
	beta_mf_all(i)=beta_all(i,imin(i));
	error_mf_all(i)=1/p*sum((y-v_mf_all(i,:)*x).^2,2);
	errorv_mf_all(i)=1/pv*sum((yv-v_mf_all(i,:)*xv).^2,2);
	if pt>0,
		errort_mf_all(i)=1/pt*sum((yt-v_mf_all(i,:)*xt).^2,2);
	end;
end;

% select the gamma that optimizes the validation error (errorv_mf_all)
[minerrorv i]=min(errorv_mf_all(1:imax));

res.gamma_mf=gamma_all(i);
res.v_mf=v_mf_all(i,:);
res.m_mf=m_mf_all(i,:);
res.n_mf1=sum(res.m_mf>0.5);

res.error_mf=error_mf_all(i);
res.errorv_mf=errorv_mf_all(i);
res.errort_mf=errort_mf_all(i);
res.beta_mf=beta_mf_all(i);

if (res.beta_mf==beta_max)
	fprintf('beta_max too small: beta_mf %6.4f, beta_max %6.4f\n', beta_mf(iruns),beta_max);
	pause
end;
if (res.gamma_mf==gamma_min)
	fprintf('gamma at minimum range boundary\n');
end;
if (res.gamma_mf==gamma_max)
	fprintf('gamma at maxium range boundary\n');
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef polyMesh
% POLYMESH Class for 3d polygon meshes.

    properties (SetAccess = public, GetAccess = public)
        % Vertices
        vertices

        % Faces
        faces
    end

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

    methods

        function obj = polyMesh(meshData)
        % POLYMESH Import of mesh data.
        % ----------------------------------------------------------------------
        % INPUT
        % [meshData]
        % Mesh input data as ply file. Currently no other formats are
        % supported.
        % ----------------------------------------------------------------------
        % OUTPUT
        % [obj]
        % Object instance of class polyMesh.
        % ----------------------------------------------------------------------
        % philipp.glira@gmail.com
        % ----------------------------------------------------------------------

        % Input parsing --------------------------------------------------------

        p = inputParser;
        p.addRequired('meshData', @(x) ischar(x));
        p.parse(meshData);
        p = p.Results;
        % Clear required inputs to avoid confusion
        clear meshData

        % Check if input file exists -------------------------------------------

        if ischar(p.meshData)
            if exist(p.meshData) ~= 2
                error('File ''%s'' does not exist!', p.meshData);
            end
        end

        % Start ----------------------------------------------------------------

        procHierarchy = {'POLYMESH' 'IMPORT'};
        msg('S', procHierarchy);

        % Import of mesh -------------------------------------------------------

        data = plyread(p.meshData);

        obj.vertices = [data.vertex.x data.vertex.y data.vertex.z];
        obj.faces    = vertcat(data.face.vertex_indices{:}) + 1; % increase IDs by one!

        % End ------------------------------------------------------------------

        procHierarchy = {'POLYMESH' 'IMPORT'};
        msg('E', procHierarchy);

        end

    end

end

author: pglira


text:
function id = generatejobid(batch, batchid)

% GENERATEJOBID generates a unique identifier for a job to be submitted to the
% batch queueing system. It maintains an internal counter to allow it to be
% called from multiple qsubfeval instances without the user having to keep
% track of the numbers.
%
% Use as
%   jobid     = generatejobid(batch)
%   batchid   = generatebatchid(batch)
%   sessionid = generatesessionid()
%
% The result is a string like
%   user_host_pid_bM_jN  % as jobid
%   user_host_pid_bM     % as batchid
%   user_host_pid        % as sessionid
% where M is the batch number and N the sequential job number (per batch).
%
% Besides specifying a batch number, it is also possible to specify the full
% batch name as string. This allows the user to override the default
% user_host_pid_bM part of the jobid. This works like
%   jobid     = generatejobid(batch, batchid)
% where for example generatejobid(1, 'freqanalysis') returns the sequence
% of strings 'freqanalysis_j001', 'freqanalysis_j002', ...
%
% See also GENERATEBATCHID, GENERATESESSIONID

% Copyright (C) 2011-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% jobNum will be a vector with one entry for each batch number
persistent jobNum

if nargin<1
  error('incorrect number of input arguments');
end

if nargin<2
  batchid = [];
end

if length(jobNum)>=batch
  % increment the previous value
  job = jobNum(batch)+1;
else
  % start with job number one
  job = 1;
end

if ~isempty(batchid)
  id = sprintf('%s_j%03d', batchid, job);
else
  id = sprintf('%s_%s_p%d_b%d_j%03d', getusername(), gethostname(), getpid(), batch, job);
end

% ensure that it can be used as filename, struct fieldname, etc.
id = fixname(id);

% remember the current job number for the current batch
jobNum(batch) = job;


text:
function disparity_in_pixels = disparity_in_pixels_cityscapes(input_disparity)

disparity_in_pixels = (double(input_disparity) - 1) / 256;

end

author: sakaridis


text:
clear;
%
addpath '..';
%
fs=100;
nt=1000;
t=[1:nt]./fs;t=t(:);
%
f0=1.5;f1=10/2.15;
opm=log2(f1/f0)*60/t(end);
fup=f0*2.^(opm/60*t);fup=fup(:);
%
f0=10-(10-1.5)/2;f1=10;
inc=(f1-f0)/nt;
fdn=[f1:-inc:f0+inc];fdn=fdn(:);
%
f=0.75/t(end);
am1=0.25*sin(2*pi*f*t-pi/4);am1=am1+0.75;
am2=0.25*sin(2*pi*f*t-pi/1.7);am2=am2+0.75;
am3=0.25*sin(2*pi*f*t-pi/8);am3=am3+0.75;
am4=0.25*sin(2*pi*f*t-pi/2.5);am4=am4+0.75;
am5=0.25*sin(2*pi*f*t-pi/5);am5=am5+0.75;
am6=0.25*sin(2*pi*f*t-pi/10);am6=am6+0.75;
%
x1=am1.*sin(1*2*pi*cumsum(fup/fs));x1=x1./std(x1);
x2=am2.*sin(2*2*pi*cumsum(fup/fs));x2=x2./std(x2);
x3=am3.*sin(3*2*pi*cumsum(fup/fs));x3=x3./std(x3);
x4=am4.*sin(1*2*pi*cumsum(fdn/fs));x4=x4./std(x4);
x5=am5.*sin(2*2*pi*cumsum(fdn/fs));x5=x5./std(x5);
x6=am6.*sin(3*2*pi*cumsum(fdn/fs));x6=x6./std(x6);
%
fp=[fdn,2*fup,fup];
xp=[x4,x2,x1];
[nt,nch]=size(xp);
%
nsens=1;
%A=randn(nsens,nch);
A=[1,1,1];
%
xe=(A*xp.').';
xc=A(ones(nt,1),:).*xp;
%
%xe=xe+0.05*randn(nt,nsens); % add noise
%
nfft=128;nwin=100;novlap=50;dflag='mean';
[Sxx_m,freq,time]=p_gram(xe,nfft,fs,nwin,novlap,dflag);
%
r=4000;
ford=1;
[xs1,bw,T] = vk2(xe,fp(:,1),fs,r,ford);
[xs2,bw,T] = vk2(xe,fp(:,2),fs,r,ford);
[xs3,bw,T] = vk2(xe,fp(:,3),fs,r,ford);
xs=[xs1,xs2,xs3];
%
figure(1);
h=plot(t,xc,'--',t,abs(xs),'r');
set(h(2:end),'linewidth',2);
xlabel('Time [Sec]');
ylabel('Amplitude');
title('Single Order Extraction');
%
rm=4000;
tol=0.001;maxit=1000;
[xm,bwm,Tm,fl,rr,it,rv] = vkm(xe,fp,fs,rm*ones(nch,1),ford,tol,maxit);
%
figure(2);
h=plot(t,xc,'--',t,abs(xm),'r');
set(h(2:end),'linewidth',2);
xlabel('Time [Sec]');
ylabel('Amplitude');
title('Multiple Order Extraction');
%
figure(3);
semilogy(rv);
xlabel('Iteration');
ylabel('Residual');
%
figure(4);
typ='lin';f_min=0;f_max=32.5;
spec_plot2(Sxx_m,freq,time,typ,f_min,f_max,0,0.8);
colorbar('off');
g=colormap('gray');colormap(flipud(g));
title('Spectrogram of orders');
%
%%


text:
function test_suite = test_graphRadius
%TEST_GRAPHRADIUS  One-line description here, please.
%
%   output = test_graphRadius(input)
%
%   Example
%   testGraphRadius
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-05-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testGraph02(testCase) %#ok<*DEFNU>

[nodes, edges] = createTestGraph02;
r = graphRadius(nodes, edges);

exp = 3;
testCase.assertEqual(exp, r);

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function complex_example_AD()
% A basic example that shows how to define the cost funtion for
% optimization problems on complex manifolds.
%
% Note that automatic differentiation for complex numbers is not supported
% for Matlab R2021a or earlier. To fully exploit the convenience of AD,
% please update to the latest version if possible. If the user cannot have
% access to Matlab R2021b or later, manopt provides an alternative way to
% deal with complex problems which requires the user to define the cost
% funtion using the basic functions listed in the folder /functions_AD or
% to define their own functions following the rules described in that file.
% See the following as an example.
%
% See also: manoptADhelp

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Xiaowen Jiang, August, 31, 2021
% Contributors: Nicolas Boumal
% Change log:
%

    % Verify that the deep learning tool box was installed
    assert(exist('dlarray', 'file') == 2, ['Deep learning tool box is '...
    'needed for automatic differentiation.\n Please install the'...
    'latest version of the deep learning tool box and \nupgrade to Matlab'...
    ' R2021b if possible.'])

    % Generate the problem data.
    n = 100;
    A = randn(n, n) + 1i*randn(n, n);
    A = .5*(A+A');

    % Create the problem structure.
    S = spherecomplexfactory(n);
    problem.M = S;

    % Define the problem cost function
    % For Matlab R2021b or later, define the problem cost function as usual
    % problem.cost  = @(X) -.5*real(X'*A*X);

    % For Matlab R2021a or earlier, translate the cost function into a
    % particular format with the basic functions in /functions_AD
    problem.cost  = @(X) -creal(cprod(cprod(ctransp(X), A), X));

    % Define the gradient and the hessian via automatic differentiation
    problem = manoptAD(problem);

    % Numerically check gradient and Hessian consistency.
    figure;
    checkgradient(problem);
    figure;
    checkhessian(problem);

    % Solve.
    [x, xcost, info] = trustregions(problem);          %#ok<ASGLU>

    % Display some statistics.
    figure;
    semilogy([info.iter], [info.gradnorm], '.-');
    xlabel('Iteration #');
    ylabel('Gradient norm');
    title(['Convergence of the trust-regions algorithm on the'...
        'complex sphere power manifold']);

end

author: NicolasBoumal

repo: manopt


text:
function srf = readBESAsrf(filename)

% This method reads surface files from BESA Research/MRI. For more
% information on the srf file format also see the file format description
% provided by the BrainVoyager support site:
% http://support.brainvoyager.com
%
% Parameters:
%     [filename]
%         In the case that the current folder is not the folder containing
%         the file it should be the full path including the name of the
%         elp file else only the name of the file should be specified.
%
% Return:
%     srf
%         Cell array containing info on surface points stored in file.
%
% Copyright (C) 2013, BESA GmbH
%
% File name: readBESAsrf.m
%
% Author: Robert Spangler
% Created: 2013-11-26

% Check default parameters
% no default parameters so far...

% Ouput parameter
srf = [];

% Try to open file.
FileID = fopen(filename, 'r');
if(FileID < 0)
	printf('Error opening file.\n');
	return
end

% Read version number
srf.VerNo = fread(FileID, 1, 'float');

% Read reserved bytes, must be '0'
srf.Reserved = fread(FileID, 1, 'int');

% Read number of vertices and triangles
srf.NoVertices = fread(FileID, 1, 'int');
srf.NoTriangles = fread(FileID, 1, 'int');

% Read mesh center (X,Y,Z), default (128.0, 128.0, 128.0)
srf.MeshCenterCoord = fread(FileID, 3, 'float');

% Read sequence of (X,Y,Z) coordinates of all vertices
srf.CoordsVertices = [];
% 1. Read NrVertices X coordinates
srf.CoordsVertices(:,1) = fread(FileID, srf.NoVertices, 'float');
% 2. Read NrVertices Y coordinates
srf.CoordsVertices(:,2) = fread(FileID, srf.NoVertices, 'float');
% 3. Read NrVertices Z coordinates
srf.CoordsVertices(:,3) = fread(FileID, srf.NoVertices, 'float');

% Read normal (X,Y,Z) components of all vertices
% Note: For historical reasons, normals do not point outward, but inward.
%       Multiply each normal by "-1.0" for custom rendering.
srf.ComponentsVertices = [];
% 1. Read NrVertices X components
srf.ComponentsVertices(:,1) = fread(FileID, srf.NoVertices, 'float');
% 2. Read NrVertices Y components
srf.ComponentsVertices(:,2) = fread(FileID, srf.NoVertices, 'float');
% 3. Read NrVertices Z components
srf.ComponentsVertices(:,3) = fread(FileID, srf.NoVertices, 'float');

% Read convex curvature color (R,G,B,A), default (0.322, 0.733, 0.98, 1.0)
srf.ConvCurvCol = fread(FileID, 4, 'float');

% Read concave curvature color (R,G,B,A), default (0.1, 0.24, 0.32, 1.0)
srf.ConcCurvCol = fread(FileID, 4, 'float');

% Read mesh color (sequence of color indices, one for each vertex)
srf.MeshCol = fread(FileID, srf.NoVertices, 'int32');

% Read nearest neighbour data for each vertex
srf.NrNeighbours = [];
srf.Neighbours = [];
for CurrVertex=1:srf.NoVertices
    srf.NrNeighbours(CurrVertex) = fread(FileID, 1, 'int32');
    CurrNeighs = fread(FileID, srf.NrNeighbours(CurrVertex), 'int32');
    srf.Neighbours = [srf.Neighbours; CurrNeighs];
end

% Read sequence of three indices to define vertices of each triangle
srf.Triangles = [];
srf.Triangles = fread(FileID, srf.NoTriangles * 3, 'int32');
srf.Triangles = reshape(srf.Triangles, 3, [])';

% Read number of triangle strip elements
srf.NoTriangleStripElems = fread(FileID, [1 1], 'int');

% Read sequence of strip elements
if srf.NoTriangleStripElems > 0
    srf.SequenceStripElements = fread(...
        FileID, [1 NrOfTriangleStripElements], 'int32');
end;

% Read name of MTC file
srf.MTCfile = []; % FIXME

% Close file.
fclose(FileID);
end


text:
function [G, nGrids, nDivs] = GridStructure(n, m)
%   GridStructure - Generate a grid deviding the first octant by polar
%   coordinates. Each element of grid is given by (m-1) integer numbers.
%   The structure is similar to that given by polar coordinates.
%
%   [G, nGrids, nDivs] = GridStructure(n, m)
%
%   Input:
%   n - population size
%   m - number of objectives
%
%   Output:
%   G - grid structure
%   nGrids - number of grids (can be larger than provided population size)
%   nDivs - number of divisions in each right angle
%
%   Example:
%   [G, nGrids, nDivs] = GridStructure(225,3)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Roman Denysiuk

% compute number of divisions
nDivs = ceil( power(n, 1/(m-1)) );

% compute number of grids
nGrids = power(nDivs, m-1);

% initialize grids
G = zeros(nGrids, m-1);

% compute grids
D = m-2;
for j = 0:D
    tmp = repmat(1:nDivs, power(nDivs, j), 1);
    G(:,j+1) = repmat( tmp(:), power(nDivs, D-j), 1);
end

end