text:
function [A, W] = fpica(X, whiteningMatrix, dewhiteningMatrix, approach, ...
    numOfIC, g, finetune, a1, a2, myy, stabilization, ...
    epsilon, maxNumIterations, maxFinetune, initState, ...
    guess, sampleSize, displayMode, displayInterval, ...
    s_verbose);
%FPICA - Fixed point ICA. Main algorithm of FASTICA.
%
% [A, W] = fpica(whitesig, whiteningMatrix, dewhiteningMatrix, approach,
%        numOfIC, g, finetune, a1, a2, mu, stabilization, epsilon,
%        maxNumIterations, maxFinetune, initState, guess, sampleSize,
%        displayMode, displayInterval, verbose);
%
% Perform independent component analysis using Hyvarinen's fixed point
% algorithm. Outputs an estimate of the mixing matrix A and its inverse W.
%
% whitesig                              :the whitened data as row vectors
% whiteningMatrix                       :can be obtained with function whitenv
% dewhiteningMatrix                     :can be obtained with function whitenv
% approach      [ 'symm' | 'defl' ]     :the approach used (deflation or symmetric)
% numOfIC       [ 0 - Dim of whitesig ] :number of independent components estimated
% g             [ 'pow3' | 'tanh' |     :the nonlinearity used
%                 'gaus' | 'skew' ]
% finetune      [same as g + 'off']     :the nonlinearity used in finetuning.
% a1                                    :parameter for tuning 'tanh'
% a2                                    :parameter for tuning 'gaus'
% mu                                    :step size in stabilized algorithm
% stabilization [ 'on' | 'off' ]        :if mu < 1 then automatically on
% epsilon                               :stopping criterion
% maxNumIterations                      :maximum number of iterations
% maxFinetune                           :maximum number of iteretions for finetuning
% initState     [ 'rand' | 'guess' ]    :initial guess or random initial state. See below
% guess                                 :initial guess for A. Ignored if initState = 'rand'
% sampleSize    [ 0 - 1 ]               :percentage of the samples used in one iteration
% displayMode   [ 'signals' | 'basis' | :plot running estimate
%                 'filters' | 'off' ]
% displayInterval                       :number of iterations we take between plots
% verbose       [ 'on' | 'off' ]        :report progress in text format
%
% EXAMPLE
%       [E, D] = pcamat(vectors);
%       [nv, wm, dwm] = whitenv(vectors, E, D);
%       [A, W] = fpica(nv, wm, dwm);
%
%
% This function is needed by FASTICA and FASTICAG
%
%   See also FASTICA, FASTICAG, WHITENV, PCAMAT

% @(#)$Id: fpica.m,v 1.7 2005/06/16 12:52:55 jarmo Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variable for stopping the ICA calculations from the GUI
global g_FastICA_interrupt;
if isempty(g_FastICA_interrupt)
    clear global g_FastICA_interrupt;
    interruptible = 0;
else
    interruptible = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values

if nargin < 3, error('Not enough arguments!'); end
[vectorSize, numSamples] = size(X);
if nargin < 20, s_verbose = 'on'; end
if nargin < 19, displayInterval = 1; end
if nargin < 18, displayMode = 'on'; end
if nargin < 17, sampleSize = 1; end
if nargin < 16, guess = 1; end
if nargin < 15, initState = 'rand'; end
if nargin < 14, maxFinetune = 100; end
if nargin < 13, maxNumIterations = 1000; end
if nargin < 12, epsilon = 0.0001; end
if nargin < 11, stabilization = 'on'; end
if nargin < 10, myy = 1; end
if nargin < 9, a2 = 1; end
if nargin < 8, a1 = 1; end
if nargin < 7, finetune = 'off'; end
if nargin < 6, g = 'pow3'; end
if nargin < 5, numOfIC = vectorSize; end     % vectorSize = Dim
if nargin < 4, approach = 'defl'; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the data

if ~isreal(X)
    error('Input has an imaginary part.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for verbose

switch lower(s_verbose)
    case 'on'
        b_verbose = 1;
    case 'off'
        b_verbose = 0;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for approach

switch lower(approach)
    case 'symm'
        approachMode = 1;
    case 'defl'
        approachMode = 2;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''approach''\n', approach));
end
if b_verbose, fprintf('Used approach [ %s ].\n', approach); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for numOfIC

if vectorSize < numOfIC
    error('Must have numOfIC <= Dimension!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the sampleSize
if sampleSize > 1
    sampleSize = 1;
    if b_verbose
        fprintf('Warning: Setting ''sampleSize'' to 1.\n');
    end
elseif sampleSize < 1
    if (sampleSize * numSamples) < 1000
        sampleSize = min(1000/numSamples, 1);
        if b_verbose
            fprintf('Warning: Setting ''sampleSize'' to %0.3f (%d samples).\n', ...
                sampleSize, floor(sampleSize * numSamples));
        end
    end
end
if b_verbose
    if  b_verbose & (sampleSize < 1)
        fprintf('Using about %0.0f%% of the samples in random order in every step.\n',sampleSize*100);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for nonlinearity.

switch lower(g)
    case 'pow3'
        gOrig = 10;
    case 'tanh'
        gOrig = 20;
    case {'gaus', 'gauss'}
        gOrig = 30;
    case 'skew'
        gOrig = 40;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''g''\n', g));
end
if sampleSize ~= 1
    gOrig = gOrig + 2;
end
if myy ~= 1
    gOrig = gOrig + 1;
end

if b_verbose,
    fprintf('Used nonlinearity [ %s ].\n', g);
end

finetuningEnabled = 1;
switch lower(finetune)
    case 'pow3'
        gFine = 10 + 1;
    case 'tanh'
        gFine = 20 + 1;
    case {'gaus', 'gauss'}
        gFine = 30 + 1;
    case 'skew'
        gFine = 40 + 1;
    case 'off'
        if myy ~= 1
            gFine = gOrig;
        else
            gFine = gOrig + 1;
        end
        finetuningEnabled = 0;
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''finetune''\n', ...
            finetune));
end

if b_verbose & finetuningEnabled
    fprintf('Finetuning enabled (nonlinearity: [ %s ]).\n', finetune);
end

switch lower(stabilization)
    case 'on'
        stabilizationEnabled = 1;
    case 'off'
        if myy ~= 1
            stabilizationEnabled = 1;
        else
            stabilizationEnabled = 0;
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''stabilization''\n', ...
            stabilization));
end

