text:
function ccvt_box ( )

%*****************************************************************************80
%
%% CCVT_BOX generates and animates a constrained CVT dataset.
%
%  Discussion:
%
%    This program is meant to be used interactively.  It's also
%    possible to prepare a simple input file beforehand and use it
%    in batch mode.
%
%    The program requests input values from the user:
%
%    * NDIM, the spatial dimension,
%    * N, the number of points to generate,
%    * NPP, the number of sample points to use on the boundary,
%    * INIT, initialize the points:
%      ** file, by reading data from file;
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * IT_MAX, the maximum number of iterations;
%    * IT_FIXED, the number of iterative steps to take
%      using a fixed set of sampling points.
%    * SAMPLE, how to conduct the sampling:
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * SAMPLE_NUM, the number of sampling points;
%    * MOVIE_NAME, the name of the movie file to create;
%
%    To indicate that no further computations are desired, it is
%    enough to input a nonsensical value, such as -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%
  DEBUG = 1;
  seed_init = 123456789;

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCVT_BOX\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Generate and animate a constrained CVT dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program is meant to be used interactively.\n' );
  fprintf ( 1, '  It is also possible to prepare a simple input \n' );
  fprintf ( 1, '  file beforehand and use it in batch mode.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * NDIM, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * NPP, the number of points to place on the boundary.\n' );
  fprintf ( 1, '  * SEED, a seed to use for random number generation,\n' );
  fprintf ( 1, '  * INIT, initialize the points:\n' );
  fprintf ( 1, '    ** file, read data from a file;\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * IT_MAX, the maximum number of iterations.\n' );
  fprintf ( 1, '  * IT_FIXED, the number of iterative steps to take\n' );
  fprintf ( 1, '    using a fixed set of sampling points.\n' );
  fprintf ( 1, '  * SAMPLE, how to conduct the sampling.\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * SAMPLE_NUM, the number of sample points;\n' );
  fprintf ( 1, '  * BATCH, the number of sampling points to create at one time;\n' );
  fprintf ( 1, '  * OUTPUT, a file into which the data is stored.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To indicate that no further computations are \n' );
  fprintf ( 1, '  desired, it is enough to input a nonsensical value, \n' );
  fprintf ( 1, '  such as -1.\n' );

  while ( 1 )

    fprintf ( 1, '  *\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '*  Ready to generate a new dataset:\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '  *\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NDIM is the spatial dimension.\n' );
    fprintf ( 1, '  (Try ''2'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    ndim = [];
    ndim = input ( '  Enter NDIM:  ' );

    fprintf ( 1, '  NDIM = %12d\n', ndim );

    if ( ndim < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of NDIM = %12d\n', ndim );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points to generate.\n' );
    fprintf ( 1, '  (Try ''25'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    n = [];
    n = input ( '  Enter N:  ' );

    fprintf ( 1, '  User input N = %12d\n', n );

    if ( n < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of N = %12d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  NPP is the number of sample points used to check\n' );
    fprintf ( 1, '  the boundary.\n' );
    fprintf ( 1, '  (Try ''1000'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    npp = [];
    npp = input ( '  Enter NPP:  ' );

    fprintf ( 1, '  User input NPP = %12d\n', npp );

    if ( npp < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of NPP = %12d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  SEED is a seed for the random number generation.\n' );
%   fprintf ( 1, '  (Try ''123456789'' if you have no preference.)\n' );
%   fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
%   seed = [];
%   seed = input ( '  Enter SEED:  ' );

    seed = seed_init;
    fprintf ( 1, '  Default SEED = %d\n', seed );

    if ( seed < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of SEED = %12d\n', seed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  INIT is the method of initializing the data:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  file       read data from a file;\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, ' \n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    init_string = [];
    init_string = input ( '  Enter INIT:  ' );

    fprintf ( 1, '  User input INIT = "%s".\n', init_string );

    if ( s_len_trim ( init_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( init_string, 'RAND'  ) )
      init = -1;
    elseif ( s_eqi ( init_string, 'RANDOM' ) )
      init_string = 'RAND';
      init = -1;
    elseif ( s_eqi ( init_string, 'UNIFORM' ) )
      init = 0;
    elseif ( s_eqi ( init_string, 'HALTON'  ) )
      init = 1;
    elseif ( s_eqi ( init_string, 'GRID'    ) )
      init = 2;
    elseif ( 0 < s_len_trim ( init_string ) )
      init = 3;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_MAX is the maximum number of iterations.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  An iteration carries out the following steps:\n' );
    fprintf ( 1, '  * the Voronoi region associated with each\n' );
    fprintf ( 1, '    generator is estimated by sampling;\n' );
    fprintf ( 1, '  * the centroid of each Voronoi region is estimated.\n' );
    fprintf ( 1, '  * the generator is replaced by the centroid.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  If "enough" sampling points are used,\n' );
    fprintf ( 1, '  and "enough" iterations are taken, this process\n' );
    fprintf ( 1, '  will converge.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''50'' if you have no preference.)\n' );
    fprintf ( 1, '  (A negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_max = [];
    it_max = input ( '  Enter IT_MAX:  ' );

    fprintf ( 1, '  User input IT_MAX = %12d\n', it_max );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of IT_MAX = %12d\n', it_max );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_FIXED is the number of consecutive iterations\n' );
    fprintf ( 1, '  to take with a fixed set of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Setting IT_FIXED to 1 means a new set of sample\n' );
    fprintf ( 1, '  points is generated on every iterative step;\n' );
    fprintf ( 1, '  Setting IT_FIXED equal to IT_MAX means a single set\n' );
    fprintf ( 1, '  of sample points is used for the entire iteration.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Any value between 1 and IT_MAX is reasonable.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try "%d" if you do not have a preference).\n', it_max );
    fprintf ( 1, '  (A 0 or negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_fixed = [];
    it_fixed = input ( '  Enter IT_FIXED:  ' );

    fprintf ( 1, '  User input IT_FIXED = %12d\n', it_fixed );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of IT_FIXED = %12d\n', it_fixed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE is the method of sampling the region:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, '\n' );

    sample_string = [];
    sample_string = input ( '  Enter SAMPLE:  ' );

    fprintf ( 1, '  User input SAMPLE = "%s".\n', sample_string );

    if ( s_len_trim ( sample_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( sample_string, 'RAND'  ) )
      sample = -1;
    elseif ( s_eqi ( sample_string, 'RANDOM' ) )
      sample = -1;
      sample_string = 'RAND';
    elseif ( s_eqi ( sample_string, 'UNIFORM' ) )
      sample = 0;
    elseif ( s_eqi ( sample_string, 'HALTON'  ) )
      sample = 1;
    elseif ( s_eqi ( sample_string, 'GRID'    ) )
      sample = 2;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM is the number of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The Voronoi regions will be explored by generating\n' );
    fprintf ( 1, '  SAMPLE_NUM points.  For each sample point, the\n' );
    fprintf ( 1, '  nearest generator is found.  Using more points\n' );
    fprintf ( 1, '  gives a better estimate of these regions.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM should be much larger than N, the\n' );
    fprintf ( 1, '  number of generators.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''10000'' if you have no preference.)\n' );
    fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
    fprintf ( 1, '\n' );

    sample_num = [];
    sample_num = input ( '  Enter SAMPLE_NUM:  ' );

    fprintf ( 1, '  User input SAMPLE_NUM = %12d\n', sample_num );

    if ( sample_num <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of SAMPLE_NUM = %12d\n', sample_num );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  BATCH is the number of sample points to create\n' );
%   fprintf ( 1, '  at one time\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  BATCH should be between 1 and SAMPLE_NUM.\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  It is FASTER to set BATCH to SAMPLE_NUM;\n' );
%   fprintf ( 1, '  setting BATCH to 1 requires the least memory.\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  (Try ''%d'' if you have no preference.)\n', ...
%     min ( 1000, sample_num ) );
%   fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
%   fprintf ( 1, '\n' );

%   batch = [];
%   batch = input ( '  Enter BATCH:  ' );
    batch = 1000;

    fprintf ( 1, '  BATCH = %12d\n', batch );

    if ( batch <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of BATCH = %12d\n', batch );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  MOVIE_NAME is the name of the file in which the movie is stored;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''movie.avi'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    movie_name = [];
    movie_name = input ( '  Enter MOVIE_NAME:  ' );

    fprintf ( 1, '  User input MOVIE_NAME = %s\n', movie_name );
%
%  If a file already exists with the name used by the movie file,
%  we need to delete it.
%
    if ( file_exist ( movie_name ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Deleting old copy of %s\n', movie_name );
      file_delete ( movie_name );
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  OUTPUT is a file into which the data is stored;\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  (Try ''cvt.txt'' if you have no preference.)\n' );
%   fprintf ( 1, '  (A blank value terminates execution).\n' );
%   fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
%   fprintf ( 1, ' \n' );

%   file_out_name = [];
%   file_out_name = input ( '  Enter OUTPUT:  ' );

    file_out_name = 'cvt.txt';

    fprintf ( 1, '  OUTPUT = "%s".\n', file_out_name );

    if ( s_len_trim ( file_out_name ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of OUTPUT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_len_trim ( file_out_name ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_BOX\n' );
      fprintf ( 1, '  The input value of OUTPUT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end
%
%  Initialize the data.
%
    seed_iter = seed_init;
    seed = seed_init;
    initialize = 1;

    if ( init == 3 )
      generator = data_read ( sample_string, ndim, n );
    else
      [ generator, seed ] = cvt_sample ( ndim, n, n, init, initialize, seed );
    end
%
%  If the initialization and sampling steps use the same random number
%  scheme, then the sampling scheme does not have to be initialized.
%
    if ( init == sample )
      initialize = 0;
    else
      initialize = 1;
    end

    seed = seed_iter;

    [ energy, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
      sample_num, seed, generator );

    initialize = 0;
%
%  Project points to the boundary.
%
    generator = mpb ( ndim, n, generator, npp );

    seed = seed_iter;
    it_num = 0;
    it_diff = 0.0;

    [ energy2, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
      sample_num, seed, generator );
%
%  Delete any old copy of the movie file.
%
    if ( file_exist ( movie_name ) )
      file_delete ( movie_name )
    end

    if ( DEBUG )
      fprintf ( 1, '\n' );
      fprintf ( 1, ...
        '                            Energy         Energy after\n' );
      fprintf ( 1, '  Iteration     Seed        of CVT         projection\n' );
      fprintf ( 1, '\n' );

      fprintf ( 1, '  %6d  %12d  %14f  %14f\n', ...
        it_num, seed_iter, energy, energy2 );
    end

    cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
      it_fixed, it_num, it_diff, energy2, sample_string, sample_num, ...
      generator, 'initial.txt' );

    points_eps ( 'initial.eps', ndim, n, generator, 'Initial generators' );

    num_frames_per_second = 10;
    aviobj = avifile ( movie_name, 'fps', num_frames_per_second );
    box = [ 0.0, 0.0; 1.0, 0.0; 1.0, 1.0; 0.0, 1.0 ]';
%
%  Start the iteration.
%
    it_diff = 0.0;
    seed_iter = seed_init;
    seed = seed_init;

    for it_num = 1 : it_max
%
%  Once every IT_FIXED steps, update the base value of SEED.
%
%  Otherwise, reset the value of SEED to SEED_BASE.
%
      if ( mod ( it_num - 1, it_fixed ) == 0 )
        seed_iter = seed;
      else
        seed = seed_iter;
      end

      av_points(1:ndim,1:n) = 0.0;
      count(1:n) = 0;

      for j = 1 : sample_num

        [ s, seed ] = cvt_sample ( ndim, sample_num, 1, sample, ...
          initialize, seed );

        nearest = find_closest ( ndim, n, 1, s, generator );

        av_points(1:ndim,nearest(1)) = av_points(1:ndim,nearest(1)) + s(1:ndim);

        count(nearest(1)) = count(nearest(1)) + 1;

      end
%
%  Replace the generators by the average of the sample points.
%
      for j = 1 : n

        if ( count(j) ~= 0 )
          generator(1:ndim,j) = av_points(1:ndim,j) / count(j);
        end

      end

      seed = seed_iter;

      [ energy, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
        sample_num, seed, generator );
%
%  Plot the raw points.
%
      scatter ( generator(1,1:n), generator(2,1:n), [], 'r', 'filled' )
      voronoi ( generator(1,1:n), generator(2,1:n) );
      line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ] );
      axis ( [ -0.05, 1.05, -0.05, 1.05 ] );
      axis square

      xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16 );

      ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16, 'Rotation', 0 );

      it_string = sprintf ( 'Step %4d', it_num );

      title ( it_string, 'FontName', 'Helvetica', 'FontWeight', ...
        'bold', 'FontSize', 16 );

      frame = getframe ( gca );
      aviobj = addframe ( aviobj, frame );
%
%  Apply Lili's projection method.
%  In this case, the energy is changed, so it must be recalculated.
%
      generator = mpb ( ndim, n, generator, npp );

      seed = seed_iter;

      [ energy2, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
        sample_num, seed, generator );

      if ( DEBUG )
        fprintf ( 1, '  %6d  %12d  %14f  %14f\n', ...
          it_num, seed_iter, energy, energy2 );
      end
%
%  Plot the projected points.
%
      scatter ( generator(1,1:n), generator(2,1:n), [], 'r', 'filled' )
      voronoi ( generator(1,1:n), generator(2,1:n) );
      line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ] );
      axis ( [ -0.05, 1.05, -0.05, 1.05 ] );
      axis square

      xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16 );

      ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16, 'Rotation', 0 );

      it_string = sprintf ( 'Step %4d+', it_num );

      title ( it_string, 'FontName', 'Helvetica', 'FontWeight', ...
        'bold', 'FontSize', 16 );

      frame = getframe ( gca );
      aviobj = addframe ( aviobj, frame );

    end

    aviobj = close ( aviobj );
%
%  Write the final data to a file.
%
    cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
      it_fixed, it_num, it_diff, energy, sample_string, sample_num, ...
      generator, 'final.txt' );

    points_eps ( 'final.eps', ndim, n, generator, 'Final generators' );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  The data was written to the file "%s".\n', ...
      file_out_name );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Press RETURN to proceed.\n' );

    pause

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final value of SEED = %d\n', seed );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function planC = getplanCDownSample(planC, optS, index)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

for i = 1: size(planC{indexS.scan}.scanArray,3)
    planC{indexS.scan}.scanInfo(1,i).grid1Units = index*planC{indexS.scan}.scanInfo(1,i).grid1Units;
    planC{indexS.scan}.scanInfo(1,i).grid2Units = index*planC{indexS.scan}.scanInfo(1,i).grid2Units;
    planC{indexS.scan}.scanInfo(1,i).sizeOfDimension2 = size(planC{indexS.scan}.scanArray,2)/index;
    planC{indexS.scan}.scanInfo(1,i).sizeOfDimension1 = size(planC{indexS.scan}.scanArray,1)/index;
end

planC{indexS.scan}.scanArray = getDownsample3(planC{indexS.scan}.scanArray,index,1);

planC{indexS.scan}.scanArraySuperior = [];
planC{indexS.scan}.scanArrayInferior = [];
planC{indexS.scan}.uniformScanInfo = [];

for i=1:length(planC{indexS.structures})
    planC{indexS.structures}(i).rasterSegments =[];
end

planC = getRasterSegs(planC, optS);

planC = setUniformizedData(planC, optS);

try
	for i=1:length(planC{indexS.dose})
        planC = clearCache(planC, i);
	end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function test_bug1832

% MEM 2gb
% WALLTIME 00:20:00
% DEPENDENCY ft_read_mri ft_volumesegment ft_prepare_headmodel ft_prepare_sourcemodel

% for the warp template example script:
% 1) make test script that checks on the units in these objects (with assert)
% 2) also go through the rest of the example script and check on units

template = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/external/spm8/templates/T1.nii'));
template.coordsys = 'spm'; % so that FieldTrip knows how to interpret the coordinate system

% segment the template brain and construct a volume conduction model (i.e. head model): this is needed
% for the inside/outside detection of voxels.
cfg = [];
template_seg = ft_volumesegment(cfg, template);
assert(isfield(template_seg, 'unit'), 'unit field in seg missing')

cfg = [];
cfg.method = 'singlesphere';
template_vol = ft_prepare_headmodel(cfg, template_seg);
assert(isfield(template_vol, 'unit'), 'unit field in vol missing')
% construct the dipole grid in the template brain coordinates
% the source units are in mm
% the negative inwardshift means an outward shift of the brain surface for inside/outside detection
cfg = [];
cfg.sourcemodel.xgrid  = -200:10:200;
cfg.sourcemodel.ygrid  = -200:10:200;
cfg.sourcemodel.zgrid  = -200:10:200;
cfg.sourcemodel.tight  = 'yes';
cfg.inwardshift = -1.5;
cfg.headmodel  = template_vol;
template_grid  = ft_prepare_sourcemodel(cfg);
assert(isfield(template_grid, 'unit'), 'unit field in grid missing')

% prepare leadfield
cfg = [];
cfg.sourcemodel = ft_convert_units(template_grid, 'cm');
cfg.headmodel = template_vol;
cfg.channel = 'EEG';
cfg.elec = ft_read_sens('standard_1020.elc');
template_grid_lf     = ft_prepare_leadfield(cfg);
assert(isfield(template_grid_lf, 'unit'), 'unit field in leadfield missing')

template_grid_lf = ft_convert_units(template_grid_lf, template_grid.unit);
assert(all(template_grid_lf.pos(:) == template_grid.pos(:)), 'leadfield positions differ from sourcemodel position');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ tvgrad ] = gradientTVnorm(f,type)
%GRADIENTTVNORM Computes the gradient of the TV-norm fucntional
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
if ndims(f)==2
    tvgrad=gradientTVnormBackward2D(f);
    return;
end

if strcmp(type,'central')
    warning('It seems that central does not give correct results. Please consider using back or forward')
    tvgrad= gradientTVnormCentral(f);
    return;
end
if strcmp(type,'backward')
    tvgrad= gradientTVnormBackward(f);
    return;
end
if strcmp(type,'forward')
    tvgrad= gradientTVnormForward(f);
    return;
end
error('Undefined type of gradient. Check spelling');
end
%% Backward differences
function tvg=gradientTVnormBackward(f)
Gx=diff(f,1,1);
Gy=diff(f,1,2);
Gz=diff(f,1,3);
tvg=zeros(size(f));
clear f
% these are not defined, but we will define them just for indexing
% readability. They should never be used.
Gx=cat(1,zeros(size(Gx(1,:,:))),Gx);
Gy=cat(2,zeros(size(Gy(:,1,:))),Gy);
Gz=cat(3,zeros(size(Gz(:,:,1))),Gz);

nrm=safenorm(Gx,Gy,Gz);

tvg(1:end,1:end,1:end)= tvg(1:end,1:end,1:end)+(Gx(1:end,1:end,1:end)+Gy(1:end,1:end,1:end)+Gz(1:end,1:end,1:end))./nrm(1:end,1:end,1:end);
tvg(1:end-1,:,:)=tvg(1:end-1,:,:)-Gx([1:end-1]+1,:,:)./nrm([1:end-1]+1,:,:);
tvg(:,1:end-1,:)=tvg(:,1:end-1,:)-Gy(:,[1:end-1]+1,:)./nrm(:,[1:end-1]+1,:);
tvg(:,:,1:end-1)=tvg(:,:,1:end-1)-Gz(:,:,[1:end-1]+1)./nrm(:,:,[1:end-1]+1);

end
%%
function tvg=gradientTVnormBackward2D(f)
Gx=diff(f,1,1);
Gy=diff(f,1,2);
tvg=zeros(size(f));
clear f
% these are not defined, but we will define them just for indexing
% readability. They should never be used.
Gx=cat(1,zeros(size(Gx(1,:))),Gx);
Gy=cat(2,zeros(size(Gy(:,1))),Gy);

nrm=safenorm(Gx,Gy,0);
nrmx=safenorm(Gx,Gy,0);
nrmy=safenorm(Gx,Gy,0);

tvg(1:end,1:end)= tvg(1:end,1:end)+(Gx(1:end,1:end)+Gy(1:end,1:end))./nrm(1:end,1:end);
tvg(1:end-1,:)=tvg(1:end-1,:)-Gx([1:end-1]+1,:)./nrm([1:end-1]+1,:);
tvg(:,1:end-1)=tvg(:,1:end-1)-Gy(:,[1:end-1]+1)./nrm(:,[1:end-1]+1);
end

%% Forward differences
function tvg=gradientTVnormForward(f)
Gx=diff(f,1,1);
Gy=diff(f,1,2);
Gz=diff(f,1,3);
tvg=zeros(size(f));
clear f
% these are not defined, but we will define them just for indexing
% readability. They should never be used.
Gx=cat(1,Gx,zeros(size(Gx(end,:,:))));
Gy=cat(2,Gy,zeros(size(Gy(:,end,:))));
Gz=cat(3,Gz,zeros(size(Gz(:,:,end))));
nrm=safenorm(Gx,Gy,Gz);
tvg(1:end-1,1:end-1,1:end-1)=tvg(1:end-1,1:end-1,1:end-1)-(Gx(1:end-1,1:end-1,1:end-1)+Gy(1:end-1,1:end-1,1:end-1)+Gz(1:end-1,1:end-1,1:end-1))./nrm(1:end-1,1:end-1,1:end-1);
tvg(2:end,:,:)=tvg(2:end,:,:)+Gx([2:end]-1,:,:)./nrm([2:end]-1,:,:);
tvg(:,2:end,:)=tvg(:,2:end,:)+Gy(:,[2:end]-1,:)./nrm(:,[2:end]-1,:);
tvg(:,:,2:end)=tvg(:,:,2:end)+Gz(:,:,[2:end]-1)./nrm(:,:,[2:end]-1);
end
%% Central differences
% Doesnt accound for edges of the image!!!
% https://math.stackexchange.com/questions/1612017/gradient-of-the-tv-norm-of-an-image
function tvg=gradientTVnormCentral(f)
[Gx,Gy,Gz]=gradient(f);
tvg=zeros(size(Gx));
nrm=safenorm(Gx,Gy,Gz);
tvg(2:end,:,:)   = tvg(2:end,:,:)   + Gx(1:end-1,:,:)./nrm(1:end-1,:,:);
tvg(1:end-1,:,:) = tvg(1:end-1,:,:) - Gx(2:end,:,:)  ./nrm(2:end,:,:);
tvg(:,2:end,:)   = tvg(:,2:end,:)   + Gy(:,1:end-1,:)./nrm(:,1:end-1,:);
tvg(:,1:end-1,:) = tvg(:,1:end-1,:) - Gy(:,2:end,:)  ./nrm(:,2:end,:);
tvg(:,:,2:end)   = tvg(:,:,2:end)   + Gz(:,:,1:end-1)./nrm(:,:,1:end-1);
tvg(:,:,1:end-1) = tvg(:,:,1:end-1) - Gz(:,:,2:end)  ./nrm(:,:,2:end);
tvg = tvg/2;
end

%% Utils
function nrm=safenorm(Gx,Gy,Gz)

nrm=sqrt(Gx.^2+Gy.^2+Gz.^2)+0.00000001;
nrm(nrm==0)=1;

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function [ F ] = gsp_ngwft_frame_matrix(G, g, param)
%GSP_NGWFT_FRAME_MATRIX Create the matrix of the GWFT frame
%   Usage:  F = gsp_ngwft_frame_matrix(G, g, param);
%           F = gsp_ngwft_frame_matrix(G, g);
%
%   Input parameters:
%         G     : Graph
%         g     : window
%         param : Structure of optional parameter
%   Output parameters:
%         F     : Frame
%
%   This function compute the graph windowed Fourier transform frame. In
%   might require a lot of memory. we create a matrix of size N^2*N.
%
%   *param* a Matlab structure containing the following fields:
%
%     * *param.verbose* : 0 no log, 1 print main steps, 2 print all steps.
%     By default, it is 1.

%   AUTHOR : Nathanael Perraudin
%   TESTING: test_fast_gwft

% Optional input arguments
if nargin<4, param=struct; end
if ~isfield(param, 'verbose'), param.verbose=1 ; end

if param.verbose && G.N>256
    warning('Create a big matrix, you can use other methods.');
end

% Compute the translate of g
ghat=G.U'*g;
Ftrans=sqrt(G.N)*G.U*(repmat(ghat,1,G.N).*G.U');

F = gsp_repmatline(Ftrans,1,G.N) .* ...
    repmat(repmat(1./G.U(:,1),1,G.N).*G.U,1,G.N);

% Normalization
F = F./repmat(sqrt(sum(abs(F).^2,1)),G.N,1);

end

repo: gspbox


text:
function a = log10(a)
%LOG10        Taylor logarithm  log10(a)
%

% written  05/22/09     S.M. Rump
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  a = log(a);
  if isa(a.t,'intval')
    INTLAB_STDFCTS_LOG10_ = getappdata(0,'INTLAB_STDFCTS_LOG10_');
    a.t = a.t .* infsup(INTLAB_STDFCTS_LOG10_.INF,INTLAB_STDFCTS_LOG10_.SUP);
  else
    a.t = a.t / log(10);
  end

repo: OpenMAS


text:
function [ output_meta ] = meta2sicd_iceye(hdffile)
%META2SICD_ICEYE Converts ICEYE HDF5 into a SICD metadata structure
%
% Written by: Tim Cox, NRL; Wade Schwartzkopf, NGA/Reseach
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

SECONDS_IN_A_DAY = 24*60*60;

%% Collection Info
output_meta.CollectionInfo.CollectorName = get_hdf_data(hdffile,'/','satellite_name');
output_meta.CollectionInfo.CoreName = get_hdf_data(hdffile,'/','product_name');
output_meta.CollectionInfo.CollectType='MONOSTATIC';
output_meta.CollectionInfo.RadarMode.ModeType = upper(get_hdf_data(hdffile,'/','acquisition_mode'));
output_meta.CollectionInfo.RadarMode.ModeID = get_hdf_data(hdffile,'/','product_type');
output_meta.CollectionInfo.Classification = 'UNCLASSIFIED';

%% Image Creation
output_meta.ImageCreation.Application = ['ICEYE_P_' ...
    num2str(get_hdf_data(hdffile,'/','processor_version'))];
output_meta.ImageCreation.DateTime = datenum(reshape(get_hdf_data(...
    hdffile,'/','processing_time'),1,[]),'yyyy-mm-ddTHH:MM:SS');
output_meta.ImageCreation.Profile='Prototype';

%% Image Data
switch get_hdf_data(hdffile,'/','sample_precision')
    case 'float32'
        output_meta.ImageData.PixelType='RE32F_IM32F';
    case 'int16'
        output_meta.ImageData.PixelType='RE16I_IM16I';
    otherwise
        warning('META2SICD_ICEYE:UnrecognizedSamplePrecision', 'Only float32 and int16 precisions supported.');
end
output_meta.ImageData.NumRows = uint32(get_hdf_data(hdffile,'/','number_of_range_samples'));
output_meta.ImageData.NumCols = uint32(get_hdf_data(hdffile,'/','number_of_azimuth_samples'));
output_meta.ImageData.FirstRow = 0;
output_meta.ImageData.FirstCol = 0;
output_meta.ImageData.FullImage.NumRows = output_meta.ImageData.NumRows;
output_meta.ImageData.FullImage.NumCols = output_meta.ImageData.NumCols;
coord_center = get_hdf_data(hdffile,'/','coord_center');
output_meta.ImageData.SCPPixel.Row = uint32(coord_center(1))-1;
if coord_center(2)>0 && coord_center(2)<output_meta.ImageData.NumCols
    output_meta.ImageData.SCPPixel.Col = uint32(coord_center(2))-1;
    if strcmpi(get_hdf_data(hdffile,'/','look_side'),'left')
        output_meta.ImageData.SCPPixel.Col = output_meta.ImageData.NumCols - ...
            output_meta.ImageData.SCPPixel.Col - 1;
    end
else
    % Bug in earlier version of ICEYE processor sometimes resulted in
    % center pixel outside of image in the column direction.  In this case,
    % we just choose our own.  The choice of SCP is essentially arbitrary
    % anyway.
    output_meta.ImageData.SCPPixel.Col = output_meta.ImageData.NumCols/2;
end

%% GeoData
avg_scene_height = double(get_hdf_data(hdffile,'/','avg_scene_height'));
SCP = geodetic_to_ecf(coord_center(3),coord_center(4),avg_scene_height);
if strcmpi(get_hdf_data(hdffile,'/','geo_ref_system'),'WGS84')
    output_meta.GeoData.EarthModel = 'WGS_84';  % Should always be the case
else
    warning('META2SICD_ICEYE:UnexpectedEarthModel', 'Only WGS84 expected for earth model.');
end
output_meta.GeoData.SCP.ECF.X = SCP(1);
output_meta.GeoData.SCP.ECF.Y = SCP(2);
output_meta.GeoData.SCP.ECF.Z = SCP(3);
output_meta.GeoData.SCP.LLH.Lat = coord_center(3);
output_meta.GeoData.SCP.LLH.Lon = coord_center(4);
output_meta.GeoData.SCP.LLH.HAE = avg_scene_height;
% Image corner coordinates
% We could do it this way, but we would have to swap first/last for left
% looking cases.  Instead, we will let derived_sicd_fields compute these.
% FRFC = get_hdf_data(hdffile,'/','coord_first_near');
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lat=FRFC(3);
% output_meta.GeoData.ImageCorners.ICP.FRFC.Lon=FRFC(4);
% FRLC = get_hdf_data(hdffile,'/','coord_last_near');
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lat=FRLC(3);
% output_meta.GeoData.ImageCorners.ICP.FRLC.Lon=FRLC(4);
% LRFC = get_hdf_data(hdffile,'/','coord_first_far');
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lat=LRFC(3);
% output_meta.GeoData.ImageCorners.ICP.LRFC.Lon=LRFC(4);
% LRLC = get_hdf_data(hdffile,'/','coord_last_far');
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lat=LRLC(3);
% output_meta.GeoData.ImageCorners.ICP.LRLC.Lon=LRLC(4);

%% Grid
% We will derive Grid.Col.SS/ImpRespBW later
% These two computations of Grid.Row.SS should be identical
% output_meta.Grid.Row.SS = get_hdf_data(hdffile,'/','slant_range_spacing');
output_meta.Grid.Row.SS = SPEED_OF_LIGHT/(2*get_hdf_data(hdffile,'/','range_sampling_rate'));
fc = get_hdf_data(hdffile,'/','carrier_frequency'); %Hz
output_meta.Grid.Row.KCtr=2*fc/SPEED_OF_LIGHT;
output_meta.Grid.Col.KCtr=0;
output_meta.Grid.Type='RGZERO'; %based on ICEYE spec, all SLC files are range, zero Doppler
output_meta.Grid.ImagePlane = 'SLANT'; %SLC is always slant
output_meta.Grid.Row.Sgn=-1;
output_meta.Grid.Col.Sgn=-1;
% ICEYE documentation says chirp_bandwidth should be float, but it is int
% in datasets we have seen.
output_meta.Grid.Row.ImpRespBW = 2*double(get_hdf_data(hdffile,'/','chirp_bandwidth'))/SPEED_OF_LIGHT;
output_meta.Grid.Row.DeltaKCOAPoly=0;
output_meta.Grid.Row.WgtType.WindowName = get_hdf_data(hdffile,'/','window_function_range');
output_meta.Grid.Col.WgtType.WindowName = get_hdf_data(hdffile,'/','window_function_azimuth');
% For now only have seen uniform weighted data
if strcmpi(output_meta.Grid.Row.WgtType.WindowName,'NONE')
    output_meta.Grid.Row.WgtType.WindowName = 'UNIFORM';
end
% Although no weighting is applied, ICEYE generally does not remove
% slow-time beam shape, so spectrum will not necessarily be flat.  (This is
% reflected in the STBeamComp parameter in SICD.)
if strcmpi(output_meta.Grid.Col.WgtType.WindowName,'NONE')
    output_meta.Grid.Col.WgtType.WindowName = 'UNIFORM';
end
% ImpRespWid will be computed in derived_sicd_fields later

%% Timeline
[start_t, start_frac] = datenum_w_frac(reshape(get_hdf_data(hdffile,'/','acquisition_start_utc'),1,[]));
[end_t, end_frac] = datenum_w_frac(reshape(get_hdf_data(hdffile,'/','acquisition_end_utc'),1,[]));
output_meta.Timeline.CollectStart = start_t + (start_frac/SECONDS_IN_A_DAY);
output_meta.Timeline.CollectDuration = round((end_t-start_t)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
        (end_frac-start_frac); % Handle fractional seconds
output_meta.Timeline.IPP.Set.TStart = 0;
output_meta.Timeline.IPP.Set.TEnd = output_meta.Timeline.CollectDuration;
acquisition_prf = get_hdf_data(hdffile,'/','acquisition_prf');
output_meta.Timeline.IPP.Set.IPPStart = 0;
% We don't realy know how many IPPs.  Just making data consistent
output_meta.Timeline.IPP.Set.IPPEnd = round(acquisition_prf*...
    (output_meta.Timeline.IPP.Set.TEnd-output_meta.Timeline.IPP.Set.TStart));
output_meta.Timeline.IPP.Set.IPPPoly = [0; acquisition_prf]; %assume constant PRF

%% Position
state_vector_time_str = get_hdf_data(hdffile,'/','state_vector_time_utc').';
state_vector_T  = zeros(1,size(state_vector_time_str,1));
state_vector_T_frac  = zeros(1,size(state_vector_time_str,1));
for i=1:size(state_vector_time_str,1)
    [state_vector_T(i), state_vector_T_frac(i)] = datenum_w_frac(state_vector_time_str(i,:));
end
% Need to make this polynomial relative to the start time of the collect
t = round((state_vector_T-start_t)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
        (state_vector_T_frac-start_frac); % Handle fractional seconds
output_meta.Position.ARPPoly.X = fliplr(polyfit(t(:),get_hdf_data(hdffile,'/','posX'),5))';
output_meta.Position.ARPPoly.Y = fliplr(polyfit(t(:),get_hdf_data(hdffile,'/','posY'),5))';
output_meta.Position.ARPPoly.Z = fliplr(polyfit(t(:),get_hdf_data(hdffile,'/','posZ'),5))';
% Velocity is in native metadata here, but we will derive it
% velX = get_hdf_data(hdffile,'/','velX');
% velY = get_hdf_data(hdffile,'/','velY');
% velZ = get_hdf_data(hdffile,'/','velZ');

%% RadarCollection
Pol = get_hdf_data(hdffile,'/','polarization');
output_meta.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization = [Pol(1) ':' Pol(2)];
output_meta.RadarCollection.TxPolarization = Pol(1);
TxBandwidth = double(get_hdf_data(hdffile,'/','chirp_bandwidth'));  % Hz
output_meta.RadarCollection.TxFrequency.Min = fc-TxBandwidth/2;
output_meta.RadarCollection.TxFrequency.Max = fc+TxBandwidth/2;
output_meta.RadarCollection.Waveform.WFParameters.TxFreqStart = ...
    output_meta.RadarCollection.TxFrequency.Min;
output_meta.RadarCollection.Waveform.WFParameters.TxRFBandwidth = TxBandwidth;
output_meta.RadarCollection.Waveform.WFParameters.TxPulseLength = ...
    get_hdf_data(hdffile,'/','chirp_duration');  % sec
output_meta.RadarCollection.Waveform.WFParameters.ADCSampleRate = ...
    get_hdf_data(hdffile,'/','range_sampling_rate');  % Hz
output_meta.RadarCollection.Waveform.WFParameters.RcvDemodType = 'CHIRP';
output_meta.RadarCollection.Waveform.WFParameters.RcvFMRate = 0;

%% ImageFormation
output_meta.ImageFormation.RcvChanProc=struct('NumChanProc',1,'PRFScaleFactor',1,'ChanIndex',1);
output_meta.ImageFormation.TxRcvPolarizationProc=[Pol(1) ':' Pol(2)];
output_meta.ImageFormation.ImageFormAlgo='RMA';
output_meta.ImageFormation.TStartProc=0;
output_meta.ImageFormation.TEndProc = output_meta.Timeline.CollectDuration;
output_meta.ImageFormation.TxFrequencyProc.MinProc = ...
    output_meta.RadarCollection.TxFrequency.Min;
output_meta.ImageFormation.TxFrequencyProc.MaxProc = ...
    output_meta.RadarCollection.TxFrequency.Max;
output_meta.ImageFormation.STBeamComp = 'NO';
output_meta.ImageFormation.ImageBeamComp = 'SV';
output_meta.ImageFormation.AzAutofocus='NO';
output_meta.ImageFormation.RgAutofocus='NO';

%% Radiometric
output_meta.Radiometric.BetaZeroSFPoly = get_hdf_data(hdffile,'/','calibration_factor');

%% RMA
output_meta.RMA.RMAlgoType='OMEGA_K';
output_meta.RMA.ImageType='INCA';
% These two computations of near_range are the same:
% near_range = get_hdf_data(hdffile,'/','slant_range_to_first_pixel');
near_range = get_hdf_data(hdffile,'/','first_pixel_time')*SPEED_OF_LIGHT/2;
output_meta.RMA.INCA.R_CA_SCP = near_range + ...
    (double(output_meta.ImageData.SCPPixel.Row)*output_meta.Grid.Row.SS);
output_meta.RMA.INCA.FreqZero=fc;
ss_zd_s = get_hdf_data(hdffile,'/','azimuth_time_interval');
[zd_first, zd_first_frac]=datenum_w_frac(reshape(get_hdf_data(hdffile,'/','zerodoppler_start_utc'),1,[]));
[zd_last, zd_last_frac]=datenum_w_frac(reshape(get_hdf_data(hdffile,'/','zerodoppler_end_utc'),1,[]));
if strcmpi(get_hdf_data(hdffile,'/','look_side'),'left')
    ss_zd_s = -ss_zd_s;
    zd_left = zd_last;
    zd_left_frac = zd_last_frac;
else
    zd_left = zd_first;
    zd_left_frac = zd_first_frac;
end
% Zero doppler time of SCP relative to collect start
zd_t_scp = round((zd_left-start_t)*SECONDS_IN_A_DAY) + ... % Convert days to seconds
    (zd_left_frac-start_frac) + ... % Handle fractional seconds
    (double(output_meta.ImageData.SCPPixel.Col) * ss_zd_s);
% Reference time for Doppler polynomials according to ICEYE documentation
tref = get_hdf_data(hdffile,'/','first_pixel_time') + ...
    double(get_hdf_data(hdffile,'/','number_of_range_samples'))/...
    (2*get_hdf_data(hdffile,'/','range_sampling_rate'));
% Compute DRateSFPoly
pos_coefs = [output_meta.Position.ARPPoly.X(end:-1:1) ...
    output_meta.Position.ARPPoly.Y(end:-1:1) ...
    output_meta.Position.ARPPoly.Z(end:-1:1)];
% Velocity is derivate of position.
vel_coefs=pos_coefs(1:end-1,:).*repmat(((size(pos_coefs,1)-1):-1:1)',[1 3]);
vel_x = polyval(vel_coefs(:,1), zd_t_scp);
vel_y = polyval(vel_coefs(:,2), zd_t_scp);
vel_z = polyval(vel_coefs(:,3), zd_t_scp);
vm_ca_sq = vel_x.^2 + vel_y.^2 + vel_z.^2; % Magnitude of the velocity squared
r_ca = [output_meta.RMA.INCA.R_CA_SCP; 1]; % Polynomial representing range as a function of range distance from SCP
% Shift 1D polynomial to account for SCP (since we picked SCP to match
% reference point, this shouldn't do anything)
dop_rate_coeffs = get_hdf_data(hdffile,'/','doppler_rate_coeffs');
% Prior to ICEYE_P_1.14 processor, absolute value of Doppler rate was
% provided, not true Doppled rate.  Since Doppler rate is always negative
% will will flip sign appropriately.
dop_rate_coeffs = -sign(dop_rate_coeffs(1))* dop_rate_coeffs;
dop_rate_poly_rg_shifted=polyshift(dop_rate_coeffs, ...
    (2*output_meta.RMA.INCA.R_CA_SCP/SPEED_OF_LIGHT-tref));  % This offset should be zero because of how we have chosen SCP
% Scale 1D polynomial to from Hz/s^n to Hz/m^n
dop_rate_poly_rg_scaled=dop_rate_poly_rg_shifted.*...
    (2/SPEED_OF_LIGHT).^(0:(length(dop_rate_poly_rg_shifted)-1)).';
output_meta.RMA.INCA.DRateSFPoly = - conv(dop_rate_poly_rg_scaled,r_ca) * ... % Multiplication of two polynomials is just a convolution of their coefficients
    SPEED_OF_LIGHT / (2 * fc * vm_ca_sq(1)); % Assumes a SGN of -1
% Fields dependent on Doppler rate
% This computation of SS is actually better than the claimed SS
% (azimuth_ground_spacing) in many ways, because this makes all of the metadata
% internally consistent.  This must be the sample spacing exactly at SCP
% (which is the definition for SS in SICD), if the other metadata from
% which is it computed is correct and consistent. Since column SS can vary
% slightly over a RGZERO image, we don't know if the claimed sample spacing
% in the native metadata is at our chosen SCP, or another point, or an
% average across image or something else.
% output_meta.Grid.Col.SS = get_hdf_data(hdffile,'/','azimuth_ground_spacing');
output_meta.Grid.Col.SS = sqrt(vm_ca_sq(1)) * abs(ss_zd_s) * ...
    output_meta.RMA.INCA.DRateSFPoly(1,1);
% Convert to azimuth spatial bandwidth (cycles per meter)
dop_bw = get_hdf_data(hdffile,'/','total_processed_bandwidth_azimuth');
output_meta.Grid.Col.ImpRespBW = dop_bw*abs(ss_zd_s)/output_meta.Grid.Col.SS;
output_meta.RMA.INCA.TimeCAPoly = [zd_t_scp; ss_zd_s/output_meta.Grid.Col.SS];
% Doppler Centroid
dc_estimate_coeffs = get_hdf_data(hdffile,'/','dc_estimate_coeffs');
dc_estimate_time_utc = get_hdf_data(hdffile,'/','dc_estimate_time_utc').';
[dc_t, dc_t_frac, dc_zd_t]  = deal(zeros(1,size(dc_estimate_time_utc,1)));
for i=1:size(dc_estimate_time_utc,1)
    [dc_t(i), dc_t_frac(i)] = datenum_w_frac(dc_estimate_time_utc(i,:));
    dc_zd_t(i) = round((dc_t(i)-start_t)*SECONDS_IN_A_DAY) + ... % Convert days to seconds
        (dc_t_frac(i)-start_frac);  % Handle fractional seconds
end
M=49; % Arbitrary.  Seems as good as any.
f_dc=zeros(M,numel(dc_zd_t));
diff_t_rg = get_hdf_data(hdffile,'/','first_pixel_time') + linspace(0, ...
    double(get_hdf_data(hdffile,'/','number_of_range_samples'))/...
    get_hdf_data(hdffile,'/','range_sampling_rate'), M) - tref;
for n=1:numel(dc_zd_t)
    % Sampled values of Doppler centroid
    f_dc(:,n)=polyval(dc_estimate_coeffs(end:-1:1,n),diff_t_rg);
end
t_coa = repmat(dc_zd_t(:)', M, 1)+(f_dc./dop_rate_poly_rg_scaled(1,1));

% For SICD, we need polynomial as a function of range/azimuth distance from SCP
% in meters.
t_rg_scp = get_hdf_data(hdffile,'/','first_pixel_time') + ...
    double(output_meta.ImageData.SCPPixel.Row)/...
    get_hdf_data(hdffile,'/','range_sampling_rate');
range_scp_m = repmat((diff_t_rg(:) + tref - t_rg_scp) * (SPEED_OF_LIGHT/2), 1, numel(dc_zd_t));
azimuth_scp_m = repmat(output_meta.Grid.Col.SS * (dc_zd_t(:)' - zd_t_scp) / ss_zd_s, M, 1);

% Least squares fit (section 2.3 in paper)
% Compute the following continuous function for doppler centroid
% f_dc(azimuth_scp_m,range_scp_m)=...
%    x(1,1) + x(2,1)*range_scp_m + x(3,1)*range_scp_m*range_scp_m+...
%    x(1,2)*azimuth_scp_m + x(1,3)*azimuth_scp_m*azimuth_scp_m+...
%    x(2,2)*azimuth_scp_m*range_scp_m;
% A*x = b
% We add x(2,3),x(3,2), and x(3,3), even though they were not in the paper.
a=[ones(numel(f_dc),1) ...
   range_scp_m(:) ...
   range_scp_m(:).^2 ...
   azimuth_scp_m(:) ...
   range_scp_m(:).*azimuth_scp_m(:) ...
   (range_scp_m(:).^2).*azimuth_scp_m(:) ...
   azimuth_scp_m(:).^2 ...
   range_scp_m(:).*(azimuth_scp_m(:).^2) ...
   (range_scp_m(:).^2).*(azimuth_scp_m(:).^2)];
[b_dc,b_coa]=deal(zeros(size(a,2),1));
for i=1:size(a,2)
   b_dc(i)=sum(f_dc(:).*a(:,i)); % Doppler centroid
   b_coa(i)=sum(t_coa(:).*a(:,i)); % Center of aperture time
end
A=zeros(size(a,2));
for i=1:size(a,2)
    for j=1:size(a,2)
        A(i,j)=sum(a(:,i).*a(:,j));
    end
end
old_warning_state=warning('off','MATLAB:nearlySingularMatrix');
x=A\b_dc; % MATLAB often flags this as badly scaled, but results still appear valid
x2=A\b_coa;
warning(old_warning_state);
output_meta.RMA.INCA.DopCentroidPoly=reshape(x,3,3);
output_meta.Grid.Col.DeltaKCOAPoly=...
    output_meta.RMA.INCA.DopCentroidPoly*ss_zd_s/output_meta.Grid.Col.SS;
output_meta.Grid.TimeCOAPoly=reshape(x2,3,3);
% Spotlight adjustments
if strcmp(output_meta.CollectionInfo.RadarMode.ModeType,'SPOTLIGHT')
    output_meta.Grid.TimeCOAPoly = output_meta.Grid.TimeCOAPoly(1);
    % This field required to compute Grid.Col.DeltaKCOAPoly, but not
    % allowed for spotlight in SICD.
    output_meta.RMA.INCA = rmfield(output_meta.RMA.INCA, 'DopCentroidPoly');
else % This field also not allowed for spotlight in SICD.
    output_meta.RMA.INCA.DopCentroidCOA=true;
end

%% GeoData
% Now that sensor model fields have been populated, we can populate
% GeoData.SCP more precisely.
ecf = point_image_to_ground([output_meta.ImageData.SCPPixel.Row;output_meta.ImageData.SCPPixel.Col],output_meta);
output_meta.GeoData.SCP.ECF.X=ecf(1);
output_meta.GeoData.SCP.ECF.Y=ecf(2);
output_meta.GeoData.SCP.ECF.Z=ecf(3);
llh=ecf_to_geodetic([output_meta.GeoData.SCP.ECF.X output_meta.GeoData.SCP.ECF.Y output_meta.GeoData.SCP.ECF.Z]);
output_meta.GeoData.SCP.LLH.Lat=llh(1);
output_meta.GeoData.SCP.LLH.Lon=llh(2);
output_meta.GeoData.SCP.LLH.HAE=llh(3);

%% SCPCOA
output_meta = derived_sicd_fields(output_meta);

end

% Reads HDF5 data
function value = get_hdf_data(hid_t,path,data_name)
gid = H5G.open(hid_t,path);
data = H5D.open(gid,data_name);
value = H5D.read(data);
% Versions of MATLAB prior to 2020a resulted in single-element cell arrays
% for some ICEYE HDF5 fields.
if iscell(value)
    value = value{1};
end
H5G.close(gid);
H5D.close(data);
end

% MATLAB's datenum function won't handle precise times down under a
% millisecond, because 1) It won't accept a format with more than 3 .FFF in
% the string description of the date format and 2) the resulting serial
% date number is stored in days from 00-JAN-0000 and just doesn't have
% enough bits to handle fractional seconds to the level we want.  Here we
% handle the fractional seconds separately so we can read date with the
% precision we need.
function [datenum_s, datenum_frac] = datenum_w_frac(datestring)
    datenum_s = datenum(datestring,'yyyy-mm-ddTHH:MM:SS');
    datenum_frac = str2double(regexp(datestring,'\.\d*','match'));
    if isnan(datenum_frac), datenum_frac = 0; end
end

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

author: ngageoint


text:
function spm_dem_search_trajectory(DEM)
% plots visual search in extrinsic and intrinsic coordinates
% FORMAT spm_dem_search_trajectory(DEM)
%
% DEM - {DEM} structures from visual search simulations
%
% hidden causes and states
%==========================================================================
% x    - hidden states:
%   o(1) - oculomotor angle
%   o(2) - oculomotor angle
%   x(1) - relative amplitude of visual hypothesis 1
%   x(2) - relative amplitude of visual hypothesis 2
%   x(3) - ...
%
% v    - hidden causes
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception - x)
%   g(2) - oculomotor angle (proprioception - y)
%   g(3) - retinal input - channel 1
%   g(4) - retinal input - channel 2
%   g(5) - ...
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dem_search_trajectory.m 4595 2011-12-19 13:06:22Z karl $

% Preliminaries
%--------------------------------------------------------------------------
clf, global STIM
N  = length(DEM);
S  = spm_read_vols(STIM.U);

% Stimulus
%======================================================================
Dx = STIM.U.dim(1)/2;
Dy = STIM.U.dim(2)/2;
a  = [];
q  = [];
c  = [];

subplot(2,2,1); hold off
image((S + 1)*32), axis image, hold on

for i = 1:N

    % i-th saccade - position
    %----------------------------------------------------------------------
    pU = DEM{i}.pU.x{1}(1:2,:)*16;
    T  = length(pU);

    % eye movements in extrinsic coordinates
    %======================================================================
    subplot(2,2,1)
    plot(pU(2,:) + Dy,pU(1,:) + Dx,'r-','LineWidth',2)
    plot(pU(2,T) + Dy,pU(1,T) + Dx,'ro','MarkerSize',16)

    % Free energy
    %======================================================================
    F(i) = DEM{i}.F;

end

% Free energy
%--------------------------------------------------------------------------
subplot(2,2,2)
plot(F)
title('Free-energy','FontSize',16)
xlabel('saccade')
axis square

author: spm


text:
function test_suite = test_steepest_descent
    clear all;
    initTestSuite;
end

function test_1
    A = [3 2; 2 6];
    X = [ 2 1 0 4;
         -2 2 0 4];
    B = A * X;
    solver = SPX_SteepestDescent(A, B);
    x = solver.solve();
    %solver.printResults();
    verifyTrue(testCase, solver.hasConverged());
end

function test_2
    A  = [
    2 -1 0;
    -1 2 -1;
    0 -1 2];
    X = randi(10, 3,3);
    B = A * X;
    solver = SPX_SteepestDescent(A, B);
    % This one doesn't converge in 30 iterations
    solver.MaxIterations = 50;
    x = solver.solve();
    verifyTrue(testCase, solver.hasConverged());
end

author: indigits

repo: sparse-plex


text:
function results = VBA_test (subfolder)
    import matlab.unittest.TestSuite;

    vba_info = VBA_version();
    target = [vba_info.path filesep 'tests'];
    if nargin > 0
        target = [target filesep subfolder];
    end
    suite = TestSuite.fromFolder(target, 'IncludingSubfolders', true);
    results = suite.run ();
end

author: MBB-team

repo: VBA-toolbox


text:
classdef (Abstract) HypothesiserX < BaseX
% HypothesiserX Abstract class
%
% Summary of HypothesiserX:
% This is the base class for all TrackingX hypothesisers.
% Any custom defined hypothesiser should be derived from this HypothesiserX base class.
%
% HypothesiserX Properties:
%   None
%
% HypothesiserX Methods:
%   + HypothesiserX - Constructor method
%   ~ hypothesise - Performs the generation an evaluation of measurement
%                   hypotheses
%
% (+) denotes puplic properties/methods
% (~) denotes abstract properties/methods
%
% February 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
    end

    methods (Abstract)
        hypothesise(this);
    end
    methods
        function this = HypothesiserX(varargin)
        % HYPOTHESISERX Constructor method
        %
        % DESCRIPTION:
        % * HypothesiserX() returns a "HypothesiserX" object handle

            this@BaseX();
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [ y,opts] = lrn( x,N,kappa,alpha,beta,dzdy,opts)
%LRN Summary of this function goes here
%   Detailed explanation goes here

    if isfield(opts,'use_nntoolbox')&&opts.use_nntoolbox==1 %
        if ~isfield(opts,'layer')||length(opts.layer)<opts.current_layer||~isfield(opts.layer{opts.current_layer},'lrn_nntb')
            opts.layer{opts.current_layer}.lrn_nntb = nnet.internal.cnn.layer.LocalMapNorm2D( 'lrn_nntb', N,alpha,beta,kappa);
            if opts.use_gpu
                opts.layer{opts.current_layer}.lrn_nntb = setupForGPUPrediction(opts.layer{opts.current_layer}.lrn_nntb);
            else
                opts.layer{opts.current_layer}.lrn_nntb =setupForHostPrediction(opts.layer{opts.current_layer}.lrn_nntb);
            end
        end
        lrn_nntb=opts.layer{opts.current_layer}.lrn_nntb ;

        if isempty(dzdy)
            y=lrn_nntb.forward(x);
        else
            y=lrn_nntb.forward(x);
            y=lrn_nntb.backward(x,y,dzdy,[]);
        end

       return;
    end

    sz=size(x);
    if (length(sz)>2)
        channel_dim=3; Index={':',':'};%cnn
    else
        channel_dim=1; Index={};%mlp
    end

    y=zeros(size(x),'like',x);
    L=y;
    x2_cumu=cumsum(x.^2,channel_dim);

    n_channels=size(x,channel_dim);
    for n=1:n_channels
       nStart=max(n-floor((N-1)/2),1)-1;
       nEnd=min(n+ceil((N-1)/2),n_channels);
       Index_new=[Index,n,':'];
       Index_End=[Index,nEnd,':'];
       Index_Start=[Index,nStart,':'];
       cumu_start=0;
       if nStart~=0
           cumu_start=x2_cumu(Index_Start{:});
       end
       L(Index_new{:})=kappa+alpha*(x2_cumu(Index_End{:}) -cumu_start);

    end

    if isempty(dzdy)
       y=x./L.^(beta);
    else

        temp_cumu=cumsum(dzdy.*x./L.^(1+beta),channel_dim);
        LL=zeros(size(L),'like',x);
        for n=1:size(x,channel_dim)
            nStart=max(n-floor((N-1)/2),1)-1;
            nEnd=min(n+ceil((N-1)/2),n_channels);
            Index_new=[Index,n,':'];
            Index_End=[Index,nEnd,':'];
            Index_Start=[Index,nStart,':'];
            cumu_start=0;
           if nStart~=0
               cumu_start=temp_cumu(Index_Start{:});
           end
           LL(Index_new{:})=x(Index_new{:}).*(temp_cumu(Index_End{:})-cumu_start);
        end

       y=dzdy./L.^(beta)-2*alpha*beta.*LL;
   end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function [soln,eqn,info] = Poisson3CR(node,elem,bdFlag,pde,option,varargin)
%% POISSON3CR Poisson equation: Crouzeix-Raviart element in 3-D.
%
% u = POISSON3CR(node,elem,bdFlag,pde,option) produces the Crouzeix and Raviart
%   linear nonconforming  finite element approximation of the Poisson
%   equation in 3-D
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
% [soln,eqn,info] = POISSON3CR(node,elem,pde,bdFlag,option)
%
% The usage is the same as Poisson3. CR element on a tetrahedron is
% summarized in Poisson3CRfemrate for detail.
%
% Example
%
%   Poisson3CRfemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('bdFlag','var'), bdFlag = []; end
if ~exist('option','var'), option = []; end
NT = size(elem,1);

%% Diffusion coefficient
t = cputime;  % record assembling time
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                   % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)
    [lambda,weight] = quadpts3(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
	   pxyz = lambda(p,1)*node(elem(:,1),:) ...
		 	+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:) ...
            + lambda(p,4)*node(elem(:,4),:);
        K = K + weight(p)*pde.d(pxyz);           % d is a function
   end
end

%% Construct data structure
[elem2face,face] = dof3face(elem);
[Dphi,volume] = gradbasis3(node,elem);
NF = size(face,1);
Ndof = NF;

%% Assemble stiffness matrix
A = sparse(Ndof,Ndof);
for i = 1:4
    for j = i:4
        % local to global index map
		ii = double(elem2face(:,i));
		jj = double(elem2face(:,j));
        % local stiffness matrix.
        Aij = 9*dot(Dphi(:,:,i),Dphi(:,:,j),2).*volume;
        if ~isempty(pde.d)
            Aij = K.*Aij;
        end
        if(j==i)
         A = A + sparse(ii,jj,Aij,Ndof,Ndof);
        else
         A = A + sparse([ii;jj],[jj;ii],[Aij; Aij],Ndof,Ndof);
        end
    end
end
clear K Aij

%% Assemble the right hand side
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    [lambda,weight] = quadpts3(option.fquadorder);
    phi = 1-3*lambda;                 % linear bases
	nQuad = size(lambda,1);
    bt = zeros(NT,4);
    for p = 1:nQuad
		% quadrature points in the x-y-z coordinate
		pxyz = lambda(p,1)*node(elem(:,1),:) ...
			 + lambda(p,2)*node(elem(:,2),:) ...
			 + lambda(p,3)*node(elem(:,3),:)...
             + lambda(p,4)*node(elem(:,4),:);
		fp = pde.f(pxyz);
        for i = 1:4
            bt(:,i) = bt(:,i) + weight(p)*phi(p,i)*fp;
        end
    end
    bt = bt.*repmat(volume,1,4);
    b = accumarray(elem2face(:),bt(:),[Ndof 1]);
end
clear pxyz bt

%% Get boundary condition
[AD,b,u,freeFace,isPureNeumann] = getbd3CR(b);

%% Record assembling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeFace), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
if isPureNeumann
    option.solver = 'mg';
end
solver = option.solver;
% solve
switch solver
    case 'direct'
    tic;
    u(freeFace) = AD(freeFace,freeFace)\b(freeFace);
    residual = norm(b - AD*u);
    info = struct('solverTime',toc,'itStep',0,'error',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'error',[],'flag',3,'stopErr',[]);
    case 'mg'
        option.x0 = u;
        option.solver = 'CG';
        if nargin>=6
            HB = varargin{1};
        else
            HB = [];
        end
        [u,info] = mg(AD,b,elem,option,face,HB);
%         option.freeDof = freeFace;
%         [u(freeFace),info] = mg(AD(freeFace,freeFace),b(freeFace),elem,option,face,HB);
    case 'amg'
        option.solver = 'CG';
        [u(freeFace),info] = amg(AD(freeFace,freeFace),b(freeFace),option);
end
% post-process for pure Neumann problem
if isPureNeumann
    patchVolume = accumarray(elem2face(:),repmat(volume/4,4,1), [NF 1]);
    uc = sum(u.*patchVolume)/sum(volume);
    u = u - uc;   % int u = 0
end

%% Compute Du
dudx = u(elem2face(:,1)).*Dphi(:,1,1)+u(elem2face(:,2)).*Dphi(:,1,2) ...
     + u(elem2face(:,3)).*Dphi(:,1,3)+u(elem2face(:,4)).*Dphi(:,1,4);
dudy = u(elem2face(:,1)).*Dphi(:,2,1)+u(elem2face(:,2)).*Dphi(:,2,2) ...
     + u(elem2face(:,3)).*Dphi(:,2,3)+u(elem2face(:,4)).*Dphi(:,2,4);
dudz = u(elem2face(:,1)).*Dphi(:,3,1)+u(elem2face(:,2)).*Dphi(:,3,2) ...
     + u(elem2face(:,3)).*Dphi(:,3,3)+u(elem2face(:,4)).*Dphi(:,3,4);
Du = -3*[dudx, dudy, dudz];

%% Output
if nargout == 1
    soln = u;
else
    soln = struct('u',u,'Du',Du);
    eqn = struct('A',AD,'b',b,'face',face,'freeFace',freeFace,'Lap',A);
    info.assembleTime = assembleTime;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbd3CR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [AD,b,u,freeFace,isPureNeumann] = getbd3CR(b)
%% GETBD3CR Boundary conditions for Poisson equation: Crouzeix-Raviart element in 3D.

    u = zeros(Ndof,1);
    %% Initial check
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      % index of Robin faces in bdFlag
    if any(idxR)
        isRobin = false(Ndof,1);
        isRobin(elem2face(idxR)) = true;
        Robin = face(isRobin,:);  % Robin faces
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        v12 = node(Robin(:,2),:)-node(Robin(:,1),:);
        v13 = node(Robin(:,3),:)-node(Robin(:,1),:);
        area = 0.5*sqrt(abs(sum(mycross(v12,v13,2).^2,2)));
        center = (node(Robin(:,1),:) + node(Robin(:,2),:)+node(Robin(:,3),:))/3;
        ii = find(isRobin);
        ss = pde.g_R(center).*area; % exact for linear g_R
        A = A + sparse(ii,ii,ss,Ndof,Ndof);
    end

    % Dirichlet boundary nodes: fixedFace
    fixedFace = [];
    freeFace = [];
    if ~isempty(bdFlag) % find boundary faces
        idxD = (bdFlag(:)==1); % all Dirichlet faces in bdFlag
        isFixedFace = false(NF,1);
        isFixedFace(elem2face(idxD)) = true;
        fixedFace = find(isFixedFace);
        freeFace = ~isFixedFace;
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        % no bdFlag, only pde.g_D is given in the input
        s = accumarray(elem2face(:), 1, [NF 1]);
        fixedFace = find(s == 1);
        freeFace = (s == 2);
    end
    % Modify the matrix for different boundary conditions
    % Dirichlet boundary condition
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedFace,fixedFace)=I, AD(fixedFace,freeFace)=0, AD(freeFace,fixedFace)=0.
    if ~isempty(fixedFace)
        bdidx = zeros(Ndof,1);
        bdidx(fixedFace) = 1;
        Tbd = spdiags(bdidx,0,Ndof,Ndof);
        T = spdiags(1-bdidx,0,Ndof,Ndof);
        AD = T*A*T + Tbd;
    end
    % pure Neumann boundary condition
    isPureNeumann = false;
    if isempty(fixedFace) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        AD = A;
        AD(1,1) = AD(1,1) + 1e-6;
%         fixedFace = 1;
%         freeFace = 2:Ndof;    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Robin boundary condition
    if isempty(fixedFace) && ~isempty(Robin)
        AD = A;
    end

    %% Part 2: Find boundary faces and modify the right hand side b
    % Find boundary faces: Neumann
    Neumann = [];
    if ~isempty(bdFlag)  % bdFlag specifies different bd conditions
        idxN = (bdFlag(:) == 2);      % all Neumann faces in bdFlag
        isNeumann = elem2face(idxN | idxR); % index of Neumann faces
        Neumann = face(isNeumann,:);      % Neumann faces
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag is given, only pde.g_N or pde.g_R is given in the input
        s = accumarray(elem2face(:), 1, [NF 1]);
        isNeumann = (s==1);         % index of boundary faces
        Neumann = face(s == 1,:);   % boundary faces are Neumann faces
    end
    % Neumann boundary condition
    if ~isempty(Neumann) && ~isempty(pde.g_N) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 3;   % default order exact for linear gN
        end
        [lambdagN,weightgN] = quadpts(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        gf = zeros(size(Neumann,1),1);
        v12 = node(Neumann(:,2),:)-node(Neumann(:,1),:);
        v13 = node(Neumann(:,3),:)-node(Neumann(:,1),:);
        area = 0.5*sqrt(abs(sum(mycross(v12,v13,2).^2,2)));
        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:)...
                 + lambdagN(pp,3)*node(Neumann(:,3),:);
            gNp = pde.g_N(ppxy);
            gf = gf+ weightgN(pp)*gNp;
        end
        gf = gf.*area;
        b(isNeumann) = b(isNeumann) +gf;
    end
    % The case with non-empty Neumann faces but g_N=0 or g_N=[] corresponds to
    % the zero flux boundary condition on Neumann faces and no modification of
    % A,u,b is needed.

    % Dirichlet boundary condition
    if ~isPureNeumann && ~isempty(fixedFace) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && all(pde.g_D == 0))    % nonzero g_D
        if isnumeric(pde.g_D)
            u(fixedFace) = pde.g_D(fixedFace);
        else % pde.g_D is a function handle
            center = (node(face(fixedFace,1),:)+node(face(fixedFace,2),:)+node(face(fixedFace,3),:))/3;
            u(fixedFace) = pde.g_D(center);
        end
        b = b - A*u;
        b(fixedFace) = u(fixedFace);
    end
    % The case with non-empty Dirichlet nodes but g_D=0 or g_D=[] corresponds
    % to the zero Dirichlet boundary condition and no modification of u,b is
    % needed.

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b); % compatilbe condition: sum(b) = 0
%         b(1) = 0;
    end
    end % end of getbd3CR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % end of Poisson3CR

repo: ifem


text:
function [vw, OK] = meshLoad(vw, mshFileName, displayFlag)
% Load mesh from a file and add it to the 3D gray view.
%
%  [vw, OK] = meshLoad(vw, mshFileName, [displayFlag=0]);
%
% When we load the mesh, we also compute the vertexGrayMap on the fly.
% This defines the mapping from the mesh vertices into gray matter for
% that particular data set (defined in vw).
%
% If the mesh filename is provided as a directory, will prompt the user for
% a mesh file, starting in that directory. If it is set as 'firstmesh', it
% will load the first mesh file in the view's current mesh directory.
%
% If the optional 'displayFlag' arg is set to 1, it will  open a mrMesh
% window for the mesh. [by default, it's set to 0: don't display, just
% load.]
%
% ras, 11/07: added display flag.

if notDefined('vw'),            vw = getCurView;                      end
if notDefined('mshFileName'),   mshFileName = viewGet(vw, 'meshdir'); end
if notDefined('displayFlag'),	displayFlag = 0;					  end

OK = 1;

if ismember(lower(mshFileName), {'firstmesh' 'newest' 'mostrecent'})
	meshDir = viewGet(vw, 'MeshDir');
	w = dir( fullfile(meshDir, '*.mat') );
	[meshFiles, I]= setdiff({w.name}, {'MeshSettings.mat' 'MeshAngles.mat'});

	if isempty(meshFiles)
		error('No mesh files in the mesh directory: %s\n', meshDir);
	else
		if ismember(lower(mshFileName), {'newest' 'mostrecent'})
			% find most recent mesh file
			[dates, order] = sortrows( datevec({w(I).date}) ); %#ok<ASGLU>
			mshFileName = fullfile(meshDir, meshFiles{ order(end) });
		else
			% find first in list
			mshFileName = fullfile(meshDir, meshFiles{1});
		end
	end
end

msh = mrmReadMeshFile(mshFileName);
if(isempty(msh)), OK = 0; return; end

% Updates the mesh mesh to a new format directly readable by the mex files
% GB 02/13/05
msh = meshFormat(msh);

% Create the appropriate map from mesh vertices to gray matter.  This can
% differ from session to session because we acquire different slices.
%
if ispref('VISTA', 'autoComputeV2GMap')
	autoComputeV2G = getpref('VISTA', 'autoComputeV2GMap');
else
	autoComputeV2G = false;
end

if isempty(meshGet(msh,'vertexGrayMap')) || autoComputeV2G
    vertexGrayMap = mrmMapVerticesToGray(...
        meshGet(msh, 'initialvertices'),...
        viewGet(vw, 'nodes'),...
        viewGet(vw, 'mmPerVox'),...
        viewGet(vw, 'edges'));

    msh = meshSet(msh, 'vertexgraymap', vertexGrayMap);
end

%% Also compute the grayVertexMap if it is a requested preference
if(ispref('VISTA','autoComputeG2VMap'))
    if (getpref('VISTA','autoComputeG2VMap'))
        disp('Computing g2vmap for mesh');
        if isempty(meshGet(msh,'grayVertexMap'))
            grayVertexMap = mrmMapGrayToVertices(...
                viewGet(vw, 'nodes'),...
                meshGet(msh, 'initialvertices'),...
                viewGet(vw, 'mmPerVox'));
            msh = meshSet(msh, 'grayvertexmap', grayVertexMap);
        end
    end
end

msh = meshSet(msh,'id',-1);

if(~isfield(msh, 'fibers')), msh = meshSet(msh, 'fibers', []); end

if displayFlag==1
	msh = meshVisualize(msh);
	mrmSet(msh, 'cursoroff');  % a preference ... feel free to change back
end

vw = viewSet(vw, 'addandselectmesh', msh);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function blkStruct = slblocks
%SLBLOCKS Defines the block library for a specific Toolbox or Blockset.

% Mark Beale, 11-31-97
% Copyright 1992-2002 The MathWorks, Inc.
% $Revision: 1.12 $

blkStruct.Name = sprintf('voderiv');
blkStruct.OpenFcn = '';
blkStruct.IsFlat  = 0;% Is this library "flat" (i.e. no subsystems)?

% blkStruct.MaskDisplay = '';

Browser(1).Library = 'voderiv';
Browser(1).Name    = 'Fractional Variable Order Derivative Toolkit';
Browser(1).IsFlat  = 1;

blkStruct.Browser = Browser;

author: Sable

repo: mcbench-benchmarks


text:
function d = DownDyadHi(x,qmf)
% DownDyadHi -- Hi-Pass Downsampling operator (periodized)
%  Usage
%    d = DownDyadHi(x,f)
%  Inputs
%    x    1-d signal at fine scale
%    f    filter
%  Outputs
%    y    1-d signal at coarse scale
%
%  See Also
%    DownDyadLo, UpDyadHi, UpDyadLo, FWT_PO, iconv
%
	d = iconv( MirrorFilt(qmf),lshift(x));
	n = length(d);
	d = d(1:2:(n-1));

%
% Copyright (c) 1993. Iain M. Johnstone
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner


text:
function [logval val]=iwdensity(sigma,S,alpha,n)

% function function [logval val]=iwdensity(sigma,S,alpha,n)
% computes the density (value and log value) of the inverse Wishart distribution
% inputs:  - matrix 'sigma': 'true' variance-covariance matrix of VAR residuals, for the original Minnesota prior
%          - matrix 'S': scale matrix for sigma
%          - integer 'alpha': degrees of freedom for sigma
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
% outputs: - scalar 'logval': log of density value
%          - scalar 'val': density value

% calculates the density (a.2.10)

temp(1,1)=-(alpha*n/2)*log(2);
temp(2,1)=-bear.mgamma(alpha/2,n);
temp(3,1)=(alpha/2)*log(det(S));
temp(4,1)=-((alpha+n+1)/2)*log(det(sigma));
temp(5,1)=-0.5*trace(sigma\S);

logval=sum(temp);
val=exp(logval);

author: european-central-bank

repo: BEAR-toolbox


text:
function status = iscolvec(x)

status = ~isempty(x) && isvector(x) && size(x,2) == 1;

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef UF9 < PROBLEM
% <multi> <real> <large/none>
% Unconstrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 3;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = [0,0,zeros(1,obj.D-2)-2];
            obj.upper    = [1,1,zeros(1,obj.D-2)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            J1 = 4 : 3 : obj.D;
            J2 = 5 : 3 : obj.D;
            J3 = 3 : 3 : obj.D;
            Y  = X - 2*repmat(X(:,2),1,obj.D).*sin(2*pi*repmat(X(:,1),1,obj.D)+repmat(1:obj.D,size(X,1),1)*pi/obj.D);
            PopObj(:,1) = 0.5*(max(0,1.1*(1-4*(2*X(:,1)-1).^2))+2*X(:,1)).*X(:,2)   + 2*mean(Y(:,J1).^2,2);
            PopObj(:,2) = 0.5*(max(0,1.1*(1-4*(2*X(:,1)-1).^2))-2*X(:,1)+2).*X(:,2) + 2*mean(Y(:,J2).^2,2);
            PopObj(:,3) = 1-X(:,2)                                                  + 2*mean(Y(:,J3).^2,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,3);
            R(R(:,1)>(1-R(:,3))/4 & R(:,1)<(1-R(:,3))*3/4,:) = [];
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            a = linspace(0,1,20)';
            x = a*a';
            y = a*(1-a');
            z = (1-a)*ones(size(a'));
            z(x>(1-z)/4 & x<(1-z)*3/4) = nan;
            R = {x,y,z};
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function inputemu(varargin)
%INPUTEMU   Java-Based Mouse/Keyboard Emulator
%   INPUTEMU(ACTION,PARAM) emulates the human input via mouse and
%   keyboard. This utility uses Java java.awt.Robot class.
%
%   INPUTEMU(MOUSE_ACTION,[X Y WHEEL]) performs mouse emulation with 5
%   action types:
%
%      'move'|'none'|'wheel' - No mouse click (default)
%      'normal'              - Click left mouse button
%      'extend'              - Shift-click left mouse button
%      'alternate'           - Control-click left mouse button
%      'open'                - Double click any mouse button
%
%   The mouse cursor is first moved to the coordinate [X Y] (with respect
%   to the lower-left corner of the primary screen), followed by button
%   click action (if specified), then it turns the mouse wheel by WHEEL
%   notches. The number of "notches" to move the mouse wheel negative
%   values indicate movement up/away from the user, positive values
%   indicate movement down/towards the user. If PARAM (or any of its
%   elements) is missing, the default position is the current cursor
%   position and WHEEL = 0.
%
%   INPUTEMU('wheel',WHEEL) may be used to turn wheel without moving mouse
%   cursor.
%
%   INPUTEMU(KEYBOARD_ACTION,'text') performs keyboard emulation with
%   4 action types:
%
%      'key_normal'  - Normal keypress (shift-key pressed as needed)
%      'key_ctrl'    - CONTROL-key keypress
%      'key_alt'     - ALT-key keypress
%      'key_win'     - WIN-key keypress
%
%   The 'text' to be typed may contain special keys as escape characters
%   with '\' prefix:
%
%      '\BACKSPACE'   '\TAB'         '\ENTER'       '\SHIFT'
%      '\CTRL'        '\ALT'         '\PAUSE'       '\CAPSLOCK'
%      '\ESC'         '\PAGEUP'      '\PAGEDOWN'    '\END'
%      '\HOME'        '\LEFT'        '\UP'          '\RIGHT'
%      '\DOWN'        '\PRINTSCREEN' '\INSERT'      '\DELETE'
%      '\WINDOWS'     '\NUMLOCK'     '\SCROLLLOCK'
%
%   These escape characters are NOT case sensitive while regular characters
%   are. For regular backslash, use '\\' unless it is the only character
%   then '\' may be used.
%
%   In addition to above action types, there are 8 low-level actions to
%   specify button or key to be down (pressed) or up (released).
%
%      'left_down'/'left_up' | 'drag_on'/'drag_off' - Left mouse button
%      'right_down'/'right_up'                      - Right mouse button
%      'middle_down'/'middle_up'                    - Middle mouse button
%      'key_down'/'key_up'                          - Keyboard key
%
%   INPUTEMU(CMDS) can be used to perform multiple commands. CMDS is an
%   2-by-N cell array where the n-th column {ACTION_n;PARAM_n} defines the
%   n-th input action. CMDS may be given in an N-by-2 array (if N~=2)
%
%   INPUTEMU(CMDS,T) performs the CMDS action sequence with T second
%   update period.
%
%   INPUTEMU(CMDS,Tints) where time is a N-element vector, specifies
%   update interval (in s) for individual action. time(1) specifies the lag
%   before performing the first action.
%
%   See also jmouseemu and inputlog.

% Version 1.0 (Aug. 31, 2010)
% Written by: Takeshi Ikuma
% Created: Aug. 31, 2010
% Revision History:
%  - (Aug. 31, 2010) : initial release

% To Do List
% * Support for pause action
% * Allow a choice between block/non-block modes
% * Internationalization
%   - support for user specified locale / OemChars def
%   - support for special keys not on the US keyboard
%   - support for non-European alphabet support (possible?)

% Java environment check
error(javachk('jvm'));
error(javachk('awt'));

% default (i.e., US|us locale) OEM character map
import java.awt.event.KeyEvent
defaultOemChars = {... % {'ch w/o shift' 'ch w/shift' Java virtual keycode}
   '`' '~'   KeyEvent.VK_BACK_QUOTE
   '1' '!'   KeyEvent.VK_1
   '2' '@'   KeyEvent.VK_2
   '3' '#'   KeyEvent.VK_3
   '4' '$'   KeyEvent.VK_4
   '5' '%'   KeyEvent.VK_5
   '6' '^'   KeyEvent.VK_6
   '7' '&'   KeyEvent.VK_7
   '8' '*'   KeyEvent.VK_8
   '9' '('   KeyEvent.VK_9
   '0' ')'   KeyEvent.VK_0
   '-' '_'   KeyEvent.VK_MINUS
   '=' '+'   KeyEvent.VK_EQUALS
   '[' '{'   KeyEvent.VK_OPEN_BRACKET
   ']' '}'   KeyEvent.VK_CLOSE_BRACKET
   '\' '|'   KeyEvent.VK_BACK_SLASH
   ';' ':'   KeyEvent.VK_SEMICOLON
   '''' '"'  KeyEvent.VK_QUOTE
   ',' '<'   KeyEvent.VK_COMMA
   '.' '>'   KeyEvent.VK_PERIOD
   '/' '?'   KeyEvent.VK_SLASH
};

% parse input arguments
[cmds,deltaT,idx_mouse,oemChars,msg] = parseInput(defaultOemChars,nargin,varargin);
if ~isempty(msg), error(msg); end

% expand upper-level actions
[actions,params,dt,msg] = convertToLowLevelActions(cmds,deltaT,idx_mouse,oemChars);
if ~isempty(msg), error(msg); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

import java.awt.Robot
import java.awt.event.InputEvent

% initialize the java.robot
robot = Robot();

for n = 1:size(actions)
	% cue the action
	switch actions{n}
		case 'move'
			robot.mouseMove(params{n}(1),params{n}(2));
		case 'wheel'
			robot.mouseWheel(params{n});
		case 'left_down'
			robot.mousePress(InputEvent.BUTTON1_MASK);
		case 'left_up'
			robot.mouseRelease(InputEvent.BUTTON1_MASK);
		case 'right_down'
			robot.mousePress(InputEvent.BUTTON3_MASK);
		case 'right_up'
			robot.mouseRelease(InputEvent.BUTTON3_MASK);
		case 'middle_down'
			robot.mousePress(InputEvent.BUTTON2_MASK);
		case 'middle_up'
			robot.mouseRelease(InputEvent.BUTTON2_MASK);
		case 'key_down'
			robot.keyPress(params{n});
		case 'key_up'
			robot.keyRelease(params{n});
	end

	% delay
	robot.delay(dt(n));
end

% wait until robot's done
robot.waitForIdle();

end % jinputemu

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse Input Arguments
function [cmds,deltaT,idx_mouse,oemChars,msg] = parseInput(defaultOemChars,n,argin)

cmds = {};
deltaT = [];
idx_mouse = [];
oemChars = defaultOemChars;

msg = nargchk(1,2,n); %#ok
if ~isempty(msg), return; end

if ~iscell(argin{1}) % single-command mode
	switch n
		case 1
			cmds = [argin(1);cell(1)];
		case 2
			cmds = argin;
		otherwise
			cmds = {}; % never gets here
	end
else % multiple-command mode
	cmds = argin{1};
	if n>1
		deltaT = argin{2};
	end
end

% get dimension
dim = size(cmds);

% use the first size-2 dimension as the command
idx = find(dim==2,1,'first');
if ~isempty(idx) && idx~=1
	idx = [idx 1:idx-1 idx+1:ndims(cmds)];
	cmds = permute(cmds,idx);
	dim(:) = dim(idx);
end

% check cmds
if dim(1)~=2
	msg = 'One of the sizes of the CMDS matrix must must be 2.';
	return;
end

mouse_actions = {'none','move','wheel','normal','extend','alternate',...
	'open','drag_on','drag_off','left_down','left_up','middle_down',...
	'middle_up','right_down','right_up'};
keyboard_actions = {'key_normal','key_alt','key_ctrl','key_win','key_up','key_down'};

if any(cellfun(@(action)~isempty(action) && ~ischar(action),cmds(1,:)))
	msg = 'Command actions (CMDS{1,:}) must be a string of characters.';
	return;
end

idx_mouse = cellfun(@(action)any(strcmpi(action,mouse_actions)),cmds(1,:));
idx_keyboard = cellfun(@(action)any(strcmpi(action,keyboard_actions)),cmds(1,:));

if any(~any(idx_mouse|idx_keyboard))
	msg = 'At least one command action is invalid.';
	return;
end

if any(cellfun(@(x)~isempty(x) && (~isnumeric(x) || ~any(numel(x)==[1 2 3]) || any(isinf(x)) || any(isnan(x))),cmds(2,idx_mouse)))
	msg = 'Position must be a 2- or 3-element vector with finite values';
	return;
end

if any(cellfun(@(x)~ischar(x),cmds(2,idx_keyboard)))
	msg = 'Text must be a character string';
	return;
end

% format mouse parameters
scrnsz = get(0,'ScreenSize');
ptrloc = get(0,'PointerLocation');
for n = find(idx_mouse)
	switch numel(cmds{2,n})
		case 0
			cmds{2,n} = [ptrloc 0];
		case 1
			if strcmpi(cmds{1,n},'wheel')
				cmds{2,n} = [ptrloc cmds{2,n}];
			else
				cmds{2,n} = [cmds{2,n} ptrloc(2) 0];
			end
		case 2
			cmds{2,n} = [cmds{2,n}(:)' 0];
	end
	if all(cmds{2,n}(1:2)==ptrloc)
		cmds{2,n}(1:2) = nan; % indicating no change
	else
      % move origin from lower-left corner to upper-left corner
      cmds{2,n}(2) = scrnsz(4)-cmds{2,n}(2);

      ptrloc = cmds{2,n}(2);  % update current mouse location
	end

end

% check time
if isempty(deltaT)
	deltaT = zeros([1 dim(2:end)]);
elseif any(numel(deltaT)~=1 || ~isnumeric(deltaT) || deltaT<0 || isnan(deltaT) || isinf(deltaT))
	msg = 'deltaT must be a non-negative scalar.';
elseif numel(deltaT)==1
	deltaT = repmat(round(deltaT*1000),[1 dim(2:end)]);
elseif numel(deltaT)~=numel(cmds)/2
	msg = 'deltaT must be a scalar or match the number of commands.';
else
   deltaT = round(deltaT*1000);
end

end % parseinput()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% expand all actions to include only 'XXX_down','XXX_up','move','wheel'
function [actions,params,dt,msg] = convertToLowLevelActions(cmds,deltaT,idx_mouse,oemChars)

import java.awt.event.KeyEvent

actions = {};
params = {};
dt = [];
msg = '';

Ncmds = numel(deltaT);

% determine total number of actions
Nactions = zeros(Ncmds,1);
for n = 1:Ncmds
	if idx_mouse(n) % mouse command
      % determine number of LL actions needed for each mouse command
		Nactions(n) = sum([...
			cmds{2,n}(1)>=0 || cmds{2,n}(2)>=0% move
			cmds{2,n}(3)~=0 % wheel
			~any(strcmpi(cmds{1,n},{'none','move','wheel'})) % click action
			strcmpi(cmds{1,n},'normal')   % both down & up
			3*any(strcmpi(cmds{1,n},{'open','extend','alternate'}))]); % down & up & modifier key down & up

		if strcmpi(cmds{1,n},'drag_on')
			cmds{1,n} = 'left_down';
		elseif strcmpi(cmds{1,n},'drag_off')
			cmds{1,n} = 'left_up';
		end
	else % keyboard command
		[cmds{1,n},cmds{2,n},msg] = parseKeyboardText(cmds{1,n},cmds{2,n},oemChars);
		if ~isempty(msg), return; end
		Nactions(n) = size(cmds{1,n},1);
	end
end

% create actions & params
Ntotal = sum(Nactions);
actions = cell(Ntotal,1);
params = cell(Ntotal,1);
dt = zeros(Ntotal,1);
I = 1;
for n = 1:Ncmds
	dt(I) = deltaT(n);
	if idx_mouse(n)
		if ~isnan(cmds{2,n}(1))
			actions{I} = 'move';
			params{I} = cmds{2,n}([1 2]);
			I = I + 1;
		end
		switch lower(cmds{1,n})
			case 'normal'
				actions{I} = 'left_down';
				I = I + 1;
				actions{I} = 'left_up';
				I = I + 1;
			case 'extend'
				actions{I} = 'key_down';
				params{I} = KeyEvent.VK_SHIFT;
				I = I+1;
				actions{I} = 'left_down';
				I = I + 1;
				actions{I} = 'left_up';
				I = I + 1;
				actions{I} = 'key_up';
				params{I} = KeyEvent.VK_SHIFT;
				I = I+1;
			case 'alternate'
				actions{I} = 'key_down';
				params{I} = KeyEvent.VK_CONTROL;
				I = I+1;
				actions{I} = 'left_down';
				I = I + 1;
				actions{I} = 'left_up';
				I = I + 1;
				actions{I} = 'key_up';
				params{I} = KeyEvent.VK_CONTROL;
				I = I+1;
			case 'open'
				actions{I} = 'left_down';
				I = I + 1;
				actions{I} = 'left_up';
				I = I + 1;
				actions{I} = 'left_down';
				I = I + 1;
				actions{I} = 'left_up';
				I = I + 1;
			case {'left_down','left_up',...
					'middle_down','middle_up',...
					'right_down','right_up'} % single button action
				actions{I} = lower(cmds{1,n});
				I = I + 1;
		end
		if cmds{2,n}(3)~=0
			actions{I} = 'wheel';
			params{I} = cmds{2,n}(3);
         I = I + 1;
		end
	else % keyboard action
		idx = I:(I+Nactions(n)-1);
		actions(idx) = cmds{1,n};
		params(idx) = cmds{2,n};
		I = I + Nactions(n);
	end
end

end % convertToLowLevelActions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% expand keyboard action to a sequence of 'key_down' and 'key_up'
function	[actions,params,msg] = parseKeyboardText(cmd,text,oemChars)

import java.awt.event.KeyEvent

actions = {};
params = {};
msg = '';

% virtual key codes for escape keys
spKeys = [{'\\'      KeyEvent.VK_BACK_SLASH
		'\BACKSPACE'   KeyEvent.VK_BACK_SPACE
		'\TAB'         KeyEvent.VK_TAB
		'\ENTER'       KeyEvent.VK_ENTER
		'\SHIFT'       KeyEvent.VK_SHIFT
		'\CTRL'        KeyEvent.VK_CONTROL
		'\ALT'         KeyEvent.VK_ALT
		'\PAUSE'       KeyEvent.VK_PAUSE
		'\CAPSLOCK'    KeyEvent.VK_CAPS_LOCK
		'\ESC'         KeyEvent.VK_ESCAPE
		'\PAGEUP'      KeyEvent.VK_PAGE_UP
		'\PAGEDOWN'    KeyEvent.VK_PAGE_DOWN
		'\END'         KeyEvent.VK_END
		'\HOME'        KeyEvent.VK_HOME
		'\LEFT'        KeyEvent.VK_LEFT
		'\UP'          KeyEvent.VK_UP
		'\RIGHT'       KeyEvent.VK_RIGHT
		'\DOWN'        KeyEvent.VK_DOWN
      '\PRINTSCREEN' KeyEvent.VK_PRINTSCREEN
		'\INSERT'      KeyEvent.VK_INSERT
		'\DELETE'      KeyEvent.VK_DELETE
		'\WINDOWS'     KeyEvent.VK_WINDOWS
		'\NUMLOCK'     KeyEvent.VK_NUM_LOCK
		'\SCROLLLOCK'  KeyEvent.VK_SCROLL_LOCK}
      [strcat({'\F'},num2str((1:12)','%02d')) num2cell(KeyEvent.VK_F1+(0:11)')
      strcat({'\F'},num2str((13:24)','%02d')) num2cell(KeyEvent.VK_F13+(0:11)')]];

% scan the text and determine regular vs. special key and total number of characters
chCodes = zeros(numel(text),1);
idx = 0;

% ID type of each character in text:
%    regular char -> ASCII code
%    escape char -> negative of row index to spKeys followed by zeros
if numel(text)==1 % single character (to account for '\' case)
   chCodes = text;
else % string of characters
   while ~isempty(text)
      if text(1)=='\' % first character = special character
         tok = '';
         rest = text;
      else % first character = regular character
         [tok,rest] = strtok(text,'\');
      end
      chCodes(idx+(1:numel(tok))) = tok; % regular characters
      idx = idx+numel(tok);
      if isempty(rest), break; end

      I = find(cellfun(@(sch)strncmpi(rest,sch,numel(sch)),spKeys(:,1)));
      if isempty(I)
         msg = 'Invalid key.';
         return;
      end
      chCodes(idx+1) = -I;
      idx = idx + numel(spKeys{I});
      text = rest(numel(spKeys{I})+1:end);
   end
end

% remove the zero fillers for escape key string
chCodes = chCodes(chCodes~=0);

% total number of characters to type
Nchars = numel(chCodes);

% determine the virtual key codes for each key & need for shift key press
keycodes = zeros(Nchars,1);
shifts = false(Nchars,1);
I = chCodes>0;

%  - regular keys first
vkc = char(chCodes(I));
s = false(size(vkc));

%  * upper case
idx = vkc>='A' & vkc<='Z';
s(idx) = true;

%  * lower case
idx = vkc>='a' & vkc<='z';
vkc(idx) = vkc(idx)-'a'+'A';

%  * OEM keys
[idx1,J1] = ismember(vkc,oemChars(:,1)); % no shift
[idx2,J2] = ismember(vkc,oemChars(:,2)); % with shift
vkc(idx1) = cell2mat(oemChars(J1(idx1),3));
vkc(idx2) = cell2mat(oemChars(J2(idx2),3));
s(idx2) = true;

%  - update regular keys
keycodes(I) = vkc;
shifts(I) = s;

%  - then escape keys (no shift)
I = ~I;
keycodes(I) = cell2mat(spKeys(-chCodes(I),2));

%determine how to press keys
idx_type = strcmpi(cmd,{'key_normal','key_ctrl','key_alt','key_win'});
mod = any(idx_type(2:4)); % true to press modifier key during
press = ~strcmpi(cmd,'key_up'); % true to press key
release = ~strcmpi(cmd,'key_down'); % true to release key
shifts = [shifts(1);diff(shifts)]; % if 1, shift down, if -1 shift up (for each key)
release_shift = sum(shifts)>0; % if last character needs shift down, release it at the end

%determine total # of key down/up actions to perform
Nacts = 2*mod ... % if needs modiefier key, must press & release (2)
   + Nchars*(press+release) ... % for each character, needs 2 (if down&up) or 1 (if down or up)
   + sum(shifts~=0) ... % for each shift press/release actions
	+ release_shift; % to release the shift at the end

% define the keyboard action command cell
actions = cell(Nacts,1);
params = cell(Nacts,1);
I = 1; % action index
if mod % press the modifier key
	actions{I} = 'key_down';
	switch find(idx_type,1)
		case 2 %'key_ctrl'
			mod_param = KeyEvent.VK_CONTROL;
		case 3 %'key_alt'
			mod_param = KeyEvent.VK_ALT;
		case 4 %'key_win'
			mod_param = KeyEvent.VK_WINDOWS;
      otherwise % should never get here
			mod_param = [];
	end
	params{I} = mod_param;
	I = I + 1;
end
for n = 1:Nchars; % for each character
	if shifts(n)~=0 % shift key action needed
		if shifts(n)>0 % 1:press shift key
			actions{I} = 'key_down';
      else           % -1:release shift key
			actions{I} = 'key_up';
		end
		params{I} = KeyEvent.VK_SHIFT;
		I = I + 1;
	end
	if press
		actions{I} = 'key_down';
		params{I} = keycodes(n);
		I = I + 1;
	end
	if release
		actions{I} = 'key_up';
		params{I} = keycodes(n);
		I = I+1;
	end
end
if release_shift
	actions{I} = 'key_up';
	params{I} = KeyEvent.VK_SHIFT;
	I = I + 1;
end
if mod % release the modifier key
	actions{I} = 'key_up';
	params{I} = mod_param;
end

end % parseKeyboardText

% Copyright (c)2010, Takeshi Ikuma
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%   * Redistributions of source code must retain the above copyright
%   notice, this list of conditions and the following disclaimer. *
%   Redistributions in binary form must reproduce the above copyright
%   notice, this list of conditions and the following disclaimer in the
%   documentation and/or other materials provided with the distribution.
%   * Neither the names of its contributors may be used to endorse or
%   promote products derived from this software without specific prior
%   written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


text:
function test_failed=test_phaselock
%TEST_PHASELOCK  Test phaselock and phaseunlock
%

test_failed=0;

disp(' ===============  TEST_PHASELOCK ================');

% set up parameters
L=420;
f=tester_rand(L,1);
g=pgauss(L);
a=10;b=30;M=L/b;

c = dgt(f,g,a,M);
cp1 = ref_phaselock(c,a);
cp2 = phaselock(c,a,'lt',[0 1]);

% compare original phaselock with mine for rectangular case
res=norm(cp1-cp2,'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASELOCK REF RECT %0.5g %s\n'],res,fail);

% comparisons for non-separable case
c_big = dgt(f,g,a,2*M);
c_quin = dgt(f,g,a,M,'lt',[1 2]);

c_bigp = phaselock(c_big,a);
c_quinp= phaselock(c_quin,a,'lt',[1 2]);

% compare the quincunx lattice with twice transform on twice as many
% chanels
res=norm(c_bigp(1:2:end,1)-c_quinp(:,1),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASELOCK QUIN 1    %0.5g %s\n'],res,fail);

res=norm(c_bigp(2:2:end,2)-c_quinp(:,2),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASELOCK QUIN 2    %0.5g %s\n'],res,fail);

% testing of phaseunlock routine
res=norm(c_big - phaseunlock(c_bigp,a),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASEUNLOCK RECT    %0.5g %s\n'],res,fail);

res=norm(c_quin - phaseunlock(c_quinp,a,'lt',[1 2]),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASEUNLOCK QUIN    %0.5g %s\n'],res,fail);

cfi = dgtreal(f,g,a,M);
cti = dgtreal(f,g,a,M,'timeinv');

res=norm(cfi - phaseunlockreal(cti,a,M),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASEUNLOCKREAL     %0.5g %s\n'],res,fail);

res=norm(phaselockreal(cfi,a,M) - cti,'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASELOCKREAL       %0.5g %s\n'],res,fail);

res=norm(cfi - phaseunlockreal(cti,a,M,'precise'),'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASEUNLOCKREAL PREC %0.5g %s\n'],res,fail);

res=norm(phaselockreal(cfi,a,M,'precise') - cti,'fro');
[test_failed,fail]=ltfatdiditfail(res,test_failed);
fprintf(['PHASELOCKREAL PREC  %0.5g %s\n'],res,fail);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function P = optpower(epsil, n, R);
% compute the smallest required power to communicate at rate R with p.e. epsil and blocklength n
K = @(A) A./sqrt(2) .* sqrt(A.^2 + 2)./(A.^2 + 1) *log2(exp(1)) * norminv(epsil, 0, 1);

normapx = @(A) 1/2 * log2(1 + A.^2) + K(A)/sqrt(n) + 1/2 *log2(n)/n;
A0 = sqrt(2^(2*R) - 1);
A1 = fzero(@(A) normapx(A) - R, A0);
P = A1^2;

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre