text:
function [param]=pp2chp(input)
% [param]=pp2chp(input)
% LQ controller for 2nd order processes.
% This function computes parameters of the controller (r0, r1, q0, q1, p0, p1).
% The dynamic behavoiour of the closed-loop is defined by LQ criterion
% Output of the controller is calculated follows:
%
%            r0 + r1*z^-1               q0 + q1*z^-1
% U(z^-1) = -------------- * W(z^-1) - -------------- * Y(z^-1)
%            p0 + p1*z^-1               p0 + p1*z^-1
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2
% Gs(z^-1) = -----------------------
%             1 + a1*z^-1 + a2*z^-2
%
% Input:
%   input(1:4) ... [a1 b1 a2 b2]
%   input(5) ... type of reference signal (1-step, 2-ramp, 3-sin)
%   input(6) ... frequency [Hz] (used if reference signal is sin wave)
%   input(7) ... sample time (used if reference signal is sin wave)
%   input(8) ... penalisation of u (constant fi in LQ criterion)
% Output: param ... controller parameters  [r0; r1; q0; q1; p0; p1];

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
rs_type = input(5);
rs_freq = input(6);
T0 = input(7);
fi = input(8);

% LQ: AP+BQ = D where D determined by spectral factorization
% A(z^-1)*fi*A(z) + B(Z^-1)*B(z) = D(z^-1)*delta*D(z)
%          M(z^-1)*M(z)          = D(z^-1)*delta*D(z)
m0 = fi*(1+a1^2+a2^2) + b1^2 + b2^2;
m1 = fi*a1*(1+a2) + b1*b2;
m2 = fi*a2;

la = m0/2 - m2 + sqrt( (m0/2+m2)^2 - m1^2 );
delta = (la + sqrt(la^2-4*m2^2)) / 2;
d2 = m2 / delta;
d1 = m1 / (delta+m2);
d0 = 1;
d3 = 0;

% FBFW controller: Y=BR/(AP+BQ)*W
% conditions: 1) AP+BQ=D
%             2) BR+FS=D  where W=H/F and S is any polynomial
% 1st condition:
%  A =  1 + a1*z^-1 + a2*z^-2   B =      b1*z^-1 + b2*z^-2
%  P = p0 + p1*z^-1             Q = q0 + q1*z^-1
% system of linear equations:
% [ 1  0   0  0]   [p0]   [d0]
% [a1  1  b1  0]   [p1]   [d1]
% [a2 a1  b2 b1] * [q0] = [d2]
% [ 0  a2  0 b2]   [q1]   [d3]
p0 = d0;
p1 = (d0*b2^2*a1-d0*b2*b1*a2-b2^2*d1+b1*b2*d2-b1^2*d3)/(-b2^2+a1*b1*b2-a2*b1^2);
q0 = -(d0*a1^2*b2-d0*a1*b1*a2-d0*a2*b2-d1*a1*b2+d1*b1*a2+b2*d2-b1*d3)/(-b2^2+a1*b1*b2-a2*b1^2);
q1 = -(d0*a2*a1*b2-d0*a2^2*b1-a2*b2*d1+b1*a2*d2+d3*b2-d3*b1*a1)/(-b2^2+a1*b1*b2-a2*b1^2);

% 2nd condition:
switch (rs_type)
case 1,     %step:
    % B = b1*z^-1 + b2*z^-2   F =  1 - z^-1
    % R = r0                  S = s0 + s1*z^-1 + s2*z^-2
    % [ 0    1    0   0]   [r0]   [d0]
    % [b1   -1    1   0]   [s0]   [d1]
    % [b2    0   -1   1] * [s1] = [d2]
    % [ 0    0    0  -1]   [s2]   [d3]
    r0 = (d0+d1+d2+d3)/(b1+b2);
    r1 = 0;
case 2,     %ramp
    % B =      b1*z^-1 + b2*z^-2   F =  1 -  2*z^-1 + z^-2
    % R = r0 + r1*z^-1             S = s0 + s1*z^-1
    % [ 0   0   1   0]   [r0]   [d0]
    % [b1   0  -2   1]   [r1]   [d1]
    % [b2  b1   1  -2] * [s0] = [d2]
    % [ 0  b2   0   1]   [s1]   [d3]
    r0 = (2*d0*b1+3*d0*b2+d1*b1+2*d1*b2+b2*d2-b1*d3)/(b1+b2)^2;
    r1 = -(d0*b1+2*d0*b2+d1*b2-b1*d2-2*b1*d3-d3*b2)/(b1+b2)^2;
case 3,     %sin(om*t/T0)
    om = 2*pi*rs_freq*T0;
    % B =      b1*z^-1 + b2*z^-2   F = 1 -  2*cos(om)*z^-1 + z^-2
    % R = r0 + r1*z^-1             S = s0 + s1*z^-1
    % [ 0   0           1           0]   [r0]   [d0]
    % [b1   0  -2*cos(om)           1]   [r1]   [d1]
    % [b2  b1           1  -2*cos(om)] * [s0] = [d2]
    % [ 0  b2           0           1]   [s1]   [d3]
    r0 = (2*d0*b1*cos(om)+2*d0*b2*cos(2*om)+d0*b2+d1*b1+2*d1*b2*cos(om)+b2*d2-b1*d3)/(b1^2+2*b1*b2*cos(om)+b2^2);
    r1 = (-d0*b1-2*d0*b2*cos(om)-b2*d1+b1*d2+2*d3*b1*cos(om)+d3*b2)/(b1^2+2*b1*b2*cos(om)+b2^2);
end;

param=[r0; r1; q0; q1; p0; p1];

author: Sable

repo: mcbench-benchmarks


text:
function pvpmod(x)
% PVPMOD             - evaluate parameter/value pairs
% pvpmod(x) assigns the value x(i+1) to the parameter defined by the
% string x(i) in the calling workspace. This is useful to evaluate
% <varargin> contents in an mfile, e.g. to change default settings
% of any variable initialized before pvpmod(x) is called.
%
% (c) U. Egert 1998

%############################################
% this loop is assigns the parameter/value pairs in x to the calling
% workspace.

if ~isempty(x)
   for i = 1:2:size(x,2)
      assignin('caller', x{i}, x{i+1});
   end;
end;

%############################################


text:
function [topics, w] = createsyntheticdata (nu, eta, K, W, L)

  % Get the number of documents to generate.
  D = length(L);

  % Generate the word proportions for each topic from the prior.
  topics = zeros(W,K);
  for k = 1:K
    topics(:,k) = dirichletrnd(eta);
  end

  % Generate the documents.
  w = cell(D,1);
  for d = 1:D

    % Generate the topic proportions from the prior.
    t = dirichletrnd(nu)';

    % Generate the topic associated with each word.
    k = sample(t,L(d));

    % Generate the word tokens.
    w{d} = int32(sample_vector(topics,k));
  end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
classdef (Abstract) mme_load < mp.mm_element

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

%     properties
%     end

    methods
        function name = name(obj)
            name = 'load';
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function result = isintersect_3dof(SHAPE, LINE)
% isintersect(SHAPE, LINE)
% this function check whether we intersect the shape or not
% SHAPE could be polygon or line
% LINE is always line
%

vert_num = size(SHAPE(:,1),1);
result = 0;

m = (LINE(2, 2) - LINE(1, 2)) /  (LINE(2, 1) - LINE(1, 1));
b =  LINE(1, 2) - m * LINE(1, 1);
radius = 0.1;

for k = 1:vert_num-1
    % y = m * x + b

    m_obs = (SHAPE(k, 2) - SHAPE(k+1, 2)) / (SHAPE(k, 1) - SHAPE(k+1, 1));
    b_obs = SHAPE(k, 2) - m_obs * SHAPE(k, 1);

    % consider this lines ???

    if (m_obs + radius > m) && (m_obs - radius < m) && (b_obs + radius > b) && (b_obs - radius < b)
        result = 1;
        return;
    end

    % ???

    x_intersection = (b - b_obs)/(m_obs - m);

    x_max = max (SHAPE(k, 1),SHAPE(k+1, 1));
    x_min = min (SHAPE(k, 1), SHAPE(k+1, 1));

    y_max = max (SHAPE(k, 2), SHAPE(k+1, 2));
    y_min = min (SHAPE(k, 2), SHAPE(k+1, 2));

    x_edge_max = max(LINE(1, 1), LINE(2, 1));
    x_edge_min = min(LINE(1, 1), LINE(2, 1));

    y_edge_max = max(LINE(1, 2), LINE(2, 2));
    y_edge_min = min(LINE(1, 2), LINE(2, 2));

    y_intersection = m_obs * x_intersection + b_obs;

    if ((x_intersection >= (x_min-radius) && x_intersection <= (x_max + radius)) ...
            && (x_intersection >= (x_edge_min - radius) && x_intersection <= (x_edge_max+radius)) ...
            && y_intersection<= y_edge_max + radius && y_intersection >= y_edge_min - radius ...
            && y_intersection<= y_max+radius && y_intersection >= y_min - radius)
        result = 1;
        return;
    end
end

author: olzhas


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    y_pred = X * theta; % X.shape=(m,3) theta.shape=(3,1)
    theta = theta - alpha * 1/m * (X'*(y_pred-y));  % y.shape=(m,1);

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: imLogM


text:
function coordsys=m_coord(coordsys)
% M_COORD Initializes the coordinate system for varous conversions.
%
%   M_COORD('set') tells you the current coordinate system
%   M_COORD('get') gives you the current possibilities
%   M_COORD(SYSTEM) sets the coordinate system to SYSTEM.
%
%   Currently the available coordinate systems are:
%              'geographic' (usual lat/long)
%              'geomagnetic' (referenced to magnetic poles)

% Rich Pawlowicz (rich@ocgy.ubc.ca) nOV/2001
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

global MAP_COORDS

if nargin==0, coordsys='usage'; end

coordsys=lower(coordsys);

coords=mc_coords('name');

switch coordsys

   case 'get'
      disp(' ');
      disp('Available coordinate systems are:');
      for k=1:length(coords.name)
        disp(['     ' coords.name{k}]);
      end

   case 'set'
      if isempty(MAP_COORDS)
         disp('No coordinate system initialized');
         m_coord('usage');
      else
         if nargout==0
           disp(MAP_COORDS.name);
         else
	       coordsys=MAP_COORDS;
         end
      end

   case 'usage'
      disp(' ');
      disp('Possible calling options are:');
      disp('  ''usage''                    - this list');
      disp('  ''set''                      - list of coordinate systems');
      disp('  ''get''                      - get current coordinate (if defined)');
      disp('  ''system''                   - initialize coordinate system\n');

   otherwise
     k=strcmp(coordsys,lower(coords.name));
     MAP_COORDS=mc_coords('parameters',coords.name{k});

end
    % Check to see if a non-lat/long coordinate system is being used.

author: goGPS-Project


text:
function value = year_length_months_greek ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_GREEK returns the number of months in a Greek year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  if ( year_is_embolismic_greek ( y ) )
    value = 13;
  else
    value = 12;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym cosh (@var{x})
%% Symbolic cosh function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = cosh (x)
%%   @result{} y = (sym) cosh(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = cosh(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('cosh', x);
end

%!error cosh (sym(1), 2)
%!assert (isequaln (cosh (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = cosh(x);
%! f2 = cosh(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = cosh(A);
%! f2 = cosh(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = cosh (d);
%! f = cosh (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
classdef net_model_dc < mp.net_model & mp.form_dc

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

    properties
        va = [];
        z  = [];
    end

    methods
        %% constructor
        function obj = net_model_dc()
            obj@mp.net_model();
            obj.element_classes = { ...
                @mp.nme_bus_dc, @mp.nme_gen_dc, @mp.nme_load_dc, ...
                    @mp.nme_branch_dc, @mp.nme_shunt_dc };

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if obj.node is empty) in BUILD() and DISPLAY(),
            %%              after object construction, but before object use.
        end

        function obj = def_set_types(obj)
            def_set_types@mp.net_model(obj);        %% call parent first
            obj.set_types.va = 'VOLTAGE VARS (va)';
            obj.set_types.z  = 'NON-VOLTAGE VARS (z)';
        end

        function obj = build_params(obj, nm, dm)
            %% call parent to build individual element parameters
            build_params@mp.net_model(obj, nm, dm);

            %% aggregate parameters from individual elements
            obj.B = obj.stack_matrix_params('B', 1);
            obj.K = obj.stack_matrix_params('K', 0);
            obj.p = obj.stack_vector_params('p');
        end

        function obj = port_inj_soln(obj)
            %% compute port injections
            obj.soln.gp = obj.port_inj_power(obj.soln.x);
        end
    end     %% methods
end         %% classdef


function vl_testder(g,x,dzdy,dzdx,delta,tau)

if nargin < 5
  delta = 1e-3 ;
end

if nargin < 6
  tau = [] ;
end

dzdy = gather(dzdy) ;
dzdx = gather(dzdx) ;

y = gather(g(x)) ;
dzdx_=zeros(size(dzdx));
for i=1:numel(x)
  x_ = x ;
  x_(i) = x_(i) + delta ;
  y_ = gather(g(x_)) ;
  factors = dzdy .* (y_ - y)/delta ;
  dzdx_(i) = dzdx_(i) + sum(factors(:)) ;
end
vl_testsim(dzdx, dzdx_, tau);


text:
% Copyright 2014 Jason Heeris, jason.heeris@gmail.com
%
% This file is part of the gammatone toolkit, and is licensed under the 3-clause
% BSD license: https://github.com/detly/gammatone/blob/master/COPYING
function test_ERBFilterBank()

    erb_space_inputs = { ...
        100, 11025,  10, sin(2*pi*220*[0:22050/100]'/22050); ...
         20, 22050,  10, square(2*pi*150*[0:44100/200]'/44100); ...
         20, 44100,  40, square(2*pi*12000*[0:88200/400]'/88200); ...
        100, 11025, 1000, sawtooth(2*pi*10100*[0:22050/100]'/22050, 0.5); ...
        500, 80000,  200, sawtooth(2*pi*3333*[0:160000/400]'/160000, 0.5); ...
    };

    erb_filter_inputs = { ...
        44100, [22050; 2205; 220], square(2*pi*220*[0:44100/200]'/44100); ...
        16000, [8000; 7000; 6000; 5000; 4000; 3000; 2000; 1000], square(2*pi*2000*[0:16000/50]'/16000); ...
        16000, [16000; 8000; 1], square(2*pi*880*[0:16000/50]'/16000); ...
    };

    num_tests = size(erb_space_inputs)(1) ...
                + size(erb_filter_inputs)(1);

    erb_filterbank_inputs = {};

    erb_filterbank_results = {};

    % This will ONLY generate tests that use the centre frequency inputs

    % ERBSpace generated inputs
    for tnum=1:size(erb_space_inputs)(1)
        [f_low, f_high, num_f, wave] = deal(erb_space_inputs{tnum,:});
        fs = f_high*2;
        f_arr = ERBSpace(f_low, f_high, num_f);
        fcoefs = MakeERBFilters(fs, f_arr, 0);
        erb_filterbank_inputs(tnum, :) = {fcoefs, wave};
    end

    % MakeERBFilters generated inputs
    for tnum=1:size(erb_filter_inputs)
        [fs, f_arr, wave] = deal(erb_filter_inputs{tnum,:});
        fcoefs = MakeERBFilters(fs, f_arr, 0);
        offset = size(erb_space_inputs)(1);
        erb_filterbank_inputs(offset+tnum, :) = {fcoefs, wave};
    end

    for tnum=1:num_tests
        fcoefs = erb_filterbank_inputs{tnum, 1};
        wave = erb_filterbank_inputs{tnum, 2};
        erb_filterbank_results(tnum, :) = ERBFilterBank(wave, fcoefs);
    end

    results_file = fullfile('..', 'tests', 'data', 'test_filterbank_data.mat');
    save(results_file, 'erb_filterbank_inputs', 'erb_filterbank_results');
end

author: detly

repo: gammatone

save_path: github-repos/MATLAB/detly-gammatone


text:
function im = ifft2c(d)
% Function performs a centered ifft2
im = ifftshift(ifft2(ifftshift(d)));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

rowsz = size(X, 1);
colsz = size(X, 2);
mu = mean(X);
sigma = std(X);
mu_matrix = (ones(size(X))) * ((ones(colsz, 1) * mu) .* eye(colsz));
sigma_matrix = ((ones(colsz, 1) * (1./sigma)) .* eye(colsz));
X_norm = (X-mu_matrix) * sigma_matrix;

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

end

author: zlotus


text:
function L = log_marg_prob_node(CPD, self_ev, pev)
% LOG_MARG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m)) for node i (linear_gaussian)
% L = log_marg_prob_node(CPD, self_ev, pev)
%
% This differs from log_prob_node because we integrate out the parameters.
% self_ev{m} is the evidence on this node in case m.
% pev{i,m} is the evidence on the i'th parent in case m
% We assume there is <= 1 case.

ncases = length(self_ev);

if ncases==0
  L = 0;
  return;
elseif ncases==1
  y = self_ev{1};
  x = cat(1, pev{:}); % column vector
  f = 1-x'*inv(x*x' + CPD.prior.n)*x;
  alpha = CPD.prior.alpha;
  L = log_student_pdf(y, x'*CPD.prior.theta, f*alpha/CPD.prior.beta, 2*alpha);
else
  error('can''t handle batch data');
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Add additional initialization if you want to.
% You can use this space to pre-compute the trajectory to avoid
% repeatedly computing the same trajectory in every call of the
% "trajectory_generator" function

% Generate trajectory
disp('Generating Trajectory ...');
trajectory_generator([], [], map, path);

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function rs = rms(s)

%tstoolbox/@signal/rms
%   Syntax:
%     * rs = rms(s)
%
%   Calculate root mean square value for signal along dimension 1.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

c = rms(s.core); 	% call real working routine for parent core object
rs = signal(c, s);

if ndim(rs) > 0
	rs.xaxes = s.xaxes(2:end);
end

rs = addhistory(rs,  ['Calculated root mean square value along dimension 1'] );
rs = addcommandlines(rs, 's = rms(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function value = p01_i4 ( action, name, value )

%*****************************************************************************80
%
%% P01_I4 sets or gets I4 parameters for problem 01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P01_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P01_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P01_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P01_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P01_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P01_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P01_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P01_I4 - Fatal error!' );

  end

  return
end


text:
function data_set = create_interaction_variables(data_set,vars,range_nways,separator,max_varname_length)
% create_interaction_variables
%
% for Matlab R13+
% version 1.1 (April 2012)
% (c) Brian Weidenbaum
% website: http://www.BrianWeidenbaum.com/.
%
%
% OUTPUT: your dataset (or a dataset based on your matrix), updated with new, aptly-named *unique* interaction variables,
% ranging from at least 2 to any number the user specifies
%
% INPUTS (** = OPTIONAL)
% input name: (input datatype/s) -- description
%
% data_set: (dataset OR matrix) -- the data you want to alter
%
% **vars: (cell array of chars/numbers, OR vector of numbers, OR 'ALL') --
% default: 'ALL'
% the names of the variables you want to interact; alternatively, the column numbers of the variables you want to interact
% OR, you can just say 'ALL' to include all variables automatically
%
% **range_nways (vector OR 'MAX') --
% default: 2
% the range of numbers of variables to include in the interaction terms generated by this function
% alternatively, just type 'MAX' to use 2 variables to the maximum possible number of vars
%
% **separator (char) --
% default: '_'
% the separator string you want to use to divide the
% variable names that contributed to a new interaction variable.  Default
% is '_'.  E.g., by default, an interaction between Var1 and Var2 will be
% named 'Var1_Var2'.
%
% **max_varname_length (number) --
% default: 63
% the maximum length of the newly created interaction terms' variable names.
% Any dynamically generated variable names (e.g. 'Var1_Var2') that exceed this number will be excluded from the new dataset.
% You should set max_varname_length according to the database you plan to use with your data--
% e.g., if you only want to use this data in MATLAB, you should set max_varname_length=63 (the maximum length supported by the dataset class,
% but if you plan to export your data to Oracle, you should set it to around 30
%
%
%
% EXAMPLES
%
% You have a dataset with 3 variables: a, b, and c.
% You want to create interaction terms, up to 3 ways, for all of your variables.
% You type:
% new_dataset = create_interaction_variables(old_dataset,'all','max');
% new_dataset will contain:
% a.*b, named 'a_b'
% a.*c, named 'a_c'
% b.*c, named 'b_c'
% and a.*b.*c, named 'a_b_c',
% plus all original variables.
% It will NOT contain b.*a, c.*a, etc because these are not unique combos.
%
% You have a dataset with 3 variables: a, b, and c.
% You want to create interaction terms, up to 2 ways, only for columns 2 and 3.
% You type:
% new_dataset = create_interaction_variables(old_dataset,[2 3],2);
% new_dataset will contain only b.*c,
% plus all original variables.
%
% You have a dataset with 3 variables: a, b, and c.
% You want to create interaction terms, up to 2 ways, only for 'a' and 'c'.
% You type:
% new_dataset = create_interaction_variables(old_dataset,{'a','c'},2);
% new_dataset will contain only a.*c,
% plus all original variables.
%
% You have a dataset with 3 variables: a, b, and c.
% You want to create interaction terms, ONLY 3 ways, only for all vars.
% You type:
% new_dataset = create_interaction_variables(old_dataset,'all',3);
% new_dataset will contain only a .* b .* c,
% plus all original variables.
%
%
%
% CHANGE LOG
% Changes between 1.0 and 1.1:
% -added 'separator' parameter, enabling users to create custom named vars
% -changed 'max_nways' parameter to 'range_nways', giving user more control
% over the types of interactions created
% -enforced maximum max_varname_length as 63
% -reduced min number of arguments to 1, setting defaults for other 4 params
% -minor performance tweaks via vectorizing some inner function loops
%
%

	%PHASE ONE: USER INPUT VALIDATION

    % set default values for all params except data_set
    if nargin==1
       vars='all';
       range_nways=2;
       separator='_';
       max_varname_length = 63; %max length supported by MATLAB dataset class
    elseif nargin==2
       range_nways=2;
       separator='_';
       max_varname_length = 63;
    elseif nargin==3
        separator='_';
        max_varname_length = 63;
    elseif nargin==4
        max_varname_length = 63;
    end

    %if data_set param is a regular matrix, convert it to a dataset
	%Otherwise, user will be unable to know which columns correspond to which interactions
	if ~strcmp(class(data_set),'dataset')
	   data_set = matrix2dataset(data_set);
	   %this data_set will include default variable names, like 'Var1', 'Var2', etc
    end

    %check if user input 'all' or 'ALL' for vars param
    %if so, convert vars = data_set.Properties.VarNames;
    if ischar(vars)
        if strcmpi(vars,'ALL')
            vars = data_set.Properties.VarNames;
        else
            error('If you are using a char array as your vars parameter, you can only set it as ''all'' or ''ALL''. ');
        end
    end

    nvars = length(vars);
    if nvars<2
       error('You must specify at least 2 variable names/column numbers inside the vars parameter.');
    end

    %validate range_nways
    if ischar(range_nways)
        range_nways = 2:nvars;
    elseif isnumeric(range_nways)
        if length(range_nways)<1
           error('The max_nways parameter should contain at least one number.');
        end
        % 2 <= max vars per interaction <= total number of variable names passed in
        % set any nways <2 =2 or >nvars = nvwars and drop all but one
        % cannot have any repeating numbers, so in end just replace with unique
        if sum(range_nways<2)>0
           disp(['At least one number in range_nways is < 2.' ...
               '  Since you must have at least 2 variables per interaction, all numbers < 2 have been replaced with 2.']);
        end
        range_nways(range_nways<2)=2;
         if sum(range_nways>nvars)>0
            disp(['At least one number in range_nways is > ' num2str(nvars)...
                ' (the total number of variables you are trying to create interactions from).' ...
                '  All such numbers will be dropped from range_nways and replaced with a single ' num2str(nvars) '.']);
         end
        range_nways(range_nways>nvars)=nvars;
        range_nways=unique(range_nways);
        %force rnways to be a row vector
        range_nways = reshape(range_nways,1,length(range_nways));
    else
        error('range_nways parameter MUST be a vector of numbers, or ''ALL''.');
    end

	%check if each variable name/column number in vars param exists inside the dataset
	%also, if necessary, convert column numbers to variable names, for future use with eval strings
    if iscell(vars)
      if ~ischar(vars{1}) %if the varnames are NOT chars
        vars = columns2names(vars);
      else %vars are chars, so check if each name is member of vars
            for i=1:length(vars)
                if ~ismember(vars{i},data_set.Properties.VarNames)
                    disp([vars{i} ' is not a valid column and will not be included in the interactions.' ]);
                    nvars=nvars-1;
                end
            end
            vars(~ismember(vars,data_set.Properties.VarNames))=[];
      end
    else %its for sure a vector of numbers
        vars = columns2names(vars);
    end

    %validate separator
    if ~ischar(separator)
       error('Seperator needs to be a char array.');
    elseif length(separator)>60
        error('Your separator term is too long.');
    elseif ismember('%',separator)
        error('Your separator contains one or more illegal characters.');
    end

    %validate max_varname_length
    if max_varname_length>63
       error('Max varname length exceeds 63, which is the highest number supported by the MATLAB dataset class.');
    end

	%END INPUT VALIDATION

	%PHASE TWO: CREATE ALL THE INTERACTION VARIABLES
	all_evalstrs = get_evalstrs;
	for i=1:length(all_evalstrs)
	%trycatch prevents multiplying chars by numbers
		try
			eval(all_evalstrs{i});
		catch e
            disp(e.message);
		end
    end

    %PHASE THREE: PROFIT

	%
	%
	% INNER FUNCTIONS
	%
	%

	%inner function that translates varname INDICES to actual varnames cellarray of chars
	%for use in eval() function
    function char_varnames = columns2names(columns)
        %first detect if the columns are a cell array; if so, xform to vector
        if iscell(columns)
            columns = cell2mat(columns);
        end

        %check if all column numbers are between 1 and total nvars in data_set
        for ii=1:length(columns)
         if ((columns(ii)<1) || (columns(ii)>size(data_set,2)))
            disp(['Column number ' num2str(columns(ii)) ' is not a valid column and will not be included in the interactions.' ]);
            nvars=nvars-1;
         end
        end
        columns(columns(:)<1 | columns(:)>size(data_set,2))=[];
        char_varnames = data_set.Properties.VarNames(columns);
    end %inner function

	% inner function that creates all UNIQUE combinations of indices for 2:maxvars_perinteraction
	% adds them all to a cell array
    function all_indices = get_interaction_indices()
        all_indices={};
        ct=0;
        %for each number of vars per interaction...
        %get all unique combos of indices
        %e.g. if 3 vars, and 3 max vars per interaction
        %first, get 2 var combos:
        %12,13,23
        %then get 3 var combo: 123
        %add all that to cell array and return it
        for nways = range_nways

            %need to dynamically create allcomb args
            %if it's choose2, there are 2 args: 1:nvars,1:nvars
            %if it's choose3, there are 3 args: 1:nvars,1:nvars,1:nvars ...etc
            allcomb_argstr = repmat('1:nvars,',1,nways);
            %remove final ','
            allcomb_argstr  = allcomb_argstr(1:end-1);

            %allcomb returns cartesian product
            all_combos = eval(['allcomb(' allcomb_argstr ');']);

            %now choose all appropriate combos
            %so we only have unique combinations (eg not 2 1 and 1 2)
            logical_str='';
            for choose_columns= nways:-1:2
               logical_str=[logical_str, '(all_combos(:,' , num2str(choose_columns) , ')>all_combos(:,' , num2str(choose_columns-1) , '))&' ];
            end
            %remove final '&'
            logical_str = logical_str(1:end-1);

            % below gives the combo indices
            good_combos = all_combos(eval(logical_str),:);

            %now add each row in good_combos to all_indices
            ngood = size(good_combos,1);
            for ii=1:ngood
                ct =ct+1;
                all_indices{ct} = good_combos(ii,:);
            end %adding rows to all_indices
        end %loop that goes thru each possible number of vars per interaction
    end %get_interaction_indices inner fx

    % inner function that gets the 'data_set.var1_var2_varN=data_set.var1.*data_set.var2.*datset.varN' version of indices
    function str = indices2str(indices)
        % goal: if give two indices [1 2], return the following string:
        %   'data_set.var1_var2 = data_set.var1.*data_set.var2;'

        n_indices = length(indices);
        %left of equals sign should be: 'data_set.var1_var2'
        left_str = 'data_set.';

        %add each varname separated by _
        for i=1:n_indices
            left_str = [left_str vars{indices(i)} separator]; %seperator var
        end %building varname string
        %remove the final separator
        left_str = left_str(1:end-length(separator));
        %check if varname is too big
        if length(left_str)-length('data_set.') > max_varname_length
           error(['Name of new variable: ''' left_str(length('data_set.')+1:end) ...
               ''' exceeded max varname length of ' num2str(max_varname_length) ', and will not be included in the new dataset.']);
        end

        %now for the righthand side of equals, which should be:
        %'data_set.var1.*data_set.var2;'
        right_str = '';
        for i=1:n_indices
            right_str = [right_str 'data_set.' vars{indices(i)} '.*'];
        end %building varname string
        %removing the final '.*'
        right_str=[right_str(1:end-2) ';'];

        str= [left_str '=' right_str];
    end %indices2str inner fx

    %will give cellarray where each cell contains a good evalstr
    function all_evalstrs = get_evalstrs
        %first get all the indices
        inds = get_interaction_indices;
        ncombos= length(inds);
        all_evalstrs = cell(ncombos,1);

        %go thru each combo of indices and try to create an eval str
        %need try catch bc of potential max varname length problems
        for ii = 1:ncombos
            try
                all_evalstrs{ii}=indices2str(inds{:,ii});
            catch e
                all_evalstrs{ii}='';
                disp(e.message);
            end
        end
    end %get all evalstrs

    function A = allcomb(varargin)
        % ALLCOMB - All combinations
        %    B = ALLCOMB(A1,A2,A3,...,AN) returns all combinations of the elements
        %    in A1, A2, ..., and AN. B is P-by-N matrix is which P is the product
        %    of the number of elements of the N inputs.
        %    Empty inputs yields an empty matrix B of size 0-by-N. Note that
        %    previous versions (1.x) simply ignored empty inputs.
        %
        %    Example:
        %       allcomb([1 3 5],[-3 8],[0 1]) ;
        %         1  -3   0
        %         1  -3   1
        %         1   8   0
        %         ...
        %         5  -3   1
        %         5   8   0
        %         5   8   1
        %
        %    ALLCOMB(A1,..AN,'matlab') causes the first column to change fastest.
        %    This is more consistent with matlab indexing. Example:
        %    allcomb(1:2,3:4,5:6,'matlab') %->
        %      1   3   5
        %      2   3   5
        %      1   4   5
        %      ...
        %      2   4   6
        %
        %    This functionality is also known as the cartesian product.
        %
        %    See also NCHOOSEK, PERMS, NDGRID
        %    and COMBN, KTHCOMBN (Matlab Central FEX)

        % for Matlab R13+
        % version 2.2 (jan 2012)
        % (c) Jos van der Geest
        % email: jos@jasen.nl

        % History
        % 1.1 (feb 2006), removed minor bug when entering empty cell arrays;
        %     added option to let the first input run fastest (suggestion by JD)
        % 1.2 (jan 2010), using ii as an index on the left-hand for the multiple
        %     output by NDGRID. Thanks to Jan Simon, for showing this little trick
        % 2.0 (dec 2010). Bruno Luong convinced me that an empty input should
        % return an empty output.
        % 2.1 (feb 2011). A cell as input argument caused the check on the last
        %      argument (specifying the order) to crash.
        % 2.2 (jan 2012). removed a superfluous line of code (ischar(..))

        error(nargchk(1,Inf,nargin)) ;

        % check for empty inputs
        q = ~cellfun('isempty',varargin) ;
        if any(~q),
            warning('ALLCOMB:EmptyInput','Empty inputs result in an empty output.') ;
            A = zeros(0,nargin) ;
        else

            ni = sum(q) ;

            argn = varargin{end} ;

            if ischar(argn) && (strcmpi(argn,'matlab') || strcmpi(argn,'john')),
                % based on a suggestion by JD on the FEX
                ni = ni-1 ;
                ii = 1:ni ;
                q(end) = 0 ;
            else
                % enter arguments backwards, so last one (AN) is changing fastest
                ii = ni:-1:1 ;
            end

            if ni==0,
                A = [] ;
            else
                args = varargin(q) ;
                if ~all(cellfun('isclass',args,'double')),
                    error('All arguments should be arrays of doubles') ;
                end
                if ni==1,
                    A = args{1}(:) ;
                else
                    % flip using ii if last column is changing fastest
                    [A{ii}] = ndgrid(args{ii}) ;
                    % concatenate
                    A = reshape(cat(ni+1,A{:}),[],ni) ;
                end
            end
        end
    end %allcomb inner fx

	function output = matrix2dataset(matrix,ifx_varnames)
	% matrix2dataset
	%
	% INPUTS
	% matrix (N x M matrix): the matrix version of the data you want to convert
	% varnames (OPTIONAL cell vector of chars): the names of the columns, in order
	%
	% OUTPUT:
	% a dataset that contains all the data from your matrix, named according to
	% your varnames
	% assumes that each column in matrix is one dataset column-- eg no dataset
	% column should consist of >1 matrix column
	%
	%
	% EXAMPLE:
	% This function is great for when you have converted an old dataset
	% into a matrix of doubles for performing some calculations, but now you
	% want to update your original dataset with new values or new columns that
	% you have created from your calculations.
	%
	% In that situation, you would type the following:
	% updated_dataset = matrix2dataset(matrix_from_old_dataset, old_dataset.Properties.VarNames);
	%
	% If you don't have any varnames, you can choose not to use the varnames parameter
	% In that case, you will get the default dataset names: Var1, Var2,..VarN
	%

		n_vars = size(matrix,2);
		if nargin==2
			if ~(length(ifx_varnames)==n_vars)
				error('You must have one variable name for each column in the matrix parameter.');
			end
		end

		str_vars = '';
		for ii=1:n_vars
		   str_vars = [str_vars, 'matrix(:,', num2str(ii), '),' ];
		end

		if nargin == 2
			ds_arguments = [str_vars, '''VarNames'',','ifx_varnames'];
		elseif nargin==1
			ds_arguments = str_vars(1:end-1);
		end
		output = eval(['dataset(' ds_arguments ');']);
	end%innerfx

	%END OF INNER FUNCTIONS

end%main fx


text:
function get_basemap_parts(resolution)
    % load lake and boundaries
    % resolution is integer from 1 to 5, with 5 being full resolution
    %
    % data from:
    %  https://www.ngdc.noaa.gov/mgg/shorelines/
    %  http://www.soest.hawaii.edu/pwessel/gshhg/index.html
    %
    % this routine expectes them to be in shapefile format
    %
    % GSHHG "Early processing and assmbly of the shoreline data" processing details:
    % Wessel, P., and W. H. F. Smith, A Global Self-consistent, Hierarchical, High-resolution Shoreline Database, J. Geophys. Res., 101, 8741-8743, 1996i
    switch resolution
        case 5
            resolution='f'; % "full"
        case 4
            resolution='h'; % "high"
        case 3
            resolution='i'; % "intermediate"
        case 2
            resolution='l'; % "low"
        otherwise
            resolution='c'; % "crude"
    end

    %% define the map features
    % each MapFeature is something that can be overlain on the main map
    %
    disp('MainInteractiveMap construction');
    obj.Features = MapFeature('coast', @load_coast_and_borders, [],...
        struct('Tag','mainmap_coastline',...
        'DisplayName', 'Coastline/Borders',...
        'HitTest','off','PickableParts','none',...
        'LineWidth',1.0,...
        'Color',[0.1 0.1 0.1])...
        );

    obj.Features(2) = MapFeature('volcanoes', @load_volcanoes, [],...
        struct('Tag','mainmap_volcanoes',...
        'Marker','^',...
        'DisplayName','Volcanoes',...
        'LineWidth', 1.5,...
        'MarkerSize', 6,...
        'LineStyle','none',...
        'MarkerFaceColor','w',...
        'MarkerEdgeColor','r')...
        );

    obj.Features(3) = MapFeature('plates', @load_plates, [],...
        struct('Tag','mainmap_plates',...
        'DisplayName','plate boundaries',...
        'LineWidth', 3.0,...
        'Color',[.2 .2 .5])...
        );

    obj.Features(4) = MapFeature('faults', @load_faults, [],...
        struct('Tag','mainmap_faultlines',...
        'DisplayName','main faultine',...
        'LineWidth', 3.0,...
        'Color','b')...
        );
    %{
            obj.Features(5) = MapFeature('wells', @load_wells, [],...
                struct('Tag','mainmap_wells',...
                    'DisplayName','Wells',...
                    'Marker','d',...
                    'LineWidth',1.5,...
                    'MarkerSize',6,...
                    'LineStyle','none',...
                    'MarkerFaceColor','k',...
                    'MarkerEdgeColor','k')...
                );
            obj.Features(6) = MapFeature('minor_faults', @load_minorfaults, [],...
                struct('Tag','mainmap_faults',...
                    'DisplayName','faults',...
                    'LineWidth',0.2,...
                    'Color','k')...
                );

        %}
        figure;
        %TODO: have this on
        base='/Users/reyesc/Downloads/gshhg-shp-2';
        load_rivers(base,resolution);
        load_lakes(base,resolution);
        %load_boundaries(base,resolution);
        load_faults();
        load_continents(base, resolution);
        load_borders(base, resolution);
        load_plates();
end

function load_borders(base, resolution)
    ZG = ZmapGlobal.Data;
    % load National boundaries
    fn=fullfile(base,'WDBII_shp',resolution,['WDBII_border_' resolution '_L1']);
    S=shaperead(fn,'UseGeoCoords',true);
    ZG.borders=S;
end

function load_continents(base, resolution)
    ZG = ZmapGlobal.Data;
    % load continent boundaries (no antarctica)
    fn=fullfile(base,'GSHHS_shp',resolution,['GSHHS_' resolution '_L1']);
    S=shaperead(fn,'UseGeoCoords',true);
    ZG.coastline=S;
end

function load_faults()
    ZG = ZmapGlobal.Data;
    % load Fault traces
    % from SHARE The European Database of Seismogenic Faults:
    % http://diss.rm.ingv.it/share-edsf/SHARE_WP3.2_Downloads.html
    fn='/Users/reyesc/Downloads/DB-CFS-ESRI-shapefile/Crustal_fault_sources_TOP';
    S=shaperead(fn,'UseGeoCoords',true);
    ZG.faults=S;
end

function load_lakes(base, resolution)
    ZG = ZmapGlobal.Data;
    % load lake boundaries
    fn=fullfile(base,'GSHHS_shp',resolution,['GSHHS_' resolution '_L2']);
    S=shaperead(fn,'UseGeoCoords',true);
    ZG.lakes=S;
end

function load_rivers(base, resolution)
    ZG = ZmapGlobal.Data;
    % load river-lake boundaries
    fn=fullfile(base,'WDBII_shp',resolution,['WDBII_river_' resolution '_L01']);
    S=shaperead(fn,'UseGeoCoords',true);

    % load perm. major river boundaries
    fn=fullfile(base,'WDBII_shp',resolution,['WDBII_river_' resolution '_L02']);
    S=[S ; shaperead(fn,'UseGeoCoords',true)];

    % load addl. major rivers
    fn=fullfile(base,'WDBII_shp',resolution,['WDBII_river_' resolution '_L03']);
    S=[S ; shaperead(fn,'UseGeoCoords',true)];

    % load addl. rivers
    fn=fullfile(base,'WDBII_shp',resolution,['WDBII_river_' resolution '_L04']);
    S=[S ; shaperead(fn,'UseGeoCoords',true)];
    ZG.rivers=S;
end

function load_plates()
    ZG=ZmapGlobal.Data;
    tmp=load('resources/plates.mat','data','metadata');
    disp('Plate Boundary Metadata:');
    disp(tmp.metadata)
    ZG.plates=tmp.data;
end

author: CelsoReyes


text:
function hout=SupTitle(str)
% =======================================================================
% SUPTITLE Puts a title above all subplots
% =======================================================================
%	SUPTITLE('text') adds text to the top of the figure
%	above all subplots (a "super title"). Use this function
%	after all subplot commands.
%
% Drea Thomas 6/15/95 drea@mathworks.com
% Rev. by Guillaume Maze on 2011-02-08: also modify figure title
%
% Warning: If the figure or axis units are non-default, this
% will break.
% =======================================================================

% Rev. by Guillaume Maze on 2011-02-08: To modify figure's name
% Is their already a suptitle:
if ~isempty(findall(gcf,'tag','suptitleText'))
	oldsuptitle = get(findall(gcf,'tag','suptitleText'),'string');
else
	oldsuptitle = '';
end

% Parameters used to position the supertitle.

% Amount of the figure window devoted to subplots
plotregion = .92;

% Y position of title in normalized coordinates
titleypos  = .96;

% Fontsize for supertitle
aux = get(gca,'Title'); aux.FontSize;
fs = aux.FontSize+4;

% Fudge factor to adjust y spacing between subplots
fudge=1;

haold = gca;
figunits = get(gcf,'units');

% Get the (approximate) difference between full height (plot + title
% + xlabel) and bounding rectangle.

	if (~strcmp(figunits,'pixels')),
		set(gcf,'units','pixels');
		pos = get(gcf,'position');
		set(gcf,'units',figunits);
    else
		pos = get(gcf,'position');
	end
	ff = (fs-4)*1.27*5/pos(4)*fudge;

        % The 5 here reflects about 3 characters of height below
        % an axis and 2 above. 1.27 is pixels per point.

% Determine the bounding rectange for all the plots

% h = findobj('Type','axes');

% findobj is a 4.2 thing.. if you don't have 4.2 comment out
% the next line and uncomment the following block.

h  = findobj(gcf,'Type','axes');  % Change suggested by Stacy J. Hills

h2 = findobj(gcf,'tag','footnote');
h  = setdiff(h,h2);

% If you don't have 4.2, use this code instead
%ch = get(gcf,'children');
%h=[];
%for i=1:length(ch),
%  if strcmp(get(ch(i),'type'),'axes'),
%    h=[h,ch(i)];
%  end
%end

max_y=0;
min_y=1;

oldtitle =0;
for i=1:length(h),
	if (~strcmp(get(h(i),'Tag'),'suptitle')),
		pos=get(h(i),'pos');
		if (pos(2) < min_y), min_y=pos(2)-ff/5*3;end;
		if (pos(4)+pos(2) > max_y), max_y=pos(4)+pos(2)+ff/5*2;end;
    else
		oldtitle = h(i);
	end
end

if max_y > plotregion,
	scale = (plotregion-min_y)/(max_y-min_y);
	for i=1:length(h),
		pos = get(h(i),'position');
		pos(2) = (pos(2)-min_y)*scale+min_y;
		pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
		set(h(i),'position',pos);
	end
end

np = get(gcf,'nextplot');
set(gcf,'nextplot','add');
% if (oldtitle),
% 	delete(oldtitle);
% end
ha=axes('pos',[0 1 1 1],'visible','off','Tag','suptitle');
ht=text(.5,titleypos-1,str);set(ht,'horizontalalignment','center','fontsize',fs,'Tag','suptitleText','fontweight','bold','Interpreter','None');
set(gcf,'nextplot',np);
axes(haold);
if nargout,
	hout=ht;
end

% Rev. by Guillaume Maze on 2011-02-08: Also Modify figure title
%stophere
% check if figure's name is already from suptitle string ? :
if ~isempty(oldsuptitle) && strcmp(oldsuptitle,get(gcf,'name'))
	overw = true; % Yes, so we over-write
elseif isempty(get(gcf,'name'))
	overw = true;
else
	overw = false; % No, this is probably a custom name, so we don't change it
end
if overw
	figname = str;
	set(gcf,'name',figname);
end
end%if simple case, parent is a figure handle

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function tests = plotXTest
  tests = functiontests(localfunctions);

  clc
end

%     2d outline, filled case
%     3d outlien, filled case
%     with LS or edgecolor, color options etc.

function teardownOnce(tc)
    close all
end

function plotpoint_test(tc)

    % simple
    points = rand(2,5);
    clf; plot_point(points);
    tc.verifyEqual(length(get(gca, 'Children')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    lines = findobj(gca, 'Type', 'line');
    for i=1:5
        tc.verifyEqual(length(lines(i).XData), 1);
        tc.verifyEqual(length(lines(i).YData), 1);
        tc.verifyEqual(lines(i).LineStyle, 'none');
        tc.verifyEqual(lines(i).MarkerSize, 6);
        tc.verifyEqual(lines(i).MarkerFaceColor, 'none');
        tc.verifyEqual(lines(i).Marker, 'square');
    end

    % markers specified
    clf; plot_point(points, 'rd');
    tc.verifyEqual(length(get(gca, 'Children')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    lines = findobj(gca, 'Type', 'line');
    for i=1:5
        tc.verifyEqual(lines(i).LineStyle, 'none');
        tc.verifyEqual(lines(i).MarkerSize, 6);
        tc.verifyEqual(lines(i).MarkerFaceColor, 'none');
        tc.verifyEqual(lines(i).Marker, 'diamond');
        tc.verifyEqual(lines(i).Color, [1 0 0]);
    end

    % markers specified solid
    clf; plot_point(points, 'rd', 'solid');
    tc.verifyEqual(length(get(gca, 'Children')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    lines = findobj(gca, 'Type', 'line');
    for i=1:5
        tc.verifyEqual(lines(i).LineStyle, 'none');
        tc.verifyEqual(lines(i).MarkerSize, 6);
        tc.verifyEqual(lines(i).MarkerFaceColor, [1 0 0]);
        tc.verifyEqual(lines(i).Marker, 'diamond');
        tc.verifyEqual(lines(i).Color, [1 0 0]);
    end

    % sequential labels
    clf; plot_point(points, 'sequence');
    tc.verifyEqual(length(get(gca, 'Children')), 10);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    set = [1:5];
    for i=1:5
        set = setdiff(set, str2num(labels(i).String));
    end
    tc.verifyEmpty(set, 'Not all labels found');

    % specified labels
    L = {'A', 'B', 'C', 'D', 'E'};
    clf; plot_point(points, 'label', L);
    tc.verifyEqual(length(get(gca, 'Children')), 10);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    set = [1:5];
    for i=1:5
        set = setdiff(set, double(strip(labels(i).String))-'A'+1);
    end
    tc.verifyEmpty(set, 'Not all labels found');

    % printf labels
    clf; plot_point(points, 'printf', {'label=%d', 21:25});
    tc.verifyEqual(length(get(gca, 'Children')), 10);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    set = [21:25];
    for i=1:5
        l = strip(labels(i).String);
        tc.verifyEqual(l(1:6), 'label=');
        set = setdiff(set, str2num(l(7:end)));
    end
    tc.verifyEmpty(set, 'Not all labels found');

    % specify font size
    clf; plot_point(points, 'sequence', 'textsize', 23);
    tc.verifyEqual(length(get(gca, 'Children')), 10);

    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    for i=1:5
        tc.verifyEqual(labels(i).FontSize, 23);
        tc.verifyEqual(labels(i).FontWeight, 'normal');
    end

    % specify font color
    clf; plot_point(points, 'sequence', 'textcolor', 'g');
    tc.verifyEqual(length(get(gca, 'Children')), 10);

    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    for i=1:5
        tc.verifyEqual(labels(i).Color, [0 1 0]);
        tc.verifyEqual(labels(i).FontWeight, 'normal');
    end

    % specify font color
    clf; plot_point(points, 'sequence', 'textcolor', [0.2 0.3 0.4]);
    tc.verifyEqual(length(get(gca, 'Children')), 10);

    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    for i=1:5
        tc.verifyEqual(labels(i).Color, [0.2 0.3 0.4]);
        tc.verifyEqual(labels(i).FontWeight, 'normal');
    end

    % specify font weight
    clf; plot_point(points, 'sequence', 'bold');
    tc.verifyEqual(length(get(gca, 'Children')), 10);

    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    for i=1:5
        tc.verifyEqual(labels(i).FontWeight, 'bold');
    end

    % specify font size, weight, color
    clf; plot_point(points, 'sequence', 'textsize', 23, 'bold', 'textcolor', [0.2 0.3 0.4]);
    tc.verifyEqual(length(get(gca, 'Children')), 10);

    tc.verifyEqual(length(findobj(gca, 'Type', 'line')), 5);
    tc.verifyEqual(length(findobj(gca, 'Type', 'text')), 5);
    labels = findobj(gca, 'Type', 'text');
    for i=1:5
        tc.verifyEqual(labels(i).FontSize, 23);
        tc.verifyEqual(labels(i).FontWeight, 'bold');
        tc.verifyEqual(labels(i).Color, [0.2 0.3 0.4]);
    end
end

function plotpoly_test(tc)

    opt.animate = false;
    % 2d and 3D

    opt.tag = [];
    opt.axis = [];

    P = [1 3 5; 1 5 1];

    clf; plot_poly(P);
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    lines = findobj(gca, 'Type', 'line');
    tc.verifyEqual(length(lines), 1);
    tc.verifyEqual(length(lines(1).XData), 4);
    tc.verifyEqual(length(lines(1).YData), 4);
    tc.verifyEqual(lines(1).Marker, 'none');
    tc.verifyEqual(lines(1).LineStyle, '-');

    clf; plot_poly(P, 'g--');
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    lines = findobj(gca, 'Type', 'line');
    tc.verifyEqual(length(lines), 1);
    tc.verifyEqual(length(lines(1).XData), 4);
    tc.verifyEqual(length(lines(1).YData), 4);
    tc.verifyEqual(lines(1).Marker, 'none');
    tc.verifyEqual(lines(1).LineStyle, '--');
    tc.verifyEqual(lines(1).Color, [0 1 0]);

    clf; plot_poly(P, 'edgecolor', 'r');
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    lines = findobj(gca, 'Type', 'line');
    tc.verifyEqual(length(lines), 1);
    tc.verifyEqual(length(lines(1).XData), 4);
    tc.verifyEqual(length(lines(1).YData), 4);
    tc.verifyEqual(lines(1).Marker, 'none');
    tc.verifyEqual(lines(1).LineStyle, '-');
    %tc.verifyEqual(lines(1).Color, [1 0 0]);

    clf; plot_poly(P, 'edgecolor', [0.2 0.3 0.4]);
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    lines = findobj(gca, 'Type', 'line');
    tc.verifyEqual(length(lines), 1);
    tc.verifyEqual(length(lines(1).XData), 4);
    tc.verifyEqual(length(lines(1).YData), 4);
    tc.verifyEqual(lines(1).Marker, 'none');
    tc.verifyEqual(lines(1).LineStyle, '-');
    %tc.verifyEqual(lines(1).Color, [0.2 0.3 0.4]);

    clf; plot_poly(P, 'tag', 'bob');
    lines = findobj(gca, 'Type', 'line');
    %tc.verifyEqual(lines(1).Tag, 'bob');  % no tag for line type

    %-------- patch mode

    clf; plot_poly(P, 'fillcolor', 'r');
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    patch = findobj(gca, 'Type', 'patch');
    tc.verifyEqual(length(patch), 1);

    tc.verifyEqual(patch.FaceColor, [1 0 0]);
    tc.verifyEqual(patch.FaceAlpha, 1);
    tc.verifyEqual(patch.EdgeColor, [0 0 0]);
    tc.verifyEqual(patch.LineStyle, '-');
    tc.verifySize(patch.Faces, [1 3]);
    tc.verifySize(patch.Vertices, [3 2]);

    clf; plot_poly(P, 'fillcolor', 'g', 'alpha', 0.5);
    tc.verifyEqual(length(get(gca, 'Children')), 1);
    patch = findobj(gca, 'Type', 'patch');
    tc.verifyEqual(length(patch), 1);

    tc.verifyEqual(patch.FaceColor, [0 1 0]);
    tc.verifyEqual(patch.FaceAlpha, 0.5);
    tc.verifyEqual(patch.EdgeColor, [0 0 0]);
    tc.verifyEqual(patch.LineStyle, '-');
    tc.verifySize(patch.Faces, [1 3]);
    tc.verifySize(patch.Vertices, [3 2]);

    clf; plot_poly(P, 'fillcolor', 'b', 'tag', 'bob');
    patch = findobj(gca, 'Type', 'patch');
    %tc.verifyEqual(patch.Tag, 'bob');  % no tag
end

function plot_sphere_test(tc)
end

function plot_box_test(tc)
end

function plot_arrow_test(tc)
end

function plot_homline_test(tc)
end

% plot_ellipse
function ellipse2d_test(tc)
    clf

    E = diag([9 4]);

    plot_ellipse(E)
    plot_ellipse(E, [4 0], 'r--');
    plot_ellipse(E, [0 4], 'edgecolor', 'g');
    plot_ellipse(E, [4 4], 'fillcolor', 'g');
    plot_ellipse(E, [0 8 0.5], 'edgecolor', 'r', 'fillcolor', 'c', 'alpha', 0.5, 'LineWidth', 3);
    plot_ellipse(E, [4 8], 'b--', 'LineWidth', 3);
    plot_ellipse(E, 'confidence', 0.5);

    h = plot_ellipse(E);
    plot_ellipse(E, [1 1], 'alter', h);
    tc.verifyError(@() plot_ellipse(E, 'alter', 7), 'SMTB:plot_ellipse:badarg');
end

function ellipse2d_animate_test(tc)
    clf
    axis([-4 4 -4 4]);

    E = diag([9 4]);

    h = plot_ellipse(E, 'g');
    for x = circle([0 0], 1)
        plot_ellipse(E, x, 'alter', h)
        pause(0.1)
    end

        clf
    axis([-4 4 -4 4]);
    h = plot_ellipse(E, [0 0 0.5], 'edgecolor', 'r', 'fillcolor', 'c', 'LineWidth', 3);

    for x = circle([0 0], 1)
        plot_ellipse(E, x, 'alter', h)
        pause(0.1)
    end

end

% plot_ellipse
function ellipse3d_test(tc)
    clf

    E = diag([9 4 6]);

    plot_ellipse(E)

    clf
    plot_ellipse(E, 'edgecolor', 'g');

    clf
    plot_ellipse(E, 'fillcolor', 'g');

    clf
    plot_ellipse(E, 'fillcolor', 'g', 'shadow');

    clf
    plot_ellipse(E, 'fillcolor', 'g', 'edgecolor', 'r', 'LineWidth', 2);

    plot_ellipse(E, [0 8 1], 'edgecolor', 'r', 'fillcolor', 'c');
    plot_ellipse(E, [4 8 1], 'LineWidth', 3);

    axis equal
end

function ellipse3d_animate_test(tc)
    clf
    axis([-4 4 -4 4 -4 4]);

    E = diag([9 4 6]);

    h = plot_ellipse(E, 'g');
    for x = circle([0 0], 1)
        plot_ellipse(E, [x; 0], 'alter', h)
        pause(0.1)
    end
end

author: petercorke

repo: spatialmath-matlab


text:
function [ad] = PL2AdSpanBySource(filename, source, channelInSource, startCount, endCount)
% PL2AdSpanBySource(filename, source, channelInSource, startCount, endCount): read a span of a/d data from a .pl2 file
%                returns a/d values in millivolts
%
% examples:
%
% ad = PL2AdSpanBySource(filename, source, channelInSource, startCount, endCount)
% ad = PL2AdSpanBySource(filename, 'FP', 1, 10000)
% ad = PL2AdSpanBySource(filename, 8, 1, 1000, 2000)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   source - source name or source numeric id
%   channel -  1-based channel number within source
%           To print the list of analog channels, use these commands:
%
%           pl2 = PL2GetFileIndex(filename);
%           PL2Print(pl2.AnalogChannels);
%
%   startCount - index of first sample to fetch
%   endCount - index of last sample to fetch
%
% OUTPUT:
%   ad.Values - array of a/d values in millivolts
%   ad.ADFreq - digitization frequency for this channel (samples per second)

ad.Values = [];
ad.ADFreq = 0;

if nargin ~= 5
    error 'expected 5 input arguments';
end

filename = internalPL2ResolveFilename(filename);
pl2 = PL2GetFileIndex(filename);
channelNumber = internalPL2ResolveChannelBySource(pl2.AnalogChannels, source, channelInSource);
ad = Pl2AdSpan(filename, channelNumber, startCount, endCount);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function Z = ne(X,Y)
%NE Not equal (~=) for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@ne,X,Y);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
report_this_filefun(mfilename('fullpath'));

av4 = [];
ni = 200
newt3 = newcat;
figure
fi = gcf;
for i = 2.0:0.1:3.1
    l = newt3(:,6) >= i;
    newt2 = newt3(l,:);
    bwithde
    drawnow
    figure_w_normalized_uicontrolunits(fi)
    hold on
    plot(av2(:,2),av2(:,1))
    drawnow

end


text:
% EEG_TOPOPLOT - plot scalp map
%
% eeg_topoplot( vals, chanlocs, 'key', 'val');
%
% Input:
%   vals     - values, one per channel
%   chanlocs - channel structure, same size as vals
%
% Optional inputs:
%   'colormap'   - colormap. Possible colormaps are 'blueredyellow', ...
%                'yellowredblue', 'bluered' or any Matlab colormap ('cool',
%                'jet', 'hsv', ...). It can also be a text file 'xxx.txt'.
%                The text file must contain 3 columns and idealy 64 rows
%                defining the colors in RGB format.
%   'maplimits'  - can be [min max]. This help defines the color scale for
%                maps.
%   'electrodes' - can be 'on' to show electrode dots, 'off', or
%               'labels' to show electrode labels. Default is 'on'
%   'dotsize'    - size of electrode dots. Default is 5.
%   'shading'    - 'flat','interp'  {default: 'interp'}
%   'exclude'    - labels or indices of electrodes not to be plotted. From the
%                compiled files, these must be entered using underscores
%                for separators (e.g., "cz_pz").
%   'sphspline'  - can be 'on' or 'off'. If 'on' spherical splines are used
%                for interpolation of the scalp map. If 'off' standard
%                planar inverse distance interpolation is used.
%   'shrink'     - shrink electrode positions (default is 0.75 to be able to
%                plot electrode at the head limit if spherical interpolation
%                is set and 0.95 for planar 2-D interpolation).
%
% References for spline interpolation:
%   [1] Perrin, F., Pernier, J., Bertrand, O., & Echallier, J. F.
%       (1989). Spherical splines for scalp potential and current
%       density mapping. Electroencephalography and Clinical
%       Neurophysiology, 72, 184-187
%   [2] Ferree, T. C. (2000). Spline Interpolation of the Scalp EEG.
%       Retrieved March 26, 2006, from
%       www.egi.com/Technotes/SplineInterpolation.pdf
%
% limitation: does not plot anything below the upper part of the head

% Copyright (C) Arnaud Delorme, SCCN, INC, 2010
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function eeg_topoplot(values, chanlocs, varargin);

g = [];
for index = 1:2:length(varargin)
    g = setfield(g, varargin{index}, varargin{index+1});
end
if ~isfield(g, 'electrodes'), g.electrodes = 'on'; end
if ~isfield(g, 'colormap'),   g.colormap   = jet;  end
if ~isfield(g, 'maplimits'),  g.maplimits  = [];   end
if ~isfield(g, 'headrad'),    g.headrad    = [];   end
if ~isfield(g, 'sphspline'),  g.sphspline  = 'on'; end
if ~isfield(g, 'shading'),    g.shading    = 'interp'; end
if ~isfield(g, 'contour'),    g.contour    = 'on'; end
if ~isfield(g, 'dotsize'),    g.dotsize    = 5;    end
if ~isfield(g, 'mark'),       g.mark       = [];   end
if ~isfield(g, 'exclude'),    g.exclude    = [];   end
if ~isfield(g, 'linewidth'),  g.linewidth  = 2;    end
if ~isfield(g, 'shrink'),     g.shrink     = 1;    end
if ischar(g.dotsize), g.dotsize = str2num(g.dotsize); end
if any(values == 0)
    inds = find(values == 0);
    if ~isempty( [ chanlocs(inds).theta ])
        g.contour = 'off';
        g.sphspline = 'off';
    end
end

% exclude electrodes
% ------------------
if ~isempty(g.exclude)
    chanlocs(g.exclude) = [];
    values(g.exclude)   = [];
end

% find channel coordinates
% ------------------------
emptyvals = cellfun('isempty', { chanlocs.theta });
th = [ chanlocs.theta ];
rd = [ chanlocs.radius ];
[y x] = pol2cart(th/180*pi, rd); x=-x;
x = x*g.shrink;
y = y*g.shrink;
newvalues            = values;
newvalues(emptyvals) = [];
labls = { chanlocs.labels };
labls(emptyvals) = [];

if strcmpi(g.sphspline, 'on')
    % default head radius
    % -------------------
    g.headrad = 0.5;

    % spherical plotting
    % ------------------
    xelec = [ chanlocs.X ];
	yelec = [ chanlocs.Y ];
	zelec = [ chanlocs.Z ];

    dist = sqrt(xelec.^2+yelec.^2+zelec.^2);
	xelec = xelec./dist;
	yelec = yelec./dist;
	zelec = zelec./dist;

    if g.shrink ~= 1
        [th phi rad] = cart2sph(xelec, yelec, zelec);
        phi = (phi-pi/2)*g.shrink+pi/2;
        [xelec, yelec, zelec] = sph2cart(th, phi, rad);
    end;

	[xsph, ysph, zsph, valsph] = spheric_spline(xelec,yelec,zelec,newvalues);
    surf(-ysph/2,xsph/2,zsph/2,double(valsph), 'edgecolor', 'none'); view([0 0 1]);hold on;
    shading(g.shading);
    top = max(abs(valsph(:)))*1000;

    if strcmpi(g.contour, 'on')
    	[c h] = contour3(-ysph/2, xsph/2, valsph+top/10, 5); view([0 0 1]);
        set(h, 'cdata', [], 'edgecolor', 'k')
    end

	% coordinates for electrodes
	% --------------------------
    xelec(find(zelec < 0)) = [];
    yelec(find(zelec < 0)) = [];
    x = yelec/2;
    y = xelec/2;

else
    % default head radius
    % -------------------
    if isempty(g.headrad);
        g.headrad = max(sqrt(x.^2+y.^2));
    end

    % data points for 2-D data plot
    % -----------------------------
    pnts = linspace(0,2*pi,200/0.25*(g.headrad.^2));
    xx = sin(pnts)*g.headrad;
    yy = cos(pnts)*g.headrad;

	% make grid and add circle
	% ------------------------
    gridres = 30;
	coords = linspace(-g.headrad, g.headrad, gridres);
	ay = repmat(coords,  [gridres 1]);
	ax = repmat(coords', [1 gridres]);
	for ind=1:length(xx)
        [tmp closex] = min(abs(xx(ind)-coords));
        [tmp closey] = min(abs(yy(ind)-coords));
        ax(closex,closey) = xx(ind);
        ay(closex,closey) = yy(ind);
	end
	xx2 = sin(pnts)*(g.headrad-0.01);
	yy2 = cos(pnts)*(g.headrad-0.01);
	for ind=1:length(xx)
        [tmp closex] = min(abs(xx2(ind)-coords));
        [tmp closey] = min(abs(yy2(ind)-coords));
        ax(closex,closey) = xx(ind);
        ay(closex,closey) = yy(ind);
	end

	% linear interpolation and removal of values outside circle
	% ---------------------------------------------------------
    a = griddata(x, y, newvalues, -ay, ax, 'v4');
	aradius = sqrt(ax.^2 + ay.^2);
	indoutcircle = find(aradius(:) > g.headrad+0.01);
	a(indoutcircle) = NaN;
	surf(ay, ax, a, 'edgecolor', 'none'); view([0 0 1]); hold on;
    shading(g.shading);
    top = max(values)*1.5;

	% plot level lines
	% ----------------
    if strcmpi(g.contour, 'on')
        [c h] = contour3(ay, ax, a, 5);
        set(h, 'cdata', [], 'edgecolor', 'k')
    end
end

% plot electrodes as dots
% -----------------------
if strcmpi(g.electrodes, 'on') || strcmpi(g.electrodes, 'labels')
    rad = sqrt(x.^2 + y.^2);
    x(find(rad > g.headrad)) = [];
    y(find(rad > g.headrad)) = [];
    plot3( -x, y, ones(size(x))*top, 'k.', 'markersize', g.dotsize);
    for i = g.mark,      plot3( -x(i), y(i), double(top), 'y.', 'markersize', 4*g.dotsize); plot3( -x(i), y(i), double(top), 'r.', 'markersize', 2*g.dotsize); end
    if strcmpi(g.electrodes, 'labels')
        for index = 1:length(x)
            text( -x(index)+0.02, y(index), double(top), labls{index});
        end
    end
else
    % invisible electrode that avoid plotting problem (no surface, only
    % contours)
    plot3( -x, y, -ones(size(x))*top, 'k.', 'markersize', 0.001);
end

% plot dipoles if any
% -------------------
if ~isempty(g.dipole)
    hold on;
    for index = 1:size(g.dipole,1)
        g.dipole(index,:)   = g.dipole(index,:)*0.5;
        g.dipole(index,3:5) = g.dipole(index,3:5)/norm(g.dipole(index,3:end))*0.2;
        if ~any(g.dipole(index,:))
            fprintf('Note: dipole contains 0 - not plotted\n')
        elseif sum(g.dipole(index,3:4).^2) <= 0.00001
            fprintf('Note: dipole is length 0 - not plotted\n')
        elseif sum(g.dipole(index,1:2).^2) > g.headrad
            fprintf('Note: dipole is outside plotting area - not plotted\n')
        else
            hh = plot3( -g.dipole(index, 2), g.dipole(index, 1), top, '.');
            set(hh, 'color', 'k', 'markersize', 30);
            hh = line( -[g.dipole(index, 2) g.dipole(index, 2)+g.dipole(index, 4)]', ...
                [g.dipole(index, 1) g.dipole(index, 1)+g.dipole(index, 3)]',[top top]);
            set(hh, 'color', 'k', 'linewidth', 30/7);
        end
    end
end

% special colormaps
% -----------------
if ischar(g.colormap)
    if ~isempty(strmatch(g.colormap, { 'hsv' 'jet' 'gray' 'hot' 'cool' 'bone' ...
            'copper', 'pink' 'flag' 'prism' }, 'exact'))
    else % read text file
        g.colormap = load('-ascii', g.colormap);
    end
end;
colormap(g.colormap);

if ~isempty(g.maplimits)
    if ~ischar(g.maplimits) && ~isempty(g.maplimits) && ~isnan(g.maplimits(1))
        caxis(g.maplimits);
    end
end

% main circle
% -----------
radiuscircle = 0.5;
pnts   = linspace(0,2*pi,200);
xc     = sin(pnts)*radiuscircle;
yc     = cos(pnts)*radiuscircle;
sf     = 1; % scaling factor
plot3(xc*sf,yc*sf,ones(size(xc))*top, 'k', 'linewidth', g.linewidth); hold on;

% ears & nose
% -----------
rmax  = 0.5;
base  = rmax-.0046;
basex = 0.18*rmax;                   % nose width
tip   = 1.15*rmax;
tiphw = .04*rmax;                    % nose tip half width
tipr  = .01*rmax;                    % nose tip rounding
q = .04; % ear lengthening
EarX  = [.497-.005  .510  .518  .5299 .5419  .54    .547   .532   .510   .489-.005]; % rmax = 0.5
EarY  = [q+.0555 q+.0775 q+.0783 q+.0746 q+.0555 -.0055 -.0932 -.1313 -.1384 -.1199];

plot3(EarX*sf,EarY*sf,ones(size(EarX))*top,'color','k','LineWidth',g.linewidth)    % plot left ear
plot3(-EarX*sf,EarY*sf,ones(size(EarY))*top,'color','k','LineWidth',g.linewidth)   % plot right ear
plot3([basex;tiphw;0;-tiphw;-basex]*sf,[base;tip-tipr;tip;tip-tipr;base]*sf,top*ones(size([basex;tiphw;0;-tiphw;-basex])),'color','k','LineWidth',g.linewidth);

% axis limits
% -----------
axis off;
set(gca, 'ydir', 'normal');
axis equal
ylimtmp = max(g.headrad, 0.58);
ylim([-ylimtmp ylimtmp]);

% ----------------
% spherical spline
% ----------------
function [x, y, z, Res] = spheric_spline( xelec, yelec, zelec, values);

SPHERERES = 40;
[x,y,z] = sphere(SPHERERES);
x(1:(length(x)-1)/2,:) = [];
y(1:(length(x)-1)/2,:) = [];
z(1:(length(x)-1)/2,:) = [];

Gelec = computeg(xelec,yelec,zelec,xelec,yelec,zelec);
Gsph  = computeg(x,y,z,xelec,yelec,zelec);

% equations are
% Gelec*C + C0  = Potential (C unknown)
% Sum(c_i) = 0
% so
%             [c_1]
%      *      [c_2]
%             [c_ ]
%    xelec    [c_n]
% [x x x x x]         [potential_1]
% [x x x x x]         [potential_ ]
% [x x x x x]       = [potential_ ]
% [x x x x x]         [potential_4]
% [1 1 1 1 1]         [0]

% compute solution for parameters C
% ---------------------------------
meanvalues = mean(values);
values = values - meanvalues; % make mean zero
C = pinv([Gelec;ones(1,length(Gelec))]) * [values(:);0];

% apply results
% -------------
Res = zeros(1,size(Gsph,1));
for j = 1:size(Gsph,1)
    Res(j) = sum(C .* Gsph(j,:)');
end
Res = Res + meanvalues;
Res = reshape(Res, size(x));

% compute G function
% ------------------
function g = computeg(x,y,z,xelec,yelec,zelec)

unitmat = ones(length(x(:)),length(xelec));
EI = unitmat - ((repmat(x(:),1,length(xelec)) - repmat(xelec,length(x(:)),1)).^2 +...
                (repmat(y(:),1,length(xelec)) - repmat(yelec,length(x(:)),1)).^2 +...
                (repmat(z(:),1,length(xelec)) - repmat(zelec,length(x(:)),1)).^2)/2;

g = zeros(length(x(:)),length(xelec));
m = 4; % 3 is linear, 4 is best according to Perrin's curve
for n = 1:7
    L = legendre(n,EI);
    g = g + ((2*n+1)/(n^m*(n+1)^m))*squeeze(L(1,:,:));
end
g = g/(4*pi);

% find electrode indices
% ----------------------
function allinds = elecind( str, chanlocs, values );

    findmax = 0;
    findmin = 0;
    if ~iscell(str)
         if strmatch(str, 'max', 'exact'), findmax = 1; end
         if strmatch(str, 'min', 'exact'), findmin = 1; end;
         indunderscore = [ 0 find( str == '_' ) length(str)+1 ];
    else indunderscore = [1:length(str)+1];
    end

    % find maximum or minimum
    % -----------------------
    if findmax, [tmp allinds] = max(values); return; end
    if findmin, [tmp allinds] = min(values); return; end

    % find indices for labels
    % -----------------------
    labels = lower({ chanlocs.labels });
    for i = 1:length(indunderscore)-1
        if ~iscell(str)
             tmpstr = str(indunderscore(i)+1:indunderscore(i+1)-1);
        else tmpstr = str{i};
        end
        tmpind = strmatch(lower(tmpstr), labels, 'exact');
        if isempty(tmpind)
            if str2num(tmpstr) > 0
                tmpind = str2num(tmpstr);
            else
                error(sprintf('Could not find channel "%s"', tmpstr));
            end
        end
        allinds(i) = tmpind;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [nodes2, edges2] = grRemoveNodes(nodes, edges, rmNodes)
%GRREMOVENODES Remove several nodes in a graph
%
%   usage:
%   [NODES2 EDGES2] = grRemoveNodes(NODES, EDGES, NODES2REMOVE)
%   remove the nodes with indices NODE2REMOVE from array NODES, and also
%   remove edges containing the nodes NODE2REMOVE.
%
%   Example
%     nodes = [...
%         10 10; 20 10; 30 10; ...
%         10 20; 20 20; 30 20];
%     edges = [...
%         1 2; 1 4; 1 5; ...
%         2 3; 2 5; 2 6; ...
%         3 6; 4 5; 5 6];
%     toRemove = [3 4];
%     [nodes2 edges2] = grRemoveNodes(nodes, edges, toRemove);
%     drawGraph(nodes2, edges2);
%     axis equal; axis([0 40 0 30]);
%
%   See also
%     grRemoveEdges
%

%   -----
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY
%   10/02/2004 doc
%   07/03/2014 rewrite using clearer algorithm

%% edges processing

% remove all edges connected to one of the nodes to remove
edges2 = edges(~any(ismember(edges, rmNodes), 2), :);

% change edges information, due to the node index shift
for i = 1:length(rmNodes)
    inds = edges2 > (rmNodes(i) - i + 1);
    edges2(inds) = edges2(inds) - 1;
end

%% nodes processing

% number of nodes
N   = size(nodes, 1);
NR  = length(rmNodes);
N2  = N-NR;

% allocate memory
nodes2 = zeros(N2, 2);

% process the first node
nodes2(1:rmNodes(1)-1,:) = nodes(1:rmNodes(1)-1,:);

for i = 2:NR
    inds = rmNodes(i-1)+1:rmNodes(i)-1;
    if isempty(inds)
        continue;
    end
    nodes2(inds - i + 1, :) = nodes(inds, :);
end

% process the last node
nodes2(rmNodes(NR)-NR+1:N2, :) = nodes(rmNodes(NR)+1:N, :);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
classdef TestBinaryDescriptorMatcher
    %TestBinaryDescriptorMatcher

    properties (Constant)
        fields = {'queryIdx', 'trainIdx', 'imgIdx', 'distance'};
    end

    methods (Static)
        function test_matcher
            %TODO: something in this test is causing memory corruption and
            %TODO: crashing MATLAB on repeated runs.. most likely an OpenCV bug!
            if true
                error('mexopencv:testskip', 'todo');
            end

            X = randi([0,255], [50,32], 'uint8');
            Y = randi([0,255], [10,32], 'uint8');
            matcher = cv.BinaryDescriptorMatcher();
            matcher.add(X);
            matcher.train();

            m = matcher.match(Y);
            validateattributes(m, {'struct'}, {'vector'});
            assert(all(ismember(TestBinaryDescriptorMatcher.fields, ...
                fieldnames(m))));

            k = 3;
            m = matcher.knnMatch(Y, k);
            validateattributes(m, {'cell'}, ...
                {'vector', 'numel',size(Y,1)});
            cellfun(@(mm) validateattributes(mm, {'struct'}, ...
                {'vector', 'numel',k}), m);
            cellfun(@(mm) assert(all(ismember(...
                TestBinaryDescriptorMatcher.fields, fieldnames(mm)))), m);

            maxDist = 125.0;
            m = matcher.radiusMatch(Y, maxDist);
            validateattributes(m, {'cell'}, ...
                {'vector', 'numel',size(Y,1)});
            cellfun(@(mm) validateattributes(mm, ...
                {'struct'}, {'vector'}), m);
            cellfun(@(mm) assert(all(ismember(...
                TestBinaryDescriptorMatcher.fields, fieldnames(mm)))), m);
            assert(all(cellfun(@(mm) all([mm.distance] <= maxDist), m)));
        end

        function test_matcher_mask
            X = randi([0,255], [50,32], 'uint8');
            Y = randi([0,255], [10,32], 'uint8');
            matcher = cv.BinaryDescriptorMatcher();
            mask = rand(size(Y,1), 1) > 0.5;

            m = matcher.match(Y, X, 'Mask',mask);
            validateattributes(m, {'struct'}, {'vector'});

            m = matcher.knnMatch(Y, X, 3, 'Mask',mask);
            validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});

            m = matcher.radiusMatch(Y, X, 2.0, 'Mask',mask);
            validateattributes(m, {'cell'}, {'vector', 'numel',size(Y,1)});
        end

        function test_closest_image
            %TODO: something in this test is causing memory corruption and
            %TODO: crashing MATLAB on repeated runs.. most likely an OpenCV bug!
            if true
                error('mexopencv:testskip', 'todo');
            end

            % compute descriptors for all training images
            obj = cv.BinaryDescriptor('NumOfOctave',1);
            matcher = cv.BinaryDescriptorMatcher();
            query = [];
            files = dir(fullfile(mexopencv.root(),'test','shape0*.png'));
            if isempty(files)
                error('mexopencv:testskip', 'missing data');
            end
            idxQuery = 1;
            idxTrain = 2:numel(files);
            for i=1:numel(files)
                img = cv.imread(fullfile(mexopencv.root(),'test',files(i).name), 'Flags',0);
                whos img
                [~,desc] = obj.detectAndCompute(img);
                assert(~isempty(desc));
                if i == idxQuery
                    query = desc;
                else
                    matcher.add(desc);
                end
            end

            % find matches between query and train set
            matcher.train();
            m = matcher.match(query);

            % naive approach to find nearest image
            [~,ord] = sort([m.distance]);         % retain best 100 (min dist)
            m = m(ord(1:min(10,end)));
            ii = idxTrain(mode([m.imgIdx]) + 1);  % most frequest image index
            query = files(idxQuery).name;
            closest = files(ii).name;
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function progressbar(varargin)
% Description:
%   progressbar() provides an indication of the progress of some task using
% graphics and text. Calling progressbar repeatedly will update the figure and
% automatically estimate the amount of time remaining.
%   This implementation of progressbar is intended to be extremely simple to use
% while providing a high quality user experience.
%
% Features:
%   - Can add progressbar to existing m-files with a single line of code.
%   - Supports multiple bars in one figure to show progress of nested loops.
%   - Optional labels on bars.
%   - Figure closes automatically when task is complete.
%   - Only one figure can exist so old figures don't clutter the desktop.
%   - Remaining time estimate is accurate even if the figure gets closed.
%   - Minimal execution time. Won't slow down code.
%   - Randomized color. When a programmer gets bored...
%
% Example Function Calls For Single Bar Usage:
%   progressbar               % Initialize/reset
%   progressbar(0)            % Initialize/reset
%   progressbar('Label')      % Initialize/reset and label the bar
%   progressbar(0.5)          % Update
%   progressbar(1)            % Close
%
% Example Function Calls For Multi Bar Usage:
%   progressbar(0, 0)         % Initialize/reset two bars
%   progressbar('A', '')      % Initialize/reset two bars with one label
%   progressbar('', 'B')      % Initialize/reset two bars with one label
%   progressbar('A', 'B')     % Initialize/reset two bars with two labels
%   progressbar(0.3)          % Update 1st bar
%   progressbar(0.3, [])      % Update 1st bar
%   progressbar([], 0.3)      % Update 2nd bar
%   progressbar(0.7, 0.9)     % Update both bars
%   progressbar(1)            % Close
%   progressbar(1, [])        % Close
%   progressbar(1, 0.4)       % Close
%
% Notes:
%   For best results, call progressbar with all zero (or all string) inputs
% before any processing. This sets the proper starting time reference to
% calculate time remaining.
%   Bar color is choosen randomly when the figure is created or reset. Clicking
% the bar will cause a random color change.
%
% Demos:
%     % Single bar
%     m = 500;
%     progressbar % Init single bar
%     for i = 1:m
%       pause(0.01) % Do something important
%       progressbar(i/m) % Update progress bar
%     end
%
%     % Simple multi bar (update one bar at a time)
%     m = 4;
%     n = 3;
%     p = 100;
%     progressbar(0,0,0) % Init 3 bars
%     for i = 1:m
%         progressbar([],0) % Reset 2nd bar
%         for j = 1:n
%             progressbar([],[],0) % Reset 3rd bar
%             for k = 1:p
%                 pause(0.01) % Do something important
%                 progressbar([],[],k/p) % Update 3rd bar
%             end
%             progressbar([],j/n) % Update 2nd bar
%         end
%         progressbar(i/m) % Update 1st bar
%     end
%
%     % Fancy multi bar (use labels and update all bars at once)
%     m = 4;
%     n = 3;
%     p = 100;
%     progressbar('Monte Carlo Trials','Simulation','Component') % Init 3 bars
%     for i = 1:m
%         for j = 1:n
%             for k = 1:p
%                 pause(0.01) % Do something important
%                 % Update all bars
%                 frac3 = k/p;
%                 frac2 = ((j-1) + frac3) / n;
%                 frac1 = ((i-1) + frac2) / m;
%                 progressbar(frac1, frac2, frac3)
%             end
%         end
%     end
%
% Author:
%   Steve Hoelzer
%
% Revisions:
% 2002-Feb-27   Created function
% 2002-Mar-19   Updated title text order
% 2002-Apr-11   Use floor instead of round for percentdone
% 2002-Jun-06   Updated for speed using patch (Thanks to waitbar.m)
% 2002-Jun-19   Choose random patch color when a new figure is created
% 2002-Jun-24   Click on bar or axes to choose new random color
% 2002-Jun-27   Calc time left, reset progress bar when fractiondone == 0
% 2002-Jun-28   Remove extraText var, add position var
% 2002-Jul-18   fractiondone input is optional
% 2002-Jul-19   Allow position to specify screen coordinates
% 2002-Jul-22   Clear vars used in color change callback routine
% 2002-Jul-29   Position input is always specified in pixels
% 2002-Sep-09   Change order of title bar text
% 2003-Jun-13   Change 'min' to 'm' because of built in function 'min'
% 2003-Sep-08   Use callback for changing color instead of string
% 2003-Sep-10   Use persistent vars for speed, modify titlebarstr
% 2003-Sep-25   Correct titlebarstr for 0% case
% 2003-Nov-25   Clear all persistent vars when percentdone = 100
% 2004-Jan-22   Cleaner reset process, don't create figure if percentdone = 100
% 2004-Jan-27   Handle incorrect position input
% 2004-Feb-16   Minimum time interval between updates
% 2004-Apr-01   Cleaner process of enforcing minimum time interval
% 2004-Oct-08   Seperate function for timeleftstr, expand to include days
% 2004-Oct-20   Efficient if-else structure for sec2timestr
% 2006-Sep-11   Width is a multiple of height (don't stretch on widescreens)
% 2010-Sep-21   Major overhaul to support multiple bars and add labels
%

persistent progfig progdata lastupdate

% Get inputs
if nargin > 0
    input = varargin;
    ninput = nargin;
else
    % If no inputs, init with a single bar
    input = {0};
    ninput = 1;
end

% If task completed, close figure and clear vars, then exit
if input{1} == 1
    if ishandle(progfig)
        delete(progfig) % Close progress bar
    end
    clear progfig progdata lastupdate % Clear persistent vars
    drawnow
    return
end

% Init reset flag
resetflag = false;

% Set reset flag if first input is a string
if ischar(input{1})
    resetflag = true;
end

% Set reset flag if all inputs are zero
if input{1} == 0
    % If the quick check above passes, need to check all inputs
    if all([input{:}] == 0) && (length([input{:}]) == ninput)
        resetflag = true;
    end
end

% Set reset flag if more inputs than bars
if ninput > length(progdata)
    resetflag = true;
end

% If reset needed, close figure and forget old data
if resetflag
    if ishandle(progfig)
        delete(progfig) % Close progress bar
    end
    progfig = [];
    progdata = []; % Forget obsolete data
end

% Create new progress bar if needed
if ishandle(progfig)
else % This strange if-else works when progfig is empty (~ishandle() does not)

    % Define figure size and axes padding for the single bar case
    height = 0.03;
    width = height * 8;
    hpad = 0.02;
    vpad = 0.25;

    % Figure out how many bars to draw
    nbars = max(ninput, length(progdata));

    % Adjust figure size and axes padding for number of bars
    heightfactor = (1 - vpad) * nbars + vpad;
    height = height * heightfactor;
    vpad = vpad / heightfactor;

    % Initialize progress bar figure
    left = (1 - width) / 2;
    bottom = (1 - height) / 2;
    progfig = figure(...
        'Units', 'normalized',...
        'Position', [left bottom width height],...
        'NumberTitle', 'off',...
        'Resize', 'off',...
        'MenuBar', 'none' );

    % Initialize axes, patch, and text for each bar
    left = hpad;
    width = 1 - 2*hpad;
    vpadtotal = vpad * (nbars + 1);
    height = (1 - vpadtotal) / nbars;
    for ndx = 1:nbars
        % Create axes, patch, and text
        bottom = vpad + (vpad + height) * (nbars - ndx);
        progdata(ndx).progaxes = axes( ...
            'Position', [left bottom width height], ...
            'XLim', [0 1], ...
            'YLim', [0 1], ...
            'Box', 'on', ...
            'ytick', [], ...
            'xtick', [] );
        progdata(ndx).progpatch = patch( ...
            'XData', [0 0 0 0], ...
            'YData', [0 0 1 1] );
        progdata(ndx).progtext = text(0.99, 0.5, '', ...
            'HorizontalAlignment', 'Right', ...
            'FontUnits', 'Normalized', ...
            'FontSize', 0.7 );
        progdata(ndx).proglabel = text(0.01, 0.5, '', ...
            'HorizontalAlignment', 'Left', ...
            'FontUnits', 'Normalized', ...
            'FontSize', 0.7 );
        if ischar(input{ndx})
            set(progdata(ndx).proglabel, 'String', input{ndx})
            input{ndx} = 0;
        end

        % Set callbacks to change color on mouse click
        set(progdata(ndx).progaxes, 'ButtonDownFcn', {@changecolor, progdata(ndx).progpatch})
        set(progdata(ndx).progpatch, 'ButtonDownFcn', {@changecolor, progdata(ndx).progpatch})
        set(progdata(ndx).progtext, 'ButtonDownFcn', {@changecolor, progdata(ndx).progpatch})
        set(progdata(ndx).proglabel, 'ButtonDownFcn', {@changecolor, progdata(ndx).progpatch})

        % Pick a random color for this patch
        changecolor([], [], progdata(ndx).progpatch)

        % Set starting time reference
        if ~isfield(progdata(ndx), 'starttime') || isempty(progdata(ndx).starttime)
            progdata(ndx).starttime = clock;
        end
    end

    % Set time of last update to ensure a redraw
    lastupdate = clock - 1;

end

% Process inputs and update state of progdata
for ndx = 1:ninput
    if ~isempty(input{ndx})
        progdata(ndx).fractiondone = input{ndx};
        progdata(ndx).clock = clock;
    end
end

% Enforce a minimum time interval between graphics updates
myclock = clock;
if abs(myclock(6) - lastupdate(6)) < 0.01 % Could use etime() but this is faster
    return
end

% Update progress patch
for ndx = 1:length(progdata)
    set(progdata(ndx).progpatch, 'XData', ...
        [0, progdata(ndx).fractiondone, progdata(ndx).fractiondone, 0])
end

% Update progress text if there is more than one bar
if length(progdata) > 1
    for ndx = 1:length(progdata)
        set(progdata(ndx).progtext, 'String', ...
            sprintf('%1d%%', floor(100*progdata(ndx).fractiondone)))
    end
end

% Update progress figure title bar
if progdata(1).fractiondone > 0
    runtime = etime(progdata(1).clock, progdata(1).starttime);
    timeleft = runtime / progdata(1).fractiondone - runtime;
    timeleftstr = sec2timestr(timeleft);
    titlebarstr = sprintf('%2d%%    %s remaining', ...
        floor(100*progdata(1).fractiondone), timeleftstr);
else
    titlebarstr = ' 0%';
end
set(progfig, 'Name', titlebarstr)

% Force redraw to show changes
drawnow

% Record time of this update
lastupdate = clock;

% ------------------------------------------------------------------------------
function changecolor(h, e, progpatch) %#ok<INUSL>
% Change the color of the progress bar patch

% Prevent color from being too dark or too light
colormin = 1.5;
colormax = 2.8;

thiscolor = rand(1, 3);
while (sum(thiscolor) < colormin) || (sum(thiscolor) > colormax)
    thiscolor = rand(1, 3);
end

set(progpatch, 'FaceColor', thiscolor)

% ------------------------------------------------------------------------------
function timestr = sec2timestr(sec)
% Convert a time measurement from seconds into a human readable string.

% Convert seconds to other units
w = floor(sec/604800); % Weeks
sec = sec - w*604800;
d = floor(sec/86400); % Days
sec = sec - d*86400;
h = floor(sec/3600); % Hours
sec = sec - h*3600;
m = floor(sec/60); % Minutes
sec = sec - m*60;
s = floor(sec); % Seconds

% Create time string
if w > 0
    if w > 9
        timestr = sprintf('%d week', w);
    else
        timestr = sprintf('%d week, %d day', w, d);
    end
elseif d > 0
    if d > 9
        timestr = sprintf('%d day', d);
    else
        timestr = sprintf('%d day, %d hr', d, h);
    end
elseif h > 0
    if h > 9
        timestr = sprintf('%d hr', h);
    else
        timestr = sprintf('%d hr, %d min', h, m);
    end
elseif m > 0
    if m > 9
        timestr = sprintf('%d min', m);
    else
        timestr = sprintf('%d min, %d sec', m, s);
    end
else
    timestr = sprintf('%d sec', s);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [y] = l1median_VaZh_z(X, maxiter, tol, zerotol, medIn )
%To calculate the L1_median of X  with the l1median_VaZh method
%INput:
%     X       ---the number of the records is the number of rows
%     maxiter ---the max number of iteration
%     tol     ---the tolerance of iteration
%     zerotol ---the tolerance of vector to zero
%     medIn   ---a row vector representing initial median of X
%output:
%     y       ---the median of X

    %To discuss the parameters
    if (nargin > 5)
        error ('Too many input arguments.') ;
    elseif (nargin < 5)
        medIn = median(X);
        if(nargin < 4)
            zerotol = 1e-15 ;
            if (nargin < 3)
                tol = 1e-9 ;
                if (nargin < 2)
                    maxiter = 200 ;
                    if (nargin < 1)
                        error ('Data matrix X is missing.') ;
                    end
                end
            end
        end
    end

    [rn,vn] = size(X);
    iterdis = 1;
    iter = 0;
    y = medIn;

    %Begin the iteration
    while (iterdis > 0) && (iter < maxiter)
        Tnum=zeros(1,vn);
        R = zeros(1,vn);
        Tden=0;
        yita = 0;
        for i=1:rn
            dist = norm(X(i,:)-y);
            if dist >= zerotol
                Tnum = Tnum + X(i,:)/dist;
                Tden=Tden + 1/dist;
                R = R + (X(i,:)-y)/dist;
            else
                yita = 1;
            end
        end

        if Tden == 0
            T=0;
        else
           T=Tnum/Tden;
        end

        if norm(R) == 0
            r = 0;
        else
            r = min(1,yita/norm(R));
        end

        Ty = (1-r)*T + r*y;

        iterdis = norm((Ty-y),1) - tol*norm(y,1);
        iter = iter + 1;
        y=Ty;

    end

    %End the iteration

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS