text:
function KeplerUniversalVsSTK()
%% Constants
mu = 398600.4418;

%% Polar GEO Altitude Orbit Propagation Comparison
load('PolarGEO.mat');
r0 = repmat(PosVel(1,1:3)',[1 size(PosVel,1)]);
v0 = repmat(PosVel(1,4:6)',[1 size(PosVel,1)]);
[r,v] = keplerUniversal(r0,v0,0:60:(size(PosVel,1)-1)*60,mu);
PlotStateVectorErrors(r,v,PosVel,0:60:(size(PosVel,1)-1)*60,'GEO Altitude, Polar Inclination Propagation');

%% LEO Orbit Propagation Comparison
load('InclinedLEO.mat');
r0 = repmat(PosVel(1,1:3)',[1 size(PosVel,1)]);
v0 = repmat(PosVel(1,4:6)',[1 size(PosVel,1)]);
[r,v] = keplerUniversal(r0,v0,0:60:(size(PosVel,1)-1)*60,mu);
PlotStateVectorErrors(r,v,PosVel,0:60:(size(PosVel,1)-1)*60,'LEO Inclined Propagation');

%% Hyperbolic LEO Altitude Orbit Propagtion Comparison
load('HyperbolicLEO.mat');
r0 = repmat(PosVel(1,1:3)',[1 size(PosVel,1)]);
v0 = repmat(PosVel(1,4:6)',[1 size(PosVel,1)]);
[r,v] = keplerUniversal(r0,v0,0:60:(size(PosVel,1)-1)*60,mu);
PlotStateVectorErrors(r,v,PosVel,0:60:(size(PosVel,1)-1)*60,'Hyperbolic Orbit Propagation');

%% Parabolic LEO Altitude Orbit Propagation comparison
load('ParabolicLEO.mat');
r0 = repmat(PosVel(1,1:3)',[1 size(PosVel,1)]);
v0 = repmat(PosVel(1,4:6)',[1 size(PosVel,1)]);
[r,v] = keplerUniversal(r0,v0,0:60:(size(PosVel,1)-1)*60,mu);
PlotStateVectorErrors(r,v,PosVel,0:60:(size(PosVel,1)-1)*60,'Parabolic Orbit Propagation');

end

function PlotStateVectorErrors(r,v,PosVel,t,titleText)
MaxPositionError = abs(r-PosVel(:,1:3)');
MaxVelocityError = abs(v-PosVel(:,4:6)');
figure('color',[1 1 1]);
subplot(4,1,1); plot(t./60,MaxPositionError); title(titleText); ylabel('Pos Err (km)');
subplot(4,1,2); plot(t./60,MaxVelocityError); ylabel('Vel Err (km/s)');
subplot(4,1,3); plot(t./60,abs(sqrt(sum(r.^2,1)) - sqrt(sum(PosVel(:,1:3)'.^2,1)))); ylabel('Pos Mag Err (km)'); xlabel('Time (m)');
subplot(4,1,4); plot(t./60,abs(sqrt(sum(v.^2,1)) - sqrt(sum(PosVel(:,4:6)'.^2,1)))); ylabel('Vel Mag Err (km/s)'); xlabel('Time (m)');
end


text:
function [theta,P,res]=JPLAY(X,Y,G,L,k,d,sigma,alfa,beta,gama,rho,maxiter,eta,epsilon)

% Joint & Progressive Learning Strategy (J-PLAY)

% Usage:
%       [theta,P] = JPLAY(X,Y,G,L,k,d,sigma,alfa,beta,gama,maxiter,eta)

%% Input:
%       X      -Input data: d*N (d: dimension, N: the number of sample)
%       Y      -Label matrix: c*N (c: the number of class),
%                             e.g. if class==1 then y=[1,0,0,...,0]'
%                                  if class==2 then y=[0,1,0,...,0]'
%       G      -Adjacency matrix: N*N
%       L      -Laplacian matrix: N*N
%       k      -The number of nearest neighbor
%       d      -Reduced dimension
%     sigma    -Gaussian kernel parameter used in constructing graph
%     alfa     -Reconstruction loss parameter (default: 1)
%     beta     -Manifold regularization parameter (default: 0.1)
%     gama     -Ridge regression parameter  (default: 0.1)
%    maxiter   -Maximum number of iterations (default: 1000)
%     eta      -Eta in AutoRULe (default: 0.1)

%% Output:
%     theta    -A series of coupled-projection matrices
%       P      -Property-labeled projection matrix

% (c) Danfeng Hong, Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Germany.
%                   Singnal Processing in Earth Oberservation (SiPEO), Technical University of Munich (TUM), Germany.
%     Naoto Yokoya, Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Germany.
%                   Singnal Processing in Earth Oberservation (SiPEO), Technical University of Munich (TUM), Germany.
%                   Research Center for Advanced Science and Technology, The Univerisity of Tokyo, Japan.
%     danfeng.hong@dlr.de
%     naoto.yokoya@dlr.de; yokoya@sal.rcast.u-tokyo.ac.jp

if nargin < 2
    error('Please specify label information!')
end;

if nargin < 3
    error('Please construct adjacency matrix!')
end;

if nargin < 4
    error('Please construct Laplacian matrix!')
end;

if nargin < 5
    error('Please specify the number of nearest neighbor!')
end;

if nargin < 6
    error('Please specify subspace dimensions!')
end;

if nargin < 7
    error('Please specify Gaussian kernel parameter used in constructing graph!')
end;

if nargin < 8 || isempty(alfa)
    alfa=1;
end;

if nargin < 9 || isempty(beta)
    beta=0.1;
end;

if nargin < 10 || isempty(gama)
    gama=0.1;
end;

if nargin < 11 || isempty(maxiter)
    maxiter=1000;
end;

if nargin < 12 || isempty(eta)
    eta=0.1;
end;

%% Parameters Setting
% epsilon = 1e-3; % Tolerance error 1e-3Houston2018
iter=1;
stop = false;
res=zeros(1,maxiter); % Residuals
Num=length(d); % The number of layers

%% JPL: Initialization-step
theta=cell(1,1);
X0=X;
for i=1:Num
       theta0=DR_LPP(X0,k,d(i),sigma,G);
        theta_init=AutoRULe(theta0'*X0,X0,L,eta,theta0',maxiter);
        X0=theta_init*X0;
        theta{1,i}=theta_init;
end

%% JPL: Fine-tuning parameters
while ~stop && iter < maxiter+1

    %% Solve W
    D=X;
    for i=1:Num
        D=theta{1,i}*D;
    end
    P=(alfa*(Y*D'))/(alfa*(D*D')+gama*eye(size(D*D')));

    %% Solve the group of theta
    for j=1:Num

        %give W
        Pi=P;
        if j<Num
            for z=Num:-1:j+1
                Pi=Pi*theta{1,z};
            end
        end

        %give H
        H=X;
        for m=1:j
            H=theta{1,m}*H;
        end

        %give X
        Xi=X;
%         XWi=XW;
        if j>1
           for n=1:j-1
               Xi=theta{1,n}*Xi;
%                XWi=theta{1,n}*XWi;
           end
        end

        theta{1,j}=Theta_ADMM(Y,Pi,H,Xi,L,alfa,beta,theta{1,j},rho,maxiter);
    end

    %% Compute the vaules of objection function
     ErrorTerm=X;
%      ErrorTerm1=XW;
     ManifoldTerm=0;
%      ManifoldTerm1=0;
     ReconstructionTerm=0;
     for r=1:Num
         ReconstructionTerm=ReconstructionTerm+(norm(ErrorTerm-theta{1,r}'*theta{1,r}*ErrorTerm,'fro')^2);
         ErrorTerm=theta{1,r}*ErrorTerm;
%          ErrorTerm1=theta{1,r}*ErrorTerm1;
         ManifoldTerm=ManifoldTerm+trace(ErrorTerm*L*ErrorTerm');
%          ManifoldTerm1=ManifoldTerm1+trace(ErrorTerm1*L1*ErrorTerm1');
     end
     res(1,iter)=0.5*alfa*(norm(Y-P*ErrorTerm,'fro')^2)+0.5*beta*(ManifoldTerm)+0.5*gama*(norm(P,'fro')^2)+0.5*(ReconstructionTerm);

    %% Check the convergence condition
    if iter>1
       r_Obj=abs(res(1,iter)-res(1,iter-1))/res(1,iter-1);
       if r_Obj<epsilon
            stop = true;
            fprintf(' i = %f,res_Obj= %f\n',iter,r_Obj);
            break;
       end

       if mod(iter,10) == 1
           fprintf(' i = %f,res_Obj= %f\n',iter,r_Obj);
       end
    end

    iter=iter+1;
end
end

author: BehnoodRasti


text:
function [lf] = ft_compute_leadfield(dippos, sens, headmodel, varargin)

% FT_COMPUTE_LEADFIELD computes a forward solution for a dipole in a a volume
% conductor model. The forward solution is expressed as the leadfield
% matrix (Nchan*3), where each column corresponds with the potential or field
% distributions on all sensors for one of the x,y,z-orientations of the
% dipole.
%
% Use as
%   [lf] = ft_compute_leadfield(dippos, sens, headmodel, ...)
% with input arguments
%   dippos    = position dipole (1*3 or Ndip*3)
%   sens      = structure with gradiometer or electrode definition
%   headmodel = structure with volume conductor definition
%
% The headmodel represents a volume conductor model, its contents
% depend on the type of model. The sens structure represents a sensor
% array, i.e. EEG electrodes or MEG gradiometers.
%
% It is possible to compute a simultaneous forward solution for EEG and MEG
% by specifying sens and grad as two cell-arrays, e.g.
%   sens       = {senseeg, sensmeg}
%   headmodel  = {voleeg,  volmeg}
% This results in the computation of the leadfield of the first element of
% sens and headmodel, followed by the second, etc. The leadfields of the
% different imaging modalities are subsequently concatenated.
%
% Additional input arguments can be specified as key-value pairs, supported
% optional arguments are
%   'reducerank'      = 'no' or number
%   'normalize'       = 'no', 'yes' or 'column'
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or 1xN vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%   'backproject'     = 'yes' (default) or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace
%
% The leadfield weight may be used to specify a (normalized)
% corresponding surface area for each dipole, e.g. when the dipoles
% represent a folded cortical surface with varying triangle size.
%
% Depending on the specific input arguments for the sensor and volume, this
% function will select the appropriate low-level EEG or MEG forward model.
% The leadfield matrix for EEG will have an average reference over all the
% electrodes.
%
% The supported forward solutions for MEG are
%   infinite homogenous medium
%   single sphere (Cuffin and Cohen, 1977)
%   multiple spheres with one sphere per channel (Huang et al, 1999)
%   realistic single shell using superposition of basis functions (Nolte, 2003)
%   leadfield interpolation using a precomputed sourcemodel
%   boundary element method (BEM)
%
% The supported forward solutions for EEG are
%   infinite homogenous medium
%   infinite halfspace homogenous medium
%   single sphere
%   multiple concentric spheres (up to 4 spheres)
%   leadfield interpolation using a precomputed sourcemodel
%   boundary element method (BEM)
%
% See also FT_PREPARE_VOL_SENS, FT_HEADMODEL_ASA, FT_HEADMODEL_BEMCP,
% FT_HEADMODEL_CONCENTRICSPHERES, FT_HEADMODEL_DIPOLI, FT_HEADMODEL_HALFSPACE,
% FT_HEADMODEL_INFINITE, FT_HEADMODEL_LOCALSPHERES, FT_HEADMODEL_OPENMEEG,
% FT_HEADMODEL_SINGLESHELL, FT_HEADMODEL_SINGLESPHERE,
% FT_HEADMODEL_HALFSPACE

% Copyright (C) 2004-2016, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if iscell(sens) && iscell(headmodel) && numel(sens)==numel(headmodel)
  % this represents combined EEG, ECoG and/or MEG
  % use recursion to compute all leadfields
  lf = cell(1, numel(sens));
  for i=1:length(sens)
    lf{i} = ft_compute_leadfield(dippos, sens{i}, headmodel{i}, varargin{:});
  end
  lf = cat(1, lf{:});
  return;
end

% get the optional input arguments
reducerank      = ft_getopt(varargin, 'reducerank'); % default is handled below
backproject     = ft_getopt(varargin, 'backproject', 'yes');
normalize       = ft_getopt(varargin, 'normalize' , 'no');
normalizeparam  = ft_getopt(varargin, 'normalizeparam', 0.5);
weight          = ft_getopt(varargin, 'weight');
chanunit        = ft_getopt(varargin, 'chanunit');   % this is something like V, T, or T/m
dipoleunit      = ft_getopt(varargin, 'dipoleunit'); % this is something like nA*m

if any(strcmp(varargin(1:2:end), 'unit'))
  ft_error('the ''unit'' option is not supported any more, please use ''chanunit''');
end
if any(strcmp(varargin(1:2:end), 'units'))
  ft_error('the ''units'' option is not supported any more, please use ''chanunit''');
end

if ~isstruct(sens) && size(sens, 2)==3
  % definition of electrode positions only, restructure it
  sens = struct('elecpos', sens);
end

% ft_prepare_vol_sens should be called prior to ft_compute_leadfield
% to ensure that the sens and headmodel are up to date, since the backward
% compatibility check should not be performed for each dipole location
% sens       = ft_datatype_sens(sens);
% headmodel  = ft_datatype_headmodel(headmodel);

% determine whether it is EEG or MEG
iseeg = ft_senstype(sens, 'eeg');
ismeg = ft_senstype(sens, 'meg');

% determine the default for this option
if isempty(reducerank)
  if iseeg
    reducerank = 'no';    % for EEG
  elseif ismeg && ft_headmodeltype(headmodel, 'infinite')
    reducerank = 'no';    % for MEG with a magnetic dipole, e.g. a HPI coil
  else
    reducerank = 'yes';   % for MEG with a current dipole in a volume conductor
  end
end

% multiple dipoles can be represented either as a 1x(N*3) vector or as a
% as a Nx3 matrix, i.e. [x1 y1 z1 x2 y2 z2] or [x1 y1 z1; x2 y2 z2]
Ndipoles = numel(dippos)/3;
if all(size(dippos)==[1 3*Ndipoles])
  dippos = reshape(dippos, 3, Ndipoles)';
end

if isfield(headmodel, 'unit') && isfield(sens, 'unit') && ~strcmp(headmodel.unit, sens.unit)
  ft_error('inconsistency in the units of the volume conductor and the sensor array');
end

if ismeg && iseeg
  % this is something that could be implemented relatively easily
  ft_error('simultaneous EEG and MEG not supported');

elseif ~ismeg && ~iseeg
  ft_error('the input does not look like EEG, nor like MEG');

elseif ismeg
  switch ft_headmodeltype(headmodel)

    case 'singlesphere'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % MEG single-sphere volume conductor model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      coilpos = sens.coilpos; % position of each coil
      coilori = sens.coilori; % orientation of each coil

      if isfield(headmodel, 'o')
        % shift dipole and magnetometers to origin of sphere
        dippos  = dippos  - repmat(headmodel.o, Ndipoles, 1);
        coilpos = coilpos - repmat(headmodel.o, size(coilpos, 1), 1);
      end

      if Ndipoles>1
        % loop over multiple dipoles
        lf = zeros(size(coilpos, 1), 3*Ndipoles);
        for i=1:Ndipoles
          lf(:, (3*i-2):(3*i)) = meg_leadfield1(dippos(i, :), coilpos, coilori);
        end
      else
        % only single dipole
        lf = meg_leadfield1(dippos, coilpos, coilori);
      end

      if isfield(sens, 'tra')
        % this appears to be the modern complex gradiometer definition
        % construct the channels from a linear combination of all magnetometers
        lf = sens.tra * lf;
      end

    case 'localspheres'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % MEG multiple overlapping sphere volume conductor model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      ncoils = length(sens.coilpos);

      if size(headmodel.r, 1)~=ncoils
        ft_error('number of spheres is not equal to the number of coils')
      end

      if size(headmodel.o, 1)~=ncoils
        ft_error('number of spheres is not equal to the number of coils');
      end

      lf = zeros(ncoils, 3*Ndipoles);
      for coil=1:ncoils
        for dip=1:Ndipoles
          % shift dipole and magnetometer coil to origin of sphere
          tmppos  = dippos(dip, :) - headmodel.o(coil, :);
          coilpos = sens.coilpos(coil, :) - headmodel.o(coil, :);
          tmp = meg_leadfield1(tmppos, coilpos, sens.coilori(coil, :));
          lf(coil, (3*dip-2):(3*dip)) = tmp;
        end
      end

      if isfield(sens, 'tra')
        % this appears to be the modern complex gradiometer definition
        % construct the channels from a linear combination of all magnetometers
        lf = sens.tra * lf;
      end

    case 'neuromag'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % use external Neuromag toolbox for forward computation
      % this requires that "megmodel" is initialized, which is done in PREPARE_VOL_SENS
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % compute the forward model for all channels
      % tmp1 = ones(1, Ndipoles);
      % tmp2 = 0.01*dippos'; %convert to cm
      % lf = megfield([tmp2 tmp2 tmp2], [[1 0 0]'*tmp1 [0 1 0]'*tmp1 [0 0 1]'*tmp1]);
      for dip=1:Ndipoles
        R = 0.01*dippos(dip, :)'; % convert from cm to m
        Qx = [1 0 0];
        Qy = [0 1 0];
        Qz = [0 0 1];
        lf(:, (3*(dip-1)+1)) = megfield(R, Qx);
        lf(:, (3*(dip-1)+2)) = megfield(R, Qy);
        lf(:, (3*(dip-1)+3)) = megfield(R, Qz);
      end
      % select only those channels from the forward model that are part of the gradiometer definition
      lf = lf(headmodel.chansel, :);

    case 'singleshell'
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % use code from Guido Nolte for the forward computation
      % this requires that "meg_ini" is initialized, which is done in PREPARE_VOL_SENS
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % the dipole position and orientation should be combined in a single matrix
      % furthermore, here I want to compute the leadfield for each of the
      % orthogonal x/y/z directions
      dippar = zeros(Ndipoles*3, 6);
      for i=1:Ndipoles
        dippar((i-1)*3+1, :) = [headmodel.forwpar.scale*dippos(i, :) 1 0 0]; % single dipole with unit strength, x-orientation
        dippar((i-1)*3+2, :) = [headmodel.forwpar.scale*dippos(i, :) 0 1 0]; % single dipole with unit strength, y-orientation
        dippar((i-1)*3+3, :) = [headmodel.forwpar.scale*dippos(i, :) 0 0 1]; % single dipole with unit strength, z-orientation
      end
      % compute the leadfield for each individual coil
      lf = meg_forward(dippar, headmodel.forwpar);
      % the leadfield is computed for cm units, convert it to the desired units
      lf = lf*headmodel.forwpar.scale^2;
      if isfield(sens, 'tra')
        % compute the leadfield for each gradiometer (linear combination of coils)
        lf = sens.tra * lf;
      end

    case 'openmeeg'
        % OpenMEEG lead field already computed in ft_prepare_leadfield;
        % load here so any post-processing options (e.g. normalization) may
        % be applied
        lf = ft_getopt(varargin, 'lf');

    case {'infinite_magneticdipole', 'infinite'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % magnetic dipole instead of electric (current) dipole in an infinite vacuum
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      coilpos = sens.coilpos; % position of each coil
      coilori = sens.coilori; % orientation of each coil

      if Ndipoles>1
        % loop over multiple dipoles
        lf = zeros(size(coilpos, 1), 3*Ndipoles);
        for i=1:Ndipoles
          lf(:, (3*i-2):(3*i)) = magnetic_dipole(dippos(i, :), coilpos, coilori);
        end
      else
        % only single dipole
        lf = magnetic_dipole(dippos, coilpos, coilori);
      end

      if isfield(sens, 'tra')
        % construct the channels from a linear combination of all magnetometer coils
        lf = sens.tra * lf;
      end

    case {'infinite_currentdipole'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % current dipole in an infinite homogenous conducting medium
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      coilpos = sens.coilpos; % position of each coil
      coilori = sens.coilori; % orientation of each coil

      if Ndipoles>1
        % loop over multiple dipoles
        lf = zeros(size(coilpos, 1), 3*Ndipoles);
        for i=1:Ndipoles
          lf(:, (3*i-2):(3*i)) = current_dipole(dippos(i, :), coilpos, coilori);
        end
      else
        % only single dipole
        lf = current_dipole(dippos, coilpos, coilori);
      end

      if isfield(sens, 'tra')
        % construct the channels from a linear combination of all magnetometer coils
        lf = sens.tra * lf;
      end

    otherwise
      ft_error('unsupported volume conductor model for MEG');
  end % switch type for MEG

elseif iseeg
  switch ft_headmodeltype(headmodel)

    case 'multisphere'
      % Based on the approximation of the potential due to a single dipole in
      % a multishell sphere by three dipoles in a homogeneous sphere, code
      % contributed by Punita Christopher. Note that this one should not get
      % confused with the MEG localspheres model.

      Nelec    = size(sens.elecpos, 1);
      Nspheres = length(headmodel.r);

      % the center of the spherical volume conduction model does not have
      % to be in the origin, therefore shift the spheres, the electrodes
      % and the dipole
      if isfield(headmodel, 'o')
        center = headmodel.o;
      else
        center = [0 0 0];
      end

      % sort the spheres from the smallest to the largest
      % furthermore, the radius should be one (?)
      [radii, indx] = sort(headmodel.r/max(headmodel.r));
      sigma = headmodel.cond(indx);
      r = (sens.elecpos-repmat(center, Nelec, 1))./max(headmodel.r);
      dippos = dippos./max(headmodel.r);

      if Ndipoles>1
        % loop over multiple dipoles
        lf = zeros(Nelec, 3*Ndipoles);
        for i=1:Ndipoles
          rq = dippos(i, :) - center;
          % compute the potential for each dipole ortientation
          % it would be much more efficient to change the punita function
          q1 = [1 0 0]; lf(:, (3*i-2)) = multisphere(Nspheres, radii, sigma, r, rq, q1);
          q1 = [0 1 0]; lf(:, (3*i-1)) = multisphere(Nspheres, radii, sigma, r, rq, q1);
          q1 = [0 0 1]; lf(:, (3*i )) = multisphere(Nspheres, radii, sigma, r, rq, q1);
        end
      else
        % only single dipole
        lf = zeros(Nelec, 3);
        rq = dippos - center;
        % compute the potential for each dipole ortientation
        % it would be much more efficient to change the punita function
        q1 = [1 0 0] ; lf(:, 1) = multisphere(Nspheres, radii, sigma, r, rq, q1);
        q1 = [0 1 0] ; lf(:, 2) = multisphere(Nspheres, radii, sigma, r, rq, q1);
        q1 = [0 0 1] ; lf(:, 3) = multisphere(Nspheres, radii, sigma, r, rq, q1);
      end

    case {'singlesphere', 'concentricspheres'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % EEG spherical volume conductor model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      % FIXME, this is not consistent between spherical and BEM
      % sort the spheres from the smallest to the largest
      [headmodel.r, indx] = sort(headmodel.r);
      headmodel.cond = headmodel.cond(indx);

      Nspheres = length(headmodel.cond);
      if length(headmodel.r)~=Nspheres
        ft_error('the number of spheres in the volume conductor model is ambiguous');
      end

      if isfield(headmodel, 'o')
        % shift the origin of the spheres, electrodes and dipole
        sens.elecpos = sens.elecpos - repmat(headmodel.o, size(sens.elecpos, 1), 1);
        dippos = dippos - repmat(headmodel.o, Ndipoles, 1);
      end

      switch Nspheres
        case 1
          funnam = 'eeg_leadfield1';
        case 2
          headmodel.r = [headmodel.r(1) headmodel.r(2) headmodel.r(2) headmodel.r(2)];
          headmodel.cond = [headmodel.cond(1) headmodel.cond(2) headmodel.cond(2) headmodel.cond(2)];
          funnam = 'eeg_leadfield4';
        case 3
          headmodel.r = [headmodel.r(1) headmodel.r(2) headmodel.r(3) headmodel.r(3)];
          headmodel.cond = [headmodel.cond(1) headmodel.cond(2) headmodel.cond(3) headmodel.cond(3)];
          funnam = 'eeg_leadfield4';
        case 4
          headmodel.r = [headmodel.r(1) headmodel.r(2) headmodel.r(3) headmodel.r(4)];
          headmodel.cond = [headmodel.cond(1) headmodel.cond(2) headmodel.cond(3) headmodel.cond(4)];
          funnam = 'eeg_leadfield4';
        otherwise
          ft_error('more than 4 concentric spheres are not supported')
      end

      lf = zeros(size(sens.elecpos, 1), 3*Ndipoles);
      for i=1:Ndipoles
        lf(:, (3*i-2):(3*i)) = feval(funnam, dippos(i, :), sens.elecpos, headmodel);
      end

    case {'bem', 'dipoli', 'asa', 'bemcp'}
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % EEG boundary element method volume conductor model
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      lf = eeg_leadfieldb(dippos, sens.elecpos, headmodel);

    case 'openmeeg'
        % OpenMEEG lead field already computed in ft_prepare_leadfield;
        % load here so any post-processing options (e.g. normalization) may
        % be applied
        lf = ft_getopt(varargin, 'lf');

    case {'infinite_currentdipole' 'infinite'}
      lf = eeg_infinite_dipole(dippos, sens.elecpos, headmodel);

    case 'halfspace'
      lf = eeg_halfspace_dipole(dippos, sens.elecpos, headmodel);

    case 'infinite_monopole'
      lf = eeg_infinite_monopole(dippos, sens.elecpos, headmodel);

    case 'halfspace_monopole'
      lf = eeg_halfspace_monopole(dippos, sens.elecpos, headmodel);

    case 'slab_monopole'
      lf = eeg_slab_monopole(dippos, sens.elecpos, headmodel);

    case 'simbio'
      ft_hastoolbox('simbio', 1);
      % note that the electrode information is contained in the headmodel (thanks to ft_prepare_vol_sens)
      lf = leadfield_simbio(dippos, headmodel);

    case 'metufem'
      p3 = zeros(Ndipoles * 3, 6);
      for i = 1:Ndipoles
        p3((3*i - 2) : (3 * i), 1:3) = [dippos(i, :); dippos(i, :); dippos(i, :)];
        p3((3*i - 2) : (3 * i), 4:6) = [1 0 0; 0 1 0; 0 0 1];
      end
      lf = metufem('pot', p3', 'interp');

    case 'metubem'
      session = headmodel.session;
      p3 = zeros(Ndipoles * 3, 6);
      for i = 1:Ndipoles
        p3((3*i - 2) : (3 * i), 1:3) = [dippos(i, :); dippos(i, :); dippos(i, :)];
        p3((3*i - 2) : (3 * i), 4:6) = [1 0 0; 0 1 0; 0 0 1];
      end
      [lf, session] = bem_solve_lfm_eeg(session, p3);

    case 'fns'
      % note that the electrode information is contained in the headmodel
      % tolerance = 1e-8;
      lf = leadfield_fns(dippos, headmodel);

    case 'interpolate'
      % note that the electrode information is contained in the headmodel
      lf = leadfield_interpolate(dippos, headmodel);
      % the leadfield is already correctly referenced, i.e. it represents the
      % channel values rather than the electrode values. Prevent that the
      % referencing is done once more.
      sens.tra = speye(length(headmodel.filename));

    otherwise
      ft_error('unsupported volume conductor model for EEG');

  end % switch type for EEG

  % the forward model potential is computed on the electrodes relative to
  % an unknown reference, not on the channels. Therefore the data has to be
  % explicitly referenced here.
  if isfield(sens, 'tra')
    % apply the correct montage to the leadfield
    lf = sens.tra*lf;
  else
    % compute average reference for EEG leadfield
    for i=1:size(lf,2)
      lf(:,i) = lf(:,i) - mean(lf(:,i));
    end
  end

end % iseeg or ismeg

% optionally apply leadfield rank reduction
switch reducerank
  case 'yes'
    reducerank = 2;
  case 'no'
    reducerank = 3;
  otherwise
    % assume that it is specified as a number, keep it like this
end

if reducerank<size(lf,2)
  % decompose the leadfield
  for ii=1:Ndipoles
    tmplfd=lf(:, (3*ii-2):(3*ii));
    [u, s, v] = svd(tmplfd);
    r = diag(s);
    s(:) = 0;
    for j=1:reducerank
      s(j, j) = r(j);
    end

    if istrue(backproject)
      % recompose the leadfield with reduced rank
      lf(:, (3*ii-2):(3*ii)) = u * s * v';
    else
      % if not backprojected, the new leadfield has a different dimension
      if ii==1
        newlf    = zeros(size(lf,1), Ndipoles*reducerank);
        origrank = size(lf,2)./Ndipoles;
      end
      newlf(:, reducerank*(ii-1) + (1:reducerank)) = lf(:, origrank*(ii-1) + (1:origrank))*v(:,1:reducerank);
    end
  end

  if ~istrue(backproject)
    lf = newlf;
  end
  clear newlf;
end

% optionally apply leadfield normalization
switch normalize
  case 'yes'
    for ii=1:Ndipoles
      tmplf = lf(:, (3*ii-2):(3*ii));
      if normalizeparam==0.5
        % normalize the leadfield by the Frobenius norm of the matrix
        % this is the same as below in case normalizeparam is 0.5
        nrm = norm(tmplf, 'fro');
      else
        % normalize the leadfield by sum of squares of the elements of the leadfield matrix to the power "normalizeparam"
        % this is the same as the Frobenius norm if normalizeparam is 0.5
        nrm = sum(tmplf(:).^2)^normalizeparam;
      end
      if nrm>0
        tmplf = tmplf ./ nrm;
      end
      lf(:, (3*ii-2):(3*ii)) = tmplf;
    end
  case 'column'
    % normalize each column of the leadfield by its norm
    for ii=1:Ndipoles
      tmplf = lf(:, (3*ii-2):(3*ii));
      for j=1:size(tmplf, 2)
        nrm = sum(tmplf(:, j).^2)^normalizeparam;
        tmplf(:, j) = tmplf(:, j)./nrm;
      end
      lf(:, (3*ii-2):(3*ii)) = tmplf;
    end
end

% optionally apply a weight to the leadfield for each dipole location
if ~isempty(weight)
  for i=1:Ndipoles
    lf(:, 3*(i-1)+1) = lf(:, 3*(i-1)+1) * weight(i); % the leadfield for the x-direction
    lf(:, 3*(i-1)+2) = lf(:, 3*(i-2)+1) * weight(i); % the leadfield for the y-direction
    lf(:, 3*(i-1)+3) = lf(:, 3*(i-3)+1) * weight(i); % the leadfield for the z-direction
  end
end

if ~isempty(chanunit) || ~isempty(dipoleunit)
  assert(strcmp(headmodel.unit,  'm'), 'unit conversion only possible for SI input units');
  assert(strcmp(sens.unit, 'm'), 'unit conversion only possible for SI input units');
end

if ~isempty(chanunit)
  assert(all(strcmp(sens.chanunit, 'V') | strcmp(sens.chanunit, 'V/m') | strcmp(sens.chanunit, 'T') | strcmp(sens.chanunit, 'T/m')), 'unit conversion only possible for SI input units');
  % compute conversion factor and multiply each row of the matrix
  scale = cellfun(@ft_scalingfactor, sens.chanunit(:), chanunit(:));
  lf = bsxfun(@times, lf, scale(:));
  % prior to this conversion, the units might be  (T/m)/(A*m) for planar gradients or   (V/m)/(A*m) for bipolar EEG
  % after this conversion, the units will be     (T/cm)/(A*m)                      or (uV/mm)/(A*m)
end

if ~isempty(dipoleunit)
  scale = ft_scalingfactor('A*m', dipoleunit); % compue the scaling factor from A*m to the desired dipoleunit
  lf    = lf/scale;                         % the leadfield is expressed in chanunit per dipoleunit, i.e. chanunit/dipoleunit
end

author: spm


text:
function classifier = lapsvm(options,data)
% {lapsvm} trains a Laplacian SVM classifier (dual solved with Libsvm).
%
%      classifier = lapsvm(options,data)
%
%      options: a structure with the following fields
%               options.gamma_A: regularization parameter (ambient space)
%               options.gamma_I: regularization parameter (intrinsic norm)
%
%               [optional fields]
%               options.UseBias: {0,1} i.e. use or not a bias (fx=Kalpha+b)
%      data: a structure with the following fields
%            data.X: a M-by-D matrix of M D-dimensional training examples
%            data.K: a M-by-M kernel Gram matrix of M training examples
%            data.Y: a M-by-1 label vector in {-1,0,+1}, where 0=unlabeled
%            data.L: a M-by-M matrix of the graph Laplacian
%
%      classifier: structure of the trained classifier (see the
%                  'saveclassfier' function)
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it
%         * based on the code of Vikas Sindhwani, vikas.sindhwani@gmail.com

tic
if ~isfield(options,'UseBias'),           options.UseBias=1; end

C=1;
parameters =    [4 ... % kernel_type (4=Gram Matrix)
                1 ... % deg
                1 ... % gamma
                0 ... % coef0
                C ... % C
                40.00 ... % cache
                0.001 ... % eps
                0 ... % svm_type
                0.5 ... % nu
                0.1 ... % p
                1 ... % shrinking
                ];

I=eye(size(data.K,1));
lab=(data.Y~=0);
l=nnz(lab);

if isempty(data.L) || options.gamma_I==0
    G=I/(2*options.gamma_A);
else
    G=(2*options.gamma_A*I + 2*options.gamma_I*data.L*data.K)\I;
end

Gram=data.K(lab,:)*G(:,lab);
Ylab=data.Y(lab);

[beta, svs, b, nsv, nlab] = mexGramSVMTrain(Gram', Ylab', parameters);

if nlab(1)==-1
    beta=-beta;
else
    b=-b;
end

betaz=zeros(l,1);
betaz(svs)=beta';
alpha=G(:,lab)*betaz;
sec=toc;

classifier = saveclassifier('lapsvm',1:length(data.Y),alpha, ...
                            data.X,b*options.UseBias,options,sec);

author: viggin

repo: domain-adaptation-toolbox


text:
function Population = subNSGAII(Population,N)
% The environmental selection of NSGA-II

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    %% Non-dominated sorting
    CC = min([N,length(Population)]);
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,CC);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:CC-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function f_x = ParFor5(in1)
%PARFOR5
%    F_X = PARFOR5(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    31-Mar-2020 16:52:33

r = in1(:,1);
s = in1(:,13);
u = in1(:,19);
ut = in1(:,20);
ux = in1(:,21);
uxx = in1(:,22);
uy = in1(:,23);
uyy = in1(:,24);
z = in1(:,7);
t2 = u.^2;
f_x = r.*-4.845016406907234+s.*6.96612441301113+ut.*5.457984690783633e-1-ux.*1.530748504418966e-2+uxx.*2.135644249305187-uy.*3.570878181824355e-2+uyy.*5.02823932467436-z.*1.323504040227272e2+r.*u.*6.191151455437648-s.*u.*9.105578580027213-t2.*u.*1.077295023277402e2+u.*ut.*5.65024831920482e-2+u.*ux.*1.830355084143775e-2-u.*uxx.*4.132096419489244+u.*uy.*4.442896435307375e-2-u.*uyy.*6.867807677976089+u.*z.*1.83605552916415e2+4.028421339788474e1;

author: dynamicslab

repo: SINDy-PI


text:
function  [im_out, Par] = TWSC_Sigma_RW(Par)
im_out     =   Par.nim;
% parameters for noisy image
[h,  w, ch] = size(im_out);
Par.h = h;
Par.w = w;
Par.ch = ch;
Par = SearchNeighborIndex( Par );
% original noisy image to patches
NY = Image2Patch( Par.nim, Par );
Par.Sigma = sqrt(mean(Par.nSig.^2));
for ite  =  1 : Par.Outerloop
    % iterative regularization
    im_out = im_out + Par.delta * (Par.nim - im_out);
    % image to patches
    Y = Image2Patch( im_out, Par );
    % estimate local noise variance, par.lambdals is put here since the MAP and Bayesian rules
    if Par.lambda1 ~= 0
        if ite == 1
            SigmaRow = zeros(size(NY));
        else
            SigmaRow = (NY-Y).^2;
        end
    end
    SigmaCol = Par.lambda2 * sqrt(abs(repmat(Par.Sigma^2, 1, size(Y,2)) - mean((NY - Y).^2))); % Estimated Local Noise Level
    % estimation of noise variance
    if mod(ite-1, Par.Innerloop)==0
        Par.nlsp = max(Par.nlspgap, Par.nlsp - Par.nlspgap);
        % searching non-local patches
        blk_arr = Block_Matching_RW( Y, Par );
    end
    % Weighted Sparse Coding
    Y_hat = zeros(Par.ps2ch, Par.maxrc, 'double');
    W_hat = zeros(Par.ps2ch, Par.maxrc, 'double');
    for i = 1:Par.lenrc
        index = blk_arr(:, i);
        nlY = Y( : , index );
        DC = mean(nlY, 2);
        nDCnlY = bsxfun(@minus, nlY, DC);
        % Compute W2
        W2 = 1 ./ (SigmaCol(index) + eps);
        % update D
        [D, S, ~] = svd( full(nDCnlY), 'econ' );
        % update S
        S = sqrt(max( diag(S).^2 - length(index) * SigmaCol(index(1))^2, 0 ));
        if Par.lambda1 == 0
            % update weight for sparse coding
            Wsc = bsxfun( @rdivide, SigmaCol(index).^2, S + eps );
            % update C by soft thresholding
            B = D' * nDCnlY;
            C = sign(B) .* max( abs(B) - Wsc, 0 );
            % update Y
            nDCnlYhat = D * C;
        else
            W1 = exp( - Par.lambda1*mean(SigmaRow(:, index), 2)); % max?
            S = diag(S);
            % min |Z|_1 + |W1(Y-DSC)W2|_F,2  s.t.  C=Z
            C = TWSC_ADMM( nDCnlY, D, S, W1, W2, Par );
            % update Y
            nDCnlYhat = D * S * C;
        end
        % add back DC components
        nlYhat = bsxfun(@plus, nDCnlYhat, DC);
        % aggregation
        Y_hat(:, index) = Y_hat(:, index) + bsxfun(@times, nlYhat, W2);
        W_hat(:, index) = W_hat(:, index) + repmat(W2, [Par.ps2ch, 1]);
    end
    % Reconstruction
    im_out = PGs2Image(Y_hat, W_hat, Par);
end
return;

author: csjunxu


text:
% Script demonstrating the latent Dirichlet allocation model on synthetic
% data.

% Script parameters.
K      = 3;   % The number of topics.
D      = 25;  % The number of documents (images).
L      = 16;  % The length of each document.
height = 4;   % The height of the images (in pixels).
width  = 4;   % The width of the images (in pixels).
nr     = 4;   % Number of rows of documents to display.
nc     = 6;   % Number of columns of documents to display.
ns     = 3;   % Number of samples to display.

% The size of the vocabulary (W) is equal to the number of pixels.
W = height*width;
L = repmat(L,1,D);

% Set a uniform prior distribution over the topics.
nu = ones(1,K)/K;

% Set a uniform prior distribution over the words.
eta = ones(1,W)/W;

% Randomly generate the topics and the documents.
[topics w] = createsyntheticdata(nu,eta,K,W,L);
wc = zeros(D,W);
for d = 1:D
  wc(d,:) = hist(w{d},1:W)';
end

% Display the topics.
figure(1);
set(gcf,'MenuBar','none');
set(gcf,'NumberTitle','off');
set(gcf,'Name','True topics');
set(gcf,'Color','white');
pos = get(gcf,'Position');
set(gcf,'Position',[pos(1) pos(2) 228 65]);
colormap(gray);
maxval = max(topics(:));
for k = 1:K
  subplot(1,K,k);
  imagesc(reshape(topics(:,k),height,width),[0 maxval]);
  axis off
end

% Show some of the documents that were generated.
figure(2);
set(gcf,'MenuBar','none');
set(gcf,'NumberTitle','off');
set(gcf,'Name','Data');
set(gcf,'Color','white');
pos = get(gcf,'Position');
set(gcf,'Position',[pos(1) pos(2) 222 355]);
colormap(gray);
maxval = max(wc(:));
d = 0;
for c = 1:nc
  for r = 1:nr
    d = d + 1;
    subplot(nc,nr,d);
    imagesc(reshape(wc(d,:),height,width),[0 maxval]);
    axis off
  end
end
drawnow

% Compute the mean field estimate.
tic;
fprintf('Finding solution to the mean field lower bound.\n');
[xi0 gamma0 Phi0]  = mfldainit(W,K,w);
[LnZ xi gamma Phi] = mflda(nu,eta,w,xi0,gamma0,Phi0);
fprintf('Optimization took %i seconds.\n', round(toc));

% Show the estimated topics.
figure(3);
set(gcf,'MenuBar','none');
set(gcf,'NumberTitle','off');
set(gcf,'Name','Topics (mean field)');
set(gcf,'Color','white');
pos = get(gcf,'Position');
set(gcf,'Position',[pos(1) pos(2) 197 226]);
colormap(gray);

% Repeat for every sample to display.
for i = 1:ns
  beta = zeros(W,K);
  for k = 1:K
    beta(:,k) = dirichletrnd(xi(:,k));
  end
  maxval = max(beta(:));

  for k = 1:K
    subplot(ns,K,K*(i-1)+k);
    imagesc(reshape(beta(:,k),height,width),[0 maxval]);
    axis off
    title(' ');
  end
  subplot(ns,K,K*(i-1)+1);
  title(sprintf('Sample #%i',i));
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function metoffice_plot_fa(Q)

[D,N] = size(Q.X);
M = size(Q.W,2);

E2 = reshape(Q.CovX, [D^2, N]);
E = sqrt(E2(1:(D+1):end,:));

hax = tserrorplot(Q.X',2*E');
ylim_centered(hax);
ylim_scale(hax, 1.2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
classdef (InferiorClasses = {?SO3FunHarmonic,?SO3FunBingham,?SO3FunCBF, ...
    ?SO3FunComposition,?SO3FunHandle,?SO3FunHomochoric,?SO3FunRBF, ...
    ?SO3VectorFieldHandle}) SO3VectorFieldHarmonic < SO3VectorField
% a class representing a function on the rotation group

properties
  SO3F
end

properties(Dependent = true)
  x
  y
  z
  bandwidth
  SLeft
  SRight
  antipodal
  isReal
end

methods

  function SO3VF = SO3VectorFieldHarmonic(SO3F, varargin)
    % initialize a rotational vector field

    if nargin == 0, return; end
    if isa(SO3F,'SO3FunHarmonic')
      SO3VF.SO3F = SO3F(:);
      return
    end
    if isa(SO3F,'SO3VectorFieldHarmonic')
      SO3VF.SO3F = SO3F.SO3F;
      return
    end
    if isa(SO3F,'SO3VectorField')
      SO3VF = SO3VectorFieldHarmonic.quadrature(SO3F,varargin{:});
      return
    end
    error('Input should be of type SO3FunHarmonic or SO3VectorField.')

  end

  function bw = get.bandwidth(SO3VF), bw = SO3VF.SO3F.bandwidth; end
  function SO3VF = set.bandwidth(SO3VF,bw), SO3VF.SO3F.bandwidth = bw; end

  function x = get.x(SO3VF), x = SO3VF.SO3F(1); end
  function y = get.y(SO3VF), y = SO3VF.SO3F(2); end
  function z = get.z(SO3VF), z = SO3VF.SO3F(3); end
  function SO3VF = set.x(SO3VF, x), SO3VF.SO3F(1) = x; end
  function SO3VF = set.y(SO3VF, y), SO3VF.SO3F(2) = y; end
  function SO3VF = set.z(SO3VF, z), SO3VF.SO3F(3) = z; end

  function SLeft = get.SLeft(SO3VF), SLeft = SO3VF.SO3F.SLeft; end
  function SRight = get.SRight(SO3VF), SRight = SO3VF.SO3F.SRight; end
  function SO3VF = set.SLeft(SO3VF,SLeft), SO3VF.SO3F.SLeft = SLeft; end
  function SO3VF = set.SRight(SO3VF,SRight), SO3VF.SO3F.SRight = SRight; end

  function a = get.antipodal(SO3VF), a = SO3VF.SO3F.antipodal; end
  function SO3VF = set.antipodal(SO3VF,a), SO3VF.SO3F.antipodal = a; end
  function r = get.isReal(SO3VF), r = SO3VF.SO3F.isReal; end
  function SO3VF = set.isReal(SO3VF,r), SO3VF.SO3F.isReal = r; end

end

methods(Static = true)
  SO3VF = quadrature(f, varargin)
  SO3VF = approximation(f, varargin)
end

end

author: mtex-toolbox

repo: mtex

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


text:
%%
% Copyright (c) 2019-present, Mahmoud Afifi
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
%
% Please, cite the following paper if you use this code:
%
% Mahmoud Afifi and Michael S. Brown. What else can fool deep learning?
% Addressing color constancy errors on deep neural network performance.
% ICCV, 2019
%
% Email: mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%%

classdef PCAFeature
    properties
        weights
        bias
    end
    methods
        function feature = encode(obj,hist)
		% Generates a compacted feature of a given histogram tensor.
            feature = (reshape(hist,1,[]) - obj.bias') *obj.weights;
        end
    end
end

author: mahmoudnafifi


text:
function determ = wilk21_determinant ( n )

%*****************************************************************************80
%
%% WILK21_DETERMINANT computes the determinant of the WILK21 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  d = zeros ( n, 1 );

  for i = 1 : n
    d(i) = round ( abs ( i - ( n + 1 ) / 2 ) );
  end

  determ_nm1 = d(n);

  if ( n == 1 )
    determ = determ_nm1;
    return
  end

  determ_nm2 = determ_nm1;
  determ_nm1 = d(n-1) * d(n) - 1.0;

  if ( n == 2 )
    determ = determ_nm1;
    return
  end

  for i = n - 2 : -1 : 1

    determ = d(i) * determ_nm1 - determ_nm2;

    determ_nm2 = determ_nm1;
    determ_nm1 = determ;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% brute_force_tune
% Code to test the performance of various tuning parameters
% Works sorta like RANSAC I guess?
% Adam Werries 2016, see Apache 2.0 license.

k_max = 50;
% Specify ranges
accel_bias_PSD = logspace(-10,-4,100);
gyro_bias_PSD = logspace(-10,-4,100);
% Repeat arrays
accel_bias_PSD = repmat(accel_bias_PSD, [1 k_max]);
gyro_bias_PSD = repmat(gyro_bias_PSD, [1 k_max]);
% Generate random selections of each vector
num_items = length(accel_bias_PSD);
accel_bias_i = randperm(num_items);
gyro_bias_i = randperm(num_items);

rms_error_filter = Inf*ones(1,num_items);
max_error_filter = Inf*ones(1,num_items);
parfor i = 1:num_items
    fprintf('Iteration: %d, ABias: %08.5e, GBias: %08.5e\n', i, accel_bias_PSD(accel_bias_i(i)), gyro_bias_PSD(gyro_bias_i(i)));
    temp_conf = LC_KF_config;
    temp_conf.accel_bias_PSD = accel_bias_PSD(accel_bias_i(i));
    temp_conf.gyro_bias_PSD = gyro_bias_PSD(gyro_bias_i(i));
    [out_profile,out_IMU_bias_est,out_KF_SD] = Loosely_coupled_INS_GNSS(init_cond, filter_time, epoch, lla, gps, imu, temp_conf, est_IMU_bias);
    xyz = out_profile(:,2:4);
    if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
        llh = ecef2lla(xyz);
        [x,y] = deg2utm(llh(:,1),llh(:,2));
        x = x-min_x;
        y = y-min_y;
%         h = -llh(:,3);
%         distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2 + (ground_truth_full(:,3)-h).^2).^0.5;
        distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
        rms_error_filter(i) = rms(distance);
        max_error_filter(i) = max(distance);
    end
end

[minmax, i] = min(max_error_filter);
fprintf('\nBest max: %08.4f, rms is %08.4f\n', minmax, rms_error_filter(i));
fprintf('Best iteration for max: %d, ABias: %08.5e, GBias: %08.5e\n', i, accel_bias_PSD(accel_bias_i(i)), gyro_bias_PSD(gyro_bias_i(i)));
[minrms, i] = min(rms_error_filter);
fprintf('Best rms: %08.4f, max is %08.4f\n', minrms, max_error_filter(i));
fprintf('Best iteration for rms: %d, ABias: %08.5e, GBias: %08.5e\n', i, accel_bias_PSD(accel_bias_i(i)), gyro_bias_PSD(gyro_bias_i(i)));
[minrms, i] = min((rms_error_filter+max_error_filter)/2);
fprintf('Best average of RMS and max: %08.4f, rms is  %08.4f, max is %08.4f\n', minrms, rms_error_filter(i), max_error_filter(i));
fprintf('Best iteration for rms: %d, ABias: %08.5e, GBias: %08.5e\n', i, accel_bias_PSD(accel_bias_i(i)), gyro_bias_PSD(gyro_bias_i(i)));

author: awerries

repo: kalman-localization


text:
options.P = 3;
filters = morlet_filter_bank_2d_pyramid(options);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
classdef RegularizedSuperEllipseExponentPlotter < handle

    properties (Access = private)
        vademecum
        mesh
        errorF
        qAveraged
        qRegularized
        qComputer
    end

    methods (Access = public)

        function obj = RegularizedSuperEllipseExponentPlotter()
            obj.createQcomputer();
            obj.loadVademecum();
            obj.createMesh();
            obj.createErrorFunction();
            obj.computeNumericalAveragedValue();
            obj.computeQregularized();
            obj.plotQforSymmetricHoles();
            obj.plotSuperEllipseOptimalExponent();
        end

    end

    methods (Access = private)

        function createQcomputer(obj)
            s.m1 = [];
            s.m2 = [];
            obj.qComputer = SmoothingExponentComputerOptimal(s);
        end

        function loadVademecum(obj)
            obj.vademecum = VademecumReader();
        end

        function createMesh(obj)
            xi = obj.vademecum.xiV;
            rho = obj.vademecum.rhoV;
            obj.mesh = obj.obtainMesh(xi,rho);
        end

        function computeNumericalAveragedValue(obj)
            s.vademecum = obj.vademecum;
            p = PonderatedOptimalSuperEllipseComputer(s);
            p.compute();
            obj.qAveraged = p.qMean;
        end

        function error = computeError(obj,a,b,r,qS)
            q = obj.computeQ(a,b,r,qS);
            error = obj.errorF(q);
        end

        function q = computeQ(obj,a,b,r,q)
             obj.qComputer.setParamsValues(a,b,r,q)
             xi  = obj.vademecum.xiV;
             rho = obj.vademecum.rhoV;
             q = obj.qComputer.computeQ(xi,rho);
        end

        function [alpha,beta,rhoQmin,qSoptMin] = computeBestParameters(obj)
            s.errorComputer = @(a,b,c,d) obj.computeError(a,b,c,d);
            r = RegularizedExponentBestCoeffComputer(s);
            [alpha,beta,rhoQmin,qSoptMin] = r.compute();
        end

        function createErrorFunction(obj)
            s.type              = 'SIMPLE';
            s.mesh              = obj.mesh;
            s.backgroundMesh    = obj.mesh;
            s.globalConnec      = obj.mesh.connec;
            s.npnod = obj.mesh.nnodes;
            int = Integrator.create(s);
            int.computeLHS();
            p = 2;
            normF = @(x) ((int.computeL2norm(x)).^(p))^(1/p);
            obj.errorF = @(x) normF(abs(x - obj.qAveraged))/normF(obj.qAveraged);
        end

        function computeQregularized(obj)
            [alpha,beta,rhoQmin,qSoptMin] = obj.computeBestParameters();
            q = obj.computeQ(alpha,beta,rhoQmin,qSoptMin);
            obj.qRegularized = q;
        end

        function plotQforSymmetricHoles(obj)
            xi     = obj.vademecum.xiV;
            isSym  = abs(xi-pi/4) < 0.01;
            isSym  = isSym(2:end);
            rhoSym = obj.vademecum.rhoV(isSym);
            qAsym  = obj.qAveraged(isSym);
            qRsym  = obj.qRegularized(isSym);
            f = figure();
            hold on
            h{1} = plot(rhoSym,qAsym,'-+');
            h{2} = plot(rhoSym,qRsym,'-+');
            xlabel('$\rho$','Interpreter','Latex')
            ylabel('$q(\xi = \pi/4,\rho)$','Interpreter','Latex')
            legN = '$\textrm{Numerical smoothing exponent} \, q_N$';
            legA = '$\textrm{Analytical smoothing exponent} \, q_A$';
            legend({legN,legA},'Interpreter','Latex','Location','Best');
            outPutPath = '/home/alex/git-repos/MicroStructurePaper/';
            outputName = [outPutPath,'qMaxM1M2AnalyticalNumerical'];
            printer = plotPrinter(f,h);
            printer.print(outputName);
        end

        function plotSuperEllipseOptimalExponent(obj)
            s.title = 'Proposed ';
            s.fileName = '/home/alex/git-repos/MicroStructurePaper/ProposedOptimal';
            s.rhoV = rho;
            s.xiV = xi;
            s.value = qA;
            s.qMean = qA;
            p = SuperEllipseExponentPlotter(s);
            p.plot();
        end

        function m = obtainMesh(obj,x,y)
            s.coord = [x,y];
            s.connec = obj.obtainConnec(x,y);
            m = Mesh().create(s);
        end

        function connec = obtainConnec(obj,x,y)
            connec = delaunay(x,y);
            connec = obj.obtainQualityElements(connec,x,y);
        end

    end

    methods (Access = private, Static)

        function connec = obtainQualityElements(connec,x,y)
            s.coord = [x,y];
            s.connec = connec;
            m = Mesh().create(s);
            qua = m.computeElementQuality';
            isQ = qua > 0.02;
            connec = connec(isQ,:);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function mObject=createMobject
%   Create Mobjcet for morphing
%   mObject=createMobject

%   Designed and coded by Hideki Kawahara
%   25/February/2005
%   14/October/2005 Added creator information

mObject.date = datestr(now);
mObject.pwd = pwd;
mObject.waveform = [];
mObject.samplingFrequency = 44100; % default frequency
mObject.F0 = [];
mObject.vuv = [];
mObject.spectrogram = [];
mObject.aperiodicityIndex = [];
mObject.frameUpdateInterval = 1; % default frame is 1ms
mObject.anchorTimeLocation = [];
mObject.maximumFrequencyPoints = 9; % default max frequency anchor points
mObject.anchorFrequency = [];
mObject.F0extractionConditions = [];
mObject.SpectrumExtractionConditions = [];
mObject.creatorInformation = which('createMobject');

author: HidekiKawahara


text:
%The spherical K-means algorithm
%(C) 2007-2011 Nguyen Xuan Vinh
%Contact: vinh.nguyenx at gmail.com   or vinh.nguyen at monash.edu
%Reference:
%   [1] Xuan Vinh Nguyen: Gene Clustering on the Unit Hypersphere with the
%       Spherical K-Means Algorithm: Coping with Extremely Large Number of Local Optima. BIOCOMP 2008: 226-233
%Usage: Normalize the data set to have mean zero and unit norm

function b=normalize_norm_mean(a)
[n dim]=size(a);
for i=1:n
    a(i,:)=(a(i,:)-mean(a(i,:)))/norm(a(i,:)-mean(a(i,:)),2);
end
b=a;

author: Sable

repo: mcbench-benchmarks


text:
function x = cc_abscissas ( n )

%*****************************************************************************80
%
%% CC_ABSCISSAS computes the Clenshaw Curtis abscissas.
%
%  Discussion:
%
%    The interval is [ -1, 1 ].
%
%    The abscissas are the cosines of equally spaced angles between
%    180 and 0 degrees, including the endpoints.
%
%      X(I) = cos ( ( ORDER - I ) * PI / ( ORDER - 1 ) )
%
%    except for the basic case ORDER = 1, when
%
%      X(1) = 0.
%
%    If the value of ORDER is increased in a sensible way, then
%    the new set of abscissas will include the old ones.  One such
%    sequence would be ORDER(K) = 2*K+1 for K = 0, 1, 2, ...
%
%    When doing interpolation with Lagrange polynomials, the Clenshaw Curtis
%    abscissas can be a better choice than regularly spaced abscissas.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Charles Clenshaw, Alan Curtis,
%    A Method for Numerical Integration on an Automatic Computer,
%    Numerische Mathematik,
%    Volume 2, Number 1, December 1960, pages 197-205.
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%    Joerg Waldvogel,
%    Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules,
%    BIT Numerical Mathematics,
%    Volume 43, Number 1, 2003, pages 1-18.
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
  if ( n == 1 )
    x(1) = 0.0;
    return
  end

  for i = 1 : n
    theta(i) = ( n - i ) * pi ...
             / ( n - 1 );
  end

  x(1:n) = cos ( theta(1:n) );

  return
end


text:
function [oAllanDev] = calculateADEV(tau,sPeriod,readings)
%Overlapping Allan (ADEV) function that uses phase error or time error values
%input argument is readings already grouped as tau values. tau is the
%desired gate or tau time needed for calculating overlap.
%sPeriod is used to determine the overlap

N = numel(readings); %get the reading count
n = tau / sPeriod; %calculate averaging factor
n = floor(n); %make sure this is an integer
const = 1/(2*(N - (2*n))*tau^2); %calculate the const mult 1/(2*(N - (2*n))*tau^2)
%sum from i=1 to N-2n (Xi+2m - 2Xi+m + Xi)^2
sum = 0; %variable to store summation calculation

%loop for performing summation
for i = 1:(N-(2*n))
   sum = sum + (readings(i+(2*n)) - (2*readings(i+n)) + readings(i))^2; %previos sum + (yi+1 - yi)^2
end

oAllanDev = sqrt(const*sum); %square root for Allan Dev


text:
function DCM = spm_dcm_tfm_data(DCM)
% gets cross-spectral density data-features using a wavelet transform
% FORMAT DCM = spm_dcm_tfm_data(DCM)
% DCM    -  DCM structure
% requires
%
%    DCM.xY.Dfile        - name of data file
%    DCM.M.U             - channel subspace
%    DCM.options.trials  - trial to evaluate
%    DCM.options.Tdcm    - time limits
%    DCM.options.Fdcm    - frequency limits
%    DCM.options.D       - Down-sampling
%
% sets
%
%    DCM.xY.pst     - Peristimulus Time [ms] sampled
%    DCM.xY.dt      - sampling in seconds [s] (down-sampled)
%    DCM.xY.U       - channel subspace
%    DCM.xY.y       - cross spectral density over channels
%    DCM.xY.csd     - cross spectral density over channels
%    DCM.xY.erp     - event-related average over channels
%    DCM.xY.It      - Indices of time bins
%    DCM.xY.Ic      - Indices of good channels
%    DCM.xY.Hz      - Frequency bins
%    DCM.xY.code    - trial codes evaluated
%    DCM.xY.scale   - scalefactor applied to data
%    DCM.xY.Rft     - Wavelet number or ratio of frequency to time
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_tfm_data_nopad.m 6110 2014-07-21 09:36:13Z karl $

% Set defaults and Get D filename
%-------------------------------------------------------------------------

% channel indices (excluding bad channels)
%--------------------------------------------------------------------------
DCM       = spm_dcm_erp_data(DCM,0);
Ic        = DCM.xY.Ic;
It        = DCM.xY.It;
Nb        = length(It);

% ensure spatial modes have been computed (see spm_dcm_csd)
%-------------------------------------------------------------------------
if ~isfield(DCM.M,'U')
    DCM.M.U = spm_dcm_eeg_channelmodes(DCM.M.dipfit,DCM.options.Nmodes);
end
if size(DCM.M.U,1) ~= length(Ic);
    DCM.M.U = spm_dcm_eeg_channelmodes(DCM.M.dipfit,DCM.options.Nmodes);
end
Nm          = size(DCM.M.U,2);

% options
%--------------------------------------------------------------------------
try, trial = DCM.options.trials; catch, trial = D.nconditions; end
try, Rft   = DCM.options.Rft;    catch, Rft   = 8;             end

% get frequency range
%--------------------------------------------------------------------------
try
    Hz1     = DCM.options.Fdcm(1);          % lower frequency
    Hz2     = DCM.options.Fdcm(2);          % upper frequency
catch
    pst     = DCM.xY.pst(end) - DCM.xY.pst(1);
    Hz1     = max(ceil(2*1000/pst),4);
    if Hz1 < 8;
        Hz2 = 48;
    else
        Hz2 = 128;
    end
end

% Frequencies
%--------------------------------------------------------------------------
Hz  = fix(Hz1:Hz2);                         % Frequencies
Nf  = length(Hz);                           % number of frequencies
Ne  = length(trial);                        % number of ERPs

% Cross spectral density for each trial type
%==========================================================================
erp   = cell(1,Ne);
csd   = cell(1,Ne);
for e = 1:Ne;

    % evoked response
    %----------------------------------------------------------------------
    Nt    = DCM.xY.nt(e);
    Y     = zeros(Nb,Nm,Nt);
    P     = zeros(Nb,Nf,Nm,Nm);
    Q     = zeros(Nb,Nf,Nm,Nm);
    for k = 1:Nt
        Y(:,:,k) = full(DCM.xY.y{e}(:,:,k)*DCM.M.U);
    end

    % ERP
    %----------------------------------------------------------------------
    erp{e} = mean(Y,3);

    % induced response
    %----------------------------------------------------------------------
    for k = 1:Nt

        fprintf('\nevaluating condition %i (trial %i)',e,k)
        G     = spm_morlet(Y(:,:,k) - erp{e},Hz*DCM.xY.dt,Rft);
        for i = 1:Nm
            for j = 1:Nm
                P(:,:,i,j) = (G(:,:,i).*conj(G(:,:,j)));
            end
        end
        Q = Q + P;
    end

    % normalise induced responses in relation to evoked responses
    %--------------------------------------------------------------------------
    Vm    = mean(mean(squeeze(var(Y,[],3))));
    Vs    = mean(diag(squeeze(mean(squeeze(mean(Q))))));
    Q     = Vm*Q/Vs;

    % store
    %----------------------------------------------------------------------
    csd{e} = Q;

end

% place cross-spectral density in xY.y
%==========================================================================
try
    scale = DCM.xY.scale;
catch
    scale = 1;
end

DCM.xY.erp   = spm_unvec(spm_vec(erp)*(scale^1),erp);
DCM.xY.csd   = spm_unvec(spm_vec(csd)*(scale^2),csd);
DCM.xY.y     = DCM.xY.csd;
DCM.xY.U     = DCM.M.U;
DCM.xY.scale = scale;
DCM.xY.Hz    = Hz;
DCM.xY.Rft   = Rft;

return

% plot responses
%--------------------------------------------------------------------------
spm_dcm_tfm_response(DCM.xY,DCM.xY.pst,DCM.xY.Hz);


text:
function k = rbfardKernDiagCompute(kern, x)

% RBFARDKERNDIAGCOMPUTE Compute diagonal of RBFARD kernel.
% FORMAT
% DESC computes the diagonal of the kernel matrix for the automatic relevance determination radial basis function kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : a vector containing the diagonal of the kernel matrix
% computed at the given points.
%
% SEEALSO : rbfardKernParamInit, kernDiagCompute, kernCreate, rbfardKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

k = repmat(kern.variance, size(x, 1), 1);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [] = gsp_plot_jtv_filter(G, filters, filtertype, param)
%GSP_PLOT_JTV_FILTER  Plot a time-vertex filterbank
%   Usage:  gsp_plot_jtv_filter(G,filters);
%           gsp_plot_jtv_filter(G,filters,param);
%
%   Input parameters:
%       G          : Time-Vertex graph structure
%       filters    : Cell array of time-vertex filters
%       filtertype : Filter domain (ts,js,ts-array,js-array)
%       param      : Structure of optional parameters
%   Output parameters:
%       none
%
%   Example:::
%
%         alpha = [0.1 0.5 1 2];
%         G = gsp_sensor(100);
%         G = gsp_jtv_graph(G,100,1);
%         G = gsp_estimate_lmax(G);
%         [g, filtertype] = gsp_jtv_design_wave(G, alpha);
%         param.domain='time-spectral';
%         gsp_plot_jtv_filter(G, g, filtertype,param);
%         param.domain='joint-spectral';
%         gsp_plot_jtv_filter(G, g, filtertype,param);
%
%
%   Additional parameters
%   ---------------------
%
%   * *param.npoints*  : Number of points where the filters are evaluated if eigenvalues not available (default 100).
%   * *param.show_sum* : Extra plot showing the sum of the squared magnitudes of the filters (default 1 if there is multiple filters).
%   * *param.verbose*  : Verbosity level (1 display the warning - 0 no log) (default 1).
%   * *param.title*    : Cell array of title for subplots. (default 1:Nf)
%   * *param.domain*   : Visualize the spectrum in 'time-spectral' or 'joint-spectral' (default param.domain='joint-spectral')
%

% Author: Francesco Grassi, Nathanael Perraudin
% Date   : September 2016

% Read input parameters
if nargin < 4
    param = struct;
end

if nargin<3
    error('Invalid number of arguments: GSP_PLOT_JTV_FILTER needs the type of time-vertex filter.')
end

Nf=numel(filters);

if ~isfield(param,'show_sum'), param.show_sum = Nf>1; end
if ~isfield(param,'npoints'),  param.npoints = 100; end
if ~isfield(param,'title'),    param.title=num2cell(1:Nf); end
if ~isfield(param,'domain'),   param.domain='joint-spectral'; end
if ~isfield(param,'fftshift'), param.fftshift = 1;end
%% Define axis

if isfield(G,'e')
    lambdas = G.e;
else
    lambdas = linspace(0,G.lmax,param.npoints);
end

switch filtertype
    case  {'ts','ts-array'}
        v = gsp_jtv_ta(G);
    case  {'js','js-array'}
        v = gsp_jtv_fa(G);
    otherwise
        error('Unknown filtertype');
end

%% Evaluating filter

fid = gsp_jtv_filter_evaluate(filters,filtertype,lambdas,v,param);

if param.show_sum
    Nf = Nf+1;
    test_sum = sum(fid.^2,3);
    fid(:,:,Nf) = test_sum;
    param.title{Nf} = 'Sum squared coeff';
end

%% Plot

switch param.domain
    case 'time-spectral'
        [v,xlab] = gsp_jtv_ta(G);
        fid = real(fid);
    case 'joint-spectral'
        [v,xlab] = gsp_jtv_fa(G,param.fftshift);
        fid = abs(fid);
        if param.fftshift
            fid = fftshift(fid,2);
        end
    otherwise
        error('Unknown domain');
end

figure
r = max(1,floor(sqrt(Nf)));
c = max(1,ceil(Nf/r));
for ii=1:Nf
    subplot(r,c,ii)
    imagesc(v,lambdas,fid(:,:,ii))
    xlabel(xlab);
    ylabel('\lambda');
    title(param.title{ii})
    axis xy
end

end

repo: gspbox


text:
function [datac,datafit,Amps,freqs]=rmlinesmovingwinc(data,movingwin,tau,params,p,plt,f0,fscanbw)
% fits significant sine waves to data (continuous data) using overlapping windows.
%
% Usage: [datac,datafit]=rmlinesmovingwinc(data,movingwin,tau,params,p,plt)
%
%  Inputs:
% Note that units of Fs, fpass have to be consistent.
%       data        (data in [N,C] i.e. time x channels/trials or as a single vector) - required.
%       movingwin         (in the form [window winstep] i.e length of moving
%                                                 window and step size)
%                                                 Note that units here have
%                                                 to be consistent with
%                                                 units of Fs - required
%       tau      parameter controlling degree of smoothing for the amplitudes - we use the
%       function 1-1/(1+exp(-tau*(x-Noverlap/2)/Noverlap) in the region of overlap to smooth
%       the sinewaves across the overlap region. Noverlap is the number of points
%       in the overlap region. Increasing tau leads to greater overlap smoothing,
%       typically specifying tau~10 or higher is reasonable. tau=1 gives an almost
%       linear smoothing function. tau=100 gives a very steep sigmoidal. The default is tau=10.
%       params      structure containing parameters - params has the
%       following fields: tapers, Fs, fpass, pad
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%                    Note that T has to be equal to movingwin(1).
%
%	        Fs 	        (sampling frequency) -- optional. Defaults to 1.
%               fpass       (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%	    p		    (P-value to calculate error bars for) - optional.
%	    Defaults to 0.05/Nwin where Nwin is length of window which
%	    corresponds to a false detect probability of approximately 0.05.
%       plt         (y/n for plot and no plot respectively) - default no
%                   plot.
%       f0          frequencies at which you want to remove the
%                   lines - if unspecified the program uses the f statistic
%                   to determine appropriate lines.
%
%  Outputs:
%       datafit        (fitted sine waves)
%       datac          (cleaned up data)
if nargin < 2; error('Need data and window parameters'); end;
if nargin < 4 || isempty(params); params=[]; end;
if nargin < 8, fscanbw = [];  end

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params); % set defaults for params
clear err trialave
if nargin < 6; plt='n'; end;
%
% Window,overlap and frequency information
%
data=change_row_to_column(data);
[N,C]=size(data);
Nwin=round(Fs*movingwin(1)); % number of samples in window
Nstep=round(movingwin(2)*Fs); % number of samples to step through
Noverlap=Nwin-Nstep; % number of points in overlap
%
% Sigmoidal smoothing function
%
if nargin < 3 || isempty(tau); tau=10; end; % smoothing parameter for sigmoidal overlap function
x=(1:Noverlap)';
smooth=1./(1+exp(-tau.*(x-Noverlap/2)/Noverlap)); % sigmoidal function
smooth=repmat(smooth,[1 C]);
%
% Start the loop
%
if nargin < 5 || isempty(p); p=0.05/Nwin; end % default for p value
if nargin < 7 || isempty(f0); f0=[]; end; % empty set default for f0 - uses F statistics to determine the frequencies
params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers
winstart=1:Nstep:N-Nwin+1;
nw=length(winstart);
datafit=zeros(winstart(nw)+Nwin-1,C);
Amps=cell(1,nw);
freqs=cell(1,nw);
for n=1:nw;
   indx=winstart(n):winstart(n)+Nwin-1;
   datawin=data(indx,:);
   [datafitwin,as,fs]=fitlinesc(datawin,params,p,'n',f0,fscanbw);
   Amps{n}=as;
   freqs{n}=fs;
   datafitwin0=datafitwin;
   if n>1; datafitwin(1:Noverlap,:)=smooth.*datafitwin(1:Noverlap,:)+(1-smooth).*datafitwin0(Nwin-Noverlap+1:Nwin,:);end;
   datafit(indx,:)=datafitwin;
end;
datac=data(1:size(datafit,1),:)-datafit;
if strcmp(plt,'y');
    [S,f]=mtspectrumsegc(data,movingwin(1),params);
    [Sc,fc]=mtspectrumsegc(datac,movingwin(1),params);
    plot(f,10*log10(S),fc,10*log10(Sc));
end;

author: VisLab

repo: EEG-Clean-Tools


text:
%IMWRITE  Saves an image to a specified file
%
%     cv.imwrite(filename, img)
%     cv.imwrite(filename, imgs)
%     cv.imwrite(filename, img, 'OptionName', optionValue, ...)
%     success = cv.imwrite(...)
%
% ## Input
% * __filename__ Name of the file.
% * __img__ Image to be saved.
% * __imgs__ A cell-array of images, to be saved as one multi-page image
%   (currently only supported by Tiff encoder).
%
% ## Output
% * __success__ optional output flag, true on success, false otherwise. If not
%   requested, the function throws an error on fail.
%
% ## Options
% * __FlipChannels__ in case the input is color image, flips the color order
%   from MATLAB's RGB/RGBA to OpenCV's BGR/BGRA order. default true
%
% The following format-specific save parameters are currently supported:
%
% * __JpegQuality__ For JPEG, it can be a quality from 0 to 100 (the higher is
%   the better). Default value is 95.
% * __JpegProgressive__ Enable JPEG features, 0 or 1, default is false.
% * __JpegOptimize__ Enable JPEG features, 0 or 1, default is false.
% * __JpegResetInterval__ JPEG restart interval, 0 - 65535, default is 0
%   (no restart).
% * __JpegLumaQuality__ Separate luma quality level, 0 - 100, default is 0
%   (don't use).
% * __JpegChromaQuality__ Separate chroma quality level, 0 - 100, default is 0
%   (don't use).
% * __PngCompression__ For PNG, it can be the compression level from 0 to 9.
%   A higher value means a smaller size and longer compression time. If
%   specified, `PngStrategy` is changed to `Default`. Default value is 1
%   (best speed setting).
% * __PngStrategy__ For PNG; used to tune the compression algorithm. These
%   flags will be modify the way of PNG image compression and will be passed
%   to the underlying zlib processing stage. The strategy parameter only
%   affects the compression ratio but not the correctness of the compressed
%   output even if it is not set appropriately. One of:
%   * __Default__ Use this value for normal data.
%   * __Filtered__ Use this value for data produced by a filter (or predictor).
%     Filtered data consists mostly of small values with a somewhat random
%     distribution. In this case, the compression algorithm is tuned to
%     compress them better. The effect of `Filtered` is to force more Huffman
%     coding and less string matching; it is somewhat intermediate between
%     `Default` and `HuffmanOnly`.
%   * __HuffmanOnly__ Use this value to force Huffman encoding only
%     (no string match).
%   * __RLE__ (default) Use this value to limit match distances to one
%     (run-length encoding). `RLE` is designed to be almost as fast as
%     `HuffmanOnly`, but give better compression for PNG image data.
%   * __Fixed__ Using this value prevents the use of dynamic Huffman codes,
%     allowing for a simpler decoder for special applications.
% * __PngBilevel__ Binary level PNG, 0 or 1, controls packing of pixels per
%   bytes. If false, PNG files pack pixels of bit-depths 1, 2, and 4 into
%   bytes as small as possible. default is false.
% * __PxmBinary__ For PPM, PGM, or PBM, it can be a binary format flag, 0 or 1,
%   to specify ASCII or binary encoding. default is true.
% * __ExrType__ override EXR storage type. Note that the EXR encoder only
%   accepts `single` images. One of:
%   * __Half__ store as half precision (FP16), see cv.convertFp16
%   * __Float__ store as single precision (FP32), this is the default.
% * __WebpQuality__ For WEBP, it can be a quality from 1 to 100 (the higher is
%   the better). By default (without any parameter) and for quality above 100
%   the lossless compression is used.
% * __PamTupleType__ For PAM, sets the TUPLETYPE field to the corresponding
%   string value that is defined for the format. One of:
%   * __Null__
%   * __BlackWhite__
%   * __Grayscale__
%   * __GrayscaleAlpha__
%   * __RGB__
%   * __RGBA__
%
% For advanced uses, you can directly pass a vector of parameters:
%
% * __Params__ Format-specific save parameters encoded as pairs:
%   `[paramId_1, paramValue_1, paramId_2, paramValue_2, ...]`.
%
% The function cv.imwrite saves the image to the specified file. The image
% format is chosen based on the filename extension (see cv.imread for the list
% of extensions). Only 8-bit (or 16-bit unsigned `uint16` in case of PNG,
% JPEG 2000, and TIFF) single-channel or 3-channel (with RGB channel order)
% images can be saved using this function. If the format, depth or channel
% order is different, use cv.cvtColor to convert it before saving. Or, use the
% universal cv.FileStorage I/O functions to save the image to XML or YAML
% format.
%
% (If the chosen encoder does not support the depth of the input image, the
% image will be implicitly cast to 8-bit).
%
% If the image cannot be saved (because of IO errors, improper permissions,
% unsupported or invalid format), the function throws an error.
%
% It is possible to store PNG images with an alpha channel using this function.
% To do this, create 8-bit (or 16-bit) 4-channel image RGBA, where the alpha
% channel goes last. Fully transparent pixels should have alpha set to 0,
% fully opaque pixels should have alpha set to 255/65535.
% Note: WEBP format can also store images with alpha channel.
% See `rgba_png_demo.m` sample for an example.
%
% See also: cv.imread, cv.imencode, imwrite, dicomwrite, hdrwrite
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function create_all_menus(obj, force)
    % create menus for main zmap figure
    % create_all_menus() - will create all menus, if they don't exist
    % create_all_menus(force) - will delete and recreate menus if force is true
    import zmaptopo.pltopo

    h = findobj(obj.fig, 'Tag', 'mainmap_menu_divider');
    if ~exist('force', 'var')
        force = false;
    end
    if ~isempty(h) && force
        delete(h); h = [];
    end
    if isempty(h)
        add_menu_divider(obj.fig, 'mainmap_menu_divider');
    end
    create_overlay_menu();

    catmen = obj.catalog_menu(force);
    create_random_data_simulations_menu(catmen);

    create_decluster_menu(findobj(obj.fig, 'Label', 'Catalog', '-and', 'type', 'uimenu'));

    add_sampling_menu(obj);
    create_ztools_menu();

    uimenu('Label', '|>', 'Enable', 'off');%
    funcsymbol = char(402);
    create_time_analysis_menu();
    create_map_analysis_menu();
    create_xsec_analysis_menu();
    create_3d_analysis_menu();

    uimenu('Label', '<|>', 'Enable', 'off');
    create_results_menu();
    uimenu('Label', '<|', 'Enable', 'off');

    % modify the file menu to add ZMAP stuff
    hFileMenu = findall(obj.fig, 'tag', 'figMenuFile');
    if force
        if ~isempty(hFileMenu)
            delete(findobj(hFileMenu,'Label','Get/Load Catalog'))
        end
        copymenus(findobj(obj.fig, 'Label', 'Get/Load Catalog'), hFileMenu);
    elseif isempty(hFileMenu)
        copymenus(findobj(obj.fig, 'Label', 'Get/Load Catalog'), hFileMenu);
    end
    addPreferencesMenuItem();
    addQuitMenuItem();
    addAboutMenuItem();

    if ZmapGlobal.Data.debug
        mainhelp = findall(obj.fig, 'Tag', 'figMenuHelp');
        uimenu(mainhelp, 'Label', 'Export ZmapMainWindow to workspace',...
            'Separator', 'on', 'MenuSelectedFcn', @export_me);
    end

    function export_me(~,~)
        prompt = 'Variable name to use';
        name = 'Export ZmapMainWindow to desktop';
        numlines = 1;
        defaultanswer = {'zmw'};
        name = inputdlg(prompt, name, numlines, defaultanswer);
        if ~isempty(name)
            assignin('base', name{:}, obj);
        end
    end

    %% time analysis_menu
    function create_time_analysis_menu()
        submenu = uimenu('Label', [funcsymbol '(t)']);
        uimenu(submenu, 'Label', 'functions evaluated through time', 'enable', 'off');
        uimenu(submenu, 'Label', 'Analyze time series ...',...
            'Separator', 'on',...
            'MenuSelectedFcn', @analyze_time_series_cb);

        function analyze_time_series_cb(~,~)
            % pick which time series we are investigating
            if ~isempty(obj.shape)
                items = ["Selected Events (IN polygon)", "Unselected Events (OUTSIDE polygon)"];
                items_data = {@()obj.catalog, @()obj.rawcatalog.subset(~obj.shape.isinterior(obj.rawcatalog.X, obj.rawcatalog.Y))};
            else
                items = ["Selected Events"]; %#ok<NBRAK>
                items_data = {@()obj.catalog};
            end
            if ~isempty(obj.XSectionTitles)
                items(end+1 : end + numel(obj.XSectionTitles)) = strcat("XSEC: ", string(obj.XSectionTitles));
                data_start = numel(items_data);
                for i = numel(obj.XSectionTitles) : -1 : 1
                    items_data(data_start+i) = {@()obj.xscats(obj.XSectionTitles{i}) };
                end
            end
            items(end+1) = "FULL (raw) Catalog";
            items_data(end+1) = {@()obj.rawcatalog};
            [selection, ok] = listdlg('PromptString', 'Select catalog to analyze',...
                'SelectionMode', 'single',...
                'ListString', items);
            if ok
                % c = items_data{selection};
                ctp = CumTimePlot(items_data{selection}() );
                ctp.plot();
            end
        end
    end

    %% map-view analysis menu
    % analyze items according to spacing in a horizontal plane
    function create_map_analysis_menu()
        submenu = uimenu('Label', [funcsymbol '(x,y)']); % 'Map');

        uimenu(submenu, 'Label', 'functions evaluated over a planar grid', 'enable', 'off');

        import XYfun.* % the map functions exist in the XYfun package
        % AB menu
        bvalgrid.AddMenuItem(submenu, @()obj.map_zap, 'Separator', 'on');
        bvalmapt.AddMenuItem(submenu, @()obj.map_zap);
        bdepth_ratio.AddMenuItem(submenu, @()obj.map_zap);

        % P menu
        bpvalgrid.AddMenuItem(submenu, @()obj.map_zap, 'Separator', 'on');
        rcvalgrid_a2.AddMenuItem(submenu, @()obj.map_zap, 'Enable', 'off');

        % Rate Change menu
        comp2periodz.AddMenuItem(submenu, @()obj.map_zap, 'Separator', 'on');
        uimenu(submenu, 'Label', 'Calculate a z-value map', 'Enable', 'off', 'MenuSelectedFcn', @(~,~)inmakegr(obj.catalog));

        % Quarry menu : detect quarry contamination
        findquar.AddMenuItem(submenu, @()obj.map_zap, 'Separator', 'on');

        %uimenu(submenu, 'Label', 'Map stress tensor', 'MenuSelectedFcn', @(~,~)stressgrid());
        stressgrid.AddMenuItem(submenu, @()obj.map_zap, 'Separator', 'on');
        %{
            uimenu(tmp, 'Label', 'Load...', 'Enable', 'off', 'MenuSelectedFcn',  @(~,~)rcvalgrid_a2('lo'));
            uimenu(submenu, 'Label', 'Load a z-value movie (map-view)', 'MenuSelectedFcn', @(~,~)loadmovz());
        %}
    end

    %% xsec-view analysis menu
    % analyze items according to position & depth along strike.

    function create_xsec_analysis_menu()
        submenu = uimenu('Label', XSection.MainMenuLabel);
        uimenu(submenu, 'Label', 'functions evaluated along a cross-section', 'enable', 'off');
        uimenu(submenu, 'Separator', 'on',...
            'Label', 'Define a cross-section', 'MenuSelectedFcn', @obj.cb_xsection, 'Tag', 'CreateXsec');

        import XZfun.* % the cross-section functions exist in the XZfun package

        magrcros.AddMenuItem(submenu, @()obj.xsec_zap, 'Separator','on');
        rc_cross_a2.AddMenuItem(submenu, @()obj.xsec_zap);

        uimenu(submenu, 'enable', 'off', 'Label', 'Calc a b-value cross-section', 'MenuSelectedFcn', @(~,~)nlammap(@()obj.xsec_zap));

        bcross.AddMenuItem(submenu, @()obj.xsec_zap, 'Enable', 'off');
        bcrossVt2.AddMenuItem(submenu, @()obj.xsec_zap, 'Enable', 'off');
        calc_Omoricross.AddMenuItem(submenu, @()obj.xsec_zap, 'Enable', 'off');
        calc_across.AddMenuItem(submenu, @()obj.xsec_zap, 'Enable', 'off');
        cross_stress.AddMenuItem(submenu, @()obj.xsec_zap, 'Enable', 'off');
    end

    %% 3D-view analysis menu
    % analyze items in a volume

    function create_3d_analysis_menu()
        submenu = uimenu('Label', [funcsymbol '(x,y,z)']); %'3D-Vol');
        uimenu(submenu, 'Label', 'functions over a 3-D grid', 'enable', 'off');
        uimenu(submenu, 'Separator', 'on', 'Label', 'Nothing here yet', 'Enable', 'off');

        import XYZfun.* % the cross-section functions exist in the XYZfun package

        bgrid3dB.AddMenuItem(submenu, @()obj.map_zap, 'Enable', 'off');
        zgrid3d.AddMenuItem(submenu, @()obj.map_zap, 'Enable', 'off');

        %{
            uimenu(submenu, 'Label', 'Load a 3D b-value grid', 'MenuSelectedFcn', @(~,~)myslicer('load'));
        %}

    end

    function create_results_menu()
        ResultsDisplay.Hdisplay.add_menus();
        set(findobj(allchild(gcf), 'flat', 'Tag', 'lookmenu'), 'Enable', 'off');
    end

    %% individual menus to create

    function create_overlay_menu()
        h = findobj(obj.fig, 'Tag', 'mainmap_menu_overlay');
        axm = obj.map_axes;
        if ~isempty(h) && exist('force', 'var') && force
            delete(h); h = [];
        end
        if ~isempty(h)
            return
        end

        % Make the menu to change symbol size and type
        %
        mapoptionmenu = uimenu(obj.fig, 'Label', 'Display', 'Tag', 'mainmap_menu_overlay');

        uimenu(mapoptionmenu, 'Label', '3-D view',...
            'MenuSelectedFcn', @obj.set_3d_view); % callback was plot3d

        uimenu(mapoptionmenu, 'Label', 'Change Z aspect ratio', 'Tag', 'z_aspectratio_menu',...
            'MenuSelectedFcn', @obj.set_z_ratio, 'Enable', 'off'); % callback was plot3d

        if ~iscartesian(obj.refEllipsoid)
            uimenu(mapoptionmenu, 'Label', 'Set aspect ratio by latitude',...
                'MenuSelectedFcn', @set_aspects_cb, ...
                'Checked', char(ZmapGlobal.Data.lock_aspect));

            if ZmapGlobal.Data.lock_aspect
                daspect(obj.get_all_map_axes(), [1 cosd(mean(axm.YLim)) 10]);
            end
        end
        m = uimenu(mapoptionmenu, 'Label', 'Toggle Lat/Lon Grid',...
            'MenuSelectedFcn', @toggle_grid,...
            'checked', char(ZmapGlobal.Data.mainmap_grid));

        if iscartesian(obj.refEllipsoid)
             axis(obj.get_all_map_axes, 'equal');
             m.Label = 'Toggle X/Y Grid';
        end

        grid(axm, char(ZmapGlobal.Data.mainmap_grid));

        % choose what to plot by
        for j = 1:numel(obj.rawcatalog.FieldnamesForColorby)
            myfn = obj.rawcatalog.FieldnamesForColorby{j};
            um(j) = uimenu(mapoptionmenu, 'Label', ['Color by ' myfn], 'MenuSelectedFcn', @(s,v)set_colorby(s, v, myfn) ,...
                'Checked', tf2onoff(strcmp(obj.colorField, myfn)));
            if j == 1
                um(j).Separator = 'on';
            end
        end
        uimenu(mapoptionmenu, 'Label', 'Edit Main Map Symbols', 'MenuSelectedFcn', @manage_mainmap_symbols);
        uimenu(mapoptionmenu, 'Label', 'Edit Current Map Symbols', 'MenuSelectedFcn', @manage_symbols_for_current_map);
        % add_symbol_menu(axm, mapoptionmenu, 'Map Symbols');
        ovmenu = uimenu(mapoptionmenu, 'Label', 'Layers');
        uimenu(ovmenu, 'Label', 'Show unselected events', 'MenuSelectedFcn', @toggle_unselected_events,...
            'Checked', "on");

        try
            MapFeature.foreach(obj.Features, 'addToggleMenu', ovmenu, axm)
        catch ME
            warning(ME.message)
        end

        uimenu(ovmenu, 'Label', 'Plot stations + station names',...
            'Separator', 'on',...
            'MenuSelectedFcn', @(~,~)plotstations(axm));

        lemenu = uimenu(mapoptionmenu, 'Label', 'Legend by ...  ', 'Enable', 'off');

        uimenu(lemenu, 'Label', 'Change legend breakpoints',...
            'MenuSelectedFcn', @change_legend_breakpoints);
        legend_types = {'Legend by time', 'tim';...
            'Legend by depth', 'depth';...
            'Legend by magnitudes', 'mag';...
            'Mag by size, Depth by color', 'mad' %;...
            % 'Symbol color by faulting type', 'fau'...
            };

        for i = 1:size(legend_types, 1)
            m = uimenu(lemenu, 'Label', legend_types{i, 1},...
                'MenuSelectedFcn', @(s,v)cb_plotby(s, v, legend_types{i, 2}) );
            if i == 1
                m.Separator = 'on';
            end
        end
        clear legend_types

        uimenu(mapoptionmenu, 'Label', 'Change font size ...', 'Enable', 'off',...
            'MenuSelectedFcn', @change_map_fonts);

        uimenu(mapoptionmenu, 'Label', 'Change background colors',...
            'MenuSelectedFcn', @(~,~)setcol, 'Enable', 'off'); %

        uimenu(mapoptionmenu,...
            'Label', ['Mark large event with M > ' num2str(ZmapGlobal.Data.CatalogOpts.BigEvents.MinMag)],...
            'MenuSelectedFcn', @cb_plot_large_quakes);

        uimenu(mapoptionmenu,...
            'Separator', 'on',...
            'Label', 'Close all result tabs', 'MenuSelectedFcn', @clear_result_tabs);

        uimenu(mapoptionmenu, 'label', 'Redraw',...
            'Separator', 'on',...
            'MenuSelectedFcn', @obj.cb_redraw);

        uimenu(mapoptionmenu, 'label', 'Open printable figure', ...
            'MenuSelectedFcn', @(~,~)make_printable_figure_copy(obj.fig));

        function set_aspects_cb(src, ev)
            callbacks.toggle_aspectratio(src, ev, obj.get_all_map_axes(), "SetGlobal");
        end

        function clear_result_tabs(~, ~)
            delete(obj.maingroup.Children(obj.maingroup.Children ~= obj.maintab))
        end
        function manage_symbols_for_current_map(src, ev)
            %

            % look for plots, but only direct children of the selected tab.
            % this way we don't access axes such as the x-section plots
            ax = findobj(obj.maingroup.SelectedTab.Children, 'flat', 'Type', 'axes');
            SymbolManager.cb(src, ev, ax);
        end

        function manage_mainmap_symbols(src, ev)
            obj.make_map_active('notify');
            SymbolManager.cb(src, ev, axm);
        end

        function toggle_unselected_events(src,~)
            h_events = findobj(obj.map_axes, 'Tag', 'all events');
            if h_events.Visible == "on"
                h_events.Visible = 'off';
                src.Checked = 'off';
            else
                h_events.Visible = 'on';
                src.Checked = 'on';
            end

        end

        function set_colorby(~, ~, val)
            obj.colorField = val;
            obj.mainEventProps.ColorBy = val;
            obj.make_map_active('notify')

            % update menus
            for jj = 1:numel(um)
                myfn = obj.rawcatalog.FieldnamesForColorby{jj};
                um(jj).Checked = tf2onoff(strcmp(obj.colorField, myfn));
            end
            h = findobj(obj.fig, 'Type', 'colorbar', '-and', 'Parent', obj.fig);

            delete(h)
            obj.plotmainmap();
        end

        function cb_plotby(~,~, s)
            ZG = ZmapGlobal.Data;
            ZG.mainmap_plotby = s;
            watchon;
            watchoff;
        end

        function cb_plot_large_quakes(src,~)
            ZG = ZmapGlobal.Data;
            [~, ~, ZG.CatalogOpts.BigEvents.MinMag] = smart_inputdlg('Choose magnitude threshold',...
                struct('prompt', 'Mark events with M > ? ', 'value', ZG.CatalogOpts.BigEvents.MinMag));
            src.Label = "Mark large event with M > " + ZG.CatalogOpts.BigEvents.MinMag;
            obj.bigEvents = obj.catalog.subset(obj.catalog.Magnitude > ZG.CatalogOpts.BigEvents.MinMag);
            obj.CatalogManager.ChangeFilter('big events', @(c) c.Magnitude > ZmapGlobal.Data.CatalogOpts.BigEvents.MinMag);
            set(findobj(obj.fig, 'Tag', 'big events'), 'DisplayName', src.Label);
        end

        function toggle_grid(src, ~)
            src.Checked = toggleOnOff(src.Checked);
            arrayfun(@(x)grid(x, src.Checked) , obj.get_all_map_axes())
            drawnow nocallbacks
        end

    end

    function create_ztools_menu()
        h = findobj(obj.fig, 'Tag', 'mainmap_menu_ztools');
        if ~isempty(h) && exist('force', 'var') && force
            delete(h); h = [];
        end
        if ~isempty(h)
            return
        end
        submenu = uimenu('Label', 'ZTools', 'Tag', 'mainmap_menu_ztools');

        uimenu(submenu, 'Label', 'Plot Current Map into Projection',...
            'MenuSelectedFcn', @(~,~)create_projectedmap_from_mainmap(obj.fig));

        create_topo_map_menu(submenu);
        uimenu(submenu,'Label', 'Random data simulations has moved to CATALOG menu', 'Enable', 'off');
        % create_random_data_simulations_menu(submenu);

        % create_histogram_menu(submenu);
        create_explore_catalog_menu(submenu);

        uimenu(submenu, 'Label', 'Misfit calculation', 'MenuSelectedFcn', @(~,~)inmisfit(),...
            'Enable', 'off'); %FIXME: misfitcalclulation poorly documented, not sure what it is comparing.

    end
    function create_topo_map_menu(parent)
        submenu   =  uimenu(parent, 'Label', 'Plot topographic map');
        uimenu(submenu, 'Label', 'Open a Web Map Display', 'MenuSelectedFcn', @(~,~)webmap_of_catalog(obj.catalog, true));
        uimenu(submenu, 'Label', 'Plot Topography on main map', 'MenuSelectedFcn', @add_topography_to_main_map);
        uimenu(submenu, 'Label', 'Plot Swiss Topography on main map', 'MenuSelectedFcn', @(s,v)add_topography_to_main_map(s, v, 'CH') );

        % FIXME the following need to be found and fixed
        %{
        uimenu(submenu, 'Label', 'Open DEM GUI', 'MenuSelectedFcn', @(~,~)zmaptopo.prepinp());
        uimenu(submenu, 'Label', '3 arc sec resolution (USGS DEM)', 'MenuSelectedFcn', @(~,~)pltopo('lo3'));
        uimenu(submenu, 'Label', '30 arc sec resolution (GLOBE DEM)', 'MenuSelectedFcn', @(~,~)pltopo('lo1'));
        uimenu(submenu, 'Label', '30 arc sec resolution (GTOPO30)', 'MenuSelectedFcn', @(~,~)pltopo('lo30'));
        uimenu(submenu, 'Label', '2 deg resolution (ETOPO 2)', 'MenuSelectedFcn', @(~,~)pltopo('lo2'));
        uimenu(submenu, 'Label', '5 deg resolution (ETOPO 5, Terrain Base)', 'MenuSelectedFcn', @(~,~)pltopo('lo5'));
        uimenu(submenu, 'Label', 'Your topography (mydem, mx, my must be defined)', 'MenuSelectedFcn', @(~,~)pltopo('yourdem'));
        uimenu(submenu, 'Label', 'Help on plotting topography', 'MenuSelectedFcn', @(~,~)pltopo('genhelp'));
        %}
    end

    function add_topography_to_main_map(~, ~, code)
        htopo = findobj(obj.map_axes, '-regexp', 'Tag', 'topographic_map.*');
        delete(htopo);

        if exist('code', 'var')
            zmaptopo.add_topo(obj.map_axes, 'locale', code);
            % zmaptopo.add_topo(obj.map_axes, 'locale', code, 'ShadedOnly', true);
        end
        % now plot the world topography too.
        zmaptopo.add_topo(obj.map_axes, 'locale', 'world');
        % zmaptopo.add_topo(obj.map_axes, 'locale', 'world', 'ShadedOnly', true);
    end

    function create_random_data_simulations_menu(parent)
        submenu  =   uimenu(parent, 'Label', 'Create Synthetic/Permuted catalog', 'Separator', 'on');
        uimenu(submenu, 'label', 'Create permutated catalog (also new b-value)...', 'MenuSelectedFcn', @cb_create_permutated);
        uimenu(submenu, 'label', 'Create synthetic catalog...', 'MenuSelectedFcn', @cb_create_synhthetic_cat);

        % uimenu(submenu, 'Label', 'Evaluate significance of b- and a-values', 'MenuSelectedFcn', @(~,~)brand());
        %  uimenu(submenu, 'Label', 'Calculate a random b map and compare to observed data', 'MenuSelectedFcn', @(~,~)brand2());
        uimenu(submenu, 'Label', 'Info on synthetic catalogs', 'MenuSelectedFcn', @(~,~)web(['file:' ZmapGlobal.Data.hodi '/zmapwww/syntcat.htm']));

        function cb_create_permutated(~,~)
            % will replace existing primary catalog
            [rand_catalog, ok] = syn_invoke_random_dialog(obj.catalog);

            if ok && ~isempty(rand_catalog)
                ZmapMainWindow(figure, rand_catalog);
            end
        end

        function cb_create_synhthetic_cat(~,~)
            % will replace existing primary catalog
            [syn_cat, ok] = syn_invoke_dialog(obj.catalog);
            if ok
                ZmapMainWindow(figure, syn_cat);
            end
        end

    end

    function create_explore_catalog_menu(parent)
        uimenu(parent, 'Separator', 'on',...
            'Label', 'Explore Catalog', 'MenuSelectedFcn', @explore_catalog);

        function explore_catalog(~,~)
            t = findOrCreateTab(obj.fig, obj.maingroup, obj.maingroup, 'Exploration', 'deleteable');
            ax = findobj(t.Children, 'Type', 'axes');
            if isempty(ax)
                ax = axes(t);
            end
            cep = CatalogExplorationPlot(ax, @get_catalog);
            cep.scatter('explore catalog')
            legend(ax, 'show', 'Location', 'northeast');
            title(ax, sprintf('Exploring : catalog: %s , with %d events',...
                strrep(obj.catalog.Name, '_', '\_'), obj.catalog.Count));
            bringToForeground(ax);
        end

        function c = get_catalog()
            c = obj.catalog;
        end

    end

    function create_decluster_menu(parent)
        submenu = parent;% uimenu(parent, 'Label', 'Decluster the catalog');

        uimenu(submenu, 'Label', 'Decluster (Reasenberg)', 'MenuSelectedFcn', @cb_reasen,...
            'Separator', 'on');
        uimenu(submenu, 'Label', 'Decluster (Gardner & Knopoff)',...
            ...'Enable', 'off',...
            'MenuSelectedFcn', @cb_declus_inp);
        % uimenu(submenu, 'Label', 'Decluster (Zaliapin)', 'MenuSelectedFcn', @cb_zaliapin);
    end
    function cb_declus_inp(~,~)
        [out, eMethod] = declus_inp(obj.catalog);
        if isempty(eMethod)
            return
        end
        % Return Decluster results to the main window
        disp('Summary of Gardiner-Knopoff clusters using window:' + string(eMethod))
        disp('Cluster #     Number of Events (inc. mainshock)')
        summary(categorical(out.allClusterIdx(out.allClusterIdx ~= 0)))

        assignin('base', 'gk_decluster_output', out);
        disp(out.description)

        choice = questdlg({'Replace Main Catalog with:','',...
            sprintf('\\bfDeclustered\\rm : background + mainshocks [%d evts]', out.declusteredCatalog.Count),...
            sprintf('\\bfMainshocks\\rm: only mainshocks [%d evts]', sum(out.mainshockClusterIdx ~= 0)),...
            sprintf('\\bfClusters\\rm: only fore- and aftershocks) [%d evts]',out.aftershockCatalog.Count),...
            '','\itClose dialog to keep original catalog\rm',...
            '','Decluster results have been written to workspace as \itgk\_decluster\_output\rm'},...
            'Declust: close to keep orig cat',... % dialog title
            'Declustered','Mainshocks', 'Clusters',... % up to 3 buttons
            struct('Default','Declustered','Interpreter','tex'));  % Default choice
        switch choice
            case 'Declustered'
                obj.rawcatalog = out.declusteredCatalog;
                %replaceMainCatalog(out.declusteredCatalog)
            case 'Mainshocks'
                obj.rawcatalog = obj.catalog.subset(out.mainshockClusterIdx ~=0);
                obj.rawcatalog.Name = "GK Declust ("+string(eMethod)+") " + obj.rawcatalog.Name
                %replaceMainCatalog(obj.catalog.subset(out.mainshockClusterIdx ~=0))
            case 'Clusters'
                obj.rawcatalog = out.aftershockCatalog;
                %replaceMainCatalog(out.aftershockCatalog)
        end
        figure(obj.fig)
        obj.replot_all()
    end

    function cb_reasen(~,~)
        rdc = ReasenbergDeclusterClass(obj.catalog, "CalcFinishedFcn", @update_window_with_declust);

        function update_window_with_declust()
            if ~isempty(rdc.declusteredCatalog)
                msg.dbdisp('replacing the catalog')
                obj.rawcatalog = rdc.declusteredCatalog;
                obj.CatalogManager.RawCatalog = obj.rawcatalog;
                obj.replot_all()
            else
                errordlg('Empty declustered catalog, Main window will not be updated')
            end
        end
    end

    function cb_zaliapin(~,~)
        rdc = ZaliapinDeclusterClass(obj.catalog, "CalcFinishedFcn", @update_window_with_declust);

        function update_window_with_declust()
            if ~isempty(rdc.declusteredCatalog)
                msg.dbdisp('replacing the catalog')
                obj.rawcatalog = rdc.declusteredCatalog;
                obj.CatalogManager.RawCatalog = obj.rawcatalog;
                obj.replot_all()
            else
                errordlg('Empty declustered catalog, Main window will not be updated')
            end
        end
    end
end

author: CelsoReyes