text:
function value = r8_gamr ( x )

%*****************************************************************************80
%
%% R8_GAMR evaluates the reciprocal gamma function of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the value of the reciprocal gamma
%    function at X.
%
  if ( x <= 0.0 && r8_aint ( x ) == x )

    value = 0.0;

  elseif ( abs ( x ) <= 10.0 )

    value = 1.0 / r8_gamma ( x );

  else

    [ alngx, sgngx ] = r8_lgams ( x );
    value = sgngx * exp ( - alngx );

  end

  return
end


text:
%filename='test2d_micro';
filename = 'MicroQuad';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'chomog_alphabeta'};
weights=[1];
constraint = {'volumeConstraint'};
constraint_case = 'EQUALITY';
%incrementFactor = 1;
designVariable = 'Density';
%designVariable = 'LevelSet';
filterType = 'PDE';
fracRadius = 0.51;
%optimizer = 'IPOPT';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

optimizerUnconstrained = 'PROJECTED GRADIENT';
line_search_initiator = 'INCREASING LAST STEP';
incrementFactor = 1.95;

%optimizerUnconstrained = 'SLERP';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target=1;
optimality_final = 0.2*1e-3;
constr_final =1e-12;

Vfrac_initial = 0.8;
optimality_initial = 0.2*1e-3;
constr_initial = 1e-12;

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;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta =[1 1 0]';

%nsteps = 10;

% For all tests
plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 200;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function loadAllDirs()
	% Forces all directories to be loaded, including external programs normally not loaded.
	% Biafra Ahanonu
	% started: 2020.05.09 [18:39:52]
	% inputs
		%
	% outputs
		%

	% changelog
		%
	% TODO
		%

	%========================
	% DESCRIPTION
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		ciapkg.loadBatchFxns('loadEverything');
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function geometry_test060 ( )

%*****************************************************************************80
%
%% TEST060 tests PLANE_IMP_TRIANGLE_INT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST060\n' );
  fprintf ( 1, '  PLANE_IMP_TRIANGLE_INT_3D finds the\n' );
  fprintf ( 1, '    intersection points of an implicit plane\n' );
  fprintf ( 1, '    and a triangle.\n' );

  a = 1.0;
  b = -2.0;
  c = -3.0;
  d = 6.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The implicit plane: A*X + B*Y + C*Z + D = 0.\n' );
  fprintf ( 1, '  A,B,C,D = %f  %f  %f  %f\n', a, b, c, d );

  for i = 1 : 4

    if ( i == 1 )
      t(1:dim_num,1:3) = [ ...
         3.0,  0.0, -7.0; ...
        13.0, -4.0, -1.0; ...
         5.0,  1.0, -2.0 ]';
    elseif ( i == 2 )
      t(1:dim_num,1:3) = [ ...
         3.0,  0.0, -7.0; ...
        13.0, -4.0, -1.0; ...
         9.0,  3.0,  8.0 ]';
    elseif ( i == 3 )
      t(1:dim_num,1:3) = [ ...
        -6.0,  0.0,  0.0; ...
         0.0,  3.0,  0.0; ...
         0.0,  0.0,  2.0 ]';
    elseif ( i == 4 )
      t(1:dim_num,1:3) = [ ...
        -4.0,  1.0,  0.0; ...
         0.0,  6.0, -2.0; ...
         0.0,  0.0,  1.0 ]';
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Case %d\n', i );

    r8mat_transpose_print ( dim_num, 3, t, '  Triangle vertices:' );

    [ num_int, pint ] = plane_imp_triangle_int_3d ( a, b, c, d, t );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of intersection points is %d\n', num_int );

    r8mat_transpose_print ( dim_num, num_int, pint, '  Intersection points:' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [aps, recall, precision] = calcAP(track_confs, tp_cell, fp_cell, num_vids, num_track_per_class, num_track_thr, defaultTrackThr)

fprintf('computing AP\n');
recall = cell(1,num_track_thr);
precision = cell(1,num_track_thr);
aps = cell(1,num_track_thr);
confs = [track_confs{:}];
[~, ind] = sort(confs,'descend');
for o = 1:num_track_thr
    tp_all = [];
    fp_all = [];
    for v = 1:num_vids
        tp_all = [tp_all(:); tp_cell{v}{o}'];
        fp_all = [fp_all(:); fp_cell{v}{o}'];
    end

    tp_all = tp_all(ind)';
    fp_all = fp_all(ind)';

    % compute precision/recall
    tp = cumsum(tp_all);
    fp = cumsum(fp_all);
    recall{o} = (tp/num_track_per_class)';
    precision{o} = (tp./(fp+tp))';
    aps{o} = VOCap(recall{o},precision{o})*100;
end

fprintf('-------------\n');
ap = aps{1};
for t = 2:length(aps)
    ap = ap + aps{t};
end
ap = ap ./ length(aps);
fprintf('Mean AP:\t\t %0.2f%%\n',mean(ap));
fprintf(' = = = = = = = = \n');
for t = 1:length(aps)
    ap = aps{t};
    fprintf('Mean AP@%0.2f:\t %0.2f%%\n',defaultTrackThr(t),mean(ap));
end
fprintf(' = = = = = = = = \n');

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function theta = polygonNormalAngle(points, ind)
%POLYGONNORMALANGLE Compute the normal angle at a vertex of the polygon
%
%   THETA = polygonNormalAngle(POLYGON, IND);
%   where POLYGON is a set of points, and IND is index of a point in
%   polygon. The function compute the angle of the normal cone localized at
%   this vertex.
%   If IND is a vector of indices, normal angle is computed for each vertex
%   specified by IND.
%
%   Example
%   % creates a simple polygon
%   poly = [0 0;0 1;-1 1;0 -1;1 0];
%   % compute normal angle at each vertex
%   theta = polygonNormalAngle(poly, 1:size(poly, 1));
%   % sum of all normal angle of a non-intersecting polygon equals 2xpi
%   % (can be -2xpi if polygon is oriented clockwise)
%   sum(theta)
%
%   See also:
%   polygons2d, polygonOuterNormal, normalizeAngle
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2005-11-30
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% number of points
np = size(points, 1);

% number of angles to compute
nv = length(ind);

theta = zeros(nv, 1);

for i = 1:nv
    % current vertex
    curr = points(ind(i), :);

    % previous and next vertices
    prev = points(mod(ind(i)-2, np)+1, :);
    next = points(mod(ind(i), np)+1, :);

    theta(i) = angle3Points(prev, curr, next) - pi;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function value = r4_round_i4 ( x )

%*****************************************************************************80
%
%% R4_ROUND_I4 rounds an R4 to the nearest integral value, returning an I4.
%
%  Discussion:
%
%    In MATLAB, it is essentially true that there is little difference between
%    this function and R4_ROUND, because we store our integers in what amounts
%    to a real variable.
%
%  Example:
%
%        X        R4_ROUND_I4
%
%      1.3         1
%      1.4         1
%      1.5         1 or 2
%      1.6         2
%      0.0         0
%     -0.7        -1
%     -1.1        -1
%     -1.6        -2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value.
%
%    Output, integer R4_ROUND_I4, the rounded value.
%
  if ( x < 0.0 )
    value = - floor ( - x + 0.5 );
  else
    value =   floor ( + x + 0.5 );
  end

  return
end


text:
function hRD = plugin_rd_rxtx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.pluto.plugin_rd('Rx & Tx');
AnalogDevices.pluto.add_io(hRD,'Rx & Tx');

author: analogdevicesinc


text:
function y = vl_nnglobalpool(x, varargin)
%VL_NNGLOBALPOOL CNN global poolinng.
%   Y = VL_NNGLOBALPOOL(X) applies the pooling operator to all
%   spatial locations of the data X. X is a SINGLE array of dimension
%   H x W x C x N where (H,W) are the height and width of the map stack,
%   C is the number of feature channels and N the number of of images
%   in the batch.
%
%   DZDX = VL_NNGLOBALPOOL(X, POOL, DZDY) computes the derivatives of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.
%
%   VL_NNGLOBALPOOL(..., 'option', value, ...) takes the following option:
%
%   `method`:: 'avg'
%     Specify method of pooling. It can be either 'max' (retain max value
%     over all spatial locations per channel) or 'avg' (compute the average
%     value over all spatial locations per channel).
%
%   The output a is a SINGLE array of dimensions 1 x 1 x C x N.
%
%   The derivative DZDY has the same dimension of the output Y and
%   The derivative DZDX has the same dimension as the input X.
%
% Copyright (C) 2016 Samuel Albanie and Andrea Vedaldi
% Licensed under The MIT License [see LICENSE.md for details]

  opts.method = 'avg' ;
  [opts, dzdy] = vl_argparsepos(opts, varargin) ;

  if nargin <= 1 || isempty(dzdy)
    switch opts.method
      case 'avg', y = mean(mean(x, 1), 2) ;
      case 'max', y = max(max(x, [], 1), [], 2) ;
      otherwise, error('Pooling method %s not recognized', opts.method) ;
    end
  else
    base = 1 / (size(x,1) * size(x,2)) * ones(size(x), 'like', x) ;
    y = bsxfun(@times, base, dzdy{1}) ;
  end

repo: MFT


text:
function [varargout]=patchPointDist(V,F,Vp,dL)

methodOpt='near-norm';

switch methodOpt
    case 'near-norm'

        [Nn,Vn]=patchNormal(F,V); %Get face normals and normal vector origins
        [Dn,indMin]=minDist(Vp,Vn); %Get distances to origins and find closest faces

        %Attempt to find orthogonal projection to closest face
        N=Nn(indMin,:); %Order the normals
        W=Vp-Vn(indMin,:);
        Dm=(dot(N,W,2));
        D=abs(Dm);
        Dp=Dm(:,ones(1,3)).*N;
        Vpd=Vp-Dp;

        %Test if projection is valid
        TR = triangulation(F,V);
        B = cartesianToBarycentric(TR,indMin,Vpd);
        logicOutside=any(B<-dL,2);

        %Replace invalid projections with face centre points
        D(logicOutside)=Dn(logicOutside); %Replace distances
        Dp(logicOutside,:)=Vp(logicOutside,:)-Vn(indMin(logicOutside),:); %Replace difference vectors

    otherwise
        error('Wrong method option chosen');
end

switch nargout
    case 1
        varargout{1}=D;
    case 2
        varargout{1}=D;
        varargout{2}=Dp;
    case 3
        varargout{1}=D;
        varargout{2}=Dp;
        varargout{3}=logicOutside;
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [mmHg] = psi2mmHg(psi)
% Convert units of pressure from pounds per square inch to millimeters of mercury.
% Chad Greene 2012
mmHg = psi*51.7149;

author: Sable

repo: mcbench-benchmarks


text:
% [HEIGHT] = lpyrHt(INDICES)
%
% Compute height of Laplacian pyramid with given its INDICES matrix.
% See buildLpyr.m

% Eero Simoncelli, 6/96.

function [ht] =  lpyrHt(pind)

% Don't count lowpass residual band
ht = size(pind,1)-1;

repo: SelfExSR


text:
function [ c_vec1, c_vec2 ] = real_cross( p_vec1, p_vec2, pcross_real,...
                                            eta_c, min_realvar, max_realvar)
%   Applies SBX over two vectors of double.

nreal = length(p_vec1);

epsilon = 1.0e-14 ;
% ncross = 0 ;

c_vec1 = zeros(1,nreal);
c_vec2 = zeros(1,nreal);

if(rand(1) <= pcross_real)
% if(randomperc() <= pcross_real) % SLOW !!!
    % ncross = ncross + 1 ;
    if(nreal < 10)
        [c_vec1, c_vec2] = sbx_looped(p_vec1, p_vec2, ...
                                    c_vec1, c_vec2, ...
                                    epsilon, eta_c, ...
                                    min_realvar, max_realvar);
    else
        [c_vec1, c_vec2] = sbx_vectorized(p_vec1, p_vec2, ...
                                    c_vec1, c_vec2, ...
                                    epsilon, eta_c, ...
                                    min_realvar, max_realvar);
    end
else
    c_vec1 = p_vec1 ;
    c_vec2 = p_vec2 ;
end
end

function [c_vec1, c_vec2] = sbx_looped(p_vec1, p_vec2, ...
                                            c_vec1, c_vec2, ...
                                            epsilon, eta_c, ...
                                            min_realvar, max_realvar)
% This is the original implementation
nreal = length(p_vec1);
for i = 1:nreal
   if (rand(1) <= 0.5)
   % if (randomperc() <= 0.5) % SLOW !!!
       if (abs(p_vec1(i) - p_vec2(i)) > epsilon)
            if (p_vec1(i) < p_vec2(i))
                y1 = p_vec1(i) ;
                y2 = p_vec2(i) ;
            else
                y1 = p_vec2(i) ;
                y2 = p_vec1(i) ;
            end
            yl = min_realvar(i);
            yu = max_realvar(i);
            r = rand(1) ;
            % r = randomperc() ; % SLOW !!!
            beta_ = 1.0 + (2.0 * (y1 - yl) / (y2 - y1));
            alpha_ = 2.0 - (beta_ ^ (-1.0 * (eta_c + 1.0)));
            if (r <= (1.0 / alpha_))
                    betaq = (r * alpha_) ^ (1.0/(eta_c + 1.0));
            else
                    betaq = (1.0 / (2.0 - r * alpha_)) ^ (1.0 / (eta_c + 1.0));
            end
            c1 = 0.5 * ((y1 + y2) - (betaq * (y2 - y1)));
            beta_ = 1.0 + (2.0 * (yu - y2)/(y2 - y1));
            alpha_ = 2.0 - (beta_ ^ (-1.0 * (eta_c + 1.0)));
            if (r <= (1.0 / alpha_))
                    betaq = (r * alpha_) ^ (1.0 / (eta_c + 1.0));
            else
                    betaq = (1.0 / (2.0 - r * alpha_)) ^ (1.0 / (eta_c + 1.0));
            end
            c2 = 0.5 * ((y1 + y2) + betaq * (y2 - y1));
            if (c1 < yl)
                c1 = yl;
            end
            if (c2 < yl)
                c2 = yl;
            end
            if (c1 > yu)
                c1 = yu;
            end
            if (c2 > yu)
                c2 = yu;
            end
            if (rand(1) <= 0.5)
            % if (randomperc() <= 0.5) % SLOW !!!
                c_vec1(i) = c2;
                c_vec2(i) = c1;
            else
                c_vec1(i) = c1;
                c_vec2(i) = c2;
            end
        else
            c_vec1(i) = p_vec1(i) ;
            c_vec2(i) = p_vec2(i) ;
        end
    else
        c_vec1(i) = p_vec1(i) ;
        c_vec2(i) = p_vec2(i) ;
    end
end
%
end

function [c_vec1, c_vec2] = sbx_vectorized(p_vec1, p_vec2, ...
                                            c_vec1, c_vec2, ...
                                            epsilon, eta_c, ...
                                            min_realvar, max_realvar)
% This is the vectorized version of the above code, this will generally
% give you 2-times speed up than the above code, especially you will
% observe even more speed up when the length of the variable gets larger.
% However, short variable length may make it slower (like in Osyczka's
% problems)
nreal = length(p_vec1);
randv1lthalf = rand(1,nreal) <  0.5 ;
% randv1lthalf = randompercv(1,nreal) <  0.5 ; % SLOW !!!
absdiffgteps = abs(p_vec1 - p_vec2) > (zeros(1,nreal) * epsilon) ;
xover_index = randv1lthalf & absdiffgteps ;
abs_xover_index = (1:length(xover_index));
abs_xover_index = abs_xover_index(xover_index);
p1 = p_vec1(xover_index);
p2 = p_vec2(xover_index);
len = length(p1);
if(len > 0)
    eta_cv = ones(1, len) * eta_c ;
    rv = rand(1,len);
    % rv = randompercv(1,len); % SLOW !!!
    randv2lthalf = rand(1,len) < 0.5 ;
    % randv2lthalf = randompercv(1,len) < 0.5 ; % SLOW !!!
    %
    pveclt = p1 < p2 ;
    y1v = p1 .* pveclt + p2 .* (~pveclt);
    y2v = p1 .* (~pveclt) + p2 .* pveclt;
    ylv = min_realvar(xover_index).' ;
    yuv = max_realvar(xover_index).' ;
    %
    beta_v = 1.0 + (2.0 .* (y1v - ylv) ./ (y2v - y1v));
    alpha_v = 2.0 - (beta_v .^ (-1.0 .* (eta_cv + 1.0)));
    rltalpha = rv <= (1.0 ./ alpha_v);
    % bsxfun here ?
    betaqv = ((rv .* alpha_v) .^ (1.0 ./ (eta_cv + 1.0))) .* ...
                rltalpha + ((1.0 ./ (2.0 - rv .* alpha_v)) .^ ...
                    (1.0 ./ (eta_cv + 1.0))) .* (~rltalpha);
    c1v = 0.5 .* ((y1v + y2v) - (betaqv .* (y2v - y1v)));
    %
    beta_v = 1.0 + (2.0 .* (yuv - y2v) ./ (y2v - y1v));
    alpha_v = 2.0 - (beta_v .^ (-1.0 .* (eta_cv + 1.0)));
    rltalpha = rv <= (1.0 ./ alpha_v);
    % bsxfun here ?
    betaqv = ((rv .* alpha_v) .^ (1.0 ./ (eta_cv + 1.0))) .* ...
                rltalpha + ((1.0 ./ (2.0 - rv .* alpha_v)) .^ ...
                    (1.0 ./ (eta_cv + 1.0))) .* (~rltalpha);
    c2v = 0.5 .* ((y1v + y2v) + betaqv .* (y2v - y1v));
    %
    c1ltyl = c1v < ylv ;
    c1gtyu = c1v > yuv ;
    c1v = (c1ltyl .* ylv) + (c1v .* (~c1ltyl));
    c1v = (c1gtyu .* yuv) + (c1v .* (~c1gtyu));
    c2ltyl = c2v < ylv ;
    c2gtyu = c2v > yuv ;
    c2v = (c2ltyl .* ylv) + (c2v .* (~c2ltyl));
    c2v = (c2gtyu .* yuv) + (c2v .* (~c2gtyu));
    %
    p1 = c2v .* randv2lthalf + c1v .* (~randv2lthalf);
    p2 = c1v .* randv2lthalf + c2v .* (~randv2lthalf);
end
%
c_vec1 = c_vec1' ;
c_vec1(abs_xover_index.',1) = p1' ;
c_vec1((~xover_index).',1) = p_vec1(~xover_index).' ;
c_vec1 = c_vec1';
%
c_vec2 = c_vec2' ;
c_vec2(abs_xover_index.',1) = p2' ;
c_vec2((~xover_index).',1) = p_vec2(~xover_index).' ;
c_vec2 = c_vec2' ;
%
end

author: chudur-budur


text:
function nnchecknumgrad(nn, x, y)
    epsilon = 1e-6;
    er = 1e-7;
    n = nn.n;
    for l = 1 : (n - 1)
        for i = 1 : size(nn.W{l}, 1)
            for j = 1 : size(nn.W{l}, 2)
                nn_m = nn; nn_p = nn;
                nn_m.W{l}(i, j) = nn.W{l}(i, j) - epsilon;
                nn_p.W{l}(i, j) = nn.W{l}(i, j) + epsilon;
				rand('state',0)
                nn_m = nnff(nn_m, x, y);
                rand('state',0)
                nn_p = nnff(nn_p, x, y);
                dW = (nn_p.L - nn_m.L) / (2 * epsilon);
                e = abs(dW - nn.dW{l}(i, j));

                assert(e < er, 'numerical gradient checking failed');
            end
        end
    end
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
classdef TargetParameters < handle

    properties (Access = public)
        constr_tol
        optimality_tol
        Vfrac
        epsilon
        epsilon_velocity
        epsilon_perimeter
        epsilon_isotropy
        stressNormExponent
        iStep
    end

end


text:
function [imgc, x1, x2, y1, y2] = removeFrame(img, framecolor)
%
% Removes a frame from a picture
%
% [imgc, xmin, xmax, ymin, ymax] = removeFrame(img, framecolor);

img = double(img);
[nrows ncols cc] = size(img);

if nargin == 1
    cf = squeeze(median(img(:,5,:),1)); % color frame (a vertical column of uniform color)
else
    cf = framecolor;
end

frame = (abs(img(:,:,1)-cf(1))+abs(img(:,:,2)-cf(2))+abs(img(:,:,3)-cf(3)))<20; % find frame pixels

bx = sum(frame,1)>(nrows*.9);
by = sum(frame,2)>(ncols*.9);

if sum(bx)+sum(by)>5
    mx = round(ncols/2);
    my = round(nrows/2);

    x1 = 2+max([-1 max(find(bx(1:mx)))]);
    x2 = -2+min([ncols+2 mx+min(find(bx(mx:ncols)))]);
    y1 = 2+max([-1 max(find(by(1:my)))]);
    y2 = -2+min([nrows+2 my+min(find(by(my:nrows)))]);

    imgc = uint8(img(y1:y2, x1:x2, :));
else
    x1 = 1;
    y1 = 1;
    x2 = ncols;
    y2 = nrows;
    imgc = img;
end

author: CSAILVision

repo: LabelMeToolbox


text:
function qoi_quad ( )

%*****************************************************************************80
%
%% QOI_QUAD seeks the expected value of a quantity of interest.
%
%  Discussion:
%
%    Our base value for DELTA_BASE is 0.01.
%
%    Our quantity of interest Q is the time at which the solution achieves
%    the value 0.99.
%
%    Our parameter U is uniformly distributed in [-1,+1], and determines
%    our actual value of DELTA by
%
%      DELTA = 2^U * DELTA_BASE.
%
%    We seek the expected value of Q, which we define as:
%
%      E(Q) = Integral ( -1 <= U <= +1 ) Q(Y(DELTA(U))) rho(U) dU
%
%    where
%
%      rho(U) = 1/2;
%      DELTA(U) = 2^U * DELTA_BASE;
%      Y(DELTA()) is the solution of the combustion ODE for an initial
%        condition of DELTA;
%      Q(Y) is the (first) time T* at which the solution Y(T*) = 0.99.
%
%    We approximate the integral using a Clenshaw-Curtis quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
  delta_base = 0.01;
  n_max = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QOI_QUAD\n' );
  fprintf ( 1, '  Use quadrature to estimate the expected value of our\n' );
  fprintf ( 1, '  quantity of interest Q, the time at which the\n' );
  fprintf ( 1, '  combustion solution reaches 0.99.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using Clenshaw-Curtis quadrature of orders 1 through %d\n', n_max );
%
%  Set up an option to have ODE45 return when the solution reached 0.99.
%
  options = odeset ( 'EVENTS', @event_function );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N      Estimated Q\n' );
  fprintf ( 1, '\n' );

  for n = 1 : n_max

    [ x, w ] = clenshaw_curtis_compute ( n );

    f = zeros ( n, 1 );

    for j = 1 : n
%
%  Uniformly select
%    U, the logarithm base 2, between -1 and +1, of
%    F, a factor between 1/2 and 2, which multiplies
%    DELTA_BASE, the base value, which gives us
%    DELTA, the initial size of the flame.
%
      u = x(j);
      factor = 2^u;
      delta = factor * delta_base;
%
%  Set the starting point.
%
      t_start = 0.0;
      y_start(1,1) = delta;
%
%  Get the stopping point.
%
      t_stop = 250.0;
%
%  Call ODE45 to solve the problem, stopping immediately if the event
%  (y=0.99) is observed.
%
      t_span = [ t_start, t_stop ];

      [ t, y ] = ode45 ( @flame_fun, t_span, y_start, options );
%
%  Our function must be weighted by rho(u), which is uniformly 1/2.
%
      f(j) = 0.5 * t(end);

    end
%
%  Compute the integral estimate.
%
    q = w' * f;

    fprintf ( 1, '  %2d  %12g\n', n, q );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QOI_QUAD\n' );
  fprintf ( 1, '  Normal end of execution\n' );

  return
end


text:
% @TTENSOR
%
% Files
%   disp      - Command window display of a ttensor.
%   display   - Command window display of a ttensor.
%   double    - Convert ttensor to double array.
%   end       - Last index of indexing expression for ttensor.
%   full      - Convert a ttensor to a (dense) tensor.
%   innerprod - Efficient inner product with a ttensor.
%   isequal   - True if each component of two ttensor's is numerically equal.
%   mtimes    - Implement scalar multiplication for a ttensor.
%   mttkrp    - Matricized tensor times Khatri-Rao product for ttensor.
%   ndims     - Return the number of dimensions for a ttensor.
%   norm      - Norm of a ttensor.
%   nvecs     - Compute the leading mode-n vectors for a ttensor.
%   permute   - Permute dimensions for a ttensor.
%   size      - Size of a ttensor.
%   subsasgn  - Subscripted reference for a ttensor.
%   subsref   - Subscripted reference for a ttensor.
%   ttensor   - Tensor stored as a Tucker operator (decomposed).
%   ttm       - Tensor times matrix for ttensor.
%   ttv       - Tensor times vector for ttensor.
%   uminus    - Unary minus for ttensor.
%   uplus     - Unary plus for ttensor.

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function calpak_test70 ( )

%*****************************************************************************80
%
%% TEST70 tests YMDF_INC_COMMON, YMDF_NEXT_COMMON, YMDF_PREV_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2013
%
%  Author:
%
%    John Burkardt
%
  days = 10.25;
  dhi = 1;
  dlo = 1;
  fhi = 0.0;
  flo = 0.0;
  mhi = 1;
  mlo = 1;
  seed = 123456789;
  yhi = 1960;
  ylo = 1970;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST70\n' );
  fprintf ( 1, '  For the Common calendar:\n' );
  fprintf ( 1, '  YMDF_INC_COMMON increments a date by days;\n' );
  fprintf ( 1, '  YMDF_NEXT_COMMON computes the next day,\n' );
  fprintf ( 1, '  YMDF_PREV_COMMON computes the previous day.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMDF      Tomorrow       Yesterday      +10.25 days\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ y1, m1, d1, f1, seed ] = ymdf_uniform_common ( ylo, mlo, dlo, flo, ...
      yhi, mhi, dhi, fhi, seed );

    s1 = ymdf_to_s_common ( y1, m1, d1, f1 );

    [ y2, m2, d2, f2 ] = ymdf_next_common ( y1, m1, d1, f1 );
    s2 = ymdf_to_s_common ( y2, m2, d2, f2 );

    [ y3, m3, d3, f3 ] = ymdf_prev_common ( y1, m1, d1, f1 );
    s3 = ymdf_to_s_common ( y3, m3, d3, f3 );

    [ y4, m4, d4, f4 ] = ymdf_inc_common ( y1, m1, d1, f1, days );
    s4 = ymdf_to_s_common ( y4, m4, d4, f4 );

    fprintf ( 1, '  %15s  %15s  %15s  %15s\n', s1, s2, s3, s4 );

  end

  return
end


text:
function calculateDDM(basePlanCFile,ddmTemplatePlanCFileName,...
    registeredPlanCFilesDir,vfDir,saveFileLocationWithDDM)

% Linux
basePlanCFile = '/lab/deasylab1/Data/RTOG0617/CERR_files_tcia/rider_template/RIDER-1225316081_First_resampled_1x1x3.mat';
ddmTemplatePlanCFileName = '0617-548359_09-09-2000-30377';
registeredPlanCFilesDir = '/lab/deasylab1/Data/RTOG0617/CERR_files_tcia/dose_mapping_original_plans'; % RTOG0617
vfDir = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/ddm/RTOG0617_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
saveFileLocationWithDDM = '/lab/deasylab1/Data/RTOG0617/registrations_pericardium/ddm/RTOG0617_to_RIDER_1225316081_First_template_ddm.mat';

% Windows
% basePlanCFile = 'L:/Data/RTOG0617/CERR_files_tcia/rider_template/RIDER-1225316081_First_resampled_1x1x3.mat';
% ddmTemplatePlanCFileName = '0617-548359_09-09-2000-30377';
% registeredPlanCFilesDir = 'L:/Data/RTOG0617/CERR_files_tcia/dose_mapping_original_plans'; % RTOG0617
% vfDir = 'L:/Data/RTOG0617/registrations_pericardium/ddm/RTOG0617_to_RIDER_1225316081_First_template'; %RIDER RIDER-1225316081
% saveFileLocationWithDDM = 'L:/Data/RTOG0617/registrations_pericardium/ddm/RTOG0617_to_RIDER_1225316081_First_template_ddm.mat';
%
% basePlanCFile = strrep(basePlanCFile,'/','\');
% ddmTemplatePlanCFileName = strrep(ddmTemplatePlanCFileName,'/','\');
% registeredPlanCFilesDir = strrep(registeredPlanCFilesDir,'/','\');
% vfDir = strrep(vfDir,'/','\');
% saveFileLocationWithDDM = strrep(saveFileLocationWithDDM,'/','\');

[~,basePlanCFileName] = fileparts(basePlanCFile);

% registeredPlanCFilesDir = 'L:\Data\RTOG0617\registrations_pericardium\ddm\RTOG0617_to_RIDER_1225316081_First_template'

% registeredPlanCFilesDir = '/lab/deasylab2/Ishita/Stenosis_ultracentral/00216080_unzipped';
% baseTemplateFile = 'Baseline';
% ddmTemplateFile = '6month';

% baseTemplateFile = fullfile(registeredPlanCFilesDir,[basePlanCFile,'.mat']) ;
% ddmTemplateFile = fullfile(registeredPlanCFilesDir,[ddmTemplatePlanCFile,'.mat']) ;

dirS = dir(registeredPlanCFilesDir);
dirS([dirS.isdir]) = [];
fileNamC = {dirS.name};

% x,y,z coordinates for DDM calculation
planC = loadPlanC(basePlanCFile,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(basePlanCFile,planC);
scanNum = 1;
indexS = planC{end};
siz = getUniformScanSize(planC{indexS.scan}(scanNum));
[xBaseV,yBaseV,zBaseV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
[xBase3M,yBase3M,zBase3M] = meshgrid(xBaseV,yBaseV,zBaseV);
% xBaseV = xBase3M(:);
% yBaseV = yBase3M(:);
% zBaseV = zBase3M(:);

% filesC = fullfile(registeredPlanCFilesDir,fileNamC);

% if ~exist('vfDir','var')
%     vfDir = fullfile(registeredPlanCFilesDir,'registered');
% end

numFiles = length(fileNamC);
xDeformM = [];
yDeformM = [];
zDeformM = [];
for iFile = 1:numFiles
    fName = fileNamC{iFile};
    fullFilePath = fullfile(registeredPlanCFilesDir,fName);
    planC = loadPlanC(fullFilePath,tempdir);
    indexS = planC{end};
    scanNum = 1;
    [xUnifV,yUnifV,zUnifV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));

    if strcmp(basePlanCFileName,strtok(fName,'.'))
        continue;
    end
    vfBaseName = [basePlanCFileName,'~',strtok(fName,'.'),'_vf.mat'];
    vfBaseFile = fullfile(vfDir,vfBaseName);
    if ~exist(vfBaseFile,'file')
        continue;
    end

    load(vfBaseFile)
    xBaseDeform3M = xBase3M + vf(:,:,:,1);
    yBaseDeform3M = yBase3M + vf(:,:,:,2);
    zBaseDeform3M = zBase3M + vf(:,:,:,3);

    if strcmp(strtok(fName,'.'),ddmTemplatePlanCFileName)
        xDeformV = xBaseDeform3M(:);
        yDeformV = yBaseDeform3M(:);
        zDeformV = zBaseDeform3M(:);
    else
        vfDdmName = [strtok(fName,'.'),'~',ddmTemplatePlanCFileName,'_vf.mat'];
        vfDdmFile = fullfile(vfDir,vfDdmName);
        if ~exist(vfDdmFile,'file')
            continue;
        end
        load(vfDdmFile)

        xFieldV = [xUnifV(1)-1e-6,xUnifV(2)-xUnifV(1),xUnifV(end)+1e-6];
        yFieldV = [yUnifV(end)-1e-6,yUnifV(1)-yUnifV(2),yUnifV(1)+1e-6];

        xDeformV = xBaseDeform3M(:) + finterp3(xBaseDeform3M(:),yBaseDeform3M(:),zBaseDeform3M(:),...
            flip(vf(:,:,:,1),1),xFieldV,yFieldV,zUnifV);
        yDeformV = yBaseDeform3M(:) + finterp3(xBaseDeform3M(:),yBaseDeform3M(:),zBaseDeform3M(:),...
            flip(vf(:,:,:,2),1),xFieldV,yFieldV,zUnifV);
        zDeformV = zBaseDeform3M(:) + finterp3(xBaseDeform3M(:),yBaseDeform3M(:),zBaseDeform3M(:),...
            flip(vf(:,:,:,3),1),xFieldV,yFieldV,zUnifV);
    end

    xDeformM(:,end+1) = xDeformV;
    yDeformM(:,end+1) = yDeformV;
    zDeformM(:,end+1) = zDeformV;

end

% % Pairwise distance
% numDeforms = size(xDeformM,2);
% pairsM = nchoosek(1:numDeforms,2);
% numPairs = size(pairsM,1);
% distM = zeros(size(xDeformM,1),numPairs);
% for i = 1:numPairs
%     xSquareV = diff(xDeformM(:,pairsM(i,:)),1,2).^2;
%     ySquareV = diff(yDeformM(:,pairsM(i,:)),1,2).^2;
%     zSquareV = diff(zDeformM(:,pairsM(i,:)),1,2).^2;
%     distSquareV = xSquareV + ySquareV + zSquareV;
%     distM(:,i) = distSquareV.^0.5;
% end
% medianDistV = nanmedian(distM,2);

xMeanV = median(xDeformM,2,'omitnan');
yMeanV = median(yDeformM,2,'omitnan');
zMeanV = median(zDeformM,2,'omitnan');
xSquareV = bsxfun(@minus,xDeformM,xMeanV);
ySquareV = bsxfun(@minus,yDeformM,yMeanV);
zSquareV = bsxfun(@minus,zDeformM,zMeanV);

medianDistV = median((xSquareV.^2 + ySquareV.^2 + zSquareV.^2).^0.5,2,...
    'omitnan')/sqrt(2);

dist3M = reshape(medianDistV,siz);

% Add as dose to base planC
planC = loadPlanC(basePlanCFile,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(basePlanCFile,planC);
indexS = planC{end};
scanType = 'UniformCT';
assocScanUID = planC{indexS.scan}(scanNum).scanUID;
planC = dose2CERR(single(dist3M),'','DDM','',...
'',scanType,'','No',assocScanUID,planC);
planC = save_planC(planC,[],'passed',saveFileLocationWithDDM);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
clear;
gcp;    % start a local cluster

filename = 'demoSue2x.tif';
if ~exist(filename,'file');
    url = 'https://www.dropbox.com/s/36xdfd28eone0hj/demoSue2x.tif?dl=1';
    fprintf('downloading the file...');
    outfilename = websave(filename,url);
    fprintf('done. \n');
end

addpath(genpath('../../../ca_source_extraction-master'));  % add packages to matlab path
addpath(genpath('../../../NoRMCorre-master'));
addpath(genpath('../../../ca_source_extraction'));  % add packages to matlab path
addpath(genpath('../../../NoRMCorre'));

%% read file and determine dynamic range

Y = read_file(filename);
[d1,d2,T] = size(Y);    % dimensions of file
Y = Y - min(Y(:));      % remove negative offset

minY = quantile(Y(1:1e7),0.0005);
maxY = quantile(Y(1:1e7),1-0.0005);

%%  view data
figure;play_movie({Y},{'raw data'},minY,maxY);

%% perform motion correction (start with rigid)
% parameters motion correction
% 'd1','d2': size of FOV
% 'bin_width': how often to update the template
% 'max_shift': maximum allowed rigid shift

options_rg = NoRMCorreSetParms('d1',size(Y,1),'d2',size(Y,2),'bin_width',100,'max_shift',15);

[M_rg,shifts_rg,template_rg] = normcorre_batch(Y,options_rg);

%% view data
tsub = 5;   % downsampling factor (only for display purposes)
Y_sub = downsample_data(Y,'time',tsub);
M_rgs = downsample_data(M_rg,'time',tsub);

play_movie({Y_sub,M_rgs},{'raw data','rigid'},minY,maxY);

%% perform non-rigid motion correction
% parameters motion correction
% 'd1','d2': size FOV movie
% 'grid_size','overlap_pre': parameters regulating size of patch (size patch ~ (grid_size + 2*overlap_pre))
% 'mot_uf': upsampling factor of the grid for shift application
% 'bin_width': how often to update the template
% 'max_shift': maximum allowed rigid shift
% 'max_dev': maximum deviation allowed for each patch from the rigid shift value

options_nr = NoRMCorreSetParms('d1',size(Y,1),'d2',size(Y,2),...
                    'grid_size',[48,48],'mot_uf',4,'overlap_pre',[16,16],...
                    'bin_width',100,'max_shift',15,'max_dev',8);

[M_nr,shifts_nr,template_nr] = normcorre_batch(Y,options_nr,template_rg);

%% view (downsampled) data

M_nrs = downsample_data(M_nr,'time',tsub);
play_movie({Y_sub,M_rgs,M_nrs},{'raw data','rigid','pw-rigid'},minY,maxY);

%% compute some metrics for motion correction quality assessment

[cY,mY,vY] = motion_metrics(Y,options_rg.max_shift);
[cM_rg,mM_rg,vM_rg] = motion_metrics(M_rg,options_rg.max_shift);
[cM_nr,mM_nr,vM_nr] = motion_metrics(M_nr,options_rg.max_shift);

%% plot shifts

shifts_r = squeeze(cat(3,shifts_rg(:).shifts));
shifts_n = cat(ndims(shifts_nr(1).shifts)+1,shifts_nr(:).shifts);
shifts_n = reshape(shifts_n,[],ndims(Y)-1,T);
shifts_x = squeeze(shifts_n(:,2,:))';
shifts_y = squeeze(shifts_n(:,1,:))';

patch_id = 1:size(shifts_x,2);
str = strtrim(cellstr(int2str(patch_id.')));
str = cellfun(@(x) ['patch # ',x],str,'un',0);

figure;
    ax1 = subplot(311); plot(1:T,cY,1:T,cM_rg,1:T,cM_nr); legend('raw data','rigid','non-rigid'); title('correlation coefficients','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[])
    ax2 = subplot(312); plot(shifts_x); hold on; plot(shifts_r(:,2),'--k','linewidth',2); title('displacements along x','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[])
    ax3 = subplot(313); plot(shifts_y); hold on; plot(shifts_r(:,1),'--k','linewidth',2); title('displacements along y','fontsize',14,'fontweight','bold')
            xlabel('timestep','fontsize',14,'fontweight','bold')
    linkaxes([ax1,ax2,ax3],'x')

%% plot metrics
figure;
    ax(1) = subplot(2,3,1); imagesc(mY,[minY,maxY]);  axis equal; axis tight; axis off; title('mean raw data','fontsize',14,'fontweight','bold')
    ax(2) = subplot(2,3,2); imagesc(mM_rg,[minY,maxY]); axis equal; axis tight; axis off; title('mean rigid corrected','fontsize',14,'fontweight','bold')
    ax(3) = subplot(2,3,3); imagesc(mM_nr,[minY,maxY]); axis equal; axis tight; axis off; title('mean non-rigid corrected','fontsize',14,'fontweight','bold')
    subplot(2,3,4); plot(1:T,cY,1:T,cM_rg,1:T,cM_nr); legend('raw data','rigid','non-rigid'); title('correlation coefficients','fontsize',14,'fontweight','bold')
    subplot(2,3,5); scatter(cY,cM_rg); hold on; plot([0.9*min(cY),1.05*max(cM_rg)],[0.9*min(cY),1.05*max(cM_rg)],'--r'); axis square;
        xlabel('raw data','fontsize',14,'fontweight','bold'); ylabel('rigid corrected','fontsize',14,'fontweight','bold');
    subplot(2,3,6); scatter(cM_rg,cM_nr); hold on; plot([0.95*min(cM_rg),1.05*max(cM_nr)],[0.95*min(cM_rg),1.05*max(cM_nr)],'--r'); axis square;
        xlabel('rigid corrected','fontsize',14,'fontweight','bold'); ylabel('non-rigid corrected','fontsize',14,'fontweight','bold');
    linkaxes(ax,'xy')

%% now perform source extraction by splitting the FOV in patches

sizY = size(M_nr);
patch_size = [30,30];                   % size of each patch along each dimension (optional, default: [32,32])
overlap = [8,8];                        % amount of overlap in each dimension (optional, default: [4,4])

patches = construct_patches(sizY(1:end-1),patch_size,overlap);
K = 4;                                            % number of components to be found
tau = 4;                                          % std of gaussian kernel (half size of neuron)
p = 0;                                            % order of autoregressive system (p = 0 no dynamics, p=1 just decay, p = 2, both rise and decay)

options = CNMFSetParms(...
    'd1',sizY(1),'d2',sizY(2),...
    'temporal_iter',2,...                       % number of block-coordinate descent steps
    'ssub',1,...                                % downsample in space
    'tsub',2,...                                % downsample in time
    'merge_thr',0.8,...                         % merging threshold
    'gSig',tau,...
    'gnb',2,...                                 % number of background components
    'spatial_method','regularized'...
    );

%% run CNMF algorithm on patches and combine
tic;
[A,b,C,f,S,P,RESULTS,YrA] = run_CNMF_patches(M_nr,K,patches,tau,p,options);
[ROIvars.rval_space,ROIvars.rval_time,ROIvars.max_pr,ROIvars.sizeA,keep] = classify_components(M_nr,A,C,b,f,YrA,options);
toc

%% a simple GUI
Cn = correlation_image_max(M_nr);
Coor = plot_contours(A,Cn,options,1); close;
run_GUI = false;
if run_GUI
    GUIout = ROI_GUI(A,options,Cn,Coor,keep,ROIvars);
    options = GUIout{2};
    keep = GUIout{3};
end

%% view contour plots of selected and rejected components (optional)
throw = ~keep;
figure;
    ax1 = subplot(121); plot_contours(A(:,keep),Cn,options,0,[],Coor,1,find(keep)); title('Selected components','fontweight','bold','fontsize',14);
    ax2 = subplot(122); plot_contours(A(:,throw),Cn,options,0,[],Coor,1,find(throw));title('Rejected components','fontweight','bold','fontsize',14);
    linkaxes([ax1,ax2],'xy')
%% inspect components

plot_components_GUI(M_nr,A(:,keep),C(keep,:),b,f,Cn,options);

%% refine temporal components
A_keep = A(:,keep);
C_keep = C(keep,:);
P.p = 0;
[C2,f2,P2,S2,YrA2] = update_temporal_components(reshape(M_nr,[],T),A_keep,b,C_keep,f,P,options);

%% detrend fluorescence and extract DF/F values

options.df_window = 1000;
[F_dff,F0] = detrend_df_f(A_keep,b,C2,f2,YrA2,options);

%% deconvolve data

nNeurons = size(F_dff,1);
C_dec = zeros(size(F_dff));
S = zeros(size(F_dff));
kernels = cell(nNeurons,1);
min_sp = 3;    % find spikes resulting in transients above min_sp x noise level

for i = 1:nNeurons
    [C_dec(i,:),S(i,:),kernels{i}] = deconvCa(F_dff(i,:), [], min_sp, true, false, [], 20, [], 0);
end

%% plot a random component

i = randi(nNeurons);

figure;plot(1:T,F_dff(i,:),'--k'); hold all; plot(1:T,C_dec(i,:),'r','linewidth',2);
    spt = find(S(i,:));
    if spt(1) == 1; spt(1) = []; end
    hold on; scatter(spt,repmat(-0.25,1,length(spt)),'m*')
    title(['Component ',num2str(i)]);

    legend('Fluorescence DF/F','Deconvolved','Spikes')

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function masks = char2img( strings, h, pad )
% Convert ascii text to a binary image using pre-computed templates.
%
% Input strings can only contain standard characters (ascii character
% 32-126). First time char2img() is ever called with a given height h,
% txt2img() is used to create a template for each ascii character. All
% subsequent calls to to char2img() with the given height are very fast as
% the pre-computed templates are used and no display/screen capture is
% needed.
%
% USAGE
%  masks = char2img( strings, h, [pad] )
%
% INPUTS
%  strings  - {1xn} text string or cell array of text strings to convert
%  h        - font height in pixels
%  pad      - [0] amount of extra padding between chars
%
% OUTPUTS
%  masks    - {1xn} binary image masks of height h for each string
%
% EXAMPLE
%  mask=char2img('hello world',50); im(mask{1})
%  mask=char2img(num2str(pi),35); im(mask{1})
%
% See also TXT2IMG
%
% Piotr's Image&Video Toolbox      Version 2.65
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% load or create character templates (or simply store persistently)
persistent chars h0
if(isempty(h0) || h0~=h)
  pth=fileparts(mfilename('fullpath'));
  fName=sprintf('%s/private/char2img_h%03i.mat',pth,h); h0=h;
  if(exist(fName,'file')), load(fName); else
    chars=char(32:126); chars=num2cell(chars);
    chars=txt2img(chars,h,{'Interpreter','none'});
    save(fName,'chars');
  end
end
% add padding to chars
if(nargin<3 || isempty(pad)), pad=0; end
charsPad=chars; if(pad), pad=ones(h,pad,'uint8');
  for i=1:length(chars), charsPad{i}=[pad chars{i} pad]; end; end
% create actual string using templates
if(~iscell(strings)), strings={strings}; end
n=length(strings); masks=cell(1,n);
for s=1:n, str=strings{s};
  str(str<32 | str>126)=32; str=uint8(str-31);
  M=[charsPad{str}]; masks{s}=M;
end

author: zouchuhang

repo: LayoutNet


text:
function ei_values_test ( )

%*****************************************************************************80
%
%% EI_VALUES_TEST demonstrates the use of EI_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EI_VALUES_TEST:\n' );
  fprintf ( 1, '  EI_VALUES stores values of\n' );
  fprintf ( 1, '  the exponential integral function EI(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X          EI(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = ei_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function at = auct(crit,y,z,tt)
%AUCS Compute area under curve for survival model at given time
%
%  Description
%    A = AUCT(CRIT,Z) Compute are under curve for survival model at
%    given time using criteria vector CRIT (where larger value
%    means larger risk of incidence), observed time vector Y,
%    censoring indicator matrix Z at times TT (0=event, 1=censored)
%    and a time vector TT.
%
%  Reference
%    L. E. Chambless, C. P. Cummiskey, and G. Cui (2011). Several
%    methods to assess improvement in risk prediction models:
%    Extension to survival analysis. Statistics in Medicine
%    30(1):22-38.

% Copyright (C) 2012 Ernesto Ulloa, Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

ip=inputParser;
ip.addRequired('crit',@(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
ip.addRequired('y', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('z', @(x) isreal(x) && all(isfinite(x(:))))
ip.addRequired('tt', @(x) isreal(x) && all(isfinite(x(:))))
ip.parse(crit,y,z,tt)

for i=1:size(tt,2)
  comp=bsxfun(@times,bsxfun(@and,y(:,i)<=tt(i),1-z(:,i)),bsxfun(@or,bsxfun(@and,y(:,i)<=tt(i),z(:,i)),y(:,i)>=tt(i))');
  conc=bsxfun(@times,bsxfun(@gt,crit(:,i),crit(:,i)'),comp);
  at(i,1)=sum(conc(:))./sum(comp(:));
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Number of measurements
M = 64;
% Sparsity level
K = 8;
% Number of signals
S = 4;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
Y = Phi * X;

solver = spx.pursuit.joint.BasisPursuit(Phi);
result = solver.solve_linf_l1(Y);
% Solution vectors
Z = result.Z;

solver = spx.pursuit.joint.BasisPursuit(Phi);
result = solver.solve_l1_l1(Y);
% Solution vectors
Z = result.Z;

for s=1:S
    mf.new_figure(sprintf('MMV signal: %d', s));
    subplot(411);
    stem(X(:, s), '.');
    title('Sparse vector');
    subplot(412);
    stem(Z(:, s), '.');
    title('Recovered sparse vector');
    subplot(413);
    stem(abs(X(:, s) - Z(:, s)), '.');
    title('Recovery error');
    subplot(414);
    stem(Y(:, s), '.');
    title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function out = dbv(in)

out = 20 * log10(abs(in));

author: lukeweston

repo: SimpleFMCWRadar


text:
function flowplot(qmethod,sol,By,Bx,A,xy,xyp,x,y,bound,spc,fig)
%flowplot   plots flow data on square-shaped domain
%   flowplot(qmethod,sol,By,Bx,A,xy,xyp,x,y,bound,spc,fig);
%   input
%          qmethod    mixed method
%          sol        flow solution vector
%          By         velocity  y-derivative matrix
%          Bx         velocity x-derivative matrix
%          A          vector diffusion matrix
%          xy         velocity nodal coordinate vector
%          xyp        pressure nodal coordinate vector
%          x          vector of x-axis interpolation points
%          y          vector of y-axis interpolation points
%          bound      boundary vertex vector
%          spc        uniform/nonuniform streamline switch
%          fig        figure number
%
%   IFISS function: DJS; 30 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nvtx=length(xy); nu=2*nvtx; np=length(xyp);
Asv=A(1:nvtx,1:nvtx);
%
% compute auxilliary quantites
u=sol(1:nu);p=sol(nu+1:end);
f=[By,-Bx]*u;
[Asv,fsv] = nonzerobc(Asv,f,xy,bound);
phi=Asv\fsv;
%
%% plot pressure
if qmethod==2
   xx=x(1:2:end); yy=y(1:2:end);
elseif qmethod==3
   p=p(1:3:end); xx=x(1:end); yy=y(1:end);
else
   xx=x(1:end); yy=y(1:end);
end
% interpolate to a cartesian product mesh
[X,Y]=meshgrid(xx,yy);
xysol = griddata(xyp(:,1),xyp(:,2),p,X,Y);
figure(fig)
subplot(122), mesh(X,Y,xysol),axis('square')
title('pressure field')
%
%% plot velocity
[X,Y]=meshgrid(x,y);
xysol = griddata(xy(:,1),xy(:,2),phi,X,Y);
if spc==1,
   subplot(121),contour(X,Y,xysol,20),axis('square')
   title('Streamlines: uniform');
   if all([min(x),max(x),min(y),max(y)] == [-1,1,-1,1]), squarex, end
elseif spc==2,
   v=[-15:1.3:6]';v=exp(v); list=sort([-v;v]);
   subplot(121),contour(X,Y,xysol,list),axis('square')
   title('Streamlines: selected');
   if all([min(x),max(x),min(y),max(y)] == [-1,1,-1,1]), squarex, end
end
axis('off')
subplot(111)
return


text:
%% parLimNat
% Below is a demonstration of the features of the |parLimNat| function

%% Syntax
% |[xx,S]=parLimNat(xx_c,[xx_min xx_max],x);|

%% Description
% The |parLimNat| function can be used to constrain parameters from
% [-inf,inf] to the range [xx_min xx_max] with the values xx_c at its
% centre. The constraining can be used in combination with optimization
% routinges that do not naturally handle parameter constraints.

%% Examples

clear; close all; clc;

%%
% PLOT SETTINGS
fontSize=25;
fontSize2=15;
markerSize=45;
lineWidth1=2;
lineWidth2=1;

%% Example: Constraining parameters (normal centre)
xx_c=5;
xx_min=0;
xx_max=10;
x=linspace(xx_c-10,xx_c+10,1000);
[xx,S]=parLimNat(xx_c,[xx_min xx_max],x);

cFigure; hold on; grid on;
title('Constraining using parLimNat','FontSize',fontSize);
xlabel('"free" x','FontSize',fontSize); ylabel('constrained x','FontSize',fontSize);

plot(x,xx,'r-','LineWidth',lineWidth1);
plot(xx_min,xx_min,'k.','markerSize',markerSize);
plot(xx_max,xx_max,'k.','markerSize',markerSize);
plot(xx_c,xx_c,'k.','markerSize',markerSize);

text(xx_c+0.5,xx_c,'xx_c = centre','Interpreter','none','FontSize',fontSize2);
text(xx_max,xx_max-0.5,'xx_max = upper bound','Interpreter','none','FontSize',fontSize2);
text(xx_min-2.5,xx_min+0.5,'xx_min = lower bound','Interpreter','none','FontSize',fontSize2);

axis tight; axis equal;
set(gca,'FontSize',fontSize);
drawnow;

%% Example: Constraining parameters (out of centre centre)
xx_c=2;
xx_min=0;
xx_max=10;
x=linspace(xx_c-10,xx_c+15,100);
[xx,S]=parLimNat(xx_c,[xx_min xx_max],x);

cFigure; hold on; grid on;
title('Constraining using parLimNat','FontSize',fontSize);
xlabel('"free" x','FontSize',fontSize); ylabel('constrained x','FontSize',fontSize);

plot(x,xx,'r-','LineWidth',lineWidth1);
plot(xx_min,xx_min,'k.','markerSize',markerSize);
plot(xx_max,xx_max,'k.','markerSize',markerSize);
plot(xx_c,xx_c,'k.','markerSize',markerSize);

text(xx_c+0.5,xx_c,'xx_c = centre','Interpreter','none','FontSize',fontSize2);
text(xx_max,xx_max-0.5,'xx_max = upper bound','Interpreter','none','FontSize',fontSize2);
text(xx_min-2.5,xx_min+0.5,'xx_min = lower bound','Interpreter','none','FontSize',fontSize2);

axis tight; axis equal;
set(gca,'FontSize',fontSize);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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


text:
function stp
global st
st=true;


text:
function test_failed = test_uwpfbt(verbose)
%TEST_WFBTPR
%
% Checks perfect reconstruction of the general wavelet transform of different
% filters
%
disp('========= TEST UWPFBT ============');
global LTFAT_TEST_TYPE;
tolerance = 1e-8;
if strcmpi(LTFAT_TEST_TYPE,'single')
   tolerance = 2e-6;
end

test_failed = 0;
if(nargin>0)
   verbose = 1;
else
   verbose = 0;
end

type = {'dec'};
ext = {'per'};

J = 4;
%! Mild tree
wt1 = wfbtinit({'db10',6,'full'});
wt1 = wfbtremove(1,1,wt1,'force');
wt1 = wfbtremove(2,1,wt1,'force');

%! Hardcore tree
wt2 = wfbtinit({'db3',1});
wt2 = wfbtput(1,1,'mband1',wt2);
wt2 = wfbtput(2,2,'mband1',wt2);
wt2 = wfbtput(3,3,'mband1',wt2);
wt2 = wfbtput(3,1,'db10',wt2);
wt2 = wfbtput(4,1,'dgrid2',wt2);
wt2 = wfbtput(5,1,'db2',wt2);

% wt2 = wfbtinit();
% wt2 = wfbtput(0,0,{'db',4},wt2);
% wt2 = wfbtput(1,0,{'algmband',1},wt2);
% wt2 = wfbtput(1,1,{'hden',3},wt2);
% wt2 = wfbtput(2,0,{'dgrid',2},wt2);
% wt2 = wfbtput(2,1,{'dgrid',2},wt2);

test_filters = {

               {'algmband1',J} % 3 filters, uniform, crit. sub.
               {'algmband2',J} % 4 filters, uniform, crit. sub.
               {'db10',J}
               %{{'hden',3},J} % 3 filters, non-uniform, no crit. sub. no correct
               {'dgrid1',J} % 4 filters. sub. fac. 2
               wt1
               wt2
               };

scaling = {'scale','sqrt','noscale'};
scalingInv = scaling(end:-1:1);

interscaling = {'intscale','intsqrt','intnoscale'};
interscalingInv = interscaling(end:-1:1);

for scIdx = 1:numel(scaling)
for iscIdx = 1:numel(interscaling)
%testLen = 4*2^7-1;%(2^J-1);
testLen = 53;
f = tester_rand(testLen,1);

for extIdx=1:length(ext)
   extCur = ext{extIdx};

   for typeIdx=1:length(type)
     for tt=1:length(test_filters)
        actFilt = test_filters{tt};
         if verbose, if(~isstruct(actFilt))fprintf('J=%d, filt=%s, ext=%s, inLen=%d \n',actFilt{2},actFilt{1},extCur,length(f)); else disp('Custom'); end; end;

        [c,info] = uwpfbt(f,actFilt,scaling{scIdx},interscaling{iscIdx});
        fhat = iuwpfbt(c,actFilt,scalingInv{scIdx},interscalingInv{iscIdx});

            err = norm(f-fhat,'fro');
            [test_failed,fail]=ltfatdiditfail(err,test_failed,tolerance);
            if(~verbose)
              if(~isstruct(actFilt))
                  fprintf('J=%d, %5.5s, ext=%4.4s, %s, %s, L=%d, err=%.4e %s \n',actFilt{2},actFilt{1},extCur,scaling{scIdx},interscaling{iscIdx},size(f,1),err,fail);
              else
                  fprintf('Custom, %s, %s, err=%.4e %s\n',scaling{scIdx},interscaling{iscIdx},err,fail);
              end;
            end
            if strcmpi(fail,'FAILED')
               if verbose
                 if(~isstruct(actFilt)) fprintf('err=%d, filt=%s, ext=%s, inLen=%d \n',err,actFilt{1},extCur,testLen);else disp('Fail. Custom'); end;
                 figure(1);clf;stem([f,fhat]);
                 figure(2);clf;stem([f-fhat]);
                 break;
               end
            end

            fhat2 = iuwpfbt(c,info);
            err = norm(f-fhat2,'fro');
            [test_failed,fail]=ltfatdiditfail(err,test_failed,tolerance);

            if(~isstruct(actFilt))
                  fprintf('INFO J=%d, %5.5s, ext=%4.4s, %s, %s, L=%d, err=%.4e %s \n',actFilt{2},actFilt{1},extCur,scaling{scIdx},interscaling{iscIdx},size(f,1),err,fail);
            else
                  fprintf('INFO Custom, %s, %s, err=%.4e %s\n',scaling{scIdx},interscaling{iscIdx},err,fail);
            end;

            if test_failed && verbose, break; end;

     end
     if test_failed && verbose, break; end;
   end
   if test_failed && verbose, break; end;
end
end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function A = calc_area(I)
% The derivative uses central difference and only applies to the interior
% points
der_x = calc_der_x(I(2:end-1,:,:));
der_y = calc_der_y(I(:,2:end-1,:));
G11 = sum(der_x.^2, 3) + 1;
G12 = sum(der_x.*der_y, 3);
G22 = sum(der_y.^2, 3) + 1;
A = sum(sum(sqrt(G11.*G22 - G12.^2)));
% A = sum(sum(G11.*G22 - G12.^2));

function der_x = calc_der_x(I)
der_x = (1/2) * (I(:,3:end,:) - I(:,1:end-2,:));

function der_y = calc_der_y(I)
der_y = (1/2) * (I(3:end,:,:) - I(1:end-2,:,:));

author: zhouyuanzxcv

repo: Hyperspectral


text:
function test_suite=test_check_external()
% tests for cosmo_check_external
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_check_external_nifti()
    % nifti comes with CoSMoMVPA, so should always be available if the path
    % is set properly

    warning_state=cosmo_warning();
    orig_path=path();

    warning_state_resetter=onCleanup(@()cosmo_warning(warning_state));
    path_resetter=onCleanup(@()path(orig_path));

    % ensure path is set; disable warnings by cosmo_set_path
    cosmo_warning('off');
    cosmo_set_path();

    assertTrue(cosmo_check_external('nifti'))
    assertEqual(cosmo_check_external({'nifti','nifti'}),[true;true]);

    % test list
    externals=cosmo_check_external('-list');
    assert(~isempty(strmatch('nifti',externals,'exact')));

    % test tic/toc
    cosmo_check_external('-tic');
    assert(isempty(cosmo_check_external('-toc')));

    % must cite CoSMoMVPA, but not NIFTI
    c=cosmo_check_external('-cite');
    assert(~isempty(findstr(c,'N. N. Oosterhof')));
    assert(~isempty(findstr(c,'CoSMoMVPA')));
    assert(isempty(findstr(c,'NIFTI toolbox')));

    % after checking for nifti it must be present
    cosmo_check_external('nifti');
    assertEqual(cosmo_check_external('-toc'),{'nifti'});
    c=cosmo_check_external('-cite');
    assert(~isempty(findstr(c,'NIFTI toolbox')));

    % removing nifti from the path makes it absent
    rmpath(fileparts(which('load_nii')));
    assertFalse(cosmo_check_external('nifti',false));

    % throws exception when second argument is true or absent
    assertExceptionThrown(@()cosmo_check_external('nifti',true),'');
    assertExceptionThrown(@()cosmo_check_external('nifti'),'');

function test_check_external_mocov()
    has_mocov=~isempty(which('mocov'));
    assertEqual(has_mocov,cosmo_check_external('mocov',false));

function test_check_external_command()
    commands={'foo','basdfds','disp'};
    n=numel(commands);
    for k=1:n
        command=commands{k};
        has_command=~isempty(which(command));

        arg=['!' command];
        if has_command
            assertTrue(cosmo_check_external(arg));
            assertTrue(cosmo_check_external(arg,true));
            assertTrue(cosmo_check_external(arg,false));
        else
            assertExceptionThrown(@()cosmo_check_external(arg),'');
            assertExceptionThrown(@()cosmo_check_external(arg,true),'');
            assertFalse(cosmo_check_external(arg,false));
        end
    end

function test_check_external_exceptions()
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_check_external(varargin{:}),'');
    aet('unknown package');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
%
% Goal: Create canonical spherical harmonic bases
%
% Li Shen
% 04/11/2002 - create
% 10/15/2002 - rename and modify
% 11/03/2008 - renamed by Sungeun Kim.

function Z = calculate_SPHARM_basis(vs, degree)

[PHI,THETA] = cart2sph(vs(:,1),vs(:,2),vs(:,3));
ind = find(PHI<0);
PHI(ind) = PHI(ind)+2*pi;
THETA = pi/2-THETA;
vertnum = size(THETA,1);

Z = spharm_basis(degree,THETA,PHI);

return;

function Z = spharm_basis(max_degree,theta,phi)

Z = []; vnum = size(theta,1);

% save calculations for efficiency
for k = 0:(2*max_degree)
    fact(k+1) = factorial(k);
end
for m = 0:max_degree
    exp_i_m_phi(:,m+1) = exp(i*m*phi);
    sign_m(m+1) = (-1)^(m);
end

for n = 0:max_degree

	% P = legendre(n,X) computes the associated Legendre functions of degree n and
	% order m = 0,1,...,n, evaluated at X. Argument n must be a scalar integer
	% less than 256, and X must contain real values in the domain -1<=x<=1.
	% The returned array P has one more dimension than X, and each element
	% P(m+1,d1,d2...) contains the associated Legendre function of degree n and order
	% m evaluated at X(d1,d2...).

    Pn = legendre(n,cos(theta'))';

    posi_Y = [];
    nega_Y = [];

    m= 0:n;
    v = sqrt(((2*n+1)/(4*pi))*(fact(n-m+1)./fact(n+m+1)));
    v = v(ones(1,vnum),:).*Pn(:,m+1).*exp_i_m_phi(:,m+1);
    posi_Y(:,m+1) = v; % positive order;
    nega_Y(:,n-m+1) = sign_m(ones(1,vnum),m+1).*conj(v); % negative order

    Z(:,end+1:end+n) = nega_Y(:,1:n);
    Z(:,end+1:end+n+1) = posi_Y(:,1:(n+1));
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function  [fig_handle,tot_image,image_XYYZ,image_YX] = DrawCellsOnAnat(I,ax)
if exist('ax','var') && isobject(ax)
       axes(ax);
       fig_handle = gcf;
else
    fig_handle = figure('Position',[600,50,458.5,608],'color',[1 1 1],...
        'Name',['Fish#' num2str(I.i_fish)]);
    set(gca,'Position',[0.01, 0.01, 0.98, 0.98]); % right ~subplot
end

%% load from Input struct that is generated by 'LoadCurrentFishForAnatPlot.m'
names = fieldnames(I); % cell of strings
for i = 1:length(names),
    eval([names{i} ' = I.',names{i},';']);
end

if ~isRefAnat, % raw images
    k_zres_ratio = 20; % 8/0.406 = 19.704
    radius_xy = 5;%7;
    width_z = 10;
    thickness_z = 1;
else % registered to ZBrain
    k_zres_ratio = 2.5; % 2/0.798 = 2.506
    radius_xy = 3;
    width_z = 5;
    thickness_z = 3;
    z_range_ventral = 40;% round(min(CellXYZ(:,3))*k_zres_ratio);
end

%% Setup
[s1,s2] = size(cIX);
if s2>s1,
    cIX = cIX';
end
[s1,s2] = size(gIX);
if s2>s1,
    gIX = gIX';
end

% down-sample
if ~isPopout,
    displaymax = 8000;
    if length(cIX) > displaymax,
        skip = round(length(cIX)/displaymax);
        cIX = cIX(1:skip:end,:);
        gIX = gIX(1:skip:end,:);
    end
end

% set up anatomical images
% anat_YX = zeros(size(anat_yx));
% anat_YZ = zeros(size(anat_yz));
% anat_ZX = zeros(size(anat_zx));
anat_YX = anat_yx/4;
anat_YZ = anat_yz/4;
anat_ZX = anat_zx/4;
anat_ZX = flipud(anat_ZX);
dimv_yx = size(anat_YX);
dimv_yz = size(anat_YZ);
dimv_zx = size(anat_ZX);

% adjust size for YZ and ZX view, and combine all 3 view into total-image
% (including white borders as indicated in dim_totimage)

anat_yz2=zeros(dimv_yz(1),dimv_yz(2)*k_zres_ratio,3);
anat_zx2=zeros(dimv_zx(1)*k_zres_ratio,dimv_zx(2),3);

% Show cropped version for normalized fish
% (only for plotting for now, because drawing on GUI relies on the true
% coordinates)

if isRefAnat && isPopout,
    % crop lengthwise (most of tail), and scale k_zres
    y_range = 81:980;%81:990;% 81:1104; % max 1406
    z_range = z_range_ventral:345; % max 138*k_zres
    dimv_yx3 = size(anat_YX(y_range,:,:));
    dimv_yz3 = [length(y_range),length(z_range),3];
    dimv_zx3 = [length(z_range),size(anat_YX,2),3];

    dim_totimage = [dimv_yx3(1)+dimv_zx3(1)+10,dimv_yx3(2)+dimv_yz3(2)+10,3];
else
    dim_totimage = [dimv_yx(1)+dimv_zx(1)*k_zres_ratio+10,dimv_yx(2)+dimv_yz(2)*k_zres_ratio+10,3];
end
tot_image=ones(dim_totimage);

% create masks (~linearized indices), to draw cells on image (later)
circlemaskIX = MakeCircularMask(radius_xy,dimv_yx);
yzmaskIX = MakeSquareMask(width_z,thickness_z,dimv_yz);
zxmaskIX = MakeSquareMask(thickness_z,width_z,dimv_zx);

%% Draw each cell on map (all 3 views)
% main loop
cIX_abs = absIX(cIX);
M_xyz = CellXYZ(cIX_abs,:);

if ~isempty(cIX)
    anat_YX = DrawMasksInRGB(anat_YX,M_xyz(:,[1,2]),circlemaskIX,clrmap,gIX,clr_alpha);
    anat_YZ = DrawMasksInRGB(anat_YZ,M_xyz(:,[1,3]),yzmaskIX,clrmap,gIX,clr_alpha/2);
    anat_ZX = DrawMasksInRGB(anat_ZX,M_xyz(:,[3,2]),zxmaskIX,clrmap,gIX,clr_alpha/2);
end

%% Draw masks on anat projections
if isShowMasks,
    nMasks = length(Msk_IDs);
    % set params
    rng(1);
    %     cmap2 = rand(nMasks,3);
    cmap2 = hsv(nMasks);%jet(nMasks);
    outline_radius = 3;
    msk_alpha = 0.1;%0.15 %0.25 for all projections
    white_alpha = 0.15;%0.25 % 0 for all projections
    for i = 1:nMasks,
        clr = cmap2(i,:);
        % get masks from database
        if Msk_IDs(i)==0,
            msk = full(newMask);
        else
            msk = full(MASKs.MaskDatabase(:,Msk_IDs(i)));
        end
        mask_3D = reshape(msk, [MASKs.height, MASKs.width, MASKs.Zs]);

        % Y-X
        masks_XY = max(mask_3D,[],3);
        masks_fit = logical(imresize(masks_XY,[size(anat_yx,1),size(anat_yx,2)]));
        if isShowMskOutline,
            masks_fit = imdilate(edge(masks_fit), strel('disk',outline_radius));
        end
        anat_YX = DrawMasksInRGB(anat_YX,[],masks_fit,clr,1,msk_alpha,white_alpha);

        % Y-Z
        masks_YZ = squeeze(max(mask_3D,[],2));
        masks_fit = logical(imresize(masks_YZ,[size(anat_yz,1),size(anat_yz,2)]));
        if isShowMskOutline,
            masks_fit = imdilate(edge(masks_fit), ones(5,2));%strel('line',outline_radius,90));
        end
        anat_YZ = DrawMasksInRGB(anat_YZ,[],masks_fit,clr,1,msk_alpha,white_alpha);

        % Z-X
        masks_ZX = squeeze(max(mask_3D,[],1))';  % notice the transpose
        masks_fit = logical(imresize(masks_ZX,[size(anat_zx,1),size(anat_zx,2)]));
        if isShowMskOutline,
            masks_fit = imdilate(edge(masks_fit), ones(2,5));%strel('line',outline_radius,0));
        end
        anat_ZX = DrawMasksInRGB(anat_ZX,[],masks_fit,clr,1,msk_alpha,white_alpha);
    end
end

%% rescale (low-res) z dimension
for k=1:3
    anat_yz2(:,:,k) = imresize(anat_YZ(:,:,k), [dimv_yz(1), dimv_yz(2)*k_zres_ratio],'nearest');
    anat_zx2(:,:,k) = imresize(anat_ZX(:,:,k), [dimv_zx(1)*k_zres_ratio, dimv_zx(2)],'nearest');
end

% Draw scale bar
if isRefAnat
    if ~isPopout,
        y_range = 1:size(anat_YX,1);
    end
    k_um = 0.798; % ZBrain: x/y/z = 0.798/0.798/2um
    scalebar_len = round(50/k_um);
    scalebar_x_range = (600-scalebar_len+1):600;
    scalebar_y = y_range(end)-30-y_range(1)+1;
    scalebar_y_range = (scalebar_y-3):(scalebar_y+3);

    image_YX = anat_YX(y_range,:,:);
    image_YX(scalebar_y_range,scalebar_x_range,:) = 1;
else
    k_um = 0.406; % xy_res
    scalebar_len = round(50/k_um);

    x_end = size(anat_YX,2)-30;
    scalebar_x_range = (x_end-scalebar_len+1):x_end;
    scalebar_y = size(anat_YX,1)-30+1;
    scalebar_y_range = (scalebar_y-3):(scalebar_y+3);

    image_YX = anat_YX;
    image_YX(scalebar_y_range,scalebar_x_range,:) = 1;
end

% compile 3 projections to image panel
if isRefAnat && isPopout,
    tot_image(dimv_zx3(1)+11:end,1:dimv_yx3(2),:) = image_YX;%anat_YX(y_range,:,:);
    tot_image(dimv_zx3(1)+11:end,dimv_yx3(2)+11:end,:) = anat_yz2(y_range,z_range,:);
    tot_image(1:dimv_zx3(1),1:dimv_zx3(2),:) = flipud(anat_zx2(z_range,:,:));
    image_XYYZ = tot_image(dimv_zx3(1)+11:end,1:end,:);
else
    tot_image(dimv_zx(1)*k_zres_ratio+11:end,1:dimv_yx(2),:) = image_YX;
    tot_image(dimv_zx(1)*k_zres_ratio+11:end,dimv_yx(2)+11:end,:) = anat_yz2;
    tot_image(1:dimv_zx(1)*k_zres_ratio,1:dimv_zx(2),:) = flipud(anat_zx2);
    image_XYYZ = tot_image(dimv_zx(1)*k_zres_ratio+11:end,1:end,:);
end

% OPTIONAL BRIGHTENING
% tot_image = tot_image*2;

tot_image(tot_image(:)>1) = 1;
tot_image(tot_image(:)<0) = 0;

image(image_XYYZ); % image(tot_image);
axis image;axis off

set(gcf,'Color','w');

end

author: xiuyechen

repo: FishExplorer


text:
function wss_dist= comp_wss(cleanFile, enhancedFile);
% ----------------------------------------------------------------------
%
%     Weighted Spectral Slope (WSS) Objective Speech Quality Measure
%
%     This function implements the Weighted Spectral Slope (WSS)
%     distance measure originally proposed in [1].  The algorithm
%     works by first decomposing the speech signal into a set of
%     frequency bands (this is done for both the test and reference
%     frame).  The intensities within each critical band are
%     measured.  Then, a weighted distances between the measured
%     slopes of the log-critical band spectra are computed.
%     This measure is also described in Section 2.2.9 (pages 56-58)
%     of [2].
%
%     Whereas Klatt's original measure used 36 critical-band
%     filters to estimate the smoothed short-time spectrum, this
%     implementation considers a bank of 25 filters spanning
%     the 4 kHz bandwidth.
%
%   Usage:  wss_dist=comp_wss(cleanFile.wav, enhancedFile.wav)
%
%         cleanFile.wav - clean input file in .wav format
%         enhancedFile  - enhanced output file in .wav format
%         wss_dist      - computed spectral slope distance
%
%  Example call:  ws =comp_wss('sp04.wav','enhanced.wav')
%
%  References:
%
%     [1] D. H. Klatt, "Prediction of Perceived Phonetic Distance
%	    from Critical-Band Spectra: A First Step", Proc. IEEE
%	    ICASSP'82, Volume 2, pp. 1278-1281, May, 1982.
%
%     [2] S. R. Quackenbush, T. P. Barnwell, and M. A. Clements,
%	    Objective Measures of Speech Quality.  Prentice Hall
%	    Advanced Reference Series, Englewood Cliffs, NJ, 1988,
%	    ISBN: 0-13-629056-6.
%
%  Authors: Bryan L. Pellom and John H. L. Hansen (July 1998)
%  Modified by: Philipos C. Loizou  (Oct 2006)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%
% ----------------------------------------------------------------------
if nargin~=2
    fprintf('USAGE: WSS=comp_wss(cleanFile.wav, enhancedFile.wav)\n');
    fprintf('For more help, type: help comp_wss\n\n');
    return;
end

alpha= 0.95;

[data1, Srate1, Nbits1]= wavread(cleanFile);
[data2, Srate2, Nbits2]= wavread(enhancedFile);
if ( Srate1~= Srate2) | ( Nbits1~= Nbits2)
    error( 'The two files do not match!\n');
end

len= min( length( data1), length( data2));
data1= data1( 1: len)+eps;
data2= data2( 1: len)+eps;

wss_dist_vec= wss( data1, data2,Srate1);
wss_dist_vec= sort( wss_dist_vec);
wss_dist= mean( wss_dist_vec( 1: round( length( wss_dist_vec)*alpha)));

function distortion = wss(clean_speech, processed_speech,sample_rate)

% ----------------------------------------------------------------------
% Check the length of the clean and processed speech.  Must be the same.
% ----------------------------------------------------------------------

clean_length      = length(clean_speech);
processed_length  = length(processed_speech);

if (clean_length ~= processed_length)
  disp('Error: Files  musthave same length.');
  return
end

% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------

winlength   = round(30*sample_rate/1000); 	   % window length in samples
skiprate    = floor(winlength/4);		   % window skip in samples
max_freq    = sample_rate/2;	   % maximum bandwidth
num_crit    = 25;		   % number of critical bands

USE_FFT_SPECTRUM = 1;		   % defaults to 10th order LP spectrum
n_fft       = 2^nextpow2(2*winlength);
n_fftby2    = n_fft/2;		   % FFT size/2
Kmax        = 20;		   % value suggested by Klatt, pg 1280
Klocmax     = 1;		   % value suggested by Klatt, pg 1280

% ----------------------------------------------------------------------
% Critical Band Filter Definitions (Center Frequency and Bandwidths in Hz)
% ----------------------------------------------------------------------

cent_freq(1)  = 50.0000;   bandwidth(1)  = 70.0000;
cent_freq(2)  = 120.000;   bandwidth(2)  = 70.0000;
cent_freq(3)  = 190.000;   bandwidth(3)  = 70.0000;
cent_freq(4)  = 260.000;   bandwidth(4)  = 70.0000;
cent_freq(5)  = 330.000;   bandwidth(5)  = 70.0000;
cent_freq(6)  = 400.000;   bandwidth(6)  = 70.0000;
cent_freq(7)  = 470.000;   bandwidth(7)  = 70.0000;
cent_freq(8)  = 540.000;   bandwidth(8)  = 77.3724;
cent_freq(9)  = 617.372;   bandwidth(9)  = 86.0056;
cent_freq(10) = 703.378;   bandwidth(10) = 95.3398;
cent_freq(11) = 798.717;   bandwidth(11) = 105.411;
cent_freq(12) = 904.128;   bandwidth(12) = 116.256;
cent_freq(13) = 1020.38;   bandwidth(13) = 127.914;
cent_freq(14) = 1148.30;   bandwidth(14) = 140.423;
cent_freq(15) = 1288.72;   bandwidth(15) = 153.823;
cent_freq(16) = 1442.54;   bandwidth(16) = 168.154;
cent_freq(17) = 1610.70;   bandwidth(17) = 183.457;
cent_freq(18) = 1794.16;   bandwidth(18) = 199.776;
cent_freq(19) = 1993.93;   bandwidth(19) = 217.153;
cent_freq(20) = 2211.08;   bandwidth(20) = 235.631;
cent_freq(21) = 2446.71;   bandwidth(21) = 255.255;
cent_freq(22) = 2701.97;   bandwidth(22) = 276.072;
cent_freq(23) = 2978.04;   bandwidth(23) = 298.126;
cent_freq(24) = 3276.17;   bandwidth(24) = 321.465;
cent_freq(25) = 3597.63;   bandwidth(25) = 346.136;

bw_min      = bandwidth (1);	   % minimum critical bandwidth

% ----------------------------------------------------------------------
% Set up the critical band filters.  Note here that Gaussianly shaped
% filters are used.  Also, the sum of the filter weights are equivalent
% for each critical band filter.  Filter less than -30 dB and set to
% zero.
% ----------------------------------------------------------------------

min_factor = exp (-30.0 / (2.0 * 2.303));       % -30 dB point of filter

for i = 1:num_crit
  f0 = (cent_freq (i) / max_freq) * (n_fftby2);
  all_f0(i) = floor(f0);
  bw = (bandwidth (i) / max_freq) * (n_fftby2);
  norm_factor = log(bw_min) - log(bandwidth(i));
  j = 0:1:n_fftby2-1;
  crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
  crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
end

% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Weighted Spectral
% Slope Measure
% ----------------------------------------------------------------------

num_frames = clean_length/skiprate-(winlength/skiprate); % number of frames
start      = 1;					% starting sample
window     = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));

for frame_count = 1:num_frames

   % ----------------------------------------------------------
   % (1) Get the Frames for the test and reference speech.
   %     Multiply by Hanning Window.
   % ----------------------------------------------------------

   clean_frame = clean_speech(start:start+winlength-1);
   processed_frame = processed_speech(start:start+winlength-1);
   clean_frame = clean_frame.*window;
   processed_frame = processed_frame.*window;

   % ----------------------------------------------------------
   % (2) Compute the Power Spectrum of Clean and Processed
   % ----------------------------------------------------------

    if (USE_FFT_SPECTRUM)
       clean_spec     = (abs(fft(clean_frame,n_fft)).^2);
       processed_spec = (abs(fft(processed_frame,n_fft)).^2);
    else
       a_vec = zeros(1,n_fft);
       a_vec(1:11) = lpc(clean_frame,10);
       clean_spec     = 1.0/(abs(fft(a_vec,n_fft)).^2)';

       a_vec = zeros(1,n_fft);
       a_vec(1:11) = lpc(processed_frame,10);
       processed_spec = 1.0/(abs(fft(a_vec,n_fft)).^2)';
    end

   % ----------------------------------------------------------
   % (3) Compute Filterbank Output Energies (in dB scale)
   % ----------------------------------------------------------

   for i = 1:num_crit
      clean_energy(i) = sum(clean_spec(1:n_fftby2) ...
		            .*crit_filter(i,:)');
      processed_energy(i) = sum(processed_spec(1:n_fftby2) ...
			        .*crit_filter(i,:)');
   end
   clean_energy = 10*log10(max(clean_energy,1E-10));
   processed_energy = 10*log10(max(processed_energy,1E-10));

   % ----------------------------------------------------------
   % (4) Compute Spectral Slope (dB[i+1]-dB[i])
   % ----------------------------------------------------------

   clean_slope     = clean_energy(2:num_crit) - ...
		     clean_energy(1:num_crit-1);
   processed_slope = processed_energy(2:num_crit) - ...
		     processed_energy(1:num_crit-1);

   % ----------------------------------------------------------
   % (5) Find the nearest peak locations in the spectra to
   %     each critical band.  If the slope is negative, we
   %     search to the left.  If positive, we search to the
   %     right.
   % ----------------------------------------------------------

   for i = 1:num_crit-1

       % find the peaks in the clean speech signal

       if (clean_slope(i)>0) 		% search to the right
	  n = i;
          while ((n<num_crit) & (clean_slope(n) > 0))
	     n = n+1;
 	  end
	  clean_loc_peak(i) = clean_energy(n-1);
       else				% search to the left
          n = i;
	  while ((n>0) & (clean_slope(n) <= 0))
	     n = n-1;
 	  end
	  clean_loc_peak(i) = clean_energy(n+1);
       end

       % find the peaks in the processed speech signal

       if (processed_slope(i)>0) 	% search to the right
	  n = i;
          while ((n<num_crit) & (processed_slope(n) > 0))
	     n = n+1;
	  end
	  processed_loc_peak(i) = processed_energy(n-1);
       else				% search to the left
          n = i;
	  while ((n>0) & (processed_slope(n) <= 0))
	     n = n-1;
 	  end
	  processed_loc_peak(i) = processed_energy(n+1);
       end

   end

   % ----------------------------------------------------------
   %  (6) Compute the WSS Measure for this frame.  This
   %      includes determination of the weighting function.
   % ----------------------------------------------------------

   dBMax_clean       = max(clean_energy);
   dBMax_processed   = max(processed_energy);

   % The weights are calculated by averaging individual
   % weighting factors from the clean and processed frame.
   % These weights W_clean and W_processed should range
   % from 0 to 1 and place more emphasis on spectral
   % peaks and less emphasis on slope differences in spectral
   % valleys.  This procedure is described on page 1280 of
   % Klatt's 1982 ICASSP paper.

   Wmax_clean        = Kmax ./ (Kmax + dBMax_clean - ...
		 	    clean_energy(1:num_crit-1));
   Wlocmax_clean     = Klocmax ./ ( Klocmax + clean_loc_peak - ...
				clean_energy(1:num_crit-1));
   W_clean           = Wmax_clean .* Wlocmax_clean;

   Wmax_processed    = Kmax ./ (Kmax + dBMax_processed - ...
			        processed_energy(1:num_crit-1));
   Wlocmax_processed = Klocmax ./ ( Klocmax + processed_loc_peak - ...
			            processed_energy(1:num_crit-1));
   W_processed       = Wmax_processed .* Wlocmax_processed;

   W = (W_clean + W_processed)./2.0;

   distortion(frame_count) = sum(W.*(clean_slope(1:num_crit-1) - ...
		       processed_slope(1:num_crit-1)).^2);

   % this normalization is not part of Klatt's paper, but helps
   % to normalize the measure.  Here we scale the measure by the
   % sum of the weights.

   distortion(frame_count) = distortion(frame_count)/sum(W);

   start = start + skiprate;

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [ icc, ccc ] = st_to_cc_index ( nst, ist, jst, ncc, n )

%*****************************************************************************80
%
%% ST_TO_CC_INDEX creates CC indices from ST data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NST, the number of ST elements.
%
%    Input, integer IST(NST), JST(NST), the ST rows and columns.
%
%    Input, integer NCC, the number of CC elements.
%
%    Input, integer N, the number of columns.
%
%    Output, integer ICC(NCC), the CC rows.
%
%    Output, integer CCC(N+1), the compressed CC columns.
%

%
%  Sort the elements.
%
  [ jst2, ist2 ] = i4vec2_sort_a ( nst, jst, ist );
%
%  Get the unique elements.
%
  [ jcc, icc ] = i4vec2_sorted_uniquely ( nst, jst2, ist2, ncc );
%
%  Compress the column index.
%
  ccc = zeros(n+1,1);

  ccc(1) = 1;
  jlo = 1;
  for i = 1 : ncc
    jhi = jcc(i);
    if ( jhi ~= jlo )
      ccc(jlo+1:jhi) = i;
      jlo = jhi;
    end
  end
  jhi = n + 1;
  ccc(jlo+1:jhi) = ncc + 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [W, U, mu] = get_svds(dWU, Nrank)

[Wall, Sv, Uall] = svd(gather_try(dWU), 0);
[~, imax] = max(abs(Wall(:,1)));
Uall(:,1) = -Uall(:,1) * sign(Wall(imax,1));
Wall(:,1) = -Wall(:,1) * sign(Wall(imax,1));

%     [~, imin] = min(diff(Wall(:,1), 1));
%     [~, imin] = min(Wall(:,1));
%     dmax(k) = - (imin- 20);

%     if dmax(k)>0
%         dWU((dmax(k) + 1):nt0, :,k) = dWU(1:nt0-dmax(k),:, k);
%         Wall((dmax(k) + 1):nt0, :)  = Wall(1:nt0-dmax(k),:);
%     else
%         dWU(1:nt0+dmax(k),:, k) = dWU((1-dmax(k)):nt0,:, k);
%         Wall(1:nt0+dmax(k),:) = Wall((1-dmax(k)):nt0,:);
%     end

Wall = Wall * Sv;

Sv = diag(Sv);
mu = sum(Sv(1:Nrank).^2).^.5;
Wall = Wall/mu;

W = Wall(:,1:Nrank);
U = Uall(:,1:Nrank);

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort


text:
function nvar = p00_nvar ( problem, option )

%*****************************************************************************80
%
%% P00_NVAR sets the number of variables for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer OPTION, the option index.
%
%    Output, integer NVAR, the number of variables.
%
  if (  problem == 1 )
    nvar = p01_nvar ( option );
  elseif ( problem == 2 )
    nvar = p02_nvar ( option );
  elseif ( problem == 3 )
    nvar = p03_nvar ( option );
  elseif ( problem == 4 )
    nvar = p04_nvar ( option );
  elseif ( problem == 5 )
    nvar = p05_nvar ( option );
  elseif ( problem == 6 )
    nvar = p06_nvar ( option );
  elseif ( problem == 7 )
    nvar = p07_nvar ( option );
  elseif ( problem == 8 )
    nvar = p08_nvar ( option );
  elseif ( problem == 9 )
    nvar = p09_nvar ( option );
  elseif ( problem == 10 )
    nvar = p10_nvar ( option );
  elseif ( problem == 11 )
    nvar = p11_nvar ( option );
  elseif ( problem == 12 )
    nvar = p12_nvar ( option );
  elseif ( problem == 13 )
    nvar = p13_nvar ( option );
  elseif ( problem == 14 )
    nvar = p14_nvar ( option );
  elseif ( problem == 15 )
    nvar = p15_nvar ( option );
  elseif ( problem == 16 )
    nvar = p16_nvar ( option );
  elseif ( problem == 17 )
    nvar = p17_nvar ( option );
  elseif ( problem == 18 )
    nvar = p18_nvar ( option );
  elseif ( problem == 19 )
    nvar = p19_nvar ( option );
  elseif ( problem == 20 )
    nvar = p20_nvar ( option );
  else
    printf ( 1, '\n' );
    printf ( 1, 'P00_NVAR - Fatal error!\n' );
    printf ( 1, '  Unrecognized problem number = %d\n', problem );
    error ( 'P00_NVAR - Fatal error!' );
  end

  return
end


text:
function b = r8cbb_to_r8ge ( n1, n2, ml, mu, a )

%*****************************************************************************80
%
%% R8CBB_TO_R8GE copies a R8CBB matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8CBB storage format is for a compressed border banded matrix.
%    Such a matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    The R8CBB format is the same as the DBB format, except that the banded
%    matrix A1 is stored in compressed band form rather than standard
%    banded form.  In other words, we do not include the extra room
%    set aside for fill in during pivoting.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (ML+MU+1)*N1 entries of A, using the obvious variant
%    of the LINPACK general band format.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+MU+1)+(J-1)*(ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= n1.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks.
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N1-1.
%
%    Input, real A((ML+MU+1)*N1+2*N1*N2+N2*N2), the R8CBB matrix.
%
%    Output, real B(N1+N2,N1+N2), the R8GE matrix.
%
  for i = 1 : n1
    for j = 1 : n1

      if ( mu+ml < (j-i) | ml < (i-j) )
        b(i,j) = 0.0;
      else
        ij = (i-j+mu+1)+(j-1)*(ml+mu+1);
        b(i,j) = a(ij);
      end

    end
  end

  for i = 1 : n1
    for j = n1+1 : n2
      ij = (ml+mu+1)*n1+(j-n1-1)*n1+i;
      b(i,j) = a(ij);
    end
  end

  for i = n1+1 : n2
    for j = 1 : n1+n2
      ij = (ml+mu+1)*n1+n2*n1+(j-1)*n2+(i-n1);
      b(i,j) = a(ij);
    end
  end

  return
end