text:
function Zf = filt2(Z,res,lambda,filtertype)
% filt2 performs a highpass, lowpass, bandpass, or bandstop 2D gaussian filter on gridded data such as
% topographic, atmospheric, oceanographic, or any kind of geospatial data. This function is designed to
% make it easy to remove features longer or shorter than a given characteristic wavelength. The
% input grid can contain NaNs!
%
%% Syntax
%
%  Zf = filt2(Z,res,lambda,filtertype)
%
%% Description
%
% Zf = filt2(Z,res,lambda,filtertype) filters 2D dataset Z that has resolution res,
% to an approximate wavelength lambda.  If the filtertype is 'lp' or 'hp' for lowpass
% or highpass, lambda must be a scalar value.  If the filtertype is 'bp' or 'bs' for
% bandpass or bandstop, lambda must be a two-element array of the two cutoff wavelengths.
%
%% Explanation of this type of filter
% For an explanation of this filter with pictures, type
%
%   showdemo filt2_documentation
%
% into your command window.
%
%% Example
% Consider a 100 km by 100 km elevation dataset that has a resolution of 200 m.  It has
% some long 25 km wavelength features aligned with north/south direction, some short
% ~5 km features oriented diagonally, and a considerable amount of random noise.  There's
% also a block of missing data. Here's what your datset looks like:
%
% res = 0.2; % 200 m resolution
% x = 0:res:100; % eastings from 0 to 100 km
% y = 0:res:100; % northings from 0 to 100 km
% [X,Y] = meshgrid(x,y);
%
% % Z contains 25 km features, ~5 km diagonal features, and noise:
% Z = cos(2*pi*X/25)+cos(2*pi*(X+Y)/7)+randn(size(X));
%
% % Z also has some missing data:
% Z(100:120,100:120) = nan;
%
% subplot(1,3,1)
% imagesc(x,y,Z);
% axis xy image
% caxis([-1 1])
% title ' original Z matrix '
% xlabel ' eastings (km) '
% ylabel ' northings (km) '
%
% % Get the lowpass-filtered version of Z:
% Zlow = filt2(Z,res,15,'lp');
%
% subplot(1,3,2)
% imagesc(x,y,Zlow);
% axis xy image
% caxis([-1 1])
% title ' 15 km lowpass filtered Z matrix '
% xlabel ' eastings (km) '
% ylabel ' northings (km) '
%
% % Get the highpass-filtered version of Z:
% Zhi = filt2(Z,res,15,'hp');
%
% subplot(1,3,3)
% imagesc(x,y,Zhi);
% axis xy image
% caxis([-1 1])
% title ' 15 km highpass filtered Z matrix '
% xlabel ' eastings (km) '
% ylabel ' northings (km) '
%
%% Author Info
% This function was written by Chad A. Greene of the University of Texas Institute for
% Geophysics in November 2016; however, all I did was repackage Carlos Adrian Vargas Aguilera's
% superb ndnanfilter function, which can be found here: http://www.mathworks.com/matlabcentral/fileexchange/20417.
% Many thanks to Carlos for his well-thought-out code and clear documentation.
%
% See also conv2, imgaussfilt, and imfilter.

%% Input checks

narginchk(4,4)
% assert(license('test','image_toolbox')==1,'Error: I''m sorry, the filt2 function requires the Image Processing Toolbox.')
assert(ismatrix(Z)==1,'Input error: Z must be a 2d matrix.')
%assert(isscalar(res)==1,'Input error: res must be a scalar value.')
assert(ismember(lower(filtertype),{'lp','hp','bp','bs'}),'Input error: filtertype must be ''hp'', ''lp'', ''bp'', or ''bs''.')
if lambda<=(2*max(res))
   warning('Nyquist says the wavelength should exceed two times the resolution of the dataset, which is an unmet condition based on these inputs. I''ll give you some numbers, but I would''t trust ''em if I were you.')
end

if ismember(lower(filtertype),{'bp','bs'})
   assert(numel(lambda)==2,'Input error: Wavelength lambda must be a two-element array for a bandpass filter.')
else
   assert(isscalar(lambda)==1,'Input error: Wavelength lambda must be a scalar for lowpass or bandpass filters.')
end

%% Design filter:

% 2*pi*sigma is the wavelength at which the amplitude is multiplied by a factor of about 0.6 (more exactly, exp(-0.5))
sigma = (lambda(1)./res) /(2*pi);
%f = fspecial('gaussian',2*ceil(2.6*sigma)+1,sigma);
% WJP implementation without image processing toolbox, subfunc at bottom of
% this function
f = gaussian2D([2*ceil(2.6*sigma(1))+1 2*ceil(2.6*sigma(end))+1], sigma);

%% Now filter the data

switch lower(filtertype)
   %WJP incorrect usage originally but happened to work if imfilter is available
   case 'lp'
      Zf = ndnanfilter(Z,f,size(f),[],[],{'replicate'}); % ndnanfilter is Carlos Adrian Vargas Aguilera's excellent function, which is included as a subfunction below.

   case 'hp'
      Zf = Z - ndnanfilter(Z,f,size(f),[],[],{'replicate'});

   case 'bp'
      Zf =  filt2(filt2(Z,res,max(lambda),'hp'),res,min(lambda),'lp');

   case 'bs'
      Zf = filt2(Z,res,max(lambda),'lp') - filt2(Z,res,min(lambda),'hp');

   otherwise
      error('No such filter type.')
end

end

function [Y,W] = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT,WNAN)
% NDNANFILTER   N-dimensional zero-phase digital filter, ignoring NaNs.
%
%   Syntax:
%         Y = ndnanfilter(X,HWIN,F);
%         Y = ndnanfilter(X,HWIN,F,DIM);
%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT);
%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT);
%         Y = ndnanfilter(X,HWIN,F,DIM,WINOPT,PADOPT,WNAN);
%     [Y,W] = ndnanfilter(...);
%
%   Input:
%     X      - Data to be filtered with/without NaNs.
%     HWIN   - Window function handle (or name) or numeric multidimensional
%              window to be used (without NaNs). See WINDOW for details.
%              Default:   @rectwin  or 'rectwin' (moving average).
%     F      - A vector specifying the semi-width of the window for each
%              dimension. The final window's width will be 2*F+1.
%              Default: 3 (i.e. a 1-dimensional window of width 6).
%     DIM    - If F is a single scalar, the window will be applied through
%              this dimension; otherwise, this will be ignored.
%              Default: columns (or the first non-singleton dimension).
%     WINOPT - Cell array specifying optional arguments for the window
%              function HWIN (in addition to the width).
%              Default: {} (window's defaults).
%     PADOPT - Cell array specifying the optional arguments for the
%              PADARRAY MATLAB's function (in addition to the array X and
%              the padsize: 2*F+1). If the function is not found, data is
%              padded with zeros or the specified value: try {mean(X(:))}
%              for example.
%              Default: {'replicate'} (repeats border elements of X).
%              Default: {0} (pads with zeros if PADARRAY not found).
%     WNAN   - Integer indicating NaNs treatment and program behaviour!:
%              0: Filters data and interpolates NaNs         (default).
%              1: Filters data but do not interpolates NaNs
%              2: "Do not filters data" but interpolates NaNs!
%              See the NOTEs below
%
%   Output:
%     Y      - Filtered X data (same size as X!).
%     W      - N-dimensional window with central symmetry generated by a
%              special subfunction called NDWIND. See the description below
%              for details.
%
%   Description:
%     This function applies a N-dimensional convolution of X with W, using
%     the MATLAB's IMFILTER or CONVN function. One important aspect of the
%     function is the generation of the N-dimensional window (W) from the
%     specified function and width, which cannot be done with MATLAB's
%     functions. Besides, unlike MATLAB's FILTER, FILTER2 and IMFILTER,
%     NaNs elements are taken into account (ignored).
%
%     The N-dimensional window is generated from rotating the 1-dimensional
%     output of the HWIN function, through each of the N-dimensions, and
%     then shrinking it through each of its axes in order to fit the
%     specified semi-widths (F). This is done in the included subfunction
%     named NDWIND. In this way, the window has central symmetry and do not
%     produce a phase shift on X data.
%
%     By default, the edges are padded with the values of X at the borders
%     with the PADARRAY MATLAB's function. In this way, the edges are
%     treated smoothly. When PADARRAY is not found, the program performs
%     zero-padding.
%
%   Notes:
%     * The use of semi-widths F's is to force the generated window to be
%       even and, therefore, the change of phase is null.
%     * The window function HWIN should output an even function, otherwise,
%       it won't generate an error but the user should be aware that this
%       program will consider only the last half of it.
%     * The function window should return a monotonically decreasing
%       result, this restriction is because I try to avoid the use of FZERO
%       function, for example, to find the expanding/shrinking factors.
%     * If the user has an already generated window, it can be used in HWIN
%       instead of a function handle or name.
%     * Accepts empty value for any input. When X is empty, the program can
%       be used as a N-dimensional window generator.
%     * NaNs elements surrounded by no-NaNs elements (which will depend on
%       window width) are the ones that will be interpolated. The others
%       are leaved untouched.
%     * When WNAN=2, the programs acts like an NAN-interpolat/GAP-filling,
%       leaving untouched the no-NaNs elements but the filtering is
%       perfomed anyway. I recomend the default behaviour (WNAN=0) in order
%       to keep the filtered data in the workspace, and then use the code
%       at the end of this function to get/remove the interpolated NaNs
%       (see the example).
%     * The program looks for the IMFILTER and PADARRAY functions from the
%       Image Processing Toolbox. If not found, then CONVN is used instead
%       (slower) and pads with zeros or the given value. In this latter
%       case, if border elements are NaNs, the window won't work properly.
%
%   Example:
%     FWIN = 'hamming';
%     F = [13 8];
%     N = 100;
%     Pnoise = 0.30;
%     PNaNs  = 0.20;
%     X = peaks(N);                                     % original
%     Y = X + ((rand(size(X))-0.5)*2)*max(X(:))*Pnoise; % add noise
%     Y(round(1 + (N^2-1).*rand(N^2*PNaNs,1))) = NaN;   % add NaNs
%     [Z0,W] = ndnanfilter(Y,FWIN,F);                   % filters
%     Z1 = Z0; Z2 = Y; inan = isnan(Y);
%     Z1(inan) = NaN;
%     Z2(inan) = Z0(inan);
%     subplot(231), imagesc(X), clim = caxis; axis equal tight
%                   title('Original data')
%     subplot(232), imagesc(Y),  caxis(clim), axis equal tight
%                   title('Data + NOISE + NaNs')
%     subplot(234), imagesc(Z0), caxis(clim), axis equal tight
%                   title('FILTERS + NaNs interpolation')
%     subplot(235), imagesc(Z1), caxis(clim), axis equal tight
%                   title('FILTERS ignoring NaNs')
%     subplot(236), imagesc(Z2), caxis(clim), axis equal tight
%                   title('GAP-filling with interpolated NaNs')
%     subplot(233), imagesc(-F(1):F(1),-F(2):F(2),W), axis equal tight,
%                    title([upper(FWIN) ' 2D window']), view(2)
%
%   See also: FILTER, FILTER2 and CONVN; WINDOW from the Signal Processing
%   Toolbox; and FWIND1, FWIND2, FSPECIAL, IMFILTER and PADARRAY from the
%   Image Processing Toolbox.

%   Copyright 2008 Carlos Adrian Vargas Aguilera
%   $Revision: 1.2 $  $Date: 2008/06/30 18:00:00 $

%   Written by
%   M.S. Carlos Adrian Vargas Aguilera
%   Physical Oceanography PhD candidate
%   CICESE
%   Mexico, 2008
%   nubeobscura@hotmail.com
%
%   Download from:
%   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
%   tType=author&objectId=1093874

