text:
function [pos] = surface_shift(pos, tri, amount)

% SURFACE_SHIFT inflates or deflates a triangulated surface by moving the
% vertices outward or inward along their normals.
%
% Use as
%   pos = surface_inflate(pos, tri, amount)
% where pos and tri describe the surface.
%
% See also SURFACE_NORMALS, SURFACE_ORIENTATION, SURFACE_INSIDE,
% SURFACE_NESTING

if isempty(tri)
  propos = elproj(pos);   % projection to 2D
  tri = delaunay(propos); % creates delaunay triangulation of 2D plane, which will be used for the the 3D case
end

nrm = surface_normals(pos, tri); % compute normals of surface
switch surface_orientation(pos, tri, nrm)
  case 'outward'
    % move along the normals
    nrm = +nrm;
  case 'inward'
    % move opposite to the normals
    nrm = -nrm;
  otherwise
    ft_warning('cannot determine the orientation of the surface');
end

% move the vertices with the specified amount along the normals
pos = pos + amount*nrm;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [rr,ar]=v_lpcrf2rr(rf,p);
%V_LPCRR2AR convert reflection coefs to autocorrelation coefs [RR,AR]=(RF,P)
%
% Inputs:  rf(:,n+1)  reflection coefficients: one row per frame
%          p          specifies number of rr coefficients to calculate (default=n)
% Outputs: rr(:,p+1)  autocorrelation coefficients
%          ar(:,n+1)  AR filter coefficients

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_lpcrf2rr.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(rf);
p0=p1-1;
if p0
   a = rf(:,2);
   rr=[ones(nf,1) -a zeros(nf,p0-1)];
   e = (a.^2-1);
   for n = 2:p0
      k=rf(:,n+1);
      rr(:,n+1) =k.*e - sum(rr(:,n:-1:2).*a,2);
      a = [a+k(:,ones(1,n-1)).*a(:,n-1:-1:1) k];
      e = e.*(1-k.^2);
   end
   ar = [ones(nf,1) a];
   r0=sum(rr.*ar,2).^(-1);
   rr=rr.*r0(:,ones(1,p1));
   if nargin>1 && ~isempty(p)
      if p<p0
         rr(:,p+2:p1)=[];
      else
         rr=[rr zeros(nf,p-p0)];
         af=-ar(:,p1:-1:2);
         for i=p0+1:p
            rr(:,i+1)=sum(af.*rr(:,i-p0+1:i),2);
         end
      end
   end
else
   rr=ones(nf,1);
   ar=rr;
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function demo_snirf()

% Syntax:
%    demo_snirf()
%
% Description:
%    Convert sample .nirs files included in the Homer3 project to SNIRF files (*.snirf) and
%    compare the saved and loaded SnirfClass objects for each coberted file to make sure
%    saving and loading of the .snirf is working correctly.
%

% Find the Homer3 Examples folder and cd to it
examplesDir = findexamplesdir();

% Delete any previously generated .snirf files to make sure to start from scratch
DeleteDataFiles(examplesDir, '.snirf', 'standalone');

% Start with .nirs files
nirsfiles = mydir([examplesDir, '*.nirs']);

fprintf('\n')

% Go through all the sample .nirs files
for ii=1:length(nirsfiles)

    % Convert sample .nirs file to SNIRF file (*.snirf)
    [pname,fname,ext] = fileparts([examplesDir, nirsfiles(ii).name]);
    pname = filesepStandard(pname);
    fprintf('Converting %s to %s\n', [pname, fname, ext], [pname, fname, '.snirf']);
    [snirf_saved, snirf_loaded] = snirf_load_save([pname, nirsfiles(ii).name]);

    % Compare the saved and loaded SnirfClass objects, using overloaded == operator
    if snirf_saved == snirf_loaded
        fprintf('Saved and loaded SnirfClass objects for %s are equal\n', [fname,'.snirf']);
    else
        fprintf('ERROR: Saved and loaded SnirfClass objects for %s are NOT equal\n', [fname,'.snirf']);
    end
    fprintf('\n');

end

author: BUNPC


text:
function Z = inormal(varargin)
% Applies a rank-based inverse normal transformation.
%
% Usage: Z = inormal(X)
%            inormal(X,c)
%            inormal(X,method)
%            inormal(X,...,quanti)
%
% Inputs:
% X      : Original data. Can be a vector or an array.
% c      : Constant to be used in the transformation.
%          Default c=3/8 (Blom).
% method : Method to choose c. Accepted values are:
%              'Blom'   (c=3/8),
%              'Tukey'  (c=1/3),
%              'Bliss', (c=1/2)  and
%              'Waerden' or 'SOLAR' (c=0).
% quanti : All data guaranteed to be quantitative and
%          without NaN?
%          This can be a true/false. If true, the function
%          runs much faster if X is an array.
%          Default is false.
%
% Outputs:
% Z      : Transformed data.
%
% References:
% * Van der Waerden BL. Order tests for the two-sample
%   problem and their power. Proc Koninklijke Nederlandse
%   Akademie van Wetenschappen. Ser A. 1952; 55:453-458
% * Blom G. Statistical estimates and transformed
%   beta-variables. Wiley, New York, 1958.
% * Tukey JW. The future of data analysis.
%   Ann Math Stat. 1962; 33:1-67.
% * Bliss CI. Statistics in biology. McGraw-Hill,
%   New York, 1967.
%
% _____________________________________
% Anderson M. Winkler
% Yale University / Institute of Living
% Aug/2011 (first version)
% Jun/2014 (this version)
% http://brainder.org

% Accept inputs & defaults
c0 = 3/8;  % Default (Blom, 1958)
quanti = false;
if nargin == 1,
    c = c0;
elseif nargin > 1 && ischar(varargin{2}),
    switch lower(varargin{2}),
        case 'blom'
            c = 3/8;
        case 'tukey'
            c = 1/3;
        case 'bliss'
            c = 1/2;
        case 'waerden'
            c = 0;
        case 'solar'
            c = 0; % SOLAR is the same as Van der Waerden
        otherwise
            error('Method %s unknown. Use either ''Blom'', ''Tukey'', ''Bliss'', ''Waerden'' or ''SOLAR''.',varargin{2});
    end
elseif nargin > 1 && isscalar(varargin{2}),
    c = varargin{2};  % For a user-specified value for c
end
if nargin == 3,
    quanti = varargin{3};
    if isempty(varargin{2}),
        c = c0;
    end
end
X = varargin{1};

% If the trait is quantitative, avoid the loop
if quanti,
    % Get the rank for each value
    [~,iX] = sort(X);
    [~,ri] = sort(iX);

    % Do the actual transformation
    N = size(X,1);
    p = ((ri-c)/(N-2*c+1));
    Z = sqrt(2)*erfinv(2*p-1);

else
    Z = nan(size(X));
    for x = 1:size(X,2),

        % Remove NaNs
        XX = X(:,x);
        ynan = ~isnan(XX);
        XX = XX(ynan);

        % Get the rank for each value
        [~,iX] = sort(XX);
        [~,ri] = sort(iX);

        % Do the actual transformation
        N = size(XX,1);
        p = ((ri-c)/(N-2*c+1));
        Y = sqrt(2)*erfinv(2*p-1);

        % Check for repeated values
        [U,~,IC] = unique(XX);
        if numel(U) < N,
            sIC = sort(IC);
            dIC = diff(vertcat(sIC,1));
            U = unique(sIC(~dIC));
            for u = 1:numel(U),
                Y(IC == U(u)) = mean(Y(IC == U(u)));
            end
        end

        % Put the NaNs back
        Z(ynan,x) = Y;
    end
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function [ a1, a2 ] = r8vec2_sort_a ( n, a1, a2 )

%*****************************************************************************80
%
%% R8VEC2_SORT_A ascending sorts an R8VEC2.
%
%  Discussion:
%
%    An R8VEC2 is two R8VEC's.
%
%    An R8VEC is a vector of R8 values.
%
%    Each item to be sorted is a pair (I,J), with the I
%    and J values stored in separate vectors A1 and A2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of items of data.
%
%    Input, real A1(N), A2(N), the data to be sorted.
%
%    Output, real A1(N), A2(N), the sorted data.
%
  if ( n <= 1 )
    return
  end
%
%  Initialize.
%
  i = 0;
  indx = 0;
  isgn = 0;
  j = 0;
%
%  Call the external heap sorter.
%
  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( n, indx, isgn );
%
%  Interchange the I and J objects.
%
    if ( 0 < indx )

      [ a1(i), a1(j) ] = r8_swap ( a1(i), a1(j) );
      [ a2(i), a2(j) ] = r8_swap ( a2(i), a2(j) );
%
%  Compare the I and J objects.
%
    elseif ( indx < 0 )

      isgn = r8vec2_compare ( n, a1, a2, i, j );

    elseif ( indx == 0 )

      break

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function t = subsasgn(t,s,b)
%SUBSASGN Subscripted assignement for ktensor.
%
%   See also KTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

switch s(1).type
    case '.'
        switch s(1).subs
            case 'lambda'
                if length(s) == 1
                    t = ktensor(b, t.u);
                else
                    newlambda = subsasgn(t.lambda, s(2:end), b);
                    t = ktensor(newlambda, t.u);
                end
            case {'u','U'}
                if length(s) == 1
                    t = ktensor(t.lambda, b);
                else
                    tmpu = subsasgn(t.u, s(2:end), b);
                    t = ktensor(t.lambda, tmpu);
                end
            otherwise
                error(['No such field: ', s(1).subs]);
        end
    case '()'
        error('Cannot change individual entries in a ktensor.')
    case '{}'
        new_s(1).type = '.';
        new_s(1).subs = 'u';
        new_s(2:length(s)+1) = s;
        t = subsasgn(t, new_s, b);
    otherwise
        error('Invalid subsasgn.');
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text: any(cellfun(@exist,{'lsl_pull_sample_buf', 'lsl_loadlib_', 'lsl_push_sample'})==3)

author: PatternRecognition

repo: OpenBMI


text:
% Builds and solves a simple least-squares problem using cvx

echo on

n = 100;
A = randn(2*n,n);
b = randn(2*n,1);

cvx_begin
   variable x1(n)
   minimize( norm( A*x1-b ) )
cvx_end
v1 = cvx_optval;

cvx_begin
   variable x2(n)
   minimize( sum_square( A*x2-b ) )
cvx_end
v2 = cvx_optval;

% The differences between the unsquared and squared versions should be small
norm( x1 - x2 ) / norm( x1 ) %#ok
abs( v1 - sqrt(v2) ) / v1 %#ok

echo off

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function ltfatmex(varargin)
%LTFATMEX   Compile Mex/Oct interfaces
%   Usage:  ltfatmex;
%           ltfatmex(...);
%
%   `ltfatmex` compiles the C backend in order to speed up the execution of
%   the toolbox. The C backend is linked to Matlab and Octave through Mex
%   and Octave C++ interfaces.
%   Please see INSTALL-Matlab or INSTALL-Octave for the requirements.
%
%   The action of `ltfatmex` is determined by one of the following flags:
%
%     'compile'  Compile stuff. This is the default.
%
%     'clean'    Removes the compiled functions.
%
%     'test'     Run some small tests that verify that the compiled
%                functions work.
%
%   The target to work on is determined by on of the following flags.
%
%   General LTFAT:
%
%     'lib'      Perform action on the LTFAT C library.
%
%     'mex'      Perform action on the mex / oct interfaces.
%
%     'pbc'      Perform action on the PolyBoolClipper code for use with MULACLAB
%
%     'auto'     Choose automatically which targets to work on from the
%                previous ones based on the operation system etc. This is
%                the default.
%
%   Block-processing framework related:
%
%     'playrec'  Perform action on the playrec code for use with real-time
%                block streaming framework.
%
%     'java'     Perform compilation of JAVA classes into the bytecode.
%                The classes makes the GUI for the blockproc. framework.
%
%   Other:
%
%      'verbose' Print action details.
%
%      'debug'   Build a debug version. This will disable compiler
%                optimizations and include debug symbols.
%

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: NA
%   REFERENCE: NA

% Verify that comp_pgauss is in path
if ~exist('comp_pgauss','file')
  disp(' ');
  disp('--- LTFAT - The Linear Time Frequency Analysis toolbox. ---');
  disp(' ')
  disp('To start the toolbox, call LTFATSTART as the first command.');
  disp(' ');
  return;
end;

bp=mfilename('fullpath');
bp=bp(1:end-length(mfilename));

definput.flags.target={'auto','lib','mex','pbc','playrec','java','blockproc'};
% This has to be also defined
definput.flags.comptarget={'release','debug'};
definput.flags.command={'compile','clean','test'};
definput.flags.verbosity={'noverbose','verbose'};
definput.flags.mexapi={'auto','legacy'};
definput.keyvals.mexfile = [];
definput.keyvals.jobs = [];

s = which('ltfatarghelper');
if strcmpi(mexext,s(end-numel(mexext)+1:end))
   % We must avoid calling corrupt ltfatarghelper.mexext
   % The following must not fail
   try
       ltfatarghelper({},struct(),{});
   catch
       if any(strcmpi('verbose',varargin))
           fprintf('Removing corrupt ltfatarghelper.%s',mexext);
       end
       delete(s);
       % Now there is just a ltfatarghelper.m
   end
end

[flags,kv]=ltfatarghelper({},definput,varargin);

% Remember the current directory.
curdir=pwd;

% Compile backend lib?
do_lib  = flags.do_lib || flags.do_auto;
% Compile MEX/OCT interfaces?
do_mex  = flags.do_mex || flags.do_auto;
% Compile MEX PolyBoolClipper.mex
do_pbc  = flags.do_pbc || flags.do_auto;
% Compile MEX playrec.mex... using Portaudio library.
% (relevant only for the bloc processing framework)
do_playrec  = flags.do_playrec || flags.do_blockproc;
% Compile Java classes containing GUI for the bloc proc. framework.
do_java  = flags.do_java || flags.do_blockproc;

if isoctave
	extname='oct';
    ext='oct';
else
    extname='mex';
    ext=mexext;
end;

fftw_lib_names = {'fftw3', 'fftw3f' };

% Check if we are on Windows
if ispc
    makefilename='Makefile_mingw';
    make_exe = 'mingw32-make';
    sharedExt = 'dll';
    %fftw_lib_names = {'fftw3', 'fftw3f' };
    % The pre-compiled Octave for Windows comes only in 32bit version (3.6.4)
    % We use different Makefiles
    if isoctave
      makefilename='Makefile_mingwoct';
    end
end;

% Check if we are on Unix-type system
if isunix
   makefilename='Makefile_unix';
   make_exe = 'make';
   sharedExt = 'so';
end;

% Check if we are on Mac
if ismac
   makefilename='Makefile_mac';
   make_exe = 'make';
   sharedExt = 'dylib';
end;

clear mex;
% -------------- Handle cleaning --------------------------------
if flags.do_clean

  if do_lib
    disp('========= Cleaning libltfat ===============');
    cd([bp,'lib',filesep,'ltfatcompat']);
    callmake(make_exe,makefilename,'target','clean',flags.verbosity);
    %[status,result]=system([make_exe, ' -f ',makefilename,' clean']);
    %disp('Done.');
  end;

  if do_mex
    fprintf('========= Cleaning %s interfaces =========\n', upper(extname));
     cd([bp,extname]);
     callmake(make_exe,makefilename,'target','clean','ext',ext,flags.verbosity);
     %[status,result]=system([make_exe, ' -f ',makefilename,' clean',...
     %                ' EXT=',ext]);
  end;

  if do_pbc
     % Use the correct makefile
     makefilename_pbc = makefilename;
     if isoctave
       if ~strcmpi(makefilename(end-2:end),ext)
          makefilename_pbc = [makefilename,ext];
       end
    end

    disp('========= Cleaning PolyBoolClipper ====================');
    cd([bp,'thirdparty',filesep,'polyboolclipper']);
    clear mex;
    callmake(make_exe,makefilename_pbc,'target','clean','ext',mexext,flags.verbosity);
    %[status,result]=system([make_exe, ' -f ',makefilename,' clean',' EXT=',mexext]);
  end;

  if do_playrec
     % Use the correct makefile
     if isoctave
       if ~strcmpi(makefilename(end-2:end),ext)
          makefilename = [makefilename,ext];
       end
    end

    disp('========= Cleaning PLAYREC ================');
    cd([bp,'thirdparty',filesep,'Playrec']);
    clear mex;
    %[status,result]=system([make_exe, ' -f ',makefilename,' clean',' EXT=',mexext]);
    callmake(make_exe,makefilename,'target','clean','ext',mexext,flags.verbosity);
  end;

  if do_java
    disp('========= Cleaning JAVA ================');
    cd([bp,'blockproc',filesep,'java']);
    %[status,result]=system([make_exe,' clean']);
    callmake(make_exe,[],'target','clean',flags.verbosity);
  end;

  cd(curdir);
end;

% -------------- Handle compiling  --------------------------------

if flags.do_compile
  if do_lib
    disp('========= Compiling libltfat ==============');
    cd([bp,'lib',filesep,'ltfatcompat']);
    clear mex;

    dfftw = ['-l',fftw_lib_names{1}];
    sfftw = ['-l',fftw_lib_names{2}];
    if ispc && ~isoctave
        fftw_lib_found_names = searchfor(bp,fftw_lib_names,sharedExt);
        if ~isempty(fftw_lib_found_names)
           dfftw = ['-l:',fftw_lib_found_names{1}];
           sfftw = ['-l:',fftw_lib_found_names{2}];
       end
    end
      % DFFTW and SFFTW are not used in the unix_makefile
      [status,result] = callmake(make_exe,makefilename,'matlabroot','arch',...
                       'dfftw',dfftw,'sfftw',sfftw,flags.verbosity,...
                       'comptarget',flags.comptarget,'jobs',kv.jobs);
      if(~status)
        disp('Done.');
      else
        error('Failed to build LTFAT libs:\n %s',result);
      end
    %end;
  end;

  if do_mex
    fprintf('========= Compiling %s interfaces ========\n', upper(extname));
    clear mex;
    cd([bp,extname]);

    dfftw = ['-l',fftw_lib_names{1}];
    sfftw = ['-l',fftw_lib_names{2}];
    if ~isoctave
        fftw_lib_found_names = searchfor(bp,fftw_lib_names,sharedExt);
        if ~isempty(fftw_lib_found_names)
            if ~ismac
               dfftw = ['-l:',fftw_lib_found_names{1}];
               sfftw = ['-l:',fftw_lib_found_names{2}];
            else
                % We need a full path here.
               dfftw = [binDirPath(),filesep,fftw_lib_found_names{1}];
               sfftw = [binDirPath(),filesep,fftw_lib_found_names{2}];
            end
       end
    end

    [status,result] = callmake(make_exe,makefilename,'matlabroot','arch',...
                      'ext',ext,'dfftw',dfftw,'sfftw',sfftw,...
                      flags.verbosity,'target',kv.mexfile,...
                      'comptarget',flags.comptarget,'jobs',kv.jobs,...
                      flags.mexapi);

    if(~status)
      disp('Done.');
    else
      error('Failed to build %s interfaces: %s \n',upper(extname),result);
    end
  end;

  if do_pbc
    makefilename_pbc = makefilename;
    if isoctave
       if ~strcmpi(makefilename(end-2:end),ext)
          makefilename_pbc = [makefilename,ext];
       end
    end

    disp('========= Compiling PolyBoolClipper ===================');
    % Compile PolyBoolClipper mex file for use with mulaclab
    cd([bp,'thirdparty',filesep,'polyboolclipper']);
    clear mex;
    [status,result] = callmake(make_exe,makefilename_pbc,'matlabroot','arch',...
                      'ext',ext,flags.verbosity);

    if(~status)
      disp('Done.');
    else
      error('Failed to build PlyBoolClipper:\n %s',result);
    end
  end;
  if do_playrec
    disp('========= Compiling PLAYREC ===============');
    cd([bp,'thirdparty',filesep,'Playrec']);
    clear mex;
    % Compile the Playrec (interface to portaudio) for the real-time block-
    % stream processing

     portaudioLib = '-lportaudio';

     binArchPath = binDirPath();
       playrecRelPath = ['thirdparty',filesep,'Playrec'];

       foundPAuser = [];
       if ispc
          foundPAuser = dir([bp,playrecRelPath,filesep,'*portaudio*',sharedExt,'*']);
       end

       foundPAmatlab = [];
       if ~isoctave
          % Check if portaudio library is present in the Matlab installation
          foundPAmatlab = dir([binArchPath,filesep,'*portaudio*',sharedExt,'*']);
       end

       if ~isempty(foundPAuser)
          if numel(foundPAuser)>1
             error('Ambiguous portaudio libraries in %s. Please leave just one.',playrecRelPath);
          end
          foundPAuser = foundPAuser(1).name;

       elseif ~isempty(foundPAmatlab)
          if numel(foundPAmatlab)>1
             if ispc
                %This should not happen on Windows
                %Use the first one on Linux
                error('Ambiguous portaudio libraries in %s.',binArchPath);
             end
          end
             foundPAmatlab = foundPAmatlab(1).name;
       else
          if ispc && isoctave || ispc
          error(['Portaudio not found. Please download Portaudio http://www.portaudio.com\n',...
                 'and build it as a shared library and copy it to the\n',...
                 '%s directory. \n'],playrecPath);
          end
       end

    if isoctave
       if ~strcmpi(makefilename(end-2:end),ext)
          makefilename = [makefilename,ext];
       end
    end

    doPAuser = ~isempty(foundPAuser);
    doPAmatlab = ~isempty(foundPAmatlab) && ~doPAuser;

    if doPAmatlab
       if ismac
          % Full path is needed on MAC since
          % clang does not understand -l: prefix.
          portaudioLib = [binArchPath,filesep,foundPAmatlab];
       else
          portaudioLib = ['-l:',foundPAmatlab];
       end
       fprintf('    ...using %s from Matlab installation.\n',foundPAmatlab);
    elseif doPAuser
        portaudioLib = ['-l:',foundPAuser];
        fprintf('   ...using %s from ltfat%s%s.\n',...
                  foundPAuser,filesep,playrecRelPath);
    end

    [status,result] = callmake(make_exe,makefilename,'matlabroot','arch',...
                      'ext',mexext,'portaudio',portaudioLib,'extra','HAVE_PORTAUDIO',...
                      flags.verbosity);
    if(~status)
      disp('Done.');
    else
      error('Failed to build PLAYREC:\n %s',result);
    end
  end;

  if do_java
    disp('========= Compiling JAVA classes ===================');
    % Compile the JAVA classes
    cd([bp,'blockproc',filesep,'java']);
    clear mex;
    [status,result] = callmake(make_exe,'Makefile',flags.verbosity);
    if(~status)
      disp('Done.');
    else
      error('Failed to build JAVA classes:\n %s',result);
    end
  end;
end;

% -------------- Handle testing ---------------------------------------

if flags.do_test

  if do_mex

    fprintf('========= Testing %s interfaces ==========\n', extname);
    fprintf('1.: Test if comp_pgauss.%s was compiled: ',ext);
    fname=['comp_pgauss.',ext];
    if exist(fname,'file')
      disp('SUCCESS.');
    else
      disp('FAILED.');
    end;

    fprintf('2.: Test if pgauss executes:              ');
    pgauss(100);
    % If the execution of the script makes it here, we know that pgauss
    % did not crash the system, so we can just print success. Same story
    % with the following entries.
    disp('SUCCESS.');

    fprintf('3.: Test if fftreal executes:             ');
    fftreal(randn(10,1),10);
    disp('SUCCESS.');

    fprintf('4.: Test if dgt executes:                 ');
    dgt(randn(12,1),randn(12,1),3,4);
    disp('SUCCESS.');

  end;

end;

% Jump back to the original directory.
cd(curdir);

function status = filesExist(filenames)
   if(~iscell(filenames))
      filenames={filenames};
   end
   for ii=1:length(filenames)
      filename = filenames{ii};
      if(~exist(filename,'file'))
         error('%s: File %s not found.',mfilename,filename);
      end
   end

function found_files=searchfor(bp,files,sharedExt)

found_names = {};
      if ispc
         for ii=1:numel(files)
            % Search the ltfat/mex lib
            L = dir([bp,'mex',filesep,'*',files{ii},'*.',sharedExt]);
            if isempty(L)
                error(['%s: %s could not be found in ltfat/mex subdir.',...
                       ' Please download the FFTW dlls and install them.'],...
                      upper(mfilename),files{ii});
            end
            found_files{ii} = L(1).name;
            fprintf('   ...using %s from ltfat/mex.\n',L(1).name);
         end
      elseif isunix
          for ii=1:numel(files)
             L = dir([binDirPath(),filesep,'*',files{ii},'*.',sharedExt,'*']);

             if isempty(L)
                 error('%s: Matlab FFTW libs were not found. Strange.',...
                      upper(mfilename));
             end

             found_files{ii} = L(1).name;

             fprintf('   ...using %s from Matlab installation.\n',...
                     found_files{ii});
          end

      end;

function path=binDirPath()
path = [matlabroot,filesep,'bin',filesep,computer('arch')];

function [status,result]=callmake(make_exe,makefilename,varargin)
%CALLMAKE
%   Usage:  callmake(make_exe,makefilename);
%           callmake(make_exe,makefilename,'matlabroot',matlabroot,...);
%
%   `callmake(make_exe,makefilename)` is a platform independent wrapper for
%   calling the make command `make_exe` on `makefilename` file. When
%   `makefilename` is missing or is empty, the default `Makefile` file is
%   used.
%
%   `callmake(...,'target',target)` used `target` from the makefile.
%
%   Flags:
%
%       matlabroot:   Pass MATLABROOT=matlabroot variable to the makefile.
%
%       arch:         Pass ARCH=computer('arch') variable to the makefile.
%
%   Key-value parameters:
%
%       ext:          Pass EXT variable to the makefile.
%
%       portaudio:    Pass PORTAUDIO variable to the makefile.
%
%       dfftw:        Pass DFFTW variable to the makefile.
%
%       sfftw:        Pass SFFTW variable to the makefile.

  if nargin < 2 || isempty(makefilename)
     systemCommand = make_exe;
  else
     systemCommand = [make_exe, ' -f ',makefilename];
  end
  definput.flags.matlabroot={'none','matlabroot'};
  definput.flags.arch={'none','arch'};
  definput.keyvals.ext=[];
  definput.keyvals.dfftw=[];
  definput.keyvals.sfftw=[];
  definput.keyvals.target=[];
  definput.keyvals.comptarget=[];
  definput.keyvals.portaudio=[];
  definput.keyvals.extra=[];
  definput.keyvals.jobs =[];
  definput.flags.targetapi={'auto','legacy'};
  definput.flags.verbosity={'noverbose','verbose'};
  [flags,kv]=ltfatarghelper({},definput,varargin);

  if flags.do_matlabroot
     systemCommand = [systemCommand, ' MATLABROOT=','"',matlabroot,'"'];
  end

  if flags.do_arch
     systemCommand = [systemCommand, ' ARCH=',computer('arch')];
  end

  if ~isempty(kv.ext)
     systemCommand = [systemCommand, ' EXT=',kv.ext];
  end

  if ~isempty(kv.dfftw)
     systemCommand = [systemCommand, ' DFFTW=',kv.dfftw];
  end

  if ~isempty(kv.sfftw)
     systemCommand = [systemCommand, ' SFFTW=',kv.sfftw];
  end

  if ~isempty(kv.portaudio)
     systemCommand = [systemCommand, ' PORTAUDIO=',kv.portaudio];
  end

  if ~isempty(kv.comptarget) && ~strcmpi(kv.comptarget,'release')
     systemCommand = [systemCommand, ' COMPTARGET=',kv.comptarget];
  end

  if ~isempty(kv.extra)
     systemCommand = [systemCommand, ' ',kv.extra,'=1'];
  end

  if ~isempty(kv.jobs)
      if isnumeric(kv.jobs)
        systemCommand = [systemCommand,sprintf(' -j%d',kv.jobs)];
      elseif ischar(kv.jobs)
        systemCommand = [systemCommand,sprintf(' -j%s',kv.jobs)];
      end
  end

  if ~isoctave && ~flags.do_legacy
      matlabversion = version('-release');
      if str2double(matlabversion(1:4)) >= 2018
          systemCommand = [systemCommand,' POST2018a=1'];
      end
  end

  if flags.do_verbose
      fprintf('Calling:\n    %s\n\n',systemCommand);
  end

  if ~isempty(kv.target)
       if ~iscell(kv.target), kv.target = {kv.target}; end
       for ii = 1:numel(kv.target)
           [status,result]=system([systemCommand,' ',kv.target{ii}]);
       end
  else
        [status,result]=system(systemCommand);
  end

  if flags.do_verbose && ~isoctave
     disp(result);
  end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% ALS Completion
% as described in
%
%   Michael Steinlechner, Riemannian optimization for high-dimensional tensor completion,
%   Technical report, March 2015, revised December 2015.
%   To appear in SIAM J. Sci. Comput.
%

%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt
function [X,cost,test,stats] = completion_als( A_Omega, Omega, A_Gamma, Gamma, X, opts )

    if ~isfield( opts, 'maxiter');  opts.maxiter = 100;     end
    if ~isfield( opts, 'tol');      opts.tol = 1e-6;        end
    if ~isfield( opts, 'reltol');   opts.reltol = 1e-6;     end

	n = X.size;
	r = X.rank;
    d = X.order;

	cost = zeros(2*opts.maxiter,1);
	test = zeros(2*opts.maxiter,1);

    norm_A_Omega = norm( A_Omega );
    norm_A_Gamma = norm( A_Gamma );

    X = orthogonalize( X, 1 );

    t = tic;
    stats.time = [0];
    stats.conv = false;

	for i = 1:opts.maxiter

        % ===================
        % FORWARD SWEEP:
        % ===================
        fprintf(1,'Currently optimizing core: ')
        for mu = 1:d-1
            fprintf(1,'%i ', mu)
            X.U{mu} = solve_least_squares( A_Omega, Omega, X, mu );
            X = orth_at( X, mu, 'left' );
        end
		cost(2*i-1) = sqrt(2*func(A_Omega, X, Omega)) / norm_A_Omega;

        if cost(2*i-1) < opts.tol
            disp(sprintf('CONVERGED AFTER %i HALF-SWEEPS. Rel. residual smaller than %0.3g', ...
                          2*i-1, opts.tol))
            stats.conv = true;
            cost = cost(1:2*i-1,1);
            stats.time = [stats.time stats.time(end)+toc(t)];
            test(2*i-1) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
            test = test(1:2*i-1,1);
            break
        end

        if i > 1
            reltol = abs(cost(2*i-1) - cost(2*i-2)) / cost(2*i-1);
            if reltol < opts.reltol
                disp(sprintf('No more progress in gradient change, but not converged after %i half-sweeps. ABORTING!. \nRelative change is smaller than %0.3g', ...
                              i, opts.reltol))
                stats.conv = false;
                cost = cost(1:2*i-1,1);
                stats.time = [stats.time stats.time(end)+toc(t)];
                test(2*i-1) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
                test = test(1:2*i-1,1);
                break
            end
        end

        stats.time = [stats.time stats.time(end)+toc(t)];
        test(2*i-1) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
        t = tic;

        fprintf(1,'\nFinished forward sweep.\n    Cost: %e\n    Test: %e\n', cost(2*i-1), test(2*i-1) );
        % ===================
        % BACKWARD SWEEP:
        % ===================
        fprintf(1,'Currently optimizing core: ')
        for mu = d:-1:2
            fprintf(1,'%i ', mu)
            X.U{mu} = solve_least_squares( A_Omega, Omega, X, mu );
            X = orth_at( X, mu, 'right' );
        end

		cost(2*i) = sqrt(2*func(A_Omega, X, Omega)) / norm_A_Omega;

        if cost(2*i) < opts.tol
            disp(sprintf('CONVERGED AFTER %i HALF-SWEEPS. Rel. residual smaller than %0.3g', ...
                          2*i, opts.tol))
            stats.conv = true;
            cost = cost(1:2*i,1);
            stats.time = [stats.time stats.time(end)+toc(t)];
            test(2*i) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
            test = test(1:2*i,1);
            break
        end

        if i > 1
            reltol = abs(cost(2*i) - cost(2*i-1)) / cost(2*i);
            if reltol < opts.reltol
                disp(sprintf('No more progress in gradient change, but not converged after %i half-sweeps. ABORTING!. \nRelative change is smaller than %0.3g', ...
                              2*i, opts.reltol))
                stats.conv = false;
                cost = cost(1:2*i,1);
                stats.time = [stats.time stats.time(end)+toc(t)];
                test(2*i) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
                test = test(1:2*i,1);
                break
            end
        end

        stats.time = [stats.time stats.time(end)+toc(t)];
        test(2*i) = sqrt(2*func(A_Gamma, X, Gamma)) / norm_A_Gamma;
        t = tic;
        fprintf(1,'\nFinished backward sweep.\n    Cost: %e\n    Test: %e\n', cost(2*i), test(2*i) );

        disp('_______________________________________________________________')
    end

    % This is to match original shape of stats.time, since we artificially start w/ [0]
    % for consistency in how we count time
    stats.time = stats.time(2:end);

end

function res = func(A_Omega, X, Omega)
	res = 0.5*norm( A_Omega - X(Omega) )^2;
end

function res = solve_least_squares( A_Omega, Omega, X, mu )

    n = X.size;
    d = X.order;
    r = X.rank;

    [jmu,idx] = sort(Omega(:,mu),'ascend');
    Omega = Omega(idx,:);
    A_Omega = A_Omega(idx);

    C = cell(1,d);
    for i=1:d
        C{i} = permute( X.U{i}, [1 3 2]);
    end
    res = zeros( size(C{mu}) );

    %B = zeros(size(Omega,1), r(mu)*r(mu+1));

    %imu = 1;
    %for sample = 1:size(Omega,1)

    %    L = 1;
    %    for i = 1:mu-1
    %        L = L * C{i}(:,:,Omega(sample,i));
    %    end

    %    R = 1;
    %    for i = d:-1:mu+1
    %        R = C{i}(:,:,Omega(sample,i)) * R;
    %    end
    %
    %    %B(sample,:) = kron(R',L);
    %    B(sample,:) = reshape( L'*R', 1, r(mu)*r(mu+1) );
    %end

    B = als_solve_mex( n, r, C, Omega', mu)';

    for i = 1:X.size(mu)
        idx = find(jmu == i);

        if isempty(idx)
            error('No samples for this slice!')
        end
        res(:,:,i) = reshape(B(idx,:)\A_Omega(idx), r(mu), r(mu+1));
    end

    res = permute( res, [1 3 2] );

end

author: NicolasBoumal

repo: manopt


text:
function transportRxns=getTransportRxns(model)
% getTransportRxns
%   Retrieves the transport reactions from a model
%
%   model           a model structure
%
%   transportRxns   logical array with true if the corresponding
%                   reaction is a transport reaction
%
%   Transport reactions are defined as reactions involving (at least) one
%   metabolite name in more than one compartment.
%
%   Usage: transportRxns=getTransportRxns(model)

transportRxns=false(numel(model.rxns),1);

for i=1:numel(model.rxns)
    %Get the involved metabolites in each reaction
    mets=model.metNames(model.S(:,i)~=0);
    transportRxns(i)=numel(mets)~=numel(unique(mets));
end
end

author: SysBioChalmers

repo: RAVEN


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% wecSimMCR
% WEC-Sim multiple condition run executable
clear mcr imcr i j k l1 l2 m n name nseed kkk len numConditions
clear body waves simu output pto constraint ptoSim

evalc('wecSimInputFile');

if isempty(simu.mcrMatFile) == 0
    load(simu.mcrMatFile);
else
    kkk=0;
    mcr.header = {'waves.height','waves.period'};
    if isempty(simu.mcrExcelFile) == 0
        mcr.waveSS = xlsread(simu.mcrExcelFile);
        mcr.waveSS(isnan(mcr.waveSS))=0;
        for i=2:length(mcr.waveSS(:,1))
            for j=2:length(mcr.waveSS(1,:))
                if (mcr.waveSS(i,j)>0)
                    kkk = kkk+1;
                    mcr.cases(kkk,1) = mcr.waveSS(i,1);
                    mcr.cases(kkk,2) = mcr.waveSS(1,j);
                end
            end
        end
    else
        for i=1:length(waves.height)
            for j=1:length(waves.period)
                kkk = kkk+1;
                mcr.cases(kkk,1) = waves.height(i);
                mcr.cases(kkk,2) = waves.period(j);
            end
        end
    end

    numConditions=2;
    if length(waves.phaseSeed)>1
        numConditions=numConditions+1;
        mcr.header{numConditions} = 'waves.phaseSeed';
        len = length(mcr.cases(:,1));
        for nseed=1:length(waves.phaseSeed)
            mcr.cases(len*(nseed-1)+1:len*(nseed-1)+len,1:numConditions-1) = mcr.cases(1:len,1:numConditions-1);
            mcr.cases(len*(nseed-1)+1:len*(nseed-1)+len,    numConditions) = waves.phaseSeed(nseed);
        end
    end

    if exist('pto','var')
        for n=1:size(pto,2)
            if (length(pto(n).damping)>1 || length(pto(n).stiffness)>1)
                numConditions=numConditions+2;
                name = sprintf('pto(%i).damping', n);
                mcr.header{numConditions-1} = name;
                name = sprintf('pto(%i).stiffness', n);
                mcr.header{numConditions  } = name;

                len = length(mcr.cases(:,1)); kkk = 0;
                for l2=1:length(pto(n).stiffness)
                    for l1=1:length(pto(n).damping)
                        kkk=kkk+1;
                        mcr.cases(len*(kkk-1)+1:len*(kkk-1)+len,1:numConditions-2) = mcr.cases(1:len,1:numConditions-2);
                        mcr.cases(len*(kkk-1)+1:len*(kkk-1)+len,  numConditions-1) = pto(n).damping(l1);
                        mcr.cases(len*(kkk-1)+1:len*(kkk-1)+len,    numConditions) = pto(n).stiffness(l2);
                    end
                end
            end
        end; clear i j k l1 l2 m n name nseed kkk len numConditions
    end
end

%% Execute wecSimMCR
% Run WEC-Sim
warning('off','MATLAB:DELETE:FileNotFound'); delete('mcrCase*.mat')
for imcr=1:length(mcr.cases(:,1))
    wecSim;
    if exist('userDefinedFunctionsMCR.m','file') == 2
        userDefinedFunctionsMCR;
    end

    %% Store hydrodata in memory for reuse in future runs.
    if simu.reloadH5Data == 0 && imcr == 1        % Off->'0', On->'1', (default = 0)
        for ii = 1:simu.numHydroBodies
            hydroData(ii) = body(ii).hydroData;
        end
    end
end; clear imcr ans hydroData

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

positive = find(y == 1);
negative = find(y == 0);

plot(X(positive, 1), X(positive, 2), 'k+', 'LineWidth', 2, 'MarkerSize', 7);
plot(X(negative, 1), X(negative, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

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

hold off;

end

repo: MachineLearning


text:
classdef AxAyComputerFromVolumeAndR < handle

    properties (Access = private)
        theta
        r
        cx
        cy
    end

    methods (Access = public)

        function obj = AxAyComputerFromVolumeAndR(cParams)
            obj.init(cParams)
        end

        function [ax,ay] = compute(obj)
            h  = obj.computeH(obj.theta,obj.r);
            Tx = obj.computeTx(obj.r,h,obj.cx);
            Ty = obj.computeTy(obj.r,h,obj.cx);
            ax = obj.computeAx(Tx,Ty,obj.cx);
            ay = obj.computeAy(Tx,Ty,obj.cx);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.theta = 1 - cParams.volume;
            obj.r     = cParams.r;
            obj.cx    = cParams.cx;
            obj.cy    = cParams.cy;
        end

    end

    methods (Access = private, Static)

        function Tx = computeTx(r,h,cx)
            n = r*h*(1+h);
            d = (1+r*h)*cx;
            Tx = n/d;
        end

        function Ty = computeTy(r,h,cx)
            n = cx*h*(1+r*h);
            d = (1+h);
            Ty = n/d;
        end

        function ax = computeAx(Tx,Ty,cx)
            ax = (cx - Ty)/(1-Tx*Ty);
        end

        function ay = computeAy(Tx,Ty,cx)
            ay = (1-Tx*cx)/(1-Tx*Ty);
        end

        function h = computeH(theta,r)
            n = -(1+r)*theta + sqrt(theta^2*(r-1)^2 + 4*r);
            d = 2*r*(1+theta);
            h = n/d;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function M = obj_read(filename)
% Read Wavefront OBJ-formatted data from disk
% FORMAT M = obj_read(filename)
%
% filename - OBJ-formatted file name
% M        - data structure
%__________________________________________________________________________
%
% Wavefront OBJ Format Specification:
% https://en.wikipedia.org/wiki/Wavefront_.obj_file
%__________________________________________________________________________
% Copyright (C) 2017 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: obj_read.m 7390 2018-08-13 09:51:20Z guillaume $

fid = fopen(filename,'rt');
if fid == -1
    error('Cannot open %s.',filename);
end

M = struct('vertices',[],'faces',[]);

while true
    l = fgetl(fid);
    if ~ischar(l), break; end
    if numel(l) < 1 || isempty(strtrim(l)) || l(1) == '#', continue; end
    switch l(1)
        case 'v'
            switch l(2)
                case 't'
                    % texture coordinates, in (u, v [,w]) coordinates
                    t = sscanf(l(2:end),'%f %f %f');
                case 'n'
                    % vertex normals in (x,y,z) form
                    n = sscanf(l(2:end),'%f %f %f');
                case 'p'
                    % Parameter space vertices in (u [,v] [,w]) form
                    p = sscanf(l(2:end),'%f %f %f');
                otherwise
                    v = sscanf(l(2:end),'%f %f %f');
                    if numel(v) > 3, v = v(1:3); end
                    M.vertices(size(M.vertices,1)+1,:) = v;
            end
        case 'f'
            f = sscanf(l(2:end),'%d %d %d');
            if numel(f) ~= 3
                f = sscanf(l(2:end),'%d/%d'); % '%d/%d %d/d %d/%d'
                if numel(f) ~= 6
                    f = sscanf(l(2:end),'%d//%d'); % '%d//%d %d//d %d//%d'
                    if numel(f) ~= 6
                        f = sscanf(l(2:end),'%d/%d/%d'); % '%d/%d/%d %d/%d/%d %d/%d/%d'
                        if numel(f) == 9
                            f = f([1 4 7]);
                        else
                            fprintf('Not a triangle.\n');
                            continue;
                        end
                    else
                        f = f([1 3 5]);
                    end
                else
                    f = f([1 3 5]);
                end
            end
            i = find(f<0);
            if isempty(i), f(i) = size(M.vertices,1) + f(i); end
            M.faces(size(M.faces,1)+1,:) = f;
        case 'o'
            fprintf('Ignoring named objects.\n');
        case 'g'
            fprintf('Ignoring polygon groups.\n');
        case 's'
            fprintf('Ignoring smooth shading.\n');
        otherwise
            if strncmp('mtllib',l,6) || strncmp('usemtl',l,6)
                fprintf('Ignoring materials.\n');
            else
                fprintf('Ignoring line starting with %c.\n',l(1));
            end
    end
end

fclose(fid);

author: spm


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu=1./m.*sum(X)

%for i=1:n
%	sigma2(i)=1./m.*sum((X-mu(i)).^2);
%end
sigma2=var(X).*(m-1)./m;

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

end

author: yhyap

repo: machine-learning-coursera


text:
function [B0, beta0, phi0, S0, alpha0]=nwprior(ar,arvar,lambda1,lambda3,lambda4,n,m,p,k,q,prior,priorexo)

% function [B0 beta0 phi0 S0 alpha0]=nwprior(ar,arvar,lambda1,lambda3,lambda4,n,m,p,k,q,prior,bex)
% returns prior values from hyperparameters, for the normal-Wishart prior
% inputs:  - scalar 'ar': prior value of the autoregressive coefficient on own first lag (defined p 15 of technical guide)
%          - vector 'arvar': residual variance of individual AR models estimated for each endogenous variable
%          - scalar 'lambda1': overall tightness hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda2': cross-variable weighting hyperparameter(defined p 16 of technical guide)
%          - scalar 'lambda3': lag decay hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda4': exogenous variable tightness hyperparameter (defined p 17 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
%          - integer 'prior': value to determine which prior applies to the model
% outputs: - matrix 'B0': the non-vectorised form of beta0
%          - vector 'beta0': vector of prior values for beta (defined in 1.3.4)
%          - matrix 'phi0': prior covariance matrix for the VAR coefficients in the case of a normal-Wishart prior (defined in 1.4.7)
%          - matrix 'S0': prior scale matrix for sigma (defined in 1.4.11)
%          - integer 'alpha0': prior degrees of freedom for sigma (defined in 1.4.11)

% start with beta0, defined in (1.3.4)
beta0=zeros(q,1);

idx = 1:n;
if isscalar(ar)
    beta0((idx-1)*k+idx,1) = ar;
else
    beta0((idx-1)*k+idx,1) = ar(idx,1);
end

% if a prior for the exogenous variables is selected put it in here:
for ii=1:n
    for jj=1:m
        beta0(k*ii-m+jj)=priorexo(ii,jj);
    end
end
% unvectorize (reshape) the vector to obtain the matrix B0
B0=reshape(beta0,k,n);

% next compute phi0, the variance-covariance matrix of beta, defined in (1.4.7)

% set first phi0 as a k*k matrix of zeros
phi0=zeros(k,k);

% set the variance for coefficients on lagged values, using (1.4.5)
for ii=1:n
    for jj=1:p
        phi0((jj-1)*n+ii,(jj-1)*n+ii)=(1/arvar(ii,1))*(lambda1/jj^lambda3)^2;
    end
end

% set the variance for exogenous variables, using (1.4.6)
for ii=1:m
    phi0(k-m+ii,k-m+ii)=(lambda1*lambda4(1,ii))^2;
end

% now compute alpha0 from (1.4.12)
alpha0=n+2;

% and finally compute S0, depending on which choice has been made for the prior ((1.4.13) or identity)
if prior==21
    S0=(alpha0-n-1)*diag(arvar);
elseif prior==22
    S0=eye(n);
else
end

author: european-central-bank

repo: BEAR-toolbox


text:
function surfacedata = ellipsoidsurface

a = 9;
b = 3;
c = 1;

a2 = a*a;
a4 = a2*a2;
b2 = b*b;
b4 = b2*b2;
c2 = c*c;
c4 = c2*c2;

surfacedata = struct('phi',@phi,'gradient', @gradient, 'unitoutnormal',...
    @unitoutnormal, 'project', @project, 'initmesh',@initmesh, 'meancurvature', @meancurvature, ...
    'u',@u, 'f',@f,'gra_u',@gra_u);

function [node, elem] = initmesh
M = load('meshdata/ellipsoid3394.mat');
node = M.node;
elem = M.elem;
end

function z = phi(p)
% level set function

x = p(:,1); y = p(:,2); z = p(:,3);
z = x.^2/a2 + y.^2/b2 + z.^2/c2 - 1.0;
end

function n = unitoutnormal(p)
%
n = gradient(p);
l = sqrt(sum(n.^2,2));
n = n./repmat(l,1,3);
end

function n = gradient(p)
% gradient of phi

x = p(:,1); y = p(:,2); z = p(:,3);
n = [2*x/a2, 2*y/b2, 2*z/c2];
end

function [node,dist] = project(p)
% projection function

s = sign(phi(p));

node = p;

normalAtNode = gradient(node);
valueAtNode = phi(node);
node = node - valueAtNode*ones(1,3).*normalAtNode./(dot(normalAtNode,normalAtNode,2)*ones(1,3));

vector = (-s)*ones(1,3).*(node - p);

d = s.*sqrt(dot(vector,vector,2));

normalAtNode = gradient(node);

node = p - d*ones(1,3).*normalAtNode./(sqrt(dot(normalAtNode,normalAtNode,2))*ones(1,3));

valueAtNode = phi(node);

normalAtNode = gradient(node);

vector = (-s)*ones(1,3).*(node - p);

d = s.*sqrt(dot(vector,vector,2));

e1 = normalAtNode./(sqrt(dot(normalAtNode, normalAtNode,2))*ones(1,3))-vector./(sqrt(dot(vector,vector,2))*ones(1,3));
error=sqrt(valueAtNode.^2./(dot(normalAtNode, normalAtNode,2))+dot(e1,e1,2));

k=1;
while max(abs(error)) > 1e-6 && k<200

    k=k+1;

    node = node - valueAtNode*ones(1,3).*normalAtNode./(dot(normalAtNode,normalAtNode,2)*ones(1,3));

    vector = -s*ones(1,3).*(node - p);
    d = s.*sqrt(dot(vector,vector,2));
    normalAtNode = gradient(node);

    node = p - d*ones(1,3).*normalAtNode./(sqrt(dot(normalAtNode,normalAtNode,2))*ones(1,3));

    valueAtNode = phi(node);

    normalAtNode = gradient(node);

    vector = (-s)*ones(1,3).*(node - p);

    d = s.*sqrt(dot(vector,vector,2));
    e1 = normalAtNode./(sqrt(dot(normalAtNode, normalAtNode,2))*ones(1,3))-vector./(sqrt(dot(vector,vector,2))*ones(1,3));
    error=sqrt(valueAtNode.^2./(dot(normalAtNode, normalAtNode,2))+dot(e1,e1,2));
end

if nargout == 2
    dist = d;
end
end

function mc = meancurvature(p)

x = p(:,1); y = p(:,2); z = p(:,3);
x2 = x.^2; y2 = y.^2; z2 = z.^2;

s1 = c4*a2*y2 + b4*a2*z2 + c4*b2*x2 + a4*b2*z2 + b4*c2*x2 + a4*c2*y2;
s2 = b4*c4*x2 + a4*c4*y2 + a4*b4*z2;
s3 = x2/a4+y2/b4+z2/c4;

mc = 0.5* s1./(s2.*sqrt(s3));
end

function z = f(p)
p = project(p);
x = p(:,1); y = p(:,2); z = p(:,3);
x2 = x.^2; y2 = y.^2; z2 = z.^2;
d1 = c4*y2 + b4*z2;
d2 = c2*y2 + b2*z2;
d3 = b4*c2 + c4*b2;
d4 = c4 * b4 * (d2 * a4 + d1 * a2 + x2 * d3);
d5 = (a4 * d1 + b4 * c4 * x2);
z = a2*( d4.*x.*cos(x) +  a2*d1.*d5.*sin(x)) ./ d5.^2;
end

function z = u(p)
% exact solution
p = project(p);
z =sin(p(:,1));
end

function z = gra_u(p)

p = project(p);
x = p(:,1); y = p(:,2); z = p(:,3);
x2 = x.^2; y2 = y.^2; z2 = z.^2;

d1 = x2/a4 + y2/b4 + z2/c4 ;

z1 = (1 - x2 ./(a4*d1));
z2 = -x.*y./(a2*b2*d1);
z3 = -x.*z./(a2*c2*d1);
z = cos(x)*ones(1,3).*[z1,z2,z3];
end

end

repo: ifem


text:
function dq = line2dquat(v,r)

% LINE2DQUAT  Transforms a line position expressed in vector notation into
%            its dual quaternion representation.
%
%     DQ = LINE2DQUAT(V,R) transforms the line position, specified by the
%       line orientation V and the coordinates of any point of the line, R,
%       into its dual quaternion representation. V and R must have the same
%       size. V  (resp. R) is either a vector of size 3 or an array of size
%       3*N (column i represents the orientation (resp. coordinates of a
%       point) of Line i) where N is the number of lines. DQ is either a
%       vector of size 8, either an array of size (8*N) depending on V
%       format. Each column of DQ represents the dual quaternion
%       representation of the corresponding line position.
%
% See also POS2DQUAT, VEL2DQUAT, LINEVEL2DQUAT, DQUAT2LINE

sv = size(v);
sr = size(r);
if sv == [1 3]
    v = v';
    sv = size(v);
end
if sr == [1 3]
    r = r';
    sr = size(r);
end

% check that v and r have the same size
if sv ~= sr
    error('DualQuaternion:line2dquat:sizesDoNotMatch',...
        'Arrays v and r should be the same size. Size of v matrix is [%d %d] while size of r matrix is [%d %d]',sv(1),sv(2),sr(1),sr(2));
end
% if the format is wrong
if sv(1) ~= 3
    error('DualQuaternion:line2dquat:wrongsize',...
        '%d rows in the V and R array. It should be 3. ',sv(1));
end

% normalization of the axis vector (if necessary)
n = length(v(1,:));
n2 = sum(v.^2).^0.5;
if max(n2)~=1 || min(n2)~=1
    n2 = repmat(n2,3,1);
    v =v./n2;
end

% construction of the line dual quaternion
dq = zeros(8,n);
dq(2:4,:) = v;
dq(6:8,:) = cross(r,v);

author: Sable

repo: mcbench-benchmarks


text:
function [FM, xf] = ft_preproc_online_filter_apply(FM, x)

% function [FM, xf] = ft_preproc_online_filter_apply(FM, x)
%
% Passes signal x (channels times samples) through the filter,
% returns updated filter model (delay states) and filtered signal.
%
% See also FT_PREPROC_ONLINE_FILTER_INIT

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

[dimX, numX] = size(x);

if dimX > numX
  % explicit algorithm - faster for many channels, 1 sample (=fMRI)
  xf = zeros(size(x));

  for k=1:numX;
    z_old = FM.z;
    z0 = x(:,k) - z_old*FM.A2;
    xf(:,k) = z_old * FM.B2 + z0*FM.B1;
    FM.z(:,2:end) = z_old(:,1:end-1);
    FM.z(:,1) = z0;
  end
else
  % use built-in MATLAB stuff - faster for many samples, few channels
  [xf, z] = filter(FM.B, FM.A, x, FM.z',2);
  FM.z = z';
end


text:
function test_int_margin_test ( )

%*****************************************************************************80
%
%% TEST_INT_MARGIN_TEST tests the TEST_INT_MARGIN library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp (  );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INT_MARGIN_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TEST_INT_MARGIN library.\n' );

  test_int_margin_test01 ( );
  test_int_margin_test02 ( );
  test_int_margin_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST_INT_MARGIN_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp (  );

  return
end


text:
function C = augment(A, alpha)
%AUGMENT  Augmented system matrix.
%         AUGMENT(A, ALPHA) is the square matrix
%         [ALPHA*EYE(m) A; A' ZEROS(n)] of dimension m+n, where A is m-by-n.
%         It is the symmetric and indefinite coefficient matrix of the
%         augmented system associated with a least squares problem
%         minimize NORM(A*x-b).  ALPHA defaults to 1.
%         Special case: if A is a scalar, n say, then AUGMENT(A) is the
%                       same as AUGMENT(RANDN(p,q)) where n = p+q and
%                       p = ROUND(n/2), that is, a random augmented matrix
%                       of dimension n is produced.
%         The eigenvalues of AUGMENT(A,ALPHA) are given in terms of the
%         singular values s(i) of A (where m>n) by
%           ALPHA/2 +/- SQRT( s(i)^2*ALPHA^2 + 1/4 ),  i=1:n  (2n eigenvalues),
%           ALPHA,  (m-n eigenvalues).
%         If m < n then the first expression provides 2m eigenvalues and the
%         remaining n-m eigenvalues are zero.
%
%         See also SPAUGMENT.

%         References:
%         G. H. Golub and C. F. Van Loan, Matrix Computations, third
%            Edition, Johns Hopkins University Press, Baltimore, Maryland,
%            1996; sec. 5.6.4.
%         N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%            Second edition, Society for Industrial and Applied Mathematics,
%            Philadelphia, PA, 2002; sec. 20.5.

[m, n] = size(A);
if nargin < 2, alpha = 1; end

if max(m,n) == 1
   n = A;
   p = round(n/2);
   q = n - p;
   A = randn(p,q);
   m = p; n = q;
end

C = [alpha*eye(m) A; A' zeros(n)];

author: opencobra

repo: cobratoolbox


text:
function mv = mv_tuningChangeMap(mv);
%
% mv = mv_tuningChangeMap(mv);
%
% For multi-voxel UI, create a map in which the value for each voxel
% measures the change in stimulus "tuning" (i.e., the N-dimensional
% vector of response amplitudes to each of the N selected stimulus
% conditions) between that voxel and its neighbors.
%
% Note that the determination of what a voxel's neighbors are depends on
% the view:
%   Inplane views: the neighbors are neighboring 6-connected voxels in the
%                  inplane data, independent of whether these neighbors
%                  are in the gray, white matter, csf, or other;
%   Gray views:    determined using the gray graph.
%
% Flat and Volume views are not yet supported.
%
% ras, 11/05.
if ieNotDefined('mv'), mv = get(gcf, 'UserData'); end

switch mv.roi.viewType
    case 'Inplane',
        nVoxels = size(mv.coords, 2);

        % get amplitudes for each voxel
        amps = mv_amps(mv);

        % initialize map vals
        vals = zeros(1, nVoxels);

        % create an 'offset' matrix describing where neighbors in
        % 6-connected data would be -- this will be useful for finding
        % neighbor coords in the main loop:
        offsets = [-1 0 0; 1 0 0; 0 -1 0; 0 1 0; 0 0 -1; 0 0 1]';

        %%%%%main loop
        hwait = mrvWaitbar(0, 'Computing Tuning Change Map...');
        for v = 1:nVoxels
            % find indices I of neighboring voxels
            pt = mv.coords(:, v);
            neighbors = repmat(pt, [1 6]) - offsets;
            [found I] = intersectCols(mv.coords, neighbors);

            if isempty(I)
                % no neighbors contained in the multi-voxel data:
                % set value for this point to -1:
                vals(v) = -1;
            else
                % get amplitudes for this voxel, neighbors
                A = amps(v, :); % amplitudes for this voxel
                B = amps(I, :); % amplitudes for neighbors

                % compute mean Euclidean distance between A and columns
                % in B:
                diff = [B - repmat(A, [1 size(B,2)])];
                dist = sqrt(sum(diff.^2));
                vals(v) = mean(dist);
            end

            mrvWaitbar(v/nVoxels, hwait);
        end
        close(hwait);

        % Create a map volume with the tuning values
        mrGlobals; loadSession;
        hI = initHiddenInplane(mv.params.dataType, mv.params.scans(1));
        mapvol = zeros(dataSize(hI));
        mapvol(roiIndices(hI, mv.coords)) = vals;

        % export as map
        hI.map = cell(1, numScans(hI));
        hI.map{mv.params.scans(1)} = mapvol;
        hI.mapName = 'Tuning_Change_Map';
        saveParameterMap(hI, [], 1);

    case 'Gray',
        nVoxels = size(mv.coords, 2);

        % get amplitudes for each voxel
        amps = mv_amps(mv);

        % initialize map vals
        vals = zeros(1, nVoxels);

        % create an 'offset' matrix describing where neighbors in
        % 6-connected data would be -- this will be useful for finding
        % neighbor coords in the main loop:
        offsets = [-1 0 0; 1 0 0; 0 -1 0; 0 1 0; 0 0 -1; 0 0 1]';

        %%%%%main loop
        hwait = mrvWaitbar(0, 'Computing Tuning Change Map...');
        for v = 1:nVoxels
            % find indices I of neighboring voxels
            pt = mv.coords(:, v);
            neighbors = repmat(pt, [1 6]) - offsets;
            [found I] = intersectCols(mv.coords, neighbors);

            if isempty(I)
                % no neighbors contained in the multi-voxel data:
                % set value for this point to -1:
                vals(v) = -1;
            else
                % get amplitudes for this voxel, neighbors
                A = amps(v, :); % amplitudes for this voxel
                B = amps(I, :); % amplitudes for neighbors

                % compute mean Euclidean distance between A and columns
                % in B:
                diff = [B - repmat(A, [size(B,1) 1])];
                dist = sqrt(sum(diff.^2, 2));
                vals(v) = mean(dist(:));
            end

            mrvWaitbar(v/nVoxels, hwait);
        end
        close(hwait);

        % Create a map volume with the tuning values
        mrGlobals; loadSession;
        hG = initHiddenGray(mv.params.dataType, mv.params.scans(1));
        mapvol = zeros(dataSize(hG));
        mapvol(roiIndices(hG, mv.coords)) = vals;

        % export as map
        hG.map = cell(1, numScans(hG));
        hG.map{mv.params.scans(1)} = mapvol;
        hG.mapName = 'Tuning_Change_Map';
        saveParameterMap(hG, [], 1, 0);

    otherwise, % not yet supported
        error('Sorry, this view type is not yet supported.');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [model, specificData, coreRxnAbbr, modelGenerationReport] = growthMediaToModel(model, specificData, param, coreRxnAbbr, modelGenerationReport)
%
% USAGE:
%   [model, specificData, coreRxnAbbr, modelGenerationReport] = growthMediaToModel(model, specificData, param, coreRxnAbbr, modelGenerationReport)
%
% INPUTS:
%  model.rxns:
%  model.rxnNames:
%  model.lb:
%  model.ub:
%  specificData.mediaData:
%  specificData.mediaData.mets:
%  specificData.mediaData.rxns:
%  specificData.exoMet:
%
%  param.metabolomicsBeforeExtraction:
%  param.debug:
%  param.workingDirectory:
%  param.printLevel:
%  param.TolMinBoundary:
%  param.TolMaxBoundary:
%  param.relaxOptions:
%  coreRxnAbbr:
%  modelGenerationReport:
%
% OUTPUTS:
%  model.rxns:
%  model.rxnNames:
%  model.lb:
%  model.ub:
%  specificData.mediaData:
%  specificData.mediaData.mediaData:
%  specificData.mediaData.mets:
%  specificData.mediaData.Properties:
%  specificData.mediaData.rxns:
%  specificData.exoMet:
%  specificData.exoMet.exoMet:
%  specificData.exoMet.ismedia:
%  coreRxnAbbr:
%  modelGenerationReport:
%
% EXAMPLE:
%
% NOTE:
%
% Author(s): Aga W.

if param.metabolomicsBeforeExtraction && param.debug
    save([param.workingDirectory filesep '10.a.debug_prior_to_media_constraints.mat'])
elseif param.debug
    save([param.workingDirectory filesep '22.a.debug_prior_to_media_constraints.mat'])
end

%% 10/22.a. Set growth media constraints (if provided)

if isfield(specificData, 'mediaData') && ~isempty(specificData.mediaData)

    if param.printLevel > 0
        disp('--------------------------------------------------------------')
        disp(' ')
        disp('Adding growth media information...')
    end

    % Constrain media reactions
    if ~ismember('rxns', specificData.mediaData.Properties.VariableNames) && ismember('mets', specificData.mediaData.Properties.VariableNames)
        modelTemp = findSExRxnInd(model);
        allExRxns = model.rxns(modelTemp.ExchRxnBool);
        for i = 1:length(specificData.mediaData.mets)
            if any(contains(allExRxns, strcat('_', specificData.mediaData.mets(i))))
                specificData.mediaData.rxns(i) = allExRxns(contains(allExRxns, strcat('_', specificData.mediaData.mets(i))));
            end
        end
        [model, rxnsConstrained, rxnBoundsCorrected] = constrainRxns(model, specificData, param, 'mediaDataConstraints', param.printLevel);
    elseif ismember('rxns', specificData.mediaData.Properties.VariableNames)
        [model, rxnsConstrained, rxnBoundsCorrected] = constrainRxns(model, specificData, param, 'mediaDataConstraints', param.printLevel);
    end

    % Identify media reactions in metabolomic data
    if isfield(specificData,'exoMet')
        specificData.exoMet.ismedia = false(length(specificData.exoMet.mean), 1);
        specificData.exoMet.ismedia(ismember(specificData.exoMet.rxns, specificData.mediaData.rxns)) = 1;
    end

    if param.printLevel > 0
        disp(' ')
    end
    clear modelTemp

    if ~isempty(rxnsConstrained)
        if param.printLevel > 1
            fprintf('%s\n\n','Growth media bounds were set on the following reactions:')
            rxnsConstrainedBool = ismember(model.rxns, rxnsConstrained);
            table(model.rxns(rxnsConstrainedBool), model.rxnNames(rxnsConstrainedBool),...
                model.lb(rxnsConstrainedBool), model.ub(rxnsConstrainedBool), ...
                printRxnFormula(model, 'rxnAbbrList', model.rxns(rxnsConstrainedBool), ...
                'printFlag', false), 'VariableNames', {'rxnsConstrained', ...
                'Name', 'lb', 'ub', 'equation'})
            %printConstraints(model, options.TolMinBoundary, options.TolMaxBoundary, ismember(model.rxns, rxnsConstrained))
        end
    end

    if ~isempty(rxnBoundsCorrected)
        if param.printLevel > 1
            fprintf('%s\n\n','Bounds were corrected by constrainRxns on the following reactions:')
            printConstraints(model, param.TolMinBoundary, param.TolMaxBoundary, ismember(model.rxns, rxnBoundsCorrected))
        end
    end

    if any(model.lb > model.ub)
        error('lower bounds greater than upper bounds')
    end
end

% Check feasibility
if isfield(specificData, 'mediaData') && ~isempty(specificData.mediaData)
    sol = optimizeCbModel(model);
    if  sol.stat ~=1
        fprintf('%s\n','Infeasible after application of growth media constraints. Trying relaxation...')
        %options.relaxOptions.rxns = unique([options.exoMet.rxns; options.mediaData.rxns]);
        [solution, modelTemp] = relaxedFBA(model, param.relaxOptions);
        if solution.stat==1
            fprintf('%s\n','.. relaxation worked.')
            model = modelTemp;
        else
            error('Infeasible after application of media constraints and relaxation failed.')
        end
    elseif sol.stat ==1
        if param.printLevel>0
            disp(' ')
            fprintf('%s\n\n','Feasible after application of media constraints.')
        end
    end
end


text:
function test_suite=test_moxunit_isa_test_skipped_exception
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_moxunit_isa_test_skipped_exception_basics()
    try
        moxunit_throw_test_skipped_exception('foo');
        assert(false,'should not come here');
    catch
        e=lasterror();
    end

    assert_matches(true,e.identifier);
    assert_matches(false,e.identifier(2:end));
    assert_matches(false,'skipped');

function assert_matches(tf,x)
    e=struct();
    e.identifier=x;

    assertEqual(tf,moxunit_isa_test_skipped_exception(e))

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
black=[0 0 0];white=[255 255 255];

Screen('FillRect', w, black, [wRect(3)/2-(s_width/2) wRect(1) wRect(3)/2+(s_width/2) wRect(1)+s_height]); % up
Screen('FillRect', w, black, [wRect(3)/2-(s_width/2) wRect(4)-s_height wRect(3)/2+(s_width/2) wRect(4)]); % down
Screen('FillRect', w, black, [wRect(1)  wRect(4)/2-(s_height/2)  wRect(1)+s_width   wRect(4)/2+(s_height/2)]); % left
Screen('FillRect',w,black,[wRect(3)-s_width wRect(4)/2-(s_height/2) wRect(3) wRect(4)/2+(s_height/2)]); % right
fixSize =80;
[X,Y] = RectCenter(wRect);
PixofFixationSize = 80;
FixCross = [X-1,Y-PixofFixationSize,X+1,Y+PixofFixationSize;X-PixofFixationSize,Y-1,X+PixofFixationSize,Y+1];
Screen('FillRect', w, black, FixCross');


text:
function [dpsi, deps] = nut2000_lp (t)

% low precison nutation based on iau 2000a

% this function evaluates a short nutation series and returns approximate
% values for nutation in longitude and nutation in obliquity for a given
% tdb julian date. in this mode, only the largest 13 terms of the iau 2000a
% nutation series are evaluated.

% input

%  t = tdb time in julian centuries since j2000.0

% output

%  dpsi = nutation in longitude in arcseconds

%  deps = nutation in obliquity in arcseconds

% note: in low-accuracy mode, max error in dpsi < 0.05 arcsec,
% max error in deps < 0.02 arcsec, average error about 1/4 of max.

% ported from NOVAS 3.0

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

% largest 13 terms of iau 2000a nutation series, with precision
% of coefficients truncated

x = [0.0, 0.0, 0.0, 0.0, 1.0, -17.2064,-0.01747,  9.2052,  0.00091; ...
     0.0, 0.0, 2.0,-2.0, 2.0, -1.3171, -0.00017,  0.5730, -0.00030; ...
     0.0, 0.0, 2.0, 0.0, 2.0, -0.2276, -0.00002,  0.0978, -0.00005; ...
     0.0, 0.0, 0.0, 0.0, 2.0,  0.2075,  0.00002, -0.0897,  0.00005; ...
     0.0, 1.0, 0.0, 0.0, 0.0,  0.1476, -0.00036,  0.0074, -0.00002; ...
     0.0, 1.0, 2.0,-2.0, 2.0, -0.0517,  0.00012,  0.0224, -0.00007; ...
     1.0, 0.0, 0.0, 0.0, 0.0,  0.0711,  0.00001, -0.0007,  0.00000; ...
     0.0, 0.0, 2.0, 0.0, 1.0, -0.0387, -0.00004,  0.0201,  0.00000; ...
     1.0, 0.0, 2.0, 0.0, 2.0, -0.0301,  0.00000,  0.0129, -0.00001; ...
     0.0,-1.0, 2.0,-2.0, 2.0,  0.0216, -0.00005, -0.0096,  0.00003; ...
     0.0, 0.0, 2.0,-2.0, 1.0,  0.0128,  0.00001, -0.0069, -0.00000; ...
    -1.0, 0.0, 2.0, 0.0, 2.0,  0.0123,  0.00000, -0.0053,  0.00000; ...
    -1.0, 0.0, 0.0, 2.0, 0.0,  0.0157,  0.00000, -0.0001,  0.00000];

% transpose x matrix

x = x';

% ----------------------------------------------------
% remaining terms all have amplitudes < 0.01 arcsecond
% ----------------------------------------------------

% computation of fundamental arguments

[el, elp, f, d, om] = funarg (t);

dpsi = 0.0d0;

deps = 0.0d0;

% sum nutation series terms

for i = 13:-1:1

    arg = x(1, i) * el ...
        + x(2, i) * elp ...
        + x(3, i) * f ...
        + x(4, i) * d ...
        + x(5, i) * om;

    dpsi = (x(6, i) + x(7, i) * t) * sin(arg) + dpsi;

    deps = (x(8, i) + x(9, i) * t) * cos(arg) + deps;

end

% add in out-of-phase component of principal (18.6-year) term
% (to avoid small but long-term bias in results)

dpsi = dpsi + 0.0033d0 * cos(om);

deps = deps + 0.0015d0 * sin(om);

author: goGPS-Project


text:
function [ok, pns_norm, pns_comp, t_axis]=calcPNS(obj,hardware,doPlots)
% calculate PNS using safe model implementation by Szczepankiewicz and Witzel
% assumes safe_pns_prediction package has been downloaded and installed in
% Matlab path. See http://github.com/filip-szczepankiewicz/safe_pns_prediction
%
% returns pns levels due to respective axes (normalized to 1 and not to 100%)
%
% inputs:
%     hardware - hardware specifications. see safe_example_hw() from
%                the safe_pns_prediction package. Alternatively a text file
%                in the .asc format (Siemens) can be passed, e.g. for Prisma
%                it is MP_GPA_K2309_2250V_951A_AS82.asc (we leave it as an
%                exercise to the interested user to find were these files
%                can be acquired from);
%     doPlots  - optional parameter (defaluts to true)

if nargin < 3
    doPlots=true;
end

% acquire the entire gradient wave form
gw=obj.waveforms_and_times();
if doPlots
    figure;
    plot(gw{1}(1,:),gw{1}(2,:),gw{2}(1,:),gw{2}(2,:),gw{3}(1,:),gw{3}(2,:)); % plot the entire gradient shape
    title('gradient wave form, in T/m');
end

% find beginning and end times and resample GWs to a regular sampling raster
tf=[];
tl=[];
for i=1:3
    if size(gw{i},2)>0
        tf(end+1)=gw{i}(1,1);
        tl(end+1)=gw{i}(1,end);
    end
end
nt_min=floor(min(tf)/obj.gradRasterTime+eps);
nt_max=ceil(max(tl)/obj.gradRasterTime-eps);
% shift raster positions to the centers of the raster periods
nt_min = nt_min + 0.5;
nt_max = nt_max - 0.5;
if (nt_min<0.5)
    nt_min=0.5
end
t_axis=(nt_min:nt_max)*obj.gradRasterTime;
gwr=zeros(length(t_axis),3);
for i=1:3
    if size(gw{i},2)>0
        gwr(:,i)=interp1(gw{i}(1,:),gw{i}(2,:),t_axis,'linear',0);
    end
end

if ischar(hardware)
    % this loads the parameters from the provided text file
    asc=mr.Siemens.readasc(hardware);
    hardware=asc_to_hw(asc);
end

% use the Szczepankiewicz' and Witzel's implementation
[pns_comp,res]=safe_gwf_to_pns(gwr/obj.sys.gamma, NaN*ones(length(t_axis),1), obj.gradRasterTime, hardware); % the RF vector is unused in the code inside but it is zeropaded and exported ...
% use the exported RF vector to detect and undo zerpopadding
pns_comp=0.01*pns_comp(~isfinite(res.rf),:)';
% calc pns_norm and the final ok/not_ok
pns_norm=vecnorm(pns_comp);
ok=all(pns_norm<1);
% ready
if doPlots
    % plot results
    figure;safe_plot(pns_comp'*100, obj.gradRasterTime);
end

end

% local utility functions

function hw = asc_to_hw(asc)
% function hw = asc_to_hw(asc)
%
% SAFE model parameters for the asc structure as read from the asc file.
% See comments for units.
%
% Maxim Zaitsev 08/10/2019

if isfield(asc,'asCOMP') && isfield(asc.asCOMP,'tName')
    hw.name          = asc.asCOMP(1).tName;
else
    hw.name          = 'unknown';
end
%hw.look_ahead    =  1.0; % MZ: this is not a real hardware parameter but a coefficient, with which the final result is multiplied

hw.x.tau1        = asc.flGSWDTauX(1);  % ms
hw.x.tau2        = asc.flGSWDTauX(2);  % ms
hw.x.tau3        = asc.flGSWDTauX(3);  % ms
hw.x.a1          = asc.flGSWDAX(1);
hw.x.a2          = asc.flGSWDAX(2);
hw.x.a3          = asc.flGSWDAX(3);
hw.x.stim_limit  = asc.flGSWDStimulationLimitX;  % T/m/s
hw.x.stim_thresh = asc.flGSWDStimulationThresholdX;  % T/m/s
hw.x.g_scale     = asc.asGPAParameters.sGCParameters.flGScaleFactorX;

hw.y.tau1        = asc.flGSWDTauY(1);  % ms
hw.y.tau2        = asc.flGSWDTauY(2);  % ms
hw.y.tau3        = asc.flGSWDTauY(3);  % ms
hw.y.a1          = asc.flGSWDAY(1);
hw.y.a2          = asc.flGSWDAY(2);
hw.y.a3          = asc.flGSWDAY(3);
hw.y.stim_limit  = asc.flGSWDStimulationLimitY;  % T/m/s
hw.y.stim_thresh = asc.flGSWDStimulationThresholdY;  % T/m/s
hw.y.g_scale     = asc.asGPAParameters.sGCParameters.flGScaleFactorY;

hw.z.tau1        = asc.flGSWDTauZ(1);  % ms
hw.z.tau2        = asc.flGSWDTauZ(2);  % ms
hw.z.tau3        = asc.flGSWDTauZ(3);  % ms
hw.z.a1          = asc.flGSWDAZ(1);
hw.z.a2          = asc.flGSWDAZ(2);
hw.z.a3          = asc.flGSWDAZ(3);
hw.z.stim_limit  = asc.flGSWDStimulationLimitZ;  % T/m/s
hw.z.stim_thresh = asc.flGSWDStimulationThresholdZ;  % T/m/s
hw.z.g_scale     = asc.asGPAParameters.sGCParameters.flGScaleFactorZ;
end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function RMP = learnRMP(Problem,SubPopulation)

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

   Population = [SubPopulation{:}];
    for i = 1 : length(Problem.SubD)
        subpops(i).data = [];
        vars(i)         = Problem.SubD(i);
    end
    for i = 1 : length(Population)
        subpops(Population(i).dec(end)).data = [subpops(Population(i).dec(end)).data;Population(i).dec(1:Problem.SubD(1))];
    end
    RMP = learnRMP_sub(subpops,vars);
end

function rmpMatrix = learnRMP_sub(subpops,vars)
% There are two inputs. subpops(i).data corresponds to the population
% corresponding to the ith task; vars(i) is the number of design variables
% of the ith task.

    numtasks  = length(subpops);
    maxdims   = max(vars);
    rmpMatrix = eye(numtasks);
    % Add noise and Build probabilistic models
    for i = 1 : numtasks
        probmodel(i).nsamples = size(subpops(i).data,1);
        nrandsamples          = floor(0.1*probmodel(i).nsamples);
        randMat               = rand(nrandsamples,maxdims);
        probmodel(i).mean     = mean([subpops(i).data;randMat]);	% Univariate distribution mean
        probmodel(i).stdev    = std([subpops(i).data;randMat]);     % Univariate distribution standard deviation
    end
    for i = 1 : numtasks
        for j = i+1 : numtasks
            popdata(1).probmatrix = ones(probmodel(i).nsamples,2);
            popdata(2).probmatrix = ones(probmodel(j).nsamples,2);
            dims = min([vars(i),vars(j)]);
            for k = 1 : probmodel(i).nsamples
                for l = 1 : dims
                    popdata(1).probmatrix(k,1) = popdata(1).probmatrix(k,1)*pdf('Normal',subpops(i).data(k,l),probmodel(i).mean(l),probmodel(i).stdev(l));
                    popdata(1).probmatrix(k,2) = popdata(1).probmatrix(k,2)*pdf('Normal',subpops(i).data(k,l),probmodel(j).mean(l),probmodel(j).stdev(l));
                end
            end
            for k = 1 : probmodel(j).nsamples
                for l = 1 : dims
                    popdata(2).probmatrix(k,1) = popdata(2).probmatrix(k,1)*pdf('Normal',subpops(j).data(k,l),probmodel(i).mean(l),probmodel(i).stdev(l));
                    popdata(2).probmatrix(k,2) = popdata(2).probmatrix(k,2)*pdf('Normal',subpops(j).data(k,l),probmodel(j).mean(l),probmodel(j).stdev(l));
                end
            end
            rmpMatrix(i,j) = max([0,fminbnd(@(x)loglik(x,popdata,numtasks),0,1)+normrnd(0,0.01)]);  %fminbnd(@(x)loglik(x,popdata,numtasks),0,1)
            rmpMatrix(i,j) = min(rmpMatrix(i,j),1);
            rmpMatrix(j,i) = rmpMatrix(i,j);
        end
    end
end

function f = loglik(rmp,popdata,ntasks)
    f = 0;
    for i = 1 : 2
        for j = 1 : 2
            if i == j
                popdata(i).probmatrix(:,j) = popdata(i).probmatrix(:,j)*(1-(0.5*(ntasks-1)*rmp/ntasks));
            else
                popdata(i).probmatrix(:,j) = popdata(i).probmatrix(:,j)*0.5*(ntasks-1)*rmp/ntasks;
            end
        end
        f = f + sum(-log(sum(popdata(i).probmatrix,2)));
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function []=strsdispols(decimaldates1,stringdates1,strshocks_estimates,endo,pref,IRFt,strctident)

if IRFt==1||IRFt==2||IRFt==3
    identified=size(endo,1); % fully identified
    labels=endo; % simply use the name of the endogenous variables as labels
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions (+ IV)
    identified=size(strctident.signreslabels_shocks,1); % count the labels provided in the sign res sheet (+ IV)
    labels=strctident.signreslabels_shocks; % signreslabels, only identified shocks
elseif IRFt==5
    identified=1; % one IV shock
    labels{identified,1}=strcat('IV Shock (',strctident.Instrument,')'); % one label for the IV shock
end

if pref.plot==1
% create shock figure
strshocks=figure('Tag','BEARresults');
set(strshocks,'Color',[0.9 0.9 0.9]);
set(strshocks,'name','structural shocks');
ncolumns=ceil(identified^0.5);
nrows=ceil(identified/ncolumns);
for ii=1:identified
subplot(nrows,ncolumns,ii);
hold on
plot(decimaldates1,strshocks_estimates(ii,:),'Color',[0.4 0.4 1],'LineWidth',2);
plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
hold off
minband=min(strshocks_estimates(ii,:));
maxband=max(strshocks_estimates(ii,:));
space=maxband-minband;
Ymin=minband-0.2*space;
Ymax=maxband+0.2*space;
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'YLim',[Ymin,Ymax],'FontName','Times New Roman');
title(labels{ii,1},'FontName','Times New Roman','FontSize',10,'FontWeight','normal','interpreter','none');
end
end

if pref.results==1
% finally, save in excel
% create the cell that will be saved on excel
horzspace=repmat({''},1,identified);
strshockcell=[{'structural shocks'} horzspace;{''} horzspace;{''} labels';stringdates1 num2cell(strshocks_estimates')];

% write in excel
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strshockcell,'structshocks','B2');
end


text:
function structColorGUI(command,varargin)
%function structColorGUI(command,varargin)
%GUI to change structure color
%
%APA, 01/25/08, based on structColorGUI
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS

indexS = planC{end};

switch upper(command)

    case 'INIT'

        units = 'pixels';

        x = 340; y = 700;

        hFig = figure('Menu','None','Position',[100,100,x,y],'Name','Reassign Color', 'units', 'pixels',...
            'NumberTitle', 'off', 'resize', 'off', 'Tag', 'structColorGUI','WindowButtonDownFcn','structColorGUI(''BTDWN'')',...
            'WindowButtonUpFcn','structColorGUI(''BTUP'')');

        statusFrame = uicontrol(hFig, 'units',units ,'Position', [10 10 x-18 y-18], 'style', 'frame','HitTest','off');
        frameColor = get(statusFrame,'backgroundcolor');
        ud.frameColor = frameColor;

        %frameColor = [0.9 0.9 0.9];

        uicontrol(hFig, 'units',units,'Position',[x/2-100,670,200,20], 'style', 'text', 'String', 'Reassign Structure Color',...
            'backgroundcolor', frameColor,'fontweight', 'bold','fontsize',10);

        uicontrol(hFig, 'units',units,'Position',[12,650,300,20], 'style', 'text', 'String', '(Hold color, Drag, Drop over another to swap)',...
            'backgroundcolor', frameColor,'fontweight', 'normal','fontsize',8);

        uicontrol(hFig, 'units',units,'Position',[10+40,630,150,20], 'style', 'text', 'String', 'Current Color',...
            'backgroundcolor', frameColor,'fontweight', 'bold','fontsize',9);

        uicontrol(hFig, 'units',units,'Position',[10+180,630,90,20], 'style', 'text', 'String', 'New Color',...
            'backgroundcolor', frameColor,'fontweight', 'bold','fontsize',9);

        numRows = 28;

        for i=1:numRows
            ud.handles.strNum(i) = uicontrol(hFig, 'units',units,'Position',[10+3    y-18-50-i*20   20   15], 'style', 'text',...
                'String', num2str(i), 'backgroundcolor', frameColor,'Visible','Off');
            ud.handles.strName(i)  = uicontrol(hFig, 'units',units,'Position',[10+40   y-18-50-i*20   150  15], 'style', 'text',...
                'String', 'struct Name', 'backgroundcolor', frameColor,'Visible','Off');
            ud.handles.strSelectColor(i)  = uicontrol(hFig, 'units',units,'Position',[10+240   y-18-50-i*20   40  15], 'style', 'push',...
                'String', 'Pick', 'backgroundcolor', frameColor,'Visible','Off', 'callback',['structColorGUI(''pickColor'',''' num2str(i) ''')']);
            ud.handles.newName(i) = uicontrol(hFig, 'units',units,'Position',[10+210  y-18-50-i*20   20  15], 'style', 'text', 'backgroundcolor', stateS.optS.colorOrder(i,:),...
                'horizontalAlignment', 'center', 'userdata', i, 'Tag', 'newStructName','Visible','Off','enable','off');
            ud.handles.newNamePos{i} = get(ud.handles.newName(i),'position');
        end

        ud.numStructRows = numRows;

        %Get structure numbers associated with the passed scan
        scanNum = str2num(varargin{1});
        [assocScansV, relStructNumV] = getStructureAssociatedScan(1:length(planC{indexS.structures}), planC);
        assocStructs = find(assocScansV == scanNum);
        ud.assocStructs = assocStructs;

        for j =  1:length(ud.assocStructs)
            ud.struct(j).newName = '';
        end

        ud.handles.structsSlider = uicontrol(hFig, 'units',units,'Position',[x-18-10 65 20 y-10-125], 'style', 'slider',...
            'enable', 'off','tag','structSlider','callback','structColorGUI(''sliderRefresh'')','backgroundcolor', ud.frameColor);

%         ud.handles.undo = uicontrol(hFig, 'units',units, 'Position', [30 15 60 20], 'style', 'pushbutton', ...
%             'string', 'Undo', 'callback', 'structColorGUI(''undo'')');

        ud.handles.reset = uicontrol(hFig, 'units',units, 'Position', [100 15 60 20], 'style', 'pushbutton', ...
            'string', 'Reset', 'callback', 'structColorGUI(''reset'')');

        ud.handles.save = uicontrol(hFig, 'units',units, 'Position', [170 15 60 20], 'style', 'pushbutton', ...
            'string', 'Save', 'callback', 'structColorGUI(''save'')');

        ud.handles.exit = uicontrol(hFig, 'units',units, 'Position', [240 15 60 20], 'style', 'pushbutton', ...
            'string', 'Exit', 'callback', 'structColorGUI(''exit'')');

        %Store structure colors in ud
        strColorM = stateS.optS.colorOrder;
        numRepeat = ceil(length(ud.assocStructs)/size(strColorM,1));
        strColorM = repmat(strColorM,[numRepeat 1]);
        for j =  1:length(ud.assocStructs)
            strColorM(j,:) = planC{indexS.structures}(ud.assocStructs(j)).structureColor;
        end
        ud.strColorM = strColorM;
        ud.strColorOldM = strColorM;

        ud.handles.currentColor = [];

        structLen = length(ud.assocStructs);
        if (structLen/28) < 2
            max = 2;
        elseif (structLen/28) > 2
            max = 3;
        elseif (structLen/28) > 3
            max = 4;
        end

        try
            value = ud.sliderValue;
        catch
            value = max;
            ud.slider.oldValue = value;
        end
        if max == 2
            sliderstep = 1;
        else
            sliderstep = 1/max;
        end
        set(ud.handles.structsSlider,'min',1,'max',max,'value',max,'sliderstep',[sliderstep sliderstep]);

        %Log
        ud.swapHistory = [];

        set(hFig,'Userdata',ud);

        structColorGUI('refresh')

    case 'REFRESH'

        hFig = findobj('Tag','structColorGUI');

        ud = get(hFig,'Userdata');

        %structLen = length(planC{indexS.structures});
        structLen = length(ud.assocStructs);

        if structLen > ud.numStructRows
            % Enable Slider
            set(ud.handles.structsSlider,'Enable','On');
        end

        if structLen > 28
            set(ud.handles.structsSlider,'enable', 'on');

            % initialize structure slider
            set(hFig,'Userdata',ud);

            structColorGUI('sliderRefresh','init');

        else
            if length(ud.assocStructs)> 28
                numStruct = 28;
            else
                numStruct = length(ud.assocStructs);
            end

            for i=1:numStruct
                structNum  = i;
                structName = planC{indexS.structures}(ud.assocStructs(structNum)).structureName;
                set(ud.handles.strName(structNum), 'string', structName, 'visible', 'on','foregroundcolor',ud.strColorOldM(i,:));
                set(ud.handles.strSelectColor(structNum),'visible', 'on')
                set(ud.handles.strNum(structNum) , 'visible','on');
                set(ud.handles.newName(structNum) , 'visible','on','backgroundcolor',ud.strColorM(i,:));
            end

            %Provide Additional color-Selections
            if numStruct<28
                for i=numStruct+1:28
                    set(ud.handles.newName(i) , 'visible','on','backgroundcolor',ud.strColorM(i,:));
                end
            end

        end

    case 'SLIDERREFRESH'
        ud = get(findobj('Tag','structColorGUI'),'Userdata');

        value = round(get(ud.handles.structsSlider,'value'));
        set(ud.handles.structsSlider,'value',value)
        ud.slider.oldValue = value;
        set(findobj('Tag','structColorGUI'),'Userdata',ud);
        assocStructs = ud.assocStructs;

        numStruct = getStructDispLen(value,assocStructs);

        for i = 1:length(numStruct)
            structNum  = assocStructs(numStruct(i));
            structName = planC{indexS.structures}(structNum).structureName;
            set(ud.handles.strName(i), 'string', structName, 'visible', 'on','foregroundcolor',ud.strColorOldM(numStruct(i),:));
            set(ud.handles.strSelectColor(i),'visible', 'on')
            set(ud.handles.strNum(i) , 'visible','on','String',structNum);
            set(ud.handles.newName(i) ,'visible','on','backgroundcolor',ud.strColorM(numStruct(i),:));
        end

        for j = length(numStruct)+1:28
            set(ud.handles.strName(j), 'visible','off');
            set(ud.handles.strSelectColor(j),'visible', 'off')
            set(ud.handles.strNum(j) , 'visible','off');
            set(ud.handles.newName(j), 'visible','off');
        end

        %Provide Additional color-Selections
        for i = length(numStruct)+1:28
            set(ud.handles.newName(i) , 'visible','on','backgroundcolor',stateS.optS.colorOrder(i,:));
        end

    case 'BTDWN'
        % hFig = findobj('Tag','structColorGUI');
        hFig = gcf;
        ud = get(hFig,'Userdata');
        %Get the color object (static text) clicked
        currentPt = get(hFig,'currentPoint');
        hObject = [];
        for i=1:length(ud.handles.newName)
            pos = get(ud.handles.newName(i),'position');
            if currentPt(1)>=pos(1) && currentPt(1)<=pos(1)+pos(3) && currentPt(2)>=pos(2) && currentPt(2)<=pos(2)+pos(4)
                hObject = ud.handles.newName(i);
                break
            end
        end
        if ~isempty(hObject)
            ud.handles.currentColor = hObject;
            ud.handles.currentColorPos = get(hObject,'position');
            indColor = find(ud.handles.newName==hObject);
            set(ud.handles.strName,'FontWeight','normal')
            set(ud.handles.strName(indColor),'FontWeight','bold')
            set(ud.handles.strName(indColor),'backgroundcolor', 'y')
            set(hFig,'Userdata',ud);
            set(hFig,'WindowButtonMotionFcn','structColorGUI(''BTMTN'')')
        end

        %Set the structure name string font to boldface

    case 'BTMTN'
        hFig = gcf;
        ud = get(hFig,'Userdata');
        set(ud.handles.strName,'fontWeight','normal')
        set(ud.handles.strName,'backgroundcolor', ud.frameColor)
        movingColorInd = find(ud.handles.newName == ud.handles.currentColor);
        set(ud.handles.strName(movingColorInd),'fontWeight','bold')
        set(ud.handles.strName(movingColorInd),'backgroundcolor', 'y')
        %Get the structure object under cursor
        %         hFig = findobj('Tag','structColorGUI');
        currentPt = get(hFig,'currentPoint');
        set(ud.handles.currentColor,'position',[currentPt 20 15])
        %Set the structure name string font to boldface
        %Get the structure object under cursor
        for i=1:length(ud.handles.newName)
            colorPos = ud.handles.newNamePos{i};
            if (currentPt(1) >= colorPos(1)) && (currentPt(1) <= colorPos(1)+colorPos(3)) && (currentPt(2) >= colorPos(2)) && (currentPt(2) <= colorPos(2)+colorPos(4))
                baseColorInd = i;
                set(ud.handles.strName(baseColorInd),'fontWeight','bold')
                set(ud.handles.strName(baseColorInd),'backgroundcolor', 'y')
                drawnow
                break
            end
        end

    case 'BTUP'
        % hFig = findobj('Tag','structColorGUI');
        hFig = gcf;
        ud = get(hFig,'Userdata');
        set(ud.handles.strName,'fontWeight','normal')
        set(ud.handles.strName,'backgroundcolor', ud.frameColor)
        if isempty(ud.handles.currentColor)
            return;
        end
        %Get the structure object under cursor
        movingColorInd = find(ud.handles.newName == ud.handles.currentColor);
        currentPt = get(hFig,'currentPoint');
        if strcmpi(get(ud.handles.structsSlider ,'enable'),'off')
            sliderVal = 1;
            sliderMax = 1;
        else
            sliderVal = get(ud.handles.structsSlider ,'value');
            sliderMax = get(ud.handles.structsSlider ,'max');
        end
        val = sliderMax - sliderVal + 1;
        movingColorInd = (val-1)*ud.numStructRows + movingColorInd;
        baseColorInd = [];
        for i=1:length(ud.handles.newName)
            colorPos = ud.handles.newNamePos{i};
            if (currentPt(1) >= colorPos(1)) && (currentPt(1) <= colorPos(1)+colorPos(3)) && (currentPt(2) >= colorPos(2)) && (currentPt(2) <= colorPos(2)+colorPos(4))
                baseColorInd = (val-1)*ud.numStructRows + i;
                break
            end
        end

        %Swap the structure colors
        if ~isempty(baseColorInd)
            baseColorTmp = ud.strColorM(baseColorInd,:);
            ud.strColorM(baseColorInd,:) = ud.strColorM(movingColorInd,:);
            ud.strColorM(movingColorInd,:) = baseColorTmp;
            ud.swapHistory = [ud.swapHistory; movingColorInd baseColorInd];
        end

        set(ud.handles.currentColor,'position',ud.handles.currentColorPos)
        set(hFig,'WindowButtonMotionFcn','');

        set(hFig,'Userdata',ud);
        if strcmp(get(ud.handles.structsSlider,'Enable'),'on')
            structColorGUI('sliderrefresh')
        else
            structColorGUI('refresh')
        end

    case 'PICKCOLOR'
        % hFig = findobj('Tag','structColorGUI');
        hFig = gcf;
        ud = get(hFig,'Userdata');
        %Get structure index
        if strcmpi(get(ud.handles.structsSlider ,'enable'),'off')
            sliderVal = 1;
            sliderMax = 1;
        else
            sliderVal = get(ud.handles.structsSlider ,'value');
            sliderMax = get(ud.handles.structsSlider ,'max');
        end
        val = sliderMax - sliderVal + 1;
        structIndex = (val-1)*ud.numStructRows + str2num(varargin{1});
        drawnow;
        newColorVal = uisetcolor;
        drawnow;
        prevColorVal = ud.strColorM(structIndex,:);
        ud.strColorM(structIndex,:) = newColorVal;
        set(hFig,'Userdata',ud);
        if strcmp(get(ud.handles.structsSlider,'Enable'),'on')
            structColorGUI('sliderrefresh')
        else
            structColorGUI('refresh')
        end

    case 'UNDO'
        hFig = findobj('Tag','structColorGUI');
        ud = get(hFig,'Userdata');
        if ~isempty(ud.swapHistory)
            baseColorInd = ud.swapHistory(end,1);
            movingColorInd = ud.swapHistory(end,2);
            %Swap the structure colors
            baseColorTmp = ud.strColorM(baseColorInd,:);
            ud.strColorM(baseColorInd,:) = ud.strColorM(movingColorInd,:);
            ud.strColorM(movingColorInd,:) = baseColorTmp;
            ud.swapHistory(end,:) = [];
            set(hFig,'Userdata',ud);
            structColorGUI('refresh')
        end

    case 'RESET'
        hFig = findobj('Tag','structColorGUI');
        ud = get(hFig,'Userdata');

        ud.strColorM = ud.strColorOldM;

        %Clear Log
        ud.swapHistory = [];

        set(hFig,'Userdata',ud);

        structColorGUI('refresh')

    case 'SAVE'
        hFig = findobj('Tag','structColorGUI');
        ud = get(hFig,'Userdata');

        sliderVis = get(ud.handles.structsSlider,'enable');

        if strcmpi(sliderVis,'off')

            %Store color in planC
            for j =  1:length(ud.assocStructs)
                planC{indexS.structures}(ud.assocStructs(j)).structureColor = ud.strColorM(j,:);
            end
        else
            oldValue = ud.slider.oldValue;
            oldStructNum = getStructDispLen(oldValue,ud.assocStructs);

            for jj = 1: length(oldStructNum)
                planC{indexS.structures}(ud.assocStructs(oldStructNum(jj))).structureColor = ud.strColorM(oldStructNum(jj),:);
            end
        end
        ud.strColorOldM = ud.strColorM;

        stateS.structsChanged = 1;
        CERRRefresh

    case 'EXIT'
        hFig = findobj('Tag','structColorGUI');
        delete(hFig)

end

function numStruct = getStructDispLen(value,structsV)
global planC

indexS = planC{end};

ud = get(findobj('Tag','structColorGUI'),'Userdata');
max = get(ud.handles.structsSlider,'max');

if max == 2
    if value == 2
        numStruct = 1:28;
    elseif value == 1
        numStruct = 29:length(structsV);
    end
elseif max == 3
    if value == 3
        numStruct = 1:28;
    elseif value == 2
        numStruct = 29:56;
    elseif value == 1
        numStruct = 57:length(structsV);
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR