text:
%PRINT2EPS  Prints figures to eps with improved line styles
%
% Examples:
%   print2eps filename
%   print2eps(filename, fig_handle)
%   print2eps(filename, fig_handle, options)
%
% This function saves a figure as an eps file, with two improvements over
% MATLAB's print command. First, it improves the line style, making dashed
% lines more like those on screen and giving grid lines their own dotted
% style. Secondly, it substitutes original font names back into the eps
% file, where these have been changed by MATLAB, for up to 11 different
% fonts.
%
%IN:
%   filename - string containing the name (optionally including full or
%              relative path) of the file the figure is to be saved as. A
%              ".eps" extension is added if not there already. If a path is
%              not specified, the figure is saved in the current directory.
%   fig_handle - The handle of the figure to be saved. Default: gcf.
%   options - Additional parameter strings to be passed to print.

% Copyright (C) Oliver Woodford 2008-2013

% The idea of editing the EPS file to change line styles comes from Jiro
% Doke's FIXPSLINESTYLE (fex id: 17928)
% The idea of changing dash length with line width came from comments on
% fex id: 5743, but the implementation is mine :)

% 14/11/2011: Fix a MATLAB bug rendering black or white text incorrectly.
%             Thanks to Mathieu Morlighem for reporting the issue and
%             obtaining a fix from TMW.
% 08/12/11: Added ability to correct fonts. Several people have requested
%           this at one time or another, and also pointed me to printeps
%           (fex id: 7501), so thank you to them. My implementation (which
%           was not inspired by printeps - I'd already had the idea for my
%           approach) goes slightly further in that it allows multiple
%           fonts to be swapped.
% 14/12/11: Fix bug affecting font names containing spaces. Thanks to David
%           Szwer for reporting the issue.
% 25/01/12: Add a font not to be swapped. Thanks to Anna Rafferty and Adam
%           Jackson for reporting the issue. Also fix a bug whereby using a
%           font alias can lead to another font being swapped in.
% 10/04/12: Make the font swapping case insensitive.
% 26/10/12: Set PaperOrientation to portrait. Thanks to Michael Watts for
%           reporting the issue.
% 26/10/12: Fix issue to do with swapping fonts changing other fonts and
%           sizes we don't want, due to listeners. Thanks to Malcolm Hudson
%           for reporting the issue.
% 22/03/13: Extend font swapping to axes labels. Thanks to Rasmus Ischebeck
%           for reporting the issue.

function print2eps(name, fig, varargin)
options = {'-depsc2'};
if nargin < 2
    fig = gcf;
elseif nargin > 2
    options = [options varargin];
end
% Construct the filename
if numel(name) < 5 || ~strcmpi(name(end-3:end), '.eps')
    name = [name '.eps']; % Add the missing extension
end
% Find all the used fonts in the figure
font_handles = findall(fig, '-property', 'FontName');
fonts = get(font_handles, 'FontName');
if ~iscell(fonts)
    fonts = {fonts};
end
% Map supported font aliases onto the correct name
fontsl = lower(fonts);
for a = 1:numel(fonts)
    f = fontsl{a};
    f(f==' ') = [];
    switch f
        case {'times', 'timesnewroman', 'times-roman'}
            fontsl{a} = 'times-roman';
        case {'arial', 'helvetica'}
            fontsl{a} = 'helvetica';
        case {'newcenturyschoolbook', 'newcenturyschlbk'}
            fontsl{a} = 'newcenturyschlbk';
        otherwise
    end
end
fontslu = unique(fontsl);
% Determine the font swap table
matlab_fonts = {'Helvetica', 'Times-Roman', 'Palatino', 'Bookman', 'Helvetica-Narrow', 'Symbol', ...
                'AvantGarde', 'NewCenturySchlbk', 'Courier', 'ZapfChancery', 'ZapfDingbats'};
matlab_fontsl = lower(matlab_fonts);
require_swap = find(~ismember(fontslu, matlab_fontsl));
unused_fonts = find(~ismember(matlab_fontsl, fontslu));
font_swap = cell(3, min(numel(require_swap), numel(unused_fonts)));
fonts_new = fonts;
for a = 1:size(font_swap, 2)
    font_swap{1,a} = find(strcmp(fontslu{require_swap(a)}, fontsl));
    font_swap{2,a} = matlab_fonts{unused_fonts(a)};
    font_swap{3,a} = fonts{font_swap{1,end}(1)};
    fonts_new(font_swap{1,a}) = {font_swap{2,a}};
end
% Swap the fonts
if ~isempty(font_swap)
    fonts_size = get(font_handles, 'FontSize');
    if iscell(fonts_size)
        fonts_size = cell2mat(fonts_size);
    end
    M = false(size(font_handles));
    % Loop because some changes may not stick first time, due to listeners
    c = 0;
    update = zeros(1000, 1);
    for b = 1:10 % Limit number of loops to avoid infinite loop case
        for a = 1:numel(M)
            M(a) = ~isequal(get(font_handles(a), 'FontName'), fonts_new{a}) || ~isequal(get(font_handles(a), 'FontSize'), fonts_size(a));
            if M(a)
                set(font_handles(a), 'FontName', fonts_new{a}, 'FontSize', fonts_size(a));
                c = c + 1;
                update(c) = a;
            end
        end
        if ~any(M)
            break;
        end
    end
    % Compute the order to revert fonts later, without the need of a loop
    [update, M] = unique(update(1:c));
    [M, M] = sort(M);
    update = reshape(update(M), 1, []);
end
% Set paper size
old_pos_mode = get(fig, 'PaperPositionMode');
old_orientation = get(fig, 'PaperOrientation');
set(fig, 'PaperPositionMode', 'auto', 'PaperOrientation', 'portrait');
% MATLAB bug fix - black and white text can come out inverted sometimes
% Find the white and black text
white_text_handles = findobj(fig, 'Type', 'text');
M = get(white_text_handles, 'Color');
if iscell(M)
    M = cell2mat(M);
end
M = sum(M, 2);
black_text_handles = white_text_handles(M == 0);
white_text_handles = white_text_handles(M == 3);
% Set the font colors slightly off their correct values
set(black_text_handles, 'Color', [0 0 0] + eps);
set(white_text_handles, 'Color', [1 1 1] - eps);
% Print to eps file
print(fig, options{:}, name);
% Reset the font colors
set(black_text_handles, 'Color', [0 0 0]);
set(white_text_handles, 'Color', [1 1 1]);
% Reset paper size
set(fig, 'PaperPositionMode', old_pos_mode, 'PaperOrientation', old_orientation);
% Correct the fonts
if ~isempty(font_swap)
    % Reset the font names in the figure
    for a = update
        set(font_handles(a), 'FontName', fonts{a}, 'FontSize', fonts_size(a));
    end
    % Replace the font names in the eps file
    font_swap = font_swap(2:3,:);
    try
        swap_fonts(name, font_swap{:});
    catch
        warning('swap_fonts() failed. This is usually because the figure contains a large number of patch objects. Consider exporting to a bitmap format in this case.');
        return
    end
end
% Fix the line styles
try
    fix_lines(name);
catch
    warning('fix_lines() failed. This is usually because the figure contains a large number of patch objects. Consider exporting to a bitmap format in this case.');
end
return

function swap_fonts(fname, varargin)
% Read in the file
fh = fopen(fname, 'r');
if fh == -1
    error('File %s not found.', fname);
end
try
    fstrm = fread(fh, '*char')';
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);

% Replace the font names
for a = 1:2:numel(varargin)
    fstrm = regexprep(fstrm, [varargin{a} '-?[a-zA-Z]*\>'], varargin{a+1}(~isspace(varargin{a+1})));
end

% Write out the updated file
fh = fopen(fname, 'w');
if fh == -1
    error('Unable to open %s for writing.', fname2);
end
try
    fwrite(fh, fstrm, 'char*1');
catch ex
    fclose(fh);
    rethrow(ex);
end
fclose(fh);
return

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function DeleteDataFiles(varargin)

% Syntax:
%
%   DeleteDataFiles()
%   DeleteDataFiles(dirname)
%   DeleteDataFiles(dirname, datafiles0)
%   DeleteDataFiles(dirname, format)
%   DeleteDataFiles(dirname, datafiles0, options)
%   DeleteDataFiles(dirname, format, options)
%
% Description:
%
%   Delete all .<format extension> files in group folder. DeleteDataFiles will find all
%   .<format extension> data acquisition files in the group folder. If dirname is not supplied it'll
%   treat the current working directory as the group folder.
%
% Examples:
%
%   1. Delete all .snirf files in the Homer3 Examples directory
%
%       DeleteDataFiles('C:\jdubb\workspaces\Homer3\DataTree\AcquiredData\Snirf\Examples', '.snirf')
%
%   2. Delete all .nirs files files in the current group folder
%
%       DeleteDataFiles(pwd, '.nirs)
%
%

global supportedFormats
if isempty(supportedFormats)
    supportedFormats = {
    '.snirf',0;
    '.nirs',0;
    };
end

% Set argument defaults
dirname = filesepStandard(pwd);
format = supportedFormats{1};
datafiles0 = [];
options = 'delete';

% Parse arguments
if nargin==1
    dirname = filesepStandard(varargin{1},'full');
elseif nargin==2
    dirname = filesepStandard(varargin{1},'full');
    if ischar(varargin{2})
        format = varargin{2};
    else
        datafiles0 = varargin{2};
    end
elseif nargin==3
    dirname = filesepStandard(varargin{1}, 'full');
    if ischar(varargin{2})
        format = varargin{2};
    else
        datafiles0 = varargin{2};
    end
    options = varargin{3};
end
if isempty(datafiles0)
    datafiles0 = DataFilesClass(dirname, format, 'standalone').files;
end

% Get final list of data files
datafiles = mydir(dirname);
if iscell(datafiles0)
    for ii = 1:length(datafiles0)
        datafiles(ii) = mydir([dirname, datafiles0{ii}]);
    end
elseif ischar(datafiles0)
    datafiles = mydir([dirname, datafiles0]);
elseif isa(datafiles0, 'FileClass')
    datafiles = datafiles0;
end

% Delete data files
for ii = 1:length(datafiles)
    if datafiles(ii).isdir
        continue;
    end
    if strcmp(options, 'delete')
        fprintf('Deleting %s\n', [datafiles(ii).rootdir, '/', datafiles(ii).name]);
        delete([datafiles(ii).rootdir, '/', datafiles(ii).name]);
        delete([datafiles(ii).rootdir, '/*_events.tsv']);
    elseif strcmp(options, 'move')
        fprintf('Moving %s to %s\n', [datafiles(ii).rootdir, '/', datafiles(ii).name], [datafiles(ii).rootdir, '/', datafiles(ii).name, '.old']);
        movefile([datafiles(ii).rootdir, '/', datafiles(ii).name], [datafiles(ii).rootdir, '/', datafiles(ii).name, '.old']);
    elseif strcmp(options, 'restore')
        [pname, fname] = fileparts(datafiles(ii).name);
        fprintf('Restoring %s to %s\n', [datafiles(ii).rootdir, '/', datafiles(ii).name], [datafiles(ii).rootdir, '/', pname, '/', fname]);
        movefile([datafiles(ii).rootdir, '/', datafiles(ii).name], [datafiles(ii).rootdir, '/', pname, '/', fname]);
    end
    pause(0.25);
end

author: BUNPC


text:
function geometry_test0765 ( )

%*****************************************************************************80
%
%% TEST0765 tests POLYGON_AREA_2D_2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  test_num = 2;
  area_exact_test = [ 2.0, 6.0 ];
  n_test = [ 4, 8 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0765\n' );
  fprintf ( 1, '  For a polygon in 2D:\n' );
  fprintf ( 1, '  POLYGON_AREA_2D_2 computes the area.\n' );

  for test = 1 : 2

    n = n_test(test);
    area_exact = area_exact_test(test);

    if ( test == 1 )

      v = [ ...
        1.0, 0.0; ...
        2.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 1.0 ]';

    elseif ( test == 2 )

      v = [ ...
        0.0, 0.0; ...
        3.0, 0.0; ...
        3.0, 3.0; ...
        2.0, 3.0; ...
        2.0, 1.0; ...
        1.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 2.0 ]';

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of polygonal vertices = %d\n', n );

    r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices:' );

    area = polygon_area_2d_2 ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Exact area is        %f\n', area_exact );
    fprintf ( 1, '  The computed area is %f\n', area );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
k=1/(2*sqrt(log(2)));
x=0:30/1999:30;
a1=1.0; t1=5; w1=1.3;
y1=a1*exp(-((x-t1)/(w1*k)).^2);
a2=1.0; t2=15.0; w2=1.3;
y2=a2*(1+4*((x-t2)/w2).^2).^(-1);
a3=1.0; t3=25;
y3=a3*ones(size(x))./(1+exp(-3.0*(x-t3)));
y=y1+y2+y3;
noise=rand(size(y));
y=y+(noise-0.5)*0.01;

figure
subplot(3,1,1)
plot(x,y);
title('noisy data');

subplot(3,1,2)
plot(x(1:end-1)+(x(2)-x(1)),diff(y)/(x(2)-x(1)));
title('diff function');

subplot(3,1,3)
plot(x,filter(-smooth_diff(10),1,y)/(x(2)-x(1)));
title('smoothed differentiation');

author: Sable

repo: mcbench-benchmarks


text:
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 9 )

    zemu = pi /  2.0;

    aj(1:m) = 0.0;

    bj(1:m) = 0.5;

  end

  return
end


text:
filename='Gripping_quad_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.75;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function minmax_install
% function minmax_install
% Installation by building the C-mex files for min/max selection package
%
% Author Bruno Luong <brunoluong@yahoo.com>
% Last update: 29-Jun-2009 built inplace functions

thisfile = mfilename('fullpath');
path = fileparts(thisfile);
oldpath = cd(path);

arch=computer('arch');
mexopts = {'-v' '-O' ['-' arch]};
% 64-bit platform
if ~isempty(strfind(computer(),'64'))
    mexopts(end+1) = {'-largeArrayDims'};
end

if ispc() && datenum(version('-date')) < datenum('January 11, 2014')
    compiler = getmexopts('COMPILER');
    islcc = strcmpi(compiler,'lcc');
    % Define the C-symbol for LCC compiler
    if islcc
        mexopts(end+1) = {'-D_LCC'};
    end
end

% Internal representation of mxArray
try
    buildInternal_mxArrayDef('Internal_mxArray.h');
catch
    if ispc()
        cpcmd = 'copy';
    else
        cpcmd = ' cp';
    end
    cmd = [cpcmd ' Internal_mxArray_2010B.h Internal_mxArray.h'];
    system(cmd);
end

% Inplace tool
mex(mexopts{:},'inplacecolumnmex.c');
mex(mexopts{:},'releaseinplace.c');

% Mex MIN/MAX functions
mex(mexopts{:},'minkmex.c');
mex(mexopts{:},'maxkmex.c');

cd(oldpath);

end

author: andrewssobral

repo: lrslibrary


text:
function [Words,Simbols] = huffcodes2bin(huffcodes)
%HUFFCODES2BIN   Convert huffcodes to binary representation
%   [W,S] = HUFFCODES2BIN(HC) returns the Huffman representation returned
%      by the function NORM2HUFF to the corresponding binary string of '0' and '1'
%   W is a cell array of binary strings,
%   S is a double array containing the corresponding simbols
%
%   REMARK: the first bit of each binary strings B is B(end)
%
%   See also HUFF2NORM, NORM2HUFF

%   $Author: Giuseppe Ridino' $
%   $Revision: 1.0 $  $Date: 25-May-2004 14:26:00 $

% get number of simbols
Nsimbols = nnz(huffcodes);

% initialize output
Words = cell(Nsimbols,1);
% Simbols = cell(Nsimbols,1);

% gte code words
CodeWords = find(not(huffcodes==0));

% convert words
for index=1:Nsimbols,
	Words{index} = Double2BinStr(CodeWords(index));
end

% return simbols
Simbols = full(huffcodes(CodeWords));

% #############################################################
function BinString = Double2BinStr(Word)
BinString = dec2bin(Word);
if not(BinString(1)=='1'),
	error('wrong codeword');        % most significat bit must always be "1"
else,
	BinString = BinString(2:end);   % remove stop bit
end


text:
function val = gamma_multivariate_ln(x,p);
% function val = gamma_multivariate_ln(x,p);
%
% x: array(1,K)
% p: scalor
%
% x must be more than (p-1)/2
% x should be more than p/2
%
% Gamma_p(x) = pi^(p(p-1)/4) prod_(j=1)^p Gamma(x+(1-j)/2)
% log Gamma_p(x) = p(p-1)/4 log pi + sum_(j=1)^p log Gamma(x+(1-j)/2)

K = length(x);
gammaln_val = gammaln(repmat(x,p,1)+0.5*(1-repmat([1:p]',1,K))); % p by K
val = p*(p-1)*0.25 * log(pi) + sum(gammaln_val,1);

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

repo: REST


text:
function [tap] = hanning(n, str)

%HANNING   Hanning window.
%   HANNING(N) returns the N-point symmetric Hanning window in a column
%   vector.  Note that the first and last zero-weighted window samples
%   are not included.
%
%   HANNING(N,'symmetric') returns the same result as HANNING(N).
%
%   HANNING(N,'periodic') returns the N-point periodic Hanning window,
%   and includes the first zero-weighted window sample.
%
%   NOTE: Use the HANN function to get a Hanning window which has the
%          first and last zero-weighted samples.
%
%   See also BARTLETT, BLACKMAN, BOXCAR, CHEBWIN, HAMMING, HANN, KAISER
%   and TRIANG.
%
%   This is a drop-in replacement to bypass the signal processing toolbox

% Copyright (c) 2010, Jan-Mathijs Schoffelen, DCCN Nijmegen
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin==1,
  str = 'symmetric';
end

switch str,
case 'periodic'
   % Includes the first zero sample
   tap = [0; hanningX(n-1)];
case 'symmetric'
   % Does not include the first and last zero sample
   tap = hanningX(n);
end

function tap = hanningX(n)

% compute taper
N   = n+1;
tap = 0.5*(1-cos((2*pi*(1:n))./N))';

% make symmetric
halfn = floor(n/2);
tap( (n+1-halfn):n ) = flipud(tap(1:halfn));

author: spm


text:
function [TB] = B2TB(B)
% Convert computery things from bytes to terabytes.
% Chad A. Greene 2012
TB = B*2^-40 ;


text:
function varargout = EnvironmentalSelection(varargin)
% The environmental selection of PM-MOEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Delete duplicated solutions and non-dominated sorting
    if nargin == 3
        [PopDec,PopObj,N] = deal(varargin{:});
        [PopDec,uni] = unique(PopDec,'rows','stable');
        PopObj = PopObj(uni,:);
        [FrontNo,MaxFNo] = NDSort(PopObj,N);
    elseif nargin == 4
        [Population,Dec,Mask,N] = deal(varargin{:});
        PopObj = Population.objs;
        [FrontNo,MaxFNo] = NDSort(PopObj,Population.cons,N);
    end
    Next = FrontNo <= MaxFNo;

    %% Truncate the solutions in the last front
    Last = find(FrontNo==MaxFNo);
    if nargin == 3
        Del = Truncation(double(PopDec(Last,:)),sum(Next)-N);
    elseif nargin == 4
        Del = Truncation(PopObj(Last,:),sum(Next)-N);
    end
    Next(Last(Del)) = false;

    %% Population for next generation
    FrontNo = FrontNo(Next);
    if nargin == 3
        PopDec    = PopDec(Next,:);
        PopObj    = PopObj(Next,:);
        varargout = {PopDec,PopObj,FrontNo};
    elseif nargin == 4
        Population = Population(Next);
        Dec        = Dec(Next,:);
        Mask       = Mask(Next,:);
        varargout  = {Population,Dec,Mask,FrontNo};
    end
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef PoissonDistributionX < ProbabilityDistributionX
% PoissonDistributionX class
%
% Summary of PoissonDistributionX:
% This is a class implementation of a Poisson distribution
%
% PoissonDistributionX Properties:
%   + Mean - Mean of the Poisson distribution
%   + Covar - Variance of the Poisson distribution
%
% PoissonDistributionX Methods:
%   + PoissonDistributionX  - Constructor method
%   + reset - Reset the distribution with a new number of random variables
%   + random - Draw random samples from the multivariate normal distribution
%   + pdf - Evaluate the density of the multivariate normal distribution
%
% (+) denotes puplic properties/methods
%
% See also TransitionModelX, MeasurementModelX and ControlModelX template classes

    properties (Dependent)
        % Mean: scalar
        %   Mean of the Poisson distribution
        Mean

        % Covar: scalar
        %   Variance of the Poisson distribution
        Covar
    end

    properties (Access = private)
        % Underlying poisson distribution
        pd_ = makedist('Poisson');
    end

    methods
        function this = PoissonDistributionX(varargin)
        % PoissonDistributionX Construct a Poisson distribution
        %
        % Parameters
        % ----------
        % mean: non-negative scalar
        %   Mean of the Poisson distribution
            this.reset(varargin{:});
        end

        function samples = random(this, numSamples)
        % Samples Draw random samples from the Gaussian
        %
        % Parameters
        % ----------
        % numSamples: scalar
        %   The number of samples to be drawn from the mixture
        %
        % Returns
        % -------
        % samples: (NumVariables x numSamples) matrix
        %   The set of samples drawn from the mixture.

            assert(numSamples >= 1);

            % Incorporate actual mean
            samples = this.pd_.random(numSamples,1)';
        end

        function prob = pdf(this, samples)
        % pdf Evaluate the density of the poisson distribution
        %   PROB = PDF(OBJ, SAMPLES) returns the density of the poisson distribution,
        %   evaluated at each column of SAMPLES.
        %
            prob = this.pd_.pdf(samples);
        end

        function reset(this, varargin)
            if(nargin >= 2)
                mean = varargin{1};

                % Reset number of random variables. This will also take care of
                % input validation.
                this.NumVariables = size(mean,1);

                % Re-initialize the mean and covariance
                this.Mean = mean;
            end
        end
    end

    methods
        function meanValue = get.Mean(this)
        %get.Mean Getter for Mean property
            meanValue = this.pd_.mean();
        end

        function set.Mean(this, meanValue)
        %set.Mean - Setter for Mean property
            this.pd_ = makedist('Poisson','lambda',  meanValue);
        end

        function covariance = get.Covar(this)
        %get.Covariance Getter for Covariance property
            covariance = this.pd_.var();
        end
    end

end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
% MTSP_GA Multiple Traveling Salesmen Problem (M-TSP) Genetic Algorithm (GA)
%   Finds a (near) optimal solution to the M-TSP by setting up a GA to search
%   for the shortest route (least distance needed for the salesmen to travel
%   to each city exactly once and return to their starting locations)
%
% Summary:
%     1. Each salesman travels to a unique set of cities and completes the
%        route by returning to the city he started from
%     2. Each city is visited by exactly one salesman
%
% Input:
%     XY (float) is an Nx2 matrix of city locations, where N is the number of cities
%     DMAT (float) is an NxN matrix of city-to-city distances or costs
%     NSALESMEN (scalar integer) is the number of salesmen to visit the cities
%     MINTOUR (scalar integer) is the minimum tour length for any of the salesmen
%     POPSIZE (scalar integer) is the size of the population (should be divisible by 8)
%     NUMITER (scalar integer) is the number of desired iterations for the algorithm to run
%     SHOWPROG (scalar logical) shows the GA progress if true
%     SHOWRESULT (scalar logical) shows the GA results if true
%
% Output:
%     OPTROUTE (integer array) is the best route found by the algorithm
%     OPTBREAK (integer array) is the list of route break points (these specify the indices
%         into the route used to obtain the individual salesman routes)
%     MINDIST (scalar float) is the total distance traveled by the salesmen
%
% Route/Breakpoint Details:
%     If there are 10 cities and 3 salesmen, a possible route/break
%     combination might be: rte = [5 6 9 1 4 2 8 10 3 7], brks = [3 7]
%     Taken together, these represent the solution [5 6 9][1 4 2 8][10 3 7],
%     which designates the routes for the 3 salesmen as follows:
%         . Salesman 1 travels from city 5 to 6 to 9 and back to 5
%         . Salesman 2 travels from city 1 to 4 to 2 to 8 and back to 1
%         . Salesman 3 travels from city 10 to 3 to 7 and back to 10
%
% Example:
%     n = 35;
%     xy = 10*rand(n,2);
%     nSalesmen = 5;
%     minTour = 3;
%     popSize = 80;
%     numIter = 5e3;
%     a = meshgrid(1:n);
%     dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),n,n);
%     [optRoute,optBreak,minDist] = mtsp_ga(xy,dmat,nSalesmen,minTour, ...
%         popSize,numIter,1,1);
%
% Example:
%     n = 50;
%     phi = (sqrt(5)-1)/2;
%     theta = 2*pi*phi*(0:n-1);
%     rho = (1:n).^phi;
%     [x,y] = pol2cart(theta(:),rho(:));
%     xy = 10*([x y]-min([x;y]))/(max([x;y])-min([x;y]));
%     nSalesmen = 5;
%     minTour = 3;
%     popSize = 80;
%     numIter = 1e4;
%     a = meshgrid(1:n);
%     dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),n,n);
%     [optRoute,optBreak,minDist] = mtsp_ga(xy,dmat,nSalesmen,minTour, ...
%         popSize,numIter,1,1);
%
% Example:
%     n = 35;
%     xyz = 10*rand(n,3);
%     nSalesmen = 5;
%     minTour = 3;
%     popSize = 80;
%     numIter = 5e3;
%     a = meshgrid(1:n);
%     dmat = reshape(sqrt(sum((xyz(a,:)-xyz(a',:)).^2,2)),n,n);
%     [optRoute,optBreak,minDist] = mtsp_ga(xyz,dmat,nSalesmen,minTour, ...
%         popSize,numIter,1,1);
%
% See also: tsp_ga, mtspf_ga, mtspo_ga, mtspof_ga, mtspofs_ga, mtspv_ga, distmat
%
% Author: Joseph Kirk
% Email: jdkirk630@gmail.com
% Release: 1.5
% Release Date: 11/07/11
function varargout = mtsp_ga(xy,dmat,nSalesmen,minTour,popSize,numIter,showProg,showResult)

% Process Inputs and Initialize Defaults
nargs = 8;
for k = nargin:nargs-1
    switch k
        case 0
            xy = 10*rand(40,2);
        case 1
            N = size(xy,1);
            a = meshgrid(1:N);
            dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),N,N);
        case 2
            nSalesmen = 5;
        case 3
            minTour = 3;
        case 4
            popSize = 80;
        case 5
            numIter = 5e3;
        case 6
            showProg = 1;
        case 7
            showResult = 1;
        otherwise
    end
end

% Verify Inputs
[N,dims] = size(xy);
[nr,nc] = size(dmat);
if N ~= nr || N ~= nc
    error('Invalid XY or DMAT inputs!')
end
n = N;

% Sanity Checks
nSalesmen = max(1,min(n,round(real(nSalesmen(1)))));
minTour = max(1,min(floor(n/nSalesmen),round(real(minTour(1)))));
popSize = max(8,8*ceil(popSize(1)/8));
numIter = max(1,round(real(numIter(1))));
showProg = logical(showProg(1));
showResult = logical(showResult(1));

% Initializations for Route Break Point Selection
nBreaks = nSalesmen-1;
dof = n - minTour*nSalesmen;          % degrees of freedom
addto = ones(1,dof+1);
for k = 2:nBreaks
    addto = cumsum(addto);
end
cumProb = cumsum(addto)/sum(addto);

% Initialize the Populations
popRoute = zeros(popSize,n);         % population of routes
popBreak = zeros(popSize,nBreaks);   % population of breaks
popRoute(1,:) = (1:n);
popBreak(1,:) = rand_breaks();
for k = 2:popSize
    popRoute(k,:) = randperm(n);
    popBreak(k,:) = rand_breaks();
end

% Select the Colors for the Plotted Routes
pclr = ~get(0,'DefaultAxesColor');
clr = [1 0 0; 0 0 1; 0.67 0 1; 0 1 0; 1 0.5 0];
if nSalesmen > 5
    clr = hsv(nSalesmen);
end

% Run the GA
globalMin = Inf;
totalDist = zeros(1,popSize);
distHistory = zeros(1,numIter);
tmpPopRoute = zeros(8,n);
tmpPopBreak = zeros(8,nBreaks);
newPopRoute = zeros(popSize,n);
newPopBreak = zeros(popSize,nBreaks);
if showProg
    pfig = figure('Name','MTSP_GA | Current Best Solution','Numbertitle','off');
end
for iter = 1:numIter
    % Evaluate Members of the Population
    for p = 1:popSize
        d = 0;
        pRoute = popRoute(p,:);
        pBreak = popBreak(p,:);
        rng = [[1 pBreak+1];[pBreak n]]';
        for s = 1:nSalesmen
            d = d + dmat(pRoute(rng(s,2)),pRoute(rng(s,1)));
            for k = rng(s,1):rng(s,2)-1
                d = d + dmat(pRoute(k),pRoute(k+1));
            end
        end
        totalDist(p) = d;
    end

    % Find the Best Route in the Population
    [minDist,index] = min(totalDist);
    distHistory(iter) = minDist;
    if minDist < globalMin
        globalMin = minDist;
        optRoute = popRoute(index,:);
        optBreak = popBreak(index,:);
        rng = [[1 optBreak+1];[optBreak n]]';
        if showProg
            % Plot the Best Route
            figure(pfig);
            for s = 1:nSalesmen
                rte = optRoute([rng(s,1):rng(s,2) rng(s,1)]);
                if dims > 2, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));
                else plot(xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); end
                title(sprintf('Total Distance = %1.4f, Iteration = %d',minDist,iter));
                hold on
            end
            hold off
        end
    end

    % Genetic Algorithm Operators
    randomOrder = randperm(popSize);
    for p = 8:8:popSize
        rtes = popRoute(randomOrder(p-7:p),:);
        brks = popBreak(randomOrder(p-7:p),:);
        dists = totalDist(randomOrder(p-7:p));
        [ignore,idx] = min(dists); %#ok
        bestOf8Route = rtes(idx,:);
        bestOf8Break = brks(idx,:);
        routeInsertionPoints = sort(ceil(n*rand(1,2)));
        I = routeInsertionPoints(1);
        J = routeInsertionPoints(2);
        for k = 1:8 % Generate New Solutions
            tmpPopRoute(k,:) = bestOf8Route;
            tmpPopBreak(k,:) = bestOf8Break;
            switch k
                case 2 % Flip
                    tmpPopRoute(k,I:J) = tmpPopRoute(k,J:-1:I);
                case 3 % Swap
                    tmpPopRoute(k,[I J]) = tmpPopRoute(k,[J I]);
                case 4 % Slide
                    tmpPopRoute(k,I:J) = tmpPopRoute(k,[I+1:J I]);
                case 5 % Modify Breaks
                    tmpPopBreak(k,:) = rand_breaks();
                case 6 % Flip, Modify Breaks
                    tmpPopRoute(k,I:J) = tmpPopRoute(k,J:-1:I);
                    tmpPopBreak(k,:) = rand_breaks();
                case 7 % Swap, Modify Breaks
                    tmpPopRoute(k,[I J]) = tmpPopRoute(k,[J I]);
                    tmpPopBreak(k,:) = rand_breaks();
                case 8 % Slide, Modify Breaks
                    tmpPopRoute(k,I:J) = tmpPopRoute(k,[I+1:J I]);
                    tmpPopBreak(k,:) = rand_breaks();
                otherwise % Do Nothing
            end
        end
        newPopRoute(p-7:p,:) = tmpPopRoute;
        newPopBreak(p-7:p,:) = tmpPopBreak;
    end
    popRoute = newPopRoute;
    popBreak = newPopBreak;
end

if showResult
% Plots
    figure('Name','MTSP_GA | Results','Numbertitle','off');
    subplot(2,2,1);
    if dims > 2, plot3(xy(:,1),xy(:,2),xy(:,3),'.','Color',pclr);
    else plot(xy(:,1),xy(:,2),'.','Color',pclr); end
    title('City Locations');
    subplot(2,2,2);
    imagesc(dmat(optRoute,optRoute));
    title('Distance Matrix');
    subplot(2,2,3);
    rng = [[1 optBreak+1];[optBreak n]]';
    for s = 1:nSalesmen
        rte = optRoute([rng(s,1):rng(s,2) rng(s,1)]);
        if dims > 2, plot3(xy(rte,1),xy(rte,2),xy(rte,3),'.-','Color',clr(s,:));
        else plot(xy(rte,1),xy(rte,2),'.-','Color',clr(s,:)); end
        title(sprintf('Total Distance = %1.4f',minDist));
        hold on;
    end
    subplot(2,2,4);
    plot(distHistory,'b','LineWidth',2);
    title('Best Solution History');
    set(gca,'XLim',[0 numIter+1],'YLim',[0 1.1*max([1 distHistory])]);
end

% Return Outputs
if nargout
    varargout{1} = optRoute;
    varargout{2} = optBreak;
    varargout{3} = minDist;
end

    % Generate Random Set of Break Points
    function breaks = rand_breaks()
        if minTour == 1 % No Constraints on Breaks
            tmpBreaks = randperm(n-1);
            breaks = sort(tmpBreaks(1:nBreaks));
        else % Force Breaks to be at Least the Minimum Tour Length
            nAdjust = find(rand < cumProb,1)-1;
            spaces = ceil(nBreaks*rand(1,nAdjust));
            adjust = zeros(1,nBreaks);
            for kk = 1:nBreaks
                adjust(kk) = sum(spaces == kk);
            end
            breaks = minTour*(1:nBreaks) + cumsum(adjust);
        end
    end
end


text:
% initialize VAR (once)
clear all;close all;clc
% clearvars -except 'CellResp' 'const'; clc

data_dir = GetCurrentDataDir();
save_dir = GetCurrentDataDir();
i_fish = 1;

%% load data
disp(['load fish ' num2str(i_fish)]);

filename = ['Data_F' num2str(i_fish) '.mat'];
load(fullfile(data_dir,filename),'data');
% [CellResp,data,dimCR] = LoadFileFromParts(data_dir,filename,'CellResp');

names = fieldnames(data); % cell of strings
for i = 1:length(names),
    eval([names{i} '= data.',names{i} ';']);
end

%% Initialize VAR
VAR = [];
%%
i_ClusGroup = 1;
VAR(i_fish).ClusGroupName{i_ClusGroup} = 'selection';

i_Cluster = 1;
cIX_abs = data.absIX(1:10:end);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = '1/10 of validated cells';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

%
i_ClusGroup = 2;
VAR(i_fish).ClusGroupName{i_ClusGroup} = 'init';

i_Cluster = 1;
cIX_abs = (1:data.numcell_full)';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = 'all raw cells (ROI''s)';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

i_Cluster = 2;
cIX_abs = data.absIX;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = 'validated cells';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

%%
i_fish = 2;
i_ClusGroup = 1;
i_Cluster = 1;
cIX_abs = (1:data.numcell_full)';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = 'all raw cells (ROI''s)';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

i_ClusGroup = 1;
i_Cluster = 2;
cIX_abs = data.absIX;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = 'validated cells';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

i_ClusGroup = 2;
i_Cluster = 1;
cIX_abs = data.absIX(1:10:end);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).name = '1/10 of validated cells';
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).cIX_abs = cIX_abs;
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).gIX = ones(length(cIX_abs),1);
VAR(i_fish).ClusGroup{i_ClusGroup}(i_Cluster).numK = 1;

%%
save(fullfile(data_dir,'VAR_current.mat'),'VAR');

%%

% varience/std for reps for each cell
% if i_fish==2 || i_fish==3 || i_fish==6,
%     period_real = period/2;
% else
%     period_real = period;
% end
% nrep_real = floor((size(CR,2)-shift)/period_real);
% while period_real*nrep_real+shift>size(CR,2),
%     nrep_real = nrep_real-1;
% end
% CRZ_3D = reshape(CRZ(:,1+shift:period_real*nrep_real+shift),nCells,period_real,[]);
%% updated method, weighing both std between each rep and (summed with) std of 1st half & 2nd half of experiment - 1/8/15
% CRZ = CONST.M_array.CellResp;
% if i_fish==2 || i_fish==3 || i_fish==6,
%     period_real = CONST.M_array.period/2;
% else
%     period_real = CONST.M_array.period;
% end
% CRZ_3D = reshape(CRZ,size(CRZ,1),period_real,[]);
% divide = round(size(CRZ_3D,3)/2);
% CRZ_std1 = std(CRZ_3D(:,:,1:divide),0,3);
% CRZ_std2 = std(CRZ_3D(:,:,divide+1:end),0,3);
% temp1 = mean(CRZ_std1,2);
% temp2 = mean(CRZ_std2,2);
%
% temp12 = horzcat(temp1,temp2);
% temp = mean(temp12,2)+std(temp12,0,2);
% [~,I] = sort(temp);
% M = temp(I);
% figure;plot(M)
%
% figure;imagesc(CRZ(I,:))
%
% nCells = size(CRZ,1);

%% find low variance / stimulus-locked cells
CRZ_std = std(CRZ_3D,0,3);
temp = mean(CRZ_std,2);

% find mean-std thres: 0.5
[~,I] = sort(temp);
M = temp(I);
figure;plot(M)
%%
i_last = length(VAR(i_fish).Class);
M_perc = [0.025,0.1,0.3];
for j = 1:length(M_perc);
    thres = M(round(nCells*M_perc(j)));
    cIX = find(temp<thres);
    i = j+i_last;
    VAR(i_fish).Class(i).round = 0;
    VAR(i_fish).Class(i).name = ['perc < ' num2str(M_perc(j)*100) '%'];
    %     VAR(i_fish).Class(i).notes = ['perc < ' num2str(M_perc(j)*100) '%'];
    VAR(i_fish).Class(i).cIX = cIX;
    VAR(i_fish).Class(i).gIX = ones(length(cIX),1);
    VAR(i_fish).Class(i).numel = length(cIX);
    VAR(i_fish).Class(i).numK = 1;
    VAR(i_fish).Class(i).datatype = 'std';
end

%% shift CR?
% shift = 161;
% nrep=floor(size(CR,2)/period)-1;
%
% skiplist=[];
% IX_rep=setdiff(1:nrep, skiplist);
% IX=zeros(period*length(IX_rep),1);
%
% for i=1:length(IX_rep)
%     IX(period*(i-1)+1:period*i)=period*(IX_rep(i)-1)+1+shift:period*(IX_rep(i))+shift;
% end
%
% % cell_resp=cell_resp(:,1:nrep*period);
% CRA=mean(reshape(CR(:,IX),[nCells,period,length(IX_rep)]),3);
% CRAZ = zscore(CRA')';

author: xiuyechen

repo: FishExplorer


text:
% SP_EXTERIOR_DERIVAITVE: computes the exterior derivative as a matrix with size
%  given by the dimension of two consecutive spaces in the De Rham sequence.
%
%   diff_op = sp_exterior_derivative (space1, space2);
%
% INPUT:
%
%   space1:  domain space of the exterior derivative (number of columns)
%   space2:  image space of the exterior derivative (number of rows)
%
% OUTPUT:
%
%   diff_op: sparse matrix representation of the differential operator.
%
% Copyright (C) 2020-2023 Bernard Kapidani, Rafael Vazquez
%
%    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/>.

function diff_op = sp_exterior_derivative (space1, space2)

  assert (strcmpi(space1.transform, 'grad-preserving'), ...
    'The first space cannot be the one for integral-preserving splines (or n-forms)')

  ndim = numel (space1.knots);
  if (ndim == 1)
    grad_curl = 'grad';
    assert (space1.degree == space2.degree+1, 'The degrees are not compatible')
    assert (numel(space1.knots{1}) == numel(space2.knots{1})+2, 'The knot vectors are not compatible')
  elseif (ndim == 2)
    if (strcmpi (space2.transform, 'curl-preserving'))
      grad_curl = 'grad';
      deg_shift = {[1 0], [0 1]};
      knt_shift = {[2 0], [0 2]};
    elseif (strcmpi (space2.transform, 'div-preserving'))
      grad_curl = 'curl';
      deg_shift = {[0 1], [1 0]};
      knt_shift = {[0 2], [2 0]};
    else
      error ('The second space should be either curl-preserving or div-preserving')
    end
    for idim = 1:ndim
      assert (all(space1.degree == space2.scalar_spaces{idim}.degree+deg_shift{idim}), 'The degrees are not compatible')
      assert (all(cellfun(@numel,space1.knots) == (cellfun(@numel, space2.scalar_spaces{idim}.knots)+knt_shift{idim})), ...
        'The knot vectors are not compatible')
    end
  elseif (ndim == 3)
    grad_curl = 'grad';
    deg_shift = {[1 0 0], [0 1 0], [0 0 1]};
    knt_shift = {[2 0 0], [0 2 0], [0 0 2]};
    assert (strcmpi(space2.transform, 'curl-preserving'), ...
      'The second space should be the one for curl-conforming splines')
    for idim = 1:ndim
      assert (all(space1.degree == space2.scalar_spaces{idim}.degree+deg_shift{idim}), 'The degrees are not compatible')
      assert (all(cellfun(@numel,space1.knots) == (cellfun(@numel, space2.scalar_spaces{idim}.knots)+knt_shift{idim})), ...
        'The knot vectors are not compatible')
    end
  end

  diff_ops = op_geom_exterior (space1.knots, space1.degree, grad_curl);
  diff_op = diff_ops{1};
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [vol_tmp,H] = iconebeamdemo2007(dir,flprefix,theta,interp,filter,d,N,step,Ro,range_slices,cut)

%   author: Gianni Schena Sept. 2004 (schena@univ.trieste.it - I welcome suggestions !)
%   Reconstructs image/s for cone beam (CB) geometry according to
%   the classical Feldkamp et al. algorithm. Individual slices are stored in a volume
%
%   FDK10_CB reconstructs the 3D volume image vol from planar projection
%   data. The routine assumes that the center of rotation
%   is the center point of the projections.
%
%   THETA describes the angles (in degrees) at which the projections
%   were taken.  It it is a vector containing angle between
%   projections. THETA is a vector and must contain angles with
%   equal spacing between them.
%
%   The routine uses the filtered backprojection.The filter is designed directly
%   in the frequency domain and then multiplied by the FFT of the
%   projections.  The projections are zero-padded to a power of 2
%   before filtering to prevent spatial domain aliasing and to
%   speed up the FFT.
%
%   INTERP specifies the type of interpolation to use in the
%   backprojection.  The available options are listed in order
%   of increasing accuracy and computational complexity:
%
%      'nearest' - nearest neighbor interpolation
%      'linear'  - linear interpolation (suggested)
%
%   FILTER specifies the filter to use for frequency domain filtering.
%   FILTER is a string that specifies any of the following standard
%   filters:
%
%   'Ram-Lak'     The cropped Ram-Lak or ramp filter (default).  The
%                 frequency response of this filter is |f|.  Because
%                 this filter is sensitive to noise in the projections,
%                 one of the filters listed below may be preferable.
%   'Shepp-Logan' The Shepp-Logan filter multiplies the Ram-Lak
%                 filter by a sinc function.
%   'Cosine'      The cosine filter multiplies the Ram-Lak filter
%                 by a cosine function.
%   'Hamming'     The Hamming filter multiplies the Ram-Lak filter
%                 by a Hamming window.
%   'Hann'        The Hann filter multiplies the Ram-Lak filter by
%                 a Hann window.
%
%   D is a scalar in the range (0,1] that modifies the filter by
%   rescaling its frequency axis.  The default is 1.  If D is less
%   than 1, the filter is compressed to fit into the frequency range
%   [0,D], in normalized frequencies; all frequencies above D are set
%   to 0.
%
%   N is a scalar that specifies the number of rows and columns in the
%   reconstructed 2D slice images.  If N is not specified, the size is determined
%
%   [vol,H] = FDKX... (...) returns also the frequency response of the filter
%   in the vector H.
%
%   Ro is the orbital radius i.e. the distance source to object
%   Ds2dct is the distance between source and planar detector
%   Ds2dct = Ro + Do2dct = Ds2obj + Do2dct
%   pxlsz is the size of the (square) pixel of the detector dct
%   Class Support : All input arguments must be of class double.
%   The output arguments are of class double.

%   References:
%   A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
%   Imaging", IEEE Press 1988. Chapter 3 (available by downloading)

%[p,theta,filter,d,interp,N] = parse_inputs(varargin{:});
thetag=theta;
theta = pi*theta/180; % from degrees to radiants

A=read_1rad(1.,dir,flprefix); % read in image frame
figure , imshow(A,[]); title(' first projection') % disp first image file

%global flat , global dark

% get size of the data files
[nr_p , nc_p] = size(A) ; % get the number of rows and colums of the data files
%
range_slices=range_slices(range_slices>-nr_p/2 & range_slices<nr_p/2);
% remove slices out of range
%  N is a scalar that specifies the number of rows and columns in the
%  reconstructed 2D slice images.  If N is not specified, the size is determined
%  If the user didn't specify the size of the reconstruction, so
%  deduce it from the length of projections

% nr_p; % sinogram size i.e. one size of the projection

if N==0  ||  isempty(N) % if N is not given then ...
    N = 2*floor( nc_p /(2*sqrt(2)) ); % take the default value
end

imgDiag = 2*ceil(N/sqrt(2))+1;  % largest distance through image.

%np = length(theta);
W=pix_weight(nr_p,nc_p,Ro,0,0); % weights to account for a planar detector
%size(W)

% Design the filter H
len=nc_p;
H = designFilter(filter, len, d);
LH=length(H); % length filter in the line vector H

% Define the x & y axes for the reconstructed image so that the origin
xax = (1:N)-ceil(N/2); x = repmat(xax, N, 1);  y= repmat(xax', 1, N);
%y = rot90(x);   % x coordinates, the y coordinates are rot90(x)
mask= ((x-0).^2+(y-0).^2) <= (N/2).^2 ; % mask the inner circle of diam N
mask=true(N); % no mask : all  the pixels are to be reconstructed
maskV=mask(:); % mask in a vector

costheta = cos(theta); sintheta = sin(theta);
ctrIdx = ceil(len/2);     % index of the center of the projections (on the hor axis)

% range slices  to be recons
if exist('range_slices','var') & ~isempty(range_slices)
    slices=range_slices; % slices to be recons
    if isempty( find(slices==0) ) % se non c'e' la slice centrale aggiungila
        slices=[slices ,0];  slices=sort(slices) ;
    end
else
    % if no slices then recns the slice z=0 only
    slices=[0];% [-floor(nc_p/2):step:+floor(nc_p/2)];
end
cnt_slc=find(slices==0);
slices

%vol = repmat(0,[N,N,length(slices)]) ;% allocate 3D volume matrix
vol_tmp=repmat(single(0),[(N.^2),length(slices)]); % allocation
szproj=[nr_p,nc_p;]; % size(proj);
%

% Filtered Back-projection - i.e. from proj pixel (u,v) to voxel (x,y,z)
if strcmp(interp, 'nearest neighbor')
    x=x(maskV); y=y(maskV); % mask !

    nw_slices = real(sort(complex(slices))) ; % sort and reordering for faster reconstruction
   % eg. +5 -5 + 6 -6 ...it allows to change the sign of v rather than recalculating v

    ZFlag= abs(nw_slices(2:end)) == nw_slices(1:(end-1)) ; % logical variable
    ZFlag=real([0, ZFlag] ); % true if +z is followed by -z and one can exploit symmetry

    for ipos=1:length(slices) % right positioning the re-orded slices in the volume
        iposz_v(ipos)=find( nw_slices(ipos)== slices) ;  % slice correspondence
    end
    % data type and array allocation statements
    proj=repmat(0,nr_p,nc_p);
    tmp=zeros(nnz(mask),1); tmp2=zeros(N.^2,1); % allocate memory to temporary arrays
   % tmp1
    v0=zeros([size(x)]); v1=v0; u1=v1; % u=u1; v=u;
    good_u=logical(v0); idx=zeros(size(v0)); % allocation

    for i=1:length(theta) % loop over the projections
       % i
        dgri=thetag(i); % theta in degrees

        proj=read_1rad(i,dir,flprefix); % read in the i-th radiograph
        proj = radio2proj(proj,W,H,len); % tranforms radiog. into projection and multiply by W
        %%%
        proj(:,length(H))=0;  % Zero pad projections for
        for ir = 1:nr_p; % filtering row  by row  with row vector H
            proj(ir,:)= real (ifft ( fft(proj(ir,:)).*H ) );  % fast frequency domain filtering
        end
        proj(:,len+1:end) = [] ;   % Truncate the filtered projections
        %%%
        s =  x.*costheta(i) + y.*sintheta(i); t = -x.*sintheta(i) + y.*costheta(i); % Goddard & Trepanier
        % s = x.*costheta(i) + y.*sintheta(i); % also in matlab iradon
        R_s = Ro-s; u=round( Ro*t./(R_s) + ctrIdx) ; % u is the 1st projection pixel coordinate
        good_u = (u>0 & u<(nc_p+1)) ; % only valid pixels -i.e. within the projection frame
        lengoodu=length(good_u);
        u1=u(good_u);%
        inv_R_s_quad=1./(R_s(good_u).^2); Ro_dev_R_s=Ro./(R_s(good_u)); % for generally .* is faster than ./
        iz=0;
        id_1_term=(u1-1)*szproj(1); % does not depend up on v
        for z=nw_slices % for all the (re-ordered) slices
            iz=iz+1;  % z slice under reconstruction
           if ZFlag(iz)==1; % z symmetry (floting point is faster than logical operators!)
          %  if ( iz>=2 & ( abs(nw_slices(iz)) == nw_slices(iz-1) ) ) % z symmetry
                v0=-v0;   % allows calc time saving
            else %  v is the 2nd proj. pixel coord. (u,v)
                v0=Ro_dev_R_s.*z ;  %
            end
            v1=round( v0+nr_p/2 ) ;  % coord. correction
            %v1=v(good_u);  % logical operations i.e. only valid pixels
            idx = id_1_term + v1 ; % idx = sub2ind(szproj,u1, v1);
            tmp = proj(idx); % just in case
            tmp = tmp(:).*inv_R_s_quad; % back-projection ... from pixels to voxels ...!
            tmp1(good_u)=tmp; tmp1(lengoodu)=0;%
            tmp2(maskV)=tmp1; % prepare for image reconst. within the mask
            iposz=iposz_v(iz); % position of the current slice iz in the volume at iposz
            vol_tmp(:,iposz)=vol_tmp(:,iposz)+tmp2; % sum up back-projections without reshaping
        end % loop over z
    end % end loop over theta for the nearest neighbout method
    vol_tmp=reshape(vol_tmp,[N,N,length(slices)]) ; % reshape 1D vector into 2D image for each z
    % shows central slice only , in figure 10
    figure(10), imshow(fliplr(vol_tmp(:,:,cnt_slc)'),[]); % visualize central slice after vol reconstruction
    stringa = [' NN back-prj. cntrl slice , radiog. degrees: ' , num2str(dgri)];
    title(stringa);   colormap(gray(256));
    display(' done cone beam back projection  ! - nearest neighbour method')
    %***********************************************************
end

if  strcmp(interp, 'bilinear')
    x=x(mask); y=y(mask); % mask !
    temp_1=zeros(nnz(mask(:)),1); temp_2=zeros(nnz(mask(:)),1);
    tmp2=zeros(N.^2,1); % viene retroproiettato
    for i=1:length(theta) % loop over projections
    %   i
        dgri=thetag(i); % angolo corrente in gradi
        proj=read_1rad(i,dir,flprefix); % read in the i-th projection
        proj = radio2proj(proj,W,H,len);

        proj(:,length(H))=0;  % Zero pad projections
        for ir = 1:nr_p; % filtering row  by row
            proj(ir,:)=real (ifft ( fft(proj(ir,:)).*H ) );  % fast frequency domain filtering
        end % done for each row
        proj(:,len+1:end) = [];   % Truncate the filtered projections
        proj=proj(:);
        proj1=[proj(2:end,:); proj(end,:)]; proj1=proj1(:);
        % imshow(proj,[ ]) ;  title(' log proj');
        % proj. is transposed with respect to the original radiograph!
        s =  x.*costheta(i) + y.*sintheta(i); t = -x.*sintheta(i) + y.*costheta(i); % Goddard&Trepanier
        R_s=Ro-s; u=Ro*t./R_s; u=u+ctrIdx; fu=floor(u(:)); %
        iz=0;
        good=( (fu>0 & fu<nc_p+1) );  lengoodu=length(good);
        for z= slices
            iz=iz+1  ;%z % slice under reconstruction
            v=Ro*z./R_s + nr_p/2;  % coord. correction
            fv=floor(v(:));
           % good=( (fu>0 & fu<=nc_p) & (fv>0 & fv<=nr_p) );
            fu1=fu(good); fv1=fv(good);
            idx = (fu1-1)*szproj(1)+fv1 ; %idx = sub2ind(size(proj),fu, fv);
            temp_1(good) = (u(good)-fu1).*proj1(idx) + (fu1+1-u(good)).*proj(idx);
            proj=[proj(:,(2:end)),proj(:,end)];   proj1=[proj(2:end,:);proj(end,:)];
            temp_2(good) = (u(good)-fu1).*proj1(idx) + (fu1+1-u(good)).*proj(idx);
            R_s1=R_s(:);
            tmp =( (v(good)-fv1) .* temp_2(good)  + (fv1+1-v(good)) .* temp_1(good) ) ./ (R_s1(good).^2);
            tmp1(good)=tmp;  tmp1( lengoodu)=0; tmp2(mask(:))=tmp1;
            vol_tmp(:,iz)=vol_tmp(:,iz)+tmp2(:); % $$$$$$$$$$$$$
        end % loop over z
    end % end loop over theta for bilinear method
    vol_tmp=reshape(vol_tmp',[N,N,length(slices)]) ; % reshape $$$$$$$$$$
    % shows central slice in figure 10
    figure(10), imshow(fliplr(vol_tmp(:,:,cnt_slc)'),[]); % visualize central slice after vol reconstruction
    stringa = [' Bi-LNR bk prj. rad. ' , num2str(i-1)];  title(stringa)
    colormap(gray(256));    display(' done back projection  ! - Bilinear method')
end % closee else linear

% Filtered Back-projection - i.e. from proj pixel (u,v) to voxel (x,y,z)
if strcmp(interp, 'n_n_demo')
    x=x(maskV); y=y(maskV); % mask !

    nw_slices = real(sort(complex(slices))) ; % sort and reordering for faster reconstruction
   % eg. +5 -5 + 6 -6 ...it allows to change the sign of v rather than recalculating v

    ZFlag= abs(nw_slices(2:end)) == nw_slices(1:(end-1)) ; % logical variable
    ZFlag=real([0, ZFlag] ); % true if +z is followed by -z and one can exploit symmetry

    for ipos=1:length(slices) % right positioning the re-orded slices in the volume
        iposz_v(ipos)=find( nw_slices(ipos)== slices) ;  % slice correspondence
    end
    % data type and array allocation statements
    proj=repmat(0,nr_p,nc_p);
    tmp=zeros(nnz(mask),1); tmp2=zeros(N.^2,1); % allocate memory to temporary arrays
   % tmp1
    v0=zeros([size(x)]); v1=v0; u1=v1; % u=u1; v=u;
    good_u=logical(v0); idx=zeros(size(v0)); % allocation

    for i=1:length(theta) % loop over the projections
       % i
        dgri=thetag(i); % theta in degrees

        proj=read_1rad(i,dir,flprefix); % read in the i-th radiograph
        proj = radio2proj(proj,W,H,len); % tranforms radiog. into projection and multiply by W
        %%%
        proj(:,length(H))=0;  % Zero pad projections for
        for ir = 1:nr_p; % filtering row  by row  with row vector H
            proj(ir,:)= real (ifft ( fft(proj(ir,:)).*H ) );  % fast frequency domain filtering
        end
        proj(:,len+1:end) = [] ;   % Truncate the filtered projections
        %%%
        s =  x.*costheta(i) + y.*sintheta(i); t = -x.*sintheta(i) + y.*costheta(i); % Goddard & Trepanier
        % s = x.*costheta(i) + y.*sintheta(i); % also in matlab iradon
        R_s = Ro-s; u=round( Ro*t./(R_s) + ctrIdx) ; % u is the 1st projection pixel coordinate
        good_u = (u>0 & u<(nc_p+1)) ; % only valid pixels -i.e. within the projection frame
        lengoodu=length(good_u);
        u1=u(good_u);%
        inv_R_s_quad=1./(R_s(good_u).^2); Ro_dev_R_s=Ro./(R_s(good_u)); % for generally .* is faster than ./
        iz=0;
        id_1_term=(u1-1)*szproj(1); % does not depend up on v
        for z=nw_slices % for all the (re-ordered) slices
            iz=iz+1;  % z slice under reconstruction
           if ZFlag(iz)==1; % z symmetry (floting point is faster than logical operators!)
          %  if ( iz>=2 & ( abs(nw_slices(iz)) == nw_slices(iz-1) ) ) % z symmetry
                v0=-v0;   % allows calc time saving
            else %  v is the 2nd proj. pixel coord. (u,v)
                v0=Ro_dev_R_s.*z ;  %
            end
            v1=round( v0+nr_p/2 ) ;  % coord. correction
            %v1=v(good_u);  % logical operations i.e. only valid pixels
            idx = id_1_term + v1 ; % idx = sub2ind(szproj,u1, v1);
            tmp = proj(idx); % just in case
            tmp = tmp(:).*inv_R_s_quad; % back-projection ... from pixels to voxels ...!
            tmp1(good_u)=tmp; tmp1(lengoodu)=0;%
            tmp2(maskV)=tmp1; % prepare for image reconst. within the mask
            iposz=iposz_v(iz); % position of the current slice iz in the volume at iposz
            vol_tmp(:,iposz)=vol_tmp(:,iposz)+tmp2; % sum up back-projections without reshaping
            figure(8), imshow(reshape(vol_tmp,[N,N,length(slices)]),[ ]) ;
            title('back-projection in progress');
        end % loop over z
    end % end loop over theta for the nearest neighbout method
    vol_tmp=reshape(vol_tmp,[N,N,length(slices)]) ; % reshape 1D vector into 2D image for each z
    % shows central slice only , in figure 10
    figure(10), imshow(fliplr(vol_tmp(:,:,cnt_slc)'),[ ]); % visualize central slice after vol reconstruction
    stringa = [' NN back-prj. cntrl slice , radiog. degrees: ' , num2str(dgri)];
    title(stringa);   colormap(gray(256));
    display(' done cone beam back projection  ! - nearest neighbour method')
    %***********************************************************
end

vol_tmp= vol_tmp*Ro.^2*pi/(2*length(theta));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%%  Sub-Function:  designFilter
%%%
function filt = designFilter(filter, len, d)
% Returns the Fourier Transform of the filter which will be
% used to filter the projections
%
% INPUT ARGS:   filter - either the string specifying the filter
%               len    - the length of the projections
%               d      - the fraction of frequencies below the nyquist
%                        which we want to pass
%
% OUTPUT ARGS:  filt   - the filter to use on the projections

order = max(64,2^nextpow2(2*len));

% First create a ramp filter - go up to the next highest
% power of 2.

filt = 2*( 0:(order/2) )./order;
w = 2*pi*(0:size(filt,2)-1)/order;   % frequency axis up to Nyquist

switch lower(filter)
    case 'ram-lak'
        % Do nothing
    case 'shepp-logan'
        % be careful not to divide by 0:
        filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)));
    case 'cosine'
        filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d));
    case 'hamming'
        filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d));
    case 'hann'
        filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2;
    otherwise
        filter
        error('Invalid filter selected.');
end

filt(w>pi*d) = 0;                      % Crop the frequency response
filt = [filt , filt(end-1:-1:2)];    % Symmetry of the filter
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = pix_weight(nr,nc,R,Nofs,Eofs)
% weight for the pixels of the projetions/radiographs
% R is the orbital radius i.e. source to object distance
% off set ... Nofs=Nord off set ; Eofs= East ..
% nr, nc number of rows and columns of W

if nargin <4 % missing off set values
    Nofs=0;Eofs=0; % no off set then
end

if abs(Nofs) > floor(nc/2) | abs(Eofs) > floor(nr/2)
    display(' off set overflow in pix_weight')
end
% grid !
[X,Y]=meshgrid([[floor(nc/2):-1:0],[1:(nc-floor(nc/2)-1)]],...
    [[floor(nr/2):-1:0],[1:(nr-floor(nr/2)-1)]]);
%W=R.* ( (X.^2 + Y.^2 + R^2).^-0.5 );
W=R* ( ((X-Nofs) .^2 + (Y-Eofs).^2 + R^2).^-0.5 );
%figure, imshow(W,[]); title(' pxls weights')
%figure, surf(W); title(' pixel weighting matrix ')
%colormap(jet)
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function y = embed(x, mask)
% %function y = embed(x, mask)
% %	embed x in nonzero elements of a logical mask
% %   x is a 1D array and mask is 2D
%
% if max(size(mask)) <= 1
%     y = x; % peculiar
% else
%     good = find(mask(:) ~= 0);
%     np = length(good);
%     if size(x, 1) ~= np
%         error 'bad input size' %
%     end
%     xdim = size(x);
%     y = zeros(prod(size(mask)), prod(xdim(2:end)));
%     y(good,:) = x;
%     y = reshape(y, size(mask), xdim(2:end));
% end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function proj = radio2proj(radio,W,H,len)
% the routine corrects the radiog. for flat and dark images
% weights the radiographs
% log process the radio to generate a projection
global flat , global dark
%
% num=(radio-dark); den=(flat-dark);
% num(num<=0)=1.E-6; den(den<=0)=1.E-4;
% proj=num./den; % correzione per dark & flat
% proj(proj>1)=1; proj(proj<=0)=1E-6;
% proj=-log(proj); % da radio 2 proj
%size(radio)
proj=radio;
proj=proj.*W; % weighting projection
%return
% filtring projection
% proj(:,length(H))=0;  % Zero pad projections
% %proj = fft(proj);    % p holds fft of projections
% for ic = 1:size(proj,1)
%     % filter column by column i.e. row sino by row sino
%     %proj(:,ic) = proj(:,ic).*H(:); % fast frequency domain filtering
%     proj(ic,:)=real (ifft ( fft(proj(ic,:)).*H ) );
% end
% %proj = real(ifft(proj));     % p is the filtered projections
% proj(:,len+1:end) = [];   % Truncate the filtered projections
%figure(6),  imshow(proj,[ ]); title(' proiezione filtrata ')
return


text:
function [W,H] = NNDSVD(A, k)
nSample = size(A, 1);
nFeature = size(A, 2);
k = k + 1;
% --------------------------------------------
[U, S, V] = svds(A, k);

% --------------------------------------------
W = zeros(nSample, k);
H = zeros(k, nFeature);
W(:,1) = sqrt(S(1,1))*U(:,1);
H(1,:) = sqrt(S(1,1))*V(:,1)';
for i = 2:k
    x = U(:,i);
    y = V(:,i);
    xp = pos(x); xn = neg(x);
    yp = pos(y); yn = neg(y);
    xpnorm = norm(xp); ypnorm = norm(yp); mp = xpnorm*ypnorm;
    xnnorm = norm(xn); ynnorm = norm(yn); mn = xnnorm*ynnorm;
    if mp > mn
        u = xp/xpnorm;
        v = yp/ypnorm;
        sigma = mp;
    else
        u = xn/xnnorm;
        v = yn/ynnorm;
        sigma = mn;
    end
    W(:,i) = sqrt(S(i,i)*sigma)*u;
    H(i,:) = sqrt(S(i,i)*sigma)*v';
end
W = W(:, 2:k);
H = H(2:k, :);
end

function [xp] = pos(x)
xp = (x >= 0).*x;
end

function [xn] = neg(x)
xn = (x < 0).*(-x);
end


text:
function makeDatFromCTF(filesearch,datfilename,typego,typebutton)

% makeDatFromCTF
% filesearch is a string that if used with dir, can identify all files
% necssary to include. Make sure you're in the right directory. It should
% lead to files with columns of continuous STIM channel output
% typego lists the types of stim that one should expect button responses
% to.
% typebutton is the button response one should get. 0 if not buttonpress

if (nargin < 3)
    typego = [6 7 8 9];
    typebutton = [2 0 0 0];
end

respvals = [33554432 65536 262144 327680 393216 458752 524288 589824];
resptype = [2 1 4 5 6 7 8 9];
maxtypeval = 589824; % above this should be button presses only.
RESPLENGTH = 50;

fs = 1200;
buttonwait = fs*1.5;
numHeaderLines = 19;
files = dir(filesearch);

allStim = [];
for i = 1:length(files)
    fprintf('Getting Data from %s ... \n',files(i).name);
    [time,stim] = textread(files(i).name,'%f %f','delimiter','\t','headerlines',2);
    allStim = [allStim;stim];
    clear time stim;
end

fid = fopen(datfilename,'W');

for i = 1:numHeaderLines
    fprintf(fid,'\n');
end

fprintf(fid,'Trial\tResp\tType\tCorrect\tLatency\tStim/Resp\n');
fprintf(fid,'-----\t----\t----\t-------\t-------\t---------\n');
fprintf(fid,'1\t0\t1\t1\t1000\tStim\n');

i = 1;
trialnum = 2;
while i <= length(allStim)
    nextJump = 1;
    tmp = allStim(i:min(i+RESPLENGTH,length(allStim)));
    if (tmp(1) > 10 && ~isempty(find(respvals == max(tmp))) && max(tmp) <= maxtypeval)
        type = resptype(find(respvals == max(allStim(i:i+RESPLENGTH))));

        correct = 1;
        latency = 0;
        buttonPressed = 0;

        if (~isempty(find(typego == type))) % s2 occured, look for button press
            validButton = typebutton(find(typego == type));
            j = i + RESPLENGTH;
            correct = 0;

            while j < i + buttonwait % see if a button was pressed during the wait period
                tmp = allStim(j:min(j+RESPLENGTH,length(allStim)));
                if (~isempty(tmp) && tmp(1) > 10 && ~isempty(find(respvals == max(tmp))))
                    latency = round((j - i) * 1000 / fs);
                    buttonPressed = resptype(find(respvals == max(tmp)));
                    j = i + buttonwait;

                end
                j = j + 1;
            end

            if (buttonPressed == validButton)
                correct = 1;
            end

        end

        if (~buttonPressed) % get Latency if not button pressed
            j = i + RESPLENGTH;
            cont = 1;
            while (cont && j < length(allStim))
                j = j + 1;
                tmp = allStim(j:min(j+RESPLENGTH,length(allStim)));
                if (isempty(tmp))
                    cont = 0;
                elseif (tmp(1) > 10 && ~isempty(find(respvals == max(tmp))) && max(tmp) <= maxtypeval)
                    cont = 0;
                end
            end
            latency = round((j - i) * 1000 / fs);
            nextJump = (j - i) - 1;
        else
            nextJump = round(latency * fs / 1000) + RESPLENGTH;
        end
        fprintf(fid,'%d\t%d\t%d\t%d\t%d\tStim\n',trialnum,buttonPressed,type,correct,latency);
        trialnum = trialnum+1;
    end
    i = i + nextJump;
end

fclose(fid);

fprintf('***********************************************\n');
fprintf('%s created successfully: %d events found\n\n',datfilename,trialnum);

author: PatternRecognition

repo: OpenBMI


text:
function calpak_test46 ( )

%*****************************************************************************80
%
%% CALPAK_TEST46 tests NOW_TO_JED, NOW_TO_YJF_COMMON, NOW_TO_YMDF_COMMON, NOW_TO_YMDHMS_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST46\n' );
  fprintf ( 1, '  For the current time and date, "NOW", \\n' );
  fprintf ( 1, '  NOW_TO_JED returns the JED,\n' );
  fprintf ( 1, '  NOW_TO_YJF_COMMON the YJF date,\n' );
  fprintf ( 1, '  NOW_TO_YMDF_COMMON returns the YMDF date,\n' );
  fprintf ( 1, '  NOW_TO_YMDHMS_COMMON the YMDHMS date.\n' );
  fprintf ( 1, '\n' );

  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '  "NOW" is %s\n', s );

  jed = now_to_jed ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NOW_TO_JED_COMMON:    Now is: %11.1f\n', jed );

  [ y, j, f ] = now_to_yjf_common ( );
  s = yjf_to_s_common ( y, j, f );

  fprintf ( 1, '  NOW_TO_YJF_COMMON:    Now is: %s\n', s );

  [ y, m, d, f ] = now_to_ymdf_common ( );
  s = ymdf_to_s_common ( y, m, d, f );

  fprintf ( 1, '  NOW_TO_YMDF_COMMON:   Now is: %s\n', s );

  [ y, m, d, h, n, second ] = now_to_ymdhms_common ( );
  s = ymdhms_to_s_common ( y, m, d, h, n, second );

  fprintf ( 1, '  NOW_TO_YMDHMS_COMMON: Now is: %s\n', s );

  return
end


text:
function W = compute_mesh_weight(vertices,faces,type,options)

% compute_mesh_weight - compute a weight matrix
%
%   W = compute_mesh_weight(vertices,faces,type,options);
%
%   W is sparse weight matrix and W(i,j)=0 is vertex i and vertex j are not
%   connected in the mesh.
%
%   todo:
%       1. validate spring weights
%       2. add Mean_curvature weights
%   type is either
%       'combinatorial': W(i,j)=1 is vertex i is conntected to vertex j.
%       'distance': W(i,j) = 1/d_ij^2 where d_ij is distance between vertex
%           i and j.
%       'spring': W(i,j) = 1/d_ij where d_ij is distance between vertex
%           i and j.
%       'conformal' or 'dcp': W(i,j) = cot(alpha_ij)+cot(beta_ij) where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j). Refer to Skeleton
%           Extraction by Mesh Extraction_08, and Intrinsic Parameterizations of Surface Meshes_02.
%       'Laplace-Beltrami': W(i,j) = (cot(alpha_ij)+cot(beta_ij))/2 where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j). Refer to Computing discrete minimal surfaces
%           andtheir conjugates_93, Lemma 2 of On the convergence of metric and geometric properties of
%           polyhedral surfaces_06, and Characterizing Shape Using Conformal Factors_08, and ,
%       'Mean_curvature',W(i,j) = (1/area_i)*(cot(alpha_ij)+cot(beta_ij))/2 where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j), area_i is the
%           area of vertex i's Voroni vicinity. Refer to ??
%       'mvc': W(i,j) = [tan(/_kij/2)+tan(/_jil/2)]/d_ij where /_kij and /_jil
%           are angles at i
%
%   If options.ring is offered, the computation of it can be avoided.
%
%   Add spring and mvc weight (JJCAO, 2009)
%   Add options.ring (JJCAO, 2009)
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
[vertices,faces] = check_face_vertex(vertices,faces);

n = max(max(faces));

if isfield(options, 'verb')
    verb = options.verb;
else
    verb = n>5000;
end

if nargin<3
    type = 'dcp';
end

switch lower(type)
    case 'combinatorial'
        W = triangulation2adjacency(faces);
    case 'distance'
        W = my_euclidean_distance_2(triangulation2adjacency(faces),vertices);
        W(W>0) = 1./W(W>0);
    case 'spring'
        W = sqrt(my_euclidean_distance(triangulation2adjacency(faces),vertices));
        W(W>0) = 1./W(W>0);
    case {'conformal','dcp'} % conformal laplacian
        W = compute_mesh_weight_dcp(vertices, faces, verb);
    case 'laplace-beltrami' %
        W = compute_mesh_weight_dcp(vertices, faces, verb)*0.5;
    case 'mvc'% mvc laplacian
        if isfield(options, 'rings')
            rings = options.rings;
        else
            rings = compute_vertex_face_ring(faces);
        end
        W = sparse(n,n);
        for i = 1:n
            if verb
                progressbar(i,n);
            end
            for b = rings{i}
                % b is a face adjacent to a
                bf = faces(:,b);
                % compute complementary vertices
                if bf(1)==i
                    v = bf(2:3);
                elseif bf(2)==i
                    v = bf([1 3]);
                elseif bf(3)==i
                    v = bf(1:2);
                else
                    error('Problem in face ring.');
                end
                j = v(1); k = v(2);
                vi = vertices(:,i);
                vj = vertices(:,j);
                vk = vertices(:,k);
                % angles
                alpha = myangle(vi-vk,vi-vj);
                % add weight
                W(i,j) = W(i,j) + tan( 0.5*alpha )/sqrt(sum((vi-vj).^2));
                W(i,k) = W(i,k) + tan( 0.5*alpha )/sqrt(sum((vi-vk).^2));
            end
        end
    otherwise
        error('Unknown type.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = compute_mesh_weight_dcp(verts, faces, verb)
n = size(verts,2);
W = sparse(n,n);
% new implementation
for i = 1:size(faces,2)
    c = faces(:,i);

    v1 = verts(:,c(1));
    v2 = verts(:,c(2));
    v3 = verts(:,c(3));

    cot1 = dot(v2-v1,v3-v1)/norm(cross(v2-v1,v3-v1));
    cot2 = dot(v3-v2,v1-v2)/norm(cross(v3-v2,v1-v2));
    cot3 = dot(v1-v3,v2-v3)/norm(cross(v1-v3,v2-v3));

    W(c(2),c(3)) = W(c(2),c(3)) + cot1;
    W(c(3),c(2)) = W(c(3),c(2)) + cot1;
    W(c(3),c(1)) = W(c(3),c(1)) + cot2;
    W(c(1),c(3)) = W(c(1),c(3)) + cot2;
    W(c(1),c(2)) = W(c(1),c(2)) + cot3;
    W(c(2),c(1)) = W(c(2),c(1)) + cot3;
end

%%
% old implementation
% for i = 1:n
%     if verb
%         progressbar(i,n);
%     end
%     for b = rings{i}
%         % b is a face adjacent to a
%         bf = faces(:,b);
%         % compute complementary vertices
%         if bf(1)==i
%             v = bf(2:3);
%         elseif bf(2)==i
%             v = bf([1 3]);
%         elseif bf(3)==i
%             v = bf(1:2);
%         else
%             error('Problem in face ring.');
%         end
%         j = v(1); k = v(2);
%         vi = verts(:,i);
%         vj = verts(:,j);
%         vk = verts(:,k);
%
%         u = vk-vi; v = vk-vj;
%         W(i,j) = W(i,j) + dot(u,v)/norm(cross(u,v));
%         u = vj-vi; v = vj-vk;
%         W(i,k) = W(i,k) + dot(u,v)/norm(cross(u,v));
%     end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function beta = myangle(u,v)
du = sqrt( sum(u.^2) );
dv = sqrt( sum(v.^2) );
du = max(du,eps); dv = max(dv,eps);
beta = acos( sum(u.*v) / (du*dv) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = my_euclidean_distance_2(A,vertex)
% square euclidean distance
if size(vertex,1)<size(vertex,2)
    vertex = vertex';
end

[i,j,s] = find(sparse(A));
d = sum( (vertex(i,:) - vertex(j,:)).^2, 2);
W = sparse(i,j,d);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = my_euclidean_distance(A,vertex)
% euclidean distance
if size(vertex,1)<size(vertex,2)
    vertex = vertex';
end

[i,j,s] = find(sparse(A));
d = sum( (vertex(i,:) - vertex(j,:)).^2, 2).^0.5;
W = sparse(i,j,d);

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function OUT = MovSum(DATA,window)
% =======================================================================
% Moving sum of the vector (or matrix) DATA (T obs x N variables). If
% DATA is a matrix moving sum is computed down each column.
% =======================================================================
% OUT = MovAvg(DATA,window)
% -----------------------------------------------------------------------
% INPUT
%   DATA: T observations x N variables
%   window: window of the moving sum
%------------------------------------------------------------------------
% OUPUT
%   OUT: T observations x N variables matrix (the first windows-1
%       obseravations are NaN)
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

if nargin<2,                error('Not enough input.'),          end
if window<=0,               error('window must be positive.'),   end
if (window~=floor(window)), error('window must be an integer.'), end

if min(size(DATA))==1,
    DATA = DATA(:); % forces DATA to be a column vector
end

[nobs,nvar] = size(DATA);
if window>nobs,
    error('window must not be greater than the length of DATA.')
end

temp=[];
for row=1:(nobs-window+1),
    temp = [temp; sum(DATA(row:(row+window-1),:))];
end

OUT = temp;
OUT = [nan(window-1,nvar); OUT]; % add nans to make conformable to original

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [L, Lratio, df] = L_Ratio(Fet, ClusterSpikes, m)

% [L, Lratio] = L_Ratio(Fet, ClusterSpikes)
%
% L-ratio
% Measure of cluster quality
%
% Inputs:   Fet:           N by D array of feature vectors (N spikes, D dimensional feature space)
%           ClusterSpikes: Index into Fet which lists spikes from the cell whose quality is to be evaluated.
%           m:             squared mahalanobis distances, default is to
%                          calculate them directly

% find # of spikes in this cluster
if nargin < 3
	nSpikes = size(Fet,1);
else
	nSpikes = size(m,1);
end

nClusterSpikes = length(ClusterSpikes);
if nClusterSpikes < size(Fet,2)
    L = nan;
    Lratio = nan;
    df = nan;
    warning('more features and spikes, L-Ratio is NaN...')
    return
end
% mark spikes which are not cluster members
NoiseSpikes = setdiff(1:nSpikes, ClusterSpikes);

%%%%%%%%%%% compute mahalanobis distances %%%%%%%%%%%%%%%%%%%%%
if nargin < 3
	m = mahal(Fet, Fet(ClusterSpikes,:));
end

mCluster = m(ClusterSpikes); % mahal dist of spikes in the cluster
mNoise = m(NoiseSpikes); % mahal dist of all other spikes

df = size(Fet,2);

L = sum(1-chi2cdf(m(NoiseSpikes),df));
Lratio = L/nClusterSpikes;


text:
function root = BesselJ_RootsCyl(starts)
% BesselJ_RootsCyl(starts) finds the roots of the equation J'_1(x) = 0
% where J_1 is the first order Bessel function of the first kind.
%
% starts is the number of starting points in the search for
% roots.  The larger it is, the more roots are returned in the root array.
% The default is 20, which returns 6 roots.  starts = 100 returns 32 roots.
% Generally the number of roots is just under 1/3 the number of starting
% points.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

if(nargin==0)
    starts = 20;
end

f = @(x)(0.5*(besselj(0,x) - besselj(2,x)));

% Get a good list of starting points
y = 0:0.1:starts;
dj1 = f(y);
starts = (find(dj1(1:end-1).*dj1(2:end)<0)+1)*0.1;

for s=1:length(starts)
    root(s) = fzero(f,starts(s));
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function answer = strend (str1, str2)
    case_sensitive = true;
    answer = strendi (str1, str2, case_sensitive);
end

%!test
%! myassert(strend('private', 'haha/private'));
%! myassert(~strend('Private', 'haha/private'));
%! myassert(~strend('private', 'haha/private2'));
%! myassert(~strend('packed:rcond:complexNotSupported', ...
%!     'packed:chol:complexNotSupported'));
%! out = strend('private', {'haha/private', 'haha/wrong', ''});
%! myassert(out, [true, false, false])


text:
clear all, close all, s = 25; randn('seed',s), rand('seed',s)
dev = @(x,z) norm(x(:)-z(:))/max([1,norm(x(:)),norm(z(:))]);

n = 194;
x = 5*rand(n,1); y = sin(x); xs = linspace(0,20,1e4)';

mean = {@meanConst}; hyp.mean = 0.2;
sn = 0.2; hyp.lik = log(sn);
ell = 0.2; sf = 2.3; hyp.cov = [log(ell); log(sf)];
v = 0; cov = {@covPPiso,v};

tic
[nlZ,dnlZ,post] = gp(hyp,[],mean,cov,[],x,y);
t = toc;

tic
  K = sparse(feval(cov{:},hyp.cov,x));
  m = feval(mean{:},hyp.mean,x);
  sn2 = sn*sn;
  L = chol(K/sn2+speye(n))';                                % fast column access
  alpha = L'\(L\((y-m)/sn2));
  nlZs = (y-m)'*alpha/2 + sum(log(diag(L))) + n*log(2*pi*sn2)/2;
  dnlZs = hyp;

  nz = 2*nnz(L)-n; r = zeros(nz,1); c = zeros(nz,1); v = zeros(nz,1);    % alloc
  k = 0; Z = zeros(n,n);
  for i=n:-1:1
    li = L(:,i); dii = li(i); jj = find(li); li = li(jj); nj = numel(jj);
    z = zeros(nj,1);
    for ii=nj:-1:1
      j = jj(ii);
      if nj>1, zij = li(2:nj)'*Z(jj(2:nj),j); else zij = 0; end
      zij = (i==j)/dii^2 - zij/dii;
      Z(i,j) = zij; Z(j,i) = zij;
      z(ii) = zij;
    end
    r(k+1  ) = i; c(k+1  ) = i;        v(k+1  ) = z(1    ); k = k+1;
    idx = 1:nj-1;
    r(k+idx) = i; c(k+idx) = jj(1+idx); v(k+idx) = z(1+idx); k = k+nj-1;
    c(k+idx) = i; r(k+idx) = jj(1+idx); v(k+idx) = z(1+idx); k = k+nj-1;
  end
  Q = sparse(r,c,v/sn^2,n,n);

%   Z = solve_chol(L',eye(n)); Z = Z.*(K>0);
%   [r,c,v] = find(K); for k=1:numel(r), v(k) = Z(r(k),c(k)); end
%   Q = sparse(r,c,v/sn^2,n,n);

  for i = 1:numel(hyp.cov)
    dK = sparse(feval(cov{:},hyp.cov,x,[],i));
    dnlZs.cov(i) = sum(sum(Q.*dK))/2 - (alpha'*dK*alpha)/2;
  end
  dnlZs.lik = sn2*(trace(Q)-alpha'*alpha);
  for i = 1:numel(hyp.mean)
    dnlZs.mean(i) = -feval(mean{:},hyp.mean,x,i)'*alpha;
  end
ts = toc;

dev(post.alpha,alpha)+dev(nlZ,nlZs) + ...
dev(dnlZ.cov,dnlZs.cov)+dev(dnlZ.lik,dnlZs.lik)+dev(dnlZ.mean,dnlZs.mean)
fprintf('times %1.2f, %1.2f\n',t,ts)

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function pdf = hypergeometric_pdf ( x, n, m, l )

%*****************************************************************************80
%
%% HYPERGEOMETRIC_PDF evaluates the Hypergeometric PDF.
%
%  Discussion:
%
%    PDF(X)(N,M,L) = C(M,X) * C(L-M,N-X) / C(L,N).
%
%    PDF(X)(N,M,L) is the probability of drawing X white balls in a
%    single random sample of size N from a population containing
%    M white balls and a total of L balls.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the desired number of white balls.
%    0 <= X <= N, usually, although any value of X can be given.
%
%    Input, integer N, the number of balls selected.
%    0 <= N <= L.
%
%    Input, integer M, the number of white balls in the population.
%    0 <= M <= L.
%
%    Input, integer L, the number of balls to select from.
%    0 <= L.
%
%    Output, real PDF, the probability of exactly K white balls.
%

%
%  Special cases.
%
  if ( x < 0 )

    pdf = 1.0;

  elseif ( n < x )

    pdf = 0.0;

  elseif ( m < x )

    pdf = 0.0;

  elseif ( l < x )

    pdf = 0.0;

  elseif ( n == 0 )

    if ( x == 0 )
      pdf = 1.0;
    else
      pdf = 0.0;
    end

  else

    c1 = binomial_coef_log (   m,   x );
    c2 = binomial_coef_log ( l-m, n-x );
    c3 = binomial_coef_log ( l,   n   );

    pdf_log = c1 + c2 - c3;

    pdf = exp ( pdf_log );

  end

  return
end


text:
function tests = test_spm_dcm_fmri_check
% Unit Tests for spm_dcm_fmri_check
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_dcm_fmri_check.m 6790 2016-04-28 14:30:27Z guillaume $

tests = functiontests(localfunctions);

% -------------------------------------------------------------------------
function test_DCM(testCase)
% Simply tests the function doesn't crash with one DCM

import matlab.unittest.constraints.*

data_path = get_data_path();
dcm_file  = fullfile(data_path,'DCM_s1_m1.mat');

DCM = spm_dcm_fmri_check(dcm_file, true);

testCase.fatalAssertThat(DCM, HasField('diagnostics'));
testCase.verifyThat(DCM.diagnostics, HasElementCount(3));

% -------------------------------------------------------------------------
function test_GCM(testCase)
% Simply tests the function doesn't crash with a GCM array

import matlab.unittest.constraints.*

data_path = get_data_path();
gcm_file  = fullfile(data_path,'GCM_simulated.mat');

GCM = spm_dcm_fmri_check(gcm_file, true);

testCase.verifyThat(GCM, IsOfClass('cell'));
testCase.fatalAssertThat(GCM{1}, HasField('diagnostics'));
testCase.verifyThat(GCM{1}.diagnostics, HasElementCount(3));

% -------------------------------------------------------------------------
function data_path = get_data_path()
data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region', 'models');


text:
function M = make_convn_mat(F, sz, shape, pad)
%MAKE_CONVN_MAT   Flexible N-D convolution matrix
%   M = MAKE_CONVN_MAT(F, SZ[, SHAPE[, PAD]]) returns the convolution
%   matrix for the matrix F.  SZ gives the size of the array that the
%   convolution should be applied to.  The returned matrix M is sparse.
%   If X is of size SZ and SHAPE is 'full', then reshape(M * X(:), SZ +
%   size(F) - 1) is the same as convn(X, F).
%
%   The optional parameter SHAPE controls which parts of the convolution
%   result to return in M:
%     * 'full':      (default) returns the full N-D convolution, i.e. the
%                    behavior is identical to convn(X, F, 'full').
%     * 'same':      returns the central part of the convolution that
%                    is the same size as X, i.e. the behavior is
%                    identical to convn(X, F, 'same').
%     * 'sameswap':  identical to 'same' except that rounding needed for
%                    even-sized filters is performed the opposite way.
%     * 'valid':     returns only the part of the result that can be
%                    computed without assuming zero-padded arrays,
%                    i.e. the behavior is identical to convn(X, F,
%                    'valid').
%
%   The optional parameter PAD controls whether to return an M that pads
%   the result of the convolution with zeros (default no padding):
%     * 'full':      returns the result padded to the output size of a
%                    full N-D convolution (works with SHAPE set to 'same'
%                    and 'valid').
%     * 'same':      returns the result padded to the central part of the
%                    convolution that is the same size as X (works with
%                    SHAPE set to 'valid').
%     * 'sameswap':  identical to 'same' except that rounding needed for
%                    even-sized filters is performed the opposite way.
%
%   See also CONVMTXN.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date: 2007-03-27 14:09:11 -0400 (Tue, 27 Mar 2007) $
%   $Revision: 252 $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Default to full matrix
  if (nargin < 3)
    shape = 'full';
  end

  ndims = length(sz);

  % Border sizes for 'same' and 'sameswap'
  Fsize_lo_2 = ceil((size(F) - 1) / 2);
  Fsize_hi_2 = floor((size(F) - 1) / 2);

  % Border sizes for 'valid'
  Fsize = size(F) - 1;

  switch(shape)
    case 'same'
      % Mark valid and invalid pixels (i.e. the ones within and outside
      % of the part to be returned)
      valid = true(sz+size(F)-1);

      for d = 1:ndims
        for e = 1:ndims
          sub{e} = ':';
        end

        sub{d} = 1:Fsize_lo_2(d);
        valid(sub{:}) = false;
        sub{d} = size(valid, d)-Fsize_hi_2(d)+1:size(valid, d);
        valid(sub{:}) = false;
      end

      if (nargin > 3 && strcmp(pad, 'full'))
        % If we're padding to 'full' set the coefficients on the border
        % to zero
        M = convmtxn(F, sz, valid);
      else
        % If we're *not* padding, then suppress the rows of M that
        % correspond to the border
        M = convmtxn(F, sz);
        M = M(valid, :);
      end

   case 'sameswap'
     % Mark valid and invalid pixels (i.e. the ones within and outside
     % of the part to be returned), but round the other way
     valid = true(sz+size(F)-1);

     for d = 1:ndims
       for e = 1:ndims
         sub{e} = ':';
       end

       sub{d} = 1:Fsize_hi_2(d);
       valid(sub{:}) = false;
       sub{d} = size(valid, d)-Fsize_lo_2(d)+1:size(valid, d);
       valid(sub{:}) = false;
     end

     if (nargin > 3 && strcmp(pad, 'full'))
       % If we're padding to 'full' set the coefficients on the border
       % to zero
       M = convmtxn(F, sz, valid);
     else
       % If we're *not* padding, then suppress the rows of M that
       % correspond to the border
       M = convmtxn(F, sz);
       M = M(valid, :);
     end

   case 'valid'
     % Mark valid and invalid pixels (i.e. the ones within and outside
     % of the part to be returned)
     valid = true(sz+size(F)-1);

     for d = 1:ndims
       for e = 1:ndims
         sub{e} = ':';
       end

       sub{d} = 1:Fsize(d);
       valid(sub{:}) = false;
       sub{d} = size(valid, d)-Fsize(d)+1:size(valid, d);
       valid(sub{:}) = false;
     end

     if (nargin > 3)
       % If we're padding, then figure out the area to be padded

       switch (pad)
         case 'same'
           % Mark valid and invalid pixels (i.e. the ones within and outside
           % of the part to be padded)
           pad_valid = true(sz+size(F)-1);

           for d = 1:ndims
             for e = 1:ndims
               sub{e} = ':';
             end

             sub{d} = 1:Fsize_lo_2(d);
             pad_valid(sub{:}) = false;
             sub{d} = size(valid, d)-Fsize_hi_2(d)+1:size(valid, d);
             pad_valid(sub{:}) = false;
           end

           % Set coefficients on the border to zero
           M = convmtxn(F, sz, valid);

           % Suppress rows of M outside of the padded area
           M = M(pad_valid, :);

         case 'sameswap'
           % Mark valid and invalid pixels (i.e. the ones within and outside
           % of the part to be padded), but round the other way
           pad_valid = true(sz+size(F)-1);

           for d = 1:ndims
             for e = 1:ndims
               sub{e} = ':';
             end

             sub{d} = 1:Fsize_hi_2(d);
             pad_valid(sub{:}) = false;
             sub{d} = size(valid, d)-Fsize_lo_2(d)+1:size(valid, d);
             pad_valid(sub{:}) = false;
           end

           % Set coefficients on the border to zero
           M = convmtxn(F, sz, valid);

           % Suppress rows of M outside of the padded area
           M = M(pad_valid, :);

         otherwise
           % Padding to 'full'; only set coefficients on the border to zero
           M = convmtxn(F, sz, valid);
       end
     else
       % No padding; suppress all rows on the border
       M = convmtxn(F, sz);
       M = M(valid, :);
     end

    otherwise
      % Full convolution; return everything
      M = convmtxn(F, sz);

  end

author: kristinbranson

repo: JAABA


text:
function [valuesSelected,costVal]=knapsack01DP(v,w,W)
%%KNAPSACK01DP Solve the 0-1 knapsack problem using dynamic programming.
%        The 0-1 knapsack problem is to find a subset sel of the items in v
%        that solves the optimization problem
%        max_{sel} sum(v(sel))
%        such that sum(w(sel))<=W
%        This can be thought of as fitting the most things in a knapsack if
%        the v's are all 1's, the w's are the sizes of things and W is the
%        capacity of the knapsack.
%
%INPUTS: v An nX1 or 1Xn vector of positive real values to maximize.
%        w An nX1 or 1Xn vector of positive integer capacities required for
%          each value.
%        W The positive integer maximum allowable capacity that can be
%          filled.
%
%OUTPUTS: valuesSelected The indices of the items in v (and w) that are
%                  selected. If the problem is infeasible, then an empty
%                  matrix is returned.
%          costVal The sum of the selected items in v. 0 is returned if the
%                  problem is infeasible.
%
%The knapsack problem is NP-hard. The complexity of this algorithm
%scaled with W. It is O(n*W) The algorithm is sovled via dynamic
%programming. A development of the recursion in such a dynamic programming
%approach is given in Chapter 5.4 of [1]. However, this is a very common
%method for solving this problem and thus solutions can be found in a
%number of basic textbooks on algorithms.
%
%EXAMPLE:
% W=7;
% v=[2.1;3.125;66;4.3;6.2];
% w=[1;2;3;4;5];
% [valuesSelected,costVal]=knapsack01DP(v,w,W)
%One will get vales 3, 2, and 1 selected.
%
%REFERENCES:
%[1] L. A. Wolsey, Integer Programming. New York: Wiley-Interscience, 1998.
%
%January 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(v);

V=zeros(n+1,W+1);
V(0+1,:)=0;
keep=zeros(n,W+1);

for i=1:n
    for wIdx=0:W
        if(wIdx-w(i)>=0)
            VTest=V(i-1+1,wIdx-w(i)+1);
        else
            VTest=-Inf;
        end

        if((w(i)<=wIdx)&&(v(i)+VTest)>V(i-1+1,wIdx+1))
            V(i+1,wIdx+1)=v(i)+VTest;
            keep(i,wIdx+1)=1;
        else
            V(i+1,wIdx+1)=V(i-1+1,wIdx+1);
            keep(i,wIdx+1)=0;
        end
    end
end

K=W;

valuesSelected=zeros(n,1);
numSelected=0;

for i=n:-1:1
    if(keep(i,K+1)==1)
        numSelected=numSelected+1;
        valuesSelected(numSelected)=i;
        K=K-w(i);
    end
end

valuesSelected=valuesSelected(1:numSelected);

costVal=V(n+1,W+1);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary