text:
%kvcocASF 'Alternating Sequential Filter of type Closing(Opening(Closing()))'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vcocASF.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% Integer: n 'Size Struct Elem. ', default: 5: 'size of structure element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvcocASF({i, str}, {'i','';'str','';'n',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vcocASF - Alternating Sequential Filter of type Closing(Opening(Closing()))
%
%  DESCRIPTION
% .I vcocASF
% will do N successive close/open/closes by the same structuring element.
%
% Non-sequential filtering can be selected, so that no iteration is done. Otherwise, there will be iteration from 1 to n.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vcocASF -i feath.xv -str b.str -n 5 -o outimage.xv
%
% Will do 5 successive close/open/closes (vcocASF) of the input image "feath.xv" by the structuring element "b.str" and the resulting image will be written to  "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vcocASF
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvcocASF(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvcocASF(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str', '__input';'n', 5;'o', '__output'};
maxval={0,0,2,0};
minval={0,0,2,0};
istoggle=[0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vcocASF"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function S=afni_niml_readsimple(fn,full)
% Reads an AFNI NIML file and returns a 'simple' struct
%
% S=AFNI_NIML_READSIMPLE(FN) reads the AFNI NIML ASCII file FN, and returns
% a 'simple' struct S with the relevant data fields. (This is unlike
% the function AFNI_NIML_READ, which returns all data fields and resembles
% the original file content of FN better).
%
% If the file is present, S will contain the following fields:
%  .data           PxN data for P nodes and N columns (values per node).
%  .node_indices   Px1 indices of P nodes that data refers to (base 0)
%  .history        String with history information
%  .stats          1xN cell with strings of the statistics of the data
%                  (for example, Z or T score).
%  .labels         1xN cell with strings of the labels of the data columns
%  .dset_type      String with the data set type
%
% S=AFNI_NIML_READSIMPLE(FN,true) returns a struct S but with a
% full data matrix (rather than sparse) for all nodes; node indices not
% defined in the file are set to zero. S does not contain a field
% .node_indices.
%
% Please note that this function is *VERY EXPERIMENTAL*, and has only been
% tested with functional data NIML files.
%
% NNO Jan 2010 <n.oosterhof@bangor.ac.uk>

D=afni_niml_read(fn);

if iscell(D)
    if numel(D)==1
        D=D{1};
    else
        error('Cell with multiple elements is not supported');
    end
end

if ~isstruct(D) || ~isfield(D,'dset_type')
    error('Unrecognized input');
end

S=struct();
S.dset_type=D.dset_type;

nodecount=numel(D.nodes);
for k=1:nodecount
    Nk=D.nodes{k};
    if ~isfield(Nk,'name')
        error('Missing name in node element %d',k);
    end

    switch(Nk.name)
        case 'INDEX_LIST'
            S.node_indices=Nk.data;

        case 'SPARSE_DATA'
            S.data=Nk.data;

        case 'AFNI_atr'
            if ~isfield(Nk,'atr_name')
                error('Field %s in node %d has missing atr_name\n', ...
                            Nk.name, D);
            end

            switch Nk.atr_name
                case 'HISTORY_NOTE'
                    S.history=get_string_cell_data(Nk.data);

                case 'COLMS_STATSYM'
                    S.stats=get_string_cell_data(Nk.data);

                case 'COLMS_LABS'
                    S.labels=get_string_cell_data(Nk.data);

            end

        otherwise
            S.misc.(Nk.name)=Nk.data;
    end
end

if nargin>1 && full && isfield(S,'node_indices')
    data=zeros(max(S.node_indices)+1,size(S.data,2));
    data(S.node_indices+1,:)=S.data;
    S.data=data;
    S=rmfield(S,'node_indices');
end

function labels=get_string_cell_data(data)
    while iscell(data) && numel(data)==1 && ~ischar(data)
        data=data{1};
    end

    labels=split_string(data,';');
    if numel(labels{end})==0        % as above
        labels=labels(1:(end-1));
    end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function sparse_grid_total_poly_test ( )

%*****************************************************************************80
%
%% SPARSE_GRID_TOTAL_POLY_TEST tests the SPARSE_GRID_TOTAL_POLY library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_TOTAL_POLY_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the SPARSE_GRID_TOTAL_POLY library.\n' );

  rules_test ( );
  point_quality_test ( );
  weight_quality_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_TOTAL_POLY_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function bf_save_path(BF,path, overwrite)
% Saves BF data in a mat file
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_save_path.m 7703 2019-11-22 12:06:29Z guillaume $

if nargin == 2 && exist(path, 'file')
    save(path, '-struct', 'BF', '-append');
else
    save(path, '-struct', 'BF', '-v7.3');
end

author: spm


text:
% ========================================================================
% Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, TIP,2020
% algorithm Version 1.0
% Copyright(c) 2020, Hui Li, Kede Ma, Yongwei Yong and Lei Zhang
% All Rights Reserved.
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
% Please refer to the following paper:
% H. Li et al., "Fast Multi-Scale Structural Patch Decomposition for Multi-Exposure Image Fusion, 2020" In press
% IEEE Transactions on Image Processing
% Please kindly report any suggestions or corrections to xiaohui102788@126.com
%----------------------------------------------------------------------

function [fI3, i_meant,aa,N1] = scale_interm(imgSeqColor,r)

[h,w,n]=size(imgSeqColor);
N = boxfilter(ones(h, w), r);

tem=ones(h, w);
tem(:,2:2:w)=0;
tem(2:2:h,:)=0;
N1= boxfilter(tem, r);

p=4;

[WD, Cmax,i_mean2]= weight_cal_detail(imgSeqColor,p,r);
WD=WD.*repmat(Cmax,[1 1 n]);

% WD= weight_cal_gray(imgSeqColor,[0 0 0 1 1 0 0]);

% WB= weight_cal_gray(imgSeqColor,[0 0 0 0 0 1 0]);
% WD=WB;

F_temp2_detail=zeros(h,w,n);
% F_temp2_base=zeros(h,w,n);

%% approximate aggregation through averaging(mean filter) the weight map

i_meant=zeros(ceil(h/2),ceil(w/2),2);
% tic
for i = 1:n

    aa=i_mean2(:,:,i).*tem;
    i_meant(:,:,i)=aa(1:2:h,1:2:w);
    %      i_meant2(:,:,i)=i_meant(:,:,i).*(1-edg(:,:,i));

    %     W_D1=boxfilter(aa.*WD(:,:,i), r)./ N;
    %     W_D2=boxfilter(tem.*WD(:,:,i), r)./ N;

    W_D1=boxfilter(i_mean2(:,:,i).*WD(:,:,i), r)./ N;
    W_D2=boxfilter(WD(:,:,i), r)./ N;

    F_temp2_detail(:,:,i)=W_D2.*(imgSeqColor(:,:,i))-W_D1;

end

fI3=sum(F_temp2_detail,3);
% toc

%figure,imshow(mat2gray(fI3))

end

author: thfylsty


text:
function test_failed=test_zak

disp(' ===============  TEST_ZAK ================');

Lr=[1,12,12,12,12,12,15,15,15];
Kr=[1, 1, 2, 3, 4, 6, 3, 5,15];

ref_failed=0;
inv_failed=0;

W=1;

for ii=1:length(Lr)

  L=Lr(ii);
  K=Kr(ii);

  f=tester_rand(L,W)+i*rand(L,W)-.5-i*.5;

  ccref=ref_zak(f,K);
  cc=zak(f,K);
  r=izak(cc);

  res=ccref-cc;

  nres=norm(res(:));
  ninv=norm(f-r);

  s=sprintf('REF L:%3i K:%3i %0.5g',L,K,nres);
  disp(s)

  [ref_failed,fail]=ltfatdiditfail(nres,ref_failed);
  % if nres>10e-10
    % disp('FAILED');
    % ref_failed=ref_failed+1;
  % end;

  s=sprintf('INV L:%3i K:%3i %0.5g',L,K,ninv);
  disp(s)

  [inv_failed,fail]=ltfatdiditfail(ninv,inv_failed);
  % if nres>10e-10
    % disp('FAILED');
    % inv_failed=inv_failed+1;
  % end;

end;

test_failed=ref_failed+inv_failed;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [r,normr,nre,s] = reorth(Q,r,normr,index,alpha,method)

%REORTH   Reorthogonalize a vector using iterated Gram-Schmidt
%
%   [R_NEW,NORMR_NEW,NRE] = reorth(Q,R,NORMR,INDEX,ALPHA,METHOD)
%   reorthogonalizes R against the subset of columns of Q given by INDEX.
%   If INDEX==[] then R is reorthogonalized all columns of Q.
%   If the result R_NEW has a small norm, i.e. if norm(R_NEW) < ALPHA*NORMR,
%   then a second reorthogonalization is performed. If the norm of R_NEW
%   is once more decreased by  more than a factor of ALPHA then R is
%   numerically in span(Q(:,INDEX)) and a zero-vector is returned for R_NEW.
%
%   If method==0 then iterated modified Gram-Schmidt is used.
%   If method==1 then iterated classical Gram-Schmidt is used.
%
%   The default value for ALPHA is 0.5.
%   NRE is the number of reorthogonalizations performed (1 or 2).

% References:
%  Aake Bjorck, "Numerical Methods for Least Squares Problems",
%  SIAM, Philadelphia, 1996, pp. 68-69.
%
%  J.~W. Daniel, W.~B. Gragg, L. Kaufman and G.~W. Stewart,
%  ``Reorthogonalization and Stable Algorithms Updating the
%  Gram-Schmidt QR Factorization'', Math. Comp.,  30 (1976), no.
%  136, pp. 772-795.
%
%  B. N. Parlett, ``The Symmetric Eigenvalue Problem'',
%  Prentice-Hall, Englewood Cliffs, NJ, 1980. pp. 105-109

%  Rasmus Munk Larsen, DAIMI, 1998.

% Check input arguments.
% warning('PROPACK:NotUsingMex','Using slow matlab code for reorth.')
if nargin<2
  error('Not enough input arguments.')
end
[n k1] = size(Q);
if nargin<3 | isempty(normr)
%  normr = norm(r);
  normr = sqrt(r'*r);
end
if nargin<4 | isempty(index)
  k=k1;
  index = [1:k]';
  simple = 1;
else
  k = length(index);
  if k==k1 & index(:)==[1:k]'
    simple = 1;
  else
    simple = 0;
  end
end
if nargin<5 | isempty(alpha)
  alpha=0.5; % This choice garanties that
             % || Q^T*r_new - e_{k+1} ||_2 <= 2*eps*||r_new||_2,
             % cf. Kahans ``twice is enough'' statement proved in
             % Parletts book.
end
if nargin<6 | isempty(method)
   method = 0;
end
if k==0 | n==0
  return
end
if nargout>3
  s = zeros(k,1);
end

normr_old = 0;
nre = 0;
while normr < alpha*normr_old | nre==0
  if method==1
    if simple
      t = Q'*r;
      r = r - Q*t;
    else
      t = Q(:,index)'*r;
      r = r - Q(:,index)*t;
    end
  else
    for i=index,
      t = Q(:,i)'*r;
      r = r - Q(:,i)*t;
    end
  end
  if nargout>3
    s = s + t;
  end
  normr_old = normr;
%  normr = norm(r);
  normr = sqrt(r'*r);
  nre = nre + 1;
  if nre > 4
    % r is in span(Q) to full accuracy => accept r = 0 as the new vector.
    r = zeros(n,1);
    normr = 0;
    return
  end
end

author: andrewssobral

repo: lrslibrary


text:
function dists = eucl(crds1,crds2)
% EUCL: Calculates the euclidean distances among a set of points, or between a
%       reference point and a set of points, or among all possible pairs of two
%       sets of points, in P dimensions.  Returns a single distance for two points.
%
%     Syntax: dists = eucl(crds1,crds2)
%
%        crds1 = [N1 x P] matrix of point coordinates.  If N=1, it is taken to
%                   be the reference point.
%        crds2 = [N2 x P] matrix of point coordinates.  If N=1, it is taken to
%                   be the reference point.
%        -----------------------------------------------------------------------
%        dists = [N1 x N1] symmetric matrix of pairwise distances (if only crds1
%                   is specified);
%                [N1 x 1]  col vector of euclidean distances (if crds1 & ref
%                   are specified);
%                [1 x N2]  row vector of euclidean distances (if ref & crds2
%                   are specified);
%                [N1 x N2] rectangular matrix of pairwise distances (if crds1
%                   & crds2 are specified);
%                [1 x 1]   scalar (if crds1 is a [2 x P] matrix or ref1 & ref2
%                   are specified);
%

% RE Strauss, 5/4/94
%   10/28/95 - output row (rather than column) vector for the (reference
%               point)-(set of points) case; still outputs column vector for the
%               (set of points)-(reference point) case.
%   10/30/95 - for double for-loops, put one matrix-col access in outer loop
%               to increase speed.
%   10/12/96 - vectorize inner loop to increase speed.
%    6/12/98 - allow for P=1.
%   11/11/03 - initialize dists to NaN for error return.

  if nargin == 0, help eucl; return; end
  dists = NaN;

  if (nargin < 2)                     % If only crds1 provided,
    [N,P] = size(crds1);
    if (N<2)
      error('  EUCL: need at least two points');
    end

    crds1 = crds1';                     % Transpose crds
    dists = zeros(N,N);                 % Calculate pairwise distances

    for i = 1:N-1
      c1 = crds1(:,i) * ones(1,N-i);
      if (P>1)
        d = sqrt(sum((c1-crds1(:,(i+1:N))).^2));
      else
        d = abs(c1-crds1(:,(i+1:N)));
      end
      dists(i,(i+1:N)) = d;
      dists((i+1:N),i) = d';
    end
    if (N==2)                            % Single distance for two points
      dists = dists(1,2);
    end

  else                                % If crds1 & crds2 provided,
    [N1,P1] = size(crds1);
    [N2,P2] = size(crds2);
    if (P1~=P2)
      error('  EUCL: sets of coordinates must be of same dimension');
    else
      P = P1;
    end

    crds1 = crds1';                     % Transpose crds
    crds2 = crds2';

    if (N1>1 && N2>1)                    % If two matrices provided,
      dists = zeros(N1,N2);               % Calc all pairwise distances between them
      for i = 1:N1
        c1 = crds1(:,i) * ones(1,N2);
        if (P>1)
          d = sqrt(sum((c1-crds2).^2));
        else
          d = abs(c1-crds2);
        end
        dists(i,:) = d;
      end
    end

    if (N1==1 && N2==1)                  % If two vectors provided,
      dists = sqrt(sum((crds1-crds2).^2));  % Calc scalar
    end

    if (N1>1 && N2==1)                   % If matrix && reference point provided,
      crds1 = crds1 - (ones(N1,1)*crds2')'; % Center points on reference point
      if (P>1)                              % Calc euclidean distances in P-space
         dists = sqrt(sum(crds1.^2))';
      else
         dists = abs(crds1)';
      end
    end;                                    % Return column vector

    if (N1==1 && N2>1)                   % If reference point && matrix provided,
      crds2 = crds2 - (ones(N2,1)*crds1')'; % Center points on reference point
      if (P>1)                              % Calc euclidean distances in P-space
        dists = sqrt(sum(crds2.^2));
      else
        dists = abs(crds2);
      end
    end;                                    % Return row vector
  end

  return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function X = parseargs(X,varargin)
%PARSEARGS - Parses name-value pairs
%
% Behaves like setfield, but accepts multiple name-value pairs and provides
% some additional features:
% 1) If any field of X is an cell-array of strings, it can only be set to
%    one of those strings.  If no value is specified for that field, the
%    first string is selected.
% 2) Where the field is not empty, its data type cannot be changed
% 3) Where the field contains a scalar, its size cannot be changed.
%
% X = parseargs(X,name1,value1,name2,value2,...)
%
% Intended for use as an argument parser for functions which multiple options.
% Example usage:
%
% function my_function(varargin)
%   X.StartValue = 0;
%   X.StopOnError = false;
%   X.SolverType = {'fixedstep','variablestep'};
%   X.OutputFile = 'out.txt';
%   X = parseargs(X,varargin{:});
%
% Then call (e.g.):
%
% my_function('OutputFile','out2.txt','SolverType','variablestep');

% The various #ok comments below are to stop MLint complaining about
% inefficient usage.  In all cases, the inefficient usage (of error, getfield,
% setfield and find) is used to ensure compatibility with earlier versions
% of MATLAB.

% Copyright 2006-2010 The MathWorks, Inc.

remaining = nargin-1; % number of arguments other than X
count = 1;
fields = fieldnames(X);
modified = zeros(size(fields));
% Take input arguments two at a time until we run out.
while remaining>=2
    fieldname = varargin{count};
    fieldind = find(strcmp(fieldname,fields));
    if ~isempty(fieldind)
        oldvalue = getfield(X,fieldname); %#ok
        newvalue = varargin{count+1};
        if iscell(oldvalue)
            % Cell arrays must contain strings, and the new value must be
            % a string which appears in the list.
            if ~iscellstr(oldvalue)
                error(sprintf('All allowed values for "%s" must be strings',fieldname));  %#ok
            end
            if ~ischar(newvalue)
                error(sprintf('New value for "%s" must be a string',fieldname));  %#ok
            end
            if isempty(find(strcmp(oldvalue,newvalue))) %#ok
                error(sprintf('"%s" is not allowed for field "%s"',newvalue,fieldname));  %#ok
            end
        elseif ~isempty(oldvalue)
            % The caller isn't allowed to change the data type of a non-empty property,
            % and scalars must remain as scalars.
            if ~strcmp(class(oldvalue),class(newvalue))
                error(sprintf('Cannot change class of field "%s" from "%s" to "%s"',...
                    fieldname,class(oldvalue),class(newvalue))); %#ok
            elseif numel(oldvalue)==1 & numel(newvalue)~=1 %#ok
                error(sprintf('New value for "%s" must be a scalar',fieldname));  %#ok
            end
        end
        X = setfield(X,fieldname,newvalue); %#ok
        modified(fieldind) = 1;
    else
        error(['Not a valid field name: ' fieldname]);
    end
    remaining = remaining - 2;
    count = count + 2;
end
% Check that we had a value for every name.
if remaining~=0
    error('Odd number of arguments supplied.  Name-value pairs required');
end

% Now find cell arrays which were not modified by the above process, and select
% the first string.
notmodified = find(~modified);
for i=1:length(notmodified)
    fieldname = fields{notmodified(i)};
    oldvalue = getfield(X,fieldname); %#ok
    if iscell(oldvalue)
        if ~iscellstr(oldvalue)
            error(sprintf('All allowed values for "%s" must be strings',fieldname)); %#ok
        elseif isempty(oldvalue)
            error(sprintf('Empty cell array not allowed for field "%s"',fieldname)); %#ok
        end
        X = setfield(X,fieldname,oldvalue{1}); %#ok
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% StackExchange Signal Processing Q63549
% https://dsp.stackexchange.com/questions/63549
% Computationally Efficient Ways of Determining If Pixels Are Clumped
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     06/01/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

imageFileName001 = 'Image001.png';
imageFileName002 = 'Image002.png';

%% Generate Data

mA = imread(imageFileName001);
mB = imread(imageFileName002);

mA = ~logical(mA(:, :, 1));
mB = ~logical(mB(:, :, 1));

% Removing outer frame
mA(1, :) = false;
mA(end, :) = false;
mA(:, 1) = false;
mA(:, end) = false;

mB(1, :) = false;
mB(end, :) = false;
mB(:, 1) = false;
mB(:, end) = false;

% Extracting the properties of the Connected Components
% Can use half of 'MajorAxisLength' for the Radius.
% Then from 'PixelIdxList' do the computation
%{
objValA = CalcImageObjValImgBinaryImageProps(mA);
objValB = CalcImageObjValImgBinaryImageProps(mB);
%}

% Doing the Same, just with a nice trick to get the circle parameters by
% Binray Image Distnace Transform of the inverted image.
objValA = CalcImageObjValImgBinaryImageDistanceTransform(mA);
objValB = CalcImageObjValImgBinaryImageDistanceTransform(mB);

%% Display Results

hFigure     = figure('Position', [100, 100, 280, 320]); %<! [x, y, width, height]
hAxes       = axes('Units', 'pixels', 'Position', [17, 10, 256, 256]); %<! [x, y, width, height]
hImageObj   = imagesc(mA);
set(get(hAxes, 'Title'), 'String', {['Input Image A'], ['Objective Value - ', num2str(objValA)]}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'DataAspectRatio', [1, 1, 1]);
set(hAxes, 'XTick', [], 'YTick', [], 'XTickLabel', [], 'YTickLabel', []);
% set(hAxes, 'LooseInset', get(hAxes, 'TightInset'));
set(hAxes, 'LooseInset', [0.05, 0.05, 0.05, 0.05]);

hFigure     = figure('Position', [100, 100, 280, 320]); %<! [x, y, width, height]
hAxes       = axes('Units', 'pixels', 'Position', [17, 10, 256, 256]); %<! [x, y, width, height]
hImageObj   = imagesc(mB);
set(get(hAxes, 'Title'), 'String', {['Input Image B'], ['Objective Value - ', num2str(objValB)]}, ...
    'FontSize', fontSizeTitle);
set(hAxes, 'DataAspectRatio', [1, 1, 1]);
set(hAxes, 'XTick', [], 'YTick', [], 'XTickLabel', [], 'YTickLabel', []);
% set(hAxes, 'LooseInset', get(hAxes, 'TightInset'));
set(hAxes, 'LooseInset', [0.05, 0.05, 0.05, 0.05]);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [x, y, z] = pol2cart(th, r, z)
%POL2CART   Transform polar to Cartesian coordinates for CHEBFUN2 objects.
%
%   [X,Y] = POL2CART(TH,R) transforms corresponding elements of data stored in
%   polar coordinates (angle TH, radius R) to Cartesian coordinates X,Y.  The
%   arrays TH and R must the same size (or either can be scalar).  TH must be in
%   radians.
%
%   [X,Y,Z] = POL2CART(TH,R,Z) transforms corresponding elements of data stored
%   in cylindrical coordinates (angle TH, radius R, height Z) to Cartesian
%   coordinates X,Y,Z. The arrays TH, R, and Z must be the same size (or any of
%   them can be scalar).  TH must be in radians.
%
% See also SPH2CART.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Polar coordinates:
x = r .* cos( th );
y = r .* sin( th );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function som_stats_plot(csS,plottype,varargin)

%SOM_STATS_PLOT Plots of data set statistics.
%
% som_stats_plot(csS, plottype, [argID, value, ...])
%
%  som_stats_plot(csS,'stats')
%  som_stats_plot(csS,'stats','p','vert','color','r')
%
%  Input and output arguments ([]'s are optional):
%   csS         (cell array) of statistics structs
%               (struct) a statistics struct
%   plottype    (string) some of the following
%                        'hist'   histogram
%                        'box'    min, max, mean, and std shown as a boxplot
%                        'stats'  both histogram (with black) and the boxplot
%   [argID, (string) See below. The values which are unambiguous can
%    value] (varies) be given without the preceeding argID.
%
% Here are the valid argument IDs and corresponding values. The values which
% are unambiguous (marked with '*') can be given without the preceeding argID.
%   'counts'      *(string) 'c' (for counts, the default) or 'p' (for percentages)
%   'color'        (vector) size 1 x 3, color to be used
%                  (string) a color string
%   'title'        (string) 'on' (default) or 'off'
%   'orientation' *(string) 'horiz' or 'vert' (default): orientation for the
%                           bin values (horizontally or vertically)
%
% See also  SOM_STATS, SOM_STATS_TABLE, SOM_TABLE_PRINT, SOM_STATS_REPORT.

% Contributed to SOM Toolbox 2.0, December 31st, 2001 by Juha Vesanto
% Copyright (c) by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 311201

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% arguments

% statistics
if isstruct(csS), csS = {csS}; end

% default values
useprob   = 0;
color     = [0 0 1];
showtitle = 1;
horiz     = 0;

% varargin
i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
     % argument IDs
     case 'counts',      i=i+1; useprob = strcmp(varargin{i}(1),'p');
     case 'color',       i=i+1; color = varargin{i};
     case 'title',       i=i+1; showtitle = strcmp(varargin{i},'on');
     case 'orientation', i=i+1; horiz = strcmp(varargin{i},'horiz');
     % unambiguous values
     case {'horiz','vert'}, horiz = strcmp(varargin{i},'horiz');
     case {'c','p'}, useprob = strcmp(varargin{i}(1),'p');
     otherwise argok=0;
    end
  elseif isstruct(varargin{i}) && isfield(varargin{i},'type'),
    argok = 0;
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_stats_plot) Ignoring invalid argument #' num2str(i+2)]);
  end
  i = i+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%% action

ss = ceil(sqrt(length(csS))); ss = [ss, ceil(length(csS)/ss)];

for j = 1:length(csS),
    sS = csS{j};
    subplot(ss(1),ss(2),j);
    switch plottype,
    case 'stats',
        cla, hold on
        Counts = sS.hist.counts;
        if useprob, for i=1:size(Counts,2), Counts(:,i) = Counts(:,i)/sum(Counts(:,i)); end, end
        hist_plot(sS.hist.bins,sS.hist.binlabels,Counts,color);
        box_plot(sS.min,sS.max,sS.mean,sS.std,[0 0 0]);
    case 'hist',
        cla, hold on
        Counts = sS.hist.counts;
        if useprob, for i=1:size(Counts,2), Counts(:,i) = Counts(:,i)/sum(Counts(:,i)); end, end
        hist_plot(sS.hist.bins,sS.hist.binlabels,Counts,color);
    case 'box',
        cla
	box_plot(sS.min,sS.max,sS.mean,sS.std,color);
    end
    if showtitle, title(sprintf('%s (valid: %d/%d)',sS.name,sS.nvalid,sS.ntotal)); end
    if ~horiz, view(90,-90); end
    a = axis; a(1) = sS.min; a(2) = sS.max; axis(a);
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%% subfunctions

function hist_plot(bins,binlabels,Counts,color)

    if nargin<4, color = jet(size(Counts,2)); end
    h = bar(bins,Counts);
    for j=1:length(h), set(h(j),'facecolor',color(j,:),'edgecolor','none'); end
    a = axis; a(3:4) = [0 max(Counts(:))]; axis(a);
    set(gca,'XTick',bins,'XTickLabel',binlabels);
    return;

function vstr = numtostring(v,d)

    nearzero = (abs(v)/(max(v)-min(v)) < 10.^-d);
    i1 = find(v > 0 & nearzero);
    i2 = find(v < 0 & nearzero);
    vstr = strrep(cellstr(num2str(v,d)),' ','');
    vstr(i1) = {'0.0'};
    vstr(i2) = {'-0.0'};
    return;

function box_plot(mi,ma,me,st,Color)

    if nargin < 5, Color = jet(length(mi)); end
    a = axis;
    y = linspace(a(3),a(4),length(mi)+2); y = y(2:end);
    d = (y(2)-y(1))/20;
    for i=1:length(mi),
        h1 = line([mi(i) ma(i)],[y(i) y(i)]);
        h2 = line([mi(i) mi(i) NaN ma(i) ma(i)],[y(i)-d y(i)+d NaN y(i)-d y(i)+d]);
        h3 = line([me(i)-st(i) me(i)+st(i)],[y(i) y(i)]);
        h4 = line([me(i) me(i)],[y(i)-2*d y(i)+2*d]);
        set([h1 h2 h3 h4],'color',Color(i,:));
        set([h1 h2],'linewidth',1);
        set([h3 h4],'linewidth',3);
    end
    return;

author: ilarinieminen

repo: SOM-Toolbox


text:
function [voxel] = spm_vb_get_Ab(Y,slice)
% Get A and b quantities - average prediction errors from AR model
% FORMAT [voxel] = spm_vb_get_Ab(Y,slice)
%
% Y      - [T x N] time series
% slice  - data structure (see spm_vb_glmar)
%
% voxel(n).A
% voxel(n).b
%
% The above quantities are estimated using pre-computed
% cross-covariance matrices
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny and Nelson Trujillo-Barreto
% $Id: spm_vb_get_Ab.m 6079 2014-06-30 18:25:37Z spm $

k = slice.k;
N = slice.N;

for n=1:N
    % Equation 63 of paper VB1 but implemented
    % efficiently using cross-covariance method described in paper VB3
    if isfield(slice.I,'A2_tilde')
        A2_tilde  = slice.I.A2_tilde(:,:,n);
    else
        A2_tilde  = reshape(slice.I.S*slice.a2{n}(:),k,k);
    end
    if isfield(slice.I,'A3a_tilde')
        A3a_tilde = slice.I.A3a_tilde(:,:,n);
    else
        A3a_tilde = -reshape(slice.I.R1*slice.ap_mean(:,n),k,k);
    end
    voxel(n).A    = slice.I.xtx+A2_tilde+A3a_tilde+A3a_tilde';

    % Equation 64 of paper VB1 but implemented
    % efficiently using cross-covariance method described in paper VB3
    b2_tilde   = -slice.I.rxy(:,:,n)'*slice.ap_mean(:,n);
    b3_tilde   = -slice.I.Gxy(:,:,n)'*slice.ap_mean(:,n);
    b4_tilde   = slice.I.D(:,:,n)*slice.a2{n}(:);
    voxel(n).b = slice.I.gxy(:,n)+b2_tilde+b3_tilde+b4_tilde;
end


text:
function colors = getColors(n)
hf = figure('visible','off');
ha = axes('parent',hf);
colors = get(ha,'colororder');
colors = repmat(colors,ceil(n/size(colors,1)),1);
colors = colors(1:n,:);
delete(hf)

author: MBB-team

repo: VBA-toolbox


text:
function varargout = source2sparse(varargin)

% This function is a backward compatibility wrapper. It allows existing
% MATLAB scripts that do not use the new FieldTrip ft_xxx function naming
% scheme to work with recent versions of the FieldTrip toolbox.
%
% Please look in ft_xxx for the help of the function that you are looking
% for, where xxx is the name of the function that you were looking for.

% Copyright (C) 2009, Robert Oostenveld

prefix    = 'ft_';
funname   = mfilename;
warning([upper(mfilename) ' is only a compatibility wrapper, which will soon be removed. Please instead call ' upper(prefix) upper(funname) '.']);
funhandle = str2func([prefix funname]);
[varargout{1:nargout}] = funhandle(varargin{:});


text:
function [L,E,EMAP] = crouzeix_raviart_cotmatrix(V,F)
  % CROUZEIX_RAVIART_COTMATRIX Compute the Crouzeix-Raviart cotangent
  % Laplacian matrix where we use test functions define at edge midpoints.
  %
  % See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
  % Harmon, Zorin, Grinspun 2007]
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by element-size list of triangle indices
  % Outputs:
  %   L  #E by #E edge-based sparse cotangent matrix
  %   E  #E by 2 list of edges
  %
  % Examples:
  %   % mesh in (V,F)
  %   [Lcr,E,EMAP] = crouzeix_raviart_cotmatrix(V,F);
  %   [Mcr,E,EMAP] = crouzeix_raviart_massmatrix(V,F);
  %   [Ucr,~] = eigs(Lcr,Mcr,5,'sm');
  %   % Convert between edge values and vertex values
  %   E2V = sparse(E(:),repmat(1:size(E,1),1,2)',1,size(V,1),size(E,1));
  %   E2V = bsxfun(@rdivide,E2V,sum(E2V,2));
  %   V2E = sparse(E(:),repmat(1:size(E,1),1,2)',1,size(V,1),size(E,1))';
  %   V2E = bsxfun(@rdivide,V2E,sum(V2E,2));
  %   tsurf(F,[V(:,1:2) E2V*Ucr(:,end-1)])
  %
  %   % Display discontinous solution
  %   FF = reshape(1:numel(F),size(F));
  %   VV = V(F,:);
  %   EMAP = reshape(EMAP,size(F));
  %   A = sparse( ...
  %     [FF FF FF], ...
  %     EMAP(:,[1 2 3 2 3 1 3 1 2]), ...
  %     repmat(-[1 1 1 -1 -1 -1 -1 -1 -1],size(F,1),1), ...
  %     size(VV,1), ...
  %     size(E,1));
  %   tsurf(FF,[VV(:,1:2) A*Ucr(:,end-1)])
  %
  %
  %
  % See also: edge_laplacian, is_boundary_edge, crouzeix_raviart_massmatrix,
  %   cotmatrix
  %

  %% check for non-manifold edges
  %S = statistics(V,F,'Fast',true);
  %if S.num_nonmanifold_edges > 0
  %  error(sprintf('There are %d non-manifold edges',num_nonmanifold_edges));
  %end

  % number of vertices
  n = size(V,1);
  % number of elements
  m = size(F,1);
  % simplex size
  ss = size(F,2);
  switch ss
  case 3
    % Compute cotangents: seems to be 0.5*C
    C = 2*cotangent(V,F);

    allE = [F(:,[2 3]);F(:,[3 1]);F(:,[1 2])];
    % Map each face-edge to a unique edge
    F2E = reshape(1:3*m,m,3);
    % Lij = -2 cot aij
    %
    % o
    % |\
    % | \
    % |  \
    % |   \
    % i    \
    % |     \
    % |      \
    % |αij    \
    % o----j---o
    %
    %
    LI =  F2E(:, [1 2 3 2 3 1        1 2 3 2 3 1]);
    LJ =  F2E(:, [2 3 1 1 2 3        1 2 3 2 3 1]);
    LV = 2*[-C(:,[3 1 2 3 1 2]) C(:,[3 1 2 3 1 2])];

    % Map duplicate edges to first instance
    [E,~,EMAP] = unique(sort(allE,2),'rows');

    assert(all(size(LI)==size(LJ)));
    assert(all(size(LI)==size(LV)));
    L = sparse(EMAP(LI),EMAP(LJ),LV,size(E,1),size(E,1));
  case 4
    % tets
    T = F;
    C = -2*cotangent(V,T);
    allF = [ ...
      T(:,2) T(:,4) T(:,3); ...
      T(:,1) T(:,3) T(:,4); ...
      T(:,1) T(:,4) T(:,2); ...
      T(:,1) T(:,2) T(:,3); ...
      ];
    % Map each element-face to a unique face
    T2F = reshape(1:4*m,m,4);
    % Lij = -2 lij cot αij
    LI = T2F(:,[1 4 4 4 2 3 2 1 2 3 3 1        1 4 4 4 2 3 2 1 2 3 3 1]);
    LJ = T2F(:,[2 1 2 3 3 1 1 4 4 4 2 3        1 4 4 4 2 3 2 1 2 3 3 1]);
    LV = [-C(:,[3 4 5 6 1 2 3 4 5 6 1 2]) C(:,[3 4 5 6 1 2 3 4 5 6 1 2])];
    % Map duplicate facets to first instance
    [F,~,FMAP] = unique(sort(allF,2),'rows');

    assert(all(size(LI)==size(LJ)));
    assert(all(size(LI)==size(LV)));
    L = sparse(FMAP(LI),FMAP(LJ),LV,size(F,1),size(F,1));
    E = F;
    EMAP = FMAP;

  otherwise
    error(['Simplex size ' num2str(ss) ' unsupported']);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function option = mfemoption(option)

if ~isfield(option,'elemType')
    option.elemType = 'RT0';
end

if ~isfield(option,'maxIt')
    option.maxIt = 4;
end

if ~isfield(option,'maxN')
    option.maxN = 1e5;
end

if ~isfield(option,'L0')
    option.L0 = 0;
end

if ~isfield(option,'refType')
    option.refType = 'red';
end

if ~isfield(option,'printlevel')
    option.printlevel = 1;
end

if ~isfield(option,'plotflag')
    option.plotflag = 1;
end

if ~isfield(option,'rateflag')
    option.rateflag = 1;
end

if ~isfield(option,'dispflag')
    option.dispflag = 1;
end

if ~isfield(option,'tol')
    option.tol = 1e-8;
end

repo: ifem


text:
function a = fourier_sine_inverse ( n )

%*****************************************************************************80
%
%% FOURIER_SINE_INVERSE returns the inverse of the FOURIER_SINE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = fourier_sine ( n );

  return
end


text:
function [pts] = annotation_load( fname, data_name )

if strcmp(data_name,'helen') == 1
    pts = helen_annotation_load( fname );
elseif strcmp(data_name,'lfpw') == 1
    pts = lfpw_annotation_load( fname );
elseif strcmp(data_name,'w300') == 1
    pts = w300_annotation_load( fname );
else
    disp(['dataset name: ' data_name ' is invalid!']);
end

end

author: tntrung


text:
function [params] = setAttnParams(params)
  if params.attnFunc>0
    params.numSrcHidVecs = params.srcMaxLen-1;

    if params.attnGlobal % global
      params.numAttnPositions = params.numSrcHidVecs;
    else % local
      params.numAttnPositions = 2*params.posWin + 1;
    end
  else
    params.numSrcHidVecs = 0;
  end
end

author: lmthang


text:
function [mfcMVN, meanV, varV] = SideMVNcore(mfcRaw)

meanAcc = 0;    varAcc = 0;
for j=1:length(mfcRaw)
    nFr(j) = size(mfcRaw{j},1);
    meanAcc = meanAcc + sum(mfcRaw{j});
    varAcc = varAcc + sum(mfcRaw{j}.^2);
end
meanV = meanAcc / sum(nFr);
varV = sqrt( varAcc/sum(nFr) - meanV.^2 );
precision = 1./varV;

% Normalize the mean and variance
for j=1:length(mfcRaw)
    mfcMVN{j} = bsxfun(@minus, mfcRaw{j}, meanV);
    mfcMVN{j} = bsxfun(@times, mfcMVN{j}, precision);
end

author: singaxiong

repo: SignalGraph


text:
%MAIN
%
% This script is used for testing the continuous, finite horizon, LQR
% solver that I'm working on.
%
%
clear; clc;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                       User Set Parameters                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

nSim = 25;
testPerturbation = 0.1;  %Initial position error amplitude

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                       Draw the dynamical system                         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

n = 13;
xLim = [-4,1.5];
yLim = [-4,1];
x = linspace(xLim(1),xLim(2),n);
y = linspace(yLim(1),yLim(2),n);

[xx,yy] = ndgrid(x,y);

[dxx, dyy] = dynamics(xx,yy,zeros(n,n));

figure(1); clf; hold on;
title('Dynamical System - Simulations');
quiver(xx,yy,dxx,dyy);
xlabel('x');
ylabel('y');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Compute a reference trajectory                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

options = odeset();
options.Events = @(t,z)eventFunc(t,z,xLim,yLim);
tSpan = [0,-10];

control = @(t) -2*sin(0.4*t-0.1)+0.6*cos(0.9*t+0.1);

userFunc = @(t,z)rhs(t,z,control(t));
z0 = [0;0];

sol = ode45(userFunc,tSpan,z0,options);

nTime = 70;
tSol = linspace(sol.x(end),sol.x(1),nTime);
zSol = deval(sol,tSol);
xSol = zSol(1,:);
ySol = zSol(2,:);
uSol = control(tSol);
tSpan = [tSol(1),tSol(end)];

plot(xSol,ySol,'r-','LineWidth',3);
plot(0,0,'ro','MarkerSize',10,'LineWidth',2)
axis(1.2*[xLim,yLim]); axis equal;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                      Store the reference trajectory                     %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%We will store each trajectory as a polynomial. For low-order fits, say
%less than 10th order, Matlab's polyval and polyfit are good. For
%higher-order fit's, it is best to use different method, based on
%barycentric interpolation. Google chebyfun.

nFit = 5;  %Order of polynomial fitting
xFit = polyfit(tSol,xSol,nFit);
yFit = polyfit(tSol,ySol,nFit);
uFit = polyfit(tSol,uSol,nFit);

figure(2); clf;
subplot(3,1,1); hold on;
plot(tSol,xSol,'k.','MarkerSize',15)
plot(tSol,polyval(xFit,tSol),'r-')
xlim(tSpan);
ylabel('x')
title('Polynomial approximation of trajectory')
subplot(3,1,2); hold on;
plot(tSol,ySol,'k.','MarkerSize',15)
plot(tSol,polyval(yFit,tSol),'r-')
xlim(tSpan);
ylabel('y')
subplot(3,1,3); hold on;
plot(tSol,uSol,'k.','MarkerSize',15)
plot(tSol,polyval(uFit,tSol),'r-')
xlim(tSpan);
ylabel('u')
xlabel('t')

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Compute LQR Gains along trajectory                   %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%Function handle for getting linearized system dynamics
linSys = @(t)getLinTraj(t,xFit,yFit,uFit);

Q = eye(2);  % Running cost on state
R = 1;       % Running cost on input
F = eye(2);  % Terminal cost on state
tol = 1e-6;  % Accuracy of ricatti propagation

Soln = trajectoryLqr(tSol,linSys,Q,R,F,tol);

K = reshape([Soln.K],2,nTime);
kxFit = polyfit(tSol,K(1,:),nFit);
kyFit = polyfit(tSol,K(2,:),nFit);

figure(3); clf;
subplot(2,1,1); hold on;
plot(tSol,K(1,:),'k.','MarkerSize',15)
plot(tSol,polyval(kxFit,tSol),'r-')
xlim(tSpan);
ylabel('Kx')
title('Polynomial approximation of lqr gains')
subplot(2,1,2); hold on;
plot(tSol,K(2,:),'k.','MarkerSize',15)
plot(tSol,polyval(kyFit,tSol),'r-')
xlim(tSpan);
ylabel('Ky')
xlabel('t')

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Compare to infinite-horizon LQR gains                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

ssSoln(nTime).t = 0;
ssSoln(nTime).K = zeros(1,2);
ssSoln(nTime).S = zeros(2);
ssSoln(nTime).E = zeros(2,1);
for i=1:nTime
    tNow = tSol(i);
    [A,B] = linSys(tNow);
    [KK,SS,EE] = lqr(A,B,Q,R);
    ssSoln(i).t = 0;
    ssSoln(i).K = KK;
    ssSoln(i).S = SS;
    ssSoln(i).E = EE;
end

ssK = reshape([ssSoln.K],2,nTime);
figure(4); clf;
subplot(2,1,1); hold on;
plot(tSol,K(1,:),'k.','MarkerSize',15)
plot(tSol,ssK(1,:),'bo','MarkerSize',10)
legend('finite','infinite')
xlim(tSpan);
ylabel('Kx')
title('Compare finite vs infinite horizon gains')
subplot(2,1,2); hold on;
plot(tSol,K(2,:),'k.','MarkerSize',15)
plot(tSol,ssK(2,:),'bo','MarkerSize',10)
legend('finite','infinite')
xlim(tSpan);
ylabel('Ky')
xlabel('t')

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Lyapunov Stability Analysis                         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

% Let V(x) = x'Sx, where S is the cost-to-go matrix from LQR solution:

%%%% TODO %%%%

% Plan:
%
%   1) Represent the S matrix, which is time-varying, as a polynomial.
%
%   2) Compute the derivative of S wrt time.
%
%   3) Try the Lyapunov function:
%           V(x,t) = x'*S(t)*x     % x = perturbation from nominal
%           dV(x,t) = 2*x*f(x,t)*S(t)*dS(t)   %Chain Rule  (check math...)
%
%   4) Verify stability by checking signs of lyapunov function and
%      derivative?
%
%   5) Not sure if this whole approach is correct - check against papers...
%

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Run a bunch of simulations                      %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

