text:
function xPosVel=makeConstSpeedLevelCurvedEarthTraj(llhStart,headingStart,speed,timeOffsets,useRhumb,trajParams)
%%MAKECONSTSPEEDLEVELCURVEDEARTHTRAJ Create a non-maneuvering, constant
%   (local tangent plane) speed, constant ellipsoidal height trajectory
%   over an ellipsoidal Earth, with position and velocity sampled at
%   specified times. By default, this is a geodesic path, but a rhumb
%   trajectory (constant heading) can also be used. Note that rhumb lines
%   will get stuck at the geographic poles.
%
%INPUTS: llhStart The [latitude;longitude;ellipsoidal height] of the
%                 starting point (time=0) of the trajectory.
%    headingStart The initial heading of the trajectory in radians East of
%                 North.
%           speed The positive speed of the target. This will typically be
%                 in meters.
%     timeOffsets A length numPts list of times where the position and
%                 velocity of the target should be obtained. This will
%                 typically be in seconds.
%        useRhumb If this is true, a rhumb (constant heading) trajectory is
%                 created as opposed to a geodesic (straightest)
%                 trajectory. The default if omitted or an empty matrix is
%                 passed is false.
%      trajParams An optional structure containing values that change how
%                 the trajectory is generated. Possible field and their
%                 defaults (if omitted or empty matrices are passed) are:
%                 a The semi-major axis of the reference ellipsoid. The
%                   default is Constants.WGS84SemiMajorAxis.
%                 f The flattening factor of the reference ellipsoid. The
%                   default is Constants.WGS84Flattening.
%                 useHeightApprox If true, an approximation for targets
%                   with nonzero height is used that is much faster than a
%                   numerical integration-based technique. The default is
%                   true.
%                 numSteps4Circ If useHeightApprox is false, then numerical
%                   integration is performed and this is the number of
%                   Runge-Kutta steps that would be taken when
%                   circumnavigating the equator. The default is 2000,
%                   though 6000 is often better.
%
%OUTPUTS: xPosVel The 6XnumPts set of [position;velocity] vectors of the
%                 target at the specified times.
%
%This function obtains positions using directRhumbProbGen or
%directGeodeticProbGen and then converts the local heading to ECEF, which
%when given magnitude speed is an approximation to the instantaneous
%velocity vector.
%
%October 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(useRhumb))
    useRhumb=false;
end

useHeightApprox=true;
numSteps4Circ=2000;
a=Constants.WGS84SemiMajorAxis;
f=Constants.WGS84Flattening;

if(nargin>5&&~isempty(trajParams))
    if(isfield(trajParams,'a'))
        a=trajParams.a;
    end

    if(isfield(trajParams,'f'))
        f=trajParams.f;
    end

    if(isfield(trajParams,'useHeightApprox'))
        useHeightApprox=trajParams.useHeightApprox;
    end

    if(isfield(trajParams,'numSteps4Circ'))
        numSteps4Circ=trajParams.numSteps4Circ;
    end
end

numPts=length(timeOffsets);
height=llhStart(3);
latLonStart=llhStart(1:2);

xPosVel=zeros(6,numPts);
if(useRhumb)
    %The instantaneous local ENU velocity never changes, so find it outside
    %the loop.
    vENULocal=geogHeading2ENUVec(headingStart,0)*speed;

    for curOffset=1:numPts
        curDist=timeOffsets(curOffset)*speed;
        latLonEnd=directRhumbProbGen(latLonStart,headingStart,curDist,height,useHeightApprox,a,f,numSteps4Circ);
        vECEF=ENUVec2ECEFVec(latLonEnd,vENULocal,a,f);
        posCart=ellips2Cart([latLonEnd;height],a,f);
        xPosVel(:,curOffset)=[posCart;vECEF];
    end
else
    %A geodesic trajectory.
    for curOffset=1:numPts
        curDist=timeOffsets(curOffset)*speed;
        [latLonEnd,azEnd]=directGeodeticProbGen(latLonStart,headingStart,curDist,height,useHeightApprox,a,f,numSteps4Circ);

        vENULocal=geogHeading2ENUVec(azEnd,0)*speed;
        vECEF=ENUVec2ECEFVec(latLonEnd,vENULocal,a,f);
        posCart=ellips2Cart([latLonEnd;height],a,f);
        xPosVel(:,curOffset)=[posCart;vECEF];
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function rez  = template_learning(rez, tF, st3)

wPCA  = rez.wPCA;
iC = rez.iC;
ops = rez.ops;

xcup = rez.xcup;
ycup = rez.ycup;

wroll = [];
tlag = [-2, -1, 1, 2];
for j = 1:length(tlag)
    wroll(:,:,j) = circshift(wPCA, tlag(j), 1)' * wPCA;
end

%% split templates into batches
rmin = 0.6;
nlow = 100;
n0 = 0;
use_CCG = 0;

Nchan = rez.ops.Nchan;
Nk = size(iC,2);
yunq = unique(rez.yc);

ktid = int32(st3(:,2)) + 1;
tmp_chan = iC(1, :);
ss = double(st3(:,1)) / ops.fs;

dmin = rez.ops.dmin;
ycenter = (min(rez.yc) + dmin-1):(2*dmin):(max(rez.yc)+dmin+1);
dminx = rez.ops.dminx;
xcenter = (min(rez.xc) + dminx-1):(2*dminx):(max(rez.xc)+dminx+1);
[xcenter, ycenter] = meshgrid(xcenter, ycenter);
xcenter = xcenter(:);
ycenter = ycenter(:);

Wpca = zeros(6, Nchan, 1000, 'single');
nst = numel(ktid);
hid = zeros(nst,1 , 'int32');

% ycup = rez.yc;

tic

for j = 1:numel(ycenter)
    if rem(j,5)==1
        fprintf('time %2.2f, GROUP %d/%d, units %d \n', toc, j, numel(ycenter), n0)
    end

    y0 = ycenter(j);
    x0 = xcenter(j);
    xchan = (abs(ycup - y0) < dmin) & (abs(xcup - x0) < dminx);
    itemp = find(xchan);

    if isempty(itemp)
        continue;
    end
    tin = ismember(ktid, itemp);
    pid = ktid(tin);
    data = tF(tin, :, :);

    if isempty(data)
        continue;
    end
%     size(data)

    ich = unique(iC(:, itemp));
%     ch_min = ich(1)-1;
%     ch_max = ich(end);

    if numel(ich)<1
        continue;
    end

    nsp = size(data,1);
    dd = zeros(nsp, 6, numel(ich), 'single');
    for k = 1:length(itemp)
        ix = pid==itemp(k);
        % how to go from channels to different order
        [~,ia,ib] = intersect(iC(:,itemp(k)), ich);
        dd(ix, :, ib) = data(ix,:,ia);
    end

    kid = run_pursuit(dd, nlow, rmin, n0, wroll, ss(tin), use_CCG);

    [~, ~, kid] = unique(kid);
    nmax = max(kid);
    for t = 1:nmax
%         Wpca(:, ch_min+1:ch_max, t + n0) = gather(sq(mean(dd(kid==t,:,:),1)));
        Wpca(:, ich, t + n0) = gather(sq(mean(dd(kid==t,:,:),1)));
    end

    hid(tin) = gather(kid + n0);
    n0 = n0 + nmax;
end
Wpca = Wpca(:,:,1:n0);
toc
%%
rez.W = zeros(61,0, 3, 'single');
rez.U = zeros(ops.Nchan,0,3, 'single');
rez.mu = zeros(1,0, 'single');
for  t = 1:n0
    dWU = wPCA * gpuArray(Wpca(:,:,t));
    [w,s,u] = svdecon(dWU);
    wsign = -sign(w(21,1));
    rez.W(:,t,:) = gather(wsign * w(:,1:3));
    rez.U(:,t,:) = gather(wsign * u(:,1:3) * s(1:3,1:3));
    rez.mu(t) = gather(sum(sum(rez.U(:,t,:).^2))^.5);
    rez.U(:,t,:) = rez.U(:,t,:) / rez.mu(t);
end

%%
rez.ops.wPCA = wPCA;

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
%% Demo to show the results of MCG
clear all;close all;home;

% Read an input image
I = imread(fullfile(root_dir, 'demos','101087.jpg'));

tic;
% Test the 'fast' version, which takes around 5 seconds in mean
[candidates_scg, ucm2_scg] = im2mcg(I,'fast');
toc;

tic;
% Test the 'accurate' version, which tackes around 30 seconds in mean
[candidates_mcg, ucm2_mcg] = im2mcg(I,'accurate');
toc;

%% Show UCM results (dilated for visualization)
figure;
subplot(1,3,1)
imshow(I), title('Image')

subplot(1,3,2)
imshow(imdilate(ucm2_scg,strel(ones(3))),[]), title('Fast UCM (SCG)')

subplot(1,3,3)
imshow(imdilate(ucm2_mcg,strel(ones(3))),[]), title('Accurate UCM (MCG)')

%% Show Object Candidates results and bounding boxes
% Candidates in rank position 11 and 12
id1 = 11; id2 = 12;

% Get the masks from superpixels and labels
mask1 = ismember(candidates_mcg.superpixels, candidates_mcg.labels{id1});
mask2 = ismember(candidates_mcg.superpixels, candidates_mcg.labels{id2});

% Bboxes is a matrix that contains the four coordinates of the bounding box
% of each candidate in the form [up,left,down,right]. See folder bboxes for
% more function to work with them

% Show results
figure;
subplot(1,3,1)
imshow(I), title('Image')
subplot(1,3,2)
imshow(mask1), title('Candidate + Box')
hold on
plot([candidates_mcg.bboxes(id1,4) candidates_mcg.bboxes(id1,4) candidates_mcg.bboxes(id1,2) candidates_mcg.bboxes(id1,2) candidates_mcg.bboxes(id1,4)],...
     [candidates_mcg.bboxes(id1,3) candidates_mcg.bboxes(id1,1) candidates_mcg.bboxes(id1,1) candidates_mcg.bboxes(id1,3) candidates_mcg.bboxes(id1,3)],'r-')
subplot(1,3,3)
imshow(mask2), title('Candidate + Box')
hold on
plot([candidates_mcg.bboxes(id2,4) candidates_mcg.bboxes(id2,4) candidates_mcg.bboxes(id2,2) candidates_mcg.bboxes(id2,2) candidates_mcg.bboxes(id2,4)],...
     [candidates_mcg.bboxes(id2,3) candidates_mcg.bboxes(id2,1) candidates_mcg.bboxes(id2,1) candidates_mcg.bboxes(id2,3) candidates_mcg.bboxes(id2,3)],'r-')

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
%%% void show_4dvideo(4dmatrix)
%
function show_4dvideo(video)
  for i = 1 : size(video,4)
    frame = video(:,:,:,i);
    imshow(frame,[],'InitialMagnification','fit');
    disp(i);
    pause(0.01);
  end
end

author: andrewssobral

repo: lrslibrary


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: binopdf_01ok.m
% The regular binopdf returns NaN whenever x=0 and p=0 or whenever
%   x=n and p=1.  The correct value for binopdf to return in both cases is 1.0.
% This corrects that bug.
% This also avoids calling binopdf when it would give an inappropriate warning
%   about taking the log of zero.
% This could more easily be fixed within binopdf.m, but that for MathWorks
%
% 022400 tdr created
% 033100 tdr added check for binopdf returning NaN or Inf, which it may for n > 1000, e.g., x=500, n=1030, p=0.1

function y = binopdf_01ok(x,n,p)

if nargin < 3,
    error('Requires three input arguments (x,n,p)');
end

[errorcode x n_bar p] = distchck(3,x,n,p);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

y = zeros(size(x));

k = find(~((p == 0) | (p == 1)));
if any(k),
   y(k) = binopdf(x(k),n,p(k));
   if ~isfinite(y),
      error('NaN or Inf value returned from binopdf.m')
   end;
end

k = find((x ~= 0 & p == 0) | (x ~= n & p == 1));
if any(k),
   y(k) = 0;
end

k = find((x == 0 & p == 0) | (x == n & p == 1));
if any(k),
   y(k) = 1;
end

author: Sable

repo: mcbench-benchmarks


text:
function [A,b,P]=gausscppSolve1(A,b)
%function [A,b,P]=gausscppSolve1(A,b)
%forward gaussian elimination with partial pivoting

[m n]=size(A);
if (m~=n), error('matrix must be square!'); end
%initialize pivoting vector
P = (1:m)';
for i=1:(n-1) %sweep through rows
   %maximum of i-th column below diagonal
   [mx imx]= max(abs(A(i:m,i)));
   %row f and column c where |afc| is maximum for the submatrix A(i:m,i)
   f = imx + i - 1;
   c = i;
   %exchange i-th and f-th row for all the system
   dummy = A(i,:);
   A(i,:) = A(f,:);
   A(f,:) = dummy;
   dummy = b(i,:);
   b(i,:) = b(f,:);
   b(f,:) = dummy;
   dummy = P(i);
   P(i) = P(f);
   P(f) = dummy;
   ip1=i+1;
   for j=ip1:m %sweep the remaining rows
     mij=A(j,i)/A(i,i);
     %A(j,:)=A(j,:) - mij* A(i,:);
     A(j,ip1:n)=A(j,ip1:n) - mij* A(i,ip1:n);
     b(j,:)=b(j,:) - mij*b(i,:);
     %keep the pivot for further LU uses
     A(j,i)=mij;
   end
end

author: opencobra

repo: cobratoolbox


text:
function x = buncompress( xR, x, sx )
narginchk(2,3);

if ~isa( xR, 'double' ) && ~isa( xR, 'sparse' ),
    error( 'First argument must be a structure matrix.' );
elseif size( x.basis_, 2 ) ~= size( xR, 1 ),
    error( 'Structure matrix incompatible with vector.' );
elseif nargin < 3 || isempty( sx ),
    sx = size( xR, 2 );
elseif ~cvx_check_dimlist( sx, false ),
    error( 'Third argument must be a size matrix.' );
elseif prod( sx ) ~= isempty( xR ) * prod( x.size_ ) + ~isempty( xR ) * size( xR, 2 ),
    error( 'Incompatible size matrix.' );
end

if isempty( xR ),
    x = cvx( sx, x.basis_ );
else
    x = cvx( sx, x.basis_ * xR );
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%  switching_divided_difference_filter - filter for nonlinear regime-swiching models
%
%  ::
%
%
%    [loglik,Incr,retcode,Filters]=switching_divided_difference_filter(...
%     syst,y,U,z,options)
%
%  Args:
%
%     - **syst** [struct]: structure containing:
%
%           - **PAI00** [vector]: initial probability distributions of regimes
%
%           - **a** [cell]: initial conditions in each regime
%
%           - **Qfunc** [function handle]: transition matrix generator
%
%           - **ff** [function handle]: ft=ff(rt,xt,et), where rt is the
%           regime, xt is the vector of state variables and et the vector of
%           shocks
%
%           - **P** [cell]: initial covariance matrix of the states in each
%           regime
%
%           - **H** [cell]: Measurement error covariance matrices in each regime
%
%           - **SIGeta** [cell]: Covariance matrix of structural shocks.
%
%     - **y** [matrix]: ny x T matrix of data
%
%     - **U** [[]|matrix]: ndx x T matrix of exogenous data
%
%     - **z** [function handle|logical|vector]: linear connection of the
%     observables to the state.
%
%     - **include_in_likelihood** [logical]: selector of increments to include
%     in the likelihood calculation
%
%     - **options** [struct]: structure with various options
%
%  Returns:
%     :
%
%     - **loglik** [scalar]: log likelihood
%
%     - **Incr** [vector]: increments of elements going into the likelihood
%
%     - **retcode** [{0}|integer]: flag for problems.
%
%     - **Filters** [struct]: Filtered, updated and smoothed variables
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function bc = specific_bc(xbd,ybd)
%zero_bc   Reference problems 5.3 and 5.4 boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies streamfunction associated with enclosed flow
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=0*xbd;
return

author: johannesgerer

repo: jburkardt-m


text:
filename='RVE_Square_Triangle';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'chomog_alphabeta','perimeter'};
weights=[1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 0.5 0]';
micro.beta =[1 0.5 0]';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%************************************************************************
%% NTcorr: corrector step for the NT direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%************************************************************************

function [dX,dy,dZ] = NTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
    dX,dZ,coeff,L,X,Z)

global matfct_options solve_ok

printlevel = par.printlevel;
%%
[rhs,EinvRc] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ);
m = length(rp); ncolU = size(coeff.mat12,2);
rhs = [rhs; zeros(m+ncolU-length(rhs),1)];
%%
if strcmp(matfct_options,'chol') || strcmp(matfct_options,'spchol') ...
        || strcmp(matfct_options,'ldl') || strcmp(matfct_options,'spldl')
    [xx,resnrm,solve_ok] = symqmr(coeff,rhs,L,[],[],printlevel);
    if (solve_ok<=0) && (printlevel)
        fprintf('\n  warning: symqmr fails: %3.1f.',solve_ok);
    end
else
    [xx,resnrm,solve_ok] = mybicgstab(coeff,rhs,L,[],[],printlevel);
    if (solve_ok<=0) && (printlevel)
        fprintf('\n  warning: bicgstab fails: %3.1f.',solve_ok);
    end
end
if (printlevel>=3); fprintf('%2.0d ',length(resnrm)-1); end
%%
[dX,dy,dZ] = NTdirfun(blk,At,par,Rd,EinvRc,xx,m);
%%************************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function mask = mri_fillholes(mask, dim)
% MRI_FILLHOLES: Detect background, then fill everything that is not background

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2011

% Initialize two accumulators, for the two directions
acc1 = false(size(mask));
acc2 = false(size(mask));
n = size(mask,dim);
% Process in required direction
switch dim
    case 1
        for i = 2:n
            acc1(i,:,:) = acc1(i-1,:,:) | mask(i,:,:);
        end
        for i = n-1:-1:1
            acc2(i,:,:) = acc2(i+1,:,:) | mask(i,:,:);
        end
    case 2
        for i = 2:n
            acc1(:,i,:) = acc1(:,i-1,:) | mask(:,i,:);
        end
        for i = n-1:-1:1
            acc2(:,i,:) = acc2(:,i+1,:) | mask(:,i,:);
        end
    case 3
        for i = 2:n
            acc1(:,:,i) = acc1(:,:,i-1) | mask(:,:,i);
        end
        for i = n-1:-1:1
            acc2(:,:,i) = acc2(:,:,i+1) | mask(:,:,i);
        end
end
% Combine two accumulators
mask = acc1 & acc2;

author: brainstorm-tools


text:
function [Population,FrontNo] = Select(Population,FrontNo,y)
% And one offspring to the population and delete the worst one

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Identify the solutions in the last front
    Population = [Population,y];
    PopObj     = Population.objs;
    [N,M]      = size(PopObj);
    FrontNo    = UpdateFront(PopObj,FrontNo);

    %% Delete the worst solution
    if max(FrontNo) > 1
        Distance  = pdist2(PopObj,PopObj);
        Distance(logical(eye(N))) = inf;
        LastFront = find(FrontNo==max(FrontNo));
        [~,worst] = min(min(Distance(LastFront,:),[],2));
        worst     = LastFront(worst);
    else
        deltaS = inf(1,N);
        if M == 2
            [~,rank] = sortrows(PopObj);
            for i = 2 : N-1
                deltaS(rank(i)) = (PopObj(rank(i+1),1)-PopObj(rank(i),1)).*(PopObj(rank(i-1),2)-PopObj(rank(i),2));
            end
        elseif N > 1
            deltaS = CalHV(PopObj,max(PopObj,[],1)*1.1,1,10000);
        end
        [~,worst] = min(deltaS);
    end
    FrontNo = UpdateFront(PopObj,FrontNo,worst);
    Population(worst) = [];
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function str = strCell2EnumStr(str_cell, separator)
% SYNTAX:
%   [str] = strcell2enumstr(str_cell, <separator == ' '>);
%
% INPUT:
%   str_cell  = cell array of strings
%   separator = <optional> contains the separator string (white space as default)
%
% OUTPUT:
%   str = single string of the element list separated by the "separator"
%
% DESCRIPTION:
%   Convert to a character array the input cell of strings

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------
    if nargin == 1
        separator = ' ';
    end
    if ischar(str_cell)
        str_cell = {str_cell};
    end
    if ~isempty(str_cell)
        str = sprintf('%d: %s', 0, str_cell{1});
        for i = 2 : numel(str_cell)
            str = sprintf('%s%s%d: %s', str, separator, i-1, str_cell{i});
        end
    else
        str = '';
    end
end

author: goGPS-Project


text:
function z = g2z(g);

% Z = g2z(G)
%
% Hybrid G to Impedance transformation
% only for 2-by-2 matrices
%
% martie 27

if g(1,1) == 0
  disp('correspondent admittance matrix non-existent');
else
z(1,1) = 1/g(1,1);
z(1,2) = -g(1,2)/g(1,1);
z(2,1) = g(2,1)/g(1,1);
z(2,2) = g(2,2) - g(1,2)*g(2,1)/g(1,1);
end;


text:
function F = structBlkDiag(S,field)
% F = structBlkDiag(S,field)
%
% This function creates a block diagional matrix from a field of a struct
% array
%

n = length(S);
cmd = ['F' ' = blkdiag('];
for i=1:n
   cmd = [cmd 'S(' num2str(i) ').' field]; %#ok<AGROW>
   if i < n
      cmd = [cmd ', ']; %#ok<AGROW>
   else
       cmd = [cmd ');'];%#ok<AGROW>
   end
