text:
function pic_serial ( )

%*****************************************************************************80
%
%% PIC_SERIAL applies the Particle in Cell method to an electrostatic problem.
%
%  Discussion:
%
%    Simple Electrostatic Particle In Cell (PIC) code in MATLAB
%    Flow of solar wind around a charged plate
%
%    For more, visit http://www.particleincell.com/2010/es-pic-method/
%    and http://www.particleincell.com/2011/particle-in-cell-example/
%
%  Modified:
%
%    27 July 2013
%
%  Parameters:
%
%    Local, real A(NX*NY,NX*NY), the finite difference stencil matrix.
%
%    Local, real ACC, the acceleration experienced by a particle.
%
%    Local, real AMU, the atomic mass unit.
%
%    Local, integer BOX(2,2), the indices of the lower left and upper right
%    corners of the internal obstruction.
%
%    Local, real CHG(NX,NY), the charge distribution.
%
%    Local, real DEBYE, the Debye length.
%
%    Local, real DEN(NX,NY), the charge density.
%
%    Local, real DH, the "diameter" of a single cell.
%
%    Local, real DT, the time step.
%
%    Local, real E(2), the electric field experienced by a particle.
%
%    Local, real EFX(NX,NY), the X component of the electric field.
%
%    Local, real EFY(NX,NY), the Y component of the electric field.
%
%    Local, real EPS0, the permittivity of free space.
%
%    Local, real F, the Lorentz force experienced by a particle.
%
%    Local, real FLUX, the flux of entering particles.
%
%    Local, real K, the Boltzmann constant.
%    Oddly enough, this value seems never to be used.
%
%    Local, real M, the ion mass for molecular oxygen.
%
%    Local, real MP_Q, the macro-particle charge.
%
%    Local, integer N0, the average particle density per cubic meter.
%
%    Local, integer NN, the total number of nodes.
%
%    Local, integer NP, the number of particles.
%
%    Local, integer NP_INSERT, the number of (computational) particles inserted
%    per time step.
%
%    Local, real NPT, the number of real particles created per time step.
%
%    Local, integer NX, the number of nodes in the X direction.
%
%    Local, integer NY, the number of nodes in the Y direction.
%
%    Local, integer PART_MAX, the maximum number of particles.
%    Here, this is set to 20,000.
%
%    Local, real PART_V(NP,2), the X and Y velocities for each particle.
%
%    Local, real PART_X(NP,2), the X and Y coordinates for each particle.
%
%    Local, real PHI(NX,NY), the potential.
%
%    Local, real PHI_P, the wall potential.
%
%    Local, real PHI0, the reference potential.
%
%    Local, real QE, the elementary charge.
%
%    Local, integer SEED, a seed for the random number generator.
%
%    Local, real SPWT, the specific weight, real particles per macroparticle.
%
%    Local, integer STEP, the current time step.
%
%    Local, integer STEP_NUM, the number of time steps.
%
%    Local, real TE, the electron temperature in eV.
%
%    Local, real TI, the ion velocity in eV.
%
%    Local, real V_DRIFT, the ion injection velocity, 7 km/s.
%
%    Local, real VTH, the thermal velocity.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PIC_SERIAL\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Apply the Particle in Cell (PIC) method to\n' );
  fprintf ( 1, '  an electrostatic problem that models the flow of\n' )
  fprintf ( 1, '  the "solar wind" around an electrically charged plate.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  It is worth considering how to turn this serial program\n' );
  fprintf ( 1, '  into one that uses parallel programming techniques.\n' );
%
%  Initialize the random number generator.
%
  seed = 123456789;
% rng ( seed );
%
%  Set physical constants.
%
  eps0 = 8.854e-12;
  qe = 1.602e-19;
  k = 1.381e-23;
  amu = 1.661e-27;
  m = 32.0 * amu;
%
%  Set problem parameters.
%
  n0 = 1.0e12;
  phi0 = 0.0;
  te = 1.0;
  ti = 0.1;
  v_drift = 7000.0;
  phi_p = -5.0;
%
%  Calculate plasma parameters.
%
  debye = sqrt ( eps0 * te / ( n0 * qe ) );
  vth = sqrt ( 2.0 * qe * ti / m );
%
%  Set the simulation domain.
%
  nx = 16;
  ny = 10;
  nn = nx * ny;
  step_num = 200;
  dh = debye;
  np_insert = ( ny - 1 ) * 15;
%
%  Compute other quantities.
%
  dt = 0.1 * dh / v_drift;
  width = ( nx - 1 ) * dh;
  height = ( ny - 1 ) * dh;
%
%  Indices of the corners of the obstruction.
%
  ox1 = floor ( nx / 3 );
  ox2 = floor ( nx / 3 ) + 2;
  oy1 = 1;
  oy2 = floor ( ny / 2 );
%
%  Indices of the corners of the obstruction, in an array.
%
  box(1,1) = floor ( nx / 3 );
  box(1,2) = floor ( nx / 3 ) + 2;
  box(2,1) = 1;
  box(2,2) = floor ( ny / 2 );
%
%  Calculate the specific weight.
%
  flux = n0 * v_drift * height;
  npt = flux * dt;
  spwt = npt / np_insert;
  mp_q = 1.0;
  part_max = 20000;
%
%  Allocate particle arrays.
%
  part_x = zeros(part_max,2);
  part_v = zeros(part_max,2);
%
%  Set the finite difference stencil matrix.
%
  A = stencil ( nx, ny, dh, box );

  phi = ones ( nx, ny ) * phi0;

  np = 0;
%
%   Time loop
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Time  Particles\n' );
  fprintf ( 1, '\n' );

  for step = 1 : step_num
%
%  1. CALCULATE THE CHARGE DENSITY.
%
    chg = zeros(nx,ny);

    for p = 1 : np

      fi = 1.0 + part_x(p,1) / dh;
      i = floor ( fi );
      hx = fi - i;

      fj = 1.0 + part_x(p,2) / dh;
      j = floor ( fj );
      hy = fj - j;

      chg(i,  j)   = chg(i,  j)   + ( 1.0 - hx ) * ( 1.0 - hy );
      chg(i+1,j)   = chg(i+1,j)   +         hx   * ( 1.0 - hy );
      chg(i,  j+1) = chg(i,  j+1) + ( 1.0 - hx ) *         hy;
      chg(i+1,j+1) = chg(i+1,j+1) +         hx   *         hy;

    end
%
%  Calculate the density.
%  Along the boundaries, the calculated density must be doubled.
%  Add a density floor for plotting, and to help the solver.
%
    den(1:nx,1:ny) = spwt * mp_q * chg(1:nx,1:ny) / ( dh * dh );

    den(1,:)  = 2.0 * den(1,:);
    den(nx,:) = 2.0 * den(nx,:);
    den(:,1)  = 2.0 * den(:,1);
    den(:,ny) = 2.0 * den(:,ny);

    den = den + 10000.0;
%
%  2. CALCULATE THE POTENTIAL.
%
    phi = eval_2dpot_gs ( nx, ny, phi, A, den, n0, phi0, te, phi_p, box, ...
      eps0, qe );
%
%  3. CALCULATE THE ELECTRIC FIELD.
%  Use central differences at internal nodes,
%  forward or backward differences at boundaries.
%
    efx = zeros(nx,ny);
    efy = zeros(nx,ny);

    efx(2:nx-1,:) = ( phi(1:nx-2,:) - phi(3:nx,:) ) / ( 2.0 * dh );
    efy(:,2:ny-1) = ( phi(:,1:ny-2) - phi(:,3:ny) ) / ( 2.0 * dh );

    efx(1,:)  = ( phi(1,:)    - phi(2,:) )  / dh;
    efx(nx,:) = ( phi(nx-1,:) - phi(nx,:) ) / dh;
    efy(:,1)  = ( phi(:,1)    - phi(:,2) )  / dh;
    efy(:,ny) = ( phi(:,ny-1) - phi(:,ny) ) / dh;
%
%  4. GENERATE NEW PARTICLES
%
    if ( part_max - np <= np_insert  )
%     np_insert = part_max - np;
    end
%
%  The new particles have coordinates randomly chosen within the first X layer,
%  and any Y layer.
%
    part_x(np+1:np+np_insert,1) = rand ( np_insert, 1 ) * dh;
    part_x(np+1:np+np_insert,2) = rand ( np_insert, 1 ) * height;
%
%  Sample Maxwellian in x and y, add drift velocity in x.
%
    part_v(np+1:np+np_insert,1) = v_drift ...
      + ( - 1.5 + rand(np_insert,1) + rand(np_insert,1) ...
      + rand(np_insert,1) ) * vth;

    part_v(np+1:np+np_insert,2) = 0.5 ...
      * ( - 1.5 + rand(np_insert,1) + rand(np_insert,1) ...
      + rand(np_insert,1) ) * vth;

    np = np + np_insert;
%
%  5. MOVE ALL THE PARTICLES
%
    p = 1;

    while ( p <= np )

      fi = 1.0 + part_x(p) / dh;
      i = floor ( fi );
      hx = fi - i;

      fj = 1.0 + part_x(p,2) / dh;
      j = floor ( fj );
      hy = fj - j;

      e = [ 0.0, 0.0 ];
      e =     [ efx(i,j),     efy(i,j)     ] * ( 1.0 - hx ) * ( 1.0 - hy );
      e = e + [ efx(i+1,j),   efy(i+1,j)   ] *         hx   * ( 1.0 - hy );
      e = e + [ efx(i,j+1),   efy(i+1,j)   ] * ( 1.0 - hx ) *         hy;
      e = e + [ efx(i+1,j+1), efy(i+1,j+1) ] *         hx   *         hy;
%
%  Compute the Lorentz force and the corresponding acceleration.
%  Then update the particle velocity and position.
%
      f = qe * e;
      acc = f / m;
      part_v(p,:) = part_v(p,:) + acc * dt;
      part_x(p,:) = part_x(p,:) + part_v(p,:) * dt;
%
%  The bottom boundary is reflective.
%
      if ( part_x(p,2) < 0 )
        part_x(p,2) = - part_x(p,2);
        part_v(p,2) = - part_v(p,2);
      end
%
%  Is the particle inside the plate?
%
      in_box = ( box(1,1) <= i && i < box(1,2) && ...
                 box(2,1) <= j && j < box(2,2) );
%
%  Particle is absorbed if it passes left, right or top boundaries,
%  or is inside the plate.
%  Kill the particle by replacing it with the last particle.
%
      if ( part_x(p,1) < 0.0 || width <= part_x(p,1) || height <= part_x(p,2) || in_box )
        part_x(p,:) = part_x(np,:);
        part_v(p,:) = part_v(np,:);
        np = np - 1;
        p = p - 1;
      end

      p = p + 1;

    end
%
%  6. PLOT RESULTS
%
    if ( mod ( step, 25 ) == 0 || step == step_num )

      figure ( 1 );
      clf
      hold on
      contourf ( den' );
      colorbar;
      patch ( [ ox1 ox2 ox2 ox1 ], [ oy1, oy1, oy2, oy2 ], 'w' )
      title ( sprintf ( 'Density on step %i', step ), 'Fontsize', 16 );
      hold off;

      figure ( 2 );
      clf
      hold on
      contourf ( phi' );
      patch ( [ ox1 ox2 ox2 ox1 ], [ oy1, oy1, oy2, oy2 ], 'w' )
      title ( sprintf ( 'Potential on step %i', step ), 'Fontsize', 16 );
      colorbar;
      hold off

      drawnow;

    end

    fprintf ( 1, '  %6i  %9i\n', step, np );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PIC_SERIAL\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function phi = eval_2dpot_gs ( nx, ny, phi, A, den, n0, phi0, te, phi_p, ...
  box, eps0, qe )

%*****************************************************************************80
%
%% EVAL_2DPOT_GS computes the 2D potential using Gauss-Seidel iteration.
%
%  Discussion:
%
%    For more information, refer to
%      http://www.particleincell.com/2010/es-pic-method/
%    and
%      http://www.particleincell.com/2011/particle-in-cell-example/
%
%  Modified:
%
%    05 July 2013
%
%  Parameters:
%
%    Input, integer NX, NY, the number of nodes in the X and Y directions.
%
%    Input, real PHI(NX,NY), the current potential function.
%
%    Input, real A(NX*NY,NX*NY), the finite difference stencil matrix.
%
%    Input, real DEN(NX,NY), the charge density.
%
%    Input, real PHI0, the reference potential.
%
%    Input, real TE, the electron temperature in eV.
%
%    Input, real EPS0, the permittivity of free space.
%
%    Input, real QE, the elementary charge.
%
%    Output, real PHI(NX,NY), the updated potential function.
%
  nn = nx * ny;
%
%  Solver tolerance.
%
  tol = 0.1;
%
%  Convert the density and potential into column vectors.
%
  den = reshape ( den, numel(den), 1 );
  phi = reshape ( phi, numel(phi), 1 );
%
%  Carry out the Gauss-Seidel iteration.
%
  for it = 1 : 2000
%
%  Recalculate the right hand side, adding Boltzmann term for the electrons.
%
    b = den - n0 * exp ( ( phi - phi0 ) / te );
    b = - b * qe / eps0;
%
%  Set the boundaries.
%  Zero electric field on y = 0, y = L, x = L.
%  Fixed potential on X = 0.
%
    b(1:nx) = 0.0;
    b(nn-nx+1:nn) = 0.0;
    b(nx:nx:nn) = 0.0;
    b(1:nx:nn) = phi0;
%
%  Set the potential on the fixed nodes.
%
    for j = box(2,1) : box(2,2)
      b([box(1,1):box(1,2)]+(j-1)*nx) = ones ( box(1,2)-box(1,1)+1, 1 ) * phi_p;
    end
%
%  Apply the Gauss-Seidel update to the current solution estimate.
%
    for i = 1 : nn
      phi(i) = ( b(i) - A(i,1:i-1)  * phi(1:i-1) ...
                      - A(i,i+1:nn) * phi(i+1:nn) ) / A(i,i);
    end
%
%  Compute the residual.
%
    if ( mod ( it, 10 ) == 0 )
      res = norm ( b - A * phi );
      if ( res <= tol )
        phi = reshape ( phi, nx, ny );
        return;
      end
    end

  end
%
%  Check if the solver converged.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EVAL_2DPOT_GS - Warning\n' );
  fprintf ( 1, '  The Gauss-Seidel iteration did not converge.\n' );
  fprintf ( 1, '  Residual norm = %g\n', res );

  phi = reshape ( phi, nx, ny );

  return
end
function A = stencil ( nx, ny, dh, box )

%*****************************************************************************80
%
%% STENCIL sets up the finite difference stencil matrix.
%
%  Modified:
%
%    05 July 2013
%
%  Parameters:
%
%    Input, integer NX, NY, the number of nodes in the X and Y directions.
%
%    Local, real DH, the "diameter" of a single cell.
%
%    Input, integer BOX(2,2), contains the coordinates of the lower left
%    and upper right corners of the box.
%
%    Output, real A(NX*NY,NX*NY), the finite difference stencil matrix.
%
  A = zeros ( nx*ny, nx*ny );
%
%  For internal nodes, here are the node numberings for the
%  north, west, central, east, and south locations:
%
%             u-nx
%      u-1    u      u+1
%             u+nx
%
  for j = 2 : ny - 1
    for i = 2 : nx - 1

      u = ( j - 1 ) * nx + i;

      A(u,u)    = -4.0 / ( dh * dh );
      A(u,u-1)  =  1.0 / ( dh * dh );
      A(u,u+1)  =  1.0 / ( dh * dh );
      A(u,u-nx) =  1.0 / ( dh * dh );
      A(u,u+nx) =  1.0 / ( dh * dh );

    end
  end
%
%  Neumann boundary on y=0
%
  j = 1;
  for i = 1 : nx
    u = ( j - 1 ) * nx + i;
    A(u,u)    = -1.0 / dh;
    A(u,u+nx) =  1.0 / dh;
  end
%
%  Neumann boundary on y=height
%
  j = ny;
  for i = 1 : nx
    u = ( j - 1 ) * nx + i;
    A(u,u-nx) =  1.0 / dh;
    A(u,u)    = -1.0 / dh;
  end
%
%  Neumann boundary on x=width
%
  i = nx;
  for j = 1 : ny
    u = ( j - 1 ) * nx + i;
    A(u,:) = zeros ( 1, nx * ny );
    A(u,u-1) =  1.0 / dh;
    A(u,u)   = -1.0 / dh;
  end
%
%  Dirichlet boundary on x=0
%
  i = 1;
  for j = 1 : ny
    u = ( j - 1 ) * nx + i;
    A(u,:) = zeros ( 1, nx * ny );
    A(u,u) = 1.0;
  end
%
%  Dirichlet boundary on nodes corresponding to the plate
%
  for j = box(2,1) : box(2,2)
    for i = box(1,1) : box(1,2)
      u = ( j - 1 ) * nx + i;
      A(u,:) = zeros ( 1, nx * ny );
      A(u,u) = 1.0;
    end
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
%% Testing ground for quad mesh adaptive refinement
% This is the testing utilizing VEM's polytopal structure
% without using the tree-structure to maintain the mesh data.
%
% Outer loop: # of vertices
% inner loop:
%
% Step 1: find the center/centroid P of marked elements
% Remark: current version using naive way of locating centroid
%
% Step 2: quad-sect using P, find duplicate nodes
% Step 3: form new elements replacing the marked elements using cell
% Step 4: add vertices to the neighbor elements through edges
%
% Two sets of mesh are kept:
% 1. For refinement elem: the first 4 vertices are always the square clockwise,
% after those are the extra vertices for polygon.
% 2. For computation elemRefine: re-ordered vertices mesh.

clear; close all;

%%
[node,elemOrig] = squarequadmesh([0,1,0,1],0.25);
figure(1);
showmesh(node,elemOrig);
findnode(node);
set(gcf,'color','w','Position', [100 600 500 500])

% local ordering should always be kept as
% 2 -- 3
% |    |
% 1 -- 4

T = auxstructurequad(elemOrig);

elem = num2cell(elemOrig,2); % convert elem matrix to cell

elemVertexNumber = cellfun('length',elem);

N = size(node,1);

%%
% Nv = min(elemVertexNumber):max(elemVertexNumber)
Nv=4;

% find polygons with Nv vertices
idx = find(elemVertexNumber == Nv); % index of elements with Nv vertices
NT = length(idx); % the number of elements having the same number of vertices
% vertex index and coordinates
vertex = cell2mat(elem(idx));

x = reshape(node(vertex,1),NT,Nv);
y = reshape(node(vertex,2),NT,Nv);

xmax = max(x,[],2); % the rightmost vertices' x-coords for each element
ymax = max(y,[],2); % the uppermost vertices' y-coords for each element
xmin = min(x,[],2); % the leftmost vertices' x-coords for each element
ymin = min(y,[],2); % the lowermost vertices' y-coords for each element

xc = 0.5*(xmax + xmin);
yc = 0.5*(ymax + ymin);

x_nb = circshift(x,[0,-1]);
y_nb = circshift(y,[0,-1]);

%% marking elements
idxMarkedElem = idx(randi(NT, [floor(NT/3), 1]));
% idxMarkedElem = idx([1,3]);
NmarkedElem = length(idxMarkedElem);

NnodeNew = 5*NmarkedElem;
% number of the new nodes per Nv consisting duplicated nodes

% o -- 2 -- o
% |    |    |
% 1 -- 5 -- 3
% |    |    |
% o -- 4 -- o
% o: original vertices
% 1 to 5: local indices of newly added vertices
% the order is consistent with neighbor ordering

nodeNew(1:5:NnodeNew,:) = [xmin(idxMarkedElem,1), yc(idxMarkedElem,1)]; % left
nodeNew(2:5:NnodeNew,:) = [xc(idxMarkedElem,1), ymax(idxMarkedElem,1)]; % top
nodeNew(3:5:NnodeNew,:) = [xmax(idxMarkedElem,1), yc(idxMarkedElem,1)]; % right
nodeNew(4:5:NnodeNew,:) = [xc(idxMarkedElem,1), ymin(idxMarkedElem,1)]; % bottom
nodeNew(5:5:NnodeNew,:) = [xc(idxMarkedElem,1), yc(idxMarkedElem,1)]; % center

[nodeNewU, idxN, idxU] = unique(nodeNew,'rows','legacy');
% nodeNewU = nodeNew(idxN,:) and nodeNew = nodeNewU(idxU,:)

idxNodeNew = N+1:N+NnodeNew; % unique indices of new nodes

%% generate new elements on marked elements
elem2nodeNew = reshape(idxNodeNew(idxU), [5, NmarkedElem])';
% elem2nodeNew(i,:): i-th marked elem's newly added vertices
% elem2nodeNew(:,1): center of the marked elements
% elem2nodeNew(:,2:5): counterclockwisely ordered new edge vertices from
% the left edge locally

node = [node; nodeNewU];

elemMarked = cell2mat(elem(idxMarkedElem));

% localNewVertex = [2 1 3; 5 1 2; 4 1 5; 3 1 4]; %not very helpful

elemNew1 = horzcat(elemMarked(:,1), elem2nodeNew(:,[1 5 4]));
elemNew2 = horzcat(elem2nodeNew(:,1), elemMarked(:,2), elem2nodeNew(:,[2 5]));
elemNew3 = horzcat(elem2nodeNew(:,[5 2]), elemMarked(:,3), elem2nodeNew(:,3));
elemNew4 = horzcat(elem2nodeNew(:,[4 5 3]), elemMarked(:,4));

NelemNew = 4*NmarkedElem;

%% change elements' vertices on unrefined neighbor elements
% no need to check whether a vertex is added twice b/c
% the marked elements will get over-written in the next routine

elemMarkedNB = T.neighbor(idxMarkedElem,:);
% idxNb = ~(elemMarkedNB==idxMarkedElem);
elemMarkedNB = num2cell(elemMarkedNB,2);

for i = 1:size(elemMarked,1)
    for j = 1:size(elemMarkedNB{i},2)
        elem{elemMarkedNB{i}(:,j)} = ...
         [elem{elemMarkedNB{i}(:,j)}, elem2nodeNew(i,j)];
    end
end

%% vertically concat all new elements into a cell
elemNew = reshape([elemNew1'; elemNew2'; elemNew3'; elemNew4'], [4, NelemNew])';

elem(idxMarkedElem) = mat2cell(elemNew, repmat(4,[1,NmarkedElem]));

% elemTemp = vertcat(elem{:});
% above doesn't work if elem is of different size when nested

outTemp = cellfun(@(x) num2cell(cat(2,x),2),elem,'UniformOutput',false);
elemRefine = cellflat(outTemp)';

% celldisp(elemRefine)

%% reorder the polygon with vertices > 4
% to-do: rewrite the for loop to cellfun
elemVertexNumberRefine = cellfun('length',elemRefine);
idxReorder = find(elemVertexNumberRefine > 4);

for i = 1:length(idxReorder)
    % compute the center (simple mean suffices) and sorting atan2
    elemTemp = elemRefine{idxReorder(i)};
    center = mean(node(elemTemp,:),1);
    vert2center = node(elemTemp,:) - center;
    [~, idxOrdered] = sort(atan2(vert2center(:,1), vert2center(:,2)));
    elemRefine{idxReorder(i)} = elemTemp(idxOrdered);
end

Nrefine = size(node,1); % update number of nodes

%%
figure(2);
showmeshpoly(node,elemRefine);
plot(xc,yc,'o');
findnode(node);
set(gcf,'color','w','Position', [650 600 500 500])

repo: ifem


text:
function batch_register_scans_with_struct_priors(baseScanFileNameC,...
    movScansFileNameC, basePriorStructNameC, movPriorStructNameC, ...
    movStrsToDeformC, movDoseToDeformC, registeredDir, paired_scans_flag)
% function batch_register_scans_with_struct_priors(baseScanFileNameC,...
%     movScansFileNameC, basePriorStructNameC, movPriorStructNameC, ...
%     strsToDeformC, registeredDir, paired_scans_flag)
%
% This script registers movScansFileNameC to baseScanFileNameC using
% basePriorStructNameC and movPriorStructNameC structure priors for initial alignment.
% Additionally, the structurnames specified in strsToDeformC are warped on to the base
% scan and the deformed images and segmwntation are written to registeredDir.
%
% APA, 12/30/2021

% RIDER template
% baseScan = '/lab/deasylab1/Data/RTOG0617/CERR_files_tcia/rider_template/RIDER-1225316081_First_resampled_1x1x3.mat';

% Directory containing CERR files
%movScanDir = '/lab/deasylab1/Data/RTOG0617/CERR_files_tcia/dose_mapping_original_plans'; % RTOG0617
%movScanDir = '/lab/deasylab1/Data/Lung_PORT/cerr_files'; % PORT
%movScanDir = '/lab/deasylab1/Maria/Heart_pCT_NSCLC/CERR/All_ph_1_pts_GTVcorrect/ATRIA_VENTR_PERI_DONE/N_241/0617structs_QAed_MT'; % OLD Conventional
% movScanDir = '/lab/deasylab1/Data/LA-NSCLC_Durva_N113/cerr_files_planning_only'; % NEW Conventional
%movScanDir = '/lab/deasylab1/Maria/Lung_EarlyStage_MSK_PFS_pCT_CT_PETCT/CERR_DCM_OMT_Validation/cerr_files_all_fx_corrected'; % MSK early stage Lung (fractionation corrected)

% Directory to write registered scans
%registeredDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/CERR_files_RTOG0617_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
%registeredDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/CERR_files_PORT_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
%registeredDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/CERR_files_OLD_CONVN_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
% registeredDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/CERR_files_NEW_CONVN_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
%registeredDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/CERR_files_MSK_EarlyStage_to_RIDER_1225316081_First_template'; % MSK early stage Lung (fractionation corrected)

if ~exist(registeredDir,'dir')
    mkdir(registeredDir)
end

% Filename to store dice scores
%diceSaveFileName = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/RTOG0617_peri_dice_RIDER.mat'; %RTOG0617
%diceSaveFileName = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/PORT_peri_dice_RIDER.mat'; %PORT
%diceSaveFileName ='/lab/deasylab1/Data/RTOG0617/registrations_pericardium/OLD_CONVN_peri_dice_RIDER.mat'; %OLD_CONVN
% diceSaveFileName ='/lab/deasylab1/Data/RTOG0617/registrations_pericardium/NEW_CONVN_peri_dice_RIDER.mat'; %NEW_CONVN
%diceSaveFileName = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/MSK_EarlyStage_peri_dice_RIDER.mat'; %EarlyStage
%
% strsToDeformC = {'DL_HEART_MT','DL_HEART','DL_AORTA', 'DL_LA', 'DL_LV', 'DL_RA', 'DL_RV', 'DL_IVC',...
%     'DL_SVC', 'DL_PA', 'DL_ATRIA', 'DL_PERICARDIUM', 'DL_VENTRICLES'};
% strsToDeformC = {'DL_HEART_MT','DL_PERICARDIUM'};
% strsToDeformC = {'DL_HEART','DL_PERICARDIUM'};

% movDirS = dir(movScanDir);
% movDirS([movDirS.isdir]) = [];
% movScanC = fullfile(movScanDir,{movDirS.name});

% registerToAtlasMultipleScans(baseScan,movScanC,registeredDir,...
%     strNameToWarp,initPlmCmdFile,refinePlmCmdFile)

%hpool = parpool(12);

% % Load base planC
% planC = loadPlanC(baseScan,tempdir);
% planC = updatePlanFields(planC);
% planC = quality_assure_planC(baseScan,planC);
% indexS = planC{end};
%
% % Generate mask by excluding GTV
% indTumor = findTumorIndex(1,planC);
% baseMask3M = [];
% if length(indTumor) == 1
%     baseMask3M = ~getStrMask(indTumor,planC);
% end
% strC = {planC{indexS.structures}.structureName};
% pericardInd = getMatchingIndex('DL_PERICARDIUM',strC,'exact');
% basePeriMask3M = getStrMask(pericardInd,planC);
% % [minr, maxr, minc, maxc, mins, maxs,baseMask3M] = compute_boundingbox(baseMask3M,3);
% se = strel('sphere',5);
% basePeriMaskDilate3M = imdilate(basePeriMask3M, se);
% diffPeriMask3M = xor(basePeriMaskDilate3M,basePeriMask3M);
% planC{indexS.scan}(1).scanArray(diffPeriMask3M) = 5000;

for iBase = 1:length(baseScanFileNameC)

    if paired_scans_flag
        movScansFileToRegisterC = movScansFileNameC(iBase);
        movStrsToDeformToUseC = movStrsToDeformC(iBase);
        movPriorStructNameToUseC = movPriorStructNameC(iBase);
        movDoseToDeformToUseC = movDoseToDeformC(iBase);
    else
        movScansFileToRegisterC = movScansFileNameC;
        movStrsToDeformToUseC = movStrsToDeformC;
        movPriorStructNameToUseC = movStructNameC;
        movDoseToDeformToUseC = movDoseToDeformC;
    end

    baseScan = baseScanFileNameC{iBase};

    % Load base planC
    planC = loadPlanC(baseScan,tempdir);
    planC = updatePlanFields(planC);
    planC = quality_assure_planC(baseScan,planC);
    indexS = planC{end};

    % Mask to restrict cost function
    baseMask3M = [];

    baseLandmarkStrC = basePriorStructNameC{iBase}; %{'DL_PERICARDIUM'};
    baseScanNum = 1;
    strBaseC = {planC{indexS.structures}.structureName};
    baseLandmarkListM = [];
    xyzBaseM = zeros(length(baseLandmarkStrC),3);
    for iStr = 1:length(baseLandmarkStrC)
        baseLandmarkStr = baseLandmarkStrC{iStr};
        strIndV = getMatchingIndex(baseLandmarkStr,strBaseC,'exact');
        assocScanNumV = getStructureAssociatedScan(strIndV,planC);
        baseLandmarkInd = strIndV(assocScanNumV == baseScanNum);
        [baseX,baseY,baseZ] = calcIsocenter(baseLandmarkInd,'COM',planC);
        xyzBaseV = [baseX,baseY,baseZ];
        xyzBaseM(iStr,:) = xyzBaseV;
        [xBaseV,yBaseV,zBaseV] = getScanXYZVals(planC{indexS.scan}(baseScanNum));
        baseY =  -yBaseV(1)+(yBaseV(1)-baseY);
        baseZ = -zBaseV(end) + zBaseV(end) - baseZ;
        baseLandmarkListM(iStr,:,1) = [iStr-1, [-baseX,-baseY,baseZ]*10];
    end

    allPlanC = {planC};

    for movNum = 1:length(movScansFileToRegisterC)

        movScanFileName = movScansFileToRegisterC{movNum};

        %[~,fName] = fileparts(movScanFileName);
        %newFileName = fullfile(registeredDir,[fName,'.mat']);
        %if exist(newFileName,'file')
        %    continue;
        %end

        try

            %     baseStrNumV = [];
            %     movStrNumV = [];

            %     % Load base planC
            %     planC = loadPlanC(baseScan,tempdir);
            %     planC = updatePlanFields(planC);
            %     planC = quality_assure_planC(baseScan,planC);
            %     indexS = planC{end};
            %
            %     % Generate mask by excluding GTV
            %     indTumor = findTumorIndex(1,planC);
            %     if length(indTumor) == 1
            %         baseMask3M = ~getStrMask(indTumor,planC);
            %     end

            % Load moving planC as planD
            planD = loadPlanC(movScanFileName,tempdir);
            planD = updatePlanFields(planD);
            planD = quality_assure_planC(movScanFileName,planD);
            indexSD = planD{end};

            % Register all image representations to each other starting from the
            % transformation generated previously
            %numScans = length(planC{indexS.scan});
            %     numScans = 1;
            %     allStrNumV = 1:length(planD{indexSD.structures}); % -1 to omit the "noise" structure
            %     assocScanV = getStructureAssociatedScan(allStrNumV,planD);
            %
            %     outBspFile = '';

            %inputCmdFile = '';
            %     numScans = 1;
            scanNumV = 1;
            for iScan = 1:length(scanNumV)

                scanNum = scanNumV(iScan);

                % Registration settings
                registration_tool = 'PLASTIMATCH';
                plmSettingsFile = '/lab/deasylab1/Data/RTOG0617/plastimatch_registration_settings/plm_reg_settings_with_landmark_stiffness_pericard.txt';
                %plmSettingsFile = '//vpensmph/deasylab1/Data/RTOG0617/plastimatch_registration_settings/plm_reg_settings_with_landmark_stiffness_pericard.txt';
                %plmSettingsFile = 'L:/Data/RTOG0617/plastimatch_registration_settings/plm_reg_settings_with_landmark_stiffness_pericard.txt';
                %baseMask3M = [];
                movMask3M = [];
                threshold_bone = -800;
                inBspFile = '';
                outBspFile = '';
                %tmpDirPath = '/lab/deasylab1/Data/RTOG0617/cerr_registration_temp_dir';
                tmpDirPath = '/cluster/home/aptea/segSessions/registration_temp';
                %tmpDirPath = '//vpensmph/deasylab1/Data/RTOG0617/cerr_registration_temp_dir';
                %tmpDirPath = 'L:/Data/RTOG0617/cerr_registration_temp_dir';
                algorithm = 'BSPLINE PLASTIMATCH';

                %         % Generate mask by excluding GTV
                %         indTumor = findTumorIndex(scanNum,planD);
                %         if length(indTumor) == 1
                %             movMask3M = ~getStrMask(indTumor,planD);
                %         end
                %         strC = {planD{indexSD.structures}.structureName};
                %         pericardInd = getMatchingIndex('DL_PERICARDIUM',strC,'exact');
                %         movPeriMask3M = getStrMask(pericardInd,planD);
                %         %         [minr, maxr, minc, maxc, mins, maxs,movMask3M] = compute_boundingbox(movMask3M,3);
                %         movPeriMaskDilate3M = imdilate(movPeriMask3M, se);
                %         diffPeriMask3M = xor(movPeriMaskDilate3M,movPeriMask3M);
                %         planD{indexSD.scan}(1).scanArray(diffPeriMask3M) = 5000;

                % Generate Landmarks
                movLandmarkStrC = movPriorStructNameToUseC{movNum}; %{'DL_PERICARDIUM'};

                % Register planD to planC
                baseScanNum = scanNum;
                movScanNum  = scanNum;
                %strBaseC = {planC{indexS.structures}.structureName};
                strMovC = {planD{indexSD.structures}.structureName};
                landmarkListM = [];
                xyzMovM = nan(length(movLandmarkStrC),3);
                for iStr = 1:length(movLandmarkStrC)
                    movLandmarkStr = movLandmarkStrC{iStr};
                    strIndV = getMatchingIndex(movLandmarkStr,strMovC,'exact');
                    assocScanNumV = getStructureAssociatedScan(strIndV,planD);
                    movLandmarkInd = strIndV(assocScanNumV == movScanNum);
                    if length(movLandmarkInd) ~= 1
                        continue;
                    end
                    [movX,movY,movZ] = calcIsocenter(movLandmarkInd,'COM',planD);
                    xyzMovV = [movX,movY,movZ];
                    xyzMovM(iStr,:) = xyzMovV;
                    [xMovV,yMovV,zMovV] = getScanXYZVals(planD{indexSD.scan}(movScanNum));
                    movY =  -yMovV(1)+(yMovV(1)-movY);
                    movZ = -zMovV(end) + zMovV(end) - movZ;
                    landmarkListM(iStr,:,1) = baseLandmarkListM(iStr,:,1);
                    landmarkListM(iStr,:,2) = [iStr-1, [-movX,-movY,movZ]*10];
                end

                % Don't use landmarks (update it to use landmarks in future)
                landmarkListM = [];

                %         initialTranslationXyzV = xyzMovV-xyzBaseV;
                initialTranslationXyzM = xyzMovM - xyzBaseM;
                initialTranslationXyzV = nanmean(initialTranslationXyzM,1);
                initialTranslationXyzV(1) = initialTranslationXyzV(1);
                initialTranslationXyzV(2) = -initialTranslationXyzV(2);
                initialTranslationXyzV(3) = -initialTranslationXyzV(3);

                plnC = allPlanC{1};

                % Call registration wrapper
                plnC = register_scans(plnC, baseScanNum, planD, movScanNum,...
                    algorithm, registration_tool, tmpDirPath, baseMask3M,...
                    movMask3M, threshold_bone, plmSettingsFile, inBspFile,...
                    outBspFile, landmarkListM,initialTranslationXyzV);

                strCreationScanNum = scanNum;

                % Warp scan (already warped in register_scans.m)
                deformS = plnC{indexS.deform}(end);
                %planC = warp_scan(deformS,movScanNum,planD,planC);

                strsToDeformC = movStrsToDeformToUseC{movScanNum};

                movStructNumsV = [];
                strNamC = {planD{indexSD.structures}.structureName};
                for iStr = 1:length(strsToDeformC)
                    movStructNumsV = [movStructNumsV,...
                        getMatchingIndex(strsToDeformC{iStr},strNamC,'exact')];
                end

                %planC = warp_structures(deformS,strCreationScanNum,movStructNumsV,planD,planC);
                plnC = warp_structures(deformS,strCreationScanNum,movStructNumsV,planD,plnC,tmpDirPath);

                doseCreationScanNum = strCreationScanNum;
%                 numDoses = length(planD{indexSD.dose});
                doseToDeformC = movDoseToDeformToUseC{movScanNum}; % update to get the correct dose num
                if ~isempty(doseToDeformC)
                    doseNum = length(planD{indexSD.dose});
                    %for doseNum = 1:numDoses
                    %planC = warp_dose(deformS,doseCreationScanNum,doseNum,planD,planC);
                    plnC = warp_dose(deformS,doseCreationScanNum,doseNum,planD,plnC,tmpDirPath);
                    %end
                end

                % delete the output nrrd file
                if strcmpi(algorithm,'DEMONS PLASTIMATCH')
                    delete(outBspFile)
                end

            end

            if exist(inBspFile,'file')
                delete(inBspFile)
            end

            % Save base and moving scans
            [~,fName] = fileparts(movScanFileName);
            newFileName = fullfile(registeredDir,fName);
            save_planC(plnC,[],'passed',newFileName);

        catch

        end % end try-catch

    end  % Iterate over moving scans

end % Iterate over base scans

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function newval=meshremap(fromval,elemid,elembary,toelem,nodeto)
%
% newval=meshremap(fromval,elemid,elembary,toelem,nodeto)
%
% Redistribute nodal values from the source mesh to the target mesh so that
% the sum of each property on each mesh is the same
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%	 fromval: values defined at the source mesh nodes, the row or column
%	          number must be the same as the source mesh node number, which
%	          is the same as the elemid length
%	 elemid: the IDs of the target mesh element that encloses the nodes of
%            the source mesh nodes; a vector of length of src mesh node
%            count; elemid and elembary can be generated by calling
%
%           [elemid,elembary]=tsearchn(node_target, elem_target, node_src);
%
%           note that the mapping here is inverse to that in meshinterp()
%
%	 elembary: the bary-centric coordinates of each source mesh nodes
%	         within the target mesh elements, sum of each row is 1, expect
%	         3 or 4 columns (or can be N-D)
%    toelem: the element list of the target mesh
%    nodeto: the total number of target mesh nodes
%
%
% output:
%	 newval: a 2D array with rows equal to the target mesh nodes (nodeto),
%            and columns equals to the value numbers defined at each source
%            mesh node
% example:
%
%    [n1,f1,e1]=meshabox([0 0 0],[10 20 5],1); % src mesh
%    [n2,f2,e2]=meshabox([0 0 0],[10 20 5],2); % target mesh
%    [id, ww]=tsearchn(n2,e2,n1);              % project src to target mesh
%    value_src=n1(:,[2 3 1]);             % create dummy values at src mesh
%    newval=meshremap(value_src,id,ww,e2,size(n2,1)); % map to target
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(fromval,1)==1)
    fromval=fromval(:);
end

if(size(fromval,2)==length(elemid))
    fromval=fromval.';
end

newval=zeros(nodeto,size(fromval,2));

idx=~isnan(elemid);
fromval=fromval(idx,:);
elembary=elembary(idx,:);
idx=elemid(idx);

nodeval=repmat(fromval,1,1,size(elembary,2)).*repmat(permute(elembary,[1,3,2]),1,size(fromval,2),1);

for i=1:size(elembary,2)
    [ix,iy]=meshgrid(toelem(idx,i),1:size(fromval,2));
    nval=nodeval(:,:,i).';
    newval=newval + accumarray([ix(:),iy(:)],nval(:), size(newval));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ t_stop, y_stop ] = p33_stop ( neqn )

%*****************************************************************************80
%
%% P33_STOP returns the stopping point for problem p33.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
    0.756245; ...
    2.67294 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function  [oldlb,oldub] = presolve_probebinary(A,b,c,lb,ub,binary_variables);

goon = 1;

AL0A  = (A>0).*A;
AG0A  = (A<0).*A;

At = A';
%used = full(any(A(:,find(changed_bounds)),2));
isbinary  = ismembc(1:length(lb),binary_variables);
isinteger = ismembc(1:length(lb),binary_variables);

goon = all(lb<=ub);

oldub = ub;
oldlb = lb;

for i = 1:length(binary_variables)

    % Probe up
    lb = oldlb;
    ub = oldub;
    if lb(binary_variables(i))~=ub(binary_variables(i))
        lb(binary_variables(i)) = 1;
        bi_dn = AL0A*lb+AG0A*ub;
        if any(bi_dn>b)
            oldlb(binary_variables(i)) = 0;
            oldub(binary_variables(i)) = 0;
        end
    end
    % Probe down
    lb = oldlb;
    ub = oldub;
    if lb(binary_variables(i))~=ub(binary_variables(i))
        ub(binary_variables(i)) = 0;
        bi_dn = AL0A*lb+AG0A*ub;
        if any(bi_dn>b)
            oldlb(binary_variables(i)) = 1;
            oldub(binary_variables(i)) = 1;
        end
    end
end

fix_this = [];
fix_value = [];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function geometry_test0201 ( )

%*****************************************************************************80
%
%% TEST0201 tests CYLINDER_POINT_DIST_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  test_num = 6;

  dist_test = [ 3.0, 0.5, 5.0, 8.0, 1.0, 0.25 ];
  p_test = [ ...
      4.0,   0.5,  0.0; ...
     -0.5,  -1.0,  0.0; ...
      4.0,   6.0,  0.0; ...
      0.75, -10.0, 0.0; ...
      0.0,   0.0,  0.0; ...
      0.25,  1.75, 0.0 ]';
  p1 = [ 0.0, -2.0, 0.0 ];
  p2 = [ 0.0,  2.0, 0.0 ];
  r = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0201\n' );
  fprintf ( 1, '  CYLINDER_POINT_DIST_3D computes the distance\n' );
  fprintf ( 1, '  to a cylinder.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radius R = %f\n', r );
  fprintf ( 1, '  Center of bottom disk = %10f  %10f  %10f\n', p1(1:dim_num) );
  fprintf ( 1, '  Center of top disk =    %10f  %10f  %10f\n', p2(1:dim_num) );

  for test = 1 : test_num

    p(1:dim_num) = p_test(1:dim_num,test);

    fprintf ( 1, '\n' );
    fprintf ( 1, '  P =    %10f  %10f  %10f\n', p(1:dim_num) );

    dist = cylinder_point_dist_3d ( p1, p2, r, p );

    fprintf ( 1, '  DIST (computed) = %f\n', dist );
    fprintf ( 1, '  DIST (exact) =    %f\n', dist_test(test) );

  end

  return
end


text:
classdef MPSOD < ALGORITHM
% <multi/many> <real/integer>
% Multi-objective particle swarm optimization algorithm based on
% decomposition

%------------------------------- Reference --------------------------------
% C. Dai, Y. Wang, and M. Ye, A new multi-objective particle swarm
% optimization algorithm based on decomposition, Information Sciences,
% 2015, 325: 541-557.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate the weight vectors
            [W,Problem.N] = UniformPoint(Problem.N,Problem.M);
            W = W./repmat(sqrt(sum(W.^2,2)),1,size(W,2));
            T = ceil(Problem.N/10);

            %% Detect the neighbours of each solution
            B = pdist2(W,W);
            [~,B] = sort(B,2);
            B = B(:,1:T);

            %% Generate random population
            Population = Problem.Initialization(2*Problem.N);
            Z          = min(Population.objs,[],1);
            Population = Classification(Problem,Population,W,Z);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                [Parent,Pbest,Gbest] = MatingSelection(Population.objs,B,W,Z);
                Offspring  = Operator(Problem,Population(Parent),Population(Pbest),Population(Gbest));
                Z          = min([Z;Offspring.objs],[],1);
                Population = Classification(Problem,[Population,Offspring],W,Z);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [beta,r,J,Sigma,mse,errorparam,robustw] = nlinmultifit(x_cell, y_cell, mdl_cell, beta0, options)
%NLINMULTIFIT Nonlinear least-squares regression of multiple data sets
%
%	A wrapper function for NLINFIT which allows simulatenous fitting of
%	multiple data sets with shared fitting parameters. See example below.
%
%	Unlike different solutions (using fminsearch or similar functions)
%	this approach enables simple estimation of model predictions and
%	their confidence intervals, as well as confidence intervals on the
%	fitted parameters (using the built-in NLPREDCI and NLPARCI functions).
%
%	KNOWN ISSUE:
%		In this implementation, different data sets are weighted according
%		to their relative lengths. Special care should be taken when those
%		lengths are considerably different from one another to avoid
%		biased results towards one data set in particular.
%
%	INPUT:
% 		x_cell,y_cell: Cell arrays containing the x,y vectors of the fitted
%					   data sets.
% 		mdl_cell: Cell array containing model functions for each data set.
% 		beta0: Vector containing initial guess of the fitted parameters.
% 		options: Structure containing control parameters for NLINFIT (see
%				 help file on NLINFIT for more details).
%
%	OUTPUT:
%		beta,r,J,Sigma,mse,errorparam,robustw: Direct output from NLINFIT.
%
%	EXAMPLE:
% 		% Generate X vectors for both data sets
% 		x1 = 0:0.1:10;
% 		x2 = 0:0.2:10;
%
% 		% Generate Y data with some noise
% 		y1 = cos(2*pi*0.5*x1).*exp(-x1/5) + 0.1*randn(size(x1));
% 		y2 = 0.5 + 2*exp(-(x2/5)) + 0.1*randn(size(x2));
%
% 		% Define fitting functions and parameters, with identical
%		% exponential decay for both data sets
% 		mdl1 = @(beta,x) cos(2*pi*beta(1)*x).*exp(-x/beta(2));
% 		mdl2 = @(beta,x) beta(4) + beta(3)*exp(-(x/beta(2)));
%
% 		% Prepare input for NLINMULTIFIT and perform fitting
% 		x_cell = {x1, x2};
% 		y_cell = {y1, y2};
% 		mdl_cell = {mdl1, mdl2};
% 		beta0 = [1, 1, 1, 1];
% 		[beta,r,J,Sigma,mse,errorparam,robustw] = ...
%					nlinmultifit(x_cell, y_cell, mdl_cell, beta0);
%
% 		% Calculate model predictions and confidence intervals
% 		[ypred1,delta1] = nlpredci(mdl1,x1,beta,r,'covar',Sigma);
% 		[ypred2,delta2] = nlpredci(mdl2,x2,beta,r,'covar',Sigma);
%
% 		% Calculate parameter confidence intervals
% 		ci = nlparci(beta,r,'Jacobian',J);
%
% 		% Plot results
% 		figure;
% 		hold all;
% 		box on;
% 		scatter(x1,y1);
% 		scatter(x2,y2);
% 		plot(x1,ypred1,'Color','blue');
% 		plot(x1,ypred1+delta1','Color','blue','LineStyle',':');
% 		plot(x1,ypred1-delta1','Color','blue','LineStyle',':');
% 		plot(x2,ypred2,'Color',[0 0.5 0]);
% 		plot(x2,ypred2+delta2','Color',[0 0.5 0],'LineStyle',':');
% 		plot(x2,ypred2-delta2','Color',[0 0.5 0],'LineStyle',':');
%
%	AUTHOR:
%		Chen Avinadav
%		mygiga (at) gmail
%

	num_curves = length(x_cell);
	if length(y_cell) ~= num_curves || length(mdl_cell) ~= num_curves
		error('Invalid input to NLINMULTIFIT');
	end

	x_vec = [];
	y_vec = [];
	mdl_vec = '@(beta,x) [';
	mdl_ind1 = 1;
	mdl_ind2 = 0;
	for ii = 1:num_curves
		if length(x_cell{ii}) ~= length(y_cell{ii})
			error('Invalid input to NLINMULTIFIT');
		end
		if size(x_cell{ii},2) == 1
			x_cell{ii} = x_cell{ii}';
		end
		if size(y_cell{ii},2) == 1
			y_cell{ii} = y_cell{ii}';
		end
		x_vec = [x_vec, x_cell{ii}];
		y_vec = [y_vec, y_cell{ii}];
		mdl_ind2 = mdl_ind2 + length(x_cell{ii});
		mdl_vec = [mdl_vec, sprintf('mdl_cell{%d}(beta,x(%d:%d)), ', ii, mdl_ind1, mdl_ind2)];
		mdl_ind1 = mdl_ind1 + length(x_cell{ii});
	end
	mdl_vec = [mdl_vec(1:end-2), '];'];
	mdl_vec = eval(mdl_vec);

	if nargin == 4
		[beta,r,J,Sigma,mse,errorparam,robustw] = nlinfit(x_vec, y_vec, mdl_vec, beta0);
	else
		[beta,r,J,Sigma,mse,errorparam,robustw] = nlinfit(x_vec, y_vec, mdl_vec, beta0, options);
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function s = rectangular_plate_analytical_solution (a, b, h, E, nu, q0)
%
% INPUT
%   a:  length of the rectangular section
%   b:  height of the rectangular section
%   h:  thickness of the plate
%   E:  Young modulus
%   nu: Poisson ratio
%   q0: distributed load
%
% OUTPUT
%   displ: maximum displacement

x=a/2;
y=b/2;
%pi = %pi;
D=E*h^3/(12*(1-nu^2));
ks = 16*q0*a^4*b^4/(pi^6*D);
km = 16*q0*a^4*b^4/(pi^4);

s=0;
mx=0;
my=0;
n = 1;
nMax = 1001;
while n<=nMax
    m = 1;
    while m<=nMax
        den = (n^2*a^2 + m^2*b^2)^2;
        sinA = sin(m*pi*x/a);
        sinB = sin(n*pi*y/b);
        knm = sinA*sinB/(m*n*den);
        s = s + ks*knm;
        mx = mx + km*knm*((m/a)^2 + nu*(n/b)^2);
        my = my + km*knm*(nu*(m/a)^2 + (n/b)^2);
        m = m+2;
    end;
    n = n+2;
end

x=0;
y=0;
kmxy = -(1-nu)*16*q0*a^3*b^3/(pi^4);
mxy = 0;
n=1;
while n<=nMax
    m = 1;
    while m<=nMax
        den = (n^2*a^2 + m^2*b^2)^2;
        cosA = cos(m*pi*x/a);
        cosB = cos(n*pi*y/b);
        mxy = mxy + kmxy*cosA*cosB/den;
        m = m+2;
    end;
    n = n+2;
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%% SCENARIO BUILDER CLASS (scenarioBuilder.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The scenario builder is a tool designed to provide a set of tools to aid
% in the design of object scenarios in the global scene.

classdef scenarioBuilder
    % BASIC SCENARIO DEFINITION
    % This class handles the methods of generating different scenarios of
    % a defined number of objects/agents.

    % SCENARIO PROPERTIES
    properties
        % META DATA
        objects;
        name = 'Default Scenario.';
        % GLOBAL DATA
        % ALL GLOBAL DATA IS GIVEN IN THE EAST-NORTH-UP (ENU) FRAME OF
        % REFERENCE
        globalTriad = [1,0,0;0,1,0;0,0,1];      % XYZ (Earth) frame of reference
        positions;                              % Global positions
        velocities;                             % Global velocities
        quaternions;                            % The quaternion moving from the global XYZ
    end
    % PUBLIC METHODS
    methods
        % CONSTRUCTION METHOD
        function obj = scenarioBuilder(varargin)
            % The scenario builder is designed to define a scenario for a
            % given number of objects. An instance of the scenario builder
            % class must be created first to ensure the

            % Access the OMAS-common properties
            addpath('environment/common');  % Access to common tools

            % Parse inputs against the builder
            [obj] = obj.configurationParser(obj,varargin);

            % Value checking
            assert(isnumeric(obj.objects),'Expecting a numeric value for field "objects"');
            assert(ischar(obj.name),'Expecting a string for field "name"');

            if obj.objects > 0
                % Initialise the scenario parameters
                obj.positions   = zeros(3,obj.objects);
                obj.velocities  = zeros(3,obj.objects);
                obj.quaternions = repmat([1;0;0;0],[1 obj.objects]);
            end
        end

        %% SCENARIO GENERATORS
        % PLANAR CONCENTRIC DISK %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = planarDisk(obj,varargin)
            % This function generates a configuration representing a given
            % number of objects distributed about a disk defined by a
            % central axis vector.
            % INPUTS:
            % pointA      - The first axis reference point.
            % pointB      - The second axis reference point; the ring center.
            % radius      - The scalar radius of the ring
            % velocity    - The velocity multiplier for the object set.
            % OUTPUTS:
            % definition  - The scenario definition

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'pointA',[0;0;-1],...
                                   'pointB',[0;0;0],...
                                   'zeroAngle',0,...                       % The rotation of the first node
                                   'velocity',0,...
                                   'scale',1);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(isnumeric(config.pointA) && numel(config.pointA) == 3,'The reference point must be a 3D Cartesian point.');
            assert(isnumeric(config.pointB) && numel(config.pointB) == 3,'The centroid point must be a 3D Cartesian point.');
            assert(isnumeric(config.zeroAngle) && numel(config.zeroAngle) == 1,'The zero angle must be a scalar.');
            assert(isnumeric(config.velocity) && numel(config.velocity) == 1,'The object velocity must be a scalar.');
            assert(isnumeric(config.scale) && numel(config.scale) == 1,'The object padding scale must be a scalar.');

            % ////////////// DESIGN THE DISK DISTRIBUTION /////////////////
            % Update scenario label
            obj.name = 'Distributed planar disk.';
            % Define object number
            obj.objects = config.objects;
            % If there is only one object, simply place it at the center
            if obj.objects == 1
                obj.positions(:,1) = config.pointB;
                obj.velocities(:,1) = [1;0;0]*config.velocity;
                obj.quaternions(:,1) = [1;0;0;0];
                return
            end
            remainingObjects = obj.objects;
            layer = 1;
            while layer < obj.objects && remainingObjects > 0
                % NEW LAYER RADIUS
                layerRadius = layer*config.scale;                          % New circle radius
                criticalLayerAngle = 2*asin(1/(2*layer));                  % New maximum nodal angle
                layerMax = floor(2*pi/criticalLayerAngle);                 % Maximum number on the new layer
                % DISTRIBUTE OBJECTS OVER THE LAYERS
                if layerMax > remainingObjects
                    objectsInLayer = remainingObjects;
                else
                    objectsInLayer = layerMax;
                end
                % BUILD THE NEXT LAYER
                [layerConfig] = obj.planarRing(...
                    'objects',objectsInLayer,...
                    'pointA',config.pointA,...
                    'pointB',config.pointB,...
                    'zeroAngle',config.zeroAngle,...
                    'radius',layerRadius);

                % EXTRACT THE LAYER GLOBAL PROPERTIES
                first = 1 + obj.objects - remainingObjects;
                last  = (first - 1) + objectsInLayer;
                indexVector = first:last;
                obj.positions(:,indexVector) = layerConfig.positions(:,1:objectsInLayer);
                obj.velocities(:,indexVector) = layerConfig.velocities(:,1:objectsInLayer);
                obj.quaternions(:,indexVector) = layerConfig.quaternions(:,1:objectsInLayer);
                % INCREMENT LAYER
                remainingObjects = remainingObjects - objectsInLayer;
                layer = layer + 1;
            end

        end
        % PLANAR RING OF STATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = planarRing(obj,varargin)
            % This function generates a configuration representing a given
            % number of objects destributed around a ring defined by a
            % central axis vector.
            % INPUTS:
            % pointA     - The first axis reference point.
            % pointB     - The second axis reference point; the ring center.
            % mod_radius - The scalar radius of the ring
            % velocity   - The velocity multiplier for the object set.
            % OUTPUTS:
            % definition - The scenario definition

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'pointA',[0;0;-1],...
                                   'pointB',[0;0;0],...
                                   'radius',10,...
                                   'zeroAngle',0,...                       % The rotation of the first node
                                   'offsetAngle',pi/5,...
                                   'velocity',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % Update scenario label
            obj.name = 'Planar ring.';
            % Define object number
            obj.objects = config.objects;

            % HAND NODAL ANGLE TO ringAngle FUNCTION
            [ obj ] = obj.planarAngle(...
                'objects',config.objects,...
                'pointA',config.pointA,...
                'pointB',config.pointB,...
                'radius',config.radius,...
                'zeroAngle',config.zeroAngle,...
                'offsetAngle',(2*pi)/config.objects,... % Define an angle for equal distribution
                'velocity',config.velocity);
        end
        % CONSTANT OFFSET ANGLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = planarAngle(obj,varargin)
            % GENERATE A SCENARIO DEFINITION
            % INPUTS:
            % pointA        - The first axis reference point.
            % pointB        - The second axis reference point; the ring center.
            % mod_radius    - The scalar radius of the ring
            % velMultiplier - The velocity multiplier for the object set.
            % offsetAngle   - The constant radial angle
            % OUTPUTS:
            % definition    - The scenario deifnition

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'pointA',[0;0;-1],...
                                   'pointB',[0;0;0],...
                                   'radius',10,...
                                   'zeroAngle',0,...                       % The rotation of the first node
                                   'offsetAngle',pi/5,...
                                   'velocity',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(isnumeric(config.pointA) && numel(config.pointA) == 3,'The reference point must be a 3D Cartesian point.');
            assert(isnumeric(config.pointB) && numel(config.pointB) == 3,'The center point must be a 3D Cartesian point.');
            assert(isnumeric(config.radius) && numel(config.radius) == 1,'The arc radius must be a scalar.');
            assert(isnumeric(config.zeroAngle) && numel(config.zeroAngle) == 1,'The initial radial angle must be a scalar.');
            assert(isnumeric(config.offsetAngle) && numel(config.offsetAngle) == 1,'The offset radial angle must be a scalar.');
            assert(isnumeric(config.velocity) && numel(config.velocity) == 1,'The object velocity must be a scalar.');

            % Update scenario label
            obj.name = sprintf('Planar offset angle of %srad.',num2str(config.offsetAngle));
            % Define object number
            obj.objects = config.objects;

            % GENERATE THE NODE(OBJECT) CARTAESIAN POSITIONS
            % GET THE AXIS VECTOR PROPERTIES
            localXAxis = config.pointB - config.pointA;
            unit_axisVector = unit(localXAxis);
            % GET THE RADIAL VECTOR
            perpVector = cross(unit_axisVector,[1;0;0]);
            if sum(perpVector) == 0
                perpVector = cross(unit_axisVector,[0;1;0]);
            end
            perpVector = unit(perpVector); % Re-normalise the perpendicular vector
            % SCALE THE UNIT RADIAL VECTOR
            perpVector = config.radius*perpVector;

            % GET THE NODE POINT SET
            nodalAngle = config.zeroAngle + pi/2; % pi/2 aligns the first object with the x-axis
            for node = 1:config.objects
                % ROTATE THE RADIAL VECTOR TO DEFINE NODAL POSITIONS
                radialVector = OMAS_geometry.rodriguesRotation(perpVector,unit_axisVector,nodalAngle);
                % NODAL POSITIONS
                globalPosition = config.pointB + radialVector;
                % NODAL VELOCITIES
                unitVelocity = - unit(radialVector);        % -ve sign to make concentric velocities +ve
                % NODAL ORIENTATIONS
                % We desire to orientate each object towards the center of
                % the ring. This is done by aligning the local body XY plane
                % With the plane of the ring. The local X axis is then
                % aligned with the cocentric global velocity.

                % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                localTriad = zeros(3);
                localTriad(:,1) = unit(unitVelocity);             % GET THE UNIT DIRECTION VECTOR
                localTriad(:,3) = unit(unit_axisVector);          % UNIT RING AXIS
                localTriad(:,2) = unit(cross(localTriad(:,1),localTriad(:,3)));

                % STORE THE OBJECTS IN SCENARIO DEFINITION
                obj.positions(:,node)   = globalPosition;
                obj.velocities(:,node)  = config.velocity*unitVelocity;
                obj.quaternions(:,node) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);

                % Object data is concatinated vertically
                % Increment the nodal position by the offset angle
                nodalAngle = nodalAngle - config.offsetAngle;
            end
        end
        % SPHERICAL EQUAL DISTRIBUTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = regularSphere(obj,varargin)
            % This function is designed to generate points equally
            % distributed about the circumference of a sphere.

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'radius',10,...
                                   'center',[0;0;0],...
                                   'velocity',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(isnumeric(config.radius) && numel(config.radius) == 1,'The sphere radius must be a scalar.');
            assert(isnumeric(config.center) && numel(config.center) == 3,'The sphere center must be a 3D Cartesian vector.');
            assert(isnumeric(config.velocity) && numel(config.velocity) == 1,'The object velocity must be a scalar.');

            % ATTEMPT TO ADD THE ICOSAHEDRALS LIBRARY
            try
                libraryPath = strcat(cd,'\scenarios\icosahedrals');
                addpath(libraryPath);
            catch
                error('Unable to load the icosahedral generator function, has it been moved?');
            end

            % UPDATE SCENARIO LABEL
            obj.name = 'Equally spaced spherical distribution.';
            obj.objects = config.objects;

            % GET SCALED ICOSAHEDRAL SPHERE
            % This function will generate a vector of points equally
            % distributed about the circumference of a sphere.
            [pointCloud,~] = spheretri(obj.objects);
            XYZ = (pointCloud.*config.radius)';
            XYZ = XYZ + config.center;                                      % Offset and scale the sphere

            % SORT THE POINTS
            XYZ = unique(XYZ','rows');                                     % Remove repeat entries
            XYZ = sortrows(XYZ,1);
            XYZ = XYZ';
            nodalPositions = XYZ();
            nodalPositions = nodalPositions(:,1:obj.objects);
            % ALLOCATE RANDOM SET TO OBJECT GLOBAL POSITIONS
            for index = 1:obj.objects
                % DETERMINE THE NODAL POSITIONS
                nodalPosition = config.center + nodalPositions(:,index);

                errorFlag = 1;
                while errorFlag == 1
                    % DEFINE GLOBAL VELOCITIES
                    unit_radii = unit(config.center - nodalPosition);
                    obj.velocities(:,index) = config.velocity*unit_radii;     % Cocentric velocity
                    % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                    localTriad = zeros(3);
                    localTriad(:,1) = unit_radii;                              % The body axis X direction
                    localTriad(:,3) = unit(cross(localTriad(:,1),obj.globalTriad(:,1)));      % The body axis Y direction
                    localTriad(:,2) = unit(cross(localTriad(:,1),localTriad(:,3)));           % The body axis Z direction
                    % DETERMINE IF A QUATERNION ROTATION EXISTS
                    try
                        % ASSIGN POSITION
                        obj.positions(:,index) = nodalPosition;
                        % ROTATION
                        obj.quaternions(:,index) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);
                        errorFlag = 0;
                    catch
                        % IF ROTATION IS INVALID, RESAMPLE THE POSITION
                        [nodalPosition,~] = datasample(XYZ,1,2,'Replace',false);
                    end
                end
            end
        end
        % COCENTRIC HELICAL STATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = helix(obj,varargin)

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'pointA',[0;0;0],...
                                   'pointB',[0;0;10],...
                                   'radius',10,...
                                   'offsetAngle',pi/5,...
                                   'velocities',0,...
                                   'plot',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(ischar(config.file),'The provided file path must be a string.');
            assert(isnumeric(config.pointA) && numel(config.pointA) == 3,'The origin point must be a 3D Cartesian point.');
            assert(isnumeric(config.pointB) && numel(config.pointB) == 3,'The terminal point must be a 3D Cartesian point.');
            assert(isnumeric(config.radius) && numel(config.radius) == 1,'The helical radius must be a scalar.');
            assert(isnumeric(config.offsetAngle) && numel(config.offsetAngle) == 1,'The helical angle must be a scalar.');
            assert(isnumeric(config.velocity) && numel(config.velocity) == 1,'The object velocity must be a scalar.');

            % DEFINE THE HELIX AXIS
            helixAxis = config.pointB - config.pointA;
            norm_helixAxis = norm(helixAxis);
            unit_helixAxis = helixAxis/norm_helixAxis;

            % THE UNIT SPACING
            objectNumber = obj.objects;
            axialSeparation = norm_helixAxis/objectNumber;

            % REFERENCE POSITION FOR FIRST NODE
            pointA_i = config.pointA - axialSeparation*unit_helixAxis;
            for i = 1:objectNumber
                % AXIAL PARAMETERS
                pointB_i = pointA_i + axialSeparation*unit_helixAxis;
                radialAngle = config.offsetAngle*(i-1);

                % DEFINE CONFIGURATION FOR OBJECT "i"
                [ ~,planarConfig ] = obj.planarAngle('objects',1,...
                                                     'pointA',pointA_i,...
                                                     'pointB',pointB_i,...
                                                     'radius',config.radius,...
                                                     'zeroAngle',radialAngle,...                       % The rotation of the first node
                                                     'offsetAngle',0);
                % DEFINE THE CONFIGURATION OF THE AGENT
                obj.positions(:,i) = planarConfig.position(:,1);
                obj.velocities(:,i) = planarConfig.velocity(:,1);
                obj.quaternions(:,i) = planarConfig.quaternion(:,1);
                % SHIFT UP THE AXIS VECTOR
                pointA_i = pointB_i;
            end

            obj.name = 'Cocentric helix.';
            % SAVE THE RESULTING SCENARIO
            scenarioConfig = struct('objects',obj.objects,...
                                    'name',obj.name,...
                                    'position',obj.positions,...
                                    'velocity',obj.velocities,...
                                    'quaternion',obj.quaternions);
        end
        % COLINEAR DISTRIBUTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = line(obj,varargin)
            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,...
                                   'pointA',[0;0;0],...
                                   'pointB',[0;0;10],...
                                   'heading',[1;0;0],...
                                   'velocities',1,...
                                   'plot',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(isnumeric(config.pointA) && numel(config.pointA) == 3,'Please provide a 3D point to the start of the line.');
            assert(isnumeric(config.pointB) && numel(config.pointB) == 3,'Please provide a 3D point to the end of the line.');
            assert(isnumeric(config.heading) && numel(config.heading) == 3,'Please provide a 3D vector defining the object heading');
            assert(isnumeric(config.velocities) && numel(config.velocities) == 1,'Please provide a scalar initial speed along the heading vector');

            % Label the scenario
            obj.name = 'Apposing line configuration';
            obj.objects = config.objects;

            % DISTRIBUTE THE OBJECTS EQUALLY ALONG THE AXIS
            if isnumeric(config.objects)
                nObjects = config.objects;
            else
                nObjects = numel(config.objects);
            end

            headingReference = unit(config.heading);    % Ensure the heading is a unit vector
            lineLength = norm(config.pointB - config.pointA);
            lineAxis = (config.pointB - config.pointA)/lineLength;
            % Linearly interpolate
            segments = linspace(0,lineLength,nObjects);

            % Design the parameters for each object
            for node = 1:nObjects
                % OBJECT POSITION
                obj.positions(:,node)  = config.pointA + segments(node)*lineAxis;
                % OBJECT VELOCITY
                obj.velocities(:,node) = config.velocities*config.heading;
                % ORIENTATION
                % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                localTriad = zeros(3);
                localTriad(:,1) = headingReference;                        % GET THE UNIT DIRECTION VECTOR
                localTriad(:,3) = unit_mex(obj.globalTriad(:,3)); % GLOBAL Z-AXIS
                localTriad(:,2) = unit_mex(cross(localTriad(:,1),localTriad(:,3)));
                % RESOLVE ROTATION QUATERNION
                obj.quaternions(:,node) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);
            end
        end

        % RANDOM SPHERE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = randomSphere(obj,varargin)
            % This fuction is designed to generate random positions
            % distributed about the circumference of a sphere.

            % DEFAULT CONFIGURATION
            defaultConfig = struct('objects',0,....
                                   'radius',10,...
                                   'center',[0;0;0],...
                                   'velocity',0);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % INPUT SANITY CHECK
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');
            assert(isnumeric(config.radius) && numel(config.radius) == 1,'Sphere radius must be a scalar.');
            assert(isnumeric(config.center) && numel(config.ceneter) == 3,'Sphere center must be a 3D Cartesian vector.');
            assert(isnumeric(config.velocity) && numel(config.velocity) == 1,'Agent velocity must be a scalar.');

            % ATTEMPT TO ADD THE ICOSAHEDRALS LIBRARY
            try
                libraryPath = strcat(cd,'\scenarios\icosahedrals');
                addpath(libraryPath);
            catch
                error('Unable to load the icosahedral generator function, has it been moved?');
            end

            % GET SCALED ICOSAHEDRAL SPHERE
            % This function will generate a vector of points equally
            % distributed about the circumference of a sphere.
            [pointCloud,~] = spheretri(obj.objects);
            XYZ = (pointCloud.*config.radius)';
            XYZ = XYZ + config.center;                                     % Offset and scale the sphere
            XYZ = unique(XYZ','rows');                                     % Remove repeat entries
            XYZ = XYZ';

            % RANDOMLY SAMPLE FROM THE POINT SET
            [nodalPositions,~] = datasample(XYZ,obj.objects,2,'Replace',false);   % Take one sample column from X coordinates

            % ALLOCATE RANDOM SET TO OBJECT GLOBAL POSITIONS
            for index = 1:config.objects
                % DETERMINE THE NODAL POSITIONS
                nodalPosition = config.center + nodalPositions(:,index);

                errorFlag = 1;
                while errorFlag == 1
                    % DEFINE GLOBAL VELOCITIES
                    unit_radii = unit(sphereCenter - nodalPosition);
                    obj.velocities(:,index) = config.velocity*unit_radii;     % Cocentric velocity
                    % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                    localTriad = zeros(3);
                    localTriad(:,1) = unit_radii;                              % The body axis X direction
                    localTriad(:,3) = unit(cross(localTriad(:,1),obj.globalTriad(:,1)));      % The body axis Y direction
                    localTriad(:,2) = unit(cross(localTriad(:,1),localTriad(:,3)));           % The body axis Z direction
                    % DETERMINE IF A QUATERNION ROTATION EXISTS
                    try
                        % ASSIGN POSITION
                        obj.positions(:,index) = nodalPosition;
                        % ROTATION
                        obj.quaternions(:,index) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);
                        errorFlag = 0;
                    catch
                        % IF ROTATION IS INVALID, RESAMPLE THE POSITION
                        [nodalPosition,~] = datasample(XYZ,1,2,'Replace',false);
                    end
                end
            end
            % SAVE THE RESULTING SCENARIO
            scenarioConfig = struct('objects',obj.objects,...
                                    'name',obj.name,...
                                    'position',obj.positions,...
                                    'velocity',obj.velocities,...
                                    'quaternion',obj.quaternions);
            % UPDATE SCENARIO LABEL
            obj.name = 'Random cocentric sphere';
        end
        % DEFAULT RANDOM STATE GENERATOR - NORMAL %%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = random(obj,varargin)
            % This function is a simple placeholder to facilitate the use
            % of .random for quick access to normally distributed state sets.
            [ obj ] = obj.randomNormal(varargin);
        end
        % RANDOM (NORMAL DISTRIBUTION) SCENARIO %%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = randomNormal(obj,varargin)
            % This function generates a random state vector set using a normal distribution.
            % assuming the order of the states is [x;y;z;u;v;w;phi;theta;psi;p;q;r]
            % INPUT:
            % obj          - The scenario definition object
            % OUTPUT:
            % randomStates - The state vector set [states x objects]
            % obj          - The returned scenario object

            % DEFAULT CONFIGURATION
            defaultConfig = struct(...
                'objects',0,....
                'velocity',0,...
                'positionGain',10,...
                'velocityGain',10,...
                'poseGain',pi);
            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % Input sanity check
            assert(isnumeric(config.objects) && config.objects > 0,'Expecting a scalar number of objects greater than zero.');

            % DEFINE A RANDOM STATE PARAMETERS
            obj.name = 'Normal-randomised';
            % The number of objects
            obj.objects = config.objects;

            stateNum = 9;
            positionTuple = 1:3;
            velocityTuple = 4:6;
            angleTuple = 7:9;
            stateGains = [config.positionGain,config.velocityGain,config.poseGain];
            normalObjectStates = zeros(stateNum,1);

            for index = 1:obj.objects
                % GENERATE RANDOM NUMBER SET
                maxR = max(randn(stateNum,1));
                normalObjectStates(:,index) = randn(stateNum,1)/maxR; % returns an n-by-n matrix of random numbers.

                % ALLOCATE STATE VALUES
                normalObjectStates(positionTuple,index) = stateGains(1)*normalObjectStates(positionTuple,index);
                normalObjectStates(velocityTuple,index) = stateGains(2)*normalObjectStates(velocityTuple,index);

                % NODAL POSITIONS
                obj.positions(:,index)  = normalObjectStates(positionTuple,index);
                % NODAL VELOCITIES
                obj.velocities(:,index) = normalObjectStates(velocityTuple,index);

                % NODAL ORIENTATIONS
                % We desire to orientate each object towards the center of
                % the ring. This is done by aligning the local body XY plane
                % With the plane of the ring. The local X axis is then
                % aligned with the cocentric global velocity.

                % ALLOCATE ANGULAR STATE VALUES
                normalObjectStates(angleTuple,index) = stateGains(3)*normalObjectStates(angleTuple,index);
                normalObjectStates(angleTuple(2:end),index) = 0;           % Remove pitch and yaw from random states (for stability)

                % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                localTriad = zeros(3);
                localTriad(:,1) = unit(normalObjectStates(velocityTuple,index));      % The body axis X direction
                localTriad(:,3) = unit(cross(localTriad(:,1),obj.globalTriad(:,1)));  % The body axis Y direction
                localTriad(:,2) = unit(cross(localTriad(:,1),localTriad(:,3)));       % The body axis Z direction
                obj.quaternions(:,index) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);
            end

        end
        % RANDOM (UNIFORM DISTRIBUTION) SCENARIO %%%%%%%%%%%%%%%%%%%%%%%%%%
        function [ obj ] = randomUniform(obj,varargin)
            % This function generates a random state vector set using a uniform distribution.
            % assuming the order of the states is [x;y;z;u;v;w;phi;theta;psi;p;q;r]
            % INPUT:
            % obj          - The scenario definition object
            % OUTPUT:
            % obj          - The returned scenario object

            % DEFAULT CONFIGURATION
            defaultConfig = struct(...
                'objects',0,....
                'velocity',0,...
                'positionGain',10,...
                'velocityGain',10,...
                'poseGain',pi);

            % PARSE CONFIGURATIONS
            [config] = obj.configurationParser(defaultConfig,varargin);

            % DEFINE A RANDOM STATE PARAMETERS
            obj.name = 'Uniformly-randomised';
            % The number of objects
            obj.objects = config.objects;

            stateNum = 9;
            positionTuple = 1:3;
            velocityTuple = 4:6;

            angleTuple = 7:9;
            stateGains = [config.positionGain,config.velocityGain,config.poseGain];
            uniformObjectStates = zeros(stateNum,1);

            for index = 1:obj.objects
                % GENERATE RANDOM NUMBER SET
                uniformObjectStates(:,index) = (rand([stateNum,1])-0.5).*2;       % returns an n-by-n matrix of random numbers.
                % ALLOCATE LINEAR STATE VALUES
                uniformObjectStates(positionTuple,index) = stateGains(1)*uniformObjectStates(positionTuple,index);
                uniformObjectStates(velocityTuple,index) = stateGains(2)*uniformObjectStates(velocityTuple,index);

                % NODAL POSITIONS
                obj.positions(:,index)  = uniformObjectStates(positionTuple,index);
                % NODAL VELOCITIES
                obj.velocities(:,index) = uniformObjectStates(velocityTuple,index);

                % NODAL ORIENTATIONS
                % We desire to orientate each object towards the center of
                % the ring. This is done by aligning the local body XY plane
                % With the plane of the ring. The local X axis is then
                % aligned with the cocentric global velocity.

                % ALLOCATE ANGULAR STATE VALUES
                uniformObjectStates(angleTuple,index) = stateGains(3)*uniformObjectStates(angleTuple,index);
                uniformObjectStates(angleTuple(2:end),index) = 0; % Remove pitch and yaw from random states (for stability)

                % GET THE LOCAL AXES (ROTATED IN THE GLOBAL FRAME)
                localTriad = zeros(3);
                localTriad(:,1) = unit(uniformObjectStates(velocityTuple,index));     % The body axis X direction
                localTriad(:,3) = unit(cross(localTriad(:,1),obj.globalTriad(:,1)));  % The body axis Y direction
                localTriad(:,2) = unit(cross(localTriad(:,1),localTriad(:,3)));       % The body axis Z direction
                obj.quaternions(:,index) = OMAS_geometry.getAnalyticalTriadRotation(obj.globalTriad,localTriad);
            end
        end
    end
    methods (Static)
        % PLOT SCENARIO FROM INITIALISED OBJECT-INDEX
        function [figureHandle] = plotObjectIndex(objectIndex)
            % This function is designed to plot the configured objectIndex
            % using their global properties and thier simulation
            % object-types

            global plotnum

            % DETERMINE PLOT PROPERTIES
            if ~exist('plotnum','var') || isempty(plotnum)
                plotnum = 1;    % Default to first plot
            end

            % GENERATE THE FIGURE
            figureHandle = figure(plotnum);
            axis vis3d;
            view([45,45]);
            xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
            hold on; grid on;
            ax = gca;
            set(ax,'FontSize',12,'fontWeight','bold');

            % OBJECT COUNTERS
            agents = 0; obstacles = 0; waypoints = 0;
            for index = 1:numel(objectIndex)
                % GET THE OBSTACLES GLOBAL PARAMETERS
                GLB = objectIndex{index}.GetGLOBAL();
                objectName  = objectIndex{index}.name;
                objectID    = objectIndex{index}.objectID;

                % GET THE ROTATION MATRIX GOING FROM BODY-GLOBAL
                [R_q] = OMAS_geometry.quaternionToRotationMatrix(GLB.quaternion);
                % PLOT THE OBJECT ORIENTATION TRIAD
                OMAS_graphics.drawTriad(figureHandle,GLB.position,R_q');
                % PLOT THE POSITIONS
                scatter3(GLB.position(1),GLB.position(2),GLB.position(3),'r');
                % PLOT THE VELOCITY VECTORS
                quiv = quiver3(GLB.position(1),GLB.position(2),GLB.position(3),...
                               GLB.velocity(1),GLB.velocity(2),GLB.velocity(3),'c');
                quiv.AutoScaleFactor = 1;
                % DETERMINE REPRESENTATION IF THE OBJECT HAS GEOMETRY
                if size(objectIndex{index}.GEOMETRY.vertices,1) < 1
                    % REPRESENT AS SPHERE WITH DEFINED RADIUS
                    [geometry] = OMAS_graphics.defineSphere(GLB.position,GLB.radius);
                    vertexData = geometry.vertices;
                    faceData = geometry.faces;
                    entityFaceAlpha = 0.3;
                    entityLineWidth = 0.05;
                    entityEdgeAlpha = 0.1;                                 % Show representative shapes with higher alpha
                else
                    vertexData = objectIndex{index}.GEOMETRY.vertices*R_q + GLB.position';
                    faceData   = objectIndex{index}.GEOMETRY.faces;
                    entityFaceAlpha = 0.8;
                    entityLineWidth = 1;
                    entityEdgeAlpha = 0.8; %0.2;
                end
                % REPRESENT GEOMETRY AS A PATCH
                entityHandle = patch(ax,...
                                    'Vertices',vertexData,...
                                    'Faces',faceData,...
                                    'EdgeColor','k',...
                                    'EdgeAlpha',entityEdgeAlpha,...
                                    'FaceAlpha',entityFaceAlpha,...
                                    'FaceLighting','gouraud',...
                                    'LineWidth',entityLineWidth);
                % PLOT REPRESENTATION
                switch GLB.type
                    case OMAS_objectType.agent
                        set(entityHandle,'FaceColor','b');
                        agents = agents + 1;
                    case OMAS_objectType.obstacle
                        set(entityHandle,'FaceColor','r');
                        obstacles = obstacles + 1;
                    case OMAS_objectType.waypoint
                        set(entityHandle,'FaceColor','g');
                        waypoints = waypoints + 1;
                    otherwise
                        set(entityHandle,'FaceColor','m');
                end
                % ADD ANNOTATION
                annotationText = sprintf('    %s [ID:%s]',objectName,num2str(objectID));
                text(GLB.position(1),GLB.position(2),GLB.position(3),char(annotationText));
            end
            % ADD TITLE
            titleStr = sprintf('Test scenario: %.0f agents, %.0f obstacles and %.0f waypoints.',agents,obstacles,waypoints);
            title(titleStr);
            hold off;
        end
        % PARSE A GENERIC INPUT SET AGAINST A DEFAULT CONFIG STRUCTURE
        function [config] = configurationParser(defaultConfig,scenarioParameters)
            % This function is designed to parse a generic set of user
            % inputs and allow them to be compared to a default input
            % structure. This should be called from a get_scenario file.

            % Call the generic parameter overrider
            [config] = GetParameterOverrides_recursive(defaultConfig,scenarioParameters);
        end
    end
end

repo: OpenMAS


text:
function sqi = psqi(signal,fs,btest,btot)
%pSQI Power of QRS SQI
%
% Returns the relative power on band P(5-20Hz)/P(5-45Hz). Operates in 1D or
% 2D vectors for some speeding up.
%
% Reference:
% Li, Q., Mark, R. G., Clifford, G. D., & Li. (2008). Robust heart rate
% estimation from multiple asynchronous noisy sources using signal quality
% indices and a Kalman filter. Physiol. Meas., 29(1), 15–32.
% http://doi.org/10.1088/0967-3334/29/1/002
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%
% Output:
%   sqi:            resulting sSQI for segment
%
% Fetal Extraction Toolbox, version 1.0, February 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Fernando Andreotti
% Dresden University of Technology, Institute of Biomedical Engineering
% fernando.andreotti@mailbox.tu-dresden.de
%
% Last updated : 09-03-2014
%
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

signal(isnan(signal)) = 0;

xdft = fft(detrend(signal),[],1);
xdft = xdft(1:floor(size(signal,1)/2+1),:);
xdft(2:end-1,:) = 2*xdft(2:end-1,:);
psdest = 1/(size(signal,1)*fs)*abs(xdft).^2;
freq = 0:fs/size(signal,1):fs/2;

% plot(freq,psdest);
% xlabel('Hz');
% grid on;
% title('Single-Sided Amplitude Spectrum of S(t)')
% xlabel('f (Hz)')
% ylabel('|P1(f)|')

if nargin < 3
    btest = [5 15];
    btot = [5 40];
end

pband = sum(psdest(freq>=btest(1)&freq<=btest(2),:));
ptot = sum(psdest(freq>=btot(1)&freq<=btot(2),:));
sqi = 1-(pband./ptot);

end

author: fernandoandreotti


text:
%  INTERNAL FUNCTION: change the parameter names from name_chain_state to name(chain,state)
%
%  ::
%
%    pname=param_name_to_param_texname(pname,chain_names)
%
%  Args:
%
%     - **pname** [char|cellstring: names of the parameters to change
%
%     - **chain_names** [cellstring]: names of the markov chains
%
%  Returns:
%     :
%
%     - **pname** [char|cellstring]: names of the changed parameters
%
%  See also:
%     - parser.param_texname_to_param_name
%
%

author: jmaih


text:
% EX_ARTICLE_SECTION_512: short example with a non-isoparametric approach, and a non-NURBS geometry.
%
% Example to solve the problem
%
%    - div ( grad (u)) = (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2)  in Omega
%                    u = 0                                                      on Gamma
%
% with                Omega = (1 < x^2+y^2 < 4) & (x > 0) & (y > 0)
% and exact solution      u = (x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x))
%
% This solves the example of Section 5.1.2 in the article
%
% C. De Falco, A. Reali, R. Vazquez
% GeoPDEs: a research tool for IsoGeometric Analysis of PDEs
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

geometry = geo_load ({@ring_polar_map, @ring_polar_map_der});

knots = kntuniform ([9 9], [4 4], [3 3]);
[qn, qw] = msh_set_quad_nodes (knots, msh_gauss_nodes ([5 5]));
msh = msh_cartesian (knots, qn, qw, geometry);

space = sp_bspline (knots, [4 4], msh);

mat = op_gradu_gradv_tp (space, space, msh);
rhs = op_f_v_tp (space, msh, @(x, y) (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2));

drchlt_dofs = [];
for iside = 1:4
  drchlt_dofs = union (drchlt_dofs, space.boundary(iside).dofs);
end
int_dofs = setdiff (1:space.ndof, drchlt_dofs);

u = zeros (space.ndof, 1);
u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);

sp_to_vtk (u, space, geometry, [20 20], 'laplace_solution.vts', 'u')
err = sp_l2_error (space, msh, u, @(x,y)(x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x)))

%!demo
%! ex_article_section_512;

%!test
%! geometry = geo_load ({@ring_polar_map, @ring_polar_map_der});
%! knots = kntuniform ([9 9], [4 4], [3 3]);
%! [qn, qw] = msh_set_quad_nodes (knots, msh_gauss_nodes ([5 5]));
%! msh = msh_cartesian (knots, qn, qw, geometry);
%! space = sp_bspline (knots, [4 4], msh);
%! mat = op_gradu_gradv_tp (space, space, msh, @(x, y) ones (size (x)));
%! rhs = op_f_v_tp (space, msh, @(x, y) (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2));
%! drchlt_dofs = [];
%! for iside = 1:4
%! drchlt_dofs = union (drchlt_dofs, space.boundary(iside).dofs);
%! end
%! int_dofs = setdiff (1:space.ndof, drchlt_dofs);
%! u = zeros (space.ndof, 1);
%! u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);
%! err = sp_l2_error (space, msh, u, @(x,y)(x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x)));
%! assert (msh.nel, 64)
%! assert (space.ndof, 144)
%! assert (err, 2.82510645627626e-07, 1e-16)


text:
function vis_topo(varargin)
% handles:
% handles.ics(1-8) has which ICs to plot

% get the updated stream buffer
W = evalin('base','W');
sphere = evalin('base','sphere');
Winv = inv(W*sphere);
handles = varargin{3};

it = mod(get(varargin{1},'TasksExecuted')-1,8)+1;
hstr = ['axesIC' int2str(it)];
hand = get(handles.(hstr),'children');
[map, cmin, cmax] = topoplotUpdate(Winv(:,handles.ics(it)), handles.chanlocs,'electrodes','off','gridscale',32);
set(hand(end),'CData',map);
set(handles.(hstr),'CLim',[cmin cmax]);

repo: REST


text:
function [ t, rank ] = ksubset_colex_successor ( k, n, t, rank )

%*****************************************************************************80
%
%% KSUBSET_COLEX_SUCCESSOR computes the K subset colex successor.
%
%  Discussion:
%
%    In the original code, there is a last element with no successor.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer K, the number of elements each K subset must
%    have.  1 <= K <= N.
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input/output, integer T(K), describes a K subset.  T(I) is the
%    I-th element.  The elements must be listed in DESCENDING order.
%    On input, T describes a K subset.
%    On output, T describes the next K subset in the ordering.
%    If the input T was the last in the ordering, then the output T
%    will be the first.
%
%    Input/output, integer RANK, the rank.
%    If RANK = -1 on input, then the routine understands that this is
%    the first call, and that the user wishes the routine to supply
%    the first element in the ordering, which has RANK = 0.
%    In general, the input value of RANK is increased by 1 for output,
%    unless the very last element of the ordering was input, in which
%    case the output value of RANK is 0.
%

%
%  Return the first element.
%
  if ( rank == -1 )
    for i = 1 : k
      t(i) = k + 1 - i;
    end
    rank = 0;
    return
  end
%
%  Check.
%
  ierror = ksubset_colex_check ( k, n, t );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KSUBSET_COLEX_SUCCESSOR - Fatal error!\n' );
    fprintf ( 1, '  The input array is illegal.\n' );
    fprintf ( 1, '  IERROR = %d\n', ierror );
    error ( 'KSUBSET_COLEX_SUCCESSOR - Fatal error!' );
  end

  for i = k - 1 : -1 : 1
    if ( t(k+1-i) + 1 < t(k-i) )
      t(k+1-i) = t(k+1-i) + 1;
      rank = rank + 1;
      return
    end
  end

  if ( t(1) < n )
    t(1) = t(1) + 1;
    for i = 1 : k - 1
      t(k+1-i) = i;
    end
    rank = rank + 1;
    return
  end
%
%  The last K subset was input.
%  Return the first one.
%
  for i = 1 : k
    t(i) = k + 1 - i;
  end

  rank = 0;

  return
end


text:
function [FF,I] = orient_outward(V,F,C)
  % ORIENT_OUTWARD Use a simple heuristic to maintain or flip each
  % manifold/orientable patch of a mesh. Assumes that independent patches have
  % already been oriented consistently up to sign.
  %
  % [FF,I] = orient_outward(V,F,C)
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   C  #F list of component ids
  % Outputs:
  %   FF  #F by 3 list of new, potentially flipped triangle indices
  %   I  #C list of bools whether patch was flipped
  %
  % See also: bfs_orient, manifold_patches
  %
  [FF,C] = bfs_orient(F);
  I = false(max(C),1);
  for c = 1:max(C)
    N = normalizerow(normals(V,FF(C==c,:)));
    A = doublearea(V,FF(C==c,:));
    BC = barycenter(V,FF(C==c,:));
    BCmean = A'*BC/sum(A);
    BC = bsxfun(@minus,BC,BCmean);
    ndot = sum(bsxfun(@times,A,sum(N.*BC,2)));
    if ndot<0
      FF(C==c,:) = fliplr(FF(C==c,:));
      I(c) = true;
    else
      I(c) = false;
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
% The script is an example for the function "inset"
% figure 1 is a plot of 'zero order bessel function of the first kind' in the range of x=0 to x=50.
% figure 2 is a close-up of this function around its first zero.
% figure 3 shows the function in the main plot, and the close-up in the
% inset plot.
%
% Moshe Lindner, August 2010 (C).

close all
x1=0:50;
y1=besselj(0,x1);
x2=2.35:.001:2.45;
y2=besselj(0,x2);
fig1=figure(1);
plot(x1,y1,'b','linewidth',2)
hold on
plot(x1,0*x1,':k')
set(gca,'fontsize',15)
title ('bessel function')
xlabel('X')
ylabel('Y')
fig2=figure(2);
plot(x2,y2,'r')
hold on
plot(x2,0*x2,':k')
title ('close-up')
[h_m h_i]=inset(fig1,fig2);

set(h_i,'xtick',2.35:.025:2.45,'xlim',[2.35,2.45])


text:
% Get matlab type corresponding to each nrrd type
% Date: November 2017
% Author: Gaetan Rensonnet
function matlabdatatype = nrrd_getMatlabDataType(metaType)

% Determine the datatype
switch (metaType)
    case {'signed char', 'int8', 'int8_t'}
        matlabdatatype = 'int8';

    case {'uchar', 'unsigned char', 'uint8', 'uint8_t'}
        matlabdatatype = 'uint8';

    case {'short', 'short int', 'signed short', 'signed short int', ...
            'int16', 'int16_t'}
        matlabdatatype = 'int16';

    case {'ushort', 'unsigned short', 'unsigned short int', 'uint16', ...
            'uint16_t'}
        matlabdatatype = 'uint16';

    case {'int', 'signed int', 'int32', 'int32_t'}
        matlabdatatype = 'int32';

    case {'uint', 'unsigned int', 'uint32', 'uint32_t'}
        matlabdatatype = 'uint32';

    case {'longlong', 'long long', 'long long int', 'signed long long', ...
            'signed long long int', 'int64', 'int64_t'}
        matlabdatatype = 'int64';

    case {'ulonglong', 'unsigned long long', 'unsigned long long int', ...
            'uint64', 'uint64_t'}
        matlabdatatype = 'uint64';

    case {'float'}
        matlabdatatype = 'single';

    case {'double'}
        matlabdatatype = 'double';

    case {'block'}
        error('Data type ''block'' is currently not supported.\n');

    otherwise
        assert(false, 'Unknown datatype')
end
end


text:
function p=v_lin2pcma(x,m,s)
%V_LIN2PCMA Convert linear PCM to A-law P=(X,M,S)
%	pcma = v_lin2pcma(lin) where lin contains a vector
%	or matrix of signal values.
%	The input values will be converted to integer
%	A-law pcm vlues in the range 0 to 255 and the XORed with m
%	(default m=85).
%
%	Input values are multiplied by the scale factor s:
%
%		   s		Input Range
%
%		   1		+-4096
%		2017.396342	+-2.03033976 (default)
%		4096		+-1
%
%	Input values outside the selected range will be clipped.
%
%	The default value of the scale factor is 2017.396342 which equals
%	sqrt((1120^2 + 2624^2)/2). This factor follows ITU standard G.711 and
%	the sine wave with PCM-A values [225 244 244 225 97 116 116 97]
%	has a mean square value of unity corresponding to 0 dBm0.
%
%	See also PCMA2LIN, LIN2PCMA, LIN2PCMU

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lin2pcma.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
  s=2017.396342;
  if nargin<2 m=85; end
end

y=x*pow2(s,-6);
y=(abs(y+63)-abs(y-63))/2;
q=floor((y+64)/64);
[a,e]=log2(abs(y));
d = (e+abs(e))/2;
p=128*q+16*d+floor(pow2(a,e-d+5));
if m p=bitxor(p,m); end;

author: ImperialCollegeLondon

repo: sap-voicebox