%   1.0    Release (2008/06/23 10:30:00)
%   1.1    Fixed Bug adding an extra dimension of unitary width.
%   1.2    Fixed Bug with ynan.

% Use the IMFILTER function? (faster than CONVN):
yimfilter = (exist('imfilter','file')==2);

% Use the PADARRAY function (or zero padding):
ypadarray = (exist('padarray','file')==2);

% Check inputs and sets defaults of principal arguments:
if nargin<3 || nargin>7
 error('Filtern:IncorrectNumberOfInputs',...
  'At least three inputs are needed and less than 7.')
end
if isempty(HWIN)
 HWIN = 'rectwin';
end
if isempty(F)
 F = 3;
end
N = length(F);
S = size(X);
% Secondary arguments:
if N && (nargin<4 || isempty(DIM))
 DIM = find(S~=1,1);   %  DIM = min(find(S~=1));
 if isempty(DIM), DIM = 1; end
end
if nargin<5 || isempty(WINOPT)
 WINOPT = {};
end
if nargin<6 || isempty(PADOPT)
 if ypadarray
  PADOPT = {'replicate'};
 else
  PADOPT = {0};
 end
elseif ~ypadarray && ~isnumeric(PADOPT{1})
 PADOPT = {0};
end
if nargin<7 || isempty(WNAN)
 WNAN = 0;
end

% Selects the 1-dimensional filter or set a row vector:
if N==1
 a = zeros(1,DIM);
 a(DIM) = F;
 F = a;
 clear a
end

% Checks if the window input is a function or an array:
if ~isa(HWIN,'function_handle') && ~ischar(HWIN)
 W = HWIN;
else
 W = [];
end

% If no input data but two outputs then generates the window only:
if isempty(X)
 Y = [];
 if nargout==2 && ~isempty(W)
  W = ndwind(HWIN,F,WINOPT{:});
 end
 return
end

% Generates the window:
if isempty(W)
 W = ndwind(HWIN,F,WINOPT{:});
end

% Check for NaN's:
inan = isnan(X);
ynan = any(inan(:));                       % Bug fixed 30/jun/2008
if ynan
 X(inan) = 0;
else
 factor = sum(W(:));
end

% Filtering:
if yimfilter                                % Use IMFILTER (faster)
 if ~isa(X,'double')
  X = double(X);
 end
 if ~isa(W,'double')
  W = double(W);
 end
 if ynan
  Y = imfilter(X,W       ,PADOPT{:},'conv');
 else
  Y = imfilter(X,W/double(factor),PADOPT{:},'conv');
 end
else                                        % Use CONVN
 % Sets F and S of equal sizes.
 F = reshape(F,1,N);
 Nx = numel(S);
 if N<Nx
  F(N+1:Nx) = 0;
 elseif N>Nx
  S(Nx+1:N) = 1;
 end
 F2 = 2*F;
 % Pads the borders:
 if ypadarray
  ind    = padarray(false(S),F2,true     );    % Index of the padding.
  Y      = padarray(X       ,F2,PADOPT{:},'both');
 elseif length(PADOPT{1})==1
  ind2 = cell(N,1);
  for n = 1:N
   ind2{n} = F2(n) + (1:S(n)).';
  end
  ind          = true(2*F2+S);
  Y            = repmat(PADOPT{1},2*F2+S);
  ind(ind2{:}) = false;
  Y(ind2{:})   = X;
 else % No padding at all
  Y    = X;
  ind  = false(S);
  warning('Ndnanfilter:PaddingOption','Do not perfom any padding.')
 end
 % Convolutes both arrays:
 if ynan
  Y = convn(Y,W       ,'same');
 else
  Y = convn(Y,W/factor,'same');
 end
 %  Eliminates the padding:
 Y(ind) = [];
 Y      = reshape(Y,S);
end

% Estimates the averages when NaNs are present:
if ynan
 if yimfilter
  factor       = imfilter(double(~inan),W,PADOPT{:},'conv');
 else
  if ypadarray
   factor      = padarray(~inan,F2,PADOPT{:});
  elseif length(PADOPT{1})==1 % (won't work properly with NaNs at borders)
   factor          = ind;
   factor(ind2{:}) = ~inan;
  else
   factor = ~inan;
  end
  factor      = convn(factor,W,'same');
  factor(ind) = [];
  factor      = reshape(factor,S);
 end
 Y = Y./factor;
end

% What about NaNs?:
if     WNAN == 1       % Leave NaNs elements untouched!
 Y(inan) = NaN;
elseif WNAN == 2       % Leave no-NaNs elements untouched!!!
 X(inan) = Y(inan);
 Y = X;
end
end

function W = ndwind(HWIN,F,varargin)
% NDWIND Generate a N-Dimensional zero-phase window.
%
%   Syntax:
%     W = ndwind(HWIN,F);
%     W = ndwind(HWIN,F,OPT);
%
%   Input:
%     HWIN - Window function handle. See WINDOW for details. By default
%            uses: @rectwin (a rectangular window).
%     F    - A vector specifying the semiwidth of the window for each
%            dimension. The window's width will be 2*F+1. By default uses:
%            3 (i.e. a window of width 6).
%     OPT  - Cell array specifying optional arguments for the window
%            function. By default uses: {[]} (window's defaults).
%
%   Output:
%     W    - N-Dimensional window with central symmetry.
%
%   Description:
%     In the axes of each dimension, W has a 1-D window defined as
%              feval(HWIN,2*F(n)+1), n = 1,...,N.
%     That is, they are defined by the same window function but have
%     different widths. So, this program creates another widther window (at
%     least 201 points), with the same definition, and finds how much the
%     former windows should be expanded in order to fit the latter one.
%
%     Afterwards, the coordinates of every point are expanded accordingly
%     and the value of the window in those points are found by linear
%     interpolation with the bigger window.
%
%     In resume, it is like rotating this big window through every
%     dimension and then shrinking it through each of its axes to fix the
%     specified widths.
%
%   Notes:
%     * Because of the use of the semi-widths F's, all the generated
%       windows are even. Therefore the change of phase is null.
%     * The window function HWIN should output an even function, otherwise,
%       it won't generate an error but this program will consider only the
%       last half of it.
%     * The window should be monotonically decreasing.
%     * Instead of the handle window, it can be given as a string:
%       'hamming' instead of @hamming, for example.
%     * Uses the MATLAB's function FUNC2STR.
%
%   Example:
%     W = ndwind(@hamming,[3 2])
%     % Results:
%     W =
%
%              0         0    0.0800         0         0
%              0    0.1417    0.3100    0.1417         0
%              0    0.3966    0.7700    0.3966         0
%         0.0800    0.5400    1.0000    0.5400    0.0800
%              0    0.3966    0.7700    0.3966         0
%              0    0.1417    0.3100    0.1417         0
%              0         0    0.0800         0         0
%
%
%   See also: WINDOW from the Signal Processing Toolbox; and FWIND1,
%   FWIND2, and FSPECIAL from the Image Processing Toolbox.

%   Copyright 2008 Carlos Adrian Vargas Aguilera
%   $Revision: 1.1 $  $Date: 2008/06/26 19:30:00 $

%   Written by
%   M.S. Carlos Adrian Vargas Aguilera
%   Physical Oceanography PhD candidate
%   CICESE
%   Mexico, 2008
%   nubeobscura@hotmail.com
%
%   Download from:
%   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
%   tType=author&objectId=1093874

%   1.0    Release (2008/06/23 10:30:00)
%   1.1    Fixed Bug adding an extra dimension of unitary width.

% Check inputs:
if nargin<1 || isempty(HWIN)
 HWIN = 'rectwin';
end
if nargin<2 || isempty(F)
 F = 3;
end

% Rectangular wind?:
if isa(HWIN,'function_handle')
 HWIN = func2str(HWIN);
end
if strcmpi(HWIN,'rectwin')
 W = ones([2*F(:).'+1 1]);
 return
end

% Generate the BIG window (only the last half):
FBIG         = max([100; F(:)]);
BIGw         = feval(HWIN,2*FBIG+1,varargin{:});
BIGw(1:FBIG) = [];       % Deletes the first half.
rBIGw        = 0:FBIG;   % Window argument (distance).

% Axial windows widths:
N  = numel(F);
F  = reshape(F,1,N);
F  = [F 0];             % BUG fixed by adding an extra dimension.
N  = N+1;
F2 = 2*F+1;

% Pre-allocates the final window and the expanded axis:
W  = zeros(F2);
An = cell(N,1);
Ae = An;

% Generates the index and expanded axes:
for n = 1:N

 % Generate temporally the window in the n-axis:
 wn = feval(HWIN,F2(n),varargin{:});

 % Finds the expansion factors (Note: the window should tends to zero):
 if F(n)
  piv = wn(end);
  ind = (BIGw == piv);
  if ~any(ind)
   ind1 = (BIGw >= piv); ind1 = length(ind1(ind1));
   ind2 = (BIGw <= piv); ind2 = length(ind2(~ind2))+1;
   if ind2>FBIG+1
    r = rBIGw(ind1);
   else
    r = interp1(BIGw([ind1 ind2]), rBIGw([ind1 ind2]),piv);
   end
  else
   r = rBIGw(ind);
  end
  Ef = r/F(n);
 else
  Ef = 1;
 end

 % Reversed index and expanded n-axis (for the following grid):
 An{n} = (F(n):-1:0);
 Ae{n} = An{n}*Ef;

end

% Estimates the expanded distances outside the axes (only at the 1st
% quarter):
% Note: In a 2-Dimensional matrix, by the 1st quarter of a matrix I mean
% the first 1/4 piece of the matrix after you divided it throuh the middle
% row and column. In N-dimensions it would be the 1st 1/2^N part.
gride4      = cell(N,1);
[gride4{:}] = ndgrid(Ae{:});
R4          = sqrt(sum(reshape([gride4{:}],prod(F+1),N).^2,2));

% Generates the window and linear index in the 1st quarter:
grid4     = cell(N,1);
[grid4{:}]= ndgrid(An{:});
in        = (R4<=rBIGw(end));           % Looks for elements inside window.
W4        = zeros(F+1);                 % 1st quarter of the window.
W4(in)    = interp1(rBIGw,BIGw,R4(in)); % Interpolates the window values.
for n=1:N                               % Linear index on the 1st quarter.
 grid4{n} = flip(grid4{n}+1,n);
end
ind4      = sub2ind(F2,grid4{:});

% Index of permutations to fill the N-D window:
np = 2^N-1;
ip = zeros(1,np);
for n = 1:N
 ini  = 2^(n-1);
 step = ini*2;
 ip(ini:step:np) = n;
end

% Fills the N-D window by flipping W4 and the index:
ones4       = false(F2);    % Avoids using new FALSE function
ones4(ind4) = true;
W(ones4)    = W4;
for kp = ip
 W4         = flip(W4,kp);
 ones4      = flip(ones4,kp);
 W(ones4)   = W4;
end
end

function h = gaussian2D(siz, std)

    % create the grid of (x,y) values
    siz = (siz-1)./2;
    [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));

    std = max(std);
    sr = siz(2)/siz(1);
    if sr > 1
        sy = sr^2; sx = 1;
    else
        sx = sr^2; sy = 1;
    end

    % analytic function
    h = exp(-(sx*x.*x + sy*y.*y)/(2*std*std));

    % truncate very small values to zero
    h(h<eps*max(h(:))) = 0;

    % normalize filter to unit L1 energy
    sumh = sum(h(:));
    if sumh ~= 0
        h = h/sumh;
    end
end

author: CHLNDDEV


text:
%kfliptoreal 'flips phase of complex object into a real one '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros fliptoreal.pane file
%
% Parameters:
% Toggle: ac 'Make Real (Add conjugate)', default: 0: 'adds the conjugate image to make it real'
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kfliptoreal(i, {'ac',0;'i','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% fliptoreal - flips phase of complex object into a real one
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kfliptoreal(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,..] = kfliptoreal(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'ac', 0;'i', '__input';'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[1,0,0];
was_set=istoggle * 0;
paramtype={'Toggle','InputFile','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 'fliptoreal"  -k'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Test Difference Routines

% Comparing Different Strategies

% 1) DERIVEST Suite
% 2) MKL djacobi
% 3) adiff

%% Test 1
clc
fun = @(x) 2*x(1) + 3*x(2);
x0 = [1;1];

tic
% g1 = gradest(fun,x0)
toc

tic
[g2,stat] = mklJac(fun,x0)
toc

tic
g3 = autoJac(fun,x0)
toc

%%
clc
f = symJac(fun)

%% Test 2
clc
fun = @(x) 2.1*x(1) + 3*x(2)^2 + 4.2*sqrt(x(3));
x0 = [1;1;2];

tic
% g1 = gradest(fun,x0)
toc

tic
[g2,stat] = mklJac(fun,x0)
toc

tic
g3 = autoJac(fun,x0)
toc

%% Test 3
clc
nlcon = @(x) [8 - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 - x(1) + x(2) - x(3) + x(4);
              10 - x(1)^2 - 2*x(2)^2 - x(3)^2 - 2*x(4)^2 + x(1) + x(4);
              5 - 2*x(1)^2 - x(2)^2 - x(3)^2 - 2*x(1) + x(2) + x(4)];
x0 = [1;1;1;1];

tic
% g1 = jacobianest(nlcon,x0)
toc

tic
[jac,stat] = mklJac(nlcon,x0,3)
toc

tic
g3 = autoJac(nlcon,x0)
toc

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
clear all; close all; clc;

t=linspace(0,10,1000);
x=sin(t)+3*cos(t/3)+4*sin(t/4);

plot(t,x,'b');
grid on;
xlabel('time(s)'); ylabel('x(t)');
hold on;

Ts=0.1;
num_of_space=10/Ts;

t_DT=linspace(0,9,num_of_space);
x_DT=sin(t_DT)+3*cos(t_DT/3)+4*sin(t_DT/4);
% alpha=1;
alpha=t_DT(2)-t_DT(1);
    for ii=1:1:length(t_DT)
      if x_DT(ii)>=0
          rectangle('Position',[t_DT(ii),0,alpha,x_DT(ii)*(1/alpha)],'EdgeColor','r');
      else
          rectangle('Position',[t_DT(ii),1*x_DT(ii)*(1/alpha),alpha,abs(x_DT(ii)*(1/alpha))],'EdgeColor','r');
      end

    end

% axis([-1 11 -0.5 2.5])

% stem(t_DT,x_DT)

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
% OCEANSITES_DEMO Plot all time series variables that have a certain
% 'standard_name' in a specified bounding box and time period

% WHOI GI-CAT catalog server, which harvests data from multiple
% THREDDS Servers:
open_url='http://geoport.whoi.edu/gi-cat/services/opensearch';

% can search for any text, not just variable 'standard_names':
%var='relative_humidity'; %free text search string
var='sea_water_temperature';
%var='sea_water_salinity';

% time_var='time'; %time coordinate variable
bbox =[-180 180 -70 80];  % [lon_min lon_max lat_min lat_max]

% specify a certain search range.  Like this...
start=[1990 1 1 0 0 0]; % specified start
stop =[2000 1 1 0 0 0]; % specified stop

% or like this....
start=now_utc-28;  % last 28 days
stop=now_utc;          % now

%%%%%%%%%%%%%%%%%%%% end of user defined input %%%%%%%%%%%%%%%%
% opensearch query
q.endpoint=open_url;
q.bbox=bbox;
q.time_start=datestr(start);% convert to ISO
q.time_end=datestr(stop);% convert to ISO
q.string_text=var

disp(['Querying ' q.endpoint ' via OpenSearch']);
[links,params]=opensearch(q);
dap=links2dap(links); % find only the OPeNDAP links
char(dap)
disp('Accessing data via OPeNDAP');
for i=1:length(dap);
    figure(i);
    nc=ncgeodataset(dap{i});
    vars=nc.variables;
    for j=1:length(vars); %loop through variables to find standard_names
        vart=nc.geovariable(vars{j});
        std_name=vart.attribute('standard_name');
        if strcmp(std_name,var),
            jd= vart.timewindowij(start, stop);
            t=vart.data(jd.index);  %extact these indices from dataset
            plot(jd.time,t);datetick;grid;...
            ylabel(sprintf('%s [%s]',var,vart.attribute('units')),...
                'interpreter','none');...
            title(nc.attribute('title'));...
            drawnow
        end
    end
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function jac = p40_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P40_JAC evaluates the jacobian for problem p40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  e = p40_param ( 'GET', 'E', [] );

  jac(1,1) = ( 2.0 * y(1) - t ) / e;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mv = mv_removeOutliers(mv,criterion,thresh,plotFlag);
%
% mv = mv_removeOutliers([mv,criterion,thresh,plotFlag]);
%
% Remove outlier voxels, whose extreme values lie
% beyond a certain threshold. This threshold
% can either be a certain # of standard deviations
% away from the mean, or greater/less than a certain
% absolute value.
%
% criterion is 0 (count standard deviations  criterion) or 1
% (use abs value criterion). thresh is the threshold # of S.D.s or
% threshold value to use.. Omitting either brings up a dialog.
%
% after removing voxels, the mv struct is re-initialized.
%
% ras 05/05.
if ieNotDefined('mv')
    mv = get(gcf,'UserData');
end

plotFlag = 0;

if ieNotDefined('thresh') | ieNotDefined('criterion')
	ui(1).string = 'Threshold value:';
	ui(1).fieldName = 'thresh';
	ui(1).style = 'edit';
	ui(1).value = '1';

	ui(2).string = 'Criterion:';
	ui(2).fieldName = 'criterion';
	ui(2).list = {'Standard Deviations from Mean' 'Absolute % Signal'};
	ui(2).style = 'popup';
	ui(2).value = 1;

    ui(3).string = 'Report on response distributions';
    ui(3).fieldName = 'plotFlag';
    ui(3).style = 'checkbox';
    ui(3).value = 0;

	resp = generalDialog(ui,'Remove Outliers');
    thresh = str2num(resp.thresh);
    criterion = cellfind(ui(2).list,resp.criterion)-1;
    plotFlag = resp.plotFlag;
end

amps = mv_amps(mv);

% find outliers, based on specified criterion
switch criterion
    case 0, % remove voxels w/ mean > thresh std. deviations from
        amps = permute(nanmean(amps),[3 2 1]);
        sigma = nanstd(abs(amps(:)));
        mu = nanmean(abs(amps));
        muTotal = nanmean(abs(amps(:))); % overall mean across voxels
        outliers = find(mu>muTotal+thresh*sigma);
    case 1, % tSeries extends greater than thresh pct. signal
        [rows cols] = find(abs(mv.tSeries) > thresh);
        outliers = unique(cols);
end

% first, check that not ALL voxels are outliers...
if length(outliers)==size(mv.roi.coords,2)
    msg = ['Warning: Given the selected criteria, ALL voxels '...
           'were found to be outliers! No action taken.'];
    mrMessage(msg);
    return
end

% report the # and location of outliers found
fprintf('%i Outlier Voxels Found. \n',length(outliers));
% fprintf('Locations: \n ');
% for i = 1:length(outliers)
%     fprintf('\t %s \n', num2str(mv.roi.coords(:,i)));
% end
% fprintf('\n')

% keep only voxels that aren't outliers
keep = setdiff(1:size(mv.tSeries,2),outliers);
mv.tSeries = mv.tSeries(:,keep);

% also set roi coords to reflect only kept voxels
mv.coords = mv.coords(:,keep);
mv.roi.coords = mv.roi.coords(:,keep);

% plot distributions if selected
if plotFlag==1
    h = figure;
    subplot(211);
    if criterion > 0,
        amps = permute(nanmean(mv.voxAmps),[3 2 1]);
    end
    hist(amps(:),100);
    title('Before')
end

% recompute voxData matrix
mv.voxData = er_voxDataMatrix(mv.tSeries,mv.trials,mv.params);

% recompute voxAmps matrix
mv.voxAmps = er_voxAmpsMatrix(mv.voxData,mv.params);

% plot distributions if selected
if plotFlag==1
    figure(h);
    subplot(212);
    newAmps = permute(nanmean(mv.voxAmps),[3 2 1]);
    hist(newAmps(:),100);
    title('After')
    xlabel('Amplitude, % Signal')
    ylabel('# Voxels')

    figure(mv.ui.fig);
end

% if a UI exists, set as user data
if isfield(mv.ui,'fig') & ishandle(mv.ui.fig)
    set(mv.ui.fig,'UserData',mv);
    multiVoxelUI; % refresh UI
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function S = matsignt(T)
%MATSIGNT    Matrix sign function of a triangular matrix.
%            S = MATSIGN(T) computes the matrix sign function S of the
%            upper triangular matrix T using a recurrence.

%            Called by SIGNM.

if ~isequal(T,triu(T)), error('Matrix must be upper triangular.'), end

n = length(T);

S = diag( sign( diag(real(T)) ) );
for p = 1:n-1
   for i = 1:n-p

      j = i+p;
      d = T(j,j) - T(i,i);

      if S(i,i) ~= -S(j,j)  % Solve via S^2 = I if we can.

         % Get S(i,j) from S^2 = I.
         k = i+1:j-1;
         S(i,j) = -S(i,k)*S(k,j) / (S(i,i)+S(j,j));

      else

         % Get S(i,j) from S*T = T*S.
         s = T(i,j)*(S(j,j)-S(i,i));
         if p > 1
            k = i+1:j-1;
            s = s + T(i,k)*S(k,j) - S(i,k)*T(k,j);
         end
         S(i,j) = s/d;

      end

   end
end

author: opencobra

repo: cobratoolbox


text:
function make3plots(X1, Y1, Y2, Y3)
%CREATEFIGURE1(X1,Y1,Y2,Y3)
%  X1:  vector of x data
%  Y1:  vector of y data
%  Y2:  vector of y data
%  Y3:  vector of y data

%  Auto-generated by MATLAB on 14-Feb-2007 22:38:36

% Create figure
figure1 = figure;

% Create subplot
subplot1 = subplot(3,1,1,'Parent',figure1,'YGrid','on','XGrid','on');
box('on');
hold('all');

% Create plot
plot(X1,Y1,'Parent',subplot1,'LineWidth',2,'Color',[0 0 0]);

% Create ylabel
ylabel('Plunge (in)','FontWeight','demi','FontSize',12);

% Create title
title('Closed-Loop Response','FontWeight','demi','FontSize',12);

% Create subplot
subplot2 = subplot(3,1,2,'Parent',figure1,'YGrid','on','XGrid','on');
box('on');
hold('all');

% Create plot
plot(X1,Y2,'Parent',subplot2,'LineWidth',2,'Color',[0 0 0]);

% Create ylabel
ylabel('Pitch (deg)','FontWeight','demi','FontSize',12);

% Create subplot
subplot3 = subplot(3,1,3,'Parent',figure1,'YGrid','on','XGrid','on');
box('on');
hold('all');

% Create plot
plot(X1,Y3,'Parent',subplot3,'LineWidth',2,'Color',[0 0 0]);

% Create xlabel
xlabel('Time (s)','FontWeight','demi','FontSize',12);

% Create ylabel
ylabel('TE Pos (deg)','FontWeight','demi','FontSize',12);

author: Sable

repo: mcbench-benchmarks


text:
function rtk=udrcvbias(rtk,tt)

global glc
VAR_HWBIAS=1^2; PRN_HWBIAS=1e-6;

for i=1:glc.NFREQGLO

    j=rtk.il+i;

    if rtk.x(j)==0
        rtk=initx(rtk,1e-6,VAR_HWBIAS,j);
    elseif rtk.nfix>=rtk.opt.minfix&&rtk.sol.ratio>rtk.opt.thresar(1)
        rtk=initx(rtk,rtk.xa(j),rtk.Pa(j,j),j);
    else
        rtk.P(j,j)=rtk.P(j,j)+PRN_HWBIAS^2*abs(tt);
    end

end

return

repo: GINav


text:
function [cfg] = ft_multiplotER(cfg, varargin)

% FT_MULTIPLOTER plots the event-related potentials or event-related fields verus
% time, or the oscillatory activity (power or coherence) versus frequency. Multiple
% datasets can be overlayed. The plots are arranged according to their location
% specified in the layout.
%
% Use as
%   ft_multiplotER(cfg, data)
% or
%   ft_multiplotER(cfg, data, data2, ..., dataN)
%
% The data can be an event-related potential or field produced by
% FT_TIMELOCKANALYSIS, a power spectrum produced by FT_FREQANALYSIS or a coherence
% spectrum produced by FT_FREQDESCRIPTIVES.
%
% If you specify multiple datasets they should contain the same channels, etc.
%
% The configuration can have the following parameters:
%   cfg.parameter     = field to be plotted on y-axis, for example 'avg', 'powspctrm' or 'cohspctrm' (default is automatic)
%   cfg.maskparameter = field in the first dataset to be used for marking significant data
%   cfg.maskstyle     = style used for masking of data, 'box', 'thickness' or 'saturation' (default = 'box')
%   cfg.maskfacealpha = mask transparency value between 0 and 1
%   cfg.xlim          = 'maxmin' or [xmin xmax] (default = 'maxmin')
%   cfg.ylim          = 'maxmin', 'maxabs', 'zeromax', 'minzero', or [ymin ymax] (default = 'maxmin')
%   cfg.gradscale     = number, scaling to apply to the MEG gradiometer channels prior to display
%   cfg.magscale      = number, scaling to apply to the MEG magnetometer channels prior to display
%   cfg.channel       = Nx1 cell-array with selection of channels (default = 'all'), see FT_CHANNELSELECTION for details
%   cfg.refchannel    = name of reference channel for visualising connectivity, can be 'gui'
%   cfg.baseline      = 'yes', 'no' or [time1 time2] (default = 'no'), see FT_TIMELOCKBASELINE or FT_FREQBASELINE
%   cfg.trials        = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.axes          = string, 'yes' or 'no' whether to draw x- and y-axes for each graph (default = 'yes')
%   cfg.box           = string, 'yes' or 'no' whether to draw a box around each graph (default = 'no')
%   cfg.showlabels    = 'yes' or 'no' (default = 'no')
%   cfg.showoutline   = 'yes' or 'no' (default = 'no')
%   cfg.showscale     = 'yes' or 'no' (default = 'yes')
%   cfg.showcomment   = 'yes' or 'no' (default = 'yes')
%   cfg.comment       = string of text (default = date + limits)
%                       Add 'comment' to graph (according to COMNT in the layout)
%   cfg.limittext     = add user-defined text instead of cfg.comment, (default = cfg.comment)
%   cfg.fontsize      = font size of comment and labels (default = 8)
%   cfg.interactive   = 'yes' or 'no', make the plot interactive (default = 'yes')
%                       In a interactive plot you can select areas and produce a new
%                       interactive plot when a selected area is clicked. Multiple areas
%                       can be selected by holding down the SHIFT key.
%   cfg.renderer      = 'painters', 'zbuffer', ' opengl' or 'none' (default = [])
%   cfg.colorgroups   = 'sequential', 'allblack', 'labelcharN' (N = Nth character in label), 'chantype' or a vector
%                       with the length of the number of channels defining the groups (default = 'sequential')
%   cfg.linestyle     = linestyle/marker type, see options of the PLOT function (default = '-')
%                       can be a single style for all datasets, or a cell-array containing one style for each dataset
%   cfg.linewidth     = linewidth in points (default = 0.5)
%   cfg.linecolor     = color(s) used for plotting the dataset(s) (default = 'brgkywrgbkywrgbkywrgbkyw')
%                       alternatively, colors can be specified as Nx3 matrix of RGB values
%   cfg.directionality = '', 'inflow' or 'outflow' specifies for connectivity measures whether the
%                       inflow into a node, or the outflow from a node is plotted. The (default) behavior
%                       of this option depends on the dimord of the input data (see below).
%   cfg.layout        = specify the channel layout for plotting using one of the supported ways (see below).

% The following options for the scaling of the EEG, EOG, ECG, EMG, MEG and NIRS channels
% is optional and can be used to bring the absolute numbers of the different
% channel types in the same range (e.g. fT and uV). The channel types are determined
% from the input data using FT_CHANNELSELECTION.
%   cfg.eegscale      = number, scaling to apply to the EEG channels prior to display
%   cfg.eogscale      = number, scaling to apply to the EOG channels prior to display
%   cfg.ecgscale      = number, scaling to apply to the ECG channels prior to display
%   cfg.emgscale      = number, scaling to apply to the EMG channels prior to display
%   cfg.megscale      = number, scaling to apply to the MEG channels prior to display
%   cfg.gradscale     = number, scaling to apply to the MEG gradiometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.magscale      = number, scaling to apply to the MEG magnetometer channels prior to display (in addition to the cfg.megscale factor)
%   cfg.nirsscale     = number, scaling to apply to the NIRS channels prior to display
%   cfg.mychanscale   = number, scaling to apply to the channels specified in cfg.mychan
%   cfg.mychan        = Nx1 cell-array with selection of channels
%   cfg.chanscale     = Nx1 vector with scaling factors, one per channel specified in cfg.channel
%
% For the plotting of directional connectivity data the cfg.directionality option
% determines what is plotted. The default value and the supported functionality
% depend on the dimord of the input data. If the input data is of dimord
% 'chan_chan_XXX', the value of directionality determines whether, given the
% reference channel(s), the columns (inflow), or rows (outflow) are selected for
% plotting. In this situation the default is 'inflow'. Note that for undirected
% measures, inflow and outflow should give the same output. If the input data is of
% dimord 'chancmb_XXX', the value of directionality determines whether the rows in
% data.labelcmb are selected. With 'inflow' the rows are selected if the
% refchannel(s) occur in the right column, with 'outflow' the rows are selected if
% the refchannel(s) occur in the left column of the labelcmb-field. Default in this
% case is '', which means that all rows are selected in which the refchannel(s)
% occur. This is to robustly support linearly indexed undirected connectivity
% metrics. In the situation where undirected connectivity measures are linearly
% indexed, specifying 'inflow' or 'outflow' can result in unexpected behavior.
%
% The layout defines how the channels are arranged and what the size of each
% subplot is. You can specify the layout in a variety of ways:
%  - you can provide a pre-computed layout structure (see prepare_layout)
%  - you can give the name of an ascii layout file with extension *.lay
%  - you can give the name of an electrode file
%  - you can give an electrode definition, i.e. "elec" structure
%  - you can give a gradiometer definition, i.e. "grad" structure
% If you do not specify any of these and the data structure contains an
% electrode or gradiometer structure, that will be used for creating a
% layout. If you want to have more fine-grained control over the layout
% of the subplots, you should create your own layout file.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
% If you specify this option the input data will be read from a *.mat
% file on disk. This mat files should contain only a single variable named 'data',
% corresponding to the input structure. For this particular function, the
% data should be provided as a cell-array.
%
% See also FT_MULTIPLOTTFR, FT_SINGLEPLOTER, FT_SINGLEPLOTTFR, FT_TOPOPLOTER,
% FT_TOPOPLOTTFR, FT_PREPARE_LAYOUT

% Undocumented local options:
% cfg.layoutname
% cfg.preproc
% cfg.orient = landscape/portrait

% Copyright (C) 2003-2006, Ole Jensen
% Copyright (C) 2007-2011, Roemer van der Meij & Jan-Mathijs Schoffelen
% Copyright (C) 2012-2017, F.C. Donders Centre

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPERS NOTE: This code is organized in a similar fashion for multiplot/singleplot/topoplot
% and for ER/TFR and should remain consistent over those 6 functions.
% Section 1: general cfg handling that is independent from the data
% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data
% Section 3: select the data to be plotted and determine min/max range
% Section 4: do the actual plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Section 1: general cfg handling that is independent from the data

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

Ndata = length(varargin);
for i=1:Ndata
  % check if the input data is valid for this function
  varargin{i} = ft_checkdata(varargin{i}, 'datatype', {'timelock', 'freq'});
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'cohrefchannel', 'refchannel'});
cfg = ft_checkconfig(cfg, 'renamed', {'hlim', 'xlim'});
cfg = ft_checkconfig(cfg, 'renamed', {'matrixside',  'directionality'});
cfg = ft_checkconfig(cfg, 'renamed', {'vlim', 'ylim'});
cfg = ft_checkconfig(cfg, 'renamed', {'zparam', 'parameter'});
cfg = ft_checkconfig(cfg, 'renamed', {'graphcolor', 'linecolor'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedback', 'inflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'directionality', 'feedforward', 'outflow'});
cfg = ft_checkconfig(cfg, 'renamedval', {'zlim', 'absmax', 'maxabs'});
cfg = ft_checkconfig(cfg, 'unused',  {'cohtargetchannel'});
% cfg = ft_checkconfig(cfg, 'deprecated', {'xparam'});

% set the defaults
cfg.baseline       = ft_getopt(cfg, 'baseline', 'no');
cfg.trials         = ft_getopt(cfg, 'trials', 'all', 1);
cfg.xlim           = ft_getopt(cfg, 'xlim', 'maxmin');
cfg.ylim           = ft_getopt(cfg, 'ylim', 'maxmin');
cfg.comment        = ft_getopt(cfg, 'comment', date);
cfg.limittext      = ft_getopt(cfg, 'limittext', 'default');
cfg.axes           = ft_getopt(cfg, 'axes', 'yes');
cfg.showlabels     = ft_getopt(cfg, 'showlabels', 'no');
cfg.showoutline    = ft_getopt(cfg, 'showoutline', 'no');
cfg.showscale      = ft_getopt(cfg, 'showscale',   'yes');
cfg.showcomment    = ft_getopt(cfg, 'showcomment', 'yes');
cfg.box            = ft_getopt(cfg, 'box', 'no');
cfg.fontsize       = ft_getopt(cfg, 'fontsize', 8);
cfg.fontweight     = ft_getopt(cfg, 'fontweight');
cfg.interactive    = ft_getopt(cfg, 'interactive', 'yes');
cfg.renderer       = ft_getopt(cfg, 'renderer'); % let MATLAB decide on the default
cfg.orient         = ft_getopt(cfg, 'orient', 'landscape');
cfg.maskparameter  = ft_getopt(cfg, 'maskparameter');
cfg.colorgroups    = ft_getopt(cfg, 'colorgroups', 'condition');
cfg.linecolor      = ft_getopt(cfg, 'linecolor', 'brgkywrgbkywrgbkywrgbkyw');
cfg.linestyle      = ft_getopt(cfg, 'linestyle', '-');
cfg.linewidth      = ft_getopt(cfg, 'linewidth', 0.5);
cfg.maskstyle      = ft_getopt(cfg, 'maskstyle', 'box');
cfg.maskfacealpha  = ft_getopt(cfg, 'maskfacealpha', 1);
cfg.channel        = ft_getopt(cfg, 'channel', 'all');
cfg.directionality = ft_getopt(cfg, 'directionality', '');
cfg.figurename     = ft_getopt(cfg, 'figurename');
cfg.preproc        = ft_getopt(cfg, 'preproc');
cfg.tolerance      = ft_getopt(cfg, 'tolerance', 1e-5);
cfg.frequency      = ft_getopt(cfg, 'frequency', 'all'); % needed for frequency selection with TFR data
cfg.latency        = ft_getopt(cfg, 'latency', 'all'); % needed for latency selection with TFR data, FIXME, probably not used

% check for linestyle being a cell-array
if ischar(cfg.linestyle)
  cfg.linestyle = repmat({cfg.linestyle}, 1, Ndata);
end
% check it's length, and lengthen it if does not have enough styles in it
if (length(cfg.linestyle) > 1) && (length(cfg.linestyle) < Ndata )
  ft_error('either specify cfg.linestyle as a cell-array with one cell for each dataset, or only specify one linestyle')
elseif (length(cfg.linestyle) == 1)
  cfg.linestyle = repmat(cfg.linestyle, 1, Ndata);
end

% this is needed for the figure title and correct labeling of linecolor later on
if isfield(cfg, 'dataname') && ~isempty(cfg.dataname)
  dataname = cfg.dataname;
elseif isfield(cfg, 'inputfile') && ~isempty(cfg.inputfile)
  dataname = cfg.inputfile;
elseif nargin>1
  dataname = arrayfun(@inputname, 2:nargin, 'UniformOutput', false);
else
  dataname = {};
end

%% Section 2: data handling, this also includes converting bivariate (chan_chan and chancmb) into univariate data

for i=1:Ndata
  dtype{i}   = ft_datatype(varargin{i});
  hastime(i) = isfield(varargin{i}, 'time');
  hasfreq(i) = isfield(varargin{i}, 'freq');
end

% check if the input has consistent datatypes
if ~all(strcmp(dtype, dtype{1})) || ~all(hastime==hastime(1)) || ~all(hasfreq==hasfreq(1))
  ft_error('different datatypes are not allowed as input');
end
dtype   = dtype{1};
hastime = hastime(1);
hasfreq = hasfreq(1);

% Set x/y/parameter defaults according to datatype and dimord
switch dtype
  case 'timelock'
    xparam = 'time';
    if isfield(varargin{1}, 'trial')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'trial');
    elseif isfield(varargin{1}, 'individual')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'individual');
    elseif isfield(varargin{1}, 'avg')
      cfg.parameter = ft_getopt(cfg, 'parameter', 'avg');
    end
  case 'freq'
    if hastime && hasfreq
      xparam = 'time';  % further down the code will compute the average over selected frequencies
    else
      xparam = 'freq';
    end
    cfg.parameter = ft_getopt(cfg, 'parameter', 'powspctrm');
  case 'comp'
    % not supported
  otherwise
    % not supported
end

% check whether rpt/subj is present and remove if necessary
dimord = getdimord(varargin{1}, cfg.parameter);
dimtok = tokenize(dimord, '_');
hasrpt = any(ismember(dimtok, {'rpt' 'subj'}));

if ~hasrpt
  assert(isequal(cfg.trials, 'all') || isequal(cfg.trials, 1), 'incorrect specification of cfg.trials for data without repetitions');
else
  assert(~isempty(cfg.trials), 'empty specification of cfg.trials for data with repetitions');
end

% parse cfg.channel
if isfield(cfg, 'channel') && isfield(varargin{1}, 'label')
  cfg.channel = ft_channelselection(cfg.channel, varargin{1}.label);
elseif isfield(cfg, 'channel') && isfield(varargin{1}, 'labelcmb')
  cfg.channel = ft_channelselection(cfg.channel, unique(varargin{1}.labelcmb(:)));
end

% apply baseline correction
if ~strcmp(cfg.baseline, 'no')
  tmpcfg = keepfields(cfg, {'baseline', 'baselinetype', 'baselinewindow', 'demean', 'parameter', 'channel'});
  for i=1:Ndata
    % keep mask-parameter if it is set
    if ~isempty(cfg.maskparameter)
      tempmask = varargin{i}.(cfg.maskparameter);
    end
    if strcmp(dtype, 'timelock') && strcmp(xparam, 'time')
      varargin{i} = ft_timelockbaseline(tmpcfg, varargin{i});
    elseif strcmp(dtype, 'freq') && strcmp(xparam, 'time')
      varargin{i} = ft_freqbaseline(tmpcfg, varargin{i});
    elseif strcmp(dtype, 'freq') && strcmp(xparam, 'freq')
      ft_error('Baseline correction is not supported for spectra without a time dimension');
    else
      ft_warning('Baseline correction not applied, please set xparam');
    end
    % put mask-parameter back if it is set
    if ~isempty(cfg.maskparameter)
      varargin{i}.(cfg.maskparameter) = tempmask;
    end
  end
end

% channels SHOULD be selected here, as no interactive action produces a new multiplot
tmpcfg = keepfields(cfg, {'channel', 'showcallinfo', 'trials'});
if hasrpt
  tmpcfg.avgoverrpt = 'yes';
else
  tmpcfg.avgoverrpt = 'no';
end
if hastime && hasfreq
  tmpcfg.avgoverfreq = 'yes';       % average over the selected frequencies
  tmpcfg.frequency = cfg.frequency; % not to be confused with cfg.xlim or cfg.ylim
  tmpcfg.keepfreqdim = 'no';
else
  tmpcfg.avgoverfreq = 'no';
end
tmpvar = varargin{1};
[varargin{:}] = ft_selectdata(tmpcfg, varargin{:});
% restore the provenance information
[cfg, varargin{:}] = rollback_provenance(cfg, varargin{:});

if isfield(tmpvar, cfg.maskparameter) && ~isfield(varargin{1}, cfg.maskparameter)
  % the mask parameter is not present after ft_selectdata, because it is
  % not included in all input arguments. Make the same selection and copy
  % it over
  tmpvar = ft_selectdata(tmpcfg, tmpvar);
  varargin{1}.(cfg.maskparameter) = tmpvar.(cfg.maskparameter);
end

clear tmpvar tmpcfg dimord dimtok hastime hasfreq hasrpt

% ensure that the preproc specific options are located in the cfg.preproc
% substructure, but also ensure that the field 'refchannel' remains at the
% highest level in the structure. This is a little hack by JM because the field
% refchannel can relate to connectivity or to an EEg reference.

if isfield(cfg, 'refchannel'), refchannelincfg = cfg.refchannel; cfg = rmfield(cfg, 'refchannel'); end
cfg = ft_checkconfig(cfg, 'createsubcfg',  {'preproc'});
if exist('refchannelincfg', 'var'), cfg.refchannel  = refchannelincfg; end

if ~isempty(cfg.preproc)
  % preprocess the data, i.e. apply filtering, baselinecorrection, etc.
  fprintf('applying preprocessing options\n');
  if ~isfield(cfg.preproc, 'feedback')
    cfg.preproc.feedback = cfg.interactive;
  end
  for i=1:Ndata
    varargin{i} = ft_preprocessing(cfg.preproc, varargin{i});
  end
end

% handle the bivariate case
dimord = getdimord(varargin{1}, cfg.parameter);
if startsWith(dimord, 'chan_chan_') || startsWith(dimord, 'chancmb_')
  % convert the bivariate data to univariate and call this plotting function with univariate input
  cfg.originalfunction = 'ft_multiplotER';
  cfg.trials = 'all'; % trial selection has been taken care off
  bivariate_common(cfg, varargin{:});
  return
end

% apply channel-type specific scaling
fn = fieldnames(cfg);
tmpcfg = keepfields(cfg, fn(endsWith(fn, 'scale') | startsWith(fn, 'mychan') | strcmp(fn, 'channel') | strcmp(fn, 'parameter')));
for i=1:Ndata
  varargin{i} = chanscale_common(tmpcfg, varargin{i});
end

%% Section 3: select the data to be plotted and determine min/max range

% Read or create the layout that will be used for plotting
tmpcfg = keepfields(cfg, {'layout', 'rows', 'columns', 'commentpos', 'scalepos', 'elec', 'grad', 'opto', 'showcallinfo'});
cfg.layout = ft_prepare_layout(tmpcfg, varargin{1});

% Take the subselection of channels that is contained in the layout, this is the same in all datasets
[selchan, sellay] = match_str(varargin{1}.label, cfg.layout.label);

% Get physical min/max range of x, i.e. time or frequency
if strcmp(cfg.xlim, 'maxmin')
  % Find maxmin throughout all varargins:
  xmin = [];
  xmax = [];
  for i=1:Ndata
    xmin = nanmin([xmin varargin{i}.(xparam)]);
    xmax = nanmax([xmax varargin{i}.(xparam)]);
  end
else
  xmin = cfg.xlim(1);
  xmax = cfg.xlim(2);
end

% Get the index of the nearest bin, this is the same in all datasets
xminindx = nearest(varargin{1}.(xparam), xmin);
xmaxindx = nearest(varargin{1}.(xparam), xmax);
xmin = varargin{1}.(xparam)(xminindx);
xmax = varargin{1}.(xparam)(xmaxindx);
selx = xminindx:xmaxindx;
xval = varargin{1}.(xparam)(selx);

% Get physical y-axis range, i.e. of the parameter to be plotted
if ~isnumeric(cfg.ylim)
  % Find maxmin throughout all varargins
  ymin = [];
  ymax = [];
  for i=1:Ndata
    % Select the channels in the data that match with the layout and that are selected for plotting
    dat = varargin{i}.(cfg.parameter)(selchan,selx);
    ymin = min([ymin min(min(min(dat)))]);
    ymax = max([ymax max(max(max(dat)))]);
  end
  if strcmp(cfg.ylim, 'maxabs') % handle maxabs, make y-axis center on 0
    ymax = max([abs(ymax) abs(ymin)]);
    ymin = -ymax;
  elseif strcmp(cfg.ylim, 'zeromax')
    ymin = 0;
  elseif strcmp(cfg.ylim, 'minzero')
    ymax = 0;
  end
else
  ymin = cfg.ylim(1);
  ymax = cfg.ylim(2);
end

% Gather the data from all input data structures
datamatrix = zeros(Ndata, length(selchan), length(selx));
for i=1:Ndata
  datamatrix(i,:,:) = varargin{i}.(cfg.parameter)(selchan, selx);
end

if ~isempty(cfg.maskparameter)
  % one value for each channel-time point
  maskmatrix = varargin{1}.(cfg.maskparameter)(selchan, selx);
else
  % create an Nx0 matrix
  maskmatrix = zeros(length(selchan), 0);
end

chanX      = cfg.layout.pos(sellay, 1);
chanY      = cfg.layout.pos(sellay, 2);
chanWidth  = cfg.layout.width(sellay);
chanHeight = cfg.layout.height(sellay);
chanLabel  = cfg.layout.label(sellay);

%% Section 4: do the actual plotting

% determine the coloring of channels/conditions
linecolor = linecolor_common(cfg, varargin{:});

cla
hold on

if ischar(linecolor)
  set(gca, 'ColorOrder', char2rgb(linecolor))
elseif isnumeric(linecolor)
  set(gca, 'ColorOrder', linecolor)
end

% Plot the data
for m=1:length(selchan)
  mask = maskmatrix(m, :);
  if strcmp(cfg.maskstyle, 'difference')
    % combine the conditions in a single plot, highlight the difference
    yval = squeeze(datamatrix(:,m,:));
    % Clip out of bounds y values:
    yval(yval > ymax) = ymax;
    yval(yval < ymin) = ymin;
    ft_plot_vector(xval, yval, 'width', chanWidth(m), 'height', chanHeight(m), 'hpos', chanX(m), 'vpos', chanY(m), 'hlim', [xmin xmax], 'vlim', [ymin ymax], 'color', linecolor, 'style', cfg.linestyle{1}, 'linewidth', cfg.linewidth, 'axis', cfg.axes, 'highlight', mask, 'highlightstyle', cfg.maskstyle, 'facealpha', cfg.maskfacealpha);
  else
    % loop over the conditions, plot them on top of each other
    for i=1:Ndata
      yval = squeeze(datamatrix(i,m,:));
      % clip out of bounds y values:
      yval(yval > ymax) = ymax;
      yval(yval < ymin) = ymin;
      % select the color for the channel/condition
      if strcmp(cfg.colorgroups, 'condition')
        color = linecolor(i,:);
      else
        color = linecolor(m,:);
      end
      ft_plot_vector(xval, yval, 'width', chanWidth(m), 'height', chanHeight(m), 'hpos', chanX(m), 'vpos', chanY(m), 'hlim', [xmin xmax], 'vlim', [ymin ymax], 'color', color, 'style', cfg.linestyle{i}, 'linewidth', cfg.linewidth, 'axis', cfg.axes, 'highlight', mask, 'highlightstyle', cfg.maskstyle, 'facealpha', cfg.maskfacealpha);
    end
  end
end % for number of channels

% plot the layout, labels and outline
ft_plot_layout(cfg.layout, 'box', istrue(cfg.box), 'label', istrue(cfg.showlabels), 'outline', istrue(cfg.showoutline), 'point', 'no', 'mask', 'no', 'fontsize', cfg.fontsize, 'labelyoffset', 1.4*median(cfg.layout.height/2), 'labelalignh', 'center', 'chanindx', find(~ismember(cfg.layout.label, {'COMNT', 'SCALE'})) );

% write comment
if istrue(cfg.showcomment)
  % Add the colors of the different conditions/datasets to the comment
  colorLabels = [];
  if Ndata > 1 && strcmp(cfg.colorgroups, 'condition')
    for i=1:Ndata
      if ischar(linecolor)
        colorLabels = [colorLabels '\n' dataname{i} '='         linecolor(i)     ];
      elseif isnumeric(linecolor)
        colorLabels = [colorLabels '\n' dataname{i} '=' num2str(linecolor(i, :)) ];
      end
    end
  end
  cfg.comment = [cfg.comment colorLabels];

  k = find(strcmp('COMNT', cfg.layout.label));
  if ~isempty(k)
    limittext = cfg.limittext;
    if ~strcmp(limittext, 'default')
      comment = limittext;
    else
      comment = cfg.comment;
      comment = sprintf('%0s\nxlim=[%.3g %.3g]', comment, xmin, xmax);
      comment = sprintf('%0s\nylim=[%.3g %.3g]', comment, ymin, ymax);
    end
    ft_plot_text(cfg.layout.pos(k, 1), cfg.layout.pos(k, 2), sprintf(comment), 'FontSize', cfg.fontsize, 'FontWeight', cfg.fontweight);
    % plot an invisible box, the text itself is not sufficient to get the automatic scaling of the figures axes to include COMNT
    xy(1) = cfg.layout.pos(k, 1) - cfg.layout.width(k, 1)/2;
    xy(2) = cfg.layout.pos(k, 1) + cfg.layout.width(k, 1)/2;
    xy(3) = cfg.layout.pos(k, 2) - cfg.layout.height(k, 1)/2;
    xy(4) = cfg.layout.pos(k, 2) + cfg.layout.height(k, 1)/2;
    ft_plot_box(xy, 'edgecolor', 'none');
  end
end

% Plot scales
if istrue(cfg.showscale)
  l = find(strcmp(cfg.layout.label, 'SCALE'));
  if ~isempty(l)
    x = cfg.layout.pos(l,1);
    y = cfg.layout.pos(l,2);
    plotScales([xmin xmax], [ymin ymax], x, y, chanWidth(1), chanHeight(1), cfg)
  end
end

axis tight
axis off
hold off

% Make the axis a little wider when boxes are shown
if strcmp(cfg.box, 'yes')
  abc = axis;
  axis(abc + [-1 +1 -1 +1]*mean(abs(abc))/10)
end

% Set orientation for printing if specified
if ~isempty(cfg.orient)
  orient(gcf, cfg.orient);
end

% set the figure window title
if ~isempty(dataname)
  set(gcf, 'Name', sprintf('%d: %s: %s', double(gcf), mfilename, join_str(', ', dataname)));
else
  set(gcf, 'Name', sprintf('%d: %s', double(gcf), mfilename));
end
set(gcf, 'NumberTitle', 'off');

% Set renderer if specified
if ~isempty(cfg.renderer)
  set(gcf, 'renderer', cfg.renderer)
end

% Make the figure interactive
if strcmp(cfg.interactive, 'yes')
  % add the cfg/data/channel information to the figure under identifier linked to this axis
  ident                 = ['axh' num2str(round(sum(clock.*1e6)))]; % unique identifier for this axis
  set(gca, 'tag', ident);
  info                  = guidata(gcf);
  info.(ident).x        = cfg.layout.pos(:, 1);
  info.(ident).y        = cfg.layout.pos(:, 2);
  info.(ident).label    = cfg.layout.label;
  info.(ident).dataname = dataname;
  info.(ident).cfg      = cfg;
  info.(ident).varargin = varargin;
  guidata(gcf, info);

  set(gcf, 'WindowButtonUpFcn',  {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonUpFcn'});
  set(gcf, 'WindowButtonDownFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonDownFcn'});
  set(gcf, 'WindowButtonMotionFcn', {@ft_select_channel, 'multiple', true, 'callback', {@select_singleplotER}, 'event', 'WindowButtonMotionFcn'});
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous varargin
ft_postamble provenance
ft_postamble savefig

% add a menu to the figure, but only if the current figure does not have subplots
menu_fieldtrip(gcf, cfg, false);

if ~ft_nargout
  % don't return anything
  clear cfg
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotScales(hlim, vlim, hpos, vpos, width, height, cfg)

% the placement of all elements is identical
placement = {'hpos', hpos, 'vpos', vpos, 'width', width, 'height', height, 'hlim', hlim, 'vlim', vlim};

ft_plot_box([hlim vlim], placement{:}, 'edgecolor', 'k');

if hlim(1)<=0 && hlim(2)>=0
  ft_plot_line([0 0], vlim, placement{:}, 'color', 'k');
  ft_plot_text(0, vlim(1), '0  ', placement{:}, 'rotation', 90, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'middle', 'FontSize', cfg.fontsize);
end

if vlim(1)<=0 && vlim(2)>=0
  ft_plot_line(hlim, [0 0], placement{:}, 'color', 'k');
  ft_plot_text(hlim(1), 0, '0  ', placement{:}, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'middle', 'FontSize', cfg.fontsize);
end

ft_plot_text(hlim(1), vlim(1), [num2str(hlim(1), 3) ' '], placement{:}, 'rotation', 90, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top',    'FontSize', cfg.fontsize);
ft_plot_text(hlim(2), vlim(1), [num2str(hlim(2), 3) ' '], placement{:}, 'rotation', 90, 'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom', 'FontSize', cfg.fontsize);
ft_plot_text(hlim(1), vlim(1), [num2str(vlim(1), 3) ' '], placement{:}, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'bottom', 'FontSize', cfg.fontsize);
ft_plot_text(hlim(1), vlim(2), [num2str(vlim(2), 3) ' '], placement{:}, 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'top',    'FontSize', cfg.fontsize);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION which is called after selecting channels in case of cfg.interactive='yes'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function select_singleplotER(label, varargin)
% fetch cfg/data based on axis indentifier given as tag
ident       = get(gca,'tag');
info        = guidata(gcf);
cfg         = info.(ident).cfg;
datvarargin = info.(ident).varargin;
if ~isempty(label)
  cfg = removefields(cfg, 'inputfile');   % the reading has already been done and varargin contains the data
  cfg.baseline = 'no';                    % make sure the next function does not apply a baseline correction again
  cfg.channel = label;
  cfg.dataname = info.(ident).dataname;   % put data name in here, this cannot be resolved by other means
  cfg.trials = 'all';                     % trial selection has already been taken care of
  fprintf('selected cfg.channel = {%s}\n', join_str(', ', cfg.channel));
  % ensure that the new figure appears at the same position
  f = figure('position', get(gcf, 'Position'));
  ft_singleplotER(cfg, datvarargin{:});
end

author: spm


text:
function bool = iseven(number)
%ISEVEN true for even integer
%
%   Usage: bool = iseven(number)
%
%   Input parameters:
%       number  - number (or vector/matrix with numbers) to be tested
%
%   Output parameters:
%       bool    - true if the number is even, false else
%
%   See also: isodd

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

%% ===== Computation =====================================================
% Create answer
bool = false( size(number) );
% Look for even numbers, use bitget to overcome a mod() bug, see
% https://bit.ly/1wcNYBI
bool(bitget(number,1)==0) = true;

author: sfstoolbox

repo: sfs-matlab


text:
function[invmat]=matinv(mat,str)
%MATINV  Fast inversion of arrays of small matrices.
%
%   MATINV is a low-level function called by POLYSMOOTH.
%
%   Let MAT be an array of K different M x M matrices A1,A2,...,AK.
%   INV=MATINV(MAT) then returns an array of N inverse matrices.
%
%   If MAT has dimensions K1 x K2 x .... M x M, then MATINV returns an
%   array of the same size containing the inverses of the M x M matrices.
%
%   For example, MAT could be 10 x 10 x 4 x 4, in which case the inverses
%   of one hundred 4 x 4 matrices are found.
%
%   MAT can have any dimensionality so long as the matrices to be inverted
%   occupy the last two dimensions.  The last dimension is interpreted as
%   "columns" and the second to last as "rows."
%
%   Note that MATINV only works matrices with M=2 through M=8.
%   ____________________________________________________________
%
%   Algorithms
%
%   MATINV can use either of two different algorithms.  This is specified
%   with INV=MATINV(MAT,STR).
%
%   MATINV(MAT,'direct') uses algebraic expressions for 2 x 2 and 3 x 3
%   matrix inverses together with Boltz's block diagonal recursion formula.
%   For details, see the following links
%
%        http://mathworld.wolfram.com/MatrixInverse.html
%        http://en.wikipedia.org/wiki/Invertible_matrix.
%
%   MATINV(MAT,'loop') uses Matlab's INV function together with a
%   straightforward loop.
%
%   The direct algorithm, which is the default, can be much faster when
%   MAT is large and the dimension to be inverted is small.
%   ____________________________________________________________
%
%   See also MATMULT.
%
%   'matinv --t' runs some tests.
%
%   Usage: inv=matinv(mat);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2020 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(mat, '--t')
    matinv_test,return
end

ndims=lnsd(mat);
K=size(mat,ndims);

if K~=size(mat,ndims-1)
    error('The last two dimensions of MAT must be the same for it to be invertible.')
end

if nargin==1
    str='direct';
end

if K>12
    error('MATINV number of dimensions of matrix should be no more than M=12.')
end

sizemat=size(mat);
if length(sizemat)>2
    %[mat,index]=matinv_strip(mat);

    if ~isempty(strfind(str,'dir'))
        invmat=matinv_direct(mat,ndims,K);
    elseif ~isempty(strfind(str,'loo'))
        invmat=matinv_loop(mat);
    end

else
    invmat=inv(mat);
end

function[mat,index]=matinv_strip(mat)

sizemat=size(mat);
mat=reshape(mat,prod(sizemat(1:end-2)),sizemat(end-1),sizemat(end));
bool=~isnan(sum(sum(mat,3),2));
index=find(bool);
mat=mat(index,:,:);

function[invmat]=matinv_loop(mat)

sizemat=size(mat);
mat=reshape(mat,prod(sizemat(1:end-2)),sizemat(end-1),sizemat(end));
invmat=0*mat;

%[lastmsg,lastid]=lastwarn;
warning('off','MATLAB:illConditionedMatrix');
warning('off','MATLAB:nearlySingularMatrix');
%warning('off','MATLAB:singularMatrix');
for i=1:size(mat,1)
     invmat(i,:,:)=inv(squeeze(mat(i,:,:)));
end

invmat=reshape(invmat,sizemat);
warning('on','MATLAB:illConditionedMatrix');
warning('on','MATLAB:nearlySingularMatrix');
warning('on','MATLAB:singularMatrix');

function[invmat]=matinv_direct(mat,ndims,K)

%if K==1
%    invmat=1./mat;
if K==2
    invmat=matinv_twoxtwo(mat,ndims);
elseif K==3
    invmat=matinv_threexthree(mat,ndims);
else
    invmat=matinv_block(mat,ndims,K);
end

function[invmat]=matinv_twoxtwo(mat,ndims)
ac=vindex(mat,1,ndims);
bd=vindex(mat,2,ndims);

a=vindex(ac,1,ndims-1);
c=vindex(ac,2,ndims-1);
b=vindex(bd,1,ndims-1);
d=vindex(bd,2,ndims-1);

deta=a.*d-b.*c;
a=a./deta;
b=b./deta;
c=c./deta;
d=d./deta;

ac=0*ac;
bd=0*bd;

ac=vindexinto(ac,d,1,ndims-1);
ac=vindexinto(ac,-c,2,ndims-1);

bd=vindexinto(bd,-b,1,ndims-1);
bd=vindexinto(bd,a,2,ndims-1);

invmat=zeros(size(mat));
invmat=vindexinto(invmat,ac,1,ndims);
invmat=vindexinto(invmat,bd,2,ndims);

function[mat]=matinv_threexthree(mat,ndims)

c1=vindex(mat,1,ndims);
c2=vindex(mat,2,ndims);
c3=vindex(mat,3,ndims);

a11=vindex(c1,1,ndims-1);
a21=vindex(c1,2,ndims-1);
a31=vindex(c1,3,ndims-1);

a12=vindex(c2,1,ndims-1);
a22=vindex(c2,2,ndims-1);
a32=vindex(c2,3,ndims-1);

a13=vindex(c3,1,ndims-1);
a23=vindex(c3,2,ndims-1);
a33=vindex(c3,3,ndims-1);

b11=a22.*a33-a23.*a32;
b21=a23.*a31-a21.*a33;
b31=a21.*a32-a22.*a31;

b12=a13.*a32-a12.*a33;
b22=a11.*a33-a13.*a31;
b32=a12.*a31-a11.*a32;

b13=a12.*a23-a13.*a22;
b23=a13.*a21-a11.*a23;
b33=a11.*a22-a12.*a21;

deta=a11.*b11+a12.*b21+a13.*b31;
%This looks funny because b21 is minus

b11=b11./deta;
b21=b21./deta;
b31=b31./deta;

b12=b12./deta;
b22=b22./deta;
b32=b32./deta;

b13=b13./deta;
b23=b23./deta;
b33=b33./deta;

clear deta

c1=vindexinto(c1,b11,1,ndims-1);
c1=vindexinto(c1,b21,2,ndims-1);
c1=vindexinto(c1,b31,3,ndims-1);

c2=vindexinto(c2,b12,1,ndims-1);
c2=vindexinto(c2,b22,2,ndims-1);
c2=vindexinto(c2,b32,3,ndims-1);

c3=vindexinto(c3,b13,1,ndims-1);
c3=vindexinto(c3,b23,2,ndims-1);
c3=vindexinto(c3,b33,3,ndims-1);

mat=vindexinto(mat,c1,1,ndims);
mat=vindexinto(mat,c2,2,ndims);
mat=vindexinto(mat,c3,3,ndims);

function[invmat]=matinv_block(mat,ndims,K)

%This keeps the two matrices about the same size, which seems the
%fastest option in tests

i1=1:floor(K/2);
i2=floor(K/2)+1:K;

%i1,i2

ac=vindex(mat,i1,ndims);
bd=vindex(mat,i2,ndims);

a=vindex(ac,i1,ndims-1);
c=vindex(ac,i2,ndims-1);
b=vindex(bd,i1,ndims-1);
d=vindex(bd,i2,ndims-1);

%vsize(a,b,c,d)

%Recursion
ainv=matinv(a);

dcab=matinv(d-matmult(matmult(c,ainv,ndims-1),b,ndims-1));

bnew=-matmult(matmult(ainv,b,ndims-1),dcab,ndims-1);
anew=ainv+matmult(matmult(-bnew,c,ndims-1),ainv,ndims-1);
cnew=-matmult(matmult(dcab,c,ndims-1),ainv,ndims-1);
dnew=dcab;

ac=vindexinto(ac,anew,i1,ndims-1);
ac=vindexinto(ac,cnew,i2,ndims-1);

bd=vindexinto(bd,bnew,i1,ndims-1);
bd=vindexinto(bd,dnew,i2,ndims-1);

invmat=zeros(size(mat));
invmat=vindexinto(invmat,ac,i1,ndims);
invmat=vindexinto(invmat,bd,i2,ndims);

function[]=matinv_test
disp('Testing that direct and looping algorithms match for non-singular matrices.')

rng(0);

for n=2:12

    mat=randn(n);

    invmat=inv(mat);
    invmat2=matinv(mat);
    disp(['MATINV testing case of ' int2str(n) 'x' int2str(n) ' matrices.'])
    reporttest('MATINV with single matrix',aresame(invmat,invmat2,1e-3))
    mat=randn(10000,n,n);
    tic
    invmat=matinv(mat,'loop');
    t1=toc;

    tic
    invmat2=matinv(mat,'direct');
    t2=toc;

    tol=1e-1;
    bool=false(size(mat,1),1);
    for i=1:size(mat)
        bool(i)=rcond(squeeze(mat(i,:,:)))>=tol;
    end

    vindex(invmat,invmat2,find(bool),1);

    reporttest('MATINV with 10000 random matrices',aresame(invmat,invmat2,1e-2))
    disp(['MATINV was ' num2str(t1./t2) ' times faster than INV.'])
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%% Ex. 2 For loop: Utility of the dummy index

b = 3;
for k = 1:5
 b^k
end
%Output:
% 3
% 9
% 27
% 81
% 243

%Remark: The outputs are 3^1, 3^2, 3^3, 3^4, and 3^5. the value of "k" keeps changing as
% we go through the loop

author: TheAlgorithms

repo: MATLAB-Octave


text:
function geometry_test061 ( )

%*****************************************************************************80
%
%% TEST061 tests PLANE_NORMAL_BASIS_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST061\n' );
  fprintf ( 1, '  PLANE_NORMAL_BASIS_3D, given a plane in\n' );
  fprintf ( 1, '    point, normal form (P,N), finds two unit\n' );
  fprintf ( 1, '    vectors Q and R that "lie" in the plane\n' );
  fprintf ( 1, '    and are mutually orthogonal.\n' );

  for test = 1 : test_num

    [ pp(1:3,1), seed ] = r8vec_uniform_01 ( 3, seed );

    [ normal(1:3,1), seed ] = r8vec_uniform_01 ( 3, seed );
    t = norm ( normal );
    normal = normal / t;

    [ pq, pr ] = plane_normal_basis_3d ( pp, normal );

    if ( test == 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Data for test 1:\n' );
      fprintf ( 1, '\n' );
      r8vec_print ( 3, pp, '  Point PP:' );
      r8vec_print ( 3, normal, '  Normal vector N:' );
      r8vec_print ( 3, pq, '  Vector PQ:' );
      r8vec_print ( 3, pr, '  Vector PR:' );
    end

    b(1,1) = normal(1:3,1)' * normal(1:3,1);
    b(1,2) = normal(1:3,1)' * pq(1:3,1);
    b(1,3) = normal(1:3,1)' * pr(1:3,1);

    b(2,1) = pq(1:3,1)' * normal(1:3,1);
    b(2,2) = pq(1:3,1)' * pq(1:3,1);
    b(2,3) = pq(1:3,1)' * pr(1:3,1);

    b(3,1) = pr(1:3,1)' * normal(1:3,1);
    b(3,2) = pr(1:3,1)' * pq(1:3,1);
    b(3,3) = pr(1:3,1)' * pr(1:3,1);

    r8mat_print ( 3, 3, b, '  Dot product matrix:' );

  end

  return
end


text:
function [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,options)
% spike inference using a constrained deconvolution approach:
%      min      sum(sp)
%    c,sp,b,c1
%      subject to: sp >= 0
%                   b >= 0
%                  G*c = sp
%                   c1 >= 0
%           ||y-b-c - c_in|| <= sn*sqrt(T)

%   Variables:
%   y:      raw fluorescence data (vector of length(T))
%   c:      denoised calcium concentration (Tx1 vector)
%   b:      baseline concentration (scalar)
%  c1:      initial concentration (scalar)
%   g:      discrete time constant(s) (scalar or 2x1 vector)
%  sn:      noise standard deviation (scalar)
%  sp:      spike vector (Tx1 vector)

%   USAGE:
%   [c,b,c1,g,sn,sp] = constrained_foopsi(y,b,c1,g,sn,OPTIONS)
%   The parameters b,cin,g,sn can be given or else are estimated from the data

%   OPTIONS: (stuct for specifying options)
%         p: order for AR model, used when g is not given (default 2)
%    method: methods for performing spike inference
%   available methods: 'dual' uses dual ascent
%                       'cvx' uses the cvx package available from cvxr.com (default)
%                      'lars' uses the least regression algorithm
%                     'spgl1' uses the spgl1 package available from
%                     math.ucdavis.edu/~mpf/spgl1/  (can be faster)
%   bas_nonneg:   flag for setting the baseline lower bound. if 1, then b >= 0 else b >= min(y)
%   noise_range:  frequency range over which the noise power is estimated. Default [Fs/4,Fs/2]
%   noise_method: method to average the PSD in order to obtain a robust noise level estimate
%   lags:         number of extra autocovariance lags to be considered when estimating the time constants
%   resparse:     number of times that the solution is resparsened (default 0). Currently available only with methods 'cvx', 'spgl'
%   fudge_factor: scaling constant to reduce bias in the time constant estimation (default 1 - no scaling)

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

defoptions.p = 2;
defoptions.method = 'cvx';
defoptions.bas_nonneg = 1;              % nonnegativity option for baseline estimation
defoptions.noise_range = [0.25,0.5];    % frequency range over which to estimate the noise
defoptions.noise_method = 'logmexp';    % method for which to estimate the noise level
defoptions.lags = 5;                    % number of extra lags when computing the AR coefficients
defoptions.resparse = 0;                % number of times to re-sparse solution
defoptions.fudge_factor = 1;            % fudge factor for time constants

if nargin < 6
    options = defoptions;
    if nargin < 5
        sn = [];
        if nargin < 4
            g = [];
            if nargin < 3
                c1 = [];
                if nargin < 2
                    b = [];
                end
            end
        end
    end
end

if ~isfield(options,'p');  options.p = defoptions.p;  end
if ~isfield(options,'method'); options.method = defoptions.method; end
if ~isfield(options,'bas_nonneg'); options.bas_nonneg = defoptions.bas_nonneg; end
if ~isfield(options,'noise_range'); options.noise_range = defoptions.noise_range; end
if ~isfield(options,'noise_method'); options.noise_method = defoptions.noise_method; end
if ~isfield(options,'lags'); options.lags = defoptions.lags; end
if ~isfield(options,'resparse'); options.resparse = defoptions.resparse; end
if ~isfield(options,'fudge_factor'); options.fudge_factor = defoptions.fudge_factor; end

method = options.method;
if isempty(b);
    bas_est = 1;
else
    bas_est = 0;
end
if isempty(c1)
    c1_est = 1;
else
    c1_est = 0;
end

y = y(:);
T = length(y);
y_full = y;
mis_data = isnan(y);
E = speye(T);
E(mis_data,:) = [];

if any(mis_data)
    y_full(mis_data) = interp1(find(~mis_data),y(~mis_data),find(mis_data));
end

if isempty(sn)
    sn = GetSn(y_full,options.noise_range,options.noise_method);
end
if isempty(g)
    g = estimate_time_constants(y_full,options.p,sn,options.lags);
    while max(abs(roots([1,-g(:)']))>1) && options.p < 5
        warning('No stable AR(%i) model found. Checking for AR(%i) model \n',options.p,options.p+1);
        options.p = options.p + 1;
        g = estimate_time_constants(y,options.p,sn,options.lags);
    end
    if options.p == 5
        g = 0;
    end
    %fprintf('Stable AR(%i) model found \n',options.p);
    % re-adjust time constant values
    rg = roots([1;-g(:)]);
    if ~isreal(rg); rg = real(rg) + .001*randn(size(rg)); end
    rg(rg>1) = 0.95 + 0.001*randn(size(rg(rg>1)));
    rg(rg<0) = 0.15 + 0.001*randn(size(rg(rg<0)));
    pg = poly(options.fudge_factor*rg);
    g = -pg(2:end);
end
if options.bas_nonneg  % lower bound for baseline
    b_lb = 0;
else
    b_lb = min(y);
end

if strcmpi(method,'dual'); method = 'dual';
elseif strcmpi(method,'cvx'); method = 'cvx';
elseif strcmpi(method,'lars'); method = 'lars';
elseif strcmpi(method,'spgl1'); method = 'spgl1';
else fprintf('Invalid choice of method. Using CVX \n'); method = 'cvx';
end

if strcmpi(method,'dual') && any(mis_data)
    warning('Dual method does not support missing data. Switching to CVX');
    method = 'cvx';
end

if options.resparse > 0 && (strcmpi(method,'dual') || strcmpi(method,'lars'))
    warning('Resparsening is not supported with chosen method. Switching to CVX');
    method = 'cvx';
end

pathCell = regexp(path, pathsep, 'split');
g = g(:);
G = spdiags(ones(T,1)*[-g(end:-1:1)',1],-length(g):0,T,T);
gd = max(roots([1,-g']));  % decay time constant for initial concentration
gd_vec = gd.^((0:T-1)');

switch method
    case 'dual'
         v = G'*ones(T,1);
        thr = sn*sqrt(T-sum(mis_data));
        if bas_est; b = 0; end
        if c1_est; c1 = 0; end
        myfun = @(Ald) lagrangian_temporal_gradient(Ald,thr^2,y(~mis_data)-b-c1*gd_vec(~mis_data),bas_est,c1_est);
        c = [G\max(G*y,0);zeros(bas_est);zeros(c1_est)];
        options_dual = optimset('GradObj','On','Display','Off','Algorithm','interior-point','TolX',1e-8);
        ld_in = 10;
        [ld,~,flag] = fmincon(myfun,ld_in,[],[],[],[],0,[],[],options_dual);
        if (flag == -2) || (flag == -3)
            warning('Problem seems unbounded or infeasible. Try a different method.');
        end
        if bas_est; b = c(T+bas_est); end
        if c1_est; c1 = c(end); end
        c = c(1:T);
        sp = G*c;
    case 'cvx'
        onPath = ~isempty(which('cvx_begin'));
        if onPath
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                [c(:,rep),bas(rep),cin(rep)] = cvx_foopsi(y,b,c1,sn,b_lb,g,w_,~mis_data);
                sp(:,rep) = G*c(:,rep);
                w_ = 1./(max(sp(:,rep),0) + 1e-8);
            end
            sp(sp<1e-6) = 0;
            c = G\sp;
            b = bas;
            c1 = cin;
        else
            error('CVX does not appear to be on the MATLAB path. It can be downloaded from cvxr.com \n');
        end
    case 'lars'
         Ginv = E*[full(G\speye(T)),ones(T,bas_est),gd_vec*ones(1,c1_est)];
         if bas_est; b = 0; end
         if c1_est; c1 = 0; end
         [~, ~, spikes, ~, ~] = lars_regression_noise(y(~mis_data)-b_lb*bas_est - b - c1*gd_vec(~mis_data), Ginv, 1, sn^2*(T-sum(mis_data)));
         sp = spikes(1:T);
         b = (spikes(T+bas_est)+b_lb)*bas_est + b*(1-bas_est);
         c1 = spikes(end)*c1_est + c1*(1-c1_est);
         c = G\sp;
    case 'spgl1'
        onPath = ~isempty(which('spgl1'));
        if onPath
            Gx = @(x,mode) G_inv_mat(x,mode,T,g,gd_vec,bas_est,c1_est,E);
            c = zeros(T,1+options.resparse);
            sp = zeros(T,1+options.resparse);
            bas = zeros(1+options.resparse,1);
            cin = zeros(1+options.resparse,1);
            w_ = ones(T,1);
            for rep = 1:options.resparse+1
                if bas_est; b = 0; w_ = [w_;1e-10]; end
                if c1_est; c1 = 0; w_ = [w_;1e-10]; end
                options_spgl = spgSetParms('project',@NormL1NN_project ,'primal_norm', @NormL1NN_primal,'dual_norm',@NormL1NN_dual,'verbosity',0,'weights',w_);
                [spikes,r,~,~] = spg_bpdn( Gx, y(~mis_data)-b_lb*bas_est - (1-bas_est)*b-(1-c1_est)*c1*gd_vec(~mis_data), sn*sqrt(T-sum(mis_data)),options_spgl);
                c(:,rep) = G\spikes(1:T); %Gx([spikes(1:T);0],1);                                  %% calcium signal
                bas(rep) = b*(1-bas_est) + bas_est*spikes(T+bas_est)+b_lb*bas_est;       %% baseline
                cin(rep) = c1*(1-c1_est) + c1_est*spikes(end);
                sp(:,rep) = spikes(1:T);                                           %% spiking signal
                w_ = 1./(spikes(1:T)+1e-8);
            end
            b = bas;
            c1 = cin;
            %sn = norm(r)/sqrt(T);
        else
            error('SPGL1 does not appear to be on the MATLAB path. It can be downloaded from math.ucdavis.edu/~mpf/spgl1 \n');
        end
end

    function sn = GetSn(Y,range_ff,method)
        % estimate noise level with a power spectral density method
        L=length(Y);
%         if ~isempty(which('pmtm'))
%             [psd_Y,ff] = pmtm(Y,5/2,1000,1);
%         end
        if ~isempty(which('pwelch'));
            [psd_Y,ff]=pwelch(Y,round(L/8),[],1000,1);
        else
            xdft = fft(Y);
            xdft = xdft(1:round(L/2)+1);
            psd_Y = (1/L) * abs(xdft).^2;
            ff = 0:1/L:1/2;
            psd_Y(2:end-1) = 2*psd_Y(2:end-1);
        end
        ind=ff>range_ff(1);
        ind(ff>range_ff(2))=0;
        switch method
            case 'mean'
                sn=sqrt(mean(psd_Y(ind)/2));
            case 'median'
                sn=sqrt(median(psd_Y(ind)/2));
            case 'logmexp'
                sn = sqrt(exp(mean(log(psd_Y(ind)/2))));
        end
    end

    function g = estimate_time_constants(y,p,sn,lags)
        % estimate time constants from the sample autocovariance function

        lags = lags + p;
        if ~isempty(which('xcov')) %signal processing toolbox
            xc = xcov(y,lags,'biased');
        else
            ynormed = (y - mean(y));
            xc = nan(lags + 1, 1);
            for k = 0:lags
                xc(k + 1) = ynormed(1 + k:end)' * ynormed(1:end - k);
            end
            xc = [flipud(xc(2:end)); xc] / numel(y);
        end
        xc = xc(:);
        A = toeplitz(xc(lags+(1:lags)),xc(lags+(1:p))) - sn^2*eye(lags,p);
        try
            g = pinv(A)*xc(lags+2:end);
        catch
            g = 0;
        end
    end

    function [f,grad] = lagrangian_temporal_gradient(Al,thr,y_raw,bas_flag,c1_flag)
        options_qp = optimset('Display','Off','Algorithm','interior-point-convex');
        H = [speye(T),ones(T,bas_flag),gd_vec*ones(1,c1_flag);...
            ones(bas_flag,T),T*ones(bas_flag),(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag);...
            (gd_vec*ones(1,c1_flag))',(1-gd^T)/(1-gd)*ones(c1_flag,bas_flag),(1-gd^(2*T))/(1-gd^2)*ones(c1_flag,c1_flag)];
        Ay = [y_raw;sum(y_raw)*ones(bas_flag);gd_vec'*y_raw*ones(c1_flag)];
        c = quadprog(2*Al(1)*H,[v;zeros(bas_flag+c1_flag,1)]-2*Al(1)*Ay,[-G,sparse(T,bas_flag+c1_flag);sparse(bas_flag+c1_flag,T),-speye(bas_flag+c1_flag)]...
            ,[sparse(T,1);-b_lb*ones(bas_flag);zeros(c1_flag)],[],[],[],[],c,options_qp);
        f = v'*c(1:T);
        grad = [sum((c(1:T)-y_raw + c(T+bas_flag)*bas_flag + c(end)*gd_vec*c1_flag).^2)-thr];
        f = f + Al(:)'*grad;
    end

    function b = G_inv_mat(x,mode,NT,gs,gd_vec,bas_flag,c1_flag,Emat)
        if mode == 1
            b = filter(1,[1;-gs(:)],x(1:NT)) + bas_flag*x(NT+bas_flag) + c1_flag*gd_vec*x(end);
            b = Emat*b;
           %b = G\x(1:NT) + x(NT+bas_flag)*bas_flag + x(end)*c1_flag;
        elseif mode == 2
            x = Emat'*x;
            b = [flipud(filter(1,[1;-gs(:)],flipud(x)));ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)];
           %b = [G'\x;ones(bas_flag,1)*sum(x);ones(c1_flag,1)*(gd_vec'*x)] ;
        end
    end

end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
classdef TestPutText
    %TestPutText

    methods (Static)
        function test_simple
            im = 255*ones(64,64,3,'uint8');
            out = cv.putText(im, 'foo', [5,30]);
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_options
            im = zeros([128,128,3],'uint8');
            out = cv.putText(im, 'foo', [40,70], ...
                'FontFace','HersheyComplex', 'FontStyle','Italic',...
                'Color',[255 0 0], 'Thickness',2, 'LineType','AA');
            validateattributes(out, {class(im)}, {'size',size(im)});
        end

        function test_error_argnum
            try
                cv.putText();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [out,x] =  symmetric_bc( x, nx)

    if(x < 0)
   border = nx - 1;
       xx = -x;
        n  = mod((xx/border) , 2);

        if ( n ) x = border - (mod( xx ,border) );
        else x = mod(xx , border);
        out = 1;
        end
        elseif ( x >= nx )
         border = nx - 1;
         n = mod((x/border) , 2);
        end
        if ( n ) x = border - ( mod(x , border) );
        else x = mod(x , border);
        out = 1;

        end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm