text:
function model = Learning_MLE_Basis_Stoc( Seqs, model, alg )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Learning Hawkes processes via maximum likelihood estimation
% Different regularizers (low-rank, sparse, group sparse) of parameters and
% their combinations are considered, which are solved via StochasticADMM.
%
% Reference:
% Xu, Hongteng, Mehrdad Farajtabar, and Hongyuan Zha.
% "Learning Granger Causality for Hawkes Processes."
% International Conference on Machine Learning (ICML). 2016.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June. 10, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial
Aest = model.A;
muest = model.mu;
% KerFint = struct('dGK', []);
% KerF = struct('gij', []);

if alg.LowRank
    UL = zeros(size(Aest));
    ZL = Aest;
end

if alg.Sparse
    US = zeros(size(Aest));
    ZS = Aest;
end

if alg.GroupSparse
    UG = zeros(size(Aest));
    ZG = Aest;
end

D = size(Aest, 1);

if alg.storeLL
    model.LL = zeros(alg.epoch, 1);
end
if alg.storeErr
    model.err = zeros(alg.epoch, 3);
end

tic;
for o = 1:alg.epoch % the number of epoch for stochastic opt

    rho = alg.rho * (1.1^o);

    idx_seq = randperm(length(Seqs));
    if length(Seqs)<alg.seqbatch;
        N_seq_batch = 1;
        Ns = length(Seqs);
    else
        N_seq_batch = floor(length(Seqs)/alg.seqbatch);
        Ns = alg.seqbatch;
    end

    for n = 1:N_seq_batch

        % initial para
        Amu = zeros(D, 1);
        Bmu = Amu;

        CmatA = zeros(size(Aest));
        AmatA = CmatA;
        BmatA = CmatA;
        if alg.LowRank
            BmatA = BmatA + rho*(UL-ZL);
            AmatA = AmatA + rho;
        end
        if alg.Sparse
            BmatA = BmatA + rho*(US-ZS);
            AmatA = AmatA + rho;
        end
        if alg.GroupSparse
            BmatA = BmatA + rho*(UG-ZG);
            AmatA = AmatA + rho;
        end

        % E-step: evaluate the responsibility using the current parameters
        for ns = 1:Ns
            c = idx_seq(Ns*(n-1)+ns); % the id of current sequence

            if ~isempty(Seqs(c).Time)
                Time = Seqs(c).Time;
                Event = Seqs(c).Mark;
                Tstart = Seqs(c).Start;

                if isempty(alg.Tmax)
                    Tstop = Seqs(c).Stop;
                else
                    Tstop = alg.Tmax;
                    indt = Time < alg.Tmax;
                    Time = Time(indt);
                    Event = Event(indt);
                end

                Nc = length(Time);
%                 if o == 1
%                     KerFint(c).dGK = zeros(alg.historyL, length(model.landmark), Nc-1);
%                     KerF(c).gij = zeros(alg.historyL, length(model.landmark), Nc-1);
%                 end

                idx_event = randperm(Nc);
                if Nc<alg.eventbatch
                    N_event_batch = 1;
                    Ne = Nc;
                else
                    N_event_batch = alg.eventbatch;
                    Ne = floor(Nc/alg.eventbatch);
                end
                TimeNew = [Tstart, Time];
                for m = 1:N_event_batch
                    for ne = 1:Ne
                        i = idx_event(Ne*(m-1)+ne);
                        ui = Event(i);
                        ti = Time(i);
                        ti1 = TimeNew(i);

                        Amu = Amu + ti - ti1;

                        % calculate intensity
                        lambdai = muest(ui);
                        pii = muest(ui);
                        pij = [];

                        if i>1
                            start = max([1, i-alg.historyL]);
                            tj = Time(start:i-1);
                            uj = Event(start:i-1);

                            %if o==1 || isempty(KerF(c).gij)
                                dt = ti - tj;
                                dt1 = ti1 - tj;
                                gij = Kernel(dt, model);

                                %KerF(c).gij(1:length(dt),:,i-1) = gij;

                                GK = Kernel_Integration(dt, model);
                                GK1 = Kernel_Integration(dt1, model);
                                dGK = GK - GK1;
                                %KerFint(c).dGK(1:length(dt),:,i-1) = dGK;
                            %end

                            auiuj = Aest(uj, :, ui);
                            %pij = auiuj .* KerF(c).gij(1:length(uj),:,i-1);%
                            pij = auiuj .* gij;
                            lambdai = lambdai + sum(pij(:));
                        end
                        pii = pii./lambdai;

                        if i>1
                            pij = pij./lambdai;
                            if ~isempty(pij) && sum(pij(:))>0
                                for j = 1:length(uj)
                                    uuj = uj(j);
                                    CmatA(uuj,:,ui) = CmatA(uuj,:,ui) - pij(j,:);

                                    BmatA(uuj,:,:) = BmatA(uuj,:,:) + ...
                                        double(Aest(uuj,:,:)>0).*...
                                        repmat( dGK(j,:), [1,1,D] );
                                        %repmat( KerFint(c).dGK(j,:,i-1), [1,1,D] );
                                end
                            end
                        end
                        Bmu(ui) = Bmu(ui) + pii;

                    end
                end

%                 Amu = Amu + Tstop - Tstart;
%
%                 dT = Tstop - Time;
%                 GK = Kernel_Integration(dT, model);
%
%                 Nc = length(Time);
%                 idx_event = randperm(Nc);
%
%                 if Nc<alg.eventbatch
%                     N_event_batch = 1;
%                     Ne = Nc;
%                 else
%                     N_event_batch = alg.eventbatch;
%                     Ne = floor(Nc/alg.eventbatch);
%                 end
%
%                 for m = 1:N_event_batch
%                     for ne = 1:Ne
%                         i = idx_event(Ne*(m-1) + ne);
%                         ui = Event(i);
%
%                         BmatA(ui,:,:) = BmatA(ui,:,:)+...
%                             double(Aest(ui,:,:)>0).*repmat( GK(i,:), [1,1,D] );
%
%                         ti = Time(i);
%                         lambdai = muest(ui);
%                         pii = muest(ui);
%                         pij = [];
%
%                         if i>1
%                             start = max([1, i-alg.historyL]);
%                             tj = Time(start:i-1);
%                             uj = Event(start:i-1);
%
%                             dt = ti - tj;
%                             gij = Kernel(dt, model);
%                             auiuj = Aest(uj, :, ui);
%                             pij = auiuj .* gij;
%                             lambdai = lambdai + sum(pij(:));
%                         end
%
%                         pii = pii./lambdai;
%
%                         if i>1
%                             pij = pij./lambdai;
%                             if ~isempty(pij) && sum(pij(:))>0
%                                 for j = 1:length(uj)
%                                     uuj = uj(j);
%                                     CmatA(uuj,:,ui) = CmatA(uuj,:,ui) - pij(j,:);
%                                 end
%                             end
%                         end
%
%                         Bmu(ui) = Bmu(ui) + pii;
%                     end
%                 end
            else
                warning('Sequence %d is empty!', c)
            end
        end

        % M-step: update parameters
        mu = Bmu./Amu;
        if alg.Sparse==0 && alg.GroupSparse==0 && alg.LowRank==0
            A = -CmatA./BmatA;%( -BA+sqrt(BA.^2-4*AA*CA) )./(2*AA);
            A(isnan(A))=0;
            A(isinf(A))=0;
        else
            A = ( -BmatA + sqrt(BmatA.^2 - 4*AmatA.*CmatA) )./(2*AmatA);
            A(isnan(A))=0;
            A(isinf(A))=0;
        end

        % check convergence
        Err=sum(sum(sum(abs(A-Aest))))/sum(abs(Aest(:)));
        Aest = A;
        muest = mu;
        model.A = Aest;
        model.mu = muest;
        fprintf('Epoch=%d, EventBatch=%d, RelErr=%f, Time=%0.2fsec\n',...
                o, n, Err, toc);

        if Err<alg.thres || (o==alg.epoch && n==N_seq_batch)
            break;
        end
    end

    % calculate Loglikelihood
    if alg.storeLL
        Loglike = Loglike_Basis( Seqs, model, alg );
        model.LL(o) = Loglike;
    end
    % calculate error
    if alg.storeErr
        Err = zeros(1,3);
        Err(1) = norm(model.mu(:) - alg.truth.mu(:))/norm(alg.truth.mu(:));
        Err(2) = norm(model.A(:) - alg.truth.A(:))/norm(alg.truth.A(:));
        Err(3) = norm([model.mu(:); model.A(:)]-[alg.truth.mu(:); alg.truth.A(:)])...
            /norm([alg.truth.mu(:); alg.truth.A(:)]);
        model.err(o,:) = Err;
    end

    if alg.LowRank
        threshold = alg.alphaLR/rho;
        ZL = SoftThreshold_LR( Aest+UL, threshold );
        UL = UL + (Aest-ZL);
    end

    if alg.Sparse
        threshold = alg.alphaS/rho;
        ZS = SoftThreshold_S( Aest+US, threshold );
        US = US + (Aest-ZS);
    end

    if alg.GroupSparse
        threshold = alg.alphaGS/rho;
        ZG = SoftThreshold_GS( Aest+UG, threshold );
        UG = UG + (Aest-ZG);
    end

end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function y = Fault_decision_system(u)
global st14

% Using model defined by the structure st14 to make decision about system fault

u(3)=1e-10*u(3).*(u(end).^6);u=u(1:3);
sker=st14.x2sup+(abs(u))'*abs(u)*ones(st14.Nlsup,1)-2*st14.xsup*abs(u);
y=(st14.w)'*exp(-sker./(2*(st14.sigma).^2))+st14.b;

if y>=0.8
    y=1;
else
    y=0;
end
end

author: Sable

repo: mcbench-benchmarks


text:
function [reg_spectrum, chi2_regNNLS] = do_regNNLS(decay_matrix, measurements, sigma, mu)

%-------------------------------------------------------------------
% function [s_NNLS, chi2_NNLS] = do_regNNLS(num_t2_vals, decay_matrix, measurements, sigma, mu)
%
% * Runs through regularized NNLS (energy)
%
% ~~~ Charmaine Chia (May 3, 2005) ~~~
%-------------------------------------------------------------------

% generate minimum energy regularization
num_t2_vals = size(decay_matrix,2);
H = eye(num_t2_vals);
H_temp = mu*H;
this_decay_matrix = [decay_matrix; H_temp];
this_data = [measurements'; zeros(num_t2_vals,1)];

%-------------------------------------------------------------------
% --- EAO: addition to fit 2D spectrum (T2* times, with freq. offsets)
%-------------------------------------------------------------------
fit_delf = 1;

if fit_delf == 1

    % compute default tolerance used by lsqnonneg (as per help lsqnonneg )
    default_tolx = 10*max(size(this_decay_matrix))*norm(this_decay_matrix,1)*eps;
    % initialize fitting options structure
    opts = optimset('TolX',default_tolx);

    [reg_spectrum, bla, blo, exitflag, output] = lsqnonneg(this_decay_matrix,  this_data);

    % if num of iterations exceeded
    while ~exitflag
        % increase the tolerance
        opts = optimset(opts,'TolX',opts.TolX*10);
        % re-fit
        [reg_spectrum, bla, blo, exitflag] = lsqnonneg(this_decay_matrix, this_data, opts);
    end

else
    reg_spectrum = lsqnonneg(this_decay_matrix, this_data);
end
%-------------------------------------------------------------------

chi2_regNNLS = sum((this_decay_matrix*reg_spectrum - this_data).^2)/sigma^2;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [ nxy, bxy, fxy ] = pce_legendre_linear_assemble ( n, p )

%*****************************************************************************80
%
%% PCE_LEGENDRE_LINEAR_ASSEMBLE assembles a particular stochastic Galerkin matrix.
%
%  Discussion:
%
%    We wish to analyze a stochastic PDE of the form:
%
%      -div A(X,Y) grad U(X,Y) = F(X)
%
%    where
%
%      U is an unknown scalar function,
%      A is a given diffusion function,
%      F is a given right hand side function,
%      X represents dependence on spatial variables,
%      Y represents dependence on stochastic variables.
%
%    We let X be a space of finite element functions generated by piecewise linear
%    functions associated with a particular triangular dissection of the unit square.
%
%    We let Y be the space of polynomials over R^N with total degree at most P.
%
%    We seek solutions U in X cross Y using a polynomial chaos expansion approach.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of factors in the probability space.
%
%    Input, integer P, the maximum degree of the basis functions for the
%    probability space.
%
%    Output, integer NXY, the order of the matrix BXY.  NXY = NX * NY.
%
%    Output, real sparse BXY(NXY,NXY), the matrix.
%
%    Output, real FXY(NXY), the right hand side.
%
%  Local Parameters:
%
%    Local, real AREA_X, the area of the current element.
%
%    Local, integer BASIS_X_NUM, the number of finite element basis functions (3).
%
%    Local, integer BASIS_Y_DEGREE(1:N), describes a particular basis function
%    in Y, in terms of the degrees of its components.
%
%    Local, integer BASIS_Y_H(*), auxilliary "external" memory needed by
%    SUBCOMP_NEXT.
%
%    Local, logical BASIS_Y_MORE, is set FALSE before the first call to
%    SUBCOMP_NEXT.  SUBCOMP_NEXT returns the next subcomposition, and sets the
%    value of BASIS_Y_MORE to TRUE if there are even more subpositions that can
%    be produced, or FALSE if there are no more.
%
%    Local, integer BASIS_Y_T(*), auxilliary "external" memory needed by
%    SUBCOMP_NEXT.
%
%    Local, real BVEC(1:N,1:POINT_X_NUM), stores the KL expansion coefficients
%    of orders 1 through N at the spatial points POINT_X_P(1:2,1:POINT_X_NUM).
%
%    Local, real BZERO(1:POINT_X_NUM), stores the KL expansion coefficients
%    of order 0 at the spatial points POINT_X_P(1:2,1:POINT_X_NUM).
%
%    Local, real CLK(1:QUAD_X_NUM), contains the integral of the KL expansion
%    function for given values of the L-th probability basis function and the
%    K-th probability test function evaluated at the quadrature points in an
%    element.
%
%    Local, integer ELEMENT_X, the element being considered.
%
%    Local, integer ELEMENT_X_NUM, the number of elements in the finite element
%    problem.
%
%    Local, integer ELEMENT_X_NODE(1:3,1:ELEMENT_X_NUM);
%    ELEMENT_X_NODE(I,J) is the global index of local node I in element J.
%
%    Local, integer ELEMENT_X1_NUM, the number of (pairs) of elements in the
%    first spatial dimension in the finite element problem.
%
%    Local, integer ELEMENT_X2_NUM, the number of (pairs) of elements in the
%    second spatial dimension in the finite element problem.
%
%    Local, real FVEC(1:QUAD_X_NUM), stores the values of the function F which
%    is the right hand side of the PDE, at the quadrature points in an element.
%
%    Local, integer I, the index of the finite element test function
%    PHI(1:NX)(X).
%
%    Local, integer J, the index of the finite element basis function
%    PHI(1:NX)(X).
%
%    Local, integer K, the index of the probability space test function
%    PSI(1:NY)(Y).
%
%    Local, integer L, the index of the probability space basis function
%    PSI(1:NY)(Y).
%
%    Local, integer NODE_X_NUM, the number of nodes.
%
%    Local, real NODE_X_P(1:2,1:NODE_X_NUM), the coordinates of nodes.
%
%    Local, integer NODE_X1_NUM, the number of nodes in the first space
%    direction.
%
%    Local, integer NODE_X2_NUM, the number of nodes in the second space
%    direction.
%
%    Local, integer NX, the dimension of the finite element space.
%    NX = NODE_X_NUM for a finite element problem involving a scalar variable U.
%
%    Local, integer NY, = ( N + P)! / N! / P! = the dimension of the space Y,
%    the space of polynomials over R^N of total degree at most P.
%
%    Local, real PHI(1:BASIS_X_NUM,1:QUAD_X_NUM), the finite element basis
%    functions, evaluated at all the quadrature points in a particular element.
%
%    Local, real PHI_DX1(1:BASIS_X_NUM,1:QUAD_X_NUM), the derivative of the
%    finite element basis functions with respect to the first spatial variable,
%    evaluated at all the quadrature points in a particular element.
%
%    Local, real PHI_DX2(1:BASIS_X_NUM,1:QUAD_X_NUM), the derivative of the
%    finite element basis functions with respect to the second spatial variable,
%    evaluated at all the quadrature points in a particular element.
%
%    Local, real PHYS_X_P(1:2,1:QUAD_X_NUM), the "physical" coordinates of the
%    quadrature points in the current element.
%
%    Local, integer QUAD_X, the index of the current X quadrature point.
%
%    Local, integer QUAD_X_NUM, the order of the X quadrature rule.
%
%    Local, real QUAD_X_P(1:2,1:QUAD_X_NUM), the points for the X quadrature
%    rule.
%
%    Local, real QUAD_X_W(1:QUAD_X_NUM), the weights for the X quadrature rule.
%
%    Local, real T3(1:2,1:3), the coordinates of the nodes that define the
%    current element.
%
%    Local, real TABLE(1:P+1,1:P+1), a table of the values of integrals of the
%    1D basis functions in Z.  TABLE(D1+1,D2+1)
%    = Integral ( -1 <= Z <= +1 ) Z * PSI(D1)(Z) * PSI(D2)(Z) dZ.
%
%    Local, integer TEST_X_NUM, the number of finite element test functions (3).
%
%    Local, integer TEST_Y_DEGREE(1:N), describes a particular test function in
%    Y, in terms of the degrees of its components.
%
%    Local, integer TEST_Y_H(*), auxilliary "external" memory needed by
%    SUBCOMP_NEXT.
%
%    Local, logical TEST_Y_MORE, is set FALSE before the first call to
%    SUBCOMP_NEXT.  SUBCOMP_NEXT returns the next subcomposition, and sets the
%    value of TESST_Y_MORE to TRUE if there are even more subpositions that can
%    be produced, or FALSE if there are no more.
%
%    Local, integer TEST_Y_T(*), auxilliary "external" memory needed by
%    SUBCOMP_NEXT.
%
  FALSE = 0;
%
%  Some setup for the finite element calculation.
%
  basis_x_num = 3;
  test_x_num = 3;
  node_x1_num = 11;
  node_x2_num = 11;
  node_x_num = node_x1_num * node_x2_num;
  nx = node_x_num;

  node_x_p = grid_nodes_01 ( node_x1_num, node_x2_num );

  element_x1_num = node_x1_num - 1;
  element_x2_num = node_x2_num - 1;
  element_x_num = 2 * element_x1_num * element_x2_num;
  element_x_node = grid_t3_element ( element_x1_num, element_x2_num );

  quad_x_num = 13;
  [ quad_x_p, quad_x_w ] = triangle_rule_13 ( );
%
%  Some setup for the probability space.
%
  ny = i4_choose ( n + p, n );

  e = 1;
  table = legendre_linear_product ( p, e );

  for i = 1 : p + 1
    for j = 1 : p + 1
      if ( abs ( table(i,j) ) < 10000 * eps )
        table(i,j) = 0.0;
      end
    end
  end
%
%  Define the order of the matrix BXY,
%  Define BXY as a sparse matrix;
%  Define FXY as a column vector.
%
  nxy = nx * ny;
  bxy = sparse ( [], [], [], nxy, nxy );
  fxy(1:nxy,1) = 0.0;
%
%  LOOP L:
%  Generate the L-th basis function PSI(D,Y) in Y space.
%
  basis_y_degree(1:n) = 0;
  basis_y_more = FALSE;
  basis_y_h = [];
  basis_y_t = [];
  basis_y_n2 = [];
  basis_y_more2 = [];
  l = 0;

  while ( 1 )

    l = l + 1;

    [ basis_y_degree, basis_y_more, basis_y_h, basis_y_t, basis_y_n2, basis_y_more2 ] = ...
      subcomp_next ( p, n, basis_y_degree, basis_y_more, basis_y_h, basis_y_t, basis_y_n2, ...
      basis_y_more2 );
%
%  LOOP K:
%  Generate the K-th test function PSI(D,Y) in Y space.
%
    test_y_degree(1:n) = 0;
    test_y_more = FALSE;
    test_y_h = [];
    test_y_t = [];
    test_y_n2 = [];
    test_y_more2 = [];
    k = 0;

    while ( 1 )

      k = k + 1;

      [ test_y_degree, test_y_more, test_y_h, test_y_t, test_y_n2, test_y_more2 ] = ...
        subcomp_next ( p, n, test_y_degree, test_y_more, test_y_h, test_y_t, test_y_n2, ...
        test_y_more2 );

      for element_x = 1 : element_x_num
%
%  For ALL the quadrature points in this element, evaluate:
%  * PHI the basis functions;
%  * PHI_DX1 and PHI_DX2, basis function derivatives;
%  * PHYS_X_P, the physical coordinates of the quadrature points;
%  * BZERO and BVEC, the coefficient functions in the KL expansion for CLK;
%  * CLK, the KL function integrated against the probability basis and test functions;
%  * FVEC, the right hand side of the PDE.
%
        t3(1:2,1:3) = node_x_p(1:2,element_x_node(1:3,element_x));

        area_x = abs ( triangle_area_2d ( t3 ) );

        [ phi, phi_dx1, phi_dx2 ] = basis_mn_t3 ( t3, quad_x_num, quad_x_p );

        phys_x_p(1:2,1:quad_x_num) = reference_to_physical_t3 ( t3, quad_x_num, quad_x_p );

        [ bzero, bvec ] = b_evaluator ( n, quad_x_num, phys_x_p );

        if ( l == k )
          clk(1:quad_x_num) = bzero(1:quad_x_num);
        else
          clk(1:quad_x_num) = 0.0;
        end

        for i2 = 1 : n

          factor = 1.0;
          for i3 = 1 : n
            if ( i3 ~= i2 )
              if ( basis_y_degree(i3) ~= test_y_degree(i3) )
                factor = 0.0;
              end
            end
          end

          clk(1:quad_x_num) = clk(1:quad_x_num) ...
            + bvec(i2,1:quad_x_num) * table(basis_y_degree(i2)+1,test_y_degree(i2)+1)...
            * factor;

        end

        fvec = f_evaluator ( quad_x_num, phys_x_p );
%
%  Integrate over the finite element space X.
%
        for quad_x = 1 : quad_x_num
%
%  LOOP J:
%  All finite element basis functions J.
%  But we actually only look at those which are nonzero in this element.
%
          for basis_x = 1 : basis_x_num

            j = element_x_node(basis_x,element_x);
%
%  LOOP I:
%  Finite element test function I.
%  But we actually only look at those which are nonzero in this element.
%
            for test_x = 1 : test_x_num

              i = element_x_node(test_x,element_x);

              ik = ( k - 1 ) * node_x_num + i;
              jl = ( l - 1 ) * node_x_num + j;
%
%  If I is a boundary node, the equation associated with its test function
%  is replaced by a simple equation that enforces a zero Dirichlet condition.
%
              if ( i <= node_x1_num || ...
                        node_x1_num * ( node_x2_num - 1) + 1 <= i || ...
                        mod ( i, node_x1_num ) == 1 || ...
                        mod ( i, node_x1_num ) == 0 )

                if ( ik == jl )
                  bxy(ik,jl) = 1.0;
                  fxy(ik) = 0.0;
                end
%
%  If I is not a boundary node, the equation associated with its test function
%  generates a finite element equation.
%
              else

                bxy(ik,jl) = bxy(ik,jl) + quad_x_w(quad_x) * area_x ...
                  * clk(quad_x) ...
                  * ( phi_dx1(test_x,quad_x) * phi_dx1(basis_x,quad_x) ...
                    + phi_dx2(test_x,quad_x) * phi_dx2(basis_x,quad_x) );
%
%  The only nonzero result for the right hand side FXY occurs when K = 1,
%  since this is the single test function in Y which is the product of N
%  copies of the constant Legendre polynomial.
%
%  (This is true as long as there is no Y dependence in F(), and
%  assuming we are using Legendre polynomials.)
%
                if ( k == 1 )
                  psi0_integral = 1 / sqrt ( 2^n );
                  fxy(ik) = fxy(ik) + quad_x_w(quad_x) * area_x ...
                    * fvec(quad_x) * phi(test_x,quad_x) * psi0_integral;
                end

              end

            end
          end
        end

      end

      if ( test_y_more == FALSE )
        break
      end

    end

    if ( basis_y_more == FALSE )
      break
    end

  end

  return
end
function [ bzero, bvec ] = b_evaluator ( n, nx, x )

%*****************************************************************************80
%
%% B_EVALUATOR evaluates KL expansion coefficients for the diffusion coefficient.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of factors in the probability space.
%
%    Input, integer NX, the number of points.
%
%    Input, real X(1:2,1:NX), the point coordinates.
%
%    Output, real BZERO(1:NX), the zero-th order coefficient evaluated at X.
%
%    Output, real BVEC(1:N,1:NX), the coefficients of orders 1 to N,
%    evaluated at X.
%
  bzero(1:nx) = 1.0;
  bvec = zeros ( n, nx );
  for i = 1 : n
    bvec(i,1:nx) = sin ( i * x(1,1:nx) );
  end

  return
end
function [ phi, dphidx, dphidy ] = basis_mn_t3 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_T3: all bases functions at N points for a T3 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices of a triangle.
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%    The sides of the triangle DO NOT have to lie along a coordinate
%    axis.
%
%    The routine evaluates the basis functions associated with each vertex,
%    and their derivatives with respect to X and Y.
%
%  Physical Element T3:
%
%            3
%           / \
%          /   \
%         /     \
%        /       \
%       1---------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the vertices of the triangle.  It is common to list
%    these points in counter clockwise order.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(2,N), the coordinates of the evaluation points.
%
%    Output, real PHI(3,N), the basis functions at the evaluation points.
%
%    Output, real DPHIDX(3,N), DPHIDY(3,N), the basis derivatives
%    at the evaluation points.
%
%  Local parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%
  area = t(1,1) * ( t(2,2) - t(2,3) ) ...
       + t(1,2) * ( t(2,3) - t(2,1) ) ...
       + t(1,3) * ( t(2,1) - t(2,2) );

  phi(1,1:n) =     (  ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) )     ...
                    - ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) );
  dphidx(1,1:n) =   - ( t(2,3) - t(2,2) );
  dphidy(1,1:n) =     ( t(1,3) - t(1,2) );

  phi(2,1:n) =     (  ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) )     ...
                    - ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) );
  dphidx(2,1:n) =   - ( t(2,1) - t(2,3) );
  dphidy(2,1:n) =     ( t(1,1) - t(1,3) );

  phi(3,1:n) =     (  ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) )     ...
                    - ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) );
  dphidx(3,1:n) =   - ( t(2,2) - t(2,1) );
  dphidy(3,1:n) =     ( t(1,2) - t(1,1) );
%
%  Normalize.
%
  phi(1:3,1:n)    = phi(1:3,1:n) / area;
  dphidx(1:3,1:n) = dphidx(1:3,1:n) / area;
  dphidy(1:3,1:n) = dphidy(1:3,1:n) / area;

  return
end
function [ a, more, h, t ] = comp_next ( n, k, a, more, h, t )

%*****************************************************************************80
%
%% COMP_NEXT computes the compositions of the integer N into K parts.
%
%  Discussion:
%
%    A composition of the integer N into K parts is an ordered sequence
%    of K nonnegative integers which sum to N.  The compositions (1,2,1)
%    and (1,1,2) are considered to be distinct.
%
%    The routine computes one composition on each call until there are no more.
%    For instance, one composition of 6 into 3 parts is
%    3+2+1, another would be 6+0+0.
%
%    On the first call to this routine, set MORE = FALSE.  The routine
%    will compute the first element in the sequence of compositions, and
%    return it, as well as setting MORE = TRUE.  If more compositions
%    are desired, call again, and again.  Each time, the routine will
%    return with a new composition.
%
%    However, when the LAST composition in the sequence is computed
%    and returned, the routine will reset MORE to FALSE, signaling that
%    the end of the sequence has been reached.
%
%    This routine originally used a SAVE statement to maintain the
%    variables H and T.  I have decided (based on wasting an
%    entire morning trying to track down a problem) that it is safer
%    to pass these variables as arguments, even though the user should
%    never alter them.  This allows this routine to safely shuffle
%    between several ongoing calculations.
%
%  Example:
%
%    The 28 compositions of 6 into three parts are:
%
%      6 0 0,
%      5 1 0,
%      5 0 1,
%      4 2 0,
%      4 1 1,
%      4 0 2,
%      3 3 0,
%      3 2 1,
%      3 1 2,
%      3 0 3,
%      2 4 0,
%      2 3 1,
%      2 2 2,
%      2 1 3,
%      2 0 4,
%      1 5 0,
%      1 4 1,
%      1 3 2,
%      1 2 3,
%      1 1 4,
%      1 0 5,
%      0 6 0,
%      0 5 1,
%      0 4 2,
%      0 3 3,
%      0 2 4,
%      0 1 5,
%      0 0 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2007
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis and Herbert Wilf
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Second Edition,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the integer whose compositions are desired.
%
%    Input, integer K, the number of parts in the composition.
%
%    Input, integer A(K), the previous composition.  On the first call,
%    with MORE = FALSE, set A = [].  Thereafter, A should be the
%    value of A output from the previous call.
%
%    Input, logical MORE.  The input value of MORE on the first
%    call should be FALSE, which tells the program to initialize.
%    On subsequent calls, MORE should be TRUE, or simply the
%    output value of MORE from the previous call.
%
%    Input, integer H, T, two internal parameters needed for the
%    computation.  The user may need to initialize these before the
%    very first call, but these initial values are not important.
%    The user should not alter these parameters once the computation
%    begins.
%
%    Output, integer A(K), the next composition.
%
%    Output, logical MORE, will be TRUE unless the composition
%    that is being returned is the final one in the sequence.
%
%    Output, integer H, T, the updated values of the two internal
%    parameters.
%
  if ( ~more )

    t = n;
    h = 0;
    a(1) = n;
    a(2:k) = 0;

  else

    if ( 1 < t )
      h = 0;
    end

    h = h + 1;
    t = a(h);
    a(h) = 0;
    a(1) = t - 1;
    a(h+1) = a(h+1) + 1;

  end

  more = ( a(k) ~= n );

  return
end
function fvec = f_evaluator ( n, x )

%*****************************************************************************80
%
%% F_EVALUATOR evaluates the finite element right hand side function F.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of point.
%
%    Input, real X(2,N), the point coordinates.
%
%    Output, real FVEC(N), the value of F at the points.
%
  fvec(1:n) = sin ( pi * x(1,1:n) ) .* sin ( pi * x(2,1:n) );

  return
end
function node_xy = grid_nodes_01 ( x_num, y_num )

%*****************************************************************************80
%
%% GRID_NODES_01 returns an equally spaced grid of nodes in the unit square.
%
%  Example:
%
%    X_NUM = 5
%    Y_NUM = 3
%
%    NODE_XY =
%    ( 0, 0.25, 0.5, 0.75, 1, 0,   0.25, 0.5, 0.75, 1,   0, 0.25, 0.5, 0.75, 1;
%      0, 0,    0,   0,    0, 0.5, 0.5,  0.5, 0.5,  0.5, 1, 1.0,  1.0, 1.0,  1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, Y_NUM, the number of nodes in the X and Y directions.
%
%    Output, real NODE_XY(2,X_NUM*Y_NUM), the coordinates of the nodes.
%
  node_num = x_num * y_num;

  node_xy(1:2,1:node_num) = 0.0;

  for i = 1 : x_num
    node_xy(1,i:x_num:i+(y_num-1)*x_num) = ( i - 1 ) / ( x_num - 1 );
  end

  for j = 1 : y_num
    node_xy(2,1+(j-1)*x_num:j*x_num) = ( j - 1 ) / ( y_num - 1 );
  end

  return
end
function [ element_node ] = grid_t3_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T3_ELEMENT produces a grid of pairs of 3 node triangles.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  2,  5;
%         6,  5,  2;
%         2,  3,  6;
%         7,  6,  3;
%         3,  4,  7;
%         8,  7,  4;
%         5,  6,  9;
%        10,  9,  6;
%         6,  7, 10;
%        11, 10,  7;
%         7,  8, 11;
%        12, 11,  8.
%
%  Grid:
%
%    9---10---11---12
%    |\ 8 |\10 |\12 |
%    | \  | \  | \  |
%    |  \ |  \ |  \ |
%    |  7\|  9\| 11\|
%    5----6----7----8
%    |\ 2 |\ 4 |\ 6 |
%    | \  | \  | \  |
%    |  \ |  \ |  \ |
%    |  1\|  3\|  5\|
%    1----2----3----4
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    2 * NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(3,2*NELEMX*NELEMY), the nodes that form
%    each element.
%

%
%  Node labeling:
%
%    NW--NE
%     |\ |
%     | \|
%    SW--SE
%
  element = 0;
  element_node = zeros ( 3, 2 * nelemx * nelemy );
  for j = 1 : nelemy
    for i = 1 : nelemx

      sw = i     + ( j - 1 ) * ( nelemx + 1 );
      se = i + 1 + ( j - 1 ) * ( nelemx + 1 );
      nw = i     +   j       * ( nelemx + 1 );
      ne = i + 1 +   j       * ( nelemx + 1 );

      element = element + 1;

      element_node(1,element) = sw;
      element_node(2,element) = se;
      element_node(3,element) = nw;

      element = element + 1;

      element_node(1,element) = ne;
      element_node(2,element) = nw;
      element_node(3,element) = se;

    end
  end

  return
end
function value = i4_choose ( n, k )

%*****************************************************************************80
%
%% I4_CHOOSE computes the binomial coefficient C(N,K).
%
%  Discussion:
%
%    The value is calculated in such a way as to avoid overflow and
%    roundoff.  The calculation is done in integer arithmetic.
%
%    The formula used is:
%
%      C(N,K) = N! / ( K! * (N-K)! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    ML Wolfson, HV Wright,
%    Algorithm 160:
%    Combinatorial of M Things Taken N at a Time,
%    Communications of the ACM,
%    Volume 6, Number 4, April 1963, page 161.
%
%  Parameters:
%
%    Input, integer N, K, are the values of N and K.
%
%    Output, integer VALUE, the number of combinations of N
%    things taken K at a time.
%
  mn = min ( k, n - k );

  if ( mn < 0 )

    value = 0;

  elseif ( mn == 0 )

    value = 1;

  else

    mx = max ( k, n - k );
    value = mx + 1;

    for i = 2 : mn
      value = ( value * ( mx + i ) ) / i;
    end

  end

  return
end
function table = legendre_linear_product ( p, e )

%*****************************************************************************80
%
%% LEGENDRE_LINEAR_PRODUCT computes a linearly weighted Legendre product table.
%
%  Discussion:
%
%    Let L(i)(X) represent the Legendre polynomial of degree i.
%
%    For polynomial chaos applications, it is of interest to know the
%    value of the integrals of products of X with every possible pair
%    of basis functions.  That is, we'd like to form
%
%      Tij = Integral ( -1 <= X <= +1 ) X^E * L(i)(X) * L(j)(X) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polyonomial factors.
%    0 <= P.
%
%    Input, integer E, the exponent of X in the integrand.
%    0 <= E.
%
%    Output, real TABLE(P+1,P+1), the table of integrals.  TABLE(I,J)
%    represents the weighted integral of X^E * L(i+1)(X) * L(j+1)(X).
%
  table(1:p+1,1:p+1) = 0.0;

  order = p + 1 + floor ( ( e + 1 ) / 2 );
  [ x_table, w_table ] = legendre_quadrature_rule ( order );

  for k = 1 : order

    x = x_table(k);
    l_table = legendre_polynomial ( p, x );
%
%  The following formula is an outer product in L_TABLE.
%
    if ( e == 0 )
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) *       l_table(1:p+1)' * l_table(1:p+1);
    else
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) * x^e * l_table(1:p+1)' * l_table(1:p+1);
    end

  end

  return
end
function l = legendre_polynomial ( p, x )

%*****************************************************************************80
%
%% LEGENDRE_POLYNOMIAL evaluates the Legendre polynomials L(0:P)(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Daniel Zwillinger, editor,
%    CRC Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, integer P, the highest evaluation degree.
%
%    Input, real X, the evaluation point.
%
%    Output, real L(1:P+1), the Legendre polynomials of order 0 through P at X.
%
  if ( p < 0 )
    l = [];
    return
  end
%
%  Allocate space.
%
  l(1:p+1) = 0.0;
%
%  Apply recursion.
%
  l(1) = 1.0;

  if ( 1 <= p )

    l(2) = x;

    for i = 2 : p

      l(i+1) = ( ( 2 * i - 1 ) * x * l(i)   ...
               - (     i - 1 ) *     l(i-1) ) ...
               / (     i     );

    end

  end
%
%  Normalize.
%
  l(1:p+1) = l(1:p+1) .* sqrt ( ( 1 : 2 : 2*p+1 ) / 2 );

  return
end
function [ xtab, weight ] = legendre_quadrature_rule ( order )

%*****************************************************************************80
%
%% LEGENDRE_QUADRATURE_RULE computes a Gauss-Legendre quadrature rule.
%
%  Discussion:
%
%    The integration interval is [ -1, 1 ].
%
%    The weight function is w(x) = 1.0.
%
%    The integral to approximate:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2008
%
%  Author:
%
%    FORTRAN77 original version by Philip Davis, Philip Rabinowitz
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%    ORDER must be greater than 0.
%
%    Output, real XTAB(ORDER), the abscissas of the rule.
%
%    Output, real WEIGHT(ORDER), the weights of the rule.
%    The weights are positive, symmetric, and should sum to 2.
%
  xtab = zeros ( order, 1 );
  weight = zeros ( order, 1 );

  if ( order < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_RULE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of ORDER = %d\n', order );
    error ( 'LEGENDRE_RULE - Fatal error!' );
  end

  e1 = order * ( order + 1 );

  m = floor ( ( order + 1 ) / 2 );

  for i = 1 : floor ( ( order + 1 ) / 2 )

    mp1mi = m + 1 - i;

    t = ( 4 * i - 1 ) * pi / ( 4 * order + 2 );

    x0 = cos(t) * ( 1.0 - ( 1.0 - 1.0 / ( order ) ) / ( 8 * order * order ) );

    pkm1 = 1.0;
    pk = x0;

    for k = 2 : order
      pkp1 = 2.0 * x0 * pk - pkm1 - ( x0 * pk - pkm1 ) / k;
      pkm1 = pk;
      pk = pkp1;
    end

    d1 = order * ( pkm1 - x0 * pk );

    dpn = d1 / ( 1.0 - x0 * x0 );

    d2pn = ( 2.0 * x0 * dpn - e1 * pk ) / ( 1.0 - x0 * x0 );

    d3pn = ( 4.0 * x0 * d2pn + ( 2.0 - e1 ) * dpn ) / ( 1.0 - x0 * x0 );

    d4pn = ( 6.0 * x0 * d3pn + ( 6.0 - e1 ) * d2pn ) / ( 1.0 - x0 * x0 );

    u = pk / dpn;
    v = d2pn / dpn;
%
%  Initial approximation H:
%
    h = - u * ( 1.0 + 0.5 * u * ( v + u * ( v * v - d3pn / ( 3.0 * dpn ) ) ) );
%
%  Refine H using one step of Newton's method:
%
    p = pk + h * ( dpn + 0.5 * h * ( d2pn + h / 3.0 ...
      * ( d3pn + 0.25 * h * d4pn ) ) );

    dp = dpn + h * ( d2pn + 0.5 * h * ( d3pn + h * d4pn / 3.0 ) );

    h = h - p / dp;

    xtemp = x0 + h;

    xtab(mp1mi) = xtemp;

    fx = d1 - h * e1 * ( pk + 0.5 * h * ( dpn + h / 3.0 ...
      * ( d2pn + 0.25 * h * ( d3pn + 0.2 * h * d4pn ) ) ) );

    weight(mp1mi) = 2.0 * ( 1.0 - xtemp * xtemp ) / fx / fx;

  end

  if ( mod ( order, 2 ) == 1 )
    xtab(1) = 0.0;
  end
%
%  Shift the data up.
%
  nmove = floor ( ( order + 1 ) / 2 );
  ncopy = order - nmove;

  for i = 1 : nmove
    iback = order + 1 - i;
    xtab(iback) = xtab(iback-ncopy);
    weight(iback) = weight(iback-ncopy);
  end
%
%  Reflect values for the negative abscissas.
%
  for i = 1 : order - nmove
    xtab(i) = - xtab(order+1-i);
    weight(i) = weight(order+1-i);
  end

  return
end
function phy = reference_to_physical_t3 ( t, n, ref )

%*****************************************************************************80
%
%% REFERENCE_TO_PHYSICAL_T3 maps a reference point to a physical point.
%
%  Discussion:
%
%    Given the vertices of an order 3 physical triangle and a point
%    (XSI,ETA) in the reference triangle, the routine computes the value
%    of the corresponding image point (X,Y) in physical space.
%
%    Note that this routine may also be appropriate for an order 6
%    triangle, if the mapping between reference and physical space
%    is linear.  This implies, in particular, that the sides of the
%    image triangle are straight and that the "midside" nodes in the
%    physical triangle are halfway along the sides of
%    the physical triangle.
%
%    The T3 reference element is suggested by the following diagram:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the coordinates of the vertices.  The vertices are assumed
%    to be the images of (0,0), (1,0) and (0,1) respectively.
%
%    Input, integer N, the number of points to transform.
%
%    Input, real REF(2,N), the coordinates of points in the reference space.
%
%    Output, real PHY(2,N), the coordinates of the corresponding points in the
%    physical space.
%
  phy = zeros ( 2, n );

  for i = 1 : 2

    phy(i,1:n) = t(i,1) * ( 1.0 - ref(1,1:n) - ref(2,1:n) ) ...
               + t(i,2) *        ref(1,1:n)                ...
               + t(i,3) *                     ref(2,1:n);
  end

  return
end
function [ a, more, h, t, n2, more2 ] = subcomp_next ( n, k, a, more, h, t, ...
  n2, more2 )

%*****************************************************************************80
%
%% SUBCOMP_NEXT computes the next subcomposition of N into K parts.
%
%  Discussion:
%
%    A composition of the integer N into K parts is an ordered sequence
%    of K nonnegative integers which sum to a value of N.
%
%    A subcomposition of the integer N into K parts is a composition
%    of M into K parts, where 0 <= M <= N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer whose subcompositions are desired.
%
%    Input, integer K, the number of parts in the subcomposition.
%
%    Input, integer A(K), the parts of the subcomposition.
%
%    Input, logical MORE, set to FALSE by the user to start the computation.
%
%    Input, integer H, T, N2, MORE2, internal parameters needed for the
%    computation.  The user may need to initialize these before the
%    very first call, but these initial values are not important.
%    The user should not alter these parameters once the computation
%    begins.
%
%    Output, integer A(K), the parts of the subcomposition.
%
%    Output, logical MORE, set to FALSE by the routine to terminate
%    the computation.
%
%    Output, integer H, T, N2, MORE2, the updated values of the two internal
%    parameters.
%

%
%  The first computation.
%
  if ( ~more )

    n2 = 0;
    a(1:k) = 0;
    more2 = 0;
    h = 0;
    t = 0;

    more = 1;
%
%  Do the next element at the current value of N.
%
  elseif ( more2 )

    [ a, more2, h, t ] = comp_next ( n2, k, a, more2, h, t );

  else

    more2 = 0;
    n2 = n2 + 1;

    [ a, more2, h, t ] = comp_next ( n2, k, a, more2, h, t );

  end
%
%  Termination occurs if MORE2 = FALSE and N2 = N.
%
  if ( ~more2 && n2 == n )
    more = 0;
  end

  return
end
function area = triangle_area_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA_2D computes the area of a triangle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real AREA, the absolute area of the triangle.
%
  area = 0.5 * abs ( ...
      t(1,1) * ( t(2,2) - t(2,3) ) ...
    + t(1,2) * ( t(2,3) - t(2,1) ) ...
    + t(1,3) * ( t(2,1) - t(2,2) ) );

  return
end
function [ x, w ] = triangle_rule_13 ( )

%*****************************************************************************80
%
%% TRIANGLE_RULE_13 sets a 13 point quadrature rule for triangles.
%
%  Discussion:
%
%    The Integration region is points (X,Y) such that
%
%      0 <= X,
%      0 <= Y, and
%      X + Y <= 1.
%
%    Graph:
%
%      ^
%    1 | *
%      | |\
%    Y | | \
%      | |  \
%    0 | *---*
%      +------->
%        0 X 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gilbert Strang, George Fix,
%    An Analysis of the Finite Element Method,
%    Prentice Hall, 1973, page 184,
%    ISBN: 096140888X,
%    LC: TA335.S77.
%
%  Parameters:
%
%    Output, real X(2,13), the abscissas.
%
%    Output, real W(13), the weights.
%
  a = 0.479308067841923;
  b = 0.260345966079038;
  c = 0.869739794195568;
  d = 0.065130102902216;
  e = 0.638444188569809;
  f = 0.312865496004875;
  g = 0.048690315425316;
  h = 1.0 / 3.0;
  t = 0.175615257433204;
  u = 0.053347235608839;
  v = 0.077113760890257;
  w = -0.149570044467670;

  x(1:2,1:13) = [ a, b, b, c, d, d, e, e, f, f, g, g, h;
                  b, a, b, d, c, d, f, g, e, g, e, f, h ];

  w(1:13) = [ t, t, t, u, u, u, v, v, v, v, v, v, w ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function autoShowClassMap(varargin)
% This function is used to provde all land classification maps for each year
% Results for LCMAP
% Version 1.2  Acquire map from annual classfication map (04/09/2015)
% Version 1.1  Add disturbance in the annual map (04/01/2015);
% Version 1.0  No disturbance class in the cover map (11/06/2015)
% Tools
%
%
% Specific parameters
% ------------------------
%   'CCDCDir'     Directory of input data.  Default is the path to
%                        the current folder.
%   'StartYear'          Start year.
%   'EndYear'            End year.
%   'Year'               Year. Only a certain date's map will be outputted
%                        when inputing this.
%   'Month'              7
%   'Day'                1
%
%
% Examples:
% ------------------------
% autoShowClassMap('StartYear',2010,'EndYear',2015)
% will give the classifiaction maps at July 1 between 2010 and 2015.
%
% autoShowClassMap('StartYear',2000,'EndYear',2015,'Month',12,'Day',1)
% will give the classifiaction maps at December 1 between 2010 and 2015.
%
% autoShowClassMap('Year',2000,'Month',7,'Day',1)
% will give the classifiaction map at July 1, 2000.
%
   %% get parameters from inputs
    % where the all Landsat zipped files are
    dir_cur = pwd;
    p = inputParser;
    p.FunctionName = 'paras';

    addParameter(p,'CCDCDir',dir_cur);

    % default date is 1th, July.
    addParameter(p,'StartYear',0);
    addParameter(p,'EndYear',0);
    addParameter(p,'Year',0);
    addParameter(p,'Month',7);
    addParameter(p,'Day',1);

    % request user's input
    parse(p,varargin{:});
    dir_cur=p.Results.CCDCDir;

    start_year=p.Results.StartYear;
    end_year=p.Results.EndYear;
    oneyear=p.Results.Year;

    mm=p.Results.Month;
    dd=p.Results.Day;

    % give certain year
    if oneyear>0
        start_year = oneyear;
        end_year = oneyear;
    end

    % optional
%     addpath('~/ccdc');

    cd(dir_cur);

    % INPUTS:
    all_yrs = start_year:end_year; % all of years for producing maps
    v_input = ccdc_Inputs;

    % dimension and projection of the image
    nrows = v_input.ijdim(1);
    ncols = v_input.ijdim(2);
    jiDim = [ncols,nrows];
    jiUL = v_input.jiul;
    res = v_input.resolu;
    zc = v_input.zc;
    % number of coefficients
    num_c = v_input.num_c;
    % number of bands
    nbands = v_input.nbands;
    % max number of maps
    max_n = length(all_yrs);
    % all julidan dates
%     mm = 7;
%     dd = 1;
    jul_d = datenummx(all_yrs,mm,dd);
    jul_start = datenummx(all_yrs,1,1);
    jul_end = datenummx(all_yrs,12,31);

    % produce land cover map
    CoverMap = 255*ones(nrows,ncols,max_n,'uint8'); % Trends categories (0~11)
    % produce land cover map QA (unsupervised emsemble margin)
    CoverQAMap = 255*ones(nrows,ncols,max_n,'uint8'); % Trends categories (0~100)

    % make Predict folder for storing predict images
    n_map=v_input.name_map;% 'CCDCMap';
    if isempty(dir(n_map))
        mkdir(n_map);
    end

    % cd to the folder for storing recored structure
    cd(v_input.name_rst);

    imf = dir(fullfile(pwd,'record_change*')); % folder names
    num_line = size(imf,1);
    line_pt = 0;
    %%
    for line = 1:num_line
        if 100*(line/num_line) - line_pt > 1
            fprintf('Processing %.0f percent\n',ceil(100*(line/num_line)));
            line_pt = 100*(line/num_line);
        end

        % load one line of time series models
        load(fullfile(pwd,imf(line).name));

        % postions
        pos = [rec_cg.pos];

        % continue if there is no model available
        l_pos = length(pos);
        if l_pos == 0
            continue
        end

        % matrix of each component
        % start time
        t_start = [rec_cg.t_start];
        % end time
        t_end = [rec_cg.t_end];
        % break time
        t_break = [rec_cg.t_break];
        % class
        class = [rec_cg.class];
        % classification QA
        class_qa = [rec_cg.classQA];

% %         % put disturbed class to grass/shrub
% %         class(class == 10) = 7;

        % For Condition & Cover Map
        for i = 1:l_pos
            % get row and col
            [I,J] = ind2sub(jiDim,pos(i));

            % initialize pixels have at least one model
            if sum(CoverMap(J,I,:) == 255) == max_n
                % write land cover to CoverMap
                CoverMap(J,I,:) = 0;
                % write land cover to CoverQAMap
                CoverQAMap(J,I,:) = 0;
            end

            % give disturbed class for each year
            % year (band) the curve belongs to
            n_band = jul_d >= t_start(i) & (jul_d <= t_end(i) | jul_d < t_break(i));
            % write land cover to CoverMap
            CoverMap(J,I,n_band) = class(i);
            % write land cover to CoverQAMap
            CoverQAMap(J,I,n_band) = class_qa(i);

            % give next land cover category for gaps
            if i > 1
                if pos(i) == pos(i-1) % same location
                    n_dist = jul_d < t_start(i) & jul_d >= t_break(i-1);
                    % find the first nonzero value from class(:,i);
                    idn = find(class(:,i)>0);
                    if ~isempty(idn)
                        % write land cover (distubed) to CoverMap
                        CoverMap(J,I,n_dist) = class(idn(1),i);
                    end
                end
            end
        end
    end
    cd ..

    enviwrite_bands(fullfile(v_input.l_dir,n_map,'CoverMap'),CoverMap,'uint8',res,jiUL,'bsq',zc,all_yrs);
    clear CoverMap;
    enviwrite_bands(fullfile(v_input.l_dir,n_map,'CoverQAMap'),CoverQAMap,'uint8',res,jiUL,'bsq',zc,all_yrs);
    clear CoverQAMap;

% %     % write ENVI files for ARD
% %     % Cover Map
% %     ARD_enviwrite_bands([v_input.l_dir,'/',n_map,'/CoverMap'],CoverMap,'uint8','bsq',all_yrs);
% %     clear CoverMap
% %     ARD_enviwrite_bands([v_input.l_dir,'/',n_map,'/CoverQAMap'],CoverQAMap,'uint8','bsq',all_yrs);
% %     clear CoverQAMap

end
% %% Add change event
% change_map = double(enviread('CCDCMap/ChangeMap'));
% % minimum mapping unit
% num_obj = 5;
%
% % spatial filtering
% for i = 1:size(change_map,3)
%     i
%     tmp = change_map(:,:,i);
%     tmp(tmp > 0) = 1;
%     segm_tmp=bwlabeln(tmp,8);
%     L = segm_tmp;
%     s = regionprops(L,'area');
%     area = [s.Area];
%
%     % filter out cloud object < than num_cldoj pixels
%     idx = find(area >= num_obj);
%     tmp(ismember(L,idx)==0) = 0;
%     copy_map = change_map(:,:,i);
%     copy_map(tmp==0) = 0;
%     change_map(:,:,i) = copy_map;
% end
% % write filtered change map
% ARD_enviwrite_bands([v_input.l_dir,'/',n_map,'/FChangeMap'],change_map,'uint16','bsq',all_yrs,'example_img');
%
% %% change analysis
% max_n = 30;
% map = double(enviread('CCDCMap/CoverMap1_4'));
% change_stat = [];
%
% for i = 1:max_n-1
%     i
%     change_map = map(:,:,i)*100 + map(:,:,i+1);
%     idstb = map(:,:,i+1)-map(:,:,i);
%     change_stat = [change_stat;change_map(idstb~=0)];
% end
%
% % update class number
% all_class = unique(change_stat);
% % update number of class
% n_class = length(all_class);
% % calculate proportion based # for each class
% number = hist(change_stat,all_class);
% prct = number/sum(number);

% %% add disturbance (3) to the classification map
% for i = 1:max_n
%     change_tmp = change_map(:,:,i);
%     class_tmp = cover_map(:,:,i);
%     class_tmp(change_tmp>0) = 3;
%     cover_map(:,:,i) = class_tmp;
% end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
 function ob = Gdelaysum1(varargin)
%function ob = Gdelaysum1([options])
%|
%| Construct Gdelaysum1 object that performs weighted sums of a delayed signal.
%| This is useful for model-based reconstruction of THz images.
%| See Gdelaysum1_test() at end for example usage.
%|
%| y[n;m] = sum_{k=0}^{Nx-1} h[n - d[k;m]] x[k], n=0,...,Ny-1, m=0,...,Nm-1
%| for 0 <= n - d[k] <= Nh - 1
%|
%| required
%|	'Ny'	[1]		output signal length
%|	'delay'	[Nx Nm]		delays
%|
%| options
%|	'h'	[Nh 1]		impulse response, default [1]
%|	'nthread' [1]		# threads, default 1
%|
%| out
%|	ob	[Ny*Nm Nx]	Fatrix object
%|
%| Copyright 2006-8-25, Jeff Fessler, University of Michigan

if nargin == 1 && streq(varargin{1}, 'test'), Gdelaysum1_test, return, end
if nargin < 4, ir_usage, end

% defaults
arg.Ny = [];
arg.delay = [];
arg.h = 1;
arg.nthread = 1;
arg.class = 'fatrix2';
arg.chat = false;

% options
arg = vararg_pair(arg, varargin);

if isempty(arg.delay), fail 'delay required', end
if isempty(arg.Ny), fail 'Ny required', end

arg.h = single(arg.h);
arg.delay = single(arg.delay);
arg.nthread = int32(arg.nthread);
arg.chat = int32(arg.chat);

[arg.Nx arg.Nm] = size(arg.delay);
arg.Nh = length(arg.h);

arg.str_forw_mex = 'delaysum1,forw';
arg.str_back_mex = 'delaysum1,back';
if arg.nthread > 1
	arg.str_forw_mex = 'delaysum1,forw,thr';
	arg.str_back_mex = 'delaysum1,back,thr';
end

switch arg.class
case 'Fatrix'
	arg.dim = [arg.Ny * arg.Nm, arg.Nx];
	ob = Fatrix(arg.dim, arg, ...
		'abs', @Gdelaysum1_abs, 'power', @Gdelaysum1_power, ...
		'forw', @Gdelaysum1_forw_Fatrix, ...
		'back', @Gdelaysum1_back_Fatrix);
case 'fatrix2'

	odim = arg.Ny*arg.Nm;
	forw = @(arg, x) delaysum1_mex(arg.str_forw_mex, arg.h, arg.delay, ...
		arg.nthread, single(x), int32(arg.Ny), arg.chat);
	back = @(arg, y) delaysum1_mex(arg.str_back_mex, arg.h, arg.delay, ...
		arg.nthread, single(y), arg.chat);
	ob = fatrix2('idim', arg.Nx, 'odim', odim, 'arg', arg, ...
		'abs', @Gdelaysum1_abs, 'power', @Gdelaysum1_power, ...
		'forw', forw, 'back', back);
otherwise
	fail 'bug'
end

% Gdelaysum1_abs(): |A|
function ob = Gdelaysum1_abs(ob)
ob.arg.h = abs(ob.arg.h);

% Gdelaysum1_power(): A .^ p
function ob = Gdelaysum1_power(ob, p)
ob.arg.h = ob.arg.h .^ p;

% Gdelaysum1_forw_Fatrix(): y = A * x
% in
%	x	[Nx L]
% out
%	y	[Ny*Nm L]
%
function y = Gdelaysum1_forw_Fatrix(arg, x)

LL = size(x, 2);
y = zeros(arg.Ny*arg.Nm, LL, 'single');
for ll=1:LL
	tmp = single(x(:,ll));
	tmp = delaysum1_mex(arg.str_forw_mex, arg.h, arg.delay, ...
		arg.nthread, tmp, int32(arg.Ny), arg.chat);
	y(:,ll) = tmp(:);
end

% Gdelaysum1_back_Fatrix(): x = A' * y
% in
%	y	[Ny*Nm L]
% out
%	x	[Nx L]
%
function x = Gdelaysum1_back_Fatrix(arg, y)

LL = size(y, 2);
x = zeros(arg.Nx, LL, 'single');
for ll=1:LL
	tmp = single(y(:,ll));
	tmp = delaysum1_mex(art.str_back_mex, arg.h, arg.delay, ...
		arg.nthread, tmp, arg.chat);
	x(:,ll) = tmp;
end

% Gdelaysum1_test
function Gdelaysum1_test

delay = [10 20 40; 50 51 52];
h = [4 3 -2 1]';
Ny = 80;
ats = {'nthread', jf('ncore')};
A = Gdelaysum1('Ny', Ny, 'h', h, 'delay', delay, ats{:});

x = [1 0]';
y1 = A * x;
y1 = reshapee(y1, Ny, []);
if im
	clf, plot(y1, '-o')
end

%mask = true(A.arg.Nx,1);
%Fatrix_test_basic(A, mask)
fatrix2_tests(A)
test_adjoint(A);

if 1 % compare Fatrix and fatrix2
	Ny = 2^11;
	Nx = 2^10;
	Nm = 2^9;
	Nh = 2^8;
	rng(0)
	delay = rand(Nx, Nm);
	x = rand(Nx, 1);
	args = {'Ny', Ny, 'h', (1:Nh)', 'delay', delay, ats{:}};
	A1 = Gdelaysum1(args{:}, 'class', 'Fatrix');
	A2 = Gdelaysum1(args{:}, 'class', 'fatrix2');
	cpu etic
	y1 = A1 * x;
	cpu etoc 'Fatrix'
	cpu etic
	y2 = A2 * x;
	cpu etoc 'fatrix2'
	jf_equal(y1, y2)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function setpaths(addremove)

%
% USAGE:
%
%   % 1. Add all search paths for this repo
%   setpaths
%
%   % 2. Same as 1.. Add all search paths for this repo
%   setpaths(1)
%
%   % 3. Add all search paths for this repo while removing search
%   %    paths of all similar workspaces
%   setpaths(2)
%
%   % 4. Remove all search paths for this repo
%   setpaths(0)
%
currdir = pwd;
global logger

try

    warning('off','MATLAB:rmpath:DirNotFound');

    appname = 'Homer3';

    % Parse arguments
    if ~exist('addremove','var')
    addremove = 1;
    end

    % Add libraries on which Homer3 depends
    d = addDependenciesSearchPaths();

    % Start logger only after adding library paths. Logger is in the Utils libary.
    logger = Logger('setpaths');

    % Create list of possible known similar apps that may conflic with current
    % app
    appNameExclList = {'Homer3','DataTree','Homer2_UI','brainScape','ResolveCommonFunctions'};
    appNameInclList = {'AtlasViewerGUI'};
    exclSearchList  = {'.git','.idea','Data','Docs','*_install','*.app','submodules'};

    appThis         = filesepStandard_startup(pwd);
    appThisPaths    = findDotMFolders(appThis, exclSearchList);

    if addremove == 0
        if ~isempty(which('deleteNamespace.m'))
            deleteNamespace(appname);
        end
        removeSearchPaths(appThis);
        return;
    elseif addremove == 2
        appNameExclList = [appNameExclList, appNameInclList];
        appNameInclList = {};
    end

    appExclList = {};
    appInclList = {};

    % Find all root folders of apps to exclude from search paths
    for ii = 1:length(appNameExclList)
        foo = which([appNameExclList{ii}, '.m'],'-all');

        % Before giving up and concluding the path does not exist in list of matlab
        % search paths, try appending the word 'Class' to path name.
        if isempty(foo)
            foo = which([appNameExclList{ii}, 'Class.m'],'-all');
        end

        for jj = 1:length(foo)
            p = filesepStandard_startup(fileparts(foo{jj}));
            if pathscompare_startup(appThis, p)
                continue
            end
            printMethod(sprintf('Exclude paths for %s\n', p));
            appExclList = [appExclList; p]; %#ok<AGROW>
        end
    end

    % Find all root folders of apps to include in search paths
    for ii = 1:length(appNameInclList)
        foo = which([appNameInclList{ii}, '.m'],'-all');
        for jj = 1:length(foo)
            if jj > 1
                p = filesepStandard_startup(fileparts(foo{jj}));
                appExclList = [appExclList; p]; %#ok<AGROW>
                printMethod(sprintf('Exclude paths for %s\n', p));
            else
                p = filesepStandard_startup(fileparts(foo{jj}));
                appInclList = [appInclList; p]; %#ok<AGROW>
                printMethod(sprintf('Include paths for %s\n', p));
            end
        end
    end

    % Remove all search paths for all other apps except for current one, to
    % make sure that we use only search from the current app for download shared
    % libraries (i.e, submodules).
    for ii = 1:length(appExclList)
        removeSearchPaths(appExclList{ii})
    end
    for ii = 1:length(appInclList)
        removeSearchPaths(appInclList{ii})
    end

    addSearchPaths(appThisPaths);

    if ~isempty(which('setNamespace.m'))
        setNamespace(appname);
    end

    % Add back all search paths for all other apps except for current app
    for ii = 1:length(appInclList)
        % This app's path has already been added
        if pathscompare_startup(appInclList{ii}, appThis)
            continue;
        end
        foo = findDotMFolders(appInclList{ii}, exclSearchList);
        addSearchPaths(foo);
    end

    if  ~isempty(which('setpaths_proprietary.m'))
        setpaths_proprietary(addremove);
    end

    warning('on','MATLAB:rmpath:DirNotFound');

    PrintSystemInfo(logger, ['Homer3'; d]);
    logger.CurrTime('Setpaths completed on ');
    logger.Close();
    cd(currdir);

catch ME

    printStack(ME)
    if exist('logger','var')
        logger.Close();
    end
    cd(currdir);
    rethrow(ME)

end

% ---------------------------------------------------
function d = dependencies()
d = {};
submodules = parseGitSubmodulesFile(pwd);
temp = submodules(:,1);
for ii = 1:length(temp)
    [~, d{ii,1}] = fileparts(temp{ii});
end

% ---------------------------------------------------
function setpermissions(appPath)
if isunix() || ismac()
    global logger
    if ~isempty(strfind(appPath, '/bin')) %#ok<*STREMP>
        cmd = sprintf('chmod 755 %s/*\n', appPath);
        logger.Write(cmd);
        files = dir([appPath, '/*']);
        if ~isempty(files)
            system(cmd);
        end
    end
end

% ----------------------------------------------------
function addSearchPaths(appPaths)
if ischar(appPaths)
    p = genpath(appPaths);
    if ispc
        delimiter = ';';
    else
        delimiter = ':';
    end
    appPaths = str2cell(p, delimiter);
end
for kk = 1:length(appPaths)
    if strfind(appPaths{kk}, '.git')
        continue;
    end
    addpath(appPaths{kk}, '-end');
    setpermissions(appPaths{kk});
end
printMethod(sprintf('ADDED search paths for app %s\n', appPaths{1}));

% ----------------------------------------------------
function removeSearchPaths(app)
p = path;
if ispc()
    delimiter = ';';
elseif ismac() || isunix()
    delimiter = ':';
end
[~,appname] = fileparts(fileparts(app));
r = version('-release');
msg = sprintf('Removing search paths for %s ...', appname);
h = waitbar(0, msg);
p = str2cell_startup(p, delimiter);
for kk = 1:length(p)
    if mod(kk,100)==0
        waitbar(kk/length(p), h);
    end
    if ~isempty(strfind(lower(p{kk}), 'matlab')) && ~isempty(strfind(p{kk}, r))
        continue;
    end
    if ~isempty(strfind(filesepStandard_startup(p{kk}), app))
        rmpath(p{kk});
    end
end
close(h);
printMethod(sprintf('REMOVED search paths for app %s\n', app));

% ----------------------------------------------------
function   d = addDependenciesSearchPaths()
if exist([pwd, '/Utils/submodules'],'dir')
    addpath([pwd, '/Utils/submodules'],'-end');
end
d = dependencies();
for ii = 1:length(d)
    rootpath = findFolder(pwd, d{ii});
    rootpath(rootpath=='\') = '/';
    if ispathvalid_startup([rootpath, '/Shared'],'dir')
        rootpath = [rootpath, '/Shared'];
    end
    if ~exist(rootpath,'dir')
        printMethod(sprintf('ERROR: Could not find required dependency %s\n', d{ii}));
        continue;
    end
    addSearchPaths(rootpath);
end
if exist([pwd, '/Utils/submodules'],'dir')
    addpath([pwd, '/Utils/submodules'],'-end');
end

% -----------------------------------------------------------------------------
function [C,k] = str2cell(str, delimiters, options)

% Option tells weather to keep leading whitespaces.
% (Trailing whitespaces are always removed)
if ~exist('options','var')
    options = '';
end

if ~strcmpi(options, 'keepblanks')
    str = strtrim(str);
end
str = deblank(str);

if ~exist('delimiters','var') || isempty(delimiters)
    delimiters{1} = sprintf('\n');
elseif ~iscell(delimiters)
    foo{1} = delimiters;
    delimiters = foo;
end

% Get indices of all the delimiters
k=[];
for kk=1:length(delimiters)
    k = [k, find(str==delimiters{kk})];
end
j = find(~ismember(1:length(str),k));

% The following line seems to hurt performance a little bit. It was
% meant to preallocate to speed things up but it does not seem to do that.
% C = repmat({blanks(max(diff([k,length(str)])))}, length(k)+1, 1);
C = {};
ii=1; kk=1;
while ii<=length(j)
    C{kk} = str(j(ii));
    ii=ii+1;
    jj=2;
    while (ii<=length(j)) && ((j(ii)-j(ii-1))==1)
        C{kk}(jj) = str(j(ii));
        jj=jj+1;
        ii=ii+1;
    end
    C{kk}(jj:end)='';
    kk=kk+1;
end
C(kk:end) = [];

% -------------------------------------------------------------------------
function printMethod(msg)
global logger
if isa(logger', 'Logger')
    try
        logger.Write(msg);
    catch
        fprintf(msg);
    end
else
    fprintf(msg);
end

% -------------------------------------------------------------------------
function submodules = parseGitSubmodulesFile(repo)
submodules = cell(0,3);

if ~exist('repo','var') || isempty(repo)
    repo = pwd;
end
currdir = pwd;
if repo(end) ~= '/' && repo(end) ~= '\'
    repo = [repo, '/'];
end

filename = [repo, '.gitmodules'];
if ~exist(filename, 'file')
    return;
end
cd(repo);

fid = fopen(filename, 'rt');
strs = textscan(fid, '%s');
strs = strs{1};
kk = 1;
for ii = 1:length(strs)
    if strcmp(strs{ii}, '[submodule')
        jj = 1;
        while ~strcmp(strs{ii+jj}, '[submodule')
            if ii+jj+2>length(strs)
                break;
            end
            if strcmp(strs{ii+jj}, 'path')
                submodules{kk,2} = [pwd, '/', strs{ii+jj+2}];
            end
            if strcmp(strs{ii+jj}, 'path')
                submodules{kk,3} = strs{ii+jj+2};
            end
            if strcmp(strs{ii+jj}, 'url')
                submodules{kk,1} = strs{ii+jj+2};
            end
            jj = jj+1;
        end
        kk = kk+1;
    end
end
fclose(fid);
cd(currdir);

author: BUNPC


text:
function init_params = init_coeffs(x,y)
% INIT_COEFFS Function to generate the initial parameters for the 4
% parameter dose response curve.
% Requires an array of doses and an array of responses
% This function is used by sigmoid.m and ec50.m
%
% Copyright 2004 Carlos Evangelista
% send comments to CCEvangelista@aol.com
% Version 1.0    01/07/2004

parms=ones(1,4);
parms(1)=min(y);
parms(2)=max(y);
parms(3)=(min(x)+max(x))/2;
sizey=size(y);
sizex=size(x);
if (y(1)-y(sizey))./(x(2)-x(sizex))>0
    parms(4)=(y(1)-y(sizey))./(x(2)-x(sizex));
else
    parms(4)=1;
end
init_params=parms;


text:
function [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(varargin)
% [bnet, Qnodes, Fnodes, Onode] = mk_motif_hhmm(...)
%
% Make the following HHMM
%
%    S2 <----------------------> S1
%    |                           |
%    |                           |
%   M1 -> M2 -> M3 -> end        B1 -> end
%
% where Mi represents the i'th letter in the motif
% and B is the background state.
% Si chooses between running the motif or the background.
% The Si and B states have self loops (not shown).
%
% The transition params are defined to respect the above topology.
% The background is uniform; each motif state has a random obs. distribution.
%
% Optional params:
% motif_length  - required, unless we specify motif_pattern
% motif_pattern - if specified, we make the motif submodel deterministically
%                  emit this pattern
% background    - if specified, we make the background submodel
%                  deterministically emit this (makes the motif easier to see!)

args = varargin;
nargs = length(args);

% extract pattern, if any
motif_pattern = [];
for i=1:2:nargs
  switch args{i},
   case 'motif_pattern', motif_pattern = args{i+1};
  end
end

% set defaults
motif_length = length(motif_pattern);
background_char = [];

% get params
for i=1:2:nargs
  switch args{i},
   case 'motif_length', motif_length = args{i+1};
   case 'background', background_char = args{i+1};
  end
end

chars = ['a', 'c', 'g', 't'];
Osize = length(chars);

Qsize = [2 motif_length];
Qnodes = 1:2;
D = 2;
transprob = cell(1,D);
termprob = cell(1,D);
startprob = cell(1,D);

% startprob{d}(k,j), startprob{1}(1,j)
% transprob{d}(i,k,j), transprob{1}(i,j)
% termprob{d}(k,j)

% LEVEL 1

startprob{1} = zeros(1, 2);
startprob{1} = [1 0]; % always start in the background model

% When in the background state, we stay there with high prob
% When in the motif state, we immediately return to the background state.
transprob{1} = [0.8 0.2;
		1.0 0.0];

% LEVEL 2
startprob{2} = 'leftstart'; % both submodels start in substate 1
transprob{2} = zeros(motif_length, 2, motif_length);
termprob{2} = zeros(2, motif_length);

% In the background model, we only use state 1.
transprob{2}(1,1,1) = 1; % self loop
termprob{2}(1,1) = 0.2; % prob transition to end state

% Motif model
transprob{2}(:,2,:) = mk_leftright_transmat(motif_length, 0); % no self loops
termprob{2}(2,end) = 1.0; % last state immediately terminates

% OBS LEVEl

obsprob = zeros([Qsize Osize]);
if isempty(background_char)
  % uniform background model
  %obsprob(1,1,:) = normalise(ones(Osize,1));
  obsprob(1,1,:) = normalise(rand(Osize,1));
else
  % deterministic background model (easy to see!)
  m = find(chars==background_char);
  obsprob(1,1,m) = 1.0;
end

if ~isempty(motif_pattern)
  % initialise with true motif (cheating)
  for i=1:motif_length
    m = find(chars == motif_pattern(i));
    obsprob(2,i,m) = 1.0;
  end
else
  obsprob(2,:,:) = mk_stochastic(rand(motif_length, Osize));
end

if 0
  Oargs = {'CPT', obsprob};
else
  % We use a minent prior for the emission distribution for the states in the motif model
  % (but not the background model). This encourages nearly deterministic distributions.
  % We create an index matrix  (where M = motif length)
  %  [2 1
  %   2 2
  %   ...
  %   2 M]
  % and then convert this to a list of integers, which
  % specifies when to use the minent prior (Q1=2 specifies motif model).
  M = motif_length;
  ndx = [2*ones(M,1) (1:M)'];
  pcases = subv2ind([2 motif_length], ndx);
  Oargs = {'CPT', obsprob, 'prior_type', 'entropic', 'entropic_pcases', pcases};
end

[bnet, Qnodes, Fnodes, Onode] = mk_hhmm('Qsizes', Qsize, 'Osize', Osize, 'discrete_obs', 1, ...
	       'Oargs', Oargs, 'Ops', Qnodes(1:2), ...
	       'startprob', startprob, 'transprob', transprob, 'termprob', termprob);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [disp_row, disp_col, scale_ind] = optimize_scores(scores_fs, iterations)

% Maximizes the continuous convolution response (classification scores).
% Find the size of the output.
[sz1, sz2, num_scales] = size(scores_fs);
output_sz = [sz1 sz2];

% Do the grid search step by finding the maximum in the sampled response
% for each scale.
sampled_scores = sample_fs(scores_fs);
[max_resp_row, max_row] = max(sampled_scores, [], 1);
[init_max_score, max_col] = max(max_resp_row, [], 2);
max_row_perm = permute(max_row, [2 3 1]);
col = max_col(:)';
row = max_row_perm(sub2ind(size(max_row_perm), col, 1:size(sampled_scores,3)));

% Shift and rescale the coordinate system to [-pi, pi]
trans_row = mod(row - 1 + floor((output_sz(1)-1)/2), output_sz(1)) - floor((output_sz(1)-1)/2);
trans_col = mod(col - 1 + floor((output_sz(2)-1)/2), output_sz(2)) - floor((output_sz(2)-1)/2);
init_pos_y = permute(2*pi * trans_row / output_sz(1), [1 3 2]);
init_pos_x = permute(2*pi * trans_col / output_sz(2), [1 3 2]);

% Set the current maximum to the sampled one
max_pos_y = init_pos_y;
max_pos_x = init_pos_x;

% construct grid
ky = -ceil((output_sz(1) - 1)/2) : floor((output_sz(1) - 1)/2);
kx = (-ceil((output_sz(2) - 1)/2) : floor((output_sz(2) - 1)/2))';

% pre-compute complex exponential
exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

ky2 = ky.*ky;
kx2 = kx.*kx;

iter = 1;
while iter <= iterations
    % Compute gradient
    ky_exp_ky = bsxfun(@times, ky, exp_iky);
    kx_exp_kx = bsxfun(@times, kx, exp_ikx);
    y_resp = mtimesx(exp_iky, scores_fs, 'speed');
    resp_x = mtimesx(scores_fs, exp_ikx, 'speed');
    grad_y = -imag(mtimesx(ky_exp_ky, resp_x, 'speed'));
    grad_x = -imag(mtimesx(y_resp, kx_exp_kx, 'speed'));

    % Compute Hessian
    ival = 1i * mtimesx(exp_iky, resp_x, 'speed');
    H_yy = real(-mtimesx(bsxfun(@times, ky2, exp_iky), resp_x, 'speed') + ival);
    H_xx = real(-mtimesx(y_resp, bsxfun(@times, kx2, exp_ikx), 'speed') + ival);
    H_xy = real(-mtimesx(ky_exp_ky, mtimesx(scores_fs, kx_exp_kx, 'speed'), 'speed'));
    det_H = H_yy .* H_xx - H_xy .* H_xy;

    % Compute new position using newtons method
    max_pos_y = max_pos_y - (H_xx .* grad_y - H_xy .* grad_x) ./ det_H;
    max_pos_x = max_pos_x - (H_yy .* grad_x - H_xy .* grad_y) ./ det_H;

    % Evaluate maximum
    exp_iky = exp(bsxfun(@times, 1i * ky, max_pos_y));
    exp_ikx = exp(bsxfun(@times, 1i * kx, max_pos_x));

    iter = iter + 1;
end

% Evaluate the Fourier series at the estimated locations to find the
% corresponding scores.
max_score = real(mtimesx(mtimesx(exp_iky, scores_fs, 'speed'), exp_ikx, 'speed'));

% check for scales that have not increased in score
ind = max_score < init_max_score;
max_score(ind) = init_max_score(ind);
max_pos_y(ind) = init_pos_y(ind);
max_pos_x(ind) = init_pos_x(ind);

% Find the scale with the maximum response
[max_scale_response, scale_ind] = max(max_score(:));

% Scale the coordinate system to output_sz
disp_row = (mod(max_pos_y(1,1,scale_ind) + pi, 2*pi) - pi) / (2*pi) * output_sz(1);
disp_col = (mod(max_pos_x(1,1,scale_ind) + pi, 2*pi) - pi) / (2*pi) * output_sz(2);
end

repo: MFT


text:
function chooseBestModels_HN_OldWithExp(pathExperimentsBinary,pathExperimentsTime,nExp,fSetNames,nameOutcomes,metric,maxOrder)
% -------------------------------------------------------------------------
% function chooseBestModels_HN(pathWORK,nExp,fSetName,nameOutcomes,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function requires user input to choose the set types and model orders
% providing the best parsimonious models for all outcomes analyzed in the
% HN study. See ref.[1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to the directory containing all experiments.
%                     --> Ex: '/myProject/WORKSPACE/CV-BASED_RESULTS'
% 2. nExp: Numerical value specifying the number of experiments to analyze.
%          --> Ex: 10
% 3. fSetName: Cell of strings specifying the name of the type of feature
%              set analyzed.
%              --> Ex: {'PET','CT','SEPARATE','FUSED'}
% 4. nameOutcomes: Cell of strings specifying the outcome names to analyze.
%                  --> Ex: {'Failure','Locoregional','Distant','Death'}
% 5. maxOrder: Integer specifying the maximal multivariable model order.
%              --> Ex: 10
% -------------------------------------------------------------------------
% OUTPUTS: Final prediction models saved in a folder named 'FINAL_MODELS/nameOutcome/fSetName'
%          in the given experiment folder.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathExperimentsBinary)
nOutcomes = length(nameOutcomes);
nType = length(fSetNames);

fprintf('\n')
for exp = 1:nExp
    cd(pathExperimentsTime), pathExperimentTime = fullfile(pathExperimentsTime,['Experiment',num2str(exp)]);
    cd(pathExperimentTime), mkdir('FINAL_MODELS'), cd('FINAL_MODELS'), pathFinalModelsTime = pwd;
    cd(pathExperimentsBinary), pathExperimentBinary = fullfile(pathExperimentsBinary,['Experiment',num2str(exp)]);
    cd(pathExperimentBinary), mkdir('FINAL_MODELS'), cd('FINAL_MODELS'), pathFinalModelsBinary = pwd;
    for o = 1:nOutcomes
        cd(pathFinalModelsTime), mkdir(nameOutcomes{o}), cd(fullfile(pathFinalModelsTime,nameOutcomes{o})), pathOutcomeTime = pwd;
        cd(pathFinalModelsBinary), mkdir(nameOutcomes{o}), cd(fullfile(pathFinalModelsBinary,nameOutcomes{o})), pathOutcomeBinary = pwd;
        plotPredictionResults_HN(fullfile(pathExperimentBinary,'RESULTS'),nameOutcomes{o},fSetNames,metric,maxOrder)
        fprintf('\n======== DISPLAYING PREDICTION RESULTS FOR "EXPERIMENT %u", "%s" OUTCOME, "%s" METRIC  ========\n',exp,nameOutcomes{o},metric)
        for f = 1:nType
            cd(pathOutcomeTime), mkdir(fSetNames{f}), cd(fullfile(pathOutcomeTime,fSetNames{f})), pathFinalTime = pwd;
            cd(pathOutcomeBinary), mkdir(fSetNames{f}), cd(fullfile(pathOutcomeBinary,fSetNames{f})), pathFinalBinary = pwd;
            while 1
                order = input(['Which model order of the ',fSetNames{f},' feature set provides the best parsimonious model? \n' ...
                               '--> Type a number between 1 to ',num2str(maxOrder),' and press ENTER \n' ...
                               'ANSWER: ']);
                fprintf('\n')
                if isnumeric(order) && order <= 10 && order >= 1
                    break
                end
            end
            cd(fullfile(pathExperimentBinary,'RESULTS'))
            results = load(['RESULTS_',fSetNames{f},'_',nameOutcomes{o}]); results = struct2cell(results); results = results{1};
            finalModel = results.(['Order',num2str(order)]);
            finalModel.Order = order;
            finalModel.outcome = nameOutcomes{o};
            cd(pathFinalBinary), save('finalModel','finalModel')
            cd(pathFinalTime), save('finalModel','finalModel')
        end
        close all
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function mov = make3Davi(varargin)
% ::
%
%    mov = make3Davi([opt] Options_Structure)
%
% Makes an avi movie file called head3d[x].avi or whatever you specify.
%
% :Options:
%
%   **O.name:**
%        'output_name.avi';
%
%   **O.fps:**
%        frames per second
%
%   **O.length:**
%        length of movie in s
%
%   **O.H:**
%        cluster handles in vector   (for adjusting transparency with time)
%
%   **O.timecourse{i}:**
%        cell array of time courses for each cluster
%
%   **O.timeres:**
%        resolution, in s, of timecourse data
%
%   **O.azOffset:**
%        azimuth value to offset, positive = move clockwise
%
%   **O.elOffset:**
%        elevation value to move through, positive = inf to superior
%
%   **O.zoom:**
%        zoom value to end up with
%
%   **O.add2movie:**
%        add to existing movie - enter mov structure in this field
%
%   **O.closemovie:**
%        1 or 0, close the movie afterward or not.
%
% Blank fields for az, el, zoom, timecourse indicate that these functions should not be performed
% This script spirals up, right, and in 36 degrees
%
% Notes: my indeo5 one wouldn't work in media player
%        also only seems to work if you add clusters before head isosurfaces
%        for Mac OS X / UNIX, choose "no compression"
%
% start with the image in the location you want to zoom in on
% but with no zoom.
%
% The movie FINISHES at the current axis locations
%
% :Example:  "Surface tour"
% ::
%
%    O = struct('name','all4_union_bg_surf.avi','fps',10,'length',6,'azOffset',180,'zoom',1,'elOffset',90,'closemovie',0); mov = make3davi(O);
%    view(90,0);lightfollowview;O = struct('name','all4_union_bg_surf.avi','add2movie',mov,'fps',10,'length',6,'azOffset',90,'zoom',1,'elOffset',-90,'closemovie',0); mov = make3davi(O);
%    view(0,90);lightfollowview;O = struct('name','all4_union_bg_surf.avi','add2movie',mov,'fps',10,'length',6,'azOffset',-90,'zoom',1,'elOffset',90,'closemovie',0); mov = make3davi(O);
%    mov = close(mov);
%
% ..
%    By Tor Wager, 10/07/01
% ..

% ..
%    set default values and inputs
% ..
myName = 'head3d1.avi'; i = 1;                      % set name to first unused file #
while exist(myName) == 2, myName = ['head3d' num2str(i) '.avi'];, i = i+1;,end
O.name = myName;
O.fps = 10;
O.length = 8;
O.H = [];
O.timecourse = [];
O.timeres = 1;
O.azOffset = 36;                                    % positive starts from left, moves r
                                                    % spins brain clockwise
O.elOffset = 20;                                    % positive starts from below, moves up
O.zoom = 2;                                         % final image will be this many times closer than at start
O.add2movie = [];                                   % existing movie structure to add to.  [] = new movie file.
O.closemovie = 1;                                   % close movie after finishing.  0 to add more frames later.

if nargin > 0
    inO = varargin{1};
    if isfield(inO,'name'), O.name = inO.name;,end
    if isfield(inO,'fps'), O.fps = inO.fps;,end
    if isfield(inO,'length'), O.length = inO.length;,end
    if isfield(inO,'H'), O.H = inO.H;,end
    if isfield(inO,'timecourse'), O.timecourse = inO.timecourse;,end
    if isfield(inO,'timeres'), O.timeres = inO.timeres;,end
    if isfield(inO,'azOffset'), O.azOffset = inO.azOffset;,end
    if isfield(inO,'elOffset'), O.elOffset = inO.elOffset;,end
    if isfield(inO,'zoom'), O.zoom = inO.zoom;,end
    if isfield(inO,'add2movie'), O.add2movie = inO.add2movie;,end
    if isfield(inO,'closemovie'), O.closemovie = inO.closemovie;,end
end

numFrames = O.fps .* O.length;
interpValue = O.fps ./ O.timeres;

% -------------------------------------------
% * check values to make sure everything's ok
% -------------------------------------------

if ~isempty(O.timecourse) & length(O.H) ~= length(O.timecourse), error('Number of handles and timecourses must be equal.'),end
if interpValue ~= round(interpValue), error('Frames per second must be a multiple of time res - adjust either one.'),end
if any(~ishandle(O.H)),
    warning('At least one cluster handle is not a valid object. Searching for handles by Tag name cluster[x]')
    for i = 1:length(O.H)
        O.H(i) = findobj('Tag',['cluster' num2str(i)]);
    end
    if any(~ishandle(O.H)), error('Cannot find valid cluster handles.'),end
end

% -------------------------------------------
% * create transparency values from timecourses
% -------------------------------------------
% normalize all timecourses to alpha values, from 0 to 1, by dividing by max value
if ~isempty(O.timecourse)
  for i = 1:length(O.timecourse),
    %O.timecourse{i} = O.timecourse{i} - min(O.timecourse{i});       % prevent negative values, set min to 0
    %min(O.timecourse{i})
    %O.timecourse{i} = O.timecourse{i} ./ max(abs(O.timecourse{i})); % max is one
    numSeconds = ceil(length(O.timecourse{i}) .* O.timeres);
    %numFrames = numSeconds .* O.fps;                               % adjust length
    alphaVal{i} = interp(O.timecourse{i},interpValue);              % interpolate to number of frames
    %alphaVal{i} = alphaVal{i} - min(alphaVal{i});                   % prevent negative values, set min to 0
    %alphaVal{i} = alphaVal{i} ./ max(abs(alphaVal{i}));             % max is one
    alphaVal{i} = [ alphaVal{i} zeros(1,max([0 numFrames - length(alphaVal{i})])) ];   % pad with zeros

    if any(alphaVal{i} < 0 | alphaVal{i} > 1),
        alphaVal{i},min(alphaVal{i}),max(alphaVal{i})
        error('Script error: alpha values must be between 0 and 1.')
    end
  end
end

% -------------------------------------------
% * set initial view and parameters
% -------------------------------------------
%set(gcf,'Color','w')
axis off
[az,el] = view;
if ~isempty(O.azOffset),az = az - O.azOffset;,az2add = O.azOffset ./ numFrames;,end
if ~isempty(O.elOffset),el = el - O.elOffset;,el2add = O.elOffset ./ numFrames;,end
view(az,el)
lightFollowView
if ~isempty(O.zoom), myZoom = 1 + ((O.zoom-1) ./ numFrames);, end

% -------------------------------------------
% * initialize movie file, if necessary
% -------------------------------------------

if isempty(O.add2movie)
    try
        %mov = avifile(O.name,'Quality',75,'Compression','Indeo5','Fps',O.fps);
        mov = avifile(O.name,'Quality',75,'Compression','none','Fps',O.fps);
    catch
        error('Cannot open movie file.  Probably already open: Try changing filename.')
    end
else
    mov = O.add2movie;
end

try

H = gca;

% -------------------------------------------
% * rotate, zoom, etc. and add frames
% -------------------------------------------
for i = 1:numFrames

	if ~isempty(O.azOffset),az = az + az2add;,end
    if ~isempty(O.elOffset),el = el + el2add;,end
	if ~isempty(O.azOffset) | ~isempty(O.elOffset)
        view(az,el);
        lightFollowView
    end
	if ~isempty(O.zoom),camzoom(myZoom);,end
    if ~isempty(O.timecourse)
        for j = 1:length(alphaVal)
            set(O.H(j),'FaceAlpha',alphaVal{j}(i))
        end
        myTime = i ./ O.fps;
        if myTime == round(myTime)
            title(['Time = ' num2str(i ./ O.fps) '.0 s'],'Color','w','FontSize',18)
        else
            title(['Time = ' num2str(i ./ O.fps) ' s'],'Color','w','FontSize',18)
        end

    end
    drawnow
	try
        mov = addframe(mov,H);
    catch
        error('Cannot write frame.  Failed to set stream format??')
    end

end

catch, lasterr, warning('There was an error making this movie file.  No frames written.')
    mov = close(mov);
end

% -------------------------------------------
% * close the movie
% -------------------------------------------

if O.closemovie, mov = close(mov);, end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

%%
clear; close all; clc;

%% Syntax
% |[Vm]=rigidTransformVertices(V,T,v);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [rms, nees] = compute_rms_nees_ave( rmsall, neesall )
% compute rms and nees from multiple tests

nsteps = size(rmsall.position, 2);
ntests = size(rmsall.position, 1);

% compute rms
rmspos_sum = sum(rmsall.position.^2, 1);
rms.position = sqrt(rmspos_sum/ntests);

rmsori_sum = sum(rmsall.orientation.^2, 1);
rms.orientation = sqrt(rmsori_sum/ntests);

% compute nees
neesall_pose = neesall.pose(:, 2:end);
neesall_orientation = neesall.orientation(:, 2:end);

neespose_sum = sum(neesall_pose, 1);
nees.pose = neespose_sum/ntests;

neesorientation_sum = sum(neesall_orientation, 1);
nees.orientation = neesorientation_sum/ntests;

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
classdef SetPowerSrcActiveStateAction < AbstractEventAction
    %SetPowerSrcActiveStateAction Summary of this class goes here
    %   Detailed explanation goes here

    properties
        src AbstractLaunchVehicleElectricalPowerSrcSnk
        activeStateToSet(1,1) logical = false;
    end

    properties(Constant)
        emptyVarArr = AbstractOptimizationVariable.empty(0,1);
    end

    methods
        function obj = SetPowerSrcActiveStateAction(src, activeStateToSet)
            if(nargin > 0)
                obj.src = src;
                obj.activeStateToSet = activeStateToSet;
            end

            obj.id = rand();
        end

        function newStateLogEntry = executeAction(obj, stateLogEntry)
            newStateLogEntry = stateLogEntry;
            stageStates = newStateLogEntry.stageStates;

            pwrSrcStates = AbstractLaunchVehicleElectricalPowerSrcState.empty(1,0);
            for(i=1:length(stageStates)) %#ok<*NO4LP>
                pwrSrcStates = horzcat(pwrSrcStates, stageStates(i).powerSrcStates); %#ok<AGROW>
            end

            pwrSrcs = AbstractLaunchVehicleElectricalPowerSrcSnk.empty(1,0);
            for(i=1:length(pwrSrcStates))
                pwrSrcs = horzcat(pwrSrcs, pwrSrcStates(i).getEpsSrcComponent()); %#ok<AGROW>
            end

            pwrSrcInd = find(pwrSrcs == obj.src,1,'first');
            pwrSrcState = pwrSrcStates(pwrSrcInd);

            pwrSrcState.setActiveState(obj.activeStateToSet);
        end

        function initAction(obj, initialStateLogEntry)
            %nothing
        end

        function name = getName(obj)
            if(obj.activeStateToSet)
                tf = 'Active';
            else
                tf = 'Inactive';
            end

            name = sprintf('Set Power Source State (%s = %s)',obj.src.getName(),tf);
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;
        end

        function tf = usesPwrSink(obj, powerSink)
            tf = false;
        end

        function tf = usesPwrSrc(obj, powerSrc)
            tf = [obj.src] == powerSrc;
        end

        function tf = usesPwrStorage(obj, powerStorage)
            tf = false;
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = obj.emptyVarArr;
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
            [~, powerSrcs] = lv.getPowerSrcsListBoxStr();

            if(not(isempty(powerSrcs)))
%                 addActionTf = lvd_EditActionSetPwrSrcStateGUI(action, lv);

                output = AppDesignerGUIOutput({false});
                lvd_EditActionSetPwrSrcStateGUI_App(action, lv, output);
                addActionTf = output.output{1};
            else
                addActionTf = false;
                warndlg('There are no power sources on the vehicle.  Create one first.','Cannot Create Action','modal');
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot