text:
function [L, Nc] = discretize(S, N)
%DISCRETIZE   Discretize a SPINOP2.
%   [L, NC] = DISCRETIZE(S, N) uses a Fourier spectral method in coefficient
%   space to discretize the SPINOP2 S with N grid points in each direction. L is
%   the linear part, a N^2xN^2 diagonal matrix stored as a NxN matrix, and NC is
%   the differentiation term of the nonlinear part (and hence is linear).
%
% Remark: DISCRETIZE will fail to discretize SPINOP2 objects which are not of
%         the right form. See HELP/SPINOP2.
%
% See also SPINOP2.

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

%% Set-up:

% Get the domain DOM, the linear part LFUN, the nonlinear part NFUN, and the
% number of variables NVARS from S:
dom = S.domain;
funcL = S.lin;
nVars = nargin(funcL);

% Get the variables of the workspace:
func = functions(funcL);
wrk = func.workspace{1};
names = fieldnames(wrk);
if ( isempty(names) == 0 )
    lengthNames = size(names, 1);
    for k = 1:lengthNames
        eval(sprintf('%s = wrk.(names{k});', names{k}));
    end
end

% Create a CHEBOPPREF object with TRIGSPEC discretization:
pref = cheboppref();
pref.discretization = @trigspec;

%% Discretize the linear part:

% Second-order Fourier differentiation matrix with TRIGSPEC (sparse diagonal
% matrix):
D2 = trigspec.diffmat(N,2)*(2*pi/(dom(2) - dom(1)))^2;
if ( mod(N,2) == 0 )
    D2 = fftshift(D2);
else
    D2 = ifftshift(D2);
end

% Look for 'laplacian'/'lap', 'biharmonic'/'biharm', 'triharmonic'/'triharm',
% 'quadharmonic'/'quadharm' or 'quintharmomic'/'quintharm':
strL = func2str(funcL);
isLap = isempty(strfind(strL,'laplacian')) && isempty(strfind(strL,'lap'));
isLap = ~isLap;
isBih = isempty(strfind(strL,'biharmonic')) && isempty(strfind(strL,'biharm'));
isBih = ~isBih;
isTrih = isempty(strfind(strL,'triharmonic')) ...
    && isempty(strfind(strL,'triharm'));
isTrih = ~isTrih;
isQuadh = isempty(strfind(strL,'quadharmonic')) ...
    && isempty(strfind(strL,'quadharm'));
isQuadh = ~isQuadh;
isQuinth = isempty(strfind(strL,'quintharmonic')) ...
    && isempty(strfind(strL,'quintharm'));
isQuinth = ~isQuinth;

% NxN identity matrix for the Kronecker products:
I = eye(N);

% Construct the Laplacian operator -- needed for all the operators:
if ( isLap || isBih || isTrih || isQuadh || isQuinth )

    % Compute the N^2xN^2 Laplacian with KRON:
    lapmat = kron(I, D2) + kron(D2, I);

    % Create a NxN matrix with the diagonal of the N^2xN^2 Laplacian:
    lapmat = reshape(full(diag(lapmat)), N, N);

else
    lapmat = 0;
end

% The linear part has a B*biharmonic(u) term:
if ( isBih == 1 )

    % Pointwise multiplication since we only store the diagonal elements:
    bihmat = lapmat.^2;

else
    bihmat = 0;
end

% The linear part has a C*triharmonic(u) term:
if ( isTrih == 1 )

    % Pointwise multiplication since we only store the diagonal elements:
    trihmat = lapmat.^3;

else
    trihmat = 0;
end

% The linear part has a D*quadharmonic(u) term:
if ( isQuadh == 1 )

    % Pointwise multiplication since we only store the diagonal elements:
    quadhmat = lapmat.^4;

else
    quadhmat = 0;
end

% The linear part has a E*quintharmonic(u) term:
if ( isQuinth == 1 )

    % Pointwise multiplication since we only store the diagonal elements:
    quinthmat = lapmat.^5;

else
    quinthmat = 0;
end

% Convert to a string and initialize L:
strL = func2str(funcL);
L = [];

% Get the constants A in front of the Laplacians:
str = strrep(strL, 'laplacian', '');
str = strrep(str, 'lap', '');
str = strrep(str, 'biharmonic', '0*');
str = strrep(str, 'biharm', '0*');
str = strrep(str, 'triharmonic', '0*');
str = strrep(str, 'triharm', '0*');
str = strrep(str, 'quadharmonic', '0*');
str = strrep(str, 'quadharm', '0*');
str = strrep(str, 'quintharmonic', '0*');
str = strrep(str, 'quintharm', '0*');
func = eval(str);
inputs = cell(1, nVars);
for k = 1:nVars
   inputs{k} = 1;
end
A = feval(func, inputs{:});

% Get the constants B in front of the biharmonic operators:
str = strrep(strL, 'laplacian', '0*');
str = strrep(str, 'lap', '0*');
str = strrep(str, 'biharmonic', '');
str = strrep(str, 'biharm', '');
str = strrep(str, 'triharmonic', '0*');
str = strrep(str, 'triharm', '0*');
str = strrep(str, 'quadharmonic', '0*');
str = strrep(str, 'quadharm', '0*');
str = strrep(str, 'quintharmonic', '0*');
str = strrep(str, 'quintharm', '0*');
func = eval(str);
B = feval(func, inputs{:});

% Get the constants C in front of the triharmonic operators:
str = strrep(strL, 'laplacian', '0*');
str = strrep(str, 'lap', '0*');
str = strrep(str, 'biharmonic', '0*');
str = strrep(str, 'biharm', '0*');
str = strrep(str, 'triharmonic', '');
str = strrep(str, 'triharm', '');
str = strrep(str, 'quadharmonic', '0*');
str = strrep(str, 'quadharm', '0*');
str = strrep(str, 'quintharmonic', '0*');
str = strrep(str, 'quintharm', '0*');
func = eval(str);
C = feval(func, inputs{:});

% Get the constants D in front of the quadharmonic operators:
str = strrep(strL, 'laplacian', '0*');
str = strrep(str, 'lap', '0*');
str = strrep(str, 'biharmonic', '0*');
str = strrep(str, 'biharm', '0*');
str = strrep(str, 'triharmonic', '0*');
str = strrep(str, 'triharm', '0*');
str = strrep(str, 'quadharmonic', '');
str = strrep(str, 'quadharm', '');
str = strrep(str, 'quintharmonic', '0*');
str = strrep(str, 'quintharm', '0*');
func = eval(str);
D = feval(func, inputs{:});

% Get the constants E in front of the quintharmonic operators:
str = strrep(strL, 'laplacian', '0*');
str = strrep(str, 'lap', '0*');
str = strrep(str, 'biharmonic', '0*');
str = strrep(str, 'biharm', '0*');
str = strrep(str, 'triharmonic', '0*');
str = strrep(str, 'triharm', '0*');
str = strrep(str, 'quadharmonic', '');
str = strrep(str, 'quadharm', '');
str = strrep(str, 'quintharmonic', '');
str = strrep(str, 'quintharm', '');
func = eval(str);
E = feval(func, inputs{:});

% Compute L:
for k = 1:nVars
    L = [L; A(k)*lapmat + B(k)*bihmat + C(k)*trihmat + D(k)*quadhmat + ...
        E(k)*quinthmat];
end

%% Disretize the differentiation term of the nonlinear part:

% We only support no differentiation, i.e., NC = 1.
% [TODO]: Support diff order > 1.
Nc = 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function calpak_test79 ( )

%*****************************************************************************80
%
%% CALPAK_TEST79 tests YMDF_TO_WEEKDAY_HEBREW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST79\n' );
  fprintf ( 1, '  For the HEBREW calendar:\n' );
  fprintf ( 1, '  YMDF_TO_WEEKDAY_HEBREW\n' );
  fprintf ( 1, '  returns the day of the week.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED   YMDF           Day of the week\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_hebrew ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0 )
      break
    end

    if ( jed_epoch <= jed1 )

      jed2 = jed_to_next_noon ( jed1 );

      [ y1, m1, d1, f1 ] = jed_to_ymdf_hebrew ( jed2 );
      s1 = ymdf_to_s_hebrew ( y1, m1, d1, f1 );

      w2 = ymdf_to_weekday_hebrew ( y1, m1, d1, f1 );
      s2 = weekday_to_name_hebrew ( w2 );

      fprintf ( 1, '  %11.2f  %20s  %d  %20s\n', jed2, s1, w2, s2 );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_bug1599(datainfo, version)

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_sourceanalysis ft_inverse_lcmv

% fixedori is set correctly in ft_inverse_lcmv, the only problem is
% that ft_sourceanalysis seems to not bother about this when using a
% predefined filter

if nargin<1
  datainfo = ref_datasets;
end
if nargin<2
  version = 'latest';
end

%% Try to reproduce
% load your very favourite data
k = 10; % ctf data
inputfile = fullfile(datainfo(k).origdir,version,'timelock',datainfo(k).type,['timelock_trl_' datainfo(k).datatype '.mat']);
load(inputfile);

% make grid'n'vol
vol = [];
vol.o = [0 0 4];
vol.r = 12;
vol.unit = 'cm';
vol.type = 'singlesphere';

sourcemodel = [];
sourcemodel.resolution = 2.5;

% compute filter
cfg                 = [];
cfg.headmodel       = vol;
cfg.sourcemodel     = sourcemodel;
cfg.method          = 'lcmv';
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.fixedori   = 'no';
filter              = ft_sourceanalysis(cfg, timelock);

cfg               = [];
cfg.headmodel     = vol;
cfg.method        = 'lcmv';
% cfg.rawtrial    = 'yes';
cfg.sourcemodel   = keepfields(filter, {'pos', 'filter', 'filterdimord', 'label'});
cfg.lcmv.fixedori = 'yes';
oriyes            = ft_sourceanalysis(cfg, timelock);
cfg.lcmv.fixedori = 'no';
orino             = ft_sourceanalysis(cfg, timelock);

% cfg is of course different, see above, so remove it
orino             = rmfield(orino, 'cfg');
oriyes            = rmfield(oriyes, 'cfg');

if ~isfield(oriyes.avg, 'ori')
  error('ori not added when it should be');
end

if isfield(orino.avg, 'ori')
  error('ori added when it should not be');
end

fprintf('test_bug1599 is all fine\n');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function f = flops_mul(n,m,c)
% FLOPS_MUL    Flops for real matrix multiplication.
% FLOPS_MUL(a,b) returns the number of flops for a*b.
% FLOPS_MUL(n,m,c) returns the number of flops for ones(n,m)*ones(m,c).
% This is more accurate than the flopcount returned by matlab.

if nargin == 2
  a = n;
  b = m;
  n = rows(a);
  m = cols(a);
  c = cols(b);
end
f = n*c*(2*m-1);

author: andrewssobral

repo: lrslibrary


text:
% GPS Signal Capture Script
% This requires LeCroy's ActiveDSO Interface, available at http://www.lecroy.com
% No additional MathWorks products are required.
% The scope state has been setup via front panel controls:
% 50 MSPS rate, 48 M words record length, 0.1V/div, BWL 20 MHz
% and a single shot capture has been made.
% Copywrite 2002-2011 The Mathworks, Inc.
% Dick Benson
%% Create ActiveX Control to establish communication with scope
hDSO = actxcontrol('lecroy.activedsoctrl.1',[130 80 800 600]);
MakeConnection(hDSO,'TCPIP: 192.168.0.13'); % static IP set for this scope

%  inspect(hDSO)      %What do we know about it?
%  methodsview(hDSO)  %What can we do with it, and how?

% returns unsigned int8, need to remove offset, and recast to int8
y = int8(int16((GetByteWaveform(hDSO,'C1',48e6,0)))-128); % 48 million samples from channel 1
figure; plot(y(1:1e5)); % sanity check
Disconnect(hDSO);
clear hDSO
s=[6 10 30 21 2 29 5]; % satellites that should be in the data set
% save the capture
save GPS_long_3_X y s

author: Sable

repo: mcbench-benchmarks


text:
%MDL_LWR Create model of Kuka LWR manipulator
%
% MDL_LWR is a script that creates the workspace variable KR5 which
% describes the kinematic characteristics of a Kuka KR5 manipulator using
% standard DH conventions.
%
% Also define the workspace vectors:
%   qz        all zero angles
%
% Notes::
% - SI units of metres are used.
%
% Reference::
% - Identifying the Dynamic Model Used by the KUKA LWR: A Reverse Engineering Approach
%   Claudio Gaz Fabrizio Flacco Alessandro De Luca
%   ICRA 2014
%
% See also mdl_kr5, mdl_irb140, mdl_puma560, SerialLink.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% MODEL: Kuka, LWR, 7DOF, standard_DH

d1 =0.4; d2 = 0.39;

%All link lengths and offsets are measured in m
clear L
%            theta    d           a       alpha
links = [
	    Link([0        0           0       pi/2])
		Link([0        0           0      -pi])
		Link([0        d1          0      -pi/2])
		Link([0        0           0       pi/2])
		Link([0        d2          0       pi/2])
		Link([0        0           0       -pi/2])
		Link([0        0           0       0])
	];

LWR=SerialLink(links, 'name', 'Kuka LWR');

qz = [0 0 0 0 0 0 0];

author: petercorke

repo: robotics-toolbox-matlab


text:
function bvec = bvec_next_grlex ( n, bvec )

%*****************************************************************************80
%
%% BVEC_NEXT_GRLEX generates the next binary vector in GRLEX order.
%
%  Discussion:
%
%    N = 3
%
%    Input      Output
%    -----      ------
%    0 0 0  =>  0 0 1
%    0 0 1  =>  0 1 0
%    0 1 0  =>  1 0 0
%    1 0 0  =>  0 1 1
%    0 1 1  =>  1 0 1
%    1 0 1  =>  1 1 0
%    1 1 0  =>  1 1 1
%    1 1 1  =>  0 0 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension.
%
%    Input, integer BVEC(N), the binary vector whose successor is desired.
%
%    Output, integer BVEC(N), the successor to the input vector.
%

%
%  Initialize locations of 0 and 1.
%
  if ( bvec(1) == 0 )
    z = 1;
    o = 0;
  else
    z = 0;
    o = 1;
  end
%
%  Moving from right to left, search for a "1", preceded by a "0".
%
  for i = n : -1 : 2
    if ( bvec(i) == 1 )
      o = i;
      if ( bvec(i-1) == 0 )
        z = i - 1;
        break
      end
    end
  end
%
%  BVEC = 0
%
  if ( o == 0 )
    bvec(n) = 1;
%
%  01 never occurs.  So for sure, B(1) = 1.
%
  elseif ( z == 0 )
    s = sum ( bvec(1:n) );
    if ( s == n )
      bvec(1:n) = 0;
    else
      bvec(1:n-s-1) = 0;
      bvec(n-s:n) = 1;
    end
%
%  Found the rightmost "01" string.
%  Replace it by "10", and shift following 1's to end.
%
  else
    bvec(z) = 1;
    bvec(o) = 0;
    s = sum ( bvec(o+1:n) );
    bvec(o+1:n-s) = 0;
    bvec(n+1-s:n) = 1;
  end

  return
end


text:
function [Contain]=BuildMPS(A, b, Aeq, beq, cost, L, U, PbName, varargin)
%
% function Contain=BuildMPS(A, b, Aeq, beq, cost, L, U, PbName); OR
%          Contain=BuildMPS(..., Param1, Value1, ...);
%
% Build ascii fixed-width MPS matrix string that contains linear
% programming (LP) problem:
%
% Minimizing (for x in R^n): f(x) = cost'*x, subject to
%       A*x <= b        (LE)
%       Aeq*x = beq     (EQ)
%       L <= x <= U     (BD).
%
% Also supported is integer/mixte programming problem similar to the above,
% where a subset of components of x is restricted to be integer (N set) or
% binary set {0,1}.
%
% INPUTS:
%   A: (m x n) matrix
%   b: (m x 1) matrix
%   Aeq: (k x n) matrix
%   beq: (k x 1) matrix
%   cost: (n x 1) matrix
%   L: (1 x n), (n x 1) or (1 x 1)
%   U: (1 x n), (n x 1) or (1 x 1)
%
% Remark: To disable constraint(s) (LE, EQ, BD), please use empty []
%         for corresponding input matrix/rhs parameters.
%
% Optional:
%   - PbName is a string of problem name, default value is 'GENERIC'.
% Other Params:
%    'EleNames', 'EqtNames', or 'VarNames'
%     Cells contain string of respectively
%        (LE) equations, (EQ) equations, or variable names
%   - 'EleNameFun', 'EqtNameFun', or 'VarNameFun'
%     Corresponding Value are function handles that return
%     Equation/Variable name from equation/Variable number
%       Example: > VarNameFun=@(m) char('x'+(m-1));
%     These functions will NOT be used if names of equations/variables
%     are defined.
%   - Param is 'MPSfilename': output MPS file to be saved
%     No saving if MPSfilename is undefined.
%   - 'I', 'Int', 'Integer', 'Integers'
%       Array that stores the indexes that defines the set of integer
%       variables (>=0). The indexes must belong to [1,..., n] and
%       correspond to the column of A, Aeq.
%   - 'B', 'Bin', 'Binary', 'Binaries'
%       Array that stores the index that defines the set of binary
%       variables {0,1}. Indexes follow the same convention as with integer
%       case.
%   - 'QUAD', 'Q': a structure with following fields
%         'Q': (n x n) matrix
%            The lower triangle of Q is assumed to be the transpose of the
%            upper triangle (in other word, Q must be symmetric and we use
%            only the upper part).
%         'g': (n x 1) vector
%         'bquad': scalar
%         'name' (optional): string, name of the constraint.
%                if qs.name is 'COST' then the quadratic term in the
%                cost function to be minimized (see below)
%         'type' (optional): must contains the string 'QLE'
%                (This field is reserved for future for extension of BuildMPS)
%
%      QUAD parameters is used to enforce an additional quadratic
%      constraint on the unknown x of the type:
%           0.5*x'*Q*x + g'*x <= bquad        (QLE)
%
%      Provide as many QUAD parameters as the number of constraints to be
%      meet.
%      Spatial case: if name is 'COST' then it corresponds to a quadratic
%      term of the cost function
%           f(x) = cost'*x + 0.5*x'*Q*x
%      For this case, the 'g' and 'bquad' fields will be ignored.
%
% OUTPUT:
%   Contain: char matrix of the MPS format description of LP/IP problem.
%
% RESTRICTION:
%   Only single column rhs (b and beq) is supported.
%
% The MPS (Mathematical Programming System) file format was introduced by
% IBM in 1970s, but has also been accepted by most subsequent linear
% programming codes. To learn about MPS format, please see:
%   http://lpsolve.sourceforge.net/5.5/mps-format.htm
%
% See also: SaveMPS
%
% Usage example:
%
%     A = [1 1 0; -1 0 -1];
%     b = [5; -10];
%     L = [0; -1; 0];
%     U = [4; +1; +inf];
%     Aeq = [0 -1 1];
%     beq = 7;
%     cost = [1 4 9];
%     VarNameFun = @(m) (char('x'+(m-1))); % returning varname 'x', 'y' 'z'
%
%     Qle = [2 1 0;
%          1 2 0;
%          0 0 1];
%     g = [0; 0; -3];
%     bquad = 100;
%     quad_le = struct('Q', Qle, ...
%                      'g', g, ...
%                      'bquad', bquad);
%
%     Qcost = speye(3);
%     quad_cost = struct('Q', Qcost, ...
%                        'name', 'cost'),
%     Contain = BuildMPS(A, b, Aeq, beq, cost, L, U, 'Pbtest', ...
%                        'VarNameFun', VarNameFun, ...
%                        'EqtNames', {'Equality'}, ...
%                        'Q', quad_le, 'Q', quad_cost, ...
%                        'Integer', [1], ... % first variable 'x' integer
%                        'MPSfilename', 'Pbtest.mps');
%
% Author: Bruno Luong
% update: 15-Jul-2008: sligly improved number formatting
%         25-Aug-2009: Improvement in handling sparse matrix
%         03-Sep-2009: integer/binary variables
%         02-May-2010: quadratic term

if nargin<8 || isempty(PbName)
    PbName='GENERIC';
end

%
% Columns indices of MPS fields
%
idx1=02:03;
idx2=05:12;
idx3=15:22;
idx4=25:36;
idx5=40:47;
idx6=50:61;
idxlist={idx1 idx2 idx3 idx4 idx5 idx6};

%
% Default returned value if error occurs
%
Contain=[]; %#ok
OK = 0;  %#ok

%
% Get the size of the input matrices
%
[neq nvar]=size(Aeq);
[nle sizeA2]=size(A);

if neq==0 % Aeq is empty, i.e., no equality constraint
    nvar=sizeA2;
    Aeq=zeros(0,nvar);
elseif nle==0 % A is empty, i.e., no LE constraint
    sizeA2=nvar;
    A=zeros(0,nvar);
end

%
% Default values for naming functions (nested functions)
%
elenamefun = @elename;
eqtnamefun = @eqtname;
varnamefun = @varname;
MPSfilename = ''; % MPSfilename

% default empty integer and binary sets
iset = [];
bset = [];

%
% Parse options (varargin)
%
parseoptions(varargin{:});

% Number of quadratic constraints
nquadle = length(quadle);

if ~exist('elenames','var')
    elenames=arrayfun(elenamefun, (1:nle), 'UniformOutput', false);
end
if ~exist('eqtnames','var')
    eqtnames=arrayfun(eqtnamefun, (1:neq), 'UniformOutput', false);
end
if ~exist('varnames','var')
    varnames=arrayfun(varnamefun, (1:nvar), 'UniformOutput', false);
end

if nargin<6 || isempty(L)
    L=-inf(1,nvar);
elseif isscalar(L) % extend L if it's a scalar input
    Lval=L;
    L=zeros(1,nvar);
    L(:)=Lval;
else % BUG corrected, reshape L in row
    L = reshape(L,1,[]);
end
if nargin<7 || isempty(U)
    U=+inf(1,nvar);
elseif isscalar(U) % extend U if it's a scalar input
    Uval=U;
    U=zeros(1,nvar);
    U(:)=Uval;
else % BUG corrected, reshape U in row
    U = reshape(U,1,[]);
end

%
% Dimension check
%
if length(beq)~=neq || length(b)~=nle || ...
   length(cost)~=nvar || ...
   length(L)~=nvar || length(U)~=nvar || ...
   sizeA2~=nvar
    error('BuildMPS:DimensionsUnMatched', ...
          'BuildMPS: dimensions do not match');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set problem name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_name=setfields([],0,'NAME');
l_name=setfields(l_name,3,PbName);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set equations in ROWS and COST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_rows=setfields([],0,'ROWS');

l_cost=setfields([],1,'N',2,'COST');

l_rows_eq=emptyline(neq);
for m=1:neq
    l_rows_eq(m,:)=setfields(l_rows_eq(m,:),1,'E',2,eqtnames{m});
end

l_rows_le=emptyline(nle+nquadle);
for m=1:nle
    l_rows_le(m,:)=setfields(l_rows_le(m,:),1,'L',2,elenames{m});
end
for m=1:nquadle
    l_rows_le(nle+m,:)=setfields(l_rows_le(nle+m,:),1,'L',2,quadle(m).name);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set coefficients of constraint equations in COLUMNS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put all the linear constraint terms of quadratic in a matrix
quadle_g = cat(2,quadle(:).g).'; % (nquadle x nvar)
CostAeq = [cost(:).';
           Aeq;
           A;
           quadle_g]; % CostAeq is sparse if any is sparse
MustWrite = (CostAeq ~= 0);
NWrite = sum(MustWrite,1);
NLines = sum(ceil(NWrite/2));

l_columns=setfields([],0,'COLUMNS');
l_columnsbody=emptyline(NLines);

c=0;
for n=1:nvar % Loop over variables
    var=varnames{n};
    field=3;
    eqtn = find(MustWrite(:,n)); % subset of (1:1+neq+nle+nquadle)
    for m=eqtn(:).' % 1:1+neq+nle+nquadle % Loop over eqt
        if m==1
            colname='COST';
            val = cost(n);
        elseif m<=1+neq
            colname=eqtnames{m-1};
            val=Aeq(m-1,n);
        elseif m<=1+neq+nle
            colname=elenames{m-(1+neq)};
            val=A(m-(1+neq),n);
        else
            iquad = m-(1+neq+nle);
            colname=quadle(iquad).name;
            val=quadle_g(iquad,n);
        end
        if field==3
            c=c+1;
            l_columnsbody(c,:)=setfields(l_columnsbody(c,:),...
                2,var,...
                field,colname, ...
                field+1,val);
            field=5;
        else % field==5
            l_columnsbody(c,:)=setfields(l_columnsbody(c,:),...
                field,colname, ...
                field+1,val);
            field=3;
        end
    end % for-loop eqt
end % for-loop variables
l_columnsbody(c+1:end,:)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set equation RHS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
quadle_bquad = [quadle(:).bquad]; % (1 x nquadle)

rhs=[beq(:); b(:); quadle_bquad(:)];
MustWrite = (rhs ~= 0);
NWrite = sum(MustWrite);
NLines = ceil(NWrite/2);

l_rhs=setfields([],0,'RHS');
l_rhsbody=emptyline(NLines);
c=0;
field=3;
eqt = find(MustWrite); % subset of (1:neq+nle)
for m=eqt(:).' % 1:neq+nle+nquadle % Loop over eqt
    if m<=neq
        colname=eqtnames{m};
        val=rhs(m);
    elseif m<=neq+nle
        colname=elenames{m-neq};
        val=rhs(m);
    else
        colname=quadle(m-(neq+nle)).name;
        val=rhs(m);
    end
    if field==3
        c=c+1;
        l_rhsbody(c,:)=setfields(l_rhsbody(c,:),...
            2,'RHS',...
            field,colname, ...
            field+1,val);
        field=5;
    else
        l_rhsbody(c,:)=setfields(l_rhsbody(c,:),...
            field,colname, ...
            field+1,val);
        field=3;
    end
end % for-loop eqt
l_rhsbody(c+1:end,:)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set bound constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_bound=setfields([],0,'BOUNDS');

VarType=zeros(size(U));
%
% Var types (local definition)
%
VarType(:)=0; % real
VarType(iset) = 1; % integer
VarType(bset) = 2; % binary

% Force lower/upper bound for integer variables to be integer as well
L(iset) = max(ceil(L(iset)),0); % integer lower bound cannot be negative
U(iset) = floor(U(iset));

% Values not used, but we set for clarity
L(bset) = 0;
U(bset) = 1;

upinf=(U==inf);
loinf=(L==-inf);
lonz=(L~=0) & ~loinf;

BoundType=zeros(size(U));

%
% Bound types (local definition)
%
BoundType(:) = 3; % Default, 0<=x, real variable
BoundType(upinf & loinf) = 1; % free, real
BoundType(upinf & lonz) = 2; % lo<=x (lo ~= 0), real
BoundType(~upinf & lonz) = 4; % lo<=x<=up, integer or real
BoundType(~upinf & loinf) = 5; % x<=up, real
BoundType(~upinf & ~loinf & ~lonz) = 6; % 0<=x<=up, integer or real
BoundType(upinf & VarType==1) = 7; %  lo<=x, integer
BoundType(bset) = 8; % binary, x = 0 or 1

NLines = sum(ismember(BoundType,[1 2 6 7 8])) + ...
         sum(ismember(BoundType,[4 5]))*2;
l_boundbody=emptyline(NLines);
c=0;
for n=1:nvar
    var=varnames{n};
    lo=L(n);
    up=U(n);
    vtype = VarType(n);
    if (vtype==2) % Type 8, binary variables
        c=c+1;
        l_boundbody(c,:)=setfields(l_boundbody(c,:),...
            1, 'BV', ...
            2, 'BND1', ...
            3, var, ...
            4, 1); % Field 4 must be 1.0 or blank
    elseif (up==inf)
        if (lo==-inf) % Type 1, Free real variable, one line
            c=c+1;
            l_boundbody(c,:)=setfields(l_boundbody(c,:),...
                1, 'FR', ...
                2, 'BND1', ...
                3, var, ...
                4, 0);
        elseif (lo~=0) || (vtype==1) % Type 2, or Type 7 lo<=x, one line
            c = c+1;
            if vtype==1 % integer, Type 7
                LOstr = 'LI';
            else % real, Type 2
                LOstr = 'LO';
            end
            l_boundbody(c,:)=setfields(l_boundbody(c,:),...
                1, LOstr, ...
                2, 'BND1', ...
                3, var, ...
                4, lo);
        % else 0<=x<=inf: Type3, real variable nothing to write
        end
    else % up<inf
        if lo>-inf
            if lo~=0 % Type 4, lo<=x<=up
                c=c+1;
                if vtype==1 % integer
                    LOstr = 'LI';
                else % real
                    LOstr = 'LO';
                end
                l_boundbody(c,:)=setfields(l_boundbody(c,:),...
                    1, LOstr, ...
                    2, 'BND1', ...
                    3, var, ...
                    4, lo);
            %else % 0<=x<=up % Type 6
            end
        else % if lo==-inf % Type 5, x<=up
            c=c+1;
            l_boundbody(c,:)=setfields(l_boundbody(c,:),...
                1, 'MI', ...
                2, 'BND1', ...
                3, var, ...
                4, 0);
        end
        % Common Type 4, 5, or 6
        % Type 6 is 0<=x<=up
        c=c+1;
        if vtype==1 % integer
            HIstr = 'UI';
        else % real
            HIstr = 'UP';
        end
        l_boundbody(c,:)=setfields(l_boundbody(c,:),...
            1, HIstr, ...
            2, 'BND1', ...
            3, var, ...
            4, up);
    end
end % for-loop on variable
l_boundbody(c+1:end,:)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quad section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for m=1:nquadle
    % Get the current quad structure
    qs = quadle(m);
    l_quad = setfields([],0,'QSECTION');
    l_quad = setfields(l_quad,3,qs.name);

    % Use only the upper-triangular part of Q
    [i j Qij] = find(triu(qs.Q));
    NLines = length(Qij);
    l_quadbody = emptyline(NLines);

    for n=1:NLines % Loop over non-zeros elements
        vari=varnames{i(n)};
        varj=varnames{j(n)};
        l_quadbody(n,:) = setfields(l_quadbody(n,:), ...
            2,vari,...
            3,varj, ...
            4,Qij(n));
    end

    quadle(m).qsection = [l_quad;
                          l_quadbody]; %#ok

end % for-loop on quadratic constraints

if nquadle>1
    % concatenate together all the qsections
    l_allquad = cat(1,quadle(:).qsection);
else
    % empty line
    l_allquad = l_name([],:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the last card
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_end=setfields([],0,'ENDATA');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Concatenate together all parts of mps format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Contain=[l_name; ...
         l_rows; ...
         l_cost; ...
         l_rows_eq; ...
         l_rows_le; ...
         l_columns; ...
         l_columnsbody; ...
         l_rhs; ...
         l_rhsbody; ...
         l_bound; ...
         l_boundbody; ...
         l_allquad; ...
         l_end];

if ~isempty(MPSfilename)
    %
    % Save the Contain in MPSfilename
    %
    OK = SaveMPS(MPSfilename, Contain);
    if ~OK % Something is wrong during saving
        warning('BuildMPS:SavingFailure', ...
                ['BuildMPS: Cannot save ' MPSfilename]);
    end
else % Nothing to save
    OK = 1;
end

% return % Uncomment the RETURN statement causes M-lint to crash on 2009A
% There is no instructions from now on, juts nested functions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nested functions: BE AWARE, the functions have access to local
% variables of BuildMPS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Generate n empty lines of MPS data
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function l=emptyline(n)
        if nargin<1 || isempty(n)
            n=1;
        end
        l=char(zeros(n,61));
        l(:)=' ';
    end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Convert to string at the fixed length of 12
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function str=num2fixedlengthstr(num, maxlength, roundingflag)
        % function str=num2fixedlengthstr(num); OR
        % str=num2fixedlengthstr(..., maxlength, roundingflag);
        %
        % Convert double NUM to decimal string having MAXLENGTH [12] as maximum
        % length. Smart conversion with accurate result despite length constraint.
        %
        % ROUNDINGFLAG: 0 or [1]
        %   0: truncate fracional part (quicker)
        %   1: rounding fracional part (more accurate).
        %
        % Last update: 15/Aug/2008, remove leading "0" when the string starts
        %              as "0.xxxx"
        %
        if nargin<2
            maxlength=12;
        end

        if nargin<3
            roundingflag=1; % rounding by default
        end

        if num>=0
            fracNDigits=maxlength;
        else
            fracNDigits=maxlength-1;
        end
        % "%G" format:
        % ANSI specification X3.159-1989: "Programming Language C,"
        % ANSI, 1430 Broadway, New York, NY 10018.
        str=num2str(num,['%0.' num2str(fracNDigits) 'G']);
        %
        % Try to compact the string data to fit inside the field length
        %
        while length(str)>maxlength
            if regexp(str,'^0\.') % delete the leading 0 in "0.xxx"
                str(1)=[];
                continue;
            end
            [istart iend]=regexp(str,'[+-](0)+'); % +/- followed by multiples 0
            if ~isempty(istart) % Remove zero in xxxE+000yy or xxxE-000yy
                str(istart+1:iend)=[];
                continue
            else
                [istart iend]=regexp(str,'E[+]');
                if ~isempty(istart) % Remove "+" char in xxxE+yyy
                    str(iend)=[];
                    continue
                end
            end
            idot=find(str=='.',1,'first');
            if ~isempty(idot)
                iE=find(str=='E',1,'first');
                if roundingflag % rounding fraction part
                    % Calculate the Length of the fractional part
                    % Adjust its number of digits and start over again
                    if ~isempty(iE) % before the mantissa
                        fracNDigits=maxlength-length(str)+iE-idot-1;
                        str=num2str(num,['%0.' num2str(fracNDigits) 'E']);
                    else %if idot<=maxlength+1 % no manissa
                        fracNDigits=maxlength-idot;
                        str=num2str(num,['%0.' num2str(fracNDigits) 'f']);
                    end
                    roundingflag=0; % won't do rounding again
                    continue % second pass with new string
                else
                    % truncate the fractional part
                    if ~isempty(iE) % before the mantissa
                        str(maxlength-length(str)+iE:iE-1)=[];
                        return;
                    else %if idot<=maxlength+1 % no mantissa
                        str(maxlength+1:end)=[];
                        return;
                    end
                end
            end
            % it should not never go here, unless BUG
            error('BuildMPS: cannot convert %0.12e to string\n',num);
        end % while loop

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Set the field of an MPS line by value
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function l=setfield(l,field,var)

        if isnumeric(var) % numerical data, convert to string
            var=num2fixedlengthstr(var); % convert to 12-length string
        end

        if isempty(l)
            l=emptyline;
        end
        if ~isempty(field) && field>0
            idx=idxlist{field};
        else
            idx=1:61;
        end
        if length(var)>length(idx)
            var=var(1:length(idx));
        else
            idx=idx(1:length(var));
        end
        l(idx)=var;

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Set multiple fields of an MPS line by values
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function l=setfields(l, varargin)
        for k=1:2:length(varargin)
            l=setfield(l, varargin{k}, varargin{k+1});
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate equation name for (LE) constraint
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function name=elename(m)
        name=['LE' num2str(m)];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate equation name for (EQ) constraint
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function name=eqtname(m)
        name=['EQ' num2str(m)];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate variable name
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function name=varname(n)
        name=['X' num2str(n)];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate equation name for (QLE) constraint
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function name=quadlename(m)
        name=['QLE' num2str(m)];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Parse a pair of Name/Value option
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function parseoption(strname, value)
        if ischar(strname)
            strname = strtrim(lower(strname));
            switch strname
                case 'elenames',
                    if ~iscell(value) || length(value)~=nle || ...
                        ~all(cellfun(@ischar, value))
                        error('BuildMPS:IncorrectEleNames', ...
                    'BuildMPS: EleNames must be cell of %d strings', nle);
                    end
                    elenames = value;
                case 'eqtnames',
                    if ~iscell(value) || length(value)~=neq || ...
                            ~all(cellfun(@ischar, value))
                        error('BuildMPS:IncorrectEqtNames', ...
                    'BuildMPS: EqtNames must be cell of %d strings', neq);
                    end
                    eqtnames = value;
                case 'varnames',
                    if ~iscell(value) || length(value)~=nvar || ...
                            ~all(cellfun(@ischar, value))
                        error('BuildMPS:IncorrectVarNames', ...
                    'BuildMPS: VarNames must be cell of %d strings', nvar);
                    end
                    varnames = value;
                case 'varnamefun',
                    if ischar(value)
                        value=str2func(value);
                    end
                    if ~isa(value,'function_handle')
                        error('BuildMPS:IncorrectVarNameFun', ...
                              'BuildMPS: VarNameFun must be a function');
                    end
                    varnamefun = value;
                case 'eqtnamefun',
                    if ischar(value)
                        value=str2func(value);
                    end
                    if ~isa(value,'function_handle')
                        error('BuildMPS:IncorrectEqtNameFun', ...
                              'BuildMPS: EqtNameFun must be a function');
                    end
                    eqtnamefun = value;
                case 'elenamefun',
                    if ischar(value)
                        value=str2func(value);
                    end
                    if ~isa(value,'function_handle')
                        error('BuildMPS:IncorrectEleNameFun', ...
                              'BuildMPS: EleNameFun must be a function');
                    end
                    elenamefun = value;
                case 'mpsfilename',
                    if ~ischar(value)
                        error('BuildMPS:IncorrectMPSfilename', ...
                              'BuildMPS: MPSfilename must be a string');
                    end
                    MPSfilename = value;
                case  {'i' 'int' 'integer' 'integers'},
                    iset = value(:);
                    if any(iset<1 | iset>nvar)
                        error('Integer set contains invalid index');
                    end
                case {'b' 'bin' 'binary' 'binaries'},
                    bset = value(:);
                    if any(bset<1 | bset>nvar)
                        error('Binary set contains invalid index');
                    end
                case {'quad' 'q'},
                    qcounter = length(quadle)+1;
                    % Basic check of quad structure
                    if isstruct(value)
                        qs = value;
                        if ~isfield(qs,'Q') || ~isequal(size(qs.Q),[nvar nvar])
                            error('Missing or invalid <Q> field in QUAD structure');
                        end
                        if ~isfield(qs,'g') || isempty(qs.g)
                            qs.g = zeros(nvar,1);
                        elseif isequal(size(qs.g), [1 nvar])
                            % reshape in column
                            qs.g = qs.g(:);
                        elseif ~isequal(size(qs.g), [nvar 1])
                            error('Invalid <g> field in QUAD');
                        end
                        if ~isfield(qs,'bquad') || isempty(qs.bquad)
                            qs.bquad = 0;
                        end
                        if ~isscalar(qs.bquad)
                            error('Missing or invalid <bquad> field in QUAD structure');
                        end
                        if ~isfield(qs,'type')
                            qs.type = 'QLE';
                        end
                        if ~strcmpi(qs.type ,'QLE')
                            error('Invalid <type> field in QUAD structure');
                        end
                        if ~isfield(qs,'name') || isempty(qs.name)
                            qs.name = quadlename(qcounter);
                        elseif strcmpi(qs.name,'COST') %
                            qs.name = 'COST'; % force to be upper case
                            % The linear term for the functional must be
                            % provides in the 5th parameter
                            % so we set 'g' to zero
                            qs.g(:) = 0;
                            qs.bquad(:) = 0;
                        end
                    else
                        if isempty(value)
                            return % ignore empty argument
                        end
                        error('Invalid input QUAD (must be a structure)');
                    end
                    quadle(qcounter) = orderfields(qs);
                otherwise
                    warning('BuildMPS:UnknownParams', ...
                        ['BuildMPS: Unknown parameter ' strname]);
            end
        else
            error('BuildMPS:IncorrectCall', ...
                  'BuildMPS: options must be pair of Name/Value');
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Parse options
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function parseoptions(varargin)
        if mod(nargin,2)
            error('BuildMPS:IncorrectCall', ...
                  'BuildMPS: options must be pair of Name/Value');
        end

        % default empty quadle
        quadle = struct('Q', {}, ...
                        'g', {}, ...
                        'bquad', {}, ...
                        'type', {}, ...
                        'name', {} ...
                        );

        quadle = orderfields(quadle);

        %
        % Loop over pair of Name/Value option
        %
        for ivararg=1:2:nargin
            parseoption(varargin{ivararg},varargin{ivararg+1});
        end

    end

end % BuildMPS

author: opencobra

repo: cobratoolbox


text:
function plot_srcColorbar(cmin,cmax)

colormap(abs((gui_Colormap_ReadPovRay('SRC-Discrete.pov'))))
set(gca,'CLim',[cmin cmax]);
colorbar('Location','EastOutside',...
    'YLim',[cmin cmax]);

author: CelsoReyes


text:
function cout = int(cin, delta)

%tstoolbox/@core/int
%   Syntax:
%     * cout = int(cin, delta)
%
%   Input Arguments:
%     * cin - core object
%     * delta - time period between two data samples
%
%   numerical integration along dimension 1 when data was sampled
%   equidistantly with samplerate = 1/delta
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

tmp = cumsum(data(cin),1)*delta;
tmp = [zeros(1, dlens(cin,2)) ; tmp];
cout = core(tmp);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function  [ConSPix ConEdge]= find_connect_superpixel(labels, K, height ,width )
%%
% obtain the neighbour relationship of the super-pixels
% Input:
%         labels:    the super-pixel label obtained from SLIC
%         K:         the number of super-pixels
%         height:    the height of the image
%         width:     the width of the image
% Output:
%         ConPix:   the one layer neighbour relationship
%%%%=====================================================
ConSPix=zeros(K,K);
%the one outerboundary super
for i=1:height-1
    for j=1:width-1
        if labels(i,j)~=labels(i,j+1)
            ConSPix(labels(i,j) ,labels(i,j+1) )=1;
        end
        if labels(i,j)~=labels(i+1,j)
            ConSPix(labels(i,j) ,labels(i+1,j) )=1;
        end
    end
    if labels(i,j+1)~=labels(i+1,j+1)
        ConSPix(labels(i,j+1) ,labels(i+1,j+1 ) )=1;
    end
end
for j=1:width-1
    if labels(height,j)~=labels(height,j+1)
        ConSPix( labels(height,j),labels(height,j+1) )=1;
    end
end
for i=1:height-1
    for j=1:width-1
        if labels(i,j)~=labels(i+1,j+1)
            ConSPix( labels(i,j),labels(i+1,j+1) )=1;
        end
    end
end
for i=1:height-1
    for j=2:width
        if labels(i,j)~=labels(i+1,j-1)
            ConSPix( labels(i,j),labels(i+1,j-1) )=1;
        end
    end
end
 ConSPix = ConSPix + ConSPix';
 ConSPix(ConSPix>0)=1;
 [edges_x edges_y] = find(triu(ConSPix)==1);
 ConEdge = [edges_x edges_y];
 %[x y] = meshgrid(1:K,1:K);
 %ConEdge = [x(:) y(:)];

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function F = airy(K, F, scale, pref)
%AIRY   Airy function of a CHEBFUN.
%   AIRY(F) returns the Airy function Ai(F) of a CHEBFUN F.
%
%   AIRY(K, F) returns various Airy functions specified by K:
%     0 - (default) is the same as airy(Z)
%     1 - returns the derivative, Ai'(Z)
%     2 - returns the Airy function of the second kind, Bi(Z)
%     3 - returns the derivative, Bi'(Z)
%
%   AIRY(K, F, SCALE) returns a scaled AIRY(K, F) specified by SCALE:
%     0 - (default) is that same as AIRY(K, Z)
%     1 - returns airy(K, F) scaled by EXP(2/3*F.^(3/2)) for K = 0, 1,
%         and scaled by EXP(-ABS(2/3.*REAL(F.^(3/2)))) for K = 2, 3.
%
% See also BESSELJ.

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

% Parse the inputs as described in help text:
if ( nargin == 1 )
    F = K;
    K = 0;
    pref = chebfunpref();
end
if ( nargin < 3 )
    scale = 0;
    pref = chebfunpref();
end
if ( nargin == 3 )
    if ( isa(scale, 'chebfunpref') )
        pref = scale;
        scale = 0;
    else
        pref = chebfunpref;
    end
end

for k = 1:numel(F)
    F(k) = columnAiry(K, F(k), scale, pref);
end

end

function g = columnAiry(K, f, scale, pref)

% The standard Airy function:
g = compose(f, @(x) airy(K, x), [], pref);

if ( scale == 0 )
    % Standard case (no scaling).

elseif ( (scale == 1) && ((K == 0) || (K == 1)) )
    % Scaled with k = 1, 2:
    scl = exp(2/3*f.^(3/2));
    g = scl.*g;

elseif ( (scale == 1) && ((K == 2) || (K == 3)) )
    % Scaled with k = 2, 3:
    scl = exp(-abs(2/3*real(f.^(3/2))));
    g = scl.*g;

else
    % Invalid parameter sequence:
    error('CHEBFUN:CHEBFUN:airy:params', 'Invalid paramter selection.');

end

end


text:
function [Zc mu] = myCenter(Z)
%--------------------------------------------------------------------------
% Syntax:       Zc = myCenter(Z);
%               [Zc mu] = myCenter(Z);
%
% Inputs:       Z is an (d x n) matrix containing n samples of a
%               d-dimensional random vector
%
% Outputs:      Zc is the centered version of Z
%
%               mu is the (d x 1) sample mean of Z
%
% Description:  This function returns the centered (i.e., zero mean)
%               version of the input samples
%
%               NOTE: Z = Zc + repmat(mu,1,n);
%
% Author:       Brian Moore
%               brimoor@umich.edu
%
% Date:         April 26, 2015
%--------------------------------------------------------------------------

% Compute sample mean
mu = mean(Z,2);

% Subtract mean
Zc = bsxfun(@minus,Z,mu);


text:
%--- help for trnd ---
%
% TRND   Random arrays from Student's t distribution.
%    R = TRND(V) returns an array of random numbers chosen from Student's t
%    distribution with V degrees of freedom.  The size of R is the size of
%    V.
%
%    R = TRND(V,M,N,...) or R = TRND(V,[M,N,...]) returns an M-by-N-by-...
%    array.
%
%    The t distribution with one degree of freedom is also known as the
%    Cauchy distribution.
%
%    See also TCDF, TINV, TPDF, TSTAT, NCTRND, RANDOM.
%
%    Reference page in Doc Center
%       doc trnd
%
%

author: jmaih


text:
function []=ps_weed(all_da_flag,no_weed_adjacent,no_weed_noisy)
%PS_WEED weeds out neighboring PS and save those kept to new version
%   PS_weed(all_da_flag,no_weed_adjacent,no_weed_noisy)
%
%   Andy Hooper, June 2006
%
%   ===================================================================
%   09/2006 AH: create all workspace files directly
%   09/2006 AH: drop noisy pixels
%   09/2006 AH: add small baselines
%   01/2007 AH: drop pixels with duplicate lon/lat
%   05/2007 AH: optional weeding of pixels with zero elevation added
%   03/2009 AH: delete scla_smooth mat files
%   02/2010 AH: change smoothing of arcs to time domain
%   02/2010 AH: option to threshold on max arc noise in any ifg
%   02/2010 AH: Speed up of noise weeding
%   02/2010 AH: Leave out ifgs in drop_ifg_index from noise calculation
%   01/2011 AH: Correct arc DEM error estimation for small baselines
%   06/2011 AH: Add weed_neighbours parm (default still 'y')
%   12/2012 AH: weed very small heights if weed_zero_elevation='y'
%   09/2015 DB: Clean the command line output, store number of PS left.
%   09/2015 DB: Fix bug when no PS is left after weeding 0 elevation.
%   09/2015 AH: use matlab triangulation if triangle program not installed
%   01/2016 DB: Replace save with stamps_save which checks for var size when
%               saving
%   10/2017 DB: if inc is present also do weeding on it.
%   ===================================================================
logit;
logit('Weeding selected pixels...')

if nargin<1
    all_da_flag=0;
end

%weed_alpha=getparm('weed_alpha',1);
time_win=getparm('weed_time_win',1);
weed_standard_dev=getparm('weed_standard_dev',1);
weed_max_noise=getparm('weed_max_noise',1);
weed_zero_elevation=getparm('weed_zero_elevation',1);
weed_neighbours=getparm('weed_neighbours',1);
drop_ifg_index=getparm('drop_ifg_index');
small_baseline_flag=getparm('small_baseline_flag',1);

if nargin<2
    if strcmpi(weed_neighbours,'y')
        no_weed_adjacent=0;
    else
        no_weed_adjacent=1;
    end
end
if nargin<3
    if weed_standard_dev>=pi && weed_max_noise>=pi
        no_weed_noisy=1;
    else
        no_weed_noisy=0;
    end
end

load psver

psname=['ps',num2str(psver)];
pmname=['pm',num2str(psver)];
phname=['ph',num2str(psver)];
selectname=['select',num2str(psver)];
hgtname=['hgt',num2str(psver),'.mat'];
laname=['la',num2str(psver),'.mat'];
incname=['inc',num2str(psver),'.mat'];
bpname=['bp',num2str(psver),'.mat'];
psothername=['ps_other'];
psothername=['pm_other'];
selectothername=['select_other'];
hgtothername=['hgt_other'];
laothername=['la_other'];
incothername=['inc_other'];
bpothername=['bp_other'];

ps=load(psname);
ifg_index=setdiff([1:ps.n_ifg],drop_ifg_index);

sl=load(selectname);

if exist([phname,'.mat'],'file')
    phin=load(phname);
    ph=phin.ph;
    clear phin
else
    ph=ps.ph;
end

day=ps.day;
bperp=ps.bperp;
master_day=ps.master_day;

if isfield(sl,'keep_ix')
    ix2=sl.ix(sl.keep_ix);
    K_ps2=sl.K_ps2(sl.keep_ix);
    C_ps2=sl.C_ps2(sl.keep_ix);
    coh_ps2=sl.coh_ps2(sl.keep_ix);
else
    ix2=sl.ix2;
    K_ps2=sl.K_ps2;
    C_ps2=sl.C_ps2;
    coh_ps2=sl.coh_ps2;
end

ij2=ps.ij(ix2,:);
xy2=ps.xy(ix2,:);
ph2=ph(ix2,:);
lonlat2=ps.lonlat(ix2,:);

pm=load(pmname);
ph_patch2=pm.ph_patch(ix2,:); % use original patch phase, with PS left in
if isfield(sl,'ph_res2')
    ph_res2=sl.ph_res2(sl.keep_ix,:);
else
    ph_res2=[];
end
clear pm

clear sl

clear ph
if isfield(ps,'ph')
    ps=rmfield(ps,'ph');
end
ps=rmfield(ps,{'xy','ij','lonlat','sort_ix'});

if all_da_flag~=0
    pso=load(psothername);
    slo=load(selectothername);
    ix_other=slo.ix_other;
    n_ps_other=sum(ix_other);
    K_ps_other2=pso.K_ps_other(ix_other);
    C_ps_other2=pso.C_ps_other(ix_other);
    coh_ps_other2=pso.coh_ps_other(ix_other);
    ph_res_other2=pso.ph_res_other(ix_other,:);
    ij2=[ij2;pso.ij_other(ix_other,:)];
    xy2=[xy2;pso.xy_other(ix_other,:)];
    ph2=[ph2;pso.ph_other(ix_other,:)];
    lonlat2=[lonlat2;pso.lonlat_other(ix_other,:)];
    clear pso slo

    pmo=load(pmothername);
    ph_patch_other2=pmo.ph_patch_other(ix_other,:);
    clear pm

    K_ps2=[K_ps2;K_ps_other2];
    C_ps2=[C_ps2;C_ps_other2];
    coh_ps2=[coh_ps2;coh_ps_other2];
    ph_patch2=[ph_patch2;ph_patch_other2];
    ph_res2=[ph_res2;ph_res_other2];
else
    n_ps_other=0;
end

if exist(hgtname,'file')
    ht=load(hgtname);
    hgt=ht.hgt(ix2);
    clear ht
    if all_da_flag~=0
        hto=load(hgtothername);
        hgt=[hgt;hto.hgt_other(ix_other)];
        clear hto
    end
end

n_ps_low_D_A=length(ix2);
n_ps=n_ps_low_D_A + n_ps_other;
ix_weed=logical(ones(n_ps,1));
logit([num2str(n_ps_low_D_A),' low D_A PS, ',num2str(n_ps_other),' high D_A PS']);

if no_weed_adjacent==0
	step_name='INITIALISE NEIGHBOUR MATRIX';
    fprintf([step_name '\n'])

    ij_shift=ij2(:,2:3)+repmat([2,2]-min(ij2(:,2:3)),n_ps,1);
	neigh_ix=zeros(max(ij_shift(:,1))+1,max(ij_shift(:,2))+1);
    miss_middle=logical(ones(3));
    miss_middle(2,2)=0;

	for i=1:n_ps
        neigh_this=neigh_ix(ij_shift(i,1)-1:ij_shift(i,1)+1,ij_shift(i,2)-1:ij_shift(i,2)+1);
        neigh_this(neigh_this==0&miss_middle)=i;
        neigh_ix(ij_shift(i,1)-1:ij_shift(i,1)+1,ij_shift(i,2)-1:ij_shift(i,2)+1)=neigh_this;

        if i/100000==floor(i/100000)
            logit([num2str(i),' PS processed'],2)
            save log step_name i
            !sync
        end
    end

	step_name='FIND NEIGHBOURS';
    fprintf([step_name '\n'])

    neigh_ps=cell(n_ps,1);
	for i=1:n_ps
        my_neigh_ix=neigh_ix(ij_shift(i,1),ij_shift(i,2));
        if my_neigh_ix~=0
            neigh_ps{my_neigh_ix}=[neigh_ps{my_neigh_ix},i];
        end
        if i/100000==floor(i/100000)
            logit([num2str(i),' PS processed'],2)
            save log step_name i
            !sync
        end

    end

    clear neigh_ix

	step_name='SELECT BEST';
	fprintf([step_name '\n'])

	for i=1:n_ps
        if ~isempty(neigh_ps{i})
            same_ps=i;
            i2=1;
            while i2 <= length(same_ps)
                ps_i=same_ps(i2);
                same_ps=[same_ps,neigh_ps{ps_i}];
                neigh_ps{ps_i}=[];
                i2=i2+1;
            end
            same_ps=unique(same_ps);
            [dummy,high_coh]=max(coh_ps2(same_ps));
            low_coh_ix=logical(ones(size(same_ps)));
            low_coh_ix(high_coh)=0;
            ix_weed(same_ps(low_coh_ix))=0;
        end
        if i/100000==floor(i/100000)
            logit([num2str(i),' PS processed'],2)
            save log step_name i
            !sync
        end

	end
    logit([num2str(sum(ix_weed)),' PS kept after dropping adjacent pixels']);

end

% output how many PS are left after weeding zero elevations out
if strcmpi(weed_zero_elevation,'y') & exist('hgt','var')
    sea_ix=hgt<1e-6;
    ix_weed(sea_ix)=false;
    logit([num2str(sum(ix_weed)),' PS kept after weeding zero elevation']);
end
xy_weed=xy2(ix_weed,:);
% update PS inofmration
n_ps=sum(ix_weed);

%% Remove dupplicated points
% Some non-adjacent pixels are allocated the same lon/lat by DORIS.
% If duplicates occur, the pixel with the highest coherence is kept.
ix_weed_num=find(ix_weed);
[dummy,I]=unique(xy_weed(:,2:3),'rows');
dups=setxor(I,[1:sum(ix_weed)]'); % pixels with duplicate lon/lat

for i=1:length(dups)
    dups_ix_weed=find(xy_weed(:,2)==xy_weed(dups(i),2)&xy_weed(:,3)==xy_weed(dups(i),3));
    dups_ix=ix_weed_num(dups_ix_weed);
    [dummy,I]=max(coh_ps2(dups_ix));
    ix_weed(dups_ix([1:end]~=I))=0; % drop dups with lowest coh
end

if ~isempty(dups)
    xy_weed=xy2(ix_weed,:);
    logit(sprintf('%d PS with duplicate lon/lat dropped\n',length(dups)'))
end

% update PS inofmration
n_ps=sum(ix_weed);
ix_weed2=true(n_ps,1);

%% Weedign noisy pixels
ps_std=zeros(n_ps,1);
ps_max=zeros(n_ps,1);

% include a check to see if any points are actually left.
if n_ps~=0
    if no_weed_noisy==0

        step_name='DROP NOISY';
        arch=computer('arch');
        if strcmpi(arch(1:3),'win')
            use_triangle='n';
        else
            tripath=system('which triangle >& /dev/null');
            if tripath==0
                use_triangle='y';
            else
                use_triangle='n';
            end
        end

        if use_triangle=='y'
            nodename=['psweed.1.node'];
            fid=fopen(nodename,'w');
            fprintf(fid,'%d 2 0 0\n',n_ps);

            for i=1:n_ps
                fprintf(fid,'%d %f %f\n',i,xy_weed(i,2),xy_weed(i,3));
            end

            fclose(fid);

            system('triangle -e psweed.1.node > triangle_weed.log');

            fid=fopen('psweed.2.edge','r');
            header=str2num(fgetl(fid));
            N=header(1);
            edgs=zeros(N,4);
            for i=1:N
                edgs(i,:)=str2num(fgetl(fid));
            end
            fclose(fid);
            edgs=edgs(:,2:3);
        else
            xy_weed=double(xy_weed);
            tri=delaunay(xy_weed(:,2),xy_weed(:,3));
            tr=triangulation(tri,xy_weed(:,2),xy_weed(:,3));
            edgs=edges(tr);
        end
        n_edge=size(edgs,1);

        ph_weed=ph2(ix_weed,:).*exp(-j*(K_ps2(ix_weed)*bperp'));  % subtract range error
        ph_weed=ph_weed./abs(ph_weed);
        if ~strcmpi(small_baseline_flag,'y')
            ph_weed(:,ps.master_ix)=exp(j*(C_ps2(ix_weed)));  % add master noise
        end
        edge_std=zeros(n_edge,1);
        edge_max=zeros(n_edge,1);
        dph_space=(ph_weed(edgs(:,2),:).*conj(ph_weed(edgs(:,1),:)));
        dph_space=dph_space(:,ifg_index);
        n_use=length(ifg_index);
        for i=1:length(drop_ifg_index)
            if strcmpi(small_baseline_flag,'y')
                logit(sprintf('%s-%s dropped from noise estimation',datestr(ps.ifgday(drop_ifg_index(i),2)),datestr(ps.ifgday(drop_ifg_index(i),2))));
            else
                logit(sprintf('%s dropped from noise estimation',datestr(day(drop_ifg_index(i)))));
            end
        end

        if ~strcmpi(small_baseline_flag,'y')
          logit(sprintf('Estimating noise for all arcs...'))
          dph_smooth=zeros(n_edge,n_use,'single');
          dph_smooth2=zeros(n_edge,n_use,'single');
          for i1=1:n_use
            time_diff=(day(ifg_index(i1))-day(ifg_index))';
            weight_factor=exp(-(time_diff.^2)/2/time_win^2);
            weight_factor=weight_factor/sum(weight_factor);

            dph_mean=sum(dph_space.*repmat(weight_factor,n_edge,1),2);
            dph_mean_adj=angle(dph_space.*repmat(conj(dph_mean),1,n_use)); % subtract weighted mean
            G=[ones(n_use,1),time_diff'];
            m=lscov(double(G),double(dph_mean_adj)',weight_factor); % weighted least-sq to find best-fit local line
            dph_mean_adj=angle(exp(j*(dph_mean_adj-(G*m)'))); % subtract first estimate
            m2=lscov(double(G),double(dph_mean_adj)',weight_factor); % weighted least-sq to find best-fit local line
            dph_smooth(:,i1)=dph_mean.*exp(j*(m(1,:)'+m2(1,:)')); % add back weighted mean
             weight_factor(i1)=0; % leave out point itself
             dph_smooth2(:,i1)=sum(dph_space.*repmat(weight_factor,n_edge,1),2);
           end
           dph_noise=angle(dph_space.*conj(dph_smooth));
           dph_noise2=angle(dph_space.*conj(dph_smooth2));
           ifg_var=var(dph_noise2,0,1);
           K=lscov(bperp(ifg_index),double(dph_noise)',1./double(ifg_var))'; % estimate arc dem error
           dph_noise=dph_noise-K*bperp(ifg_index)';
           clear dph_space dph_smooth dph_smooth2 dph_noise2
           edge_std=std(dph_noise,0,2);
           edge_max=max(abs(dph_noise),[],2);
           clear dph_noise
        else
           ifg_var=var(dph_space,0,1);
           K=lscov(bperp(ifg_index),double(dph_space)',1./double(ifg_var))'; % estimate arc dem error
           dph_space=dph_space-K*bperp(ifg_index)';
           edge_std=std(angle(dph_space),0,2);
           edge_max=max(abs(angle(dph_space)),[],2);
           clear dph_space
        end

        logit(sprintf('Estimating max noise for all pixels...'))
        ps_std=inf(n_ps,1,'single');
        ps_max=inf(n_ps,1,'single');
        for i=1:n_edge
           ps_std(edgs(i,:))=min([ps_std(edgs(i,:)),[edge_std(i);edge_std(i)]],[],2);
           ps_max(edgs(i,:))=min([ps_max(edgs(i,:)),[edge_max(i);edge_max(i)]],[],2);
        end
        ix_weed2=ps_std<weed_standard_dev&ps_max<weed_max_noise;
        ix_weed(ix_weed)=ix_weed2;
        n_ps=sum(ix_weed);

        logit([num2str(n_ps),' PS kept after dropping noisy pixels']);
    end
end

%% Keep information about number of PS left.
if exist('no_ps_info.mat','file')~=2
   stamps_step_no_ps = zeros([5 1 ]);       % keep for the first 5 steps only
else
   load('no_ps_info.mat');
   % reset as we are currently re-processing
   stamps_step_no_ps(4:end)=0;
end
if n_ps==0
   fprintf('***No PS points left. Updating the stamps log for this****\n')
   % update the flag indicating no PS left in step 3
   stamps_step_no_ps(4)=1;
end
save('no_ps_info.mat','stamps_step_no_ps')

%% Saving the results
weedname=['weed',num2str(psver)];
% save(weedname,'ix_weed','ix_weed2','ps_std','ps_max','ifg_index')
stamps_save(weedname,ix_weed,ix_weed2,ps_std,ps_max,ifg_index)

coh_ps=coh_ps2(ix_weed);
K_ps=K_ps2(ix_weed);
C_ps=C_ps2(ix_weed);
ph_patch=ph_patch2(ix_weed,:);
if ~isempty(ph_res2)
    ph_res=ph_res2(ix_weed,:);
else
    ph_res=ph_res2;
end

pmname=['pm',num2str(psver+1)];
% save(pmname,'ph_patch','ph_res','coh_ps','K_ps','C_ps')
stamps_save(pmname,ph_patch,ph_res,coh_ps,K_ps,C_ps)

clear ph_patch ph_res coh_ps K_ps C_ps ph_patch2 ph_res2 coh_ps2 K_ps2 C_ps2

ph2=ph2(ix_weed,:);
ph=ph2;
phname=['ph',num2str(psver+1)];
% save(phname,'ph')
stamps_save(phname,ph)

clear ph

xy2=xy2(ix_weed,:);
ij2=ij2(ix_weed,:);
lonlat2=lonlat2(ix_weed,:);
ps.xy=xy2;
ps.ij=ij2;
ps.lonlat=lonlat2;
ps.n_ps=size(ph2,1);
psname=['ps',num2str(psver+1)];
save(psname,'-struct','ps');
clear ps xy2 ij2 lonlat2

if exist(hgtname,'file')
    hgt=hgt(ix_weed);
    % save(['hgt',num2str(psver+1),'.mat'],'hgt');
    stamps_save(['hgt',num2str(psver+1),'.mat'],hgt);

    clear hgt
end

if exist(laname,'file')
    la=load(laname);
    la=[la.la(ix2)];
    if all_da_flag~=0
        lao=load(laothername);
        la=[la;lao.la_other(ix_other)];
        clear lao
    end
    la=la(ix_weed);
%    save(['la',num2str(psver+1),'.mat'],'la');
    stamps_save(['la',num2str(psver+1),'.mat'],la);
    clear la
end

if exist(incname,'file')
    inc=load(incname);
    inc=[inc.inc(ix2)];
    if all_da_flag~=0
        inco=load(incothername);
        inc=[inc;inco.inc_other(ix_other)];
        clear inco
    end
    inc=inc(ix_weed);
    stamps_save(['inc',num2str(psver+1),'.mat'],inc);
    clear inc
end

if exist(bpname,'file')
    bp=load(bpname);
    bperp_mat=[bp.bperp_mat(ix2,:)];
    clear bp
    if all_da_flag~=0
        bpo=load(bpothername);
        bperp_mat=[bperp_mat;bpo.bperp_other(ix_other,:)];
        clear bpo
    end
    bperp_mat=bperp_mat(ix_weed,:);
    % save(['bp',num2str(psver+1),'.mat'],'bperp_mat');
    stamps_save(['bp',num2str(psver+1),'.mat'],bperp_mat);
end

if exist(['scla_smooth',num2str(psver+1),'.mat'],'file')
    delete(['scla_smooth',num2str(psver+1),'.mat'])
end
if exist(['scla',num2str(psver+1),'.mat'],'file')
    delete(['scla',num2str(psver+1),'.mat'])
end
if exist(['scla_smooth_sb',num2str(psver+1),'.mat'],'file')
    delete(['scla_smooth_sb',num2str(psver+1),'.mat'])
end
if exist(['scla_sb',num2str(psver+1),'.mat'],'file')
    delete(['scla_sb',num2str(psver+1),'.mat'])
end
if exist(['aps',num2str(psver+1),'.mat'],'file')
    delete(['aps',num2str(psver+1),'.mat'])
end
if exist(['scn',num2str(psver+1),'.mat'],'file')
    delete(['scn',num2str(psver+1),'.mat'])
end

setpsver(psver+1)
logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% Begin initialization code - DO NOT EDIT
function varargout = chanoperGUI(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
        'gui_Singleton',  gui_Singleton, ...
        'gui_OpeningFcn', @chanoperGUI_OpeningFcn, ...
        'gui_OutputFcn',  @chanoperGUI_OutputFcn, ...
        'gui_LayoutFcn',  [] , ...
        'gui_Callback',   []);
if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
        gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

%--------------------------------------------------------------------------
function chanoperGUI_OpeningFcn(hObject, eventdata, handles, varargin)

handles.output = [];
example{1}  = 'ch@ = ch6 - ch5 label HEOG';
example{2}  = 'chan@ = ch6 - ch5 label = HEOG';
example{3}  = 'ch@ = 0.5*chan7 + 0.5*chan9 label LINKM';
example{4}  = 'ch@ = mahaleeg(ch1,ch2) label MAHAL01';
example{5}  = 'CH@ = mahaleeg(ch40) label SELFMAHAL';
example{6}  = 'ch@ = (ch11+ch12+ch13)/3 label = ROI temporal';
example{7}  = 'CHAN@ = (ch11 + ch12 + ch13)/3 label ROI temporal';
example{8}  = 'ch@ = abs(ch110) label E110 rectified ';

try
        ERPLAB = varargin{1};

        if iserpstruct(ERPLAB)
                nchan      = ERPLAB.nchan; % Total number of channels
                typedata   = 'ERP';
                datastr    = 'ERPset';
                formtype   = 'erpchanformulas';
                example{5} = 'ch@ = mgfperp(ERP) label MGFPower';
        else
                example{3} = 'ch@ = chinterpol';
                nchan    = ERPLAB.nbchan; % Total number of channels
                typedata = 'EEG';
                datastr  = 'dataset';
                formtype = 'eegchanformulas';
        end
catch
        ERPLAB.chanlocs = [];
        listch   = '';
        nchan    = 1;
        formtype = [];
        typedata = 'unspecific data';
        datastr  = 'who-knows-what';
end
try
        def = varargin{2};
        formulas = def{1};
        wchmsgon = def{2};
catch
        formulas = [];
        wchmsgon = 1;
end

handles.nchan    = nchan;
% handles.listname = [];
handles.example  = example;
handles.exacounter = 0;
handles.typedata = typedata;
handles.formtype = formtype;

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Channel Operation GUI for ' typedata])

% formulas = erpworkingmemory(formtype);

if isempty(formulas)
        set(handles.editor,'String','');
else
        set(handles.editor, 'String', formulas)
end

% For undo & redo
% handles.formulas{1} = get(handles.editor, 'String');
% handles.undocount = 1;
% set(handles.pushbutton_undo,'Enable','off')
% set(handles.pushbutton_redo,'Enable','off')

%
% Prepare List of current Channels
%
listch=[];
if isempty(ERPLAB.chanlocs)
        for e=1:nchan
                ERPLAB.chanlocs(e).labels = ['Ch' num2str(e)];
        end
end
listch = cell(1,nchan);
for ch =1:nchan
        listch{ch} = [num2str(ch) ' = ' ERPLAB.chanlocs(ch).labels ];
end

set(handles.listboxchan1,'String', listch)
label1 = '<left>Send file rather than individual equations';
label2 = '<left>(creates compact history)';
set(handles.checkbox_sendfile2history, 'string',[label1 '
' label2])
handles.listch = listch;

%
% Mode buttons
%
set(handles.button_recursive, 'String', sprintf('Modify existing %s (recursive updating)', datastr));
set(handles.button_no_recu, 'String', sprintf('Create new %s (independent transformations)',datastr));

%
% Gui memory
%
chanopGUI = erpworkingmemory('chanopGUI');

if isempty(chanopGUI)
        set(handles.button_recursive,'Value', 1); % default is Modify existing ERPset (recursive updating)
        set(handles.button_savelist, 'Enable','off')
        set(handles.chkeeplocs,'Value', 1) %default is Preserve Channel locations

        %
        % File List
        %
        set(handles.edit_filelist,'String','');
        set(handles.checkbox_sendfile2history,'Value',0)
        handles.listname = [];
else
        if chanopGUI.emode==0
                set(handles.button_recursive,'Value', 1);
                set(handles.button_no_recu,'Value', 0);
        else
                set(handles.button_recursive,'Value', 0);
                set(handles.button_no_recu,'Value', 1);
        end
        if chanopGUI.hmode==0
                set(handles.checkbox_sendfile2history,'Value', 0);
        else
                set(handles.checkbox_sendfile2history,'Value', 1);
        end
        if chanopGUI.keeplocs ==1
            set(handles.chkeeplocs,'Value',1);
        else
            set(handles.chkeeplocs,'Value',0);
        end
        listname = chanopGUI.listname;
        set(handles.edit_filelist,'String', listname );
        handles.listname = listname; % JLC Sept 1, 2012
end

% wchmsgon = erpworkingmemory('wchmsgon');

if isempty(wchmsgon) || wchmsgon==0
        set(handles.chwarning,'Value', 0)
elseif wchmsgon==1
        set(handles.chwarning,'Value', 1)
else
        error('Oops...checkbox_warning memory failed')
end

handles.locs = ERPLAB.chanlocs;

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
helpbutton

drawnow
uiwait(handles.gui_chassis);

%--------------------------------------------------------------------------
function varargout = chanoperGUI_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure

varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

%--------------------------------------------------------------------------
function editor_Callback(hObject, eventdata, handles)

set(handles.edit_filelist, 'String','');
compacteditor(hObject, eventdata, handles);

% handles = editorbackup(hObject, eventdata, handles);
% formulas  = get(handles.editor, 'String');
% undocount = handles.undocount;
% undocount = undocount + 1;
% handles.formulas{undocount} = formulas;
% handles.undocount = undocount;

handles.listname  = [];

% Update handles structure
guidata(hObject, handles);

%--------------------------------------------------------------------------
function editor_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function pushbutton_help_Callback(hObject, eventdata, handles)
% doc pop_eegchanoperator
web https://github.com/lucklab/erplab/wiki/EEG-and-ERP-Channel-Operations -browser

%--------------------------------------------------------------------------
function pushbutton_RUN_Callback(hObject, eventdata, handles)
listname    = handles.listname;
compacteditor(hObject, eventdata, handles);
formulalist = get(handles.editor,'String');
wchmsgon    = get(handles.chwarning,'Value');
keeplocs    = get(handles.chkeeplocs,'Value');

if strcmp(formulalist,'')
        msgboxText =  'You have not written any formula!';
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(msgboxText, title);
        return
end
if size(formulalist,2)>256
        msgboxText = ['Formulas length exceed 256 characters.\n\n'...
                'Be sure to press [Enter] after you have entered each formula.'];
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(sprintf(msgboxText), title);
        return
end

%
% Check formulas
%
if get(handles.button_recursive,'Value')
        editormode = 0;
else
        editormode = 1;
end

typedata = lower(handles.typedata);
[option, recall, goeson] = checkformulas(cellstr(formulalist), [typedata 'chanoperGUI'], editormode);

if goeson==0
        return
end
if isempty(listname) && get(handles.checkbox_sendfile2history,'Value')==1
        BackERPLABcolor = [1 0.9 0.3];    % yellow
        question = ['Equations at editor window have not been saved yet.\n\n'...
                'What would you like to do?'];
        title = 'WARNING: Save List of edited chans';
        oldcolor = get(0,'DefaultUicontrolBackgroundColor');
        set(0,'DefaultUicontrolBackgroundColor',BackERPLABcolor)
        button = questdlg(sprintf(question), title,'Save and run','Run without saving', 'Cancel','Run without saving');
        set(0,'DefaultUicontrolBackgroundColor',oldcolor)

        if strcmpi(button,'Save and run')
                fullname = savelist(hObject, eventdata, handles);
                listname = fullname;
                handles.output = {listname, wchmsgon}; % sent filenam string)
        elseif strcmpi(button,'Run without saving')
                handles.output = {cellstr(formulalist), wchmsgon}; % sent like a cell string (with formulas)
        elseif strcmpi(button,'Cancel') || strcmpi(button,'')
                handles.output   = [];
                handles.listname = [];
                % Update handles structure
                guidata(hObject, handles);
                return
        end
elseif isempty(listname) && get(handles.checkbox_sendfile2history,'Value')==0
        handles.output = {cellstr(formulalist), wchmsgon}; % sent like a cell string (with formulas)

elseif ~isempty(listname) && get(handles.checkbox_sendfile2history,'Value')==1
        handles.output = {listname, wchmsgon}; % sent filename string

elseif ~isempty(listname) && get(handles.checkbox_sendfile2history,'Value')==0
        handles.output = {cellstr(formulalist), wchmsgon}; % sent like a cell string (with formulas)
end

% handles.output
% formtype = handles.formtype;
% erpworkingmemory(formtype, formulalist);

%
% memory for Gui
%
chanopGUI.emode = editormode;
chanopGUI.hmode = get(handles.checkbox_sendfile2history,'Value');
chanopGUI.listname  = listname;
chanopGUI.keeplocs = keeplocs;
chanopGUI.chanlocs = handles.locs;
disp(chanopGUI.chanlocs)
erpworkingmemory('chanopGUI', chanopGUI);

try
    handles.output{3} = keeplocs;
catch
        disp('Problem setting location preference')
end

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function cancel_Callback(hObject, eventdata, handles)

handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function eraser_Callback(hObject, eventdata, handles)
% %
% % for undo
% %
% formulas  = get(handles.editor, 'String');
% undocount = handles.undocount;
% handles.formulas{undocount} = formulas;
% undocount = undocount + 1;
% handles.undocount = undocount;
formulas  = get(handles.editor, 'String');
if isempty(formulas)
        return
end

% handles = editorbackup(hObject, eventdata, handles);
set(handles.editor, 'String','')
handles.output = [];
disp('Formulas were erased.')
set(handles.button_savelist, 'Enable','off')
handles.listname = [];
% Update handles structure
guidata(hObject, handles);

%--------------------------------------------------------------------------
function help_Callback(hObject, eventdata, handles)

fprintf('\n\n\n\n\n')
help pop_chanoperator

%--------------------------------------------------------------------------
function button_saveaslist_Callback(hObject, eventdata, handles)

compacteditor(hObject, eventdata, handles);
fulltext = strtrim(get(handles.editor,'String'));

if size(fulltext,2)>256
        msgboxText =  ['Formulas length exceed 256 characters.\n\n'...
                'Be sure to press [Enter] after you have entered each formula.'];
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(sprintf(msgboxText), title);
        return
end
if ~strcmp(fulltext,'')
        fullname = savelist(hObject, eventdata, handles);
        if isempty(fullname)
                return
        end
        set(handles.edit_filelist, 'String', fullname )
        set(handles.button_savelist, 'Enable', 'on')
        handles.listname = fullname;
        % Update handles structure
        guidata(hObject, handles);
else
        set(handles.button_saveaslist,'Enable','off')
        msgboxText =  'You have not written any formula yet!';
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(msgboxText, title);
        set(handles.button_saveaslist,'Enable','on')
        return
end

%--------------------------------------------------------------------------
function button_loadlist_Callback(hObject, eventdata, handles)
[filename, filepath] = uigetfile({'*.txt';'*.*'},'Select a formulas-file');

if isequal(filename,0)
        disp('User selected Cancel')
        return
else
        fullname = fullfile(filepath, filename);
        disp(['pop_chanoperation(): For formulas-file, user selected ', fullname])
end

set(handles.edit_filelist,'String',fullname);
fid_formula = fopen( fullname );

try
        formcell    = textscan(fid_formula, '%s','delimiter', '\r');
        formulas    = char(formcell{:});
catch
        serr = lasterror;
        msgboxText =  ['Please, check your file: \n';
                fullname '\n'...
                serr.message];
        title = 'ERPLAB: pop_chanoperation() error:';
        errorfound(sprintf(msgboxText), title);
        return
end
if size(formulas,2)>256
        msgboxText =  ['Formulas length exceed 256 characters.\n\n'...
                'Be sure to press [Enter] after you have entered each formula.'];
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(sprintf(msgboxText), title);
        return
end

compacteditor(hObject, eventdata, handles);
set(handles.editor,'String',formulas);
fclose(fid_formula);
set(handles.button_savelist, 'Enable','on')
handles.listname = fullname;

% Update handles structure
guidata(hObject, handles);

%--------------------------------------------------------------------------
function listboxchan1_Callback(hObject, eventdata, handles)

numchan   = get(hObject, 'Value');
if isempty(numchan)
        return
end

linet    = get(handles.editor, 'Value');

if nnz(linet) == 0
        linet = 1;
end

formulas = cellstr(get(handles.editor, 'String'));
formulas{linet} = [formulas{linet} 'ch' num2str(numchan)];
set(handles.editor, 'String', char(formulas));

%--------------------------------------------------------------------------
function listboxchan1_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function button_example_Callback(hObject, eventdata, handles)

% %
% % for undo
% %
% formulas  = get(handles.editor, 'String');
% undocount = handles.undocount;
% handles.formulas{undocount} = formulas;
% undocount = undocount + 1;
% handles.undocount = undocount;

nchan      = handles.nchan;
example    = handles.example;
exacounter = handles.exacounter;
exacounter = exacounter + 1;
text       = cellstr(get(handles.editor, 'String'));

if get(handles.button_no_recu,'Value')
        prechar = 'n';
else
        prechar = '';
end
if isempty([text{:}]) || exacounter>length(example)
        exacounter = 1;
end
if length(text)==1 && strcmp(text{1}, '')
        exacurr = char(regexprep(example{exacounter},'@',num2str(nchan+exacounter)));
        text{1} = [prechar exacurr];
else
        exacurr = char(regexprep(example{exacounter},'@',num2str(nchan+exacounter)));
        text{end+1} = [prechar exacurr];
end

set(handles.editor, 'String', char(text));
% handles    = editorbackup(hObject, eventdata, handles);
handles.exacounter = exacounter;
% Update handles structure
guidata(hObject, handles);

%--------------------------------------------------------------------------
function fullname = savelist(hObject, eventdata, handles)
fulltext = char(get(handles.editor,'String'));

if isempty(fulltext)
        return
end

fullnamepre = get(handles.edit_filelist,'String');

%
% Save OUTPUT file
%
[filename, filepath, filterindex] = uiputfile({'*.txt';'*.dat';'*.*'},'Save formulas-file as', fullnamepre);

if isequal(filename,0)
        disp('User selected Cancel')
        fullname = [];
        return
else
        [px, fname, ext] = fileparts(filename);

        if strcmp(ext,'')
                if filterindex==1 || filterindex==3
                        ext   = '.txt';
                else
                        ext   = '.dat';
                end
        end

        fname = [ fname ext];
        fullname = fullfile(filepath, fname);

        fid_list   = fopen( fullname , 'w');

        for i=1:size(fulltext,1)
                fprintf(fid_list,'%s\n', fulltext(i,:));
        end

        fclose(fid_list);
        set(handles.button_savelist, 'Enable','on')
end

%--------------------------------------------------------------------------
function compacteditor(hObject, eventdata, handles)
texteditor = strtrim(get(handles.editor,'String'));

if isempty(texteditor)
        return
end

formul = cellstr(texteditor);
nfl    = length(formul);
k = 1;
formulalist = cell(1);

% removes blank lines
for i=1:nfl
        if ~strcmp(formul{i},'')
                % removes blank lines
                formulalist{k} = regexprep(formul{i},'[^\s](label)[^\s]',' $1 ', 'ignorecase');
                k = k+1;
        end
end

formulalist = strtrimx(formulalist); % special white space removing
set(handles.editor,'String', char(formulalist));
drawnow

%--------------------------------------------------------------------------
function check_aoperations(hObject, eventdata, handles)
texteditor = strtrim(get(handles.editor,'String'));
if isempty(texteditor)
        return
end

formul = cellstr(texteditor);
nfl = length(formul);
k   = 1;
formulalist = cell(1);

for i=1:nfl
        if ~strcmp(formul{i},'')
                formulalist{k} = formul{i};
                k = k+1;
        end
end

formulalist = strtrimx(formulalist);
set(handles.editor,'String', char(formulalist));

%--------------------------------------------------------------------------
function button_savelist_Callback(hObject, eventdata, handles)
compacteditor(hObject, eventdata, handles);
fulltext = strtrim(get(handles.editor,'String'));

if size(fulltext,2)>256
        msgboxText = ['Formulas length exceed 256 characters.\n\n'...
                'Be sure to press [Enter] after you have entered each formula.'];
        title = 'ERPLAB: chanoperGUI few inputs';
        errorfound(sprintf(msgboxText), title);
        return
end
if ~isempty(fulltext)

        fullname = get(handles.edit_filelist, 'String');

        if ~strcmp(fullname,'')
                fid_list   = fopen( fullname , 'w');

                for i=1:size(fulltext,1)
                        fprintf(fid_list,'%s\n', fulltext(i,:));
                end

                fclose(fid_list);
                handles.listname = fullname;
                % Update handles structure
                guidata(hObject, handles);
                disp(['Saving equation list at ' fullname ''])
        else
                button_saveaslist_Callback(hObject, eventdata, handles)
                return
        end
else
        set(handles.button_saveaslist,'Enable','off')
        msgboxText =  'You have not written any formula yet!';
        title      = 'ERPLAB: chanoperGUI few inputs';
        errorfound(msgboxText, title);
        set(handles.button_saveaslist,'Enable','on')
        return
end

%--------------------------------------------------------------------------
function checkbox_sendfile2history_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_filelist_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_filelist_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function button_clearfile_Callback(hObject, eventdata, handles)
set(handles.edit_filelist,'String','');
set(handles.button_savelist, 'Enable', 'off')

%--------------------------------------------------------------------------
function button_recursive_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
        set(handles.button_no_recu,'Value',0)
        val = testsyntaxtype(hObject, eventdata, handles, 'recu');

        if val==0;
                set(handles.button_recursive, 'Value', 0)
                set(handles.button_no_recu, 'Value', 1)
        end

        %%%       handles = editorbackup(hObject, eventdata, handles);
        %%% Update handles structure
        %%guidata(hObject, handles);
else
        set(hObject,'Value',1)
end

%--------------------------------------------------------------------------
function button_no_recu_Callback(hObject, eventdata, handles)
if get(hObject,'Value')
        set(handles.button_recursive,'Value',0)
        val = testsyntaxtype(hObject, eventdata, handles, 'norecu');
        if val==0
                set(handles.button_recursive, 'Value', 1)
                set(handles.button_no_recu, 'Value', 0)
        end
        %%%       handles = editorbackup(hObject, eventdata, handles);
        %%% Update handles structure
        %%%guidata(hObject, handles);
else
        set(hObject,'Value',1)
end

%--------------------------------------------------------------------------
function val = testsyntaxtype(hObject, eventdata, handles, whocall)
val = 1;
formulaArray = get(handles.editor,'String');

if isempty(formulaArray)
        return
else
        formulaArray = strtrim(cellstr(formulaArray));
end
nformulas = length(formulaArray);

[expspliter parts] = regexp(formulaArray, '=','match','split');
ask4fix = 1;
wantfix = 0;
newnumchan = 1;

for t=1:nformulas
        fcomm = formulaArray{t};
        tokcommentb  = regexpi(fcomm, '^#', 'match');  % comment symbol (June 3, 2013)

        if isempty(tokcommentb) % skip comment symbol
                pleft  = regexpi(parts{t}{1}, '(\s*nch[an]*\d+)', 'tokens');
                plcom  = regexpi(parts{t}{1}, '(\s*ch[an]*\d+)', 'tokens');

                if isempty(pleft) &&  ~isempty(plcom) && strcmpi(whocall,'norecu')
                        if ask4fix
                                BackERPLABcolor = [1 0.9 0.3];    % yellow
                                question = ['For non recursive mode, left side of equation\nmust be define as a new channel.\n'...
                                            'For instance, nchan1 = ...\n\n'...
                                            'Do you want that ERPLAB corrects the syntax for you?'];
                                title = 'WARNING: Syntax is not proper for non recursive mode';
                                oldcolor = get(0,'DefaultUicontrolBackgroundColor');
                                set(0,'DefaultUicontrolBackgroundColor',BackERPLABcolor)
                                button = questdlg(sprintf(question), title,'Cancel','No', 'Yes','Yes');
                                set(0,'DefaultUicontrolBackgroundColor',oldcolor)

                                if strcmpi(button,'Yes')
                                        ask4fix = 0;
                                        wantfix = 1;
                                elseif strcmpi(button,'Cancel')
                                        val = 0; % cancel
                                        break
                                else
                                        ask4fix = 0;
                                        wantfix = 0;
                                end
                                %else
                                %      wantfix =1;
                        end
                elseif ~isempty(pleft) && strcmpi(whocall,'recu')
                        if ask4fix
                                BackERPLABcolor = [1 0.9 0.3];    % yellow
                                question = ['For recursive mode, left side of equation cannot\nbe define as a new channel.\n'...
                                            'For instance, you must write chan1 = ...\n\n'...
                                            'Do you want that ERPLAB corrects the syntax for you?'];
                                title = 'WARNING: Syntax is not proper for recursive mode';
                                oldcolor = get(0,'DefaultUicontrolBackgroundColor');
                                set(0,'DefaultUicontrolBackgroundColor',BackERPLABcolor)
                                button = questdlg(sprintf(question), title,'Cancel','No', 'Yes','Yes');
                                set(0,'DefaultUicontrolBackgroundColor',oldcolor)

                                if strcmpi(button,'Yes')
                                        ask4fix = 0;
                                        wantfix =1;
                                elseif strcmpi(button,'Cancel')
                                        val = 0; % cancel
                                        break
                                else
                                        ask4fix = 0;
                                        wantfix = 0;
                                end
                                %else
                                %      wantfix =1;
                        end
                        %else
                        %      wantfix = 0;
                end
                if wantfix && (~isempty(pleft) || ~isempty(plcom))% fixed  (June 3, 2013): JLC
                        fprintf('WARNING: equation %s ', formulaArray{t})
                        if strcmpi(whocall,'recu') % for recursive mode delete the n in nch
                                formulaArray{t} = sprintf('%s = %s', strtrim(regexprep(parts{t}{1}, '^n*','','ignorecase')), strtrim(parts{t}{2}));
                        else
                                formulaArray{t} = sprintf('%s = %s', strtrim(regexprep(parts{t}{1}, '^n*ch(\D*)(\d*)',['nch$1' num2str(newnumchan)],'ignorecase')), strtrim(parts{t}{2}));
                                newnumchan = newnumchan+1;
                        end
                        fprintf('was changed to equation %s \n', formulaArray{t})
                end
        end
end
if val==1
        set(handles.editor,'String', char(formulaArray));
end

%--------------------------------------------------------------------------
function chwarning_Callback(hObject, eventdata, handles)
%
% if get(hObject,'Value')
%       %
%       % Gui memory
%       %
%       erpworkingmemory('wchmsgon', 1);
% else
%       %
%       % Gui memory
%       %
%       erpworkingmemory('wchmsgon', 0);
% end

%--------------------------------------------------------------------------
function pushbutton_export_chan_list_Callback(hObject, eventdata, handles)
list_of_chans = get(handles.listboxchan1,'String');
nloch = length(list_of_chans);

for i=1:nloch
        fprintf('%s\n',list_of_chans{i});
end

fprintf('\n\nEnd of list.\n\n')

%--------------------------------------------------------------------------
function pushbutton_removechan_Callback(hObject, eventdata, handles)
typedata = handles.typedata;
chan2del = deletechanGUI(typedata);
nchan    = handles.nchan;
keeplocs = get(handles.chkeeplocs,'Value');

if ~isempty(chan2del)
        chan2del = unique_bc2(cell2mat(chan2del));
        if max(chan2del)>nchan && nnz(isinf(chan2del))==0
                msgboxText = 'You are specifying non-existing channels!';
                title = 'ERPLAB: chanoperGUI  remove channels';
                errorfound(sprintf(msgboxText), title);
                return
        end
        if min(chan2del)<1 && nnz(isinf(chan2del))==0
                msgboxText = 'You are specifying freaky non-existing channels!';
                title = 'ERPLAB: chanoperGUI  remove channels';
                errorfound(sprintf(msgboxText), title);
                return
        end
        if nnz(isnan(chan2del))>0  || nnz(isinf(chan2del))>0
                msgboxText = ['Channel problem\n'...
                        'You are specifying freaky non-existing channels!'];
                title = 'ERPLAB: chanoperGUI  remove channels';
                errorfound(sprintf(msgboxText), title);
                return
        end
        if strcmpi(typedata, 'EEG')
                eqtn = sprintf('deletechan(%s)', vect2colon(chan2del, 'Delimiter','off'));
        else
                eqtn = sprintf('delerpchan(%s)', vect2colon(chan2del, 'Delimiter','off'));
        end

        wchmsgon    = get(handles.chwarning,'Value');
        handles.output = {{eqtn}, wchmsgon}; % sent like a cell string (with formulas)
        try
            handles.output{3} = keeplocs;
        catch
            disp('Problem setting location preference')
        end
        % Update handles structure
        guidata(hObject, handles);
        uiresume(handles.gui_chassis);
else
        return
end

%--------------------------------------------------------------------------
function pushbutton_rerefwizard_Callback(hObject, eventdata, handles)
% set to non-recursive mode
set(handles.button_recursive, 'Value', 0)
set(handles.button_no_recu, 'Value', 1)

% open reference wizard
formulalist = rerefassistantGUI(handles.nchan, handles.listch);
if isempty(formulalist)
        return
end

compacteditor(hObject, eventdata, handles);
formulas  = get(handles.editor, 'String');

if get(handles.button_no_recu,'Value')
        formulalist = cellstr([formulalist{:}]);
        for t=1:length(formulalist)
                [expspliter parts] = regexp(formulalist, '=','match','split');
                formulalist{t} = sprintf('%s = %s', strtrim(regexprep(parts{t}{1}, '[^n]*ch','nch','ignorecase')), strtrim(parts{t}{2}));
        end
        formulalist = char(formulalist);
end

if isempty(formulas)
        set(handles.editor,'String',formulalist);
else
        formulalist = cellstr(formulalist);
        formulas = cellstr(formulas);
        formt = [formulas ;formulalist];
        formt = char(formt);
        set(handles.editor,'String',formt);
end

% %--------------------------------------------------------------------------
% function pushbutton_undo_Callback(hObject, eventdata, handles)
%
% undocount = handles.undocount;
% if undocount>1
%       undocount = undocount - 1;
%       formulas  = handles.formulas{undocount};
%       handles.undocount = undocount;
%       set(handles.editor, 'String',char(formulas));
%       set(handles.pushbutton_redo,'Enable','on')
%
%       % formulas = get(handles.edit_filelist, 'String');
%       % Update handles structure
%       guidata(hObject, handles);
% end
% if undocount==1
%       set(handles.pushbutton_undo,'Enable','off')
%       set(handles.pushbutton_redo,'Enable','on')
% end

%--------------------------------------------------------------------------
% function pushbutton_redo_Callback(hObject, eventdata, handles)
%
% nformu    = length(handles.formulas);
% undocount = handles.undocount;
% if undocount<nformu
%       undocount = undocount + 1;
%       formulas  = handles.formulas{undocount};
%       handles.undocount = undocount;
%       set(handles.editor, 'String',char(formulas));
%       % formulas = get(handles.edit_filelist, 'String');
%       % Update handles structure
%       guidata(hObject, handles);
% end
% if undocount==nformu
%       set(handles.pushbutton_redo,'Enable','off')
%       set(handles.pushbutton_undo,'Enable','on')
% end

%--------------------------------------------------------------------------
% function handles = editorbackup(hObject, eventdata, handles)
%
% %
% % for undo
% %
% formulas = get(handles.editor, 'String');
% if isempty(formulas)
%       return
% end
%
% lastform = handles.formulas(end);
%
% if strcmp(char(formulas),char(lastform))
%       return
% end
%
% undocount = length(handles.formulas);
% undocount = undocount + 1;
% handles.formulas{undocount} = formulas;
% set(handles.pushbutton_undo,'Enable','on')
% set(handles.pushbutton_redo,'Enable','off')
%
% handles.undocount = undocount;
% % Update handles structure
% %guidata(hObject, handles);
% return

%--------------------------------------------------------------------------
function pushbutton_RUN_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_loadlist_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_savelist_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_saveaslist_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_clearfile_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function checkbox_sendfile2history_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function eraser_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_example_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function help_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_recursive_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function button_no_recu_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function cancel_CreateFcn(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)

if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')
        % The GUI is still in UIWAIT, us UIRESUME
        uiresume(handles.gui_chassis);
else
        % The GUI is no longer waiting, just close it
        delete(handles.gui_chassis);
end

% --- Executes on button press in chkeeplocs.
function chkeeplocs_Callback(hObject, eventdata, handles)
% hObject    handle to chkeeplocs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of chkeeplocs

% --- Executes during object creation, after setting all properties.
function chwarning_CreateFcn(hObject, eventdata, handles)
% hObject    handle to chwarning (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab