text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  internal function

%  'xform_nii.m' is an internal function called by "load_nii.m", so
%  you do not need run this program by yourself. It does simplified
%  NIfTI sform/qform affine transform, and supports some of the
%  affine transforms, including translation, reflection, and
%  orthogonal rotation (N*90 degree).
%
%  For other affine transforms, e.g. any degree rotation, shearing
%  etc. you will have to use the included 'reslice_nii.m' program
%  to reslice the image volume. 'reslice_nii.m' is not called by
%  any other program, and you have to run 'reslice_nii.m' explicitly
%  for those NIfTI files that you want to reslice them.
%
%  Since 'xform_nii.m' does not involve any interpolation or any
%  slice change, the original image volume is supposed to be
%  untouched, although it is translated, reflected, or even
%  orthogonally rotated, based on the affine matrix in the
%  NIfTI header.
%
%  However, the affine matrix in the header of a lot NIfTI files
%  contain slightly non-orthogonal rotation. Therefore, optional
%  input parameter 'tolerance' is used to allow some distortion
%  in the loaded image for any non-orthogonal rotation or shearing
%  of NIfTI affine matrix. If you set 'tolerance' to 0, it means
%  that you do not allow any distortion. If you set 'tolerance' to
%  1, it means that you do not care any distortion. The image will
%  fail to be loaded if it can not be tolerated. The tolerance will
%  be set to 0.1 (10%), if it is default or empty.
%
%  Because 'reslice_nii.m' has to perform 3D interpolation, it can
%  be slow depending on image size and affine matrix in the header.
%
%  After you perform the affine transform, the 'nii' structure
%  generated from 'xform_nii.m' or new NIfTI file created from
%  'reslice_nii.m' will be in RAS orientation, i.e. X axis from
%  Left to Right, Y axis from Posterior to Anterior, and Z axis
%  from Inferior to Superior.
%
%  NOTE: This function should be called immediately after load_nii.
%
%  Usage: [ nii ] = xform_nii(nii, [tolerance], [preferredForm])
%
%  nii	- NIFTI structure (returned from load_nii)
%
%  tolerance (optional) - distortion allowed for non-orthogonal rotation
%	or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
%	if it is default or empty.
%
%  preferredForm (optional)  -  selects which transformation from voxels
%	to RAS coordinates; values are s,q,S,Q.  Lower case s,q indicate
%	"prefer sform or qform, but use others if preferred not present".
%	Upper case indicate the program is forced to use the specificied
%	tranform or fail loading.  'preferredForm' will be 's', if it is
%	default or empty.	- Jeff Gunter
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function nii = xform_nii(nii, tolerance, preferredForm)

   %  save a copy of the header as it was loaded.  This is the
   %  header before any sform, qform manipulation is done.
   %
   nii.original.hdr = nii.hdr;

   if ~exist('tolerance','var') | isempty(tolerance)
      tolerance = 0.1;
   elseif(tolerance<=0)
      tolerance = eps;
   end

   if ~exist('preferredForm','var') | isempty(preferredForm)
      preferredForm= 's';				% Jeff
   end

   %  if scl_slope field is nonzero, then each voxel value in the
   %  dataset should be scaled as: y = scl_slope * x + scl_inter
   %  I bring it here because hdr will be modified by change_hdr.
   %
   if nii.hdr.dime.scl_slope ~= 0 & ...
	ismember(nii.hdr.dime.datatype, [2,4,8,16,64,256,512,768]) & ...
	(nii.hdr.dime.scl_slope ~= 1 | nii.hdr.dime.scl_inter ~= 0)

      nii.img = ...
	nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

      if nii.hdr.dime.datatype == 64

         nii.hdr.dime.datatype = 64;
         nii.hdr.dime.bitpix = 64;
      else
         nii.img = single(nii.img);

         nii.hdr.dime.datatype = 16;
         nii.hdr.dime.bitpix = 32;
      end

      nii.hdr.dime.glmax = max(double(nii.img(:)));
      nii.hdr.dime.glmin = min(double(nii.img(:)));

      %  set scale to non-use, because it is applied in xform_nii
      %
      nii.hdr.dime.scl_slope = 0;

   end

   %  However, the scaling is to be ignored if datatype is DT_RGB24.

   %  If datatype is a complex type, then the scaling is to be applied
   %  to both the real and imaginary parts.
   %
   if nii.hdr.dime.scl_slope ~= 0 & ...
	ismember(nii.hdr.dime.datatype, [32,1792])

      nii.img = ...
	nii.hdr.dime.scl_slope * double(nii.img) + nii.hdr.dime.scl_inter;

      if nii.hdr.dime.datatype == 32
         nii.img = single(nii.img);
      end

      nii.hdr.dime.glmax = max(double(nii.img(:)));
      nii.hdr.dime.glmin = min(double(nii.img(:)));

      %  set scale to non-use, because it is applied in xform_nii
      %
      nii.hdr.dime.scl_slope = 0;

   end

   %  There is no need for this program to transform Analyze data
   %
   if nii.filetype == 0 & exist([nii.fileprefix '.mat'],'file')
      load([nii.fileprefix '.mat']);	% old SPM affine matrix
      R=M(1:3,1:3);
      T=M(1:3,4);
      T=R*ones(3,1)+T;
      M(1:3,4)=T;
      nii.hdr.hist.qform_code=0;
      nii.hdr.hist.sform_code=1;
      nii.hdr.hist.srow_x=M(1,:);
      nii.hdr.hist.srow_y=M(2,:);
      nii.hdr.hist.srow_z=M(3,:);
   elseif nii.filetype == 0
      nii.hdr.hist.rot_orient = [];
      nii.hdr.hist.flip_orient = [];
      return;				% no sform/qform for Analyze format
   end

   hdr = nii.hdr;

   [hdr,orient]=change_hdr(hdr,tolerance,preferredForm);

   %  flip and/or rotate image data
   %
   if ~isequal(orient, [1 2 3])

      old_dim = hdr.dime.dim([2:4]);

      %  More than 1 time frame
      %
      if ndims(nii.img) > 3
         pattern = 1:prod(old_dim);
      else
         pattern = [];
      end

      if ~isempty(pattern)
         pattern = reshape(pattern, old_dim);
      end

      %  calculate for rotation after flip
      %
      rot_orient = mod(orient + 2, 3) + 1;

      %  do flip:
      %
      flip_orient = orient - rot_orient;

      for i = 1:3
         if flip_orient(i)
            if ~isempty(pattern)
               pattern = flipdim(pattern, i);
            else
               nii.img = flipdim(nii.img, i);
            end
         end
      end

      %  get index of orient (rotate inversely)
      %
      [tmp rot_orient] = sort(rot_orient);

      new_dim = old_dim;
      new_dim = new_dim(rot_orient);
      hdr.dime.dim([2:4]) = new_dim;

      new_pixdim = hdr.dime.pixdim([2:4]);
      new_pixdim = new_pixdim(rot_orient);
      hdr.dime.pixdim([2:4]) = new_pixdim;

      %  re-calculate originator
      %
      tmp = hdr.hist.originator([1:3]);
      tmp = tmp(rot_orient);
      flip_orient = flip_orient(rot_orient);

      for i = 1:3
         if flip_orient(i) & ~isequal(tmp(i), 0)
            tmp(i) = new_dim(i) - tmp(i) + 1;
         end
      end

      hdr.hist.originator([1:3]) = tmp;
      hdr.hist.rot_orient = rot_orient;
      hdr.hist.flip_orient = flip_orient;

      %  do rotation:
      %
      if ~isempty(pattern)
         pattern = permute(pattern, rot_orient);
         pattern = pattern(:);

         if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
		hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            tmp = reshape(nii.img(:,:,:,1), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,1) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            tmp = reshape(nii.img(:,:,:,2), [prod(new_dim) hdr.dime.dim(5:8)]);
            tmp = tmp(pattern, :);
            nii.img(:,:,:,2) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
               tmp = reshape(nii.img(:,:,:,3), [prod(new_dim) hdr.dime.dim(5:8)]);
               tmp = tmp(pattern, :);
               nii.img(:,:,:,3) = reshape(tmp, [new_dim       hdr.dime.dim(5:8)]);
            end

         else
            nii.img = reshape(nii.img, [prod(new_dim) hdr.dime.dim(5:8)]);
            nii.img = nii.img(pattern, :);
            nii.img = reshape(nii.img, [new_dim       hdr.dime.dim(5:8)]);
         end
      else
         if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792 | ...
		hdr.dime.datatype == 128 | hdr.dime.datatype == 511

            nii.img(:,:,:,1) = permute(nii.img(:,:,:,1), rot_orient);
            nii.img(:,:,:,2) = permute(nii.img(:,:,:,2), rot_orient);

            if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
               nii.img(:,:,:,3) = permute(nii.img(:,:,:,3), rot_orient);
            end
         else
            nii.img = permute(nii.img, rot_orient);
         end
      end
   else
      hdr.hist.rot_orient = [];
      hdr.hist.flip_orient = [];
   end

   nii.hdr = hdr;

   return;					% xform_nii

%-----------------------------------------------------------------------
function [hdr, orient] = change_hdr(hdr, tolerance, preferredForm)

   orient = [1 2 3];
   affine_transform = 1;

   %  NIFTI can have both sform and qform transform. This program
   %  will check sform_code prior to qform_code by default.
   %
   %  If user specifys "preferredForm", user can then choose the
   %  priority.					- Jeff
   %
   useForm=[];					% Jeff

   if isequal(preferredForm,'S')
       if isequal(hdr.hist.sform_code,0)
           error('User requires sform, sform not set in header');
       else
           useForm='s';
       end
   end						% Jeff

   if isequal(preferredForm,'Q')
       if isequal(hdr.hist.qform_code,0)
           error('User requires qform, qform not set in header');
       else
           useForm='q';
       end
   end						% Jeff

   if isequal(preferredForm,'s')
       if hdr.hist.sform_code > 0
           useForm='s';
       elseif hdr.hist.qform_code > 0
           useForm='q';
       end
   end						% Jeff

   if isequal(preferredForm,'q')
       if hdr.hist.qform_code > 0
           useForm='q';
       elseif hdr.hist.sform_code > 0
           useForm='s';
       end
   end						% Jeff

   if isequal(useForm,'s')
      R = [hdr.hist.srow_x(1:3)
           hdr.hist.srow_y(1:3)
           hdr.hist.srow_z(1:3)];

      T = [hdr.hist.srow_x(4)
           hdr.hist.srow_y(4)
           hdr.hist.srow_z(4)];

      if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
         hdr.hist.old_affine = [ [R;[0 0 0]] [T;1] ];
         R_sort = sort(abs(R(:)));
         R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
         hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

         if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
         end
      end

   elseif isequal(useForm,'q')
      b = hdr.hist.quatern_b;
      c = hdr.hist.quatern_c;
      d = hdr.hist.quatern_d;

      if 1.0-(b*b+c*c+d*d) < 0
         if abs(1.0-(b*b+c*c+d*d)) < 1e-5
            a = 0;
         else
            error('Incorrect quaternion values in this NIFTI data.');
         end
      else
         a = sqrt(1.0-(b*b+c*c+d*d));
      end

      qfac = hdr.dime.pixdim(1);
      i = hdr.dime.pixdim(2);
      j = hdr.dime.pixdim(3);
      k = qfac * hdr.dime.pixdim(4);

      R = [a*a+b*b-c*c-d*d     2*b*c-2*a*d        2*b*d+2*a*c
           2*b*c+2*a*d         a*a+c*c-b*b-d*d    2*c*d-2*a*b
           2*b*d-2*a*c         2*c*d+2*a*b        a*a+d*d-c*c-b*b];

      T = [hdr.hist.qoffset_x
           hdr.hist.qoffset_y
           hdr.hist.qoffset_z];

      %  qforms are expected to generate rotation matrices R which are
      %  det(R) = 1; we'll make sure that happens.
      %
      %  now we make the same checks as were done above for sform data
      %  BUT we do it on a transform that is in terms of voxels not mm;
      %  after we figure out the angles and squash them to closest
      %  rectilinear direction. After that, the voxel sizes are then
      %  added.
      %
      %  This part is modified by Jeff Gunter.
      %
      if det(R) == 0 | ~isequal(R(find(R)), sum(R)')

         %  det(R) == 0 is not a common trigger for this ---
         %  R(find(R)) is a list of non-zero elements in R; if that
         %  is straight (not oblique) then it should be the same as
         %  columnwise summation. Could just as well have checked the
         %  lengths of R(find(R)) and sum(R)' (which should be 3)
         %
         hdr.hist.old_affine = [ [R * diag([i j k]);[0 0 0]] [T;1] ];
         R_sort = sort(abs(R(:)));
         R( find( abs(R) < tolerance*min(R_sort(end-2:end)) ) ) = 0;
         R = R * diag([i j k]);
         hdr.hist.new_affine = [ [R;[0 0 0]] [T;1] ];

         if det(R) == 0 | ~isequal(R(find(R)), sum(R)')
            msg = [char(10) char(10) '   Non-orthogonal rotation or shearing '];
            msg = [msg 'found inside the affine matrix' char(10)];
            msg = [msg '   in this NIfTI file. You have 3 options:' char(10) char(10)];
            msg = [msg '   1. Using included ''reslice_nii.m'' program to reslice the NIfTI' char(10)];
            msg = [msg '      file. I strongly recommand this, because it will not cause' char(10)];
            msg = [msg '      negative effect, as long as you remember not to do slice' char(10)];
            msg = [msg '      time correction after using ''reslice_nii.m''.' char(10) char(10)];
            msg = [msg '   2. Using included ''load_untouch_nii.m'' program to load image' char(10)];
            msg = [msg '      without applying any affine geometric transformation or' char(10)];
            msg = [msg '      voxel intensity scaling. This is only for people who want' char(10)];
            msg = [msg '      to do some image processing regardless of image orientation' char(10)];
            msg = [msg '      and to save data back with the same NIfTI header.' char(10) char(10)];
            msg = [msg '   3. Increasing the tolerance to allow more distortion in loaded' char(10)];
            msg = [msg '      image, but I don''t suggest this.' char(10) char(10)];
            msg = [msg '   To get help, please type:' char(10) char(10) '   help reslice_nii.m' char(10)];
            msg = [msg '   help load_untouch_nii.m' char(10) '   help load_nii.m'];
            error(msg);
         end

      else
         R = R * diag([i j k]);
      end					% 1st det(R)

   else
      affine_transform = 0;	% no sform or qform transform
   end

   if affine_transform == 1
      voxel_size = abs(sum(R,1));
      inv_R = inv(R);
      originator = inv_R*(-T)+1;
      orient = get_orient(inv_R);

      %  modify pixdim and originator
      %
      hdr.dime.pixdim(2:4) = voxel_size;
      hdr.hist.originator(1:3) = originator;

      %  set sform or qform to non-use, because they have been
      %  applied in xform_nii
      %
      hdr.hist.qform_code = 0;
      hdr.hist.sform_code = 0;
   end

   %  apply space_unit to pixdim if not 1 (mm)
   %
   space_unit = get_units(hdr);

   if space_unit ~= 1
      hdr.dime.pixdim(2:4) = hdr.dime.pixdim(2:4) * space_unit;

      %  set space_unit of xyzt_units to millimeter, because
      %  voxel_size has been re-scaled
      %
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,1,0));
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,2,1));
      hdr.dime.xyzt_units = char(bitset(hdr.dime.xyzt_units,3,0));
   end

   return;					% change_hdr

%-----------------------------------------------------------------------
function orient = get_orient(R)

   orient = [];

   for i = 1:3
      switch find(R(i,:)) * sign(sum(R(i,:)))
      case 1
         orient = [orient 1];		% Left to Right
      case 2
         orient = [orient 2];		% Posterior to Anterior
      case 3
         orient = [orient 3];		% Inferior to Superior
      case -1
         orient = [orient 4];		% Right to Left
      case -2
         orient = [orient 5];		% Anterior to Posterior
      case -3
         orient = [orient 6];		% Superior to Inferior
      end
   end

   return;					% get_orient

%-----------------------------------------------------------------------
function [space_unit, time_unit] = get_units(hdr)

   switch bitand(hdr.dime.xyzt_units, 7)	% mask with 0x07
   case 1
      space_unit = 1e+3;		% meter, m
   case 3
      space_unit = 1e-3;		% micrometer, um
   otherwise
      space_unit = 1;			% millimeter, mm
   end

   switch bitand(hdr.dime.xyzt_units, 56)	% mask with 0x38
   case 16
      time_unit = 1e-3;			% millisecond, ms
   case 24
      time_unit = 1e-6;			% microsecond, us
   otherwise
      time_unit = 1;			% second, s
   end

   return;					% get_units

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% predicting numerical values using Logistic Regression
clear all;
format long
disp('===== Logistic Regression ====');
disp('Reading featur vector');

for feat = 1:3
    featurs = csvread('data\forWeka_featuresonly.csv');
    num_data = 50%size(featurs,1); %5000;
    disp(sprintf('Number of datapoints %d',num_data))

    possiblefeaturizations =  {'bernouli', 'tfidf','multinomial'};
    %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
    featurization  = possiblefeaturizations{feat}

    featurs = featurs(:,2:size(featurs,2));
    if strcmp(featurization,'multinomial')
        %just pass
    elseif strcmp(featurization,'bernouli')
        featurs = double(featurs>0);
    elseif strcmp(featurization,'tfidf')
        occurance = (featurs>0);
        idf = log(size(featurs,1)./sum(occurance));
        featurs = featurs.*repmat( idf, size(featurs,1),1);
    end

    size_training = floor(.8*num_data);

    trainingset = featurs(1:size_training,:);
    testset = featurs((size_training+1):num_data,:);

    disp('Splitting up data into training/test sets');
    [num,txt,raw] = xlsread('data\final104.xls');

    % reading the description of each shoe
    descriptions = raw(2:size(raw,1),2);
    style_ratings = num(1:size(num,1),1);
    comfort_ratings = num(1:size(num,1),4);
    overal_ratings = num(1:size(num,1),5);

    % only take m data points
    m=num_data;
    descriptions = descriptions(1:m);
    style_ratings = style_ratings(1:m);
    comfort_ratings = comfort_ratings(1:m);
    overal_ratings = overal_ratings(1:m);

    responsevals = [style_ratings, comfort_ratings, overal_ratings];

    responsevals_training = responsevals(1:size_training,:);
    responsevals_test = responsevals((size_training+1):num_data,:);

    disp('Logistic Regression');
    %

    tic;

    predictions = [];
    actual = [];
    for i =1:3
        a = responsevals_training(:,i);
        b = responsevals_test(:,i)';
        [COEFF,SCORE] = princomp(trainingset);
        regresscoeff = mnrfit(SCORE,a);
        C2 = (regresscoeff'*(testset'));
        predictions = [predictions, C2'];
        actual = [actual, b'];
    end

    MSE = mean(sum(((predictions-actual).^2)'))
    toc;

end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
% DMRG eigenvalue solver.
%   function [x,theta,testdata]=dmrg_eig(A, tol, varargin)
% Tries to solve the lowest eigenvalue (ground state) problem A*x=theta*x
% using the DMRG iteration.
% With default parameters, reconstructs the genuine 2-site DMRG for one
% targeted state, without enrichments.
% However, it can also work as the "block 1-site DMRG" from ref. [*], where
% at least 2 states must be targeted, but only one block at a time can be
% considered. See parameters 'b' and 'numblocks' below.
%
% A must be a Hermitian matrix in the tt_matrix format,
% tol is the relative truncation/stopping threshold,
% varargin may contain a sequence of tuning parameters of the form
% 'parameter1_name', parameter1_value, 'parameter2_name', parameter2_value
% and so on. Available parameters and default values are the following.
%   o b: Number of eigenpairs targeted (default: 1)
%   o numblocks: Number of TT blocks included into the current supercore.
%       Can be either 2 for 2-site dmrg, or 1 for 1-site dmrg (default: 2)
%   o nswp: maximal number of DMRG sweeps (default: 20)
%   o max_full_size: if the local system size (rx(i)*n(i)*rx(i+1)) is less
%       than max_full_size, the local problem is solved via eig, otherwise via
%       lobpcg (iteratively) (default: 50)
%   o local_iters: maximal number of the lobpcg iterations (default: 100)
%   o usemex: whether shall we use the optimized MEX library eig3d_primme
%       instead of lobpcg. It safely switches off automatically if
%       eig3d_primme is not found in the MATLAB path, as well as on complex
%       data. To obtain eig3d_primme, you need to compile it in the
%       TT-Toolbox/fmex directory, please follow instructions there (default: true)
%   o resid_damp: residual gap for the local solvers, i.e. the real
%       threshold for e.g. lobpcg is tol/(sqrt(d)*resid_damp) (default: 2)
%   o trunc_norm: truncation strategy: either the standard Frobenius norm
%       ('fro'), or the residual ('resid') (default: 'fro')
%   o rmax: maximal TT rank (bond dimension) allowed (default: Inf)
%   o tol_exit: stopping tolerance (default: tol)
%   o verb: verbosity level: silent (0), sweep info (1), block info (2) or
%       the full debug returned in the testdata output (3) (default: 1)
%   o kickrank: size of the random enrichment (default: 0)
%   o x0: Initial guess in the tt_tensor format (default: random rank-2)
%
% Output data are the solution x in the tt_tensor format,
% the eigenvalue theta, and
% (optionally) the history for debug (or convergence tests) purposes testdata.
%
%
% Implementation: Sergey Dolgov (sergey.v.dolgov@gmail.com)
%
% [*] Dolgov, Khoromskij, Oseledets, Savostyanov,
% "Computation of extreme eigenvalues in higher dimensions using block
% tensor train format", Comp. Phys. Comm. 2014, http://dx.doi.org/10.1016/j.cpc.2013.12.017
%

function [x,theta,testdata]=dmrg_eig(A, tol, varargin)
% Number of eigenstates
b = 1;
% Threshold on the local size between direct-iterative local solving
max_full_size=50;
% Number of local iters
local_iters = 100;
% Whether shall we use the local solver from the MEX library eig3d_primme
usemex = true;
% Number of sweeps
nswp=20;
% Residual gap for the local solver
resid_damp = 2;
% Truncation strategy
%trunc_norm = 'resid';
trunc_norm = 'fro';
% Max rank
rmax = Inf;
% Exit tol
tol_exit = tol;
% Verb
verb=1;
% Random enrichment rank
kickrank = 0;
% Number of blocks to consider: 2 or 1
numblocks = 2;
% Initial guess
x=[];

for i=1:2:length(varargin)-1
    switch lower(varargin{i})
        case 'nswp'
            nswp=varargin{i+1};
        case 'rmax'
            rmax=varargin{i+1};
        case 'x0'
            x=varargin{i+1};
        case 'verb'
            verb=varargin{i+1};
        case 'resid_damp'
            resid_damp=varargin{i+1};
        case 'trunc_norm'
            trunc_norm=varargin{i+1};
        case 'kickrank'
            kickrank=varargin{i+1};
        case  'max_full_size'
            max_full_size=varargin{i+1};
        case  'local_iters'
            local_iters=varargin{i+1};
        case 'usemex'
            usemex=varargin{i+1};
        case 'tol_exit'
            tol_exit=varargin{i+1};
        case 'b'
            b = varargin{i+1};
        case 'numblocks'
            numblocks = varargin{i+1};
        otherwise
            error('Unknown tuning parameter "%s"', varargin{i});
    end;
end;

% Extract the dimensions
if (~isa(A, 'tt_matrix'))
    error('A must be given in a tt_matrix class');
end;
d = A.d;
n = A.n;
% Initialize from random initial state if not passed otherwise
if (isempty(x))
    x = tt_rand(n, d, b, -1);
else
    if (~isa(x, 'tt_tensor'))
        error('x0 must be given in a tt_tensor class');
    end;
    if (x.d~=d)
        error('Inconsistent dim(A) and dim(x0)');
    end;
    if (any(x.n~=n))
        error('Inconsistent size(A) and size(x0)');
    end;
end;

% Disable MEX if it does not exist
if (usemex)&&(exist('eig3d_primme', 'file')<2)
    warning('MEX local solver is not found, disabled');
    usemex = false;
end;

% More housekeeping
ra = A.r;
rx = x.r;
if (rx(d+1)>1)
    % Leave only the first component from the initial guess
    x = x*eye(rx(d+1),1);
    rx(d+1) = 1;
end;
crA = core2cell(A);
crx = core2cell(x);

% Threshold for local problems
real_tol = (tol/sqrt(d))/resid_damp;

% Interfaces
phixax = cell(d+1,1); phixax{1}=1; phixax{d+1}=1;

% This is some convergence output for test purposes
testdata = cell(3,1);
testdata{1} = zeros(d-numblocks+1, nswp); % times
testdata{2} = zeros(d-numblocks+1, nswp, b); % evs
testdata{3} = zeros(d-numblocks+1, nswp); % local res or dx
t_dmrg_eig = tic;

% Presetup: compute the initial projections
for i=d:-1:2
    % QR (Gauge conditions) right-to-left
    ux = reshape(crx{i}, rx(i), n(i)*rx(i+1));
    [ux,vx]=qr(ux.', 0);
    cr2 = reshape(crx{i-1}, rx(i-1)*n(i-1), rx(i));
    cr2 = cr2*vx.';
    rx(i) = size(ux,2);
    crx{i} = reshape(ux.', rx(i), n(i), rx(i+1));
    crx{i-1} = reshape(cr2, rx(i-1), n(i-1), rx(i));

    % Pass the reduction to the left
    phixax(i) = rightreduce_matrix(phixax(i+1), crx{i}, crA(i), crx{i}, rx(i),n(i),rx(i+1), 1,ra(i),ra(i+1), rx(i),n(i),rx(i+1));
    if (usemex)&&(~isreal(phixax{i}))
        warning('Complex data detected, turning MEX local solver off');
        usemex = false;
    end;
end;
% Initial guess for EVs
theta = rightreduce_matrix(phixax(2), crx{1}, crA(1), crx{1}, rx(1),n(1),rx(1+1), 1,ra(1),ra(1+1), rx(1),n(1),rx(1+1));
theta = squeeze(theta{1});
if (usemex)&&(~isreal(theta))
    warning('Complex data detected, turning MEX local solver off');
    usemex = false;
end;
theta = diag(theta);
if (b>1)
    crx{1} = randn(rx(1), n(1), rx(2), b);
end;

% DMRG sweeps
swp = 1;
i = 1; % Start from the first block
dir = 1; % Current direction: forward(+1) or backward(-1)
% Meausure the errors and # of operations
max_dx = 0;
max_res = 0;
max_matvecs = 0;

while (swp<=nswp)||(dir>0)
    % Extract the matrix parts, accelerate a plenty of iterations with them
    % Phi1: 1, rx'1, rx1, ra1
    % Phi2: ra2, rx'2, rx2, 1
    Phi1 = phixax(i);
    A1 = crA(i);
    if (numblocks==1)
        % one-block DMRG, everything is ready
        Phi2 = phixax(i+1);
        rx1 = rx(i);
        nloc = n(i);
        rx2 = rx(i+1);
        ra1 = ra(i);
        ra2 = ra(i+1);
        % sol_prev. It is also an initial guess.
        sol_prev = reshape(crx{i}, rx(i)*n(i)*rx(i+1), b);
    else
        % Two-block DMRG, we have to merge some of two blocks into one.
        Phi2 = phixax(i+2);
        A2 = crA(i+1);
        if ((rx(i)*n(i))<=(n(i)*n(i+1)))&&((rx(i)*n(i))<=(n(i+1)*rx(i+2)))
            % Merge Phi1 and A1
            Phi1 = reshape(Phi1{1}, rx(i)*rx(i), ra(i));
            A1 = reshape(A1{1}, ra(i), n(i)*n(i)*ra(i+1));
            Phi1 = Phi1*A1;
            Phi1 = reshape(Phi1, rx(i), rx(i), n(i), n(i)*ra(i+1));
            Phi1 = permute(Phi1, [1,3,2,4]);
            Phi1 = {reshape(Phi1, rx(i)*n(i), rx(i)*n(i), ra(i+1))};
            A1 = A2;
            rx1 = rx(i)*n(i);
            nloc = n(i+1);
            rx2 = rx(i+2);
            ra1 = ra(i+1);
            ra2 = ra(i+2);
        elseif ((n(i)*n(i+1)<=(rx(i)*n(i))))&&((n(i)*n(i+1))<=(n(i+1)*rx(i+2)))
            % Merge A1 and A2
            A1 = reshape(A1{1}, ra(i)*n(i)*n(i), ra(i+1));
            A2 = reshape(A2{1}, ra(i+1), n(i+1)*n(i+1)*ra(i+2));
            A1 = A1*A2;
            A1 = reshape(A1, ra(i)*n(i), n(i), n(i+1), n(i+1)*ra(i+2));
            A1 = permute(A1, [1,3,2,4]);
            A1 = {reshape(A1, ra(i), n(i)*n(i+1), n(i)*n(i+1), ra(i+2))};
            rx1 = rx(i);
            nloc = n(i)*n(i+1);
            rx2 = rx(i+2);
            ra1 = ra(i);
            ra2 = ra(i+2);
        else
            % Merge A2 and Phi2
            Phi2 = reshape(Phi2{1}, ra(i+2), rx(i+2)*rx(i+2));
            A2 = reshape(A2{1}, ra(i+1)*n(i+1)*n(i+1), ra(i+2));
            Phi2 = A2*Phi2;
            Phi2 = reshape(Phi2, ra(i+1)*n(i+1), n(i+1), rx(i+2), rx(i+2));
            Phi2 = permute(Phi2, [1,3,2,4]);
            Phi2 = {reshape(Phi2, ra(i+1), n(i+1)*rx(i+2), n(i+1)*rx(i+2))};
            rx1 = rx(i);
            nloc = n(i);
            rx2 = n(i+1)*rx(i+2);
            ra1 = ra(i);
            ra2 = ra(i+1);
        end;
        % sol_prev. It is also an initial guess.
        if (dir>0)
            sol_prev = reshape(crx{i}, rx(i)*n(i)*rx(i+1), b);
            sol_prev = sol_prev.';
            sol_prev = reshape(sol_prev, b*rx(i)*n(i), rx(i+1));
            sol_prev = sol_prev*reshape(crx{i+1}, rx(i+1), n(i+1)*rx(i+2));
            sol_prev = reshape(sol_prev, b, rx(i)*n(i)*n(i+1)*rx(i+2));
            sol_prev = sol_prev.';
        else
            sol_prev = reshape(crx{i+1}, rx(i+1), n(i+1)*rx(i+2)*b);
            sol_prev = reshape(crx{i}, rx(i)*n(i), rx(i+1))*sol_prev;
            sol_prev = reshape(sol_prev, rx(i)*n(i)*n(i+1)*rx(i+2), b);
        end;
    end;

    % Initial residual
    res_prev = norm(local_matvec(sol_prev, rx1,nloc,rx2,b, rx1,nloc,rx2, Phi1, A1, Phi2, 1,ra1,ra2)-sol_prev*diag(theta));

    if (rx1*nloc*rx2<max_full_size) % Full solution
        %      |     |    |
        % B = Phi1 - A1 - Phi2
        %      |     |    |
        Bxx = assemble_local_matrix(Phi1, A1, Phi2, 1,ra1,ra2, rx1,nloc,rx2, rx1,nloc,rx2);
        Bxx = (Bxx+Bxx')*0.5; % Ensure the symmetry. At least now...
        [sol,L]=eig(Bxx);
        L = diag(L);
        L = real(L); % exclude possible i*1e-15 noise
        [~,ind] = sort(L, 'ascend');
        theta = L(ind(1:b));
        sol = sol(:,ind(1:b));
        num_matvecs = 1;
    else % Iterative solution
        if (usemex)
            Phi1m = reshape(Phi1{1}, rx1, rx1, ra1);
            Phi2m = reshape(Phi2{1}, ra2*rx2, rx2);
            Phi2m = Phi2m.';
            Phi2m = reshape(Phi2m, rx2, ra2, rx2);
            A1m = reshape(A1{1}, ra1, nloc, nloc, ra2);
            [sol,theta, num_matvecs]=eig3d_primme(Phi1m,A1m,Phi2m,real_tol,b,sol_prev,local_iters);
        else
            [sol,theta,~,lambdaHistory]=lobpcg(sol_prev,@(x)local_matvec(x, rx1,nloc,rx2,[], rx1,nloc,rx2, Phi1, A1, Phi2, 1,ra1,ra2),  real_tol, local_iters);
            num_matvecs = numel(lambdaHistory);
        end;
    end;

    % count the number of MatVecs
    max_matvecs = max(max_matvecs, num_matvecs);

    if (~strcmp(trunc_norm, 'fro'))
        % We need the new residual for the corresp. trunc. strategy
        res_new = norm(local_matvec(sol, rx1,nloc,rx2,b, rx1,nloc,rx2, Phi1, A1, Phi2, 1,ra1,ra2)-sol*diag(theta));
    end;

    % L2-norm convergence check.
    dx = norm(sol*(sol'*sol_prev)-sol_prev);
    max_dx = max(max_dx, dx);
    max_res = max(max_res, res_prev);

    if ((dir>0)&&(i<d))||((dir<0)&&(i>1))||(numblocks==2) % SVD and enrichment
        if (dir>0)
            if (numblocks==2)
                % In 2-block DMRG, we have to separate modes
                sol = reshape(sol, rx(i)*n(i), n(i+1)*rx(i+2)*b);
            else
                sol = reshape(sol, rx(i)*n(i), rx(i+1)*b);
            end;
        else
            sol = sol.';
            if (numblocks==2)
                % In 2-block DMRG, we have to separate modes
                sol = reshape(sol, b*rx(i)*n(i), n(i+1)*rx(i+2));
            else
                sol = reshape(sol, b*rx(i), n(i)*rx(i+1));
            end;
        end;

        [ux,s,vx]=svd(sol, 'econ');
        s = diag(s);

        if (strcmp(trunc_norm, 'fro')) % Fro-norm truncation
            r = my_chop2(s, real_tol*resid_damp*norm(s));
        else
            % Residual truncation
            % start from the old rank
            r = min([rx(i),rx(i+1),numel(s)]);
            cursol = ux(:,1:r)*diag(s(1:r))*vx(:,1:r)';
            if (dir>0)
                cursol = reshape(cursol, rx1*nloc*rx2, b);
            else
                cursol = reshape(cursol, b, rx1*nloc*rx2);
                cursol = cursol.';
            end;
            res = norm(local_matvec(cursol, rx1,nloc,rx2,b, rx1,nloc,rx2, Phi1, A1, Phi2, 1,ra1,ra2)-cursol*diag(theta));
            if (res<max(real_tol, res_new)*resid_damp)
                drank = -1; % rank is overestimated; decrease
            else
                drank = 1; % residual is large; increase the rank
            end;
            while (r>0)&&(r<=numel(s))
                cursol = ux(:,1:r)*diag(s(1:r))*vx(:,1:r)';
                if (dir>0)
                    cursol = reshape(cursol, rx1*nloc*rx2, b);
                else
                    cursol = reshape(cursol, b, rx1*nloc*rx2);
                    cursol = cursol.';
                end;
                res = norm(local_matvec(cursol, rx1,nloc,rx2,b, rx1,nloc,rx2, Phi1, A1, Phi2, 1,ra1,ra2)-cursol*diag(theta));
                if (drank>0)
                    if (res<max(real_tol, res_new)*resid_damp)
                        break;
                    end;
                else
                    if (res>=max(real_tol, res_new)*resid_damp)
                        break;
                    end;
                end;
                r = r+drank;
            end;
            if (drank<0)
                r=r+1;
            end;
        end;

        r = min(r, numel(s));
        r = min(r, rmax);

        if (verb==2)
            fprintf('=dmrg_eig= swp=%d, block=%d, dx=%3.3e, r=%d\n', swp, i, dx, r);
        end;
    end;

    if (dir>0)&&(i<d) % Forward sweep
        ux = ux(:,1:r);
        vx = conj(vx(:,1:r))*diag(s(1:r));

        if (kickrank>0)
            % Enrich the solution
            zx = randn(rx(i)*n(i), kickrank);
            % Concatenate the bases and reinforce the gauges.
            % In future: insert symmetries here?
            [ux,rv]=qr([ux,zx], 0);
            rv = rv(:,1:r);
            vx = vx*rv.';
        end;

        r = size(ux, 2);

        if (numblocks==1)
            % Drop the [r x r] factor to the next core -- White's prediction
            cr2 = crx{i+1};
            cr2 = reshape(cr2, rx(i+1), n(i+1)*rx(i+2));
            vx = reshape(vx, rx(i+1), b*r);
            cr2 = vx.'*cr2;
            cr2 = reshape(cr2, b, r*n(i+1)*rx(i+2));
            cr2 = cr2.';
        else
            % Replace the next core by vx
            cr2 = vx.';
            cr2 = reshape(cr2, r, n(i+1), rx(i+2), b);
        end;

        % Stuff them back
        rx(i+1) = r;
        crx{i} = reshape(ux, rx(i), n(i), rx(i+1));
        crx{i+1} = reshape(cr2, rx(i+1), n(i+1), rx(i+2), b);

        % Compute new reductions
        phixax(i+1) = leftreduce_matrix(phixax(i), ux, crA(i), ux, rx(i),n(i),rx(i+1), 1,ra(i),ra(i+1), rx(i),n(i),rx(i+1));
    elseif (dir<0)&&((i>1)||(numblocks==2)) % Backward sweep
        vx = conj(vx(:,1:r));
        ux = ux(:,1:r)*diag(s(1:r));

        if (kickrank>0)
            % Enrich the solution
            if (numblocks==1)
                zx = randn(n(i)*rx(i+1), kickrank);
            else
                zx = randn(n(i+1)*rx(i+2), kickrank);
            end;
            % Concatenate the bases and reinforce the gauges.
            % In future: insert symmetries here?
            [vx,rv]=qr([vx,zx], 0);
            rv = rv(:,1:r);
            ux = ux*rv.';
        end;

        r = size(vx, 2);

        if (numblocks==1)
            % Drop the [r x r] factor to the next core -- White's prediction
            cr2 = crx{i-1};
            cr2 = reshape(cr2, rx(i-1)*n(i-1), rx(i));
            ux = reshape(ux, b, rx(i)*r);
            ux = ux.';
            ux = reshape(ux, rx(i), r*b);
            cr2 = cr2*ux;
            cr2 = reshape(cr2, rx(i-1), n(i-1), r, b);
            % Stuff them back
            rx(i) = r;
            vx = vx.';
            crx{i} = reshape(vx, rx(i), n(i), rx(i+1));
            crx{i-1} = reshape(cr2, rx(i-1), n(i-1), rx(i), b);
            % Compute new reductions
            phixax(i) = rightreduce_matrix(phixax(i+1), vx, A1, vx, rx(i),n(i),rx(i+1), 1,ra(i),ra(i+1), rx(i),n(i),rx(i+1));
        else
            % Replace both cores
            vx = vx.';
            vx = reshape(vx, r, n(i+1), rx(i+2));
            ux = reshape(ux, b, rx(i)*n(i)*r);
            ux = ux.';
            rx(i+1) = r;
            crx{i+1} = vx;
            crx{i} = reshape(ux, rx(i), n(i), r, b);
            % Compute new reductions
            phixax(i+1) = rightreduce_matrix(phixax(i+2), vx, crA(i+1), vx, rx(i+1),n(i+1),rx(i+2), 1,ra(i+1),ra(i+2), rx(i+1),n(i+1),rx(i+2));
        end;
    else
        % Just stuff back the last core
        sol = reshape(sol, rx(i), n(i), rx(i+1), b);
        crx{i} = sol;
    end;

    if (verb>2)
        % Report the debug data if necessary
        if (dir>0)
            testdata{1}(i,swp) = toc(t_dmrg_eig);
            testdata{2}(i,swp,:) = theta;
            if (strcmp(trunc_norm, 'fro'))
                testdata{3}(i,swp) = dx;
            else
                testdata{3}(i,swp) = res_prev;
            end;
        else
            testdata{1}(d-numblocks+2-i,swp) = toc(t_dmrg_eig);
            testdata{2}(d-numblocks+2-i,swp,:) = theta;
            if (strcmp(trunc_norm, 'fro'))
                testdata{3}(d-numblocks+2-i,swp) = dx;
            else
                testdata{3}(d-numblocks+2-i,swp) = res_prev;
            end;
        end;
    end;

    i = i+dir;

    % Check for the end of the sweep
    if ((dir>0)&&(i>d-numblocks+1))||((dir<0)&&(i<1))
        % Report all
        if (verb>0)
            fprintf('=dmrg_eig= sweep %d, max_dx: %3.3e, max_res: %3.3e, erank: %g, theta: %3.15e, mv: %d\n', swp, max_dx, max_res, sqrt(rx(1:d)'*(n.*rx(2:d+1))/sum(n)), sum(theta), max_matvecs);
        end;

        % Check the stops
        if (strcmp(trunc_norm, 'fro'))
            if (max_dx<tol_exit)&&(verb<3)&&(dir>0)
                break;
            end;
        else
            if (max_res<tol_exit)&&(verb<3)&&(dir>0)
                break;
            end;
        end;

        swp = swp+1;
        % Go backward
        dir = -dir;
        i = i+dir;
        % Clear the errors
        max_dx = 0;
        max_res = 0;
        max_matvecs = 0;
    end;
end;

crx{d} = reshape(crx{d}, rx(d), n(d), b);
x = cell2core(tt_tensor, crx);

end

% Accumulates the left reduction W{1:k}'*A{1:k}*X{1:k}
function [WAX2] = leftreduce_matrix(WAX1, w, A, x, rw1,n,rw2, Ra,ra1,ra2, rx1,m,rx2)
% Left WAX has the form of the first matrix TT block, i.e. [rw, rx, ra]
WAX2 = WAX1;
wc = reshape(w, rw1, n*rw2);
xc = reshape(x, rx1*m, rx2);
for k=1:Ra
    WAX2{k} = reshape(WAX2{k}, rw1, rx1*ra1(k));
    WAX2{k} = wc'*WAX2{k}; % size n rw2 x rx1 ra1
    WAX2{k} = reshape(WAX2{k}, n, rw2*rx1*ra1(k));
    WAX2{k} = WAX2{k}.';
    WAX2{k} = reshape(WAX2{k}, rw2*rx1, ra1(k)*n);
    tmp = reshape(A{k}, ra1(k)*n, m*ra2(k));
    WAX2{k} = WAX2{k}*tmp; % size rw2 rx1 m ra2
    WAX2{k} = reshape(WAX2{k}, rw2, rx1*m*ra2(k));
    WAX2{k} = WAX2{k}.';
    WAX2{k} = reshape(WAX2{k}, rx1*m, ra2(k)*rw2);
    WAX2{k} = xc.'*WAX2{k}; % size rx2, ra2 rw2
    WAX2{k} = reshape(WAX2{k}, rx2*ra2(k), rw2);
    WAX2{k} = WAX2{k}.';
end;
end

% Accumulates the right reduction W{k:d}'*A{k:d}*X{k:d}
function [WAX1] = rightreduce_matrix(WAX2, w, A, x, rw1,n,rw2, Ra,ra1,ra2, rx1,m,rx2)
% Right WAX has the form of the last matrix TT block, i.e. [ra, rw, rx]
WAX1 = WAX2;
wc = reshape(w, rw1, n*rw2);
wc = conj(wc);
xc = reshape(x, rx1*m, rx2);
for k=1:Ra
    WAX1{k} = reshape(WAX1{k}, ra2(k)*rw2, rx2);
    WAX1{k} = xc*WAX1{k}.'; % size rx1 m x ra2 rw2
    WAX1{k} = reshape(WAX1{k}, rx1, m*ra2(k)*rw2);
    WAX1{k} = WAX1{k}.';
    WAX1{k} = reshape(WAX1{k}, m*ra2(k), rw2*rx1);
    tmp = reshape(A{k}, ra1(k)*n, m*ra2(k));
    WAX1{k} = tmp*WAX1{k}; % size ra1(k)*n, rw2*rx1
    WAX1{k} = reshape(WAX1{k}, ra1(k), n*rw2*rx1);
    WAX1{k} = WAX1{k}.';
    WAX1{k} = reshape(WAX1{k}, n*rw2, rx1*ra1(k));
    WAX1{k} = wc*WAX1{k}; % size rw1, rx1 ra1
    WAX1{k} = reshape(WAX1{k}, rw1*rx1, ra1(k));
    WAX1{k} = WAX1{k}.';
end;
end

% A matrix-vectors product for the matrix in the 3D TT (WAX1-A-WAX2), and
% full vectors of size (rx1*m*rx2) x b. Returns (rw1*n*rw2) x b
function [w]=local_matvec(x, rx1,m,rx2,b, rw1,n,rw2, WAX1, A, WAX2, Ra,ra1,ra2)
xc = reshape(x, rx1*m*rx2, []);
if (isempty(b))
    b = size(xc, 2);
end;
w = zeros(rw1*n*rw2, b);
xc = xc.';
xc = reshape(xc, b*rx1*m, rx2);
for k=1:Ra
    tmp = reshape(WAX2{k}, ra2(k)*rw2, rx2);
    wk = xc*tmp.';
    wk = reshape(wk, b*rx1, m*ra2(k)*rw2);
    wk = wk.';
    wk = reshape(wk, m*ra2(k), rw2*b*rx1);
    tmp = reshape(A{k}, ra1(k)*n, m*ra2(k));
    wk = tmp*wk;
    wk = reshape(wk, ra1(k)*n*rw2*b, rx1);
    wk = wk.';
    wk = reshape(wk, rx1*ra1(k), n*rw2*b);
    tmp = reshape(WAX1{k}, rw1, rx1*ra1(k));
    wk = tmp*wk;
    wk = reshape(wk, rw1*n*rw2, b);
    w = w+wk;
end;
end

% Builds the full (rw1*n*rw2) x (rx1*m*rx2) matrix from its TT blocks
function [B,sparseflag]=assemble_local_matrix(WAX1, A, WAX2, Ra,ra1,ra2, rw1,n,rw2, rx1,m,rx2)
% Check the sparsity of the matrix blocks
sparseflag = true;
for k=1:Ra
    if (~issparse(A{k}))
        sparseflag=false;
    end;
end;
if (sparseflag)
    B = sparse(rw2*rw1*n, rx2*rx1*m); % reverse order !!!
    % The reverse order is needed since usually the n x m part is large and
    % sparse, so let it be the senior dimension.
    % Note that currently only canonical sparse matrices are allowed
    for k=1:Ra
        tmp = reshape(WAX2{k}, rw2, rx2);
        tmp = sparse(tmp);
        Bk = reshape(WAX1{k}, rw1, rx1);
        Bk = sparse(Bk);
        Bk = kron(Bk, tmp); % mind endiannes
        Bk = kron(A{k}, Bk); % mind endiannes
        B = B+Bk;
    end;
else
    % There are dense blocks, everything is dense, and in the natural index
    % order
    B = zeros(rw1*n*rw2, rx1*m*rx2);
    for k=1:Ra
        Bk = reshape(WAX1{k}, rw1*rx1, ra1(k));
        tmp = reshape(A{k}, ra1(k), n*m*ra2(k));
        if (issparse(tmp))
            % Don't mix sparse if we are already full
            tmp = full(tmp);
        end;
        Bk = Bk*tmp;
        Bk = reshape(Bk, rw1, rx1, n, m*ra2(k));
        Bk = permute(Bk, [1,3,2,4]);
        Bk = reshape(Bk, rw1*n*rx1*m, ra2(k));
        tmp = reshape(WAX2{k}, ra2(k), rw2*rx2);
        Bk = Bk*tmp;
        Bk = reshape(Bk, rw1*n, rx1*m, rw2, rx2);
        Bk = permute(Bk, [1,3,2,4]);
        Bk = reshape(Bk, rw1*n*rw2, rx1*m*rx2);
        B = B+Bk;
    end;
end;
end

author: oseledets

repo: TT-Toolbox


text:
function [u,p,edge,A,eqn,info] = StokesP2P1(node,elem,pde,bdFlag,option)
%% STOKESP2P1 Stokes equation: P2-P1 Taylor-Hood elements.
%
%   [u,p] = STOKESP2P1(node,elem,pde,bdFlag) use quadratic and piceswise
%   linear elements to approximate velocity u and pressure p, repectively.
%
%       -div(mu*grad u) + grad p = f in \Omega,
%                        - div u = 0  in \Omega,
%   with
%       Dirichlet boundary condition        u = g_D  on \Gamma_D,
%       Neumann boundary condition du/dn - np = g_N  on \Gamma_N.
%
%   It is a choice of option.fem in Stokes. Please read Stokes for more
%   information on the input and output.
%
% See also Stokes, Poisson, StokesP2P1
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

tic;
%% Construct Data Structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1);  Nu = N+size(edge,1);   Np = N;

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix for Laplace operator
% generate sparse pattern
ii = zeros(21*NT,1); jj = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end
% quadrature points
if ~isfield(pde,'nu'), pde.nu = []; end
if ~isfield(option,'quadorder')
    % constant viscosity
    option.quadorder = 2;        % default order
    if ~isempty(pde.nu) && isnumeric(pde.nu) % numerical viscosity
        option.quadorder = 3;    % exact for linear diffusion coefficient
    end
end
[lambda, w] = quadpts(option.quadorder);
nQuad = size(lambda,1);
% compute non-zeros
sA = zeros(21*NT,nQuad);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    index = 0;
    for i = 1:6
        for j = i:6
            Aij = 0;
            if isempty(pde.nu) || isnumeric(pde.nu)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
            if ~isempty(pde.nu) && (pde.nu~=1)
                Aij = pde.nu*Aij;
            end
            Aij = Aij.*area;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Nu,Nu);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Nu,Nu);
A = A + AU + AU';
A = blkdiag(A,A);
clear Aij ii jj sA

%% Assemble the matrix for divergence operator
Dx = sparse(Np,Nu);
Dy = sparse(Np,Nu);
[lambda, w] = quadpts(2); % (div(P2), P1) is P2
nQuad = size(lambda,1);
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,1) = (4*lambda(p,1)-1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (4*lambda(p,2)-1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (4*lambda(p,3)-1).*Dlambda(:,:,3);
    Dphip(:,:,4) = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
    Dphip(:,:,5) = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
    Dphip(:,:,6) = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    for i = 1:6
        for j = 1:3
            Dxij = 0;
            Dyij = 0;
            Dxij = Dxij + w(p)*Dphip(:,1,i).*lambda(p,j);
            Dyij = Dyij + w(p)*Dphip(:,2,i).*lambda(p,j);
            Dx = Dx + sparse(elem(:,j),double(elem2dof(:,i)),Dxij.*area,Np,Nu);
            Dy = Dy + sparse(elem(:,j),double(elem2dof(:,i)),Dyij.*area,Np,Nu);
        end
    end
end
B = -[Dx Dy];
clear Dxij Dyij Dx Dy Dphip

%% Assemble right hand side
f1 = zeros(Nu,1);
f2 = zeros(Nu,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,weight] = quadpts(option.fquadorder);
    % basis values at quadrature points
    phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
    phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
    phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
    phi(:,4) = 4*lambda(:,2).*lambda(:,3);
    phi(:,5) = 4*lambda(:,3).*lambda(:,1);
    phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    nQuad = size(lambda,1);
    ft1 = zeros(NT,6);
    ft2 = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        % function values at quadrature points
        fp = pde.f(pxy);
        % evaluate fp outside.
        for j = 1:6
            ft1(:,j) = ft1(:,j) + fp(:,1).*phi(p,j)*weight(p);
            ft2(:,j) = ft2(:,j) + fp(:,2).*phi(p,j)*weight(p);
        end
    end
    ft1 = ft1.*repmat(area,1,6);
    ft2 = ft2.*repmat(area,1,6);
    f1 = accumarray(elem2dof(:),ft1(:),[Nu 1]);
    f2 = accumarray(elem2dof(:),ft2(:),[Nu 1]);
end
clear phi ft1 ft2

%% Boundary Conditions
[AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P1;

%% Record assembeling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(ufreeDof), return; end
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if length(f)+length(g) <= 1e3  % Direct solver for small size systems
        option.solver = 'direct';
    else          % Multigrid-type  solver for large size systems
        option.solver = 'asmg';
    end
end
solver = option.solver;

% solve
switch solver
    case 'direct'
        tic
        bigA = [AD, B'; ...
                B, sparse(Np,Np)];
        bigF = [f; g];
        bigu = [u; p];
        bigFreeDof = [ufreeDof; 2*Nu+pDof];
        bigu(bigFreeDof) = bigA(bigFreeDof,bigFreeDof)\bigF(bigFreeDof);
        u = bigu(1:2*Nu);
        p = bigu(2*Nu+1:end);
        residual = norm(bigF - bigA*bigu);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
%         option.tol = Np^(-2);
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                        u(ufreeDof),p,elem,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A(ufreeDof,ufreeDof),B(:,ufreeDof),f(ufreeDof),g,...
                                          u,p,node,elem,bdFlag,ufreeDof,option);
end

%% Post-process
if length(pDof) ~= Np % p is unique up to a constant
    % impose the condition int(p)=0
    c = sum(mean(p(elem),2).*area)/sum(area);
    p = p - c;
end

%% Output information
eqn = struct('A',AD,'B',BD,'f',f,'g',g,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdStokesP2P1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,BD,f,g,u,p,ufreeDof,pDof] = getbdStokesP2P1
    %% Boundary condition of Stokes equation: P2-P0 elements

    %% Initial set up
%     f = [f1; f2];    % set in Neumann boundary condition
    g = zeros(Np,1);
    u = zeros(2*Nu,1);
    p = zeros(Np,1);
    ufreeDof = (1:Nu)';
    pDof = (1:Np)';

    if ~exist('bdFlag','var'), bdFlag = []; end
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Find Dirichlet dof and modify the matrix
    % Find Dirichlet boundary dof: fixedDof and pDof
    isFixedDof = false(Nu,1);
    if ~isempty(bdFlag)       % case: bdFlag is not empty
        elem2edge = elem2dof(:,4:6)-N;
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;   % nodes of all D-edges
        isFixedDof(N + find(isDirichlet')) = true;% dof on D-edges
        fixedDof = find(isFixedDof);
        ufreeDof = find(~isFixedDof);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        fixedDof = bdDof;
        isFixedDof(fixedDof) = true;
        ufreeDof = find(~isFixedDof);
    end
    if isempty(fixedDof) % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        fixedDof = 1;
        ufreeDof = (2:Nu)';    % eliminate the kernel by enforcing u(1) = 0;
    end

    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % AD(fixedDof,fixedDof)=I, AD(fixedDof,ufreeDof)=0, AD(ufreeDof,fixedDof)=0.
    % BD(:,fixedDof) = 0 and thus BD'(fixedDof,:) = 0.
    bdidx = zeros(2*Nu,1);
    bdidx([fixedDof; Nu+fixedDof]) = 1;
    Tbd = spdiags(bdidx,0,2*Nu,2*Nu);
    T = spdiags(1-bdidx,0,2*Nu,2*Nu);
    AD = T*A*T + Tbd;
    BD = B*T;

    %% Part 2: Find boundary edges and modify the right hand side f and g
    % Find boundary edges: Neumann and Robin
    Neumann = []; Robin = []; %#ok<*NASGU>
    if ~isempty(bdFlag)
        isNeumann(elem2edge((bdFlag(:)==2)|(bdFlag(:) == 3))) = true;
        isRobin(elem2edge(bdFlag(:)==3)) = true;
        Neumannidx = find(isNeumann);
        Neumann   = edge(isNeumann,:);
        Robin     = edge(isRobin,:);
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        Neumann = edge(bdDof>N,:);
        if ~isempty(pde.g_R)
            Robin = Neumann;
        end
    end

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        [lambda,w] = quadpts1(3);
        nQuad = size(lambda,1);
        % quadratic bases (1---3---2)
        bdphi(:,1) = (2*lambda(:,1)-1).*lambda(:,1);
        bdphi(:,2) = (2*lambda(:,2)-1).*lambda(:,2);
        bdphi(:,3) = 4*lambda(:,1).*lambda(:,2);
        % length of edge
        ve = node(Neumann(:,1),:) - node(Neumann(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        % update RHS
        gex = zeros(size(Neumann,1),2);   % x-component
        gey = zeros(size(Neumann,1),2);   % y-component
        for pp = 1:nQuad
            pxy = lambda(pp,1)*node(Neumann(:,1),:)+lambda(pp,2)*node(Neumann(:,2),:);
            gp = pde.g_N(pxy);
            gex(:,1) = gex(:,1) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,1);
            gex(:,2) = gex(:,2) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,2);
            gey(:,1) = gey(:,1) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,1);
            gey(:,2) = gey(:,2) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,2);
            f1(N+Neumannidx) = f1(N+Neumannidx) + w(pp)*edgeLength.*gp(:,1)*bdphi(pp,3); % interior bubble
            f2(N+Neumannidx) = f2(N+Neumannidx) + w(pp)*edgeLength.*gp(:,2)*bdphi(pp,3); % interior bubble
        end
        f1(1:N) = f1(1:N) + accumarray(Neumann(:), gex(:),[N,1]);
        f2(1:N) = f2(1:N) + accumarray(Neumann(:), gey(:),[N,1]);
    end
    f = [f1; f2];
    % The case non-empty Neumann but g_N=[] corresponds to the zero flux
    % boundary condition on Neumann edges and no modification is needed.

    % Dirichlet boundary conditions
    if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
        u1 = zeros(Nu,1);
        u2 = zeros(Nu,1);
        idx = (fixedDof > N);              % index of edge dof
        uD = pde.g_D(node(fixedDof(~idx),:));  % bd value at vertex dofs
        u1(fixedDof(~idx)) = uD(:,1);
        u2(fixedDof(~idx)) = uD(:,2);
        bdEdgeIdx = fixedDof(idx)-N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:)+node(edge(bdEdgeIdx,2),:))/2;
        uD = pde.g_D(bdEdgeMid);         % bd values at middle points of edges
        u1(fixedDof(idx)) = uD(:,1);
        u2(fixedDof(idx)) = uD(:,2);
        u = [u1; u2]; % Dirichlet bd condition is built into u
        f = f - A*u;  % bring affect of nonhomgenous Dirichlet bd condition to
        g = g - B*u;  % the right hand side
        g = g - mean(g);
        f(fixedDof) = u1(fixedDof);
        f(fixedDof+Nu) = u2(fixedDof);
    end
    % The case non-empty Dirichlet but g_D=[] corresponds to the zero Dirichlet
    % boundary condition and no modification is needed.

    % modfiy pressure dof for pure Dirichlet
    if isempty(Neumann)
        pDof = (1:Np-1)';
    end

    ufreeDof = [ufreeDof; Nu+ufreeDof];
    end % end of function getbdStokesP2P1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

repo: Aero-matlab


text:
function Labels = mri_getlabels_svreg()
% ATLAS     : BrainSuite SVREG (Brainsuite1, USCBrain)
% REFERENCE : http://brainsuite.org/uscbrain-description/

% @=============================================================================
% 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.
% =============================================================================@

Labels = {...
       0, 'Background',                                         [    0    0    0]; ...
       1, 'CSF',                                                [  255    0    0]; ...
       2, 'GM',                                                 [    0  255    0]; ...
       3, 'WM',                                                 [  255  255    0]; ...
       4, 'GM/CSF',                                             [    0    0  255]; ...
       5, 'GM/WM',                                              [  255    0  255]; ...
       6, 'CSF/other',                                          [  255  255  255]; ...
      16, 'Scalp',                                              [  122  157  221]; ...
      17, 'Skull',                                              [  176  176  176]; ...
      18, 'Space',                                              [    0    0  255]; ...
      19, 'Brain',                                              [  128    0   64]; ...
     100, 'Frontal lobe R',                                     [  149   46  145]; ...
     101, 'Frontal lobe L',                                     [  242   56  240]; ...
     110, 'Lateral frontal lobe R',                             [  175  236   13]; ...
     111, 'Lateral frontal lobe L',                             [   80   19  242]; ...
     120, 'Superior frontal gyrus R',                           [    0  255  127]; ...
     121, 'Superior frontal gyrus L',                           [  255    0  128]; ...
     130, 'Middle frontal gyrus R',                             [  128  255  255]; ...
     131, 'Middle frontal gyrus L',                             [  127    0    0]; ...
     140, 'Inferior frontal gyrus R',                           [  175  236   13]; ...
     141, 'Inferior frontal gyrus L',                           [   80   19  242]; ...
     142, 'Pars opercularis R',                                 [  255  255    0]; ...
     143, 'Pars opercularis L',                                 [    0    0  255]; ...
     144, 'Pars triangularis R',                                [  127    0  128]; ...
     145, 'Pars triangularis L',                                [  128  255  127]; ...
     146, 'Pars orbitalis R',                                   [  255  255    0]; ...
     147, 'Pars orbitalis L',                                   [    0    0  255]; ...
     150, 'Pre-central gyrus R',                                [  128  255    0]; ...
     151, 'Pre-central gyrus L',                                [  127    0  255]; ...
     160, 'Pole/orbital frontal lobe R',                        [  175  236   13]; ...
     161, 'Pole/orbital frontal lobe L',                        [   80   19  242]; ...
     162, 'Transvers frontal gyrus R',                          [  255  255    0]; ...
     163, 'Transvers frontal gyrus L',                          [    0    0  255]; ...
     164, 'Gyrus rectus R',                                     [  255    0  255]; ...
     165, 'Gyrus rectus L',                                     [    0  255    0]; ...
     166, 'Middle orbito-frontal gyrus R',                      [  128  255    0]; ...
     167, 'Middle orbito-frontal gyrus L',                      [  127    0  255]; ...
     168, 'Anterior orbito-frontal gyrus R',                    [  255    0  128]; ...
     169, 'Anterior orbito-frontal gyrus L',                    [    0  255  127]; ...
     170, 'Posterior orbito-frontal gyrus R',                   [  128  255  255]; ...
     171, 'Posterior orbito-frontal gyrus L',                   [  127    0    0]; ...
     172, 'Lateral orbitofrontal gyrus R',                      [  128    0  255]; ...
     173, 'Lateral orbitofrontal gyrus L',                      [  127  255    0]; ...
     180, 'Mesial frontal lobe R',                              [  175  236   13]; ...
     181, 'Mesial frontal lobe L',                              [   80   19  242]; ...
     182, 'Paracentral lobule R',                               [  255    0    0]; ...
     183, 'Paracentral lobule L',                               [    0  255  255]; ...
     184, 'Cingulate gyrus R',                                  [  255  255  128]; ...
     185, 'Cingulate gyrus L',                                  [    0    0  127]; ...
     186, 'Subcallosal gyrus R',                                [  255    0    0]; ...
     187, 'Subcallosal gyrus L',                                [    0  255  255]; ...
     200, 'Parietal lobe R',                                    [  117  236   61]; ...
     201, 'Parietal lobe L',                                    [  171  242   49]; ...
     220, 'Lateral parietal lobe R',                            [  175  236   13]; ...
     221, 'Lateral parietal lobe L',                            [   80   19  242]; ...
     222, 'Post-central gyrus R',                               [  128  127    0]; ...
     223, 'Post-central gyrus L',                               [  127  128  255]; ...
     224, 'Supramarginal gyrus R',                              [    0  128  127]; ...
     225, 'Supramarginal gyrus L',                              [  255  127  128]; ...
     226, 'Angular gyrus R',                                    [  255  255  127]; ...
     227, 'Angular gyrus L',                                    [    0    0  128]; ...
     228, 'Superior parietal gyrus R',                          [    0  128  255]; ...
     229, 'Superior parietal gyrus L',                          [  255  127    0]; ...
     240, 'Mesial parietal lobe R',                             [  175  236   13]; ...
     241, 'Mesial parietal lobe L',                             [   80   19  242]; ...
     242, 'Pre-cuneus R',                                       [    0  255    0]; ...
     243, 'Pre-cuneus L',                                       [  255    0  255]; ...
     300, 'Temporal lobe R',                                    [   50  132  236]; ...
     301, 'Temporal lobe L',                                    [   93  239  242]; ...
     310, 'Temporal pole R',                                    [  255    0    0]; ...
     311, 'Temporal pole L',                                    [    0  255  255]; ...
     320, 'Lateral/superior temporal lobe R',                   [  128  255  127]; ...
     321, 'Lateral/superior temporal lobe R',                   [  127    0  128]; ...
     322, 'Superior temporal gyrus R',                          [  128  255  127]; ...
     323, 'Superior temporal gyrus L',                          [  127    0  128]; ...
     324, 'Transverse temporal gyrus R',                        [  255    0  128]; ...
     325, 'Transverse temporal gyrus L',                        [    0  255  127]; ...
     326, 'Middle temporal gyrus R',                            [  255    0  255]; ...
     327, 'Middle temporal gyrus L',                            [    0  255    0]; ...
     328, 'Inferior temporal gyrus R',                          [  128  255  255]; ...
     329, 'Inferior temporal gyrus L',                          [  127    0    0]; ...
     330, 'Fusiforme gyrus R',                                  [  255  255    0]; ...
     331, 'Fusiforme gyrus L',                                  [    0    0  255]; ...
     340, 'Infero/mesial temporal lobe R',                      [    0    0  255]; ...
     341, 'Infero/mesial temporal lobe L',                      [  255  255    0]; ...
     342, 'Parahippocampal gyrus R',                            [    0    0  255]; ...
     343, 'Parahippocampal gyrus L',                            [  255  255    0]; ...
     344, 'Hippocampus R',                                      [  255  128  127]; ...
     345, 'Hippocampus L',                                      [    0  127  128]; ...
     346, 'Amygdala R',                                         [  255  128    0]; ...
     347, 'Amygdala L',                                         [    0  127  255]; ...
     370, 'Hippocampus R',                                      [  255  128  127]; ...
     371, 'Hippocampus L',                                      [    0  127  128]; ...
     400, 'Occipital lobe R',                                   [  255  252  122]; ...
     401, 'Occipital lobe L',                                   [  242  225   42]; ...
     420, 'Lateral occipital lobe R',                           [  175  236   13]; ...
     421, 'Lateral occipital lobe L',                           [   80   19  242]; ...
     422, 'Superior occipital gyrus R',                         [  128  127    0]; ...
     423, 'Superior occipital gyrus L',                         [  127  128  255]; ...
     424, 'Middle occipital gyrus R',                           [  128  255    0]; ...
     425, 'Middle occipital gyrus L',                           [  127    0  255]; ...
     440, 'Infero/mesial occipital lobe R',                     [  175  236   13]; ...
     441, 'Infero/mesial occipital lobe L',                     [   80   19  242]; ...
     442, 'Inferior occipital gyrus R',                         [  127    0  255]; ...
     443, 'Inferior occipital gyrus L',                         [  128  255    0]; ...
     444, 'Lingual gyrus R',                                    [  127    0  127]; ...
     445, 'Lingual gyrus L',                                    [  128  255  128]; ...
     446, 'Cuneus R',                                           [  255  128  255]; ...
     447, 'Cuneus L',                                           [    0  127    0]; ...
     500, 'Insula R',                                           [  127  255  128]; ...
     501, 'Insula L',                                           [  128    0  127]; ...
     600, 'Central gray structures R',                          [  175  236   13]; ...
     601, 'Central gray structures L',                          [   80   19  242]; ...
     610, 'Basal ganglia R',                                    [  175  236   13]; ...
     611, 'Basal ganglia L',                                    [   80   19  242]; ...
     612, 'Caudate nucleus R',                                  [  127  128  255]; ...
     613, 'Caudate nucleus L',                                  [  128  127    0]; ...
     614, 'Putamen R',                                          [  255  255    0]; ...   % OLD: 630
     615, 'Putamen L',                                          [    0    0  255]; ...   % OLD: 631
     616, 'Pallidum R',                                         [  128    0  127]; ...   % OLD: 640
     617, 'Pallidum L',                                         [  127  255  128]; ...   % OLD: 641
     620, 'Accumbens R',                                        [    0  255  255]; ...
     621, 'Accumbens L',                                        [  255    0    0]; ...
     630, 'Claustrum R',                                        [  255    0    0]; ...
     631, 'Claustrum L',                                        [    0  255  255]; ...
     640, 'Thalamus R',                                         [    0  255  255]; ...   % OLD: 660
     641, 'Thalamus L',                                         [  255    0    0]; ...   % OLD: 661
     650, 'Basal forebrain R',                                  [    0  128    0]; ...
     651, 'Basal forebrain L',                                  [  255  127  255]; ...
     660, 'Lateral geniculate nucleus R',                       [  128    0    0]; ...
     661, 'Lateral geniculate nucleus L',                       [  127  255  255]; ...
     662, 'Medial geniculate nucleus R',                        [  127  127  255]; ...
     663, 'Medial geniculate nucleus L',                        [  128  128    0]; ...
     670, 'Superior colliculus R',                              [  255  128    0]; ...
     671, 'Superior colliculus L',                              [    0  127  255]; ...
     680, 'Inferior colliculus R',                              [  255    0  128]; ...
     681, 'Inferior colliculus L',                              [    0  255  127]; ...
     690, 'Mamillary body R',                                   [  127  255  127]; ...
     691, 'Mamillary body L',                                   [  128    0  128]; ...
     700, 'Ventricular system R',                               [  175  236   13]; ...
     701, 'Ventricular system L',                               [   80   19  242]; ...
     710, 'Pineal gland',                                       [  255  255    0]; ...
     720, 'Lateral ventricle R',                                [  255  255    0]; ...
     721, 'Lateral ventricle L',                                [    0    0  255]; ...
     740, 'Third ventricle',                                    [    0    0  127]; ...
     760, 'Fourth ventricle',                                   [    0    0  255]; ...
     780, 'Cerebral aqueduct',                                  [  175  236   13]; ...
     800, 'Brainstem',                                          [  127    0  127]; ...
     850, 'Corpus Callosum',                                    [  255    0    0]; ...
     900, 'Cerebellum',                                         [  128   64   64]; ...
    4000, 'Central sulcus R',                                   [   18    0  242]; ...
    4001, 'Central sulcus L',                                   [   18    0  242]; ...
    4002, 'Precentral sulcus R',                                [   35  235   61]; ...
    4003, 'Precentral sulcus L',                                [   35  235   61]; ...
    4004, 'Superior frontal sulcus R',                          [  245  125  249]; ...
    4005, 'Superior frontal sulcus L',                          [  245  125  249]; ...
    4006, 'Inferio frontal sulcus R',                           [    0  255  255]; ...
    4007, 'Inferio frontal sulcus L',                           [    0  255  255]; ...
    4008, 'Ascending branch of the Sylvian fissure R',          [  236   55  173]; ...
    4009, 'Ascending branch of the Sylvian fissure L',          [  236   55  173]; ...
    4010, 'Horizontal branch of the Sylvian fissure R',         [    0  122  244]; ...
    4011, 'Horizontal branch of the Sylvian fissure L',         [    0  122  244]; ...
    4012, 'Diagonal sulcus R',                                  [  192  128  255]; ...
    4013, 'Diagonal sulcus L',                                  [  192  128  255]; ...
    4014, 'Lateral orbital sulcus R',                           [   17  238  193]; ...
    4015, 'Lateral orbital sulcus L',                           [   17  238  193]; ...
    4016, 'Frontomarginal sulcus R',                            [  125  190  255]; ...
    4017, 'Frontomarginal sulcus L',                            [  125  190  255]; ...
    4018, 'Cingulate sulcus R',                                 [  170  253  251]; ...
    4019, 'Cingulate sulcus L',                                 [  170  253  251]; ...
    4020, 'Paracentral sulcus R',                               [   52  180   78]; ...
    4021, 'Paracentral sulcus L',                               [   52  180   78]; ...
    4022, 'Superior supra orbital sulcus R',                    [  255  255    0]; ...
    4023, 'Superior supra orbital sulcus L',                    [  255  255    0]; ...
    4024, 'Inferior supra orbital sulcus R',                    [  255  255    0]; ...
    4025, 'Inferior supra orbital sulcus L',                    [  255  255    0]; ...
    4026, 'Olfactory or medial orbital sulcus R',               [  125  190  255]; ...
    4027, 'Olfactory or medial orbital sulcus L',               [  125  190  255]; ...
    4028, 'H shaped sulci, mesial R',                           [  255    0    0]; ...
    4029, 'H shaped sulci, mesial L',                           [  255    0    0]; ...
    4030, 'H shaped sulci, lateral R',                          [  255    0    0]; ...
    4031, 'H shaped sulci, lateral L',                          [  255    0    0]; ...
    4032, 'H shaped sulci, transverse R',                       [  255    0    0]; ...
    4033, 'H shaped sulci, transverse L',                       [  255    0    0]; ...
    4034, 'Sylvian fissure terminal split R',                   [  255  255    0]; ...
    4035, 'Sylvian fissure terminal split L',                   [  255  255    0]; ...
    4036, 'Superior temporal sulcus R',                         [   17  238  193]; ...
    4037, 'Superior temporal sulcus L',                         [   17  238  193]; ...
    4038, 'Inferior temporal sulcus R',                         [  245  152   90]; ...
    4039, 'Inferior temporal sulcus L',                         [  245  152   90]; ...
    4040, 'Occipito temporal sulcus R',                         [   92    0  185]; ...
    4041, 'Occipito temporal sulcus L',                         [   92    0  185]; ...
    4042, 'Collateral or medial occipito temporal sulcus R',    [  217    0  255]; ...
    4043, 'Collateral or medial occipito temporal sulcus L',    [  217    0  255]; ...
    4044, 'Transverse temporal sulcus R',                       [   35  235   61]; ...
    4045, 'Transverse temporal sulcus L',                       [   35  235   61]; ...
    4046, 'Circular sulcus R',                                  [  217    0  255]; ...
    4047, 'Circular sulcus L',                                  [  217    0  255]; ...
    4048, 'Postcentral sulcus R',                               [   35  235   61]; ...
    4049, 'Postcentral sulcus L',                               [   35  235   61]; ...
    4050, 'Intraparietal sulcus R',                             [  245  125  249]; ...
    4051, 'Intraparietal sulcus L',                             [  245  125  249]; ...
    4052, 'Primary sulcus of Jensen R',                         [    0  255  255]; ...
    4053, 'Primary sulcus of Jensen L',                         [    0  255  255]; ...
    4054, 'Secondary sulcus of Jensen R',                       [    0  255  255]; ...
    4055, 'Secondary sulcus of Jensen L',                       [    0  255  255]; ...
    4056, 'Parieto occipital sulcus R',                         [  128  128    0]; ...
    4057, 'Parieto occipital sulcus L',                         [  128  128    0]; ...
    4058, 'Subparietal sulcus R',                               [    0  255  255]; ...
    4059, 'Subparietal sulcus L',                               [    0  255  255]; ...
    4060, 'Calcarine sulcus R',                                 [   18    0  242]; ...
    4061, 'Calcarine sulcus L',                                 [   18    0  242]; ...
    4062, 'Calcarine sulcus terminal T R',                      [  255    0  255]; ...
    4063, 'Calcarine sulcus terminal T L',                      [  255    0  255]; ...
    4064, 'Transverse occipital sulcus superior R',             [   17  238  193]; ...
    4065, 'Transverse occipital sulcus superior L',             [   17  238  193]; ...
    4066, 'Superior lateral occipital sulcus R',                [   92    0  185]; ...
    4067, 'Superior lateral occipital sulcus L',                [   92    0  185]; ...
    4068, 'Inferior lateral occipital sulcus R',                [   92    0  185]; ...
    4069, 'Inferior lateral occipital sulcus L',                [   92    0  185]; ...
    4070, 'Anterior occipital sulcus R',                        [    0  122  244]; ...
    4071, 'Anterior occipital sulcus L',                        [    0  122  244]; ...
};

% Duplicate labels < 1000:  GM=1000+label, WM=2000+label
iDup = find([Labels{:,1}] < 1000);
Labels = cat(1, Labels, ...
    cat(2, cellfun(@(c)plus(1000, c), Labels(iDup,1), 'UniformOutput', 0), cellfun(@(c)cat(2, 'GM ', c), Labels(iDup,2), 'UniformOutput', 0), Labels(iDup,3)), ...
    cat(2, cellfun(@(c)plus(2000, c), Labels(iDup,1), 'UniformOutput', 0), cellfun(@(c)cat(2, 'WM ', c), Labels(iDup,2), 'UniformOutput', 0), Labels(iDup,3)));

author: brainstorm-tools


text:
function res = bf_write_gifti(BF, S)
% Writes out beamformer results as GIfTI meshes
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_write_gifti.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    normalise         = cfg_menu;
    normalise.tag     = 'normalise';
    normalise.name    = 'Global normalisation';
    normalise.help    = {'Normalise image values by the mean'};
    normalise.labels  = {
        'No'
        'Each image separately'
        'Across images'
        }';
    normalise.values  = {
        'no'
        'separate'
        'all'
        }';
    normalise.val = {'no'};

    space         = cfg_menu;
    space.tag     = 'space';
    space.name    = 'Image space';
    space.help    = {'Specify image space'};
    space.labels  = {
        'MNI'
        'Native'
        'MNI-aligned'
        }';
    space.values  = {
        'mni'
        'native'
        'aligned'
        }';
    space.val = {'mni'};

    visualise = cfg_menu;
    visualise.tag = 'visualise';
    visualise.name = 'Visualise the outputs';
    visualise.help = {'Visualise the outputs in the graphics window.'};
    visualise.labels = {'inflated', 'yes', 'no'};
    visualise.values = {2, 1, 0};
    visualise.val = {0};

    res      = cfg_branch;
    res.tag  = 'gifti';
    res.name = 'GIfTI';
    res.val  = {normalise, space, visualise};

    return
elseif nargin < 2
    error('Two input arguments are required');
end

if ~isfield(BF.sources, 'mesh')
    error('Only mesh source space is supported');
end

scale = ones(1, numel(BF.output.image));
switch S.normalise
    case  'separate'
        for i = 1:numel(BF.output.image)
            val = BF.output.image(i).val;
            scale(i) = 1./mean(abs(val(~isnan(val))));
        end
    case  'all'
        val = spm_vec({BF.output.image(:).val});
        scale = scale./mean(abs(val(~isnan(val))));
end

switch S.space
    case 'mni'
        source = BF.sources.mesh.canonical;
    case 'aligned'
        source = BF.sources.mesh.individual;
        source.vert =  spm_eeg_inv_transform_points(BF.data.transforms.toMNI_aligned, source.vert);
    case 'native'
        source = BF.sources.mesh.individual;
        source.vert =  spm_eeg_inv_transform_points(BF.data.transforms.toNative, source.vert);
end

save(gifti(source), [BF.data.D.fname '.surf.gii']);

nimages = numel(BF.output.image);

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nimages , 'Writing out images'); drawnow;
if nimages  > 100, Ibar = floor(linspace(1, nimages ,100));
else Ibar = 1:nimages; end

for i = 1:nimages
    fname = fullfile(pwd, [BF.output.image(i).label '.gii']);

    G = gifti;
    G.private.metadata(1).name = 'SurfaceID';
    G.private.metadata(1).value = [BF.data.D.fname '.surf.gii'];

    G.cdata = scale(i)*BF.output.image(i).val;

    G.cdata = G.cdata(:);

    save(G, fname, 'ExternalFileBinary');

    res.files{i, 1} = fname;

    if S.visualise
        Fgraph  = spm_figure('GetWin', fname); figure(Fgraph); clf
        a = axes;
        if S.visualise == 2
            H = spm_mesh_render('Disp', spm_mesh_inflate(gifti([BF.data.D.fname '.surf.gii'])), 'Parent', a);
            H = spm_mesh_render('Overlay', H, gifti(fname));
        else
            H = spm_mesh_render('Disp',gifti([BF.data.D.fname '.surf.gii']), 'Parent', a);
            H = spm_mesh_render('Overlay', H, gifti(fname));
        end
         spm_mesh_render('Colormap', H, jet(256));
    end

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

author: spm


text:
close all; clear all; clc;

% Create the directory for storing results
spx.fs.ensure_dir('bin');

h = spx.data.motion.Hopkins155;
% pre-load all examples
h.load_all_examples();
h.describe();
examples = h.get_2_3_motions();
ne = length(examples);
fprintf('Number of 2, 3 motions: %d\n', ne);
min_angles = zeros(1, ne);
max_angles = zeros(1, ne);
all_ranks = zeros(3, ne);
for i=1:ne
    fprintf('Example (%d): ', i);
    example = examples{i};
    fprintf('%s, %d motions, %d points, %d frames;', example.name, example.num_motions, example.num_points, example.num_frames);
    % The dataset
    X = example.X;
    % average norm
    norms1 = spx.norm.norms_l2_cw(X);
    fprintf('\n');
    %X = spx.la.affine.homogenize(X, 1);
    norms2 = spx.norm.norms_l2_cw(X);
    m1 = mean(norms1);
    m2 = mean(norms2);
    d = (m2 - m1) * 100 / m1 ;
    %fprintf('Mean norm: before: %0.2f, after: %0.2f, diff: %0.2f %%\n', m1, m2, d);
    % now compute the bases for 4 dimensional approximations
    % bases  = spx.la.svd.low_rank_bases(X, example.counts, 2);
    [bases, ranks]  = spx.la.svd.vidal_rank_bases(X, example.counts, 0.1);
    all_ranks(1:example.num_motions, i) = ranks;
    fprintf('Ranks: ');
    fprintf('%d ', ranks);
    fprintf('\n');
    % bases  = spx.la.spaces.bases(X, example.counts);
    angles = spx.la.spaces.smallest_angles_deg(bases);
    off_diag_angles = spx.matrix.off_diag_upper_tri_elements(angles);
    off_diag_angles = sort(off_diag_angles);
    minimum_angle = min(off_diag_angles);
    maximum_angle = max(off_diag_angles);
    fprintf('Angles between subspaces (degrees): \n');
    disp(angles);
    % disp(off_diag_angles');
    % fprintf('Minimum angle: %.2f degrees\n\n', minimum_angle);
    % fprintf('Maximum angle: %.2f degrees\n\n', maximum_angle);
    fprintf('\n');
    min_angles(i) = minimum_angle;
    max_angles(i) = maximum_angle;
end
fprintf ('Minimum angle: %0.2f deg\n', min(min_angles));
fprintf ('Maximum angle: %0.2f deg\n', max(max_angles));
save ('bin/principal_angles','min_angles', 'max_angles', 'all_ranks');

author: indigits

repo: sparse-plex


text:
function [TB] = EB2TB(EB)
% Convert computery things from exabytes to terabytes.
% Chad A. Greene 2012
TB = EB*1048576 ;

author: Sable

repo: mcbench-benchmarks


text:
function varargout = plot_spline(P,C)
  % PLOT_SPLINE Plot a cubic Bezier spline.
  %
  % [pe,p] = plot_spline(P,C)
  %
  % Inputs:
  %   P  #P by dim list of control point locations
  %   C  #C by 4 list of indices into P of cubic Bezier curves
  % Outputs:
  %   pe  plot handles for UI
  %   p  plot handle for curves
  % Example:
  %   cellfun(@(pe) arrayfun(@(p) set(p,'Color','r'),pe),plot_spline(P,C))
  assert(max(C(:))<=size(P,1));
  assert(min(C(:))>=1);
  p = {};
  pe = {};
  ish = ishold;
  for c = 1:size(C,1)
    [pe{c},p{c}] = plot_cubic(P(C(c,:),:));
    hold on;
  end
  hold off;
  if ish
    hold on
  end
  if nargout>=1
    varargout{1} = pe;
    if nargout >= 2
      varargout{2} = p;
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [ p, t ] = distmesh_2d ( fd, fh, h0, box, iteration_max, pfix, ...
  varargin )

%*****************************************************************************80
%
%% DISTMESH_2D is a 2D mesh generator using distance functions.
%
%  Example:
%
%    Uniform Mesh on Unit Circle:
%
%      fd = inline('sqrt(sum(p.^2,2))-1','p');
%      [p,t] = distmesh_2d(fd,@huniform,0.2,[-1,-1;1,1],100,[]);
%
%    Rectangle with circular hole, refined at circle boundary:
%
%      fd = inline('ddiff(drectangle(p,-1,1,-1,1),dcircle(p,0,0,0.5))','p');
%      fh = inline('min(4*sqrt(sum(p.^2,2))-1,2)','p');
%      [p,t] = distmesh_2d(fd,fh,0.05,[-1,-1;1,1],500,[-1,-1;-1,1;1,-1;1,1]);
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    09 June 2012
%
%  Author:
%
%    Per-Olof Persson
%    Modifications by John Burkardt
%
%  Reference:
%
%    Per-Olof Persson, Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, function FD, signed distance function d(x,y).
%
%    Input, function FH, scaled edge length function h(x,y).
%
%    Input, real H0, the initial edge length.
%
%    Input, real BOX(2,2), the bounding box [xmin,ymin; xmax,ymax].
%
%    Input, integer ITERATION_MAX, the maximum number of iterations.
%    The iteration might terminate sooner than this limit, if the program decides
%    that the mesh has converged.
%
%    Input, real PFIX(NFIX,2), the fixed node positions.
%
%    Input, VARARGIN, aditional parameters passed to FD and FH.
%
%    Output, real P(N,2), the node positions.
%
%    Output, integer T(NT,3), the triangle indices.
%
%  Local parameters:
%
%    Local, real GEPS, a tolerance for determining whether a point is "almost" inside
%    the region.  Setting GEPS = 0 makes this an exact test.  The program currently
%    sets it to 0.001 * H0, that is, a very small multiple of the desired side length
%    of a triangle.  GEPS is also used to determine whether a triangle falls inside
%    or outside the region.  In this case, the test is a little tighter.  The centroid
%    PMID is required to satisfy FD ( PMID ) <= -GEPS.
%
  dptol = 0.001;
  ttol = 0.1;
  Fscale = 1.2;
  deltat = 0.2;
  geps = 0.001 * h0;
  deps = sqrt ( eps ) * h0;
  iteration = 0;
  triangulation_count = 0;
%
%  1. Create the initial point distribution by generating a rectangular mesh
%  in the bounding box.
%
  [ x, y ] = meshgrid ( box(1,1) : h0           : box(2,1), ...
                        box(1,2) : h0*sqrt(3)/2 : box(2,2) );
%
%  Shift the even rows of the mesh to create a "perfect" mesh of equilateral triangles.
%  Store the X and Y coordinates together as our first estimate of "P", the mesh points
%  we want.
%
  x(2:2:end,:) = x(2:2:end,:) + h0 / 2;
  p = [ x(:), y(:) ];
%
%  Instead of a regular mesh, you can initialize P with random values here.
%
%  2. Remove mesh points that are outside the region,
%  then satisfy the density constraint.
%
%  Keep only points inside (or almost inside) the region, that is, FD(P) < GEPS.
%
  p = p( feval ( fd, p, varargin{:} ) < geps, : );
%
%  Set R0, the relative probability to keep a point, based on the mesh density function.
%
  r0 = 1 ./ feval ( fh, p, varargin{:} ).^2;
%
%  Apply the rejection method to thin out points according to the density.
%
  p = [ pfix; p(rand(size(p,1),1) < r0 ./ max ( r0 ),: ) ];
%
%  The following obscure commands are written so that:
%  * We ALWAYS keep the fixed points at the beginning of the P array.
%  * We remove any points which are duplicates of these points.
%  That way, we do not later allow the fixed points to move,
%  because we know they are at the beginning of the array,
%  and DELAUNAYN stops complaining about duplicate points.
%  JVB, 09 June 2012.
%
  [ q, i, j ] = unique ( p, 'rows', 'first' );
  k = unique ( i );
  p = p(k,:);

  N = size ( p, 1 );
%
%  If ITERATION_MAX is 0, we're almost done.
%  For just this case, do the triangulation, then exit.
%  Setting ITERATION_MAX to 0 means that we can see the initial mesh
%  before any of the improvements have been made.
%
  if ( iteration_max <= 0 )
    t = delaunayn ( p );
    triangulation_count = triangulation_count + 1;
    return
  end

  pold = inf;

  while ( iteration < iteration_max )

    iteration = iteration + 1;

    if ( mod ( iteration, 10 ) == 0 )
      fprintf ( 1, '  %d iterations, %d triangulations\n', ...
        iteration, triangulation_count );
    end
%
%  3. Retriangulation by the Delaunay algorithm.
%
%  Was there large enough movement to retriangulate?
%
%  If so, save the current positions, get the list of
%  Delaunay triangles, compute the centroids, and keep
%  the interior triangles (whose centroids are within the region).
%
    if ( ttol < max ( sqrt ( sum ( ( p - pold ).^2, 2 ) ) / h0 ) )
      N = size ( p, 1 );
      pold = p;
      t = delaunayn ( p );
      triangulation_count = triangulation_count + 1;
      pmid = ( p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:) ) / 3;
      t = t( feval ( fd, pmid, varargin{:} ) < -geps, : );
%
%  4. Describe each bar by a unique pair of nodes.
%
      bars = [ t(:,[1,2]); t(:,[1,3]); t(:,[2,3]) ];
      bars = unique ( sort ( bars, 2 ), 'rows' );
%
%  5. Graphical output of the current mesh.
%
      trimesh ( t, p(:,1), p(:,2), zeros(N,1), 'EdgeColor', 'b', 'Linewidth', 2 )
      title ( sprintf ( 'Iteration %d', iteration ) );
      view(2), axis equal, axis off, drawnow
%
%  Put a "pause" command here if you'd like to see each new mesh.
%
    end
%
%  6. Move mesh points based on bar lengths L and forces F.
%
%  Make a list of the bar vectors and lengths.
%  Set L0 to the desired lengths, F to the scalar bar forces,
%  and FVEC to the x, y components of the bar forces.
%
%  At the fixed positions, reset the force to 0.
%
    barvec = p(bars(:,1),:) - p(bars(:,2),:);
    L = sqrt ( sum ( barvec.^2, 2 ) );
    hbars = feval ( fh, (p(bars(:,1),:)+p(bars(:,2),:))/2, varargin{:} );
    L0 = hbars * Fscale * sqrt ( sum(L.^2) / sum(hbars.^2) );
    F = max ( L0 - L, 0 );
    Fvec = F ./ L * [1,1] .* barvec;
    Ftot = full ( sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2) );
    Ftot(1:size(pfix,1),:) = 0;
    p = p + deltat * Ftot;
%
%  7. Bring outside points back to the boundary.
%
%  Use the numerical gradient of FD to project points back to the boundary.
%
    d = feval ( fd, p, varargin{:} );
    ix = d > 0;
    dgradx = ( feval(fd,[p(ix,1)+deps,p(ix,2)],varargin{:}) - d(ix) ) / deps;
    dgrady = ( feval(fd,[p(ix,1),p(ix,2)+deps],varargin{:}) - d(ix) ) / deps;
    p(ix,:) = p(ix,:) - [ d(ix) .* dgradx, d(ix) .* dgrady ];
%
%  8. Termination criterion: All interior nodes move less than dptol (scaled).
%
    if ( max ( sqrt ( sum ( deltat * Ftot ( d < -geps,:).^2, 2 ) ) / h0 ) < dptol )
      break;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% MAIN.m  --  Particle Swarm Optimization
%
% This script is used to run particle swarm optimization on several test
% functions.
%

quadraticBowl = @(x)( sum(x.^2,1) );
quadraticBowlSkew2d = @(x)( 0.3*x(1,:).^2 + 4.2*x(2,:).^2 );
noisyBowl = @(x)( sum(x.^2 + 0.1*randn(size(x)),1) );

problem.options.populationCount = 15;   % Number of particles in the search
problem.options.w = 0.3;  % Particle damping coefficient
problem.options.pg = 0.7;   % Global search coefficient
problem.options.pl = 0.5;   % Local search coefficient
problem.options.maxIter = 20;  % Number of iterations

problem.xLow = -ones(2,1);
problem.xUpp = ones(2,1);

problem.objFun = quadraticBowl;
% problem.objFun = quadraticBowlSkew2d;
% problem.objFun = noisyBowl;

soln = particleSwarmOptimization(problem);

%%%% Plotting:
figure(1);
for i=1:problem.options.maxIter
    plotIterData(soln,i);   %Only works for 2D objective function
    pause(0.5);
end

author: MatthewPeterKelly

repo: dscTutorials


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%Relationship between FS coefficients

%b(k)=a(k)+a(-k)
%c(k)=j[a(k)-a(-k)]
syms t ;
x=exp(-t);
t0=0;
T=5;
w=2*pi/T;
k=-6:6;
n=1:6 ;
b=(2/T)*int(x*cos(n*w*t),t,t0,t0+T);
b=eval(b)
c=(2/T)*int(x*sin(n*w*t),t,t0,t0+T);
c=eval(c)
a=(1/T)*int(x*exp(-j*k*w*t), t,t0,t0+T);
a=eval(a);
for i=1:6
bb(7-i)=a(14-i)+a(i);
cc(7-i)=j*(a(14-i)-a(i));
end
bb
cc

% a(k)=0.5[b(k)-jc(k)]
an=(1/2)*(b-j*c)
ak(1:6)=a(8:13)