text:
function varargout = process_matlab_eval2( varargin )
% PROCESS_MATLAB_EVAL2: Evaluate a matlab script.

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Run Matlab command';
    sProcess.FileTag     = 'matlab';
    sProcess.Category    = 'Filter2';
    sProcess.SubGroup    = 'Other';
    sProcess.Index       = 902;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/TutUserProcess#Alternative';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'results', 'timefreq', 'raw', 'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isPaired    = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = -1;
    sProcess.processDim  = [];
    % === Matlab command
    sProcess.options.matlab.Comment = 'Type your Matlab script below: ';
    sProcess.options.matlab.Type    = 'textarea';
    sProcess.options.matlab.Value   = ['% Input variables: DataA, DataB, TimeVector' 10  ...
                                       '% Output variables: Data, Comment, Condition', 10 ...
                                       'Data = DataA - DataB;' 10 ...
                                       'Comment = ''New file'';' 10 ...
                                       'Condition = ''NewCondition'';'];
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data', 'raw'};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function sOutput = Run(sProcess, sInputA, sInputB) %#ok<DEFNU>
    % Prepare input
    DataA = sInputA.A;
    DataB = sInputB.A;
    TimeVector = sInputA.TimeVector;
    Comment    = sInputA.Comment;
    Condition  = [];

    % Evaluate Matlab code
    try
        eval(sProcess.options.matlab.Value);
    catch
        % Get error message
        e = lasterror();
        if ~isempty(e)
            % Remove HTML tags
            errMsg = str_striptag(e.message);
            % Strip "Error:"
            errMsg = strrep(errMsg, 'Error: ', '');
            % Strip first line
            iNewLine = find(errMsg == 10);
            if ~isempty(iNewLine)
                errMsg = errMsg(iNewLine(1)+1:end);
            end
            errMsg = ['Error while executing Matlab command:' 10 errMsg];
        else
            errMsg = 'Error while executing Matlab command.';
        end
        % Report error
        bst_report('Error', sProcess, sInputA, errMsg);
        sOutput = [];
        return;
    end

    % Check if time was modified with only a few sensors selected
    if isfield(sProcess.options, 'sensortypes') && isfield(sProcess.options.sensortypes, 'Value') && ~isempty(sProcess.options.sensortypes.Value) && (length(TimeVector) ~= length(sInputA.TimeVector))
        sInputA = [];
        bst_report('Error', sProcess, sInputA, ...
            ['This process can modify the time definition only if all the channels are processed at once.' 10 ...
             'Either you leave the option "Sensor types or names" empty to process all the channels at once,' 10 ...
             'or your keep the number of time points of the original file.']);
        return;
    end

    % Output structure
    sOutput = sInputA;
    % Output condition name
    sOutput.A         = Data;
    sOutput.Condition = Condition;
    sOutput.Comment   = Comment;
end

author: brainstorm-tools


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
    = cornk_run(fid, data, K, L, pc, tc, opts)
% This program simulates the CORN-K algorithm. This definition is slightly
% different from the algorithm in LHG11, as we define a percentage of all
% experts here, rather than a specified number of experts. NOTE.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%    = cornk_run(fid, data, K, L, pc, tc, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
% exp_ret: experts' return
%
% data: market sequence vectors
% fid: handle for write log file
% K: maximum window size
% L: splits into L parts, in each K, useless in CORN-U, L=1
% pc: top-k percentage of K*L
% tc: transaction cost rate parameter
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%            = cornk_run(fid, data, K, L, pc, tc, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m] = size(data);

% Variables for return, start with uniform weight
% cumprod_ret = 1;
% daily_ret = 1;
% weight = ones(nStocks, 1)/nStocks;

cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% Variables for expert
exp_ret = ones(K, L);
exp_w = ones(K*L, m)/m;

% print file head
fprintf(fid, '-------------------------------------\n');
if (~opts.quiet_mode)
    fprintf(fid, 'Parameters [K=%d, L=%d, pc=%f, tc=%f\n]', K, L, pc, tc);
    fprintf(fid, 'day\t Daily Return\t Total return\n');
end
if (opts.progress)
	progress = waitbar(0,'Executing Algorithm...');
end
for t = 1:1:n,
    % Calculate t's portfolio
    if (t >=2)
        [day_weight, exp_w] ...
            = cornk_kernel(data(1:t-1, :), K, L, pc, exp_ret, exp_w);
    end

    % Normalize the constraint
    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    % Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % Normalize the portfolio
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Cal t's experts return
    for k=1:K,
        for l=1:L,
            exp_ret(k, l) = exp_ret(k, l)*data(t, :)*exp_w((k-1)*L+l, :)';
        end
    end

    % Debug information
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cum_ret);
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cum_ret);
        end
    end
    if (opts.progress)
		if mod(t, 50) == 0
			waitbar((t/n));
		end
	end
end

% Debug Information
fprintf(fid, 'CORN topK (K:%d, L:%d, pc:%.2f, tc:%.4f), Final return: %.2f\n', ...
    K, L, pc, tc, cum_ret);

fprintf(fid, 'CORN topK, Experts return:\n');
fprintf(fid, '%f\t', exp_ret);
fprintf(fid, '\n');
fprintf(fid, '-------------------------------------\n');

fprintf(1, 'CORN topK (K:%d, L:%d, pc:%.2f, tc:%.4f), Final return: %.2f\n', ...
    K, L, pc, tc, cum_ret);
if (~opts.quiet_mode)
    fprintf(1, 'CORN topK, Experts return:\n');
    fprintf(1, '%f\t', exp_ret);
    fprintf(1, '\n');
end
fprintf(fid, '-------------------------------------\n');
	if (opts.progress)
		close(progress);
	end
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function visData = sc_vis_planar_structure(modelPlane, modelReg)
%
%     sc_vis_planar_structure:
%     visualize the extracted planar structure

%%
visData = [];

end

repo: SelfExSR


text:
tic
nt0 = ops.nt0;

if ~isempty(ops.chanMap)
    load(ops.chanMap);
    chanMapConn = chanMap(connected>1e-6);
else
    chanMapConn = 1:ops.Nchan;
end
batch_path = fullfile(root, 'batches');
if ~exist(batch_path, 'dir')
    mkdir(batch_path);
end
NchanTOT = ops.NchanTOT;

d = dir(fullfile(root, fname));
ops.sampsToRead = floor(d.bytes/NchanTOT/2);

NT          = 128*1024+ ops.ntbuff;
NTbuff      = NT + 4*ops.ntbuff;
Nbatch      = ceil(d.bytes/2/NchanTOT /(NT-ops.ntbuff));

% load data into patches, filter, compute covariance, write back to
% disk

fprintf('Time %3.0fs. Loading raw data... \n', toc);
fid = fopen(fullfile(root, fname), 'r');
ibatch = 0;
Nchan = ops.Nchan;

Nchans = ops.Nchan;
ts = [1:1:nt0]';

clear stimes
% for iNN = 1:size(rez.W,2)
%     stimes{iNN} = rez.st3pos(rez.st3pos(:,2)==iNN,1);
% end
stimes = gtimes;

Wraw = zeros(nt0, Nchans, numel(stimes));

while 1
    ibatch = ibatch + 1;

    offset = max(0, 2*NchanTOT*((NT - ops.ntbuff) * (ibatch-1) - 2*ops.ntbuff));
    if ibatch==1
        ioffset = 0;
    else
        ioffset = ops.ntbuff;
    end
    fseek(fid, offset, 'bof');
    buff = fread(fid, [NchanTOT NTbuff], '*int16');

    if isempty(buff)
        break;
    end
    nsampcurr = size(buff,2);
    if nsampcurr<NTbuff
        buff(:, nsampcurr+1:NTbuff) = repmat(buff(:,nsampcurr), 1, NTbuff-nsampcurr);
    end

    offset = (NT-ops.ntbuff)*(ibatch-1)-64 - 40;

    buff = buff(chanMapConn,:)';
    %
    for iNN = 1:numel(stimes)
        st = stimes{iNN} - offset;
        st(st<0) = [];
        st(st>NT-ops.ntbuff) = [];

        if ~isempty(st)
            inds = repmat(st', nt0, 1) + repmat(ts, 1, numel(st));

            Wraw(:,:,iNN) = Wraw(:,:,iNN) + ...
                squeeze(sum(reshape(buff(inds, :), nt0, numel(st), Nchans),2));
        end
    end

end
fclose(fid);

for iNN = 1:numel(stimes)
     Wraw(:,:,iNN) = Wraw(:,:,iNN)/numel(stimes{iNN});
end
fprintf('Time %3.2f. Mean waveforms computed... \n', toc);

author: cortex-lab

repo: KiloSort

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


text:
function [rd, boardName] = hdlcoder_ref_design_customization
% Reference design plugin registration file
% 1. The registration file with this name inside of a board plugin folder
%    will be picked up
% 2. Any registration file with this name on MATLAB path will also be picked up
% 3. The registration file returns a cell array pointing to the location of
%    the reference design plugins
% 4. The registration file also returns its associated board name
% 5. Reference design plugin must be a package folder accessible from
%    MATLAB path, and contains a reference design definition file

%   Copyright 2013-2014 The MathWorks, Inc.

rd = {...
    'AnalogDevices.fmcomms2.zed.plugin_rd_rx', ...
    'AnalogDevices.fmcomms2.zed.plugin_rd_tx', ...
    'AnalogDevices.fmcomms2.zed.plugin_rd_rxtx', ...
    };

boardName = 'AnalogDevices FMCOMMS2/3 ZED';

end

author: analogdevicesinc


text:
% PutSegyHeader : Writes SEG-Y header to disk.
% PutSegyHeader(segyid,SegyHeader)
%
  % (C) 2001-2016, Thomas Mejer Hansen, thomas.mejer.hansen@gmail.com
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
function SegyHeader=PutSegyHeader(segyid,SegyHeader);

NA=0;

if isfield(SegyHeader,'Rev')==0
  SegyHeader.Rev=GetSegyHeaderBasics;
end

if (isfield(SegyHeader,'DataSampleFormat')==0)&(isfield(SegyHeader,'SegyFormatRevisionNumber')==0)
  SegyHeader.DataSampleFormat=5; % '5'->4-byte IEEE floating point
  SegyHeader.SegyFormatRevisionNumber=256;
  SegymatVerbose([mfilename,' : Using datasample format : ',SegyHeader.Rev(2).DataSampleFormat(SegyHeader.DataSampleFormat).name],1)
  PrintSegyFormatRevsisionNumber(SegyHeader,1);
end

% SET Revsion number. Always use revision 1, if not otherwise specified.
if (isfield(SegyHeader,'SegyFormatRevisionNumber')==0)
  SegyHeader.SegyFormatRevisionNumber=256;
  PrintSegyFormatRevsisionNumber(SegyHeader,1);
end

% SET DATA SAMPLE FORMAT IF NOT SPECIFIED
% Rev0->IBM FLOATING POINT
% Rev1->IEEE FLOATING POINT
if (isfield(SegyHeader,'DataSampleFormat')==0)
  if SegyHeader.SegyFormatRevisionNumber==0
    SegyHeader.DataSampleFormat=1; % '1'->IBM floating point
    SegymatVerbose([mfilename,' : Using datasample format : ',SegyHeader.Rev(1).DataSampleFormat(SegyHeader.DataSampleFormat).name],1)
  else
    SegyHeader.DataSampleFormat=5; % '5'->4-byte IEEE floating point
    SegymatVerbose([mfilename,' : Using datasample format : ',SegyHeader.Rev(2).DataSampleFormat(SegyHeader.DataSampleFormat).name],1)
  end
end

if isfield(SegyHeader,'TextualFileHeader')
  % MAKE SURE LENGTH IF Text.. is 3200
  if length(SegyHeader.TextualFileHeader)<3200
    dummyTXT=SegyHeader.TextualFileHeader;
    SegyHeader.TextualFileHeader=zeros(3200,0);
    SegyHeader.TextualFileHeader(1:length(dummyTXT))=dummyTXT;
  end
else
  SegyHeader.TextualFileHeader=sprintf('%3200s','SEGY READER (tmh@gfy.ku.dk)');
end

if ~isfield(SegyHeader,'Job'), SegyHeader.Job=NA; end
if ~isfield(SegyHeader,'Line'), SegyHeader.Line=NA; end
if ~isfield(SegyHeader,'Reel'), SegyHeader.Reel=NA; end
if ~isfield(SegyHeader,'DataTracePerEnsemble'), SegyHeader.DataTracePerEnsemble=NA;end
if ~isfield(SegyHeader,'AuxiliaryTracePerEnsemble'), SegyHeader.AuxiliaryTracePerEnsemble=0;end
if ~isfield(SegyHeader,'dt'), SegyHeader.dt=4000; end
if ~isfield(SegyHeader,'dtOrig'), SegyHeader.dtOrig=NA; end
if ~isfield(SegyHeader,'ns'), SegyHeader.ns=NA; end
if ~isfield(SegyHeader,'nsOrig'), SegyHeader.nsOrig=NA; end
if ~isfield(SegyHeader,'DataSampleFormat'), SegyHeader.DataSampleFormat=5; end % '5'->4-byte IEEE floating point
if ~isfield(SegyHeader,'EnsembleFold'), SegyHeader.EnsembleFold=NA; end
if ~isfield(SegyHeader,'TraceSorting'), SegyHeader.TraceSorting=NA; end
if ~isfield(SegyHeader,'VerticalSumCode'), SegyHeader.VerticalSumCode=NA; end
if ~isfield(SegyHeader,'SweepFrequencyStart');SegyHeader.SweepFrequencyStart=NA; end
if ~isfield(SegyHeader,'SweepFrequencyEnd');SegyHeader.SweepFrequencyEnd=NA; end
if ~isfield(SegyHeader,'SweepLength');SegyHeader.SweepLength=NA; end
if ~isfield(SegyHeader,'SweepType');SegyHeader.SweepType=NA; end
if ~isfield(SegyHeader,'SweepChannel');SegyHeader.SweepChannel=NA; end
if ~isfield(SegyHeader,'SweepTaperlengthStart');SegyHeader.SweepTaperlengthStart=NA; end
if ~isfield(SegyHeader,'SweepTaperLengthEnd');SegyHeader.SweepTaperLengthEnd=NA; end
if ~isfield(SegyHeader,'TaperType');SegyHeader.TaperType=NA; end
if ~isfield(SegyHeader,'CorrelatedDataTraces');SegyHeader.CorrelatedDataTraces=NA; end
if ~isfield(SegyHeader,'BinaryGain');SegyHeader.BinaryGain=NA; end
if ~isfield(SegyHeader,'AmplitudeRecoveryMethod');SegyHeader.AmplitudeRecoveryMethod=NA; end
if ~isfield(SegyHeader,'MeasurementSystem');SegyHeader.MeasurementSystem=1; end %1-Meters, 2-Feet
if ~isfield(SegyHeader,'ImpulseSignalPolarity');SegyHeader.ImpulseSignalPolarity=NA; end
if ~isfield(SegyHeader,'VibratoryPolarityCode');SegyHeader.VibratoryPolarityCode=NA; end
% 3261-3500 UNASSIGNED
if ~isfield(SegyHeader,'Unassigned1');SegyHeader.Unassigned1=NA*ones(1,120); end

if ~isfield(SegyHeader,'FixedLengthTraceFlag');SegyHeader.FixedLengthTraceFlag=1; end % ALL TRACES HAVE THE SAME dt AND ns
if ~isfield(SegyHeader,'NumberOfExtTextualHeaders'); SegyHeader.NumberOfExtTextualHeaders=0; end % WE DO NOT YET USE TEXTURAL HEADERS as can be done according to draft 6
% 3506-3600 UNASSIGNED
if ~isfield(SegyHeader,'Unassigned2');SegyHeader.Unassigned2=NA*ones(1,47); end

% IF EXTENDED TEXTUAL FILE HEADERS EXISTS THEN WRITE THEM TO DISK

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WRITE TO DISK
%

fseek(segyid,0,'bof');

fwrite(segyid,SegyHeader.TextualFileHeader(1:3200),'uchar');         % 1-3200
fwrite(segyid,SegyHeader.Job,'int32');                           % 3204

fwrite(segyid,SegyHeader.Line,'int32');                          % 3208
fwrite(segyid,SegyHeader.Reel,'int32');                          % 3212
fwrite(segyid,SegyHeader.DataTracePerEnsemble,'int16');        % 3214
fwrite(segyid,SegyHeader.AuxiliaryTracePerEnsemble,'int16');   % 3216
fwrite(segyid,SegyHeader.dt,'uint16');                          % 3218
fwrite(segyid,SegyHeader.dtOrig,'uint16');                      % 3220
fwrite(segyid,SegyHeader.ns,'uint16');                          % 3222
fwrite(segyid,SegyHeader.nsOrig,'uint16');                      % 3224
fwrite(segyid,SegyHeader.DataSampleFormat,'int16');            % 3226
fwrite(segyid,SegyHeader.EnsembleFold,'int16');
fwrite(segyid,SegyHeader.TraceSorting,'int16');               % 3228
fwrite(segyid,SegyHeader.VerticalSumCode,'int16');            % 3230

fwrite(segyid,SegyHeader.SweepFrequencyStart,'int16');        % 3232
fwrite(segyid,SegyHeader.SweepFrequencyEnd,'int16');          % 3234
fwrite(segyid,SegyHeader.SweepLength,'int16');                % 3236
fwrite(segyid,SegyHeader.SweepType,'int16');                  % 3238
fwrite(segyid,SegyHeader.SweepChannel,'int16');               % 3240
fwrite(segyid,SegyHeader.SweepTaperlengthStart,'int16');               % 3242
fwrite(segyid,SegyHeader.SweepTaperLengthEnd,'int16');               % 3244
fwrite(segyid,SegyHeader.TaperType,'int16');               % 3246

fwrite(segyid,SegyHeader.CorrelatedDataTraces,'int16');               % 3248

fwrite(segyid,SegyHeader.BinaryGain,'int16');               % 3250
fwrite(segyid,SegyHeader.AmplitudeRecoveryMethod,'int16');               % 3252

fwrite(segyid,SegyHeader.MeasurementSystem,'int16');               % 3254

fwrite(segyid,SegyHeader.ImpulseSignalPolarity,'int16');               % 3256
fwrite(segyid,SegyHeader.VibratoryPolarityCode,'int16');               % 3258

% 3261-3500 UNASSIGNED1 => (120 int32 = 240 bytes)
fwrite(segyid,SegyHeader.Unassigned1,'int16');               % 3260
%SegymatVerbose(ftell(segyid),10);
%fseek(segyid,3500,'bof');
fwrite(segyid,SegyHeader.SegyFormatRevisionNumber,'uint16');   % 3500
fwrite(segyid,SegyHeader.FixedLengthTraceFlag,'int16');        % 3502
fwrite(segyid,SegyHeader.NumberOfExtTextualHeaders,'uint16');        % 3504

% 3506-3600 UNASSIGNED2 => 94/2=47 int16
fwrite(segyid,SegyHeader.Unassigned2,'int16');               % 3506

%
%
%
if SegyHeader.NumberOfExtTextualHeaders>0
  n=SegyHeader.NumberOfExtTextualHeaders;
  SegymatVerbose(sprintf('Writing %d Extended Textual File Headers',n));
  if isfield(SegyHeader,'ExtTextualHeader')
    % MAKE SURE LENGTH IF Text.. is 3200*n
    if length(SegyHeader.ExtTextualHeader)<(3200*n)
      dummyTXT=SegyHeader.ExtTextualHeader;
      SegyHeader.ExtTextualHeader=zeros(n*3200,0);
      SegyHeader.ExtTextualHeader(1:length(dummyTXT))=dummyTXT;
    end
  else
    SegyHeader.ExtTextualHeader=sprintf('%3200s','SEGY READER (tmh@gfy.ku.dk)');
    for i=2:n,
      SegyHeader.ExtTextualHeader=[SegyHeader.ExtTextualHeader,sprintf('%3200s','SEGY READER (tmh@gfy.ku.dk)')];
    end
  end
  % WRITE TEXTUAL FILE HEADER
  fwrite(segyid,SegyHeader.ExtTextualHeader(1:(n*3200)),'uchar');         % 1:(n*3200)
end

Revision=SegyHeader.SegyFormatRevisionNumber;
if Revision>0, Revision=1; end
Format=SegyHeader.Rev(Revision+1).DataSampleFormat(SegyHeader.DataSampleFormat).name;
SegymatVerbose([mfilename,' - SegyRevision ',sprintf('%0.4g',Revision),', ',Format],2)

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function R = SeparateAngles(X,deep,is_real)
% SeparateAngles: Frequency angular partitioning
%  Usage:
%     R = SeparateAngles(X,deep)
%  Inputs:
%    X     m by m matrix (m = 2*2^(j+1)).
%          The input is the Fourier transform of an object obtained
%          after windowing with a bandpass Meyer window which
%          isolates frequencies in the range 2^(j-1) <= |k| <=
%          2^j. We interpret the array as samples at the frequency
%          points (k1,k2), -2^(j+1) <= k1,k2 < 2^(j+1)
%
%   deep   2^deep is the number of angular wedges in each of the
%          directional  panels NESW (North, East, South, West).
%  Outputs:
%    R      Array of smoothly localized Fourier samples
%           size = 4 * m * 2^j *  2*2^(j+1)/m, m = 2^deep
%
%     (i)   first index -- cardinal point 1 = W, 2 = E, 3 = N, 4 = S.
%     (ii)  second index -- indexes orientations according to the
%              organization below; e.g. (1,3) would correspond to
%              the wedge marked with a star.
%
%                                j
%             ------------------------------------->
%            |
%            |  1                                4
%            |
%            |
%            |
%            |  2                                3
%            |
%       i    |
%            |
%            |  3*                                2
%            |
%            |
%            |
%            |  4                                1
%            \/
%
%
%                     West                      East
%
%              North and South are obtained by transposition.
%
%     (iii) third and fourth indices-- frequency points, regular
%           samples along parallel slanted lines (shape of a parallelogram)
%
% Description
%   Smoothly localizes the FT near angular wedges
%   Operates by interpolating the FT along vertical and horizontal
%   directions
% See Also
%   Adj_SeparateAngles, AtA, AtA_Toeplitz
%
% By Emmanuel Candes, 2003-2004

  if nargin < 3
    is_real = 0;
  end

  n  = size(X,1);
  n2 = n/2;
  boxlen = n/2^deep;
  boxcnt = 2^deep;

  R = zeros(4,boxcnt,n2/2,2*boxlen);

  [ix,w] = DetailMeyerWindow([n2/4 n2/2],3);
  lx = reverse(-ix);

  m      =  0:(boxcnt - 1);
  ym     =  (m-boxcnt/2).*boxlen + boxlen/2;
  slope  =  -ym./n2;

  % recall that fft_mid0 operates along columns
  % Columns: West and East

  F = ifft_mid0(X).*sqrt(n);

  for r = 1:(n2/2),
    k = lx(r); t = n2 + k + 1;
    shift  =  ym + slope.*(k+n2);
    alpha =  -k./n2;
    w = MakeSineWindow(boxlen,alpha);
    y = Evaluate_FT(F(:,t),shift,boxlen,0,w);
    R(1,:,r,:)  =  reshape(y,2*boxlen,boxcnt).';
    rsym = n2/2+1-r;
    if (is_real)
      R(2,:,rsym,:) = conj(fliplr(squeeze(R(1,:,r,:))));
    else
      t = n2 - k + 1;
      shift = -shift + 1;
      y = Evaluate_FT(F(:,t),shift,boxlen,0,w);
      R(2,:,rsym,:) = reshape(y,2*boxlen,boxcnt).';
    end
  end

  % Rows: North and South
  F = ifft_mid0(X.').*sqrt(n);
  for r = 1:(n2/2),
    k = lx(r);  t = n2 + k + 1;
    shift  =  ym + slope.*(k+n2);
    alpha =  -k./n2;
    w = MakeSineWindow(boxlen,alpha);
    y = Evaluate_FT(F(:,t),shift,boxlen,0,w);
    R(3,:,r,:) = reshape(y,2*boxlen,boxcnt).';

    rsym = n2/2+1-r;
    if (is_real)
      R(4,:,rsym,:) = conj(fliplr(squeeze(R(3,:,r,:))));
    else
      t = n2 - k + 1;
      shift = -shift + 1;
      y = Evaluate_FT(F(:,t),shift,boxlen,0,w);
      R(4,:,rsym,:) = reshape(y,2*boxlen,boxcnt).';
    end
  end

author: thomaskuestner


text:
clear; clc; close all;
dataset = 'grid';
img = imread(strcat('../noncvx-lowrank(ICDM 2015)/data/images/', dataset, '.jpg'));

if(size(img, 2) > size(img, 1))
    img = permute(img, [2,1,3]);
end

imgSize = size(img);

img = double(img)/255;
imgMean = mean(img(:));
img = img - imgMean;
img = img/std(img(:));

img = reshape(img, imgSize(1), prod(imgSize)/imgSize(1));

missRatio = 0.95;
noisRatio = 0.05;

mask = (rand(size(img)) > missRatio);
G = randn(size(img))*noisRatio;
traD = img + G;
traD = traD.*mask;
traD = sparse(traD);

clear maxk missRatio noisRatio mask G;

para.maxIter = 5000;
para.tol = 1e-6;
para.maxR = 50 ;

%% ---------------------------------------------------------------
mask = (rand(size(img)) > 0.9);
tstD = img.*mask;
tstD = sparse(tstD);

[tstRow, tstCol, tstVal] = find(tstD);

para.test.row  = tstRow;
para.test.col  = tstCol;
para.test.data = tstVal;
para.test.m = size(tstD, 1);
para.test.n = size(tstD, 2);

clear tstRow tstCol tstVal tstD;

%% ---------------------------------------------------------------
lambdaMax = 4;
gridLambda = lambdaMax*(0.8).^(0:9);

gridNMSE = zeros(2, size(gridLambda,2 ));
gridRank = zeros(1, size(gridLambda,2 ));
for g = 1:length(gridLambda)
    lambda = gridLambda(g);

    % [U, S, V] = SoftImpute(traD, lambda, para );
    [U, S, V] = AISImpute(traD, lambda, para );
    gridRank(g) = nnz(S);

    X = U*S*V';
    gridNMSE(1, g) = sqrt(norm(X - img, 'fro')^2/numel(img));

    [ U, S, V ] = PostProcess(traD, U, V, S);
    X = U*S*V';

    gridNMSE(2, g) = sqrt(norm(X - img, 'fro')^2/numel(img));

    if(g > 1 && gridNMSE(2, g) > gridNMSE(2, g - 1))
        break;
    end
end

gridNMSE = gridNMSE(2,1:g);
[~, lambda] = min(gridNMSE);
gndRank = gridRank(lambda);
para.maxR = ceil(gndRank*1.2);

lambda = gridLambda(lambda);

clear gridNMSE gridRank g X U S V gridLambda;

%% active --------------------------------------------------------
method = 1;
t = tic;
[U, S, V, out{method}] = ActiveSubspace(traD, lambda, para );
% [U, S, V, out{1}]=mc_alt(traData', lambda, para);
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
hold on;
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
hold on;
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('active');

%% boost ---------------------------------------------------------
method = 2;
t = tic;
[U, S, V, out{method}] = Boost( traD, lambda, para);
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
hold on;
plot(out{method}.Time, out{method}.RMSE);
figure(2);
hold on;
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('boost');

%% TR ------------------------------------------------------------
method = 3;
t = tic;
[U, S, V, out{method}] = MMBS( traD, lambda, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('TR');

%% ALT-Impute ----------------------------------------------------
method = 4;
t = tic;
[U, S, V, out{method}] = SoftImputeALS( traD, lambda, para.maxR, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('ALT-Impute');

%% SSGD ----------------------------------------------------------
method = 5;
t = tic;
[U, S, V, out{method}] = SSGD( traD, lambda, gndRank, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('SSGD');

%% LMaFit --------------------------------------------------------

gridNMSE = zeros(1, 10);
for g = 1:10
    [U, S, V] = FixedRank( traD, g, para );

    X = U*S*V';
    gridNMSE(g) = sqrt(norm(X - img, 'fro')^2/numel(img));

    if(g > 1 && gridNMSE(g) > gridNMSE(g - 1))
        break;
    end
end

gridNMSE = gridNMSE(1:g);
[~, rnk] = min(gridNMSE);

clear gridNMSE X g U S V;

method = 6;
t = tic;
[U, S, V, out{method}] = FixedRank( traD, rnk, para );
Time(method) = toc(t);

X = U*S*V';
RMSE(method) = out{method}.RMSE(end);
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('LMaFit');

%% APG -----------------------------------------------------------
method = 7;
t = tic;
[U, S, V, out{method}] = APGMatComp( traD, lambda, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('APG');

%% Soft-Impute ---------------------------------------------------
method = 8;
t = tic;
[U, S, V, out{method}] = SoftImpute( traD, lambda, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('Soft-Impute');

%% AIS-Impute ----------------------------------------------------
method = 9;
t = tic;
[U, S, V, out{method}] = AISImpute( traD, lambda, para );
Time(method) = toc(t);

[ U, S, V ] = PostProcess(traD, U, V, S);
X = U*S*V';
RMSE(method) = sqrt(norm(X - img, 'fro')^2/numel(img));
clear U S V t;

figure(1);
plot(out{method}.Time, out{method}.RMSE);
figure(2);
semilogy(out{method}.Time, out{method}.obj - min(out{method}.obj));
figure;
X = X + imgMean;
X = reshape(X, imgSize);
imshow(X, []);
title('AIS-Impute');

% clear gndRank img imgMean lambda lambdaMax method rnk traD X mask;

author: HKUST-KnowComp

repo: FMG

save_path: github-repos/MATLAB/HKUST-KnowComp-FMG


text:
function varargout = create_main_inputfile(GID)
fid = fopen('main_inputfile.m','w+');  % open file

fprintf(fid, ...
 ['function varargout = main_inputfile(varargin)\n' ...
  '%% Task(menu).name   - name of task. Task(1).name is ''File'', Task(2).name is\n' ...
  '%%                     ''Model'', Task(8).name is ''Help''.\n' ...
  '%% Task(menu).string - For each menu item (or any sub-items), there are 4\n' ...
  '%%                     parameters that needed to be specified.\n' ...
  '%%                     (1) index -- There are up to 3 cascading levels of\n' ...
  '%%                     menu items. Each level can have up to 9 menu items.\n' ...
  '%%                     For example, "100" represents menu item 1 of level 1.\n' ...
  '%%                     "236" is item 2 selected from level 1 and item 3 from\n' ...
  '%%                     level 2 and item 6 from level 3.\n' ...
  '%%                     (2) label -- menu item label corresponding to index\n' ...
  '%%                     (3) type  -- menu item type; type ''doc1'' if the\n' ...
  '%%                     action is to open html or pdf file; ''run1'' to execute\n' ...
  '%%                     a matlab function; etc.\n' ...
  '%%                     (4) operands -- for type ''doc1'', this would be the\n' ...
  '%%                     html or pdf file; for type ''run1'', this would be the\n' ...
  '%%                     name of the function and its arguments.\n' ...
  '%%                     You can enter the menu items in any order\n' ...
  ]);
fprintf(fid, ...
    ['clear all\n' ...
     'if nargin > 0\n' ...
     '%% update gui4gui_inputfile with modified "Task" from the gui4gui input GUI\n' ...
     '    GID = varargin{1};\n' ...
     '    %% insert any additional modifications here\n' ...
     '    save gui4gui_inputfile GID;\n' ...
     'else\n' ...
     '%% the gui4gui package comes with this mfile pre-loaded with the setup\n' ...
     '%% for the MODE code. Some of these data are universal for all codes if\n' ...
     '%% you prepare PDF/HTML files and label them with the same names, such as\n' ...
     '%% to call the tutorial file Tutorial.ppt or name the abstract Abstract.pdf.\n' ...
     '%% You may also edit this file and change any contents that you wish.\n' ...
     '%% However, you should not remove any fields in Task. If you don''t need\n' ...
     '%% a specific menu item, just remove its content, like this\n' ...
     '%% Task(1).label = {};\n' ...
     '%%\n' ...
     '%% Define number of menus (File, Model, Articles, Tutorial, Examples, Run, Code, Help)\n' ...
     '%%\n']);
fprintf(fid,'Nmenus = %d;\n',GID.Nmenus);
fprintf(fid,'\n');
Nmenus = GID.Nmenus;
for menu=1:Nmenus
fprintf(fid, ...
    'Task(%d).name = ''%s'';\n', menu,GID.task(menu).name);
fprintf(fid, ...
    'Task(%d).string = { ...\n',menu);
n = length(GID.task(menu).string(:,1));
for k=1:n
  fprintf(fid,' %d ''%s''  ''%s''  ''%s'' ...\n',GID.task(menu).string{k,1:4});
end
fprintf(fid,'};\n\n');
end

fprintf(fid, ...
['for menu=1:Nmenus\n' ...
 '    task(menu).name = Task(menu).name;\n' ...
 '    for n=1:length(Task(menu).string(:))/4\n' ...
 '        task(menu).string(n,1) = Task(menu).string(1+(n-1)*4);\n' ...
 '        task(menu).string(n,2) = Task(menu).string(2+(n-1)*4);\n' ...
 '        task(menu).string(n,3) = Task(menu).string(3+(n-1)*4);\n' ...
 '        task(menu).string(n,4) = Task(menu).string(4+(n-1)*4);\n' ...
 '    end\n' ...
 'end\n']);

fprintf(fid, ...
['  GID.guiName = ''%s'';\n' ...
 '  GID.Nmenus = Nmenus;\n' ...
 '  GID.frontImage = ''%s'';\n' ...
 '  GID.bgcolor = ''%s'';\n' ...
 '  GID.fontsize = %d;\n' ...
 '  GID.fontweight = ''%s'';\n' ...
 '  GID.position = [%4.2f %4.2f %4.2f %4.2f];\n' ...
 '  GID.resize = ''%s'';\n' ...
 '  GID.task = task;\n'], ...
 GID.guiName, GID.frontImage, GID.bgcolor, GID.fontsize, ...
 GID.fontweight, GID.position, GID.resize);

fprintf(fid,'\n');
fprintf(fid, ...
['%% The buildMainGUI builder needs fields label, type and ops\n' ...
 '%% The addLabel function can extract them from GID''s string field\n' ...
 '  GID = addLabel(GID);\n']);

fprintf(fid,'\n');
fprintf(fid, ...
['%% returns the final GID to the caller.\n' ...
 '  varargout{1} = GID;\n' ...
 '%% Also save GID to the following file.\n' ...
 '  save main_inputfile GID;\n' ...
 'end   %% if\n']);

author: Sable

repo: mcbench-benchmarks


text:
function [Y_train, label_train, Y_test, label_test] = pickTrainTest_2(dataset, N_train_c)
    data_fn = fullfile('data', strcat(dataset, '.mat'))
    load(data_fn);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ================== block: random projection ==========================
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     fprintf('random projection...')
%     d_new = 3000;
%     A = randn(d_new, size(Y,1))        ;
%     Y = A*Y;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ------------------end of block: random projection ----------------------------
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Y = normc(Y);
    d = size(Y,1);
    if ~exist('Y_range', 'var')
        Y_range = label_to_range(label);
    end

    C = numel(Y_range) - 1;
    N_total = Y_range(C+1);
    N_train = C*N_train_c;
    N_test = N_total - N_train;

    Y_train = zeros(d, N_train);
    Y_test = zeros(d, N_test);
    label_train = zeros(1, N_train);
    label_test = zeros(1, N_test);
    %%
    cur_train = 0;
    cur_test = 0;
    for c = 1: C
        Yc = get_block_col(Y, c, Y_range);
        N_total_c = size(Yc, 2);
        N_test_c = N_total_c - N_train_c;
        label_train(:, cur_train + 1: cur_train + N_train_c) = c*ones(1, N_train_c);
        label_test(:, cur_test + 1: cur_test + N_test_c) = c*ones(1, N_test_c);

        idx = randperm(N_total_c);

        Y_train(:, cur_train + 1: cur_train + N_train_c) = Yc(:, idx(1: N_train_c));
        Y_test(:, cur_test + 1: cur_test + N_test_c) = Yc(:, idx(N_train_c + 1: end));

        cur_train = cur_train + N_train_c;
        cur_test = cur_test + N_test_c;
    end

    % d_new = 3000;
    % if size(Y,1) > d_new
    %     fprintf('pca...')
    %     [W,frac] = pcam(Y',d_new);
    %     Y_train = normc((Y_train'*W)');
    %     Y_test = normc((Y_test'*W)');
    %     fprintf('DONE\n');
    % else
        Y_train = normc(Y_train);
        Y_test = normc(Y_test);
    % end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%%begin
% We will begin by constructing a HAL object with default properties

hal = HAL()

% Let's plot it in the zero angle configuration
qz = zeros(1,7);
hal.plot(qz); campos([6.5 6 5]);
% Note with a neutral wrist, the x-axis is in the direction of the
% forearm, and the z-axis is the elbow axis.

% We can see already that as a subclass, HAL inherits lots of
% functionality SerialLink objects. Forward kinematics is the same too
Th = hal.fkine(qz)

% Methods that are overloaded are ikine and islimit. Check out the help
% or doc files for these methods.
%
% Extending these is the method reachable, which combines the two
% inverse kinematic solutions and checks for joint limits. Poses that
% are unreachable are ouput as a NaN row.
[q1, q2] = hal.ikine(Th)
q = hal.reachable(q1, q2)

% You can use the method h2fsu to convert hand frames into forearm,
% swivel and upper arm frames (hence h2fsu). You can also use the hand
% point only - in which case the swivel angle must be explictly specified,
% or h2fsu uses the hand-to-goal method.
[Tf, Ts, Tu] = hal.h2fsu(Th(1:3,4))

% h2fsu in this case has used the hand-to-goal method of resolving the
% swivel angle. The elbow axis is made perpendicular to the goal vector
hal.goal
dot(hal.goal, Tf(1:3,3))
% The dot product is near-zero which says they are perpendicular.

% Let's see what it looks like. We will use plot3d this time, but first
% we must add the collision models to our HAL object.
[q1, q2] = hal.ikine(Tf)
q = hal.reachable(q1, q2)
hal.addCM;
hal.plot3d(q); campos([6.5 6 5]);

author: Allopart

repo: rbpf-gmapping


text:
function [xMin,fMin,exitCode]=steepestDescent(f,x0,epsilon,deltaTestDist,delta,lineSearchParams,maxIter)
%%STEEPESTDESCENT Perform unconstrained nonlinear optimization using
%              the steepest descent algorithm. The algorithm performs
%              unconstrained minimization of a nonlinear function without
%              having to provide a Hessian matrix.
%
%INPUTS: f A handle to the function (and its gradient) over which the
%          minimization is to be performed. The function [fVal,gVal]=f(x)
%          takes the NX1 x vector and returns the real scalar function
%          value fVal and gradient gVal at the point x.
%       x0 The NX1-dimensional point from which the minimization starts.
%  epsilon The parameter determining the accuracy of the desired solution
%          in terms of the gradient. The function terminates when
%          norm(g) < epsilon*max([1, norm(x)])
%          where g is the gradient. The default if omitted or an empty
%          matrix is passed is 1e-6.
% deltaTestDist The number of iterations back to use to compute the
%          decrease of the objective function if a delta-based convergence
%          test is performed. If zero, then no delta-based convergence
%          testing is done. The default if omitted or an empty matrix is
%          passed is zero.
%    delta The delta for the delta convergence test. This determines the
%          minimum rate of decrease of the objective function. Convergence
%          is determined if (f'-f)<=delta*f, where f' is the value of the
%          objective function f deltaTestDist iterations ago,and f is the
%          current objective function value. The default if this parameter
%          is omitted or an empty matrix is passed is 0.
% lineSearchParams An optional structure whose members specify tolerances
%          for the line search. The parameters are described as in the
%          lineSearch function, except if -1 is passed instead of a
%          parameter structure, then no line search is performed.
%  maxIter The maximum number of iterations to use for the algorithm. The
%          default if this parameter is omitted or an empty matrix is
%          passed is 100.
%
%OUTPUTS: xMin The value of x at the minimum point found. If exitCode is
%              negative, then this might be an empty matrix.
%         fMin The cost function value at the minimum point found. If
%              exitCode is negative, then this might be an empty matrix.
%     exitCode A value indicating the termination condition of the
%              algorithm. Nonnegative values indicate success; negative
%              values indicate some type of failure. Possible values are:
%                  0 The algorithm termiated successfully based on the
%                    gradient criterion.
%                  1 The algorithm terminated successfully based on the
%                    accuracy criterion.
%                 -1 The maximum number of overall iterations was reached.
%                 -2 A non-finite value was encoutered outside the line
%                    search.
%              Other negative values correspond to a failure in lineSearch
%              and correspond to the exitCode returned by the lineSearch
%              function.
%
%The algorithm is implemented based on the description of steepest descent
%given in Chapter 1.2 of [1].
%
%EXAMPLE:
%The example is that used in the lineSearch file.
% f=@(x)deal((x(1)+x(2)-3)*(x(1)+x(2))^3*(x(1)+x(2)-6)^4,... %The function
%            [(-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)));
%            (-6+x(1)+x(2))^3*(x(1)+x(2))^2*(54+8*x(1)^2+x(2)*(-45+8*x(2))+x(1)*(-45+16*x(2)))]);%And the gradient as the second return.
% %Note that the deal function is used to make an anonymous function have
% %two outputs.
% x0=[0.5;0.25];
% [xMin,fMin,exitCode]=steepestDescent(f,x0)
%The optimum point found is such that sum(xMin) is approximately
%1.73539450 with a minimum function value of approximately -2.1860756.
%
%REFERENCES:
%[1] D. P. Bertsekas, Nonlinear Programming, 3rd ed. Belmont, MA: Athena
%    Science, 2016.
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

    if(nargin<3||isempty(epsilon))
        epsilon=1e-6;
    end

    if(nargin<4||isempty(deltaTestDist))
        deltaTestDist=0;
    end

    if(nargin<5||isempty(delta))
        delta=0;
    end

    %The default line search parameter options
    useLineSearch=true;
    if(nargin>=6&&~isempty(lineSearchParams))
        if(isnumeric(lineSearchParams)&&lineSearchParams==-1)
            %If one just wishes to blindly take steps.
            useLineSearch=false;
        elseif(~isstruct(lineSearchParams)||(isnumeric(lineSearchParams)&&lineSearchParams~=-1))
            error('Unknown lineSearchParams value given')
        end
    else
        lineSearchParams=[];
    end

    if(nargin<7||isempty(maxIter))
        maxIter=100;
    end

    x=x0;
    [fValCur,gradF]=f(x0);

    if(deltaTestDist>0)
        pastFVals=zeros(deltaTestDist,1);
        pastFVals(1)=fValCur;
    end

    for cutIter=1:maxIter
         %Steepest descent direction.
         d=-gradF;

        if(useLineSearch)
            %Perform a line search in the given descent direction.
            [xNew,fValNew,gradFCur,~,~,exitCode]=lineSearch(f,x,d,[],[],lineSearchParams);
            if(isempty(xNew))
                %Just return the last valid values.
                xMin=x;
                fMin=fValCur;
                return;
            end
            xCur=xNew;
            fValCur=fValNew;
        else%Take a step without a line search.
            xCur=x+d;
            [fValCur,gradFCur]=f(xCur);
        end

        if(any(~isfinite(xCur(:))))
            exitCode=-2;
            xMin=[];
            fMin=[];
            return;
        end

        %Check for convergence based on the gradient.
        if(norm(gradFCur)<epsilon*max([1,norm(xCur)]))
            xMin=xCur;
            fMin=fValCur;
            exitCode=0;
            return;
        end

        %Check for convergence based on the actual function value.
        if(deltaTestDist~=0)
            if(pastFVals(end)-fValCur<=delta*fValCur)
                xMin=xCur;
                fMin=fValCur;
                exitCode=1;
                return;
            end
            pastFVals=circshift(pastFVals,[1,0]);
            pastFVals(1)=fValCur;
        end

        x=xCur;
        gradF=gradFCur;
    end

    %The maximum number of iterations elapsed without convergence
    xMin=xCur;
    fMin=fValCur;
    exitCode=-1;
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


function example_bar_plot()
test_data =[18     0; 20     0;   21     2;    30    14;    35    34;    40    57;    45    65;    50    46;    55     9;    60     2;    65     1;    70     0];

% Create figure
figure1 = figure('Color',[1 1 1]);

subplot(1,2,1)

hb=barh(test_data(:,1),test_data(:,2),'DisplayName','Test Data');

ylabel('parameter [units]');
xlabel('#');
legend('show','Location','northwest');
subplot(1,2,2)

hb=bar(test_data(:,1),test_data(:,2),'DisplayName','Test Data');

xlabel('parameter [units]');
ylabel('#');
legend('show','Location','northwest');

xdata=test_data(:,1);
barWidth=test_getBarWidthInAbsolutUnits(hb);

x_l=xdata-barWidth/2;
x_u=xdata+barWidth/2;
max_y=max(test_data(:,2))*1.2;
x=[];
y=[];
for i=1:length(x_l)
    x = [x , x_l(i),x_l(i),nan,x_u(i),x_u(i),nan];
    y = [y,       0,max_y ,nan,0     ,max_y ,nan];

end
hold on
plot(x,y,'r');

matlab2tikz('figurehandle',figure1,'filename','example_v_bar_plot.tex' ,'standalone', true);

    function BarWidth=test_getBarWidthInAbsolutUnits(h)
        % astimates the width of a bar plot
        XData_bar=get(h,'XData');
        length_bar = length(XData_bar);
        BarWidth= get(h, 'BarWidth');
        if length_bar > 1
            BarWidth = min(diff(XData_bar))*BarWidth;
        end


text:
function [D] = spm_DEM_eval_diff(x,v,qp,M,bilinear)
% evaluates derivatives for DEM schemes
% FORMAT [D]       = spm_DEM_eval_diff(x,v,qp,M,bilinear)
% v{i} - casual states
% x(i) - hidden states
% qp - conditional density of parameters
%  qp.p{i} - parameter deviates for i-th level
%  qp.u(i) - basis set
%  qp.x(i) - expansion point ( = prior expectation)
% M        - model structure
% bilinear - optional flag to suppress second-order derivatives
%
% D        - derivatives
%   D.dgdv
%   ...
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_eval_diff.m 6734 2016-03-02 12:02:46Z peter $

% check for evaluation of bilinear terms
%--------------------------------------------------------------------------
try
    bilinear;
catch
    bilinear = 1;
end

% get dimensions
%==========================================================================
nl    = size(M,2);                       % number of levels
ne    = sum(spm_vec(M.l));               % number of e (errors)
nx    = sum(spm_vec(M.n));               % number of x (hidden states)
np    = sum(spm_vec(M.p));               % number of p (parameters)
ny    = M(1).l;                          % number of y (inputs)
nc    = M(end).l;                        % number of c (prior causes)

% initialise cell arrays for hierarchical structure
%--------------------------------------------------------------------------
df.dv = cell(nl - 1,nl - 1);
df.dx = cell(nl - 1,nl - 1);
df.dp = cell(nl - 1,nl - 1);
dg.dv = cell(nl    ,nl - 1);
dg.dx = cell(nl    ,nl - 1);
dg.dp = cell(nl    ,nl - 1);

for i = 1:(nl - 1)
    dg.dv{i + 1,i} = sparse(M(i).m,M(i).m);
    dg.dx{i + 1,i} = sparse(M(i).m,M(i).n);
    dg.dp{i + 1,i} = sparse(M(i).m,M(i).p);
    dg.dv{i    ,i} = sparse(M(i).l,M(i).m);
    dg.dx{i    ,i} = sparse(M(i).l,M(i).n);
    dg.dp{i    ,i} = sparse(M(i).l,M(i).p);
    df.dv{i    ,i} = sparse(M(i).n,M(i).m);
    df.dx{i    ,i} = sparse(M(i).n,M(i).n);
    df.dp{i    ,i} = sparse(M(i).n,M(i).p);
end

if bilinear
    for i = 1:(nl - 1)
        dg.dvp{i}      = cell(M(i).p,1);
        dg.dxp{i}      = cell(M(i).p,1);
        df.dvp{i}      = cell(M(i).p,1);
        df.dxp{i}      = cell(M(i).p,1);
        [dg.dvp{i}{:}] = deal(dg.dv);
        [dg.dxp{i}{:}] = deal(dg.dx);
        [df.dvp{i}{:}] = deal(df.dv);
        [df.dxp{i}{:}] = deal(df.dx);
    end
end

% Derivatives at each hierarchical level
%==========================================================================

% inline function for evaluating projected parameters
%--------------------------------------------------------------------------
h     = @(f,x,v,q,u,p) f(x,v,spm_unvec(spm_vec(p) + u*q,p));
for i = 1:(nl - 1)

    % states level i
    %----------------------------------------------------------------------
    xvp = {x{i},v{i},qp.p{i},qp.u{i},M(i).pE};

    % 1st and 2nd partial derivatives (states)
    %----------------------------------------------------------------------
    if bilinear && np
        try
            [dgdxp, dgdx] = spm_diff(h,M(i).gx,xvp{:},4,'q');
            [dgdvp, dgdv] = spm_diff(h,M(i).gv,xvp{:},4,'q');
            [dfdxp, dfdx] = spm_diff(h,M(i).fx,xvp{:},4,'q');
            [dfdvp, dfdv] = spm_diff(h,M(i).fv,xvp{:},4,'q');
        catch
            [dgdxp, dgdx] = spm_diff(h,M(i).g,xvp{:},[2 4],'q');
            [dgdvp, dgdv] = spm_diff(h,M(i).g,xvp{:},[3 4],'q');
            [dfdxp, dfdx] = spm_diff(h,M(i).f,xvp{:},[2 4],'q');
            [dfdvp, dfdv] = spm_diff(h,M(i).f,xvp{:},[3 4],'q');
        end
    else
        try
            dgdx = h(M(i).gx,xvp{:});
            dgdv = h(M(i).gv,xvp{:});
            dfdx = h(M(i).fx,xvp{:});
            dfdv = h(M(i).fv,xvp{:});
        catch
            dgdx = spm_diff(h,M(i).g,xvp{:},2);
            dgdv = spm_diff(h,M(i).g,xvp{:},3);
            dfdx = spm_diff(h,M(i).f,xvp{:},2);
            dfdv = spm_diff(h,M(i).f,xvp{:},3);
        end
    end

    % 1st-order partial derivatives (parameters)
    %----------------------------------------------------------------------
    try
        dfdp   = h(M(i).fp,xvp{:});
        dgdp   = h(M(i).gp,xvp{:});
    catch
        dfdp   = spm_diff(h,M(i).f,xvp{:},4);
        dgdp   = spm_diff(h,M(i).g,xvp{:},4);
    end

%     % check which dervatives need to be evaluated
%     %====================================================================
%     D(i).dgdv = nnz(dgdv) + nnz(spm_vec(dgdvp));
%     D(i).dgdx = nnz(dgdx) + nnz(spm_vec(dgdxp));
%     D(i).dfdv = nnz(dfdv) + nnz(spm_vec(dfdvp));
%     D(i).dfdx = nnz(dfdx) + nnz(spm_vec(dfdxp));

    % Constant terms (linking causes over levels)
    %----------------------------------------------------------------------
    dg.dv{i + 1,i} = -speye(M(i).m,M(i).m);

    % place 1st derivatives in array
    %----------------------------------------------------------------------
    dg.dx{i,i} = dgdx;
    dg.dv{i,i} = dgdv;
    df.dx{i,i} = dfdx;
    df.dv{i,i} = dfdv;
    df.dp{i,i} = dfdp;
    dg.dp{i,i} = dgdp;

    % place 2nd derivatives in array
    %----------------------------------------------------------------------
    if bilinear && np
        for j = 1:length(dgdxp)
            dg.dxp{i}{j}{i,i} = dgdxp{j};
            dg.dvp{i}{j}{i,i} = dgdvp{j};
            df.dxp{i}{j}{i,i} = dfdxp{j};
            df.dvp{i}{j}{i,i} = dfdvp{j};
        end
    end
end

% concatenate hierarchical forms
%==========================================================================
D.dgdv  = spm_cat(dg.dv);
D.dgdx  = spm_cat(dg.dx);
D.dfdv  = spm_cat(df.dv);
D.dfdx  = spm_cat(df.dx);
D.dfdp  = spm_cat(df.dp);
D.dgdp  = spm_cat(dg.dp);

% fixed derivatives w.r.t. prediction errors and states
%--------------------------------------------------------------------------
D.dfdy  =  sparse(nx,ny);
D.dfdc  =  sparse(nx,nc);
D.dedy  =  spm_speye(ne,ny);
D.dedc  = -spm_speye(ne,nc,nc - ne);

% bilinear terms if required
%--------------------------------------------------------------------------
if bilinear
    D.dgdvp = {};
    D.dgdxp = {};
    D.dfdvp = {};
    D.dfdxp = {};
    for i = 1:length(dg.dvp)
        for j = 1:length(dg.dvp{i})
            D.dgdvp{end + 1} = spm_cat(dg.dvp{i}{j});
            D.dgdxp{end + 1} = spm_cat(dg.dxp{i}{j});
            D.dfdvp{end + 1} = spm_cat(df.dvp{i}{j});
            D.dfdxp{end + 1} = spm_cat(df.dxp{i}{j});
        end
    end
end

author: spm


text:
%% Principal Component Analysis and Partial Least Squares
% Principal Component Analysis (PCA) and Partial Least Squares (PLS) are
% widely used tools. This code is to show their relationship through the
% Nonlinear Iterative PArtial Least Squares (NIPALS) algorithm.

%% The Eigenvalue and Power Method
% The NIPALS algorithm can be derived from the Power method to solve the
% eigenvalue problem. Let x be the eigenvector of a square matrix, A,
% corresponding to the eignvalue s:
%
% $$Ax=sx$$
%
% Modifying both sides by A iteratively leads to
%
% $$A^kx=s^kx$$
%
% Now, consider another vectro y, which can be represented as a linear
% combination of all eigenvectors:
%
% $$y=\sum_i^n b_ix_i=Xb$$
%
% where
%
% $$X=\left[x_1\,\,\, \cdots\,\,\, x_n \right]$$
%
% and
%
% $$b = \left[b_1\,\,\, \cdots\,\,\, b_n \right]^T$$
%
% Modifying y by A gives
%
% $$Ay=AXb=XSb$$
%
% Where S is a diagnal matrix consisting all eigenvalues. Therefore, for
% a large enough k,
%
% $$A^ky=XS^kb\approx \alpha x_1$$
%
% That is the iteration will converge to the direction of x_1, which is the
% eigenvector corresponding to the eigenvalue with the maximum module.
% This leads to the following Power method to solve the eigenvalue problem.

A=randn(10,5);
% sysmetric matrix to ensure real eigenvalues
B=A'*A;
%find the column which has the maximum norm
[dum,idx]=max(sum(A.*A));
x=A(:,idx);
%storage to judge convergence
x0=x-x;
%convergence tolerant
tol=1e-6;
%iteration if not converged
while norm(x-x0)>tol
    %iteration to approach the eigenvector direction
    y=A'*x;
    %normalize the vector
    y=y/norm(y);
    %save previous x
    x0=x;
    %x is a product of eigenvalue and eigenvector
    x=A*y;
end
% the largest eigen value corresponding eigenvector is y
s=x'*x;
% compare it with those obtained with eig
[V,D]=eig(B);
[d,idx]=max(diag(D));
v=V(:,idx);
disp(d-s)
% v and y may be different in signs
disp(min(norm(v-y),norm(v+y)))

%% The NIPALS Algorithm for PCA
% The PCA is a dimension reduction technique, which is based on the
% following decomposition:
%
% $$X=TP^T+E$$
%
% Where X is the data matrix (m x n) to be analysed, T is the so called
% score matrix (m x a), P the loading matrix (n x a) and E the residual.
% For a given tolerance of residual, the number of principal components, a,
% can be much smaller than the orginal variable dimension, n.
% The above power algorithm can be extended to get T and P by iteratively
% subtracting A (in this case, X) by x*y' (in this case, t*p') until the
% given tolerance satisfied. This is the so called NIPALS algorithm.

% The data matrix with normalization
A=randn(10,5);
meanx=mean(A);
stdx=std(A);
X=(A-meanx(ones(10,1),:))./stdx(ones(10,1),:);
B=X'*X;
% allocate T and P
T=zeros(10,5);
P=zeros(5);
% tol for convergence
tol=1e-6;
% tol for PC of 95 percent
tol2=(1-0.95)*5*(10-1);
for k=1:5
    %find the column which has the maximum norm
    [dum,idx]=max(sum(X.*X));
    t=A(:,idx);
    %storage to judge convergence
    t0=t-t;
    %iteration if not converged
    while norm(t-t0)>tol
        %iteration to approach the eigenvector direction
        p=X'*t;
        %normalize the vector
        p=p/norm(p);
        %save previous t
        t0=t;
        %t is a product of eigenvalue and eigenvector
        t=X*p;
    end
    %subtracing PC identified
    X=X-t*p';
    T(:,k)=t;
    P(:,k)=p;
    if norm(X)<tol2
        break
    end
end
T(:,k+1:5)=[];
P(:,k+1:5)=[];
S=diag(T'*T);
% compare it with those obtained with eig
[V,D]=eig(B);
[D,idx]=sort(diag(D),'descend');
D=D(1:k);
V=V(:,idx(1:k));
fprintf('The number of PC: %i\n',k);
fprintf('norm of score difference between EIG and NIPALS: %g\n',norm(D-S));
fprintf('norm of loading difference between EIG and NIPALS: %g\n',norm(abs(V)-abs(P)));

%% The NIPALS Algorithm for PLS
% For PLS, we will have two sets of data: the independent X and dependent
% Y. The NIPALS algorithm can be used to decomposes both X and Y so that
%
% $$X=TP^T+E,\,\,\,\,Y=UQ^T+F,\,\,\,\,U=TB$$
%
% The regression, U=TB is solved through least sequares whilst the
% decompsition may not include all components. That is why the approach is
% called partial least squares. This algorithm is implemented in the PLS
% function.

%% Example: Discriminant PLS using the NIPALS Algorithm
% From Chiang, Y.Q., Zhuang, Y.M and Yang, J.Y, "Optimal Fisher
% discriminant analysis using the rank decomposition", Pattern Recognition,
% 25 (1992), 101--111.

% Three classes data, each has 50 samples and 4 variables.
x1=[5.1 3.5 1.4 0.2; 4.9 3.0 1.4 0.2; 4.7 3.2 1.3 0.2; 4.6 3.1 1.5 0.2;...
    5.0 3.6 1.4 0.2; 5.4 3.9 1.7 0.4; 4.6 3.4 1.4 0.3; 5.0 3.4 1.5 0.2; ...
    4.4 2.9 1.4 0.2; 4.9 3.1 1.5 0.1; 5.4 3.7 1.5 0.2; 4.8 3.4 1.6 0.2; ...
    4.8 3.0 1.4 0.1; 4.3 3.0 1.1 0.1; 5.8 4.0 1.2 0.2; 5.7 4.4 1.5 0.4; ...
    5.4 3.9 1.3 0.4; 5.1 3.5 1.4 0.3; 5.7 3.8 1.7 0.3; 5.1 3.8 1.5 0.3; ...
    5.4 3.4 1.7 0.2; 5.1 3.7 1.5 0.4; 4.6 3.6 1.0 0.2; 5.1 3.3 1.7 0.5; ...
    4.8 3.4 1.9 0.2; 5.0 3.0 1.6 0.2; 5.0 3.4 1.6 0.4; 5.2 3.5 1.5 0.2; ...
    5.2 3.4 1.4 0.2; 4.7 3.2 1.6 0.2; 4.8 3.1 1.6 0.2; 5.4 3.4 1.5 0.4; ...
    5.2 4.1 1.5 0.1; 5.5 4.2 1.4 0.2; 4.9 3.1 1.5 0.2; 5.0 3.2 1.2 0.2; ...
    5.5 3.5 1.3 0.2; 4.9 3.6 1.4 0.1; 4.4 3.0 1.3 0.2; 5.1 3.4 1.5 0.2; ...
    5.0 3.5 1.3 0.3; 4.5 2.3 1.3 0.3; 4.4 3.2 1.3 0.2; 5.0 3.5 1.6 0.6; ...
    5.1 3.8 1.9 0.4; 4.8 3.0 1.4 0.3; 5.1 3.8 1.6 0.2; 4.6 3.2 1.4 0.2; ...
    5.3 3.7 1.5 0.2; 5.0 3.3 1.4 0.2];
x2=[7.0 3.2 4.7 1.4; 6.4 3.2 4.5 1.5; 6.9 3.1 4.9 1.5; 5.5 2.3 4.0 1.3; ...
    6.5 2.8 4.6 1.5; 5.7 2.8 4.5 1.3; 6.3 3.3 4.7 1.6; 4.9 2.4 3.3 1.0; ...
    6.6 2.9 4.6 1.3; 5.2 2.7 3.9 1.4; 5.0 2.0 3.5 1.0; 5.9 3.0 4.2 1.5; ...
    6.0 2.2 4.0 1.0; 6.1 2.9 4.7 1.4; 5.6 2.9 3.9 1.3; 6.7 3.1 4.4 1.4; ...
    5.6 3.0 4.5 1.5; 5.8 2.7 4.1 1.0; 6.2 2.2 4.5 1.5; 5.6 2.5 3.9 1.1; ...
    5.9 3.2 4.8 1.8; 6.1 2.8 4.0 1.3; 6.3 2.5 4.9 1.5; 6.1 2.8 4.7 1.2; ...
    6.4 2.9 4.3 1.3; 6.6 3.0 4.4 1.4; 6.8 2.8 4.8 1.4; 6.7 3.0 5.0 1.7; ...
    6.0 2.9 4.5 1.5; 5.7 2.6 3.5 1.0; 5.5 2.4 3.8 1.1; 5.5 2.4 3.7 1.0; ...
    5.8 2.7 3.9 1.2; 6.0 2.7 5.1 1.6; 5.4 3.0 4.5 1.5; 6.0 3.4 4.5 1.6; ...
    6.7 3.1 4.7 1.5; 6.3 2.3 4.4 1.3; 5.6 3.0 4.1 1.3; 5.5 2.5 5.0 1.3; ...
    5.5 2.6 4.4 1.2; 6.1 3.0 4.6 1.4; 5.8 2.6 4.0 1.2; 5.0 2.3 3.3 1.0; ...
    5.6 2.7 4.2 1.3; 5.7 3.0 4.2 1.2; 5.7 2.9 4.2 1.3; 6.2 2.9 4.3 1.3; ...
    5.1 2.5 3.0 1.1; 5.7 2.8 4.1 1.3];
x3=[6.3 3.3 6.0 2.5; 5.8 2.7 5.1 1.9; 7.1 3.0 5.9 2.1; 6.3 2.9 5.6 1.8; ...
    6.5 3.0 5.8 2.2; 7.6 3.0 6.6 2.1; 4.9 2.5 4.5 1.7; 7.3 2.9 6.3 1.8; ...
    6.7 2.5 5.8 1.8; 7.2 3.6 6.1 2.5; 6.5 3.2 5.1 2.0; 6.4 2.7 5.3 1.9; ...
    6.8 3.0 5.5 2.1; 5.7 2.5 5.0 2.0; 5.8 2.8 5.1 2.4; 6.4 3.2 5.3 2.3; ...
    6.5 3.0 5.5 1.8; 7.7 3.8 6.7 2.2; 7.7 2.6 6.9 2.3; 6.0 2.2 5.0 1.5; ...
    6.9 3.2 5.7 2.3; 5.6 2.8 4.9 2.0; 7.7 2.8 6.7 2.0; 6.3 2.7 4.9 1.8; ...
    6.7 3.3 5.7 2.1; 7.2 3.2 6.0 1.8; 6.2 2.8 4.8 1.8; 6.1 3.0 4.9 1.8; ...
    6.4 2.8 5.6 2.1; 7.2 3.0 5.8 1.6; 7.4 2.8 6.1 1.9; 7.9 3.8 6.4 2.0; ...
    6.4 2.8 5.6 2.2; 6.3 2.8 5.1 1.5; 6.1 2.6 5.6 1.4; 7.7 3.0 6.1 2.3; ...
    6.3 3.4 5.6 2.4; 6.4 3.1 5.5 1.8; 6.0 3.0 4.8 1.8; 6.9 3.1 5.4 2.1; ...
    6.7 3.1 5.6 2.4; 6.9 3.1 5.1 2.3; 5.8 2.7 5.1 1.9; 6.8 3.2 5.9 2.3; ...
    6.7 3.3 5.7 2.5; 6.7 3.0 5.2 2.3; 6.3 2.5 5.0 1.9; 6.5 3.0 5.2 2.0; ...
    6.2 3.4 5.4 2.3; 5.9 3.0 5.1 1.8];
%Split data set into training (1:25) and testing (26:50)
idxTrain = 1:25;
idxTest = 26:50;

% Combine training data with normalization
X = [x1(idxTrain,:);x2(idxTrain,:);x3(idxTrain,:)];
% Define class indicator as Y
Y = kron(eye(3),ones(25,1));
% Normalization
xmean = mean(X);
xstd = std(X);
ymean = mean(Y);
ystd = std(Y);
X = (X - xmean(ones(75,1),:))./xstd(ones(75,1),:);
Y = (Y - ymean(ones(75,1),:))./ystd(ones(75,1),:);
% Tolerance for 90 percent score
tol = (1-0.9) * 25 * 4;
% Perform PLS
[T,P,U,Q,B] = pls(X,Y,tol);
% Results
fprintf('Number of components retained: %i\n',size(B,1))
% Predicted classes
X1 = (x1(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
X2 = (x2(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
X3 = (x3(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
Y1 = X1 * (P*B*Q');
Y2 = X2 * (P*B*Q');
Y3 = X3 * (P*B*Q');
Y1 = Y1 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
Y2 = Y2 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
Y3 = Y3 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
% Class is determined from the cloumn which is most close to 1
[dum,classid1]=min(abs(Y1-1),[],2);
[dum,classid2]=min(abs(Y2-1),[],2);
[dum,classid3]=min(abs(Y3-1),[],2);
bar(1:25,classid1,'b');
hold on
bar(26:50,classid2,'r');
bar(51:75,classid3,'g');
hold off

% The results show that most samples are classified correctly.


text:
classdef NSGAIIARSBX < ALGORITHM
% <multi> <real/integer> <constrained/none>
% NSGA-II with adaptive rotation based simulated binary crossover

%------------------------------- Reference --------------------------------
% L. Pan, W. Xu, L. Li, C. He, and R. Cheng, Adaptive simulated binary
% crossover for rotated multi-objective optimization, Swarm and
% Evolutionary Computation, 2021, 60: 100759.
%--------------------------------------------------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N);
            B  = eye(Problem.D);
            m  = 0.5*(Problem.upper - Problem.lower);
            ps = 0.5;

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = ARSBX(Problem,Population(MatingPool),{B,m,ps});
                [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N);
                [B,m,ps,Population] = UpdateParameter(Problem,Population);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef SVORLin < Algorithm
    %SVORLin Linear Support Vector for Ordinal Regression (Implicit constraints)
    %   This class derives from the Algorithm Class and implements the
    %   linear SVORIM method. This class uses SVORIM implementation by
    %   W. Chu et al (http://www.gatsby.ucl.ac.uk/~chuwei/svor.htm)
    %
    %   SVORLin methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %     [1] W. Chu and S. S. Keerthi, Support Vector Ordinal Regression,
    %         Neural Computation, vol. 19, no. 3, pp. 792–815, 2007.
    %         http://10.1162/neco.2007.19.3.792
    %     [2] P.A. Gutiérrez, M. Pérez-Ortiz, J. Sánchez-Monedero,
    %         F. Fernández-Navarro and C. Hervás-Martínez
    %         Ordinal regression methods: survey and experimental study
    %         IEEE Transactions on Knowledge and Data Engineering, Vol. 28. Issue 1
    %         2016
    %         http://dx.doi.org/10.1109/TKDE.2015.2457911
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html
    properties
        description = 'Support Vector for Ordinal Regression (Implicit constraints / Linear)';
        parameters = struct('C', 0.1);
    end

    properties (Access = private)
        algorithmMexPath = fullfile(fileparts(which('Algorithm.m')),'SVORIM');
    end

    methods

        function obj = SVORLin(varargin)
            %SVORLin constructs an object of the class SVORLin and sets its default
            %   characteristics
            %   OBJ = SVORLin() builds SVORLin object
        end

        function [projectedTrain,predictedTrain] = privfit(obj,train,parameters)
            %PRIVFIT trains the model for the SVORLin method with TRAIN data and
            %vector of parameters PARAMETERS.
            if isempty(strfind(path,obj.algorithmMexPath))
                addpath(obj.algorithmMexPath);
            end
            [alpha, thresholds, projectedTrain] = svorim([train.patterns train.targets],1,parameters.C,0,0,1);
            predictedTrain = obj.assignLabels(projectedTrain, thresholds);
            model.projection = alpha;
            model.thresholds = thresholds;
            model.parameters = parameters;
            model.train = train.patterns;
            obj.model = model;
            projectedTrain = projectedTrain';
            if ~isempty(strfind(path,obj.algorithmMexPath))
                rmpath(obj.algorithmMexPath);
            end
        end

        function [projected, predicted] = privpredict(obj, test)
            %PREDICT predicts labels of TEST patterns labels. The object needs to be fitted to the data first.
            kernelMatrix = computeKernelMatrix(obj.model.train',test','linear',1);
            projected = obj.model.projection*kernelMatrix;

            predicted = SVORLin.assignLabels(projected, obj.model.thresholds);
            projected = projected';
        end
    end
    methods (Static = true)
        function predicted = assignLabels(projected, thresholds)
            numClasses = size(thresholds,2)+1;
            %TEST assign the labels from projections and thresholds
            project2 = repmat(projected, numClasses-1,1);
            project2 = project2 - thresholds'*ones(1,size(project2,2));

            % Asignation of the class
            % f(x) = max {Wx-bk<0} or Wx - b_(K-1) > 0
            wx=project2;

            % The procedure for that is the following:
            % We assign the values > 0 to NaN
            wx(wx(:,:)>0)=NaN;

            % Then, we choose the biggest one.
            [maximum,predicted]=max(wx,[],1);

            % If a max is equal to NaN is because Wx-bk for all k is >0, so this
            % pattern belongs to the last class.
            predicted(isnan(maximum(:,:)))=numClasses;

            predicted = predicted';
        end

    end

end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
function [oa, pa, K, CM] = USFE_MSTV(HSI, Tr, Te, dim, Trees)
[FE_MSTV]= MSTV_Xu(HSI, dim);
[acc_Mean,acc_std,CM]=RF_ntimes_overal(FE_MSTV,Tr,Te,Trees);
pa=acc_Mean(1:dim,1);
oa=acc_Mean(dim+2,1);
K=acc_Mean(dim+3,1);

author: BehnoodRasti


text:
function snc = vol2talairachVolume(vw,skipTalFlag)
% Function to give the mean Talairach and MNI coordinates of an ROI.  If no
% output argument is given, a message box will pop up.  Otherwise, the
% coordinates (a vector of 3 numbers) will be output in snc.
%
%    snc = vol2talairachVolume(vw,[skipTalFlag = 0])
%
% skipTalFlag tells the function to ignore the Talairach space.  This is
% useful for scripting across subjects where a Talairach normalization has
% not been created and you just want the spatialNorm (MNI).
%
% Note, the coordinate that you are given is for the current ROI and is
% based on the mean MNI coordinate (after transforming all coordinates to
% MNI).  If you want the peak activation within an ROI, then find that
% point, make a point ROI there, and then call this function for that ROI.
%

global mrSESSION;

if notDefined('skipTalFlag'), skipTalFlag = 0; end

[talairach, spatialNorm] = loadTalairachXform(mrSESSION.subject,[],skipTalFlag);

roiNum = viewGet(vw, 'curROI');
coords = vw.ROIs(roiNum).coords';
name = vw.ROIs(roiNum).name;

if(~isempty(talairach))
    talairach = volToTalairach(coords,talairach.vol2Tal);
    t = round(mean(talairach,1));
    msg = sprintf('ROI name:\t%s\nTalairach:\t%d, %d, %d',name,t(1),t(2),t(3));
else
    msg = [];
end

if(~isempty(spatialNorm))
    if(isfield(spatialNorm,'invLUT'))
        % We need to flip the coords from vAnat space to nifti space
        anatSz = viewGet(vw, 'anat size');
        c = [coords(:,3) anatSz(2)-coords(:,2) anatSz(3)-coords(:,1)];
        sz = size(squeeze(spatialNorm.invLUT.coordLUT(:,:,:,1)));
        inds = sub2ind(sz, c(:,1), c(:,2), c(:,3));
        ssCoords = [spatialNorm.invLUT.coordLUT(inds) spatialNorm.invLUT.coordLUT(inds+prod(sz)) spatialNorm.invLUT.coordLUT(inds+prod(sz)*2)];
        snc = round(mean(ssCoords,1));  % gets mean of all ROI coordinates
    elseif(isfield(spatialNorm,'voxToTemplateLUT'))
        % Support old-style xform:
        sz = size(squeeze(spatialNorm.voxToTemplateLUT(1,:,:,:)));
        inds = sub2ind(sz, coords(:,1), coords(:,2), coords(:,3));
        snc = round(mean(spatialNorm.voxToTemplateLUT(:,inds)',1));
    end
    [junk,tName] = fileparts(spatialNorm.sn.VG.fname);

    msg = [msg sprintf('\nSpatial Norm (%s):\t%d, %d, %d\n', tName, snc)];
    msg = [msg sprintf('MNI2Tal on SN:\t%d, %d, %d', round(mni2tal(snc)))];
    %msg = [msg tName];
end

if nargout<1  % only do message box if no output is called for
    msgbox(msg, 'Talairach');
end
fprintf('\n%s\n',msg);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
u0 = double(imread('58468899e2951a4b87b00897.png'));
u0 = mean(u0,3);

[M N] = size(u0);
MN = M*N;
I=reshape([1:M*N],M,N);
east=[I(:,2:end), I(:,end)];
north=[I(2:end,:); I(end,:)];
D1 = sparse(I,east,1,MN,MN) -speye(MN,MN);
D2 = sparse(I,north,1,MN,MN) -speye(MN,MN);
D = [D1 ; D2];

h0 = [ 1 2 1 ];
h = h0
for i=1:5 %% level of blur
    h = conv(h0,h);
end
size(h)
h = h'*h;
h = h/sum(sum(h));

stddev=5;
g = filter2(h,u0,'valid');
[Ms Ns]=size(g);
g = g + stddev*randn(Ms,Ns);

figure(1);
imagesc(u0); colormap(gray); %drawnow();
figure(2);
imagesc(g); colormap(gray); drawnow();

u=zeros(size(u0));
L2 = 8;
tau = 1; %% max value for tau = 1;
sig = 1/tau/L2;
lambda = .2;

p = zeros(2*MN,1);
%% primal-dual
for i=1:1000
    um = u;
    u = u-tau*(conv2(h,filter2(h,u,'valid')-g,'full') + reshape(D'*p,M,N));
    um = 2*u-um;
    p = p + sig*D*um(:);
    no= max(1,hypot(p(1:MN),p(MN+1:end))/lambda);
    p = p./[no;no];
    if mod(i,20)==0
        i
        figure(3);
        imagesc(u);
        colormap(gray); drawnow();
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function cycle_brent_test ( )

%*****************************************************************************80
%
%% CYCLE_BRENT_TEST tests the CYCLE_BRENT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLE_BRENT_PRB\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the CYCLE_BRENT library.\n' );

  cycle_brent_test01 ( );
  cycle_brent_test02 ( );
  cycle_brent_test03 ( );
  cycle_brent_test04 ( );
  cycle_brent_test05 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLE_BRENT_PRB\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function net = gtm(dim_latent, nlatent, dim_data, ncentres, rbfunc, ...
   prior)
%GTM	Create a Generative Topographic Map.
%
%	Description
%
%	NET = GTM(DIMLATENT, NLATENT, DIMDATA, NCENTRES, RBFUNC), takes the
%	dimension of the latent space DIMLATENT, the number of data points
%	sampled in the latent space NLATENT, the dimension of the data space
%	DIMDATA, the number of centres in the RBF model NCENTRES, the
%	activation function for the RBF RBFUNC and returns a data structure
%	NET. The parameters in the RBF and GMM sub-models are set by calls to
%	the corresponding creation routines RBF and GMM.
%
%	The fields in NET are
%	  type = 'gtm'
%	  nin = dimension of data space
%	  dimlatent = dimension of latent space
%	  rbfnet = RBF network data structure
%	  gmmnet = GMM data structure
%	  X = sample of latent points
%
%	NET = GTM(DIMLATENT, NLATENT, DIMDATA, NCENTRES, RBFUNC, PRIOR),
%	sets a Gaussian zero mean prior on the parameters of the RBF model.
%	PRIOR must be a scalar and represents the inverse variance of the
%	prior distribution.  This gives rise to a weight decay term in the
%	error function.
%
%	See also
%	GTMFWD, GTMPOST, RBF, GMM
%

%	Copyright (c) Ian T Nabney (1996-2001)

net.type = 'gtm';
% Input to functions is data
net.nin = dim_data;
net.dim_latent = dim_latent;

% Default is no regularisation
if nargin == 5
   prior = 0.0;
end

% Only allow scalar prior
if isstruct(prior) | size(prior) ~= [1 1]
   error('Prior must be a scalar');
end

% Create RBF network
net.rbfnet = rbf(dim_latent, ncentres, dim_data, rbfunc, ...
   'linear', prior);

% Mask all but output weights
net.rbfnet.mask = rbfprior(rbfunc, dim_latent, ncentres, dim_data);

% Create field for GMM output model
net.gmmnet = gmm(dim_data, nlatent, 'spherical');

% Create empty latent data sample
net.X = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function out = reconSpyr(coeff, filter)

% Reconstruct the image from steerable pyramid
% Input:
%        coeff:     the cell structured pyramid
%       filter:     filter, typically sp3.mat, sp1.mat, sp5.mat
%
% Output:
%          out:     reconstructed image

load(filter,'lo0filt','hi0filt','lofilt','bfilts');

res = reconSpyrLevs(coeff(2:length(coeff)),lofilt, bfilts);
temp = upConv(res, lo0filt,[1 1]);

highpass = upConv( coeff{1}, hi0filt,[1 1]);

out=highpass+temp;


text:
% pop_interp() - interpolate data channels
%
% Usage: EEGOUT = pop_interp(EEG, badchans, method);
%
% Inputs:
%     EEG      - EEGLAB dataset
%     badchans - [integer array] indices of channels to interpolate.
%                For instance, these channels might be bad.
%                [chanlocs structure] channel location structure containing
%                either locations of channels to interpolate or a full
%                channel structure (missing channels in the current
%                dataset are interpolated).
%     method   - [string] method used for interpolation (default is 'spherical').
%                'invdist' uses inverse distance on the scalp
%                'spherical' uses superfast spherical interpolation.
%                'spacetime' uses griddata3 to interpolate both in space
%                and time (very slow and cannot be interupted).
% Output:
%     EEGOUT   - data set with bad electrode data replaced by
%                interpolated data
%
% Author: Arnaud Delorme, CERCO, CNRS, 2009-

% Copyright (C) Arnaud Delorme, CERCO, 2009, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG com] = pop_interp(EEG, bad_elec, method)

    com = '';
    if nargin < 1
        help pop_interp;
        return;
    end;

    if nargin < 2
        disp('Warning: interpolation can be done on the fly in studies');
        disp('         this function will actually create channels in the dataset');
        disp('Warning: do not interpolate channels before running ICA');
        disp('You may define channel location to interpolate in the channel');
        disp('editor and declare such channels as non-data channels');

        enablenondat = 'off';
        if isfield(EEG.chaninfo, 'nodatchans')
            if ~isempty(EEG.chaninfo.nodatchans)
                enablenondat = 'on';
            end;
        end;
        cb_nondat   = [ 'tmpchaninfo = EEG.chaninfo; [chanlisttmp chanliststr] = pop_chansel( { tmpchaninfo.nodatchans.labels } );' ...
                      'if ~isempty(chanlisttmp),' ...
                      '   set(gcbf, ''userdata'', EEG.chaninfo.nodatchans(chanlisttmp));' ...
                      '   set(findobj(gcbf, ''tag'', ''chanlist''), ''string'', chanliststr);' ...
                      'end;' ...
                      'clear chanlisttmp chanliststr tmpchaninfo;' ];

        cb_otherdat = [ 'tmpanswer = inputdlg2({ ''Dataset index'' }, ''Choose dataset'', 1, { '''' });' ...
                      'if ~isempty(tmpanswer),' ...
                      '   tmpanswernum = round(str2num(tmpanswer{1}));' ...
                      '   if ~isempty(tmpanswernum),' ...
                      '       if tmpanswernum > 0 & tmpanswernum < length(ALLEEG),' ...
                      '           tmpchans1 = ALLEEG(tmpanswernum).chanlocs;' ...
                      '           tmpchans2 = EEG.chanlocs;' ...
                      '           tmpchanlist = setdiff( { tmpchans1.labels }, { tmpchans2.labels } );' ...
                      '           if ~isempty(tmpchanlist),' ...
                      '               set(gcbf, ''userdata'', ALLEEG(tmpanswernum).chanlocs);' ...
                      '               tmpchanlist(2,:) = { '' '' };' ...
                      '               set(findobj(gcbf, ''tag'', ''chanlist''), ''string'', [ tmpchanlist{:} ]);' ...
                      '           else,' ...
                      '               warndlg2(''No new channels in this dataset'');' ...
                      '           end;' ...
                      '       else,' ...
                      '           warndlg2(''Wrong index'');' ...
                      '       end;' ...
                      '   end;' ...
                      'end;' ...
                      'clear tmpanswer tmpanswernum;' ];

        uilist = { { 'Style' 'text' 'string' 'What channel(s) to interpolate' 'fontweight' 'bold' } ...
                   { 'style' 'text' 'string' 'none' 'tag' 'chanlist' } ...
                   { 'style' 'pushbutton' 'string' 'Select from non-data channels' 'callback' 'pop_interp(''nondatchan'',gcbf);' 'enable' enablenondat } ...
                   { 'style' 'pushbutton' 'string' 'Select from other dataset' 'callback' 'pop_interp(''selectchan'',gcbf);'} ...
                   { 'style' 'pushbutton' 'string' 'Use list of other dataset' 'callback' 'pop_interp(''uselist'',gcbf);'} ...
                   { } ...
                   { 'style' 'text'  'string' 'Interpolation method'} ...
                   { 'style' 'popupmenu'  'string' 'Spherical|Planar (slow)'  'tag' 'method' } ...
                   };

        geom = { 1 1 1 1 1 1 [1.8 1] };
        [res userdata tmp restag ] = inputgui( 'uilist', uilist, 'title', 'Interpolate channel(s) -- pop_interp()', 'geometry', geom, 'helpcom', 'pophelp(''pop_interp'')');
        if isempty(res) | isempty(userdata), return; end;

        if restag.method == 1
             method = 'spherical';
        else method = 'invdist';
        end;
        bad_elec = userdata.chans;

        com = sprintf('EEG = pop_interp(EEG, %s, ''%s'');', userdata.chanstr, method);
        if ~isempty(findstr('nodatchans', userdata.chanstr))
            eval( [ userdata.chanstr '=[];' ] );
        end;

    elseif isstr(EEG)
        command = EEG;
        clear EEG;
        fig = bad_elec;
        userdata = get(fig, 'userdata');

        if strcmpi(command, 'nondatchan')
            global EEG;
            tmpchaninfo = EEG.chaninfo;
            [chanlisttmp chanliststr] = pop_chansel( { tmpchaninfo.nodatchans.labels } );
            if ~isempty(chanlisttmp),
                userdata.chans   = EEG.chaninfo.nodatchans();
                userdata.chanstr = [ 'EEG.chaninfo.nodatchans([' num2str(chanlisttmp) '])' ];
                set(fig, 'userdata', userdata);
                set(findobj(fig, 'tag', 'chanlist'), 'string', chanliststr);
            end;
        else
            global ALLEEG EEG;
            tmpanswer = inputdlg2({ 'Dataset index' }, 'Choose dataset', 1, { '' });
            if ~isempty(tmpanswer),
                tmpanswernum = round(str2num(tmpanswer{1}));
                if ~isempty(tmpanswernum),
                    if tmpanswernum > 0 & tmpanswernum < length(ALLEEG),
                        TMPEEG = ALLEEG(tmpanswernum);

                        tmpchans1 = TMPEEG.chanlocs;
                        if strcmpi(command, 'selectchan')
                            chanlist = pop_chansel( { tmpchans1.labels } );
                        else
                            chanlist = 1:length(TMPEEG.chanlocs); % use all channels
                        end

                        % look at what new channels are selected
                        tmpchans2 = EEG.chanlocs;
                        [tmpchanlist chaninds] = setdiff( { tmpchans1(chanlist).labels }, { tmpchans2.labels } );
                        if ~isempty(tmpchanlist),
                            if length(chanlist) == length(TMPEEG.chanlocs)
                                userdata.chans   = TMPEEG.chanlocs;
                                userdata.chanstr = [ 'ALLEEG(' tmpanswer{1} ').chanlocs' ];
                            else
                                userdata.chans   = TMPEEG.chanlocs(chanlist(sort(chaninds)));
                                userdata.chanstr = [ 'ALLEEG(' tmpanswer{1} ').chanlocs([' num2str(chanlist(sort(chaninds))) '])' ];
                            end;
                            set(fig, 'userdata', userdata);
                            tmpchanlist(2,:) = { ' ' };
                            set(findobj(gcbf, 'tag', 'chanlist'), 'string', [ tmpchanlist{:} ]);
                        else
                            warndlg2('No new channels selected');
                        end;
                    else
                        warndlg2('Wrong index');
                    end;
                end;
            end;
        end;
        return;
    end;

    EEG = eeg_interp(EEG, bad_elec, method);

author: PatternRecognition

repo: OpenBMI


text:
function grid_level = index_level_own ( level, level_max, dim_num, ...
  point_num, grid_index, grid_base )

%*****************************************************************************80
%
%% INDEX_LEVEL_OWN: determine first level at which given index is generated.
%
%  Discussion:
%
%    We are constructing a sparse grid based on a 1D OWN rule (Gauss Legendre
%    or Gauss Hermite).  The grid is built up of product grids,
%    with a characteristic LEVEL.
%
%    We are concerned with identifying points in this product grid which
%    have actually been generated previously, on a lower value of LEVEL.
%
%    This routine determines the lowest value of LEVEL at which each of
%    the input points would be generated.
%
%    In 1D, given LEVEL, the number of points is ORDER = 2**(LEVEL+1) + 1,
%    (except that LEVEL = 0 implies ORDER = 1%), the BASE is (ORDER-1)/2,
%    and the point INDEX values range from -BASE to +BASE.
%
%    The values of INDEX and BASE allow us to determine the abstract
%    properties of the point.  In particular, if INDEX is 0, the corresponding
%    abscissa is 0, the special "nested" value we need to take care of.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer LEVEL, the level at which these points were
%    generated.  LEVEL_MIN <= LEVEL <= LEVEL_MAX.
%
%    Input, integer LEVEL_MAX, the maximum level.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points to be tested.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), the indices of the
%    points to be tested.
%
%    Input, integer GRID_BASE(DIM_NUM), the "base", which is essentially
%    the denominator of the index.
%
%    Output, integer GRID_LEVEL(POINT_NUM), the value of LEVEL at
%    which the point would first be generated.  This will be the same as
%    the input value of LEVEL, unless the point has an INDEX of 0 and
%    a corresponding BASE that is NOT zero.
%
  if ( dim_num == 1 )
    level_min = level_max;
  else
    level_min = 0;
  end
%
%  If a point has a DIM-th component whose INDEX is 0, then the
%  value of LEVEL at which this point would first be generated is
%  less than LEVEL, unless the DIM-th component of GRID_BASE is 0.
%
  for point = 1 : point_num

    grid_level(point) = max ( level, level_min );

    for dim = 1 : dim_num
      if ( grid_index(dim,point) == 0 )
        grid_level(point) = max ( grid_level(point) - grid_base(dim), level_min );
      end
    end

  end

  return
end


text:
function table = p_power_product ( p, e )

%*****************************************************************************80
%
%% P_POWER_PRODUCT: power products for Legendre polynomial P(n,x).
%
%  Discussion:
%
%    Let P(n,x) represent the Legendre polynomial of degree i.
%
%    For polynomial chaos applications, it is of interest to know the
%    value of the integrals of products of powers of X with every possible pair
%    of basis functions.  That is, we'd like to form
%
%      Tij = Integral ( -1 <= X <= +1 ) X^E * P(i,x) * P(j,x) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polyonomial factors.
%    0 <= P.
%
%    Input, integer E, the exponent of X in the integrand.
%    0 <= E.
%
%    Output, real TABLE(P+1,P+1), the table of integrals.
%
  table(1:p+1,1:p+1) = 0.0;

  order = p + 1 + floor ( ( e + 1 ) / 2 );
  [ x_table, w_table ] = p_quadrature_rule ( order );

  for k = 1 : order

    x = x_table(k);
    l_table = p_polynomial_value ( 1, p, x );
%
%  The following formula is an outer product in L_TABLE.
%
    if ( e == 0 )
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) *       l_table(1:p+1)' * l_table(1:p+1);
    else
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) * x^e * l_table(1:p+1)' * l_table(1:p+1);
    end

  end

  return
end


text:
function name = dfgetfitname()

% $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:38 $
% Copyright 2003-2004 The MathWorks, Inc.

count=dfgetset('fitcount');
if isempty(count)
    count = 1;
end
taken = 1;
while taken
    name=sprintf('fit %i', count);
    if isempty(find(getfitdb,'name',name))
        taken = 0;
    else
        count=count+1;
    end
end
dfgetset('fitcount',count+1);

author: zouchuhang

repo: LayoutNet


text:
function [listDuplicates] = check4DuplicatesInList(list)
% This function checks for duplicate entries in a list.
%
% INPUT
% list              List of e.g. metabolite abbr
%
% OUTPUT
% listDuplicates    List of duplicated entries. Second (or more) occurance
%                   of the duplicate is provided.
%
% Ines Thiele, 09/2021

listDuplicates = [];
cnt = 1;
[D,IA,ID]= duplicates((list));
Di = find(IA);
if isempty(Di)
    fprintf('No duplicate metabolites exists.\n')
else
    for i = 1 : length(Di)
        fprintf([list{Di(i)} ' appears more than once.\n'])
        listDuplicates{cnt,1} = list{Di(i)};
        listDuplicates{cnt,2} = num2str(Di(i));
        cnt = cnt +1;
    end
end

author: opencobra

repo: cobratoolbox


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function [H_LC, P, min_cost, final_result, results] = optimizeIoU(opt, X, Y, intrinsic, display)
    % prepare for random initilization (not use for now)
    for i = 1:1
        theta_x = optimvar('theta_x', 1, 1,'LowerBound',-180,'UpperBound',180); % 1x1
        theta_y = optimvar('theta_y', 1, 1,'LowerBound',-180,'UpperBound',180); % 1x1
        theta_z = optimvar('theta_z', 1, 1,'LowerBound',-180,'UpperBound',180); % 1x1
        T = optimvar('T', 1, 3,'LowerBound',-0.5,'UpperBound',0.5);
%         theta_x = 0;
%         theta_y = 10;
%         theta_z = 20;
%         T = [ 1 2 3];  % 1x3
%         costIoU(theta_x, theta_y, theta_z, T, X, Y, intrinsic);

        prob = optimproblem;
        f = fcn2optimexpr(@costIoU, theta_x, theta_y, theta_z, T, X, Y, intrinsic);
        prob.Objective = f;

        noise = (rand([1,6]) - 0.5) * 2; % -1 to 1
        if i==1
            x0.theta_x = opt(1);
            x0.theta_y = opt(2) ;
            x0.theta_z = opt(3);
            x0.T = [0 0 0];
        elseif i==2
            x0.theta_x = 82.0122;
            x0.theta_y = -0.0192 ;
            x0.theta_z = 87.7953;
            x0.T = [0.0228
                   -0.2070
                   -0.0783];
        else
            x0.theta_x = opt(1) + noise(1) * 10;
            x0.theta_y = opt(2) + noise(2) * 10;
            x0.theta_z = opt(3) + noise(3) * 10;
            x0.T = [0 0 0] + noise(4:6) * 0.1;
        end

    %             options = optimoptions('fmincon', 'MaxIter',5e2,'Display','iter', 'TolX', 1e-12, 'FunctionTolerance', 1e-8, 'MaxFunctionEvaluations', 3e4);
        options = optimoptions('fmincon', 'MaxIter',5e2, 'TolX', 1e-12, 'Display','off', 'FunctionTolerance', 1e-8, 'MaxFunctionEvaluations', 3e4);
        % options.UseParallel = true;
        % showproblem(prob);
        % [sol, fval, exitflag, output] = solve(prob, x0, 'Options',options);
        [sol, fval, ~, ~] = solve(prob, x0, 'Options', options);
        R_final = rotx(sol.theta_x) * roty(sol.theta_y) * rotz(sol.theta_z);
        H_LC = eye(4);
        H_LC(1:3, 1:3) = R_final;
        H_LC(1:3, 4) = sol.T';

    %             disp('H_LC: ')
    %             disp('     R:')
    %             disp(H_LC(1:3, 1:3))
    %             disp('     T:')
    %             disp(-inv(H_LC(1:3, 1:3))*H_LC(1:3, 4))
        P = intrinsic * [eye(3) zeros(3,1)] * H_LC;
        results(i).total_cost = fval;
        results(i).RMSE = sqrt(fval/size(Y,2));
        results(i).H = H_LC;
        results(i).P = P;
        if display
            disp('-- H_LC: ')
            disp('------- R:')
            disp(H_LC(1:3, 1:3))
            disp('------- T:')
            disp(-inv(H_LC(1:3, 1:3))*H_LC(1:3, 4))
            disp('-- cost:')
            disp(results(i).total_cost)
            disp('-- RMSE:')
            disp(results(i).RMSE)
        end
    end

    [min_cost, k] = min([results(:).total_cost]);
    [max_cost, ~] = max([results(:).total_cost]);
    [min_RMSE, ~] = min([results(:).RMSE]);
    final_result.RMSE = min_RMSE;
    final_result.std = std([results(:).total_cost]);
    final_result.min_cost = min_cost;
    final_result.max_cost = max_cost;
    final_H = results(k).H;
    final_P = results(k).P;
    final_result.H = final_H;
    final_result.P = final_P;
    if display
%         disp("std:")
%         disp(final_result.std)
%         disp("minimum cost:")
%         disp(min_cost)
%         disp('H_LC: ')
%         disp(' R:')
%         disp(results(k).H(1:3, 1:3))
%         disp(' RPY (XYZ):')
%         disp(rad2deg(rotm2eul(results(k).H(1:3, 1:3), "XYZ")))
%         disp(' T:')
%         disp(-inv(results(k).H(1:3, 1:3))*results(k).H(1:3, 4))
%         disp("max cost:")
%         disp(max_cost)
%         disp('H_LC: ')
%         disp(' R:')
%         disp(results(j).H(1:3, 1:3))
%         disp(' RPY (XYZ):')
%         disp(rad2deg(rotm2eul(results(j).H(1:3, 1:3), "XYZ")))
%         disp(' T:')
%         disp(-inv(results(j).H(1:3, 1:3))*results(j).H(1:3, 4))
    end
end

author: UMich-BipedLab


text:
% forcerow.m
% function to force a vector to be a single row
%
% Brian Birge
% Rev 2.0 - vectorized
% 7/1/98

function[out]=forcerow(in)
len=prod(size(in));
out=reshape(in,[1,len]);

author: HuangCongQing


text: %% S2AxisFieldTri

author: mtex-toolbox

repo: mtex

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