text:
function helpstring = hlp_HelpToString(func,args,maxwrapcols)
% Automatically create Help Text from the arg specification for a given
% function. This will only generate help for arguments that are reportable
% (i.e., any  arg_norep()  arguments will be ignored). This program
% can only be applied to functions which use Christian Kothe's argument
% input/output specification (http://sccn.ucsd.edu/wiki/BCILAB).
%
% Inputs:
%
%   func:   the handle to the function for which to generate helptext
%   args:   cell array of mandatory argument inputs for the function
%   maxwrapcols: the max number of columns before wrapping. Use Inf for
%                no wrapping
% Output:
%
%   helpstring: formatted help text which can be pasted as a function
%               header.
%
% Example:
%
% helpstring = hlp_HelpToString(@pre_prepData,struct([]))
%
% helpstring =
%
% Optional                 Information
% -------------------------------------------------------------------------------------------------------------------------------
% VerbosityLevel:          Verbosity level. 0 = no output, 1 = text, 2 = graphical
%                          Possible values: 0,1,2
%                          Default value  : 0
%                          Input Data Type: real number (double)
%
% backupOriginalData:      Keep a nonnormalized copy of the data
%                          Input Data Type: boolean
%
% SelectComponents:        Select components to analyze
%                          Input Data Type: string
%
%     | VerbosityLevel:    Verbosity level. 0 = no output, 1 = text, 2 = graphical
%                          Possible values: 0,1,2
%                          Default value  : 0
%                          Input Data Type: real number (double)
%
%     | ComponentsToKeep:  Select components to analyze
%                          This should be a cell array of strings containing the IDs of components you wish to keep
%                          Input Data Type: boolean
%
% EpochTimeRange:          [Min Max] Epoch time range (sec)
%                          If blank, use all original epoch length
%                          Input Data Type: real number (double)
%
% NewSamplingRate:         New sampling rate
%                          Data will be down/upsampled using a zero-phase filter (see 'help resample')
%                          Input Data Type: real number (double)
% ...
%
%
% See Also: arg_report()
%
% Author: Tim Mullen 2011, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

if nargin<3
    maxwrapcols = 100;
end

% get the full argument spec for this function
rep = arg_report('rich',func,args);

helpstring = [{'Input','Information'}; parseHelpText(rep,0,maxwrapcols)];

helpstring = cprintf([helpstring cell(size(helpstring,1), 1)],'-la',true,'-it',true,'-E','','-Lcs','-');
helpstring = [repmat('-',[1 size(helpstring,2)]); helpstring];

function output = parseHelpText(rep,indentlevel,maxwrapcols)
% This functon generates formatted help text for an argument (and its children).
%
% Inputs:
%         rep:              The rich report for an argument
%         indentlevel:      The depth level for this argument
%         maxwrapcols:      The max number of columns before wrapping. Use Inf for
%                           no wrapping
%
% Outputs:
%         output:           Formatted help text

% indent delimiter
delim = '    ';

if nargin < 3
    maxwrapcols = 100;
end
if nargin < 2
    indentlevel = 0;
end

output = {};

% iterate through all arguments and generate help text
for i=1:length(rep)

    if strcmpi(rep(i).names{1},'arg_selection')
        continue;
    end

    % extract the name of the entry
    name = rep(i).names{fastif(length(rep(i).names)==1,1,2)};

    % extract the help text for this entry
    % apply textwrapping if necessary
    helptext = {};
    rep(i).help = strtrim(rep(i).help);
    for j=1:length(rep(i).help)
        helptext = [helptext; textwrap([],rep(i).help(j),maxwrapcols)];
    end

    % add the data range and default value to helptext
    range = rep(i).range;
    def   = rep(i).value;
    if isempty(def), def='n/a'; end

%     if ~strcmpi(def,'__arg_mandatory__') %~isempty(range)

    def = fastif(strcmpi(def,'__arg_mandatory__'),'MANDATORY INPUT',def);

    if isempty(range), range = 'Unrestricted'; end

    if iscell(rep(i).range)
        if ~isempty(rep(i).range), range = cell2str(rep(i).range); end %['''' strrep(cell2str(rep(i).range),',',''',''') ''''];
        def = cell2str(def);
        helptext = [helptext; ...
                    sprintf('Possible values: %s', range); ...
                    sprintf('Default value  : %s', def)];
    elseif isnumeric(rep(i).range)
        if ~isempty(rep(i).range), range = ['[' num2str(rep(i).range) ']']; end
        helptext = [helptext; ...
                    sprintf('Input Range  : %s', range); ...
                    sprintf('Default value: %s', fastif(ischar(def),def,num2str(def)))];
    else
        if ~isempty(rep(i).range), range = char(rep(i).range); end
        helptext = [helptext; ...
                    sprintf('Possible values: %s', range); ...
                    sprintf('Default value  : ''%s''', def)];
    end

%     else
%         helptext = [helptext; ...
%                         'MANDATORY INPUT'];
%     end

    % add the input type to help text
    type  = simplifyDataType(rep(i).type);
    helptext = [helptext; sprintf('Input Data Type: %s',type)];

    % construct help text for this entry
%     output = cprintf({[name ':'], arg.(name)},'-d','\t','-cla',true,'-O',@(x) cprintf(x,'-E','','-d','\t','-la',true), ...
%                          '-L',repmat('\t',[1 indentlevel]));
%     output = [output; {name,cprintf(char(helptext),'-cla',true,'-O',@(x) cprintf(x,'-E','','-d','\t','-la',true),'-L',repmat('\t',[1 indentlevel])) } ];

    % create blank leading "row"
    output = [output; {'',''}];

    % fill row with [Name:   Help Text]
    output = [output; {[name ':'],char(helptext)} ];

%     % check if this has children and if so, create underline
%     % TM: omitted
%     if ismember_bc(char(rep(i).head),{'arg_sub','arg_subswitch','arg_subtoggle'})
%         output = [output; {repmat('-',1,length(name)+1), ''}];
%     end

    % now parse children and insert helptext as needed
    if ~isempty(rep(i).alternatives)
        for k=2:length(rep(i).alternatives)
            for child=1:length(rep(i).alternatives{k})
                ret    =  parseHelpText(rep(i).alternatives{k}(child),indentlevel,maxwrapcols);
                if ~isempty(ret)
                    for q=1:size(ret,1)
                        % insert indented child argument block with pipe
                        % character prepended
                        output = [output; {[repmat(delim,[1 indentlevel+1]) fastif(isempty(ret{q,1}),'','| ') ret{q,1}] ret{q,2}}];
                    end
                end
            end
        end
    end

    if ~isempty(rep(i).children)
        % parse children help text
        for child=1:length(rep(i).children)
            ret    =  parseHelpText(rep(i).children(child),indentlevel,maxwrapcols);
            if ~isempty(ret)
                for q=1:size(ret,1)
                    % insert indented child argument block with pipe
                    % character prepended
                    output = [output; {[repmat(delim,[1 indentlevel+1]) fastif(isempty(ret{q,1}),'','| ') ret{q,1}] ret{q,2}}];
                end
            end
        end
    end

end

% convert cell array of numbers or other data to a comma-delimited string
% numeric contents are converted to form a,b,c
% char contents are wrapped with quotes in form 'a','b','c'
% all other datatypes are converted to form a,b,c
function c = cell2str(C)

c = '';

if ~iscell(C)
    C = {C};
end

for k=1:length(C)
    if isnumeric(C{k})
        c = [c num2str(C{k}) ','];
    elseif ischar(C{k})
        c = [c '''' C{k} ''','];
    else
        c = [c char(C{k}) ','];
    end
end

c(end) = []; % delete trailing comma

% convert coded datatypes into more human-readable expressions
function plaintext = simplifyDataType(datatype)

    switch datatype
        case 'logical'
            plaintext = 'boolean';
        case 'char'
            plaintext = 'string';
        case 'int8'
            plaintext = 'integer (int8)';
        case 'uint8'
            plaintext = 'integer (uint8)';
        case 'int16'
            plaintext = 'integer (int16)';
        case 'uint16'
            plaintext = 'integer (uint16)';
        case 'int32'
            plaintext = 'integer (int32)';
        case 'int64'
            plaintext = 'integer (int64)';
        case 'uint64'
            plaintext = 'integer (uint64)';
        case 'denserealsingle'
            plaintext = 'real number (single)';
        case 'denserealdouble'
            plaintext = 'real number (double)';
        case 'densecomplexsingle'
            plaintext = 'complex number (single)';
        case 'densecomplexdouble'
            plaintext = 'complex number (double)';
        case 'sparserealsingle'
            plaintext = 'real number (sparse single)';
        case 'sparserealdouble'
            plaintext = 'real number (sparse double)';
        case 'sparsecomplexsingle'
            plaintext = 'complex number (sparse single)';
        case 'sparsecomplexdouble'
            plaintext = 'complex number (double)';
        case 'cellstr'
            plaintext = 'cell array of strings (cellstr)';
        case 'object'
            plaintext = 'any matlab object';
        case 'expression'
            plaintext = 'any evaluable Matlab expression.';
        otherwise
            plaintext = datatype;
    end

repo: REST


text:
%rate = mimo_iso_conv(nn,P,error,tx,rx,K) computes the converse bound on
%the maximal channel coding rate over a quasi-static MIMO Rician fading
%channel with isotropic inputs.
%
%The function takes the following inputs:
%
%nn: blocklength, scalar/vector;
%P: total input power (linear scale), scalar;
%error: block error probability, scalar;
%tx: number of transmit antennas, scalar;
%rx: number of receive antennas, scalar;
%K:  rician K factor, scalar; the default value is 0.

function rate_c = mimo_iso_conv(nn,P,error,tx,rx,K)
if (nargin < 6) || isempty(K)
	K = 0;
end
loop_h=round(1000/error);

P=P/tx;

rate_c=[];

min_tx_rx = min(tx,rx);

eig_vec=zeros( min_tx_rx, loop_h);

for ij = 1:loop_h
    h = (randn(tx,rx)+1i*randn(tx,rx))*sqrt(0.5/(K+1)) + sqrt(K/(K+1));
    %eigenvalues of H^{H}*Q*H, where Q is a scaled identity matrix P*I
    eig_vec(:,ij)= P* abs(svd( h)).^2;
end

%nn=[60:10:200 220:20:500 550:50:1000];

for n=nn
    %information density under P
    S_n = sum( n*log(1 + eig_vec)) + n*min_tx_rx  - sum(ncx2rnd(2*n, 2*n ./ eig_vec).* eig_vec./ (1 + eig_vec))/2;

    S_n_sort=sort(S_n);

    tau = (error + 1/loop_h ): 1/loop_h : 2*error;

    gamma_n_array = S_n_sort(round(tau*loop_h));

    %lower-bounding beta_{error} using the standard inequality
    y_a = gamma_n_array/n - log(tau-error)/n;
    %optimizing over gamma
    y_min=min(y_a);

    rate =y_min./log(2);
    rate_c=[rate_c,rate];
end

%

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
% GUI Layout Toolbox
% Version 2.3.4 (R2018b) 31-January-2019
%
% Panels
%   uix.Panel          - arrange a single element inside a standard panel
%   uix.CardPanel      - show one element from a list
%   uix.BoxPanel       - arrange a single element inside a panel with a boxed title
%   uix.TabPanel       - arrange elements in a panel with tabs for selecting which is visible
%   uix.ScrollingPanel - arrange a single element inside a scrollable panel
%
% Boxes
%   uix.HBox           - arrange elements horizontally in a single row
%   uix.VBox           - arrange elements vertically in a single column
%   uix.HBoxFlex       - arrange elements horizontally with draggable dividers
%   uix.VBoxFlex       - arrange elements vertically with draggable dividers
%   uix.HButtonBox     - arrange buttons horizontally in a single row
%   uix.VButtonBox     - arrange buttons vertically in a single column
%
% Grids
%   uix.Grid           - arrange elements in a two-dimensional grid
%   uix.GridFlex       - arrange elements in a two-dimensional grid with draggable dividers
%
% Other
%   uix.Empty          - create an empty space
%   uix.tracking       - track anonymized usage data

%  Copyright 2009-2019 The MathWorks, Inc.
%  $Revision: 1774 $ $Date: 2018-12-06 06:03:23 +0000 (Thu, 06 Dec 2018) $

author: goGPS-Project


text:
function [fig] = plot_acceleration_offline(accel_history, ...
    time_history, accel_norm_max, fontsize, color)

% plot_acceleration_offline - Plot the min/avg/max acceleration of the
% agents

x0 = 10;
y0 = 10;
width = 600;
height = 150;

t_steps = length(accel_history(:,1));
min_accel = zeros(t_steps,1);
max_accel = zeros(t_steps,1);
avg_accel = zeros(t_steps,1);

for k = 1:t_steps
    Accel_k = accel_history(k,:);
    Accel_k = reshape(Accel_k,3,[]);
    Accel_norm_k = sqrt(sum(Accel_k.^2,1));
    min_accel(k) = min(Accel_norm_k);
    max_accel(k) = max(Accel_norm_k);
    avg_accel(k) = mean(Accel_norm_k);
end

fig = figure('Name','Offline swarm acceleration','NumberTitle','off');
hold on;
grid on;
box on;

% Plot envelope
err_bar(1,:,:) = [max_accel-avg_accel, avg_accel-min_accel];
if ~isempty(color)
    line_props.col = {color};
    mseb(time_history(1:t_steps)',avg_accel',err_bar,line_props);
else
    mseb(time_history(1:t_steps)',avg_accel',err_bar);
end

% Plot reference value
reference = yline(0,'--','LineWidth',1.5);
reference.Color = [0.25 0.25 0.25];
if ~isempty(accel_norm_max)
    threshold = yline(accel_norm_max,'-.','LineWidth',1.5);
    threshold.Color = [0.25 0.25 0.25];
end

xlabel('Time [s]','fontsize',fontsize);
ylabel('Acceleration [m/s^2]','fontsize',fontsize);
if ~isempty(accel_norm_max)
    h = legend('Average', 'Reference', 'Threshold','fontsize',fontsize);
else
    h = legend('Average', 'Reference','fontsize',fontsize);
end
set(h, 'location', 'northeastoutside');

set(fig,'units','pixels','position',[x0,y0,width,height]);
set(fig,'PaperPositionMode','auto');

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [error_x,error_y] = stokespost_q1q1_bc(aez,fezx,fezy,elerrorx,elerrory,xy,ev,ebound)
%stokespost_q1q1_bc   postprocesses Poisson error estimator
%   [error_x,error_y] = stokespost_q1q1_bc(aez,fezx,fezy,elerrorx,elerrory,xy,ev,ebound);
%   input
%          aez                  elementwise Poisson problem matrices
%          fezx,fezy            elementwise rhs vectors
%          elerrorx, elerrory   elementwise error estimate (without BC imposition)
%          xy                   vertex coordinate vector
%          ev                   element mapping matrix
%          ebound               element edge boundary matrix
%   output
%          error_x, error_y     component elementwise error estimate
%
%   calls function localbc_xy
%   IFISS function: DJS; 8 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      x=xy(:,1); y=xy(:,2);
      nel=length(ev(:,1));
      lev=[ev,ev(:,1)];
	  error_x=elerrorx;  error_y=elerrory;
%
% recompute contributions from elements with Dirichlet boundaries
      nbde=length(ebound(:,1));
      ebdy = zeros(nel,1);
      edge = zeros(nel,1);
% isolate boundary elements
      for el = 1:nbde
      ee = ebound(el,1);
      ebdy(ee) = ebdy(ee)+1; edge(ee)=ebound(el,2);
      end
%
% two edge elements
      k2=find(ebdy==2);
      nel2b=length(k2);
% loop over two edge elements
      for el = 1:nel2b
      el2e=k2(el);
      kk=find(ebound(:,1) == el2e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
	  ae=squeeze(aez(el2e,1:5,1:5));
      fex=fezx(el2e,:)'; fey=fezy(el2e,:)';
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el2e,:)); yl=y(lev(el2e,:));
	  [bae,fex,fey] = localbc_xy(ae,fex,fey,edges,xl,yl);
% solve local problems
      errx=bae\fex;  erry=bae\fey;
	  error_x(el2e,1) = errx'*ae*errx; error_y(el2e,1) = erry'*ae*erry;
	  end
% end of element loop
%
% one edge elements
      k1=find(ebdy==1);
      nel1b=length(k1);
% loop over one edge elements
      for el = 1:nel1b
      el1e=k1(el);
      kk=find(ebound(:,1) == el1e);
      edges=ebound(kk,2);
% set up original matrix and RHS vector
      fex=fezx(el1e,:)'; fey=fezy(el1e,:)';
	  ae=squeeze(aez(el1e,1:5,1:5));
% set up local coordinates and impose interpolated error as Dirichlet bc
      xl=x(lev(el1e,:)); yl=y(lev(el1e,:));
	  [bae,fex,fey] = localbc_xy(ae,fex,fey,edges,xl,yl);
% solve local problems
      errx=bae\fex;  erry=bae\fey;
	  error_x(el1e,1) = errx'*ae*errx; error_y(el1e,1) = erry'*ae*erry;

      end
% end of element loop
%
      err_x = sqrt(sum(error_x)); error_x = sqrt(error_x);
	  err_y = sqrt(sum(error_y)); error_y = sqrt(error_y);
      fprintf('estimated velocity error (in energy):  (%10.6e,%10.6e) \n',err_x,err_y)
 return

author: johannesgerer

repo: jburkardt-m


text:
function varargout=mod(varargin)
%MOD (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        if ~isa(varargin{2},'double')
            error('MOD is currently only supported for DOUBLE second argument');
        end

        x = varargin{1};
        y = varargin{2};
        % Some boring code for scalarization
        if prod(size(x))==1 & prod(size(y))>1
            x = x*ones(size(y));
        elseif prod(size(y))==1 & prod(size(x))>1
            y = y*ones(size(x));
        end
        if ~all(size(x) == size(y))
            error('Matrix dimensions must agree.');
        end
        dim = size(x);
        x = reshape(x,prod(dim),1);
        y = reshape(y,prod(dim),1);
        z = [];
        % Create one variable for each element
        for i = 1:length(x)
            xi = extsubsref(x,i);
            yi = extsubsref(y,i);
            inarg = {xi,yi};
            z = [z;yalmip('define',mfilename,inarg{:})];
        end
        z = reshape(z,dim);
        varargout{1} = z;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph

        % Description using epigraphs
        t = varargin{2};
        x = varargin{3};
        y = varargin{4};

        % t = mod(x,y), i.e. t = x - n*y, n = floor(x/y)
        n = intvar(1,1);
        F = (t == x - y*n)  + ((x/y) - 1 <= n <= (x/y));

        varargout{1} = F;
        varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','integer');
        varargout{3} = [x(:);y(:)];

    otherwise
        error('Strange type on first argument in SDPVAR/MOD');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [axis, theta] = rotation3dAxisAndAngle(mat)
%ROTATION3DAXISANDANGLE Determine axis and angle of a 3D rotation matrix.
%
%   [AXIS, ANGLE] = rotation3dAxisAndAngle(MAT)
%   Where MAT is a 4-by-4 matrix representing a rotation, computes the
%   rotation axis (containing the points that remain invariant under the
%   rotation), and the rotation angle around that axis.
%   AXIS has the format [DX DY DZ], constrained to unity, and ANGLE is the
%   rotation angle in radians.
%
%   Note: this method use eigen vector extraction. It would be more precise
%   to use quaternions, see:
%   http://www.mathworks.cn/matlabcentral/newsreader/view_thread/160945
%
%
%   Example
%     origin = [1 2 3];
%     direction = [4 5 6];
%     line = [origin direction];
%     angle = pi/3;
%     rot = createRotation3dLineAngle(line, angle);
%     [axis angle2] = rotation3dAxisAndAngle(rot);
%     angle2
%     angle2 =
%           1.0472
%
%   See also
%   transforms3d, vectors3d, angles3d, eulerAnglesToRotation3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-08-11,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% extract the linear part of the rotation matrix
A = mat(1:3, 1:3);

% extract eigen values and eigen vectors
[V, D] = eig(A - eye(3));

% we need the eigen vector corresponding to eigenvalue==1
[dummy, ind] = min(abs(diag(D)-1)); %#ok<ASGLU>

% extract corresponding eigen vector
vector = V(:, ind)';

% compute rotation angle
t = [A(3,2)-A(2,3) , A(1,3)-A(3,1) , A(2,1)-A(1,2)];
theta = atan2(dot(t, vector), trace(A)-1);

% If angle is negative, invert both angle and vector direction
if theta<0
    theta  = -theta;
    vector = -vector;
end

% try to get a point on the line
% seems to work, but not sure about stability
[V, D] = eig(mat-eye(4)); %#ok<ASGLU>
origin = V(1:3,4)'/V(4, 4);

% create line corresponding to rotation axis
axis = [origin vector];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ lam, mu ] = cycle_floyd ( f, x0 )

%*****************************************************************************80
%
%% CYCLE_FLOYD finds a cycle in an iterated mapping using Floyd's method.
%
%  Discussion:
%
%    Suppose we a repeatedly apply a function f(), starting with the argument
%    x0, then f(x0), f(f(x0)) and so on.  Suppose that the range of f is finite.
%    Then eventually the iteration must reach a cycle.  Once the cycle is reached,
%    succeeding values stay within that cycle.
%
%    Starting at x0, there is a "nearest element" of the cycle, which is
%    reached after MU applications of f.
%
%    Once the cycle is entered, the cycle has a length LAM, which is the number
%    of steps required to first return to a given value.
%
%    This function uses Floyd's method to determine the values of MU and LAM,
%    given F and X0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Knuth,
%    The Art of Computer Programming,
%    Volume 2, Seminumerical Algorithms,
%    Third Edition,
%    Addison Wesley, 1997,
%    ISBN: 0201896842,
%    LC: QA76.6.K64.
%
%  Parameters:
%
%    Input, integer F( integer I ), the name of the function
%    to be analyzed.
%
%    Input, integer X0, the starting point.
%
%    Output, integer LAM, the length of the cycle.
%
%    Output, integer MU, the index in the sequence starting
%    at X0, of the first appearance of an element of the cycle.
%
  tortoise = f ( x0 );
  hare = f ( tortoise );

  while ( tortoise ~= hare )
    tortoise = f ( tortoise );
    hare = f ( f ( hare ) );
  end

  mu = 0;
  tortoise = x0;

  while ( tortoise ~= hare )
    tortoise = f ( tortoise );
    hare = f ( hare );
    mu = mu + 1;
  end

  lam = 1
  hare = f ( tortoise )
  while ( tortoise ~= hare )
    hare = f ( hare );
    lam = lam + 1;
  end

  return
end


text:
function a = r8vec_sort_bubble_d ( n, a )

%*****************************************************************************80
%
%% R8VEC_SORT_BUBBLE_D descending sorts an R8VEC using bubble sort.
%
%  Discussion:
%
%    Bubble sort is simple to program, but inefficient.  It should not
%    be used for large arrays.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), an unsorted array.
%
%    Output, real A(N), the array has been sorted.
%
  for i = 1 : n-1
    for j = i+1 : n
      if ( a(i) < a(j) )
        t    = a(i);
        a(i) = a(j);
        a(j) = t;
      end
    end
  end

  return
end


text:
function [rx,re] = addToGraph(x,e,rx,re)

global Map

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

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

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function e=v_rotro2eu(m,r)
%V_ROTRO2EU converts a 3x3 rotation matrix into the corresponding euler angles
% inverse of v_roteu2ro
%     M        a string of n characters from the set determining the order of rotation axes
%              as listed below. Note that the control characters 'rdoOaA' may occur anywhere in the string:
%                'x','y','z'    rotate around the given axis by the corresponding angle
%                               given in e()
%                '1','2','3'    90 degree rotation around x,y or z axis; doesn't use a value from e()
%                '4','5','6'    180 degree rotation around x,y or z axis; doesn't use a value from e()
%                '7','8','9'    270 degree rotation around x,y or z axis; doesn't use a value from e()
%                'r','d'        all angles are given in radians or degrees  [default='r']
%                'R','D'        all angles are given in radians or degrees and are negated
%             'o','O','a','A'   selects whether to rotate the object or the coordinate axes and
%                               whether the rotation axes remain fixed in space for consecutive
%                               rotations (extrinsic) or else move with each rotation (intrinsic).
%                                  'o' = object-extrinsic [default]
%                                  'O' = object-intrinsic
%                                  'a' = axes-extrinsic
%                                  'A' = axes-intrinsic
%
%     R(3,3,...)   Input rotation matrix (or array of matrices)
%
% Outputs:
%
%     E(K,...)  K Euler angles in radians (or degrees if 'd' or 'D' specified) per quaternion where K is the number of 'xyz' characters in M.
%               A positive rotation is clockwise if looking along the +ve axis away from the origin or anti-clockwise if 'R' or 'D' is given.
%               The x, y, z axes form a right-handed triple.
%
% The string M specifies the seqeunce of axes about which the rotations are performed. There are 12
% possible 3-character sequences that avoid consecutive repetitions. These are 'Euler angles' if
% there is a repeated axis or 'Tait-Bryan angles' if not. Common choices are:
% (1) 'zxz' the most common Euler angle set
% (2) 'xyz' corresponds to 'roll, pitch, yaw' for an aeroplane heading in the x direction with y to
%     the right and z down. The intrinsic equivalent is 'Ozyx' corresponding to 'yaw, pitch, roll'.
% (3) 'z1z1z' involves 5 rotations, in which all the non-fixed rotations are around the z axis.
%
% The Euler angles are not, in general, unique. In particular:
%  (1) v_roteu2ro('zxz',[a b c]) = v_roteu2ro('zxz',[a+pi -b c+pi])
%  (2) v_roteu2ro('xyz',[a b c]) = v_roteu2ro('xyz',[a+pi pi-b c+pi])

%      Copyright (C) Mike Brookes 2007-2020
%      Version: $Id: v_rotro2eu.m 11260 2020-07-18 20:07:58Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent mch mvch rtr rtci rtsi scai w6 th6 x6 pmw
if isempty(mvch)
    mch=''; % cached rotation code string
    mvch=[0;0;1;1;0;0;1];
    rtr=[1 4 7 2 5 8 3 6 9]; % indices to transpose a vectorized 3x3 matrix
    rtci=[2 3 5 6 8 9; 3 1 6 4 9 7; 1 2 4 5 7 8]';
    rtsi=[3 2 6 5 9 8; 1 3 4 6 7 9; 2 1 5 4 8 7]';
    scai=[0 0 0 1; 0 0 0 2; 0 0 0 3; 1 -1 0 1; 1 -1 0 2; 1 -1 0 3; 0 0 -1 1; 0 0 -1 2; 0 0 -1 3; -1 1 0 1; -1 1 0 2; -1 1 0 3]'; % [sin; -sin; cos; xyz] for fixed rotations
    w6=ones(6,1); %
    th6=3*w6;
    x6=[2 1 2 1 2 1]'; % Index for sin components
    pmw=[1; -1];
end

if ~nargout
    v_rotro2qr(r);
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert the m string
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ischar(m) && strcmp(m,mch) % check if the rotation code string is cached
        mv=mvch;
    else
        mv=v_roteucode(m); % else decode the string
        mch=m; % and save the result in the cache for next time.
        mvch=mv;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now calculate euler angles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nm=size(mv,2)-1; % number of rotation codes
    sz=size(r);
    r=reshape(r,9,[]); % vectorize the rotation matrices
    nr=size(r,2);                       % number of rotation matrices
    if mv(end)<0
        r=r(rtr,:);                     % transpose rotation matrix
    end
    e=zeros(mv(2,end),nr);              % initialize array of euler angles
    ef=mv(end-3);
    for i=nm:-1:1                       % process rotations in reverse order
        mvi=mv(:,i);
        mi=mvi(1);
        if mi<=3                        % rotation around x,y or z
            if mvi(6)~=0               % skip if this rotation is redundant
                [si,ci,ri,ti]=v_atan2sc(mvi(7)*r(mvi(4),:),mvi(7)*r(mvi(5),:));
                e(mv(2,i),:)=ti*mvi(6)/ef;  % save the euler angle
                si=mvi(6)*pmw*si;          % make -si available and correct sign
                r(rtci(:,mi),:)=ci(w6,:).*r(rtci(:,mi),:)-si(x6,:).*r(rtsi(:,mi),:); % apply reverse rotation
            end
        else % fixed rotation
            ai=scai(4,mi);                      % axis of rotation: 1, 2 or 3
            r(rtci(:,ai),:)=scai(th6,mi).*r(rtci(:,ai),:)-scai(x6,mi).*r(rtsi(:,ai),:); % apply reverse rotation
        end
    end
    if nr>1                                     % if there was >1 input matrix
        e=reshape(e,[size(e,1) sz(3:end)]);     % restore the shape of the Euler angle array
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function ComparingSymmetry

incPhi = pi/180;

fileName = 'StressSymmetryTraction';
sPnormT = computeExperiment(incPhi,fileName);

fileName = 'StressSymmetryCompression';
sPnormC = computeExperiment(-incPhi,fileName);

end

function sPnorm = computeExperiment(incPhi,fileName)

txi = pi/2 - 0.2;%1083;
rho = 0.9;
q = 4;

phi = 0+incPhi;
pNorm = 'max';
print = false;
hMesh = 0.01;
hasToCaptureImage = true;

mx = SuperEllipseParamsRelator.mx(txi,rho,q);
my = SuperEllipseParamsRelator.my(txi,rho,q);

s.mx       = mx;
s.my       = my;
s.q        = q;
s.phi      = phi;
s.pNorm    = pNorm;
s.print    = print;
s.hMesh    = hMesh;
s.fileName = fileName;
s.hasToCaptureImage = false;
sN = StressNormSuperEllipseComputer(s);
sPnorm = sN.compute();
sN.printStress();

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [xVi, mask] = spm_est_non_sphericity(SPM)
% Non-sphericity estimation using ReML
% FORMAT [xVi, mask] = spm_est_non_sphericity(SPM)
%
% Required fields of SPM structure (see spm_spm):
% SPM.xY.VY  - nScan x 1 struct array of file handles
% SPM.xX     - structure containing design matrix information
% SPM.xX.W   - optional whitening/weighting matrix
% SPM.xVi    - structure describing intrinsic non-sphericity
% SPM.xM     - structure containing masking information
%
% Return xVi from SPM.xVi with extra fields:
% xVi.V      - estimated non-sphericity, trace(V) = rank(V)
% xVi.h      - hyperparameters  xVi.V = xVi.h(1)*xVi.Vi{1} + ...
% xVi.Cy     - spatially whitened <Y*Y'> (used by ReML to estimate h)
%
% mask       - logical array of voxels within analysis mask
%__________________________________________________________________________
%
% In a first pass, voxels over which non-sphericity will be estimated are
% selected using an 'effects of interest' F-contrast (can be specified in
% SPM.xVi.Fcontrast) and critical threshold taken from SPM defaults
% stats.<modality>.UFp.
% The sample covariance matrix (xVi.Cy) is then estimated by pooling over
% these voxels, assuming V is constant over them.
% Finally, SPM will invoke ReML to estimate hyperparameters (xVi.h) of an
% array of non-sphericity components (xVi.Vi), providing a high precise
% estimate of the non-sphericity matrix (xVi.V).
%__________________________________________________________________________
% Copyright (C) 1994-2019 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Guillaume Flandin
% $Id: spm_est_non_sphericity.m 7577 2019-04-24 08:59:56Z guillaume $

SVNid = '$Rev: 7577 $';

%-Say hello
%--------------------------------------------------------------------------
spm('FnBanner',mfilename,SVNid);

%-Get data, design, mask and variance components details
%--------------------------------------------------------------------------
VY           = SPM.xY.VY;
M            = VY(1).mat;
DIM          = VY(1).dim;
YNaNrep      = spm_type(VY(1).dt(1),'nanrep');

xX           = SPM.xX;
nScan        = size(xX.X,1);
if ~isfield(xX,'W')
    xX.W     = speye(nScan,nScan);
end

xM           = SPM.xM;

xVi          = SPM.xVi;

%-Compute Hsqr and F-threshold under i.i.d.
%--------------------------------------------------------------------------
xX.xKXs      = spm_sp('Set',spm_filter(xX.K,xX.W*xX.X));
xX.xKXs.X    = full(xX.xKXs.X);
xX.pKX       = spm_sp('x-',xX.xKXs);

if isfield(xVi,'Fcontrast')
    Fcname   = 'User-specified contrast';
    xCon     = spm_FcUtil('Set',Fcname,'F','c',xVi.Fcontrast,xX.xKXs);
else
    Fcname   = 'effects of interest';
    iX0      = [xX.iB xX.iG];
    xCon     = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
end

if ~isempty(xCon(1).c)
    X1o      = spm_FcUtil('X1o', xCon(1),xX.xKXs);
    Hsqr     = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
    trMV     = spm_SpUtil('trMV',X1o);
else
    % Force all voxels to enter non-sphericity
    trMV     = 1;
    Hsqr     = Inf;
end
trRV         = spm_SpUtil('trRV',xX.xKXs);

%-Threshold for voxels entering non-sphericity estimates
%--------------------------------------------------------------------------
try
    modality = lower(spm_get_defaults('modality'));
    UFp      = spm_get_defaults(['stats.' modality '.ufp']);
catch
    UFp      = 0.001;
end
xVi.UFp      = UFp;
UF           = spm_invFcdf(1 - UFp,[trMV,trRV]);

%==========================================================================
%-         P O O L E D   V A R I A N C E   E S T I M A T I O N
%==========================================================================

%-Get explicit mask(s)
%==========================================================================
mask = true(DIM);
for i = 1:numel(xM.VM)
    if ~(isfield(SPM,'xVol') && isfield(SPM.xVol,'G'))
        %-Assume it fits entirely in memory
        C = spm_bsplinc(xM.VM(i), [0 0 0 0 0 0]');
        v = true(DIM);
        [x1,x2] = ndgrid(1:DIM(1),1:DIM(2));
        for x3 = 1:DIM(3)
            M2  = inv(M\xM.VM(i).mat);
            y1 = M2(1,1)*x1+M2(1,2)*x2+(M2(1,3)*x3+M2(1,4));
            y2 = M2(2,1)*x1+M2(2,2)*x2+(M2(2,3)*x3+M2(2,4));
            y3 = M2(3,1)*x1+M2(3,2)*x2+(M2(3,3)*x3+M2(3,4));
            v(:,:,x3) = spm_bsplins(C, y1,y2,y3, [0 0 0 0 0 0]') > 0;
        end
        mask = mask & v;
        clear C v x1 x2 x3 M2 y1 y2 y3
    else
        if spm_mesh_detect(xM.VM(i))
            v = full(xM.VM(i).private.cdata) > 0;
        else
            v = spm_mesh_project(gifti(SPM.xVol.G), xM.VM(i)) > 0;
        end
        mask = mask & v(:);
        clear v
    end
end

Cy        = 0;                                  %-<Y*Y'> spatially whitened
Cm        = mask;                               %-mask of selected voxels

%-Split data into chunks
%==========================================================================
chunksize = floor(spm_get_defaults('stats.maxmem') / 8 / nScan);
nbchunks  = ceil(prod(DIM) / chunksize);
chunks    = min(cumsum([1 repmat(chunksize,1,nbchunks)]),prod(DIM)+1);

spm_progress_bar('Init',nbchunks,'Hyperparameter estimation','Chunks');

for i=1:nbchunks
    chunk = chunks(i):chunks(i+1)-1;

    %-Report progress
    %======================================================================
    if i > 1, fprintf(repmat(sprintf('\b'),1,72)); end                  %-#
    fprintf('%-40s: %30s', sprintf('Chunk %3d/%-3d',i,nbchunks),...
                           '...processing');                            %-#

    %-Get data & construct analysis mask
    %----------------------------------------------------------------------
    Y       = zeros(nScan,numel(chunk));
    cmask   = mask(chunk);
    for j=1:nScan
        if ~any(cmask), break, end                 %-Break if empty mask

        Y(j,cmask) = spm_data_read(VY(j),chunk(cmask));%-Read chunk of data

        cmask(cmask) = Y(j,cmask) > xM.TH(j);      %-Threshold (& NaN) mask
        if xM.I && ~YNaNrep && xM.TH(j) < 0        %-Use implicit mask
            cmask(cmask) = abs(Y(j,cmask)) > eps;
        end
    end
    cmask(cmask) = any(diff(Y(:,cmask),1));        %-Mask constant data
    mask(chunk)  = cmask;
    Cm(chunk)    = cmask;
    if ~any(cmask), continue, end
    Y       = Y(:,cmask);                          %-Data within mask

    %-Remove filter confounds
    %----------------------------------------------------------------------
    KWY      = spm_filter(xX.K,xX.W*Y);

    %-Ordinary Least Squares estimation
    %----------------------------------------------------------------------
    beta    = xX.pKX*KWY;                          %-Parameter estimates
    if any(cmask)
        res = spm_sp('r',xX.xKXs,KWY);             %-Residuals
    else
        res = zeros(nScan,0);
    end
    ResSS   = sum(res.^2);                         %-Residual SSQ
    clear res

    %-F-threshold & accumulate spatially whitened Y*Y'
    %----------------------------------------------------------------------
    j       = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV;
    Cm(chunk(cmask)) = j;
    q       = nnz(j);
    if q
        q   = spdiags(sqrt(trRV./ResSS(j)'),0,q,q);
        Y   = Y(:,j)*q;
        Cy  = Cy + Y*Y';
    end

    %-Report progress
    %======================================================================
    fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...done');             %-#
    spm_progress_bar('Set',i);
end

fprintf('\n');                                                          %-#
spm_progress_bar('Clear');

s = nnz(Cm);                                    %-Number of selected voxels
if ~s
    error('Please check your data: There are no significant voxels.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Estimate sample covariance matrix Y*Y' from voxels in Cm
% Cy        = 0;
% nbchunks  = ceil(s / chunksize);
% chunks    = min(cumsum([1 repmat(chunksize,1,nbchunks)]),s+1);
% vx        = find(Cm);
% for i=1:nbchunks
%     chunk = chunks(i):chunks(i+1)-1;
%     n     = numel(chunk);
%     Y     = zeros(nScan,n);
%     for j=1:nScan
%         Y(j,:) = spm_data_read(VY(j),vx(chunk));
%     end
%     % store ResSS above to prevent recomputing it here:
%     %--------------------------------------------------
%     KY    = spm_filter(xX.K,Y);
%     ResSS = spm_sp('r',xX.xKXs,KY);
%     ResSS = sum(ResSS.^2);
%     %--------------------------------------------------
%     q     = spdiags(sqrt(trRV./ResSS'),0,n,n);
%     Y     = Y*q;
%     Cy    = Cy + Y*Y';
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cy = Cy / s;                                    %-Sample covariance matrix

%==========================================================================
%-                     R e M L   E S T I M A T I O N
%==========================================================================

%-ReML estimate of residual correlations through hyperparameters (h)
%--------------------------------------------------------------------------
str    = sprintf('Temporal non-sphericity (%d voxels)',s);
fprintf('%-40s: %30s\n',str,'...ReML estimation');                      %-#

%-ReML for separable designs and covariance components
%--------------------------------------------------------------------------
if isstruct(xX.K)
    m     = length(xVi.Vi);
    h     = zeros(m,1);
    V     = sparse(nScan,nScan);
    for i = 1:length(xX.K)

        % extract blocks from bases
        %------------------------------------------------------------------
        q     = xX.K(i).row;
        p     = [];
        Qp    = {};
        for j = 1:m
            if nnz(xVi.Vi{j}(q,q))
                Qp{end + 1} = xVi.Vi{j}(q,q);
                p           = [p j];
            end
        end

        % design space for ReML (with confounds in filter)
        %------------------------------------------------------------------
        Xp     = xX.X(q,:);
        try
            Xp = [Xp xX.K(i).X0];
        end

        % ReML
        %------------------------------------------------------------------
        fprintf('%-30s\n',sprintf('  ReML Block %i',i));
        [Vp,hp] = spm_reml(Cy(q,q),Xp,Qp);
        V(q,q)  = V(q,q) + Vp;
        h(p)    = hp;
    end
else
    [V,h] = spm_reml(Cy,xX.X,xVi.Vi);
end

%-Normalize non-sphericity and save hyperparameters
%--------------------------------------------------------------------------
V      = V*nScan/trace(V);
xVi.h  = h;
xVi.V  = V;                                     % Save non-sphericity xVi.V
xVi.Cy = Cy;                                    % spatially whitened <Y*Y'>

author: spm


text:
function x = hypersphere_to_cartesian ( m, n, c, r, theta )

%*****************************************************************************80
%
%% HYPERSPHERE_TO_CARTESIAN: hypersphere to Cartesian coordinate transform.
%
%  Discussion:
%
%    We allow the trivial case M = 1; in that case alone, the value R
%    must be assumed to have a sign.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%    1 <= M.
%
%    Input, integer N, the number of points to transform.
%
%    Input, real C(M,1), the center of the hypersphere.
%
%    Input, real R(N,1), the radius of the points on the hypersphere.
%    Except for the trivial case M = 1, R is assumed nonnegative.
%
%    Input, real THETA(M-1,N), the coordinate angles of the points,
%    measured in radians.
%
%    Output, real X(M,N), the Cartesian coordinates of the points.
%
  x = zeros ( m, n );
%
%  Make R a row vector.
%
  r = ( r(:) )';
%
%  Handle special case of M = 1.
%
  if ( m == 1 )
    x(1:1,1:n) = repmat ( c(1:1,1:1), 1, n ) ...
               + repmat ( r(1:1,1:n), 1, 1 );
    return
  end

  x(1:m,1:n) = repmat ( r(1:1,1:n), m, 1 );

  for i1 = 1 : m - 1
    x(i1,1:n) = x(i1,1:n) .* cos ( theta(i1,1:n) );
    for i2 = i1 + 1 : m
      x(i2,1:n) = x(i2,1:n) .* sin ( theta(i1,1:n) );
    end
  end
%
%  Add the center.
%
  x(1:m,1:n) = x(1:m,1:n) + repmat ( c(1:m,1:1), 1, n );

  return
end


text:
function fig = RlcResponse()
% This is the machine-generated representation of a Handle Graphics object
% and its children.  Note that handle values may change when these objects
% are re-created. This may cause problems with any callbacks written to
% depend on the value of the handle at the time the object was saved.
% This problem is solved by saving the output as a FIG-file.
%
% To reopen this object, just type the name of the M-file at the MATLAB
% prompt. The M-file and its associated MAT-file must be on your path.
%
% NOTE: certain newer features in MATLAB may not have been saved in this
% M-file due to limitations of this format, which has been superseded by
% FIG-files.  Figures which have been annotated using the plot editor tools
% are incompatible with the M-file/MAT-file format, and should be saved as
% FIG-files.
warning('off','MATLAB:dispatcher:InexactMatch')
load RlcResponse

h0 = figure('Units','normalized', ...
	'Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
	'FileName','C:\MATLABR11\Circuits\RlcResponse.m', ...
	'HitTest','off', ...
	'PaperPosition',[18 180 576.0000000000001 432.0000000000002], ...
	'PaperUnits','points', ...
	'Position',[0.09 0.04952380952380953 0.8 0.8], ...
	'Tag','Fig1', ...
	'ToolBar','none');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 1 0.501960784313725], ...
	'ListboxTop',0, ...
	'Position',[0.3284493284493285 0.04715447154471544 0.4932844932844932 0.09593495934959348], ...
	'Style','frame', ...
	'Tag','Frame1');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 0 0], ...
	'Callback','rlcparallel', ...
	'FontSize',10, ...
	'FontWeight','bold', ...
	'ForegroundColor',[0 0 0.627450980392157], ...
	'ListboxTop',0, ...
	'Position',[0.3401785714285714 0.05952380952380953 0.1794642857142857 0.06785714285714287], ...
	'String','Parallel Circuit', ...
	'Tag','Pushbutton1', ...
	'TooltipString','Parallel circuit response');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 0 0], ...
	'Callback','rlcseries', ...
	'FontSize',10, ...
	'FontWeight','bold', ...
	'ForegroundColor',[0 0 0.627450980392157], ...
	'ListboxTop',0, ...
	'Position',[0.5311355311355309 0.06025293586269195 0.1790123456790123 0.06811201445347785], ...
	'String','Series Circuit', ...
	'Tag','Pushbutton2', ...
	'TooltipString','Series circuit response');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[0 0.501960784313725 1], ...
	'FontName','Script MT Bold', ...
	'FontSize',14, ...
	'ForegroundColor',[1 0 0], ...
	'ListboxTop',0, ...
	'Position',mat1, ...
	'String','    Hadi Saadat', ...
	'Style','text', ...
	'Tag','StaticText7', ...
	'TooltipString','www.msoe.edu/~saadat');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 0 0], ...
	'FontName','Agaramond', ...
	'FontSize',16, ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Max',0, ...
	'Position',[0.06227106227106227 0.04552845528455284 0.0525030525030525 0.04715447154471544], ...
	'String','OE', ...
	'Style','text', ...
	'Tag','StaticText5', ...
	'TooltipString','Milwaukee School of Engineering');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 0 0], ...
	'FontName','Agaramond', ...
	'FontSize',16, ...
	'FontWeight','demi', ...
	'ForegroundColor',[1 1 1], ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[0.06227106227106227 0.08943089430894309 0.0525030525030525 0.05203252032520326], ...
	'String','MS', ...
	'Style','text', ...
	'Tag','StaticText5', ...
	'TooltipString','Milwuakee School of Engineering');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[1 0 0], ...
	'Callback','close(gcbf)', ...
	'FontSize',10, ...
	'FontWeight','bold', ...
	'ForegroundColor',[0 0 0.627450980392157], ...
	'ListboxTop',0, ...
	'Position',[0.7219508886175551 0.06025293586269195 0.08262108262108259 0.06811201445347785], ...
	'String','Exit', ...
	'Tag','Pushbutton3', ...
	'TooltipString','Exit');
h1 = uicontrol('Parent',h0, ...
	'Units','normalized', ...
	'BackgroundColor',[0 0.501960784313725 1], ...
	'FontSize',9, ...
	'FontWeight','bold', ...
	'ForegroundColor',[1 0 0], ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[0.1233211233211233 0.09105691056910568 0.0293040293040293 0.04390243902439023], ...
	'String','(c)', ...
	'Style','text', ...
	'Tag','StaticText6', ...
	'TooltipString','Copyright 1999');
h1 = axes('Parent',h0, ...
	'Box','on', ...
	'CameraUpVector',[0 1 0], ...
	'CameraUpVectorMode','manual', ...
	'Color',[1 1 1], ...
	'ColorOrder',mat2, ...
	'HitTest','off', ...
	'Position',[0.0702 0.1385 0.876 0.7662], ...
	'Tag','Axes1', ...
	'Visible','off', ...
	'XColor',[0 0 0], ...
	'XLim',[0 4], ...
	'XLimMode','manual', ...
	'YColor',[0 0 0], ...
	'YLim',[-3.5 1.5], ...
	'YLimMode','manual', ...
	'ZColor',[0 0 0]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'EraseMode','none', ...
	'Tag','Axes1Line18', ...
	'XData',[0.2 0.2 1.1 1.1], ...
	'YData',[0.6 1 1 0.55]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'EraseMode','none', ...
	'Tag','Axes1Line17', ...
	'XData',[0.2 0.2 1.1 1.1], ...
	'YData',[0.4 0 0 0.45]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line16', ...
	'XData',mat3, ...
	'YData',mat4);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line15', ...
	'XData',mat5, ...
	'YData',mat6);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line14', ...
	'XData',mat7, ...
	'YData',mat8);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line13', ...
	'XData',[0.8 0.8], ...
	'YData',[0 0.3]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line12', ...
	'XData',[0.8 0.8], ...
	'YData',[0.7 1]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line11', ...
	'XData',[1.05 1.15], ...
	'YData',[0.45 0.45]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line10', ...
	'XData',[1.05 1.15], ...
	'YData',[0.55 0.55]);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[1.995918367346939 1.55443234836703 17.32050807568877], ...
	'Tag','Axes1Text56', ...
	'VerticalAlignment','bottom');
set(get(h2,'Parent'),'Title',h2);
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.025 0.5 0], ...
	'String','I_s', ...
	'Tag','Axes1Text55');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.375 0.5 0], ...
	'String','R', ...
	'Tag','Axes1Text54');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.75 0.5 0], ...
	'String','L', ...
	'Tag','Axes1Text53');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[1.015 0.36 0], ...
	'String','C', ...
	'Tag','Axes1Text52');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.425 0.825 0], ...
	'String','\downarrow', ...
	'Tag','Axes1Text51');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.725 0.825 0], ...
	'String','\downarrow', ...
	'Tag','Axes1Text50');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[1.025 0.825 0], ...
	'String','\downarrow', ...
	'Tag','Axes1Text49');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.175 0.485 0], ...
	'String','\uparrow', ...
	'Tag','Axes1Text48');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'Position',[0.725 0.65 0], ...
	'String','\downarrow', ...
	'Tag','Axes1Text47');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[0.65 0.65 0], ...
	'String','I_0', ...
	'Tag','Axes1Text46');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[0.925 0.5 0], ...
	'String','V_0', ...
	'Tag','Axes1Text45');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[0.925 0.725 0], ...
	'String','+', ...
	'Tag','Axes1Text44');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[0.925 0.35 0], ...
	'String','-', ...
	'Tag','Axes1Text43');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[1.2 0.5 0], ...
	'String','v', ...
	'Tag','Axes1Text42');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[1.15 0.9 0], ...
	'String','+', ...
	'Tag','Axes1Text41');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[1.15 0.1 0], ...
	'String','-', ...
	'Tag','Axes1Text40');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.35 0.9 0], ...
	'String','i_R', ...
	'Tag','Axes1Text39');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.65 0.9 0], ...
	'String','i_L', ...
	'Tag','Axes1Text38');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0.95 0.9 0], ...
	'String','i_C', ...
	'Tag','Axes1Text37');
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'EraseMode','none', ...
	'Tag','Axes1Line9', ...
	'XData',mat9, ...
	'YData',mat10);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'EraseMode','none', ...
	'Tag','Axes1Line8', ...
	'XData',[2.2 2.2 3.1 3.1], ...
	'YData',[0.4 0 0 0.45]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line7', ...
	'XData',mat11, ...
	'YData',mat12);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line6', ...
	'XData',mat13, ...
	'YData',mat14);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line5', ...
	'XData',mat15, ...
	'YData',mat16);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line4', ...
	'XData',mat17, ...
	'YData',mat18);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line3', ...
	'XData',[3 3.1 3.1], ...
	'YData',[1 1 0.55]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line2', ...
	'XData',[3.05 3.15], ...
	'YData',[0.45 0.45]);
h2 = line('Parent',h1, ...
	'Color',[0 0 1], ...
	'Tag','Axes1Line1', ...
	'XData',[3.05 3.15], ...
	'YData',[0.55 0.55]);
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2 0.5 0], ...
	'String','V_s', ...
	'Tag','Axes1Text36');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.375 1.15 0], ...
	'String','R', ...
	'Tag','Axes1Text35');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.75 1.15 0], ...
	'String','L', ...
	'Tag','Axes1Text34');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[3.025 0.35 0], ...
	'String','C', ...
	'Tag','Axes1Text33');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.6 1.1 0], ...
	'String','\rightarrow', ...
	'Tag','Axes1Text32');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.175 0.5 0], ...
	'String','\pm', ...
	'Tag','Axes1Text31');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'Position',[2.875 1.125 0], ...
	'String','\rightarrow', ...
	'Tag','Axes1Text30');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[2.85 1.275 0], ...
	'String','I_0', ...
	'Tag','Axes1Text29');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[2.9 0.5 0], ...
	'String','V_0', ...
	'Tag','Axes1Text28');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[2.95 0.725 0], ...
	'String','+', ...
	'Tag','Axes1Text27');
h2 = text('Parent',h1, ...
	'Color',[1 0 0], ...
	'FontSize',9, ...
	'Position',[2.95 0.3 0], ...
	'String','-', ...
	'Tag','Axes1Text26');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[3.175 0.5 0], ...
	'String','v_C', ...
	'Tag','Axes1Text25');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[3.15 0.9 0], ...
	'String','+', ...
	'Tag','Axes1Text24');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[3.15 0.1 0], ...
	'String','-', ...
	'Tag','Axes1Text23');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.425 0.85 0], ...
	'String','v_R', ...
	'Tag','Axes1Text22');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.3 0.875 0], ...
	'String','+', ...
	'Tag','Axes1Text21');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.55 0.875 0], ...
	'String','-', ...
	'Tag','Axes1Text20');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.625 1.2 0], ...
	'String','i', ...
	'Tag','Axes1Text19');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.825 0.85 0], ...
	'String','v_L', ...
	'Tag','Axes1Text18');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.7 0.875 0], ...
	'String','+', ...
	'Tag','Axes1Text17');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[2.95 0.875 0], ...
	'String','-', ...
	'Tag','Axes1Text16');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',14, ...
	'Position',[0.5 1.95 0], ...
	'String',mat19, ...
	'Tag','Axes1Text15');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -0.25 0], ...
	'String',mat20, ...
	'Tag','Axes1Text14');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 1.5 0], ...
	'String',mat21, ...
	'Tag','Axes1Text13');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -0.525 0], ...
	'String',mat22, ...
	'Tag','Axes1Text12');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -0.775 0], ...
	'String',mat23, ...
	'Tag','Axes1Text11');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -1 0], ...
	'String',mat24, ...
	'Tag','Axes1Text10');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -1.25 0], ...
	'String',mat25, ...
	'Tag','Axes1Text9');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -1.65 0], ...
	'String',mat26, ...
	'Tag','Axes1Text8');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -1.95 0], ...
	'String',mat27, ...
	'Tag','Axes1Text7');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -2.175 0], ...
	'String',mat28, ...
	'Tag','Axes1Text6');
h2 = text('Parent',h1, ...
	'Color',[0 0 0.502], ...
	'FontSize',9, ...
	'Position',[0 -2.35 0], ...
	'String',mat29, ...
	'Tag','Axes1Text5');
h2 = text('Parent',h1, ...
	'Color',[1 0 0.502], ...
	'FontName','Roman', ...
	'FontSize',11, ...
	'Position',[0 -2.95 0], ...
	'String',mat30, ...
	'Tag','Axes1Text4');
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[1.995918367346939 -3.686625194401245 17.32050807568877], ...
	'Tag','Axes1Text3', ...
	'VerticalAlignment','cap', ...
	'Visible','off');
set(get(h2,'Parent'),'XLabel',h2);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','center', ...
	'Position',[-0.1346938775510204 -1.011664074650078 17.32050807568877], ...
	'Rotation',90, ...
	'Tag','Axes1Text2', ...
	'VerticalAlignment','baseline', ...
	'Visible','off');
set(get(h2,'Parent'),'YLabel',h2);
h2 = text('Parent',h1, ...
	'Color',[0 0 0], ...
	'HandleVisibility','off', ...
	'HorizontalAlignment','right', ...
	'Position',[-0.326530612244898 2.114307931570762 17.32050807568877], ...
	'Tag','Axes1Text1', ...
	'Visible','off');
set(get(h2,'Parent'),'ZLabel',h2);
if nargout > 0, fig = h0; end

author: Sable

repo: mcbench-benchmarks


text:
% The COBRAToolbox: testDynamicModelFieldModification.m
%
% Purpose:
%     - Test functionality in removeFieldEntriesForType, extendModelFields,
%     getModelFieldsForType
%
% Author:
%     - Original file: Thomas Pfau

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testDynamicModelFieldModification.m'));
cd(fileDir);

%Lets create an empty model.
model = createModel();

%Lets add a gene:
model.genes{end+1,1} = 'gene1';

%extend all relevant fields (the original size was 0.
%get all Model fields for genes
fprintf('Testing getModelFieldsForType and extendModelFieldsForType ...\n');
[matchingGeneFields,dimensions] = getModelFieldsForType(model,'genes','fieldSize',0);
assert(dimensions == 2);
assert(isequal(matchingGeneFields,{'rxnGeneMat'}));
model = extendModelFieldsForType(model,'genes');
%This should result in an empty rxnGEneMatrix of size 0x1
assert(all(size(model.rxnGeneMat) == [0,1]));

%Recreate the model.
model = createModel();

%lets add a couple of reactions
reactions = {'R1',{'A','B'},[-1 1],0,1000;...
             'R2',{'B','C'},[-1 2],0,1000;...
             'R3',{'C','D'},[-1 1],0,1000;...
             'R4',{'E','F'},[-1 1],0,1000;...
             'R5',{'F','G'},[-1 1],0,1000;...
             'R6',{'G','E'},[-1 1],0,1000;...
             'R7',{'G','D'},[-1 1],0,1000;...
             'R8',{'H','C'},[-1 1],0,1000};
%Add Reactions
for i = 1:size(reactions,1)
    %All reactions are irreversible
    model = addReaction(model,reactions{i,1},'metaboliteList',reactions{i,2},...
          'stoichCoeffList',reactions{i,3},'lowerBound',reactions{i,4},'upperBound',reactions{i,5},...
          'printLevel',-1);
end
%Now, the size of model.mets is exactly the same as model.rxns (which could
%be a problematic situation).
%so, lets see, which fields we get for mets:
[matchingMetFields,dimensions] = getModelFieldsForType(model,'mets');
%This should be:
assert(isempty(setxor(matchingMetFields,{'mets','b','metNames','csense','S'})) && numel(matchingMetFields) == 5); %No duplicates!
%And S should have dimension 1, and only 1.
assert(all(dimensions(ismember(matchingMetFields,'S')) == 1)); %all to make sure that duplicate appearences are covered.

%Lets also check this for reactions
[matchingRxnFields,dimensions] = getModelFieldsForType(model,'rxns');
assert(isempty(setxor(matchingRxnFields,{'rxns','lb','ub','c','rxnGeneMat','rules','grRules','S'})) && numel(matchingRxnFields) == 8); %No duplicates!
%And S should have dimension 1, and only 1.
assert(all(dimensions(ismember(matchingRxnFields,'S')) == 2)); %all to make sure that duplicate appearences are covered.

%Finally for genes:
[matchingGeneFields,dimensions] = getModelFieldsForType(model,'genes');
assert(isempty(setxor(matchingGeneFields,{'genes','rxnGeneMat'})) && numel(matchingGeneFields) == 2); %No duplicates!
%And S should have dimension 1, and only 1.
assert(all(dimensions(ismember(matchingGeneFields,'rxnGeneMat')) == 2)); %all to make sure that duplicate appearences are covered.

%Now, add a metabolite
model.mets{end+1,1} = 'NewMet';
%and see, which fields we get when asking for all met fields one smaller
%than mets
[matchingMetFields,dimensions] = getModelFieldsForType(model,'mets','fieldSize',length(model.mets)-1);
%This should be all except mets:
assert(isempty(setxor(matchingMetFields,{'b','metNames','csense','S'})) && numel(matchingMetFields) == 4); %No duplicates!
%And S should have dimension 1, and only 1.
assert(all(dimensions(ismember(matchingMetFields,'S')) == 1)); %all to make sure that duplicate appearences are covered.

%now, lets add an unknown field
model.something = zeros(8,1); %This field will not be returned, as there are other fields (rxns) of the same size, so its unclear, what it belongs to.
[matchingMetFields,dimensions] = getModelFieldsForType(model,'mets','fieldSize',length(model.mets)-1);
assert(~any(ismember(matchingMetFields,'something')));

%If this would be size 9, it should be returned if all equal sized fields
%are requested:
model.something = zeros(1,9); %This field will not be returned, as there are other fields (rxns) of the same size, so its unclear, what it belongs to.
[matchingMetFields,dimensions] = getModelFieldsForType(model,'mets');
assert(any(ismember(matchingMetFields,'something')));
assert(dimensions(ismember(matchingMetFields,'something')) == 2);

% As a size 8 field it woudl be returned as a rxn field:
model.something = zeros(8,1); %This field will not be returned, as there are other fields (rxns) of the same size, so its unclear, what it belongs to.
[matchingRxnFields,dimensions] = getModelFieldsForType(model,'rxns');
assert(any(ismember(matchingRxnFields,'something')));
assert(dimensions(ismember(matchingRxnFields,'something')) == 1);

%So, lets extend the model to update all of those metabolite derived fields
model = extendModelFieldsForType(model,'mets');
[matchingMetFields,dimensions] = getModelFieldsForType(model,'mets');
assert(isempty(setxor(matchingMetFields,{'mets','b','metNames','csense','S'})) && numel(matchingMetFields) == 5);
assert(length(model.something) == 8) %This field was not modified, as it could not be clearly associated with mets.

fprintf('Testing removeFieldEntriesForType ...\n');
%Also, test the removeFieldEntries method now
%To be able to check whether the right fields are removed, we will create a
%random s matrix:
model.S = rand(size(model.S));
removedMets = [3 4 5 6];
model2 = removeFieldEntriesForType(model,removedMets,'mets', length(model.mets));
reducedSMatrix = model.S;
reducedSMatrix(removedMets,:) = [];
assert(isequal(model2.S, reducedSMatrix));
%Assert that the right mets are removed and that the others are retained.
assert(isempty(intersect(model2.mets,model.mets(removedMets))) && isempty(setxor(setdiff(model.mets,model.mets(removedMets)),model2.mets)))

%Now, remove a metabolite
model = removeMetabolites(model,model.mets{end});
% and try this again (same sized rxns and mets)
model2 = removeFieldEntriesForType(model,removedMets,'mets', length(model.mets));
reducedSMatrix = model.S;
reducedSMatrix(removedMets,:) = [];
assert(isequal(model2.S, reducedSMatrix));
%Assert that the right mets are removed and that the others are retained.
assert(isempty(intersect(model2.mets,model.mets(removedMets))) && isempty(setxor(setdiff(model.mets,model.mets(removedMets)),model2.mets)))

%Also give it a try with reactions:
% and try this again (same sized rxns and mets)
removedRxns = [3 4 5];
%Also modify lb, ub and rules for this case
model.lb = rand(size(model.rxns));
model.ub = rand(size(model.rxns));
model = changeGeneAssociation(model,model.rxns{1},'A or B');
model = changeGeneAssociation(model,model.rxns{3},'B and C');
model = changeGeneAssociation(model,model.rxns{6},'D or E');
model2 = removeFieldEntriesForType(model,removedRxns,'rxns', length(model.rxns));
[rxnPres,rxnPos] = ismember(model2.rxns,model.rxns);
assert(isequal(model2.S, model.S(:,rxnPos(rxnPres))));
%Assert that the right mets are removed and that the others are retained.
assert(isempty(intersect(model2.rxns,model.rxns(removedRxns))) && isempty(setxor(setdiff(model.rxns,model.rxns(removedRxns)),model2.rxns)))
assert(all(cellfun(@(x,y) isequal(x,y), model2.rules, model.rules(rxnPos(rxnPres)))));
assert(all(model2.lb == model.lb(rxnPos(rxnPres))));
assert(all(model2.ub == model.ub(rxnPos(rxnPres))));

%And finally give ti a try for genes:
%First create the GR rules.
model = creategrRulesField(model);
removedGenes = [2 3];
model2 = removeFieldEntriesForType(model,removedGenes,'genes', length(model.genes));
assert(isempty(model2.rules{3})); % Both genes got removed.
assert(isequal(model2.grRules{1},'A')); % B got removed, so only A is retained.
assert(isequal(model2.rules{6},'x(2) | x(3)')); % The gene positions got changed.

%And test remobval of gene 4
model2 = removeFieldEntriesForType(model,4,'genes', length(model.genes));
assert(isequal(model2.rules{6},'x(4)')); % 4 got removed, and 5 renamed
[e] = parseBoolean(model2.grRules{6});
[e2] = parseBoolean('E');
assert(isequal(e,e2));
%Compare irrespective of actual format.

%Switch back
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [ymax, yu, xpos, xu] = max(s)

%tstoolbox/@signal/max
%   Syntax:
%     * [maximum, yunit, xpos, xunit] = max(s)
%
%   Give information about maximum of scalar signal s.
%
%   Example:
%disp('maximum of signal : ')
%disp(['y = ' num2str(m) ' ' label(yunit(s))]);
%disp(['x = ' num2str(xpos) ' ' label(a)]);
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

if ndim(s) ~= 1
	help(mfilename)
	return
end

[m, ind] = max(data(s));
a = getaxis(s, 1);
xpos = first(a) + delta(a) * (ind-1);	% determine position of maximum
% disp('maximum of signal : ')
% disp(['y = ' num2str(m) ' ' label(yunit(s))]);
% disp(['x = ' num2str(xpos) ' ' label(a)]);
xu = unit(a);
yu = yunit(s);
ymax = m;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [inF] = prepare_dcm(A,B,C,D,kernelGroup)
% prepares matrices and indices vectors for DCM (without HRF)
% function [options] = prepare_fullDCM(A,B,C,D,TR,microDT,homogeneous)
% IN:
%   - A: binary matrix indicating where the connections are
%   - B: cell-array of binary matrices of modulatory effects
%   - C: binary matrix of input-state coupling
%   - D: cell-array of binay matrices for gating effects
% OUT:
%   - inF: the optional input structure to @f_dcm4fmri

%check for empty matrices

if nargin<5
    kernelGroup=1;
end

nu = max([numel(B) size(C,2)]);
nx1 = max([size(A,1) cellfun(@(cel) size(cel,1),B) size(C,1) cellfun(@(cel) size(cel,1),D) ]);
nx2 = max([size(A,2) cellfun(@(cel) size(cel,2),B) size(C,1) cellfun(@(cel) size(cel,2),D) ]);

if isempty(A)
    A = zeros(nx1,nx2);
end
for i=1:nu
    if isempty(B) || numel(B)<i || isempty(B{i})
        B{i} = zeros(nx1,nx2);
    end
end
if isempty(C)
    C = zeros(nx1,nu);
end
for i=1:nx2
    if isempty(D) || numel(D)<i || isempty(D{i})
        D{i} = zeros(nx1,nx2);
    end
end

%prepare indices
[indA,indB,indC,indD,indself] = find_dcm(A,B,C,D,kernelGroup);

%save
inF.indself = indself;
inF.A = A;
inF.indA = indA;
inF.B = B;
inF.indB = indB;
inF.C = C;
inF.indC = indC;
inF.D = D;
inF.indD = indD;

[inF.dA,inF.dB,inF.dC,inF.dD] = get_dMatdvec(inF);

function [dA,dB,dC,dD] = get_dMatdvec(inF)
if ~isempty(inF.indA)
    dA = dMatdvec(inF.A);
else
    dA = [];
end
dB = cell(length(inF.B),1);
for i=1:length(inF.B)
    if ~isempty(inF.indB{i})
        dB{i} = dMatdvec(inF.B{i});
    else
        dB{i} = [];
    end
end
if ~isempty(inF.indC)
    dC = dMatdvec(inF.C);
else
    dC = [];
end
dD = cell(length(inF.D),1);
for i=1:length(inF.D)
    if ~isempty(inF.indD{i})
        dD{i} = dMatdvec(inF.D{i});
    else
        dD{i} = [];
    end
end

function C = dMatdvec(A)
A = ~~A;
ind = find(A~=0);
n = numel(A);
ni = length(ind);
C = zeros(n,ni);
for i=1:length(ind)
    C(ind(i),i) = 1;
end

function [indA,indB,indC,indD,indself] = find_dcm(A,B,C,D,kernelGroup)
ia = find(A~=0);
if ~isempty(ia)
    indA = 1:length(ia);
else
    indA = [];
end
ib = [];
indB = cell(length(B),1);
for i=1:length(B)
    tmp = find(B{i}~=0);
    if ~isempty(tmp)
        indB{i} = length(ia)+length(ib)+1:...
            length(ia)+length(ib)+length(tmp);
    else
        indB{i} = [];
    end
    ib = [ib;tmp];
end
ic = find(C~=0);
if ~isempty(ic)
    indC = length(ia)+length(ib)+1:length(ia)+length(ib)+length(ic);
else
    indC = [];
end
id = [];
indD = cell(length(D),1);
for i=1:length(D)
    tmp = find(D{i}~=0);
    if ~isempty(tmp)
        indD{i} = length(ia)+length(ib)+length(ic)+length(id)+1:...
            length(ia)+length(ib)+length(ic)+length(id)+length(tmp);
    else
        indD{i} = [];
    end
    id = [id;tmp];
end
indself = (length(ia)+length(ib)+length(ic)+length(id))+kernelGroup;

author: MBB-team

repo: VBA-toolbox


text:
function pass = test_printTree(~)
%%TEST_PRINTTREE    Do TREEVAR computations, check that print and printTree works.

%% Basic computation:
u = treeVar();
v = cos(u);
w = sin(u);
t = v + w;

%% Call both static and non-static versions of the method
pass(1) = doesNotCrash(@() treeVar.printTree(t.tree));
pass(2) = doesNotCrash(@() print(t));
%% Introducing differentiation
u = treeVar();
myfun = @(u) 2 + diff(u,2);
t = myfun(u);
pass(3) = doesNotCrash(@() treeVar.printTree(t.tree));
pass(4) = doesNotCrash(@() print(t));
end

function pass = doesNotCrash(fn)
try
    fn();
    pass = true;
catch ME %#ok<NASGU>
    pass = false;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function X2=henon(X1)
%HENON   Henon map (a 2nd-order discrete system)
%
%             x(n+1) = 1 - a*x(n)^2 + y(n)
%             y(n+1) = b*x(n)
%
%        In this demo, a = 1.4, b = 0.3,
%        Initial conditions: x(0) = 0, y(0) = 0
%        Reference values are:  LE1 = 0.418, LE2 = -1.621, LD = 1.26
%
%        The reference values are from the following references:
%
%        [1] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano,
%           "Determining Lyapunov Exponents from a Time Series,"
%            Physica D, Vol. 16, pp. 285-317, 1985.
%
%        [2] Keith Briggs, "An Improved Method for Estimating Liapunov
%            Exponents of Chaotic Time Series," Phys. Lett. A, Vol. 151,
%            pp. 27-32, Nov. 1990.

%        by Steve Wai Kam SIU, Jun. 29, 1998.

%Parameters
a=1.4;
b=0.3;

%Rearrange input data in desired format
%Note: the input data is a column vector
x=X1(1);
y=X1(2);
Q=[X1(3),X1(5);
   X1(4),X1(6)];

%Henon map
X=1-a*(x^2)+y;
Y=b*x;

%Linearized system
J=[-2*a*x, 1;
        b, 0];

%Variational equation
F=J*Q;

%Output data
X2=[X; Y; F(:)];


text:
function gomakegr() % autogenerated function wrapper
    % This M-file displays a message, then calls makegrid to produce a grid and
    % calculate Cumulative Number curves.
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun(mfilename('fullpath'));

    clf
    set(gca,'visible','off');
    gh = gcf;

    txt4  = text(...
        'Position',[0.0 0.65 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'Please wait until cursor changes to a CROSS');

    txt4B = text(...
        'Position',[0.0 0.58 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'and select the edges of a polygon    ');
    txt4 = text(...
        'Position',[0.0 0.51 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'on the map (using the left mouse buton).');
    txt4B  = text(...
        'Position',[0.0 0.36 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'Mac Users: Use the character "p" on your      ');
    txt4  = text(...
        'Position',[0.0 0.29 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'keyboard. Use the right mouse button to select ');
    txt4B  = text(...
        'Position',[0.0 0.22 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String', 'the final point (or charcter "l")           ');

    pause(0.1)
    makegrid

end

author: CelsoReyes