text:
function [ cluster, cluster_center, cluster_population, cluster_energy, ...
  it_num ] = kmeans_02 ( dim_num, point_num, cluster_num, it_max, point, ...
  cluster_center )

%*****************************************************************************80
%
%% KMEANS_02 applies the K-Means algorithm.
%
%  Discussion:
%
%    The routine attempts to divide POINT_NUM points in
%    DIM_NUM-dimensional space into CLUSTER_NUM clusters so that the within
%    cluster sum of squares is minimized.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    Original FORTRAN77 by John Hartigan, Manchek Wong.
%    FORTRAN90 version by John Burkardt.
%
%  Reference:
%
%    John Hartigan, Manchek Wong,
%    Algorithm AS 136:
%    A K-Means Clustering Algorithm,
%    Applied Statistics,
%    Volume 28, Number 1, 1979, pages 100-108.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, integer IT_MAX, the maximum number of iterations.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates
%    of the points.
%
%    Input, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Output, integer CLUSTER(POINT_NUM), the cluster each
%    point belongs to.
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number
%    of points in each cluster.
%
%    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the
%    within-cluster sum of squares.
%
%    Output, integer IT_NUM, the number of iterations taken.
%
  it_num = 0;

  if ( cluster_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KMEANS_02 - Fatal error!\n' );
    fprintf ( 1, '  CLUSTER_NUM < 1.\n' );
    error ( 'KMEANS_02 - Fatal error!' )
  end

  if ( point_num <= cluster_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KMEANS_02 - Fatal error!\n' );
    fprintf ( 1, '  POINT_NUM <= CLUSTER_NUM.\n' );
    error ( 'KMEANS_02 - Fatal error!' )
  end

  if ( dim_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KMEANS_02 - Fatal error!\n' );
    fprintf ( 1, '  DIM_NUM < 1.\n' );
    error ( 'KMEANS_02 - Fatal error!' )
  end

  if ( point_num < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KMEANS_02 - Fatal error!\n' );
    fprintf ( 1, '  POINT_NUM < 1.\n' );
    error ( 'KMEANS_02 - Fatal error!' )
  end
%
%  For each point I, find its two closest centers, CLUSTER(I) and CLUSTER2(I).
%  Assign it to CLUSTER(I).
%
  for i = 1 : point_num

    cluster(i) = 1;
    cluster2(i) = 2;

    for il = 1 : 2
      dt(il) = sum ( ...
        ( point(1:dim_num,i) - cluster_center(1:dim_num,il) ).^2 );
    end

    if ( dt(2) < dt(1) )
      cluster(i) = 2;
      cluster2(i) = 1;
      temp = dt(1);
      dt(1) = dt(2);
      dt(2) = temp;
    end

    for l = 3 : cluster_num

      db = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l) ).^2 );

      if ( db < dt(1) )
        dt(2) = dt(1);
        cluster2(i) = cluster(i);
        dt(1) = db;
        cluster(i) = l;
      elseif ( db < dt(2) )
        dt(2) = db;
        cluster2(i) = l;
      end

    end

  end
%
%  Update cluster centers to be the average of points contained within them.
%
  cluster_population(1:cluster_num) = 0;
  cluster_center(1:dim_num,1:cluster_num) = 0.0;

  for i = 1 : point_num
    l = cluster(i);
    cluster_population(l) = cluster_population(l) + 1;
    cluster_center(1:dim_num,l) = cluster_center(1:dim_num,l) ...
      + point(1:dim_num,i);
  end
%
%  Check to see if there is any empty cluster.
%
  for l = 1 : cluster_num

    if ( cluster_population(l) == 0 )
      j = 1 + floor ( point_num * rand ( ) );
      cluster_center(1:dim_num) = point(1:dim_num,j)
    else
      cluster_center(1:dim_num,l) = cluster_center(1:dim_num,l) ...
        / cluster_population(l);
    end
%
%  Initialize AN1, AN2, ITRAN and NCP
%  AN1(L) = CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
%  AN2(L) = CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
%  ITRAN(L) = 1 if cluster L is updated in the quick-transfer stage,
%           = 0 otherwise
%  In the optimal-transfer stage, NCP(L) stores the step at which
%  cluster L is last updated.
%  In the quick-transfer stage, NCP(L) stores the step at which
%  cluster L is last updated plus POINT_NUM.
%
    an2(l) = cluster_population(l) / ( cluster_population(l) + 1 );

    if ( 1 < cluster_population(l) )
      an1(l) = cluster_population(l) / ( cluster_population(l) - 1 );
    else
      an1(l) = Inf;
    end

    itran(l) = 1;
    ncp(l) = -1;

  end

  indx = 0;
  ifault = 2;
  it_num = 0;

  d = zeros(point_num);
  live = zeros(cluster_num);

  while ( it_num < it_max )

    it_num = it_num + 1;
%
%  In this stage, there is only one pass through the data.   Each
%  point is re-allocated, if necessary, to the cluster that will
%  induce the maximum reduction in within-cluster sum of squares.
%
    [ cluster_center, cluster, cluster2, cluster_population, an1, an2, ...
      ncp, d, itran, live, indx ] = kmeans_02_optra ( dim_num, point_num, ...
      cluster_num, point, cluster_center, cluster, cluster2, cluster_population, ...
      an1, an2, ncp, d, itran, live, indx );
%
%  Stop if no transfer took place in the last POINT_NUM optimal transfer steps.
%
    if ( indx == point_num )
      ifault = 0;
      break
    end
%
%  Each point is tested in turn to see if it should be re-allocated
%  to the cluster to which it is most likely to be transferred,
%  CLUSTER2(I), from its present cluster, CLUSTER(I).   Loop through the
%  data until no further change is to take place.
%
    [ cluster_center, cluster, cluster2, cluster_population, an1, ...
      an2, ncp, d, itran, indx ] = kmeans_02_qtran ( dim_num, point_num, ...
      cluster_num, point, cluster_center, cluster, cluster2, cluster_population, ...
      an1, an2, ncp, d, itran, indx );
%
%  If there are only two clusters, there is no need to re-enter the
%  optimal transfer stage.
%
    if ( cluster_num == 2 )
      ifault = 0;
      break
    end
%
%  NCP has to be set to 0 before entering OPTRA.
%
    ncp(1:cluster_num) = 0;

  end

  if ( ifault == 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'KMEANS_02 - Warning!\n' );
    fprintf ( 1, '  Maximum number of iterations reached\n' );
    fprintf ( 1, '  without convergence.\n' );
  end
%
%  Compute the within-cluster sum of squares for each cluster.
%
  cluster_center(1:dim_num,1:cluster_num) = 0.0;

  for i = 1 : point_num
    cluster_center(1:dim_num,cluster(i)) = ...
      cluster_center(1:dim_num,cluster(i)) + point(1:dim_num,i);
  end

  for j = 1 : dim_num
    cluster_center(j,1:cluster_num) = cluster_center(j,1:cluster_num) ...
      ./ cluster_population(1:cluster_num);
  end
%
%  Compute the cluster energies.
%
  cluster_energy(1:cluster_num) = 0.0;

  for i = 1 : point_num

    j = cluster(i);

    cluster_energy(j) = cluster_energy(j) + sum ( ...
      ( point(1:dim_num,i) - cluster_center(1:dim_num,j) ).^2 );

  end

  return
end
function [ cluster_center, cluster, cluster2, cluster_population, an1, an2, ...
  ncp, d, itran, live, indx ] = kmeans_02_optra ( dim_num, point_num, ...
  cluster_num, point, cluster_center, cluster, cluster2, cluster_population, ...
  an1, an2, ncp, d, itran, live, indx )

%*****************************************************************************80
%
%% KMEANS_02_OPTRA carries out the optimal transfer stage.
%
%  Discussion:
%
%    Each point is re-allocated, if necessary, to the cluster that
%    will induce a maximum reduction in the within-cluster sum of
%    squares.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    Original FORTRAN77 by John Hartigan, Manchek Wong.
%    FORTRAN90 version by John Burkardt.
%
%  Reference:
%
%    John Hartigan, Manchek Wong,
%    Algorithm AS 136:
%    A K-Means Clustering Algorithm,
%    Applied Statistics,
%    Volume 28, Number 1, 1979, pages 100-108.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of
%    the points.
%
%    Input, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Input, integer CLUSTER(POINT_NUM), the cluster
%    each point belongs to.
%
%    Input, integer CLUSTER2(POINT_NUM), the cluster
%    to which each point is most likely to be transferred to.
%
%    Input, integer CLUSTER_POPULATION(CLUSTER_NUM),
%    the number of points in each cluster.
%
%    Input, real AN1(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
%
%    Input, real AN2(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
%
%    Input, integer NCP(CLUSTER_NUM), ?
%
%    Input, real D(POINT_NUM), ?
%
%    Input, integer ITRAN(CLUSTER_NUM),
%    1 if cluster L is updated in the quick-transfer stage,
%    0 otherwise.  Reset to 0 on output.
%
%    Input, integer LIVE(CLUSTER_NUM), ?
%
%    Input, integer INDX, ?
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Output, integer CLUSTER(POINT_NUM), the cluster
%    each point belongs to.
%
%    Output, integer CLUSTER2(POINT_NUM), the cluster
%    to which each point is most likely to be transferred to.
%
%    Output, integer CLUSTER_POPULATION(CLUSTER_NUM),
%    the number of points in each cluster.
%
%    Output, real AN1(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
%
%    Output, real AN2(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
%
%    Output, integer NCP(CLUSTER_NUM), ?
%
%    Output, real D(POINT_NUM), ?
%
%    Output, integer ITRAN(CLUSTER_NUM),
%    1 if cluster L is updated in the quick-transfer stage,
%    0 otherwise.  Reset to 0 on output.
%
%    Output, integer LIVE(CLUSTER_NUM), ?
%
%    Output, integer INDX, ?
%

%
%  If cluster L is updated in the last quick-transfer stage, it
%  belongs to the live set throughout this stage.   Otherwise, at
%  each step, it is not in the live set if it has not been updated
%  in the last POINT_NUM optimal transfer steps.
%
  for l = 1 : cluster_num
    if ( itran(l) == 1 )
      live(l) = point_num + 1;
    end
  end

  for i = 1 : point_num

    indx = indx + 1;
    l1 = cluster(i);
    l2 = cluster2(i);
    ll = l2;
%
%  If point I is the only member of cluster L1, no transfer.
%
    if ( 1 < cluster_population(l1) )
%
%  If L1 has been updated in this stage, re-compute D(I).
%
      if ( ncp(l1) ~= 0 )
        d(i) = an1(l1) * sum ( ...
          ( point(1:dim_num,i) - cluster_center(1:dim_num,l1) ).^2 );
      end
%
%  Find the cluster with minimum R2.
%
      r2 = an2(l2) * sum ( ...
        ( point(1:dim_num,i) - cluster_center(1:dim_num,l2) ).^2 );

      for l = 1 : cluster_num
%
%  If LIVE(L1) <= I, then L1 is not in the live set.   If this is
%  true, we only need to consider clusters that are in the live set
%  for possible transfer of point I.
%
%  Otherwise, we need to consider all possible clusters.
%
        if ( ( i < live(l1) | i < live(l) ) & ...
             l ~= l1 & ...
             l ~= ll )

          rr = r2 / an2(l);

          dc = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l) ).^2 );

          if ( dc < rr )
            r2 = dc * an2(l);
            l2 = l;
          end

        end

      end
%
%  If no transfer is necessary, L2 is the new CLUSTER2(I).
%
      if ( d(i) <= r2 )

        cluster2(i) = l2;

      else
%
%  Update cluster centers, LIVE, NCP, AN1 and AN2 for clusters L1 and
%  L2, and update CLUSTER(I) and CLUSTER2(I).
%
        indx = 0;
        live(l1) = point_num + i;
        live(l2) = point_num + i;
        ncp(l1) = i;
        ncp(l2) = i;
        al1 = cluster_population(l1);
        alw = al1 - 1.0;
        al2 = cluster_population(l2);
        alt = al2 + 1.0;

        cluster_center(1:dim_num,l1) = ( cluster_center(1:dim_num,l1) * al1 ...
          - point(1:dim_num,i) ) / alw;

        cluster_center(1:dim_num,l2) = ( cluster_center(1:dim_num,l2) * al2 ...
          + point(1:dim_num,i) ) / alt;

        cluster_population(l1) = cluster_population(l1) - 1;
        cluster_population(l2) = cluster_population(l2) + 1;
        an2(l1) = alw / al1;

        if ( 1.0 < alw )
          an1(l1) = alw / ( alw - 1.0 );
        else
          an1(l1) = Inf;
        end

        an1(l2) = alt / al2;
        an2(l2) = alt / ( alt + 1.0 );
        cluster(i) = l2;
        cluster2(i) = l1;

      end

    end

    if ( indx == point_num )
      return
    end

  end
%
%  ITRAN(L) = 0 before entering QTRAN.
%
  itran(1:cluster_num) = 0;
%
%  LIVE(L) has to be decreased by POINT_NUM before re-entering OPTRA.
%
  live(1:cluster_num) = live(1:cluster_num) - point_num;

  return
end
function [ cluster_center, cluster, cluster2, cluster_population, an1, ...
  an2, ncp, d, itran, indx ] = kmeans_02_qtran ( dim_num, point_num, ...
  cluster_num, point, cluster_center, cluster, cluster2, cluster_population, ...
  an1, an2, ncp, d, itran, indx )

%*****************************************************************************80
%
%% KMEANS_02_QTRAN carries out the quick transfer stage.
%
%  Discussion:
%
%    For each point I, CLUSTER(I) and CLUSTER2(I) are switched, if necessary,
%    to reduce within-cluster sum of squares.  The cluster centers are
%    updated after each step.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2009
%
%  Author:
%
%    Original FORTRAN77 by John Hartigan, Manchek Wong.
%    FORTRAN90 version by John Burkardt.
%
%  Reference:
%
%    John Hartigan, Manchek Wong,
%    Algorithm AS 136:
%    A K-Means Clustering Algorithm,
%    Applied Statistics,
%    Volume 28, Number 1, 1979, pages 100-108.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer CLUSTER_NUM, the number of clusters.
%
%    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates
%    of the points.
%
%    Input, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Input, integer CLUSTER(POINT_NUM), the cluster
%    each point belongs to.
%
%    Input, integer CLUSTER2(POINT_NUM), the cluster to which
%    each point is most likely to be transferred to.
%
%    Input, integer CLUSTER_POPULATION(CLUSTER_NUM),
%    the number of points in each cluster.
%
%    Input, real AN1(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
%
%    Input, real AN2(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
%
%    Input, integer NCP(CLUSTER_NUM), ?
%
%    Input, real D(POINT_NUM), ?
%
%    Input, integer ITRAN(CLUSTER_NUM),
%    1 if cluster L is updated in the quick-transfer stage,
%    0 otherwise.
%
%    Input, integer INDX, is set to 0 if any updating occurs.
%
%    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM),
%    the cluster centers.
%
%    Output, integer CLUSTER(POINT_NUM), the cluster
%    each point belongs to.
%
%    Output, integer CLUSTER2(POINT_NUM), the cluster to which
%    each point is most likely to be transferred to.
%
%    Output, integer CLUSTER_POPULATION(CLUSTER_NUM),
%    the number of points in each cluster.
%
%    Output, real AN1(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
%
%    Output, real AN2(CLUSTER_NUM),
%    CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
%
%    Output, integer NCP(CLUSTER_NUM), ?
%
%    Output, real D(POINT_NUM), ?
%
%    Output, integer ITRAN(CLUSTER_NUM),
%    1 if cluster L is updated in the quick-transfer stage,
%    0 otherwise.
%
%    Output, integer INDX, is set to 0 if any updating occurs.
%

%
%  In the optimal transfer stage, NCP(L) indicates the step at which
%  cluster L is last updated.   In the quick transfer stage, NCP(L)
%  is equal to the step at which cluster L is last updated plus POINT_NUM.
%
  count = 0;
  step = 0;

  while ( 1 )

    for i = 1 : point_num

      count = count + 1;
      step = step + 1;
      l1 = cluster(i);
      l2 = cluster2(i);
%
%  If point I is the only member of cluster L1, no transfer.
%
      if ( 1 < cluster_population(l1) )
%
%  If NCP(L1) < STEP, no need to re-compute distance from point I to
%  cluster L1.   Note that if cluster L1 is last updated exactly POINT_NUM
%  steps ago, we still need to compute the distance from point I to
%  cluster L1.
%
        if ( step <= ncp(l1) )
          d(i) = an1(l1) * sum ( ...
            ( point(1:dim_num,i) - cluster_center(1:dim_num,l1) ).^2 );
        end
%
%  If STEP >= both NCP(L1) and NCP(L2) there will be no transfer of
%  point I at this step.
%
        if ( step < ncp(l1) | step < ncp(l2) )

          r2 = d(i) / an2(l2);

          dd = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l2) ).^2 );
%
%  Update cluster centers, NCP, CLUSTER_POPULATION, ITRAN, AN1 and AN2
%  for clusters L1 and L2.   Also update CLUSTER(I) and CLUSTER2(I).
%
%  Note that if any updating occurs in this stage, INDX is set back to 0.
%
          if ( dd < r2 )

            count = 0;
            indx = 0;
            itran(l1) = 1;
            itran(l2) = 1;
            ncp(l1) = step + point_num;
            ncp(l2) = step + point_num;
            al1 = cluster_population(l1);
            alw = al1 - 1.0;
            al2 = cluster_population(l2);
            alt = al2 + 1.0;

            cluster_center(1:dim_num,l1) = ...
              ( cluster_center(1:dim_num,l1) * al1 ...
              - point(1:dim_num,i) ) / alw;

            cluster_center(1:dim_num,l2) = ...
              ( cluster_center(1:dim_num,l2) * al2 ...
              + point(1:dim_num,i) ) / alt;

            cluster_population(l1) = cluster_population(l1) - 1;
            cluster_population(l2) = cluster_population(l2) + 1;
            an2(l1) = alw / al1;

            if ( 1.0 < alw )
              an1(l1) = alw / ( alw - 1.0 );
            else
              an1(l1) = Inf;
            end

            an1(l2) = alt / al2;
            an2(l2) = alt / ( alt + 1.0 );
            cluster(i) = l2;
            cluster2(i) = l1;

          end

        end

      end
%
%  If no re-allocation took place in the last POINT_NUM steps, return.
%
      if ( count == point_num )
        return
      end

    end

  end

  return
end


text:
function varargout=lpczz2ss(varargin)
%
% For calling details please see v_lpczz2ss.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpczz2ss.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpczz2ss.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpczz2ss(varargin{:});
else
    v_lpczz2ss(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [CV,CT,TV,TT,w,SV,SF] = solid(V,F,varargin)
  % SOLID Clean a mesh using the generalized winding number pipeline with some
  % default settings [Jacobson et al. 2013]: "Shoot for the moon!!!!"
  %
  % [CV,CT] = solid(V,F)
  % [CV,CT] = solid(V,F,'ParameterName',ParameterValue,..)
  %
  % Inputs:
  %   V  #V by 3 input mesh positions
  %   F  #F by 3 list of triangle indices into V
  %   Optional:
  %     'SurfaceOnly' followed by whether to return CF  #CF by 3 list of
  %       surface triangle indices instead of CT, and remove rows unrefernced
  %       by CF from CV. Equivalent to running:
  %           [CV,CT] = winding_number_clean(V,F);
  %           CF = boundary_faces(CT);
  %           [CV,I] = remove_unreferenced(TV,CF)
  %           CF = I(CF);
  %       {false}.
  % Outputs:
  %   CV  #CV by 3 output mesh positions
  %   CT  #CT by 3 list of tet indices into CV
  %   TV  #TV by 3 list of tetrahedral mesh vertex positions
  %   TT  #TT by 3 list of tetrahedral mesh tet indices
  %   w  #TT list of winding numbers
  %   SV  #SV by 3 list of self-intersected "clean", but non-solid mesh indices
  %   SF  #SF by 3 list of triangle indices into SV
  %

  surface_only = false;
  % parse optional input parameters
  % Map of parameter names to variable names
  params_to_variables = containers.Map( {'SurfaceOnly'},{'surface_only'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin))
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % shoot for the moon with winding number on default settings
  [SV,SF] = clean(V,F,'Quiet',true,'Single',true);
  tetgen_flags_alts = {'-Y','-Y -T1e-16','','-T1e-16'};
  TT = [];
  for t = 1:numel(tetgen_flags_alts)
    tetgen_flags = tetgen_flags_alts{t};
    try
      [TV,TT] = cdt(SV,SF,'TetgenFlags',tetgen_flags,'Quiet',true);
      fprintf('cdt succeeded with "%s"...\n',tetgen_flags);
      break;
    catch
      fprintf('cdt failed with "%s"...\n',tetgen_flags);
    end
  end
  if isempty(TT)
    error('All tetgen flags alternatives failed');
  end
  w = winding_number(V,F,barycenter(TV,TT));
  vol = volume(TV,TT);
  avg_w = w'*(vol/sum(vol));
  if avg_w < 0
    fprintf('inside out?...\n');
    w = -w;
    avg_w = -avg_w;
  end
  CT = TT(w>avg_w,:);

  if surface_only
    CT = boundary_faces(CT);
  end
  [CV,I] = remove_unreferenced(TV,CT);
  CT = I(CT);

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function rec = compute_recovery_rate(x,x0)

% compute_recovery_rate - compute the number of recovered coefficients
%
%   rec = compute_recovery_rate(x,x0);
%
%   x is the original signal with 0/1 coefficients
%   x0 is the recovered signal.
%
%   rec=1 means perfect recovery and rec=0 means no recovery (bad).
%
%   Copyright (c) 2007 Gabriel Peyre

% works for positive coefficients 0/1

s = sum(x);
x = double( x(:)>0 );

[tmp,I] = sort(x0(:));
x0 = x0 * 0; x0(I(end-s+1:end)) = 1;
rec = sum(x0.*x)/s;

author: gpeyre

repo: matlab-toolboxes


text:
function [Population,Fitness] = EnvironmentalSelection(Population,N)
% The environmental selection of SPEA2+SDE

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Calculate the fitness of each solution
    Fitness = CalFitness(Population.objs);

    %% Environmental selection
    Next = Fitness < 1;
    if sum(Next) < N
        [~,Rank] = sort(Fitness);
        Next(Rank(1:N)) = true;
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
    end
    % Population for next generation
    Population = Population(Next);
    Fitness    = Fitness(Next);
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    N = size(PopObj,1);

    %% Calculate the shifted distance between each two solutions
    Distance = inf(N);
    for i = 1 : N
        SPopObj = max(PopObj,repmat(PopObj(i,:),N,1));
        for j = [1:i-1,i+1:N]
            Distance(i,j) = norm(PopObj(i,:)-SPopObj(j,:));
        end
    end

    %% Truncation
    Del = false(1,N);
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function lik = lik_multinom(varargin)
%LIK_MULTINOM  Create a multinomial likelihood structure
%
%  Description
%    LIK = LIK_MULTINOM creates multinomial likelihood for multi-class
%    count data. The observed numbers in each class with C classes is
%    given as 1xC vector.
%
%    The likelihood is defined as follows:
%                              __ n                __ C
%      p(y|f^1, ..., f^C, z) = || i=1 [ gamma(N+1) || c=1 p_i^c^(y_i^c)/gamma(y_i^c+1)]
%
%    where p_i^c = exp(f_i^c)/ (sum_c=1^C exp(f_i^c)) is the succes
%    probability for class c, which is a function of the latent variable
%    f_i^c for the corresponding class and N=sum(y) is the number of trials.
%
%  See also
%    GP_SET, LIK_*
%
% Copyright (c) 2010 Jaakko Riihimäki, Pasi Jylänki
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2010 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_MULTINOM';
  ip.addOptional('lik', [], @isstruct);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'Multinom';
    lik.nondiagW = true;
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Multinom')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_multinom_pak;
    lik.fh.unpak = @lik_multinom_unpak;
    lik.fh.ll = @lik_multinom_ll;
    lik.fh.llg = @lik_multinom_llg;
    lik.fh.llg2 = @lik_multinom_llg2;
    lik.fh.llg3 = @lik_multinom_llg3;
    lik.fh.tiltedMoments = @lik_multinom_tiltedMoments;
    lik.fh.predy = @lik_multinom_predy;
    lik.fh.invlink = @lik_multinom_invlink;
    lik.fh.recappend = @lik_multinom_recappend;
  end

end

function [w,s,h] = lik_multinom_pak(lik)
%LIK_MULTINOM_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_MULTINOM_PAK(LIK) takes a likelihood structure LIK and
%    returns an empty verctor W. If Multinom likelihood had
%    parameters this would combine them into a single row vector
%    W (see e.g. lik_negbin). This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%
%  See also
%    LIK_MULTINOM_UNPAK, GP_PAK

  w = []; s = {};h=[];
end

function [lik, w] = lik_multinom_unpak(lik, w)
%LIK_MULTINOM_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_MULTINOM_UNPAK(W, LIK) Doesn't do anything.
%
%    If Multinom likelihood had parameters this would extracts them
%    parameters from the vector W to the LIK structure. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%
%  See also
%    LIK_MULTINOM_PAK, GP_UNPAK

  lik=lik;
  w=w;
end

function ll = lik_multinom_ll(lik, y, f, z)
%LIK_MULTINOM_LL  Log likelihood
%
%  Description
%    LL = LIK_MULTINOM_LL(LIK, Y, F) takes a likelihood structure
%    LIK, class counts Y (NxC matrix), and latent values F (NxC
%    matrix). Returns the log likelihood, log p(y|f,z). This
%    subfunction is needed when using Laplace approximation or
%    MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria
%    (DIC, WAIC) computations.
%
%  See also
%    LIK_MULTINOM_LLG, LIK_MULTINOM_LLG3, LIK_MULTINOM_LLG2, GPLA_E

  f=reshape(f,size(y));
  expf = exp(f);
  p = expf ./ repmat(sum(expf,2),1,size(expf,2));
  N = sum(y,2);

  ll = sum(gammaln(N+1) - sum(gammaln(y+1),2) + sum(y.*log(p),2) );

end

function llg = lik_multinom_llg(lik, y, f, param, z)
%LIK_MULTINOM_LLG    Gradient of the log likelihood
%
%  Description
%    LLG = LIK_MULTINOM_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F. Returns
%    the gradient of the log likelihood with respect to PARAM. At
%    the moment PARAM can be 'param' or 'latent'. This subfunction
%    is needed when using Laplace approximation or MCMC for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG2, LIK_MULTINOM_LLG3, GPLA_E

  f=reshape(f,size(y));
  C = size(y,2);
  expf2 = exp(f);
  N=sum(y, 2);
  pi2 = (N*ones(1,C)).*expf2./(sum(expf2, 2)*ones(1,C));
  pi_vec=pi2(:);
  llg = y(:)-pi_vec;

end

function [pi_vec, pi_mat] = lik_multinom_llg2(lik, y, f, param, z)
%LIK_MULTINOM_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_MULTINOM_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F. Returns
%    the Hessian of the log likelihood with respect to PARAM. At
%    the moment PARAM can be only 'latent'. LLG2 is a vector with
%    diagonal elements of the Hessian matrix (off diagonals are
%    zero). This subfunction is needed when using Laplace
%    approximation or EP for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG, LIK_MULTINOM_LLG3, GPLA_E

% multinom:
  [n,nout]=size(y);
  N = sum(y,2)*ones(1,nout);
  f=reshape(f,n,nout);

  expf2 = exp(f);
  pi2 = expf2./(sum(expf2, 2)*ones(1,nout));
  pi_vec=pi2(:).*N(:);

  pi_mat=zeros(nout*n, n);
  for i1=1:nout
    pi_mat((1+(i1-1)*n):(nout*n+1):end)=pi2(:,i1).*sqrt(N(:,i1));
  end
  %     D = diag(pi_vec);
  %     llg2 = -D + pi_mat*pi_mat';

end

function [dw_mat] = lik_multinom_llg3(lik, y, f, param, z)
%LIK_MULTINOM_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_MULTINOM_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, class labels Y, and latent values F and
%    returns the third gradients of the log likelihood with
%    respect to PARAM. At the moment PARAM can be only 'latent'.
%    LLG3 is a vector with third gradients. This subfunction is
%    needed when using Laplace approximation for inference with
%    non-Gaussian likelihoods.
%
%  See also
%    LIK_MULTINOM_LL, LIK_MULTINOM_LLG, LIK_MULTINOM_LLG2, GPLA_E, GPLA_G

  [n,nout] = size(y);
  f2 = reshape(f,n,nout);

  N=sum(y, 2);
  expf2 = exp(f2);
  pi2 = expf2./(sum(expf2, 2)*ones(1,nout));
  pi_vec=pi2(:);

  dw_mat=zeros(nout,nout,nout,n);

  for cc3=1:nout
    for ii1=1:n

      pic=pi_vec(ii1:n:(nout*n));
      for cc1=1:nout
        for cc2=1:nout

          % multinom third derivatives
          cc_sum_tmp=0;
          if cc1==cc2 && cc1==cc3 && cc2==cc3
            cc_sum_tmp=cc_sum_tmp+pic(cc1);
          end
          if cc1==cc2
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc3);
          end
          if cc2==cc3
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc2);
          end
          if cc1==cc3
            cc_sum_tmp=cc_sum_tmp-pic(cc1)*pic(cc2);
          end
          cc_sum_tmp=cc_sum_tmp+2*pic(cc1)*pic(cc2)*pic(cc3);

          dw_mat(cc1,cc2,cc3,ii1)=cc_sum_tmp.*N(ii1);
        end
      end
    end
  end

end

function [logM_0, m_1, sigm2hati1] = lik_multinom_tiltedMoments(lik, y, i1, S2_i, M_i, z)
    %LIK_COXPH_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
    %
    %  Description
    %    [M_0, M_1, M2] = LIK_COXPH_TILTEDMOMENTS(LIK, Y, I, S2,
    %    MYY, Z) takes a likelihood structure LIK, class labels
    %    Y, index I and cavity variance S2 and
    %    mean MYY. Returns the zeroth moment M_0, mean M_1 and
    %    variance M_2 of the posterior marginal (see Rasmussen and
    %    Williams (2006): Gaussian processes for Machine Learning,
    %    page 55). This subfunction is needed when using EP for
    %    inference with non-Gaussian likelihoods.
    %
    %  See also
    %    GPEP_E

    error('tiltedMoment has not been implemented for multinom likelihood');

end

function [lpy, Ey, Vary] = lik_multinom_predy(lik, Ef, Varf, yt, zt)
%LIK_MULTINOM_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_MULTINOM_PREDY(LIK, EF, VARF YT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | y) = \int p(yt | f) p(f|y) df.
%    This requires also the incedence counts YT. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%    [LPY, EY, VARY] = LIK_MULTINOM_PREDY(LIK, EF, VARF, YT) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  N=sum(yt,2);
  S=10000;
  [ntest, nout]=size(yt);
  pi=zeros(ntest,nout);
  lpy=zeros(ntest,nout);
  Ey=zeros(ntest,nout);
  Vary=zeros(size(Varf));
  Ef=reshape(Ef(:),ntest,nout);
  [notused,notused,c] =size(Varf);
  if c>1
    mcmc=false;
  else
    mcmc=true;
    Varf=reshape(Varf(:), ntest, nout);
  end
  for i1=1:ntest
    if mcmc
      Sigm_tmp = (Varf(i1,:));
      f_star=bsxfun(@plus, Ef(i1,:), bsxfun(@times, sqrt(Sigm_tmp), ...
        randn(S,nout)));
    else
      Sigm_tmp=(Varf(:,:,i1)'+Varf(:,:,i1))./2;
      f_star=mvnrnd(Ef(i1,:), Sigm_tmp, S);
    end

    tmp = exp(f_star);
    tmp = tmp./(sum(tmp, 2)*ones(1,size(tmp,2)));

    if nargout > 1
        Ey(i1,:) = N(i1).*mean(tmp);
        for z1 = 1:nout;
          for z2 = 1:nout
            for z3=1:S
              Var_tmp(:,:,z3) = (diag(tmp(z3,:)) - tmp(z3,:)'*tmp(z3,:));
            end
            if mcmc
              Vary(i1+(0:nout-1)*ntest,:) = diag(N(i1).*mean(Var_tmp,3));
            else
              Vary(:,:,i1) = N(i1).*mean(Var_tmp,3);
            end
          end
        end
    end
    lpy=[];
    if ~isempty(yt)
      ytmp = repmat(yt(i1,:),S,1);
      lpy(i1,:) = log(mean( mnpdf(ytmp,tmp) ));
    end
  end
  lpy=lpy(:);
  Ey=Ey(:);
  Vary=Vary(:);
end

function p = lik_multinom_invlink(lik, f, z)
%LIK_MULTINOM_INVLINK Returns values of inverse link function
%
%  Description
%    P = LIK_MULTINOM_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_MULTINOM_LL, LIK_MULTINOM_PREDY
p = multinominv(f).*z;
end

function reclik = lik_multinom_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_MULTINOM_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    reclik.type = 'Multinom';
    reclik.nondiagW = true;

    % Set the function handles
    reclik.fh.pak = @lik_multinom_pak;
    reclik.fh.unpak = @lik_multinom_unpak;
    reclik.fh.ll = @lik_multinom_ll;
    reclik.fh.llg = @lik_multinom_llg;
    reclik.fh.llg2 = @lik_multinom_llg2;
    reclik.fh.llg3 = @lik_multinom_llg3;
    reclik.fh.tiltedMoments = @lik_multinom_tiltedMoments;
    reclik.fh.predy = @lik_multinom_predy;
    reclik.fh.invlink = @lik_multinom_invlink;
    reclik.fh.recappend = @lik_multinom_recappend;
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function Y = prtUtilY(varargin)
% prtUtilY     Quickly generate class labels for the PRT.
%
% Syntax: Y = prtUtilY(nH0,nH1,nH2,...)
%
% Inputs:
%   nH0 - The number of H0 samples
%   nH1 - The number of H1 samples
%       ...
%
% Outputs:
%   Y - A PRT compliant labeled vector
%
% Example:
%   Y = prtUtilY(100,100);
%   Y = prtUtilY([],100,100);
%   Y = prtUtilY(100,0,100);
%
% See also: prtDataGen*

Y = [];

if nargin == 1
    if length(varargin{1}) > 1

        temp = mat2cell(varargin{1}(:),ones(length(varargin{1}),1),1);
        Y = prtUtilY(temp{:});
    else
        Y = zeros(varargin{1},1);
    end
    return
end

for iInput = 1:nargin
    nHi = varargin{iInput};
    if ~isempty(nHi)
        if nHi > 0
            Y(end+1:end+nHi,1) = iInput-1;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [Y, xY] = spm_summarise(V,xY,fhandle,keepNaNs)
% Summarise data within a Region of Interest
% FORMAT [Y, xY] = spm_summarise(V,xY,fhandle)
% V       - [1 x n] vector of mapped image volumes to read (from spm_vol)
%           Or a char array of filenames
% xY      - VOI structure (from spm_ROI)
%           Or a VOI_*.mat (from spm_regions) or a mask image filename
%           Or the keyword 'all' to summarise all voxels in the images
%           Or a [3 x m] matrix of voxel coordinates {mm}
% fhandle - function handle to be applied on image data within VOI
%           Must transform a [1 x m] array into a [1 x p] array
%           Default is Identity (returns raw data, vectorised into rows).
%           Can also use keyword 'litres' to compute the total volume,
%           within the region of interest, for a tissue segment image.
%
% Y       - [n x p] data summary
% xY      - (updated) VOI structure
%__________________________________________________________________________
%
% Example:
% spm_summarise('beta_0001.nii',...
%               struct('def','sphere', 'spec',8, 'xyz',[10 20 30]'),...
%               @mean)
%__________________________________________________________________________
% Copyright (C) 2010-2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin, Ged Ridgway
% $Id: spm_summarise.m 7384 2018-07-31 13:36:15Z guillaume $

%-Argument checks
%--------------------------------------------------------------------------
if nargin < 1 || isempty(V)
    [V, sts] = spm_select([1 Inf], 'image', 'Specify Images');
    if ~sts, error('Must select 1 or more images'), end
end
if iscellstr(V), V = char(V); end
if ischar(V), V = spm_data_hdr_read(V); end
spm_check_orientations(V);

if nargin < 2 || isempty(xY), xY = struct; end
if ischar(xY)
    if strcmpi(xY, 'all')
        xY = struct('def', 'all');
    elseif any(regexpi(xY, '\.mat$'))
        try
            load(xY,'xY'); % VOI_*.mat file. Warns if .mat has no xY ...
            xY = rmfield(xY,'XYZmm'); % ...  error if .mat has no xY
        catch
            xY = struct; % GUI specification in spm_ROI
        end
    else % assume mask image filename
        xY = struct('def','mask', 'spec',xY);
    end
elseif isnumeric(xY) && any(size(xY, 1) == [3 4])
    xY = struct('XYZmm', xY(1:3, :));
elseif isstruct(xY) && isfield(xY,'fname')
    xY = struct('def','mask', 'spec',xY);
elseif ~isstruct(xY)
    error('Incorrect xY specified')
end
if ~isfield(xY,'XYZmm'), [xY, xY.XYZmm] = spm_ROI(xY,V(1)); end

if nargin < 3 || isempty(fhandle), fhandle = @(x) x; end
if ischar(fhandle) && strcmp(fhandle, 'litres')
    vsz     = abs(det(V(1).mat));  % voxel size in mm^3
    fhandle = @(x) sum(x) * vsz / 1e6;
end
if ischar(fhandle), fhandle = str2func(fhandle); end

% Undocumented option in case anyone wants to keep (e.g. to check for) NaNs
if nargin < 4, keepNaNs = false; end
if keepNaNs
    dropNaNs = @(x) x;
else
    dropNaNs = @(x) x(~isnan(x));
end

%-Summarise data
%--------------------------------------------------------------------------
XYZ = round(V(1).mat \ [xY.XYZmm; ones(1, size(xY.XYZmm, 2))]);

% Run on first volume to determine p, and transpose if column vector
Y = fhandle(dropNaNs(spm_data_read(V(1), 'xyz', XYZ)));
if ndims(Y) > 2
    error('Function must return a [1 x p] array')
elseif size(Y, 1) ~= 1
    if size(Y, 2) == 1
        Y = Y';
    else
        error('Function returned a [%d x %d] array instead of [1 x p]', ...
            size(Y, 1), size(Y, 2))
    end
end

% Preallocate space and then run on remaining volumes
Y(2:numel(V), :) = 0;
for i = 2:numel(V)
    Y(i, :) = fhandle(dropNaNs(spm_data_read(V(i), 'xyz', XYZ)));
end

author: spm


text:
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 5;

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

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d:', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12g', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [N,rho,P]=atmosParam4GasTemp(gasTable,T,rhow)
%%ATMOSPARAM4GASTEMP  Get basic parameters for atmospheric refractivity,
%                     density, and pressure given a table of constitudent
%                     gasses in the atmosphere, the temperature and the
%                     absolute humidity. The model is best suited for
%                     altitudes below 90km as ionozed parameters, such as
%                     anomolous oxygen, are not used.
%
%INPUTS: gasTable An NX2 cell array where gasTable{i,1} is a string
%                 describing the ith constituent atmospheric element and
%                 gasTable{i,2} is the number density of the element in
%                 particles per cubic meter. For a list of constituent
%                 elements that can be handled, see the documentation for
%                 the Constants.gasProp method. Unknown constituent
%                 elements that are passed will be ignored.
%               T The temperature in degrees Kelvin.
%            rhow An optional parameter specifying the mass density of
%                 water vapor at the point in question in units of
%                 kilograms per cubic meter. If omitted, the air is assumed
%                 to be dry (rhow=0). The total density of the air is
%                 assumed to be the sum of the dry air density and rhow.
%                 Alternatively, this parameter can be omitted and 'H2O'
%                 can be given as one of the constituent elements in
%                 gasTable.
%
%OUTPUTS: N The refractivity of the atmosphere. In this model, N is always
%           real. N=10^6*(n-1) where n is the index of refraction. This is
%           generally valid for frequencies from L-band (1GHz) to 10 GHz
%           (the middle of X-band).
%       rho The atmospheric density at the point in question in units of
%           kilograms per cubic meter.
%         P The atmospheric pressure at the point in question in units of
%           Newtons per square meter (Pascals). It assumes that the gasses
%           can be treated as ideal gasses.
%
%The refractive index is then found using the dry and wet air densities
%using the formula of [1], which should be valid for frequencies between
%1GHz and 10GHz. It ignores the lossiness of the atmosphere.
%
%REFERENCES:
%[1] J. M. Aparicio and S. Laroche, "An evaluation of the expression of
%    the atmospheric refractivity for GPS signals," Journal of Geophysical
%    Research, vol. 116, no. D11, 16 Jun. 2011.
%
%September 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3)
   rhow=0;
end

numGasses=size(gasTable,1);
numberDensities=[gasTable{:,2}];
totalGasNumberDensity=sum(numberDensities);

%rhod will be the total mass density of dry air at the point in kg/m^3.
rhod=0;
for curGas=1:numGasses
    AMU=Constants.gasProp(gasTable{curGas,1});
    %If an unknown gas is given, then ignore it.
    if(~isempty(AMU))
        %The number density has units of particles per cubic meter.
        %Multiplied by the atomic mass of the gas, we have atomic mass
        %units (AMU) per cubic meter. Multiplied by the value of the atomic
        %mass unit in kilograms, we get kilograms per cubic meter.
        rhod=rhod+Constants.atomicMassUnit*numberDensities(curGas)*AMU;
    else
        %This is so that the number density of the unknown constitutent is
        %also ignored when computing the pressure, so that the results are
        %consistent with the computation of the density.
        numberDensities(curGas)=0;
    end
end

rho=rhod+rhow;

%To use the ideal gas law to find the air pressure, the number of water
%molecules per cubic meter of the atmosphere is needed. This is obtained
%using the molar mass of water (H2O) and Avogadro's constant
Av=Constants.AvogadroConstant;
%The molar mass of water (H2O) in atomic mass units (grams per mole).
HAMU=Constants.elementAMU(1);
OAMU=Constants.elementAMU(8);
MMWater=HAMU*2+OAMU;

%The number of atoms of water per cubic meter. The 1000 transforms grams to
%kilograms.
NH2O=rhow/(1000*Av*MMWater);

%The total number density of the gasses in the atmosphere. That is, the
%number of atoms per cubic meter.
NTotal=totalGasNumberDensity+NH2O;
kB=Constants.BoltzmannConstant;
P=NTotal*kB*T;

%T is the temperature in Kelvin.
tau=273.15/T-1;
N0=(222.682+0.069*tau)*rhod+(6701.605+6385.886*tau)*rhow;
N=N0*(1+10^(-6)*N0/6);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function x_hat=phit_fp(y,m,n,ops)
% function x_hat=phit_fp(y,m,n,ops)
% PHIT_FP projects a length m signal onto the transpose of an mxn measurment matrix.
% Input:
%       y   : length m signals (for example a vectorized image)
%       m   : the number of measurements
%       n   : the signal length
%       ops : currently unused argument
%Output:
%       x   : result of projecting y onto phi'.
    x_hat=zeros(n,1);
    remaining_rows=n;
    K=4096;
    iters=ceil(n/K);
    for i=0:iters-1
        rng(i);
        row_num=min(K,remaining_rows);
        phi_columns=randn(m,row_num);
        column_norms=sqrt(sum(abs(phi_columns).^2,1));
        phi_columns=bsxfun(@rdivide,phi_columns,column_norms);
        phit_rows=phi_columns';
        x_hat(K*i+1:K*i+row_num)=phit_rows*y;
        remaining_rows=remaining_rows-K;
    end
end

author: ricedsp


text:
% Hessian of the measurement function in BOT-demo.

% Copyright (C) 2007 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function dY = bot_d2h_dx2(x,s)
  % Space for Hessians. Note that we need a Hessian for
  % each dimension in the measurement space, that is we need
  % a Hessian for each sensor in this case.
  dY = zeros(size(s,2),size(x,1),size(x,1));

  % Loop through sensors.
  for i=1:size(s,2)
    % Derivative twice wrt. x
    dx2 = -2*(x(1)-s(1,i)) / ((x(1)-s(1,i))^2+(x(2)-s(2,i))^2)^2;
    % Derivative twice wrt. y
    dy2 = -2*(x(2)-s(2,i)) / ((x(1)-s(1,i))^2+(x(2)-s(2,i))^2)^2;
    % Derivative wrt. x and y
    dxdy = ((x(2)-s(2,i))^2-(x(1)-s(1,i))^2) / ((x(1)-s(1,i))^2+(x(2)-s(2,i))^2)^2;
    dh = [dx2  dxdy 0 0;...
	  dxdy dy2  0 0;...
	  0    0    0 0;...
          0    0    0 0];
    dY(i,:,:) = dh;
  end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function [ output_meta ] = meta2sicd_aimidb( aimidb_struct )
%META2SICD_AIMIDB Converts metadata stored in the AIMIDB NITF TRE into a
% SICD metadata format
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if any(~isspace(aimidb_struct.COUNTRY)) % Could be blank
    output_meta.CollectionInfo.CountryCode=aimidb_struct.COUNTRY;
end
output_meta.ImageFormation.SegmentIdentifier=aimidb_struct.CURRENT_SEGMENT;
try
    output_meta.Timeline.CollectStart=datenum(aimidb_struct.ACQUISITION_DATE,'yyyymmddHHMMSS');
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [y,outside,leads] = spm_eeg_wrap_momfit_vbecd(P,M,U)
% A cost function/wrapper to sit between non-linear optimisation spm_nlsi_gn.m
% and dipole fit routine spm_cfg_eeg_momentfit.m
% FORMAT [y,outside,leads] = spm_eeg_wrap_momfit_vbecd(P,M,U)
% sens and vol structures should be passed in M, where
%   sens = M.Setup.forward.sens;
%   vol  = M.Setup.forward.vol;
% P contains a list of the free parameters (assuming all position
%   parameters come first (in triplets) followed by all moment paameters
%   (also in triplets)
% U is unused
% At the momnent reduces the rank of the MEG leadfield 2 dimensions.
% leads are the lead fields of the dipoles fit
%__________________________________________________________________________
% Copyright (C) 2020 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_eeg_wrap_momfit_vbecd.m 7764 2020-01-02 15:34:03Z spm $

x = U.u; % input , unused

sens = M.Setup.forward.sens;
vol  = M.Setup.forward.vol;

siunits = M.Setup.forward.siunits;
if ~siunits
    warning('Data not in SI units, scaling (and therefore priors) maybe wrong');
end
chanunits = M.Setup.forward.chanunits;

mompars = P; % one moment parameter per source
Ndips   = length(mompars);
allpos  = M.pos; % source positions
allori  = M.ori; % source orientations

if (size(allpos,1)~=size(allori,1))|| (size(allpos,1)~=length(mompars))
    error('There must be equal number of moments, positions and orientations');
end

Ndips = length(mompars);

allmom = allori.*repmat(spm_vec(mompars),1,3);

if ft_senstype(sens, 'meg')
    RANK = 2; % restricting rank of MEG data, could change this in future
else
    RANK = 3;
end

y = 0;
outside = 0;
leads = zeros(Ndips,3,numel(sens.label));
for i=1:Ndips

    pos = allpos(i,:);

    mom = allmom(i,:); %% in nAm

    if siunits
        [tmp] = ft_compute_leadfield(1e-3*pos, sens, vol, 'reducerank',RANK, 'dipoleunit', 'nA*m', 'chanunit', chanunits);

    else
        [tmp] = ft_compute_leadfield(pos, sens, vol, 'reducerank',RANK);

    end

    gmn = tmp;
    leads(i,:,:) = gmn';

    y = y + gmn * mom';

end % for i


text:
function [pHr, pHAdjustment] = realpH(pHa, temp, is)
% Apparent glass electrode pH is not the same as real pH for thermodynamic calculations.
%
% Given the experimental glass electrode measurement of pH, this function returns
% the pH to be used for thermodynamic calculations, `pHc = -log10[H+]`,
% by subtracting the effect of the ion atmosphere around H+ which
% reduces its activity coefficient below unity.
% See `p49 Alberty 2003`.
%
% USAGE:
%
%    [pHr, pHAdjustment] = realpH(pHa, temp, is)
%
% INPUTS:
%    pHa:             apparent pH, measured by glass electrode experimentally
%    temp:            experimentally measured temperature
%    is:              estimate of ionic strength
%
% OUTPUTS:
%    pHr:             real pH to be used for thermodynamic calculations
%    pHAdjustment:    adjustment to pH
%
% .. Author: -Ronan M.T. Fleming

gibbscoeff = 1.10708 - (1.54508*temp)/10^3 + (5.95584*temp^2)/10^6; % p48 Alberty

%Adjust pH using Extended Debye-Huckle equation
pHAdjustment = ( (gibbscoeff*(is^(0.5))) / (log(10)*(1+1.6*(is^(0.5)))) );

% fprintf('%s\t%f\n','pH adjustment',phAdjustment);
pHr = pHa - pHAdjustment;

author: opencobra

repo: cobratoolbox


text:
classdef prtClassLogisticDiscriminant < prtClass & prtClassBig
    % prtClassLogisticDiscriminant  Logistic Discriminant classifier
    %
    %    CLASSIFIER = prtClassLogisticDiscriminant returns a LogisticDiscriminant classifier
    %
    %    CLASSIFIER = prtClassLogisticDiscriminant(PROPERTY1, VALUE1, ...) constructs a
    %    prtClassLogisticDiscriminant object CLASSIFIER with properties as specified by
    %    PROPERTY/VALUE pairs.
    %
    %    A prtClassLogisticDiscriminant object inherits all properties from the abstract class
    %    prtClass. In addition is has the following properties:
    %
    %   wTolerance       - The convergance tolerance of the weights
    %   irlsStepSize     - Step size used in training. Can be set to a
    %                      double, or 'hessian'. If 'hessian', IRLS is
    %                      solved using the Hessian to estimate steps.
    %   maxIter          - maximum IRLS iterations
    %   nIterations      - number of iterations used, set during training
    %   wInitTechnique   - Technique to initialize weights, can be set to
    %                      'FLD', 'randn', and 'manual'
    %   manualInitialW   - The values the weights are initialized to if
    %                      wInitTechnique is set to 'manual'
    %   wTolerance       - Convergence tolerance on weight vector
    %   handleNonPosDefR - What to do when R is non-positive definte, can
    %                      be set to 'regularize' or 'exit'. When set to
    %                      regularize, the classifier will attempt to
    %                      regularize the matrix. When set to exit the
    %                      classifier will exit.
    %
    %   w                - The regression weights, estimated during training
    %                      w(1) corresponds to the DC bias and w(2:end)
    %                      corresponds to the weights for the features
    %
    %    For more information on LogisticDiscriminant classifiers, refer to the
    %    following URL:
    %
    %    http://en.wikipedia.org/wiki/Logistic_regression
    %
    %    A prtClassLogisticDiscriminant object inherits the TRAIN, RUN,
    %    CROSSVALIDATE and KFOLDS methods from prtAction. It also inherits
    %    the PLOT method from prtClass.
    %
    %    Example:
    %
    %     TestDataSet = prtDataGenUnimodal;           % Create some test and
    %     TrainingDataSet = prtDataGenUnimodal;       % training data
    %     classifier = prtClassLogisticDiscriminant;  % Create a classifier
    %     classifier = classifier.train(TrainingDataSet);       % Train
    %     classified = run(classifier, TestDataSet);            % Test
    %     subplot(2,1,1);
    %     classifier.plot;
    %     subplot(2,1,2);
    %     [pf,pd] = prtScoreRoc(classified,TestDataSet);
    %     h = plot(pf,pd,'linewidth',3);
    %     title('ROC'); xlabel('Pf'); ylabel('Pd');
    %
    %    See also prtClass, prtClassLogisticDiscriminant, prtClassBagging,
    %    prtClassMap, prtClassCap, prtClassBinaryToMaryOneVsAll,
    %    prtClassDlrt, prtClassPlsda, prtClassFld, prtClassRvm,
    %    prtClassLogisticDiscriminant,  prtClass

    properties (SetAccess=private)
        name = 'Logistic Discriminant'  % Logistic Discriminant

        nameAbbreviation = 'LogDisc'  % LogDisc

        isNativeMary = false;  % True
    end

    % properties (SetAccess = protected)
    properties
        % w
        %   w is a DataSet.nDimensions + 1 x 1 vector of projection weights
        %   learned during LogDisc.train(DataSet)

        w = [];  % Regression weights

        % nIterations
        %   Number of iterations used in training.  This is set to a number
        %   between 1 and maxIter during training.

        nIterations = nan;  % The number of iterations used in training

        wPriorSigma = nan;  % Prior standard deviation on weight vector, use NAN to implement maximum likelihood estimates
    end
    properties
        % irlsStepSize
        %   irlsStepSize can be the string 'hessian', or a double value
        %   (typically << 1).  If irlsStepSize is 'hessian', IRLS is solved
        %   using the Hessian to estimate steps.  This can be numerically
        %   unstable.  Otherwise, training takes steps in the direction of
        %   the gradient with a step size of irlsStepSize*gradient.  Default
        %   value is 0.05.

        irlsStepSize = .05; % The stepsize

        % maxIter
        %   Maximum number of iterations to allow before exiting without
        %   convergence.

        maxIter = 500;  % Maxmimuum number of iterations

        % wInitTechnique
        %   wInitTechnique specifies how training should initialize the w
        %   vector.  Possible values are 'FLD', 'randn', and 'manual'.

        wInitTechnique = 'FLD';  % Weight initialization technique

        % manualInitialW
        %   manualInitialW is used to set the initial w value when
        %   wInitTechnique is 'manual'.  manualInitialW must be a vector
        %   and of length(TrainDataSet.nDimensions + 1)

        manualInitialW = []; % The value of the initial weights if initialized manually

        % wTolerance
        %   Convergence tolerance on w.  When norm(w - wOld) < wTolerance,
        %   convergence is reached, and training exits.

        wTolerance = 1e-2;  % The convergance tolerance of the weights

        % handleNonPosDefR
        %   It is possible to obtain non-positive definite re-weighted
        %   least-squares matrices in the optimization process.  Often this
        %   indicates convergence.  handleNonPosDefR can be one of 'exit'
        %   or 'regularize'.  'exit' tells training to exit with the
        %   current weight vector, and 'regularize' attempts to
        %   diagonal-load the R matrix to acheive a well conditioned
        %   matrix.  Often regularizing is a losing battle.

        handleNonPosDefR = 'exit';  % The action taken when R is non-positive definite

        sgdPassesThroughTheData = 10;
        sgdLearningRate = @(t)((sqrt(t) + 10).^(-0.9));
        sgdRegularization = 0.1;
        sgdWeightTolerance = 1e-6;

    end

    methods

        function Obj = prtClassLogisticDiscriminant(varargin)
            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end

        function Obj = set.irlsStepSize(Obj,val)
            if ~prtUtilIsPositiveScalar(val) && ~strcmpi(val,'hessian')
                error('prt:prtClassLogisticDiscriminant:irlsStepSize','irlsStepSize must be a positive scalar (or the string ''hessian'')');
            end
            Obj.irlsStepSize = val;
        end

        function Obj = set.maxIter(Obj,val)
            if ~prtUtilIsPositiveScalarInteger(val)
                error('prt:prtClassLogisticDiscriminant:maxIter','maxIter must be a positive scalar integer');
            end
            Obj.maxIter = val;
        end

        function Obj = set.wInitTechnique(Obj,val)
            if ~isa(val,'char') || ~any(strcmpi(val,{'fld','randn','manual'}))
                error('prt:prtClassLogisticDiscriminant:wInitTechnique','wInitTechnique must be a string matching one of ''fld'', ''randn'', or ''manual''');
            end
            Obj.wInitTechnique = val;
        end

        function Obj = set.manualInitialW(Obj,val)
            if ~isnumeric(val) && ~isvector(val) && ~isempty(val)
                error('prt:prtClassLogisticDiscriminant:manualInitialW','manualInitialW must be a numeric vector');
            end
            Obj.manualInitialW = val;
        end

        function Obj = set.wTolerance(Obj,val)
            if ~prtUtilIsPositiveScalar(val)
                error('prt:prtClassLogisticDiscriminant:wTolerance','wTolerance must be a positive scalar');
            end
            Obj.wTolerance = val;
        end

        function Obj = set.handleNonPosDefR(Obj,val)
            if ~isa(val,'char') || ~any(strcmpi(val,{'exit','regularize'}))
                error('prt:prtClassLogisticDiscriminant:handleNonPosDefR','handleNonPosDefR must be a string matching one of ''exit'' or ''regularize''');
            end
            Obj.handleNonPosDefR = val;
        end

    end

    methods (Access=protected, Hidden = true)

        function Obj = trainAction(Obj,DataSet)
            %Obj = trainAction(Obj,DataSet)

            if ~DataSet.isBinary
                error('prtClassLogisticDiscriminant:nonBinaryTraining','Input dataSet for prtClassLogisticDiscriminant.train must be binary');
            end

            %Helper functions:
            sigmaFn = @(x) 1./(1 + exp(-x));
            rCondDiag = @(vec) min(vec)/max(vec);

            switch lower(Obj.wInitTechnique)
                case 'fld'
                    Fld = prtClassFld;
                    Fld = Fld.train(DataSet);
                    Obj.w = [1;Fld.w]; %append zero for offset
                case 'randn'
                    Obj.w = randn(DataSet.nFeatures+1,1);
                case 'manual'
                    assert(isvector(Obj.manualInitialW) & numel(Obj.manualInitialW) == DataSet.nFeatures + 1,'manualInitialW must be a vector and have %d elements',DataSet.nFeatures + 1);
                    assert(isequal(size(Obj.manualInitialW), [DataSet.nFeatures + 1,1]), 'manualInitialW must be a %d x 1 vector',DataSet.nFeatures + 1);
                    Obj.w = Obj.manualInitialW;
                otherwise
                    error('Invalid value for Options.wInitTechnique; wInitTechnique must be one of {''FLD'',''randn'',''manual''}');
            end

            DataSet = DataSet.retainLabeled; % Only use labeled data to train

            x = DataSet.getObservations;
            x = cat(2,ones(size(x,1),1),x);  %append "ones"
            %             y = DataSet.getTargets;
            %             if ~isequal(DataSet.uniqueClasses,[0;1])
            %                 bm = DataSet.getTargetsAsBinaryMatrix;
            %                 y(logical(bm(:,1))) = 0;
            %                 y(logical(bm(:,2))) = 1;
            %             end
            y = DataSet.getBinaryTargetsAsZeroOne;

            yOut = sigmaFn((x*Obj.w)')';
            rVec = yOut.*(1-yOut);

            converged = 0;
            nIter = 1;

            while ~converged
                Obj.nIterations = nIter;

                %these next lines are for numerical instabilities; these are detected
                %using the rCond of rVec or any nans in rVec
                if rCondDiag(rVec) < eps*2 || any(isnan(rVec))
                    %Numerical instability options are normalize, exit, or stepsize.
                    switch lower(Obj.handleNonPosDefR)
                        case 'regularize'
                            warning('prt:generateLogDisc:stepSize','rcond(R) < eps; attempting to diagonally load R');
                            diagAdd = 1e-5*max(rVec);
                            while(rCondDiag(rVec) < eps*2)
                                rVec = rVec + diagAdd;
                            end
                        case 'exit'
                            warning('prt:generateLogDisc:stepSize','rcond(R) < eps; Exiting; Try reducing classifier.irlsStepSize');
                            return;
                        otherwise
                            error('Invalid Obj.handleNonPosDefR field');
                    end
                end
                if isnan(Obj.wPriorSigma)
                    priorSigma = inf;
                else
                    priorSigma = Obj.wPriorSigma;
                end
                if isa(Obj.irlsStepSize,'char') && strcmpi(Obj.irlsStepSize,'hessian')
                    % Bishop equations:
                    %         wOld = w;
                    %         z = x*wOld - (R^-1)*(yOut - y);
                    %         w = (x'*R*x)^-1*x'*R*z;
                    %         %re-calculate yOut and R
                    %         yOut = sigmaFn((x*w)')';
                    %         R = diag(yOut.*(1-yOut));

                    %Non-matrix R Bishop equations (memory saving equations):
                    wOld = Obj.w;
                    z = x*wOld - bsxfun(@times,rVec.^-1,(yOut - y));
                    Obj.w = (x'*bsxfun(@times,rVec,x) + eye(size(x,2))*1./priorSigma)^-1*x'*bsxfun(@times,rVec,z);
                    %re-calculate yOut and R
                    yOut = sigmaFn((x*Obj.w)')';
                    rVec = yOut.*(1-yOut);
                elseif isa(Obj.irlsStepSize,'char')
                    error('String irlsStepSize must be ''hessian''');
                elseif isa(Obj.irlsStepSize,'double')
                    % Raw update equations:
                    %         wOld = w;
                    %         H = x'*R*x;
                    %         dE = x'*(yOut - y);
                    %         w = w - H^-1*dE*Options.irlsStepSize;  %move by stepsize * the hessian
                    %         %re-calculate yOut and R
                    %         yOut = sigmaFn((x*w)')';
                    %         R = diag(yOut.*(1-yOut));

                    %Non-matrix R Raw equations (memory saving equations):
                    wOld = Obj.w;
                    H = x'*bsxfun(@times,rVec,x) + eye(size(x,2))*1./priorSigma;

                    if rcond(H) < eps
                        warning('prt:generateLogDisc:stepSize','rcond(H) < eps; Exiting; Try reducing Options.irlsStepSize');
                        return;
                    end

                    dE = x'*(yOut - y);
                    Obj.w = Obj.w - H^-1*dE*Obj.irlsStepSize;  %move by stepsize * the hessian
                    %re-calculate yOut and R
                    yOut = sigmaFn((x*Obj.w)')';
                    rVec = yOut.*(1-yOut);
                end

                if norm(Obj.w - wOld) < Obj.wTolerance
                    converged = 1;
                end

                if nIter > Obj.maxIter
                    warning('prt:generateLogDisc:maxIter',sprintf('nIterations (%d) > maxIterations; exiting',nIter)); %#ok
                    return;
                end
                nIter = nIter + 1;
            end
        end

        function self = trainActionBig(self,ds)

            nTriesForNonEmptyBlock = 1000;
            useVariableLearningRate = isa(self.sgdLearningRate,'function_handle');
            converged = false;

            nMaxIterations = self.sgdPassesThroughTheData*ds.getNumBlocks();

            for iter = 1:nMaxIterations

                % Try to load a block
                for blockLoadTry = 1:nTriesForNonEmptyBlock
                    cBlockDs = ds.getRandomBlock;
                    if ~isempty(cBlockDs)
                        break
                    end
                end
                if blockLoadTry == nTriesForNonEmptyBlock
                    warning('Exiting after %d empty blocks were consecutlivly found.', nTriesForNonEmptyBlock);
                    break
                end

                % cBlockDs is now our dataset

                if iter==1
                    cW = zeros(cBlockDs.nFeatures+1,1);
                    wOld = inf;

                end
                cX = cat(2,ones(cBlockDs.nObservations,1), cBlockDs.X);

                % Get Y as -1, 1
                cY = cBlockDs.getTargetsAsBinaryMatrix;
                cY = cY(:,end); % In case of m-ary data we also do the last one.
                cY(cY == 0) = -1;

                %cW = cW + self.irlsStepSize*sum(bsxfun(@times,cY./(1 + exp(cY.*(cX*cW))),cX),1)';

                if useVariableLearningRate
                    cLearningRate = self.sgdLearningRate(iter);
                else
                    cLearningRate = self.sgdLearningRate;
                end
                cStep = (self.sgdRegularization*cW - mean(bsxfun(@times,cY./(1 + exp(cY.*(cX*cW))),cX),1)');
                cW = cW - cLearningRate*cStep;

                %cW = cW - cLearningRate*G^(-1/2)*cStep;
                cChange = norm(cW - wOld)/norm(cW);
                if cChange < self.sgdWeightTolerance
                    converged = true;
                    break
                end

                wOld = cW;

                plotOnIter = 1;
                if plotOnIter && ~mod(iter,plotOnIter)
                    plot(cW)
                    title(sprintf('iteration=%d - change=%0.3f',iter,cChange));
                    drawnow;
                end
            end

            if ~converged
                warning('prt:prtClassLogisticDiscriminant:notConverged','Convergence was not reached in the maximum number of allowed iterations.');
            end

            self.w = cW;
        end

        function ds = runAction(self,ds)
            ds.X = runFast(self,ds.X);
        end

        function y = runActionFast(self, x)
            x = cat(2,ones(size(x,1),1),x);
            y = 1./(1 + exp(-((x*self.w)')'));
        end
    end
end


text:
function [maps,stats] = SurveyFiringMaps(units,varargin)

%SurveyFiringMaps - Compute and plot firing maps for all subsessions.
%
%  USAGE
%
%    [maps,stats] = SurveyFiringMaps(units,<options>)
%
%    units          optional list of units, i.e. [electrode group, cluster]
%                   pairs; set cluster to -1 to process all clusters
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'minv'        minimum instantaneous velocity (default = 0)
%     'pixel'       size of the video pixel in cm (no default value)
%     'show'        set to 'off' to compute but not plot data (default = 'on')
%    =========================================================================
%
%  OUTPUT
%
%    The outputs are the same as for MapStats, except for map.z which is
%    replaced by map.rate.
%
%  SEE
%
%    See also FiringMap, PlotColorMap.

% Copyright (C) 2009-2013 by Anne Cei, Michaël Zugaro
%
% 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.

% Default values
minv = [];
pixel = [];
show = 'on';

% No unit list provided?
if ~(nargin >= 1 && isimatrix(units)),
	units = GetUnits;
	varargin = {units,varargin{:}};
else
	all = units(:,2) == -1;
	groups = unique(units(all,1));
	units(all,:) = [];
	units = unique([units;GetUnits(groups)],'rows');
end
nUnits = size(units,1);

if mod(length(varargin),2) ~= 0,
	error('Incorrect number of parameters (type ''help SurveyFiringMaps'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+1) ' is not a property (type ''help SurveyFiringMaps'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'minv',
			minv = varargin{i+1};
			if ~isdscalar(minv,'>=0'),
				error('Incorrect value for property ''minv'' (type ''help SurveyFiringMaps'' for details).');
			end
		case 'pixel',
			pixel = varargin{i+1};
			if ~isdscalar(pixel,'>0'),
				error('Incorrect value for property ''pixel'' (type ''help SurveyFiringMaps'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring_FMAT(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help SurveyFiringMaps'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SurveyFiringMaps'' for details).']);
	end
end

% Get positions
positions = GetPositions;
if isempty(positions), warning('No positions found for current subsession'); return; end
if ~isempty(minv),
	if isempty(pixel),
		error(['Missing pixel size for minimum velocity (type ''help SurveyFiringMaps'' for details).']);
	end
	x = GetPositions('coordinates','real','pixel',pixel);
	v = LinearVelocity(x,30);
	[~,in] = Threshold(v,'>',minv,'min',10,'max',2);
	positions = positions(in,:);
end

% Get start/stop events for each subsession
start = GetEvents('beginning of .*');
stop = GetEvents('end of .*');
nSubsessions = length(start);

if strcmp(show,'on'),
	titles = GetEventTypes('beginning of .*');
	status = Hide('status');
	Hide('on');
	figureList = [];
end

try
	n = 0;
	% Loop through units
	for j = 1:nUnits,
		if strcmp(show,'on'), figureList = [figureList figure]; end
		spikes = GetSpikes(units(j,:));
		% Loop through subsessions
		for i = 1:nSubsessions,
			p = Restrict(positions,[start(i) stop(i)]);
			s = Restrict(spikes,[start(i) stop(i)]);
			if nargin == 0,
				map = FiringMap(p(:,1:3),s,'nbins',[250 250],'smooth',5,'mintime',0);
			else
				[map,st] = FiringMap(p(:,1:3),s,'nbins',[250 250],'smooth',5,'mintime',0);
				st.unit = units(j,:);
				st.subsession = i;
				n = n + 1;
				stats(n) = st;
				maps(n) = map;
			end
			if strcmp(show,'on'),
				SquareSubplot(nSubsessions,i);
				PlotColorMap(map.rate,map.time,'ydir','reverse','bar','off','cutoffs',[0 10],'gamma',2);
				title = titles{i}(14:end);
				SplitTitle([title ' (' int2str(units(j,1)) '-' int2str(units(j,2)) ')'],round(150/sqrt(nSubsessions)));
			end
		end
	end
catch err
	if strcmp(show,'on'),
		Hide(status);
		Hide(figureList,status);
	end
	err.rethrow;
end

if strcmp(show,'on'),
	Hide(status);
	Hide(figureList,status);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function aboxes_out = merge_detected_bboxes(aboxes)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

num_methods  = length(aboxes);
num_classes  = length(aboxes{1});
num_imgs     = length(aboxes{1}{1});

aboxes_out   = cell(num_classes, 1);

for j = 1:num_classes, aboxes_out{j} = cell(num_imgs,1); end

for i = 1:num_imgs
    for j = 1:num_classes
        aboxes_this_img_this_cls = cell(num_methods, 1);
        for d = 1:num_methods
            aboxes_this_img_this_cls{d} = aboxes{d}{j}{i};
        end
        aboxes_out{j}{i} = cell2mat(aboxes_this_img_this_cls);
    end
end

aboxes_out = {aboxes_out};

end

author: gidariss

repo: mrcnn-object-detection


text:
%% HELLO_FSU uses the fsuClusterMatlab command to run the HELLO "computation".
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    '~/matlab/hello_parfor' allows us to specify an output directory;
%
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%
%    'p" indicates this is a "parfor" job;
%
%    'w' indicates that the MATLAB session should wait for completion;
%
%    workers is the number of processors requested;
%
%    @hello_fun indicates the function we want to execute;
%
%    { n } contains input arguments to the function.
%
%    The function "hello_fun" must correspond to a MATLAB M-file "hello_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed
%    by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2011
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  workers = 4;

  tic

  fsuClusterMatlab ( '~/matlab/hello_parfor', [], 'p', 'w', workers, ...
    @hello_fun, { n } )

  toc


text:
function [v2gMap, sqDist] = mrmMapVerticesToGray(vertexCoords, grayNodes, mmPerVox, grayEdges,distThresh)
%
% [v2gMap sqDist] = mrmMapVerticesToGray(vertices, grayNodes, mmPerVox, [grayEdges]);
%
% Finds a map between the mesh vertices and gray nodes.
% To see the coordinates of the gray matter node (or nodes) nearest to
% mesh vertex I, we can use:
%
%      nearestNodeCoords = grayNodes(1:3, v2gMap(:,I))
%
%
% If the grayEdges are also passed in, then we'll compute the mapping for
% all gray layers. Layer 1 mapping will be in the first row, and layers
% >1 in subsequent columns. As always, vertices with no mapping will be set
% to 0.
%
% The vertexGrayMap now has multiple rows to represent mappings from
% each vertex (the columns) to all gray layers. This mapping is computed in
% two steps:
%
% 1. For each vertex, find the nearest layer 1 gray node, where 'nearest'
% is simply the smallest 3d Euclidian distance, with a threshold to leave
% vertices that are too far from any layer 1 nodes unmapped (zero entry in
% vertexGrayMap). This threshold is currently set to 4mm- edit
% mrmMapVerticesToGray if you want to change this. This layer-1 mapping
% forms the first column of the vertexGrayMap.
%
% 2. For each vertex (row), find all the nodes from layers >1 that should
% map to this vertex. This is a somewhat ill-posed problem, so we had to
% make a choice about which solution we wanted. We decided to use the gray
% layer connections to find all the layer >1 nodes that form a direct
% connection to the nearest layer 1 node that we mapped in step 1. As you
% discovered, this gives us a very redundant sampling of the gray nodes.
% The number of nodes >1 that are most directly connected to a given layer
% 1 node can be very high, especially in regions of high convexity. The
% algorithm for doing this essentially starts at the highest layers and
% maps them down to the most direct (ie. shortest path) layer 1 node. In a
% region of GM at the end of a thick finger of WM, this can map many nodes
% to one layer 1 node (apparently 84 in one case). The vertexGrayMap is an
% array, so the number of columns represents the most densely connected
% case. If you look at it, there are probably only a few vertices with as
% many as 84 connections. Most vertices probably have lots of zeros in
% those higher columns. In fact, in regions of high concavity, you should
% see some with no layers >1 (ie. a row with all zeros except for the first
% column).
%
% We realize that this choice of implementation might not be what everyone
% wants. But it works pretty well for our purposes. Of course, one downside
% to the highly redundant mapping is an effective 'smoothing' of the data,
% especially when you choose to average all layers (in your 3d window
% preferences settings). What we do is look at our data collapsed across
% layers as well as combining all layers, just to make sure we aren't
% missing something.

%
% HISTORY:
%  2003.09.16 RFD (bob@white.stanford.edu)
%  2004.05.21 RFD: now uses Dan Merget's nearpoints to find a proper
%  mapping. No more hacks! Aside from being a bit faster, the result should
%  be substantially more accurate.
%  2005.07.26 RFD: We now return a mapping for all gray layers, if
%  grayEdges is passed in.
%  2005.08.05 RFD: all gray layer mapping is turned off by default, until
%  we figure out how to make it go faster.
%  2005.10.26 GB: Tentative speed improvement in the main loop
%  2006.06.01 RFD: added comments above. These were prompted by an email
%  exchange with David Ress.
%  2007.05.23 ARW : Gray layer mapping seems to be available. Which is
%  nice. Added some minor mods to return sqDist and detect dist thresholds.
verbose = prefsVerboseCheck;

if notDefined('mmPerVox'), error('Voxel size (mm per vox) is required.'); end;
if (notDefined('distThresh'))
    if (ispref('VISTA','defaultSurfaceWMMapDist'))
		if verbose,
	        disp('Setting distThresh to the one in VISTA preferences');
		end
        distThresh = getpref('VISTA','defaultSurfaceWMMapDist');
    else
       if verbose, disp('Setting distThresh to 2'); end
       distThresh = 3;
    end
end

vertexCoords = double(vertexCoords);
grayNodes = double(grayNodes);

if notDefined('grayEdges')
    grayEdges=[];
end

grayEdges = double(grayEdges);

prefs = mrmPreferences;
if(strcmp(prefs.layerMapMode,'all'))
    mapToAllLayers = true;
else
    mapToAllLayers = false;
end

% This is now a real distance threshold, in mm. For each mesh vertex, the
% algorithm will find the nearest gray coord. If the nearest coord is >
% distThresh from the vertex, then that vertex gets no mapping ('0').

% The gray coordinates are in voxels in the vAnatomy file.  This scales
% them into real physical (mm) coordinates.  And transposes them.
grayCoords = grayNodes([1,2,3], :);
grayCoords = [grayCoords(1,:).*mmPerVox(1); ...
        grayCoords(2,:).*mmPerVox(2); ...
        grayCoords(3,:).*mmPerVox(3) ]';

% Transposes these mesh coordinates, which were already built in real
% physical coordinates.
% Major comments needed here.
vertexCoords = vertexCoords' + 1;

% Mask out non-layer 1 nodes so that they are not found
if (~strcmp(prefs.layerMapMode,'any'))
	if verbose,     disp('Masking out layers > 1');   end
    grayCoords(grayNodes(6,:)~=1,:) = -9999;
end

[v2gMap, sqDist] = nearpoints(double(vertexCoords'), double(grayCoords'));

if verbose,
	fprintf('Excluding mesh nodes further than %d away from the boundary \n', ...
			distThresh);
end

v2gMap(sqDist > (distThresh^2)) = 0;

v2gMap = int32(v2gMap);

if(mapToAllLayers && exist('grayEdges','var') && ~isempty(grayEdges))
    layer1Nodes = v2gMap;
    curValid = double(layer1Nodes) > 0;
    numLayers = max(grayNodes(6,:));
    n = length(curValid);
    curLayerNum = 2;

    % GB 2005.10.26
    % We map from higher layers to lower layers,
    % selecting the closest connected node in the lower layer. This way we
    % will avoid multiply-connected

    noNeighborsDown = 0;
    noNeighbor = 0;
    noNeighbor2 = 0;

    v2gMapTemp = zeros(numLayers,size(grayNodes,2));

    h = mrvWaitbar(0, 'Computing between layers connections...');
    progress = 0;
    for iLayer = numLayers:-1:curLayerNum
        midprogress = 0;

        curNodes = find(grayNodes(6,:) == iLayer);
        for iNode = curNodes
            % disp(iNode) % debugging
            offset = grayNodes(5,iNode);
            numConnected = grayNodes(4,iNode);
            if numConnected == 0
                noNeighbor = noNeighbor + 1;
                continue %mrmMapVerticesToGray
            end
            neighbors = grayEdges(offset:offset + numConnected - 1);
            neighbors(grayNodes(6,neighbors) ~= (iLayer - 1)) = [];
            distance = sum(grayNodes(1:3,neighbors).^2,1);
            [value,nearest] = min(distance);
            if ~isempty(nearest)
                v2gMapTemp(iLayer,iNode) = neighbors(nearest);
            elseif length(nearest) > 1
                neighbors = neighbors(nearest);
                distanceIndex = abs(neighbors - iNode);
                [value,nearest] = min(distanceIndex);
                v2gMapTemp(iLayer,iNode) = neighbors(nearest);
            else
                neighbors = grayEdges(offset:offset + numConnected - 1);
                neighbors(grayNodes(6,neighbors) < iLayer - 1) = [];
                nextNeighbors = [];
                neighborsDown = [];
                iter = 0; % debugging
                while isempty(neighborsDown) && ~isequal(neighbors,nextNeighbors)
                    iter = iter+1; %debugging
                    if iter> 1000,
                       disp foo
                       break
                    end
                    neighbors = union(nextNeighbors,neighbors);
                    nextNeighbors = [];
                    for iNeighbor = 1:length(neighbors)
                        offset = grayNodes(5,neighbors(iNeighbor));
                        numConnected = grayNodes(4,neighbors(iNeighbor));
                        nextNeighbors = [nextNeighbors grayEdges(offset:offset+numConnected-1)];
                    end
                    nextNeighbors = unique(nextNeighbors);
                    neighborsDown = nextNeighbors(grayNodes(6,nextNeighbors) == (iLayer - 1));
                end

                if isempty(neighborsDown)
                    noNeighbor2 = noNeighbor2 + 1;
                else
                    distance = sum(grayNodes(1:3,neighborsDown).^2,1);
                    [value,nearest] = min(distance);
                    if length(nearest) == 1
                        v2gMapTemp(iLayer,iNode) = neighborsDown(nearest(1));
                    else
                        neighborsDown = neighborsDown(nearest);
                        distanceIndex = abs(neighborsDown - iNode);
                        [value,nearest] = min(distanceIndex);
                        v2gMapTemp(iLayer,iNode) = neighborsDown(nearest);
                    end
                end
                noNeighborsDown = noNeighborsDown + 1;
            end
        end
        progress = progress + 1/(numLayers - curLayerNum + 1);
        mrvWaitbar(progress/2,h);
    end

    indices = 2*ones(1,size(v2gMap,2));
    indexNodes = cell(1,size(grayNodes,2));
    for curNode = find(curValid)
        indexNodes{v2gMap(1,curNode)} = [indexNodes{v2gMap(1,curNode)} curNode];
    end

    h = mrvWaitbar(1/2, h, 'Creating connections table...');
    progress = 0;
    for iNode = 1:size(grayNodes,2)
        if iNode > (progress + 1/10)*size(grayNodes,2)
            progress = progress + 1/10;
            mrvWaitbar(0.5 + progress/2);
        end

        curLayer = grayNodes(6,iNode);
        curNode = iNode;
        for iLayer = curLayer:-1:2
            if curNode == 0
                break
            end
            curNode = v2gMapTemp(iLayer,curNode);
        end
        if (curNode == 0) || (curNode == iNode)
            continue
        end

        indexNode = indexNodes{curNode};
        if isempty(indexNode)
            continue
        end

        v2gMap(indices(indexNode),indexNode) = iNode;
        indices(indexNode) = indices(indexNode) + 1;

    end
    close(h);

end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [y] = gt(x, val)

y = cell(size(x));
for k = 1:numel(x)
  y{k} = x{k}>val;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [A, b, x, ProbInfo] = PRtomo(varargin)
% IRtomo Generates data for X-ray tomographic reconstruction problems
%
% [A, b, x, ProbInfo] = PRtomo
% [A, b, x, ProbInfo] = PRtomo(n)
% [A, b, x, ProbInfo] = PRtomo(n, options)
% [A, b, x, ProbInfo] = PRtomo(options)
%
% This function uses the "line model" to create a 2D X-ray tomography test
% problem with an N-times-N pixel domain, using p rays for each angle in
% the vector theta.
%
% The severity of the problem is almost unaffected by the problem parameters.
%
% Input:
%  n : size of the image; a scalar such that the size is n x n.
%      Default: n = 256.
%  options : Structure containing the following optional fields:
%    phantomImage - user supplied test image of size n x n, of type numeric,
%             2-D only, or character string indicating
%             'shepplogan' : Shepp-Logan phantom
%             'smooth' : a smooth image
%             'binary' : a binary image
%             'threephases' : random image with pixel values 0, 0.5, 1
%                arranged in domains
%             'threephasessmooth' : similar to threephases, but the
%                domains have smoothly varying pixel values and there is
%                a smooth background
%             'fourphases' : similar to 'binary' but with three phases
%                separated by (thin) structures that form the fourth phase
%             'grains' : a random image with Voronoi cells
%             'ppower' : a random image with patterns of nonzero pixels
%             Default: 'shepplogan'.
%             This image is then stored in the output vector x.
%    CTtype - string that defines the type of CT problem:
%             'parallel'  : parallel beam geometry (default),
%             'fancurved' : fan beam with curved detector.
%    sm     - logical; if true (default )then A is a sparse matrix,
%             otherwise it is a function handle.
%    angles - vector of projection angles, in degrees.
%             Default: 0:1:179 for parallel beam, 0:2:358 for fan beam.
%    p      - number of rays for each source angle.
%             Default: p = round(sqrt(2)*n).
%    d      - Parallel beam only: scalar denoting the distance from the
%             first ray to the last; default: d = p-1.
%    R      - Fan beam only: the distance from the source to the center
%             of the phantom is R*N.  Default: R = 2.
%    span   - Fan beam only: scalar that determines the angular span of
%             the rays, in degrees. The default value is defined such
%             that from the source at (0,R*N) the first and last rays
%             hit the image corners (-n/2,n/2) and (n/2,n/2).
%
% Output:
%  A : Sparse matrix or function handle for forward/adjoint problem.
%  b : Vector with projection data (the sinogram).
%  x : Vector with image (i.e., exact image with stacked columns).
%  ProbInfo : structure containing some information about problem
%      problemType - kind of test problem (in this case: 'tomography')
%      xType       - solution type (in this case 'image2D')
%      bType       - data type (in this case 'image2D')
%      xSize       - size of image x
%      bSize       - size of sinogram b
%
% See also: PRblur, PRdiffusion, PRinvinterp2, PRnmr, PRseismic,
% PRspherical, PRtomo, PRnoise, PRshowb, PRshowx, radon

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD Licence. A separate license file should be provided as part
% of the package.

% Initialization: set the default image size.
default_N = 256;

% Check for acceptable number of optional input arguments.
switch length(varargin)
    case 0
        N = default_N; options = [];
    case 1
        if isa(varargin{1}, 'double')
            N = varargin{1}; options = [];
        else
            N = default_N; options = varargin{1};
        end
    case 2
        if isa(varargin{1}, 'double')
            N = varargin{1}; options = varargin{2};
        else
            N = varargin{2}; options = varargin{1};
        end
    otherwise
        error('Too many input parameters')
end

% Set default values for options.
defaultopt = struct('phantomImage', 'shepplogan', 'CTtype', 'parallel', ...
    'sm', true, 'angles', 0:1:179, 'p', NaN, 'R', 2, 'd', NaN, 'span', NaN);

% If input is 'defaults,' return the default options in A.
if nargin == 1 && nargout <= 1 && strcmp(varargin,'defaults')
    A = defaultopt;
    return;
end

if isempty(options)
    options = defaultopt;
end

options = PRset(defaultopt, options);

phantom    = PRget(options, 'phantomImage', [], 'fast');
CTtype     = PRget(options, 'CTtype',       [], 'fast');
angles     = PRget(options, 'angles',       [], 'fast');
sm         = PRget(options, 'sm',           [], 'fast');
p          = PRget(options, 'p',            [], 'fast');
R          = PRget(options, 'R',            [], 'fast');
d          = PRget(options, 'd',            [], 'fast');
span       = PRget(options, 'span',         [], 'fast');

% If a phantom image is given, then this defines N and the other parametes.
if isnumeric(phantom)
    % Make sure user input image is a matrix
    if ~ismatrix(phantom)
        error('Expected user supplied phantom image to be a 2-D array of type numeric')
    else
        N = size(phantom,1);
        if size(phantom,2)~=N, error('phantomImage image must be square'); end
        x = double(phantom(:));
    end
else
    x = phantomgallery(phantom,N);
    x = x(:);
end
if isnan(p), p = round(sqrt(2)*N); end

switch CTtype
    case 'parallel'
        if isnan(d), d = p-1; end
        A = paralleltomo(N,angles,p,d,0,sm);
    case 'fancurved'
        if isnan(span), span = 2*atand(1/(2*R-1)); end
        A = fancurvedtomo(N,angles,p,R,span,0,sm);
   otherwise
        error('Type of CT problem not provided')
end
if sm
    b = A*x;
else
    b = A(x,'notransp');
end

% Providing information about the test problem.
ProbInfo.problemType = 'tomography';
ProbInfo.xType = 'image2D';
ProbInfo.bType = 'image2D';
ProbInfo.xSize = [N,N];
ProbInfo.bSize = [p,length(angles)];

repo: IRtools


text:
function calpak_test24 ( )

%*****************************************************************************80
%
%% CALPAK_TEST24 tests JED_TO_YMDF_PERSIAN and YMDF_TO_JED_PERSIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST24\n' );
  fprintf ( 1, '  For the Persian calendar:\n' );
  fprintf ( 1, '  JED_TO_YMDF_PERSIAN: JED -> YMDF.\n' );
  fprintf ( 1, '  YMDF_TO_JED_PERSIAN: YMDF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YMDF                JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_persian ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, m2, d2, f2 ] = jed_to_ymdf_persian ( jed1 );

      s2 = ymdf_to_s_numeric ( y2, m2, d2, f2 );

      jed3 = ymdf_to_jed_persian ( y2, m2, d2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
function Y = SurfStatSmooth( Y, surf, FWHM );

%Smooths surface data by repeatedly averaging over edges.
%
% Usage: Y = SurfStatSmooth( Y, sv, FWHM );
%
% Y        = n x v or n x v x k matrix of surface data, v=#vertices;
%            n=#observations; k=#variates, or memory map of same.
% surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
% or
% surf.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume).
% FWHM     = approximate FWHM of Gaussian smoothing filter, in mesh units.
%
% Note that if the data is memory mapped, then the data is overwriten by
% the smoothed data.

niter=ceil(FWHM^2/(2*log(2)));

if isnumeric(Y)
    [n,v,k]=size(Y);
    isnum=true;
else
    Ym=Y;
    s=Ym.Format{2};
    if length(s)==2
        s=s([2 1]);
        k=1;
    else
        s=s([3 1 2]);
        k=s(3);
    end
    n=s(1);
    v=s(2);
    isnum=false;
end

edg=SurfStatEdg(surf);

Y1=accumarray(edg(:,1),2,[v 1])'+accumarray(edg(:,2),2,[v 1])';

if n>1
    fprintf(1,'%s',[num2str(n) ' x ' num2str(k) ' surfaces to smooth, % remaining: 100 ']);
end
n10=floor(n/10);
for i=1:n
    if rem(i,n10)==0
        fprintf(1,'%s',[num2str(100-i/n10*10) ' ']);
    end
    for j=1:k
        if isnum
            Ys=squeeze(Y(i,:,j));
            for iter=1:niter
                Yedg=Ys(edg(:,1))+Ys(edg(:,2));
                Ys=(accumarray(edg(:,1),Yedg',[v 1]) + ...
                    accumarray(edg(:,2),Yedg',[v 1]))'./Y1;
            end
            Y(i,:,j)=Ys;
        else
            if length(s)==2
                Y=Ym.Data(1).Data(:,i);
            else
                Y=Ym.Data(1).Data(:,j,i);
            end
            for iter=1:niter
                Yedg=Y(edg(:,1))+Y(edg(:,2));
                Y=(accumarray(edg(:,1),Yedg',[v 1]) + ...
                    accumarray(edg(:,2),Yedg',[v 1]))'./Y1;
            end
            if length(s)==2
                Ym.Data(1).Data(:,i)=Y;
            else
                Ym.Data(1).Data(:,j,i)=Y;
            end
        end
    end
end
if n>1
    fprintf(1,'%s\n','Done');
end
if ~isnum
    Y=Ym;
end

return
end

author: brainstorm-tools


text:
function [u,s] = cubRsol(ul,ur,delta)
%
% A wrapper for a scalar Riemannsolver
% to fit with genFrontTrack(...). Uses the flux function
% f=scalflux(u);
%
flux='cub';

[u,s]=ScalarRiemannsol(ul,ur,delta,flux);

end

repo: Aero-matlab


text:
function sF = norm(sVF)
% pointwise norm of the vectorfield
%
% Syntax
%   norm(sVF)
%
% Output
%   sF - S2FunTri
%

sF = S2FunTri(sVF.tri,norm(sVF.values));

end

author: mtex-toolbox

repo: mtex

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


text:
function r8po_print ( n, a, title )

%*****************************************************************************80
%
%% R8PO_PRINT prints a R8PO matrix.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of an SPO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the R8PO matrix.
%
%    Input, string TITLE, a title to be printed.
%
  r8po_print_some ( n, a, 1, 1, n, n, title );

  return
end