end

eval(cmd);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%++++++++++++++++++++++++++++++++++++++++++++
%PURPOSE: TO DETERMINE THE OPTIMUM VALUE OF
%OVER-RELAXATION (OMEGA) FOR APPLICATION
%IN SOR GS, PICKED OVER MINIMUM VALUE OF
%SOR/GS ITERATION.
%++++++++++++++++++++++++++++++++++++++++++++
function [T]=optimum(om1,om2)
N=200;
eps=1e-6;%error
L=2;
dx=L/N;
n2=25;  %n^2=hP/kA
tb=400+273.15;
ta=34+273.15;
kira=0;

global omega eps N tb ta
warning off

[gs,T]=GS;

omega=om1;
while omega<om2
[sor,T]=SOR;
kira=kira+1;
wye(kira)=sor/gs;
exx(kira)=omega;
omega=omega+0.01;
end

plot(exx,wye);
ylabel('#SOR/#GS RATIO');
xlabel('OMEGA');
grid


text:
% Codes for CVPR-15 work `Face Alignment by Coarse-to-Fine Shape Searching'
% Any question please contact Shizhan Zhu: zhshzhutah2@gmail.com
% Released on July 25, 2015

function reg_model = trainLR_averagedHalfBagging(X,y,iter,regressorInfo)

% train regressor using all the samples only once

assert(size(X,1) == size(y,1));
m = size(X,1);
REG = zeros(size(X,2)+1,size(y,2));

for i = 1:regressorInfo.times
    train_ind = randperm(m,ceil(m/2));
    REG = REG + (1/regressorInfo.times) * ...
        regressorInfo.trainMethod(X(train_ind,:),y(train_ind,:),regressorInfo.lambda(iter));
end;
reg_model = REG;

end

author: zhusz


text:
function [U]=routh(R1,R2)
%these function make next row of routh hurwits
%R2 , R1 must be vector array or cell vector but both format must be same
%length R2 should be same with R1
%these function use match.m for same dimension
%output is a vector
nu=0;
ce=0;
if isnumeric(R1) && isnumeric(R2)
    a=[];
    nu=1;
elseif iscell(R1) && iscell(R2)
    a=sym;
    ce=1;
else
    nu=0;
    ce=0;
end
m=length(R1);
n=length(R2);
if m<n
    R1=match(R1,R2);
end

if n<m
    R2=match(R2,R1);
end
k=2;

if nu==1
    a(1,1)=R1(1,1);
    a(2,1)=R2(1,1);
elseif ce==1
    a(1,1)=R1{1,1};
    a(2,1)=R2{1,1};
end

for f=1:m
    if k~=f
        if nu==1
            a(1,2)=R1(1,k);
            a(2,2)=R2(1,k);
            U(1,f)=(-det(a))/R2(1,1);
        elseif ce==1
            a(1,2)=R1{1,k};
            a(2,2)=R2{1,k};
            U(1,f)=(-det(a))/R2{1,1};
        end
    end
    if k~=m
        k=k+1;
    end
end

if length(U)<m
    U=match(U,R1);
end


text:
% -------------------------------------------------------------------------------------------------
function [net, stats] = experiment_triplet(imdb_video, varargin)
%EXPERIMENT
%   main function - creates a network and trains it on the dataset indexed by imdb_video.
%
%   Luca Bertinetto, Jack Valmadre, Joao Henriques, 2016

% the enhancers are produced by the training net changed with each iteration
% delete pairwise loss.
% -------------------------------------------------------------------------------------------------
    % Default parameters (set the experiment-specific ones in run_experiment)
    opts.net.type = 'alexnet';
    opts.net.conf = struct(); % Options depend on type of net.
    opts.pretrain = false; % Location of model file set in env_paths.
    opts.init.scale = 1;
    opts.init.weightInitMethod = 'xavierimproved';
    opts.init.initBias = 0.1;
    opts.expDir = 'data'; % where to save the trained net
    opts.numFetchThreads = 12; % used by vl_imreadjpg when reading dataset
    opts.validation = 0.1; % fraction of imbd reserved to validation
    opts.exemplarSize = 127; % exemplar (z) in the paper
    opts.instanceSize = 255 - 2*8; % search region (x) in the paper
    opts.loss.type = 'simple';
    opts.loss.rPos = 16; % pixel with distance from center d > rPos are given a negative label
    opts.loss.rNeg = 0; % if rNeg != 0 pixels rPos < d < rNeg are given a neutral label
    opts.loss.labelWeight = 'balanced';
    opts.numPairs =  5.32e4; % Number of example pairs per epoch, if empty, then equal to number of videos.
    opts.randomSeed = 0;
    opts.shuffleDataset = false; % do not shuffle the data to get reproducible experiments
    opts.frameRange = 100; % range from the exemplar in which randomly pick the instance
    opts.gpus = [];
    opts.prefetch = false; % Both get_batch and cnn_train_dag depend on prefetch.
    opts.train.numEpochs = 50;
    opts.train.learningRate = logspace(-2, -5, opts.train.numEpochs);
    opts.train.weightDecay = 5e-4;
    opts.train.batchSize = 8; % we empirically observed that small batches work better
    opts.train.profile = false;
    % add parameter
%     opts.train.exemplarSize = opts.exemplarSize; % exemplar (z) in the paper
%     opts.train.instanceSize = opts.instanceSize; % search region (x) in the paper
%     opts.loss.weight  = 0.01;
    % Data augmentation settings
    opts.subMean = false;
    opts.colorRange = 255;
    opts.augment.translate = true;
    opts.augment.maxTranslate = 4;
    opts.augment.stretch = true;
    opts.augment.maxStretch = 0.05;
    opts.augment.color = true;
    opts.augment.grayscale = 0; % likelihood of using grayscale pair
    % Override default parameters if specified in run_experiment

    % Get environment-specific default paths.
    opts = env_paths_training(opts);
    opts.train.gpus = opts.gpus;
    opts.train.prefetch = opts.prefetch;

    opts = vl_argparse(opts, varargin);
% -------------------------------------------------------------------------------------------------
    % Get ImageNet Video metadata
    if isempty(imdb_video)
        fprintf('loading imdb video...\n');
        imdb_video = load(opts.imdbVideoPath);
        imdb_video = imdb_video.imdb_video;
    end

    % Load dataset statistics
    [rgbMean_z, rgbVariance_z, rgbMean_x, rgbVariance_x] = load_stats(opts);
    if opts.shuffleDataset
        s = RandStream.create('mt19937ar', 'Seed', 'shuffle');
        opts.randomSeed = s.Seed;
    end

    opts.train.expDir = opts.expDir;

    rng(opts.randomSeed); % Re-seed before calling make_net.

    % -------------------------------------------------------------------------------------------------
    net = make_net(opts);
    % -------------------------------------------------------------------------------------------------

    [imdb_video, imdb] = choose_val_set(imdb_video, opts);

    [resp_sz, resp_stride] = get_response_size(net, opts);
    % We want an odd number so that we can center the target in the middle
    assert(all(mod(resp_sz, 2) == 1), 'resp. size is not odd');

    [net, derOutputs, label_inputs_fn] = setup_loss(net, resp_sz, resp_stride, opts.loss);

%     enhancer = zeros(a_feat_sz(1),a_feat_sz(2),256,opts.train.batchSize,'single');
% %     imout_t(end)=1;
% %     if numel(opts.train.gpus) >= 1
% %         imout_t = gpuArray(imout_t);
% %     end
%     t_labels = ones(opts.train.batchSize,1);
%     neg_labels = -ones(opts.train.batchSize,1);
    batch_fn = @(db, batch) get_batch(db, batch, ...
                                        imdb_video, ...
                                        opts.rootDataDir, ...
                                        numel(opts.train.gpus) >= 1, ...
                                        struct('exemplarSize', opts.exemplarSize, ...
                                               'instanceSize', opts.instanceSize, ...
                                               'frameRange', opts.frameRange, ...
                                               'subMean', opts.subMean, ...
                                               'colorRange', opts.colorRange, ...
                                               'stats', struct('rgbMean_z', rgbMean_z, ...
                                                               'rgbVariance_z', rgbVariance_z, ...
                                                               'rgbMean_x', rgbMean_x, ...
                                                               'rgbVariance_x', rgbVariance_x), ...
                                               'augment', opts.augment, ...
                                               'prefetch', opts.train.prefetch, ...
                                               'numThreads', opts.numFetchThreads), ...
                                        label_inputs_fn);

    opts.train.derOutputs = derOutputs;
    opts.train.randomSeed = opts.randomSeed;
    % -------------------------------------------------------------------------------------------------
    % Start training
    [net, stats] = cnn_train_dag(net, imdb, batch_fn, opts.train);
    % -------------------------------------------------------------------------------------------------
end

% -----------------------------------------------------------------------------------------------------
function [rgbMean_z, rgbVariance_z, rgbMean_x, rgbVariance_x] = load_stats(opts)
% Dataset image statistics for data augmentation
% -----------------------------------------------------------------------------------------------------
    stats = load(opts.imageStatsPath);
    % Subtracted if opts.subMean is true
    if ~isfield(stats, 'z')
        rgbMean = reshape(stats.rgbMean, [1 1 3]);
        rgbMean_z = rgbMean;
        rgbMean_x = rgbMean;
        [v,d] = eig(stats.rgbCovariance);
        rgbVariance_z = 0.1*sqrt(d)*v';
        rgbVariance_x = 0.1*sqrt(d)*v';
    else
        rgbMean_z = reshape(stats.z.rgbMean, [1 1 3]);
        rgbMean_x = reshape(stats.x.rgbMean, [1 1 3]);
        % Set data augmentation statistics, used if opts.augment.color is true
        [v,d] = eig(stats.z.rgbCovariance);
        rgbVariance_z = 0.1*sqrt(d)*v';
        [v,d] = eig(stats.x.rgbCovariance);
        rgbVariance_x = 0.1*sqrt(d)*v';
    end
end

% -------------------------------------------------------------------------------------------------
function net = make_net(opts)
% -------------------------------------------------------------------------------------------------

    net = make_siam_adapt2(opts);
%     net = make_siameseFC(opts);

    % Save the net graph to disk.
    inputs = {'exemplar', [opts.exemplarSize*[1 1] 3 opts.train.batchSize], ...
              'instance', [opts.instanceSize*[1 1] 3 opts.train.batchSize]};
    net_dot = net.print(inputs, 'Format', 'dot');
    if ~exist(opts.expDir)
        mkdir(opts.expDir);
    end
    f = fopen(fullfile(opts.expDir, 'arch.dot'), 'w');
    fprintf(f, net_dot);
    fclose(f);
end

% -------------------------------------------------------------------------------------------------
function [resp_sz, resp_stride] = get_response_size(net, opts)
% -------------------------------------------------------------------------------------------------

    sizes = net.getVarSizes({'exemplar', [opts.exemplarSize*[1 1] 3 256], ...
                             'instance', [opts.instanceSize*[1 1] 3 256]});
    resp_sz = sizes{net.getVarIndex('score')}(1:2);
    rfs = net.getVarReceptiveFields('exemplar');
    resp_stride = rfs(net.getVarIndex('score')).stride(1);
    assert(all(rfs(net.getVarIndex('score')).stride == resp_stride));
end

% -------------------------------------------------------------------------------------------------
function [net, derOutputs, inputs_fn] = setup_loss(net, resp_sz, resp_stride, loss_opts)
% Add layers to the network, specifies the losses to minimise, and
% constructs a function that returns the inputs required by the loss layers.
% -------------------------------------------------------------------------------------------------
net.addLayer('tri_select', select_pairs(), ...
                 {'score', 'eltwise_label'}, ...
                 {'tri_pairs'}, ...
                 {});
    net.addLayer('objective_softmaxlog', ...
                 dagnn.Loss('loss', 'softmaxlog'), ...
                 {'tri_pairs','tri_pairs_label'}, 'objective_softmaxlog');
   net.vars(end).precious = 1;
%     %xcorr ac  correlation between exemplar and positive instance
%     net.addLayer('xcorr_ac', XCorr(), ...
%                  {'a_feat', 'c_feat'}, ...
%                  {'xcorr_ac_out'}, ...
%                  {});
%     add_adjust_layer(net, 'adjust_ac', 'xcorr_ac_out', 'score_ac', ...
%                  {'adjust_f', 'adjust_b'}, 1e-3, 0, 0, 1);
%     %xcorr ad   correlation between exemplar and negative instance
%     net.addLayer('xcorr_ad', XCorr(), ...
%                  {'a_feat', 'd_feat'}, ...
%                  {'xcorr_ad_out'}, ...
%                  {});
%     add_adjust_layer(net, 'adjust_ad', 'xcorr_ad_out', 'score_ad', ...
%                  {'adjust_f', 'adjust_b'}, 1e-3, 0, 0, 1);
%
% %     net.addLayer('concat', dagnn.Concat(), ...
% %                  {'score_ac', 'score_ad'}, ...
% %                  {'concat_out'}, ...
% %                  {});
%     % softmaxloss between score_ac and score_ad
%     net.addLayer('objective_softmax', ...
%                  TripletLoss2(), ...
%                  {'score_ac', 'score_ad'}, 'objective_softmax');
% %                  net.layers(end).block.opts = [...
% %         net.layers(end).block.opts, ...
% %         {'instanceWeights', loss_opts.weight}];
% %     add_adjust_layer(net, 'adjust_loss', 'objective_softmax', 'loss_softmax',  ...
% %                  {'adjust_f_sm', 'adjust_b_sm'}, 1e-3, 0, 1, 0);
%
% %     % create label and weights for logistic loss
% %     net.addLayer('objective', ...
% %                  dagnn.Loss('loss', 'logistic'), ...
% %                  {'score', 'eltwise_label'}, 'objective');
% %     % adding weights to loss layer

    [pos_eltwise, instanceWeight] = create_labels(...
        resp_sz, loss_opts.labelWeight, ...
        loss_opts.rPos/resp_stride, loss_opts.rNeg/resp_stride);
    neg_eltwise = [];   % no negative pairs at the moment

    ind_pos = pos_eltwise>0;
        ind_neg = pos_eltwise<0;
        ind_pos = ind_pos(:);
        ind_neg = ind_neg(:);
        n_pos = sum(ind_pos);
        n_neg = sum(ind_neg);
        ind_select = ones(n_pos,n_neg)>0;
        n_pairs = sum(ind_select(:));
        instanceWeight = ones(1,n_pairs,'single')/n_pairs;
    net.layers(end).block.opts = [...
        net.layers(end).block.opts, ...
        {'instanceWeights', instanceWeight}];

%         net.addLayer('sum_loss', ...
%                  dagnn.Sum('numInputs',2), ...
%                  {'objective', 'objective_softmax'}, 'sum_loss');

%         net.addLayer('sum_loss', ...
%                  weightSum('numInputs',2,'ws',[1,loss_opts.weight]), ...
%                  {'objective', 'objective_softmax'}, 'sum_loss');

%      add_norm_weight_sum_layer(net, 'sum_loss', {'objective', 'objective_softmax'}, {'sum_loss'}, ...
%                  {'ws'}, [0.9,0.1], 1);

%     % loss between a_feat and c_feat
%     net.addLayer('objective_ac', ...
%                  dagnn.Loss('loss', 'logistic'), ...
%                  {'score_ac', 'pos_label'}, 'objective_ac');
%     add_adjust_layer(net, 'adjust', 'score_ac', 'adjust_score_ac',  ...
%                  {'adjust_f_ac', 'adjust_b_ac'}, 1e-3, 0, 1, 0);
%     % loss between a_feat and d_feat
%     net.addLayer('objective_ad', ...
%                  dagnn.Loss('loss', 'logistic'), ...
%                  {'score_ad', 'neg_label'}, 'objective_ad');
%     add_adjust_layer(net, 'adjust', 'score_ad', 'adjust_score_ad',  ...
%                  {'adjust_f_ad', 'adjust_b_ad'}, 1e-3, 0, 1, 0);
%     net.layers(end).block.opts = [...
%         net.layers(end).block.opts, ...
%         {'instanceWeights', loss_opts.weight}];

%     derOutputs = {'sum_loss', 1};
%     derOutputs = {'objective', 1, 'objective_softmax',1};
    derOutputs = { 'objective_softmaxlog',1};

    inputs_fn = @(labels, obj_sz_z, obj_sz_x) get_label_inputs_simple(...
        labels, obj_sz_z, obj_sz_x, pos_eltwise, neg_eltwise, n_pairs);

    net.addLayer('errdisp', centerThrErr(), {'score', 'label'}, 'errdisp');
    net.addLayer('errmax', MaxScoreErr(), {'score', 'label'}, 'errmax');
end

% -------------------------------------------------------------------------------------------------
function inputs = get_label_inputs_simple(labels, obj_sz_z, obj_sz_x, pos_eltwise, neg_eltwise, n_pairs)
% GET_LABEL_INPUTS_SIMPME returns the network inputs that specify the labels.
%
% labels -- Label of +1 or -1 per image pair, size [1, n].
% obj_sz_z -- Size of exemplar box, dims [2, n].
% obj_sz_x -- Size of instance box, dims [2, n].
% -------------------------------------------------------------------------------------------------

%     pos = (labels > 0);
%     neg = (labels < 0);
%
%     resp_sz = size(pos_eltwise);
%     eltwise_labels = zeros([resp_sz, 1, numel(labels)], 'single');
%     eltwise_labels(:,:,:,pos) = repmat(pos_eltwise, [1 1 1 sum(pos)]);
%     eltwise_labels(:,:,:,neg) = repmat(neg_eltwise, [1 1 1 sum(neg)]);
        label0 = ones(1,n_pairs,1,numel(labels),'single');
    inputs = {'label', labels, ...
              'eltwise_label', pos_eltwise,...
              'tri_pairs_label',label0};
%           inputs = {'label', labels};
end

% -------------------------------------------------------------------------------------------------
function [imdb_video, imdb] = choose_val_set(imdb_video, opts)
% Designates some examples for validation.
% It modifies imdb_video and constructs a dummy imdb.
% -------------------------------------------------------------------------------------------------
    TRAIN_SET = 1;
    VAL_SET = 2;

    % set opts.validation to validation and the rest to training.
    size_dataset = numel(imdb_video.id);
    size_validation = round(opts.validation * size_dataset);
    size_training = size_dataset - size_validation;
    imdb_video.set = uint8(zeros(1, size_dataset));
    imdb_video.set(1:size_training) = TRAIN_SET;
    imdb_video.set(size_training+1:end) = VAL_SET;

    %% create imdb of indexes to imdb_video
    % train and val from disjoint video sets
    imdb = struct();
    imdb.images = struct(); % we keep the images struct for consistency with cnn_train_dag (MatConvNet)
    imdb.id = 1:opts.numPairs;
    n_pairs_train = round(opts.numPairs * (1-opts.validation));
    imdb.images.set = uint8(zeros(1, opts.numPairs)); % 1 -> train
    imdb.images.set(1:n_pairs_train) = TRAIN_SET;
    imdb.images.set(n_pairs_train+1:end) = VAL_SET;
end

% -------------------------------------------------------------------------------------------------
function inputs = get_batch(db, batch, imdb_video, data_dir, use_gpu, sample_opts, label_inputs_fn)
% Returns the inputs to the network.
% -------------------------------------------------------------------------------------------------

    [imout_z, imout_x, labels, sizes_z, sizes_x] = vid_get_random_batch(...
        db, imdb_video, batch, data_dir, sample_opts);
    if use_gpu
        imout_z = gpuArray(imout_z);
        imout_x = gpuArray(imout_x);
    end
    % Constructs full label inputs from output of vid_get_random_batch.
    label_inputs = label_inputs_fn(labels, sizes_z, sizes_x);
    inputs = [{'exemplar', imout_z, 'instance', imout_x}, label_inputs];
end

author: shenjianbing

repo: TripletTracking


text:
function [Nv, VX, VY, VZ, K, EToV] = MeshReaderGambit3D(FileName)

% function [Nv, VX, VY, VZ, K, EToV] = MeshReaderGambit3D(FileName)
% Purpose  : Read in basic grid information to build grid
% NOTE     : gambit *.neu format is assumed

Fid = fopen(FileName, 'rt');

% read intro
for i=1:6
  line = fgetl(Fid);
end

% fine number of nodes and number of elements
dims = fscanf(Fid, '%d');

Nv = dims(1); K = dims(2);

for i=1:2
  line = fgetl(Fid);
end

% read node coordinates
xyz = fscanf(Fid, '%lf', [4, Nv]);
xyz = xyz(2:4, :);
VX = xyz(1,:); VY = xyz(2,:); VZ = xyz(3,:);

for i=1:3
  line = fgetl(Fid);
end

% read element to node connectivity
EToV = zeros(K, 4);
for k = 1:K
  line   = fgetl(Fid);
  tmpcon = sscanf(line, '%lf');
  EToV(k,1:4) = tmpcon(4:7);
end
return

repo: Aero-matlab


text:
function [Lc,Mc,Rc] = compressfactors(L,M,R)

Lc = {};
Rc = {};
Mc = {};
taken = zeros(1,length(M));
for i = 1:length(M)
    Lsameas = [];
    Rsameas = [];
    if ~taken(i)
        for j = i+1:length(M)
            if isequal(M{i},M{j})
                if isequal(L{i},L{j})
                    Lsameas = [Lsameas j];
                    taken(j) = 1;
                elseif isequal(R{i},R{j})
                    Rsameas = [Rsameas j];
                    taken(j) = 1;
                end
            end
        end
        if isempty(Rsameas) & ~isempty(Lsameas)
            Lc{end+1} = L{i};
            Mc{end+1} = M{i};
            Rc{end+1} = R{i};
            for j = 1:length(Lsameas)
                Rc{end} = Rc{end} + R{Lsameas(j)};
            end
        elseif isempty(Lsameas) & ~isempty(Rsameas)
            Lc{end+1} = L{i};
            Mc{end+1} = M{i};
            Rc{end+1} = R{i};
            for j = 1:length(Rsameas)
                Lc{end} = Lc{end} + L{Rsameas(j)};
            end
        elseif ~taken(i)
            Lc{end+1} = L{i};
            Mc{end+1} = M{i};
            Rc{end+1} = R{i};
        end
    end
end


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcFlipZ(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

if ~isempty(handles.V.Segs)
    choice = questdlg('Segmentation mask is detected, transform operation will reset them, preceed?','Mask Reset', ...
                      'No, go save mask','Yes','No, go save mask');
    if isempty(choice)
        warndlg('Transform is cancelled.');
        return;
    end
    % Handle response
    switch choice
        case 'No, go save mask'
            warndlg('Save your mask before transform.');
            return;
    end

    handles.Mask=handles.Mask*0;
    handles.V.Segs=[];
end

% close 3D slicer
global Figure_handles
if isfield(Figure_handles,'MU_display2')
    slicer_display_handles=guidata(Figure_handles.MU_display2);
    if Figure_handles.MU_display == slicer_display_handles.Parent
        close(Figure_handles.MU_display2);
    end
end

handles.TMatrix=flipdim(handles.TMatrix,3);
handles.Mask=flipdim(handles.Mask,3);

MergeM=get(handles.Matrix_name_edit,'String');
set(handles.Matrix_name_edit,'String',[MergeM '_fpz']);

% update current display matrix
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function params = mrParamsCorAnal(mr);
% Put up a dialog to get params for a coherence analysis.
%
% params = mrParamsCorAnal(mr);
%
% This is an attempt to emulate the style of mrLoadRet end+1.0,
% in which analyses have related GUIs to get params.
% So, corAnal has a corresponding corAnalGUI to get cor anal
% params. This is the equivalent mrVista 2.0 implementation.
%
% Thought on naming conventions: for other analyses, it might
% be good to also start with mrParams -- mr since the params are
% operating on mr objects, and params, because that is the output
% of this dialog.
%
% Note: an mr struct must be entered, since it's needed for
% the default set of frames (all).
%
% ras, 07/05.

%%%%% default params (if omitted from params struct):
nCycles = 6;   % just a guess that many scans use
noiseBand = 0; % use all available spectral bands for analysis
saveDir = fileparts(mr.path);
detrend = 1;

% put up dialog to get params
dlg(1).fieldName = 'nCycles';
dlg(1).style = 'edit';
dlg(1).string = 'Number of Cycles in Time Series';
dlg(1).value = num2str(nCycles);

dlg(end+1).fieldName = 'detrend';
dlg(end).style = 'checkbox';
dlg(end).string = 'Detrend Time Series?';
dlg(end).value = detrend;

dlg(end+1).fieldName = 'noiseBand';
dlg(end).style = 'edit';
dlg(end).string = 'Noise Band Option';
dlg(end).value = num2str(noiseBand);

dlg(end+1).fieldName = 'frames';
dlg(end).style = 'edit';
dlg(end).string = 'Frames to Analyze';
dlg(end).value = sprintf('1:%i',size(mr.data,4));

dlg(end+1).fieldName = 'saveDir';
dlg(end).style = 'edit';
dlg(end).string = 'Save directory';
dlg(end).value = saveDir;

params = generalDialog(dlg,sprintf('Coherence Analysis: %s',mr.name));

% if canceled, exit quietly
if isempty(params), disp('User Aborted corAnal.'); return;   end

% parse params
params.nCycles = str2num(params.nCycles);
params.noiseBand = str2num(params.noiseBand);
params.frames = str2num(params.frames);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef TP3 < PROBLEM
% <multi> <real> <large/none> <robust>
% Test problem for robust multi-objective optimization
% delta --- 0.05 --- Maximum disturbance degree
% H     ---   50 --- Number of disturbances

%------------------------------- Reference --------------------------------
% A. Gaspar-Cunha, J. Ferreira, and G. Recio, Evolutionary robustness
% analysis for multi-objective optimization: benchmark problems, Structural
% and Multidisciplinary Optimization, 2014, 49: 771-793.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        delta;      % Maximum disturbance degree
        H;          % Number of disturbances
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.delta,obj.H] = obj.ParameterSet(0.05,50);
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj(:,1) = 1 - PopDec(:,1).^2;
            g = 1 + 10*mean(PopDec(:,2:end),2);
            PopObj(:,2) = sin(pi/2*PopDec(:,1)).*g;
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(~,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = sin(pi/2*sqrt(1-R(:,1)));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case {'Mean_IGD','Mean_HV','Worst_IGD','Worst_HV'}
                    score = feval(metName,Population,obj);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Perturb solutions multiple times
        function PopX = Perturb(obj,PopDec,N)
            if nargin < 3; N = obj.H; end
            Delta = repmat(obj.delta.*(obj.upper-obj.lower),N*size(PopDec,1),1);
            w     = UniformPoint(N,obj.D,'Latin');
            Dec   = 2*Delta.*w(reshape(repmat(1:end,size(PopDec,1),1),1,[]),:) + repmat(PopDec,N,1) - Delta;
            Dec   = obj.CalDec(Dec);
            PopX  = SOLUTION(Dec,obj.CalObj(Dec),obj.CalCon(Dec));
        end
    end
end


text:
function demo3 (C, sym, name)
%DEMO3: Cholesky update/downdate
%
% Example:
%   demo3 (C, 1, 'name of system')
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

clf
subplot (2,2,1) ; cspy (C) ;
title (name, 'FontSize', 16, 'Interpreter', 'none') ;
drawnow

[m n] = size (C) ;
if (m ~= n | ~sym)                                                          %#ok
    return ;
end

b = rhs (n) ;
fprintf ('chol then update/downdate ') ;
print_order (0) ;

tic ;
[L,p] = cs_chol (C) ;
t = toc ;
fprintf ('\nchol  time: %8.2f\n', t) ;

subplot (2,2,2) ; cspy (L) ; title ('L') ;
drawnow

tic ;
x = b (p) ;
x = cs_lsolve (L,x) ;
x = cs_ltsolve (L,x) ;
x (p) = x ;
t = toc ;
fprintf ('solve time: %8.2f\n', t) ;

fprintf ('original: ') ;
print_resid (C, x, b) ;

k = fix (n/2) ;
w = L(k,k) * sprand (L (:,k)) ;

parent = cs_etree (C (p,p)) ;

tic ;
L2 = cs_updown (L, w, parent, '+') ;
t1 = toc ;
fprintf ('update:   time: %8.2f\n', t1) ;

subplot (2,2,3) ; cspy (L2) ; title ('updated L') ;
subplot (2,2,4) ; cspy (L-L2) ; title ('L - updated L') ;
drawnow

tic ;
x = b (p) ;
x = cs_lsolve (L2,x) ;
x = cs_ltsolve (L2,x) ;
x (p) = x ;
t = toc ;

w2 = sparse (n,1) ;
w2 (p) = w ;                    % w2 = P'*w
wt = cs_transpose (w2) ;
ww = cs_multiply (w2,wt) ;
E = cs_add (C, ww, 1, 1) ;      % E = C + w2*w2' ;

fprintf ('update:   time: %8.2f (incl solve) ', t1+t) ;
print_resid (E, x, b) ;

tic
[L,p2] = cs_chol (E) ;
x = b (p2) ;
x = cs_lsolve (L,x) ;
x = cs_ltsolve (L,x) ;
x (p2) = x ;
t = toc ;
fprintf ('rechol:   time: %8.2f (incl solve) ', t) ;
print_resid (E, x, b) ;

tic ;
L3 = cs_updown (L2, w, parent, '-') ;
t1 = toc ;
fprintf ('downdate: time: %8.2f\n', t1) ;

tic ;
x = b (p) ;
x = cs_lsolve (L3,x) ;
x = cs_ltsolve (L3,x) ;
x (p) = x ;
t = toc ;
fprintf ('downdate: time: %8.2f (incl solve) ', t1+t) ;
print_resid (C, x, b) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ EEGbase ] = eeg_baseline(EEGdata,baseline_samples)

% eeg_baseline - remove mean of baseline period from all EEG epochs
%
% [ EEGbase ] = eeg_baseline(EEGdata,baseline_samples)
%
% It is assumed that EEGdata is an MxNxP matrix, with M data samples (ie,
% time), N channels and P epochs (if P=1, EEGdata is MxN).
%
% baseline_samples is an array of data sample indices (of M) that define
% the baseline period.  If baseline_samples is a scalar, the function
% assumes the baseline index array is 1:baseline_samples.
%
% The function calculates the baseline mean for each epoch as follows:
%
% baseline_mean = mean(EEGdata(baseline_samples,:,epoch));
%

% Copyright (C) 2004  Darren L. Weber
%
% 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:50 $
% Created:  10/2004, copyright 2004 Darren.Weber_at_radiology.ucsf.edu
% Modified: 05/2005, Darren.Weber_at_radiology.ucsf.edu
%                    baseline_samples can now be an array or a scalar input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

version = '[$Revision: 1.1 $]';
fprintf('\nEEG_BASELINE [v%s]\n',version(12:16)); tic

if ~exist('baseline_samples','var'),
    error('no baseline_samples specified');
end
if isempty(baseline_samples),
    error('no baseline_samples specified');
end
if length(baseline_samples) == 1,
    baseline_samples = 1:baseline_samples;
end

number_samples = size(EEGdata,1);

for epoch = 1:size(EEGdata,3),

    baseline_mean = mean(EEGdata(baseline_samples,:,epoch));

    baseline_mean = repmat(baseline_mean,number_samples,1);

    EEGbase(:,:,epoch) = EEGdata(:,:,epoch) - baseline_mean;

end

t = toc; fprintf('...done (%6.2f sec).\n\n',t);
return

author: PatternRecognition

repo: OpenBMI


text:
function xdot=dpendulum(t,x)
%% Double Pendulum ODE
%
% ivp=[gamma0; dtgamma0; alpha0; dtalpha0; beta0; lambda; omega; psi; eta];

beta0=x(5); lambda=x(6); omega=x(7); psi=x(8); eta=x(9);
xdot=zeros(9,1); % a column vector

C = cos(beta0)*cos(x(3)-x(1))-sin(beta0)*sin(x(3)-x(1));
S = sin(beta0)*cos(x(3)-x(1))+cos(beta0)*sin(x(3)-x(1));

xdot(1) = x(2);
xdot(2) = ((omega^2)*sin(x(1))-psi*(S*(x(4)^2+C*x(2)^2*eta)+...
    C*(lambda^2)*sin(x(3))))/(-1+(C^2)*eta*psi);
xdot(3) = x(4);
xdot(4) = (S*eta*(x(2)^2+C*x(4)^2*psi)-C*eta*(omega^2)*sin(x(1))+...
    (lambda^2)*sin(x(3)))/(-1+(C^2)*eta*psi);


text:
classdef TableExporter < handle

    properties (Access = public)
        writeVariableNames = true;
        writeRowNames = false;
        delimiter = ',';
    end

    methods (Access = public)

        function testExportTableForFileNames(obj,fileNames)

            labelingStrategy = DefaultLabelingStrategy();
            tableCreator = TableCreator();
            table = tableCreator.loadFeaturesTableFiltered(fileNames, labelingStrategy);
            featureSelector = FeatureSelector();
            bestFeatureIdxs = featureSelector.getBestNFeatures(9);
            table = table(:,[bestFeatureIdxs end]);
            obj.exportTable(table);
        end

        function exportTable(obj,table,fileName)
            writetable(table,fileName,'Delimiter',obj.delimiter,'WriteVariableNames',obj.writeVariableNames,'WriteRowNames',obj.writeRowNames);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
% 梯形公式的加速--外推技巧
clear; clc;
f = @(x) sqrt(x.^3);
a = 0;
b = 1;
tol = 1e-7;
maxit = 100; % 最大迭代步数

T(1,1) = (b-a)/2 * (f(a) + f(b));
fprintf(' T(1,1)=%.8f \n',T(1,1));
k = 1;  err = 1.0;
for k = 1 : maxit
    h = (b-a)/(2^(k-1));
    T(k+1,1) = T(k,1)/2 + (h/2) * sum(f(a+0.5*h+[0:2^(k-1)-1]*h));
    fprintf(' T(%d,1)=%.8f ',k+1,T(k+1,1));
    for j = 1 : k
        T(k+1,j+1) = (4^j * T(k+1,j) - T(k,j) )/(4^j-1);
        fprintf(' T(%d,%d)=%.8f ',k+1,j+1,T(k+1,j+1));
        err = abs(T(k+1,j+1)-T(k+1,j));
        if (err<tol), break, end % check convergence
    end
    fprintf('\n')
    if (err<tol), break, end % check convergence
end

repo: NumericalAnalysis


text:
function varargout = clickpoint(varargin)
% VL_CLICKPOINT  Select a point by clicking
%   SEL = VL_CLICKPOINT(V) let the user vl_click a point in the current
%   figure and returns the index v of the closest point (in Euclidean
%   norm) in the collection V. The 2xK matrix V has a a column for
%   each point.
%
%   The user can abort the operation by pressing any key. In this case
%   the function returns the empty matrix.
%
%   VL_CLICKPOINT(V,N) selects N points in a row. The user can stop the
%   selection at any time by pressing any key. In this case the
%   partial selection is returned. This can be used in combination
%   with N=inf to get an arbitrary number of points.
%
%   VL_CLICKPOINT() accepts the following options:
%
%   PlotMarker:: [0]
%     Put a marker as points are selected. The markers are deleted on
%     exiting the function.
%
%   See also: VL_CLICK(), VL_HELP().
[varargout{1:nargout}] = vl_clickpoint(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama