text:
function [filt, B, A] = ft_preproc_highpassfilter(dat, Fs, Fhp, N, type, dir, instabilityfix, df, wintype, dev, plotfiltresp, usefftfilt)

% FT_PREPROC_HIGHPASSFILTER applies a high-pass filter to the data and thereby removes
% the low frequency components in the data
%
% Use as
%   [filt] = ft_preproc_highpassfilter(dat, Fsample, Fhp, N, type, dir, instabilityfix)
% where
%   dat             data matrix (Nchans X Ntime)
%   Fs              sampling frequency in Hz
%   Fhp             filter frequency in Hz
%   order           optional filter order, default is 6 (but) or dependent on frequency band and data length (fir/firls)
%   type            optional filter type, can be
%                     'but'       Butterworth IIR filter (default)
%                     'firws'     FIR filter with windowed sinc
%                     'fir'       FIR filter using MATLAB fir1 function
%                     'firls'     FIR filter using MATLAB firls function (requires MATLAB Signal Processing Toolbox)
%                     'brickwall' frequency-domain filter using forward and inverse FFT
%   dir             optional filter direction, can be
%                     'onepass'                   forward filter only
%                     'onepass-reverse'           reverse filter only, i.e. backward in time
%                     'onepass-zerophase'         zero-phase forward filter with delay compensation (default for firws, linear-phase symmetric FIR only)
%                     'onepass-reverse-zerophase' zero-phase reverse filter with delay compensation
%                     'onepass-minphase'          minimum-phase converted forward filter (non-linear, only for firws)
%                     'twopass'                   zero-phase forward and reverse filter (default, except for firws)
%                     'twopass-reverse'           zero-phase reverse and forward filter
%                     'twopass-average'           average of the twopass and the twopass-reverse
%   instabilityfix  optional method to deal with filter instabilities
%                     'no'       only detect and give error (default)
%                     'reduce'   reduce the filter order
%                     'split'    split the filter in two lower-order filters, apply sequentially
%   df              optional transition width (firws)
%   wintype         optional window type (firws), can be
%                     'hamming' (default)    maximum passband deviation 0.0022 [0.22%], stopband attenuation -53dB
%                     'hann'                 maximum passband deviation 0.0063 [0.63%], stopband attenuation -44dB
%                     'blackman'             maximum passband deviation 0.0002 [0.02%], stopband attenuation -74dB
%                     'kaiser'
%   dev             optional max passband deviation/stopband attenuation (only for firws with kaiser window, default = 0.001 [0.1%, -60 dB])
%   plotfiltresp    optional, 'yes' or 'no', plot filter responses (only for firws, default = 'no')
%   usefftfilt      optional, 'yes' or 'no', use fftfilt instead of filter (only for firws, default = 'no')
%
% Note that a one- or two-pass filter has consequences for the strength of the filter,
% i.e. a two-pass filter with the same filter order will attenuate the signal twice as
% strong.
%
% Further note that the filter type 'brickwall' filters in the frequency domain,
% but may have severe issues. For instance, it has the implication that the time
% domain signal is periodic. Another issue pertains to that frequencies are
% not well defined over short time intervals; particularly for low frequencies.
%
% See also PREPROC

% Copyright (c) 2003-2022, Robert Oostenveld, Arjen Stolk, Andreas Widmann,
% Jan-Mathijs Schoffelen
%
% 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$

% determine the size of the data
[nchans, nsamples] = size(dat);

% set the default filter order later
if nargin<4 || isempty(N)
  N = [];
end

% set the default filter type
if nargin<5 || isempty(type)
  type = 'but';
end

% set the default filter direction
if nargin<6 || isempty(dir)
  if strcmp(type, 'firws')
    dir = 'onepass-zerophase';
  else
    dir = 'twopass';
  end
end

% set the default method to deal with filter instabilities
if nargin<7|| isempty(instabilityfix)
  instabilityfix = 'no';
end

% Set default transition width later
if nargin < 8 || isempty(df)
  df = [];
end

% Set default window type
if nargin < 9 || isempty(wintype)
  wintype = 'hamming';
end

% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 10 || isempty(dev)
  if strcmp(wintype, 'kaiser')
    dev = 0.001;
  else
    dev = [];
  end
end

% Set default passband deviation/stopband attenuation for Kaiser window
if nargin < 11 || isempty(plotfiltresp)
  plotfiltresp = 'no';
end

% Set default filter function
if nargin < 12 || isempty(usefftfilt)
  usefftfilt = false;
else
  % convert to boolean value
  usefftfilt = istrue(usefftfilt);
end

% Filtering does not work on integer data
if ~isa(dat, 'double') && ~isa(dat, 'single')
  dat = cast(dat, 'double');
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

% Nyquist frequency
Fn = Fs/2;

% demean the data before filtering
meandat = nanmean(dat,2);
dat = bsxfun(@minus, dat, meandat);

% compute filter coefficients
switch type
  case 'but'
    if isempty(N)
      N = 6;
    end
    [B, A] = butter(N, max(Fhp)/Fn, 'high');

  case 'firws'
    % Input arguments
    if length(Fhp) ~= 1
      ft_error('One cutoff frequency required.')
    end

    % Filter order AND transition width set?
    if ~isempty(N) && ~isempty(df)
      ft_warning('firws:dfOverridesN', 'Filter order AND transition width set - transition width setting will override filter order.')
    elseif isempty(N) && isempty(df) % Default transition width heuristic
      df = fir_df(Fhp, Fs);
    end

    % Compute filter order from transition width
    [foo, maxDf] = fir_df(Fhp, Fs); %#ok<ASGLU>
    isOrderLow = false;
    if ~isempty(df)
      if df > maxDf
        ft_error('Transition band too wide. Maximum transition width is %.2f Hz.', maxDf)
      end
      [N, dev] = firwsord(wintype, Fs, df, dev);
    else % Check filter order otherwise
      [df, dev] = invfirwsord(wintype, Fs, N, dev);
      if df > maxDf
        nOpt = firwsord(wintype, Fs, maxDf, dev);
        ft_warning('firws:filterOrderLow', 'Filter order too low. For better results a minimum filter order of %d is recommended. Effective cutoff frequency might deviate from requested cutoff frequency.', nOpt)
        isOrderLow = true;
      end
    end

    % Window
    if strcmp(wintype, 'kaiser')
      beta = kaiserbeta(dev);
      win = windows('kaiser', N + 1, beta);
    else
      win = windows(wintype, N + 1);
    end

    % Impulse response
    B = firws(N, Fhp / Fn, 'high', win);
    A = 1;

    % Convert to minimum phase
    if strcmp(dir, 'onepass-minphase')
      B = minphaserceps(B);
    end

    % Twopass filtering
    if strncmp(dir, 'twopass', 7)
      pbDev = (dev + 1)^2 - 1;
      sbAtt = 20 * log10(dev^2);
      order = 2 * (length(B) - 1);
      isTwopass = true;
    else
      pbDev = dev;
      sbAtt = 20 * log10(dev);
      order = length(B) - 1;
      isTwopass = false;
    end

    % Reporting
    ft_info once
    ft_info('Highpass filtering data: %s, order %d, %s-windowed sinc FIR\n', dir, order, wintype);
    if ~isTwopass && ~isOrderLow % Do not report shifted cutoffs
      ft_info('  cutoff (-6 dB) %g Hz\n', Fhp);
      tb = [max([Fhp - df / 2 0]), min([Fhp + df / 2 Fn])]; % Transition band edges
      ft_info('  transition width %.1f Hz, stopband 0-%.1f Hz, passband %.1f-%.0f Hz\n', df, tb, Fn);
    end
    if ~isOrderLow
      ft_info('  maximum passband deviation %.4f (%.2f%%), stopband attenuation %.0f dB\n', pbDev, pbDev * 100, sbAtt);
    end

  case 'fir'
    if isempty(N)
      N = 3*fix(Fs / Fhp);
      if rem(N,2)==1,   N=N+1;    end
    end
    if N > floor( (size(dat,2) - 1) / 3)
      N=floor(size(dat,2)/3) - 2;
      if rem(N,2)==1,   N=N+1;    end
    end
    B = fir1(N, max(Fhp)/Fn, 'high');
    A = 1;

  case 'firls' % from NUTMEG's implementation
    % Deprecated: see bug 2453
    ft_warning('The filter type you requested is not recommended for neural signals, only proceed if you know what you are doing.')
    if isempty(N)
      N = 3*fix(Fs / Fhp);
      if rem(N,2)==1,   N=N+1;    end
    end
    if N > floor( (size(dat,2) - 1) / 3)
      N=floor(size(dat,2)/3) - 2;
      if rem(N,2)==1,   N=N+1;    end
    end
    f = 0:0.001:1;
    if rem(length(f),2)~=0
      f(end)=[];
    end
    z = zeros(1,length(f));
    [val,pos1] = min(abs(Fs*f/2 - Fhp));
    pos2 = length(f);
    z(pos1:pos2) = 1;
    A = 1;
    B = firls(N,f,z); % requires MATLAB signal processing toolbox

  case 'brickwall'
    ax = (0:(size(dat,2)-1))./(Fs/size(dat,2)); % frequency axis

    a     = ones(1, size(dat,2));
    fbin1 = nearest(ax, Fhp);
    fbin2 = nearest(ax, Fs-Fhp);

    a(1:(fbin1-1))   = 0;
    a((fbin2+1):end) = 0;

    f    = fft(dat,[],2);             % FFT
    f    = f.*a(ones(size(dat,1),1),:); % brickwall
    filt = real(ifft(f,[],2));        % iFFT

  otherwise
    ft_error('unsupported filter type "%s"', type);
end

% Plot filter responses
if strcmp(plotfiltresp, 'yes')
  plotfresp(B, A, [], Fs, dir)
end

if ~isequal(type, 'brickwall')
  try
    filt = filter_with_correction(B,A,dat,dir,usefftfilt);
  catch
    switch instabilityfix
      case 'no'
        rethrow(lasterror);
      case 'reduce'
        ft_warning('off','backtrace');
        ft_warning('filter instability detected - reducing the %dth order filter to an %dth order filter', N, N-1);
        ft_warning('on','backtrace');
        filt = ft_preproc_highpassfilter(dat,Fs,Fhp,N-1,type,dir,instabilityfix);
      case 'split'
        N1 = ceil(N/2);
        N2 = floor(N/2);
        ft_warning('off','backtrace');
        ft_warning('filter instability detected - splitting the %dth order filter in a sequential %dth and a %dth order filter', N, N1, N2);
        ft_warning('on','backtrace');
        filt = ft_preproc_highpassfilter(dat ,Fs,Fhp,N1,type,dir,instabilityfix);
        filt = ft_preproc_highpassfilter(filt,Fs,Fhp,N2,type,dir,instabilityfix);
      otherwise
        ft_error('incorrect specification of instabilityfix');
    end % switch
  end
end


text:
%Recursive Zonal Equal Area Sphere Partitioning: Utilities
%
% Recursive Zonal Equal Area (EQ) Sphere Partitioning Toolbox.
% Release 1.10 2005-06-01
%
%Functions
%=========
%
%  area_of_cap            Area of spherical cap
%  area_of_collar         Area of spherical collar
%  area_of_ideal_region   Area of one region of an EQ partition
%  area_of_sphere         Area of sphere
%  cart2polar2            Convert Cartesian to spherical polar coordinates on S^2
%  euc2sph_dist           Convert Euclidean to spherical distance
%  euclidean_dist         Euclidean distance between two points
%  fatcurve               Create a parameterized cylindrical surface
%  haslight               Check if axis handle has a light attached
%  ideal_collar_angle     Ideal angle for spherical collars of an EQ partition
%  illustration_options   Options for illustrations of EQ partitions
%  partition_options      Options for EQ partition
%  polar2cart             Convert spherical polar to Cartesian coordinates
%  sph2euc_dist           Convert spherical to Euclidean distance
%  spherical_dist         Spherical distance between two points on the sphere
%  sradius_of_cap         Spherical radius of spherical cap of given area
%  volume_of_ball         Volume of the unit ball

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Function changed name from e2s to euc2sph_dist
% Function changed name from s2e to sph2euc_dist
% Function changed name from s2x to polar2cart
% Function changed name from x2s2 to cart2polar2
% Add new function fatcurve
% Add new function haslight
% Clean up descriptions
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-13 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

repo: gspbox


text:
function varargout = CoilDesignPanel(varargin)
% CoilDESIGNPANEL M-file for CoilDesignPanel.fig
%      CoilDESIGNPANEL, by itself, creates a new CoilDESIGNPANEL or raises the existing
%      singleton*.
%
%      H = CoilDESIGNPANEL returns the handle to a new CoilDESIGNPANEL or the handle to
%      the existing singleton*.
%
%      CoilDESIGNPANEL('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CoilDESIGNPANEL.M with the given input arguments.
%
%      CoilDESIGNPANEL('Property','Value',...) creates a new CoilDESIGNPANEL or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before CoilDesignPanel_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to CoilDesignPanel_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help CoilDesignPanel

% Last Modified by GUIDE v2.5 17-Apr-2014 19:26:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @CoilDesignPanel_OpeningFcn, ...
    'gui_OutputFcn',  @CoilDesignPanel_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before CoilDesignPanel is made visible.
function CoilDesignPanel_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to CoilDesignPanel (see VARARGIN)
global VCtl;
global VObj;
global VMco;

handles.Simuh=varargin{1};
% Load coil element list
try
    handles.CoilElemListStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'CoilElem' filesep 'CoilElem.xml']);
    guidata(hObject,handles);
    Root=DoConvStruct2Tree(handles.CoilElemListStruct);
    handles.CoilElemListTree=uitree('v0', 'Root',Root,'SelectionChangeFcn', {@ChkCoilElemListTreeNode,handles});
    set(handles.CoilElemListTree.getUIContainer,'Units','normalized');
    set(handles.CoilElemListTree.getUIContainer,'Position',[0.0,0.4,0.2,0.58]);
    handles.Attrh1=[];
catch ME
    errordlg('Coil element list can not be loaded!');
    close(hObject);
    return;
end

%Load tabs for coil simulation
try
    handles.CoilSimuAttrStruct=DoParseXML([handles.Simuh.MRiLabPath filesep 'Macro' filesep 'CoilElem' filesep 'CoilSimuAttr.xml']);
catch ME
    errordlg('CoilSimuAttr.xml file is missing or can not be loaded!');
    close(hObject);
    return;
end

handles.Setting_tabgroup=uitabgroup(handles.Setting_uipanel);
guidata(handles.CoilDesignPanel_figure,handles);
for i=1:length(handles.CoilSimuAttrStruct.Children)
    eval(['handles.' handles.CoilSimuAttrStruct.Children(i).Name '_tab=uitab( handles.Setting_tabgroup,' '''title'',' '''' handles.CoilSimuAttrStruct.Children(i).Name ''',''Units'',''normalized'');']);
    eval(['DoEditValue(handles,handles.' handles.CoilSimuAttrStruct.Children(i).Name '_tab,handles.CoilSimuAttrStruct.Children(' num2str(i) ').Attributes,2,[0.2,0.15,0.0,0.15,0.15,0.05]);']);
    handles=guidata(handles.CoilDesignPanel_figure);
end

% Create grid and display structure
VMco.xdimres=str2num(get(handles.Attrh2.('XDimRes'),'String'));
VMco.ydimres=str2num(get(handles.Attrh2.('YDimRes'),'String'));
VMco.zdimres=str2num(get(handles.Attrh2.('ZDimRes'),'String'));

Mxdims=size(VObj.Rho);
[VMco.xgrid,VMco.ygrid,VMco.zgrid]=meshgrid((-(Mxdims(2)-1)/2)*VObj.XDimRes:VMco.xdimres:((Mxdims(2)-1)/2)*VObj.XDimRes,...
    (-(Mxdims(1)-1)/2)*VObj.YDimRes:VMco.ydimres:((Mxdims(1)-1)/2)*VObj.YDimRes,...
    (-(Mxdims(3)-1)/2)*VObj.ZDimRes:VMco.zdimres:((Mxdims(3)-1)/2)*VObj.ZDimRes);
VMco.axes=handles.CoilSen_axes;
handles.IV=struct(...
    'xslice',(handles.Simuh.SV.Slice - (Mxdims(2)-1)/2)* VObj.XDimRes,...
    'yslice',(handles.Simuh.CV.Slice - (Mxdims(1)-1)/2)* VObj.YDimRes,...
    'zslice',(handles.Simuh.AV.Slice - (Mxdims(3)-1)/2)* VObj.ZDimRes ...
    );
view(handles.CoilSen_axes,3);

if verLessThan('matlab','8.5')
    % Code to run in MATLAB R2015a and earlier here
    %Open current loaded Coil XML file
    if strcmp(get(handles.Simuh.CoilTx_tab,'Visible'),'on')
        if isfield(handles.Simuh,'CoilTxXMLFile')
            CoilXMLFile = handles.Simuh.CoilTxXMLFile;
        end
    end

    if strcmp(get(handles.Simuh.CoilRx_tab,'Visible'),'on')
        if isfield(handles.Simuh,'CoilRxXMLFile')
            CoilXMLFile = handles.Simuh.CoilRxXMLFile;
        end
    end
else
    % Code to run in MATLAB R2015a and later here
    tabtitle=get(get(handles.Simuh.CoilSel_tabgroup,'SelectedTab'),'Title');
    if strcmp(tabtitle,'Tx')
        if isfield(handles.Simuh,'CoilTxXMLFile')
            CoilXMLFile = handles.Simuh.CoilTxXMLFile;
        end
    elseif strcmp(tabtitle,'Rx')
        if isfield(handles.Simuh,'CoilRxXMLFile')
            CoilXMLFile = handles.Simuh.CoilRxXMLFile;
        end
    end
end

if exist('CoilXMLFile','var')
    import javax.swing.*
    import javax.swing.tree.*;
    handles.CoilStruct=DoParseXML(CoilXMLFile);
    handles.CoilXMLFile=CoilXMLFile;
    DoWriteXML2m(handles.CoilStruct,[handles.CoilXMLFile(1:end-3) 'm']);
    %     guidata(hObject, handles);
    Root=DoConvStruct2Tree(handles.CoilStruct);
    handles.CoilTree=uitree('V0','Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
    set(handles.CoilTree.getUIContainer,'Units','normalized');
    set(handles.CoilTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
    handles.CoilTreeModel=DefaultTreeModel(Root);
    handles.CoilTree.setModel(handles.CoilTreeModel);
    try
        DoPlotCoilSen(handles);
        set(handles.Undock_pushbutton,'Enable','on');
        set(handles.Display_pushbutton,'Enable','on');
        set(handles.Forwards_pushbutton,'Enable','on');
        set(handles.Backwards_pushbutton,'Enable','on');
        handles=guidata(hObject);
    catch me
        error_msg{1,1}='ERROR!!! Displaying B1 field aborted.';
        error_msg{2,1}=me.message;
        errordlg(error_msg);
    end

end

% Choose default command line output for CoilDesignPanel
handles.input=varargin;
handles.output=hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes CoilDesignPanel wait for user response (see UIRESUME)
% uiwait(handles.CoilDesignPanel_figure);

function ChkCoilElemListTreeNode(tree, value, handles)

handles=guidata(handles.CoilDesignPanel_figure);
SelNode=tree.SelectedNodes;
SelNode=SelNode(1);
Level=SelNode.getValue;
if ~ischar(Level)
    LevelChg=diff(Level,1);
    Level([LevelChg; -1]==1)=[];
end

Node='handles.CoilElemListStruct';
if Level=='0'
    eval(['handles.CoilElemNode=' Node ';']);
else
    for i=1:length(Level)
        Node=[Node '.Children(' num2str(Level(i)) ')'];
    end
    eval(['handles.CoilElemNode=' Node ';']);
end
guidata(handles.CoilDesignPanel_figure, handles);

% --- Outputs from this function are returned to the command line.
function varargout = CoilDesignPanel_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = hObject;

% --------------------------------------------------------------------
function LoadCoil_uimenu_Callback(hObject, eventdata, handles)
% hObject    handle to LoadCoil_uimenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function LoadCoilFile_uimenu_Callback(hObject, eventdata, handles)
% hObject    handle to LoadCoilFile_uimenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import javax.swing.*
import javax.swing.tree.*;

[filename,coilpathname,filterindex]=uigetfile({'Coil*.xml','XML-files (*.xml)'},'MultiSelect','off');
if filename~=0
    if ~isempty(handles.Attrh1)
        Attrh1name=fieldnames(handles.Attrh1);
        for i=1:length(Attrh1name)
            delete(handles.Attrh1.(Attrh1name{i}));
        end
        handles.Attrh1=[];
    end
    handles.CoilStruct=DoParseXML(fullfile(coilpathname,filename));
    handles.CoilXMLFile=fullfile(coilpathname,filename);
    guidata(hObject, handles);
    Root=DoConvStruct2Tree(handles.CoilStruct);
    handles.CoilTree=uitree('v0', 'Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
    set(handles.CoilTree.getUIContainer,'Units','normalized');
    set(handles.CoilTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
    handles.CoilTreeModel=DefaultTreeModel(Root);
    handles.CoilTree.setModel(handles.CoilTreeModel);
    cla(handles.CoilSen_axes);
    clear global VCco;
    % Add searchig path
    path(path,coilpathname);
else
    errordlg('No Coil is loaded!');
    return;
end
guidata(hObject, handles);

function ChkNode(tree, value, handles)

handles=guidata(handles.CoilDesignPanel_figure);
SelNode=tree.SelectedNodes;
SelNode=SelNode(1);
Level=SelNode.getValue;
if ~ischar(Level)
    LevelChg=diff(Level,1);
    Level([LevelChg; -1]==1)=[];
end
Node='handles.CoilStruct';
if Level=='0'
    eval(['handles.CoilNode=' Node ';']);
else
    for i=1:length(Level)
        Node=[Node '.Children(' num2str(Level(i)) ')'];
    end
    eval(['handles.CoilNode=' Node ';']);
end
handles.CoilNodeLvl=Node;
handles.SelNode=SelNode;

if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:length(Attrh1name)
        delete(handles.Attrh1.(Attrh1name{i}));
    end
    handles.Attrh1=[];
end

set(handles.CoilAttri_uipanel,'Title',['''' char(SelNode.getName) ''' Coil Attribute'],'Unit','normalized');
if SelNode.isLeaf | SelNode.isRoot
    if strcmp(char(SelNode.getName),'Specials')
        DoEditValue(handles,handles.CoilAttri_uipanel,handles.CoilNode.Attributes,1,[0.25,0.25,0.0,0.1,0.1,0.05]);
    else
        DoEditValue(handles,handles.CoilAttri_uipanel,handles.CoilNode.Attributes,1,[0.25,0.25,0.0,0.1,0.1,0.05]);
    end
else
    guidata(handles.CoilDesignPanel_figure, handles);
end

% --- Executes on button press in DelNode_pushbutton.
function DelNode_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to DelNode_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.CoilNode.Name
        case 'MRiLabCoil'
            errordlg('MRiLabCoil Root can not be deleted!');
            return;
    end

    eval([handles.CoilNodeLvl '=[];']);
    nN=handles.SelNode.getNextSibling;
    while ~isempty(nN)
        Level=nN.getValue;
        Level(end)=[];
        nN.setValue(Level);
        nN=nN.getNextSibling;
    end
    guidata(hObject, handles);
    pause(0.1);

    % remove node
    handles.CoilTree.setSelectedNode(handles.SelNode.getParent);
    handles.CoilTreeModel.removeNodeFromParent(handles.SelNode);

    DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

% --- Executes on button press in AddNode_pushbutton.
function AddNode_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to AddNode_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import javax.swing.*
import javax.swing.tree.*;

try
    ChkCoilElemListTreeNode(handles.CoilElemListTree, [], handles); % In case CoilElem.xml updated
    handles=guidata(handles.CoilDesignPanel_figure);
    switch handles.CoilNode.Name
        case {'MRiLabCoil'}
            % add node
            SelNode=handles.CoilElemListTree.SelectedNodes;
            SelNode=SelNode(1);
            AddedNode=SelNode.clone;
            switch char(AddedNode.getName)
                case {'CoilElem'}
                    errordlg([char(AddedNode.getName) ' can not be added!']);
                    return;
            end
            handles.CoilTreeModel.insertNodeInto(AddedNode,handles.SelNode,handles.SelNode.getChildCount());
            handles.CoilTree.setSelectedNode(AddedNode); % expand to show added child
            handles.CoilTree.setSelectedNode(handles.SelNode); % insure additional nodes are added to parent

            nP=AddedNode.getPreviousSibling;
            if ~isempty(nP)
                Level=nP.getValue;
                AddedNode.setValue([Level; Level(end)+1]);
                eval([handles.CoilNodeLvl '.Children(' num2str(Level(end)+1) ')=handles.CoilElemNode;']);
                DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
            else
                nP=AddedNode.getParent;
                Level=nP.getValue;
                AddedNode.setValue([Level; 1]);
                eval([handles.CoilNodeLvl '.Children=handles.CoilElemNode;']);
                DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
            end
        otherwise
            errordlg(['Coil Element can not be added under ' handles.CoilNode.Name '!']);
            return;
    end
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

% --- Executes on button press in UpdateCoilXML_pushbutton.
function UpdateCoilXML_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to UpdateCoilXML_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% update Coil XML file
if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:2:length(Attrh1name)
        switch get(handles.Attrh1.(Attrh1name{i+1}),'Style')
            case 'edit'
                handles.CoilNode.Attributes((i+1)/2).Value=get(handles.Attrh1.(Attrh1name{i+1}),'String');
            case 'checkbox'
                if get(handles.Attrh1.(Attrh1name{i+1}),'Value')
                    handles.CoilNode.Attributes((i+1)/2).Value='^1';
                else
                    handles.CoilNode.Attributes((i+1)/2).Value='^0';
                end
            case 'popupmenu'
                handles.CoilNode.Attributes((i+1)/2).Value(2)=num2str(get(handles.Attrh1.(Attrh1name{i+1}),'Value'));
        end
    end
    eval([handles.CoilNodeLvl '=handles.CoilNode;']);
    DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
    % update associated m function
    DoWriteXML2m(handles.CoilStruct,[handles.CoilXMLFile(1:end-3) 'm']);
end

% update CoilSimuAttr.xml
tabs=get(handles.Setting_tabgroup,'Children');
for j=1:length(handles.CoilSimuAttrStruct.Children)
    Attrh2=get(tabs(j),'Children');
    for i=1:2:length(Attrh2)
        if ~iscell(get(Attrh2(end-i),'String'))
            handles.CoilSimuAttrStruct.Children(j).Attributes((i+1)/2).Value=get(Attrh2(end-i),'String');
        else
            handles.CoilSimuAttrStruct.Children(j).Attributes((i+1)/2).Value(2)=num2str(get(Attrh2(end-i),'Value'));
        end
    end
end
DoWriteXML(handles.CoilSimuAttrStruct,[handles.Simuh.MRiLabPath filesep 'Macro' filesep 'CoilElem' filesep 'CoilSimuAttr.xml']);

guidata(hObject, handles);

% --- Executes on button press in CopSelNode_pushbutton.
function CopNode_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to CopSelNode_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.CoilNode.Name
        case 'MRiLabCoil'
            errordlg('MRiLabCoil Root can not be copied!');
            return;
    end
    handles.CopSelNode=handles.SelNode.clone;
    handles.CopCoilNode=handles.CoilNode;
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

% --- Executes on button press in PasNode_pushbutton.
function PasNode_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to PasNode_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

import javax.swing.*
import javax.swing.tree.*;

try
    switch handles.CoilNode.Name
        case {'MRiLabCoil'}
            % paste node
            handles.CoilTreeModel.insertNodeInto(handles.CopSelNode,handles.SelNode,handles.SelNode.getChildCount());
            handles.CoilTree.setSelectedNode(handles.CopSelNode); % expand to show added child
            handles.CoilTree.setSelectedNode(handles.SelNode); % insure additional nodes are added to parent

            nP=handles.CopSelNode.getPreviousSibling;
            if ~isempty(nP)
                Level=nP.getValue;
                handles.CopSelNode.setValue([Level; Level(end)+1]);
                eval([handles.CoilNodeLvl '.Children(' num2str(Level(end)+1) ')=handles.CopCoilNode;']);
                DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
            else
                nP=handles.CopSelNode.getParent;
                Level=nP.getValue;
                handles.CopSelNode.setValue([Level; 1]);
                eval([handles.CoilNodeLvl '.Children=handles.CopCoilNode;']);
                DoWriteXML(handles.CoilStruct,handles.CoilXMLFile);
            end
        otherwise
            errordlg(['Coil Element can not be added under ' handles.CoilNode.Name '!']);
            return;
    end
    handles.CopSelNode=[];
    handles.CopCoilNode=[];
    guidata(hObject, handles);
catch me
    error_msg{1,1}='ERROR!!! Node operation aborted.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);

end

% --- Executes on button press in CoilExecute_pushbutton.
function CoilExecute_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to CoilExecute_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

set(hObject,'Enable','off');
set(hObject,'String','Calc...');
pause(0.01);
try
    DoPlotCoilSen(handles);
    set(handles.Undock_pushbutton,'Enable','on');
    set(handles.Display_pushbutton,'Enable','on');
    set(handles.Forwards_pushbutton,'Enable','on');
    set(handles.Backwards_pushbutton,'Enable','on');
catch me
    error_msg{1,1}='Error occurred at calculating B1/E1 field.';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
    set(handles.Undock_pushbutton,'Enable','off');
    set(handles.Display_pushbutton,'Enable','off');
    set(handles.Forwards_pushbutton,'Enable','off');
    set(handles.Backwards_pushbutton,'Enable','off');
end
set(hObject,'Enable','on');
set(hObject,'String','Execute');

% --------------------------------------------------------------------
function SavePanel_uipushtool_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to SavePanel_uipushtool (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
DoSaveSnapshot(handles.CoilDesignPanel_figure);

% --- Executes when user attempts to close CoilDesignPanel_figure.
function CoilDesignPanel_figure_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to CoilDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

clear global VCco;
% Hint: delete(hObject) closes the figure
delete(hObject);

% --- Executes on button press in Forwards_pushbutton.
function Forwards_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Forwards_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global VObj;
TAttr=get(handles.Attrh2.('Plane'),'String');
switch TAttr{get(handles.Attrh2.('Plane'),'Value')}
    case 'XY'
        handles.IV.zslice=handles.IV.zslice+VObj.ZDimRes;
    case 'XZ'
        handles.IV.yslice=handles.IV.yslice+VObj.YDimRes;
    case 'YZ'
        handles.IV.xslice=handles.IV.xslice+VObj.XDimRes;
end
DoUpdateSlice(handles.CoilSen_axes,handles.F,handles.IV,'Coil');
guidata(hObject, handles);

% --- Executes on button press in Backwards_pushbutton.
function Backwards_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Backwards_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global VObj;
TAttr=get(handles.Attrh2.('Plane'),'String');
switch TAttr{get(handles.Attrh2.('Plane'),'Value')}
    case 'XY'
        handles.IV.zslice=handles.IV.zslice-VObj.ZDimRes;
    case 'XZ'
        handles.IV.yslice=handles.IV.yslice-VObj.YDimRes;
    case 'YZ'
        handles.IV.xslice=handles.IV.xslice-VObj.XDimRes;
end
DoUpdateSlice(handles.CoilSen_axes,handles.F,handles.IV,'Coil');
guidata(hObject, handles);

% --------------------------------------------------------------------
function Cursor_uitoggletool_ClickedCallback(hObject, eventdata, handles)
% hObject    handle to Cursor_uitoggletool (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if ~isfield(handles,'dcm')
    handles.dcm=datacursormode(handles.CoilDesignPanel_figure);
    datacursormode on
    set(handles.dcm,'updatefcn',{@Cursor_update,handles});
else
    if strcmp(get(handles.dcm,'Enable'),'on')
        set(handles.dcm,'Enable','off');
    else
        set(handles.dcm,'Enable','on');
    end
end
guidata(hObject, handles);

function output_txt = Cursor_update(obj,event_obj,handles)
% Display the position of the data cursor
% obj          Currently not used (empty)
% event_obj    Handle to event object
% output_txt   Data cursor text string (string or cell array of strings).

global VObj;
global VCtl;
global VMco;

handles=guidata(handles.CoilDesignPanel_figure);
pos = get(event_obj,'Position');
pos(1)=round((pos(1)+VCtl.ISO(1)*VObj.XDimRes)/VMco.xdimres);
pos(2)=round((pos(2)+VCtl.ISO(2)*VObj.YDimRes)/VMco.ydimres);
pos(3)=round((pos(3)+VCtl.ISO(3)*VObj.ZDimRes)/VMco.zdimres);

output_txt = {['X: ',num2str(pos(1),4)],...
    ['Y: ',num2str(pos(2),4)],...
    ['Z: ',num2str(pos(3),4)],...
    ['Value: ',num2str(handles.F(pos(2),pos(1),pos(3)),4)]};

% --- Executes on button press in Undock_pushbutton.
function Undock_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Undock_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global VMco;
TAttr=get(handles.Attrh2.('Plane'),'String');
switch TAttr{get(handles.Attrh2.('Plane'),'Value')}
    case 'XY'
        zmax=max(VMco.zgrid(1,1,:));
        zmin=min(VMco.zgrid(1,1,:));
        zslice=handles.IV.zslice;
        if zslice<zmax & zslice>zmin
            zslice=floor((zslice-zmin)/VMco.zdimres)+1;
            figure;
            quiver(VMco.xgrid(:,:,zslice),VMco.ygrid(:,:,zslice),handles.Fx(:,:,zslice),handles.Fy(:,:,zslice));
            title('Field direction in XY plane');
            set(gca,'YDir','reverse');
            xlabel('X');  ylabel('Y');
            axis image;
        end
    case 'XZ'
        ymax=max(max(max(VMco.ygrid)));
        ymin=min(min(min(VMco.ygrid)));
        yslice=handles.IV.yslice;
        if yslice<ymax & yslice>ymin
            yslice=floor((yslice-ymin)/VMco.ydimres)+1;
            figure;
            quiver(VMco.xgrid(yslice,:,:),VMco.zgrid(yslice,:,:),handles.Fx(yslice,:,:),handles.Fz(yslice,:,:));
            title('Field direction in XZ plane');
            set(gca,'YDir','reverse');
            xlabel('X');  ylabel('Z');
            axis image;
        end
    case 'YZ'
        xmax=max(max(max(VMco.xgrid)));
        xmin=min(min(min(VMco.xgrid)));
        xslice=handles.IV.xslice;
        if xslice<xmax & xslice>xmin
            xslice=floor((xslice-xmin)/VMco.xdimres)+1;
            figure;
            quiver(VMco.ygrid(:,xslice,:),VMco.zgrid(:,xslice,:),handles.Fy(:,xslice,:),handles.Fz(:,xslice,:));
            title('Field direction in YZ plane');
            set(gca,'YDir','reverse');
            xlabel('Y');  ylabel('Z');
            axis image;
        end
end

% --------------------------------------------------------------------
function NewCoil_uimenu_Callback(hObject, eventdata, handles)
% hObject    handle to NewCoil_uimenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function NewCoilFile_uimenu_Callback(hObject, eventdata, handles)
% hObject    handle to NewCoilFile_uimenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

prompt = {'Enter NEW coil name:','Enter NEW coil note:'};
dlg_title = 'Input for creating new coil';
num_lines = 1;
def = {'Coil_0ChHead','0 Chanel Head Coil'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
if isempty(answer)
    warndlg('No new coil file is created!');
    return
end
if isempty(answer{1})
    warndlg('No new coil file is created!');
    return
end
if ~strcmp(answer{1}(1:5),'Coil_')
    warndlg('Coil name must have a prefix ''Coil_'', please rename your coil!');
    return
end
coilpath=uigetdir(pwd,'Specify a saving path for new coil file.');
if coilpath==0
    warndlg('You have to specify a saving path!');
    return;
end

% Create new coil file
mkdir([coilpath filesep answer{1}]); % make folder
copyfile([handles.Simuh.MRiLabPath filesep 'Config' filesep 'Coil' filesep 'Head' filesep 'Coil_1ChHead' filesep 'Coil_1ChHead.xml'],...
    [coilpath filesep answer{1} filesep answer{1} '.xml']); % copy coil temple
% Add searchig path
path(path,[coilpath filesep answer{1}]);
% Load new coil file
import javax.swing.*;
import javax.swing.tree.*;
if ~isempty(handles.Attrh1)
    Attrh1name=fieldnames(handles.Attrh1);
    for i=1:length(Attrh1name)
        delete(handles.Attrh1.(Attrh1name{i}));
    end
    handles.Attrh1=[];
end
handles.CoilStruct=DoParseXML([coilpath filesep answer{1} filesep answer{1} '.xml']);
for i=1:length(handles.CoilStruct.Attributes)
    if strcmp(handles.CoilStruct.Attributes(i).Name,'Name')
        handles.CoilStruct.Attributes(i).Value = answer{1};
    elseif strcmp(handles.CoilStruct.Attributes(i).Name,'Notes')
        handles.CoilStruct.Attributes(i).Value = answer{2};
    end
end
DoWriteXML(handles.CoilStruct,[coilpath filesep answer{1} filesep answer{1} '.xml']); %update new coil file
handles.CoilXMLFile=[coilpath filesep answer{1} filesep answer{1} '.xml'];
DoWriteXML2m(handles.CoilStruct,[handles.CoilXMLFile(1:end-3) 'm']);
guidata(hObject, handles);
Root=DoConvStruct2Tree(handles.CoilStruct);
handles.CoilTree=uitree('V0','Root',Root,'SelectionChangeFcn', {@ChkNode,handles});
set(handles.CoilTree.getUIContainer,'Units','normalized');
set(handles.CoilTree.getUIContainer,'Position',[0.205,0.4,0.2,0.58]);
handles.CoilTreeModel=DefaultTreeModel(Root);
handles.CoilTree.setModel(handles.CoilTreeModel);
cla(handles.CoilSen_axes);
clear global VCco;
guidata(hObject, handles);

% --- Executes on mouse press over figure background, over a disabled or
% --- inactive control, or over an axes background.
function CoilDesignPanel_figure_WindowButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to CoilDesignPanel_figure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

tempaxes = gca;
try
    if strcmp(get(gcf,'Selectiontype'),'alt')
        figure 			% Create a new figure
        ax_new = axes;
        copyaxes(tempaxes, ax_new);
    end
catch me
end

% --- Executes on button press in Display_pushbutton.
function Display_pushbutton_Callback(hObject, eventdata, handles)
% hObject    handle to Display_pushbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

Conts=get(handles.Attrh2.FieldType,'String');
SelCont=Conts{get(handles.Attrh2.FieldType,'Value')};
switch SelCont
    case 'B1Field'
        Conts=get(handles.Attrh2.Mode,'String');
        SelCont=Conts{get(handles.Attrh2.Mode,'Value')};
        assignin('base', ['B1_' SelCont], handles.F);
        MU_Matrix_Display(['B1_' SelCont],'Magnitude');
    case 'E1Field'
        assignin('base', 'E1_Magnitude', handles.F);
        MU_Matrix_Display('E1_Magnitude','Magnitude');
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function Offspring = OperatorGA(Problem,Parent,Parameter)
%OperatorGA - Crossover and mutation operators of genetic algorithm.
%
%   Off = OperatorGA(Pro,P) uses genetic operators to generate offsprings
%   for problem Pro based on parents P. If P is an array of SOLUTION
%   objects, then Off is also an array of SOLUTION objects. While if P is a
%   matrix of decision variables, then Off is also a matrix of decision
%   variables, i.e., the offsprings are not evaluated. P is split into two
%   subsets P1 and P2 with the same size, where each object or row of P1
%   and P2 is used to generate two offsprings. Different operators are used
%   for different encoding schemes.
%
%   Off = OperatorGA(Pro,P,{proC,disC,proM,disM}) specifies the parameters
%   of operators, where proC is the probability of crossover, disC is the
%   distribution index of simulated binary crossover, proM is the
%   expectation of the number of mutated variables, and disM is the
%   distribution index of polynomial mutation.
%
%   Example:
%       Offspring = OperatorGA(Problem,Parent)
%       Offspring = OperatorGA(Problem,Parent.decs,{1,20,1,20})
%
%   See also OperatorGAhalf

%------------------------------- Reference --------------------------------
% [1] K. Deb, K. Sindhya, and T. Okabe, Self-adaptive simulated binary
% crossover for real-parameter optimization, Proceedings of the Annual
% Conference on Genetic and Evolutionary Computation, 2007, 1187-1194.
% [2] K. Deb and M. Goyal, A combined genetic adaptive search (GeneAS) for
% engineering design, Computer Science and informatics, 1996, 26: 30-45.
% [3] L. Davis, Applying adaptive algorithms to epistatic domains,
% Proceedings of the International Joint Conference on Artificial
% Intelligence, 1985, 162-164.
% [4] D. B. Fogel, An evolutionary approach to the traveling salesman
% problem, Biological Cybernetics, 1988, 60(2): 139-144.
%------------------------------- 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".
%--------------------------------------------------------------------------

    if nargin > 2
        [proC,disC,proM,disM] = deal(Parameter{:});
    else
        [proC,disC,proM,disM] = deal(1,20,1,20);
    end
    if isa(Parent(1),'SOLUTION')
        evaluated = true;
        Parent    = Parent.decs;
    else
        evaluated = false;
    end
    Parent1   = Parent(1:floor(end/2),:);
    Parent2   = Parent(floor(end/2)+1:floor(end/2)*2,:);
    Offspring = zeros(2*size(Parent1,1),size(Parent1,2));
    Type      = arrayfun(@(i)find(Problem.encoding==i),1:5,'UniformOutput',false);
    if ~isempty([Type{1:2}])    % Real and integer variables
        Offspring(:,[Type{1:2}]) = GAreal(Parent1(:,[Type{1:2}]),Parent2(:,[Type{1:2}]),Problem.lower([Type{1:2}]),Problem.upper([Type{1:2}]),proC,disC,proM*length([Type{1:2}])/size(Parent1,2),disM);
    end
    if ~isempty(Type{3})        % Label variables
        Offspring(:,Type{3}) = GAlabel(Parent1(:,Type{3}),Parent2(:,Type{3}),Problem.lower(Type{3}),Problem.upper(Type{3}),proC,proM*length(Type{3})/size(Parent1,2));
    end
    if ~isempty(Type{4})        % Binary variables
        Offspring(:,Type{4}) = GAbinary(Parent1(:,Type{4}),Parent2(:,Type{4}),proC,proM*length(Type{4})/size(Parent1,2));
    end
    if ~isempty(Type{5})        % Permutation variables
        Offspring(:,Type{5}) = GApermutation(Parent1(:,Type{5}),Parent2(:,Type{5}),proC);
    end
    if evaluated
        Offspring = Problem.Evaluation(Offspring);
    end
end

function Offspring = GAreal(Parent1,Parent2,lower,upper,proC,disC,proM,disM)
% Genetic operators for real and integer variables

    %% Simulated binary crossover
    [N,D] = size(Parent1);
    beta  = zeros(N,D);
    mu    = rand(N,D);
    beta(mu<=0.5) = (2*mu(mu<=0.5)).^(1/(disC+1));
    beta(mu>0.5)  = (2-2*mu(mu>0.5)).^(-1/(disC+1));
    beta = beta.*(-1).^randi([0,1],N,D);
    beta(rand(N,D)<0.5) = 1;
    beta(repmat(rand(N,1)>proC,1,D)) = 1;
    Offspring = [(Parent1+Parent2)/2+beta.*(Parent1-Parent2)/2
                 (Parent1+Parent2)/2-beta.*(Parent1-Parent2)/2];

    %% Polynomial mutation
    Lower = repmat(lower,2*N,1);
    Upper = repmat(upper,2*N,1);
    Site  = rand(2*N,D) < proM/D;
    mu    = rand(2*N,D);
    temp  = Site & mu<=0.5;
    Offspring       = min(max(Offspring,Lower),Upper);
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                      (1-(Offspring(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp = Site & mu>0.5;
    Offspring(temp) = Offspring(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                      (1-(Upper(temp)-Offspring(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
end

function Offspring = GAlabel(Parent1,Parent2,lower,upper,proC,proM)
% Genetic operators for label variables

    %% Uniform crossover
    [N,D] = size(Parent1);
    k     = rand(N,D) < 0.5;
    k(repmat(rand(N,1)>proC,1,D)) = false;
    Offspring1    = Parent1;
    Offspring2    = Parent2;
    Offspring1(k) = Parent2(k);
    Offspring2(k) = Parent1(k);
    Offspring     = [Offspring1;Offspring2];

    %% Bitwise mutation
    Site = rand(2*N,D) < proM/D;
    Rand = round(unifrnd(repmat(lower,2*N,1),repmat(upper,2*N,1)));
    Offspring(Site) = Rand(Site);
end

function Offspring = GAbinary(Parent1,Parent2,proC,proM)
% Genetic operators for binary variables

    %% Uniform crossover
    [N,D] = size(Parent1);
    k     = rand(N,D) < 0.5;
    k(repmat(rand(N,1)>proC,1,D)) = false;
    Offspring1    = Parent1;
    Offspring2    = Parent2;
    Offspring1(k) = Parent2(k);
    Offspring2(k) = Parent1(k);
    Offspring     = [Offspring1;Offspring2];

    %% Bit-flip mutation
    Site = rand(2*N,D) < proM/D;
    Offspring(Site) = ~Offspring(Site);
end

function Offspring = GApermutation(Parent1,Parent2,proC)
% Genetic operators for permutation variables

    %% Order crossover
    [N,D]     = size(Parent1);
    Offspring = [Parent1;Parent2];
    k = randi(D,1,2*N);
    for i = 1 : N
        if rand < proC
            Offspring(i,k(i)+1:end)   = setdiff(Parent2(i,:),Parent1(i,1:k(i)),'stable');
            Offspring(i+N,k(i)+1:end) = setdiff(Parent1(i,:),Parent2(i,1:k(i)),'stable');
        end
    end

    %% Slight mutation
    k = randi(D,1,2*N);
    s = randi(D,1,2*N);
    for i = 1 : 2*N
        if s(i) < k(i)
            Offspring(i,:) = Offspring(i,[1:s(i)-1,k(i),s(i):k(i)-1,k(i)+1:end]);
        elseif s(i) > k(i)
            Offspring(i,:) = Offspring(i,[1:k(i)-1,k(i)+1:s(i)-1,k(i),s(i):end]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function TF = raisedCos (f,Tb,r)
% Raised Cosine Filter
% f = Frequency Matrix
% Tb = Bit Period
% r = Rolloff Factor

W0 = 1/(2*Tb);
W = (r+1)*W0;
x = abs(f)*(1/(W-W0));
x = x + ((W - 2*W0)/(W-W0))*ones(size(x));
TF = cos((pi/4)*x).^2;
xx = find(abs(f) < (2*W0 - W));
for ll = xx, TF(ll) = 1;end
xx = find(abs(f) > W);
for ll = xx, TF(ll) = 0;end

author: Sable

repo: mcbench-benchmarks


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function[x,r] = mfJacobi(x,b,para,steps)
%
% Jacobi smoother for multigrid, see mfvcycle
%==============================================================================

function[x,r] = mfJacobi(x,b,para,steps)

if nargin == 0;
  MGsolver;
  x = 'endOfMinimalExample';
  return;
end;

% n = length(b);
[r,D] = mfAy(x,para);
r     = b - r;

for i=1:steps,
  ss = D\r;
  x = x + para.MGomega*ss;   %x = rmnspace(x,para.Z);
  r = b - mfAy(x,para);      %r = b - rmnspace(mfAu(x,para),para.Z);
%     his(i) = norm(r)/norm(b);
%    figure(2); plot(r); pause
end;

%  figure(1); clf; plot(his);
%  mfilename, keyboard
%==============================================================================


text:
function [objFuncValue, dv, deltaV1, deltaV2, deltaV1R, deltaV2R, xfrOrbit, deltaV1NTW, deltaV2NTW] = rendezvousObjFunc(x, iniOrbit, finOrbit, gmuXfr, weights, onlyOptBurn1)
%rendezvousObjFunc Summary of this function goes here
%   Detailed explanation goes here
    time1 = x(1);
    time2 = time1 + x(2);

    iniOrbBodyInfo = getBodyInfoStructFromOrbit(iniOrbit);
    [rVect, vVect] = getStateAtTime(iniOrbBodyInfo, time1, iniOrbit(8));
    [sma, ecc, inc, raan, arg, iniTA] = getKeplerFromState(rVect,vVect,iniOrbit(8));
    iniOrbit(1) = sma;
    iniOrbit(2) = ecc;
    iniOrbit(3) = inc;
    iniOrbit(4) = raan;
    iniOrbit(5) = arg;
    iniOrbit(6) = computeMeanFromTrueAnom(iniTA, ecc);
    iniOrbit(7) = time1;

    iniOrbBodyInfo.mean = computeMeanFromTrueAnom(iniTA, ecc);
    iniOrbBodyInfo.epoch = time1;
    % iniOrbit(7) = [];
    % iniOrbit(6) = [];

    finOrbBodyInfo = getBodyInfoStructFromOrbit(finOrbit);
    [rVect, vVect] = getStateAtTime(finOrbBodyInfo, time2, finOrbit(8));
    [sma, ecc, inc, raan, arg, finTA] = getKeplerFromState(rVect,vVect,finOrbit(8));
    finOrbit(1) = sma;
    finOrbit(2) = ecc;
    finOrbit(3) = inc;
    finOrbit(4) = raan;
    finOrbit(5) = arg;
    finOrbit(6) = computeMeanFromTrueAnom(finTA, ecc);
    finOrbit(7) = time2;
    % finOrbit(7) = [];
    % finOrbit(6) = [];

    x(1) = iniTA;
    x(2) = finTA;
    x(3) = time2 - time1;
    [dv, deltaV1, deltaV2, deltaV1R, deltaV2R, xfrOrbit, deltaV1NTW, deltaV2NTW] = twoBurnOrbitChangeObjFunc(x, iniOrbit, finOrbit, gmuXfr);

    if(onlyOptBurn1)
        dv = norm(deltaV1);
    end

    dvWt = weights(1);
    timeWt = weights(2);
    timeNorm = computePeriod(mean([abs(iniOrbit(1)),abs(finOrbit(1))]), gmuXfr); %time norm is the period of the average of the initial and final orbit SMAs

    objFuncValue = dvWt*dv + timeWt*(time2 - time1)/timeNorm;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = p09_autonomous ( )

%*****************************************************************************80
%
%% P09_AUTONOMOUS reports whether problem P09 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [spectogram_y_lims, s, f] = cutSpectogramFrequency(s,f, handles)

    upper_freq = find(f>=handles.data.settings.HighFreq*1000,1);
    if isempty(upper_freq)
        upper_freq = length(f);
    end
    lower_freq = find(f>=handles.data.settings.LowFreq*1000,1);

    % Extract the region within the frequency range
    f = f(lower_freq:upper_freq,:);
    s = s(lower_freq:upper_freq,:);
    spectogram_y_lims = [ min(f), max(f)];
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
function z = entropy(x)
% Compute entropy z=H(x) of a discrete variable x.
% Input:
%   x: a integer vectors
% Output:
%   z: entropy z=H(x)
% Written by Mo Chen (sth4nth@gmail.com).
n = numel(x);
[~,~,x] = unique(x);
Px = accumarray(x, 1)/n;
Hx = -dot(Px,log2(Px));
z = max(0,Hx);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function poly = EulerCycles(F)
% retrieve Euler cycles from an not oriented edge list

% convert to consecute numbers
% now it remains to identify 1->2, 3->4, 5->6 and so on
[~,A] = sort(F(:));
B(A) = 1:numel(F);
B = reshape(B,[],2);

% index for B, i.e., B(iB(n)) == n
[~,iB] = sort(B(:));

numF = size(F,1);
polyB = zeros(numF,1); % these are positions in B
tourStart = 1;

if isempty(iB)
  poly = {};
  return;
end

currentiB = iB(1);  % starting vertex

for ipoly = 1:numF

  polyB(ipoly) = currentiB;
  currentB = B(currentiB)-1;
  iB(2*fix(currentB/2)+(1:2)) = 0; % do not visit this again

  if currentiB > numF
    nextB = B(currentiB - numF);
  else
    nextB = B(currentiB + numF);
  end

  if mod(nextB,2), nextB = nextB + 1; else nextB = nextB - 1; end

  currentiB = iB(nextB);

  if currentiB == 0 %start new cycle
    [~,~,currentiB] = find(iB,1);
    tourStart(end+1) = ipoly+1; %#ok<AGROW>
  end

end

% entries of poly should be indeces of vertices
poly = F(polyB).';

numTours = numel(tourStart)-1;
polys = cell(1,numTours);
for k=1:numTours
  polys{k} = [poly(tourStart(k):tourStart(k+1)-1),poly(tourStart(k))];
end

if numTours == 1
  poly = polys;
else
  % order such that largest Tour is the first one
  [~,order] = sort(diff(tourStart),'descend');
  poly = polys(order);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function dick_grid_display ( ng, xy )

%*****************************************************************************80
%
%% DISK_GRID_DISPLAY displays grid points inside a disk.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NG, the number of grid points inside the disk.
%
%    Input, real XY(2,NG), the grid points.
%
  scatter ( xy(1,:), xy(2,:), 'b.' );
  axis equal
  title ( sprintf ( '%d grid points inside a disk', ng ) )
  grid on

  return
end


text:
function writeInfluenceMatrixBixel(IM, structNum, planC, outFile)
%"writeInfluenceMatrixBixel"
%   Write the bixel information to binary file.
%
%   All beams in one file, for each beam:
%       float    a_g, a_t, a_c  - gantry, table and collimator angles
%       float    dx_b, dy_b     - bixel dimensions
%       float    dx, dy, dz     - voxel dimensions
%       int      Nx, Ny, Nz     - dose cube dimensions (number of voxels)
%       int      Npb            - number of pencil beams (bixels) used for this field
%       float    DoseScalefactor- conversion factor to the absolute dose, for the best resolution, this will be max(Dij_entry)/max(short) = max(Dij_entry)/(2^15-1)
%
%       For each of Npb bixels:
%           PB header:
%           float    energy         - energy
%           float    spot_x, spot_y - position of the bixel w/ respect to the central axis of the field
%           int      Nvox           - number of voxels with non-zero dose from this beamlet
%
%           For each of Nvox voxels:
%               int      VoxelNumber    - voxel ID in the cube (between 1 and Nx*Ny*Nz)
%               short    Value          - multiply this by DoseScalefactor to get the dose deposited by this beamlet to the voxel VoxelNumber, assuming the beamlet weight of 1
%
% Repeat for each beam.
%
% JRA 3/11/04
%
%Usage:
%   function writeInfluenceMatrixBixel(IM, structNum, planC, outFile)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

%Open file
if nargin < 4   %default
    outFile = 'c:\tmp\default.txt'
else
    outFile
end
fid = fopen(outFile,'w+');

%Rows of inflM are voxels, columns are PBs.
disp('Constructing influence matrix...')
[inflM] = getGlobalInfluenceM(IM, structNum);
disp('Done.')

%Find indices of voxels with nonzero dose.
indexV = find(any(inflM,2));

%Get total number of non zero voxels for each bixel.
nnzVoxelsPerBixel = full(sum(inflM ~= 0, 1));

%Get beamNumber for each beamlet.
beamN = [IM.beamlets(structNum, :).beamNum];

%Index map from influence beamlets to beam style beamlets
for i = 1:length(IM.beams)
    beamletsPerBeam([IM.beamlets(structNum, :).beamNum] == i) = find([IM.beamlets(structNum, :).beamNum] == i);
end

%For each beam...
for beamNum = 1:length(IM.beams)
	%Generate header info.
	a_g = IM.beams(beamNum).gantryAngle;
	a_t = IM.beams(beamNum).couchAngle;
    if isempty(a_t), a_t = 0; end
	a_c = IM.beams(beamNum).collimatorAngle;
    if isempty(a_c), a_c = 0; end

	dx_b = IM.beams(beamNum).beamletDelta_x;
	dy_b = IM.beams(beamNum).beamletDelta_y;

	dx = planC{indexS.scan}.uniformScanInfo.grid2Units;
	dy = planC{indexS.scan}.uniformScanInfo.grid1Units;
	dz = planC{indexS.scan}.uniformScanInfo.sliceThickness;

	siz = getUniformizedSize(planC)
	Nx = siz(1);
	Ny = siz(2);
	Nz = siz(3);

	Npb = length(find([IM.beamlets(structNum, :).beamNum] == beamNum));

    %Find bixels belonging to this beam.
    myBixels = find([IM.beamlets(structNum, :).beamNum] == beamNum);

    %Scale determined by maximum value contributed to influence matrix from all bixels belonging to this beam.
  	dScale = max(max(inflM(:, myBixels)))/(2^15-1);

	%Write header info:
	fwrite(fid, [a_g a_t a_c dx_b dy_b dx dy dz],'float32');
	fwrite(fid, [Nx Ny Nz Npb], 'int32');
    fwrite(fid, dScale, 'float32');

	%For every beamlet in this beam...
	for i = myBixels
        %Prepare data that will be written.
        bE     = IM.beams(beamNum).beamEnergy;
        spot_x = IM.beams(beamNum).xPBPosV(beamletsPerBeam(i));
        spot_y = IM.beams(beamNum).yPBPosV(beamletsPerBeam(i));
        nVox   = nnzVoxelsPerBixel(i);
        voxNum = find(inflM(:,i));
        %Convert to KonRad Cube.
        voxI = KonRadCube(voxNum, siz);
        doses = full(inflM(voxNum, i)) / dScale;

        %Write it.
        write  = fwrite(fid, bE,'float32');
        write  = fwrite(fid, spot_x,'float32');
        write  = fwrite(fid, spot_y,'float32');
        write  = fwrite(fid, nVox,'uint32');

        %Use seeking+skips to interleave values.
        pos = ftell(fid);
        fseek(fid, -2, 'cof');
        write = fwrite(fid, voxI,'int32', 2);
        fseek(fid, pos, 'bof');
        write = fwrite(fid, doses,'ushort', 4);

        if mod(i,200) == 0
            disp(['Processed ' num2str(i) ' bixels for beam ' num2str(beamNum) '.']);
        end
    end
end
clear inflM;

disp('Done.');
fclose(fid);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function F = vertcat(varargin)

F = horzcat(varargin{:});

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function pcircles(fhdl, center, radius, varargin)
%
% PCIRCLE(fhdl, center, radius, varargin)
% A simple function for drawing a circle on the selected figure using
% polylines.
%
% INPUT:
% fhdl:         figure handler
% center:       center of the circle in [x0 y0] format
% radius:       radius of the circle (in pixels)
% varargin:     colors, plot symbols and line types

format = '-w';  % default
if nargin == 4
    format = varargin{1};
end

x0 = center(1);
y0 = center(2);
nseg = ceil(2*pi*radius);
theta = 0 : (2 * pi / nseg) : (2 * pi);
x = radius * cos(theta) + x0;
y = radius * sin(theta) + y0;

figure(fhdl), hold;
plot(x0, y0, '.r', 'LineWidth',5);
plot(x, y, format);


text:
%USING_HG2 Determine if the HG2 graphics pipeline is used
%
%   tf = using_hg2(fig)
%
%IN:
%   fig - handle to the figure in question.
%
%OUT:
%   tf - boolean indicating whether the HG2 graphics pipeline is being used
%        (true) or not (false).

function tf = using_hg2(fig)
try
    tf = ~graphicsversion(fig, 'handlegraphics');
catch
    tf = false;
end

author: Cloud-CV

repo: object-proposals


text:
function set_colorbar(conf)
%SET_COLORBAR draws a color bar to the plot
%
%   Usage: set_colorbar(conf)
%
%   Input parameters:
%       conf        - configuration struct (see SFS_config)
%
%   SET_COLORBAR(conf) draws a color bar on the figure and sets the map to the
%   color specified in conf.plot.colormap.
%
%   See also: plot_sound_field, set_colormap

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameter =====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);
isargstruct(conf);

%% ===== Configuration ===================================================
p.caxis = conf.plot.caxis;
p.usedb = conf.plot.usedb;
p.colormap = conf.plot.colormap;

%% ===== Plotting ========================================================
% Change color map (default: gray)
set_colormap(p.colormap);

% Set the limits of the colormap and add a colorbar
if length(p.caxis)==2
    caxis(p.caxis);
end
if p.usedb
    h = colorbar;
    ylabel(h,'Amplitude (dB)');
    % Get the font size and name of the figure and adjust the colorbar
    fsize = get(gca,'FontSize');
    fname = get(gca,'FontName');
    set(h,'FontSize',fsize);
    set(h,'FontName',fname);
    temp = get(h,'Ylabel');
    set(temp,'FontSize',fsize);
    set(temp,'FontName',fname);
end

author: sfstoolbox

repo: sfs-matlab


text:
function lattice_rule_test07 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST07 tests LATTICE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994, page 18.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST07\n' );
  fprintf ( 1, '  LATTICE applies a lattice rule to integrate\n' );
  fprintf ( 1, '  a function over the unit hypercube.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  The lattice rule order M will vary.\n' );

  z(1:dim_num) = [ 1, 2 ];
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  i4vec_print ( dim_num, z, '  The lattice generator vector:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         I         M    EXACT       ESTIMATE    ERROR\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    m = prime ( 3 * i );

    quad = lattice ( dim_num, m, z, @f_01_2d );

    exact = e_01_2d ( dim_num, a, b );

    error = abs ( exact - quad );

    fprintf ( 1, '  %8d  %8d  %10.6f  %10.6f  %10.6e\n', i, m, exact, quad, error );

  end

  return
end


text:
function point_num = cfn_e_size ( dim_num, level_max )

%*****************************************************************************80
%
%% CFN_E_SIZE: Closed Fully Nested, Exponential Growth.
%
%  Discussion:
%
%    This calculation assumes that an exponential growth rule is being used,
%    that is, that the 1D rules have orders 1, 3, 7, 15, 31, and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the maximum value of LEVEL.
%
%    Output, integer POINT_NUM, the total number of unique
%    points in the grids.
%

%
%  Special case.
%
  if ( level_max < 0 )
    point_num = 0;
    return
  end

  if ( level_max == 0 )
    point_num = 1;
    return
  end
%
%  Construct the vector that counts the new points in the 1D rule.
%
  new_1d = zeros ( level_max+1, 1 );

  new_1d(0+1) = 1;
  new_1d(1+1) = 2;

  j = 1;
  for l = 2 : level_max
    j = j * 2;
    new_1d(l+1) = j;
  end
%
%  Count the number of points by counting the number of new points
%  associated with each level vector.
%
  level_1d = zeros ( dim_num, 1 );

  point_num = 0;

  for level = 0 : level_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );

      point_num = point_num + prod ( new_1d(level_1d(1:dim_num)+1) );

      if ( ~more )
        break
      end

    end

  end

  return
end


text:
function I = gcmi_mixture_cd(x, y, Ym)
% GCMI_MIXTURE_CD Gaussian-Copula Mutual Information between a continuous and a
%         discrete variable in bits calculated from a Gaussian mixture
%
%   The  Gaussian mixture is fit using robust measures of location (median) and
%   scale (median absolute deviation) for each class.
%   I = gcmi_mixture_cd(x,y,Ym) returns the MI between the (possibly multidimensional)
%   continuous variable x and the discrete variable y.
%   Rows of x correspond to samples, columns to dimensions/variables.
%   (Samples first axis)
%   y should contain integer values in the range [0 Ym-1] (inclusive).
%
%   See also: GCMI_MODEL_CD

% ensure samples first axis for vectors
if isvector(x)
    x = x(:);
end
if ndims(x)~=2
    error('gcmi_mixture_cd: input array should be 2d')
end

if isvector(y)
    y = y(:);
else
    error('gcmi_mixture_cd: only univariate discrete variable supported');
end

Ntrl = size(x,1);
Nvar = size(x,2);

if size(y,1) ~= Ntrl
    error('gcmi_mixture_cd: number of trials do not match');
end

% check for repeated values
for xi=1:Nvar
    if length(unique(x(:,xi)))./Ntrl < 0.9
        warning('Input x has more than 10% repeated values.')
        break
    end
end

% check values of discrete variable
if min(y)~=0 || max(y)~=(Ym-1) || any(round(y)~=y)
    error('Values of discrete variable y are not correct')
end

% copula normalise each class
% shift and rescale to match loc and scale of raw data
% this provides a robust way to fit a guassian mixture
groupdat = cell(Ym,1);
yval = cell(Ym,1);
for yi=1:Ym
    % class conditional data
    idx = y==(yi-1);
    ydat = x(idx,:);
    cydat = copnorm(ydat);
    % robust measure of s.d. under Gaussian assumption from median absolute deviation
    cyscaled = bsxfun(@times, cydat, 1.482602218505602*mad(ydat,1));
    % robust measure of loc from median
    cyscaled = bsxfun(@plus, cyscaled, median(ydat));
    groupdat{yi} = cyscaled;
    yval{yi} = (yi-1)*ones(size(ydat,1),1);
end
cx = cell2mat(groupdat);
newy = cell2mat(yval);
% Gaussian mixture MI
I = mi_mixture_gd(cx,newy,Ym);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% This library is coded for the double pendulum example
% Last Updated: 2019/07/30
% Coded By: K.Kahirman

function [Data,Sym_Struct]=SINDyLib(X,dX,iter,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order)
%% First get the size of the X matrix, determin the data length and the number of variables we have.
[Data_Length,Variable_Number]=size(X);
[~,Variable_Number_dX]=size(dX);
[~,Variable_Number_u]=size(u);
%Also create the symbolic variable
Symbol=sym('z',[Variable_Number,1]);
Symbol_dX=sym('dz',[Variable_Number,1]);
Symbol_u=sym('u',[Variable_Number_u,1]);

%Now according the Highest Polynomial Order entered, we will calculate the data matrix.
Data=[];
Index=1;

%% First calculate the polynomial term

% Form basis vector
Basis=[X(:,1) X(:,2) X(:,1)-X(:,2) X(:,1)-2*X(:,2) 2*X(:,1)-X(:,2) 2*X(:,1)-2*X(:,2)];
Basis_Sym=[Symbol(1) Symbol(2) Symbol(1)-Symbol(2) Symbol(1)-2*Symbol(2) 2*Symbol(1)-Symbol(2) 2*Symbol(1)-2*Symbol(2)];

% Add the trigonometric form
for i=1:size(Basis,2)
   Data(:,Index)=sin(Basis(:,i));
   Sym_Struct{1,Index}=sin(Basis_Sym(1,i));
   Index=Index+1;
end

for i=1:size(Basis,2)
   Data(:,Index)=cos(Basis(:,i));
   Sym_Struct{1,Index}=cos(Basis_Sym(1,i));
   Index=Index+1;
end

for i=3:size(Basis,2)
   Data(:,Index)=cos(Basis(:,i)).^2;
   Sym_Struct{1,Index}=cos(Basis_Sym(1,i))^2;
   Index=Index+1;
end

% Adding following terms will reduce the noise robustness
% for i=3:size(Basis,2)
%    Data(:,Index)=sin(Basis(:,1)).*cos(Basis(:,i));
%    Sym_Struct{1,Index}=sin(Basis_Sym(1,1))*cos(Basis_Sym(1,i));
%    Index=Index+1;
% end
%
% for i=3:size(Basis,2)
%    Data(:,Index)=cos(Basis(:,1)).*cos(Basis(:,i));
%    Sym_Struct{1,Index}=cos(Basis_Sym(1,1))*cos(Basis_Sym(1,i));
%    Index=Index+1;
% end
%
% for i=3:size(Basis,2)
%    Data(:,Index)=sin(Basis(:,2)).*cos(Basis(:,i));
%    Sym_Struct{1,Index}=sin(Basis_Sym(1,2))*cos(Basis_Sym(1,i));
%    Index=Index+1;
% end
%
% for i=3:size(Basis,2)
%    Data(:,Index)=cos(Basis(:,2)).*cos(Basis(:,i));
%    Sym_Struct{1,Index}=cos(Basis_Sym(1,2))*cos(Basis_Sym(1,i));
%    Index=Index+1;
% end

% Add polynomial term
Data(:,Index)=X(:,3);
Sym_Struct{1,Index}=Symbol(3);
Index=Index+1;

Data(:,Index)=X(:,4);
Sym_Struct{1,Index}=Symbol(4);
Index=Index+1;

for i=3:size(Basis,2)
   Data(:,Index)=X(:,3).*sin(Basis(:,i));
   Sym_Struct{1,Index}=Symbol(3)*sin(Basis_Sym(1,i));
   Index=Index+1;
end

for i=3:size(Basis,2)
   Data(:,Index)=X(:,4).*sin(Basis(:,i));
   Sym_Struct{1,Index}=Symbol(4)*sin(Basis_Sym(1,i));
   Index=Index+1;
end

for i=3:size(Basis,2)
   Data(:,Index)=X(:,3).^2.*sin(Basis(:,i));
   Sym_Struct{1,Index}=Symbol(3)^2*sin(Basis_Sym(1,i));
   Index=Index+1;
end

for i=3:size(Basis,2)
   Data(:,Index)=X(:,4).^2.*sin(Basis(:,i));
   Sym_Struct{1,Index}=Symbol(4)^2*sin(Basis_Sym(1,i));
   Index=Index+1;
end

% Add dx term
Data(:,Index)=dX(:,1);
Sym_Struct{1,Index}=Symbol_dX(iter);
Index=Index+1;

for i=3:size(Basis,2)
   Data(:,Index)=dX(:,1).*cos(Basis(:,i)).^2;
   Sym_Struct{1,Index}=Symbol_dX(iter)*cos(Basis_Sym(1,i)).^2;
   Index=Index+1;
end

% Add constant
%Order zero:
Data(:,Index)=ones(Data_Length,1);
Sym_Struct{1,Index}=1;

author: dynamicslab

repo: SINDy-PI


text:
function f = ackleysfcn(x)
% Ackley's Function.

if strcmp(x,'init')
    f.PopInitRange = [-32.768; 32.768] ;
    f.KnownMin = [0,0] ; % For plotting only
else
    x = reshape(x,1,[]) ;
    a = 20 ;
    b = 0.2 ;
    n = size(x,2) ;
    c = 2*pi ;
    f = -a*exp(-b*sqrt(1/n*(x*x'))) - exp(1/n*(cos(c*x)*cos(c*x'))) + ...
        a + exp(1) ;
end


text:
function p2Vals=floorPow2(nums)
%%FLOORPOW2   Find the largest powers of 2 that are less than or equal
%             to the non-negative, real values in vals.
%
%INPUTS:  nums  An array or matrix of non-negative, real values.
%
%OUTPUTS:p2Vals Values corresponding to the power of 2 less than or
%               equal to the values in nums having the same datatype as
%               nums (float, int32, etc.).
%
%This function checks for the data type of the passed value. If the
%value passed returns true from the function isfloat then an algorithm
%that manipulates the mantissa and exponent of a floating point number is
%used. If the data type is not a float, then it is assumed to be an integer
%and binary shift operations are used to find the next highest power of
%two.
%
%The implementation here could be easily translated in C, since the
%floating point operations correspond to functions in math.h and the
%bitshift function for integer values just implements the << and >>
%operators in C.
%
%April 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

assert(all(nums>=0))

if(isfloat(nums))
    %log2 with two return values is the frexp function in C.
    [F,E]=log2(nums);

    %pow2 with two arguments is the ldexp function in C.
    p2Vals=pow2(0.5*(F>0),E);
else
    p2Vals=zeros(size(nums),'like',nums);
    numVals=size(nums(:));

    for curVal=1:numVals
        val=ones(1,1,'like',nums);
        temp=nums(curVal);

        if(temp==0)
            p2Vals(curVal)=0;
            continue;
        end

        while(temp~=1)
            temp=bitshift(temp,-1);
            val=bitshift(val,1);
        end

        p2Vals(curVal)=val;
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function x = r8lt_sl ( n, a, b, job )

%*****************************************************************************80
%
%% R8LT_SL solves a R8LT system.
%
%  Discussion:
%
%    The R8LT storage format is used for an M by N lower triangular matrix,
%    and sets aside storage even for the entries that must be zero.
%
%    No factorization of the lower triangular matrix is required.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the R8LT matrix.
%
%    Input, real B(N), the right hand side.
%    On output, the solution vector.
%
%    Input, integer JOB, is 0 to solve the untransposed system,
%    nonzero to solve the transposed system.
%
%    Output, real X(N), the solution vector.
%
  x(1:n) = b(1:n);

  if ( job == 0 )

    for j = 1 : n
      x(j) = x(j) / a(j,j);
      x(j+1:n) = x(j+1:n) - a(j+1:n,j)' * x(j);
    end

  else

    for j = n : -1 : 1
      x(j) = x(j) / a(j,j);
      x(1:j-1) = x(1:j-1) - a(j,1:j-1) * x(j);
    end

  end

  return
end


text:
function square_arbq_rule_test03 ( degree, n, header )

%*****************************************************************************80
%
%% SQUARE_ARBQ_RULE_TEST03 gets a rule and creates GNUPLOT input files.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree exactness
%    of the quadrature rule.  0 <= DEGREE <= 50.
%
%    Input, integer N, the number of nodes to be used by the rule.
%
%    Input, string HEADER, an identifier for the filenames.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_ARBQ_RULE_TEST03\n' );
  fprintf ( 1, '  Get a quadrature rule for the symmetric square.\n' );
  fprintf ( 1, '  Set up GNUPLOT graphics input.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );
%
%  Retrieve a symmetric quadrature rule.
%
  [ x, w ] = square_arbq ( degree, n );
%
%  Create files for input to GNUPLOT.
%
  square_arbq_gnuplot ( n, x, header );

  return
end


text:
function PlotImu6(aedat, numBins, startTime, endTime)

%{
%}

if isfield(aedat.data, 'imu6')
    data = aedat.data.imu6;
else
    return
end

if ~exist('numBins', 'var') || (exist('numBins', 'var') && numBins == 0)
	numBins = data.numEvents;
end

if exist('startTime', 'var') && startTime > 0
    startTime = startTime * 1e6;
else
    startTime = aedat.info.firstTimeStamp;
end

if exist('endTime', 'var') && endTime > 0
    endTime = endTime * 1e6;
else
    endTime = aedat.info.lastTimeStamp;
end

durationUs = double(endTime - startTime);
durationOfBinUs = durationUs / numBins;
durationOfBinS = durationOfBinUs / 1000000;

timeBinBoundariesUs = double(startTime) : durationOfBinUs : double(endTime);
timeBinCentresS = (timeBinBoundariesUs(1 : end - 1) + durationOfBinUs / 2) / 1000000;

if numBins == data.numEvents
    accelX = data.accelX;
    accelY = data.accelY;
    accelZ = data.accelZ;
    gyroX  = data.gyroX;
    gyroY  = data.gyroY;
    gyroZ  = data.gyroZ;
    temperature   = data.temperature;
else
    accelX = zeros(numBins, 1);
    accelY = zeros(numBins, 1);
    accelZ = zeros(numBins, 1);
    gyroX  = zeros(numBins, 1);
    gyroY  = zeros(numBins, 1);
    gyroZ  = zeros(numBins, 1);
    temperature = zeros(numBins, 1);

    for bin = 1 : numBins
		firstTimeStampIndex = find(data.timeStamp >= timeBinBoundariesUs(bin), 1, 'first');
		lastTimeStampIndex = max(firstTimeStampIndex, find(data.timeStamp < timeBinBoundariesUs(bin + 1), 1, 'last'));
		if ~isempty(firstTimeStampIndex) && ~isempty(lastTimeStampIndex)
			accelX(bin) = mean(data.accelX(firstTimeStampIndex : lastTimeStampIndex));
			accelY(bin) = mean(data.accelZ(firstTimeStampIndex : lastTimeStampIndex));
			accelZ(bin) = mean(data.accelY(firstTimeStampIndex : lastTimeStampIndex));
			gyroX(bin) = mean(data.gyroX(firstTimeStampIndex : lastTimeStampIndex));
			gyroY(bin) = mean(data.gyroY(firstTimeStampIndex : lastTimeStampIndex));
			gyroZ(bin) = mean(data.gyroZ(firstTimeStampIndex : lastTimeStampIndex));
			temperature(bin) = mean(data.temperature(firstTimeStampIndex : lastTimeStampIndex));
		end
    end
end

figure
legendLocal = {};
hold all
plot(timeBinCentresS, accelX, '-')
legendLocal = [legendLocal 'accelX'];
plot(timeBinCentresS, accelY, '-')
legendLocal = [legendLocal 'accelY'];
plot(timeBinCentresS, accelZ, '-')
legendLocal = [legendLocal 'accelZ'];
xlabel('Time (s)')
ylabel('Acceleration (g)')
legend(legendLocal)

figure
legendLocal = {};
hold all
plot(timeBinCentresS, gyroX)
legendLocal = [legendLocal 'gyroX'];
plot(timeBinCentresS, gyroY)
legendLocal = [legendLocal 'gyroY'];
plot(timeBinCentresS, gyroZ)
legendLocal = [legendLocal 'gyroZ'];
xlabel('Time (s)')
ylabel('Angular velocity (deg/s)')
legend(legendLocal)

figure
plot(timeBinCentresS, temperature)
xlabel('Time (s)')
ylabel('temperature (C)')

author: panpanfei


text:
% Lightspeed Toolbox.
% Efficient operations for Matlab programming.
% Version 2.8   04-Jan-2017
% By Tom Minka
%
% Matrix algebra
%   repmat           - Fast replacement for matlab's repmat.
%   xrepmat          - Matlab's original repmat.
%   setnonzeros      - Fast creation of sparse matrix.
%   row_sum          - Sum for each row.  Faster than 'sum'.
%   scale_rows       - Scale each row of a matrix.
%   scale_cols       - Scale each column of a matrix.
%   solve_triu       - Left division by upper triangular matrix.
%   solve_tril       - Left division by lower triangular matrix.
%   sqdist           - Squared Euclidean and Mahalanobis distance.
%   isposdef         - Check for positive-definiteness.
%   logdet           - log(determinant) for positive definite matrix.
%   cholproj         - Projected Cholesky factorization.
%   inv_posdef       - Invert positive definite matrix.
%
% Statistics
%   mvnormpdf        - Multivariate normal density.
%   mvnormpdfln      - Log of multivariate normal density.
%   normcdf          - Normal cumulative distribution.
%   normcdfln        - Log of normal cumulative distribution.
%   normcdflogit     - Logit of normal cumulative distribution.
%   invnormcdf       - Normal quantile function.
%   wishpdf          - Wishart probability density function.
%   wishpdfln        - Log of Wishart probability density function.
%   sample           - Sample from categorical distribution.
%   sample_vector    - Sample from multiple categorical distributions.
%   sample_hist      - Sample from multinomial distribution.
%   randbinom        - Sample from binomial distribution.
%   randnorm         - Sample from multivariate normal.
%   randgamma        - Sample from Gamma distribution.
%   randbeta         - Sample from Beta distribution.
%   randwishart      - Sample from Wishart distribution.
%   randomseed       - Get or set the random seed.
%   int_hist         - Histogram of integer values.
%
% Utility
%   logsumexp        - Sum in the log domain.
%   logmulexp        - Matrix multiply in the log domain.
%   ndgridmat        - Matrix of grid points.
%   ind2subv         - Subscript vector from linear index.
%   subv2ind         - Linear index from subscript vector.
%   gammaln          - Fast replacement for matlab's gammaln.
%   digamma          - Derivative of gammaln.
%   trigamma         - Derivative of digamma.
%   tetragamma       - Derivative of trigamma.
%   ndsum            - Sum over multiple dimensions.
%   ndmax            - Maximum over multiple dimensions.
%   ndlogsumexp      - Sum over multiple dimensions in the log domain.
%   maxdiff          - Maximum difference between structs or arrays.
%   sameobject       - Test if two variables correspond to the same object.
%   find_sameobject  - Find an object in a cell array.
%   toJava           - Convert to Java representation.
%   fromJava         - Convert from Java to Matlab.
%   glob             - Filename expansion via wildcards.
%   globstrings      - String matching via wildcards.
%
% Argument lists
%   argfilter        - Remove unwanted arguments from a key/value list.
%   makestruct       - Cell-friendly alternative to STRUCT.
%   setfields        - Set multiple fields of a structure.
%   struct2arglist   - Convert structure to cell array of fields/values.
%
% Mutation
%   mutable          - Convert to a mutable object.
%   immutable        - Convert to an ordinary (immutable) object.
%
% Set operations
%   ismember_sorted  - True for member of sorted set.
%   match            - Location of matches in a set.
%   match_sorted     - Location of matches in a sorted set.
%   setdiff_sorted   - Set difference between sorted sets.
%   intersect_sorted - Set intersection between sorted sets.
%   union_sorted     - Set union of sorted sets.
%   union_sorted_rows - Set union of sorted sets of row vectors.
%   duplicated       - Find duplicated rows in a matrix.
%
% Readability
%   rows             - Number of rows.
%   cols             - Number of columns.
%   col_sum          - Sum for each column.
%   setdiag          - Modify the diagonals of a matrix.
%   finddiag         - Index of elements on diagonals.
%   argmax           - Index of maximum element.
%   argmin           - Index of minimum element.
%
% Flop counting
%   flops            - Read/write flop counter.
%   addflops         - Add to flop counter.
%   flops_chol       - Flops for Cholesky decomposition.
%   flops_col_sum    - Flops for column sums.
%   flops_det        - Flops for matrix determinant.
%   flops_digamma    - Flops for gammaln, digamma, and trigamma.
%   flops_div        - Flops for division.
%   flops_exp        - Flops for exponential.
%   flops_inv        - Flops for matrix inversion.
%   flops_mul        - Flops for real matrix multiplication.
%   flops_normpdfln  - Flops for normpdfln.
%   flops_pow        - Flops for raising to real power.
%   flops_randnorm   - Flops for randnorm.
%   flops_row_sum    - Flops for row sums.
%   flops_sample     - Flops for sample(p,n).
%   flops_solve      - Flops for matrix left division.
%   flops_solve_tri  - Flops for triangular left division.
%   flops_spadd      - Flops for sparse matrix addition.
%   flops_spmul      - Flops for sparse matrix multiplication.
%   flops_sqrt       - Flops for square root.
%
% Stand alone programs
%   matfile          - Read/write MAT files.
%   tests/test_flops - Compare time versus flops for various math operations.
%
% Graphics utilities
%  see graphics/Contents.m
%
% Demos
%   tests/test_repmat,
%   tests/test_solve_tri, ...

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
function [l, L_rf, L_sf, L_obs, L_n, N] = retroProjLmk(Rob,Sen,Obs,Opt)

% RETROPROJLMK  Retro project landmark.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

switch Sen.type

    % camera pinHole
    case {'pinHole'}
        % type of lmk to init
        switch Opt.init.initType
            case {'idpPnt'}
                % INIT LMK OF TYPE: Inverse depth point
                [l, L_rf, L_sf, L_k, L_c, L_obs, L_n] = ...
                    retroProjIdpPntFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y, ...
                    Opt.init.idpPnt.nonObsMean) ;

                N = Opt.init.idpPnt.nonObsStd^2 ;

            case {'hmgPnt'}
                % INIT LMK OF TYPE: Homogeneous point
                [l, L_rf, L_sf, L_k, L_c, L_obs, L_n] = ...
                    retroProjHmgPntFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y, ...
                    Opt.init.idpPnt.nonObsMean) ;

                N = Opt.init.idpPnt.nonObsStd^2 ;

            case {'ahmPnt'}
                % INIT LMK OF TYPE: Anchored Homogeneous point
                [l, L_rf, L_sf, L_k, L_c, L_obs, L_n] = ...
                    retroProjAhmPntFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y, ...
                    Opt.init.idpPnt.nonObsMean) ;

                N = Opt.init.idpPnt.nonObsStd^2 ;

            case {'fhmPnt'}
                % INIT LMK OF TYPE: Framed Homogeneous point
                [l, L_rf, L_sf, L_k, L_c, L_obs, L_n] = ...
                    retroProjFhmPntFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y, ...
                    Opt.init.idpPnt.nonObsMean) ;

                N = Opt.init.idpPnt.nonObsStd^2 ;

            case {'plkLin'}
                % INIT LMK OF TYPE: Plucker line
                [hm, HM_obs] = seg2hmgLin(Obs.meas.y);
                [l, L_rf, L_sf, L_k, L_hm, L_n] = ...
                    retroProjPlkLinFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    hm, ...
                    Opt.init.plkLin.nonObsMean) ;
                L_obs = L_hm*HM_obs;

                N = diag(Opt.init.plkLin.nonObsStd.^2) ;

            case 'aplLin'
                % INIT LMK OF TYPE: Anchored Plucker line
                [hm, HM_obs] = seg2hmgLin(Obs.meas.y);
                [l, L_rf, L_sf, L_k, L_hm, L_n] = ...
                    retroProjAplLinFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    hm, ...
                    Opt.init.plkLin.nonObsMean) ;
                L_obs = L_hm*HM_obs;

                N = diag(Opt.init.plkLin.nonObsStd.^2) ;

            case 'idpLin'
                % INIT LMK TYPE: Inverse-depth line.
                nMean = Opt.init.idpPnt.nonObsMean*[1;1]; % non-measured prior
                nStd  = Opt.init.idpPnt.nonObsStd*[1;1];

                [l, L_rf, L_sf, L_k, L_obs, L_n] = ...
                    retroProjIdpLinFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Obs.meas.y, ...
                    nMean) ;

                N = diag(nStd.^2);

            case 'hmgLin'
                % INIT LMK TYPE: Homogenous points line.
                nMean = Opt.init.idpPnt.nonObsMean*[1;1]; % non-measured prior
                nStd  = Opt.init.idpPnt.nonObsStd*[1;1];

                [l, L_rf, L_sf, L_k, L_obs, L_n] = ...
                    retroProjHmgLinFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Obs.meas.y, ...
                    nMean) ;

                N = diag(nStd.^2);

            case 'ahmLin'
                % INIT LMK TYPE: Anchored Homogenous points line.
                nMean = Opt.init.idpPnt.nonObsMean*[1;1]; % non-measured prior
                nStd  = Opt.init.idpPnt.nonObsStd*[1;1];

                [l, L_rf, L_sf, L_k, L_obs, L_n] = ...
                    retroProjAhmLinFromPinHoleOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Obs.meas.y, ...
                    nMean) ;

                N = diag(nStd.^2);

            otherwise
                error('??? Unknown landmark type ''%s'' for initialization.',Opt.init.initType)
        end

        % Pin hole camera with depth information
    case 'pinHoleDepth'

        switch Opt.init.initType
            case 'eucPnt'
                % INIT LMK OF TYPE: Euclidean point
                [l, L_rf, L_sf, L_k, L_c, L_obs] = ...
                    retroProjEucPntFromPhdOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y) ;

                L_n = zeros(3,0);
                N   = [] ;

            otherwise
                error('??? Unknown landmark type ''%s'' for initialization.',Opt.init.initType)
        end

    % Omnidirectional camera
    case {'omniCam'}
        % type of lmk to init
        switch Opt.init.initType
            case {'ahmPnt'}
                % INIT LMK OF TYPE: Anchored Homogeneous point
                [l, L_rf, L_sf, L_k, L_c, L_obs, L_n] = ...
                    retroProjAhmPntFromOmniCamOnRob( ...
                    Rob.frame, ...
                    Sen.frame, ...
                    Sen.par.k, ...
                    Sen.par.c, ...
                    Obs.meas.y, ...
                    Opt.init.idpPnt.nonObsMean) ;

                N = Opt.init.idpPnt.nonObsStd^2 ;
            otherwise
                error('??? Unknown landmark type ''%s'' for initialization.',Opt.init.initType)
        end

    otherwise % -- Sen.type
        % Print an error and exit
        error('??? Unknown sensor type ''%s''.',Sen.type);
end % -- Sen.type

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb