text:
% KNT_DERHAM: construct the knot vectors for discrete B-spline spaces in the De Rham diagram.
%
%  USAGE:
%
%     [knots, degree] = knt_derham (knots_h1, degree_h1, output_space)
%
% INPUTS:
%
%     knots_h1:     knot vectors for the H^1-conforming space, cell array of size (1 x ndim)
%     degree_h1:    degree for the H^1 space, vector of size (1 x ndim)
%     output_space: string, with one of the following: 'H1', 'Hcurl', 'Hdiv', 'L2'
%
% OUTPUT
%
%    The output are the knot vectors and the degrees for the chosen output space
%    Their size depend on the space:
%
% For H^1 and L^2 (scalar spaces)
%
%     NAME               TYPE          SIZE       DESCRIPTION
%     knots              cell-array  (1 x ndim)
%     degree             vector      (1 x ndim)
%
% For H(curl) and H(div) (vectorial spaces) each component of the
% cell-array contains the same information for one component of the space.
% That is, the knot vector is a cell-array of cell-arrays.
%
%     NAME               TYPE          SIZE
%     knots              cell-array  (1 x ndim)
%     degree             cell-array  (1 x ndim)
%
% Copyright (C) 2010, 2015 Rafael Vazquez
%
%    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/>.

function [knots, degree] = knt_derham (knots_h1, degree_h1, output_space)

  if (nargin < 3)
    output_space = 'Hcurl';
  end

  if (nargout > 2)
    error (['Wrong number of output parameters. The use of knt_derham has changed since version 2.1.' ...
    'Please read the help, or the examples in GeoPDEs, to understand its usage'])
  end

  if (~iscell (knots_h1))
    knots_h1 = {knots_h1};
  end

  ndim = numel (knots_h1);
  if (ndim ~= numel (degree_h1))
    error ('knt_derham: wrong input parameters. The dimension of the knot vector and the degree do not match')
  end

  for ii = 1:ndim
    knots_l2{ii} = knots_h1{ii}(2:end-1);
  end

  if (strcmpi (output_space, 'H1'))
    knots = knots_h1;
    degree = degree_h1;
  elseif (strcmpi (output_space, 'L2'))
    knots = knots_l2;
    degree = degree_h1 - 1;
  elseif (strcmpi (output_space, 'Hcurl'))
    if (ndim < 2 || ndim > 3)
      error ('knt_derham: Wrong dimension to compute the H(curl) space')
    end
    knots = cell (ndim, 1);
    degree = cell (ndim, 1);
    for idim = 1:ndim
      knots{idim} = knots_h1;
      degree{idim} = degree_h1;
      knots{idim}{idim} = knots_l2{idim};
      degree{idim}(idim) = degree_h1(idim) - 1;
    end
  elseif (strcmpi (output_space, 'Hdiv'))
    if (ndim < 2 || ndim > 3)
      error ('knt_derham: Wrong dimension to compute the H(div) space')
    end
    knots = cell (ndim, 1);
    degree = cell (ndim, 1);
    for idim = 1:ndim
      knots{idim} = knots_l2;
      degree{idim} = degree_h1 - 1;
      knots{idim}{idim} = knots_h1{idim};
      degree{idim}(idim) = degree_h1(idim);
    end
  end

%   if (numel (knots_h1) == 2)
%     knots_u1 = {knots_l2{1} knots_h1{2}};
%     knots_u2 = {knots_h1{1} knots_l2{2}};
%
%     if (nargout == 2)
%       varargout = {knots_u1 knots_u2};
%     elseif (nargout == 4 && nargin == 2)
%       degree1 = [degree(1) - 1, degree(2)];
%       degree2 = [degree(1), degree(2) - 1];
%       varargout = {knots_u1 knots_u2 degree1 degree2};
%     elseif (nargout == 3 && nargin == 1)
%       varargout = {knots_u1 knots_u2 knots_l2};
%     elseif (nargout == 6 && nargin == 2)
%       degree1 = [degree(1) - 1, degree(2)];
%       degree2 = [degree(1), degree(2) - 1];
%       degree_l2 = [degree(1) - 1, degree(2) - 1];
%       varargout = {knots_u1 knots_u2 knots_l2 degree1 degree2 degree_l2};
%     else
%       error ('knt_derham: wrong number of input or output parameters');
%     end
%
%   elseif (numel (knots_h1) == 3)
%     knots_u1 = {knots_l2{1} knots_h1{2} knots_h1{3}};
%     knots_u2 = {knots_h1{1} knots_l2{2} knots_h1{3}};
%     knots_u3 = {knots_h1{1} knots_h1{2} knots_l2{3}};
%
%     if (nargout == 3)
%       varargout = {knots_u1 knots_u2 knots_u3};
%     elseif (nargout == 6 && nargin == 2)
%       degree1 = [degree(1) - 1, degree(2), degree(3)];
%       degree2 = [degree(1), degree(2) - 1, degree(3)];
%       degree3 = [degree(1), degree(2), degree(3) - 1];
%       varargout = {knots_u1 knots_u2 knots_u3 degree1 degree2 degree3};
%     else
%       error ('knt_derham: wrong number of input or output parameters');
%     end
%   end
%
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%ISIZE Size of image
%
% N = ISIZE(IM,D) is the size of the D'th dimension of IM.
%
% [W,H] = ISIZE(IM) is the image width W and height H.
%
% WH = ISIZE(IM) is the image size WH = [W H].
%
% [W,H,P] = ISIZE(IM) is the image width W,  height H and and number of
% planes P.  Even if the image has only two dimensions P will be one.
%
% Notes::
% - A simple convenience wrapper on the MATLAB function SIZE.
%
% See also SIZE.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [o1,o2,o3] = isize(im, idx)

    if nargin == 2
        % isize(im, dim)
        o1 = size(im, idx);
    else
        % isize(im)
        s = size(im);
        if nargout == 1
            o1 = [s(2) s(1)];  % width height
        elseif nargout == 2
            o1 = s(2);  % width
            o2 = s(1);  % height, number of rows
        elseif nargout == 3
            o1 = s(2);  % width
            o2 = s(1);  % height, number of rows
            if ndims(im) == 2
                o3 = 1;
            else
                o3 = s(3);
            end
        end
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [x,info,Ai,Bi,BTi,Res,Pro,isFreeDof] = mg(A,b,elem,option,varargin)
%% MG multigrid-type solvers
%
% x = MG(A,b,elem) attempts to solve the system of linear equations A*x =
% b for x using geometric multigrid solvers. Inside mg, an coarsening algorithm
% is applied. See coarsen for the coarsening algorithm on bisection grids.
%
% The method is designed for the system from several finite element
% descritzations of elliptic equations on a grid whose topology
% is given by the array elem.
%
% - 2D: P0, P1, P2, P3, CR, WG
% - 3D: P0, P1, P2, CR, WG
%
% The algorithm is based on fast auxiliary space preconditioner (FASP).
% We first transfer a given element to P1 element and then build V-cycle
% for P1 element. The whole cycle is used as a preconditioner in CG to
% solve the linear system.
%
% Reference: J. Xu, The auxiliary space method and optimal multigrid
% preconditioning techniques for unstructured grids, Computing. 56 (1996)
% 215?235.
%
% The type of the finite element can be build into the size of the problem.
%
% - NT: number of elements         P0 element
% - N : number of vertices         P1 element
% - NE: number of edges            CR element in 2D
% - NF: number of faces            CR element in 3D
% - N + NE: sum of vertices and edges    P2 element
% - N + 2*NE + NT:                       P3 element in 2D
%
% For Dirichlet problems, usually the matrix equation can be restricted to
% free dofs and results a smaller matrix in which the link of the size of
% the matrix and the type of elements is missing. In this case, use
% option.freeDof to provide a logic array for free dofs and the length of
% option.freeDof can be used to determine the type.
%
% So for Dirichlet problems, use
% - mg(A,b,elem,option)      with option.freeDof is given
% - mg(AD,b,elem)       where AD is a larger matrix containing boundary dof
%
% For elements involving dof on edges and faces, additional mesh structure
% such as |edge|, |face| can be provided. Otherwise, the subroutine will
% generate them using |elem|.
%
% For 3-D *adaptive meshes*, HB information is needed, and should be
% listed as the first parameter in varargin.
%
%   - mg(A,b,elem,option,HB)            3-D linear P1 element
%   - mg(A,b,elem,option,HB,edge)       3-D quadratic P2 element
%   - mg(A,b,elem,option,HB,face)       3-D non-conforming CR P1 element
%
% For 3-D meshes obtained by uniformrefine3 or uniformbisect3, no HB is
% needed and simple mg(A,b,elem) is okay.
%
% In some applications, multigrid solvers are used as a build-in block of
% other larger systems. Then use
%
% [~,~,Ai,Bi,BBi,Res,Pro,isFreeDof] = mg(A,f,elem,setupOption);
%
% with setupOption.solver = 'NO' to get hierarichal structure and then use
%
%   e = mg(A,r,elem,option,Ai,Bi,BBi,Res,Pro,isFreeDof);
%
% This will save the setting time. See diapreStokes, tripremixPoisson.
%
% When testing the solvers, more options can be provided.
%
%   x = MG(A,b,elem,options) specifies options in the following list.
%   - option.x0: the initial guess. Default setting x0 = 0.
%   - option.tol: the tolerance of the convergence. Default setting 1e-8.
%   - option.maxIt: the maximum number of iterations. Default setting 200.
%   - option.N0: the size of the coarest grid. Default setting 500.
%   - option.mu: smoothing steps. Default setting 1
%   - option.coarsegridsolver: solver used in the coarest grid. Default
%     setting: direct solver.
%   - option.freeDof: free d.o.f
%   - option.solver: various cycles and Krylov space methods
%       * 'NO'     only setup the transfer matrix
%       * 'Vcycle'      V-cycle MultiGrid Method
%       * 'Wcycle'      W-cycle MultiGrid Method
%       * 'Fcycle'      Full cycle Multigrid Method
%       * 'cg'     MG-Preconditioned Conjugate Gradient
%       * 'minres' MG-Preconditioned Minimal Residual Method
%       * 'gmres'  MG-Preconditioned Generalized Minimal Residual Method
%       * 'bicg'   MG-Preconditioned BiConjugate Gradient Method
%       * 'bicgstable' MG-Preconditioned BiConjugate Gradient Stabilized Method
%       * 'bicgstable1' MG-Preconditioned BiConjugate Gradient Stabilized Method
%       The default setting is 'cg' which works well for SPD matrices. For
%       non-symmetric matrices, try 'gmres' and for symmetric but indefinite
%       matrices, try 'minres' or 'bicg' sequences.
%       The string option.solver is not case sensitive.
%   - option.preconditioner:  multilevel preconditioners including:
%       * 'V'   V-cycle MultiGrid used as a Preconditioner
%       * 'W'   W-cycle MultiGrid used as a Preconditioner
%       * 'F'   Full cycle Multigrid used as a Preconditioner
%       * 'bpx' BPX-Preconditioner
%   - option.printlevel: the level of screen print out
%       * 0: no output
%       * 1: name of solver and convergence information (step, err, time)
%       * 2: convergence history (err in each iteration step)
%
%   [x,info] = MG(A,b,elem) also returns information of the solver
%   - info.flag:
%       * 0: mg converged to the desired tolerance tol within maxIt iterations
%       * 1: mg iterated maxIt times but did not converge.
%       * 2: direct solver
%   - info.itStep: the iteration number at which x was computed.
%   - info.time: the cpu time to get x
%   - info.err: the approximate relative error in the energy norm in
%   err(:,1) and the relative residual norm(b-A*x)/norm(b) in err(:,2). If
%   flag is 0, then max(err(end,:)) <= tol.
%   - info.stopErr: the error when iteration stops
%
%   Example:
%
%   [node,elem] = squaremesh([0 1 0 1],0.5);
%   for k = 1:8
%     [node,elem] = uniformrefine(node,elem);
%   end
%   pde.f = inline('p(:,1).*p(:,2)','p');
%   pde.g_D = inline('zeros(size(p,1),1)','p');
%   option.solver = 'notsolve';
%   [soln,eqn] = Poisson(node,elem,[],pde,option);
%   fprintf('\n Number of unknowns: %8.0u\n',length(eqn.b))
%   tic; display('Direct solver'); u = eqn.A\eqn.b; toc;
%   tic; x = mg(eqn.A,eqn.b,elem); toc;
%   format shorte
%   fprintf('Difference between direct and mg solvers %0.2g \n',norm(u-x));
%
% See also mgMaxwell, mgstokes
%
% Documentation in Help browser ifem mgdoc
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

t = cputime;
%% Size of systems
Nborig = size(b,1);                % number of dof
Ndof = Nborig;
nb = size(b,2);                    % number of bs
N = max(elem(:));                  % number of nodes
NT = size(elem,1);                 % number of elements
dim = size(elem,2)-1;
if N > NT       % 2-D quad mesh
    dim = 2;
end

%% Options
% Assign default values to unspecified parameters
if ~exist('option','var')
    option = [];
end
option = mgoptions(option,Nborig);    % parameters
x0 = option.x0;
if size(x0,2) == 1
    x0 = repmat(x0,1,nb);
end
N0 = option.N0;
tol = option.tol;
maxIt = option.solvermaxit;
mu = option.smoothingstep;
smoothingRatio = option.smoothingratio; % for variable smoothing
solver = option.solver;
preconditioner = option.preconditioner;
coarsegridsolver = option.coarsegridsolver;
printlevel = option.printlevel;
setupflag = option.setupflag;
if nargin > 8      % with the multilevel structure in the input
    setupflag = 0;
end

%% Set up multilevel structure
if setupflag == true
%% eliminate isolated dof
if isfield(option,'freeDof') % freeDof is given
   if islogical(option.freeDof)
      Ndof = max(length(option.freeDof),Nborig);
   else
      error('Provide logic array for option.freeDof');
   end
   isFreeDof = false(Ndof,1);
   isFreeDof(option.freeDof) = true;
   isFixDof = true(Ndof,1);
   isFixDof(isFreeDof) = false;
else % find free dofs and eliminate isolated dofs
    Ndof = Nborig;
    degreeA = sum(spones(A));  % degree
    isFreeDof = false(size(A,1),1);
    isFreeDof(degreeA>1) = true;
    isFixDof = ~isFreeDof;
    isFixDof(degreeA == 0) = false;
end
if Nborig > sum(isFreeDof) % truncate the matrix
    xD = zeros(Nborig,nb);
    Nfix = sum(isFixDof);
    ADinv = spdiags(1./diag(A(isFixDof,isFixDof)),0,Nfix,Nfix);
    xD(isFixDof,:) = ADinv*b(isFixDof,:);
    A = A(isFreeDof,isFreeDof);
end
if Nborig > sum(isFreeDof) % truncate the rhs
    b = b(isFreeDof,:);
end
if size(x0,1) > sum(isFreeDof) % truncate the initial guess
    x0 = x0(isFreeDof,:);
    option.x0 = x0;
end
isFreeNode = isFreeDof(1:N); % the first N element is for P1

%% Additional mesh structure for different elements
% list of elements: P0,P2,P3,CR,HB,WG
if Ndof > N % other than P1 element
    NE = N + NT -1; % estimate of NE by Euler formula
    if dim == 3
        NF = 2*NT;  % estimate of NF
    end
    % get additional data structure
    if Ndof ~= NT   % not piecewise constant element
        if dim == 2
            if nargin > 4
                edge = varargin{1};
                NE = size(edge,1);
            else
                [tempvar,edge] = dofedge(elem);
            end
        elseif dim == 3
            if nargin > 5 % additional data structure is in the input
                if size(varargin{1},2) == 2 % edge
                    edge = varargin{1};
                    NE = size(edge,1);
                elseif size(varargin{1},2) == 3 % face
                    face = varargin{1};
                    NF = size(face,1);
                end
            else % if no additional edge/face is input, generate it
                 [tempvar,edge] = dof3edge(elem); %#ok<*ASGLU>
                 [tempvar,face] = dof3face(elem);
                 NE = size(edge,1);
                 NF = size(face,1);
            end
        end
    end
    % transfer operators from P1 to the current element
    isFreeNode = true(N,1);
    if Ndof == NT % piecewise constant element
       if dim == 2
          P1toP0 = sparse([1:NT;1:NT;1:NT]',elem,ones(3*NT,1)/3,NT,N);
       elseif dim == 3
          P1toP0 = sparse([1:NT;1:NT;1:NT;1:NT]',elem,ones(4*NT,1)/4,NT,N);
       end
       auxPro = P1toP0(isFreeDof,1:N);
       isFreeNode = [];
    end
    if Ndof == N + NE % quadratic element
        isFreeNode = isFreeDof(1:N);
        P1toP2 = sparse([(1:N)'; N+(1:NE)'; N+(1:NE)'], ...
                        [(1:N)'; double(edge(:))],...
                        [ones(N,1); 0.5*ones(2*NE,1)]',Ndof,N);
        auxPro = P1toP2(isFreeDof,isFreeNode);
    end
    if (dim == 2) && (Ndof == N + 2*NE +NT) % cubic element in 2-D
        isFreeNode = isFreeDof(1:N);
        P1toP3 = sparse([(1:N)'; N+2*(1:NE)'-1; N+2*(1:NE)'-1; N+2*(1:NE)';...
                          N+2*(1:NE)';N+2*NE+(1:NT)';N+2*NE+(1:NT)';N+2*NE+(1:NT)'], ...
                        [(1:N)'; double(edge(:));double(edge(:));elem(:);],...
                        [ones(N,1); 2/3*ones(NE,1);1/3*ones(NE,1);...
                         1/3*ones(NE,1);2/3*ones(NE,1);1/3*ones(3*NT,1)]',Ndof,N);
         auxPro = P1toP3(isFreeDof,isFreeNode);
    end
    if (dim == 2) && (Ndof == NE)  % 2-D CR nonconforming element
        isFreeNode(edge(isFixDof,:)) = false;
        P1toCR = sparse([1:NE;1:NE]',double(edge(:)),0.5*ones(2*NE,1),NE,N);
        auxPro = P1toCR(isFreeDof,isFreeNode);
    end
    if (dim == 3) && (Ndof == NF) % 3-D CR nonconforming element
        isFreeNode(face(isFixDof,:)) = false;
        P1toCR = sparse([1:NF;1:NF;1:NF]',double(face(:)),ones(3*NF,1)/3,NF,N);
        auxPro = P1toCR(isFreeDof,isFreeNode);
    end
    if (dim == 2) && (Ndof == NT + NE) % 2-D weak Galerkin element (P0,P0,RT0) element
        fixEdgeDof = find(isFixDof) - NT;
        isFreeNode(edge(fixEdgeDof,:)) = false;
        P1toWG = sparse([repmat((1:NT)',3,1); repmat((NT+1:Ndof)',2,1)], ... % i
                        [elem(:);             double(edge(:))], ...   % j
                        [ones(3*NT,1)/3;      ones(2*NE,1)/2], NT+NE, N);
        auxPro = P1toWG(isFreeDof,isFreeNode);
    end
    if (dim == 3) && (Ndof == NT + NF) % 3-D weak Galerkin element (P0,P0,RT0) element
        fixFaceDof = find(isFixDof) - NT;
        isFreeNode(face(fixFaceDof,:)) = false;
        P1toWG = sparse([repmat((1:NT)',4,1); repmat((NT+1:Ndof)',3,1)], ... % i
                        [elem(:);           double(face(:))], ...   % j
                        [ones(4*NT,1)/4;    ones(3*NF,1)/3], NT+NF, N);
        auxPro = P1toWG(isFreeDof,isFreeNode);
    end
    if (nargin > 4) && (dim == 2) && ischar(varargin{1}) && strcmp(varargin{1},'HB') % HB basis element
        isFreeNode = isFreeDof(1:N);
        P1toPk = speye(Ndof,N);
        auxPro = P1toPk(isFreeDof,isFreeNode);
    end
end

%% Hierarchical Structure of Mesh
if dim == 2  % 2-D
    [HB, NL, level] = HBstructure(elem,N0);
end
if dim == 3  % 3-D
    if nargin > 4
        HBmesh = varargin{end}; % need HB for bisection refinement
        [HB, NL, level] = HBstructure3(elem,N0,HBmesh);
    else % no HBmesh is given. only work for the red uniform refinement
        [HB, NL, level] = HBstructure3(elem,N0);
    end
end

%% Transfer operators between multilevel meshes for P1 element
% standard prolongation and restriction operator for P1 element
[Pro,Res] = transferoperator(HB,NL,isFreeNode);
% add one more level from P1 to the current element
if Ndof > N
    if ~exist('auxPro','var')
        disp('The current element is not supported by mg');
    else
        Pro{level} = auxPro;
        Res{level+1} = Pro{level}';
        level = level + 1;
    end
end
clear HB auxPro

%% Matrices in each level
Ai = cell(level,1);
Ai{level} = A;
if level == 1
    Bi{1} = tril(Ai{1});
    BTi{1} = transpose(Bi{1});
    Res = []; Pro = [];
end
for j = level:-1:2
    Ai{j-1} = Res{j}*Ai{j}*Pro{j-1};           % Ac = Res*Af*Pro
    switch option.smoother
        case 'GS'
            Bi{j} = tril(Ai{j});        % Forward Gauss-Seidel   B = D+L
            BTi{j} = triu(Ai{j});       % Backward Gauss-Seidel BT = D+U
        case 'JAC'                      % Jacobi iteration  B = BT = D;
            Bi{j} = spdiags(diag(Ai{j}),0,size(Ai{j},1),size(Ai{j},1));
            BTi{j} = Bi{j};
    end
    if option.smoothingparameter~=1     % scaling of the smoother
        Bi{j} = Bi{j}/option.smoothingparameter;
        BTi{j} = BTi{j}/option.smoothingparameter;
    end
end
end % end for setup

%% Only set up the transfer operators
if strcmp(solver,'NO')
    x = x0; flag = 0; itStep = 0; err = 0; time = cputime-t;
    info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));
    return
end

%% No need of set up
if setupflag == false
    Ai = varargin{1};
    Bi = varargin{2};
   BTi = varargin{3};
   Res = varargin{4};
   Pro = varargin{5};
   if length(varargin)>=6
      isFreeDof = varargin{6};
   else
      isFreeDof = true(Ndof,1);
   end
   level = length(Ai);
   isFixDof = ~isFreeDof;
   xD = zeros(Ndof,nb);
   if any(isFixDof) && Ndof == length(isFreeDof) % larger system
        Nfix = sum(isFixDof);
        ADinv = spdiags(1./diag(A(isFixDof,isFixDof)),0,Nfix,Nfix);
        xD(isFixDof,:) = ADinv*b(isFixDof,:);
        b = b(isFreeDof,:);
        A = A(isFreeDof,isFreeDof);
        x0 = x0(isFreeDof,:);
   end
end
if condest(Ai{1}) > 1e16 % Ai{1} is singular
    Ai{1} = Ai{1} + 1e-12*speye(size(Ai{1}));
%     coarsegridsolver = 'pcg';
end

%% No coarsening or coarsened nodes is small
if Ndof <= N    % linear element
    if level == 1
        if Ndof < N0 % small size
            solver = 'DIRECT';
        else
            [x,info] = amg(A,b,option);
            if any(isFixDof) && Nborig > size(x,1) % a larger system
                xD(isFreeDof,:) = x;
                x = xD;
            end
            Ai{1} = A; Bi{1} = tril(A); BTi{1} = Bi'; Res = []; Pro = [];
            fprintf('No hierarchical structure of the mesh is found. AMG is used. \n');
            return
        end
    end
end

%% Direct solver
if strcmp(solver,'DIRECT')
    x = A\b;                       % use direct solver
    flag = 2; itStep = 0; err = norm(b-A*x)/norm(b); time = cputime - t;
    %% Modify x to include fix dof
    if any(isFixDof) && Nborig > size(x,1)
        xD(isFreeDof,:) = x;
        x = xD;
    end
    if printlevel >= 1
        fprintf('Direct solver \n')
        fprintf('#dof: %8.0u,  #nnz: %8.0u,  iter: %2.0u,  err = %8.4e,  time = %4.2g s\n',...
                 Ndof, nnz(A), itStep, max(err(end,:)), time)
        if level == 1
            fprintf('No coarsening of the grid. \n')
        end
    end
    info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));
    Ai{1} = A; Bi{1} = tril(A); BTi{1} = Bi'; Res = []; Pro = [];
    return
end

%% Krylov iterative methods use Multigrid-type Preconditioners
% set up preconditioner
if ~strcmp(solver(2:end),'CYCLE')
    switch preconditioner
        case 'V'
            prefunc = @vcycle;
            if printlevel >= 1
                fprintf('\n Multigrid V-cycle Preconditioner with ')
            end
        case 'W'
            prefunc = @wcycle;
            if printlevel >= 1
                fprintf('\n Multigrid W-cycle Preconditioner with ')
            end
        case 'F'
            prefunc = @fcycle;
            if printlevel >= 1
                fprintf('\n Multigrid Full Cycle Preconditioner with ')
            end
        case 'BPX'
            % modify smoother
            for j = level:-1:2
                Di{j} = diag(Ai{j});
            end
            prefunc = @bpx;
            if printlevel >= 1
                fprintf('BPX Preconditioner with ')
            end
    end
end
% initial set up
k = 1;
x = x0;
r = b - A*x;
nb = max(sqrt(sum(b.^2,1)));
err = zeros(maxIt,2);
if nb > eps  % nb is non-zero
    err(1,:) = max(sqrt(sum(r.^2,1)))/nb;
else
    err(1,:) = max(sqrt(sum(r.^2,1)));
end
% solvers
switch solver
    case 'CG'
        if printlevel >= 1
            fprintf('Conjugate Gradient Method\n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            % compute Br by MG
            Br = prefunc(r);
            % update tau, beta, and p
            rho = dot(Br,r);  % e'*ABA*e approximates e'*A*e
            if k == 1
                p = Br;
            else
                beta = rho./rho_old;
                p = Br + beta.*p;
            end
            % update alpha, x, and r
            Ap = A*p;
            alpha = rho./dot(Ap,p);
            r = r - alpha.*Ap;
            x = x + alpha.*p;
            rho_old = rho;
            k = k + 1;
            % compute err for the stopping criterion
        %     err(k,1) = alpha*sqrt(p'*Ap/(x'*A*x)); % increamental error in energy norm
            err(k,1) = max(sqrt(abs(rho./dot(x,b)))); % approximate relative error in energy norm
            if nb > eps
                err(k,2) = max(sqrt(sum(r.^2,1)))/nb; % relative error of the residual in L2-norm
            else
                err(k,2) = max(sqrt(sum(r.^2,1)));
            end
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MGCG iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'VCYCLE'
        if printlevel >= 1
            fprintf('\n Multigrid Vcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = vcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            % Step 1: Form residual r
            r = r - A*Br;
            err(k,1) = max(sqrt(abs(dot(Br,r)/dot(x,b)))); % approximate relative error in energy norm
            if nb > eps
                err(k,2) = max(sqrt(sum(r.^2,1)))/nb; % relative error of the residual in L2-norm
            else
                err(k,2) = max(sqrt(sum(r.^2,1)));
            end
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Vcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'WCYCLE'
        if printlevel >= 1
            fprintf('\n Multigrid Wcycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Vcylce MG
            Br = wcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            err(k,1) = max(sqrt(abs(dot(Br,r)/dot(x,b)))); % approximate relative error in energy norm
            % Step 1: Form residual r
            r = r - A*Br;
            if nb > eps
                err(k,2) = max(sqrt(sum(r.^2,1)))/nb; % relative error of the residual in L2-norm
            else
                err(k,2) = max(sqrt(sum(r.^2,1)));
            end
            if printlevel >= 2
                fprintf('#dof: %8.0u,  #nnz: %8.0u, MG Wcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'FCYCLE'
        if printlevel >= 1
            fprintf('\n Multigrid Full Cycle Iteration \n')
        end
        while (max(err(k,:)) > tol) && (k <= maxIt)
            k = k + 1;
            % Step 2: Compute Br by one Full cycle MG
            Br = fcycle(r);
            % Step 3: Correct the solution
            x = x + Br;
            err(k,1) = max(sqrt(abs(dot(Br,r)/dot(x,b)))); % approximate relative error in energy norm
            % Step 1: Form residual r
            r = r - A*Br;
            if nb > eps
                err(k,2) = max(sqrt(sum(r.^2,1)))/nb; % relative error of the residual in L2-norm
            else
                err(k,2) = max(sqrt(sum(r.^2,1)));
            end
            if printlevel >= 2
                fprintf('#dof: %8.0u, #nnz: %8.0u, MG Fcycle iter: %2.0u, err = %8.4e\n',...
                         Ndof, nnz(A), k-1, max(err(k,:)));
            end
        end
        err = err(1:k,:);
        itStep = k-1;
    case 'MINRES'
        if printlevel >= 1
            fprintf('Minimum Residual Method \n')
        end
        [x,flag,err,itStep] = minres(A,b,tol,maxIt,prefunc,[],x0);
    case 'GMRES'
        if printlevel >= 1
            fprintf('General Minimum Residual Method\n')
        end
        if isfield(option,'restart')
            restart = option.restart;
        else
            restart = min(N,10);
        end
        [x,flag,err,itStep] = gmres(A,b,restart,tol,maxIt,prefunc,[],x0);
        itStep = (itStep(1)-1)*restart + itStep(2);
    case 'BICG'
        if printlevel >= 1
            fprintf('BiConjugate Gradient Method\n')
        end
        [x,flag,err,itStep] = bicg(A,b,tol,maxIt,prefunc,[],x0);
    case 'BICGSTAB'
        if printlevel >= 1
            fprintf('Stablilized BiConjugate Gradient Method\n')
        end
        [x,flag,err,itStep] = bicgstab(A,b,tol,maxIt,prefunc,[],x0);
    case 'BICGSTAB1'
        if printlevel >= 1
            fprintf('Stablilized BiConjugate Gradient Method\n')
        end
        [x,flag,err,itStep] = bicgstab1(A,b,tol,maxIt,prefunc,[],x0);
end

%% Modify x to include fix dof
if Nborig > size(x,1) % a larger system
    xD(isFreeDof,:) = x;
    x = xD;
end

%% Output
if k > maxIt
    flag = 1;
else
    flag = 0;
end
time = cputime - t;
if printlevel >= 2
    fprintf('#dof: %8.0u, level: %2.0u,   coarse grid %2.0u, #nnz: %8.0u\n',...
              size(b,1), level, size(Ai{1},1), nnz(Ai{1}))
end
if printlevel >= 1
    fprintf('#dof: %8.0u,  #nnz: %8.0u, smoothing: (%1.0u,%1.0u), iter: %2.0u,   err = %8.2e,   time = %4.2g s\n',...
                 size(b,1), nnz(A), mu, mu, itStep, max(err(end,:)), time)
end
if (flag == 1) && (printlevel>0)
   fprintf('NOTE: the iterative method does not converge! \n');
end
info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions vcycle, wcycle, fcycle, bpx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Vcycle MG
    function Br = vcycle(r,J)  % solve equations Ae = r in each level
    if nargin<=1
        J = level;
    end
    ri = cell(J,1);            % residual in each level
    ei = cell(J,1);            % correction in each level
    mui = cell(J,1);           % variable smoothing steps
    ri{J} = r;
    mui{J} = mu;
    for i = J:-1:2
        ei{i} = Bi{i}\ri{i};   % pre-smoothing
        for s = 1:mui{i}-1     % extra mu-1 steps smoothing
            if mod(s,2)        % switch between B and BT to symmetrize smoother
                ei{i} = ei{i} + BTi{i}\(ri{i}-Ai{i}*ei{i});
            else
                ei{i} = ei{i} + Bi{i}\(ri{i}-Ai{i}*ei{i});
            end
        end
        ri{i-1} = Res{i}*(ri{i} - Ai{i}*ei{i});
        mui{i-1} = ceil(mui{i}*smoothingRatio);
    end
    switch coarsegridsolver
        case 'direct'
            ei{1} = Ai{1}\ri{1};   % direct solver in the coarest level
        case 'pcg'
            D = spdiags(diag(Ai{1}),0,size(Ai{1},1),size(Ai{1},1));
            [ei{1},flag] = pcg(Ai{1},ri{1},1/size(Ai{1},1),1000,D);
        case 'amg'
            amgoption.printlevel = 0;
            ei{1} = amg(Ai{1},ri{1},amgoption);
    end
    for i = 2:J
        ei{i} = ei{i} + Pro{i-1}*ei{i-1};
        for s = 1:mui{i}       % post-smoothing
            if mod(s,2)        % switch between B and BT to symmetrize smoother
                ei{i} = ei{i} + BTi{i}\(ri{i}-Ai{i}*ei{i});
            else
                ei{i} = ei{i} + Bi{i}\(ri{i}-Ai{i}*ei{i});
            end
        end
    end
    Br = ei{J};
    end

%% Wcycle MG
    function e = wcycle(r,J)    % solve equations Ae = r in level J
    if nargin<=1
        J = level;
    end
    if J == 1
        switch coarsegridsolver
            case 'direct'
                e = Ai{J}\r;   % direct solver in the coarest level
            case 'pcg'
                D = spdiags(diag(Ai{J}),0,size(Ai{J},1),size(Ai{J},1));
                [e,flag] = pcg(Ai{J},r,1/size(Ai{J},1),1000,D);
            case 'amg'
                amgoption.printlevel = 0;
                e = amg(Ai{J},r,amgoption.printlevel);
        end
        return
    end
    % fine grid pre-smoothing
    e = Bi{J}\r;        % pre-smoothing
    for s = 1:mu-1      % extra mu-1 steps smoothing
        if mod(s,2)
            e = e + Bi{J}\(r-Ai{J}*e);
        else
            e = e + BTi{J}\(r-Ai{J}*e);
        end
    end
    % restriction
    rc = Res{J}*(r - Ai{J}*e);
    % coarse grid correction twice
    ec = wcycle(rc,J-1);
    ec = ec + wcycle(rc - Ai{J-1}*ec,J-1);
    % prolongation
    e = e + Pro{J-1}*ec;
    % fine grid post-smoothing
    for s = 1:mu
        if mod(s,2)
            e = e + Bi{J}\(r-Ai{J}*e);
        else
            e = e + BTi{J}\(r-Ai{J}*e);
        end
    end
    end

%% Fcycle MG
    function Br = fcycle(r)
    ri = cell(level,1);            % residual in each level
    ei = cell(level,1);            % correction in each level
    ri{level} = r;
    for i = level:-1:2
        ei{i} = vcycle(ri{i},i);   % pre-smoothing
        for s = 1:mu-1             % extra smoothing steps
            ei{i} = ei{i} + vcycle(ri{i}-Ai{i}*ei{i},i); % pre-smoothing
        end
        ri{i-1} = Res{i}*(ri{i} - Ai{i}*ei{i});
    end
    switch coarsegridsolver
        case 'direct'
            ei{1} = Ai{1}\ri{1};   % direct solver in the coarest level
        case 'pcg'
            D = spdiags(diag(Ai{1}),0,size(Ai{1},1),size(Ai{1},1));
            [ei{1},flag] = pcg(Ai{1},ri{1},1/size(Ai{1},1),1000,D);
        case 'amg'
            amgoption.printlevel = 0;
            ei{1} = amg(Ai{1},ri{1},amgoption);
    end
    for i = 2:level
        ei{i} = ei{i} + Pro{i-1}*ei{i-1};
        for s = 1:mu   % post-smoothing
            ei{i} = ei{i} + vcycle(ri{i}-Ai{i}*ei{i},i); % post-smoothing
        end
    end
    Br = ei{level};
    end

%% BPX preconditioner
    function Br = bpx(r)
    ri = cell(level,1);            % record residual in each level
    ei = cell(level,1);            % record err in each level
    % compute Br by BPX
    ri{level} = r;
    for i = level:-1:2
        ei{i} = ri{i}./Di{i};       % Jacobi smoothing
        ri{i-1} = Res{i}*ri{i};     % restriction of the residual
    end
    switch coarsegridsolver
        case 'direct'
            ei{1} = Ai{1}\ri{1};   % direct solver in the coarest level
        case 'pcg'
            D = spdiags(diag(Ai{1}),0,size(Ai{1},1),size(Ai{1},1));
            [ei{1},flag] = pcg(Ai{1},ri{1},1/size(Ai{1},1),1000,D);
        case 'amg'
            amgoption.printlevel = 0;
            ei{1} = amg(Ai{1},ri{1},amgoption);
    end
    for i=2:level
        ei{i} = ei{i} + Pro{i-1}*ei{i-1};  % prolongation of the correction
    end
    Br = ei{level};
    end
end %% end of mg

repo: ifem


text:
function F = bst_window( Method, L, R )
% BST_WINDOW: Generate a window of length L, ot the given type: hann, hamming, blackman, parzen, tukey
%
% USAGE:  F = bst_window( Method, L, R )
%
% References:
%    Formulas documented on Wikipedia: http://en.wikipedia.org/wiki/Window_function

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2013-2014

% Parse inputs
if (nargin < 3) || isempty(R)
    R = [];
end
if (nargin < 2) || isempty(L)
    L = 64;
end
if (nargin < 1) || isempty(Method)
    Method = 'hann';
end
% Calculate normalized time vector
t = (0:L-1)' ./ (L-1);
% Switch according to windowing method
switch (lower(Method))
    case 'hann'
        F = 0.5 - 0.5 * cos(2*pi*t);
    case 'hamming'
        F = 0.54 - 0.46 * cos(2*pi*t);
    case 'blackman'
        F = 0.42 - 0.5 * cos(2*pi*t) + 0.08 * cos(4*pi*t);
    case 'tukey'
        % Tukey window = tapered cosine window, cosine lobe of width aL/2
        if isempty(R)
            R = 0.5;
        end
        a = R;
        % If a=0: square function
        if (a <= 0)
            F = ones(L,1);
        % If a=1: Hann window
        elseif (a >= 1)
            F = bst_window('hann', L);
        % Else: Function in three blocks
        else
            % Define three blocks
            len = floor(a*(L-1)/2) + 1;
            t1 = t(1:len);
            t3 = t(L-len+1:end);
            % Window is defined in three sections: taper, constant, taper
            F = [ 0.5 * (1 + cos(pi * (2*t1/a - 1)));  ...
                  ones(L-2*len,1); ...
                  0.5 * (1 + cos(pi * (2*t3/a - 2/a + 1)))];
        end
    case 'parzen'
        % Reference:
        % Harris FJ, On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform,
        % Proceedings of IEEE, Vol. 66, No. 1, January 1978   [Equation 37]

        % Time indices
        t = -(L-1)/2 : (L-1)/2;
        i1 = find(abs(t) <= (L-1)/4);  %   0 <= |n| <= N/4
        i2 = find(abs(t) >  (L-1)/4);  % N/4 <  |n| <  N/2

        % Definition of the Parzen window: 2 parts
        F = zeros(length(t), 1);
        F(i1) = 1 - 6 .* (t(i1)/L*2).^2 .* (1 - abs(t(i1))/L*2);
        F(i2) = 2 .* (1 - abs(t(i2))/L*2) .^3;

    otherwise
        error(['Unsupported windowing method: "' Method '".']);
end

author: brainstorm-tools


text:
function test_failed=test_wmdct
% Test the algorithm using LONG windows.

which comp_dwiltiii
which comp_idwiltiii

disp(' ===============  TEST_WMDCT ================');

Lr=[4, 6, 8,12,16,12,18,32,30];
Mr=[2, 3, 2, 3, 4, 2, 3, 4, 3];

test_failed=0;

for ii=1:length(Lr);
  for W=1:3
    for ftype=1:2
      for wtype=1:2
	L=Lr(ii);
	M=Mr(ii);

	a=M;

	if wtype==1
	  % Full length window
	  g=pgauss(L);
	  gd=wildual(g,M);
          wtype='LONG';
	else
	  g=firwin('sqrthann',2*M,'2');
	  gd=g;
          wtype='FIR ';
	end;

	if ftype==1
	  % Complex-valued test case
	  f=tester_crand(L,W);
	  S='CMPLX';
	else
	  % Real-valued tes
	  f=tester_rand(L,W);
	  S='REAL ';
	end;

	c=wmdct(f,g,M,L);

	a=M;

	c2=ref_dwiltiii(f,g,a,M);
	r=iwmdct(c,gd);

	res=norm(c(:)-c2(:));

        [test_failed,fail]=ltfatdiditfail(res,test_failed);
	s=sprintf('REF  %s %s L:%3i W:%2i a:%3i M:%3i %0.5g %s',S,wtype,L,W,a,M,res,fail);
	disp(s);

	rdiff=f-r;
	res=norm(rdiff(:));

        [test_failed,fail]=ltfatdiditfail(res,test_failed);

	s=sprintf('REC  %s %s L:%3i W:%2i a:%3i M:%3i %0.5g %s',S,wtype,L,W,a,M,res,fail);
	disp(s);

        g=wilorth(M,L);
        c=wmdct(f,g,M);
        r=iwmdct(c,g);
        rdiff=f-r;

	res=norm(rdiff(:));

        [test_failed,fail]=ltfatdiditfail(res,test_failed);
	s=sprintf('ORTH %s %s L:%3i W:%2i a:%3i M:%3i %0.5g %s',S,wtype,L,W,a,M,res,fail);
	disp(s);

      end;
    end;
  end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% eeg_interp_scalp_script - Script to interpolate scalp potentials over mesh
%
% This is part of the development and testing of
% functions to implement:
%
% Oostendorp T, Oosterom A, & Huiskamp G (1989),
% Interpolation on a triangulated 3D surface.
% Journal of Computational Physics, 80: 331-343.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:51 $

% Licence:  GNU GPL, no implied or express warranties
% History:  04/2002, Darren.Weber_at_radiology.ucsf.edu
%           - needs verification and conversion to function
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

% run R. Oostendorp's script to load simulated data
cd d:\matlab\cvs\eeg_toolbox\lapint
%sphere_load;
%cd ..

% If this script has already been run and the data
% saved, then just load the saved data and return
if isequal(exist('eeg_interp_scalp_mesh_data.mat'),2),
    fprintf('\nLoading saved data.\n');
    load 'eeg_interp_scalp_mesh_data';
    return
end

p = eeg_toolbox_defaults('create');
% Load electrode coordinates 'Cartesian'
p = elec_open(p);
% Load potential values at each electrode
p = eeg_open(p);
% Load a tesselation dataset (scalp, 'skull', cortex)
p = mesh_open(p);

% Find the nearest scalp vertices to each
% electrode vertex
p = mesh_plot(p);
close all

% Extract tesselations from p struct
scalpvert = p.mesh.data.vertices{3};
scalpface = p.mesh.data.faces{3};
elecvert  = p.mesh.data.vertices{4};
elecface  = p.mesh.data.faces{4};

% Extract electric potential at time t for electrodes
Velec = p.volt.data;

clear p;

% This command retains the shape of scalp, while reducing
% the number of faces from ~4000 to 1000.  In doing so, it
% reduces the resolution of the scalp mesh and we need to
% recompute the nearest scalp vertices for each electrode
% vertex below
%[scalpface, scalpvert] = reducepatch(scalpface,scalpvert,1000);

% find the scalp vertex indices for each electrode vertex.
% In this example, this generates replicate indices when
% the reducepatch command above is executed
Scalpindex = dsearchn(scalpvert,elecvert)';

% Calculate the Laplacian matrix
lap = mesh_laplacian(scalpvert,scalpface);

% Calculate interpolation matrix, based on Laplacian
[Lint, keepindex, repindex] = mesh_laplacian_interp(lap,Scalpindex);

% Interpolate potential to scalp
time = 100;
if isempty(repindex),
    Vscalp = Lint * Velec(time,:)';
else
    Vscalp = Lint * Velec(time,keepindex)';
end

%save eeg_interp_scalp_mesh_data

return

% Plot scalp with voltages at all interpolated vertices

fig = figure;
Hp = patch('vertices',scalpvert,'faces',scalpface,'FaceVertexCdata',Vscalp,'facecolor','interp','edgecolor','none');
axis off tight vis3d
lighting gouraud

map = eeg_colormap('Red/Blue/White');
colormap(map)

[az,el] = view;
lit = lightangle(az,el);

set(Hp,'AmbientStrength',.7);
set(Hp,'DiffuseStrength',.3);
set(Hp,'SpecularStrength',0);
set(Hp,'SpecularColorReflectance',0);

%colorbar

% Animation of the timeseries & save movie
set(gcf,'BackingStore','off');
set(Hp,'EraseMode','normal');
figure(fig);
% Define movie region as figure size
rect = get(fig,'Position');
rect(1:2) = [0 0];

clear M

start = 120;
finish = 160;
for t=start:finish,
    Vscalp = Lint * Velec(t,:)';
    set(Hp,'FaceVertexCdata',Vscalp);
    drawnow;
    M(t-(start-1)) = getframe(gcf,rect);
end
% Play the movie in another figure
playtimes = 1;
fps = 15;
figure; axis off; movie(gcf,M,playtimes,fps,rect);
movie2avi(M,'scalp_interpolation.avi','quality',100);

return

author: PatternRecognition

repo: OpenBMI


text:
function value = hex_digit_to_i4 ( c )

%*****************************************************************************80
%
%% HEX_DIGIT_TO_I4 converts a hexadecimal digit to an I4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the hexadecimal digit, '0'
%    through '9', or 'A' through 'F', or also 'a' through 'f'
%    are allowed.
%
%    Output, integer VALUE, the corresponding integer, or -1 if C was illegal.
%
  if ( '0' <= c && c <= '9' )
    value = c - '0';
  elseif ( 'a' <= c && c <= 'f' )
    value = 10 + c - 'a';
  elseif ( 'A' <= c && c <= 'F' )
    value = 10 + c - 'A';
  elseif ( c == ' ' )
    value = 0;
  else
    value = -1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = non_domination_sort_mod(x, M, V)

%% function f = non_domination_sort_mod(x, M, V)
% This function sort the current popultion based on non-domination. All the
% individuals in the first front are given a rank of 1, the second front
% individuals are assigned rank 2 and so on. After assigning the rank the
% crowding in each front is calculated.

%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.

% [N, m] = size(x);
% clear m;
N = size(x,1); % modified by zzb
x_temp = zeros(N,M + V + 1); % modified by zzb
x_temp(:,1:M + V) = x(:,1:M + V); % modified by zzb
% Initialize the front number to 1.
front = 1;
% There is nothing to this assignment, used only to manipulate easily in
% MATLAB.
% F(front).f = [];
% individual = [];
F = struct('f',cell(1,N)); % modified by zzb
individual = struct('n',cell(1,N),'p',cell(1,N)); % modified by zzb

%% Non-Dominated sort.
% The initialized population is sorted based on non-domination. The fast
% sort algorithm [1] is described as below for each

% ?for each individual p in main population P do the following
%   ?Initialize Sp = []. This set would contain all the individuals that is
%     being dominated by p.
%   ?Initialize np = 0. This would be the number of individuals that domi-
%     nate p.
%   ?for each individual q in P
%       * if p dominated q then
%           ?add q to the set Sp i.e. Sp = Sp ? {q}
%       * else if q dominates p then
%           ?increment the domination counter for p i.e. np = np + 1
%   ?if np = 0 i.e. no individuals dominate p then p belongs to the first
%     front; Set rank of individual p to one i.e prank = 1. Update the first
%     front set by adding p to front one i.e F1 = F1 ? {p}
% ?This is carried out for all the individuals in main population P.
% ?Initialize the front counter to one. i = 1
% ?following is carried out while the ith front is nonempty i.e. Fi != []
%   ?Q = []. The set for storing the individuals for (i + 1)th front.
%   ?for each individual p in front Fi
%       * for each individual q in Sp (Sp is the set of individuals
%         dominated by p)
%           ?nq = nq-1, decrement the domination count for individual q.
%           ?if nq = 0 then none of the individuals in the subsequent
%             fronts would dominate q. Hence set qrank = i + 1. Update
%             the set Q with individual q i.e. Q = Q ? q.
%   ?Increment the front counter by one.
%   ?Now the set Q is the next front and hence Fi = Q.
%
% This algorithm is better than the original NSGA ([2]) since it utilize
% the informatoion about the set that an individual dominate (Sp) and
% number of individuals that dominate the individual (np).

for i = 1 : N
    % Number of individuals that dominate this individual
    individual(i).n = 0;
    % Individuals which this individual dominate
    individual(i).p = [];
    for j = 1 : N
        dom_less = 0;
        dom_equal = 0;
        dom_more = 0;
        for k = 1 : M
            if (x_temp(i,V + k) < x_temp(j,V + k))
                dom_less = dom_less + 1;
            elseif (x_temp(i,V + k) == x_temp(j,V + k))
                dom_equal = dom_equal + 1;
            else
                dom_more = dom_more + 1;
            end
        end
        if dom_less == 0 && dom_equal ~= M
            individual(i).n = individual(i).n + 1;
        elseif dom_more == 0 && dom_equal ~= M
            individual(i).p = [individual(i).p j];
        end
    end
    if individual(i).n == 0
        x_temp(i,M + V + 1) = 1;
        F(front).f = [F(front).f i];
    end
end
% Find the subsequent fronts
while ~isempty(F(front).f)
   Q.temp = []; % modified by zzb
   for i = 1 : length(F(front).f)
       if ~isempty(individual(F(front).f(i)).p)
            for j = 1 : length(individual(F(front).f(i)).p)
            	individual(individual(F(front).f(i)).p(j)).n = ...
                	individual(individual(F(front).f(i)).p(j)).n - 1;
                if individual(individual(F(front).f(i)).p(j)).n == 0
                    x_temp(individual(F(front).f(i)).p(j),M + V + 1) = front + 1;
                    Q.temp = [Q.temp individual(F(front).f(i)).p(j)]; % modified by zzb
                end
            end
       end
   end
   front =  front + 1;
   F(front).f = Q.temp; % modified by zzb
end

%% Crowding distance
%The crowing distance is calculated as below
% ?For each front Fi, n is the number of individuals.
%   ?initialize the distance to be zero for all the individuals i.e. Fi(dj ) = 0,
%     where j corresponds to the jth individual in front Fi.
%   ?for each objective function m
%       * Sort the individuals in front Fi based on objective m i.e. I =
%         sort(Fi,m).
%       * Assign infinite distance to boundary values for each individual
%         in Fi i.e. I(d1) = ? and I(dn) = ?
%       * for k = 2 to (n-1)
%           ?I(dk) = I(dk) + (I(k + 1).m ? I(k ? 1).m)/fmax(m) - fmin(m)
%           ?I(k).m is the value of the mth objective function of the kth
%             individual in I

% Find the crowding distance for each individual in each front
[~,index_of_fronts] = sort(x_temp(:,M + V + 1));
sorted_based_on_front = x_temp; % Initialize by x % modified by zzb
for i = 1 : length(index_of_fronts)
    sorted_based_on_front(i,:) = x_temp(index_of_fronts(i),:);
end
current_index = 0;

% for front = 1 : (length(F) - 1) % The last front of F is []
for front_count = 1 : (front - 1) % modified by zzb
    % objective = [];
    % distance = 0;
    % y = [];
    y = zeros(length(F(front_count).f),2*M + V + 1); % modified by zzb
    previous_index = current_index + 1;
    for i = 1 : length(F(front_count).f) % modified by zzb
        %y(i,:) = sorted_based_on_front(current_index + i,:);
        y(i,1:M + V + 1) = sorted_based_on_front(current_index + i,:); % modified by zzb
    end
    current_index = current_index + i;
    % Sort each individual based on the objective
    % sorted_based_on_objective = []; % modified by zzb
    for i = 1 : M
        [~, index_of_objectives] = sort(y(:,V + i)); % modified by zzb
        % sorted_based_on_objective = [];
        sorted_based_on_objective = y;  % Initialize by y  % modified by zzb
        for j = 1 : length(index_of_objectives)
            sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);
        end
        f_max = sorted_based_on_objective(length(index_of_objectives), V + i);
        f_min = sorted_based_on_objective(1, V + i);
        y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i) = Inf;
        y(index_of_objectives(1),M + V + 1 + i) = Inf;
         for j = 2 : length(index_of_objectives) - 1
            next_obj = sorted_based_on_objective(j + 1,V + i);
            previous_obj = sorted_based_on_objective(j - 1,V + i);
            if (f_max - f_min == 0)
                y(index_of_objectives(j),M + V + 1 + i) = Inf;
            else
                y(index_of_objectives(j),M + V + 1 + i) = ...
                     (next_obj - previous_obj)/(f_max - f_min);
            end
         end
    end
    % distance = [];
    % distance(:,1) = zeros(length(F(front_count).f),1); % modified by zzb
    distance = zeros(length(F(front_count).f),1); % modified by zzb
    for i = 1 : M
       % distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);% modified by zzb
       distance = distance + y(:,M + V + 1 + i);% modified by zzb
    end
    y(:,M + V + 2) = distance;
    % y = y(:,1 : M + V + 2);
    % z(previous_index:current_index,:) = y;
    z(previous_index:current_index,:) = y(:,1 : M + V + 2);
end
% f = z();
f = z;
%% References
% [1] *Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and T. Meyarivan*, |A Fast
% Elitist Multiobjective Genetic Algorithm: NSGA-II|, IEEE Transactions on
% Evolutionary Computation 6 (2002), no. 2, 182 ~ 197.
%
% [2] *N. Srinivas and Kalyanmoy Deb*, |Multiobjective Optimization Using
% Nondominated Sorting in Genetic Algorithms|, Evolutionary Computation 2
% (1994), no. 3, 221 ~ 248.

author: sfvsfv

repo: Mathematical-modeling


text:
%% ALS optimization for CP and Tucker tensor decompositions

%% Alternating least squares for PARAFAC/CANDECOMP
% The function |parafac_als| computes an estimate of the best rank-R
% PARAFAC model of a tensor X using an alternating least-squares
% algorithm.  The input X can be a tensor, sptensor, ktensor, or
% ttensor. The result P is a ktensor.
rand('state',0);
X = sptenrand([5 4 3], 10)
%%
P = parafac_als(X,2)
%%
P = parafac_als(X,2,struct('dimorder',[3 2 1]))
%%
P = parafac_als(X,2,struct('dimorder',[3 2 1],'init','nvecs'))
%%
U0 = {rand(5,2),rand(4,2),[]}; %<-- Initial guess for factors of P
P = parafac_als(X,2,struct('dimorder',[3 2 1],'init',{U0}))
%% Alternating least squares for Tucker model
% The function |tucker_als| computes the best rank(R1,R2,..,Rn)
% approximation of tensor X, according to the specified dimensions in
% vector R.  The input X can be a tensor, sptensor, ktensor, or
% ttensor.  The result returned in T is a ttensor.
X = sptenrand([5 4 3], 10)
%%
T = tucker_als(X,2)        %<-- best rank(2,2,2) approximation
%%
T = tucker_als(X,[2 2 1])  %<-- best rank(2,2,1) approximation
%%
T = tucker_als(X,2,struct('dimorder',[3 2 1]))
%%
T = tucker_als(X,2,struct('dimorder',[3 2 1],'init','eigs'))
%%
U0 = {rand(5,2),rand(4,2),[]}; %<-- Initial guess for factors of T
T = tucker_als(X,2,struct('dimorder',[3 2 1],'init',{U0}))

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%  INTERNAL FUNCTION
%
%  see also mygrid
%

author: jmaih