if b_verbose & stabilizationEnabled
    fprintf('Using stabilized algorithm.\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some other parameters
myyOrig = myy;
% When we start fine-tuning we'll set myy = myyK * myy
myyK = 0.01;
% How many times do we try for convergence until we give up.
failureLimit = 5;

usedNlinearity = gOrig;
stroke = 0;
notFine = 1;
long = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for initial state.

switch lower(initState)
    case 'rand'
        initialStateMode = 0;
    case 'guess'
        if size(guess,1) ~= size(whiteningMatrix,2)
            initialStateMode = 0;
            if b_verbose
                fprintf('Warning: size of initial guess is incorrect. Using random initial guess.\n');
            end
        else
            initialStateMode = 1;
            if size(guess,2) < numOfIC
                if b_verbose
                    fprintf('Warning: initial guess only for first %d components. Using random initial guess for others.\n', size(guess,2));
                end
                guess(:, size(guess, 2) + 1:numOfIC) = ...
                    rand(vectorSize,numOfIC-size(guess,2))-.5;
            elseif size(guess,2)>numOfIC
                guess=guess(:,1:numOfIC);
                fprintf('Warning: Initial guess too large. The excess column are dropped.\n');
            end
            if b_verbose, fprintf('Using initial guess.\n'); end
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''initState''\n', initState));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the value for display mode.

switch lower(displayMode)
    case {'off', 'none'}
        usedDisplay = 0;
    case {'on', 'signals'}
        usedDisplay = 1;
        if (b_verbose & (numSamples > 10000))
            fprintf('Warning: Data vectors are very long. Plotting may take long time.\n');
        end
        if (b_verbose & (numOfIC > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    case 'basis'
        usedDisplay = 2;
        if (b_verbose & (numOfIC > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    case 'filters'
        usedDisplay = 3;
        if (b_verbose & (vectorSize > 25))
            fprintf('Warning: There are too many signals to plot. Plot may not look good.\n');
        end
    otherwise
        error(sprintf('Illegal value [ %s ] for parameter: ''displayMode''\n', displayMode));
end

% The displayInterval can't be less than 1...
if displayInterval < 1
    displayInterval = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if b_verbose, fprintf('Starting ICA calculation...\n'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SYMMETRIC APPROACH
if approachMode == 1,

    % set some parameters more...
    usedNlinearity = gOrig;
    stroke = 0;
    notFine = 1;
    long = 0;

    A = zeros(vectorSize, numOfIC);  % Dewhitened basis vectors.
    if initialStateMode == 0
        % Take random orthonormal initial vectors.
        B = orth (randn (vectorSize, numOfIC));
    elseif initialStateMode == 1
        % Use the given initial vector as the initial state
        B = whiteningMatrix * guess;
    end

    BOld = zeros(size(B));
    BOld2 = zeros(size(B));

    % This is the actual fixed-point iteration loop.
    for round = 1:maxNumIterations + 1,
        if round == maxNumIterations + 1,
            fprintf('No convergence after %d steps\n', maxNumIterations);
            fprintf('Note that the plots are probably wrong.\n');
            if ~isempty(B)
                % Symmetric orthogonalization.
                B = B * real(inv(B' * B)^(1/2));

                W = B' * whiteningMatrix;
                A = dewhiteningMatrix * B;
            else
                W = [];
                A = [];
            end
            return;
        end

        if (interruptible & g_FastICA_interrupt)
            if b_verbose
                fprintf('\n\nCalculation interrupted by the user\n');
            end
            if ~isempty(B)
                W = B' * whiteningMatrix;
                A = dewhiteningMatrix * B;
            else
                W = [];
                A = [];
            end
            return;
        end

        % Symmetric orthogonalization.
        B = B * real(inv(B' * B)^(1/2));

        % Test for termination condition. Note that we consider opposite
        % directions here as well.
        minAbsCos = min(abs(diag(B' * BOld)));
        minAbsCos2 = min(abs(diag(B' * BOld2)));

        if (1 - minAbsCos < epsilon)
            if finetuningEnabled & notFine
                if b_verbose, fprintf('Initial convergence, fine-tuning: \n'); end;
                notFine = 0;
                usedNlinearity = gFine;
                myy = myyK * myyOrig;
                BOld = zeros(size(B));
                BOld2 = zeros(size(B));

            else
                if b_verbose, fprintf('Convergence after %d steps\n', round); end

                % Calculate the de-whitened vectors.
                A = dewhiteningMatrix * B;
                break;
            end
        elseif stabilizationEnabled
            if (~stroke) & (1 - minAbsCos2 < epsilon)
                if b_verbose, fprintf('Stroke!\n'); end;
                stroke = myy;
                myy = .5*myy;
                if mod(usedNlinearity,2) == 0
                    usedNlinearity = usedNlinearity + 1;
                end
            elseif stroke
                myy = stroke;
                stroke = 0;
                if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
                    usedNlinearity = usedNlinearity - 1;
                end
            elseif (~long) & (round>maxNumIterations/2)
                if b_verbose, fprintf('Taking long (reducing step size)\n'); end;
                long = 1;
                myy = .5*myy;
                if mod(usedNlinearity,2) == 0
                    usedNlinearity = usedNlinearity + 1;
                end
            end
        end

        BOld2 = BOld;
        BOld = B;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Show the progress...
        if b_verbose
            if round == 1
                fprintf('Step no. %d\n', round);
            else
                fprintf('Step no. %d, change in value of estimate: %.3g \n', round, 1 - minAbsCos);
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Also plot the current state...
        switch usedDisplay
            case 1
                if rem(round, displayInterval) == 0,
                    % There was and may still be other displaymodes...
                    % 1D signals
                    icaplot('dispsig',(X'*B)');
                    drawnow;
                end
            case 2
                if rem(round, displayInterval) == 0,
                    % ... and now there are :-)
                    % 1D basis
                    A = dewhiteningMatrix * B;
                    icaplot('dispsig',A');
                    drawnow;
                end
            case 3
                if rem(round, displayInterval) == 0,
                    % ... and now there are :-)
                    % 1D filters
                    W = B' * whiteningMatrix;
                    icaplot('dispsig',W);
                    drawnow;
                end
            otherwise
        end

        switch usedNlinearity
            % pow3
            case 10
                B = (X * (( X' * B) .^ 3)) / numSamples - 3 * B;
            case 11
                % optimoitu - epsilonin kokoisia eroja
                % t�m� on optimoitu koodi, katso vanha koodi esim.
                % aikaisemmista versioista kuten 2.0 beta3
                Y = X' * B;
                Gpow3 = Y .^ 3;
                Beta = sum(Y .* Gpow3);
                D = diag(1 ./ (Beta - 3 * numSamples));
                B = B + myy * B * (Y' * Gpow3 - diag(Beta)) * D;
            case 12
                Xsub=X(:, getSamples(numSamples, sampleSize));
                B = (Xsub * (( Xsub' * B) .^ 3)) / size(Xsub,2) - 3 * B;
            case 13
                % Optimoitu
                Ysub=X(:, getSamples(numSamples, sampleSize))' * B;
                Gpow3 = Ysub .^ 3;
                Beta = sum(Ysub .* Gpow3);
                D = diag(1 ./ (Beta - 3 * size(Ysub', 2)));
                B = B + myy * B * (Ysub' * Gpow3 - diag(Beta)) * D;

                % tanh
            case 20
                hypTan = tanh(a1 * X' * B);
                B = X * hypTan / numSamples - ...
                    ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / numSamples * ...
                    a1;
            case 21
                % optimoitu - epsilonin kokoisia
                Y = X' * B;
                hypTan = tanh(a1 * Y);
                Beta = sum(Y .* hypTan);
                D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
                B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;
            case 22
                Xsub=X(:, getSamples(numSamples, sampleSize));
                hypTan = tanh(a1 * Xsub' * B);
                B = Xsub * hypTan / size(Xsub, 2) - ...
                    ones(size(B,1),1) * sum(1 - hypTan .^ 2) .* B / size(Xsub, 2) * a1;
            case 23
                % Optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                hypTan = tanh(a1 * Y);
                Beta = sum(Y .* hypTan);
                D = diag(1 ./ (Beta - a1 * sum(1 - hypTan .^ 2)));
                B = B + myy * B * (Y' * hypTan - diag(Beta)) * D;

                % gauss
            case 30
                U = X' * B;
                Usquared=U .^ 2;
                ex = exp(-a2 * Usquared / 2);
                gauss =  U .* ex;
                dGauss = (1 - a2 * Usquared) .*ex;
                B = X * gauss / numSamples - ...
                    ones(size(B,1),1) * sum(dGauss)...
                    .* B / numSamples ;
            case 31
                % optimoitu
                Y = X' * B;
                ex = exp(-a2 * (Y .^ 2) / 2);
                gauss = Y .* ex;
                Beta = sum(Y .* gauss);
                D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
                B = B + myy * B * (Y' * gauss - diag(Beta)) * D;
            case 32
                Xsub=X(:, getSamples(numSamples, sampleSize));
                U = Xsub' * B;
                Usquared=U .^ 2;
                ex = exp(-a2 * Usquared / 2);
                gauss =  U .* ex;
                dGauss = (1 - a2 * Usquared) .*ex;
                B = Xsub * gauss / size(Xsub,2) - ...
                    ones(size(B,1),1) * sum(dGauss)...
                    .* B / size(Xsub,2) ;
            case 33
                % Optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                ex = exp(-a2 * (Y .^ 2) / 2);
                gauss = Y .* ex;
                Beta = sum(Y .* gauss);
                D = diag(1 ./ (Beta - sum((1 - a2 * (Y .^ 2)) .* ex)));
                B = B + myy * B * (Y' * gauss - diag(Beta)) * D;

                % skew
            case 40
                B = (X * ((X' * B) .^ 2)) / numSamples;
            case 41
                % Optimoitu
                Y = X' * B;
                Gskew = Y .^ 2;
                Beta = sum(Y .* Gskew);
                D = diag(1 ./ (Beta));
                B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;
            case 42
                Xsub=X(:, getSamples(numSamples, sampleSize));
                B = (Xsub * ((Xsub' * B) .^ 2)) / size(Xsub,2);
            case 43
                % Uusi optimoitu
                Y = X(:, getSamples(numSamples, sampleSize))' * B;
                Gskew = Y .^ 2;
                Beta = sum(Y .* Gskew);
                D = diag(1 ./ (Beta));
                B = B + myy * B * (Y' * Gskew - diag(Beta)) * D;

            otherwise
                error('Code for desired nonlinearity not found!');
        end
    end

    % Calculate ICA filters.
    W = B' * whiteningMatrix;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Also plot the last one...
    switch usedDisplay
        case 1
            % There was and may still be other displaymodes...
            % 1D signals
            icaplot('dispsig',(X'*B)');
            drawnow;
        case 2
            % ... and now there are :-)
            % 1D basis
            icaplot('dispsig',A');
            drawnow;
        case 3
            % ... and now there are :-)
            % 1D filters
            icaplot('dispsig',W);
            drawnow;
        otherwise
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFLATION APPROACH
if approachMode == 2

    B = zeros(vectorSize);

    % The search for a basis vector is repeated numOfIC times.
    round = 1;

    numFailures = 0;

    while round <= numOfIC,
        myy = myyOrig;
        usedNlinearity = gOrig;
        stroke = 0;
        notFine = 1;
        long = 0;
        endFinetuning = 0;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Show the progress...
        if b_verbose, fprintf('IC %d ', round); end

        % Take a random initial vector of lenght 1 and orthogonalize it
        % with respect to the other vectors.
        if initialStateMode == 0
            w = randn (vectorSize, 1);
        elseif initialStateMode == 1
            w=whiteningMatrix*guess(:,round);
        end
        w = w - B * B' * w; %orthogonalization - B is the previous w
        w = w / norm(w);

        wOld = zeros(size(w));
        wOld2 = zeros(size(w));

        % This is the actual fixed-point iteration loop.
        %    for i = 1 : maxNumIterations + 1
        i = 1;
        gabba = 1;
        while i <= maxNumIterations + gabba
            if (usedDisplay > 0)
                drawnow;
            end
            if (interruptible & g_FastICA_interrupt)
                if b_verbose
                    fprintf('\n\nCalculation interrupted by the user\n');
                end
                return;
            end

            % Project the vector into the space orthogonal to the space
            % spanned by the earlier found basis vectors. Note that we can do
            % the projection with matrix B, since the zero entries do not
            % contribute to the projection.
            w = w - B * B' * w;
            w = w / norm(w);

            if notFine
                if i == maxNumIterations + 1
                    if b_verbose
                        fprintf('\nComponent number %d did not converge in %d iterations.\n', round, maxNumIterations);
                    end
                    round = round - 1;
                    numFailures = numFailures + 1;
                    if numFailures > failureLimit
                        if b_verbose
                            fprintf('Too many failures to converge (%d). Giving up.\n', numFailures);
                        end
                        if round == 0
                            A=[];
                            W=[];
                        end
                        return;
                    end
                    % numFailures > failurelimit
                    break;
                end
                % i == maxNumIterations + 1
            else
                % if notFine
                if i >= endFinetuning
                    wOld = w; % So the algorithm will stop on the next test...
                end
            end
            % if notFine

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Show the progress...
            if b_verbose, fprintf('.'); end;

            % Test for termination condition. Note that the algorithm has
            % converged if the direction of w and wOld is the same, this
            % is why we test the two cases.
            if norm(w - wOld) < epsilon | norm(w + wOld) < epsilon
                if finetuningEnabled & notFine
                    if b_verbose, fprintf('Initial convergence, fine-tuning: '); end;
                    notFine = 0;
                    gabba = maxFinetune;
                    wOld = zeros(size(w));
                    wOld2 = zeros(size(w));
                    usedNlinearity = gFine;
                    myy = myyK * myyOrig;

                    endFinetuning = maxFinetune + i;

                else
                    numFailures = 0;
                    % Save the vector
                    B(:, round) = w;

                    % Calculate the de-whitened vector.
                    A(:,round) = dewhiteningMatrix * w;
                    % Calculate ICA filter.
                    W(round,:) = w' * whiteningMatrix;

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Show the progress...
                    if b_verbose, fprintf('computed ( %d steps ) \n', i); end

                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    % Also plot the current state...
                    switch usedDisplay
                        case 1
                            if rem(round, displayInterval) == 0,
                                % There was and may still be other displaymodes...
                                % 1D signals
                                temp = X'*B;
                                icaplot('dispsig',temp(:,1:numOfIC)');
                                drawnow;
                            end
                        case 2
                            if rem(round, displayInterval) == 0,
                                % ... and now there are :-)
                                % 1D basis
                                icaplot('dispsig',A');
                                drawnow;
                            end
                        case 3
                            if rem(round, displayInterval) == 0,
                                % ... and now there are :-)
                                % 1D filters
                                icaplot('dispsig',W);
                                drawnow;
                            end
                    end
                    % switch usedDisplay
                    break; % IC ready - next...
                end
                %if finetuningEnabled & notFine
            elseif stabilizationEnabled
                if (~stroke) & (norm(w - wOld2) < epsilon | norm(w + wOld2) < ...
                        epsilon)
                    stroke = myy;
                    if b_verbose, fprintf('Stroke!'); end;
                    myy = .5*myy;
                    if mod(usedNlinearity,2) == 0
                        usedNlinearity = usedNlinearity + 1;
                    end
                elseif stroke
                    myy = stroke;
                    stroke = 0;
                    if (myy == 1) & (mod(usedNlinearity,2) ~= 0)
                        usedNlinearity = usedNlinearity - 1;
                    end
                elseif (notFine) & (~long) & (i > maxNumIterations / 2)
                    if b_verbose, fprintf('Taking long (reducing step size) '); end;
                    long = 1;
                    myy = .5*myy;
                    if mod(usedNlinearity,2) == 0
                        usedNlinearity = usedNlinearity + 1;
                    end
                end
            end

            wOld2 = wOld;
            wOld = w;

            switch usedNlinearity
                % pow3
                case 10
                    w = (X * ((X' * w) .^ 3)) / numSamples - 3 * w;
                case 11
                    EXGpow3 = (X * ((X' * w) .^ 3)) / numSamples;
                    Beta = w' * EXGpow3;
                    w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
                case 12
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    w = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2) - 3 * w;
                case 13
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    EXGpow3 = (Xsub * ((Xsub' * w) .^ 3)) / size(Xsub, 2);
                    Beta = w' * EXGpow3;
                    w = w - myy * (EXGpow3 - Beta * w) / (3 - Beta);
                    % tanh
                case 20
                    hypTan = tanh(a1 * X' * w);
                    w = (X * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / numSamples;
                case 21
                    hypTan = tanh(a1 * X' * w);
                    Beta = w' * X * hypTan;
                    w = w - myy * ((X * hypTan - Beta * w) / ...
                        (a1 * sum((1-hypTan .^2)') - Beta));
                case 22
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    hypTan = tanh(a1 * Xsub' * w);
                    w = (Xsub * hypTan - a1 * sum(1 - hypTan .^ 2)' * w) / size(Xsub, 2);
                case 23
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    hypTan = tanh(a1 * Xsub' * w);
                    Beta = w' * Xsub * hypTan;
                    w = w - myy * ((Xsub * hypTan - Beta * w) / ...
                        (a1 * sum((1-hypTan .^2)') - Beta));
                    % gauss
                case 30
                    % This has been split for performance reasons.
                    u = X' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    w = (X * gauss - sum(dGauss)' * w) / numSamples;
                case 31
                    u = X' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    Beta = w' * X * gauss;
                    w = w - myy * ((X * gauss - Beta * w) / ...
                        (sum(dGauss)' - Beta));
                case 32
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    u = Xsub' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    w = (Xsub * gauss - sum(dGauss)' * w) / size(Xsub, 2);
                case 33
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    u = Xsub' * w;
                    u2=u.^2;
                    ex=exp(-a2 * u2/2);
                    gauss =  u.*ex;
                    dGauss = (1 - a2 * u2) .*ex;
                    Beta = w' * Xsub * gauss;
                    w = w - myy * ((Xsub * gauss - Beta * w) / ...
                        (sum(dGauss)' - Beta));
                    % skew
                case 40
                    w = (X * ((X' * w) .^ 2)) / numSamples;
                case 41
                    EXGskew = (X * ((X' * w) .^ 2)) / numSamples;
                    Beta = w' * EXGskew;
                    w = w - myy * (EXGskew - Beta*w)/(-Beta);
                case 42
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    w = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
                case 43
                    Xsub=X(:,getSamples(numSamples, sampleSize));
                    EXGskew = (Xsub * ((Xsub' * w) .^ 2)) / size(Xsub, 2);
                    Beta = w' * EXGskew;
                    w = w - myy * (EXGskew - Beta*w)/(-Beta);

                otherwise
                    error('Code for desired nonlinearity not found!');
            end

% TEST % % % % % % % % % % % % % % % % % % % % % % % % %
interval = 1;
if mod(i,interval)==0
    modifyA =1; %for now...
    nx = 20;
    ny = 20;
    sigfix = 3.7/4; % sig of PSF...

    if exist ('modifyA','var')
        if modifyA == 1
            %                     fprintf ('Modifying matrix conevrgence of w...')
            Atemp = dewhiteningMatrix * w;
%             signum = 1;
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 signum = -1;
%             end
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 Atemp = -Atemp;
%             end
            absA = abs(Atemp);
            maxA = max(absA);

            gf = (gauss2d([nx ny],[nx/2,ny/2],sigfix,maxA));
            Aconv = (conv2(gf, reshape (Atemp, nx, ny),'same')).^5;
            Atemp = reshape(Aconv,ny*ny,1);

%             [x_mu, y_mu, sig] = fitgauss2d(reshape(absA, nx, ny),sigfix); %fit to abs...
%             Atemp = (Atemp + signum * reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1))/2;
            % Atemp = reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1);
        end
    end
    w = whiteningMatrix*Atemp;
end

% TEST % % % % % % % % % % % % % % % % % % % % % % % % %

% %%% TEST version old %%%%%%%%%
% if mod(i,5)==0
%     modifyA =1; %for now...
%     nx = 20;
%     ny = 20;
%     sigfix = 3.7/4; % sig of PSF...
%
%     if exist ('modifyA','var')
%         if modifyA == 1
%             %                     fprintf ('Modifying matrix conevrgence of w...')
%             Atemp = dewhiteningMatrix * w;
%             if abs(min(Atemp))>max(Atemp) %or sum(Atemp)<0...?
%                 Atemp = -Atemp;
%             end
%             [x_mu, y_mu, sig] = fitgauss2d(reshape(Atemp, nx, ny),sigfix);
%             maxA = max(abs(Atemp));
%             Atemp = (Atemp + reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1))/2;
%             % Atemp = reshape(gauss2d([nx ny],[x_mu,y_mu],sig,maxA),nx*ny,1);
%         end
%     end
%     w = whiteningMatrix*Atemp;
%
% %%% TEST version old %%%%%%%%%

% Normalize the new w.
            w = w / norm(w);
            i = i + 1;
        end
        round = round + 1;
    end
    if b_verbose, fprintf('Done.\n'); end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Also plot the ones that may not have been plotted.
    if (usedDisplay > 0) & (rem(round-1, displayInterval) ~= 0)
        switch usedDisplay
            case 1
                % There was and may still be other displaymodes...
                % 1D signals
                temp = X'*B;
                icaplot('dispsig',temp(:,1:numOfIC)');
                drawnow;
            case 2
                % ... and now there are :-)
                % 1D basis
                icaplot('dispsig',A');
                drawnow;
            case 3
                % ... and now there are :-)
                % 1D filters
                icaplot('dispsig',W);
                drawnow;
            otherwise
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In the end let's check the data for some security
if ~isreal(A)
    if b_verbose, fprintf('Warning: removing the imaginary part from the result.\n'); end
    A = real(A);
    W = real(W);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunction
% Calculates tanh simplier and faster than Matlab tanh.
function y=tanh(x)
y = 1 - 2 ./ (exp(2 * x) + 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Samples = getSamples(max, percentage)
Samples = find(rand(1, max) < percentage);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function error_frobenius = r8mat_is_inverse_left ( m, n, a, b )

%*****************************************************************************80
%
%% R8MAT_IS_INVERSE_LEFT determines if one matrix is the left inverse of another.
%
%  Discussion:
%
%    This routine returns the Frobenius norm of B * A - I.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns
%    of the matrix.
%
%    Input, real A(M,N), the matrix.
%
%    Input, real B(N,M), the matrix to be checked as a left inverse of A.
%
%    Output, real ERROR_FROBENIUS, the Frobenius norm
%    of the difference matrix B * A - I, which would be exactly zero
%    if B was the exact left inverse of A and computer arithmetic were exact.
%
  c(1:n,1:n) = b(1:n,1:m) * a(1:m,1:n);

  for i = 1 : n
    c(i,i) = c(i,i) - 1.0;
  end

  error_frobenius = r8mat_norm_fro ( n, n, c );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y=multimidfilter(x,m)
a=x;
for k=1 : m
    b=medfilt1(a, 5);
    a=b;
end
y=b;

author: bastamon


text:
function [ tet_index, face, step_num ] = tet_mesh_search_delaunay ( ...
  node_num, node_xyz, tet_order, tet_num, tet_node, tet_neighbor, p )

%*****************************************************************************80
%
%% TET_MESH_SEARCH_DELAUNAY searches a Delaunay tet mesh for a point.
%
%  Discussion:
%
%    The algorithm "walks" from one tetrahedron to its neighboring tetrahedron,
%    and so on, until a tetrahedron is found containing point P, or P is found
%    to be outside the convex hull.
%
%    The algorithm computes the barycentric coordinates of the point with
%    respect to the current tetrahedron.  If all 4 quantities are positive,
%    the point is contained in the tetrahedron.  If the I-th coordinate is
%    negative, then P lies on the far side of edge I, which is opposite
%    from vertex I.  This gives a hint as to where to search next.
%
%    For a Delaunay tet mesh, the search is guaranteed to terminate.
%    For other meshes, a cycle may occur.
%
%    Note the surprising fact that, even for a Delaunay tet mesh of
%    a set of nodes, the nearest node to P need not be one of the
%    vertices of the tetrahedron containing P.
%
%    The code can be called for tet meshes of any order, but only
%    the first 4 nodes in each tetrahedron are considered.  Thus, if
%    higher order tetrahedrons are used, and the extra nodes are intended
%    to give the tetrahedron a polygonal shape, these will have no effect,
%    and the results obtained here might be misleading.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2009
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Barry Joe,
%    GEOMPACK - a software package for the generation of meshes
%    using geometric algorithms,
%    Advances in Engineering Software,
%    Volume 13, pages 325-331, 1991.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the coordinates of
%    the nodes.
%
%    Input, integer TET_ORDER, the order of the tetrahedrons.
%
%    Input, intege TET_NUM, the number of tetrahedrons.
%
%    Input, integer TET_NODE(TET_ORDER,TET_NUM),
%    the nodes that make up each tetrahedron.
%
%    Input, integer TET_NEIGHBOR(4,TET_NUM), the
%    tetrahedron neighbor list.
%
%    Input, real P(3), the coordinates of a point.
%
%    Output, integer TET_INDEX, the index of the tetrahedron
%    where the search ended.  If a cycle occurred, then TET_INDEX = -1.
%
%    Output, integer FACE, indicates the position of the point P in
%    face TET_INDEX:
%    0, the interior or boundary of the tetrahedron;
%    -1, outside the convex hull of the tet mesh, past face 1;
%    -2, outside the convex hull of the tet mesh, past face 2;
%    -3, outside the convex hull of the tet mesh, past face 3.
%    -4, outside the convex hull of the tet mesh, past face 4.
%
%    Output, integer STEP_NUM, the number of steps taken.
%
  persistent tet_index_save;

  if ( length ( tet_index_save ) == 0 )
    tet_index_save = -1;
  end
%
%  If possible, start with the previous successful value of TET_INDEX.
%
  if ( tet_index_save < 1 | tet_num < tet_index_save )
    tet_index = floor ( ( tet_num + 1 ) / 2 );
  else
    tet_index = tet_index_save;
  end

  step_num = -1;
  face = 0;

  while ( 1 )

    step_num = step_num + 1;

    if ( tet_num < step_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TET_MESH_SEARCH_DELAUNAY - Fatal error!\n' );
      fprintf ( 1, '  The algorithm seems to be cycling.\n' );
      tet_index = -1;
      face = -1;
      error ( 'TET_MESH_SEARCH_DELAUNAY - Fatal error!' );
    end

    alpha = tetrahedron_barycentric ( node_xyz(1:3,tet_node(1:4,tet_index)), p );
%
%  If the barycentric coordinates are all positive, then the point
%  is inside the tetrahedron and we're done.
%
    if ( 0.0 <= alpha(1) & 0.0 <= alpha(2) & 0.0 <= alpha(3) & 0.0 <= alpha(4) )
      break
    end
%
%  At least one barycentric coordinate is negative.
%
%  If there is a negative barycentric coordinate for which there exists an
%  opposing tetrahedron neighbor closer to the point, move to that tetrahedron.
%
    if ( alpha(1) < 0.0 & 0 < tet_neighbor(1,tet_index) )
      tet_index = tet_neighbor(1,tet_index);
      continue
    elseif ( alpha(2) < 0.0 & 0 < tet_neighbor(2,tet_index) )
      tet_index = tet_neighbor(2,tet_index);
      continue
    elseif ( alpha(3) < 0.0 & 0 < tet_neighbor(3,tet_index) )
      tet_index = tet_neighbor(3,tet_index);
      continue
    elseif ( alpha(4) < 0.0 & 0 < tet_neighbor(4,tet_index) )
      tet_index = tet_neighbor(4,tet_index);
      continue
    end
%
%  All negative barycentric coordinates correspond to vertices opposite
%  faces on the convex hull.
%
%  Note the face and exit.
%
    if ( alpha(1) < 0.0 )
      face = -1;
      break
    elseif ( alpha(2) < 0.0 )
      face = -2;
      break
    elseif ( alpha(3) < 0.0 )
      face = -3;
      break
    elseif ( alpha(4) < 0.0 )
      face = -4;
      break
    end

  end

  tet_index_save = tet_index;

  return
end


text:
function [parent, post] = cs_etree (A, mode)                                %#ok
%CS_ETREE elimination tree of A or A'*A.
%   parent = cs_etree (A) returns the elimination tree of A.
%   parent = cs_etree (A,'col') returns the elimination tree of A'*A.
%   parent = cs_etree (A,'sym') is the same as cs_etree(A).
%   For the symmetric case (cs_etree(A)), only triu(A) is used.
%
%   [parent,post] = cs_etree(...) also returns a postorder of the tree.
%
%   Example:
%       Prob = UFget ('HB/bcsstk01') ; A = Prob.A ;
%       parent = cs_etree (A) ; treeplot (parent) ;
%
%   See also ETREE, TREEPLOT.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_etree mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clc;
clear all;
close all;
A=imread('cameraman.tif');
figure,imshow(uint8(A))
title('Original Image');
A=double(A);
[s1 s2]=size(A);
% bs=input('Enter the block sizes for division of the image: '); % Block Size
bs=256;

% Haar
temp=double(zeros(size(A)));
for y=1:bs:s1-bs+1
    for x=1:bs:s2-bs+1
        croppedImage = A((y:y+bs-1),(x:x+bs-1));
        t=getHaarTransform(croppedImage,bs);
        temp((y:y+bs-1),(x:x+bs-1))=t;
    end
end
figure,imshow(uint8(temp))

% Inverse Haar
temp1=double(zeros(size(A)));
for y=1:bs:s1-bs+1
    for x=1:bs:s2-bs+1
        croppedImage = temp((y:y+bs-1),(x:x+bs-1));
        t=getInvHaarTransform(croppedImage,bs);
        temp1((y:y+bs-1),(x:x+bs-1))=t;
    end
end
figure,imshow(uint8(temp1))

author: Sable

repo: mcbench-benchmarks


text:
function varargout = gzipencode(varargin)
%
% output = gzipencode(input)
%    or
% [output, info] = gzipencode(input)
%
% Compress a string or numerical array using the GZIP-compression
%
% This function depends on JVM in MATLAB or, can optionally use the ZMat
% toolbox (http://github.com/fangq/zmat)
%
% Copyright (c) 2012, Kota Yamaguchi
% URL: https://www.mathworks.com/matlabcentral/fileexchange/39526-byte-encoding-utilities
%
% Modified by: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%      input: the original data, can be a string, a numerical vector or array
%
% output:
%      output: the decompressed byte stream stored in a uint8 vector; if info is
%            given, output will restore the original data's type and dimensions
%
% examples:
%      [bytes, info]=gzipencode(eye(10));
%      orig=gzipdecode(bytes,info);
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(nargin==0)
    error('you must provide at least 1 input');
end

if(exist('zmat','file')==2 || exist('zmat','file')==3)
    [varargout{1:nargout}]=zmat(varargin{1},1,'gzip');
    return;
elseif(isoctavemesh)
    error('You must install the ZMat toolbox (http://github.com/fangq/zmat) to use this function in Octave');
end

error(javachk('jvm'));

if(ischar(varargin{1}))
    varargin{1}=uint8(varargin{1});
end

input=typecast(varargin{1}(:)','uint8');

buffer = java.io.ByteArrayOutputStream();
gzip = java.util.zip.GZIPOutputStream(buffer);
gzip.write(input, 0, numel(input));
gzip.close();

varargout{1} = typecast(buffer.toByteArray(), 'uint8')';

if(nargout>1)
    varargout{2}=struct('type',class(varargin{1}),'size',size(varargin{1}),'method','gzip','status',0);
end


text:
function out= proc_average(epo, varargin)
%PROC_AVERAGE - Classwise calculated averages
%
%Synopsis:
% EPO= proc_average(EPO, <OPT>)
% EPO= proc_average(EPO, CLASSES)
%
%Arguments:
% EPO -      data structure of epoched data
%            (can handle more than 3-dimensional data, the average is
%            calculated across the last dimension)
% OPT struct or property/value list of optional arguments:
%  .Policy - 'mean' (default), 'nanmean', or 'median'
%  .Std    - if true, standard deviation is calculated also
%  .Stats  - if true, additional statistics are calculated, including the
%            standard error of the mean, the t score, the p-value for the null
%            Hypothesis that the mean is zero, and the "signed log p-value"
%  .Classes - classes of which the average is to be calculated,
%            names of classes (strings in a cell array), or 'ALL' (default)
% 'Bonferroni' - if true, Bonferroni corrected is used to adjust p-values
%                and their logarithms
% 'Alphalevel' - if provided, a binary indicator of the significance to the
%                alpha level is returned for each feature in fv_rval.sigmask
%
% For compatibility PROC_AVERAGE can be called in the old format with CLASSES
% as second argument (which is now set via opt.Classes):
% CLASSES - classes of which the average is to be calculated,
%           names of classes (strings in a cell array), or 'ALL' (default)
%
%Returns:
% EPO     - updated data structure with fields
%  .x     - classwise means
%  .N     - vector of epochs per class across which average was calculated
%  .std   - standard deviation, if requested (opt.Std==1), format as epo.x
%  .se    - contains the standard error of the mean, if opt.Stats==1
%  .tstat     - Student t statistics of the difference, if opt.Stats==1
%  .df    - degrees of freedom of the t distribution (one sample test)
%  .p     - p value of null hypothesis that the mean is zero,
%           derived from t Statistics using two-sided test, if opt.Stats==1
%           If opt.Bonferroni==1, the p-value is multiplied by
%           epo.corrfac and cropped at 1.
%  .sgnlogp - contains the signed log10 p-value, if opt.Stats==1
%           if opt.Bonferroni==1, the p-value is multiplied by
%           epo.corrfac, cropped, and then logarithmized
%  .sigmask - binary array indicating significance at alpha level
%             opt.Alphalevel, if opt.Stats==1 and opt.Alphalevel > 0
%  .corrfac - Bonferroni correction factor (number of simultaneous tests),
%             if opt.Bonferroni==1
%  .crit    - 'significance' threshold of t statistics with respect to
%             level alpha
%
% Benjamin Blankertz
% 09-2012 stefan.haufe@tu-berlin.de
% 10-2015 Daniel Miklody

props= {  'Policy'   'mean' '!CHAR(mean nanmean median)';
          'Classes' 'ALL'   '!CHAR';
          'Std'      0      '!BOOL';
          'Stats'      0    '!BOOL';
          'Bonferroni' 0    '!BOOL';
          'Alphalevel' []   'DOUBLE'};

if nargin==0,
  out = props; return
end

misc_checkType(epo, 'STRUCT(x clab y)');
if nargin==2&&(iscellstr(varargin{1})||ischar(varargin{1}))
  opt.Classes = varargin{:};
else
  opt= opt_proplistToStruct(varargin{:});
end
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);
epo = misc_history(epo);

%% delegate a special case:
if isfield(epo, 'yUnit') && isequal(epo.yUnit, 'dB'),
  out= proc_dBAverage(epo, varargin{:});
  return;
end

%%
classes = opt.Classes;

if ~isfield(epo, 'y'),
  warning('no classes label found: calculating average across all epochs');
  nEpochs= size(epo.x, ndims(epo.x));
  epo.y= ones(1, nEpochs);
  epo.className= {'all'};
end

if isequal(opt.Classes, 'ALL'),
  classes= epo.className;
end
if ischar(classes), classes= {classes}; end
if ~iscell(classes),
  error('classes must be given cell array (or string)');
end
nClasses= length(classes);

if max(sum(epo.y,2))==1,
  warning('only one epoch per class - nothing to average');
  out= proc_selectClasses(epo, classes);
  out.N= ones(1, nClasses);
  return;
end

out= epo;
%  clInd= find(ismember(epo.className, classes));
%% the command above would not keep the order of the classes in cell 'ev'
evInd= cell(1,nClasses);
for ic= 1:nClasses,
  clInd= find(ismember(epo.className, classes{ic},'legacy'));
  evInd{ic}= find(epo.y(clInd,:));
end

sz= size(epo.x);
out.x= zeros(prod(sz(1:end-1)), nClasses);
if opt.Std,
  out.std= zeros(prod(sz(1:end-1)), nClasses);
end
if opt.Stats,
  out.se = zeros(prod(sz(1:end-1)), nClasses);
  out.p = zeros(prod(sz(1:end-1)), nClasses);
  out.tstat = zeros(prod(sz(1:end-1)), nClasses);
  out.sgnlogp = zeros(prod(sz(1:end-1)), nClasses);
end
out.y= eye(nClasses);
out.className= classes;
out.N= zeros(1, nClasses);
epo.x= reshape(epo.x, [prod(sz(1:end-1)) sz(end)]);
for ic= 1:nClasses,
  switch(lower(opt.Policy)),  %% alt: feval(opt.Policy, ...)
   case 'mean',
    out.x(:,ic)= mean(epo.x(:,evInd{ic}), 2);
   case 'nanmean',
    out.x(:,ic)= nanmean(epo.x(:,evInd{ic}), 2);
   case 'median',
    warning('median computation will be handled by proc_percentiles in the future');
    out.x(:,ic)= median(epo.x(:,evInd{ic}), 2);
   otherwise,
    error('unknown policy');
  end
  if opt.Std,
    if strcmpi(opt.Policy,'nanmean'),
      out.std(:,ic)= nanstd(epo.x(:,evInd{ic}), 0, 2);
    else
      out.std(:,ic)= std(epo.x(:,evInd{ic}), 0, 2);
    end
  end
  out.N(ic)= length(evInd{ic});
  if opt.Stats,
    if strcmpi(opt.Policy,'nanmean'),
      [H out.p(:, ic) ci stats] = ttest(epo.x(:,evInd{ic}), [], [], [], 2);
      out.se(:,ic)= stats.sd/sqrt(out.N(ic));
    else
      [H out.p(:, ic) ci stats] = ttest(epo.x(:,evInd{ic}), [], [], [], 2);
      out.se(:,ic)= stats.sd/sqrt(out.N(ic));
    end
    out.tstat(:, ic) = stats.tstat;
    out.df(ic) = stats.df(1);
    if ~isempty(opt.Alphalevel)
      out.crit(ic) = stat_calcTCrit(opt.Alphalevel, stats.df(1));
    end
  end
end

out.x= reshape(out.x, [sz(1:end-1) nClasses]);
if opt.Std,
  out.std= reshape(out.std, [sz(1:end-1) nClasses]);
end

if opt.Stats,
  out.tstat= reshape(out.tstat, [sz(1:end-1) nClasses]);
  out.se = reshape(out.se, [sz(1:end-1) nClasses]);
  out.p = reshape(out.p, [sz(1:end-1) nClasses]);
  if opt.Bonferroni
    out.corrfac = prod(sz(1:end-1));
    out.p = min(out.p*out.corrfac, 1);
  end
  out.sgnlogp = -log10(out.p).*sign(out.x);
  if ~isempty(opt.Alphalevel)
    out.alphalevel = opt.Alphalevel;
    out.sigmask = out.p < opt.Alphalevel;
  end
end

out.indexedByEpochs = {};

author: bbci


text:
function [ SimilarityMatrix ] = similarityNeighbor( x, n, range)
%SIMILARITYNEIGHBOR Summary of this function goes here
%   Detailed explanation goes here

    sz = size(x,1);
    SimilarityMatrix = eye(sz);

    i = 1:sz-n;
    SimilarityMatrix(sub2ind([sz, sz], i+n,i)) = 1;
    SimilarityMatrix(sub2ind([sz, sz], i,i+n)) = 1;

    % invalidate the illegal values from the mask (if at least one element is
    % not present in the mask set similarity to 0)
%     if(numel(mask)~=0)
%         invalidInds = sum(mask(:,range),2) < numel(range);
%
%         SimilarityMatrix(invalidInds,:) = 0;
%         SimilarityMatrix(:,invalidInds) = 0;
%     end

    DiagMask = ones(size(x, 1)) - eye(size(x,1));
    SimilarityMatrix = SimilarityMatrix .* DiagMask;
    SimilarityMatrix = SimilarityMatrix + eye(size(x, 1));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function rotatedData = rotatePoints(alignmentVector, originalData)

% rotatedData = rotatePoints(alignmentVector, originalData) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     Rotate the 'originalData' in the form of Nx2 or Nx3 about the origin by aligning the x-axis with the alignment vector
%
%       Rdata = rotatePoints([1,2,-1], [Xpts(:), Ypts(:), Zpts(:)]) - rotate the (X,Y,Z)pts in 3D with respect to the vector [1,2,-1]
%
%       Rotating using spherical components can be done by first converting using [dX,dY,dZ] = cart2sph(theta, phi, rho);  alignmentVector = [dX,dY,dZ];
%
% Example:
%   %% Rotate the point [3,4,-7] with respect to the following:
%   %%%% Original associated vector is always [1,0,0]
%   %%%% Calculate the appropriate rotation requested with respect to the x-axis.  For example, if only a rotation about the z-axis is
%   %%%% sought, alignmentVector = [2,1,0] %% Note that the z-component is zero
%   rotData = rotatePoints(alignmentVector, [3,4,-7]);
%
%     Author: Shawn Arseneau
%     Created: Feb.2, 2006
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    alignmentDim = numel(alignmentVector);
    DOF = size(originalData,2); %---- DOF = Degrees of Freedom (i.e. 2 for two dimensional and 3 for three dimensional data)

    if alignmentDim~=DOF
        error('Alignment vector does not agree with originalData dimensions');
    end
    if DOF<2 || DOF>3
        error('rotatePoints only does rotation in two or three dimensions');
    end

    if DOF==2  % 2D rotation...
        [rad_theta, rho] = cart2pol(alignmentVector(1), alignmentVector(2));
        deg_theta = -1 * rad_theta * (180/pi);
        ctheta = cosd(deg_theta);  stheta = sind(deg_theta);

        Rmatrix = [ctheta, -1.*stheta;...
                   stheta,     ctheta];
        rotatedData = originalData*Rmatrix;

    else    % 3D rotation...
        [rad_theta, rad_phi, rho] = cart2sph(alignmentVector(1), alignmentVector(2), alignmentVector(3));
        rad_theta = rad_theta * -1;
        deg_theta = rad_theta * (180/pi);
        deg_phi = rad_phi * (180/pi);
        ctheta = cosd(deg_theta);  stheta = sind(deg_theta);
        Rz = [ctheta,   -1.*stheta,     0;...
              stheta,       ctheta,     0;...
              0,                 0,     1];                  %% First rotate as per theta around the Z axis
        rotatedData = originalData*Rz;

        [rotX, rotY, rotZ] = sph2cart(-1* (rad_theta+(pi/2)), 0, 1);          %% Second rotation corresponding to phi
        rotationAxis = [rotX, rotY, rotZ];
        u = rotationAxis(:)/norm(rotationAxis);        %% Code extract from rotate.m from MATLAB
        cosPhi = cosd(deg_phi);
        sinPhi = sind(deg_phi);
        invCosPhi = 1 - cosPhi;
        x = u(1);
        y = u(2);
        z = u(3);
        Rmatrix = [cosPhi+x^2*invCosPhi        x*y*invCosPhi-z*sinPhi     x*z*invCosPhi+y*sinPhi; ...
                   x*y*invCosPhi+z*sinPhi      cosPhi+y^2*invCosPhi       y*z*invCosPhi-x*sinPhi; ...
                   x*z*invCosPhi-y*sinPhi      y*z*invCosPhi+x*sinPhi     cosPhi+z^2*invCosPhi]';

        rotatedData = rotatedData*Rmatrix;
    end


text:
function [B,L,N,P] = gp_bwboundaries(BW,CONN)
  % GP_BWBOUNDARIES Find boundaries of regions and holes in a logical image
  %
  % B = gp_bwboundaries(BW)
  % [B,L,N] = gp_bwboundaries(BW,CONN)
  %
  % Inputs:
  %   BW  h by w logical BWage
  % Outputs:
  %   B  cell array containing where each element B{b} contains a #B{c} by 2
  %     list of row and column indices of pixels in the bth boundary
  %   L  h by w BWage of ids of regions and holes (0 for "unlabeled" false
  %     pixels touching the border)
  %   N  number of regions (first N elements of B are regions, the rest are
  %     "holes")
  %
  % See also: gp_bwlabel
  %

  if nargin<2
    CONN = 8;
  end

  A4 = fd_laplacian(size(BW))~=0;
  A8 = ((A4*1*A4)>1) + A4;
  A = {[],[]};
  switch CONN
  case 4
    A{1} = A8;
    A{2} = A4;
  case 8
    A{1} = A4;
    A{2} = A8;
  end

  LN = gp_bwlabel(BW);
  LP = gp_bwlabel(~BW);
  % Combine
  L = LN + (LN==0).*(LP+max(LN(:)));

  % BFS from outside
  O = reshape(full(sparse([2:size(BW,1)-1 1:size(BW,1):numel(BW)-size(BW,1) size(BW,1):size(BW,1):numel(BW) numel(BW)-size(BW,1)+1:numel(BW)-1],1,1)),size(BW));

  K = {O&BW,false(size(BW))};
  outside = [];
  inside = [];

  outside = union(outside,setdiff(unique(L.*O.*~BW),0));

  while true
    G = ismember(L,union(inside,outside));
    %BWshow(G);
    %pause
    K{1} = K{1} | (reshape((A{1}*1*G(:))~=0,size(BW)) & BW & ~G);
    %BWshow(matrixnormalize(K{1}-K{2}))
    %pause
    inside = union(inside,setdiff(unique(L.*(K{1}|K{2})),0));
    if isempty(setdiff(L,union(inside,outside)))
      break;
    end
    G = ismember(L,union(inside,outside));
    K{2} = K{2} | (reshape((A{2}*1*G(:))~=0,size(BW)) & ~G);
    outside = setdiff(union(outside,L.*((~BW & reshape((A{1}*1*G(:))~=0,size(BW))))),0);
    %BWshow(matrixnormalize(K{1}-K{2}))
    %pause
    if isempty(setdiff(L,union(inside,outside)))
      break;
    end
  end

  B = {{},{}};
  for c = 1:max(L(:))
    for e = 1:2
      BWc = K{e} & L==c;
      [Lc,ncc] = gp_bwlabel(BWc,8);
      for cc = 1:ncc
        BWcc = Lc == cc;
        % Super inefficient search
        Acc = A{3-e};
        Acc(~BWcc,:) = 0;
        Acc(:,~BWcc) = 0;
        Acc = Acc - diag(diag(Acc));
        [s,t] = find(Acc,1,'first');
        P = [];
        p = s;
        while sum(Acc(:))>0
          P = [P;p];
          Acc(p,:) = 0;
          [~,p] = max(Acc(:,p));
        end
        [LI,LJ] = ind2sub(size(BW),P);
        if ~isempty(LI)
          %
          B{e} = {B{e}{:} [LI([1:end 1]) LJ([1:end 1])]}';
        end
      end
    end
  end
  N = numel(B{1});
  B = {B{1}{:} B{2}{:}}';

  L = L.*~ismember(L,L.*(O&~BW));
  [~,~,L] = unique(L);
  L = reshape(L,size(BW));

end

author: alecjacobson

repo: gptoolbox


text:
function [uniformScan, firstSliceZValue] = scanUniformize(scanStruct, scanArray, scanInfo, dim3Spacing, tMin, tMax, keepSlice, optS, hBar)
%"scanUniformize"
%   Interpolates the input CT data values to a grid of values
%   uniform in the z direction (and therefore in each direction).
%   The scan in the passed scanStruct is used.
%
%   Format: Scan array input in (rows,cols,k) and output in
%   (rows, cols, k).
%
%   keepSlice is either 'last' or 'first', depending on which
%   slice should be kept. The uniform z-value calculations are
%   made beginning from that slice. Default value is 'first'.
%
%   note: this program updates a waitbar between the values of
%         tMin and tMax. Suggested use is to create a waitbar
%         before calling this function and to close it afterwards.
%
%Latest modifications:
% 14 Aug 02, V H Clark, creation
% Jan 02, JOD, added 8-bit support for reduced memory.
% 23 Apr 03, JOD, added back 16-bit support and fixed bug.
% 09 Apr 03, JOD, added hBar to parameter list.
% 18 Feb 05, JRA, Added scanStruct parameter for multiscan support.
%
%Usage:
%   function [uniformScan, firstSliceZValue] = scanUniformize(scanStruct, scanArray, scanInfo, dim3Spacing, tMin, tMax, keepSlice, optS, hBar)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if nargin < 6
  keepSlice = 'first';
end

tDelta = tMax - tMin;

%find desired z values, zi:

firstZValue = scanInfo(1).zValue;
lastZValue  = scanInfo(end).zValue;

if strcmp(keepSlice,'first')
  ziValues = [firstZValue : dim3Spacing : lastZValue];
elseif strcmp(keepSlice,'last')
  ziValuesRev = [lastZValue : -dim3Spacing : firstZValue];
  ziValues = ziValuesRev(end:-1:1);
end
firstSliceZValue = ziValues(1);

zValues  = [scanInfo(:).zValue];
helperIndex = 1;

xSize = scanInfo(1).sizeOfDimension2; %doesn't matter which slice is used as the scanInfo index -- they should all be the same.
ySize = scanInfo(1).sizeOfDimension1;
numOfElts = xSize*ySize;

xUnits = scanInfo(1).grid1Units;
yUnits = scanInfo(1).grid2Units;

xMax = xUnits*(xSize - 1);
yMax = yUnits*(ySize - 1);

if strcmpi(optS.uniformizedDataType,'uint8')
  uniformScan = zeros(ySize, xSize, length(ziValues),'uint8');
elseif strcmpi(optS.uniformizedDataType,'uint16') %assume uint16
  uniformScan = zeros(ySize, xSize, length(ziValues),'uint16');
end

if strcmpi(optS.uniformizedDataType,'uint8')
  CTMin = scanStruct.uniformScanInfo.minCTValue;
  CTMax = scanStruct.uniformScanInfo.maxCTValue;
  CTScale = 255 / (CTMax - CTMin);
elseif strcmpi(optS.uniformizedDataType,'uint16')
  CTMin = scanStruct.uniformScanInfo.minCTValue;
  CTMax = scanStruct.uniformScanInfo.maxCTValue;
  CTScale = 65535 / (CTMax - CTMin);
end

waitbar(tMin, hBar);
for k = 1 : length(ziValues)
  zi = ziValues(k);

  %Find nearest slices to zi value.
  [cranialSlice, caudalSlice, helperIndex] = findSurroundingSlices(zi, zValues, helperIndex);

  sliceWidth = abs(zValues(cranialSlice) - zValues(caudalSlice));
  if sliceWidth ~= 0
      cranialWeight = 1 - abs(zValues(cranialSlice) - zi) / sliceWidth;
      caudalWeight  = 1 - abs(zValues(caudalSlice) - zi) / sliceWidth;
  else
      cranialWeight = 1;
      caudalWeight  = 0;
  end

  %Use simple linear interpolation between the two slices, by weight.
  interpSlice = double(scanArray(:,:,cranialSlice)) * cranialWeight + double(scanArray(:,:,caudalSlice)) * caudalWeight;
  if strcmpi(optS.uniformizedDataType,'uint8')
      uniformScan(:,:,k) = uint8((interpSlice - CTMin) * CTScale);
  elseif  strcmpi(optS.uniformizedDataType,'uint16')
      uniformScan(:,:,k) = uint16((interpSlice - CTMin) * CTScale);
  end
  waitbar(tMin + (k/length(ziValues))*tDelta, hBar);
end

return

function [cranialSlice, caudalSlice, hi] = findSurroundingSlices(zi, zValues, hi)
% Find the cranial and caudal slice numbers.
% assumes that z values increase with index number.
% hi is a helper index that must be equal to or lower than the cranial slice number desired.
%    (it makes this function faster)

done = 0;
while ~done
  if (zi < zValues(hi))
    caudalSlice  = hi;
    cranialSlice = hi-1;
    hi = cranialSlice;
    done = 1;
  elseif (zValues(hi) == zi)
    cranialSlice = hi;
    caudalSlice = hi;
    done = 1;
  else
    hi = hi+1;
  end
end

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [aVals,aJacob,aHess,papt]=aSpiralSimp(xPoints,numRetDims)
%%ASPIRALSIMP The drift function for a non-ballistic spiraling target
%         motion model in 3 dimensions, formulated in a manner such that it
%         is simple to make a spiraling target follow a nominal trajectory.
%         This model is probably bad to design into a target tracking
%         algorithm, but might be good for designing spiraling trajectories
%         to use to test target tracking algorithms.
%
%INPUTS: xPoints The target state at time t. It consists of position
%                (3 elements), instantaneous velocity (3 elements), the
%                velocity vector (3 elements) of the overall direction of
%                motion of the spiraling model (ground speed) and the
%                spiral rate of the target. Thus xDim=10. If x is an
%                xDim X numStates matrix, then the spiraling model is
%                evaluated for all of the state vectors. The last 4
%                elements in xPoints are taken to be constants and thus
%                aVal entries will be 0 for them. If only aVals is
%                requested on the output, then xPoints can be a matrix of
%                points.
%     numRetDims If numRetDims=6, then the last 4 elements in xPoints (the
%                velocity vector of the overall spiraling model and the
%                spiral rate) are taken to be constants and the returned
%                aVals Vec is 6X1. Otherwise, numRetDims can be 10 and
%                aVals is 10X1. The default if omitted or an empty matrix
%                is passed is 6.
%
%OUTPUTS: aVals The 6X1 (or 10X1 depending on numRetDims) flat-Earth time-
%               derivative of the state. If xPoints was a matrix of N
%               points, then this will be 6XN (or 10XN).
%        aJacob The 6X6 (or 10X10) matrix of partial derivatives of aVals
%               such that aJacob(:,k) is the partial derivative of
%               aVals(:,k) with respect to
%               xPoints(k). xPoints can't be a matrix if this output is
%               requested.
%         aHess The 6X6X6 (or 10X10X10) matrix of second derivatives of
%               aVals such that aHess(:,k1,k2) is the second partial
%               derivative of aVals with respect to xPoints(k1) and
%               xPoints(k2).
%          papt The 6X1 or 10X1 derivative with resect to time of aVals.
%               This is all zeros, because the model is time invariant.
%
%A derivation of the simplified flat-Earth spiraling dynamic model is given
%in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "Simulating aerial targets in 3D accounting for the
%    Earth's curvature," Journal of Advances in Information Fusion, vol.
%    10, no. 1, Jun. 2015.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<3||isempty(numRetDims))
        numRetDims=6;
    end

    if(numRetDims<6||numRetDims>10)
       error('numRetDims must be between 6 and 10.')
    end

    numPoints=size(xPoints,2);
    aVals=zeros(numRetDims,numPoints);

    vl=xPoints(7:9,:);
    omega=xPoints(10,:);

    vlMag2=sum(vl.*vl,1);
    vlMag=sqrt(vlMag2);

    aVals(1:3,:)=xPoints(4:6,:);
    vs=xPoints(4:6,:)-vl;
    ul=bsxfun(@rdivide,vl,vlMag);
    vsDot=omega.*bsxfun(@cross,ul,vs);
    aVals(4:6,:)=vsDot;

    if(nargout>1)
        rdotx=xPoints(4);
        rdoty=xPoints(5);
        rdotz=xPoints(6);

        vlx=vl(1);
        vly=vl(2);
        vlz=vl(3);

        vlx2=vlx*vlx;
        vly2=vly*vly;
        vlz2=vlz*vlz;

        vlMag2=vlMag*vlMag;
        vlMag3=vlMag2*vlMag;

        pXrdy=-(omega*vlz)/vlMag;
        pXrdz=(omega*vly)/vlMag;

        pYrdx=(omega*vlz)/vlMag;
        pYrdz=-((omega*vlx)/vlMag);

        pZrdx=-((omega*vly)/vlMag);
        pZrdy=(omega*vlx)/vlMag;

        if(numRetDims==6)
            aJacob=[0,0,0,1,        0,      0;
                    0,0,0,0,        1,      0;
                    0,0,0,0,        0,      1;
                    0,0,0,0,        pXrdy,  pXrdz;
                    0,0,0,pYrdx,    0,      pYrdz;
                    0,0,0,pZrdx,    pZrdy,  0];
        else
            pXvx=(omega*vlx*(-rdotz*vly+rdoty*vlz))/vlMag3;
            pXvy=(omega*(rdoty*vly*vlz+rdotz*(vlx2+vlz2)))/vlMag3;
            pXvz=-((omega*(rdoty*(vlx2+vly2)+rdotz*vly*vlz))/vlMag3);

            pXomega=(rdotz*vly-rdoty*vlz)/vlMag;

            pYvx=-((omega*(rdotx*vlx*vlz+rdotz*(vly2+vlz2)))/vlMag3);
            pYvy=(omega*vly*(rdotz*vlx-rdotx*vlz))/vlMag3;
            pYvz=(omega*(rdotx*(vlx2+vly2)+rdotz*vlx*vlz))/vlMag3;
            pYomega=(-rdotz*vlx+rdotx*vlz)/vlMag;

            pZvx=(omega*(rdotx*vlx*vly+rdoty*(vly2+vlz2)))/vlMag3;
            pZvy=-((omega*(rdoty*vlx*vly+rdotx*(vlx2+vlz2)))/vlMag3);
            pZvz=(omega*(-rdoty*vlx+rdotx*vly)*vlz)/vlMag3;
            pZomega=(rdoty*vlx-rdotx*vly)/vlMag;

            aJacob=[0,0,0,1,        0,      0,      0,      0,      0,      0;
                    0,0,0,0,        1,      0,      0,      0,      0,      0;
                    0,0,0,0,        0,      1,      0,      0,      0,      0;
                    0,0,0,0,        pXrdy,  pXrdz,  pXvx,   pXvy,   pXvz,   pXomega;
                    0,0,0,pYrdx,    0,      pYrdz,  pYvx,   pYvy,   pYvz,   pYomega;
                    0,0,0,pZrdx,    pZrdy,  0,      pZvx,   pZvy,   pZvz,   pZomega;
                    zeros(4,10)];
        end

        if(nargout>2)
            aHess=zeros(numRetDims,numRetDims,numRetDims);

            vlMag5=vlMag3*vlMag2;

            %Second deriavtives are all 0 if only
            if(numRetDims>6)
                pXrdyvx=(omega*vlx*vlz)/vlMag3;
                pXrdyvy=(omega*vly*vlz)/vlMag3;
                pXrdyvz=-((omega*(vlx2+vly2))/vlMag3);
                pXrdyOmega=-(vlz/vlMag);

                pXrdzvx=-((omega*vlx*vly)/vlMag3);
                pXrdzvy=(omega*(vlx2+vlz2))/vlMag3;
                pXrdzvz=-((omega*vly*vlz)/vlMag3);
                pXrdzOmega=vly/vlMag;

                pXvxrdx=0;
                pXvxrdy=pXrdyvx;
                pXvxrdz=pXrdzvx;
                pXvxvx=-((omega*(rdotz*vly-rdoty*vlz)*(-2*vlx2+vly2+vlz2))/vlMag5);
                pXvxvy=-((omega*vlx*(3*rdoty*vly*vlz+rdotz*(vlx2-2*vly2+vlz2)))/vlMag5);
                pXvxvz=(omega*vlx*(3*rdotz*vly*vlz+rdoty*(vlx2+vly2-2*vlz2)))/vlMag5;
                pXvxOmega=(-rdotz*vlx*vly+rdoty*vlx*vlz)/vlMag3;

                pXvyrdx=0;
                pXvyrdy=pXrdyvy;
                pXvyrdz=pXrdzvy;
                pXvyvx=pXvxvy;
                pXvyvy=(omega*(-3*rdotz*vly*(vlx2+vlz2)+rdoty*vlz*(vlx2-2*vly2+vlz2)))/vlMag5;
                pXvyvz=(omega*(rdoty*vly*(vlx2+vly2-2*vlz2)-rdotz*vlz*(vlx2-2*vly2+vlz2)))/vlMag5;
                pXvyOmega=(rdoty*vly*vlz+rdotz*(vlx2+vlz2))/vlMag3;

                pXvzrdx=0;
                pXvzrdy=pXrdyvz;
                pXvzrdz=pXrdzvz;
                pXvzvx=pXvxvz;
                pXvzvy=pXvyvz;
                pXvzvz=(omega*(3*rdoty*(vlx2+vly2)*vlz-rdotz*vly*(vlx2+vly2-2*vlz2)))/vlMag5;
                pXvzOmega=(-rdoty*(vlx2+vly2)-rdotz*vly*vlz)/vlMag3;

                pXOmegardx=0;
                pXOmegardy=pXrdyOmega;
                pXOmegardz=pXrdzOmega;
                pXOmegavx=pXvxOmega;
                pXOmegavy=pXvyOmega;
                pXOmegavz=pXvzOmega;
                pXOmegaOmega=0;

                pYrdxvx=-((omega*vlx*vlz)/vlMag3);
                pYrdxvy=-((omega*vly*vlz)/vlMag3);
                pYrdxvz=(omega*(vlx2+vly2))/vlMag3;
                pYrdxOmega=vlz/vlMag;

                pYrdzvx=-((omega*(vly2+vlz2))/vlMag3);
                pYrdzvy=(omega*vlx*vly)/vlMag3;
                pYrdzvz=(omega*vlx*vlz)/vlMag3;
                pYrdzOmega=-(vlx/vlMag);

                pYvxrdx=pYrdxvx;
                pYvxrdy=0;
                pYvxrdz=pYrdzvx;
                pYvxvx=(omega*(3*rdotz*vlx*(vly2+vlz2)-rdotx*vlz*(-2*vlx2+vly2+vlz2)))/vlMag5;
                pYvxvy=(omega*vly*(3*rdotx*vlx*vlz+rdotz*(-2*vlx2+vly2+vlz2)))/vlMag5;
                pYvxvz=(omega*(-rdotx*vlx*(vlx2+vly2-2*vlz2)+rdotz*vlz*(-2*vlx2+vly2+vlz2)))/vlMag5;
                pYvxOmega=(-rdotx*vlx*vlz-rdotz*(vly2+vlz2))/vlMag3;

                pYvyrdx=pYrdxvy;
                pYvyrdy=0;
                pYvyrdz=pYrdzvy;
                pYvyvx=pYvxvy;
                pYvyvy=(omega*(rdotz*vlx-rdotx*vlz)*(vlx2-2*vly2+vlz2))/vlMag5;
                pYvyvz=-((omega*vly*(3*rdotz*vlx*vlz+rdotx*(vlx2+vly2-2*vlz2)))/vlMag5);
                pYvyOmega=(vly*(rdotz*vlx-rdotx*vlz))/vlMag3;

                pYvzrdx=pYrdxvz;
                pYvzrdy=0;
                pYvzrdz=pYrdzvz;
                pYvzvx=pYvxvz;
                pYvzvy=pYvyvz;
                pYvzvz=(omega*(-3*rdotx*(vlx2+vly2)*vlz+rdotz*vlx*(vlx2+vly2-2*vlz2)))/vlMag5;
                pYvzOmega=(rdotx*(vlx2+vly2)+rdotz*vlx*vlz)/vlMag3;

                pYOmegardx=pYrdxOmega;
                pYOmegardy=0;
                pYOmegardz=pYrdzOmega;
                pYOmegavx=pYvxOmega;
                pYOmegavy=pYvyOmega;
                pYOmegavz=pYvzOmega;
                pYOmegaOmega=0;

                pZrdxvx=(omega*vlx*vly)/vlMag3;
                pZrdxvy=-((omega*(vlx2+vlz2))/vlMag3);
                pZrdxvz=(omega*vly*vlz)/vlMag3;
                pZrdxOmega=-(vly/vlMag);

                pZrdyvx=(omega*(vly2+vlz2))/vlMag3;
                pZrdyvy=-((omega*vlx*vly)/vlMag3);
                pZrdyvz=-((omega*vlx*vlz)/vlMag3);
                pZrdyOmega=vlx/vlMag;

                pZvxrdx=pZrdxvx;
                pZvxrdy=pZrdyvx;
                pZvxrdz=0;
                pZvxvx=(omega*(-3*rdoty*vlx*(vly2+vlz2)+rdotx*vly*(-2*vlx2+vly2+vlz2)))/vlMag5;
                pZvxvy=(omega*(rdotx*vlx*(vlx2-2*vly2+vlz2)-rdoty*vly*(-2*vlx2+vly2+vlz2)))/vlMag5;
                pZvxvz=-((omega*vlz*(3*rdotx*vlx*vly+rdoty*(-2*vlx2+vly2+vlz2)))/vlMag5);
                pZvxOmega=(rdotx*vlx*vly+rdoty*(vly2+vlz2))/vlMag3;

                pZvyrdx=pZrdxvy;
                pZvyrdy=pZrdyvy;
                pZvyrdz=0;
                pZvyvx=pZvxvy;
                pZvyvy=(omega*(3*rdotx*vly*(vlx2+vlz2)-rdoty*vlx*(vlx2-2*vly2+vlz2)))/vlMag5;
                pZvyvz=(omega*vlz*(3*rdoty*vlx*vly+rdotx*(vlx2-2*vly2+vlz2)))/vlMag5;
                pZvyOmega=(-rdoty*vlx*vly-rdotx*(vlx2+vlz2))/vlMag3;

                pZvzrdx=pZrdxvz;
                pZvzrdy=pZrdyvz;
                pZvzrdz=0;
                pZvzvx=pZvxvz;
                pZvzvy=pZvyvz;
                pZvzvz=-((omega*(rdoty*vlx-rdotx*vly)*(vlx2+vly2-2*vlz2))/vlMag5);
                pZvzOmega=(-rdoty*vlx*vlz+rdotx*vly*vlz)/vlMag3;

                pZOmegardx=pZrdxOmega;
                pZOmegardy=pZrdyOmega;
                pZOmegardz=0;
                pZOmegavx=pZvxOmega;
                pZOmegavy=pZvyOmega;
                pZOmegavz=pZvzOmega;
                pZOmegaOmega=0;

                aHess(:,:,4)=[zeros(3,10);
                              0,0,0,0,0,0,  pXvxrdx,   pXvyrdx,   pXvzrdx,   pXOmegardx;
                              0,0,0,0,0,0,  pYvxrdx,   pYvyrdx,   pYvzrdx,   pYOmegardx;
                              0,0,0,0,0,0,  pZvxrdx,   pZvyrdx,   pZvzrdx,   pZOmegardx;
                              zeros(4,10)];
                aHess(:,:,5)=[zeros(3,10);
                              0,0,0,0,0,0,  pXvxrdy,   pXvyrdy,   pXvzrdy,   pXOmegardy;
                              0,0,0,0,0,0,  pYvxrdy,   pYvyrdy,   pYvzrdy,   pYOmegardy;
                              0,0,0,0,0,0,  pZvxrdy,   pZvyrdy,   pZvzrdy,   pZOmegardy;
                              zeros(4,10)];
                aHess(:,:,6)=[zeros(3,10);
                              0,0,0,0,0,0,  pXvxrdz,   pXvyrdz,   pXvzrdz,   pXOmegardz;
                              0,0,0,0,0,0,  pYvxrdz,   pYvyrdz,   pYvzrdz,   pYOmegardz;
                              0,0,0,0,0,0,  pZvxrdz,   pZvyrdz,   pZvzrdz,   pZOmegardz;
                              zeros(4,10)];
                aHess(:,:,7)=[zeros(3,10);
                              0,0,0,0,          pXrdyvx,  pXrdzvx,  pXvxvx,   pXvyvx,   pXvzvx,   pXOmegavx;
                              0,0,0,pYrdxvx,    0,        pYrdzvx,  pYvxvx,   pYvyvx,   pYvzvx,   pYOmegavx;
                              0,0,0,pZrdxvx,    pZrdyvx,  0,        pZvxvx,   pZvyvx,   pZvzvx,   pZOmegavx;
                              zeros(4,10)];
                aHess(:,:,8)=[zeros(3,10);
                              0,0,0,0,          pXrdyvy,  pXrdzvy,  pXvxvy,   pXvyvy,   pXvzvy,   pXOmegavy;
                              0,0,0,pYrdxvy,    0,        pYrdzvy,  pYvxvy,   pYvyvy,   pYvzvy,   pYOmegavy;
                              0,0,0,pZrdxvy,    pZrdyvy,  0,        pZvxvy,   pZvyvy,   pZvzvy,   pZOmegavy;
                              zeros(4,10)];
                aHess(:,:,9)=[zeros(3,10);
                              0,0,0,0,          pXrdyvz,  pXrdzvz,  pXvxvz,   pXvyvz,   pXvzvz,   pXOmegavz;
                              0,0,0,pYrdxvz,    0,        pYrdzvz,  pYvxvz,   pYvyvz,   pYvzvz,   pYOmegavz;
                              0,0,0,pZrdxvz,    pZrdyvz,  0,        pZvxvz,   pZvyvz,   pZvzvz,   pZOmegavz;
                              zeros(4,10)];
               aHess(:,:,10)=[zeros(3,10);
                              0,0,0,0,             pXrdyOmega,  pXrdzOmega,  pXvxOmega,   pXvyOmega,   pXvzOmega,   pXOmegaOmega;
                              0,0,0,pYrdxOmega,    0,           pYrdzOmega,  pYvxOmega,   pYvyOmega,   pYvzOmega,   pYOmegaOmega;
                              0,0,0,pZrdxOmega,    pZrdyOmega,  0,           pZvxOmega,   pZvyOmega,   pZvzOmega,   pZOmegaOmega;
                              zeros(4,10)];
            end
            if(nargout>3)
                papt=zeros(10,1);
            end
        end
    end
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 fx1 = p09_fx1 ( x )

%*****************************************************************************80
%
%% P09_FX1 evaluates the derivative of the function for problem 9.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX1, the first derivative of the function at X.
%
  x2 = x - 6.25;

  if ( x2 < - 0.25 )
    fx1 = 0.75;
  elseif ( x2 < 0.25 )
    fx1 = 2.0;
  else
    fx1 = 0.75;
  end

  return
end


text:
function xyzf_header_write ( output_filename, output_unit, point_num, ...
  face_num, face_data_num )

%*****************************************************************************80
%
%% XYZF_HEADER_WRITE writes the header of an XYL file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the output file.
%
%    Input, integer OUTPUT_UNIT, the output file unit number.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, integer FACE_DATA_NUM, the number of face items.
%
  fprintf ( output_unit, '#  %s\n', output_filename );
  fprintf ( output_unit, '#  created by xyz_io::xyzf_header_write.m\n' );
  fprintf ( output_unit, '#\n');
  fprintf ( output_unit, '#  Number of points =     %d\n', point_num );
  fprintf ( output_unit, '#  Number of faces =      %d\n', face_num );
  fprintf ( output_unit, '#  Number of face items = %d\n', face_data_num );
  fprintf ( output_unit, '#\n' );

  return
end


text:
  function Cdiff1_test1(C)
%|function Cdiff1_test1(C)
%| private helper for Cdiff1_test
%| tests abs and power and adjoints

Cf = C(:,:);

% abs
Ca = abs(C); Caf = Ca(:,:);
jf_equal(Caf, abs(Cf))

% squared
C2 = C.^2; C2f = C2(:,:);
jf_equal(C2f, Cf.^2)

test_adjoint(C);
test_adjoint(Ca);
test_adjoint(C2);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [tfVMVariates] = vmrand(fMu, fKappa, varargin)

% vmrand - FUNCTION Draw random variates from the Von Mises circular distribution
%
% Usage: [tfVMVariates] = vmrand(fMu, fKappa)
%        [tfVMVariates] = vmrand(..., M)
%        [tfVMVariates] = vmrand(..., M, N, P, ...)
%        [tfVMVariates] = vmrand(..., [M N P ...])
%
% This function uses an envelope-rejection method based on a wrapped Cauchy
% distribution to draw random variates from an arbitrary Von Mises
% distribution, first proposed in [1].
%
% 'fMu' and 'fKappa' are the mean and variance parameter of the Von Mises
% distribution over [-pi, pi).  'tVMVariates' will be a tensor containing
% random variates drawn from the defined distribution.  If 'fMu' and
% 'fKappa' are non-scalar, then they must be the same size.  In this case
% 'tVMVariates' will be the same size.  If 'fMu' and 'fKappa' are scalar,
% then the number of variates returned can be specified as extra arguments.
%
% If only single dimension 'M' is provided, then the return argument
% 'tfVMVariates' will be M x M.
%
% This implementation is vectorised, and requires O(7.5*N) space.
%
% References:
% [1] D. J. Best and N. I. Fisher, 1979. "Efficient Simulation of the von
% Mises Distribution", Journal of the Royal Statistical Society. Series C
% (Applied Statistics), Vol. 28, No. 2, pp. 152-157.

% Author: Dylan Muir <muir@hifo.uzh.ch>
% Created: 19th June, 2012

% -- Check arguments

if (nargin < 2)
   help vmrand;
   error('VMRAND:Usage', '*** vmrand: Incorrect usage');
end

% -- Check sizes

vbScalarArgs = [isscalar(fMu) isscalar(fKappa)];

if (nnz(vbScalarArgs) == 1)
   if vbScalarArgs(1)
      % - fMu is scalar
      fMu = repmat(fMu, size(fKappa));
   else
      % - fKappa is scalar
      fKappa = repmat(fKappa, size(fMu));
   end

   % - Set return sizes
   vnTensorSize = size(fMu);

elseif (nnz(vbScalarArgs == 0))
   % - Two non-scalar arguments
   if (~isequal(size(fMu), size(fKappa)))
      error('VMRAND:UnequalSizes', ...
         '*** vmrand: ''fMu'' and ''fKappa must be the same size.');
   else
      vnTensorSize = size(fMu);
   end

elseif (~isempty(varargin))
   % - Get argument sizes from varargin (be forgiving)
   varargin = cellfun(@(c)(reshape(c, 1, [])), varargin, 'UniformOutput', false);
   vnTensorSize = [varargin{:}];

   % - Take Matlab semantics to make square matrices
   if (isscalar(vnTensorSize))
      vnTensorSize = vnTensorSize([1 1]);
   end

   fKappa = repmat(fKappa, vnTensorSize);

else
   % - Return a scalar variate
   vnTensorSize = [1 1];
end

% -- Check values

if (any(fKappa < 0))
   error('VMRAND:InvalidArguments', ...
      '*** vmrand: ''fKappa'' must be >= 0');
end

if (any(vnTensorSize < 1))
   error('VMRAND:InvalidArguments', ...
      '*** vmrand: Tensor size dimensions must be positive');
end

% -- Preallocate data

tfVMVariates = nan(vnTensorSize);
tfZ = nan(vnTensorSize);
tfF = nan(vnTensorSize);
tfC = nan(vnTensorSize);

% -- Short-cut fKappa == 0

tbUniform = fKappa == 0;
tfVMVariates(tbUniform) = rand(nnz(tbUniform), 1) * 2*pi - pi;
tbDraw = ~tbUniform;
tbAccept = tbUniform;

% -- Pre-compute what we can

if (all(vbScalarArgs))
   tfTau = sqrt(4 .* fKappa(1).^2 + 1) + 1;
   tfRho = (tfTau - sqrt(2 .* tfTau)) ./ (2 .* fKappa(1));
   tfR = repmat((1 + tfRho.^2) ./ (2 .* tfRho), vnTensorSize);

else
   tfTau = sqrt(4 .* fKappa.^2 + 1) + 1;
   tfRho = (tfTau - sqrt(2 .* tfTau)) ./ (2 .* fKappa);
   tfR = (1 + tfRho.^2) ./ (2 .* tfRho);
   clear tfTau tfRho;   % - To save some space
end

% -- Draw random variates

while (nnz(tbDraw > 0))
   % - Draw partial variates and estimate wrapped Cauchy distribution envelope
   nNumToDraw = nnz(tbDraw);
   tfZ(tbDraw) = cos(pi .* rand(nNumToDraw, 1));
   tfF(tbDraw) = (1 + tfR(tbDraw) .* tfZ(tbDraw)) ./ (tfR(tbDraw) + tfZ(tbDraw));
   tfC(tbDraw) = fKappa(tbDraw) .* (tfR(tbDraw) - tfF(tbDraw));

   % - Filter variates
   vfRand2 = rand(nNumToDraw, 1);
   tbAccept(tbDraw) = (tfC(tbDraw) .* (2 - tfC(tbDraw)) - vfRand2) > 0;
   vbRecheck = ~tbAccept(tbDraw);
   if (any(vbRecheck))
      tbAccept(tbDraw & ~tbAccept) = (log(tfC(tbDraw & ~tbAccept) ./ vfRand2(vbRecheck)) + 1 - tfC(tbDraw & ~tbAccept)) >= 0;
   end

   % - Construct final variates
   nNumToAccept = nnz(tbDraw & tbAccept);
   tfVMVariates(tbDraw & tbAccept) = sign(rand(nNumToAccept, 1) - 0.5) .* acos(tfF(tbDraw & tbAccept));

   % - Mark as being accepted (don't need to try again)
   tbDraw(tbAccept) = false;
end

% -- Shift variates to mean

tfVMVariates = tfVMVariates + fMu;
tfVMVariates(tfVMVariates > pi) = -2*pi + tfVMVariates(tfVMVariates > pi);
tfVMVariates(tfVMVariates < -pi) = 2*pi + tfVMVariates(tfVMVariates < -pi);

% --- END of vmrand.m ---


text:
function [cIX_,gIX_,num] = ArtifactAnalysis(cIX,gIX,CellXYZ,absIX)
U = unique(gIX);
cIX_abs = absIX(cIX);
%%
% Mean = zeros(length(U),3);
% Median = zeros(length(U),3);
% Midpoint = zeros(length(U),3);
Std = zeros(length(U),3);
% TF_clus = zeros(length(U),1);
TF = zeros(length(gIX),1);
for i = 1:length(U),
   IX = gIX==U(i);
   YXZ = CellXYZ(cIX_abs(IX),:);
   for j = 1:3,
%        Mean(i,j) = mean(YXZ(:,j));
%        Median(i,j) = median(YXZ(:,j));
%        Midpoint(i,j) = mean([max(YXZ(:,j)),min(YXZ(:,j))]);
       Std(i,j) = std(YXZ(:,j));
   end
   tf1 = Std(i,1)<0.5;
   tf2 = Std(i,2)<0.5;
   tf3 = Std(i,3)<0.5;
   tf = tf1 || tf2 || tf3;
%    TF_clus(i) = tf;
   TF(IX) = tf;
end
%%
cIX_ = cIX(find(TF));
gIX_ = gIX(find(TF));
num = length(unique(gIX_));
end

author: xiuyechen

repo: FishExplorer


text:
function cplxPot
% complex potential flow
%    using MATLAB analytical solutions
%
%   $Ekkehard Holzbecher  $Date: 2006/05/31 $
%--------------------------------------------------------------------------
% Baseflow
H = 10;             % thickness [L]
h0 = 5;             % reference piezometric head [L]
K = 5.e-5;          % hydraulic conductivity [L/T]
Qx0 = 0;            % baseflow in x-direction [L^2/T]
Qy0 = 0;            % baseflow in y-direction [L^2/T]

% Wells
xwell = [150 250];          % x-coordinates well position [L]
ywell = [0 0];              % y-coordinates well position [L]
Qwell = 1.e-4*[1 -1];       % pumping / recharge rates [L^3/T]
R = [1 1];                  % well radius [L]

% Mesh
xmin = 0;           % minimum x-position of mesh [L]
xmax = 400;         % maximum x-position of mesh [L]
ymin = -100;        % minimum y-position of mesh [L]
ymax = 100;         % maximum y-position of mesh [L]

% Reference point position in mesh
iref = 1; jref = 1;

% Graphical output options
gsurfh = 0;         % piezometric head surface plot
gcontf = 16;        % no. filled contour lines (=0: none)
gquiv = 0;          % arrow field plot
gflowp_fit = 0;     % flowpaths forward in time
gflowp_bit = 0;     % no. flowpaths backward in time (=0: none)
gflowp_dot = 0;     % flowpaths with dots indicating speed
gstream = 10;       % streamfunction plot

%----------------------------------------execution-------------------------------
xvec = linspace(xmin,xmax,50);
yvec = linspace(ymin,ymax,50);
[x,y] = meshgrid (xvec,yvec);                      % mesh

phi = -Qx0*x - Qy0*y;
psi = -Qx0*y + Qy0*x;
for i = 1:size(xwell,2)
    r = sqrt((x-xwell(i)).*(x-xwell(i))+(y-ywell(i)).*(y-ywell(i)));
    phi = phi + (Qwell(i)/(2*pi))*log(r);   % potential
    psi = psi + (Qwell(i)/(2*pi))*atan2((y-ywell(i)),(x-xwell(i)));
end
if h0 > H
    phi0 = -phi(iref,jref) + K*H*h0 - 0.5*K*H*H;
else
    phi0 = -phi(iref,jref) + 0.5*K*h0*h0;       % reference potential
end
hc = 0.5*H+(1/K/H)*(phi+phi0);                  % head confined
hu = sqrt ((2/K)*(phi+phi0));                   % head unconfined
phicrit = phi0 + 0.5*K*H*H;                     % transition confined / unconfined
confined = (phi>=phicrit);                      % confined / unconfined indicator
h = confined.*hc+~confined.*hu;                 % head

%---------------------------------------display messages-------------------
if all(all(confined))
    display ('aquifer confined');
else
    if all(all(~confined))
        display ('aquifer unconfined');
    else
        display ('aquifer partially confined and unconfined');
    end
end
if any(any(h<0))
    display ('aquifer falls partially dry');
    h = max(0, h);
end
[u,v] = gradient (-phi);

%--------------------------------------graphical output--------------------
if gsurfh
    figure; surf (x,y,h);                             % surface
end
figure;
if gcontf                                             % filled contours
    colormap(winter);
    contourf (x,y,h,linspace(5-max(max(h-5)),5+max(max(h-5)),gcontf),'w');
    colorbar; hold on;
end
if gquiv
    quiver (x,y,u,v,'y'); hold on;                    % arrow field
end
if gflowp_fit                                         % flowpaths
    xstart = []; ystart = [];
    for i = 1:100
        if v(1,i) > 0 xstart = [xstart xvec(i)];...
                ystart = [ystart yvec(1)]; end
        if v(100,i) < 0 xstart = [xstart xvec(i)];...
                ystart = [ystart yvec(100)]; end
        if u(i,1) > 0 xstart = [xstart xvec(1)];...
                ystart = [ystart yvec(i)]; end
        if u(i,100) < 0 xstart = [xstart xvec(100)];...
                ystart = [ystart yvec(i)]; end
    end
    h = streamline (x,y,u,v,xstart,ystart);
    set (h,'Color','r');
end
if gflowp_bit
    xstart = x0 + R*cos(2*pi*[1:1:gflowp_bit]/gflowp_bit);
    ystart = y0 + R*sin(2*pi*[1:1:gflowp_bit]/gflowp_bit);
    h = streamline (x,y,-u,-v,xstart,ystart);
    set (h,'Color','y')
end
if gflowp_dot
    [verts averts] = streamslice(x,y,u,v,gflowp_dot);
    sc = 10/mean(mean(sqrt(u.*u+v.*v)));
    iverts = interpstreamspeed(x,y,u,v,verts,sc);
    h = streamline(iverts);
    set (h,'Marker','.','Color','y','MarkerSize',18)
end
if gstream
    h = contour (x,y,psi,gstream,'k','LineWidth',1);
end


text:
function y=agd(x)

% Inverse Gudermannian function
y=atanh(sin(x));

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function hdr = analyzeWrite(imArray, fileName, mmPerVox, notes, origin)
%  hdr = analyzeWrite(imArray, fileName, [mmPerVox], [notes], [origin]);
%
% Saves Analyze7.5 image file from a matlab array 'imArray'.
% fileName should be a root filename- no '.img' or '.hdr' suffix.
% mmPerVox is the row,col,slice voxel size, in mm. Defaults to
% [1 1 1].
%
% If imarray is int16 or int32 then it will be saved as such. If it is
% logical it will be saved as uint1 (1-bit). Any other matlab type will
% be saved as 'float'.
%
% Returns hdr - Analyse header modified after image write.
%
% REQUIRES:
%   fmris_write_analyze from VISTASOFT 'filters' module.
%
% SEE ALSO: loadAnalyze, analyze2mrGray, mrGray2Analyze
%
% HISTORY:
%   2002.04.01 RFD (bob@white.stanford.edu) wrote it.
%   2003.07.09 RFD added origin option.
%   2003.12.04 RFD: minor edits to allow for 4-d data.
%   2004.01.06 RFD: replaced 'V = spm_vol(fileName);' with code to set up
%   the V struct manually. This fixes a bug where an existing image
%   transform 'mat' file would be read in a replace our carefully
%   calculated scale factors. (This bug would not affect those who don't
%   use the spm-style '.mat' transform files.)
%   2004.07.27 RFD: Removed SPM99 dependency. We now use a modified version
%   of Worsley's fmris_write_analyze.
%   2004.10.19 AG: Added version check for Matlab7 compatibility
%   datestr(now) format changed, used date instead
%
% BOB (c) Stanford VISTASOFT, 2002

if (~exist('mmPerVox','var') | isempty(mmPerVox))
    if(ndims(imArray==4)) mmPerVox = [1 1 1 1];
    else mmPerVox = [1 1 1]; end
end

% datestr(now) format changed, used date instead
if (~exist('notes','var') | isempty(notes))
    matlabVersion=ver('Matlab');
    majorVersion = str2num(matlabVersion.Version(1));
    if (majorVersion>=7)
        notes = ['Created by saveAnalyze on ',date];
    else
        notes = ['Created by saveAnalyze on ',datestr(now)];
    end
end

imDim = size(imArray);
if (~exist('origin','var') | isempty(origin))
    origin = (imDim(1:3)+1)./2;
end

d.file_name = [fileName '.img'];
d.data = imArray;
d.vox = mmPerVox;
d.vox_units = 'mm';
d.vox_offset = 0;
d.calib_units = 'min^-1';
d.origin = origin;
d.descrip = notes;
switch(class(imArray))
    case {'uint8','int16','int32'}
        d.precision = class(imArray);
    case {'logical'}
        d.precision = 'uint1';
    otherwise
        d.precision = 'float';
end
disp(['saving data as ' d.precision '...']);
d = fmris_write_analyze(d);
hdr = d;

% s = spm_hwrite(fileName, imDim, mmPerVox, 1, spm_type('int16'), 0, origin, notes);
% %V = spm_vol(fileName);
% V.fname = fileName;
% V.dim = [imDim spm_type('int16')];
% % Offset is from center of image
% offset = -mmPerVox.*origin;
% V.mat = [mmPerVox(1) 0 0 offset(1) ; 0 mmPerVox(2) 0 offset(2) ; 0 0 mmPerVox(3) offset(3) ; 0 0 0 1];
% V.pinfo = [1 0 0]';
% V.descrip = notes;
%
% hdr = spm_write_vol(V, imArray);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft