text:
function [f,A,b,sdcone] = optiReadSDPA(filename,dense,print)
%OPTIREADSDPA  Read an SDPA (.dat-s or .dat) file and convert to Matlab Values
%
% [f,A,b,sdcone] = optiReadSDPA(filename) reads the sparse file specified by
% filename and return the semidefinite problem in OPTI format. If you do
% not specify a file a extension, it will default to .dat-s.
%
% [f,A,b,sdcone] = optReadSDPA(filename,dense) allows the user to specify
% whether the file is dense or not. dense = 1 indicates a dense SDPA file,
% and the file extensions will default to .dat if not provided. The default
% is sparse.
%
% [f,A,b,sdcone] = optiReadSDPA(filename,dense,print) allows print out of
% reading progress.
%
%
% *sdcone will be a cell array, where each cell is a sparse matrix of the
% form [C(:) A0(:) A1(:) ... ], containing symmetric matrices.

%   Copyright (C) 2013 Jonathan Currie (I2C2)
%#ok<*AGROW>

if(nargin < 2), dense = 0; end
if(nargin < 3), print = 0; end

% Make sure we have file extension!
if(isempty(strfind(filename,'.')))
    filename = [filename '.dat-s'];
end

%Open File
fid = fopen(filename);
if(fid == -1)
    error('Error opening file %s',filename);
end

%States
sNCON = 1;
sNBLOCKS = 2;
sSIZES = 3;
sRHS = 4;
sDNDATA = 5;
sSPDATA = 6;

%Default Internals
eof = 0;
state = sNCON; cind = 0; bind = 1; rind = 1;
f = []; A = []; b = []; sdcone = [];

if(print)
    ind = strfind(filename,filesep);
    if(isempty(ind))
        file = filename;
    else
        file = filename(ind(end)+1:end);
    end
    if(dense)
        fprintf('Reading Dense SDPA File ''%s''...\n',file);
    else
        fprintf('Reading Sparse SDPA File ''%s''...\n',file);
    end
end
try
    %Main Loop
    while(~eof)
        %Read DAT File String
        str = fgetl(fid);
        %Check for comment, skip if found
        if(str(1) == '*' || str(1) == '"')
            continue;
        end
        %Replace { } ( ) [ ] , ; : with spaces
        str = regexprep(str,{'[',']',',','{','}','(',')'},' ');
        %Decode based on state
        switch(state)
            case sNCON
                M = sscanf(str,'%d',1);
                %Check for errors
                if(M < 1 || M > 1e9), error('The number of constraints is negative or > 1e9'); end
                %Print what we found
                if(print), fprintf('Variables:              %6d\n',M); end
                state = sNBLOCKS;
            case sNBLOCKS
                N = sscanf(str,'%d',1);
                %Check for errors
                if(N < 1 || N > 1e9), error('The number of blocks is negative or > 1e9'); end
                %Print what we found
                if(print), fprintf('Number of Blocks:       %6d\n',N); end
                state = sSIZES;
            case sSIZES
                blocksize = sscanf(str,'%d');
                %Check length = number of blocks
                if(length(blocksize) ~= N), error('The number of blocks does not match the number of entries of block sizes'); end
                if(any(blocksize == 0) || any(blocksize > 1e9)), error('A block size is 0 or > 1e9'); end
                %Print Total Dims
                if(print)
                    fprintf('  Linear:               %6d\n',length(blocksize(blocksize < 0)));
                    fprintf('  Semidefinite:         %6d\n',length(blocksize(blocksize > 0)));
                    fprintf('Total Block Dimension:  %6d\n',sum(abs(blocksize)));
                end
                %Allocate problem data
                sdcone = cell(sum(blocksize > 0),1); j = 1;
                for i = 1:N
                    if(blocksize(i) > 0)
                        if(dense)
                            sdcone{j} = zeros(blocksize(i)^2,M+1);
                        else
                            sdcone{j} = sparse(blocksize(i)^2,M+1); %not the greatest...
                        end
                        j = j + 1;
                    end
                end
                %Create linear part
                ind = blocksize < 0;
                if(any(ind))
                    if(dense)
                        A = zeros(sum(abs(blocksize(ind))),M);
                    else
                        A = sparse(sum(abs(blocksize(ind))),M);
                    end
                    b = zeros(sum(abs(blocksize(ind))),1);
                end
                %Create index vector
                blockindex = zeros(N,1);  lp = 1; sdp = 1;
                for i = 1:N
                    if(blocksize(i) < 0)
                        blockindex(i) = lp; lp  = lp + 1;
                    else
                        blockindex(i) = sdp; sdp = sdp + 1;
                    end
                end
                state = sRHS;
            case sRHS
                f = sscanf(str,'%f'); f = -f;
                %Check length = number of constraints
                if(length(f) ~= M), error('The number of RHS elements does not match the number of constraints'); end
                if(dense)
                    state = sDNDATA;
                else
                    state = sSPDATA;
                end
            %Sparse Data File
            case sSPDATA
                data = sscanf(str,'%d %d %d %d %f');
                if(~isempty(data)) %could be {}
                    %Check for C matrix
                    if(data(1) == 0)
                        %Check block for linear
                        if(blocksize(data(2)) < 0)
                            if(data(3) ~= data(4)), error('A diagonal entry in C is not on the diagonal!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:data(2)-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)));
                            b(stIndex + data(3)) = -data(5);
                        else %semidef
                            sdcone{blockindex(data(2))}(data(3) + (data(4) - 1)*blocksize(data(2)),1) = data(5);
                            %if off-diagonal, place in transposed position (build symmetric matrices)
                            if(data(3) ~= data(4))
                                sdcone{blockindex(data(2))}(data(4) + (data(3) - 1)*blocksize(data(2)),1) = data(5);
                            end
                        end
                    else
                        %Check block for linear
                        if(blocksize(data(2)) < 0)
                            if(data(3) ~= data(4)), error('A diagonal entry in A is not on the diagonal!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:data(2)-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)));
                            A(stIndex + data(3),data(1)) = data(5);
                        else %semidef
                            sdcone{blockindex(data(2))}(data(3) + (data(4) - 1)*blocksize(data(2)),data(1)+1) = -data(5);
                            %if off-diagonal, place in transposed position (build symmetric matrices)
                            if(data(3) ~= data(4))
                                sdcone{blockindex(data(2))}(data(4) + (data(3) - 1)*blocksize(data(2)),data(1)+1) = -data(5);
                            end
                        end
                    end
                end
            %Dense Data File
            case sDNDATA
                data = sscanf(str,'%f');
                if(~isempty(data)) %could be {}
                    %Check for C matrix
                    if(cind == 0)
                        %Check block for linear
                        if(blocksize(bind) < 0)
                            if(length(data) ~= abs(blocksize(bind))), error('Data read in dense entry C (diagonal) is not the correct length!'); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:bind-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)))+1;
                            b(stIndex:stIndex+abs(blocksize(bind))-1) = -data;
                        else
                            if(length(data) == blocksize(bind))
                                sdcone{blockindex(bind)}(rind:blocksize(bind):blocksize(bind)^2) = data;
                                rind = rind + 1;
                            elseif(length(data) == blocksize(bind)^2)
                                %Create matrix and transpose it
                                temp = reshape(data,blocksize(bind),blocksize(bind))';
                                sdcone{blockindex(bind)}(:,1) = temp(:);
                                rind = rind + blocksize(bind);
                            else
                                error('Data read in dense entry C matrix is not the correct length. Expected first row of a matrix, or the entire matrix (all rows and columns, row major).');
                            end
                        end
                    else
                        %Check block for linear
                        if(blocksize(bind) < 0)
                            if(length(data) ~= abs(blocksize(bind))), error('Data read in dense entry A%d (diagonal) is not the correct length!',cind); end
                            %Find start index by determining sizes of previous LP blocks encountered
                            bsizes = blocksize(1:bind-1);
                            stIndex = abs(sum(bsizes(bsizes < 0)))+1;
                            A(stIndex:stIndex+abs(blocksize(bind))-1,cind) = data;
                        else
                            if(length(data) == blocksize(bind))
                                sdcone{blockindex(bind)}(rind + (cind)*blocksize(bind)^2:blocksize(bind):blocksize(bind)^2 + (cind)*blocksize(bind)^2) = -data;
                                rind = rind + 1;
                            elseif(length(data) == blocksize(bind)^2)
                                %Create matrix and transpose it
                                temp = reshape(data,blocksize(bind),blocksize(bind))';
                                sdcone{blockindex(bind)}(:,cind+1) = -temp(:);
                                rind = rind + blocksize(bind);
                            else
                                error('Data read in dense entry A%d matrix is not the correct length. Expected first row of a matrix, or the entire matrix (all rows and columns, row major).',cind);
                            end
                        end
                    end
                    %Increment pointers
                    if(blocksize(bind) < 0 || rind > blocksize(bind))
                        rind = 1;
                        bind = bind + 1;
                    end
                    if(bind > N)
                        bind = 1;
                        cind = cind + 1;
                    end
                end
        end
        eof = feof(fid);
    end
catch ME
    fclose(fid);
    rethrow(ME);
end

%Ensure sparse return args
if(~isempty(sdcone))
    for i = 1:length(sdcone)
        if(~issparse(sdcone{i}))
            sdcone{i} = sparse(sdcone{i});
        end
    end
end
if(~isempty(A))
    if(~issparse(A)), A = sparse(A); end
end
if(print)
    nz = 0;
    if(~isempty(sdcone))
        for i = 1:length(sdcone)
            nz = nz + nnz(sdcone{i});
        end
    end
    fprintf('Total Semidefinite NZs: %6d\n',nz);
    fprintf('Finished Reading SDPA File\n\n');
end
%Close
fclose(fid);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function l_curr = gpu_forward_conv(l_prev, weights)
    tic;
    s = size(l_prev);
    num_prev_row = s(1);
    num_prev_col = s(2);
    num_prev_ch = s(3);

    ss = size(weights);
    num_curr_ch = ss(4);

    %disp('zero-pad the l_prev')
    padded_l_prev = zeros(num_prev_row+2, num_prev_col+2, num_prev_ch);
    padded_l_prev(2:num_prev_row+1, 2:num_prev_col+1, 1:num_prev_ch) = l_prev;
    padded_gpu_l_prev = gpuArray(padded_l_prev);
    %disp('start looping')
    parfor ch = 1:num_curr_ch
        l_curr(:,:,ch) = convn(padded_gpu_l_prev, gpuArray(squeeze(weights(:,:,:,ch))),'valid');
    end
    l_curr = gather(l_curr);
    %l_curr(l_curr<0) = 0; % this is the relu
    toc
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function options = pnopt_optimset( varargin )
% pnopt_optimset : Set options for PNOPT
%
% options = pnopt_optimset( 'param1', val1, 'param2', val2,... ) creates an options
%   structure in which the named parameters have the specified values. Default
%   values are used for parameters that are not specified.
%
% pnopt_optimset with no input arguments and no output arguments displays all
%   parameter names and their default values.
%
%   $Revision: 0.8.0 $  $Date: 2012/12/01 $
%
  % Print out possible values of options.
  if nargin == 0 && nargout == 0
    fprintf('Default parameters for PNOPT:\n');
    fprintf('     curv_param : [  positive scalar |     0.9 ] \n');
    fprintf('          debug : [          logical |       0 ] \n');
    fprintf('     desc_param : [  positive scalar |  0.0001 ] \n');
    fprintf('        display : [ positive integer |      10 ] \n');
    fprintf('      Lbfgs_mem : [ positive integer |      50 ] \n');
    fprintf('  backtrack_mem : [ positive integer |      10 ] \n');
    fprintf('       maxfunEv : [ positive integer |    5000 ] \n');
    fprintf('        maxIter : [ positive integer |     500 ] \n');
    fprintf('         method : [           string |   Lbfgs ] \n');
    fprintf(' subprob_solver : [           string |  sparsa ] \n');
    fprintf('           ftol : [  positive scalar |    1e-9 ] \n');
    fprintf('      optim_tol : [  positive scalar |    1e-6 ] \n');
    fprintf('           xtol : [  positive scalar |    1e-9 ] \n');
    fprintf('\n');
    return;
  end

  Names = [
    'curv_param     '
    'debug          '
    'desc_param     '
    'display        '
    'Lbfgs_mem      '
    'backtrack_mem  '
    'maxfunEv       '
    'maxIter        '
    'method         '
    'sparsa_options '
    'subprob_solver '
    'Tfocs_opts     '
    'ftol           '
    'optim_tol      '
    'xtol           '
    'useMex         '
    ];

  [m,n] = size(Names); %#ok<NASGU>
  names = lower(Names);

  % Combine all leading options structures o1, o2, ... in l1Set(o1,o2,...).
  options = [];
  for j = 1:m
    eval(['options.' Names(j,:) '= [];']);
  end
  i = 1;
  while i <= nargin
    arg = varargin{i};
    if ischar(arg), break; end
    if ~isempty(arg)                      % [] is a valid options argument
      if ~isa(arg,'struct')
        error(sprintf(['Expected argument %d to be a string parameter name ' ...
          'or an options structure\ncreated with SetPNoptOptions.'], i)); %#ok<SPERR>
      end
      for j = 1:m
        if any(strcmp(fieldnames(arg),deblank(Names(j,:))))
          eval(['val = arg.' Names(j,:) ';']);
        else
          val = [];
        end
        if ~isempty(val)
          eval(['options.' Names(j,:) '= val;']);
        end
      end
    end
    i = i + 1;
  end

  % A finite state machine to parse name-value pairs.
  if rem(nargin-i+1,2) ~= 0
    error('Arguments must occur in name-value pairs.');
  end
  expectval = 0;                          % start expecting a name, not a value
  while i <= nargin
    arg = varargin{i};

    if ~expectval
      if ~ischar(arg)
         error('Expected argument %d to be a string parameter name.', i);
      end

      lowArg = lower(arg);
      j = strmatch(lowArg,names);
      if isempty(j)                       % if no matches
        error('Unrecognized parameter name ''%s''.', arg);
      elseif length(j) > 1                % if more than one match
        % Check for any exact matches (in case any names are subsets of others)
        k = strmatch(lowArg,names,'exact');
        if length(k) == 1
          j = k;
        else
          msg = sprintf('Ambiguous parameter name ''%s'' ', arg);
          msg = [msg '(' deblank(Names(j(1),:))]; %#ok<AGROW>
          for k = j(2:length(j))'
            msg = [msg ', ' deblank(Names(k,:))]; %#ok<AGROW>
          end
          error('%s).', msg);
        end
      end
      expectval = 1;                      % we expect a value next

    else
      eval(['options.' Names(j,:) '= arg;']);
      expectval = 0;

    end
    i = i + 1;
  end

  if expectval
    error('Expected value for parameter ''%s''.', arg);
  end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function status = test_modal_weighting(modus)
%TEST_MODAL_WEIGHTING evaluates the modal weighting implementation
%
%   Usage: status = test_modal_weighting(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual
%
%   Output parameters:
%       status  - true or false
%
%   TEST_MODAL_WEIGHTING(modus) tests the implementation of the modal weighting
%   window for different local sound field synthesis techniques.

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

status = false;

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Configuration ==================================================
window_type = {...
    'rect'; ...
    'max-rE'; ...
    'kaiser'; ...
    'tukey'; ...
    };
% For plotting
conf = SFS_config;
conf.dimension = '2.5D';
conf.secondary_sources.number = 16;
conf.secondary_sources.geometry = 'circular';
conf.secondary_sources.size = 3;
conf.resolution = 300;
conf.plot.usedb = true;
conf.t0 = 'source';
X = [-2 2];
Y = [-2 2];
Z = 0;
src = 'pw';

%% ===== Computation and Plotting ========================================
if modus % graphical mode
    for window = window_type'
        conf.modal_window = window{:};
        update_title = false;
        base_title = ['Window: ' window{:}];
        switch conf.modal_window
        case {'rect', 'max-rE'}
            parameter = 1;
        otherwise
            parameter = [0 0.5 1];
            update_title = true;
        end
        for pp = 1:length(parameter)
            conf.modal_window_parameter = parameter(pp);
            sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
            if update_title
                figure_title = [base_title ...
                                ', Parameter: ' ...
                                num2str(parameter(pp))];
            else
                figure_title = base_title;
            end
            title(figure_title);
        end
    end
else % numerical mode
    for window = window_type'
        conf.modal_window = window{:};
        conf.modal_window_parameter = 0.5;
        for order = [1 2]
            for ndtft = [2*order+1 3*order+1]
                [win,Win,Phi] = ...
                    modal_weighting(order,ndtft,conf);
                [win_ref,Win_ref,Phi_ref] = ...
                    ref_modal_weighting(order,ndtft,conf);
                if any([abs(win-win_ref) ...
                        abs(Win-Win_ref) ...
                        abs(Phi-Phi_ref)])>1e14
                    return
                end
            end
        end
    end
end

status = true;

end

function [win,Win,Phi] = ref_modal_weighting(order,ndtft,conf)
    % Returns reference values for some orders as calculated by the Matlab
    % implementation of commit ...
    switch conf.modal_window
    case 'rect'
        if order==1 && ndtft==2*order+1
            win = [1 1];
            Win = [1 0 0];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 1];
            Win = [0.75 0.25 -0.25 0.25];
            Phi = [0 1.570796326794897 3.141592653589793 4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 1 1];
            Win = [1 0 0 0 0];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 1 1];
            Win = [0.714285714285714 ...
                   0.257419676543548 ...
                  -0.178139943388210 ...
                   0.063577409701804 ...
                   0.063577409701804 ...
                  -0.178139943388210 ...
                   0.257419676543548];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case 'max-rE'
        if order==1 && ndtft==2*order+1
            win = [1 0.707106781186548];
            Win = [0.804737854124365 0.097631072937817 0.097631072937817];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.707106781186548];
            Win = [0.603553390593274 ...
                   0.250000000000000 ...
                  -0.103553390593274 ...
                   0.250000000000000];
            Phi = [0 ...
                   1.570796326794897 ...
                   3.141592653589793 ...
                   4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 0.866025403784439 0.500000000000000];
            Win = [0.746410161513775 ...
                   0.145243228056937 ...
                  -0.018448308813825 ...
                  -0.018448308813825 ...
                   0.145243228056937];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 0.866025403784439 0.500000000000000];
            Win = [0.533150115366983 ...
                   0.265342154409152 ...
                  -0.040912347323205 ...
                   0.008995135230562 ...
                   0.008995135230562 ...
                  -0.040912347323205 ...
                   0.265342154409152];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case {'kaiser', 'kaiser-bessel'}
        % conf.modal_window_parameter == 0.5
        if order==1 && ndtft==2*order+1
            win = [1 0.581816879161903];
            Win = [0.721211252774602 0.139394373612699 0.139394373612699];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.581816879161903];
            Win = [0.540908439580952 ...
                   0.250000000000000 ...
                  -0.040908439580952 ...
                   0.250000000000000];
            Phi = [0 ...
                   1.570796326794897 ...
                   3.141592653589793 ...
                   4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 0.883766863470753 0.581816879161903];
            Win = [0.786233497053062 ...
                   0.120959694808693 ...
                  -0.014076443335224 ...
                  -0.014076443335224 ...
                   0.120959694808693];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 0.883766863470753 0.581816879161903];
            Win = [0.561595355037902 ...
                   0.263300911786297 ...
                  -0.063101577944049 ...
                   0.019002988638802 ...
                   0.019002988638802 ...
                  -0.063101577944049 ...
                   0.263300911786297];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case 'tukey'
        % conf.modal_window_parameter == 0.5
        if order==1 && ndtft==2*order+1
            win = [1 0.75];
            Win = [0.833333333333333 0.083333333333333 0.083333333333333];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.75];
            Win = [0.625 0.250 -0.125 0.250];
            Phi = [0 1.570796326794897 3.141592653589793 4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 1 0.5];
            Win = [0.800000000000000 ...
                   0.161803398874989 ...
                  -0.061803398874989 ...
                  -0.061803398874989 ...
                   0.161803398874989];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 1 0.5];
            Win = [0.571428571428571 ...
                   0.289208381394450 ...
                  -0.049430105116435 ...
                  -0.025492561992301 ...
                  -0.025492561992301 ...
                  -0.049430105116435 ...
                   0.289208381394450];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
function P = khatrirao(varargin)
%KHATRIRAO Khatri-Rao product of matrices.
%
%   KHATRIRAO(A,B) computes the Khatri-Rao product of matrices A and
%   B that have the same number of columns.  The result is the
%   column-wise Kronecker product
%   [KRON(A(:,1),B(:,1)) ... KRON(A(:,n),B(:,n))]
%
%   KHATRIRAO(A1,A2,...) computes the Khatri-Rao product of
%   multiple matrices that have the same number of columns.
%
%   KHATRIRAO(C) computes the Khatri-Rao product of
%   the matrices in cell array C.
%
%   KHATRIRAO(...,'r') computes the Khatri-Rao product in reverse
%   order.
%
%   Examples
%   A = rand(5,2); B = rand(3,2); C = rand(2,2);
%   khatrirao(A,B) %<-- Khatri-Rao of A and B
%   khatrirao(B,A,'r') %<-- same thing as above
%   khatrirao({C,B,A}) %<-- passing a cell array
%   khatrirao({A,B,C},'r') %<-- same as above
%
%   See also TENSOR, 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

%% Error checking on input and set matrix order
% Note that this next if/else check forces A to be a cell array.
if ischar(varargin{end}) && varargin{end} == 'r'
    if nargin == 2 && iscell(varargin{1})
        % Input is a single cell array
        A = varargin{1};
    else
        % Input is a sequence of matrices
        A = {varargin{1:end-1}};
    end
    matorder = length(A):-1:1;
else
    if nargin == 1 && iscell(varargin{1})
        % Input is a single cell array
        A = varargin{1};
    else
        % Input is a sequence of matrices
        A = varargin;
    end
    matorder = 1:length(A);
end

%% Error check on matrices and compute number of rows in result

% N = number of columns (must be the same for every input)
N = size(A{1},2);

% After loop, M = number of rows in the result
M = 1;

for i = matorder
    if ndims(A) ~= 2
        error('Each argument must be a matrix');
    end
    if (N ~= size(A{i},2))
        error('All matrices must have the same number of columns.')
    end
    M = M * size(A{i},1);
end

%% Computation

% Preallocate
P = zeros(M,N);

% Loop through all the columns
for n = 1:N
    % Loop through all the matrices
    ab = A{matorder(1)}(:,n);
    for i = matorder(2:end)
       % Compute outer product of nth columns
       ab = A{i}(:,n) * ab(:).';
    end
    % Fill nth column of P with reshaped result
    P(:,n) = ab(:);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%% PTspec2DUIcontrol - ui controls for spectral analyses plots

% ----------------------------------------------------------------------------------
% "THE BEER-WARE LICENSE" (Revision 42):
% <brian.white@queensu.ca> wrote this file. As long as you retain this notice you
% can do whatever you want with this stuff. If we meet some day, and you think
% this stuff is worth it, you can buy me a beer in return. -Brian White
% ----------------------------------------------------------------------------------

if ~isempty(fnameMaster)

%%% tooltips
TooltipString_specRun=['Run current spectral configuration'];
TooltipString_cmap=['Choose from a selection of colormaps'];
TooltipString_smooth=['Choose amount of smoothing'];
TooltipString_user=['Choose the variable you wish to plot'];
TooltipString_sub100=['Zoom data to show sub 100Hz details',...
    newline, 'Typically used to see propwash or mid-throttle vibration in e.g. Gyro/Pterm/PIDerror'];

%%%
PTcolormap;
% define
smat=[];%string
ampmat=[];%spec matrix
amp2d2=[];%spec 2d
freq2d2=[];% freq

clear posInfo.Spec2Pos
cols=[0.05 0.48];
rows=[0.69 0.395 0.1];
k=0;
for c=1 : size(cols,2)
    for r=1 : size(rows,2)
        k=k+1;
        if c == 1
            posInfo.Spec2Pos(k,:)=[cols(c) rows(r) 0.38 0.25];
        else
            posInfo.Spec2Pos(k,:)=[cols(c) rows(r) 0.38 0.25];
        end
    end
end

vPosSpec2d = .037;

posInfo.fileListWindowSpec=[.898 .7+vPosSpec2d .088 .20];
posInfo.TermListWindowSpec=[.898 .55+vPosSpec2d .088 .14];

posInfo.computeSpec=            [.896 .52+vPosSpec2d .0455 .026];
posInfo.resetSpec=              [.942 .52+vPosSpec2d .0455 .026];
posInfo.spectrogramButton2=     [.896 .495+vPosSpec2d .0915 .026];
posInfo.spectrogramButton3=    [.896 .47+vPosSpec2d .0915 .026];
posInfo.saveFig2=               [.896 .445+vPosSpec2d .0455 .026]; % .896 .495 .092 .026
posInfo.saveSettings2=          [.942 .445+vPosSpec2d .0455 .026];

posInfo.smooth_select =         [.895 .405+vPosSpec2d .0955 .04];
posInfo.Delay =                 [.895 .379+vPosSpec2d .0955 .04];

posInfo.plotRspec =                 [.90 .37+vPosSpec2d .03 .025];
posInfo.plotPspec =                 [.925 .37+vPosSpec2d .03 .025];
posInfo.plotYspec =                 [.95 .37+vPosSpec2d .03 .025];

posInfo.checkboxPSD =            [.90 .35+vPosSpec2d .03 .02];
posInfo.RPYcomboSpec =           [.932 .35+vPosSpec2d .055 .02];

posInfo.climMax1_text = [.894 .322+vPosSpec2d .026 .022];
posInfo.climMax1_input = [.918 .322+vPosSpec2d .022 .022];
posInfo.climMax2_text = [.943 .322+vPosSpec2d .026 .022];
posInfo.climMax2_input = [.968 .322+vPosSpec2d .022 .022];

climScale1=[0 ; -50 ];
climScale2=[0.5 ; 20];

PTspecfig2=figure(3);
set(PTspecfig2, 'units','normalized','outerposition',[.1 .1 .75 .8])
PTspecfig2.NumberTitle='off';
PTspecfig2.Name= ['PIDtoolbox (' PtbVersion ') - Spectral Analyzer'];
PTspecfig2.InvertHardcopy='off';
set(PTspecfig2,'color',bgcolor);

dcm_obj2 = datacursormode(PTspecfig2);
set(dcm_obj2,'UpdateFcn',@PTdatatip);

specCrtlpanel = uipanel('Title','select files (max 10)','FontSize',fontsz,...
              'BackgroundColor',[.95 .95 .95],...
              'Position',[.89 .31+vPosSpec2d .105 .61]);

guiHandlesSpec2.computeSpec = uicontrol(PTspecfig2,'string','Run','fontsize',fontsz,'TooltipString', [TooltipString_specRun],'units','normalized','outerposition',[posInfo.computeSpec],...
    'callback','PTplotSpec2D;');
guiHandlesSpec2.computeSpec.ForegroundColor=[colRun];

guiHandlesSpec2.resetSpec = uicontrol(PTspecfig2,'string','Reset','fontsize',fontsz,'TooltipString', ['Reset Spectral Tool'],'units','normalized','outerposition',[posInfo.resetSpec],...
    'callback',' for k = 1 : 6, delete(subplot(''position'',posInfo.Spec2Pos(k,:))), end; set(PTspecfig2, ''pointer'', ''arrow'');');
guiHandlesSpec2.resetSpec.ForegroundColor=[cautionCol];

guiHandlesSpec2.saveFig2 = uicontrol(PTspecfig2,'string','Save Fig','fontsize',fontsz,'TooltipString',[TooltipString_saveFig],'units','normalized','ForegroundColor',[saveCol],'outerposition',[posInfo.saveFig2],...
    'callback','guiHandlesSpec2.saveFig2.FontWeight=''bold'';PTsaveFig;guiHandlesSpec2.saveFig2.FontWeight=''normal'';');

guiHandlesSpec2.saveSettings2 = uicontrol(PTspecfig2,'string','Save Settings','fontsize',fontsz, 'TooltipString',['Save current settings to PTB defaults' ], 'units','normalized','outerposition',[posInfo.saveSettings2],...
    'callback','guiHandlesSpec2.saveSettings2.FontWeight=''bold'';PTsaveSettings; guiHandlesSpec2.saveSettings2.FontWeight=''normal'';');
guiHandlesSpec2.saveSettings2.ForegroundColor=[saveCol];

% create string list for SpecSelect
sA={'Gyro','Gyro prefilt','Dterm','Dterm prefilt','Pterm','PID error','Set point','PIDsum'};

guiHandlesSpec2.SpecList = uicontrol(PTspecfig2,'Style','listbox','string',[sA],'max',3,'min',1, 'fontsize',fontsz, 'TooltipString',[TooltipString_user],'units','normalized','outerposition', [posInfo.TermListWindowSpec], 'callback', 'if length(guiHandlesSpec2.SpecList.Value) > 2, guiHandlesSpec2.SpecList.Value=1; end;');
guiHandlesSpec2.SpecList.Value =[1 2];

guiHandlesSpec2.FileSelect = uicontrol(PTspecfig2,'Style','listbox','string',[fnameMaster],'max', 10, 'min', 1, 'fontsize',fontsz,'TooltipString',[TooltipString_user],'units','normalized','outerposition', [posInfo.fileListWindowSpec], 'callback', 'if length(guiHandlesSpec2.FileSelect.Value) > 10, guiHandlesSpec2.FileSelect.Value=1; end;');

guiHandlesSpec2.smoothFactor_select = uicontrol(PTspecfig2,'style','popupmenu','string',{'smoothing low' 'smoothing low-med' 'smoothing medium' 'smoothing med-high' 'smoothing high'},'fontsize',fontsz,'TooltipString', [TooltipString_smooth], 'units','normalized','outerposition',[posInfo.smooth_select],...
     'callback','@selection2;PTplotSpec2D;');

guiHandlesSpec2.spectrogramButton2 = uicontrol(PTspecfig2,'string','Freq x Throttle','fontsize',fontsz,'TooltipString', ['Opens Freq x Throttle Spectrogram in New Window'], 'units','normalized','outerposition',[posInfo.spectrogramButton2],...
    'callback','PTspecUIcontrol;');
guiHandlesSpec2.spectrogramButton2.ForegroundColor=[colorA];

 guiHandlesSpec2.spectrogramButton3 = uicontrol(PTspecfig2,'string','Freq x Time','fontsize',fontsz,'TooltipString', ['Opens Freq x Time Spectrogram in New Window'], 'units','normalized','outerposition',[posInfo.spectrogramButton3],...
     'callback','PTfreqTimeUIcontrol;');
 guiHandlesSpec2.spectrogramButton3.ForegroundColor=[colorB];

 guiHandlesSpec2.Delay = uicontrol(PTspecfig2,'style','popupmenu','string',{'filter delay', 'SP-gyro delay', 'SP smoothing delay', 'phase shift'},'fontsize',fontsz,'TooltipString', ['Select which Delay Display'], 'units','normalized','outerposition',[posInfo.Delay],...
     'callback','PTplotSpec2D;');

guiHandlesSpec2.plotR =uicontrol(PTspecfig2,'Style','checkbox','String','R','fontsize',fontsz,'TooltipString', ['Plot Roll '],...
    'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.plotRspec], 'callback', 'for k = 1 : 6, delete(subplot(''position'',posInfo.Spec2Pos(k,:))), end; set(PTspecfig2, ''pointer'', ''arrow'');');

guiHandlesSpec2.plotP =uicontrol(PTspecfig2,'Style','checkbox','String','P','fontsize',fontsz,'TooltipString', ['Plot Pitch '],...
    'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.plotPspec], 'callback', 'for k = 1 : 6, delete(subplot(''position'',posInfo.Spec2Pos(k,:))), end; set(PTspecfig2, ''pointer'', ''arrow'');');

guiHandlesSpec2.plotY =uicontrol(PTspecfig2,'Style','checkbox','String','Y','fontsize',fontsz,'TooltipString', ['Plot Yaw '],...
    'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.plotYspec], 'callback', 'for k = 1 : 6, delete(subplot(''position'',posInfo.Spec2Pos(k,:))), end; set(PTspecfig2, ''pointer'', ''arrow'');');

guiHandlesSpec2.checkboxPSD =uicontrol(PTspecfig2,'Style','checkbox','String','PSD','fontsize',fontsz,'TooltipString', ['Power Spectral Density'],...
    'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.checkboxPSD],'callback', 'PTplotSpec2D;');
guiHandlesSpec2.checkboxPSD.Value = 1;

guiHandlesSpec2.RPYcomboSpec =uicontrol(PTspecfig2,'Style','checkbox','String','Single Panel','fontsize',fontsz,'TooltipString', ['Plot RPY in same panel '],...
    'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.RPYcomboSpec],'callback', 'PTplotSpec2D;');

guiHandlesSpec2.climMax1_text = uicontrol(PTspecfig2,'style','text','string','Y min','fontsize',fontsz,'TooltipString',['Y min'],'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.climMax1_text]);
guiHandlesSpec2.climMax1_input = uicontrol(PTspecfig2,'style','edit','string',[num2str(climScale1(guiHandlesSpec2.checkboxPSD.Value+1, 1))],'fontsize',fontsz,'TooltipString',['Y min'],'units','normalized','outerposition',[posInfo.climMax1_input],...
     'callback','@textinput_call2; climScale1(guiHandlesSpec2.checkboxPSD.Value+1, 1)=str2num(guiHandlesSpec2.climMax1_input.String);PTplotSpec2D;');

 guiHandlesSpec2.climMax2_text = uicontrol(PTspecfig2,'style','text','string','Y max','fontsize',fontsz,'TooltipString',['Y max'],'units','normalized','BackgroundColor',bgcolor,'outerposition',[posInfo.climMax2_text]);
guiHandlesSpec2.climMax2_input = uicontrol(PTspecfig2,'style','edit','string',[num2str(climScale2(guiHandlesSpec2.checkboxPSD.Value+1, 1))],'fontsize',fontsz,'TooltipString',['Y max'],'units','normalized','outerposition',[posInfo.climMax2_input],...
     'callback','@textinput_call2; climScale2(guiHandlesSpec2.checkboxPSD.Value+1, 1)=str2num(guiHandlesSpec2.climMax2_input.String);PTplotSpec2D;');

try guiHandlesSpec2.SpecList.Value = [defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-term1'))) defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-term2')))]; catch, guiHandlesSpec2.SpecList.Value = [1 2], end
try guiHandlesSpec2.smoothFactor_select.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-smoothing'))), catch, guiHandlesSpec2.smoothFactor_select.Value = 3, end
try guiHandlesSpec2.Delay.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-delay'))), catch, guiHandlesSpec2.Delay.Value = 1, end
try guiHandlesSpec2.plotR.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-plotR'))), catch, guiHandlesSpec2.plotR.Value = 1, end
try guiHandlesSpec2.plotP.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-plotP'))), catch, guiHandlesSpec2.plotP.Value = 1, end
try guiHandlesSpec2.plotY.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-plotY'))), catch, guiHandlesSpec2.plotY.Value = 1, end
try guiHandlesSpec2.RPYcomboSpec.Value = defaults.Values(find(strcmp(defaults.Parameters, 'spec2D-SinglePanel'))), catch, guiHandlesSpec2.RPYcomboSpec.Value = 0, end

FilterDelayDterm={};
SPGyroDelay=[];
Debug01={};
Debug02={};
gyro_phase_shift_deg=[];
dterm_phase_shift_deg=[];
for k = 1 : Nfiles
    Fs=1000/A_lograte(k);% yields more consistent results (mode(diff(tta)));
    maxlag=int8(round(30000/Fs)); %~30ms delay

    clear d pg g1 g1 s1 g2 s2  g3 s3
    try
        pg = smooth(T{k}.debug_0_(tIND{k}),50);
    catch
        pg = 0;
    end
    g1 = smooth(T{k}.gyroADC_0_(tIND{k}),50);
    s1 = smooth(T{k}.setpoint_0_(tIND{k}),50);

    g2 = smooth(T{k}.gyroADC_1_(tIND{k}),50);
    s2 = smooth(T{k}.setpoint_1_(tIND{k}),50);

    g3 = smooth(T{k}.gyroADC_2_(tIND{k}),50);
    s3 = smooth(T{k}.setpoint_2_(tIND{k}),50);

    [c,lags] = xcorr(g1,pg,maxlag);
    d = lags(find(c==max(c),1));
    d = d * (Fs / 1000);
    if d<.1,  Debug01{k} = ' '; else Debug01{k} = num2str(d);end

    [c,lags] = xcorr(s1,pg,maxlag);
    d = lags(find(c==max(c),1));
    d = d * (Fs / 1000);
    if d<.1,  Debug02{k} = ' '; else Debug02{k} = num2str(d);end

    [c,lags] = xcorr(g1,s1,maxlag);
    d = lags(find(c==max(c),1));
    d = d * (Fs / 1000);
    if d<.1, SPGyroDelay(k,1) = 0, else, SPGyroDelay(k,1) = d; end

    [c,lags] = xcorr(g2,s2,maxlag);
    d = lags(find(c==max(c),1));
    d = d * (Fs / 1000);
    if d<.1, SPGyroDelay(k,2) = 0, else, SPGyroDelay(k,2) = d; end

    [c,lags] = xcorr(g3,s3,maxlag);
    d = lags(find(c==max(c),1));
    d = d * (Fs / 1000);
    if d<.1, SPGyroDelay(k,3) = 0, else, SPGyroDelay(k,3) = d; end

    clear d d1 d2
    d1 = smooth(T{k}.axisDpf_0_(tIND{k}),50);
    d2 = smooth(T{k}.axisD_0_(tIND{k}),50);
    [c,lags] = xcorr(d2,d1,maxlag);
    d = lags(find(c==max(c)));
    d=d * (Fs / 1000);
    if d<.1, FilterDelayDterm{k} = ' '; else FilterDelayDterm{k} = num2str(d); end

    if ~isempty(str2num(Debug01{k})) && ~isempty(SPGyroDelay(k,1))
        [gyro_phase_shift_deg(k,1)] = round(PTphaseShiftDeg(str2num(Debug01{k}), 1000/(SPGyroDelay(k,1)) ));
    end
    if ~isempty(str2num(FilterDelayDterm{k})) && ~isempty(SPGyroDelay(k,1))
        [dterm_phase_shift_deg(k,1)] = round(PTphaseShiftDeg(str2num(FilterDelayDterm{k}), 1000/(SPGyroDelay(k,1)) ));
    end
end

else
     warndlg('Please select file(s)');
end

% functions
function selection2(src,event)
    val = c.Value;
    str = c.String;
    str{val};
end

function getList2(hObj,event)
v=get(hObj,'value')
end

function textinput_call2(src,eventdata)
str=get(src,'String');
    if isempty(str2num(str))
        set(src,'string','0');
        warndlg('Input must be numerical');
    end
end

repo: PIDtoolbox


text:
function ray = createRay(varargin)
%CREATERAY Create a ray (half-line), from various inputs.
%
%   RAY = createRay(POINT, ANGLE)
%   POINT is a N*2 array giving starting point of the ray, and ANGLE is the
%   orientation of the ray.
%
%   RAY = createRay(X0, Y0, ANGLE)
%   Specify ray origin with 2 input arguments.
%
%   RAY = createRay(P1, P2)
%   Create a ray starting from point P1 and going in the direction of point
%   P2.
%
%   Ray is represented in a parametric form: [x0 y0 dx dy]
%   x = x0 + t*dx
%   y = y0 + t*dy;
%   for all t>0
%
%   Example
%   origin  = [3 4];
%   theta   = pi/6;
%   ray = createRay(origin, theta);
%   figure(1); clf; hold on;
%   axis([0 10 0 10]);
%   drawRay(ray);
%
%   See also
%   rays2d, createLine, points2d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2007-10-18
% Copyright 2007-2022 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas

if length(varargin)==2
    p0 = varargin{1};
    arg = varargin{2};
    if size(arg, 2)==1
        % second input is the ray angle
        ray = [p0 cos(arg) sin(arg)];
    else
        % second input is another point
        ray = [p0 arg-p0];
    end

elseif length(varargin)==3
    x = varargin{1};
    y = varargin{2};
    theta = varargin{3};
    ray = [x y cos(theta) sin(theta)];

else
    error('Wrong number of arguments in ''createRay'' ');
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%
% In order to obtain a simple and effective source code for
% Fingerprint Recognition System please visit
% http://utenti.lycos.it/matlab/beginner.htm
%
% There you will be able to make a small donation
% ( just 2 Euros, less than 2.8 US Dollars ) using
% the secure connection to PayPal.
%
% The source code implements FingerCode-based fingerprint
% verification.
%
% References:
%
%   A. K. Jain, S. Prabhakar, and S. Pankanti, "A Filterbank-based Representation for
%   Classification and Matching of Fingerprints", International Joint Conference on
%   Neural Networks (IJCNN), pp. 3284-3285, Washington DC, July 10-16, 1999.
%   http://www.cse.msu.edu/~prabhaka/publications.html
%
%   "Fingerprint Classification and Matching Using a Filterbank", Salil Prabhakar
%   A DISSERTATION Submitted to Michigan State University in partial fulfillment
%   of the requirements for the degree of DOCTOR OF PHILOSOPHY, Computer
%   Science & Engineering, 2001
%   http://biometrics.cse.msu.edu/SalilThesis.pdf
%
% For any question please email me luigi.rosa@tiscali.it
%
% Luigi ROSA
% Via Centrale 35
% 67042 Civita Di Bagno
% L'Aquila - ITALY
% mobile +39 3207214179
% email luigi.rosa@tiscali.it
%
%
%

author: Sable

repo: mcbench-benchmarks


text:
% created: Zoya Bylinskii, Aug 2014

% This finds the KL-divergence between two different saliency maps when
% viewed as distributions: it is a non-symmetric measure of the information
% lost when saliencyMap is used to estimate fixationMap.

function score = KLdiv(saliencyMap, fixationMap)
% saliencyMap is the saliency map
% fixationMap is the human fixation map

map1 = im2double(imresize(saliencyMap, size(fixationMap)));
map2 = im2double(fixationMap);

% make sure map1 and map2 sum to 1
if any(map1(:))
    map1 = map1/sum(map1(:));
end

if any(map2(:))
    map2 = map2/sum(map2(:));
end

% compute KL-divergence
score = sum(sum(map2 .* log(eps + map2./(map1+eps))));

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function value = smach ( job )

%*****************************************************************************80
%
%% SMACH computes machine parameters of floating point arithmetic.
%
%  Discussion:
%
%    This routine is for testing only.  It is not required by LINPACK.
%
%    If there is trouble with the automatic computation of these quantities,
%    they can be set by direct assignment statements.
%    We assume the computer has
%
%      B = base of arithmetic
%      T = number of base B digits
%      L = smallest possible exponent
%      u = largest possible exponent
%
%    then
%
%      EPS = B**(1-T)
%      TINY = 100.0D+00 *B**(-L+T)
%      HUGE = 0.01D+00 *B**(U-T)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer JOB:
%    1: requests EPSILON;
%    2: requests TINY;
%    3: requests HUGE.
%
%    Output, real VALUE, the requested value.
%
  s = 1.0;

  while ( 1 )

    tiny = s;
    s = s / 2.0;

    if ( s * 1.0 == 0.0 | s == 0.0 )
      break
    end

  end

  tiny = ( tiny / eps ) * 100.0;
  huge = 1.0 / tiny;

  if ( job == 1 )
    value = eps;
  elseif ( job == 2 )
    value = tiny;
  elseif ( job == 3 )
    value = huge;
  else
    xerbla ( 'SMACH', 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = daub6_scale ( n, x )

%*****************************************************************************80
%
%% DAUB6_SCALE recursively evaluates the DAUB6 scaling function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the recursion level.
%
%    Input, real X, the point at which the function is to be evaluated.
%
%    Output, real Y, the estimated value of the function.
%
  c = [  0.3326705529500826E+00; ...
         0.8068915093110925E+00; ...
         0.4598775021184915E+00; ...
       - 0.1350110200102545E+00; ...
       - 0.08544127388202666E+00; ...
         0.03522629188570953E+00 ];

  c = c * sqrt ( 2.0 );

  if ( 0 < n )
    y = c(1) * daub6_scale ( n - 1, 2 * x     ) ...
      + c(2) * daub6_scale ( n - 1, 2 * x - 1 ) ...
      + c(3) * daub6_scale ( n - 1, 2 * x - 2 ) ...
      + c(4) * daub6_scale ( n - 1, 2 * x - 3 ) ...
      + c(5) * daub6_scale ( n - 1, 2 * x - 4 ) ...
      + c(6) * daub6_scale ( n - 1, 2 * x - 5 );
  else
    y = ( 0 <= x & x < 1 );
  end

  return
end


text:
% SB1_EXAMPLECLASSIFY   Example of Sparse Bayes Classification
%
% SB1_EXAMPLECLASSIFY(N,KERNEL,WIDTH,MAXITS)
%
% INPUT ARGUMENTS:
%
%       N       Number of training points (up to 250)
%       KERNEL  Kernel function to use (see SB1_KERNELFUNCTION)
%       WIDTH   Kernel length scale parameter
%       MAXITS  Maximum number of iterations to run for
%
% NOTES: The data is taken from the book "Pattern Recognition for
%        Neural Networks" (Ripley). Running SB1_EXAMPLECLASSIFY with
%        no arguments will result in sensible defaults.
%
%
% Copyright 2009 :: Michael E. Tipping
%
% This file is part of the SPARSEBAYES baseline implementation (V1.10)
%
% Contact the author: m a i l [at] m i k e t i p p i n g . c o m
%
function SB1_ExampleClassify(N,kernel_,width,maxIts)

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

% Set verbosity of output (0 to 4)
setEnvironment('Diagnostic','verbosity',3);
% Set file ID to write to (1 = stdout)
setEnvironment('Diagnostic','fid',1);

%
% Set default values for data and model
%
useBias	= true;
%
rand('state',1)
if nargin==0
  % Some acceptable defaults
  %
  N		= 100;
  kernel_	= 'gauss';
  width		= 0.5;
  maxIts	= 1000;
end
monIts		= round(maxIts/10);
N		= min([250 N]); % training set has fixed size of 250
Nt		= 1000;
%
% Load Ripley's synthetic training data (see reference in doc)
%
load synth.tr
synth	= synth(randperm(size(synth,1)),:);
X	= synth(1:N,1:2);
t	= synth(1:N,3);

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

COL_data1	= 'k';
COL_data2	= 0.75*[0 1 0];
COL_boundary50	= 'r';
COL_boundary75	= 0.5*ones(1,3);
COL_rv		= 'r';

%
% Plot the training data
%
figure(1)
whitebg(1,'w')
clf
h_c1	= plot(X(t==0,1),X(t==0,2),'.','MarkerSize',18,'Color',COL_data1);
hold on
h_c2	= plot(X(t==1,1),X(t==1,2),'.','MarkerSize',18,'Color',COL_data2);
box	= 1.1*[min(X(:,1)) max(X(:,1)) min(X(:,2)) max(X(:,2))];
axis(box)
set(gca,'FontSize',12)
drawnow

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

%
% Set up initial hyperparameters - precise settings should not be critical
%
initAlpha	= (1/N)^2;
% Set beta to zero for classification
initBeta	= 0;
%
% "Train" a sparse Bayes kernel-based model (relevance vector machine)
%
[weights, used, bias, marginal, alpha, beta, gamma] = ...
    SB1_RVM(X,t,initAlpha,initBeta,kernel_,width,useBias,maxIts,monIts);

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

%
% Load Ripley's test set
%
load synth.te
synth	= synth(randperm(size(synth,1)),:);
Xtest	= synth(1:Nt,1:2);
ttest	= synth(1:Nt,3);
%
% Compute RVM over test data and calculate error
%
PHI	= SB1_KernelFunction(Xtest,X(used,:),kernel_,width);
y_rvm	= PHI*weights + bias;
errs	= sum(y_rvm(ttest==0)>0) + sum(y_rvm(ttest==1)<=0);
SB1_Diagnostic(1,'RVM CLASSIFICATION test error: %.2f%%\n', errs/Nt*100)

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

%
% Visualise the results over a grid
%
gsteps		= 50;
range1		= box(1):(box(2)-box(1))/(gsteps-1):box(2);
range2		= box(3):(box(4)-box(3))/(gsteps-1):box(4);
[grid1 grid2]	= meshgrid(range1,range2);
Xgrid		= [grid1(:) grid2(:)];
size(Xgrid)
%
% Evaluate RVM
%
PHI		= SB1_KernelFunction(Xgrid,X(used,:),kernel_,width);
y_grid		= PHI*weights + bias;

% apply sigmoid for probabilities
p_grid		= 1./(1+exp(-y_grid));

%
% Show decision boundary (p=0.5) and illustrate p=0.25 and 0.75
%
[c,h05]		= ...
    contour(range1,range2,reshape(p_grid,size(grid1)),[0.5],'-');
[c,h075]	= ...
    contour(range1,range2,reshape(p_grid,size(grid1)),[0.25 0.75],'--');
set(h05, 'Color',COL_boundary50,'LineWidth',3);
set(h075,'Color',COL_boundary75,'LineWidth',2);
%
% Show relevance vectors
%
h_rv	= plot(X(used,1),X(used,2),'o','LineWidth',2,'MarkerSize',10,...
	       'Color',COL_rv);
%
legend([h_c1 h_c2 h05 h075(1) h_rv],...
       'Class 1','Class 2','Decision boundary','p=0.25/0.75','RVs',...
       'Location','NorthWest')
%
hold off
title('RVM Classification of Ripley''s synthetic data','FontSize',14)

author: epfl-lasa


text:
function obj = updateColorbar(obj,colorbarIndex)

    % title: ...[DONE]
    % titleside: ...[DONE]
    % titlefont: ...[DONE]
    % thickness: ...[DONE]
    % thicknessmode: ...[DONE]
    % len: ...[DONE]
    % lenmode: ...[DONE]
    % x: ...[DONE]
    % y: ...[DONE]
    % autotick: ...[DONE]
    % nticks: ...[DONE]
    % ticks: ...[DONE]
    % showticklabels: ...[DONE]
    % tick0: ...[DONE]
    % dtick: ...[DONE]
    % ticklen: ...[DONE]
    % tickwidth: ...[DONE]
    % tickcolor: ...[DONE]
    % tickangle: ...[NOT SUPPORTED IN MATLAB]
    % tickfont: ...[DONE]
    % exponentformat: ...[DONE]
    % showexponent: ...[NOT SUPPORTED IN MATLAB]
    % xanchor: ...[DONE]
    % yanchor: ...[DONE]
    % bgcolor: ...[DONE]
    % outlinecolor: ...[DONE]
    % outlinewidth: ...[DONE]
    % borderwidth: ...[NOT SUPPORTED IN MATLAB]
    % bordercolor: ...[NOT SUPPORTED IN MATLAB]
    % xpad: ...[DONE]
    % ypad: ...[DONE]

    %-FIGURE STRUCTURE-%
    figureData = get(obj.State.Figure.Handle);

    %-PLOT DATA STRUCTURE- %
    try
        colorbarData = get(obj.State.Colorbar(colorbarIndex).Handle);
    catch
        disp('could not get colorbar data');
    end

    %-STANDARDIZE UNITS-%
    colorbarUnits = colorbarData.Units;
    set(obj.State.Colorbar(colorbarIndex).Handle,'Units','normalized');

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

    %-variable initialization-%
    if isHG2
        outlineColor = [0 0 0];
    else
        if colorbarData.Position(4) > colorbarData.Position(3)
            outlineColor = 255*colorbarData.YColor;
        else
            outlineColor = 255*colorbarData.XColor;
        end
    end

    outlineColor = sprintf('rgb(%f,%f,%f)', outlineColor);
    lineWidth = colorbarData.LineWidth*obj.PlotlyDefaults.AxisLineIncreaseFactor;
    tickLength = min(obj.PlotlyDefaults.MaxTickLength,...
        max(colorbarData.TickLength(1)*colorbarData.Position(3)*obj.layout.width,...
        colorbarData.TickLength(1)*colorbarData.Position(4)*obj.layout.height));

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

    %-colorbar placement-%
    colorbar.x = colorbarData.Position(1);
    colorbar.y = colorbarData.Position(2);
    colorbar.len = colorbarData.Position(4);
    colorbar.thickness = colorbarData.Position(3);

    colorbar.xpad = obj.PlotlyDefaults.MarginPad;
    colorbar.ypad = obj.PlotlyDefaults.MarginPad;
    colorbar.xanchor = 'left';
    colorbar.yanchor = 'bottom';

    colorbar.outlinewidth = lineWidth;
    colorbar.outlinecolor = outlineColor;
    colorbar.exponentformat = obj.PlotlyDefaults.ExponentFormat;
    colorbar.thicknessmode = 'fraction';
    colorbar.lenmode = 'fraction';

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

    %-tick setings-%
    colorbar.tickcolor = outlineColor;
    colorbar.tickfont.color = outlineColor;
    colorbar.tickfont.size = colorbarData.FontSize;
    colorbar.tickfont.family = matlab2plotlyfont(colorbarData.FontName);
    colorbar.ticklen = tickLength;
    colorbar.tickwidth = lineWidth;

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

    %-get colorbar title and labels-%
    colorbarTitle = colorbarData.Label;

    if isHG2
        colorbarTitleData = get(colorbarTitle);
        colorbarYLabel = colorbarTitle;
        colorbarYLabelData = get(colorbarTitle);
        colorbarXLabelData.String = [];
    else
        colorbarTitleData = get(colorbarTitle);
        colorbarXLabel = colorbarData.XLabel;
        colorbarXLabelData = get(colorbarXLabel);
        colorbarYLabel = colorbarData.YLabel;
        colorbarYLabelData = get(colorbarYLabel);
    end

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

    %-STANDARDIZE UNITS FOR TITLE-%
    titleunits = colorbarTitleData.Units;
    titlefontunits = colorbarTitleData.FontUnits;
    ylabelunits = colorbarYLabelData.Units;
    ylabelfontunits = colorbarYLabelData.FontUnits;
    set(colorbarTitle,'Units','data');
    set(colorbarYLabel,'Units','data');
    set(colorbarYLabel,'FontUnits','points');
    if ~isHG2
        xlabelunits = colorbarXLabelData.Units;
        xlabelfontunits = colorbarXLabelData.FontUnits;
        set(colorbarTitle,'FontUnits','points');
        set(colorbarXLabel,'Units','data');
        set(colorbarXLabel,'FontUnits','points');
    end

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

    %-colorbar title settings-%
    isTitle = true;

    if ~isempty(colorbarTitleData.String)
        titleString = colorbarTitleData.String;
        titleInterpreter = colorbarTitleData.Interpreter;

        if colorbarTitleData.Rotation == 90
            titleSide = 'right';
        else
            titleSide = 'top';
        end

        titleFontSize = 1.20 * colorbarTitleData.FontSize;
        titleFontColor = sprintf('rgb(%f,%f,%f)', 255*colorbarTitleData.Color);
        titleFontFamily = matlab2plotlyfont(colorbarTitleData.FontName);

    elseif ~isempty(colorbarXLabelData.String)
        titleString = colorbarXLabelData.String;
        titleInterpreter = colorbarXLabelData.Interpreter;

        titleSide = 'right';
        titleFontSize = 1.20 * colorbarXLabelData.FontSize;
        titleFontColor = sprintf('rgb(%f,%f,%f)', 255*colorbarXLabelData.Color);
        titleFontFamily = matlab2plotlyfont(colorbarXLabelData.FontName);

    elseif ~isempty(colorbarYLabelData.String)
        titleString = colorbarYLabelData.String;
        titleInterpreter = colorbarYLabelData.Interpreter;

        titleSide = 'bottom';
        titleFontSize = 1.20 * colorbarYLabelData.FontSize;
        titleFontColor = sprintf('rgb(%f,%f,%f)', 255*colorbarYLabelData.Color);
        titleFontFamily = matlab2plotlyfont(colorbarYLabelData.FontName);

    else
        isTitle = false;
    end

    if isTitle
        colorbar.title = parseString(titleString, titleInterpreter);
        colorbar.titleside = titleSide;
        colorbar.titlefont.size = titleFontSize;
        colorbar.titlefont.color = titleFontColor;
        colorbar.titlefont.family = titleFontFamily;
    end

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

    %-REVERT UNITS FOR TITLE-%
    set(colorbarTitle,'Units',titleunits);
    set(colorbarTitle,'FontUnits',titlefontunits);
    set(colorbarYLabel,'Units',ylabelunits);
    set(colorbarYLabel,'FontUnits',ylabelfontunits);

    if ~isHG2
        set(colorbarXLabel,'Units',xlabelunits);
        set(colorbarXLabel,'FontUnits',xlabelfontunits);
    end

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

    %-tick labels-%
    tickValues = colorbarData.Ticks;
    tickLabels = colorbarData.TickLabels;
    showTickLabels = true;

    if isHG2
        if isempty(tickValues)
            showTickLabels = false;
            colorbar.ticks = '';

        elseif isempty(tickLabels)
            colorbar.tickvals = tickValues;

        else
            colorbar.tickvals = tickValues;
            colorbar.ticktext = tickLabels;
        end

        if showTickLabels
            colorbar.showticklabels = showTickLabels;

            switch colorbarData.AxisLocation
                case 'in'
                    colorbar.ticklabelposition = 'inside';
                case 'out'
                    colorbar.ticklabelposition = 'outside';
            end

            switch colorbarData.TickDirection
                case 'in'
                    colorbar.ticks = 'inside';
                case 'out'
                    colorbar.ticks = 'outside';
            end

        end

    else
        colorbar = setTicksNotHG2(colorbar, colorbarData);
    end

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

    %-colorbar bg-color-%
    if ~isHG2
        if ~ischar(colorbarData.Color)
            bgColor = 255*colorbarData.Color;
        else
            bgColor = 255*figureData.Color;
        end

        obj.layout.plot_bgcolor = sprintf('rgb(%f,%f,%f', bgColor);
    end

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

    %-ASSOCIATED DATA-%
    if isfield(colorbarData.UserData,'dataref')
        colorbarDataIndex = colorbarData.UserData.dataref;
    else
        colorbarDataIndex = findColorbarData(obj,colorbarIndex,colorbarData);
    end

    obj.data{colorbarDataIndex}.colorbar = colorbar;
    obj.data{colorbarDataIndex}.showscale = true;

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

    %-REVERT UNITS-%
    set(obj.State.Colorbar(colorbarIndex).Handle,'Units',colorbarUnits);

    %-------------------------------------------------------------------------%
end

function colorbar = setTicksNotHG2(colorbar, colorbarData)
    verticalOrientation = colorbar.len > colorbar.thickness;

    if verticalOrientation
        if isempty(colorbarData.YTick)
            %-show tick labels-%
            colorbar.ticks = '';
            colorbar.showticklabels = false;
        else
            %-tick direction-%
            switch colorbarData.TickDir
                case 'in'
                    colorbar.ticks = 'inside';
                case 'out'
                    colorbar.ticks = 'outside';
            end

            if strcmp(colorbarData.YTickLabelMode,'auto')
                %-autotick-%
                colorbar.autotick = true;
                %-numticks-%
                colorbar.nticks = length(colorbarData.YTick) + 1; %nticks = max ticks (so + 1)
            else
                %-show tick labels-%
                if isempty(colorbarData.YTickLabel)
                    colorbar.showticklabels = false;
                else
                    %-autotick-%
                    colorbar.autotick = false;
                    %-tick0-%
                    colorbar.tick0 = str2double(colorbarData.YTickLabel(1,:));
                    %-dtick-%
                    colorbar.dtick = str2double(colorbarData.YTickLabel(2,:)) - str2double(colorbarData.YTickLabel(1,:));
                end
            end
        end
    else
        if isempty(colorbarData.XTick)
            %-show tick labels-%
            colorbar.ticks = '';
            colorbar.showticklabels = false;
        else
            %-tick direction-%
            switch colorbarData.TickDir
                case 'in'
                    colorbar.ticks = 'inside';
                case 'out'
                    colorbar.ticks = 'outside';
            end

            if strcmp(colorbarData.XTickLabelMode,'auto')
                %-autotick-%
                colorbar.autotick = true;
                %-numticks-%
                colorbar.nticks = length(colorbarData.XTick) + 1;
            else
                %-show tick labels-%
                if isempty(colorbarData.XTickLabel)
                    colorbar.showticklabels = false;
                else
                    %-autotick-%
                    colorbar.autotick = false;
                    %-tick0-%
                    colorbar.tick0 = str2double(colorbarData.XTickLabel(1,:));
                    %-dtick-%
                    colorbar.dtick = str2double(colorbarData.XTickLabel(2,:)) - str2double(colorbarData.XTickLabel(1,:));
                end
            end
        end
    end
end

author: plotly


text:
function dat= func_powerspectrum (data, band, varargin)
% func_powerspectrum :
% This function calculates the power spectrum
%
% Example:
%      dat= func_powerspectrum(dat, [8 15], {'win',[];'step',100})
%
% Input:
%      dat  - data structure of continuous or epoched data
%      band - frequency band
%
% Options:
%      win  - window for FFT, default ones(dat.fs, 1)
%      N    - window width for FFT -> square window, default dat.fs
%      step - step for window (= # of overlapping samples), default N/2
%      db_scaled - boolean, if true values are db scaled (10*log10),
%                  default true
% Retuns:
%     dat - Data structure of power specturm result
%%
% data
dat=data;
ban=band;
opt=opt_cellToStruct(varargin{:});
epo=struct('win',[],'N',[],'step',[],'scale',[]);

if isempty(dat)
    warning('Warning! data is empty');
end
if isempty(ban)
    warning('Band is not exist.');
end

if ~isfield(opt,'win')
   epo.win=dat.fs;
else
    epo.win=opt.win;
end

if ~isfield(opt,'N')
   epo.N=dat.fs;
else
    epo.N=opt.N;
end

if ~isfield(opt,'step')
   epo.step=dat.fs/2;
else
   epo.N=opt.N;
end

if ~isfield(opt,'scale')
   epo.scale='db';
else
    epo.scale=opt.scale;
end

[T, nEvents , nChans]= size(dat.x);

if length(epo.win)==1
    if epo.win>T
        warning('window legth is higher than signal')
    end
    epo.win=ones(epo.win,1);
end
N=length(epo.win);
normwin=norm(epo.win);
Freq=(0:N)/2*dat.fs/N;

%%
XX= zeros(N, nChans*nEvents);
nWindows= 1 + max(0, floor((T-N)/epo.step));
iv= 1:min(N, T);
Win= repmat(epo.win(:), [1 nChans*nEvents]);
bInd= band(1): ban(2);

%%calculate file

switch(lower(epo.scale)),
    case 'db',
        for iw= 1:nWindows,
            XX= XX + abs(fft(dat.x(iv,:).*Win, N)).^2;
            iv= iv + epo.step;
        end
        XX = XX/(nWindows*normwin^2);
        dat.x= reshape( 10*log10( XX(bInd,:)+eps ), [length(bInd), nChans, nEvents]);
        dat.yUnit= 'dB';
    case 'power',
        for iw= 1:nWindows,
            XX= XX + abs(fft(dat.x(iv,:).*Win, N).^2);
            iv= iv + epo.step;
        end
        dat.x= reshape(XX(bInd,:)/(nWindows*normwin^2), [length(bInd), nChans, nEvents]);
        dat.yUnit= 'power';
end

end

author: PatternRecognition

repo: OpenBMI


text:
function B = subsref(A, S)
B = builtin('subsref', A, S);


text:
%% Percolation Cluster
clf
clc, clear
threshold = .63;
%
ax = axes('units','pixels','position',[1 1 500 400],'color','k');
text('units', 'pixels', 'position', [50,255,0],...
    'string','BioNB','color','w','fontname','helvetica','fontsize',100)
text('units', 'pixels', 'position', [120,120,0],...
    'string','441','color','w','fontname','helvetica','fontsize',100)
initial = getframe(gca);
[a,b,c]=size(initial.cdata);
z=zeros(a,b);
cells = double(initial.cdata(:,:,1)==255);
visit = z ;
sum = z;
imh = image(cat(3,z,cells,z));
set(imh, 'erasemode', 'none')
%return
for i=1:1000
    sum(2:a-1,2:b-1) = cells(2:a-1,1:b-2) + cells(2:a-1,3:b) + ...
    cells(1:a-2, 2:b-1) + cells(3:a,2:b-1) + ...
    cells(1:a-2,1:b-2) + cells(1:a-2,3:b) + ...
    cells(3:a,1:b-2) + cells(3:a,3:b);

    pick = rand(a,b);
    %edges only
    %cells = (cells & (sum<8)) | ((sum>=1) & (pick>=threshold) & (visit==0)) ;
    cells = cells | ((sum>=1) & (pick>=threshold) & (visit==0)) ;
    visit = (sum>=1) ;%& (pick<threshold) ;

    set(imh, 'cdata', cat(3,z,cells,z) )
    drawnow
end
return
figure(2)
image(cat(3,z,cells,z))

author: hliangzhao

repo: Mathematical-Model-Implementation


text:
function   v = eigenvector_centrality_und(CIJ)
%EIGENVECTOR_CENTRALITY_UND      Spectral measure of centrality
%
%   v = eigenvector_centrality_und(CIJ)
%
%   Eigenector centrality is a self-referential measure of centrality:
%   nodes have high eigenvector centrality if they connect to other nodes
%   that have high eigenvector centrality. The eigenvector centrality of
%   node i is equivalent to the ith element in the eigenvector
%   corresponding to the largest eigenvalue of the adjacency matrix.
%
%   Inputs:     CIJ,        binary/weighted undirected adjacency matrix.
%
%   Outputs:      v,        eigenvector associated with the largest
%                           eigenvalue of the adjacency matrix CIJ.
%
%   Reference: Newman, MEJ (2002). The mathematics of networks.
%
%   Contributors:
%   Xi-Nian Zuo, Chinese Academy of Sciences, 2010
%   Rick Betzel, Indiana University, 2012
%   Mika Rubinov, University of Cambridge, 2015

%   MODIFICATION HISTORY
%   2010/2012: original (XNZ, RB)
%   2015: ensure the use of leading eigenvector (MR)

n = length(CIJ);
if n < 1000
    [V,D] = eig(CIJ);
else
    [V,D] = eigs(sparse(CIJ));
end
[~,idx] = max(diag(D));
ec = abs(V(:,idx));
v = reshape(ec, length(ec), 1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef WheelEncoder < handle

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties
        type

        radius
        length
        ticks_per_rev

        ticks

        total_distance
    end

    methods
        function obj = WheelEncoder(type, radius, length, ticks_per_rev)
            obj.radius = radius;
            obj.length = length;
            obj.type = type;
            obj.ticks_per_rev = ticks_per_rev;
            obj.ticks = 0;
            obj.total_distance = 0;
        end

        function update_ticks(obj, wheel_velocity, dt)
            obj.ticks = obj.ticks + obj.distance_to_ticks(wheel_velocity*dt);
        end

        function reset_ticks(obj)
            obj.ticks = 0;
        end

        function ticks = distance_to_ticks(obj, distance)
            obj.total_distance = obj.total_distance + distance;
            ticks = round((obj.total_distance*obj.ticks_per_rev)/(2*pi));
            obj.total_distance = obj.total_distance - obj.ticks_to_distance(ticks);
        end

        function distance = ticks_to_distance(obj, ticks)
            distance = (ticks*2*pi)/obj.ticks_per_rev;
        end
    end
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function prob_test030 ( )

%*****************************************************************************80
%
%% TEST030 tests CHI_CDF, CHI_CDF_INV, CHI_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST030\n' );
  fprintf ( 1, '  For the Chi PDF:\n' );
  fprintf ( 1, '  CHI_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  CHI_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  CHI_PDF evaluates the PDF.\n' );

  a = 1.0;
  b = 2.0;
  c = 3.0;

  check = chi_check ( a, b, c );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST030 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );
  fprintf ( 1, '  PDF parameter C = %14f\n', c );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed  ] = chi_sample ( a, b, c, seed );

    pdf = chi_pdf ( x, a, b, c );

    cdf = chi_cdf ( x, a, b, c );

    x2 = chi_cdf_inv ( cdf, a, b, c );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function apAll = evalRP_multicut(expidxs, bSave, subset)

fprintf('evalRP_multicut()\n');

if (nargin < 2)
    bSave = true;
end

if (nargin < 3)
    subset = [];
end

fontSize = 18;
% figure(100); clf; hold on;
legendName = cell(0);
set(0,'DefaultAxesFontSize', fontSize)
set(0,'DefaultTextFontSize', fontSize)

prefix = 'pck-head-rpc-';
[~,parts] = util_get_parts24();
nJoints = 14;
pck_thresh = 0.5;

% clusSizeThresh = 15;

num_exps = length(expidxs);
imgidxs_missing_all = cell(num_exps, 1);
keypointsAll_all = cell(num_exps, 1);
annolist_gt_all = cell(num_exps, 1);
annolist_bbox_all = cell(num_exps, 1);

for ii = 1:num_exps
    fprintf('**********************************************************************************\n');
    expidx = expidxs(ii);

    % load experiment parameters
    p = exp_params(expidx);

    load(p.testGT);
    annolist_gt = annolist;%(1:300);

    mkdir_if_missing(p.latexDir);
    mkdir_if_missing(p.plotsDir);

    assert(length(p.pidxs) == nJoints);

    if (isfield(p,'clusSizeThresh'))
        clusSizeThresh = p.clusSizeThresh;
    else
        clusSizeThresh = 1;
    end

    if (isfield(p,'disjointGTassignment'))
        disjointGTassignment = p.disjointGTassignment;
    else
        disjointGTassignment = false;
    end

    % load keypoints
    if (isfield(p,'multicut') && p.multicut == true)
        [keypointsAll,imgidxs_missing] = det2keypointsImgMulticutMulti(expidx,parts,annolist_gt,clusSizeThresh);
    elseif (isfield(p,'oracleDetBBox') && p.oracleDetBBox == true && ~isfield(p,'poseClusters'))
        annolist_bbox = compute_gt_bbox(annolist_gt,p.pidxs,parts,p.bbox_offset);
        [keypointsAll,imgidxs_missing] = det2keypointsImgMulticutBBoxMulti(expidx,parts,annolist_bbox,disjointGTassignment);
    elseif (isfield(p,'neckDetBBox') && p.neckDetBBox == true && ~isfield(p,'poseClusters'))
        try
            assert(false);
            load([p.outputDir '/annolist_bbox'],'annolist_bbox');
        catch
            annolist_bbox = compute_det_bbox(expidx,annolist_gt,p.neck_nms_thresh,p.neck_det_thresh,p.bbox_offset);
            if (~exist(p.outputDir,'dir')),mkdir(p.outputDir);end
            save([p.outputDir '/annolist_bbox'],'annolist_bbox');
        end
        [keypointsAll,imgidxs_missing] = det2keypointsImgMulticutBBoxMulti(expidx,parts,annolist_bbox,disjointGTassignment);
    elseif (isfield(p,'poseClusters') && p.poseClusters == true)
        load(p.testGTnopad);
        annolist_gt = annolist(1:length(annolist_gt));
        [keypointsAll,imgidxs_missing] = det2keypointsImgMulticutBBoxMultiPoses(expidx,parts,length(annolist_gt));
    end

    imgidxs_missing_all{ii} = imgidxs_missing;
    keypointsAll_all{ii} = keypointsAll;
    annolist_gt_all{ii} = annolist_gt;
    %annolist_bbox_all{ii} = annolist_bbox;
end

imgidxs_missing = [];

for ii = 1:num_exps
    imgidxs_missing = [imgidxs_missing; imgidxs_missing_all{ii}];
end

imgidxs_missing = unique(imgidxs_missing);

%im_idx = 7;
%imgidxs_missing = [1:(im_idx-1) (im_idx+1):1758];

fprintf('final images missing %d\n', length(imgidxs_missing));

for ii = 1:num_exps
    expidx = expidxs(ii);
    p = exp_params(expidx);

    keypointsAll = keypointsAll_all{ii};
    annolist_gt = annolist_gt_all{ii};
    %annolist_bbox = annolist_bbox_all{ii};

    keypointsAll16 = keypointsAll;
    if isempty(subset)
        subset = true(length(annolist_gt), 1);
        subset(imgidxs_missing) = 0;
    end
    annolist_gt = annolist_gt(subset);
    keypointsAll = keypointsAll(subset);

    assert(length(annolist_gt) == length(keypointsAll));
    [scoresAll, labelsAll, nGTall] = assignGTMulticutMulti(keypointsAll,annolist_gt,p.pidxs,parts,pck_thresh);

%     groups = {0,2,4,5,7,9,12,14,16,17,19,21,22,23};
    partLabels = {'rankle','rknee','rhip','lhip','lknee','lankle','rwrist','relbow','rshoulder','lshoulder','lelbow','lwrist','neck','tophead'};

    apAll = zeros(length(p.pidxs),1);
%     figure(100); clf;

    for i = 1:length(p.pidxs)
        scores = [];
        labels = [];
        for imgidx = 1:length(annolist_gt)
            scores = [scores; scoresAll{i}{imgidx}];
            labels = [labels; labelsAll{i}{imgidx}];
        end
        [precision,recall,sorted_scores,sortidx,sorted_labels] = getRPC(scores,labels,sum(nGTall(i,:)));
        ap = VOCap(recall,precision)*100;
        apAll(i) = ap;
        legendName = sprintf('%s, AP: %1.1f%%', p.name, ap);
        fprintf('%s\n',legendName);

%         plotRPC(precision,recall,p.colorName,'-',partLabels{i});
%         legend(legendName,'Location','NorthWest');
%         if (bSave)
%             print(gcf, '-dpng', [p.plotsDir '/' prefix partLabels{i} '-' 'expidx' num2str(expidx) '.png']);
%             printpdf([p.plotsDir '/' prefix partLabels{i} '-' 'expidx' num2str(expidx) '.pdf']);
%         end
    end

    num_images = length(labelsAll{1});
    num_parts = length(labelsAll);
    simple_scores = zeros(num_images, 1);
    for k = 1:num_images
        num_hits = 0;
        num_all = 0;
        for i = 1:num_parts
            lbl = labelsAll{i}{k};
            num_all = num_all + length(lbl);
            num_hits = num_hits + sum(lbl);
        end

        if num_all ~= 0
            simple_scores(k) = num_hits/num_all;
        end
    end

    mkdir_if_missing(fullfile(p.exp_dir, 'data'));
    save(fullfile(p.exp_dir, 'data', 'scores'), 'simple_scores', 'subset', 'keypointsAll16');

    if (bSave)
        tableFilename = [p.latexDir '/ap-expidx' num2str(expidx) '.tex'];
        apAll(15,:) = mean(apAll(7:14));
        apAll(16,:) = mean(apAll);
        [row, header] = genTable(apAll,p.name);
        fid = fopen(tableFilename,'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',row{1});fclose(fid);
        fid = fopen([p.latexDir '/ap-header.tex'],'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',header);fclose(fid);
        save([p.latexDir '/ap-expidx' num2str(expidx) '.mat'],'apAll');
    end
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function pot = cg_can_to_mom(pot)
% CG_CAN_TO_MOM Convert a CG potential from canonical to moment form, if necessary.
% pot = cg_can_to_mom(pot)

if pot.subtype ~= 'm'
  for i=1:pot.dsize
    pot.mom{i} = cpot_to_mpot(pot.can{i});
  end
  pot.subtype = 'm';
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% bootstrapping trial

%% make test Data matrix (from data examined by hand)
Data = cell(1,3);

i_fish = 8;
gIX = FC{i_fish}.gIX;
IX = find(gIX <100);
Data{1}.cIX = FC{i_fish}.cIX(IX);
Data{1}.xyz_norm = FC{i_fish}.xyz_norm(IX,:);

i_fish = 9;
Data{2}.cIX = FC{i_fish}.cIX;
Data{2}.xyz_norm = FC{i_fish}.xyz_norm;

i_fish = 10;
gIX = FC{i_fish}.gIX;
IX = find(gIX <100);
Data{3}.cIX = FC{i_fish}.cIX(IX);
Data{3}.xyz_norm = FC{i_fish}.xyz_norm(IX,:);

%%

%%
thres_prc = 5;
numDim = 3;
numFish = 3;

TF = zeros(numFish,numFish,numDim);

for i_ref = 1:numFish,
    numRand = 100;
    means = zeros(numRand,3);
    for i = 1:numRand,
        IX = ceil(100*rand(1,numRand));
        means(i,:) = mean(Data{i_ref}.xyz_norm(IX,:),1);
    end
    for i_test = 1:numFish,
        if i_test ~= i_ref,
            for i_dim = 1:numDim,
                mean_test = mean(Data{i_test}.xyz_norm(:,i_dim),1);
                lowerbound = prctile(means(:,i_dim),thres_prc/2);
                upperbound = prctile(means(:,i_dim),100-thres_prc/2);
                TF(i_ref,i_test,i_dim) = mean_test>lowerbound & mean_test<upperbound;
            end
        else
            TF(i_ref,i_test,:) = NaN;
        end
    end
end

author: xiuyechen

repo: FishExplorer


text:
% TEST_MAXWELL_RING_H_DRCHLT: data function for Dirichlet boundary condition.

function h = test_maxwell_ring_h_drchlt (x, y, ind)

  [theta, r] = cart2pol (x, y);
  h = zeros (size(x));
  switch (ind)
    case 1
      h = -sin(theta) .* sin(y) + cos(theta) .* sin(x);
    case 2
      h = sin(theta) .* sin(y) - cos(theta) .* sin(x);
    case 3
      h = -sin(y);
    case 4
      h = sin(x);
    otherwise
      error ('h_drchlt: unknown reference number')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [Q, maxIter] = powerMethodMatComp( U, V, spa, R, maxIter, tol)
% low rank (U*V') + sparse term (data)

if(isempty(R))
    maxIter = 0;
    Q = randn(size(U,1), 1);
    return;
end

Y = (R'*spa')';
Y = U*(V'*R) + Y;
[Q, ~] = qr(Y, 0);
err = zeros(maxIter, 1);
for i = 1:maxIter
    % Y = A*(A'*Q);
    AtQ = (Q'*spa)';
    AtQ = V*(U'*Q) + AtQ;
    Y = (AtQ'*spa')';
    Y = U*(V'*AtQ) + Y;
    [iQ, ~] = qr(Y, 0);

    err(i) = norm(iQ(:,1) - Q(:,1), 2);
    Q = iQ;

    if(err(i) < tol)
        break;
    end
end

maxIter = i;

end

author: HKUST-KnowComp

repo: FMG

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


text:
function newmodel = addMultipleReactions(model,rxnIDs,metList,Stoichiometries,varargin)
% Add multiple Reactions to the model. In contrast to addReaction, this
% function allows the addition of multiple reactions at once but requires
% all used metabolites to be already present in the model.
%
% USAGE:
%
%    newmodel = addMultipleReactions(model,rxnIDs,metList,Stoichiometries,varargin)
%
% INPUTS:
%    model:             The model to add the Metabolite batch to.
%    rxnIDs:            The IDs of the reactions that shall be added (cell Array).
%    metList:           A list of metabolites for which stoichiometric coefficients are given
%    Stoichiometries:   A Stoichiometric matrix of dimension
%                       numel(metList) x numel(rxnID).
% OPTIONAL INPUTS:
%    varargin:          fieldName, Value pairs. The given fields will be set
%                       according to the values. Only defined COBRA fields may be
%                       used. The following fields will be ignored (as they
%                       are dependent on the existing model structure)
%
%                        - `S`:  this is being resolved by the metList/Stoichiometries combination
%                        - `rxnGeneMat`: This depends on the original model structure, and is thus nor considered.
%
%                       You can provide GPR rules in different ways to
%                       this function:
%
%                        - By providing a grRules array in the varargin
%                        - By providing a genes field in the varargin and a
%                          corresponding rules field
%                        - By providing a rules field, that directly refers
%                          to the genes in the model.
%
%                       Examples of parameters are:
%
%                        - 'rxnNames',{'RuBisCo','Transketloase'}
%                        - 'rxnKEGGID,{'R00024','R01641'}
%                        - 'lb',[0,-1000]
%                        - or any field name associated with reactions (except
%                          for S and rxnGeneMat) as defined in the COBRA
%                          Model field definitions.
%                       It is also possible to add a printout to the
%                       function by using the 'printLevel', 1 argument.
%
%
% OUTPUTS:
%
%    newmodel:     The model structure with the additional reactions.
%
% EXAMPLE:
%
%    % To add the following reactions, with lower and upper bounds:
%    % ExA: A <-> ; ATob: A <-> B ; BToC: B <-> C
%    model = addMultipleReactions(model, {'ExA','ATob','BToC'}, {'A','b','c'}, [1 -1 0; 0,1,-1;0,0,1], 'lb', [-50,30,1], 'ub', [0,60,15])
%    % To add them with GPRs in text form:
%    % ExA: A <-> ; ATob: A <-> B ; BToC: B <-> C
%    model = addMultipleReactions(model, {'ExA','ATob','BToC'}, {'A','b','c'}, [1 -1 0; 0,1,-1;0,0,1], 'grRules', {'G1 or G2', 'G3 and G4',''})
%
%
%    % To add them with the same GPRs in logical format assuming that model.genes is {'G1';'G2';'G3';'G4'}:
%    % ExA: A <-> ; ATob: A <-> B ; BToC: B <-> C
%    model = addMultipleReactions(model, {'ExA','ATob','BToC'}, {'A','b','c'}, [1 -1 0; 0,1,-1;0,0,1], 'rules', {'x(1) | x(2)', 'x(3) and x(4)',''})
%    % To add them with the same GPRs in logical format without assuming anything for the model:
%    % ExA: A <-> ; ATob: A <-> B ; BToC: B <-> C
%    model = addMultipleReactions(model, {'ExA','ATob','BToC'}, {'A','b','c'}, [1 -1 0; 0,1,-1;0,0,1], 'rules', {'x(3) | x(2)', 'x(4) & x(1)',''}, 'genes', {'G4';'G2';'G1';'G3'})
%
%    % Print the above:
%    model = addMultipleReactions(model, {'ExA','ATob','BToC'}, {'A','b','c'}, [1 -1 0; 0,1,-1;0,0,1], 'printLevel', 1, 'rules', {'x(3) | x(2)', 'x(4) & x(1)',''}, 'genes', {'G4';'G2';'G1';'G3'})

[metPres,metPos] = ismember(metList,model.mets);
if any(metPos == 0)
    error('The following Metabolites are not part of the model:\n%s',strjoin(metList(metPos==0)));
end

% check rxnIDs
if ~iscell(rxnIDs)
    error('rxnIDs has to be a cell array of strings!')
end

%[bool,existing] = checkIDsForTypeExist(model,ids,basefield)
% Check whether the given IDs exist for the given type of base field. Base
% fields include rxns, mets, genes, comps, proteins, ctrs, evars. ctrs/mets
% as well as rxns/evars will be considered as a combined field.
% OUTPUT:
%    bool:                Boolean vector true if any of the IDs exist.
%    existing:            Unique set of pre-existing IDs for this base field.
[bool,existing] = checkIDsForTypeExist(model,rxnIDs,'rxns');

if bool
    if exist('varargin','var')
        warning('The following reaction IDs are already IDs of variables in the model:\n%s', strjoin(existing,'\n'));
    else
        %try to recover rather than fail ungracefully, passes
        %testBatchAddition and testMergeTwoModels!

        % C = setdiff(A,B) for vectors A and B, returns the values in A that
        % are not in B with no repetitions. C will be sorted.
        rxnIDs = setdiff(rxnIDs,existing);
        if isempty(rxnIDs)
            warning('All of the reaction IDs provided are already IDs of reactions in the model, so nothing to add.:\n%s', strjoin(existing,'\n'));
            newmodel = model;
            return
        else
            duplicateRxnIDBool = ismember(existing,rxnIDs);
            %pare down the stoichiometries
            Stoichiometries = Stoichiometries(:,duplicateRxnIDBool);
            warning('The following reaction IDs are already IDs of reactions in the model and will not be added:\n%s', strjoin(existing,'\n'));
        end
    end
end

if numel(unique(metList)) < numel(metList)
    error('Duplicate Metabolite ID detected.');
end

printLevel = 0;

nRxns = numel(model.rxns);

%We extract those fields which are either associated with the rxns field
%(2nd or 3rd column contains 'rxns' in the definitions), and we also look up which fields are in the model and
%associated with the rxns field (from the sizes and definitions)
fieldDefs = getDefinedFieldProperties();
fieldDefs = fieldDefs(cellfun(@(x) strcmp(x,'rxns'), fieldDefs(:,2)) | cellfun(@(x) strcmp(x,'rxns'), fieldDefs(:,3)));
modelRxnFields = getModelFieldsForType(model,'rxns');
model.rxns= [model.rxns;columnVector(rxnIDs)];

for field = 1:2:numel(varargin)
    cfield = varargin{field};
    values = varargin{field+1};
    %Anything thats not a model field or not a specialised field is
    %ignored.
    if strcmp(cfield,'printLevel')
        printLevel = varargin{field+1};
        continue;
    end
    if any(ismember({'rules','genes'},cfield))
        continue; %Don't throw a warning here.
    end
    if any(ismember({'S','rxnGeneMat'},cfield)) || (~any(ismember(fieldDefs(:,1),cfield)) && ~any(ismember(modelRxnFields,cfield)))
        warning('Field %s is excluded.',cfield);
        continue;
    end
    % if we reach this point we can transpose it savely
    values = columnVector(values);
    if any(size(values) ~= size(columnVector(rxnIDs)))
        if size(values,2) ~=1
            error('The supplied field %s has dimension %ix%i. But should be a %ix1 vector',cfield,size(values),numel(rxnIDs));
        else
            error('The supplied field %s has %i elements but it should have %i elements, one for each element in rxnIDs.',cfield,numel(values),numel(rxnIDs));
        end
    end
    if ~isfield(model,cfield)
        model = createEmptyFields(model,cfield);
        model.(cfield)((end-numel(varargin{field+1})+1):end) = columnVector(values);
    else
        model.(cfield) = [model.(cfield);columnVector(values)];
    end
end
newmodel = extendModelFieldsForType(model,'rxns','originalSize',nRxns);
%Now, we have extended the S matrix. So lets fill it.
newmodel.S(metPos,(nRxns+1):end) = Stoichiometries;

if isfield(model,'grRules') && ~any(ismember(varargin(1:2:end),'rules')) %There is a grRules field, and no Rules are provided.
    rulesToUpdate = ~cellfun(@isempty, model.grRules(nRxns+1:end));
    if any(rulesToUpdate)
        %We have non Empty grRules in the input. Parse them.
        newgrRules = nRxns + find(rulesToUpdate);
        for i = 1:numel(newgrRules)
            newmodel = changeGeneAssociation(newmodel,newmodel.rxns{newgrRules(i)},newmodel.grRules{newgrRules(i)});
        end
    end
end

if any(ismember(varargin(1:2:end),'rules'))
    %We have a rules field. First, reset any grRules
    if isfield(model,'grRules')
        newmodel.grRules(nRxns+1:end) = {''};
    end
    %Now, we have to check, if there is a gene list provided:
    if any(ismember(varargin(1:2:end),'genes'))
        %If genes are provided, we have to A: find the overlap, be
        %merge the genes and C adjust the rules.
        nGenes = length(newmodel.genes);
        genePos = 2 * find(ismember(varargin(1:2:end),'genes'));
        rulesPos = 2 * find(ismember(varargin(1:2:end),'rules'));
        genes = varargin{genePos};
        rules = varargin{rulesPos};
        if ~any(size(rules) == numel(rxnIDs)) %something must fit
            error('Size of field rules does not fit to the rxnList size');
        end
        %Find all genes (and their positions) which are already in the model.
        [genePres,genePos] = ismember(genes,newmodel.genes);
        additionalGenes = sum(~genePres);
        newmodel = addGenes(newmodel,genes(~genePres));
        genePos(~genePres) = nGenes+(1:additionalGenes);
        %now, replace the ids.
        genePos = cellfun(@(x) strcat('x(', num2str(x),')'), num2cell(genePos),'Uniformoutput',false);
        rules = regexprep(rules,'x\(([0-9]+)\)','${genePos{str2num($1)}}');
        newmodel.rules(nRxns+1:length(newmodel.rxns)) = rules;
    else
        rulesPos = 2 * find(ismember(varargin(1:2:end),'rules'));
        if ~isfield(newmodel,'rules') %If rules does not exist, create it.
            newmodel = generateRules(newmodel);
        end
        if ~any(size(varargin{rulesPos}) == numel(rxnIDs)) %something must fit
            error('Size of field rules does not fit to the rxnList size');
        end
        newmodel.rules((nRxns+1):end) = columnVector(varargin{rulesPos});
    end

    %Update the corresponding grRules fields and rxnGeneMat fields.
    %create the grRules
    rulesToUpdate = nRxns + find(~cellfun(@isempty, model.rules(nRxns+1:end)));
    if ~isempty(rulesToUpdate)
        if isfield(newmodel,'grRules')
            cGrRules = newmodel.grRules(nRxns+1:end);
            cGrRules = strrep(cGrRules,'|','or');
            cGrRules = strrep(cGrRules,'&','and');
            cGrRules = regexprep(cGrRules,'x\(([0-9]+)\)','${model.genes{str2num($1)}}');
            newmodel.grRules(nRxns+1:end) = cGrRules;
        end
        %Also update the rxnGeneMat, if present
        if isfield(model,'rxnGeneMat')
            nGenes = length(newmodel.genes);
            for i = 1:numel(rulesToUpdate)
                assoc = false(1,nGenes);
                %get all gene positions
                pos = regexp(newmodel.rules{rulesToUpdate(i)},'x\((?<pos>[0-9]+)\)','names');
                genePos = cellfun(@str2num, {pos.pos});
                assoc(genePos) = true;
                if ~issparse(newmodel.rxnGeneMat)
                    newmodel.rxnGeneMat=sparse(newmodel.rxnGeneMat);
                end
                newmodel.rxnGeneMat(rulesToUpdate(i),:) = sparse(assoc);
            end
        end
    end

end
if printLevel > 0
    fprintf('addMultipleReactions: Adding the following reactions to the model:\n');
    printRxnFormula(newmodel,rxnIDs);
end

author: opencobra

repo: cobratoolbox


text:
function value = r8_log_2 ( x )

%*****************************************************************************80
%
%% R8_LOG_2 returns the logarithm base 2 of |X|.
%
%  Discussion:
%
%    R8_LOG_2 ( X ) = Log ( |X| ) / Log ( 2.0 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose base 2 logarithm is desired.
%    X should not be 0.
%
%    Output, real VALUE, the logarithm base 2 of the absolute
%    value of X.  It should be true that |X| = 2**R8_LOG_2.
%
  if ( x == 0.0 )
    value = -inf;
  else
    value = log ( abs ( x ) ) / log ( 2.0 );
  end

  return
end


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function checkStartup()
  OCL_CASADI_SETUP = getenv('OCL_CASADI_SETUP');

  if strcmp(OCL_CASADI_SETUP, 'true')
    ocl_casadi_setup_completed = true;
  else
    ocl_casadi_setup_completed = false;
  end

  if isempty(ocl_casadi_setup_completed) || ~ocl_casadi_setup_completed
    disp('Running OpenOCL setup procedure. This may required your input, and may take a while at the first time.')
    ocl.utils.startup();
    setenv('OCL_CASADI_SETUP', 'true');
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function HDF5_GroupClose(fileobj, gid, fid)

H5G.close(gid);
if ischar(fileobj)
    H5F.close(fid);
end


text:
function candoit = canGetPartialGradient(problem)
% Checks whether the partial gradient can be computed for a given problem.
%
% function candoit = canGetPartialGradient(problem)
%
% Returns true if the partial gradient of the cost function can be computed
% given the problem description, false otherwise.
%
% See also: getPartialGradient canGetPartialEuclideanGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 28, 2016.
% Contributors:
% Change log:

    candoit = (isfield(problem, 'partialgrad') && isfield(problem, 'ncostterms')) || ...
              canGetPartialEuclideanGradient(problem);

    if isfield(problem, 'partialgrad') && ~isfield(problem, 'ncostterms')
        warning('manopt:partialgrad', ...
               ['If problem.partialgrad is specified, indicate the number n\n' ...
                'of terms in the cost function with problem.ncostterms = n.']);
    end

end

author: NicolasBoumal

repo: manopt


text:
function neqn = p11_neqn ( )

%*****************************************************************************80
%
%% P11_NEQN returns the number of equations for problem P11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 10;

  return
end