text:
function varargout = process_test_parametric2p( varargin )
% PROCESS_TEST_PARAMETRIC2P: Parametric two-sample tests (dependent=paired).

% @=============================================================================
% 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, Dimitrios Pantazis, 2008-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Parametric test: Paired';
    sProcess.Category    = 'Stat2';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 102;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',  'results',  'timefreq',  'matrix'};
    sProcess.OutputTypes = {'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 2;
    sProcess.isPaired    = 1;
    sProcess.isSeparator = 1;

    % === GENERIC EXTRACT OPTIONS
    % Label
    sProcess.options.extract_title.Comment    = 'Select data to test:';
    sProcess.options.extract_title.Type       = 'label';
    % Options
    sProcess = process_extract_values('DefineExtractOptions', sProcess);
    % DISABLE ABSOLUTE VALUE
    sProcess.options.isabs.Value = 0;
    sProcess.options.isnorm.Value = 0;
    sProcess.options.isabs.Hidden = 1;
    sProcess.options.isnorm.Hidden = 1;

    % === OUTPUT COMMENT
    sProcess.options.Comment.Comment = 'Comment (empty=default): ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';

    % === TEST: title
    sProcess.options.test_title.Comment    = '
Test statistic:';
    sProcess.options.test_title.Type       = 'label';
    % === TEST: type
    sProcess.options.test_type.Comment = {['Paired Student''s t-test        (A-B)~N(m,v)
' ...
                                           't = mean(A-B) / std(A-B) * sqrt(n)      df=n-1']; ...
                                          'ttest_paired'};
    sProcess.options.test_type.Type    = 'radio_label';
    sProcess.options.test_type.Value   = 'ttest_paired';
    % === TAIL FOR THE TEST STATISTIC
    sProcess.options.tail.Comment  = {'One-tailed (-)', 'Two-tailed', 'One-tailed (+)', ''; ...
                                      'one-', 'two', 'one+', ''};
    sProcess.options.tail.Type     = 'radio_linelabel';
    sProcess.options.tail.Value    = 'two';
end

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

%% ===== RUN =====
function sOutput = Run(sProcess, sInputsA, sInputsB) %#ok<DEFNU>
    sOutput = process_test_parametric2('Run', sProcess, sInputsA, sInputsB);
end

author: brainstorm-tools


text:
function studyreg = readLeksellStudyregFile(filename)
%"readLeksellStudyregFile"
%   Reads a Leksell studyreg file using decodeLeksellData and places the
%   fields into a datastructure with meaningful names.  This structure
%   shows the meaning of the transformation matrices in the Leksell Study
%   files so that it is clear which modality the matrices are going to and
%   from.  However, this information is not used when importing to CERR
%   because the target modality is always LGP (the Leksell coordinate
%   system).
%
%JRA 6/13/05
%
%LM: KRK, 05/29/07, added additional documentation
%
%Usage:
%   studyreg = readLeksellStudyregFile(filename)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

fid = fopen(filename, 'r', 'b');

data = decodeLeksellData(fid);

fclose(fid);

%If data does not have at least two elements it is effectively empty.
if length(data) < 2
    studyreg = [];
    return;
end

for i=1:length(data) - 1
    rawStudyreg = data{i};

	studyreg(i).original_modality       = rawStudyreg{1}; %MR/CT/etc
	studyreg(i).target_modality         = rawStudyreg{2}; %Leksell
  	studyreg(i).registration_value      = rawStudyreg{3}; %original modalities reg. value

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [g, gdata, gprior] = mlpgrad(net, x, t)
%MLPGRAD Evaluate gradient of error function for 2-layer network.
%
%	Description
%	G = MLPGRAD(NET, X, T) takes a network data structure NET  together
%	with a matrix X of input vectors and a matrix T of target vectors,
%	and evaluates the gradient G of the error function with respect to
%	the network weights. The error funcion corresponds to the choice of
%	output unit activation function. Each row of X corresponds to one
%	input vector and each row of T corresponds to one target vector.
%
%	[G, GDATA, GPRIOR] = MLPGRAD(NET, X, T) also returns separately  the
%	data and prior contributions to the gradient. In the case of multiple
%	groups in the prior, GPRIOR is a matrix with a row for each group and
%	a column for each weight parameter.
%
%	See also
%	MLP, MLPPAK, MLPUNPAK, MLPFWD, MLPERR, MLPBKP
%

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

% Check arguments for consistency
errstring = consist(net, 'mlp', x, t);
if ~isempty(errstring);
  error(errstring);
end
[y, z] = mlpfwd(net, x);
delout = y - t;

gdata = mlpbkp(net, x, z, delout);

[g, gdata, gprior] = gbayes(net, gdata);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function K = K_all2K(K_all)
%% Transform between K and K_all

M = size(K_all,2);
K = zeros(1,M);
for j = 1:M
    K(j) = max(K_all(:,j));
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function [Jacobian, Joint] = ARM_2DOF(L1, L2, theta1, theta2, draw, orgX, orgY)
%% Function Configuration
if (nargin == 5)
    orgX = 0;
    orgY = 0;
end

%% Modelling of the Manipulator
X(1)= orgX+ L1*cosd(theta1);
Y(1)= orgY+ L1*sind(theta1);
X(2)= X(1)+ L2*cosd(theta1 + theta2);
Y(2)= Y(1)+ L2*sind(theta1 + theta2);

% Storing the Joint Locations of the Manipulator
Joint = [orgX, orgY; X(1), Y(1); X(2), Y(2)];

% Storing the Jacobian matrix for the current Manipulator configuration
Jacobian = [-L1*sind(theta1) - L2*sind(theta1+theta2), -L2*sind(theta1+theta2);
    L1*cosd(theta1)+ L2*cosd(theta1+theta2), L2*cosd(theta1+theta2)];

% If it was mentioned in the parameter 'draw' then draw the manipulator.
if draw == true
    plot([orgX,X(1)],[orgY,Y(1)]);
    hold on
    plot([X(1),X(2)],[Y(1),Y(2)]);
    axis([-(L1 + L2), (L1 + L2), -(L1 + L2), (L1 + L2)]);
    hold off;
    title('2 DOF Planar Manipulator');
    legend('1st Link of Arm', '2nd Link of Arm');
    drawnow;
end

end

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


text:
% pop_rmbase() - remove channel baseline means from an epoched or
%                continuous EEG dataset. Calls rmbase().
% Usage:
%   >> OUTEEG = pop_rmbase( EEG ); % pop up an interactive arg entry window
%   >> OUTEEG = pop_rmbase( EEG, timerange, pointrange); % call rmbase()
%
% Graphic interface:
%    "Baseline latency range" - [edit box] Latency range for the baseline in ms.
%                               Collects the 'timerange' command line input.
%                               Empty or [] input -> Use whole epoch as baseline
%    "Baseline points vector" - [edit box] Collects the 'pointrange' command line
%                               option (below). (Overwritten by 'timerange' above).
%                               Empty or [] input -> Use whole epoch as baseline
% Inputs:
%   EEG        - Input dataset
%   timerange  - [min_ms max_ms] Baseline latency range in milliseconds.
%                                Empty or [] input -> Use whole epoch as baseline
%   pointrange - [min:max]       Baseline points vector (overwritten by timerange).
%                                Empty or [] input -> Use whole epoch as baseline
% Outputs:
%   OUTEEG     - Output dataset
%
% Note: If dataset is continuous, channel means are removed separately
%       for each continuous data region, respecting 'boundary' events
%       marking boundaries of excised or concatenated data portions.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: rmbase(), eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, 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

% 01-25-02 reformated help & license -ad

function [EEG, com] = pop_rmbase( EEG, timerange, pointrange);

com ='';
if nargin < 1
    help pop_rmbase;
    return;
end;
if isempty(EEG(1).data)
    disp('pop_rmbase(): cannot remove baseline of an empty dataset'); return;
end;
if nargin < 1
	help pop_rmbase;
	return;
end;
if nargin < 2 & EEG(1).trials > 1
	% popup window parameters
	% -----------------------
    promptstr    = {'Baseline latency range (min_ms max_ms) ([] = whole epoch):',...
         		   strvcat('Else, baseline points vector (ex:1:56) ([] = whole epoch):', ...
                           '(overwritten by latency range above).') };
	inistr     = { [num2str(EEG(1).xmin*1000) ' 0'], '' }; % latency range in ms
	result       = inputdlg2( promptstr, 'Epoch baseline removal -- pop_rmbase()', ...
                                  1,  inistr, 'pop_rmbase');
	size_result  = size( result );
	if size_result(1) == 0 return; end;

	% decode parameters
	% -----------------
	if numel(result) < 2 | ((isempty(result{1}) | strcmp(result{1},'[]') ) ...
                           & (isempty(result{2}) | strcmp(result{2},'[]')))
       timerange = [num2str(EEG(1).xmin*1000) num2str(EEG(1).xmax*1000)]; % whole epoch latency range
       fprintf('pop_rmbase(): using whole epoch as baseline.\n');

       % fprintf('pop_rmbase(): baseline limits must be specified.\n');
       % return; end;
    else
      timerange = eval( [ '[' result{1} ']' ] );
      pointrange = eval( [ '[' result{2} ']' ] );
    end
elseif nargin < 2 & EEG(1).trials == 1
	% popup window parameters
	% -----------------------
    resp = questdlg2(strvcat('Remove mean of each data channel'), 'pop_rmbase', 'Cancel', 'Ok', 'Ok');
    if strcmpi(resp, 'Cancel'), return; end;
    timerange = [];
    pointrange = [1:EEG(1).pnts];
end;

% process multiple datasets
% -------------------------
if length(EEG) > 1
    [ EEG com ] = eeg_eval( 'pop_rmbase', EEG, 'warning', 'on', 'params', ...
                            { timerange pointrange } );
    return;
end;

if exist('pointrange') ~= 1 && ~isempty(timerange)
    if (timerange(1) < EEG.xmin*1000) & (timerange(2) > EEG.xmax*1000)
        error('pop_rmbase(): Bad time range');
    end;
    pointrange = round((timerange(1)/1000-EEG.xmin)*EEG.srate+1):round((timerange(2)/1000-EEG.xmin)*EEG.srate);
end;

if isempty(timerange)
    timerange = [ EEG(1).xmin*1000 EEG(1).xmax*1000];
end;

if exist('pointrange') ~= 1 || isempty(pointrange)
    if ~isempty(timerange) && (timerange(1) < EEG.xmin*1000) & (timerange(2) > EEG.xmax*1000)
        error('pop_rmbase(): Bad time range');
    end;
    pointrange = round((timerange(1)/1000-EEG.xmin)*EEG.srate+1):round((timerange(2)/1000-EEG.xmin)*EEG.srate);
end;

if ~isempty(pointrange) && ((min(pointrange) < 1) || (max( pointrange ) > EEG.pnts))
   error('pop_rmbase(): Wrong point range');
end;

fprintf('pop_rmbase(): Removing baseline...\n');
%
% Respect excised data boundaries if continuous data
% ---------------------------------------------------
if EEG.trials == 1 & ~isempty(EEG.event) ...
                     & isfield(EEG.event, 'type') ...
                        & isstr(EEG.event(1).type)
    tmpevent = EEG.event;
	boundaries = strmatch('boundary', {tmpevent.type});
	if ~isempty(boundaries) % this is crashing
        fprintf('Pop_rmbase(): finding continuous data discontinuities\n');
        boundaries = round([ tmpevent(boundaries).latency ] -0.5-pointrange(1)+1);
        boundaries(find(boundaries>=pointrange(end)-pointrange(1))) = [];
        boundaries(find(boundaries<1)) = [];
        boundaries = [0 boundaries pointrange(end)-pointrange(1)];
        for index=1:length(boundaries)-1
            tmprange = [boundaries(index)+1:boundaries(index+1)];
            EEG.data(:,tmprange) = rmbase( EEG.data(:,tmprange), length(tmprange), ...
                                               [1:length(tmprange)]);
        end;
    else
        EEG.data = rmbase( EEG.data(:,:), EEG.pnts, pointrange );
    end;
else
    EEG.data = rmbase( EEG.data(:,:), EEG.pnts, pointrange );
end;

EEG.data = reshape( EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
EEG.icaact = [];

if ~isempty(timerange)
	com = sprintf('%s = pop_rmbase( %s, [%s]);', inputname(1), inputname(1), ...
			num2str(timerange));
else
	com = sprintf('%s = pop_rmbase( %s, [], %s);', inputname(1), inputname(1), ...
			vararg2str({pointrange}));
end;
return;

author: PatternRecognition

repo: OpenBMI


text:
function A = connect_by_inherit_neigh(pts, spls, corresp, neigh)
% build a connection matrix of downsamples (spls) by inherit neighbors of samples (pts) they
% correspond.
% There are two kinds of ways to use the neighbors, seeing the code.
%
% pts: samples
% spls: downsamples
% corresp: correspondence between pts and spls, array of |pts|*1
% neigh: neighbors of pts, a |pts|*? cell, the neighbors are sorted in ascending order by distance
% A: connection matrix of downsamples (spls)
%
% @author: JJCAO
% Changed from Andrea Tagliasacchi's code
% @data:   2010-5-7
% @version 2.0

%% visual debug conditions
SHOW_CONNECT_PROGRESS = false;
SHOW_RESULTS = false;

if SHOW_CONNECT_PROGRESS || SHOW_RESULTS
    close all;
    figure(1); movegui('northwest');set(gcf,'color','white');hold on;
    plot3( pts(:,1), pts(:,2), pts(:,3), '.r', 'markersize', 1);
    plot3( spls(:,1), spls(:,2), spls(:,3), '.g', 'markersize', 20);
    axis off; axis equal;set(gcf,'Renderer','OpenGL');
end
if ~iscell(neigh)
    tmp = cell(size(neigh,1),1);
    for i = 1:size(neigh,1)
        tmp{i}=neigh(i,:);
    end
    neigh = tmp;  clear tmp;
end
%pause(5);
%change to use the
%%
A = zeros( length(spls), length(spls) );
for pIdx=1:length(pts)
    ns =neigh{pIdx};
    pc = corresp(pIdx);
    if pc == 0
        warning('some points have no correspondence');continue;
    end

    for nIdx=1:length(ns)
        nc = corresp(ns(nIdx));
        if nc == 0
            warning('some points have no correspondence');continue;
        end
        if nc~=pc
            A(pc,nc) = A(pc,nc) + 1;
            A(nc,pc) = A(nc,pc) + 1;
            if SHOW_CONNECT_PROGRESS
                figure(1); idx = [pc, nc];
                line( spls(idx,1),spls(idx,2),spls(idx,3), 'LineWidth', 2,'Color','b');
            end
            break;%new method, generates less edges
%           %continue;old method, generates more edges.
        end
    end
end

%% if there are isolate points, connect it with its nearest neighbors.
isopts = zeros(1,0);
for i=1:size(A,1)
    A(i,i) = 1;
    if length( find( A(i,:)>0) ) == 1
        isopts(1, end+1) = i;
        warning('there are isolate points: %d', i);
    end
end
if isempty(isopts)
    spls_kdtree = kdtree_build( spls );
    for i=isopts
        neighs = kdtree_k_nearest_neighbors( spls_kdtree, spls(i,:), 2)';
        for j = neighs
            A(i,j) = 1;
            A(j,i) = 1;
        end
    end
    kdtree_delete( spls_kdtree );
end

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
%kcpRGBNTSCtoUVWstar 'converts an image from RGB NTSC to U*V*W* '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cpRGBNTSCtoUVWstar.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kcpRGBNTSCtoUVWstar(i, {'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cpRGBNTSCtoUVWstar - converts an image from RGB NTSC to U*V*W*
%
%  DESCRIPTION
% call ckfromRGBNTSC with the parameter -labstar to convert an image from RGB NTSC to U*V*W*
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
% All references for the COLOR toolbox are on the file $COLOR/Readme.html
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
import matlab.unittest.plugins.DiagnosticsRecordingPlugin
import matlab.unittest.plugins.TestReportPlugin

runner = matlab.unittest.TestRunner.withNoPlugins;
runner.addPlugin(DiagnosticsRecordingPlugin);
runner.addPlugin(TestReportPlugin.producingPDF('Report.pdf',...
    'IncludingPassingDiagnostics',true,'IncludingCommandWindowText',true));
runner.ArtifactsRootFolder = pwd;

suite = testsuite('HardwarePerformanceTests');

if ~exist('logs', 'dir')
    mkdir('logs')
end

results = runner.run(suite);
t = table(results);
disp(t);

exit(any([results.Failed]));

author: analogdevicesinc


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z).*(1-sigmoid(z));

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

end

repo: MachineLearningEx


text:
function writePLY(filename,V,F,mode)
  % WRITEPLY wrapper for write_ply
  %
  % Inputs:
  %   filename  path to output .ply file
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of mesh triangle indices
  %   mode  followed by
  %      'ascii'                ASCII text data
  %      {'binary_little_endian'}   binary data, little endian
  %      'binary_big_endian'      binary data, big endian
  %
  %

  if nargin < 4
    mode = 'binary_little_endian';
  end

  fid = fopen(filename,'wt');
  fprintf(fid,'ply\nformat %s 1.0\n',mode);
  fprintf(fid,'element vertex %d\n',size(V,1));
  fprintf(fid,'property double x\n');
  fprintf(fid,'property double y\n');
  fprintf(fid,'property double z\n');
  fprintf(fid,'element face %d\n',size(F,1));
  fprintf(fid,'property list int int vertex_indices\n');
  fprintf(fid,'end_header\n');
  FF = [size(F,2)*ones(size(F,1),1) F-1];
  switch mode
  case 'ascii'
    % do nothing
    fprintf(fid,'%0.17f %0.17f %0.17f\n',V');
    format = [repmat('%d ',1,size(FF,2)) '\n'];
    fprintf(fid,format,FF');
  case {'binary_little_endian','binary_big_endian'}
    fclose(fid);
    switch mode
    case 'binary_little_endian'
      fid = fopen(filename,'a','ieee-le');
    case 'binary_big_endian'
      fid = fopen(filename,'a','ieee-be');
    end
    fwrite(fid,V','double');
    fwrite(fid,FF','int');
  otherwise
    error('Unsupported format: %s',mode);
  end
  fclose(fid);

  %% OLD 200x slower way (possible to produce slightly smaller files and control
  %% over precision)
  %write_ply(V,F,filename,mode);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function T = CreateDatabase(TrainDatabasePath)
% Align a set of face images (the training set T1, T2, ... , TM )
%
% Description: This function reshapes all 2D images of the training database
% into 1D column vectors. Then, it puts these 1D column vectors in a row to
% construct 2D matrix 'T'. Each column of 'T' is a training image, which has been reshaped into a 1D vector.
% Also, P is the total number of MxN training images and C is the number of
% classes.
%
%
% Argument:     TrainDatabasePath      - Path of the training database
%
% Returns:      T                      - A 2D matrix, containing all 1D image vectors.
%                                        The length of 1D column vectors is MN and 'T' will be a MNxP 2D matrix.
%
% See also: STRCMP, STRCAT, RESHAPE

% Original version by Amir Hossein Omidvarnia, October 2007
%                     Email: aomidvar@ece.ut.ac.ir

%%%%%%%%%%%%%%%%%%%%%%%% File management
TrainFiles = dir(TrainDatabasePath);
Train_Number = 0;

for i = 1:size(TrainFiles,1)
    if not(strcmp(TrainFiles(i).name,'.')|strcmp(TrainFiles(i).name,'..')|strcmp(TrainFiles(i).name,'Thumbs.db'))
        Train_Number = Train_Number + 1; % Number of all images in the training database
    end
end

%%%%%%%%%%%%%%%%%%%%%%%% Construction of 2D matrix from 1D image vectors
T = [];
for i = 1 : Train_Number

    % I have chosen the name of each image in databases as a corresponding
    % number. However, it is not mandatory!
    str = int2str(i);
    str = strcat('\',str,'.jpg');
    str = strcat(TrainDatabasePath,str);

    img = imread(str);
    img = rgb2gray(img);

    [irow icol] = size(img);

    temp = reshape(img',irow*icol,1);   % Reshaping 2D images into 1D image vectors
    T = [T temp]; % 'T' grows after each turn
end

T = double(T);

author: Sable

repo: mcbench-benchmarks


text:
function deformation_field = PTKSolveMatchedImagesForFluidRegistration(image_to_transform, reference_image, reporting)
    % PTKSolveMatchedImagesForFluidRegistration.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    MimImageUtilities.MatchSizes(reference_image, image_to_transform);

    if ~isequal(image_to_transform.VoxelSize, reference_image.VoxelSize)
        reporting.Error('SolveMatchedImagesForFluidRegistration:UnequalVoxelSize', 'SolveMatchedImagesForFluidRegistration requires images to have the same voxel size');
    end

    % Convert to distance transforms
    dt_float = PTKRunForEachComponentAndCombine(@MimImageUtilities.GetNormalisedDT, image_to_transform, image_to_transform, reporting);
    dt_ref = PTKRunForEachComponentAndCombine(@MimImageUtilities.GetNormalisedDT, reference_image, reference_image, reporting);

    % Solve to find the deformation field between these images
    [deformation_field, ~] = Solve(dt_float, dt_ref, reporting);
end

function [deformation_field, deformed_image] = Solve(image_1, image_2, reporting)
    deformation_field = image_2.BlankCopy;
    deformed_image = image_2.BlankCopy;
    voxel_size = image_2.VoxelSize;
    if ~isequal(voxel_size, image_1.VoxelSize)
        reporting.Error('PTKGetFluidDeformationField:DifferentVoxelSizes', 'Images must have the same voxel size');
    end
    if ~isequal(image_1.ImageSize, image_2.ImageSize)
        reporting.Error('PTKGetFluidDeformationField:DifferentVoxelSizes', 'Images must have the same image size');
    end

    if PTKSoftwareInfo.ShowRegistrationConvergence
        display_value = 'iter';
    else
        display_value = 'off';
    end
    iterations = 200;
    % X corresponds to image columns and Y image rows
    options = npRegSet(...
        'Display', display_value, ...
        'Regularizer', 'curvature', ... % elastic | {fluid} | diffusion | curvature
        'BoundaryCond', 'Neumann', ...
        'VoxSizeX', voxel_size(2),'VoxSizeY', voxel_size(1), 'VoxSizeZ', voxel_size(3), ...
        'maxiter', iterations, ...
        'RegularizerFactor', 1,...
        'BodyForceTol', PTKSoftwareInfo.RegistrationBodyForceTol, ... % Termination tolerance on the body force [1e-2]
        'BodyForceDiffTol', PTKSoftwareInfo.RegistrationBodyForceDiffTol, ... % Termination tolerance on the difference between successive body force estimates [1e-2]..
        'SimMeasPercentDiffTol', 0.001, ... % Termination tolerance on the percentage difference between successive similarity measure estimates [ positive scalar {1e-2} ]
        'FixedPointMaxFlowDistance', 5 ...
);

%         'SimilarityMeasure', 'SSD', ... % [ {SSD} | NCC | CR | MI | NMI ]
%UDiffTol - Fixed point iteration termination tolerance on the maximum sum
%       of squared differences between successive displacement field
%       estimates [ positive scalar {1e-2} ]
%UDiffTol - Fixed point iteration termination tolerance on the maximum sum
%       of squared differences between successive velocity field estimates
%           [ positive scalar {1e-2} ]
%SimMeasTol - Fixed point iteration termination tolerance on the similarity
%       measure [ positive scalar {1e-2} ]
%SimMeasDiffTol - Fixed point iteration termination tolerance on the
%       difference between successive similarity measure estimates
%           [ positive scalar {1e-2} ]
%SimMeasPercentDiffTol - Fixed point iteration termination tolerance on
%       the percentage difference between successive similarity measure
%       estimates [ positive scalar {1e-2} ]
%FixedPointMaxFlowDistance - Maximum distance (in voxels) that the
%       deformation field is allowed to flow during each fixed point
%       iteration [ positive scalar {5.0} ]
%RegridTol - Tolerance on the Jacobian of the deformation field, below which
%       regridding will take place [ positive scalar {0.0025} ]
%Mu - Lame constant for use with elastic or fluid regularizers
%           [ positive scalar {1.0} ]
%Lambda - Lame constant for use with elastic or fluid regularizers
%           [ nonnegative scalar {0.0} ]
    [deformed_image_raw, deformation_field_raw, exit_flag] = npReg(image_2.RawImage, image_1.RawImage, options);
    deformed_image.ChangeRawImage(deformed_image_raw);

    if (exit_flag ~= 1)
        reporting.ShowWarning('PTKGetFluidDeformationField:RegistrationDidNotConverge', 'npReg registration did not converge.', []);
    end

    translation = (image_2.Origin - image_1.Origin).*image_2.VoxelSize;
    deformation_field_raw(:,:,:,1) = deformation_field_raw(:,:,:,1) + translation(1);
    deformation_field_raw(:,:,:,2) = deformation_field_raw(:,:,:,2) + translation(2);
    deformation_field_raw(:,:,:,3) = deformation_field_raw(:,:,:,3) + translation(3);

    deformation_field.ChangeRawImage(deformation_field_raw);
end

author: tomdoel

repo: pulmonarytoolkit


text:
%
% Replaces model parameters (from a vector) into an object.
% Used for optimizing a parse, with fixed type-level parameters,
% to a new image
%
% Replace variables:
%  image blur, pixel noise, and affine parameters
%  stroke position, token scale, and token shape
%
function vec_to_model_fit_token(theta,M,list_sid)

    if ~exist('list_sid','var')
       list_sid = 1:M.ns;
    end

    count = 1;

    % pixel noise
    M.epsilon = theta(count);
    count = count + 1;

    % image blur
    M.blur_sigma = theta(count);
    count = count + 1;

    % affine warp
    x = count:count+3;
    M.A = vec(theta(x));
    count = count + 4;

    for i=list_sid

        % stroke position
        x = count:count+1;
        M.S{i}.pos_token = vec(theta(x))';
        count = count + 2;

        % sub-stroke shapes (token)
        nn = numel(M.S{i}.shapes_token);
        x = count:count+nn-1;
        M.S{i}.shapes_token = reshape(theta(x),size(M.S{i}.shapes_token));
        count = count + nn;

        % sub-stroke scales (token)
        nn = numel(M.S{i}.invscales_token);
        x = count:count+nn-1;
        M.S{i}.invscales_token = reshape(theta(x),size(M.S{i}.invscales_token));
        count = count + nn;

        % eval_spot_token
        if ~isempty(M.S{i}.R) && strcmp(M.S{i}.R.type,'mid')
            M.S{i}.R.eval_spot_token = theta(count);
            count = count + 1;
        end
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% Example illustrating the performance of different preconditioners for
% the Gauss-Newton system in LDDMM. The example is described in detail in
% Section 4.2 of the paper:
%
% @article{MangRuthotto2017,
%   Title = {A {L}agrangian {G}auss--{N}ewton--{K}rylov solver for mass- and intensity-preserving diffeomorphic image registration},
%   Year = {2017},
%   Journal = {SIAM Journal on Scientific Computing},
%   Author = {A. Mang, L. Ruthotto},
% }
%
% =========================================================================
clear all; clc;
setup2Ddisc2CData

%% run affine pre-registration
imgModel('reset','imgModel','splineInterMex','regularizer','moments','theta',.1);
        trafo('set','trafo','affine2D');
        distance('set','distance','SSD');

alpha = [6e2 0];
parametric = 0;
pad  = 0.5;
nt = 0;
N    = 3;
mV     = @(m) ceil(1*m);
minLevel = 5;
maxLevel = 5;

%% run multilevel LDDMM registration

% 1) setup grid for velocities (padded)
omegaV = omega; omegaV(1:2:end) = omegaV(1:2:end)-pad;  omegaV(2:2:end) = omega(2:2:end)+pad;

% 2) setup regularizer
regularizer('reset','regularizer','mfDiffusionCC','alpha',alpha,'nt',nt,'HessianShift',1e-2); % stationary velocity

NPIRpara    = optPara('NPIR-GN');
NPIRpara.maxIter = 40;
NPIRpara.scheme = @GaussNewtonLDDMM;
[vc,~,wc,his] = MLLDDMM(ML,'minLevel',minLevel,'maxLevel',maxLevel,...
    'omegaV',omegaV,'mV',mV,'N',N,'parametric',parametric,'NPIRpara',NPIRpara,'plots',1);

%% generate plots to explore conditioning / sparsity
xc = getNodalGrid(omega,ML{minLevel}.m);
[T,R] = imgModel('coefficients',ML{minLevel}.T,ML{minLevel}.R,omega);
Rc    = imgModel(R,omega,center(xc,ML{minLevel}.m));

mV = ML{minLevel}.m;
fctn = @(vc) LDDMMobjFctn(T,Rc,omega,ML{minLevel}.m,0*vc,center(xc,ML{minLevel}.m),omegaV,mV,N,vc);

close all
regularizer('set','regularizer','mbDiffusionCC')
[Sc,dS,d2Smb] = regularizer(vc,omegaV,mV);
% get objective function
[J0,p0,dJ0,Hmb0] = fctn(0*vc);
[Jc,pc,dJc,Hmbc] = fctn(vc);

regularizer('set','regularizer','mfDiffusionCC')
[J0,p0,dJ0,Hmf0] = fctn(0*vc);
[Jc,pc,dJc,Hmfc] = fctn(vc);

%%
fig = figure(); clf;
fig.Name =sprintf('%s: Sparsity of Hessian',mfilename);

subplot(1,2,1)
spy(Hmb0);
title('H(v0)');

subplot(1,2,2);
spy(Hmbc);
title('H(vOpt');
%% explore PCG convergence at first iteration
hd = prod((omega(2:2:end)-omega(1:2:end))./ML{minLevel}.m);
Hmb0 = Hmb0 + hd* regularizer('get','HessianShift')* speye(size(Hmb0));
D   = diag(Hmb0); % D is diagonal
L   = tril(Hmb0); % Symmetric Gauss Seidel Preconditioning,

PCjac = @(x) D.\x;
PCsgs = @(x) full(L\(D.*(L'\x)));
PCnone = @(x) x;

[~,iter0cg,relres0cg,resvec0cg] = spectralPrecondPCG(-dJ0(:), Hmf0, 250, 1e-10,'prec',PCnone);
[~,iter0jac,relres0jac,resvec0jac] = spectralPrecondPCG(-dJ0(:), Hmf0, 250, 1e-10,'prec',PCjac);
[~,iter0sgs,relres0sgs,resvec0sgs] = spectralPrecondPCG(-dJ0(:), Hmf0, 250, 1e-10,'prec',PCsgs);
[~,iter0sp,relres0sp,resvec0sp] = spectralPrecondPCG(-dJ0(:), Hmf0, 250, 1e-10);

%%
figConv = figure(); clf
figConv.Name = 'PCG convergence';

subplot(1,2,1);
semilogy(resvec0cg/resvec0cg(1));
hold on;
semilogy(resvec0jac/resvec0jac(1));
semilogy(resvec0sgs/resvec0sgs(1));
semilogy(resvec0sp/resvec0sp(1));
title('first GN iteration');
ax = [0 250 1e-10 10];
axis(ax)
legend('CG','PCG-Jac','PCG-SGS','PCG-Spec','Location','SouthWest')

%% explore PCG convergence at final iteration
hd = prod((omega(2:2:end)-omega(1:2:end))./ML{minLevel}.m);
Hmb0 = Hmb0 + hd* regularizer('get','HessianShift')* speye(size(Hmb0));
D   = diag(Hmbc); % D is diagonal
L   = tril(Hmbc); % Symmetric Gauss Seidel Preconditioning,

linSol = @spectralPrecondPCG;
PCjac = @(x) D.\x;
PCsgs = @(x) full(L\(D.*(L'\x)));
PCnone = @(x) x;

[~,itercg,relrescg,resveccg]    = linSol(-dJc(:), Hmfc, 250, 1e-10,'prec',PCnone);
[~,iterjac,relresjac,resvecjac] = linSol(-dJc(:), Hmfc, 250, 1e-10,'prec',PCjac);
[~,itersgs,relressgs,resvecsgs] = linSol(-dJc(:), Hmfc, 250, 1e-10,'prec',PCsgs);
[~,itersp,relressp,resvecsp]    = linSol(-dJc(:), Hmfc, 250, 1e-10);

%%
fig = figure(figConv);
subplot(1,2,2);
semilogy(resveccg/resveccg(1));
hold on;
semilogy(resvecjac/resvecjac(1));
semilogy(resvecsgs/resvecsgs(1));
semilogy(resvecsp/resvecsp(1));
title('final iteration');

legend('CG','PCG-Jac','PCG-SGS','PCG-Spec','Location','SouthWest')
axis(ax)


text:
 function ob = subsasgn(ob, sub, arg)
%function ob = subsasgn(ob, sub, arg)
% method for "ob.sub = arg" etc.
% caution: reassigning dimensions (.idim .odim .imask .omask) is dangerous!

if ~streq(sub(1).type, '.')
	fail('only "ob.field = value" is supported')
end

% trick: for struct.sub = object
% e.g., A.arg.field = A;
if isstruct(ob)
%	ob = subsasgn(ob, sub, arg); % no, causes recursion
%	ob = builtin('subsasgn', ob, sub, arg);	% samit basu dislikes this

	% recurse a.b.c.d... to avoid builtin:
	if numel(sub) > 1
		ob.(sub(1).subs) = subsasgn(ob.(sub(1).subs), sub(2:end), arg);
	else
		ob.(sub(1).subs) = arg;
	end

% trick: just convert it to a structure and then back!
else
	cl = class(ob);
	ob = struct(ob);
	name = sub.subs;

	% note: ob.arg.mask = new_mask may cause errors.
	if numel(sub) > 1 % ob.arg.name = ... or ob.meth.name =
		ob.(sub(1).subs) = subsasgn(ob.(sub(1).subs), sub(2:end), arg);

	elseif isfield(ob, name) % reassigning private object data is dangerous
		fail('too dangerous to subsasgn "%s"', name)

%	elseif isfield(ob.meth, name) % user may reassign method handles
%		ob.meth = subsasgn(ob.meth, sub, arg);

%	elseif isfield(ob.arg, name) % user may reassign arg internals
%		ob.arg = subsasgn(ob.arg, sub, arg);

	else
		fail('%s is not a field', name)
	end

	ob = class(ob, cl);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [money, accum_bet, acc, central_money, bk_money, nTrans] = calc_money_bks_H_A_opCls(aux, bet_type, money, ...
    bet, increase_bet, bookies_home_code_yrs, bookies_away_code_yrs, bk_money, accumulator)

% Game by game, place bets on each bookie
% Check if there is enough money in a bookie and place a bet.
% Keep the history of the central money to see the zeros

% lim =  money + bet * 12;

loss = 0.03; %loss on every transaction

central_money = money;

z = 1;
c = 1;
nTrans = 0;

bet_type = bet_type(~isnan(bet_type));
aux = aux(~isnan(bet_type), :);
bookies_home_code_yrs = bookies_home_code_yrs(~isnan(bet_type));
bookies_away_code_yrs = bookies_away_code_yrs(~isnan(bet_type));

for gm = 1 : size(aux, 1)

    if bet > 3000
        bet = 3000;
    end

    accum_bet(z) = bet; %#ok

    if bet_type(gm) == 1

        id_bk = bookies_home_code_yrs(gm);
        % Add money to the bookie if we are below the limit
        if bk_money(id_bk) < bet
            bk_money(id_bk) = bk_money(id_bk) + bet;
            central_money(c+1) = central_money(c) - bet;
            % Transfer money from another bookie if we ran out of central
            % money
            if central_money(c+1) <= bet
               [~, id_max] = max(bk_money);
               bk_money(id_max) = bk_money(id_max) - (bet + (bet * loss)); % we are losing a margin with the movement
               central_money(c+1) = central_money(c+1) + bet;
               nTrans = nTrans + 1;
            end
            c = c + 1;
        end

        if aux(gm,1) == 1
            mony = (aux(gm,5) - 1) * bet; % bet won
            bk_money(id_bk) = bk_money(id_bk) + mony;
            acc(z) = 1;
        else
            mony = -1 * bet; % bet lost
            bk_money(id_bk) = bk_money(id_bk) + mony;
            acc(z) = 0;
        end

    elseif bet_type(gm) == 2

        id_bk = bookies_away_code_yrs(gm);
        % Add money to the bookie if we are below the limit
        if bk_money(id_bk) < bet
            bk_money(id_bk) = bk_money(id_bk) + bet;
            central_money(c+1) = central_money(c) - bet;

            % Transfer money from another bookie if we ran out of central
            % money
            if central_money(c+1) <= bet
               [~, id_max] = max(bk_money);
               bk_money(id_max) = bk_money(id_max) - (bet + (bet * loss)); % we are losing a margin with the movement
               central_money(c+1) = central_money(c+1) + bet;
               nTrans = nTrans + 1;
            end

            c = c + 1;
        end

        if aux(gm,1) == 2
            mony = (aux(gm,7) - 1) * bet; % bet won
            bk_money(id_bk) = bk_money(id_bk) + mony;
            acc(z) = 1;
        else
            mony = -1 * bet; % bet lost
            bk_money(id_bk) = bk_money(id_bk) + mony;
            acc(z) = 0;
        end

    elseif bet_type(gm) == 3

        if aux(gm,1) == 0
            mony = (aux(gm,6) - 1) * bet; % bet won
            acc(z) = 1;
        else
            mony = -1 * bet; % bet lost
            acc(z) = 0;
        end

    else
        continue; % not a game to bet, go to next game

    end

    money(z+1) = money(z) + mony;

    bet = money(z+1) * accumulator;
%     if money(z+1) > lim && increase_bet == 1
%
%         bet = bet * 2;
%
%         lim = lim * 2;
%
%     end
    z = z + 1;

%     if z == 4000
%         disp('here')
%     end
end

repo: BeatTheBookie


text:
function haveNewerRelease = tapas_check_for_new_release(verbose)
    %% Check if a new release is available.
    %
    % Input
    %       verbose         -- Show information on command prompt:
    %                               0   no information
    %                               1   only if there is a new release (or not)
    %                               2   as 1, but if there is a new release also
    %                                   the release notes for the lase release
    %                               3   as 1, but if there is a new release also
    %                                   the release notes for all newer
    %                                   releases.
    %
    % Output
    %       haveNewerRelease  -- True if there is a newer release on github

    % muellmat@ethz.ch
    % copyright (C) 2020
    %

    if nargin < 1
        verbose = 3;
    end

    [online_version,data] = tapas_check_online_version();
    if isempty(online_version) % Could not communicate with server
        haveNewerRelease = false;
        if verbose
            fprintf(1, 'Could not reach github to get current TAPAS version!\n');
        end
        return;
    end
    offline_version = tapas_get_current_version();
    % both versions are strings in the form '3.3.0'.
    haveNewerRelease = tapas_compare_versions(online_version,offline_version);

    if verbose % Show info
        if ~haveNewerRelease
            fprintf(1, 'Your TAPAS version is up-to-date.\n\n');
        else
            fprintf(1, ['There is a new TAPAS release available (installed %s /'...
                +'newest %s)!\n\n'],offline_version,online_version);
            if verbose > 1 % Show release notes of the new releases
                try % in case the api changed and the structure of the struct is
                    % different
                    n_data = numel(data);
                    fprintf(1, 'Release notes:\n')
                    if verbose == 2
                        n_data = 1; % Just iterate once through loop
                    end
                    for i_data = 1:n_data
                        str_version = data(i_data).tag_name;
                        str_version = strrep(str_version,'v','');
                        if tapas_compare_versions(str_version,offline_version)
                            fprintf(1, '===== TAPAS v%s =====\n',str_version)
                            body = data(i_data).body;
                            body = tapas_remove_problematic_escape_characters(body);
                            fprintf(1,body);
                            fprintf(1,'\n======================\n')
                        end
                    end
                catch
                    fprintf(1, ['Cannot print release notes. Maybe the github API'...
                        ' has changed.\n']);
                end
            end
        end
    end

end

function str = tapas_remove_problematic_escape_characters(str)
    %% Removing problematic escape characters
    %
    % Input
    %       str     -- String with problematic escape characters
    %
    % Output
    %       str     -- String without problematic escape characters

    % matlab does not like some escape sequences from github
    str = strrep(str,'\','');
end

author: translationalneuromodeling

repo: tapas


text:
function [out] = baseflow_5(p1,p2,S,Smax,dt)
%baseflow_5

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Non-linear scaled outflow from a reservoir
% Constraints:  f <= S/dt
% @(Inputs):    p1   - base outflow rate [mm/d]
%               p2   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               dt   - time step size [d]

out = min(S/dt,p1*((max(S,0)/Smax)^p2));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function indx = i4vec_indexed_heap_d ( n, a, indx )

%*****************************************************************************80
%
%% I4VEC_INDEXED_HEAP_D creates a descending heap from an indexed I4VEC.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%    An indexed I4VEC is an I4VEC of data values, and an I4VEC of N indices,
%    each referencing an entry of the data vector.
%
%    The function adjusts the index vector INDX so that, for 1 <= J <= N/2,
%    we have:
%      A(INDX(2*J))   <= A(INDX(J))
%    and
%      A(INDX(2*J+1)) <= A(INDX(J))
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the size of the index array.
%
%    Input, integer A(*), the data vector.
%
%    Input, integer INDX(N), the index array.
%    Each entry of INDX must be a valid index for the array A.
%
%    Output, integer INDX(N), the indices have been reordered into
%    a descending heap.
%

%
%  Only nodes N/2 down to 1 can be "parent" nodes.
%
  for i = floor ( n / 2 ) : - 1 : 1
%
%  Copy the value out of the parent node.
%  Position IFREE is now "open".
%
    key = indx(i);
    ifree = i;

    while ( 1 )
%
%  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
%  IFREE.  (One or both may not exist because they exceed N.)
%
      m = 2 * ifree;
%
%  Does the first position exist?
%
      if ( n < m )
        break
      end
%
%  Does the second position exist?
%
      if ( m + 1 <= n )
%
%  If both positions exist, take the larger of the two values,
%  and update M if necessary.
%
        if ( a(indx(m)) < a(indx(m+1)) )
          m = m + 1;
        end

      end
%
%  If the large descendant is larger than KEY, move it up,
%  and update IFREE, the location of the free position, and
%  consider the descendants of THIS position.
%
      if ( a(indx(m)) <= a(key) )
        break
      end

      indx(ifree) = indx(m);
      ifree = m;

    end
%
%  Once there is no more shifting to do, KEY moves into the free spot IFREE.
%
    indx(ifree) = key;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%HOUGHPOINT2LINE  Calculates coordinates of line segment corresponded by point in Hough space
%
%     line = cv.HoughPoint2Line(houghPoint, srcImgInfo)
%     line = cv.HoughPoint2Line(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __houghPoint__ Point in Hough space `[x,y]`.
% * __srcImgInfo__ The source (input) image of Hough transform.
%
% ## Output
% * __line__ Coordinates of line segment corresponded by point in Hough space,
%   a 4-element integer vector `[vx,vy, ux,uy]`.
%
% ## Options
% * __AngleRange__ The part of Hough space where point is situated. See
%   cv.FastHoughTransform, default `ARO_315_135`.
% * __MakeSkew__ Specifies to do or not to do image skewing. See
%   cv.FastHoughTransform, default 'Deskew'.
% * __Rules__ Specifies strictness of line segment calculating. This specifies
%   the degree of rules validation. This can be used, for example, to choose
%   a proper way of input arguments validation. Default 'IgnoreBorders':
%   * __Strict__ Validate each rule in a proper way.
%   * __IgnoreBorders__ Skip validations of image borders.
%
% The function calculates coordinates of line segment corresponded by point in
% Hough space.
%
% ### Notes
% - If `Rules` parameter set to 'Strict' then returned line cut along the
%   border of source image.
% - If `Rules` parameter set to 'IgnoreBorders' then in case of point, which
%   belongs the incorrect part of Hough image, returned line will not
%   intersect source image.
%
% See also: cv.FastHoughTransform, cv.HoughLines, hough, houghlines, houghpeaks
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Used for printing and displaying image
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

function [ varargout ] = drawFigFrames( varargin )
% evaluate function according to the number of inputs and outputs

    if nargout(varargin{1}) > 0
        [varargout{1:nargout(varargin{1})}] = feval(varargin{:});
    else
        feval(varargin{:});
    end
end

function [fig_info] = init_fig(frame_sz, varargin)
    addpath_export_fig();

    screen_sz = get(0,'ScreenSize');
    fig_info.fg_h = figure('units', 'pixels', 'position', ...
        [screen_sz([3 4])-frame_sz([2 1])-40 frame_sz([2 1])], ...
        'paperpositionmode', 'auto', varargin{:});
    fig_info.ax_h = axes('Parent',fig_info.fg_h);
end

function h = imshow_fig(frame, fig_info, varargin)
    h = imshow(frame, 'Parent',fig_info.ax_h, ...
        'InitialMagnification','fit', 'Border','tight', varargin{:});
%     truesize;
end

function h = imagesc_fig(frame, fig_info, varargin)
    h = imagesc(frame, 'Parent',fig_info.ax_h, varargin{:});
    adjust_ax(fig_info);
end

function h = image_fig(frame, fig_info, varargin)
    h = image(frame, 'Parent', fig_info.ax_h, varargin{:});
    adjust_ax(fig_info);
end

function holdonoff(on_off, fig_info)
    hold(fig_info.ax_h, on_off);
end

function fg_h = get_fig_h(fig_info)
    fg_h = fig_info.fg_h;
end

function ax_h = get_ax_h(fig_info)
    ax_h = fig_info.ax_h;
end

function export_fig_only(im_filepath, fig_info)
    adjust_ax(fig_info);
%     set(fig_info.ax_h, 'Units','pixels','Position',[200 200 c r]);
    export_fig(im_filepath, '-a1', fig_info.fg_h);
%     set(gca,'Position',[0 0 1 1]);  %# Modify axes size
end

function export_fig_and_clear(im_filepath, fig_info)
    export_fig_only(im_filepath, fig_info);
    delete(get(fig_info.ax_h, 'children'));
end

function adjust_ax(fig_info)
    set(fig_info.ax_h, 'Units','normalized', ...
        'position', [0 0 1 1], 'visible', 'off');
end

function close_fig(fig_info)
    delete(fig_info.ax_h);
    close(fig_info.fg_h);
end

function addpath_export_fig(extern_dir)
    % get export_fig in the path
    if exist('extern_dir', 'var') == 0
        extern_dir = 'C:\Users\Administrator\Documents\MATLAB\Rigor_untouched\extern_src';
    end
    addpath(fullfile(extern_dir, 'utils', 'export_fig'));
end

author: Cloud-CV

repo: object-proposals


text:
function Y = ttm(X,V,varargin)
%TTM Sparse tensor times matrix.
%
%   Y = TTM(X,A,N) computes the n-mode product of the sptensor X with
%   a dense matrix A; i.e., X x_N A.  The integer N specifies the
%   dimension (or mode) of X along which A should be multiplied.  If
%   size(A) = [J,I], then X must have size(X,N) = I.  The result will
%   will be a (dense) tensor or sptensor of the same order and size as X
%   except that size(Y,N) = J.
%
%   Y = TTM(X,{A,B,C,...}) computes the n-mode product of the sptensor
%   X with a sequence of matrices in the cell array.  The n-mode
%   products are computed sequentially along all dimensions (or modes)
%   of X. The cell array contains ndims(X) matrices.
%
%   Y = TTM(X,{A,B,C,...},DIMS) computes the sequence tensor-matrix
%   products along the dimensions specified by DIMS.
%
%   Y = TTM(...,'t') performs the same computations as above except
%   the matrices are transposed.
%
%   Examples
%   X = sptenrand([5 3 4 2], 10);
%   A = rand(4,5); B = rand(4,3); C = rand(3,4); D = rand(3,2);
%   Y = ttm(X, A, 1)         %<-- computes X times A in mode-1
%   Y = ttm(X, {A,B,C,D}, 1) %<-- same as above
%   Y = ttm(X, A', 1, 't')   %<-- same as above
%   Y = ttm(X, {A,B,C,D}, [1 2 3 4]) %<-- 4-way multiply
%   Y = ttm(X, {D,C,B,A}, [4 3 2 1]) %<-- same as above
%   Y = ttm(X, {A,B,C,D})            %<-- same as above
%   Y = ttm(X, {A',B',C',D'}, 't')   %<-- same as above
%   Y = ttm(X, {C,D}, [3 4])     %<-- X times C in mode-3 & D in mode-4
%   Y = ttm(X, {A,B,C,D}, [3 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, [1 2 4])   %<-- 3-way multiply
%   Y = ttm(X, {A,B,C,D}, [1 2 4]) %<-- same as above
%   Y = ttm(X, {A,B,D}, -3)        %<-- same as above
%   Y = ttm(X, {A,B,C,D}, -3)      %<-- same as above
%
%   See also SPTENSOR, TENSOR/TTM.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Check the number of arguments
if (nargin < 2)
    error('TTM requires at least two arguments.');
end

%% Create 'n' and 'tflag' arguments from varargin
n = 1:ndims(X);
tflag = '';
if numel(varargin) == 1
    if ischar(varargin{1})
        tflag = varargin{1};
    else
        n = varargin{1};
    end
elseif numel(varargin) == 2
    n = varargin{1};
    tflag = varargin{2};
end

%% Handle cell array
if iscell(V)
    % Copy n into dims
    dims = n;
    % Check that the dimensions are valid
    [dims,vidx] = tt_dimscheck(dims,ndims(X),numel(V));
    % Calculate individual products
    Y = ttm(X, V{vidx(1)}, dims(1), tflag);
    for i = 2 : numel(dims)
        Y = ttm(Y, V{vidx(i)}, dims(i), tflag);
    end
    % All done
    return;
end

%% Check the second argument
if ndims(V) ~= 2
    error('tensor/ttm: 2nd argument must be a matrix.');
end

%% Flip V is transposed
if tflag == 't'
    V = V';
end

%% Check n
if numel(n) ~= 1 || (n < 0) || (n > ndims(X))
    error('Dimension N must be between 1 and NDIMS(X).');
end

%% Compute the product

% Check that sizes match!
if size(X,n) ~= size(V,2)
    error('Size mismatch on V');
end

% Compute the new size
siz = size(X);
siz(n) = size(V,1);

% Compute Xn'
Xnt = sptenmat(X,n,'t');

% Extract the dimensions
rdims = Xnt.rdims;
cdims = Xnt.cdims;

% Convert to sparse matrix and do the multiplication; result is generally a
% dense matrix
Z = double(Xnt) * V';

if nnz(Z) <= 0.5 * prod(siz)
    % Final result is a *sparse* tensor
    Ynt = sptenmat(Z, rdims, cdims, siz);
    Y = sptensor(Ynt);
else
    % Final result is a *dense* tensor
    Ynt = tenmat(Z, rdims, cdims, siz);
    Y = tensor(Ynt);
end

end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


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

function MU_update_dim_control(handles, dimSize)

% refresh dimension tabgroup
if isfield(handles,'MDimension_tabgroup')
    tabs=get(handles.MDimension_tabgroup,'Children');
    for i=1:length(tabs)
        delete(get(tabs(i),'Children'));
    end
    delete(handles.MDimension_tabgroup);
    handles = rmfield(handles,'MDimension_tabgroup');
end

% initialize dimension tabgroup
if  numel(dimSize)> 2
    handles.MDimension_tabgroup=uitabgroup(handles.MDimension_uipanel);
    for i=3:numel(dimSize)

        handles.(['MDim' num2str(i) '_tab'])=uitab( handles.MDimension_tabgroup,'title',['Dim' num2str(i)] ,'Units','normalized');
        handles.(['Dim' num2str(i) '_edit'])=uicontrol(handles.(['MDim' num2str(i) '_tab']),'Style', 'edit','Units','normalized','BackgroundColor',[1 1 1],...
            'Position', [0.81 0 0.1 1],'TooltipString',['Matrix Dimension ' num2str(i)],'string',1,'Enable','off');
        handles.(['Dim' num2str(i) '_text'])=uicontrol(handles.(['MDim' num2str(i) '_tab']),'Style', 'text','Units','normalized','BackgroundColor',[1 1 1],...
            'Position', [0.91 0 0.08 1],'string',['/' num2str(dimSize(i))]);
        if dimSize(i) > 1
            set(handles.(['Dim' num2str(i) '_edit']),'Enable','on');
            handles.(['Dim' num2str(i) '_slider'])=uicontrol(handles.(['MDim' num2str(i) '_tab']),'Style', 'slider','Units','normalized','BackgroundColor',[1 1 1],...
                'Position', [0.01 0 0.8 1],'TooltipString',['Matrix Dimension ' num2str(i)],...
                'Value',1,'Min',1,'Max',dimSize(i),'SliderStep',[1/dimSize(i), 4/dimSize(i)]);
            set(handles.(['Dim' num2str(i) '_slider']),'Callback',{@MU_linkSliderEditDim,handles.(['Dim' num2str(i) '_slider']),handles.(['Dim' num2str(i) '_slider']),handles.(['Dim' num2str(i) '_edit']),i});
            set(handles.(['Dim' num2str(i) '_edit'])  ,'Callback',{@MU_linkSliderEditDim,handles.(['Dim' num2str(i) '_edit']),handles.(['Dim' num2str(i) '_slider']),handles.(['Dim' num2str(i) '_edit']),i});
        end
    end
end

% turn off uitab warning
warning('off');
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [ a, seed ] = r8pp_random ( n, seed )

%*****************************************************************************80
%
%% R8PP_RANDOM randomizes a R8PP matrix.
%
%  Discussion:
%
%    The R8PP storage format is appropriate for a symmetric positive
%    definite matrix.  Only the upper triangle of the matrix is stored,
%    by successive partial columns, in an array of length (N*(N+1))/2,
%    which contains (A11,A12,A22,A13,A23,A33,A14,...,ANN)
%
%    The matrix is computed by setting a "random" upper triangular
%    Cholesky factor R, and then computing A = R'*R.
%    The randomness is limited by the fact that all the entries of
%    R will be between 0 and 1.  A truly random R is only required
%    to have positive entries on the diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A((N*(N+1))/2), the R8PP matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  a(1:(n*(n+1))/2) = 0.0;

  for i = n : -1 : 1
%
%  Set row I of R.
%
    for j = i : n
      ij = i + ( j * ( j - 1 ) ) / 2;
      [ a(ij), seed ] = r8_uniform_01 ( seed );
    end
%
%  Consider element J of row I, last to first.
%
    for j = n : -1 : i
%
%  Add multiples of row I to lower elements of column J.
%
      ij = i + ( j * ( j - 1 ) ) / 2;

      for k = i+1 : j
        kj = k + ( j * ( j - 1 ) ) / 2;
        ik = i + ( k * ( k - 1 ) ) / 2;
        a(kj) = a(kj) + a(ik) * a(ij);
      end
%
%  Reset element J.
%
      ii = i + ( i * ( i - 1 ) ) / 2;
      a(ij) = a(ii) * a(ij);

    end
  end

  return
end


text:
function D = tsss_spm_momentspace(S)
% Switch a dataset to SSS space using virtual montage
% FORMAT D = spm_eeg_crop(S)
%
% S        - input struct
%  fields of S:
%   D          - MEEG object or filename of M/EEG mat-file with data after
%                TSSS tool
%   condthresh - threshold on condition number for regularisation
%
% Output:
% D        - MEEG object (also written on disk)
%
% Reference: Vrba J, Taulu S, Nenonen J, Ahonen A. Signal space separation
% beamformer. Brain Topogr. 2010 Jun;23(2):128-33.
%__________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: tsss_spm_momentspace.m 7703 2019-11-22 12:06:29Z guillaume $

SVNrev = '$Rev: 7703 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','TSSS momenspace transformation'); spm('Pointer','Watch');

if ~isfield(S, 'condthresh'),       S.condthresh   = 80;           end

D = spm_eeg_load(S.D);

if ~isfield(D, 'SSS')
    error('Run the TSSS tool first');
end

D = montage(D, 'clear');

SSS = D.SSS;

[SN_new, sss_indices , nmodes] = basis_condition_adjustment(SSS(1).SN, size(SSS.SNin, 2), S.condthresh);
pSN = pinv(SN_new);
pSN = pSN(1:nmodes, :);

labelnew = {};
for i = 1:nmodes
    labelnew{i} = ['moment' num2str(i)];
end

mont = [];
mont.labelorg = D.chanlabels(D.indchantype('MEGANY'))';
mont.labelnew = labelnew;
mont.tra = pSN;

if ~isempty(S.addchannels)
    chantypeorg = D.chantype(D.indchannel(S.addchannels));
    chanunitorg = D.units(D.indchannel(S.addchannels));
    for c = 1:numel(S.addchannels)
        mont.labelorg(end+1) = S.addchannels(c);
        mont.labelnew(end+1) = S.addchannels(c);
        mont.tra(end+1, end+1) = 1;
    end
end

D = montage(D, 'add', mont);
D = chantype(D, strmatch('moment', D.chanlabels), 'MEG');
D = units(D, strmatch('moment', D.chanlabels), 'fT');

if ~isempty(S.addchannels)
    D = chantype(D, D.indchannel(S.addchannels), chantypeorg);
    D = units(D, D.indchannel(S.addchannels), chanunitorg);
end

D = badchannels(D, D.indchantype('MEGANY'), 0);

save(D);

spm('FigName','TSS momenspace transformation: done'); spm('Pointer','Arrow');

author: spm


text:
function sd_estimate = sdest3(x)
 wc=dwt3(x,'sym4','mode','sym');
s=wc.dec{2,2,2};
sd_estimate = mad(s(:),1)/0.6745;

author: thomaskuestner


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%backwardstep_flow   Reference problem 5.2 inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies backward step flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=0*xbd; bcy=0*xbd;
k=find(xbd==-1); bcx(k)=4*ybd(k).*(1-ybd(k));
return


text:
% Copyright (C) 2008   Sylvain Pelissier   <sylvain.pelissier@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, see <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn {Function File} {[@var{yout}] =} bohmanwin(@var{xin},@var{h}),@var{p},@var{q})
%	Upsample, filter and downsample a signal.
% @seealso{rectwin,  bartlett}
% @end deftypefn

function yout = upfirdn(xin,h,p,q)

if(nargin < 2)
  error('usage : yout = upfirdn(xin,h,p,q)');
end

if(nargin < 3)
	p = 1;
	q = 1;
end

if(nargin < 4)
	q = 1;
end

if(floor(p) ~= p || floor(q) ~= q || p < 1 || q < 1)
	error('p and q must be positive integer');
end

yout = upsample(xin,p);
yout = convn(yout, h).*p; % original was filter(h, 1, yout);
% the scaling with p is needed as per github issue 2085, causing the output
% to be scaled by the value of p, with this change, the compat/matlab
% versions will give an output that is about equal (scaled with about 0.9993)
yout = downsample(yout,q);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ferns,hsPr] = fernsClfChangeFernDepth( data, hs, ferns, Snew, varargin )
% Train random fern classifier.
%
% See "Fast Keypoint Recognition in Ten Lines of Code" by Mustafa Ozuysal,
% Pascal Fua and Vincent Lepetit, CVPR07.
%
% Dimensions:
%  M - number ferns
%  S - fern depth
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  [ferns,hsPr] = fernsClfTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .S        - [10] fern depth (ferns are exponential in S)
%   .M        - [50] number of ferns to train
%   .thrr     - [0 1] range for randomly generated thresholds
%   .bayes    - [1] if true combine probs using bayes assumption
%   .ferns    - [] if given reuse previous ferns (recompute pFern)
%
% OUTPUTS
%  ferns    - learned fern model w the following fields
%   .fids     - [MxS] feature ids for each fern for each depth
%   .thrs     - [MxS] threshold corresponding to each fid
%   .pFern    - [2^SxHxM] learned log probs at fern leaves
%   .bayes    - if true combine probs using bayes assumption
%   .inds     - [NxM] cached indices for original training data
%   .H        - number classes
%  hsPr     - [Nx1] predicted output labels
%
% EXAMPLE
%  N=5000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  fernPrm=struct('S',4,'M',50,'thrr',[-1 1],'bayes',1);
%  tic, [ferns,hsPr0]=fernsClfTrain(xs0,hs0,fernPrm); toc
%  tic, hsPr1 = fernsClfApply( xs1, ferns ); toc
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also fernsClfApply, fernsInds
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2010 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

% get additional parameters and check dimensions
dfs={'thrr',[0 1],'tmp',''};
[thrr,~]=getPrmDflt(varargin,dfs,1);

[M,Sold] = size(ferns.fids);
[N,F]=size(data); assert(length(hs)==N);
H=max(hs); assert(all(hs>0)); assert(Snew<=20);

if Snew == Sold,
elseif Snew < Sold,
  Sremove = Sold - Snew;

  % remove classifiers from ferns
  ferns.thrs(:,Snew+1:Sold) = [];
  ferns.fids(:,Snew+1:Sold) = [];

  % remove corresponding bits from inds
  ferns.inds = bitshift(ferns.inds,-Sremove);

  % accumulate counts into parent nodes
  ferns.counts = reshape(sum(reshape(ferns.counts,[2^Snew,2^Sremove,H,M]),2),[2^Snew,H,M]);

  % convert fern leaf class counts into probabilities
  if( ferns.bayes<=0 )
    norm = 1./sum(ferns.counts,2);
    ferns.pFern = bsxfun(@times,ferns.counts,norm);
  else
    norm = 1./sum(ferns.counts,1);
    ferns.pFern = bsxfun(@times,ferns.counts,norm);
    ferns.pFern=log(ferns.pFern);
  end
else

  Sadd = Snew - Sold;

  % generate some new classifiers for each fern
  ferns.thrs(:,Sold+1:Snew) = rand(M,Sadd)*(thrr(2)-thrr(1))+thrr(1);
  ferns.fids(:,Sold+1:Snew) = uint32(floor(rand(M,Sadd)*F+1));

  % new lower bits to inds
  indsadd=fernsInds(data,ferns.fids(:,Sold+1:Snew),ferns.thrs(:,Sold+1:Snew));
  ferns.inds = bitshift(ferns.inds,Sadd) + indsadd;

  % re-histogram
  ferns.counts = zeros(2^Snew,H,M);
  edges = 1:2^Snew;
  for h=1:H,
    inds1=ferns.inds(hs==h,:);
    for m=1:M,
      ferns.counts(:,h,m)=histc(inds1(:,m),edges);
    end
  end
  ferns.counts = ferns.counts + ferns.bayes;

  % convert fern leaf class counts into probabilities
  if( ferns.bayes<=0 )
    norm = 1./sum(ferns.counts,2);
    ferns.pFern = bsxfun(@times,ferns.counts,norm);
  else
    norm = 1./sum(ferns.counts,1);
    ferns.pFern = bsxfun(@times,ferns.counts,norm);
    ferns.pFern=log(ferns.pFern);
  end
end

if(nargout==2), hsPr=fernsClfApply([],ferns,inds); end

end

author: kristinbranson

repo: JAABA


text:
function [ scmap_all,poly ] = get_sticks_segmentation( p, im, joints )

stride = 4; %p.stride;
half_stride = stride/2;
scale_factor = p.scale_factor;
sz = 17;

scmap_height = ceil(size(im, 1) * scale_factor / stride);
scmap_width = ceil(size(im, 2) * scale_factor / stride);

joint_pairs = [1 2; 2 3; 6 5; 4 5; 7 8; 8 9; 12 11; 11 10; 13 14];
limb_size_coefs = [1.0 1.0 1.0 1.0 0.8 0.8 0.8 0.8 1.0];

num_sticks = size(joint_pairs, 1) + 1;

scmap_all = zeros(scmap_height, scmap_width, num_sticks);

for k = 1:size(joint_pairs, 1)
    scmap = zeros(scmap_height, scmap_width);
    limb_sz = sz * limb_size_coefs(k);

    jnt1 = joints(joint_pairs(k, 1), :);
    jnt2 = joints(joint_pairs(k, 2), :);
    if ~isnan(jnt1(1)) && ~isnan(jnt2(1))

        % stick heatmap
        diff = jnt2-jnt1;
        if norm(diff) > 1.0
            perp = [-diff(2) diff(1)];
            perp = perp/norm(perp);
            % construct stick polygon
            poly = zeros(5, 2, 'double');
            poly(1,:) = jnt1-perp*limb_sz;
            poly(2,:) = jnt1+perp*limb_sz;
            poly(3,:) = jnt2+perp*limb_sz;
            poly(4,:) = jnt2-perp*limb_sz;
            poly(5,:) = poly(1,:);

            for j = 1:scmap_height
                for i = 1:scmap_width

                    crd = [i-1, j-1]*stride;
                    if p.res_net
                        crd = crd + half_stride;
                    end
                    crd = single(crd)/scale_factor;
                    [in, on] = inpolygon(crd(1),crd(2),poly(:,1),poly(:,2));
                    scmap(j, i) = in | on;
                end
            end

            if k ~= 9
                for j = 1:scmap_height
                    for i = 1:scmap_width

                        crd = [i-1, j-1]*stride;
                        if p.res_net
                            crd = crd + half_stride;
                        end
                        crd = single(crd)/scale_factor;

                        if norm(crd-jnt1) <= limb_sz
                            scmap(j, i) = 1;
                        end
                        if norm(crd-jnt2) <= limb_sz
                            scmap(j, i) = 1;
                        end
                    end
                end
            end
        end
    end

    scmap_all(:,:,k) = scmap(:,:);
end

% torso
jnt = zeros(4, 2);
joints = int32(joints);
jnt1 = double(joints(3, :));
jnt2 = double(joints(4, :));
jnt3 = double(joints(9, :));
jnt4 = double(joints(10, :));

if ~isnan(jnt1(1)) && ~isnan(jnt2(1)) && ~isnan(jnt3(1)) && ~isnan(jnt4(1))
    points = zeros(1, 2);

    index = 1;

    if all(jnt1 == jnt2)
        jnt2(1) = jnt1(1) + 1;
    end
    diff12 = normalise(jnt2-jnt1);
    points(index, :) = jnt2 + diff12*sz;
    index = index + 1;
    points(index, :) = jnt1 - diff12*sz;
    index = index + 1;

    if all(jnt1 == jnt3)
        jnt3(2) = jnt1(2) - 1;
    end
    diff13 = normalise(jnt3-jnt1);
    points(index, :) = jnt3 + diff13*sz;
    index = index + 1;
    points(index, :) = jnt1 - diff13*sz;
    index = index + 1;

    if norm(jnt3 - jnt4) <= sz*1.5
        if all(jnt4 == jnt3)
            jnt4(1) = jnt3(1) + 1;
        end

        diff34 = normalise(jnt4-jnt3);
        points(index, :) = jnt4 + diff34*sz;
        index = index + 1;
        points(index, :) = jnt3 - diff34*sz;
        index = index + 1;
    end

    if all(jnt2 == jnt4)
        jnt4(2) = jnt2(2) - 1;
    end
    diff24 = normalise(jnt4-jnt2);
    points(index, :) = jnt4 + diff24*sz;
    index = index + 1;
    points(index, :) = jnt2 - diff24*sz;
    index = index + 1;

    I = convhull(points(:,1), points(:, 2));
    poly = points(I, :);

    scmap = zeros(scmap_height, scmap_width);
    for j = 1:scmap_height
        for i = 1:scmap_width

            crd = [i-1, j-1]*stride;
            if p.res_net
                crd = crd + half_stride;
            end
            crd = single(crd)/scale_factor;
            [in, on] = inpolygon(crd(1),crd(2),poly(:,1),poly(:,2));
            scmap(j, i) = in | on;
        end
    end

    scmap_all(:,:,size(scmap_all, 3)) = scmap(:,:);
end

end

function out = normalise(vec)
vec = double(vec);
n = norm(vec);
if n <= 1
end
out = vec/norm(vec);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function latexcmd_example(varargin)

% Usage: latexcmd_example(outfilename)
%
% Run this command and then look at this file and the resulting LaTeX file
% (latexcmd_example_out.tex) to better understand how latexcmd can be used.
%
% If you find any errors, please let me know! (peder at axensten dot se)
%
% Example:
% latexcmd_example('testout.tex')
%
% Copyright (C) Peder Axensten (peder at axensten dot se), 2006.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	outfile=	'latexcmd_example_out';		% Default file name
	if(nargin > 0)
		if((nargin == 1) && ischar(varargin{1}))
			% File name is given.
			outfile= varargin{1};

		% Run multiple times, for timing purposes.
		elseif((nargin == 1) && isnumeric(varargin{1}))
			iterations=	varargin{1};
			outfile=	tempname();
			tic;
			for n= 1:iterations
				latexcmd_example(outfile);
				fprintf(1, 'Iteration %d (of %d)\r', n, iterations);
			end
			fprintf(1, 'Finished %d iterations in  %f seconds.\n', iterations, toc);
			delete(outfile);
		else
			error('Usage: latexcmd_example(outfilename)');
		end
	end

	latexcmd(outfile);
	outfile=	['>' outfile];

	% Symbolic Toolbox.
	% Exporting symbolic variables works in the same way as numeric variables.
	if(exist('sym', 'file') == 2)
		% Create 3 symbolic variables.
		x=			sym('x');
		y=			sym('y');
		z=			sym('z');
		% Calculate the Hessian, a 3x3 symbolic array.
		symbolic=	jacobian(jacobian(x^4 + y^3*x + z^2*x^2));
		% Get a symbolic scalar, for later purposes.
		symitem=	symbolic(1, 1);
	else
		symbolic=	'No symbolic toolbox...';
		symitem=	'No symb.';
	end
	latexcmd(outfile, symbolic);

	% Name an expression.
		latexcmd(outfile, '\mypi', 3.1415+0.000092);

	% Check various container forms (array, cell array, character string, structure).
		% Zero sized (this is for testing purposes).
		stringzero=				'';
		arrayzero=				[];
		cellzero=				{};
		structzero.first=		[];
		latexcmd(outfile, stringzero, arrayzero, cellzero, structzero);

		% One item sized (this is for testing purposes).
		stringone=				'A';
		arrayone=				 1.1;
		cellone=				{1.1};
		structone.first=		1.1;		% Numeric.
		structone.second=		'a.b';		% String.
		structone.third=		symitem;	% Symbolic.
		latexcmd(outfile, stringone, arrayone, cellone, structone);

		% Three items sized.
		stringnorm=				'Normal sized string';
		arraynorm=				[1.1  1.2      1.3];
		cellnorm=				{1.1 '1.2' symitem};	% Numeric, string, and symbolic items.
		structnorm(1).first=	1.1;		% Numeric.
		structnorm(1).second=	1.2;		% Numeric.
		structnorm(1).third=	1.3;		% Numeric.
		structnorm(2).first=	2.1;		% Numeric.
		structnorm(2).second=	'b.b';		% String.
		structnorm(2).third=	symitem;	% Symbolic.
		latexcmd(outfile, stringnorm, arraynorm, cellnorm, structnorm);

	% Check different formats.
		% Create array of values.
		nums=		    [-0, -9.8765e-113, -1.234e-7, -0.5-0.125i, -23.4, -1.234e7, -9.8765e113, -Inf];
		nums=			[nums(end:-1:1), NaN, 1+ NaN*i];
		nums=			[nums, +0,  9.8765e-113,  1.234e-7,  0.5+0.125i,  23.4,  1.234e7,  9.8765e113,  Inf];
		% Various output formats, look at the definition of '\numforms' to see what happens!
		numformsformat=	'%g & %3s & %3S & %r & %R & %12.1E & %12.2e & %d';
		numformsarr=	eval(['{''' regexprep(numformsformat, '\s*&\s*', ''',''') '''}']);
		numforms=		transpose([nums; nums; nums; nums; nums; nums; nums; nums; nums]);

		% First 'numformstr' is exported as a string, the second decides the numerical format of 'numforms'.
		latexcmd(outfile, '\numformstr', numformsarr, numformsformat, numforms, '%g');

	% Check "special" numbers.
		% Complex
		numscomplex=	transpose([-1-i, -i, i, 1+i; ...
						 	-Inf-Inf*i, 1-Inf*i, 1+Inf*i, Inf+Inf*i; ...	% Matlab makes real part a NaN!??
						  	Inf+i, -Inf-i, NaN, ...
											  	NaN*i ...					% Matlab makes real part a NaN!??
						]);

		% Polynomial (actually a vector of six polynomials).
		polynomialA=	[ 	 0,  1, -1+i,  2, -2; ...
					 		-2,  0,  1-i, -1,  2; ...
					  		 2, -2,  0+i,  1, -1; ...
					 		-1,  2,   -2,  0,  1-i; ...
					  		 1, -1,    2, -2,  0; ...
					  		 0,  0,    0,  0,  0];
		polynomialB=	polynomialA;

		latexcmd(outfile, numscomplex, '%poly:\alpha', polynomialA, ...
					'%3S & %.2e & %.1f & %d & %R', '%poly:t', polynomialB);

	% Statistics AND using the prefix option.
		arraydd=				[1.1, 1.2,     1.3; 2.1, 2.2, 2.3; 3.1,  3.2,  3.3];
		latexcmd(outfile, arraydd, '@prefix', '%stat', arraydd);

end


text:
% ir_example_ct_lir1
% example of local impulse response (LIR) in CT
% Copyright 2012-07-28, Jeff Fessler, University of Michigan

if ~isvar('A'), printm 'setup geometry, image, sinogram'
	f.down = 4;
	ig = image_geom('nx', 512, 'fov', 30, 'down', f.down);
	ig.mask = ig.circ > 0;
	sg = sino_geom('par', 'nb', ig.nx, 'na', ig.nx, ...
			'dr', ig.dx, 'strip_width', 'd');
	sg.plot(ig);

	ell = [0 0 10 10 0 0.2;
                 0  6 3 3 0  0.2;
                 6  0 3 3 0  0.1
                 0 -6 3 3 0 -0.1;
                -6  0 3 3 0 -0.2;
		];
        xtrue = ellipse_im(ig, ell, 'oversample', 2);

%	A = Gtomo2_dscmex(sg, ig);
	A = Gtomo2_wtmex(sg, ig, 'nthread', jf('ncore'));
	ytrue = A * xtrue;
	wi = exp(-ytrue); % ideal Poisson weighting

	im plc 2 3
	clim = [0 0.4];
	im(1, xtrue, 'x', clim), cbar
	im(2, ytrue, 'y: ideal sinogram'), cbar
	im(3, wi, 'w: ideal weighting'), cbar
%	ir_savefig ir_example_ct_lir1_x_y_w
prompt
end

if ~isvar('yp'), printm 'yp'
	ix = round(ell(2:end, 1) / ig.dx + (ig.nx-1)/2);
	iy = round(ell(2:end, 2) / ig.dy + (ig.ny-1)/2);
	xp = ig.zeros;
	for ii=1:numel(ix)
		xp = xp + ig.unitv(ix(ii), iy(ii));
	end
	f.amp = 0.1;
	xp = f.amp * xp ;
	im(2, xtrue + xp)

	yp = A * xp;
	climp = [0 f.amp];
	im(3, xp, climp), cbar
	im(4, yp)
end

if ~isvar('fbp'), printm 'fbp'
	tmp = fbp2(sg, ig);
	fbp = fbp2(yp, tmp, 'window', 'boxcar,0.8');
	im(2, fbp, 'FBP', climp), cbar
prompt
end

if ~isvar('fw.fbp'), printm 'fw.fbp'
	ox = [-9:9]; oy = ox;
	fw.fun = @(x, ii) fwhm2(x(ix(ii)+ox, iy(ii)+oy));
	for ii=1:numel(ix)
		fw.fbp(ii) = fw.fun(fbp, ii);
	end
	pr fw.fbp
end

if ~isvar('kappa'), printm 'kappa: try to make resolution approximately uniform'
	kappa = sqrt( div0(A' * wi, A' * ones(size(wi))) );
	im(4, kappa), cbar
prompt
end

% use local psf to help select beta
if ~isvar('R1'), printm 'R1, R2'
	f.l2b = 0;
	W = Gdiag(wi);
	R1 = Reg1(ig.mask * kappa(end/2+1,end/2+1), 'beta', 2^f.l2b); % usual
%	qpwls_psf(A, R1, 1, ig.mask, W, 'loop', 1); % choose beta
	R2 = Reg1(kappa, 'beta', 2^f.l2b); % kappa
	qpwls_psf(A, R2, 1, ig.mask, W, 'loop', 1); % choose beta
prompt
end

%	psf = qpwls_psf(A, R1, 1, ig.mask, W);
%	init = conv2(psf, xp, 'same');
	init = fbp;
	im(init, climp)

if ~isvar('xpwls1'), printm 'pwls1'
	f.niter = 300;
	xpwls1 = pwls_pcg1(init(ig.mask), A, W, yp(:), R1, 'niter', f.niter);
	xpwls1 = ig.embed(xpwls1);
	im(3, xpwls1, 'PWLS R1'), cbar
end

if ~isvar('xpwls2'), printm 'pwls2'
	xpwls2 = pwls_pcg1(init(ig.mask), A, W, yp(:), R2, 'niter', f.niter);
	xpwls2 = ig.embed(xpwls2);
	im(4, xpwls2, 'PWLS R2'), cbar
end

if 1
	for ii=1:numel(ix)
		fw.pwls1(ii) = fw.fun(xpwls1, ii);
		fw.pwls2(ii) = fw.fun(xpwls2, ii);
	end
	pr fw.fbp
	pr fw.pwls1
	pr fw.pwls2
end

if 1
	climp = [-0.004 0.05];
%	climp = []
	im plc 1 3
	ax = [18 110 20 120];
	im(1, fbp, 'FBP', climp), axis(ax)%, cbar h
	fun = @(fw, ii) text(ix(ii), iy(ii)+15, ...
		sprintf('%3.1f', fw(ii)), 'horiz', 'center', 'color', 'green');
	for ii=1:numel(ix), fun(fw.fbp, ii); end
	im(2, xpwls1, 'PWLS standard R', climp), axis(ax)%, cbar h
	for ii=1:numel(ix), fun(fw.pwls1, ii); end
	im(3, xpwls2, 'PWLS modified R', climp), axis(ax)%, cbar h
	for ii=1:numel(ix), fun(fw.pwls2, ii); end

%	ir_savefig eps_c ir_example_ct_lir1_fwhm
end


text:
filename='ImprovedBridge_hexahedra';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'IPOPT'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [F] = Filter2D(Norder,Nc,sp)

% function [F] = Filter2D(Norder,sp)
% Purpose : Initialize 2D filter matrix of order sp and cutoff Nc

Globals2D;

filterdiag = ones((Norder+1)*(Norder+2)/2,1);
alpha = -log(eps);

% build exponential filter
sk = 1;
for i=0:Norder
  for j=0:Norder-i
    if (i+j>=Nc)
      filterdiag(sk) = exp(-alpha*((i+j - Nc)/(Norder-Nc))^sp);
    end
    sk = sk+1;
  end
end

F = V*diag(filterdiag)*invV;
return;

repo: Aero-matlab


text:
function coeffs=terms2MultiDimPolyMat(termMat)
%%TERMS2MULTIDIMPOLYMAT Given a multidimensional polynomial represented as
%                    a 2D matrix of terms (monomials), convert it into a
%                    hypermatrix of coefficients suitable for use with
%                    functions such as polyValMultiDim and where convn can
%                    be used to multiply two such polynomials.
%
%INPUTS: termMat An (n+1)XnumTerms matrix such that
%                termMat(:,i)=[c,a1,a2,...,an] is a monomial term where c
%                is the value of of the monomial coefficient and the
%                monomial is
%                x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1). The ordering
%                of the terms in termMat does not matter.
%
%OUTPUTS: coeffs A hypermatrix of the coefficients for the multivariate
%               polynomial. These are arranged such that
%               coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%               x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%               number of indices coeffs takes is equal to the
%               dimensionality of x (not counting singleton dimensions at
%               the end of coeffs). Note that this ordering is the reverse
%               that used in the 1D polyval function that is built into
%               Matlab. The number of elements for each index in coeffs is
%               the maximum order of that dimension +1.
%
%This function is the opposite of multiDimHyperMat2Terms. See
%multiDimHyperMat2Terms for usage examples.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numTerms=size(termMat,2);
numIdx=size(termMat,1)-1;

numDims=zeros(1,numIdx);

if(isscalar(numDims))
    numDims=[numDims,1];
end

for idx=1:numIdx
    numDims(idx)=max(termMat(idx+1,:))+1;
end

coeffs=zeros(numDims);
for curTerm=1:numTerms
    idx=nDim2Index(numDims,termMat(2:end,curTerm)+1);
    coeffs(idx)=coeffs(idx)+termMat(1,curTerm);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function ival = s_to_i4 ( s )

%*****************************************************************************80
%
%% S_TO_I4 reads an integer value from a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string to be examined.
%
%    Output, integer IVAL, the integer value read from the string.
%
  sgn = 1;
  state = 0;
  ival = 0;

  i = 0;

  while ( i < s_len_trim ( s ) )

    i = i + 1;
    c = s(i);

    if ( state == 0 )

      if ( c == ' ' )

      elseif ( c == '-' )
        state = 1;
        sgn = -1;
      elseif ( c == '+' )
        state = 1;
        sgn = +1;
      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return;
      end
%
%  Have read the sign, now expecting the first digit.
%
    elseif ( state == 1 )

      if ( c == ' ' )

      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return
      end
%
%  Have read at least one digit, expecting more.
%
    elseif ( state == 2 )

      if ( '0' <= c & c <= '9' )
        ival = 10 * ival + c - '0';
      else
        ival = sgn * ival;
        return;
      end

    end

  end
%
%  If we read all the characters in the string, see if we're OK.
%
  if ( state ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'S_TO_I4 - Fatal error!\n' );
    fprintf ( '  Did not read enough information to define an integer!\n' );
    return;
  end

  ival = sgn * ival;

  return
end


text:
function chrpak_test115 ( )

%*****************************************************************************80
%
%% TEST115 tests S_REPLACE_CH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
  s_old = 'No pennies now.';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST115\n' );
  fprintf ( 1,'  S_REPLACE_CH replaces one character by another;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    C1  C2  Original String  Modified String\n' );
  fprintf ( 1, '\n' );

  c1 = 'n';
  c2 = 't';

  s = s_replace_ch ( s_old, c1, c2 );

  fprintf ( 1, '     %c   %c   %s  %s\n', c1, c2, s_old, s );

  return
end


text:
function DoToHDF5(Simuh)

% convert MRiLab output to HDF5 file which works with Gadgetron
global VCtl
global VSig
global VCoi
global VObj

vctl.protocolName               = VCtl.SeriesName;
vctl.systemVendor               = 'MRiLab';
vctl.systemModel                = '2013a';
vctl.systemFieldStrength_T      = single(VCtl.B0);
vctl.receiverChannels           = uint16(VCoi.RxCoilNum);
vctl.institutionName            = 'WIMR';
vctl.stationName                = 'L1122-E';
vctl.H1resonanceFrequency_Hz    = uint32(VCtl.B0 * VObj.Gyro/(2*pi));
vctl.ESMatrixSizeX              = uint16(VCtl.ResFreq);
vctl.ESMatrixSizeY              = uint16(VCtl.ResPhase);
vctl.ESMatrixSizeZ              = uint16(VCtl.SliceNum);
vctl.ESFOVX                     = single(VCtl.FOVFreq  * 1000);
vctl.ESFOVY                     = single(VCtl.FOVPhase * 1000);
vctl.ESFOVZ                     = single(VCtl.FOVSlice * 1000);
vctl.RSMatrixSizeX              = uint16(VCtl.ResFreq);
vctl.RSMatrixSizeY              = uint16(VCtl.ResPhase);
vctl.RSMatrixSizeZ              = uint16(VCtl.SliceNum);
vctl.RSFOVX                     = single(VCtl.FOVFreq  * 1000);
vctl.RSFOVY                     = single(VCtl.FOVPhase * 1000);
vctl.RSFOVZ                     = single(VCtl.FOVSlice * 1000);
vctl.trajectory                 = VCtl.TrajType;
vctl.TR                         = single(VCtl.TR);
vctl.TE                         = single(VCtl.TE);
vctl.BW                         = single(VCtl.BandWidth);
vctl.outputFile                 = [Simuh.OutputDir filesep 'Series' num2str(Simuh.ScanSeriesInd) '.h5'];

Sx=single(sum(reshape(VSig.Sx, length(VSig.Sx)/VObj.TypeNum, VObj.TypeNum),2));
Sy=single(sum(reshape(VSig.Sy, length(VSig.Sy)/VObj.TypeNum, VObj.TypeNum),2));
Kx=single(VSig.Kx);
Ky=single(VSig.Ky);
Kz=single(VSig.Kz);
vsig.S                          = [Sx'; Sy'];
vsig.K                          = [Kx; Ky; Kz];
vsig.echoNumber                 = uint32(VCtl.TEPerTR);
vsig.firstPhaseNumber           = uint32(VCtl.FirstPhNum);
vsig.secondPhaseNumber          = uint32(VCtl.SecondPhNum);
vsig.readoutNumber              = uint32(length(VSig.Sx)/(VCtl.TEPerTR*VCtl.FirstPhNum*VCtl.SecondPhNum*VCoi.RxCoilNum*VObj.TypeNum));

% convert MRiLab output to HDF5 file
DoMatToHDF5(vsig,vctl);

end


text:
function [hd_record,hd_estimates]=panel6hd(Xi,theta_gibbs,D_record,strshocks_record,It,Bu,Ymat,N,n,m,p,k,T,d,HDband)

% preliminary tasks
% first create the hd_record and temp cells
hd_record=cell(n*N,n*N+1);
temp=cell(n*N,2);

% then initiate the Gibbs algorithm
for ii=1:It-Bu

% recover theta for the current iteration (one column for each period)
theta_iter=reshape(theta_gibbs(:,ii,:),d,T);

% recover D for the current iteration (one page for each period)
D_iter=[];
for jj=1:T
D_iter(:,:,jj)=reshape(D_record(:,ii,jj),N*n,N*n);
end

% recover the series of period-specific orthogonal IRFs
[IRFcell]=bear.panel6hdsim(theta_iter,D_iter,Xi,N,n,m,p,T,k);

% recover the structural disturbances
ETA=[];
for jj=1:N*n
ETA=[ETA;strshocks_record{jj,1}(ii,:)];
end
ETA=ETA';

% then compute the historical decomposition
   % loop over variables
   for jj=1:n*N
      % loop over shocks
      for kk=1:n*N
      % loop over shocks
      vshocks=ETA(:,kk);
         % loop over time periods
         for ll=1:T
         % initiate the vectors of IRFs and shocks
         virf=[];
            % then loop over IRF periods
            for mm=1:ll
            % create the vector of IRF coefficients
            virf=[virf IRFcell{mm,ll}(jj,kk)];
            end
         % compute then the contribution of shock kk for variable jj at period ll
         hd_record{jj,kk}(ii,ll)=virf*flipud(vshocks(1:ll,1));
         end
      end
   end

% then go for next Gibbs iteration
end

% compute the contributions of deterministic variables
% loop over rows of temp/hd_record
for ii=1:n*N
% fill the Ytot matrix in temp
% initial condition
temp{ii,1}=hd_record{ii,1};
   % sum over the remaining columns of hd_record
   for jj=2:n*N
   temp{ii,1}=temp{ii,1}+hd_record{ii,jj};
   end
% fill the Y matrix in temp
temp{ii,2}=repmat(Ymat(:,ii)',It-Bu,1);
% fill the Yd matrix in hd_record
hd_record{ii,N*n+1}=temp{ii,2}-temp{ii,1};
% go for next variable
end

% finally, obtain point esimates and credibility intervals
[hd_estimates]=bear.hdestimates(hd_record,N*n,T,HDband);

author: european-central-bank

repo: BEAR-toolbox