controller = @(t,z)stabilizingController(t,z(1,:),z(2,:),xFit,yFit,uFit,kxFit,kyFit);
userFunc = @(t,z)rhs(t,z,controller(t,z));
zNom = [xSol(1);ySol(1)];

xSim = zeros(nTime,nSim);
ySim = zeros(nTime,nSim);
uSim = zeros(nTime,nSim);
failFlag = false(nSim,1);
for i=1:nSim
    try
        z0 = zNom + testPerturbation*randn(2,1);
        [~, yout] = ode45(userFunc,tSol,z0);
        xSim(:,i) = yout(:,1);
        ySim(:,i) = yout(:,2);
        uSim(:,i) = controller(tSol,yout');
    catch ME
        xSim(:,i) = z0(1);
        ySim(:,i) = z0(2);
        uSim(:,i) = 0;
        failFlag(i) = true;
    end
end

% Plot the simulated trajectories:
figure(1); hold on;  %Plot on top of dynamics vector field
for i=1:nSim
    plot(xSim(:,i),ySim(:,i),'k-','LineWidth',1);
end

%%%% Plot trajectories against time:
figure(5); clf;

for i=1:nSim
    if ~failFlag(i)
        subplot(3,1,1); hold on;
        plot(tSol,xSim(:,i),'k-','LineWidth',1)
        subplot(3,1,2); hold on;
        plot(tSol,ySim(:,i),'k-','LineWidth',1)
        subplot(3,1,3); hold on;
        plot(tSol,uSim(:,i),'k-','LineWidth',1)
    end
end

subplot(3,1,1); hold on;
plot(tSol,xSol,'r-','LineWidth',3)
xlim(tSpan);
ylabel('x')
title('Compare simulations against reference trajectory')
subplot(3,1,2); hold on;
plot(tSol,ySol,'r-','LineWidth',3)
xlim(tSpan);
ylabel('y')
subplot(3,1,3); hold on;
plot(tSol,uSol,'r-','LineWidth',3)
xlim(tSpan);
ylabel('u')
xlabel('t')

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%    Computate the backward reachable set   --  Convex Hull Method        %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

nSamples = 500;  %How many simulations to run?
nPointsPerSample = 10;  %Used for computing reachable set

% What is "close enough" to the final state?
xRadFinal = 0.015;
yRadFinal = 0.015;

%Compute a random set of initial conditions
IC = randDisk(nSamples,[xRadFinal,yRadFinal]);

%Simulation runs backwards in time
tSpanFlip = linspace(tSpan(2),tSpan(1),nPointsPerSample);   %Run system backwards in time

%Create a data structure for storing the point cloud:
ReachCvx.t = ones(nSamples,1)*tSpanFlip;  % Time
ReachCvx.x = zeros(nSamples,nPointsPerSample);  % x position
ReachCvx.y = zeros(nSamples,nPointsPerSample);  % y position
ReachCvx.c(nPointsPerSample).x = []; %x contour
ReachCvx.c(nPointsPerSample).y = []; %x contour

%Run the simulation:
for i=1:nSamples
    idxLow = (i-1)*nPointsPerSample+1;
    idxUpp= i*nPointsPerSample;
    [~, yout] = ode45(userFunc,tSpanFlip,IC(:,i));
    ReachCvx.x(i,:) = yout(:,1);
    ReachCvx.y(i,:) = yout(:,2);
end

%Compute the convex hull at each time step
for i=1:nPointsPerSample
    k = convhull(ReachCvx.x(:,i),ReachCvx.y(:,i));
    ReachCvx.c(i).t = ReachCvx.t(k,i);
    ReachCvx.c(i).x = ReachCvx.x(k,i);
    ReachCvx.c(i).y = ReachCvx.y(k,i);
end

%Plot the reachable sets against time:
figure(6); clf; hold on;
colors = jet(nPointsPerSample);
for i=1:nPointsPerSample
    %Convex Hull
    plot3(ReachCvx.c(i).x, ReachCvx.c(i).y, ReachCvx.c(i).t,'color',colors(i,:));
    %Point Cloud:
    plot3(ReachCvx.x(:,i), ReachCvx.y(:,i), ReachCvx.t(:,i),'.','MarkerSize',10,'color',colors(i,:));
end

%plot the nominal trajectory
plot3(xSol,ySol,tSol,'k-','LineWidth',3);

title('Convex Hull approximation of the reachable set');
xlabel('x'); ylabel('y'); zlabel('time');
view(3);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%    Computate the backward reachable set   --  Front Propagation         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

nSamples = 250;  %How many simulations to run?
nPointsPerSample = 10;  %Used for computing reachable set

% What is "close enough" to the final state?
xRadFinal = 0.015;
yRadFinal = 0.015;

%Compute initial conditions on the circle around the final set:
angle = linspace(0,2*pi,nSamples+1); angle(end) = [];
IC = [cos(angle); sin(angle)];
IC(1,:) = IC(1,:)*xRadFinal;
IC(2,:) = IC(2,:)*yRadFinal;

%Simulation runs backwards in time
tSpanFlip = linspace(tSpan(2),tSpan(1),nPointsPerSample);   %Run system backwards in time

%Create a data structure for storing the point cloud:
ReachFrt.t = ones(nSamples,1)*tSpanFlip;  % Time
ReachFrt.x = zeros(nSamples,nPointsPerSample);  % x position
ReachFrt.y = zeros(nSamples,nPointsPerSample);  % y position

%Run the simulation:
for i=1:nSamples
    idxLow = (i-1)*nPointsPerSample+1;
    idxUpp= i*nPointsPerSample;
    [~, yout] = ode45(userFunc,tSpanFlip,IC(:,i));
    ReachFrt.x(i,:) = yout(:,1);
    ReachFrt.y(i,:) = yout(:,2);
end

%Plot the reachable sets against time:
figure(7); clf; hold on;
colors = jet(nPointsPerSample);
idx = [1:size(ReachFrt.x,1),1];  %Make the contours wrap around
for i=1:nPointsPerSample
    %Front contour
    plot3(ReachFrt.x(idx,i), ReachFrt.y(idx,i), ReachFrt.t(idx,i),'color',colors(i,:));
    %Point Cloud:
    plot3(ReachFrt.x(:,i), ReachFrt.y(:,i), ReachFrt.t(:,i),'k.','MarkerSize',15);
    %Fill the central region:
    h = patch(ReachFrt.x(idx,i), ReachFrt.y(idx,i), ReachFrt.t(idx,i),colors(i,:));
end
plot3(xSol,ySol,tSol,'k-');

hp = findobj(gcf,'type','patch');
set(hp,'facealpha',0.7)
set(hp,'edgealpha',0.1)

title('Front Propagation approximation of the reachable set');
xlabel('x'); ylabel('y'); zlabel('time');
view(3);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%              Computation with the Level Set Toolbox                     %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

%This section of code requires the levelSetToolbox to operate:
%         http://www.cs.ubc.ca/~mitchell/ToolboxLS/
if exist('odeCFLset','file')
    %This section of code is still under development. As of 11/20/2014 the
    %following two scripts should run, but they are slow. It also seems
    %that there is some numerical dissapation, that is causing inaccurate
    %results for any reasonable level of discretization.

    %%%% reachableForward
    %%%% reachableForwardRel
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function inum=myRound(dnum)

if dnum>0
    inum=fix(dnum+0.5);
else
    inum=fix(dnum-0.5);
end

return

repo: GINav


text:
function [x, c, funVal, ValueL]=mcLogisticC(A, y, z, opts)
%
%%
% Function mcLeastR:
%      Logistic Loss for Multi-class (task) Learning
%             with the (group) L2/L1-norm Constraint
%
%% Problem
%
%  min  - sum_{il} weight_{il} log( p_{il} )
%  s.t.  sum_j ||x^j||_q <=z
%
%  p_{il}= 1 / (1+ exp(-y_i (x_i' * a_i + c_l) ) ) denotes the probability
%  weight_{il} is the weight for the i-th sample in the l-th classifier
%                is a m x k matrix
%  c_l is the intercept for the l-th classfier, and is a 1xk vector
%  x_i denotes the i-th column of x
%  x^j denotes the j-th row of x
%  a_i' denotes the i-th row of A
%
%  y_i (either 1 or -1) is the response
%
%  In this implementation, we assume weight_{il}=1/(mk)
%
%  When y are of binary values (1 and -1), this problem is the
%  multi-class classificaiton problem, which learns the projection
%  matrix x, by sharing information across k classification tasks (with the
%  Lq/L1-norm Regularization)
%
%  For the case that the multi tasks have separate data
%  matrix, please refer to the functions:
%           mtLeastR and mtLogisticR.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size m x k)
%  z -        L_1/L_q norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x k
%  c-         The obtained intercept if size 1 x k
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 19, 2010.
%
%% Related papers
%
% [1]  Jun Liu, Shuiwang Ji, and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, UAI, 2009
%
% [2]  Jun Liu, Lei Yuan, Songcan Chen and Jieping Ye, Multi-Task Feature Learning
%      Via Efficient L2,1-Norm Minimization, Technical Report ASU, 2009.
%
%% Related functions:
%
%  sll_opts, initFactor, pathSolutionLeast
%  mcLeastR, mcLeastC, mcLogisticR, ep21d
%
%%

%% Verify and initialize the parameters
%%
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
end

[m,n]=size(A);

k=size(y,2);

if (size(y,1) ~=m)
    error('\n Check the length of y!\n');
end

if (z<=0)
    error('\n z should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group Property

% Initialize q
if (~isfield(opts,'q'))
    q=2; opts.q=2;
else
    q=opts.q;
    if (q~=2)
        error('\n The current program only supports q=2');
    end
end

%% Starting point initialization

% % % compute AT y
% % if (opts.nFlag==0)
% %     ATy =A'* y;
% % elseif (opts.nFlag==1)
% %     ATy= A'* y -  mu' * sum(y, 1);    ATy=ATy./repmat(nu, 1, k);
% % else
% %     invNu=y./repmat(nu, 1, k);        ATy=A'*invNu-mu' * sum(invNu, 1);
% % end

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(p_flag,1);               % the total number of the positive samples
m2=m-m1;                        % the total number of the negative samples
% m1 and m2 are 1 x k vectors

% initialize a starting point
if opts.init==2
    x=zeros(n,k); c=zeros(1,k);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if ( size(x,1)~=n && size(x,2)~=k )
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,k);
    end

    if isfield(opts,'c0')
        c=opts.c0;

        if ( length(c)~=k )
            error('\n Check the input .c0');
        end
    else
        c=log(m1./m2);
    end
end

% compute Ax= A * x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
end

%% The main program

%% The Armijo Goldstein line search schemes

if (opts.lFlag==0)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/(m*k); % the intial guess of the Lipschitz continuous gradient

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;         ccp=zeros(1,k);

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ repmat(sc,m,1));

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y) * (1 - prob)
        b= -y.*(1-prob) / (m * k);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g= A'* b -  mu' * sum(b, 1);  g=g./repmat(nu, 1, k);
        else
            invNu=b./repmat(nu, 1, k);        g=A'*invNu-mu' * sum(invNu, 1);
        end

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L; c= sc- gc/L;

            % Lq/L1-norm regularized projection
            [x, lambda, zf_step]=ep21d(v, n, k, z, lambda0);
            lambda0=lambda;

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute Ax= A * x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./repmat(nu, 1, k);  mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./repmat(nu, 1, k);
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ repmat(c,m,1));

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

            r_sum=(norm(v,'fro')^2 + norm(c-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (c-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc> )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;  ccp=c-cp;

        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% adpative line search

if (opts.lFlag==1)

    lambda0=0;
    % a guess of the root

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/(m*k); % the intial guess of the Lipschitz continuous gradient

    gamma=1;
    % we shall set the value of gamma = L,
    % and L is appropriate for the starting point

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,k);
    cp=c;         ccp=zeros(1,k);

    for iterStep=1:opts.maxIter
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp; sc=c + beta* ccp;
                As=Ax + beta* (Ax-Axp);
            else
                alpha= (-1+ sqrt(5)) / 2;  beta=0;
                s=x; sc= c;
                As=Ax;
            end

            % aa= - diag(y) * (A * s + sc)
            aa=- y.*(As+ repmat(sc,m,1));

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_s= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

            % compute prob=[p_1;p_2;...;p_m]
            prob=1./( 1+ exp(aa) );

            % b= - diag(y) * (1 - prob)
            b= -y.*(1-prob) / (m * k);

            gc=sum(b); % the gradient of c

            % compute g= AT b, the gradient of x
            if (opts.nFlag==0)
                g=A'*b;
            elseif (opts.nFlag==1)
                g= A'* b -  mu' * sum(b, 1);  g=g./repmat(nu, 1, k);
            else
                invNu=b./repmat(nu, 1, k);        g=A'*invNu-mu' * sum(invNu, 1);
            end

            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L; cnew= sc- gc/L;

            % projection
            [xnew, lambda, zf_step]=ep21d(v, n, k, z, lambda0);
            lambda0=lambda;

            v=xnew-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A * xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Axnew+ repmat(cnew,m,1));

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= sum(sum ( log( exp(-bb) +  exp(aa-bb) ) + bb ) ) / (m * k);

            r_sum=(norm(v,'fro')^2 + norm(cnew-sc,2)^2) / 2;
            l_sum=fun_x - fun_s - sum(sum(v.* g)) - (cnew-sc)* gc';

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + <v, g> + <c ,gc>
            %                           + L/2 * (<v,v> + <c-sc,c-sc> )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;
        % store values for L

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;     Axp=Ax;
        x=xnew;   Ax=Axnew;
        cp=c; c=cnew;
        % update x and Ax with xnew and Axnew

        xxp=x-xp;  ccp=c-cp;

        funVal(iterStep)=fun_x;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=norm(xp,'fro');    norm_xxp=norm(xxp,'fro');
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

