text:
clear();

CONVOLUTION_SHAPE_FULL  = 1;
CONVOLUTION_SHAPE_SAME  = 2;
CONVOLUTION_SHAPE_VALID = 3;

maxThr = 1e-9;

tic();
for numElementsSignal = 8:21

    vI = rand(numElementsSignal, 1);

    for numElementsKernel = 1:7

        vH = rand(numElementsKernel, 1);

        for convShape = 1:3

            switch(convShape)
                case(CONVOLUTION_SHAPE_FULL)
                    numElementsOut  = numElementsSignal + numElementsKernel - 1;
                    convShapeString = 'full';
                case(CONVOLUTION_SHAPE_SAME)
                    numElementsOut  = numElementsSignal;
                    convShapeString = 'same';
                case(CONVOLUTION_SHAPE_VALID)
                    numElementsOut  = numElementsSignal - numElementsKernel + 1;
                    convShapeString = 'valid';
            end

            vORef   = conv2(vI, vH, convShapeString);
            mK      = GenerateToeplitzConvMatrixEx(vH, numElementsSignal, convShape);
            % mKK      = full(CreateConvMtx1D(vH, numElementsSignal, convShape));
            vO      = reshape(mK * vI, numElementsOut, 1);

            disp([' ']);
            disp(['Validating solution for the following parameters:']);
            disp(['Signal Size - [', num2str(numElementsSignal), ' x 1]']);
            disp(['Kernel Size - [', num2str(numElementsKernel), ' x 1]']);
            disp(['Convolution Shape - ', convShapeString]);

            vE = vO - vORef;
            maxAbsDev = max(abs(vE(:)));
            if(maxAbsDev >= maxThr)
                disp([' ']);
                disp(['Validation Failed']);
                disp([' ']);
            end
            assert(maxAbsDev < maxThr);

        end
    end
end

toc();

author: RoyiAvital

repo: StackExchangeCodes


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Laplace Transform properties

% Time shifting

syms t s
t0=2;

Le=cos(t-t0)*heaviside(t-t0);
Left= laplace(Le)

X=laplace(cos(t),s);
Right=exp(-s*t0)*X

author: Sable

repo: mcbench-benchmarks


text:
function [XNorm mu stddev] = normalizeFeatures(X)

% This function is used to normalize the features. It makes the mean of
% each feature 0. Formula for normalizing feature taking all examples into
% consideration is : [X(particular feature) - mean(of feature)]/ standard
% deviation of feature.

mu = mean(X);
stddev = std(X);

XNorm = bsxfun(@minus,X,mu); % subtract mean of each feature from original value
XNorm = bsxfun(@rdivide,XNorm,stddev); % divide by standard deviation


text:
function errFolders = db_delete_subjects( iSubjects )
% DB_DELETE_SUBJECTS: Delete some subjects from the brainstorm database.
%
% USAGE:  errFolders = db_delete_subjects( iSubjects )
%
% INPUT:
%    - iSubjects : indices of the subejcts to delete
%
% OUTPUT:
%    - errFolders : Cell-array with full paths to folders that could not be deleted

% @=============================================================================
% 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, 2012-2022

ProtocolInfo     = bst_get('ProtocolInfo');
ProtocolSubjects = bst_get('ProtocolSubjects');
errFolders = {};

% Remove invalid indices and default subject
iInvalid = find(iSubjects <= 0);
if ~isempty(iInvalid)
    disp('DELETE> Cannot delete default subject.');
    iSubjects(iInvalid) = [];
end

% For each subject
for i = 1:length(iSubjects)
    SubjectFile = ProtocolSubjects.Subject(iSubjects(i)).FileName;
    % === DELETE STUDIES ===
    % Find all the studies that are associated with the current brainstormsubject file
    [sStudies, iStudies] = bst_get('StudyWithSubject', SubjectFile, 'intra_subject', 'default_study');
    if ~isempty(sStudies)
        % Delete folders
        errFolders = cat(2, errFolders, db_delete_studies(iStudies));
        % Delete the studies folder for the subject
        subjDataDir = bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(bst_fileparts(sStudies(1).FileName)));
        if (file_delete(subjDataDir, 1, 1) ~= 1)
            errFolders{end+1} = subjDataDir;
            return;
        end
    end
    % === DELETE SUBJECT ===
    % Remove subject's directory
    subjAnatDir = bst_fullfile(ProtocolInfo.SUBJECTS, bst_fileparts(SubjectFile));
    if (file_delete(subjAnatDir, 1, 1) ~= 1)
        errFolders{end+1} = subjAnatDir;
        return;
    end
end
% Remove subjects from database
ProtocolSubjects.Subject(iSubjects) = [];
bst_set('ProtocolSubjects', ProtocolSubjects);
% Update tree
panel_protocols('UpdateTree');
% Save database
db_save();

author: brainstorm-tools


text:
function outVec = vecInsert(elemToInsert,position,vecIn)

%Inserts a vector of elements into the given positions.  elemToInsert is a
%vector that can be as small as 1 element and as large as desired.
%Position is an array of positions that can be as small as one element and
%as large as the number of elements you wish to insert.
%
%If position is smaller than the number of elements, this function will
%assume that the remaining elements will be inserted sequentially.
%
%If position is larger than the number of elements, this function will insert the
%last element of elemToInsert at the remaining given positions.

N = max(length(elemToInsert),length(position));
tmpVec = zeros(length(vecIn) + N,1);
tmpVec(:,1) = NaN;
if length(position) == length(elemToInsert)
    tmpVec(position) = elemToInsert;
elseif length(position) < length(elemToInsert)
    %More elements to insert than given positions
    tmpVec(position,1) = elemToInsert(1:length(position));
    tmpVec(position(end) + 1:position(end) + (length(elemToInsert) - length(position))) = elemToInsert(length(position)+1:end);
else
    %More positions than elements to insert
    tmpVec(position(1:length(elemToInsert))) = elemToInsert;
    tmpVec(position(length(elemToInsert)+1:end)) = elemToInsert(end);
end

%put vecIn back in
tmpVec(find(isnan(tmpVec(:)))) = vecIn;

outVec = tmpVec;


text:
function[varargout]=deg180(varargin)
%DEG180  Converts degrees to the range [-180,180].
%
%   [D1,D2,...,DN]=DEG180(D1,D2,...,DN) converts the input angles, which
%   are measured in degrees, to the range [-180, 180].
%
%   DEG180 also works if D1,D2,...,DN are cell arrays of numerical arrays.
%
%   See also JDEG2RAD, JRAD2DEG, DEG180, DEGUNWRAP.
%
%   'deg180 --t' runs a test.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2007--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    deg180_test,return
end

varargout=varargin;
for i=1:nargin;
    if ~iscell(varargin{1})
        varargout{i}=deg360(varargin{i});
        bool=varargout{i}>180;
        varargout{i}(bool)=varargout{i}(bool)-360;
    else
        for j=1:length(varargin{1})
            varargout{i}{j,1}=deg360(varargin{i}{j});
            bool=varargout{i}{j}>180;
            varargout{i}{j}(bool)=varargout{i}{j}(bool)-360;
        end
    end
    %varargout{i}=angle(exp(sqrt(-1)*varargin{i}*2*pi/360))*360/2/pi;
    %varargout{i}(~isfinite(varargin{i}))=varargin{i}(~isfinite(varargin{i}));
    %    varargout{i}=jrad2deg(jdeg2rad(varargin{i})); Same but slower
end

function[]=deg180_test
thi=[359 181 nan inf];
tho=[-1 -179 nan inf];
tol=1e-10;
reporttest('DEG180 simple',aresame(deg180(thi),tho,tol))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% Data file PILZ3
% Forced damped vibrations of a mechanical
% system with one degree of freedom, contained
% a reverse pendulum.
  s     = 1; % degree of freedom
  L     = '2.5*qt1^2 - 1/2*c*q1^2 - 3.9*cos(10*q1)'; % Lagrangian
  QN{1} = '-k*qt1 + 10*sin(p*t)'; % generalized non potential force
  Tend  = 20;    % upper bound of integration
  qj0   = 0.05;  % initial coordinate
  qtj0  = 0;     % initial velocity
  eps   = 1e-10; % desirable accuracy
  np    = 3;     % number of parameters
  P{1}  = 'c';   % spring stiffness
  P{2}  = 'k';   % coefficient of damping
  P{3}  = 'p';   % disturbance frequency


text:
classdef ANSGAIII < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation> <constrained/none>
% Adaptive NSGA-III

%------------------------------- Reference --------------------------------
% H. Jain and K. Deb, An evolutionary many-objective optimization algorithm
% using reference-point based non-dominated sorting approach, part II:
% Handling constraints and extending to an adaptive approach, IEEE
% Transactions on Evolutionary Computation, 2014, 18(4): 602-622.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate the reference points and random population
            % All the reference points
            [Z,Problem.N] = UniformPoint(Problem.N,Problem.M);
            Z = sortrows(Z);
            % Distance between two consecutive reference points for the adaption
            interval = Z(1,end) - Z(2,end);
            % Initial population
            Population = Problem.Initialization();
            % Ideal point
            Zmin = min(Population(all(Population.cons<=0,2)).objs,[],1);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,sum(max(0,Population.cons),2));
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Zmin       = min([Zmin;Offspring(all(Offspring.cons<=0,2)).objs],[],1);
                Population = EnvironmentalSelection([Population,Offspring],Problem.N,Z,Zmin);
                Z          = Adaptive(Population.objs,Z,Problem.N,interval);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function pot = dpot(domain, sizes, T)
% DPOT Make a discrete (sparse) potential.
% pot = dpot(domain, sizes, T, spar)
%
% sizes(i) is the size of the i'th domain element.
% T defaults to all 1s.

%assert(length(sizes) == length(domain));

pot.domain = domain(:)'; % so we can see it when we display
if nargin < 3
  pot.T = myones(sizes);
  %pot.T = ones(1,prod(sizes)); % 1D vector
else
   if isempty(T)
      pot.T = [];
   else
      if issparse(T)
         pot.T = T;
      else
         pot.T = myreshape(T, sizes);
      end
   end
end
pot.sizes = sizes(:)';
pot = class(pot, 'dpot');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [R Rx Ry Rz] = angles2rotmat(angles)
% [R Rx Ry Rz] = angles2rotmat(angles)
%
% Convert 3 euler angles into a rotation matrix
%
% angles is a 3x1 vector in radians
% angles(1) - pitch - rotation about x or LR (gamma)
% angles(2) - yaw   - rotation about y or AP (beta)
% angles(3) - roll  - rotation about z or SI (alpha)
% R = Rz*Ry*Rx;
%
% See also: rotmat2angles
% Ref: Craig, Intro to Robotics
%
% $Id: angles2rotmat.m,v 1.4 2011/03/02 00:04:12 nicks Exp $

%
% angles2rotmat.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.4 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

R  = [];
Rx = [];
Ry = [];
Rz = [];
if(nargin ~= 1)
  fprintf('R = angles2rotmat(angles)\n');
  return;
end

gamma = angles(1);
beta  = angles(2);
alpha = angles(3);

Rx = zeros(3,3);
Rx(1,1) = +1;
Rx(2,2) = +cos(gamma);
Rx(2,3) = -sin(gamma);
Rx(3,2) = +sin(gamma);
Rx(3,3) = +cos(gamma);

Ry = zeros(3,3);
Ry(1,1) = +cos(beta);
Ry(1,3) = +sin(beta);
Ry(2,2) = +1;
Ry(3,1) = -sin(beta);
Ry(3,3) = +cos(beta);

Rz = zeros(3,3);
Rz(1,1) = +cos(alpha);
Rz(1,2) = -sin(alpha);
Rz(2,1) = +sin(alpha);
Rz(2,2) = +cos(alpha);
Rz(3,3) = +1;

R = Rz*Ry*Rx;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [] = outSVMData(f)
    fid = fopen('f.txt', 'w');
	for ii = 1 : size (f, 1)
	    fprintf (fid, '%f ', -1);
	    for jj = 1 : size (f, 2)
		    fprintf (fid, '%d:', jj);
		    fprintf (fid, '%f ', f(ii, jj));
		end
		fprintf (fid, '\n');
	end
	fclose (fid);
return;

author: dsoellinger


text:
function value = hermite_integral ( n )

%*****************************************************************************80
%
%% HERMITE_INTEGRAL returns the value of a Hermite polynomial integral.
%
%  Discussion:
%
%    H(n) = Integral ( -oo < x < +oo ) x^n exp(-x^2) dx
%
%    H(n) is 0 for n odd.
%
%    H(n) = (n-1)!! * sqrt(pi) / 2^(n/2) for n even.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the integral.
%    0 <= N.
%
%    Output, real VALUE, the value of the integral.
%
  if ( n < 0 )

    value = - r8_huge ( );

  elseif ( mod ( n, 2 ) == 1 )

    value = 0.0;

  else

    value = i4_factorial2 ( n - 1 ) * sqrt ( pi ) / 2.0^( n / 2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [post nlZ dnlZ] = infKL_sprox(hyp, mean, cov, lik, x, y)
% PG-SVI

if hyp.is_cached==1
	global cache_post;
	global cache_nlz;
	global cache_idx;

	post=cache_post(cache_idx);
	nlZ=cache_nlz(cache_idx);
	if nargout>2
		warning('to be implemented\n');
		dnlZ = NaN;
	end
	return
end

snu2=hyp.snu2;
n=size(x,1);

% GP prior
K = feval(cov{:}, hyp.cov, x);                  % evaluate the covariance matrix
m = feval(mean{:}, hyp.mean, x);                      % evaluate the mean vector
K=snu2*eye(n)+K;

lik_name = func2str(lik{1});
mini_batch_size=hyp.mini_batch_size;
assert (mini_batch_size>0)
mini_batch_num=ceil(n/mini_batch_size);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%init value
post_m=hyp.init_m;%k=0
tW = zeros(n,1);%k=-1
post_v=diag(hyp.init_V);%k=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

iter = 0;%iteration
pass=0;%pass
max_pass=hyp.max_pass;
beta = hyp.learning_rate;
r = 1/(beta+1);
index=1:n;
while pass<max_pass
	if mini_batch_size<n
		index=randperm(n);
	end
	offset=0;
	mini_batch_counter=0;
	pass=pass+1;
	while mini_batch_counter<mini_batch_num
		mini_batch_counter=mini_batch_counter+1;
		iter=iter+1;
		%mini batch
		tmp_idx = mini_batch_counter*mini_batch_size;
		idx=index( (tmp_idx-mini_batch_size+1):min(tmp_idx,n) );

		weight=double(n)/size(x(idx,:),1);

		if hyp.stochastic_approx==1
			[ll, gf, gv] = sampling_E(y(idx), post_m(idx), post_v(idx), lik, hyp.sample_size, hyp.lik);
		else
			switch lik_name
			case {'laplace','likLaplace','poisson','bernoulli_logit','likLogistic'}
				[ll, gf, gv] = E_log_p(lik_name, y(idx), post_m(idx), post_v(idx), hyp.lik);
			otherwise
				[ll,gf,d2f,gv] = likKL(post_v(idx), lik, hyp.lik, y(idx), post_m(idx));
			end
		end

		% pseudo observation
		pseudo_y = m + K(:,idx)*(weight*gf) - post_m;
		tW = r.*tW;
		tW(idx) = tW(idx)+(1-r).*((-2*weight)*gv);%tW^{k}, where W=-2*gv
		sW = sqrt(abs(tW)) .* sign(tW);
		L = chol(eye(n)+sW*sW'.*K); %L = chol(sW*K*sW + eye(n));

		%use this following line if we approximate r^{k} .* tW.^{k-1} by tW.{k}
		post_m = post_m + (1-r).*(pseudo_y - K*(sW.*(L\(L'\(sW.*pseudo_y)))));%m^{k+1}
		T = L'\(repmat(sW,1,n).*K); %T  = L'\(sW*K);
		post_v = diag(K) - sum(T.*T,1)'; % v = diag(inv(inv(K)+diag(W))); %v^{k+1}

		if isfield(hyp,'save_iter') && hyp.save_iter==1
			global cache_nlz_iter
			global cache_iter

			alpha=K\(post_m-m);
			nlZ2=compute_nlz(lik, hyp, sW, K, m, alpha, post_m, y);
			cache_iter=[cache_iter; iter];
			cache_nlz_iter=[cache_nlz_iter; nlZ2];
		end

	end
	alpha=K\(post_m-m);
	nlZ=compute_nlz(lik, hyp, sW, K, m, alpha, post_m, y);

	if isfield(hyp,'save_iter') && hyp.save_iter==1
		if pass==1
			global num_iters_at_pass;
			num_iters_at_pass=iter;
		end
		fprintf('pass:%d) at %d iter %.4f %.4f\n', pass, iter, nlZ, nlZ2);
	else
		fprintf('pass:%d) %.4f\n', pass, nlZ);
	end

	if hyp.is_save==1
		global cache_post;
		global cache_nlz;

		post.sW = sW;                                             % return argument
		post.alpha = alpha;
		post.L = L;

		cache_post=[cache_post; post];
		cache_nlz=[cache_nlz; nlZ];
	end
end

alpha=K\(post_m-m);
post.sW = sW;                                             % return argument
post.alpha = alpha;
post.L = L;                                              % L'*L=B=eye(n)+sW*K*sW

nlZ=compute_nlz(lik, hyp, sW, K, m, alpha, post_m, y);
fprintf('final: %.4f\n', nlZ);

if nargout>2
  warning('to be implemented\n');
  dnlZ = NaN;
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i = 1 : size(X, 1)
    minDist = 10000000;
    minDistIdx = 0;
    for j = 1 : K
        dist = sqrt(sum(power((X(i,:) - centroids(j,:)),2)));
        if dist < minDist
            minDist = dist;
            minDistIdx = j;
        end
    end
    idx(i,1) = minDistIdx;
end

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

end

author: UtkarshPathrabe


text:
function [stat,msg] = myweb_nocheck(html_file)

if ismac,
  [stat,msg] = unix(['open ' html_file]);
elseif isunix,
  browsers = {'xdg-open','google-chrome','firefox','opera','konqueror','epiphany'};
  browser = '';
  for i = 1:numel(browsers),
    [a,b] = unix(sprintf('which %s',browsers{i}));
    b = strtrim(b);
    if a == 0 && ~isempty(b),
      browser = b;
    end
  end
  if isempty(browser),
    stat = 1;
    msg = sprintf('Could not find a web browser');
    return
  end
  [stat,msg] = unix(sprintf('%s %s',browser,html_file));
elseif ispc,
  html_file = fullfile(html_file);
  [stat,msg] = dos(['cmd.exe /c rundll32 url.dll,FileProtocolHandler "' html_file '"']);
end

author: kristinbranson

repo: JAABA


text:
function [trl, event] = ft_trialfun_brainvision_segmented(cfg)

% FT_TRIALFUN_BRAINVISION_SEGMENTED creates trials for a Brain Vision Analyzer
% dataset that was segmented in the BVA software.
%
% Use as
%   cfg          = [];
%   cfg.dataset  = 'filename.vhdr';
%   cfg.trialfun = 'ft_trialfun_brainvision_segmented';
%   cfg  = ft_definetrial(cfg);
%   data = ft_preprocessing(cfg);
%
% Optionally, you can specify:
%   cfg.stimformat = 'S %d'
%
% The stimformat instruct this function to parse stimulus triggers according to
% the specific format. The default is 'S %d'. The cfg.stimformat always
% needs to contain exactly one %d code. The trigger values parsed in this way
% will be stored in columns 4 and upwards of the output 'trl' matrix, and
% after FT_PREPROCESSING will end up in data.trialinfo.
%
% A BrainVision dataset consists of three files: an .eeg, .vhdr, and a .vmrk
% file. The option cfg.dataset should refer to the .vhdr file.
%
% See also FT_DEFINETRIAL, FT_PREPROCESSING

% Copyright (C) 2014, Robert Oostenveld, DCCN
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

hdr = ft_read_header(cfg.dataset);
event = ft_read_event(cfg.dataset);

% for backward compatibility
cfg = ft_checkconfig(cfg, 'renamed', {'trigformat', 'stimformat'});

% set the defaults
cfg.stimformat = ft_getopt(cfg, 'stimformat', 'S %d');

sel = find(strcmp({event.type}, 'New Segment'));
ntrial = length(sel);
fprintf('found %d segments in the data\n', ntrial);

sample = [event.sample];
if any(diff(sample)<0)
  % events should be ordered according to the sample in the file at which they occur
  ft_warning('reordering events based on their sample number');
  [sample, order] = sort(sample);
  event = event(order);
end

% for Brain Vision Analyzer the event type and value are strings
type  = {event.type};
value = {event.value};

begsample = nan(1,ntrial);
endsample = nan(1,ntrial);
offset    = nan(1,ntrial);
stim      = cell(1,ntrial);

for i=1:ntrial
  begevt = sel(i);
  if i<ntrial
    endevt = sel(i+1);
  else
    endevt = length(event);
  end

  % specify at which data sample the trial begins
  begsample(i) = sample(sel(i));

  % find the sample on which the trial is time-locked
  t0 = find(strcmp(type((begevt+1):endevt), 'Time 0'));
  t0 = t0(1); % take the first one in case there are multiple
  offset(i) = -(sample(begevt+t0) - begsample(i));

  % construct a list of all stimulus events in each segment
  stim{i} = find(strcmp(type((begevt+1):endevt), 'Stimulus'))+begevt;
end

% the endsample of each trial aligns with the beginsample of the next one
endsample(1:end-1) = begsample(2:end) - 1;
% the last endsample corresponds to the end of the file
endsample(end)     = hdr.nSamples*hdr.nTrials;

% add the stimulus events to the output, if possible
numstim = cellfun(@length, stim);
if (numstim > 0)
  if all(numstim==numstim(1))
    for i=1:length(stim)
      for j=1:numstim(i)
        if isempty(value{stim{i}(j)})
          ft_error('missing a stimulus type definition in the related *.vmrk file');
        end
        stimvalue  = sscanf(value{stim{i}(j)}, cfg.stimformat);
        stimsample = sample(stim{i}(j));
        stimtime   = (stimsample - begsample(i) + offset(i))/hdr.Fs; % relative to 'Time 0'
        if isempty(stimvalue)
          ft_error('the upper case letter of the stimulus value does not match with definition of "cfg.stimformat"');
        end
        trialinfo(i,2*(j-1)+1) = stimvalue;
        trialinfo(i,2*(j-1)+2) = stimtime;
      end % j
    end % i
  else
    ft_warning('the trials have a varying number of stimuli, not adding them to the "trl" matrix');
    trialinfo = [];
  end
else
  ft_warning('the trials have no stimuli, no "trialinfo" will be added to the "trl" matrix');
  trialinfo = [];
end

% combine the sample information and the trial information
trl = cat(2, [begsample(:) endsample(:) offset(:)], trialinfo);

author: spm


text:
function cost = snd_runningcosts(k, x, u, varargin) %this shows the objective funciton
    u0_ref  = cell2mat(varargin(1));

    if nargin >= 5
        battery = cell2mat(varargin(2));
    end
        %@@UPDATE tou price
        A = battery.lifeParam(1,1);
        b = battery.lifeParam(1,2);
        capacity = battery.capacity;
        totalprice = battery.totalprice;
        coeff = totalprice/(2*A*(capacity^b));
    cost = ( (u(1)-u0_ref(1,k))^2 + (u(2)-u0_ref(2,k))^2 ) ;

    if u(2)*x(1)>=0   % x(1) is the cumulative kWh
        cost = cost + coeff*( abs(x(1)+u(2))^b - abs(x(1)^b ) );
    else %0.0001*(1/u(2))^2 +
        cost = cost + coeff*( abs(u(2))^b );
    end
end

author: juchengquan


text:
function [vol, M, dcminfo, mr_parms] = load_dicom_fl(flist)
% [vol, M, dcminfo, mr_parms] = load_dicom_fl(flist)
%
% Loads a volume from the dicom files in flist.
%
% The volume dimensions are arranged such that the
% readout dimension is first, followed by the phase-encode,
% followed by the slices (this is not implemented yet)
%
% M is the 4x4 vox2ras transform such that
% vol(i1,i2,i3), xyz1 = M*[i1 i2 i3 1] where the
% indicies are 0-based.
%
% mr_parms = [tr flipangle te ti]
%
% Does not handle multiple frames correctly yet.
%

%
% load_dicom_fl.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.11 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

vol=[];
M=[];

if(nargin ~= 1)
  fprintf('[vol, M] = load_dicom_fl(flist)\n');
  return;
end

nfiles = size(flist,1);

tic
fprintf('Loading dicom info foreach file \n');
for n = 1:nfiles
  fname = deblank(flist(n,:));
  % fprintf('n = %d/%d, %s   %g\n',n,nfiles,fname,toc);
  tmpinfo = dicominfo(fname);
  if(isempty(tmpinfo))
    fprintf('ERROR: reading %s\n',fname);
    return;
  end
  if(n > 1)
    % Check that the nth series number agrees with the first
    if(tmpinfo.SeriesNumber ~= dcminfo0(1).SeriesNumber)
      fprintf('ERROR: series number inconsistency (%s)\n',fname);
      return;
    end
  end
  tmpinfo.fname = fname;
  dcminfo0(n) = tmpinfo;
end

% Sort by slice location %
dcminfo = sort_by_sliceloc(dcminfo0);

% Slice direction cosine %
sdc = dcminfo(nfiles).ImagePositionPatient-dcminfo(1).ImagePositionPatient;
sdc = sdc /sqrt(sum(sdc.^2));

% Distance between slices %
dslice = sqrt(sum((dcminfo(2).ImagePositionPatient-dcminfo(1).ImagePositionPatient).^2));

% Matrix of direction cosines %
Mdc = zeros(3,3);
Mdc(:,1) = dcminfo(1).ImageOrientationPatient(1:3);
Mdc(:,2) = dcminfo(1).ImageOrientationPatient(4:6);
Mdc(:,3) = sdc;

% Voxel resolution %
delta = [dcminfo(1).PixelSpacing; dslice];
D = diag(delta);

% XYZ of first voxel in first slice %
P0 = dcminfo(1).ImagePositionPatient;

% Change Siemens to be RAS %
% GE is also LPS - change it to be RAS, too - ebeth %
Manufacturer = dcminfo(1).Manufacturer;
if(strcmpi(Manufacturer,'Siemens') | strcmpi(Manufacturer,'ge medical systems'))
  % Change to RAS
  Mdc(1,:) = -Mdc(1,:);
  Mdc(2,:) = -Mdc(2,:);
  P0(1)    = -P0(1);
  P0(2)    = -P0(2);
end

% Compute vox2ras transform %
M = [Mdc*D P0; 0 0 0 1];
if (0&strcmpi(Manufacturer,'Siemens'))
  % Correcting for P0 being at corner of
  % the first voxel instead of at the center
  M = M*[[eye(3) [0.5 0.5 0]']; 0 0 0 1];  %'
end

% Pre-allocate vol. Note: column and row designations do
% not really mean anything. The "column" is the fastest
% dimension. The "row" is the next fastest, etc.
ndim1 = dcminfo(1).Columns;
ndim2 = dcminfo(1).Rows;
ndim3 = nfiles;
vol = zeros(ndim1,ndim2,ndim3);

fprintf('Loading data from each file.\n');
for n = 1:nfiles
  %fprintf('n = %d, %g\n',n,toc);
  fname = dcminfo(n).fname;
  x = dicomread(fname);
  if(isempty(x))
    fprintf('ERROR: could not load pixel data from %s\n',fname);
    return;
  end
  % Note: dicomread will transposed the image. This is supposed
  % to help. Anyway, to make the vox2ras transform agree with
  % the volume, the image is transposed back.
  vol(:,:,n) = x'; %'
end

% Reorder dimensions so that ReadOut dim is first %
if(0 & ~strcmpi(dcminfo(1).PhaseEncodingDirection,'ROW'))
  % This does not work
  fprintf('INFO: permuting vol so that ReadOut is first dim\n');
  vol = permute(vol,[2 1 3]);
  Mtmp = M;
  M(:,1) = Mtmp(:,2);
  M(:,2) = Mtmp(:,1);
end

% Lines below correct for the Z-offset in GE machines - ebeth %
% We're told ge machines recenter along superior/inferior axis but
% don't update c_ras - but now c_s should be zero.
if(strcmpi(Manufacturer,'ge medical systems'))
  % Lines below correct for the Z-offset in GE machines
  firstZ = dcminfo(1).ImagePositionPatient(3);
  lastXYZ = M*[size(vol)';1]; %'
  % size(imvol) = number of slices in all 3 dirs
  lastZ = lastXYZ(3);
  offsetZ = (lastZ + firstZ)/2.0;
  % Z0 = Z + offsetZ;  [XYZ1]' = M*[CRS1]', need to add to M(3,4)(?)
  M(3,4) = M(3,4) - offsetZ;
end

% if(strcmpi(Manufacturer,'ge medical systems'))
%   M(3,4) = 0; % Wow - that was actually completely wrong!
% end

% Pull out some info from the header %
if(isfield(dcminfo(1),'FlipAngle')) FlipAngle = pi*dcminfo(1).FlipAngle/180;
else FlipAngle = 0;
end
if(isfield(dcminfo(1),'EchoTime')) EchoTime = dcminfo(1).EchoTime;
else EchoTime = 0;
end
if(isfield(dcminfo(1),'RepetitionTime')) RepetitionTime = dcminfo(1).RepetitionTime;
else RepetitionTime = 0;
end
InversionTime = 0;
mr_parms = [RepetitionTime FlipAngle EchoTime InversionTime];

return;

%----------------------------------------------------------%
function dcminfo2 = sort_by_sliceloc(dcminfo)
  nslices = length(dcminfo);
  sliceloc = zeros(nslices,1);
  for n = 1:nslices
    sliceloc(n) = dcminfo(n).SliceLocation;
  end

  [tmp ind] = sort(sliceloc);
  dcminfo2 = dcminfo(ind);
return;

%----------------------------------------------------------%


text:
function xyPoints=getEllipseHullPoints(z,A,gammaVal,numPoints,invertA)
%%GETELLIPSHULLPOINTS  Consider an ellipsoid in 3D (x,y,z) coordinates. For
%       every fixed z coordinate, if the x-y plane at the z coordinate
%       intersects the ellipsoid, it cuts an ellipse. For example, the
%       function projEllipse2ZPlane gets the parameters of the ellipse of
%       intersection. This function gets points marking the outer limits of
%       all of those ellipses in the x-y plane (the hull of the ellipses).
%
%INPUTS: z A 3XN vector corresponding to the centers of the N ellipses for
%          which points should be obtained.
%        A A 3X3XN set of N positive definite matrices that specify the
%          size and shape of the ellipsoids, where a point zp is on the ith
%          ellipsoid if
%          (zp-z(:,i))'*A(:,:,i)*(zp-z(:,i))=gammaVal (if invertA is true,
%          then replace A(:,:,i) with inv(A(:,:,i))).
% gammaVal An optional parameter specifying the size of the ellipse/
%          ellipsoid. If omitted or an empty matrix is passed, then
%          gammaVal=18.8049 is used. This is approximately the value for a
%          99.97% confidence region if A are inverse covariance matrices of
%          a Gaussian distribution. gammaVal must be positive gammaVal must
%          be positive.
% numPoints An optional parameter specifying how many points should be
%          generated. The default if omitted or an empty matrix is passed
%          is 500.
%  invertA If this is true, then A is inverted before use. The default if
%          omitted or an empty matrix is passed is false.
%
%OUTPUTS: xyPoints A 2XnumPointsXN set of the points for each of the
%                  ellipsoid projection hulls.
%
%The 3D ellipsoid (with the coordinate system shited to put the origin at
%the center) is defined as x'*R*x=gammaVal. We can rewrite this as
%[r*u,z]*R*[r*u;z]=gammaVal
%where r is a positive scalar and u is a 2X1 unit vector. The matrix R can
%be broken up accordingly as
%R=[Rxy, rz;
%   rz', rzz];
%We want the maximum extent of the ellipsoid in the x-y plane for each
%direction (as specified by u). This will provide the hull of the ellipsoid
%points. Thus, we perform the optimization:
% maximize (over r,z) r^2
% such that r^2*u'*Rxy*u+2*r*z*rz'*u+z^2*rzz=gammaVal
%where the constraint is just from rewriting the expression for the
%definition of an ellipsoid. Writing the Lagrangian, (with lambda as the
%Lagrangian parameter) taking derivatives and setting them equal to 0, we
%get two additional equations. Solving those two equations with the
%original constraint, we get expressions for r, z, and lambda. Choosing the
%solution such that r is positive and nonzero, one gets the expressions
%implemented in this function.
%
%December 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

if(nargin<4||isempty(numPoints))
    numPoints=500;
end

if(nargin<3||isempty(gammaVal))
    gammaVal=18.8049;
end

if(invertA)
    A=applyFunToEachMatrix(@inv,A);
end

numEllipse=size(z,2);

theta=linspace(-pi,pi,numPoints);
u=[cos(theta);
   sin(theta)];
xyPoints=zeros(2,numPoints,numEllipse);
for curEllipse=1:numEllipse
    Rxyrzz=A(1:2,1:2,curEllipse)*A(3,3,curEllipse);
    rz=A(1:2,3,curEllipse);
    rzzGamma=A(3,3,curEllipse)*gammaVal;
    zXYCur=z(1:2,curEllipse);

    for k=1:numPoints
        uCur=u(:,k);

        r=sqrt(rzzGamma/(uCur'*Rxyrzz*uCur-(rz'*uCur)^2));
        xyPoints(:,k,curEllipse)=zXYCur-r*uCur;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Y, X] = svmlread(fname)
% SVMLREAD - Read a data file generated by SVM light
%
%   Y = SVMLREAD(FNAME)
%   FNAME gives the name of an output file generated by SVM light. It
%   may contain predicted labels, coefficients alpha, or an input
%   (example) file with class values and features. From this file the
%   data in the first column (class labels or alphas) is extrated and
%   returned in Y.
%   [Y, X] = SVMLREAD(FNAME), where FNAME is the name of an input file
%   with class values and features, returns both the vector of class
%   labels Y and the matrix of examples X. Each line of X corresponds to
%   a line in the file.
%   Attention: this may take a while...
%
%   See also SVML, SVM_LEARN, SVM_CLASSIFY, SVMLOPT, SVMLWRITE
%

%
% Copyright (c) by Anton Schwaighofer (2001)
% $Revision: 1.6 $ $Date: 2002/02/19 12:26:07 $
% mailto:anton.schwaighofer@gmx.net
%
% This program is released unter the GNU General Public License.
%

error(nargchk(1, 1, nargin));

X = [];
Y = [];

f = fopen(fname, 'rt');
if (f<0),
  error(sprintf('Unable to open file %s', fname));
end

i = 0;
% fprintf('Scanning ');
while ~feof(f),
  s = fgetl(f);
  [Yi, count, errmsg, nextind] = sscanf(s, '%f', 1);
  % read the class label resp. anything else that is in the first column
  if (count==1),
    i = i+1;
    Y(i,1) = Yi;
    [Xi, count] = sscanf(s(nextind:end), ' %i:%f');
    % scan for the feature:value pairs
    if (rem(count,2)==0) & (count~=0),
      % if they really come in pairs, then accept
      ind = 2:2:count;
      if isempty(X),
        maxCol = max(Xi(ind-1));
        approxSparsity = (count/2)/maxCol;
        % a rough estimate of the sparsity, based on the first line of
        % data
        if approxSparsity>0.5,
          approxSparsity = 1;
          X = zeros(maxCol, 1000);
        else
          X = spalloc(maxCol, 1000, round(1000*maxCol*approxSparsity));
          % allocate for 1000 data points (lines) beforehand
          % We store everything *columnwise* and transpose afterwards,
          % this greatly improves performance
        end
      end
      X(Xi(ind-1),i) = Xi(ind);
%       if (rem(i,100)==0),
%         fprintf(' %i', i);
%       end
    end
  end
end
% fprintf(' done.\n');
if ~isempty(X),
  X = X(:,1:i)';
  sparsity = length(find(X))/prod(size(X));
  if sparsity<0.5,
    X = sparse(X);
    % remove any surplus lines & convert to sparse a second time for
    % optimal memory usage
  end
end

fclose(f);

author: viggin

repo: domain-adaptation-toolbox


text:
function ppm = getPPM(MRSIStruct)
    ppm = MRSIStruct.ppm;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [ relpos, psize3, delta_midy ] = gen2dOffset_v2( box1, box2, issplitting )
%GEN2DOFFSET_V2 compute the rotation and offsets between two sibling boxes
% INPUT: box1, box2 - 3d box of base box with absolute position
%
% OUTPUT: relpos - 13d vector (rp of box2 relative to box1), including:
%                    (1) relative orientation - 1 bit degree
%                    (3) offsets - 2 bits
%                    (4) 16 classes (4 classes along each axis) - 8 bits
%                    (5) attachments (4 possible cases) - 4 bits
%                    (6) alignments (aligned along each axis) - 2 bit

% box1 info
pcent = box1(1:3);
pfront = box1(4:6);
pup = box1(7:9);
psize = box1(10:12);
paxes = cross(pfront,pup);
paxes = paxes/norm(paxes,2);

% box2 info
ccent = box2(1:3);
cfront = box2(4:6);
cup = box2(7:9);
csize = box2(10:12);
caxes = cross(cfront,cup);
caxes = caxes/norm(caxes,2);

transform = genTransMat(pfront,pup);
ncent = ccent-pcent;
ncent = transform*ncent; % child's center position in parent's local frame
nfront = cfront;
nfront = transform*nfront; % child's front direction in parent's local frame
nfront = nfront/norm(nfront,2);
nup = cup;
nup = transform*nup; % child's up direction in parent's local frame
nup = nup/norm(nup,2);
naxes = cross(nfront,nup);
naxes = naxes/norm(naxes,2);

% 4 corner points of parent box under its own local frame
plfront = [1;0;0];
plaxes = [0;0;1];
ppoints(:,1) = psize(1)*plfront/2 + psize(3)*plaxes/2;
ppoints(:,2) = psize(1)*plfront/2 - psize(3)*plaxes/2;
ppoints(:,3) = -psize(1)*plfront/2 - psize(3)*plaxes/2;
ppoints(:,4) = -psize(1)*plfront/2 + psize(3)*plaxes/2;

% 4 corner points of child box under its parent's local frame
cpoints(:,1) = ncent + csize(1)*nfront/2 + csize(3)*naxes/2;
cpoints(:,2) = ncent + csize(1)*nfront/2 - csize(3)*naxes/2;
cpoints(:,3) = ncent - csize(1)*nfront/2 - csize(3)*naxes/2;
cpoints(:,4) = ncent - csize(1)*nfront/2 + csize(3)*naxes/2;

% the relative positions
pmaxx = max(ppoints(1,:));
pmaxy = max(ppoints(3,:));
pminx = min(ppoints(1,:));
pminy = min(ppoints(3,:));
cmaxx = max(cpoints(1,:));
cmaxy = max(cpoints(3,:));
cminx = min(cpoints(1,:));
cminy = min(cpoints(3,:));

delta_midy = 0;
if(issplitting==1)
    pmidy = (pmaxy+pminy)/2;
    cmidy = (cmaxy+cminy)/2;
    delta_midy = pmidy-cmidy;
    pmaxy = cmaxy;
    pminy = cminy;
end
psize3 = pmaxy-pminy;

% relative position class along x axis
delta_x = [pminx-cmaxx;pminx-cminx;pmaxx-cmaxx;pmaxx-cminx];
[~,classxindex] = min(abs(delta_x));
offsets(1) = delta_x(classxindex);

% relative position class along y axis
delta_y = [pminy-cmaxy;pminy-cminy;pmaxy-cmaxy;pmaxy-cminy];
[~,classyindex] = min(abs(delta_y));
offsets(2) = delta_y(classyindex);

% alignment info
attachindex = 1;
attacheps = 0.011; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(abs(offsets(1))<attacheps)
    attachindex = attachindex + 1;
end
if(abs(offsets(2))<attacheps)
    attachindex = attachindex + 2;
end
attachvec = zeros(4,1);
attachvec(attachindex) = 1;

% form the relpos vector
rot_eps = 0.01;
alignment = [1;0;0;0;0];
if(abs(nfront(1)-1)<rot_eps)
    alignment = [0;1;0;0;0]; % nfront = [1,0,0];
end
if(abs(nfront(1)+1)<rot_eps)
    alignment = [0;0;1;0;0]; % nfront = [-1,0,0];
end
if(abs(nfront(3)-1)<rot_eps)
    alignment = [0;0;0;1;0]; % nfront = [0,0,1];
end
if(abs(nfront(3)+1)<rot_eps)
    alignment = [0;0;0;0;1]; % nfront = [0,0,-1];
end

if(nfront(1)==0)
    if(nfront(3)>0)
        degree = -0.5;
    else
        degree = 0.5;
    end
else
    degree = atand(nfront(3)/nfront(1));
    if(nfront(1)<0)
        if(degree<0)
            degree = degree-180;
        else
            degree = degree+180;
        end
    end
    degree = mod(degree,360)/180-1;
end
if(abs(degree)>1)
        o=0;
    end

classindex = classxindex + (classyindex-1)*4;
classvec = zeros(16,1);
classvec(classindex) = 1;
relpos = [degree;offsets(1);offsets(2);classvec;attachvec;alignment];

% tmpcent = pcent;
% tmpfsize = (max(cmaxx,pmaxx)-min(cminx,pminx));
% tmpasize = (max(cmaxy,pmaxy)-min(cminy,pminy));
% tmpobb = [tmpcent;plfront;0;1;0;tmpfsize;1;tmpasize];
% draw3dOBB_v2(box1,'r');
% draw3dOBB_v2(box2,'g');
% draw3dOBB_v2(tmpobb,'y');

end

repo: GRAINS


text:
function out = input2(text,def)
    %
    %	out = input2(text,default)
    %
    %	Compliments of Charles Stark Draper Labs
    if ischar(def) == 1
        eval(['y = input('' ',text,' [',def,'] '',''s'');'])
        if isempty(y)
            out = def;
        else
            out = y;
        end
    else
        eval(['y = input('' ',text,' [',num2str(def),'] '',''s'');'])
        if isempty(y)
            out = def;
        else
            out = eval(y);
        end
    end

author: CelsoReyes


text:
function obj = modelPreprocessMovie(obj)
	% Controller for pre-processing movies, mainly aimed at calcium imaging data. Wrapper for modelPreprocessMovieFunction.
	% Biafra Ahanonu
	% started 2013.11.09 [10:46:23]

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% verify that there are folders present
	if isempty(obj.inputFolders)
		uiwait(ciapkg.overloaded.msgbox('Please run modelAddNewFolders to add folders to the object','Success','modal'));
		return
	end
	% check that Miji exists, if not, have user enter information
	% try
	% 	Miji
	% 	MIJ.exit
	% catch
	modelAddOutsideDependencies('miji');
	if obj.guiEnabled==1
		if isempty(obj.foldersToAnalyze)
			scnsize = get(0,'ScreenSize');
			[fileIdxArray, ok] = listdlg('ListString',obj.fileIDNameArray,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','which folders to analyze?');
		else
			fileIdxArray = obj.foldersToAnalyze;
		end
	else
		if isempty(obj.foldersToAnalyze)
			fileIdxArray = 1:length(obj.fileIDNameArray);
		else
			fileIdxArray = obj.foldersToAnalyze;
		end
	end

	options.fileFilterRegexp = 'concat_.*.h5';
	folderListInfo = {obj.inputFolders{fileIdxArray}};
	options.datasetName = obj.inputDatasetName;

	% controllerPreprocessMovie2('folderListPath',folderListInfo,'fileFilterRegexp',options.fileFilterRegexp,'datasetName',options.datasetName,'frameList',[]);
	obj.modelPreprocessMovieFunction('folderListPath',folderListInfo,'fileFilterRegexp',options.fileFilterRegexp,'datasetName',options.datasetName,'frameList',[]);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function varargout = drawBox(box, varargin)
%DRAWBOX Draw a box defined by coordinate extents
%
%   drawBox(BOX)
%   Draws a box defined by its extent: BOX = [XMIN XMAX YMIN YMAX].
%
%   drawBox(..., NAME, VALUE)
%   Specifies drawing parameters using parameter name and value. See plot
%   function for syntax.
%
%   drawBox(AX, ...)
%   Specifies the handle of the axis to draw on.
%
%   Example
%     % define some points, compute their box, display everything
%     points = [10 30; 20 50; 20 20; 30 10;40 30;50 20];
%     box = pointSetBounds(points);
%     figure; hold on;
%     drawPoint(points, 's');
%     drawBox(box);
%     axis([0 60 0 60]);
%
%   See Also:
%   drawOrientedBox, drawRect
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/12/2003.
%

%   HISTORY
%   2010-02-22 creation
%   2011-04-01 add support for drawing option, fix bug for several boxes
%   2011-10-11 add management of axes handle

% extract handle of axis to draw on
if isAxisHandle(box)
    ax = box;
    box = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% default values
xmin = box(:,1);
xmax = box(:,2);
ymin = box(:,3);
ymax = box(:,4);

nBoxes = size(box, 1);
r = zeros(nBoxes, 1);

% iterate on boxes
for i = 1:nBoxes
    % exract min and max values
    tx(1) = xmin(i);
    ty(1) = ymin(i);
    tx(2) = xmax(i);
    ty(2) = ymin(i);
    tx(3) = xmax(i);
    ty(3) = ymax(i);
    tx(4) = xmin(i);
    ty(4) = ymax(i);
    tx(5) = xmin(i);
    ty(5) = ymin(i);

    % display polygon
    r(i) = plot(ax, tx, ty, varargin{:});
end

% format output
if nargout > 0
    varargout = {r};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function wt = wfbtremove(d,kk,wt,varargin)
%WFBTREMOVE Remove node(s) from the filterbank tree
%   Usage:  wt = wbftremove(d,kk,wt);
%           wt = wfbtremove(d,kk,wt,'force');
%
%   Input parameters:
%           d   : Level in the tree (0 - root).
%           kk  : Index of the node at level *d* (starting at 0) or array
%                 of indexes.
%           wt  : Wavelet filterbank tree structure (as returned from
%                 |wfbtinit|).
%
%   Output parameters:
%           wt : Modified filterbank structure.
%
%   `wfbtremove(d,kk,wt)` removes existing node at level *d* and index *kk*
%   from the filterbank tree structure *wt*. The function fails if the
%   node has any children (it is not a leaf node).
%
%   `wfbtremove(d,k,wt,'force')` does the same, but any childern of the
%   node are removed too.
%
%   Examples:
%   ---------
%
%   The following example shows magnitude frequency responses of filterbank
%   tree before and after prunning.:::
%
%      % Create a full filterbank tree usinf 'db10' basic filterbank.
%      wt1 = wfbtinit({'db10',4,'full'});
%      % Remove a subtree starting by root's high-pass filter. Force flag
%      % is used because we are removing a non-leaf node.
%      wt2 = wfbtremove(1,1,wt1,'force');
%
%      % Create identical filterbanks
%      [g1,a1] = wfbt2filterbank(wt1,'freq');
%      [g2,a2] = wfbt2filterbank(wt2,'freq');
%
%      % Plot the frequency responses
%      subplot(2,1,1);
%      filterbankfreqz(g1,a1,1024,'plot','posfreq','linabs');
%      subplot(2,1,2);
%      filterbankfreqz(g2,a2,1024,'plot','posfreq','linabs');
%
%

% AUTHOR: Zdenek Prusa

complainif_notenoughargs(nargin,3,'WFBTREMOVE');

definput.flags.force = {'noforce','force'};
flags=ltfatarghelper({},definput,varargin);

if isempty(wt.nodes)
   error('%s: Tree is empty.',mfilename);
end

for k=kk
   [nodeNo,nodeChildIdx] = depthIndex2NodeNo(d,k,wt);
   if(nodeNo==0)
       % removing root
       rootNo = find(wt.parents==0);
       % check for any children of the root
       if any(wt.children{rootNo}~=0) && ~flags.do_force
            error(['%s: Deleting root node. To delete the whole tree ',...
                   'use FORCE option.'],mfilename,d,k);
       else
           wt = nodeSubtreeDelete(rootNo,wt);
           continue;
       end
   end

   % check if node exists
   childrenIdx = find(wt.children{nodeNo}~=0);
   found = find(childrenIdx==nodeChildIdx,1);

   if(isempty(found))
        error('%s: Such node (depth=%d, idx=%d) does not exist.',mfilename,d,k);
   end

   nodeToDelete = wt.children{nodeNo}(nodeChildIdx);
   % Check if it is a leaf (terminal node)
   if any(wt.children{nodeToDelete}~=0) && ~flags.do_force
       error(['%s: Deleting a non-leaf node. To delete whole subtree use ',...
              'FORCE option.'],mfilename);
   else
       wt = nodeSubtreeDelete(nodeToDelete,wt);
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%INTLAB interval hessian toolbox
%
%Hessian constructors
%  hessianinit  - Initialization of dependent variables
%  horzcat      - Horizontal concatenation          [ , ]
%  vertcat      - Vertical concatenation            [ ; ]
%  subsasgn     - Subscripted assignment A(i,:) = 1
%  subsref      - Subscripted reference r = A(3,4)
%  hessian      - Double to gradient
%
%Display of hessians and interval hessians (rigorous)
%  display      - Command window display of hessian
%  disp         - Display function for pop-up windows in debugger
%  realimag     - Real and Imaginary part separately
%  infsup       - Display infimum and supremum
%  midrad       - Display midpoint and radius
%  disp_        - Display in "_" notation
%
%Hessian arithmetic operations
%  plus         - Plus                              +
%  uplus        - Unary plus                        +
%  minus        - Minus                             -
%  uminus       - Unary minus                       -
%  mtimes       - Matrix multiply                   *
%  times        - Elementwise multiply              .*
%  mrdivide     - Slash or right division           /
%  mldivide     - Backslash or left division        \
%  rdivide      - Elementwise right division        ./
%  ldivide      - Elementwise left division         .\
%  mpower       - Matrix power                      ^
%  power        - Elementwise power                 .^
%  intersect    - Intersection
%
%Other hessian operations
%  abs          - Absolute value
%  inf          - Infimum
%  sup          - Supremum
%  mid          - Midpoint
%  rad          - Radius
%  diam         - Diameter
%  real         - Real part
%  imag         - Imaginary part
%  trace        - Trace
%  sum          - Sum
%  prod         - Product
%  ctranspose   - Complex conjugate transpose       '
%  transpose    - Transpose                         .'
%
%Utility routines
%  isnan        - True for Not a Number
%  isreal       - Hessian is real
%  isintval     - Hessian is intval
%  isfinite     - Interval is finite
%  isinf        - Interval is infinite
%  isempty      - Hessian is empty in Matlab sense, i.e. []
%  emptyintersect - detect empty intersections
%  issparse     - Hessian has sparse structure
%  find         - find indices of nonzero elements
%  all          - Determine if all array elements are nonzero
%  any          - Determine if any array elements are nonzero
%  logical      - Convert hessian values to logical
%  nnz          - number of nonzero elements
%  end          - determine last index
%
%Structural operations
%  full         - Convert to full hessian
%  sparse       - Convert to sparse hessian
%  band         - Extract band
%  diag         - Extract diagonal
%  tril         - Extract lower triangular
%  triu         - Extract upper triangular
%  bandwidth    - Bandwidth
%  length       - Length
%  size         - Size
%  dim          - Dimension of square matrix
%  reshape      - Reshape
%  repmat       - Duplicate arrays
%
%Hessian trigonometric functions (rigorous, real and complex)
%  sin          - Sine
%  cos          - Cosine
%  tan          - Tangent
%  cot          - Cotangent
%  sec          - Secant
%  csc          - Cosecant
%  asin         - Inverse sine
%  acos         - Inverse cosine
%  atan         - Inverse tangent
%  acot         - Inverse cotangent
%  asec         - Inverse secant
%  acsc         - Inverse cosecant
%  sinh         - Hyperbolic sine
%  cosh         - Hyperbolic cosine
%  tanh         - Hyperbolic tangent
%  coth         - Hyperbolic cotangent
%  asinh        - Inverse hyperbolic sine
%  acosh        - Inverse hyperbolic cosine
%  atanh        - Inverse hyperbolic tangent
%  acoth        - Inverse hyperbolic cotangent
%
%Hessian exponential functions (rigorous, real and complex)
%  exp          - Exponential
%  log          - Natural logarithm
%  log10        - Logarithm to base 10
%  sqr          - Square
%  sqrt         - Square root
%
%Hessian comparison of ".x" part
%  eq           - Equal                             ==
%  ne           - Not equal                         ~=
%  gt           - Greater than                      >
%  ge           - Greater than or equal             >=
%  lt           - Less than                         <
%  le           - Less than or equal                <=
%
%Verification routines and auxiliary
%  typeadj      - Type adjustment
%  typeof       - Type for type adjustment
%
%Initialization of INTLAB hessian package and system variables
%  hessianinit  - Initialization and definition of INTLAB switches, also
%                   initialization of dependent variables
%
%Demonstration, samples
%  demohessian - Some examples for using INTLAB hessian package
%
%
%The hessian package uses forward mode of automatic differentiation. For
%an introduction to forward differentation, see
%  L.B. Rall: Automatic Differentiation: Techniques and Applications,
%    Lecture Notes in Computer Science 120, Springer, 1981.
%

% written  04/04/04     S.M. Rump  INTLAB Version 5
%
% Copyright (c) Siegfried M. Rump, head of the Institute for Reliable Computing,
%               Hamburg University of Technology

repo: OpenMAS


text:
function [x, cost, info, options] = conjugategradient(problem, x, options)
% Conjugate gradient minimization algorithm for Manopt.
%
% function [x, cost, info, options] = conjugategradient(problem)
% function [x, cost, info, options] = conjugategradient(problem, x0)
% function [x, cost, info, options] = conjugategradient(problem, x0, options)
% function [x, cost, info, options] = conjugategradient(problem, [], options)
%
% Apply the conjugate gradient minimization algorithm to the problem
% defined in the problem structure, starting at x0 if it is provided
% (otherwise, at a random point on the manifold). To specify options whilst
% not specifying an initial guess, give x0 as [] (the empty matrix).
%
% In most of the examples bundled with the toolbox (see link below), the
% solver can be replaced by the present one if need be.
%
% The outputs x and cost are the best reached point on the manifold and its
% cost. The struct-array info contains information about the iterations:
%   iter : the iteration number (0 for the initial guess)
%   cost : cost value
%   time : elapsed time in seconds
%   gradnorm : Riemannian norm of the gradient
%   stepsize : norm of the last tangent vector retracted
%   beta : value of the beta parameter (see options.beta_type)
%   linesearch : information logged by options.linesearch
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below this.
%   maxiter (1000)
%       The algorithm terminates if maxiter iterations have been executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%   minstepsize (1e-10)
%       The algorithm terminates if the linesearch returns a displacement
%       vector (to be retracted) smaller in norm than this value.
%   beta_type ('H-S')
%       Conjugate gradient beta rule used to construct the new search
%       direction, based on a linear combination of the previous search
%       direction and the new (preconditioned) gradient. Possible values
%       for this parameter are:
%           'S-D', 'steep' for beta = 0 (preconditioned steepest descent)
%           'F-R' for Fletcher-Reeves's rule
%           'P-R' for Polak-Ribiere's modified rule
%           'H-S' for Hestenes-Stiefel's modified rule
%           'H-Z' for Hager-Zhang's modified rule
%       See Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods" for a description of these rules in the Euclidean case and
%       for an explanation of how to adapt them to the preconditioned case.
%       The adaption to the Riemannian case is straightforward: see in code
%       for details. Modified rules take the max between 0 and the computed
%       beta value, which provides automatic restart, except for H-Z which
%       uses a different modification.
%   orth_value (Inf)
%       Following Powell's restart strategy (Math. prog. 1977), restart CG
%       (that is, make a -preconditioned- gradient step) if two successive
%       -preconditioned- gradients are "too" parallel. See for example
%       Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods", page 12. An infinite value disables this strategy. See in
%       code formula for the specific criterion used.
%   linesearch (@linesearch_adaptive or @linesearch_hint)
%       Function handle to a line search function. The options structure is
%       passed to the line search too, so you can pass it parameters. See
%       each line search's documentation for info. Another available line
%       search in manopt is @linesearch, in /manopt/linesearch/linesearch.m
%       If the problem structure includes a line search hint, then the
%       default line search used in @linesearch_hint.
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (3)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent.
%   storedepth (2)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb. If the
%       caching features of Manopt are not used, this is irrelevant. For
%       the CG algorithm, a store depth of 2 should always be sufficient.
%
%
% See also: steepestdescent trustregions manopt/solvers/linesearch manopt/examples

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors: Nicolas Boumal
% Change log:
%
%   March 14, 2013, NB:
%       Added preconditioner support : see Section 8 in
%       https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
%
%   Sept. 13, 2013, NB:
%       Now logging beta parameter too.
%
%	Nov. 7, 2013, NB:
%       The search direction is not normalized before it is passed to the
%       linesearch anymore. This way, it is up to the designers of the
%       linesearch to decide whether they want to use the norm of the
%       search direction in their algorithm or not. There are reasons
%       against it, but practical evidence that it may help too, so we
%       allow it. The default linesearch_adaptive used does exploit the
%       norm information. The base linesearch does not. You may select it
%       by setting options.linesearch = @linesearch;
%
%	Nov. 29, 2013, NB:
%       Documentation improved: options are now explicitly described.
%       Removed the Daniel rule for beta: it was not appropriate for
%       preconditioned CG and I could not find a proper reference for it.

% Verify that the problem description is sufficient for the solver.
if ~canGetCost(problem)
    warning('manopt:getCost', ...
        'No cost provided. The algorithm will likely abort.');
end
if ~canGetGradient(problem)
    warning('manopt:getGradient', ...
        'No gradient provided. The algorithm will likely abort.');
end

% Set local defaults here
localdefaults.minstepsize = 1e-10;
localdefaults.maxiter = 1000;
localdefaults.tolgradnorm = 1e-6;
localdefaults.storedepth = 2;
% Changed by NB : H-S has the "auto restart" property.
% See Hager-Zhang 2005/2006 survey about CG methods.
% Well, the auto restart comes from the 'max(0, ...)', not so much from the
% reason stated in Hager-Zhang I believe. P-R also has auto restart.
localdefaults.beta_type = 'H-S';
localdefaults.orth_value = Inf; % by BM as suggested in Nocedal and Wright

% Depending on whether the problem structure specifies a hint for
% line-search algorithms, choose a default line-search that works on
% its own (typical) or that uses the hint.
if ~canGetLinesearch(problem)
    localdefaults.linesearch = @linesearch_adaptive;
else
    localdefaults.linesearch = @linesearch_hint;
end

% Merge global and local defaults, then merge w/ user options, if any.
localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

% for convenience
inner = problem.M.inner;
lincomb = problem.M.lincomb;

% Create a store database
storedb = struct();

timetic = tic();

% If no initial point x is given by the user, generate one at random.
if ~exist('x', 'var') || isempty(x)
    x = problem.M.rand();
end

% Compute objective-related quantities for x
[cost grad storedb] = getCostGrad(problem, x, storedb);
gradnorm = problem.M.norm(x, grad);
[Pgrad storedb] = getPrecon(problem, x, grad, storedb);
gradPgrad = inner(x, grad, Pgrad);

% Iteration counter (at any point, iter is the number of fully executed
% iterations so far)
iter = 0;

% Save stats in a struct array info and preallocate,
% see http://people.csail.mit.edu/jskelly/blog/?x=entry:entry091030-033941
stats = savestats();
info(1) = stats;
info(min(10000, options.maxiter+1)).iter = [];

% Initial linesearch memory
lsmem = [];

if options.verbosity >= 2
    fprintf(' iter\t                cost val\t     grad. norm\n');
end

% Compute a first descent direction (not normalized)
desc_dir = lincomb(x, -1, Pgrad);

% Start iterating until stopping criterion triggers
while true

    % Display iteration information
    if options.verbosity >= 2
        fprintf('%5d\t%+.16e\t%.8e\n', iter, cost, gradnorm);
    end

    % Start timing this iteration
    timetic = tic();

    % Run standard stopping criterion checks
    [stop reason] = stoppingcriterion(problem, x, options, info, iter+1);

    % Run specific stopping criterion check
    if ~stop && abs(stats.stepsize) < options.minstepsize
        stop = true;
        reason = 'Last stepsize smaller than minimum allowed. See options.minstepsize.';
    end

    if stop
        if options.verbosity >= 1
            fprintf([reason '\n']);
        end
        break;
    end

    % The line search algorithms require the directional derivative of the
    % cost at the current point x along the search direction.
    df0 = inner(x, grad, desc_dir);

    % If we didn't get a descent direction: restart, i.e., switch to the
    % negative gradient. Equivalent to resetting the CG direction to a
    % steepest descent step, which discards the past information.
    if df0 >= 0

        % Or we switch to the negative gradient direction.
        if options.verbosity >= 3
            fprintf(['Conjugate gradient info: got an ascent direction '...
                     '(df0 = %2e), reset to the (preconditioned) '...
                     'steepest descent direction.\n'], df0);
        end
        % Reset to negative gradient: this discards the CG memory.
        desc_dir = lincomb(x, -1, Pgrad);
        df0 = -gradPgrad;

    end

    % Execute line search
    [stepsize newx storedb lsmem lsstats] = options.linesearch(...
                 problem, x, desc_dir, cost, df0, options, storedb, lsmem);

    % Compute the new cost-related quantities for x
    [newcost newgrad storedb] = getCostGrad(problem, newx, storedb);
    newgradnorm = problem.M.norm(newx, newgrad);
    [Pnewgrad storedb] = getPrecon(problem, x, newgrad, storedb);
    newgradPnewgrad = inner(newx, newgrad, Pnewgrad);

    % Apply the CG scheme to compute the next search direction.
    %
    % This paper https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
	% by Hager and Zhang lists many known beta rules. The rules defined
    % here can be found in that paper (or are provided with additional
    % references), adapted to the Riemannian setting.
	%
    if strcmpi(options.beta_type, 'steep') || ...
       strcmpi(options.beta_type, 'S-D')              % Gradient Descent

        beta = 0;
        desc_dir = lincomb(x, -1, Pnewgrad);

    else

        oldgrad = problem.M.transp(x, newx, grad);
        orth_grads = inner(newx, oldgrad, Pnewgrad)/newgradPnewgrad;

        % Powell's restart strategy (see page 12 of Hager and Zhang's
        % survey on conjugate gradient methods, for example)
        if abs(orth_grads) >= options.orth_value,
            beta = 0;
            desc_dir = lincomb(x, -1, Pnewgrad);

        else % Compute the CG modification

            desc_dir = problem.M.transp(x, newx, desc_dir);

            if strcmp(options.beta_type, 'F-R')  % Fletcher-Reeves
                beta = newgradPnewgrad / gradPgrad;

            elseif strcmp(options.beta_type, 'P-R')  % Polak-Ribiere+
                % vector grad(new) - transported grad(current)
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                ip_diff = inner(newx, Pnewgrad, diff);
                beta = ip_diff/gradPgrad;
                beta = max(0, beta);

            elseif strcmp(options.beta_type, 'H-S')  % Hestenes-Stiefel+
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                ip_diff = inner(newx, Pnewgrad, diff);
                beta = ip_diff / inner(newx, diff, desc_dir);
                beta = max(0, beta);

            elseif strcmp(options.beta_type, 'H-Z') % Hager-Zhang+
                diff = lincomb(newx, 1, newgrad, -1, oldgrad);
                Poldgrad = problem.M.transp(x, newx, Pgrad);
                Pdiff = lincomb(newx, 1, Pnewgrad, -1, Poldgrad);
                deno = inner(newx, diff, desc_dir);
                numo = inner(newx, diff, Pnewgrad);
                numo = numo - 2*inner(newx, diff, Pdiff)*...
                                       inner(newx, desc_dir, newgrad)/deno;
                beta = numo/deno;

                % Robustness (see Hager-Zhang paper mentioned above)
                desc_dir_norm = problem.M.norm(newx, desc_dir);
                eta_HZ = -1/(desc_dir_norm * min(0.01, gradnorm));
                beta = max(beta,  eta_HZ);

            else
                error(['Unknown options.beta_type. ' ...
                       'Should be steep, S-D, F-R, P-R, H-S or H-Z.']);
            end
            desc_dir = lincomb(newx, -1, Pnewgrad, beta, desc_dir);
        end

    end

    % Make sure we don't use too much memory for the store database.
    storedb = purgeStoredb(storedb, options.storedepth);

    % Update iterate info
    x = newx;
    cost = newcost;
    grad = newgrad;
    Pgrad = Pnewgrad;
    gradnorm = newgradnorm;
    gradPgrad = newgradPnewgrad;

    % iter is the number of iterations we have accomplished.
    iter = iter + 1;

    % Log statistics for freshly executed iteration
    stats = savestats();
    info(iter+1) = stats; %#ok<AGROW>

end

info = info(1:iter+1);

if options.verbosity >= 1
    fprintf('Total time is %f [s] (excludes statsfun)\n', info(end).time);
end

% Routine in charge of collecting the current iteration stats
    function stats = savestats()
        stats.iter = iter;
        stats.cost = cost;
        stats.gradnorm = gradnorm;
        if iter == 0
            stats.stepsize = nan;
            stats.time = toc(timetic);
            stats.linesearch = [];
            stats.beta = 0;
        else
            stats.stepsize = stepsize;
            stats.time = info(iter).time + toc(timetic);
            stats.linesearch = lsstats;
            stats.beta = beta;
        end
        stats = applyStatsfun(problem, x, storedb, options, stats);
    end

end

author: andrewssobral

repo: lrslibrary


text:
function [TorF, vstr, rdate] = have_feature_sdp_pf()
%HAVE_FEATURE_SDP_PF  Detect availability/version info for SDP_PF
%
%   Feature detection function implementing 'sdp_pf' tag for HAVE_FEATURE
%   to detect availability/version of SDP_PF, a MATPOWER extension for
%   applications of semi-definite programming relaxations of power flow
%   equations (https://github.com/MATPOWER/mx-sdp_pf/).
%
%   See also HAVE_FEATURE.

%   MATPOWER
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

TorF = have_feature('yalmip') && exist('mpoption_info_sdp_pf', 'file') == 2;
if TorF
    v = sdp_pf_ver('all');
    vstr = v.Version;
    rdate = v.Date;
else
    vstr = '';
    rdate = '';
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function validatedir(status, mess, messid, filename)
    % check success of directory creation
    if (status == 0)
        if(~strcmp(messid, 'MATLAB:MKDIR:DirectoryExists'))
            error(['Error saving %s folder: ' mess ', ' messid ...
                   '. Please contact support@plot.ly for assistance.'], ...
                   filename);
        end
    end
end

author: plotly


text:
function Animation(ivp, duration, fps, movie)
%--------------------------------------------------------------------------
% To Animate the Bead on a rotating hoop, Equation is solved by MATLAB ode45, and
% running a loop over the time step; for every time step the position
% coordinates of the bead are updated. Phase plane plot and time Vs.
% displacement plot is also plotted
%--------------------------------------------------------------------------
nframes=duration*fps;           % Number of Frames
% Solving the Equation using ode45
sol=ode45(@Equation,[0 duration], ivp);
t = linspace(0,duration,nframes);
y = deval(sol,t);
% Position and Angular velocity of the bead
thita = y(1,:)';
dthita = y(2,:)';
R = ivp(5);                     % Radius of the Hoop
w0 = ivp(7) ;                   % Grequency of rotation of the hoop/ring
%
% Figure starts
fh = figure ;
set(fh,'name','Bead on Rotating Hoop','numbertitle','off','color', 'k','menubar','none') ;
% Plot for the Hoop (Ring/Circle)
phi = linspace(0,2*pi,nframes) ;
c = [ 0 0 0] ;                  % Center of the hoop
u = [1 0 0] ;                   % Normal vectors lying on the plane of circle
v = [0 0 1] ;
px = c(1)+R*cos(phi).*u(1)+R*sin(phi).*v(1) ;
py = c(2)+R*cos(phi).*u(2)+R*sin(phi).*v(2) ;
pz = c(3)+R*cos(phi).*u(3)+R*sin(phi).*v(3) ;
hoop = plot3(px,py,pz,'Color','r','Linewidth',3) ;
% Bead position
beadX = R*sin(thita).*cos(w0*t') ;
beadY = R*sin(thita).*sin(w0*t') ;
beadZ = -R*cos(thita) ;
hold on ;
bead =  plot3(beadX(1),beadY(1),beadZ(1),'MarkerSize',50,'Marker','.','Color','b');
% Trajectory of the bead on rotating hoop
path = plot3(beadX(1),beadY(1),beadZ(1),'Color','w','Linewidth',0.5) ;
% Figure Settings
title('Bead on a Rotating Hoop','Color','w','Fontsize',10);
range = 1.3*R;
campos([0 R 0]) ;
axis([-range range -range range -range range]) ;
axis square ;
axis off ;
rotate3d ;
stop = uicontrol('style','toggle','string','stop','background','white');
%
% Animation starts
for i=1:length(t)-1
 % Animation Plot
    if (ishandle(fh)==1)
        T = [cos(w0*t(i))  -sin(w0*t(i)) 0 ;
             sin(w0*t(i))   cos(w0*t(i)) 0 ;
             0          0       1] ;        % Rotation Matrix
       pn = T*[px;py;pz] ;                  % New Position of the Hoop
       pn = pn' ;
       set(hoop,'XData',pn(:,1),'YData',pn(:,2),'ZData',pn(:,3)) ;
       set(bead,'XData',beadX(i),'YData',beadY(i),'ZData',beadZ(i)) ;
      % drawnow;
        % Trajectory of the Bead
       PP(i,:) = [beadX(i,1) beadY(i,1) beadZ(i,1)];
       set(path,'XData',PP(:,1),'YData',PP(:,2),'ZData',PP(:,3));
            if get(stop,'value')==0
                drawnow;
            elseif get(stop,'value')==1
                break
            end
    end
    F(i) = getframe(fh) ;
    % Pause for length of time step
    if movie == false
            pause(t(i+1)-t(i));
    end
end
if movie == true
        movie2avi(F,'BeadOnRotatingHoop.avi','compression','Cinepak','fps',fps)
end
set(stop,'style','pushbutton','string','close','callback','close(gcf)');

author: Sable

repo: mcbench-benchmarks


function d2img = draw_2Dskeleton_ICVL(tline, coord_pixel)

    db_path = '/home/gyeongsikmoon/workspace/Data/Hand_pose_estimation/ICVL/Testing/Depth/';
    jointNum =  16;
     cubicSz = 200;
    imgWidth = 320;
    imgHeight = 240;

    coord_pixel = squeeze(coord_pixel);
    rgb_img = zeros(imgHeight,imgWidth,3);
    refDepths = zeros(1,jointNum);
    line_width = 4;

    splitted = strsplit(tline);
    img_name = splitted{1};
    bin_name = strcat(db_path,img_name(1:size(img_name,2)-3),'bin');

    for jid = 1:jointNum
        refDepths(1,jid) = str2num(splitted{(jid-1)*3+4});
    end
    refDepth = (min(refDepths(:)) + max(refDepths(:)))/2;

    fp_bin = fopen(bin_name,'r');
    img = fread(fp_bin,[imgWidth imgHeight],'float');
    img = permute(img,[2,1]);
    img(img==0) = refDepth+cubicSz/2;
    fclose(fp_bin);

    img(img>refDepth+cubicSz/2) = refDepth + cubicSz/2;
    img(img<refDepth-cubicSz/2) = refDepth - cubicSz/2;
    img = img - refDepth;
    img = img/(cubicSz/2);

    img = img + 1;
    img = img/2;

    rgb_img(:,:,1) = img*255/255;
    rgb_img(:,:,2) = img*240/255;
    rgb_img(:,:,3) = img*204/255;

    f = figure;
    set(f, 'visible', 'off');
    imshow(rgb_img);
    hold on;

    plot([coord_pixel(1,1),coord_pixel(1,2)],[coord_pixel(2,1),coord_pixel(2,2)],'Color',[255/255,153/255,153/255],'LineWidth',line_width); %wrist to thumb root
    plot([coord_pixel(1,2),coord_pixel(1,3)],[coord_pixel(2,2),coord_pixel(2,3)],'Color',[255/255,102/255,102/255],'LineWidth',line_width) %thumb root to thumb mid
    plot([coord_pixel(1,3),coord_pixel(1,4)],[coord_pixel(2,3),coord_pixel(2,4)],'Color',[255/255,51/255,51/255],'LineWidth',line_width) %thumb mid to thumb tip

    plot([coord_pixel(1,1),coord_pixel(1,5)],[coord_pixel(2,1),coord_pixel(2,5)],'Color',[153/255,255/255,153/255],'LineWidth',line_width) %wrist to index root
    plot([coord_pixel(1,5),coord_pixel(1,6)],[coord_pixel(2,5),coord_pixel(2,6)],'Color',[102/255,255/255,102/255],'LineWidth',line_width) %index root to index mid
    plot([coord_pixel(1,6),coord_pixel(1,7)],[coord_pixel(2,6),coord_pixel(2,7)],'Color',[51/255,255/255,51/255],'LineWidth',line_width) %index mid to index tip

    plot([coord_pixel(1,1),coord_pixel(1,8)],[coord_pixel(2,1),coord_pixel(2,8)],'Color',[255/255,204/255,153/255],'LineWidth',line_width) %wrist to middle root
    plot([coord_pixel(1,8),coord_pixel(1,9)],[coord_pixel(2,8),coord_pixel(2,9)],'Color',[255/255,178/255,102/255],'LineWidth',line_width) %middle root to middle mid
    plot([coord_pixel(1,9),coord_pixel(1,10)],[coord_pixel(2,9),coord_pixel(2,10)],'Color',[255/255,153/255,51/255],'LineWidth',line_width) %middle mid to middle tip

    plot([coord_pixel(1,1),coord_pixel(1,11)],[coord_pixel(2,1),coord_pixel(2,11)],'Color',[153/255,204/255,255/255],'LineWidth',line_width) %wrist to ring root
    plot([coord_pixel(1,11),coord_pixel(1,12)],[coord_pixel(2,11),coord_pixel(2,12)],'Color',[102/255,178/255,255/255],'LineWidth',line_width) %ring root to ring mid
    plot([coord_pixel(1,12),coord_pixel(1,13)],[coord_pixel(2,12),coord_pixel(2,13)],'Color',[51/255,153/255,255/255],'LineWidth',line_width) %ring mid to ring tip

    plot([coord_pixel(1,1),coord_pixel(1,14)],[coord_pixel(2,1),coord_pixel(2,14)],'Color',[255/255,153/255,255/255],'LineWidth',line_width) %wrist to pinky root
    plot([coord_pixel(1,14),coord_pixel(1,15)],[coord_pixel(2,14),coord_pixel(2,15)],'Color',[255/255,102/255,255/255],'LineWidth',line_width) %pinky root to pinky mid
    plot([coord_pixel(1,15),coord_pixel(1,16)],[coord_pixel(2,15),coord_pixel(2,16)],'Color',[255/255,51/255,255/255],'LineWidth',line_width) %pinky mid to pinky tip

    colorList = [
    230/255 230/255 0/255;

    255/255,153/255,153/255;
    255/255 102/255 102/255;
    255/255 51/255 51/255;

    153/255,255/255,153/255;
    102/255,255/255,102/255;
    51/255,255/255,51/255;

    255/255,204/255,153/255;
    255/255 178/255 102/255;
    255/255 153/255 51/255;

    153/255,204/255,255/255;
    102/255 178/255 255/255;
    51/255 153/255 255/255;

    255/255,153/255,255/255;
    255/255,102/255,255/255;
    255/255,51/255,255/255];
    scatter(coord_pixel(1,:),coord_pixel(2,:),100,colorList,'filled');

    set(gca,'Units','normalized','Position',[0 0 1 1]);  %# Modify axes size
    set(gcf,'Units','pixels','Position',[200 200 2*imgWidth 2*imgHeight]);  %# Modify figure size

    frame = getframe(gcf);
    framedata = frame.cdata;

    xmin = min(coord_pixel(1,:));
    xmax = max(coord_pixel(1,:));
    ymin = min(coord_pixel(2,:));
    ymax = max(coord_pixel(2,:));

    len = max(xmax-xmin+1,ymax-ymin+1) + 70;
    xcenter = (xmin + xmax)/2;
    ycenter = (ymin + ymax)/2;

    xmin = max(round(xcenter - len/2),1);
    xmax = min(round(xmin + len),imgWidth);
    ymin = max(round(ycenter - len/2),1);
    ymax = min(round(ymin + len),imgHeight);
    framedata = framedata(2*ymin:2*ymax,2*xmin:2*xmax,:);

    hold off;
    close(f);

    d2img = framedata;
end


text:
function [features] = collect(conf, imgs, scale, filters, verbose)

if nargin < 5
    verbose = 0;
end

num_of_imgs = numel(imgs);
feature_cell = cell(num_of_imgs, 1); % contains images' features
num_of_features = 0;

if verbose
    fprintf('Collecting features from %d image(s) ', num_of_imgs)
end
feature_size = [];

h = [];
for i = 1:num_of_imgs
    h = progress(h, i / num_of_imgs, verbose);
    sz = size(imgs{i});
    if verbose
        fprintf(' [%d x %d]', sz(1), sz(2));
    end

    F = extract(conf, imgs{i}, scale, filters);
    num_of_features = num_of_features + size(F, 2);
    feature_cell{i} = F;

    assert(isempty(feature_size) || feature_size == size(F, 1), ...
        'Inconsistent feature size!')
    feature_size = size(F, 1);
end
if verbose
    fprintf('\nExtracted %d features (size: %d)\n', num_of_features, feature_size);
end
clear imgs % to save memory
features = zeros([feature_size num_of_features], 'single');
offset = 0;
for i = 1:num_of_imgs
    F = feature_cell{i};
    N = size(F, 2); % number of features in current cell
    features(:, (1:N) + offset) = F;
    offset = offset + N;
end

author: thomas-koehler

repo: SupER


text:
function h = explode_tetramesh(T,V,steps)
  % EXPLODE_TETRAMESH animated, interactive tetramesh plot of tets of a tet
  % mesh being pulled away from each
  %
  % h = explode_tetramesh(T,V)
  %
  % Inputs:
  %   T  #T by 4 list of tetrahedra indices
  %   V  #V by dim by #{frames | 1} list of vertex positions
  %   steps  # of explosion animation steps {100}
  % Outputs:
  %   h  handle to plot
  %
  % See also: tsurf, trisurf, tetramesh, animated_trisurf, animated_tetramesh
  %

  % Remap vertices and tets so that all corners are separated
  V = [ ...
    V(T(:,1),:); ...
    V(T(:,2),:); ...
    V(T(:,3),:); ...
    V(T(:,4),:); ...
  ];
  T = [ ...
    0*size(T,1) + (1:size(T,1)); ...
    1*size(T,1) + (1:size(T,1)); ...
    2*size(T,1) + (1:size(T,1)); ...
    3*size(T,1) + (1:size(T,1)); ...
    ]';
  factor = 4;
  if ~exist('steps','var')
    steps = 100;
  end
  s = linspace(1,factor,steps);
  % centroid
  c = mean(V);
  VV = zeros(size(V,1),size(V,2),steps);
  for ii = 1:steps
    sii = s(ii);
    % dialate globally
    Vii = bsxfun(@plus,sii*V,-sii*c+c);
    % contract locally
    % centroids of each tet
    Cii = ( ...
      Vii(T(:,1),:) + ...
      Vii(T(:,2),:) + ...
      Vii(T(:,3),:) + ...
      Vii(T(:,4),:))/4;
    Vii = Vii/sii+repmat(Cii,4,1)*(1-1/sii);
    VV(:,:,ii) = Vii;
  end
  h = animated_tetramesh(T,VV);
  axis equal;
  A = reshape(axis,2,[]);
  % grow by factor about center
  dA = [A(1,:)-A(2,:);A(2,:)-A(1,:)];
  axis(reshape(A+dA,1,prod(size(A))));
  % don't resize mesh
  axis manual;

end

author: alecjacobson

repo: gptoolbox


text:
function test_suite = test_survival_aft

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_SURVIVAL_AFT
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictions(testCase)
  verifyVarsEqual(testCase, getName(), {'pmu'}, {@(x)x(:,1)}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
  verifyVarsEqual(testCase, getName(), {'pmu'}, {@(x)x(:,2)}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
  verifyVarsEqual(testCase, getName(), {'pmu'}, {@(x)x(:,3)}, ...
    'RelTolElement', 0.05, 'RelTolRange', 0.01)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function y = cumprod( x, dim )

%   Disciplined geometric programming information for PROD:
%      CUMPROD(X,DIM) is a vectorized version of multiplication,
%      so in most cases it would be incompatible with DCPs. Therefore it
%      has not been implemented to support DCPs. DGPs however support
%      products more liberally. When PROD is used in a DCP, elements in
%      each subvector must satisfy the corresponding combination rule
%      for multiplication (see TIMES). For example, suppose that X looks
%      like this:
%         X = [ log-convex log-concave log-affine  ;
%               log-affine log-concave log-concave ]
%      Then CUMPROD(X,1) would be legal, but CUMPROD(X,2) would not,
%      because the top row contains the product of log-convex and
%      log-concave terms, in violation of the DGP ruleset.

error( nargchk( 1, 2, nargin ) );

%
% Basic argument check
%

sx = size( x );
if nargin < 2,
    dim = cvx_default_dimension( sx );
elseif ~cvx_check_dimension( dim ),
    error( 'Second argument must be a positive integer.' );
end

%
% Determine sizes
%

sx = [ sx, ones( 1, dim - length( sx ) ) ];
nx = sx( dim );

%
% Quick exit for empty arrays
%

if any( sx == 0 ),
    y = ones( sy );
    return
end

%
% Type check
%

persistent remap_1 remap_2 remap_3 remap_0
if isempty( remap_3 ),
    remap_0 = cvx_remap( 'zero' );
    remap_1 = cvx_remap( 'constant' );
    remap_2 = cvx_remap( 'log-convex' );
    remap_3 = cvx_remap( 'log-concave' );
end
vx = cvx_reshape( cvx_classify( x ), sx );
t0 = any( reshape( remap_0( vx ), sx ), dim );
t1 = all( reshape( remap_1( vx ), sx ), dim );
t2 = all( reshape( remap_2( vx ), sx ), dim ) | ...
     all( reshape( remap_3( vx ), sx ), dim );
t3 = t2 & t0;
ta = ( t1 | t3 ) + 2 * ( t2 & ~t3 );
nu = sort( ta(:) );
nu = nu([true;diff(nu)~=0]);
nk = length( nu );

%
% Quick exit for easy case
%

if nx == 1 && nu(1) > 0,
    y = x;
    return
end

%
% Permute and reshape, if needed
%

perm = [];
if nk > 1 || ( any( nu > 1 ) && nx > 1 ),
    if dim > 1 && any( sx( 1 : dim - 1 ) > 1 ),
        perm = [ dim, 1 : dim - 1, dim + 1 : length( sx ) ];
        x    = permute( x,  perm );
        ta   = permute( ta, perm );
        sx   = sx( perm );
        dim  = 1;
    end
    nv = prod( sx ) / nx;
    x  = reshape( x, nx, nv );
end

%
% Perform the computations
%

if nk > 1,
    y = cvx( [ nx, nv ], [] );
end
for k = 1 : nk,

    if nk == 1,
        xt = x;
    else
        tt = ta == nu( k );
        xt = cvx_subsref( x, ':', tt );
    end

    switch nu( k ),
        case 0,
            error( 'Disciplined convex programming error:\n   Invalid computation: cumprod( {%s} )', cvx_class( xt, true, true ) );
        case 1,
            yt = cvx( cumprod( cvx_constant( xt ), dim ) );
        case 2,
            yt = exp( cumsum( log( xt ), dim ) );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    if nk == 1,
        y = yt;
    else
        y = cvx_subsasgn( y, ':', tt, yt );
    end

end

%
% Reverse the reshaping and permutation steps
%

y = reshape( y, sx );
if ~isempty( perm ),
    y = ipermute( y, perm );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function L = hash2landmark(H)
% L = hash2landmark(H)
%  Convert a set of <time hash> pairs ready from store
%  into a set of 4-entry landmarks <t1 f1 f2 dt>.
%  If H is 3 cols, first col (song ID) is discarded.
% 2008-12-29 Dan Ellis dpwe@ee.columbia.edu

% Hash value is 20 bits: 8 bits of F1, 6 bits of F2-F1, 6 bits of delta-T

%H = uint32(L(:,1));
%F1 = rem(round(L(:,2)),2^7);
%F2 = rem(round(L(:,4)),2^7);
%DT = rem(abs(round(L(:,3) - L(:,1))), 2^6);
%H = [H,uint32(F1*(2^13)+F2*(2^6)+DT)];

if size(H,2) == 3
  H = H(:,[2 3]);
end

H1 = H(:,1);
H2 = double(H(:,2));
F1 = floor(H2/(2^12));
H2 = H2 - (2^12)*F1;
F1 = F1 + 1;
DF = floor(H2/(2^6));
H2 = H2 - (2^6)*DF;
if DF > 2^5
  DF = DF-2^6;
end
F2 = F1+DF;

DT = H2;

L = [H1,F1,F2,DT];


text:
function [vin, vout] = EdgeExtract(img, imgf)
    S = im2double(img);
    vin = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
    img_filtered = im2double(imgf);
    S = im2double(img_filtered);
    vout = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
end

author: jimmy-ren