text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Launch code for sinusoidal input                 %
%       code by Fabrizio Conso, university of pavia, student       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                launch code example                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                       %
 format short;
 f_s=1;            % normalized sampling frequency
 f_in=0.4;         % normalized signal frequency
 f_bw=5;                    % DAC's bandwidth
 sr=2;                      % DAC's slew-rate
 nbit=12;                   % converter bits
 N=2^12;                    % fft points
 t=[1:(1/f_s):N*(1/f_s)];   % discrete time array
 bw=0.5;
 fin=antismear(f_in,f_s,N);
 for i=1:N         % cycle for collect ADC outputs
     input=0.5+0.5*sin(2*pi*fin*t(i));
     counter2=Approx_2(input,nbit,f_bw,sr,f_s);
     counter=Approx(input,nbit,f_bw,sr,f_s);
     y(i)=counter;
     z(i)=counter2;
 end

 z=z-mean(z);                % eliminating DC tone
 y=y-mean(y);                % eliminating DC tone
 w=gausswin(N);              % window
 f=fin/f_s;                  % Normalized signal frequency
 fB=N*(bw/f_s);              % Base-band frequency bins
 [snr2,ptot2]=calcSNR(z(1:N),f,fB,w',N);
 ptot2=ptot2-max(ptot2);
 [snr,ptot]=calcSNR(y(1:N),f,fB,w',N);
 ptot=ptot-max(ptot);

figure(1);
clf;
plot(linspace(0,f_s/2,N/2), ptot(1:N/2), 'r');
grid on;
title('Output PSD without correction')
xlabel('Frequency [Hz]')
ylabel('PSD [dB]')
axis([0.3 f_s/2 -120 0]);

figure(2);
clf;
plot(linspace(0,f_s/2,N/2), ptot2(1:N/2), 'r');
grid on;
title('Output PSD with correction')
xlabel('Frequency [Hz]')
ylabel('PSD [dB]')
axis([0.3 f_s/2 -120 0]);
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
% [M_F, V_F] = GPPDE(X_Y, Y, S2_Y, THETA, X_G, G, ALPHA, D, X_F)

% Copyright (c) 2010 Jaakko Luttinen

% TODO: Each equation as a struct/cell: {X, y, K, alpha, D}
%
% X : input
% y : observations
% K : observation covariance matrix
% alpha : coefficients in (differential) equation
% D : differential operators in (differential) equation
%
% Then, noisy function observations are:
% { X_y, y, s2_y, 1, 0 }
% Noiseless derivative observations:
% { X_dy, dy, 1e-8, 1, 1 }
% E.g., heat equation:
% { X_g, zeros(..), 1e-8, [1; -a], [1 0; 0 2] }
%
% Actually, you should allow s2 to be a covariance matrix.  That way you
% could use other unknown functions in equations.. Or maybe even more
% generally, s2 could be a covariance function, thus its hyperparameters
% could be learnt simultaneously. :)
%
% And you can give several of these.
%
% And same for the predictions!

% TODO: Learn hyperparameters.

function [m_f, V_f] = gppde(varargin)

%
% Generate full joint covariance matrix for the equations
%

% Identity matrix for data noise
I_y = speye(size(X_y,1));
alpha0 = 1;
D0 = zeros(1, size(X_y,2));

% Variables for PDE numerical inaccuracy
I_g = speye(length(g));
s2_g = 1e-6;

% Joint data
yg = [y; g];

%
% DEBUG STUFF
%
% $$$ func = @(x) covfunc(x(:)', X_g(1,:), theta, alpha0,D0, alpha, D);
% $$$ mycheckgrad(func,

% Joint prior covariance matrix for observations and PDE
K_g = covfunc(X_g, X_g, theta, alpha, D, alpha, D) + s2_g*I_g;
K_y = covfunc(X_y, X_y, theta, alpha0,D0, alpha0,D0) + s2_y*I_y;
K_y_g = covfunc(X_y, X_g, theta, alpha0,D0, alpha, D);
K_yg = [K_y, K_y_g; K_y_g', K_g];

% $$$ figure
% $$$ imagesc(K_yg)

L_yg = chol(K_yg, 'lower');

% Posterior covariance for predictive function values
K_y_f = covfunc(X_y, X_f, theta, alpha0,D0, alpha0,D0);
K_g_f = covfunc(X_g, X_f, theta, alpha,D, alpha0,D0);
K_yg_f = [K_y_f; K_g_f];
K_f = covfunc(X_f, X_f, theta, alpha0,D0, alpha0,D0);
Z = linsolve_tril(L_yg, K_yg_f);
V_f = K_f - Z' * Z;
m_f = Z' * linsolve_tril(L_yg, yg);

function K = covfunc(X1, X2, theta, alpha1, D1, alpha2, D2)

M = length(alpha1);
N = length(alpha2);
P = size(X1,2);

% For now, distances with unit length scales
D = bsxfun(@minus, ...
           reshape(X1, [size(X1,1),1,P]), ...
           reshape(X2, [1,size(X2,1),P]));

% Unit scale
K0 = theta(1)^2 * exp(-0.5/theta(2)^2*sum(D.^2,3));

K = 0;
for m=1:M
  for n=1:N
    % [c,E] = tmp(D1(m,:)+D2(n,:));

    % Evaluate exponents (E) and coefficients (c)
    b = D1(m,:)+D2(n,:) ;
    A = ntuples(floor(b/2));
    B = repmat(b, [size(A,1),1]);
    E = bsxfun(@minus, b, 2*A) ;
% $$$     c = (-1)^(sum(b)) * (-1).^(sum(A,2)) .* prod(npairsk(B, A), 2) % debug test
    c = (-1)^(sum(D1(m,:))) * (-1).^(sum(A,2)) .* prod(npairsk(B, A), 2) ;

    for l=1:length(c)
      Z = prod(bsxfun(@power, D, reshape(E(l,:), [1,1,P])), 3);
      Z = c(l) * Z;
      thetaexp = sum(0.5*(E(l,:)+b));
      K = K + alpha1(m)*alpha2(n) * theta(2)^(-2*thetaexp) * Z;
    end
  end
end
K = K.*K0;

% $$$ function [c, E] = term_coeff(d)
% $$$
% $$$ Z = ntuples(floor(d/2))
% $$$ D = repmat(d, [size(Z,1),1]);
% $$$ E = bsxfun(@minus, d, 2*Z);
% $$$ c = (-1).^(sum(Z,2)) .* prod(npairsk(D, Z), 2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function logpdf = logmvtpdf(t,mu,f,Sigma);
% function logpdf = logmvtpdf(t,mu,f,Sigma);
% log pdf of multivariate t-student dist.
% t : D by N

[d,n] = size(t);

c = gammaln((d+f)*0.5) - (d*0.5)*log(f*pi) - gammaln(f*0.5) - 0.5*detln(Sigma);
diff = t - repmat(mu,1,n); % d by n
logpdf = c - (f+d)*0.5 * log(1 + sum(diff.*(inv(f*Sigma)*diff),1));

% Local Variables: ***
% mode: matlab ***
% End: ***

repo: REST


text:
%% Pose Estimation
%
% In this sample, we learn to exploit calib3d module to create some 3D effects
% in images from a calibrated camera.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d7/d53/tutorial_py_pose.html>
% * <https://docs.opencv.org/3.2.0/dc/d43/tutorial_camera_calibration_square_chess.html>
%

%%
% In the camera calibration demo, we have found the camera matrix, distortion
% coefficients etc. Given a pattern image, we can utilize the above
% information to calculate its pose, or how the object is situated in space,
% like how it is rotated, how it is displaced etc. For a planar object, we can
% assume Z=0, such that, the problem now becomes how camera is placed in space
% to see our pattern image. So, if we know how the object lies in the space,
% we can draw some 2D diagrams in it to simulate the 3D effect.
%
% Our problem is, we want to draw our 3D coordinate axis (X, Y, Z axes) on our
% chessboard's first corner. X axis in blue color, Y axis in green color and
% Z axis in red color. So in-effect, Z axis should feel like it is
% perpendicular to our chessboard plane.
%

%%
% First, let's load the camera matrix and distortion coefficients from the
% previous calibration result.
fname = fullfile(tempdir(), 'calibration_chessboard.yml');
assert(exist(fname, 'file') == 2, ...
    'calibration result not found, run calibration_demo.m first');
fs = cv.FileStorage(fname);
display(fs)

%%
% Then as in previous case, we create object points
% (planar 3D points of chessboard corners)
[X,Y] = ndgrid(1:fs.board_width, 1:fs.board_height);
objPts = ([X(:), Y(:)] - 1) * fs.square_size;
objPts(:,3) = 0;

%%
% Create axis points. Axis points are points in 3D space for drawing the axis.
% We draw axis of length 3 (units will be in terms of chess square size since
% we calibrated based on that size). So our X axis is drawn from (0,0,0) to
% (3,0,0), so for Y axis. For Z axis, it is drawn from (0,0,0) to (0,0,-3).
% Negative denotes it is drawn towards the camera.
p3d = [3 0 0; 0 3 0; 0 0 -3] * fs.square_size;

%%
% Now, as usual, we loop over each image
files = cv.glob(fullfile(mexopencv.root(), 'test', 'left*.jpg'));
N = numel(files);
for i=1:N
    % load image
    img = cv.imread(files{i}, 'Color',true);

    % search for 9x6 grid
    [imgPts, found] = cv.findChessboardCorners(img, ...
        [fs.board_width, fs.board_height]);
    imgPts = cat(1, imgPts{:});
    if found
        % if found, refine it with subcorner pixels.
        gray = cv.cvtColor(img, 'RGB2GRAY');
        imgPts = cv.cornerSubPix(gray, imgPts, 'WinSize',[11 11], ...
            'Criteria',struct('type','Count+EPS', 'maxCount',30, 'epsilon',0.1));
    end

    % calculate the rotation and translation vectors
    [rvecs, tvecs] = cv.solvePnP(objPts, imgPts, fs.camera_matrix, ...
        'DistCoeffs',fs.distortion_coefficients);

    % use transformation matrices to project "axis points" to the image plane.
    % In simple words, we find the points on image plane corresponding to each
    % of (3,0,0),(0,3,0),(0,0,3) in 3D space.
    p2d = cv.projectPoints(p3d, rvecs, tvecs, fs.camera_matrix, ...
        'DistCoeffs',fs.distortion_coefficients);

    % draw lines from the first corner to each of projected axis points.
    % Notice that each axis is 3 squares long.
    clrs = [0 0 1 0; 0 1 0 0; 1 0 0 0] * 255;
    img = cv.line(img, imgPts([1 1 1],:), p2d, 'Colors',clrs, 'Thickness',3);
    imshow(img), title(sprintf('%02d / %02d', i, N))
    pause(1)
end

%%
% We can also modify the axis points above to draw a cube instead of axis
% vectors. We draw ground floor in green, pillars in blue, top layer in red.
%
% This can be extended to render more complicated objects such as for
% augmented reality applications.
%
p3d = [0 0 0; 0 3 0; 3 3 0; 3 0 0;
    0 0 -3; 0 3 -3; 3 3 -3; 3 0 -3] * fs.square_size;
p2d = cv.projectPoints(p3d, rvecs, tvecs, fs.camera_matrix, ...
    'DistCoeffs',fs.distortion_coefficients);
img = cv.imread(files{N}, 'Color',true);
img = cv.drawContours(img, p2d(1:4,:), 'Color',[0 255 0], 'Thickness','Filled');
img = cv.line(img, p2d(1:4,:), p2d(5:8,:), 'Color',[0 0 255], 'Thickness',3);
img = cv.drawContours(img, p2d(5:8,:), 'Color',[255 0 0], 'Thickness',3);
figure, imshow(img), title(sprintf('%02d / %02d', N, N))

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function i4vec_permute_test ( )

%*****************************************************************************80
%
%% I4VEC_PERMUTE_TEST tests I4VEC_PERMUTE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
  n = 12;
  b = 0;
  c = n;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_PERMUTE_TEST\n' );
  fprintf ( 1, '  I4VEC_PERMUTE reorders an integer vector\n' );
  fprintf ( 1, '  according to a given permutation.\n' );
  fprintf ( 1, '  Using initial random number seed = %d\n', seed );

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  A, before rearrangement:' );

  p = perm1_uniform ( n, seed );

  i4vec_print ( n, p, '  Permutation vector P:' );

  a = i4vec_permute ( n, a, p );

  i4vec_print ( n, a, '  A, after rearrangement:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function month_name = month_to_month_name_republican ( m )

%*****************************************************************************80
%
%% MONTH_TO_MONTH_NAME_REPUBLICAN returns the name of a Republican month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the month index.
%
%    Output, string MONTH_NAME, the month name.
%
  name = [ ...
    'Vendemaire    '; 'Brumaire      '; 'Frimaire      '; 'Nivose        '; ...
    'Pluviose      '; 'Ventose       '; 'Germinal      '; 'Floreal       '; ...
    'Prairial      '; 'Messidor      '; 'Thermidor     '; 'Fructidor     '; ...
    'Sansculottides' ];

  if ( m < 1 || 13 < m )

    month_name = '??????????????';

  else

    month_name = name(m,:);

  end

  return
end


text:
function [nconds,cforecast_record,cforecast_estimates]=...
    panel2cf(N,n,m,p,k,q,cfconds,cfshocks,cfblocks,data_endo_a,data_exo_a,data_exo_p,It,Bu,Fperiods,const,beta_gibbs,D_record,gamma_record,CFt,Fband)

% initiate the cell recording the Gibbs sampler draws
cforecast_record={};
cforecast_estimates={};

% because conditional forecasts can be computed for many (potentially all) units, loop over units
for ii=1:N
% check wether there are any conditions on unit ii
temp=cfconds(:,:,ii);
nconds(ii,1)=numel(temp(cellfun(@(x) any(~isempty(x)),temp)));

   % if there are conditions
   if nconds(ii,1)~=0
   % prepare the elements for conditional forecast estimation, depending on the type of conditional forecasts
   temp1=cfconds(:,:,ii);
      if CFt==1
      temp2={};
      temp3=[];
      elseif CFt==2
      temp2=cfshocks(:,:,ii);
      temp3=cfblocks(:,:,ii);
      end
   % run the Gibbs sampler for unit ii
   cforecast_record(:,:,ii)=bear.cforecast(data_endo_a(:,:,ii),data_exo_a,data_exo_p,It,Bu,Fperiods,temp1,temp2,temp3,CFt,const,beta_gibbs,D_record,gamma_record,n,m,p,k,q);
   % then obtain point estimates and credibility intervals
   cforecast_estimates(:,:,ii)=bear.festimates(cforecast_record(:,:,ii),n,Fperiods,Fband);

   % if there are no conditions, return empty elements
   elseif nconds(ii,1)==0
   cforecast_record(:,:,ii)=cell(n,1);
   cforecast_estimates(:,:,ii)=cell(n,1);
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function value = r4_besi1e ( x )

%*****************************************************************************80
%
%% R4_BESI1E: exponentially scaled Bessel function I of order 1 of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the exponentially scaled Bessel function I
%    of order 1 of X.
%
  persistent ai1cs
  persistent ai12cs
  persistent bi1cs
  persistent ntai1
  persistent ntai12
  persistent nti1
  persistent xmin
  persistent xsml

  if ( isempty ( nti1 ) )

    ai1cs = [ ...
     -0.02846744181881479E+00, ...
     -0.01922953231443221E+00, ...
     -0.00061151858579437E+00, ...
     -0.00002069971253350E+00, ...
      0.00000858561914581E+00, ...
      0.00000104949824671E+00, ...
     -0.00000029183389184E+00, ...
     -0.00000001559378146E+00, ...
      0.00000001318012367E+00, ...
     -0.00000000144842341E+00, ...
     -0.00000000029085122E+00, ...
      0.00000000012663889E+00, ...
     -0.00000000001664947E+00, ...
     -0.00000000000166665E+00, ...
      0.00000000000124260E+00, ...
     -0.00000000000027315E+00, ...
      0.00000000000002023E+00, ...
      0.00000000000000730E+00, ...
     -0.00000000000000333E+00, ...
      0.00000000000000071E+00, ...
     -0.00000000000000006E+00 ]';
    ai12cs = [ ...
      0.02857623501828014E+00, ...
     -0.00976109749136147E+00, ...
     -0.00011058893876263E+00, ...
     -0.00000388256480887E+00, ...
     -0.00000025122362377E+00, ...
     -0.00000002631468847E+00, ...
     -0.00000000383538039E+00, ...
     -0.00000000055897433E+00, ...
     -0.00000000001897495E+00, ...
      0.00000000003252602E+00, ...
      0.00000000001412580E+00, ...
      0.00000000000203564E+00, ...
     -0.00000000000071985E+00, ...
     -0.00000000000040836E+00, ...
     -0.00000000000002101E+00, ...
      0.00000000000004273E+00, ...
      0.00000000000001041E+00, ...
     -0.00000000000000382E+00, ...
     -0.00000000000000186E+00, ...
      0.00000000000000033E+00, ...
      0.00000000000000028E+00, ...
     -0.00000000000000003E+00 ]';
    bi1cs = [ ...
     -0.001971713261099859E+00, ...
      0.40734887667546481E+00, ...
      0.034838994299959456E+00, ...
      0.001545394556300123E+00, ...
      0.000041888521098377E+00, ...
      0.000000764902676483E+00, ...
      0.000000010042493924E+00, ...
      0.000000000099322077E+00, ...
      0.000000000000766380E+00, ...
      0.000000000000004741E+00, ...
      0.000000000000000024E+00 ]';

    nti1 = r4_inits ( bi1cs, 11, 0.1 * r4_mach ( 3 ) );
    ntai1 = r4_inits ( ai1cs, 21, 0.1 * r4_mach ( 3 ) );
    ntai12 = r4_inits ( ai12cs, 22, 0.1 * r4_mach ( 3 ) );
    xmin = 2.0 * r4_mach ( 1 );
    xsml = sqrt ( 8.0 * r4_mach ( 3 ) );

  end

  y = abs ( x );

  if ( x == 0.0 )
    value = 0.0;
  elseif ( y <= xmin )
    value = 0.0;
  elseif ( y <= xsml )
    value = 0.5 * x;
    value = exp ( - y ) * value;
  elseif ( y <= 3.0 )
    value = x * ( 0.875 ...
      + r4_csevl ( y * y / 4.5 - 1.0, bi1cs, nti1 ) );
    value = exp ( - y ) * value;
  elseif ( y <= 8.0 )
    value = ( 0.375 + r4_csevl ( ( 48.0 / y - 11.0 ) / 5.0, ...
      ai1cs, ntai1) ) / sqrt ( y );
    if ( x < 0.0 )
      value = - value;
    end
  else
    value = ( 0.375 + r4_csevl ( 16.0 / y - 1.0, ai12cs, ntai12 ) ) ...
      / sqrt ( y );
    if ( x < 0.0 )
      value = - value;
    end
  end

  return
end


text:
function handleExecErrorsWarnings(varargin)
%handleExecErrors Summary of this function goes here
%   Detailed explanation goes here
    global warn_error_label_handles warn_error_slider_handle;
    warnErrorLbls = warn_error_label_handles;
    hSlider = warn_error_slider_handle;

    execErrors = getExecutionErrors();
    execWarns = getExecutionWarnings();

    numWarnError = length(execErrors) + length(execWarns);

    if(~isempty(varargin))
        updateSlider = varargin{1};
    else
        updateSlider = true;
    end

    if(numWarnError <= 6)
        set(hSlider,'Max',1.0);
        set(hSlider,'Min',0.0);
        set(hSlider,'Value',1.0);
        set(hSlider,'Enable','off');
    else
        if(updateSlider)
            set(hSlider,'Min',0.0);
            set(hSlider,'Max',numWarnError-6);
            set(hSlider,'Value',numWarnError-6);
            set(hSlider,'SliderStep',[1/(numWarnError-6),1.0]);
            set(hSlider,'Enable','on');
        end
    end

    lblOffset = round(get(hSlider,'Max') - get(hSlider,'Value'));
    set(hSlider,'Value',round(get(hSlider,'Value')));

    lblUseCnt = 0;
    for(i=1+lblOffset:length(execErrors))
        if(lblUseCnt>=6)
            break;
        end
        errorStr = execErrors{i};
        styleLabelAsError(warnErrorLbls(lblUseCnt+1), errorStr);
        drawnow;
        lblUseCnt = lblUseCnt+1;
    end

    if(lblOffset >= length(execErrors))
        lblOffset = lblOffset - length(execErrors);
    else
        lblOffset = 0;
    end

    for(i=1+lblOffset:length(execWarns)) %#ok<*NO4LP>
        if(lblUseCnt>=6)
            break;
        end
        warnStr = execWarns{i};
        styleLabelAsWarn(warnErrorLbls(lblUseCnt+1), warnStr);
        drawnow;
        lblUseCnt = lblUseCnt+1;
    end

    if(numWarnError == 0)
        styleLabelAsValid(warnErrorLbls(1), 'No errors or warnings found.');
        drawnow;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Test file for ADCHEBFUN times

function pass = test_power

% Initialise pass vector
pass = zeros(3, 5);

% Tolerance for Taylor testing
tolOrder = 1e-2;
tolDiff = 1e-14;

% The function handle we're working with:
func = @power;

% Compare values
[err, lin] = adchebfun.valueTestingBinary(func);

% Confirm that the error returned is zero
pass(1, :) = ( err < tolDiff );

% Taylor testing
[order1, order2] = adchebfun.taylorTestingBinary(func);

% We expect all elements of ORDER1 to be close to 1. Since the methods being
% tested in this all nonlinear, ORDER2 should have entries close to 2.

% Check whether we get the expected results for linear and nonlinear
% operations.
pass(2,:) = ( (max(abs(order1 - 1)) < tolOrder) & ...
              (max(abs(order2 - 2)) < tolOrder) );

% Linearity checking. All operations should be detected as nonlinear.
pass(3,:) = ( lin == 0 );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [feature_matrix, labels] = ...
    get_speaker_features_and_labels(db, field_name)
    % Extracts features and labels from all the samples corresponding
    % to the specified speakers.
    % db - Database of samples generated by gendb()
    % speakers - List of speakers for which to extract features
    % feature_matrix - Feature matrix for input to Diffusion Maps
    % labels - Label for each column of feature_matrix

    unique_field_values = get_field(db, field_name);
    if (strcmp(field_name, 'digit'))
       unique_field_values = cell2mat(unique_field_values);
       unique_field_values(:, 2) = '*';
       unique_field_values = mat2cell(unique_field_values, ...
           ones(size(unique_field_values,1), 1));
       unique_field_values = unique(unique_field_values);
    end

    feature_matrix = [];
    labels = [];
    total_entries = 0; % counter of handled entries

    progressbar;
    NUM_OF_UNIQUE_VALUES = length(unique_field_values);
    for k = 1:NUM_OF_UNIQUE_VALUES
        speaker_db = filterdb(db, field_name, unique_field_values{k});

        % extract features for all speaker samples
        func = @calc_mfcc;
%         func = @calc_stft;
        [features, num_of_success] = extract_features(speaker_db, func);

        if num_of_success > 0
            labels = [labels; ones(num_of_success, 1) * k];
            feature_matrix = [feature_matrix features];
            total_entries = total_entries + num_of_success;
        end

        progressbar(k/NUM_OF_UNIQUE_VALUES);
    end

    fprintf('Handled %d entries\n', total_entries);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%kkcolorspace 'Change the color space attribute on a data object'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kcolorspace.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% MultiChoice: type 'Image Colorspace', default: 1: 'Type of colorspace to be applied to the input image'
%    Choices are:
%   1: 'Greyscale'
%   2: 'RGB Colorspace'
%   3: 'CMY Colorspace'
%   4: 'YIQ Colorspace'
%   5: 'HSV Colorspace'
%   6: 'HLS Colorspace'
%   7: 'IHS Colorspace'
%   8: 'XYZ Colorspace'
%   9: 'UVW Colorspace'
%   10: 'UCSUVW Colorspace'
%   11: 'UCSSOW Colorspace'
%   12: 'UCSLab Colorspace'
%   13: 'UCSLuv Colorspace'
%   14: 'User Defined Colorspace'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkcolorspace(i, {'i','';'type',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kcolorspace - Change the color space attribute on a data object
%
%  DESCRIPTION
% .I kcolorspace
% is used to change the setting of the KCOLOR_COLORSPACE attribute on a data
% object to any of the allowed values. Note that changing the attribute using
% kcolorspace does not modify the data (say by actually performing an RGB to
% HLS converssion if the KCOLOR_COLORSPACE is changed from RGB to HLS). Instead,
% changing the attribute will cause the data in the object to be interpreted
% as to already be in the new color space.
%
%
%
%  EXAMPLES
% Consider an W=512 H=512 E=3 data object, with each of the WxH planes obtained
% from a color scanner using R, G, and B filters. Simply inserting these data
% planes into a "blank" data object will cause the data to be treated as
% three planes of greyscale data. To cause the data to be interpreted as RGB
% data (for those operations that perform color processing, including the
% display tools), use kcolorspace to change the setting to RGB.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkcolorspace(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkcolorspace(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'type', 1;'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function centroid = polyhedronCentroid(vertices, faces) %#ok<INUSD>
%POLYHEDRONCENTROID Compute the centroid of a 3D convex polyhedron.
%
%   CENTRO = polyhedronCentroid(V, F)
%   Computes the centroid (center of mass) of the polyhedron defined by
%   vertices V and faces F.
%   The polyhedron is assumed to be convex.
%
%   Example
%     % Creates a polyhedron centered on origin, and add an arbitrary
%     % translation
%     [v, f] = createDodecahedron;
%     v2 = bsxfun(@plus, v, [3 4 5]);
%     % computes the centroid, that should equal the translation vector
%     centroid = polyhedronCentroid(v2, f)
%     centroid =
%         3.0000    4.0000    5.0000
%
%
%   See also
%   meshes3d, meshVolume, meshSurfaceArea, polyhedronMeanBreadth
%

% ------
% Author: David Legland
% E-mail: david.legland@nantes.inra.fr
% Created: 2012-04-05, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% compute set of elementary tetrahedra
DT = delaunayTriangulation(vertices);
T = DT.ConnectivityList;

% number of tetrahedra
nT  = size(T, 1);

% initialize result
centroid = zeros(1, 3);
vt = 0;

% Compute the centroid and the volume of each tetrahedron
for i = 1:nT
    % coordinates of tetrahedron vertices
    tetra = vertices(T(i, :), :);

    % centroid is the average of vertices.
    centi = mean(tetra);

    % compute volume of tetrahedron
    vol = det(tetra(1:3,:) - tetra([4 4 4],:)) / 6;

    % add weighted centroid of current tetraedron
    centroid = centroid + centi * vol;

    % compute the sum of tetraedra volumes
    vt = vt + vol;
end

% compute by sum of tetrahedron volumes
centroid = centroid / vt;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Test file for singfun/plus.m
% This has been adapted from bndfun test for plus.
function pass = test_plus(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(666);
x = 2 * rand(100, 1) - 1;

% A random number to use as an arbitrary additive constant.
alpha = -0.194758928283640 + 0.075474485412665i;

%%
% Check operation in the case of empty arguments.
f = singfun();
data.exponents = [-1, 0];
g = singfun(@(x) 1./(1+x), data, pref);
pass(1) = (isempty(f + f) && isempty(f + g) && isempty(g + f));

% Addition of smooth SINGFUNs should not return a SINGFUN
f = singfun(@(x) sin(x));
g = singfun(@(x) cos(x));
pass(2) = ~isa(f+g, 'singfun');

% SMOOTHFUN + SINGFUN
f = smoothfun.constructor(@(x) sin(x));
g = singfun(@(x) cos(x));
pass(3) = isa(f + g, 'smoothfun') && isa(g + f, 'smoothfun');
%%
% Check addition with scalars.

fh = @(x) 1./((1+x).*(1-x));
data.exponents = [-1, -1];
f = singfun(fh, data, pref);
pass(4:5) = test_add_function_to_scalar(f, fh, alpha, x);

%%
% Check addition of two singfun objects.

fh = @(x) zeros(size(x));
data.exponents = [];
data.singType = {'none', 'none'};
f = singfun(fh, data, pref);
pass(6:7) = test_add_function_to_function(f, fh, f, fh, x);

fh = @(x) sin(pi*x)./(1-x);
f = singfun(fh, [], pref);

gh = @(x) cos(pi*x)./(1-x);
g = singfun(gh, [], pref);
pass(8:9) = test_add_function_to_function(f, fh, g, gh, x);

gh = @(x) cos(1e2*x);
g = singfun(gh, [], pref);
pass(10:11) = test_add_function_to_function(f, fh, g, gh, x);

gh = @(t) sinh(t*exp(2*pi*1i/6));
g = singfun(gh, [], pref);
pass(12:13) = test_add_function_to_function(f, fh, g, gh, x);

%%
% Check that direct construction and PLUS give comparable results.

f = singfun(@(x) x, [], pref);
g = singfun(@(x) cos(x) - 1, [], pref);
h1 = f + g;
h2 = singfun(@(x) x + cos(x) - 1, [], pref);
tol = 10*max(get(h1, 'vscale')*eps, ...
    get(h2, 'vscale')*eps);
pass(14) = norm(feval(h1, x) - feval(h2, x), inf) < tol;

% Check that plus handles "small" results appropriately in the case of a
% non-integer exponent difference.
warnState = warning('off', 'CHEBFUN:SINGFUN:plus:exponentDiff');

pref_fixed = pref;
pref_fixed.fixedLength = 256;
op = @(x) ((x + 1)/2).^pi;
data.exponents = [pi - 3, 0];
f = singfun(op, data, pref);
g = singfun(op, [], pref_fixed);
h = f - g;
pass(15) = get(h, 'ishappy') && (length(h) < 1024);

warning(warnState);

end

% Test the addition of a SINGFUN F, specified by Fh, to a scalar C using
% a grid of points X in [a  b] for testing samples.
function result = test_add_function_to_scalar(f, fh, c, x)
    g1 = f + c;
    g2 = c + f;
    result(1) = isequal(g1, g2);
    g_exact = @(x) fh(x) + c;
    result(2) = norm(feval(g1, x) - g_exact(x), inf) < 1e3*eps;
end

% Test the addition of two SINGFUN objects F and G, specified by FH and
% GH, using a grid of points X in [-1  1] for testing samples.
function result = test_add_function_to_function(f, fh, g, gh, x)
    h1 = f + g;
    h2 = g + f;
    result(1) = isequal(h1, h2);
    h_exact = @(x) fh(x) + gh(x);
    result(2) = norm(feval(h1, x) - h_exact(x), inf) <= 2e3*eps;
end


text:
function kos_map = obtain_kos_map(lines_refined,img_size)
    lines = lines_refined;
    init_map = double(zeros(img_size(1),img_size(2)));
    shortest = min(lines(5,:));
    longest = max(lines(5,:));
    lines(5,:) = (lines(5,:)-shortest) ./ (longest-shortest);
    for lidx = 1 : size(lines,2)
        y = sort(lines([1,2],lidx));
        x = sort(lines([3,4],lidx));
        expand_factor = 0.2; % expand 20%
        x_min = round( x(1)*(1-expand_factor) );
        x_max = round( x(2)*(1+expand_factor) );
        y_min = round( y(1)*(1-expand_factor) );
        y_max = round( y(2)*(1+expand_factor) );
        [x_min,x_max,y_min,y_max] = handle_cross_boundary(x_min,x_max,y_min,y_max,img_size);
        s = lines(5,lidx);
        init_map(x_min:x_max,y_min:y_max) = init_map(x_min:x_max,y_min:y_max) + s;
    end
    w = ones(5,5) / 25;
    kos_map = imfilter(mat2gray(init_map),w);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
vI = cat(3, 85, 80, 217) / 255;

vPhotoshopValues = [30; 98; 51; 73; 5; 53];

vCoeffValues = (vPhotoshopValues - 50) ./ 50;
vO = ApplyBlackWhiteFilter(vI, vCoeffValues);
vO = max(min(vO, 1), 0);

round(255 * vI(:))
round(255 * vO(:))

author: RoyiAvital

repo: StackExchangeCodes


text:
%
% NAME: mrAlign (ver 3)
% AUTHOR: DJH
% DATE: 8/2001
% PURPOSE:
%	Main routine for menu driven matlab program to view and analyze
%	volumes of anatomical and functional MRI data.
% HISTORY:
%	Dozens of people contributed to the original mrAlign. It was written
%   using scripts in matlab 4 without any data structures. This is an
%   attempt to make a fresh start.
%
disp('mrAlign (version 3)');

% Global Variables

global HOMEDIR
HOMEDIR = pwd;

% Matlab controls for menu driven program
global volinc voldec
global volslislice volslimin1 volslimax1 volslimin2 volslimax2 volslicut
global interpflag
global volcmap 	% Color map
global checkAutoRotFlag   % Flag whether auto rotation has been checked

% Windows
global sagwin obwin joywin navigwin

% Number of Sagittal slices
global numSlices

% Local Variables
volselpts = []; 		% Selected region in volume.
rvolselpts = [];
inpts = [];			    % list of alignment points in inplane
volpts = [];			% corresponding points in volume
obX = [0,0];			% Coordinates of sagittal and oblique slices
obY = [0,0];
obXM = [];			    % Coordinates of user set inplane slices
obYM = [];
lp = [];			    % pointers to the inplane lines we draw
ipThickness = -99;		% inplanes thickness (mm)
ipSkip = -99;			% amount of space skipped between inplane (mm)
curSag = -99;			% sagittal slice currently displaying
reflections = [1,1];	% Keeps track up left/right, up/down flips done by user
sagX = [0,0];
sagY = [0,0];
obMin = 0;
obMax = 0;
sagMin = 0;
sagMax = 500;
numSlices = 124;		% Number of planes of volume anatomy
obslice = [];
sagwin = [];			% Figure with sagittal view
obwin = [];			    % Figure with oblique view
sagSlice = [];			% Current sagittal image
sagPts = [];			% Samples for the current sagittal slice.
sagSize = [];			% Current sagittal size
sagCrop = [];			% Sagittal crop region
volume = [];			% Volume of data
obPts = [];			    % Locations in volume of oblique slice
obSlice = [];			% Current oblique image
obSize = [];			% Size of oblique image
obSizeOrig = [];        % real size of oblique, used for point selection

%%%%%%%%%%%% global constants and variable inits for rotation %%%%%%%%%%%

global axial coronal
axial = 1;
coronal = 2;                    % used to identify axis of rotation
cTheta = 0;                     % default rotation angles
aTheta = 0;
sagDelta = 0.03;               % angle delta in radians
sagDeltaMin = 0.003;
sagDeltaMax = 0.06;
transDelta = 10;               % trans delta in mm
transDeltaMin = 0;
transDeltaMax = 20;
aThetaSave = 0;
cThetaSave = 0;
curInplane = 0;

inpRotDelta = 20;           % inplane grid rotation increment (deg)
inpRotDeltaMin = .5;
inpRotDeltaMax = 35;

inOutDelta = 5;		     % in/out translation in mm
inOutDeltaMin = 0;
inOutDeltaMax = 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mrLoadRet stuff
% The goal here is to make mrAlign a standalone program so that you don't
% need to be running mrLoadRet. This stuff needs to be set up to allow
% mrAlign to run on its own.

global mrSESSION dataTYPES vANATOMYPATH
load mrSESSION
subject = mrSESSION.subject;
vAnatPath = getVAnatomyPath(subject)
numofanats = mrSESSION.inplanes.nSlices;
curSize = mrSESSION.inplanes.cropSize;
anatsz = [curSize,numofanats];
inplane_pix_size = 1./mrSESSION.inplanes.voxelSize;

% open inplane window
global INPLANE
openInplaneWindowForAlign3;
set(gcf, 'MenuBar', 'figure');
retwin = INPLANE.ui.figNum;		%mrLoadRet window
INPLANE = loadAnat(INPLANE);
INPLANE = refreshView(INPLANE);
set(gcf,'Units','Normalized','Position',[.5 .05 .45 .45]);

% open screen save window
if (~exist('Raw/Anatomy/SS','dir'))
    disp('No SS found');
    sswin=0;
else

openSSWindow;
sswin = gcf;
set(sswin,'Units','Normalized','Position',[.5 .05 .45 .45]);
end

try
    openSSWindow;
    sswin = gcf;
    set(sswin,'Units','Normalized','Position',[.5 .05 .45 .45]);
catch
    disp('Couldn''t find Screen Save, so no SS Window')
    sswin = [];
end

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

volcmap = gray(110);		% Set up color map
volcmap = [volcmap;hsv(110)];

% Sagittal window
sagwin = figure('MenuBar','none');
set(sagwin,'Units','Normalized','Position',[.05 .5 .45 .45]);
colormap(volcmap);

% Interpolated Oblique window
obwin = figure('MenuBar','none');
set(obwin,'Units','Normalized','Position',[.5 .5 .45 .45]);
colormap(volcmap);

% Navigation window
joywin = figure('MenuBar','none');
set(joywin,'Position', [0, 100, 650, 170]);
colormap(volcmap);

set(retwin, 'Name', 'Inplane');
if (sswin)
set(sswin, 'Name', 'Screen Save');
end

set(sagwin, 'Name', 'Interpolated Screen Save');
set(obwin, 'Name', 'Interpolated Inplane');
set(joywin, 'Name', 'Navigation Control');

% Make sagwin active.
figure(sagwin);

%%%%% Sagittal Buttons %%%%%

%  These are the arrows for moving from one sagittal to another
%  Rightward one step
volinc = uicontrol('Style','pushbutton','String','->','Units','normalized','Position',[.9,.95,.1,.05],'Callback','curSag=curSag+1;[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
% Leftward one step
voldec = uicontrol('Style','pushbutton','String','<-','Units','normalized','Position',[.8,.95,.1,.05],'Callback','curSag=curSag-1;[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
% Slider for moving across sagittals faster by dragging
volslislice = uicontrol('Style','slider','String','Pos','Units','normalized','Position',[.6,.95,.1,.05],'Callback','curSag=ceil(numSlices*get(volslislice,''value''));[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');

% Contrast controls
% Choose the max value set to highest display intensity
volslimax1 = uicontrol('Style','slider','String','Max','Units',...
	'normalized','Position',[.9,.6,.1,.05],'Callback',...
	'myShowImageVol(sagSlice,sagSize,max(sagSlice)*get(volslimin1,''value''),max(sagSlice)*get(volslimax1,''value''),obX,obY)');
% Choose the min value set to lowest display intensity
volslimin1 = uicontrol('Style','slider','String','Min','Units',...
	'normalized','Position',[.9,.4,.1,.05],'Callback',...
	'myShowImageVol(sagSlice,sagSize,max(sagSlice)*get(volslimin1,''value''),max(sagSlice)*get(volslimax1,''value''),obX,obY)');
set(volslimin1,'value',0);
set(volslimax1,'value',.50);

%%%%% Sagittal File Menu %%%%%
ld = uimenu('Label','File','separator','on');

% Reload volume anatomy data set.
uimenu(ld,'Label','Load Volume Anatomy','CallBack', ...
    '[volume,mmPerPix,volSize] = readVolAnat(vAnatPath); sagSize = volSize(1:2); numSlices = volSize(3); curSag = floor(numSlices/2); sagSlice = mrShowSagVol(volume,sagSize,curSag,[]);',...
    'Separator','off');
% Load AlignParams
uimenu(ld,'Label','(Re-)Load AlignParams','CallBack',...
    'load AlignParams;[obX,obY,obSize,obSizeOrig,sagPts,sagSlice,lp,obPts,obSlice]=mrReloadParams(lp,curInplane,obXM,obYM,sagSize,numSlices,volume,cTheta,aTheta,curSag,reflections,scaleFac);',...
    'Separator','on');
% Save AlignParams
uimenu(ld,'Label','Save AlignParams','CallBack',...
    'aThetaSave= aTheta; cThetaSave=cTheta;mrSaveAlignParams(obXM,obYM,subject,inplane_pix_size,ipThickness,ipSkip,curSag,curInplane,aTheta,cTheta)',...
    'Separator','on');
% Load alignment
uimenu(ld,'Label','Load Alignment (bestrotvol)','CallBack',...
    'eval(sprintf(''load bestrotvol''));disp(''Loaded: bestrotvol.mat'');',...
    'Separator','on');
% Save alignment
uimenu(ld,'Label','Save Alignment (bestrotvol)','CallBack',...
    'eval(sprintf(''save bestrotvol inpts volpts trans rot scaleFac''));disp(''Saved: bestrotvol.mat'');',...
    'Separator','on');
% Quit
uimenu(ld,'Label','Quit','CallBack', ...
    'close all; clear all;',...
    'Separator','on');

%%%%% Sagittal Inplanes Menu %%%%%
ld = uimenu('Label','Inplanes','separator','on');

% Create of the set of inplane lines initially
uimenu(ld,'Label','Setup/Refresh Inplanes','CallBack',...
	'[obX,obY,obXM,obYM,lp,ipThickness,ipSkip] =mrSetupInplanes(numofanats,obXM,obYM,lp,ipThickness,ipSkip,volume_pix_size,inplane_pix_size,curInplane);','Separator','on');
% Translate inplanes
uimenu(ld,'Label','Translate Inplanes','CallBack',...
	'[obX,obY,obXM,obYM,lp] = mrTransInplanes(numofanats,obXM,obYM,lp,curInplane);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);','Separator','on');
% Rotate inplanes
uimenu(ld,'Label','Clip Inplanes','CallBack', ...
 	'[obX,obY,obXM,obYM,lp] = mrClipInplanes(numofanats,obXM,obYM,lp,curInplane); [obSize,obSizeOrig] = mrFindObSize(obX,obY,sagSize,numSlices,curInplane);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);','Separator','on');
% Select inplane
uimenu(ld,'Label','Select Inplane','CallBack',...
    ['[obX,obY,lp,curInplane]=mrSelInplane(numofanats,obXM,obYM,lp,curInplane);[obSize,obSizeOrig] = mrFindObSize(obX,obY,sagSize,numSlices,curInplane);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);'...
    'INPLANE = viewSet(INPLANE, ''Current Slice'',curInplane);refreshView(INPLANE);'],'Accelerator','S','Separator','on');

%%%%% Sagittal Alignment Menu %%%%%

ld = uimenu('Label','Alignment','separator','on');

uimenu(ld,'Label','Compute Alignment from Many Points','CallBack', ...
 '[trans,rot]=mrDoAlignVol(inpts,volpts,scaleFac,curSize,sagSize,volume,numSlices,retwin,sagwin,obwin);','Separator','on');

uimenu(ld,'Label','Compute Alignment Automatically','CallBack', ...
   'NCoarseIter = 4; Usebestrotvol = 0; regEstRot;', 'separator', 'on');

uimenu(ld,'Label','Compute Alignment Automatically (use bestrotvol.mat as starting point)','CallBack', ...
   'NCoarseIter = 4; Usebestrotvol = 1; regEstRot;', 'separator', 'on');
% ON - in this case, if the alignment in bestrotvol is good, I guess that the
% coarse iterations can be eliminated. For now, and for safety, I leave
% it to be 4 (if bestrotvol is a previous manual alignment, then it will be
% good. But if it is a copy of an alignment computed for data acquired using
% the same slice prescription, it might be not good enough to eliminate
% the coarse iterations).

uimenu(ld,'Label','Check Alignment (current slice)','CallBack', ...
'chkImg = mrCheckAlignVol(rot,trans,scaleFac,curSize,viewGet(INPLANE, ''Current Slice''),volume,sagSize,numSlices,obwin);','Separator','on');

uimenu(ld,'Label','Check Alignment (all slices)',...
   'CallBack','regCheckRotAll', 'separator', 'on');

uimenu(ld,'Label','Check Alignment (overlay)',...
   'CallBack','OverlayRegCheck;', 'separator', 'on');

%%%%% Oblique Buttons %%%%%

figure(obwin);

%Select Points button
uicontrol('style','pushbutton','string','Select Points','units','normalized',...
	'Position',[0.0,.95,.2,.05],'CallBack',...
	'[inpts,volpts]=mrSelectPoints(inpts,volpts,retwin,obwin,viewGet(INPLANE, ''Current Slice''),obPts,obSizeOrig,reflections);');

%Undo Last Point button
uicontrol('style','pushbutton','string','Undo Last Point','units','normalized',...
	'Position',[.26,.95,.2,.05],'CallBack', ...
	'[inpts,volpts]=mrUndoAPoint(inpts,volpts);');

%Clear all points button
uicontrol('style','pushbutton','string','Clear All Points','units','normalized',...
	'Position',[.53,.95,.2,.05],'CallBack', ...
	'[inpts,volpts]=mrClearPoints(inpts,volpts);');

%Flip image left/right
uicontrol('style','pushbutton','string','Flip Right/Left','units','normalized',...
	'Position',[.26,.0,.2,.05],'CallBack', '[obSlice,reflections] = mrReflectObl(obSlice,obSize,reflections,1,1);');

%Flip image up/down
uicontrol('style','pushbutton','string','Flip Up/Down','units','normalized',...
	'Position',[.53,.0,.2,.05],'CallBack', '[obSlice,reflections] = mrReflectObl(obSlice,obSize,reflections,2,1);');

% Contrast controls
volslimax2 = uicontrol('Style','slider','String','Max','Units',...
		'normalized','Position',[.9,.6,.1,.05],'Callback',...
	'myShowImageVol(obSlice,obSize,max(obSlice)*get(volslimin2,''value''),max(obSlice)*get(volslimax2,''value''),sagX,sagY)');
volslimin2 = uicontrol('Style','slider','String','Min','Units',...
		'normalized','Position',[.9,.4,.1,.05],'Callback',...
	'myShowImageVol(obSlice,obSize,max(obSlice)*get(volslimin2,''value''),max(obSlice)*get(volslimax2,''value''),sagX,sagY)');
set(volslimin2,'value',.15);
set(volslimax2,'value',.90);

%%%%%%%%%%% Joystick Control Buttons %%%%%%%%%%%%

figure(joywin);

%%% Sagittal rotation %%%
uicontrol('Style','Text','Position',[42,140,90,14],'String','Sagittal rotation');
uicontrol('Style','pushbutton','String','<--','Position',[25,65,45,20],'CallBack', '[sagSlice,sagPts,cTheta,aTheta,obSlice,obPts]=mrRotSagVol2(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,axial,sagDelta,curSag,reflections,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','-->','Position',[100,65,45,20],'CallBack', '[sagSlice,sagPts,cTheta,aTheta,obSlice,obPts]=mrRotSagVol2(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,axial,-sagDelta,curSag,reflections,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','^','Position',[75,90,20,45],'CallBack', '[sagSlice,sagPts,cTheta,aTheta,obSlice,obPts]=mrRotSagVol2(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,coronal,sagDelta,curSag,reflections,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','v','Position',[75,15,20,45],'CallBack', '[sagSlice,sagPts,cTheta,aTheta,obSlice,obPts]=mrRotSagVol2(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,coronal,-sagDelta,curSag,reflections,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','0,0','Position',[75,65,20,20],'CallBack','aTheta=0;cTheta=0;[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
sagDeltaslider = uicontrol('Style','slider','Position',[120,100,80,20],...
			'Min',sagDeltaMin,'Max',sagDeltaMax,'Value',sagDelta,...
			'Callback','sagDelta=get(sagDeltaslider,''value'');');

%%% Translation %%%
uicontrol('Style','Text','Position',[230,140,60,14],'String','Translation');
uicontrol('Style','pushbutton','String','<--','Position',[200,65,45,20],'CallBack', '[obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,transDelta,1);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','-->','Position',[275,65,45,20],'CallBack', '[obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,transDelta,2);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','^','Position',[250,90,20,45],'CallBack','[obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,transDelta,3);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);' );
uicontrol('Style','pushbutton','String','v','Position',[250,15,20,45],'CallBack','[obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,transDelta,4);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
transDeltaslider = uicontrol('Style','slider','Position',[295,100,80,20],...
			'Min',transDeltaMin,'Max',transDeltaMax,'Value',transDelta,...
			'Callback','transDelta=get(transDeltaslider,''value'');');

%%% Rotation %%%
uicontrol('Style','Text','Position',[405,140,60,14],'String','Rotation');
uicontrol('Style','pushbutton','String','<--','Position',[375,65,45,20],'CallBack','[obXM,obYM]=mrRotInplanes(numofanats,obXM,obYM,(-1*inpRotDelta),curInplane);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','-->','Position',[450,65,45,20],'CallBack','[obXM,obYM] = mrRotInplanes(numofanats,obXM,obYM,inpRotDelta,curInplane);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
inpRotDeltaslider = uicontrol('Style','slider','Position',[415,100,80,20],...
			'Min',inpRotDeltaMin,'Max',inpRotDeltaMax,'Value', inpRotDelta,...
			'Callback','inpRotDelta=get(inpRotDeltaslider,''Value'');');

%%% In/out translation %%%
uicontrol('Style','Text','Position',[550,140,90,14],'String','In/Out');
uicontrol('Style','pushbutton','String','Out','Position',[525,65,45,20],'CallBack', '[obX,obY,obXM,obYM,lp] = mrPerpTransByButton(numofanats,obXM,obYM,lp,curInplane,inOutDelta,1);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
uicontrol('Style','pushbutton','String','In','Position',[600,65,45,20],'CallBack', '[obX,obY,obXM,obYM,lp] = mrPerpTransByButton(numofanats,obXM,obYM,lp,curInplane,inOutDelta,-1);[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);lp = mrRefreshInplanes(lp,obXM,obYM,curInplane,0);');
inOutDeltaslider = uicontrol('Style','slider','Position',[525,100,122,20],...
			'Min',inOutDeltaMin,'Max',inOutDeltaMax,'Value', inOutDelta,...
			'Callback','inOutDelta=get(inOutDeltaslider,''Value'');disp([''inOutDelta: '' num2str(inOutDelta)]);');

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

%%%% Matlab executes code from here to end, exits, then stays resident %%%%%%%

% Re-worked the logic of obtaining critical parameters -- 07.23.97 SPG,ABP
% AlignParams.mat doesn't exist

%Check if 'AlignParams' exist.
%If not, then create 'AlignParams'.
if ~check4File('AlignParams')
  mrSaveAlignParams(obXM,obYM,subject,...
    inplane_pix_size,...
    ipThickness,ipSkip,curSag,curInplane,aTheta,cTheta);
end
% Load AlignParams
mrLoadAlignParams;

%Load in the volume anatomy and display it sagwin
[volume,mmPerPix,volSize] = readVolAnat(vAnatPath);
sagSize = volSize(1:2);
numSlices = volSize(3);
%[volume, sagSize, numSlices, calc, dataRange] = mrLoadVAnatomy(voldr,subject);

% Get the volume pixel size here.
volume_pix_size = 1./mmPerPix;
%[volume_pix_size] = mrGetVolPixSize(voldr,subject);

%compile the scale factors for inplane and volume anatomies
scaleFac = [inplane_pix_size;volume_pix_size];

% First time through this will not be a parameter in VolParams.mat
if curSag < 0
	curSag = floor(numSlices/2);
end

figure(sagwin);
sagSlice = mrShowSagVol(volume,sagSize,curSag,[]);
%[sagSlice,sagPts,obSlice,obPts]=mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);

interpflag = 1; % global required to do 3dlinearinterp

% If the user has an inplane that is selected, display it.
if (curInplane ~=0)
	% draw the inplanes
	[obX,obY,obXM,obYM,lp,ipThickness,ipSkip] = mrSetupInplanes(numofanats,obXM,obYM,lp,ipThickness,ipSkip,volume_pix_size,inplane_pix_size,curInplane);
	% interpolated image size
	[obSize,obSizeOrig] = mrFindObSize(obX,obY,sagSize,numSlices,curInplane);
	% draw it
	[sagSlice,sagPts,obSlice,obPts] = mrRotSagVol(volume,obXM,obYM,obSize,sagSize,cTheta,aTheta,curSag,reflections,1,1/scaleFac(1,3),curInplane);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [rez, X] = splitAllClusters(rez, flag)
% i call this algorithm "bimodal pursuit"
% split clusters if they have bimodal projections
% the strategy is to maximize a bimodality score and find a single vector projection
% that maximizes it. If the distribution along that maximal projection crosses a
% bimodality threshold, then the cluster is split along that direction
% it only uses the PC features for each spike, stored in rez.cProjPC

ops = rez.ops;

wPCA = gather(ops.wPCA); % use PCA projections to reconstruct templates when we do splits

ccsplit = rez.ops.AUCsplit; % this is the threshold for splits, and is one of the main parameters users can change

NchanNear   = min(ops.Nchan, 32);
Nnearest    = min(ops.Nchan, 32);
sigmaMask   = ops.sigmaMask;

ik = 0;
Nfilt = size(rez.W,2);
nsplits= 0;

[iC, mask, C2C] = getClosestChannels(rez, sigmaMask, NchanNear); % determine what channels each template lives on

ops.nt0min = getOr(ops, 'nt0min', 20); % the waveforms must be aligned to this sample

[~, iW] = max(abs(rez.dWU(ops.nt0min, :, :)), [], 2); % find the peak abs channel for each template
iW = squeeze(int32(iW));

isplit = 1:Nfilt; % keep track of original cluster for each cluster. starts with all clusters being their own origin.
dt = 1/1000;
nccg = 0;

while ik<Nfilt
    if rem(ik, 100)==1
      % periodically write updates
       fprintf('Found %d splits, checked %d/%d clusters, nccg %d \n', nsplits, ik, Nfilt, nccg)
    end
    ik = ik+1;

    %
    isp = find(rez.st3(:,2)==ik); % get all spikes from this cluster
    nSpikes = numel(isp);
    if  nSpikes<300
       continue; % do not split if fewer than 300 spikes (we cannot estimate cross-correlograms accurately)
    end

    ss = rez.st3(isp,1)/ops.fs; % convert to seconds

    clp0 = rez.cProjPC(isp, :, :); % get the PC projections for these spikes
    clp0 = gpuArray(clp0(:,:));
    clp = clp0 - mean(clp0,1); % mean center them

    % now use two different ways to initialize the bimodal direction
    % the main script calls this function twice, and does both initializations
    if flag
        [u s v] = svdecon(clp');
        w = u(:,1); % initialize with the top PC
    else
        w = mean(clp0, 1)'; % initialize with the mean of NOT drift-corrected trace
        w = w/sum(w.^2)^.5; % unit-normalize
    end

    % initial projections of waveform PCs onto 1D vector
    x = gather(clp * w);
    s1 = var(x(x>mean(x))); % initialize estimates of variance for the first
    s2 = var(x(x<mean(x))); % and second gaussian in the mixture of 1D gaussians

    mu1 = mean(x(x>mean(x))); % initialize the means as well
    mu2 = mean(x(x<mean(x)));
    p  = mean(x>mean(x)); % and the probability that a spike is assigned to the first Gaussian

    logp = zeros(numel(isp), 2); % initialize matrix of log probabilities that each spike is assigned to the first or second cluster

    % do 50 pursuit iteration
    for k = 1:50
        % for each spike, estimate its probability to come from either Gaussian cluster
        logp(:,1) = -1/2*log(s1) - (x-mu1).^2/(2*s1) + log(p);
        logp(:,2) = -1/2*log(s2) - (x-mu2).^2/(2*s2) + log(1-p);

        lMax = max(logp,[],2);
        logp = logp - lMax; % subtract the max for floating point accuracy
        rs = exp(logp); % exponentiate the probabilities

        pval = log(sum(rs,2)) + lMax; % get the normalizer and add back the max
        logP(k) = mean(pval); % this is the cost function: we can monitor its increase

        rs = rs./sum(rs,2); % normalize so that probabilities sum to 1

        p = mean(rs(:,1)); % mean probability to be assigned to Gaussian 1
        mu1 = (rs(:,1)' * x )/sum(rs(:,1)); % new estimate of mean of cluster 1 (weighted by "responsibilities")
        mu2 = (rs(:,2)' * x )/sum(rs(:,2)); % new estimate of mean of cluster 2 (weighted by "responsibilities")

        s1 = (rs(:,1)' * (x-mu1).^2 )/sum(rs(:,1)); % new estimates of variances
        s2 = (rs(:,2)' * (x-mu2).^2 )/sum(rs(:,2));

        if (k>10 && rem(k,2)==1)
            % starting at iteration 10, we start re-estimating the pursuit direction
            % that is, given the Gaussian cluster assignments, and the mean and variances,
            % we re-estimate w
            StS  = clp' * (clp .* (rs(:,1)/s1 + rs(:,2)/s2))/nSpikes; % these equations follow from the model
            StMu = clp' * (rs(:,1)*mu1/s1 + rs(:,2)*mu2/s2)/nSpikes;

            w = StMu'/StS; % this is the new estimate of the best pursuit direection
            w = normc(w'); % which we unit normalize
            x = gather(clp * w);  % the new projections of the data onto this direction
        end
    end

    ilow = rs(:,1)>rs(:,2); % these spikes are assigned to cluster 1
    %    ps = mean(rs(:,1));
    plow = mean(rs(ilow,1)); % the mean probability of spikes assigned to cluster 1
    phigh = mean(rs(~ilow,2)); % same for cluster 2
    nremove = min(mean(ilow), mean(~ilow)); % the smallest cluster has this proportion of all spikes

    % did this split fix the autocorrelograms?
    [K, Qi, Q00, Q01, rir] = ccg(ss(ilow), ss(~ilow), 500, dt); % compute the cross-correlogram between spikes in the putative new clusters
    Q12 = min(Qi/max(Q00, Q01)); % refractoriness metric 1
    R = min(rir);                % refractoriness metric 2

    % if the CCG has a dip, don't do the split.
    % These thresholds are consistent with the ones from merges.
    if Q12<.25 && R<.05 % if both metrics are below threshold.
        nccg = nccg+1; % keep track of how many splits were voided by the CCG criterion
        continue;
    end

    % now decide if the split would result in waveforms that are too similar
    c1  = wPCA * reshape(mean(clp0(ilow,:),1), 3, []); %  the reconstructed mean waveforms for putatiev cluster 1
    c2  = wPCA * reshape(mean(clp0(~ilow,:),1), 3, []); %  the reconstructed mean waveforms for putative cluster 2
    cc = corrcoef(c1, c2); % correlation of mean waveforms
    n1 =sqrt(sum(c1(:).^2)); % the amplitude estimate 1
    n2 =sqrt(sum(c2(:).^2)); % the amplitude estimate 2

    r0 = 2*abs(n1 - n2)/(n1 + n2); % similarity of amplitudes

    % if the templates are correlated, and their amplitudes are similar, stop the split!!!
    if cc(1,2)>.9 && r0<.2
        continue;
    end

    % finaly criteria to continue with the split: if the split piece is more than 5% of all spikes,
    % if the split piece is more than 300 spikes, and if the confidences for assigning spikes to
    % both clusters exceeds a preset criterion ccsplit
    if nremove > .05 && min(plow,phigh)>ccsplit && min(sum(ilow), sum(~ilow))>300
       % one cluster stays, one goes
       Nfilt = Nfilt + 1;

       % the templates for the splits have been estimated from PC coefficients
       rez.dWU(:,iC(:, iW(ik)),Nfilt) = c2;
       rez.dWU(:,iC(:, iW(ik)),ik)    = c1;

       % the temporal components are therefore just the PC waveforms
       rez.W(:,Nfilt,:) = permute(wPCA, [1 3 2]);
       iW(Nfilt) = iW(ik); % copy the best channel from the original template
       isplit(Nfilt) = isplit(ik); % copy the provenance index to keep track of splits

       rez.st3(isp(ilow), 2)    = Nfilt; % overwrite spike indices with the new index
       rez.simScore(:, Nfilt)   = rez.simScore(:, ik); % copy similarity scores from the original
       rez.simScore(Nfilt, :)   = rez.simScore(ik, :); % copy similarity scores from the original
       rez.simScore(ik, Nfilt) = 1; % set the similarity with original to 1
       rez.simScore(Nfilt, ik) = 1; % set the similarity with original to 1

       rez.iNeigh(:, Nfilt)     = rez.iNeigh(:, ik); % copy neighbor template list from the original
       rez.iNeighPC(:, Nfilt)     = rez.iNeighPC(:, ik); % copy neighbor channel list from the original

       % try this cluster again
       ik = ik-1; % the cluster piece that stays at this index needs to be tested for splits again before proceeding
       % the piece that became a new cluster will be tested again when we get to the end of the list
       nsplits = nsplits + 1; % keep track of how many splits we did

    end
end

fprintf('Finished splitting. Found %d splits, checked %d/%d clusters, nccg %d \n', nsplits, ik, Nfilt, nccg)

Nfilt = size(rez.W,2); % new number of templates
Nrank = 3;
Nchan = ops.Nchan;
Params     = double([0 Nfilt 0 0 size(rez.W,1) Nnearest ...
    Nrank 0 0 Nchan NchanNear ops.nt0min 0]); % make a new Params to pass on parameters to CUDA

% we need to re-estimate the spatial profiles
[Ka, Kb] = getKernels(ops, 10, 1); % we get the time upsampling kernels again
[rez.W, rez.U, rez.mu] = mexSVDsmall2(Params, rez.dWU, rez.W, iC-1, iW-1, Ka, Kb); % we run SVD

[WtW, iList] = getMeWtW(single(rez.W), single(rez.U), Nnearest); % we re-compute similarity scores between templates
rez.iList = iList; % over-write the list of nearest templates

isplit = rez.simScore==1; % overwrite the similarity scores of clusters with same parent
rez.simScore = gather(max(WtW, [], 3));
rez.simScore(isplit) = 1; % 1 means they come from the same parent

rez.iNeigh   = gather(iList(:, 1:Nfilt)); % get the new neighbor templates
rez.iNeighPC    = gather(iC(:, iW(1:Nfilt))); % get the new neighbor channels

rez.isplit = isplit; % keep track of origins for each cluster

% figure(1)
% subplot(1,4,1)
% plot(logP(1:k))
%
% subplot(1,4,2)
% [~, isort] = sort(x);
% epval = exp(pval);
% epval = epval/sum(epval);
% plot(x(isort), epval(isort))
%
% subplot(1,4,3)
% ts = linspace(min(x), max(x), 200);
% xbin = hist(x, ts);
% xbin = xbin/sum(xbin);
%
% plot(ts, xbin)
%
% figure(2)
% plotmatrix(v(:,1:4), '.')
%
% drawnow
%
% % compute scores for splits
% ilow = rs(:,1)>rs(:,2);
% ps = mean(rs(:,1));
% [mean(rs(ilow,1)) mean(rs(~ilow,2)) max(ps, 1-ps) min(mean(ilow), mean(~ilow))]

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function ph = montage_clusters_points(ovl,clusters,XYZpts,varargin)
% ::
%
%    ph = montage_clusters_points(ovl,clusters,XYZpts,varargin)
%
% :Inputs:
%
%   **varargin:**
%        additional clusters structures
%
%   **XYZpts:**
%        XYZ mm coordinates of points to plot
%
%   **ph:**
%        point handles
%
% ..
%    Tor Wager
% ..

if isempty(ovl), ovl = which('scalped_single_subj_T1.img');, end
myc = {'b' 'r' 'y' 'g' 'm' 'c'};

if size(XYZpts,1) ~= 3, XYZpts = XYZpts';,end

XYZmm = cat(2,clusters.XYZmm);
XYZmm = [XYZmm XYZpts];

if length(varargin) > 0
    for i = 1:length(varargin)
        cl = varargin{i};
        clXYZmm = cat(2,cl.XYZmm);
        XYZmm = [XYZmm clXYZmm];
    end
end

% ----------------------------------------
% * overlay image
% ----------------------------------------
V = spm_vol(ovl);
oimg = spm_read_vols(V);
V.M = V.mat;

textx = size(oimg,2) - 50;
texty = size(oimg,1) - 6;

%[array,hdr,h,whichslices,rows,cols,figh] = readim2(ovl,'p');

% how many slices, which ones

XYZ = mm2voxel(XYZmm,V)';
whsl = unique(XYZ(3,:));
nsl = length(whsl) + 1;
rc = ceil(sqrt(nsl));
h = [];

figure; colormap gray; set(gcf,'Color','w')

% plot first cluster structure

index = plot_cluster(clusters,oimg,rc,V,whsl,myc{1},textx,texty,1,size(oimg));
ph = plot_points(XYZpts,rc,V,whsl,myc{2});

% plot additional cluster structures

if length(varargin) > 0
    for i = 1:length(varargin)
        cl = varargin{i};
        index = plot_cluster(cl,[],rc,V,whsl,myc{i+1},textx,texty,i+1,size(oimg));

    end
end

return

% sub-functions
%

function index = plot_cluster(clusters,oimg,rc,V,whsl,myc,textx,texty,clind,odims)

XYZmm = cat(2,clusters.XYZmm);
XYZ = mm2voxel(XYZmm,V)';

% surface patch method
% ----------------------------------------------------------------------------------------
vol = voxel2mask(XYZ',odims);
vol = smooth3(vol);

index = 1;
for z = whsl

    subplot(rc,rc,index);

    if ~isempty(oimg)
        set(gca,'YDir','reverse');
        imagesc(oimg(:,:,z))
        hold on; axis image; axis off
        zmm = voxel2mm([1 1 z]',V.mat);
        text(textx,texty,['z = ' num2str(zmm(3))],'Color','w')
    else
        hold on
    end

    if z>1,FVC = isocaps(vol(:,:,z-1:z),0,'zmax');
    else FVC = isocaps(vol(:,:,z:z+1),0,'zmax');
    end
    try
	patch(FVC,'EdgeColor','none','FaceColor',myc,'FaceAlpha',.7)
    catch
	patch(FVC,'EdgeColor','none','FaceColor',myc)
    end

    % plot method

    %myXYZ = XYZ(1:2,XYZ(3,:) == z);
    %plot(myXYZ(2,:),myXYZ(1,:),[myc 's'],'MarkerFaceColor',myc,'MarkerSize',3)

    index = index + 1;
    drawnow

end

subplot(rc,rc,index)
a = pwd; a = a(end-6:end);
b = num2str(clusters(1).threshold);

c = num2str(length(clusters));
text(0,clind-1,[myc ': ' clusters(1).title ' ' a ' u = ' b ', ' c ' clusters'])
axis off
axis([0 1 -1 clind])

return

function ph = plot_points(XYZmm,rc,V,whsl,myc)

XYZ = mm2voxel(XYZmm,V)';
index = 1;
phind = 1;
for z = whsl

    subplot(rc,rc,index);
    hold on
    myXYZ = XYZ(:,XYZ(3,:) == z);

    for i = 1:size(myXYZ,2)
        ph(phind) = plot3(myXYZ(2,i),myXYZ(1,i),myXYZ(3,i),[myc(1) '.'],'MarkerFaceColor',myc(1),'MarkerSize',8);
        phind = phind + 1;
    end

    index = index + 1;

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [mMLS] = calc_FitComp(mCatalog, fMinRad, fMaxRad, fRadIncr, fX, fY)
% function [mMLS] = calc_FitComp(mCatalog, fMinRad, fMaxRad, fRadIncr, fX, fY);
% -----------------------------------------------------------------------------
%
% Function to calculate the MLS fit and the Goodness of fit rate to discriminate which
% method works best exploring increasing radii
%
% Incoming variables:
%
% Author: J. Woessner
% woessner@seismo.ifg.ethz.ch
% last update: 23.01.03

% Initialize
mMlS = [];

% Start value
fRadius = fMinRad;
while fRadius <= fMaxRad
    % Create catalog
    vDistances_ = sqrt(((mCatalog(:,1)-fX)*cos(pi/180*fY)*111).^2 + ((mCatalog(:,2)-fY)*111).^2);
    % Select those in between the maximum radius
    vSel = (vDistances_ <= fRadius);
    vCheckDist = vDistances_(vSel, :);
    mRadCatalog = mCatalog(vSel, :);
    % Determine Mc-values
    [result]=sv_NodeCalcMc(mRadCatalog); % For parameters in result see sv_NodeCalcMc

    % Calculate normal and lognormal fit
    [mResult, fProbNorm, fMcNorm, vX_resNorm, fNmaxNorm, mDatPredNorm] = calc_McCdfnormal(mRadCatalog, 0.1);
    [mResult2, fProbLog, fMcLog, fMuLog, fSigmaLog, mDatPredLog, vPredBest] = calc_McCdflognormal(mCat, fBinning);
    [fProbExp, fMcExp, vX_resExp, fNmaxExp, mDatPredExp] = calc_McCdfexp(mRadCatalog, 0.1);

    % Show data fit below Mc
    vSel = (mDatPredNorm(:,2) < fMcNorm);
    mTmpNorm = mDatPredNorm(vSel,:);
    vSel = (mDatPredLog(:,2) < fMcLog);
    mTmpLog = mDatPredLog(vSel,:);
    vSel = (mDatPredExp(:,2) < fMcExp);
    mTmpExp = mDatPredExp(vSel,:);
    figure;
    plot(mTmpNorm(:,2), mTmpNorm(:,3),'k*',mTmpNorm(:,2), mTmpNorm(:,1),'-b');
    hold on;
    plot(mTmpLog(:,2), mTmpLog(:,1),'-r');
    plot(mTmpExp(:,2), mTmpExp(:,1),'-g');
    legend('Data','Normal CDF', 'Lognorm. CDF', 'Exponential func.');
    sTitlestr = ['McNorm = ' num2str(fMcNorm) ', McLog = ' num2str(fMcLog) ', McExp = ' num2str(fMcExp) ', R = ' num2str(fRadius)];
    title(sTitlestr)
    xlabel('Magnitude')
    drawnow;
    hold off;
    sPrintstr = ['Fit_kobe_1986_92_135.6_35' num2str(fRadius) 'km.eps'];
    print('-deps2c', '-tiff','-r400', sPrintstr);

    % Result array
    mMLS = [mMLS; fRadius result.fProbMcNorm result.fProbMcLog fProbExp result.fMc_max result.fMc_90 result.fMc_95...
            result.fMc_com result.fMcNorm result.fMcLog fMcExp];

    % Plot Non-cumulative distribution, original and predicted
    % Time period
    [vFMD, vNonCFMD] = calc_FMD(mRadCatalog);
    vNonCFMD = fliplr(vNonCFMD);
    fPeriod1 = max(mRadCatalog(:,3)) - min(mRadCatalog(:,3));
    %     figure_w_normalized_uicontrolunits('tag','ncumdist','Name','Best model','Units','normalized','Nextplot','add',...
    %         'Numbertitle','off','visible','on');
    figure;
    semilogy(vNonCFMD(1,:)', vNonCFMD(2,:)', '-k^',mDatPredNorm(:,2) ,mDatPredNorm(:,1).*fPeriod1,'-bo');
    hold on;
    semilogy(mDatPredLog(:,2) ,mDatPredLog(:,1).*fPeriod1,'-r*');
    semilogy(mDatPredExp(:,2) ,mDatPredExp(:,1).*fPeriod1,'-gs');
    legend('Data','Normal CDF', 'Lognorm. CDF', 'Exponential func.');
    sTitlestr = ['McNorm = ' num2str(fMcNorm) ', McLog = ' num2str(fMcLog) ', McExp = ' num2str(fMcExp) ', R = ' num2str(fRadius)];
    title(sTitlestr)
    xlabel('Magnitude')
    ylabel('Non-cumulative FMD')
    drawnow;
    sPrintstr = ['FMD_kobe_1986_92_135.6_35' num2str(fRadius) 'km.eps'];
    print('-deps2c', '-tiff','-r400', sPrintstr);
    % Increase radius
    fRadius = fRadius+fRadIncr;
end

figure;
subplot(2,1,1);
plot(mMLS(:,1), mMLS(:,8),'-bd', mMLS(:,1), mMLS(:,9),'-rd', mMLS(:,1), mMLS(:,10),'-gd');
xlabel('Radius / [km]')
ylabel('Mc')
legend('Normal CDF', 'Lognorm. CDF', 'Exponential func.');
subplot(2,1,2);
plot(mMLS(:,1), mMLS(:,2),'-bd', mMLS(:,1), mMLS(:,3),'-rd', mMLS(:,1), mMLS(:,4),'-gd');
xlabel('Radius / [km]')
ylabel('MLS')
drawnow;
sPrintstr = ['Radius_kobe_1986_92_135.6_35.eps'];
print('-deps2c', '-tiff','-r400', sPrintstr);

author: CelsoReyes


text:
function [alpha_fin,b_fin] = rsptrain(K,Y,C,T,S_min),
% Recursive Stabilization Procedure
%% Variable Description
%%
%%
m = size(K,1);
k = ceil(log(m/S_min)/log(T/(T-1)));
current_position = [0,0];
Node_done = zeros(k,T);
Current_Path = zeros(k*T,m);
not_finished = 1;
count = 0;
Path_from_root = zeros(k,1);
while not_finished,

    cp1 = current_position(1);
    cp2 = current_position(2);

    if (cp1~=k),
        active = find(Node_done(cp1+1,1:T)==0);
    else
        active = [];
    end;%if (cp1~=k)

    if ~isempty(active),

        % The part of the execution tree is not completely explored
        if cp1==0,
            tmp=[1:m];
        else
            tmp = find(Current_Path(cp2 + (cp1-1)*T,:)==1);
        end; % if cp1==0,

        % Build the training set for the node
        tab_tmp = zeros(1,m);
        train_ind_tmp = tmp([1:floor((active(1)-1)*length(tmp)/T),floor(active(1)*length(tmp)/T):length(tmp)]);
        %% If all the points are from the same class then do not go further
        is_plus = find(Y(train_ind_tmp)==1);

        if (length(is_plus)==0|(length(is_plus)==length(train_ind_tmp))),
            Node_done(cp1+1,1:T)=1;
            Current_Path([1:T] + (cp1)*T,:) = zeros(T,m);
        else
            tab_tmp(train_ind_tmp)=1;
            current_position(1) = cp1+1;
            Path_from_root(current_position(1)) = current_position(2);
            current_position(2) = active(1);
            Current_Path(current_position(2) + (current_position(1)-1)*T,:) = tab_tmp;
        end;
    else
        if (cp1==0),
            tmp = [1:m];
        else
            tmp = find(Current_Path(current_position(2)+(current_position(1)-1)*T,:)==1);
        end;
            Ylearn = Y(tmp,:);
            Clearn = (C*m)/(length(tmp));
            Klearn = K(tmp,tmp);

        if current_position(1)==k,
            [alpha] = quadsolve(Klearn,-ones(length(tmp),1),Ylearn',0,Clearn);
            alpha_bias=alpha;
            scale_f = sqrt(alpha'*Klearn*alpha);
            alpha = alpha/scale_f;
        else
            c = Klearn*(Current_Path(current_position(1)*T+1:current_position(1)*T+T,tmp))';
            c = mean(c,2);
            c = -1+c;
            H = 1/(T)*Klearn;
            [alpha,y] = quadsolve(H,c,Ylearn',0,Clearn);
            alpha_bias = alpha;
            alpha = 1/T*(alpha + sum(Current_Path(current_position(1)*T+1:current_position(1)*T+T,tmp)',2));
            scale_f = sqrt(alpha'*Klearn*alpha);
            alpha = alpha/scale_f;
        end; % if current_pos...==k-1

        if current_position ~=0,
            Current_Path(current_position(2)+(current_position(1)-1)*T,tmp) = alpha';
            Node_done(current_position(1)+1,1:T)=0;

            Node_done(current_position(1),current_position(2)) = 1;
            count=count+1;
            current_position(2) = Path_from_root(current_position(1));
            current_position(1) = current_position(1)-1;
        else,
            not_finished=0;
            alpha_fin = alpha;
            b_fin = -y/scale_f;
        end; % if current_pos...~=0

    end;% if ~isempty(active)

%    disp(sprintf('Node done: %d\n',count));

end; %while not_finished
% End function


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% Discrete Fourier Transform  properties

% symmetry

x=[1 2 3 4 5];
Xk=dft(x)

k=0:4
N=5
R=Xk(1+mod(N-k,N))
Xnk=conj(R)


text:
function C = EstimateUnitNormalCorr(Y_X, exp_window_size, shrink_factor)
%C = EstimateUnitNormalCorr(Y_X, exp_window_size, shrink_factor)
%  Estimate the correlation matrix of a joint normal sample, where each
%  element is assumed to be unit normal. Each column of Y_F is a random
%  variable and each row is a sample. Exponential weighting is applied in
%  time, and RMT filtering is applied using the "PG" algorithm.
%
% To do!!:
%   - review implementation of exponential weighting and RMT
%   - add choice for rectangular averaging?
%   - add choice for other RMT schemes?
% Code by S. Gollamudi. This version December 2009.

[T,N] = size(Y_X);
use_rmt_cleaning = 0;

%% Estimate sample correlation matrix with exponential averaging

% compute exponential weights
lambda = 2/exp_window_size;
exp_wts = (1-lambda).^(T-1:-1:0)*(lambda/(1-(1-lambda)^T));
m_exp_wts = repmat(exp_wts',1,N);

% compute the exponentially weighted sample mean
meanYX = sum(m_exp_wts.*Y_X, 1);
% remove mean from the time series Y_X
Y_X = Y_X - repmat(meanYX,T,1);

% compute the exponentially weighted standard deviation
stdYX = sqrt(sum(m_exp_wts.*(Y_X.^2), 1));

% compute the exponentially weighted correlation matrix
C = ((m_exp_wts.*Y_X)'*Y_X)./(stdYX'*stdYX);
diag_indx = 1:(N+1):(N*N);

%% Clean the correlation matrix estimate using Random Matrix Theory

if use_rmt_cleaning,

    % ratio of the effective series length and the number of variables
    Q = exp_window_size/N;

    % maximum random eigenvalue
    e_max = rmtEigenLim(Q,1);

    [eigVec,eigVal] = eig(C);
    [eigVal,Idx] = sort(diag(eigVal),'descend');
    eigVec=eigVec(:,Idx);
    % K is the # of non-random eigenvalues
    K=length(find(eigVal>e_max));
    % PG filter
    eigVal(K+1:end)=mean(eigVal(K+1:end));
    C = eigVec*diag(eigVal)*eigVec';
    % make the diagonal elements equal to one
    C(diag_indx) = 1;

end

%% Shrink to average correlation

% compute the average pairwise correlation coefficient
rho = (sum(sum(C)) - N)/(N*(N-1));
C_rho = rho*ones(N,N);
C_rho(diag_indx) = 1;

% shrink C towards C_rho
if nargin < 3, shrink_factor = 1/3; end
C = (1-shrink_factor)*C + shrink_factor*C_rho;


text:
% poolType = 'mean' or 'max'
% poolLayer = integer layer number after which we do the pooling
function layer = genNetworkFeedForward_pool(inputDim, hiddenLayerSize, outputDim, costFn, poolType, poolAfterNlayer, LastActivation4MSE)
if nargin<7
    LastActivation4MSE = 'linear';
end

layer = genNetworkFeedForward_v2(inputDim, hiddenLayerSize, outputDim, costFn, LastActivation4MSE);

poolLayer.name = poolType;
poolLayer.dim = [1 1]*layer{poolAfterNlayer}.dim(1);
poolLayer.prev = layer{poolAfterNlayer}.prev;
layer = [layer(1:poolAfterNlayer) poolLayer layer(poolAfterNlayer+1:end)];

% automatically derive the list of layers that the output of the current layer goes.
for i=1:length(layer); layer{i}.next = []; end
for i=length(layer):-1:1
    if isfield(layer{i}, 'prev')
        for j=1:length(layer{i}.prev)
            layer{i+layer{i}.prev(j)}.next(end+1) = -layer{i}.prev(j);
        end
    end
end

end

author: singaxiong

repo: SignalGraph


text:
function [x y pos trials map mapping] = getBehavEvents_linearTrack(pos,bins)
dbstop if error
scatter(pos(:,8),pos(:,10),'.')
% axis([-1 1 -1 1])
[x y]=ginput();
cx =  x(end-1:end);
y = y(end-1:end);
dist(:,1)=fastrms(abs(pos(:,8)-x(1))+abs(pos(:,10)-y(1)),120);
dist(:,2)=fastrms(abs(pos(:,8)-x(2))+abs(pos(:,10)-y(2)),120);

plot(dist,'.')
% axis([0 length(dist) 0 3])
[xx yy] = ginput();

[pks_one locs_one]=findpeaks(-dist(:,1),'MINPEAKHEIGHT',-yy);
[pks_two locs_two]=findpeaks(-dist(:,2),'MINPEAKHEIGHT',-yy);
c=1;
lastTrial = 0;
for i=1:length(locs_one)
    f = locs_one(i)-locs_two;
    ff = find(f>0);
    [a b]=min(locs_one(i)-locs_two(ff));
    if a < 500 & a > 60 & locs_two(ff(b)) > lastTrial + 60
       trials{c} = pos(locs_two(ff(b)):locs_one(i),:);
       c=1+c;
       lastTrial = locs_two(ff(b));
    end
    f = locs_two-locs_one(i);
    ff = find(f>0);
    [a b]=min(locs_two(ff)-locs_one(i));
    if a < 500 & a > 60 & locs_one(i) > lastTrial + 60
       trials{c} = pos(locs_one(i):locs_two(ff(b)),:);
       c=1+c;
       lastTrial = locs_one(i);
    end
end

% for i=1:length(locs_one)
%     two = locs_two(locs_two-locs_one(i)>0);
%     if ~isempty(two)
%     [a b]=min(abs(two-locs_one(i)));
%     if two(b)-locs_one(i) < 24000 && two(b)-locs_one(i) > 120 && locs_one(i) > lastTrial-60
%         trials{c} = pos(locs_one(i):two(b),:);
%         c=1+c;
%         lastTrial = two(b);
%     end
%     end
%     two = locs_two(locs_one(i)-locs_two>0);
%     if ~isempty(two)
%     [a b]=min(abs(two-locs_one(i)));
%     if locs_one(i)-two(b) < 24000 && locs_one(i)-two(b) > 120 && two(b) > lastTrial-60
%         trials{c} = pos(two(b):locs_one(i),:);
%         c=1+c;
%         lastTrial = locs_one(i);
%     end
%     end
% end

%% fine tune trials here
% for i=1:length(trials)
% clear t tt ttt tttt; t = trials{i};
% tt=t(:,8:10);
% d = diff(trials{i}(:,8:10));
% tt(find(d==0)+1)=nan;
% for j=1:3
% tttt(:,j) = smooth(fillmissing(tt(:,j),'pchip'),15);
% end
% subplot(3,1,1)
% plot((tttt))
% axis([0 length(tttt) -1 1])
% % hold on
% % plot(tttt(:,2));
% % hold off
% subplot(3,1,2)
% plot(diff(tttt))
% hold on
% plot(sum(abs(diff(tttt(:,[1 3])))'))
% hold off
% axis([0 length(tttt) -.015 .015])
% subplot(3,1,3)
% plot(abs(tttt(:,1)-tttt(:,3)))
% axis([0 length(tttt) 0 1.4])
% % scatter3(tttt(:,1),tttt(:,3),tttt(:,2),'.k')
% [xxx yyy]=ginput();
% if xxx(2) > length(trials{i})
%     xxx(2)=length(trials{i});
% end
% trials_new{i}=trials{i}(round(xxx(1)):round(xxx(2)),:)
% end

%% merge trials to the same length
trials_unsorted = trials;
for i=1:length(trials_unsorted)
       d = pdist(trials_unsorted{i}(:,[8 9 10]),'euclidean');
       dd = squareform(d);
       dd(dd==0)=nan;
       [mins ind] = min(dd(:));
           [row col] = ind2sub(size(dd),ind);
%            f = find(abs(row-col)==1);  % this worked but it was sloowww,
%            row = row(f);               % removed and replaced w ind2sub
%            col = col(f);
    while length(trials_unsorted{i})>bins
         [mins ind] = min(dd(:));
         [row col] = ind2sub(size(dd),ind);
       % check that row and col are adjacent
       if ((row(1)-col(1)))== -1
           trials_unsorted{i} = [trials_unsorted{i}(1:row(1)-1,:);...
                                nanmean(trials_unsorted{i}(row(1):col(1),:));...
                                trials_unsorted{i}(col(1)+1:end,:)];
           dd(row,:)=[];
           dd(:,col)=[];
       elseif ((row(1)-col(1)))== 1
           trials_unsorted{i} = [trials_unsorted{i}(1:col(1)-1,:);...
                                nanmean(trials_unsorted{i}(col(1):row(1),:));...
                                trials_unsorted{i}(row(1)+1:end,:)];
           dd(row,:)=[];
           dd(:,col)=[];
       elseif abs((row(1)-col(1)))>1 | abs((row(1)-col(1)))<1
           dd(row,col)=nan;
       end
    end
    i
end

% cluster similarity
for i=1:length(trials_unsorted)
    for j =1:length(trials_unsorted)
        for ax = 8:10
        temp(ax-7,:)=crosscorr(trials_unsorted{i}(:,[ax]),trials_unsorted{j}(:,[ax]),bins-1);
        end
        if sum(~isreal(temp(:)))>0
            disp([i j]);
        end
        a = trials_unsorted{i}(:,[8 9 10]);
        b = trials_unsorted{j}(:,[8 9 10]);
% a = [trials_unsorted{i}(:,[9]) trials_unsorted{i}(:,[10])-trials_unsorted{i}(:,[9])];
% b = [trials_unsorted{j}(:,[9]) trials_unsorted{j}(:,[10])-trials_unsorted{j}(:,[9])];
%         temp = corrcoef(reshape(a,bins*3,1),reshape(b,bins*2,1));
%           temp = crosscorr(reshape(a,bins*3,1),reshape(b,bins*,1),120);
%         cc(i,j)=temp(2);
          [cc(i,j) locs]=max(nanmean(temp));
    end
end
% remove NaN's?
n = isnan(cc(:,1));

%%
numConditions = str2num(input('How many conditions were there: ','s'));
while ~isnumeric(numConditions) || isempty(numConditions)
    numConditions = str2num(input('How many conditions were there: ','s'));
end
exitC=0;
while exitC == 0
    clust = kmeans(cc,numConditions);
    clear trials
    for i=1:numConditions
        f = find(clust==i);
        for j=1:length(f)
            trials{i}{j}=trials_unsorted{f(j)};
        end
    end
    for i=1:numConditions
        subplot(4,5,i)
        hold off
        for t = 1:length(trials{i})
            scatter3(trials{i}{t}(:,8),trials{i}{t}(:,10),trials{i}{t}(:,9),'.k')
            hold on
            scatter3(trials{i}{t}(1,8),trials{i}{t}(1,10),trials{i}{t}(1,9),'.g')
            scatter3(trials{i}{t}(end,8),trials{i}{t}(end,10),trials{i}{t}(end,9),'.r')
        end
    end
    pause(.01);
    adjust = input('merge[m]/split[s] ','s');
    if strcmp(adjust,'s')
        numConditions = numConditions+1;
    end
    if strcmp(adjust,'m')
        on = str2num(input('condition #1:','s'));
        tw = str2num(input('condition #2','s'));
        for t = 1:length(trials{tw})
           trials{on}{length(trials{on})+1}=trials{tw}{t};
        end
        trials{tw}=[];
    end
    if ~strcmp(adjust,'s') && ~strcmp(adjust,'m')
       exitC = 1;
       clf();
    end
end
% check that clustering was correct

c=1;
for tt = 1:length(trials)
%             if d == 5  %% only for rec 21st
%                 bins = 100;
%             elseif d == 4
%                 bins = 100;
%             else
%                 bins = 80;
%             end

    hold off
    map{tt}=[];
    t_conc=[];
    for t = 1:length(trials{tt})
        t_conc = [trials{tt}{t}(:,:),20*(trials{tt}{t}(:,1)-trials{tt}{t}(1,1))];
    if length(t_conc)>=bins
        while length(t_conc)>bins+1
        di = pdist(t_conc);
        s = squareform(di);
        s(find(eye(size(s))))=nan;
        [a b] = min(s(:));

        [coords blah] = find(s==a);
        t_conc(coords(1),:) = (t_conc(coords(1),:)+t_conc(coords(2),:))./2;
        t_conc(coords(2),:) = [];
        % debug
%         scatter(t_conc(:,1),t_conc(:,2))
%         pause(.01)
        end
    t_conc_all(t,:,:) = t_conc;
    end
    end
    if length(trials{tt})>0
    map{tt} = squeeze(median(t_conc_all(:,:,:),1));

    if ~isempty(trials{tt})
    subplot(8,4,c)
    c=1+c
    for t = 1:length(trials{tt})
    scatter(trials{tt}{t}(:,1),trials{tt}{t}(:,2),'.k')
    hold on
    scatter(trials{tt}{t}(1,1),trials{tt}{t}(1,2),'.g')
    scatter(trials{tt}{t}(end,1),trials{tt}{t}(end,2),'.r')
    end
    end
    end
    clear t_conc_all
    for t =1:length(trials{tt})  % all trial types (rotations)
        for p = 1:length(trials{tt}{t})
            [a b] = min(nansum(abs([trials{tt}{t}(p,1)-map{tt}(:,1),...
                trials{tt}{t}(p,8)-map{tt}(:,8),...
                trials{tt}{t}(p,9)-map{tt}(:,9),...
                trials{tt}{t}(p,10)-map{tt}(:,10),...
                (trials{tt}{t}(p,1)-trials{tt}{t}(1,1))*50-map{tt}(:,1),...  % penalty for time differences
                40*(p./length(trials{tt}{t})*length(map{tt}) - (1:length(map{tt})))'])'));     % penalty for order differences
            mapping{tt}{t}(p,:) = [map{tt}(b,1:end) b trials{tt}{t}(p,1)];
%             plot(nansum(abs([trials{tt}{t}(p,1)-map{tt}(:,1),...
%                 trials{tt}{t}(p,2)-map{tt}(:,2),...
%                 trials{tt}{t}(p,3)-map{tt}(:,3),...
%                 trials{tt}{t}(p,4)-map{tt}(:,4),...
%                 (trials{tt}{t}(p,5)-trials{tt}{t}(1,5))*20-map{tt}(:,5),...  % penalty for time differences
%                 50*(p./length(trials{tt}{t})*length(map{tt}) - (1:length(map{tt})))'])'))
%             pause
        end
    end
end

%% incorporate this below to reformat things and double check trial assignments...
% load behav
%
c=1
clear m mm tt
for i=1:length(trials)
    if ~isempty(trials{i})
    tt{c}=trials{i};
    m{c}=map{i};
    mm{c}=mapping{i};

    tt{c}=tt{c}(~cellfun('isempty',tt{c}));
    mm{c}=mm{c}(~cellfun('isempty',mm{c}));
    c=1+c;
    end
end
dbmap=m;
mapping=mm;
trials= tt;

for tt = 1:length(trials)
subplot(5,4,tt)
for t = 1:length(trials{tt})
%     scatter(map{tt}(:,1),map{tt}(:,2),'.')
scatter(trials{tt}{t}(:,8),trials{tt}{t}(:,10),'.k')
hold on

%     axis([0 550 0 550])
end
for t = 1:length(trials{tt})
%     scatter(map{tt}(:,1),map{tt}(:,2),'.')
scatter(trials{tt}{t}(1,8),trials{tt}{t}(1,10),'.g')
scatter(trials{tt}{t}(end,8),trials{tt}{t}(end,10),'.r')
hold on

%     axis([0 550 0 550])
end
title(tt);
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% DEMCMU35GPLVM3 Learn a GPLVM on CMU 35 data set.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

% Get the sequence numbers.
[Y, lbls] = lvmLoadData('cmu35WalkJog');
seq = cumsum(sum(lbls)) - [1:31];

dataSetName = 'cmu35gplvm';
experimentNo = 3;

% load data
[Y, lbls, Ytest, lblstest] = lvmLoadData(dataSetName);
origBias = mean(Y, 1);
origScale = 1./sqrt(var(Y, 1));
Y = Y - repmat(origBias, size(Y, 1), 1);
Ytest = Ytest - repmat(origBias, size(Ytest, 1), 1);
Y = Y.*repmat(origScale, size(Y, 1), 1);
Ytest = Ytest.*repmat(origScale, size(Ytest, 1), 1);

% Set up model
options = fgplvmOptions('fitc');
options.optimiser = 'conjgrad';
options.back = 'mlp';
options.backOptions = mlpOptions(10);
options.numActive = 100;
options.fixInducing = 1;
options.fixIndices = round(linspace(1, size(Y, 1), options.numActive));
latentDim = 5;

d = size(Y, 2);
model = fgplvmCreate(latentDim, d, Y, options);

% Add dynamics model.
optionsDyn = gpOptions('fitc');
optionsDyn.numActive = 100;
optionsDyn.fixInducing = 1;
optionsDyn.kern = kernCreate(model.X, {'rbf', 'white'});
optionsDyn.kern.comp{1}.inverseWidth = 0.2;
% This gives signal to noise of 0.1:5e-3 or 20:1.
optionsDyn.kern.comp{1}.variance = 0.01;
optionsDyn.kern.comp{2}.variance = 0.95;
diff = 1;
learn = 1;
optionsDyn.fixIndices = round(linspace(1, size(Y, 1)-length(seq), options.numActive));
model = fgplvmAddDynamics(model, 'gp', optionsDyn, diff, learn, seq);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model', 'origScale', 'origBias');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [X, objV] = proxF_tSVD_1(Y,rho,opts)
% proximal function of tensor nuclear norm
% Authors: G. Ely, S. Aeron, Z. Zhang, ECE, Tufts Univ. 03/16/2015

%% Input checking
parOP = false;
if ~exist('opts','var') && ~isempty(opts)
    nARG = length(opts);
    if nARG > 0
        parOP=opts{1};
    end
end

%%
sa = size(Y);
la = length(sa);
[n1,n2,n3] = size(Y);

% this returned this STILL FFT'd version
[U,S,V]=ntsvd(Y,1,parOP);

sTrueV =zeros(min(n1,n2),n3);
for i = 1:n3
    s = S(:,:,i);  % tube
    s = diag(s);
    sTrueV(:,i) = s;
end
%%
[sTrueV] = proxF_l1(sTrueV,rho);
%objV = sum(sTrueV(:));
objV = sum(sTrueV(:));
%%
for i = 1:min(n1,n2)
    for j = 1:n3
        S(i,i,j) = sTrueV(i,j);
    end
end

for i = la:-1:3
    U = ifft(U,[],i);
    S = ifft(S,[],i);
    V = ifft(V,[],i);
end

% X = tprod( tprod(U,S), ttrans_HO(V));

X = tprod( tprod(U,S), tran(V));

author: andrewssobral


text:
function [ rank, x ] = tuple_next ( m1, m2, n, rank, x )

%*****************************************************************************80
%
%% TUPLE_NEXT computes the next element of a tuple space.
%
%  Discussion:
%
%    The elements are N vectors.  Each entry is constrained to lie
%    between M1 and M2.  The elements are produced one at a time.
%    The first element is
%      (M1,M1,...,M1),
%    the second element is
%      (M1,M1,...,M1+1),
%    and the last element is
%      (M2,M2,...,M2)
%    Intermediate elements are produced in lexicographic order.
%
%  Example:
%
%    N = 2, M1 = 1, M2 = 3
%
%    INPUT        OUTPUT
%    -------      -------
%    Rank  X      Rank   X
%    ----  ---    -----  ---
%    0     * *    1      1 1
%    1     1 1    2      1 2
%    2     1 2    3      1 3
%    3     1 3    4      2 1
%    4     2 1    5      2 2
%    5     2 2    6      2 3
%    6     2 3    7      3 1
%    7     3 1    8      3 2
%    8     3 2    9      3 3
%    9     3 3    0      0 0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M1, M2, the minimum and maximum entries.
%
%    Input, integer N, the number of components.
%
%    Input, integer RANK, counts the elements.
%    On first call, set RANK to 0.  On subsequent calls, the input value of
%    RANK should be the output value of RANK from the previous call.
%
%    Input, integer X(N), the previous tuple.
%
%    Output, integer RANK, the order of the next tuple.  When there are no
%    more elements, RANK will be returned as 0.
%
%    Output, integer X(N,1), the next tuple.
%
  if ( m2 < m1 )
    rank = 0;
    return
  end

  if ( rank <= 0 )

    x(1:n,1) = m1;
    rank = 1;

  else

    rank = rank + 1;
    i = n;

    while ( 1 )

      if ( x(i) < m2 )
        x(i) = x(i) + 1;
        break
      end

      x(i) = m1;

      if ( i == 1 )
        rank = 0;
        x(1:n) = m1;
        break
      end

      i = i - 1;

    end

  end

  return
end

author: johannesgerer


text:
function dockfig(fig)
% Copyright 2008 The MathWorks, Inc.
% Programmatically dock one or all open figures.
%
% SYNTAX:
% dockfig
%    Docks the current figure.
%
% dockfig(fig)
%    Docks figure with handle FIG, and brings it to the front.
%
% dockfig('all')
%    Docks all open figures
%
% Written by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% 01/01/08

% Revision 1
% 01/03/08: Providing default behavior when no input argument is provided
% (operates on current figure).

if nargin == 0
    fig = gcf;
end
if ischar(fig) && strcmpi(fig,'all')
    fig = findall(0,'type','figure');
end
for ii = 1:numel(fig)
    set(fig(ii),'windowstyle','docked');
end
if numel(fig) == 1
    figure(fig);
end

author: Sable

repo: mcbench-benchmarks