text:
function drosPlot(model, totgenes, thisgene)

% DROSPLOT Plot a GPDISIM model/models
% FORMAT
% DESC Plot a GPDISIM model/models
% ARG model : The model to plot
% ARG totgenes : Total number of genes if plotting several single
% target models in a single figure (optional)
% ARG thisgene : Index of this gene in the above case (optional)
%
% COPYRIGHT : Antti Honkela, 2009

% SHEFFIELDML

FONTSIZE = 8;
LINEWIDTH = 1;
MARKERSIZE = 6;

if strcmp(model.type, 'cgpdisim'),
  numGenes = model.comp{1}.numGenes + 1;
else
  numGenes = model.comp{1}.numGenes;
end

selectgenes = 1:numGenes;

numPlots = length(selectgenes);

tf = model.comp{1}.annotation.tf;
genes = model.comp{1}.annotation.genes;
targets = genes(2:end);

genenames = genes;
times = model.comp{1}.t;

y = {};
yvar = {};
for k=1:length(model.comp),
  y{k} = reshape(model.comp{k}.y, [length(times), length(genes)]);
  yvar{k} = reshape(model.comp{k}.yvar, [length(times), length(genes)]);
end

for j = 1:length(model.comp)

  % Generate predictions of the functions.
  % to do this we need to compute the K_xf portions of the kernel
  % (simXrbfKernCompute does this for us).
  predt = [1:0.1:12 model.comp{j}.t']';

  if model.comp{j}.includeNoise,
    simMultiKern = model.comp{j}.kern.comp{1};
  else
    simMultiKern = model.comp{j}.kern;
  end
  if strcmp(model.type, 'cgpdisim'),
    proteinKern = kernCreate(model.comp{j}.t, 'sim');
    proteinKern.inverseWidth = simMultiKern.comp{1}.inverseWidth;
    proteinKern.decay = model.comp{j}.delta;
    proteinKern.variance = simMultiKern.comp{2}.di_variance;
    % simXrbf requires an RBF with unit variance
    inputKern = kernCreate(model.comp{j}.t, 'rbf');
    inputKern.inverseWidth = simMultiKern.comp{1}.inverseWidth;
    inputKern.variance = 1;
    K = simXrbfKernCompute(proteinKern, inputKern, ...
			   predt, model.comp{j}.t);
    K = K * simMultiKern.comp{1}.variance;
    warning('off', 'KERN:simRBFVariance');
    for i=2:simMultiKern.numBlocks
      blockK =  disimXsimKernCompute(simMultiKern.comp{i}, proteinKern, ...
				     model.comp{j}.t, predt);
      K = [K blockK'];
    end
    warning('on', 'KERN:simRBFVariance');

    ymean = reshape(ones(length(times),1)*[0, model.comp{j}.mu], ...
		    size(model.comp{j}.y));
  else
    proteinKern = kernCreate(model.comp{1}.t, 'rbf');
    proteinKern.inverseWidth = ...
	simMultiKern.comp{1}.inverseWidth;
    K = [];

    for i=1:simMultiKern.numBlocks
      K = [K; simXrbfKernCompute(simMultiKern.comp{i}, proteinKern, ...
				 model.comp{j}.t, predt)];
    end

    K = K';
    ymean = reshape(ones(length(times),1)*model.comp{j}.mu, ...
		    size(model.comp{j}.y));
  end

  predF = K*model.comp{j}.invK*(model.comp{j}.y-ymean);
  if strcmp(model.type, 'cgpdisim'),
    varF = simMultiKern.comp{1}.variance * ...
	   kernDiagCompute(proteinKern, predt) - sum(K'.* ...
						     (model.comp{j}.invK* ...
						      K'), 1)';
  else
    varF = kernDiagCompute(proteinKern, predt) - sum(K'.* ...
						     (model.comp{j}.invK* ...
						      K'), 1)';
  end

  % Predicted Gene Expressions
  if strcmp(model.type, 'cgpdisim'),
    Kxx = multiKernCompute(simMultiKern, predt, model.comp{j}.t);
    varX = real(kernDiagCompute(simMultiKern, predt) - ...
		sum(Kxx'.* (model.comp{j}.invK*Kxx'), 1)');
    meanPredX = reshape(ones(length(predt),1)*([0 model.comp{j}.B./ ...
		    model.comp{j}.D]), length(predt)*(numGenes), 1);
  else
    Kxx = multiKernCompute(simMultiKern, predt, model.comp{j}.t);
    varX = real(kernDiagCompute(simMultiKern, predt) - ...
		sum(Kxx'.* (model.comp{j}.invK*Kxx'), 1)');
    meanPredX = reshape(ones(length(predt),1)*(model.comp{j}.B./ ...
		    model.comp{j}.D), length(predt)*(numGenes), 1);
  end

  predX = meanPredX + real(Kxx*model.comp{j}.invK*(model.comp{j}.y-ymean));

  % Take out predictions at data points.
  % Use them to get the scale for the other data.
  numData = length(predX)/(numGenes);
  predExprs = reshape(predX, numData, numGenes);
  meanExprs = ones(numData, 1)*mean(predExprs);
  scaleExprs = ones(numData, 1)*(sqrt(var(predExprs))./sqrt(var(y{j})));

  predExprs(end-length(model.comp{j}.t)+1:end,:) = [];
  varExprs = reshape(varX, numData, numGenes);
  varExprs = varExprs./scaleExprs./scaleExprs;
  varExprs(end-length(model.comp{j}.t)+1:end,:) = [];
  predF(end-length(model.comp{j}.t)+1:end) = [];
  varF(end-length(model.comp{j}.t)+1:end) = [];
  predt(end-length(model.comp{j}.t)+1:end) = [];

  yscale = ones(length(times), 1)*sqrt(var(y{j}));
  ymean = ones(length(times),1)*mean(y{j});
  ynormal = ymean + (y{j}-ymean)./yscale;
  yvarNormal = yvar{j}./yscale./yscale;

  scalePred = sqrt(var(predExprs));

  if nargin < 2,
    figure;
    if strcmp(model.type, 'cgpdisim'),
      numRows = 3;
      numCols = numGenes - 1;
      inputCol = floor((numCols+1)/2);

      subplot(numRows, numCols, inputCol + numCols);
    else
      numRows = 2;
      numCols = numGenes;
      inputCol = floor((numCols+1)/2);

      subplot(numRows, numCols, inputCol);
    end
  else
    figure(j);
    subplot(numPlots+1, totgenes, totgenes + thisgene);
  end
  bh = fill([predt', fliplr(predt')], [predF' + 2*sqrt(varF'), fliplr(predF' - 2*sqrt(varF'))], [.8 .8 .8], 'EdgeColor', [.8, .8, .8]);
  hold on,
  lin = plot(predt, predF, '-');
  hold off
  title(sprintf('Inferred %s protein', tf), 'fontsize', FONTSIZE);
  set(bh, 'lineWidth', LINEWIDTH);
  set(lin, 'lineWidth', 1.5*LINEWIDTH);
  %set(lin, 'markersize', 20);
  set(gca, 'fontname', 'arial', 'fontsize', FONTSIZE, 'xlim', [min(model.comp{j}.t) ...
                                                    max(model.comp{j} ...
                                                    .t)]);
    set(gca, 'YTick', []);
    set(gca, 'XTickLabel', []);
    axis tight
    v = axis;
    axis([v(1)-.05*(v(2)-v(1)), v(2)+.05*(v(2)-v(1)), ...
	  v(3)-.05*(v(4)-v(3)), v(4)+.05*(v(4)-v(3))])
%  set(gca, 'ylim', [-0.1 0.4]);

  for index = 1:numPlots
    if nargin >= 2,
      if index==1,
	subplot(numPlots+1, totgenes, thisgene);
      else
	subplot(numPlots+1, totgenes, index*totgenes+thisgene);
      end
    else
      if strcmp(model.type, 'cgpdisim'),
	if index == 1,
	  subplot(numRows, numCols, inputCol);
	else
	  subplot(numRows, numCols, (numRows-1)*numCols + index - 1);
	end
      else
	subplot(numRows, numCols, (numRows-1)*numCols + index);
      end
    end
    bh = fill([predt', fliplr(predt')], [predExprs(:,selectgenes(index))' + 2*sqrt(varExprs(:,selectgenes(index))'), fliplr(predExprs(:,selectgenes(index))' - 2*sqrt(varExprs(:,selectgenes(index))'))], [.8, .8, .8], 'EdgeColor', [.8, .8, .8]);
    hold on,
    lin = plot(predt, predExprs(:,selectgenes(index)), '-');
    lin = [lin plot(times, y{j}(:,selectgenes(index)), 'rx')];
    lin1 = errorbar(times, y{j}(:,selectgenes(index)), 2*sqrt(yvar{j}(:,selectgenes(index))), ...
                    'rx');
    hold off
    if strcmp(model.type, 'cgpdisim') && selectgenes(index) == 1,
      titleText = sprintf('%s mRNA (input)', tf);
%      lin = [lin plot(predt, predDi, 'm-')];
    else
      titleText = sprintf('%s mRNA', genenames{selectgenes(index)});
    end
    title(titleText, 'fontsize', FONTSIZE);
    set(bh, 'lineWidth', LINEWIDTH);
    set(lin, 'lineWidth', 1.5*LINEWIDTH);
    set(lin, 'markersize', MARKERSIZE);
    set(lin1, 'lineWidth', LINEWIDTH);
    set(gca, 'fontname', 'arial', 'fontsize', FONTSIZE, 'xlim', [min(predt) ...
                        max(predt)]);
    set(gca, 'YTick', []);
    if nargin < 2,
      if strcmp(model.type, 'cgpdisim') && index == 1,
	set(gca, 'XTickLabel', []);
      else
	xlabel('Time (h)')
      end
    else
      if index ~= numPlots,
	set(gca, 'XTickLabel', []);
      else
	xlabel('Time (h)')
      end
    end
    axis tight
    v = axis;
    axis([v(1)-.05*(v(2)-v(1)), v(2)+.05*(v(2)-v(1)), ...
	  v(3)-.05*(v(4)-v(3)), v(4)+.05*(v(4)-v(3))])
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function runningOctave = isOctave()

% Returns whether this code has been run using Octave or MATLAB

try
    OCTAVE_VERSION;
    runningOctave = 1;
catch
    runningOctave = 0;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Simulating Transducer Field Patterns Example
%
% This example demonstrates the use of k-Wave to compute the field pattern
% generated by a curved single element transducer in two dimensions. It
% builds on the Monopole Point Source In A Homogeneous Propagation Medium
% Example.
%
% author: Bradley Treeby
% date: 10th December 2009
% last update: 24th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 216;           % number of grid points in the x (row) direction
Ny = 216;           % number of grid points in the y (column) direction
dx = 50e-3/Nx;    	% grid point spacing in the x direction [m]
dy = dx;            % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_power = 1.5;   % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

