text:
function handle = draw_waypoints(fig_handle, wpt_list, radius, handle, mode)
% DRAW_WAYPOINTS - plot waypoints

    if wpt_list(1,4)==-9999 % check to see if Dubins paths
        XX = [wpt_list(:,1)];
        YY = [wpt_list(:,2)];
        ZZ = [wpt_list(:,3)];
    else
        XX = [];
        YY = [];
        ZZ = [];
        for i=2:size(wpt_list,1)
            dubinspath = compute_dubins_param(wpt_list(i-1,:),wpt_list(i,:),radius);
            [tmpX,tmpY,tmpZ] = points_along_dubins_path(dubinspath,0.1);
            XX = [XX; tmpX];
            YY = [YY; tmpY];
            ZZ = [ZZ; tmpZ];
        end
%         ZZ = wpt_list(i,3)*ones(size(XX));
    end

    if isempty(handle)
        handle = plot3(fig_handle.Children, YY, XX, -ZZ, 'b');
    else
        set(handle,'XData', YY, 'YData', XX, 'ZData', -ZZ);
        drawnow
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X,Y,Z] = points_along_dubins_path(dubinspath,Del)
% Points Along Dubins Path -
%   Find points along Dubin's path separted by Del (to be used in
%   collision detection)

  % points along start circle
  th1 = mod(atan2(dubinspath.ps(2)-dubinspath.cs(2),dubinspath.ps(1)-dubinspath.cs(1)),2*pi);
  th2 = mod(atan2(dubinspath.w1(2)-dubinspath.cs(2),dubinspath.w1(1)-dubinspath.cs(1)),2*pi);
  if dubinspath.lams>0
      if th1>=th2
        th = [th1:Del:2*pi,0:Del:th2];
      else
        th = [th1:Del:th2];
      end
  else
      if th1<=th2
        th = [th1:-Del:0,2*pi:-Del:th2];
      else
        th = [th1:-Del:th2];
      end
  end
  X = [];
  Y = [];
  Z = [];
  for i=1:length(th)
    X = [X; dubinspath.cs(1)+dubinspath.R*cos(th(i))];
    Y = [Y; dubinspath.cs(2)+dubinspath.R*sin(th(i))];
    Z = [Z; dubinspath.cs(3)];
  end

  % points along straight line
  sig = 0;
  while sig<=1
      X = [X; (1-sig)*dubinspath.w1(1) + sig*dubinspath.w2(1)];
      Y = [Y; (1-sig)*dubinspath.w1(2) + sig*dubinspath.w2(2)];
      Z = [Z; (1-sig)*dubinspath.w1(3) + sig*dubinspath.w2(3)];
      sig = sig + Del;
  end

  % points along end circle
  th2 = mod(atan2(dubinspath.pe(2)-dubinspath.ce(2),dubinspath.pe(1)-dubinspath.ce(1)),2*pi);
  th1 = mod(atan2(dubinspath.w2(2)-dubinspath.ce(2),dubinspath.w2(1)-dubinspath.ce(1)),2*pi);
  if dubinspath.lame>0
      if th1>=th2
        th = [th1:Del:2*pi,0:Del:th2];
      else
        th = [th1:Del:th2];
      end
  else
      if th1<=th2
        th = [th1:-Del:0,2*pi:-Del:th2];
      else
        th = [th1:-Del:th2];
      end
  end
  for i=1:length(th)
    X = [X; dubinspath.ce(1)+dubinspath.R*cos(th(i))];
    Y = [Y; dubinspath.ce(2)+dubinspath.R*sin(th(i))];
    Z = [Z; dubinspath.ce(3)];
  end
end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function Psi0 = linard2VardistPsi0Compute(linard2kern, vardist)

% LINARD2VARDISTPSI0COMPUTE Description

% VARGPLVM

A = linard2kern.inputScales;
Psi0 = sum(A.*sum((vardist.means.*vardist.means) + vardist.covars,1));

%Psi0covs = sum(A.*sum(vardist.covars,1));

%psi00 = kernCompute(linard2kern,vardist.means);
%ok = 0;
%for n=1:size(vardist.means,1)
%    ok = ok + trace(diag(A)*(vardist.means(n,:)'*vardist.means(n,:))) + trace(diag(A)*diag(vardist.covars(n,:)));
%end
%ok
%pause
%Psi0 = Psi0means + Psi0covs;

author: SheffieldML

repo: vargplvm


text:
% Demo: learn a sinc. Run and compare all algorithms using their default
% parameters.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

close all
clear
rng('default'); rng(1)

%% PARAMETERS

N = 1000; % number of training data
N_test = 500; % number of test data
SNR = 20; % SNR in dB

%% PROGRAM

% generate data
x = randn(N,1);
x_test = linspace(min(x),max(x),N_test)';
y_ref = sinc([x;x_test]);
y = y_ref + sqrt(10^(-SNR/10)*var(y_ref))*randn(N+N_test,1);
y_test = y_ref(N+1:N+N_test);

% get list of kernel adaptive filters in 'lib' folder
fdir = fileparts(which('kafbox_template.m'));
files = dir(fullfile(fdir,'*.m'));
[~,algorithms] = cellfun(@fileparts, {files.name}, 'UniformOutput',0);
for i=length(algorithms):-1:1
    if ~exist(algorithms{i},'class')
        algorithms(i) = []; % remove files that do not represent classes
    end
end

% perform online learning for each algorithm
fprintf('\n')
num_alg = length(algorithms);
titles = cell(num_alg,1);
MSE = zeros(num_alg,1);
Y_est = zeros(N_test,num_alg);
for algo_ind=1:num_alg
    t1 = tic;
    algorithm = algorithms{algo_ind};
    fprintf('%2d. %9s: ',algo_ind,upper(algorithm));
    titles{algo_ind} = strrep(upper(algorithm),'_','\_');

    kaf = feval(algorithm);
    for i=1:N
        if ~mod(i,floor(N/10)), fprintf('.'); end
        kaf.train(x(i),y(i));
    end
    y_est = kaf.evaluate(x_test);
    Y_est(:,algo_ind) = y_est;
    MSE(algo_ind) = mean((y_test-y_est).^2);

    fprintf(' %.2fs. MSE=%3.2fdB\n',toc(t1),10*log10(MSE(algo_ind)))
end

%% OUTPUT

% plot results in different "leagues"
[MSE_sorted,ind] = sort(MSE,'descend');
num_fig = ceil(num_alg/5);

remaining = num_alg;
titles{num_alg+1} = 'data';
for fig_ind=num_fig:-1:1
    figure; hold all
    plot(x,y(1:N),'.')

    rm = rem(remaining,5);
    num_in_league = (rm==0)*5 + rm;
    % plot the results for the num_in_league worst results
    league_inds = num_alg-remaining+num_in_league:-1:num_alg-remaining+1;
    for i=league_inds
        plot(x_test,Y_est(:,ind(i)),'LineWidth',2)
    end
    title(sprintf('League %d',fig_ind))
    legend(titles([num_alg+1; ind(league_inds)]))

    axis([min(x)-0.5 max(x)+0.5 min(y)-0.5 max(y)+0.5]);
    remaining = remaining - num_in_league;
end

repo: kafbox


text:
function str = opti_PreReqCheck(mode,verb)
%Check for VC++ 2010 / VC++ 2012 on x86 and x64 systems

if(nargin < 2), verb = 1; end
if(nargin < 1), mode = 'VS2012'; end

ROOTKEY = 'HKEY_LOCAL_MACHINE';
switch(lower(mode))
    case {'2012','vs2012'}
        %VC++ 2012 Keys
        KEY32 = 'SOFTWARE\Wow6432Node\Microsoft\DevDiv\VC\Servicing\11.0\RuntimeMinimum';
        KEY64 = 'SOFTWARE\Microsoft\DevDiv\VC\Servicing\11.0\RuntimeMinimum1';
        rtstr = 'VC++ 2012';
        rtredist86 = '';
        rtrunt86 = '';
        rtredist64 = '';
        rtrunt64 = '';
        instkey = 'Install';
    case {'2010','vs2010'}
        %VC++ 2010 Keys
        KEY32 = 'SOFTWARE\Microsoft\VisualStudio\10.0\VC\';
        KEY64 = 'SOFTWARE\Wow6432Node\Microsoft\VisualStudio\10.0\VC\';
        rtstr = 'VC++ 2010';
        rtredist86 = 'VCRedist\x86';
        rtrunt86 = 'Runtimes\x86';
        rtredist64 = 'VCRedist\x64';
        rtrunt64 = 'Runtimes\x64';
        instkey = 'Installed';
    otherwise
        error('Unknown Visual C++ / Intel C++ Version to check for! Only VC++ 2010 and 2012 supported');
end

mver = ver('MATLAB');
if(verb), fprintf('- Checking operating system...\n'); end
switch(mexext)
    case 'mexw32'
        if(verb), fprintf('MATLAB %s 32bit (Windows x86) detected\n',mver.Release); end
    case 'mexw64'
        if(verb), fprintf('MATLAB %s 64bit (Windows x64) detected\n',mver.Release); end
    otherwise
        error('OPTI Toolbox is compiled only for Windows systems - sorry!');
end

if(verb), fprintf('\n- Checking for required pre-requisites...\n'); end

str = [];
switch(mexext)
    %If mexw32 - could be 32bit windows or 64bit windows with 32bit matlab, must check both!
    case 'mexw32'
        %Check 32bit Location, Redist
        try
            if(winqueryreg(ROOTKEY,[KEY32 rtredist86],instkey))
                str = [rtstr ' x86 Redistributable (Win32) Found'];
            end
        catch %#ok<*CTCH>
            %Check 32bit Location, Runtime
            try
                if(winqueryreg(ROOTKEY,[KEY32 rtrunt86],instkey))
                    str = [rtstr ' x86 Runtime (Win32) Found'];
                end
            catch
                %Check 64bit Location, Redist
                try
                    if(winqueryreg(ROOTKEY,[KEY64 rtredist86],instkey))
                        str = [rtstr ' x86 Redistributable (Win64) Found'];
                    end
                catch
                    %Check 64bit Location, Runtime
                    try
                        if(winqueryreg(ROOTKEY,[KEY64 rtrunt86],instkey))
                            str = [rtstr ' x86 Redistributable (Win64) Found'];
                        end
                    catch
                        error(prereqerror('x86',rtstr));
                    end
                end
            end
        end
    case 'mexw64'
        %Check Redist
        try
            if(winqueryreg(ROOTKEY,[KEY64 rtredist64],instkey))
                str = [rtstr ' x64 Redistributable Found'];
            end
        catch
            %Check Runtime
            try
                if(winqueryreg(ROOTKEY,[KEY64 rtrunt64],instkey))
                    str = [rtstr ' x64 Runtime Found'];
                end
            catch
                error(prereqerror('x64',rtstr));
            end
        end
end

if(isempty(str))
    error('Could not determine the required pre-requisites, ensure you are running Windows x86 or x64!');
else
    if(verb), fprintf('%s\n',str); end
end

function str = prereqerror(build,rtstr)

switch(rtstr)
    case 'VC++ 2012'
        l = 'http://www.microsoft.com/en-us/download/details.aspx?id=30679';
        m = 'Microsoft';
    case 'VC++ 2010'
        switch(build)
            case 'x64'
                l = 'http://www.microsoft.com/download/en/details.aspx?id=14632';
            otherwise
                l = 'http://www.microsoft.com/download/en/details.aspx?id=5555';
        end
        m = 'Microsoft';
    case 'Intel C++ XE 2013'
        l = 'http://software.intel.com/en-us/articles/redistributable-libraries-for-intel-c-and-visual-fortran-composer-xe-2013-for-windows';
        m = 'Intel';
end

str = [sprintf(['Cannot find the %s ' rtstr ' ' build ' Redistributable!\n\n'],m)...
	   sprintf('Please close MATLAB and download the redistributable from %s:\n',m)...
       l...
       sprintf('\n\nThen reinstall OPTI Toolbox.')];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% lssep;
% separate close local minimizers
% and maybe add a few global points
%

nsep=0;

while nsep<nmin,
  % find intervals where the monotonicity behavior of both
  % adjacent intervals is opposite
  down=[flist(2:s)<flist(1:s-1)];
  sep= ( [1,1,down] & [0,up,0] & [down,1,1] );
  sep= sep | ( [1,1,up] & [0,down,0] & [up,1,1] );
  ind=find(sep);
  if isempty(ind), break; end;

  aa=0.5*(alist(ind)+alist(ind-1));	% interval midpoints
  if length(aa)>nloc,
    % select nloc best interval midpoints
    ff=min(flist(ind),flist(ind-1));
    [ff,ind]=sort(ff);
    aa=aa(ind(1:nloc));
  end;
  for alp=aa,
    if prt>2, disp(['separate minimizer at ',num2str(alp)]); end;
    % new function value
    falp=feval(func,data,x+alp*p);
    alist=[alist,alp];flist=[flist,falp];
    nsep=nsep+1;
    if nsep>=nmin, break; end;
  end;
  lssort;
end;

% instead of unnecessary separation, add some global points
for times=1:nmin-nsep,
  lsnew;			% extrapolation or split
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


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: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [u,res,iter] = deflationshifted(elem,A,b,M,isfreeNode,tol)
global ZH AH;
[HB, NL, level]      = HBstructure(elem,10);
% isfreeNode           = false(size(node,1),1);
% isfreeNode(freeNode) = 1;
[Pro]                = transferoperator(HB,NL,isfreeNode);
ZH                   = Pro{level-1};
AH                   = ZH'*A*ZH;
%invAH = inv(AH);
%P = (speye(Ni)-AD(freeNode,freeNode)*ZH*invAH*ZH');
u1 = ZH*(AH\(ZH'*b));
[u,~,res,iter]  = Pgmres_deflation(A,b,[],tol,size(b,1),M);
iter = max(iter);
u2 = u - ZH*(AH\(ZH'*(A*u)));
u0 = u1+u2;
u  = u0;
end

repo: ifem


text:
function t = s2t(s, epsilon);

% T = s2t(S, EPSILON)
%
% Scattering to Transmission Transformation
%
% since 2 is the only allowed number of ports for a T-matrix,
% S, T are matrices of size [2,2,F],
% where F is the number of frequencies
%
% EPSILON is used to produce an approximation of the corresponding S-matrix
% in exceptional cases
%
% 3-d version 4 nov 1999

if nargin < 2 epsilon = 1e-14; end;
[n,i] = min(abs(s(2,1,:)));
while n <= epsilon
    s(2,1,i) = s(2,1,i)+rand*epsilon;
    [n,i] = min(abs(s(2,1,:)));
end;

t(1,1,:) = 1./s(2,1,:);
t(1,2,:) = -s(2,2,:)./s(2,1,:);
t(2,1,:) = s(1, 1,:)./s(2,1,:);
t(2,2,:) = s(1, 2,:) - s(1, 1,:).*s(2,2,:)./s(2,1,:);

author: Sable

repo: mcbench-benchmarks


text:
function [features, background] = extractBlackBlobs(im)
% EXTRACTIMAGEBLOBS  Find the dark blobs in the images

f = hessiandet(im) ;
ok = f(3,:) > 0.0006 ;
idx = sub2ind(size(im), round(f(2,ok)), round(f(1,ok))) ;
features = false([size(im,1) size(im,2)]) ;
features(idx) = true ;
background = ~imdilate(features, strel('disk', 5, 0)) ;

b = 5 ;
features([1:b, end-b:end],:) = 0 ;
features(:,[1:b, end-b:end]) = 0 ;
background([1:b, end-b:end],:) = 0 ;
background(:,[1:b, end-b:end]) = 0 ;

function f = hessiandet(im)
% HESSIANDET  Basic Hessian detector
%   F = HESSIANDET(IM) runs a basic implementation of the Hessian
%   detector on the gray-scale image IM.

ims = imsmooth(im,2.5) ;

d2 = [1 -2 1];
d = [-1 0 1]/2 ;
im11 = conv2(1,d2,ims,'same') ;
im22 = conv2(d2,1,ims,'same') ;
im12 = conv2(d,d,ims,'same') ;
score = im11.*im22 - im12.*im12 ;

points = (islocalmax(score,1) .* islocalmax(score,2)) + ...
         (islocalmin(score,1) .* islocalmin(score,2)) ;
points = points .* (abs(score) > 0.0006) ;

scores = score(find(points)) ;
[i,j] = find(points) ;
f = [j(:),i(:),scores(:)]' ;

if 0
  clf ;subplot(1,2,1); imagesc(im) ; hold on ; colormap gray ;
  plot(f(1,:),f(2,:),'r.') ;
  subplot(1,2,2) ; imagesc(abs(score).^.25) ;
  keyboard
end

function m = islocalmax(x,dim)
m  = (circshift(x,1,dim) < x) & (circshift(x,-1,dim) < x) ;

function m = islocalmin(x,dim)
m = (circshift(x,1,dim) > x) & (circshift(x,-1,dim) > x) ;

author: vedaldi

repo: practical-cnn


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Q=SELECT_CLOSEST_CONFIGURATION(ROBOT, Qinv, CONF)
% Returns the joint coordinates Q that are closest with the axes configuration
% vector CONF, given a set of solutions of the inverse kinematic problem Qinv.
% For 6DOF or less manipulators, the variable CONF={CF1, CF4, CF6, CFX}
% specifies univoquely only one of the solutions. The first axis has
% priority over the rest.
%
% See also:
%   COMPUTE_CONFIGURATION, SELECT_CONFIGURATION , GET_CONF_DATA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function q=select_closest_configuration(robot, qinv, conf)

q=qinv(:,1);

distance=64;

for i=1:size(qinv,2),
    confi=compute_configuration(robot, qinv(:,i));

    d = abs(conf(1)-confi(1)) + abs(conf(2)-confi(2)) + abs(conf(3)-confi(3));

    if d < distance
       q =  qinv(:,i);
       distance=d;
    end
end

if d ~= 0
    disp('WARNING: RAPID/select_closest_configuration: No solutions complies with the specified configuration ');
end

repo: ARTE


text:
function vGlobal=getGlobalVectors(vLocal,uList)
%%GETGLOBALVECTORS Change a collection of local vectors into global
%                  vectors using the local coordinate axes. This multiplies
%                  the components of the vectors by the corresponding
%                  coordinate axis vectors.
%
%INPUTS: vLocal A numDimsXN matrix of N local vectors that are to be
%               converted.
%         uList A numDimsXnumDimsXN matrix of orthonormal unit coordinate
%               axes that are associated with the local coordinate system
%               in which the vectors are expressed. If all N vectors use
%               the same local coordinate system, then a numDimsxnumDimsx1
%               matrix can be passed instead.
%
%OUTPUTS:  vGlobal The vectors in the global coordinate system.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    numVel=size(vLocal,2);

    if(size(uList,3)==1)
        uList=repmat(uList,[1,1,numVel]);
    end

    numDims=size(vLocal,1);
    vGlobal=zeros(numDims,numVel);
    for curVel=1:numVel
        for curDim=1:numDims
            vGlobal(:,curVel)=vGlobal(:,curVel)+vLocal(curDim,curVel)*uList(:,curDim,curVel);
        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 comp_nmf_online_algorithms_test()
% Test file for stochastic/online NMF.
%
% Created by H.Kasai on Apr. 17, 2018

    clc;
    clear;
    close all;

    %% set options
    max_epoch                       = 100;
    verbose                         = 2;
    batch_size                      = 50;
    permute_on                      = true;
    svrmu_inneriter                 = 5;
    max_h_repeat                    = 5;
    calc_sol                        = 1;
    lambda                          = 1;

    dataset                         = 'CBCL';
    %dataset                         = 'COIL';
    %dataset                         = 'Synthetic';

    F                               = 100;
    K                               = 10;
    N                               = 500;
    noise_level                     = 0.0;
    outlier_rho                     = 0.0; % 0.9

    dic_display_flag                = 1;
    plot_flag                       = 1;
    clustering_flag                 = 1;
    denoise_display_flag            = 1;

    %% algorithms
    if outlier_rho == 0
        % batch
        nmf_hals_flag                   = 0;
        nmf_acc_hals_flag               = 0;
        rnmf_flag                       = 0;

        % spg
        spg_nmf_flag                    = 0;
        spg_acc_nmf_flag                = 0;
        apg_nmf_flag                    = 0;
        spg_ls_nmf_flag                 = 0;
        spg_precon_ls_nmf_flag          = 0;

        % online mu
        ronmf_flag                      = 0;
        onmf_flag                       = 0;
        onmf_acc_flag                   = 0;
        inmf_flag                       = 0;
        inmf_online_flag                = 0;
        asag_mu_nmf_flag                = 0;

        % smu
        smu_nmf_flag                    = 0;
        smu_acc_nmf_flag                = 0;
        smu_ls_nmf_flag                 = 0;

        % svrmu
        svrmu_nmf_flag                  = 0;
        svrmu_acc_nmf_flag              = 1;
        svrmu_ls_nmf_flag               = 0;
        svrmu_ls_acc_nmf_flag           = 0;
        svrmu_precon_ls_nmf_flag        = 0;
        svrmu_precon_ls_acc_nmf_flag    = 0;
        rsvrmu_acc_nmf_flag             = 0;
        svrmu_acc_adaptive_nmf_flag     = 0;
    else
        % batch
        nmf_hals_flag                   = 0;
        nmf_acc_hals_flag               = 0;
        rnmf_flag                       = 1;

        % spg
        spg_nmf_flag                    = 0;
        spg_acc_nmf_flag                = 0;
        apg_nmf_flag                    = 0;
        spg_ls_nmf_flag                 = 0;
        spg_precon_ls_nmf_flag          = 0;

        % online mu
        ronmf_flag                      = 1;
        onmf_flag                       = 0;
        onmf_acc_flag                   = 0;
        inmf_flag                       = 0;
        inmf_online_flag                = 0;
        asag_mu_nmf_flag                = 0;

        % smu
        smu_nmf_flag                    = 0;
        smu_acc_nmf_flag                = 0;
        smu_ls_nmf_flag                 = 0;

        % svrmu
        svrmu_nmf_flag                  = 0;
        svrmu_acc_nmf_flag              = 0;
        svrmu_ls_nmf_flag               = 0;
        svrmu_ls_acc_nmf_flag           = 0;
        svrmu_precon_ls_nmf_flag        = 0;
        svrmu_precon_ls_acc_nmf_flag    = 0;
        rsvrmu_acc_nmf_flag             = 1;
        svrmu_acc_adaptive_nmf_flag     = 0;
    end

    %% generate/load data
    fprintf('Loading data [%s] ...', dataset);
    if strcmp(dataset, 'PIE') || strcmp(dataset, 'AR') || strcmp(dataset, 'COIL') || strcmp(dataset, 'CBCL')

        if strcmp(dataset, 'PIE')
            img_in_dim  = [32, 32];
            img_out_dim = img_in_dim;
            batch_size  = 100; % up to N
            dic_display_dim = [4, 5];

        elseif strcmp(dataset, 'AR')
            img_in_dim  = [28, 20];
            img_out_dim = [28, 28];
            batch_size  = 100; % up to N
            dic_display_dim = [10, 12];

        elseif strcmp(dataset, 'COIL')
            img_in_dim  = [32, 32];
            img_out_dim = img_in_dim;
            batch_size  = 100; % up to N
            dic_display_dim = [4, 5];

        elseif strcmp(dataset, 'CBCL')
            img_in_dim  = [19, 19];
            img_out_dim = img_in_dim;
            batch_size  = 100; % up to N
            dic_display_dim = [7, 7];

        end

        [N, F, K, Vo, V, ~, classes] = load_dataset(dataset, '../../../data', outlier_rho);

        oriimg_display_flag = 1;

        if isempty(classes)
            clustering_flag = 0;
        end

    elseif strcmp(dataset, 'Synthetic')
        [Vo, Ro, No, V] = generate_syntheticdata(F, N, K, noise_level, outlier_rho);

        oriimg_display_flag     = 0;
        clustering_flag         = 0;
        denoise_display_flag    = 0;
        dic_display_flag        = 0;
    end

    dim = F * N;
    fprintf(' done\n');

    %% set initial data
    if 0
        x_init.W = rand(F, K);
        x_init.H = rand(K, N);
    elseif 1
        [W_init, H_init] = Qiao_SVD_Init(V, K);
        x_init.W = W_init;
        x_init.H = H_init;
    else
        [W_init, H_init] = NNDSVD(V, K, 0);
        x_init.W = W_init;
        x_init.H = H_init;
    end

    if outlier_rho == 0
        x_init.R = zeros(F, N);
    else
        x_init.R = rand(F, N);
    end

    %% calc solution
    if calc_sol

        clear options;
        options.max_epoch = max_epoch;
        options.x_init = x_init;
        options.verbose = 0;
        options.verbose = verbose;

        if ~outlier_rho
            fprintf('Calculating f_opt by HALS ...\n');
            options.alg = 'hals';
            [w_sol, infos_sol] = als_nmf(V, K, options);
        else
            fprintf('Calculating f_opt by R-NMF ...\n');
            options.lambda = lambda;
            [w_sol, infos_sol] = robust_mu_nmf(V, K, options);
        end

        f_opt = infos_sol.cost(end);
        fprintf('Done.. f_opt: %.16e\n', f_opt);
    else
        f_opt = -Inf;
    end

    %% execute algorithms
    names = cell(1);
    sols = cell(1);
    infos = cell(1);
    costs = cell(1);
    alg_idx = 0;

    %% Batch
    if nmf_hals_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.max_epoch = max_epoch;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.alg = 'hals';

        [w_nmf_hals, infos_nmf_hals] = als_nmf(V, K, options);

        names{alg_idx} = 'NMF HALS';
        sols{alg_idx} = w_nmf_hals;
        infos{alg_idx} = infos_nmf_hals;
        costs{alg_idx} = nmf_cost(Vo, w_nmf_hals.W, w_nmf_hals.H, zeros(F, N)) * 2 / dim;
    end

    if nmf_acc_hals_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.max_epoch = max_epoch;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.alg = 'acc_hals';

        [w_nmf_acc_hals, infos_nmf_acc_hals] = als_nmf(V, K, options);

        names{alg_idx} = 'NMF ACC HALS';
        sols{alg_idx} = w_nmf_acc_hals;
        infos{alg_idx} = infos_nmf_acc_hals;
        costs{alg_idx} = nmf_cost(Vo, w_nmf_acc_hals.W, w_nmf_acc_hals.H, zeros(F, N)) * 2 / dim;
    end

    if rnmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.max_epoch = max_epoch;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;

        [w_rnmf, infos_rnmf] = robust_mu_nmf(V, K, options);

        names{alg_idx} = 'R-NMF(Batch)';
        sols{alg_idx} = w_rnmf;
        infos{alg_idx} = infos_rnmf;
        costs{alg_idx} = nmf_cost(Vo, w_rnmf.W, w_rnmf.H, zeros(F, N)) * 2 / dim;
    end

    %% Online MU
    if ronmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;

        [w_ronmf, infos_ronmf] = robust_online_mu_nmf(V, K, options);

        names{alg_idx} = 'R-ONMF';
        sols{alg_idx} = w_ronmf;
        infos{alg_idx} = infos_ronmf;
        costs{alg_idx} = nmf_cost(Vo, w_ronmf.W, w_ronmf.H, zeros(F, N)) * 2 / dim;
    end

    if onmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        %options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;

        [w_onmf, infos_onmf] = online_mu_nmf(V, K, options);

        names{alg_idx} = 'ONMF';
        sols{alg_idx} = w_onmf;
        infos{alg_idx} = infos_onmf;
        costs{alg_idx} = nmf_cost(Vo, w_onmf.W, w_onmf.H, zeros(F, N)) * 2 / dim;
    end

    if onmf_acc_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;

        [w_onmf_acc, infos_onmf_acc] = acc_online_mu_nmf(V, K, options);

        names{alg_idx} = 'ONMF-ACC';
        sols{alg_idx} = w_onmf_acc;
        infos{alg_idx} = infos_onmf_acc;
        costs{alg_idx} = nmf_cost(Vo, w_onmf_acc.W, w_onmf_acc.H, zeros(F, N)) * 2 / dim;
    end

    if inmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.max_epoch = max_epoch;
        options.online = 0; % online mode
        options.x_init = x_init;
        options.verbose = 1;
        options.f_opt = f_opt;
        options.batch_size = batch_size;

        [w_inmf, infos_inmf] = incremental_mu_nmf(V, K, options);

        names{alg_idx} = 'INMF';
        sols{alg_idx} = w_inmf;
        infos{alg_idx} = infos_inmf;
        costs{alg_idx} = nmf_cost(Vo, w_inmf.W, w_inmf.H, zeros(F, N)) * 2 / dim;
    end

    if inmf_online_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.max_epoch = 1;
        options.online = 1; % online mode
        options.max_inneriter = 100;
        options.x_init = x_init;
        options.verbose = 2;
        options.f_opt = f_opt;
        options.batch_size = batch_size;

        [w_online_inmf, infos_online_inmf] = incremental_mu_nmf(V, K, options);

        names{alg_idx} = 'INMF (Online)';
        sols{alg_idx} = w_online_inmf;
        infos{alg_idx} = infos_online_inmf;
        costs{alg_idx} = nmf_cost(Vo, w_online_inmf.W, w_online_inmf.H, zeros(F, N)) * 2 / dim;
    end

    %%
    if asag_mu_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = 1;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.permute_on = permute_on;
        options.lambda = 0.5;
        options.f_opt = f_opt;

        [w_asag_mu_nmf, infos_asag_mu_nmf] = asag_mu_nmf(V, K, options);

        names{alg_idx} = 'ASAG-MU';
        sols{alg_idx} = w_asag_mu_nmf;
        infos{alg_idx} = infos_asag_mu_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_asag_mu_nmf.W, w_asag_mu_nmf.H, zeros(F, N)) * 2 / dim;
    end

    %% SMU
    if smu_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 0;

        [w_smu_nmf, infos_smu_nmf] = smu_nmf(V, K, options);

        names{alg_idx} = 'SMU';
        sols{alg_idx} = w_smu_nmf;
        infos{alg_idx} = infos_smu_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_smu_nmf.W, w_smu_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if smu_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;

        [w_smu_nmf_acc, infos_smu_nmf_acc] = smu_nmf(V, K, options);

        names{alg_idx} = 'SMU-ACC';
        sols{alg_idx} = w_smu_nmf_acc;
        infos{alg_idx} = infos_smu_nmf_acc;
        costs{alg_idx} = nmf_cost(Vo, w_smu_nmf_acc.W, w_smu_nmf_acc.H, zeros(F, N)) * 2 / dim;
    end

    if smu_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.ls = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;

        [w_smu_nmf_acc, infos_smu_nmf_acc] = smu_nmf(V, K, options);

        names{alg_idx} = 'SMU-LS';
        sols{alg_idx} = w_smu_nmf_acc;
        infos{alg_idx} = infos_smu_nmf_acc;
        costs{alg_idx} = nmf_cost(Vo, w_smu_nmf_acc.W, w_smu_nmf_acc.H, zeros(F, N)) * 2 / dim;
    end

    if smu_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.ls = 1;
        options.precon = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;

        [w_smu_nmf_acc, infos_smu_nmf_acc] = smu_nmf(V, K, options);

        names{alg_idx} = 'SMU-Precon-LS';
        sols{alg_idx} = w_smu_nmf_acc;
        infos{alg_idx} = infos_smu_nmf_acc;
        costs{alg_idx} = nmf_cost(Vo, w_smu_nmf_acc.W, w_smu_nmf_acc.H, zeros(F, N)) * 2 / dim;
    end

    %% SPG
    if spg_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 0;
        options.ls = 0;

        [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, K, options);

        names{alg_idx} = 'SPG';
        sols{alg_idx} = w_spg_nmf;
        infos{alg_idx} = infos_spg_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_spg_nmf.W, w_spg_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if spg_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 1;
        options.ls = 0;
        options.h_repeat = max_h_repeat;

        [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, K, options);

        names{alg_idx} = 'SPG-ACC';
        sols{alg_idx} = w_spg_nmf;
        infos{alg_idx} = infos_spg_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_spg_nmf.W, w_spg_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if apg_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 0;
        options.ls = 0;
        %options.h_repeat = max_h_repeat;
        options.W_sub_mode = 'Nesterov';

        [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, K, options);

        names{alg_idx} = 'APG';
        sols{alg_idx} = w_spg_nmf;
        infos{alg_idx} = infos_spg_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_spg_nmf.W, w_spg_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if spg_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 0;
        options.ls = 1;

        [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, K, options);

        names{alg_idx} = 'SPG-LS';
        sols{alg_idx} = w_spg_nmf;
        infos{alg_idx} = infos_spg_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_spg_nmf.W, w_spg_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if spg_precon_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;
        options.max_epoch = max_epoch;
        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.f_opt = f_opt;
        options.accel = 0;
        options.ls = 1;
        options.precon = 1;

        [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, K, options);

        names{alg_idx} = 'SPG-Precon-LS';
        sols{alg_idx} = w_spg_nmf;
        infos{alg_idx} = infos_spg_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_spg_nmf.W, w_spg_nmf.H, zeros(F, N)) * 2 / dim;
    end

    %% SVRMU
    if svrmu_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.x_init.R = zeros(F, N); % enfore non-robust mode
        options.verbose = verbose;
        options.permute_on = permute_on;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 0;
        options.ratio = 1; % 1: original 1<: adaptive
        options.permute_on = permute_on;
        options.robust = false;
        %options.stepsize_ratio = 1;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        %[w_svrmu_nmf, infos_svrmu_nmf] = svrmu_nmf_old(V, K, options);
        [w_svrmu_nmf, infos_svrmu_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU';
        sols{alg_idx} = w_svrmu_nmf;
        infos{alg_idx} = infos_svrmu_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_nmf.W, w_svrmu_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-ACC';
        sols{alg_idx} = w_svrmu_acc_nmf;
        infos{alg_idx} = infos_svrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_nmf.W, w_svrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 0;
        options.ls = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-LS';
        sols{alg_idx} = w_svrmu_acc_nmf;
        infos{alg_idx} = infos_svrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_nmf.W, w_svrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_ls_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 1;
        options.ls = 1;
        options.rep_mode = 'fix';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-LS-ACC';
        sols{alg_idx} = w_svrmu_acc_nmf;
        infos{alg_idx} = infos_svrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_nmf.W, w_svrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_precon_ls_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 0;
        options.ls = 1;
        options.precon = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-Precon-LS';
        sols{alg_idx} = w_svrmu_acc_nmf;
        infos{alg_idx} = infos_svrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_nmf.W, w_svrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_precon_ls_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 1;
        options.ls = 1;
        options.precon = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-Precon-LS-ACC';
        sols{alg_idx} = w_svrmu_acc_nmf;
        infos{alg_idx} = infos_svrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_nmf.W, w_svrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if svrmu_acc_adaptive_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.permute_on = permute_on;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 1;
        %options.rep_mode = 'fix';
        options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.ratio = 0.3; % 1: original 1<: adaptive
        options.permute_on = permute_on;
        options.robust = false;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_svrmu_acc_adaptive_nmf, infos_svrmu_acc_adaptive_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'SVRMU-ACC-Adaptive';
        sols{alg_idx} = w_svrmu_acc_adaptive_nmf;
        infos{alg_idx} = infos_svrmu_acc_adaptive_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_svrmu_acc_adaptive_nmf.W, w_svrmu_acc_adaptive_nmf.H, zeros(F, N)) * 2 / dim;
    end

    if rsvrmu_acc_nmf_flag
        alg_idx = alg_idx + 1;
        clear options;
        options.lambda = lambda;

        options.batch_size = batch_size;
        options.x_init = x_init;
        options.verbose = verbose;
        options.permute_on = permute_on;
        options.fast_calc = 0;
        options.f_opt = f_opt;
        options.repeat_inneriter = svrmu_inneriter;
        options.accel = 0;
        options.ls = 1;
        options.rep_mode = 'fix';
        %options.rep_mode = 'adaptive';
        options.w_repeat = 1;
        options.h_repeat = max_h_repeat;
        options.ratio = 1;
        options.permute_on = permute_on;
        options.robust = true;

        options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));

        [w_rsvrmu_acc_nmf, infos_rsvrmu_acc_nmf] = svrmu_nmf(V, K, options);

        names{alg_idx} = 'R-SVRMU-LS';
        sols{alg_idx} = w_rsvrmu_acc_nmf;
        infos{alg_idx} = infos_rsvrmu_acc_nmf;
        costs{alg_idx} = nmf_cost(Vo, w_rsvrmu_acc_nmf.W, w_rsvrmu_acc_nmf.H, zeros(F, N)) * 2 / dim;
    end

    alg_total = alg_idx;

    %% plot
    if plot_flag
        %display_graph_icassp('epoch','cost', names, sols, infos);
        %display_graph_icassp('grad_calc_count','cost', names, sols, infos);
        %display_graph_icassp('time','cost', names, sols, infos);

        display_graph('numofgrad','cost', names, sols, infos);
        display_graph('time','cost', names, sols, infos);
        if f_opt ~= -Inf
            display_graph('numofgrad','optimality_gap', names, sols, infos);
            display_graph('time','optimality_gap', names, sols, infos);
        end
    end

    % display original dic
    if oriimg_display_flag && ~isempty(classes)
        for ii = 1 : K
            k_index = find(classes==ii);
            col_img = Vo(:,k_index(1));
            col_img_reshape = reshape(col_img, img_in_dim);
            col_img_resize = imresize(col_img_reshape, img_out_dim);
            img(:,ii) = col_img_resize(:);
        end
        plot_dictionnary(img, [], dic_display_dim);
    end

    clear results;

    for alg_idx = 1 : alg_total

        mse = calc_mse(Vo, sols{alg_idx}.W, sols{alg_idx}.H);
        psnr = calc_psnr(Vo, sols{alg_idx}.W, sols{alg_idx}.H);
        fprintf('# %s:\t MSE:%e,\t PSNR:%e,\t time:%e [sec]\n', names{alg_idx}, mse, psnr, infos{alg_idx}.time(end));

        if clustering_flag
            cluster_num = length(unique(classes));
            Coeff = sols{alg_idx}.H;
            M = Coeff';
            num_samples = size(M,1);
            perm_idx = randperm(num_samples);
            center = zeros(cluster_num, K);
            for i = 1 : cluster_num
                center(i,:) = M(perm_idx(i),:);
            end

            % do k-menas
            k_means_result = k_means(M, center, cluster_num);

            purity = calc_purity(classes, k_means_result);
            nmi = calc_nmi(classes, k_means_result);

            fprintf('# %s:\t Purity:%e,\t NMI: %e\n', names{alg_idx}, purity, nmi);
        end

    end

    if dic_display_flag
        for alg_idx = 1 : alg_total
            clear img;
            for ii = 1 : K
                col_img = sols{alg_idx}.W(:,ii);
                col_img_reshape = reshape(col_img, img_in_dim);
                col_img_resize = imresize(col_img_reshape, img_out_dim);
                img(:,ii) = col_img_resize(:);
            end

            plot_dictionnary(img, [], dic_display_dim);

        end
    end

    if denoise_display_flag && ~isempty(classes)
        for alg_idx = 1 : alg_total

            Rec = sols{alg_idx}.W * sols{alg_idx}.H;

            for ii = 1 : K
                k_index = find(classes==ii);
                col_img = Rec(:,k_index(1));
                col_img_reshape = reshape(col_img, img_in_dim);
                col_img_resize = imresize(col_img_reshape, img_out_dim);
                img(:,ii) = col_img_resize(:);
            end

            plot_dictionnary(img, [], dic_display_dim);

        end
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
%  This scriptfile ask for several input parameters that can be setup
%  at the beginning of each session. The default values are the
%  extrema in the catalog
%
%a = org;        % resets the main catalogue "a" to initial state

%special version of the inpu script for Mapseis, it does not open a
%parameter selection window is the catalog should be already filtered.

report_this_filefun(mfilename('fullpath'));

%  default values
t0b = min(a.Date);
teb = max(a.Date);
tdiff = (teb - t0b)*365;

if ~exist('par1','var')
    %  if tdiff>10                 %select bin length respective to time in catalog
    %     par1 = ceil(tdiff/100);
    %  elseif tdiff<=10 & tdiff>1
    %     par1 = 0.1;
    %  elseif tdiff<=1
    %     par1 = 0.01;
    %  end
    par1 = 30;
end

minmag = max(a.Magnitude) -0.2;
dep1 = 0.3*max(a.Depth);
dep2 = 0.6*max(a.Depth);
dep3 = max(a.Depth);
minti = min(a.Date);
maxti  = max(a.Date);
minma = min(a.Magnitude);
maxma = max(a.Magnitude);
mindep = min(a.Depth);
maxdep = max(a.Depth);

think;
sele_sub;

author: CelsoReyes


text:
function diplayFinalModelsCT_HN(pathWORK,outcomes)
% -------------------------------------------------------------------------
% function diplayFinalModelsCT_HN(pathWORK,outcomes)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function displays the final models for all analyzed outcomes in the
% HN study for the 'CT' feature set only, in terms of: multivariable model
% response, model variables, and prediction performance estimation. See
% ref.[1] for more details. Goal: comparison with the "Radiomics signature"
% of (Aerts et al., Nat Commun, 2014)
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the HN WORKSPACE directory.
% - outcomes: Structure specifying the status (1 or 0) for different
%             outcomes in HN cancer. Contains: outcomes.Failure,
%             outcomes.Locoregional, outcomes.Distant, outcomes.Death. See
%             ref.[1] for more details.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
nameOutcomes = fieldnames(outcomes);
nOutcomes = length(nameOutcomes);

for o = 1:nOutcomes
    cd([pathWORK,'/FINAL_MODELS/',nameOutcomes{o},'/CTonly'])
    load('finalModel'), load('coeff'), load('modelCI'), load('response') % Variables 'finalModel', 'coeff',  'modelCI', and 'response' now present in MATLAB workspace
    order = size(finalModel.Data,2);
    plotSigmoidalResponse(response,outcomes.(nameOutcomes{o}),modelCI,upper(nameOutcomes{o}))
    fprintf(['\n\n --> THE FINAL MULTIVARIABLE MODEL (CT FEATURE SET) FOR ''%s'' OUTCOME IS:\n\n'...
             '               g(x) =               \n'],upper(nameOutcomes{o}))
    for i = 1:order
        fprintf([num2str(coeff(i)),' X ',finalModel.Name{i},'\n'])
        fprintf('                    +               \n')
    end
    fprintf(['                   ',num2str(coeff(end)),'\n'])
    fprintf('\nWITH CORRESPONDING PREDICTION PERFORMANCE ESTIMATION:\n')
    fprintf(['AUC = ',num2str(roundsd(finalModel.AUC632,ceil(log10(finalModel.AUC632/roundsd(finalModel.SE_AUC632,1))))),' ± ',num2str(roundsd(finalModel.SE_AUC632,1)),'\n'])
    fprintf(['Sensitivity = ',num2str(roundsd(finalModel.Sensitivity632,ceil(log10(finalModel.Sensitivity632/roundsd(finalModel.SE_Sensitivity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Sensitivity632,1)),'\n'])
    fprintf(['Specificity = ',num2str(roundsd(finalModel.Specificity632,ceil(log10(finalModel.Specificity632/roundsd(finalModel.SE_Specificity632,1))))),' ± ',num2str(roundsd(finalModel.SE_Specificity632,1)),'\n'])
    fprintf('\n\n')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function view = recomputeSSImage(view,numGrays,numColors,clipMode)
%
% function view = recomputeSSImage(view,numGrays,numColors,clipMode)
%
% rmk 09/17/98 based on recomputeImage

% Initialize images
anatIm=[];

% Get anatClip and overlayClip from sliders
anatClip = getAnatClip(view);

% Get anatomy image
anatIm = view.anat;

% Rescale anatIm to [1:numGrays], anatClip determines the range
% of anatomy values that gets mapped to the available grayscales.
% If anatClip=[0,1] then there is no clipping and the entire
% range anatomy values is scaled to the range of available gray
% scales.
minVal = min(anatIm(:));
maxVal = max(anatIm(:));
anatClipMin = min(anatClip)*(maxVal-minVal) + minVal;
anatClipMax = max(anatClip)*(maxVal-minVal) + minVal;
if (length(anatIm(:))~=0)
    anatIm=rescale2(anatIm,[anatClipMin,anatClipMax],[1,numGrays]);
end

% Replace NaNs
indices = find(isnan(anatIm));
anatIm(indices) = 1;

% Finally, set the view.ui.image field
view.ui.image = uint8(anatIm-1);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [result_features]=back_find_high_node(W,C,nROI,w,midw_lasso,IDX,opt_t)

% This function aims to find features used in the classification, suitable
% for the dHOFC network construction method. And it is also
% suitable for the LOOCV or the 10-fold cross validation. The found
% features can be applied on the visualization software to present the
% important high-order nodes (a set of links in each high-order node).
% Input:
%         W: window length
%         C: number of cluster
%         nROI: number of ROIs
%         w: weight of each selected features;
%         midw_lasso: the feature selection indexes of lasso;
%         IDX: the index in clustering;
%         opt_t: the selected feat in each outer LOOCV;
% Output:
%        result_features: cell array, consisting of occurrence of the features,
%        node number, the cluster matrix and weight of the cluster;

% Written by Zhen Zhou, zzstefan@email.unc.edu
% IDEA lab, https://www.med.unc.edu/bric/ideagroup
% Department of Radiology and BRIC, University of North Carolina at Chapel Hill

%
% load dHOFC_loocv_middle.mat;
%load dHOFC_kfold_middle.mat;
fprintf('Begin contributing feature identification\n');

for i =1:size(midw_lasso,1)
    for j=1:size(midw_lasso,2)
        which_C=ceil(opt_t(i,j)/length(W));
        which_W=mod(opt_t(i,j),length(W));
        if which_W ==0
            which_W=length(W);
        end
        index{i,j}=IDX{which_W,which_C};
        high_index{i,j}=find(midw_lasso{i,j}~=0);
    end
end
%% index means the corresponding cluster number of each of the 6670 lines in the original low-order FC;
%% high_index ,each LOOCV ,there is one combination of hyper-parameter constructed brain network being choosen,
%% find the cluster number after feature selection

index=index(:);
high_index=high_index(:);
for i=1:length(high_index)
    for j=1:length(high_index{i})
        tmp{i,j}=find(index{i}==high_index{i}(j));
    end
    %length_test(i)=length(find(~cellfun(@isempty,tmp(i,:))));
end

% B=unique(length_test);
% for i=1:length(B)
%     A(i)=length(find(length_test==B(i)));
% end

%% find the frequency of each element in the cell array
tmp=cellfun(@(x) x',tmp,'UniformOutput',false);
temp=tmp(:);

[a1,b,c]=unique(cellfun(@char,temp,'un',0));
lo=histc(c,1:max(c));
loo=lo(:)>1;
out=[temp(b(loo)),num2cell(lo(loo))];
out(1,:)=[];
new_out=out;

%% find the weight of each cluster and the averaged accuracy would be the mean weight for each cluster.
for k=1:length(new_out)
    weight=[];
    test{k}=cellfun(@(x)isequal(x,new_out{k}),tmp,'un',0);
    [location{k}(:,1),location{k}(:,2)]=find(cell2mat(test{k}));
    for j=1:length(location{k})
        weight=[weight,w{location{k}(j,1)}(location{k}(j,2))];
    end
    W(k)=mean(weight);
end

for j=1:length(new_out)
    inner_tmp=new_out{j};
    for i=1:length(inner_tmp)
        [first{j}(i,1),first{j}(i,2)]=find_elements(116,new_out{j,1}(i));
    end
    new_out(j,3)=first(j);
end

%%according to the connection between ROI, create the 0-1 matrix for the
%%convenience of drawing in brainnet viewer

%Matrix=create_matrix(first);

for m=1:length(first)
    node_matrix{m}=zeros(nROI,nROI);
    for n=1:length(first{m})
        node_matrix{m}(first{m}(n,1),first{m}(n,2))=1;
        node_matrix{m}(first{m}(n,2),first{m}(n,1))=1;
    end
end

% nROI=116;
% node_matrix=eye(nROI,nROI);
% for m=1:length(test)
%
%         node_matrix(test{m}(:,1),test{m}(:,2))=1;
%         node_matrix(test{m}(:,2),test{m}(:,1))=1;
%
% end

result_features(:,1)=new_out(:,2);
result_features(:,2)=new_out(:,3);
result_features(:,3)=node_matrix';
result_features(:,4)=num2cell(W);
%result_features=sortrows(result_features,1,'descend');%% sort according to the frequency of each cluster ,descend
result_features=sortrows(result_features,-1);

fprintf('End contributing feature identification\n');

author: zzstefan

repo: BrainNetClass


text:
function varargout = xyy2xyz(varargin)
%XYY2XYZ Convert chromaticity coordinates to XYZ tristimulus values.
%   [X,Y,Z] = XYY2XYZ(x,y,Y) converts the chromaticity coordinates (x
%   and y) and the Y tristimulus value to the tristimulus values X, Y,
%   and Z. Y, x, y, X, Y, and Z are all arrays with the same size.
%
%   XYZ = XYY2XYZ(xyY) performs the conversion on a Px3 input matrix
%   where P(:,1) contains the x values, P(:,2) contains the y values,
%   and P(:,3) contains the Y values. The output matrix is also Px3.

%   Written by Steve Eddins to accompany Digital Image Processing Using
%   MATLAB, 3rd edition, Gatesmark Press, 2020,
%   http://imageprocessingplace.com.
%
%   Copyright 2019 The MathWorks, Inc.
%   License: https://github.com/mathworks/matlab-color-tools/blob/master/license.txt

if nargin == 3
   style = "separate";
   x = varargin{1};
   y = varargin{2};
   Y = varargin{3};
elseif nargin == 1
   style = "3-column";
   xyY = varargin{1};
   x = xyY(:,1);
   y = xyY(:,2);
   Y = xyY(:,3);
end

X = Y .* x ./ y;
Z = Y .* (1 - x - y) ./ y;

% Handle the numerical edge case where y is 0.
X(y == 0) = 0;
Z(y == 0) = 0;

if style == "3-column"
   varargout{1} = [X Y Z];
else
   varargout{1} = X;
   varargout{2} = Y;
   varargout{3} = Z;
end

author: dipum

repo: dipum-toolbox

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


text:
function roi = dtiRoiMakePlane(cornerCoords, roiName, roiColor)
%
% roi = dtiRoiMakePlane(cornerCoords, name, color)
%
% e.g.:
% roi = dtiRoiMakePlane([-60,-30,-60; 60,-30,80], 'coronal_plane', 'c')
%
% For a more complex sample, edit this function to see the code below the
% return.
%
% 2008.10.01 RFD: wrote it.

if(~exist('roiName','var') || isempty(roiName))
    roiName = 'plane';
end
if(~exist('roiColor','var') || isempty(roiColor))
    roiColor = 'r';
end
lc = min(cornerCoords);
uc = max(cornerCoords);
[x,y,z] = ndgrid([lc(1):uc(1)], [lc(2):uc(2)], [lc(3):uc(3)]);

roi = dtiNewRoi(roiName, roiColor, [x(:),y(:),z(:)]);

return;

% To run this on a bunch of subjects:
bd = '/biac3/wandell4/data/reading_longitude/dti_y4';
[subList,subCodes,subDirs] = findSubjects(fullfile(bd,'/*'),'dti06');
x = [-70 70];
z = [-60 80];
y = [-30 -31 -34 -30 ];
for(ii=1:numel(subList))
    roi = dtiRoiMakePlane([x(1),y(ii),z(1); x(2),y(ii),z(2)], 'coronal', 'c');
    roiFileName = fullfile(subDirs{ii},'ROIs','IPSproject','coronal');
    dtiWriteRoi(roi,roiFileName);
end


text:
classdef NoDirectiveChanged < MBeautifier.DirectiveChange
    methods
        function obj = NoDirectiveChanged()
            obj@MBeautifier.DirectiveChange(MBeautifier.DirectiveChangeType.NONE);
        end
    end
end

author: davidvarga

repo: MBeautifier


text:
function Population = archive(Population,N)
% Select feasible and non-dominated solutions by using SPEA2-CDP

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

    %% Select feasible solutions
    fIndex     = all(Population.cons <= 0,2);
    Population = Population(fIndex);
    if isempty(Population)
        return
    elseif size(Population,2) > N
        Fitness = CalFitness(Population.objs,Population.cons);
        Next = Fitness < 1;
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
        Population = Population(Next);
    end
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [proj_mat,D] = tomoproj2d(im,angles)
%TOMOPROJ2D   [proj_mat,D] = tomoproj2d(im,angles)
%   unit of angles: Degree;
%   projection direction: When angle == 0, X-ray passes through vertical
%   (up-down) direction.
%
%   Phymhan
%   02-Aug-2013 14:07:06

%Pad image
[im_pad,D] = impad(im);
%Calculate projection
num_proj = length(angles);
proj_mat = zeros(num_proj,D);
for k = 1:num_proj
    im_rot = imrotate(im_pad,-angles(k),'bilinear','crop');
    proj_mat(k,:) = sum(im_rot,1);
end
end


text:
function RightImpactConstraints(nlp, src, tar, bounds, varargin)
    plant = nlp.Plant;

    % no need to be time-continuous
    removeConstraint(nlp,'tContDomain');

    % first call the class method (calling impact model since it no longer
    % applies if we have a custom function)
    plant.rigidImpactConstraint(nlp, src, tar, bounds, varargin{:});

    % the relabeling/periodicity of joint coordiante is no longer valid
    % (this only affects position peridicity, velocity still applies)
    removeConstraint(nlp,'xDiscreteMapRightImpact');

    % Readding Periodicity (ignoring first 6 coordinates)
    R = plant.R;
    x = plant.States.x;
    xn = plant.States.xn;
    x_diff = R*x-xn;
    x_map = SymFunction(['xPartialDiscreteMap' plant.Name],x_diff(7:end),{x,xn});

    addNodeConstraint(nlp, x_map, {'x','xn'}, 'first', 0, 0, 'Linear');
end

author: ayonga

repo: frost-dev

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


text:
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function csd = CSD(lfp)

%CSD - Compute current source density.
%
%  USAGE
%
%    csd = CSD(lfp)
%
%    lfp            local field potential samples
%
%  SEE
%
%    See also PlotCSD.

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% 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.

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help CSD'' for details).');
end

t = lfp(:,1);
y = lfp(:,2:end);
y = y - repmat(mean(y),length(t),1);
d = -diff(y,2,2);
csd = [t d];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% change in body orientation
function [data,units] = compute_dtheta(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  if trx(fly).nframes <= 1,
    data{i} = [];
  else
    data{i} = modrange(diff(trx(fly).theta_mm,1,2),-pi,pi)./trx(fly).dt;
  end
end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
function [ fea, out ] = ex_spanner( varargin )
%EX_SPANNER STL geometry import and stress calculation of a spanner.
%
%   [ FEA, OUT ] = EX_SPANNER( VARARGIN ) Example to import a STL CAD
%   geometry and calculate displacements on a fixed spanner. The load
%   force may be distributed in the tangential load direction with the
%   force fraction parameter FRAC.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       E           scalar {190e3}         Modulus of elasticity [N/mm^2]
%       nu          scalar {0.29}          Poissons ratio
%       force       scalar {1000}          Load force [N]
%       frac        scalar {0}             Fraction of stress against pulling direction
%       sfun        string {sflag1}        Shape function for displacements
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'E',        190e3;
            'nu',       0.29;
            'force',    1000;
            'frac',     0;
            'sfun',     'sflag1';
            'iplot',    1;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});

% Define scaling factor m to mm.
USE_METERS = true;
s = double( ~USE_METERS + USE_METERS*1e-3 );

% Import (and scale) grid.
fea.sdim = {'x','y','z'};
fea.geom = impexp_stl( 'spanner1.stl', 'import', [], 'solid', 2, 'extend', 0 );
fea.grid = gridgen( fea, 'hmax', 5, 'gridgen', 'robust', 'fid', opt.fid );
fea.grid.p = fea.grid.p*s;

% Add linear elasticity physics mode and define material parameters.
fea = addphys( fea, @linearelasticity );
fea.phys.el.sfun            = { opt.sfun opt.sfun opt.sfun };
fea.phys.el.eqn.coef{1,end} = { opt.nu };
fea.phys.el.eqn.coef{2,end} = { opt.E/s^2 };

% Set all boundaries to no load per default.
n_bdr  = max(fea.grid.b(3,:));
bc_sel = cell(3,n_bdr);
[bc_sel{:}] = deal(0);

% Fix all displacements on mandible boundaries.
i_fix   = [13 14];
i_force = [3 1];
faxis   = [1 3];
[bc_sel{:,i_fix}] = deal(1);
fea.phys.el.bdr.coef{5} = bc_sel;

% Apply force for x > 140 mm.
force = opt.force/(6*s*80*s);
fea.phys.el.bdr.coef{7}{faxis(1),i_force(1)} = ['-',num2str((1-opt.frac)*force),'*(y>140*',num2str(s),')'];
fea.phys.el.bdr.coef{7}{faxis(2),i_force(2)} = [num2str(opt.frac*force),'*(y>140*',num2str(s),')'];

% Parse and solve problem.
fea = parsephys(fea);
fea = parseprob(fea);
fea.sol.u = solvestat( fea, 'fid', opt.fid, 'icub', 1+str2num(strrep(opt.sfun,'sflag','')) );

% Postprocessing.
if( opt.iplot>0 )
  subplot(1,2,1)
  postplot( fea, 'surfexpr', ['sqrt(u^2+v^2+w^2)/',num2str(s)] )
  view([30 20])
  title('Total displacement (mm)')

  subplot(1,2,2)
  DSCALE = 5;
  dp = zeros(size(fea.grid.p));
  for i=1:3
    dp(i,:) = DSCALE*evalexpr( fea.dvar{i}, fea.grid.p, fea );
  end
  fea_disp.grid   = fea.grid;
  fea_disp.grid.p = fea_disp.grid.p + dp;
  plotgrid( fea_disp )
  title(['Displacement plot'])
  view([30 20])
end

% Error checking.
u = fea.sol.u(unique(fea.eqn.dofm{1}(:)))/s;
v = fea.sol.u(unique(fea.eqn.ndof(1)+fea.eqn.dofm{2}(:)))/s;
w = fea.sol.u(unique(sum(fea.eqn.ndof(1:2))+fea.eqn.dofm{3}(:)))/s;
out.disp = sqrt( u.^2 + v.^2 + w.^2 );
out.pass = nan;
if( ~(got.frac || got.E || got.nu || got.force) )
  out.pass = abs( max(out.disp) - 2.5 )/2.5 < 0.1;
end

if( nargout==0 )
  clear fea out
else
  fea.grid.p = fea.grid.p/s;
end

author: precise-simulation

repo: featool-multiphysics


text:
function [oz] = cm32oz(cm3)
% Convert volume from cubic centimeters to US liquid ounces.
% Chad Greene 2012
oz = cm3*0.033814022701;


text:
function GUI = sessionGUI_loadROI(rois, viewer);
% Load a saved ROI into a mrViewer UI.
%
% GUI = sessionGUI_loadROI(<rois=get from GUI>, <mrViewer handle>);
%
% rois can be specified as the name of an ROI in the relevant ROI dir
% (inplane or volume, depending on the GUI's roi type), a path to an
% ROI file, an index into the ROI list being displayed, a handle
% to a uicontrol with the ROI list, or a cell array of several of these
% specifications. By default, it uses the GUI.controls.roi listbox.
%
% This function will load the ROIs into the current mrViewer UI specified
% by GUI.settings.viewer. An alternate viewer can be specified by passing
% the handle to the main UI figure. If no viewer is open, this will open
% one.
%
%
% ras, 07/06.
mrGlobals2;

if notDefined('rois'), rois = GUI.controls.roi; end
if notDefined('viewer')
    if isempty(GUI.viewers)
        viewer = sessionGUI_viewInplane;
	else
	    viewer = GUI.viewers(GUI.settings.viewer);
	end
end

% parse the format of the rois arg, getting it into the format
% of a cell array of file paths
if iscell(rois)
    % check that each entry is a path
    for i = 1:length(rois)
        rois{i} = sessionGUI_roiPath(rois{i});
    end

elseif isnumeric(rois) & all(mod(rois, 1)==0)
    % index into ROIs list
    for i = 1:length(rois)
        roiList{i} = sessionGUI_roiPath(rois(i));
    end
    rois = roiList;

elseif ishandle(rois)
    % get from uicontrol
    rois = get(rois, 'Value');
    for i = 1:length(rois)
        roiList{i} = sessionGUI_roiPath(rois(i));
    end
    rois = roiList;

elseif ischar(rois)
    rois = {sessionGUI_roiPath(rois)};

else
    error('Invalid ROI specification. ')

end

% Now that we have a cell array of paths, load each one into the current
% viewer:
ui = get(viewer, 'UserData');
for i = 1:length(rois)
    ui = mrViewLoad(ui, rois{i}, 'roi');
end

% select the last one in the list
mrViewROI('select', ui, length(ui.rois));

return


text:
function [V,net,genFlag] = TrainGrowingGasNet(V,temp1,net,scale,params,Problem,wholeObj,genFlag,zmin)

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

% This function is written by Qiqi Liu

    %% Parameters
    N = params.N;
    MaxIt = params.MaxIt;
    L = params.L;
    epsilon_b = params.epsilon_b;
    epsilon_n = params.epsilon_n;
    alpha = params.alpha;
    delta = params.delta;
    T = params.T;

    t = net.t;
    w = net.w;
    E = net.E;
    C = net.C;
    nx = net.nx;
    ageSumBefore = net.ageSumBefore;
    flag = net.flag;

    % select the corner solutions using GNG
    output=[];
    fm =[];
    C = net.C;
    w= net.w;
    for i = 1:size(w,1)
        neighbor = find(C(i,:)==1);

        ageSum(i,:) = sum(t(i,neighbor,:),2);
        if ageSum(i,:) == ageSumBefore(i,:)
            flag(i,:) = flag(i,:) + 1;
        end

    end

    cw = w(output,:);
    d = pdist2(cw,temp1);

    ageSumBefore = ageSum;
    maxN = 1.5;

    gen    = ceil(Problem.FE/Problem.N);
    maxgen = ceil(Problem.maxFE/Problem.N);
    if gen <= round(0.9*maxgen)
        maxIter = 1;
        maxPZ = maxN;
        if size(w,1) == round(maxN*N)
            %      PlotResults(w, C)
            % plot(w(:,1),w(:,2),'k*')
            % hold on
            % plot(temp1(:,1),temp1(:,2),'rs')
            [~,rankFlag] = sort(flag,'descend');
            r = rankFlag(1:round(maxN*N)-N);
            %     [~,index] = max(flag,[],1);
            %     r = index;
            C(r, :) = [];
            C(:, r) = [];
            t(r, :) = [];
            t(:, r) = [];
            w(r, :) = [];
            E(r) = [];
            ageSumBefore(r,:) = [];
            flag(r,:) = [];
            flag = zeros(N,1);

        end
    else
        if size(w,1) < round(maxN*N) && isempty(genFlag)
            maxPZ = maxN;
            maxIter = 1;
        else
            maxPZ = 1;
            maxIter = 0;
        end
        if size(w,1) == round(maxN*N)
            maxPZ = 1;
            maxIter = 0;
            genFlag = gen;
        end
    end
    % plot(w(:,1),w(:,2),'k*')

    if isempty(genFlag)
        for iter = 1:maxIter
            for kk = 1:size(temp1,1)
                nx = nx + 1;
                x = temp1(kk,:);
                w = w./sum(w,2);
                d = pdist2(x, w);
                [~, SortOrder] = sort(d);
                s1 = SortOrder(1);
                s2 = SortOrder(2);

                % Aging: the age of all neighbours of s1 is increased by 1

                t(s1, :) = t(s1, :) + 1;
                t(:, s1) = t(:, s1) + 1;

                % Add Error
                E(s1) = E(s1) + d(s1)^2;

                % Adaptation
                w(s1,:) = w(s1,:) + epsilon_b*(x-w(s1,:));
                Ns1 = find(C(s1,:)==1);
                for j=Ns1
                    w(j,:) = w(j,:) + epsilon_n*(x-w(j,:));
                end

                % Create Link
                C(s1,s2) = 1;
                C(s2,s1) = 1;
                t(s1,s2) = 0;
                t(s2,s1) = 0;

                % Remove Old Links
                C(t>T) = 0;
                nNeighbor = sum(C);
                AloneNodes = (nNeighbor==0);
                if ~isempty(find(AloneNodes == true))
                    %         AloneNodes
                end
                C(AloneNodes, :) = [];
                C(:, AloneNodes) = [];
                t(AloneNodes, :) = [];
                t(:, AloneNodes) = [];
                w(AloneNodes, :) = [];
                E(AloneNodes) = [];
                ageSumBefore(AloneNodes,:) = [];
                flag(AloneNodes,:) = [];

                % Add New Nodes

                if mod(nx, L) == 0 && size(w,1) < round(maxPZ*N)
                    [~, q] = max(E);
                    [~, f] = max(C(:,q).*E);
                    r = size(w,1) + 1;
                    w(r,:) = (w(q,:) + w(f,:))/2;
                    C(q,f) = 0;
                    C(f,q) = 0;
                    C(q,r) = 1;
                    C(r,q) = 1;
                    C(r,f) = 1;
                    C(f,r) = 1;
                    t(r,:) = 0;
                    t(:, r) = 0;
                    E(q) = alpha*E(q);
                    E(f) = alpha*E(f);
                    E(r) = E(q);
                    ageSumBefore(r,:) = 0;
                    flag(r,:) = 0;
                end

                % Decrease Errors
                E = delta*E;
            end
        end

        w = w./sum(w,2);
        net.w = w;
        net.E = E;
        net.C = C;
        net.t = t;
        net.nx = nx;
        net.ageSumBefore = ageSumBefore;
        net.flag = flag;
    end

    if isempty(genFlag)
        output = [];
        for i = 1:size(w,1)
            neighbor = find(C(i,:)==1);
            [~,minInd] = min(w([i neighbor],:),[],1);
            [~,maxInd] = max(w([i neighbor],:),[],1);
            if ~isempty(find(minInd==1))
                output = [output i];
            elseif ~isempty(find(maxInd==1))
                output = [output i];
            end

        end

        for t = 1:size(output,2)
            s = output(:,t);
            x = mean(w(find(C(s,:)==1),:),1);
            w(s,:) = w(s,:) - 1*(x-w(s,:));
        end

        t = any(w<0,2);
        l = find(t == true);
        if ~isempty(l)
            for tt = 1:size(l,1)
                w(l(tt,1),any(w(l(tt,1),:)<0,1)) = 0;
            end
        end

        V = w;
        V = V.*repmat(scale,size(V,1),1);
    end

    if gen == genFlag
        N = size(wholeObj,1);
        wholeObj1 = (wholeObj - repmat(zmin,N,1));
        wholeObj = wholeObj1./scale;
        temp2 = wholeObj1./sum(wholeObj1,2);

        output = [];
        for i = 1:size(w,1)
            neighbor = find(C(i,:)==1);
            [~,minInd] = min(w([i neighbor],:),[],1);
            [~,maxInd] = max(w([i neighbor],:),[],1);
            if ~isempty(find(minInd==1))
                output = [output i];
            elseif ~isempty(find(maxInd==1))
                output = [output i];
            end

        end

        for t = 1:size(output,2)
            s = output(:,t);
            x = mean(w(find(C(s,:)==1),:),1);
            w(s,:) = w(s,:) - 1*(x-w(s,:));
        end

        t = any(w<0,2);
        l = find(t == true);
        if ~isempty(l)
            for tt = 1:size(l,1)
                w(l(tt,1),any(w(l(tt,1),:)<0,1)) = 0;
            end
        end

        V = [w.*repmat(scale,size(w,1),1);temp2];
    end
end


text:
function [scr_dev_pav,scr_eval_pav] = pav_calibrate_scores_dev_eval(scr_dev,scr_eval,key_dev,ndx_eval,prior,score_offset,addzeros)
% A function to calibrate a set of scores (eval) using another set
% as reference (dev).  The function requires a key for the dev
% scores (supervised) and an ndx for the eval scores
% (unsupervised).  It trains a PAV transformation and then applies
% the transformation to both sets of scores.  Comparing the act_dcf
% and min_dcf values at the operating point gives an indication of
% the success of the calibration (assuming the score distribution
% is the same in the two sets --- if not, don't bother calibrating).
% Inputs:
%   scr_dev: The development scores.
%   scr_eval: The eval scores.
%   key_dev: The Key for the development scores i.e. must have
%     'tar' and 'non' fields.
%   ndx_eval: The Ndx for the eval scores i.e. must have a
%     'trialmask' field.
%   prior: The effective target prior.
%   score_offset: This is used to make the transform monotonically
%     increasing by tilting each flat portion of the PAV output.
%     The tilt is controlled by this value, and making it zero
%     results in no tilt.
%   addzeros: If true, missing scores (required by key but not
%     present in scr) are added (with value zero).
% Outputs:
%   scr_dev_pav: The calibrated development scores.
%   scr_eval_pav: The calibrated eval scores.

assert(nargin==7)
assert(isa(scr_dev,'Scores'))
assert(isa(scr_eval,'Scores'))
assert(isa(key_dev,'Key'))
assert(isa(ndx_eval,'Ndx'))
assert(scr_dev.validate())
assert(scr_eval.validate())
assert(key_dev.validate())
assert(ndx_eval.validate())

[tar,non] = scr_dev.get_tar_non(key_dev);
logprint(Logger.Info,'training calibration\n')
pav_trans = pav_calibration(tar,non,score_offset);

logprint(Logger.Info,'calibrating scores\n')
scr_dev_pav = scr_dev.align_with_ndx(key_dev);
scr_eval_pav = scr_eval.align_with_ndx(ndx_eval);

scr_dev_pav = scr_dev_pav.transform(pav_trans);
scr_eval_pav = scr_eval_pav.transform(pav_trans);

logprint(Logger.Info,'displaying calibration results\n')
res = Results(scr_dev_pav,key_dev);
logprint(Logger.Info,'norm_act_dcf = %g, norm_min_dcf = %g, prbep = %g\n',res.get_norm_act_dcf(prior),res.get_norm_min_dcf(prior),res.get_prbep());

if addzeros
    scr_dev_pav = scr_dev_pav.set_missing_to_value(key_dev,0.0);
    scr_eval_pav = scr_eval_pav.set_missing_to_value(ndx_eval,0.0);
end

assert(scr_dev_pav.validate())
assert(scr_eval_pav.validate())

end

author: nesl


text:
function [f,A,b,sdcone,xtype] = optiReadSDPA(filename,dense,print)
%OPTIREADSDPA  Read an SDPA (.dat-s or .dat) file and convert to Matlab Values
%
% [f,A,b,sdcone] = optiReadSDPA(filename) reads the sparse file specified by
% filename and return the semidefinite problem in OPTI format. If you do
% not specify a file a extension, it will default to .dat-s.
%
% [f,A,b,sdcone] = optReadSDPA(filename,dense) allows the user to specify
% whether the file is dense or not. dense = 1 indicates a dense SDPA file,
% and the file extensions will default to .dat if not provided. The default
% is sparse.
%
% [f,A,b,sdcone] = optiReadSDPA(filename,dense,print) allows print out of
% reading progress.
%
% [f,A,b,sdcone,xtype] = optiReadSDPA(..) returns variable integrality.
%
% *sdcone will be a cell array, where each cell is a sparse matrix of the
% form [C(:) A0(:) A1(:) ... ], containing symmetric matrices.

%   Copyright (C) 2013 Jonathan Currie (IPL)
%#ok<*AGROW>

if(nargin < 2), dense = 0; end
if(nargin < 3), print = 0; end

% Make sure we have file extension!
if(isempty(strfind(filename,'.')))
    filename = [filename '.dat-s'];
end

%Open File
fid = fopen(filename);
if(fid == -1)
    error('Error opening file %s',filename);
end

%States
sNCON = 1;
sNBLOCKS = 2;
sSIZES = 3;
sRHS = 4;
sDNDATA = 5;
sSPDATA = 6;
sINTEGER = 7;

%Default Internals
eof = 0;
state = sNCON; cind = 0; bind = 1; rind = 1;
f = []; A = []; b = []; sdcone = []; xtype = [];

if(print)
    ind = strfind(filename,filesep);
    if(isempty(ind))
        file = filename;
    else
        file = filename(ind(end)+1:end);
    end
    if(dense)
        fprintf('Reading Dense SDPA File ''%s''...\n',file);
    else
        fprintf('Reading Sparse SDPA File ''%s''...\n',file);
    end
end
try
    %Main Loop
    while(~eof)
        %Read DAT File String
        str = fgetl(fid);
        %Check for empty
        if(~ischar(str) && str == -1)
            eof = 1;
            continue;
        end
        %Check for comment, skip if found
        if(state ~= sINTEGER && (str(1) == '*' || str(1) == '"'))
            if(~isempty(strfind(str,'*INTEGER')))
                state = sINTEGER;
            end
            continue;
        end
        %Check for equals (Sonja files)
        idx = strfind(str,'=');
        if(~isempty(idx))
            str = str(1:idx-1);
        end
        %Replace { } ( ) [ ] , ; : with spaces
        try
            str = regexprep(str,{'[',']',',','{','}','(',')'},' ');
        catch
            error('Error in SDPA File with "%s"',str);
        end
        %Decode based on state
        switch(state)
            case sNCON
                M = sscanf(str,'%d',1);
                xtype = repmat('C',1,M);
                %Check for errors
                if(M < 1 || M > 1e9), error('The number of constraints is negative or > 1e9'); end
                %Print what we found
                if(print), fprintf('Variables:              %6d\n',M); end
                state = sNBLOCKS;
            case sNBLOCKS
                N = sscanf(str,'%d',1);
                %Check for errors
                if(N < 1 || N > 1e9), error('The number of blocks is negative or > 1e9'); end
                %Print what we found
                if(print), fprintf('Number of Blocks:       %6d\n',N); end
                state = sSIZES;
            case sSIZES
                blocksize = sscanf(str,'%d');
                %Check length = number of blocks
                if(length(blocksize) ~= N), error('The number of blocks does not match the number of entries of block sizes'); end
                if(any(blocksize == 0) || any(blocksize > 1e9)), error('A block size is 0 or > 1e9'); end
                %Print Total Dims
                if(print)
                    fprintf('  Linear:               %6d\n',length(blocksize(blocksize < 0)));
                    fprintf('  Semidefinite:         %6d\n',length(blocksize(blocksize > 0)));
                    fprintf('Total Block Dimension:  %6d\n',sum(abs(blocksize)));
                end
                %Allocate problem data
                sdcone = cell(sum(blocksize > 0),1); j = 1;
                for i = 1:N
                    if(blocksize(i) > 0)
                        if(dense)
                            sdcone{j} = zeros(blocksize(i)^2,M+1);
                        else
                            sdcone{j} = sparse(blocksize(i)^2,M+1); %not the greatest...
                        end
                        j = j + 1;
                    end
                end
                %Create linear part
                ind = blocksize < 0;
                if(any(ind))
                    if(dense)
                        A = zeros(sum(abs(blocksize(ind))),M);
                    else
                        A = sparse(sum(abs(blocksize(ind))),M);
                    end
                    b = zeros(sum(abs(blocksize(ind))),1);
                end
                %Create index vector
                blockindex = zeros(N,1);  lp = 1; sdp = 1;
                for i = 1:N
                    if(blocksize(i) < 0)
                        blockindex(i) = lp; lp  = lp + 1;
                    else
                        blockindex(i) = sdp; sdp = sdp + 1;
                    end
                end
                state = sRHS;
            case sRHS
                f = sscanf(str,'%f'); f = -f;
                %Check length = number of constraints
                if(length(f) ~= M), error('The number of RHS elements does not match the number of constraints'); end
                if(dense)
                    state = sDNDATA;
                else
                    state = sSPDATA;
                end
            %Sparse Data File
            case sSPDATA
                data = sscanf(str,'%d %d %d %d %f');
                if(~isempty(data)) %could be {}
                    %Check for C matrix
                    if(data(1) == 0)
                        %Check block for linear
                        if(blocksize(data(2)) < 0)
                            if(data(3) ~= data(4)), error('A diagonal entry in C is not on the diagonal!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:data(2)-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)));
                            b(stIndex + data(3)) = -data(5);
                        else %semidef
                            sdcone{blockindex(data(2))}(data(3) + (data(4) - 1)*blocksize(data(2)),1) = data(5);
                            %if off-diagonal, place in transposed position (build symmetric matrices)
                            if(data(3) ~= data(4))
                                sdcone{blockindex(data(2))}(data(4) + (data(3) - 1)*blocksize(data(2)),1) = data(5);
                            end
                        end
                    else
                        %Check block for linear
                        if(blocksize(data(2)) < 0)
                            if(data(3) ~= data(4)), error('A diagonal entry in A is not on the diagonal!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:data(2)-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)));
                            A(stIndex + data(3),data(1)) = data(5);
                        else %semidef
                            sdcone{blockindex(data(2))}(data(3) + (data(4) - 1)*blocksize(data(2)),data(1)+1) = -data(5);
                            %if off-diagonal, place in transposed position (build symmetric matrices)
                            if(data(3) ~= data(4))
                                sdcone{blockindex(data(2))}(data(4) + (data(3) - 1)*blocksize(data(2)),data(1)+1) = -data(5);
                            end
                        end
                    end
                end
            %Dense Data File
            case sDNDATA
                data = sscanf(str,'%f');
                if(~isempty(data)) %could be {}
                    %Check for C matrix
                    if(cind == 0)
                        %Check block for linear
                        if(blocksize(bind) < 0)
                            if(length(data) ~= abs(blocksize(bind))), error('Data read in dense entry C (diagonal) is not the correct length!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:bind-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)))+1;
                            b(stIndex:stIndex+abs(blocksize(bind))-1) = -data;
                        else
                            if(length(data) == blocksize(bind))
                                sdcone{blockindex(bind)}(rind:blocksize(bind):blocksize(bind)^2) = data;
                                rind = rind + 1;
                            elseif(length(data) == blocksize(bind)^2)
                                %Create matrix and transpose it
                                temp = reshape(data,blocksize(bind),blocksize(bind))';
                                sdcone{blockindex(bind)}(:,1) = temp(:);
                                rind = rind + blocksize(bind);
                            else
                                error('Data read in dense entry C matrix is not the correct length. Expected first row of a matrix, or the entire matrix (all rows and columns, row major).');
                            end
                        end
                    else
                        %Check block for linear
                        if(blocksize(bind) < 0)
                            if(length(data) ~= abs(blocksize(bind))), error('Data read in dense entry A%d (diagonal) is not the correct length!',cind); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:bind-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)))+1;
                            A(stIndex:stIndex+abs(blocksize(bind))-1,cind) = data;
                        else
                            if(length(data) == blocksize(bind))
                                sdcone{blockindex(bind)}(rind + (cind)*blocksize(bind)^2:blocksize(bind):blocksize(bind)^2 + (cind)*blocksize(bind)^2) = -data;
                                rind = rind + 1;
                            elseif(length(data) == blocksize(bind)^2)
                                %Create matrix and transpose it
                                temp = reshape(data,blocksize(bind),blocksize(bind))';
                                sdcone{blockindex(bind)}(:,cind+1) = -temp(:);
                                rind = rind + blocksize(bind);
                            else
                                error('Data read in dense entry A%d matrix is not the correct length. Expected first row of a matrix, or the entire matrix (all rows and columns, row major).',cind);
                            end
                        end
                    end
                    %Increment pointers
                    if(blocksize(bind) < 0 || rind > blocksize(bind))
                        rind = 1;
                        bind = bind + 1;
                    end
                    if(bind > N)
                        bind = 1;
                        cind = cind + 1;
                    end
                end
            case sINTEGER
                if(isempty(xtype))
                    error('Your integer declaration must come last!');
                end
                idx = sscanf(str,'*%d');
                xtype(idx) = 'I';

        end
        eof = feof(fid);
    end
catch ME
    fclose(fid);
    rethrow(ME);
end

%Ensure sparse return args
if(~isempty(sdcone))
    for i = 1:length(sdcone)
        if(~issparse(sdcone{i}))
            sdcone{i} = sparse(sdcone{i});
        end
    end
end
if(~isempty(A))
    if(~issparse(A)), A = sparse(A); end
end
if(print)
    nz = 0;
    if(~isempty(sdcone))
        for i = 1:length(sdcone)
            nz = nz + nnz(sdcone{i});
        end
    end
    fprintf('Total Semidefinite NZs: %6d\n',nz);
    fprintf('Finished Reading SDPA File\n\n');
end
%Close
fclose(fid);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot