text:
function params = stblfit(X,varargin)
%PARAMS = STBLFIT(X) returns an estimate of the four parameters in a
% fit of the alpha-stable distribution to the data X.  The output
% PARAMS is a 4 by 1 vector which contains the estimates of the
% characteristic exponent ALPHA, the skewness BETA, the scale GAMMA and
% location DELTA.
%
%PARAMS = STBLFIT(X,METHOD) Specifies the algorithm used to
% estimate the parameters.  The choices for METHOD are
%        'ecf' - Fits the four parameters to the empirical characteristic
%               function estimated from the data.  This is the default.
%               Based on Koutrouvelis (1980,1981), see [1],[2] below.
% 'percentile' - Fits the four parameters using various
%                percentiles of the data X.  This is faster than 'ecf',
%                however studies have shown it to be slightly less
%                accurate in general.
%                Based on McCulloch (1986), see [2] below.
%
%PARAMS = STBLFIT(...,OPTIONS) specifies options used in STBLFIT.  OPTIONS
% must be an options stucture created with the STATSET function.  Possible
% options are
%       'Display' - When set to 'iter', will display the values of
%                   alpha,beta,gamma and delta in each
%                   iteration.  Default is 'off'.
%       'MaxIter' - Specifies the maximum number of iterations allowed in
%                   estimation. Default is 5.
%       'TolX'    - Specifies threshold to stop iterations. Default is
%                   0.01.
%
%   See also: STBLRND, STBLPDF, STBLCDF, STBLINV
%
% References:
%   [1] I. A. Koutrouvelis (1980)
%       "Regression-Type Estimation of the Paramters of Stable Laws.
%       JASA, Vol 75, No. 372
%
%   [2] I. A. Koutrouvelis (1981)
%       "An Iterative Procedure for the estimation of the Parameters of
%       Stable Laws"
%       Commun. Stat. - Simul. Comput. 10(1), pages 17-28
%
%   [3] J. H. McCulloch (1986)
%       "Simple Consistent Estimators of Stable Distribution Parameters"
%       Cummun. Stat. Simul. Comput. 15(4)
%
%   [4] A. H. Welsh (1986)
%       "Implementing Empirical Characteristic Function Procedures"
%       Statistics & Probability Letters Vol 4, pages 65-67

% ==== Gather additional options
dispit = false;
maxiter = 5;
tol = .01;
if ~isempty(varargin)
    if isstruct(varargin{end})
        opt = varargin{end};
        try
            dispit = opt.Display;
        catch ME
            error('OPTIONS must be a structure created with STATSET');
        end
        if ~isempty(opt.MaxIter)
            maxiter = opt.MaxIter;
        end
        if ~isempty(opt.TolX)
            tol = opt.TolX;
        end
    end
end

if strcmp(dispit,'iter')
    dispit = true;
    fprintf('    iteration\t    alpha\t    beta\t  gamma\t\t  delta\n');
    dispfmt = '%8d\t%14g\t%8g\t%8g\t%8g\n';
end

% === Find which method.
if any(strcmp(varargin,'percentile'))
    maxiter = 0; % This is McCulloch's percentile method
end

% ==== Begin estimation =====
N = numel(X); % data size

% function handle to compute empirical char. functions
I = sqrt(-1);
phi = @(theta,data) 1/numel(data) * sum( exp( I * ...
    reshape(theta,numel(theta),1) *...
    reshape(data,1,numel(data)) ) , 2);

% Step 1 - Obtain initial estimates of parameters using McCulloch's method
%          then standardize data
[alpha beta] = intAlpBet(X);
[gam delta ] = intGamDel(X,alpha,beta);

if gam==0
    % Use standard deviation as initial guess
    gam = std(X);
end

s = (X - delta)/gam;

if dispit
    fprintf(dispfmt,0,alpha,beta,gam,delta);
end

% Step 2 - Iterate until convergence
alphaold = alpha;
deltaold = delta;
diffbest = inf;
for iter = 1:maxiter

    % Step 2.1 - Regress against ecf to refine estimates of alpha & gam
    %            After iteration 1, use generalized least squares
    if iter <= 2
        K = chooseK(alpha,N);
        t = (1:K)*pi/25;
        w = log(abs(t));
    end

    y = log( - log( abs(phi(t,s)).^2 ) );

    if iter == 1  % use ordinary least squares regression
        ell = regress(y,[w' ones(size(y))]);
        alpha = ell(1);
        gamhat = (exp(ell(2))/2)^(1/alpha);
        gam = gam * gamhat;
    else          % use weighted least squares regression
        sig = charCov1(t ,N, alpha , beta, 1);
        try
            ell = lscov([w' ones(size(y))],y,sig);
        catch   % In case of badly conditioned covariance matrix, just use diagonal entries
            try
                ell = lscov([w' ones(size(y))],y,eye(K).*(sig+eps));
            catch
                break
            end
        end
        alpha = ell(1);
        gamhat = (exp(ell(2))/2)^(1/alpha);
        gam = gam * gamhat;
    end

    % Step 2.2 - Rescale data by estimated scale, truncate
    s = s/gamhat;
    alpha = max(alpha,0);
    alpha = min(alpha,2);
    beta = min(beta,1);
    beta = max(beta,-1);
    gam = max(gam,0);

    % Step 2.3 - Regress against ecf to refine estimates of beta, delta
    %            After iteration 1, use generalized least squares
    if iter <=  2
        L = chooseL(alpha,N);
        % To ensure g is continuous, find first zero in real part of ecf
        A = efcRoot(s);
        u = (1:L)*min(pi/50,A/L);
    end

    ecf = phi(u,s);
    U = real(ecf);
    V = imag(ecf);
    g = atan2(V,U);
    if iter == 1  % use ordinary least squares
        ell = regress(g, [u',  sign(u').*abs(u').^alpha]);
        beta =  ell(2)/tan(alpha*pi/2) ;
        delta = delta + gam* ell(1) ;
    else         % use weighted least squares regression
        sig = charCov2(u ,N, alpha , beta, 1);
        try
            ell = lscov([u',  sign(u').*abs(u').^alpha],g,sig);
        catch % In case of badly conditioned covariance matrix, use diagonal entries
            try
                ell = lscov([u',  sign(u').*abs(u').^alpha],g,eye(L).*(sig+eps));
            catch
                break
            end
        end
        beta =  ell(2)/tan(alpha*pi/2) ;
        delta = delta + gam* ell(1) ;
    end

    % Step 2.4 Remove estimated shift
    s = s - ell(1);

    % display
    if dispit
        fprintf(dispfmt,iter,alpha,beta,gam,delta);
    end

    % Check for blow-up
    if any(isnan([alpha, beta, gam, delta]) | isinf([alpha, beta, gam, delta]))
        break
    end

    % Step 2.5 Check for convergence, keep track of parameters with
    % smallest 'diff'
    diff = (alpha - alphaold)^2 + (delta - deltaold)^2;
    if abs(diff) < diffbest
        bestparams = [alpha; beta; gam; delta];
        diffbest = diff;
        if diff < tol
            break;
        end
    end

    alphaold = alpha;
    deltaold = delta;

end

% Pick best
if maxiter > 0 && iter >= 1
    alpha = bestparams(1);
    beta = bestparams(2);
    gam = bestparams(3);
    delta = bestparams(4);
end

% Step 3 - Truncate if necessary
alpha = max(alpha,0);
alpha = min(alpha,2);
beta = min(beta,1);
beta = max(beta,-1);
gam = max(gam,0);

params = [alpha; beta; gam; delta];

end % End stblfit

%===============================================================
%===============================================================

function [alpha beta] = intAlpBet(X)
% Interpolates Tables found in MuCulloch (1986) to obtain a starting
% estimate of alpha and beta based on percentiles of data X

% Input tables
nuA = [2.439 2.5 2.6 2.7 2.8 3.0 3.2 3.5 4.0 5.0 6.0 8.0 10 15 25];
nuB = [0 .1 .2 .3 .5 .7 1];
[a b] = meshgrid( nuA , nuB );
alphaTab=  [2.000 2.000 2.000 2.000 2.000 2.000 2.000;...
             1.916 1.924 1.924 1.924 1.924 1.924 1.924;...
             1.808 1.813 1.829 1.829 1.829 1.829 1.829;...
             1.729 1.730 1.737 1.745 1.745 1.745 1.745;...
             1.664 1.663 1.663 1.668 1.676 1.676 1.676;...
             1.563 1.560 1.553 1.548 1.547 1.547 1.547;...
             1.484 1.480 1.471 1.460 1.448 1.438 1.438;...
             1.391 1.386 1.378 1.364 1.337 1.318 1.318;...
             1.279 1.273 1.266 1.250 1.210 1.184 1.150;...
             1.128 1.121 1.114 1.101 1.067 1.027 0.973;...
             1.029 1.021 1.014 1.004 0.974 0.935 0.874;...
             0.896 0.892 0.887 0.883 0.855 0.823 0.769;...
             0.818 0.812 0.806 0.801 0.780 0.756 0.691;...
             0.698 0.695 0.692 0.689 0.676 0.656 0.595;...
             0.593 0.590 0.588 0.586 0.579 0.563 0.513]';
betaTab=  [ 0.000 2.160 1.000 1.000 1.000 1.000 1.000;...
            0.000 1.592 3.390 1.000 1.000 1.000 1.000;...
            0.000 0.759 1.800 1.000 1.000 1.000 1.000;...
            0.000 0.482 1.048 1.694 1.000 1.000 1.000;...
            0.000 0.360 0.760 1.232 2.229 1.000 1.000;...
            0.000 0.253 0.518 0.823 1.575 1.000 1.000;...
            0.000 0.203 0.410 0.632 1.244 1.906 1.000;...
            0.000 0.165 0.332 0.499 0.943 1.560 1.000;...
            0.000 0.136 0.271 0.404 0.689 1.230 2.195;...
            0.000 0.109 0.216 0.323 0.539 0.827 1.917;...
            0.000 0.096 0.190 0.284 0.472 0.693 1.759;...
            0.000 0.082 0.163 0.243 0.412 0.601 1.596;...
            0.000 0.074 0.147 0.220 0.377 0.546 1.482;...
            0.000 0.064 0.128 0.191 0.330 0.478 1.362;...
            0.000 0.056 0.112 0.167 0.285 0.428 1.274]';

% Calculate percentiles
Xpcts = prctile(X,[95 75 50 25 5]);
nuAlpha = (Xpcts(1) - Xpcts(5))/(Xpcts(2) - Xpcts(4));
nuBeta = (Xpcts(1) + Xpcts(5) - 2*Xpcts(3))/(Xpcts(1) - Xpcts(5));
% Bring into range
if nuAlpha < 2.4390
    nuAlpha = 2.439 + 1e-12;
elseif nuAlpha > 25
    nuAlpha = 25 - 1e-12;
end

s = sign(nuBeta);

% Get alpha
alpha = interp2(a,b,alphaTab,nuAlpha,abs(nuBeta));

% Get beta
beta = s * interp2(a,b,betaTab,nuAlpha,abs(nuBeta));

% Truncate beta if necessary
if beta>1
    beta = 1;
elseif beta < -1
    beta =-1;
end

end

function [gam delta] = intGamDel(X,alpha,beta)
% Uses McCulloch's Method to obtain scale and location of data X given
% estimates of alpha and beta.

% Get percentiles of data and true percentiles given alpha and beta;
Xpcts = prctile(X,[75 50 25]);

% If alpha is very close to 1, truncate to avoid numerical instability.
warning('off','stblcdf:ScaryAlpha');
warning('off','stblpdf:ScaryAlpha');
if abs(alpha - 1) < .02
    alpha = 1;
end

% With the 'quick' option, these are equivalent to McCulloch's tables
Xquart = stblinv([.75 .25],alpha,beta,1,0,'quick');
Xmed = stblinv(.5,alpha,beta,1,-beta*tan(pi*alpha/2),'quick');

% Obtain gamma as ratio of interquartile ranges
gam = (Xpcts(1) - Xpcts(3))/(Xquart(1) - Xquart(2));

% Obtain delta using median of shifted data and estimate of gamma
zeta = Xpcts(2) - gam * Xmed;
delta = zeta - beta*gam*tan(alpha*pi/2);

end

function K = chooseK(alpha,N)
    % Interpolates Table 1 in [1] to calculate optimum K given alpha and N

    % begin parameters into correct ranges.
    alpha = max(alpha,.3);
    alpha = min(alpha,1.9);
    N = max(N,200);
    N = min(N,1600);
    a = [1.9, 1.5: -.2: .3];
    n = [200 800 1600];
    [X Y] = meshgrid(a,n);
    Kmat = [ 9   9   9 ; ...
            11  11  11 ; ...
            22  16  14 ; ...
            24  18  15 ; ...
            28  22  18 ; ...
            30  24  20 ; ...
            86  68  56 ; ...
            134 124 118  ];
     K = round(interp2(X,Y,Kmat',alpha,N,'linear'));

end

function L = chooseL(alpha,N)
    % Interpolates Table 2 in [1] to calculate optimum L given alpha and N

    alpha = max(alpha,.3);
    alpha = min(alpha,1.9);
    N = max(N,200);
    N = min(N,1600);
    a = [1.9, 1.5, 1.1:-.2:.3];
    n = [200 800 1600];
    [X Y] = meshgrid(a,n);
    Lmat = [ 9  10  11 ; ...
            12  14  15 ; ...
            16  18  17 ; ...
            14  14  14 ; ...
            24  16  16 ; ...
            40  38  36 ; ...
            70  68  66 ];
    L = round(interp2(X,Y,Lmat',alpha,N,'linear'));

end

function A = efcRoot(X)
% An iterative procedure to find the first positive root of the real part
% of the empirical characteristic function of the data X. Based on [4].

N = numel(X);
U = @(theta) 1/N * sum( cos(   ...
    reshape(theta,numel(theta),1) *...
      reshape(X,1,N)       ) , 2 );  % Real part of ecf
m = mean(abs(X));
A = 0;
val = U(A);
iter1 = 0;
while abs(val) > 1e-3 && iter1 < 10^4
    A = A + val/m;
    val = U(A);
    iter1 = iter1 + 1;
end

end

function sig = charCov1(t ,N, alpha , beta,gam)
% Compute covariance matrix of y = log (- log( phi(t) ) ), where phi(t) is
% ecf of alpha-stable random variables. Based on Theorem in [2].

    K = length(t);
    w = tan(alpha*pi/2);
    calpha = gam^alpha;

    Tj = repmat( t(:) , 1 , K);
    Tk = repmat( t(:)'  , K , 1);
    Tjalpha = abs(Tj).^alpha;
    Tkalpha = abs(Tk).^alpha;
    TjxTk = abs(Tj .* Tk);
    TjpTk = Tj + Tk ;
    TjpTkalpha = abs(TjpTk).^alpha;
    TjmTk = Tj - Tk ;
    TjmTkalpha = abs(TjmTk).^alpha;

    A = calpha*( Tjalpha + Tkalpha - TjmTkalpha);
    B = calpha * beta *...
        (-Tjalpha .* sign(Tj) * w ...
        + Tkalpha .* sign(Tk) * w ...
        + TjmTkalpha .* sign(TjmTk) * w) ;
    D = calpha * (Tjalpha + Tkalpha - TjpTkalpha);
    E = calpha * beta *...
        ( Tjalpha .* sign(Tj) * w ...
        + Tkalpha .* sign(Tk) * w ...
        - TjpTkalpha .* sign(TjpTk) * w);

    sig = (exp(A) .* cos(B) + exp(D).*cos(E) - 2)./...
          ( 2 * N * gam^(2*alpha) * TjxTk.^alpha);

end

function sig = charCov2(t ,N, alpha , beta, gam)
% Compute covariance matrix of z = Arctan(imag(phi(t))/real(phi(t)),
% where phi(t) is ecf of alpha-stable random variables.
% Based on Theorem in [2].
    K = length(t);
    w = tan(alpha*pi/2);
    calpha = gam^alpha;

    Tj = repmat( t(:) , 1 , K);
    Tk = repmat( t(:)'  , K , 1);
    Tjalpha = abs(Tj).^alpha;
    Tkalpha = abs(Tk).^alpha;
    TjpTk = Tj + Tk ;
    TjpTkalpha = abs(TjpTk).^alpha;
    TjmTk = Tj - Tk ;
    TjmTkalpha = abs(TjmTk).^alpha;

    B = calpha * beta *...
        (-Tjalpha .* sign(Tj) * w ...
        + Tkalpha .* sign(Tk) * w ...
        + TjmTkalpha .* sign(TjmTk) * w) ;
    E = calpha * beta *...
        ( Tjalpha .* sign(Tj) * w ...
        + Tkalpha .* sign(Tk) * w ...
        - TjpTkalpha .* sign(TjpTk) * w);
    F = calpha * (Tjalpha + Tkalpha);
    G = -calpha * TjmTkalpha;
    H = -calpha * TjpTkalpha;

    sig = exp(F) .*(exp(G) .* cos(B) - exp(H) .* cos(E))/(2*N);

end

author: Sable

repo: mcbench-benchmarks


text:
function [ccg,ic] = chXcorr(hc_L,hc_R,fs,varargin)
%CHXCORR Calculate cross-correlograms with a wide range of options.
%
%   CCG = IOSR.AUDITORY.CHXCORR(HC_L,HC_R,FS) cross-correlates the input
%   2-D matrices HC_L and HC_R over 10ms frame with a maximum lag of 1ms.
%   It is assumed that the number of frequency channels is min(size(HC_L))
%   and hence HC_L and HC_R can be in either orientation. The
%   cross-correlograms consist of cross-correlations for every frame and
%   frequency channel. CCG has dimensions [lag,frequency,frame]. The
%   function calculates running cross-correlations for every sample and
%   integrates these cross-correlations over each frame. The number of
%   frames frame_count is calculated thus:
%
%       frame_count = ...
%           floor((max(size(hc_L))-maxlag-1)/frame_length);
%
%   The underlying cross-correlation algorithm is based on that proposed by
%   Faller & Merimaa [1]. In this implmentation, the time constant of the
%   backward infinite exponential window is given by tau (in samples).
%
%   CCG = IOSR.AUDITORY.CHXCORR(HC_L,HC_R,FS,'PARAMETER',VALUE) allows a
%   number of options to be specified. The options are:
%
%   ({} indicates the default value)
%
%   'frame_length'   : {round(0.01*fs)} | scalar
%       The length of frames used to calculate for integrating
%       cross-correlations.
%   'noverlap'       : {1} | scalar
%       The number of frames over which to integrate the
%       cross-correlations. Note that the frame count is reduced
%       accordingly.
%   'maxlag'         : {round(0.001*fs)} | scalar
%       The maximum lag of the cross-correlation.
%   'tau'            : {round(0.01*fs)} | scalar
%       The time constant of the exponential window used to calculate
%       running cross-correlations. It is recommended that if norm_flag = 1
%       then tau >> 1. As can be seen below, as tau -> 1 then
%       [aL(m,i,n?1), aR(m,i,n?1)] -> 0, and hence c(m,i,n) -> 1.
%   'inhib'          : {[]} | array
%       Specificies an array with which to multiply the cross-correlations
%       before they are integrated. The value defaults to an empty array,
%       meaning that no inhibition will be applied.
%   'ic_t'           : {0} | scalar
%       Specifies the interaural coherence (IC) threshold. Only samples for
%       which the IC exceeds this threshold will be used to integrate
%       cross-correlations. The algorithm calculates Interaural Coherence
%       (IC) according to [1]. The value should be in the range [0,1].
%   'norm_flag'      : {0} | scalar
%       Specifies whether the cross-correlograms are calculated using
%       normalised cross-correlations. A non-zero value indicates that
%       normalised cross-correlations are used.
%   'inhib_mode'     : {'subtract'} | 'multiply'
%       Specify how the inhibition is applied. The default 'subtract' will
%       subtract inhib from the running cross-correlations (negative values
%       are set to zero); 'multiply' will multiply inhib with the running
%       cross-correlations.
%
%   [CCG,IC] = IOSR.AUDITORY.CHXCORR(...) returns the calculated IC to the
%   matrix IC. Although the matrix returned is the same size as hc_L, IC is
%   only calculated for samples 1:frame_count*frame_length, other values
%   will be set to 0.
%
%   Algorithm
%
%   The running normalised cross-correlation is calculated as [1]:
%
%   C(m,i,n) = c(m,i,n) / sqrt( aL(m,i,n) * aR(m,i,n)
%
%   where
%
%   c(m,i,n) = (1/tau) * HC_L(i, max(n+m,n)) * HC_R(i, max(n-m,n)) + ...
%       (1-1/tau) * c(m,i,n-1),
%
%   aL(m,i,n) = (1/tau) * (HC_L(i, max(n+m,n)))^2 + ...
%       (1-1/tau) * aL(m,i,n-1),
%
%   aR(m,i,n) = (1/tau) * (HC_R(i, max(n+m,n)))^2 + ...
%       (1-1/tau) * aR(m,i,n-1),
%
%   i is the frequency index, m is the lag index, and n is the sample
%   index. The running (non-normalised) cross-correlation is calculated is
%   simply c(m,i,n).
%
%   The interaural coherence is
%
%   IC(i,n) = max(C(m,i,n)) % (i.e. over m)
%
%   The cross-correlogram is calculated as the sum of cross- correlations
%   in a given frame:
%
%   CCG(m,i,j) = sum(C(m,i,J),3)
%
%   where
%
%   J = (j-1) * frame_length + 1 : j * frame_length
%
%   References
%
%   [1] C. Faller and J. Merimaa, "Source localization in complex listening
%       situations: Selection of binaural cues based on interaural
%       coherence", The Journal of the Acoustical Society of America, vol.
%       116, pp.3075-3089, Nov. 2004.
%
%   Further Reading
%
%   C. Hummersone, R. Mason, and T. Brookes, "A comparison of computational
%       precedence models for source separation in reverberant
%       environments", The Journal of the Audio Engineering Society, vol.
%       61(7/8), pp.508-520, July 2013.

%   Copyright 2016 University of Surrey.

    assert(nargin>=3, 'iosr:chXcorr:nargin', 'Number of input arguments must be greater than or equal to three.')

    if isparameter(varargin,'inhib_mode') && ~isparameter(varargin,'inhib')
        warning('iosr:instIld:inhibMode','''inhib_mode'' specified, but no inhibition array ''inhib''.');
    end

    % Check source file is compiled
    iosr.general.checkMexCompiled('-largeArrayDims',fullfile(fileparts(mfilename('fullpath')),'chXcorr_c.c'))

    options = struct(...
        'frame_length',round(0.01*fs),...
        'noverlap',1,...
        'maxlag',round(0.001*fs),...
        'tau',round(0.01*fs),...
        'inhib',[],...
        'ic_t',0,...
        'norm_flag',0,...
        'inhib_mode','subtract');

    % read parameter/value inputs
    if nargin > 3 % if parameters are specified
        % read the acceptable names
        optionNames = fieldnames(options);
        % count arguments
        nArgs = length(varargin);
        if round(nArgs/2)~=nArgs/2
           error('iosr:chXcorr:nameValuePair','CHXCORR needs propertyName/propertyValue pairs')
        end
        % overwrite defults
        for pair = reshape(varargin,2,[]) % pair is {propName;propValue}
           IX = strcmpi(pair{1},optionNames); % find match parameter names
           if any(IX)
              % do the overwrite
              options.(optionNames{IX}) = pair{2};
           else
              error('iosr:chXcorr:unknownOption','%s is not a recognized parameter name',pair{1})
           end
        end
    end

    % assign options to variables
    frame_length = options.frame_length;
    noverlap = options.noverlap;
    maxlag = options.maxlag;
    tau = options.tau;
    inhib_mode = options.inhib_mode;
    norm_flag = options.norm_flag;
    ic_t = options.ic_t;
    inhib = options.inhib;

    % check inputs
    assert(all(size(hc_L)==size(hc_R)), 'iosr:chXcorr:invalidInput', '''hc_L'' and ''hc_R'' must be the same size')
    assert(round(frame_length)==frame_length && isscalar(frame_length) && frame_length>0, 'iosr:chXcorr:invalidFrame', ...
        '''frame_length'' must be an integer greater than zero')
    assert(round(noverlap)==noverlap && isscalar(noverlap) && noverlap>0, 'iosr:chXcorr:invalidNoverlap', ...
        '''noverlap'' must be an integer greater than zero')
    assert(round(maxlag)==maxlag && isscalar(maxlag) && maxlag>0, 'iosr:chXcorr:invalidMaxlag', ...
        '''maxlag'' must be an integer greater than zero')
    assert(isscalar(tau) && tau>=1, 'iosr:chXcorr:invalidTau', '''tau'' must be a scalar greater than or equal to one')
    assert(isscalar(norm_flag), 'iosr:chXcorr:invalidNorm', '''norm_flag'' must be a scalar')
    assert(isscalar(ic_t) && ic_t>=0 && ic_t<=1, 'iosr:chXcorr:invalidIct', '''ic_t'' must be a scalar in the range [0,1]')
    assert(ischar(inhib_mode), 'iosr:chXcorr:invalidInhibMode', '''inhib_mode'' must be a char array (string)')

    % Calculate frame count
    frame_count = floor(max(size(hc_L))/(frame_length));
    frame_count = frame_count-noverlap+1;

    % Calculate number of frequency channels
    numchans = min(size(hc_L));
    numsamples = max(size(hc_L));

    % Check orientation of HC and inhib data (i.e. that frequency runs across the rows)
    dims = size(hc_L);
    hc_L = check_input(hc_L,2,numchans);
    hc_R = check_input(hc_R,2,numchans);

    % set a flag if data has been transposed in this way
    if dims(1)~=size(hc_L,1)
        rot = true;
    else
        rot = false;
    end

    % set inhibition mode ID
    switch inhib_mode
        case 'multiply'
            inhib_mode_ID = 1;
            if isempty(inhib)
                inhib = ones(size(hc_L));
            end
        case 'subtract'
            inhib_mode_ID = 2;
            if isempty(inhib)
                inhib = zeros(size(hc_L));
            end
        otherwise
            error('iosr:chXcorr:unknownInhibMode','''inhib_mode'' must be set to ''multiply'' or ''subtract''')
    end

    inhib = check_input(inhib,2,numchans);

    % Append HC and inhibition data with zeros for cross-correlation
    hc_L = [hc_L; zeros(maxlag+1,numchans)];
    hc_R = [hc_R; zeros(maxlag+1,numchans)];
    inhib = [inhib; zeros(maxlag+1,numchans)];

    assert(all(size(inhib)==size(hc_L)), 'iosr:chXcorr:invalidInhib', '''inhib'' must be a matrix the same size as ''hc_L'' or ''hc_R''')

    % Calculate cross-correlograms
    [ccg,ic] = iosr.auditory.chXcorr_c(hc_L,hc_R,frame_count,frame_length,noverlap,maxlag,tau,inhib,ic_t,norm_flag,inhib_mode_ID);

    % Correct orientation of IC data, if data was transposed, and crop to remove appended zeros
    ic = ic(1:numsamples,:);
    if rot
        ic = ic';
    end

end

function output = check_input(input,dim,target)
%CHECK_INPUT check input is correct orientation

    if size(input,dim)~=target
        output = input';
        assert(size(output,dim)==target, 'iosr:chXcorr:invalidInputs', 'Input invalid')
    else
        output = input;
    end

end

function set = isparameter(input,parameter)
%ISPARAMETER check for input parameter

    set = any(strcmpi(input(cellfun(@ischar,input)),parameter));

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
% EX_STOKES_SQUARE_NDL: solve the Stokes problem in the unit square with generalized Nedelec elements of the second family.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_square.txt';

% Type of boundary conditions for each side of the domain
problem_data.drchlt_sides = 1:4;
problem_data.nmnn_sides = [];

% Physical parameters
problem_data.viscosity = @(x, y) ones (size (x));

% Force term
fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
problem_data.f  = @(x, y) cat(1, ...
                reshape (fx (x,y), [1, size(x)]), ...
                reshape (fy (x,y), [1, size(x)]));

% Boundary terms
uxex = @(x, y) (sin(x) .* cos(y));
uyex = @(x, y) (-sin(y) .* cos(x));

problem_data.h = @(x, y, iside) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

% Exact solution, to compute the errors
problem_data.velex = @(x, y) cat(1, ...
                  reshape (uxex (x,y), [1, size(x)]), ...
                  reshape (uyex (x,y), [1, size(x)]));

problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);

problem_data.gradvelex = @test_stokes_square_bc_graduex;

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.element_name = 'ndl';   % Element type for discretization
method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
method_data.nsub         = [10 10];  % Number of subdivisions
method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule

% Penalization parameter for Nitsche's method
factor = 10;
method_data.Cpen = factor*(min(method_data.degree)+1);

% 3) CALL TO THE SOLVER
[geometry, msh, space_v, vel, space_p, press] = ...
                       solve_stokes (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) COMPARISON WITH EXACT SOLUTION
error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex)
[error_h1_v, error_l2_v] = ...
   sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex)

% 4.2) EXPORT TO PARAVIEW
output_file = 'SQUARE_NDL_Deg3_Reg2_Sub10';

fprintf ('The result is saved in the files %s \n and %s \n \n', ...
           [output_file '_vel'], [output_file '_press']);
vtk_pts = {linspace(0, 1, 20), linspace(0, 1, 20)};
sp_to_vtk (press, space_p, geometry, vtk_pts, [output_file '_press'], 'press')
sp_to_vtk (vel,   space_v, geometry, vtk_pts, [output_file '_vel'  ], {'velocity', 'divergence'}, {'value', 'divergence'})

% 4.3) PLOT IN MATLAB
[eu, F] = sp_eval (vel, space_v, geometry, vtk_pts);
[X,  Y] = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));

figure()
subplot (1,2,2)
eu2 = problem_data.velex (X, Y);
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
axis equal
title('Exact solution')
subplot (1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
axis equal
title('Computed solution')

%!demo
%! ex_stokes_square_ndl

%!test
%! problem_data.geo_name = 'geo_square.txt';
%! problem_data.drchlt_sides = 1:4;
%! problem_data.nmnn_sides = [];
%! problem_data.viscosity = @(x, y) ones (size (x));
%! fx = @(x, y) (6 * x + y .* cos(x .* y) + 2 * cos(y) .* sin(x));
%! fy = @(x, y) (x .* cos(x .* y) - 2 * cos(x) .* sin(y));
%! problem_data.f  = @(x, y) cat(1, ...
%!                 reshape (fx (x,y), [1, size(x)]), ...
%!                 reshape (fy (x,y), [1, size(x)]));
%! uxex = @(x, y) (sin(x) .* cos(y));
%! uyex = @(x, y) (-sin(y) .* cos(x));
%! problem_data.h = @(x, y, iside) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.velex = @(x, y) cat(1, ...
%!                   reshape (uxex (x,y), [1, size(x)]), ...
%!                   reshape (uyex (x,y), [1, size(x)]));
%! problem_data.pressex = @(x, y) (3 * x.^2 + sin(x .* y) - 1.239811742000564725943866);
%! problem_data.gradvelex = @test_stokes_square_bc_graduex;
%! method_data.element_name = 'ndl';   % Element type for discretization
%! method_data.degree       = [ 3  3];  % Degree of the splines (pressure space)
%! method_data.regularity   = [ 2  2];  % Regularity of the splines (pressure space)
%! method_data.nsub         = [10 10];  % Number of subdivisions
%! method_data.nquad        = [ 5  5];  % Points for the Gaussian quadrature rule
%! factor = 10;
%! method_data.Cpen = factor*(min(method_data.degree)+1);
%! [geometry, msh, space_v, vel, space_p, press] = ...
%!                        solve_stokes (problem_data, method_data);
%! error_l2_p = sp_l2_error (space_p, msh, press, problem_data.pressex);
%! [error_h1_v, error_l2_v] = ...
%!    sp_h1_error (space_v, msh, vel, problem_data.velex, problem_data.gradvelex);
%! assert (msh.nel, 100)
%! assert (space_p.ndof, 169)
%! assert (space_v.ndof, 644)
%! assert (error_l2_p, 1.96345892858352e-08, 1e-15)
%! assert (error_h1_v, 9.54443704246158e-08, 1e-15)
%! assert (error_l2_v, 1.54091389875648e-09, 1e-15)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function a = acosh(a)
%ACOSH        Gradient inverse hyperbolic cosine acosh(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 03/22/04     S.M. Rump  correction for complex input
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  wng = warning;
  warning off

  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = full(a.x(:));
  ax = 1 ./ ( sqrt( ax-1 ) .* sqrt( ax+1 ) );
  a.x = acosh(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  warning(wng)

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function test_suite = test_triangleArea3d
%TEST_TRIANGLEAREA3D  Test case for the file triangleArea3d
%
%   Test case for the file triangleArea3d

%   Example
%   test_triangleArea3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-08-23,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>

% rectangle triangle with sides 20 and 30 -> area is 300.
p1 = [10 20 30];
p2 = [30 20 30];
p3 = [10 50 30];

exp = 300;

area = triangleArea3d(p1, p2, p3);
testCase.assertEqual(exp, area);
area = triangleArea3d(p2, p3, p1);
testCase.assertEqual(exp, area);
area = triangleArea3d(p3, p1, p2);
testCase.assertEqual(exp, area);

area = triangleArea3d(p3, p2, p1);
testCase.assertEqual(exp, area);
area = triangleArea3d(p1, p3, p2);
testCase.assertEqual(exp, area);
area = triangleArea3d(p2, p1, p3);
testCase.assertEqual(exp, area);

area = triangleArea3d([p1; p2; p3]);
testCase.assertEqual(exp, area);

function test_XZ(testCase)

% rectangle triangle with sides 20 and 30 -> area is 300.
p1 = [10 30 20];
p2 = [30 30 20];
p3 = [10 30 50];
exp = 300;

area = triangleArea3d(p1, p2, p3);
testCase.assertEqual(exp, area);

area = triangleArea3d([p1; p2; p3]);
testCase.assertEqual(exp, area);

function test_YZ(testCase)

% rectangle triangle with sides 20 and 30 -> area is 300.
p1 = [30 10 20];
p2 = [30 30 20];
p3 = [30 10 50];
exp = 300;

area = triangleArea3d(p1, p2, p3);
testCase.assertEqual(exp, area);

area = triangleArea3d([p1; p2; p3]);
testCase.assertEqual(exp, area);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function pass = test_arithmetic( pref )
% Check the Chebfun2v constructor for simple arithmetic operations.
% Alex Townsend, March 2013.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3 * pref.cheb2Prefs.chebfun2eps;
j = 1;

% These function chosen so that scl does not change.
f = @(x,y) cos(x); f=chebfun2v(f,f);
g = @(x,y) sin(y); g=chebfun2v(g,g);
% exact answers.
plus_exact = @(x,y) cos(x) + sin(y);
plus_exact=chebfun2v(plus_exact, plus_exact);
minus_exact = @(x,y) cos(x) - sin(y);
minus_exact=chebfun2v(minus_exact, minus_exact);
mult_exact = @(x,y) cos(x).*sin(y);
mult_exact=chebfun2v(mult_exact, mult_exact);

pass(j) = norm(f + g - plus_exact) < tol; j=j+1;
pass(j) = norm(f - g - minus_exact) < tol; j=j+1;
pass(j) = norm(f.*g - mult_exact) < tol; j=j+1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_issue2193

% WALLTIME 00:10:00
% MEM 3gb
% DEPENDENCY ft_apply_montage ft_prepare_montage

load(dccnpath('/home/common/matlab/fieldtrip/data/test/issue2193.mat'));
data = ft_preprocessing(cfg, data_segmented);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function d = ddiff ( d1, d2 )

%*****************************************************************************80
%
%% DDIFF returns the signed distance to a region that is the difference of two regions.
%
%  Discussion:
%
%    The author comments that this is not the true signed distance function for
%    the difference set, in particular, around corners.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real D1, D2, the signed distances to region 1 and 2.
%
%    Output, real D, the signed distance to the region formed by
%    removing from region 1 its intersection with region 2.
%
  d = max ( d1, -d2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function benford_test ( )

%*****************************************************************************80
%
%% BENFORD_TEST tests BENFORD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BENFORD_TEST\n' );
  fprintf ( 1, '  BENFORD(I) is the Benford probability of the\n' );
  fprintf ( 1, '  initial digit sequence I.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I     BENFORD(I)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 9
    fprintf ( 1, '  %2d  %12f\n', i, benford(i) );
  end

  return
end


text:
function outfigdir1 = PlotSampleFramesWrapperRotate(behavior,rootdatadir,experiment_name,scoresfilestr,mainfly,otherflies,ts,ts_overlay,varargin)

annfilestr = 'movie.ufmf.ann';
colorpos = [.7,0,0];
colorneg = [0,0,.7];
border = 16;
bg_thresh = 10/255;
wmah = .5;
frac_a_back = 1;
dist_epsilon = .1;
figpos_example = [97 886 1800 650];
targettype = 'fly';
nframes_sample_bkgd = 20;

[annfilestr,colorpos,colorneg,border,...
  bg_thresh,wmah,frac_a_back,dist_epsilon,...
  figpos,hfig_base,outfigdir,doplotellipseoverlay,predictions,...
  targettype,nframes_sample_bkgd,...
  fg_thresh,ncolors_reference,intensityoff,...
  leftovers] = ...
  myparse_nocheck(varargin,...
  'annfilestr',annfilestr,...
  'colorpos',colorpos,...
  'colorneg',colorneg,...
  'border',border,...
  'bg_thresh',bg_thresh,...
  'wmah',wmah,...
  'frac_a_back',frac_a_back,...
  'dist_epsilon',dist_epsilon,...
  'figpos',figpos_example,...
  'hfig_base',0,...
  'outfigdir','.',...
  'doplotellipseoverlay',false,...
  'predictions',[],...
  'targettype',targettype,...
  'nframes_sample_bkgd',nframes_sample_bkgd,...
  'fg_thresh',.7,...
  'ncolors_reference',256,...
  'intensityoff',0);

flipim = false;
switch targettype,
  case 'fly',
    moviefilestr = 'movie.ufmf';
    trxfilestr = 'registered_trx.mat';
    trx = Trx('trxfilestr',trxfilestr,'moviefilestr',moviefilestr,'perframedir','perframe',leftovers{:});
  case 'mouse'
    moviefilestr = 'movie.seq';
    trxfilestr = 'trx.mat';
    trx = Trx('trxfilestr',trxfilestr,'moviefilestr',moviefilestr,'perframedir','perframe',leftovers{:});
  case 'larva'
    moviefilestr = 'movie.mov';
    trxfilestr = 'trx.mat';
    trx = Trx('trxfilestr',trxfilestr,'moviefilestr',moviefilestr,'perframedir','perframe',leftovers{:});
    flipim = true;
end
expdir = fullfile(rootdatadir,experiment_name);
trx.AddExpDir(expdir,'openmovie',false);
moviename = fullfile(expdir,moviefilestr);
global STOREDREADFRAME
global STOREDHEADERINFO
% readframe = STOREDREADFRAME;
% headerinfo = STOREDHEADERINFO;
[readframe,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
STOREDREADFRAME = readframe;
STOREDHEADERINFO = headerinfo;

if isempty(predictions),
  if ~exist(fullfile(expdir,scoresfilestr),'file'),
    fprintf('PREDICTIONS NOT LOADED FOR REAL FOR %s! FIX THIS\n',behavior);
    predictions = cell(1,trx.nflies);
    for i = 1:trx.nflies,
      predictions{i} = true(trx(i).nframes);
    end
  else
    scoresdata = load(fullfile(expdir,scoresfilestr));
    predictions = cellfun(@(x) x > 0,scoresdata.allScores.scores,'UniformOutput',false);
  end
end

if ~isempty(annfilestr) && exist(fullfile(rootdatadir,experiment_name,annfilestr),'file'),
  annfilename = fullfile(rootdatadir,experiment_name,annfilestr);
  [bkgdim,fg_thresh] = read_ann(annfilename,'background_center','n_bg_std_thresh');
  bkgdim = reshape(bkgdim,[headerinfo.nc,headerinfo.nr])'/255;
  fg_thresh = fg_thresh / 255;
else
  sampleims = nan(headerinfo.nr,headerinfo.nc,1,nframes_sample_bkgd);
  samplets_bkgd = round(linspace(1,headerinfo.nframes,nframes_sample_bkgd+2));
  samplets_bkgd = samplets_bkgd(2:end-1);
  for i = 1:nframes_sample_bkgd,
    imcurr = im2double(readframe(samplets_bkgd(i)));
    if i == 1 && size(imcurr,3) > 1,
      sampleims = repmat(sampleims,[1,1,3,1]);
    end
    sampleims(:,:,:,i) = imcurr;
  end
  bkgdim = median(sampleims,4);
end

sigma_bkgd = (fg_thresh-bg_thresh);

hfig = hfig_base+1;

[hfig,~,pfore] = PlotSampleFramesColorRotate(trx,readframe,bkgdim,predictions,mainfly,otherflies,ts,...
  'colorpos',colorpos,...
  'colorneg',colorneg,...
  'hfig',hfig,'figpos',figpos,...
  'colormap',{@jet,@jet},...
  'fg_thresh',fg_thresh,...
  'bg_thresh',bg_thresh,...
  'sigma_bkgd',sigma_bkgd,...
  'wmah',wmah,...
  'frac_a_back',frac_a_back,...
  'dist_epsilon',dist_epsilon,...
  'border',border,...
  'ncolors_reference',ncolors_reference,...
  'flipim',flipim,...
  'intensityoff',intensityoff);

shortname = regexprep(experiment_name,'^GMR_','');
shortname = regexprep(shortname,'_A[DE]_01_.*$','');

outfigdir1 = fullfile(outfigdir,sprintf('Example%s_%s_fly%02d_Frames%05dto%05d',behavior,shortname,mainfly,ts_overlay(1),ts_overlay(end)));
if ~exist(outfigdir1,'dir'),
  mkdir(outfigdir1);
end

SaveFigLotsOfWays(hfig,fullfile(outfigdir1,sprintf('Example%s_%s_fly%02d_frame%02d',behavior,shortname,mainfly,ts(1))));
for i = 1:size(pfore,3),
  imwrite(pfore(:,:,i),fullfile(outfigdir1,sprintf('Example%s_pfore_%s_fly%02d_frame%02d.png',behavior,shortname,mainfly,ts(i))));
end
%
% hfig = hfig_base+2;
% [hfig,~,pfore] = OverlayFrames(trx,readframe,bkgdim,predictions,mainfly,otherflies,ts_overlay,...
%   'colorpos',colorpos,'colorneg',colorneg,'border',border,...
%   'hfig',hfig,'figpos',figpos,...
%   'colormap',{@jet,@jet},...
%   'fg_thresh',fg_thresh,...
%   'bg_thresh',bg_thresh,...
%   'sigma_bkgd',sigma_bkgd,...
%   'wmah',wmah,...
%   'frac_a_back',frac_a_back,...
%   'dist_epsilon',dist_epsilon,...
%   'ncolors_reference',ncolors_reference,...
%   'flipim',flipim,...
%   'intensityoff',intensityoff);
%
% SaveFigLotsOfWays(hfig,fullfile(outfigdir1,sprintf('Example%s_Overlay_%s_fly%02d_Frames%05dto%05d',behavior,shortname,mainfly,ts_overlay(1),ts_overlay(end))));
% imwrite(pfore,fullfile(outfigdir1,sprintf('Example%s_Overlay_pfore_%s_fly%02d_frame%02d.png',behavior,shortname,mainfly,ts(i))));

if doplotellipseoverlay,
  hfig = hfig_base+3;
  if strcmp(targettype,'larva'),
    tmp = load(fullfile(expdir,trxfilestr),'trx');
    trx1 = tmp.trx;
    [hfig] = OverlayFrames_Skeletons(trx,trx1,readframe,bkgdim,predictions,mainfly,otherflies,ts,...
      'colorpos',colorpos,'colorneg',colorneg,'border',border,...
      'hfig',hfig,'figpos',figpos,...
      'colormap',{@jet,@jet},...
      'fg_thresh',fg_thresh,...
      'bg_thresh',bg_thresh,...
      'sigma_bkgd',sigma_bkgd,...
      'wmah',wmah,...
      'frac_a_back',frac_a_back,...
      'dist_epsilon',dist_epsilon,...
      'border',border,...
      'ncolors_reference',ncolors_reference,...
      'flipim',flipim);
  else
    [hfig] = OverlayFrames_Ellipses(trx,readframe,bkgdim,predictions,mainfly,otherflies,ts,...
      'colorpos',colorpos,'colorneg',colorneg,'border',border,...
      'hfig',hfig,'figpos',figpos,...
      'colormap',{@jet,@jet},...
      'fg_thresh',fg_thresh,...
      'bg_thresh',bg_thresh,...
      'sigma_bkgd',sigma_bkgd,...
      'wmah',wmah,...
      'frac_a_back',frac_a_back,...
      'dist_epsilon',dist_epsilon,...
      'border',border,...
      'ncolors_reference',ncolors_reference);
  end
  SaveFigLotsOfWays(hfig,fullfile(outfigdir1,sprintf('Example%s_OverlayEllipses_%s_fly%02d_Frames%05dto%05d',behavior,shortname,mainfly,ts_overlay(1),ts_overlay(end))));
end

author: kristinbranson

repo: JAABA


text:
function dwped = dtiInitPhaseDim(dwRawPed)
%
%  dwped = dtiInitPhaseDim(dwParams,dwRaw)
%
% Make sure there is a valid phase-encode dir, as this is crucial for
% eddy-current correction. If the phase_dim is not valid in the dwRaw or
% not set in dwParams.phaseEncode prompt the user for it.
%
% If the phase encode direction is not correct in dwRaw.phase_encode. Here
% we ask the user to tell us what it is. We could read this from the dicom
% if we were smart about it. It would be in the dicom header - we would
% have to know where the raw dicoms are and get them out if they're zipped.
% (dicominfo will read the dicom header).
%
%
% INPUTS
%       (dwRawPed) - passed in from dtiInit
% RETURNS
%       dwped - with a valid phase_dim
%
% Web Resources
%       mrvBrowseSVN('dtiInitPhaseDim');
%
% (C) Stanford VISTA, 8/2011 [lmp]
%

%%
%
% Set it to output what was passed in and check it
dwped = dwRawPed;

% Handle the case where the value is passed in as a string
if ischar(dwped)
    dwped = str2double(dwped);
end

% Check the value passed in: If the phase encode direction is not a valid
% value prompt the user until they give a valid value (default is 2 = 'col'
% =  'A/P').
while (dwped < 1 || dwped > 3)
    prompt = sprintf('Phase-encode dir is currently [ %d ], but must be 1 2 or 3. \n (1 = L/R "row", 2 = A/P "col"). \n Enter new value:\n', dwRawPed);
    resp = inputdlg(prompt, 'Set phase encode direction', 1, {'2'});
    if ~isempty(resp)
        dwped = round(str2double(resp{1}));
    elseif isempty(resp)
        error('User aborted! [dtiInitPhaseDim]');
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef CantileverBeamMeshCreator < handle

    properties (Access = public)
        connec
        coords
    end

    properties (Access = private)
        dim
        length
        height
    end

    methods (Access = public)

        function obj = CantileverBeamMeshCreator(cParams)
            obj.init(cParams)
        end

        function mesh = create(obj, xdiv, ydiv)
            switch obj.dim
                case '2D'
                    mesh = obj.create2Dcantilever(xdiv, ydiv);
                case '3D'
                    mesh = obj.create3Dcantilever(xdiv, ydiv);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.dim    = cParams.dim;
            obj.length = cParams.length;
            obj.height = cParams.height;
        end

        function mesh = create2Dcantilever(obj, xdiv, ydiv)
            obj.computeCoords2D(xdiv, ydiv);
            obj.computeConnec2D(xdiv, ydiv);
            mesh = obj.createMesh();
        end

        function mesh = create3Dcantilever(obj, xdiv, ydiv)
            obj.computeCoords3D(xdiv, ydiv);
            obj.computeConnec3D(xdiv, ydiv);
            mesh = obj.createMesh();
        end

        function coords = computeCoords2D(obj, xdiv, ydiv)
            x = linspace(0, obj.length, xdiv+1);
            y = linspace(0, obj.height, ydiv+1);
            [X,Y,Z] = meshgrid(x,y,0);
            fvc = surf2patch(X,Y,Z,'triangles');
            fvc.vertices(:,3) = []; % 2D
            coords = fvc.vertices;
            obj.coords = coords;
        end

        function computeConnec2D(obj, xdiv, ydiv)
            conn = [];
            for j = 0:1:xdiv-1
                for i = 1:1:ydiv
                    node1 = j*(ydiv+1) + i;
                    node2 = node1 + 1;
                    node3 = node1 + (ydiv+1);
                    node4 = node2 + (ydiv+1);
                    elem = [node1, node2, node4, node3];
                    conn = [conn; elem];
                end
            end
            obj.connec = conn;
        end

        function computeCoords3D(obj, xdiv, ydiv)
            x = linspace(0, obj.length,    xdiv+1);
            y = linspace(0, obj.height, ydiv+1);
            z = linspace(0, obj.height, ydiv+1);
            [X,Y,Z] = meshgrid(x,y,z);
            npnod = size(X,1)*size(X,2)*size(X,3);
            Xr = reshape(X, npnod,1);
            Yr = reshape(Y, npnod,1);
            Zr = reshape(Z, npnod,1);
            coor = [Xr, Yr, Zr];
            obj.coords = coor;
        end

        function computeConnec3D(obj, xdiv, ydiv)
            conn = [];
%             for z = 0:1:ydiv-1
%                 for j = 0:1:ydiv-1
%                     for i = 1:1:xdiv
%                         addZ = (xdiv+1)*(ydiv+1)*z;
%                         addZ1 = (xdiv+1)*(ydiv+1)*(z+1);
%                         node1 = j*(xdiv+1)+i + addZ;
%                         node2 = j*(xdiv+1)+i+1 + addZ;
%                         node3 = (j+1)*(xdiv+1)+i + addZ;
%                         node4 = (j+1)*(xdiv+1)+i+1 + addZ;
%                         node5 = j*(xdiv+1)+i + addZ1;
%                         node6 = j*(xdiv+1)+i+1 + addZ1;
%                         node7 = (j+1)*(xdiv+1)+i + addZ1;
%                         node8 = (j+1)*(xdiv+1)+i+1 + addZ1;
%                         elem = [node1, node2, node3, node4, ...
%                                 node5, node6, node7, node8];
%                         conn = [conn; elem];
%                     end
%                 end
%             end
            for z = 0:1:ydiv
                for j = 0:1:xdiv-1
                    for i = 1:1:ydiv
                        addZ  = (xdiv+1)*(ydiv+1)*z;
                        node1 = j*(ydiv+1) + i;
                        node2 = node1 + 1;
                        node3 = node1 + (ydiv+1);
                        node4 = node2 + (ydiv+1);
                        node5 = j*(ydiv+1) + i + addZ;
                        node6 = node5 + 1;
                        node7 = node5 + (ydiv+1);
                        node8 = node6 + (ydiv+1);
                        elem = [node1, node2, node4, node3, ...
                                node5, node6, node7, node8];
                        conn = [conn; elem];
                    end
                end
            end
            obj.connec = conn;
        end

        function mesh = createMesh(obj)
            m.coord  = obj.coords;
            m.connec = obj.connec;
            mesh = Mesh(m);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [H,Q,M,W]=matExpIntFam(T,A,B,Qc)
%%MATEXPINTFAM This function explicitly solves a number of real matrix
%       integral problems that arise when solving for optimal linear
%       regulators. The integrals are:
%       H=int_0^T expm(A*s)*B ds
%       Q=int_0^T expm(A'*s)*Qc*expm(A*s) ds
%       M=int_0^T expm(A'*s)*Qc*H(s) ds
%       W=int_0^T H(s)'*Qc*H(s) ds
%       where s is the scalar parameter of integration, expm is the matrix
%       exponential function int refers to a definite integral, and the
%       function
%       H(s) is defined to be
%       H(s)=(int_0^sexpm(A*tau) dtau)*B
%
%INPUTS: T The real, positive, scalar upper bound of the integrals.
%        A A real nXn matrix.
%        B A real nXp matrix. If this parameter is omitted or an empty
%          matrix is passed, then it is assumed that p=0, in which case H,
%          M, and W will be empty matrices.
%       Qc A symmetric positive (semi)definite nXn matrix. If omitted or an
%       empty matrix is passed, then it is assumed that Qc=eye(n,n);
%
%OUTPUTS: H,Q,M,W The values of the desired matrix integrals, which are
%         defined above
%
%Expressions for obtaining the solutions to the integrals by evaluating a
%single matrix exponential are given in [1]. The type of linear regulator
%problem addressed by these integrals is given in [2].
%
%REFERENCES:
%[1] C. F. Van Loan, "Computing Integrals Involving the Matrix
%    Exponential," IEEE Transactions on Automatic Control, vol. AC-23, no.
%    3, pp. 395-404, Jun. 1967.
%[2] E. S. Armstrong and A. K. Caglayan, "An algorithm for the weighting
%    matrices in the sampled-data optimal linear regulator problem,"
%    National Aeronautics and Space Administration, Langley Research
%    Center, Hampton, VA, Tech. Rep. NASA TN D-8372, Dec. 1976.
%
%September 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

if(nargin<4||isempty(Qc))
    Qc=eye(n,n);
end

if(isempty(B)||nargin<2)
    p=0;
    B=zeros(n,p);
else
    p=size(B,2);
end

C=[-A.',        eye(n,n), zeros(n,n+p);
    zeros(n,n), -A.',     Qc,   zeros(n,p);
    zeros(n,2*n),         A,    B;
    zeros(p,3*n+p)];

CExp=expm(C*T);

F3=CExp((2*n+1):(3*n),(2*n+1):(3*n));
G2=CExp((n+1):(2*n),(2*n+1):(3*n));
G3=CExp((2*n+1):(3*n),(3*n+1):(3*n+p));
H2=CExp((n+1):(2*n),(3*n+1):(3*n+p));
K1=CExp(1:n,(3*n+1):(3*n+p));

H=G3;
Q=F3.'*G2;
M=F3.'*H2;

temp=B.'*F3.'*K1;
W=temp+temp.';

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


text:
function ndx = read_mat(infilename)
% Reads a struct from a mat file and constructs an Ndx object.
% Inputs:
%   infilename: The name for the mat file to read.
% Outputs:
%   ndx: An Ndx object created from the information in the mat
%     file.

assert(nargin==1)
assert(isa(infilename,'char'))

load(infilename,'ndx');
ndx = Ndx(ndx);

assert(ndx.validate())

author: nesl


text:
function [thrs,cntR,sumR,cntP,sumP,V] = edgesEvalImg( E, G, varargin )
% Calculate edge precision/recall results for single edge image.
%
% Enhanced replacement for evaluation_bdry_image() from BSDS500 code:
%  http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/
% Uses same format and is fully compatible with evaluation_bdry_image.
% Given default prms results are *identical* to evaluation_bdry_image.
%
% In addition to performing the evaluation, this function can optionally
% create a visualization of the matches and errors for a given edge result.
% The visualization of edge matches V has the following color coding:
%  green=true positive, blue=false positive, red=false negative
% If multiple ground truth labels are given the false negatives have
% varying strength (and true positives can match *any* ground truth).
%
% This function calls the mex file correspondPixels. Pre-compiled binaries
% for some systems are provided in /private, source for correspondPixels is
% available as part of the BSDS500 dataset (see link above). Note:
% correspondPixels is computationally expensive and very slow in practice.
%
% USAGE
%  [thrs,cntR,sumR,cntP,sumP,V] = edgesEvalImg( E, G, [prms] )
%
% INPUTS
%  E          - [h x w] edge probability map (may be a filename)
%  G          - file containing a cell of ground truth boundaries
%  prms       - parameters (struct or name/value pairs)
%   .out        - [''] optional output file for writing results
%   .thrs       - [99] number or vector of thresholds for evaluation
%   .maxDist    - [.0075] maximum tolerance for edge match
%   .thin       - [1] if true thin boundary maps
%
% OUTPUTS
%  thrs       - [Kx1] vector of threshold values
%  cntR,sumR  - [Kx1] ratios give recall per threshold
%  cntP,sumP  - [Kx1] ratios give precision per threshold
%  V          - [hxwx3xK] visualization of edge matches
%
% EXAMPLE
%
% See also edgesEvalDir
%
% Structured Edge Detection Toolbox      Version 3.01
% Code written by Piotr Dollar, 2014.
% Licensed under the MSR-LA Full Rights License [see license.txt]

% get additional parameters
dfs={ 'out','', 'thrs',99, 'maxDist',.0075, 'thin',1 };
[out,thrs,maxDist,thin] = getPrmDflt(varargin,dfs,1);
if(any(mod(thrs,1)>0)), K=length(thrs); thrs=thrs(:); else
  K=thrs; thrs=linspace(1/(K+1),1-1/(K+1),K)'; end

% load edges (E) and ground truth (G)
if(all(ischar(E))), E=double(imread(E))/255; end
G=load(G); G=G.groundTruth; n=length(G);
for g=1:n, G{g}=double(G{g}.Boundaries); end

% evaluate edge result at each threshold
Z=zeros(K,1); cntR=Z; sumR=Z; cntP=Z; sumP=Z;
if(nargout>=6), V=zeros([size(E) 3 K]); end
for k = 1:K
  % threshhold and thin E
  E1 = double(E>=max(eps,thrs(k)));
  if(thin), E1=double(bwmorph(E1,'thin',inf)); end
  % compare to each ground truth in turn and accumualte
  Z=zeros(size(E)); matchE=Z; matchG=Z; allG=Z;
  for g = 1:n
    [matchE1,matchG1] = correspondPixels(E1,G{g},maxDist);
    matchE = matchE | matchE1>0;
    matchG = matchG + double(matchG1>0);
    allG = allG + G{g};
  end
  % compute recall (summed over each gt image)
  cntR(k) = sum(matchG(:)); sumR(k) = sum(allG(:));
  % compute precision (edges can match any gt image)
  cntP(k) = nnz(matchE); sumP(k) = nnz(E1);
  % optinally create visualization of matches
  if(nargout<6), continue; end; cs=[1 0 0; 0 .7 0; .7 .8 1]; cs=cs-1;
  FP=E1-matchE; TP=matchE; FN=(allG-matchG)/n;
  for g=1:3, V(:,:,g,k)=max(0,1+FN*cs(1,g)+TP*cs(2,g)+FP*cs(3,g)); end
  V(:,2:end,:,k) = min(V(:,2:end,:,k),V(:,1:end-1,:,k));
  V(2:end,:,:,k) = min(V(2:end,:,:,k),V(1:end-1,:,:,k));
end

% if output file specified write results to disk
if(isempty(out)), return; end; fid=fopen(out,'w'); assert(fid~=1);
fprintf(fid,'%10g %10g %10g %10g %10g\n',[thrs cntR sumR cntP sumP]');
fclose(fid);

end

author: pdollar

repo: edges

save_path: github-repos/MATLAB/pdollar-edges


text:
%kClassify 'classifies a series of images according to a number of color vectors '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros Classify.pane file
%
% Parameters:
% InputFile: i1 'Images (along elements)', required: 'Series of images (also 3D) stacked along elment direction'
% InputFile: i2 'Colors', required: 'Color vectors (along X) stacked along element direcition'
% Double: t 'DAPI threshold', default: 10: 'Element 0 is assumed to contain the DAPI channel. Pixels below threshold are ignored.'
% Double: minq 'Min Quality', default: 0.97: 'Defines the minimal quality for single channel matching. If not reached, method reverts to two channels'
% Integer: co 'Color Optimize Iterations', default: 1: 'If selected, some iterations for finding the optimal color vectors are performed'
% Double: mincq 'Min Color Quality', default: 0.8: 'Defines the minimum quality for adjusting the colors.'
% OutputFile: o1 'Classified', required: 'Resulting classification image'
% OutputFile: o2 'Quality Img', required: 'A classification quality is given for every pixel'
% OutputFile: o3 'Output Colors', optional: 'The (eventually optimized) colors'
%
% Example: [o1, o2, o3] = kClassify({i1, i2}, {'i1','';'i2','';'t',10;'minq',0.97;'co',1;'mincq',0.8;'o1','';'o2','';'o3',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% Classify - classifies a series of images according to a number of color vectors
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ama] = loadama(filename)

% LOADAMA read an inverted A-matrix and associated geometry information
% from an ama file that was written by Tom Oostendorp's DIPOLI
%
% Use as
%   [ama] = loadama(filename)
%
% See also LOADTRI, LOADMAT

% Copyright (C) 2005, Robert Oostenveld
% Copyright (C) 2021, Thom Oostendorp
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

fid = fopen(filename, 'rb', 'ieee-le');
version = fread(fid, 1, 'int');

% TFO 2021-08-10
% version 11 and higher of dipoli uses double precision so all 'float's have been replaced by 'float64's
switch version
  case 10
    floatFormat='float';
  case 11
    floatFormat='float64';
  otherwise
    ft_error('%s is either not an inverted A matrix, or neither version 10 nor 11', filename);
end

mode = fread(fid, 1, 'int');
ngeo = fread(fid, 1, 'int');

totpnt = 0;
tottri = 0;
nrow   = 0;

% read the boundaries
geo  = [];
for i=1:ngeo
  geo(i).name    = char(fread(fid, [1 80], 'uchar'));
  geo(i).npos    = fread(fid, 1, 'int');
  geo(i).pos     = fread(fid, [3 geo(i).npos], floatFormat)';
  geo(i).ntri    = fread(fid, 1, 'int');
  geo(i).tri     = fread(fid, [3 geo(i).ntri], 'int')' + 1;  % MATLAB indexing starts at 1
  geo(i).sigmam  = fread(fid, 1, floatFormat);
  geo(i).sigmap  = fread(fid, 1, floatFormat);
  geo(i).geocon  = fread(fid, ngeo, 'int');
  geo(i).deflat  = fread(fid, ngeo, floatFormat);
  totpnt = totpnt + geo(i).npos;
  tottri = tottri + geo(i).ntri;
end

% read the electrodes
if mode~=1
  elec.name    = char(fread(fid, [1 80], 'uchar'));
  elec.npos    = fread(fid, 1, 'int');
  for i=1:(elec.npos+1)
    elec.el(i).tri  = fread(fid, 1, 'int') + 1; % MATLAB indexing starts at 1
    % TFO 2021-08-10
    % In the 64-bit version, each field is padded to a multiple of 8 bytes, so we need to read 4 dummy bytes
    if version>10
      dum = fread(fid, 4, 'char');
    end
    elec.el(i).la   = fread(fid, 1, floatFormat);
    elec.el(i).mu   = fread(fid, 1, floatFormat);
    elec.el(i).name = char(fread(fid, [1 10], 'char'));
    % the ELECTRODE c-structure is padded to word boundaries, i.e. to 4 bytes
    dum = fread(fid, 2, 'char');
    % TFO 2021-08-10
    % In the 64-bit version, each field is padded to a multiple of 8 bytes, so we need to read 4 more dummy bytes
    if version>10
      dum = fread(fid, 4, 'char');
    end
  end
  elec.vertex  = fread(fid, 1, 'int');
  elec.surface = fread(fid, 1, 'int');
  nrow = nrow + elec.npos;
else
  elec = [];
end

% read the gradiometers
if mode~=0
  ft_error('gradiometers not yet implemented');
else
  grad = [];
end

% read the inverted A-matrix
bi = fread(fid, [totpnt nrow], floatFormat)';

% read the isolated source compartment information, if present
iso_sur    = fread(fid, 1, 'int') + 1;  % Matlab indexing starts at 1
inner_only = fread(fid, 1, 'int');
if iso_sur~=0
  iso_totpnt = geo(iso_sur).npos;
  iso_b      = fread(fid, [iso_totpnt iso_totpnt], floatFormat)';
else
  iso_b = [];
end

fclose(fid);
clear floatFormat;
clear dum;
clear i;

% put all local variables into a structure, this is a bit unusual programming style
% the output structure is messy, but contains all relevant information
tmp = whos;
ama = [];
for i=1:length(tmp)
  if isempty(strmatch(tmp(i).name, {'tmp', 'fid', 'ans', 'handles'}))
    ama = setfield(ama, tmp(i).name, eval(tmp(i).name));
  end
end


text:
function Gxy = femGreenKernel(X,Y,green,k)
%+========================================================================+
%|                                                                        |
%|              OPENFEM - LIBRARY FOR FINITE ELEMENT METHOD               |
%|           openFem is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : femGreenKernel.m                              |
%|    #    |   VERSION    : 0.55                                          |
%|   _#_   |   AUTHOR(S)  : M. Aussal & F. Alouges & M. Averseng          |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.05.2019                                    |
%| ( === ) |   SYNOPSIS   : Usefull green kernel functions                |
%|  `---'  |                                                              |
%+========================================================================+

% Security
if (size(X,2) ~= 3) || (size(Y,2) ~= 3)
    error('femGreenKernel.m : unavailable case')
end
if isempty(k)
    k = 0;
end

% Distances between particles
Rxy = sqrt( ...
    (X(:,1) - Y(:,1)).^2 + ...
    (X(:,2) - Y(:,2)).^2 + ...
    (X(:,3) - Y(:,3)).^2 );

% For empty wave-number
if isempty(k)
    k = 0;
end

% Green kernel definition
if strcmp(green,'[1/r]')
    Gxy = 1./Rxy;

elseif strcmp(green,'dx[1/r]')
    Gxy = - 1 ./ (Rxy.^3);

elseif strcmp(green,'dy[1/r]')
    Gxy = 1 ./ (Rxy.^3);

elseif strcmp(green,'[exp(ikr)/r]')
    Gxy = exp(1i*k*Rxy)./Rxy;

elseif strcmp(green,'dx[exp(ikr)/r]')
    Gxy = (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) ./ (Rxy.^2);

elseif strcmp(green,'dy[exp(ikr)/r]')
    Gxy = - (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) ./ (Rxy.^2);

elseif strcmp(green(1:end-1),'gradx[1/r]')
    j = str2double(green(end));
    Gxy = - (X(:,j)-Y(:,j)) ./ (Rxy.^3);

elseif strcmp(green(1:end-1),'grady[1/r]')
    j = str2double(green(end));
    Gxy = (X(:,j)-Y(:,j)) ./ (Rxy.^3);

elseif strcmp(green(1:end-1),'gradx[exp(ikr)/r]')
    j = str2double(green(end));
    Gxy = (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) .* ...
        (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green(1:end-1),'grady[exp(ikr)/r]')
    j = str2double(green(end));
    Gxy = - (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) .* ...
        (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green,'[log(r)]')
    Gxy = log(Rxy);

elseif strcmp(green,'[H0(kr)]')
    Gxy = besselh(0,k*Rxy);

elseif strcmp(green(1:end-1),'gradx[log(r)]')
    j = str2double(green(end));
    Gxy = (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green(1:end-1),'grady[log(r)]')
    j = str2double(green(end));
    Gxy = - (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green(1:end-1),'gradx[H0(kr)]')
    j   = str2double(green(end));
    Gxy = - k * besselh(1,k*Rxy) .* (X(:,j)-Y(:,j)) ./ Rxy;

elseif strcmp(green(1:end-1),'grady[H0(kr)]')
    j   = str2double(green(end));
    Gxy = k * besselh(1,k*Rxy) .* (X(:,j)-Y(:,j)) ./ Rxy;

elseif strcmp(green(1:end-2),'[ij/r+rirj/r^3]')
    i = str2double(green(end-1));
    j = str2double(green(end));
    Gxy = (i==j)./Rxy + (X(:,i)-Y(:,i)).*(X(:,j)-Y(:,j))./(Rxy.^3);

elseif strcmp(green(1:end-3),'[rirjrk/r^5]')
    i = str2double(green(end-2));
    j = str2double(green(end-1));
    k = str2double(green(end));
    Gxy = (X(:,i)-Y(:,i)).*(X(:,j)-Y(:,j)).*(X(:,k)-Y(:,k))./(Rxy.^5);

else
    error('Error in femGreenKernel.m : unknown green kernel')
end

% Singularity
if strcmp(green,'[exp(ikr)/r]')
    Gxy(Rxy<1e-12) = 0 + 1i*k;
elseif strcmp(green,'[H0(kr)]')
    gamma         = 0.5772156649;
    Gxy(Rxy<1e-12) = 1 + 1i*(2/pi*(gamma+log(k/2)));
else
    Gxy(Rxy<1e-12) = 0;
end
end


text:
  function xs = alg_art1(x, At, y, varargin)
%|function xs = alg_art1(x, At, y, [options])
%| classical ART algorithm, aka, Kaczmarz algorithm; tries to solve y=Ax
%|
%| in
%|	x	[np 1]		initial guess, possibly empty
%|	Caution: x must be in the range of A' for convergence!
%|	At	[np nd]		(hermitian) *transpose* of system matrix
%|	y	[nb na]		measurement
%|
%| option
%|	wi	[nb na]		weights
%|	niter			# of iterations
%|	isave			default: [] 'last'
%|	pixmin
%|	pixmax
%|	anorms	[nd 1]		see below
%|	eps
%|
%| out
%|	xs	[np niter]	iterates
%|
%| Copyright 2006-4-2, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test'), alg_art1_test, return, end
if nargin < 3, ir_usage, end
if isempty(x), x = zeros(nrow(A),1); end

% defaults
arg.niter = 1;
arg.isave = [];
arg.anorms = [];
arg.eps = eps;
%arg.pixmax = inf;
%arg.pixmin = -inf;

% options
arg = vararg_pair(arg, varargin);

arg.isave = iter_saver(arg.isave, arg.niter);

[nb na] = size(y);
starts = subset_start(na);

% For WLS, premultiply y and postmultiply At by W^{1/2}
%Wh = spdiag(sqrt(wi(:)), 'nowarn');
%y = Wh * y(:);
%At = At * Wh;

% weighted row norms
if isempty(arg.anorms)
	arg.anorms = sum(At.^2); % | e_i' A |^2
end

iglist = col(outer_sum(1:nb, (starts-1)*nb));
iglist = iglist(arg.anorms(iglist) ~= 0);

adenom = arg.anorms + eps; % trick:

np = length(x);
xs = zeros(np, length(arg.isave));
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = x;
end

ticker(mfilename, 1, arg.niter)

for iter=1:arg.niter
	ticker(mfilename, iter, arg.niter)

	for ii=iglist'
		ai = At(:,ii);
		step = (y(ii) - ai' * x) / adenom(ii);
		x = x + ai * step;

%		todo: try following approach to see if faster
%		[j ignore ai] = find(At(:,ii));
%		xj = x(j);
%		step = (y(ii) - ai' * xj) / adenom(ii);
%		x(j) = xj - step * ai;
	end

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function f = imag(f)
%IMAG   Imaginary part of a CHEBTECH.
%   IMAG(F) is the imaginary part of F.
%
% See also REAL, ISREAL, CONJ.

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

% Compute imaginary part of the coefficients:
f.coeffs = imag(f.coeffs);

if ( ~any(f.coeffs(:)) )
    % Input was real, so output a zero CHEBTECH:
    f = f.make(zeros(1, size(f.coeffs, 2)));
    f.ishappy = 1;
end

end


text:
function Y=AveEmbedding(dataset,parameters)
dataStr=['./../dataset/',dataset,'-STC2.mat'];
load(dataStr);
parameters.vocSize = size_vocab;
% Step 1. Generate word vector sets
CR_E = randi([-25,25],parameters.wordDim,parameters.vocSize)/100;
disp(strcat('Number of weights E:',num2str(size(CR_E))));
vocab_emb_length = length(vocab_emb_Word2vec_48(1,:));
if vocab_emb_length > size_vocab
    error(['Error, and the size fo vocab_emb is:',vocab_emb_length])
end
CR_E(1:parameters.wordDim,vocab_emb_Word2vec_48_index) = vocab_emb_Word2vec_48(1:parameters.wordDim,1:vocab_emb_length);
% Step 2. Compute TF-IDF
if parameters.weightMode
    fea_All=tf_idf(fea_All);
end
% Step 3. Average Embedding
textSize = length(fea_All(:,1));
fea_vector =[];
for i=1:textSize
    tmp_fea_vector_weight = repmat(fea_All(i,find(fea_All(i,:)>0)),parameters.wordDim,1);
    tmp_fea_vector_matrix = CR_E(:,find(fea_All(i,:)>0)) .* tmp_fea_vector_weight;
    tmp_fea_vector = sum(tmp_fea_vector_matrix,2);
    fea_vector(i,:) = tmp_fea_vector';
end
% Step 4. Normalize features
Y = normalize(fea_vector);
end

author: jacoxu


text:
function r = issparse(c)
%ISSPARSE     Returns 1 if c is sparse
%
%  r = issparse(c);
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if c.complex
    r = issparse(c.mid);
  else
    r = issparse(c.inf);
  end


text:
 function x = dtft2_adj(X, omega, N1, N2, n_shift, useloop)
%function x = dtft2_adj(X, omega, N1, N2, n_shift, useloop)
%|
%| Compute adjoint of 2D DTFT for spectrum X at frequency locations omega
%| in
%|	X	[M,L]		2D DTFT values
%|	omega	[M,2]		frequency locations (radians)
%|	n_shift [2,1]		use [0:N-1]-n_shift (default [0 0])
%|	useloop			1 to reduce memory use (slower)
%| out
%|	x	[N1,N2,L]	signal values
%|
%| Requires enough memory to store M * (N1*N2) size matrices (for testing)
%|
%| Copyright 2001-9-17, Jeff Fessler, University of Michigan

if nargin == 1 && streq(X, 'test'), dtft2_adj_test, return, end
if nargin < 2, ir_usage(), end

if ~isvar('n_shift') || isempty(n_shift), n_shift = [0 0]; end
if ~isvar('useloop') || isempty(useloop), useloop = 0; end

[nn1, nn2] = ndgrid([0:(N1-1)]-n_shift(1), [0:(N2-1)]-n_shift(2));

if useloop % loop way: slower but less memory

	M = length(omega);
	x = zeros(N1,N2,ncol(X)); % [N1 N2 M]
	t1 = 1i * nn1;
	t2 = 1i * nn2;
	for ii=1:M
		x = x + exp(omega(ii,1)*t1 + omega(ii,2)*t2) * X(ii,:);
	end

else % non-loop way
	x = exp(1i*(nn1(:)*omega(:,1)' + nn2(:)*omega(:,2)')) * X; % [N1*N2 L]
	x = reshape(x, [N1 N2 numel(x)/N1/N2]); % [N1 N2 L]
end

function dtft2_adj_test()
N1 = 4; N2 = 6;
n_shift = [2 1];
% test with uniform frequency locations:
o1 = 2*pi*[0:(N1-1)]'/N1;
o2 = 2*pi*[0:(N2-1)]'/N2;
[o1, o2] = ndgrid(o1, o2);
X = o1 + o2; % test spectrum
om = [o1(:) o2(:)];
xd = dtft2_adj(X(:), om, N1, N2, n_shift);
xl = dtft2_adj(X(:), om, N1, N2, n_shift, 1);
printm('loop max %% difference = %g', max_percent_diff(xl,xd))
Xp = X .* reshape(exp(-1i * om * n_shift(:)), size(X));
xf = ifft2(Xp) * N1 * N2;
printm('ifft max %% difference = %g', max_percent_diff(xf,xd))


text:
%--- help for intercept_column ---
%
%  intercept_column - builds a scalar splanar object from a
%  vectorized splanar object. The pointer argument points the
%  element in the vector to be used.
%

author: jmaih


text:
% Shark Demo
%
% Copyright (c) by Lorenzo Torresani, Stanford University
%
% A demo of Non-Rigid Structure From Motion on artificial shark sequence
%
%
% The 3D reconstruction technique is based on the following paper:
%
%  Lorenzo Torresani, Aaron Hertzmann and Christoph Bregler,
%     Learning Non-Rigid 3D Shape from 2D Motion, NIPS 16, 2003
%  http://cs.stanford.edu/~ltorresa/projects/learning-nr-shape/
%
%
% Function em_sfm implements the algorithms "EM-Gaussian" and "EM-LDS" described
% in the paper
%
% I recommend that you try to compile the CMEX code for the function computeH:
% type 'mex computeH.c' in the Matlab Command Window ('mex computeH.c -l matlb' under Unix)
%

% loads the matrix P3_gt containing the ground thruth data: P3_gt([t t+T t+2*T],:) contains the 3D coordinates of the J points at time t
% (T is the number of frames, J is the number of points)
load('jaws.mat');
[T, J] = size(P3_gt); T = T/3;

% 2D motion resulting from orthographic projection (Eq (1))
p2_obs = P3_gt(1:2*T, :);

% runs the non-rigid structure from motion algorithm
use_lds = 1;
max_em_iter = 60;
tol = 0.0001;
K = 2; % number of deformation shapes
Zcoords_gt = P3_gt(2*T+1:3*T,:) - mean(P3_gt(2*T+1:3*T,:),2)*ones(1,J);
Zdist = max(Zcoords_gt,[],2) - min(Zcoords_gt,[],2); % size of the 3D shape along the Z axis for each time frame
MD = zeros(T,J);

[P3, S_hat, V, RO, Tr, Z] = em_sfm(p2_obs, MD, K, use_lds, tol, max_em_iter);

%% Compares it with ground truth.
% Note that there are still 2 unresolvable ambiguities:
% 1. depth direction (i.e. the shape could be "flipped" along the Z axis) -> we test both possibilities
% 2. Z translation                                                        -> we subtract the mean of the Z coords to evaluate reconstruction results
Zcoords_em = P3(2*T+1:3*T,:) - mean(P3(2*T+1:3*T,:),2)*ones(1,J);

Zerror1 = mean( mean(abs(Zcoords_em - Zcoords_gt), 2)./Zdist );
Zerror2 = mean( mean(abs(-Zcoords_em - Zcoords_gt), 2)./Zdist );

if Zerror2 < Zerror1,
   avg_zerror = 100*Zerror2;
   P3(2*T+1:3*T,:) = -(P3(2*T+1:3*T,:) - mean(P3(2*T+1:3*T,:),2)*ones(1,J));
else
   avg_zerror = 100*Zerror1;
   P3(2*T+1:3*T,:) = P3(2*T+1:3*T,:) - mean(P3(2*T+1:3*T,:),2)*ones(1,J);
end
fprintf('Average reconstruction error in Z: %f%%\n', avg_zerror);

vis_reconstruction(P3_gt, P3);

author: TadasBaltrusaitis

repo: OpenFace


text:
function value = p05_exact ( dim_num )

%*****************************************************************************80
%
%% P05_EXACT returns the exact integral for problem 05.
%
%  Discussion:
%
%    The exact value is given only for DIM_NUM = 1, 2, 3, 4 or 5.
%    For other cases, the value R8_HUGE is returned instead.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral,
%    or R8_HUGE if the exact value is not known.
%
  if ( dim_num == 1 )

    value = log ( 3.0 );

  elseif ( dim_num == 2 )

    value = 5.0 * log ( 5.0 ) - 6.0 * log ( 3.0 );

  elseif ( dim_num == 3 )

    value = 0.5 * ( 49.0 * log ( 7.0 ) ...
      - 75.0 * log ( 5.0 ) + 27.0 * log ( 3.0 ) );

  elseif ( dim_num == 4 )

    value = 225.0 * log ( 3.0 ) + 125.0 * log ( 5.0 ) ...
      - 686.0 * log ( 7.0 ) / 3.0;

  elseif ( dim_num == 5 )

    value = ( ...
      - 65205.0 * log ( 3.0 ) ...
      - 6250.0 * log ( 5.0 ) ...
      + 24010.0 * log ( 7.0 ) ...
      + 14641.0 * log ( 11.0 ) ) / 24.0;

  else

    value = r8_huge ( );

  end

  return
end


text:
function trans = createTranslation(varargin)
%CREATETRANSLATION Create the 3*3 matrix of a translation.
%
%   TRANS = createTranslation(DX, DY);
%   Returns the translation corresponding to DX and DY.
%   The returned matrix has the form :
%   [1 0 TX]
%   [0 1 TY]
%   [0 0  1]
%
%   TRANS = createTranslation(VECTOR);
%   Returns the matrix corresponding to a translation by the vector [x y].
%
%
%   See also:
%   transforms2d, transformPoint, createRotation, createScaling
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2004.
%

%   HISTORY
%   22/04/2009: rename as createTranslation

% process input arguments
if isempty(varargin)
    tx = 0;
    ty = 0;
elseif length(varargin)==1
    var = varargin{1};
    tx = var(1);
    ty = var(2);
else
    tx = varargin{1};
    ty = varargin{2};
end

% create the matrix representing the translation
trans = [1 0 tx ; 0 1 ty ; 0 0 1];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [nbcd,bcdof] = BoundaryConditions(sdof,bc)

%--------------------------------------------------------------------------
% Purpose :
%         To get the arrested degree's of freedom for the beam depending on
% type of the boundary conditions
%
% Synopsis :
%          [nbcd,bcdof] = BoundaryConditions(sdof,bc)
%
% Variable Description:
% INPUT parameters:
%           sdof : system degrees of freedom
%           bc : boundary condition type
%
% OUTPUT PARAMETERS :
%           bcdof : boundary degrees of freedom
%           nbcd : number of boundary conditions
%
%--------------------------------------------------------------------------

 if bc == 'c-c'     % clamped-clamped beam
    bcdof = [1 2 sdof-1 sdof] ;
    nbcd = length(bcdof) ;

 elseif bc == 'c-f' % clamped-free beam
    bcdof = [1 2] ;
    nbcd = length(bcdof) ;

 elseif bc == 'c-s'     % clamped-supported beam
     bcdof = [1 2 sdof-1] ;
     nbcd = length(bcdof) ;

 elseif bc == 's-s'     % supported-supported beam
     bcdof = [ 1 sdof-1] ;
     nbcd = length(bcdof) ;

 end

author: Sable

repo: mcbench-benchmarks


text:
function triangulation_test08 ( )

%*****************************************************************************80
%
%% TEST08 tests R8TRIS2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2005
%
%  Author:
%
%    John Burkardt
%
  node_num = 9;
  dim_num = 2;

  node_xy = [ ...
       0.0, 0.0; ...
       0.0, 1.0; ...
       0.2, 0.5; ...
       0.3, 0.6; ...
       0.4, 0.5; ...
       0.6, 0.4; ...
       0.6, 0.5; ...
       1.0, 0.0; ...
       1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  R8TRIS2 computes the Delaunay triangulation of\n' );
  fprintf ( 1, '    a set of nodes in 2D.\n' );
%
%  Set up the Delaunay triangulation.
%
  [ triangle_num, triangle_node, triangle_neighbor ] = r8tris2 ( ...
    node_num, node_xy );

  triangulation_order3_print ( node_num, triangle_num, node_xy, ...
    triangle_node, triangle_neighbor );

  return
end


text:
function model = rmSearchFit_oneGaussianNonlinear(model, data, params, wProcess, t)
% rmSearchFit_oneGaussianNonlinear - wrapper for 'fine' one Gaussian fit
%
% model = rmSearchFit_oneGaussianNonlinear(model, data, params, wProcess, t);
%
% 2008/01 SOD: split of from rmSearchFit.
% 2010/02 SOD: cleanup.
% 2015/02 JW: branched from rmSearchFit_oneGaussian, now includes
%             non-linear model (see Kay et al, 2013, on Compressive
%             Spatial Summation)

% fminsearch options
searchOptions = params.analysis.fmins.options;
expandRange   = params.analysis.fmins.expandRange;

% convert to double just in case
params.analysis.X = double(params.analysis.X);
params.analysis.Y = double(params.analysis.Y);

params.analysis.allstimimages_unconvolved = double(params.analysis.allstimimages_unconvolved);
data = double(data);

% get starting upper and lower range and reset TolFun
% (raw rss computation (similar to norm) and TolFun adjustments)
model.s = model.s_major;
[range, TolFun] = rmSearchFit_range(params,model,data);

% GLU 2021-10-14: not sute if best option, but I am going to restrict here
%                 the exponent
%                 If we are here we know that we are fitting the nonlinear css
if params.analysis.fixcssexp ~= 0
   range.start(4,:) = params.analysis.fixcssexp * ones(size(range.start(4,:)));
   range.lower(4,:) = params.analysis.fixcssexp * ones(size(range.start(4,:)));;
   range.upper(4,:) = params.analysis.fixcssexp * ones(size(range.start(4,:)));;
end

% amount of negative fits
nNegFit  = 0;
vethresh = params.analysis.fmins.vethresh;
trends   = t.trends;
t_id     = t.dcid+1;

%-----------------------------------
% Go for each voxel
%-----------------------------------
progress = 0;tic;
for ii = 1:numel(wProcess),
    % progress monitor (10 dots)
    if floor(ii./numel(wProcess)*10)>progress,
        % print out estimated time left
        if progress==0,
            esttime = toc.*10;
            if floor(esttime./3600)>0,
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d hours.\n',...
                    mfilename,numel(wProcess),ceil(esttime./3600));
            else
                fprintf(1,'[%s]:Estimated processing time: %d voxels: %d minutes.\n',...
                    mfilename,numel(wProcess),ceil(esttime./60));
            end;
            fprintf(1,'[%s]:Nonlinear optimization (x,y,sigma, exponent):',mfilename);
        end;
        fprintf(1,'.');drawnow;
        progress = progress + 1;
    end;

    % volume index
    vi = wProcess(ii);
    vData = data(:,ii);

    % reset tolFun: Precision of evaluation function.
    % We define RMS improvement relative to the initial raw 'no-fit' data
    % RMS. So, 1 means stop if there is less than 1% improvement on the fit:
    % searchOptions = optimset(searchOptions,'tolFun',optimget(params.analysis.fmins.options,'tolFun')./100.*rawrss);
    % optimset and optimget are a little slow so:
    searchOptions.TolFun = TolFun(ii);

    % actual fitting routine
    if searchOptions.MaxIter>0
        outParams = ...
            fmincon(@(x) rmModelSearchFit_oneGaussianNonlinear(x,vData,...
                                                            params.analysis.X,...
                                                            params.analysis.Y,...
                                                            params.analysis.allstimimages_unconvolved,...
                                                            params.analysis.Hrf,...
                                                            params.analysis.scans,...
                                                            trends),...
                     range.start(:,vi),...  % X0. Vector of 5 [x,y,sigma,exp,?]
                     [],[], ...             % A, B
                     [],[], ...             % Aeq, Beq
                     range.lower(:,vi), ... % LB
                     range.upper(:,vi),...  % UB
                     @(x) distanceCon(x,range.start(:,vi),range.step(:,vi).*expandRange),... % NONLCON
                     searchOptions); % OPTIONS
    else
        outParams = range.start(:,vi);
    end

    %[ outParams range.lower(:,vi) range.start(:,vi) range.upper(:,vi)]

    % make RF, prediction and get rss,b
    Xv = params.analysis.X-outParams(1);
    Yv = params.analysis.Y-outParams(2);
    n  = outParams(4);
    rf = exp( (Yv.*Yv + Xv.*Xv) ./ (-2.*(outParams(3).^2)) );
    pred = (params.analysis.allstimimages_unconvolved * rf).^ n;
    for scan = 1:numel(params.stim)
        inds = params.analysis.scans == scan;
        pred(inds) = filter(params.analysis.Hrf{scan}, 1, pred(inds));
    end

    X  = [pred trends];
    b    = pinv(X)*vData;
    rss  = norm(vData-X*b).^2;

    % store results only if the first beta is positive, somehow fmincon
    % outputs negative fits. If the fit is negative keep old (grid) fit. We
    % do adjust the rss, so it won't be accidentally counted as a 'good'
    % fit.
    if b(1)>0,
        model.x0(vi)         = outParams(1);
        model.y0(vi)         = outParams(2);
        model.s(vi)          = outParams(3);
        model.s_major(vi)    = outParams(3);
        model.s_minor(vi)    = outParams(3);
        model.s_theta(vi)    = 0;
        model.exponent(vi)   = outParams(4);
        model.rss(vi)        = rss;
        model.b([1 t_id],vi) = b;
    else
        % change the percent variance explained to be just under the
        % current vethresh. So it counts as a 'coarse'-fit but can still be
        % included in later 'fine'-fits
        model.rss(vi)  = (1-max((vethresh-0.01),0)).*model.rawrss(vi);
        nNegFit = nNegFit + 1;
    end
end

% end time monitor
et  = toc;
if floor(et/3600)>0,
    fprintf(1,'Done [%d hours].\n',ceil(et/3600));
else
    fprintf(1,'Done [%d minutes].\n',ceil(et/60));
end;
fprintf(1,'[%s]:Removed negative fits: %d (%.1f%%).\n',...
    mfilename,nNegFit,nNegFit./numel(wProcess).*100);
return;

%-----------------------------------
% make sure that the pRF can only be moved "step" away from original
% position "startParams" - for the one Gaussian model
function [C, Ceq]=distanceCon(x,startParams,step)
Ceq = [];
dist = x([1 2])-startParams([1 2]);
C = norm(dist) - step;
return;
%-----------------------------------


text:
function D = spm_eeg_ft_artefact_visual(S)
% Function for interactive artefact rejection using Fieldtrip
%
% Disclaimer: this code is provided as an example and is not guaranteed to work
% with data on which it was not tested. If it does not work for you, feel
% free to improve it and contribute your improvements to the MEEGtools toolbox
% in SPM (http://www.fil.ion.ucl.ac.uk/spm)
%
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_ft_artefact_visual.m 5674 2013-10-09 10:00:26Z vladimir $

[Finter,Fgraph,CmdLine] = spm('FnUIsetup', 'Fieldtrip visual artefact rejection',0);

if nargin == 0
    S = [];
end

try
    D = S.D;
catch
    D = spm_select(1, '\.mat$', 'Select EEG mat file');
    S.D = D;
end

if ischar(D)
    try
        D = spm_eeg_load(D);
    catch
        error(sprintf('Trouble reading file %s', D));
    end
end

if isequal(D.type, 'continuous')
    error('Artefact detection can only be applied to epoched data');
end

trlind = find(~badtrials(D, ':'));

data = ftraw(D, ':', ':', trlind);

data.trialinfo = [1:length(trlind)]';

cfg=[];

if isfield(S, 'method')
    cfg.method = S.method;
else
    cfg.method =  spm_input('What method?','+1', 'm', 'summary|channel|trial', strvcat('summary', 'channel', 'trial'));
end

if isfield(S, 'latency')
    cfg.latency = S.latency;
else
    cfg.latency = 1e-3*spm_input('PST ([start end] in ms):', '+1', 'r', num2str(1e3*[data.time{1}(1) data.time{1}(end)]), 2);
end

cfg.keepchannel = 'no';

if isfield(S, 'channel')
    cfg.channel = S.channel;
else
    [junk, chanind] = spm_eeg_modality_ui(D, 0, 1);
    cfg.channel = data.label(chanind);
end

data = ft_rejectvisual(cfg, data);

% Figure out based on the output of FT function what trials and channels to
% reject
trlsel = ones(1, length(trlind));
trlsel(data.trialinfo (:, 1)) = 0;

D = badtrials(D, trlind, trlsel);

badchan = setdiff(cfg.channel, data.label);
if ~isempty(badchan)
    badchanind = spm_match_str(D.chanlabels, badchan);
    D = badchannels(D, badchanind, 1);
end

if ~isfield(S, 'save') || S.save
    save(D);
end

author: spm