text:
function value = lens_half_w_area_2d ( r, w )

%*****************************************************************************80
%
%% LENS_HALF_W_AREA_2D returns the area of a circular half lens in 2D.
%
%  Discussion:
%
%    A circular half lens is formed by drawing a circular arc, and joining its endpoints.
%    This particular half lens is described by the "width" of the region.  In other words,
%    it is the portion of the circle under water if the width
%    of the water surface is W.  There are two possible values for this
%    area, A and (PI*R**2-A).  The routine returns the smaller of the
%    two values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real W, the width of the half lens.
%
%    Output, real VALUE, the area of the half lens.
%
  if ( w <= 0.0E+00 )

    value = 0.0E+00;

  elseif ( 2.0E+00 * r <= w )

    value = 0.5E+00 * pi * r * r;

  else

    half_width = 0.5E+00 * w;
    h = r - sqrt ( r * r - half_width * half_width );
    angle = 2.0E+00 * atan2 ( half_width, r - h );
    sector = r * r * angle / 2.0E+00;
    triangle = ( r - h ) * half_width;
    value = sector - triangle;

  end

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% Artificial example motivated by Gary Christensen
%==============================================================================

viewPara = { 'viewImage','viewImage2D','colormap',flipud(gray(256))};
imgPara  = {'imgModel','splineInter','regularizer','moments','theta',1};
traPara  = {'trafo','affine2D'};
disPara  = {'distance','SSD'};
regPara  = {'regularizer','mbHyperElastic','alpha',1,...
  'alphaLength',1000,'alphaArea',0,'alphaVolume',500};

expfile = jpgs2data('','disc.jpg','c.jpg', ...
  'omega',[0,1,0,1],'m',[128,128],...
  'viewPara',viewPara,'imgPara',imgPara);

save(expfile,'-append','viewPara','imgPara','traPara','disPara','regPara');
checkDataFile
%==============================================================================


text:
function value = hexagon_contains_point_2d ( v, p )

%*****************************************************************************80
%
%% HEXAGON_CONTAINS_POINT_2D finds if a point is inside a hexagon in 2D.
%
%  Discussion:
%
%    This test is only valid if the hexagon is convex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real V(2,6), the vertics, in counterclockwise order.
%
%    Input, real P(2), the point to be tested.
%
%    Output, logical VALUE, is TRUE if X is in the hexagon.
%
  n = 6;
%
%  A point is inside a convex hexagon if and only if it is "inside"
%  each of the 6 halfplanes defined by lines through consecutive
%  vertices.
%
  for i = 1 : n

    j = mod ( i, n ) + 1;

    if (  v(1,i) * ( v(2,j) - p(2  ) ) ...
        + v(1,j) * ( p(2  ) - v(2,i) ) ...
        + p(1  ) * ( v(2,i) - v(2,j) ) < 0.0 )

      value = 0;
      return

    end

  end

  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function st_io_test03 ( )

%*****************************************************************************80
%
%% ST_IO_TEST03 tests ST_SORT_A.
%
%  Discussion:
%
%    The matrix is:
%
%      11  12   0   0  15
%      21  22   0   0   0
%       0   0  33   0  35
%       0   0   0  44   0
%      51   0  53   0  55
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2014
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 5;
  nst = 11;
  ast = [ 51.0, 12.0, 11.0, 33.0, 15.0, 53.0, 55.0, 22.0, 35.0, 44.0, 21.0 ]';
  ist = [ 5, 1, 1, 3, 1, 5, 5, 2, 3, 4, 2 ];
  jst = [ 1, 2, 1, 3, 5, 3, 5, 2, 5, 4, 1 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'ST_IO_TEST03\n' );
  fprintf ( 1, '  ST_SORT_A sorts an ST matrix by columns.\n' );

  i_min = min ( ist );
  i_max = max ( ist );
  j_min = min ( jst );
  j_max = max ( jst );

  st_header_print ( i_min, i_max, j_min, j_max, m, n, nst );

  st_print ( m, n, nst, ist, jst, ast, '  Matrix data before sorting:' );

  [ ist, jst, ast ] = st_sort_a ( m, n, nst, ist, jst, ast );

  st_print ( m, n, nst, ist, jst, ast, '  Matrix data after sorting:' );

  return
end


text:
function test_tutorial_beamformingextended

% MEM 5gb
% WALLTIME 00:30:00
% DEPENDENCY ft_read_mri ft_redefinetrial ft_freqanalysis ft_volumesegment ft_appenddata ft_selectdata ft_prepare_singleshell ft_sourceanalysis ft_prepare_leadfield ft_prepare_headmodel ft_prepare_sourcemodel ft_plot_headmodel ft_plot_sens ft_plot_mesh ft_sourceinterpolate ft_sourceplot

datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/sensor_analysis');
mridir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer_extended');
templatedir  = dccnpath('/home/common/matlab/fieldtrip/template/sourcemodel');

load(fullfile(datadir, 'subjectK.mat'));

% Time windows of interest
data = ft_appenddata([], data_left, data_right);
cfg = [];
cfg.toilim = [-0.8 1.1];
cfg.minlength = 'maxperlen';
data = ft_redefinetrial(cfg, data);

cfg = [];
cfg.toilim = [-0.8 0];
data_bsl = ft_redefinetrial(cfg, data);

cfg.toilim = [0.3 1.1];
data_exp = ft_redefinetrial(cfg, data);

cfg = [];
data_cmb = ft_appenddata(cfg, data_bsl, data_exp);
% code the trial: 0 = baseline, 1 = experimental condition
data_cmb.trialinfo = [zeros(length(data_bsl.trial), 1); ones(length(data_exp.trial), 1)];

%% calculating the cross spectral density matrix
cfg = [];
cfg.method        = 'mtmfft';
cfg.output        = 'fourier'; % add hint: why fourier?
cfg.tapsmofrq     = 15;
cfg.foi           = 55;
cfg.keeptrials    = 'yes';
freq_cmb = ft_freqanalysis(cfg, data_cmb);

cfg = [];
cfg.trials = freq_cmb.trialinfo == 0;
freq_bsl = ft_selectdata(cfg, freq_cmb);
freq_bsl.cumtapcnt = freq_cmb.cumtapcnt(cfg.trials);
freq_bsl.cumsumcnt = freq_cmb.cumsumcnt(cfg.trials);
cfg.trials = freq_cmb.trialinfo == 1;
freq_exp = ft_selectdata(cfg, freq_cmb);
freq_exp.cumtapcnt = freq_cmb.cumtapcnt(cfg.trials);
freq_exp.cumsumcnt = freq_cmb.cumsumcnt(cfg.trials);

%%foward model and lead field

mri = ft_read_mri(fullfile(mridir, 'subjectK.mri'));
% cfg = [];
% [segmentedmri] = ft_volumesegment(cfg, mri);
%
oldsegmented = load(fullfile(mridir, 'segmentedmri.mat'));
segmentedmri = oldsegmented.segmentedmri;
%
% % check whether the segmentation gives results of more than 99% consistency
% assert(max(abs(oldsegmented.segmentedmri.gray(:)-segmentedmri.gray(:))) < .01, 'Gray matter segmentation differs from stored data')
% assert(max(abs(oldsegmented.segmentedmri.csf(:)-segmentedmri.csf(:))) < .01, 'CSF segmentation differs from stored data')
% assert(max(abs(oldsegmented.segmentedmri.white(:)-segmentedmri.white(:))) < .01, 'White matter segmentation differs from stored data')
% % transformation should be absolutely identical
% assert(isequal(oldsegmented.segmentedmri.transform, segmentedmri.transform), 'Transform differs from stored data')
%
% %save segmentedmri segmentedmri
%
% % add anatomical information to the segmentation
% segmentedmri.transform = mri.transform;
% segmentedmri.anatomy   = mri.anatomy;
% % call ft_sourceplot
% cfg = [];
% cfg.funparameter = 'gray';
% ft_sourceplot(cfg,segmentedmri);

% create ht head model from the segmented brain surface
cfg = [];
cfg.method = 'singleshell';
hdm = ft_prepare_headmodel(cfg, segmentedmri);

template = load(fullfile(templatedir, 'standard_sourcemodel3d8mm'));
% inverse-warp the subject specific grid to the template grid
cfg = [];
cfg.sourcemodel.warpmni   = 'yes';
cfg.sourcemodel.template  = template.sourcemodel;
cfg.sourcemodel.nonlinear = 'yes'; % use non-linear normalization
cfg.mri            = mri;
sourcemodel        = ft_prepare_sourcemodel(cfg);

figure;
hold on;
% note that when calling different plotting routines, all objects that we plot
% need to be in the same unit and coordinate space, here, we need to transform
% the head model to 'cm'
ft_plot_headmodel(ft_convert_units(hdm, freq_cmb.grad.unit), 'edgecolor', 'none');
alpha 0.4;
ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:));
ft_plot_sens(freq_cmb.grad);

cfg         = [];
cfg.sourcemodel    = sourcemodel;
cfg.headmodel = hdm;
cfg.channel = {'MEG'};
cfg.grad    = freq_cmb.grad;
sourcemodel_lf     = ft_prepare_leadfield(cfg, freq_cmb);

%% contrasting source activity
cfg                   = [];
cfg.frequency         = freq_cmb.freq;
cfg.grad              = freq_cmb.grad;
cfg.method            = 'dics';
cfg.keeptrials        = 'yes';
cfg.channel           = 'MEG';
cfg.sourcemodel       = sourcemodel_lf;
cfg.headmodel         = hdm;
cfg.keeptrials        = 'yes';
cfg.dics.lambda       = '5%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.fixedori     = 'yes';
cfg.dics.realfilter   = 'yes';
source  = ft_sourceanalysis(cfg, freq_cmb);

% beam pre- and poststim by using the common filter
cfg.sourcemodel.label    = source.avg.label;
cfg.sourcemodel.filter   = source.avg.filter;
source_bsl  = ft_sourceanalysis(cfg, freq_bsl);
source_exp  = ft_sourceanalysis(cfg, freq_exp);

source_diff = source_exp;
source_diff.avg.pow = (source_exp.avg.pow ./ source_bsl.avg.pow) - 1;
source_diff.pos = template.sourcemodel.pos;
source_diff.dim = template.sourcemodel.dim;

% note that the exact directory is user-specific
templatefile = dccnpath('/home/common/matlab/fieldtrip/external/spm8/templates/T1.nii');
template_mri = ft_read_mri(templatefile);
template_mri.coordsys = 'spm';

cfg              = [];
cfg.parameter    = 'pow';
cfg.interpmethod = 'nearest';
source_diff_int  = ft_sourceinterpolate(cfg, source_diff, template_mri);

cfg               = [];
cfg.method        = 'slice';
cfg.funparameter  = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim   = [0.0 1.2];
cfg.opacitylim    = [0.0 1.2];
cfg.opacitymap    = 'rampup';
ft_sourceplot(cfg,source_diff_int);

cfg.method = 'ortho';
cfg.atlas  = dccnpath('/home/common/matlab/fieldtrip/template/atlas/aal/ROI_MNI_V4.nii');
ft_sourceplot(cfg,source_diff_int);

cfg.method = 'surface';
cfg.projmethod     = 'nearest';
cfg.surffile       = 'surface_white_both.mat';
cfg.surfdownsample = 10;
ft_sourceplot(cfg,source_diff_int);

%% coherence beaming

data = ft_appenddata([], data_left, data_right);
cfg                 = [];
cfg.toilim          = [-1 -0.0025];
cfg.minlength       = 'maxperlen'; % this ensures all resulting trials are equal length
data_stim           = ft_redefinetrial(cfg, data);

cfg                 = [];
cfg.output          = 'powandcsd';
cfg.method          = 'mtmfft';
cfg.taper           = 'dpss';
cfg.tapsmofrq       = 5;
cfg.foi             = 20;
cfg.keeptrials      = 'yes';
cfg.channel         = {'MEG' 'EMGlft' 'EMGrgt'};
cfg.channelcmb      = {'MEG' 'MEG'; 'MEG' 'EMGlft'; 'MEG' 'EMGrgt'};
freq_csd            = ft_freqanalysis(cfg, data_stim);

cfg                 = [];
cfg.method          = 'dics';
cfg.refchan         = 'EMGlft';
cfg.frequency       = 20;
cfg.headmodel       = hdm;
cfg.sourcemodel     = sourcemodel;
source_coh_lft      = ft_sourceanalysis(cfg, freq_csd);

source_coh_lft.pos = template.sourcemodel.pos;
source_coh_lft.dim = template.sourcemodel.dim;

% note that the exact directory is user-specific
templatefile = dccnpath('/home/common/matlab/fieldtrip/external/spm8/templates/T1.nii');
template_mri = ft_read_mri(templatefile);
template_mri.coordsys = 'spm';

cfg              = [];
cfg.parameter    = 'coh';
cfg.interpmethod = 'nearest';
source_coh_int   = ft_sourceinterpolate(cfg, source_coh_lft, template_mri);

cfg              = [];
cfg.method       = 'ortho';
cfg.funparameter = 'coh';
cfg.funcolormap = 'jet';

cfg.funcolorlim   = [00 .15];
cfg.opacitylim    = [00 .15];

cfg.maskparameter = cfg.funparameter;
cfg.opacitymap    = 'rampup';

cfg.atlas         = dccnpath('/home/common/matlab/fieldtrip/template/atlas/aal/ROI_MNI_V4.nii');

ft_sourceplot(cfg, source_coh_int);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Computes the Euclidean distance between two n-dimensional point-sets
%
% Input:
%   g1,g2: input point sets
%
% Output:
%   Dist: output distance matrix
%

function Dist = dist_matrix(g1,g2)

l1 = size(g1,1);
l2 = size(g2,1);

Dist = zeros(l1,l2);
for it_dim = 1:size(g1,2)
    Dist = Dist + ...
        ((g1(:,it_dim)*ones(1,l2) - ones(l1,1)*g2(:,it_dim)').^2);
end
Dist = sqrt(Dist);

author: Sable

repo: mcbench-benchmarks


text:
function soln = smoothJointTrajectory(problem)

% Joint limits
qLow = problem.qLow;
qUpp = problem.qUpp;
dqMax = problem.dqMax;   %Joint speed limit

% Waypoints
tNode = problem.tNode;   %time
qNode = problem.qNode;   %angle

%Order of interpolating polynomial in each segment;
nGrid = problem.nGrid;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                      Build problem matricies                            %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

nDecVar = sum(nGrid);  %Number of decision variables
nSegment = length(nGrid);

% Rate constraint and objective matrix along each segment:
for i=1:length(nGrid)
    S(i) = chebyshevSegment(nGrid(i),tNode([i,i+1]),dqMax); %#ok<SAGROW>
end
A = structBlkDiag(S,'A');
H = structBlkDiag(S,'H');
D = structBlkDiag(S,'D');
DD = structBlkDiag(S,'DD');

b = []; t = [];  % Slow memory allocation - ok since few iterations
for i=1:length(nGrid)
    b = [b; S(i).b];  %#ok<AGROW>
    t = [t; S(i).t];  %#ok<AGROW>
end

%%%% Boundary Values:
nCstBc = 2*nSegment + 2*(nSegment+1);  %angle + rate

Aeq = zeros(nCstBc,nDecVar); beq = zeros(nCstBc,1);
finalIdx = cumsum(nGrid);
startIdx = 1 + [0, finalIdx(1:(end-1))];

cstIdx = 0;
for i=1:nSegment   %Angle at start of segment
    cstIdx = cstIdx + 1;
    Aeq(cstIdx,startIdx(i)) = 1;  beq(cstIdx) = qNode(i);
end
for i=1:nSegment   %Angle at end of segment
    cstIdx = cstIdx + 1;
    Aeq(cstIdx,finalIdx(i)) = 1;  beq(cstIdx) = qNode(i+1);
end

cstIdx = cstIdx + 1; Aeq(cstIdx,:) = D(1,:);  %zero initial velocity
cstIdx = cstIdx + 1; Aeq(cstIdx,:) = D(end,:);  %zero final velocity

% Defect constraint on rate at segment boundaries
for i=1:(nSegment-1)
    cstIdx = cstIdx + 1;
    Aeq(cstIdx,:) = D(startIdx(i+1),:) - D(finalIdx(i),:);
end

for i=1:(nSegment-1)
    cstIdx = cstIdx + 1;
    Aeq(cstIdx,:) = DD(startIdx(i+1),:) - DD(finalIdx(i),:);
end

%%%% Options:
options = optimset(...
    'Display','iter',...  % {'iter','final'}
    'Algorithm','interior-point-convex');

%%%% Build Problem:
problem.H = (H+H')/2;   %Correct for numerically introduced asymmetry
problem.f = zeros(nDecVar,1);
problem.Aineq = A;
problem.bineq = b;
problem.Aeq = Aeq;
problem.beq = beq;
problem.lb = qLow*ones(nDecVar,1);
problem.ub = qUpp*ones(nDecVar,1);
problem.x0 = [];
problem.options = options;
problem.solver = 'quadprog';

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             Solve Problem                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

[x, fVal, exitFlag, output] = quadprog(problem);
for i=1:nSegment
    SS.grid.t = S(i).t';
    tSpan = SS.grid.t([1,end]);
    SS.grid.q = x(startIdx(i):finalIdx(i))';
    SS.grid.dq = chebyshevDerivative(SS.grid.q,tSpan);
    SS.interp.t = linspace( tSpan(1), tSpan(2), 5*length(SS.grid.q) );
    [q,dq,ddq] = chebyshevInterpolate(SS.grid.q,SS.interp.t,tSpan);
    SS.interp.q = q;
    SS.interp.dq = dq;
    SS.interp.ddq = ddq;
    soln.segment(i) = SS;
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                             Return Solution                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

soln.info = output;
soln.info.exitFlag = exitFlag;
soln.info.fVal = fVal;

soln.grid.t = t';
soln.grid.q = x';

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%QMINCON Resolve redundancy in robots by avoiding joint limits
%
% A popular way to resolve redundant robots is to keep joints away from
% their mechanical limits to allow freer motion. This function will do
% that process. Requires fmincon from the optimization toolbox.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% This file requires file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% Syntax:
%  (1) [qstar, error, exitflag] = robot.qmincon(q)
%
% Outputs:
%  qstar    : Optimised joint angles, mxrobot.n where m = size(q,2)
%  error    : The error measurement (value of objective function)
%  exitflag : The exitflag direct from fmincon
%
% Inputs:
%  q : Joint configuration(s), may be an mxrobot.n matrix of m poses
%
% See also fmincon ikcon ikunc

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license

function [qstar, error, exitflag] = qmincon(robot, q)

M = size(q,1);
n = robot.n;

qstar = zeros(M,n);
error = zeros(M,1);
exitflag = zeros(M,1);

opt = optimoptions('fmincon', ...
    'Algorithm', 'active-set', ...
    'Display', 'off');

lb = robot.qlim(:,1);
ub = robot.qlim(:,2);

x_m = 0; % Little trick for setting x0 in first iteration of loop

for m = 1:M
    q_m = q(m,:);

    J = robot.jacobn(q(m,:));
    N = null(J);

    if isempty(N)
        error(pHRIWARE('error', 'Robot is not redundant'));
    end

    f = @(x) sumsqr((2*(N*x + q_m') - ub - lb)./(ub-lb));

    x0 = zeros(size(N,2), 1) + x_m;

    A = [N; -N];
    b = [ub-q_m'; q_m'-lb];

    [x_m, err_m, ef_m] = fmincon(f,x0,A,b,[],[],[],[],[],opt);

    qstar(m,:) = q(m,:) + (N*x_m)';
    error(m) = err_m;
    exitflag(m) = ef_m;
end

end

function s = sumsqr(A)
    s = sum(A.^2);
end

author: Allopart

repo: rbpf-gmapping


text:
function checkSpatialFT(MRSIStruct)
    if(~getFlags(MRSIStruct, 'spatialFT'))
        error('Please fourier transform along the spatial dimension')
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function slider = mrvSlider(pos, name, varargin);
%
% slider = mrvSlider([pos], [name], [properties]);
%
% Create a slider with an attached edit field,  both of which
% specify the same value,  and a text label describing the
% parameter being set. Often used in mrVista UIs.
%
% pos: position [lrhc_x lrhc_y x_size y_size] of box containing
%      slider,  edit field,  and labels.
% name: name of slider (will be default text in the label).
%
% properties: specify in pairs,  the name of the property and
% its value. Here are the properties an mrvSlider can have:
%    'parent': parent figure or uipanel for slider. Default gcf.
%    'range': [min max] values for the slider. Default [0 1].
%    'val': initial value for the slider. Default is [min] value.
%
%    'callback': callback function for the object. When the edit field
%          or sliders are modified by the user,  the control will
%          first set the other control to agree with the value that
%          was set. (E.g., if you type a value, it will first set the
%          slider to match that value.)
%		   Then it will evaluate the cb string entered here.
%          Default cb is empty.
%          For the callback, the current value of the slider will be
%          evaluated to the temporary variable 'val', so no get functions
%          need to be added. For example:
%          'callback', 'colormap( gray( val ))' will cause the slider
%          to set the figure's colormap to the [val] number of grays after
%          a value is selected.
%
%    'units': units for slider objects. Default 'normalized',  relative
%          to parent object.
%    'intFlag': flag for whether slider value should be constrained
%          to be an integer or not. Default 0.
%    'flexFlag': flag for whether slider should be flexible for
%          values entered outside the specified range. If 0,  when
%          a value is entered into the edit field outside the min
%          and max values set by the slider,  the value is automatically
%          reset to be within the bounds. If 1,  will allow the higher
%          value,  and expand the slider's range accordingly. Default
%          is 0,  reset out-of-bound values.
%    'sigFigs': # of significant figures to show in the edit field.
%           Default is 2. If intFlag is 1,  this is overridden and 0
%           sig figs are always shown.
%    'maxLabelFlag': flag to put a small label to the side of the
%           slider indicating the maximum value. Default is 0,
%           don't add this.
%    'color': background color of the object. Default is the
%           same value as the parent object.
%    'fontSize': size of fonts for the edit / text controls.
%           Default is 12 point.
%
% Returns a slider struct with handles to the slider,  edit,  and
% label(s),  the callback strings for the slider and edit controls,
% and the values of the properties above.
%
% Properties of the mrvSlider can be set later on using
% mrvSliderSet.m.
%
% ras,  07/05/05.
% ras,  09/06 -- imported into VISTASOFT repository from mrVista 2 tools.
if nargin<2,  help(mfilename); error('Not enough args.');            end

%%%%% Default properties
parent = gcf;
range = [0 1];
cb = '';
units = 'normalized';
intFlag = 0;
flexFlag = 0;
sigFigs = 2;
maxLabelFlag = 0;
fontName = 'Helvetica';
fontSize = 12;
sliderStep = [.03 .15];
visible = 'on';

%%%%% parse the property flags
for i = 1:2:length(varargin)
    prop = lower(varargin{i});
    switch prop
        case 'parent',  parent = varargin{i+1};
        case 'range',  range = varargin{i+1};
        case {'val', 'value'},  val = varargin{i+1};
        case {'cb', 'callback'},  cb = varargin{i+1};
        case 'units',  units = varargin{i+1};
        case 'intflag',  intFlag = varargin{i+1};
        case 'flexflag',  flexFlag = varargin{i+1};
        case 'sigfigs',  sigFigs = varargin{i+1};
        case 'maxlabelflag',  maxLabelFlag = varargin{i+1};
        case 'color',  color = varargin{i+1};
		case 'fontname',  fontName = varargin{i+1};
		case 'fontsize',  fontSize = varargin{i+1};
		case 'sliderstep', sliderStep = varargin{i+1};
        case 'visible',  visible = varargin{i+1};
        otherwise,  fprintf('mrvSlider: %s unrecognized property.\n', prop);
    end
end

% default initial val
if ~exist('val', 'var') | isempty(val),  val = range(1);   end

% default initial val
if ~exist('color', 'var')|isempty(color),
    if isequal(get(parent, 'Type'), 'figure')
        color = get(parent, 'Color');
    else
        color = get(parent, 'BackgroundColor');
    end
end

% range check
if range(2)<=range(1),     range(2) = range(1)+eps;         end

if flexFlag % overrides input ranges if flexFlag
    range(1) = min([range(1), val]);
    range(2) = max([range(2), val]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make slider
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slider callback string :
sliderCb = 'val = get(gcbo, ''Value'');';
if intFlag==1
	sliderCb = [sliderCb ' val = round(val); '];
end
sliderCb = sprintf('%s \n slider = get(gcbo, ''UserData'');', sliderCb);
sliderCb = sprintf('%s \n mrvSliderSet(slider, ''Val'', val);', sliderCb);
if ~isempty(cb)
    sliderCb = sprintf('%s \n %s;', sliderCb, cb);
end
sliderCb = sprintf('%s \n clear val slider; ', sliderCb);

% slider position within box
sliderPos = pos;
sliderPos(2) = pos(2) + 0.5*pos(4);
sliderPos(4) = 0.5*pos(4);

% create slider control
sliderHandle = uicontrol('Style', 'slider', ...
                        'Parent', parent, ...
                        'Units', units, ...
                        'Position', sliderPos, ...
                        'fontName', fontName, ...
						'FontSize', fontSize, ...
                        'ForegroundColor', color, ...
                        'BackgroundColor', [1 1 1], ...
                        'Min', range(1), ...
                        'Max', range(2), ...
                        'Val', val, ...
                        'SliderStep', sliderStep, ...
						'Visible', visible, ...
                        'Callback', sliderCb);
if intFlag==1
    %% make slider step across integer values
	sigFigs = 0;
	sliderStep = [.03 .15];

    if diff(range)>=1
		denom = max(1, abs(diff(range)));
        sliderStep(1) = 1/ denom; % should always page 1
	end

    set(sliderHandle, 'SliderStep', sliderStep);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if selected,  add a little label for the max slider val
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if maxLabelFlag==1
    sliderPos(3) = sliderPos(3) * 0.8;
    set(sliderHandle, 'Position', sliderPos);
    maxLabelPos = sliderPos;
    maxLabelPos(1) = sliderPos(1) + sliderPos(3);
    maxLabelPos(3) = sliderPos(3) * 0.2/0.8;
	maxTxt = sprintf( sprintf('%%.%if', sigFigs), range(2) );
    maxLabelHandle = uicontrol('Style', 'text', ...
                        'Parent', parent, ...
                        'fontName', fontName, ...
                        'FontSize', fontSize-2, ...
                        'Units', units, ...
                        'Position', maxLabelPos, ...
                        'BackgroundColor', 'w', ...
                        'Visible', visible, ...
                        'HorizontalAlignment', 'left', ...
                        'String', maxTxt);
else
    maxLabelHandle = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create text label
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% label pos within box
labelPos = pos;
labelPos(3) = 0.7*pos(3);
labelPos(4) = 0.5*pos(4);

% Create a text label w/ the slider name
labelHandle = uicontrol('Style', 'text', ...
                        'Parent', parent, ...
                        'fontName', fontName, ...
                        'FontSize', fontSize, ...
                        'Units', units, ...
                        'Position', labelPos, ...
                        'BackgroundColor', color, ...
                        'Visible', visible, ...
                        'HorizontalAlignment', 'left', ...
                        'String', [name ':']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create edit field
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% edit pos within box
editPos = pos;
editPos(1) = pos(1) + 0.5*pos(3);
editPos(3) = 0.5*pos(3);
editPos(4) = 0.5*pos(4);

% Edit field callback string:
editCb = 'val = str2num(get(gcbo, ''String''));';
if intFlag==1
	editCb = [editCb ' val = round(val); '];
end
editCb = sprintf('%s \n slider = get(gcbo, ''UserData'');', editCb);
editCb = sprintf('%s \n mrvSliderSet(slider, ''Val'', val);', editCb);
if ~isempty(cb)
    editCb = sprintf('%s \n %s;', editCb, cb);
end
editCb = sprintf('%s \n clear val slider;', editCb);

% make edit control
editHandle = uicontrol('Style', 'edit', ...
                        'Parent', parent, ...
                        'Units', units, ...
                        'fontName', fontName, ...
                        'FontSize', fontSize, ...
                        'Position', editPos, ...
                        'String', num2str(val), ...
                        'Visible', visible, ...
                        'BackgroundColor', color, ...
                        'Callback', editCb);

% add handles to struct
slider.name = name;
slider.sliderHandle = sliderHandle;
slider.editHandle = editHandle;
slider.labelHandle = labelHandle;
slider.maxLabelHandle = maxLabelHandle;
slider.range = range;
slider.intFlag = intFlag;
slider.flexFlag = flexFlag;
slider.sigFigs = sigFigs;
slider.maxLabelFlag = maxLabelFlag;
slider.sliderCb = sliderCb;
slider.editCb = editCb;

% set as user data to controls
set(editHandle, 'UserData', slider);
set(sliderHandle, 'UserData', slider);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% computes the pose vector v from a homogeneous transform A
function v=t2v(A)
  v = [A(1:2,3); atan2(A(2,1),A(1,1))];
end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
%% This function will generate the simulation data of an ODE function.
%You could determine the noise level by input variable "Noise". If you do
%not want any noise, set noise to zero. Please indicate whether your ODE
%function have control input, if the answer is yes, please set the
%"Control" as 1.

% Last Update: 2019/04/21
% Coded By: K

function [d_Data,Data]=Get_Sim_Data(ODE,state0,u,tspan,Noise,Control,Shuffle)
%% Get the size of the state and control
[N1,M1]=size(state0);
[N2,M2]=size(u);

%% Get simulation data by simulating the system using ODE113

% Determine the left hand side derivative
if Control==1
    y_list(1,:)=state0;
    d_y_list(1,:)=ODE(0,y_list(1,:),u(1,:));
    for i=2:length(u)
        [t_1,y_1] = ode15s(@(t_1,y_1)ODE(t_1,y_1,u(i-1,:)),tspan(1,i-1:i),state0);
        y_list(i,:)=y_1(end,:);
        d_y_list(i,:)=ODE(0,y_list(i,:),u(i,:));
        state0=y_list(i,:)';
    end
else
    %Simulate the system ODE
    [t,y]=ode15s(@(t,y)ODE(t,y),tspan,state0);
    y_list=y;
    % Get the derivative data
    d_y_list=ODE(0,y_list')';
end

%% Add some noise to the system
for i=1:N1
    Data(:,i)=y_list(:,i)+Noise*randn(size(y_list(:,i)));
end
%
for i=1:M2
    d_Data(:,i)=d_y_list(:,i)+Noise*randn(size(d_y_list(:,i)));
end

%% Shuffle the data
if Shuffle==1
    Sequence=randperm(size(Data,1));
    Data=Data(Sequence,:);
    d_Data=d_Data(Sequence,:);
end

author: dynamicslab

repo: SINDy-PI


text:
function [x0, exitflag] = vertex_untangle(tri, x, idx0)
% VERTEX_UNTANGLE  Find 2D coordinates for a free vertex surrounded
% by a polygon of fixed counter-clockwise sorted vertices such that the
% free vertex is untangled.
%
% This function is an enhanced implementation of the 2D case (i.e. each
% simplex is a triangle) of the "Optimization-based mesh untangling"
% method proposed by Freitag and Plassmann (2000).
%
% Note that the original paper contains some errors in the formulation of
% the linear programming problem. The correct problem formulation is
%
%     max b^T \pi (equivalenty: min -b^T \pi)
%     such that -A^T \pi <= c
%
% where b = [0 0 1], A^T = [ax ay -2*ones(N, 1)], N=number of
% triangles, and ax, ay, c as given in the paper.
%
% We check whether the free vertex is already untangled. In that case, we
% do not relocate it (this is an untangling algorithm, not a mesh
% improvement or mesh smoothing algorithm).
%
% X0 = vertex_untangle(TRI, X);
%
%   Let TRI, X describe a closed fan triangular mesh. That is, we have a
%   central or free vertex connected to N neighbours. The neighbours are
%   connected between them to form a closed perimeter. The perimeter is
%   assumed to be oriented in counter-clockwise orientation (that is, if a
%   solution exists, the areas of the triangles will all be positive).
%
%   We want to find the coordinates X0 for the free vertex, such that the
%   edges that connect vertices don't overlap. The N neighbours remain
%   fixed. If the free vertex was already producing any overlaps (i.e. the
%   mesh was entangled), this can be seen as an untangling algorithm.
%
%   TRI is a 3-column matrix. Each row contains the 3 nodes that form one
%   triangular facet in the mesh. The orientation provided by the triangles
%   will be assumed by the algorithm to the counter-clockwise, i.e. the
%   algorithm will produce a solution where all triangles have positive
%   signed areas.
%
%   X is a 2-column matrix. X(i, :) contains the xy-coordinates of the
%   i-th node in the mesh. Note that Freitag and Plassmann (2000) also
%   provided the solution for the 3D case (i.e. each simplex is a
%   tetrahedron), but we haven't implemented it here yet.
%
%   X0 is a 2-vector with the optimal coordinates of the central (a.k.a.
%   free) vertex, in a linear programming sense,
%
% [X0, EXITFLAG] = vertex_untangle(...)
%
%   EXITFLAG is the exit condition of the algorithm.
%
%     2: Free vertex was not tangled, no relocation performed.
%     1: Linear programming algorithm converged to a relocation solution.
%     0: Linear programming algorithm stopped because of maximum number of
%        iterations reached.
%    -1: Linear programming algorithm found no valid solution exists.
%
% L. A. Freitag and P. Plassmann, "Local optimization-based simplicial mesh
% untangling and improvement", International Journal for Numerical Methods
% in Engineering, 49(1):109-125, 2000.
%
% See also: sphtri_untangle, linprog, surfreorient.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.2.3
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(3, 3);
nargoutchk(0, 3);

if (size(tri, 2) ~= 3)
    error('TRI must have 3 columns')
end
if (size(x, 2) ~= 2)
    error('X must have 2 columns')
end
if (~isscalar(idx0))
    error('IDX must be a scalar')
end

% number of vertices in the neighbourhood
N = size(x, 1);

% adjacency matrix for the mesh. Edges are directed, thus the matrix is not
% symmetric. We do it this way to preserve the counter-clockwise
% orientation of the neighbourhood, i.e. the positive sign of the areas
d = sparse(N, N);
d(sub2ind([N, N], tri(:, 1), tri(:, 2))) = 1;
d(sub2ind([N, N], tri(:, 2), tri(:, 3))) = 1;
d(sub2ind([N, N], tri(:, 3), tri(:, 1))) = 1;

% number of neighbours for each vertex
nn = full(sum((d | d') ~= 0));

% valid neighbourhood check: the central vertex is connected to every other
% vertex
if (nn(idx0) ~= size(x, 1)-1)
    error(['Invalid neighbourhood: central vertex connected to ' ...
        num2str(nn(idx0)) ' vertices instead of ' num2str(size(x, 1) - 1)])
end
% valid neighbourhood check: non-central vertices must be connected to the
% central vertex, and to another two vertices each
if (any(nn([1:idx0-1 idx0+1:end]) ~= 3))
    error('Invalid neighbourhood: at least one boundary vertex not connected to another 3 vertices')
end

% get sorted list of vertices that form the perimeter
v = sort_perim_vertices(d, idx0);

% % DEBUG: plot boundary edges, marking the first segment
% subplot(2, 2, 1)
% hold off
% gplot(d, x)
% hold on
% plot(x(v([1:end 1]), 1), x(v([1:end 1]), 2), 'r')
% plot(x(v(1), 1), x(v(1), 2), 'ro')
% plot(x(v(2), 1), x(v(2), 2), 'rp')

% compute the optimal coordinates
% [x1, triarea, exitflag, output] ...
[x1, ~, exitflag] ...
    = linear_programming_coordinates(x, idx0, v);

% % DEBUG: plot solution of current orientation
% subplot(2, 2, 2)
% aux = x;
% aux(idx0, :) = x1;
% hold off
% gplot(d, aux)
% hold on
% plot(aux(v([1:end 1]), 1), aux(v([1:end 1]), 2), 'r')
% plot(aux(v(1), 1), aux(v(1), 2), 'ro')
% plot(aux(v(2), 1), aux(v(2), 2), 'r*')

% % DEBUG: plot solution of opposite orientation
% subplot(2, 2, 4)
% aux = x;
% aux(idx0, :) = x2;
% hold off
% gplot(d, aux)
% hold on
% plot(aux(v([1:end 1]), 1), aux(v([1:end 1]), 2), 'r')
% plot(aux(v(1), 1), aux(v(1), 2), 'ro')
% plot(aux(v(2), 1), aux(v(2), 2), 'r*')

%% choose the best solution

% if unsuccessful, we return NaN as a solution
if (exitflag < 0)

    x0 = nan(1, 2);
    exitflag = -1;

% else, we return the new location
else

    x0 = x1;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Auxiliary functions

% d: the sparse adjacency matrix with the directed connections between the
% perimeter vertices
%
% idx0: index of the central vertex
%
% v: ordered list of perimeter vertices in counter-clockwise orientation
% (so that the signed area is positive)
function v = sort_perim_vertices(d, idx0)

% clear up the connections from the perimeter to the central vertex
d(idx0, :) = 0;
d(:, idx0) = 0;

% Dijkstra's shortest paths from one arbitrary vertex in the perimeter to
% all the others
if (idx0 == 1)
    va = 2;
else
    va = 1;
end
[l, p] = dijkstra(d, va);

% furtherst vertex from origin vertex
l(isinf(l)) = 0;
[~, vb] = max(l);

% path from va to vb
v = graphpred2path(p, vb);

% invert path
v = v(end:-1:1);

end

% set up the Freitag and Plassmann (2000) linear programming problem and
% compute the solution. Note that the formulation in the original paper is
% wrong, as has been corrected here.
%
% x: 2D vertex coordinates
%
% v0: index of the free vertex
%
% v: indices of the perimeter vertices, sorted in counter-clockwise
% orientation (so that the signed area is positive)
%
% lb, ub: lower and upper bounds for the linear programming algorithm
function [x0, triarea, exitflag, output] ...
    = linear_programming_coordinates(x, v0, v)

% for convenience, express the boundary as a list of edges
e = [v(1:end)' v([2:end 1])'];

% number of triangles
Ntri = size(e, 1);

% terms for the linear programming problem
xi = x(e(:, 1), 1);
xj = x(e(:, 2), 1);
yi = x(e(:, 1), 2);
yj = x(e(:, 2), 2);
ax = yi - yj;
ay = xj - xi;
c = xi .* yj - xj .* yi;
A = [ax'; ay'; -2*ones(1, Ntri)];

% compute triangle areas (this formula comes from expanding the determinant
% form of a triangle's area)
triarea = 0.5 * (ax * x(v0, 1) + ay * x(v0, 2) + c);

% if all the areas are positive, then there's no tangling, and we won't
% relocate the free vertex (this is not a smoothing algorithm, thus the
% free vertex will only be moved if it's tangled)
if all(triarea > 0)
    x0 = x(v0, :);
    exitflag = 2;
    output = [];
    return
end

% the solution cannot be outside the convex hull of the boundary vertices.
% Thus, it cannot be outside of a rectangular box containing the boundary
% vertices. We can use this as the lower and upper bounds for the linear
% programming solution
lb = min(x(v, :));
ub = max(x(v, :));

% the third component of the linear programming solution vector = minimum
% area. We bound the area value between 0 and infinite. We could use as the
% upper bound the convex hull's area, but it's not worth wasting time
% computing it.
lb = [lb 0]';
ub = [ub Inf]';

% otherwise, we solve linear programming problem to try to untangle the
% free vertex
options = optimset('Display', 'off', 'Simplex', 'on', 'LargeScale', 'off');
[x0, ~, exitflag, output] = linprog([0 0 -1]', -A', c, [], [], ...
    lb, ub, [], options);

% the two first elements are the target coordinates of the central vertex
x0 = x0(1:2)';

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [psi] = dynpcm22psi(dynpcm2)
% Convert pressure from dynes per square centimeter to pounds per square
% inch
% Chad Greene 2012
psi = dynpcm2*0.0000145038;


text:
function [ hdr ] = Load_BV_hdr( file )
%LOAD_BP_HDR Summary of this function goes here
%   Detailed explanation goes here
global BMI
[fid, message]=fopen(fullfile([file '.vhdr']));
if fid == -1
    [fid, message] = fopen(fullfile(pathname, lower(filename)));
    if fid == -1
        error(message)
    end;
end

hdr.DataFile= getEntry(fid, 'DataFile=', 0, [file '.eeg']);
getEntry(fid, '[Common Infos]');
hdr.DataFile= getEntry(fid, 'DataFile=', 0, [file '.eeg']);
hdr.MarkerFile= getEntry(fid, 'MarkerFile=', 0, [file '.vmrk']);
hdr.DataFormat= getEntry(fid, 'DataFormat=', 0);
hdr.DataOrientation= getEntry(fid, 'DataOrientation=', 0);
hdr.DataType= getEntry(fid, 'DataType=', 0);
getEntry(fid, '[Common Infos]');
hdr.NumberOfChannels= str2num(getEntry(fid, 'NumberOfChannels='));
hdr.DataPoints= getEntry(fid, 'DataPoints=', 0, '0');
getEntry(fid, '[Common Infos]');
hdr.SamplingInterval= str2num(getEntry(fid, 'SamplingInterval='));

getEntry(fid, '[Binary Infos]');
hdr.BinaryFormat= getEntry(fid, 'BinaryFormat=', 0);
hdr.UseBigEndianOrder= getEntry(fid, 'UseBigEndianOrder=', 0);

getEntry(fid, '[Channel Infos]');
hdr.chan=cell(1,hdr.NumberOfChannels);
for chan=1:hdr.NumberOfChannels
    temp=getEntry(fid,['Ch' num2str(chan) '='],0);
    chInfo=strsplit(temp,',');
    chInfo_=chInfo(~cellfun('isempty',chInfo));
    hdr.ChannelResolution(chan)=str2num(chInfo_{2});
    hdr.chan{chan}=chInfo_{1};
end

hdr.fs= 1000000/hdr.SamplingInterval;
%make a check if hdr.fs is a whole number
if(hdr.fs ~= round(hdr.fs))
  warning('file_readBVheader: hdr.fs was not a whole number: %f',hdr.fs);
  hdr.fs= round(hdr.fs);
  hdr.SamplingInterval=1000000/hdr.fs;
end
fclose(fid);
end

function [entry, str]= getEntry(fid, keyword, mandatory, default_value)

if ~exist('mandatory','var'), mandatory=1; end
if ~exist('default_value','var'), default_value=[]; end
entry= 1;

if keyword(1)=='[',
    fseek(fid, 0, 'bof');
end
ok= 0;
while ~ok && ~feof(fid),
    str= fgets(fid);
    ok= strncmp(keyword, str, length(keyword));
end
if ~ok,
    if mandatory,
        error(sprintf('keyword <%s> not found', keyword));
    else
        entry= default_value;
        return;
    end
end
if keyword(end)=='=',
    entry= deblank(str(length(keyword)+1:end));
end

end

author: PatternRecognition

repo: OpenBMI


text:
function kern = linKernParamInit(kern)

% LINKERNPARAMINIT LIN kernel parameter initialisation.
% The linear kernel (LIN) is the simple inner product
% kernel. Sampling from this kernel produces linear functions.
%
% k(x_i, x_j) = sigma2 * x_i'*x_j
%
% There is one parameter, sigma2, which is stored in the field
% kern.variance.
%
% SEEALSO : linardKernParamInit
%
% FORMAT
% DESC initialises the linear
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

kern.variance = 1;
kern.nParams = 1;

kern.transforms.index = 1;
kern.transforms.type = optimiDefaultConstraint('positive');

kern.isStationary = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function blend_test04 ( )

%*****************************************************************************80
%
%% BLEND_TEST04 checks out BLEND_IJ_0D1 and BLEND_IJ_1D1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
  m1 = 5;
  m2 = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLEND_TEST04\n' );
  fprintf ( 1, '  BLEND_IJ_0D1 interpolates data in a table,\n' );
  fprintf ( 1, '  from corner data.\n' );
  fprintf ( 1, '  BLEND_IJ_1D1 interpolates data in a table,\n' );
  fprintf ( 1, '  from edge data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The table is %d rows by %d columns.\n', m1, m2 );

  x = zeros ( m1, m2 );
%
%  Load data in the corners only.
%
  i = 1;
  j = 1;
  r = (  i - 1 ) / ( m1 - 1 );
  s = (  j - 1 ) / ( m2 - 1 );
  x(i,j) = cubic_rs ( r, s, 1 );

  i = m1;
  j = 1;
  r = (  i - 1 ) / ( m1 - 1 );
  s = (  j - 1 ) / ( m2 - 1 );
  x(i,j) = cubic_rs ( r, s, 1 );

  i = 1;
  j = m2;
  r = (  i - 1 ) / ( m1 - 1 );
  s = (  j - 1 ) / ( m2 - 1 );
  x(i,j) = cubic_rs ( r, s, 1 );

  i = m1;
  j = m2;
  r = (  i - 1 ) / ( m1 - 1 );
  s = (  j - 1 ) / ( m2 - 1 );
  x(i,j) = cubic_rs ( r, s, 1 );

  x = blend_ij_0d1 ( x, m1, m2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Values interpolated by BLEND_IJ_0D1:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m1
    fprintf ( 1, '  %10f  %10f  %10f  %10f\n', x(i,1:m2) );
  end
%
%  Load data in the edges.
%
  j = 1;
  for i = 1 : m1
    r = (  i - 1 ) / ( m1 - 1 );
    s = (  j - 1 ) / ( m2 - 1 );
    x(i,j) = cubic_rs ( r, s, 1 );
  end

  j = m2;
  for i = 1 : m1
    r = (  i - 1 ) / ( m1 - 1 );
    s = (  j - 1 ) / ( m2 - 1 );
    x(i,j) = cubic_rs ( r, s, 1 );
  end

  i = 1;
  for j = 2 : m2 - 1
    r = (  i - 1 ) / ( m1 - 1 );
    s = (  j - 1 ) / ( m2 - 1 );
    x(i,j) = cubic_rs ( r, s, 1 );
  end

  i = m1;
  for j = 2 : m2 - 1
    r = (  i - 1 ) / ( m1 - 1 );
    s = (  j - 1 ) / ( m2 - 1 );
    x(i,j) = cubic_rs ( r, s, 1 );
  end

  x = blend_ij_1d1 ( x, m1, m2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Values interpolated by BLEND_IJ_1D1:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m1
    fprintf ( 1, '  %10f  %10f  %10f  %10f\n', x(i,1:m2) );
  end
%
%  Compare with BLEND_RS_1D1
%
  for i = 1 : m1
    r = (  i - 1 ) / ( m1 - 1 );
    for j = 1 : m2
      s = (  j - 1 ) / ( m2 - 1 );
      x(i,j) = blend_rs_1dn ( r, s, 1, @cubic_rs );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data blended by BLEND_RS_1DN:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m1
    fprintf ( 1, '  %10f  %10f  %10f  %10f\n', x(i,1:m2) );
  end
%
%  Load all data.
%
  for i = 1 : m1
    for j = 1 : m2
      r = (  i - 1 ) / ( m1 - 1 );
      s = (  j - 1 ) / ( m2 - 1 );
      x(i,j) = cubic_rs ( r, s, 1 );
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact data:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : m1
    fprintf ( 1, '  %10f  %10f  %10f  %10f\n', x(i,1:m2) );
  end

  return
end


text:
function output = calloptidsdp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

model = yalmip2optidsdp(interfacedata);
if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end

if options.savedebug
    save dsdpdebug model
end

solvertime = tic;
[y,fvals,exitflag,stats,X] = dsdp(model.f,model.A,model.b,model.lb,model.ub,model.sdcone,model.y0,model.ops);
solvertime = toc(solvertime);

% Create dual variable in internal format
if options.saveduals
    if any(K.l)
        top = 1;
        D_struc = X{1};
    else
        top = 0;
        D_struc = [];
    end
    if any(K.s)
        for j = 1:length(K.s)
            n = K.s(j);
            Z = triu(ones(n));
            i = find(Z);
            Z(i) = X{j+top};
            Z = Z + Z' - diag(diag(Z));
            D_struc = [D_struc;Z(:)];
        end
    end
else
    D_struc = [];
end

x = y;  % Our notation do not coincide ...
switch exitflag
    case 1
        switch stats.pdflag
            case 1
                problem = 0;
            case 3
                problem = 2;
            case 4
                problem = 1;
            otherwise
                problem = 4;
        end
    case -6
        problem = 6;
    case {-2,-8,-9}
        problem = 4;
    case 7
        problem = 16;
    case {-3,-27}
        problem = 3;
    otherwise
        problem = -1;
end

if options.savesolveroutput
	solveroutput.y = y;
    solveroutput.fvals = fvals;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
    solveroutput.X = X;
else
	solveroutput = [];
end

if options.savesolverinput
	solverinput.model = model;
else
	solverinput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function model = edgesTrainDL(dlPara, varargin )
% Train structured edge detector with deep contour features.
% The code is modified from Structured Edge Detection Toolbox

% get default parameters
dfs={'imWidth',32, 'gtWidth',16, 'nPos',5e5, 'nNeg',5e5, 'nImgs',inf, ...
  'nTrees',10, 'fracFtrs',0.25, 'minCount',1, 'minChild',8, ...
  'maxDepth',64, 'discretize','pca', 'nSamples',256, 'nClasses',2, ...
  'split','gini', 'nOrients',4, 'grdSmooth',0, 'chnSmooth',2, ...
  'simSmooth',8, 'normRad',4, 'shrink',2, 'nCells',5, 'rgbd',0, ...
  'stride',2, 'multiscale',0, 'sharpen',2, 'nTreesEval',4, ...
  'nThreads',4, 'nms',0, 'seed',1, 'useParfor',0, 'modelDir','models/', ...
  'modelFnm','model', 'bsdsDir','BSR/BSDS500/data/'};
opts = getPrmDflt(varargin,dfs,1);
if(nargin==0), model=opts; return; end
% opts.nBatch = dlPara.nBatch;
% opts.nOutputNum = dlPara.nOutputNum;
opts.patch_mean = dlPara.patch_mean;
opts.modelDir = dlPara.modelDir;
opts.selected_dims = dlPara.selected_dims;
% if forest exists load it and return
cd(fileparts(mfilename('fullpath')));
forestDir = [opts.modelDir '/forest/'];
forestFn = [forestDir opts.modelFnm];
if(exist([forestFn '.mat'], 'file'))
  load([forestFn '.mat']); return; end
if(~exist(forestDir,'dir')), mkdir(forestDir); end
% compute constants and store in opts
nTrees=opts.nTrees; nCells=opts.nCells; shrink=opts.shrink;
opts.nPos=round(opts.nPos); opts.nNeg=round(opts.nNeg);
opts.nTreesEval=min(opts.nTreesEval,nTrees);
opts.stride=max(opts.stride,shrink);
imWidth=opts.imWidth; gtWidth=opts.gtWidth;
% imWidth=round(max(gtWidth,imWidth)/shrink/2)*shrink*2;
% opts.imWidth=imWidth; opts.gtWidth=gtWidth;
% nChnsGrad=(opts.nOrients+1)*2; nChnsColor=3;
% if(opts.rgbd==1), nChnsColor=1; end
% if(opts.rgbd==2), nChnsGrad=nChnsGrad*2; nChnsColor=nChnsColor+1; end
% nChns = nChnsGrad+nChnsColor;

% opts.nChnFtrs = imWidth*imWidth*nChns/shrink/shrink;
% opts.nSimFtrs = (nCells*nCells)*(nCells*nCells-1)/2*nChns;
nChns = length(opts.selected_dims);
if(nChns == 0)
    I = zeros([1 1 3], 'single');
    deepFeat = edgesChnsDL2(I, opts);
    nChns = size(deepFeat, 3);
    opts.selected_dims = [1:nChns]';
end
opts.nChns = nChns;
% opts.nTotFtrs = imWidth*imWidth*nChns/shrink/shrink; disp(opts);
opts.nTotFtrs = round(imWidth/shrink)*round(imWidth/shrink)*nChns; disp(opts);
% generate stream for reproducibility of model
stream=RandStream('mrg32k3a','Seed',opts.seed);

% train nTrees random trees (can be trained with parfor if enough memory)
if(opts.useParfor),
    parfor i=1:nTrees,
        trainTree(opts,stream,i);
    end
else
    for i=1:nTrees,
        trainTree(opts,stream,i);
    end;
end

% merge trees and save model
model = mergeTrees( opts );
if(~exist(forestDir,'dir')), mkdir(forestDir); end
save([forestFn '.mat'], 'model', '-v7.3');

end

function model = mergeTrees( opts )
% accumulate trees and merge into final model
nTrees=opts.nTrees; gtWidth=opts.gtWidth;
treeFn = [opts.modelDir '/tree/' opts.modelFnm '_tree'];
for i=1:nTrees
  t=load([treeFn int2str2(i,3) '.mat'],'tree'); t=t.tree;
  if(i==1), trees=t(ones(1,nTrees)); else trees(i)=t; end
end
nNodes=0; for i=1:nTrees, nNodes=max(nNodes,size(trees(i).fids,1)); end
% merge all fields of all trees
model.opts=opts; Z=zeros(nNodes,nTrees,'uint32');
model.thrs=zeros(nNodes,nTrees,'single');
model.fids=Z; model.child=Z; model.count=Z; model.depth=Z;
model.segs=zeros(gtWidth,gtWidth,nNodes,nTrees,'uint8');
for i=1:nTrees, tree=trees(i); nNodes1=size(tree.fids,1);
  model.fids(1:nNodes1,i) = tree.fids;
  model.thrs(1:nNodes1,i) = tree.thrs;
  model.child(1:nNodes1,i) = tree.child;
  model.count(1:nNodes1,i) = tree.count;
  model.depth(1:nNodes1,i) = tree.depth;
  model.segs(:,:,1:nNodes1,i) = tree.hs-1;
end
% remove very small segments (<=5 pixels)
segs=model.segs; nSegs=squeeze(max(max(segs)))+1;
parfor i=1:nTrees*nNodes, m=nSegs(i);
  if(m==1), continue; end; S=segs(:,:,i); del=0;
  for j=1:m, Sj=(S==j-1); if(nnz(Sj)>5), continue; end
    S(Sj)=median(single(S(convTri(single(Sj),1)>0))); del=1; end
  if(del), [~,~,S]=unique(S); S=reshape(S-1,gtWidth,gtWidth);
    segs(:,:,i)=S; nSegs(i)=max(S(:))+1; end
end
model.segs=segs; model.nSegs=nSegs;
% store compact representations of sparse binary edge patches
nBnds=opts.sharpen+1; eBins=cell(nTrees*nNodes,nBnds);
eBnds=zeros(nNodes*nTrees,nBnds);
parfor i=1:nTrees*nNodes
  if(model.child(i) || model.nSegs(i)==1), continue; end %#ok<PFBNS>
  E=gradientMag(single(model.segs(:,:,i)))>.01; E0=0;
  for j=1:nBnds, eBins{i,j}=uint16(find(E & ~E0)'-1); E0=E;
    eBnds(i,j)=length(eBins{i,j}); E=convTri(single(E),1)>.01; end
end
eBins=eBins'; model.eBins=[eBins{:}]';
eBnds=eBnds'; model.eBnds=uint32([0; cumsum(eBnds(:))]);
end

function trainTree( opts, stream, treeInd )
% Train a single tree in forest model.

% location of ground truth
trnImgDir = [opts.bsdsDir '/images/train/'];
trnDepDir = [opts.bsdsDir '/depth/train/'];
trnGtDir = [opts.bsdsDir '/groundTruth/train/'];
trnFeatFir = [opts.modelDir '/feat/'];
if(~exist(trnFeatFir,'dir')), mkdir(trnFeatFir); end
imgIds=dir(trnImgDir); imgIds=imgIds([imgIds.bytes]>0);
imgIds={imgIds.name}; ext=imgIds{1}(end-2:end);
nImgs=length(imgIds); for i=1:nImgs, imgIds{i}=imgIds{i}(1:end-4); end

% extract commonly used options
imWidth=opts.imWidth; imRadius=floor(imWidth/2);
gtWidth=opts.gtWidth; gtRadius=gtWidth/2;
nChns=opts.nChns;
nTotFtrs=opts.nTotFtrs;
selected_dims = opts.selected_dims;
rgbd=opts.rgbd;
nPos=opts.nPos; nNeg=opts.nNeg; shrink=opts.shrink;

% finalize setup
treeDir = [opts.modelDir '/tree/'];
treeFn = [treeDir opts.modelFnm '_tree'];
if(exist([treeFn int2str2(treeInd,3) '.mat'],'file'))
  fprintf('Reusing tree %d of %d\n',treeInd,opts.nTrees); return; end
fprintf('\n-------------------------------------------\n');
fprintf('Training tree %d of %d\n',treeInd,opts.nTrees); tStart=clock;

% set global stream to stream with given substream (will undo at end)
streamOrig = RandStream.getGlobalStream();
set(stream,'Substream',treeInd);
RandStream.setGlobalStream( stream );

% collect positive and negative patches and compute features
fids=sort(randperm(nTotFtrs,round(nTotFtrs*opts.fracFtrs)));
k = nPos+nNeg; nImgs=min(nImgs,opts.nImgs);
ftrs = zeros(k,length(fids),'single');
labels = zeros(gtWidth,gtWidth,k,'uint8'); k = 0;
tid = ticStatus('Collecting data',30,1);
for i = 1:nImgs
  % get image and compute channels
  gt=load([trnGtDir imgIds{i} '.mat']); gt=gt.groundTruth;
  I=imread([trnImgDir imgIds{i} '.' ext]); siz=size(I);
  if(rgbd), D=single(imread([trnDepDir imgIds{i} '.png']))/1e4; end
  if(rgbd==1), I=D; elseif(rgbd==2), I=cat(3,single(I)/255,D); end
  p=zeros(1,4); p([2 4])=mod(4-mod(siz(1:2),4),4);
  if(any(p)), I=imPad(I,p,'symmetric'); end
%   [chnsReg,chnsSim] = edgesChns(I,opts);
%   deepFeat = chnsReg(:, :, 4);
  tic;
if exist([trnFeatFir imgIds{i} '_deepFeat.mat'], 'file');
    load([trnFeatFir imgIds{i} '_deepFeat.mat']);
else
  deepFeat = edgesChnsDL2(I, opts);
  save([trnFeatFir imgIds{i} '_deepFeat.mat'], 'deepFeat');
end
  deepFeat = deepFeat(:, :, selected_dims);
  toc;
  % sample positive and negative locations
  nGt=length(gt); xy=[]; k1=0; B=false(siz(1),siz(2));
  B(shrink:shrink:end,shrink:shrink:end)=1;
  B([1:imRadius end-imRadius:end],:)=0;
  B(:,[1:imRadius end-imRadius:end])=0;
  for j=1:nGt
    M=gt{j}.Boundaries; M(bwdist(M)<gtRadius)=1;
    [y,x]=find(M.*B); k2=min(length(y),ceil(nPos/nImgs/nGt));
    rp=randperm(length(y),k2); y=y(rp); x=x(rp);
    xy=[xy; x y ones(k2,1)*j]; k1=k1+k2; %#ok<AGROW>
    [y,x]=find(~M.*B); k2=min(length(y),ceil(nNeg/nImgs/nGt));
    rp=randperm(length(y),k2); y=y(rp); x=x(rp);
    xy=[xy; x y ones(k2,1)*j]; k1=k1+k2; %#ok<AGROW>
  end
  if(k1>size(ftrs,1)-k), k1=size(ftrs,1)-k; xy=xy(1:k1,:); end
  % crop patches and ground truth labels
  psReg=zeros(round(imWidth/shrink),round(imWidth/shrink),nChns,k1,'single');
  lbls=zeros(gtWidth,gtWidth,k1,'uint8');
%   psSim=psReg;
    ri=round(imRadius/shrink);
    rg=gtRadius;
%     ftrs1 = zeros(k1, nTotFtrs);
  for j=1:k1, xy1=xy(j,:); xy2=xy1/shrink;
    psReg(:,:,:,j)=deepFeat(round(xy2(2))-ri+1:round(xy2(2))+ri,round(xy2(1))-ri+1:round(xy2(1))+ri,:);
%     psSim(:,:,:,j)=chnsSim(xy2(2)-ri+1:xy2(2)+ri,xy2(1)-ri+1:xy2(1)+ri,:);
%     ftrs1(j, :) = deepFeat(round(xy2(2)), round(xy2(1)), :);
    t=gt{xy1(3)}.Segmentation(xy1(2)-rg+1:xy1(2)+rg,xy1(1)-rg+1:xy1(1)+rg);
    if(all(t(:)==t(1))), lbls(:,:,j)=1; else [~,~,t]=unique(t);
      lbls(:,:,j)=reshape(t,gtWidth,gtWidth); end
  end
%   if(0), figure(1); montage2(squeeze(psReg(:,:,1,:))); drawnow; end
%   if(0), figure(2); montage2(lbls(:,:,:)); drawnow; end
  % compute features and store
    ftrs1=reshape(psReg,[],k1)';
%   ftrs1=[reshape(psReg,[],k1)' stComputeSimFtrs(psSim,opts)];
%   ftrs1=reshape(ftrs1,[],k1)';

  ftrs(k+1:k+k1,:)=ftrs1(:,fids);
  labels(:,:,k+1:k+k1)=lbls;
  k=k+k1;
  if(k==size(ftrs,1)), tocStatus(tid,1); break; end
  tocStatus(tid,i/nImgs);
end
if(k<size(ftrs,1)), ftrs=ftrs(1:k,:); labels=labels(:,:,1:k); end

% train structured edge classifier (random decision tree)
pTree=struct('minCount',opts.minCount, 'minChild',opts.minChild, ...
  'maxDepth',opts.maxDepth, 'H',opts.nClasses, 'split',opts.split);
t=labels; labels=cell(k,1); for i=1:k, labels{i}=t(:,:,i); end
pTree.discretize=@(hs,H) discretize(hs,H,opts.nSamples,opts.discretize);
tree=forestTrain(ftrs,labels,pTree); tree.hs=cell2array(tree.hs);
tree.fids(tree.child>0) = fids(tree.fids(tree.child>0)+1)-1;
if(~exist(treeDir,'dir')), mkdir(treeDir); end
save([treeFn int2str2(treeInd,3) '.mat'],'tree'); e=etime(clock,tStart);
fprintf('Training of tree %d complete (time=%.1fs).\n',treeInd,e);
RandStream.setGlobalStream( streamOrig );

end

function ftrs = stComputeSimFtrs( chns, opts )
% Compute self-similarity features (order must be compatible w mex file).
w=opts.imWidth/opts.shrink; n=opts.nCells; if(n==0), ftrs=[]; return; end
nSimFtrs=opts.nSimFtrs; nChns=opts.nChns; m=size(chns,4);
inds=round(w/n/2); inds=round((1:n)*(w+2*inds-1)/(n+1)-inds+1);
chns=reshape(chns(inds,inds,:,:),n*n,nChns,m);
ftrs=zeros(nSimFtrs/nChns,nChns,m,'single');
k=0; for i=1:n*n-1, k1=n*n-i; i1=ones(1,k1)*i;
  ftrs(k+1:k+k1,:,:)=chns(i1,:,:)-chns((1:k1)+i,:,:); k=k+k1; end
ftrs = reshape(ftrs,nSimFtrs,m)';
end

function [hs,segs] = discretize( segs, nClasses, nSamples, type )
% Convert a set of segmentations into a set of labels in [1,nClasses].
persistent cache; w=size(segs{1},1); assert(size(segs{1},2)==w);
if(~isempty(cache) && cache{1}==w), [~,is1,is2]=deal(cache{:}); else
  % compute all possible lookup inds for w x w patches
  is=1:w^4; is1=floor((is-1)/w/w); is2=is-is1*w*w; is1=is1+1;
  kp=is2>is1; is1=is1(kp); is2=is2(kp); cache={w,is1,is2};
end
% compute n binary codes zs of length nSamples
nSamples=min(nSamples,length(is1)); kp=randperm(length(is1),nSamples);
n=length(segs); is1=is1(kp); is2=is2(kp); zs=false(n,nSamples);
for i=1:n, zs(i,:)=segs{i}(is1)==segs{i}(is2); end
zs=bsxfun(@minus,zs,sum(zs,1)/n); zs=zs(:,any(zs,1));
if(isempty(zs)), hs=ones(n,1,'uint32'); segs=segs{1}; return; end
% find most representative segs (closest to mean)
[~,ind]=min(sum(zs.*zs,2)); segs=segs{ind};
% apply PCA to reduce dimensionality of zs
U=pca(zs'); d=min(5,size(U,2)); zs=zs*U(:,1:d);
% discretize zs by clustering or discretizing pca dimensions
d=min(d,floor(log2(nClasses))); hs=zeros(n,1);
for i=1:d, hs=hs+(zs(:,i)<0)*2^(i-1); end
[~,~,hs]=unique(hs); hs=uint32(hs);
if(strcmpi(type,'kmeans'))
  nClasses1=max(hs); C=zs(1:nClasses1,:);
  for i=1:nClasses1, C(i,:)=mean(zs(hs==i,:),1); end
  hs=uint32(kmeans2(zs,nClasses,'C0',C,'nIter',1));
end
% optionally display different types of hs
for i=1:0, figure(i); montage2(cell2array(segs(hs==i))); end
end

repo: DeepContour


text:
function [xi,w]=orthoPolyZerosFromRecur(n,a,b,c,mu0)
%%ORTHOPOLYZEROSFROMRECUR Consider a family of polynomials p_n(x) such that
%              p_n(x) is nth order and the polynomials are orthonormal with
%              respect to some weighting function w(x). This means that
%              integral_lowL^upL w(x)*p_n(x)*p_m(x) dx =1 if m=n and is
%              zero otherwise. All such polynomials can be written in
%              recursive form as
%              p_i(x)=(a(i)*x+b(i))*p_{(i-1)}(x)-c(i)*p_{i-2}(x)
%              starting with p_{-1}(x)=0 and p_0(x)=1. This function finds
%              the zeros of the nth polynomial given expressions for the
%              recursion of the polynomials. Parameters for common
%              orthogonal polynomials are given below. The zeros of such
%              orthonormal polynomials can be used as quadrature points for
%              1D integration from lowL to upL of all polynomials up to
%              order 2*n-1 times the weighting function w(x). The
%              associated quadrature weights for such integrals can also be
%              returned by this function if two outputs are requested.
%
%INPUTS: n The order of the polynomial whose zeros are desired.
%    a,b,c There are two formulations for the input. In the first one,
%          a,b, and c are vectors or function handles such that a(i), b(i),
%          and c(i) give the appropriate values for the recursive
%          formulation of the function values in
%          p_i(x)=(a(i)*x+b(i))*p_{(i-1)}(x)-c(i)*p_{i-2}(x)
%          where the values of i used in a and b go from 1 to n, and
%          indices used in  c(i) go from 2 to n. c(1) is not used.
%          Alternatively, if c is omitted or an empty matrix is passed,
%          then a and b are assumed to be coefficients in the three-term
%          expansion of the form
%          b(i)*p_i(x)=(x-a(i))*p_{(i-1)}(x)-b(i-1)*p_{i-2}(x)
%          and can be vectors or function handles. Indices of a run from 1
%          to n and b run from 1 to (n-1).
%      mu0 This parameter is only needed if the second output of this
%          function is desired. The quadrature weights must be scaled by
%          mu0=integral_lowL^upL w(x) dx. This is that mu0.
%
%OUTPUTS: xi The 1Xn set of zeros of the given set of orthonormal
%            polynomials.
%          w The set of n cubature weights associated with the orthonormal
%            polynomials. This requires the input mu0 be specified.
%
%When xi and w are used for quadrature integration, one can write
%integral_lowL^upL w(x)*f(x) dx=sum_{i=1}^n w(i)*f(xi(i))
%where the equality holds for all polynomials up to degree 2*n-1. This is
%the fundamental theorem of Gaussian quadrature. Thus, finding the zeros of
%polynomials that are orthonormal over a region of interest with respect to
%a particular weighting function can be useful.
%
%The algorithm is that of [1]. Because an eigenvalue decomposition is used,
%the results can be numerically unstable for large n.
%
%Examples of values for common polynomials, are given in [2]. Some are:
%1) Hermite polynomials. These are orthonormal with respect to exp(-x^2)
%   integrated over (-Inf,Inf). The recursion from [3] is
% mu0=sqrt(pi);
% a=@(i)2;
% b=@(i)0;
% c=@(i)2*(i-1);
%
%2) Legendre polynomials. These are orthonormal with respect to w(x)=1 for
%   integration over [-1,+1]. The recursion from [3] is
% mu0=2;
% a=@(i)(2*(i-1)+1)/i;
% b=@(i)0;
% c=@(i)(i-1)/i;
%However, the function GaussLegendrePoints1D should provide
%higher-precision zeros than this function.
%
%3) Laguerre polynomials. These are orthonormal with respect to
%   w(x)=exp(-x) for integration over [0,Inf). The recursion from [3] is
% mu0=1;
% a=@(i)-1/i;
% b=@(i)(2*(i-1)+1)/i;
% c=@(i)(i-1)/i;
%
%4) Associated Laguerre polynomials. These are orthonormal with respect to
%   w(x)=x^c1*exp(-x) for integration over [0,Inf). When considering
%   c1>-1, the recursion from [3] is
% mu0=gamma(1+c1);
% a=@(i)(-1/i);
% b=@(i)(2*i-1+c1)/i;
% c=@(i)(i-1+c1)/i;
%
%5) Gegenbauer polynomials. These are orthonormal with respect to
%   w(x)=(1-x^2)^(c1-1/2) for integration over [-1,1]. When considering
%   c1>-(1/2), the recursion from [3] is
% mu0=(sqrt(pi)*gamma((1/2)+c1))/gamma(1+c1);
% a=@(i)2*(i+c1-1)/i;
% b=@(i)0;
% c=@(i)(i+2*c1-2)/i;
%
%6) Modified Hermite polynomials that are orthogonal with respect to the
%   standard normal N(0,1) distribution integrated over (-Inf,Inf). These
%   are for w(x)=1/sqrt(2*pi)*exp(-x^2/2). The recursion [3] is
% mu0=1;
% a=@(i)1;
% b=@(i)0;
% c=@(i)i-1;
%
%REFERENCES:
%[1] G. H. Golub and J. H. Welsh, "Calculation of Gauss quadrature rules,"
%    Mathematics of Computation, vol. 23, pp. 221-230, 1969.
%[2] Weisstein, Eric W. "Orthogonal Polynomials." From MathWorld--A Wolfram
%    Web Resource. http://mathworld.wolfram.com/OrthogonalPolynomials.html
%[3] Abramowitz, M. and Stegun, I. A. (Eds.). "Orthogonal Polynomials."
%    Table 22.7 in Ch. 22 in Handbook of Mathematical Functions with
%    Formulas, Graphs, and Mathematical Tables, 9th printing. New York:
%    Dover, pp. 782, 1972.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

J=zeros(n,n);
if(nargin<4||isempty(c))
    %If alpha and beta are directly given as a and b.
    for i=1:(n-1)
        J(i,i)=a(i);
        J(i,i+1)=b(i);
        J(i+1,i)=b(i);
    end
    J(n,n)=a(n);
else
    %The first row
    %Rows 1 through n-1
    for i=1:(n-1)
        bi=b(i);
        ai=a(i);
        aiNext=a(i+1);
        ciNext=c(i+1);

        %From Equation 2.2 in [1].
        alpha=-bi/ai;
        beta=sqrt(ciNext/(ai*aiNext));

        J(i,i)=alpha;
        J(i,i+1)=beta;
        J(i+1,i)=beta;
    end
    %The final row
    i=n;
    bi=b(i);
    ai=a(i);
    alpha=-bi/ai;
    J(i,i)=alpha;
end

[V,D]=eig(J);

xi=diag(D)';
w=abs(V(1,:)').^2*mu0;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function linplus_test44 ( )

%*****************************************************************************80
%
%% TEST44 tests R8LT_DET, R8LT_INVERSE, R8LT_MXM, R8LT_RANDOM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST44\n' );
  fprintf ( 1, '  For a matrix in lower triangular storage,\n' );
  fprintf ( 1, '  R8LT_DET computes the determinant.\n' );
  fprintf ( 1, '  R8LT_INVERSE computes the inverse.\n' );
  fprintf ( 1, '  R8LT_MXM computes matrix products.\n' );
  fprintf ( 1, '  R8LT_RANDOM sets a random value.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  [ a, seed ] = r8lt_random ( n, n, seed );

  r8lt_print ( n, n, a, '  Matrix A:' );
%
%  Compute the determinant.
%
  det = r8lt_det ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Determinant is %f\n', det );
%
%  Compute the inverse matrix.
%
  a_inverse = r8lt_inverse ( n, a );

  r8lt_print ( n, n, a_inverse, '  Inverse matrix B:' );
%
%  Check
%
  c = r8lt_mxm ( n, a, a_inverse );

  r8lt_print ( n, n, c, '  Product A * B:' );

  return
end


text:
function [err,Info, opt] = New_HEAD (opt)
%
%   [err,Info, Opt] = New_HEAD (Opt)
%
%Purpose:
%   A function for the likes of Greg Detre.
%   It makes the creation of an AFNI dataset
%   from a 3 or 4 dimensional matrix in matlab a breeze.
%
%
%
%Input Parameters:
%   opt is an options structure with the following fields:
%     .prefix: A string containg the AFNI dataset's prefix
%              .prefix and one of .dimen or .master must
%              be used.
%     .dimen: A string or array containing the dimensions of
%             the dataset (at least 3 dimensions need be present).
%     The following four options are from 3dUndump. See 3dUndump -help
%     .view:    '+orig' (default) or '+acpc' or '+tlrc'
%               If you choose +tlrc, the resultant volume
%               will fit the Talairach box.
%     .master: Name of an existing AFNI dataset that would
%              provide the needed matrix size, orientation etc.
%              Do not mix .master with .dimen
%     .datum:  'byte', 'short' (default), or 'float'
%     .orient: Orientation of data in matrix. Default in 'RAI'
%
%     .tr: a float specifying the TR in seconds.
%          The presence of such a field necessitates a 4 dimensional
%     .scale: 0/1 This is only meaningful with 'short' data where it
%             defaults to 1. It defaults to 0 for other types. It is
%             used when storing float data as shorts to minimize disk
%             use while preserving numeric percision.
%     .Overwrite: y/[n] allow header to be created even if one with
%                 similar name is found on disk
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%   Info: The header structure, given the options specified
%   Opt: A modified version of the input Opt. It will be passed
%        along with the data array to WriteBrik for writing the
%        dataset to disk. (See Examples below)
%
%Examples:
%     Run New_HEAD('test') or see Test_New_HEAD.m for examples.
%
%More Info :
%     AFNI programs 3drefit and 3dAttribute are your friends
%     BrikLoad, BrikInfo, WriteBrik
%     README.attributes
%
%     The function makes heavy use of 3dUndump
%
%This function requires AFNI binaries postdating Feb. 12 2007
%
%
%     Author : Ziad Saad   saadz@mail.nih.gov
%     Date : Fri Feb 9 16:29:51 EST 2007
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'New_HEAD';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
Info = [];

%check on return parameters
if (nargout == 2),
   fprintf(1, 'Notice %s:\nYou are not specifying a return argument for Opt.\n', FuncName);
   fprintf(1, ' Be aware that some fields in Opt might have been modified inside\n');
   fprintf(1, ' %s but these changes are not reflected in the Opt structure you\n', FuncName);
   fprintf(1, ' will be passing to WriteBrik.\n');
   fprintf(1, ' If this makes no sense to you, just use:\n');
   fprintf(1, '    [err, Info, Opt] = New_HEAD(Opt);\n');
end

%work options
if (ischar(opt)), %perhaps test mode
   Test_New_HEAD;
   return;
end

%standardize input
if (isfield(opt,'dimen')),
   if (ischar(opt.dimen)),
      v = str2num(opt.dimen);
      rmfield(opt, 'dimen');
      opt.dimen = v;
   end
end

%create command for -master option or non-tlrc views
if (~isfield(opt,'view') | ~strcmp(opt.view,'+tlrc') | isfield(opt,'master')),
   tmp_suf = '___NeW_hEaD_';
   sopt = '3dUndump -head_only';
   if (isfield(opt,'dimen')),
      sopt = sprintf('%s -dimen %s', sopt, num2str(opt.dimen(1:3)));
   end
   if (isfield(opt,'orient')),
      sopt = sprintf('%s -orient %s', sopt, opt.orient);
   end
   if (isfield(opt,'master')),
      [Status, mPrefix, mView] = PrefixStatus (opt.master);
      sopt = sprintf('%s -master %s', sopt, opt.master);
   else
      mView = '+orig';
   end
   if (~isfield(opt,'datum')), opt.datum = 'short'; end

   if (isfield(opt,'datum')),
      if (~strcmp(opt.datum,'short') & ~strcmp(opt.datum, 'byte') & ~strcmp(opt.datum,'float')),
         fprintf(1,'Error %s:\ndatum option must be one of ''byte'', ''short'', ''float''\nI have %s\n', FuncName, opt.datum);
         return;
      end
      sopt = sprintf('%s -datum %s', sopt, opt.datum);
   end

   if (~isfield(opt,'view')),
      opt.view = '+orig';
   end

   if (isfield(opt,'prefix')),
      [Status, Prefix, View] = PrefixStatus (opt.prefix);
      if (Status < 1 & strcmp(View,opt.view)),
         if (isfield(opt,'overwrite') & ~strcmp(opt.overwrite,'y')),
            fprintf(1,'Error %s:\nLooks like %s exists already.\n',...
                     FuncName, opt.prefix);
            return;
         end
      end
      opt.prefix = Prefix;
      ohead = sprintf('%s%s', tmp_suf, opt.prefix);
      sopt = sprintf('%s -prefix %s%s', sopt, ohead);
   else
      fprintf(1,'Error %s:\nNeed a .prefix option.\n', FuncName);
      return;
   end

   [e,w] = unix(sopt);
   if (e),
      fprintf(1,'Error %s:\nHeader creating command %s failed.\nSee this function''s help and 3dUndump -help\n3dUndump''s output was:\n%s\n', ...
               FuncName, sopt, w);
      New_HEAD_CLEAN(tmp_suf);
      return;
   end

   [err,Info]  = BrikInfo(sprintf('%s%s', ohead,mView));
   New_HEAD_CLEAN(tmp_suf);
else %have +tlrc and no master
   mView = '+tlrc';
   [e,d] = unix('which afni');
   if (e),
      fprintf(1,'Error %s:\nFailed to find afni!\n', FuncName);
      return;
   end
   [e,pt] = GetPath(d);
   templ = sprintf('%s/TT_N27+tlrc.HEAD', pt);
   [e,Info] = BrikInfo(templ);
   if (e),
      fprintf(1,'Error %s:\nFailed to get header of %s!\n', FuncName, templ);
      return;
   end
   %recalculate the origin (res. of template is 1x1x1mm)
   oOrig = Info.ORIGIN;
   oDelta = Info.DELTA;
   oDimen = Info.DATASET_DIMENSIONS;
   %Calculate dimension ratio and adjust so that volume fits in box.
   rat = oDimen(1:3)./opt.dimen(1:3);
   Info.DELTA = Info.DELTA .* rat;
   Info.DATASET_DIMENSIONS(1:3) = opt.dimen(1:3);
   %Now shift the origin so that the final volume still fits in the same box (edge to edge) of TLRC box
   dOrig = 1./2*(Info.DELTA-oDelta);
   Info.ORIGIN = Info.ORIGIN+dOrig;
   %Have to deal with datum
   if (isfield(opt,'datum')),
      if (strcmp(opt.datum,'byte')),
         Info.BRICK_TYPES(1) = 0;
      elseif (strcmp(opt.datum,'short')),
         Info.BRICK_TYPES(1) = 1;
      elseif (strcmp(opt.datum,'float')),
         Info.BRICK_TYPES(1) = 3;
      else
         fprintf(1,'Error %s: Bad data type %s\n', FuncName, opt.datum);
         return;
      end
   else
      opt.datum = 'short'
      Info.BRICK_TYPES(1) = 1;
   end
end

%the scaling option
if (isfield(opt,'scale')),
   if (opt.scale & ~strcmp(opt.datum,'short')),
      fprintf(1,'Warning %s:\n .scale option is only for ''short'' type data.\n Resetting it to 0.\n', FuncName);
      opt.scale = 0;
   end
else
   if (strcmp(opt.datum,'short')),
      opt.scale = 1;
   else
      opt.scale = 0;
   end
end

%take care of prefix business
Info.RootName = sprintf('%s%s', opt.prefix, opt.view);

%take care of view
if (strcmp(mView, opt.view) == 0), %different
   if (strcmp(opt.view,'+orig')) Info.SCENE_DATA(1) = 0;
   elseif (strcmp(opt.view,'+acpc')) Info.SCENE_DATA(1) = 1;
   elseif (strcmp(opt.view,'+tlrc')) Info.SCENE_DATA(1) = 2;
   else
      fprintf(1,'Error %s:\nBad view %s\n', FuncName, opt.view);
      return;
   end
end

%take care of 4th dimen
if (isfield(opt,'dimen')),
   if (length(opt.dimen) > 3),
      Info.DATASET_RANK(2) = opt.dimen(4);
      Info.BRICK_TYPES = Info.BRICK_TYPES(1).*ones(1,opt.dimen(4));
      Info.BRICK_FLOAT_FACS  = Info.BRICK_FLOAT_FACS(1) .*ones(1, opt.dimen(4));
   end
end

%take care of TR
if (isfield(opt,'tr')),
   Info.TAXIS_NUMS(1) = Info.DATASET_RANK(2);
   Info.TAXIS_NUMS(2) = 0; %no time offset for slices at the moment
   Info.TAXIS_NUMS(3) = 77002; %units in seconds for tr (below)

   Info.TAXIS_FLOATS(1) = 0;  %time origin 0
   Info.TAXIS_FLOATS(2) = opt.tr;   %TR in units of Info.TAXIS_NUMS(3)
   Info.TAXIS_FLOATS(3) = 0; %duration of acquisition
   Info.TAXIS_FLOATS(4) = 0; %no time offset please
   Info.TAXIS_FLOATS(5) = 0; %no time offset please

   Info.TAXIS_OFFSETS = zeros(1,Info.TAXIS_NUMS(1));  %no bloody time offset
end

err = 0;
return;

function New_HEAD_CLEAN(sss)
   unix(sprintf('rm -f %s*.HEAD >& /dev/null', sss));
return;


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the BEAM-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BeamCurve()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx = 256;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny = 256;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
N = 1.5*Nx;        % Number of Lagrangian Pts. (2x resolution of Eulerian grid)
a = 0.5;         % Length of beam (only in x-coordinate)
struct_name = 'BeamCurve'; % Name for .vertex, .spring, .beam, .target, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,Lx,a);

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis square;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .beam file!
k_Beam = 5e11; C = 0.0;
print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1.75e8;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', 2 );

    fprintf(target_fid, '%d %1.16e\n', 1, k_Target);
    fprintf(target_fid, '%d %1.16e\n', N, k_Target);

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,Lx,a)

% The immsersed structure is a curved line %
ds = a/(N-1);

xLag(1) = -a/2;%Lx/2-a/2;
yLag(1) = 0.0;
for i=2:N

    xLag(i) = xLag(i-1) + ds;
    x = xLag(i);
    yLag(i) = -2*( (x)^2 - (a/2)^2 );

end

xLag = xLag + Lx/2;
yLag = yLag + Lx/2;

author: nickabattista


text:
function [ArcData, Ang] = ComputeArc(V1,V2,Dist,NumPoints,DirectionFlag)
%COMPUTEARC Computes Arc between two 3D vectors at a specified distance

%Direction of angle, 0 indicates the shortest direction, one indicates to
%go the other direction which results in an angle > 180
if ~exist('DirectionFlag','var')
    DirectionFlag = 0;
end

%Make vectors unit vectors
V1 = V1./norm(V1);
V2 = V2./norm(V2);

%get angle between vectors
Ang = acosd(dot(V1,V2));

if DirectionFlag
    Ang = -1*(360-Ang); %make neg to go the other direction
end

%cross V1 and V2 to get axis or rotation
V3 = cross(V1,V2);
V3 = V3./norm(V3);

%allocate ArcDta
ArcData = zeros(NumPoints,3);

%Set up angles to rotate through
Angs = linspace(0,Ang,NumPoints);

for i=1:NumPoints
    %create rotation matrix to perform rotations
    T = RotateAboutAxis(V3,Angs(i));
    ArcData(i,:) = (T*V1)'*Dist;
end

%set angle back to positive
if DirectionFlag
    Ang = -1*Ang;
end

author: ngageoint


text:
function b = cpbsl ( abd, lda, n, m, b )

%*****************************************************************************80
%
%% CPBSL solves a complex hermitian positive definite band system.
%
%  Discussion:
%
%    The system matrix must have been factored by CPBCO or CPBFA.
%
%    A division by zero will occur if the input factor contains
%    a zero on the diagonal.  Technically this indicates
%    singularity but it is usually caused by improper subroutine
%    arguments.  It will not occur if the subroutines are called
%    correctly and INFO == 0.
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call cpbco(abd,lda,n,rcond,z,info)
%
%      if (rcond is too small .or. info /= 0) then
%        error
%      end if
%
%      do j = 1, p
%        call cpbsl(abd,lda,n,c(1,j))
%      end do
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex ABD(LDA,N), the output from CPBCO or CPBFA.
%
%    Input, integer LDA, the leading dimension of ABD.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer M, the number of diagonals above the main diagonal.
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%

%
%  Solve hermitian(R) * Y = B.
%
  for k = 1 : n
    lm = min ( k - 1, m );
    la = m + 1 - lm;
    lb = k - lm;
    t = abd(la:la+lm-1,k)' * transpose ( b(lb:lb+lm-1) );
    b(k) = ( b(k) - t ) / abd(m+1,k);
  end
%
%  Solve R * X = Y.
%
  for k = n : -1 : 1
    lm = min ( k - 1, m );
    la = m + 1 - lm;
    lb = k - lm;
    b(k) = b(k) / abd(m+1,k);
    t = -b(k);
    b(lb:lb+lm-1) = b(lb:lb+lm-1) + t * transpose ( abd(la:la+lm-1,k) );
  end

  return
end


text:
function dat = alignment(algo,dat)

  [x y]=get_xy(dat);

  if size(x,1)~=size(x,2)
      x=x*x';
  end;
  if size(y,1)~=size(y,2)
      y=y*y';
  end;

  lss=sum(sum( x .* y)) / sqrt(sum(sum(x.^2)) * sum(sum(y.^2)) );

  dat=data([get_name(dat) ' -> alignment='  num2str(lss,4) ],[],lss);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function res = removeObsRobpca(data,i,k,Hsets_min_i,same,factor_ind,csteps)

%REMOVEOBSROBPCA is an auxiliary function to perform cross-validation with ROBPCA,
% RPCR, RSIMPLS, RSIMCA (see cvRobpca.m, cvRpcr.m, cvRsimpls.m).
%
% Input:
%    data : the data set
%       i : the observation that is removed, index with respect to the whole data set.
%       k : the number of principal components that has to be calculated.
%  Hsets_min_i : contains H0_min_i, H1_min_i and Hfreq_min_i as first, second and third row respectively.
%                The h-subsets are implemented by means of indices of the observations in data_min_i, which is
%                is the original data set minus sample i.
%  same : structure :
%           same.value : indicates whether some part of this algorithm can be skipped (= 1) or not ( = 0).
%           same.res : if same.value = 1, then some additional information is needed.
%  factor : optional. default  = 0, then the original consistency factor is used.
%             Else factor = 1, the consistency factor is adapted to the kmax approach.
%  csteps : optional structure:
%            csteps.value : 1 (default, then csteps are performed within robpca), 0 (no csteps are performed within robpca)
%            csteps.number : the number of csteps that need to be performed (default = 2).
%
% Output:
% res is the result structure. It contains:
%   Pk_min_i  : update of the loadingmatrix for a certain k when observation i is deleted.
%   Lk_min_i  : update of the eigenvalues for a certain k when observation i is deleted.
%   muk_min_i : update of the center for a certain k when observation i is deleted.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S. Engelen

rot = [];
center = [];
P1 = [];

n = size(data,1);
p = size(data,2);
h = size(Hsets_min_i,2) + 1;

if nargin < 6
    factor_ind = 0;
    csteps.number = 2;
    csteps.value = 1;
elseif nargin == 6
    csteps.number = 2;
    csteps.value = 1;
end

data_min_i= removal(data,i,0);
mu0 = mean(data);
center = mu0;

H0_min_j = Hsets_min_i(1,:);
H1_min_t = Hsets_min_i(2,:);

if ~same.value
    [PH0_min_i,TH0,LH0_min_i,rH0_min_i,centerX,mu1trafo]=kernelEVD(data_min_i(H0_min_j,:));
    LH0_min_i = diag(LH0_min_i);

    % calculation of the projection
    % res.T2tilde = T2tilde;
    center = mu1trafo;
    rot = PH0_min_i;
    T2tilde = (data_min_i - repmat(mu1trafo,n-1,1))*PH0_min_i;
    T2tilde = T2tilde(:,1:k);
    rot = rot(:,1:k);

    % defining the outputstructure res:
    res.PH0_min_i = PH0_min_i;
    res.LH0_min_i = LH0_min_i;
    res.mu1trafo = mu1trafo;
else
    % defining the input structure input
    res = same.res;
    PH0_min_i = res.PH0_min_i;
    LH0_min_i = res.LH0_min_i;
    mu1trafo = res.mu1trafo;
    T2tilde = (data_min_i - repmat(mu1trafo,n-1,1))*PH0_min_i;
    center = mu1trafo;
    rot = PH0_min_i;
    rot = rot(:,1:k);
    T2tilde = T2tilde(:,1:k);
end

mah = libra_mahalanobis(T2tilde,zeros(1,k),'invcov',1./diag(LH0_min_i(1:k,1:k)));
oldobj = prod(diag(LH0_min_i(1:k,1:k)));P4 = eye(k);

if csteps.value
    oldobj = prod(diag(LH0_min_i(1:k,1:k)));
    for j = 1:csteps.number
        [mahsort, indsort] = sort(mah);
        dataH1 = T2tilde(indsort(1:(h-1)),:);
        [P,T,L,r3,Xm,clmX] = classSVD(dataH1);
        obj = prod(L);
        T2tilde = (T2tilde - repmat(clmX,n-1,1))*P;
        center = center + clmX*rot';
        rot = rot*P;
        mah = libra_mahalanobis(T2tilde,zeros(1,size(T2tilde,2)),'invcov',1./L);
        P4 = P4*P;
        if abs(oldobj - obj) > 1.e-12
            oldobj = obj;
        else
            break
        end
    end
else
    obj = oldobj;
end

% extra reweighting:
if k~=r3
    XRc= data_min_i-repmat(center,n-1,1);
    Xtilde = T2tilde*rot';
    Rdiff = XRc-Xtilde;
    for i=1:n
        odh(i,1)=norm(Rdiff(i,:));
    end
    [m,s]=unimcd(odh.^(2/3),h);
    cutoffodh = sqrt(norminv(0.975,m,s).^3);
    indexset = find(odh<=cutoffodh)';
    [P,Threw,Lrew,rrew,Xmrew,clmX]=kernelEVD(data_min_i(indexset,:));
    center = clmX;
    rot = P(:,1:k);
end
T2tilde = (data_min_i - repmat(center,n-1,1))*rot;

% Perform mcdcov on some H-subsets:
[res_min_i,raw_min_i] = mcdcov(T2tilde,'Hsets',Hsets_min_i,'h',h-1,'ntrial',250,'plots',0,'factor',1);
% perform the last part of ROBPCA:
if raw_min_i.objective < obj
    z = res_min_i;
else
    sortmah = sort(mah);
    if h==n
        factor=1;
    else
        if factor_ind == 1
            factor = sortmah(h-1)/chi2inv((h-1)/size(T2tilde,1),size(T2tilde,2)/2); % adjusted factor.
        else
            factor = sortmah(h-1)/chi2inv((h-1)/size(T2tilde,1),size(T2tilde,2)); % non adjusted factor
        end
    end
    mah = mah/factor;
    weights = mah <= chi2inv(0.975,size(T2tilde,2));
    [center_noMCD,cov_noMCD] = weightmecov(T2tilde,weights);
    mah = libra_mahalanobis(T2tilde,center_noMCD,'cov',cov_noMCD);
    z.flag = (mah <= chi2inv(0.975,size(T2tilde,2)));
    z.center = center_noMCD;
    z.cov = cov_noMCD;
end

covf=z.cov;
centerf=z.center;

% The final PC:
[P6tilde,L6]=eig(covf);
[L6,I]=greatsort(diag(L6));
P6tilde=P6tilde(:,I);
T_min_i=(T2tilde-repmat(centerf,n-1,1))*P6tilde;

Pk_min_i=rot(:,1:k)*P6tilde;
Lk_min_i = L6;

res.Pk_min_i = Pk_min_i;
res.Lk_min_i = Lk_min_i;
res.muk_min_i = center + centerf*rot';

author: marianux

repo: ecg-kit

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


text:
function kern = simwhiteKernParamInit(kern)

% SIMWHITEKERNPARAMINIT SIM-WHITE kernel parameter initialisation.
% FORMAT
% DESC initialises the SIM-White (Single Input Motif - White)
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : David Luengo, 2009

% KERN

if kern.inputDimension > 1
  error('SIM-WHITE kernel only valid for one-D input.')
end

kern.nParams = 3;
kern.decay = 1;
kern.variance = 1;
kern.sensitivity = 1;

kern.delay = 0;
kern.initVal = 1;

kern.transforms.index = [1 2];
kern.transforms.type = optimiDefaultConstraint('positive');

kern.isStationary = false;
kern.positiveTime = true;


text:
function [cIX,gIX] = MakeFoxels(cIX,gIX,M_0,cIX_reg,isWkmeans,clusParams,absIX,i_fish)
% Obtain foxels (functional voxels)

% if isAutoclusWithAllCells,
%     cIX_reg = (1:size(M_0,1))';
% else
%     cIX_reg = cIX;
% end

%% set params
% default:
% thres_reg2 = 0.5;
thres_reg = 0.5;
% thres_merge = 0.5;
% thres_cap = 0.3;
thres_minsize = 10;
numK1 = 20;
avrFxsize = 10;
numK2 = round(length(cIX)/avrFxsize/numK1);

isShowProgress = 1;

% optional override:
if exist('clusParams','var'),
    if ~isempty(clusParams),
        thres_reg = clusParams.reg1; % correlation coeff
        thres_minsize = clusParams.minSize; % number of cells
        numK1 = clusParams.k1;

        isShowProgress = 0;
    end
else

end

    %% 1. k-means (2-step)
    % 1.1. first level k-means
    k1Start = tic;
    if isWkmeans,
        if isShowProgress,
            disp(['kmeans step 1: k = ' num2str(numK1)]);
        end
        gIX = Kmeans_Direct(M_0(cIX,:),numK1);
    else
        [gIX, numK1] = SqueezeGroupIX(gIX);
    end
    k1Time = toc(k1Start);

    % 1.2. divide each of the above clusters again
    if isShowProgress,
        disp(['kmeans step 2: k = ' num2str(numK2)]);
    end
    k2Start = tic;

    gIX_old = gIX;

    % temporarily suppress kmeans warnings
    wid = 'stats:kmeans:FailedToConverge';
    orig_warn = warning('off',wid);

    for i = 1:numK1,
        IX = find(gIX_old == i);
        M_sub = M_0(cIX(IX),:);
        rng('default');
        if numK2<length(IX),
            gIX_sub = kmeans(M_sub,numK2,'distance','correlation');
        else
            gIX_sub = kmeans(M_sub,length(IX),'distance','correlation');
        end
        gIX(IX) = (i-1)*numK2+gIX_sub;
    end
    % restore warnings
    warning(orig_warn);

    k2Time = toc(k2Start);

    %% 2. Regression with the centroid of each cluster
    if isShowProgress,
        disp('regression to make foxels');
    end
    reg0Start = tic;

    Reg = FindCentroid_Direct(gIX,M_0(cIX,:));
    [cIX,gIX,numFoxels] = AllCentroidRegression_SizeThres_direct(M_0(cIX_reg,:),cIX_reg,thres_reg,Reg,thres_minsize/2);

    reg0Time = toc(reg0Start);

%     if exist('absIX','var') && exist('i_fish','var'),
%         clusgroupID = 2;
%         SaveCluster_Direct(cIX,gIX,absIX,i_fish,'foxels',clusgroupID);
%     end

    %% (3. report timing)
    disp(['(Size) nFox:' num2str(numFoxels)]);
    disp(['(Time) k1:' num2str(k1Time) ...
        ' k2:' num2str(k2Time) ...
        ' reg0:' num2str(reg0Time)]);
end

author: xiuyechen

repo: FishExplorer


text:
close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 5;
% Number of signals
S = 4;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
Y = Phi.apply(X);
%  MMV Thresholding solver instance
solver = spx.pursuit.joint.Thresholding(Phi, K);
% Solve the sparse recovery problem
result = solver.solve(Y);
% Solution vector
Z = result.Z;

for s=1:S
    mf.new_figure(sprintf('MMV signal: %d', s));
    subplot(411);
    stem(X(:, s), '.');
    title('Sparse vector');
    subplot(412);
    stem(Z(:, s), '.');
    title('Recovered sparse vector');
    subplot(413);
    stem(abs(X(:, s) - Z(:, s)), '.');
    title('Recovery error');
    subplot(414);
    stem(Y(:, s), '.');
    title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function test_bug1427

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY dataset2files ft_read_header ft_read_data

% the Long64ChannelWithEvents fails because it consists of muliple segments
% filepath = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1427/Long64ChannelWithEvents.mff');

filepath = dccnpath('/home/common/matlab/fieldtrip/data/test/bug1427/NS500Sine6Hz.mff');

filename = dir(filepath);
filename = {filename(~[filename.isdir]).name};

% determine the type for any of the files
typ = cell(size(filename));
for i=1:length(filename)
  typ{i} = ft_filetype(fullfile(filepath, filename{i}));
end

sel = find(strcmp(typ, 'unknown'));
if ~isempty(sel)
  error('the filetype detection failed on %s',  fullfile(filepath, filename{sel(1)}))
end

% read the header for any of the files
hdr = cell(size(filename));
for i=1:length(filename)
  disp(filename{i})
  hdr{i} = ft_read_header(fullfile(filepath, filename{i}));
end

% read some data for any of the files
dat = cell(size(filename));
for i=1:length(filename)
  disp(filename{i})
  dat{i} = ft_read_data(fullfile(filepath, filename{i}), 'begsample', 1, 'endsample', 1000);
end

author: fieldtrip

repo: fieldtrip