% define a curved transducer element
source.p_mask = makeCircle(Nx, Ny, 61, 61, 60, pi/2);

% define a time varying sinusoidal source
source_freq = 0.25e6;       % [Hz]
source_mag = 0.5;           % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);

% filter the source to remove high frequencies not supported by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);

% create a display mask to display the transducer
display_mask = source.p_mask;

% create a sensor mask covering the entire computational domain using the
% opposing corners of a rectangle
sensor.mask = [1, 1, Nx, Ny].';

% set the record mode capture the final wave-field and the statistics at
% each sensor point
sensor.record = {'p_final', 'p_max', 'p_rms'};

% assign the input options
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

% add the source mask onto the recorded wave-field
sensor_data.p_final(source.p_mask ~= 0) = 1;
sensor_data.p_max(source.p_mask ~= 0) = 1;
sensor_data.p_rms(source.p_mask ~= 0) = 1;

% plot the final wave-field
figure;
subplot(1, 3, 1), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_final, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Final Wave Field');

% plot the maximum recorded pressure
subplot(1, 3, 2), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_max, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('Maximum Pressure');

% plot the rms recorded pressure
subplot(1, 3, 3), imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data.p_rms, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
title('RMS Pressure');
scaleFig(2, 1);

repo: Aero-matlab


text:
function dd=tritrig(a,b,c,A,B,C)
if A==90 | B==90
    'error, right angle can only be C'

end
x=0;
if C==90
    if a~=0 & b~=0
        c=sqrt(a^2+b^2);
        A=asind(a/c);
        B=asind(b/c);
        x=1;

    end
    if a~=0 & c~=0 & x==0
        b=sqrt(c^2-a^2);
        A=asind(a/c);
        B=asind(b/c);
        x=1;

    end
    if b~=0 & c~=0 & x==0
        a=sqrt(c^2-b^2);
        A=asind(a/c);
        B=asind(b/c);
        x=1;

    end
    if A~=0 & x==0
        B=90-A;
        if a~=0
            b=a/tand(A);
            c=a/sind(a);
            x=1;

        end
        if b~=0 & x==0
            a=tand(A)*b;
            c=b/cosd(A);
            x=1;

        end
        if c~=0 & x==0
            a=sind(A)*c;
            b=cosd(A)*c;
            x=1;

        end
    end
    if B~=0 & x==0
        A=90-B;
        if a~=0
            b=a/tand(A);
            c=a/sind(a);
            x=1;

        end
        if b~=0 & x==0
            a=tand(A)*b;
            c=b/cosd(A);
            x=1;

        end
        if c~=0 & x==0
            a=sind(A)*c;
            b=cosd(A)*c;
            x=1;

        end
    end
end

if a~=0 & b~=0 & c~=0 & x==0
    A=acosd((b^2+c^2-a^2)/2/b/c);
    B=asind(sind(A)*b/a);
    C=180-A-B;
    x=1;

end
if a~=0 & b~=0 & C~=0 & x==0
    c=sqrt(a^2+b^2-2*a*b*cosd(C));
    A=asind(sind(B)*a/b);
    B=asind(sind(A)*b/a);
    x=1;

end
if a~=0 & c~=0 & B~=0 & x==0
    b=sqrt(a^2+c^2-2*a*c*cosd(B));
    A=asind(sind(B)*a/b);
    C=180-B-A;
    x=1;

end
if b~=0 & c~=0 & A~=0 & x==0
    a=sqrt(b^2+c^2-2*b*c*cosd(A));
    B=asind(sind(A)*b/a);
    C=180-A-B;
    x=1;

end
if A~=0 & B~=0 & x==0
    C=180-A-B;
    if a~=0
        b=sind(B)*a/sind(A);
        c=sind(C)*a/sind(A);
        x=1;

    end
    if b~=0 & x==0
        a=sind(A)*b/sind(B);
        c=sind(C)*b/sind(B);
        x=1;

    end
    if c~=0 & x==0
        a=sind(A)*c/sind(C);
        b=sind(B)*c/sind(C);
        x=1;

    end
end
if A~=0 & C~=0 & x==0
    B=180-A-C;
    if a~=0
        b=sind(B)*a/sind(A);
        c=sind(C)*a/sind(A);
        x=1;

    end
    if b~=0 & x==0
        a=sind(A)*b/sind(B);
        c=sind(C)*b/sind(B);
        x=1;

    end
    if c~=0 & x==0
        a=sind(A)*c/sind(C);
        b=sind(B)*c/sind(C);
        x=1;

    end
end
if B~=0 & C~=0 & x==0
    A=180-B-C;
    if a~=0
        b=sind(B)*a/sind(A);
        c=sind(C)*a/sind(A);
        x=1;

    end
    if b~=0 & x==0
        a=sind(A)*b/sind(B);
        c=sind(C)*b/sind(B);
        x=1;

    end
    if c~=0 & x==0
        a=sind(A)*c/sind(C);
        b=sind(B)*c/sind(C);
        x=1;

    end
end
if x==0
    'not enough information was given'
else
a
b
c
A
B
C
end
%How to use tritrig:
%tritrig is in the form 'tritrig(a,b,c,A,B,C)' where a b c A B C are the
%angles and the lengths of a triangle. If your triangle has a right angle
%in it, it must be 'C' and the hypotenuse must be 'c'. The angles and
%lengths that you do not know must be typed as '0'
%Eg.
%tritrig(3,4,5,0,0,0)
%A =
%    36.8699
%B=
%    53.1301
%C=
%    90.0000

author: Sable

repo: mcbench-benchmarks


text:
%% DML Toolbox Functions
%
% * <analysis.html |dml.analysis|> - multivariate analysis class
% * <blogreg.html |dml.blogreg|> - Bayesian logistic regression
% * <bootstrap.html |dml.bootstrap|> - bootstrapping to determine parameter relevance
% * <circreg.html |dml.circreg|> - circular regression method
% * <corclas.html |dml.corclas|> - template matching correlation based classifier
% * <crossvalidator.html |dml.crossvalidator|> - crossvalidation class
% * <enet.html |dml.enet|> - efficient elastic net algorithm
% * <filterer.html |dml.filterer|> - filtering approach to feature selection
% * <garrote.html |dml.garrote|> - variational garrote
% * <gp.html |dml.gp|> - Gaussian process
% * <graphnet.html |dml.graphnet|> - native implementation of graphnet algorithm
% * <gridsearch.html |dml.gridsearch|> - grid search method
% * <hmm.html |dml.hmm|> - Hidden Markov model with continuous observations
% * <lds.html |dml.lds|> - linear dynamical system
% * <method.html |dml.method|> - abstract class for multivariate methods
% * <naive.html |dml.naive|> - gaussian naive Bayes classifier
% * <ndata.html |dml.ndata|> - wrapper class to make methods handle multiple datas
% * <noutput.html |dml.noutput|> - wrapper class to make methods handle multiple outpu
% * <one_against_one.html |dml.one_against_one|> - one-against-one binary classification
% * <one_against_rest.html |dml.one_against_rest|> - one-against-rest binary classification
% * <permutation.html |dml.permutation|> - permutation testing class
% * <prior.html |dml.prior|> - used to created smoothing and shrinkage priors
% * <searchlight.html |dml.searchlight|> - searchlight analysis
% * <sincos.html |dml.sincos|> - circular regression by decomposing into sine and co
% * <slda.html |dml.slda|> - shrinkage linear discriminant analysis
% * <sopls.html |dml.sopls|> - sparse orthonormalized partial least squares
% * <standardizer.html |dml.standardizer|> - takes zscores
% * <statistic.html |dml.statistic|> - test statistics
% * <svm.html |dml.svm|> - support vector machine
% * <whiten.html |dml.whiten|> - whitens the data

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cdot(x,y,z,c,dot_size,label,range);
% CDOT plots colored dots over the range of data
%
% Usage:  cdot(x,y,z,c,dot_size,label,range);
%
%   x=vector of x locations
%   y = vector of x locations
%   z = vector of z values
%   c = colormap
%   dot_size = Marker size for dots (e.g 8)
%   label = 1 to label with value, 0 otherwise
%   range = range of data for colormap
%
%  example
%    [x,y,z]=peaks(20);cdot(x(:),y(:),z(:),jet,20,1,[-8 8])
%
 [m,n]=size(c);
 ndots=length(z);
 zsave=z;
 if(exist('range')==1),
  zmin=range(1);
  zmax=range(length(range))+eps;
  iz=find(z>=zmax);
  if(~isempty(iz)),
    z(iz)=(zmax-eps)*ones(size(iz));
  end
  iz=find(z<zmin);
  if(~isempty(iz)),
    z(iz)=zmin*ones(size(iz));
  end
 else
  zmin=min(z(:));
  zmax=max(z(:))+eps;
 end
 zinc=(zmax-zmin)/m;
% set(gca,'xlim',[min(x) max(x)],'ylim',[min(y) max(y)]);
 for i=1:m;
   z1=zmin+zinc*(i-1);
   z2=zmin+zinc*i;
   ind=find(z>=z1& z<=z2);
   line(x(ind),y(ind),'linestyle','none','marker','.','markersize',...
      dot_size,'color',c(i,:));
 end
 if(label==1);
     for i=1:ndots
%      text(x(i),y(i),sprintf(' %d',zsave(i)),'fontsize',14,...
    text(double(x(i)),double(y(i)),sprintf('%5.2f',zsave(i)),...
     'HorizontalAlignment','left','VerticalAlignment','bottom')
%      bgtext(x(i),y(i),sprintf(' %d',zsave(i)),[1 1 .9],'fontsize',10,...
%       bgtext(x(i),y(i),sprintf(' %d',zsave(i)),[1 1 .9]);
    end
 end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
% NORMAL_LOGPDF
%
% Y = NORMAL_LOGPDF(X)
% Y = NORMAL_LOGPDF(X, MU, SIGMA)
%
% MU is the location parameter in (-INF, INF), default 0
% SIGMA is the scale parameter in (0, INF)

% Last modified 2010-11-12
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function y = lognormal_logpdf(x, varargin)

if nargin <= 3
  if nargin == 1
    mu = 0;
    sigma = 1;
  else
    mu = varargin{1};
    sigma = varargin{2};
  end
  y = -0.5*log(2*pi) ...
      -log(sigma) ...
      -0.5*((x-mu)./sigma).^2;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function b_im = Blur_IM( im, H )

[h w ch]  =  size(im);

if ch==3
    m1  = im(:,:,1);
    m2  = im(:,:,2);
    m3  = im(:,:,3);
    b_im(:,:,1)  =  reshape( H*m1(:), h, w );
    b_im(:,:,2)  =  reshape( H*m2(:), h, w );
    b_im(:,:,3)  =  reshape( H*m3(:), h, w );
else
    b_im  =  reshape( H*im(:), h, w );
end

author: lbasek

repo: image-denoising-benchmark


text:
% NMFLAB for Signal Processing written by A. Cichocki and R. Zdunek
% in cooperation with other members of Laboratory for Advanced Brain Signal
% Processing, BSI, RIKEN, Saitama, JAPAN

function [A,X,Distance_output]=nmf_smart_alg(Y,r,Index_norm,A_true,S, mc_on, NoAlts,restart_mc_on,AL,Y_true,type_alg_A,type_alg_X, max_restart,no_iter, alpha, alpha_BE, omegaX, omegaA, alphaS)
%
%
% Non-negative Matrix Factorization (NMF) with Extended SMART algorithms
%
% [A,X]=nmf_smart_alg(Y, r, Index_norm, A_true, S, mc_on, NoAlts, restart_mc_on, AL, Y_true, type_alg_A, type_alg_X, max_restart, no_iter,
%       alpha, alpha_BE, omegaX, omegaA, alphaS)
%       returns mixing matrix A of dimension [m by r],
%       and source matrix X of dimension [r by T], for the linear mixing model: AX = Y,
%       where Y is an observation matrix [m by T].
% Note: > m: number of sensors,
%       > r: number of sources,
%       > T: number of samples,
%
% INPUTS:
%       > Index_norm:    vector of 13 binary entries indicating which number
%                        divergence measures are turned on in the View Options,
%
%       > A_true:        true mixing matrix (only for synthetic data),
%       > S:             true source matrix (only for synthetic data),
%       > mc_on:         1 - Monte Carlo analysis enabled, 0 - otherwise,
%       > NoAlts:        number of alternating steps (only for Monte Carlo
%                        analysis and the option "Fixed Alternatings" is selected)
%       > restart_mc_on: 1 - restarts in Monte Carlo analysis are enabled,
%                        0 - otherwise,
%       > AL:            mixing matrix estimaed from the preceeding layer,
%       > Y_true:        the first layer mixed signals (mixtures),
%       > type_alg_A:    indicates the selected algorithm for computation
%                        of the mixing matrix,
%       > type_alg_X:    indicates the selected algorithm for computation of the sources,
%       > no_iter:       number of inner iterations,
%       > alpha:         parameter "alpha" in the Amari alpha-divergence
%       > alpha_BE:      parameter "alpha" in the Bose-Einstein divergence
%       > omegaX:        relaxation parameter (omega) in computation of the sources,
%       > omegaA:        relaxation parameter (omega) in computation of
%                        the mixing matrix,
%       > alphaS:        parameter of non-linear projection in computation
%                        of the sources,
%
% OUTPUTS:
%       > A:               estimated mixing matrix,
%       > X:               estimated source matrix,
%       > Distance_output: structures of different divergences measured between "Y" and estimated "AX" versus iterations,
%
% #########################################################################
A = [];
X = [];
if (nargin < 19) | isempty(alphaS) | max(size(alphaS) > 1)
   disp('Incorrect parameter alphaS');
   return
end
if (nargin < 18) | isempty(omegaA) | max(size(omegaA) > 1)
   disp('Incorrect parameter alpha');
   return
end
if (nargin < 17) | isempty(omegaX) | max(size(omegaX) > 1)
   disp('Incorrect parameter omegaX');
   return
end
if (nargin < 16) | isempty(alpha_BE) | max(size(alpha_BE) > 1)
   disp('Incorrect parameter alpha in the Bose_Einstein divergence');
   return
end
if (nargin < 15) | isempty(alpha) | max(size(alpha) > 1)
   disp('Incorrect parameter alpha');
   return
end
if (nargin < 14) | isempty(no_iter) | (no_iter < 1) | max(size(no_iter) > 1)
   disp('Incorrect number of inner iterations');
   return
end
if (nargin < 13) | isempty(max_restart)  | (max_restart < 0) | max(size(max_restart) > 1)
   disp('Number of restarts must be given correctly');
   return
end
if (nargin < 12) | isempty(type_alg_X) | (type_alg_X < 1) | max(size(type_alg_X) > 1)
   disp('Incorrect algorithm for X');
   return
end
if (nargin < 11) | isempty(type_alg_A) | (type_alg_A < 1) | max(size(type_alg_A) > 1)
   disp('Incorrect algorithm for A');
   return
end
if (nargin < 10) | isempty(Y_true)
   disp('The first layer mixed signals are unknown');
   Y_true = zeros(size(Y_true));
end
if (nargin < 9) | isempty(AL)
   disp('Mixing matrix from the preceeding layer unknown');
   AL = eye(size(Y,1));
end
if (nargin < 8) | isempty(restart_mc_on) | max(size(restart_mc_on) > 1)
   disp('Index od restarts in MC analysis unknown');
   restart_mc_on = 0;
end
if (nargin < 7) | isempty(NoAlts) | max(size(NoAlts) > 1)
   disp('Adjustable number of alternatings');
   NoAlts = [];
end
if (nargin < 6) | isempty(mc_on) | max(size(mc_on) > 1)
   disp('No Monte Carlo Analysis');
   mc_on = 0;
end
if (nargin < 5) | isempty(S)
   disp('X_true not given');
end
if (nargin < 4) | isempty(A_true)
   disp('A_true not given');
   index_fixed_A = 1;
else
   index_fixed_A = 0;
end
if (nargin < 3) | isempty(Index_norm)
   '"Index_norm" must be specified'
   return
end
if (nargin < 2) | isempty(r)
   'Rank of factorization must be given'
   return
end
if isempty(Y) | isnan(Y)
   error('No data');
   return
end
% test for negative values in Y
if min(min(Y)) < 0
    disp('Some matrix entries are changed from negative to small positive');
    Y(Y< 0) = eps;
end
if min(sum(Y,2)) == 0
    disp('Not all entries in a row can be zero');
    return
end

if (alpha == 0)
    if (type_alg_X == 1)
        type_alg_X = 2;
    end
    if (type_alg_A == 1)
        type_alg_A = 2;
    end
end

if (type_alg_A == 12) & (size(A_true,1) ~= size(Y,1))
   disp('Multilayer technique cannot be used with A fixed');
   Distance_output = [];
   return
end

Y = Y + eps;

[m,T]=size(Y);
niter_selected = 2000;     % maximum number of iterations for the selected sample (can be adjusted)
niter_sample = 30; % maximum number of iterations for each random sample
epsil_normA = 1E-4; % tolerance for alternating

% Monte Carlo and alternatings adjustment
if mc_on & ~restart_mc_on
    max_restart = 0;
end
if ~isempty(NoAlts)
    niter_selected = NoAlts;
end

% Declaration for A and X
A=zeros(m,r);
Ap = A;
X=zeros(r,T);
Ainit = A;
Xinit = X;
nr_best = -1;
Z = zeros(m,T);
KL_outer_temp = 0;
nr = 0; restart_on = 0; norm_A = 10;
m_sx = 1:m; r_sx = 1:r; T_sx = 1:T; s_dist = 0;

while (nr <= max_restart)

     % Initialize random A and X
       if ~nr & (~mc_on | restart_mc_on)
          Ainit(m_sx',r_sx) = abs(repmat(.1*sin(2*pi*.1*m_sx'),1,r) + repmat(.1*cos(2*pi*.1*r_sx),m,1) + repmat(cos(2*pi*.471*m_sx'),1,r) + repmat(sin(2*pi*.471*r_sx),m,1));
          Ainit = Ainit/max(max(Ainit));

          Xinit(r_sx',T_sx) = abs(repmat(.1*sin(2*pi*.1*r_sx'),1,T) + repmat(.1*cos(2*pi*.1*T_sx),r,1) + repmat(cos(2*pi*.471*r_sx'),1,T) + repmat(sin(2*pi*.471*T_sx),r,1));
          Xinit = Xinit/max(max(Xinit));
       else
          Ainit=rand(m,r);
          Xinit=rand(r,T);
       end

        % Normalization of initial guess
        Ainit = Ainit*diag(1./sum(Ainit,1));

        if (nr == max_restart)&(max_restart > 0)
           A = A_best;
           X = X_best;
        else
           A = Ainit;
           X = Xinit;
        end % initial guess assignment

    Yx = zeros(m,T);
    n = 0; k = 0;

while ((k <= niter_sample)&(nr < max_restart)) | ((k <= niter_selected)&(nr == max_restart)&(norm_A > epsil_normA)& isempty(NoAlts)) | ((k <= niter_selected)&(nr == max_restart)& (NoAlts > 0))

k = k + 1;

if no_iter == 1

    switch type_alg_X

        case 1 % Asymmetric alpha-divergence

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (1/alpha)*((Y./Z).^alpha - 1);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 2 % EG-DKL

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(Y./Z + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 3 % EG-RAG

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(2*Y./(Z + Y) + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 4 % EG-SAG

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(2*sqrt(Y.*Z)./(Z + Y) + eps) + (Y - Z)./Z ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 5 % EG-DJ

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = .5*log(Y./Z + eps) + .5*(Y - Z)./Z ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 6 % EG-RJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (Y - Z)./(Y + Z) ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 7 % EG-DRJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = 2*log(.5*(Y + Z)./Z) + (Z - Y)./(Y + Z);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 8 % EG-SJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(.5*(Y + Z)./Z + eps) ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 9 % EG-Tri

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (2*Y./(Y + Z)).^2 - 1;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 10 % EG-BE

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = alpha_BE*log((Y + alpha_BE*Z)./((1 + alpha_BE)*Z) + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1 + alphaS);

        case 11 % Pinv

            X = max(1E6*eps,pinv(A)*Y);

        case 12 % Fixed X

            X = S + eps;
            niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
            X = diag(sum(X,2))*X;

    end % switch for X

 switch type_alg_A

        case 1 % Asymmetric alpha-divergence

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (1/alpha)*((Y./Z).^alpha - 1);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

        case 2 % EG-DKL

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(Y./Z + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

        case 3 % EG-RAG

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(2*Y./(Z + Y) + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 4 % EG-SAG

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(2*sqrt(Y.*Z)./(Z + Y) + eps) + (Y - Z)./Z ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 5 % EG-DJ

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = .5*log(Y./Z + eps) + .5*(Y - Z)./Z ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 6 % EG-RJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (Y - Z)./(Y + Z) ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 7 % EG-DRJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = 2*log(.5*(Y + Z)./Z) + (Z - Y)./(Y + Z);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 8 % EG-SJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(.5*(Y + Z)./Z + eps) ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 9 % EG-Tri

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (2*Y./(Y + Z)).^2 - 1;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 10 % EG-BE

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = alpha_BE*log((Y + alpha_BE*Z)./((1 + alpha_BE)*Z) + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 11 % Pinv

            Ap = A;
            A = max(1E6*eps,Y*pinv(X')');
            A = A*diag(1./sum(A,1));

          case 12 % Fixed A

            niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
            A = A_true + eps;
            A = A*diag(1./(sum(A,1) + eps));

    end % switch for A

else % no_iter

   for t = 1:no_iter

        switch type_alg_X

        case 1 % Asymmetric alpha-divergence

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (1/alpha)*((Y./Z).^alpha - 1);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 2 % EG-DKL

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(Y./Z + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 3 % EG-RAG

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(2*Y./(Z + Y) + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 4 % EG-SAG

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(2*sqrt(Y.*Z)./(Z + Y) + eps) + (Y - Z)./Z ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 5 % EG-DJ

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = .5*log(Y./Z + eps) + .5*(Y - Z)./Z ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 6 % EG-RJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (Y - Z)./(Y + Z) ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 7 % EG-DRJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = 2*log(.5*(Y + Z)./Z) + (Z - Y)./(Y + Z);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 8 % EG-SJS

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = log(.5*(Y + Z)./Z + eps) ;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 9 % EG-Tri

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = (2*Y./(Y + Z)).^2 - 1;
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 10 % EG-BE

            gamma_m = sum(A,1);
            Z = A*X + eps;
            Q = alpha_BE*log((Y + alpha_BE*Z)./((1 + alpha_BE)*Z) + eps);
            X = (X.*exp(diag(omegaX./gamma_m)*(A'*Q))).^(1+alphaS);

        case 11 % Pinv

            X = max(1E6*eps,pinv(A)*Y);

        case 12 % Fixed X

            X = S + eps;
            niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
            X = diag(sum(X,2))*X;

    end % switch for X

   end % for t in X

   for t = 1:no_iter

       switch type_alg_A

        case 1 % Asymmetric alpha-divergence

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (1/alpha)*((Y./Z).^alpha - 1);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

        case 2 % EG-DKL

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(Y./Z + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

        case 3 % EG-RAG

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(2*Y./(Z + Y) + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 4 % EG-SAG

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(2*sqrt(Y.*Z)./(Z + Y) + eps) + (Y - Z)./Z ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 5 % EG-DJ

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = .5*log(Y./Z + eps) + .5*(Y - Z)./Z ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 6 % EG-RJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (Y - Z)./(Y + Z) ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 7 % EG-DRJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = 2*log(.5*(Y + Z)./Z) + (Z - Y)./(Y + Z);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

         case 8 % EG-SJS

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = log(.5*(Y + Z)./Z + eps) ;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 9 % EG-Tri

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = (2*Y./(Y + Z)).^2 - 1;
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 10 % EG-BE

            Ap = A;
            gamma_T = (sum(X,2))';
            Z = A*X + eps;
            Q = alpha_BE*log((Y + alpha_BE*Z)./((1 + alpha_BE)*Z) + eps);
            A = A.*exp(Q*X'*diag(omegaA./gamma_T));
            A = A*diag(1./(sum(A,1) + eps));

          case 11 % Pinv

            Ap = A;
            A = max(1E6*eps,Y*pinv(X')');
            A = A*diag(1./sum(A,1));

          case 12 % Fixed A

            niter_selected = 1000;     % maximum number of iterations for the selected sample (can be adjusted)
            A = A_true + eps;
            A = A*diag(1./(sum(A,1) + eps));

        end % switch for A

   end % for t in A

end % if no_iter

                if (nr == max_restart)&(mod(k,50)==0)& (~mc_on | restart_mc_on)
                    norm_A = norm(abs(A - Ap),'fro');
                    fprintf(1, 'Restart %d,  %d-th alternating step\n',nr_best+1,k);
                end

                if sum(Index_norm)
                   if (nr == max_restart) & (((k < 50) & (mod(k,5)==0)) | ((k>49) & ((mod(k,50)==0))))

                       s_dist = s_dist + 1;
                       k_select(s_dist) = k;
                       Z = A*X + eps;
                       Z = diag(1./sqrt(var(Z')))*Z;

                        dist_Fro(s_dist) = norm(Y - Z,'fro');
                        dist_KL(s_dist) = sum(sum(Y.*log(Y./Z + eps) - Y + Z));
                        dist_KL2(s_dist) = sum(sum(Z.*log(Z./Y + eps) + Y - Z));
                        dist_Pearson(s_dist) = sum(sum( ((Y - Z).^2)./Z ));
                        dist_Hellinger(s_dist) = sum(sum( (sqrt(Z) - sqrt(Y)).^2 ));
                        dist_JS_rel(s_dist) = sum(sum(2*Y.*log(2*Y./(Y + Z) + eps) + Z - Y));
                        dist_JS_rel2(s_dist) = sum(sum(2*Z.*log(2*Z./(Y + Z) + eps) - Z + Y));
                        Zy = Y + Z;
                        dist_JS(s_dist) = sum(sum(Y.*log(2*Y./Zy + eps) + Z.*log(2*Z./Zy + eps) ));
                        dist_AG_rel(s_dist) = sum(sum(Zy.*log(.5*Zy./Y + eps) + Y - Z));
                        dist_AG(s_dist) = sum(sum(.5*Zy.*log(.5*Zy./sqrt(Y.*Z) + eps)));
                        dist_J(s_dist) = sum(sum( .5*(Y - Z).*log(Y./Z + eps) ));
                        dist_Chi(s_dist) = sum(sum( ((Y + Z).*(Y - Z).^2)./(Y.*Z) ));
                        dist_Tria(s_dist) = sum(sum( ((Y - Z).^2)./(Y + Z) ));
                    end % if multiple
                end % if sum

end % while (k)

% Outer KL divergence
Z = AL*A*X;
Z_outer = norm(Z,'fro') + eps;
KL_outer = sum(sum(Y_true.*log((Y_true + eps)./(Z + eps)) - Y_true + Z))/Z_outer;

           if (nr == 0) | (KL_outer < KL_outer_temp)
              A_best = A; X_best = X; KL_outer_temp = KL_outer; nr_best = nr;
           end % multi-conditions

   nr = nr + 1;

   if nr <=max_restart
      fprintf(1, 'Restart %d, Kullback-Leibler divergence = %e\n',	nr, KL_outer);
   end

end % while (restarts)

X(X <= 0) = eps;

Distance_output = cell(length(s_dist),1);
Distance_output(1) = {[]};
Distance_output(2) = {[]};
Distance_output(3) = {[]};
Distance_output(4) = {[]};
Distance_output(5) = {[]};
Distance_output(6) = {[]};
Distance_output(7) = {[]};
Distance_output(8) = {[]};
Distance_output(9) = {[]};
Distance_output(10) = {[]};
Distance_output(11) = {[]};
Distance_output(12) = {[]};
Distance_output(13) = {[]};
Distance_output(14) = {[]};

if sum(Index_norm)
   Distance_output(1) = {k_select};
   Distance_output(2) = {dist_Fro};
   Distance_output(3) = {dist_KL};
   Distance_output(4) = {dist_KL2};
   Distance_output(5) = {dist_Pearson};
   Distance_output(6) = {dist_Hellinger};
   Distance_output(7) = {dist_JS_rel};
   Distance_output(8) = {dist_JS_rel2};
   Distance_output(9) = {dist_JS};
   Distance_output(10) = {dist_AG_rel};
   Distance_output(11) = {dist_AG};
   Distance_output(12) = {dist_J};
   Distance_output(13) = {dist_Chi};
   Distance_output(14) = {dist_Tria};
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [mimgR,mimgG] = red_channel_mean3(ops)

% numPlanes = length(ops.planesToProcess);

% build file list with red channel

if (isfield(ops, 'SubDirsRed') && ~isempty(ops.SubDirsRed))
    subDirsRed = ops.SubDirsRed;
else
    if (isfield(ops, 'expred') && ~isempty(ops.expred))
        for i = 1:length(ops.expred)
            subDirsRed{i} = sprintf('%d', ops.expred(i));
        end
    else
        warning('could not find red channel info, returning...')
        return;
    end
end
%
% build file list
fsRED = dir(fullfile(ops.RootDir, subDirsRed{1}, '*.tif'));
for k = 1:length(fsRED)
    fsRED(k).name = fullfile(ops.RootDir, subDirsRed{1}, fsRED(k).name);
end

root = ops.ResultsSavePath;
fregops =  sprintf('regops_%s_%s.mat', ops.mouse_name, ops.date);
if exist(fullfile(root, fregops), 'file')
    load(fullfile(root, fregops))
else
    ops1 = cell(ops.nplanes, 1);
    for j = 1:ops.nplanes

        fname = sprintf('regops_%s_%s_plane%d.mat', ops.mouse_name, ops.date, j);
        dat = load(fullfile(root, fname));
        ops1{j} = dat.ops;
        ops1{j}.useGPU = ops.useGPU;
    end
end

%
ntf0 = 0;
numPlanes = ops.nplanes;
%iplane0 = 1:1:ops.nplanes;

totFrames=0;
for k = 1:length(fsRED)
    %iplane0 = mod(iplane0-1, numPlanes) + 1;
    startPlane=((mod(totFrames+1,ops.nplanes*2)-1)/2)+1;

    nFr = nFramesTiff(fsRED(k).name);
    totFrames=totFrames+nFr;
    data = loadFramesBuff(fsRED(k).name, 1, nFr, 1, ops.temp_tiff);

    if ~exist('mimgR', 'var')
        [Ly, Lx, ~] = size(data);
        mimgR = zeros(Ly, Lx, ops.nplanes);
    end
    if ~exist('mimgG', 'var')
        [Ly, Lx, ~] = size(data);
        mimgG = zeros(Ly, Lx, ops.nplanes);
    end
    %

    for iPlane=1:ops.nplanes

        idx0=mod((ops.nplanes-startPlane+iPlane)*2,ops.nplanes*2);
        planesG=(idx0+1):(2*ops.nplanes):nFr;
        planesR=(idx0+2):(2*ops.nplanes):nFr;
        dataG0=data(:,:,planesG);
        dataR0=data(:,:,planesR);

        BiDiPhase = ops1{iPlane}.BiDiPhase;
        if abs(BiDiPhase) > 0
            yrange = 2:2:Ly;
            if BiDiPhase>0
                dataG0(yrange,(1+BiDiPhase):Lx,:,:) = dataG0(yrange, 1:(Lx-BiDiPhase),:,:);
                dataR0(yrange,(1+BiDiPhase):Lx,:,:) = dataR0(yrange, 1:(Lx-BiDiPhase),:,:);
            else
                dataG0(yrange,1:Lx+BiDiPhase,:,:)   = dataG0(yrange, 1-BiDiPhase:Lx,:,:);
                dataR0(yrange,1:Lx+BiDiPhase,:,:)   = dataR0(yrange, 1-BiDiPhase:Lx,:,:);
            end
        end

        [ds, ~]  = regoffKriging(dataG0, ops1{iPlane}, 0);
        %[ds, ~]  = registration_offsets(dataG0, ops1{iPlane}, 0);

        if k==1
            ds(1,:) = 0;
        end
        dataR       = ...
            register_movie(dataR0, ops1{iPlane}, ds);
        dataG     = ...
            register_movie(dataG0, ops1{iPlane}, ds);

        mimgR(:,:,iPlane) = mimgR(:,:,iPlane) + mean(dataR, 3);
        mimgG(:,:,iPlane) = mimgG(:,:,iPlane) + mean(dataG, 3);

    end

    ntf0 = ntf0 + 1;

    %iplane0 = iplane0 - nFr/ops.nchannels_red;
    fprintf('processing tiff %d/%d\n',k,length(fsRED))
end

mimgR = mimgR/ntf0;
mimgG = mimgG/ntf0;

author: cortex-lab


text:
function drawAcrobot(t,z,p)

clf; hold on;
length = p.l1+p.l2;
axis equal; axis(length*[-1,1,-1,1]);

[p1,p2] = acrobotKinematics(z,p);
pos = [[0;0],p1,p2];

plot(0,0,'ks','MarkerSize',25,'LineWidth',4)
plot(pos(1,:),pos(2,:),'Color',[0.1, 0.8, 0.1],'LineWidth',4)
plot(pos(1,:),pos(2,:),'k.','MarkerSize',50)

title(['Acrobot Animation,  t = ' num2str(t,2)])

drawnow; pause(0.01);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function sysdfod=dfod3(n,T,r)
%
% sysdfod=dfod3(n,T,r): digital fractional - order differentiator (r > 0)
%                       and integrator (r < 0) in form of the IIR filter;
%                       Recommended restriction for order r: (-1 < r < 1)
% Output: =>
% Discrete system in the form of the IIR filter of the given order 'n'
% obtained by power series expansion of the trapezoidal (Tustin) rule.
%
% Inputs: <=
% n: order of truncation (min. n = 20 is recommended) --> filter order
% T: sampling period in [sec]
% r: approximated fractional order (s^r), r is an arbitrary real number
%
% Copyright (c), 2011, Ivo Petras (ivo.petras@tuke.sk)
%
% Note: It requires a Matlab Control System Toolbox (->FILT function<-)
%
% Example: fractional half order integrator for T=0.1 sec and n = 20 :
% >> FHOI=dfod3(20, 0.1, -0.5); bode(FHOI); figure; step(FHOI);

bcN(1)=1.0; bcD(1)=1.0;
for i=1:n
  bcN(i+1)=((-1)^i)*(gamma(abs(r)+1)./(gamma(i+1).*gamma(abs(r)-i+1)));
  bcD(i+1)=gamma(abs(r)+1)./(gamma(i+1).*gamma(abs(r)-i+1));
end
if r>=0
  sysdfod=((2/T)^r)*(filt(bcN,bcD,T));
end
if r<0
  sysdfod=((2/T)^r)*(filt(bcD,bcN,T));
end
%


text:
function pf = noisepf(pf,fak,background,varargin)
% simulate diffraction counts
%
% noisepf simulates realistic diffraction counts by generating random
% samples of the Poisson distribution with mean m = alpha * pdf + bg
%
% Syntax
%   pdfn = noisepf(pdf,alpha,bg)
%
% Input
%  pf    - @PoleFigure
%  alpha - uniform radiation (double)
%  bg    - background radiation (double)
%
% Options
%  NONNEGATIV - force data to be non negative
%
% See also
% SO3Fun/calcPoleFigure

if nargin == 2, background = 0;end
if numel(fak) == 1, fak = repmat(fak,pf.numPF,1);end

for i = 1:pf.numPF
  pf.allI{i} = ...
  randp(fak(i)*pf.allI{i} + background) - background;
  if check_option(varargin,'NONNEGATIV')
    pf.allI{i}(pf.allI{i} < 0) = 0;
  end
end

author: mtex-toolbox

repo: mtex

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


text:
function piece_processing
global ai ps st hp pc paus frs datas
global highpass lowpass numerator_hp numerator_lp
if st
    stop(ai);
    delete(ai);
else
    data = getdata(ai,ps);
    %datas(1:ps)=datas(ps+1:2*ps);
    datas(1:2*ps)=datas(1*ps+1:3*ps);
    %datas(ps+1:2*ps)=data;
    datas(2*ps+1:3*ps)=data;
    pc=pc+1;
    if (mod(pc,frs)==0)&&(~paus)
        [tmp ii20]=max(data); % syncrhinize to maximum
        %datas1=datas(ii20:ps+ii20-1);
        if (~highpass)&&(~lowpass)
            datas1=datas(ps+ii20:2*ps+ii20-1);
            %datas1=datas(ps+ii20-0:2*ps+ii20-1);
            %datas1(0+1:end)
        end
        if highpass
            datas10=datas(ps+ii20-300:2*ps+ii20-1);
            datas10=filter(numerator_hp,1,datas10);
            datas1=datas10(301:end);
        end
        if lowpass
            datas10=datas(ps+ii20-300:2*ps+ii20-1);
            datas10=filter(numerator_lp,1,datas10);
            datas1=datas10(301:end);
        end
        set(hp,'YData',datas1);
        drawnow;
    end
end


text:
function value = r8_chu ( a, b, x )

%*****************************************************************************80
%
%% R8_CHU evaluates the confluent hypergeometric function of R8 arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, B, the parameters.
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the function value.
%
  persistent eps

  if ( isempty ( eps ) )
    eps = r8_mach ( 3 );
  end

  if ( x < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_CHU - Fatal error!\n' );
    fprintf ( 1, '  X < 0.\n' );
    error ( 'R8_CHU - Fatal error!' )
  end

  if ( x == 0.0 )
    if ( 1.0 <= b )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8_CHU - Fatal error!\n' );
      fprintf ( 1, '  X = 0 and 1 <= B.\n' );
      error ( 'R8_CHU - Fatal error!' )
    end
    value = r8_gamma ( 1.0 - b ) / r8_gamma ( 1.0 + a - b );
    return
  end

  if ( max ( abs ( a ), 1.0 ) * max ( abs ( 1.0 + a - b ), 1.0 ) ...
    < 0.99 * abs ( x ) )
    value = x^( - a ) * r8_chu_scaled ( a, b, x );
    return
  end
%
%  The ascending series will be used, because the descending rational
%  approximation (which is based on the asymptotic series) is unstable.
%
  if ( 0.0 <= b )
    aintb = r8_aint ( b + 0.5 );
  else
    aintb = r8_aint ( b - 0.5 );
  end
  beps = b - aintb;
  n = aintb;
  alnx = log ( x );
  xtoeps = exp ( - beps * alnx );
%
%  Evaluate the finite sum.
%
%  Consider the case b < 1.0 first.
%
  if ( n < 1 )

    sum = 1.0;

    t = 1.0;
    m = - n;
    for i = 1 : m
      xi1 = i - 1;
      t = t * ( a + xi1 ) * x / ( ( b + xi1 ) * ( xi1 + 1.0 ) );
      sum = sum + t;
    end

    sum = r8_poch ( 1.0 + a - b, - a ) * sum
%
%  Now consider the case 1 <= b.
%
  else

    sum = 0.0;
    m = n - 2;

    if ( 0 <= m )

      t = 1.0;
      sum = 1.0;

      for i = 1 : m
       xi = i;
        t = t * ( a - b + xi ) * x / ( ( 1.0 - b + xi ) * xi );
        sum = sum + t;
      end

      sum = r8_gamma ( b - 1.0 ) * r8_gamr ( a ) * x^( 1 - n ) * xtoeps * sum;

    end

  end
%
%  Next evaluate the infinite sum.
%
  if ( n < 1 )
    istrt = 1 - n;
  else
    istrt = 0;
  end

  xi = istrt;

  factor = r8_mop ( n ) * r8_gamr ( 1.0 + a - b ) * x^istrt;

  if ( beps ~= 0.0 )
    factor = factor * beps * pi / sin ( beps * pi );
  end

  pochai = r8_poch ( a, xi );
  gamri1 = r8_gamr ( xi + 1.0 );
  gamrni = r8_gamr ( aintb + xi );
  b0 = factor * r8_poch ( a, xi - beps ) * gamrni ...
    * r8_gamr ( xi + 1.0 - beps );
%
%  x^(-beps) is close to 1.0, so we must be careful in evaluating
%  the differences.
%
  if ( abs ( xtoeps - 1.0 ) <= 0.5 )

    pch1ai = r8_poch1 ( a + xi, - beps );
    pch1i = r8_poch1 ( xi + 1.0 - beps, beps );
    c0 = factor * pochai * gamrni * gamri1 * ( ...
      - r8_poch1 ( b + xi,- beps ) + pch1ai ...
      - pch1i + beps * pch1ai * pch1i );
%
%  xeps1 = (1.0 - x^(-beps))/beps = (x^(-beps) - 1.0)/(-beps)
%
    xeps1 = alnx * r8_exprel ( - beps * alnx );

    value = sum + c0 + xeps1 * b0;
    xn = n;

    for i = 1 : 1000
      xi = istrt + i;
      xi1 = istrt + i - 1;
      b0 = ( a + xi1 - beps ) * b0 * x ...
        / ( ( xn + xi1 ) * ( xi - beps ) );
      c0 = ( a + xi1 ) * c0 * x / ( ( b + xi1) * xi ) ...
        - ( ( a - 1.0 ) * ( xn + 2.0 * xi - 1.0 ) ...
        + xi * ( xi - beps ) ) * b0 ...
        / ( xi * ( b + xi1 ) * ( a + xi1 - beps ) );
      t = c0 + xeps1 * b0;
      value = value + t;
      if ( abs ( t ) < eps * abs ( value ) )
        return
      end
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_CHU - Fatal error!\n' );
    fprintf ( 1, '  No convergence in 1000 terms.\n' );
    error ( 'R8_CHU - Fatal error!' )

  end
%
%  x^(-beps) is very different from 1.0, so the straightforward
%  formulation is stable.
%
  a0 = factor * pochai * r8_gamr ( b + xi ) * gamri1 / beps;
  b0 = xtoeps * b0 / beps;

  value = sum + a0 - b0;

  for i = 1 : 1000
    xi = istrt + i;
    xi1 = istrt + i - 1;
    a0 = ( a + xi1 ) * a0 * x / ( ( b + xi1 ) * xi );
    b0 = ( a + xi1 - beps ) * b0 * x ...
      / ( ( aintb + xi1 ) * ( xi - beps ) );
    t = a0 - b0;
    value = value + t;
    if ( abs ( t ) < eps * abs ( value ) )
      return
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_CHU - Fatal error!\n' );
  fprintf ( 1, '  No convergence in 1000 terms.\n' );
  error ( 'R8_CHU - Fatal error!' )

end

author: johannesgerer

repo: jburkardt-m


text:
function gibbsKernDisplay(kern, spaceNum)

% GIBBSKERNDISPLAY Display parameters of the GIBBS kernel.
% FORMAT
% DESC displays the parameters of the Mark Gibbs's non-stationary
% kernel and the kernel type to the console.
% ARG kern : the kernel to display.
%
% FORMAT does the same as above, but indents the display according
% to the amount specified.
% ARG kern : the kernel to display.
% ARG spacing : how many spaces to indent the display of the kernel by.
%
% SEEALSO : gibbsKernParamInit, modelDisplay, kernDisplay
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

if nargin > 1
else
  spaceNum = 0;
end
spacing = repmat(32, 1, spaceNum);
spacing = char(spacing);
fprintf(spacing);
fprintf('GIBBS variance: %2.4f\n', kern.variance)
fprintf(spacing);
fprintf('GIBBS length scale function: \n')
modelDisplay(kern.lengthScaleFunc, length(spacing) + 2);


text:
function [par, parb, Pf, Pb] = rc2parv(rc,rcb)

%function [par, parb, Pf, Pb] = rc2parv(rc,rcb)
%  Transforms forward and backward reflection matrices rc and rcb
%  into parameters.

%S. de Waele, March 2003.

s = kingsize(rc);
order = s(3)-1;
dim = s(1); I = eye(dim);

par = zeros(dim,dim,order+1);
parb = zeros(dim,dim,order+1);

par(:,:,1) = I;
parb(:,:,1)= I;
if order,
	par(:,:,2) = rc(:,:,2);
	parb(:,:,2)= rcb(:,:,2);
	par_o  = par;
	parb_o = parb;
end
for p = 2:order,
   par(:,:,2:p) =  par_o(:,:,2:p) +fliptime(filterv(rc(:,:,p+1),1,parb_o(:,:,2:p)));
   par(:,:,p+1)= rc(:,:,p+1);
   parb(:,:,2:p) =  parb_o(:,:,2:p) +fliptime(filterv(rcb(:,:,p+1) ,1,par_o(:,:,2:p)));
   parb(:,:,p+1)= rcb(:,:,p+1);

   par_o  = par;
   parb_o = parb;
end %for p = 2:order,


text:
function partitions = cosmo_nchoosek_partitioner(chunks_or_ds, k, varargin)
% partitions for into nchoosek(n,k) parititions with optional grouping schemas.
%
% partitions=cosmo_nchoosek_partitioner(chunks, k, group_values1, test_group_by1,...)
%
% Input
%  - chunks          Px1 chunk indices for P samples. It can also be a
%                    dataset struct with field .sa.chunks
%  - k               When an integer, k chunks are in each test partition.
%                    When between 0 and 1, this is interpreted as
%                    round(k*nchunks) where nchunks is the number of unique
%                    chunks in chunks.
%                    A special case, mostly aimed at split-half
%                    correlations, is when k=='half'; this sets k to .5,
%                    and if k is even, it treats train_indices and
%                    test_indices as symmetrical, meaning it returns only
%                    half the number of partitions (avoiding duplicates).
%                    If k is odd then train and test indices have
%                    (k+1)/nchunks and (k-1)/nchunks elements,
%                    respectively, or vice versa.
%  - group_values*   } Intended for cross-participant or cross-condition
%  - test_group_by*  } generalizability analyses.
%                      Pairs of these determine a subsequent level group
%                      partition scheme. Each group_values can be
%                      a cell with the labels for a test group, or,
%                      if chunks is a string, the name of a sample
%                      attribute whose values are used as labels.
%                      Each test_group_by indicates which values in
%                      group_values are used as a test value in the
%                      cross-validation scheme (and all other ones are used
%                      as training value). If empty it is set to the unique
%                      values in group_values.
%
% Output:
%  - partitions      A struct with fields:
%     .train_indices } Each of these is an 1xQ cell for Q partitions, where
%     .test_indices  } Q=nchoosek(nchunks,k)). It contains all possible
%                      combinations with k test indices and (nchunks-k)
%                      training indices (except when k=='half', see above).
%                    If group_values* and test_group_by* are specified,
%                    then the number of output partitions is multiplied
%                    by the product of the number of values in
%                    test_group_by*.
%
%
% Examples:
%     % make a simple dataset with 4 chunks, 2 samples each
%     % assume two targets (i.e. conditions, say piano versus guitar)
%     ds=struct();
%     ds.samples=randn(8,99); % 8 samples, 99
%     ds.sa.targets=[1 1 1 1 2 2 2 2]';
%     ds.sa.chunks=2+[1 2 3 4 4 3 2 1]';
%     cosmo_check_dataset(ds); % sanity check
%     %
%     % take-one-chunk out partitioning
%     p=cosmo_nchoosek_partitioner(ds,1);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 2    [ 1    [ 1    [ 1
%     %||       3      3      2      2
%     %||       4      4      4      3
%     %||       5      5      5      6
%     %||       6      6      7      7
%     %||       7 ]    8 ]    8 ]    8 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 2    [ 3    [ 4
%     %||       8 ]    7 ]    6 ]    5 ] }
%     %
%     % take-two chunks out partitioning
%     p=cosmo_nchoosek_partitioner(ds,2);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 3    [ 2    [ 2    [ 1    [ 1    [ 1
%     %||       4      4      3      4      3      2
%     %||       5      5      6      5      6      7
%     %||       6 ]    7 ]    7 ]    8 ]    8 ]    8 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 1    [ 1    [ 2    [ 2    [ 3
%     %||       2      3      4      3      4      4
%     %||       7      6      5      6      5      5
%     %||       8 ]    8 ]    8 ]    7 ]    7 ]    6 ] }
%     %
%     % take-half-of-the-chunks out partitioning
%     % (this effectively gives same chunks as above)
%     p_alt=cosmo_nchoosek_partitioner(ds,.5);
%     isequal(p, p_alt)
%     %|| true
%     %
%     % do half split (for correlation measure); this leaves out
%     % mirror partitions of train and test indices
%     p=cosmo_nchoosek_partitioner(ds,'half');
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 1    [ 1    [ 1
%     %||       2      3      4
%     %||       7      6      5
%     %||       8 ]    8 ]    8 ] }
%     %|| .test_indices
%     %||   { [ 3    [ 2    [ 2
%     %||       4      4      3
%     %||       5      5      6
%     %||       6 ]    7 ]    7 ] }
%     %
%     % test on samples with chunk=3 only using take-one-fold out
%     p=cosmo_nchoosek_partitioner(ds,1,'chunks',3);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 2
%     %||       3
%     %||       4
%     %||       5
%     %||       6
%     %||       7 ] }
%     %|| .test_indices
%     %||   { [ 1
%     %||       8 ] }
%     % test on samples with chunk=[3 4] only using take-one-fold out;
%     % only samples with chunks=1 or 2 are used for training
%     p=cosmo_nchoosek_partitioner(ds,1,'chunks',{[3 4]});
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 3    [ 3
%     %||       4      4
%     %||       5      5
%     %||       6 ]    6 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 2
%     %||       8 ]    7 ] }
%     % test separately on samples with chunk=3 and samples with chunk=4;
%     % in some folds, samples with chunks=1,2,4 are used for training, in
%     % other folds samples with chunks=1,2,3 are used for training
%     p=cosmo_nchoosek_partitioner(ds,1,'chunks',{[3],[4]});
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 2    [ 1
%     %||       3      3
%     %||       4      4
%     %||       5      5
%     %||       6      6
%     %||       7 ]    8 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 2
%     %||       8 ]    7 ] }
%     %||
%
%
%     % make a slightly more complicated dataset: with three chunks,
%     % suppose there are two modalities (e.g. (1) visual and (2)
%     % auditory stimulation) which are stored in an
%     % additional field 'modality'
%     ds=struct();
%     ds.samples=randn(12,99);
%     ds.sa.chunks  =[1 1 1 1 2 2 2 2 3 3 3 3]';
%     ds.sa.targets =[1 2 1 2 1 2 1 2 1 2 1 2]';
%     ds.sa.modality=[1 1 2 2 1 1 2 2 1 1 2 2]';
%     cosmo_check_dataset(ds);
%     %
%     % take-one-chunk out, test on samples with modality=1 (and train
%     % on samples with other modalities, i.e. modality=2)
%     p=cosmo_nchoosek_partitioner(ds,1,'modality',1);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [  7    [  3    [ 3
%     %||        8       4      4
%     %||       11      11      7
%     %||       12 ]    12 ]    8 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 5    [  9
%     %||       2 ]    6 ]    10 ] }
%     %
%     % take-one-chunk out, test on samples with modality=2
%     p=cosmo_nchoosek_partitioner(ds,1,'modality',2);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [  5    [  1    [ 1
%     %||        6       2      2
%     %||        9       9      5
%     %||       10 ]    10 ]    6 ] }
%     %|| .test_indices
%     %||   { [ 3    [ 7    [ 11
%     %||       4 ]    8 ]    12 ] }
%     % take-one-chunk out, test on samples with modality=1 (and train on
%     % modality=2) and vice verse
%     p=cosmo_nchoosek_partitioner(ds,1,'modality',[]);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [  7    [  3    [ 3    [  5    [  1    [ 1
%     %||        8       4      4       6       2      2
%     %||       11      11      7       9       9      5
%     %||       12 ]    12 ]    8 ]    10 ]    10 ]    6 ] }
%     %|| .test_indices
%     %||   { [ 1    [ 5    [  9    [ 3    [ 7    [ 11
%     %||       2 ]    6 ]    10 ]    4 ]    8 ]    12 ] }
%
%     % between-subject classification: 3 chunks, 2 modalities, 5 subjects
%     ds=struct();
%     ds.samples=randn(60,99);
%     ds.sa.targets=repmat([1 2],1,30)';
%     ds.sa.chunks=repmat([1 1 1 1 2 2 2 2 3 3 3 3],1,5)';
%     ds.sa.modality=repmat([1 1 2 2],1,15)';
%     ds.sa.subject=kron(1:5,ones(1,12))';
%     cosmo_check_dataset(ds);
%     %
%     % test on subject 1, train on other subjects using take-one-chunk out
%     p=cosmo_nchoosek_partitioner(ds,1,'subject',1);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 17         [ 13         [ 13
%     %||       18           14           14
%     %||       19           15           15
%     %||        :            :            :
%     %||       58           58           54
%     %||       59           59           55
%     %||       60 ]@32x1    60 ]@32x1    56 ]@32x1 }
%     %|| .test_indices
%     %||   { [ 1    [ 5    [  9
%     %||       2      6      10
%     %||       3      7      11
%     %||       4 ]    8 ]    12 ] }
%     %
%     % test on each subject after training on each other subject
%     % in each fold, the test data is from one subject and one chunk,
%     % and the train data from all other subjects and all other chunks.
%     % since there are 5 subjects and 3 chunks, there are 15 folds.
%     p=cosmo_nchoosek_partitioner(ds,1,'subject',[]);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 17         [ 13         [ 13        ... [  5         [  1         [  1
%     %||       18           14           14               6            2            2
%     %||       19           15           15               7            3            3
%     %||        :            :            :               :            :            :
%     %||       58           58           54              46           46           42
%     %||       59           59           55              47           47           43
%     %||       60 ]@32x1    60 ]@32x1    56 ]@32x1       48 ]@32x1    48 ]@32x1    44 ]@32x1   }@1x15
%     %|| .test_indices
%     %||   { [ 1    [ 5    [  9   ... [ 49    [ 53    [ 57
%     %||       2      6      10         50      54      58
%     %||       3      7      11         51      55      59
%     %||       4 ]    8 ]    12 ]       52 ]    56 ]    60 ]   }@1x15
%     %
%     % as above, but test on modality=2 (and train on other values for
%     % modality, i.e. modality=1)
%     p=cosmo_nchoosek_partitioner(ds,1,'subject',[],'modality',2);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 17         [ 13         [ 13        ... [  5         [  1         [  1
%     %||       18           14           14               6            2            2
%     %||       21           21           17               9            9            5
%     %||        :            :            :               :            :            :
%     %||       54           50           50              42           38           38
%     %||       57           57           53              45           45           41
%     %||       58 ]@16x1    58 ]@16x1    54 ]@16x1       46 ]@16x1    46 ]@16x1    42 ]@16x1   }@1x15
%     %|| .test_indices
%     %||   { [ 3    [ 7    [ 11   ... [ 51    [ 55    [ 59
%     %||       4 ]    8 ]    12 ]       52 ]    56 ]    60 ]   }@1x15
%     %
%     % as above, but test on each modality after training on the other
%     % modality. There are 30 folds (5 subjects, 3 chunks, 2 modalities).
%     p=cosmo_nchoosek_partitioner(ds,1,'subject',[],'modality',[]);
%     cosmo_disp(p);
%     %|| .train_indices
%     %||   { [ 19         [ 15         [ 15        ... [  5         [  1         [  1
%     %||       20           16           16               6            2            2
%     %||       23           23           19               9            9            5
%     %||        :            :            :               :            :            :
%     %||       56           52           52              42           38           38
%     %||       59           59           55              45           45           41
%     %||       60 ]@16x1    60 ]@16x1    56 ]@16x1       46 ]@16x1    46 ]@16x1    42 ]@16x1   }@1x30
%     %|| .test_indices
%     %||   { [ 1    [ 5    [  9   ... [ 51    [ 55    [ 59
%     %||       2 ]    6 ]    10 ]       52 ]    56 ]    60 ]   }@1x30
%
% Notes:
%    - When k=1 and two input arguments, this function behaves equivalently
%      to cosmo_nfold_partitioner. Thus, in the most simple case
%      (nfold-partitioning), cosmo_nfold_partitioner with k=1 can be used
%      as well as this function.
%    - This function does not consider any .sa.targets (trial condition)
%      or .samples information.
%    - As shown in the examples below, this function can be used for
%      cross-modal and/or cross-participant cross-validation.
%    - For cross-validation it is recommended to balance partitions using
%      cosmo_balance_partitions.
%   - this function can be used for cross-decoding analyses. Doing so may
%     require a re-assignment of .sa.targets, and adding another sample
%     attribute to specify which samples are used for training and testing.
%     For example, consider the following dataset with six unique
%     conditions as specified in the sample attribute field .sa:
%
%       .targets    .chunks     .labels
%       1           1           'vis_dog'
%       2           1           'vis_cat'
%       3           1           'vis_frog'
%       4           1           'aud_dog'
%       5           1           'aud_cat'
%       6           1           'aud_frog'
%       1           2           'vis_dog'
%       2           2           'vis_cat'
%       :           :               :
%       6           8           'aud_frog'
%
%    This dataset has 8 chunks, each with 6 conditions: three for visual
%    stimuli of dogs, cats, and frogs, and three for auditory stimuli for
%    the same animals. The field .labels is not required, but used for a
%    human-readable description of the condition of each sample.
%
%    Suppose that one wants to do cross-decoding to see
%    if discrimination of animals generalizes between the visual and
%    auditory modalities.
%    To do so, the user has to:
%       * change the .targets field, to indicate the anima species,
%       * add another field (here 'modality') indicating which samples are
%         used for the cross-decoding
%
%    In this example, the sample attribute field .sa can be set as follows:
%
%       .targets    .chunks     .labels     .modality
%       1           1           'vis_dog'   1
%       2           1           'vis_cat'   1
%       3           1           'vis_frog'  1
%       1           1           'aud_dog'   2
%       2           1           'aud_cat'   2
%       3           1           'aud_frog'  2
%       1           2           'vis_dog'   1
%       2           2           'vis_cat'   1
%       :           :               :
%       6           8           'aud_frog'  2
%
%    so that:
%     * .modality corresponds to the visual (=1) and auditory (=2)
%       modality.
%     * .targets corresponds to the dog (=1), cat (=2), and frog (=3)
%       species.
%
%    With this re-assignment of .targets, testing on the auditory modality
%    and training on the visual modality with take-one-chunk out
%    cross-validation can be done using:
%
%       test_modality=2; % train on all other modalities (here: only 1)
%       partitions=cosmo_nchoosek_partitioner(ds,1,...
%                                   'modality',test_modality);
%
%    If test_modality is set to empty ([]), then both modalities are used
%    for training and for testing (in separate folds).
%
% See also: cosmo_nfold_partitioner, cosmo_balance_partitions
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    % if the first input is a dataset, get the chunks from it
    is_ds=isstruct(chunks_or_ds);
    if is_ds
        if isfield(chunks_or_ds,'sa') && isfield(chunks_or_ds.sa,'chunks')
            % it is a dataset
            ds=chunks_or_ds;
            chunks=ds.sa.chunks;
        else
            error('illegal struct input');
        end
    else
        chunks=chunks_or_ds;
    end

    % use helper function defined below
    partitions=nchoosek_partitioner(chunks,k);

    if nargin<3
        return;
    elseif mod(nargin,2)~=0
        error('Need even number of arguments');
    else
        nsamples=numel(chunks);
        for k=1:2:(nargin-2);
            % loop over grouping values
            group_values=varargin{k};
            test_group_by=varargin{k+1};

            % get sample attribute, if input was dataset
            if ischar(group_values)
                if is_ds
                    group_values=ds.sa.(group_values);
                else
                    error('String for group value requires dataset input');
                end
            end

            % check number of values matches
            ngroup_by=numel(group_values);
            if ngroup_by ~= nsamples
                error('group_by has %d values, expected %d', ...
                                ngroup_by, nsamples);
            end

            if isempty(test_group_by)
                test_group_by=num2cell(unique(group_values));
            elseif isnumeric(test_group_by)
                test_group_by={test_group_by};
            end

            % update partitions using helper function below
            partitions=group_by(partitions, group_values, test_group_by);
        end
    end

function partitions=group_by(partitions, group_values, test_group_by)
    % helper function to group partitions
    % the output has N times as many partitions as the input,
    % where N=numel(test_by).

    npartitions=numel(partitions.test_indices);

    assert(iscell(test_group_by));
    % see which values to split on
    ntest_group_by=numel(test_group_by);

    % allocate space for output
    train_indices=cell(1,ntest_group_by);
    test_indices=cell(1,ntest_group_by);

    % run for each unique value in group_by_values
    for m=1:ntest_group_by
        test_by=test_group_by{m};

        % allocate space for this iteration
        train_indices_cell=cell(1,npartitions);
        test_indices_cell=cell(1,npartitions);

        % some filtered partitions may be empty, so keep track of the
        % last position where a value was stored
        pos=0;
        for j=1:npartitions
            % get testing chunk indices
            p_test=partitions.test_indices{j};

            % see which ones match the group_by_value
            msk_test=cosmo_match(group_values(p_test),test_by);

            % keep just those indices
            p_test_masked=p_test(msk_test);

            % the same for training, but different from test_by
            p_train=partitions.train_indices{j};
            msk_train=~cosmo_match(group_values(p_train),test_by);
            p_train_masked=p_train(msk_train);

            if ~isempty(p_test_masked) && ~isempty(p_train_masked)
                % both training and test set are non-empty, so keep result
                pos=pos+1;
                test_indices_cell{pos}=p_test_masked;
                train_indices_cell{pos}=p_train_masked;
            end
        end

        % store the test-indices for m-th group_by_value
        test_indices{m}=test_indices_cell(1:pos);
        train_indices{m}=train_indices_cell(1:pos);
    end

    partitions=struct();

    % concatenate results for training and test indices
    partitions.train_indices=cat(2,train_indices{:});
    partitions.test_indices=cat(2,test_indices{:});

function partitions=nchoosek_partitioner(chunks,k)
% straightfoward partitioner

[unq,unused,chunk_indices]=unique(chunks);
nchunks=numel(unq);

if nchunks<2
    error(['at least two unique values for .sa.chunks are required, '...
                    'found %d.'], nchunks);
end

% deal with special 'half' case
is_symmetric=false;
if ischar(k)
    if strcmp(k,'half')
        k=.5;
        is_symmetric=true;
    else
        error('illegal string k');
    end
end
if isnumeric(k)
    if ~isscalar(k)
        error('k must be a scalar');
    end

    if k ~= round(k)
        k=round(nchunks*k);
    end

    if ~any(k==1:(nchunks-1))
        error('illegal k=%d: test class count should be between 1 and %d', ...
                    k, nchunks-1);
    end
else
    error('illegal parameter for k');
end

% little optimization: if just two chunks, the split is easy
if all(cosmo_match(chunks,[1 2]))
    chunk_msk1=chunks==1;
    chunk1_idxs=find(chunk_msk1);
    chunk2_idxs=find(~chunk_msk1);

    partitions=struct();
    if is_symmetric
        partitions.train_indices={chunk1_idxs};
        partitions.test_indices={chunk2_idxs};
    else
        partitions.train_indices={chunk1_idxs,chunk2_idxs};
        partitions.test_indices={chunk2_idxs,chunk1_idxs};
    end
    return
end

npartitions=nchoosek(nchunks,k);
combis=nchoosek(1:nchunks,k); % make all combinations

if is_symmetric && mod(nchunks,2)==0
    % when nclasses is even, return the first half of the permutations:
    % the current implementation of nchoosek results is that
    % combis(k,:) and combis(npartitions+1-k) are complementary
    % (i.e. together they make up 1:nchunks). In principle this could
    % change in the future leading to wrong results (if Mathworks, in its
    % infinite wisdom, decides to change the implementation of nchoosek),
    % so to be sure we check that the output matches what is expected.
    % The rationale is that this reduces computation time of subsequent
    % analyses by half, in particular for commonly used split half
    % correlations.
    nhalf=npartitions/2;

    check_combis=[combis(1:nhalf,:) combis(npartitions:-1:(nhalf+1),:)];

    % each row, when sorted, should be 1:nchunks
    matches=bsxfun(@eq,sort(check_combis,2),1:nchunks);
    assert(all(matches(:)),'nchoosek behaves weirdly');

    % we're good - just take the second half and update npartitions
    combis=combis(npartitions:-1:(nhalf+1),:);
    npartitions=nhalf;
end

% allocate space for output
train_indices=cell(1,npartitions);
test_indices=cell(1,npartitions);

% make all partitions
for j=1:npartitions
    combi=combis(j,:);
    sample_count=cosmo_match(chunk_indices,combi);
    test_indices{j}=find(sample_count);
    train_indices{j}=find(~sample_count);
end

partitions=struct();
partitions.train_indices=train_indices;
partitions.test_indices=test_indices;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function c=countif(vector,x)
%COUNTIF adalah fungsi untuk menghitung jumlah kehadiran bilangan dalam
%suatu array
[m,n]=size(vector);
d=0;
for a=1:m
    if vector(a)==x
    d=d+1;
    c=d;
    else
    c=d;
    end
end
c;


text:
function slSpawn(searchlight, process, processes, partitions)
% slSpawn(searchlight, process, processes, partitions)
% SEARCHLIGHT - SPAWN
% ---------------------------------------------------------
% Used by SLRUN to complete bite-sized portions of a searchlight.
%
% INPUTS
% 	searchlight - Structure initialized with slInit.
% 	process - Which of N portions to complete.
% 	processes - Number of portions to be completed.
% 	partitions - How many smaller-sized portions to break up each portion
%       into.
%
% OUTPUTS
%   N/A
%
% USAGE
%   Calling this is either handled by the SGE or layed out explicitly by
%   provided commands.  Assuming you'd like to write up a command on your
%   own, let's suppose we want to run the first portion of a searchlight
%   broken into three pieces.  Each of those portions we'll divide into 30
%   bite-sized portions to iterate over.
%       searchlight = slInit(...);
%       slSpawn(searchlight, 1, 3, 30);
%
% See also SLINIT, SLRUN, SLCOMPLETE, SVMINIT.
%
% renobowen@gmail.com [2010]
%
    isGray = strcmpi(searchlight.searchView, 'gray');
    growBy = searchlight.growBy;

    processesToRun = (((process - 1) * partitions):(process * partitions - 1)) + 1;
    for i = processesToRun
        filename = sprintf('%s/%s/searchlight%04d.mat', searchlight.path, searchlight.tmpDir, i);
        processString = sprintf('[%d/%d]', i, processes);
        if (exist(filename, 'file'))
            fprintf('%s Process already completed.  Proceeding...\n', processString);
            continue;
        end

        indices = searchlight.startInds(i):searchlight.endInds(i);

        if (isGray)
            centers = searchlight.ipInds(:, indices);
        else
            centers = searchlight.allCoords(:, indices);
        end
        nCenters = size(centers, 2);

        fprintf('%s Processing coordinate block...', processString);
        meanAccs = zeros(1, nCenters);
        for j = 1:nCenters
            [roiCoords centers(:, j)] = growROI(searchlight, centers(:, j), isGray);
            if (isempty(roiCoords))
                meanAccs(j) = 0;
            else
                meanAccs(j) = runSVM(searchlight, roiCoords);
            end
        end
        fprintf('done.\n');
        fprintf('%s Saving...', processString);
        save(filename, 'meanAccs', 'centers', 'growBy', '-v7.3');
        fprintf('done.\n');
    end
end

function [roiCoords center] = growROI(searchlight, center, isGray)
    growBy = searchlight.growBy;
    if (isGray)
        % Set nodes and edges to relevant hemisphere
        ind = center(1);
        isRightHemi = center(2);
        if (isRightHemi)
            nodes   = searchlight.vw_gray.allRightNodes;
            edges   = searchlight.vw_gray.allRightEdges;
        else
            nodes   = searchlight.vw_gray.allLeftNodes;
            edges   = searchlight.vw_gray.allLeftEdges;
        end

        % Compute coordinates of grown disc
        dist        = mrManDist(nodes, edges, double(ind), [1 1 1], 0, growBy);
        discCoords  = nodes(1:3, dist > 0);
        if (isempty(discCoords)), roiCoords = []; return; end
        roiCoords  = [discCoords(2, :); discCoords(1, :); discCoords(3, :)]; % rearranging because of fudged order

        % Recompute hemispheric index to parameter map index
        if (isRightHemi)
            [c center(1)] = intersectCols(searchlight.vw_gray.nodes(1:3, :), ...
                searchlight.vw_gray.allRightNodes(1:3, ind)); % for 3 coordinates
        else
            [c center(1)] = intersectCols(searchlight.vw_gray.nodes(1:3, :), ...
                searchlight.vw_gray.allLeftNodes(1:3, ind)); % for 3 coordinates
        end

    else
        minmax = repmat(center, [1 2]);
        minmax(:, 1) = minmax(:, 1) - growBy;
        minmax(:, 2) = minmax(:, 2) + growBy;
        if (sum(minmax(:, 1) < 1) || sum(minmax(:, 2) > searchlight.dims)), roiCoords = []; return; end % out of bounds

        cubewidth = growBy * 2 + 1;
        xcoords = reshape(repmat(minmax(1, 1):minmax(1, 2), [cubewidth^2 1]), [1 cubewidth^3]);
        ycoords = repmat(reshape(repmat(minmax(2, 1):minmax(2, 2), [cubewidth 1]), [1 cubewidth^2]), [1 cubewidth]);
        zcoords = repmat(minmax(3, 1):minmax(3, 2), [1 cubewidth^2]);

        roiCoords = [xcoords; ycoords; zcoords];
    end
end

function meanAcc = runSVM(searchlight, roiCoords)
    meanAcc = [];
    if (isempty(roiCoords)), return; end

    svmTmp = searchlight.svm;
    if (strcmp(searchlight.searchView, 'inplane'))
        coords = svmTmp.coordsInplane;
    else
        coords = svmTmp.coordsGray;
    end

    [c svmInds] = intersectCols(coords, roiCoords); % for 3 coordinates
    if (isempty(svmInds)), svmInds = 1:0; end
    svmTmp.data = svmTmp.data(:, svmInds);
    svmTmp.data = svmTmp.data(:, ~isnan(svmTmp.data(1,:)));
    svmTmp.voxel = [];

    [meanAcc models] = svmRun(svmTmp, 'options', searchlight.runOptions);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = vscale(F)
%VSCALE   Vertical scale of a CHEBFUN3V.
%   VSCL = VSCALE(F) returns the maximal vertical scale of the components of a
%   CHEBFUN3V object F as determined by evaluating F on a coarse tensor-product
%   grid.

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

for jj = 1:F.nComponents
    vscl(jj,1) = vscale(F.components{jj});
end
out = max(vscl);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [H] = rigidbody(f);

% RIGIDBODY creates the homogenous spatial transformation matrix
% for a 6 parameter rigid-body transformation
%
% Use as
%   [H] = rigidbody(f)
%
% The transformation vector f should contain the
%   x-shift
%   y-shift
%   z-shift
% followed by the
%   pitch (rotation around x-axis)
%   roll  (rotation around y-axis)
%   yaw   (rotation around z-axis)

% Copyright (C) 2000-2005, Robert Oostenveld
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% $Log: rigidbody.m,v $
% Revision 1.1  2009/01/30 04:02:11  arno
% *** empty log message ***
%
% Revision 1.4  2006/04/13 10:37:38  roboos
% added a ; to the end of a line
%
% Revision 1.3  2005/08/15 08:15:32  roboos
% reimplemented the rotate function, which contained an error (the error is in the AIR technical reference)
% changed all functions to be dependent on the rotate, translate and scale function
% all functions now behave consistenly, which also means that they are not compleetly backward compatible w.r.t. the order of the rotations
%
% Revision 1.2  2004/05/19 09:57:07  roberto
% added GPL copyright statement, added CVS log item
%

% compute the homogenous transformation matrix for the translation
T = translate(f([1 2 3]));

% compute the homogenous transformation matrix for the rotation
R = rotate(f([4 5 6]));

% compute the homogenous transformation matrix for the combination
H = T*R;

author: PatternRecognition

repo: OpenBMI


text:
function [sz]=tt_qsize(tt,s)

% For a QTT decomposition _tt_,
% k-th core of which is indexed by
% i^{k}_{1} ,..., i^{k}_{s},
% k=1,...,d,
% returns the array of mode lengths
% sz = [n^{1}_{1} ,..., n^{1}_{s};
%                 ,...,
%       n^{d}_{1} ,..., n^{d}_{s}]
%
% April 26, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%

d=size(tt,1);
sz=zeros(d,s);

for k=1:d
	szk=size(tt{k});
	szk=[szk,ones(1,s)];
	sz(k,1:s)=szk(1:s);
end

return
end

author: oseledets

repo: TT-Toolbox


text:
function triangulation_test31 ( )

%*****************************************************************************80
%
%% TEST31 tests VORONOI_POLYGON_AREA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  neighbor_num = 4;
  node_num = 5;

  area_correct = 0.5;
  center = 5;
  neighbor_index = [ 1, 2, 3, 4 ];
  node_xy = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    1.0, 1.0; ...
    0.0, 1.0; ...
    0.5, 0.5 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST31\n' );
  fprintf ( 1, '  VORONOI_POLYGON_AREA computes the area of\n' );
  fprintf ( 1, '  a finite Voronoi polygon.\n' );

  area = voronoi_polygon_area ( center, neighbor_num, neighbor_index, ...
    node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The computed area is %f\n', area );
  fprintf ( 1, '  The correct area is  %f\n', area_correct );

  return
end