text:
%GRIDGEN Grid and mesh generation for geometry objects.
%
%   [ GRID, STATS ] = GRIDGEN( SIN, VARARGIN ) Generates a grid/mesh
%   for the geometry defined by the objects in SIN by calling an
%   external grid generation algorithm. SIN can be either a valid fea
%   problem struct with geometry defined in SIN.geom.objects, a cell
%   array of multiple geometry objects, or a single geometry
%   object. Accepts the following property/value pairs are available
%   with the default grid generation algorithm.
%
%       Property    Value/{Default}           Description
%       -----------------------------------------------------------------------------------
%       gridgen     string {default}          Grid generation algorithm: default, gmsh
%                                                robust (3D), gridgen2d, or triangle (2D)
%       dprim       logical    {true}         Structrured meshing of primitives
%       hmax        scalar/arr {0.1}          Target grid size for subdomains
%       hmaxb       scalar/arr {[]}           Target grid size for boundaries
%       hmaxe       scalar/arr {[]}           Target grid size for edges (3D)
%       grading     scalar     0.3            Mesh grading/growth rate
%       eledge      scalar     1              Elements per edge
%       quad        logical    {false}        Use quad meshing (2D)
%       intb        logical    {true}         Output interior/internal boundaries
%       waitbar     scalar     {0}            Show/hide waitbar
%       fid         scalar     {1}            File identifier for output ([]=no output)
%
%   GRIDGEN specifies which grid generation algorithm to use and calls
%   the corresponding grid generation code (default, gmsh, robust
%   (3D), gridgen2d (2D), or triangle (2D)). Enabling DPRIM with the
%   default grid generator enables structured meshing of single
%   geometry object primitives (and no differing boundary/edge mesh
%   sizing).
%
%   HMAX indicates target grid cell diameters, and is either a numeric
%   scalar prescribing the grid size for the entire geometry, or an
%   array with HMAX values corresponding to individual
%   subdomains. Positive HMAX values uses the minimum mesh size for
%   shared boundaries, while negative applys the mean value.
%
%   HMAXB is analogous to HMAX but related to boundaries
%   (edges/faces). HMAXB can be a single scalar applicable to all
%   boundaries, a numeric array with entries corresponding to
%   individual boundaries.
%
%   HMAXE is analogous to HMAXB but related to 3D edges (only
%   applicable to default grid generator for 3D geometries).
%
%   GRADING specifies the rate at which smaller cells will grow to
%   larger (0 - 1), and ELEDGE optionally prescribes the minimum
%   elements per edge.
%
%   INTB toggels interior/internal boundaries on (default) or off.
%
%
%   The following additional property/value pairs are available with
%   and specific to the GMSH mesh generator.
%
%       Property    Value/{Default}            Description
%       -----------------------------------------------------------------------------------
%       hmaxp       scal/arr {[]}              Target grid size for vertices
%       nref        scalar   {0}               Number of uniform grid refinements
%       algo2       scalar   {2}               2D mesh algorithm (1=MeshAdapt, 2=Automatic,
%                                              5=Delaunay, 6=Frontal, 7=BAMG, 8=DelQuad)
%       algo3       scalar   {1}               3D mesh algorithm (1=Del, 2=New Del, 4=Front
%                                              5=Front Del, 6=Front Hex, 7=MMG3D, 9=R-tree)
%       blayer      struct   {[]}              Data struct for Gmsh 2D boundary layers
%       quad        logical  {false}           Use quad meshing (for 2D)
%       intb        logical  {true}            Output interior/internal boundaries
%       avhb        logical  {true}            Average hmax to boundaries
%       nsm         scalar   {3}               Number of (post) grid smoothing steps
%       palign      scalar   {eps}             Distance tolerance to align point objects
%       tol         scalar   {eps*1e3}         Deduplication tolerance
%       compound    logical  {true}            Use Gmsh compound boundaries
%       mshopt      cell     {}                Cell array of additional Gmsh options
%       mshall      logical  {true}            Output/save all meshed entities
%       mshver      integer  {2}               Gmsh msh file version (1/2/4)
%       verbosity   integer  {5}               Gmsh verbosity/output level
%       syscmd      string   {'default'}       Gmsh system call command
%                                              (default 'gmsh fdir/fname.geo -')
%       fname       string   {'featool_gmsh_UID'}  Gmsh imp/exp file name (root)
%       fdir        string   {tempdir}         Directory to write help files
%       clean       logical  {true}            Delete (clean) Gmsh help files
%
%   NREF (default 0) the number of post uniform grid refinement steps.
%
%   ALGO2 and ALGO3 the Gmsh 2D and 3D mesh generation algorithms.
%
%   QUAD (default 0) toggles Blossom-Quad conversion for 2D geometries.
%
%   The BLAYER flag enables boundary layers for boundaries. May
%   optionally given as array of structs where each entry corresponds
%   to a boundary layer field entry.
%
%   The AVHB logical flag toggles if internal boundaries inheriting
%   HMAX values (when HMAXB is unspecified) should be assigned the
%   smallest HMAX value from neighbouring subdomains, or the mean value
%   (default).
%
%   NSM (default 3) the number of GRIDSMOOTH smoothing steps to perform.
%
%   PALIGN sets a minimum distance over which to re-align mesh
%   vertices to point objects.
%
%   Additional Gmsh options can be provided with the cell array
%   MSHOPT. For example to set the "CharacteristicLengthMax" and "AnisoMax"
%   Gmsh options, MSHOPT could be given as
%
%       {{'Mesh', 'CharacteristicLengthMax', '1'}, {'Mesh', 'AnisoMax', '10'}}
%
%
%   The following additional property/value pairs are available with
%   and specific to the GRIDGEN2D mesh generator.
%
%       Property    Value/{Default}            Description
%       -----------------------------------------------------------------------------------
%       q           scalar   {0.65}            Target quality
%       blayer      logical  {false}           Enable boundary layers
%       avhb        logical  {true}            Average hmax to boundaries
%       nsm         scalar   {3}               Number of (post) grid smoothing steps
%       palign      scalar   {eps}             Distance tolerance to align point objects
%
%   Q (default 0.65) specifies a quality target [0 < q < 1.0].
%
%   The BLAYER flag enables boundary layers for internal boundaries
%   (holes). May optionally given as a vector where the 1st entry
%   specifies relative thickness of stretched layer.
%
%   PALIGN sets a minimum distance over which to re-align mesh
%   vertices to point objects.
%
%
%   The following additional property/value pairs are available with
%   and specific to the TRIANGLE mesh generator.
%
%       Property    Value/{Default}            Description
%       -----------------------------------------------------------------------------------
%       q           scalar   {28}              Minimum target angle (quality)
%       syscmd      string   {'default'}       Triangle system call command (default
%                                              'triangle -I -q%f -j -e -a -A %s.poly')
%       fname       string   {'featool_tri_UID'}  Triangle imp/exp file name (root)
%       fdir        string   {tempdir}         Directory to write help files
%       clean       logical  {true}            Delete (clean) Triangle help files
%
%   Q (default 28 degrees) specifies a minimum target angle (values less than 33 are
%   generally acceptable, while higher values might prevent Triangle convergence).
%
%
%   Examples:
%
%      1) Unit circle with uniform global grid size set to 0.1.
%
%      grid = gridgen( gobj_cone, 'hmax', 0.1 );
%      plotgrid( grid )
%
%      2) Unit square with a finer grid along the top boundary (using Triangle).
%
%      grid = gridgen( gobj_rectangle, 'hmax', 0.5, ...
%                      'hmaxb', [0.5 0.5 0.01 0.5], 'gridgen', 'triangle' );
%      plotgrid( grid )
%
%      3) Domain with curved boundaries meshed with quadrilaterals (using Gmsh).
%
%      geom.objects = {gobj_rectangle() gobj_circle([0 0],.6) gobj_circle([1 1],.3,'C2')};
%      geom = geom_apply_formula( geom, 'R1-C1-C2' );
%      grid = gridgen( geom, 'hmax', 0.1, 'gridgen', 'gmsh', 'quad', true, 'verbosity', 3 );
%      plotgrid( grid )
%
%      4) Two connected subdomains with a shared boundary (using Gridgen2D).
%
%      geom.objects = { gobj_polygon([-2e-3 -8e-3;0 -8e-3;0 -6e-3;0 6e-3;0 8e-3;-2e-3 8e-3]), ...
%      gobj_polygon([0 -6e-3;2e-3 -5e-3;2e-3 4e-3;0 6e-3]) };
%      hmax  = 5e-4;
%      hmaxb = hmax*ones(1,4);
%      hmaxb(9) = hmax/5;
%      grid  = gridgen( geom, 'hmax', hmax, 'hmaxb', hmaxb, 'gridgen', 'gridgen2d' );
%      plotgrid( grid )
%
%      5) Composite component with two subdomains (using Gridgen2D).
%
%      r1 = gobj_rectangle( 0, 0.11, 0, 0.12,  'R1' );
%      c1 = gobj_circle( [ 0.065 0 ],   0.015, 'C1' );
%      c2 = gobj_circle( [ 0.11 0.12 ], 0.035, 'C2' );
%      c3 = gobj_circle( [ 0 0.06 ],    0.025, 'C3' );
%      r2 = gobj_rectangle( 0.065, 0.16, 0.05, 0.07, 'R2' );
%      c4 = gobj_circle( [ 0.065 0.06 ], 0.01, 'C4' );
%      geom.objects = { r1 c1 c2 c3 r2 c4 };
%      geom = geom_apply_formula( geom, 'R1-C1-C2-C3' );
%      geom = geom_apply_formula( geom, 'R2+C4' );
%
%      grid  = gridgen( geom, 'hmax', [0.0025 0.01 0.0025], 'gridgen', 'gridgen2d' );
%      plotgrid( grid )
%
%      6) Complex geometry with several holes and subdomains (using Gridgen2D).
%
%      w = 10e-4; L = 3*w; H = 5*w;
%      p1  = gobj_polygon( [w/10 0;(L-w/4)/2 0;(L-w/4)/2 H;0 H;0 H/3], 'P1' );
%      p2  = gobj_polygon( [(L+w/4)/2 0;L 0;L H-H/3;L H;(L+w/4)/2 H], 'P2' );
%      r1  = gobj_rectangle( (L-w/4)/2, (L+w/4)/2, 0, H, 'R1' );
%      c1  = gobj_circle( [2*w/3 3*w], w/3, 'C1' );
%      c2  = gobj_circle( [2*w/3 2*w], w/3, 'C2' );
%      c3  = gobj_circle( [2*w/3 1*w], w/3, 'C3' );
%      c4  = gobj_circle( [L-w/2 4.5*w], w/8, 'C4' );
%      c5  = gobj_circle( [L-w   4.5*w], w/8, 'C5' );
%      c6  = gobj_circle( [L-w/2 4*w], w/8, 'C6' );
%      c7  = gobj_circle( [L-w   4*w], w/8, 'C7' );
%      c8  = gobj_circle( [L-w/2 3.5*w], w/8, 'C8' );
%      c9  = gobj_circle( [L-w   3.5*w], w/8, 'C9' );
%      c10 = gobj_circle( [L-w/2 3*w], w/8, 'C10' );
%      c11 = gobj_circle( [L-w   3*w], w/8, 'C11' );
%
%      geom.objects = { p1 p2 r1 c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 };
%      geom = geom_apply_formula( geom, 'P1-C1-C2-C3' );
%      geom = geom_apply_formula( geom, 'P2-C4-C5-C6-C7-C8-C9-C10-C11' );
%
%      hmax = w./[5 5 20];     % Set finer mesh size in subdomain 3.
%      hmaxb = zeros(1,21);
%      hmaxb([6 21]) = w/50;   % Set finer mesh size on the in and outlets, boundaries 6 and 21.
%      grid = gridgen( geom, 'hmax', hmax, 'hmaxb', hmaxb, 'gridgen', 'gridgen2d' );
%      plotgrid( grid )
%
%      7) Cone showing the interior y>0 and z<0 (using Gmsh).
%
%      grid = gridgen( gobj_cone, 'hmax', 0.25, 'gridgen', 'gmsh' );
%      plotgrid( grid, 'facecolor', 'none', 'edgecolor', [.7 .7 .7] )
%      plotgrid( grid, 'selcells', '(y>0)&(z<0)' )

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [feat_vector] = SPC_extract_featuresAR(data)
no_range = 0;
if(range(data)<0.01)
    no_range = 1;
    data = randn(size(data));
end

dat1=data(1:floor(end/2));
dat2=data(ceil(end/2)+1:end);
fit = ar_prediction_error(dat1,dat2,9);
AR1 = fit(1);
AR2 = fit(2);
AR3 = fit(3);
AR4 = fit(4);
AR5 = fit(5);
AR6 = fit(6);
AR7 = fit(7);
AR8 = fit(8);
AR9 = fit(9);

feat_vector = [AR1;AR2;AR3;AR4;AR5;AR6;AR7;AR8;AR9];

if(no_range == 1)
    feat_vector = NaN(size(feat_vector));
end

author: alexandrebarachant


text:
function linkText = varlink(var,name)
% link to a variable

if getMTEXpref('generatingHelpMode')
  linkText = name;
else
  %linkText = ['' name ''];
  linkText = ['' name ''];
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function daub10_scale_plot ( n )

%*****************************************************************************80
%
%% DAUB10_SCALE_PLOT plots the DAUB10 scaling function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the recursion level.
%
  x = linspace ( 0.0, 10.0, 801 );

  y = daub10_scale ( n, x );

  plot ( x, y, 'LineWidth', 2 );

  grid on
  xlabel ( '<---X--->' );
  ylabel ( '<---Y--->' );
  title ( sprintf ( 'DAUB10 Scale Function, Recursion level n = %d', n ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% im = appendimages(image1, image2)
%
% Return a new image that appends the two images side-by-side.

function im = appendimages(image1, image2)

% Select the image with the fewest rows and fill in enough empty rows
%   to make it the same height as the other image.
rows1 = size(image1,1);
rows2 = size(image2,1);

if (rows1 < rows2)
     image1(rows2,1) = 0;
else
     image2(rows1,1) = 0;
end

% Now append both images side-by-side.
im = [image1 image2];

author: Sable

repo: mcbench-benchmarks


text:
function u = cos(a)
%COS          Slope cosine cos(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = cos(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  Xxs = hull(a.r(:,indexc),a.r(:,indexr));
  u.s = - sin(Xxs) .* a.s;

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%dir=0 -> q=qab (vb=Cab*va)
%dir=1 -> q=qba (vb=Cba'*va)
function vb=quatrot(qab,va,dir)

va=[0;va(1);va(2);va(3)];
if dir==0
    q=qab;
else
    q=qab;
    q(2:4)=-q(2:4);
end
vr_a=quatmult_v001(q,va,0);
q(2:4)=-q(2:4);
vb=quatmult_v001(vr_a,q,0);
vb=vb(2:4);

repo: ekfmonoslam


text:
function determ = jordan_determinant ( n, alpha )

%*****************************************************************************80
%
%% JORDAN_DETERMINANT returns the determinant of the JORDAN matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real ALPHA, the eigenvalue of the Jordan matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = alpha^n;

  return
end


text:
function [ x, w ] = bdfp_set ( n )

%*****************************************************************************80
%
%% BDFP_SET sets weights for backward differentiation predictor quadrature.
%
%  Discussion:
%
%    A backward differentiation predictor formula is defined for a set
%    of evenly spaced abscissas X(I) with X(1) = 1 and X(2) = 0.  Assuming
%    that the values of the function to be integrated are known at the
%    abscissas, the formula is written in terms of the function value at
%    X(2), and the backward differences at X(2) that approximate the
%    derivatives there.  A backward differentiation predictor formula
%    is equivalent to an Adams-Bashforth formula of the same order.
%
%    The integral:
%
%      Integral ( 0 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * BD**(I-1) F ( 0 ),
%
%    Here, "BD**(I-1) F ( 0 )" denotes the (I-1)st backward difference
%    of F at X = 0, using a spacing of 1.  In particular,
%
%    BD**0 F(0) = F(0)
%    BD**1 F(0) = F(0) - F(-1)
%    BD**2 F(0) = F(0) - 2 * F(-1) + F(-2 )
%
%    The relationship between a backward difference predictor and the
%    corresponding Adams-Bashforth formula may be illustrated for the
%    BDF predictor of order 3:
%
%      BD**0 F(0) + 0.5 * BD**1 F(0) + 5/12 * BD**2 F(0)
%      =            F(0)
%        + 1/2  * ( F(0) -         F(1) )
%        + 5/12 * ( F(0) - 2     * F(-1) +      F(-2) )
%      =  23/12 *   F(0) - 16/12 * F(-1) + 5/12 F(-2)
%
%    which is the Adams-Bashforth formula of order 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Simeon Fatunla,
%    Numerical Methods for Initial Value Problems in Ordinary Differential
%    Equations,
%    Academic Press, 1988.
%
%  Parameters:
%
%    Input, integer N, the order.
%    1 <= N <= 19.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  n_max = 19;

  w_save(1) =                       1.0;
  w_save(2) =                       1.0 /                2.0;
  w_save(3) =                       5.0 /               12.0;
  w_save(4) =                       3.0 /                8.0;
  w_save(5) =                     251.0 /              720.0;
  w_save(6) =                      95.0 /              288.0;
  w_save(7) =                   19087.0 /            60480.0;
  w_save(8) =                    5257.0 /            17280.0;
  w_save(9) =                 1070017.0 /          3628800.0;
  w_save(10) =                  25713.0 /            89600.0;
  w_save(11) =               26842253.0 /         95800320.0;
  w_save(12) =                4777223.0 /         17418240.0;
  w_save(13) =           703604254357.0 /    2615348736000.0;
  w_save(14) =           106364763817.0 /     402361344000.0;
  w_save(15) =          1166309819657.0 /    4483454976000.0;
  w_save(16) =               25221445.0 /         98402304.0;
  w_save(17) =       8092989203533249.0 /    3201186852864.0;
  w_save(18) =         85455477715379.0 /      34237292544.0;
  w_save(19) =   12600467236042756559.0 / 5109094217170944.0;

  if ( n_max < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BDFP_SET - Fatal error!\n' );
    fprintf ( 1, '  N = %d exceeds N_MAX = %d\n', n, n_max );
    error ( 'BDFP_SET - Fatal error!' );
  end

  w(1:n) = w_save(1:n);

  for i = 1 : n
    x(i) = 1 - i;
  end

  return
end


text:
function [handles,msh] = dtiMrMeshAddFGs(handles,msh,showTheseFGs)
%Add the selected fiber groups to the mrMesh window.
%
% [handles,msh] = dtiMrMeshAddFGs(handles,msh,[showTheseFGs])
%
%   handles is from dtiFiberUI (guidata)
%
%   The show selection list is determined from data in the handles
%   structure. You can over-ride that by specifying a vector of 0s and 1s
%   in showTheseFGs. 0 means don't show, 1 means show.  The length is equal
%   to the number of fiber groups.
%
% See also:  dtiFGShowList, dtiRefreshFigure; dtiMrMesh3AxisImage
%
% HISTORY:
% 2004.09.?? Dougherty & Wandell wrote it.
% 2006.11.21 RFD: fixed a couple of little bugs that affected the rendering
% of fibers with few points
%
% (c) Stanford VISTA Team 2005

if ~exist('showTheseFGs','var')||isempty(showTheseFGs),
    showTheseFGs =  dtiFGShowList(handles);
    if(isempty(showTheseFGs)), return; end
end

% This is the list of fiber groups to show
fiberGroups = handles.fiberGroups(showTheseFGs);
defaultAlpha = 255;

stepSizeMm = sqrt(mean(sum(diff(fiberGroups(1).fibers{1},1,2).^2,1)));
skipPoints = round(2./stepSizeMm);

% Display the fiber groups
numGrps = length(fiberGroups);
totalNumActors = 0;
for grpNum=1:numGrps
    fg = fiberGroups(grpNum);
    if(fg.visible)
        fgRgba = uint8(fg.colorRgb);
        sz = size(fgRgba);
        if(sz(2)==3), fgRgba(:,4) = repmat(defaultAlpha,sz(1),1); end
        if(fg.thickness>0)
            % Positive thickness means render as tubes
            clear t;
            t.class = 'mesh';
            [id,s,t] = mrMesh(msh.host, msh.id, 'add_actor', t);
            totalNumActors = totalNumActors+1;
            msh.fiberGroupActors(totalNumActors) = t.actor;

            % This value could change to make really beautiful fibers.  But
            % for exploration, this is fine.  It makes little
            % square/rectangular extrusions.
            t.sides = 5;
            t.radius = fg.thickness;
            t.cap = 1;
            t.points = [];
            % t.points are assembled into a set of tubes separated by a
            % [999,999,999] value.  Inelegant, but OK for now.
            t.color = fgRgba(1,:);
            for ii=1:length(fg.fibers)
                if(size(fg.fibers{ii},2)>skipPoints*10)
                    % We now ensure that we always include the first and last
                    % point. This matters more with metrotrac paths that are
                    % sampled at 2mm.
                    fg.fibers{ii} = horzcat(fg.fibers{ii}(:,1:skipPoints:end-1),fg.fibers{ii}(:,end));
                end
                fg.fibers{ii}(:,end+1) = [999;999;999];
            end
            t.points = horzcat(fg.fibers{:});
            [id,s,r] = mrMesh(msh.host, msh.id, 'tube', t);
        else
            % A negative thickness means render fibers as polylines
            clear t
            t.class = 'polyline';
            t.width = abs(fg.thickness);
            if(sz(1)==1)
                [id,s,r] = mrMesh(msh.host, msh.id, 'add_actor', t);
                totalNumActors = totalNumActors+1;
                msh.fiberGroupActors(totalNumActors) = r.actor(1);
                t.actor = r.actor;
                t.color = uint8(fgRgba);
                % Work-around for fiber color bug in mrMesh. With polylines, it
                % adds a random luminance offset in the range of +/- 8 to the
                % color of each fiber to make individual polylines distinct.
                % But, it donesn't clamp the values at 0,255, so if any of the
                % RGB values are too close to 0,255 the value might wrap and
                % make rainbow colors.
                t.color(t.color(1:3)<8) = 8;
                t.color(t.color(1:3)>247) = 247;
                t.points = [];
                for(ii=1:length(fg.fibers))
                    fg.fibers{ii}(:,end) = [999;999;999];
                end
                t.points = horzcat(fg.fibers{:});
                [id,s,r] = mrMesh(msh.host, msh.id, 'set', t);
            else
                for ii=1:length(fg.fibers)
                    [id,s,r] = mrMesh(msh.host, msh.id, 'add_actor', t);
                    totalNumActors = totalNumActors+1;
                    msh.fiberGroupActors(totalNumActors) = r.actor(1);
                    t.actor = r.actor;
                    if(size(fg.fibers{ii},2)>skipPoints*10)
                        fg.fibers{ii} = horzcat(fg.fibers{ii}(:,1:skipPoints:end-1),fg.fibers{ii}(:,end));
                    end
                    t.points = fg.fibers{ii};
                    t.color = fgRgba(ii,:);
                    t.color(t.color(1:3)<8) = 8;
                    t.color(t.color(1:3)>247) = 247;
                    [id,s,r] = mrMesh(msh.host, msh.id, 'set', t);
                end

            end
        end
    end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rhs=a_rhs(t,a,dummy,phi,L)
rhs= L*a + i*phi'*( (abs(phi*a).^2).*(phi*a) );

author: dynamicslab


text:
function C=multiDimPolyMatMultiply(A,B)
%%MULTIDIMPOLYMATMULTIPLY Given 2D matrices of multivariate polynomials,
%          multiply the matrices to get a resulting matrix containing
%          multivariate polynomials. The multivariate polynomials
%          themselves are represented as hypermatrices, so the "matrices"
%          of these polynomials are represented by cell arrays of
%          hypermatrices holding the polynomial coefficients.
%
%INPUTS:   A An mXr cell array where each element contains a hypermatrix
%            representing a multivariate polynomials. The format of these
%            hypermatrices is discussed below.
%          B An rXn cell array where each element contains a hypermatrix
%            representing a multivariate polynomial.
%
%OUTPUTS: C The mXn-sized product of the matrices of multivariate
%           polynomials, A and B. C is a cell array where each entry is a
%           hypermatrix of coefficients of the multivariate polynomial.
%
%The coefficients of a hypermatrix representing a multivariate polynomial
%coeffs are arranged such that coeffs(a1,a2,a3...an) corresponds to the
%coefficient of an x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus,
%the number of indices coeffs takes is equal to the dimensionality of x
%(not counting singleton dimensions at the end of coeffs). Note that this
%ordering is the reverse that used in the 1D polyval function that is built
%into Matlab. The number of elements for each index in coeffs is the
%maximum order of that dimension +1. The indices of the multivariate
%polynomials in the different elements of the cell arrays must represent
%the same variable.
%
%This just implements the rules for basic matrix multiplication as in
%Chapter 1.11 of [1], using the appropriate functions for multiplying
%(convn) and summing (polySumMultiDim) multivariate polynomials.
%
%EXAMPLE:
%COnsider multiplying a 2X2 matrix of polynomials by a 2X1 vector of
%polynomials.
% A=cell(2,2);
% B=cell(2,1);
% %A multivariate linear polynomial with 3 variables.
% %-3+4*x1+5*x1*x2*x3
% A11=zeros(2,2,2);
% A11(0+1,0+1,0+1)=-3;
% A11(1+1,0+1,0+1)=4;
% A11(1+1,1+1,1+1)=5;
% A{1,1}=A11;
% %Both diagonals equal and just a single third-order term.
% %12*x1^3
% A12=zeros(4,1,1);
% A12(3+1,0+1,0+1)=12;
% A{1,2}=A12;
% A{2,1}=A12;
% %A second order polynomial only with x1 and x3.
% %6+x1^2-2*x3^2
% A22=zeros(3,1,3);
% A22(0+1,0+1,0+1)=6;
% A22(2+1,0+1,0+1)=1;
% A22(0+1,0+1,2+1)=-2;
% A{2,2}=A22;
% %Now the B vector.
% %-6+x1*x2+x1^2
% B1=zeros(2,3);
% B1(0+1,0+1)=-6;
% B1(1+1,1+1)=1;
% B1(2+1,0+1)=1;
% B{1}=B1;
% %-12-4*x1^2*x2^2
% B2=zeros(3,3);
% B2(0+1,0+1)=-12;
% B2(2+1,2+1)=-4;
% B{2}=B2;
% C=multiDimPolyMatMultiply(A,B)
% %One finds that C is a 2X1 cell array where the highest exponent in C{1} is
% %3 and the highest in C{2} is 4. We can evaluate the polynomial at a
% %point. Say
% x=[1;2;3];
% f=@(coeffs)polyValMultiDim(coeffs,x);
% val=cellfun(f,C)
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

m=size(A,1);
r=size(A,2);
n=size(B,2);

C=num2cell(zeros(m,n));

for j=1:n
    for k=1:r
        for i=1:m
           polyProd=convn(A{i,k},B{k,j});
           C{i,j}=polySumMultiDim(C{i,j},polyProd);
        end
    end
end
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:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% 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.

function CstTeflonPTFElossy(mws)

%'@ define material: Teflon (PTFE) (lossy)

material = invoke(mws,'material');
invoke(material,'Reset');
invoke(material,'Name','Teflon (PTFE) (lossy)');
invoke(material,'FrqType','all');
invoke(material,'Type','Normal');
invoke(material,'SetMaterialUnit','GHz','mm');
invoke(material,'Epsilon','2.1');
invoke(material,'Mue','1.0');
invoke(material,'Kappa','0.0');
invoke(material,'TanD','0.0002');
invoke(material,'TanDFreq','10.0');
invoke(material,'TanDGiven','True');
invoke(material,'TanDModel','ConstTanD');
invoke(material,'KappaM','0.0');
invoke(material,'TanDM','0.0');
invoke(material,'TanDMFreq','0.0');
invoke(material,'TanDMGiven','False');
invoke(material,'TanDMModel','ConstKappa');
invoke(material,'DispModelEps','None');
invoke(material,'DispModelMue','None');
invoke(material,'DispersiveFittingSchemeEps','General 1st');
invoke(material,'DispersiveFittingSchemeMue','General 1st');
invoke(material,'UseGeneralDispersionEps','False');
invoke(material,'UseGeneralDispersionMue','False');
invoke(material,'Rho','2200.0');
invoke(material,'ThermalType','Normal');
invoke(material,'ThermalConductivity','0.2');
invoke(material,'HeatCapacity','1.0');
invoke(material,'SetActiveMaterial','all');
invoke(material,'MechanicsType','Isotropic');
invoke(material,'YoungsModulus','0.5');
invoke(material,'PoissonsRatio','0.4');
invoke(material,'ThermalExpansionRate','140');
invoke(material,'Colour','0.94','0.82','0.76');
invoke(material,'Wireframe','False');
invoke(material,'Transparency','0');
invoke(material,'Create');
release(material);
end

repo: CST-MATLAB-API


text:
%kCTRecon 'ML/EM reconstruction of ct-slice '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros CTRecon.pane file
%
% Parameters:
% InputFile: i 'Sinogram', required: 'First Input data object'
% InputFile: t 'BeamTable', optional: 'BeamTable file in raw format'
% Integer: SPos 'Start Projection', default: 0: 'ProjectionNr to start with'
% Integer: PSize 'Projection Size', default: 768: 'integer'
% Integer: PNumber 'Projection Number', default: 1056: 'integer'
% Integer: n 'Iterations', default: 80: 'integer'
% Toggle: q 'No Quarter Shift', default: 0: 'flag'
% Integer: IX 'Size X', default: 128: 'integer'
% Integer: IY 'Size Y', default: 128: 'integer'
% Toggle: divc 'Use Div Cmp', default: 0: 'compute correction by division'
% Toggle: nn 'Without Norms', default: 0: 'flag'
% Toggle: nfn 'No first norm', default: 0: 'do not use norm in forwardprojection'
% Toggle: rn 'Norms from file', default: 0: 'read the norm from the two files'
% OutputFile: o 'Reconstructed', required: 'Resulting output data object'
% OutputFile: c 'Last Correction', required: 'Correction'
% InputFile: Over 'OverralaxTable', optional: 'optional input file'
% InputFile: s 'Start Image', optional: 'optional input file'
% OutputFile: f 'Forward Projection', optional: 'output file'
%
% Example: [o, c, f] = kCTRecon({i, t, Over, s}, {'i','';'t','';'SPos',0;'PSize',768;'PNumber',1056;'n',80;'q',0;'IX',128;'IY',128;'divc',0;'nn',0;'nfn',0;'rn',0;'o','';'c','';'Over','';'s','';'f',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% CTRecon - ML/EM reconstruction of ct-slice
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kCTRecon(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kCTRecon(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'t', '__input';'SPos', 0;'PSize', 768;'PNumber', 1056;'n', 80;'q', 0;'IX', 128;'IY', 128;'divc', 0;'nn', 0;'nfn', 0;'rn', 0;'o', '__output';'c', '__output';'Over', '__input';'s', '__input';'f', '__output'};
maxval={0,1,1,2,2,2,0,2,2,0,0,0,0,0,0,1,1,1};
minval={0,1,1,2,2,2,0,2,2,0,0,0,0,0,0,1,1,1};
istoggle=[0,1,0,0,0,0,1,0,0,1,1,1,1,0,0,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Integer','Integer','Toggle','Integer','Integer','Toggle','Toggle','Toggle','Toggle','OutputFile','OutputFile','InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=2; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function srImage = callSRMethod(slidingWindow, magFactor, srMethod)

    % Assemble name of the SR wrapper function.
    funName = sprintf('superresolve_%s', srMethod.name);

    % Call the desired SR method.
    srImage = feval(funName, slidingWindow, magFactor);

author: thomas-koehler

repo: SupER


text:
function [ freq,amp,phase,dc ] = fourier_series_real( t,x )
% function [ freq,amp,phase,dc ] = fourier_series_real( t,x )
% Fourier series of real signals.
% Written by Yoash Levron, January 2013.
%
% This function computes the fourier series of a signal x(t).
% the amplitudes of the fourier series have the same dimension
% of the original signal, so this function is useful for immediate
% computation of the actual frequency components, without
% further processing.
%
% for example, x(t) = 2 + 3*cos(2*pi*50*t) will result in
% dc value = 2
% frequencies = [50   100  150   200 ...]
% amplitudes   = [3     0      0        0     ...]
% phases         = [0     0      0        0     ...]
%
% x(t) is one cycle of an infinite cyclic signal. The function
% computes the fourier transform of that infinite signal.
% the period of the signal (T) is determined by the length
% of the input time vector, t.
% x(t) must be real (no imaginary values).
%
% The signal x(t) is represented as:
% x(t) = Adc + A1*cos(w*t + ph1) + A2*cos(2*w*t + ph2) + ...
% the function computes the amplitudes, Adc,A1,A2...
% and the phases ph1,ph2,...
%
% T = period of the signal = t(end) - t(1)
% w = basic frequency = 2*pi/T
%
% The function automatically interpolates the original signal
% to avoid aliasing. Likewise, the function automatically determines
% the number of fourier components, and truncates trailing zeros.
%
% inputs:
% t - [sec] time vector. Sample time may vary within the signal.
% x - signal vector. same length as t.
%
% outputs:
% freq - [Hz] frequencies of the fourier series, not including zero.
% amp - amplitudes vector. amp=[A1 A2 A3 ...], not including the DC component.
% phase - [rad/sec] . phases, not including the DC component.
% dc - the DC value (average of the signal).

%%%%%%%%%%% computation %%%%%%%%
rel_tol = 1e-4;  % relative tolerance, to determine trailing zero truncation

if (~isreal(x))
        clc;
        beep;
        disp('fourier_series_real Error:  x(t) must be real.');
        dc = NaN;  amp = NaN;  freq = NaN;  phase = NaN;
        return;
end

t = t-t(1);  % shifting time to zero.
T = t(end);  % period time.
N = 100;  % number of samples
if (mod(N,2) == 1)
    N = N + 1;
end
N = N/2;

ok = 0;
while (~ok)
    N = N*2;  % increase number of samples

    if (N > 10e6)
        clc;
        beep;
        disp('fourier_series_real Error: signal bandwidth seems too high.');
        disp('Try decreasing the sample time in the input time vector t,');
        disp('or increasing the relative tolerance rel_tol');
        dc = NaN;  amp = NaN;  freq = NaN;  phase = NaN;
        return;
    end
    dt = T/N;
    t1 = 0:dt:(T-dt);
    x1 = interp1(t,x,t1,'cubic',0);
    xk = (1/N)*fft(x1);

    dc = abs(xk(1));
    xkpos = xk(2:(N/2));
    xkneg = xk(end:-1:(N/2+2));

    freq = [1:length(xkpos)]/T;  % Hz
    amp = 2*abs(xkpos);
    phase = angle(xkpos);  % rad/sec

    %%% check if enough samples are used.
    %%% if not, try again, with more samples.
    Am = max(amp);
    ii = find((amp(end-10:end)/Am)>rel_tol);
    ok = isempty(ii);
end

% %%% truncate output vectors to remove trailing zeros
Am = max(amp);
ii = length(amp);
while (amp(ii) < Am*rel_tol)
    ii = ii - 1;
end
amp = amp(1:ii);
freq = freq(1:ii);
phase = phase(1:ii);

end


text:
function [g,dgdx,dgdp] = g_GLM_missingData(x,P,u,in)
X = in.X;
md = P(in.md);
b = P(in.b);
X(in.xmd) = md;
g = X*b;
n = length(g);
dgdx = [];
dgdp = zeros(length(P),n);
dgdp(in.b,:) = X';
nmd = length(in.md);
A = zeros(n,nmd);
A([0:nmd-1].*n+mod(in.xmd,n)) = b(ceil(in.xmd./n));
dgdp(in.md,:) = A';

author: MBB-team

repo: VBA-toolbox


text:
function [varargout] = ndgrid(varargin)

%   NDGRID Generation of arrays for N-D functions and interpolation.
%   [X1,X2,X3,...] = NDGRID(x1,x2,x3,...) transforms the domain
%   specified by vectors x1,x2,x3, etc. into arrays X1,X2,X3, etc. that
%   can be used for the evaluation of functions of N variables and N-D
%   interpolation.  The i-th dimension of the output array Xi are copies
%   of elements of the vector xi.
%
%   [X1,X2,...] = NDGRID(x) is the same as [X1,X2,...] = NDGRID(x,x,...).
%
%   For example, to evaluate the function  x2*exp(-x1^2-x2^2-x^3) over the
%   range  -2 < x1 < 2,  -2 < x2 < 2, -2 < x3 < 2,
%
%       [x1,x2,x3] = ndgrid(-2:.2:2, -2:.25:2, -2:.16:2);
%       z = x2 .* exp(-x1.^2 - x2.^2 - x3.^2);
%       slice(x2,x1,x3,z,[-1.2 .8 2],2,[-2 -.2])
%
%   NDGRID is like MESHGRID except that the order of the first two input
%   arguments are switched (i.e., [X1,X2,X3] = NDGRID(x1,x2,x3) produces
%   the same result as [X2,X1,X3] = MESHGRID(x2,x1,x3)).  Because of
%   this, NDGRID is better suited to N-D problems that aren't spatially
%   based, while MESHGRID is better suited to problems in cartesian
%   space (2-D or 3-D).
%
%   This is a drop-in replacement for the matlab version in elmat, which is
%   relatively slow for big grids. FIXME this function still only works up
%   to 5 dimensions
%
%   See also MESHGRID, INTERPN.

%   Copyright(C) 2010, Jan-Mathijs Schoffelen, DCCN
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: ndgrid.m 2956 2011-02-25 19:57:31Z jansch $

if nargin==0
  error('MATLAB:ndgrid:NotEnoughInputs', 'Not enough input arguments.');
end
if nargin==1, varargin = repmat(varargin,[1 max(nargout,2)]); end

ndims = numel(varargin);
switch ndims
case 2
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';

  varargout{1} = x(:, ones2);
  varargout{2} = y(ones1, :);
case 3
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};

  varargout{1} = x(:, ones2, ones3);
  varargout{2} = y(ones1, :, ones3);
  varargout{3} = z(ones1, ones2, :);
case 4
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));
  ones4 = ones(1,numel(varargin{4}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};
  xx   = zeros(1,1,1,numel(varargin{4}));
  xx(:) = varargin{4};

  varargout{1} = x(:, ones2, ones3, ones4);
  varargout{2} = y(ones1, :, ones3, ones4);
  varargout{3} = z(ones1, ones2, :, ones4);
  varargout{4} = xx(ones1, ones2, ones3, :);
case 5
  ones1 = ones(1,numel(varargin{1}));
  ones2 = ones(1,numel(varargin{2}));
  ones3 = ones(1,numel(varargin{3}));
  ones4 = ones(1,numel(varargin{4}));
  ones5 = ones(1,numel(varargin{5}));

  x   = varargin{1}(:);
  y   = varargin{2}(:)';
  z   = zeros(1,1,numel(varargin{3}));
  z(:) = varargin{3};
  xx   = zeros(1,1,1,numel(varargin{4}));
  xx(:) = varargin{4};
  yy   = zeros(1,1,1,1,numel(varargin{5}));
  yy(:) = varargin{5};

  varargout{1} = x(:, ones2, ones3, ones4, ones5);
  varargout{2} = y(ones1, :, ones3, ones4, ones5);
  varargout{3} = z(ones1, ones2, :, ones4, ones5);
  varargout{4} = xx(ones1, ones2, ones3, :,ones5);
  varargout{5} = yy(ones1, ones2, ones3, :,ones5);
otherwise
  error('this version of ndgrid supports inputs up to 5 dimensions');
  %call the ndgrid from elmat
  %FIXME this has to be done
end

author: PatternRecognition

repo: OpenBMI


text:
function c8vec_print_some ( n, x, max_print, title )

%*****************************************************************************80
%
%% C8VEC_PRINT_SOME prints some of a C8VEC.
%
%  Discussion:
%
%    The user specifies MAX_PRINT, the maximum number of lines to print.
%
%    If N, the size of the vector, is no more than MAX_PRINT, then
%    the entire vector is printed, one entry per line.
%
%    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
%    followed by a line of periods suggesting an omission,
%    and the last entry.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries of the vector.
%
%    Input, complex X(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines to print.
%
%    Input, string TITLE, a title.
%
  if ( max_print <= 0 )
    return
  end

  if ( n <= 0 )
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  if ( n <= max_print )

    for i = 1 : n
      fprintf ( 1, '  %6d:  %s\n', i, num2str ( x(i) ) );
    end

  elseif ( 3 <= max_print )

    for i = 1 : max_print-2
      fprintf ( 1, '  %6d: %s\n', i, num2str ( x(i) ) );
    end
    fprintf ( 1, '  ......  ..............\n' );
    i = n;
    fprintf ( 1, '  %6d: %s\n', i, num2str ( x(i) ) );

  else

    for i = 1 : max_print - 1;
      fprintf ( 1, '  %6d: %s\n', i, num2str ( x(i) ) );
    end
    i = max_print;
    fprintf ( 1, '  %6d: %s  ...more entries...\n', i, num2str ( x(i) ) );

  end

  return
end


text:
 function Distance=DistanceToBorder(Position,Corner_position)

%     Position=FullNeuronBehaviorDataSet.FullNeuronBehaviorMatrix{1,1}(:,2:3);
%     Corner_position=FullNeuronBehaviorDataSet.EBC_results{1,4}.QP;
%%
     for i=1:1:size(Position,1)

         Position_now=[squeeze(Position(i,:)),0];
         V1=[Corner_position(1,:),0];
         V2=[Corner_position(2,:),0];
         V3=[Corner_position(3,:),0];
         V4=[Corner_position(4,:),0];

        a = V1 - V2;
        b = Position_now - V2;
        d1 = norm(cross(a,b)) /norm(a);

        a = V2 - V3;
        b = Position_now - V3;
        d2 = norm(cross(a,b)) /norm(a);

        a = V3 - V4;
        b = Position_now - V4;
        d3 = norm(cross(a,b)) /norm(a);

        a = V4 - V1;
        b = Position_now - V1;
        d4 = norm(cross(a,b)) /norm(a);

        Distance(i)=min([d1,d2,d3,d4]);
    end

 end

author: kavli-ntnu


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% File: FRSPC
% Time: Sep. 7th, 2015.
%
% Smooth modeling for incomplete tensor by the PARAFAC decomposition.
% This algorithm is proposed in
%  "Yokota, Tatsuya, et al. "Smooth PARAFAC Decomposition for Tensor Completion." arXiv:1505.06611 (2015)."
%
% minimize || T.*Q - Z.*Q ||_F^2 + penalty_func(U),
%
%    s.t.  Z = [G; U{1}, U{2}, ..., U{N}]
%
% Inputs
% - T       : N-way incomplete tensor
% - Q       : binary tensor which represents elements are available or not (available:1, missing:0)
% - R       : Number of components
% - TV_QV   : it take 'tv' or 'qv' for selecting types of smoothing
% - rho     : N-dimensional vector which represents smoothness of individual modes
% - K       : accelerate parameter for fast optimization (small for fast <--> large for slow, typically 10)
% - SNR     : error threshold via signal-to-noise ratio
% - maxiter : maximum number of iteration
% - tol     : tolerance parameter for convergence evaluation
% - out     : 1 for image completion, 0 for the others
%
% Outputs
% - X       : Results of tensor completion X, where X(Q) = T and X(~Q) = Z
% - Z       : Results of smooth PARAFAC decomposition Z
% - G       : Results of core values G
% - U       : Results of factor matrices U
% - histo   : optimization behavior of error || T.*Q - Z.*Q ||_F^2
% - histo_R : optimization behavior of number of components R
%
% This code was implemented by Tatsuya Yokota
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X Z G U histo histo_R] = FRSPC(T,Q,R,TV_QV,rho,K,SNR,maxiter,tol,out)

  warning off;

  histo=[];
  histo_R=[];
  N  = ndims(T);
  II = size(T);
  NN = sum(Q(:));
  epsiron = 10^(-SNR/10) * sum_all(T(Q).^2);

  %% initialization
  E  = zeros(II);
  X  = zeros(II);
  Z  = zeros(II);
  ONE= zeros(II);

  X(Q) = T(Q);
  X(~Q) = sum(X(:))/NN;
  for n = 1:N
    Pp  = eye(II(n)-1,II(n));
    Pm  = eye(II(n)); Pm(1,:) = [];
    P{n} = Pp - Pm;
    PtP{n} = P{n}'*P{n};
    PtPI{n} = rho(n)*PtP{n} + eye(II(n));
    Pr{n} = inv(rho(n)*PtP{n} + eye(II(n)));
    U{n} = Pr{n}*randn(II(n),R);
    for r = 1:R
      U{n}(:,r) = U{n}(:,r)/norm(U{n}(:,r));
    end
  end

  G = 1e-16*ones(1,R);
  for r = 1:R
    Z = Z + G(r)*outerprod(U,r);
  end
  obj = sum_all((T(Q) - Z(Q)).^2);
  X(~Q) = Z(~Q);
  E = X - Z;

  %% output some figures (only for color image)
  h1 = figure();clf;
  subplot(1,2,1);cla;hold on;
  subplot(1,2,2);
  drawnow;
  if out == 1
    h2 = figure();clf;
    imagesc(uint8(X));drawnow;
    imwrite(uint8(X),['saved/' TV_QV '_iter_0.png']);
  end

  %% start main algorithm
  for iter = 1:maxiter

    [val ID] = sort(abs(G));
    for k = ID(1:min(end,K))

      ONE = G(k)*outerprod(U,k);
      Z = Z - ONE;
      E = E + ONE;

      div = 1;
      for n = 1:N

        if strcmp(TV_QV,'tv')

          u = innerprod_one_exc(E,U,k,n);
          u = u(:);
          %initilization
          a = u/norm(u);
          % main iteration for constrained version
          object = 0.5*G(k)^2*rho(n)*sum(abs(P{n}*a)) - G(k)*a'*u;
          for nn = 1:1000
            df = P{n}'*sign(P{n}*a);
            dL = (0.5*G(k)^2*rho(n)*df - G(k)*u + G(k)^2*a)/G(k)^2;
            al = [0 0.1 0.01 0.001 0.0001 0.00001];
            for ai = 1:length(al)
              a2 = a - al(ai)*dL;
              a2 = a2/norm(a2);
              score(ai) = 0.5*G(k)^2*rho(n)*sum(abs(P{n}*a2)) - G(k)*a2'*u;
            end
            [object2 ai] = min(score);
            a2 = a - al(ai)*dL;
            a2 = a2/norm(a2);
            if abs(object2 - object)/II(n) < 1e-3
              break;
            else
              a = a2;
              object = object2;
            end
          end
          lam = norm(a);
          u = a/lam;
          U{n}(:,k) = u;
          v{n} = u;
          div  = div + rho(n)*sum(abs(P{n}*u));

        elseif strcmp(TV_QV,'qv')

          u = innerprod_one_exc(E,U,k,n);
          u = u(:);
          % initialization
          a = Pr{n}*u;
          a = a/norm(a);
          % main iteration for constrained version
          object = 0.5*G(k)^2*rho(n)*a'*PtP{n}*a - G(k)*a'*u;
          mu = 0.1;
          for nn = 1:1000
            dL = (G(k)^2*rho(n)*PtPI{n}*a - G(k)*u'*a)/G(k)^2;
            al = [0 0.1 0.01 0.001 0.0001 0.00001];
            for ai = 1:length(al)
              a2 = a - al(ai)*dL;
              a2 = a2/norm(a2);
              score(ai) = 0.5*G(k)^2*rho(n)*a2'*PtP{n}*a2 - G(k)*a2'*u;
            end
            [object2 ai] = min(score);
            a2 = a - al(ai)*dL;
            if abs(object2 - object)/II(n) < 1e-3
              break;
            else
              a = a2;
              object = object2;
            end
          end

          u = a2;
          lam = norm(u);
          u = u/lam;
          U{n}(:,k) = u;
          v{n} = u;
          div  = div + rho(n)*u'*PtP{n}*u;

        else
          error('2rd input is ''tv'' or ''qv'' ');
        end

      end

      G(k) = tensor_allprod(E,v,1)/div;
      if G(k) < 0
        G(k) = -G(k);
        U{1}(:,k) = - U{1}(:,k);
      end

      ONE = G(k)*outerprod(U,k);
      E = E - ONE;
      E(~Q) = 0;
      Z = Z + ONE;
      X(~Q) = Z(~Q);

    end

    %% calculate MSE
    obj2 = sum_all(E.^2);

    %% convergence speed
    speed = abs(obj2 - obj)/abs(epsiron - obj2);

    %% checking convergence
    if obj2 < epsiron || abs(obj2-obj)/NN < tol
      break;
    else
      obj = obj2;
      if mod(iter,5)==0
        fprintf('%d:  %f :: %f :: %f :: Pid %d \n',iter,obj2/NN,epsiron/NN,speed,R);
      end
      histo(iter) = obj;
      histo_R(iter,:) = size(G);
    end

    %% output figures
    set(0,'CurrentFigure',h1);
    subplot(1,2,1);cla;hold on;
    plot((histo));
    plot((epsiron)*ones(1,length(histo)))
    grid on;
    set(gca,'YScale','log')
    title('MSE')
    subplot(1,2,2);
    plot(histo_R(:,2));
    title('number of components R')

    drawnow;
    if out == 1
      set(0,'CurrentFigure',h2);
      imagesc(uint8(X));
      title(['Number of R = ' num2str(R)]);
      drawnow;
    end

    if mod(iter,10) == 0 & out == 1
      imwrite(uint8(Z),['saved/' TV_QV '_iter_' num2str(iter) '.png']);
    end

    if mod(iter,100) == 0
      pack;
    end

  end

author: andrewssobral


text:
function airy_ai_prime_values_test ( )

%*****************************************************************************80
%
%% AIRY_AI_PRIME_VALUES_TEST demonstrates the use of AIRY_AI_PRIME_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AIRY_AI_PRIME_VALUES_TEST:\n' );
  fprintf ( 1, '  AIRY_AI_PRIME_VALUES stores values of \n' );
  fprintf ( 1, '  the Airy function A''(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           Ai''(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, aip ] = airy_ai_prime_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, aip );

  end

  return
end


text:
function [SuperFibersGroup, IDsOfFlipped, startCoords, endCoords]=dtiReorientSuperfiberGroups(SuperFibersGroup)

%Given a set of superfiber groups across the participants, checks that all start points and all end points are
%grouped respectively. Flip some groups if needed. Returned all the
%superfibers (including the reoriented ones) and a vector with "1" for
%groups which had to be reoriented (first node -> last)

%Input: 1xX array of SuperFibersGroup
%ER wrote it 08/10/2009
numsubjects=length(SuperFibersGroup);
IDsOfFlipped=[];
for s=1:numsubjects
numNodes(s)=size(SuperFibersGroup(s).fibers{1}, 2);
end

if length(unique(numNodes))>1
    error ('need to have the same number of nodes in each superfibergroup');
end

curves=zeros(3, unique(numNodes), numsubjects);
for i=1:numsubjects
    curves(:, :, i)=SuperFibersGroup(i).fibers{1};
end

%Check that the all the starting points and all the end points are grouped.
%Flip the fibers whose end point groups with the starting points of the
%other.
T = clusterdata([squeeze(curves(:, 1, :))'; squeeze(curves(:, end, :))'], 'maxclust', 2);
Tstart=T(1:end/2); Tend=T(end/2+1:end);

%Take first and last points in every fibers and cluster together. The two
%clouds should be quite separate for a tract with distinct termination
%ROIs. Looking at the labels for the start points, keep fibers
%corresponding to classlabel=1 intact. Swap the node order for the fibers
%corresponding to classlabel=2;
if ~(sum(Tstart==2)==numsubjects | sum(Tstart==1)==numsubjects)
curves(:, :, Tstart==2)=flipdim(curves(:, :, Tstart==2), 2);
display(['Flipped ' num2str(sum(Tstart==2)) ' superfibers of ' num2str(numsubjects)]);
IDsOfFlipped=Tstart==2;
end

startCoords=zeros(3, 1);
endCoords=zeros(3, 1);

for i=1:numsubjects
SuperFibersGroup(i).fibers{1}= curves(:, :, i);
startCoords=startCoords+SuperFibersGroup(i).fibers{1}(:, 1);
endCoords=endCoords+SuperFibersGroup(i).fibers{1}(:, end);
end

startCoords=startCoords./numsubjects;
endCoords=endCoords./numsubjects;


text:
function ori = byEuler(varargin)
% define orientations by Euler angles
%
% Syntax
%   ori = orientation.byEuler(phi1,Phi,phi2,CS,SS) % Bunge convention
%   ori = orientation.byEuler(alpha,beta,gamma,'ZYZ',CS,SS) % Matthies convention
%   ori = orientation.byEuler(phi1,Phi,phi2,'Kocks',CS,SS) % Kocks convention
%
% Input
%  phi1, Phi, phi2 - Euler angles in radiant
%  CS - @crystalSymmetry
%  SS - @specimenSymmetry
%
% Output
%  ori - @orientation
%
% Flags
%  Bunge, ZXZ -
%  ABG, Matthies, ZYZ -
%  Roe -
%  Kocks -
%  Canova -
%
% See also
% orientation/orientation orientation/byMiller orientation/byAxisAngle
% orientation/map

q = euler2quat(varargin{:});

ori = orientation(q,varargin{:});


text:
function  [m, bool] = maxEntConsVector(SInt, printLevel)
% USAGE:
%
%    [m, bool] = maxEntConsVector(SInt, printLevel)
%
% INPUTS:
%    SInt:          `m` x `n` matrix
%    printLevel:    verbose level
%
% OUTPUTS:
%    m:             primal solution
%    bool:          is 1 if `m > 0`
%

[mlt,nlt]=size(SInt);

d1=1e-4;
d2=1;

%pdco parameters
options = pdcoSet;
% options.FeaTol    = 1e-6; %6 March 2010 medium FeaTol
% options.OptTol    = 1e-6; %6 March 2010 medium OptTol
% options.MaxIter   = 200;
% options.Method    = 1;    % 1=Chol  2=QR (more reliable)
% options.mu0       = 0;  % 0 lets pdco decide.  1 or 10 assumes good scaling
% options.StepTol   = 0.9;
% options.StepSame  = 1;
% options.Print     = printLevel-1;
% options.wait      = 0;

x0=ones(mlt,1);
y0=ones(nlt,1);
z0=ones(mlt,1);
xsize=1;
zsize=1;

alpha=1e-6;
%    -----------------------------------------------------------------------
%    pdco.m: Primal-Dual Barrier Method for Convex Objectives (28 Apr 2012)
%    -----------------------------------------------------------------------
%            [x,y,z,inform,PDitns,CGitns,time] = ...
%       pdco(pdObj,pdMat,b,bl,bu,d1,d2,options,x0,y0,z0,xsize,zsize);
%
%     solves optimization problems of the form
%
%        minimize    phi(x) + 1/2 norm(D1*x)^2 + 1/2 norm(r)^2
%          x,r
%        subject to  A*x + D2*r = b,   bl <= x <= bu,   r unconstrained,

[x,y,z,inform,PDitns,CGitns,time] = ...
    pdco(@(x) pdcoObj(x,alpha),SInt',zeros(nlt,1),zeros(mlt,1),inf*ones(mlt,1),d1,d2,options,x0,y0,z0,xsize,zsize);

m=x;
%boolean indicating metabolites involved in the maximal consistent vector
bool=m>0;
end

function [obj,grad,Hess]=pdcoObj(x,alpha)
obj  = alpha*(x'*(log(x) -x));
grad = alpha*log(x);
Hess = alpha./x;
end

author: opencobra

repo: cobratoolbox


text:
clear all; close all; clc;
%% Initialization
%data load
% name={'1_bykim','2_dblee','3_eskim','9_prchoi','10_smkang','12_yskim'}; %session 0
name={'1_bykim','2_dblee','3_eskim','4_jhsuh','5_mgjung','6_mhlee','7_mjkim','8_oykwon','9_prchoi','10_smkang','11_spseo','12_yskim','13_jyha','14_hmyoo','15_jwlee','16_wbsim','17_hwnoh','18_dysull','19_jyhan','20_jwbae','21_syjung','22_dkhan','23_hjwon','24_nkkil','25_khshim'}; % session 1
% name={'1_bykim','2_dblee','3_eskim','4_jhsuh','5_mgjung','7_mjkim','9_prchoi','10_smkang','12_yskim','13_jyha','14_hmyoo','15_jwlee','16_wbsim','17_hwnoh','18_dysull','19_jyhan','20_jwbae','21_syjung','22_dkhan','23_hjwon'}; % session 2
session = {'session1'};
task = {'mi_off','mi_on'};
fs=100;

%pre-processing
channel_index = [1:32];
band = [8 30];
time_interval = [1000 3500];

%feature-extraction
CSPFilter=2;
%classifiaction

%% Data load and mat save
for sub=1:length(name)
    for onoff=1:2
        file3 = ['E:\Users\cvpr\Desktop\StarlabDB_2nd\',name{sub},'\'];
        BMI.EEG_DIR=['E:\Users\cvpr\Desktop\StarlabDB_2nd\',name(sub),'\',session];
        BMI.EEG_DIR=cell2mat(BMI.EEG_DIR);
        file=fullfile(BMI.EEG_DIR, task{onoff});
        marker={'1','left';'2','right'};
        [EEG.data, EEG.marker, EEG.info]=Load_EEG(file,{'device','brainVision';'marker', marker;'fs', fs});
        field={'x','t','fs','y_dec','y_logic','y_class','class', 'chan'};
        CNT=opt_eegStruct({EEG.data, EEG.marker, EEG.info}, field);
        CNTT{sub,onoff} = CNT;
    end
    %     filename1= ['mi_cnt_s1_on'];
    %     save([file3, filename1], 'CNT');
end
%% 'mat' load
% for sub=1:11
%
% %     name={'bykim','dblee','ejlee','eskim','mhlee','prchoi','sbsim','smkang','yelee','yskim'};
%     name={'bykim','eskim','prchoi','smkang','oykwon'};
%     file3 = ['E:\Users\cvpr\Desktop\StarlabDB_2nd\',name{sub},'\'];
%     filename4 = ['mi_cnt_s1_on.mat'];
%
%     cnt = load([file3, filename4]);
%     % MI 19 ch
%     cnt= prep_selectChannels(cnt.CNT, {'Name',{'FC5','FC3','FC1','C5','C3','C1','Cz','CP5','CP3','CP1','CP6','CP4','CP2','C6','C4','C2','FC6','FC4','FC2'}});
%     cnt=prep_filter(cnt, {'frequency', [7 13]});
%     new_smt{sub,1}=prep_segmentation(cnt, {'interval', [750 3500]});
%
% end
%% Pre-processing
clear CNT clear BMI ans EEG field file file3 marker sub
for NUM=1:length(CNTT)
    for onoff=1:2
        CNTch = prep_selectChannels(CNTT{NUM,onoff}, {'Index', channel_index});
        CNTchfilt =prep_filter(CNTch , {'frequency', band});
        all_SMT{NUM,onoff} = prep_segmentation(CNTchfilt, {'interval', time_interval});
        clear CNTch CNTchfilt
    end
end
%% Feature extraction and classification
%% CSP - LDA
for num = 1:length(all_SMT)
    [SMT, CSP_W, CSP_D]=func_csp(all_SMT{num,1},{'nPatterns', CSPFilter});
    FT=func_featureExtraction(SMT, {'feature','logvar'});
    [CF_PARAM]=func_train(FT,{'classifier','LDA'});
    SMT_te=func_projection(all_SMT{num,2}, CSP_W);

    FT_te=func_featureExtraction(SMT_te, {'feature','logvar'});
    [cf_out]=func_predict(FT_te, CF_PARAM);
    [loss out]=eval_calLoss(FT_te.y_dec, cf_out);
    sub_acc(:,num)=1-loss';
    clear SMT CSP_W CSP_D FT CF_PARAM SMT_te FT_Te cf_out loss out
end
    all_acc =sub_acc';
    all_mean_acc = mean(sub_acc,2);
%% FBCSP - LDA

%% BSSFO - LDA


text:
function box = mergeBoxes3d(box1, box2)
%MERGEBOXES3D Merge 3D boxes, by computing their greatest extent.
%
%   BOX = mergeBoxes3d(BOX1, BOX2);
%
%   Example
%   box1 = [5 20 5 30 10 50];
%   box2 = [0 15 0 15 0 20];
%   mergeBoxes3d(box1, box2)
%   ans =
%       0 20 0 30 0 50
%
%
%   See also
%   boxes3d, drawBox3d, intersectBoxes3d
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2010-07-26, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% unify sizes of data
if size(box1,1) == 1
    box1 = repmat(box1, size(box2,1), 1);
elseif size(box2, 1) == 1
    box2 = repmat(box2, size(box1,1), 1);
elseif size(box1,1) ~= size(box2,1)
    error('Bad size for inputs');
end

% compute extreme coords
mini = min(box1(:,1:2:end), box2(:,1:2:end));
maxi = max(box1(:,2:2:end), box2(:,2:2:end));

% concatenate result into a new box structure
box = [mini(:,1) maxi(:,1) mini(:,2) maxi(:,2) mini(:,3) maxi(:,3)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
classdef EvoXBenchProblem < PROBLEM
%EvoXBenchProblem - The base problems of Multiobjective Neural Architecture Search.
%
%   EvoXBenchProblem is a subclass of PROBLEM. In EvoXBench, offspring
%   generated by algorithms are always evaluated on the training set.When
%   saving the results, it is necessary to re-evaluate on the test set,
%   which can also obtain accurate objective values.
%   Before running the EvoXBench test suite, you need to install evoxbench,
%   the installation tutorial address is: https://github.com/EMI-Group/evoxbench.git
%
% EvoXBenchProblem properties:
%   conn            <class>     the TCP/IP object
%   true_eval       <scalar>    true is validation pthase or test phase, false is training phase
%
% EvoXBenchProblem methods:
%   delete          <public>    release resources in Python
%   InitRemoteObj   <public>    initializes the TCP/IP object
%   Setting         <public>    default settings of the problem
%   CalDec          <public>    repair invalid solutions
%   CalObj          <public>    calculate the objective values of solutions
%   Sample          <public>    sample the network architectures
%   GetOptimum      <public>    generate the optimal objective values of the problem
%   GetPF          	<public>    generate the image of the Pareto front
%   CalMetric       <public>    calculate the metric value of a population
%   DrawObj         <public>    display solutions in the objective space
%   SetValid        <public>    set to validation phase
%   SetTest         <public>    set to test phase
%
%   Example:
%
%       platemo()
%
%   displays the GUI of PlatEMO.When users need to save the population on
%   EvoXBench, the objective value of the population needs to be
%   recalculated on the test set. Since the saved results are evaluated on
%   the training set, the decision variables are accurate, but objective
%   values are less accurate.
%
%       pro = @C10MOP1;
%       problem = pro();
%       [decs, objs, cons] = platemo('algorithm',@NSGAII,'problem',pro,'N',100,'maxFE',10000)
%       population = SOLUTION(decs, objs, cons);
%       igd = problem.CalMetric('IGD', population);
%   runs NSGAII with a population size of 100 on C10MOP1 for 10000
%   evaluations and calculate the IGD value. If the users want to save the
%   population, repair invalid solutions and recalculate the objective
%   value.

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    properties (Access = private)
        conn;
        true_eval;
    end
    methods
        %% Release resources in Python
        function delete(obj)
            if ~isempty(obj.conn)
                request.operation = 'delete';
                try
                    obj.conn.writeline(jsonencode(request));
                    obj.conn.flush('output');
                catch
                end
            end
        end
        %% Initializes the TCP/IP object
        function InitRemoteObj(obj, config)
            try
                request.config = config;
                request.operation = 'create';
                obj.conn = tcpclient('127.0.0.1', 9876, 'timeout', 30);
                obj.conn.writeline(jsonencode(request));
                obj.conn.flush('output');
                response = jsondecode(obj.conn.readline());
                assert(strcmp(response.status,'ok'),'Failed to connect to the remote server.');
            catch err
                error('EvoXBench:createFailed','%s\nPlease visit https://github.com/EMI-Group/evoxbench for details.',err.message);
            end
        end
        %% Default settings of the problem
        function Setting(obj, config)
            obj.InitRemoteObj(config);
            query.operation = 'settings';
            obj.conn.writeline(jsonencode(query));
            obj.conn.flush('output');
            response = jsondecode(obj.conn.readline());
            if strcmp(response.status, 'ok')
                basicinfo = response.result;
            else
                exception = MException('TestFunction:queryFailed','Failed to query');
                throw(exception);
            end
            obj.M = basicinfo.n_obj;
            if isempty(obj.D); obj.D = basicinfo.n_var; end
            obj.lower    = reshape(basicinfo.lb, 1, []);
            obj.upper    = reshape(basicinfo.ub, 1, []);
            obj.encoding = ones(1,obj.D);
            obj.true_eval = false;
        end
        %% Repair invalid solutions
        function PopDec = CalDec(obj, PopDec)
            % clip and round, because evoxbench uses discrete encoding
            PopDec = min(PopDec, repmat(obj.upper, size(PopDec,1), 1));
            PopDec = max(PopDec, repmat(obj.lower, size(PopDec,1), 1));
            PopDec = round(PopDec);

            % deduplicate
            TotalLen = size(PopDec, 1);
            [~, ia] = unique(PopDec, 'rows');
            NonDupLen = size(ia, 1);
            DupLen = TotalLen - NonDupLen;
            DupIndex = true(TotalLen, 1);
            DupIndex(ia) = false;
            PopDec(DupIndex, :) = obj.Sample(DupLen);
        end
        %% calculate the objective values
        function PopObj = CalObj(obj, x)
            query.operation = 'query';
            query.encoding = x;
            query.true_eval = obj.true_eval;
            obj.conn.writeline(jsonencode(query));
            obj.conn.flush('output');
            response = jsondecode(obj.conn.readline());
            if strcmp(response.status, 'ok')
                PopObj = response.result;
            else
                exception = MException('EvoXBench:queryFailed','Failed to query');
                throw(exception);
            end
        end
        %% Sample the network architectures
        function R = Sample(obj, n_samples)
            query.operation = 'sample';
            query.n_samples = n_samples;
            obj.conn.writeline(jsonencode(query));
            obj.conn.flush('output');
            response = jsondecode(obj.conn.readline());
            if strcmp(response.status, 'ok')
                R = response.result;
            else
                exception = MException('EvoXBench:queryFailed','Failed to query');
                throw(exception);
            end
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj, N)
            query.operation = 'pareto_front';
            query.n_samples = N;
            obj.conn.writeline(jsonencode(query));
            obj.conn.flush('output');
            response = jsondecode(obj.conn.readline());
            if strcmp(response.status, 'ok')
                R = response.result;
            else
                exception = MException('EvoXBench:queryFailed','Failed to query');
                throw(exception);
            end
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            % -1 will get all points
            R = obj.GetOptimum(-1);
        end
        function score = CalMetric(obj,metName,Population)
        %CalMetric - calculate the metric value of a population
        %
        %   value = obj.CalMetric(Met,P) returns the metric value of a
        %   population P, where Met is a string denoting the name of a
        %   metric function. If Met is 'HV', 'IGD' or 'normalized_hv', the
        %   metric value is calculate through the Python side. Before the
        %   metric value is calculated, the test phase is set and the
        %   objective value is recalculated. Also, invalid solutions are
        %   removed.
        %
        %   Example:
        %       value = Problem.CalMetric('HV',Population);

            obj.true_eval = true;
            PopDec = obj.CalDec(Population.decs());
            PopObj = obj.CalObj(PopDec);
            PopCon = obj.CalCon(PopDec);
            invalidSolutions = any(PopObj == inf, 2);
            validSolutions = ~invalidSolutions;
            PopDec = PopDec(validSolutions, :);
            indicator = lower(metName);
            if strcmp(indicator, 'hv') || strcmp(indicator, 'igd') || strcmp(indicator, 'normalized_hv')
                query.operation = 'calc_perf_indicator';
                query.inputs = PopDec;
                query.indicator = indicator;
                obj.conn.writeline(jsonencode(query));
                obj.conn.flush('output');
                response = jsondecode(obj.conn.readline());
                if strcmp(response.status, 'ok')
                    score = response.result;
                else
                    exception = MException('EvoXBench:queryFailed','Failed to query');
                    throw(exception);
                end
            else
                Population = SOLUTION(PopDec,PopObj,PopCon);
                score = feval(metName,Population,obj.optimum);
            end
        end
        function DrawObj(obj,Population)
            %DrawObj - Display a population in the objective space.
            %
            %   To display the results of test phase in real time in the
            %   GUI, set to validation phase. Then recalculate the
            %   objective value and create a new SOLUTION object.
            %
            %	This function is usually called by the GUI.
            obj.true_eval = true;
            PopDec = obj.CalDec(Population.decs());
            PopObj = obj.CalObj(PopDec);
            PopCon = obj.CalCon(PopDec);
            Population = SOLUTION(PopDec,PopObj,PopCon);

            ax = Draw(Population.objs,{'\it f\rm_1','\it f\rm_2','\it f\rm_3'});
            if ~isempty(obj.PF)
                if ~iscell(obj.PF)
                    if obj.M == 2
                        plot(ax,obj.PF(:,1),obj.PF(:,2),'-k','LineWidth',1);
                    elseif obj.M == 3
                        plot3(ax,obj.PF(:,1),obj.PF(:,2),obj.PF(:,3),'-k','LineWidth',1);
                    end
                else
                    if obj.M == 2
                        surf(ax,obj.PF{1},obj.PF{2},obj.PF{3},'EdgeColor','none','FaceColor',[.85 .85 .85]);
                    elseif obj.M == 3
                        surf(ax,obj.PF{1},obj.PF{2},obj.PF{3},'EdgeColor',[.8 .8 .8],'FaceColor','none');
                    end
                    set(ax,'Children',ax.Children(flip(1:end)));
                end
            elseif size(obj.optimum,1) > 1 && obj.M < 4
                if obj.M == 2
                    plot(ax,obj.optimum(:,1),obj.optimum(:,2),'.k');
                elseif obj.M == 3
                    plot3(ax,obj.optimum(:,1),obj.optimum(:,2),obj.optimum(:,3),'.k');
                end
            end
            obj.true_eval = false;
        end
        %% Set to validation phase
        function SetValid(obj)
            obj.true_eval = true;
        end
        %% Set to test phase
        function SetTest(obj)
            obj.true_eval = true;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% MOVE_INDR - Removes the indices indc from the indices inda
%             and adds them to the reserve vector index list
%             if necessary.
%
% Syntax: [inda,i] = move_indr(inda,indc)
%
%      i: when i > 0, the i-th element was removed from ind{RESERVE}
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function [inda,removed_i] = move_indr(inda,indc)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLABELED = 4;

% define global variables
global g;							% partial derivatives of cost function w.r.t. alpha coefficients
global ind;							% cell array containing indices of margin, error, reserve and unlearned vectors
global max_reserve_vectors;   % maximum number of reserve vectors

removed_i = 0;
num_RVs_orig = length(ind{RESERVE});

% shift indc from inda to ind{RESERVE}
[inda,ind{RESERVE}] = move_ind(inda,ind{RESERVE},indc);

% if we need to remove some reserve vectors
if (length(ind{RESERVE}) > max_reserve_vectors)

   % sort g(ind{RESERVE})
   [g_sorted,i] = sort(g(ind{RESERVE}));

   % reserve vectors that need to be removed
   removed = i(max_reserve_vectors+1:length(i));

   % find any original reserve vectors that need to be removed
   k = find(removed <= num_RVs_orig);
   if (length(k) > 0)
      removed_i = removed(k);
   end;

   % remove the necessary reserve vectors
   ind{RESERVE}(removed) = [];

end;


text:
function [SampEn, ApEn] = EvalEntropyMetrics(rr, t_rr, m, r, HRVparams, tWin, sqi)

% [SampEn, ApEn] = EvalEntropyMetrics(rr, t_rr, m, r, HRVparams, tWin, sqi)
%
% Overview
%	Calculates Sample Entropy and Approximate Entropy values of input data.
%
% Input
%   rr          - (seconds) rr intervals
%   t_rr        - (seconds) time of the rr intervals
%   m           - pattern length
%   r           - radius of similarity
%   HRVparams   - struct of settings for hrv_toolbox analysis
%   sqi         - Signal Quality Index; Requires a matrix with at least two
%                 columns. Column 1 should be timestamps of each sqi
%                 measure, and Column 2 should be SQI on a scale from 0 to 1
%   tWin        - vector containing the starting time of each
%                                windows (in seconds)
% Output
%   SampEn      -
%   ApEn        -
%
% Written by Giulia Da Poian <giulia.dap@gmail.com>
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%

if nargin<3
    Error('Not enough input arguments!');
end
if nargin<4 || isempty(HRVparams)
    windowlength = length(data);
    SQI_th = 0.9;          % SQI threshold
    WinQuality_th = 0.20;  % Low quality windows threshold
else
    windowlength = HRVparams.windowlength;
    SQI_th = HRVparams.sqi.LowQualityThreshold;        % SQI threshold
    WinQuality_th = HRVparams.RejectionThreshold;  % Low quality windows threshold
end
if nargin<5 || isempty(tWin)
    tWin = 0;
end
if nargin <6 || isempty(sqi)
    sqi(:,1) = t_rr;
    sqi(:,2) = ones(length(t_rr),1);
end

% Preallocation (all NaN)
SampEn = nan(length(tWin),1);
ApEn = nan(length(tWin),1);

% What Sample Entropy function?
SampEnType = 'Maxim'; % Initialize default SampEn method
if windowlength < 34000; SampEnType = 'Fast'; end

% Loop through each window of RR data
for iWin = 1:length(tWin)
    if ~isnan(tWin(iWin))
        % Isolate data in this window
        sqi_win = sqi( sqi(:,1) >= tWin(iWin) & sqi(:,1) < tWin(iWin) + windowlength,:);
        nn_win = rr( t_rr >= tWin(iWin) & t_rr < tWin(iWin) + windowlength );
        lowqual_idx = find(sqi_win(:,2) < SQI_th);  % Analysis of SQI for the window
        % If enough data has an adequate SQI, perform the calculations
        if numel(lowqual_idx)/length(sqi_win(:,2)) < WinQuality_th

            nn_win = zscore(nn_win);  % normalization of the signal that replace the common
                                      % practice of expressing the tolerance as r times the std
            % 1. Sample Entropy
            switch SampEnType
                case 'Fast'
                   SampEn(iWin) = fastSampen(nn_win, m, r);
                otherwise
                   SampEn(iWin) = sampenMaxim(nn_win, m, r);
            end
            % 2. Approximate Entropy
            ApEn(iWin) =  ApproxEntropy( nn_win, m, r);
        end
    end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function W = weight_conversion(W, wcm)
% WEIGHT_CONVERSION    Conversion of weights in input matrix
%
%   W_bin = weight_conversion(W, 'binarize');
%   W_nrm = weight_conversion(W, 'normalize');
%   L = weight_conversion(W, 'lengths');
%   W_fix = weight_conversion(W, 'autofix');
%
%   This function may either binarize an input weighted connection matrix,
%   normalize an input weighted connection matrix, convert an input
%   weighted connection matrix to a weighted connection-length matrix, or
%   fix common connection problems in binary or weighted connection matrices.
%
%       Binarization converts all present connection weights to 1.
%
%       Normalization rescales all weight magnitudes to the range [0,1] and
%   should be done prior to computing some weighted measures, such as the
%   weighted clustering coefficient.
%
%       Conversion of connection weights to connection lengths is needed
%   prior to computation of weighted distance-based measures, such as
%   distance and betweenness centrality. In a weighted connection network,
%   higher weights are naturally interpreted as shorter lengths. The
%   connection-lengths matrix here is defined as the inverse of the
%   connection-weights matrix.
%
%       Autofix removes all Inf and NaN values, remove all self connections
%   (sets all weights on the main diagonal to 0), ensures that symmetric matrices
%   are exactly symmetric (by correcting for round-off error), and ensures that
%   binary matrices are exactly binary (by correcting for round-off error).
%
%   Inputs: W           binary or weighted connectivity matrix
%           wcm         weight-conversion command - possible values:
%                           'binarize'      binarize weights
%                           'normalize'     normalize weights

%                           'lengths'       convert weights to lengths
%                           'autofix'       fixes common weights problems
%
%   Output: W_          output connectivity matrix
%
%
%   Mika Rubinov, U Cambridge, 2012

%   Modification History:
%   Sep 2012: Original
%   Jan 2015: Added autofix feature.

switch wcm
    case 'binarize'
        W=double(W~=0);         % binarize
    case 'normalize'
        W=W./max(abs(W(:)));    % scale by maximal weight
    case 'lengths'
        E=find(W);
        W(E)=1./W(E);           % invert weights
    case 'autofix'
        % clear diagonal
        n = length(W);
        W(1:n+1:end)=0;

        % remove Infs and NaNs
        idx = isnan(W) | isinf(W);
        if any(any(idx));
            W(idx)=0;
        end

        % ensure exact binariness
        U = unique(W);
        if numel(U)>1
            idx_0 = abs(U)<1e-10;
            idx_1 = abs(U-1)<1e-10;
            if all(idx_0 | idx_1)
                W(idx_0)=0;
                W(idx_1)=1;
            end
        end

        % ensure exact symmetry
        if ~isequal(W,W.');
            if max(max(abs(W-W.'))) < 1e-10;
                W=(W+W).'/2;
            end
        end
    otherwise
        error('Unknown weight-conversion command.')
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function plot_setPrmAndPlot(hFig, idxPC, varargin)

% Actual parameter
p = hFig.UserData.PC{idxPC}.U.p;

for i = 1:2:numel(varargin)

    prm = varargin{i};
    val = varargin{i+1};

    % Color --------------------------------------------------------------------

    if strcmpi(prm, 'Color')
        if ischar(val)
            % Convert color defined as char to RGB value
            % Note: this is necessary, as:
            %       - Unicolor:  colors have to be defined as RGB values
            %       - Attribute: attributes have to be defined as char
            [row, ~] = find(strcmp(hFig.UserData.colors, val));
            if ~isempty(row) % a color name was selected
                p.Color = hFig.UserData.colors{row,1};
            % Temp !!!!! - start
            elseif strcmpi(val(1:2), 'A.')
                p.Color = val;
                if ~any(strcmp(varargin, 'CAxisLim'))
                    p.CAxisLim = [];
                end
            % Temp !!!!! - end
            else % 'RGB' or 'by PC'
                p.Color = val;
            end
        elseif numel(val) == 3
            p.Color = val;
        else
            plot_throwError(prm); return;
        end
    end

    % MarkerSize ---------------------------------------------------------------

    if strcmpi(prm, 'MarkerSize')
        if isnumeric(val) && isscalar(val) && val>0
            p.MarkerSize = val;
        else
            plot_throwError(prm); return;
        end
    end

    % MaxPoints ----------------------------------------------------------------

    if strcmpi(prm, 'MaxPoints')
        if isnumeric(val) && isscalar(val) && val>0
            p.MaxPoints = val;
        else
            plot_throwError(prm); return;
        end
    end

    % ColormapName -------------------------------------------------------------

    if strcmpi(prm, 'ColormapName')
        if ischar(val)
            p.ColormapName = val;
        else
            plot_throwError(prm); return;
        end
    end

    % CAxisLim -----------------------------------------------------------------

    if strcmpi(prm, 'CAxisLim')
        if isnumeric(val) && numel(val)==2
            p.CAxisLim = val;
        else
            plot_throwError(prm); return;
        end
    end

    % ToDo: error if non-valid prm name is defined

end

hFig.UserData.PC{idxPC}.U.p = p;

% Update figure and axes -------------------------------------------------------

% Axes color
hAxes = findobj('Parent', hFig, 'Type', 'Axes');
row = strcmpi(hFig.UserData.options.AxesColor, hFig.UserData.colors(:,3));
if any(row), hAxes.Color = hFig.UserData.colors{row,1}; end

% Figure color
row = strcmpi(hFig.UserData.options.FigureColor, hFig.UserData.colors(:,3));
if any(row), hFig.Color = hFig.UserData.colors{row,1}; end

% Help text
hFig.UserData.hHelptext.BackgroundColor = hFig.UserData.options.FigureColor;

% Plot! ------------------------------------------------------------------------

plot_updatePlot(hFig, idxPC);

end

author: pglira


text:
function I=wave_inv_transform(x,ld)
[m n]=size(x);
% [ld hd lr hr]=wfilters('sym8');
t=0:(length(ld)-1);
hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
Wn=zeros(n);
Wm=zeros(m);
z=length(ld);
% l=ld;
% l(end:-1:1)=ld;
% ld=l;
j=1;
for i=1:n/2

        if j+z-1>n
            a=j+z-1-n;
            Wn(i,1:a)=ld(z-a+1:end);
            Wn(i,j:end)=ld(1:z-a);
        else
             Wn(i,j:j+z-1)=ld;
        end
        j=j+2;
end
j=1;
for i=n/2+1:n

        if j+z-1>n
            a=j+z-1-n;
            Wn(i,1:a)=hd(z-a+1:end);
            Wn(i,j:end)=hd(1:z-a);
        else
             Wn(i,j:j+z-1)=hd;
        end
        j=j+2;
end
Wn;
j=1;
for i=1:m/2

         if j+z-1>m
             a=j+z-1-m;
             Wm(i,1:a)=ld(z-a+1:end);
             Wm(i,j:end)=ld(1:z-a);
         else
              Wm(i,j:j+z-1)=ld;
         end
         j=j+2;
 end
 j=1;
 for i=m/2+1:m

         if j+z-1>m
             a=j+z-1-m;
             Wm(i,1:a)=hd(z-a+1:end);
             Wm(i,j:end)=hd(1:z-a);
         else
              Wm(i,j:j+z-1)=hd;
         end
         j=j+2;
 end
 Wm;
x=double(x);
 c=Wm'*x;
 I=c*Wn;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [Q,J,F,T,I,V,E] = spinv(A,varargin)
  % SPINV Compute a sparse factorization of the Moore-Penrose pseudoinverse A⁺
  % of a sparse matrix A.  This is useful for solving
  % under-and-over-constrained systems Ax = b: where A is m by n with
  % rank r < min(m,n). Then to solve the system:
  %
  %     x = E*(V*(I'*(T'\(F'*(J'*(Q'*b))))))
  %
  % Note: Store these sparse factors and sovle using the routine above **with
  % parentheses**. Do not precompute A⁺ as it will be dense.
  %
  % Inputs:
  %   A  m by n sparse
  %   Optional:
  %     'Epsilon' followed by epsilon used to determine rank from qr {1e-15}
  % Outputs:
  %   Q  m by m sparse orthogonal matrix
  %   J  m by r sparse rectangular identity matrix
  %   F  r by r sparse permutation matrix
  %   T  r by r sparse triangular matrix (diagonal?)
  %   I  r by n sparse rectangular identity matrix
  %   V  n by n sparse orthogonal matrix
  %   E  n by n sparse permutation matrix
  %
  % Example:
  %   % over and under determined system of 3 variables
  %   A = [1 0 0;1 0 0;0 1 1;0 1 1];
  %   b = [4;5;6;7];
  %   % Compute facotrization of sparse pseudo inverse
  %   [Q,J,F,T,I,V,E] = spinv(A);
  %   % Compute minimum norm best fit
  %   xs = E*(V*(I'*(T'\(F'*(J'*(Q'*b))))));
  %   % Compare to Tikhonov regularization
  %   xt = (A'*A+1e-10*speye(3,3))\(A'*b);
  %

  % default values
  epsilon = 1e-15;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'Epsilon'}, {'epsilon'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end
  % http://www.irisa.fr/sage/wg-statlin/WORKSHOPS/LEMASSOL05/SLIDES/QR/Guyomarch.pdf

  % dimensions of A
  m = size(A,1);
  n = size(A,2);
  % QR factorization of AE
  [Q,R,E] = qr(A);
  % We have that:
  %     AE = QR
  % rank
  r = find(any(abs(R)>epsilon,2),1,'last');
  % Only keep upper, non-zero block
  R1R2 = R(1:r,:);
  J = speye(m,r);
  % We have that:
  %    AE = Q*J*R1R2
  % QR factorization of R1R2'
  [V,Tfull,F] = qr(R1R2');
  % This means:
  %    R1R2'F = V*Tfull
  % only keep upper, non-zero block
  T = Tfull(1:r,:);
  I = speye(r,n);
  % and now:
  %    R1R2'F = V*I'*T
  %    R1R2' = V*I'*T*F'
  % Now we have that
  %     AE = Q*J*F*T'*I*V'
  % We want to solve:
  %     Ax = b
  %     AEE'x = b
  %     E'x = (Q*J*F*T'*I*V')⁺b
  %     x = E*(Q*J*F*T'*I*V')⁺b
  %     x = E*(V*(I'*(T'\(F'*(J'*(Q'*b))))))
  %
end

author: alecjacobson

repo: gptoolbox


text:
function tapas_physio_save_batch_mfile(matlabbatch, fileBatchM, ...
    doRemoveEmptyDefaults)
% Saves created matlabbatch variable to .m-file for later call of
% spm_jobman or run(fileBatchM)
%
%    tapas_physio_save_batch_mfile.m(matlabbatch, fileBatchM)
%
% IN
%   matlabbatch     variable from workspace
%   fileBatchM      m-file that we save to
%   doRemoveEmptyDefaults
%                   if true, all lines in the batch that have empty
%                   (default) values will be removed
%                   default: false
%
% OUT
%
% EXAMPLE
%   tapas_physio_save_batch_mfile.m
%
%   See also

% Author:   Lars Kasper
% Created:  2022-03-30
% Copyright (C) 2022 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under
% the terms of the GNU General Public License (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.

if nargin < 3
    doRemoveEmptyDefaults = false;
end

str = gencode(matlabbatch)';

if doRemoveEmptyDefaults
    indLineRemove = tapas_physio_find_string(str, ...
        {'= \[\];', '= {};', '= {''''};', '= '''';'});
    indLineRemove = cell2mat(indLineRemove);
    str(indLineRemove) = [];
end

nLines = numel(str);
fid = fopen(fileBatchM, 'w+');

for iLine = 1:nLines
    fprintf(fid, '%s\n', str{iLine});
end
fclose(fid);

author: translationalneuromodeling

repo: tapas


text:
%CODEGENERATOR.GENMEXFRICTION Generate C-MEX-function for joint friction
%
% cGen.genmexfriction() generates a robot-specific MEX-function to compute
% the vector of joint friction.
%
% Notes::
% - Is called by CodeGenerator.genfriction if cGen has active flag genmex
% - The MEX file uses the .c and .h files generated in the directory
%   specified by the ccodepath property of the CodeGenerator object.
% - Access to generated function is provided via subclass of SerialLink
%   whose class definition is stored in cGen.robjpath.
% - You will need a C compiler to use the generated MEX-functions. See the
%   MATLAB documentation on how to setup the compiler in MATLAB.
%   Nevertheless the basic C-MEX-code as such may be generated without a
%   compiler. In this case switch the cGen flag compilemex to false.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.gengravload.

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

function [] = genmexfriction(CGen)

%% Forward kinematics up to tool center point
CGen.logmsg([datestr(now),'\tGenerating friction vector MEX-function: ']);
symname = 'friction';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genmfungravload:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

funfilename = fullfile(CGen.robjpath,[symname,'.c']);
[QD] = CGen.rob.gencoords;

% Function description header
hStruct = createHeaderStruct(CGen.rob,symname);

% Generate and compile MEX function
CGen.mexfunction(tmpStruct.(symname),...
    'funfilename',funfilename,...
    'funname',[CGen.getrobfname,'_',symname],...
    'vars',{QD},...
    'output','F',...
    'header',hStruct);

CGen.logmsg('\t%s\n',' done!');

end

%% Definition of the function description header contents
function hStruct = createHeaderStruct(rob,fname)
[~,hStruct.funName] = fileparts(fname);
hStruct.calls = '';
hStruct.shortDescription = ['Joint friction for the ',rob.name,' arm.'];
hStruct.detailedDescription = {['Given a full set of generalized joint velocities the function'],...
    'computes the friction forces/torques. Angles have to be given in radians!'};
hStruct.inputs = { ['input1:  ',int2str(rob.n),'-element vector of generalized velocities.']};
hStruct.outputs = {['F:  [',int2str(rob.n),'x1] output vector of joint friction forces/torques.']};
hStruct.references = {'Robot Modeling and Control - Spong, Hutchinson, Vidyasagar',...
    'Modelling and Control of Robot Manipulators - Sciavicco, Siciliano',...
    'Introduction to Robotics, Mechanics and Control - Craig',...
    'Modeling, Identification & Control of Robots - Khalil & Dombre'};
hStruct.authors = {'This is an autogenerated function!',...
    'Code generator written by:',...
    'Joern Malzahn (joern.malzahn@tu-dortmund.de)'};
hStruct.seeAlso = {'gravload'};
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ x, error_norm, iter, flag ]  = sor ( A, x, b, w, max_it, tol )

%*****************************************************************************80
%
%% SOR solves the linear system Ax=b using the Successive Over-Relaxation Method.
%
%  Discussion:
%
%    When the parameter W
%    is set to 1, this is equivalent to the Gauss-Seidel iteration.
%
%  Reference:
%
%    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
%      June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
%      Charles Romine, Henk van der Vorst
%    Templates for the Solution of Linear Systems: Building Blocks for
%      Iterative Methods,
%    SIAM Publications, 1993.
%
%  Parameters:
%
%    Input, real A(N,N), the symmetric positive definite matrix.
%
%    Input, real X(N), the initial guess vector.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, real W, the relaxation scalar, between 0 and 2.
%
%    Input, integer MAX_IT, the maximum number of iterations.
%
%    Input, real TOL, an error tolerance.
%
%    Output, real X(N), the solution.
%
%    Output, real ERROR_NORM, the norm of the error.
%
%    Output, integer ITER, the number of iterations performed.
%
%    Output, integer FLAG, the return flag.
%    0 = the solution was found to within the specified tolerance.
%    1 = a satisfactory solution was not found.  The iteration limit
%        was exceeded.
%

%
%  Initialization.
%
  flag = 1;
  iter = 0;

  bnrm2 = norm ( b );
  if ( bnrm2 == 0.0 )
    bnrm2 = 1.0
  end

  r = b - A * x;
  error_norm = norm ( r ) / bnrm2;
  errorhist = [ ];
  errorhist(1) = error_norm;

  if ( error_norm < tol )
    flag = 0;
    return
  end
%
%  Split the matrix.
%
  [ M, N, b ] = split ( A, b, w, 2 );

  for iter = 1 : max_it

    x_1 = x;
%
%  Update the approximation.
%
    x = M \ ( N * x + b );
%
%  Compute the error.
%
    error_norm = norm ( x - x_1 ) / norm ( x );
    errorhist(iter+1) = error_norm;
%
%  Check for convergence.
%
    if ( error_norm <= tol )
      flag = 0;
      break
    end

  end
%
%  Restore the right hand side.
%
  b = b / w;

  error_norm = errorhist;

  return
end

author: johannesgerer

repo: jburkardt-m