text:
function varargout = process_segment_brainsuite( varargin )
% PROCESS_SEGMENT_BRAINSUITE: Run the segmentation of a T1 MRI with BrainSuite.
%
% USAGE:     OutputFiles = process_segment_brainsuite('Run',     sProcess, sInputs)
%         [isOk, errMsg] = process_segment_brainsuite('Compute', iSubject, iAnatomy=[default], nVertices, isInteractive)
%                          process_segment_brainsuite('ComputeInteractive', iSubject, iAnatomy)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Segment MRI with BrainSuite';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import anatomy'};
    sProcess.Index       = 32;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/SegBrainSuite';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'import'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    % Title
    sProcess.options.doc.Comment = ['BrainSuite must be installed on the computer,
' ...
                                    'the command "bse" must be available in the system path.
' ...
                                    'Website: http://brainsuite.org

'];
    sProcess.options.doc.Type    = 'label';
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = '';
    % Option: Number of vertices
    sProcess.options.nvertices.Comment = 'Number of vertices (cortex): ';
    sProcess.options.nvertices.Type    = 'value';
    sProcess.options.nvertices.Value   = {15000, '', 0};
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % Number of vertices
    nVertices = sProcess.options.nvertices.Value{1};
    if isempty(nVertices) || (nVertices < 50)
        bst_report('Error', sProcess, [], 'Invalid number of vertices.');
        return
    end
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return;
    end
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    if isempty(iSubject)
        bst_report('Error', sProcess, [], ['Subject "' SubjectName '" does not exist.']);
        return
    end
    % Call processing function
    [isOk, errMsg] = Compute(iSubject, [], nVertices, 0);
    % Handling errors
    if ~isOk
        bst_report('Error', sProcess, [], errMsg);
    elseif ~isempty(errMsg)
        bst_report('Warning', sProcess, [], errMsg);
    end
    % Return an empty structure
    OutputFiles = {'import'};
end

%% ===== COMPUTE BRAINSUITE SEGMENTATION =====
function [isOk, errMsg] = Compute(iSubject, iAnatomy, nVertices, isInteractive)
    errMsg = '';
    isOk = 0;

    % ===== INSTALL BRAINSUITE =====
    bst_progress('text', 'Testing BrainSuite installation...');
    % Get BrainSuite path from Brainstorm preferences
    BsDir = bst_get('BrainSuiteDir');
    % Get executable names
    if ispc
        cse_exe = fullfile(BsDir, 'bin', 'cortical_extraction.cmd');
        svreg_exe = fullfile(BsDir, 'svreg', 'bin', 'svreg.exe');
    else
        cse_exe = fullfile(BsDir, 'bin', 'cortical_extraction.sh');
        svreg_exe = fullfile(BsDir, 'svreg', 'bin', 'svreg.sh');
    end
    % Check BrainSuite installation
    if isempty(BsDir) || ~file_exist(cse_exe)
        errMsg = ['BrainSuite is not installed on your computer.' 10 ...
                      'Download it from http://brainsuite.org and install it.' 10 ...
                      'Then set its installation folder in the Brainstorm options (File > Edit preferences)'];
        return
    end

    % ===== GET SUBJECT =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', iSubject);
    if isempty(sSubject)
        errMsg = 'Subject does not exist.';
        return
    end
    % Check if a MRI is available for the subject
    if isempty(sSubject.Anatomy)
        errMsg = ['No MRI available for subject "' SubjectName '".'];
        return
    end
    % Get default MRI if not specified
    if isempty(iAnatomy)
        iAnatomy = sSubject.iAnatomy;
    end

    % ===== DELETE EXISTING SURFACES =====
    % Confirm with user that the existing surfaces will be removed
    if isInteractive && ~isempty(sSubject.Surface)
        isDel = java_dialog('confirm', ['Warning: There are already surfaces in this subject.' 10 ...
            'Running BrainSuite will remove all the existing surfaces.' 10 10 ...
            'Delete the existing files?' 10 10], 'BrainSuite segmentation');
        if ~isDel
            errMsg = 'Process aborted by user.';
            return;
        end
    end

    % ===== VERIFY FIDUCIALS IN MRI =====
    % Load MRI file
    T1FileBst = sSubject.Anatomy(iAnatomy).FileName;
    sMri = in_mri_bst(T1FileBst);
    % If the SCS transformation is not defined: compute MNI transformation to get a default one
    if isempty(sMri) || ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'NAS') || ~isfield(sMri.SCS, 'LPA') || ~isfield(sMri.SCS, 'RPA') || (length(sMri.SCS.NAS)~=3) || (length(sMri.SCS.LPA)~=3) || (length(sMri.SCS.RPA)~=3) || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS, 'T') || isempty(sMri.SCS.T)
        % Issue warning
        errMsg = 'Missing NAS/LPA/RPA: Computing the MNI normalization to get default positions.';
        % Compute MNI normalization
        [sMri, errNorm] = bst_normalize_mni(T1FileBst);
        % Handle errors
        if ~isempty(errNorm)
            errMsg = [errMsg 10 'Error trying to compute the MNI normalization: ' 10 errNorm 10 ...
                'Missing fiducials: the surfaces cannot be aligned with the MRI.'];
        end
    end

    % ===== SAVE MRI AS NII =====
    bst_progress('text', 'Saving temporary files...');
    % Create temporay folder for BrainSuite output
    TmpDir = bst_get('BrainstormTmpDir', 0, 'brainsuite');
    % Save MRI in .nii format
    subjid = strrep(sSubject.Name, '@', '');
    NiiFile = bst_fullfile(TmpDir, [subjid, '.nii']);
    out_mri_nii(sMri, NiiFile);
    % If a "world transformation" was not available in the MRI in the database, it was set to a default when saving to .nii
    % Let's reload this file to get the transformation matrix, it will be used when importing the results
    if ~isfield(sMri, 'InitTransf') || isempty(sMri.InitTransf) || isempty(find(strcmpi(sMri.InitTransf(:,1), 'vox2ras')))
        % Load again the file, with the default vox2ras transformation
        [tmp, vox2ras] = in_mri_nii(NiiFile, 0, 0, 0);
        % Prepare the history of transformations
        if ~isfield(sMri, 'InitTransf') || isempty(sMri.InitTransf)
            sMri.InitTransf = cell(0,2);
        end
        % Add this transformation in the MRI
        sMri.InitTransf(end+1,[1 2]) = {'vox2ras', vox2ras};
        % Save modification on hard drive
        bst_save(file_fullpath(T1FileBst), sMri, 'v7');
    end

    % ===== 1. CORTICATL EXTRACTION =====
    bst_progress('text', '1/2: Cortical extraction...    (see command window)');
    strCall = ['"' cse_exe '" "' NiiFile(1:end-4) '"'];
    disp(['BST> System call: ' strCall]);
    status = system(strCall)
    % Error handling
    if (status ~= 0)
        errMsg = ['BrainSuite failed at step 1/2 (Cortical extraction).', 10, 'Check the Matlab command window for more information.'];
        return
    end

    % ===== 2. SVREG =====
    bst_progress('text', '2/2: SVREG...    (see command window)');
    AtlasPath = fullfile(BsDir, 'svreg', 'BrainSuiteAtlas1', 'mri');
    strCall = ['"' svreg_exe '" "' NiiFile(1:end-4) '" "' AtlasPath '"'];
    disp(['BST> System call: ' strCall]);
    status = system(strCall)
    % Error handling
    if (status ~= 0)
        errMsg = ['BrainSuite failed at step 2/2 (SVREG).', 10, 'Check the Matlab command window for more information.'];
        return
    end

    % ===== 3. SKULL EXTRACTION =====
    % cmd=[skullfinder_exe,' -i ', subbasename,'.nii.gz -o ',subbasename,'.skull.label.nii.gz -m ',subbasename,'.mask.nii.gz --scalplabel ',subbasename,'.scalp.label.nii.gz -s ',subbasename];
    % '$HOME/BrainSuite21a/bin/skullfinder' -i subject.nii.gz -o subject.skull.label.nii.gz -m subject.mask.nii.gz -s subject
    %
    % Anand 28-Oct-2022: There is a caveat. It seems that due to a bug in the executable, outer skull and scalp surfaces are swapped.
    % For now, we can write a code that swaps subject.scalp.dfs and subject.outer_skull.dfs

    % ===== IMPORT OUTPUT FOLDER =====
    % Import BrainSuite anatomy folder
    isKeepMri = 1;
    isVolumeAtlas = 1;
    errMsg = import_anatomy_bs(iSubject, TmpDir, nVertices, 0, [], isVolumeAtlas, isKeepMri);
    if ~isempty(errMsg)
        return;
    end
    % Delete temporary folder
    file_delete(TmpDir, 1, 1);
    % Return success
    isOk = 1;
end

%% ===== COMPUTE/INTERACTIVE =====
function ComputeInteractive(iSubject, iAnatomy) %#ok<DEFNU>
    % Get inputs
    if (nargin < 2) || isempty(iAnatomy)
        iAnatomy = [];
    end
    % Ask for number of vertices
    nVertices = java_dialog('input', 'Number of vertices on the cortex surface:', 'BrainSuite segmentation', [], '15000');
    if isempty(nVertices)
        return
    end
    nVertices = str2double(nVertices);
    % Open progress bar
    bst_progress('start', 'BrainSuite', 'BrainSuite MRI segmentation...');
    % Run BrainSuite
    isInteractive = 1;
    [isOk, errMsg] = Compute(iSubject, iAnatomy, nVertices, isInteractive);
    % Error handling
    if ~isOk
        bst_error(errMsg, 'BrainSuite MRI segmentation', 0);
    elseif ~isempty(errMsg)
        java_dialog('msgbox', ['Warning: ' errMsg]);
    end
    % Close progress bar
    bst_progress('stop');
end

author: brainstorm-tools


text:
function y = ga44(x)
    y = 1 - 0.1 * (sin(x(1)^2 + x(2)) - 0.1) / (x(1)^2 + x(2)^2);
end

author: LiYangSir

repo: Smart-Algorithm


text:
function boundary_segment_num = p06_boundary_segment_num ( )

%*****************************************************************************80
%
%% P06_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
%
  boundary_segment_num = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [sig1,sig2]=spreader(Isymbols,Qsymbols,seq1,seq2,N)
%Bob Gess, June 2008 for EE473 (Header data  and time vector from QPSK_TX_IQ_RX program
%written by JC and posted to Mathworks Dec 2005)

%This module spreads the data by the PRN sequence generated previously

%The plots are scaled by the data rate.  If you have a high ratio between
%the chip rate and the data rate, some of the plots become unintelligible.

%N=1e3;		        % Number of data bits(bit rate)
fs=40*2e3;		    % Sampling frequency
Fn=fs/2;            % Nyquist frequency
Ts=1/fs;	        % Sampling time = 1/fs
T=1/N;		        % Bit time
randn('state',0);   % Keeps PRBS from changing on reruns
td=[0:Ts:(N*T)-Ts]';% Time vector(data)(transpose)

tiq = [0:Ts*2:(N*T)-Ts]';% Time vector for I and Q symbols(transpose)

%Spreading data by prn sequences
sig1=Isymbols.*seq1;
sig2=Qsymbols.*seq2;

%=====================================================================
%Plots
%======================================================================
figure(2)

subplot(2,1,1)
plot(tiq,sig1)
axis([0 5/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('I Channel(one bit/symbol(phase)) Data')

subplot(2,1,2)
plot(tiq,sig2)
axis([0 5/N -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Q Channel(one bit/symbol(phase)) Data')

end

author: Sable

repo: mcbench-benchmarks


text:
function openAxisInFig(hAxis)
%function openAxisInFig(hAxis)
%
%This function opens the passed hAxis in a new figure.
%hAxis must be one the axes stored under stateS.handle.CERRAxis
%
%APA, 11/17/2006
%
% 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/>.

global stateS
doseSet = getAxisInfo(hAxis,'doseSets');

f = figure;

%copy the specified CERR axis (2 in this example) to new figure
h1 = copyobj(hAxis,f);
set(h1,'units','normalized','Position',[0.01 0.01 0.9 0.98])
axis(h1,'square')
h1Color = get(h1,'color');
h1Children = get(h1,'children');
hLabelV = findobj(h1Children,'tag','CERRAxisLabel');
if sum(h1Color)==0
    set(hLabelV,'color',[1 1 1])
else
    set(hLabelV,'color',[0 0 0])
end

%copy colorbar to new figure
if stateS.doseToggle==1 && ~isempty(doseSet)
    h2 = copyobj([stateS.handle.doseColorbar.trans],f);
    set(h2,'units','normalized','Position',[0.9 0.15 0.08 0.8])
    axis(h2,'auto')
else
    set(h1,'units','normalized','Position',[0.01 0.01 0.95 0.98])
    colormap('gray')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
pkg load image
files = dir('tmp/all/*.jpg');
for d=1:length(files)
file_name = fullfile('tmp/all', files(d).name);
fprintf('%s\n', file_name);
I = imread(file_name);
max_wh = max(size(I));
if max_wh > 1024
    I = imresize(I, 1024/max_wh);
end

b = pencil_draw(I);
R = colour_pencil_draw(I);
imwrite(b, strcat('pencil_',files(d).name), 'jpg');
imwrite(R, strcat('colour_',files(d).name), 'jpg');
end

author: fumin

repo: pencil

save_path: github-repos/MATLAB/fumin-pencil


text:
function G=gauss2dmultislice(sizevec, m, s, a)
% G=gauss2dmultislice(sizevec, m, s, a)
lm=length(m);
if length(s)<lm
    s=repmat(s,lm); %all the sama sigma
end
if length(a)<lm
    a=repmat(a,lm); %all the sama amplitude
end

G=zeros(sizevec);
for ii=1:sizevec(3)
    G(:,:,ii) = gauss2d(sizevec(1:2), m(ii,:), s(ii), a(ii));
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function hm=plotsurf(node,face,varargin)
%
% hm=plotsurf(node,face,opt)
%
% plot 3D surface meshes
%
% author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
%
% input:
%      node: node coordinates, dimension (nn,3); if node has a
%            4th column, it will be used to set the color at each node.
%      face: triangular surface face list; if face has a 4th column,
%            it will be used to separate the surface into
%            sub-surfaces and display them in different colors;
%            face can be a cell array, each element of the array represents
%            a polyhedral facet of the mesh, if an element is an array with
%            two array subelements, the first one is the node index, the
%            second one is a scalar as the group id of the facet.
%      opt:  additional options for the plotting, see plotmesh
%
% output:
%   hm: handle or handles (vector) to the plotted surfaces
%
% example:
%
%   h=plotsurf(node,face);
%   h=plotsurf(node,face,'facecolor','r');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

rngstate = rand ('state');

if(nargin>=2)
  randseed=hex2dec('623F9A9E'); % "U+623F U+9A9E"
  if(isoctavemesh) randseed=randseed+3; end
  if(~isempty(getvarfrom({'caller','base'},'ISO2MESH_RANDSEED')))
        randseed=getvarfrom({'caller','base'},'ISO2MESH_RANDSEED');
  end
  rand('state',randseed);

  if(iscell(face))
    sc=sparse(10,3); % face colormap
    sc(1:10,:)=rand(3,10)';
    len=length(face);
    newsurf=cell(1);
    % reorganizing each labeled surface into a new cell
    for i=1:len
        fc=face{i};
        if(iscell(fc) && length(fc)>=2)
            if(fc{2}+1>10)
                sc(fc{2}+1,:)=rand(1,3);
            end
            if(fc{2}+1>length(newsurf))
                newsurf{fc{2}+1}={};
            end
            newsurf{fc{2}+1}{end+1}=fc{1};
        else % unlabeled facet is tagged by 0
            if(iscell(fc))
                newsurf{1}{end+1}=cell2mat(fc);
            else
                newsurf{1}{end+1}=fc;
            end
        end
    end
    hold on;
	h=[];
    newlen=length(newsurf);

    for i=1:newlen
        if(isempty(newsurf{i})); continue; end
        try
            subface=cell2mat(newsurf{i}')';
            if(size(subface,1)>1 && ndims(subface)==2)
               subface=subface';
            end
            h=[h patch('Vertices',node,'Faces',subface,'facecolor',sc(i,:),varargin{:})];
        catch
            for j=1:length(newsurf{i})
                h=[h patch('Vertices',node,'Faces',newsurf{i}{j},'facecolor',sc(i,:),varargin{:})];
            end
        end
    end
  else
    if(size(face,2)==4)
        tag=face(:,4);
		types=unique(tag);
        hold on;
		h=[];
	    for i=1:length(types)
            if(size(node,2)==3)
                h=[h plotasurf(node,face(find(tag==types(i)),1:3),'facecolor',rand(3,1),varargin{:})];
            else
                h=[h plotasurf(node,face(find(tag==types(i)),1:3),varargin{:})];
            end
        end
    else
        h=plotasurf(node,face,varargin{:});
    end
  end
end
if(~isempty(h))
  axis equal;
  if(all(get(gca,'view')==[0 90]))
      view(3);
  end
end
if(~isempty(h) && nargout>=1)
  hm=h;
end

rand ('state',rngstate);

%-------------------------------------------------------------------------
function hh=plotasurf(node,face,varargin)
isoct=isoctavemesh;
if(size(node,2)==4)
	if(isoct && ~exist('trisurf','file'))
	    h=trimesh(face(:,1:3),node(:,1),node(:,2),node(:,3),node(:,4),'edgecolor','k',varargin{:});
	else
	    h=trisurf(face(:,1:3),node(:,1),node(:,2),node(:,3),node(:,4),varargin{:});
	end
else
	if(isoct && ~exist('trisurf','file'))
	    h=trimesh(face(:,1:3),node(:,1),node(:,2),node(:,3),'edgecolor','k',varargin{:});
	else
	    h=trisurf(face(:,1:3),node(:,1),node(:,2),node(:,3),varargin{:});
	end
end
if(exist('h','var')) hh=h; end

repo: REST


text:
function Ainv = incremental_invert(Pinv,Q,S)
%
%   Computes the new inverse of a symmetric matrix from an old submatrix inverse and the new parts.
%
%   So if   A = [P,Q;Q',S] incremental_invert gets inv(P) , Q and S and
%   returns inv(A)
%
%

if isempty(S) && isempty(Q)
    Ainv = Pinv; return
end

Qt = Q';
M = S-Qt*Pinv*Q;
Minv = pinv(M);

Ptmp = Pinv + Pinv*Q*Minv*Qt*Pinv;
Qtmp = -Pinv*Q*Minv;

Ainv = [Ptmp,Qtmp;Qtmp',Minv];

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [xi,w]=fifthOrder2DCubPoints()
%%FIFTHORDER2DCUBPOINTS Generate fifth-order cubature points
%               for integration over a 2-dimensional cube with bounds of
%               (-1,-1), (-1,1), (1,1), and (1,-1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 2XnumCubPoints set of points for the standard
%            square.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard square (4).
%
%This function implements the points given in [1] (8 points).
%
%EXAMPLE:
%We compare a 4th-order moment computed using these cubature points
%to one computed using monomialIntCube (a 5th order moment would have just
%been 0). The results are the same within typical finite precision limits.
% [xi,w]=fifthOrder2DCubPoints();
% alpha=[2;2];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntCube(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[0.68313005106397322554806924536807013272,                                          0,    0.8163265306122448979591836734693877551;
                                        0,   0.68313005106397322554806924536807013272,    0.8163265306122448979591836734693877551;
-0.68313005106397322554806924536807013272,                                          0,    0.8163265306122448979591836734693877551;
                                        0,  -0.68313005106397322554806924536807013272,    0.8163265306122448979591836734693877551;
  0.8819171036881968635005385845464201419,    0.8819171036881968635005385845464201419,    0.1836734693877551020408163265306122449;
  0.8819171036881968635005385845464201419,   -0.8819171036881968635005385845464201419,    0.1836734693877551020408163265306122449;
 -0.8819171036881968635005385845464201419,    0.8819171036881968635005385845464201419,    0.1836734693877551020408163265306122449;
 -0.8819171036881968635005385845464201419,   -0.8819171036881968635005385845464201419,    0.1836734693877551020408163265306122449];
w=M(:,3);
xi=M(:,1:2)';

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 [afixed,sqnorm,Qahat,Z] = lambda_routine1 (afloat,Qahat,ncands,fidlog)
%LAMBDA1: Integer ambiguity estimation using LAMBDA (extended version)
%
% This routine performs an integer ambiguity estimation using the
% LAMBDA-method, as developed by the Delft University of Technology,
% Mathematical Geodesy and Positioning.
%
% This version (lambda1) is intended to be used mainly for research.
% It has more options than the original code, such as output of
% intermediate results, optional number of candidates to be returned
% and so on. If you are not interested in these options, you can use
% "lambda2"
%
% Input arguments:
%    afloat: Float ambiguities (\hat{a}, must be a column!)
%    Qahat : Variance/covariance matrix of ambiguities (Q_\hat{a})
%    ncands: Number of candidates to be returned (default 2)
%    fidlog: File-id for intermediate results (default 0, no output)
%
% Output arguments:
%    afixed: Estimated integers
%    sqnorm: Distance between candidates and float ambiguity vector
%    Qzhat : Decorrelated variance/covariance matrix
%    Z     : Transformation matrix

% ----------------------------------------------------------------------
% File.....: lambda1.m
% Date.....: 19-MAY-1999
% Version..: 2.0b
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

% ----------------------------------------------------------------------
% Initialise
% ----------------------------------------------------------------------

if nargin < 3; ncands = 2; end;
if nargin < 4; fidlog = 0; end;

n      = size (Qahat,1);
afixed = zeros(n,ncands);
sqnorm = zeros(1,ncands);

if fidlog ~= 0;
   if fidlog == 1; clc; end;
   fprintf (fidlog, '\n');
   fprintf (fidlog, ' -----------------------------------------------------------------------\n');
   fprintf (fidlog, ' "LAMBDA" - Least Squares Ambiguity Decorrelation Adjustment            \n');
   fprintf (fidlog, ' Rel. 2.0b d.d. 19 MAY 1999                                             \n');
   fprintf (fidlog, ' (c) Department of Mathematical Geodesy and Positioning (MGP)           \n');
   fprintf (fidlog, '     Faculty of Civil Engineering and Geosciences                       \n');
   fprintf (fidlog, '     Delft University of Technology, The Netherlands                    \n');
   fprintf (fidlog, ' -----------------------------------------------------------------------\n');
   fprintf (fidlog, '\n');
end;

% ----------------------------------------------------------------------
% Display the input-arguments
% ----------------------------------------------------------------------

writemat (fidlog,'Original float ambiguities',afloat);
writemat (fidlog,'Original variance/covariance matrix (Qahat)',Qahat);
writemat (fidlog,'Number of candidates requested',ncands);

% ----------------------------------------------------------------------
% --- Perform tests on the input, these tests are rather extensive   ---
% --- and can be switched off (by commenting-out) if necessary       ---
% --- Tests: Is the Q-matrix symmetric?                              ---
% ---        Is the Q-matrix positive-definite?                      ---
% ---        Do the Q-matrix and ambiguity-vector have identical     ---
% ---        dimensions?                                             ---
% ---        Is the ambiguity vector a column?                       ---
% ----------------------------------------------------------------------

if ~isequal(Qahat-Qahat'<1E-12,ones(size(Qahat)));
  error ('Variance/covariance matrix is not symmetric!');
end;

if sum(eig(Qahat)>0) ~= size(Qahat,1);
  error ('Variance/covariance matrix is not positive definite!');
end;

if length(afloat) ~= size(Qahat,1);
  error (['Variance/covariance matrix and vector of ambiguities do not have' ...
	  ' identical dimensions!']);
end;

if size(afloat,2) ~= 1;
  error ('Ambiguity-vector should be a column-vector');
end;

% ----------------------------------------------------------------------
% Make estimates in 'a' between -1 and +1 by subtracting an
% integer number, store the increments in 'incr' (= shift the centre
% of the ellipsoid over the grid by an integer translation)
% ----------------------------------------------------------------------

incr   = afloat - rem(afloat,1);
afloat = rem(afloat,1);
writemat (fidlog,'shifted ambiguities - increments',[afloat incr]);

% ----------------------------------------------------------------------
% Compute the Z-transformation based on L and D of Q, ambiguities
% are transformed according to \hat{z} = Z^T\hat{a}, Q is transformed
% according to Q_\hat{z} = Z^T * Q_\hat{a} * Z
% ----------------------------------------------------------------------

[Qz,Z,L,D,z] = decorrel (Qahat,afloat);

writemat (fidlog,'Decorrelated Q-matrix',Qz);
writemat (fidlog,'Transformation matrix (Z)',Z);
writemat (fidlog,'L-matrix, as used in ''search''',L);
writemat (fidlog,'D-matrix, as used in ''search''',D);
writemat (fidlog,'Transformed reduced ambiguities',z);

% ----------------------------------------------------------------------
% If the number of requested candidates is no more than "n+1", with "n"
% the dimension of Q_\hat(a), we can compute a suitable Chi^2 such that
% we have the requested number of candidates at minimum; use an eps
% to make sure the candidates are inside the ellipsoid and not exactly
% on the border.
%
% If the requested number of canditates is larger, we have to "guess"
% the volume of the search-space, it can not be guarenteed the requested
% number of canditates will actually be found.
% ----------------------------------------------------------------------

Chi2 = chistart (D,L,afloat,ncands);
writemat (fidlog,'Chi-squared, size of the search-ellipsoid',Chi2);

% ----------------------------------------------------------------------
% Find the requested number of candidates (search)
% ----------------------------------------------------------------------

[afixed,sqnorm,ierr] = lsearch (afloat,L,D,Chi2,ncands);

if ierr == 1;
   fprintf (1,'Warning: Not enough candidates were found!!\n');
   fprintf (1,'         Requested number.: %2d\n',ncands);
end;

writemat (fidlog,'Fixed transformed ''reduced'' ambiguities',afixed);
writemat (fidlog,'Corresponding squared norms',sqnorm);

% ----------------------------------------------------------------------
% Finally, perform the back-transformation and add the increments
% ----------------------------------------------------------------------

afixed = (afixed' * inv(Z))';
afixed = round(afixed + repmat(incr,1,ncands));
writemat (fidlog,'Final fixed ambiguities',afixed);

% ----------------------------------------------------------------------
% End of routine: lambda1
% ----------------------------------------------------------------------

author: goGPS-Project


text:
classdef DehomogenizationTestsSuite < handle

    methods

        function obj = DehomogenizationTestsSuite()
            warning('off', 'MATLAB:structOnObject')
            results = runtests("DehomogenizingTests","Tag","Singularities", 'Debug', true);
            table(results)
            warning('on', 'MATLAB:structOnObject')
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function variance = cardioid_variance ( x, a, b )

%*****************************************************************************80
%
%% CARDIOID_VARIANCE returns the variance of the Cardioid PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 <= B <= 0.5.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = a;

  return
end


text:
function table = i4mat_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_READ reads information from an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  [ m, n ] = i4mat_header_read ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_READ - Fatal error!\n' );
    fprintf ( 1, '  The spatial dimension is not positive.\n' );
  end

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_READ - Fatal error!\n' );
    fprintf ( 1, '  The number of points is not positive.\n' );
  end

  table(1:m,1:n) = i4mat_data_read ( input_filename, m, n );

  return
end


text:
function test_toolbox
%TEST_TOOLBOX Tests all functionalities of the dimension reduction toolbox
%
%   test_toolbox
%
% Tests all functionalities of the dimension reduction toolbox.
%
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, Delft University of Technology

    % Generate data
    disp('Testing data generation functions...');
    datasets = {'helix', 'twinpeaks', '3d_clusters', 'intersect', 'swiss'};
    for i=1:length(datasets)
        try
            X = generate_data(datasets{i}, 500);
        catch e
            disp(e);
            warning(['Generation of data set ' datasets{i} ' failed! Press any key to continue tests...']);
            pause
        end
    end

    % Test prewhitening
    disp('Testing prewhitening...');
    try
        X = prewhiten(X);
    catch e
        disp(e);
        warning('Prewhitening failed! Press any key to continue tests...');
        pause
    end
    unscaled_X = X;
    X = X - min(X(:));
    X = X / max(X(:));

    % Test all intrinsic dimensionality estimators
    disp('Testing intrinsic dimensionality estimators...');
    techniques = {'CorrDim', 'NearNbDim', 'GMST', 'PackingNumbers', 'EigValue', 'MLE'};
    for i=1:length(techniques)
        try
            intrinsic_dim(X, techniques{i});
        catch e
            disp(e);
            warning(['Intrinsic dimensionality estimation using ' techniques{i} ' failed! Press any key to continue tests...']);
            pause
        end
    end

    % Test all unsupervised dimension reduction techniques
    disp('Testing dimensionality reduction techniques...');
    techniques = {'PCA', 'MDS', 'ProbPCA', 'FactorAnalysis', 'GPLVM', 'Sammon', 'Isomap', ...
        'LandmarkIsomap', 'LLE', 'Laplacian', 'HessianLLE', 'LTSA', 'MVU', 'CCA', 'LandmarkMVU', ...
        'FastMVU', 'DiffusionMaps', 'KernelPCA', 'GDA', 'SNE', 'SymSNE', 'tSNE', 'LPP', 'NPE', ...
        'LLTSA', 'SPE', 'Autoencoder', 'LLC', 'ManifoldChart', 'CFA'};
    for i=1:length(techniques)

        % Test the dimension reduction technique
        try
            if any(strcmpi(techniques{i}, {'GPLVM', 'CFA'}))
                [mappedX, mapping] = compute_mapping(unscaled_X, techniques{i}, 2);
            else
                [mappedX, mapping] = compute_mapping(X, techniques{i}, 2);
            end
            if any(strcmpi(techniques{i}, {'Isomap', 'LandmarkIsomap', 'LLE', 'Laplacian', 'MVU', 'CCA', 'FastMVU', 'LPP', 'NPE', 'LLTSA'}))
                [mappedX, mapping] = compute_mapping(X, techniques{i}, 2, 'adaptive');
            end
        catch e
            disp(e);
            warning(['Technique ' techniques{i} ' failed! Press any key to continue tests...']);
            pause
        end

        % Test the out-of-sample extension code
        if any(strcmpi(techniques{i}, {'PCA', 'LPP', 'NPE', 'LLTSA', 'SPCA', 'PPCA', 'FA'}))
            try
                out_of_sample(X, mapping);
            catch e
                disp(e);
                warning(['Out-of-sample extension for technique ' techniques{i} ' failed! Press any key to continue tests...']);
                pause
            end
        end

        % Test reconstruction code
        if any(strcmpi(techniques{i}, {'PCA', 'LPP', 'NPE', 'LLTSA', 'SPCA', 'PPCA', 'FA', 'Autoencoder'}))
            try
                reconstruct_data(mappedX, mapping);
            catch e
                disp(e);
                warning(['Reconstruction for technique ' techniques{i} ' failed! Press any key to continue tests...']);
                pause
            end
        end
    end

    % Test approximate out-of-sample function
    try
        out_of_sample_est(X, X, mappedX);
    catch e
        disp(e);
        warning(['Approximate out-of-sample extension failed! Press any key to continue tests...']);
        pause
    end

    % Test all supervised dimension reduction techniques
    labels = double(X(:,1) > .5) + 1;
    X = [labels X];
    techniques = {'LDA', 'NCA', 'MCML', 'LMNN'};
    for i=1:length(techniques)

        % Test the actual technique
        try
            [mappedX, mapping] = compute_mapping(X, techniques{i}, 2);
        catch e
            disp(e);
            warning(['Technique ' techniques{i} ' failed! Press any key to continue tests...']);
            pause
        end

        % Test out-of-sample extension
        try
            out_of_sample(X(:,2:end), mapping);
        catch e
            disp(e);
            warning(['Out-of-sample extension for technique ' techniques{i} ' failed! Press any key to continue tests...']);
            pause
        end

        % Test reconstruction code
        try
            reconstruct_data(mappedX, mapping);
        catch e
            disp(e);
            warning(['Reconstruction for technique ' techniques{i} ' failed! Press any key to continue tests...']);
            pause
        end

    end
    disp('All tests completed!');

repo: cluster


text:
function [out,Wc,Wo] = hankelsv(SYS)
%HANKELSV Compute Hankel singular values and grammians.
%
% [OUT,Wc,Wo] = HANKELSV(SYS) computes controllability and observability
%    grammians Wc, Wo, and the Hankel singular values OUT of an LTI
%    model SYS (created with either TF, ZPK, SS, or FRD). The model
%    SYS can be either continuous-time or discrete-time. However, only
%    in continuous-time case that SYS is allowed to be unstable.
%    The computed Hankel singular values are sorted in ascending order.
%
%    For unstable continuous-time system, state-space stable/anti-stable
%    decomposition is used instead, and OUT=[OUT_stable;OUT_anti-stable].
%    In addition, Wc={Wc_stable,Wc_anti-stable} and Wo is either.

%    The former version of HKSV employs an obsolete fashion in using the
%    Matlab function GRAM, i.e., instead of using GRAM(SYS,'c') and
%    GRAM(SYS,'o'), it uses GRAM(A,B) and GRAM(A',C'), respectively.
%    This restricts the computation of gramians to continuous-time case,
%    only, since GRAM(A,B) and GRAM(A',C') solve LYAP, not DLYAP.
%
%    This improved version also correct some other bugs in the former
%    version, for example, HKSV does not return gramians when SYS is
%    unstable, but not anti-stable.

%    Developer: Wathanyoo Khaisongkram
%    Date Developed: Oct 20, 2004
%    Email: sunboom15@yahoo.com
%    Improved version of HKSV by R. Y. Chiang & M. G. Safonov March, 1986
% -----------------------------------------------------------------------

SYS=ss(SYS);                            % convert to state-space model

% Discrete-time LTI model
if get(SYS,'Ts')~=0                     % discrete-time LTI model
   if max(abs(pole(SYS))) > 1           % unstable discrete system
      error('For discrete LTI model, system must be stable')
   end
   Wc = gram(SYS,'c');                  % controllability matrix
   Wo = gram(SYS,'o');                  % observability matrix
   out = sqrt(eig(Wc*Wo));              % Hankel singular values
   [no_use,index] = sort(out);          % index for sorting
   out = out(index);                    % sorted singular values
   return                               % end function
end

% Continuous-time LTI model
[A,B,C,D]=ssdata(SYS);                  % extract the system matrices
   mode = eig(A);                       % the eigen values of 'a'
   nrow = size(A,1);                    % the number of rows in 'a'
   nsta = length(find(real(mode) < 0)); % the number of unstable modes
if isequal(nsta,nrow),                  % completely stable
   Wc = gram(SYS,'c');                  % controllability matrix
   Wo = gram(SYS,'o');                  % observability matrix
   out = sqrt(eig(Wc*Wo));              % Hankel singular values
   [no_use,index] = sort(out);          % index for sorting
   out = out(index);                    % sorted singular values
elseif isequal(nsta,0),                 % completely unstable
   SYS=ss(-A,-B,C,D,SYS);               % change the dynamic and the input matrices to '-a' and '-b',
                                        % respectively with all LTI properties inherited from
                                        % original SYS, e.g., discrete or continuous domain.
   Wc = gram(SYS,'c');                  % controllability matrix
   Wo = gram(SYS,'o');                  % observability matrix
   out = sqrt(eig(Wc*Wo));              % Hankel singular values
   [no_use,index] = sort(out);          % index for sorting
   out = out(index);                    % sorted singular values
else, % 0 < nsta < nrow
   [SYSs,SYSu] = stabproj(SYS);         % decompose stable/anti-stable parts
   % Stable part
   Wcs = gram(SYSs,'c');                % controllability matrix
   Wos = gram(SYSs,'o');                % observability matrix
   outl = sqrt(eig(Wcs*Wos));           % Hankel singular values1
   [no_use,index] = sort(outl);         % index for sorting
   outl = outl(index);                  % sorted singular values
   % Anti-stable part
   [Au,Bu,Cu,Du]=ssdata(SYSu);          % obtain the system matrices of SYSu
   SYSu=ss(-Au,-Bu,Cu,Du,SYSu);         % change the dynamic and the input matrices to '-ar' and '-br'.
   Wcu = gram(SYSu,'c');                % controllability matrix
   Wou = gram(SYSu,'o');                % observability matrix
   outr = sqrt(eig(Wcu*Wou));           % Hankel singular values
   [no_use,index] = sort(outr);         % index for sorting
   outr = outr(index);                  % sorted singular values
   out = [outl;outr];                   % gather two cases
   Wc={Wcs,Wcu};                        % controllability gramian
   Wo={Wos,Wou};                        % observability gramian
end

% end of HKSV_MOD


text:
function incorrectGeneRules = testGeneRules(model)
% Finds gene rules that have an incorrect nomenclature.
%
% INPUT
% model             COBRA model structure
%
% OUTPUT
% incorrectGeneRules       Cell array listing entries in model.rules that
%                          have incorrect nomenclature.
%
% Almut Heinken, Oct 2019

cnt=1;
incorrectGeneRules={};

for i=1:length(model.rules)
   tf = verifyRuleSyntax(model.rules{i});
   if ~tf
       incorrectGeneRules{cnt}=i;
       cnt=cnt+1;
   end
end

end

author: opencobra

repo: cobratoolbox


text:
function mmldemo
% MMLDEMO is a guided tour about using MML.
disp(' ')
disp('  MMLDEMO')
disp(' ')
disp('  MMLDEMO demonstrates the capabilities of the MML function and how to use it.')
disp('  This demo guides through a number of examples showing the syntactical')
disp('  variations MML can be used with.')
disp(' ')
disp('  To run properly this demo requires:')
disp('  1.) Symbolic Math Toolbox 3.0.1 for MATLAB with Maple 8 Kernel or later.')
disp('  2.) To display the expression in MS Internet Explorer install MATHPLAYER.')
disp('      Download MATHPLAYER (freeware) at')
disp('      http://www.dessci.com/webmath/mathplayer')
disp(' ')
disp('[Press any key to continue with Example 1]')
pause
disp('  ')
disp('  Example 1: MML(TERM)')
disp('  ')
disp('  The following line in the MATLAB Command Window displays a mathematical')
disp('  expression in the system''s standard webbrowser. In this example the area')
disp('  function of a circle is entered for TERM. When the line below is executed')
disp('  the webbrowser opens and displays the term.')
disp('  ')
disp('  >>mml(pi/4*d^2)')
disp('  ')
disp('[Press any key to demonstrate]')
pause
disp('  ')
disp('  Now the webbrowser opens and displays the term.')
mml('pi/4*d^2')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  In the example above the term was entered as a symbolic MATLAB expression.')
disp('  Using symbolic expressions is subject to certain restrictions of MATLAB syntax.')
disp('  For example it is not valid to enter equations, like')
disp('  >>mml(a=pi/4*d^2)')
disp('  To avoid some restrictions it is possible to enter the term as a character ')
disp('  string using apostrophes, like in the following example. The character string is')
disp('  passed to the Symbolic Toolbox''s Maple Kernel for interpretation. The term has ')
disp('  to be entered in valid Maple or MATLAB syntax.')
disp('  ')
disp('  >>mml(''a=pi/4*d^2'')')
disp('  ')
disp('[Press any key to demonstrate]')
pause
disp('  ')
disp('  Now the webbrowser opens and displays the term.')
mml('a=pi/4*d^2')
disp('  ')
disp('  ')
disp('[Press any key to continue with Example 2]')
pause
disp('  ')
disp('  Example 2: MML')
disp('  ')
disp('  A useful feature of MML is the possibility of displaying the current content ')
disp('  of the global variable ANS by simply typing MML.')
disp('  To demonstrate this we assume that ANS contains i.e. a symbolic matrix.')
disp('  ')
disp('  >> [x1 x2 x3; y1 y2 y3]')
disp('  ')
disp('  ans =')
disp('  ')
disp('  [ x1, x2, x3]')
disp('  [ y1, y2, y3]')
disp('  ')
disp('  >>mml')
disp('  ')
disp('[Press any key to demonstrate]')
pause
disp('  ')
disp('  Now the webbrowser opens and displays the content of ANS, here: the matrix.')
mml('[x1 x2 x3; y1 y2 y3]')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  Example 3: Options')
disp('  ')
disp('  This example demonstrates the use of the display options ''type'' and ''nwb''.')
disp('  Both options can be used in each syntax variation of the MML function.')
disp('  If MathML code typed in the Command Window is needed in addition to displaying')
disp('  the term in the webbrowser, the option ''type'' can be used.')
disp('  ')
disp('  >>mml(sqrt(a),''type'')')
disp('  ')
disp('[Press any key to demonstrate]')
pause
disp('  ')
disp('  Now the webbrowser is opened displaying the term and the corresponding ')
disp('  MathML code is typed in the Command Window.')
disp('  ')
disp('  ')
mml('sqrt(a)','type')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  The following example shows the use of the ''type'' option with the syntax ')
disp('  of EXAMPLE 2: MML.')
disp('  ')
disp('  >> sqrt(a)')
disp('  ')
disp('  ans =')
disp('  ')
disp('  sqrt(a)')
disp('  ')
disp('  >>mml(''type'')')
disp('  ')
disp('[Press any key to demonstrate]')
pause
disp('  ')
disp('  Now the webbrowser is opened displaying the content of ANS ')
disp('  and the corresponding MathML code is typed in the Command Window.')
disp('  ')
disp('  ')
mml('sqrt(a)','type')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  If only MathML code typed in the Command Window is needed instead of displaying')
disp('  the term in the webbrowser, the option ''nwb''("No Web Browser") can be used.')
disp('  ')
disp('  >>mml(sqrt(a),''nwb'')')
disp('  ')
disp('[Press any key to demonstrate]')
disp('  ')
pause
disp('  The MathML code of the term is now typed in the Command Window.')
disp('  The webbrowser is not opened.')
disp('  ')
disp('  ')
mml('sqrt(a)','nwb')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  The following example shows the use of the ''nwb'' option with the syntax ')
disp('  of EXAMPLE 2: MML.')
disp('  ')
disp('  >> sqrt(a)')
disp('  ')
disp('  ans =')
disp('  ')
disp('  sqrt(a)')
disp('  ')
disp('  >>mml(''nwb'')')
disp('  ')
disp('[Press any key to demonstrate]')
disp('  ')
pause
disp('  Now the MathML code of the content of ANS is typed in the Command Window.')
disp('  The webbrowser is not opened.')
disp('  ')
disp('  ')
mml('sqrt(a)','nwb')
disp('  ')
disp('  ')
disp('[Press any key to continue]')
pause
disp('  ')
disp('  Example 4: "Torture Test"')
disp('  ')
disp('  The following example demonstrates a term with higher complexity. ')
disp('  The displayed term makes no mathematical sense. It has been choosen to contain')
disp('  a number of different mathematical aspects.')
disp('  ')
disp('  >>mml(''g(x+y)<=asinh(exp(15^(1/3)*pi^factorial(x)))/(a*x^2-b*x+c)+2/3-2.5*y'')')
disp('  ')
disp('[Press any key to demonstrate]')
disp('  ')
pause
mml('g(x+y)<=asinh(exp(15^(1/3)*pi^factorial(x)))/(a*x^2-b*x+c)+2/3-2.5*y')
disp('  ')
disp('  End of MMLDEMO.')
disp('  ')
disp('  See also: help MML ')


text:
classdef ACO < ALGORITHM
% <single> <permutation> <large/none>
% Ant colony optimization

%------------------------------- Reference --------------------------------
% M. Dorigo and G. D. Caro, Ant colony optimization: a new meta-heuristic,
% Proceedings of the IEEE Congress on Evolutionary Computation, 1999.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            Tau = ones(Problem.D);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Move the ants
                PopDec = zeros(Problem.N,Problem.D);
                PopDec(:,1) = randi(Problem.D,Problem.N,1);
                for i = 1 : Problem.N
                    for j = 1 : Problem.D-1
                        Remain = setdiff(1:Problem.D,PopDec(i,1:j));
                        next = RouletteWheelSelection(1,Problem.C(PopDec(i,j),Remain)./Tau(PopDec(i,j),Remain));
                        PopDec(i,j+1) = Remain(next);
                    end
                end
                Population = Problem.Evaluation(PopDec);
                % Update the pheromone matrix
                dTau = zeros(Problem.D) + 1e-6;
                for i = 1 : Problem.N
                    for j = 1 : Problem.D-1
                        dTau(PopDec(i,j),PopDec(i,j+1)) = dTau(PopDec(i,j),PopDec(i,j+1)) + 1/Population(i).obj;
                    end
                    dTau(PopDec(i,end),PopDec(i,1)) = dTau(PopDec(i,end),PopDec(i,1)) + 1/Population(i).obj;
                end
                Tau = 0.5*Tau + dTau;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [proj,geo,angles]=loadNikonProjections(filepath,geo,angles,varargin)
%[proj,angles]=loadNikonProjections(filepath,geo,angles,varargin)
%   loads Nikon uCT machine projections
%
%   loadNikonData(filepath,geo,angles) Loads a dataset given its FILEPATH,
%      a geometry GEO and angles ANGLES (loaded from readXtekctGeometry())
%
%   loadNikonData(filepath,geo,angles, OPT,VAL,...) uses options and values.
%      These are options in case you don't want to load the entire
%      dataset, but only particular sets of projections.
%      The possible options in OPT are:
%           'sampling': type of sampling. default 'equidistant' Can be:
%                  'equidistant': equidistantly sample the entire set
%                       of angles. 'num_angles' should be set for partial
%                       loading of the data.
%                  'step': sample the entire set of projections every
%                         'sampling_step' angles.
%                  'continous': Load the first 'num_angles' amount of
%                             angles only.
%
%           'num_angles': Number of total angles to load. Default all of
%                    them. Useful for 'equidistant' and 'continous' loading
%
%           'sampling_step': step to load when loading projections. Default
%           1. Useful for 'step' loading.
%

% developed by A. Biguri and W. Sun 06.07.2020

%% Parse inputs.

[angles_to_load,index]=parse_inputs(geo,angles,varargin);

% make sure its path
if filepath(end)~='\' && filepath(end)~='/'
   filepath=[filepath '/'];
end
%% get filename
% assuming TIF and 4 digits.
firstfile = dir([filepath,'/*.tif']); %
firstfile = firstfile(1).name;
filename = firstfile(1:end-8);

fprintf("Dataset in: %s \n", filepath);

%% load images
%proj=[];
l = length(angles_to_load);
proj = zeros(geo.nDetector(2),geo.nDetector(1),l,'single');
for ii=1:length(angles_to_load)
    if(~mod(ii,50))
      fprintf("Loading: %d / %d \n", ii, length(angles_to_load));
    end
     proj(:,:,ii)=single(imread([filepath,filename,num2str(index(ii),'%04d'),'.tif']));
end
%% Beer lambert
if any(proj(:))>single(geo.whitelevel)
    warning('Changing the Whitelevel value as projection data has higher value than specified in the file');
    geo.whitelevel=max(proj(:)+1);
end
proj=-log(proj/single(geo.whitelevel));
geo=rmfield(geo,'whitelevel');
%%
angles=angles_to_load;
end

function [angles,indices]=parse_inputs(geo,angles,argin)
opts=     {'sampling','num_angles','sampling_step'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error(['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option is not default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error(['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        case 'sampling'
            if default
                sampling='equidistant';
            else
                if strcmpi(val,'equidistant') || strcmpi(val,'step')|| strcmpi(val,'continuous')
                    sampling=val;
                else
                   error('sampling type not understood, must be  equidistant, step or continuous');
                end
            end
        case 'num_angles'
            if default
                num_angles=length(angles);
            else
                num_angles=val;
            end
        case 'sampling_step'
            if default
               sampling_step=1;
            else
               sampling_step=val;
            end
        otherwise
            error(['Invalid input name:', num2str(opt),'\n No such option']);
    end
end
% now lets build the angles

if strcmpi(sampling,'equidistant')
   sampling_step=round(length(angles)/num_angles);
   indices=1:sampling_step:length(angles);
   angles=angles(indices);
   return;
end
if strcmpi(sampling,'continuous')
   indices=1:num_angles;
   angles=angles(indices);
   return;
end
if strcmpi(sampling,'step')
   indices=1:sampling_step:length(angles);
   angles=angles(indices);
   return;
end
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MATLAB Code for                                                  %
%                                                                   %
%  Multi-Objective Particle Swarm Optimization (MOPSO)              %
%  Version 1.0 - Feb. 2011                                          %
%                                                                   %
%  According to:                                                    %
%  Carlos A. Coello Coello et al.,                                  %
%  "Handling Multiple Objectives with Particle Swarm Optimization," %
%  IEEE Transactions on Evolutionary Computation, Vol. 8, No. 3,    %
%  pp. 256-279, June 2004.                                          %
%                                                                   %
%  Developed Using MATLAB R2009b (Version 7.9)                      %
%                                                                   %
%  Programmed By: S. Mostapha Kalami Heris                          %
%                                                                   %
%         e-Mail: sm.kalami@gmail.com                               %
%                 kalami@ee.kntu.ac.ir                              %
%                                                                   %
%       Homepage: http://www.kalami.ir                              %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function G=CreateHypercubes(costs,ngrid,alpha)

    nobj=size(costs,1);

    empty_grid.Lower=[];
    empty_grid.Upper=[];
    G=repmat(empty_grid,nobj,1);

    for j=1:nobj

        min_cj=min(costs(j,:));
        max_cj=max(costs(j,:));

        dcj=alpha*(max_cj-min_cj);

        min_cj=min_cj-dcj;
        max_cj=max_cj+dcj;

        gx=linspace(min_cj,max_cj,ngrid-1);

        G(j).Lower=[-inf gx];
        G(j).Upper=[gx inf];

    end

end

author: sfvsfv

repo: Mathematical-modeling


text:
function [s,t]=HFMChirp(T,fStart,fEnd,tSamp,phi)
%%LFMCHIRP Generate a complex hyperbolic frequency modulated (HFM) signal.
%          It can be an upchirp or a downchirp. HFM chirps have a role in
%          radar and sonar signal processing. They are more Doppler
%          insensitive than linear FM chirps.
%
%INPUTS: T The duration of the chirp (seconds).
%   fStart The frequency at the start of the chirp (Hertz). This must be
%          >0.
%     fEnd The frequency at the end of the chirp (Hertz). This must be >0.
%    tSamp If this is scalar, it is the sampling rate of the chirp and the
%          number of samples depends on T. Otherwise, this is a numSampX1
%          array of times at which samples are desired. Samples before
%          time 0 and after time T are set to zero. If this parameter is
%          omitted or an empty matrix is passed, the Nyquist sampling
%          rate is used.
%      phi The phase offset of the chirp. If this parameter is omitted or
%          an empty matrix is passed, then the default of pi/2 is used,
%          which makes the real part of the signal start at zero.
%
%OUTPUT: s The numSampX1 complex chirp signal.
%        t The sample times used. If tSamp was not a scalar, then this is
%          just tSamp.
%
%HFM chirps are described in [1] along with a range bias model for target
%tracking that is relevant when Doppler filtering is not performed.
%
%EXAMPLE 1:
%Here is an example of an upchirp:
% T=1e-3;
% fStart=1e3;
% fEnd=50e3;
% tSamp=1/(100*2*fEnd);%100 times the Nyquist frequency.
% [s,t]=HFMChirp(T,fStart,fEnd,tSamp);
% figure(1)
% clf
% plot(t,real(s),'linewidth',2)
% h1=xlabel('t');
% h2=ylabel('s');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')

%EXAMPLE 2:
%Here is an example of a downchirp:
% T=1e-3;
% fStart=50e3;
% fEnd=1e3;
% tSamp=1/(100*2*fStart);%100 times the Nyquist frequency.
% [s,t]=HFMChirp(T,fStart,fEnd,tSamp);
% figure(2)
% clf
% plot(t,real(s),'linewidth',2)
% h1=xlabel('t');
% h2=ylabel('s');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] X. Song, P. Willett, and S. Zhou, "Range bias modeling for hyperbolic-
%    frequency-modulated waveforms in target tracking," IEEE Journal of
%    Oceanic Engineering, vol. 37, no. 4, pp. 670-679, Oct. 2012.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(tSamp))
   %Default to the Nyquist sampling rate.
   tSamp=1/(2*max(fStart,fEnd));
end

if(nargin<5||isempty(phi))
   phi=pi/2;%Phase offset so the signal starts at zero, by default.
end

%Equation 1
b=(fStart-fEnd)/(fStart*fEnd*T);

if(isscalar(tSamp))
    %If a scalar sampling rate is given.
    NPulse=fix(T/tSamp);
    t=(0:(NPulse-1))*tSamp;
else
    %If the samples themselves are given
    t=tSamp;
end

s=exp(1j*(phi+(2*pi/b)*log(1+b*fStart*t)));

sel=(t>T)&(t<0);
s(sel)=0;

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.