text:
function [ x, y, seed ] = circle_segment_sample_from_height ( r, h, n, seed )

%*****************************************************************************80
%
%% CIRCLE_SEGMENT_SAMPLE_FROM_HEIGHT samples points from a circle segment.
%
%  Discussion:
%
%    Begin with a circle of radius R.  Choose two points P1 and P2 on the
%    circle, and draw the chord P1:P2.  This chord divides the circle
%    into two pieces, each of which is called a circle segment.
%    Consider one of the pieces.  The "angle" of this segment is the angle
%    P1:C:P2, where C is the center of the circle.  Let Q be the point on
%    the chord P1:P2 which is closest to C.  The "height" of the segment
%    is the distance from Q to the perimeter of the circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%    0 < R.
%
%    Input, real H, the height of the circle segment.
%    0 <= H <= 2 * R.
%
%    Input, integer N, the number of sample points.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real X(N,1), Y(N,1), the sample points.
%
  area = circle_segment_area_from_height ( r, h );
%
%  Pick CDF's randomly.
%
  if ( seed == 0 )
    u = rand ( n, 1 );
  else
    [ u, seed ] = r8vec_uniform_01 ( n, seed );
  end
%
%  Choose points randomly by choosing ordered areas randomly.
%
  area2(1:n,1) = u(1:n,1) * area;
%
%  Each area corresponds to a height H2.  Find it.
%
  h2 = zeros ( n, 1 );
  for i = 1 : n
    h2(i,1) = circle_segment_height_from_area ( r, area2(i,1) );
  end
%
%  Determine the half-width WH of the segment for each H2.
%
  wh(1:n,1) = sqrt ( h2(1:n,1) .* ( 2.0 * r - h2(1:n,1) ) );
%
%  Choose an X position randomly in [-WH,+WH].
%
  if ( seed == 0 )
    u = rand ( n, 1 );
  else
    [ u, seed ] = r8vec_uniform_01 ( n, seed );
  end

  x(1:n,1) = ( 2.0 * u(1:n,1) - 1.0 ) .* wh(1:n,1);
%
%  Our circle center is at (0,0).  Our height of H2 is subtracted
%  from the height R at the peak of the circle.  Determine the Y
%  coordinate using this fact.
%
  y(1:n,1) = r - h2(1:n,1);

  return
end


text:
function md_fast ( nd, np, step_num, dt )

%*****************************************************************************80
%
%% MD_FAST is the main program for the molecular dynamics simulation.
%
%  Discussion:
%
%    MD_FAST is a modification of the MD code.
%
%    The compute() function used in the original version of MD has been
%    rewritten to take advantage of some of MATLAB's vector capabilities.
%    It is hoped that MD_FAST will execute faster than MD does.
%
%    Like MD, MD_FAST implements a simple molecular dynamics simulation.
%
%    The velocity Verlet time integration scheme is used.
%
%    The particles interact with a central pair potential.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2014
%
%  Author:
%
%    Original FORTRAN90 version by Bill Magro.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer ND, the spatial dimension.
%    A value of 2 or 3 is usual.
%    The default value is 3.
%
%    Input, integer NP, the number of particles.
%    A value of 1000 or 2000 is small but "reasonable".
%    The default value is 1000.
%
%    Input, integer STEP_NUM, the number of time steps.
%    A value of 500 is a small but reasonable value.
%    The default value is 100.
%
%    Input, real DT, the time step.
%    A value of 0.1 is large; the system will begin to move quickly but the
%    results will be less accurate.
%    A value of 0.0001 is small, but the results will be more accurate.
%    The default value is 0.1.
%
  if ( nargin < 1 )
    nd = 3;
  end

  if ( nargin < 2 )
    np = 1000;
  end

  if ( nargin < 3 )
    step_num = 100;
  end

  if ( nargin < 4 )
    dt = 0.1;
  end

  mass = 1.0;
%
%  Set SHOW_PLOTS to 1 if you want to see a plot each time.
%
  show_plots = 0;
%
%  Report to the user.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MD_FAST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  A molecular dynamics program.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ND, the spatial dimension, is %d\n', nd );
  fprintf ( 1, '  NP, the number of particles in the simulation is %d.\n', np );
  fprintf ( 1, '  STEP_NUM, the number of time steps, is %d.\n', step_num );
  fprintf ( 1, '  DT, the time step size, is %f seconds.\n', dt );
%
%  Set the dimensions of the box.
%
  box(1:nd) = 10.0;
%
%  Initialize the random number generator to a prescribed state.
%  RNG is the new recommended interface to the MATLAB random number generator.
%
  seed = 123456789;
  rng ( seed );
%
%  Set initial positions, velocities, and accelerations.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initializing positions, velocities, and accelerations.\n' );

  [ pos, vel, acc ] = initialize ( np, nd, box, seed );
%
%  Compute the forces and energies.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computing initial forces and energies.\n' );

  [ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );
%
%  Save the initial total energy for use in the accuracy check.
%
  e0 = potential + kinetic;
%
%  This is the main time stepping loop:
%    Compute forces and energies,
%    Update positions, velocities, accelerations.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  At each step, we report the potential and kinetic energies.\n' );
  fprintf ( 1, '  The sum of these energies should be a constant.\n' );
  fprintf ( 1, '  As an accuracy check, we also print the relative error\n' );
  fprintf ( 1, '  in the total energy.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Step      Potential       Kinetic        (P+K-E0)/E0\n' );
  fprintf ( 1, '                Energy P        Energy K       Relative Energy Error\n' );
  fprintf ( 1, '\n' );

  step_print_index = 0;
  step_print_num = 10;

  step = 0;
  fprintf ( 1, '  %8d  %14f  %14f  %14e\n', ...
    step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
  step_print_index = step_print_index + 1;
  step_print = floor ( ( step_print_index * step_num ) / step_print_num );

  tic;

  for step = 1 : step_num

    [ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );

    if ( step == step_print )
      fprintf ( 1, '  %8d  %14f  %14f  %14e\n', ...
        step, potential, kinetic, ( potential + kinetic - e0 ) / e0 );
      step_print_index = step_print_index + 1;
      step_print = floor ( ( step_print_index * step_num ) / step_print_num );
    end
%
%  Plot the particles.
%
    if ( show_plots )

      if ( nd == 2 || nd == 3 )
        clf
        if ( nd == 2 )
          scatter ( pos(1,:), pos(2,:), 10.0, 'r', 'filled' )
        elseif ( nd == 3 )
          scatter3 ( pos(1,:), pos(2,:), pos(3,:), 10.0, 'r', 'filled' )
        end
        title ( sprintf ( 'Step %d\n', step ) );
        pause ( 0.5 )

      end

    end
%
%  Update the data for the next time step.
%
    [ pos, vel, acc ] = update ( np, nd, pos, vel, force, acc, mass, dt );

  end

  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Main computation:\n' );
  fprintf ( 1, '    Wall clock time = %f seconds.\n', wtime );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MD_FAST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ f, pot, kin ] = compute ( np, nd, pos, vel, mass )

%*****************************************************************************80
%
%% COMPUTE computes the forces and energies.
%
%  Discussion:
%
%    This is a modified version of the original compute90 function.  It has
%    been rewritten in an attempt to speed it up, by taking advantage of
%    MATLAB's vector capabilities.
%
%    The potential function V(X) is a harmonic well which smoothly
%    saturates to a maximum value at PI/2:
%
%      v(x) = ( sin ( min ( x, PI2 ) ) )^2
%
%    The derivative of the potential is:
%
%      dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos ( min ( x, PI2 ) )
%            = sin ( 2.0 * min ( x, PI2 ) )
%
%    Thanks to Abouzar Moshfegh for suggesting a revision to the kinetic
%    energy computation, 05 July 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2012
%
%  Author:
%
%    Original FORTRAN90 version by Bill Magro.
%    MATLAB version by John Burkardt and Gene Cliff.
%
%  Parameters:
%
%    Input, integer NP, the number of particles.
%
%    Input, integer ND, the number of spatial dimensions.
%
%    Input, real POS(ND,NP), the position of each particle.
%
%    Input, real VEL(ND,NP), the velocity of each particle.
%
%    Input, real MASS, the mass of each particle.
%
%    Output, real F(ND,NP), the forces.
%
%    Output, real POT, the total potential energy.
%
%    Output, real KIN, the total kinetic energy.
%
  f = zeros ( nd, np );

  pot = 0.0;

  pi2 = pi / 2.0;

  for i = 1 : np
%
%  Set up an array of vectors to particle I.
%
    Ri = pos - repmat ( pos( :, i ), 1, np );
%
%  Set up an array of distances.
%
    D = sqrt ( sum ( Ri.^2 ) );

    Ri = Ri( :, ( D > 0.0 ) );

    D = D( D > 0.0 );
%
%  Truncate the potential.
%
    D2 = D .* ( D <= pi2 ) + pi2 * ( D > pi2 );
%
%  Accumulate the potential energy.
%
    pot = pot + 0.5 * sum ( sin ( D2 ).^2 );
%
%  Compute force on particle I.
%
    f( :, i) = Ri * ( sin ( 2*D2 ) ./ D )';

  end
%
%  Compute the kinetic energy.
%
  kin = 0.5 * mass * sum ( sum ( vel.^2 ) );

  return
end
function [ pos, vel, acc, seed ] = initialize ( np, nd, box, seed )

%*****************************************************************************80
%
%% INITIALIZE initializes the positions, velocities, and accelerations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 July 2008
%
%  Author:
%
%    Original FORTRAN90 version by Bill Magro.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer NP, the number of particles.
%
%    Input, integer ND, the number of spatial dimensions.
%
%    Input, real BOX(ND), specifies the maximum position
%    of particles in each dimension.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real POS(ND,NP), the position of each particle.
%
%    Output, real VEL(ND,NP), the velocity of each particle.
%
%    Output, real ACC(ND,NP), the acceleration of each particle.
%
%    Output, integer SEED, a seed for the random number generator.
%

%
%  Start by setting the positions to random numbers between 0 and 1.
%
  pos(1:nd,1:np) = rand ( nd, np );
%
%  Use these random values as scale factors to pick random locations
%  inside the box.
%
  for i = 1 : nd
    pos(i,1:np) = box(i) * pos(i,1:np);
  end
%
%  Velocities and accelerations.
%
  vel(1:nd,1:np) = randn ( nd, np );
  acc(1:nd,1:np) = 0.0;

  return
end
function timestamp ( )

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

  return
end
function [ pos, vel, acc ] = update ( np, nd, pos, vel, f, acc, mass, dt )

%*****************************************************************************80
%
%% UPDATE updates positions, velocities and accelerations.
%
%  Discussion:
%
%    The time integration is fully parallel.
%
%    A velocity Verlet algorithm is used for the updating.
%
%    x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
%    v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
%    a(t+dt) = f(t) / m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 July 2008
%
%  Author:
%
%    Original FORTRAN90 version by Bill Magro.
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer NP, the number of particles.
%
%    Input, integer ND, the number of spatial dimensions.
%
%    Input, real POS(ND,NP), the position of each particle.
%
%    Input, real VEL(ND,NP), the velocity of each particle.
%
%    Input, real F(ND,NP), the force on each particle.
%
%    Input, real ACC(ND,NP), the acceleration of each
%    particle.
%
%    Input, real MASS, the mass of each particle.
%
%    Input, real DT, the time step.
%
%    Output, real POS(ND,NP), the updated position of each particle.
%
%    Output, real VEL(ND,NP), the updated velocity of each particle.
%
%    Output, real ACC(ND,NP), the updated acceleration of each
%    particle.
%
  rmass = 1.0 / mass;

  pos(1:nd,1:np) = pos(1:nd,1:np) + vel(1:nd,1:np) * dt ...
    + 0.5 * acc(1:nd,1:np) * dt * dt;

  vel(1:nd,1:np) = vel(1:nd,1:np) ...
    + 0.5 * dt * ( f(1:nd,1:np) * rmass + acc(1:nd,1:np) );

  acc(1:nd,1:np) = f(1:nd,1:np) * rmass;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Demo file for StructBrowser
%
% H.Lahdili,(hassan.lahdili@crc.ca)
% Communications Research Centre (CRC) | Advanced Audio Systems (AAS)
% www.crc.ca | www.crc.ca/aas
% Ottawa. Canada
%
% 30/04/2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load example
StructBrowser

author: Sable

repo: mcbench-benchmarks


text:
16 May 2004 immoments added
            ploto     added
            featsel routines with crossval implemented
25 Jun 2004 fdsc      added
30 Aug 2004 proxm     homogeneous measure added (x*y)^p
            disperror, confmat: writing to files enabled
26 Sep 2004 gaussm, emclust: get right number of components
            crisp / soft label handling for these routines
30 Sep 2004 plotr     allow multiple error curves
            scalem    variance option corrected
21 Oct 2004 svc_nu, svo_nu, added
 6 Dec 2004 parzendc, parzen_map: soft labels
            nmc: bug fix, allow soft labels
            prtest: bug fix, label generation
12 Dec 2004 plotr: cell array of error structures
14 Dec 2004 setlabels, roc, normm: bug fixes to correct / prevent
            reordering of lablist
20 Jan 2005 prwarning level now default 1
            some warning levels adapted
            in particular for empty priors
 1 Feb 2005 clevalf, debugged from reusing the same testset
24 Apr 2005 labcmp added
29 Oct 2005 trying to avoid prior warnings

----------------------------

 8 Mar 2006 start datafile options
26 Sep 2006 multiple labeling, see multi-labeling
26 Sep 2006 structure field handling for user and ident field in datasets
28 Nov 2006 testauc added

18 Dec 2006 global variables like PRMEMORY, PRTRACE, GRIDSIZE, etc
            are made persistent instead of global
19 Dec 2006 plotr renamed in plote (make space for a regression plot command)
            svc svc_nu: general option for kernel computation
22 Mar 2007 datafile construct ready
19 Apr 2007 Normal_map: inversion at learning
            Gaussm, mogs, improved: noise, more trials, regularisation
17 May 2007 NAIVEBC generalised with PARZENM and GAUSSM
18 May 2007 Streamlining PRPROGRESS and CLOSEMESS
            NMC now based on normal distributions with identically shaped
            spherical classifiers. This gives a better soft output. Consequently
            the soft version of EMCLUST is improved as it uses NMC for
            initialisation.
 5 Jun 2007 TESTC extended with various criteria
            Automatic parameter optimization in some routines
            RSSCC and RBSVC added
 3 Jul 2008 PRPROGRESS partially replaced by PRWAITBAR
10 Jul 2008 ident reorganised
21 Jul 2008 createdatafile
18 Aug 2008 Image Database Retrieval GUI
 5 Jan 2009 im_patch, band2obj and bandsel added
- new trainable combiners: naivebcc and mlrc (thanks to Chunxia!)
- Leave One Set Out crossvalidation: loso
- testc accepts missing classes in case of no_priors
- implementation of multi-band time-signals as multiband 1-D images
- this script generates now proper p-code :-)
10 Mar 2009 DCSC added
18 Mar 2009 MODSELC added

April  2011 FDSC, MDSC, DISNORM added

author: marianux

repo: ecg-kit

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


text:
% [BOUNDX,BOUNDY] = LABEL_OPEN_EDGES(RGBI)
%
% Input: RGBI is a RGB image we want to label edges in.
% This program allows the user to input any number of open
% edges. Here is the interactive usage:
% label_open_edges commands:
% LEFT CLICK with the mouse to add a new point.
% RIGHT CLICK with the mouse to add a new point on an image edge.
% MIDDLE CLICK with the mouse to finish the curve.
% Press ENTER to finish the curve.
% Press <CONTROL> z to undo the last curve piece added.
% Press H for a list of commands.
% Press z to zoom out.
% Press Z to put in zoomin mode, then draw a rectangle to zoomin on.
% Press - to decrease the first edge threshold.
% Press = to increase the first edge threshold.
% Press _ to decrease the second edge threshold.
% Press + to increase the second edge threshold.
% Press t to see the image with/without the edges.
% Press <CONTROL> q to quit.

function varargout = label_open_edges(arg1)

if nargin ~= 1,
  error('LABEL_OPEN_EDGES: Usage: label_open_edges(rgbI)');
end;

if strcmp(arg1,'Move'),
  move;
  return;
elseif strcmp(arg1,'ButtonDown'),
  buttondown;
  return;
elseif strcmp(arg1,'Drag'),
  drag;
  return;
elseif strcmp(arg1,'ButtonUp'),
  buttonup;
  return;
elseif strcmp(arg1,'KeyPress'),
  keypress;
  return;
end;

fprintf('Label edges: Type ''H'' for a list of commands.\n');

global FIG_H AXIS_H IMAGE_H NROWS NCOLS;
global BOUNDX BOUNDY CURVEN LASTLINE_H CURVEPIECESTART;
global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH ISEDGESHOWN;

% rgb image
RGBI0 = im2double(arg1);
RGBI = RGBI0;
[NROWS,NCOLS,ncolors] = size(RGBI0);

[EDGEI,THRESH] = edge(rgb2gray(RGBI0),'canny');

% draw the edges on the image
RGBI = reshape(RGBI,[NROWS*NCOLS],ncolors);
RGBI(EDGEI,1) = 1; RGBI(EDGEI,2) = .75; RGBI(EDGEI,3) = 1;
RGBI = reshape(RGBI,[NROWS,NCOLS,ncolors]);

% draw the image
clf;
image(RGBI); axis('image'); hold on;

% find the closest edge to each pixel
[DISTEDGE,NNEDGEI] = bwdist(EDGEI);

% save which image we are looking at
FIG_H = gcf;
AXIS_H = gca;
IMAGE_H = findobj(AXIS_H,'type','image');

% we will use IMAGE_H's UserData as a flag to see when we are done
set(IMAGE_H,'UserData',1);

% Initialize the lines drawn
BOUNDX{1} = []; BOUNDY{1} = [];
CURVEN = 1;
LASTLINE_H = [];
CURVEPIECESTART = [];
ISEDGESHOWN = 1;

% wait for the mouse to be pushed
set(IMAGE_H,'ButtonDownFcn','label_open_edges(''ButtonDown'')');
set(FIG_H,'WindowButtonMotionFcn','label_open_edges(''Move'')');

% or for a key to be pressed
set(FIG_H,'KeyPressFcn','label_open_edges(''KeyPress'')');

waitfor(IMAGE_H,'UserData','Done');

varargout = {BOUNDX , BOUNDY};

clear global BOUNDX BOUNDY;

% mouse moved
function move

global FIG_H IMAGE_H AXIS_H;
global BOUNDX BOUNDY CURVEN LASTLINE_H;
global EDGEI RGBI DISTEDGEI NNEDGEI;
global EDGEPOINTER_H;

% find the mouse location
[imageHandle,x,y] = over_image(AXIS_H);
if imageHandle ~= IMAGE_H | imageHandle == 0, return; end;

% find the closest edge point
[yedge,xedge] = ind2sub(size(EDGEI),NNEDGEI(round(y),round(x)));
if isempty(EDGEPOINTER_H) | ~ishandle(EDGEPOINTER_H),
  EDGEPOINTER_H = plot(xedge,yedge,'b.','markersize',12,'hittest','off');
else,
  set(EDGEPOINTER_H,'XData',xedge,'YData',yedge);
end;

% mouse button pressed down
function buttondown

global FIG_H IMAGE_H AXIS_H;
global BOUNDX BOUNDY CURVEN LASTLINE_H CURVEPIECESTART;
global EDGEI RGBI DISTEDGEI NNEDGEI;

% which mouse button was pushed?
selection_type = get(FIG_H,'SelectionType');

if strcmp(selection_type,'normal'),
  b = 1;
elseif strcmp(selection_type,'extend'),
  b = 2;
elseif strcmp(selection_type,'alt'),
  b = 3;
end;

% left mouse button: input a new point
if b == 1,
  % find the point clicked
  [imageHandle,x,y] = over_image(AXIS_H);
  if imageHandle ~= IMAGE_H | imageHandle == 0, return; end;

  i = length(BOUNDX{CURVEN});
  % add the point to the current curve
  BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; x];
  BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; y];

  % plot the new point
  if isempty(LASTLINE_H),
    LASTLINE_H = plot(BOUNDX{CURVEN},BOUNDY{CURVEN},...
		      'r','erasemode','none','hittest','off');
    CURVEPIECESTART = 1;
  else,
    set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN});
    CURVEPIECESTART = [CURVEPIECESTART;i+1];
  end;

  set(FIG_H,'WindowButtonMotionFcn','label_open_edges(''Drag'')',...
	    'WindowButtonUpFcn','label_open_edges(''ButtonUp'')');
  set(IMAGE_H,'ButtonDownFcn','');

% right mouse button: input a new point on an edge
elseif b == 3,

  % input a point
  [imageHandle,x,y] = over_image(AXIS_H);
  if imageHandle ~= IMAGE_H | imageHandle == 0, return; end;

  i = length(BOUNDX{CURVEN});

  % find the closest edge point to the boundary point
  [yedge,xedge] = ind2sub(size(EDGEI),NNEDGEI(round(y),round(x)));

  % add the point to the current curve
  BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; xedge];
  BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; yedge];

  % plot the new point
  if isempty(LASTLINE_H),
    LASTLINE_H = plot(BOUNDX{CURVEN},BOUNDY{CURVEN},...
		      'r','erasemode','none','hittest','off');
    CURVEPIECESTART = 1;
  else,
    set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN});
    CURVEPIECESTART = [CURVEPIECESTART;i+1];
  end;

  set(FIG_H,'WindowButtonMotionFcn','label_open_edges(''Drag'')',...
	    'WindowButtonUpFcn','label_open_edges(''ButtonUp'')');
  set(IMAGE_H,'ButtonDownFcn','');

% Done entering the current curve
elseif b == 2,
  set(LASTLINE_H,'color',[.75,.75,1],'linewidth',4);
  CURVEN = CURVEN + 1;
  LASTLINE_H = [];
  BOUNDX{CURVEN} = [];
  BOUNDY{CURVEN} = [];
end;

% drag the mouse
function drag

global EDGEI RGBI DISTEDGEI NNEDGEI;
global FIG_H AXIS_H IMAGE_H BOUNDX BOUNDY CURVEN LASTLINE_H;

% which mouse button was pushed?
selection_type = get(FIG_H,'SelectionType');
if strcmp(selection_type,'normal'),
  b = 1;
elseif strcmp(selection_type,'extend'),
  b = 2;
elseif strcmp(selection_type,'alt'),
  b = 3;
end;

[imageHandle,x,y] = over_image(AXIS_H);

if imageHandle ~= 0 & imageHandle == IMAGE_H & (b == 1 | b == 3),
  % left mouse button: input a new point
  if b == 1,
    % set control point to the current mouse location
    BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; x];
    BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; y];
  elseif b == 3,
    % right mouse button: input a new point on an edge
    [yedge,xedge] = ind2sub(size(EDGEI),NNEDGEI(round(y),round(x)));
    BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; xedge];
    BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; yedge];
  end;
  set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN});
end;

% mouse button up
function buttonup

global EDGEI RGBI DISTEDGEI NNEDGEI;
global IMAGE_H AXIS_H FIG_H BOUNDX BOUNDY CURVEN LASTLINE_H;
selection_type = get(FIG_H,'SelectionType');
if strcmp(selection_type,'normal'),
  b = 1;
elseif strcmp(selection_type,'extend'),
  b = 2;
elseif strcmp(selection_type,'alt'),
  b = 3;
end;

[imageHandle,x,y] = over_image(AXIS_H);
if imageHandle ~= 0,
  if imageHandle == IMAGE_H,
    if b == 2, return;
    elseif b == 1,
      % set control point to the current mouse location
      BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; x];
      BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; y];
    elseif b == 3,
      % right mouse button: input a new point on an edge
      [yedge,xedge] = ind2sub(size(EDGEI),NNEDGEI(round(y),round(x)));
      BOUNDX{CURVEN} = [BOUNDX{CURVEN} ; xedge];
      BOUNDY{CURVEN} = [BOUNDY{CURVEN} ; yedge];
    end;
    set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN});
  end;
end;

set(FIG_H,'WindowButtonUpFcn','',...
	  'WindowButtonMotionFcn','label_open_edges(''Move'')');
set(IMAGE_H,'ButtonDownFcn','label_open_edges(''ButtonDown'')');

% key was pressed
function keypress

global IMAGE_H AXIS_H FIG_H;

CONTROLZ = 26;
CONTROLQ = 17;
ENTER = 13;

c = get(FIG_H,'CurrentCharacter');

if c == 'h' | c == 'H',
  le_help;
elseif c == CONTROLZ,
  undo;
elseif c == 'z',
  zoomout;
elseif c == 'Z',
  zoomin;
elseif c == '-',
  decreasethresh1;
elseif c == '=',
  increasethresh1;
elseif c == '_',
  decreasethresh2;
elseif c == '+',
  increasethresh2;
elseif c == CONTROLQ,
  finish;
elseif c == ENTER,
  closecurve;
elseif c == 't',
  toggleedgeimage;
else,
end;

function le_help

fprintf('label_open_edges commands:\n');
fprintf('LEFT CLICK with the mouse to add a new point.\n');
fprintf('RIGHT CLICK with the mouse to add a new point on an image edge.\n');
fprintf('MIDDLE CLICK with the mouse to finish the curve.\n');
fprintf('Press ENTER to close and then finish the curve.\n'),
fprintf('Press <CONTROL> z to undo the last curve piece added.\n');
fprintf('Press H for a list of commands.\n');
fprintf('Press z to zoom out.\n');
fprintf('Press Z to put in zoomin mode, then draw a rectangle to zoomin on.\n');
fprintf('Press - to decrease the first edge threshold.\n');
fprintf('Press = to increase the first edge threshold.\n');
fprintf('Press _ to decrease the second edge threshold.\n');
fprintf('Press + to increase the second edge threshold.\n');
fprintf('Press t to see the image with/without the edges.\n');
fprintf('Press <CONTROL> q to quit.\n');

% undo the last curve piece entered
function undo

global BOUNDX BOUNDY CURVEN CURVEPIECESTART LASTLINE_H;
if ~isempty(CURVEPIECESTART) & ~isempty(BOUNDX),
  BOUNDX{CURVEN} = BOUNDX{CURVEN}(1:CURVEPIECESTART(end)-1);
  BOUNDY{CURVEN} = BOUNDY{CURVEN}(1:CURVEPIECESTART(end)-1);
  CURVEPIECESTART = CURVEPIECESTART(1:end-1);
  set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN});
  drawnow;
elseif isempty(CURVEPIECESTART) & length(BOUNDX) >= 2,
  BOUNDX = BOUNDX(1:CURVEN-1);
  BOUNDY = BOUNDY(1:CURVEN-1);
else,
  fprintf('Cannot undo -- no information left.\n');
end;

% zoom out so that the whole picture is in view
function zoomout

global AXIS_H NROWS NCOLS;

axes(AXIS_H);
axis([.5,NCOLS+.5,.5,NROWS+.5]);

% zoom in on a rectangle
function zoomin

global AXIS_H;
p = getrect(AXIS_H);
axis([p(1)-.5,p(3)+p(1)+.5,p(2)-.5,p(4)+p(2)+.5]);

function decreasethresh2

global IMAGE_H NROWS NCOLS;
global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH;

THRESH(2) = THRESH(2) *.95;
if THRESH(2) <= THRESH(1),
  THRESH(1) = THRESH(2) - eps;
end;
fprintf('Threshold is now [%f,%f]\n',THRESH);
% find edges in the image
EDGEI = edge(rgb2gray(RGBI0),'canny',THRESH);

% draw the edges on the image
ncolors = size(RGBI0,3);
RGBI = reshape(RGBI0,[NROWS*NCOLS],ncolors);
RGBI(EDGEI,1) = 1; RGBI(EDGEI,2) = .75; RGBI(EDGEI,3) = 1;
RGBI = reshape(RGBI,[NROWS,NCOLS,ncolors]);
set(IMAGE_H,'CData',RGBI);

% find the closest edge to each pixel
[DISTEDGE,NNEDGEI] = bwdist(EDGEI);

function decreasethresh1

global IMAGE_H NROWS NCOLS;
global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH;

THRESH(1) = THRESH(1) *.95;
if THRESH(1) == 0,
  THRESH(1) = 10^(-5);
end;
fprintf('Threshold is now [%f,%f]\n',THRESH);
% find edges in the image
EDGEI = edge(rgb2gray(RGBI0),'canny',THRESH);

% draw the edges on the image
ncolors = size(RGBI0,3);
RGBI = reshape(RGBI0,[NROWS*NCOLS],ncolors);
RGBI(EDGEI,1) = 1; RGBI(EDGEI,2) = .75; RGBI(EDGEI,3) = 1;
RGBI = reshape(RGBI,[NROWS,NCOLS,ncolors]);
set(IMAGE_H,'CData',RGBI);

% find the closest edge to each pixel
[DISTEDGE,NNEDGEI] = bwdist(EDGEI);

function increasethresh2

global IMAGE_H NROWS NCOLS;
global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH;

THRESH(2) = THRESH(2) /.95;

% find edges in the image
EDGEIt = edge(rgb2gray(RGBI0),'canny',THRESH);

if ~any(EDGEI(:)),
  THRESH(2) = THRESH(2) * .95;
else,
  EDGEI = EDGEIt;
  fprintf('Threshold is now %f %f\n',THRESH);

  % draw the edges on the image
  ncolors = size(RGBI0,3);
  RGBI = reshape(RGBI0,[NROWS*NCOLS],ncolors);
  RGBI(EDGEI,1) = 1; RGBI(EDGEI,2) = .75; RGBI(EDGEI,3) = 1;
  RGBI = reshape(RGBI,[NROWS,NCOLS,ncolors]);
  set(IMAGE_H,'CData',RGBI);

  % find the closest edge to each pixel
  [DISTEDGE,NNEDGEI] = bwdist(EDGEI);
end;

function increasethresh1

global IMAGE_H NROWS NCOLS;
global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH;

THRESH(1) = THRESH(1) /.95;
if THRESH(1) >= THRESH(2),
  THRESH(2) = THRESH(1) + eps;
end;
fprintf('Threshold is now [%f,%f]\n',THRESH);
% find edges in the image
EDGEI = edge(rgb2gray(RGBI0),'canny',THRESH);

% draw the edges on the image
ncolors = size(RGBI0,3);
RGBI = reshape(RGBI0,[NROWS*NCOLS],ncolors);
RGBI(EDGEI,1) = 1; RGBI(EDGEI,2) = .75; RGBI(EDGEI,3) = 1;
RGBI = reshape(RGBI,[NROWS,NCOLS,ncolors]);
set(IMAGE_H,'CData',RGBI);

% find the closest edge to each pixel
[DISTEDGE,NNEDGEI] = bwdist(EDGEI);

function closecurve

global BOUNDX BOUNDY CURVEN LASTLINE_H;

if isempty(BOUNDX{CURVEN}), return; end;
%BOUNDX{CURVEN} = [BOUNDX{CURVEN};BOUNDX{CURVEN}(1)];
%BOUNDY{CURVEN} = [BOUNDY{CURVEN};BOUNDY{CURVEN}(1)];
set(LASTLINE_H,'XData',BOUNDX{CURVEN},'YData',BOUNDY{CURVEN}),
set(LASTLINE_H,'color',[.75,.75,1],'linewidth',4);
CURVEN = CURVEN + 1;
LASTLINE_H = [];
BOUNDX{CURVEN} = [];
BOUNDY{CURVEN} = [];

function finish

global IMAGE_H FIG_H;
set(FIG_H,'WindowButtonMotionFcn','');
set(IMAGE_H,'ButtonDownFcn','');
tmp = IMAGE_H;
clear global FIG_H AXIS_H IMAGE_H NROWS NCOLS;
clear global CURVEN LASTLINE_H CURVEPIECESTART;
clear global EDGEI RGBI RGBI0 DISTEDGEI NNEDGEI THRESH ISEDGESHOWN;

set(tmp,'UserData','Done');

function [imageHandle,x,y] = over_image(axesHandle)

% Return the index of which image we are over, and return a 0 if we
% aren't above an image.

imagefound = findobj(axesHandle, 'type', 'image');
if isempty(imagefound)
   imageHandle=0; x=0; y=0;
   return
end
% Make sure that the Image's Button Down & Up functions will queue
set(imagefound, 'Interruptible', 'off', 'BusyAction', 'Queue');
axHandle = get(imagefound, 'Parent');
axPosition = get(axHandle, 'Position');
axCurPt = get(axHandle, 'CurrentPoint');

% See if we are above the desired axes
imageHandle = 0;
XLim = get(axHandle, 'XLim');
YLim = get(axHandle, 'YLim');
pt = axCurPt;
x = pt(1,1); y = pt(1,2);
if x>=XLim(1) & x<=XLim(2) & y>=YLim(1) & y<=YLim(2)
  imageHandle = imagefound;
else,
  x = 0; y = 0;
end;

function toggleedgeimage

global RGBI;
global RGBI0;
global ISEDGESHOWN;
global IMAGE_H;
if ISEDGESHOWN == 1,
  set(IMAGE_H,'CData',RGBI0);
  ISEDGESHOWN = 0;
else,
  set(IMAGE_H,'CData',RGBI);
  ISEDGESHOWN = 1;
end;

author: kristinbranson

repo: JAABA


text:
function pass = test_eigsGeneralizedSys(pref)
% Test generalized eigenvalueproblems for LINOP
if ( nargin == 0 )
    pref = cheboppref();
end

dom = [-1, 1];
diffOp = operatorBlock.diff(dom);
I = operatorBlock.eye(dom);
Z = operatorBlock.zeros(dom);
ev = functionalBlock.eval(dom);
z = functionalBlock.zero(dom);

A = [diffOp^2 diffOp ; diffOp diffOp^2];
B = [I diffOp ; diffOp I];
A = linop(A);
A = addbc(A, [ev(-1) z]);
A = addbc(A, [ev(1) z]);
A = addbc(A, [z ev(-1)]);
A = addbc(A, [z ev(1)]);
B = linop(B);
e_true = -1 + 1i*pi*[-1 1 -1 1 -2 2 -2 2 -3 3 -3 3].';

% CHEBCOLLOC1
pref.discretization = @chebcolloc1;
[V, D] = eigs(A, B, 12, 0, pref);
e = diag(D);
AV = A*V;
BV = B*V;
err(1,1) = norm(e - e_true);
err(1,2) = norm(AV - BV*D);

% CHEBCOLLOC2
pref.discretization = @chebcolloc2;
[V, D] = eigs(A, B, 12, 0, pref);
e = diag(D);
AV = A*V;
BV = B*V;
err(1,3) = norm(e - e_true);
err(1,4) = norm(AV - BV*D);

% ULTRAS
pref.discretization = @ultraS;
[V, D] = eigs(A, B, 12, 0, pref);
e = diag(D);
AV = A*V;
BV = B*V;
err(1,5) = norm(e - e_true);
err(1,6) = norm(AV - BV*D);

%%
tolVals = 6e-9;
tolFuns = 4e-7;

tol = repmat([tolVals, tolFuns], 1, 3);

pass = err < 10*tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function result = isintersect(SHAPE, LINE, m_obs_array, b_obs_array, vert_num)
% isintersect(SHAPE, LINE)
% this function check whether we intersect the shape or not
% SHAPE could be polygon or line
% LINE is always line
%

result = 0;
m = (LINE(2, 2) - LINE(1, 2)) /  (LINE(2, 1) - LINE(1, 1));
b =  LINE(1, 2) - m * LINE(1, 1);
radius = 0.01;

% find min and max of a line
if LINE(1, 1) >= LINE(2, 1)
    x_edge_max = LINE(1, 1);
    x_edge_min = LINE(2, 1);