author: thomaskuestner


text:
%
%   FILE NAME: getMetsFromBiGG.m
%
%   PURPOSE: Generate data structure for BiGG metabolites
%

% Move to the target folder
cd('/Users/haowa/Box Sync/HMR3/BiGG');

% Load the text format BiGG reactions, after fixing errors with M01870
T=readtable('bigg_models_metabolites_20180424.txt','ReadVariableNames',1);
BiGGMets=table2struct(T,'ToScalar',true);

% Rename the fields according to RAVEN specification
BiGGMets.mets=BiGGMets.bigg_id;
BiGGMets.metNames=BiGGMets.name;
BiGGMets.universalID=BiGGMets.universal_bigg_id;

% Add oldids field as cell array based on info from old_bigg_ids
num=numel(BiGGMets.bigg_id);
BiGGMets.oldids=cell(num,1);
BiGGMets.oldids(:)={''};
% Add MNX field
metMNXID=regexp(BiGGMets.database_links,'MNXM\d+','match');
BiGGMets.metMNXID=cell(num,1);
BiGGMets.metMNXID(:)={''};
count=0;
for i=1:num
		if ~isempty(BiGGMets.old_bigg_ids{i})
				% Remove existing BiGG ids and universal ids from oldids
				BiGGMets.oldids{i}=setdiff(transpose(strsplit(BiGGMets.old_bigg_ids{i},'; ')),BiGGMets.mets{i});
				BiGGMets.oldids{i}=setdiff(BiGGMets.oldids{i},BiGGMets.universalID{i});
		end
		if ~isempty(metMNXID{i})
				if numel(metMNXID{i})==1
						BiGGMets.metMNXID{i}=metMNXID{i}{1};
				else
						BiGGMets.metMNXID{i}=metMNXID{i};
						count=count+1;
				end
		end
end
%count=0
numel(find(~cellfun(@isempty,BiGGMets.metMNXID)))  %ans = 9915

% Remove some fields
BiGGMets=rmfield(BiGGMets,{'bigg_id','name','universal_bigg_id','old_bigg_ids'});

save('BiGGMets.mat','BiGGMets');

author: SysBioChalmers

repo: Human-GEM


text:
function [kJ] = Btu2kJ(Btu)
% Convert energy or work from British thermal units to kilojoules.
% Chad A. Greene 2012
kJ = Btu*1.0550559 ;


text:
function Plugin_Timer

global VVar
global VCtl

if VVar.TRCount==1
    tic;
    set(VCtl.h.TimeWait_text,'String', ['Est. Time Left:  ' 'Estimating...']);
    pause(0.0001);
%     set(VCtl.h.TimeBar_text,'String', '0%');
else
    elptime=toc;
    tic;
    lefttime=double(VCtl.TRNum-VVar.TRCount+1)*elptime;
    % Convert seconds to other units
    d = floor(lefttime/86400); % Days
    lefttime = lefttime - d*86400;
    h = floor(lefttime/3600); % Hours
    lefttime = lefttime - h*3600;
    m = floor(lefttime/60); % Minutes
    lefttime = lefttime - m*60;
    s = floor(lefttime); % Seconds
    pause(0.0001);
    set(VCtl.h.TimeWait_text,'String', ['Est. Time Left:  ' num2str(h) ' : ' num2str(m) ' : ' num2str(s)]);
%     set(VCtl.h.TimeBar_text,'String',[num2str((double(VVar.TRCount)/double(VCtl.TRNum))*100, '%02.0f') '%']);

    DoUpdateBar(VCtl.h.TimeBar_axes,double(VVar.TRCount),double(VCtl.TRNum));
    pause(0.001);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function eThrottle = throttleefficiency(h,M,throttle,assumptions)
% Change in gas turbine efficiency as a function of throttle setting. Note
% that for this function throttle is shaft Preq/Pavail, NOT Treq/Tavail.
%
%   relativeEfficiency = throttleefficiency(h,M,throttle,assumptions)
%
%   See also CALCULATEPSFC.

%% Interpolation method
%{
tau = [0
0.185
0.35
0.5
0.65
0.8
1];

e = [0.56
0.746
0.86
0.947
0.988
1
0.988];

eThrottle = interp1(tau,e,throttle,'linear');
%}

%% Quadratic curve fit from interpolation data
eThrottle =-0.6713*throttle.^2 + 1.0928*throttle + 0.5626;

%% Raymer method
%{
TSFCoverTSFCfullthrottle = .1./throttle+.24./throttle.^.8+...
    .66*throttle.^.8 + .1*M.*(1./throttle+throttle);
eThrottle = max(2-TSFCoverTSFCfullthrottle,0);
%}

end


text:
classdef PTKAirwayGrowingLabelledBySegment < PTKPlugin
    % PTKAirwayGrowingLabelledBySegment. Plugin to find the segmental bronchi in each lobe
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Growing bronchi<br>labelled by segment'
        ToolTip = 'Finds the segmental bronchi in each lobe'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'DoNothing'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = false
        Visibility = 'Developer'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            growing_centreline_tree = dataset.GetResult('PTKAirwayGrowing');
            growing_centreline_tree = growing_centreline_tree.Airways;

            segments_by_bronchus = dataset.GetResult('PTKSegmentsByNearestBronchus');

            labelled_segmental_tree = segments_by_bronchus.AirwaysBySegment.Trachea;

            growing_tree_with_segmental_labels = PTKMapSegmentalParameters(labelled_segmental_tree, growing_centreline_tree, reporting);

            results = [];
            results.StartBranches = growing_tree_with_segmental_labels;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function r8r8_print ( a1, a2, title )

%*****************************************************************************80
%
%% R8R8_PRINT prints an R8R8.
%
%  Discussion:
%
%    An R8R8 is simply a pair of R8's, thought of as a single item.
%
%    A format is used which suggests a coordinate pair:
%
%  Example:
%
%    Center : ( 1.23, 7.45 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A1, A2, the coordinates of the vector.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '  %s : ( %12g, %12g )\n', title, a1, a2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function DEM_demo_contact_lens
% This demo illustrates tracking under the contact lens problem:
% The contact lens refers to the non-Gaussian uncertainty induced by
% nonlinear measurements. Here it is illustrated in terms of tracking the
% motion of a target in Cartesian coordinates, given the distance to target
% (range) and direction as measurements. The problem is to accumulate
% information over time about the target location under random fluctuations
% on the velocity (technically this is a constant acceleration model).
% Comparative evaluations are made with Extended Kalman filtering.
%
% See: X. Tian, Y. Bar-Shalom, Coordinate Conversion and Tracking for
% Very Long Range Radars. IEEE Transactions on Aerospace and Electronic
% Systems, AES-45(3):1073-1088, July 2009.
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_contact_lens.m 7679 2019-10-24 15:54:07Z spm $

% non-linear generative model
%==========================================================================

% The problem: states = [x(1) x(2)]; causes = [v(1) v(2)]
%--------------------------------------------------------------------------
x        = [1e4 1e4 -8 -2]';   % initial states: location = 10km,  10km
                               %                 velocity = -8m/s, -2m/s
V        = [1e-2 4];           % observation precision (inverse variance)
                               % standard deviation of range: 1/sqrt(V(1)) = 10m
                               % standard deviation of angle: 1/sqrt(V(2)) = .5 mrad
s        = 1;                  % smoothness of fluctuations
w        = 8;                  % precision of fluctuations in motion
W        = [128 128 w w];      % standard deviation of velocity: 1/sqrt(w) = .3536 m/s^2

% preliminaries
%--------------------------------------------------------------------------
N        = 256;                % length of sequence
t        = 1:N;                % time (seconds)
M(1).E.s = s;                  % smoothness of fluctuations
M(1).E.n = 4;                  % order of generalised coordinates
M(1).E.K = 128;                % rate of generalized gradient descent

% Level 1: Hidden states = [x(1) x(2) v(1) v(2)];
%--------------------------------------------------------------------------
f        = '[x(3); x(4); 0; 0]';
g        = '[sqrt(x(1)^2 + x(2)^2); 1000*atan(x(2)/x(1))]';
M(1).f   = inline(f,'x','v','P');
M(1).g   = inline(g,'x','v','P');
M(1).x   = x;
M(1).V   = diag(V);            % precision of observation noise
M(1).W   = diag(W);            % precision of fluctuations on hidden states

% create data
%==========================================================================
DEM    = spm_DEM_generate(M,N);

% Comparative inversions (variants of generalised filtering)
%==========================================================================

% reset initial position and bearing
%--------------------------------------------------------------------------
DEM.M(1).x = [1e3; 1e3; 0; 0];

% DEM and EKF
%--------------------------------------------------------------------------
DEM     = spm_DEM(DEM);
[EKF S] = spm_ekf(DEM.M,DEM.Y);

spm_figure('Getwin','DEM');
spm_DEM_qU(DEM.qU,DEM.pU)

% show prediction errors
%==========================================================================
spm_figure('GetWin','Figure 1');

D(1,:) = sqrt(sum((DEM.pU.x{1}([1 2],:) - DEM.qU.x{1}([1 2],:)).^2));
D(2,:) = sqrt(sum((DEM.pU.x{1}([1 2],:) - EKF([1 2],:)).^2));
D(3,:) = sqrt(sum((DEM.pU.x{1}([3 4],:) - DEM.qU.x{1}([3 4],:)).^2));
D(4,:) = sqrt(sum((DEM.pU.x{1}([3 4],:) - EKF([3 4],:)).^2));

for i = 1:N
    E         = DEM.pU.x{1}(:,i) - DEM.qU.x{1}(:,i);
    NEES(1,i) = E'*spm_inv(DEM.qU.S{i})*E;
    E         = DEM.pU.x{1}(:,i) - EKF(:,i);
    NEES(2,i) = E'*spm_inv(S{i})*E;
end

% plot errors
%--------------------------------------------------------------------------
subplot(2,2,1)
plot(t,log(D(1:2,:)))
title('log(location error (m))','Fontsize',16)
xlabel('time (secs)')
ylabel('log distance from target')
axis square
legend('DEM','EKF')

subplot(2,2,2)
plot(t,log(D(3:4,:)))
title('log(speed error (m/s))','Fontsize',16)
xlabel('time (secs)')
ylabel('log distance from target')
axis square
legend('DEM','EKF')

% plot trajectories
%--------------------------------------------------------------------------
subplot(2,2,3)
plot(DEM.pU.x{1}(1,:),DEM.pU.x{1}(2,:)     ),hold on
plot(DEM.qU.x{1}(1,:),DEM.qU.x{1}(2,:),'-.'),hold on
plot(EKF(1,:),EKF(2,:),':'),hold off
title('location (m)','Fontsize',16)
axis([8 12 8 12]*1000)
xlabel('x')
ylabel('y')
axis square
legend('true','DEM','EKF')

% plot NEES
%--------------------------------------------------------------------------
subplot(2,2,4)
plot(t,log(NEES))
title('log(NEES)','Fontsize',16)
xlabel('time (secs)')
axis square
legend('DEM','EKF')
drawnow

return

% prediction errors as a function of K
%==========================================================================
K     = -0:1/2:12;
for i = 1:length(K)

    % DEM and EKF (linear)
    %----------------------------------------------------------------------
    DEM.M(1).E.K = exp(K(i));
    DEM          = spm_DEM(DEM);

    % mean square error
    %----------------------------------------------------------------------
    kse(i,1) = mean(mean((DEM.pU.x{1} - DEM.qU.x{1}).^2));

end

% show prediction errors
%==========================================================================
spm_DEM_qU(DEM.qU,DEM.pU)

% plot trajectories
%--------------------------------------------------------------------------
subplot(2,1,2)
plot(K,log(kse))
title('log(MSE)','Fontsize',16)
xlabel('log(K)','Fontsize',12)
axis square

author: spm


text:
% Poster Figures
% Make figures for posters
%
% Author: Derek Macklin, macklin@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Created: 9/13/2011
classdef PosterFigures
    methods (Static = true);
        function run(eventName, outputDir)
            import edu.stanford.covert.cell.sim.analysis.PosterFigures

            PosterFigures.(eventName)(outputDir);
        end

        function pioneer2011(outputDir)
            import edu.stanford.covert.cell.sim.analysis.SingleCell;
            import edu.stanford.covert.cell.sim.analysis.ChromosomeSpaceTimePlot;
            import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.DiskLogger;
            import edu.stanford.covert.cell.sim.analysis.ProcessMetaboliteUsage;

            simBatchDir = '2011_09_09_19_17_51';
            simNum = 1;
            [~, ~, sim] = SimulationDiskUtil.getSimulation([simBatchDir filesep num2str(simNum)]);

            %% Gene -> RNA -> Protein -> Growth Correlation
            [~, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            SingleCell.growth(sim, simBatchDir, simNum, figHandle)
            saveas(figHandle, [outputDir filesep 'Growth.eps'], 'epsc');
            close(figHandle);

            %% Chromosomal stuff
            [simDir, ~, sim] = SimulationDiskUtil.getSimulation([simBatchDir filesep num2str(simNum)]);
            c = sim.state('Chromosome');

            simIdx = num2str(SimulationDiskUtil.getSimulationIndex(simDir));
            simTimeStamp = SimulationDiskUtil.getSimulationTimeStamp(simDir);
            stateNames = {'ProteinComplex' 'counts'
                'Chromosome' 'complexBoundSites'
                'Chromosome' 'monomerBoundSites'
                'Chromosome' 'polymerizedRegions'
                'RNAPolymerase' 'positionStrands'
                'RNAPolymerase' 'states'
                'RNAPolymerase' 'transcriptionFactorBindingProbFoldChange'
                'Time' 'values'};
            states = DiskLogger.load(simDir, stateNames, [], [], [], 'extract');

            height = 1000;
            width = 1000;

            RNAPSparseMat = ChromosomeSpaceTimePlot.makeRNAPSparseMat(states, sim);

            [ax, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            [densityBEMatrix, maxTime, maxSpace] = ChromosomeSpaceTimePlot.makeSpaceTimeDensityPlots(...
                ax, sim, states, RNAPSparseMat, ...
                ChromosomeSpaceTimePlot.chromProteins, [], width, height, simTimeStamp, simIdx);
            close(figHandle);
            [ax, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            set(figHandle, 'Color', [1 1 1]);
            ChromosomeSpaceTimePlot.plotCircularDensity(...
                ax, densityBEMatrix, maxTime, sim, false, ...
                sprintf('Circular Chromosome Density Plot\nSimulation: %s #%s', simTimeStamp, simIdx));
            set(ax, 'Position', [0.2 0.2 0.6 0.6])
            saveas(figHandle, [outputDir filesep 'CircularDensity.eps'], 'epsc');
            close(figHandle);

            selectedTime = 4 * 3600;
            [ax, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            ChromosomeSpaceTimePlot.makeRingPlot(...
                ax, sim, states, RNAPSparseMat, ...
                ChromosomeSpaceTimePlot.chromProteins, [], selectedTime, ...
                maxSpace);
            saveas(figHandle, [outputDir filesep 'RingPlot.eps'], 'epsc');
            close(figHandle);

            [ax, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            whatToShow = {'Topoisomerase IV', 'DnaA Complex', 'LuxR', 'HTH Regulator', 'Ferric uptake repressor', 'Helicase'};
            ChromosomeSpaceTimePlot.plotSpaceTimeOverlay(...
                ax, sim, states, RNAPSparseMat, ...
                ChromosomeSpaceTimePlot.chromProteins, whatToShow, ...
                sprintf('Chromosome 1\nSimulation: %s #%s', simTimeStamp, simIdx));
            set(figHandle, 'Renderer', 'Painters')
            saveas(figHandle, [outputDir filesep 'SpaceTimeOverlay.eps'], 'epsc');
            close(figHandle);
            clear('states', 'stateNames', 'RNAPSparseMat', 'whatToShow', 'height', 'width');
%
            %% Metabolite Usages
            wIDsToShow = {'ATP', 'GTP'};
            desiredInfo = 'Usages';
            [ax, figHandle] = edu.stanford.covert.cell.sim.util.PlotUtil.newAxesHandle();
            ProcessMetaboliteUsage.processSingleCellMetaboliteHandler(wIDsToShow, 'Usages', figHandle, '2011_09_09_19_17_51', simNum)
            saveas(figHandle, [outputDir filesep 'ProcessUsages.eps'], 'epsc');
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [element2edges, edge2nodes]=getEdges(elements)
%function: [element2edges, edge2nodes]=edge_numbering(elements)
%requires: deleterepeatedrows
%generates edges of (triangular) triangulation defined in elements
%elements is matrix, whose rows contain numbers of its element nodes
%element2edges returns edges numbers of each triangular element
%edge2nodes returns two node numbers of each edge
%example in 2D: [element2edges, edge2nodes]=getEdges([1 2 3; 2 4 3])
%example in 3D: [element2edges, edge2nodes]=getEdges([1 2 3 4; 1 2 3 5; 1 2 4 6])

%2D case
if (size(elements,2)==3)
    %extracts sets of edges
    edges1=elements(:,[2 3]);
    edges2=elements(:,[3 1]);
    edges3=elements(:,[1 2]);

    %as sets of their nodes (vertices)
    vertices=zeros(size(elements,1)*3,2);
    vertices(1:3:end,:)=edges1;
    vertices(2:3:end,:)=edges2;
    vertices(3:3:end,:)=edges3;

    %repeated sets of nodes (joint edges) are eliminated
    [edge2nodes,element2edges]=deleterepeatedrows(vertices);
    element2edges=reshape(element2edges,3,size(elements,1))';
end

%3D case
if (size(elements,2)==4)
    %extracts sets of edges
    edges1=elements(:,[2 3]);
    edges2=elements(:,[3 1]);
    edges3=elements(:,[1 2]);
    edges4=elements(:,[3 4]);
    edges5=elements(:,[1 4]);
    edges6=elements(:,[2 4]);

    %as sets of their nodes (vertices)
    vertices=zeros(size(elements,1)*6,2);
    vertices(1:6:end,:)=edges1;
    vertices(2:6:end,:)=edges2;
    vertices(3:6:end,:)=edges3;
    vertices(4:6:end,:)=edges4;
    vertices(5:6:end,:)=edges5;
    vertices(6:6:end,:)=edges6;

    %repeated sets of nodes (joint edges) are eliminated
    [edge2nodes,element2edges]=deleterepeatedrows(vertices);
    element2edges=reshape(element2edges,6,size(elements,1))';
end

author: Sable


text:
function [rho pval] = circ_corrcl(alpha, x)
%
% [rho pval ts] = circ_corrcc(alpha, x)
%   Correlation coefficient between one circular and one linear random
%   variable.
%
%   Input:
%     alpha   sample of angles in radians
%     x       sample of linear random variable
%
%   Output:
%     rho     correlation coefficient
%     pval    p-value
%
% References:
%     Biostatistical Analysis, J. H. Zar, p. 651
%
% PHB 6/7/2008
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html

if size(alpha,2) > size(alpha,1)
	alpha = alpha';
end

if size(x,2) > size(x,1)
	x = x';
end

if length(alpha)~=length(x)
  error('Input dimensions do not match.')
end

n = length(alpha);

% compute correlation coefficent for sin and cos independently
rxs = corr(x,sin(alpha));
rxc = corr(x,cos(alpha));
rcs = corr(sin(alpha),cos(alpha));

% compute angular-linear correlation (equ. 27.47)
rho = sqrt((rxc^2 + rxs^2 - 2*rxc*rxs*rcs)/(1-rcs^2));

% compute pvalue
pval = 1 - chi2cdf(n*rho^2,2);

author: kavli-ntnu


text:
function [elem,HB,newHB,tree]= uniformcoarsen3(elem,HB)
%% UNIFORMCOARSEN3 uniform coarsening in 3-D
%
% [elem,HB,newHB] = uniformcarsen3(elem,HB) remove all good-to-coarsen nodes.
% It is mainly used to get multilevel decomposition in multigrid methods.
% The input matrix elem stands for the fine mesh and the output one for the
% coarse mesh, whose nodal indices are shifted and thus different with the
% fine mesh.
%
% Unlike the two dimensional case, an additional matrix HB is introduced.
% The local index of HB is 2 -- 1 -- 3, i.e., HB(:,1) is the middle point
% of the edge formed by HB(:,2:3). HB(:,4) is used to store the generation
% of the node HB(:,1).
%
% In the output, newHB(:,1) are all removed nodes and newHB(:,2:3) are two
% neighboring nodes whose indices are in the coarse mesh.
%
% See also: coarsen3, bisect3, uniformcoarsen, MGP1, MGP2, MG3P1, MG3P2
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Find good-to-coarsen nodes
N = max(elem(:));
NT = size(elem,1);
generation = zeros(N,1);
generation(HB(:,1)) = HB(:,4);
valence = accumarray(elem(:),ones(4*NT,1),[N 1]);
valenceNew = accumarray(elem(:,4),ones(NT,1), [N 1]); % for newest nodes only
isGoodNode = (valence == valenceNew) & (generation > 0);
if ~any(isGoodNode)
    newHB = []; tree = [];
    return
else
    nGoodNode = sum(isGoodNode);
    newHB = zeros(nGoodNode,3);
    newHB(:,1) = find(isGoodNode);
    newHB(:,2) = HB(newHB(:,1),2);
    newHB(:,3) = HB(newHB(:,1),3);
end

%% Remove good-to-coarsen nodes
t = find(isGoodNode(elem(:,4)));     % elements containing good nodes
leftNode = HB(elem(t,4),2);          % element-wise left node of good nodes
idx = (elem(t,1)==leftNode) | (elem(t,2)==leftNode) | (elem(t,3)==leftNode);
tl = t(idx);
tr = t(~idx);
Nr = length(tl); % number of refined elements
tree = zeros(Nr,3);
% sort tl and tr by the common vertices such that tl and tr matches
if nargout == 4
    leftNode = HB(elem(tl,4),2);
    rightNode = HB(elem(tr,4),3);
    temptl = zeros(Nr,2);
    idx = (elem(tl,1) == leftNode);
    temptl(idx,:) = elem(tl(idx),[2 3]);
    idx = (elem(tl,2) == leftNode);
    temptl(idx,:) = elem(tl(idx),[1 3]);
    idx = (elem(tl,3) == leftNode);
    temptl(idx,:) = elem(tl(idx),[1 2]);
    temptl = sort(temptl,2);
    temptr = zeros(Nr,2);
    idx = (elem(tr,1) == rightNode);
    temptr(idx,:) = elem(tr(idx),[2 3]);
    idx = (elem(tr,2) == rightNode);
    temptr(idx,:) = elem(tr(idx),[1 3]);
    idx = (elem(tr,3) == rightNode);
    temptr(idx,:) = elem(tr(idx),[1 2]);
    temptr = sort(temptr,2);
    [temptl, Il] = sortrows([temptl elem(tl,4)]);
    [temptr, Ir] = sortrows([temptr elem(tr,4)]);
    tl = tl(Il);
    tr = tr(Ir);
    tree(:,3) = tr;
    tree(:,1) = tl;
    tree(:,2) = tl;
end
% coarsen tl to t
elem(tl,4) = HB(elem(tl,4),3);       % replace new node by right node

%% Sort element nodes by generations
% the newest node is the node with maxmum generation
if (length(tl)==1)
    idx = max(transpose(generation(elem(t1,:))));
else
    [tempvar,idx] = max(generation(elem(tl,:)),[],2);  %#ok<*ASGLU>
end
elem(tl((idx==1)),1:4) = elem(tl((idx==1)),[2 4 3 1]);
elem(tl((idx==2)),1:4) = elem(tl((idx==2)),[3 4 1 2]);
elem(tl((idx==3)),1:4) = elem(tl((idx==3)),[4 2 1 3]);

%% Clean and shift index
elem(tr,:) = [];
inCoarse = true(NT,1);
inCoarse(tr) = false;
elemidxMap = zeros(NT,1);
elemidxMap(inCoarse) = 1:size(elem,1);
if nargout == 4
    tree(:,1) = elemidxMap(tree(:,1));
end
HB(isGoodNode,:) = [];
indexMap = zeros(N,1);
indexMap(~isGoodNode)= 1:N-nGoodNode;
elem = indexMap(elem);
HB(:,1:3) = indexMap(HB(:,1:3));
newHB(:,2:3) = indexMap(newHB(:,2:3));


text:
function varargout = arg_toworkspace(s)
% Copy the arguments in the given Struct into the workspace of the calling function.
% Outputs... = arg_toworkspace(Struct)
%
% This function can be used to populate the function's workspace with fields of a struct, and
% optionally allows to assign the first k fields of the struct to explicitly named output variables.
% For all other fields variables with the same name as the field are created implicitly in the
% function's workspace by arg_toworkspace.
%
% Using explicit names and return values allows to work around situations where MATLAB "doesn't
% realize" that a given variable has been written to the workspace and accidentally looks up the
% identifier from the MATLAB path (e.g., a function with the same name). The latter issue is
% explained in online posts about "poofing" variables into the workspace.
%
% In:
%   Struct : an argument structure, as produced by arg_define
%
% Out:
%   Outputs... : optional output arguments that shall receive the first k fields in the Struct
%                all remaining fields are written into the workspace by their name.
%
% Examples:
%   function myfunction(varargin)
%   opts = arg_define(varargin, ...
%               arg('myarg1',...), ...
%               arg('myarg2',...), ...
%               arg('myarg3',...));
%
%   % copy all args into the workspace
%   arg_toworkspace(opts);
%
%   % alternative: accept the first two arguments (myarg1,myarg2) into variables named myfirstarg,
%   % and mysecondarg, copy the rest into the workspace by name (myarg3)
%   [myfirstarg,mysecondarg] = arg_toworkspace(opts);
%
% See also:
%   arg_define
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-09-24

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% 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 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

% assign the first k arguments to the given output variables
if nargout > 0
    fields = struct2cell(s);
    varargout(1:nargout) = fields(1:nargout);
end

% place all remaining variables in the workspace by their name
fnames = fieldnames(s);
for fn=fnames(nargout+1:end)'
    assignin('caller',fn{1},s.(fn{1})); end

repo: REST


text:
function out = ST_FitPolynomial(y,k)
% ST_FitPolynomial   Goodness of a polynomial fit to a time series
%
% Usually kind of a stupid thing to do with a time series, but it's sometimes
% somehow informative for time series with large trends.
%
%---INPUTS:
% y, the input time series.
% k, the order of the polynomial to fit to y.
%
%---OUTPUT:
% RMS error of the fit.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

doPlot = 0; % Plot stuff to screen

if nargin < 2 || isempty(k)
    k = 1; % Linear by default
end

N = length(y); % Length of the time series (number of samples)
t = (1:N)'; % Get a range for the time axis for time series y

% ------------------------------------------------------------------------------
% Fit a polynomial to the time series
% ------------------------------------------------------------------------------
% Supress the (valid!) warning from stupidly fitting a polynomial to a time series...
warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');
cf = polyfit(t,y,k);
warning('on','MATLAB:polyfit:RepeatedPointsOrRescale');

f = polyval(cf,t);
out = mean((y-f).^2); % mean RMS ERROR OF FIT

% ------------------------------------------------------------------------------
% Plot
% ------------------------------------------------------------------------------
if doPlot
    n = 10;
    errs = zeros(n,1);
    x = 1:length(y);
    for i = 1:n
        cf = polyfit(x,y',i);
        f = polyval(cf,x);
        errs(i) = sum((y'-f).^2);
    end
    f = figure('color','w'); hold on;
    plot(f,'k')
    plot(errs);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa