text:
function classifier = train_boosted_dt_mc(features, cat_features, labels, ...
    num_iterations, num_nodes, stopval, init_weights, varargin)
% Train a classifier based on boosted decision trees.  Boosting done by the
% logistic regression version of Adaboost (Adaboost.L - Collins, Schapire,
% Singer 2002).  At each
% iteration, a set of decision trees is created for each class, with
% confidences equal to 1/2*ln(P+/P-) for that class, according to the
% weighted distribution.  Final classification is based on the largest
% confidence label (possibly incorporating a prior as h0(c) =
% 1/2*ln(Pc/(1-Pc)).  Weights are assigned as
% w(i,j) = 1 / (1+exp(sum{t in iterations}[yij*ht(xi, j)])).

if length(varargin) == 1  % class names supplied
    gn = varargin{1};
    gid = zeros(size(labels));
    for c = 1:length(gn)
        ind = find(strcmp(labels, gn{c}));
        gid(ind) = c;
        if ~isempty(init_weights)
            disp([gn{c} ': ' num2str(sum(init_weights(ind)))]);
        else
            disp([gn{c} ': ' num2str(length(ind))]);
        end
    end
    ind = find(gid==0);
    gid(ind) = [];
    labels(ind) = [];
    features(ind, :) = [];
else
    [gid, gn] = grp2idx(labels);
    gn
end

classifier.names = gn;

num_classes = length(gn);
num_data = length(gid);

if ~exist('init_weights', 'var') || isempty(init_weights)
    init_weights = ones(num_data, 1)/num_data;
else
    init_weights = init_weights(:) / sum(init_weights);
end

% if no examples from a class are present, create one dummy example for
% that class with very small weight
for c = 1:numel(gn)
    if ~any(gid==c)
        disp(['warning: no examples from class ' gn(c)])
        gid(end+1) = c;
        features(end+1, :) = zeros(size(features(end, 1)));
        num_data = num_data + 1;
        init_weights(end+1) = min(init_weights)/2;
    end
end

all_conf = zeros(num_data, num_classes);
for c = 1:num_classes

    disp(['class: ' num2str(gn{c})]);
    y = (gid == c)*2-1;
    cl = [-1 1];
    nc = 2;
    w = zeros(num_data, 1);
    cw = zeros(num_classes, 1);
    for i = 1:2
        indices = find(y==cl(i));
        %count = sum(init_weights(indices));
        %w(indices) = init_weights(indices) / count / 2;
        w(indices) = init_weights(indices);

        if cl(i)==1
            %classifier.h0(c) = log(count / (1-count));
            classifier.h0(c) = 0;
        end

    end

    data_confidences = zeros(num_data, 1);
    aveconf = [];

    for t = 1:num_iterations
        % learn decision tree based on weighted distribution
        dt = treefitw(features, y, w, 1/num_data/2, 'catidx', cat_features, 'method', 'classification', 'maxnodes', num_nodes*4);
        [tmp, level] = min(abs(dt.ntermnodes-num_nodes));
        dt = treeprune(dt, 'level', level-1);

        % assign partition confidences
        pi = (strcmp(dt.classname{1},'1')) + (2*strcmp(dt.classname{2},'1'));
        ni = (strcmp(dt.classname{1},'-1')) + (2*strcmp(dt.classname{2},'-1'));
        classprob = dt.classprob;
        confidences = 1/2*(log(classprob(:, pi)) - log(classprob(:, ni)));

        % assign weights
        [class_indices, nodes, classes] = treeval(dt, features);
        data_confidences = data_confidences + confidences(nodes);

        w = 1 ./ (1+exp(y.*data_confidences));
        % was w = 1 ./ (1+exp(y.*data_confidences));
        w = w / sum(w);

%         disp(['c: ' num2str(sum(init_weights ./ (1+exp(-y.*data_confidences)))) ...
%             '  e: ' num2str(sum(init_weights .* (y.*data_confidences < 0))) ...
%             '   w: ' num2str(max(w))]);

        classifier.wcs(t, c).dt = dt;
        classifier.wcs(t, c).confidences = confidences;

        %aveconf(t) = mean(1 ./ (1+exp(-y.*data_confidences)));
        aveconf(t) = sum(1 ./ (1+exp(-y.*data_confidences)).*init_weights);
        if t>10 && (aveconf(t)-aveconf(t-10) < stopval)
            disp(num2str(aveconf))
            disp(['Stopping after ' num2str(t) ' trees'])
            break;
        end

    end

    finalconf = 1 ./ (1+exp(-y.*data_confidences)) .* init_weights;
    finalerr = (y.*data_confidences < 0);
    disp(['confidence:: mean: ' num2str(sum(finalconf)) ...
        '  pos: ' num2str(sum(finalconf(y==1))/sum(init_weights(y==1))) ...
        '  neg: ' num2str(sum(finalconf(y~=1))/sum(init_weights(y~=1)))]);
    disp(['training error:: mean: ' num2str(sum(init_weights.*finalerr)) ...
        '  pos: ' num2str(sum(init_weights(y==1).*finalerr(y==1))/sum(init_weights(y==1))) ...
        '  neg: ' num2str(sum(init_weights(y~=1).*finalerr(y~=1))/sum(init_weights(y~=1)))]);
    all_conf(:, c) = data_confidences+classifier.h0(c);

end

% compute and display training error
[tmp, assigned_label] = max(all_conf, [], 2);
conf_matrix = zeros(num_classes, num_classes);
for c = 1:num_classes
    indices = find(gid==c);
    for c2 = 1:num_classes
        conf_matrix(c, c2) = sum(init_weights(indices).*(assigned_label(indices)==c2))/sum(init_weights(indices));
    end
    disp([gn{c} ' error: ' num2str(sum(init_weights(indices).*(assigned_label(indices)~=c))/sum(init_weights(indices)))]);
end
disp('Confusion Matrix: ');
disp(num2str(conf_matrix));
disp(['total error: ' num2str(sum(init_weights.*(assigned_label~=gid)))]);


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%FxNLMS
%Marko Stamenovic
%April 28, 2016

mus = [0.01 0.05 0.1 0.5 1 2 2.5 3];% 0.5 1];
for j = 1:length(mus)
    for i = 1:100
        %%INITIALIZE VALUES%%
        % generate input signal
        muOG=mus(j); %learning rate
        M=128; %buffer size (num filter weights)
        x=randn(10000,1); %input signal
        x=x/max(x); %sample rate
        fs=8000; %number of samples of the input signal
        N=length(x); %length of input signal

        % generate known filter coefficients
        Pz=0.5*[0:127];  %linear coefficients
        %Pz=randn(128,1); %random coefficients
        ylim = max(Pz)*1.20;
        ymin = min(Pz)-.2*max(Pz);
        % generate filtered input signal == desired signal
        d=conv(Pz,x); %input signal filtered by known filter Pz (primary path)

        %% ESTIMATE SECONDARY PATH SIGNAL USING LMS %%
        %generate dummy secondary path response Sz
        Sz = Pz/2;
        %run known signal through filter
        xp=conv(Sz,x);
        %initalize Sz hat values
        Szh=zeros(M,1);

        for n=M:N
            xpvec=x(n:-1:n-M+1); %input has to be in reverse orxer has to be
            %update mu
            mu(n) = 1/(xpvec'*xpvec);
            e(n)=xp(n)-Szh'*xpvec; %update error
            %plot(e)
            Szh=Szh+mu(n)*xpvec*(e(n)); %update filter coefficient
        end
        Szh = abs(ifft(1./abs(fft(Szh))));

        %% LMS FOR MAIN ANC %%
        emean=zeros(N,1);
        %filter input by learned filter to get x prime
        xph = conv(Szh,x);
        %initalize Wz filter values
        Wz=zeros(M,1);
        %Make sure that x and d are column vectors
        x=x(:);
        d=d(:);
        %LMS
        for n=M:N
            xvec=x(n:-1:n-M+1); %input has to be in reverse orxer
            xpvec=xph(n:-1:n-M+1);
            %update mu
            mu(n) = muOG/(xvec'*xvec);
            e(n)=d(n)-Wz'*xvec; %update error
            %plot(e)
            Wz=Wz+mu(n)*xpvec*(e(n)); %update filter coefficient

            %draw the learned filter in realtime
            %         plot(Pz)
            %         hold on
            %         plot(Wz)
            %         axis([0 inf ymin ylim])
            %         title(sprintf('n=%f time=%fs error = %f',n-M, (n-M)/fs, e(n)))
            %         hold off
            %         legend('Input coefficients','Learned Coefficients')
            %         drawnow;
        end
        e=e(:);

        emean = (emean(:)+e);
    end
    emean=(emean)/i;
    if max(emean)>1e4
        for l = 1:length(emean)
            if abs(emean(l)) > 1e3
                emean(l)=1e3;
            end
        end
        eall(j,:)=emean;
    else
        try
            [eall(j,:),q]=(envelope(abs(emean),500,'peaks'));
        catch
            eall(j,:) = 1e3;
        end
    end
end
figure
for i = 1:length(mus)
    plot(10*log10(abs(eall(i,:))))
    hold on
end
title('Convergence Time in Cycles')
ylabel('Error (dB)');
xlabel('Cycles');
hleg=legend('0.01','0.05','0.1','0.5','1.0', '2.0');
htitle=get(hleg,'Title');
set(htitle,'String','mu');
% %% PLOT RESULTS %%
% figure
% subplot(2,1,1)
% plot(e)
% title('Convergence Time in Cycles')
% ylabel('Amplitude');
% xlabel('Cycles');
% legend('Error');
% subplot(2,1,2)
% stem(Pz)
% hold on
% stem(Wz, 'r*')
% title('Input Coefficients vs Learned Coefficients')
% ylabel('Amplitude');
% xlabel('Numbering of filter tap');
% legend('Input Coefficients', 'learned coefficients')

author: markostam

repo: active-noise-cancellation


text:
function fd_predator_prey ( p0, tspan, step_num )

%*****************************************************************************80
%
%% FD_PREDATOR_PREY solves a pair of predator-prey ODE's.
%
%  Discussion:
%
%    The physical system under consideration is a pair of animal populations.
%
%    The PREY reproduce rapidly; for each animal alive at the beginning of the
%    year, two more will be born by the end of the year.  The prey do not have
%    a natural death rate; instead, they only die by being eaten by the predator.
%    Every prey animal has 1 chance in 1000 of being eaten in a given year by
%    a given predator.
%
%    The PREDATORS only die of starvation, but this happens very quickly.
%    If unfed, a predator will tend to starve in about 1/10 of a year.
%    On the other hand, the predator reproduction rate is dependent on
%    eating prey, and the chances of this depend on the number of available prey.
%
%    The resulting differential equations can be written:
%
%      PREY(0) = 5000
%      PRED(0) =  100
%
%      d PREY / dT =    2 * PREY(T) - 0.001 * PREY(T) * PRED(T)
%      d PRED / dT = - 10 * PRED(T) + 0.002 * PREY(T) * PRED(T)
%
%    Here, the initial values (5000,100) are a somewhat arbitrary starting point.
%
%    The pair of ordinary differential equations that result have an interesting
%    behavior.  For certain choices of the interaction coefficients (such as
%    those given here), the populations of predator and prey will tend to
%    a periodic oscillation.  The two populations will be out of phase; the number
%    of prey will rise, then after a delay, the predators will rise as the prey
%    begins to fall, causing the predator population to crash again.
%
%    In this program, the pair of ODE's is solved with a simple finite difference
%    approximation using a fixed step size.  In general, this is NOT an efficient
%    or reliable way of solving differential equations.  However, this program is
%    intended to illustrate the ideas of finite difference approximation.
%
%    In particular, if we choose a fixed time step size DT, then a derivative
%    such as dPREY/dT is approximated by:
%
%      d PREY / dT = approximately ( PREY(T+DT) - PREY(T) ) / DT
%
%    which means that the first differential equation can be written as
%
%      PREY(T+DT) = PREY(T) + DT * ( 2 * PREY(T) - 0.001 * PREY(T) * PRED(T) ).
%
%    We can rewrite the equation for PRED as well.  Then, since we know the
%    values of PRED and PREY at time 0, we can use these finite difference
%    equations to estimate the values of PRED and PREY at time DT.  These values
%    can be used to get estimates at time 2*DT, and so on.  To get from time
%    T_START = 0 to time T_STOP = 5, we simply take STEP_NUM steps each of size
%    DT = ( T_STOP - T_START ) / STEP_NUM.
%
%    Because finite differences are only an approximation to derivatives, this
%    process only produces estimates of the solution.  And these estimates tend
%    to become more inaccurate for large values of time.  Usually, we can reduce
%    this error by decreasing DT and taking more, smaller time steps.
%
%    In this example, for instance, taking just 100 steps gives nonsensical
%    answers.  Using STEP_NUM = 1000 gives an approximate solution that seems
%    to have the right kind of oscillatory behavior, except that the amplitude
%    of the waves increases with each repetition.  Using 10000 steps, the
%    approximation begins to become accurate enough that we can see that the
%    waves seem to have a fixed period and amplitude.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Second Edition,
%    Prentice Hall, 1999,
%    ISBN: 0-13-012641-1,
%    LC: QA297.P45.
%
%  Parameters:
%
%    Input, real P0(2), the initial number of prey and predators.
%
%    Input, real TSPAN(2), the initial and final times.
%
%    Input, integer STEP_NUM, the number of steps.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD_PREDATOR_PREY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A finite difference approximate solution of a pair\n' );
  fprintf ( 1, '  of ordinary differential equations for a population\n' );
  fprintf ( 1, '  of predators and prey.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The exact solution shows wave behavior, with a fixed\n' );
  fprintf ( 1, '  period and amplitude.  The finite difference approximation\n' );
  fprintf ( 1, '  can provide a good estimate for this behavior if the stepsize\n' );
  fprintf ( 1, '  DT is small enough.\n' );

  t_start = tspan(1);
  t_stop =  tspan(2);
  dt = ( t_stop - t_start ) / step_num;

  trf = zeros (3,step_num+1);

  trf(1,1) = t_start;
  trf(2,1) = p0(1);
  trf(3,1) = p0(2);

  for i = 1 : step_num
    trf(1,i+1) = trf(1,i) + dt;
    trf(2,i+1) = trf(2,i) + dt * (    2 * trf(2,i) - 0.004 * trf(2,i) * trf(3,i) );
    trf(3,i+1) = trf(3,i) + dt * ( - 10 * trf(3,i) + 0.003 * trf(2,i) * trf(3,i) );
  end
%
%  Time plots.
%
  figure ( 1 )
  plot ( trf(1,:), trf(2,:), 'g-', ...
         trf(1,:), trf(3,:), 'r-', 'LineWidth', 3 )
  filename = sprintf ( 'trf_%d_time.png', step_num );
  title ( 'Predator Prey System Solved by Finite Differences' );
  grid on
  xlabel ( 'Time' );
  ylabel ( 'Populations' );
  print ( '-dpng', filename );
%
%  Phase plot.
%
  figure ( 2 )
  plot ( trf(2,:), trf(3,:), 'b-', 'LineWidth', 3 )
  filename = sprintf ( 'trf_%d_phase.png', step_num );
  title ( 'Predator Prey System Solved by Finite Differences' );
  grid on
  xlabel ( 'Rabbits' );
  ylabel ( 'Foxes' );
  print ( '-dpng', filename );
%
%  Write data to files.
%
  filename = sprintf ( 'trf_%d.txt', step_num );
  r8mat_write ( filename, 3, step_num + 1, trf );
  fprintf ( 1, '  T, R, F values written to "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD_PREDATOR_PREY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6e', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16e', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function timestamp ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function markedElem = mark(elem,eta,theta,method)
% MARK mark element.
%
% markedElem = mark(elem,eta,theta) mark a subset of elements by Dorfler
% marking strategy. It returns an array of indices of marked elements
% markedElem such that sum(eta(markedElem)^2) > theta*sum(eta^2).
%
% markedElem = mark(elem,eta,theta,'max') choose markedElem such that
% eta(markedElem) > theta*max(eta).
%
% markedElem = mark(elem,eta,theta,'COARSEN') choose markedElem such that
% eta(markedElem) < theta*max(eta).
%
% Copyright (C) 2008 Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1); isMark = false(NT,1);
if ~exist('method','var'), method = 'L2'; end  % default marking is L2 based
switch upper(method)
    case 'MAX'
        isMark(eta>theta*max(eta))=1;
    case 'COARSEN'
        isMark(eta<theta*max(eta))=1;
    case 'L2'
        [sortedEta,idx] = sort(eta.^2,'descend');
        x = cumsum(sortedEta);
        isMark(idx(x < theta* x(NT))) = 1;
        isMark(idx(1)) = 1;
end
markedElem = uint32(find(isMark==true));

repo: Aero-matlab


text:
function [p_min, f_min, iter]=pso(func, numInd, range, n_var, tolerance, numIter, pesoStoc)
% It finds the absolut minimum of a n variables function with the Particle
% Swarm Optimization Algorithm.
% The input parameters are:
% -func: it's the objective function's handle to minimize
% -numInd: it's the number of the swarm's elements
% -range: it's the range in which the elements must be created
% -n_var: it's the number of function's variables
% -tolerance: it's the tolerance for the stop criterion on the swarm's
% radius
% -numIter: it's the max iterations' number
% -pesoStoc: it's the swarm's movability
%
% The output parameters are:
% -p_min: the minimum point find
% -f_min: the minimum value of the function
% -iter: the number of iterations processed

range_min=range(1); % Range for initial swarm's elements
range_max=range(2);
numVar=n_var; % Number of variables
fOb=func; % Objective function
iter=1; % Number of iteration
ind = range_min + (range_max-range_min).*rand(numInd,n_var); % Initial swarm
k=pesoStoc; % weight of stocastic element

v=zeros(numInd,numVar); % Vector of swarm's velocity

radius=1000; % Initial radius for stop criterion

while iter<numIter && radius>tolerance
    for l=1:numInd
        valF(l,1)=fOb(ind(l,:)); % Fitness function for the swarm
    end
    [valF_ord,index]=sort(valF); % Sort the objective function's values for the swarm and identify the leader
    leader=ind(index(1),:);
    for l=1:size(ind,1) % Calculates the new velocity and positions for all swarm's elements
        fi=rand();
        v(l,:)=(1-(sqrt(k*fi))/2)*v(l,:)+k*fi*(leader-ind(l,:)); % Velocity
        ind(l,:)=ind(l,:)+(1-(sqrt(k*fi))/2)*v(l,:)+(1-k*fi)*(leader-ind(l,:)); % Position
    end
    radius=norm(leader-ind(index(end),:)); % Calculates the new radius
    iter=iter+1; % Increases the number of iteration
end

p_min=ind(1:20,:); % Output variables
f_min=valF_ord(1:20,:);

author: Sable

repo: mcbench-benchmarks


text:
function b = r8vm_to_r8ge ( m, n, a, b )

%*****************************************************************************80
%
%% R8VM_TO_R8GE copies a R8VM matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8VM storage format is used for an M by N Vandermonde matrix.
%    An M by N Vandermonde matrix is defined by the values in its second
%    row, which will be written here as X(1:N).  The matrix has a first
%    row of 1's, a second row equal to X(1:N), a third row whose entries
%    are the squares of the X values, up to the M-th row whose entries
%    are the (M-1)th powers of the X values.  The matrix can be stored
%    compactly by listing just the values X(1:N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(N,1), the R8VM matrix.
%
%    Output, real B(M,N), the R8GE matrix.
%
  a = a(:);

  for i = 1 : m
    for j = 1 : n
      if ( i == 1 )
        b(i,j) = 1.0;
      else
        b(i,j) = b(i-1,j) * a(j,1);
      end
    end
  end

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
%
% KERNEL/TOOLS
%
% Contents of FAIR TOOLS
%
%   contents.m          this file
%   dealOptions.m       deals persistent options for use in main modules
%   reportStatus.m      reports the current setting of the toolbox
%
%   testFAIR.m          tests the complete FAIR toolbox
%   testTools.m         test the kernel/tools folder
%
%   FAIRtestPara.m      handles the variables used for testings
%   testEnd.m           finishes    the testing environment
%
%   FAIRcheckFolder.m   generates list of files to be tested
%   FAIRcheckFiles.m    checks all files from a list
%   FAIRmake.m          convenient way of calling the matlab mex compiler
%   FAIRcell2struct.m   converts cell to struct and vice versa
%   FAIRmessage.m       generates nice message
%   FAIRinput.m         input, disabled in the testing framework
%   FAIRpause.m         pause, disabled in the testing framework
%   FAIRpath.m          shortcut to the root of the toolbox
%==============================================================================

function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
  'contents.m'
  'dealOptions.m'
  'reportStatus.m'

  'testFAIR.m'
  'testTools.m'

  'FAIRtestPara.m'
  'testEnd.m'

  'FAIRcheckFolder.m'
  'FAIRcheckFiles.m'
  'FAIRmake.m'

  'FAIRcell2struct.m'
  'FAIRmessage.m'
  'FAIRinput.m'
  'FAIRpause.m'
  'FAIRpath.m'

};
%==============================================================================


text:
function S = DSAFinit_Morgan(M,mu,N,J)

% DSAFinit_Morgan   Initialize Parameter Structure for the Delayless
%                   Subband Adaptive Filter Proposed by Morgan and Thi
% Arguments:
% M                 Length of corresponding fullband filter
% mu                Step size
% N                 Number of subbands
% J                 Update rate, typically in the range 1 to 8
%
% Subband Adaptive Filtering: Theory and Implementation
% by Lee, Gan, and Kuo, 2008
% Publisher: John Wiley and Sons, Ltd

D = N/2;                            % 2x oversampling (recommended)
OverFac = 4;                        % Overlapping factor (Try other integer values)
L = 2*OverFac*N-1;                  % Length of analysis filters, L = 2KN-1

% DFT filter bank with fractional delay (see Section 4.4.3.2)

wc = 1/N;                           % Cut-off frequency
hopt = fir1(L-1,wc);                % Generate prototype filter
[H,F] = make_bank_DFT(hopt,N);      % Complex modulation

% Assign structure fields
S.FULLcoeffs    = zeros(M,1);       % Fullband adaptive filter
S.SUBcoeffs     = zeros(M/D,N/2+1); % Adaptive subfilters (initialized to zeros)
                                    %  consider only the first N/2 + 1 subfilters
S.step          = mu;               % Step size
S.decfac        = D;                % Decimation factor
S.analysis      = H;                % Analysis filter bank
S.synthesis     = F;                % Synthesis filter bank
S.iter          = 0;                % Iteration count
S.alpha         = 1e-4;             % Small positive constant
S.AdaptStart    = L + M;            % Running effect of the analysis and adaptive filter, minimum L+M
S.UpdateRate    = round(M/J);       % J is typically in the range 1 to 8

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function [ a, R, T, T3D, params, error, shapeOrtho ] = fit_PDM_ortho_proj_to_2D( M, E, V, shape2D, f, cx, cy)
%FITPDMTO2DSHAPE Summary of this function goes here
%   Detailed explanation goes here

    params = zeros(size(E));

    hidden = false;

    % if some of the points are unavailable modify M, V, and shape2D (can
    % later infer the actual shape from this)
    if(sum(shape2D(:)==0) > 0)

        hidden = true;
        % which indices to remove
        inds_to_rem = shape2D(:,1) == 0 | shape2D(:,2) == 0;

        shape2D = shape2D(~inds_to_rem,:);

        inds_to_rem = repmat(inds_to_rem, 3, 1);

        M_old = M;
        V_old = V;

        M = M(~inds_to_rem);
        V = V(~inds_to_rem,:);

    end

    num_points = numel(M) / 3;

    m = reshape(M, num_points, 3)';
    width_model = max(m(1,:)) - min(m(1,:));
    height_model = max(m(2,:)) - min(m(2,:));

    bounding_box = [min(shape2D(:,1)), min(shape2D(:,2)),...
                    max(shape2D(:,1)), max(shape2D(:,2))];

    a = (((bounding_box(3) - bounding_box(1)) / width_model) + ((bounding_box(4) - bounding_box(2))/ height_model)) / 2;

    tx = (bounding_box(3) + bounding_box(1))/2;
    ty = (bounding_box(4) + bounding_box(2))/2;

    % correct it so that the bounding box is just around the minimum
    % and maximum point in the initialised face
    tx = tx - a*(min(m(1,:)) + max(m(1,:)))/2;
    ty = ty - a*(min(m(2,:)) + max(m(2,:)))/2;

    R = eye(3);
    T = [tx; ty];

    currShape = getShapeOrtho(M, V, params, R, T, a);

    currError = getRMSerror(currShape, shape2D);

    reg_rigid = zeros(6,1);
    regFactor = 20;
    regularisations = [reg_rigid; regFactor ./ E]; % the above version, however, does not perform as well
    regularisations = diag(regularisations)*diag(regularisations);

    red_in_a_row = 0;

    for i=1:1000

        shape3D = M + V * params;
        shape3D = reshape(shape3D, numel(shape3D) / 3, 3);

        % Now find the current residual error
        currShape = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
        currShape = currShape';

        error_res = shape2D - currShape;

        eul = Rot2Euler(R);

        p_global = [a; eul'; T];

        % get the Jacobians
        J = CalcJacobian(M, V, params, p_global);

        % RLMS style update
        p_delta = (J'*J + regularisations) \ (J'*error_res(:) - regularisations*[p_global;params]);

        [params, p_global] = CalcReferenceUpdate(p_delta, params, p_global);

        a = p_global(1);
        R = Euler2Rot(p_global(2:4));
        T = p_global(5:6);

        shape3D = M + V * params;
        shape3D = reshape(shape3D, numel(shape3D) / 3, 3);
        currShape = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
        currShape = currShape';

        error = getRMSerror(currShape, shape2D);

        if(0.999 * currError < error)
            red_in_a_row = red_in_a_row + 1;
            if(red_in_a_row == 5)
                break;
            end
        end

        currError = error;

    end

    if(hidden)
        shapeOrtho = getShapeOrtho(M_old, V_old, params, R, T, a);
    else
        shapeOrtho = currShape;
    end
    if(nargin == 7)

        Zavg = f / a;
        Xavg = (T(1) - cx) / a;
        Yavg = (T(2) - cy) / a;

        T3D = [Xavg;Yavg;Zavg];
    else
        T3D = [0;0;0];
    end

end

function [shape2D] = getShapeOrtho(M, V, p, R, T, a)

    % M - mean shape vector
    % V - eigenvectors
    % p - parameters of non-rigid shape
    % R - rotation matrix
    % T - translation vector (tx, ty)
    shape3D = getShape3D(M, V, p);
    shape2D = a * R(1:2,:)*shape3D' + repmat(T, 1, numel(M)/3);
    shape2D = shape2D';
end

function [shape2D] = getShapeOrthoFull(M, V, p, R, T, a)

    % M - mean shape vector
    % V - eigenvectors
    % p - parameters of non-rigid shape
    % R - rotation matrix
    % T - translation vector (tx, ty)
    T = [T; 0];
    shape3D = getShape3D(M, V, p);
    shape2D = a * R*shape3D' + repmat(T, 1, numel(M)/3);
    shape2D = shape2D';
end

function [shape3D] = getShape3D(M, V, params)

    shape3D = M + V * params;
    shape3D = reshape(shape3D, numel(shape3D) / 3, 3);

end

function [error] = getRMSerror(shape2Dv1, shape2Dv2)

    error = sqrt(mean(reshape(shape2Dv1 - shape2Dv2, numel(shape2Dv1), 1).^2));

end

% This calculates the combined rigid with non-rigid Jacobian
function J = CalcJacobian(M, V, p, p_global)

    n = size(M, 1)/3;

    non_rigid_modes = size(V,2);

    J = zeros(n*2, 6 + non_rigid_modes);

    % now the layour is
    % ---------- Rigid part -------------------|----Non rigid part--------|
    % dx_1/ds, dx_1/dr1, ... dx_1/dtx, dx_1/dty dx_1/dp_1 ... dx_1/dp_m
    % dx_2/ds, dx_2/dr1, ... dx_2/dtx, dx_2/dty dx_2/dp_1 ... dx_2/dp_m
    % ...
    % dx_n/ds, dx_n/dr1, ... dx_n/dtx, dx_n/dty dx_n/dp_1 ... dx_n/dp_m
    % dy_1/ds, dy_1/dr1, ... dy_1/dtx, dy_1/dty dy_1/dp_1 ... dy_1/dp_m
    % ...
    % dy_n/ds, dy_n/dr1, ... dy_n/dtx, dy_n/dty dy_n/dp_1 ... dy_n/dp_m

    % getting the rigid part
    J(:,1:6) = CalcRigidJacobian(M, V, p, p_global);

    % constructing the non-rigid part
    R = Euler2Rot(p_global(2:4));
    s = p_global(1);

    % 'rotate' and 'scale' the principal components

    % First reshape to 3D
    V_X = V(1:n,:);
    V_Y = V(n+1:2*n,:);
    V_Z = V(2*n+1:end,:);

    J_x_non_rigid = s*(R(1,1)*V_X + R(1,2)*V_Y + R(1,3)*V_Z);
    J_y_non_rigid = s*(R(2,1)*V_X + R(2,2)*V_Y + R(2,3)*V_Z);

    J(1:n, 7:end) = J_x_non_rigid;
    J(n+1:end, 7:end) = J_y_non_rigid;

end

function J = CalcRigidJacobian(M, V, p, p_global)

 	n = size(M, 1)/3;

	% Get the current 3D shape (not affected by global transform, as this
	% is how the Jacobian was derived (for derivation please see
	% ../derivations/orthoJacobian
	shape3D = GetShape3D(M, V, p);

	% Get the rotation matrix corresponding to current global orientation
	R = Euler2Rot(p_global(2:4));
	s = p_global(1);

    % Rigid Jacobian is laid out as follows
    % dx_1/ds, dx_1/dr1, dx_1/dr2, dx_1/dr3, dx_1/dtx, dx_1/dty
    % dx_2/ds, dx_2/dr1, dx_2/dr2, dx_2/dr3, dx_2/dtx, dx_2/dty
    % ...
    % dx_n/ds, dx_n/dr1, dx_n/dr2, dx_n/dr3, dx_n/dtx, dx_n/dty
    % dy_1/ds, dy_1/dr1, dy_1/dr2, dy_1/dr3, dy_1/dtx, dy_1/dty
    % ...
    % dy_n/ds, dy_n/dr1, dy_n/dr2, dy_n/dr3, dy_n/dtx, dy_n/dty

    J = zeros(n*2, 6);

    % dx/ds = X * r11  + Y * r12 + Z * r13
    % dx/dr1 =  s*(r13 * Y - r12 * Z)
    % dx/dr2 = -s*(r13 * X - r11 * Z)
    % dx/dr3 =  s*(r12 * X - r11 * Y)
    % dx/dtx = 1
    % dx/dty = 0

    % dy/ds = X * r21  + Y * r22 + Z * r23
    % dy/dr1 =  s * (r23 * Y - r22 * Z)
    % dy/dr2 = -s * (r23 * X - r21 * Z)
    % dy/dr3 =  s * (r22 * X - r21 * Y)
    % dy/dtx = 0
    % dy/dty = 1

    % set the Jacobian for x's

    % with respect to scaling factor
    J(1:n,1) = shape3D * R(1,:)';

    % with respect to angular rotation around x, y, and z axes

    % Change of x with respect to change in axis angle rotation
	dxdR = [      0,  R(1,3), -R(1,2);
            -R(1,3),       0,  R(1,1);
             R(1,2), -R(1,1),      0];

    J(1:n,2:4) = s*(dxdR * shape3D')';

    % with respect to translation
    J(1:n,5) = 1;
    J(1:n,6) = 0;

    % set the Jacobian for y's

    % with respect to scaling factor
    J(n+1:end,1) = shape3D * R(2,:)';

     % with respect to angular rotation around x, y, and z axes

    % Change of y with respect to change in axis angle rotation
    dydR = [      0,  R(2,3), -R(2,2);
            -R(2,3),       0,  R(2,1);
             R(2,2), -R(2,1),      0];

     J(n+1:end,2:4) = s*(dydR * shape3D')';

     % with respect to translation
     J(n+1:end,5) = 0;
     J(n+1:end,6) = 1;

end

% This updates the parameters based on the updates from the RLMS
function [non_rigid, rigid] = CalcReferenceUpdate(params_delta, current_non_rigid, current_global)

    rigid = zeros(6, 1);
    % Same goes for scaling and translation parameters
    rigid(1) = current_global(1) + params_delta(1);
    rigid(5) = current_global(5) + params_delta(5);
    rigid(6) = current_global(6) + params_delta(6);

    % for rotation however, we want to make sure that the rotation matrix
    % approximation we have
	% R' = [1, -wz, wy
	%       wz, 1, -wx
	%       -wy, wx, 1]
    % is a legal rotation matrix, and then we combine it with current
    % rotation (through matrix multiplication) to acquire the new rotation

	R = Euler2Rot(current_global(2:4));

    wx = params_delta(2);
    wy = params_delta(3);
    wz = params_delta(4);

    R_delta = [1, -wz, wy;
               wz, 1, -wx;
               -wy, wx, 1];

	% Make sure R_delta is orthonormal
	R_delta = OrthonormaliseRotation(R_delta);

    % Combine rotations
	R_final = R * R_delta;

	% Extract euler angle
	euler = Rot2Euler(R_final);

	rigid(2:4) = euler;

    if(length(params_delta) > 6)
        % non-rigid parameters can just be added together
        non_rigid = params_delta(7:end) +  current_non_rigid;
    else
        non_rigid = current_non_rigid;
    end

end

function R_ortho = OrthonormaliseRotation(R)

    % U * V' is basically what we want, as it's guaranteed to be
    % orthonormal
    [U, ~, V] = svd(R);

    % We also want to make sure no reflection happened

    % get the orthogonal matrix from the initial rotation matrix
    X = U*V';

    % This makes sure that the handedness is preserved and no reflection happened
    % by making sure the determinant is 1 and not -1
    W = eye(3);
    W(3,3) = det(X);
    R_ortho = U*W*V';
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function poly2 = resamplePolylineByLength(poly, step)
%RESAMPLEPOLYLINEBYLENGTH Resample a polyline with a fixed sampling step.
%
%   RES = resamplePolyline(POLY, STEP)
%   Resample the input polyline POLY by distributing new vertices on the
%   original polyline such that the (curvilinear) distance between the new
%   vertices is approximately equal to STEP.
%
%   Example
%     poly = [0 10;0 0;10 0; 10 10; 20 10;20 0];
%     figure; drawPolyline(poly, 'k');
%     poly2 = resamplePolylineByLength(poly, 4);
%     hold on;
%     drawPolyline(poly2, 'm');
%     drawPoint(poly2, 'mo');
%     axis equal; axis([-10 30 -10 20]);
%     legend('Original polyline', 'Resampled polyline');
%
%   See also
%     polygons2d, drawPolyline, resamplePolygon
%

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

% parametrisation of the curve
s = parametrize(poly);

% compute the number of points for sampling the polygon
% (equal to the number of segments plus one)
Lmax = s(end);
n = round(Lmax / step) + 1;

% distribute N points equally spaced
pos = linspace(0, Lmax, n);

poly2 = zeros(n, size(poly, 2));
for i = 1:n
    % index of surrounding vertices before and after
    ind0 = find(s <= pos(i), 1, 'last');
    ind1 = find(s >= pos(i), 1, 'first');

    if ind0 == ind1
        % get position of a vertex in input polyline
        poly2(i, :) = poly(ind0, :);
        continue;
    end

    % position of surrounding vertices
    pt0 = poly(ind0, :);
    pt1 = poly(ind1, :);

    % weights associated to each neighbor
    l0 = pos(i) - s(ind0);
    l1 = s(ind1) - pos(i);

    % linear interpolation of neighbor positions
    if (l0 + l1) > Lmax * 1e-12
        poly2(i, :) = (pt0 * l1 + pt1 * l0) / (l0 + l1);
    else
        % if neighbors are too close, do not use interpolation
        poly2(i, :) = pt0;
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function newcat2=plotmi(var1, newcat2, mi)
    %plot misfit (?)
    % TODO make this work with the new catalogs
    report_this_filefun(mfilename('fullpath'));

    global  mif2 mif1
    global tmp % REALLY? global tmp?  "tmp" is 1:nEvents
    % cumu2 mi2
    figNumber=findobj('Type','Figure','-and','Name','Misfit ');
    figure(figNumber);
    delete(findobj(figNumber,'Type','axes'));
    rect = [0.15,  0.15, 0.75, 0.65];
    axes('position',rect)
    ax=gca;
    nEvents=newcat2.Count;
    tmp=1:nEvents;
    sixSlices=round(0 : nEvents/5 : nEvents);
    sixSlices(1)=1;

    var2=var1;

    misfitAngle = mi(:,2);
    X = 1:nEvents;
    xtitle=sprintf('Number of Eqs (sorted by %s)',lower(var1));
    switch (var1)
        case {'Longitude','Latitude','Magnitude','Depth'}
            % plot_by_lon(); %by lon
            plot_by_field(var1);
        case 'Time'
            plot_by_time(); % by date
        case 'Strike'
            plot_by_strike(); % along strike
        case 'Default'
            option_7(); %unsorted
        otherwise
            error('unknown choice for plotmi');
    end

    grid('on')
    set(ax,'Color',color_bg);
    set(ax,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    ylabel('Cumulative Misfit ','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m);
    xlabel(xtitle,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)

    if ~ strcmp(var1 ,'Default')
        hold on
        for i=1:6
            plot(ax,tmp(sixSlices(i)),cumu2(sixSlices(i)),'xr');
            str=['  ',num2str(newcat2(sixSlices(i),var2))];
            te=text(tmp(sixSlices(i)),cumu2(sixSlices(i)),str);
            set(te,'FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.s)
        end
    else

    end

    function plot_by_field(name)
        % assumes that misfit matrix (mi) has same number of rows as
        % number of earthquakes in catalog
        [~,is] = sort(newcat2.(name));
        newcat2.sort(name); % sort the catalog itself by this field

        cumu2=cumsum(misfitAngle(is));
        plot(1:nEvents , cumu2 , 'o');
        xtitle=sprintf('Number of Eqs (sorted by %s)',lower(name));
    end

    function plot_by_time()
        [~,is] = sort(newcat2.Date);
        newcat2.sort(Date);
        cumu2=cumsum(misfitAngle(is));
        pl = plot(tmp,cumu2,'o');
        xtitle='Number of Eqs (sorted by time)';
    end

    function plot_by_strike()
        % [~,is] = sort(newcat2(:,15));
        [~,is] = sort(newcat2(:,end));
        newa2 = newcat2.subset(is) ;
        cumu2=cumsum(misfitAngle(is));
        pl = plot(newa2(:,16)-18.6,cumu2,'o');
        xtitle='Number of Eqs (sorted along strike)';
        var2=15;
    end

    function option_7()
        mi2 = mi ;
        cumu2=cumsum(mi2(:,2));
        pl = plot(tmp,cumu2,'o');
        xtitle='Number of Eqs ';
    end
end

author: CelsoReyes


text:
function [out] = smoothThreshold_temperature_logistic(T,Tt,r)
%smoothThreshold_temperature_logistic Logisitic smoother for temperature threshold functions.

% Copyright (C) 2018 Wouter J.M. Knoben
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

%   Smooths the transition of threshold functions of the form:
%
%   Snowfall = { P, if T <  Tt
%              { 0, if T >= Tt
%
%   By transforming the equation above to Sf = f(P,T,Tt,r):
%   Sf = P * 1/ (1+exp((T-Tt)/r))
%
%   Inputs:
%   P       : current precipitation
%   T       : current temperature
%   Tt      : threshold temperature below which snowfall occurs
%   r       : [optional] smoothing parameter rho, default = 0.01
%
%   NOTE: this function only outputs the multiplier. This needs to be
%   applied to the proper flux utside of this function.

% Check for inputs and use defaults if not provided
% NOTE: this is not very elegant, but it is more than a factor 10 faster then:
% if ~exist('r','var'); r = 0.01; end
if nargin == 2
    r = 0.01;
end

% Calculate multiplier
out = 1 ./ (1+exp((T-Tt)/(r)));

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function [ class, clsize, critvl, ntrans, ifault ] = trnsfr ( varval, class, ...
  clsize, in, ik, iv, critvl )

%*****************************************************************************80
%
%% TRNSFR transfers objects between classes to improve a criterion.
%
%  Discussion:
%
%    This routine is given a classification of objects, including the
%    number of objects in each class, and the current value of some criterion
%    which is desired to be minimized.
%
%    The routine calculates the change in criterion for all possible transfers
%    of any object from its current class to a different class.  Each transfer
%    that would result in a lowering of the criterion is executed, and the
%    related quantities are updated.
%
%    When no more advantageous transfers can be found, the routine returns.
%
%    The routine relies on a user-supplied routine, CRTRAN, to report the
%    expected change in the criterion for a given transfer, and to carry
%    out that transfer if requested.
%
%    The variables CLASS and CRITVL have been added to the argument list
%    of CRTRAN.
%
%    Also, the order of the two classes "L" and "M" was interchanged in
%    the call to CRTRAN.  The original order was counterintuitive.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by Banfield, Bassill.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Colin Banfield, LC Bassill,
%    Algorithm AS 113:
%    A transfer for non-hierarchichal classification,
%    Applied Statistics,
%    Volume 26, Number 2, 1977, pages 206-210.
%
%  Parameters:
%
%    Input, real VARVAL(IN,IV), the data values.  There are IN
%    objects, each having spatial dimension IV.
%
%    Input, integer CLASS(IN), the initial classification of
%    each object.
%
%    Input, integer CLSIZE(IK), the initial number of objects in
%    each class.
%
%    Input, integer IN, the number of objects.
%
%    Input, integer IK, the number of classes.
%
%    Input, integer IV, the number of spatial dimensions, or
%    variates, of the objects.
%
%    Input, real CRITVL, the initial value of the criterion.
%
%    Output, integer CLASS(IN), the classification of
%    each object.
%
%    Output, integer CLSIZE(IK), the number of objects in
%    each class.
%
%    Output, real CRITVL, the current value of the criterion.
%
%    Output, integer NTRANS, the number of transfers executed.
%
%    Output, integer IFAULT, error indicator.
%    0, no error detected.
%    1, the number of classes was less than 2.
%    2, the number of objects was less than the number of classes.
%
  eps = 1.0E-38;
  ntrans = 0;

  if ( ik <= 1 )
    ifault = 1;
    return
  end

  if ( in <= ik )
    ifault = 2;
    return
  end

  ifault = 0;
  i = 0;
  icount = 0;

  while ( 1 )

    i = i + 1;

    if ( in <= icount )
      break;
    end

    if ( in < i )
      i = 0;
      icount = 0;
      continue
    end

    m = class(i);
    if ( clsize(m) <= 1 )
      icount = icount + 1;
      continue
    end

    inco = - eps;
    lo = m;
%
%  Test the transfer of object I from class M to class L.
%
    for l = 1 : ik

      if ( l ~= m )

        iswitch = 1;
        inc = crtran ( varval, class, clsize, in, ik, iv, critvl, ...
          i, m, l, iswitch );
%
%  Remember the values of L and INC.
%
        if ( inc < inco )
          lo = l;
          inco = inc;
        end

      end

    end

    icount = icount + 1;
%
%  Execute the transfer of object I from class M to class LO.
%
    if ( lo ~= m )

      l = lo;
      critvl = critvl + inco;
      icount = 0;

      iswitch = 2;
       crtran ( varval, class, clsize, in, ik, iv, critvl, ...
        i, m, l, iswitch );

      ntrans = ntrans + 1;
      class(i) = l;
      clsize(l) = clsize(l) + 1;
      clsize(m) = clsize(m) - 1;

    end

  end

  return
end


text:
function Retain_Img = RemoveMinorCC(SegImg,reject_T)
    if nargin < 2
        reject_T = 0.2;
    end
    Retain_Img = zeros(size(SegImg));

    % class 2
    Img_C2 = (SegImg==2);
    Retain_C2 = Connection_Judge_3D(Img_C2,reject_T);
    Retain_C2 = imfill(Retain_C2,'hole');
    Retain_Img(Retain_C2) = 2;

    % class 1
    Img_C1 = (SegImg==1);
    Retain_C1 = Connection_Judge_3D(Img_C1,reject_T);
    Retain_C1 = imfill(Retain_C1,'hole');
    Retain_Img(Retain_C1) = 1;
end

author: yulequan

repo: HeartSeg

save_path: github-repos/MATLAB/yulequan-HeartSeg


text:
function Out = QualityIndices(I_HS,I_REF,ratio)
%--------------------------------------------------------------------------
% Quality Indices
%
% USAGE
%   Out = QualityIndices(I_HS,I_REF,ratio)
%
% INPUT
%   I_HS  : target HS data (rows,cols,bands)
%   I_REF : reference HS data (rows,cols,bands)
%   ratio : GSD ratio between HS and MS imagers
%
% OUTPUT
%   Out.cc   : CC
%   Out.sam  : SAM
%   Out.rmse : RMSE
%   Out.ergas: ERGAS
%
%--------------------------------------------------------------------------
[rows,cols,bands] = size(I_REF);

%Remove border from the analysis
I_HS  = I_HS(ratio+1:rows-ratio,ratio+1:cols-ratio,:);
I_REF = I_REF(ratio+1:rows-ratio,ratio+1:cols-ratio,:);

cc = CC(I_HS,I_REF);
Out.cc = mean(cc);
Out.cc_std = std(cc);
[angle_SAM,map] = SAM(I_HS,I_REF);
Out.sam = angle_SAM;
Out.rmse = RMSE(I_HS,I_REF);
Out.ergas = ERGAS(I_HS,I_REF,ratio);

disp(['CC   : ', num2str(Out.cc), ' (std: ',num2str(Out.cc_std),')']);
disp(['SAM  : ' num2str(Out.sam)]);
disp(['RMSE : ' num2str(Out.rmse)]);
disp(['ERGAS: ' num2str(Out.ergas)]);

author: zhouyuanzxcv

repo: Hyperspectral


text:
filename_bag='experimento_from_min_sgo.bag'
filename_mat='experimento_from_min_sgo.mat'
topic = '/robot/joint_states'
bagMsgs=rosbag(filename_bag);
bagselect = select(bagMsgs,'Topic',topic);

msgStructs = readMessages(bagselect);
pathq=[];
pathqd=[];
for i = 1:size(msgStructs,1)
    msg = msgStructs{i};
    q=msg.Position;
    q_s=q(2:8);
    pathq = [pathq q_s];

    qd=msg.Velocity;
    qd_s=qd(2:8);
    pathqd = [pathqd qd_s];
end
figure, plot(pathq')
save(filename_mat)
'ended'

repo: ARTE


text:
function [imgPyr, maskPyr, scaleImgPyr] = sc_create_pyramid(img, mask, optS)

% SC_CREAT_IMG_PYRAMID
%
% Create image pyramid with linear or log scale for coarse to fine image
% completion
%
% Input:
%   - img:  input image with hole
%   - mask: Hole mask
%   - optS: options
% Output:
%   - imgPyr:      Image pyramid
%   - maskPyr:     Mask pyramid
%   - scaleImgPyr: Image dimensions in each level

% Image size in the high-resolution image
[imgHeight, imgWidth, nCh] = size(img);
img = sc_init_coarsest_level(img, logical(mask));

% =========================================================================
% Create pyramid: scale
% =========================================================================
scaleImgPyr = sc_create_scale_pyramid(imgHeight, imgWidth, optS);

% =========================================================================
% Create pyramid: mask
% =========================================================================
maskPyr = sc_create_image_pyramid(mask, scaleImgPyr, 'mask', optS);

% =========================================================================
% Create pyramid: image
% =========================================================================
imgPyr = sc_create_image_pyramid(img, scaleImgPyr, 'image', optS);

% =========================================================================
% Recover image boundary
% =========================================================================

% nCh = 3;
% for iLvl = 1: optS.numPyrLvl
%     maskCur  = maskPyr{iLvl};
%     bdRegion = maskCur < 0.99 & maskCur > 0.1;
%     bdRegionC = bdRegion(:,:,ones(nCh,1));
%
%     imgCur   = imgPyr{iLvl};
%     imgCurBd = bsxfun(@rdivide, imgCur, maskCur);
%     imgCur(bdRegionC) = imgCurBd(bdRegionC);
%     imgPyr{iLvl} = imgCur;
% end

% Initialize the coarsest level
imgPyr{optS.numPyrLvl} = sc_init_coarsest_level(imgPyr{optS.numPyrLvl}, maskPyr{optS.numPyrLvl});

% Convert to single type
imgPyr = cellfun(@im2single, imgPyr, 'UniformOutput', false);

end

function img = sc_init_coarsest_level(img, mask)

% Get the inital solution
[~, idMap] = bwdist(~mask, 'euclidean');

% Intepolate only in the interior to avoid dark values near the image borders
maskInt = mask;
maskInt(1,:) = 0;   maskInt(end,:) = 0;
maskInt(:,1) = 0;   maskInt(:,end) = 0;

for ch = 1: 3
    imgCh = img(:,:,ch);
    imgCh = imgCh(idMap);
    img(:,:,ch) = roifill(imgCh, maskInt);
end

end

function imgPyr = sc_create_image_pyramid(img, scaleImgPyr, imageType, optS)

% h = fspecial('gaussian', 5, 1);

% Initialize image pyramid
imgPyr  = cell(optS.numPyrLvl, 1);

% The finest level
imgPyr{1} = img;

%
for iLvl = 2: optS.numPyrLvl
    imgHCurLvl = scaleImgPyr{iLvl}.imgSize(1);
    imgWCurLvl = scaleImgPyr{iLvl}.imgSize(2);

    % Previous layer
    imgCur = imgPyr{iLvl - 1};

    % Anti-alising by blurring
%     imgCur   = imfilter(imgCur, h, 'same', 'replicate', 'conv');

    % Resampling
    imgPyr{iLvl} = imresize(imgCur,  [imgHCurLvl, imgWCurLvl], optS.resampleKernel);
end

if(strcmp(imageType, 'mask'))
    % Convert resampled masks into logical type
    for iLvl = 1: optS.numPyrLvl
        imgPyr{iLvl} = imgPyr{iLvl} > 0.5;
    end
end

end

function scaleImgPyr = sc_create_scale_pyramid(imgHeight, imgWidth, optS)

% Compute the coarsest image scale
imgSizeMin   = min(imgHeight, imgWidth);
coarestScale = optS.coarestImgSize/imgSizeMin;

% Compute the scale in each layer in the image pyramid
if(optS.useLogScale)      % use log scale
    scalePyr = 2.^linspace(0, log2(coarestScale), optS.numPyrLvl);
else                      % use linear scale
    scalePyr = linspace(1, coarestScale, optS.numPyrLvl);
end

% Image size in each layer
imgHPyr = round(imgHeight *scalePyr);
imgWPyr = round(imgWidth  *scalePyr);

% Initialize scales
scaleImgPyr = cell(optS.numPyrLvl, 1);

% Finest level
scaleImgPyr{1}.imgScale = 1;
scaleImgPyr{1}.imgSize = [imgHeight, imgWidth];

% Downsampled images
for k = 2: optS.numPyrLvl
    scaleImgPyr{k}.imgScale = scalePyr(k);
    scaleImgPyr{k}.imgSize  = [imgHPyr(k), imgWPyr(k)];
end

end

repo: StructCompletion


text:
function w = rotation_mat_vector_3d ( a, v )

%*****************************************************************************80
%
%% ROTATION_MAT_VECTOR_3D applies a marix rotation to a vector in 3d.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(3,3), the matrix defining the rotation.
%
%    Input, real V(3), the vector to be rotated.
%
%    Output, real W(3), the rotated vector.
%
  dim_num = 3;
%
%  We want V to be a ROW vector.
%  If MATLAB idiocies make it look like a COLUMN vector, please fix that.
%
  if ( size ( v, 2 ) == 1 )
    v = v';
  end

  w(1:dim_num) = a(1:dim_num,1:dim_num) * v(1:dim_num)';

  return
end


text:
function [ds, g, mystd, d, d2, c, c2, mi, b, eigv, eigval] = compare_subjects(varargin)
% This function compares a set of images to one another and does some diagnostics on the similarity among images.
% - It returns multivariate distances and dissimilarities among images
% - It works on the GLOBAL signal after standardizing each image (case 1) or the REGIONAL values in each cluster (case 2)
% - You can also enter a reference image, in which case each image will be correlated with the ref.
%
%
% :Usage:
% ::
%
%     function [ds, g, mystd, d, d2, c, c2, mi, b, eigv, eigval] = compare_subjects([img files or clusters], [mask], ...
%                                        [plot flag], [title on figure], [standardize flag], [text labels], [ref image])
%
% :Inputs:
%
%     a list of image names to compare
%
%     OR
%
%     a clusters structure, with data to compare
%     in timeseries field
%
% If a mask is entered, only voxels in the mask (e.g., with value of 1) will be used.
% You can use this option to specify brain-only or gray-matter only voxels
%
% textlab: optional text labels for each image, can be empty []
%
% If a ref image is entered, each image will be correlated with the ref,
% and values will be saved for the correlation (plot 2 will show these values)
% Useful for comparing anatomical imgs with template, etc.
%
% :Outputs:  from correls with ref image are in variable "c"
%
%   **ds:**
%        multivariate distance (sim. to Mahalanobis) for each image
%        ds is a matrix of squared distances, case numbers, and
%        expected chi2 values (in columns in this order) rows are cases
%
%   **g:**
%        global value for each image
%
%   **d:**
%        global distance from mean image
%        distance, or dissimilarity, is the average absolute deviation between images
%
%   **d2:**
%        matrix of distances among all images
%
%   **c:**
%        correlation between real valued voxels and mean image
%
%   **c2:**
%        correlations among all images (treating voxels as cases)
%
%   **mi:**
%        mutual information between images, with hist2.m
%
%   **b:**
%        principal component scores on correlation matrix for eigenvalues > 1
%
%   **eigv:**
%        eigenvectors
%
%   **eigval:**
%        eigenvalues
%
% :Examples:
% ::
%
%    % Compare normalized anatomcals with standard brain
%    P = get_filename2(['sub*\Anatomy\nscalped_ft1.img']);
%    [ds, g, mystd, d, d2, c, c2, mi] = compare_subjects(P, which('brain_avg152T1.img'), 1, 'intext_countloc', 1, [], which('avg152T1.img'));
%
% ..
%    by Tor Wager
% ..

    doplot = 1;
    dostd = 0;
    do_mi = 1;
    do_corr = 1;
    do_chi2 = 1;
    textlab = [];
    if isempty(varargin)
        error('no inputs.');
    end

    if length(varargin) > 2, doplot = varargin{3}; end
    if length(varargin) > 3, mytitle = varargin{4}; end
    if length(varargin) > 4, dostd = varargin{5}; end
    if length(varargin) > 5, textlab = varargin{6}; end
    if length(varargin) > 6, refimg = varargin{7}; end

    for i=8:length(varargin)
        if(ischar(varargin{i}))
            switch(varargin{i})
                case 'no_mi'
                    do_mi = 0;
                case 'no_corr'
                    do_corr = 0;
                case 'no_chi2'
                    do_chi2 = 0;
            end
        end
    end

    if isstruct(varargin{1})
        clusters = varargin{1};
        error('No method implemented yet.');
    else
        hP = varargin{1};

        disp(' compare subjects.m Running ');

        % Check for images whose dims do not match
        % -----------------------------------------------------------------
        V = spm_vol(hP);
        [bad, whbad] = iimg_check_volinfo(V(1), V);

        hP(whbad, :) = [];
        V(whbad) = [];

        if any(bad)
            disp('Warning!!! Some image dims do not match.  Removing subjects:');
            fprintf('%3.0f ', whbad);
            fprintf('\n');

            if ~isempty(textlab),
                fprintf('%s ', textlab{whbad});
                fprintf('\n');
                textlab(whbad) = [];
            end

        end

        % -----------------------------------------------------------------

        fprintf(1, '\nLoading volumes.\t');

        v = spm_read_vols(V);
        num_images = size(v, 4);

        if length(varargin) > 1
            mP = varargin{2};
            fprintf(1, 'masking volumes.\t')
            % ----------------------------------------------------------------------------------
            % * load mask, and mask all subjects' contrast values
            % so that we show only those voxels existing for all subjects.
            % ----------------------------------------------------------------------------------

            vm = spm_read_vols(spm_vol(mP));
            tmp = size(v);
            if any(size(vm) - tmp(1:3)),
                fprintf(1, '(reslicing mask).\t')
                [tmp, mP] = reslice_imgs(deblank(hP(1,:)), mP);
                vm = spm_read_vols(spm_vol(mP));
            end

            vm = (vm ~= 0);
            vm = +vm;
            vm(vm == 0) = NaN;

            %attempting something less memory-intensive - commented out line
            %required 3 times the size of v in memory... approx 1.2 gigs
            %v = v .* repmat(vm, [1 1 1 size(v, 4)]);
            for i=1:num_images
                v(:,:,:,i) = v(:,:,:,i) .* vm;
            end

            if length(varargin) > 6
                rimg = spm_read_vols(spm_vol(refimg));
                tmp = size(v);
                if any(size(rimg) - tmp(1:3)),
                    fprintf(1, '(reslicing ref image).\t')
                    [tmp, refimg] = reslice_imgs(deblank(hP(1,:)), refimg);
                    rimg = spm_read_vols(spm_vol(refimg));
                end
                rimg(rimg == 0) = NaN;
            end
        end

        % ----------------------------------------------------------------------------------
        % * Standardize, if asked for
        % ----------------------------------------------------------------------------------
        if dostd
            fprintf('\nScaling images to mean 0 and var 1.\t')
        end

        for i = 1:num_images
            tmp = v(:,:,:,i);
            tmp = tmp(:);
            mystd(i) = std(tmp(~isnan(tmp)));
            g(i) = mean(tmp(~isnan(tmp) & tmp ~= 0));

            if dostd
                v(:,:,:,i) = (v(:,:,:,i) - g(i)) ./ mystd(i);
            end
        end
        clear tmp

        % Metrics
        [d, d2, g, mystd] = get_dist(v);

        c = [];
        c2 = [];
        mi = [];
        if exist('rimg', 'var')
            [c, c2, mi] = get_correl(v, do_corr, do_mi, rimg);
        else
            [c, c2, mi] = get_correl(v, do_corr, do_mi);
        end

        if(do_chi2)
            [b, eigv, eigval] = pc(c2, doplot);
            [ds, S, p] = multivar_dist(b);
        else
            ds = [];
            S = [];
            p = [];
            b = [];
            eigv = [];
            eigval = [];
        end

        % plotting stuff

        if doplot
            plot_results();
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Inline Functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function plot_results()
        figure('Color', 'w'); subplot(2, 2, 1);
        bar(g);
        xlabel('Image');
        ylabel('Global for each image');
        hold on;
        plot([0 num_images], [mean(g) mean(g)], 'r-');
        if length(varargin) > 3
            title(mytitle);
        end

        subplot(2, 2, 2);
        bar(mi);
        xlabel('Image');
        if length(varargin) < 7
            varargin{7} = 'refimg';
        end

        if exist('rimg', 'var')
            ylabel('Mutual info with reference');
            title(varargin{7});
        else
            ylabel('Mutual info with mean');
        end

        subplot(2, 2, 3);
        imagesc(d2);
        colormap hot;
        xlabel('Image');
        ylabel('Image');
        title('Avg abs dist');

        if isempty(textlab)
            for j = 1:size(b, 1)
                textlab{j} = num2str(j);
            end
        end

        % mds-like (pca version) on similarities (correlations)
        subplot(2, 2, 4); hold on
        if size(b, 2) > 1
            plot(b(:,1), b(:,2), 'Color', 'w');
            xlabel('Component 1');
            ylabel('Component 2');
            for j = 1:size(b, 1)
                text(b(j, 1), b(j, 2), textlab{j}, 'Color', 'b');
            end
        else
            plot(ones(1, length(b)), b, 'Color', 'w');
            xlabel('Component 1');
            for j = 1:size(b, 1)
                text(1, b(j), textlab{j}, 'Color', 'b');
            end
            set(gca, 'XTick', [-1 1]);
        end

        title('MDS of global image values');
    end
end

function [d, d2, g, mystd] = get_dist(v)
    % get the distances from the average
    fprintf('getting distances from mean.\t')
    gmn = mean(v, 4);
    d2 = zeros(size(v, 4), size(v, 4));

    for i = 1:size(v, 4)
        dd = abs((v(:,:,:,i) - gmn));
        d(i) = nanmean(dd(:));
        gg = v(:,:,:,i);
        gg = gg(~isnan(gg));
        g(i) = mean(gg(:));
        mystd(i) = std(gg(:));

        % get the distances from all other images
        d2(i, i) = 0;
        for j = (i+1):size(v, 4)
            dd = abs((v(:,:,:,i) - v(:,:,:,j)));
            d2(i, j) = mean(dd(~isnan(dd)));
        end
    end

    %fill in lower triangle
    d2 = d2 + d2';
    if any(isnan(g)), disp('Warning!  Some images have no valid values.'); end
end

function [c, c2, mi] = get_correl(v, do_corr, do_mi, varargin)
    % get correlations with the average and with all
    fprintf(1, 'getting correlations.\t')

    if length(varargin) > 0  %#ok  % we have a ref image instead of the grand mean
        rimg = varargin{1};
        gv = rimg(:);
    end

    if ~exist('gv', 'var') || isempty(gv)
        gmn = mean(v, 4);
        gv = gmn(:);
    end

    szv = size(v);
    v = reshape(v, prod(szv(1:3)), szv(4));

    % eliminate NaNs
    wh = union(find(isnan(gv)), find(any(isnan(v), 2)));
    if ~isempty(wh)
        gv(wh,:) = [];
        v(wh,:) = [];
    end

    c = [];
    mi = [];
    for i = 1:size(v, 2),
        if(do_corr)
            cc = corrcoef(gv, v(:,i));
            c(i) = cc(1, 2);
        end

        % mutual information
        if(do_mi)
            fprintf(1, 'MI.')
            [H, mi(i)] = hist2(gv, v(:,i), 64);
        end
    end

    c2 = corrcoef(v);
end

function [b, v, d] = pc(a, doplot)
    % a is original matrix, b is principal components, v is eigenvectors
    % (weights on columns, which = weights on voxels)

    [v, d] = eig(a);
    b = (pinv(v) * a')' ./ repmat((diag(d)').^.5, size(a, 1), 1);
    b = fliplr(b);
    v = fliplr(v);

    num = min(10, sum(diag(d) >= 1));
    b = b(:,1:num);
    v = v(:,1:num);
    origd = diag(d);
    d = diag(d)';
    d = fliplr(d);

    if doplot
        figure('Color', 'w');
        bar(real(d.^.5));
        cumprct = cumsum(d ./ sum(d));
        for i = 1:length(d)
            str = sprintf('%3.2f', cumprct(i));
            text(i-.2, d(i)^.5+.2, str);
        end
        title('Eigenvalues');
    end

    d = d(1:num);

    if num == 0
        warning('No eigenvalues above 1!');
        disp(origd);
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function f16=f16(x)

Bound=[-5 5];

if nargin==0
    f16 = Bound;
else
    f16=4*x(1,:)^2-2.1*x(1,:)^4+(1/3)*x(1,:)^6+x(1,:)*x(2,:)-4*x(2,:)^2+4*x(2,:)^4;
end

author: sfvsfv

repo: Mathematical-modeling


text:
classdef LHSintegrator_DiffReactNeumann < handle

    properties (GetAccess = public, SetAccess = private)
        M
        K
    end

    properties (Access = private)
        mesh
    end

    methods (Access = public)

        function obj = LHSintegrator_DiffReactNeumann(cParams)
            obj.mesh  = cParams.mesh;
            obj.computeStiffnessMatrix(cParams);
            obj.computeMassMatrix();
        end

        function LHS = compute(obj, epsilon)
            LHS = epsilon^2*obj.K + obj.M;
        end

    end

    methods (Access = private)

        function computeStiffnessMatrix(obj,cParams)
            s      = cParams;  % For anisotropic stiffness
            s.type = cParams.stiffType;
            s.mesh = obj.mesh;
            s.fun  = P1Function.create(obj.mesh,1);
            LHS    = LHSintegrator.create(s);
            obj.K  = LHS.compute();
        end

        function computeMassMatrix(obj)
            s.type  = 'MassMatrix';
            s.mesh  = obj.mesh;
            s.fun  = P1Function.create(obj.mesh, 1);
            s.quadratureOrder = 'QUADRATICMASS';
            LHS     = LHSintegrator.create(s);
            obj.M   = LHS.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Data File EKIN1
% Oscillations of
% a mathematical pendulum
  Ek   = '1/2*m*qt^2'; % Kinetic energy
  N    = '-(k*qt + 9.81*sin(q))*qt'; % power
  q0   = '0.1';  % initial coordinate
  qt0  = '5';    % initial velocity
  Tend = 20;     % upper bound of integration
  eps  = 1e-10;  % desirable accuracy
  np   = 2;      % number of parameters
  P{1} = 'm';    % mass of the pendulum
  P{2} = 'k';    % coefficient of resistance


text:
function [Q]=pointCloudPrinDir(V)

% function [Q]=pointCloudPrinDir(V)
% ------------------------------------------------------------------------
%
% Change log:
% 2019/04/28 Marked obsolete since function is the same as pointSetPrincipalDir
% ------------------------------------------------------------------------

[Q]=pointSetPrincipalDir(V);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function checkopt(inopt, dfopt)

% checkopt -- Check options structure for unrecognized fields
%
% Usage:
%       checkopt(inopt, dfopt)
%
% Input:
%       inopt         Input options structure
%       dfopt         Default options structure
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-06-17
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if ~isempty(inopt) && (~isfield(inopt,'NoOptionCheck') || ~inopt.NoOptionCheck),
  ifnc = fieldnames(inopt);
  for k = 1:length(ifnc),
    ifn = ifnc{k};
    if ~strcmp(ifn, 'NoOptionCheck') && ~isfield(dfopt, ifn),
      warning('SPORCO:UnknownOption', 'Unknown option field %s', ifn);
    end
  end
end

return

author: thfylsty


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

repo: sap-voicebox


text:
function [ ISIstats ] = bz_ISIStats( spikes,varargin )
%ISIstats = bz_ISIStats(spikes,varargin) calculates the statistics
%inter-spike intervals for the spiketimes in spikes.
%
%   INPUTS
%       spikes      Structure with spikes.times and spikes.UID
%
%       (options)
%       'ints'        A structure with intervals in which to calculate ISIs.
%                       states.stateNAME = [start stop]
%                       Will calculate ISIsstats separately for each state
%                       (Can also 'load' from SleepState.states.mat)
%       'cellclass'     Cell array of strings - label for each cell.
%                       (Can also 'load' from CellClass.cellinfo.mat)
%       'savecellinfo'  logical (default: false) save a cellinfo file?
%       'basePath'
%       'figfolder'     a folder to save the figure in
%       'showfig'       logical (default: false) show the figure?
%       'forceRedetect' logical (default: false) to re-compute even if saved
%       'shuffleCV2'    logical (devault: false)
%       'numISIbins'    number of bins for ISI distribution (default: 60)
%       'numCV2bins'    number pf bins for CV2 distribution (default: 60)
%
%   OUTPUTS
%       ISIstats        cellinfo structure with ISI statistics
%           .summstats  summary statistics
%           .ISIhist    histograms of ISIs etc
%           .sorts      sorting indices
%           .allspikes  ISI/CV2 value for each spike
%                       (ISI is PRECEDING interval for the spike at allspikes.times)
%
%DLevenstein 2018
%% Parse the inputs
defaultstates.ALL = [-Inf Inf];

% parse args
p = inputParser;
addParameter(p,'ints',defaultstates)
addParameter(p,'savecellinfo',false,@islogical)
addParameter(p,'basePath',pwd,@isstr)
addParameter(p,'figfolder',false)
addParameter(p,'showfig',false,@islogical);
addParameter(p,'cellclass',[]);
addParameter(p,'forceRedetect',false,@islogical);
addParameter(p,'shuffleCV2',false,@islogical);
addParameter(p,'numISIbins',60,@islogical);
addParameter(p,'numCV2bins',50,@islogical);

parse(p,varargin{:})
ints = p.Results.ints;
cellclass = p.Results.cellclass;
basePath = p.Results.basePath;
SAVECELLINFO = p.Results.savecellinfo;
figfolder = p.Results.figfolder;
SHOWFIG = p.Results.showfig;
forceRedetect = p.Results.forceRedetect;
SHUFFLECV2 = p.Results.shuffleCV2;
numISIbins = p.Results.numISIbins;
numCV2bins = p.Results.numCV2bins;

%% Load the stuff
baseName = bz_BasenameFromBasepath(basePath);
cellinfofilename = fullfile(basePath,[baseName,'.ISIStats.cellinfo.mat']);

if exist(cellinfofilename,'file') && ~forceRedetect
    ISIstats = bz_LoadCellinfo(basePath,'ISIStats');
    return
end

if strcmp(cellclass,'load')
    cellclass = bz_LoadCellinfo(basePath,'CellClass');
    cellclass = cellclass.label;
end

if strcmp(ints,'load')
    ints = bz_LoadStates(basePath,'SleepState');
    ints = ints.ints;
end

%% Get the States
ints.ALL = [-Inf Inf];  %Add ALL to the options
statenames = fieldnames(ints);
numstates = length(statenames);

%% ISI and CV2 statistics
numcells = length(spikes.times);

%Calculate ISI and CV2 for allspikes
allspikes.ISIs = cellfun(@diff,spikes.times,'UniformOutput',false);
allspikes.meanISI = cellfun(@(X) (X(1:end-1)+X(2:end))./2,allspikes.ISIs,'UniformOutput',false);
allspikes.CV2 = cellfun(@(X) 2.*abs(X(2:end)-X(1:end-1))./(X(2:end)+X(1:end-1)),allspikes.ISIs ,'UniformOutput',false);
%Make sure times line up
allspikes.times = cellfun(@(X) X(2:end-1),spikes.times,'UniformOutput',false);
allspikes.ISIs = cellfun(@(X) X(1:end-1),allspikes.ISIs,'UniformOutput',false);
%%
for ss = 1:numstates
%ss=1;

%Find which spikes are during state of interest
statespikes = cellfun(@(X) InIntervals(X,ints.(statenames{ss})),...
    allspikes.times,'UniformOutput',false);
CV2 = cellfun(@(X,Y) X(Y(2:end-1)),allspikes.CV2,statespikes,'Uniformoutput',false);
ISIs = cellfun(@(X,Y) X(Y(2:end-1)),allspikes.ISIs,statespikes,'Uniformoutput',false);
normISIs = cellfun(@(X) X./mean(X),ISIs,'Uniformoutput',false);

%Summary Statistics
summstats.(statenames{ss}).meanISI = cellfun(@(X) mean(X),ISIs);
summstats.(statenames{ss}).meanrate = 1./summstats.(statenames{ss}).meanISI;
summstats.(statenames{ss}).ISICV = cellfun(@(X) std(X)./mean(X),ISIs);
summstats.(statenames{ss}).meanCV2 = cellfun(@(X) mean(X),CV2);

%Account for no spikes
summstats.(statenames{ss}).meanrate(isnan(summstats.(statenames{ss}).meanrate))=0;

if SHUFFLECV2
    % CV2 for shuffled ISIs
    numshuffle = 100;
    for sh = 1:numshuffle
        ISIs_shuffle = cellfun(@(X) shuffle(X),ISIs,'UniformOutput',false);
        CV2_shuffle = cellfun(@(X) 2.*abs(X(2:end)-X(1:end-1))./(X(2:end)+X(1:end-1)),...
            ISIs_shuffle ,'UniformOutput',false);
        meanshuffle(sh,:) = cellfun(@(X) mean(X),CV2_shuffle);
    end
    summstats.(statenames{ss}).shufflemeanCV2 = mean(meanshuffle);
    summstats.(statenames{ss}).shufflestdCV2 = std(meanshuffle);

end
%%
%Set up all the bins and matrices
ISIhist.linbins = linspace(0,10,numISIbins);
ISIhist.logbins = linspace(log10(0.001),log10(200),numISIbins);
ISIhist.(statenames{ss}).lin = zeros(numcells,numISIbins);
ISIhist.(statenames{ss}).log = zeros(numcells,numISIbins);
ISIhist.(statenames{ss}).return = zeros(numISIbins,numISIbins,numcells);
ISIhist.(statenames{ss}).meannorm = zeros(numcells,numISIbins);

CV2hist.bins = linspace(0,2,numCV2bins+1);
CV2hist.bins = CV2hist.bins(1:end-1)+0.5.*diff(CV2hist.bins([1 2]));
CV2hist.(statenames{ss}) = zeros(numcells,numCV2bins);

Jointhist.(statenames{ss}).log = zeros(numcells,numISIbins,numCV2bins);
Jointhist.(statenames{ss}).norm = zeros(numcells,numISIbins,numCV2bins);

%Calculate all the histograms: ISI, log(ISI), 1/ISI, log(1/ISI)
for cc = 1:numcells
    numspks(cc) = length(ISIs{cc});
    spikethresh = 100; %need 100 spikes to show in mean plot
    enoughspikes(cc) = numspks(cc)>spikethresh;

    if numspks(cc)<1
        continue
    end
    %Calculate ISI histograms
    ISIhist.(statenames{ss}).lin(cc,:) = hist(ISIs{cc},ISIhist.linbins);
    ISIhist.(statenames{ss}).log(cc,:) = hist(log10(ISIs{cc}),ISIhist.logbins);
    ISIhist.(statenames{ss}).loginv(cc,:) = hist(log10(1./ISIs{cc}),ISIhist.logbins);
    ISIhist.(statenames{ss}).meannorm(cc,:) = hist(log10(normISIs{cc}),ISIhist.logbins);

    CV2hist.(statenames{ss})(cc,:) = hist(CV2{cc},CV2hist.bins);
    Jointhist.(statenames{ss}).log(cc,:,:) = hist3(...
        [log10([ISIs{cc};ISIs{cc}(2:end)]),[CV2{cc};CV2{cc}(1:end-1)]],...
        {ISIhist.logbins,CV2hist.bins});
    Jointhist.(statenames{ss}).norm(cc,:,:) = hist3(...
        [log10([normISIs{cc};normISIs{cc}(2:end)]),[CV2{cc};CV2{cc}(1:end-1)]],...
        {ISIhist.logbins,CV2hist.bins});

    %Normalize histograms to number of spikes
    ISIhist.(statenames{ss}).lin(cc,:) = ISIhist.(statenames{ss}).lin(cc,:)./numspks(cc);
    ISIhist.(statenames{ss}).log(cc,:) = ISIhist.(statenames{ss}).log(cc,:)./numspks(cc);
    ISIhist.(statenames{ss}).loginv(cc,:) = ISIhist.(statenames{ss}).loginv(cc,:)./numspks(cc);
    ISIhist.(statenames{ss}).meannorm(cc,:) = ISIhist.(statenames{ss}).meannorm(cc,:)./numspks(cc);

    CV2hist.(statenames{ss})(cc,:) = CV2hist.(statenames{ss})(cc,:)./numspks(cc);

    Jointhist.(statenames{ss}).log(cc,:,:) = Jointhist.(statenames{ss}).log(cc,:,:)./numspks(cc);
    Jointhist.(statenames{ss}).norm(cc,:,:) = Jointhist.(statenames{ss}).norm(cc,:,:)./numspks(cc);

    %Calculate Return maps
    if numspks(cc)>1
    ISIhist.(statenames{ss}).return(:,:,cc) = hist3(log10([ISIs{cc}(1:end-1) ISIs{cc}(2:end)]),{ISIhist.logbins,ISIhist.logbins});
    end
    ISIhist.(statenames{ss}).return(:,:,cc) = ISIhist.(statenames{ss}).return(:,:,cc)./numspks(cc);

end

summstats.(statenames{ss}).numspikes = numspks;

%Sortings
[~,sorts.(statenames{ss}).rate]=sort(summstats.(statenames{ss}).meanrate);
[~,sorts.(statenames{ss}).ISICV]=sort(summstats.(statenames{ss}).ISICV);
[~,sorts.(statenames{ss}).CV2]=sort(summstats.(statenames{ss}).meanCV2);

%Make the cell-type specific sortings and average distributions
if ~isempty(cellclass)
    %Check for empty cell class entries
    noclass = cellfun(@isempty,cellclass);
    sorts.numclassycells = sum(~noclass);
    %cellclass(noclass)={'none'};
    classnames = unique(cellclass(~noclass));
    numclasses = length(classnames);
    for cl = 1:numclasses
        inclasscells{cl} = strcmp(classnames{cl},cellclass);

        %Mean distributions
        meandists.(statenames{ss}).(classnames{cl}).ISIdist = squeeze(nanmean(ISIhist.(statenames{ss}).log(inclasscells{cl}&enoughspikes,:),1));
        meandists.(statenames{ss}).(classnames{cl}).CV2dist = squeeze(nanmean(CV2hist.(statenames{ss})(inclasscells{cl}&enoughspikes,:),1));
        meandists.(statenames{ss}).(classnames{cl}).Jointdist = squeeze(nanmean(Jointhist.(statenames{ss}).log(inclasscells{cl}&enoughspikes,:,:),1));

        %Sorts
        sorttypes = {'rate','ISICV','CV2'};
        for tt = 1:length(sorttypes)
        sorts.(statenames{ss}).([sorttypes{tt},classnames{cl}]) = ...
            intersect(sorts.(statenames{ss}).(sorttypes{tt}),find(inclasscells{cl}),'stable');

        if cl==1
            sorts.(statenames{ss}).([sorttypes{tt},'byclass'])=[];
        end
        sorts.(statenames{ss}).([sorttypes{tt},'byclass']) = ...
            [sorts.(statenames{ss}).([sorttypes{tt},'byclass']) sorts.(statenames{ss}).([sorttypes{tt},classnames{cl}])];
        end

    end
else
    numclasses = 1; cl=1;
    classnames{1} = 'ALL';
    sorts.numclassycells = numcells;
    inclasscells{1} = true(1,numcells);
    sorts.(statenames{ss}).ratebyclass = sorts.(statenames{ss}).rate;
    sorts.(statenames{ss}).ISICVbyclass = sorts.(statenames{ss}).ISICV;
    sorts.(statenames{ss}).CV2byclass = sorts.(statenames{ss}).CV2;

        %Mean distributions
        meandists.(statenames{ss}).(classnames{cl}).ISIdist = squeeze(nanmean(ISIhist.(statenames{ss}).log(inclasscells{cl}&enoughspikes,:),1));
        meandists.(statenames{ss}).(classnames{cl}).CV2dist = squeeze(nanmean(CV2hist.(statenames{ss})(inclasscells{cl},:)&enoughspikes,1));
        meandists.(statenames{ss}).(classnames{cl}).Jointdist = squeeze(nanmean(Jointhist.(statenames{ss}).log(inclasscells{cl}&enoughspikes,:,:),1));
end

%% ACG - get its own analysis
%[ccg,t] = CCG(statespiketimes,[],<options>)

%%
if SHOWFIG | figfolder
figure
    subplot(2,2,1)
        for cl = 1:numclasses
            plot(log10(summstats.(statenames{ss}).meanrate(inclasscells{cl})),...
                log2(summstats.(statenames{ss}).ISICV(inclasscells{cl})),'.','markersize',11)
            hold on
        end
        plot(get(gca,'xlim'),log2([1 1]),'k')
        LogScale('x',10);LogScale('y',2);
        xlabel('Mean Rate (Hz)');ylabel('ISI CV')
        title(statenames{ss})
        box off

    subplot(2,2,2)
    for cl = 1:numclasses

        if SHUFFLECV2
            plot([log10(summstats.(statenames{ss}).meanrate(inclasscells{cl}));...
                log10(summstats.(statenames{ss}).meanrate(inclasscells{cl}))],...
                [summstats.(statenames{ss}).shufflemeanCV2(inclasscells{cl});...
                summstats.(statenames{ss}).meanCV2(inclasscells{cl})],...
                'color',0.8.*[1 1 1],'linewidth',0.25)
            hold on
            plot([log10(summstats.(statenames{ss}).meanrate(inclasscells{cl}));...
                log10(summstats.(statenames{ss}).meanrate(inclasscells{cl}))],...
                [summstats.(statenames{ss}).shufflemeanCV2(inclasscells{cl})-summstats.(statenames{ss}).shufflestdCV2(inclasscells{cl});...
                summstats.(statenames{ss}).shufflemeanCV2(inclasscells{cl})+summstats.(statenames{ss}).shufflestdCV2(inclasscells{cl})],...
                'color',0.8.*[1 1 1],'linewidth',2.5)
        end

        plot(log10(summstats.(statenames{ss}).meanrate(inclasscells{cl})),...
            (summstats.(statenames{ss}).meanCV2(inclasscells{cl})),'.','markersize',11)
        hold on

    end
        plot(get(gca,'xlim'),[1 1],'k')
        LogScale('x',10);
        xlabel('Mean Rate (Hz)');ylabel('ISI <CV2>')
        title(statenames{ss})
        box off

    subplot(2,3,4)
        imagesc((ISIhist.logbins),[1 sorts.numclassycells],...
            ISIhist.(statenames{ss}).log(sorts.(statenames{ss}).ratebyclass,:))
        hold on
        plot(log10(1./(summstats.(statenames{ss}).meanrate(sorts.(statenames{ss}).ratebyclass))),[1:sorts.numclassycells],'k.','LineWidth',2)
        plot(ISIhist.logbins([1 end]),sum(inclasscells{1}).*[1 1]+0.5,'r')
        LogScale('x',10)
        xlabel('ISI (s)')
        xlim(ISIhist.logbins([1 end]))
        %colorbar
      %  legend('1/Mean Firing Rate (s)','location','southeast')
        ylabel('Cell (Sorted by FR, Type)')
        %legend('1/Mean Firing Rate (s)','location','southeast')
        caxis([0 0.1])
        title('ISI Distribution (Log Scale)')

    subplot(2,3,5)
        imagesc((CV2hist.bins),[1 sorts.numclassycells],...
            CV2hist.(statenames{ss})(sorts.(statenames{ss}).ratebyclass,:))
        hold on
        plot((summstats.(statenames{ss}).meanCV2(sorts.(statenames{ss}).ratebyclass)),[1:sorts.numclassycells],'k.','LineWidth',2)
        plot(CV2hist.bins([1 end]),sum(inclasscells{1}).*[1 1]+0.5,'r')
        %LogScale('x',10)
        xlabel('CV2')
        xlim(CV2hist.bins([1 end]))
        %colorbar
        xlim([0 2])
      %  legend('1/Mean Firing Rate (s)','location','southeast')
        ylabel('Cell (Sorted by FR, Type)')
        %legend('1/Mean Firing Rate (s)','location','southeast')
        %caxis([0 0.1])
        title('CV2 Distribution')

	for cl = 1:numclasses
        subplot((numclasses.*2),3,cl.*3+(numclasses.*3))
            imagesc(ISIhist.logbins,CV2hist.bins,meandists.(statenames{ss}).(classnames{cl}).Jointdist')
            hold on
            plot(ISIhist.logbins,meandists.(statenames{ss}).(classnames{cl}).ISIdist.*20,'k')
            plot(meandists.(statenames{ss}).(classnames{cl}).CV2dist.*20-3,CV2hist.bins,'k')
            axis xy
            LogScale('x',10)
            ylabel({(classnames{cl}),'CV2'});
            xlabel('ISI (s)');

            bottomrow = meandists.(statenames{ss}).(classnames{cl}).Jointdist(:,1);
            caxis([0 max([bottomrow;0])])
	end

%     subplot(2,3,5)
%         imagesc((ISIhist.logbins),[1 sorts.numclassycells],...
%             ISIhist.(statenames{ss}).log(sorts.(statenames{ss}).ISICVbyclass,:))
%         hold on
%         plot(log10(1./(summstats.(statenames{ss}).meanrate(sorts.(statenames{ss}).ISICVbyclass))),[1:sorts.numclassycells],'k.','LineWidth',2)
%         plot(ISIhist.logbins([1 end]),sum(inclasscells{1}).*[1 1]+0.5,'r')
%         LogScale('x',10)
%         xlabel('ISI (s)')
%         xlim(ISIhist.logbins([1 end]))
%         %colorbar
%       %  legend('1/Mean Firing Rate (s)','location','southeast')
%         ylabel('Cell (Sorted by CV, Type)')
%         %legend('1/Mean Firing Rate (s)','location','southeast')
%         caxis([0 0.1])
%         title('ISI Distribution (Log Scale)')
%
%     subplot(2,3,6)
%         imagesc((ISIhist.logbins),[1 sorts.numclassycells],...
%             ISIhist.(statenames{ss}).log(sorts.(statenames{ss}).CV2byclass,:))
%         hold on
%         plot(log10(1./(summstats.(statenames{ss}).meanrate(sorts.(statenames{ss}).CV2byclass))),[1:sorts.numclassycells],'k.','LineWidth',2)
%         plot(ISIhist.logbins([1 end]),sum(inclasscells{1}).*[1 1]+0.5,'r')
%         LogScale('x',10)
%         xlabel('ISI (s)')
%         xlim(ISIhist.logbins([1 end]))
%         %colorbar
%       %  legend('1/Mean Firing Rate (s)','location','southeast')
%         ylabel('Cell (Sorted by CV2, Type)')
%         %legend('1/Mean Firing Rate (s)','location','southeast')
%         caxis([0 0.1])
%         title('ISI Distribution (Log Scale)')

%     subplot(2,2,3)
%         plot(log2(summstats.(statenames{ss}).meanrate(CellClass.pE)),...
%             log2(summstats.(statenames{ss}).meanCV2(CellClass.pE)),'k.')
%         hold on
%         plot(log2(summstats.(statenames{ss}).meanrate(CellClass.pI)),...
%             log2(summstats.(statenames{ss}).meanCV2(CellClass.pI)),'r.')
%         LogScale('xy',2)
%         xlabel('Mean Rate (Hz)');ylabel('Mean CV2')
%         title(statenames{ss})
%         box off

if figfolder
    NiceSave(['ISIstats_',(statenames{ss})],figfolder,baseName);
end

%exneurons %top/bottom 25 %ile rate/ISICV
%%
% exwindur = 4; %s
% STATEtimepoints = Restrict(lfp.timestamps,double(SleepState.ints.(statenames{ss})));
% samplewin = STATEtimepoints(randi(length(STATEtimepoints))) + [0 exwindur];
% %%
% figure
% bz_MultiLFPPlot( lfp,'spikes',spikes,'timewin',samplewin,...
%     'sortmetric',summstats.(statenames{ss}).meanCV2,...
%     'cellgroups',{CellClass.pI,CellClass.pE})

end

ISIstats.summstats = summstats;
ISIstats.ISIhist = ISIhist;
ISIstats.CV2hist = CV2hist;
ISIstats.Jointhist = Jointhist;
ISIstats.sorts = sorts;
ISIstats.UID = spikes.UID;
ISIstats.allspikes = allspikes;

if SAVECELLINFO
    save(cellinfofilename,'ISIstats')
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ ret ] =setPin1Off( t_Kuka )
%% About:
% This function is used to set off the Pin 1 of the KUKA iiwa 7 R 800.

%% Syntax:
% [ ret ] = setPin1Off( t_Kuka )

%% Arreguments:
% t_Kuka: is the TCP/IP connection

%% Return value:
% ret: a boolean variable
% true: if the message is treated successfully
% false: if an error has happened

% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='pin1off';
fprintf(t_Kuka, theCommand);
message=fgets(t_Kuka);
ret=checkAcknowledgment(message);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_snr_40db_s_16.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_snr_40_db_s_16';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'MMV, SNR=40dB s=16';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function inspect_issue1674

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_prepare_layout ft_convert_coordsys

%%

elec = [];
elec.label = {
  'Fpz'
  'T7'
  'Cz'
  'T8'
  'Oz'
  };
elec.elecpos = [
   80   0   0 % a
    0  80   0 % l
    0   0  80 % s
    0 -80   0
  -80   0   0
  ];
elec.unit = 'mm';
elec.coordsys = 'eeglab';

%%

% this was the key issue that was reported, Fpz should be at the front of the head

cfg = [];
cfg.elec = elec;
layout = ft_prepare_layout(cfg);
ft_plot_layout(layout);

%%

% this was a side issue that I detected while testing, since 'eeglab' is also 'als'
% the conversion from 'als' to 'ras' (and others) was not working as expected

elec = [];
elec.label = {
  'a'
  'l'
  's'
  };
elec.elecpos = [
  1 0 0 % a
  0 1 0 % l
  0 0 1 % s
  ];
elec.unit = 'mm';
elec.coordsys = 'als';

target = {'als' 'ali' 'ars' 'ari' 'pls' 'pli' 'prs' 'pri' 'las' 'lai' 'ras' 'rai' 'lps' 'lpi' 'rps' 'rpi' 'asl' 'ail' 'asr' 'air' 'psl' 'pil' 'psr' 'pir' 'sal' 'ial' 'sar' 'iar' 'spl' 'ipl' 'spr' 'ipr' 'sla' 'ila' 'sra' 'ira' 'slp' 'ilp' 'srp' 'irp' 'lsa' 'lia' 'rsa' 'ria' 'lsp' 'lip' 'rsp' 'rip'};

elec = {elec};

for i=2:length(target)
  disp(target{i});

  elec{i} = ft_convert_coordsys(elec{i-1}, target{i});

  % the 1st electrode is at the nose
  letter1indx = find(elec{i}.elecpos(1,:)~=0); % where did it go?
  if elec{i}.elecpos(1,letter1indx)==1         % is the axis flipped?
    letter1val = 'a';
  else
    letter1val = 'p';
  end
  % the 2nd electrode is at the left
  letter2indx = find(elec{i}.elecpos(2,:)~=0);
  if elec{i}.elecpos(2,letter2indx)==1
    letter2val = 'l';
  else
    letter2val = 'r';
  end
  % the 3rd electrode is at the top
  letter3indx = find(elec{i}.elecpos(3,:)~=0);
  if elec{i}.elecpos(3,letter3indx)==1
    letter3val = 's';
  else
    letter3val = 'i';
  end

  result = 'xxx';
  result(letter1indx) = letter1val;
  result(letter2indx) = letter2val;
  result(letter3indx) = letter3val;

  assert(isequal(result, elec{i}.coordsys));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [node,elem,bdFlag,HB] = uniformbisect3(node,elem,bdFlag,HB)
%% UNIFORMBISECT3 uniformly bisect a 3-D triangulation.
%
% [node,elem] = uniformbisect3(node,elem) divides each tetrahedron into 8 small
% tetrahedrons using longest vertex bisection.
%
% [node,elem,bdFlag,HB] = uniformbisect3(node,elem,bdFlag,HB) update HB and
% boundary conditions.
%
% Example
%
%     node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1];
%     elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7];
%     elem = label3(node,elem);
%     figure(1); subplot(1,3,1);
%     set(gcf,'Units','normal'); set(gcf,'Position',[0.25,0.25,0.5,0.3]);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     [node,elem] = uniformbisect3(node,elem);
%     figure(1); subplot(1,3,2);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%     bdFlag = setboundary3(node,elem,'Dirichlet');
%     [node,elem,bdFlag] = uniformbisect3(node,elem,bdFlag);
%     figure(1); subplot(1,3,3);
%     showmesh3(node,elem,[],'FaceAlpha',0.35); view([210 8]);
%
% See also uniformbisect, uniformrefine, bisect, bisect3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('bdFlag','var')
    bdFlag =[];
end
if ~exist('HB','var')
    HB = [];
end
[node,elem,bdFlag,HB] = bisect3(node,elem,'all',bdFlag,HB);
[node,elem,bdFlag,HB] = bisect3(node,elem,'all',bdFlag,HB);
[node,elem,bdFlag,HB] = bisect3(node,elem,'all',bdFlag,HB);

repo: ifem


text:
classdef RankTwoLaminateHomogenizerFromAllaireTestedWebPageCode

    properties
    end

    methods (Static)
        function Chomog = ChinvOld(lambda,mu,d1,d2,m1,m2,theta,epsil)

            sqrt2 = sqrt(2);

            e1x = d1(1);
            e1y = d1(2);

            e2x = d2(1);
            e2y = d2(2);

            K = (mu+lambda)/(mu*(2*mu+lambda));

            muVoigt = mu;

            A =	m1 * ( (lambda+2*mu) - 1/mu*(lambda^2*e1y^2+(lambda+2*mu)^2*e1x^2) + K*((lambda+2*mu)*e1x^2+lambda*e1y^2)^2 ) + ...
                +m2* ( (lambda+2*mu) - 1/mu*(lambda^2*e2y^2+(lambda+2*mu)^2*e2x^2) + K*((lambda+2*mu)*e2x^2+lambda*e2y^2)^2 );
            B =	m1 * ( (lambda+2*mu) - 1/mu*(lambda^2*e1x^2+(lambda+2*mu)^2*e1y^2) + K*((lambda+2*mu)*e1y^2+lambda*e1x^2)^2 )  + ...
                +m2* ( (lambda+2*mu) - 1/mu*(lambda^2*e2x^2+(lambda+2*mu)^2*e2y^2) + K*((lambda+2*mu)*e2y^2+lambda*e2x^2)^2 );
            C =	m1 * ( 4*mu - 1/mu*(2*mu)^2 + K*(4*mu*e1x*e1y)^2 )/2  + ...
                +m2* ( 4*mu - 1/mu*(2*mu)^2 + K*(4*mu*e2x*e2y)^2 )/2;
            D = 	m1 * ( 2*lambda - 1/mu*(2*lambda*(lambda+2*mu)) + K*2*((lambda+2*mu)*e1y^2+lambda*e1x^2)*((lambda+2*mu)*e1x^2+lambda*e1y^2) )/2  + ...
                +m2* ( 2*lambda - 1/mu*(2*lambda*(lambda+2*mu)) + K*2*((lambda+2*mu)*e2y^2+lambda*e2x^2)*((lambda+2*mu)*e2x^2+lambda*e2y^2) )/2;
            E =	m1 * ( -1/mu*(4*mu*(2*lambda+2*mu)*e1x*e1y) + K*2*(4*mu*e1x*e1y*((lambda+2*mu)*e1y^2+lambda*e1x^2)) )/(2*sqrt2)  + ...
                +m2* ( -1/mu*(4*mu*(2*lambda+2*mu)*e2x*e2y) + K*2*(4*mu*e2x*e2y*((lambda+2*mu)*e2y^2+lambda*e2x^2)) )/(2*sqrt2);
            F = 	m1 * ( -1/mu*(4*mu*(2*lambda+2*mu)*e1x*e1y) + K*2*(4*mu*e1x*e1y*((lambda+2*mu)*e1x^2+lambda*e1y^2)) )/(2*sqrt2)  + ...
                +m2* ( -1/mu*(4*mu*(2*lambda+2*mu)*e2x*e2y) + K*2*(4*mu*e2x*e2y*((lambda+2*mu)*e2x^2+lambda*e2y^2)) )/(2*sqrt2);

            %// Ajout du matériau mou (qui simule le vide)
            A = epsil/(1.-epsil)*(lambda+2*mu)	+ theta*A;
            B = epsil/(1.-epsil)*(lambda+2*mu)	+ theta*B;
            C = epsil/(1.-epsil)*muVoigt		+ theta*C;
            D = epsil/(1.-epsil)*lambda		+ theta*D;
            E = epsil/(1.-epsil)*0			+ theta*E;
            F = epsil/(1.-epsil)*0			+ theta*F;

            %//Première inversion
            DET = A*B*C-A*E*E-B*F*F-C*D*D+2*D*E*F;
            A1=(B*C-E*E)/DET;
            B1=(A*C-F*F)/DET;
            C1=(A*B-D*D)/DET;
            D1=(E*F-C*D)/DET;
            E1=(D*F-A*E)/DET;
            F1=(D*E-B*F)/DET;

            A = (1-theta)*A1+(lambda+2*mu)/(4*mu*(lambda+mu));
            B = (1-theta)*B1+(lambda+2*mu)/(4*mu*(lambda+mu));
            C = (1-theta)*C1+1/(muVoigt);
            D = (1-theta)*D1-lambda/(4*mu*(lambda+mu));
            E = (1-theta)*E1;
            F = (1-theta)*F1;

            %//Deuxième inversion
            DET = A*B*C-A*E*E-B*F*F-C*D*D+2*D*E*F;
            A1=(B*C-E*E)/DET;
            B1=(A*C-F*F)/DET;
            C1=(A*B-D*D)/DET;
            D1=(E*F-C*D)/DET;
            E1=(D*F-A*E)/DET;
            F1=(D*E-B*F)/DET;
            A=A1;
            B=B1;
            C=C1;
            D=D1;
            E=E1;
            F=F1;

            Chomog = [A   D	 F;
                D   B  E;
                F   E  C];

        end

    end

end


text:
function value=get(obj, key)
% Get test property
%
% Inputs:
%   obj                 MOxUnitFieldTripTestCase instance
%   key                 test property key, one of:
%                       'mem'       Maximum memory the test uses (bytes)
%                       'walltime'  Maximum time the test takes (seconds)
%                       'dccnpath'  Set to true if this function requires
%                                   dccnpath
%
% Output:
%   value               The value associated with key, or empty if the
%                       property is not set for this test case
%
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #

    properties=obj.properties;

    if ~isfield(properties,key)
        value=[];
        return;
    end

    value=properties.(key);


text:
function product=prod(X)
%PROD (overloaded)

product = 1;
for i = 1:length(X)
   pick = cell(1,1);pick{1}={i};
   product = product*subsref(X,struct('type','()','subs',pick));
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP