text:
%kmeigen 'Compute eigenvectors and eigenvalues of a square matrix '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros meigen.pane file
%
% Parameters:
% InputFile: i 'Input Matrix', required: 'Input Matrix'
% OutputFile: o1 'Eigenvectors', optional: 'Eigenvector output matrix'
% OutputFile: o2 'Eigenvalues ', optional: 'Eigenvalue output matrix'
%
% Example: [o1, o2] = kmeigen(i, {'i','';'o1','';'o2',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% meigen - Compute Eigenvectors and Eigenvalues of a Square Matrix
%
%  DESCRIPTION
% .I meigen
% computes the eigenvalues and right eigenvectors of a square matrix.
% The input matrix can be nonsymmetric.
% The input matrix is caste to KDCOMPLEX for reading, and all output
% matrices are of type KDCOMPLEX.
%
% The eigenvectors are stored as columns of the eigenvector output matrix.
% The eigenvectors are all normalized to have unit 2-norm and a purely real
% largest component.
%
% The eigenvalues are stored as a diagonal matrix with the i'th diagonal
% being the eigenvalue corresponding the the eigenvector in column i of the
% eigenvector output matrix.
%
% The eigenvalues (and thier associated eigenvectors) appear in decreasing
% order. They are explicitly sorted (while keeping track of the associated
% eigenvectors).
%
% "meigen\fR uses the LAPACK routine ZGEEV to perform the calculations.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% MATRIX:zlapack:zgeev
%
%  RESTRICTIONS
% It has been observed that
% LAPACK generally gets the eigenvalues sorted in decending order, but
% occasionally will return with some of the smaller eigenvalues out of
% order. This is not a bug in "meigen\fR, although a sort may be inserted
% at some future date to make sure that the eigenvalues and thier associated
% eigenvectors are sorted. The 1992 version of the LAPACK USERS' GUIDE doesn't
% specify an ordering for the eigenvalues.
%
%  REFERENCES
% LAPACK Users' Guide, E. Anderson et. al., SIAM 1992. ISBN 0-89871-294-7
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function pot_type = determine_pot_type(model, onodes, nodes)
% DETERMINE_POT_TYPE Determine the type of potential based on the evidence pattern.
% pot_type = determine_pot_type(model, onodes, nodes)
%
% If there are any utility nodes, pot_type = 'u'
% else
% If all hidden nodes are discrete, pot_type = 'd'.
% If all hidden nodes are continuous, pot_type = 'g' (Gaussian).
% If some hidden nodes are discrete, and some cts, pot_type = 'cg' (conditional Gaussian).
%
% nodes defaults to all nodes in graph

nnodes = length(model.node_sizes);
if nargin < 3, nodes = 1:nnodes; end

hnodes = mysetdiff(nodes, onodes);
if isfield(model, 'limid') %~isempty(model.utility_nodes)
  pot_type = 'u';
elseif isempty(myintersect(model.cnodes, hnodes))
  pot_type = 'd';
elseif mysubset(hnodes, model.cnodes)
  pot_type = 'g';
else
  pot_type = 'cg';
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function demo_fast_sc(opt_choice)
% opt_choice = 1: use epslion-L1 penalty
% opt_choice = 2: use L1 penalty

if ~exist('opt_choice', 'var')
    opt_choice = 1;
end

% natural image data
load ../data/IMAGES.mat
X = getdata_imagearray(IMAGES, 14, 10000);

% sparse coding parameters
num_bases = 128;
beta = 0.4;
batch_size = 1000;
num_iters = 100;
if opt_choice==1
    sparsity_func= 'epsL1';
    epsilon = 0.01;
elseif opt_choice==2
    sparsity_func= 'L1';
    epsilon = [];
end

Binit = [];
fname_save = sprintf('../results/sc_%s_b%d_beta%g_%s', sparsity_func, num_bases, beta, datestr(now, 30));

% run fast sparse coding
[B S stat] = sparse_coding(X, num_bases, beta, sparsity_func, epsilon, num_iters, batch_size, fname_save, Binit);

author: thomas-koehler

repo: SupER


text:
function obj = modelVarsFromFilesCheck(obj,folderNo,varargin)
	% Checks whether a specific folder has variables loaded, if not loads them.
	% Biafra Ahanonu
	% started: 2021.06.18 [20:25:49]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.06.21 [21:03:42] - Add check for objLocations
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% ========================
	% DESCRIPTION
	options.baseOption = '';
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	% ========================

	try
		try
			max(1,round(obj.objLocations{folderNo}.(obj.signalExtractionMethod)(:,1)));
			[rawSignals rawImages signalPeaks signalPeaksArray, ~, ~, rawSignals2] = modelGetSignalsImages(obj,'returnType','raw');
			% obj.nSignals{fileNum}
			skipReload = 1;
		catch
			obj.guiEnabled = 0;
			% originalFileNum =
            oldTmp = obj.foldersToAnalyze;
            obj.foldersToAnalyze = folderNo;
			obj.modelVarsFromFiles();
            obj.foldersToAnalyze = oldTmp;
			obj.fileNum = folderNo;
			obj.guiEnabled = 1;
			skipReload = 0;
		end
	catch err
		display(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		display(repmat('@',1,7))
	end

end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function degree = r8poly_degree ( na, a )

%*****************************************************************************80
%
%% R8POLY_DEGREE returns the degree of a polynomial in power sum form.
%
%  Discussion:
%
%    The power sum form of a polynomial is:
%
%      p(x) = a(0) + a(1) * x + ... + a(n-1) * x**(n-1) + a(n) * x**(n)
%
%    The degree of a polynomial is the index of the highest power
%    of X with a nonzero coefficient.
%
%    The degree of a constant polynomial is 0.  The degree of the
%    zero polynomial is debatable, but this routine returns the
%    degree as 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NA, the dimension of A.
%
%    Input, real A(1:NA+1), the coefficients of the polynomials.
%
%    Output, integer DEGREE, the degree of A.
%
  degree = na;

  while ( 0 < degree )

    if ( a(degree+1) ~= 0.0 )
      return
    end

    degree = degree - 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [outputStruct]=febioStructTemplate(varargin)

% function [outputStruct]=febioStructTemplate(inputStruct)
% ------------------------------------------------------------------------
% This function creates a template for the febio_spec structure for FEBio.
% This structure is intended to be converted to an XML file using
% |febioStruct2xml|.
%
%
% Change log:
% 2020/12/01: KMM Changed to use febio_spec 3.0 by default
% 2021/04/29: KMM Removed strain energy output as this is not available for
% all materials (e.g. viscoelastic materials).
% 2023/03/06: KMM Updating for FEBio spec 4.0.
%
% To do:
% * Check if still plot_stride is still needed
% * Check if both 0 and 1 and "symemtric" and presumably
% "assymetric" are options for the solver
% * No documentation for check_zero_diagonal
%

% ------------------------------------------------------------------------

%%

switch nargin
    case 0
        inputStruct=[];
    case 1
        inputStruct=varargin{1};
    otherwise
        error('Wrong number of input arguments');
end

%% Check for old FEBio version (febio2)

if contains(lower(getFEBioPath),'febio2')
    [outputStruct]=febioStructTemplate_v2p5(inputStruct);
elseif contains(lower(getFEBioPath),'febio3')
    [outputStruct]=febioStructTemplate_v3p0(inputStruct);
    return
end

%% Set febio_spec version
febio_spec.ATTR.version='4.0';

%% Module section
febio_spec.Module.ATTR.type='solid'; %Use default set

%% Control section

febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=10;
febio_spec.Control.step_size=0.1;

febio_spec.Control.plot_zero_state=1;
febio_spec.Control.plot_range=[0,-1];
febio_spec.Control.plot_level='PLOT_MAJOR_ITRS';
febio_spec.Control.plot_stride=1;
febio_spec.Control.output_level='OUTPUT_MAJOR_ITRS';

febio_spec.Control.adaptor_re_solve=1;

febio_spec.Control.time_stepper.ATTR.type='default';
febio_spec.Control.time_stepper.max_retries=5;
febio_spec.Control.time_stepper.opt_iter=6;
febio_spec.Control.time_stepper.dtmin=0.01;
febio_spec.Control.time_stepper.dtmax=0.1;
febio_spec.Control.time_stepper.aggressiveness=0;
febio_spec.Control.time_stepper.cutback=0.5;
febio_spec.Control.time_stepper.dtforce=0;

febio_spec.Control.solver.ATTR.type='solid';
febio_spec.Control.solver.symmetric_stiffness=1;%'symmetric';
febio_spec.Control.solver.equation_scheme='staggered';
febio_spec.Control.solver.equation_order='default';
febio_spec.Control.solver.optimize_bw=0;
febio_spec.Control.solver.lstol=0.9;

febio_spec.Control.solver.lsmin=0.01;
febio_spec.Control.solver.lsiter=5;

febio_spec.Control.solver.max_refs=25;
febio_spec.Control.solver.check_zero_diagonal=0;
febio_spec.Control.solver.zero_diagonal_tol=0;
febio_spec.Control.solver.force_partition=0;
febio_spec.Control.solver.reform_each_time_step=1;
febio_spec.Control.solver.reform_augment=0;
febio_spec.Control.solver.diverge_reform=1;
febio_spec.Control.solver.min_residual=1e-20;
febio_spec.Control.solver.max_residual=0;
febio_spec.Control.solver.dtol=0.001;
febio_spec.Control.solver.etol=0.01;
febio_spec.Control.solver.rtol=0;
febio_spec.Control.solver.rhoi=0;
febio_spec.Control.solver.alpha=1;
febio_spec.Control.solver.beta=0.25;
febio_spec.Control.solver.gamma=0.5;
febio_spec.Control.solver.logSolve=0;
febio_spec.Control.solver.arc_length=0;
febio_spec.Control.solver.arc_length_scale=0;

febio_spec.Control.solver.qn_method.ATTR.type='BFGS';
febio_spec.Control.solver.qn_method.max_ups=0; %Default 10
febio_spec.Control.solver.qn_method.max_buffer_size=0;
febio_spec.Control.solver.qn_method.cycle_buffer=1;
febio_spec.Control.solver.qn_method.cmax=1e5;

%% Globals section

febio_spec.Globals.Constants.R=8.314e-6;
febio_spec.Globals.Constants.T=298;
febio_spec.Globals.Constants.Fc=96485e-9;

%% Material section

%% LoadData section
% febio_spec.LoadData.load_controller{1}.ATTR.id=1;
% febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
% febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
% febio_spec.LoadData.load_controller{1}.extend='CONSTANT';
% febio_spec.LoadData.load_controller{1}.points.pt.VAL=[0 0; 1 1];

%% Output section

% Plot file
febio_spec.Output.plotfile.ATTR.type='febio';
febio_spec.Output.plotfile.var{1}.ATTR.type='displacement';
febio_spec.Output.plotfile.var{2}.ATTR.type='stress';
febio_spec.Output.plotfile.var{3}.ATTR.type='relative volume';
febio_spec.Output.plotfile.var{4}.ATTR.type='reaction forces';
febio_spec.Output.plotfile.var{5}.ATTR.type='contact pressure';
febio_spec.Output.plotfile.compression=0;

%% Fill in missing if input structure is provided

switch nargin
    case 0 %Make output the full default
        outputStruct=febio_spec;
    case 1 %Make output by complementing input with default
        [outputStruct]=structComplete(inputStruct,febio_spec,1); %Complement provided with default if missing
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function X = Adj_DetailCurveCoeff(C,is_real);
% Adj_DetailCurveCoeff: Adjoint of DetailCurveCoeff
%  Usage:
%    X = Adj_DetailCurveCoeff(C);
%  Inputs:
%    C    matrix of curvelet coefficients at scale 2^j
%  Outputs:
%    X    matrix of Fourier samples; jth dyadic subband
%  See Also
%   DetailCurveCoeff, Adj_SeparateAngles, Adj_Curvelet02Xform
%
% By Emmanuel Candes, 2003-2004

  C = ClockwisetoWENS(C);
  nn = size(C);
  R = zeros(nn);
  deep = log2(nn(2));

  for j = 1:size(R,1),
    for m = 1:size(R,2),
      W = squeeze(C(j,m,:,:));
      W = fft2_mid0(W)/sqrt(prod(size(W)));
      R(j,m,:,:) = W;
    end
  end

  for w=1:size(R,2)
    tmp = squeeze(R(2,w,:,:));
    R(2,w,:,:) = tmp([2:end,1], [2:end,1]);
  end
  for w=1:size(R,2)
    tmp = squeeze(R(4,w,:,:));
    R(4,w,:,:) = tmp([2:end,1], [2:end,1]);
  end

  X = Adj_SeparateAngles(Adj_SqueezeAngularFT(R),deep,is_real);

author: thomaskuestner


text:
classdef TestDrawMarkerAruco
    %TestDrawMarkerAruco

    methods (Static)
        function test_1
            img = cv.drawMarkerAruco('6x6_50', 0, 60);
            validateattributes(img, {'uint8'}, {'size',[60 60]});
        end

        function test_2
            img = cv.drawMarkerAruco('6x6_50', 0, 60, 'BorderBits',1);
            validateattributes(img, {'uint8'}, {'size',[60 60]});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [mps] = mach2mps(mach)
% Convert speed from mach (at standard temperature and pressure!) to
% meters per second.
% Chad A. Greene 2012
mps = mach*343;

author: Sable

repo: mcbench-benchmarks


text:
function net = glminit(net, prior)
%GLMINIT Initialise the weights in a generalized linear model.
%
%	Description
%
%	NET = GLMINIT(NET, PRIOR) takes a generalized linear model NET and
%	sets the weights and biases by sampling from a Gaussian distribution.
%	If PRIOR is a scalar, then all of the parameters (weights and biases)
%	are sampled from a single isotropic Gaussian with inverse variance
%	equal to PRIOR. If PRIOR is a data structure similar to that in
%	MLPPRIOR but for a single layer of weights, then the parameters are
%	sampled from multiple Gaussians according to their groupings (defined
%	by the INDEX field) with corresponding variances (defined by the
%	ALPHA field).
%
%	See also
%	GLM, GLMPAK, GLMUNPAK, MLPINIT, MLPPRIOR
%

%	Copyright (c) Ian T Nabney (1996-2001)

errstring = consist(net, 'glm');
if ~isempty(errstring);
  error(errstring);
end
if isstruct(prior)
  sig = 1./sqrt(prior.index*prior.alpha);
  w = sig'.*randn(1, net.nwts);
elseif size(prior) == [1 1]
  w = randn(1, net.nwts).*sqrt(1/prior);
else
  error('prior must be a scalar or a structure');
end

net = glmunpak(net, w);


text:
classdef dme_line3p_opf < mp.dme_line3p & mp.dme_shared_opf

%   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     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% POP_RMBASE - remove channel baseline means from an epoched or
%                continuous EEG dataset. Calls RMBASE.
% Usage:
%   >> OUTEEG = pop_rmbase( EEG ); % pop up an interactive arg entry window
%   >> OUTEEG = pop_rmbase( EEG, timerange, pointrange, chanlist); % call RMBASE
%
% Graphic interface:
%    "Baseline latency range" - [edit box] Latency range for the baseline in ms.
%                               Collects the 'timerange' command line input.
%                               Empty or [] input -> Use whole epoch as baseline
%    "Baseline points vector" - [edit box] Collects the 'pointrange' command line
%                               option (below). (Overwritten by 'timerange' above).
%                               Empty or [] input -> Use whole epoch as baseline
% Inputs:
%   EEG        - Input dataset
%   timerange  - [min_ms max_ms] Baseline latency range in milliseconds.
%                                Empty or [] input -> Use whole epoch as baseline
%   pointrange - [min:max]       Baseline points vector (overwritten by timerange).
%                                Empty or [] input -> Use whole epoch as baseline
%   chanlist   - [cell]          List of channels. Default is all.
%
% Outputs:
%   OUTEEG     - Output dataset
%
% Note: If dataset is continuous, channel means are removed separately
%       for each continuous data region, respecting 'boundary' events
%       marking boundaries of excised or concatenated data portions.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: RMBASE, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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.

% 01-25-02 reformated help & license -ad

function [EEG, com] = pop_rmbase( EEG, timerange, pointrange, chanlist)

com ='';
if nargin < 1
    help pop_rmbase;
    return;
end
if isempty(EEG(1).data)
    disp('pop_rmbase(): cannot remove baseline of an empty dataset'); return;
end
if nargin < 1
	help pop_rmbase;
	return;
end
if nargin < 4 || isempty(chanlist)
    chanlist = 1:EEG(1).nbchan;
end
if nargin < 2, timerange = [];  end
if nargin < 3, pointrange = []; end
if nargin < 2
    % popup window parameters
    % -----------------------
    defaultbase = [num2str(EEG(1).times(1)) ' 0'];
    if EEG(1).times(1) >= 0
        defaultbase = '[ ]';
    end
    if EEG(1).trials > 1
        uilist = { { 'style' 'text' 'string' 'Baseline latency range ([min max] in ms) ([ ] = whole epoch):' } ...
            { 'style' 'edit' 'string'  defaultbase } ...
            { 'style' 'text' 'string' 'Or remove baseline points vector (ex:1:56):' } ...
            { 'style' 'edit' 'string' '' } ...
            { 'style' 'text' 'string' 'Note: press Cancel if you do not want to remove the baseline' } };
        uigeom = { [3 1] [3 1] [1] };
    else
        uilist = { { 'style' 'text' 'string' 'Removing the mean of each data channel (press cancel to skip)' } };
        uigeom = { [1] };
    end

    % add channel selection
    % --------------
    if ~isempty(EEG(1).chanlocs)
        tmpchanlocs = EEG(1).chanlocs;
    else
        tmpchanlocs = [];
        for index = 1:EEG(1).nbchan
            tmpchanlocs(index).labels = int2str(index);
            tmpchanlocs(index).type = '';
        end
    end
    cb_type = 'pop_chansel(get(gcbf, ''userdata''), ''field'', ''type'',   ''handle'', findobj(''parent'', gcbf, ''tag'', ''chantypes''));';
    cb_chan = 'pop_chansel(get(gcbf, ''userdata''), ''field'', ''labels'', ''handle'', findobj(''parent'', gcbf, ''tag'', ''channels''));';
    enableChan = fastif(length(EEG) > 1, 'off', 'on') ;
    uilist = { uilist{:} ...
        { 'style' 'text'       'string' 'Channel type(s)' 'enable' enableChan} ...
        { 'style' 'edit'       'string' '' 'tag' 'chantypes' 'enable' enableChan}  ...
        { 'style' 'pushbutton' 'string' '...'  'callback' cb_type  'enable' enableChan} ...
        { 'style' 'text'       'string' 'OR channel(s) (default all)'  'enable' enableChan} ...
        { 'style' 'edit'       'string' '' 'tag' 'channels'  'enable' enableChan}  ...
        { 'style' 'pushbutton' 'string' '...' 'callback' cb_chan  'enable' enableChan}
        };
    uigeom = { uigeom{:} [2 1.5 0.5] [2 1.5 0.5] };
    [result, usrdat, sres2, sres] = inputgui( 'uilist', uilist, 'geometry', uigeom, 'title', 'Baseline removal - pop_rmbase()', 'helpcom', 'pophelp(''pop_rmbase'');', 'userdata', tmpchanlocs);
    if isempty(result), return; end

    % decode parameters
    % -----------------
    if EEG(1).trials > 1
        if numel(result) < 2 || ((isempty(result{1}) || strcmp(result{1},'[]') ) ...
                && (isempty(result{2}) || strcmp(result{2},'[]')))
            timerange = [ EEG(1).times(1) EEG(1).times(end) ]; % whole epoch latency range
            pointrange = [];
            fprintf('pop_rmbase(): using whole epoch as baseline.\n');
        else
            timerange  = eval( [ '[' result{1} ']' ] );
            pointrange = eval( [ '[' result{2} ']' ] );
        end
    end
    if ~isempty(sres.chantypes)
        chanlist = eeg_decodechan(EEG.chanlocs, parsetxt(sres.chantype), 'type');
    elseif ~isempty(sres.channels)
        chanlist = eeg_decodechan(EEG(1).chanlocs, sres.channels);
    end

end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 2
        [ EEG, com ] = eeg_eval( 'pop_rmbase', EEG, 'warning', 'on', 'params', { timerange pointrange } );
    else
        [ EEG, com ] = eeg_eval( 'pop_rmbase', EEG, 'params', { timerange pointrange } );
    end
    return;
end

flag_timerange = 1; % provide timerange as input (so use in history)
if ~isempty(timerange)
    if timerange(1) < EEG.times(1) || timerange(end) > EEG.times(end)
        error('pop_rmbase(): Bad time range');
    end
    pointrange = find( EEG.times >= timerange(1) & EEG.times <= timerange(2));
elseif ~isempty(pointrange)
    flag_timerange = 0;
    if pointrange(1) < 1, pointrange(1) = 1; end
    if pointrange(end) > EEG.pnts, pointrange(end) = EEG.pnts; end
else
    pointrange = [1:EEG.pnts];
end

%
% Respect excised data boundaries if continuous data
% ---------------------------------------------------
fprintf('pop_rmbase(): Removing baseline...\n');
if EEG.trials == 1 && ~isempty(EEG.event) ...
                     && isfield(EEG.event, 'type') ...
                        && ischar(EEG.event(1).type)
    tmpevent = EEG.event;
    boundaries = eeg_findboundaries(tmpevent);
	if ~isempty(boundaries) % this is crashing
        fprintf('Pop_rmbase(): finding continuous data discontinuities\n');
        boundaries = round([ tmpevent(boundaries).latency ] -0.5-pointrange(1)+1);
        boundaries(boundaries>=pointrange(end)-pointrange(1)) = [];
        boundaries(boundaries<1) = [];
        boundaries = [0 boundaries pointrange(end)-pointrange(1)+1];
        for index=1:length(boundaries)-1
            tmprange = [boundaries(index)+1:boundaries(index+1)];
            if length(tmprange) > 1
                EEG.data(chanlist,tmprange) = rmbase( EEG.data(:,tmprange), length(tmprange), ...
                                                   [1:length(tmprange)]);
            elseif length(tmprange) == 1
                EEG.data(chanlist,tmprange) = 0;
            end
        end
    else
        EEG.data(chanlist,:) = rmbase( EEG.data(chanlist,:), EEG.pnts, pointrange );
    end
else
    for indc = chanlist
        tmpmean  = mean(double(EEG.data(indc,pointrange,:)),2);
        EEG.data(indc,:,:) = EEG.data(indc,:,:) - repmat(tmpmean, [1 EEG.pnts 1]);
    end
%    EEG.data = rmbase( reshape(EEG.data, EEG.nbchan, EEG.pnts*EEG.trials), EEG.pnts, pointrange );
end

EEG.data = reshape( EEG.data, EEG.nbchan, EEG.pnts, EEG.trials);
EEG.icaact = [];

if isequal(chanlist, [1:EEG.nbchan]), chanlist = []; end
if flag_timerange, pointrange = []; else, timerange = []; end
if isempty(chanlist)
    com = sprintf('EEG = pop_rmbase( EEG, %s);',vararg2str({timerange pointrange}));
else
    com = sprintf('EEG = pop_rmbase( EEG, %s);',vararg2str({timerange pointrange chanlist}));
end
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function res = chanlabels(this, varargin)
% Method for getting/setting the channel labels
% FORMAT res = chanlabels(this, ind, label)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: chanlabels.m 1373 2008-04-11 14:24:03Z spm $

if nargin == 3
    ind = varargin{1};
    label = varargin{2};

    if iscell(label) && length(label)>1
        if ~isempty(ind) && length(ind)~=length(label)
            error('Indices and values do not match');
        end

        if length(label)>1
            for i = 1:length(label)
                for j = (i+1):length(label)
                    if strcmp(label{i}, label{j})
                        error('All labels must be different');
                    end
                end
            end
        end

    end
end

res = getset(this, 'channels', 'label', varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef CEC2020_F4 < PROBLEM
% <single> <real>
% Expanded Rosenbrock's plus Griewangk's function

%------------------------------- Reference --------------------------------
% C .T. Yue, K. V. Price, P. N. Suganthan, J. J. Liang, M. Z. Ali, B. Y.
% Qu, N. H. Awad, and P. P Biswas, Problem definitions and evaluation
% criteria for the CEC 2020 special session and competition on single
% objective bound constrained numerical optimization, Zhengzhou University,
% China and Nanyang Technological University, Singapore, 2019.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2020.mat'),'Data');
            obj.O = Data{4}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 10
                obj.D   = 5;
                obj.Mat = Data{4}.M_5;
            elseif obj.D < 15
                obj.D   = 10;
                obj.Mat = Data{4}.M_10;
            elseif obj.D < 20
                obj.D   = 15;
                obj.Mat = Data{4}.M_15;
            else
                obj.D   = 20;
                obj.Mat = Data{4}.M_20;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = 0.05*Z*obj.Mat';
            Z = Y + 1;
            temp   = 100*(Z.^2-Z(:,[2:end,1])).^2 + (Z-1).^2;
            PopObj = 1900 + sum(temp.^2/4000-cos(temp)+1,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%CALCBLURRINESS  Calculate image blurriness
%
%     blurriness = cv.calcBlurriness(img)
%
% ## Input
% * __img__ input image.
%
% ## Output
% * __blurriness__ calculated blurriness.
%
% See also: cv.Sobel
%


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment1B
close all;
global robot
global parameters
global hfigures

%STOMP PARAMETERS
%conversion from cost to Prob factor
parameters.lambda = .4;
parameters.lambda_obstacles = .2;
%height of the obstacle
parameters.yo = 2.5;
%cost function starts at this distance
%must be below 0.3 for the 4 DOF robot
parameters.epsilon = 0.2;
%multiply noise by this facto
%parameters.noise_k = 5;
%parameters.noise_sigma_null_space = 0.01;
parameters.alpha=0.02;
parameters.time_step=0.01;

%number of waypoints
parameters.N = 12;
%number of particles
parameters.K = 1;
%repeat experiment number of times
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment1B_K1_N30.mat';
parameters.animate = 0;

parameters.obstacles = [];
%LINE 1
x1 = -1.5;
y1 = .5; %m
x2 = 0;
y2 = 2; %m
phi = 3*pi/4;
p0 = [x1 y1 0]';
pf = [x2 y2 0]';
T0 = build_T_4dof(p0, phi);
parameters.obstacles{1}.line = [p0 pf];
parameters.obstacles{1}.T0 = T0;

%repeat the experiment E times
random_manips=[];
Gout = [];
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = stomp_null_space(robot);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end

function T = build_T_4dof(p, phi)
T = [cos(phi) -sin(phi) 0 p(1);
     sin(phi) cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

 function T = build_T_sawyer(p, phi)
T = [1  -sin(phi) 0 p(1);
     0  cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

repo: ARTE


text:
function [xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, beamGeometry]=...
planCToPB(planC, beamNumber, bResx, bResy);

indexS=planC{end};
% assume the words Isocenter Coordinates are enclosed by quotes
% look for the end quotes to denote the begining of the actual numbers
s=planC{indexS.beamGeometry}(beamNumber).file{1};
ind=max(strfind(s, '"'));
isocenter=str2num(s(ind+1:end));
indiceJaws = 2;

% Note: In the RTOG files export by Pinnacle, the isocenter is defined at
% the first three field of
% planC{indexS.beamGeometry}(beamNumber).file
% Not the the first field only.
if (length(isocenter) == 1)
isocenter(1)=str2num(s(ind+1:end));
isocenter(2)=str2num(planC{indexS.beamGeometry}(beamNumber).file{2});
isocenter(3)=str2num(planC{indexS.beamGeometry}(beamNumber).file{3});
indiceJaws = 4;
end

gantryAngle=planC{indexS.beamGeometry}(beamNumber).gantryAngle;
isoDistance=planC{indexS.beamGeometry}(beamNumber).nominalIsocenterDistance;
beamEnergy = planC{indexS.beamGeometry}(beamNumber).beamEnergyMeV;

% if planC{indexS.beamGeometry}(beamNumber).file has only
% 3 entries, then this is a square field defined by collimator
  s=planC{indexS.beamGeometry}(beamNumber).file{indiceJaws};
  ind=max(strfind(s, '"'));
  xjaws=str2num(s(ind+1:end));
  if(length(xjaws)==1),
      xjaws = [xjaws/2 xjaws/2];
  end

  s=planC{indexS.beamGeometry}(beamNumber).file{indiceJaws+1};
  ind=max(strfind(s, '"'));
  yjaws=str2num(s(ind+1:end));
  if(length(yjaws)==1),
      yjaws = [yjaws/2 yjaws/2];
  end

  xjaws(1)=-xjaws(1);
  % by symmetry with the below(line 45),
  % I think I should swap the positive and negative y jaws) ??
  yjaws(3)=yjaws(1);
  yjaws(1)=-yjaws(2);
  yjaws(2)=yjaws(3);
  yjaws(3)=[];

if(length(planC{indexS.beamGeometry}(beamNumber).file)<7),
  % no MLC or block shapes
  % field is formed by the collimator jaws
  MLC = 0;
%   for i=1:yjaws(2)-yjaws(1),
%      x1(i)=xjaws(1);
%      xend(i)=xjaws(2);
%      vertices_y(i)=yjaws(1)-1+i;
%      vertices_ywidth(i) = 1;
%    end
   beamletInput.miny=yjaws(1);
  beamletInput.maxy=yjaws(2);
%   input = [];
%   input(:,1) = [x1'; xend'];
%   input(:,2) = [vertices_y'; vertices_y'];
  beamletInput.xjaws=xjaws;
  beamletInput.yjaws=yjaws;

else
  MLC = 1;
  s=planC{indexS.beamGeometry}(beamNumber).file{7};
  ind=max(strfind(s, '"'));
  numPairs=str2num(s(ind+1:end));
  for i=1:numPairs,
    input(i, :)=str2num(planC{indexS.beamGeometry}(beamNumber).file{7+i});
  end

  % +Y towards head, -Y towards feet (looking from BEV)
  input(:, 2)=-1*input(:,2);

 % now, calculate the intersect between the jaws and the MLC/block
  % taking into account the collimator rotation
  % jaws and collimator are both rotated EXACTLY the same way, so this is
  % just stupid!

  midPointx = xjaws(2) - (xjaws(2) - xjaws(1))/2;
  midPointy = yjaws(2) - (yjaws(2) - yjaws(1))/2;

 ind = find(input(:,1)>=midPointx & input(:,1)>xjaws(2));
  if(~isempty(ind)),
    input(ind, 1) = xjaws(2);
  end
  ind = find(input(:,1)<midPointx & input(:,1)<xjaws(1));
  if(~isempty(ind)),
    input(ind, 1) = xjaws(1);
  end

  ind = find(input(:,2)>=midPointy & input(:,2)>yjaws(2));
  if(~isempty(ind)),
    input(ind, 2) = yjaws(2);
  end
  ind = find(input(:,2)<midPointy & input(:,2)<yjaws(1));
  if(~isempty(ind)),
    input(ind, 2) = yjaws(1);
  end
end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % input is list of 2xn points that form the corners of the polygon
  % want to sort this into rectangles, with lower-left corners (x, y) and
  % width and height (w, h)

  % assume plan is MLC

  beamGeometry.isoDistance = isoDistance;
  beamGeometry.beamEnergy = beamEnergy;
  beamGeometry.isocenter=isocenter;
  beamGeometry.gantryAngle=gantryAngle;
  beamGeometry.collimatorAngle=planC{indexS.beamGeometry}(beamNumber).collimatorAngle;
  beamGeometry.couchAngle=planC{indexS.beamGeometry}(beamNumber).couchAngle;

  beamN=planC{indexS.beamGeometry}(beamNumber);
  if(MLC==0),
      % different function if the field shaping is done just by the
      % collimator jaws
      nBeamletx = ceil(diff(beamletInput.xjaws)/bResx)+1;
      nBeamlety = ceil(diff(beamletInput.yjaws)/bResy)+1;
      xvals = linspace(beamletInput.xjaws(1), beamletInput.xjaws(2), nBeamletx);
      yvals = linspace(beamletInput.yjaws(1), beamletInput.yjaws(2), nBeamlety);
      for i=1:length(yvals)-1,
          beamlets{i} = xvals;
          beamletsy{i} = yvals(i);
          beamletsDeltay{i} = yvals(i+1) - yvals(i);
      end

  else
      [x1, xend, vertices_y, vertices_ywidth]=cornersMLC(input);
      % redo setting y-vertices based on jaws
      if(vertices_y(1)<yjaws(1))
          vertices_y(1) = yjaws(1);
          vertices_ywidth(1) = vertices_y(2) -vertices_y(1);
      end
      if(vertices_y(end)+vertices_ywidth(end)>yjaws(2))
          vertices_ywidth(end) = yjaws(2) - vertices_y(end);
      end

      beamletInput.input=input;
      beamletInput.miny=min(input(:,2));
      beamletInput.maxy=max(input(:,2));
      beamletInput.x1=x1;
      beamletInput.xend=xend;
      beamletInput.vertices_y=vertices_y;
      beamletInput.vertices_ywidth=vertices_ywidth;
      beamletInput.xjaws=xjaws;
      beamletInput.yjaws=yjaws;

      numBeamlets=ceil((xend-x1)/bResx)+1;
      %numBeamlets_y=ceil(length(vertices_y)/bResy);

      numcopies = ceil(vertices_ywidth/bResy);
      % still need to deal with "fractional" number of copies

      for i=1:length(x1),
          beamlets{i}=linspace(x1(i), xend(i), numBeamlets(i));
      end

      cntr=1;

      beamletsOrig = beamlets;
      beamlets = [];

      for i=1:length(beamletsOrig),
          if(numcopies(i)==1),
              y=vertices_y(i);
              beamletsy{cntr} = y;
              ynext = vertices_y(i)+vertices_ywidth(i);
              beamletsDeltay{cntr} = ynext - y;
              beamlets{cntr} = beamletsOrig{i};
              cntr = cntr+1;
          else
              for k=1:numcopies(i)
                  y=vertices_y(i) + bResy*(k-1);
                  beamletsy{cntr} = y;
                  ynext=y+bResy;
                  beamletsDeltay{cntr} = ynext -y;
                  beamlets = insertCellEntry(beamlets, cntr, beamletsOrig{i});
                  cntr=cntr+1;
              end
          end
      end

  end

  xPosV = [];
  yPosV = [];
  beamlet_delta_x = [];
  beamlet_delta_y = [];

  for i = 1:length(beamlets)
      n = length(beamlets{i})-1;
      xPosV = [xPosV beamlets{i}(1:n)];
      beamlet_delta_x = [beamlet_delta_x diff(beamlets{i})];
      yPosV = [yPosV ones(1, n)*beamletsy{i}];
      beamlet_delta_y = [beamlet_delta_y ones(1, n)*beamletsDeltay{i}];
  end

  % need to move the position to the center of the voxel?
  xPosV = xPosV + beamlet_delta_x/2;
  yPosV = yPosV + beamlet_delta_y/2;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function [x1, xend, vertices_y, vertices_ywidth]=cornersMLC(input);

  input(:,1) = round(input(:,1)*10)/10;
  % assume 0.5 cm leaf width as the narrowest
  input(:,2) = round(input(:,2)*2)/2;

  input = checkInputPoints(input);
  miny=min(input(:,2));
  maxy=max(input(:,2));

  yvalues=sort(unique(input(:,2)));

  % Have to deal with the case that the 1st and/or last MLC positions
  % are caused by collimator jaws not MLC

  % still assuming 1 cm leaf widths, so

  %if((mod(miny, 1)~=0) | mod(maxy, 1)~=0)

  % miny1=floor(miny);
  % maxy1=ceil(maxy);
  % input(input(:,2)==miny, 2)=miny1;
  % input(input(:,2)==maxy, 2)=maxy1;
  %
  % miny=miny1;
  % maxy=maxy1;

  % numRect=maxy-miny+1;
  %
  % for i=1:numRect
  %     vertices{i}=sort(input(input(:,2)==i+miny-1, 1));
  %     vertices_y(i)=i+miny-1;
  % end

  numRect=length(yvalues);

  for i=1:numRect
      vertices{i}=sort(input(input(:,2)==yvalues(i), 1));
      vertices_y(i)=yvalues(i);
  end

  if(length(vertices{1}>2)),
      vertices{1} = unique(vertices{1});
  end
  if(length(vertices{end}>2)),
      vertices{end} = unique(vertices{end});
  end

  for i = 2:length(vertices)-1,
      if(length(vertices{i})>4)
          vertices{i} = unique(vertices{i});
      end
  end

  % for i = 2:length(vertices)-1,
  %     if(length(vertices{i}==2)),
  %         % either left or right side MLC positions not given because leaf
  %         % same as leaf above and below
  %

  x1(1)=vertices{1}(1);
  xend(1)=vertices{1}(end);

  %find the points in each line which are also in the next line
  % because the vertices are sorted, 1:2 will be left, and 3:4 will be right
  % vertices

  % june30, 2005, if it's a straight vertical line, there may be only 2
  % vertices!

  for i=2:numRect-1
      %     if(length(vertices{i} == 2)),
      %         x1(i) = vertices{i}(1);
      %         xend(i) = vertices{i}(2);
      %     else

      t=vertices{i}(1:2);
      if(t(1)==t(2)),
          x1(i)=t(1);
      else
          x1(i)=t(t~=x1(i-1));
      end
      t=vertices{i}(3:4);
      if(t(1)==t(2)),
          xend(i)=t(1);
      else
          xend(i)=t(t~=xend(i-1));
      end
      % end

  end

  vertices_ywidth = diff(vertices_y);
  vertices_y=vertices_y(1:end-1);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------

% Check that 'mcg_root' has been set
if ~exist(mcg_root,'dir')
    error('Error installing the package, try updating the value of mcg_root in the file "mcg_root.m"')
end

% Check that 'mcg_root' has the needed folder
if ~exist(fullfile(mcg_root,'lib'),'dir')
    error('Error installing the package, the folder "lib" not found, have you compiled it? See build.m')
end
if ~exist(fullfile(mcg_root,'src'),'dir')
    error('Error installing the package, the folder "src" not found')
end

% Install own lib
addpath(mcg_root);
addpath(fullfile(mcg_root,'lib'));
addpath(fullfile(mcg_root,'scripts'));
addpath(fullfile(mcg_root,'datasets'));
addpath(genpath(fullfile(mcg_root,'src')));

%% Check that the needed functions are compiled
% Included in our code
needed_files = {'mex_assess_one_sel','mex_base_perimeters','mex_fast_features',...
                'mex_fast_intersections', 'mex_fast_reduction', 'mex_box_reduction',...
                'mex_get_tree_cands', 'mex_prune_tree_to_regions',...
                'mex_max_margin', 'mex_hole_filling',...
                'mex_intersect_hierarchies','mex_cands2masks','mex_cands2labels','mex_ucm2hier',...
                'paretofront',...                  % Included from paretofront
                'mexRF_train', 'mexRF_predict',... % Included from RF_Reg_C
                'convConst','gradientMex','imPadMex','imResampleMex','rgbConvertMex',... % Included from piotr_toolbox
                'edgesDetectMex',...               % Included from structured_forest
                'buildW', 'mex_contour_sides', 'ucm_mean_pb',... % Included from BSR
                };
for ii=1:length(needed_files)
    if exist(needed_files{ii})~=3 %#ok<EXIST>
        error(['The needed function (' needed_files{ii} ') not found. Have you built the package properly?'])
    end
end

%% Check that the databases are available and show a warning if not
dbs        = {'pascal2012',    'SBD',   'COCO'};
im_folders = {'JPEGImages', 'images', 'images'};
all_ok = 1;
for ii=1:length(dbs)
    db = dbs{ii};
    if ~exist(fullfile(database_root_dir(db), im_folders{ii}),'dir')
        all_ok = 0;
        disp(['WARNING: Database ' db ' (folder ' im_folders{ii} ') not found in ' database_root_dir(db)])
    end
end

if ~all_ok
    disp('-- You can disable this warning in install.m --')
end

%% Clear
clear ii needed_files;
disp('-- Successful installation of MCG. Enjoy! --');

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function y = atanminuspihalf(x)
%ATANMINUSPIHALF Inverse tangent minus pi/2.
%
%   ATANMINUSPIHALF(X) is ATAN(X)-PI/2 calculated in a way that is numerically
%   better when X is large and positive.
%
%   See also TANPLUSPIHALF.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 14:55:25 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   y = zeros(size(x));

   k = x <= 0;
   y(k) = atan(x(k)) - (pi / 2);

   k = x > 0;
   y(k) = atan(-1 ./ x(k));

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function value = p27_exact ( dim_num )

%*****************************************************************************80
%
%% P27_EXACT returns the exact integral for problem 27.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  r = 0.0;
  r = p27_r8 ( 'G', 'R', r );

  c = [];
  c = p27_r8vec ( 'G', 'C', dim_num, c );

  value = 2.0^dim_num * ...
    cos ( 2.0 * pi * r + 0.5 * sum ( c(1:dim_num) ) ) * ...
    prod ( sin ( 0.5 * c(1:dim_num) ) ./ c(1:dim_num ) );

  return
end


text:
% ********************************************************************
% * mex File compiling code for Random Forest (for windows)
% * mex interface to Andy Liaw et al.'s C code (used in R package randomForest)
% * Added by Abhishek Jaiantilal ( abhishek.jaiantilal@colorado.edu )
% * License: GPLv2
% * Version: 0.02
% ********************************************************************/

function compile_windows

%need to do tricks for making Makfile run on windows as one needs cygwin.
%so instead use mex to compile everything up.

%system('del *.mexw32;');

mex src/cokus.cpp src/reg_RF.cpp src/mex_regressionRF_train.cpp   -DMATLAB -output mexRF_train
mex src/cokus.cpp src/reg_RF.cpp src/mex_regressionRF_predict.cpp   -DMATLAB  -output mexRF_predict

fprintf('\n Mex`s compiled correctly\n')

fprintf('As this package has precompiled mex files ignore this file i.e. compile_windows.m\n')

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function pass = test_grad( )
% Test curl

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Check gradient of an empty diskfun is an empty diskfunv.
u = diskfun;
f = grad(u);
pass(1) = isempty(f) & isa(f,'diskfunv');

% Check gradient of the zero function is zero
f = diskfun(@(x,y) 0*x);
pass(2) = norm(grad(f)) < tol;

%
% Check gradient of a non-zero function gives the correct result
%

% Example 1
u = grad(diskfun(@(x,y,z) x.^4.*(y-y.^2)) );
% Exact gradient
f = diskfun(@(x,y) 4*x.^3.*(y-y.^2));
g = diskfun(@(x,y) x.^4.*(1-2*y));
exact = diskfunv(f,g);
pass(3) = norm(u-exact) < tol;

% Example 2
u = grad(diskfun(@(x,y) cos(4*x).*sin(x.*y)));
% Exact gradient
f = diskfun(@(x,y) -4*sin(4*x).*sin(x.*y)+cos(4*x).*cos(x.*y).*y);
g = diskfun(@(x,y) cos(4*x).*cos(x.*y).*x);
exact = diskfunv(f,g);
pass(4) = norm(u-exact) < 2*tol;

% Example 3
u = grad(diskfun(@(x,y) exp(-3*(x.^2+(y+.2).^2)) ));
% Exact gradient
f = diskfun(@(x,y) -6.*x.*exp(-3*(x.^2+(y+.2).^2)));
g = diskfun(@(x,y) -6.*(y+.2).*exp(-3*(x.^2+(y+.2).^2)));
exact = diskfunv(f,g);
pass(5) = norm(u-exact) < 10*tol;

% Check that the grad and gradient give the same result.
gradientu = gradient(diskfun(@(x,y) exp(-3*(x.^2+(y+.2).^2)) ));
pass(6) = isequal(u,gradientu);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% StackExchange Mathematics Q3872982
% https://math.stackexchange.com/questions/3872982
% Solve Linear Least Squares with Squared L2 Norm Regularization (Tikhonov
% / Ridge Regression) with Non Negativity Constraint Using FISTA
% References:
%   1.
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     26/12/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 0;2165;42; %<! Set to 0 for Random

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numRows         = 100;
numCols         = 20;
paramLambda     = 0.1;

% Solvers parameters
numIterations       = 5000;
stepSizeGd          = 5e-5;
stepSizeMomentum    = 0.8;
stepSizeAccel       = 0.4;

%% Load / Generate Data

mA = randn(numRows, numCols);
vB = randn(numRows, 1);

hObjFun = @(vX) 0.5 * sum((mA * vX - vB) .^ 2) + (0.5 * paramLambda * sum(vX .^ 2));

mAA = mA.' * mA;
vAb = mA.' * vB;

hG = @(vX) (mAA * vX) - vAb + (paramLambda * vX); %<! Gradient of the Objection Function
hP = @(vX) max(vX, 0); %<! Projection onto the constraint set

solverIdx       = 0;
cMethodString   = {};

mObjFunValMse   = zeros([numIterations, 1]);
mSolMse         = zeros([numIterations, 1]);

%% Solution by CVX

solverString = 'CVX';

% cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
% cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
% cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    minimize( 0.5 * sum_square(mA * vX - vB) + (0.5 * paramLambda * square_pos(norm(vX))) );
    subject to
        vX >= 0;
cvx_end

runTime = toc(hRunTime);

% vX = mX(:);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The ', solverString, ' Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

sCvxSol.vXCvx     = vX;
sCvxSol.cvxOptVal = hObjFun(vX);

%% Solution by Projected Gradient Descent

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Method'];

hRunTime = tic();

[vX, mX] = ProjectedGd(zeros(numCols, 1), hG, hP, numIterations, stepSizeGd);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with Momentum

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with Momentum'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdMomentum(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd, stepSizeMomentum);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with Nesterov Acceleration

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with Nesterov Acceleration'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdAccel(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd, stepSizeAccel);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Solution by Projected Gradient Descent with FISTA

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Descent with FISTA Acceleration'];

hRunTime = tic();

[vX, ~, mX] = ProjectedGdFista(zeros(numCols, 1), zeros(numCols, 1), hG, hP, numIterations, stepSizeGd);
% [vX, mX] = SolveLsFista(zeros(numCols, 1), mA, vB, paramLambda, numIterations, stepSizeGd);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

[mObjFunValMse, mSolMse] = UpdateAnalysisData(mObjFunValMse, mSolMse, mX, hObjFun, sCvxSol, solverIdx);

%% Display Results

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);

hAxes       = subplot(2, 1, 1);
hLineSeries = plot(1:numIterations, 10 * log10(mObjFunValMse));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['Objective Function Value vs. Optimal Value (CVX)'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ 10 \log_{10} {\left( \left| f \left( x \right) - f \left( {x}_{CVX} \right) \right| \right)}^{2} $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend(cLegendString);

hAxes       = subplot(2, 1, 2);
hLineSeries = plot(1:numIterations, 10 * log10(mSolMse));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['Solution Error Norm'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'Iteration Number', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', '$ 10 \log_{10} \left( {\left\| x - {x}_{CVX} \right\|}_{2}^{2} \right) $', ...
    'FontSize', fontSizeAxis, 'Interpreter', 'latex');
set(hAxes, 'XLim', [1, numIterations]);
hLegend = ClickableLegend(cLegendString);

if(generateFigures == ON)
    % saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

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

author: RoyiAvital

repo: StackExchangeCodes


text:
function cfactor = SC_AmplitudeFactor(blk, page, edgewt, sccalib)
% Amplitude Factor: ce_i with Edge Response Function (edgewt)
% Reference: Improved scatter correction using adaptive scatter kernel superposition
% Author: Yi Du (yi.du@hotmail.com)
% Date: 2021-05-24

%% group number
ngroup = length(sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel);

% Amplitude factor groups
cfactor = [];

for ii=1:ngroup
    tmp = sccalib.CalibrationResults.ObjectScatterModels.ObjectScatterModel{ii}.ObjectScatterFit;
    % Amplitude Factor
    % unit: mm - > cm
    A = str2double(tmp.A.Text) / 10;
    % unitless
    alpha = str2double(tmp.alpha.Text);
    beta = str2double(tmp.beta.Text);

    % fill holes
    term = (page + eps)./(blk + eps);
    logterm = -log(term);
    logterm(logterm<0) = NaN;
    logterm = single(inpaint_nans(double(logterm), 2));

    % amplitude factor wi edge response function as well
    cfactor(:,:,ii) = A .*edgewt .* (term).^(alpha) .* ( logterm ).^(beta);
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
% AssignmentToIndex Convert assignment to index.
%
%   I = AssignmentToIndex(A, D) converts an assignment, A, over variables
%   with cardinality D to an index into the .val vector for a factor.
%   If A is a matrix then the function converts each row of A to an index.
%
%   See also IndexToAssignment.m and SampleFactors.m

function I = AssignmentToIndex(A, D)

D = D(:)'; % ensure that D is a row vector
if (any(size(A) == 1)),
    I = cumprod([1, D(1:end - 1)]) * (A(:) - 1) + 1;
else
    I = sum(repmat(cumprod([1, D(1:end - 1)]), size(A, 1), 1) .* (A - 1), 2) + 1;
end;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function C = coeffs3(f, m, n, p)
% COEFFS3   Trivariate Cheybshev-Fourier-Fourier expansion coefficients of F.
%   C = COEFFS3(F) returns the tensor of trivariate coefficients.  The
%   coefficients are arranged so they correspond the spherical coordinates
%   radial-azimuthal-polar.
%
%   X = COEFFS3(F, M, N, P) is the same as above but it returns
%   coefficients as an M x N x P tensor.

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

if ( nargin == 1 )
    % Return the coefficients
    C = f.coeffs;
else
    if ( nargin == 2 )
        n = m;
        p = m;
    elseif (nargin == 3 )
        p = n;
    end

    %% Copy the diff alias when we solve it below

    [mf, nf, pf] = size(f);
    F = f.coeffs;

    G = zeros(m,n,pf);

    for k = 1:pf
        G(:,:,k) = chebtech2.alias(trigtech.alias(F(:,:,k).',n).',m);
    end

    G = permute(G,[3,2,1]);
    C = zeros(p,n,m);

    for k = 1:m
       C(:,:,k) = trigtech.alias(G(:,:,k),p);
    end

    C = permute(C,[3,2,1]);
end


text:
function f = rhs ( node_num, node_xy, time )

%*****************************************************************************80
%
%% RHS gives the right-hand side of the differential equation.
%
%  Discussion:
%
%    We assume that the equation to be solved is
%
%      dUdT - Laplacian U + K * U = F
%
%    with
%
%      K = 0,
%
%    and
%
%      F = (2*pi*pi-1)*sin(pi*x)*sin(pi*y)*exp(-t).
%
%    The exact solution is:
%
%      U = sin(pi*x) * sin(pi*y) * exp(-t)
%
%  Modified:
%
%    07 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of the points.
%
%    Input, real TIME, the current time (assumed to be
%    the initial time).
%
%    Output, real F(NODE_NUM), the value of the right hand side.
%
  f(1:node_num) = ( 2.0 * pi * pi - 1.0 ) ...
     * sin ( pi * node_xy(1,1:node_num) ) ...
    .* sin ( pi * node_xy(2,1:node_num) ) * exp ( - time );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef (Abstract) math_model < mp.element_container & opt_model
%MP.MATH_MODEL  MATPOWER mathematical model abstract base class.
%   ?
%
%   MP.MATH_MODEL provides properties and methods related to the specific
%   problem specification being solved (e.g. power flow, continuation
%   power flow, optimal power flow, etc.) ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-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
        aux_data    %% struct of auxiliary data relevant to the model,
                    %% e.g. can be passed to model constraint functions
    end

    methods
        function tag = task_tag(obj)
            error('mp.math_model/task_tag: must be implemented in subclass');
        end

        function name = task_name(obj)
            error('mp.math_model/task_name: must be implemented in subclass');
        end

        function tag = form_tag(obj)
            error('mp.math_model/form_tag: must be implemented in subclass');
        end

        function name = form_name(obj)
            error('mp.math_model/form_name: must be implemented in subclass');
        end

        function obj = build(obj, nm, dm, mpopt)
            %% 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:  Initialize MP_IDX_MANAGER fields here, if needed,
            %%              after object construction, but before object use.
            if isempty(obj.var)         %% only if not already initialized
                obj.init_set_types();
            end

            %% create element objects for each class with data
            obj.elements = mp.mapped_array();
            for c = obj.element_classes
                mme = c{1}();       %% element constructor
                if dm.online(mme.name)      %% dm element exists
                    obj.elements.add_elements(mme, mme.name);
                end
            end

            obj.add_aux_data(nm, dm, mpopt);
            obj.add_vars(nm, dm, mpopt);
            obj.add_constraints(nm, dm, mpopt);
            obj.add_costs(nm, dm, mpopt);
        end

        function display(obj)
            fprintf('MATH MODEL CLASS : %s\n', class(obj));
            fprintf('    TASK NAME               : %s\n', obj.task_name());
            fprintf('    TASK TAG                : %s\n', obj.task_tag());
            fprintf('    FORMULATION NAME        : %s\n', obj.form_name());
            fprintf('    FORMULATION TAG         : %s\n', obj.form_tag());
            display@opt_model(obj)

            %% elements
            fprintf('\nELEMENTS\n')
            fprintf('========\n')
            fprintf('  name              class\n');
            fprintf(' ----------------  --------------------\n');
            for k = 1:length(obj.elements)
                mme = obj.elements{k};
                fprintf('  %-13s     %s\n', mme.name, class(mme));
            end
        end

        function obj = add_aux_data(obj, nm, dm, mpopt)
            %% create aux_data struct
            obj.aux_data = obj.build_aux_data(nm, dm, mpopt);
        end

        function ad = build_base_aux_data(obj, nm, dm, mpopt)
            %% get model variables
            vvars = nm.model_vvars();
            zvars = nm.model_zvars();
            vars = {vvars{:} zvars{:}};
            vals = cellfun(@(x)nm.params_var(x), vars, 'UniformOutput', false);

            %% get node types
            [ref, pv, pq, by_elm] = nm.node_types(nm, dm);

            %% create aux_data struct
            ad = cell2struct(...
                {vals{:}, {}, ref, pv, pq, ...
                    length(ref), length(pv), length(pq), by_elm}, ...
                {vars{:}, 'var_map', 'ref', 'pv', 'pq', ...
                    'nref', 'npv', 'npq', 'node_type_by_elm'}, 2);
        end

        function obj = add_vars(obj, nm, dm, mpopt)
            obj.add_system_vars(nm, dm, mpopt);

            %% each element adds its OPF variables
            for k = 1:length(obj.elements)
                obj.elements{k}.add_vars(obj, nm, dm, mpopt);
            end
        end

        function obj = add_system_vars(obj, nm, dm, mpopt)
        end

        function obj = add_constraints(obj, nm, dm, mpopt)
            obj.add_system_constraints(nm, dm, mpopt);

            %% each element adds its OPF variables
            for k = 1:length(obj.elements)
                obj.elements{k}.add_constraints(obj, nm, dm, mpopt);
            end
        end

        function obj = add_system_constraints(obj, nm, dm, mpopt)
            %% node balance constraints
            obj.add_node_balance_constraints(nm, dm, mpopt);
        end

        function obj = add_node_balance_constraints(obj, nm, dm, mpopt)
        end

        function obj = add_costs(obj, nm, dm, mpopt)
            obj.add_system_costs(nm, dm, mpopt);

            %% each element adds its OPF variables
            for k = 1:length(obj.elements)
                obj.elements{k}.add_costs(obj, nm, dm, mpopt);
            end
        end

        function obj = add_system_costs(obj, nm, dm, mpopt)
        end

        function opt = solve_opts(obj, nm, dm, mpopt)
            opt = struct();
        end

        function nm_vars = update_nm_vars(obj, mmx, nm)
            %% Returns a struct with the network model variables as fields.
            %% The ad.var_map cell array is used to track mappings
            %% of math model variables back to network model variables.
            %% Each entry is a cell array:
            %%  {nm_var_type, nm_i1, nm_iN, nm_idx, mm_i1, mm_iN, mm_idx}
            %% where
            %%  nm_var_type - network model variable type (e.g. va, vm, zr, zi)
            %%  nm_i1 - starting index for network model variable type
            %%  nm_iN - ending index for network model variable type
            %%  nm_idx - vector of indices for network model variable type
            %%  mm_i1 - starting index for math model variable
            %%  mm_iN - ending index for math model variable
            %%  mm_idx - vector of indices for math model variable
            %% Uses either i1:iN (if i1 is not empty) or idx as the indices,
            %% unless both are empty, in which case it uses ':'

            ad = obj.aux_data;
            vvars = nm.model_vvars();
            zvars = nm.model_zvars();
            vars = {vvars{:} zvars{:}};
            nm_vars = cell2struct( ...
                cellfun(@(x)ad.(x), vars, 'UniformOutput', false), ...
                vars, 2);

            %% copy mm state values to nm state based on ad.var_map
            for k = 1:length(ad.var_map)
                d = ad.var_map{k};

                %% right hand side from mmx
                if ~isempty(d{5})       %% use i1:iN
                    rhs = mmx(d{5}:d{6});
                elseif isempty(d{7})    %% use :
                    rhs = mmx;
                else                    %% use idx
                    rhs = mmx(d{7});
                end

                %% left hand side to nm_vars
                name = d{1};
                if ~isempty(rhs)
                    if ~isempty(d{2})       %% use i1:iN
                        nm_vars.(name)(d{2}:d{3}) = rhs;
                    elseif isempty(d{4})    %% use :
                        nm_vars.(name) = rhs;
                    else                    %% use idx
                        nm_vars.(name)(d{4}) = rhs;
                    end
                end
            end
        end

        function dm = data_model_update(obj, nm, dm, mpopt)
            %% each element updates its data model
            for k = 1:length(obj.elements)
                obj.elements{k}.data_model_update(obj, nm, dm, mpopt);
            end
        end

        function nm = network_model_x_soln(obj, nm)
            %% convert solved state from math model to network model soln
            [nm.soln.v, nm.soln.z, nm.soln.x] = ...
                obj.convert_x_m2n(obj.soln.x, nm);
        end
    end     %% methods
end         %% classdef


text:
function [packet, state]= bbci_control_BGUIERP(cfy_out, state, event, opt)
%same as BBCI_CONTROL_ERPSPELLER - just need another packet.
% 02-2012 Javier Pascual

    if(strcmp(opt.kind, 'MI')),
        packet= {['i:' opt.cfy], cfy_out};
    else,

        if isempty(state),
          state.counter= zeros(1, opt.nClasses);
          state.output= zeros(1, opt.nClasses);
        end

        this_cue= 1 + mod(event.desc-11, 10);
        state.counter(this_cue)= state.counter(this_cue) + 1;
        state.output(this_cue)= state.output(this_cue) + cfy_out;

        %fprintf('cfy_out = %.2f\n', cfy_out);
        %s = sprintf('output %d(%d) = %.2f %.2f %.2f %.2f %.2f %.2f\n',sum(state.counter),opt.nClasses*opt.nSequences,state.output(1),state.output(2),state.output(3),state.output(4),state.output(5),state.output(6) );
        %fprintf(s);
        %msgbox(s);

        if sum(state.counter) >= opt.nClasses*opt.nSequences,
          idx= find(state.counter>0);  % avoid divide by zero
          state.output(idx)= state.output(idx) ./ state.counter(idx);
          [max_score, selected_class]= min(state.output);

          %fprintf('  -> selected class = %d\n',selected_class );

          packet= {'i:erp_output', selected_class};
          state.counter(:)= 0;
          state.output(:)= 0;
        else
          packet= [];
        end;
    end;
end

author: bbci