else
    x_edge_max = LINE(2, 1);
    x_edge_min = LINE(1, 1);
end

if LINE(1, 2) >= LINE(2, 2)
    y_edge_max = LINE(1, 2);
    y_edge_min = LINE(2, 2);
else
    y_edge_max = LINE(2, 2);
    y_edge_min = LINE(1, 2);
end

for k = 1:vert_num
    % y = m * x + b

    m_obs = m_obs_array(k);
    b_obs = b_obs_array(k);

    % consider this lines ???

%     if (m_obs + radius > m) && (m_obs - radius < m) && (b_obs + radius > b) && (b_obs - radius < b)
%        result = 1;
%        return;
%     end

    % I couldn't find better solution to deal this min max detection =)
    % this is the fastest one
    if SHAPE(k, 1) >= SHAPE(k+1, 1)
        x_max = SHAPE(k, 1);
        x_min = SHAPE(k+1, 1);
    else
        x_max = SHAPE(k+1, 1);
        x_min = SHAPE(k, 1);
    end

    if SHAPE(k, 2) >= SHAPE(k+1, 2)
        y_max = SHAPE(k, 2);
        y_min = SHAPE(k+1, 2);
    else
        y_max = SHAPE(k+1, 2);
        y_min = SHAPE(k, 2);
    end

    x_intersection = (b - b_obs)/(m_obs - m);
    y_intersection = m_obs * x_intersection + b_obs;

    if (x_intersection >= (x_min-radius) && x_intersection <= (x_max+radius)) ...
            && (x_intersection >= (x_edge_min-radius) && x_intersection <= (x_edge_max+radius)) ...
            && (y_intersection >= (y_edge_min-radius) && y_intersection <= (y_edge_max+radius)) ...
            && (y_intersection >= (y_min-radius) && y_intersection <= (y_max+radius))
        result = 1;
        return;
    end
end

author: olzhas


text:
function s = mean(f, dim)
%MEAN   Average or mean value of a SEPARABLEAPPROX.
%   MEAN(F) takes the mean in the y-direction (default), i.e.,
%          MEAN(F) = 1/(ymax-ymin) sum(F).
%
%   MEAN(F, DIM) takes the mean along the direction DIM. If DIM = 1 it is the
%   y-direction and if DIM = 2 then it is the x-direction.
%
% See also MEAN2, STD2.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( f ) )
    s = chebfun;
    return
end

if ( nargin == 1)
    % Default to the y-direction:
    dim = 1;
end
dom = f.domain;

s = sum( f, dim );
if ( dim == 1 )
    s = s / diff( dom(3:4) ); % Mean in the y direction (default)
elseif ( dim == 2 )
    s = s / diff( dom(1:2) ); % Mean in the x direction
else
    error('CHEBFUN:SEPARABLEAPPROX:mean:dim', 'Mean not in x or y direction.')
end

end


text:
% Version 1.000
%
% Code provided by Ruslan Salakhutdinov and Geoff Hinton
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

% This program fine-tunes an autoencoder with backpropagation.
% Weights of the autoencoder are going to be saved in mnist_weights.mat
% and trainig and test reconstruction errors in mnist_error.mat
% You can also set maxepoch, default value is 200 as in our paper.

maxepoch=200;
fprintf(1,'\nFine-tuning deep autoencoder by minimizing cross entropy error. \n');
fprintf(1,'60 batches of 1000 cases each. \n');

load mnistvh
load mnisthp
load mnisthp2
load mnistpo

makebatches;
[numcases numdims numbatches]=size(batchdata);
N=numcases;

%%%% PREINITIALIZE WEIGHTS OF THE AUTOENCODER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w1=[vishid; hidrecbiases];
w2=[hidpen; penrecbiases];
w3=[hidpen2; penrecbiases2];
w4=[hidtop; toprecbiases];
w5=[hidtop'; topgenbiases];
w6=[hidpen2'; hidgenbiases2];
w7=[hidpen'; hidgenbiases];
w8=[vishid'; visbiases];

%%%%%%%%%% END OF PREINITIALIZATIO OF WEIGHTS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l1=size(w1,1)-1;
l2=size(w2,1)-1;
l3=size(w3,1)-1;
l4=size(w4,1)-1;
l5=size(w5,1)-1;
l6=size(w6,1)-1;
l7=size(w7,1)-1;
l8=size(w8,1)-1;
l9=l1;
test_err=[];
train_err=[];

for epoch = 1:maxepoch

%%%%%%%%%%%%%%%%%%%% COMPUTE TRAINING RECONSTRUCTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
err=0;
[numcases numdims numbatches]=size(batchdata);
N=numcases;
 for batch = 1:numbatches
  data = [batchdata(:,:,batch)];
  data = [data ones(N,1)];
  w1probs = 1./(1 + exp(-data*w1)); w1probs = [w1probs  ones(N,1)];
  w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
  w3probs = 1./(1 + exp(-w2probs*w3)); w3probs = [w3probs  ones(N,1)];
  w4probs = w3probs*w4; w4probs = [w4probs  ones(N,1)];
  w5probs = 1./(1 + exp(-w4probs*w5)); w5probs = [w5probs  ones(N,1)];
  w6probs = 1./(1 + exp(-w5probs*w6)); w6probs = [w6probs  ones(N,1)];
  w7probs = 1./(1 + exp(-w6probs*w7)); w7probs = [w7probs  ones(N,1)];
  dataout = 1./(1 + exp(-w7probs*w8));
  err= err +  1/N*sum(sum( (data(:,1:end-1)-dataout).^2 ));
  end
 train_err(epoch)=err/numbatches;

%%%%%%%%%%%%%% END OF COMPUTING TRAINING RECONSTRUCTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% DISPLAY FIGURE TOP ROW REAL DATA BOTTOM ROW RECONSTRUCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%
fprintf(1,'Displaying in figure 1: Top row - real data, Bottom row -- reconstructions \n');
output=[];
 for ii=1:15
  output = [output data(ii,1:end-1)' dataout(ii,:)'];
 end
   if epoch==1
   close all
   figure('Position',[100,600,1000,200]);
   else
   figure(1)
   end
   mnistdisp(output);
   drawnow;

%%%%%%%%%%%%%%%%%%%% COMPUTE TEST RECONSTRUCTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[testnumcases testnumdims testnumbatches]=size(testbatchdata);
N=testnumcases;
err=0;
for batch = 1:testnumbatches
  data = [testbatchdata(:,:,batch)];
  data = [data ones(N,1)];
  w1probs = 1./(1 + exp(-data*w1)); w1probs = [w1probs  ones(N,1)];
  w2probs = 1./(1 + exp(-w1probs*w2)); w2probs = [w2probs ones(N,1)];
  w3probs = 1./(1 + exp(-w2probs*w3)); w3probs = [w3probs  ones(N,1)];
  w4probs = w3probs*w4; w4probs = [w4probs  ones(N,1)];
  w5probs = 1./(1 + exp(-w4probs*w5)); w5probs = [w5probs  ones(N,1)];
  w6probs = 1./(1 + exp(-w5probs*w6)); w6probs = [w6probs  ones(N,1)];
  w7probs = 1./(1 + exp(-w6probs*w7)); w7probs = [w7probs  ones(N,1)];
  dataout = 1./(1 + exp(-w7probs*w8));
  err = err +  1/N*sum(sum( (data(:,1:end-1)-dataout).^2 ));
  end
 test_err(epoch)=err/testnumbatches;
 fprintf(1,'Before epoch %d Train squared error: %6.3f Test squared error: %6.3f \t \t \n',epoch,train_err(epoch),test_err(epoch));

%%%%%%%%%%%%%% END OF COMPUTING TEST RECONSTRUCTION ERROR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 tt=0;
 for batch = 1:numbatches/10
 fprintf(1,'epoch %d batch %d\r',epoch,batch);

%%%%%%%%%%% COMBINE 10 MINIBATCHES INTO 1 LARGER MINIBATCH %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 tt=tt+1;
 data=[];
 for kk=1:10
  data=[data
        batchdata(:,:,(tt-1)*10+kk)];
 end

%%%%%%%%%%%%%%% PERFORM CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  max_iter=3;
  VV = [w1(:)' w2(:)' w3(:)' w4(:)' w5(:)' w6(:)' w7(:)' w8(:)']';
  Dim = [l1; l2; l3; l4; l5; l6; l7; l8; l9];

  [X, fX] = minimize(VV,'CG_MNIST',max_iter,Dim,data);

  w1 = reshape(X(1:(l1+1)*l2),l1+1,l2);
  xxx = (l1+1)*l2;
  w2 = reshape(X(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
  xxx = xxx+(l2+1)*l3;
  w3 = reshape(X(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
  xxx = xxx+(l3+1)*l4;
  w4 = reshape(X(xxx+1:xxx+(l4+1)*l5),l4+1,l5);
  xxx = xxx+(l4+1)*l5;
  w5 = reshape(X(xxx+1:xxx+(l5+1)*l6),l5+1,l6);
  xxx = xxx+(l5+1)*l6;
  w6 = reshape(X(xxx+1:xxx+(l6+1)*l7),l6+1,l7);
  xxx = xxx+(l6+1)*l7;
  w7 = reshape(X(xxx+1:xxx+(l7+1)*l8),l7+1,l8);
  xxx = xxx+(l7+1)*l8;
  w8 = reshape(X(xxx+1:xxx+(l8+1)*l9),l8+1,l9);

%%%%%%%%%%%%%%% END OF CONJUGATE GRADIENT WITH 3 LINESEARCHES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 end

 save mnist_weights w1 w2 w3 w4 w5 w6 w7 w8
 save mnist_error test_err train_err;

end

author: qiuwch

repo: DeepLearning

save_path: github-repos/MATLAB/qiuwch-DeepLearning


text:
function fact_test ( )

%*****************************************************************************80
%
%% FACT_TEST tests the use of the MEX file FACT.F
%
%  Discussion:
%
%    The file fact.F is a FORTRAN77 function which computes the factorial.
%
%    This M file "compiles" fact.F, and then shows how it can be called.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FACT_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the MEX compiler,\n' );
  fprintf ( 1, '  which allows MATLAB to call FORTRAN77 functions.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.  The file "fact.F" should,\n' );
  fprintf ( 1, '  show up here.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compile the file "fact.F".\n' );

  mex fact.F

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get a directory listing.  A new file should show up,\n' );
  fprintf ( 1, '  containing the compiled information.\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now use FACT as though it were a MATLAB M-file function.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N  (N Factorial)' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    j = fact ( i );

    fprintf ( 1, '  %2d  %10d\n', i, j );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FACT_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

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

  return
end


text:
%KLI_demo

% S. de Waele, March 2003.

clear

a = 1;
b = rc2arset([1 -.3 .4]);
varx = 1;

nobs = 200;
x = gendata(a,b,nobs);

Lmax = 50;
ASAglob_subtr_mean = 0;
arh = sig2ar(x,Lmax);
varh= mean(x.^2);

%Fit of model to data
%And Auto-KLI
%These quantities should be equal according to eq. 4.14 on page 88

ar_set = ar2arset(arh,0:Lmax);
for t = 0:Lmax
    kli_hat(t+1) = KLIndex_hat(x,ar_set{t+1},1,varh);
    kli_auto(t+1) = KLIndex(ar_set{t+1},1,varh,ar_set{t+1},1,varh,nobs);
end

%Kullback-Leibler index (KLI)
%for increasing AR model order
[dummy, set_kli] = KLIndex(arh,1,1,a,b,varx,nobs);

%Kullback-Leibler Discrepancy (KLI)
%The KLD is equal to the KLI apart from a constant
[dummy, set_kld] = KLDiscrepancy(arh,1,1,a,b,varx,nobs);

%Akaike information criterion
aic = kli_hat + 2*(0:Lmax);

plot(0:Lmax,[set_kld' set_kli' kli_hat'],0:Lmax,kli_auto,'-.',0:Lmax,aic)
xlabel('p')
ylabel('KL')
legend('KLD','KLI','KLI_{hat}','KLI_{auto}','AIC')
legend boxoff

title('Fit and error and AIC for estimated AR models.')


text:
function [v,y,w]=v_nearnonz(x,d)
%V_NEARNONZ replace each zero element with the nearest non-zero element [V,Y,W]=v_nearnonz(X,D)
%
%  Inputs:  x         input vector, matrix or larger array
%           d         dimension to apply filter along [default 1st non-singleton]
%
% Outputs:  v         v is the same size as x but with each zero entry replaced by
%                     the nearest non-zero value along dimension d
%                     elements equidistant from two non-zero entries will be taken
%                     from the higher index
%           y         y is the same size as x and gives the index along dimension d
%                     from which the corresponding entry in v was taken
%                     If there are no non-zero entries, then the corresponding
%                     elements of y will be zero.
%           w         w is the same size as x and gives the distance (+ or -) to the
%                     nearest non-zero entry in x

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_nearnonz.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e=size(x);
p=prod(e);
if nargin<2             % if no dimension given, find the first non-singleton
    d=find(e>1,1);
    if ~numel(d)
        d=1;
    end
end
k=e(d);                 % size of active dimension
q=p/k;                  % size of remainder
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
xx=z~=0;
cx=cumsum(xx);
[i,j]=find(z);
qq=cx(xx);
pos=full(sparse(qq,j,i,k,q)); % list the positions of non-zero elements in each column
mp=ceil((pos(1:end-1,:)+pos(2:end,:))*0.5); % find the mid point between consecutive non-zero elements
[i2,j2]=find(pos(2:end,:)>0);
zz=1+cumsum(full(sparse(mp(pos(2:end,:)>0),j2,1,k,q)));
y=pos(zz+repmat((0:q-1)*k,k,1));
v=z(max(y,1)+repmat((0:q-1)*k,k,1));
w=y-repmat((1:k)',1,q);
w(y==0)=0;
if d==1
    y=reshape(y,e);
    v=reshape(v,e);
    w=reshape(w,e);
else
    y=shiftdim(reshape(y,r),length(e)+1-d);
    v=shiftdim(reshape(v,r),length(e)+1-d);
    w=shiftdim(reshape(w,r),length(e)+1-d);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [varargout] = ismatrix(varargin)

% ISMATRIX is a drop-in replacement for the same function that was
% introduced in MATLAB R2010b.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% see https://github.com/fieldtrip/fieldtrip/issues/899

alternatives = which(mfilename, '-all');
if ~iscell(alternatives)
  % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
  alternatives = {alternatives};
end

keep = true(size(alternatives));
for i=1:numel(alternatives)
  keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
  keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
end
alternatives = alternatives(keep);

if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
  % remove this directory from the path
  p = fileparts(mfilename('fullpath'));
  sw = warning('off','backtrace');
  warning('Removing "%s" from your path.\nSee http://www.fieldtriptoolbox.org/faq/should_i_add_fieldtrip_with_all_subdirectories_to_my_matlab_path/', p);
  warning(sw);
  rmpath(p);
  % call the original MATLAB function
  if exist(mfilename, 'builtin')
    [varargout{1:nargout}] = builtin(mfilename, varargin{:});
  else
    [varargout{1:nargout}] = feval(mfilename, varargin{:});
  end
  return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = ismatrix(x)

% deal with the input arguments
if nargin==1
  [x] = deal(varargin{1:1});
else
  error('incorrect number of input arguments')
end

siz = size(x);
tf = numel(siz)==2 && siz(1)>=0 && siz(2)>=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};

author: spm


text:
function jed = ymdf_to_jed_saka ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_SAKA converts a Saka YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding JED.
%

%
%  Convert the calendar date to a computational date.
%
  y_prime = y + 4794 - floor ( ( 13 - m ) / 12 );
  m_prime = mod ( m + 10, 12 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 1461 * y_prime ) / 4 );

  z = floor ( m_prime / 6 );

  j2 = ( 31 - z ) * m_prime + 5 * z;

  g = floor ( ( y_prime + 184 ) / 100 );
  g = floor ( ( 3 * g ) / 4 ) - 36;

  jed = j1 + j2 + d_prime - 1348 - g - 0.5 + f;

  return
end


text:
function ctrParams = ctrInitBatchParams(varargin)
%  Initialize the parameters structure used with ctrInitBatchTrack when
%  doing multi-subject tractography using the conTrack algorithm.
%
%  ctrParams = ctrInitBatchParams([varargin])
%
% VARIABLES:
%    ctrParams ...
%      .projectName = ctrParams.projectName = ['myConTrackProj'];
%                     This variable will be used in the name of all the
%                     files that are created. E.g., outFile =
%                     ['fg_',projectName,'_',roi1,'_',roi2,'_',timeStamp,'.
%                     pdb'];
%      .logName     = ctrParams.logName- = ['myConTrackProj'];
%                     This will be used to provide a unique name for easy
%                     ID in the log directory.
%      .baseDir     = Top-level directory containing your data.
%                     The level below baseDir should have each subjects
%                     data directory.
%                     * e.g., ctrParams.baseDir = '/biac/wandell/data/dti/';
%      .dtDir       = This should be the name of the directory containing
%                     the dt6.mat file.* Relative to the subjects directory.
%                     * e.g., ctrParams.dtDir = 'dti40trilinrt';.
%      .roiDir      = Directory containing the ROIs. * Relative to the
%                     subjects directory.
%                     * e.g., ctrParams.roiDir = 'ROIs';
%      .subs        = This is the cell array that will contain the names
%                     of all the sujbect's directories that will be
%                     processed.
%                     * e.g., ctrParams.subs = {'sub1','sub2','sub3'};
%                     NOTE: This script also supports the ability to load a
%                     list of subjects from a text file. If you wish to do
%                     this simply leave the cell empty. You will be
%                     prompted to select a text file that contains a list
%                     of subjects. Please assure that this list is a simple
%                     text file with only subject names seperated by new
%                     lines or spaces. * This can also be a path to a text
%                     file with subject IDs.
%      .roi1/roi2   = These two cell arrays should contain the names of
%                     each ROI to be used in tracking. The script will
%                     track from ROI1{1} to ROI2{1} and ROI1{2} to ROI2{2}
%                     etc... Use .mat files, but DO NOT include file
%                     extensions. In case that you wish to track from
%                     multiple rois (ROI1) to the same roi (ROI2) you can
%                     just place the name of one roi in ROI2 and each roi
%                     in ROI1 will be tracked to the single roi in ROI2.
%                     * e.g., ctrParams.roi1 = {'Roi1a','Roi1b'};
%                             ctrParams.roi2 = {'Roi2a','Roi2b'};
%      .nSamples     = [50000];    % Enter the number of pathway samples to generate.
%      .maxNodes     = [240];      % Enter the max length of the samples.
%      .minNodes     = [10];       % Enter the minimum length of the samples.
%      .stepSize     = [1];        % Enter the step size.
%      .pddpdfFlag   = [0];        % 0 = Only compute if file does not already exist. 1= Always recompute.
%      .wmFlag       = [0];        % 0 = Only compute if file does not already exist. 1= Always recompute.
%      .roi1SeedFlag = ['true'];   % We usually want to equally seed both ROIs, so both flags = 'true'.
%      .roi2SeedFlag = ['true'];   % For speed you can choose not to seed the second ROI
%      .multiThread  = [0];        % Spawn many tracking jobs at once.
%      .executeSh    = [0];        % Execute the script on cuurent host immediately using an xterm.
%
% EXAMPLE USUAGE:
%       ctrParams          = ctrInitBatchParams;
%       ctrParams.baseDir  = '/Directory/Containing/Data';
%       ctrParams.dtDir    = 'dti40trilin';
%       ctrParams.roiDir   = 'ROIs'; % Relative to the subject dir. Contains the ROIs
%       ctrParams.subs     = {'subDir1','subDir2'};
%       ctrParams.ROI1     = {'Roi1a','Roi1b'}; % assumes '.mat'
%       ctrParams.ROI2     = {'Roi2a','Roi2b'}; % assumes '.mat'
%       ctrParams.nSamples = 10000;
%       [cmd info] = ctrInitBatch(ctrParams);
%
% WEB RESOURCES:
%       http://white.stanford.edu/newlm/index.php/ConTrack
%       mrvBrowseSVN('ctrInitBatchParams');
%
% See Also:
%       ctrInitBatchTrack.m , ctrInitBatchScore.m
%
%
% (C) Stanford Vista, 8/2011 [lmp]
%
%
%% Initialize ctrParams struct

ctrParams = struct;

%% I. Set Naming and Directory Structure

ctrParams.projectName = 'myConTrackProj';
ctrParams.logName     = 'myConTrackLog';
ctrParams.baseDir     = '';
ctrParams.dtDir       = '';
ctrParams.roiDir      = '';

%% II. Set Subjects and ROIs

ctrParams.subs = {};
ctrParams.roi1 = {};
ctrParams.roi2 = {};

%% III. Set Algorithm Parameters

ctrParams.nSamples     = 50000;
ctrParams.maxNodes     = 240;
ctrParams.minNodes     = 10;
ctrParams.stepSize     = 1;
ctrParams.pddpdfFlag   = 0;
ctrParams.wmFlag       = 0;
ctrParams.roi1SeedFlag = 'true';
ctrParams.roi2SeedFlag = 'true';
ctrParams.multiThread  = 0;
ctrParams.executeSh    = 0;

%% Varargin

ctrParams = mrVarargin(ctrParams,varargin);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function s = yjf_to_s_julian ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_S_JULIAN writes a Julian YJF date into a string.
%
%  Format:
%
%    AD YYYY/JJJ.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, string S, the representation of the date.
%
  if ( y < 0 )
    s = sprintf ( 'BC %d/%06.2f', -y, j + f );
  else
    s = sprintf ( 'AD %d/%06.2f', y, j + f );
  end

  return
end


text:
function opt= defopt_scalp_power(varargin)
%DEFOPT_SCALP_POWER - Properties for plot_scalp* functions, optimized for
%visualizing spectral power

props= {'Shading'                 'flat'
        'Extrapolation'           1
        'ExtrapolateToMean'       0
        'Resolution'              61
        'CLim'                    'range'
        'IvalColor'               [0.8 0.8 0.8; 0.6 0.6 0.6]
        'LineWidth'               2
        'ChannelLineStyleOrder'   {'thick' 'thin'}
        'Colormap'                cmap_hsvFade(51, [4/6 0], 1, 1)
        'Contour'                 5
        'ContourPolicy'           'strict'
        'ContourLineprop'         {'LineWidth' 0.3}
        'MarkMarkerProperties'    {'MarkerSize' 6, 'LineWidth' 1}
        'GlobalCLim'              1
        'LegendPos'               'NorthEast'};

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);

author: bbci


text:
function [ifeat,testStat,testStat_rand,featureClassifier] = TS_TopFeatures(whatData,whatTestStat,cfnParams,varargin)
% TS_TopFeatures  Top individual features for discriminating labeled time series
%
% This function compares each feature in an hctsa dataset individually for its
% ability to separate the labeled classes of time series according to a given
% test statistic.
%
% Can also compare this performance to a set of randomized null features to
% evaluate the statistical significance of the result (pooled permutation testing).
%
% Note that it can be more interpretable to run this function on filtered but
% un-normalized data (e.g., from TS_Normalize('none')), so that feature values
% correspond to raw feature outputs.
%
%---INPUTS:
%
% whatData, the hctsa data to use (input to TS_LoadData, default: 'raw')
% whatTestStat, the test statistic to quantify the goodness of each feature
%               (e.g., 'tstat','ustat'; or set 'classification' to use classifier
%                   described in cfnParams)
% cfnParams, the classification settings if using 'classification'-based selection
%
%---OPTIONAL extra inputs:
%
% 'whatPlots', can specify what output plots to produce (cell of strings), e.g.,
%               specify {'histogram','distributions','cluster','datamatrix'} to
%               produce all four possible output plots (this is the default).
% 'numTopFeatures' [40], specifies the number of top features to analyze, both in
%                   terms of the list of outputs, the histogram plots, and the
%                   cluster plot.
% 'numFeaturesDistr' [16], sets a custom number of distributions to display (can
%                   set this lower to avoid producing large numbers of figures).
% 'numNulls' [0], the number of shuffled nulls to generate (e.g., 10 shuffles pools
%               shuffles for all M features, for a total of 10*M elements in the
%               null distribution)
%
%---EXAMPLE USAGE:
%
% TS_TopFeatures('norm','ustat',struct(),'whatPlots',{'histogram','distributions',...
%           'cluster','datamatrix'},'numTopFeatures',40,'numFeaturesDistr',10,'numNulls',10);
%
% TS_TopFeatures('norm','classification',cfnParams,'whatPlots',{'histogram','distributions',...
%           'cluster','datamatrix'},'numTopFeatures',40,'numFeaturesDistr',10,'numNulls',10);
%
%---OUTPUTS:
%
% ifeat, the ordering of operations by their performance
% testStat, the test statistic (whatTestStat) for each operation
% testStat_rand, test statistics making up the null distributions
% featureClassifier, the individual-feature classifier (if using single-feature
%               classification)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%% Check inputs and set defaults
%-------------------------------------------------------------------------------
if nargin < 1 || isempty(whatData)
    whatData = 'norm';
end
% Set other defaults later when we know what we're dealing with

% Use an inputParser to control additional options as parameters:
inputP = inputParser;
% whatPlots
default_whatPlots = {'histogram','distributions','dendrogram','cluster'}; % 'datamatrix'
check_whatPlots = @(x) iscell(x) || ischar(x);
addParameter(inputP,'whatPlots',default_whatPlots,check_whatPlots);
% numTopFeatures
default_numTopFeatures = 40;
addParameter(inputP,'numTopFeatures',default_numTopFeatures,@isnumeric);
% numFeaturesDistr
default_numFeaturesDistr = 16;
addParameter(inputP,'numFeaturesDistr',default_numFeaturesDistr,@isnumeric);
% numNulls
default_numNulls = 0; % by default, don't compute an empirical null distribution
                      % by randomizing class labels
addParameter(inputP,'numNulls',default_numNulls,@isnumeric);

parse(inputP,varargin{:});

whatPlots = inputP.Results.whatPlots;
if ischar(whatPlots)
    whatPlots = {whatPlots};
end
numTopFeatures = inputP.Results.numTopFeatures;
numFeaturesDistr = inputP.Results.numFeaturesDistr;
numNulls = inputP.Results.numNulls;
clear('inputP');

%-------------------------------------------------------------------------------
%% Load the data
%-------------------------------------------------------------------------------
[TS_DataMat,TimeSeries,Operations,whatDataFile] = TS_LoadData(whatData);

%-------------------------------------------------------------------------------
%% Checks
%-------------------------------------------------------------------------------
% Does the grouping information exist?:
% Assign group labels (removing unlabeled data):
[TS_DataMat,TimeSeries] = FilterLabeledTimeSeries(TS_DataMat,TimeSeries);
[groupLabels,classLabels,groupLabelsInteger,numClasses] = TS_ExtractGroupLabels(TimeSeries);
% Give basic info about the represented classes:
TellMeAboutLabeling(TimeSeries);

%-------------------------------------------------------------------------------
% Set defaults for the test statistic:
if nargin < 2 || isempty(whatTestStat)
    if numClasses == 2
        % ustat or ttest are way faster than proper prediction models
        whatTestStat = 'ustat';
    else
        whatTestStat = 'classification';
    end
    fprintf(1,'Using ''%s'' test statistic by default.\n', whatTestStat);
end
% Set cfnParams default (if using 'classification') later after loading the time series
if strcmp(whatTestStat,'classification')
    if nargin < 3 || isempty(fieldnames(cfnParams))
        cfnParams = GiveMeDefaultClassificationParams(TimeSeries);
        cfnParams.whatClassifier = 'fast-linear';
        cfnParams = UpdateClassifierText(cfnParams);
    end
else
    cfnParams = struct();
end

%-------------------------------------------------------------------------------
% Filter down to reduced features if specified/required:
[TS_DataMat,Operations] = FilterFeatures(TS_DataMat,Operations,cfnParams);
numFeatures = size(TS_DataMat,2);
numTopFeatures = min(numTopFeatures,numFeatures);

%-------------------------------------------------------------------------------
%% Define the train/test classification rate function, fn_testStat
%-------------------------------------------------------------------------------
% chanceLevel -- what you'd expect by chance

% Check that simple stats are being applied just for pairs:
if ismember(whatTestStat,{'ustat','ranksum','ustatExact','ustatP','ranksumExact',...
                            'ttest','ttestP','tstat','svmBeta'})
    if numClasses~=2
        error(['There are %u labeled classes.\nStatistics like ''%s'' are only valid',...
                    ' for problems with 2 labeled classes'],numClasses,whatTestStat);
    end
end

% Set up the test statistic computation and relevant description text:
switch whatTestStat
    case 'classification'
        % Set up the loss function for a classifier-based metric:
        fn_testStat = GiveMeFunctionHandle(cfnParams);
        chanceLevel = 100/numClasses; % (could be a bad assumption: accuracy for equiprobable groups...)
        testStatText = sprintf('%s %s',cfnParams.classifierText,cfnParams.whatLoss);
        statUnit = cfnParams.whatLossUnits;
        whatIsBetter = 'high';
    case 'svmBeta'
        % Feature weight from a linear (in-sample) SVM
        fn_testStat = [];
        chanceLevel = NaN;
        testStatText = '(multivariate) Linear SVM (in-sample) feature weight';
        statUnit = '';
        whatIsBetter = 'abs';
    case {'ustatP','ranksumP'}
        % Approximate p-value from ranksum test
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                fn_uStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),false,true,false);
        chanceLevel = NaN;
        testStatText = 'Mann-Whitney approx p-value';
        statUnit = ' (-log10(p))';
        whatIsBetter = 'high';
    case {'ustatPsign','ranksumPsign'}
        % Approximate p-value from ranksum test, signed by the direction of difference
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                fn_uStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),false,true,true);
        chanceLevel = NaN;
        testStatText = 'Mann-Whitney approx p-value (signed)';
        statUnit = ' signed(-log10(p))';
        whatIsBetter = 'abs';
    case {'ustatExactP','ranksumExactP'}
        % Exact p-value from ranksum test
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                    fn_uStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),true,true,false);
        chanceLevel = NaN;
        testStatText = 'Mann-Whitney exact p-value';
        statUnit = ' (-log10(p))';
        whatIsBetter = 'high';
    case {'ustat','ranksum'}
        % Approximate test statistic from ranksum test
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                fn_uStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),false,false,false);
        chanceLevel = NaN; % (check if you can compute the chance level stat?)
        testStatText = 'Mann-Whitney approx test statistic';
        statUnit = '';
        whatIsBetter = 'high';
    case {'ttest','tstat'}
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                        fn_tStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),false);
        chanceLevel = 0; % chance-level t statistic is zero
        testStatText = 'Welch''s t-stat';
        statUnit = '';
        whatIsBetter = 'abs';
    case {'ttestP','tstatP'}
        fn_testStat = @(XTrain,yTrain,Xtest,yTest) ...
                        fn_tStat(XTrain(yTrain==classLabels{1}),XTrain(yTrain==classLabels{2}),true);
        chanceLevel = NaN;
        testStatText = 'Welch''s t-test p-value';
        statUnit = ' (-log10(p))';
        whatIsBetter = 'high';
    otherwise
        error('Unknown test statistic: ''%s''',whatTestStat);
end

%-------------------------------------------------------------------------------
%% Loop over all features
%-------------------------------------------------------------------------------
% Use the same data for training and testing:
fprintf(1,['Computing the performance of %u individual features to differentiate',...
                ' %u classes using a %s...\n'],...
                        height(Operations),numClasses,testStatText);
timer = tic;
testStat = giveMeStats(TS_DataMat,TimeSeries.Group,true);
fprintf(1,' Done in %s.\n',BF_TheTime(toc(timer)));

if all(isnan(testStat))
    error('Error computing test statistics for %s (maybe there are missing data?)',...
                testStatText);
end

%-------------------------------------------------------------------------------
% Give mean and that expected from random classifier (there may be a little overfitting)
if ~isnan(chanceLevel)
    fprintf(1,['Mean %s across %u features = %4.2f%s\n' ...
            '(Random guessing for %u equiprobable classes = %4.2f%s)\n'], ...
        testStatText,numFeatures,nanmean(testStat),statUnit,numClasses,chanceLevel,statUnit);
else
    fprintf(1,'Mean %s across %u features = %4.2f%s\n',...
        testStatText,numFeatures,nanmean(testStat),statUnit);
end

%-------------------------------------------------------------------------------
%% Display information about the top features (numTopFeatures)
%-------------------------------------------------------------------------------
switch whatIsBetter
case 'high'
    % Bigger numbers indicate better features:
    [testStat_sort,ifeat] = sort(testStat,'descend');
case 'abs'
    % Bigger magnitudes indicate better features:
    [~,ifeat] = sort(abs(testStat),'descend');
    testStat_sort = testStat(ifeat);
otherwise
    error('Unknown ''whatIsBetter'' option: %s',whatIsBetter)
end

isNaN = isnan(testStat_sort);
testStat_sort = testStat_sort(~isNaN);
ifeat = ifeat(~isNaN);

% List the top features:
for i = 1:numTopFeatures
    ind = ifeat(i);
    fprintf(1,'[%u] %s (%s) -- %4.2f%s\n',Operations.ID(ind),...
            Operations.Name{ind},Operations.Keywords{ind},...
            testStat_sort(i),statUnit);
end

%-------------------------------------------------------------------------------
%% Define the indices of features we will go on to visualize
% (and construct text labels for them)
%-------------------------------------------------------------------------------
featInd = ifeat(1:numTopFeatures);
makeLabelMinimal = @(x) sprintf('[%u] %s (%4.2f%s)',Operations.ID(x),Operations.Name{x},...
                        testStat(x),statUnit);
makeLabelFull = @(x) sprintf('[%u] %s (%s) (%4.2f%s)',Operations.ID(x),Operations.Name{x},...
                    Operations.Keywords{x}, testStat(x),statUnit);
topFeatureLabelsMinimal = arrayfun(@(x) makeLabelMinimal(x),featInd,'UniformOutput',false);
topFeatureLabelsFull = arrayfun(@(x) makeLabelFull(x),featInd,'UniformOutput',false);

%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
%% Plot outputs
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
%% Histogram of distribution of test statistics for labeled and null data
%-------------------------------------------------------------------------------
if ismember('histogram',whatPlots)
    % Plot the distribution of test statistics across all features:

    %-------------------------------------------------------------------------------
    % Compute null distribution
    %-------------------------------------------------------------------------------
    testStat_rand = zeros(numFeatures,numNulls);
    if numNulls > 0
        fprintf(1,'Now for %u nulls... ',numNulls);
        tic
        for j = 1:numNulls
            if j < numNulls
                fprintf(1,'%u,',j);
            else
                fprintf(1,'%u.',j);
            end
            shuffledLabels = TimeSeries.Group(randperm(height(TimeSeries)));
            testStat_rand(:,j) = giveMeStats(TS_DataMat,shuffledLabels,false);
        end
        fprintf(1,'\n%u %s statistics computed in %s.\n',numFeatures*numNulls,...
                                        testStatText,BF_TheTime(toc(timer)));

        % Pool nulls to estimate p-values
        pooledNulls = testStat_rand(:);
        pVals = arrayfun(@(x) mean(testStat(x) < pooledNulls),1:length(testStat));
        % FDR-corrected q-values:
        FDR_qVals = mafdr(pVals,'BHFDR','true');
        fprintf(1,'Estimating FDR-corrected p-values across all features by pooling across %u nulls.\n',numNulls);
        fprintf(1,'(Given strong dependences across %u features, will produce conservative p-values)\n',numFeatures);

        % Give summary:
        sigThreshold = 0.05;
        if any(FDR_qVals < sigThreshold)
            fprintf(1,['%u/%u features show better performance using %s than the null distribution\n' ...
                        'at the magical 0.05 threshold (FDR corrected).\n'],...
                            sum(FDR_qVals < 0.05),length(FDR_qVals),testStatText);
        else
            fprintf(1,['Tough day at the office, hey? No features show statistically better performance than ' ...
                    'the null distribution at a FDR of 0.05.\nDon''t you go p-hacking now, will you?\n']);
        end
    end

    %---------------------------------------------------------------------------
    % Plot histogram
    %---------------------------------------------------------------------------
    f = figure('color','w'); hold('on')
    f.Position(3:4) = [559,278];
    colors = BF_GetColorMap('spectral',5,1);
    if numNulls == 0
        % Just plot the real distribution of test statistics across all features
        h_real = histogram(testStat,'Normalization','probability',...
                    'BinMethod','auto','FaceColor',colors{4},'EdgeColor','k');
        maxH = max(h_real.Values);
    else
        % Plot both real distribution and null distribution:
        numBins = 20;
        allTestStat = [testStat(:); testStat_rand(:)];
        minMax = [min(allTestStat),max(allTestStat)];
        histEdges = linspace(minMax(1),minMax(2),numBins+1);
        h_null = histogram(testStat_rand(:),histEdges,'Normalization','probability','FaceColor',colors{1});
        h_real = histogram(testStat,histEdges,'Normalization','probability',...
                                'FaceColor',colors{5},'EdgeColor','k');
        maxH = max([max(h_real.Values),max(h_null.Values)]);
        l_meannull = plot(nanmean(testStat_rand(:))*ones(2,1),[0,maxH],'--','color',colors{1},'LineWidth',2);
    end

    % Add chance line:
    if ~isnan(chanceLevel)
        l_chance = plot(chanceLevel*ones(2,1),[0,maxH],'--','color',colors{1},'LineWidth',2);
    end

    % Add mean of real distribution:
    l_mean = plot(nanmean(testStat)*ones(2,1),[0,maxH],'--','color',colors{5},'LineWidth',2);

    % Labels:
    xlabel(sprintf('Individual %s %s across %u features',testStatText,statUnit,numFeatures))
    ylabel('Probability')

    % Legend:
    if numNulls > 0
        if ~isnan(chanceLevel)
            legend([h_null,h_real,l_chance,l_meannull,l_mean],'null','real','chance','null mean','mean')
        else
            legend([h_null,h_real,l_meannull,l_mean],'null','real','null mean','mean')
        end
        title(sprintf(['%u features significantly informative of groups',...
                        ' (FDR-corrected at 0.05)'],sum(FDR_qVals < 0.05)))
    else
        if ~isnan(chanceLevel)
            legend([h_real,l_chance,l_mean],{'real','chance','mean'});
        else
            legend([h_real,l_mean],{'real','mean'});
        end
    end
else
    testStat_rand = [];
end

%-------------------------------------------------------------------------------
%% Distributions across classes for top features
%-------------------------------------------------------------------------------
if ismember('distributions',whatPlots)
    subPerFig = 16; % subplots per figure

    % Space the figures out properly:
    numFigs = ceil(numFeaturesDistr/subPerFig);

    % Make data structure for TS_SingleFeature
    data = struct('TS_DataMat',TS_DataMat,'TimeSeries',TimeSeries,...
                'Operations',Operations);

    for figi = 1:numFigs
        if figi*subPerFig > length(ifeat)
            break % We've exceeded number of features
        end
        % Get the indices of features to plot
        r = ((figi-1)*subPerFig+1:figi*subPerFig);
        if figi==numFigs % filter down for last one
            r = r(r <= numFeaturesDistr);
        end
        featHere = ifeat(r); % features to plot on this figure
        % Make the figure
        f = figure('color','w');
        f.Position(3:4) = [1353, 857];
        % Loop through features
        for opi = 1:length(featHere)
            subplot(ceil(length(featHere)/4),4,opi);
            TS_SingleFeature(data,Operations.ID(featHere(opi)),true,false,...
                                                testStat(featHere(opi)),false);
        end
    end
end

%-------------------------------------------------------------------------------
%% Data matrix containing top features
%-------------------------------------------------------------------------------
if ismember('datamatrix',whatPlots)
    ixFeat = BF_ClusterReorder(TS_DataMat(:,featInd)','corr','average');
    dataLocal = struct('TS_DataMat',BF_NormalizeMatrix(TS_DataMat(:,featInd(ixFeat)),'maxmin'),...
                    'TimeSeries',TimeSeries,...
                    'Operations',Operations(featInd(ixFeat),:));
    TS_PlotDataMatrix(dataLocal,'colorGroups',true,'groupReorder',true);
end

%-------------------------------------------------------------------------------
%% Pre-computing pairwise dependencies
% (for annotated dendrogram or feature–feature 'cluster' plot)
%-------------------------------------------------------------------------------
if ismember('dendrogram',whatPlots) || ismember('cluster',whatPlots)
    % 1. Compute pairwise Spearman correlations between all pairs of top features

    % (If it exists already, use that; otherwise compute it on the fly)
    % clustStruct = TS_GetFromData(whatData,'op_clust');
    % if isempty(clustStruct) % doesn't exist
    %     clustStruct = struct();
    % end
    % if isfield(clustStruct,'Dij') && ~isempty(clustStruct.Dij)
    %     % pairwise distances already computed, stored in the HCTSA .mat file
    %     fprintf(1,'Loaded %s distances from %s\n',clustStruct.distanceMetric,whatDataFile);
    %     Dij = squareform(clustStruct.Dij);
    %     Dij = Dij(featInd,featInd);
    %     distanceMetric = clustStruct.distanceMetric;
    % else
    %     % Compute correlations on the fly
    %     Dij = BF_pdist(TS_DataMat(:,featInd)','abscorr');
    %     distanceMetric = 'abscorr';
    % end

    % spearmanCorrs = corr(TS_DataMat(:,featInd)','abscorr');
    distanceMetric = 'abscorr'; % absolute Spearman correlation distances
    Dij = BF_pdist(TS_DataMat(:,featInd)',distanceMetric);
end

%-------------------------------------------------------------------------------
%% Dendrogram
% An annotated dendrogram showing the similarity between the set of top features
% (like the feature–feature correlation matrix but without the matrix)
%-------------------------------------------------------------------------------
if ismember('dendrogram',whatPlots)
    colorThreshold = 0.25;
    doVertical = false;
    f = figure('color','w');

    if doVertical
        f.Position(3:4) = [950,750];
    else
        f.Position(3:4) = [750,950];
    end

    [distVec,links] = BF_DoLinkage(Dij,distanceMetric,'average');
    % Get the dendrogram reordering:
    ord = BF_LinkageOrdering(distVec,links);
    if doVertical
        h_dend = dendrogram(links,0,'Orientation','right','Reorder',ord,...
                    'ColorThreshold',colorThreshold);
    else
        h_dend = dendrogram(links,0,'Orientation','top','Reorder',ord,...
                    'ColorThreshold',colorThreshold);
    end

    % Thicken lines:
    set(h_dend,'LineWidth',2)

    % Fix colors:
    BF_RecolorDendrogram(h_dend);

    % Label:
    box('on')
    ax = gca();
    ax.TickLabelInterpreter = 'none';
    if doVertical
        ax.YDir = 'reverse'; % needs to be reversed to match the imagesc plot (if 'cluster' is selected)
        ax.YTick = 1:numTopFeatures;
        ax.YTickLabel = topFeatureLabelsFull(ord);
        xlabel(sprintf('%s distance',distanceMetric));
    else
        ax.XTick = 1:numTopFeatures;
        ax.XTickLabel = topFeatureLabelsFull(ord);
        ylabel(sprintf('%s distance',distanceMetric));
        ax.XTickLabelRotation = 40;
    end
end

%-------------------------------------------------------------------------------
%% Inter-dependence of top features
%-------------------------------------------------------------------------------
% Plot the pairwise feature–feature correlation matrix (using BF_ClusterDown)
if ismember('cluster',whatPlots)
    clusterThreshold = 0.25; % threshold at which split into clusters
    [~,cluster_Groupi] = BF_ClusterDown(Dij,'clusterThreshold',clusterThreshold,...
                        'whatDistance',distanceMetric,...
                        'objectLabels',topFeatureLabelsMinimal);
    title(sprintf('Dependencies between %u top features (organized into %u clusters)',...
                            numTopFeatures,length(cluster_Groupi)))
end

%-------------------------------------------------------------------------------
%% Save a classifier (for the single best feature) to file
%-------------------------------------------------------------------------------
featureClassifier = struct();
if isfield(cfnParams,'classifierFilename') && ~isempty(cfnParams.classifierFilename)
    theTopFeatureIndex = ifeat(1);
    topFeatureValues = TS_DataMat(:,theTopFeatureIndex);

    % Build in-sample classifier (no CV):
    numFoldsAbove = cfnParams.numFolds;
    cfnParams.numFolds = 0;

    % Train a classification model on the top feature:
    GiveMeCfn(XTrain,yTrain,XTest,yTest,cfnParams,beVerbose)
    [bestTestStat,bestMdl,whatTestStat] = GiveMeCfn(topFeatureValues,TimeSeries.Group,...
                                                    topFeatureValues,TimeSeries.Group,...
                                                    cfnParams);

    % Prepare the featureClassifier structure:
    featureClassifier.Operation.ID = Operations.ID(theTopFeatureIndex); % Sorted list of top operations
    featureClassifier.Operation.Name = Operations.Name{theTopFeatureIndex}; % Sorted list of top operations
    featureClassifier.Mdl = bestMdl; % sorted feature models
    if numFoldsAbove > 0
        % From sorted accuracy of models
        featureClassifier.CVAccuracy = testStat_sort(1);
    else
        % CV was never done
        featureClassifier.CVAccuracy = NaN;
    end
    featureClassifier.Accuracy = bestTestStat;
    featureClassifier.whatTestStat = whatTestStat;
    featureClassifier.normalizationInfo = TS_GetFromData(whatData,'normalizationInfo');
    classes = cfnParams.classLabels;

    % Save to file
    if exist(cfnParams.classifierFilename,'file')~=0
        save(cfnParams.classifierFilename,'featureClassifier','classes','-append');
        fprintf(1,'Appended individual-feature classifiers to %s\n',cfnParams.classifierFilename);
    else
        save(cfnParams.classifierFilename,'featureClassifier','classes','-v7.3');
        fprintf(1,'Saved individual-feature classifiers to %s\n',cfnParams.classifierFilename);
    end
end

%-------------------------------------------------------------------------------
% Don't display unused outputs to screen:
if nargout == 0
    clear('ifeat','testStat','testStat_rand','featureClassifier');
end

%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
function [theStatistic,stats] = fn_uStat(d1,d2,doExact,doP,doSigned)
    % Return a statistic from a Mann-Whitney U-test
    if doExact
        % use the exact method (slower)
        [p,~,stats] = ranksum(d1,d2,'method','exact');
    else
        % use the approximate method (faster)
        [p,~,stats] = ranksum(d1,d2,'method','approx');
    end
    if doP
        % -log10(p) from a Mann-Whitney U-test
        theStatistic = -log10(p);
    else
        % ranksum test statistic from a Mann-Whitney U-test
        theStatistic = stats.ranksum;
    end
    if doSigned
        % Gives the statistic the sign of the z statistic
        % (tracks wether higher or lower in the class)
        theStatistic = sign(stats.zval)*theStatistic;
    end
    % theStatistic = sign(stats.ranksum)*(-log10(p)); % (signed p-value)
end
%-------------------------------------------------------------------------------
function [theStatistic,stats] = fn_tStat(d1,d2,doP)
    % Return a statistic from a 2-sample Welch's t-test
    [~,p,~,stats] = ttest2(d1,d2,'Vartype','unequal');
    if doP
        % -log10(p) from the 2-sample Welch's t-test
        theStatistic = -log10(p);
    else
        % t-statistic from the 2-sample Welch's t-test
        theStatistic = stats.tstat;
    end
end
%-------------------------------------------------------------------------------
function [testStat,Mdl] = giveMeStats(dataMatrix,groupLabels,beVerbose)
    % Return a test statistic for each feature:
    testStat = zeros(numFeatures,1);
    if nargout==2
        Mdl = cell(numFeatures,1);
    end

    if strcmp(whatTestStat,'svmBeta')
        % Fit a linear SVM (multivariate) and return feature weights:
        % Mdl = fitcsvm(dataMatrix,dataLabels,'KernelFunction','linear');
        Mdl = fitcsvm(dataMatrix,groupLabels,'KernelFunction','linear','Weights',InverseProbWeight(groupLabels));
        testStat = Mdl.Beta; % (code assumes bigger is better)
        testStat = testStat/max(abs(testStat)); % rescale by max-abs
    else
        loopTimer = tic;
        BF_ProgressBar('new')
        for k = 1:numFeatures
            try
                if nargout == 2
                    % This is slower for the fast-linear classifier (but returns a model)
                    [testStat(k),Mdl{k}] = fn_testStat(dataMatrix(:,k),groupLabels,dataMatrix(:,k),groupLabels);
                else
                    testStat(k) = fn_testStat(dataMatrix(:,k),groupLabels,dataMatrix(:,k),groupLabels);
                end
            catch
                warning('Error computing %s test statistic for feature %u.',whatTestStat,k);
            end
            BF_ProgressBar(k/numFeatures)
        end
        BF_ProgressBar('close')
    end
end
%-------------------------------------------------------------------------------

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%-------------------------------------------------------------------------------------------------------------
% This is an implementation of the TWSC algorithm for real-world image denoising
% by processing the YUV channels insead of RGB channels, it would improve
% TWSC by 0.1dB on PSNR.
% Author:  Jun Xu, csjunxu@comp.polyu.edu.hk / nankaimathxujun@gmail.com
%          The Hong Kong Polytechnic University
%
% Please refer to the following paper if you find this code helps:
%
% @article{TWSC_ECCV2018,
% 	author = {Jun Xu and Lei Zhang and David Zhang},
% 	title = {A Trilateral Weighted Sparse Coding Scheme for Real-World Image Denoising},
% 	journal = {ECCV},
%   year = {2018}
% }
% Please see the file License.txt for the license governing this code.
%-------------------------------------------------------------------------------------------------------------
clear,clc;
load '/home/csjunxu/Github/data/NTIRE2019/ValidationNoisyBlocksSrgb.mat';

method = 'TWSCrgb2yuv';
dataset = 'SIDD_2018';
% write image directory
write_MAT_dir = ['/home/csjunxu/Github/data/NTIRE2019/'];
% Parameters
Par.ps = 6;        % patch size
Par.step = 3;      % the step of two neighbor patches
Par.win = 20;      % size of window around the patch
Par.Outerloop = 4;
Par.Innerloop = 2;
Par.nlspini = 70;
Par.display = 0;
Par.delta = 0;
Par.nlspgap = 0; %10
Par.lambda1 = 0;
% Par.lambda2 = 3;
for lambda2 = [4]
    results = zeros(size(ValidationNoisyBlocksSrgb));
    Par.lambda2 =lambda2;
    save([write_MAT_dir 'TWSC' num2str(lambda2) '.mat'], 'results');
    for i = 1:size(ValidationNoisyBlocksSrgb,1)
        for j = 1:size(ValidationNoisyBlocksSrgb,2)
            Par.nlsp = Par.nlspini;  % number of non-local patches
            % iterate over bounding boxes
            Par.nim = im2double(squeeze(ValidationNoisyBlocksSrgb(i,j,:,:,:)));
            Par.I = Par.nim;
            Par.nim = rgb2ycbcr(Par.nim);
            [h,w,ch] = size(Par.nim);
            % noise estimation
            for c = 1:ch
                Par.nSig(c) = NoiseEstimation(Par.nim(:, :, c)*255, Par.ps)/255;
            end
            [IMout, Par]  =  TWSC_RW(Par);
            IMout = ycbcr2rgb(IMout);
            fprintf('%s/%s: \n', num2str(i), num2str(j));
            %% output
            results(i,j,:,:,:) = uint8(IMout*255);
        end
        fprintf('Image %d/%d done\n', i,j);
    end
end

author: csjunxu


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% /**
%  * Copyright (C) 2016 CoDyCo
%  * @author: Daniele Pucci, Gabriele Nava
%  * Permission is granted to copy, distribute, and/or modify this program
%  * under the terms of the GNU General Public License, version 2 or any
%  * later version published by the Free Software Foundation.
%  *
%  * A copy of the license can be found at
%  * http://www.robotcub.org/icub/license/gpl.txt
%  *
%  * This program is distributed in the hope that it will be useful, but
%  * WITHOUT ANY WARRANTY; without even the implied warranty of
%  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
%  * Public License for more details
%  */
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% NOTE: THIS SCRIPT IS RUN AUTOMATICALLY WHEN THE USER STARTS THE ASSOCIATED
% SIMULINK MODEL. NO NEED TO RUN THIS SCRIPT EVERY TIME.
clearvars -except sl_synch_handles simulinkStaticGUI
clc

% Add path to local source code
addpath('./src/')

%% GENERAL SIMULATION INFO
%
% If you are simulating the robot with Gazebo, remember that it is required
% to launch Gazebo as follows:
%
%     gazebo -slibgazebo_yarp_clock.so
%
% and properly set the environmental variable YARP_ROBOT_NAME in the .bashrc.

% Simulation time in seconds. For long simulations, put an high number
% (not inf) for allowing automatic code generation
Config.SIMULATION_TIME = 70;

% Controller period [s]
Config.tStep           = 0.01;

%% PRELIMINARY CONFIGURATION
%
% DEMO_TYPE: defines the kind of demo that will be performed.
%
% 'YOGA': the robot will perform the YOGA++ demo (highly dynamic movements
%         while balancing on one foot and two feet)
%
% 'COORDINATOR': the robot can either balance on two feet or move from left
%                to right follwing a desired center-of-mass trajectory.
%
DEMO_TYPE                     = 'YOGA';

% Config.SCOPES: debugging scopes activation
Config.SCOPES_WRENCHES        = true;
Config.SCOPES_GAIN_SCHE_INFO  = true;
Config.SCOPES_MAIN            = true;
Config.SCOPES_QP              = true;

% DATA PROCESSING
%
% Save the Matlab workspace after stopping the simulation
Config.SAVE_WORKSPACE         = false;

% Verify that the integration time has been respected during the simulation
Config.CHECK_INTEGRATION_TIME = true;
Config.PLOT_INTEGRATION_TIME  = false;

% Run robot-specific configuration parameters
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/configRobot.m'));
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/configStateMachine.m'));
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/gainsAndReferences.m'));

% Deactivate/activate the internal coordinator
if strcmpi(DEMO_TYPE, 'COORDINATOR')

    Config.COORDINATOR_DEMO = true;

elseif strcmpi(DEMO_TYPE, 'YOGA')

    Config.COORDINATOR_DEMO = false;
end

author: robotology

repo: whole-body-controllers


text: ezplot('tan(x)')

author: Eurus-Holmes


text:
classdef TopOptTestsSuite < handle

    methods

        function obj = TopOptTestsSuite()
            warning('off', 'MATLAB:structOnObject')
            % testFastDisplacement, testMacro, testMicro
            results = runtests("TopOptTests","ProcedureName","testFastDisplacement", 'Debug', true);
%             results = runtests("TopOptTests","ProcedureName","testAnalyticVsRegularizedPerimeter", 'Debug', true);
%             results = runtests("TopOptTests", 'Debug', true);
            table(results)
            warning('on', 'MATLAB:structOnObject')
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function S = spm_cfg_eeg_montage
% configuration file for reading montage files
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_cfg_eeg_montage.m 6983 2017-01-06 17:25:22Z vladimir $

D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the EEG mat file.'};

montagefile = cfg_files;
montagefile.tag = 'montagefile';
montagefile.name = 'Montage file name';
montagefile.filter = 'mat';
montagefile.num = [1 1];
montagefile.help = {'Select a montage file.'};

keepothers = cfg_menu;
keepothers.tag = 'keepothers';
keepothers.name = 'Keep other channels';
keepothers.labels = {'Yes', 'No'};
keepothers.values = {1, 0};
keepothers.val = {0};
keepothers.help = {'Specify whether you want to keep channels that are not contributing to the new channels'};

montage = cfg_branch;
montage.tag = 'montage';
montage.name = 'Load montage from file';
montage.val = {montagefile, keepothers};

montagename         = cfg_entry;
montagename.tag     = 'montagename';
montagename.name    = 'Montage name';
montagename.help    = {'Specify the name of existing montage'};
montagename.strtype = 's';
montagename.num     = [1 Inf];

montagenum         = cfg_entry;
montagenum.tag     = 'montagenum';
montagenum.name    = 'Montage index';
montagenum.help    = {'Specify the index of existing montage'};
montagenum.strtype = 'w';
montagenum.num     = [1 1];

montspec = cfg_choice;
montspec.tag = 'montspec';
montspec.name = 'Montage specification';
montspec.values = {montage, montagename, montagenum};
montspec.val = {montage};
montspec.help = {'Choose how to specify the montage'};

blocksize = cfg_entry;
blocksize.tag = 'blocksize';
blocksize.name = 'Block size';
blocksize.strtype = 'r';
blocksize.val = {655360};
blocksize.num = [1 1];
blocksize.help = {'size of blocks used internally to split large files default ~20Mb'};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the montaged ',...
    'dataset (if generated). Default prefix is ''M''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'M'};

write = cfg_branch;
write.tag = 'write';
write.name = 'Write';
write.val = {montspec, blocksize, prefix};
write.help = {'Apply montage and write out to a new dataset.'};

online = cfg_branch;
online.tag = 'switch';
online.name = 'Switch';
online.val = {montspec};
online.help = {'Switch to new online montage without writing out a new file'};

add = cfg_branch;
add.tag = 'add';
add.name = 'Add';
add.val = {montage};
add.help = {'Add montage to the set of online montages without applying'};

clr = cfg_const;
clr.tag = 'clear';
clr.name = 'Clear';
clr.val  = {1};
clr.help = {'Clear all online montages and revert to the original state.'};

mode = cfg_choice;
mode.tag = 'mode';
mode.name = 'Mode';
mode.values = {write, online, add, clr};
mode.val = {write};
mode.help = {'Choose between writing a new dataset or online montage operation'};

S = cfg_exbranch;
S.tag = 'montage';
S.name = 'Montage';
S.val = {D, mode};
S.help = {'Apply a montage (linear transformation) to EEG/MEG data.'};
S.prog = @eeg_montage;
S.vout = @vout_eeg_montage;
S.modality = {'EEG'};

function out = eeg_montage(job)
% construct the S struct
S.D = job.D{1};
S.mode = char(fieldnames(job.mode));
switch S.mode
    case 'write'
        S.blocksize = job.mode.write.blocksize;
        S.prefix = job.mode.write.prefix;
        switch char(fieldnames(job.mode.write.montspec))
            case 'montagename'
                S.montage = job.mode.write.montspec.montagename;
            case 'montagenum'
                S.montage = job.mode.write.montspec.montagenum;
            case 'montage'
                S.montage = char(job.mode.write.montspec.montage.montagefile);
                S.keepothers = job.mode.write.montspec.montage.keepothers;
        end
    case 'switch'
        switch char(fieldnames(job.mode.switch.montspec))
            case 'montagename'
                S.montage = job.mode.switch.montspec.montagename;
            case 'montagenum'
                S.montage = job.mode.switch.montspec.montagenum;
            case 'montage'
                S.montage = char(job.mode.switch.montspec.montage.montagefile);
                S.keepothers = job.mode.switch.montspec.montage.keepothers;
        end
    case 'add'
        S.montage = char(job.mode.add.montage.montagefile);
        S.keepothers = job.mode.add.montage.keepothers;
end
out.D = spm_eeg_montage(S);
out.Dfname = {fullfile(out.D.path, out.D.fname)};

function dep = vout_eeg_montage(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Montaged data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Montaged Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});