text:
function [sma, ecc, inc, raan, arg, tru] = vect_getKeplerFromState(rVect,vVect,gmu, varargin)
% getKeplerFromState() returns Keplerian orbital elements when provided
% with the state (cartesian position vector, cartesian velocity vector) of
% a spacecraft or celestial body.
%
%INPUTS
% rVect - a 3x1 vector that contains the x,y,z components of the orbiting
% body's current position relative to the central body.  Units: [km]
% vVect - a 3x1 vector that contains the x,y,z components of the orbiting
% body's current velocity vector relative to the central body.  Units: [km/sec]
% muCB - the gravitational parameter of the central body.  Units: km^3/s^2
%
%OUTPUTS
% sma - semi-major axis of the orbit.  Units: [km]
% ecc - eccentricity of the orbit.  Units: dimensionless
% inc - inclination angle of the orbit.  Units: radian
% raan - Longitude of ascending node of the orbit.  Units: radian
% arg - Argument of periapse of the orbit.  Units: radian.
% tru - Current true anomaly of the spacecraft/body in the orbit.
% Units: radian
if(~isempty(varargin))
    consistencyCheck = varargin{1};
else
    consistencyCheck = true;
end

numRV = size(rVect,2);

rVect = reshape(rVect,3,numRV);
vVect = reshape(vVect,3,numRV);

if(isscalar(gmu))
    gmu = gmu * ones(1,numRV);
end

try
    [sma, ecc, inc, raan, arg, tru] = vect_getKeplerFromState_Alg_mex(rVect,vVect,gmu(:)');
catch
    [sma, ecc, inc, raan, arg, tru] = vect_getKeplerFromState_Alg(rVect,vVect,gmu(:)');
end

bool = ecc<1.0;

tru(bool)=AngleZero2Pi(tru(bool));
tru(~bool) = angleNegPiToPi(tru(~bool));

tol = 1E-6;
if(consistencyCheck)
    [rVect2,~]=vect_getStatefromKepler(sma, ecc, inc, raan, arg, tru, gmu, false);
    rDang = real(acos(dot(rVect,rVect2)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect2).^2,1)) )));
    if(all(rDang < tol))
        return;
    end

    [rVect3,~]=vect_getStatefromKepler(sma, ecc, inc, raan-rDang, arg+rDang, tru, gmu, false);
    rDang2 = real(acos(dot(rVect,rVect3)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect3).^2,1)) )));
    if(any(rDang2 < tol))
        raan(rDang2 < tol) = raan(rDang2 < tol) - rDang(rDang2 < tol);
        arg(rDang2 < tol)  = arg(rDang2 < tol)  + rDang(rDang2 < tol);
        raan = AngleZero2Pi(raan);
        arg = AngleZero2Pi(arg);

        if(all(rDang2 < tol))
            return;
        end
    end

    [rVect4,~]=vect_getStatefromKepler(sma, ecc, inc, raan+rDang, arg-rDang, tru, gmu, false);
    rDang3 = real(acos(dot(rVect,rVect4)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect4).^2,1)) )));
    if(any(rDang3 < tol))
        raan(rDang3 < tol) = raan(rDang3 < tol) + rDang(rDang3 < tol);
        arg(rDang3 < tol)  = arg(rDang3 < tol)  - rDang(rDang3 < tol);
        raan = AngleZero2Pi(raan);
        arg = AngleZero2Pi(arg);

        if(all(rDang3 < tol))
            return;
        end
    end

    [rVect5,~]=vect_getStatefromKepler(sma, ecc, inc, raan+rDang, arg+rDang, tru, gmu, false);
    rDang4 = real(acos(dot(rVect,rVect5)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect5).^2,1)) )));
    if(any(rDang4 < tol))
        raan(rDang4 < tol) = raan(rDang4 < tol) + rDang(rDang4 < tol);
        arg(rDang4 < tol)  = arg(rDang4 < tol) + rDang(rDang4 < tol);
        raan = AngleZero2Pi(raan);
        arg = AngleZero2Pi(arg);

        if(all(rDang4 < tol))
            return;
        end
    end

    [rVect6,~]=vect_getStatefromKepler(sma, ecc, inc, raan-rDang, arg-rDang, tru, gmu, false);
    rDang5 = real(acos(dot(rVect,rVect6)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect6).^2,1)) )));
    if(any(rDang5 < tol))
        raan(rDang5 < tol) = raan(rDang5 < tol) - rDang(rDang5 < tol);
        arg(rDang5 < tol)  = arg(rDang5 < tol) - rDang(rDang5 < tol);
        raan = AngleZero2Pi(raan);
        arg = AngleZero2Pi(arg);

        if(all(rDang5 < tol))
            return;
        end
    end

    [rVect7,~]=vect_getStatefromKepler(sma, ecc, inc, raan, arg, tru-rDang, gmu, false);
    rDang6 = real(acos(dot(rVect,rVect7)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect7).^2,1)) )));
    if(any(rDang6 < tol))
        tru(rDang6 < tol) = tru(rDang6 < tol) - rDang(rDang6 < tol);

        if(all(rDang6 < tol))
            return;
        end
    end

    [rVect8,~]=vect_getStatefromKepler(sma, ecc, inc, raan, arg, tru+rDang, gmu, false);
    rDang7 = real(acos(dot(rVect,rVect8)./(sqrt(sum(abs(rVect).^2,1)) .* sqrt(sum(abs(rVect8).^2,1)) )));
    if(any(rDang7 < tol))
        tru(rDang7 < tol) = tru(rDang7 < tol) + rDang(rDang7 < tol);

        if(all(rDang7 < tol))
            return;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ vX ] = ProjectOntoHalfSpace( vY, vA, valB )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = OrthogonalProjectionOntoConvexSets( cProjFun, vY, numIterations, stopThr )
%   Solves \arg \min_{x} 0.5 || x - y ||, s.t. x \in \bigcap {C}_{i} using
%   Dykstra's Projection Algorithm.
% Input:
%   - mA            -   Model Matrix.
%                       Input model matrix.
%                       Structure: Matrix (m x n).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - vX            -   Solution Vector.
%                       The solution to the optimization problem..
%                       Structure: Vector (m x 1).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  Orthogonal Projection onto a Half Space - https://math.stackexchange.com/questions/318740.
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     19/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

valR = (vA.' * vY) - valB;

if(valR > 0)
    vX = vY - ((valR / (vA.' * vA)) * vA);
else
    vX = vY;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
classdef PESAII < ALGORITHM
% <multi> <real/integer/label/binary/permutation>
% Pareto envelope-based selection algorithm II
% div --- 10 --- The number of divisions in each objective

%------------------------------- Reference --------------------------------
% D. W. Corne, N. R. Jerram, J. D. Knowles, and M. J. Oates, PESA-II:
% Region-based selection in evolutionary multiobjective optimization,
% Proceedings of the Annual Conference on Genetic and Evolutionary
% Computation, 2001, 283-290.
%------------------------------- 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
        function main(Algorithm,Problem)
            %% Parameter setting
            div = Algorithm.ParameterSet(10);

            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Population.objs,Problem.N,div);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N,div);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [A, rhs] = sp_weak_drchlt_bc (varargin)

  warning ('Using the function sp_weak_drchlt_bc_stokes')
  [A, rhs] = sp_weak_drchlt_bc_stokes (varargin{:});

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function Data = Hash2Struct(hashMap)
% This function converts object of type java.util.HashMap to matlab struct.
% Inputs:
%   - hashMap : java.util.HashMap object
% Outputs:
%   - Data    : matlab struct
% This code is based on template provided by Jan Siroky on
% http://stackoverflow.com/questions/1638018/is-there-already-a-yaml-library-parser-for-matlab

%  %======================================================================
%{
		Copyright (c) 2011
		This program is a result of a joined cooperation of Energocentrum
		PLUS, s.r.o. and Czech Technical University (CTU) in Prague.
        The program is maintained by Energocentrum PLUS, s.r.o. and
        licensed under the terms of MIT license. Full text of the license
        is included in the program release.

        Author(s):
		Jiri Cigler, Dept. of Control Engineering, CTU Prague
		Jan  Siroky, Energocentrum PLUS s.r.o.

        Implementation and Revisions:

        Auth  Date        Description of change
        ----  ---------   --------------------------------------------------
        jc    01-Mar-11   First implementation
        jc    03-Mar-11   Arrays of struct support
        jc    09-Mar-11   DateTime support extended
        jc    15-Mar-11   Support for different encodings
        jc    08-May-11   Support for yaml read of cell of structs
%}
%======================================================================

if isa(hashMap,'java.util.ArrayList')
    iterator = hashMap.iterator();
    index = 0;
    while (iterator.hasNext())
        index = index +1;
        %Data(index)= Hash2Struct(iterator.next());
        res   = Hash2Struct(iterator.next());
        if index==1
            Data = res;
        else
            if iscell(Data)
                Data{index} = res;
            elseif isequal(fieldnames(res),fieldnames(Data))
                Data(index) = res;
            else % convert array to cell
                Data = num2cell(Data);
                Data{index}= res;
            end
        end
    end
else
    Data = [];
    iterator = hashMap.keySet().iterator();
    while (iterator.hasNext())
        field = iterator.next();
        if ~isempty(field)
            d =  hashMap.get(java.lang.String(field));
            switch class(d)
                case {'java.util.LinkedHashMap'}
                    Data.(field) = Hash2Struct(d);
                case {'java.util.ArrayList'}
                    switch class(d.get(0))
                        case {'char' ,'double','java.util.Date'}
                            it = d.iterator();
                            val={};
                            while (it.hasNext())
                                val(end+1)={it.next()};
                                if isa(val{end},'java.util.Date')
                                    val{end} = DateTime(val{end});
                                end
                                if ischar(val{end})
                                    val{end} = DecodeUnicode(val{end});
                                end
                            end
                            if all(cellfun(@(x) isnumeric(x),val))
                                val = cell2mat(val);
                            end
                        case {'java.util.ArrayList'}
                            itr = d.iterator();
                            val={};
                            r=1;
                            while (itr.hasNext())
                                row=itr.next();
                                itc = row.iterator();
                                c=1;
                                while(itc.hasNext())
                                    val{r,c}=itc.next();
                                    if isa(val{r,c},'java.util.Date')
                                        val{r,c} = DateTime(val{r,c});
                                    end
                                    if ischar(val{r,c})
                                        val{r,c} = DecodeUnicode(val{r,c});
                                    end
                                    c=c+1;
                                end
                                r=r+1;
                            end
                            if all(cellfun(@(x) isnumeric(x),val))
                                val = cell2mat(val);
                            end
                        case {'java.util.LinkedHashMap'}
                            for hh=0:(d.size-1)
                                hm=d.get(hh);
                                h2sres = Hash2Struct(hm);
                                h2sres_fn = fieldnames(h2sres);

                                if isfield(Data,field)
                                    if iscell(Data.(field))
                                        Data.(field){hh+1}= h2sres;
                                    elseif isequal(h2sres_fn,fieldnames(Data.(field)))
                                        Data.(field)(hh+1)= h2sres;
                                    else % convert array to cell
                                        Data.(field) = num2cell(Data.(field));
                                        Data.(field){hh+1}= h2sres;
                                    end
                                else
                                    Data.(field)= h2sres;
                                end
                            end
                            continue;

                        otherwise
                            error('unknown  java datatype');
                    end

                    Data.(field) = val;
                otherwise
                    if isa(d,'java.util.Date')
                        d = DateTime(d);
                    end
                    if ischar(d)
                        d = DecodeUnicode(d);
                    end
                    Data.(field) = d;
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
%
% Convert a polynomial coefficient matrix to a polynomial string
%
% Syntax:  >> p = CoefMat2PolynString(F,var)
%          >> p = CoefMat2PolynString(F,var,digits)
%
% Input:      F --- (matrix) Coefficient matrix of a polynomial
%           var --- (cell)   variable names in order
%        digits --- (integer) optional, number of digits, default = 15
%
% Output:     p --- (string) polynomial as a string
%
% Example:  >> F = [2 0; 1 3; 6, 8]
%              F =
%                2     0
%                1     3
%                6     8
%           >> CoefMat2PolynString(F,{'x','y'})
%              ans =
%                  6*x^2*y + 8*y^3

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%  rise_exit exit from RISE
%
%  ::
%
%
%    rise_exit
%    rise_exit()
%
%  Args:
%
%     none
%
%  Returns:
%     :
%
%     none
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function [x,paramvec,paramindex,ons] = cond2reg(c,mspec,paramvec,paramindex)
% given a condition structure, returns regressor
%
% x is reg or set of regs to add to model
% paramvec is cell array of onsets (relative to...) for design reconstruction
% paramindex is index of which cell in paramvec to use for reconstruction of which condition
% ons is vector of trial onsets, for updating conditions, if necessary
%
% This function is not used in GA2: it's an extra function that mirrors a subfunction of
% construct_model.m
%
% THIS IS OLD! DOES IT WORK?
%
% Tor Wager

    if isfield(c,'parts')
        if c.parts ~= 1
            % can't build regressor from trial type with parts - return
            x = [];
           return
       end
    end

    % determine onset times
    % -----------------------------------------------------------------------------------
    delta = zeros(mspec.numframes,1);
    if size(c.onsets,1) == 1
        % fixed onsets for condition, but sub-trial may be variable

        if isfield(c,'baseons')     % for sub-parts of trials
            if size(c.onsets,2) == 2    % sub-trial, range is given; get random
                if length(paramvec) < paramindex,
                    paramvec{paramindex} = round(randrange(c.onsets,length(c.baseons)));
                end
                c.onsets = paramvec{paramindex};
                paramindex = paramindex + 1;
            end
            ons = c.onsets + c.baseons;

        else
            ons = c.onsets;

        end
        ons = round(ons);
        delta(1 + (ons)) = 1;

    else
        % variable (random) onsets
        if length(paramvec) < paramindex
            % build new ones if they're not input
            paramvec{paramindex} = onsetbuilder(c.onsets,mspec.numframes);
        end

        if isfield(c,'baseons')     % for sub-parts of trials
            ons = paramvec{paramindex} + c.baseons;
        else
            ons = paramvec{paramindex};
        end
        delta(1 + (paramvec{paramindex})) = 1;
        paramindex = paramindex + 1;
    end

    % convolve or shift
    % -----------------------------------------------------------------------------------

    if isfield(c,'parts')
        if c.parts == 1
            doconv = c.subcond(1).convolve;
            hrfest = c.subcond(1).hrfest;
        else
            error('This should never ever happen.')
        end
    else
        % this c IS a subpart
        doconv = c.convolve;
        hrfest = c.hrfest;
    end

    if doconv,
         x = conv(mspec.hrf,delta);
    elseif ~isempty(hrfest)
         sf{1} = delta;
         [x] = tor_make_deconv_mtx(sf,round(hrfest),1);
    else
         x = delta;
    end

    x = x(1:ceil(mspec.numframes),:);

return

% -----------------------------------------------------------------------------------
% * fill the run with as many trials as possible, given range of random lengths and delays
% -----------------------------------------------------------------------------------

function [onsets] = onsetbuilder(range,numframes)
% builds list of trials with random onset times
% range is range of onset times in condition.onsets, a 2 x 2 matrix

onsets = []; tend = 0;  % 0 is first time point in run
dlen = round(randrange(range(1,:),1));
dnaindex = 2;

while tend + dlen < numframes     % while onset of next trial is < end of run

    tlen = round(randrange(range(2,:),1));  % trial length
    if tend + dlen + tlen > numframes, break,end  % trial will not fit in run

    % add a trial
    onsets(end+1) = tend + dlen;            % specify onset after start delay of dlen
    tend = onsets(end) + tlen;              % specify ending point of trial
    dlen = round(randrange(range(1,:),1));  % re-randomize delay for next trial

end

% -----------------------------------------------------------------------------------
% * get a random number or vector within a specified range
% -----------------------------------------------------------------------------------
function rval = randrange(range,num)
    rval = rand(1,num) * (range(1,2) - range(1,1)) + range(1,1);
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [kJ] = J2kJ(J)
% Convert energy or work from joules to joules.
% Chad A. Greene 2012
kJ = J*.001;

author: Sable

repo: mcbench-benchmarks


text:
function reverse_mask(inname,outname)
% Changes 1 and greater's to 0's in an .img file, and vice versa
% NaNs are still NaNs.
%
% :Usage:
% ::
%
%     reverse_mask(inname,outname)

Q = spm_imcalc_ui(inname,outname,'abs(i1) < eps');

return


text:
function x = Ftrunc( x, prec)

% function x = Ftrunc( x, prec)
%
% rounds x to prec significant digits

  f = 10.^fix( prec - log10( abs( x)));
  x = round( x.*f) ./ f;


text:
function a = ipermute(a,order)
%IPERMUTE     Inverse permute array dimensions
%
%Implements the Matlab function "ipermute", similar syntax and semantics
%

% written  09/17/08     S.M. Rump
%

  if isa(a,'intval')
    if a.complex
      a.mid = ipermute(a.mid,order);
      if ~isequal(a.rad,0)
        a.rad = ipermute(a.rad,order);
      end
    else
      a.inf = ipermute(a.inf,order);
      a.sup = ipermute(a.sup,order);
    end
  else
    error('invalid call of permute for intervals')
  end

repo: OpenMAS


text:
%                                          [ux,ispos] = psdfactor(x,K)
% PSDFACTOR  UX'*UX Cholesky factorization
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

function [ux,ispos] = psdfactor(x,K) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
report_this_filefun(mfilename('fullpath'));

%l = isnan(tmap);
%tmap(l) = 1;

l = tmap < 100;
tmap(l) = 100;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
%  [X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

re4 = re3;

l = re4 < -3.9;
re4(l) = -3.9;
l = re4 > -1.4;
re4(l) = -1.3;

ren = interp2(X2,Y2,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-20;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[32.6 36.2],'MapLonLimit',[-118.3 -115.3])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',6);
tightmap
view([0 90])
hl = camlight ; lighting phong
set(gca,'projection','perspective');

aa = a;

for i = 1:length(maepi(:,3))
    dep = interp2(lon,lat,tmap,maepi(i,1),maepi(i,2));

    pl =plot3m(maepi(i,2),maepi(i,1),dep+175,'^k');
    set(pl,'Markersize',14,'markerfacecolor','w');

end

load usahi
[plat,plon] = extractm(stateline,'california');
pl =plot3m(plat,plon,175,'linewidth',2);
set(pl,'color','w');

%dep = interp2(lon,lat,tmap,faults(:,1),faults(:,2));
%pl = plotm(faults(:,2), faults(:,1),dep+10,'k','Linewidth',1.0);

dep = interp2(lon,lat,tmap,mainfault(:,1),mainfault(:,2));
pl = plotm(mainfault(:,2), mainfault(:,1),dep+10,'y','Linewidth',2.0);

%zdatam(handlem('allline'),10000) % keep line on surface
j = jet;
j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j ; [0.9 0.9 0.9] ];
caxis([ -4 -1.3]);

colormap(j);
axis off; set(gcf,'color','k')
ax = axis;
axis([ax(1:4) 0 5000])

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',1)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',1)
setm(gca,'parallellabel','on','Labelrotation','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',14,'Labelunits','degrees')

h5 = colorbar;
set(h5,'position',[0.7 0.21 0.01 0.2],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',14);
set(gcf,'Inverthardcopy','off');

a = aa;

author: CelsoReyes


text:
x = [-117.6 -115.6 -115.6 -117.6 -117.6]
y = [35.5 35.5 33.5 33.5 35.5]

%create a rectangular grid
xvect=[min(x):dx:max(x)];
yvect=[min(y):dy:max(y)];
gx = xvect;
gy= yvect;
tmpgri=zeros((length(xvect)*length(yvect)),2);
n=0;
for i=1:length(xvect)
    for j=1:length(yvect)
        n=n+1;
        tmpgri(n,:)=[xvect(i) yvect(j)];
    end
end
%extract all gridpoints in chosen polygon
XI=tmpgri(:,1);
YI=tmpgri(:,2);

    ll = polygon_filter(x,y, XI, YI, 'inside');
%grid points in polygon
newgri=tmpgri(ll,:);

% plot the grid points
figure_w_normalized_uicontrolunits(map)
pl = plot(newgri(:,1),newgri(:,2),'+k','era','normal');
set(pl,'MarkerSize',8,'LineWidth',1)
drawnow


text:
function [assignment, cost] = assignmentoptimal(distMatrix)
%ASSIGNMENTOPTIMAL    Compute optimal assignment by Munkres algorithm
%		ASSIGNMENTOPTIMAL(DISTMATRIX) computes the optimal assignment for the
%		given rectangular distance (or weight) matrix, for example the assignment
%		of tracks (in rows) to observations (in columns). The result is a column
%		vector containing the assigned column number in each row (or 0 if no
%		assignment could be done).
%
%		[ASSIGNMENT, COST] = ASSIGNMENTOPTIMAL(DISTMATRIX) returns the assignment
%		vector and the overall cost.
%
%		The distance matrix may contain infinite values (forbidden
%		assignments). Internally, the infinite values are set to a very large
%		finite number, so that the Munkres algorithm itself works on finite-number
%		matrices. Before returning the assignment, all assignments with infinite
%		distance are deleted (i.e. set to zero).
%
%		A description of Munkres algorithm (also called Hungarian algorithm) can
%		easily be found on the web.
%
%		Written by Markus Buehren, www.Lss.uni-stuttgart.de
%		Last modified 14.12.2004

% save original distMatrix for cost computation
originalDistMatrix    = distMatrix;

% check for negative elements
if any(distMatrix(:) < 0)
	error('All matrix elements have to be non-negative.');
end

% get matrix dimensions
[nOfRows, nOfColumns] = size(distMatrix);

% check for infinite values
finiteIndex   = isfinite(distMatrix);
infiniteIndex = find(~finiteIndex);
if ~isempty(infiniteIndex)
	% set infinite values to large finite value
	maxFiniteValue = max(max(distMatrix(finiteIndex)));
	if maxFiniteValue > 0
		infValue = abs(10 * maxFiniteValue * nOfRows * nOfColumns);
	else
		infValue = 10;
	end
	if isempty(infValue)
		% all elements are infinite
		assignment = zeros(nOfRows, 1);
		cost       = 0;
		return
	end
	distMatrix(infiniteIndex) = infValue;
end

% memory allocation
coveredColumns = zeros(1,nOfColumns);
coveredRows    = zeros(nOfRows,1);
starMatrix     = zeros(nOfRows, nOfColumns);
primeMatrix    = zeros(nOfRows, nOfColumns);

% preliminary steps
if nOfRows <= nOfColumns
	minDim = nOfRows;

	% find the smallest element of each row
	minVector = min(distMatrix,[],2);

	% subtract the smallest element of each row from the row
	distMatrix = distMatrix - repmat(minVector, 1, nOfColumns);

	% Steps 1 and 2
	for row = 1:nOfRows
		for col = find(distMatrix(row,:)==0)
			if ~coveredColumns(col)%~any(starMatrix(:,col))
				starMatrix(row, col) = 1;
				coveredColumns(col)  = 1;
				break
			end
		end
	end

else % nOfRows > nOfColumns
	minDim = nOfColumns;

	% find the smallest element of each column
	minVector = min(distMatrix);

	% subtract the smallest element of each column from the column
	distMatrix = distMatrix - repmat(minVector, nOfRows, 1);

	% Steps 1 and 2
	for col = 1:nOfColumns
		for row = find(distMatrix(:,col)==0)'
			if ~coveredRows(row)
				starMatrix(row, col) = 1;
				coveredColumns(col)  = 1;
				coveredRows(row)     = 1;
				break
			end
		end
	end
	coveredRows(:) = 0; % was used auxiliary above

end

if sum(coveredColumns) == minDim
	% algorithm finished
	assignment = buildassignmentvector__(starMatrix);
else
	% move to step 3
	[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);
end

% compute cost and remove invalid assignments
[assignment, cost] = computeassignmentcost__(assignment, originalDistMatrix, nOfRows);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function assignment = buildassignmentvector__(starMatrix)

[maxValue, assignment] = max(starMatrix, [], 2);
assignment(maxValue == 0) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, cost] = computeassignmentcost__(assignment, distMatrix, nOfRows)

rowIndex   = find(assignment);
costVector = distMatrix(rowIndex + nOfRows * (assignment(rowIndex)-1));
finiteIndex = isfinite(costVector);
cost = sum(costVector(finiteIndex));
assignment(rowIndex(~finiteIndex)) = 0;

% Step 2: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step2__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% cover every column containing a starred zero
maxValue = max(starMatrix);
coveredColumns(maxValue == 1) = 1;

if sum(coveredColumns) == minDim
	% algorithm finished
	assignment = buildassignmentvector__(starMatrix);
else
	% move to step 3
	[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);
end

% Step 3: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

zerosFound = 1;
while zerosFound

	zerosFound = 0;
	for col = find(~coveredColumns)
		for row = find(~coveredRows')
			if distMatrix(row,col) == 0

				primeMatrix(row, col) = 1;
				starCol = find(starMatrix(row,:));
				if isempty(starCol)
					% move to step 4
					[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step4__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim);
					return
				else
					coveredRows(row)        = 1;
					coveredColumns(starCol) = 0;
					zerosFound              = 1;
					break % go on in next column
				end
			end
		end
	end
end

% move to step 5
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step5__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

% Step 4: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step4__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim)

newStarMatrix          = starMatrix;
newStarMatrix(row,col) = 1;

starCol = col;
starRow = find(starMatrix(:, starCol));

while ~isempty(starRow)

	% unstar the starred zero
	newStarMatrix(starRow, starCol) = 0;

	% find primed zero in row
	primeRow = starRow;
	primeCol = find(primeMatrix(primeRow, :));

	% star the primed zero
	newStarMatrix(primeRow, primeCol) = 1;

	% find starred zero in column
	starCol = primeCol;
	starRow = find(starMatrix(:, starCol));

end
starMatrix = newStarMatrix;

primeMatrix(:) = 0;
coveredRows(:) = 0;

% move to step 2
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step2__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

% Step 5: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step5__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% find smallest uncovered element
uncoveredRowsIndex    = find(~coveredRows');
uncoveredColumnsIndex = find(~coveredColumns);
[s, index1] = min(distMatrix(uncoveredRowsIndex,uncoveredColumnsIndex));
[s, index2] = min(s);
h = distMatrix(uncoveredRowsIndex(index1(index2)), uncoveredColumnsIndex(index2));

% add h to each covered row
index = find(coveredRows);
distMatrix(index, :) = distMatrix(index, :) + h;

% subtract h from each uncovered column
distMatrix(:, uncoveredColumnsIndex) = distMatrix(:, uncoveredColumnsIndex) - h;

% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = step3__(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

author: kristinbranson

repo: JAABA


text:
function [hdr] = read_bucn_nirshdr(filename)

% READ_BUCN_NIRSHDR reads the header information of ASCII-formatted NIRS
% data acquired with the UCL-BIRKBECK machine and postprocessed by the
% Paris group. The first line contains the channel labels and the rest of
% the file contains per line a time sample. The first column specifies the
% time axis.
%
% Use as
%   [hdr] = read_bucn_nirshdr(filename)
%
% See also read_bucn_nirsdata, READ_BUCN_NIRSEVENT

% Copyright (C) 2011, Jan-Mathijs Schoffelen
%
% 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: read_bucn_nirshdr.m$

fid = fopen_or_error(filename, 'r');

% read the first line
line1 = textscan(fid, '%[^\n]',1);

% field delimiter can be space or tab
labelspc = textscan(line1{1}{1}, '%[^ ]');
labeltab = textscan(line1{1}{1}, '%[^\t]');

% let tab as a delimiter prevail
if numel(labeltab{1})>1
  label = labeltab{1};
else
  label = labelspc{1};
end
nchan = numel(label);
Fs    = str2num(strtok(strtok(label{1},'#Time.'),'Hz'));

% test whether the channel labels are non-numeric
labelnumber = cellfun(@str2num, label, 'UniformOutput', false);
labelstring = cellfun(@isempty, labelnumber, 'UniformOutput', true);
if ~any(labelstring)
  ft_error('channel labels were not found in the first line of the file');
end

% read the rest
dat = textscan(fid, '%f');
fclose(fid);

dat  = reshape(dat{1}, nchan, []);
nsmp = size(dat,2);

% create the output
hdr          = [];
hdr.Fs       = Fs;
hdr.label    = label;
hdr.nTrials  = 1;
hdr.nSamples = nsmp;
hdr.nSamplesPre = 0;
hdr.nChans   = nchan;
hdr.time     = dat(1,:); % events in the raw event file have both a sample and a time stamp

author: spm


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function area = box_area( bbox )
    up    = bbox(1);
    left  = bbox(2);
    down  = bbox(3);
    right = bbox(4);

    area = (down-up+1)*(right-left+1);
end

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
%%%C-Type implementation for imu interpolation for iphone data
%oper=new output period
%cdef= col definitions [total, Time , x , y, z]
function [ini_time]=interp_imu(INFILES, OUTFILE, cdef, oper)

% fclose all;
% clear;
% DIRNAME='C:\Documents and Settings\Yigiter\Desktop\data\CalibTest1\';
% INFILES=[[DIRNAME '07-16-00-12-acclog.bin ']; [DIRNAME '07-16-00-12-gyrolog.bin']];
% OUTFILE=[DIRNAME 'imu.bin'];
% cdef=[5 2 3 4 5;5 2 3 4 5];
% oper=2;

%i/o files
nfiles=size(cdef,1);
ndat=size(cdef,2)-1;

infids=zeros(nfiles,1);
for (in=1:nfiles)
    infids(in,1)=fopen(INFILES(in,:),'rb');
end
ofid=fopen(OUTFILE,'wb');

%first data from all files
data=zeros(ndat,nfiles);
for in=1:nfiles
    vr_a=fread(infids(in),[cdef(in,1) 1],'double');
    data(:,in)=vr_a(cdef(in,2:end),1);
end

%Latest data
ct=data(1,1);
for in=2:nfiles
    if ct<data(1,in)
        ct=data(1,in);
    end
end

%Read all files till they past the latest
flag=ones(nfiles,1);
while (sum(flag,1))
    for in=1:nfiles
        if (data(1,in)<ct)
            vr_a=fread(infids(in),[cdef(in,1) 1],'double');
            data(:,in)=vr_a(cdef(in,2:end),1);
        else
            flag(in)=0;
        end
    end
end

ini_time=ct;

%%Write the first data (useless)
oin=0;
vr_a=reshape(data(2:end,:), (ndat-1)*nfiles, 1);
fwrite(ofid, [oin;vr_a],'double');

%Start the main part
ipt=ones(nfiles,1)*ct;
idt=zeros(nfiles,1);
data_inc=zeros(ndat,nfiles);
iflag=zeros(nfiles,1);
fflag=zeros(nfiles,1);

while (~sum(fflag,1))
    %%process all files till they exceed the next sampling point
    for in=1:nfiles
        if (~iflag(in))
            dt=data(1,in)-ipt(in);
            if (idt(in)+dt<=oper)
                data_inc(:,in)=data_inc(:,in)+dt*data(:,in);
                idt(in)=idt(in)+dt;

                ipt(in)=data(1,in);
                vr_a=fread(infids(in),[cdef(in,1) 1],'double');
                if (isempty(vr_a))
                    break;  %Why does not feof() work properly! whyyyyyyy
                end
                data(:,in)=vr_a(cdef(in,2:end),1);
            else
                data_inc(:,in)=data_inc(:,in)+(oper-idt(in))*data(:,in);
                ipt(in)=ipt(in)+(oper-idt(in));
                iflag(in)=1;
            end
        end
    end

    %Check if one the files reaches the feof
    for in=1:nfiles
        fflag(in)=feof(infids(in));
    end

    %%In all files pass the next sampling point, write the results
    if sum(iflag,1)==nfiles
        %%Write data
        oin=oin+1;
        vr_a=reshape(data_inc(2:end,:), (ndat-1)*nfiles, 1)/oper;
        fwrite(ofid, [oin;vr_a],'double');

        %Prepare for the next cycle
        idt=zeros(nfiles,1);
        data_inc=zeros(ndat,nfiles);
        iflag=zeros(nfiles,1);
    end
end

%%Close all the files
for in=1:nfiles
    fclose(infids(in));
end
fclose(ofid);

return;

repo: ekfmonoslam


text:
function [in2] = ft22in2(ft2)
% Convert area from square feet to square inches.
% Chad A. Greene 2012
in2 = ft2*144;


text:
function [f,g] = autoGrad(x,useComplex,funObj,varargin)
% [f,g] = autoGrad(x,useComplex,funObj,varargin)
%
% Numerically compute gradient of objective function from function values

p = length(x);
mu = 1e-150;

if useComplex % Use Complex Differentials
    diff = zeros(p,1);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        diff(j,1) = funObj(x + mu*i*e_j,varargin{:});
    end

    f = mean(real(diff));
    g = imag(diff)/mu;
else % Use Finite Differencing
    f = funObj(x,varargin{:});
    mu = 2*sqrt(1e-12)*(1+norm(x))/norm(p);
    for j = 1:p
        e_j = zeros(p,1);
        e_j(j) = 1;
        diff(j,1) = funObj(x + mu*e_j,varargin{:});
    end
    g = (diff-f)/mu;
end

if 0 % DEBUG CODE
    [fReal gReal] = funObj(x,varargin{:});
    [fReal f]
    [gReal g]
    pause;
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
%kkresample 'Resample Object Data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kresample.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input Object'
% Toggle: center 'Place COA in middle of dimensions (unless overridden below)', default: 0: 'center'
% Integer: wc 'Width Coordinate of COA ', default: 0: 'Width coordinate for the center of action'
% Integer: hc 'Height Coordinate of COA ', default: 0: 'Height coordinate for the center of action'
% Integer: dc 'Depth Coordinate of COA ', default: 0: 'Depth coordinate for the center of action'
% Integer: tc 'Time Coordinate of COA ', default: 0: 'Time coordinate for the center of action'
% Integer: ec 'Elements Coordinate of COA ', default: 0: 'Elements coordinate for the center of action'
% OutputFile: o 'Output', required: 'Resulting output object'
%
% Example: o = kkresample(i, {'i','';'center',0;'wc',0;'hc',0;'dc',0;'tc',0;'ec',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kresample - Resample Object Data
%
%  DESCRIPTION
% .I kresample
% is used to magnify or reduce the Input data object (i) by using the
% selected interpolation method (inter).  (At this time, only pixel
% replication is supported.)  Magnification factors greater than one
% magnify the object, and magnification factors between zero and one
% reduce the object.
%
% Pixel replication and subsampling are done on a floating point indexing
% basis, where the data point index is obtained by truncating the cumulative
% floating point index.  Therefore, when the calculated replication factor
% or the subsampling step size is non-integer, the resulting integer step
% size varies during the process.
%
%  "Magnification Factors"
% Magnification factors for each dimension, Width (wmag), Height (hmag)
% Depth (dmag), Time (tmag), and Elements (emag) can be specified.  Magnification
% factors can be any value greater than zero.
%
%  "Center of Action"
% The magnification or reduction is done about a center point that can be
% specified to be anywhere in the object (even outside of the object).
% The correct way to think about this is by analogy to a magnifying glass
% where you place the magnifying glass over the point of interest and then
% choose the lens to get the magnification you want.  The center of action
% coordinates can be specified individually for each dimension (wc, hc, dc,
% tc, ec) or, for convenience, the center of action can be automatically
% computed to be in the middle of all dimensions.  This is accomplished
% by selecting the "center" flag.  Automatic centering can be overridden
% on selected dimensions by explicitly specifying the centers of action
% for those dimensions.
%
%  "Output Object Size"
% The resultant output object size can either be the same as the original
% input object size if "Resize dimensions" is set to "No" (-resize 0), or
% it can be resized to fit the magnified or reduced data set if "Resize
% dimensions" is set to "Yes" (-resize 1).
% If the user chooses not to resize the dimensionality of the output object,
% and the input object is reduced (magnification
% factor < 1.0) then unknown pixels are set to zero.  Likewise, if the
% input object is magnified (magnification factor > 1.0) then input data
% that would be mapped outside of the output object are lost.
%
%  "Data Type"
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
% If there is map and value data in the source object, and zero
% order interpolation (pixel replication) is used, the map is transferred
% to the destination object.  If any interpolation algorithm is used,
% the value data is mapped before the operation can begin and the
% destination will not have map data.
%
%  "Failure Modes"
% .cI $DATAMANIP/repos/shared/man/sections/fail_1input
% This program fails if there is mask data and the operation uses anything
% other than zero order interpolation.  In other words, if value data
% is being replicated or subsample on grid, then the mask can be
% correspondingly replicated or subsampled.  (Bilinear interpolation can not
% be performed on mask data.)
% If explicit location or time data are available, the sampling algorithm
% can not resample these data yet, and will fail.  This will be implemented
% later ....
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% kshrink, kexpand
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkresample(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkresample(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'center', 0;'wc', 0;'hc', 0;'dc', 0;'tc', 0;'ec', 0;'o', '__output'};
maxval={0,0,1,1,1,1,1,0};
minval={0,0,1,1,1,1,1,0};
istoggle=[0,1,1,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Toggle','Integer','Integer','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kresample"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function paf = graph2paf(nodes, edges, sz, channelsOnly, sigma)
%GRAPH2PAF Converts a set of edges into part affinity fields.
% Usage:
%   graph2paf(nodes, edges, sz, sigma)
%
% Args:
%   nodes: set of points (N x 2)
%   edges: indices into nodes defining directed edges (E x 2)
%   sz: grid/image size (1 x 2)
%   channelsOnly: stack all PAFs along channels (dim 3) instead of dim 4 (default: true)
%   sigma: maximum distance from edge to keep (default: 5)
%
% Returns:
%   paf: part affinity fields (sz(1) x sz(2) x 2E) or (sz(1) x sz(2) x 2 x E)
%
% See also: pts2confmaps

if nargin < 4 || isempty(channelsOnly); channelsOnly = true; end
if nargin < 5 || isempty(sigma); sigma = 5; end

% Create image coordinate grid
[XX,YY] = meshgrid(1:sz(2), 1:sz(1));

% Create PAFs for each edge
E = size(edges,1);
paf = cell(E,1);
for i = 1:E
    % Pull out edge points
    src = nodes(edges(i,2),:);
    dst = nodes(edges(i,1),:);

    % Edge length
    L = norm(dst - src, 2);

    % Unit vectors
    V = (dst - src) ./ L; % pointing along edge
    Vp = [-V(:,2), V(:,1)]; % perpendicular

    % Signed distance along edge
    D1 = sum(V .* ([XX(:) YY(:)] - src),2);

    % Absolute distance orthogonal to edge
    D2 = abs(sum(Vp .* ([XX(:) YY(:)] - src),2));

    % Vector field mask
    paf_mask = reshape(D1 >= 0 & D1 <= L & D2 <= sigma, sz);

    % Create vector field along channels (X and Y)
    paf{i} = paf_mask .* permute(V, [1 3 2]);
end

% Merge all edge PAFs
if channelsOnly
    paf = cat(3, paf{:});
else
    paf = cat(4, paf{:});
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function ll = composite_likelihood(S,data,indices)
% Computes the negative of the composite normal loglikelihood (bivariate) for a K-dimensional array
%
% USAGE:
%   [LL] = composite_likelihood(S,DATA,INDICES)
%
% INPUTS:
%   S       - K by K covariance matrix
%   DATA    - Either a K by K matrix (e.g. outer-produces) or a 1 by K vector (e.g. returns)
%   INDICES - Q by 2 array of indices to use when computing the composite likleihood
%
% OUTPUTS:
%   LL      - Composite likelihood
%
% COMMENTS:
%  This is a helper function for various multivariate GARCH models. A MEX file version of this file
%  is available which provides a large speed-up.
%
%  See also DCC, SCALAR_VEC_VECH, BEKK, RARCH

q = size(indices,1);
[m,n] = size(data);
likConst = 3.67575413281869;
ll = 0;
if m==n
    for k=1:q
        i = indices(k,1);
        j = indices(k,2);
        s11 = S(i,i);
        s12 = S(i,j);
        s22 = S(j,j);
        det = s11*s22-s12*s12;
        x11 = data(i,i);
        x12 = data(i,j);
        x22 = data(j,j);
        ll = ll + 0.5*(likConst + log(det) + (s22*x11 - 2*s12*x12 + s11*x22)/det)/q;
    end
else
    for k=1:q
        i = indices(k,1);
        j = indices(k,2);
        s11 = S(i,i);
        s12 = S(i,j);
        s22 = S(j,j);
        det = s11*s22-s12*s12;
        x11 = data(i)*data(i);
        x12 = data(i)*data(j);
        x22 = data(j)*data(j);
        ll = ll + 0.5*(likConst + log(det) + (s22*x11 - 2*s12*x12 + s11*x22)/det)/q;
    end
end

author: bashtage

repo: mfe-toolbox


text:
function res = copy(this, newname)
% Method for copying a dataset
% FORMAT res = copy(this, fname)
%
% fname can be
% - path\filename -> data copied and renamed
% - path          -> data copied only
%__________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: copy.m 5025 2012-10-31 14:44:13Z vladimir $

[p, f] = fileparts(newname);

if ~isempty(p)
    if ~exist(p,'dir'), mkdir(p); end;
else
    p = path(this);
end

if isempty(f)
    f = fname(this);
else
    f = [f '.mat'];
end

if strcmpi(fullfile(this), fullfile(p, f))
    res = this;
    return;
end

%-Copy dataset (.mat and .dat)
%--------------------------------------------------------------------------
new = clone(this, fullfile(p, f));
[r, msg] = copyfile(fnamedat(this), ...
    fnamedat(new), 'f');
if ~r
    error(msg);
    res = [];
else
    res = new;
end


text:
function mColormap = gui_Colormap_RedGrayGreen(nSize)

% If size is not specified, set it to 256
if nargin < 1
  nSize = 256;
end
nPart = floor(nSize / 2);

mRedPart = [];
mGreenPart = [];

for nCnt = 1:nPart
  mRedPart = [mRedPart; 0.75  (nCnt*(0.75/nPart)) (nCnt*(0.75/nPart))];
  mGreenPart = [mGreenPart; (0.75 - (nCnt*(0.75/nPart))) (0.75 - (nCnt*(0.75/nPart))) 0.75];
end

mColormap = [mRedPart; mGreenPart];


text:
function plot_3d_odomap(groundTruthFile, odoMapFile, ...
    xOdoMapScaling, yOdoMapScaling, zOdoMapScaling, ...
    xOdoMapTrans, yOdoMapTrans, zOdoMapTrans, ...
    xGtScaling, yGtScaling, zGtScaling)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    % load ground truth data
    [frameId, gt_x, gt_y, gt_z, gt_rx, gt_ry, gt_rz] = load_ground_truth_data(groundTruthFile);

    % load experience map data
    [id, odomap_x, odomap_y, odomap_z, odomap_yaw] = load_exp_map_data(odoMapFile);

    figure('color',[1 1 1]);
    hold on

    plot3(odomap_y * yOdoMapScaling + yOdoMapTrans, ...
        odomap_x * xOdoMapScaling + xOdoMapTrans, ...
        odomap_z * zOdoMapScaling + zOdoMapTrans, '.b', 'MarkerSize',7);
    plot3((gt_x - gt_x(1)) * xGtScaling,(gt_y - gt_y(1)) * yGtScaling, (gt_z - gt_z(1)) * zGtScaling, '.r','MarkerSize',7);
    hold off;

    view(3)
% view(0, 90);
    grid on
%     xlabel('x');
%     ylabel('y');
    xl = xlabel('x','FontSize',18);
    yl = ylabel('y','FontSize',18);
    zlabel('z','FontSize',18);
    set(xl,'Rotation',15);
    set(yl,'Rotation',-30);
%     title('3D Experience Map');
    %                     legend('Result','Truth' ,'1');
    % SynPanData
    axis([-2 16 -20 5 -1 3]);

% SynPerData
%     axis([-5 20 -20 5 -1 3]);
    %                     axis equal
%     legend('Experience map','Ground truth', 'location', '0');
    hh = legend('OM','GT');

    set(hh,'position',[.76 .82 .05 .05]);
    set(gca,'FontSize',18); % axis font
    set (gcf,'Position',[500,300,400,300], 'color','w')
    axis on
    rotate3d on
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function TS_LoadFeatures(whatDataFile,TS_DataMatToLoad,TS_CalcTimeToLoad,TS_QualityToLoad,INP_opsToLoad,INP_tsToLoad)
% TS_LoadFeatures   Load some pre-computed features into an hctsa data file
%                     if using externally computed values, e.g., from Catch22.
%
%---INPUTS:
% whatDataFile, the input/output hctsa file (already created with TS_Init)
% TS_DataMatToLoad, the feature matrix to load, as a time series-by-features
%               matrix
% TS_CalcTimeToLoad, the calculation time for each operation, same size as
%                       TS_DataMatToLoad
% TS_QualityToLoad, the quality for each operation (zero means no issues), same
%                       size as TS_DataMatToLoad
% INP_opsToLoad, the names of the hctsa operations, calculation and quality
%                   matrices to be loaded (see INP_ops.txt)
% INP_tsToLoad, the names of the time series that are in whatDataFile that
%                   we're loading
%
% whatDataFile should contain all operations specified in TS_opsToLoad

% ------------------------------------------------------------------------------
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% whatDataFile must be a filename

% Load in the data:
TS_DataMat = TS_GetFromData(whatDataFile,'TS_DataMat');
TS_CalcTime = TS_GetFromData(whatDataFile,'TS_CalcTime');
TS_Quality = TS_GetFromData(whatDataFile,'TS_Quality');

% Which operations are we loading?
if isempty(INP_opsToLoad)
    myOpIDs = 1:size(TS_DataMat,2);
else
    myOpIDs = TS_GetIDs(INP_opsToLoad,whatDataFile,'ops','Name');
end

% Which time series are we loading?
if isempty(INP_tsToLoad)
    myTSIDs = 1:size(TS_DataMat,1);
else
    myTSIDs = TS_GetIDs(INP_tsToLoad,whatDataFile,'ts','Name');
end

% Fill the data:
TS_DataMat(myTSIDs,myOpIDs) = TS_DataMatToLoad;
TS_CalcTime(myTSIDs,myOpIDs) = TS_CalcTimeToLoad;
TS_Quality(myTSIDs,myOpIDs) = TS_QualityToLoad;

% Save to the data file
save(whatDataFile,'TS_DataMat','TS_CalcTime','TS_Quality','-append');

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Waveform synthesis for the Harmonic Model + Phase Distortion (HMPD)
%
% Octave compatible
%
% This the main entry point for synthesizing a waveform using HMPD parameters.
%
% Please read the README.txt file for general information about HMPD before
% using it.
%
% Inputs
%  f0s    : [s, Hz] [Nx2] A time/data column vector, containing the
%           fundamental frequency f0 values.
%  AE     : [NxD] A matrix containing the amplitude envelope.
%           D is either opt.dftlen/2+1 (from hmpd_features_compute.m), or
%           opt.amp_order+1, depending if compression is disabled or enabled.
%  PDM    : [NxD] A matrix containing the Phase Distortion Mean.
%           D is either opt.dftlen/2+1 or opt.pdm_order+1, depending if
%           compression is disabled or enabled.
%  PDD    : [NxD] A matrix containing the Phase Distortion Deviation.
%           D is either opt.dftlen/2+1 or opt.pdd_order+1, depending if
%           compression is disabled or enabled.
%  fs     : [Hz] The sampling rate of the waveform.
%           (It has to be the same as the one used during analysis)
% [wavlen]: Length of the synthesized waveform [in number of samples].
%           If this option is omited, it is predicted from the last time instant
%           of the f0s argument.
% [opt]   : Additional options for synthesis (see code below).
%           The option opt.enc has to contain absolutely the options used
%           during analysis.
%
% Outputs
%  syn   : The samples of the synthesized waveform
%  opt   : The options used during synthesis.
%
% Copyright (c) 2012 University of Crete - Computer Science Department (UOC-CSD)
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <degottex@csd.uoc.gr>
%

function [syn, opt] = hmpd_synthesis(f0s, AE, PDM, PDD, fs, wavlen, opt)

    if nargin<7
        % Options
        opt.enc = hmpd_analysis();
        opt.amp_minphasegain = 1; % 1:minimum-phase, 0:zero-phase, -1:max-phase
        opt.pd_gain = 1;
        opt.rps_randunvoiced = false; % Randomize the RPS in unvoiced segments
                                      % Unvoiced segments are indicated by
                                      % f0=0 values.
                                      % (not used by default with HMPD)
        opt.pdv_corrthresh = 0.75;    % PDD correction threshold

        opt.defdbamp = -300;  % [dB]
        opt.usemex   = false; % Use mex fn, faster but use linear interpolation
        opt.debug    = 1;
    end
    if nargin==0; syn=opt; return; end
    if nargin==1;
        syn = hmpd_synthesis();
        syn.enc = f0s;
        return;
    end

    if opt.debug>0; disp('HMPD Vocoder: Synthesis ...'); end
    if opt.debug>1; disp(opt); end

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

    if nargin<5 || isempty(fs);     fs=44100; end
    if nargin<6 || isempty(wavlen); wavlen=round(fs*f0s(end,1)); end

    unvoiced = f0s(:,2)==0; % Remember the unvoiced indices ...
    f0s = fillf0(f0s);      % ... and replace them with interpolated f0

    [AH, APH] = hmpd_amplitude_decompress(AE, f0s, fs, opt); % Amp decoding
    Hmax = size(AH,2)-1; % Max number of harmonics (without DC)

    % Reconstruct a Relative Phase Shift (RPS) of the voice source
    RPS = hmpd_phase_decompress(PDM, PDD, f0s, fs, opt); % Phase decoding

    % If asked, force full RPS randomization in unvoiced segments
    if opt.rps_randunvoiced
        if opt.debug>0; disp('    Randomize unvoiced segments (where f0=0)');end
        idx = find(unvoiced);
        RPS(idx,:) = wrap(2.*pi.*rand(length(idx),size(RPS,2)));
    end

    % Add the VTF phase
    RPS = wrap(opt.pd_gain*RPS + opt.amp_minphasegain*APH);

    if opt.debug>0; disp(['    HM synthesis (' num2str(Hmax) ' harmonics)']);end
    deflogamp = log(db2mag(opt.defdbamp));

    % Compute instantaneous fundamental frequency along the whole recording
    f1 = interp1(f0s(:,1), f0s(:,2), (0:wavlen-1)'/fs, 'linear');
    f1 = interp1_fixnan(f1)'; % Fix bounds
    p1 = filter(1, [1 -1], (2*pi/fs)*f1);

    syn = zeros(wavlen,1); % Allocate the synthetic signal
    % For each harmonic
    for h=1:Hmax % here h is the harmonic number

        % Synthesize only indices where the frequency doesn't go above Nyquist
        idxsyn = find(f1*h<fs/2);
        if length(idxsyn)==0; break; end
        synts = (idxsyn-1)/fs;

        % Amplitude ------------------------------------------------------------
        % Interpolate on a log scale
        Y = AH(:,1+h);
        if opt.usemex
            am = interp1ordered(f0s(:,1), Y, synts, deflogamp).';
        else
            am = interp1(f0s(:,1), Y, synts, 'linear', deflogamp);
        end

        % Phase ----------------------------------------------------------------
        % Interpolate on RPS domain
        Y = unwrap(RPS(:,1+h));
        if opt.usemex
            rpst = interp1ordered(f0s(:,1), Y, synts, 0).'; % Lose the splines !
        else
            rpst = interp1(f0s(:,1), Y, synts, 'spline');
        end
        % Instantaneous phase = linear phase + relative phase shift
        ph = h.*p1(idxsyn) + rpst;

        % Harmonic synthesis ---------------------------------------------------
        syn(idxsyn) = syn(idxsyn) + 2*exp(am).*cos(ph);
    end

return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [eta Heta inner_it stop_tCG storedb] ...
                      = tCG(problem, x, grad, eta, Delta, options, storedb)
% tCG - Truncated (Steihaug-Toint) Conjugate-Gradient method
% minimize <eta,grad> + .5*<eta,Hess(eta)>
% subject to <eta,eta> <= Delta^2

% This file is part of Manopt: www.manopt.org.
% This code is an adaptation to Manopt of the original GenRTR code:
% RTR - Riemannian Trust-Region
% (c) 2004-2007, P.-A. Absil, C. G. Baker, K. A. Gallivan
% Florida State University
% School of Computational Science
% (http://www.math.fsu.edu/~cbaker/GenRTR/?page=download)
% See accompanying license file.
% The adaptation was executed by Nicolas Boumal.
% Change log:
%   NB Feb. 12, 2013:
%       We do not project r back to the tangent space anymore: it was not
%       necessary, and as of Manopt 1.0.1, the proj operator does not
%       coincide with this notion anymore.
%   NB April 3, 2013:
%       tCG now also returns Heta, the Hessian at x along eta. Additional
%       esthetic modifications.
%   NB Dec. 2, 2013:
%       If options.useRand is activated, we now make sure the preconditio-
%       ner is not used, as was originally intended in GenRTR. In time, we
%       may want to investigate whether useRand can be modifed to work well
%       with preconditioning too.
%   NB Jan. 9, 2014:
%       Now checking explicitly for model decrease at each iteration. The
%       first iteration is a Cauchy point, which necessarily realizes a
%       decrease of the model cost. If a model increase is witnessed
%       (which is theoretically impossible if a linear operator is used for
%       the Hessian approximation), then we return the previous eta. This
%       ensures we always achieve at least the Cauchy decrease, which
%       should be sufficient for convergence.

% All terms involving the trust-region radius will use an inner product
% w.r.t. the preconditioner; this is because the iterates grow in
% length w.r.t. the preconditioner, guaranteeing that we will not
% re-enter the trust-region.
%
% The following recurrences for Prec-based norms and inner
% products come from [CGT2000], pg. 205, first edition.
% Below, P is the preconditioner.
%
% <eta_k,P*delta_k> =
%          beta_k-1 * ( <eta_k-1,P*delta_k-1> + alpha_k-1 |delta_k-1|^2_P )
% |delta_k|^2_P = <r_k,z_k> + beta_k-1^2 |delta_k-1|^2_P
%
% therefore, we need to keep track of
% 1)   |delta_k|^2_P
% 2)   <eta_k,P*delta_k> = <eta_k,delta_k>_P
% 3)   |eta_k  |^2_P
%
% initial values are given by:
%    |delta_0|_P = <r,z>
%    |eta_0|_P   = 0
%    <eta_0,delta_0>_P = 0
% because we take eta_0 = 0 (if useRand = false).
%
% [CGT2000] Conn, Gould and Toint: Trust-region methods, 2000.

inner = problem.M.inner;

theta = options.theta;
kappa = options.kappa;

if ~options.useRand % and therefore, eta == 0
    Heta = problem.M.zerovec(x);
    r = grad;
    e_Pe = 0;
else % and therefore, no preconditioner
    % eta (presumably) ~= 0 was provided by the caller
    [Heta storedb] = getHessian(problem, x, eta, storedb);
    r = problem.M.lincomb(x, 1, grad, 1, Heta);
    e_Pe = inner(x, eta, eta);
end
r_r = inner(x, r, r);
norm_r = sqrt(r_r);
norm_r0 = norm_r;

% precondition the residual
if ~options.useRand
    [z storedb] = getPrecon(problem, x, r, storedb);
else
    z = r;
end

% compute z'*r
z_r = inner(x, z, r);
d_Pd = z_r;

% Initial search direction
delta  = problem.M.lincomb(x, -1, z);
if ~options.useRand % and therefore, eta == 0
    e_Pd = 0;
else % and therefore, no preconditioner
    e_Pd = inner(x, eta, delta);
end

% If the Hessian or a linear Hessian approximation is in use, it is
% theoretically guaranteed that the model value decreases monotonically
% with each iteration of tCG. Hence, there is no need to monitor the model
% value. But, when a nonlinear Hessian approximation is used (such as the
% built-in finite-difference approximation for example), the model may
% increase. It is then important to terminate the tCG iterations and return
% the previous (the best-so-far) iterate. The variable below will hold the
% model value.
model_fun = @(eta, Heta) inner(x, grad, eta) + .5*inner(x, eta, Heta);
if ~options.useRand
    model_value = 0;
else
    model_value = model_fun(eta, Heta);
end

% Pre-assume termination b/c j == end
stop_tCG = 5;

% Begin inner/tCG loop
j = 0;
for j = 1 : options.maxinner

    [Hdelta storedb] = getHessian(problem, x, delta, storedb);

    % Compute curvature
    d_Hd = inner(x, delta, Hdelta);

    alpha = z_r/d_Hd;
    % <neweta,neweta>_P =
    % <eta,eta>_P + 2*alpha*<eta,delta>_P + alpha*alpha*<delta,delta>_P
    e_Pe_new = e_Pe + 2.0*alpha*e_Pd + alpha*alpha*d_Pd;

    if options.debug > 2,
        fprintf('DBG:   (r,r)  : %e\n',r_r);
        fprintf('DBG:   (d,Hd) : %e\n',d_Hd);
        fprintf('DBG:   alpha  : %e\n',alpha);
    end

    % Check against negative curvature and trust-region radius violation.
    % If either condition triggers, we bail out.
    if d_Hd <= 0 || e_Pe_new >= Delta^2,
        % want
        %  ee = <eta,eta>_prec,x
        %  ed = <eta,delta>_prec,x
        %  dd = <delta,delta>_prec,x
        tau = (-e_Pd + sqrt(e_Pd*e_Pd + d_Pd*(Delta^2-e_Pe))) / d_Pd;
        if options.debug > 2,
            fprintf('DBG:     tau  : %e\n', tau);
        end
        eta  = problem.M.lincomb(x, 1,  eta, tau,  delta);
        Heta = problem.M.lincomb(x, 1, Heta, tau, Hdelta);
        if d_Hd <= 0,
            stop_tCG = 1;     % negative curvature
        else
            stop_tCG = 2;     % exceeded trust region
        end
        break;
    end

    % No negative curvature and eta_prop inside TR: accept it
    e_Pe = e_Pe_new;
    new_eta  = problem.M.lincomb(x, 1,  eta, alpha,  delta);
    new_Heta = problem.M.lincomb(x, 1, Heta, alpha, Hdelta);

    % Verify that the model cost decreased in going from eta to new_eta. If
    % the cost increased (which can only occur if the Hessian approximation
    % is nonlinear), then we return the previous eta (which necessarily is
    % the best reached so far, according to the model cost). Otherwise, we
    % accept the new eta and go on.
    new_model_value = model_fun(new_eta, new_Heta);
    if new_model_value > model_value
        stop_tCG = 6;
        break;
    end

    eta = new_eta;
    Heta = new_Heta;

    % Update the residual
    r = problem.M.lincomb(x, 1, r, alpha, Hdelta);

    % Compute new norm of r
    r_r = inner(x, r, r);
    norm_r = sqrt(r_r);

    % Check kappa/theta stopping criterion
    if j >= options.mininner && norm_r <= norm_r0*min(norm_r0^theta, kappa)
        % Residual is small enough to quit
        if kappa < norm_r0^theta,
            stop_tCG = 3;  % linear convergence
        else
            stop_tCG = 4;  % superlinear convergence
        end
        break;
    end

    % Precondition the residual
    if ~options.useRand
        [z storedb] = getPrecon(problem, x, r, storedb);
    else
        z = r;
    end

    % Save the old z'*r
    zold_rold = z_r;
    % Compute new z'*r
    z_r = inner(x, z, r);

    % Compute new search direction
    beta = z_r/zold_rold;
    delta = problem.M.lincomb(x, -1, z, beta, delta);

    % Update new P-norms and P-dots [CGT2000, eq. 7.5.6 & 7.5.7]
    e_Pd = beta*(e_Pd + alpha*d_Pd);
    d_Pd = z_r + beta*beta*d_Pd;

end  % of tCG loop
inner_it = j;

end

author: andrewssobral

repo: lrslibrary


text:
function Mutual_Information=MutualInformation(x,y)

% Compute Mutual Information and Joint entropy of two images
% Takes two images and returns the mutual information and  joint
% entropy. For implementing this function also download function
% joint_histogram.m in my file exchange.
%
% written by Amir Pasha Mahmoudzadeh
% Wright State University
% Biomedical Imaging Lab

h=joint_histogram(x,y); % Joint histogram for 2D
[m,n] = size(h);
J= h./(m*n); % normalized joint histogram
y_summ=sum(J); %sum of the rows of normalized joint histogram
x_summ=sum(J);%sum of columns of normalized joint histogran

sy=0;
for s=1:n;    %  column
      if( y_summ(s)==0 )

      else
         sy = sy + -(y_summ(s)*(log2(y_summ(s)))); %marginal entropy for image 1
      end
   end

sx=0;
for s=1:m;    %rows
   if( x_summ(s)==0 )

      else
         sx = sx + -(x_summ(s)*(log2(x_summ(s)))); %marginal entropy for image 2
      end
   end
joint_entropy = -sum(sum(J.*(log2(J+(J==0))))) % joint entropy
Mutual_Information = sx + sy - joint_entropy; % Mutual information
end


text:
% ------------------------------------------------------------------------
%          Run Allen Atlas Browser
% ------------------------------------------------------------------------

%% ENTER FILE LOCATION AND PROBE-SAVE-NAME

% directory of histology
processed_images_folder = 'C:\Drive\Histology\brainX\processed';

% name the saved probe points, to avoid overwriting another set of probes going in the same folder
probe_save_name_suffix = '';

% directory of reference atlas files
annotation_volume_location = 'C:\Drive\Histology\annotation_volume_10um_by_index.npy';
structure_tree_location = 'C:\Drive\Histology\structure_tree_safe_2017.csv';
template_volume_location = 'C:\Drive\Histology\template_volume_10um.npy';

% plane to view ('coronal', 'sagittal', 'transverse')
plane = 'coronal';

%% GET PROBE TRAJECTORY POINTS

% load the reference brain and region annotations
if ~exist('av','var') || ~exist('st','var') || ~exist('tv','var')
    disp('loading reference atlas...')
    av = readNPY(annotation_volume_location);
    st = loadStructureTree(structure_tree_location);
    tv = readNPY(template_volume_location);
end

% select the plane for the viewer
if strcmp(plane,'coronal')
    av_plot = av;
    tv_plot = tv;
elseif strcmp(plane,'sagittal')
    av_plot = permute(av,[3 2 1]);
    tv_plot = permute(tv,[3 2 1]);
elseif strcmp(plane,'transverse')
    av_plot = permute(av,[2 3 1]);
    tv_plot = permute(tv,[2 3 1]);
end

% create Atlas viewer figure
f = figure('Name','Atlas Viewer');

% show histology in Slice Viewer
try; figure(slice_figure_browser); title('');
catch; slice_figure_browser = figure('Name','Slice Viewer'); end
reference_size = size(tv_plot);
sliceBrowser(slice_figure_browser, processed_images_folder, f, reference_size);

% % use application in Atlas Transform Viewer
% % use this function if you have a processed_images_folder with appropriately processed .tif histology images
f = AtlasTransformBrowser(f, tv_plot, av_plot, st, slice_figure_browser, processed_images_folder, probe_save_name_suffix, plane);

% use the simpler version, which does not interface with processed slice images
% just run these two lines instead of the previous 5 lines of code
%
%  save_location = processed_images_folder;
%  f = allenAtlasBrowser(f, tv_plot, av_plot, st, save_location, probe_save_name_suffix, plane);

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
function h = ksizeHall(npd)
%
% Find kernel size according to "plug-in" method of
%     Hall, Marron, Sheather, Jones (91)
%
%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

  x = getPoints(npd);
  [N1,N2] = size(x);
  sig = std(x,0,2);                     % estimate sigma (standard)
  lamS= .7413 * iqr(x')';               % find sigma by interquartile range lam
  if (max(lamS)==0) lamS=sig; end;      % replace sigma est. if possible
  BW = 1.0592 * lamS * N2^(-1/(4+N1));
  BW = repmat(BW,[1,N2]);

  dX = repmat(permute(x,[1,3,2]),[1,N2,1]);  % compute Xi-Xj for all i,j
  for i=1:N2,
    dX(:,:,i) = (dX(:,:,i)-x)./BW;
  end;
  for i=1:N2, dX(:,i,i) = 2e22; end;
  dX = reshape(dX,[N1,N2*N2]);

%  use that to find I2 and I3
  I2=h_findI2(N2,dX,BW(:,1));      % I2 = \hat R( f^(2) ) = \int f^(2)^2 dx
  I3=h_findI3(N2,dX,BW(:,1));      % I3 = \hat R( f^(3) )

% for Gaussian Kernel, we evaluate to find:
%  RK = 1.0/(2^N1) * 1.0/pi^(N1/2.0);    % R(K) = \int K^2(x) dx
%  mu2 = 1.0;                            % \mu_i = \int x^i K(x) dx
%  mu4 = 3.0^N1;
  switch (npd.type)
    case 0, RK = 0.282095;   mu2 = 1.000000;   mu4 = 3.000000; % Gauss
    case 1, RK = 0.600000;   mu2 = 0.199994;   mu4 = 0.085708; % Epanetch
    case 2, RK = 0.250002;   mu2 = 1.994473;   mu4 = 23.299070;% Laplace
  end;

  J1 = RK/mu2^2 .* 1./I2;
  J2 = (mu4 * I3) ./ (20 * mu2) .* 1./I2;
  h  = (J1/N2).^(1.0/5) + J2.*(J1/N2).^(3.0/5);

% Let us estimate R(f^(p)) by R( \hat f^(p) )
%   (f is the original density to be est'd;  f^(p) is its pth derivative)
%   Let L be the kernel function for this second estimator, with bandwidth alpha
%
% Ip = \int f^(p)^2_\alpha(x) dx
%    = [(-1)^p/n^2] \sum_i \sum_j L^(p)_\alpha * L^(p)_\alpha
%    = [(-1)^p/(n^2 \alpha^(2p+1))]  \sum_i \sum_j (L^(p) * L^(p))( (Xi-Xj)/\alpha )
%
% Take L to be a Gaussian kernel; we then evaluate L^(p) by:
%
  % L^(p)(x) = (-1)^p H_p(x) L(x)
  % H_p(x) = x H_{p-1}(x) - (p-1)H_{p-2}(x)
  % p=2 =>
  %   H_2(x) = x H_1(x) - H_0(x) = x * x - 1
  %   =>  L^(2)(x) = (x^2 - 1) * L(x)
  % p=3 =>
  %   H_3(x) = x H_2(x) - 2*H_1(x) = x*(x^2-1) - 2*x
  %   =>  L^(3)(x) = (x^3 - 3x) * L(x)
  %

function I2 = h_findI2(n,dXa,alpha)
%%  load ksizeHSJM.mat;
%%  xInd = fix(Nquant * (dXa-Xmin) / (Xmax-Xmin));
%%  xInd = max(xInd,1); xInd = min(xInd,Nquant);
%%  s = sum( L2data(xInd) ,2);
%  s = sum( (dXa.^2 -1) .* 1/sqrt(2*pi) .* exp(-.5*dXa.^2) , 2);
  s = sum( (dXa.^2 -1) .* 1/sqrt(2*pi) .* repmat(exp(-.5*sum(dXa.^2,1)),[size(dXa,1),1]) , 2);
  I2=s./((n*(n-1))*alpha.^5);

function I3 = h_findI3(n,dXb,beta)
%%  load ksizeHSJM.mat;
%%  xInd = fix(Nquant * (dXb-Xmin) / (Xmax-Xmin));
%%  xInd = max(xInd,1); xInd = min(xInd,Nquant);
%%  s = sum( L3data(xInd) , 2);
%  s = sum( (dXb.^3 -3*dXb) .* 1/sqrt(2*pi) .* exp(-.5*dXb.^2) , 2);
  s = sum( (dXb.^3 -3*dXb) .* 1/sqrt(2*pi) .* repmat(exp(-.5*sum(dXb.^2,1)),[size(dXb,1),1]) , 2);
  I3  = -s./((n*(n-1))*beta.^7);

author: ShapeNet

repo: RenderForCNN


text:
function resized = resize_img(img)
    s = size(img);
    % scaling up to 600 if necessary
    minside = min([s(1), s(2)]);
    if minside < 600
        scale = 600/minside;
        numrow = floor(size(img, 1) * scale);
        numcol = floor(size(img, 2) * scale);
        resized= imresize(img, [numrow, numcol]);
    end
    % scaling down to 1000 if necessary
    maxside = max([s(1), s(2)]);
    if maxside > 1000
        scale = 1000/maxside;
        numrow = floor(size(img, 1) * scale);
        numcol = floor(size(img, 2) * scale);
        resized = imresize(img, [numrow, numcol]);
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function test_suite = test_independent_samples_partitioner
% tests for cosmo_independent_samples_partitioner
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_independent_samples_partitioner_tiny()
    min_class_count=2;
    max_class_count=2+min_class_count;
    rng_class_count=[min_class_count,max_class_count];

    for nclasses=[2 3]
        for test_count=0:(nclasses+1)
            seed=ceil(rand()*1e6);
            opt=struct();
            opt.test_count=test_count;
            p1=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,opt,seed);

            opt=struct();
            opt.test_count=-test_count; % use test_ratio
            p2=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,opt,seed);

            assertEqual(p1,p2);
        end
    end

function test_independent_samples_partitioner_big()
    nclasses=ceil(rand()*4+10);
    for test_count=0:(nclasses+1)
        min_class_count=ceil(rand()*10+2);
        max_class_count=min_class_count+ceil(rand()*10+2);
        rng_class_count=[min_class_count,max_class_count];

        opt=struct();
        opt.test_count=test_count;
        opt.fold_count=ceil(rand()*100+10);

        p1=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt);

        opt=struct();
        opt.test_count=-test_count; % use test_ratio
        opt.fold_count=ceil(rand()*100+10);
        p2=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt);
    end

function test_independent_samples_partitioner_big_with_seed()
    opt=struct();
    opt.fold_count=100;
    opt.test_count=1;

    nclasses=ceil(rand()*4+10);
    min_class_count=ceil(rand()*10+2);
    max_class_count=min_class_count+ceil(rand()*10+2);
    rng_class_count=[min_class_count,max_class_count];

    ds_seed=(ceil(rand()*1e6));

    % without seed they use the default seed, must be equal
    p1=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);
    p2=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);
    assertEqual(p1,p2);
    opt.seed=1;
    p3=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);
    assertEqual(p1,p3);
    opt.seed=2;
    p4=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);
    assertFalse(isequal(p1,p4));

    % with no seed they must be identical
    opt.seed=0;
    p1=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);
    p2=helper_test_independent_samples_partitioner(nclasses,...
                                rng_class_count,...
                                opt,ds_seed);

    assertFalse(isequal(p1,p2));

function p=helper_test_independent_samples_partitioner(nclasses,...
                            rng_class_count,opt,gen_seed)
    if nargin<4
        gen_seed=0;
    end

    ds=helper_generate_dataset(nclasses,...
                            rng_class_count,...
                            gen_seed);

    % compute how many samples per class in test set
    class_counts=histc(ds.sa.targets,1:nclasses);
    min_class_count=min(class_counts);
    max_class_count=max(class_counts);

    test_count=opt.test_count;
    if test_count<0
        % use test_ratio
        test_count=-test_count;
        opt=rmfield(opt,'test_count');
        opt.test_ratio=test_count/min_class_count;

    end

    train_count=min_class_count-test_count;

    has_illegal_count=train_count<=0 || test_count<=0;
    if has_illegal_count
        if ~isfield(opt,'fold_count')
            opt.fold_count=1;
        end
    else
        if isfield(opt,'fold_count')
            % given by calling funciton
            fold_count=opt.fold_count;
        else
            % use all folds available
            assert(nclasses<=3,'fold_count required with nclasses>3');
            assert(max_class_count<=5,'too many classes');
            % When not set we generate all possible folds; thereare at most
            % nchoosek(5,3)^nreps <= 20^3 = 8000 possible folds
            max_fold_count=8000; %

            % see how many possible folds based on each class
            combi_test=@(i)nchoosek(class_counts(i),test_count);
            test_fold_counts=arrayfun(combi_test,1:nclasses);
            combi_train=@(i)nchoosek(class_counts(i)-test_count,train_count);
            train_fold_counts=arrayfun(combi_train,1:nclasses);

            % take product
            fold_count=prod(test_fold_counts)*prod(train_fold_counts);
            assert(fold_count<=max_fold_count,'memory safety limit exceeded');
            opt.fold_count=fold_count;
        end
    end

    func=@()cosmo_independent_samples_partitioner(ds,opt);

    if has_illegal_count
        % not enough samples, expect an exception
        assertExceptionThrown(func,'')
        p=[];
        return
    end

    p=func();

    assertEqual(numel(p.train_indices),fold_count);
    assertEqual(numel(p.test_indices),fold_count);

    nsamples=size(ds.samples,1);

    train_count=min(class_counts)-test_count;

    for f=1:fold_count
        tr=p.train_indices{f};
        te=p.test_indices{f};
        assertTrue(isempty(intersect(tr,te)));

        assert_all_int_less_than(tr,nsamples);
        assert_all_int_less_than(te,nsamples);

        % equal number of targets in all classes, for train and test
        assert_all_hist_equal(ds.sa.targets(te),nclasses,test_count);
        assert_all_hist_equal(ds.sa.targets(tr),nclasses,train_count);
    end

    assert_all_folds_unique(p.train_indices,p.test_indices)

function ds=helper_generate_dataset(nclasses,rng_class_count,seed)
    % generate dataset where each target occurs at least min_class_count
    % times
    min_count=rng_class_count(1);
    max_count=rng_class_count(2);

    seed_arg={'seed',seed};

    delta=(max_count-min_count);
    assert(delta>0);

    % make lots of trials
    nregular=nclasses*min_count;
    nextra=ceil(delta*nclasses*cosmo_rand(seed_arg{:}));

    nsamples=nregular+nextra;

    ds=struct();
    ds.samples=rand(nsamples,2);
    ds.sa.targets=zeros(nsamples,1);
    rp=cosmo_randperm(nsamples,seed_arg{:});
    ds.sa.targets(1:nsamples)=mod(rp,nclasses)+1;
    ds.sa.chunks=(1:nsamples)';

    h=histc(ds.sa.targets,1:nclasses);
    assert(min(h)>=min_count)
    assert(max(h)<=max_count)

function assert_all_folds_unique(xs,ys)
    assert(all(min(cellfun(@min,xs))>=0));
    assert(all(min(cellfun(@min,ys))>=0));

    xs_max=max(cellfun(@max,xs));
    ys_max=max(cellfun(@max,ys));

    % value greater than max in all
    ys_mark=1+max(xs_max,ys_max);

    nfolds=numel(xs);
    merged=cell(nfolds,1);
    for f_i=1:nfolds

        xy=sort([xs{f_i}(:); (ys_mark+ys{f_i}(:))]);
        merged{f_i}=xy(:)';
    end

    % must all be same length
    c=cellfun(@numel,merged);
    assertEqual(c(1)+zeros(size(c)),c,'inputs do not have same length');

    % put in matrix
    merged_mat=cat(1,merged{:});
    s=sortrows(merged_mat);

    % look for duplicate rows
    eq_msk=bsxfun(@eq,s(1:(end-1),:),s(2:end,:));
    row_same=find(all(eq_msk,2),1);

    assertEqual(row_same,zeros(0,1),'row duplicate');

function assert_all_hist_equal(targets,nclasses,nreps)
    h_targets=histc(targets(:)',1:nclasses);
    assertEqual(h_targets,nreps+zeros(1,nclasses));

function assert_all_int_less_than(x,mx)
    assert(isnumeric(x));
    assertEqual(sort(x),x);
    assert(all(x>=1));
    assert(all(x<=mx));
    assert(all(round(x)==x));

function test_independent_samples_partitioner_mismatch_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                   cosmo_independent_samples_partitioner(varargin{:}),'');

    aet_targets=@(counts,varargin)aet(...
                 helper_generate_dataset_with_target_counts(counts{:}),...
                                                            varargin{:});

    opt=struct();
    opt.test_count=1;
    opt.fold_count=1;

    % missing target
    aet_targets({[3 3 3],[4 4]},opt);

    % not enought targets in one class
    aet_targets({[3 3 3],[2 2 1]},opt);

    opt.test_count=3;
    aet_targets({[4 3 4],[4 4 4]},opt);

    % try with ratio
    opt=rmfield(opt,'test_count');

    % missing target
    opt.test_ratio=.25;
    aet_targets({[10 10 10],[10 10]},opt);

    % too few targets
    aet_targets({[2 2 2],[4 4 4]},opt);

    % with too many folds
    aet_targets({[4 4 4],[4 4 4]},opt,'max_fold_count',0);

function test_independent_samples_partitioner_arg_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                   cosmo_independent_samples_partitioner(varargin{:}),'');
    ds=cosmo_synthetic_dataset();
    ds.sa.chunks=(1:size(ds.samples,1))';

    % missing arguments
    aet(ds);
    aet(ds,'fold_count');
    aet(ds,'fold_count',2);

    % mutually exclusive arguments
    aet(ds,'fold_count',2,'test_count',1,'test_ratio',.5);

    % not a dataset
    aet(struct,'fold_count',2,'test_count',1);

    % non-unique
    ds_bad=ds;
    ds_bad.sa.chunks(2)=ds.sa.chunks(1);
    aet(ds_bad,'fold_count',2,'test_count',1);

function ds=helper_generate_dataset_with_target_counts(varargin)
    nfeatures=2;
    nchunks=numel(varargin);
    ds_cell=cell(nchunks,1);
    for i=1:nchunks
        counts=varargin{i};
        nclasses=numel(counts);
        ds_parts=cell(nclasses,1);
        for j=1:nclasses
            nt=counts(j);
            ds=struct();

            ds.samples=randn(nt,nfeatures);
            ds.sa.targets=zeros(nt,1)+j;
            ds.sa.chunks=zeros(nt,1)+i;
            ds_parts{j}=ds;
        end

        ds_cell{i}=cosmo_stack(ds_parts);
    end

    ds=cosmo_stack(ds_cell);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function determ = gk324_determinant ( n, x )

%*****************************************************************************80
%
%% GK324_DETERMINANT returns the determinant of the GK324 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N-1), the first N-1 entries of the
%    last row.
%
%    Output, real DETERM, the determinant.
%
  determ = prod ( 1.0 - x(1:n-1) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = english_word_length_cdf_inv ( cdf )

%*****************************************************************************80
%
%% ENGLISH_WORD_LENGTH_CDF_INV inverts the English Word Length CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Henry Kucera, Winthrop Francis,
%    Computational Analysis of Present-Day American English,
%    Brown University Press, 1967.
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%    0.0 <= CDF <= 1.0.
%
%    Output, integer X, the corresponding word length for which
%    CDF(X-1) < CDF <= CDF(X)
%
  word_length_max = 27;

  pdf_vec = [ ...
    0.03160, ...
    0.16975, ...
    0.21192, ...
    0.15678, ...
    0.10852, ...
    0.08524, ...
    0.07724, ...
    0.05623, ...
    0.04032, ...
    0.02766, ...
    0.01582, ...
    0.00917, ...
    0.00483, ...
    0.00262, ...
    0.00099, ...
    0.00050, ...
    0.00027, ...
    0.00022, ...
    0.00011, ...
    0.00006, ...
    0.00005, ...
    0.00002, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001, ...
    0.00001 ];
  pdf_sum = 0.99997;

  if ( cdf < 0.0 || 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ENGLISH_WORD_LENGTH_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'ENGLISH_WORD_LENGTH_CDF_INV - Fatal error!' );
  end

  cum = 0.0;

  for j = 1 : word_length_max

    cum = cum + pdf_vec(j);

    if ( cdf <= cum / pdf_sum )
      x = j;
      return
    end

  end

  x = word_length_max;

  return
end


text:
x=zeros(2,100);
x(1,:)=linspace(-3,3,100);
x(2,:)=linspace(-3,3,100);

[X1,X2]=meshgrid(x(1,:),x(2,:));
Z=zeros(100,100);
for i=1:100
    for j=1:100
        Z(i,j)=obj2([X1(i,j),X2(i,j)]);
    end
end

plot3(X1,X2,Z)
% surf(X1,X2,Z)
% contour(X1,X2,Z)

repo: Optimization-Theory-and-Algorithm


text:
function [ sa, seed ] = wshrt ( d, n, np, seed )

%*****************************************************************************80
%
%% WSHRT returns a random Wishart variate.
%
%  Discussion:
%
%    This routine is a Wishart variate generator.
%
%    On output, SA is an upper-triangular matrix of size NP * NP,
%    written in linear form, column ordered, whose elements have a
%    Wishart(N, SIGMA) distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2014
%
%  Author:
%
%    Original FORTRAN77 version by William Smith, Ronald Hocking.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Smith, Ronald Hocking,
%    Algorithm AS 53, Wishart Variate Generator,
%    Applied Statistics,
%    Volume 21, Number 3, pages 341-345, 1972.
%
%  Parameters:
%
%    Input, real D(NP*(NP+1)/2), the upper triangular array that
%    represents the Cholesky factor of the correlation matrix SIGMA.
%    D is stored in column-major form.
%
%    Input, integer N, the number of degrees of freedom.
%    1 <= N <= NP.
%
%    Input, integer NP, the size of variables.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real SA(NP*(NP+1)/2), a sample from the Wishart distribution.
%
  k = 1;
%
%  Load SB with independent normal (0, 1) variates.
%
  nnp = floor ( ( np * ( np + 1 ) ) / 2 );
  sb = zeros ( nnp );

  while ( k <= nnp )

    [ u1, u2, seed ] = rnorm ( seed );

    sb(k) = u1;
    k = k + 1;

    if ( k <= nnp )
      sb(k) = u2;
      k = k + 1;
    end

  end
%
%  Load diagonal elements with square root of chi-square variates.
%
  sa = zeros ( nnp );

  ns = 0;

  for i = 1 : np
%
%  The original text read "DF = N - I + 1".
%  This should read "DF = NP - I + 1".
%
    df = np - i + 1;
    ns = ns + i;
    u1 = 2.0 / ( 9.0 * df );
    u2 = 1.0 - u1;
    u1 = sqrt ( u1 );
%
%  Wilson-Hilferty formula for approximating chi-square variates:
%  The original code did not take the absolute value!
%
    sb(ns) = sqrt ( df * abs ( ( u2 + sb(ns) * u1 ) ^ 3 ) );

  end

  rn = n;
  nr = 1;

  for i = 1 : np
    nr = nr + i - 1;
    for j = i : np
      ip = nr;
      nq = floor ( ( j * ( j - 1 ) ) / 2 ) + i - 1;
      c = 0.0;
      for k = i : j
        ip = ip + k - 1;
        nq = nq + 1;
        c = c + sb(ip) * d(nq);
      end
      sa(ip) = c;
    end
  end

  for i = 1 : np
    ii = np - i + 1;
    nq = nnp - np;
    for j = 1 : i
      ip = floor ( ( ii * ( ii - 1 ) ) / 2 );
      c = 0.0;
      for k = i : np
        ip = ip + 1;
        nq = nq + 1;
        c = c + sa(ip) * sa(nq);
      end
      sa(nq) = c / rn;
      nq = nq - 2 * np + i + j - 1;
    end
  end

  return
end


text:
function [Xrls,yrls,s]=RLSreg(y,X)
%Syntax: [Xrls,yrls,s]=RLSreg(y,X)
%_______________________________
%
% Calculates the Reweighted Least Squares (RLS) regression data points
% and the scale parameter from the LMS regression.
%
% Xrls is the X values matrix to be taken into account for RLS.
% yrls is the y values vector to be taken into account for RLS.
% s is RLS scale parameter.
% y is the vector of the dependent variable.
% X is the data matrix of the independent variable.
%
% Reference:
% Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
%
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
%
% University e-mail: al10@ukc.ac.uk (until December 2002)
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Sep 3, 2001.

if nargin<1 | isempty(y)==1
   error('Not enough input arguments.');
end

% Estimate the LMS values
if nargin<2 | isempty(X)==1
   LMSout=LMSreg(y);
   X=(1:length(y))';
else
   LMSout=LMSreg(y,X);
end

% p is the number of parameters to be estimated
p=size(X,2)+1;

% Calculate the residuals
r=y-LMSout;

% Estimate the preliminary scale parameter
s=LMSsca(r,0,p);

% Take into account a data point, if its residual is relatively small
w=find(abs(r/s)<=2.5);
Xrls=X(w,:);
yrls=y(w);

author: Sable

repo: mcbench-benchmarks