text:
function pass = test_nonlinSys1Breaks_C1(pref)
% Test 2x2 system (sin/cos). This is pecewiseificaion of the test
%       test_nonlinearSystem1
%
% Asgeir Birkisson, April 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

d = [-pi 0 pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Piecewise (chebcolloc1):
pref.discretization = @chebcolloc1;

A = chebop(@(x,u,v) [u - diff(v,2) + u.^2; diff(u) + sin(v)],d);
A.lbc = @(u,v) u-1;
A.rbc = @(u,v) [v-1/2; diff(v)];

u = mldivide(A, f, pref);
u1 = u{1}; u2 = u{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = chebfun(A.rbc(u1,u2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(1) = norm( chebfun(A(x, u1, u2))) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
pass(3) = norm(u1jump) < tol && norm(u2jump) < tol;

end


text:
function xList=RungeKSteps(xInit,t0,f,deltaT,numSteps,order,solutionChoice,onlyEnd)
%%RUNGEKSTEPS Perform multiple steps of Runge-Kutta propagation with a
%             fixed time interval between steps.
%
%INPUTS: xInit The xDimX1 initial value of the state (scalar or vector)
%              over which integration is being performed.
%           t0 The scalar time at which xInit is taken.
%            f f(xVal,curT) returns the derivative of xVal taken at time
%              curT.
%       deltaT The scalar size of the steps in the Runge-Kutta integration.
%     numSteps The scalar number of steps over which Runge-Kutta
%              integration is performed.
%        order The order of the Runge-Kutta method. If this parameter is
%              omitted, then the default order of 4 is used. Order can
%              range from 1 to 7.
% solutionChoice When multiple formulae are implemented, this selects which
%              one to use. Otherwise, this parameter is not used.
%              Currently, only the order=4 method has multiple solutions
%              implemented in which case omitting this parameter or setting
%              it to zero used the Dormand and Prince Algorithm, and
%              setting it to 1 uses the Fehlberg algorithm.
%      onlyEnd If this is true, only the value of x at the final step is
%              returned, not any of the intermediate values. The default if
%              this is omitted or an empty matrix is passed is false.
%
%OUTPUTS: xList If onlyEnd is false, then this is the xDimX(numSteps+1 The
%               first column of xList is xInit. The subsequent columns
%               are the state propagated forward at intervals of deltaT.
%               Otherwise, if onlyEnd is true, then only the final
%               propagated value is returned and not any of the
%               intermediate steps.
%
%This function calls the RungeKStep function to propagate forward the
%target state each step of duration deltaT. See the comments in RungeKStep
%for more information.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<8||isempty(onlyEnd))
        onlyEnd=false;
    end

    if(nargin<7||isempty(solutionChoice))
        solutionChoice=0;
    end

    if(nargin<6||isempty(order))
        order=4;
    end

    if(onlyEnd)
        xList=xInit;
        curT=t0;
        for curStep=1:numSteps
            xList=RungeKStep(xList,curT,f,deltaT,[],order,solutionChoice);
            curT=curT+deltaT;
        end
    else%If all of the steps are desired.
        xDim=size(xInit,1);
        xList=zeros(xDim,numSteps+1);
        xList(:,1)=xInit;

        curT=t0;
        for curStep=1:numSteps
            xList(:,curStep+1)=RungeKStep(xList(:,curStep),curT,f,deltaT,[],order,solutionChoice);
            curT=curT+deltaT;
        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 [error_x, error_y, fex, fey, ae] = stokespost_q1p0_p(jmpx,jmpy,els,xy,ev)
%stokespost_q1p0_p  computes Poisson error estimator for Q1-P0
%   [err_x, err_y, fex, fey, ae] = stokespost_q1p0_p(jmpx,jmpy,els,xy,ev);
%   input
%          jmpx, jmpy     component elementwise edge stress jumps
%          els            elementwise edge lengths
%          xy             vertex coordinate vector
%          ev             element mapping matrix
%   output
%          err_x, err_y   component of velocity elementwise error estimate
%          fex, fey       component elementwise rhs vectors
%          ae             elementwise Poisson problem matrices
%
%   IFISS function: DJS; 8 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      fprintf('computing local error estimator... ')
      x=xy(:,1); y=xy(:,2);
      nel=length(ev(:,1));
      error_x=zeros(nel,1);  error_y=zeros(nel,1);
%
% set up 3x3 Gauss points
      gpt=sqrt(0.6);
      s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
      s(2) =  gpt; t(2) = -gpt; wt(2)=25/81;
      s(3) =  gpt; t(3) =  gpt; wt(3)=25/81;
      s(4) = -gpt; t(4) =  gpt; wt(4)=25/81;
      s(5) =  0.0; t(5) = -gpt; wt(5)=40/81;
      s(6) =  gpt; t(6) =  0.0; wt(6)=40/81;
      s(7) =  0.0; t(7) =  gpt; wt(7)=40/81;
      s(8) = -gpt; t(8) =  0.0; wt(8)=40/81;
      s(9) =  0.0; t(9) =  0.0; wt(9)=64/81;
%
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(ev(:,ivtx));
        yl_v(:,ivtx) = y(ev(:,ivtx));
		end
        ae = zeros(nel,5,5); elerrx=zeros(5,nel);  elerry=zeros(5,nel);
        fex = zeros(nel,5); fey = zeros(nel,5);
% loop over Gauss points
         for igpt = 1:9
         sigpt=s(igpt);
         tigpt=t(igpt);
         wght=wt(igpt);
% evaluate derivatives etc
         [jac_v,invjac_v,phi_v,dphidx_v,dphidy_v] = deriv(sigpt,tigpt,xl_v,yl_v);
         [psi_v,dpsidx_v,dpsidy_v] = qderiv(sigpt,tigpt,xl_v,yl_v);
            for j = 1:5
               for i = 1:5
               ae(:,i,j) = ae(:,i,j)+wght*dpsidx_v(:,i+4).*dpsidx_v(:,j+4).*invjac_v(:);
               ae(:,i,j) = ae(:,i,j)+wght*dpsidy_v(:,i+4).*dpsidy_v(:,j+4).*invjac_v(:);
               end
            end
% end of Gauss point loop
         end
%
% include edge jumps (evaluated at the midpoint)
         for ee = 1:4
         fex(:,ee) = fex(:,ee) - jmpx(:,ee) .* els(:,ee)*(1/3);
         fey(:,ee) = fey(:,ee) - jmpy(:,ee) .* els(:,ee)*(1/3);
          end
%
% solve for local estimate
%         err=ae\fe;
%         elerr_p(ielem,1) = err'*fe;
%% sequential code
         for ielem = 1:nel
		 elerrx(:,ielem) = squeeze(ae(ielem,1:5,1:5))\(fex(ielem,1:5)');
		 elerry(:,ielem) = squeeze(ae(ielem,1:5,1:5))\(fey(ielem,1:5)');
	     end
%%
         for ivtx=1:5,
	     error_x(:) = error_x(:) + fex(:,ivtx) .* elerrx(ivtx,:)';
		 error_y(:) = error_y(:) + fey(:,ivtx) .* elerry(ivtx,:)';
		 end
%%
         fprintf('done.\n')
		 return

author: johannesgerer

repo: jburkardt-m


text:
% New SVD based initialization strategy for Non-negative Matrix Factorization
% Hanli Qiao

function [W, H] = Qiao_SVD_Init(Z, rank)
    [u, s, v, p] = ChoosingR(Z);
    W = abs(u(:,1:rank));
    H = abs(s(1:rank,:)*v');
end

function [u, s, v, p] = ChoosingR(Z)
    [u,s,v] = svd(Z);
    sum1= sum(s);
    sum2=sum(sum1);
    extract=0;
    p = 0;
    dsum=0;
    while(extract/sum2<0.90)
        p = p + 1;
        dsum=dsum+s(p,p);
        extract=dsum;
    end
end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% [mu,S] = weighted_mean_cov(x,w)
% x is n x d, w is n x 1
function [mu,S] = weighted_mean_cov(x,w)

%[n,d] = size(x);

z = sum(w);
mu = sum(bsxfun(@times,x,w(:)),1)/z;

diffs = bsxfun(@minus,x,mu);
diffs = bsxfun(@times,diffs,sqrt(w));
S = (diffs'*diffs)/z;

author: kristinbranson

repo: JAABA


text:
function [A, Z, title, key, mtype] = RBfix (filename)
%RBFIX read a possibly corrupted matrix from a R/B file
% (assembled format only).  Usage:
%
% [A Z title key mtype] = RBfix (filename)
%
% The Rutherford/Boeing format stores a sparse matrix in a file in compressed-
% column form, using 3 arrays: Ap, Ai, and Ax.  The row indices of entries in
% A(:,j) are in Ai(p1:p2) and the corresponding numerical values are Ax(p1:p2),
% where p1 = Ap(j) and p2 = Ap(j+1)-1.  The row indices ought to be sorted, and
% no duplicates should appear, but this function ignores that requirement.
% Duplicate entries are summed if they exist, and A is returned with sorted
% columns.  Symmetric matrices are stored with just their lower triangular
% parts in the file.  Normally, it is an error if entries are present in the
% upper triangular part of a matrix that is declared in the file to be
% symmetric.  This function simply ignores those entries.
%
% If CHOLMOD is installed, this function is faster and uses less memory.
%
% Example:
%
%   load west0479
%   RBwrite ('mywest', west0479, [ ], 'My west0479 file', 'west0479') ;
%   [A Z title key mtype] = RBfix ('mywest') ;
%   isequal (A, west0479)
%   title, key, mtype
%
% See also mread, RBread, RBwrite, RBreade, sparse2.

% Optionally uses the CHOLMOD sparse2 mexFunction.

% Copyright 2007, Timothy A. Davis

%-------------------------------------------------------------------------------
% read in the raw contents of the Rutherford/Boeing file
%-------------------------------------------------------------------------------

[mtype Ap Ai Ax title key nrow] = RBraw (filename) ;
mtype = lower (mtype) ;

%-------------------------------------------------------------------------------
% determine dimension, number of entries, and convert numerical entries
%-------------------------------------------------------------------------------

% number of columns
ncol = length (Ap) - 1 ;

% number of entries
nz = length (Ai) ;

% check column pointers
if (any (Ap ~= sort (Ap)) | (Ap (1) ~= 1) | (Ap (ncol+1) - 1 ~= nz))	    %#ok
    error ('invalid column pointers') ;
end

% check row indices
if ((double (max (Ai)) > nrow) | double (min (Ai)) < 1)			    %#ok
    error ('invalid row indices') ;
end

% Ax can be empty, for a p*a matrix
if (~isempty (Ax))
    if (mtype (1) == 'c')
	% Ax is real, with real/imaginary parts interleaved
	if (2 * nz ~= length (Ax))
	    error ('invalid matrix') ;
	end
	Ax = Ax (1:2:end) + (1i * Ax (2:2:end)) ;
    elseif (mtype (1) == 'i')
	Ax = double (Ax) ;
    end
    % numerical values must be of the right size
    if (nz ~= length (Ax))
	error ('invalid matrix') ;
    end
end

%-------------------------------------------------------------------------------
% create the triplet form
%-------------------------------------------------------------------------------

% construct column indices
Aj = zeros (nz,1) ;
for j = 1:ncol
    p1 = Ap (j) ;
    p2 = Ap (j+1) - 1 ;
    Aj (p1:p2) = j ;
end

%-------------------------------------------------------------------------------
% create the sparse matrix form
%-------------------------------------------------------------------------------

if (exist ('sparse2') == 3)						    %#ok
    % Use sparse2 in CHOLMOD.  It's faster, allows integer Ai and Aj, and
    % returns the Z matrix as the 2nd output argument.
    if (isempty (Ax))
	Ax = 1 ;
    end
    % numerical matrix
    [A Z] = sparse2 (Ai, Aj, Ax, nrow, ncol) ;
else
    % stick with MATLAB, without CHOLMOD.  This is slower and takes more memory.
    Ai = double (Ai) ;
    Aj = double (Aj) ;
    if (isempty (Ax))
	% pattern-only matrix
	A = spones (sparse (Ai, Aj, 1, nrow, ncol)) ;
	Z = sparse (nrow, ncol) ;
    else
	% numerical matrix
	A = sparse (Ai, Aj, Ax, nrow, ncol) ;
	% determine the pattern of explicit zero entries
	S = spones (sparse (Ai, Aj, 1, nrow, ncol)) ;
	Z = S - spones (A) ;
    end
end

% check for entries in upper part
if (any (mtype (2) == 'shz') & nnz (triu (A,1) > 0))			    %#ok
    fprintf ('entries in upper triangular part of %s matrix ignored\n', mtype);
end

% add the upper triangular part
if (mtype (2) == 's')
    A = A + tril (A,-1).' ;
    Z = Z + tril (Z,-1)' ;
elseif (mtype (2) == 'h')
    A = A + tril (A,-1)' ;
    Z = Z + tril (Z,-1)' ;
elseif (mtype (2) == 'z')
    A = A - tril (A,-1).' ;
    Z = Z + tril (Z,-1)' ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Hv] = LogisticHv(v,w,X,y)
% v(feature,1) - vector that we will multiply Hessian by
% w(feature,1)
% X(instance,feature)
% y(instance,1)

sig = 1./(1+exp(-y.*(X*w)));
Hv = X.'*(sig.*(1-sig).*(X*v));

author: zouchuhang

repo: LayoutNet


text:
function bcOut = bcReform(dom, bcInput, isIorF)
%BCREFORM    Convert a general BC from chebgui to IVP/FVP form
%
% BCOUT = BCREFORM(DOM, BCINPUT, ISIORF), where
%
%   DOM:     The domain that the problem is specified.
%   BCINPUT: The input from the BC field of CHEBGUI.
%   ISIORF:  Is equal to 1 if we're solving an initial value problem, and 2
%            if we're solving a final value problem.
%
% returns
%   BCOUT:   Boundary conditions that can be assigned to the LBC or RBC field of
%            a chebop.
%
% This method is used for converting input that usually would be assigned to
% the BC field of a chebop to form that can be assigned to the LBC or RBC field.
% It is required since in CHEBGUI, we specify conditions for ODEs on the format
%   u(0) = 1,
% which can then be imposed via
%   N.bc = @(x,u) u(0) - 1;
% On the other hand, for IVPs/FVPs, the format
%   u = 1
% is more useful, since that can then be imposed via
%   N.lbc = @(u) u;

% Find where whitespace and commas appear in the domain string, so that we can
% determine what the endpoints of the domain are:
spaceLoc = union(strfind(dom, ' '), strfind(dom, ','));

% Look at the BCINPUT string, and throw away parts of the string that we need to
% get rid of so that we convert it to an anonymous function that is suitable for
% N.lbc/rbc. Essentially, this is converting a string such as
%   u(0) = 1
% to
%   u = 1.
% We need to throw different parts of the BCINPUT string away, depending on
% whether we're solving an IVP or FVP:
if ( isIorF ==  1 )
    % Initial value problem.
    firstSpace = min(spaceLoc);
    leftDomStr = dom(2:firstSpace-1);
    % String that we're going to be throwing away:
    throwString = ['(' ,leftDomStr, ')'];
else
    % Final value problem.
    lastSpace = max(spaceLoc);
    rightDomStr = dom(lastSpace+1:end-1);
    throwString = ['(' ,rightDomStr, ')'];
end

% Clear out the parts of the string required!
bcOut = strrep(bcInput, throwString, '');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ x, seed ] = ball_unit_sample_nd ( n, seed )

%*****************************************************************************80
%
%% BALL_UNIT_SAMPLE_ND picks a random point in the unit ball in ND.
%
%  Discussion:
%
%    N-1 random Givens rotations are applied to the point ( 1, 0, 0, ..., 0 ).
%
%    The I-th Givens rotation is in the plane of coordinate axes I and I+1,
%    and has the form:
%
%     [ cos ( theta )  - sin ( theta ) ] * x(i)      = x'(i)
%     [ sin ( theta )    cos ( theta ) ]   x(i+1)      x'(i+1)
%
%    Finally, a scaling is applied to set the point at a distance R
%    from the center, in a way that results in a uniform distribution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the space.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), the random point.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  x(1) = 1.0;
  x(2:n) = 0.0;

  for i = 1 : n-1

    [ r, seed ] = r8_uniform_01 ( seed );
    random_cosine = 2.0 * r - 1.0;
    [ r, seed ] = r8_uniform_01 ( seed );
    random_sign = 2 * floor ( 2.0 * r ) - 1;
    [ r, seed ] = r8_uniform_01 ( seed );
    random_sine = random_sign * sqrt ( 1.0 - random_cosine * random_cosine );

    xi = x(i);
    x(i  ) = random_cosine * xi;
    x(i+1) = random_sine   * xi;

  end

  [ r, seed ] = r8_uniform_01 ( seed );

  r = r^( 1.0 / n );

  x(1:n) = r * x(1:n);

  return
end


text:
%Install_lightspeed
% Compiles mex files for the lightspeed library.

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

% thanks to Kevin Murphy for suggesting this routine.
% thanks to Ruben Martinez-Cantin for UNDERSCORE_LAPACK_CALL

% check that we are in the right directory
d = dir('install_lightspeed.m');
if length(d) == 0
	error('You must run install_lightspeed from within the lightspeed directory.  Use cd to change the current directory.');
end
fprintf('Compiling lightspeed 2.7 mex files...\n');

% Matlab version
v = sscanf(version,'%d.%d.%*s (R%d) %*s');
% v(3) is the R number
% could also use v(3)>=13
atleast65 = (v(1)>6 || (v(1)==6 && v(2)>=5));
atleast73 = (v(1)>7 || (v(1)==7 && v(2)>=3));
atleast75 = (v(1)>7 || (v(1)==7 && v(2)>=5));
atleast76 = (v(1)>7 || (v(1)==7 && v(2)>=6));
atleast78 = (v(1)>7 || (v(1)==7 && v(2)>=8)); % R2009a
atleast82 = (v(1)>8 || (v(1)==8 && v(2)>=2)); % R2013b
atleast83 = (v(1)>8 || (v(1)==8 && v(2)>=3)); % R2014a
if atleast73
	% largeArrayDims and mwSize were added in version 7.3 (R2006b)
	% http://www.mathworks.com/help/techdoc/rn/bqt6wtq.html
	flags = ' -largeArrayDims ';
else
	flags = ' -DmwSize=int -DmwIndex=int ';
end
if atleast83
	flags = [flags ' -silent '];
end
% To catch ANSI C violations with gcc:
%flags = [flags ' CFLAGS=''-Wall -std=c99'' '];
replace_repmat = ~atleast82;

% copy matlab's original repmat.m as xrepmat.m
if replace_repmat && exist('xrepmat.m') ~= 2
  w = fullfile(matlabroot,'toolbox','matlab','elmat','repmat.m');
  cmd = ['"' w '" xrepmat.m'];
  if ispc
    system(['copy ' cmd]);
  else
    system(['cp -rp ' cmd]);
  end
end

% these are done first to initialize mex
eval(['mex' flags '-c flops.c']);
eval(['mex' flags 'sameobject.c']);
eval(['mex' flags 'int_hist.c']);
eval(['mex' flags '-c mexutil.c']);
eval(['mex' flags '-c util.c']);

libdir = '';
if ispc || atleast83
	[compiler,options] = mexcompiler;
	libdir = options.LIBLOC;
	engmatopts = [compiler 'engmatopts.bat'];
elseif ismac
	options = struct;
	% this installer is set up for 64-bit MacOSX 10.6 with gcc-4.0
	% if you are using something else, run 'mex -v -c flops.c'
	% and use the output to change these strings
	options.COMPILER = 'gcc-4.0';
	options.COMPFLAGS = '-fno-common -no-cpp-precomp -arch x86_64 -isysroot /Developer/SDKs/MacOSX10.6.sdk -mmacosx-version-min=10.6  -fexceptions';
	options.OPTIMFLAGS = '-O -DNDEBUG';
else
	options = struct;
	options.COMPILER = 'cc';
	options.COMPFLAGS = '-fPIC';
	options.OPTIMFLAGS = '-O';
end

% Routines that use LAPACK
lapacklib = '';
blaslib = '';
lapackflags = flags;
if atleast78
	lapackflags = [lapackflags ' -DBLAS64'];
end
if ispc
  if strncmp(compiler,'MSVC',4)
    if atleast65
      % version >= 6.5
      lapacklib = fullfile(libdir,'libmwlapack.lib');
    end
  else
    lapacklib = fullfile(libdir,'libmwlapack.lib');
  end
  if atleast75
    blaslib = fullfile(libdir,'libmwblas.lib');
  end
  %%% Paste the location of libmwlapack.lib %%%
  %lapacklib = '';
  if ~exist(lapacklib,'file')
    lapacklib = 'dtrsm.c';
    fprintf('libmwlapack.lib was not found.  To get additional optimizations, paste its location into install_lightspeed.m\n');
  else
    fprintf('Using the lapack library at %s\n',lapacklib);
  end
else
  % in version 7.5, non-PC systems do not need to specify lapacklib,
  % but they must use an underscore when calling lapack routines
	% http://www.mathworks.com/help/techdoc/matlab_external/br_2m24-1.html
	lapackflags = [lapackflags ' -DUNDERSCORE_LAPACK_CALL'];
  if atleast76
    lapacklib = '-lmwlapack';
    blaslib = '-lmwblas';
  end
end
disp(['mex' lapackflags ' solve_triu.c "' lapacklib '" "' blaslib '"'])
eval(['mex' lapackflags ' solve_triu.c "' lapacklib '" "' blaslib '"']);
eval(['mex' lapackflags ' solve_tril.c "' lapacklib '" "' blaslib '"']);

if ispc
  % Windows
  %if exist('util.obj','file')
  eval(['mex' flags 'addflops.c flops.obj'])
	if atleast78
		[success,message] = mkdir('@double');
		eval(['mex' flags 'gammaln.c util.obj -outdir @double'])
	else
		eval(['mex' flags 'gammaln.c util.obj'])
	end
  eval(['mex' flags 'digamma.c util.obj'])
  eval(['mex' flags 'trigamma.c util.obj'])
  eval(['mex' flags 'tetragamma.c util.obj'])
	eval(['mex' flags 'setnonzeros.c'])
  if strncmp(compiler,'MSVC',4)
		clear random.dll randomseed randbinom randgamma sample_hist
    disp(['install_random.bat "' options.VSINSTALLDIR '" ' options.vcvarsopts]);
    system(['install_random.bat "' options.VSINSTALLDIR '" ' options.vcvarsopts]);
    eval(['mex' flags 'randomseed.c util.obj random.lib'])
    eval(['mex' flags 'randbinom.c mexutil.obj util.obj random.lib'])
    eval(['mex' flags 'randgamma.c mexutil.obj util.obj random.lib'])
    eval(['mex' flags 'sample_hist.c util.obj random.lib'])
  else
    fprintf('mexcompiler is not MSVC. The randomseed() function will have no effect.');
    eval(['mex' flags 'randomseed.c util.obj random.c'])
    eval(['mex' flags 'randbinom.c mexutil.obj util.obj random.c'])
    eval(['mex' flags 'randgamma.c mexutil.obj util.obj random.c'])
    eval(['mex' flags 'sample_hist.c util.obj random.c'])
  end
	if replace_repmat
		eval(['mex' flags 'repmat.c mexutil.obj'])
	else
		eval(['mex' flags 'repmat.c mexutil.obj -output xrepmat'])
	end
  try
    % standalone programs
    % compilation instructions are described at:
    % http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_external/ch1_im15.html#27765
    if atleast78
			%disp('lightspeed''s matfile utility is not supported for this version of Matlab, but everything else should work.');
		elseif atleast65
      % -V5 is required for Matlab >=6.5
      mex('-f',engmatopts,'matfile.c','-V5');
    else
      mex('-f',engmatopts,'matfile.c');
    end
    % uncomment the line below if you want to build test_flops.exe
    % This program lets you check the flop counts on your processor.
    % mex('-f',engmatopts,'tests/test_flops.c');
  catch
    disp('Could not install the standalone programs.');
    disp(lasterr)
  end
else
  % UNIX
  eval(['mex' flags 'addflops.c flops.o'])
	if atleast78
		[success,message] = mkdir('@double');
		eval(['mex' flags 'gammaln.c util.o -lm -outdir @double'])
	else
		eval(['mex' flags 'gammaln.c util.o -lm'])
	end
  eval(['mex' flags 'digamma.c util.o -lm'])
  eval(['mex' flags 'trigamma.c util.o -lm'])
  eval(['mex' flags 'tetragamma.c util.o -lm'])
	eval(['mex' flags 'setnonzeros.c'])
  if ismac
    % thanks to Nicholas Butko for these mac-specific lines
		clear librandom.dylib randomseed randbinom randgamma sample_hist
		cmd = [options.COMPILER ' ' options.COMPFLAGS ' ' options.OPTIMFLAGS ' -c random.c; ' options.COMPILER ' ' options.COMPFLAGS ' -dynamiclib -Wl,-install_name,"`pwd`/librandom.dylib" -o librandom.dylib random.o'];
		disp(cmd);
		system(cmd)
    eval(['mex' flags 'randomseed.c util.o librandom.dylib -lm'])
    eval(['mex' flags 'randbinom.c mexutil.o util.o librandom.dylib -lm'])
    eval(['mex' flags 'randgamma.c mexutil.o util.o librandom.dylib -lm'])
    eval(['mex' flags 'sample_hist.c util.o librandom.dylib -lm'])
  else
    % this command only works on linux
		clear librandom.so randomseed randbinom randgamma sample_hist
    cmd = [options.COMPILER ' ' options.COMPFLAGS ' ' options.OPTIMFLAGS ' -c random.c; ' options.COMPILER ' ' options.COMPFLAGS ' -shared -Wl,-E -Wl,-soname,`pwd`/librandom.so -o librandom.so random.o'];
		disp(cmd);
		system(cmd)
    eval(['mex' flags 'randomseed.c util.o librandom.so -lm'])
    eval(['mex' flags 'randbinom.c mexutil.o util.o librandom.so -lm'])
    eval(['mex' flags 'randgamma.c mexutil.o util.o librandom.so -lm'])
    eval(['mex' flags 'sample_hist.c util.o librandom.so -lm'])
  end
	if replace_repmat
		eval(['mex' flags 'repmat.c mexutil.o'])
	else
		eval(['mex' flags 'repmat.c mexutil.o -output xrepmat'])
	end
  try
    % standalone programs
    if atleast78
			disp('lightspeed''s matfile utility is not supported for this version of Matlab');
		elseif atleast65
      % -V5 is required only for Matlab >=6.5
      mex -f matopts.sh matfile.c -V5
    else
      mex -f matopts.sh matfile.c
    end
    % uncomment the line below if you want to build test_flops.exe
    % This program lets you check the flop counts on your processor.
    % mex -f matopts.sh tests/test_flops.c
  catch
    disp('Could not install the standalone programs.');
    disp(lasterr);
    fprintf('Note: if matlab cannot find matopts.sh, your installation of matlab is faulty.\nIf you get this error, don''t worry, lightspeed should still work.');
  end
end

addpath(genpath(pwd))
fprintf('Done.\n');
fprintf('Type "test_lightspeed" to verify the installation.\n');

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
% -----------------------------------------------------------------  %
% Matlab Programs included the Appendix B in the book:               %
%  Xin-She Yang, Engineering Optimization: An Introduction           %
%                with Metaheuristic Applications                     %
%  Published by John Wiley & Sons, USA, July 2010                    %
%  ISBN: 978-0-470-58246-6,   Hardcover, 347 pages                   %
% -----------------------------------------------------------------  %
% Citation detail:                                                   %
% X.-S. Yang, Engineering Optimization: An Introduction with         %
% Metaheuristic Application, Wiley, USA, (2010).                     %
%                                                                    %
% http://www.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html %
% http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0470582464.html  %
% -----------------------------------------------------------------  %
% ===== ftp://  ===== ftp://   ===== ftp:// =======================  %
% Matlab files ftp site at Wiley                                     %
% ftp://ftp.wiley.com/public/sci_tech_med/engineering_optimization   %
% ----------------------------------------------------------------   %

% Simulated Annealing (by X-S Yang, Cambridge University)            %
% Usage: B2_sa                                                       %
% For the constrained optimization, please see the file: sa_mincon.m %
% ------------------------------------------------------------------ %

disp('Simulating ... it will take a minute or so!');
% Rosenbrock's function with f*=0 at (1,1)
fstr='(1-x)^2+100*(y-x^2)^2';
% Convert into an inline function
f=vectorize(inline(fstr));
% Show the topography of the objective function
range=[-2 2 -2 2];
xgrid=range(1):0.1:range(2); ygrid=range(3):0.1:range(4);
[x,y]=meshgrid(xgrid,ygrid);
surfc(x,y,f(x,y));
% Initializing parameters and settings
T_init = 1.0;       % Initial temperature
T_min =  1e-10;     % Finial stopping temperature
F_min = -1e+100;    % Min value of the function
max_rej=5000;       % Maximum number of rejections
max_run=500;        % Maximum number of runs
max_accept = 250;   % Maximum number of accept
k = 1;              % Boltzmann constant
alpha=0.95;         % Cooling factor
Enorm=1e-8;         % Energy norm (eg, Enorm=1e-8)
guess=[2 2];        % Initial guess
% Initializing the counters i,j etc
i= 0; j = 0; accept = 0; totaleval = 0;
% Initializing various values
T = T_init;
E_init = f(guess(1),guess(2));
E_old = E_init; E_new=E_old;
best=guess;  % initially guessed values
% Starting the simulated annealling
while ((T > T_min) & (j <= max_rej) & E_new>F_min)
    i = i+1;
    % Check if max numbers of run/accept are met
    if (i >= max_run) | (accept >= max_accept)
    % Cooling according to a cooling schedule
        T = alpha*T;
        totaleval = totaleval + i;
        % reset the counters
        i = 1;  accept = 1;
    end
    % Function evaluations at new locations
      ns=guess+rand(1,2)*randn;
      E_new = f(ns(1),ns(2));
    % Decide to accept the new solution
    DeltaE=E_new-E_old;
    % Accept if improved
    if (-DeltaE > Enorm)
        best = ns; E_old = E_new;
        accept=accept+1;   j = 0;
    end
    % Accept with a small probability if not improved
    if (DeltaE<=Enorm & exp(-DeltaE/(k*T))>rand );
        best = ns; E_old = E_new;
        accept=accept+1;
    else
        j=j+1;
    end
    % Update the estimated optimal solution
    f_opt=E_old;
end
% Display the final results
disp(strcat('Obj function  :',fstr));
disp(strcat('Evaluations   :', num2str(totaleval)));
disp(strcat('Best  solution:', num2str(best)));
disp(strcat('Best objective:', num2str(f_opt)));

author: Sable

repo: mcbench-benchmarks


text:
function output = call_xpressfico_miqp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
model = yalmip2xpress(interfacedata);

if options.savedebug
    save xpressdebug model
end

solvertime = tic;
if isempty(model.extra.integer_variables) & isempty(model.extra.binary_variables) & isempty(model.extra.semicont_variables) & isempty(model.sos)
    if options.verbose
        [x,fval,exitflag,output,lambda] = xprsqp(model.H,model.f,model.A,model.b,model.rtype,model.lb,model.ub,model.ops);
    else
        evalc('[x,fval,exitflag,output,lambda] = xprsqp(model.H,model.f,model.A,model.b,model.rtype,model.lb,model.ub,model.ops);');
    end
else
    if options.verbose
        [x,fval,exitflag,output] = xprsmiqp(model.H,model.f,model.A,model.b,model.rtype,model.ctype,model.clim,model.sos,model.lb,model.ub,[],model.ops);
    else
        evalc('[x,fval,exitflag,output] = xprsmiqp(model.H,model.f,model.A,model.b,model.rtype,model.ctype,model.clim,model.sos,model.lb,model.ub,[],model.ops);        ');
    end
    lambda = [];
end
solvertime = toc(solvertime);

if ~isempty(lambda)
    D_struc = [lambda.lin];
else
    D_struc = [];
end

if length(x) == length(model.f)
    if ~isempty(model.extra.NegatedSemiVar)
        x(model.extra.NegatedSemiVar) = -x(model.extra.NegatedSemiVar);
    end
end

% Check, currently not exhaustive...
switch exitflag
    case {1}
        problem = 0;
    case {-2}
        problem = 1; % Infeasible
    case {0,-4,-5}
        problem = 3;
    case -8
        problem = -11;
    otherwise
        problem = -1;
end

% Save all data sent to solver?
if options.savesolverinput
    solverinput.f = f;
    solverinput.A = A;
    solverinput.b = b;
    solverinput.ctype = ctype;
    solverinput.rtype = rtype;
    solverinput.lb = lb;
    solverinput.ub = ub;
    solverinput.x0 = [];
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fval = fval;
    solveroutput.exitflag = exitflag;
    solveroutput.output = output;
    solveroutput.lambda = lambda;
else
    solveroutput = [];
end

if isempty(x)
    x = zeros(length(model.f),1);
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function sVF1 = plus(sVF1, sVF2)
%
% Syntax
%   sVF = sVF1 + sVF2
%   sVF = sVF + v
%
% Input
%  sVF1 - @S2VectorFieldTri
%  sVF2 - @S2VectorField
%  v    - @vector3d
%
% Output
%  sF - @S2FunTri
%

% first should be S2VectorFieldTri
if ~isa(sVF1,'S2VectorFieldTri'), [sVF1,sVF2] = deal(sVF2,sVF1); end

if isa(sVF2,'vector3d')
  sVF1.values = sVF1.values +  a;
elseif isa(sVF2,'S2VectorFieldTri')
  sVF1.values = sVF1.values + sVF2.values;
else
  sVF1.values = sVF1.values + sVF2.eval(SVF.vertices);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function sig = MakeSignal(Name,n)
% MakeSignal -- Make artificial signal
%  Usage
%    sig = MakeSignal(Name,n)
%  Inputs
%    Name   string: 'HeaviSine', 'Bumps', 'Blocks',
%            'Doppler', 'Ramp', 'Cusp', 'Sing', 'HiSine',
%            'LoSine', 'LinChirp', 'TwoChirp', 'QuadChirp',
%            'MishMash', 'WernerSorrows' (Heisenberg),
%            'Leopold' (Kronecker), 'Piece-Regular' (Piece-Wise Smooth),
%	     'Riemann','HypChirps','LinChirps', 'Chirps', 'Gabor'
%	     'sineoneoverx','Cusp2','SmoothCusp','Gaussian'
%	     'Piece-Polynomial' (Piece-Wise 3rd degree polynomial)
%    n      desired signal length
%  Outputs
%    sig    1-d signal
%
%  References
%    Various articles of D.L. Donoho and I.M. Johnstone
%
   if nargin > 1,
	t = (1:n) ./n;
   end
	if strcmp(Name,'HeaviSine'),
	    sig = 4.*sin(4*pi.*t);
	    sig = sig - sign(t - .3) - sign(.72 - t);
	elseif strcmp(Name,'Bumps'),
	    pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
	    hgt = [ 4  5   3   4  5  4.2 2.1 4.3  3.1 5.1 4.2];
	    wth = [.005 .005 .006 .01 .01 .03 .01 .01  .005 .008 .005];
	    sig = zeros(size(t));
	    for j =1:length(pos)
	       sig = sig + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
	    end
	elseif strcmp(Name,'Blocks'),
	    pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
	    hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3  (-3.1) 2.1 (-4.2)];
	    sig = zeros(size(t));
	    for j=1:length(pos)
	        sig = sig + (1 + sign(t-pos(j))).*(hgt(j)/2) ;
	    end
	elseif strcmp(Name,'Doppler'),
	    sig = sqrt(t.*(1-t)).*sin((2*pi*1.05) ./(t+.05));
	elseif strcmp(Name,'Ramp'),
	    sig = t - (t >= .37);
	elseif strcmp(Name,'Cusp'),
	    sig = sqrt(abs(t - .37));
	elseif strcmp(Name,'Sing'),
	    k = floor(n * .37);
	    sig = 1 ./abs(t - (k+.5)/n);
	elseif strcmp(Name,'HiSine'),
	    sig = sin( pi * (n * .6902) .* t);
	elseif strcmp(Name,'LoSine'),
	    sig = sin( pi * (n * .3333) .* t);
	elseif strcmp(Name,'LinChirp'),
	    sig = sin(pi .* t .* ((n .* .500) .* t));
	elseif strcmp(Name,'TwoChirp'),
	    sig = sin(pi .* t .* (n .* t)) + sin((pi/3) .* t .* (n .* t));
	elseif strcmp(Name,'QuadChirp'),
	    sig = sin( (pi/3) .* t .* (n .* t.^2));
	elseif strcmp(Name,'MishMash'),  % QuadChirp + LinChirp + HiSine
	    sig = sin( (pi/3) .* t .* (n .* t.^2)) ;
	    sig = sig +  sin( pi * (n * .6902) .* t);
	    sig = sig +  sin(pi .* t .* (n .* .125 .* t));
	elseif strcmp(Name,'WernerSorrows'),
	    sig = sin( pi .* t .* (n/2 .* t.^2)) ;
	    sig = sig +  sin( pi * (n * .6902) .* t);
	    sig = sig +  sin(pi .* t .* (n .* t));
	    pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
	    hgt = [ 4  5   3   4  5  4.2 2.1 4.3  3.1 5.1 4.2];
	    wth = [.005 .005 .006 .01 .01 .03 .01 .01  .005 .008 .005];
	    for j =1:length(pos)
	       sig = sig + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
	    end
	elseif strcmp(Name,'Leopold'),
	    sig = (t == floor(.37 * n)/n);  % Kronecker
	elseif strcmp(Name,'Riemann'),
		sqn = round(sqrt(n));
	    sig = t .* 0;  % Riemann's Non-differentiable Function
		sig((1:sqn).^2) = 1. ./ (1:sqn);
		sig = real(ifft(sig));
	elseif strcmp(Name,'HypChirps'), % Hyperbolic Chirps of Mallat's book
		alpha	= 15*n*pi/1024;
		beta    = 5*n*pi/1024;
		t  	= (1.001:1:n+.001)./n;
		f1      = zeros(1,n);
		f2      = zeros(1,n);
		f1  	= sin(alpha./(.8-t)).*(0.1<t).*(t<0.68);
		f2  	= sin(beta./(.8-t)).*(0.1<t).*(t<0.75);
		M  	= round(0.65*n);
		P 	= floor(M/4);
		enveloppe = ones(1,M); % the rising cutoff function
            enveloppe(1:P) = (1+sin(-pi/2+((1:P)-ones(1,P))./(P-1)*pi))/2;
            enveloppe(M-P+1:M) = reverse(enveloppe(1:P));
  		env 	= zeros(1,n);
  		env(ceil(n/10):M+ceil(n/10)-1) = enveloppe(1:M);
		sig     = (f1+f2).*env;
	elseif strcmp(Name,'LinChirps'), % Linear Chirps of Mallat's book
		b 	= 100*n*pi/1024;
		a 	= 250*n*pi/1024;
		t 	= (1:n)./n;
		A1 	= sqrt((t-1/n).*(1-t));
		sig	= A1.*(cos((a*(t).^2)) + cos((b*t+a*(t).^2)));
	elseif strcmp(Name,'Chirps'), % Mixture of Chirps of Mallat's book
		t 	= (1:n)./n.*10.*pi;
  		f1 	= cos(t.^2*n/1024);
		a 	= 30*n/1024;
  		t 	= (1:n)./n.*pi;
  		f2 	= cos(a.*(t.^3));
  		f2 	= reverse(f2);
		ix 	= (-n:n)./n.*20;
 		g 	= exp(-ix.^2*4*n/1024);
		i1 	= (n/2+1:n/2+n);
		i2 	= (n/8+1:n/8+n);
		j  	= (1:n)/n;
    		f3 	= g(i1).*cos(50.*pi.*j*n/1024);
		f4 	= g(i2).*cos(350.*pi.*j*n/1024);
		sig 	= f1+f2+f3+f4;
   	 enveloppe = ones(1,n); % the rising cutoff function
  	 enveloppe(1:n/8) = (1+sin(-pi/2+((1:n/8)-ones(1,n/8))./(n/8-1)*pi))/2;
  	 enveloppe(7*n/8+1:n) = reverse(enveloppe(1:n/8));
 		sig 	= sig.*enveloppe;
        elseif strcmp(Name,'Gabor'), % two modulated Gabor functions in
				     % Mallat's book
		N = 512;
   		t = (-N:N)*5/N;
        	j = (1:N)./N;
		g = exp(-t.^2*20);
		i1 = (2*N/4+1:2*N/4+N);
		i2 = (N/4+1:N/4+N);
		sig1 = 3*g(i1).*exp(i*N/16.*pi.*j);
		sig2 = 3*g(i2).*exp(i*N/4.*pi.*j);
   		sig = sig1+sig2;
	elseif strcmp(Name,'sineoneoverx'), % sin(1/x) in Mallat's book
		N = 1024;
		i1 = (-N+1:N);
		i1(N) = 1/100;
		i1 = i1./(N-1);
		sig = sin(1.5./(i1));
		sig = sig(513:1536);
	elseif strcmp(Name,'Cusp2'),
		N = 64;
		i1 = (1:N)./N;
		x = (1-sqrt(i1)) + i1/2 -.5;
		M = 8*N;
		sig = zeros(1,M);
		sig(M-1.5.*N+1:M-.5*N) = x;
		sig(M-2.5*N+2:M-1.5.*N+1) = reverse(x);
		sig(3*N+1:3*N + N) = .5*ones(1,N);
	elseif strcmp(Name,'SmoothCusp'),
		sig = MakeSignal('Cusp2');
		N = 64;
		M = 8*N;
		t = (1:M)/M;
		sigma = 0.01;
		g = exp(-.5.*(abs(t-.5)./sigma).^2)./sigma./sqrt(2*pi);
		g = fftshift(g);
		sig2 = iconv(g',sig)'/M;
   	elseif strcmp(Name,'Piece-Regular'),
		sig1=-15*MakeSignal('Bumps',n);
		t = (1:fix(n/12)) ./fix(n/12);
		sig2=-exp(4*t);
		t = (1:fix(n/7)) ./fix(n/7);
		sig5=exp(4*t)-exp(4);
		t = (1:fix(n/3)) ./fix(n/3);
		sigma=6/40;
		sig6=-70*exp(-((t-1/2).*(t-1/2))/(2*sigma^2));
		sig(1:fix(n/7))= sig6(1:fix(n/7));
		sig((fix(n/7)+1):fix(n/5))=0.5*sig6((fix(n/7)+1):fix(n/5));
		sig((fix(n/5)+1):fix(n/3))=sig6((fix(n/5)+1):fix(n/3));
		sig((fix(n/3)+1):fix(n/2))=sig1((fix(n/3)+1):fix(n/2));
		sig((fix(n/2)+1):(fix(n/2)+fix(n/12)))=sig2;
		sig((fix(n/2)+2*fix(n/12)):-1:(fix(n/2)+fix(n/12)+1))=sig2;
sig(fix(n/2)+2*fix(n/12)+fix(n/20)+1:(fix(n/2)+2*fix(n/12)+3*fix(n/20)))=...
-ones(1,fix(n/2)+2*fix(n/12)+3*fix(n/20)-fix(n/2)-2*fix(n/12)-fix(n/20))*25;
		k=fix(n/2)+2*fix(n/12)+3*fix(n/20);
		sig((k+1):(k+fix(n/7)))=sig5;
		diff=n-5*fix(n/5);
		sig(5*fix(n/5)+1:n)=sig(diff:-1:1);
		% zero-mean
		bias=sum(sig)/n;
		sig=bias-sig;
   	elseif strcmp(Name,'Piece-Polynomial'),
		t = (1:fix(n/5)) ./fix(n/5);
		sig1=20*(t.^3+t.^2+4);
		sig3=40*(2.*t.^3+t) + 100;
		sig2=10.*t.^3 + 45;
		sig4=16*t.^2+8.*t+16;
		sig5=20*(t+4);
		sig6(1:fix(n/10))=ones(1,fix(n/10));
		sig6=sig6*20;
		sig(1:fix(n/5))=sig1;
		sig(2*fix(n/5):-1:(fix(n/5)+1))=sig2;
		sig((2*fix(n/5)+1):3*fix(n/5))=sig3;
		sig((3*fix(n/5)+1):4*fix(n/5))=sig4;
		sig((4*fix(n/5)+1):5*fix(n/5))=sig5(fix(n/5):-1:1);
		diff=n-5*fix(n/5);
		sig(5*fix(n/5)+1:n)=sig(diff:-1:1);
		%sig((fix(n/20)+1):(fix(n/20)+fix(n/10)))=-ones(1,fix(n/10))*20;
		sig((fix(n/20)+1):(fix(n/20)+fix(n/10)))=ones(1,fix(n/10))*10;
		sig((n-fix(n/10)+1):(n+fix(n/20)-fix(n/10)))=ones(1,fix(n/20))*150;
		% zero-mean
		bias=sum(sig)/n;
		sig=sig-bias;
   	elseif strcmp(Name,'Gaussian'),
		sig=GWN(n,beta);
		g=zeros(1,n);
		lim=alpha*n;
		mult=pi/(2*alpha*n);
		g(1:lim)=(cos(mult*(1:lim))).^2;
		g((n/2+1):n)=g((n/2):-1:1);
		g = rnshift(g,n/2);
		g=g/norm(g);
		sig=iconv(g,sig);
       else
	    disp(sprintf('MakeSignal: I don*t recognize <<%s>>',Name))
	    disp('Allowable Names are:')
	       disp('HeaviSine'),
	       disp('Bumps'),
	       disp('Blocks'),
	       disp('Doppler'),
	       disp('Ramp'),
	       disp('Cusp'),
	       disp('Crease'),
	       disp('Sing'),
	       disp('HiSine'),
	       disp('LoSine'),
	       disp('LinChirp'),
	       disp('TwoChirp'),
	       disp('QuadChirp'),
	       disp('MishMash'),
	       disp('WernerSorrows'),
	       disp('Leopold'),
	       disp('Sing'),
	       disp('HiSine'),
	       disp('LoSine'),
	       disp('LinChirp'),
	       disp('TwoChirp'),
	       disp('QuadChirp'),
	       disp('MishMash'),
	       disp('WernerSorrows'),
	       disp('Leopold'),
	       disp('Riemann'),
	       disp('HypChirps'),
	       disp('LinChirps'),
	       disp('Chirps'),
	       disp('sineoneoverx'),
	       disp('Cusp2'),
	       disp('SmoothCusp'),
	       disp('Gabor'),
	       disp('Piece-Regular');
	       disp('Piece-Polynomial');
	       disp('Gaussian');
	end

%
% Originally made by David L. Donoho.
% Function has been enhanced.

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:39 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner


text:
function [OutputSurfaces, OutputMris] = bst_warp(destPts, srcPts, SurfaceFiles, MriFiles, OutputTag, OutputDir, isSurfaceOnly)
%BST_WARP:  Deform an anatomy (MRI+surfaces) to fit a set of landmarks.
%
% USAGE:  bst_warp(destPts, srcPts, SurfaceFiles, MriFiles, OutputTag, OutputDir, isSurfaceOnly=0)
%
% INPUTS:
%     - destPts      : landmarks in real head coordinates, i.e. from a Polhemus or
%     - srcPts       : landmarks in intial anatomy coordinates
%     - SurfaceFiles : Cell array of the surface files to be warped (full path)
%     - MriFiles     : Cell array of the MRI files to be warped (full path)
%     - OutputTag    : Tag to add at the end of the filenames of the wrapped files
%     - OutputDir    : Output directory
%     - isSurfaceOnly: If 1, do not warp the MRI

% @=============================================================================
% 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: Felix Darvas, 2005
%          Louis Hovasse, 2009
%          Francois Tadel, 2010-2022

if (nargin < 7) || isempty(isSurfaceOnly)
    isSurfaceOnly = 0;
end
if ischar(MriFiles)
    MriFiles = {MriFiles};
end
OutputSurfaces = {};
OutputMris = {};

%% ===== COMPUTE WARP PARAMETER =====
% Compute warp transform parameters
[W,A] = warp_transform(srcPts, destPts);

%% ===== WARP SURFACES =====
% Progress bar
bst_progress('start', 'Warp anatomy', 'Warping surfaces...', 0, length(SurfaceFiles));
% Loop on all surface files
for i = 1:length(SurfaceFiles)
    bst_progress('inc', 1);
    % Read input
    sSurf = load(file_fullpath(SurfaceFiles{i}));
    % Warp surface
    switch file_gettype(SurfaceFiles{i})
        case 'fibers'
            nRows = size(sSurf.Points, 1);
            Points = permute(reshape(sSurf.Points, [], 1, 3), [1 3 2]);
            Points = warp_lm(Points, A, W, srcPts) + Points;
            Points = reshape(permute(Points, [1 3 2]), nRows, [], 3);
            sSurfNew.Points = Points;
            sSurfNew.Header = sSurf.Header;
            sSurfNew.Colors = sSurf.Colors;
            sSurfNew.Scouts = sSurf.Scouts;
        case 'fem'
            sSurfNew.Elements     = sSurf.Elements;
            sSurfNew.Tissue       = sSurf.Tissue;
            sSurfNew.TissueLabels = sSurf.TissueLabels;
            sSurfNew.Tensors      = sSurf.Tensors;
            sSurfNew.Vertices     = warp_lm(sSurf.Vertices, A, W, srcPts) + sSurf.Vertices;
        otherwise
            sSurfNew.Faces    = sSurf.Faces;
            sSurfNew.Vertices = warp_lm(sSurf.Vertices, A, W, srcPts) + sSurf.Vertices;
    end
    % Add tag to comment
    sSurfNew.Comment  = [sSurf.Comment ' warped'];
    % Copy previous field
    if isfield(sSurf, 'Atlas') && ~isempty(sSurf.Atlas)
        sSurfNew.Atlas = sSurf.Atlas;
    end
    if isfield(sSurf, 'Reg') && ~isempty(sSurf.Reg)
        sSurfNew.Reg = sSurf.Reg;
    end
    if isfield(sSurf, 'History') && ~isempty(sSurf.History)
        sSurfNew.History = sSurf.History;
    end
    % History: Warp
    sSurfNew = bst_history('add', sSurfNew, 'warp', 'Surface deformed to match the head points from channel file.');
    % Output filename
    [tmp__, fileBase, fileExt] = bst_fileparts(SurfaceFiles{i});
    OutputSurfaces{end+1} = bst_fullfile(OutputDir, [fileBase, OutputTag, fileExt]);
    % Save new file
    bst_save(OutputSurfaces{end}, sSurfNew, 'v7');
end

%% ====== COMPUTE TRANSFORMATION ========
bst_progress('start', 'Warp anatomy', 'Preparing MRI...');
% Load reference MRI (first in the list)
sMriSrc = in_mri_bst(MriFiles{1});
% Transform landmarks into MRI coordinates (meters)
srcPts_mr  = cs_convert(sMriSrc, 'scs', 'voxel', srcPts);
destPts_mr = cs_convert(sMriSrc, 'scs', 'voxel', destPts);
% Compute warp transform in MR coordinates
[Wmr,Amr] = warp_transform(srcPts_mr, destPts_mr);
% Compute "inverse" warp transform in MR coordinates
[Wmr_inv,Amr_inv] = warp_transform(destPts_mr, srcPts_mr);

%% ===== WRAP FIDUCIALS =====
% % SCS transformation
% destSCS.NAS = warp_lm(sMriSrc.SCS.NAS, Amr, Wmr, srcPts_mr) + sMriSrc.SCS.NAS;
% destSCS.LPA = warp_lm(sMriSrc.SCS.LPA, Amr, Wmr, srcPts_mr) + sMriSrc.SCS.LPA;
% destSCS.RPA = warp_lm(sMriSrc.SCS.RPA, Amr, Wmr, srcPts_mr) + sMriSrc.SCS.RPA;
% NCS transformation
if isfield(sMriSrc, 'NCS') && isfield(sMriSrc.NCS, 'AC') && ~isempty(sMriSrc.NCS.AC) && ~isempty(sMriSrc.NCS.PC) && ~isempty(sMriSrc.NCS.IH)
    destNCS.AC = warp_lm(sMriSrc.NCS.AC, Amr, Wmr, srcPts_mr) + sMriSrc.NCS.AC;
    destNCS.PC = warp_lm(sMriSrc.NCS.PC, Amr, Wmr, srcPts_mr) + sMriSrc.NCS.PC;
    destNCS.IH = warp_lm(sMriSrc.NCS.IH, Amr, Wmr, srcPts_mr) + sMriSrc.NCS.IH;
    destNCS.R = [];
    destNCS.T = [];
    destNCS.Origin = [];
else
    destNCS = [];
end

%% ===== WARP MRI =====
bst_progress('start', 'Warp anatomy', 'Warping MRI...', 0, 100*length(MriFiles));
for iMri = 1:length(MriFiles)
    bst_progress('text', sprintf('Warping MRI...   [%d/%d]', iMri, length(MriFiles)));
    % Load MRI (the first one is already loaded above)
    if (iMri >= 2)
        sMriSrc = in_mri_bst(MriFiles{iMri});
    end
    % If warping of the MRI volumes is enabled
    if ~isSurfaceOnly
        % Process coordinates by blocks: Doing all at once costs too much memory, doing only 1 at a time costs too much time
        sizeMri = size(sMriSrc.Cube);
        if (length(sizeMri) > 3)
            error('No support for 4D volumes. Ask on the Brainstorm for help.');
        end
        newCube = ones(sizeMri);
        nVoxels = numel(newCube);
        BLOCK_SIZE = 10000;
        nBlocks = ceil(nVoxels / BLOCK_SIZE);
        ix0 = 1;
        for i = 1:nBlocks
            % Increment progress bar
            if (mod(i, round(nBlocks/100)) == 0)
                bst_progress('inc', 1);
            end
            % Get indices in dest volume
            ix1 = min(ix0 - 1 + BLOCK_SIZE, nVoxels);
            [xv,yv,zv] = ind2sub(sizeMri, ix0:ix1);
            rv = [xv;yv;zv]';
            % Unwarp MRI coordinates
            rv_inv = warp_lm(rv, Amr_inv, Wmr_inv, destPts_mr) + rv;
            % Round coordinates (nearest neighor interpolation)
            rv_inv = round(rv_inv);
            % Remove values that are outside the volume
            iOutside = find(sum((rv_inv < 1) | (rv_inv > repmat(sizeMri,size(rv_inv,1),1)),2) > 0);
            rv_inv(iOutside,:) = 1;
            % Get indices from xyz coordinates
            ix_inv = sub2ind(sizeMri, rv_inv(:,1), rv_inv(:,2), rv_inv(:,3));
            % Get values in initial volume
            newCube(ix0:ix1) = sMriSrc.Cube(ix_inv);
            % Set values outside of the volume to zero
            newCube(iOutside) = 0;
            % Go to next block
            ix0 = ix1 + 1;
        end
        newComment = [sMriSrc.Comment, ' warped'];
    else
        newCube = sMriSrc.Cube;
        newComment = sMriSrc.Comment;
    end

    % === SAVE NEW MRI ===
    % Create new structure
    sMriDest = sMriSrc;
    sMriDest.Cube    = newCube;
    sMriDest.Comment = newComment;
    sMriDest.NCS     = destNCS;
    % History: Copy previous field
    if isfield(sMriSrc, 'History') && ~isempty(sMriSrc.History)
        sMriDest.History = sMriSrc.History;
    end
    % History: Warp
    sMriDest = bst_history('add', sMriDest, 'warp', 'MRI deformed to match the head points from channel file.');
    % Output filename
    [tmp__, fileBase, fileExt] = bst_fileparts(MriFiles{iMri});
    OutputMris{iMri} = bst_fullfile(OutputDir, [fileBase, OutputTag, fileExt]);
    % Save new MRI
    bst_save(OutputMris{iMri}, sMriDest, 'v7');
end

bst_progress('stop');

end

%% =====================================================================================
%  ===== HELPER FUNCTIONS ==============================================================
%  =====================================================================================

%% ===== WARP TRANSFORM =====
% Calculates nonlinear transformation coefficents (see Ermer's Thesis)
% INPUT:
%    - p : Landmarks in system 1
%    - q : Landmarks in system 2
% OUTPUT:
%    - e : Warp energy
function [W,A,e] = warp_transform(p, q)
    N = size(p,1);
    px = repmat(p(:,1), 1, N);
    py = repmat(p(:,2), 1, N);
    pz = repmat(p(:,3), 1, N);
    K = sqrt((px - px').^2 + (py - py').^2 + (pz - pz').^2);

    P = [p, ones(N,1)];
    L = [K P; P' zeros(4,4)];
    D = [q - p; zeros(4,3)];
    warning off
    H = L \ D;
    warning on
    if any(isnan(H))
        H = pinv(L) * D;
    end
    W = H(1:N,:);
    A = H(N+1:end, :);
    e = sum(diag(W' * K * W));
end

%% ===== WARP LANDMARKS: VECTORIZED =====
% Performs warp transformation with linear 3D RFB (see Ermer's Thesis)
function rw = warp_lm(r, A, W, p)
    rw = r * A(1:3,1:3);
    rw = bst_bsxfun(@plus, rw, A(4,:));
    np = size(p,1);
    U = sqrt(bst_bsxfun(@minus, repmat(r(:,1),1,np), p(:,1)') .^ 2 + ...
             bst_bsxfun(@minus, repmat(r(:,2),1,np), p(:,2)') .^ 2 + ...
             bst_bsxfun(@minus, repmat(r(:,3),1,np), p(:,3)') .^ 2);
    rw = rw + U * W;
end

%% ===== WARP LANDMARKS: LOOP VERSION =====
% This function is twice slower
% function rw = warp_lm(r, A, W, p)
%     rw = r * A(1:3,1:3) + repmat(A(4,:), size(r,1), 1);
%     for i = 1:size(p,1)
%         U = sqrt(sum((r - repmat(p(i,:), size(r,1), 1)) .^ 2, 2));
%         rw = rw + U * W(i,:);
%     end
% end

author: brainstorm-tools


text:
function count = persistent_variables(n)
    if nargin == 0
        n = 1;
    end

    persistent acc;
    if isempty(acc)
        acc = 1;
    else
        acc = acc + n;
    end
    count = acc;
end

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function sparse_count_test08 ( dim_min, dim_max, level_max_min, ...
  level_max_max )

%*****************************************************************************80
%
%% SPARSE_COUNT_TEST08 tests OWN_E_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_MIN, the minimum spatial dimension to consider.
%
%    Input, integer DIM_MAX, the maximum spatial dimension to consider.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX to consider.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX to consider.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_COUNT_TEST08\n' );
  fprintf ( 1, '  OWN_E_SIZE returns the number of\n' );
  fprintf ( 1, '  distinct points in an OWN_E sparse grid made from \n' );
  fprintf ( 1, '  product grids formed from open weakly nested \n' );
  fprintf ( 1, '  quadrature rules with exponential growth, including:\n' );
  fprintf ( 1, '  * GGH_E, the Generalized Gauss-Hermite Exponential Growth Family;\n' );
  fprintf ( 1, '  * GH_E, the Gauss-Hermite Exponential Growth Family;\n' );
  fprintf ( 1, '  * LG_E, the Gauss-Legendre Exponential Growth Family;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   DIM: ' );

  for dim_num = dim_min : dim_max
    fprintf ( 1, '  %10d', dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   LEVEL_MAX\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max
    fprintf ( 1, '    %4d', level_max );
    for dim_num = dim_min : dim_max
      point_num = own_e_size ( dim_num, level_max );
      fprintf ( 1, '  %10d', point_num );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [blocks,idx] = my_im2col(I,blkSize,slidingDis)
if (slidingDis==1)
    blocks = im2col(I,blkSize,'sliding');
    idx = [1:size(blocks,2)];
    return
end

idxMat = zeros(size(I)-blkSize+1);
idxMat([[1:slidingDis:end-1],end],[[1:slidingDis:end-1],end]) = 1; % take blocks in distances of 'slidingDix', but always take the first and last one (in each row and column).
idx = find(idxMat);
[rows,cols] = ind2sub(size(idxMat),idx);
blocks = zeros(prod(blkSize),length(idx));
for i = 1:length(idx)
    currBlock = I(rows(i):rows(i)+blkSize(1)-1,cols(i):cols(i)+blkSize(2)-1);
    blocks(:,i) = currBlock(:);
end

author: thfylsty


text:
function graphics ( it, x, q, q_hat_plot, u, u_hat_plot )

%*****************************************************************************80
%
%% GRAPHICS displays the computed and optimal functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2011
%
%  Author:
%
%    Jeff Borggaard, John Burkardt, Catalin Trenchea, Clayton Webster
%
%  Parameters:
%
%    Input, integer IT, the iteration number.
%
%    Input, real X(*), the node coordinates.
%
%    Input, real Q(*), the computed control at each node.
%
%    Input, real Q_HAT_PLOT(*), the optimal control at each node.
%
%    Input, real U(*), the computed BVP solution at each node.
%
%    Input, real U_HAT_PLOT(*), the optimal BVP solution at each node.
%
  figure ( 1 )
  clf
  hold on
  plot ( x, u )
  plot ( x, u_hat_plot, 'r+' )
  grid
  title ( sprintf ( '  Target function U\\_hat (red), computed U (blue), iteration %d', it ) )
  hold off

  figure ( 2 )
  clf
  hold on
  plot ( x, q )
  plot ( x, q_hat_plot, 'r+' )
  grid
  title ( sprintf ( '  Optimal control Q\\_hat (red), computed Q (blue), iteration %d', it ) )
  hold off

  return
end


text:
function [elem] = subsref(tt,s)
% SUBSREF V2
% Evaluate elements, cores of TT-formats and fields of the TT-structure
%   A=TT{I} computes I-th core of the TT-representation
%
%   ELEM=TT(IND) computes element with index IND of the tensor TT
%
% TT-Toolbox 2.2.1, 2009-2016
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com

% Possibility to extract subtensor added by Alexey Boyko,
% Skolkovo Institute of Science and Technology
% For questions email alexey.boyko@skolkovotech.ru or a.boyko@skoltech.ru

% Example 1
% y =tt_x(2,10)
% y2=y(1,2,1,2,2,:,:,:,:,:), or, alternatively,
% y2=y(1,2,1,2,2) or
% y2=y([1,2,1,2,2]) or
% y2=y({1,2,1,2,2,':',':',':',':',':'}) or
% y2=y({1,2,1,2,2}) or
% will return 5-dimensional tt_tensor, corresponding to
% fist 5 modes\harmonics being frozen to indices 1,2,1,2,2 respectively

% Example 2
% y =tt_x(2,10)
% y2=y(1,:,1,2,2,:,:,1,:,:), or, alternatively,
% y2=y({1,':',1,2,2,':',':',1}) or
% y2=y({1,':',1,2,2,':',':',1,':',':'}) or

% will also return 5-dimensional tt_tensor, corresponding to freezing
% respective modes\harmonics

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

switch s(1).type
    case '()'
    %Evaluate element of a tensor in position s
    d=tt.d; n=tt.n; ps=tt.ps; cr=tt.core; r=tt.r;

    %ind=double(s);
    pp=s.subs;
    mn=numel(pp);

    if mn == 1
        if iscell( pp{1})
            %one argument, cell array
            pp=pp{1};
        else
            % one argument, ordinary array
            pp=num2cell(pp{1});
        end
    end
    mn=numel(pp);

    if ( mn > d )
        error('Invalid number of indices specified: given %d, need no more than %d \n',mn,d);
    end

    idxs=[];
    vals=[];
    for i=1:mn
        if isnumeric(pp{i})
            idxs=[idxs,i];
            vals=[vals,pp{i}];
        end
    end

    mn=numel(idxs);
    if(mn < d)
        elem=tt_subtensor(tt,idxs,vals);
    else

        %               if ( mn == 1 ) %Special case
        %                  ind=pp{1};
        %                  if ( is_array(ind) )
        %                    for i=1:d
        %                      pp{i}=ind(i);
        %                    end
        %                  else
        %                     pp=ind; %Supposedly it was a cell array of index subsets
        %                  end
        %               end

        elem=tt_tensor;
        elem.r=r;
        elem.d=d;
        n0=zeros(d,1);
        for i=1:d
            ind=pp{i};
            if ( strcmp(ind, ':'))
                pp{i}=1:n(i);
                ind=pp{i};
            end
            n0(i)=numel(ind);
        end
        ps1=cumsum([1;n0.*r(1:d).*r(2:d+1)]);
        cr1=zeros(ps1(d+1)-1,1);

        for i=1:d
            ind=pp{i};
            crx=cr(ps(i):ps(i+1)-1);
            crx=reshape(crx,[r(i),n(i),r(i+1)]);
            crx=crx(:,ind,:);
            cr1(ps1(i):ps1(i+1)-1)=crx(:);
        end

        if ( prod(n0) == 1 ) %single element
            v=1;
            for i=1:d
                crx=cr1(ps1(i):ps1(i+1)-1); crx=reshape(crx,[r(i),r(i+1)]);
                v=v*crx;
            end
            elem=v;
        else
            elem.n=n0;
            elem.core=cr1;
            elem.ps=ps1;
        end

    end

    case '.'
        switch s(1).subs
            case 'r'
                elem = tt.r;
                if (numel(s)>1)
                    s = s(2:end);
                    elem = subsref(elem, s);
                end;
%                 if (numel(s)>1)&&(strcmp(s(2).type,'()'))
%                     elem = tt.r(s(2).subs{1});
%                 else
%                     elem=tt.r;
%                 end;
            case 'over'
                elem = tt.over;
                if (numel(s)>1)
                    s = s(2:end);
                    elem = subsref(elem, s);
                end;
            case 'core'
                if (numel(s)>1)&&(strcmp(s(2).type,'()'))
                    elem = tt.core(s(2).subs{1});
                else
                    elem=tt.core;
                end;
            case 'd'
                elem=tt.d;
            case 'ps'
                elem = tt.ps;
                if (numel(s)>1)
                    s = s(2:end);
                    elem = subsref(elem, s);
                end;

%                 if (numel(s)>1)&&(strcmp(s(2).type,'()'))
%                     elem = tt.ps(s(2).subs{1});
%                 else
%                     elem=tt.ps;
%                 end;
            case 'n'
                elem = tt.n;
                if (numel(s)>1)
                    s = s(2:end);
                    elem = subsref(elem, s);
                end;

%                 if (numel(s)>1)&&(strcmp(s(2).type,'()'))
%                     elem = tt.n(s(2).subs{1});
%                 else
%                     elem=tt.n;
%                 end;
            otherwise
                error(['No field ', s.subs, ' is here.']);
        end

    case '{}'
        %Return the core in the old (not exactly!!! r1-n-r2 here) format
        pp=s.subs;
        mn=numel(pp);
        if ( mn > 1 )
          error('Invalid number of cores asked');
        end
        elem=core(tt,pp{1});

        if (numel(s)>1)
            s = s(2:end);
            elem = subsref(elem, s);
        end;
%         if (pp{1}~=1)
%             elem=permute(elem,[2,1,3]);
%         end;

    otherwise
        error('Invalid subsref.');
end
end

author: oseledets

repo: TT-Toolbox


text:
function c = mtexdegchar
% returns the degree character

c = getMTEXpref('degreeChar');


text:
function combineKeypoints(expidx)

fprintf('combineKeypoints()\n');

if (ischar(expidx))
    expidx = str2double(expidx);
end

% load annolist with full size images
[~,imgidxs] = getAnnolist(expidx);

p = rcnn_exp_params(expidx);
saveTo = [fileparts(p.evalTest) '/pred_backproj'];

fname = [saveTo '/keypoints_' padZeros(num2str(imgidxs(1)),5)];
load(fname,'keypoints');
keypointsAll = keypoints;

fname = [saveTo '/bboxes_' padZeros(num2str(imgidxs(1)),5)];
load(fname,'bboxes');
bboxAll = bboxes;

for i = 1:length(imgidxs)
    fprintf('.');

    fname = [saveTo '/keypoints_' padZeros(num2str(imgidxs(i)),5)];
    load(fname,'keypoints');
    keypointsAll(i) = keypoints;

    fname = [saveTo '/bboxes_' padZeros(num2str(imgidxs(i)),5)];
    load(fname,'bboxes');
    bboxAll(i) = bboxes;

    if (~mod(i, 100))
        fprintf(' %d/%d\n',i,length(imgidxs));
    end
end
fprintf(' done\n');

fnameKeypoints = [fileparts(p.evalTest) '/keypointsAll'];
save(fnameKeypoints, 'keypointsAll');

fnameBBox = [fileparts(p.evalTest) '/bboxAll'];
save(fnameBBox, 'bboxAll');
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
classdef m_47_IHM19_16p_4s < MARRMoT_model
% Class for hydrologic conceptual model: Forellenbach model

% Copyright (C) 2021 Clara Brandes, 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.
%
% Model reference
% Brandes, C. (2020). Erstellung eines Konzeptionellen
% Hochwasserabfluss-modells für das Einzugsgebiet des Forellenbachs, NP
% Bayerischer Wald. MSc Thesis. Technische Univerisität Dresden, Germany.

    properties
        % model-specific attributes
    end
    methods

        % creator method
        function obj = m_47_IHM19_16p_4s()
            obj.numStores = 4;                                             % number of model stores
            obj.numFluxes = 18;                                            % number of model fluxes
            obj.numParams = 16;

            obj.JacobPattern  = [1,0,0,0;
                                 1,1,0,0;
                                 1,1,1,0;
                                 0,1,1,1];                                 % Jacobian matrix of model store ODEs

            obj.parRanges = [   2,     5;      % 01 SIMAX, maximum interception storage [mm]
                                0.9,   1;      % 02 A, splitting coeffcient for excess precipitation [-]
                                0.4,   0.95;   % 03 FF, forest fraction [-]
                                0.05,  5;      % 04 SMPMAX, maximum storage macropores [mm]
                                0,     1;      % 05 CQMP, runoff time parameter (fast/slow runnoff) first soil layer [1/d]
                                1,     5;      % 06 XQMP, runoff scale parameter first soil layer [-]
                              400,   600;      % 07 SS1MAX, maximum soil moisture storage first soil layer [mm]
                                0.3,   0.7;    % 08 FCCS1, field capacity as fraction of maximum storage first soil layer [-]
                                0,  1000;      % 09 CFS1, maximum infiltration rate first soil layer [mm/d]
                                0,    15;      % 10 XFS1, infiltration loss exponent first soil layer [-]
                                0,     1;      % 11 CQS1, runoff time parameter for (fast/slow runnoff) first soil layer [1/d]
                                1,     5;      % 12 XQS1, runoff scale parameter first soil layer [-]
                              300,   500;      % 13 SS2MAX, maximum soil moisture storage second soil layer [mm]
                                0,     1;      % 14 CQS2, runoff time parameter for (fast/slow runnoff) second soil layer [1/d]
                                1,     5;      % 15 XQS2, runoff scale parameter second soil layer [-]
                                0.01,  5];     % 16 D, Flow delay before surface runoff [d]

            obj.StoreNames = {"S1", "S2" "S3" "S4"};                        % Names for the stores
            obj.FluxNames  = {"ei",    "pex", "pexmp",  "pexs1", "fmp",...
                              "qexmp", "qmp", "pqexsl", "fs1",   "etas1",...
                              "qs1",   "q0",  "q0r",    "qmps1", "pc",...
                              "qh",    "qs2", "qgw"};                      % Names for the fluxes

            obj.FluxGroups.Ea = [1 10];                                    % Index or indices of fluxes to add to Actual ET
            obj.FluxGroups.Q  = [13 16 17 18];                             % Index or indices of fluxes to add to Streamflow
            obj.FluxGroups.GW = -18;                                       % Index of GW runoff flow.

        end

        % INITialisation function
        function obj = init(obj)
            % parameters
            theta   = obj.theta;
            delta_t = obj.delta_t;
            SIMAX   = theta(1);                         % maximum interception storage [mm]
            SMPMAX  = theta(4);                         % maximum storage macropores [mm]
            SS1MAX  = theta(7);                         % maximum soil moisture storage first soil layer [mm]
            SS2MAX  = theta(13);                        % maximum soil moisture storage second soil layer [mm]
            D       = theta(16);                        % Flow delay before surface runoff [d]

            % initial store values
            obj.S0 = 0.8 * [SIMAX;                   % Initial interception storage
                            SMPMAX;                  % Initial macropore storage
                            SS1MAX;                  % Initial soil moisture storage 1
                            SS2MAX];                 % Initial soil moisture storage 2

            % initialise the unit hydrographs and still-to-flow vectors
            uh_q0r = uh_4_full(D,delta_t);

            obj.uhs = {uh_q0r};
        end

        % MODEL_FUN are the model governing equations in state-space formulation
        function [dS, fluxes] = model_fun(obj, S)
            % parameters
            theta   = obj.theta;
            SIMAX   = theta(1);                         % maximum interception storage [mm]
            A       = theta(2);                         % splitting coeffcient for excess precipitation [-]
            FF      = theta(3);                         % forest fraction [-]
            SMPMAX  = theta(4);                         % maximum storage macropores [mm]
            CQMP    = theta(5);                         % runoff time parameter (fast/slow runnoff) first soil layer [1/d]
            XQMP    = theta(6);                         % runoff scale parameter first soil layer [-]
            SS1MAX  = theta(7);                         % maximum soil moisture storage first soil layer [mm]
            FCCS1   = theta(8);                         % field capacity coefficient fist soil layer [-]
            CFS1    = theta(9);                         % maximum infiltration rate first soil layer [-]
            XFS1    = theta(10);                        % infiltration loss exponent first soil layer [-]
            CQS1    = theta(11);                        % runoff time parameter for (fast/slow runnoff) first soil layer [1/d]
            XQS1    = theta(12);                        % runoff scale parameter first soil layer [-]
            SS2MAX  = theta(13);                        % maximum soil moisture storage second soil layer [mm]
            CQS2    = theta(14);                        % runoff time parameter for (fast/slow runnoff) second soil layer [1/d]
            XQS2    = theta(15);                        % runoff scale parameter second soil layer [-]

            % delta_t
            delta_t = obj.delta_t;

            % unit hydrographs and still-to-flow vectors
            uhs = obj.uhs;
            uh_q0r = uhs{1};

            % stores
            S1 = S(1);
            S2 = S(2);
            S3 = S(3);
            S4 = S(4);

            % climate input
            t = obj.t;                             % this time step
            climate_in = obj.input_climate(t,:);   % climate at this step
            P  = climate_in(1);
            Ep = climate_in(2);
            T  = climate_in(3);

            % fluxes functions
            flux_ei      = evap_1(S1,Ep,delta_t);
            flux_pex     = interception_1(P,S1,SIMAX);
            flux_pexmp   = split_1(A,flux_pex);
            flux_pexs1   = split_2(A,flux_pex);
            flux_fmp     = infiltration_3(flux_pexmp,S2,SMPMAX);
            flux_qexmp   = flux_pexmp - flux_fmp;
            flux_qmp     = interflow_3(CQMP,XQMP,S2,delta_t);
            flux_pqexs1  = flux_pexs1 + flux_qexmp;
            flux_fs1     = infiltration_7(CFS1,XFS1,S3,SS1MAX,flux_pqexs1);
            flux_etas1   = evap_23(FF,FCCS1,S3,SS1MAX,Ep,delta_t);
            flux_qs1     = interflow_9(S3,CQS1,FCCS1*SS1MAX,XQS1,delta_t);
            flux_q0      = flux_pqexs1 - flux_fs1;
            flux_q0r     = route(flux_q0, uh_q0r);
            flux_qmps1   = flux_qmp + flux_qs1;
            flux_pc      = infiltration_3(flux_qmps1,S4,SS2MAX);
            flux_qh      = flux_qmps1 - flux_pc;
            flux_qs2     = interflow_3(CQS2,XQS2,S4,delta_t);
            flux_qgw     = 0.0195;

            % stores ODEs
            dS1 = P - flux_ei - flux_pex;
            dS2 = flux_fmp - flux_qmp;
            dS3 = flux_fs1 - flux_etas1 - flux_qs1;
            dS4 = flux_pc - flux_qs2;

            % outputs
            dS = [dS1 dS2 dS3 dS4];
            fluxes = [flux_ei, flux_pex, flux_pexmp, flux_pexs1,...
                      flux_fmp, flux_qexmp, flux_qmp, flux_pqexs1,...
                      flux_fs1, flux_etas1, flux_qs1, flux_q0,...
                      flux_q0r, flux_qmps1, flux_pc , flux_qh,...
                      flux_qs2, flux_qgw];
        end

        % STEP runs at the end of every timestep.
        function obj = step(obj)
            % unit hydrographs and still-to-flow vectors
            uhs = obj.uhs;
            uh_q0r = uhs{1};

            % input fluxes to the unit hydrographs
            fluxes = obj.fluxes(obj.t,:);
            flux_q0   = fluxes(12);

            % update still-to-flow vectors using fluxes at current step and
            % unit hydrographs
            uh_q0r = update_uh(uh_q0r, flux_q0);
            obj.uhs = {uh_q0r};
        end
    end
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function sParams = rmStimTemplate(vw, nScans)
% Defines the retinotopy model stimulus slot
%
%    sParams = rmStimTemplate(vw, nScans)
%
% This routine defines all of the parameters in a rm stimulus.  If you ever
% change the default, you must edit this function

n = viewGet(vw,'nscans');
if notDefined('nScans'), nScans = n; end

sParams = struct;

for ii=1:nScans

    sParams(ii).stimType   = '8Bars';
    sParams(ii).stimSize   =  14;
    sParams(ii).stimWidth  =  45;
    sParams(ii).stimStart  =  0;
    sParams(ii).stimDir    =  0;
    sParams(ii).nCycles    =  6;
    sParams(ii).nStimOnOff =  0;
    sParams(ii).nUniqueRep =  1;
    sParams(ii).prescanDuration = 8;
    sParams(ii).nDCT = 0;
    sParams(ii).hrfType =  'one gamma (Boynton style)';  % 'two gammas (SPM style)'
    sParams(ii).hrfParams = {[1.68 3 2.05],[5.4 5.2 10.8 7.35 0.35]};

    %If we are within the number scans in the dataTYPES, get frame
    %parameters from the view struct
    if ii <= n,
        sParams(ii).framePeriod = viewGet(vw,'framePeriod',ii);
        sParams(ii).nFrames     = viewGet(vw,'nFrames',ii);
    % Otherwise copy them from the previous scan
    else
        sParams(ii).framePeriod = sParams(ii-1).framePeriod;
        sParams(ii).nFrames     = sParams(ii-1).nFrames;
    end

    sParams(ii).fliprotate = [0 0 0];
    sParams(ii).imFile      = 'None';
    sParams(ii).jitterFile  = 'None';
    sParams(ii).paramsFile  = 'None';
    sParams(ii).imFilter    = 'None';

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Gs = lp_trfia(freq,A,B,C,D,E)
%
%  Computes the transfer function for systems
%      .
%    E*x  =  A*x + B*u
%      y  =  C*x + D*u
%
%  on the imaginary axis (more precisely, on the points
%  sqrt(-1)*freq(i), where i = 1,...,length(freq) ). This routine can only
%  be applied to systems, where the matrices are given explicitely.
%
%  Calling sequence:
%
%    Gs = lp_trfia(freq,A,B,C,D,E)
%
%  Input:
%
%    freq      (row) vector containing frequency points;
%    A         system matrix A (n-x-n matrix);
%    B         system matrix B (n-x-m matrix);
%    C         system matrix C (q-x-n matrix);
%    D         system matrix D (q-x-m matrix);
%              Set D = [] if D is not existing or the zero matrix!
%    E         system matrix E. Set E = [] if E is not existing or the
%              (sparse) identity matrix!
%
%  Output:
%
%    Gs        transfer function sampling matrix, which is a
%              q*m-x-length(freq) matrix);
%              Gs(:,i) contains the stacked columns of the matrix
%              D + C*(sqrt(-1)*freq(i)*E-A)^(-1)*B.
%
%  Remarks:
%
%    The vector freq can easily be computed by the function 'lp_lpfrq'.
%
%    The permutation of the matrices A and/or E before calling this
%    routine is not necessary.
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

% Input data not completely checked!

na = nargin;

if length(A)==0,
  Gs = [];
  return;
end

if na<6, E = []; end
with_E = length(E)>0;
with_D = length(D)>0;

[n,m] = size(B);
q = size(C,1);
nop = length(freq);

if ~with_D
  D = zeros(q,m);
end

Gs = zeros(m*q,nop);

j = sqrt(-1);

for i = 1:nop
  if with_E
    G0 = D + C*(((j*freq(i))*E-A)\B);
  else
    G0 = D + C*(((j*freq(i))*speye(n)-A)\B);
  end
  for k = 1:m
    Gs((k-1)*q+1:k*q,i) = G0(:,k);
  end
end


text:
% Script: call_rcgrid.m
% Script to determine rate changes in aftershock sequences on a grid with calc_rcgrid.m.
% Plots the first time slice at the end.
% The function works on the catalog newt2!!!
%
% Output variables:
% RCREL : Relative rate change matrix
%
% J.Woessner
% last update: 03.07.03

% Get input parameters
prompt  = {'Enter forecast period (number of events):','Enter maximum time of learning period (days):',...
    'Enter timesteps (days):','Longitude spacing / [deg]:','Latitude spacing / [deg]:','Radius / [deg]:'};
title   = 'Parameters ';
lines= 1;
def     = {'50','100','5','0.1','0.1','0.15'};
answer  = inputdlg(prompt,title,lines,def);
step = str2double(answer{1});
maxtime = str2double(answer{2});
timestep = str2double(answer{3});
dx = str2double(answer{4});
dy = str2double(answer{5});
r = str2double(answer{6});

% Calculate rate changes on grid
[RCREL,xx,yy] = calc_rcgrid(newt2,dx,dy,r,step,maxtime,timestep);

[existFlag,figNumber]=figure_exists('Rate change time slice',1);
newRCMapFlag=~existFlag;

% Set up figure for time slice plots
if newRCMapFlag
    rctfig = figure_w_normalized_uicontrolunits('tag','rtslice',...
        'Name','Rate change time slice',...
        'NumberTitle','off', ...
        'NextPlot','replace', ...
        'backingstore','on',...
        'Visible','on');
    rctax = axes('tag','axrctfig','NextPlot','replace','box','on');

    matdraw
    %
    uimenu('Label','Choose time slice', 'Callback','plot_timeslice')
end

% Other colormap
% cc = colormap;
% for i = 25:40
%     cc(i,:) = [1 1 1];
% end
% colormap(cc)

% get longitude / latitude
lon = a.Longitude; lat = a.Latitude;
% define grid
xmax = round(10*max(lon))/10+dx;
xmin = round(10*min(lon))/10-dx;
ymax = round(10*max(lat))/10+dy;
ymin = round(10*min(lat))/10-dy;

figure_w_normalized_uicontrolunits(rctfig)
%gcf=findobj('tag','rtslice')
%set(gcf,'Name','Rate change time slice');
set(gca,'tag','axrctfig');
hold on
xx = xmin-dx/2:dx:xmax-dx/2;
yy = ymax+dy/2:-dy:ymin+dy/2; yy = yy';
pcolor(xx,yy,RCREL(:,:,1))
shading flat
axis equal
caxis([-4 4])
colorbar

% plot_tslice=uicontrol('Style', 'pushbutton', 'String', 'Choose time slice','Units','normalized',...
%     'Position',[0.02 0.02 0.3 0.06],'Callback','plot_timeslice');
%hold off;

author: CelsoReyes


text:
function r8mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% R8MAT_PRINT prints an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end

author: johannesgerer


text:
function out=cumprod(x,dim)

if ~isempty(x)

 precision=x(1).precision;
 % can only handle up to 2 dimensions
 s=size(x);
 out=x;
 if nargin==1
  if any(s==1)
   for ii=2:numel(x)
    out(ii)=out(ii-1)*x(ii);
   end
  else
   for j=1:s(2)
    for i=2:s(1)
     out(i,j)=out(i-1,j)*x(i,j);
    end
   end
  end
 elseif nargin==2
  if dim==1
   for j=1:s(2)
    for i=2:s(1)
     out(i,j)=out(i-1,j)*x(i,j);
    end
   end
  elseif dim==2
   for j=2:s(2)
    for i=1:s(1)
     out(i,j)=out(i,j-1)*x(i,j);
    end
   end
  else
   error('mp cumprod can only handle up to 2-D arrays (11/04)');
  end
 end

else
 out=x;
end

author: opencobra

repo: cobratoolbox


text:
function varargout = lvd_EditContinuityConstraintGUI(varargin)
% LVD_EDITCONTINUITYCONSTRAINTGUI MATLAB code for lvd_EditContinuityConstraintGUI.fig
%      LVD_EDITCONTINUITYCONSTRAINTGUI, by itself, creates a new LVD_EDITCONTINUITYCONSTRAINTGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITCONTINUITYCONSTRAINTGUI returns the handle to a new LVD_EDITCONTINUITYCONSTRAINTGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITCONTINUITYCONSTRAINTGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITCONTINUITYCONSTRAINTGUI.M with the given input arguments.
%
%      LVD_EDITCONTINUITYCONSTRAINTGUI('Property','Value',...) creates a new LVD_EDITCONTINUITYCONSTRAINTGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditContinuityConstraintGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditContinuityConstraintGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditContinuityConstraintGUI

% Last Modified by GUIDE v2.5 19-Jan-2020 15:12:54

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditContinuityConstraintGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditContinuityConstraintGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditContinuityConstraintGUI is made visible.
function lvd_EditContinuityConstraintGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditContinuityConstraintGUI (see VARARGIN)

    % Choose default command line output for lvd_EditContinuityConstraintGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    constraint = varargin{1};
    setappdata(hObject, 'constraint', constraint);

    lvdData = varargin{2};
    setappdata(hObject,'lvdData',lvdData);

	populateGUI(handles, constraint, lvdData);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditContinuityConstraintGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditContinuityConstraintGUI);

function populateGUI(handles, constraint, lvdData)
    handles.constraintTypeLabel.String = constraint.getConstraintType();
    handles.constraintTypeLabel.TooltipString = constraint.getConstraintType();

    handles.scaleFactorText.String = fullAccNum2Str(constraint.getScaleFactor());

    handles.eventCombo.String = lvdData.script.getListboxStr();
    evtNum = lvdData.script.getNumOfEvent(constraint.event);
    if(not(isempty(evtNum)))
        handles.eventCombo.Value = evtNum;
    else
        handles.eventCombo.Value = 1;
    end

    evtNum = lvdData.script.getNumOfEvent(constraint.constraintEvent);
    if(not(isempty(evtNum)))
        handles.constraintEventCombo.Value = evtNum;
    else
        handles.constraintEventCombo.Value = 1;
    end

    handles.constraintEventCombo.String = lvdData.script.getListboxStr();

    handles.constActiveCheckbox.Value = constraint.active;

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditContinuityConstraintGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        constraint = getappdata(hObject, 'constraint');
        lvdData = getappdata(hObject,'lvdData');

        constraint.setScaleFactor(str2double(handles.scaleFactorText.String));

        numEvents = lvdData.script.getTotalNumOfEvents();
        if(numEvents > 0)
            eventNum = handles.eventCombo.Value;
            constraint.event = lvdData.script.getEventForInd(eventNum);

            eventNum = handles.constraintEventCombo.Value;
            constraint.constraintEvent = lvdData.script.getEventForInd(eventNum);
        else
            constraint.event = LaunchVehicleEvent.empty(1,0);
            constraint.constraintEvent = LaunchVehicleEvent.empty(1,0);
        end

        constraint.active = logical(handles.constActiveCheckbox.Value);

        varargout{1} = true;
        close(handles.lvd_EditContinuityConstraintGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.lvd_EditContinuityConstraintGUI);
    else
        msgbox(errMsg,'Invalid Constraint Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    sf = str2double(get(handles.scaleFactorText,'String'));
    enteredStr = get(handles.scaleFactorText,'String');
    numberName = 'Scale Factor';
    lb = 1E-12;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(sf, numberName, lb, ub, isInt, errMsg, enteredStr);

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditContinuityConstraintGUI);

    % --- Executes during object creation, after setting all properties.
function constraintTypeLabel_CreateFcn(hObject, eventdata, handles)
% hObject    handle to constraintTypeLabel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --- Executes on selection change in eventCombo.
function eventCombo_Callback(hObject, eventdata, handles)
% hObject    handle to eventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns eventCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from eventCombo

% --- Executes during object creation, after setting all properties.
function eventCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function scaleFactorText_Callback(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of scaleFactorText as text
%        str2double(get(hObject,'String')) returns contents of scaleFactorText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function scaleFactorText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to scaleFactorText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditContinuityConstraintGUI or any of its controls.
function lvd_EditContinuityConstraintGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditContinuityConstraintGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditContinuityConstraintGUI);
    end

% --- Executes on button press in constActiveCheckbox.
function constActiveCheckbox_Callback(hObject, eventdata, handles)
% hObject    handle to constActiveCheckbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of constActiveCheckbox

% --- Executes on selection change in constraintEventCombo.
function constraintEventCombo_Callback(hObject, eventdata, handles)
% hObject    handle to constraintEventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns constraintEventCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from constraintEventCombo

% --- Executes during object creation, after setting all properties.
function constraintEventCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to constraintEventCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [Res] = rmsr_inv(P, VecI, Mat, Pth)

% rmsr_inv replaces all descriptors from Mat (except the ones in the given descriptor vector) in the given position Pth and returns
%      the vector with higher S_err, using multiple linear regression analysis.
%
%
%	   Input:
%             P             Property vector
%             VecI          Descriptor vector
%             Mat           Descriptor matrix
%             Pth           Path to follow
%
%     Returns:
%
%             Res           Vector containing in the first place the
%                           S_err, and afterwards the corresponding descriptor vector
%
%
%
% Andrew G. Mercader
% INIFTA, La Plata, Argentina
% Created: 30 Jan 2007

if (nargin < 4)
   error('rmsr_inv requires at least 4 input variables. Type ''help rmsr_inv''.');
end

[k, n_m] = size(Mat);

Num=[1:n_m];
Vec2=VecI;
Vec2(Pth)=[];
Num(Vec2)=[];

[k,n_n]=size(Num);

Smax=0.00000000000000000000000000000000001;   % A very big number compared to a normal S necessary just to start the program
for j=1:n_n;
VecI(Pth)=Num(j);
Ser=rms_inv(P,VecI,Mat);
    if (Ser>Smax)
         Smax=Ser;
        desc=Num(j);
    end
end
VecI(Pth)=desc;
Res=[Smax, VecI];
%End of rmsr_inv


text:
function [rxnList, metList, Stats] = reactionCompare(CmodelIn, TmodelIn, scoreIn, varargin)
% This is the function front end for `autoMatchReactions` and
% the `reactionCompareGUI`. Loads `CMODEL`, `TMODEL`, `Score` and `ScoreTotal` from
% globals if they are not provided. If `rxnList` is not provided then it is
% generated with `autoMatchReactions`. If you are continuing comparison,
% include the optional inputs.
% Called by `optimalScores`, `autoMatchReactions`, `reactionCompareGUI`, calls `mergeModelsBorg`, `driveModelBorgifier`.
%
% USAGE:
%
%    [rxnList, metList, Stats] = reactionCompare(CmodelIn, TmodelIn, scoreIn, [rxnList, metList, Stats])
%
% INPUTS:
%    CmodelIn:    Comparison model
%    TmodelIn:    Template model
%    scoreIn:     Score
%
% OPTIONAL INPUTS:
%    rxnList:     Array pairs reactions in `CMODEL` with matches from `TMODEL` or
%                 declares them as new.
%    metList:     Array pairs metabolites in `CMODEL` with matches from `TMODEL`,
%                 new metabolites are given their new met number in `TMODEL`.
%    Stats:       Stats array that contains weighting information from previous
%                 scoring work.
%
% OUTPUTS:
%    rxnList:     reactions' list
%    metList:     metabolites' list
%    Stats:       Structure of information regarding the merging and also stats
%                 on the combined model.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

global CMODEL TMODEL SCORE % Declare variables
% Make the inputs variables so they can be easily accessed by downstream
% scripts. None of these variables are edited by any of these scripts.
CMODEL = CmodelIn ;
TMODEL = TmodelIn ;
SCORE = scoreIn ;

% Need scoreTotal now to create rxnList if Stats is not provided.
if nargin <= 5
%     Stats = optimalScores(CMODEL,TMODEL,SCORE) ;
    Stats = optimalScores ;
else
    Stats = varargin{3} ;
end

% Was metList supplied?
if nargin >= 5
    metList = varargin{2} ;
else
    metList = zeros(length(CMODEL.mets),1) ;
end

% How 'bout rxnList?
if nargin > 3
    rxnList = varargin{1} ;
else
    rxnList = ones(length(CMODEL.rxns),1)*-1 ;
end

%% Manual reaction comparison.

% Create information structure to pass to reactionCompareGUI
InfoBall.rxnList = rxnList ;
InfoBall.metList = metList ;
InfoBall.CmodelName = CMODEL.description ;
InfoBall.Stats = Stats ;
InfoBall.S = CMODEL.S ;

% Launch GUI.
[rxnList, metList, Stats] = reactionCompareGUI(InfoBall) ;


text:
function [ vX, mX ] = SolveLsL1ComplexCd( mA, vB, lambdaFctr, numIterations )
% ----------------------------------------------------------------------------------------------- %
%[ vX, mX ] = SolveLsL1ComplexPgm( mA, vB, lambdaFctr, numIterations )
% Solves the 0.5 * || A x - b ||_2 + \lambda || x ||_1 problem using
% Coordinate Descent Method. The model allows A, b and x to be Complex.
% Input:
%   - mA                -   Model Matrix.
%                           The model matrix.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - vB                -   Input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double' (Complex).
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  Wikipedia Coordinate Descent Method - https://en.wikipedia.org/wiki/Coordinate_descent.
% Remarks:
%   1.  Coordienat Descent is basically Steepest Descnt in L1 Norm.
% Known Issues:
%   1.  A
% TODO:
%   1.  B
% Release Notes:
%   -   1.0.000     07/11/2016
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

vANorm = sum(mA .* conj(mA), 1);
numElements = size(mA, 2); %<! Size of solution

vX  = pinv(mA) * vB; %<! Dealing with "Fat Matrix"

mX = zeros([size(vX, 1), numIterations]);
mX(:, 1) = vX;

for ii = 2:numIterations

    for jj = 1:numElements
         vA = mA(:, jj);
         colNormSqr = vANorm(jj);
         vExcCoord = [1:jj - 1, jj + 1:numElements];

         vR = vB - (mA(:, vExcCoord) * vX(vExcCoord));
         vBeta = vA' * vR;

         vX(jj) = ProxL1(vBeta, lambdaFctr) / colNormSqr;
    end

    mX(:, ii) = vX;

end

end

function [ vX ] = ProxL1( vX, lambdaFactor )

% Soft Thresholding - Complex Domain -> Keep Phase, Soft Threshold the
% Modulus

% vXAbs   = abs(vX);
% vXPhase = angle(vX);
%
% vX = max(vXAbs - lambdaFactor, 0) .* exp(1i * vXPhase);

vXAbs = abs(vX);

vX              = (vX ./ vXAbs) .* max((vXAbs - lambdaFactor), 0);
vX(vXAbs == 0)  = 0;

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [goodTF, msg] = verifyCelBodyData(celBodyData)
    %verifyCelBodyData Summary of this function goes here
    %   Detailed explanation goes here
    goodTF = true;
    msg = {};

    numTopLvlBodies = 0;
    topLevelBodies = KSPTOT_BodyInfo.empty(1,0);
    bodyIds = [];
    bodyNames = {};
    topLevelBodyIdStrs = {};
    celBodyFields = fields(celBodyData);
    for(i=1:length(celBodyFields))
        bInfo = celBodyData.(celBodyFields{i});
        bodyIds(end+1) = bInfo.id; %#ok<AGROW>
        bodyNames{end+1} = bInfo.name; %#ok<AGROW>

        if(isempty(bInfo.parentBodyInfo))
            numTopLvlBodies = numTopLvlBodies + 1;
            topLevelBodies(end+1) = bInfo; %#ok<AGROW>
            topLevelBodyIdStrs{end+1} = fullAccNum2Str(bInfo.id); %#ok<AGROW>
        end
    end

    if(numTopLvlBodies ~= 1)
        goodTF = false;

        idStr = strjoin(topLevelBodyIdStrs,', ');
        msg{end+1} = sprintf('Multiple bodies have been detected without parents.  Only one body should be the "top-level" body without a parent body.  Detected body IDs: %s', idStr);
    end

    if(numel(bodyIds) ~= numel(unique(bodyIds)))
        goodTF = false;

        msg{end+1} = 'Duplicate body IDs were found.  Check to ensure that all bodies have unique integer ID numbers.';
    end

    if(numel(bodyNames) ~= numel(unique(bodyNames)))
        goodTF = false;

        msg{end+1} = 'Duplicate body names were found.  Check to ensure that all bodies have unique name.';
    end
end