text:
function y = wavelift(x, nlevel, wname)
%WAVELIFT: Multi-level discrete two-dimension wavelet transform
%based on lifting method.
%
% c = wavelift(x, nlevel, wname) performs the follows according to the
% value of nlevel:
%   nlevel > 0:   decomposes 2-dimension matrix x up to nlevel level;
%   nlevel < 0:   does the inverse transform to nlevel level;
%   nlevel = 0:   sets c equal to x;
%
% wname is name of wavelet used for DWT or IDWT. It can be omitted.
% If so, WAVELIFT use the default Cohen-Daubechies-Feauveau (CDF) 9/7
% wavelet, which is the name 'cdf97'.Currently, WAVELIFT only support
% two kind of wavelets, i.e. cdf97 and spline 5/3 with the name 'spl53'.
% However, aided with the organized lifting structure illustrated below,
% it is adaptive to other specific lifting realizations. The only thing
% needed in most cases is only to modify the structure L and the mode to
% indicate lossy or lossless compression.
%
% WAVELIFT call another function COLWAVELIFT to perform 1-D FWT based on
% lifting method. Deliberately organized lifting structure is provided
% to COLWAVELIFT as a major parameter.
%
% The lifting structure is organized as follows:
% L: 1-by-1 structure with two fields lamdaz and K.
%   K is two-element vector [K0, K1], which is the lifting gains.
%   lamdaz is 1-by-M structure if M lifting units are used.
%     lamdaz's two fields coeff and zorder denote the transfer function
%     of every lifting units lamda(Z)
%   e.g. for a wavelet transform with 3 lifting units as
%     lamda1 = a1+a2*z, lamda2 = b1+b2*z^-1, lamda3 = c1*z^(-1)+c2*z
%     and the lifting gains K0 and K1
%   L is to be organized as
%     lamdaz = struct('coeff', {[a1, a2], [b1, b2], [c1, c2], ...
%                    'zorder', {[ 0,  1], [ 0, -1], [-1, 1 ]} );
%     L = struct('lamdaz', lamdaz, 'K', [K0, K1]);
%
% You can test WAVELIFT with following lines:
%   x=imread('E:\study\jpeg2000\images\lena.tif');
%   % see the decomposition coefficients
%   y=wavelift(x, 1, 'spl53'); % using lossless spline 5/3 wavelet
%   figure; subplot(1,2,1); imshow(x); subplot(1,2,2); imshow(mat2gray(y))
%   % see the reconstruction precision
%   yy=wavelift(x, 5); % using lossy cdf 9/7 wavelet
%   ix=wavelift(yy, -5); % inverse
%   sum(sum((double(x)-ix).^2))
%
% Reference:
%   [1] D.S.Taubman et al., JPEC2000 Image Compression: F. S. & P.,
%       Chinese Edition, section 6.4, 6.5, 10.3 and 10.4
%   [2] Pascal Getreuer, waveletcdf97.m from Matlab file Exchange website
%
% Cantact information:
%   Email/MSN messenger: wangthth@hotmail.com
%
% Tianhui Wang at Beijing, China,  Aug 5, 2006
%                  Last Revision:  Aug 6, 2006

%----------------------- input arguments checking ----------------------%
error(nargchk(2, 3, nargin));
% default decomposition/forward lifting level
if nargin < 3
    wname = 'cdf97';
end
% check nlevel
if ~isreal(nlevel) || ~isnumeric(nlevel) || round(nlevel)~=nlevel
    error('WAVELIFT:InArgErr', ['The 2nd argument shall be ' ...
        'a real and numeric integer.']);
end
% check x
if ~isreal(x) || ~isnumeric(x) || (ndims(x) > 2)
    error('WAVELIFT:InArgErr', ['The first argument must' ...
        ' be a real, numeric 2-D or 1-D matrix.']);
end
if isinteger(x)
    x = double(x);
end
% check wname
if ~ischar(wname) || ~ismember(wname, {'cdf97', 'spl53'})
    error('WAVELIFT:InArgErr', ['The last argument must be a wavelet ' ...
        'name. \nCurrently only ''cdf97'' and ''spl53'' are supported.']);
end
%------------- forming lifting structure and lifting mode --------------%
switch wname
    case 'cdf97'
       lamdaz=struct('coeff',{[-1.5861343420693648,-1.5861343420693648],...
                              [-0.0529801185718856,-0.0529801185718856],...
                              [ 0.8829110755411875, 0.8829110755411875],...
                              [ 0.4435068520511142,0.4435068520511142]},...
                   'zorder', {[0 1], [0 -1], [0 1], [0 -1]});
       L=struct('lamdaz',lamdaz,'K',[1/1.230174104914, 1.230174104914/2]);
       % the line below is another version of K used by P.Getreuer[2]
       % L = struct('lamdaz', lamdaz, 'K', ...
       %    [1.1496043988602418, 1/1.1496043988602418]);
       mode='lossy';
    case 'spl53'
       lamdaz = struct('coeff',  {[-.5, -.5], [.25 .25]}, ...
                       'zorder', {[ 0,   1 ], [ 0, -1 ]});
       L = struct('lamdaz', lamdaz, 'K', [1, 1/2]);
       mode='lossless';
end
clear lamdaz;
% set y be x for the sake of lifting processes
% and also for the case of nlevel = zero
y = x;
%-----------  decomposition/forward lifting,  when nlevel > 0  ---------%
if nlevel > 0
    for i = 1:nlevel
        sx = size(x);
        % first lift all columns of x
        [temp0, temp1] = colwavelift(x, L, 'd', mode);
        % the inverse lifting for rows of temp0 and temp1 can be
        % performed simultaneously using 1-D column lifting process
        [temp0, temp1] = colwavelift([temp0; temp1]', L, 'd', mode);
        temp = [temp0', temp1'];
        % update coefficient matrix
        y(1:sx(1), 1:sx(2)) = temp;
        % replace x with temp upper left quarter for next level FWT
        x = temp(1:ceil(sx(1)/2), 1:ceil(sx(2)/2));
        % give a warning if nlevel is too large
        if size(x,1)<=1 && size(x,2)<=1 && i~=nlevel
            warning('WAVELIFT:InArgDegrade', ['Only decompose to ' ...
                num2str(i) '-level instead of ' num2str(nlevel) ...
                ', \nas the approximation coefficients at ' num2str(i) ...
                '-level has row or/and column of length 1.']);
            break
        end
    end
%------------  reconstruction/inverse lifting,  if nlevel < 0  ---------%
else
    sx = size(x);
    % reconstruction level
    nl = -nlevel;
    while sx(1)/2^nl<=1/2 && sx(2)/2^nl<=1/2,  nl = nl-1;  end
    if nl ~= -nlevel
        warning('WAVELIFT:InArgDegrade', ['Only reconstruct to ' ...
            num2str(nl) '-level instead of ' num2str(-nlevel) ...
            ', \n as the approximation coefficients at ' num2str(nl) ...
            '-level has row or/and column of length 1.']);
    end
    % 2-D reconstruction
    for i = 1 : nl
        % find the target LL block
        sTarget = ceil(sx/2^(nl-i));
        target = y(1:sTarget(1), 1:sTarget(2));
        % perform inverse lifting for all rows using column 1-D lifting
        sLL = ceil(sTarget/2);
        temp0 = target(:, 1: sLL(2));
        temp1 = target(:, sLL(2)+1:end);
        temp = colwavelift(temp0', temp1', L, 'r', mode);
        temp = temp';
        % with the upper half of temp being the even sequences and the
        % lower half being the odd, perform inverse lifting
        % simultaneously for all columns
        temp0 = temp(1: sLL(1), :);
        temp1 = temp(sLL(1)+1 :end, :);
        temp = colwavelift(temp0, temp1, L, 'r', mode);
        % update y with the new LL block
        y(1:sTarget(1), 1:sTarget(2)) = temp;
    end
end
% EOF


text:
%% Copyright (C) 2014-2016, 2019, 2022-2023 Colin B. Macdonald
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2016 Abhinav Tripathi
%% Copyright (C) 2017 NVS Abhilash
%% Copyright (C) 2020 Fernando Alvarruiz
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @defun mat_replace (@var{A}, @var{subs}, @var{rhs})
%% Private helper routine for setting symbolic array entries.
%%
%% @end defun

function z = mat_replace(A, subs, b)

  if (length (subs) == 1 && islogical (subs{1}))
    %% A(logical) = B
    subs{1} = find (subs{1});
    if isempty (subs{1})
      z = A;
      return;
    end
  end
  %% Check when b is []
  if (isempty(b))
    switch length(subs)
      case 1
        if strcmp(subs{1}, ':')
          z = sym([]);
          return
        end
        if (isempty (subs{1}))
          z = A;
          return
        end
          if rows(A) == 1
            z = delete_col(A, subs{1});
            return
          elseif columns(A) == 1
            z = delete_row(A, subs{1});
            return
          else
            z = sym([]);
            for i=1:A.size(2)
              z = [z subsref(A, substruct ('()', {':', i})).'];
            end
            z = subsasgn (z, substruct ('()', {subs{1}}), []);
            return
          end
      case 2
        if (isempty (subs{1}) || isempty (subs{2}))
          z = A;
          return
        end
        if strcmp(subs{1}, ':')
          if strcmp(subs{2}, ':')
            z = sym(zeros(0,columns(A)));
            return
          else
            z = delete_col(A, subs{2});
            return
          end
        elseif strcmp(subs{2}, ':')
          z = delete_row(A, subs{1});
          return
        else
          error('A null assignment can only have one non-colon index.'); % Standard octave error
        end
      otherwise
        error('Unexpected subs input')
    end
  end

  if (length(subs) == 1 && strcmp(subs{1}, ':') && length(b) == 1)
    z = pycall_sympy__ ('return ones(_ins[0], _ins[1])*_ins[2],', uint64(A.size(1)), uint64(A.size(2)), sym(b));
    return

  elseif (length(subs) == 1)
    % can use a single index to grow a vector, so we carefully deal with
    % vector vs linear index to matrix (not required for access)
    [n,m] = size(A);
    if (n == 0 || n == 1)
      c = subs{1};  r = ones(size(c));
    elseif (m == 1)
      r = subs{1};  c = ones(size(r));
    else
      % linear indices into 2D array
      % Octave 8 does not raise error from ind2sub so we do it ourselves
      sz = size (A);
      i = subs{1};
      if (i > prod (sz))
        error ('%d is out of bound %d (dimensions are %dx%d)\n', i, prod (sz), sz)
      end
      [r, c] = ind2sub (size (A), i);
      % keep all the indices in a row vector
      r = reshape (r, 1, []);
      c = reshape (c, 1, []);
    end

  elseif (length(subs) == 2)
    r = subs{1};
    c = subs{2};
    [n,m] = size(A);

    if (isnumeric(r) && ((isvector(r) || isempty(r))))
      % no-op
    elseif (strcmp(r,':'))
      r = 1:n;
    elseif (islogical(r) && isvector(r) && (length(r) == n))
      I = r;
      r = 1:n;  r = r(I);
    else
      error('unknown 2d-indexing in rows')
    end

    if (isnumeric(c) && ((isvector(c) || isempty(c))))
      % no-op
    elseif (strcmp(c,':'))
      c = 1:m;
    elseif (islogical(c) && isvector(c) && (length(c) == m))
      J = c;
      c = 1:m;  c = c(J);
    else
      error('unknown 2d-indexing in columns')
    end

    [r,c] = ndgrid(r,c);
    if ~(isscalar (b) || (isvector (r) && isvector (b)) || is_same_shape (r, b))
      % vectors may have diff orientations but if we have matrices then
      % they must have the same shape (Octave/Matlab do this for double)
      error('A(I,J,...) = X: dimensions mismatch')
    end
    r = r(:);
    c = c(:);
  else
    error('unknown indexing')
  end

  z = mat_rclist_asgn(A, r, c, b);

end

function z = delete_col(A, subs)
  if isscalar (A)
    z = sym(zeros (1, 0));
  else
    cmd = { 'A, subs = _ins'
            'if isinstance(subs, Integer):'
            '    A.col_del(subs - 1)'
            '    return A,'
            'for i in sorted(subs, reverse=True):'
            '    A.col_del(i - 1)'
            'return A,' };
    z = pycall_sympy__ (cmd, A, sym(subs));
  end
end

function z = delete_row(A, subs)
  if isscalar (A)
    % no test coverage: not sure how to hit this
    z = sym(zeros (0, 1));
  else
    cmd = { 'A, subs = _ins'
            'if isinstance(subs, Integer):'
            '    A.row_del(subs - 1)'
            '    return A,'
            'for i in sorted(subs, reverse=True):'
            '    A.row_del(i - 1)'
            'return A,' };
    z = pycall_sympy__ (cmd, A, sym(subs));
  end
end

repo: octsympy


text:
function p = eval_pdf_clg(X,Y,mu,Sigma,W)
% function p = eval_pdf_clg(X,Y,mu,Sigma,W)
%
% p(c,t) = N(Y(:,t); mu(:,c) + W(:,:,c)*X(:,t), Sigma(:,:,c))

[d T] = size(Y);
[d nc] = size(mu);
p = zeros(nc,T);
for c=1:nc
  denom = (2*pi)^(d/2)*sqrt(abs(det(Sigma(:,:,c))));
  M = repmat(mu(:,c), 1, T) + W(:,:,c)*X;
  mahal = sum(((Y-M)'*inv(Sigma(:,:,c))).*(Y-M)',2);
  p(c,:) = (exp(-0.5*mahal) / denom)';
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Test file for ADCHEBFUN norm related functions

function pass = test_norm

% List of trigonometric functions to test.
normFun1 = @(u) norm(u, 2);
normFun2 = @(u) norm(u, 2).^2;
normFun3 = @(u) 1./norm(u, 2);
normFun4 = @(u) 1./norm(u, 2).^2;
normFun5 = @(u) norm(u, 2) + norm(diff(u));
normFun5 = @(u) norm(u, 2).^2 + norm(diff(u)).^3;

funcList = {normFun1, normFun2, normFun3, normFun4, normFun5};

% Call the ADCHEBFUN testUnary() method to do the tests.
pass = adchebfun.testUnary(funcList);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [V,lagsused]=covvar(data,maxlag,method)
% Long-run covariance estimation using the VAR-based method of DenHaan and Levin (1996)
%
%   USAGE:
%     V=covvar(DATA)
%     [V,LAGSUSED]=covvar(DATA,MAXLAG,METHOD)
%
%   INPUTS:
%     DATA    - T by K vector of dependent data
%     MAXLAG  - Non-negative integer containing the maximum lag length to use.  If empty or not
%                 included, MAXLAG=min(1.2*floor(T^(1/3)),floor(T/K)) is used
%     METHOD  - An integer value 1, 2, 3, 4 or 5. 2 is the DEFAULT.
%                 1 - Use lags 1 to MAXLAG to compute the VAR based covariance
%                 2 - [DEFAULT] Use up to MAXLAG but select using SIC
%                 3 - Use up to MAXLAG but select using AIC
%                 4 - Use up to MAXLAG but select using SIC and a global search
%                 5 - Use up to MAXLAG but select using AIC and a global search
%
%   OUTPUTS:
%     V        - A K by K covariance matrix estimated using the VAR based esitmator
%     LAGSUSED - A row vector indicating the vlags used in computing the VAR-based long-run
%                  covariance.  If empty, no lags selected, and the estimated covariance is
%                  identical to the usual covariance.
%
%   COMMENTS:
%    If options 4 or 5 are used, all combinations of lags between none and MAXLAG are tried (e.g. [0
%    (constant)],[1],[2],[1 2],[3],[1 3],[2 3],[1 2 3] ...).  These two options are only practically
%    viable for smallish MAXLAG (~<10 to 20)  and small K

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 5/1/2007

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T,K]=size(data);
if nargin==1
    maxlag=min(floor(T/K),floor(1.2*T^(1/3)));
    method=2;
elseif nargin==2
    method=2;
end
if isempty(maxlag)
    maxlag=min(floor(T/K),floor(1.2*T^(1/3)));
end
if isempty(method)
    method=2;
end
if ~ismember(method,1:5)
    error('MATHOD must be a scalar between 1 and 5.')
end
if floor(maxlag)~=maxlag || maxlag<0 || maxlag>floor(T/K)
    error('MAXLAG must be a non-negative integer with MAXLAG<=floor(T/K).')
end
if ndims(data)>2
    error('DATA must be a T by K matrix of data.')
end
switch method
    case {1,2,4}
        ICtype=1;
    case {3,5}
        ICtype=2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up the indices to use
if method==1
    indices{1}=1:maxlag;
elseif method==2 || method==3
    indices=cell(maxlag+1,1);
    for i=1:maxlag
        indices{i}=1:i;
    end
    indices{maxlag+1}=[];
else
    indices=cov_VAR_dec2bin(maxlag);
end

% Need to use the same amount of 'y' data for each to make SIC comparable
indep=cell(K,1);
Y=zeros(T-maxlag,K);
X=zeros(T-maxlag,K*maxlag);

for i=1:K;
    [dep,indep{i}]=newlagmatrix(data(:,i),maxlag,0);
    Y(:,i)=dep;
end
index=1;
for j=1:maxlag
    for i=1:K;
        temp=indep{i};
        X(:,index)=temp(:,j);
        index=index+1;
    end
end
X=[ones(size(X,1),1) X];

% Now that the X and Y matrices are set up, I can loop over the indices to
% run the regressions.  I need to save the SIC/AIC
N=length(indices);
IC=zeros(N,1);
T2=(T-maxlag);
for i=1:N
    P=length(indices{i});
    if ~isempty(indices{i});
        cols=repmat(indices{i}-1,K,1)*K+repmat((0:K-1)',1,P)+2;
    else
        cols=[];
    end
    regressors=X(:,[1 cols(:)']);
    B=regressors\Y;
    e=Y-X(:,[1 cols(:)'])*B;
    covE=e'*e/T2;
    if ICtype==1
        IC(i)=log(det(covE))+log(T2)/T2*(P*K^2+K);
    else
        IC(i)=log(det(covE))+2/T2*(P*K^2+K);
    end
end

% Finally pick the best
[temp,ICpos]=min(IC);

% Use the maximum amount of data to estimate the covariance
if ~isempty(indices{ICpos})
    maxlag=max(indices{ICpos});
    Y=zeros(T-maxlag,K);
    X=zeros(T-maxlag,K*maxlag);
    for i=1:K;
        [dep,indep{i}]=newlagmatrix(data(:,i),maxlag,0);
        Y(:,i)=dep;
    end
    index=1;
    for j=1:maxlag
        for i=1:K;
            temp=indep{i};
            X(:,index)=temp(:,j);
            index=index+1;
        end
    end
    X=[ones(size(X,1),1) X];
    T2=(T-maxlag);
    i=ICpos;
    P=length(indices{i});
    if ~isempty(indices{i});
        cols=repmat(indices{i}-1,K,1)*K+repmat((0:K-1)',1,P)+2;
    else
        cols=[];
    end
    regressors=X(:,[1 cols(:)']);
    B=regressors\Y;
    e=Y-X(:,[1 cols(:)'])*B;
    covE=e'*e/T2;
    B=B(2:P*K+1,:);
    B=reshape(B',[K,K,P]);
    A=(eye(K)-sum(B,3))^(-1);
    V=A*covE*A';
else
    % If empty then it's is the usual covariance estimator
    V=cov(Y)*((T-1)/T);
end
lagsused=indices{ICpos};

function indices=cov_VAR_dec2bin(maxlag)
% This is a helper function that creates the indices for the search
indices=cell(2^maxlag-1,1);
default=1:maxlag;
for p=1:(2^maxlag-1)
    rem = p;
    selector=false(1,maxlag);
    for i=maxlag-1:-1:0
        if floor(rem/(2^i))
            selector(i+1)=true;
            rem=rem-2^i;
        end
    end
    indices{p}=default(selector);
end
indices{2^maxlag}=[];

author: bashtage

repo: mfe-toolbox


text:
function y = triu(varargin)
% When matrix is empty and there's a non-zero second argument it should give [] (that's what Matlab does), not error
if nargin==2 && isempty(varargin{1})
    y = [];
else
    y = builtin('triu', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
%% Copyright (C) 2014-2017 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym ezplot (@var{f})
%% @defmethodx @@sym ezplot (@var{f1}, @var{f2})
%% @defmethodx @@sym ezplot (@var{f}, @var{dom})
%% @defmethodx @@sym ezplot (@var{f1}, @var{f2}, @var{dom})
%% @defmethodx @@sym ezplot (@dots{}, @var{N})
%% Simple plotting of symbolic expressions.
%%
%% Example parametric plot of a Lissajous Curve:
%% @example
%% @group
%% syms t
%% x = cos(3*t), y = sin(2*t)
%%   @result{} x = (sym) cos(3⋅t)
%%   @result{} y = (sym) sin(2⋅t)
%%
%% ezplot(x, y)                                 % doctest: +SKIP
%% @end group
%% @end example
%%
%% Example plotting the zero level curve of a function of two
%% variables:
%% @example
%% @group
%% syms x y
%% f = x^2 + y^2 - 1;
%% ezplot (f)                                   % doctest: +SKIP
%% @end group
%% @end example
%% Here the curve is defined implicitly by @code{f(x, y) == 0},
%% but we do not enter the @code{== 0} part.
%%
%% See help for the (non-symbolic) @code{ezplot}, which this
%% routine calls after trying to convert sym inputs to
%% anonymous functions.
%%
%% Using sym arguments for @var{dom} and @var{n} can lead to
%% ambiguity where OctSymPy cannot tell if you are specifying @var{n}
%% or @var{f2}.  For example:
%% @example
%% @group
%% syms t
%% f = sin(t);
%% N = sym(50);
%%
%% % parametric plot of f(t), N(t)
%% ezplot(f, N)                                 % doctest: +SKIP
%%
%% % plot f vs t using 50 pts
%% ezplot(f, double(N))                         % doctest: +SKIP
%% @end group
%% @end example
%%
%% The solution, as shown in the example, is to convert the sym to
%% a double.
%%
%% @seealso{ezplot, @@sym/ezplot3, @@sym/ezsurf, @@sym/function_handle}
%% @end defmethod

function varargout = ezplot(varargin)

  % first input is handle, shift
  if (ishandle(varargin{1}))
    fshift = 1;
  else
    fshift = 0;
  end

  firstsym = [];

  for i = (1+fshift):nargin
    if (isa(varargin{i}, 'sym'))
      if ( (i == 1 + fshift) || ...
           (i == 2 + fshift && isscalar(varargin{i})) ...
         )
        % This is one of the fcns to plot, so convert to handle fcn
        % The "i == 2" issscalar cond is to supports ezplot(f, sym([0 1]))

        % Each is function of one var, and its the same var for all
        % (or could be a single function of two variables)
        thissym = symvar(varargin{i});
        assert(length(thissym) <= 2, ...
          'ezplot: plotting curves: functions should have at most two inputs');
        if (isempty(thissym))
          % a number, create a constant function in a dummy variable
          % (0*t works around some Octave oddity on 3.8 and hg Dec 2014)
          thisf = inline(sprintf('%g + 0*t', double(varargin{i})), 't');
          %thisf = @(t) 0*t + double(varargin{i});  % no
        else
          % check variables match (sanity check)
          if (isempty(firstsym))
            firstsym = thissym;
          else
            assert (all (logical (thissym == firstsym)), ...
              'ezplot: all functions must be in terms of the same variables');
          end
          thisf = function_handle(varargin{i});
        end

        varargin{i} = thisf;

      else
        % plot ranges, etc, convert syms to doubles
        varargin{i} = double(varargin{i});
      end
    end
  end

  h = ezplot(varargin{:});

  if (nargout)
    varargout{1} = h;
  end

end

%%!shared hf
%%! % offscreen rendering currently (2016-06) causing crashes:
%%! % e.g., https://savannah.gnu.org/bugs/?44478
%%! hf = figure ('visible', 'off');

%!test
%! % simple
%! syms x
%! f = cos(x);
%! h = ezplot(f);
%! xx = get(h, 'xdata');
%! yy = get(h, 'ydata');
%! assert (abs(yy(end) - cos(xx(end))) <= 2*eps)
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%! % matlab misses endpoint with nodisplay
%! assert (abs(xx(end) - 2*pi) <= 4*eps)
%! assert (abs(yy(end) - cos(2*pi)) <= 4*eps)
%! end

%!test
%! % parametric
%! syms t
%! x = cos(t);
%! y = sin(t);
%! h = ezplot(x, y);
%! xx = get(h, 'xdata');
%! assert (abs(xx(end) - cos(2*pi)) <= 4*eps)

%!error <all functions must be in terms of the same variables>
%! syms x t
%! ezplot(t, x)

%!error
%! syms x t
%! ezplot(t, t*x)

%!test
%! % implicit plot of f(x,y) == 0
%! syms x y
%! f = sqrt(x*x + y*y) - 1;
%! h = ezplot(f);
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%!   xx = get (h, 'xdata');
%!   yy = get (h, 'ydata');
%! else
%!   if (isempty (get (h, 'zdata')))
%!     xx = get (h, 'xdata');
%!     yy = get (h, 'ydata');
%!   else
%!     cm = get (h, 'ContourMatrix');
%!     xx = cm(1, 2:end);
%!     yy = cm(2, 2:end);
%!     assert (cm(1, 1) == 0)
%!     assert (cm(2, 1) == length (xx))
%!   end
%! end
%! assert (abs (max (xx) - 1) <= 0.02)
%! assert (abs (max (yy) - 1) <= 0.02)

%!error
%! % implicit plot supports single function
%! syms x y
%! f = sqrt(x*x + y*y) - 1;
%! g = sqrt(x*x + y*y) - 4;
%! h = ezplot(f, g);

%!test
%! % bounds etc as syms
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%! % this number-of-points option not supported on matlab
%! syms x
%! f = cos(x);
%! h = ezplot(f, [0 2*sym(pi)], sym(42));
%! y = get(h, 'ydata');
%! assert (length(y) == 42)
%! assert (abs(y(end) - cos(4*pi)) <= 4*eps)
%! end

%!test
%! close all

%%!test
%%! close (hf);


text:
function varargout = bst_navigator( varargin )
% BST_NAVIGATOR: Navigation through the Brainstorm database.
%
% USAGE : bst_navigator( 'DbNavigation', action, iDataSets ) : Process only the selected datasets
%         bst_navigator( 'DbNavigation', action )            : Process all the loaded datasets
%         bst_navigator( 'CreateNavigatorMenu', jMenuNavigator)
%
% INPUT : action    : possible values {'NextData', 'NextSubject', 'NextCondition', 'NextResult', 'Previous...'}
%         iDataSets : indices of GlobalData.DataSet array

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2015

eval(macro_method);
end

%% ===== NAVIGATION =====
function DbNavigation( action, iDataSets )
    global GlobalData mutexNavigator timerNavigator;
    % ===== MUTEX =====
    % Use a mutex to prevent the function from being executed more than once at the same time
    if isempty(mutexNavigator) || (mutexNavigator > 1)
        % Entrance accepted
        timerNavigator = tic;
        mutexNavigator = 0;
    else
        % Entrance rejected (another call is not finished,and was call less than 1 seconds ago)
        if ~isempty(timerNavigator)
            mutexNavigator = toc(timerNavigator);
        else
            mutexNavigator = toc;
        end
        disp('Call to bst_navigator() ignored...');
        return
    end

    % ===== PARSE INPUTS =====
    % If no DataSets are specified : process all datasets
    if (nargin < 2)
        iDataSets = 1:length(GlobalData.DataSet);
    end
    % Interpret action
    action = lower(action);
    if ~isempty(strfind(action, 'next'))
        moveDirection = 1;
        action = strrep(action, 'next', '');
    elseif ~isempty(strfind(lower(action), 'previous'))
        moveDirection = -1;
        action = strrep(action, 'previous', '');
    else
        moveDirection = 0;
    end
    % Get protocol information
    ProtocolSubjects = bst_get('ProtocolSubjects');

%% ===== FOR EACH DATA SET =====
    iDSclose = [];
    for iDS = iDataSets
        % ===== CHECK DATASET =====
        % Two necessary conditions : have a SubjectFile defined, and have at least one figure displayed
        if isempty(GlobalData.DataSet(iDS).SubjectFile) || isempty(GlobalData.DataSet(iDS).Figure)
            % Closing all datasets figures
            % close([GlobalData.DataSet(iDS).Figure.hFigure]);
            continue
        end
        % Save modifications
        if GlobalData.DataSet(iDS).Measures.isModified
            panel_record('SaveModifications', iDS);
        end

        % ===== GET DATASET INFORMATION =====
        % === SUBJECT ===
        [sOldSubject, iOldSubject] = bst_get('Subject', GlobalData.DataSet(iDS).SubjectFile);
        % === STUDY ===
        if isempty(GlobalData.DataSet(iDS).StudyFile)
            isNoStudyFile = 1;
        else
            isNoStudyFile = 0;
            [sOldStudy, iOldStudy] = bst_get('Study', GlobalData.DataSet(iDS).StudyFile);
            if isempty(sOldStudy)
                isNoStudyFile = 1;
            end
        end
        % === DATA ===
        if isempty(GlobalData.DataSet(iDS).DataFile)
            isNoDataFile = 1;
        else
            isNoDataFile = 0;
            [tmp__, tmp__, iOldData] = bst_get('DataFile', GlobalData.DataSet(iDS).DataFile);
            if isempty(iOldData)
                isNoDataFile = 1;
            end
        end

        % ===== GET NEW LOCATION =====
        isSubjectFileChanged = 0;
        isStudyFileChanged   = 0;
        isDataFileChanged    = 0;
        isResultsFileChanged = 0;
        waitMsg = '';
        % Switch according to action
        switch lower(action)
            % ===== SUBJECT =====
            case 'subject'
                % === NEW SUBJECT ===
                % Get next/previous subject
                iNewSubject = iOldSubject + moveDirection;
                if (iNewSubject <= 0) || (iNewSubject > length(ProtocolSubjects.Subject))
                    continue
                end
                % Switching to next subject
                sNewSubject = bst_get('Subject', iNewSubject);
                newSubjectName = ProtocolSubjects.Subject(iNewSubject).Name;
                newSubjectFile = ProtocolSubjects.Subject(iNewSubject).FileName;
                isSubjectFileChanged = 1;
                % === NEW STUDY ===
                % If study file is defined
                if ~isNoStudyFile
                    % Get study with NEW subject AND OLD condition
                    newConditionPath = bst_fullfile(newSubjectName, sOldStudy.Condition{:});
                    [sNewStudy, iNewStudy] = bst_get('StudyWithCondition', newConditionPath);
                    % If no study was found
                    if isempty(sNewStudy)
                        continue
                    end
                    isStudyFileChanged = 1;
                else
                    iNewStudy = [];
                    sNewStudy = [];
                end
                % === NEW DATA ===
                % If a DataFile is defined
                if ~isNoDataFile
%                     % If data index is not valid anymore
%                     if (iOldData > length(sNewStudy.Data))
%                         continue
%                     end
%                     iNewData = iOldData;
                    % Cannot switch to an empty folder
                    if isempty(sNewStudy.Data)
                        continue;
                    end
                    % Always the first file of the new condition
                    iNewData = 1;
                    sNewData = sNewStudy.Data(iNewData);
                    % Cannot switch to a continuous raw file
                    if strcmpi(sNewData.DataType, 'raw')
                        continue;
                    end
                    % Display progress bar
                    waitMsg = ['Switching to subject #' num2str(iNewSubject) ': "' newSubjectName '"...' ];
                    % Ask to update the current datafile
                    isDataFileChanged = 1;
                end

            % ===== CONDITION =====
            case 'condition'
                % === NEW SUBJECT ===
                sNewSubject = sOldSubject;
                newSubjectName = sNewSubject.Name;
                newSubjectFile = sNewSubject.FileName;
                % === NEW STUDY ===
                % If study file is defined
                if ~isNoStudyFile
                    % Get all the conditions available for the current subject
                    availableConditions = bst_get('ConditionsForSubject', sOldStudy.BrainStormSubject);
                    % Find current condition in conditions list
                    iOldCond = [];
                    for i=1:length(availableConditions)
                        if isequal(availableConditions{i}, sOldStudy.Condition{1})
                            iOldCond = i;
                            break;
                        end
                    end
                    if isempty(iOldCond)
                        continue;
                    end
                    % Get next/previous condition
                    iNewCond = iOldCond + moveDirection;
                    if (iNewCond <= 0) || (iNewCond > length(availableConditions))
                        continue
                    end
                    % Get study with NEW subject AND OLD condition
                    newConditionPath = bst_fullfile(newSubjectName, availableConditions{iNewCond});
                    [sNewStudy, iNewStudy] = bst_get('StudyWithCondition', newConditionPath);
                    % If no study was found
                    if isempty(sNewStudy)
                        continue
                    end
                    isStudyFileChanged = 1;
                end
                % === NEW DATA ===
                % If a DataFile is defined
                if ~isNoDataFile
%                     % If data index is not valid anymore
%                     if (iOldData > length(sNewStudy.Data))
%                         continue
%                     end
%                     iNewData = iOldData;
                    % Cannot switch to an empty folder
                    if isempty(sNewStudy.Data)
                        continue;
                    end
                    % Always the first file of the new condition
                    iNewData = 1;
                    sNewData = sNewStudy.Data(iNewData);
                    % Cannot switch to a continuous raw file
                    if strcmpi(sNewData.DataType, 'raw')
                        continue;
                    end
                    % Display progress bar
                    waitMsg = sprintf('Switching to condition "%s"...', newConditionPath);
                    % Ask to update the current datafile
                    isDataFileChanged = 1;
                end

            % ===== DATA =====
            case 'data'
                % If DataFile or StudyFile is not defined
                % => Changing Data cannot have consequences : go to next DataSet
                if isNoStudyFile || isNoDataFile
                    continue
                end
                % === NEW STUDY ===
                sNewStudy   = sOldStudy;
                iNewStudy   = iOldStudy;
                % === NEW SUBJECT ===
                sNewSubject = sOldSubject;
                newSubjectFile = sNewSubject.FileName;

                % ===== DATA LISTS =====
                % Get standardized comments
                listComments = cellfun(@str_remove_parenth, {sOldStudy.Data.Comment}, 'UniformOutput', 0);
                % Remove empty matrices
                iEmpty = cellfun(@isempty, listComments);
                if ~isempty(iEmpty)
                    listComments(iEmpty) = {''};
                end
                % Group comments
                [tmp,tmp,iData2List] = unique(listComments);
                % If current file is in a group of trials and not the last one
                iList = iData2List(iOldData);
                % Move to the next one in the trial group
                if (moveDirection == 1) && (nnz(iData2List(iOldData+1:end) == iList) > 0)
                    iNewData = iOldData + find(iData2List(iOldData+1:end) == iList, 1, 'first');
                % Move to the pervious one in the trial group
                elseif (moveDirection == -1) && (nnz(iData2List(1:iOldData-1) == iList) > 0)
                    iNewData = find(iData2List(1:iOldData-1) == iList, 1, 'last');
                % Else go to the next/previous file in the .Data structure
                else
                    iNewData = iOldData + moveDirection;
                end

                % === NEW DATA ===
                % Is action valid ?
                if (iNewData <= 0) || (iNewData > length(sOldStudy.Data))
                    continue
                end
                sNewData = sNewStudy.Data(iNewData);
                % Cannot switch to a continuous raw file
                if strcmpi(sNewData.DataType, 'raw')
                    continue;
                end
                % Display progress bar
                waitMsg = sprintf('Switching to data file #%d: "%s"...', iNewData, sNewData.Comment);
                % Ask to update the current datafile
                isDataFileChanged = 1;

            % ===== RESULT =====
            case 'results'
                % If DataFile or StudyFile is not defined
                % => Changing Results cannot have consequences : go to next DataSet
                if isNoStudyFile || isNoDataFile
                    continue;
                end
                % Current ResultsFiles changed
                isResultsFileChanged = 1;
                % But associated data file did not change
                sNewSubject = sOldSubject;
                newSubjectFile = sNewSubject.FileName;
                sNewStudy = sOldStudy;
                sNewData  = sNewStudy.Data(iOldData);
                % Waiting message
                switch (moveDirection)
                    case 1
                        waitMsg = 'Switching to next results file...';
                    case -1
                        waitMsg = 'Switching to previous results file...';
                end
            otherwise
                bst_error(['Invalid action: "' action '"']);
                return
        end

%% ===== UPDATE GLOBAL STRUCTURES =====
        % === COLLAPSE PREVIOUS NODE ===
        if isSubjectFileChanged
            panel_protocols('CollapseAncestor', 'subject');
        elseif isStudyFileChanged
            panel_protocols('CollapseAncestor', 'study');
        elseif isDataFileChanged
            panel_protocols('CollapseAncestor', 'data');
        end

        % === DATAFILE CHANGED ===
        if isDataFileChanged
            % If other DataSets already exist for the new DataFile, unload them
            iOldDS = bst_memory('GetDataSetData', sNewData.FileName);
            if ~isempty(iOldDS)
                bst_memory('UnloadDataSets', iOldDS);
            end
            % Empty loaded data
            GlobalData.DataSet(iDS).Measures = db_template('Measures');
            % Select new data in tree
            %panel_protocols('SelectNode', [], 'data', iNewStudy, iNewData);
            panel_protocols('SelectNode', [], sNewData.FileName);
        % === STUDY CHANGED ===
        elseif isStudyFileChanged
            % Select new data in tree (for the SUBJECTS and STUDIES exploration modes)
            panel_protocols('SelectStudyNode', iNewStudy);
        % === SUBJECT CHANGED ===
        elseif isSubjectFileChanged
            % Select new data in tree (for the SUBJECTS and STUDIES exploration modes)
            panel_protocols('SelectNode', [], 'subject', iNewSubject, -1);
            panel_protocols('SelectNode', [], 'studysubject', -1, iNewSubject);
        end

        % === UPDATE DATASET ===
        GlobalData.DataSet(iDS).SubjectFile = newSubjectFile;
        % StudyFile
        if ~isNoStudyFile
            GlobalData.DataSet(iDS).StudyFile = sNewStudy.FileName;
        end
        % DataFile
        if ~isNoDataFile
            GlobalData.DataSet(iDS).DataFile = sNewData.FileName;
            % Old time window
            oldTime = GlobalData.UserTimeWindow.Time;
            % Force data reload
            bst_memory('LoadDataFile', sNewData.FileName, 1, 0);
            % If the time was modified
            if any(abs(GlobalData.DataSet(iDS).Measures.Time - GlobalData.UserTimeWindow.Time) > 1e-6)
                % Change the current time window
                GlobalData.UserTimeWindow.Time            = GlobalData.DataSet(iDS).Measures.Time;
                GlobalData.UserTimeWindow.SamplingRate    = GlobalData.DataSet(iDS).Measures.SamplingRate;
                GlobalData.UserTimeWindow.NumberOfSamples = GlobalData.DataSet(iDS).Measures.NumberOfSamples;
                % Change current time if necessary
                if (GlobalData.UserTimeWindow.CurrentTime < GlobalData.UserTimeWindow.Time(1)) || (GlobalData.UserTimeWindow.CurrentTime > GlobalData.UserTimeWindow.Time(2))
                    panel_time('SetCurrentTime', GlobalData.UserTimeWindow.Time(1));
                % Update time panel
                elseif any(abs(oldTime - GlobalData.UserTimeWindow.Time) > 1e-6)
                    panel_time('UpdatePanel');
                end
                % Close all the other datasets
                iDSclose = setdiff(1:length(GlobalData.DataSet), iDS);
            end
        else
            % Reset UserTimeWindow
            GlobalData.UserTimeWindow.Time            = [];
            GlobalData.UserTimeWindow.SamplingRate    = [];
            GlobalData.UserTimeWindow.NumberOfSamples = 0;
        end
        % === CHANNEL FILE ===
        isUpdatedChannelFile = 0;
        if ~isempty(iNewStudy)
            % Get old channel file
            oldChannelFile = GlobalData.DataSet(iDS).ChannelFile;
            % Get new channel
            sNewChannel = bst_get('ChannelForStudy', iNewStudy);
            % If ChannelFile changed
            if ~isempty(sNewChannel) && ~file_compare(sNewChannel.FileName, oldChannelFile)
                newChannelFile = sNewChannel.FileName;
                isUpdatedChannelFile = 1;
                % Update GlobalData structure
                GlobalData.DataSet(iDS).ChannelFile = newChannelFile;
                % Load new channel file
                ChannelMat = in_bst_channel(newChannelFile);
                GlobalData.DataSet(iDS).Channel         = ChannelMat.Channel;
                GlobalData.DataSet(iDS).MegRefCoef      = ChannelMat.MegRefCoef;
                GlobalData.DataSet(iDS).Projector       = ChannelMat.Projector;
                GlobalData.DataSet(iDS).Clusters        = ChannelMat.Clusters;
                GlobalData.DataSet(iDS).IntraElectrodes = ChannelMat.IntraElectrodes;
            end
        end

        % === GET NEW RESULTS ===
        % Get results loaded in DataSet
        if ~isempty(GlobalData.DataSet(iDS).Results)
            oldResultsFiles = {GlobalData.DataSet(iDS).Results.FileName};
        else
            oldResultsFiles = {};
        end
        % New results
        newResultsFiles = cell(1, length(oldResultsFiles));
        if (isDataFileChanged || isResultsFileChanged) && ~isempty(oldResultsFiles)
            % Unload all DataSet results
            bst_memory('UnloadDataSetResult', iDS, 1:length(GlobalData.DataSet(iDS).Results));

            % For each old ResultsFile, build the new ResultsFile (or [] if no equivalent)
            for i=1:length(oldResultsFiles)
                % Get old result file indice in sOldStudy
                [tmp__, tmp__, iOldResult] = bst_get('ResultsFile', oldResultsFiles{i});
                % Get old list of available results in previous sOldStudy
                [tmp__, tmp__, iOldAvailableResults] = bst_get('ResultsForDataFile', sOldStudy.Result(iOldResult).DataFile, iOldStudy);
                % Get indice of old result file in available results
                iOldResultInAvailable = find(iOldAvailableResults == iOldResult);

                % === RESULTS NAVIGATION ===
                if isResultsFileChanged
                    % Data file did not change
                    iNewAvailableResults = iOldAvailableResults;
                    % Get next or previous result
                    iNewResultInAvailable = iOldResultInAvailable + moveDirection;
                    % If requested ResultsFile is NOT available
                    if (iNewResultInAvailable <= 0) || (iNewResultInAvailable > length(iOldAvailableResults))
                        iNewResultInAvailable = [];
                    end
                % === DATA/STUDY NAVIGATION ===
                else
                    % Get all the results available for the new DataFile
                    [tmp__, tmp__, iNewAvailableResults] = bst_get('ResultsForDataFile', sNewData.FileName, iNewStudy);
                    % If equivalent ResultsFile is found for new data
                    if (iOldResultInAvailable <= length(iNewAvailableResults))
                        iNewResultInAvailable = iOldResultInAvailable;
                    else
                        iNewResultInAvailable = [];
                    end
                end

                % If a new results file was found
                if ~isempty(iNewResultInAvailable)
                    iNewResult = iNewAvailableResults(iNewResultInAvailable);
                    newResultsFiles{i} = sNewStudy.Result(iNewResult).FileName;
                else
                    newResultsFiles{i} = [];
                end
            end
        end

%% ===== UPDATE FIGURES =====
        % Process each figure : update or close it
        iFig = 1;
        isDSUnloaded = 0;
        while ~isDSUnloaded && (iDS <= length(GlobalData.DataSet)) && (iFig <= length(GlobalData.DataSet(iDS).Figure))
            Figure = GlobalData.DataSet(iDS).Figure(iFig);
            isFigureRemoved = 0;
            % Update figure appdata
            setappdata(Figure.hFigure, 'SubjectFile', GlobalData.DataSet(iDS).SubjectFile);
            if ~isNoStudyFile
                setappdata(Figure.hFigure, 'StudyFile', GlobalData.DataSet(iDS).StudyFile);
            end
            if ~isNoDataFile
                setappdata(Figure.hFigure, 'DataFile', GlobalData.DataSet(iDS).DataFile);
            end

            % Switch according to figure type
            switch(Figure.Id.Type)
                % === DATA TIME SERIES ===
                case 'DataTimeSeries'
                    if ~isDataFileChanged
                        % Nothing to do
                    elseif ~isempty(Figure.Id.Modality) && (Figure.Id.Modality(1) ~= '$')
                        % Display progress bar
                        bst_progress('start', 'Database navigator', waitMsg);
                        % Update sensors selected in the figure
                        GlobalData.DataSet(iDS).Figure(iFig).SelectedChannels = bst_figures('GetChannelsForFigure', iDS, iFig);
                        % Update TsInfo structure
                        if ~isNoDataFile
                            TsInfo = getappdata(Figure.hFigure, 'TsInfo');
                            if ~isempty(TsInfo) && isfield(TsInfo, 'FileName') && ~isempty(TsInfo.FileName)
                                TsInfo.FileName = GlobalData.DataSet(iDS).DataFile;
                                setappdata(Figure.hFigure, 'TsInfo', TsInfo);
                            end
                        end
                        % Update "isStatic" status
                        setappdata(Figure.hFigure, 'isStatic', (GlobalData.DataSet(iDS).Measures.NumberOfSamples <= 2));
                        % Reload figure
                        bst_figures('ReloadFigures', Figure.hFigure, 0);
                    else
                        isFigureRemoved = 1;
                    end

                % === RESULTS TIME SERIES ===
                case 'ResultsTimeSeries'
                    % Get the results list that is displayed in this figure
                    oldFigRes = getappdata(Figure.hFigure, 'ResultsFiles');
                    % If no or more than one results file : close figure
                    if (length(oldFigRes) ~= 1) || isempty(oldResultsFiles)
                        isFigureRemoved = 1;
                    else
                        % Get new results file for this results file
                        iNewRes = find(file_compare(oldResultsFiles, oldFigRes{1}));
                        newFigRes = newResultsFiles{iNewRes};
                        % If cannot switch
                        if isempty(newFigRes)
                            isFigureRemoved = 1;
                        % Else : update display
                        else
                            % Update figure appdata
                            setappdata(Figure.hFigure, 'ResultsFiles', {newFigRes});
                            setappdata(Figure.hFigure, 'ResultsFile', newFigRes);
                            % Remove figure's unique ID that prevents from writing over this figure
                            GlobalData.DataSet(iDS).Figure(iFig).Id.SubType = '';
                            % Display progress bar
                            bst_progress('start', 'Database navigator', waitMsg);
                            % Call display function again
                            view_scouts({newFigRes}, 'SelectedScouts');
                        end
                    end

                % === RESULTS TIME SERIES ===
                case 'Topography'
                    if isDataFileChanged
                        % Display progress bar
                        bst_progress('start', 'Database navigator', waitMsg);
                        % Call plot function again for new DataFile
                        view_topography(GlobalData.DataSet(iDS).DataFile, Figure.Id.Modality, Figure.Id.SubType);
                    end

                % === 3D FIGURES ===
                case '3DViz'
                    % Get the kind of data represented in this window
                    TessInfo = getappdata(Figure.hFigure, 'Surface');
                    % === PROCESS SURFACES ===
                    for iTess = 1:length(TessInfo)
                        % Get old surface file
                        oldSurfaceFile = TessInfo(iTess).SurfaceFile;
                        % Get new surface file
                        fieldName = ['i' TessInfo(iTess).Name];
                        if (isfield(sNewSubject, fieldName) && ~isempty(sNewSubject.(fieldName)))
                            newSurfaceFile = sNewSubject.Surface(sNewSubject.(fieldName)).FileName;
                        % If old surface type is not available for new subject
                        else
                            % Do not process this surface
                            continue
                        end

                        % === SURFACE CHANGED ===
                        if ~file_compare(oldSurfaceFile, newSurfaceFile)
                            isSurfaceChanged = 1;
                        else
                            isSurfaceChanged = 0;
                        end

                        % === DATA CHANGED ===
                        if ~isempty(TessInfo(iTess).DataSource.Type)
                            % Switch according to type of data on the surface
                            switch (TessInfo(iTess).DataSource.Type)
                                case 'Data'
                                    if isDataFileChanged
                                        % Update Surface description
                                        TessInfo(iTess).DataSource.FileName = sNewStudy.Data(iNewData).FileName;
                                        setappdata(Figure.hFigure, 'Surface', TessInfo);
                                    end
                                case 'Source'
                                    if isResultsFileChanged || isDataFileChanged
                                        % Get new results file for this results file
                                        iNewRes = find(file_compare(oldResultsFiles, TessInfo(iTess).DataSource.FileName));
                                        % If cannot switch
                                        if isempty(iNewRes) || isempty(newResultsFiles{iNewRes(1)})
                                            isFigureRemoved = 1;
                                        % Else : update display
                                        else
                                            % Make sure that only one result is found
                                            iNewRes = iNewRes(1);
                                            newFigRes = newResultsFiles{iNewRes};
                                            % Reset the maximum values for this surface
                                            TessInfo(iTess).DataMinMax = [];
                                            % Update Surface description
                                            TessInfo(iTess).DataSource.FileName = newFigRes;
                                            setappdata(Figure.hFigure, 'Surface', TessInfo);
                                        end
                                    end
                            end
                            if ~isFigureRemoved
                                % Display progress bar
                                bst_progress('start', 'Database navigator', waitMsg);
                                % View new surface / new data on surface
                                view_surface_data(TessInfo(iTess).SurfaceFile, TessInfo(iTess).DataSource.FileName, Figure.Id.Modality);
                            end

                        % === SURFACE CHANGED / NO DATA ===
                        elseif isSurfaceChanged
                            % Display progress bar
                            bst_progress('start', 'Database navigator', waitMsg);
                            % Display new surface
                            view_surface(newSurfaceFile);
                            % Remove previous surface
                            panel_surface('RemoveSurface', Figure.hFigure, iTess);
                        end
                    end
                    % === PROCESS SENSORS ===
                    if isUpdatedChannelFile
                        % Get the elements to be displayed
                        isMarkers = ~isempty(Figure.Handles.hSensorMarkers);
                        isLabels  = ~isempty(Figure.Handles.hSensorLabels);
                        % Update channels display
                        %view_channels(GlobalData.DataSet(iDS).ChannelFile, Figure.Id.Modality, isMarkers, isLabels);
                        figure_3d('ViewSensors', Figure.hFigure, isMarkers, isLabels);
                    end

                % === FREQ DISPLAYS ===
                case {'Timfreq', 'Spectrum', 'Pac', 'Image'}
                    isFigureRemoved = 1;

                % === CONNECTIVITY ===
                case 'Connect'
                    isFigureRemoved = 1;

                % Other types of windows: just close them
                otherwise
                    isFigureRemoved = 1;
            end

            % If close request for this figure
            if isFigureRemoved
                % If it is the last DataSet figure => DataSet will be unloaded
                isDSUnloaded = (length(GlobalData.DataSet(iDS).Figure) <= 1);
                % Close figure
                close(Figure.hFigure);
                % Do not increment the figure => A figure was deleted, so the next figure has the same indice
            else
                % Go to next figure
                iFig = iFig + 1;
            end
        end
    end
    % Unload some datasets
    if ~isempty(iDSclose)
    	bst_memory('UnloadDataSets', iDSclose);
    end
    % Update events lists
    panel_record('UpdatePanel');
    % Hide progress bar
    bst_progress('stop');
    % Release mutex
    mutexNavigator = [];
end

%% ====== CREATE NAVIGATOR MENU ======
function CreateNavigatorMenu(jMenuNavigator) %#ok<DEFNU>
    import java.awt.event.KeyEvent;
    import javax.swing.KeyStroke;
    import org.brainstorm.icon.*;
    % Create items
    jItemPrevSubj = gui_component('MenuItem', jMenuNavigator, [], 'Previous subject',   IconLoader.ICON_PREVIOUS_SUBJECT,   [], @(h,ev)DbNavigation('PreviousSubject'));
    jItemNextSubj = gui_component('MenuItem', jMenuNavigator, [], 'Next subject',       IconLoader.ICON_NEXT_SUBJECT,       [], @(h,ev)DbNavigation('NextSubject'));
    jMenuNavigator.addSeparator();
    jItemPrevCond = gui_component('MenuItem', jMenuNavigator, [], 'Previous condition', IconLoader.ICON_PREVIOUS_CONDITION, [], @(h,ev)DbNavigation('PreviousCondition'));
    jItemNextCond = gui_component('MenuItem', jMenuNavigator, [], 'Next condition',     IconLoader.ICON_NEXT_CONDITION,     [], @(h,ev)DbNavigation('NextCondition'));
    jMenuNavigator.addSeparator();
    jItemPrevData = gui_component('MenuItem', jMenuNavigator, [], 'Previous data file', IconLoader.ICON_PREVIOUS_DATA,      [], @(h,ev)DbNavigation('PreviousData'));
    jItemNextData = gui_component('MenuItem', jMenuNavigator, [], 'Next data file',     IconLoader.ICON_NEXT_DATA,          [], @(h,ev)DbNavigation('NextData'));
    % Set keyboards shortcuts
    jItemNextSubj.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F1), 0)); % F1
    jItemNextCond.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F2), 0)); % F2
    jItemNextData.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F3), 0)); % F3
    jItemPrevSubj.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F1), KeyEvent.SHIFT_MASK)); % SHIFT+F1
    jItemPrevCond.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F2), KeyEvent.SHIFT_MASK)); % SHIFT+F2
    jItemPrevData.setAccelerator(KeyStroke.getKeyStroke(int32(KeyEvent.VK_F3), KeyEvent.SHIFT_MASK)); % SHIFT+F3
end

author: brainstorm-tools


text:
function fiff_write_short(fid, kind, data)
%
% fiff_write_short(fid, kind, data)
%
% Writes a 16-bit integer tag to a fif file
%
%     fid           An open fif file descriptor
%     kind          Tag kind
%     data          The integers to use as data
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.1  2006/04/27 22:38:37  msh
%   Splitting an empty list now results in an empty output.
%   Added fiff_write_double and fiff_write_short
%   Write an array of floats, ints, and doubles instead of just one value.
%   Fixed help text of fiff_pick_channels.
%
%

me = 'MNE:fiff_write_short';

if nargin ~= 3
    error(me, 'Incorrect number of arguments');
end

FIFFT_SHORT = 2;
FIFFV_NEXT_SEQ = 0;
nel = numel(data);
datasize = nel * 2;
count = fwrite(fid, int32(kind), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(FIFFT_SHORT), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(datasize), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(FIFFV_NEXT_SEQ), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int16(data), 'int16');
if count ~= nel
    error(me, 'write failed');
end

return


text:
% function cl = enlarge(cl)
%
% enlarge cluster

function cl = enlarge_cluster(cl)

    XYZ = cl.XYZ;
    XYZnew = XYZ;

    for i = 1:3
        Xtmp = XYZ;
        Xtmp(i, :) = Xtmp(i, :) + 1;
        XYZnew = [XYZnew Xtmp];         % add to dim
        Xtmp(i, :) = Xtmp(i, :) - 2;
        XYZnew = [XYZnew Xtmp];      % subtract 1 from dim
    end

    XYZnew = unique(round(XYZnew)', 'rows')';
    cl.XYZ = XYZnew;
    cl.XYZmm = voxel2mm(XYZnew, cl.M);
    cl.Z = ones(1, size(cl.XYZ, 2));

    cl.numVox = size(cl.XYZmm, 2);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Trajectory = Create_Trajectory(Params)
%
% FUNCTION:
%   This function is used to draw the road that the tractor trailer truck
%   is driving on. A different representation is used for the forward and
%   reverse versions of the road.
%
%   This function defines a trajectory as a series of points and gains
%   associated with that point. Trajectory is a struct with two fields.
%
% INPUT:
%   Params = a struct of parameters for this problem. The field Traj is
%   used in this function. See Set_Parameters.m and Trajectory_Examples.m
%   for more details.
%
% OUTPUT:
%   Trajectory.States = a (2xN) cell array with the state at each of the N points
%   Trajectory.Inputs = a (2xN) cell array the nominal actuator effort
%   Trajectory.Gains = A (2xN) cell array.
%   {1,:} ==> Forward
%   {2,:} ==> Reverse

Npts = Params.Traj.Npts;

if Params.Traj.Func == false

    %Set constraint points:
    SetPts = Params.Traj.SetPts;
    Order = Params.Traj.Order;

    P = polyfit(SetPts(:,1),SetPts(:,2),Order);
    x = linspace(min(SetPts(:,1)),max(SetPts(:,1)),Npts+1);
    y = polyval(P,x);
else
   x = Params.Traj.x;
   y = Params.Traj.y;
end

%Now we need the desired trailer angle theta, which is defined as the angle
%between the trailer center line and the positive y axis.

%Vector direction between successive points
x_step = diff(x);
y_step = diff(y);
angle = atan2(y_step,x_step);  %angle from the positive x axis
th = angle-pi/2;

%Assume that the cab angle should always try to be centered
phi = zeros(size(th));

%Truncate the vectors
x=x(1:end-1);
y=y(1:end-1);

% % % %PLOT THINGS
% % %
% % % subplot(2,1,2)
% % % plot(x(1:end-1),th)
% % % title('Angle vs X Position')
% % %
% % % subplot(2,1,1)
% % % plot(x(1:end-1),y(1:end-1))
% % % axis equal
% % % title('Road')

%% Now do gain matrix calculations

Trajectory.States = cell(2,Npts);
Trajectory.Inputs = cell(2,Npts);
Trajectory.Gains = cell(2,Npts);
for i=1:Npts
   %Drive Forward
        State = [x(i); y(i); th(i); phi(i)];
        Input = [1;0];  %[speed;steering]
        Trajectory.States{1,i} = State;
        Trajectory.Inputs{1,i} = Input;
        Trajectory.Gains{1,i} =  Set_Gain_Matrix(State, Input, Params);

   %Drive Backward
        State_Rev = [x(i); y(i); th(i)+pi; phi(i)];
        Input_Rev = [-1,0];   %[speed;steering]
        Trajectory.States{2,i} = State_Rev;
        Trajectory.Inputs{2,i} = Input_Rev;
        Trajectory.Gains{2,i} =  Set_Gain_Matrix(State_Rev, Input_Rev, Params);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function om = init_indexed_name(om, set_type, name, dim_list)
%INIT_INDEXED_NAME  Initializes the dimensions for an indexed named set.
%
%   OM.INIT_INDEXED_NAME(SET_TYPE, NAME, DIM_LIST)
%
%   Initializes the dimensions for an indexed named variable, constraint
%   or cost set.
%
%   Variables, constraints and costs are referenced in OPT_MODEL in terms
%   of named sets. The specific type of named set being referenced is
%   given by SET_TYPE, with the following valid options:
%       SET_TYPE = 'var'   => variable set
%       SET_TYPE = 'lin'   => linear constraint set
%       SET_TYPE = 'nle'   => nonlinear equality constraint set
%       SET_TYPE = 'nli'   => nonlinear inequality constraint set
%       SET_TYPE = 'qdc'   => quadratic cost set
%       SET_TYPE = 'nlc'   => nonlinear cost set
%       SET_TYPE = 'cost'  => legacy cost set
%
%   Indexed Named Sets
%
%   A variable, constraint or cost set can be identified by a single NAME,
%   such as 'Pmismatch', or by a name that is indexed by one or more indices,
%   such as 'Pmismatch(3,4)'. For an indexed named set, before adding the
%   indexed variable, constraint or cost sets themselves, the dimensions of
%   the indexed set must be set by calling INIT_INDEXED_NAME, where
%   DIM_LIST is a cell array of the dimensions.
%
%   Examples:
%       %% linear constraints with indexed named set 'R(i,j)'
%       om.init_indexed_name('lin', 'R', {2, 3});
%       for i = 1:2
%         for j = 1:3
%           om.add_lin_constraint('R', {i, j}, A{i,j}, ...);
%         end
%       end
%
%   See also OPT_MODEL, ADD_VAR, ADD_LIN_CONSTRAINT, ADD_NLN_CONSTRAINT,
%            ADD_QUAD_COST, ADD_NLN_COST and ADD_LEGACY_COST.

%   MATPOWER
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% call parent method (also checks for valid type for named set)
om = init_indexed_name@opt_model(om, set_type, name, dim_list);

%% add type-specific info about this named set
switch set_type
    case 'cost'         %% cost set
        %% use column vector if single dimension
        if length(dim_list) == 1
            dim_list = {dim_list{:}, 1};
        end
        empty_cell  = cell(dim_list{:});
        om.cost.data.N.(name)   = empty_cell;
        om.cost.data.Cw.(name)  = empty_cell;
        om.cost.data.vs.(name)  = empty_cell;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function r8_sign_test ( )

%*****************************************************************************80
%
%% R8_SIGN_TEST tests R8_SIGN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2014
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;

  r8_test = [ -1.25, -0.25, 0.0, +0.5, +9.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_SIGN_TEST\n' );
  fprintf ( 1, '  R8_SIGN returns the sign of an R8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    R8  R8_SIGN(R8)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    r8 = r8_test(test);
    s = r8_sign ( r8 );
    fprintf ( 1, '  %8f  %8f\n', r8, s );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ui = mrViewAddMap(ui, newMR);
%
% Add a new map to a ui struct.
% ui = mrViewAddMap(ui,newMR);
%
% ras 09/2005
if notDefined('ui'),	ui = mrViewGet;				end
if ishandle(ui),		ui = get(ui, 'UserData');	end

% allow map to be a struct array of new maps, and recursively add them:
if length(newMR) > 1
	for i = 1:length(newMR)
		ui = mrViewAddMap(ui, newMR(i));
	end
	return
end

% integer check
if ~isa(newMR.data,'double'), newMR.data = double(newMR.data); end

if isempty(ui.maps) % first map loaded
	ui.maps = newMR;
else                % add to other maps
	% check if name exists; if so, make unique
	exNames = {ui.maps.name}';
	test = strmatch(newMR.name,exNames);
	if ~isempty(test)
		newMR.name = sprintf('%s (%i)',newMR.name,length(test)+1);
	end

	ui.maps(end+1).baseXform = []; % will initialize new entry

	% copy fields over, in a memory-efficient way (?)
	for f = fieldnames(newMR)'
		ui.maps(end).(f{1}) = newMR.(f{1});
		newMR.(f{1}) = []; % clear to save memory
	end
end
ui.maps(end).baseXform = mrBaseXform(ui.mr, ui.maps(end));

% if loading a map for the first time, open
% a new overlay for the map:
if length(ui.maps)==1
	ui = mrViewAddOverlay(ui);
end

% also add the map to any popup lists of maps which
% may exist for overlays:
for o = 1:length(ui.overlays)
	if isfield(ui.overlays(o),'mapPopup') & ...
			ishandle(ui.overlays(o).mapPopup)
		set(ui.overlays(o).mapPopup,'String',{ui.maps.name});

		for j = 1:length(ui.overlays(o).threshMap)
			str = get(ui.overlays(o).threshMap(j),'String');
			str{end+1} = ui.maps(end).name;
			set(ui.overlays(o).threshMap(j),'String',str);
		end
	end
end

% update the info panel to include the new map as an option:
if checkfields(ui,'panels','info') & ishandle(ui.panels.info)
	ui = mrViewSet(ui,'infopanel');
end

ui = mrViewRefresh(ui);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = functionPoints(disc)
%FUNCTIONPOINTS   Points at which functions are discretized.

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

pointsFun = @(n) trigtech.trigpts(n);
[varargout{1:nargout}] = valsDiscretization.points(disc, pointsFun);

end


text:
% Plot the aircraft's geometry setup.
% Example call: plotAircraftGeometry('ExperimentalCarrier', 'State_1.5alpha_12.5ms.mat', 'F:\svn\dev\matlab\tornado\T135_export');

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

function plotAircraftGeometry(aircraft_geometry_filename, state, tornado_root_directory)

currDir = pwd;

cd(tornado_root_directory);
settings=config('startup');

% Load the aircraft
cd(settings.acdir)
load(aircraft_geometry_filename); % aircraft geometry is now in variable 'geo'.

cd(tornado_root_directory)

% Generate the Lattice
latticetype=0; % Tornado lattice type
[lattice,ref]=fLattice_setup2(geo,state,latticetype);

%  Plot the geometry
geometryplot(lattice,geo,ref);
%bodyplot(body);

% back to the original directory
cd(currDir);

author: jrgenerative

repo: fixed-wing-sim


text:
function [] = plot_route(ax, route, parsed_osm)
% plot (over map) the route found by route planner
%
% usage
%   PLOT_ROUTE(ax, route, parsed_osm)
%
% input
%   ax = axes object handle.
%   route = matrix of nodes traversed along route, as returned by the
%           route_planner function.
%   parsed_osm = parsed OpenStreetMap XML data, as returned by the
%                parse_openstreetmap function.
%
% 2012.04.24 (c) Ioannis Filippidis, jfilippidis@gmail.com
%
% See also ROUTE_PLANNER, PARSE_OPENSTREETMAP.

% empty path ?
if isempty(route)
    warning('path is empty. This means that the planner found no path.')
    return
end

nodexy = parsed_osm.node.xy;
start_xy = nodexy(:, route(1, 1) );
path_xy = nodexy(:, route);
path_end = nodexy(:, route(1, end) );

held = takehold(ax);

plotmd(ax, start_xy, 'Color', 'm', 'Marker', 'o', 'MarkerSize', 15)
plotmd(ax, path_xy, 'Color', 'r', 'LineStyle', '--', 'LineWidth', 5)
plotmd(ax, path_end, 'Color', 'c', 'Marker', 's', 'MarkerSize', 15)

restorehold(ax, held)

author: johnyf

repo: openstreetmap


text:
function nd = lagrange_interp_nd_size ( m, n_1d )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_ND_SIZE sizes an M-dimensional Lagrange interpolant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N_1D(M), the order of the 1D rule to be used
%    in each dimension.
%
%    Output, integer ND, the number of points in the product grid.
%
  nd = prod ( n_1d(1:m) );

  return
end


text:
function [P, V, A] = extremePools(model, positivity, inequality)
% Calculates the extreme pools of a stoichiometric model using the vertex / facet enumeration package
%
% INPUT:
%    model:                structure with:
%
%                            * .S - `m x n` Stoichiometric matrix with integer coefficients. If no
%                              other inputs are specified it is assumed that all reactions are
%                              reversible and `S.v = 0`
%                            * .SIntRxnBool - `n x 1` boolean vector with 1 for internal reactions
%                            * .description - description
%
% OPTIONAL INPUT:
%    positivity:           {0, (1)} if `positivity == 1`, then positive orthant base
%    inequality:           {(0), 1} if `inequality == 1`, then use two inequalities rather than a single equaltiy
%
% .. Author: - lrs by David Avis, McGill University

[nMet, nRxn] = size(model.S);

if isfield(model, 'SIntRxnBool')
    A = model.S(:, model.SIntRxnBool)';
else
    A = model.S';
end

if nnz(A - round(A))
    figure
    spy(A - round(A))
    title('S-round(S)')
    error('Stoichiometric coefficients must be all integers')
end

a = zeros(size(A, 1),1);

if isfield(model, 'description')
    filename = model.description;
else
    filename = 'model';
end

if ~exist('positivity', 'var')
    positivity = 1;
end
if ~exist('inequality', 'var')
    inequality = 0;
end

suffix = '';
if positivity
    suffix = [suffix 'pos_'];
else
    suffix = [suffix 'neg_'];
end
if inequality
    suffix = [suffix 'ineq'];
else
    suffix = [suffix 'eq'];
end

% no inequalities
D = [];
d = [];

% no linear objective
f = [];

% no shell script
sh = 0;

% INPUT
% A          matrix of linear equalities A*x=(a)
% D          matrix of linear inequalities D*x>=(d)
% filename   base name of output file
%
% OPTIONAL INPUT
% positivity {0,(1)} if positivity==1, then positive orthant base
% inequality {0,(1)} if inequality==1, then use two inequalities rather than a single equaltiy
% a          boundry values for matrix of linear equalities A*x=a
% d          boundry values for matrix of linear inequalities D*x>=d
% f          linear objective for a linear optimization problem in rational arithmetic
%            minimise     f'*x
%            subject to   A*x=(a)
%                         D*x>=(d)
lrsInputHalfspace(A, D, filename, positivity, inequality, a, d, f, sh);

% pause(eps)
[status, result] = system('which lrs');
if ~isempty(result)
    % call lrs and wait until extreme pathways have been calculated
    systemCallText = ['lrs ' pwd filesep filename '_' suffix '.ine > ' pwd filesep filename '_' suffix '.ext'];
    [status, result] = system(systemCallText);
else
    error('lrs not installed or not in path')
end

[P, V] = lrsOutputReadRay([filename '_' suffix '.ext']);
P = P';
A = A';
if any(any(P * A ~= 0))
    warning('extreme pool not in nullspace of stoichiometric matrix')
end

author: opencobra

repo: cobratoolbox


text:
% example: get s
%
% load ROIfile
% for each voxel:
%   - remove session means
%   - HP filter
%   - convert to percent change
%   - fit model
%   - save fits and residual

/Users/tor/Documents/Tor_Documents/CurrentExperiments/VNL/ROIFILES
load ../vnl_condf.mat

hrf = hrf ./ max(hrf);
[LinearX] = getPredictors(vnl_condf,hrf);
LinearX = LinearX - repmat(mean(LinearX),size(LinearX,1),1);

xtxitx = pinv(LinearX);

D = dir('*vis.mat');
b_by_ymean = [];

for s = 1:length(D)

    % load the file
    load(D(s).name)
    eval(['ROI = ' D(s).name(1:end-4)])

    clear b, clear sigma

for i = 1:size(ROI.ts.indiv,2),

        if length(ROI.adjustedy) == 6020
			[y,ymean] = voistat('adjusty',ROI.ts.indiv(:,i),41,.5,7);
		else
			[y,ymean] = voistat('adjusty',ROI.ts.indiv(:,i),41,.5,8);
		end

    b(:,i) = xtxitx(:,1:length(y)) * y;
    sigma(i) = std(y - LinearX(1:length(y),:) * b(:,i));,

end

%mean(b')
%sigma

%btmp = xtxitx * ROI.adjustedy
%figure; plot(ROI.adjustedy); hold on; plot(LinearX * btmp,'r')

    B(s,:) = mean(b');
    SIGMA(s) = mean(sigma);
    SNR(s,:) = mean(b ./ repmat(sigma,6,1),2)';

end

save vnl_snr_out B SIGMA SNR b_by_ymean D

B(1:2,:) = [];
sigma(1:2) = [];
SNR(1:2,:) = [];

% mean B, 1 2 5 6 10 11 flashes
% 0.8258    1.1419    1.2200    1.2461    0.8393    0.9499

% mean SNR
% 0.4724    0.6557    0.7084    0.7278    0.5038    0.5574

% mean sigma
% 1.6684    percent change

% to get betas from UNSCALED data, to assess fits of betas by mean
% level...does the signal scale with baseline levels?

b_by_ymean = [];

for s = 3:length(D)

    % load the file
    load(D(s).name)
    eval(['ROI = ' D(s).name(1:end-4)])

    clear b, clear sigma

for i = 1:size(ROI.ts.indiv,2),

    %   ALTER VOISTAT TO GET RID OF SCALING BEFORE DOING THIS
        if length(ROI.adjustedy) == 6020
			[y,ymean] = voistat('adjusty',ROI.ts.indiv(:,i),41,.5,7);
		else
			[y,ymean] = voistat('adjusty',ROI.ts.indiv(:,i),41,.5,8);
		end

    b(:,i) = xtxitx(:,1:length(y)) * y;
    sigma(i) = std(y - LinearX(1:length(y),:) * b(:,i));,

    b_by_ymean(end+1,:) = [b(:,i)' ymean s];

end

end

figure;gscatter(b_by_ymean(:,7),b_by_ymean(:,1),b_by_ymean(:,8),'rgbycmk','o+xv^s');
t1 = unique(b_by_ymean(:,8));
for i = t1'

    tmp = b_by_ymean(b_by_ymean(:,8) == i,[1 7]);
    t2 = corrcoef(tmp);  b_by_mean_r(i) = t2(1,2);
    bb = polyfit(tmp(:,2),tmp(:,1),1); h = refcurve(bb);
    set(h,'Color',rand(1,3))
end

set(gca,'FontSize',16); xlabel('Mean scanner signal'),ylabel('Estimated response magnitude'),title('Estimated magnitude as a function of mean signal')


text:
function [x] = denormdata(xn,xmean,xstd)
%DENORMDATA  De-normalize normalized data
%
%  Description
%    X = DENORMDATA(XN,XMEAN,XSTD) de-normalize XN using
%    precomputed XMEAN and XSTD.
%
%  See also NORMDATA
%
% Copyright (c) 2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  if nargin<3
    error('Too few arguments')
  end
  x=bsxfun(@plus,bsxfun(@times,xn,xstd),xmean);

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [a,f,B,W]=ldatrace(b,w,n,c)
%LDATRACE Calculates an LDA transform to maximize trace discriminant [a,f,B,W]=(b,w,n,c)
% If a feature vector X can come from one of several class and W and B are respectively
% the within-class and between-class covariance matrices, then the generalized Fisher discriminant
% F=trace(W\B) is a measure of how well the feature vector discriminates between the classes.
% If we choose a rectangular (tall, skinny) transformation matrix, we can define a smaller
% feature vector Y=A'*X. The aim of this routine is to choose A to maximize the Fisher
% discriminant. We assume that W is positive definite and B is positive semi-definite.
% The input argument C allows the uset to pre-specify some of the columns of A.
%
% Inputs:
%     w[m,m] = within class covariance matrix of x
%     b[m,m] = between class covariance matrix of x [default = I]
%     n is the number of columns in output matrix A [default = M]
%     c[m,r] specifies the first few columns of A to be predefined values [default = null)
%
% Outputs:
%     a[m,n] is the transformation matrix: y=a'*x
%     f[n,1] gives the incremental gain in f value for successive columns of A
%     B(n,n) gives the between-class covariance matrix of y
%     W[n,n] gives the within-class covariance matrix of y

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

m=size(b,1);    % dimension of data vectors
if nargin<4
    r=0;
    if nargin<3
        n=m;
        if nargin<2
            w=eye(m);
        end
    end
else
    r=size(c,2);    % number of columns that are pre-specified
end
if r
    if n>r          % need to find additional vectors
        g=chol(w);
        v=g\null(c'*g');
        [p,l,q]=svd(v'*b*v);
        a(:,r+1:n)=v*p(:,1:n-r);
        a(:,1:r)=c;
    else
        a=c;        % no new vectors to find
    end
    if nargout>1
        ari=a/triu(qr(chol(a'*w*a))); % matrix a must be of full rank
        f=diag(ari'*b*ari);
    end
else
    [g,d]=eig(b,w,'qz');
    [ds,is]=sort(-diag(d));
    a=g(:,is(1:n));
    if nargout>1
        f=-ds(1:n);
    end
end
if nargout > 2
    B=a'*b*a;
    if nargout > 3
        W=a'*w*a;
    end

end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

temp = sigmoid(X * theta);
p = (temp >= 0.5)

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

end

author: sfvsfv

repo: Mathematical-modeling


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function displayStat(confs, objs, method, windowHandle, cpath)

load(deblank(char(objs)));

[p,fname,e] = fileparts(deblank(char(objs)));

switch lower(method)
    case 'res_t_map'
        switch deblank(char(confs.Overlay))
            case 'p-value'
                cutoff = confs.Threshold_p_value;
                signal = pvalue;
            case 't-map'
                df = length(grInfo)-1;
                cutoff = tinv(1-(confs.Threshold_p_value/2),df);
                signal = tstats;
        end
        patch_overlay(atlas_vertices, faces, signal, cutoff, deblank(char(confs.Overlay)), vtnorm, deblank(char(confs.Colormap)),windowHandle,fname);
    case 'res_pca'
        project_PCA(fvecs, eigenvecs, eigenvals, deblank(char(confs.Mesh)), confs.MaxSPHARMDegree, confs.Level, confs.Sigma, windowHandle,fname, cpath);
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function y = inv_pos( x )

%INV_POS   Internal cvx version.

y = pow_cvx( x, -1, 'pow_p' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
%---------------------------------------------------------
% Low-pass FIR digital differentiator (LPFIRDD) design   -
% via constrained quadratic programming (QP)             -
% [Copt,c]=lpfirdd(N,alpha,beta,r,idraw)                 -
% By Dr Yangquan Chen		019-07-1999                   -
% Email=<yqchen@ieee.org>; URL=http://www.crosswinds.net/~yqchen/
% --------------------------------------------------------
% LPFIRDD: only 1st order derivative estimate
% total taps=2N. c(i)=-c(i+N+1); c(N+1)=0 (central point)
%
%          -N      -N+1            -1                    N
% FIR=c(1)z   +c(2)z    +...+ c(N)z  + 0 + ... + c(2N+1)z
%
%       N
%     ------
%     \                    j   -j
%      >       Copt(j) * (z - z  )
%     /
%     ------
%      j=1
% N: Taps  (N=2, similar to sgfilter(2,2,1,1)
% alpha ~ beta: transit band of frequency
%				    (in percentage of Nyquest freq)
% r: the polynomial order. r<=N Normally, set it to 1.
%---------------------------------------------------------------
function [Copt,bd]=lpfirdd(N,alpha,beta,r,idraw)
% testing parameters
% alpha=1./pi;beta=1.5/pi;N=10;r=1;idraw=1;
if (alpha>beta)
   disp('Error in alpha (alpha<=beta)');return;
end
if ((beta>1) | (beta <0))
   disp('Error in Beta! (beta in [0,1]');return;
end
if ((alpha>1) | (alpha <0))
   disp('Error in Alpha! (Alpha in [0,1]');return;
end
% default r=1
if (r<1); r=1; end
% matrix W
W=zeros(r,N);
for ix=1:N;for jx=1:r;
      W(jx,ix)=ix^(2*jx-1);
end;end
%matrix L
L=zeros(N,1);
if (beta>alpha)
	for ix=1:N
	   L(ix)=(alpha*sin(ix*beta*pi)-beta*sin(ix*alpha*pi))/ix/ix/(beta-alpha);
	end
elseif (beta==alpha)
	for ix=1:N
	   L(ix)=(ix*alpha*pi*cos(ix*alpha*pi)-sin(ix*alpha*pi))/ix/ix;
	end
end
% matrix e
ex=zeros(r,1);ex(1)=1;
% optimal solution
Copt=W'*inv(W*W')*(ex + 2.*W*L/pi)-2.*L/pi
Copt=Copt/2;
% fr plots
if (idraw==1)
   bd=[-fliplr(Copt'),0,Copt']';
	%ad=1;sys_sg=tf(bd',ad,1./Fs);bode(sys_sg)
	Fs=12790;nL=N;nR=N;npts=1000;%w=logspace(0,4,npts);
	w=((1:npts)-1)*pi/npts;
	j=sqrt(-1);ejw=zeros(nL+nR+1,npts);
	for ix=(-nL:nR)
		ejw(ix+nL+1,:)=exp(j*ix*w);
	end
	freq=bd'*ejw;
	figure;subplot(2,1,1)
	plot(w/pi*Fs/2,(abs(freq)));grid on;
	hold on; ax=axis;ax(2)=Fs/2;axis(ax);
	xlabel('freq. (Hz)');ylabel('amplitude (dB)');
	subplot(2,1,2)
	plot(w/pi*Fs/2,180*(angle(freq))/pi );grid on;
	hold on; ax=axis;ax(2)=Fs/2;axis(ax);
   xlabel('freq. (Hz)');ylabel('phase anlge (deg.)');

	figure;subplot(2,1,1);Fs=12600; % Hz for U8
	semilogx(w*Fs/pi/2,20*log10(abs(freq)));grid on;
	hold on; ax=axis;ax(2)=Fs/2;axis(ax);
	semilogx([Fs/2,Fs/2],[ax(3),ax(4)],'o-r');grid on;
	xlabel('freq. (Hz)');ylabel('amplitude (dB)');
	subplot(2,1,2)
	semilogx(w*Fs/pi/2,180*(angle(freq))/pi );grid on;
	hold on; ax=axis;ax(2)=Fs/2;axis(ax);
	semilogx([Fs/2,Fs/2],[ax(3),ax(4)],'o-r');grid on;
	xlabel('freq. (Hz)');ylabel('phase anlge (deg.)');
end
return

author: Sable

repo: mcbench-benchmarks


text:
function [res,fopt] = nlp_ape(Ccell,rrPar,R,t,theta)
blk         = rrPar.blk;
tBound      = rrPar.translationBound;
dBound      = rrPar.depthBound;
FOV         = rrPar.FOV;
C           = Ccell{1};
N           = length(theta);
elements.A  = specialeuclideanfactory(3,1);
elements.B  = obliquefactory(1,N);
manifold    = productmanifold(elements);

problem.M   = manifold;

warning('off', 'manopt:getHessian:approx')
% Define the problem cost function and its Euclidean gradient.
problem.cost  = @(x) ape_cost(x,C);
problem.egrad = @(x) ape_egrad(x,C);

% Numerically check gradient consistency (optional).
% checkgradient(problem);
% Solve.
x0.A.R               = R;
x0.A.t               = t;
x0.B                 = theta';
options.verbosity    = 0;
options.tolgradnorm  = 1e-6;
options.maxiter      = 1000;
[xopt, fopt, output, options] = trustregions(problem,x0,options);
% [xopt, fopt, output, options] = arc(problem,x0,options);
Ropt     = xopt.A.R;
topt     = xopt.A.t;
thetaopt = xopt.B;
constraintviolationR          = norm(Ropt*Ropt'-eye(3),'fro');
constraintviolationtheta      = max(abs(thetaopt.^2 - 1));
firstorderopt                 = output(end).gradnorm;
fprintf('        MANOPT: itr: %3d, constraint violation: %3.2e, %3.2e, gradnorm: %3.2e, cost: %3.8e.\n',...
    length(output),constraintviolationR,constraintviolationtheta,firstorderopt,fopt);
if max([constraintviolationR,constraintviolationtheta]) < 1e-8 ...
        && firstorderopt < 1e-3
    % Do nothing
else
    fopt    = inf;
end
if ~check_translation(topt,tBound,FOV)
    fopt    = inf;
end
% Rnew = project2SO3(Ropt);
% rnew = Rnew(:);
% tnew = topt;
% thetanew = sign(thetaopt);
vnew     = lift_ape_v1(Ropt(:),topt,thetaopt,tBound,dBound,FOV);
Xmat     = {vnew{1} * vnew{1}';vnew{2} * vnew{2}';vnew{3} * vnew{3}';vnew{4} * vnew{4}'};
% fopt     = blktrace(blk,Xmat,Ccell);

res.X     = Xmat;
res.R     = Ropt;
res.t     = topt;
res.theta = thetaopt;
end

function f = ape_cost(x,C)
R       = x.A.R;
r       = R(:);
t       = x.A.t;
theta   = x.B;
theta   = theta(:);
v       = lift_pcr_v3(r,t,theta);
f       = v'*C*v;
end

function g = ape_egrad(x,C)
R       = x.A.R;
r       = R(:);
t       = x.A.t;
theta   = x.B;
theta   = theta(:);
N       = length(theta);
v       = lift_pcr_v3(r,t,theta);

gv      = 2 * v' * C; % 1 x n

vdr     = [sparse(1,9);...
            speye(9);...
            sparse(3,9);...
            sparse(N,9);...
            kron(theta,[speye(9);sparse(3,9)])]; % n x 9
gr      = (gv * vdr)'; % 9 x 1

vdt     = [ sparse(1,3);...
            sparse(9,3);...
            speye(3);...
            sparse(N,3);...
            kron(theta,[sparse(9,3);speye(3)])]; % n x 4
gt      = (gv * vdt)'; % 3 x 1

vdtheta = [sparse(13,N);...
            speye(N);...
            kron(speye(N),[r;t])]; % n x N
gtheta  = (gv * vdtheta); % 1 x N

g.A.R   = reshape(gr,3,3);
g.A.t   = gt;
g.B     = gtheta;
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
filename='Bridge_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
  function arg = remove_spaces(arg)
%|function arg = remove_spaces(arg)
%| replace extra spaces at ends of matrix string array with zeros.
%| also add extra '0' at end to be sure to null terminate string
%| Copyright May 2000, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

if streq(arg, 'test'), remove_spaces_test, clear arg, return, end

for ii=1:size(arg,1)
	jj = size(arg,2);
	while (arg(ii,jj) == ' ')
		arg(ii,jj) = 0;
		jj = jj - 1;
		if jj < 1
			warn('bug in %s', mfilename)
			keyboard
		end
	end
end
arg(:,end+1) = 0;

arg(arg(:,1) == 10, :) = []; % remove blank lines too

function remove_spaces_test
str = strvcat('a', 'bb ', 'c c ');
out = remove_spaces(str);
tmp = char([97 0 0 0 0; 98 98 0 0 0; 99 32 99 0 0]);
jf_equal(out, tmp)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt