text:
%% CONTRAST2_BATCH_ITHACA uses the BATCH command to run the CONTRAST2 code on the Ithaca cluster.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONTRAST2_BATCH_ITHACA\n' );
  fprintf ( 1, '  Run CONTRAST2_SCRIPT on Ithaca.\n' );
%
%  Read a black and white image and filter it locally.
%
  x = imread ( 'surfsup.tif' );
  xf = nlfilter ( x, [3,3], @contrast_enhance );
  xf = uint8 ( xf );
%
%  Set the number of workers to use.
%
  n = 4;
%
%  The BATCH command sends the script for execution.
%
  job = batch ( 'contrast2_script', ...
    'Configuration', 'ithaca_2009b', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'contrast2_fun', 'contrast_enhance', 'surfsup.tif' }, ...
    'CurrentDirectory', '/home/burkardt/matlab', ...
    'matlabpool', n );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( job );
%
%  LOAD makes the script's workspace available.
%
%  y = output #1.
%
  load ( job );
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( job );
%
%  Display the original, single filtered, and SPMD filtered versions.
%
  figure ( 1 );

  subplot ( 1, 3, 1 );
  imshow ( x );
  title ( 'Original image' );

  subplot ( 1, 3, 2 );
  imshow( xf );
  title( 'Filtered on Client' );

  subplot ( 1, 3, 3 );
  imshow ( y );
  title ( 'Filtered on 4 SPMD Workers' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CONTRAST2_BATCH_ITHACA\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function hRD = plugin_rd_rx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.pluto.plugin_rd('Rx');
AnalogDevices.pluto.add_io(hRD,'Rx');

author: analogdevicesinc


text:
function value = r4_sign ( x )

%*****************************************************************************80
%
%% R4_SIGN returns the sign of an R4.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( x < 0.0 )
    value = -1.0;
  else
    value = +1.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function y = reale(x, arg2, arg3)
%| Return real part of complex data (with error checking).
%function y = reale(x, arg2, arg3)
%|
%| y = reale(x)
%| y = reale(x, tol) [default tol is 1e-13 for double, else 1e-6]
%| y = reale(x, 'warn', 'message')
%| y = reale(x, 'error')
%| y = reale(x, 'report')
%| y = reale(x, 'prompt')
%| y = reale(x, 'disp')
%|
%| Checks that imaginary part is negligible (warns etc. if not).
%|
%| Copyright Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(x, 'test'), reale_test, return, end

com = 'error';
if isa(x, 'double')
	tol = 1e-13;
else
	tol = 1e-6;
end

if nargin > 1
	if ischar(arg2)
		com = arg2;
	elseif isnumeric(arg2)
		tol = arg2;
	end
end

switch com
case {'disp','prompt', 'report'}
	;
case 'warn'
	onlywarn = 1;
case 'error'
	onlywarn = 0;
otherwise
	fail('bad argument "%s"', com)
end

max_abs_x = max(abs(x(:)));
if max_abs_x == 0
	if any(imag(x(:)) ~= 0)
		fail 'max real 0, but imaginary!'
	else
		y = real(x);
		return
	end
end

frac = max(abs(imag(x(:)))) / max_abs_x;
if streq(com, 'report')
	printm('imaginary part %g%%', frac * 100)
	return
end

if frac > tol
	[cname line] = caller_name;
	t = sprintf('%s(%d): %s: imaginary fraction of %s [class %s] is %g', ...
		cname, line, mfilename, inputname(1), class(x), frac);
	if isvar('arg3')
		t = [t ', ' arg3];
	end
	if streq(com, 'disp')
		disp(t)

	elseif streq(com, 'prompt')
		printm('reale() called for input with imaginary part %g%%', frac * 100)
		printm('reale() called in context where a large imaginary part')
		printm('is likely an *error*.  proceed with caution!')
		t = input('proceed? [y|n]: ', 's');
		if isempty(t) || t(1) ~= 'y'
			printm('ok, aborting is probably wise!')
			error ' '
		end

	elseif onlywarn
		disp(t)
	else
		fail(t)
	end
end
y = real(x);

function reale_test
x = 7 + 1i*eps;
reale(x, 'warn');
reale(x, 'prompt');
%reale(x, 'report'); % check error reporting
%reale(x, eps/100) % check error faulting

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function model = person_init_grammar()
% Initialize the structure and filters of the person grammar model.
%   model = person_init_grammar()

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

model = train_full_person_2x_res();
model = make_person_grammar_occ_def(model);

% --------------------------------------------------------------------
function model = train_full_person_2x_res()
% Train a root filter for (roughly) fully visible people.
%   model = train_full_person_2x_res()
%
%   This root filter is at about 2x the resolution of the root filter
%   used for the mixture models because it will be divided up into
%   high resolution parts. Fully visible people are estimated to be
%   the top third most vertically elongated examples.

conf = voc_config();
cachedir = conf.paths.model_dir;

try
  load([cachedir 'person_full_person_2x']);
catch
  seed_rand();
  cls = 'person';
  note = 'full person trained at 2x resolution';
  max_num_examples = conf.training.cache_example_limit;
  fg_overlap = conf.training.fg_overlap;

  [pos, neg, impos] = pascal_data(cls, conf.pascal.year);
  % Split data by aspect ratio into n groups
  spos = split(pos, 3);
  spos = spos{3};

  try
    load([cachedir 'person_model_full_person_2x_1_1_1']);
    load([cachedir 'person_model_full_person_2x_pos_inds']);
  catch
    model = root_model(cls, spos, note, 8, [23 8]);
    % Allow root detections in the first pyramid octave
    bl = model.rules{model.start}(1).loc.blocklabel;
    model.blocks(bl).w(:) = 0;

    inds = lrsplit(model, spos);
    % Train with warpped positives and random negatives
    model = train(model, spos(inds), neg, true, true, 1, 1, ...
                  max_num_examples, fg_overlap, 0, false, 'full_person_2x_1');
    save([cachedir 'person_model_full_person_2x_pos_inds'], 'inds');
  end
  % Train with latent postives and hard negatives
  model = train(model, spos(inds), neg(1:200), false, false, 1, 20, ...
                max_num_examples, fg_overlap, 0, false, 'full_person_2x_2');

  save([cachedir cls '_full_person_2x'], 'model');
end

%-------------------------------------------------------------------------
function M = make_person_grammar_occ_def(model)
% Construct a person grammar model with 6 levels of occlusion, deformable
% parts, and an 'occluder' part.
%   M = make_person_grammar_occ_def(model)
%
% Return value
%   M         The person grammar model
%
% Arugment
%   model     The root only full person model output by
%             train_full_person_2x_res()

conf = voc_config();
cachedir = conf.paths.model_dir;

cls = 'person';
note = 'NIPS 2011 person grammar model';

%{{{ Initialize a fresh model
  M = model_create(cls, note);
  M.type = model_types.Grammar;

  % Start non-terminal
  [M, Q] = model_add_nonterminal(M);
  M.start = Q;
%}}}

%{{{ Filters
  % Get the root filter for fully visible people
  w = model_get_block(model, model.filters(1));
  % Divide it into 6 pieces
  X_f = w(1:8, :, :);     % Head+shoulders
  Y1_f = w(9:11, :, :);   % and so on ...
  Y2_f = w(12:14, :, :);
  Y3_f = w(15:17, :, :);
  Y4_f = w(18:20, :, :);
  Y5_f = w(21:23, :, :);  % ... feet
  % Occluder filter (initially all zeros)
  O_f = zeros(4, 8, size(w, 3));

  % Create the set of terminals: one for each filter
  % and its left-right mirrored counterpart
  % Notation: S_{l,r}f  =>  terminal symbol S, l = left, r = right
  [M, X_lf]  = model_add_terminal(M, 'w', X_f);
  [M, X_rf]  = model_add_terminal(M, 'mirror_terminal', X_lf);
  [M, Y1_lf] = model_add_terminal(M, 'w', Y1_f);
  [M, Y1_rf] = model_add_terminal(M, 'mirror_terminal', Y1_lf);
  [M, Y2_lf] = model_add_terminal(M, 'w', Y2_f);
  [M, Y2_rf] = model_add_terminal(M, 'mirror_terminal', Y2_lf);
  [M, Y3_lf] = model_add_terminal(M, 'w', Y3_f);
  [M, Y3_rf] = model_add_terminal(M, 'mirror_terminal', Y3_lf);
  [M, Y4_lf] = model_add_terminal(M, 'w', Y4_f);
  [M, Y4_rf] = model_add_terminal(M, 'mirror_terminal', Y4_lf);
  [M, Y5_lf] = model_add_terminal(M, 'w', Y5_f);
  [M, Y5_rf] = model_add_terminal(M, 'mirror_terminal', Y5_lf);
  [M, O_lf]  = model_add_terminal(M, 'w', O_f);
  [M, O_rf]  = model_add_terminal(M, 'mirror_terminal', O_lf);
%}}}

%{{{ Nonterminals and deformation rules for placing the terminals
  % Initial deformation model
  def_w = [0.01 0 0.01 0];

  % Symbols for the lhs of the deformation rules
  % Parts 2-6 and the occluder can move, part 1 (head) is fixed
  [M, Y1_l] = model_add_nonterminal(M);
  [M, Y2_l] = model_add_nonterminal(M);
  [M, Y3_l] = model_add_nonterminal(M);
  [M, Y4_l] = model_add_nonterminal(M);
  [M, Y5_l] = model_add_nonterminal(M);
  [M, O_l]  = model_add_nonterminal(M);
  [M, Y1_r] = model_add_nonterminal(M);
  [M, Y2_r] = model_add_nonterminal(M);
  [M, Y3_r] = model_add_nonterminal(M);
  [M, Y4_r] = model_add_nonterminal(M);
  [M, Y5_r] = model_add_nonterminal(M);
  [M, O_r]  = model_add_nonterminal(M);

  [M, rule] = model_add_def_rule(M, Y1_l, Y1_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, Y1_r, Y1_rf, 'mirror_rule', rule);

  [M, rule] = model_add_def_rule(M, Y2_l, Y2_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, Y2_r, Y2_rf, 'mirror_rule', rule);

  [M, rule] = model_add_def_rule(M, Y3_l, Y3_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, Y3_r, Y3_rf, 'mirror_rule', rule);

  [M, rule] = model_add_def_rule(M, Y4_l, Y4_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, Y4_r, Y4_rf, 'mirror_rule', rule);

  [M, rule] = model_add_def_rule(M, Y5_l, Y5_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, Y5_r, Y5_rf, 'mirror_rule', rule);

  [M, rule] = model_add_def_rule(M, O_l, O_lf, 'def_w', def_w);
  [M, rule] = model_add_def_rule(M, O_r, O_rf, 'mirror_rule', rule);
%}}}

%{{{ Add structural rules for deriving parts

  % The following rules derive either the (deformed) left-facing or
  % right-facing subtype.
  %  X  -> X_lf | X_rf
  %  Yi -> Yi_l | Yi_r
  %  O  -> O_l  | O_r

  [M, X]  = model_add_nonterminal(M);
  [M, Y1] = model_add_nonterminal(M);
  [M, Y2] = model_add_nonterminal(M);
  [M, Y3] = model_add_nonterminal(M);
  [M, Y4] = model_add_nonterminal(M);
  [M, Y5] = model_add_nonterminal(M);
  [M, O]  = model_add_nonterminal(M);

  [M, rule] = model_add_struct_rule(M, X, X_lf, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, X, X_rf, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, Y1, Y1_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, Y1, Y1_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, Y2, Y2_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, Y2, Y2_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, Y3, Y3_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, Y3, Y3_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, Y4, Y4_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, Y4, Y4_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, Y5, Y5_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, Y5, Y5_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;

  [M, rule] = model_add_struct_rule(M, O, O_l, {[0 0 0]});
  [M, rule] = model_add_struct_rule(M, O, O_r, {[0 0 0]}, 'mirror_rule', rule);
  M.blocks(rule.offset.blocklabel).learn = 0;
%}}}

%{{{ Add top-level structural rules

  % The following rules derive people of varying degrees of visibility
  % from head to toe.
  %  Q -> X O | X Y1 O | X Y1 Y2 O | ... | X Y1 Y2 Y3 Y4 O | X Y1 Y2 Y3 Y4 Y5

  [M, rule] = model_add_struct_rule(M, Q, [X O], ...
                                    {[0 0 0], [0 8 0]}, ...
                                    'detection_window', [8 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;

  [M, rule] = model_add_struct_rule(M, Q, [X Y1 O], ...
                                    {[0 0 0], [0 8 0], [0 11 0]}, ...
                                    'detection_window', [11 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;

  [M, rule] = model_add_struct_rule(M, Q, [X Y1 Y2 O], ...
                                    {[0 0 0], [0 8 0], [0 11 0], [0 14 0]}, ...
                                    'detection_window', [14 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;

  [M, rule] = model_add_struct_rule(M, Q, [X Y1 Y2 Y3 O], ...
                                    {[0 0 0], [0 8 0], [0 11 0], [0 14 0], [0 17 0]}, ...
                                    'detection_window', [17 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;

  [M, rule] = model_add_struct_rule(M, Q, [X Y1 Y2 Y3 Y4 O], ...
                                    {[0 0 0], [0 8 0], [0 11 0], [0 14 0], [0 17 0], [0 20 0]}, ...
                                    'detection_window', [20 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;

  [M, rule] = model_add_struct_rule(M, Q, [X Y1 Y2 Y3 Y4 Y5], ...
                                    {[0 0 0], [0 8 0], [0 11 0], [0 14 0], [0 17 0], [0 20 0]}, ...
                                    'detection_window', [23 8]);
  M.blocks(rule.loc.blocklabel).learn = 1;
  M.blocks(rule.loc.blocklabel).reg_mult = 1;
%}}}

model = M;
save([cachedir cls '_simple_grammar_occ_def'], 'model');

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function lorenza
% Lorenz convection model

sigma = 16; rho = 45.92; beta = 4;  % parameters
N = 1000;         % no. of time steps
span = 0.05;      % inner iteration span
AbsTol = 1.e-5;   % absolute tolerance for ODE solver
RelTol = 1.e-5;   % relative tolerance for ODE solver

H = figure; set(H,'DefaultLineLineWidth',1.0);
options = odeset('RelTol',RelTol,'AbsTol',ones(1,3)*AbsTol);
u0 = [1;1;1];
for i = 1:N
    [t,u] = ode45(@lornz,[0 span],u0,odeset,beta,rho,sigma);
    hold on;
    plot(u(:,1),u(:,2),'r');
    u0 = u(end,:);
end

title('Attractor of Lorenz System');
xlabel('Component 1'); ylabel('Component 2');
axis off; hold off;

function dydt = lornz(t,y,beta,rho,sigma)
dydt = [sigma*(y(2)-y(1)); rho*y(1)-y(2)-y(1)*y(3); y(1)*y(2)-beta*y(3)];

author: Sable

repo: mcbench-benchmarks


text:
function K = cotmatrix3(V,T)
  % COTMATRIX3 computes cotangent matrix for 3D tetmeshes, area/mass terms
  % already cancelled out: laplacian mesh operator
  %
  % This is distinctly NOT following definition that
  % appears in the appendix of: ``Interactive Topology-aware Surface
  % Reconstruction,'' by Sharf, A. et al
  % http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
  %
  % Instead it is a purely geometric construction. Find more details in Section
  % 1.1 of "Algorithms and Interfaces for Real-Time Deformation of 2D and 3D
  % shapes" [Jacobson 2013]
  %
  % ND derivation given in "A MONOTONE FINITE ELEMENT SCHEME FOR
  % CONVECTION-DIFFUSION EQUATIONS" [Xu & ZIKATANOV 1999]
  %
  % 3D derivation given in "Aspects of unstructured grids and finite-volume
  % solvers for the Euler and Navier-Stokes equations" [Barth 1992]
  %
  % K = cotmatrix(V,T)
  % Inputs:
  %   V  #V x 3 matrix of vertex coordinates
  %   T  #T x 4  matrix of indices of tetrahedral corners
  % Output:
  %   K  #V x #V matrix of cot weights
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also cotmatrix
  %
  warning('Deprecated. Call cotmatrix directly.');
  K = cotmatrix(V,T);
end

author: alecjacobson

repo: gptoolbox


text:
function drag(varargin)
%DRAG   Move around on a 2-D or 3-D plot.
%   DRAG with no arguments toggles the plot dragging state.
%   DRAG ON turns plot dragging on for the current figure.
%   DRAG XON, YON, or ZON turns plot dragging on for the x, y, or z axis only.
%   DRAG OFF turns plot dragging off in the current figure.
%   DRAG RESET resets the restore point to the current axis settings.
%
%   When plot dragging is on, you can click and drag the axes around while
%   maintaining the current level of zoom.  If DRAG RESET has not been
%   called the restore point is the original non-explored plot.  If DRAG
%   RESET has been called the restore point is the retore point that existed
%   when it was called.  Double clicking retores the axes to the current
%   restore point - the point at which plot dragging was first turned on for
%   this figure (or to the state to which the restore point was set by DRAG
%   RESET).
%
%   In a figure with multiple subplots, pressing and holding the left mouse
%   button drags around the currently active subplot, while pressing and
%   holding the right mouse button drags around all subplots simultaneously.
%   Double clicking with the left mouse button restores the currently active
%   plot to its restore point, while double clicking with the right mouse
%   button restores all subplots to their defined restore point.
%
%   DRAG(FIG,OPTION) applies the drag command to the figure specified by
%   FIG. OPTION can be any of the above arguments.
%
%   Example:
%   x = 0:0.1:5; y = exp(x);
%   figure; subplot(211); plot(x,y,'r.'); subplot(212); plot(y,x,'g.');
%   drag xon;
%
%   Ryan M. Eustice 08-15-2003
%   Woods Hole Oceanographic Institution, MS 7
%   Woods Hole, MA 02543
%   508.289.3269   ryan@whoi.edu

% History
% DATE          WHO                      WHAT
%----------    ------------------        ------------------------------
%2003-07-23    Ryan Eustice              Created & written.
%2003-08-15    Boyko Stoimenov           Modified the private function
%                                        dragreset to correctly count the
%                                        number of axes.  Also, switched
%                                        short-circuit && to & for matlab
%                                        backwards compatability.
%2004-01-29    Ryan Eustice              Fixed bug in private function
%                                        dragreset to properly count
%                                        number of axes and not count
%                                        legend or colorbar objects.

%   Note: drag uses the figure buttondown and buttonmotion functions

%
% PARSE ARGS - set fig and drag command
%
fignum = get(0,'CurrentFigure');
if isempty(fignum)
  return; % no figure
end

switch nargin
case 0 % no arg in
 switch get(fignum,'Tag')
 case 'dragon'  % toggle state
  dragoff(fignum);
 otherwise  % turn on
  figstate.motion = 'all';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 end % switch get(fignum,'Tag')
case 1 % one arg in
 switch lower(varargin{1})
 case 'on'
  figstate.motion = 'all';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'off'
  dragoff(fignum);
 case 'reset'
  dragreset(fignum);
 case 'xon'
  figstate.motion = 'xon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'yon'
  figstate.motion = 'yon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'zon'
  figstate.motion = 'zon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 otherwise
  error(sprintf('Unrecognized argument ''%s''',varargin{1}));
 end % switch varargin{1}
case 2 % two args in
 fignum = varargin{1};
 switch lower(varargin{2})
 case 'on'
  figstate.motion = 'all';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'off'
  dragoff(fignum);
 case 'reset'
  dragreset(fignum);
 case 'xon'
  figstate.motion = 'xon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'yon'
  figstate.motion = 'yon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 case 'zon'
  figstate.motion = 'zon';
  set(fignum,'UserData',figstate);
  dragon(fignum);
 otherwise
  error(sprintf('Unrecognized argument ''%s''',varargin{1}));
 end % switch varargin{2}
otherwise % too many args
 error(nargchk(0,2,nargin));
end % switch nargin

%====================================================================
% drag: private
%====================================================================
function dragon(fignum)
figstate = get(fignum,'UserData');
dragreset(fignum);
set(fignum,'Interruptible','on');
set(fignum,'DoubleBuffer','on');
set(fignum,'Tag','dragon');
set(fignum,'WindowButtonDownFcn',@startexploring);
set(fignum,'WindowButtonUpFcn',@stopexploring);

function dragoff(fignum)
set(fignum,'DoubleBuffer','off');
set(fignum,'Tag','');
set(fignum,'WindowButtonDownFcn','');
set(fignum,'WindowButtonUpFcn','');

function dragreset(fignum)
figstate = get(fignum,'UserData');
kids = get(fignum,'Children');
cc = 1;
for ii=1:length(kids)
  % check for regular axes handles by keying off the Tag field which should be
  % empty.  note that objects like legend and colorbar fill in the Tag field.
  if strcmp(get(kids(ii),'Type'),'axes') & strcmp(get(kids(ii),'Tag'),'')
    figstate.axishandle(cc) = kids(ii);
    figstate.origlimits(cc,:) = axis(kids(ii));
    cc = cc+1;
  end
end
figstate.numaxes = length(figstate.axishandle);
set(fignum,'UserData',figstate);

%====================================================================
% callbacks
%====================================================================
function startexploring(obj,eventdata)
fignum = gcbf;
figstate = get(fignum,'UserData');
switch get(fignum,'SelectionType')
case 'normal' %mouse left single-click
 figstate.whichaxes = 'current'; % only drag current axis
 axishandle = gca;
 initaxis(axishandle);
case 'open'   %mouse left or right double-click
 axishandle = gca;
 for ii=1:figstate.numaxes
   if strcmp(figstate.whichaxes,'current') & figstate.axishandle(ii) == axishandle
     % restore currently active axis limits to orig state
     axis(axishandle,figstate.origlimits(ii,:));
     break;
   elseif strcmp(figstate.whichaxes,'all')
     % restore other subplot axis limits to orig state
     axis(figstate.axishandle(ii),figstate.origlimits(ii,:));
   end
 end
case 'alt'    %mouse right single-click
 figstate.whichaxes = 'all'; % drag all axes simultaneously
 for ii=1:figstate.numaxes
   axishandle = figstate.axishandle(ii);
   initaxis(axishandle);
 end
otherwise
  error('Mouse SelectionType unknown');
end
set(fignum,'WindowButtonMotionFcn',@setnewview);
set(fignum,'Pointer','Fleur');
set(fignum,'UserData',figstate); % store figure state

function setnewview(obj,eventdata)
persistent cc;
fignum = gcbf;
figstate = get(fignum,'UserData');
switch figstate.whichaxes
case 'current'
 axishandle = gca;
 updateaxes(axishandle,figstate.motion);
case 'all'
 for ii=1:figstate.numaxes
   axishandle = figstate.axishandle(ii);
   updateaxes(axishandle,figstate.motion);
 end
otherwise
 %do nothing
end % switch figstate.whichaxes
cc = cc+1;
if mod(cc,2)
  drawnow; % force screen update every other motion cycle (lowers cpu usage)
end

function stopexploring(obj,eventdata)
fignum = gcbf;
set(fignum,'WindowButtonMotionFcn','');
set(fignum,'Pointer','arrow');

%====================================================================
% callbacks: private
%====================================================================
function pos = getcurrentposition(axishandle)
pos = get(axishandle,'CurrentPoint');
pos = mean(pos);
pos = round(pos*1000)/1000;

function initaxis(axishandle)
state.currentaxislimits = axis(axishandle);
state.startpoint = getcurrentposition(axishandle);
set(axishandle,'UserData',state); % store axis state

function updateaxes(axishandle,motion)
state = get(axishandle,'UserData');
state.endpoint = getcurrentposition(axishandle);
motionvector = state.startpoint - state.endpoint;

if length(state.currentaxislimits) == 4 % 2D plot
  switch motion
  case 'all'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2) + motionvector(1), ...
	state.currentaxislimits(3:4) + motionvector(2)];
  case 'xon'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2) + motionvector(1), ...
	state.currentaxislimits(3:4)];
  case 'yon'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2), ...
	state.currentaxislimits(3:4) + motionvector(2)];
  otherwise %do nothing
  end % switch motion
elseif length(state.currentaxislimits) == 6 % 3D plot
  switch motion
  case 'all'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2) + motionvector(1), ...
	state.currentaxislimits(3:4) + motionvector(2), ...
	state.currentaxislimits(5:6) + motionvector(3)];
  case 'xon'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2) + motionvector(1), ...
	state.currentaxislimits(3:4), ...
	state.currentaxislimits(5:6)];
  case 'yon'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2), ...
	state.currentaxislimits(3:4) + motionvector(2), ...
	state.currentaxislimits(5:6)];
  case 'zon'
   state.currentaxislimits = ...
       [state.currentaxislimits(1:2), ...
	state.currentaxislimits(3:4), ...
	state.currentaxislimits(5:6) + motionvector(3)];
  otherwise %do nothing
  end % switch motion
end
axis(axishandle,state.currentaxislimits);
set(axishandle,'UserData',state);


text:
function x = reErr(ref, sig)
%
% psnr -- Compute relative error for images
%
% Usage:
%       x = reErr(ref, sig)
%
% Input:
%       ref         Reference image
%       sig         Modified image
%
% Output:
%       x           reErr value
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe
%   Brendt Wohlberg   brendt@tmail.lanl.gov
%

[Nrows Ncols depth] = size(ref);

num = mean((ref(:)-sig(:)).^2);
den = mean((sig(:)).^2);

x = sqrt(num/den);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% GAINSANDREFERENCES compute gains matrices, references and other control
%                    related quantities for each state of the state machine.

%% --- Initialization ---

% CoM gains
Gain.KP_CoM = [50    100  5  % state ==  1  TWO FEET BALANCING
               50    100  5  % state ==  2  COM TRANSITION TO LEFT
               50    100  5  % state ==  3  LEFT FOOT BALANCING
               50    100  5  % state ==  4  YOGA LEFT FOOT
               50    100  5  % state ==  5  PREPARING FOR SWITCHING
               50    100  5  % state ==  6  LOOKING FOR CONTACT
               50    100  5  % state ==  7  TRANSITION TO INITIAL POSITION
               50    150  5  % state ==  8  COM TRANSITION TO RIGHT FOOT
               50    100  5  % state ==  9  RIGHT FOOT BALANCING
               50    100  5  % state == 10  YOGA RIGHT FOOT
               50    100  5  % state == 11  PREPARING FOR SWITCHING
               50    100  5  % state == 12  LOOKING FOR CONTACT
               50    100  5];% state == 13  TRANSITION TO INITIAL POSITION

Gain.KD_CoM = 2*sqrt(Gain.KP_CoM)/15;

% Angular momentum gains
Gain.KI_AngularMomentum = 3 ;
Gain.KP_AngularMomentum = 2*sqrt(Gain.KI_AngularMomentum)/5;

% Postural task gains
%                   %  TORSO  %%        LEFT ARM   %%       RIGHT ARM   %%         LEFT LEG           %%         RIGHT LEG           %%
Gain.KP_postural = [10   30   20, 10   10    10    8, 10   10    10    8, 30   30   20    20    100 100, 30   50   30    60    100 100  % state ==  1  TWO FEET BALANCING
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   30   20    20    100 100, 30   50   30    60    100 100  % state ==  2  COM TRANSITION TO LEFT
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   30    60    100 100, 30   30   20    20    100 100  % state ==  3  LEFT FOOT BALANCING
                    30   30   30, 10   10    10   10, 10   10    10   10,100  200  100   400    100 100,100   50   30    50    100 100  % state ==  4  YOGA LEFT FOOT
                    30   30   30,  5    5    10   10, 10   10    20   10,200  250   20    20     10  10,220  550  220   200     65 300  % state ==  5  PREPARING FOR SWITCHING
                    30   30   30, 10   10    20   10, 10   10    20   10,100  350   20   200     10 100,220  550  220   200     65 300  % state ==  6  LOOKING FOR CONTACT
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   60    30      5   5, 30   30   30    20      5   5  % state ==  7  TRANSITION TO INITIAL POSITION
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   60    30    100 100, 30   30   30    20    100 100  % state ==  8  COM TRANSITION TO RIGHT FOOT
                    10   30   20, 10   10    10    8, 10   10    10    8, 30   50   30    60    100 100, 30   30   20    20    100 100  % state ==  9  RIGHT FOOT BALANCING
                    30   30   30, 10   10    10   10, 10   10    10   10,100   50   30    50    100 100,100  200  100   100     10  10  % state == 10  YOGA RIGHT FOOT
                    30   30   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,200  250   20    20     10  10  % state == 11  PREPARING FOR SWITCHING
                    30   30   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,100  350   20   200     10 100  % state == 12  LOOKING FOR CONTACT
                    30   40   30, 10   10    10   10, 10   10    10   10,220  550  220   200     65 300,100  350   20   200     10 100];% state == 13  TRANSITION TO INITIAL POSITION

Gain.KD_postural = 0*sqrt(Gain.KP_postural(1,:))/20;

% symmetric gains
Gain.KP_postural(4,1:3)   = Gain.KP_postural(4,1:3)*3;
Gain.KP_postural(6,18:23) = Gain.KP_postural(6,18:23)*2;
Gain.KP_postural(10,1:3)  = Gain.KP_postural(10,1:3)*3;
Gain.KP_postural(13,1:3)  = Gain.KP_postural(13,1:3)*3;

%% Smoothing times

% Smoothing time gain scheduling
Config.SmoothingTimeGainScheduling = 2;

% Smoothing time CoM references
StateMachine.CoMSmoothingTime      = [1;    %% state ==  1  TWO FEET BALANCING
                                      1;    %% state ==  2  COM TRANSITION TO LEFT FOOT
                                      1;    %% state ==  3  LEFT FOOT BALANCING
                                      0.9;  %% state ==  4  YOGA LEFT FOOT
                                      2;    %% state ==  5  PREPARING FOR SWITCHING
                                      2;    %% state ==  6  LOOKING FOR CONTACT
                                      1;    %% state ==  7  TRANSITION INIT POSITION
                                      1;    %% state ==  8  COM TRANSITION TO RIGHT FOOT
                                      1;    %% state ==  9  RIGHT FOOT BALANCING
                                      0.9;  %% state == 10  YOGA RIGHT FOOT
                                      2;    %% state == 11  PREPARING FOR SWITCHING
                                      2;    %% state == 12  LOOKING FOR CONTACT
                                      5];   %% state == 13  TRANSITION INIT POSITION

% Smoothing time for joints references
StateMachine.jointsSmoothingTime   = [1;    %% state ==  1  TWO FEET BALANCING
                                      1;    %% state ==  2  COM TRANSITION TO LEFT FOOT
                                      1;    %% state ==  3  LEFT FOOT BALANCING
                                      0.9;  %% state ==  4  YOGA LEFT FOOT
                                      2;    %% state ==  5  PREPARING FOR SWITCHING
                                      2;    %% state ==  6  LOOKING FOR CONTACT
                                      1;    %% state ==  7  TRANSITION INIT POSITION
                                      1;    %% state ==  8  COM TRANSITION TO RIGHT FOOT
                                      1;    %% state ==  9  RIGHT FOOT BALANCING
                                      0.9;  %% state == 10  YOGA RIGHT FOOT
                                      2;    %% state == 11  PREPARING FOR SWITCHING
                                      2;    %% state == 12  LOOKING FOR CONTACT
                                      5];   %% state == 13  TRANSITION INIT POSITION

% scale factor smoothing time multiplies the smoothing factor during the
% Yoga (state 4 and 10). The purpose is to reduce the time necessary for
% the reference to converge to the next position, but without changing also
% the valuse stored in Sm.joints_leftYogaRef/Sm.joints_rightYogaRef
StateMachine.scaleFactorSmoothingTime = 0.9;

%% CoM delta

% To be summed to the reference CoM position
StateMachine.CoM_delta  = [% THIS REFERENCE IS USED AS A DELTA W.R.T. THE POSITION OF THE LEFT FOOT
                           0.0,  0.00,  0.0;   %% NOT USED
                           0.0,  0.00,  0.0;   %% state ==  2  COM TRANSITION TO LEFT FOOT
                           0.0,  0.005, 0.0;   %% state ==  3  LEFT FOOT BALANCING
                           0.0,  0.005, 0.0;   %% state ==  4  YOGA LEFT FOOT
                           0.0,  0.005, 0.0;   %% state ==  5  PREPARING FOR SWITCHING
                           0.02,-0.09,  0.0;   %% state ==  6  LOOKING FOR CONTACT
                           0.0,  0.00,  0.0;   %% state ==  7  TWO FEET BALANCING
                           % THIS REFERENCE IS USED AS A DELTA W.R.T. THE POSITION OF THE RIGHT FOOT
                           0.0,  0.00,  0.0;   %% state ==  8  COM TRANSITION TO RIGHT FOOT
                           0.0, -0.005,  0.0;  %% state ==  9  RIGHT FOOT BALANCING
                           0.0, -0.005, 0.0;   %% state == 10  YOGA RIGHT FOOT
                           0.0, -0.015, 0.0;   %% state == 11  PREPARING FOR SWITCHING
                           0.02, 0.02,  0.0;   %% state == 12  LOOKING FOR CONTACT
                           0.0,  0.00,  0.0];  %% NOT USED

%% Joint references
StateMachine.joints_references = [  zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 2  COM TRANSITION TO LEFT
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151];     %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 3  LEFT FOOT BALANCING
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,-0.1060,-0.1151];    %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 5  PREPARING FOR SWITCHING
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630, ...   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151];     %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 6  LOOKING FOR CONTACT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369, ...  %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277];      %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 8  COM TRANSITION TO RIGHT FOOT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369, ...  %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277];      %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 9  RIGHT FOOT BALANCING
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151, ...  %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630];      %
                                   zeros(1,ROBOT_DOF);                                %% THIS REFERENCE IS IGNORED
                                 [-0.0348,0.0779,0.0429, ...                          %% state == 11  PREPARING FOR SWITCHING
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                  -0.1493,0.8580,0.2437,0.8710, ...                   %
                                   0.0005,0.0793,-0.0014,-0.0051,0.0073,-0.1151, ...  %
                                  -0.0015,-0.1109,-0.0001,0.0003,0.0160,0.1630];      %
                                 [ 0.0864,0.0258,0.0152, ...                          %% state == 12  LOOKING FOR CONTACT
                                   0.1253,0.8135,0.3051,0.7928, ...                   %
                                   0.0563,0.6789,0.3340,0.6214, ...                   %
                                  -0.0026,0.0225,0.0093,-0.0020,0.0027,-0.0277, ...   %
                                   0.0107,-0.0741,-0.0001,-0.0120,0.0252,0.1369];     %
                                   zeros(1,ROBOT_DOF)];                               %% THIS REFERENCE IS IGNORED

% YOGA MOVESET (joint references during state 4 and 10)
q1 =        [-0.0790,0.2279, 0.4519, ...
             -1.1621,0.6663, 0.4919, 0.9947, ...
             -1.0717,1.2904,-0.2447, 1.0948, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3484,0.4008,-0.0004,-0.3672,-0.1060,-0.0875];

q2 =        [-0.0790,0.1279, 0.4519, ...
             -1.1621,0.6663, 0.4965, 0.9947, ...
             -1.0717,1.2904,-0.2493, 1.0948, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044,0.0700, ...
              0.3714,0.9599, 1.3253,-1.6594,-0.1060,-0.0614];

q3 =        [-0.0852,-0.3273,0.0821,...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092,0.2060, 0.0006,-0.1741,-0.1044,0.0700, ...
              0.3714,0.9599, 1.3253,-1.6594, 0.5000,-0.0614];

q4 =        [-0.0852,-0.4273,0.0821,...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.3473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q5 =        [-0.0790,-0.2273, 0.4519, ...
             -1.1621,0.6663, 0.4965, 0.9947, ...
             -1.0717,1.2904,-0.2493, 1.0948, ...
              0.2092, 0.4473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q6 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q7 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253, -1.6217, 0.5000,-0.0614];

q8 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q9 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 0.0107,1.3253,-0.0189, 0.5000,-0.0614];

q10 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q11 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q12 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 0.3107,1.3253,-0.0189, 0.5000,-0.0614];

q13 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

q14 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.8514, 0.0107,1.3253,-0.0189, 0.5000,-0.0614];

q15 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              1.5514, 0.3107,1.3253,-0.0189, 0.5000,-0.0614];

q16 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.2514, 0.0107,1.3253,-0.0189, 0.5000,-0.0614];

q17 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
             -0.3514, 0.3107,1.3253,-0.0189, 0.5000,-0.0614];

q18 =        [-0.0852,-0.4273,0.0821, ...
              0.1391, 1.4585,0.2464, 0.3042, ...
             -0.4181, 1.6800,0.7373, 0.3031, ...
              0.2092, 0.6473,0.0006,-0.1741,-0.1044, 0.0700, ...
              0.3514, 1.3107,1.3253,-0.0189, 0.5000,-0.0614];

StateMachine.joints_leftYogaRef  = [ 0,                                    q1;
                                     1*StateMachine.jointsSmoothingTime(4),q2;
                                     2*StateMachine.jointsSmoothingTime(4),q3;
                                     3*StateMachine.jointsSmoothingTime(4),q4;
                                     4*StateMachine.jointsSmoothingTime(4),q5;
                                     5*StateMachine.jointsSmoothingTime(4),q6;
                                     6*StateMachine.jointsSmoothingTime(4),q7;
                                     7*StateMachine.jointsSmoothingTime(4),q8;
                                     8*StateMachine.jointsSmoothingTime(4),q9;
                                     9*StateMachine.jointsSmoothingTime(4),q10;
                                    10*StateMachine.jointsSmoothingTime(4),q11;
                                    11*StateMachine.jointsSmoothingTime(4),q12;
                                    12*StateMachine.jointsSmoothingTime(4),q13;
                                    13*StateMachine.jointsSmoothingTime(4),q14;
                                    14*StateMachine.jointsSmoothingTime(4),q15;
                                    15*StateMachine.jointsSmoothingTime(4),q16;
                                    16*StateMachine.jointsSmoothingTime(4),q17;
                                    17*StateMachine.jointsSmoothingTime(4),q10;
                                    18*StateMachine.jointsSmoothingTime(4),q11;
                                    19*StateMachine.jointsSmoothingTime(4),q12;
                                    20*StateMachine.jointsSmoothingTime(4),q13;
                                    21*StateMachine.jointsSmoothingTime(4),q14;
                                    22*StateMachine.jointsSmoothingTime(4),q15;
                                    23*StateMachine.jointsSmoothingTime(4),q16;
                                    24*StateMachine.jointsSmoothingTime(4),q17;
                                    25*StateMachine.jointsSmoothingTime(4),q18];

StateMachine.joints_rightYogaRef      = StateMachine.joints_leftYogaRef;
StateMachine.joints_rightYogaRef(:,1) = [0;
                                         1*StateMachine.jointsSmoothingTime(10);
                                         2*StateMachine.jointsSmoothingTime(10);
                                         3*StateMachine.jointsSmoothingTime(10);
                                         4*StateMachine.jointsSmoothingTime(10);
                                         5*StateMachine.jointsSmoothingTime(10);
                                         6*StateMachine.jointsSmoothingTime(10);
                                         7*StateMachine.jointsSmoothingTime(10);
                                         8*StateMachine.jointsSmoothingTime(10);
                                         9*StateMachine.jointsSmoothingTime(10);
                                        10*StateMachine.jointsSmoothingTime(10);
                                        11*StateMachine.jointsSmoothingTime(10);
                                        12*StateMachine.jointsSmoothingTime(10);
                                        13*StateMachine.jointsSmoothingTime(10);
                                        14*StateMachine.jointsSmoothingTime(10);
                                        15*StateMachine.jointsSmoothingTime(10);
                                        16*StateMachine.jointsSmoothingTime(10);
                                        17*StateMachine.jointsSmoothingTime(10);
                                        18*StateMachine.jointsSmoothingTime(10);
                                        19*StateMachine.jointsSmoothingTime(10);
                                        20*StateMachine.jointsSmoothingTime(10);
                                        21*StateMachine.jointsSmoothingTime(10);
                                        22*StateMachine.jointsSmoothingTime(10);
                                        23*StateMachine.jointsSmoothingTime(10);
                                        24*StateMachine.jointsSmoothingTime(10);
                                        25*StateMachine.jointsSmoothingTime(10)];

% if the demo is not "yogaExtended", stop at the 8th move
if ~StateMachine.yogaExtended

    StateMachine.joints_leftYogaRef       = StateMachine.joints_leftYogaRef(1:8,:);
    StateMachine.joints_rightYogaRef      = StateMachine.joints_rightYogaRef(1:8,:);
    StateMachine.joints_leftYogaRef(8,1)  = 15*StateMachine.jointsSmoothingTime(4);
    StateMachine.joints_rightYogaRef(8,1) = 15*StateMachine.jointsSmoothingTime(10);
end

% MIRROR YOGA LEFT MOVESET FOR RIGHT YOGA
for i = 1:size(StateMachine.joints_rightYogaRef,1)

	StateMachine.joints_rightYogaRef(i,2:4)           = [StateMachine.joints_rightYogaRef(i,2) -StateMachine.joints_rightYogaRef(i,3) -StateMachine.joints_rightYogaRef(i,4)];
	rightArm                                          =  StateMachine.joints_rightYogaRef(i,end-15:end-12);
	StateMachine.joints_rightYogaRef(i,end-15:end-12) =  StateMachine.joints_rightYogaRef(i,end-19:end-16);
	StateMachine.joints_rightYogaRef(i,end-19:end-16) =  rightArm;
	rightLeg                                          =  StateMachine.joints_rightYogaRef(i,end-5:end);
	StateMachine.joints_rightYogaRef(i,end-5:end)     =  StateMachine.joints_rightYogaRef(i,end-11:end-6);
	StateMachine.joints_rightYogaRef(i,end-11:end-6)  =  rightLeg;
end

%% References for CoM trajectory (COORDINATOR DEMO ONLY)

% that the robot waits before starting the left-and-right
Config.noOscillationTime       = 0;
Config.directionOfOscillation  = [0;1;0];
Config.amplitudeOfOscillation  = 0.02; % [m]
Config.frequencyOfOscillation  = 0.2;  % [Hz]

author: robotology

repo: whole-body-controllers


text:
% pop_chancoresp() - define correspondances between two channel locations structures
%                    (EEG.chanlocs) automatically (by matching channel labels)
%                    else using a user input gui.
% Usage:
%   >> [chanlist1 chanlist2] = pop_chancoresp(chanstruct1, chanstruc2, 'key', 'val', ...);
%
% Inputs:
%   chanstruct1     - first (new) channel locations structure (EEG.chanlocs).
%                      For details, >> help readlocs
%   chanstruct2     - second (reference) chanlocs structure.
%
% Optional parameters:
%   'gui'           - ['on'|'off'] display gui or not ('on' -> yes)
%   'autoselect'    - ['none'|'fiducials'|'all'] automatically pair channels
%   'chaninfo1'     - EEG.chaninfo structure for first (new) EEG.chanlocs
%   'chaninfo2'     - EEG.chaninfo structure for second (reference) EEG.chanlocs
%   'chanlist1'     - [integer] selected channel to pair in the graphic interface
%                     for the first channel structure. This requires the input
%                     'chanlist2' below.
%   'chanlist2'     - [integer] selected channel to pair in the graphic interface
%                     for the first channel structure. This requires the input
%                     requires the input 'chanlist1' that must be of the same length.
% Output:
%   chanlist1       - [int vector] indices of selected channels from first (new) EEG.chanlocs
%   chanlist2       - [int vector] selected channels from second (reference) EEG.chanlocs
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2005

% Copyright (C) 2005 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [chanlistout1, chanlistout2, thirdout, outfourth] = pop_chancoresp(chans1, chans2, varargin);

    if nargin < 2
        help pop_chancoresp;
        return;
    end;
    chanlistout1 = [];
    chanlistout2 = [];

    % process sub command
    % -------------------
    if isstr(chans1)
       if strcmpi(chans1, 'pair')
           [chanlistout1, chanlistout2, thirdout, outfourth] = pair(chans2, varargin{:});
       elseif strcmpi(chans1, 'unpair')
           [chanlistout1, chanlistout2, thirdout, outfourth] = unpair(chans2, varargin{:});
       elseif strcmpi(chans1, 'clear')
           [chanlistout1, chanlistout2, thirdout, outfourth] = clearchans(chans2, varargin{:});
       elseif strcmpi(chans1, 'auto')
           [chanlistout1, chanlistout2, thirdout, outfourth] = autoselect(chans2, varargin{:});
       end;
       return;
    end;

    g = finputcheck(varargin, { 'autoselect'     'string'  {'none' 'fiducials' 'all'}   'all';
                                'chanlist1'      'integer' [1 Inf] [];
                                'chanlist2'      'integer' [1 Inf] [];
                                'chaninfo1'      ''        []      [];
                                'chaninfo2'      ''        []      [];
                                'gui'            'string'  { 'on' 'off' }  'on' } );
    if isstr(g), error(g); end;
    g.chanstruct1 = chans1;
    g.chanstruct2 = chans2;
    if length(g.chanlist1) ~= length(g.chanlist2)
        error('input arguments ''chanlist1'' and ''chanlist2'' must have the same length');
    end;

    % decode different input formats
    % ------------------------------
    if isstruct(chans1)
        if isfield(chans1, 'label') % fieldtrip
            chanstr1 = chans1.label;
            chanstr2 = chans2.label;
        else % EEGLAB
            chanstr1 = { chans1.labels };
            chanstr2 = { chans2.labels };
        end;
    else % only channel labels
        chanstr1 = chans1;
        chanstr2 = chans2;
    end;

    % convert selection to integer
    % ----------------------------
    if isempty(g.chanlist1)
        if strcmpi(g.autoselect, 'fiducials')
            % find fiducials in both channel location strustures
            % --------------------------------------------------
            naz1 = strmatch('nz', lower( chanstr1 ), 'exact'); if isempty(naz1), naz1 = strmatch('nasion', lower( chanstr1 ), 'exact'); end;
            naz2 = strmatch('nz', lower( chanstr2 ), 'exact'); if isempty(naz2), naz2 = strmatch('nasion', lower( chanstr2 ), 'exact'); end;
            lpa1 = strmatch('lpa', lower( chanstr1 ), 'exact'); if isempty(lpa1), lpa1 = strmatch('left',  lower( chanstr1 ), 'exact'); end;
            lpa2 = strmatch('lpa', lower( chanstr2 ), 'exact'); if isempty(lpa2), lpa2 = strmatch('left',  lower( chanstr2 ), 'exact'); end;
            rpa1 = strmatch('rpa', lower( chanstr1 ), 'exact'); if isempty(rpa1), rpa1 = strmatch('right', lower( chanstr1 ), 'exact'); end;
            rpa2 = strmatch('rpa', lower( chanstr2 ), 'exact'); if isempty(rpa2), rpa2 = strmatch('right', lower( chanstr2 ), 'exact'); end;
            g.chanlist1 = [ naz1 lpa1 rpa1 ];
            g.chanlist2 = [ naz2 lpa2 rpa2 ];
            if length(g.chanlist1) ~= length(g.chanlist2) | length(g.chanlist1) == 0
                disp('Warning: could not find fiducials in at least one of the channel location structure');
                g.chanlist1 = [];
                g.chanlist2 = [];
            end;
        elseif strcmpi(g.autoselect, 'all')
            % find common channels in both channel location strustures
            % --------------------------------------------------------
            chanstr2low = lower( chanstr2 );
            chanstr1low = lower( chanstr1 );
            for index = 1:length( chanstr1 )
                ind = strmatch(chanstr1low{index}, chanstr2low, 'exact' );
                if ~isempty(ind)
                    g.chanlist1(end+1) = index;
                    g.chanlist2(end+1) = ind;
                end;
            end;
        end;
    end;

    % plot
    % ----
    if strcmpi(g.gui, 'off')
        chanlistout1 = g.chanlist1;
        chanlistout2 = g.chanlist2;
        return;
    end;
    try,  g.promptstring;  catch, g.promptstring = ''; end;
    try,  g.selectionmode; catch, g.selectionmode = 'multiple'; end;
    try,  g.listsize;      catch, g.listsize = []; end;
    try,  g.initialvalue;  catch, g.initialvalue = []; end;
    try,  g.name;          catch, g.name = ''; end;
	g.chanstr1 = chanstr1;
	g.chanstr2 = chanstr2;

    fig = figure('visible', 'off');
    set(fig, 'name', 'Select corresponding channels to pair');

    % make text for list
    % ------------------
    [ g.newchanstr1 g.newchanstr2 ] = makelisttext( chanstr1, chanstr2, g.chanlist1, g.chanlist2);

    % callback for paring and unpairing
    % ---------------------------------
    cb_pair = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                'tmpval1 = get(findobj(gcbf, ''tag'', ''list1''), ''value'');' ...
                'tmpval2 = get(findobj(gcbf, ''tag'', ''list2''), ''value'');' ...
                '[tmpdat.newchanstr1{tmpval1}, tmpdat.newchanstr2{tmpval2}, tmpdat.chanlist1, tmpdat.chanlist2] = pop_chancoresp(''pair'', tmpval1, tmpval2, tmpdat.chanstr1, tmpdat.chanstr2, tmpdat.chanlist1, tmpdat.chanlist2, tmpdat.newchanstr1{tmpval1}, tmpdat.newchanstr2{tmpval2});' ...
                'set(findobj(gcbf, ''tag'', ''list1''), ''string'', tmpdat.newchanstr1);' ...
                'set(findobj(gcbf, ''tag'', ''list2''), ''string'', tmpdat.newchanstr2);' ...
                'set(gcbf, ''userdata'', tmpdat);' ...
                'clear tmpdat tmpval1 tmpval2;' ];
    cb_unpair = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                'tmpval1 = get(findobj(gcbf, ''tag'', ''list1''), ''value'');' ...
                'tmpval2 = get(findobj(gcbf, ''tag'', ''list2''), ''value'');' ...
                '[tmpdat.newchanstr1{tmpval1}, tmpdat.newchanstr2{tmpval2}, tmpdat.chanlist1, tmpdat.chanlist2] = pop_chancoresp(''unpair'', tmpval1, tmpval2, tmpdat.chanstr1, tmpdat.chanstr2, tmpdat.chanlist1, tmpdat.chanlist2, tmpdat.newchanstr1{tmpval1}, tmpdat.newchanstr2{tmpval2});' ...
                'set(findobj(gcbf, ''tag'', ''list1''), ''string'', tmpdat.newchanstr1);' ...
                'set(findobj(gcbf, ''tag'', ''list2''), ''string'', tmpdat.newchanstr2);' ...
                'set(gcbf, ''userdata'', tmpdat);' ...
                'clear tmpdat tmpval1 tmpval2;' ];
    cb_plot1 = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                 'figure; topoplot([], tmpdat.chanstruct1, ''style'', ''blank'', ''drawaxis'', ''on'', ' ...
                 '''electrodes'', ''labelpoint'', ''chaninfo'', tmpdat.chaninfo1);' ];
    cb_plot2 = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                 'figure; topoplot([], tmpdat.chanstruct2, ''style'', ''blank'', ''drawaxis'', ''on'', ' ...
                 '''electrodes'', ''labelpoint'', ''chaninfo'', tmpdat.chaninfo2);' ];
    cb_list1 = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                'tmpval = get(findobj(gcbf, ''tag'', ''list1''), ''value'');' ...
                'tmppos = find(tmpdat.chanlist1 == tmpval);' ...
                'if ~isempty(tmppos), set(findobj(gcbf, ''tag'', ''list2''), ''value'', tmpdat.chanlist2(tmppos)); end;' ...
                'clear tmpdat tmpval tmppos;' ];
    cb_list2 = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                'tmpval = get(findobj(gcbf, ''tag'', ''list2''), ''value'');' ...
                'tmppos = find(tmpdat.chanlist2 == tmpval);' ...
                'if ~isempty(tmppos), set(findobj(gcbf, ''tag'', ''list1''), ''value'', tmpdat.chanlist1(tmppos)); end;' ...
                'clear tmpdat tmpval tmppos;' ];
    cb_clear = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                 '[tmpdat.newchanstr1, tmpdat.newchanstr2, tmpdat.chanlist1, tmpdat.chanlist2] = pop_chancoresp(''clear'', tmpdat.chanstr1, tmpdat.chanstr2);' ...
                 'set(gcbf, ''userdata'', tmpdat);' ...
                 'set(findobj(gcbf, ''tag'', ''list1''), ''string'', tmpdat.newchanstr1);' ...
                 'set(findobj(gcbf, ''tag'', ''list2''), ''string'', tmpdat.newchanstr2);' ...
                 'set(gcbf, ''userdata'', tmpdat);' ...
                 'clear tmpdat;' ];
    cb_auto  = [ 'tmpdat = get(gcbf, ''userdata'');' ...
                 '[tmpdat.newchanstr1, tmpdat.newchanstr2, tmpdat.chanlist1, tmpdat.chanlist2] = pop_chancoresp(''auto'', tmpdat.chanstr1, tmpdat.chanstr2, tmpdat.newchanstr1, tmpdat.newchanstr2, tmpdat.chanlist1, tmpdat.chanlist2);' ...
                 'set(gcbf, ''userdata'', tmpdat);' ...
                 'set(findobj(gcbf, ''tag'', ''list1''), ''string'', tmpdat.newchanstr1);' ...
                 'set(findobj(gcbf, ''tag'', ''list2''), ''string'', tmpdat.newchanstr2);' ...
                 'set(gcbf, ''userdata'', tmpdat);' ...
                 'clear tmpdat;' ];
    geometry = {[1 1] [1 1] [1 1] [1 1] [1 1]};
    geomvert = [ 1 min(max(length(chanstr1),length(chanstr2)), 10) 1 1 1];
    listui = { ...
              { 'Style', 'pushbutton', 'string', 'Plot new montage', 'callback', cb_plot1 }  ...
              { 'Style', 'pushbutton', 'string', 'Plot ref montage', 'callback', cb_plot2 }  ...
              { 'Style', 'listbox', 'tag', 'list1', 'string', g.newchanstr1, 'value', 1, 'min', 1, 'max', 2, 'callback', cb_list1 } ...
              { 'Style', 'listbox', 'tag', 'list2', 'string', g.newchanstr2, 'value', 1, 'min', 1, 'max', 2, 'callback', cb_list2 } ...
              { 'Style', 'pushbutton', 'string', 'Pair channels'  , 'callback', cb_pair } ...
              { 'Style', 'pushbutton', 'string', 'Clear this pair', 'callback', cb_unpair } ...
              { 'Style', 'pushbutton', 'string', 'Clear all pairs' , 'callback', cb_clear } ...
              { 'Style', 'pushbutton', 'string', 'Auto select', 'callback', cb_auto } ...
              { 'Style', 'pushbutton', 'string', 'Cancel', 'callback', 'close(gcbf);' }  ...
              { 'Style', 'pushbutton', 'string', 'Ok'    , 'tag', 'ok', 'callback', ['set(gcbo, ''userdata'', ''ok'');'] } };

    [tmp tmp2 allobj] = supergui( fig, geometry, geomvert, listui{:} );
    set(fig, 'userdata', g);

    % decode output
    % -------------
    okbut = findobj( 'parent', fig, 'tag', 'ok');
    figure(fig);
    drawnow;
    waitfor( okbut, 'userdata');
    try,
        tmpdat = get(fig, 'userdata');
        chanlistout1 = tmpdat.chanlist1;
        chanlistout2 = tmpdat.chanlist2;
        close(fig);
        drawnow;
    end;

% unpair channels
% ---------------
function [ str1, str2, chanlist1, chanlist2 ] = unpair(ind1, ind2, chanstr1, chanstr2, chanlist1, chanlist2, str1, str2);
    if nargin > 4
        if isempty(find(chanlist2 == ind2)), disp('Channels not associated'); return; end;
        if isempty(find(chanlist1 == ind1)), disp('Channels not associated'); return; end;
    end;

    str1 = sprintf('%2d - %3s', ind1, chanstr1{ind1});
    str2 = sprintf('%2d - %3s', ind2, chanstr2{ind2});
    if nargout > 2
        tmppos = find( chanlist1 == ind1);
        chanlist1(tmppos) = [];
        chanlist2(tmppos) = [];
    end;

% pair channels
% -------------
function [ str1, str2, chanlist1, chanlist2 ] = pair(ind1, ind2, chanstr1, chanstr2, chanlist1, chanlist2, str1, str2);
    if nargin > 4
        if ~isempty(find(chanlist2 == ind2)), disp('Channel in second structure already associated'); return; end;
        if ~isempty(find(chanlist1 == ind1)), disp('Channel in first structure already associated'); return; end;
    end;

    str1 = sprintf('%2d - %3s   -> %2d - %3s', ind1, chanstr1{ind1}, ind2, chanstr2{ind2});
    str2 = sprintf('%2d - %3s   -> %2d - %3s', ind2, chanstr2{ind2}, ind1, chanstr1{ind1});
    if nargout > 2
        chanlist1 = [ chanlist1 ind1 ];
        chanlist2 = [ chanlist2 ind2 ];
    end;

% make full channel list
% ----------------------
function [ newchanstr1, newchanstr2 ] = makelisttext( chanstr1, chanstr2, chanlist1, chanlist2);
    for index = 1:length(chanstr1)
        if ismember(index, chanlist1)
            pos = find(chanlist1 == index);
            newchanstr1{index} = pair( chanlist1(pos), chanlist2(pos), chanstr1, chanstr2 );
        else
            newchanstr1{index} = sprintf('%2d - %3s', index, chanstr1{index});
        end;
    end;
    for index = 1:length(chanstr2)
        if ismember(index, chanlist2)
            pos = find(chanlist2 == index);
            [tmp newchanstr2{index}] = pair( chanlist1(pos), chanlist2(pos), chanstr1, chanstr2 );
        else
            newchanstr2{index} = sprintf('%2d - %3s', index, chanstr2{index});
        end;
    end;

% clear channel pairs
% -------------------
function [ newchanstr1, newchanstr2, chanlist1, chanlist2 ] = clearchans(chanstr1, chanstr2);

    chanlist1 = [];
    chanlist2 = [];
    [ newchanstr1, newchanstr2 ] = makelisttext( chanstr1, chanstr2, chanlist1, chanlist2);

% autoselect channel pairs
% ------------------------
function [ newchanstr1, newchanstr2, chanlist1, chanlist2 ] = autoselect(chanstr1, chanstr2, newchanstr1, newchanstr2, chanlist1, chanlist2);

    % GUI for selecting pairs
    % -----------------------
    listoptions = { 'All channels (same labels)' 'Fiducials (same labels)' 'BIOSEMI -> 10-20' ...
                    'EGI -> 10-20' };
    listui = { { 'style' 'text' 'string' [ 'How to pair channels (click to select)' 10 ] } ...
               { 'style' 'listbox' 'string' listoptions 'value' 1 } };
    results = inputgui({ [1 1] }, listui, '', 'Auto select channel pairs', [], 'normal', 2);

    % decode results
    % --------------
    if isempty(results), return; end;
    if results{1} == 1 % select all pairs
        [chanlist1 chanlist2] = pop_chancoresp(chanstr1, chanstr2, 'autoselect', 'all', 'gui', 'off');
    elseif results{1} == 2
        [chanlist1 chanlist2] = pop_chancoresp(chanstr1, chanstr2, 'autoselect', 'fiducials', 'gui', 'off');
    else
        disp('Not implemented yet'); return;
    end;
    [ newchanstr1, newchanstr2 ] = makelisttext( chanstr1, chanstr2, chanlist1, chanlist2);

author: PatternRecognition

repo: OpenBMI


text:
function make_train_lut

    hks_i = comm.KasamiSequence('SamplesPerFrame', 65,'Index',1);
    hks_q = comm.KasamiSequence('SamplesPerFrame', 65,'Index',3);
    x_i = step(hks_i);
    x_q = step(hks_q);
    t_i = x_i*2-1;
    t_q = x_q*2-1;

    fid = fopen('TB_i.m','w+');
    fprintf(fid,'function y = TB_i\n');
    fprintf(fid,'%%#codegen\n');
    fprintf(fid,'y = [\n');
    fprintf(fid,'%1d\n',t_i);
    fprintf(fid,'];\n');
    fclose(fid);
    fid = fopen('TB_q.m','w+');
    fprintf(fid,'function y = TB_q\n');
    fprintf(fid,'%%#codegen\n');
    fprintf(fid,'y = [\n');
    fprintf(fid,'%1d\n',t_q);
    fprintf(fid,'];\n');
    fclose(fid);
end


text:
function rez = learnAndSolve8b_old(rez)
% This is the main optimization. Takes the longest time and uses the GPU heavily.

ops = rez.ops;
ops.fig = getOr(ops, 'fig', 1); % whether to show plots every N batches

NrankPC = 6; % this one is the rank of the PCs, used to detect spikes with threshold crossings
Nrank = 3; % this one is the rank of the templates
rng('default'); rng(1);

% we need PC waveforms, as well as template waveforms
[wTEMP, wPCA]    = extractTemplatesfromSnippets(rez, NrankPC);

% move these to the GPU
wPCA = gpuArray(wPCA(:, 1:Nrank));
wTEMP = gpuArray(wTEMP);
wPCAd = double(wPCA); % convert to double for extra precision
ops.wPCA = gather(wPCA);
ops.wTEMP = gather(wTEMP);
nt0 = ops.nt0;
nt0min  = rez.ops.nt0min;
rez.ops = ops;
nBatches  = rez.temp.Nbatch;
NT  	= ops.NT;
Nfilt 	= ops.Nfilt;
Nchan 	= ops.Nchan;

% two variables for the same thing? number of nearest channels to each primary channel
NchanNear   = min(ops.Nchan, 32);
Nnearest    = min(ops.Nchan, 32);

% decay of gaussian spatial mask centered on a channel
sigmaMask  = ops.sigmaMask;

% spike threshold for finding missed spikes in residuals
ops.spkTh = -6; % why am I overwriting this here?

batchstart = 0:NT:NT*nBatches;

% find the closest NchanNear channels, and the masks for those channels
[iC, mask, C2C] = getClosestChannels(rez, sigmaMask, NchanNear);

% sorting order for the batches
isortbatches = rez.iorig(:);
nhalf = ceil(nBatches/2); % halfway point

% this batch order schedule goes through half of the data forward and backward during the model fitting
% and then goes through the data symmetrically-out from the center during the final pass
ischedule = [nhalf:nBatches nBatches:-1:nhalf];
i1 = [(nhalf-1):-1:1];
i2 = [nhalf:nBatches];

irounds = cat(2, ischedule, i1, i2);

niter   = numel(irounds);
if irounds(niter - nBatches)~=nhalf
    error('mismatch between number of batches'); % this check is in here in case I do somehting weird when I try different schedules
end

% these two flags are used to keep track of what stage of model fitting we're at
flag_final    =  0;
flag_resort   = 1;

% this is the absolute temporal offset in seconds corresponding to the start of the
% spike sorted time segment
t0 = ceil(rez.ops.trange(1) * ops.fs);

nInnerIter  = 60; % this is for SVD for the power iteration

% schedule of learning rates for the model fitting part
% starts small and goes high, it corresponds approximately to the number of spikes
% from the past that were averaged to give rise to the current template
pmi = exp(-1./linspace(ops.momentum(1), ops.momentum(2), niter-nBatches));

Nsum = min(Nchan,7); % how many channels to extend out the waveform in mexgetspikes
% lots of parameters passed into the CUDA scripts
Params     = double([NT Nfilt ops.Th(1) nInnerIter nt0 Nnearest ...
    Nrank ops.lam pmi(1) Nchan NchanNear ops.nt0min 1 Nsum NrankPC ops.Th(1)]);

% W0 has to be ordered like this
W0 = permute(double(wPCA), [1 3 2]);

% initialize the list of channels each template lives on
iList = int32(gpuArray(zeros(Nnearest, Nfilt)));

% initialize average number of spikes per batch for each template
nsp = gpuArray.zeros(0,1, 'double');

% this flag starts 0, is set to 1 later
Params(13) = 0;

% kernels for subsample alignment
[Ka, Kb] = getKernels(ops, 10, 1);

p1 = .95; % decay of nsp estimate in each batch

fprintf('Time %3.0fs. Optimizing templates ...\n', toc)

fid = fopen(ops.fproc, 'r');

ntot = 0;
ndrop = zeros(1,2); % this keeps track of dropped templates for debugging purposes

m0 = ops.minFR * ops.NT/ops.fs; % this is the minimum firing rate that all templates must maintain, or be dropped

for ibatch = 1:niter
    korder = irounds(ibatch); % korder is the index of the batch at this point in the schedule
    k = isortbatches(korder); % k is the index of the batch in absolute terms

    if ibatch>niter-nBatches && korder==nhalf
        % this is required to revert back to the template states in the middle of the batches
        [W, dWU] = revertW(rez);
        fprintf('reverted back to middle timepoint \n')
    end

    if ibatch<=niter-nBatches
        % obtained pm for this batch
        Params(9) = pmi(ibatch);
        pm = pmi(ibatch) * gpuArray.ones(Nfilt, 1, 'double');
    end

    % loading a single batch (same as everywhere)
    offset = 2 * ops.Nchan*batchstart(k);
    fseek(fid, offset, 'bof');
    dat = fread(fid, [NT ops.Nchan], '*int16');
    dataRAW = single(gpuArray(dat))/ ops.scaleproc;

    if ibatch==1
       % only on the first batch, we first get a new set of spikes from the residuals,
       % which in this case is the unmodified data because we start with no templates
        [dWU, cmap] = mexGetSpikes2(Params, dataRAW, wTEMP, iC-1); % CUDA function to get spatiotemporal clips from spike detections
        dWU = double(dWU);
        dWU = reshape(wPCAd * (wPCAd' * dWU(:,:)), size(dWU)); % project these into the wPCA waveforms

        W = W0(:,ones(1,size(dWU,3)),:); % initialize the low-rank decomposition with standard waves
        Nfilt = size(W,2); % update the number of filters/templates
        nsp(1:Nfilt) = m0; % initialize the number of spikes for new templates with the minimum allowed value, so it doesn't get thrown back out right away
        Params(2) = Nfilt; % update in the CUDA parameters
    end

    if flag_resort
        % this is a flag to resort the order of the templates according to best peak channel
        % this is important in order to have cohesive memory requests from the GPU RAM
        [~, iW] = max(abs(dWU(nt0min, :, :)), [], 2); % max channel (either positive or negative peak)
        iW = int32(squeeze(iW));

        [iW, isort] = sort(iW); % sort by max abs channel
        W = W(:,isort, :); % user ordering to resort all the other template variables
        dWU = dWU(:,:,isort);
        nsp = nsp(isort);
    end

    % decompose dWU by svd of time and space (via covariance matrix of 61 by 61 samples)
    % this uses a "warm start" by remembering the W from the previous iteration
    [W, U, mu] = mexSVDsmall2(Params, dWU, W, iC-1, iW-1, Ka, Kb);

    % UtU is the gram matrix of the spatial components of the low-rank SVDs
    % it tells us which pairs of templates are likely to "interfere" with each other
    % such as when we subtract off a template
    [UtU, maskU] = getMeUtU(iW, iC, mask, Nnearest, Nchan); % this needs to change (but I don't know why!)

    % main CUDA function in the whole codebase. does the iterative template matching
    % based on the current templates, gets features for these templates if requested (featW, featPC),
    % gets scores for the template fits to each spike (vexp), outputs the average of
    % waveforms assigned to each cluster (dWU0),
    % and probably a few more things I forget about
    [st0, id0, x0, featW, dWU0, drez, nsp0, featPC, vexp] = ...
        mexMPnu8(Params, dataRAW, single(U), single(W), single(mu), iC-1, iW-1, UtU, iList-1, ...
        wPCA);

    % Sometimes nsp can get transposed (think this has to do with it being
    % a single element in one iteration, to which elements are added
    % nsp, nsp0, and pm must all be row vectors (Nfilt x 1), so force nsp
    % to be a row vector.
    [nsprow, nspcol] = size(nsp);
    if nsprow<nspcol
        nsp = nsp';
    end

    % updates the templates as a running average weighted by recency
    % since some clusters have different number of spikes, we need to apply the
    % exp(pm) factor several times, and fexp is the resulting update factor
    % for each template
    fexp = exp(double(nsp0).*log(pm(1:Nfilt)));
    fexp = reshape(fexp, 1,1,[]);
    dWU = dWU .* fexp + (1-fexp) .* (dWU0./reshape(max(1, double(nsp0)), 1,1, []));

    % nsp just gets updated according to the fixed factor p1
    nsp = nsp * p1 + (1-p1) * double(nsp0);

    % \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

    if ibatch==niter-nBatches
        % if we reached this point, we need to disable secondary template updates
        % like dropping, and adding new templates. We need to memorize the state of the
        % templates at this timepoint, and set the processing mode to "extraction and tracking"

        flag_resort   = 0; % no need to resort templates by channel any more
        flag_final = 1; % this is the "final" pass

        % final clean up, triage templates one last time
        [W, U, dWU, mu, nsp, ndrop] = ...
            triageTemplates2(ops, iW, C2C, W, U, dWU, mu, nsp, ndrop);

        % final number of templates
        Nfilt = size(W,2);
        Params(2) = Nfilt;

        % final covariance matrix between all templates
        [WtW, iList] = getMeWtW(single(W), single(U), Nnearest);

        % iW is the final channel assigned to each template
        [~, iW] = max(abs(dWU(nt0min, :, :)), [], 2);
        iW = int32(squeeze(iW));

        % extract ALL features on the last pass
        Params(13) = 2; % this is a flag to output features (PC and template features)

        % different threshold on last pass?
        Params(3) = ops.Th(end); % usually the threshold is much lower on the last pass

        rez = memorizeW(rez, W, dWU, U, mu); % memorize the state of the templates
        fprintf('memorized middle timepoint \n')
    end

    if ibatch<niter-nBatches %-50
        % during the main "learning" phase of fitting a model
        if rem(ibatch, 5)==1
            % this drops templates based on spike rates and/or similarities to other templates
            [W, U, dWU, mu, nsp, ndrop] = ...
                triageTemplates2(ops, iW, C2C, W, U, dWU, mu, nsp, ndrop);
        end
        Nfilt = size(W,2); % update the number of filters
        Params(2) = Nfilt;

        % this adds new templates if they are detected in the residual
        [dWU0,cmap] = mexGetSpikes2(Params, drez, wTEMP, iC-1);

        if size(dWU0,3)>0
            % new templates need to be integrated into the same format as all templates
            dWU0 = double(dWU0);
            dWU0 = reshape(wPCAd * (wPCAd' * dWU0(:,:)), size(dWU0)); % apply PCA for smoothing purposes
            dWU = cat(3, dWU, dWU0);

            W(:,Nfilt + [1:size(dWU0,3)],:) = W0(:,ones(1,size(dWU0,3)),:); % initialize temporal components of waveforms

            nsp(Nfilt + [1:size(dWU0,3)]) = ops.minFR * NT/ops.fs; % initialize the number of spikes with the minimum allowed
            mu(Nfilt + [1:size(dWU0,3)])  = 10; % initialize the amplitude of this spike with a lowish number

            Nfilt = min(ops.Nfilt, size(W,2)); % if the number of filters exceed the maximum allowed, clip it
            Params(2) = Nfilt;

            W   = W(:, 1:Nfilt, :); % remove any new filters over the maximum allowed
            dWU = dWU(:, :, 1:Nfilt); % remove any new filters over the maximum allowed
            nsp = nsp(1:Nfilt); % remove any new filters over the maximum allowed
            mu  = mu(1:Nfilt); % remove any new filters over the maximum allowed
        end

    end

    if ibatch>niter-nBatches
        % during the final extraction pass, this keesp track of all spikes and features

        % we memorize the spatio-temporal decomposition of the waveforms at this batch
        % this is currently only used in the GUI to provide an accurate reconstruction
        % of the raw data at this time
        rez.WA(:,:,:,k) = gather(W);
        rez.UA(:,:,:,k) = gather(U);
        rez.muA(:,k) = gather(mu);

        % we carefully assign the correct absolute times to spikes found in this batch
        ioffset         = ops.ntbuff;
        if k==1
            ioffset         = 0; % the first batch is special (no pre-buffer)
        end
        toff = nt0min + t0 -ioffset + (NT-ops.ntbuff)*(k-1);
        st = toff + double(st0);

        irange = ntot + [1:numel(x0)]; % spikes and features go into these indices

        if ntot+numel(x0)>size(st3,1)
          % if we exceed the original allocated memory, double the allocated sizes
           fW(:, 2*size(st3,1))    = 0;
           fWpc(:,:,2*size(st3,1)) = 0;
           st3(2*size(st3,1), 1)   = 0;
        end

        st3(irange,1) = double(st); % spike times
        st3(irange,2) = double(id0+1); % spike clusters (1-indexing)
        st3(irange,3) = double(x0); % template amplitudes
        st3(irange,4) = double(vexp); % residual variance of this spike
        st3(irange,5) = korder; % batch from which this spike was found

        fW(:, irange) = gather(featW); % template features for this batch
        fWpc(:, :, irange) = gather(featPC); % PC features

        ntot = ntot + numel(x0); % keeps track of total number of spikes so far
    end

    if ibatch==niter-nBatches
        % allocate variables when switching to extraction phase
        st3 = zeros(1e7, 5); % this holds spike times, clusters and other info per spike

        % these next three store the low-d template decompositions
        rez.WA = zeros(nt0, Nfilt, Nrank,nBatches,  'single');
        rez.UA = zeros(Nchan, Nfilt, Nrank,nBatches,  'single');
        rez.muA = zeros(Nfilt, nBatches,  'single');

        % these ones store features per spike
        fW  = zeros(Nnearest, 1e7, 'single'); % Nnearest is the number of nearest templates to store features for
        fWpc = zeros(NchanNear, Nrank, 1e7, 'single'); % NchanNear is the number of nearest channels to take PC features from
    end

    if (rem(ibatch, 100)==1)
        % this is some of the relevant diagnostic information to be printed during training
        fprintf('%2.2f sec, %d / %d batches, %d units, nspks: %2.4f, mu: %2.4f, nst0: %d, merges: %2.4f, %2.4f \n', ...
            toc, ibatch, niter, Nfilt, sum(nsp), median(mu), numel(st0), ndrop)

        % these diagnostic figures should be mostly self-explanatory
        if ibatch==1
            figHand = figure;
        else
            figure(figHand);
        end

       if ops.fig
           make_fig(W, U, mu, nsp)
        end
    end
end
fclose(fid);
toc

% discards the unused portion of the arrays
st3 = st3(1:ntot, :);
fW = fW(:, 1:ntot);
fWpc = fWpc(:,:, 1:ntot);

% just display the total number of spikes
ntot

rez.st3 = st3;
rez.st2 = st3; % keep also an st2 copy, because st3 will be over-written by one of the post-processing steps

% the similarity score between templates is simply the correlation,
% taken as the max over several consecutive time delays
rez.simScore = gather(max(WtW, [], 3));

% the template features are stored in cProj, like in Kilosort1
rez.cProj    = fW';
% the neihboring templates idnices are stored in iNeigh
rez.iNeigh   = gather(iList);

rez.ops = ops; % update these (only rez comes out of this script)
rez.nsp = nsp;

%  permute the PC projections in the right order
rez.cProjPC     = permute(fWpc, [3 2 1]); %zeros(size(st3,1), 3, nNeighPC, 'single');
% iNeighPC keeps the indices of the channels corresponding to the PC features
rez.iNeighPC    = gather(iC(:, iW));

% this whole next block is just done to compress the compressed templates
% we separately svd the time components of each template, and the spatial components
% this also requires a careful decompression function, available somewhere in the GUI code
nKeep = min(Nchan*3,20); % how many PCs to keep
rez.W_a = zeros(nt0 * Nrank, nKeep, Nfilt, 'single');
rez.W_b = zeros(nBatches, nKeep, Nfilt, 'single');
rez.U_a = zeros(Nchan* Nrank, nKeep, Nfilt, 'single');
rez.U_b = zeros(nBatches, nKeep, Nfilt, 'single');
for j = 1:Nfilt
    % do this for every template separately
    WA = reshape(rez.WA(:, j, :, :), [], nBatches);
    WA = gpuArray(WA); % svd on the GPU was faster for this, but the Python randomized CPU version might be faster still
    [A, B, C] = svdecon(WA);
    % W_a times W_b results in a reconstruction of the time components
    rez.W_a(:,:,j) = gather(A(:, 1:nKeep) * B(1:nKeep, 1:nKeep));
    rez.W_b(:,:,j) = gather(C(:, 1:nKeep));

    UA = reshape(rez.UA(:, j, :, :), [], nBatches);
    UA = gpuArray(UA);
    [A, B, C] = svdecon(UA);
    % U_a times U_b results in a reconstruction of the time components
    rez.U_a(:,:,j) = gather(A(:, 1:nKeep) * B(1:nKeep, 1:nKeep));
    rez.U_b(:,:,j) = gather(C(:, 1:nKeep));
end

fprintf('Finished compressing time-varying templates \n')
%%

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function U = spm_compact_svd(Y,xyz,nu)
% local SVD with compact support for large matrices
% FORMAT U = spm_compact_svd(Y,xyz,nu)
% Y     - matrix
% xyz   - location
% nu    - number of vectors
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_compact_svd.m 5219 2013-01-29 17:07:07Z spm $

% get orders
%--------------------------------------------------------------------------
ns     = size(Y,1);                % number of samples
nv     = size(Y,2);                % number of components (voxels/vertices)

% get kernel (compact vectors)
%--------------------------------------------------------------------------
nc    = max(fix(nv/nu),1);         % voxels in compact support
C     = sum(Y.^2);                 % variance of Y
U     = spalloc(nv,nu,nc*nu);
J     = 1:nv;
for i = 1:nu

    % find maximum variance voxel
    %----------------------------------------------------------------------
    [v,j] = max(C);
    d     = 0;
    for k = 1:size(xyz,1)
        d  = d + (xyz(k,:) - xyz(k,j)).^2;
    end
    [d,j] = sort(d);
    try
        j = j(1:nc);
    end

    % save principal eigenvector
    %----------------------------------------------------------------------
    k        = J(j);
    u        = spm_svd(Y(:,k)');
    U(k,i)   = u(:,1);

    % remove compact support voxels and start again
    %----------------------------------------------------------------------
    J(j)     = [];
    C(j)     = [];
    xyz(:,j) = [];

end

author: spm


text:
function svgClipping(s, path)
% Adds a clipping path
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgClipping(s, path)
% Parameters:
%   s : Array of plot object handles
%   path : Clipping path nx3 or nx2.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    userdata.svg.ClippingPath = path;
    set(s(i),'UserData', userdata);
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function name = p28_name ( )

%*****************************************************************************80
%
%% P28_NAME returns the name of problem 28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'Genz #2 / Patterson #6, Product Peak';

  return
end


text:
%                                                y = makereal(x,K,cpx)
% MAKEREAL  Converts matrix in MATLAB-complex format to internal
%  SeDuMi format.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi

function y = makereal(x,K,cpx) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
function p = er_ftest(dof1, dof2, F, dof2max)
%
% p = er_ftest(dof1, dof2, F, <dof2max>)
%
% Computes p-value given F-value. p and F can be vectors.
% dof1 = dof of numerator (Number of Rows in RM)
% dof2 = dof of denominator (DOF)
%
% Ref: Numerical Rec in C, pg 229.
%
% ras, 05/05; based on Ftest in fs-fast code.

if(nargin ~= 3 & nargin ~= 4)
    msg = 'Usage: p = FTest(dof1, dof2, F, <dof2max>)';
    error(msg);
    return;
end

if(length(dof1) > 1)
    error('dof1 must be a scalar');
end
if(length(dof2) > 1)
    error('dof2 must be a scalar');
end

if(nargin == 4) dof2 = min(dof2,dof2max); end

z = dof2./(dof2 + dof1 * F);

% 08/05/04: temp debug stuff, b/c of a funny session
if any(z(:) < 0 | z(:) > 1 | isnan(z(:))) | ~isreal(z)
    fprintf('any < 0?: %i \n',any(z(:) < 0))
    fprintf('any > 1?: %i \n',any(z(:) > 1))
    fprintf('Isnan?: %i \n',any(isnan(z(:))))
    fprintf('~Isreal?: %i \n',~isreal(z))
    z = 0.5 * ones(size(z)); % temp hack so that betainc doesn't fail
end
% fprintf('DEBUG: Z ranges from %f to %f \n',min(z(:)),max(z(:)));

p = betainc(z, dof2/2, dof1/2);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%  Matlab script to read DEM 3 arc seconds file.
%
%  Creates the following variables:
%    elev: elevation data (row = lat; col = lon)
%    xv: longitude vector
%    yv: latitude vector
%    clat & clon: SouthEast corner latitude and longitude of the file

%  Written 27 Nov 96  by Guy Tytgat

report_this_filefun(mfilename('fullpath'));

[file,path] = uigetfile('*','Select a .3CD file',300,300);

if file == []
    disp('Error selecting file')
    clear file path
    return
elseif file == 0
    disp('No file selected')
    clear file path
    return
end

filename = sprintf('%s%s',path,file);
[fid,message2] = fopen(filename);
if file(7:9) ~= '3cd'
    error('Wrong file type selected')
    return
end

nrow = 1201;
clat = str2double(file(1:2));
clon = str2double(file(3:5));

if fid ~= -1
    if clat < 50
        ncol = 1201;
    elseif (50 <= clat)  &&  (clat < 70)
        ncol = 601;
    elseif (70 <= clat)  &&  (clat < 75)
        ncol = 401;
    elseif (75 <= clat)  &&  (clat < 80)
        ncol = 301;
    elseif (80 <= clat)  &&  (clat <= 90)
        ncol = 201;
    else
        error('Incorrect file type selected')
        return
    end

    elev = ones([nrow,ncol]);
    [elev, count] = fread(fid,[nrow,ncol],'short');
    if (nrow*ncol) ~= count
        disp('WARNING: Not all data points were read')
    end

    xv = -(clon+1):1/ncol:-clon;
    yv = clat:1/nrow:clat+1;
    clear count fid file filename message2 ncol nrow path
else
    clear fid file filename message2 clat clon path nrow
end

author: CelsoReyes


text:
% NMF-MU: NMF solved by Multiplicative Updates
% process_video('NMF', 'NMF-MU', 'dataset/demo.avi', 'output/demo_NMF-MU.avi');
alg_path_aux = fullfile(lrs_conf.nmf_path,'NMF-DTU-Toolbox');
addpath(genpath(alg_path_aux));

M = sparse(M);

% mm: Multiplicative update method using euclidean distance measure.
[W, H] = nmf(M,1,'mm');

L = W * H;
S = M - L;

rmpath(genpath(alg_path_aux));

author: andrewssobral

repo: lrslibrary


text:
function pass = test_conj( pref )
% Test CONJ

if ( nargin == 0 )
    pref = chebfunpref;
end

tol = 100*pref.cheb2Prefs.chebfun2eps;

f = chebfun2(@(x,y) cos(x.*y));
g = conj( f );
pass(1) = ( norm( f - g ) < tol );

f = chebfun2(@(x,y) cos(x.*y));
g = conj( 1i*f );
pass(2) = ( norm( 1i*f + g ) < 100*tol );

f1 = chebfun2(@(x,y) cos(x.*y));
f2 = chebfun2(@(x,y) sin(x + y.^2));
g = conj( f1 + 1i*f2 );
pass(3) = ( norm( f1 - 1i*f2 - g ) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%%
% Test for recovering a 2-D metric field over the a*b* plane of the L*a*b*
% domain that fits the ciede2000 perceptual metric.

addpath('../toolbox/');
addpath('../colors_functions/');
addpath('../image_blur/');
addpath('../blur_functions//');
addpath('../convolutional_wasserstein/');
% addpath('../../data/images/colors/'); % low-res images
addpath('../../data/images/colors-big/'); % high-res images

rep = '../results/color-metric/';
if not(exist(rep))
    mkdir(rep);
end

%%
% L in [0,100], A in [-85,100], B in [-107,95]
% base fixed luminance for the metric
L = 50;

%%
% helpers

mmin = @(x)min(x(:));
mmax = @(x)max(x(:));
%
eta = .001;
DiffDir = @(ab,h)deltaE2000([L ab],[L ab]+eta*[0 h])/eta;

%%
% test

ab = [3 -80];
h = [1 1];
d = DiffDir(ab,h);

%%
% fit a single metric

M = fit_metric(@(h)DiffDir(ab,h));

%%
% Display a comparison between metric and fit.

P = 100;
tlist = linspace(0,2*pi,P+1)'; tlist(end)=[];
EuclMetric = @(M,h)sqrt( h*M*h' );

D = arrayfun(@(h1,h2)DiffDir(ab,[h1 h2]), cos(tlist), sin(tlist));
D1 = arrayfun(@(h1,h2)EuclMetric(M,[h1 h2]), cos(tlist), sin(tlist));

clf; hold on;
plot( cos(tlist).*D, sin(tlist).*D, 'b' );
plot( cos(tlist).*D1, sin(tlist).*D1, 'r:' );
axis equal;

%%
% compute some image histogram

name = 'blue-7';
n = 512;
f = rescale( load_image(name, n) );
A = colorspace(['RGB->LAB'], f);
% compute histo
N = 60;
H = compute_histogram_2d(A(:,:,2),A(:,:,3),N);
delta = .0001;
% display
clf;
imageplot(log(H+delta));

%%
% Compute a field of metrics
% A in [-85,100], B in [-107,95]

% background color image
v = 100;
a = linspace(-v,v,N);
b = linspace(-v,v,N);
[B,A] = meshgrid(b,a);
L = 80*ones(N,N);
CM = colorspace('LAB->RGB', cat(3,L,A,B));
% metric
F = fit_metric_field(DiffDir, arange,brange,N);
% inverse metric
[e1,e2,l1,l2] = perform_tensor_decomp(F);
G = perform_tensor_recomp(e1,e2,1./l1,1./l2);

options.sub = 2;
clf;
h = plot_tensor_field(F, CM, options);
saveas(gcf, [rep name '-cie-metric.png'], 'png');
clf;
h = plot_tensor_field(G, CM, options);
saveas(gcf, [rep name '-cie-inverse.png'], 'png');

author: gpeyre


text:
%% Demo: fMRI searchlights LDA classifier
%
% The data used here is available from http://cosmomvpa.org/datadb.zip
%
% This example uses the following dataset:
% + 'digit'
%    A participant made finger pressed with the index and middle finger of
%    the right hand during 4 runs in an fMRI study. Each run was divided in
%    4 blocks with presses of each finger and analyzed with the GLM,
%    resulting in 2*4*4=32 t-values
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set data paths
% The function cosmo_config() returns a struct containing paths to tutorial
% data. (Alternatively the paths can be set manually without using
% cosmo_config.)
config=cosmo_config();

digit_study_path=fullfile(config.tutorial_data_path,'digit');
readme_fn=fullfile(digit_study_path,'README');
cosmo_type(readme_fn);

output_path=config.output_data_path;

% reset citation list
cosmo_check_external('-tic');

%% LDA classifier searchlight analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This analysis identified brain regions where the categories can be
% distinguished using an odd-even partitioning scheme and a Linear
% Discriminant Analysis (LDA) classifier.

data_path=digit_study_path;
data_fn=fullfile(data_path,'glm_T_stats_perblock+orig.HEAD');
mask_fn=fullfile(data_path,'brain_mask+orig.HEAD');

targets=repmat(1:2,1,16)';    % class labels: 1 2 1 2 1 2 1 2 1 2 ... 1 2
chunks=floor(((1:32)-1)/8)+1; % run labels:   1 1 1 1 1 1 1 1 2 2 ... 4 4

ds_per_run = cosmo_fmri_dataset(data_fn, 'mask', mask_fn,...
                                'targets',targets,'chunks',chunks);

% print dataset
fprintf('Dataset input:\n');
cosmo_disp(ds_per_run);

% Use the cosmo_cross_validation_measure and set its parameters
% (classifier and partitions) in a measure_args struct.
measure = @cosmo_crossvalidation_measure;
measure_args = struct();

% Define which classifier to use, using a function handle.
% Alternatives are @cosmo_classify_{svm,matlabsvm,libsvm,nn,naive_bayes}
measure_args.classifier = @cosmo_classify_lda;

% Set partition scheme. odd_even is fast; for publication-quality analysis
% nfold_partitioner is recommended.
% Alternatives are:
% - cosmo_nfold_partitioner    (take-one-chunk-out crossvalidation)
% - cosmo_nchoosek_partitioner (take-K-chunks-out  "             ").
measure_args.partitions = cosmo_oddeven_partitioner(ds_per_run);

% print measure and arguments
fprintf('Searchlight measure:\n');
cosmo_disp(measure);
fprintf('Searchlight measure arguments:\n');
cosmo_disp(measure_args);

% Define a neighborhood with approximately 100 voxels in each searchlight.
nvoxels_per_searchlight=100;
nbrhood=cosmo_spherical_neighborhood(ds_per_run,...
                        'count',nvoxels_per_searchlight);

% Run the searchlight
lda_results = cosmo_searchlight(ds_per_run,nbrhood,measure,measure_args);

% print output dataset
fprintf('Dataset output:\n');
cosmo_disp(lda_results);

% Plot the output
cosmo_plot_slices(lda_results);

% Define output location
output_fn=fullfile(output_path,'lda_searchlight+orig');

% Store results to disc
cosmo_map2fmri(lda_results, output_fn);

% Show citation information
cosmo_check_external('-cite');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [seq, init_image] = get_sequence_info(seq)

if ~isfield(seq, 'format') || isempty(seq.format)
    if isempty(seq)
        seq.format = 'vot';
    else
        seq.format = 'otb';
    end
end

seq.frame = 0;

if strcmpi(seq.format, 'otb')
    seq.image_files = seq.s_frames;
    seq = rmfield(seq, 's_frames');
    seq.init_sz = [seq.init_rect(1,4), seq.init_rect(1,3)];
    seq.init_pos = [seq.init_rect(1,2), seq.init_rect(1,1)] + (seq.init_sz - 1)/2;
    seq.num_frames = numel(seq.image_files);
    seq.rect_position = zeros(seq.num_frames, 4);
    init_image = imread(seq.image_files{1});
elseif strcmpi(seq.format, 'vot')
    [seq.handle, init_image_file, init_region] = vot('polygon');

    if isempty(init_image_file)
        init_image = [];
        return;
    end

    init_image = imread(init_image_file);

    bb_scale = 1;

    % If the provided region is a polygon ...
    if numel(init_region) > 4
        % Init with an axis aligned bounding box with correct area and center
        % coordinate
        cx = mean(init_region(1:2:end));
        cy = mean(init_region(2:2:end));
        x1 = min(init_region(1:2:end));
        x2 = max(init_region(1:2:end));
        y1 = min(init_region(2:2:end));
        y2 = max(init_region(2:2:end));
        A1 = norm(init_region(1:2) - init_region(3:4)) * norm(init_region(3:4) - init_region(5:6));
        A2 = (x2 - x1) * (y2 - y1);
        s = sqrt(A1/A2);
        w = s * (x2 - x1) + 1;
        h = s * (y2 - y1) + 1;
    else
        cx = init_region(1) + (init_region(3) - 1)/2;
        cy = init_region(2) + (init_region(4) - 1)/2;
        w = init_region(3);
        h = init_region(4);
    end

    init_c = [cy cx];
    init_sz = bb_scale * [h w];

    im_size = size(init_image);

    seq.init_pos = init_c;
    seq.init_sz = min(max(round(init_sz), [1 1]), im_size(1:2));
    seq.num_frames = Inf;
    seq.region = init_region;
else
    error('Uknown sequence format');
end

author: martin-danelljan

repo: ECO