text:
function schemeCoeffs = computeCoeffs(K, dt, L, M, S)
%COMPUTECOEFFS   Compute coefficients of an EXPINTEG.
%   SCHEMECOEFFS = COMPUTECOEFFS(K, DT, L, M, S) computes the coefficients
%   needed by the EXPINTEG K from the time-step DT, the linear part L, the
%   number of points for complex means M, and the SPINOPERATOR S.
%
% See also EXPINTEG.

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

% Note (1): At the end of the method the coefficients will be stored in
% SCHEMECOEFFS, a STRUCT object with the following fields:
%
% SCHEMECOEFFS.A -> sxs CELL-ARRAY that stores the coefficients A, functions of
%                   the phi-functions
% SCHEMECOEFFS.B -> sx1 CELL-ARRAY that stores the coefficients B, functions of
%                   the phi-functions
% SCHEMECOEFFS.C -> sx1 CELL-ARRAY that stores the numbers C
% SCHEMECOEFFS.E -> (s+1)x1 CELL-ARRAY that stores the s matrix exponentials
%                    e^{C(i)*dt*L}, i=1,..,s, and e^{*dt*L}
% SCHEMECOEFFS.U -> sx(q-1) CELL-ARRAY that stores the coefficients U, functions
%                   of the phi-functions
% SCHEMECOEFFS.V -> (q-1)x1 CELL-ARRAY that stores the coefficients V, functions
%                   of the phi-functions
%
% where s is the number of number of internal stages and q is number of steps
% (>1 for multistep methods). For informations about what these coefficients
% are, see [1].
%
% [1] H. Montanelli and N. Bootland, Solving periodic semilinear stiff PDEs in
% 1D/2D/3D with exponential integrators, submitted (2017).

% Note (2): EXPINTEG schemes are used for solving PDEs with diagonal linear
% operators, e.g., periodic PDEs in 1D/2D/3D with SPIN/SPIN2/SPIN3.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE-PROCESSING:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set-up:
s = K.stages;          % number of stages
q = K.steps;           % number of steps (>1 for multistep methods)
dim = getDimension(S); % spatial dimension (1, 2 or 3)
nVars = S.numVars;     % number of unknown functions
schemeName = K.scheme; % scheme
N = size(L, 1)/nVars;  % grid points

% Initialize the coefficients of the scheme:
A = cell(s);
B = cell(s, 1);
C = zeros(s, 1);
U = cell(s, q-1);
V = cell(q-1, 1);

% Create a contour around each eigenvalue of the linear part L:
LR = computeLR(S, dt, L, M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ETD ADAMS-BASHFORTH:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ( strcmpi(schemeName, 'abnorsett4') == 1 )

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,1} = expinteg.psiEval(1, C(1), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute V:
    V{1} = -3*phi{2} - 5*phi{3} - 3*phi{4};
    V{2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4};
    V{3} = -1/3*phi{2} - phi{3} - phi{4};

elseif ( strcmpi(schemeName, 'abnorsett5') == 1 )

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,1} = expinteg.psiEval(1, C(1), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute V:
    V{1} = -(4*phi{2} + 26/3*phi{3} + 9*phi{4} + 4*phi{5});
    V{2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5};
    V{3} = -(4/3*phi{2} + 14/3*phi{3} + 7*phi{4} + 4*phi{5});
    V{4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};

elseif ( strcmpi(schemeName, 'abnorsett6') == 1 )

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,1} = expinteg.psiEval(1, C(1), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute U:
    V{1} = -(5*phi{2} + 77/6*phi{3} + 71/4*phi{4} + 14*phi{5} + 5*phi{6});
    V{2} = 5*phi{2} + 107/6*phi{3} + 59/2*phi{4} + 26*phi{5} + 10*phi{6};
    V{3} = -(10/3*phi{2} + 13*phi{3} + 49/2*phi{4} + 24*phi{5} + 10*phi{6});
    V{4} = 5/4*phi{2} + 61/12*phi{3} + 41/4*phi{4} + 11*phi{5} + 5*phi{6};
    V{5} = -(1/5*phi{2} + 5/6*phi{3} + 7/4*phi{4} + 2*phi{5} + phi{6});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ETD RUNGE-KUTTA:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( strcmpi(schemeName, 'etdrk2') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
    end

    % Compute A:
    A{2,1} = phi{1};

    % Compute B:
    B{1} = phi{1} - phi{2};
    B{2} = phi{2};

elseif ( strcmpi(schemeName, 'etdrk4') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = psi{1,2};
    A{4,3} = 2*psi{1,2};

    % Compute B:
    B{2} = 2*phi{2} - 4*phi{3};
    B{3} = 2*phi{2} - 4*phi{3};
    B{4} = -phi{2} + 4*phi{3};

elseif ( strcmpi(schemeName, 'exprk5s8') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1/4;
    C(5) = 1/2;
    C(6) = 1/5;
    C(7) = 2/3;
    C(8) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = psi{1,2};
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{1,5} = psi{1,2};
    psi{1,6} = expinteg.psiEval(1, C(6), LR, N, dim, nVars);
    psi{1,7} = expinteg.psiEval(1, C(7), LR, N, dim, nVars);
    psi{1,8} = phi{1};
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{2,4} = expinteg.psiEval(2, C(4), LR, N, dim, nVars);
    psi{2,6} = expinteg.psiEval(2, C(6), LR, N, dim, nVars);
    psi{2,7} = expinteg.psiEval(2, C(7), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);
    psi{3,6} = expinteg.psiEval(3, C(6), LR, N, dim, nVars);
    psi{3,7} = expinteg.psiEval(3, C(7), LR, N, dim, nVars);
    psi{4,6} = expinteg.psiEval(4, C(6), LR, N, dim, nVars);
    psi{4,7} = expinteg.psiEval(4, C(7), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 2*psi{2,2};
    A{4,3} = 2*psi{2,4};
    A{5,3} = -2*psi{2,2} + 16*psi{3,2};
    A{5,4} = 8*psi{2,2} - 32*psi{3,2};
    A{6,4} = 8*psi{2,6} - 32*psi{3,6};
    A{7,4} = -(125/162)*A{6,4};
    A{6,5} = -2*psi{2,6} + 16*psi{3,6};
    A{7,5} = (125/1944)*A{6,4} - (4/3)*psi{2,7} + (40/3)*psi{3,7};
    Phi = (5/32)*A{6,4} - (25/28)*psi{2,6} + (81/175)*psi{2,7} - ...
        (162/25)*psi{3,7} + (150/7)*psi{4,6} + (972/35)*psi{4,7} + 6*phi{4};
    A{8,5} = -(16/3)*phi{2} + (208/3)*phi{3} - 40*Phi;
    A{7,6} = (3125/3888)*A{6,4} + (25/3)*psi{2,7} - (100/3)*psi{3,7};
    A{8,6} = (250/21)*phi{2} - (250/3)*phi{3} + (250/7)*Phi;
    A{8,7} = (27/14)*phi{2} - 27*phi{3} + (135/7)*Phi;

    % Compute B:
    B{6} = (125/14)*phi{2} - (625/14)*phi{3} + (1125/14)*phi{4};
    B{7} = -(27/14)*phi{2} + (162/7)*phi{3} - (405/7)*phi{4};
    B{8} = (1/2)*phi{2} - (13/2)*phi{3} + (45/2)*phi{4};

elseif ( strcmpi(schemeName, 'friedli') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 2*psi{2,2};
    A{4,2} = -26/25*phi{1} + 2/25*phi{2};
    A{4,3} = 26/25*phi{1} + 48/25*phi{2};

    % Compute B:
    B{3} = 4*phi{2} - 8*phi{3};
    B{4} = -phi{2} + 4*phi{3};

elseif ( strcmpi(schemeName, 'hochbruck-ostermann') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;
    C(5) = 1/2;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{1,5} = expinteg.psiEval(1, C(5), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 4*psi{2,2};
    A{4,2} = phi{2};
    A{5,2} = 1/4*phi{2} - phi{3} + 2*psi{2,2} - 4*psi{3,2};
    A{4,3} = A{4,2};
    A{5,3} = A{5,2};
    A{5,4} = psi{2,2} - A{5,2};

    % Compute B:
    B{4} = -phi{2} + 4*phi{3};
    B{5} = 4*phi{2} - 8*phi{3};

elseif ( strcmpi(schemeName, 'krogstad') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 4*psi{2,2};
    A{4,3} = 2*phi{2};

    % Compute B:
    B{2} = 2*phi{2} - 4*phi{3};
    B{3} = 2*phi{2} - 4*phi{3};
    B{4} = -phi{2} + 4*phi{3};

elseif ( strcmpi(schemeName, 'minchev') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 4/25*psi{1,2} + 24/25*psi{2,2};
    A{4,2} = 21/5*phi{2} - 108/5*phi{3};
    A{4,3} = 1/20*phi{1} - 33/10*phi{2} + 123/5*phi{3};

    % Compute B:
    B{2} = -1/10*phi{1} + 1/5*phi{2} - 4*phi{3} + 12*phi{4};
    B{3} = 1/30*phi{1} + 23/5*phi{2} - 8*phi{3} - 4*phi{4};
    B{4} = 1/30*phi{1} - 7/5*phi{2} + 6*phi{3} - 4*phi{4};

elseif ( strcmpi(schemeName, 'strehmel-weiner') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 2*psi{2,2};
    A{4,2} = -2*phi{2};
    A{4,3} = 4*phi{2};

    % Compute B:
    B{3} = 4*phi{2} - 8*phi{3};
    B{4} = -phi{2} + 4*phi{3};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LAWSON:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( strcmpi(schemeName, 'ablawson4') == 1 )

    % Compute the phi-functions:
    phi0 = expinteg.phiEval(0, LR, N, dim, nVars);
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,1} = expinteg.psiEval(1, C(1), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        phi0 = real(phi0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    e2z = phi0.*phi0;
    e3z = e2z.*phi0;
    e4z = e2z.*e2z;

    % Compute B:
    B{1} = 55/24*phi0;

    % Compute V:
    V{1} = -59/24*e2z;
    V{2} = 37/24*e3z;
    V{3} = -3/8*e4z;

elseif ( strcmpi(schemeName, 'lawson4') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi0 = expinteg.phiEval(0, LR, N, dim, nVars);
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi0 = real(phi0);
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{2,1} = 1/2*psi02;
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{1} = 1/6*phi0;
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERALIZED LAWSON:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( strcmpi(schemeName, 'genlawson41') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

elseif ( strcmpi(schemeName, 'genlawson42') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

    % Compute U:
    U{2,1} = -psi{2,2};
    U{3,1} = -psi{2,2} + 1/4;
    U{4,1} = -phi{2} + 1/2*psi02;

    % Compute V:
    V{1} = -phi{2} + 1/3*psi02 + 1/6;

elseif ( strcmpi(schemeName, 'genlawson43') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

    % Compute U:
    U{2,1} = -2*psi{2,2} - 2*psi{3,2};
    U{3,1} = -2*psi{2,2} - 2*psi{3,2} + 5/8;
    U{4,1} = -2*phi{2} - 2*phi{3} + 5/4*psi02;
    U{2,2} = 1/2*psi{2,2} + psi{3,2};
    U{3,2} = 1/2*psi{2,2} + psi{3,2} - 3/16;
    U{4,2} = 1/2*phi{2} + phi{3} - 3/8*psi02;

    % Compute V:
    V{1} = -2*phi{2} - 2*phi{3} + 5/6*psi02 + 1/2;
    V{2} = 1/2*phi{2} + phi{3} - 1/4*psi02 - 1/6;

elseif ( strcmpi(schemeName, 'genlawson44') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);
    psi{4,2} = expinteg.psiEval(4, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

    % Compute U:
    U{2,1} = -3*psi{2,2} - 5*psi{3,2} - 3*psi{4,2};
    U{2,2} = 3/2*psi{2,2} + 4*psi{3,2} + 3*psi{4,2};
    U{2,3} = -1/3*psi{2,2} - 1*psi{3,2} - 1*psi{4,2};
    U{3,1} = U{2,1} + 35/32;
    U{3,2} = U{2,2} - 21/32;
    U{3,3} = U{2,3} + 5/32;
    U{4,1} = -3*phi{2} - 5*phi{3} - 3*phi{4} + 35/16*psi02;
    U{4,2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4} - 21/16*psi02;
    U{4,3} = -1/3*phi{2} - phi{3} - phi{4} + 5/16*psi02;

    % Compute V:
    V{1} = -3*phi{2} - 5*phi{3} - 3*phi{4} + 35/24*psi02 + 1;
    V{2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4} - 7/8*psi02 - 2/3;
    V{3} = -1/3*phi{2} - 1*phi{3} - 1*phi{4} + 5/24*psi02 + 1/6;

elseif ( strcmpi(schemeName, 'genlawson45') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);
    psi{4,2} = expinteg.psiEval(4, C(2), LR, N, dim, nVars);
    psi{5,2} = expinteg.psiEval(5, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/6;

    % Compute U:
    U{2,1} = -4*psi{2,2} - 26/3*psi{3,2} - 9*psi{4,2} - 4*psi{5,2};
    U{2,2} = 3*psi{2,2} + 19/2*psi{3,2} + 12*psi{4,2} + 6*psi{5,2};
    U{2,3} = -4/3*psi{2,2} - 14/3*psi{3,2} - 7*psi{4,2} - 4*psi{5,2};
    U{2,4} = 1/4*psi{2,2} + 88/96*psi{3,2} + 3/2*psi{4,2} + psi{5,2};
    U{3,1} = U{2,1} + 105/64;
    U{3,2} = U{2,2} - 189/128;
    U{3,3} = U{2,3} + 45/64;
    U{3,4} = U{2,4} - 35/256;
    U{4,1} = -4*phi{2} - 26/3*phi{3} - 9*phi{4} - 4*phi{5} + 105/32*psi02;
    U{4,2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5} - 189/64*psi02;
    U{4,3} = -4/3*phi{2} - 14/3*phi{3} - 7*phi{4} - 4*phi{5} + 45/32*psi02;
    U{4,4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5} - 35/128*psi02;

    % Compute V:
    V{1} = -4*phi{2} - 26/3*phi{3} - 9*phi{4} - 4*phi{5} + 35/16*psi02 + ...
        5/3;
    V{2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5} - 63/32*psi02 - ...
        5/3;
    V{3} = -4/3*phi{2} - 14/3*phi{3} - 7*phi{4} - 4*phi{5} + ...
        15/16*psi02 + 5/6;
    V{4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5} - ...
        35/192*psi02 - 1/6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODIFIED GENERALIZED LAWSON:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( strcmpi(schemeName, 'modgenlawson41') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = phi{2} - 1/3*psi02;

elseif ( strcmpi(schemeName, 'modgenlawson42') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} =  1/2*phi{2} + phi{3} - 1/4*psi02;

    % Compute U:
    U{2,1} = -psi{2,2};
    U{3,1} = -psi{2,2} + 1/4;
    U{4,1} = -phi{2} + 1/2*psi02;

    % Compute V:
    V{1} = -1/2*phi{2} + phi{3} + 1/12*psi02;

elseif ( strcmpi(schemeName, 'modgenlawson43') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/3*phi{2} + phi{3} + phi{4} - 5/24*psi02;

    % Compute U:
    U{2,1} = -2*psi{2,2} - 2*psi{3,2};
    U{3,1} = -2*psi{2,2} - 2*psi{3,2} + 5/8;
    U{4,1} = -2*phi{2} - 2*phi{3} + 5/4*psi02;
    U{2,2} = 1/2*psi{2,2} + psi{3,2};
    U{3,2} = 1/2*psi{2,2} + psi{3,2} - 3/16;
    U{4,2} = 1/2*phi{2} + phi{3} - 3/8*psi02;

    % Compute V:
    V{1} = -phi{2} + phi{3} + 3*phi{4} + 5/24*psi02;
    V{2} = 1/6*phi{2} - phi{4} - 1/24*psi02;

elseif ( strcmpi(schemeName, 'modgenlawson44') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);
    psi{4,2} = expinteg.psiEval(4, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5} - 35/192*psi02;

    % Compute U:
    U{2,1} = -3*psi{2,2} - 5*psi{3,2} - 3*psi{4,2};
    U{2,2} = 3/2*psi{2,2} + 4*psi{3,2} + 3*psi{4,2};
    U{2,3} = -1/3*psi{2,2} - 1*psi{3,2} - 1*psi{4,2};

    U{3,1} = U{2,1} + 35/32;
    U{3,2} = U{2,2} - 21/32;
    U{3,3} = U{2,3} + 5/32;

    U{4,1} = -3*phi{2} - 5*phi{3} - 3*phi{4} + 35/16*psi02;
    U{4,2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4} - 21/16*psi02;
    U{4,3} = -1/3*phi{2} - phi{3} - phi{4} + 5/16*psi02;

    % Compute V:
    V{1} = -3/2*phi{2} + 1/2*phi{3} + 6*phi{4} + 6*phi{5} + 35/96*psi02;
    V{2} = 1/2*phi{2} + 1/3*phi{3} - 3*phi{4} - 4*phi{5} - 7/48*psi02;
    V{3} = -1/12*phi{2} - 1/12*phi{3} + 1/2*phi{4} + phi{5} + 5/192*psi02;

elseif ( strcmpi(schemeName, 'modgenlawson45') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1/2;
    C(3) = 1/2;
    C(4) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi02 = expinteg.psiEval(0, C(2), LR, N, dim, nVars);
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);
    psi{1,4} = expinteg.psiEval(1, C(4), LR, N, dim, nVars);
    psi{2,2} = expinteg.psiEval(2, C(2), LR, N, dim, nVars);
    psi{3,2} = expinteg.psiEval(3, C(2), LR, N, dim, nVars);
    psi{4,2} = expinteg.psiEval(4, C(2), LR, N, dim, nVars);
    psi{5,2} = expinteg.psiEval(5, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi02 = real(psi02);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/2;
    A{4,3} = psi02;

    % Compute B:
    B{2} = 1/3*psi02;
    B{3} = B{2};
    B{4} = 12/59*phi{2} + 50/59*phi{3} + 105/59*phi{4} + 120/59*phi{5} - ...
        60/59*phi{6} - 157/944*psi02;

    % Compute U:
    U{2,1} = -4*psi{2,2} - 26/3*psi{3,2} - 9*psi{4,2} - 4*psi{5,2};
    U{2,2} = 3*psi{2,2} + 19/2*psi{3,2} + 12*psi{4,2} + 6*psi{5,2};
    U{2,3} = -4/3*psi{2,2} - 14/3*psi{3,2} - 7*psi{4,2} - 4*psi{5,2};
    U{2,4} = 1/4*psi{2,2} + 88/96*psi{3,2} + 3/2*psi{4,2} + psi{5,2};
    U{3,1} = U{2,1} + 105/64;
    U{3,2} = U{2,2} - 189/128;
    U{3,3} = U{2,3} + 45/64;
    U{3,4} = U{2,4} - 35/256;
    U{4,1} = -4*phi{2} - 26/3*phi{3} - 9*phi{4} - 4*phi{5} + 105/32*psi02;
    U{4,2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5} - 189/64*psi02;
    U{4,3} = -4/3*phi{2} - 14/3*phi{3} - 7*phi{4} - 4*phi{5} + 45/32*psi02;
    U{4,4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5} - 35/128*psi02;

    % Compute V:
    V{1} = -116/59*phi{2} - 34/177*phi{3} + 519/59*phi{4} + ...
        964/59*phi{5} - 600/59*phi{6} + 495/944*psi02;
    V{2} = 57/59*phi{2} + 121/118*phi{3} - 342/59*phi{4} - ...
        846/59*phi{5} + 600/59*phi{6} - 577/1888*psi02;
    V{3} = -56/177*phi{2} - 76/177*phi{3} + 112/59*phi{4} + ...
        364/59*phi{5} - 300/59*phi{6} + 25/236*psi02;
    V{4} = 11/236*phi{2} + 49/708*phi{3} - 33/118*phi{4} - ...
        61/59*phi{5} + 60/59*phi{6} - 181/11328*psi02;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PREDICTOR-CORRECTOR:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif ( strcmpi(schemeName, 'pec423') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute B:
    B{2} = 1/3*phi{2} + phi{3} + phi{4};

    % Compute U:
    U{2,1} = -2*phi{2} - 2*phi{3};
    U{2,2} = 1/2*phi{2} + phi{3};

    % Compute V:
    V{1} = -phi{2} + phi{3} + 3*phi{4};
    V{2} = 1/6*phi{2} - phi{4};

elseif ( strcmpi(schemeName, 'pecec433') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;
    C(3) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/3*phi{2} + phi{3} + phi{4};

    % Compute B:
    B{3} = 1/3*phi{2} + phi{3} + phi{4};

    % Compute U:
    U{2,1} = -2*phi{2} - 2*phi{3};
    U{3,1} = -phi{2} + phi{3} + 3*phi{4};
    U{2,2} = 1/2*phi{2} + phi{3};
    U{3,2} = 1/6*phi{2} - phi{4};

    % Compute V:
    V{1} = -phi{2} + phi{3} + 3*phi{4};
    V{2} = 1/6*phi{2} - phi{4};

elseif ( strcmpi(schemeName, 'pec524') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute B:
    B{2} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};

    % Compute U:
    U{2,1} = -3*phi{2} - 5*phi{3} - 3*phi{4};
    U{2,2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4};
    U{2,3} = -1/3*phi{2} - phi{3} - phi{4};

    % Compute V:
    V{1} = -3/2*phi{2} + 1/2*phi{3} + 6*phi{4} + 6*phi{5};
    V{2} = 1/2*phi{2} + 1/3*phi{3} - 3*phi{4} - 4*phi{5};
    V{3} = -1/12*phi{2} - 1/12*phi{3} + 1/2*phi{4} + phi{5};

elseif ( strcmpi(schemeName, 'pecec534') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;
    C(3) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};

    % Compute B:
    B{3} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};

    % Compute U:
    U{2,1} = -3*phi{2} - 5*phi{3} - 3*phi{4};
    U{2,2} = 3/2*phi{2} + 4*phi{3} + 3*phi{4};
    U{2,3} = -1/3*phi{2} - phi{3} - phi{4};

    U{3,1} = -3/2*phi{2} + 1/2*phi{3} + 6*phi{4} + 6*phi{5};
    U{3,2} = 1/2*phi{2} + 1/3*phi{3} - 3*phi{4} - 4*phi{5};
    U{3,3} = -1/12*phi{2} - 1/12*phi{3} + 1/2*phi{4} + phi{5};

    % Compute V:
    V{1} = -3/2*phi{2} + 1/2*phi{3} + 6*phi{4} + 6*phi{5};
    V{2} = 1/2*phi{2} + 1/3*phi{3} - 3*phi{4} - 4*phi{5};
    V{3} = -1/12*phi{2} - 1/12*phi{3} + 1/2*phi{4} + phi{5};

elseif ( strcmpi(schemeName, 'pec625') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute B:
    B{2} = 1/5*phi{2} + 5/6*phi{3} + 7/4*phi{4} + 2*phi{5} + phi{6};

    % Compute U:
    U{2,1} = -4*phi{2} - 26/3*phi{3} - 9*phi{4} - 4*phi{5};
    U{2,2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5};
    U{2,3} = -4/3*phi{2} - 14/3*phi{3} - 7*phi{4} - 4*phi{5};
    U{2,4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};

    % Compute V:
    V{1} = -2*phi{2} - 1/3*phi{3} + 17/2*phi{4} + 16*phi{5} + 10*phi{6};
    V{2} = phi{2} + 7/6*phi{3} - 11/2*phi{4} - 14*phi{5} - 10*phi{6};
    V{3} = -1/3*phi{2} - 1/2*phi{3} + 7/4*phi{4} + 6*phi{5} + 5*phi{6};
    V{4} = 1/20*phi{2} + 1/12*phi{3} - 1/4*phi{4} - phi{5} - phi{6};

elseif ( strcmpi(schemeName, 'pecec635') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;
    C(3) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/5*phi{2} + 5/6*phi{3} + 7/4*phi{4} + 2*phi{5} + phi{6};

    % Compute B:
    B{3} = 1/5*phi{2} + 5/6*phi{3} + 7/4*phi{4} + 2*phi{5} + phi{6};

    % Compute U:
    U{2,1} = -4*phi{2} - 26/3*phi{3} - 9*phi{4} - 4*phi{5};
    U{2,2} = 3*phi{2} + 19/2*phi{3} + 12*phi{4} + 6*phi{5};
    U{2,3} = -4/3*phi{2} - 14/3*phi{3} - 7*phi{4} - 4*phi{5};
    U{2,4} = 1/4*phi{2} + 11/12*phi{3} + 3/2*phi{4} + phi{5};
    U{3,1} = -2*phi{2} - 1/3*phi{3} + 17/2*phi{4} + 16*phi{5} + 10*phi{6};
    U{3,2} = phi{2} + 7/6*phi{3} - 11/2*phi{4} - 14*phi{5} - 10*phi{6};
    U{3,3} = -1/3*phi{2} - 1/2*phi{3} + 7/4*phi{4} + 6*phi{5} + 5*phi{6};
    U{3,4} = 1/20*phi{2} + 1/12*phi{3} - 1/4*phi{4} - phi{5} - phi{6};

    % Compute V:
    V{1} = -2*phi{2} - 1/3*phi{3} + 17/2*phi{4} + 16*phi{5} + 10*phi{6};
    V{2} = phi{2} + 7/6*phi{3} - 11/2*phi{4} - 14*phi{5} - 10*phi{6};
    V{3} = -1/3*phi{2} - 1/2*phi{3} + 7/4*phi{4} + 6*phi{5} + 5*phi{6};
    V{4} = 1/20*phi{2} + 1/12*phi{3} - 1/4*phi{4} - phi{5} - phi{6};

elseif ( strcmpi(schemeName, 'pec726') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);
    phi{7} = expinteg.phiEval(7, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute B:
    B{2} = 1/6*phi{2} + 137/180*phi{3} + 15/8*phi{4} + 17/6*phi{5} + ...
        5/2*phi{6} + phi{7};

    % Compute U:
    U{2,1} = -5*phi{2} - 77/6*phi{3} - 71/4*phi{4} - 14*phi{5} - 5*phi{6};
    U{2,2} = 5*phi{2} + 107/6*phi{3} + 59/2*phi{4} + 26*phi{5} + 10*phi{6};
    U{2,3} = -10/3*phi{2} - 13*phi{3} - 49/2*phi{4} - 24*phi{5} - 10*phi{6};
    U{2,4} = 5/4*phi{2} + 61/12*phi{3} + 41/4*phi{4} + 11*phi{5} + 5*phi{6};
    U{2,5} = -1/5*phi{2} - 5/6*phi{3} - 7/4*phi{4} - 2*phi{5} - phi{6};

    % Compute V:
    V{1} = -5/2*phi{2} - 17/12*phi{3} + 83/8*phi{4} + 57/2*phi{5} + ...
        65/2*phi{6} + 15*phi{7};
    V{2} = 5/3*phi{2} + 47/18*phi{3} - 8*phi{4} - 92/3*phi{5} - ...
        40*phi{6} - 20*phi{7};
    V{3} = -5/6*phi{2} - 19/12*phi{3} + 29/8*phi{4} + 37/2*phi{5} + ...
        55/2*phi{6} + 15*phi{7};
    V{4} = 1/4*phi{2} + 31/60*phi{3} - phi{4} - 6*phi{5} - 10*phi{6} - ...
        6*phi{7};
    V{5} = -1/30*phi{2} - 13/180*phi{3} + 1/8*phi{4} + 5/6*phi{5} + ...
        3/2*phi{6} + phi{7};

elseif ( strcmpi(schemeName, 'pecec736') == 1 )

    % Compute C:
    C(1) = 0;
    C(2) = 1;
    C(3) = 1;

    % Compute the phi-functions:
    phi{1} = expinteg.phiEval(1, LR, N, dim, nVars);
    phi{2} = expinteg.phiEval(2, LR, N, dim, nVars);
    phi{3} = expinteg.phiEval(3, LR, N, dim, nVars);
    phi{4} = expinteg.phiEval(4, LR, N, dim, nVars);
    phi{5} = expinteg.phiEval(5, LR, N, dim, nVars);
    phi{6} = expinteg.phiEval(6, LR, N, dim, nVars);
    phi{7} = expinteg.phiEval(7, LR, N, dim, nVars);

    % Compute the psi-functions:
    psi{1,2} = expinteg.psiEval(1, C(2), LR, N, dim, nVars);
    psi{1,3} = expinteg.psiEval(1, C(3), LR, N, dim, nVars);

    % Take real part for diffusive problems (real eigenvalues):
    if ( isreal(L) == 1 )
        phi = cellfun(@(f) real(f), phi, 'UniformOutput', 0);
        psi = cellfun(@(f) real(f), psi, 'UniformOutput', 0);
    end

    % Compute A:
    A{3,2} = 1/6*phi{2} + 137/180*phi{3} + 15/8*phi{4} + 17/6*phi{5} + ...
        5/2*phi{6} + phi{7};

    % Compute B:
    B{3} = 1/6*phi{2} + 137/180*phi{3} + 15/8*phi{4} + 17/6*phi{5} + ...
        5/2*phi{6} + phi{7};

    % Compute U:
    U{2,1} = -5*phi{2} - 77/6*phi{3} - 71/4*phi{4} - 14*phi{5} - 5*phi{6};
    U{2,2} = 5*phi{2} + 107/6*phi{3} + 59/2*phi{4} + 26*phi{5} + 10*phi{6};
    U{2,3} = -10/3*phi{2} - 13*phi{3} - 49/2*phi{4} - 24*phi{5} - 10*phi{6};
    U{2,4} = 5/4*phi{2} + 61/12*phi{3} + 41/4*phi{4} + 11*phi{5} + 5*phi{6};
    U{2,5} = -1/5*phi{2} - 5/6*phi{3} - 7/4*phi{4} - 2*phi{5} - phi{6};
    U{3,1} = -5/2*phi{2} - 17/12*phi{3} + 83/8*phi{4} + 57/2*phi{5} + ...
        65/2*phi{6} + 15*phi{7};
    U{3,2} = 5/3*phi{2} + 47/18*phi{3} - 8*phi{4} - 92/3*phi{5} - ...
        40*phi{6} - 20*phi{7};
    U{3,3} = -5/6*phi{2} - 19/12*phi{3} + 29/8*phi{4} + 37/2*phi{5} + ...
        55/2*phi{6} + 15*phi{7};
    U{3,4} = 1/4*phi{2} + 31/60*phi{3} - phi{4} - 6*phi{5} - 10*phi{6} - ...
        6*phi{7};
    U{3,5} = -1/30*phi{2} - 13/180*phi{3} + 1/8*phi{4} + 5/6*phi{5} + ...
        3/2*phi{6} + phi{7};

    % Compute V:
    V{1} = -5/2*phi{2} - 17/12*phi{3} + 83/8*phi{4} + 57/2*phi{5} + ...
        65/2*phi{6} + 15*phi{7};
    V{2} = 5/3*phi{2} + 47/18*phi{3} - 8*phi{4} - 92/3*phi{5} - ...
        40*phi{6} - 20*phi{7};
    V{3} = -5/6*phi{2} - 19/12*phi{3} + 29/8*phi{4} + 37/2*phi{5} + ...
        55/2*phi{6} + 15*phi{7};
    V{4} = 1/4*phi{2} + 31/60*phi{3} - phi{4} - 6*phi{5} - 10*phi{6} - ...
        6*phi{7};
    V{5} = -1/30*phi{2} - 13/180*phi{3} + 1/8*phi{4} + 5/6*phi{5} + ...
        3/2*phi{6} + phi{7};

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POST-PROCESSING:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put everything in SCHEMECOEFFS:
schemeCoeffs.A = A;
schemeCoeffs.B = B;
schemeCoeffs.C = C;
schemeCoeffs.U = U;
schemeCoeffs.V = V;

% Compute the missing oefficients using the summation properties of the coeffs:
% (Unless the scheme is ABLAWSON4, LAWSON4 or ETDRK2.)
if ( ~strcmpi(schemeName, 'lawson4') && ~strcmpi(schemeName, 'ablawson4') && ...
     ~strcmpi(schemeName, 'etdrk2') )
    schemeCoeffs = computeMissingCoeffs(K, schemeCoeffs, phi, psi);
end

% Precompute the E quantities:
E = cell(s+1, 1);
for i = 1:s
    E{i} = exp(C(i)*dt*L);
end
E{s+1} = exp(dt*L);
schemeCoeffs.E = E;

% Multiply by time-step:
schemeCoeffs.A = cellfun(@(A) dt*A, schemeCoeffs.A, 'UniformOutput', 0);
schemeCoeffs.B = cellfun(@(B) dt*B, schemeCoeffs.B, 'UniformOutput', 0);
schemeCoeffs.U = cellfun(@(U) dt*U, schemeCoeffs.U, 'UniformOutput', 0);
schemeCoeffs.V = cellfun(@(V) dt*V, schemeCoeffs.V, 'UniformOutput', 0);

end

function schemeCoeffs = computeMissingCoeffs(K, schemeCoeffs, phi, psi)
%COMPUTEMISSINGCOEFFS   Compute the missing oefficients of an EXPINTEG using
%the summation properties of the coefficients.
%   SCHEMECOEFFS = COMPUTEMISSINGCOEFFS(K, SCHEMECOEFFS, PHI, PSI) uses the row
%   summation properties to compute the SCHEMECOEFFS A{i,1} and B{1} of the
%   EXPINTEG K, using and the phi- and psi-functions.

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

% Get the coefficients:
A = schemeCoeffs.A;
B = schemeCoeffs.B;
U = schemeCoeffs.U;
V = schemeCoeffs.V;

% Number of internal stages S and number of steps used Q:
s = K.stages;
q = K.steps;

% Compute the coefficients A{i,1} using the row summation property:
for i = 2:s
    A{i,1} = psi{1,i};
    for j = 2:i-1
        if ( isempty(A{i,j}) == 0 )
            A{i,1} = A{i,1} - A{i,j};
        end
    end
    for j = 1:q-1
        if ( isempty(U{i,j}) == 0 )
            A{i,1} = A{i,1} - U{i,j};
        end
    end
end

% Compute the coefficient B{1} using the row summation property:
B{1} = phi{1};
for i = 2:s
    if ( isempty(B{i}) == 0 )
        B{1} = B{1} - B{i};
    end
end
for i = 1:q-1
    if ( isempty(V{i}) == 0 )
        B{1} = B{1} - V{i};
    end
end

% Put everything in SCHEMECOEFFS:
schemeCoeffs.A = A;
schemeCoeffs.B = B;

end


text:
function test_suite = test_polygonSubcurve
%TESTPOLYGONSUBCURVE  One-line description here, please.
%   output = testPolygonSubcurve(input)
%
%   Example
%   testPolygonSubcurve
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-16,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

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

p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4];

subcurve = polygonSubcurve(square, .5, 1.5);
testCase.assertEqual(subcurve, [15 10;20 10;20 15], 'AbsTol', .01);

subcurve = polygonSubcurve(square, 3.5, .5);
testCase.assertEqual(subcurve, [10 15;10 10;15 10], 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


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

function test_init()
    ds=cosmo_synthetic_dataset();

    assertTrue(numel(setxor(fieldnames(ds),...
               {'samples','a','sa','fa'}))==0);

function test_nifti()
    ni=generate_test_nifti_struct();
    g=cosmo_synthetic_dataset('size','normal');

    tmpfn='__tmp1_.nii';
    cleaner=onCleanup(@()delete(tmpfn));
    save_nii(ni,tmpfn);
    ds=cosmo_fmri_dataset(tmpfn,'targets',g.sa.targets,...
                                'chunks',g.sa.chunks);

    assertTrue(numel(setxor(fieldnames(ds),...
               fieldnames(g)))==0);

    assertElementsAlmostEqual(ds.samples, g.samples,'relative',1e-6);
    assertEqual(ds.sa.targets, g.sa.targets)
    assertEqual(ds.sa.chunks, g.sa.chunks)

function test_io()
    ds=cosmo_synthetic_dataset('size','normal');
    tmpfn='__tmp2_.nii';
    cleaner=onCleanup(@()delete(tmpfn));

    cosmo_map2fmri(ds,tmpfn);

    es=cosmo_fmri_dataset(tmpfn);
    assertElementsAlmostEqual(ds.samples,es.samples,'relative',1e-6);
    assertEqual(ds.a.fdim,es.a.fdim);

function ni=generate_test_nifti_struct()
% Generates a struct that behaves like a NIFTI struct
    ni=struct();

    ds=cosmo_synthetic_dataset('size','normal');
    ni.img=shiftdim(cosmo_unflatten(ds),1);

    hdr=struct();

    dime=struct();
    dime.datatype=16; %single
    dime.dim=[4 3 2 5 6 1 1 1];
    dime.pixdim=[0 2 2 2 0 0 0 0];
    fns={'intent_p1','intent_p2','intent_p3','intent_code',...
        'slice_start','slice_duration','slice_end',...
        'scl_slope','scl_inter','slice_code','cal_max',...
        'cal_min','toffset'};

    dime=nifti_helper_set_all(dime,fns);
    dime.xyzt_units=10;
    hdr.dime=dime;

    hk=struct();
    hk.sizeof_hdr=348;
    hk.data_type='';
    hk.db_name='';
    hk.extents=0;
    hk.session_error=0;
    hk.regular='r';
    hk.dim_info=0;
    hdr.hk=hk;

    hist=struct();
    hist.sform_code=1;
    hist.originator=[2 1 3 3];
    hist=nifti_helper_set_all(hist,{'descrip','aux_file'},'');
    hist=nifti_helper_set_all(hist,{'qform_code','quatern_b',...
                        'quatern_d',...
                        'qoffset_x','qoffset_y','qoffset_z'});
    hist=nifti_helper_set_all(hist,{'intent_name'},'');

    hist.srow_x=[2 0 0 -1];
    hist.srow_y=[0 2 0 -1];
    hist.srow_z=[0 0 2 -1];

    hist.quatern_c=0;
    hdr.hist=hist;

    ni.hdr=hdr;

function s=nifti_helper_set_all(s, fns, v)
    % sets all fields in fns in struct s to v
    % if v is omitted it is set to 0.
    if nargin<3, v=0; end
    n=numel(fns);
    for k=1:n
        fn=fns{k};
        s.(fn)=v;
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function out = iszero(f)
%ISZERO    True for zero TRIGTECH objects.
%   ISZERO(F) returns logical TRUE is F.VALUES has only zero entries and logical
%   FALSE otherwise.

out = ~any(f.values, 1);

if ( any(out) )
    % We need this as any(NaN) = 0, which will pass the test above.
    out = out & ~any(isnan(f.values), 1);
    return
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% IM = mkGaussian(SIZE, COVARIANCE, MEAN, AMPLITUDE)
%
% Compute a matrix with dimensions SIZE (a [Y X] 2-vector, or a
% scalar) containing a Gaussian function, centered at pixel position
% specified by MEAN (default = (size+1)/2), with given COVARIANCE (can
% be a scalar, 2-vector, or 2x2 matrix.  Default = (min(size)/6)^2),
% and AMPLITUDE.  AMPLITUDE='norm' (default) will produce a
% probability-normalized function.  All but the first argument are
% optional.

% Eero Simoncelli, 6/96.

function [res] = mkGaussian(sz, cov, mn, ampl)

sz = sz(:);
if (size(sz,1) == 1)
  sz = [sz,sz];
end

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('cov') ~= 1)
  cov = (min(sz(1),sz(2))/6)^2;
end

if (exist('mn') ~= 1)
  mn = (sz+1)/2;
end

if (exist('ampl') ~= 1)
  ampl = 'norm';
end

%------------------------------------------------------------

[xramp,yramp] = meshgrid([1:sz(2)]-mn(2),[1:sz(1)]-mn(1));

if (sum(size(cov)) == 2)  % scalar
  if (strcmp(ampl,'norm'))
    ampl = 1/(2*pi*cov(1));
  end
  e = (xramp.^2 + yramp.^2)/(-2 * cov);
elseif (sum(size(cov)) == 3) % a 2-vector
  if (strcmp(ampl,'norm'))
    ampl = 1/(2*pi*sqrt(cov(1)*cov(2)));
  end
  e = xramp.^2/(-2 * cov(2)) + yramp.^2/(-2 * cov(1));
else
  if (strcmp(ampl,'norm'))
    ampl = 1/(2*pi*sqrt(det(cov)));
  end
  cov = -inv(cov)/2;
  e = cov(2,2)*xramp.^2 + (cov(1,2)+cov(2,1))*(xramp.*yramp) ...
      + cov(1,1)*yramp.^2;
end

res = ampl .* exp(e);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function closeOpenWindows(windowName)
    openWindow = findobj('Type', 'Figure', '-and', 'Name', windowName);
    if ~(isempty(openWindow))
        close(openWindow);
    end
end

author: VisLab

repo: EEG-Clean-Tools


text:
function [mg] = slug2mg(slug)
% Convert units of mass from slugs to milligrams.
% Chad A. Greene 2012
mg = slug*14.5939*1000000;

author: Sable

repo: mcbench-benchmarks


text:
function pcof = r8poly_lagrange_coef ( npol, ipol, xpol )

%*****************************************************************************80
%
%% R8POLY_LAGRANGE_COEF returns the coefficients of a Lagrange polynomial.
%
%  Discussion:
%
%    Given distinct abscissas XPOL(1:NPOL), the IPOL-th Lagrange
%    polynomial L(IPOL)(X) is defined as the polynomial of degree
%    NPOL - 1 which is 1 at XPOL(IPOL) and 0 at the NPOL - 1 other
%    abscissas.
%
%    A formal representation is:
%
%      L(IPOL)(X) = Product ( 1 <= I <= NPOL, I /= IPOL )
%       ( X - X(I) ) / ( X(IPOL) - X(I) )
%
%    However sometimes it is desirable to be able to write down
%    the standard polynomial coefficients of L(IPOL)(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NPOL, the number of abscissas.
%    NPOL must be at least 1.
%
%    Input, integer IPOL, the index of the polynomial to evaluate.
%    IPOL must be between 1 and NPOL.
%
%    Input, real XPOL(NPOL), the abscissas of the
%    Lagrange polynomials.  The entries in XPOL must be distinct.
%
%    Output, real PCOF(1:NPOL), the standard polynomial
%    coefficients of the IPOL-th Lagrange polynomial:
%      L(IPOL)(X) = SUM ( 0 <= I <= NPOL-1 ) PCOF(I+1) * X**I
%

%
%  Make sure IPOL is legal.
%
  if ( ipol < 1 || npol < ipol )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY_LAGRANGE_COEF - Fatal error!\n' );
    fprintf ( 1, '  1 <= IPOL <= NPOL is required.\n' );
    error ( 'R8POLY_LAGRANGE_COEF - Fatal error!' );
  end
%
%  Check that the abscissas are distinct.
%
  if ( ~r8vec_distinct ( npol, xpol ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY_LAGRANGE_COEF - Fatal error!\n' );
    fprintf ( 1, '  Two or more entries of XPOL are equal:\n' );
    error ( 'R8POLY_LAGRANGE_COEF - Fatal error!' );
  end

  pcof(1) = 1.0;
  pcof(2:npol) = 0.0;

  indx = 0;

  for i = 1 : npol

    if ( i ~= ipol )

      indx = indx + 1;

      for j = indx : -1 : 0

        pcof(j+1) = -xpol(i) * pcof(j+1) / ( xpol(ipol) - xpol(i) );

        if ( 0 < j )
          pcof(j+1) = pcof(j+1) + pcof(j) / ( xpol(ipol) - xpol(i) );
        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y,arg]=spmax2(row,col,val1,val2,N)
% MATLAB version
% [y,arg]=spmax2(row,col,val1,val2,N)
% is equivalent to [y,arg]=spmax(x)
% where x has N cols and is given by val1+val2 in position(row,col). If val1==-INF => entry disregarded
% -INF in y signify that column was empty.
y=-inf*ones(1,N);
arg=zeros(1,N);
val1=val1(:);
val2=val2(:);
interesting=find(val1~=-inf);
val=val1(interesting)+val2(interesting);
col=col(interesting);
row=row(interesting);
for i=1:N,
   sel=find(col==i);
   if ~isempty(sel),
      [y(i),arg(i)]=max(val(sel));
      arg(i)=row(sel(arg(i)));
   end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function cS = staurolite

%
cs_St = crystalSymmetry('2/m',[0.469 1 0.34], [90, 90, 90]*degree);
N = Miller({1,0,0},{0,1,0},{0,0,1},{1,1,0},{2,0,1},{2,0,-1},cs_St);
dist = [2.43, 1.3, 6.5, 2.45, 9.85, 9.85];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function ori = subSet(S3G,varargin)
% indexing of rotation
%
% Syntax
%   ori = subSet(SO3G,ind) %
%

ori = subSet@orientation(orientation(S3G),varargin{:});


text:
% Duffing_Bif - Bifurcation diagram for the Duffing equation.
% Copyright Springer 2013 Stephen Lynch.
% Make sure Duffing is in your directory.
clear
global Gamma;
Max=120;step=0.001;interval=Max*step;a=1;b=0;
% Ramp the amplitude up.
for n=1:Max
    Gamma=step*n;
    [t,x]=ode45('Duffing',0:(2*pi/1.25):(4*pi/1.25),[a,b]);
    a=x(2,1);
    b=x(2,2);
    rup(n)=sqrt((x(2,1))^2+(x(2,2))^2);
end
% Ramp the amplitude down.
for n=1:Max
    Gamma=interval-step*n;
    [t,x]=ode45('Duffing',0:(2*pi/1.25):(4*pi/1.25),[a,b]);
    a=x(2,1);
    b=x(2,2);
    rdown(n)=sqrt((x(2,1))^2+(x(2,2))^2);
end
hold on
rr=step:step:interval;
plot(rr,rup)
plot(interval-rr,rdown)
hold off
fsize=15;
axis([0 .12 0 2])
xlabel('\Gamma','FontSize',fsize)
ylabel('r','FontSize',fsize)
title('Bifurcation Diagram of the Duffing System. Bistable Region.')


text:
%% Setup
p = get_default_truck_trailer_params();
p.noiseLevel = 0.01;
p.trailerWheelbase = 11;
p.forwardTarget = [20 2 0.3];
xx = [0 -5 0 0];
opt = odeset('Events', @obstacleEvents);
dt = 0.2;
p.forwardGain = [0 0 0];
p.feedbackGain = [0 0 0];
%% setup viz
figure(101)
ax(1) = subplot(131);
tp1 = truck_trailer_plot(p, ax(1));
tp1.updateFig([0 -5 0 0])
title('Initial Pose')
ax(2) = subplot(132);
tp2 = truck_trailer_plot(p, ax(2));
tp2.updateFig([p.forwardTarget 0]);
title('Intermediate Target')
ax(3) = subplot(133);
tp2 = truck_trailer_plot(p, ax(3));
title('Final Target')

%% setup viz
figure(101)
ax(1) = subplot(121);
tp1 = truck_trailer_plot(p, ax(1));
tp1.updateFig([0 -5 0 0])
title('Initial Pose')
ax(2) = subplot(122);
tp2 = truck_trailer_plot(p, ax(2));
title('Desired Target')

set(gcf, 'Color', 'w');
export_fig ttproblem.png

%% setup viz
figure(101)
ax(1) = subplot(121);
tp1 = truck_trailer_plot(p, ax(1));
tp1.updateFig([20 2 0.3 0])
title('Intermediate Target #1')
ax(2) = subplot(122);
tp2 = truck_trailer_plot(p, ax(2));
tp2.updateFig([20 0 0  0])
title('Intermediate Target #2')
%%
set(gcf, 'Color', 'w');
export_fig ttinter.png

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% Local Regression and Likelihood, Figure 4.1.
% Author: Catherine Loader
%
% Local Likelihood (Poisson Regression).

load mine;
fit = locfit(extrp,frac,'family','poisson','deg',1,'alpha',0.6);
figure('Name','fig4_1: Poisson Regression');
lfplot(fit);
lfband(fit);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
filename='Chair_Tetrahedra';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';
constraint_case = 'INEQUALITY';
line_search_initiator = 'STANDARD';
showBC = true;

nsteps = 15;
Vfrac_final = 0.15;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ x, w ] = tetrahedron_arbq ( degree, n )

%*****************************************************************************80
%
%% TETRAHEDRON_ARBQ returns a quadrature rule for a tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired degree of exactness.
%    1 <= DEGREE <= 15.
%
%    Input, integer N, the number of points in the rule.
%    This value should be requested first from TETRAHEDRON_ARBQ_SIZE.
%
%    Output, real X(3,N), the quadrature nodes.
%
%    Output, real W(N), the quadrature weights.
%
  if ( degree == 1 )
    [ x, w ] = rule01 ( n );
  elseif ( degree == 2 )
    [ x, w ] = rule02 ( n );
  elseif ( degree == 3 )
    [ x, w ] = rule03 ( n );
  elseif ( degree == 4 )
    [ x, w ] = rule04 ( n );
  elseif ( degree == 5 )
    [ x, w ] = rule05 ( n );
  elseif ( degree == 6 )
    [ x, w ] = rule06 ( n );
  elseif ( degree == 7 )
    [ x, w ] = rule07 ( n );
  elseif ( degree == 8 )
    [ x, w ] = rule08 ( n );
  elseif ( degree == 9 )
    [ x, w ] = rule09 ( n );
  elseif ( degree == 10 )
    [ x, w ] = rule10 ( n );
  elseif ( degree == 11 )
    [ x, w ] = rule11 ( n );
  elseif ( degree == 12 )
    [ x, w ] = rule12 ( n );
  elseif ( degree == 13 )
    [ x, w ] = rule13 ( n );
  elseif ( degree == 14 )
    [ x, w ] = rule14 ( n );
  elseif ( degree == 15 )
    [ x, w ] = rule15 ( n );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TETRAHEDRON_ARBQ - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of DEGREE.\n' );
    error ( 'TETRAHEDRON_ARBQ - Fatal error!' );
  end

  d = sum ( w(1:n) );

  volume = sqrt ( 8.0 ) / 3.0;
  w(1:n) = w(1:n) * volume / d;

  return
end


text:
function [ pols, dersx, dersy ] = klegeypols3 ( x, y, n )

%*****************************************************************************80
%
%% KLEGEYPOLS3 evaluate scaled Legendre polynomials and derivatives.
%
%  Discussion:
%
%    This routine evaluates a sequence of scaled Legendre polynomials
%    P_n(x/y) y^n, with the parameter y in [0,1], together with their
%    derivatives with respect to the parameters x and y.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    27 June 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Input, real Y, the parameter value.
%
%    Input, integer N, the highest degree to be evaluated.
%
%    Output, real POLS(N+1), the polynomial values.
%
%    Output, real DERSX(N+1), the derivatives with respect to X.
%
%    Output, real DERSY(N+1), the derivatives with respect to Y.
%
  pkp1 = 1.0;
  pols(1) = pkp1;
  dkp1 = 0.0;
  dersx(1) = dkp1;
  ykp1 = 0.0;
  dersy(1) = ykp1;

  if ( n == 0 )
    return
  end

  pk = pkp1;
  pkp1 = x;
  pols(2) = pkp1;
  dk = dkp1;
  dkp1 = 1.0;
  dersx(2) = dkp1;
  yk = ykp1;
  ykp1 = 0.0;
  dersy(2) = ykp1;

  if ( n == 1 )
    return
  end

  for k = 1 : n - 1
    pkm1 = pk;
    pk = pkp1;
    dkm1 = dk;
    dk = dkp1;
    ykm1 = yk;
    yk = ykp1;
    pkp1 = ( ( 2.0 * k + 1.0 ) * x * pk - k * pkm1 * y * y ) / ( k + 1.0 );
    dkp1 = ( ( 2.0 * k + 1.0 ) * ( x * dk + pk ) ...
      - k * dkm1 * y * y ) / ( k + 1.0 );
    ykp1 = ( ( 2.0 * k + 1.0 ) * ( x * yk ) ...
      - k * ( pkm1 * 2.0 * y + ykm1 * y * y ) ) / ( k + 1.0 );
    pols(k+2) = pkp1;
    dersx(k+2) = dkp1;
    dersy(k+2) = ykp1;
  end

  return
end


text:
function ensure(condition)

%   ensure.m
%   $Id: ensure.m 800 2008-02-26 22:32:04Z mpf $
%
%   ----------------------------------------------------------------------
%   This file is part of SPGL1 (Spectral Projected Gradient for L1).
%
%   Copyright (C) 2007 Ewout van den Berg and Michael P. Friedlander,
%   Department of Computer Science, University of British Columbia, Canada.
%   All rights reserved. E-mail: <{ewout78,mpf}@cs.ubc.ca>.
%
%   SPGL1 is free software; you can redistribute it and/or modify it
%   under the terms of the GNU Lesser General Public License as
%   published by the Free Software Foundation; either version 2.1 of the
%   License, or (at your option) any later version.
%
%   SPGL1 is distributed in the hope that it will be useful, but WITHOUT
%   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
%   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
%   Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with SPGL1; if not, write to the Free Software
%   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
%   USA
%   ----------------------------------------------------------------------

   if ~condition
      [ST,I] = dbstack;
      msg = 'Assertion failed';
      if length(ST) > 1
         msg = sprintf('%s in %s; function %s line %d', ...
                       msg, ST(2).file, ST(2).name, ST(2).line);
      end

      errstr.message = msg;
      error(errstr);
   end

end % function ensure

author: mpf


text:
function r8vec_variance_test ( )

%*****************************************************************************80
%
%% R8VEC_VARIANCE_TEST tests R8VEC_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_VARIANCE_TEST\n' );
  fprintf ( 1, '  R8VEC_VARIANCE computes the variance of an R8VEC.\n' );;

  n = 10;
  r8_lo = -5.0;
  r8_hi = +5.0;
  seed = 123456789;
  [ a, seed ] = r8vec_uniform_ab ( n, r8_lo, r8_hi, seed );

  r8vec_print ( n, a, '  Input vector:' );

  variance = r8vec_variance ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Variance:   %f\n', variance );

  return
end


text:
function sincos_test ( )

%*****************************************************************************80
%
%% SINCOS_TEST demonstrates SPINTERP on a function of a 2D argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SINCOS_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate the use of SPINTERP to construct an\n' );
  fprintf ( 1, '  interpolant to the function z(x,y) = sin(x) + cos(y)\n' );
%
%  We need to have the spinterp program in the Matlab path.
%
  addpath ( '../spinterp' );
%
%  To alter the default spinterp options, we must call spset().
%  We want to use a grid that uses Chebyshev spacing, and we
%  want to use at least 5 levels of interpolation.
%
  OPTIONS = spset ( 'GridType', 'Chebyshev', ...
                    'MinDepth', 5 );
%
%  Now we call spvals() to set up in C the data defining
%  the sparse grid interpolant.  C is a Matlab structure.
%
  m = 2;

  box = [ 0.0, pi; ...
          0.0, pi ];

  c = spvals ( @sincos_f, m, box, OPTIONS );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the sparse interpolant structure:\n' );
  fprintf ( 1, '\n' );

  c
%
%  Just for information, print out the grid points added at each level.
%
  l_max = 4;
  for l = 0 : l_max
    fprintf ( 1, '\n' );
    label = sprintf ( '  Grid points added at level %d', l );
    x = spgrid ( l, m, OPTIONS );
    [ n, ~ ] = size ( x );
    r8mat_print ( n, m, x, label );
  end

  figure ( 1 )
  l_max = 5;
  plotgrid ( l_max, m );
  grid on
  xlabel ( '<---X--->' );
  ylabel ( '<---Y--->' );
  title ( sprintf ( 'Spinterp Chebyshev Grid of level %d', l_max ) )

  filename = 'sincos_grid.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Grid plot stored as "%s".\n', filename );
%
%  Evaluate the interpolant at some random points in the region.
%
  n = 25;
  x = pi * rand ( 1, n );
  y = pi * rand ( 1, n );
  z = spinterp ( c, x, y );

  e = max ( abs ( z - sincos_f ( x, y ) ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum approximation error at %d random points is %g\n', n, e );
%
%  Display plots of the function and interpolant.
%
  figure ( 2 )
  subplot ( 1, 2, 1 );
  ezmesh ( @sincos_f, [ 0.0, pi ] );
  title ( 'z(x,y) = sin(x) + cos(y)' );

  subplot ( 1, 2, 2 );
  ezmesh ( @(x,y) spinterp ( c, x, y ), [ 0, pi ] );
  title ( 'Sparse grid interpolant' );

  filename = 'sincos_interp.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Interpolant plot stored as "%s".\n', filename );

  rmpath ( '../spinterp' )
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SINCOS_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% new procedure: generate motor seed/map in trialRes space
% find top 100 cells correlating to motor, screen with Rh4+Rh5(+Rh6??) masks

% corr sweep (0.5-0.7) to get top %2 cells, then kmeans, save; rank by stim-lock, save
clear all; close all; clc

%% folder setup
isSaveFig = 1;
isPlotFig = 1;

outputDir = GetOutputDataDir;
saveDir1 = fullfile(outputDir,'motor_map_030317');
saveDir2 = fullfile(outputDir,'motor_map_seed_030317');
if ~exist(saveDir1, 'dir'), mkdir(saveDir1), end;
if ~exist(saveDir2, 'dir'), mkdir(saveDir2), end;

%% init

hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

setappdata(hfig,'isMotorseed',0);

%% run fish
range_fish = GetFishRange;%[1:3,5:18];
M_thres_reg = zeros(3,18);
M_numTopCorr = zeros(1,18);
M_motorseedRegs = cell(1,18);
M_compareMotorCellNumber = zeros(2,18);

for i_fish = 14 %range_fish
    ClusterIDs = [1,1];
    [~,~,~,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);

    numTopCorr = 100;
    [~,~,regressor_m_raw] = GetMotorRegressor(behavior,i_fish);

    %% convert to trialRes (converting locally, not updating hfig to isTrialRes)
    [~,regressor_m_raw_tRes] = GetTrialAvrLongTrace(hfig,regressor_m_raw);
    [~,M_0_tRes] = GetTrialAvrLongTrace(hfig,M_0);

    %%
    M_lr = [1,3];
%     M_name = {'Left','Right'};
    cIX_seed = [];
    gIX_seed = [];
    for i_lr = 1:2,
        % left:
        reg_tRes = regressor_m_raw_tRes(M_lr(i_lr),:);

        %% for each cell, find correlation coeff
        M_corr = corr(reg_tRes',M_0_tRes');%cell_resp_ave');

        [~,I] = sort(M_corr,'descend');
        cIX = I(1:numTopCorr)';
        gIX = i_lr*ones(size(cIX));

        %% screen with anat masks #222 and #233 (Rhombomore 4 and 5)
        MASKs = getappdata(hfig,'MASKs');
        CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');
        absIX = getappdata(hfig,'absIX');

        Msk_IDs = [222,223,224]; % manual input
        [cIX,gIX] = ScreenCellsWithMasks(Msk_IDs,cIX,gIX,MASKs,CellXYZ_norm,absIX);

        % enforce that seed cells for the left are in the left hemisphere,
        % and vice versa for the right
        if i_lr ==1,
            [cIX,gIX]  = DivideCellsbyHemisphere(CellXYZ_norm,absIX,cIX,gIX);
        else
            [~,~,cIX,gIX] = DivideCellsbyHemisphere(CellXYZ_norm,absIX,cIX,gIX);
        end

        num = numTopCorr;
        % need at least 10 cells to qualify as a seed
        while length(cIX)<10, % repeat with larger numTopCorr
            num = num+100;
            cIX = I(1:num)';
            gIX = i_lr*ones(size(cIX));
            [cIX,gIX] = ScreenCellsWithMasks(Msk_IDs,cIX,gIX,MASKs,CellXYZ_norm,absIX);
            M_numTopCorr(i_fish) = num;

            % screen by left/right again
            if i_lr ==1,
                [cIX,gIX]  = DivideCellsbyHemisphere(CellXYZ_norm,absIX,cIX,gIX);
            else
                [~,~,cIX,gIX] = DivideCellsbyHemisphere(CellXYZ_norm,absIX,cIX,gIX);
            end
        end

        cIX_seed = [cIX_seed;cIX];
        gIX_seed = [gIX_seed;gIX];

    end
    % this is not tRes, just loading default:
    M = UpdateIndices_Manual(hfig,cIX_seed,gIX_seed);
    Reg = FindClustermeans(gIX_seed,M);
    % save motor seed functional trace, e.g. new motor regressors
    M_motorseedRegs{i_fish} = Reg;

    %% save motor seeds in VAR, and visualize (not in tRes)
    clusgroupID = 11;
    clusIDoverride = 2;
    name = 'Motor_seed_030317';
    SaveCluster_Direct(cIX_seed,gIX_seed,absIX,i_fish,name,clusgroupID,clusIDoverride);

    %% left-right combined plot
    figure('Position',[50,100,1400,800]);
    % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
    subplot(121)
    setappdata(hfig,'isPlotBehavior',1);
    setappdata(hfig,'isStimAvr',0);
    setappdata(hfig,'isPlotLines',0);
    UpdateTimeIndex(hfig);
    DrawTimeSeries(hfig,cIX_seed,gIX_seed);

    % right plot
    ax = subplot(122)
    I = LoadCurrentFishForAnatPlot(hfig,cIX_seed,gIX_seed);
    DrawCellsOnAnat(I,ax);

    if isSaveFig,
        filename = fullfile(saveDir2, ['Fish',num2str(i_fish),'_motormap_seed']);
        saveas(gcf, filename, 'png');
        close(gcf)
    end

    %% regression (in tRes), thresholding by % of cells (instead of corr thres)

    [~,Reg_tRes] = GetTrialAvrLongTrace(hfig,Reg);

    Corr = corr(Reg_tRes',M_0_tRes');
    [corr_max,IX] = max(Corr,[],1);
    [~,I] = sort(corr_max,'descend');

    M_target_prctcell= [1,2,3];
    for i_numcell = 2,
        nCells_total = size(M_0,1);
        prctcell = M_target_prctcell(i_numcell);
        nCells_target = round(prctcell/100 * nCells_total);

        %% target cell number: counting cells in hindbrain only

        Msk_IDs = 114; % mask for full hindbrain

        % isScreenMskFromAllCells
        cIX = (1:length(absIX))';
        gIX = ones(size(cIX));
        [cIX_hb,gIX_hb] = ScreenCellsWithMasks(Msk_IDs,cIX,gIX,MASKs,CellXYZ_norm,absIX);

        I_hb = ismember(I,cIX_hb);
        cum_I_hb = cumsum(I_hb);
        lastIX = find(cum_I_hb==nCells_target);

        % cut off the desired number of cells to display
        cIX = I(1:lastIX)';
        gIX = IX(cIX)';
        M = UpdateIndices_Manual(hfig,cIX,gIX);

        M_thres_reg(i_numcell,i_fish) = corr_max(I(nCells_target));
        M_compareMotorCellNumber(1,i_fish) = length(cIX);

        %% k-means (not in tRes)
        numK = 16; % manual
        gIX = Kmeans_Direct(M,numK);

        %% Plot anat
        if isPlotFig,
            %             figure('Position',[50,100,800,1000]);
            %             I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
            %             DrawCellsOnAnat(I);
            figure('Position',[50,100,1400,800]);
            % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
            subplot(121)
            setappdata(hfig,'isPlotBehavior',1);
            setappdata(hfig,'isStimAvr',0);
            setappdata(hfig,'isPlotLines',0);
            UpdateTimeIndex(hfig);
            DrawTimeSeries(hfig,cIX,gIX);

            % right plot
            ax = subplot(122)
            I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
            DrawCellsOnAnat(I,ax);

            %% Save plot
            if isSaveFig,
                filename = fullfile(saveDir1, ['Fish',num2str(i_fish),'_motormap_',num2str(prctcell),'%']);
                saveas(gcf, filename, 'png');
                close(gcf)
            end
        end
        %%
%         %% Rank by stim-lock
%         [gIX,rankscore] = RankByStimLock_Direct(hfig,gIX);
%
%         setappdata(hfig,'rankscore',round(rankscore*100)/100);
%         setappdata(hfig,'clrmap_name','jet');
%
%         % Plot anat
%         if isPlotFig,
%             figure('Position',[50,100,800,1000]);
%             I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
%             DrawCellsOnAnat(I);
%
%
%             % Save plot
%             if isSaveFig,
%                 filename = fullfile(saveDir2, ['Fish',num2str(i_fish),'_motormap_',num2str(target_numcell),'_stimlock_anat']);
%                 saveas(gcf, filename, 'png');
%                 close(gcf)
%             end
%
%             % Plot functional traces
%             figure('Position',[50,100,800,1000]);
%             % isCentroid,isPlotLines,isPlotBehavior,isPlotRegWithTS
%             setappdata(hfig,'isPlotBehavior',1);
%             setappdata(hfig,'isStimAvr',0);
%             UpdateTimeIndex(hfig);
%             DrawTimeSeries(hfig,cIX,gIX);
%
%             % Save plot
%             isSaveFig
%             filename = fullfile(saveDir2, ['Fish',num2str(i_fish),'_motormap_',num2str(target_numcell),'_stimlock_func']);
%             saveas(gcf, filename, 'png');
%             close(gcf)
%         end
%     end
%     % reset colormap
%     setappdata(hfig,'clrmap_name','hsv_new');

    end

    %% count cell numbers regressed (in tRes) with fictive trace
    thres_reg = M_thres_reg(2,i_fish);
%     [~,~,regressor_m_raw] = GetMotorRegressor(behavior,i_fish);

    [rawmotorcorr,IX] = max(corr(regressor_m_raw_tRes',M_0_tRes'),[],1);
    ix = find(rawmotorcorr>thres_reg);
    M_compareMotorCellNumber(2,i_fish) = length(ix);
end
SaveVARwithBackup();

%%
% range_fish excludes Fish 4
M_compareMotorCellNumber(:,4) = NaN;
figure;bar(M_compareMotorCellNumber')

author: xiuyechen

repo: FishExplorer

path: github-repos/MATLAB/xiuyechen-FishExplorer/FishExplorer-c61392cf0835480d64fc03c15f1992935fdc7106/figure scripts/SensoryMotor/arc code/motormap, tRes versions, comparisons, w Clusters or Autoclus to follow/fig3_motormap_trialRes_030317.m


text:
function [K] = ku1v1(x, y, xp, yp, hyp, i)

logsigma = hyp(1);
logthetax = hyp(2);
logthetay = hyp(3);

n_x = size(x,1);
n_y = size(y,1);
n_xp = size(xp,1);
n_yp = size(yp,1);

x = repmat(x,1,n_xp);
y = repmat(y,1,n_yp);
xp = repmat(xp',n_x,1);
yp = repmat(yp',n_y,1);

switch i

case 0

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*(y+(-1).*yp);

case 1 % logsigma

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*(y+(-1).*yp);

case 2 % logthetax

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*((-1)+(1/2).*exp(1).^((-1).*logthetax).*(x+(-1).*xp).^2).*(x+ ...
  (-1).*xp).*(y+(-1).*yp);

case 3 % logthetay

K=exp(1).^(logsigma+(-1).*logthetax+(-1).*logthetay+(-1/2).*exp(1).^((-1) ...
  .*logthetax).*(x+(-1).*xp).^2+(-1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(x+(-1).*xp).*((-1)+(1/2).*exp(1).^((-1).*logthetay).*(y+(-1) ...
  .*yp).^2).*(y+(-1).*yp);

otherwise

        K = zeros(n_x, n_xp);
end

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function vis_benchmark(timing)
% display results of benchmarking (timing)

fprintf('\n\n------------------------------\n');
fnames = fieldnames(timing);
ttot = 0; % total sum of timing
for k=1:length(fnames);
    tsec = timing.(fnames{k});
    if isnan(tsec), continue; end;
    fprintf('%s\t: %0.5g\n',fnames{k},tsec);
    ttot = ttot+tsec;
end
fprintf('Total\t: %0.5g\n',ttot);
fprintf('------------------------------\n\n');

repo: REST


text:
function output = isfinite(input)

output = cellfun(@isfinite, input, 'uniformoutput', false);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Evaluate cost function for closest representative search given coefficients
%function ft=essential_distMinAnglePair_ft(t,m1,p1,c1,m2,p2,c2)
%Evaluates the cost function used by essential_distMinAnglePair to find the
%closest representative in the equivalence class of a QREM
%If m2,p2,c2 are omitted or empty, get value of a single term
function ft=essential_distMinAnglePair_ft(t,m1,p1,c1,m2,p2,c2)
flagSingleTerm=false;
if ~exist('m2','var') || isempty(m2)
    flagSingleTerm=true;
end

if flagSingleTerm
    ft=acos((m1*sin(t+p1)+c1-1)/2)^2;
else
    ft=acos((m1*sin(t+p1)+c1-1)/2)^2+acos((m2*sin(t+p2)+c2-1)/2)^2;
end

author: NicolasBoumal

repo: manopt


text:
function b = iswholenum(x)
% Reason for this function is that isinteger() looks at the class of the number
% rather it's actual value. So isinteger(4) will return true but integer(4.0) will
% be false even though 4.0 = 4. This is a problem.
b = false;
if ~ismember(class(x), {'double','single','int8','uint8','int16','uint16','int32','uint32'});
    return;
end
b = (x == floor(x));