text:
function [p] = gui_mesh_open(p,command,parent)

% gui_mesh_open - GUI to mesh_open
%
% Usage: [p] = gui_mesh_open(p,[command],[parent])
%
% p        - struct, generated by 'eeg_toolbox_defaults'
% command  - either 'init' or 'load'
% parent   - a handle to the gui that calls this gui, used
%            to update the UserData of the parent. The p
%            struct is returned to the parent in this way.
%
% See command line version, mesh_open
%

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

% Licence:  GNU GPL, no express or implied warranties
% History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('p','var'),
 [p] = eeg_toolbox_defaults;
elseif isempty(p),
 [p] = eeg_toolbox_defaults;
end

if ~exist('command','var'),
  command = 'init';
elseif isempty(command),
  command = 'init';
end

command = lower(command);

switch command,

  case 'init',

    if exist('parent','var'),
      MESHOpen = INIT(p,parent);
    else
      MESHOpen = INIT(p);
    end

  case 'plot',

    MESHOpen = get(gcbf,'Userdata');

    MESHOpen.p = mesh_open(MESHOpen.p);

    set(MESHOpen.gui,'Userdata',MESHOpen);

    if isequal(get(MESHOpen.handles.Bhold,'Value'),1),

     [p] = gui_updateparent(MESHOpen,0);
     [p] = gui_mesh_plot(p,'init',MESHOpen.gui);

    else

     [p] = gui_updateparent(MESHOpen);
      if isfield(MESHOpen,'parent'),
        if ~isempty(MESHOpen.parent),
          gui_mesh_plot(MESHOpen.p,'init',MESHOpen.parent.gui);
        end
      else
        gui_mesh_plot(MESHOpen.p,'init','');
      end
      close gcbf;
    end

  case 'save',

    MESHOpen = get(gcbf,'Userdata');
    gui_mesh_write(MESHOpen.p);

    if isequal(get(MESHOpen.handles.Bhold,'Value'),1),
     [p] = gui_updateparent(MESHOpen,0);
    else
     [p] = gui_updateparent(MESHOpen);
      close gcbf;
    end

  case 'return',

    MESHOpen = get(gcbf,'Userdata');

    set(MESHOpen.gui,'Pointer','watch');

    MESHOpen.p = mesh_open(MESHOpen.p);

    set(MESHOpen.gui,'Userdata',MESHOpen);

    set(MESHOpen.gui,'Pointer','arrow');

    if isequal(get(MESHOpen.handles.Bhold,'Value'),1),
     [p] = gui_updateparent(MESHOpen,0);
    else
     [p] = gui_updateparent(MESHOpen);
      close gcbf;
    end

  otherwise,

    close gcbf;

end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [MESHOpen] = INIT(p,parent),

% GUI General Parameters

GUIwidth  = 500;
GUIheight = 120;

version = '$Revision: 1.1 $';
name = sprintf('MESH File Open [v %s]\n',version(11:15));

GUI = figure('Name',name,'Tag','MESH_OPEN',...
  'NumberTitle','off',...
  'MenuBar','none','Position',[1 1 GUIwidth GUIheight]);
movegui(GUI,'center');

MESHOpen.gui = GUI;

Font.FontName   = 'Helvetica';
Font.FontUnits  = 'Pixels';
Font.FontSize   = 12;
Font.FontWeight = 'normal';
Font.FontAngle  = 'normal';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mesh Data Selection and Parameters

G.Title_data = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .75 .17 .2],...
  'TooltipString','Mesh/Tesselation Formats (see, ''help mesh_open'')',...
  'String','Mesh Type:','HorizontalAlignment','left');

switch lower(p.mesh.type),
  case 'brainstorm',          meshType = 1;
  case 'emse',                meshType = 2;
  case 'ascii',               meshType = 3;
  case 'freesurfer_ascii',    meshType = 4;
  case 'freesurfer_surf',     meshType = 5;
  case 'freesurfer_curv',     meshType = 6;
  case 'freesurfer_overlay',  meshType = 7;
  otherwise,                  meshType = 1;
end

meshtypes = {'BrainStorm' 'EMSE' 'ASCII' ...
    'FreeSurfer_ascii' 'FreeSurfer_surf' 'FreeSurfer_curv' 'FreeSurfer_overlay' };
G.PmeshType = uicontrol('Tag','PmeshType','Parent',GUI,'Style','popupmenu',...
  'Units','Normalized',Font,  ...
  'Position',[.20 .75 .28 .2],...
  'String',meshtypes,'Value',meshType,...
  'TooltipString','Mesh/Triangulation Formats (see, ''help mesh_open'', FS = FreeSurfer)',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'MESHOpen.p.mesh.type = popupstr(MESHOpen.handles.PmeshType);',...
  'set(gcbf,''Userdata'',MESHOpen); clear MESHOpen;'));

G.BmeshReplace = uicontrol('Tag','BmeshReplace','Parent',GUI,'Style','checkbox',...
  'Units','Normalized',Font,...
  'Position',[.40 .75 .18 .2],...
  'String','Replace','Value',1,...
  'TooltipString','Reload current mesh, if any?',...
  'Visible','off',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'MESHOpen.p.mesh.reload = get(MESHOpen.handles.BmeshReplace,''Value'');',...
  'set(gcbf,''Userdata'',MESHOpen); clear MESHOpen;'));

G.Title_path = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .50 .17 .2],...
  'String','Path','HorizontalAlignment','left');
G.EmeshPath = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
  'Position',[.20 .50 .58 .2], 'String',p.mesh.path,...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'MESHOpen.p.mesh.path = get(MESHOpen.handles.EmeshPath,''String'');',...
  'set(gcbf,''Userdata'',MESHOpen); clear MESHOpen;'));

G.Title_file = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
  'Position',[.01 .25 .17 .2],...
  'String','File','HorizontalAlignment','left');
G.EmeshFile = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
  'Position',[.20 .25 .58 .2], 'String',p.mesh.file,...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'MESHOpen.p.mesh.file = get(MESHOpen.handles.EmeshFile,''String'');',...
  'set(gcbf,''Userdata'',MESHOpen); clear MESHOpen;'));

Font.FontWeight = 'bold';

% BROWSE: Look for the data
browsecommand = strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'cd(MESHOpen.p.mesh.path);',...
  '[file, path] = uigetfile(',...
  '{''*.dat;*.txt'', ''ASCII Tesselations (*.dat,*.txt)'';', ...
  ' ''*.asc;*.tri'', ''FreeSurfer Tesselations (*.asc,*.tri)'';', ...
  ' ''*.orig;*.smoothwm;*.white;*.pial'', ''FreeSurfer Surfaces'';', ...
  ' ''*.wfr'', ''EMSE Tesselations (*.wfr)'';', ...
  ' ''*.mat'', ''BrainStorm Tesselations (*.mat)'';', ...
  ' ''*.*'',   ''All Files (*.*)''},', ...
  '''Select Mesh/Tesselation File'');',...
  'if ~isequal(path,0), MESHOpen.p.mesh.path = path; end;',...
  'if ~isequal(file,0), MESHOpen.p.mesh.file = file; end;',...
  'set(MESHOpen.handles.EmeshPath,''String'',MESHOpen.p.mesh.path);',...
  'set(MESHOpen.handles.EmeshFile,''String'',MESHOpen.p.mesh.file);',...
  'if     findstr(file,''.mat''), value = 1; ',...
  'elseif findstr(file,''.wfr''), value = 2; ',...
  'elseif findstr(file,''.dat''), value = 3; ',...
  'elseif findstr(file,''.txt''), value = 3; ',...
  'elseif findstr(file,''.asc''), value = 3; ',...
  'elseif findstr(file,''.tri''), value = 4; ',...
  'elseif findstr(file,''.orig''),     value = 5; ',...
  'elseif findstr(file,''.smoothwm''), value = 5; ',...
  'elseif findstr(file,''.white''),    value = 5; ',...
  'elseif findstr(file,''.pial''),     value = 5; ',...
  'else, value = 3; ',...
  'end;',...
  'set(MESHOpen.handles.PmeshType,''Value'',value); ',...
  'MESHOpen.p.mesh.type = popupstr(MESHOpen.handles.PmeshType);',...
  'set(gcbf,''Userdata'',MESHOpen); clear value file path MESHOpen;');
G.BmeshFile = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized',Font, ...
  'Position',[.01 .01 .17 .2], 'String','BROWSE',...
  'BackgroundColor',[0.8 0.8 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback', browsecommand );

% PLOT: Load & plot the data!
G.Bplot = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.20 .01 .18 .2],...
  'String','PLOT','BusyAction','queue',...
  'TooltipString','Plot the mesh and return the file parameters.',...
  'BackgroundColor',[0.0 0.5 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'p = gui_mesh_open(MESHOpen.p,''plot'');',...
  'clear MESHOpen;'));

% SAVE AS: Load & convert the data!
G.Bsave = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.40 .01 .18 .2],...
  'String','SAVE AS','BusyAction','queue',...
  'TooltipString','Load and convert the mesh.',...
  'BackgroundColor',[0.0 0.0 0.5],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'p = gui_mesh_open(MESHOpen.p,''save'');',...
  'clear MESHOpen;'));

% Return file parameters
G.Breturn = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.60 .01 .18 .2],...
  'String','RETURN','BusyAction','queue',...
  'TooltipString','Close and return the file parameters.',...
  'BackgroundColor',[0.75 0.0 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'p = gui_mesh_open(MESHOpen.p,''return'');',...
  'clear MESHOpen;'));

% Cancel, do not return parameters
G.Bcancel = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.80 .01 .18 .2],...
  'String','CANCEL','BusyAction','queue',...
  'TooltipString','Close, do not return the p struct.',...
  'BackgroundColor',[0.75 0.0 0.0],...
  'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
  'Callback',strcat('MESHOpen = get(gcbf,''Userdata'');',...
  'p = gui_mesh_open(MESHOpen.p,''cancel'');',...
  'clear MESHOpen;'));

% Help
G.Bhelp = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
  'Position',[.80 .25 .18 .2],'String','Help','BusyAction','queue',...
  'BackgroundColor',[1 1 0],...
  'ForegroundColor',[0 0 0], 'HorizontalAlignment', 'center',...
  'Callback','doc mesh_open;');

% Hold GUI Open checkbox
G.Bhold = uicontrol('Parent',GUI,'Style','checkbox','Units','Normalized', Font, ...
  'Position',[.80 .50 .18 .2],'String','Hold GUI','BusyAction','queue',...
  'TooltipString','MESH File Load GUI remains open after ''Plot'' and ''Return'' commands.',...
  'Value',p.hold,'HorizontalAlignment', 'center');

% Store userdata
if exist('parent','var'), MESHOpen.parent.gui = parent; end
MESHOpen.handles = G;
MESHOpen.p = p;
set(GUI,'HandleVisibility','callback');
set(GUI,'Userdata',MESHOpen);

return

author: PatternRecognition

repo: OpenBMI


text:
function [traj] = updateTraj(traj,Rot,v,x,omega_b,a_b,i)
%record trajectory
angles = rotm2eul(Rot);
traj.Rot(:,:,i) = Rot;
traj.phi(i) = angles(3);
traj.theta(i) = angles(2);
traj.psi(i) = angles(1);
traj.v(:,i) = v;
traj.x(:,i) = x;
traj.omega_b(:,i) = omega_b;
traj.a_b(:,i) = a_b;
end

author: mbrossar


text:
function [p,npix] = histroi(f,c,r)
%HISTROI Computes the histogram of an ROI in an image.
%   [P,NPIX] = HISTROI(F,C,R) computes the histogram, P, of a polygonal
%   region of interest (ROI) in image F. The polygonal region is defined
%   by the column and row coordinates of its vertices, which are
%   specified (sequentially) in vectors C and R, respectively. All
%   pixels of F must be >= 0. Parameter NPIX is the number of pixels in
%   the polygonal region. For consistency with Toolbox function imhist,
%   histogram p is unnormalized. To normalize it, let p = p/npix.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Generate the binary mask image.
B = roipoly(f,c,r);

% Compute the histogram of the pixels in the ROI.
p = imhist(f(B));

% Obtain the number of pixels in the ROI if requested in the output. All
% ROI pixels have value 1.
if nargout > 1
   npix = sum(B(:));
end

author: dipum

repo: dipum-toolbox

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


text:
function midrad(p)
%MIDRAD       Display of polynom in "midrad" mode
%

% written  07/22/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  olddisp = intvalinit('display',0);
  intvalinit('displaymidrad',0);
  display(p,inputname(1));
  intvalinit(olddisp,0);

repo: OpenMAS


text:
function [yinterp,ypinterp] = ntrp1210(tinterp,t,y,tnew,ynew,klast,phi,psi,idxNonNegative)
%NTRP1210  Interpolation helper function for RKN1210.
%
%   YINTERP = NTRP1210(TINTERP,T,Y,TNEW,YNEW,KLAST,PHI,PSI,IDX) uses data
%   computed in RKN1210 to approximate the solution at time TINTERP. TINTERP
%   may be a scalar or a row vector.
%
%   [YINTERP,YPINTERP] = NTRP1210(TINTERP,T,Y,TNEW,YNEW,KLAST,PHI,PSI,IDX)
%   returns also the derivative of the polynomial approximating the solution.
%
%   IDX has indices of solution components that must be non-negative. Negative
%   YINTERP(IDX) are replaced with zeros and the derivative YPINTERP(IDX) is
%   set to zero.
%
%   See also RKN1210, DEVAL.

% References;
% [1] "Interpolating Runge-Kutta-Nyström Methods of High Order", C.
%     Tsitouras, G. Papageorgiou, Intern. J. Computer Math., vol 47,
%     pp. 209-217 (1993).

% Please report bugs and inquiries to:
%
% Name       : Rody P.S. Oldenhuis
% E-mail     : oldenhuis@gmail.com
% Licence    : 2-clause BSD (See Licence.txt)

% If you find this work useful, please consider a donation:
% https://www.paypal.me/RodyO/3.5

% TODO: make magic happen here

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function test_all(test_dir)
% @=============================================================================
% 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.
% =============================================================================@

%% ===== START =====
% Check inputs
if (nargin < 1) || isempty(test_dir)
    error('You must specify an empty directory in input.');
end
% Create test folder
if ~isdir(test_dir)
    if ~mkdir(test_dir)
        error(['Could not create folder: ' test_dir]);
    end
end
% Start profiling
% profile on
% Reset brainstorm
bst_reset();
% Prepare test folders
BrainstormHomeDir = fileparts(fileparts(fileparts(which(mfilename))));
BrainstormDbDir = fullfile(test_dir, 'brainstorm_db');
if isdir(BrainstormDbDir)
    error(['You should delete the test database folder first: ' BrainstormDbDir]);
end
mkdir(BrainstormDbDir);
% Start Brainstorm with GUI
bst_startup(BrainstormHomeDir, 0, BrainstormDbDir);

% Initialize random generator
if exist('rng', 'file')
    rng('default');
else
    sRnd = RandStream.getDefaultStream;
    sRnd.reset();
end

% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'ProtocolTest';
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0, BrainstormDbDir);
% Start a new report
bst_report('Start');

%% ===== CREATE SUBJECT =====
% Create subject
SubjectName = 'Subject01';
[sSubject, iSubject] = db_add_subject(SubjectName);
% Set anatomy
sTemplates = bst_get('AnatomyDefaults');
iTemplate = find(strcmpi('ICBM152_2019', {sTemplates.Name}));
db_set_template(iSubject, sTemplates(iTemplate), 0);
% Re-compute the MNI transformation
bst_process('CallProcess', 'process_mni_affine', [], [], ...
    'subjectname', SubjectName);
% Create downsampled cortex
CortexFile = 'Subject01\tess_cortex_pial_low.mat';
CortexLowFile = tess_downsize(CortexFile, 500, 'reducepatch');
% Create a folder
ConditionName = 'simulate';
iStudySim = db_add_condition(SubjectName, ConditionName);

%% ===== DISPLAY ANATOMY =====
% Get the updated subject structure
sSubject = bst_get('Subject', SubjectName);
% Get MRI file and surface files
MriFile    = sSubject.Anatomy(sSubject.iAnatomy).FileName;
CortexFile = sSubject.Surface(sSubject.iCortex).FileName;
HeadFile   = sSubject.Surface(sSubject.iScalp).FileName;
% Display MRI
hFigMri1 = view_mri(MriFile);
hFigMri2 = view_mri_3d(MriFile, [], [], 'NewFigure');
hFigMri3 = view_mri_slices(MriFile, 'x', 20);
% Close figures
bst_report('Snapshot', hFigMri1, [], 'view_mri');
bst_report('Snapshot', hFigMri2, [], 'view_mri_3d');
bst_report('Snapshot', hFigMri3, [], 'view_mri_slices');
close([hFigMri1 hFigMri2 hFigMri3]);
% Display scalp and cortex
hFigSurf = view_surface(HeadFile);
hFigSurf = view_surface(CortexFile, [], [], hFigSurf);
hFigMriSurf = view_mri(MriFile, CortexFile);
% Figure configuration
iTess = 2;
panel_surface('SetShowSulci',     hFigSurf, iTess, 1);
panel_surface('SetSurfaceColor',  hFigSurf, iTess, [1 0 0]);
panel_surface('SetSurfaceSmooth', hFigSurf, iTess, 0.5, 0);
panel_surface('SetSurfaceTransparency', hFigSurf, iTess, 0.8);
figure_3d('SetStandardView', hFigSurf, 'left');
pause(0.5);
% Close figures
bst_report('Snapshot', hFigSurf, [], 'view_surface');
bst_report('Snapshot', hFigMriSurf, [], 'view_mri + surface');
close([hFigSurf hFigMriSurf]);

%% ===== DATA SIMULATION =====
% Process: Simulate generic signals
nfiles = 10;
sfreq = 1000;
sSimScout = bst_process('CallProcess', 'process_simulate_matrix', [], [], ...
    'subjectname', SubjectName, ...
    'condition',   ConditionName, ...
    'samples',     nfiles * sfreq, ...
    'srate',       sfreq, ...
    'matlab',      ['ferp = 1; fnoise = 60;' 10 ...
                    'Data(1,:) = sin(pi + ferp*2*pi*t);' 10 ...
                    'Data(Data < 0) = 0;' 10 ...
                    'Data(1,:) = Data(1,:) + sin(fnoise*2*pi*t + pi*rand(1)) + 4 * (0.5 - rand(1,length(t)));' ...
                    'Data(2,:) =             sin(fnoise*2*pi*t + pi*rand(1)) + 4 * (0.5 - rand(1,length(t)));']);
% Process: Set channel file
bst_process('CallProcess', 'process_import_channel', sSimScout, [], ...
    'usedefault',   'Colin27: 10-20 19');  % Colin27: 10-20 19
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sSimScout, [], ...
    'sourcespace', 1, ...  % Cortex surface
    'eeg',         3, ...  % OpenMEEG BEM
    'openmeeg',    struct(...
         'BemSelect',    [0, 0, 1], ...
         'BemCond',      [1, 0.0125, 1], ...
         'BemNames',     {{'Scalp', 'Skull', 'Brain'}}, ...
         'BemFiles',     {{}}, ...
         'isAdjoint',    0, ...
         'isAdaptative', 1, ...
         'isSplit',      0, ...
         'SplitLength',  4000));
% Process: Simulate recordings from scouts
sSimData = bst_process('CallProcess', 'process_simulate_recordings', sSimScout, [], ...
    'scouts',      {'Desikan-Killiany', {'lateraloccipital L', 'medialorbitofrontal R'}}, ...
    'savesources', 0);
% Add time markers
DataMat.Events = db_template('event');
DataMat.Events.label   = 'sin+';
DataMat.Events.color   = [0 1 0];
DataMat.Events.times   = round((.5:1:nfiles) * sfreq) ./ sfreq;
DataMat.Events.epochs  = ones(size(DataMat.Events.times));
DataMat.Events.channels = [];
DataMat.Events.notes    = [];
bst_save(file_fullpath(sSimData.FileName), DataMat, 'v6', 1);
% Save as an EGI raw file
RawFile = bst_fullfile(test_dir, 'run01.raw');
export_data(sSimData.FileName, [], RawFile, 'EEG-EGI-RAW');

%% ===== DISPLAY SENSORS =====
% View sensors
hFig = view_surface(HeadFile);
hFig = view_channels(sSimData.ChannelFile, 'EEG', 1, 1, hFig);
% Hide sensors
hFig = view_channels(sSimData.ChannelFile, 'EEG', 0, 0, hFig);
% Edit good/bad channel for current file
gui_edit_channel(sSimData.ChannelFile);
% Unload everything
bst_memory('UnloadAll', 'Forced');

%% ===== ACCESS THE RECORDINGS =====
% Process: Create link to raw file
sFilesRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',    SubjectName, ...
    'datafile',       {RawFile, 'EEG-EGI-RAW'}, ...
    'channelreplace', 0, ...
    'channelalign',   0, ...
    'evtmode',        'value');
% Process: Set channel file
bst_process('CallProcess', 'process_import_channel', sFilesRaw, [], ...
    'usedefault',   'Colin27: 10-20 19');  % Colin27: 10-20 19
% Process: Refine registration
bst_process('CallProcess', 'process_headpoints_refine', sFilesRaw, []);
% Process: Project electrodes on scalp
bst_process('CallProcess', 'process_channel_project', sFilesRaw, []);
% Process: Snapshot: Sensors/MRI registration
bst_process('CallProcess', 'process_snapshot', sFilesRaw, [], ...
    'target',   1, ...  % Sensors/MRI registration
    'modality', 4, ...  % EEG
    'orient',   1, ...  % left
    'Comment',  'MEG/MRI Registration');

%% ===== PRE-PROCESSING =====
% Process: Power spectrum density (Welch)
sFilesPsdBefore = bst_process('CallProcess', 'process_psd', sFilesRaw, [], ...
    'timewindow',  [], ...
    'win_length',  1, ...
    'win_overlap', 50, ...
    'sensortypes', 'EEG', ...
    'edit', struct(...
         'Comment',         'Power', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0));
% Process: Sinusoid removal: 60Hz 120Hz 180Hz 300Hz
sFilesClean = bst_process('CallProcess', 'process_notch', sFilesRaw, [], ...
    'freqlist',    60, ...
    'sensortypes', 'EEG', ...
    'read_all',    0);
% Process: Band-pass:0.5-40Hz
sFilesClean = bst_process('CallProcess', 'process_bandpass', sFilesClean, [], ...
    'sensortypes', 'EEG', ...
    'highpass',    0.5, ...
    'lowpass',     80, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0);

% Process: Set bad channels
bst_process('CallProcess', 'process_channel_setbad', sFilesClean, [], ...
    'sensortypes', 'Fp2');
% Process: Re-reference EEG
bst_process('CallProcess', 'process_eegref', sFilesClean, [], ...
    'eegref',      'AVERAGE', ...
    'sensortypes', 'EEG');

% Process: Power spectrum density (Welch)
sFilesPsdAfter = bst_process('CallProcess', 'process_psd', sFilesClean, [], ...
    'timewindow',  [], ...
    'win_length',  1, ...
    'win_overlap', 50, ...
    'sensortypes', 'EEG', ...
    'edit', struct(...
         'Comment',         'Power', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0));
% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', [sFilesPsdBefore, sFilesPsdAfter], [], ...
    'target',   10, ...  % Frequency spectrum
    'Comment',  'Power spectrum density');

%% ===== ARTIFACT CLEANING =====
% TODO: SSP, ICA, bad segments, other detections, add ECG/EOG to signal

%% ===== IMPORT RECORDINGS =====
% Process: Import MEG/EEG: Events
sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFilesClean, [], ...
    'subjectname', SubjectName, ...
    'condition',   '', ...
    'eventname',   'sin+', ...
    'timewindow',  [], ...
    'epochtime',   [-0.100, 0.500], ...
    'createcond',  0, ...
    'ignoreshort', 1, ...
    'usectfcomp',  1, ...
    'usessp',      1, ...
    'freq',        [], ...
    'baseline',    [-0.1, -0.0017]);
% Display raster plot
hFigRaster = view_erpimage({sFilesEpochs.FileName}, 'erpimage', 'EEG');
bst_report('Snapshot', hFigRaster, sFilesEpochs(1).FileName, 'ERP image');
close(hFigRaster);

%% ===== AVERAGE =====
% Process: Average: By trial group (folder average)
sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
    'avgtype',    5, ...  % By trial groups (folder average)
    'avg_func',   7, ...  % Arithmetic average + Standard error
    'weighted',   0, ...
    'keepevents', 1);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   4, ...  % EEG
    'Comment',    'Evoked response');
% Process: Snapshot: Recordings topography (contact sheet)
bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
    'target',         7, ...  % Recordings topography (contact sheet)
    'modality',       4, ...  % EEG
    'contact_time',   [0, 0.500], ...
    'contact_nimage', 15, ...
    'Comment',        'Evoked response');

%% ===== DISPLAY ERP =====
% View averages
hFigEeg = view_timeseries(sFilesAvg.FileName, 'EEG');
hFigTp1 = view_topography(sFilesAvg.FileName, 'EEG', '2DSensorCap');
hFigTp2 = view_topography(sFilesAvg.FileName, 'EEG', '3DSensorCap');
hFigTp3 = view_topography(sFilesAvg.FileName, 'EEG', '2DDisc');
hFigTp4 = view_topography(sFilesAvg.FileName, 'EEG', '2DLayout');
hFigTp5 = view_topography(sFilesAvg.FileName, 'EEG', '3DElectrodes');
figure_3d('SetStandardView', hFigTp5, 'left');
% Set time: 90ms
panel_time('SetCurrentTime', 0.250);
% Set filters: 40Hz low-pass, no high-pass
panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0);
% View selected sensors
SelectedChannels = {'O1', 'T5'};
bst_figures('SetSelectedRows', SelectedChannels);
hFigSel = view_timeseries(sFilesAvg.FileName, [], SelectedChannels);
% Select time window
figure_timeseries('SetTimeSelectionManual', hFigSel, [0.070, 0.130]);
% Show sensors on 2DSensorCap topography
figure_3d('ViewSensors', hFigTp1, 1, 1);
% Display time contact sheet for a figure
hContactFig = view_contactsheet( hFigTp2, 'time', 'fig', [], 12, [0 0.120] );
% Screen captures
bst_report('Snapshot', hFigSel, sFilesEpochs(1).FileName, 'Selected channels');
bst_report('Snapshot', hFigTp1, sFilesEpochs(1).FileName, '2DSensorCap');
bst_report('Snapshot', hFigTp2, sFilesEpochs(1).FileName, '3DSensorCap');
bst_report('Snapshot', hFigTp3, sFilesEpochs(1).FileName, '2DDisc');
bst_report('Snapshot', hFigTp4, sFilesEpochs(1).FileName, '2DLayout');
bst_report('Snapshot', hFigTp5, sFilesEpochs(1).FileName, '3DElectrodes');
bst_report('Snapshot', hContactFig, sFilesEpochs(1).FileName, 'Contact sheet');
% Colormap tests
bst_colormaps('SetColormapName', 'eeg', 'jet');
bst_colormaps('SetColormapName', 'eeg', 'cmap_rbw');
bst_colormaps('SetColormapAbsolute', 'eeg', 1);
bst_colormaps('SetMaxMode', 'eeg', 'local');
bst_colormaps('SetDisplayColorbar', 'eeg', 0);
bst_colormaps('RestoreDefaults', 'eeg');
gui_edit_channelflag(sFilesAvg.FileName);
bst_memory('UnloadAll', 'Forced');
close(hContactFig);

%% ===== SOURCE ANALYSIS =====
% Copy headmodel from simulation folder to data folder
sHeadmodel = bst_get('HeadModelForStudy', iStudySim);
db_add(sFilesAvg.iStudy, sHeadmodel.FileName);
% No noise model
bst_process('CallProcess', 'process_noisecov', sFilesAvg, [], ...
    'sensortypes',    'EEG', ...
    'target',         1, ...  % Noise covariance
    'identity',       1);     % Identity matrix
% Process: Compute sources [2018]
sAvgSrc = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'sLORETA: EEG', ...
         'InverseMethod',  'minnorm', ...
         'InverseMeasure', 'sloreta', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       0, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'reg', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'fixed', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'EEG'}}));
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sAvgSrc, [], ...
    'target',    8, ...  % Sources (one time)
    'modality',  4, ...  % EEG
    'orient',    3, ...  % top
    'time',      0.230, ...
    'threshold', 60, ...
    'Comment',   'Average sources');

%% ===== SCOUTS =====
% Process: Scouts time series: V1 L
sScout = bst_process('CallProcess', 'process_extract_scout', sAvgSrc, [], ...
    'timewindow',     [], ...
    'scouts',         {'Brodmann', {'V1 L', 'V1 R'}}, ...
    'scoutfunc',      1, ...  % Mean
    'isflip',         1, ...
    'isnorm',         1, ...
    'concatenate',    1, ...
    'save',           1, ...
    'addrowcomment',  1, ...
    'addfilecomment', 1);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sScout, [], ...
    'target',     5, ...  % Recordings time series
    'modality',   4, ...  % EEG
    'Comment',    'Scouts');

% ===== SOURCE ANALYSIS: VOLUME =====
% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesAvg, [], ...
    'sourcespace', 2, ...  % MRI volume
    'volumegrid',  struct(...
         'Method',        'adaptive', ...
         'nLayers',       17, ...
         'Reduction',     3, ...
         'nVerticesInit', 4000, ...
         'Resolution',    0.005, ...
         'FileName',      []), ...
    'eeg',         2);     % 3-shell sphere
% Process: Compute sources [2018]
sAvgSrcVol = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
    'output',  1, ...  % Kernel only: shared
    'inverse', struct(...
         'Comment',        'Dipoles: EEG', ...
         'InverseMethod',  'gls', ...
         'InverseMeasure', 'performance', ...
         'SourceOrient',   {{'free'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'reg', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'rms', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'EEG'}}));
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sAvgSrcVol, [], ...
    'target',    8, ...  % Sources (one time)
    'modality',  1, ...  % MEG (All)
    'orient',    3, ...  % top
    'time',      0, ...
    'threshold', 0, ...
    'Comment',   'Dipole modeling');
% Process: Dipole scanning
sDipScan = bst_process('CallProcess', 'process_dipole_scanning', sAvgSrcVol, [], ...
    'timewindow', [-0.040, 0.100], ...
    'scouts',     {});
% Process: Snapshot: Dipoles
bst_process('CallProcess', 'process_snapshot', sDipScan, [], ...
    'target',    13, ...  % Dipoles
    'orient',    3, ...   % top
    'threshold', 90, ...
    'Comment',   'Dipole scanning');
% Process: Snapshot: Dipoles
bst_process('CallProcess', 'process_snapshot', sDipScan, [], ...
    'target',    13, ...  % Dipoles
    'orient',    1, ...   % left
    'threshold', 90, ...
    'Comment',   'Dipole scanning');

%% ===== TIME-FREQUENCY =====
% Process: Time-frequency (Morlet wavelets)
sFilesTf = bst_process('CallProcess', 'process_timefreq', sFilesEpochs, [], ...
    'sensortypes', 'EEG', ...
    'edit',        struct(...
         'Comment',         'Avg,Power,1-80Hz', ...
         'TimeBands',       [], ...
         'Freqs',           [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80], ...
         'MorletFc',        1, ...
         'MorletFwhmTc',    3, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'average', ...
         'RemoveEvoked',    0, ...
         'SaveKernel',      0), ...
    'normalize',   'none');  % None: Save non-standardized time-frequency maps
% Process: Event-related perturbation (ERS/ERD): [-100ms,-1ms]
sFilesTfNorm = bst_process('CallProcess', 'process_baseline_norm', sFilesTf, [], ...
    'baseline',  [-0.1, -0.001], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 0);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sFilesTfNorm, [], ...
    'target',   14, ...  % Time-frequency maps
    'Comment',  'Not normalized');

%% ===== STOP AND REPORTING =====
ReportFile = bst_report('Save', []);
HtmlFile = bst_report('Export', ReportFile, test_dir);
web(HtmlFile);

% Start Brainstorm with no GUI
% brainstorm nogui
% brainstorm stop

% Stop profiling
% profile viewer
% profiview
% coveragerpt

author: brainstorm-tools


text:
function [est_map,cost,iter] = regionpixel(data,arr,col,map,te,c,d,A1,A2,A3,niter,cal)
% data = source image
% arr, col = pixel position
% map = initial field map (=0) at (arr,col)th pixel
% te = time
% c,d = cos(F_fw*te), sin(F_fw*te)
% A1 = according to Appendix A
% A2 = inv(A1'*A1)*A1'
% A3 = [1 exp(j*F_fw*te)]
% max = maximum for delta_field map
% niter = iteration number

%% initialize

iter = 0;
nsets = size(data,3); % number of te  = 3
delta_map = 11;
if cal == 1
    % cost function at (arr,col)th pixel
    cost = zeros(2,niter*2); % cost(1,:) = field map, cost(2,:) = cost value
end
% source images at (arr,col)th pixel
orig_s = reshape(data(arr,col,:),nsets,1);

%% Iteration code

% while (abs(delta_map)>1 && iter<niter && abs(delta_map)<230)
while (abs(delta_map)>1 && iter<niter)

    iter = iter+1; % number of iteration

    %% Apendix A

    s =  orig_s.*(exp(-j*2*pi*map*te).'); % shat (m is freq)
    shat = [real(s);imag(s)];
    % shat = [s1_R s2_R s3_R s1_I s2_I s3_I]
    % A2 = inv(A1'*A1)*A1'
    rowhat = A2*shat; % rowhat vector = [water_R water_I fat_R fat_I]

   if cal == 1
        %% cost function of Appendix A
        water = rowhat(1) + j*rowhat(2);
        fat = rowhat(3) + j*rowhat(4);
        cost(1,2*iter-1) = map;
        cost(2,2*iter-1) = norm(orig_s - A3*[water;fat].*(exp(j*2*pi*map*te).'));
        % cost(1,:) = field map values
        % cost(2,:) = cost function valuses
   end
    %% Apendix B

    % c_1n = 1, d_1n = 0 for all t_n, c_2n = c, d_2n = d
    f1 = -rowhat(1) -rowhat(3)*c +rowhat(4)*d;
    f2 = -rowhat(2) -rowhat(3)*d -rowhat(4)*c;

    shathat = shat + [f1;f2]; % according to  appendix B

    g = 2*pi*[te';te'].*[f2;-f1];

    B = [g A1];
    y = inv(B'*B)*B'*shathat; % y(1)is delta_fieldmap

    % update the fieldmap = prefieldmap + delta_field map
    delta_map = y(1);
    if abs(delta_map) > 30
        map = map + delta_map/abs(delta_map)*30;
    else
        map = map + delta_map;
    end

    if cal == 1
        % cost function of Appendix B
        cost(1,2*(iter)) = map;
        cost(2,2*(iter)) = norm(orig_s - A3*[water;fat].*(exp(j*2*pi*map*te).'));
    end
end
if cal == 1
    est_map = cost(1,iter*2-1);
else
    est_map = map;
    cost = 0;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% dicomNumSlices
%
% Usage: numSlices = dicomNumSlices;
%
% By Davie
% 2008/02/19
%
% From Bob:
% Most dicom-to-nifti functions (including my niftiFromDicom, which calls
% my dicomLoadAllSeries) guess the number of slices by loading all of the
% files and then counting the slices. For 4d sequences like DTI and fMRI,
% where you have a set of volumes acquired many times, you also need to
% look at the DICOM SliceLocation field, which tells you the slice's
% position in the volume (in physical space units). The number of unique
% slice locations is the 'official' DICOM indication of the number of
% slices.
%
% Why a new function?:
% Number of slices can be computed by dicomLoadAllSeries, but I couldn't
% get it to work for some reason. So here is a simple (but less general)
% hack.
%
% This will take a long time to run (which is why I added printing status)

function numSlices = dicomNumSlices

d = dir('*.dcm*');

sliceLocations = zeros(length(d));

for ii=1:length(d)
    curFile=d(ii).name;
    info=dicominfo(curFile);
    sliceLocations(ii)=info.SliceLocation;
    fprintf('\n Loaded slice %d',ii);
    clear info
end

numSlices = unique(sliceLocations);
numSlices = length(numSlices)-1; % subtract 1, b/c there is a 0 location

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function t=tquant(n, p)
%TQUANT  Quantiles of Student's t distribution
%
%  TQUANT(n, p) is the p-quantile of a t distributed random variable
%  with n degrees of freedom; that is, TQUANT(n, p) is the value below
%  which 100p percent of the t distribution with n degrees of freedom
%  lies.

%  Modified 13-Jul-06
%  Author: Tapio Schneider
%          tapio@gps.caltech.edu

%  References:
%  L. Devroye, 1986: "Non-Uniform Random Variate Generation", Springer
%
%  M. Abramowitz and I. A. Stegun, 1964: "Handbook of Mathematical
%     Functions"
%
%  See also: tcdf.m in the Matlab Statistics Toolbox (evaluates
%     cumulative distribution function of Student's t)

  if (n ~= round(n) | n < 1)
    error('Usage: TQUANT(n,p) - Degrees of freedom n must be positive integer.')
  end

  if (p<0 | p>1)
    error('Usage: TQUANT(n,p) - Probability p must be in [0,1].')
  elseif p == 1
    t   = Inf;
    return
  elseif p == 0
    t   = -Inf;
    return
  end

  if n == 1
    % Cauchy distribution (cf. Devroye [1986, pp. 29 and 450])
    t   = tan(pi*(p-.5));
  elseif p >= 0.5
    % positive t-values (cf. M. Abramowitz and I. A. Stegun [1964,
    % Chapter 26])
    b0  = [0, 1];
    f   = inline('1 - betainc(b, n/2, .5)/2 - p', 'b', 'n', 'p');
    opt = optimset('Display', 'off');
    b   = fzero(f, b0, opt, n, p);
    % old calling sequence (on Windows/Mac with older Matlab versions):
    %b   = fzero(f, b0, eps, 0, n, p);

    t   = sqrt(n/b-n);
  else
    % negative t-values
    b0  = [0, 1];
    f   = inline('betainc(b, n/2, .5)/2 - p', 'b', 'n', 'p');
    %opt = optimset('Display', 'off'); % does not work on Windows/Mac
    %b   = fzero(f, b0, opt, n, p);
    % old calling sequence (for Windows/Mac compatibility):
    b   = fzero(f, b0, eps, 0, n, p);
    t   = -sqrt(n/b-n);
  end

author: Sable

repo: mcbench-benchmarks


text:
function [TW] = Btuph2TW(Btuph)
% Convert power from British thermal units per hour to terawatts.
% Chad A. Greene 2012
TW = Btuph*2.930710702e-13;


text:
%%*****************************************************************************
%% checkdepconst: compute AAt to determine if the
%%             constraint matrices Ak are linearly independent.
%%
%% [At,b,y,idxB,neardepconstr,feasible,AAt] = checkdepconstr(blk,At,b,y,rmdepconstr);
%%
%% rmdepconstr = 1, if want to remove dependent columns in At
%%             = 0, otherwise.
%%
%% idxB = indices of linearly independent columns of original At.
%% neardepconstr = 1 if there is nearly dependent columns in At
%%            = 0, otherwise.
%% Note: the definition of "nearly dependent" is dependent on the
%%       threshold used to determine the small diagonal elements in
%%       the LDLt factorization of A*At.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************************

   function  [At,b,y,idxB,neardepconstr,feasible,AAt] = ...
              checkdepconstr(blk,At,b,y,rmdepconstr);

   global existlowrank printlevel
   global nnzmatold
%%
%% compute AAt
%%
   m = length(b);
   AAt = sparse(m,m);  numdencol = 0; UU = [];
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s');
         m1 = size(At{p,1},2); m2 = m - m1;
         if (m2 > 0)
            if (m2 > 0.3*m); AAt = full(AAt); end
            dd = At{p,3};
            lenn = size(dd,1);
            idxD = [0; find(diff(dd(:,1))); lenn];
            ss = [0, cumsum(pblk{3})];
            if (existlowrank)
               AAt(1:m1,1:m1) = AAt(1:m1,1:m1) + At{p,1}'*At{p,1};
               for k = 1:m2
                  idx = [ss(k)+1 : ss(k+1)];
                  idx2 = [idxD(k)+1: idxD(k+1)];
                  ii = dd(idx2,2)-ss(k); %% undo cumulative indexing
                  jj = dd(idx2,3)-ss(k);
                  len = pblk{3}(k);
                  Dk = spconvert([ii,jj,dd(idx2,4); len,len,0]);
                  tmp = svec(pblk,At{p,2}(:,idx)*Dk*At{p,2}(:,idx)');
                  tmp2 = AAt(1:m1,m1+k) + (tmp'*At{p,1})';
                  AAt(1:m1,m1+k) = tmp2;
                  AAt(m1+k,1:m1) = tmp2';
               end
            end
            DD = spconvert([dd(:,2:4); sum(pblk{3}),sum(pblk{3}),0]);
            VtVD = (At{p,2}'*At{p,2})*DD;
            VtVD2 = VtVD'.*VtVD;
            for k = 1:m2
               idx0 = [ss(k)+1 : ss(k+1)];
               %%tmp = VtVD(idx0,:)' .* VtVD(:,idx0);
               tmp = VtVD2(:,idx0);
               tmp = tmp*ones(length(idx0),1);
               tmp3 = AAt(m1+[1:m2],m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
               AAt(m1+[1:m2],m1+k) = tmp3;
            end
         else
            AAt = AAt + abs(At{p,1})'*abs(At{p,1});
         end
      else
         decolidx = checkdense(At{p,1}');
         %% checkdense punts if the matrix has too many dense columns
         %% in this case we make the whole matrix dense
         if ~isempty(decolidx);
            n2 = size(At{p,1},1);
            dd = ones(n2,1);
            len= length(decolidx);
            dd(decolidx) = zeros(len,1);
            AAt = AAt + (abs(At{p,1})' *spdiags(dd,0,n2,n2)) *abs(At{p,1});
            tmp = At{p,1}(decolidx,:)';
            UU = [UU, tmp];
            numdencol = numdencol + len;
         else
            if (nnz(At{p,1})>0.2*numel(At{p,1}))
               Atp=full(abs(At{p,1}));
            else
               Atp=abs(At{p,1});
            end
            AAt = AAt + Atp'*Atp;
         end
      end
   end
   if (numdencol > 0) & (printlevel)
      fprintf('\n number of dense column in A = %d',numdencol);
   end
   numdencol = size(UU,2);
%%
%%
%%
   feasible = 1; neardepconstr = 0;
   if ~issparse(AAt); AAt = sparse(AAt); end
   nnzmatold = mexnnz(AAt);
   rho = 1e-15;
   diagAAt = diag(AAt);
   mexschurfun(AAt,rho*max(diagAAt,1));
   [L.R,indef,L.perm] = chol(AAt,'vector');
   L.d = full(diag(L.R)).^2;
   if (indef)
      msg = 'AAt is not pos. def.';
      idxB = [1:m];
      neardepconstr = 1;
      if (printlevel); fprintf('\n checkdepconstr: %s',msg); end
      return;
   end
%%
%% find independent rows of A
%%
   dd = zeros(m,1);
   idxB = [1:m]';
   dd(L.perm) = abs(L.d);
   idxN = find(dd < 1e-13*mean(L.d));
   ddB = dd(setdiff([1:m],idxN));
   ddN = dd(idxN);
   if ~isempty(ddN) & ~isempty(ddB) & (min(ddB)/max(ddN) < 10)
      %% no clear separation of elements in dd
      %% do not label constraints as dependent
      idxN = [];
   end
   if ~isempty(idxN)
      neardepconstr = 1;
      if (printlevel)
         fprintf('\n number of nearly dependent constraints = %1.0d',length(idxN));
      end
      if (numdencol==0)
         if (rmdepconstr)
            idxB = setdiff([1:m]',idxN);
            if (printlevel)
               fprintf('\n checkdepconstr: removing dependent constraints...');
            end
            [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);
   	    tol = 1e-8;
            if (resnorm > sqrt(tol))
               idxB = [1:m]';
               neardepconstr = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: basis rows cannot be reliably identified,');
                  fprintf(' abort removing nearly dependent constraints');
               end
               return;
            end
            tmp = W'*b(idxB) - b(idxN);
            nnorm = norm(tmp)/max(1,norm(b));
            if (nnorm > tol)
               feasible = 0;
               if (printlevel)
                  fprintf('\n checkdepconstr: inconsistent constraints exist,');
                  fprintf(' problem is infeasible.');
               end
            else
               feasible = 1;
               for p = 1:size(blk,1)
                  At{p,1} = At{p,1}(:,idxB);
               end
               b = b(idxB);
               y = y(idxB);
               AAt = AAt(idxB,idxB);
            end
	 else
            if (printlevel)
               fprintf('\n To remove these constraints,');
               fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
            end
         end
      else
         if (printlevel)
            fprintf('\n warning: the sparse part of AAt may be nearly singular.');
         end
      end
   end
%%*****************************************************************************
%% findcoeffsub:
%%
%% [W,resnorm] = findcoeffsub(blk,At,idXB,idXN);
%%
%% idXB = indices of independent columns of At.
%% idxN = indices of   dependent columns of At.
%%
%% AB = At(:,idxB); AN = At(:,idxN) = AB*W
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

   function [W,resnorm] = findcoeffsub(blk,At,idxB,idxN);

   AB = []; AN = [];
   for p = 1:size(blk,1)
      AB = [AB; At{p,1}(:,idxB)];
      AN = [AN; At{p,1}(:,idxN)];
   end
   [m,n] = size(AB);
%%
%%-----------------------------------------
%% find W so that AN = AB*W
%%-----------------------------------------
%%
   [L,U,P,Q] = lu(sparse(AB));
   rhs  = P*AN;
   Lhat = L(1:n,:);
   W = Q*( U \ (Lhat \ rhs(1:n,:)));
   resnorm = norm(AN-AB*W,'fro')/max(1,norm(AN,'fro'));
%%*****************************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function chrpak_test026 ( )

%*****************************************************************************80
%
%% TEST026 tests CH_TO_ROT13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST026\n' );
  fprintf ( 1, '  CH_TO_ROT13 "encodes" characters using ROT13\n' );
  fprintf ( 1, '  (and digits using ROT5).\n' );
  fprintf ( 1, '  A second application of the function returns the\n' );
  fprintf ( 1, '  original character.\n' );

  s1 = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ';
  for i = 1 : length ( s1 )
    s2(i) = ch_to_rot13 ( s1(i) );
    s3(i) = ch_to_rot13 ( s2(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '            CH  : %s\n', s1 );
  fprintf ( 1, '      ROT13(CH) : %s\n', s2 );
  fprintf ( 1, 'ROT13(ROT13(CH)): %s\n', s3 );

  s1 = 'CH_TO_ROT13 "encodes" characters using ROT13';
  for i = 1 : length ( s1 )
    s2(i) = ch_to_rot13 ( s1(i) );
    s3(i) = ch_to_rot13 ( s2(i) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '            CH  : %s\n', s1 );
  fprintf ( 1, '      ROT13(CH) : %s\n', s2 );
  fprintf ( 1, 'ROT13(ROT13(CH)): %s\n', s3 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef prtRvUniform < prtRv
    % prtRvUniform  Uniform random variable
    %
    %   RV = prtRvUniform creates a prtRvUniform object with empty
    %   upperBounds and lowerBounds. The upperBounds and lowerBounds must
    %   be set either directly, or by calling the MLE method. upperBounds
    %   and lowerBounds specify the range of the uniform variable.
    %
    %   RV = prtRvUniform(PROPERTY1, VALUE1,...) creates a prtRvUniform
    %   object RV with properties as specified by PROPERTY/VALUE pairs.
    %
    %   A prtRvUniform object inherits all properties from the prtRv class.
    %   In addition, it has the following properties:
    %
    %   upperBounds - 1 x nDims double vector specifying the upper bound of
    %                 the region with uniform density
    %   lowerBounds - 1 x nDims double vector specifying the lower bound of
    %                 the region with uniform density
    %
    %  A prtRvUniform object inherits all methods from the prtRv class.
    %  The MLE  method can be used to estimate the distribution parameters
    %  from data.
    %
    %  Example:
    %
    %  dataSet = prtDataGenUnimodal;        % Load a dataset consisting of
    %                                       % 2 features
    %  dataSet = retainFeatures(dataSet,1); % Retain only the first feature
    %
    %  RV = prtRvUniform;                   % Create a prtRvUniform object
    %  RV = RV.mle(dataSet);                % Compute the bounds
    %                                       % form the data
    %  RV.plotPdf                           % Plot the pdf
    %
    %   See also: prtRv, prtRvMvn, prtRvGmm, prtRvMultinomial,
    %   prtRvVq, prtRvKde

    properties (SetAccess = private)
        name = 'Uniform Random Variable'
        nameAbbreviation = 'RVUnif';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties
        upperBounds  % The lower bounds of the random variable
        lowerBounds  % The upper bounds of the random variable
    end

    properties (Hidden = true, Dependent = true)
        nDimensions
    end

    methods
        % The Constructor
        function R = prtRvUniform(varargin)
            R = constructorInputParse(R,varargin{:});
        end

        function R = mle(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            R.upperBounds = max(X,[],1);
            R.lowerBounds = min(X,[],1);
        end

        function vals = pdf(R,X)

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            X = R.dataInputParse(X); % Basic error checking etc
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)

            isIn = true(size(X,1),1);
            for iDim = 1:size(X,2)
                isIn = isIn & X(:,iDim) >= R.lowerBounds(iDim) & X(:,iDim) <= R.upperBounds(iDim);
            end

            vals = zeros(size(X,1),1);
            vals(isIn) = 1/R.area();
        end

        function vals = logPdf(R,X)
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = log(pdf(R,X));
        end

        function vals = cdf(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            vals = prod(max(bsxfun(@minus,bsxfun(@min,X,R.upperBounds),R.lowerBounds),0),2)./R.area();
        end

        function a = area(R)
            a = prod(R.upperBounds - R.lowerBounds);
        end

        function vals = draw(R,N)
            if nargin < 2
                N = 1;
            end

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = bsxfun(@plus,bsxfun(@times,rand(N,R.nDimensions),R.upperBounds - R.lowerBounds),R.lowerBounds);
        end

        function val = get.nDimensions(R)
            val = length(R.upperBounds);
        end
    end

    methods (Hidden = true)
        function [val, reasonStr] = isValid(R)
            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            badOrder = ~all(R.upperBounds>R.lowerBounds);

            val = ~isempty(R.upperBounds) && ~isempty(R.lowerBounds) && ~badOrder;

            if val
                reasonStr = '';
            else
                badUpper = isempty(R.upperBounds);
                badLower = isempty(R.lowerBounds);

                if badUpper && ~badLower
                    reasonStr = 'because upperBounds has not been set';
                elseif ~badUpper && badLower
                    reasonStr = 'because lowerBounds has not been set';
                elseif badUpper && badLower
                    reasonStr = 'because upperBounds and lowerBounds have not been set';
                elseif badOrder
                    reasonStr = 'because at least one entry of lowerBounds is greater than the corresponding entry of upperBounds';
                else
                    reasonStr = 'because of an unknown reason';
                end
            end
        end

        function val = plotLimits(R)
            [isValid, reasonStr] = R.isValid;
            if isValid

                range = R.upperBounds-R.lowerBounds;

                val = zeros(2*R.nDimensions,1);
                val(1:2:end) = R.lowerBounds-range/10;
                val(2:2:end) = R.upperBounds+range/10;

            else
                error('prtRvUniform:plotLimits','Plotting limits can not be determined for this RV. It is not yet valid %s.',reasonStr)
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% pop_writelocs() - load a EGI EEG file (pop out window if no arguments).
%
% Usage:
%   >> EEG = pop_writelocs(chanstruct);             % a window pops up
%   >> EEG = pop_writelocs(chanstruct, filename, 'key', val, ...);
%
% Inputs:
%   chanstruct     - channel structure. See readlocs()
%   filename       - Electrode location file name
%   'key',val      - same as writelocs()
%
% Author: Arnaud Delorme, CNL / Salk Institute, 17 Dec 2002
%
% See also: writelocs()

% Copyright (C) Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function com = pop_writelocs(chans, filename, varargin);

com = '';
if nargin < 1
   help pop_writelocs;
   return;
end;

if isfield(chans, 'shrink')
    chans = rmfield(chans, 'shrink');
    disp('Warning: shrink factor ignored');
end;

disp('WARNING: ELECTRODE COORDINATES MUST BE WITH NOSE ALONG THE +X DIMENSION TO BE EXPORTED')
disp('         IF NOT, THE EXPORTED FILE COORDINATES MAY BE INNACURATE')

% get infos from readlocs
% -----------------------
[chanformat listcolformat] = readlocs('getinfos');
chanformat(end)    = [];
listcolformat(end) = []; % remove chanedit
chanformat(end)    = [];
listcolformat(end) = []; % remove chanedit
indformat  = [];
for index = 1:length(chanformat),
    if ~isstr(chanformat(index).importformat)
        indformat = [ indformat index ];
    end;
    if isempty(chanformat(index).skipline), chanformat(index).skipline = 0; end;
end;
listtype   = { chanformat(indformat).type };
formatinfo = { chanformat(indformat).importformat };
formatskip = [ chanformat(indformat).skipline ];

%[listtype formatinfo listcolformat formatskip] = readlocs('getinfoswrite');

listtype{end+1} = 'custom';
formatinfo{end+1} = {};
formatskip = [ formatskip 0];

if nargin < 2
   updatefields = [ 'tmpdata = get(gcf, ''userdata'');' ...
                    'tmpobj = findobj(gcf, ''tag'', ''list2'');' ...
                    'set(tmpobj, ''string'', strvcat(tmpdata{2}));' ...
                    'clear tmpobj tmpdata;' ];
   addfieldcom = [ 'tmpdata = get(gcbf, ''userdata'');' ...
                   'tmpobj = findobj(gcf, ''tag'', ''list1'');' ...
                   'tmpdata{2}{end+1} = tmpdata{1}{get(tmpobj, ''value'')};' ...
                   'set(gcbf, ''userdata'', tmpdata);' ...
                   updatefields ];
   rmfieldcom  = [ 'tmpdata = get(gcbf, ''userdata'');' ...
                   'tmpobj = findobj(gcbf, ''tag'', ''list2'');' ...
                   'try, tmpdata{2}(get(tmpobj, ''value'')) = [];' ...
                   '    set(tmpobj, ''value'', 1);' ...
                   'catch, end;' ...
                   'set(gcbf, ''userdata'', tmpdata);' ...
                   updatefields ];
   filetypecom = [ 'tmpdata = get(gcf, ''userdata'');' ...
                   'tmpobj = findobj(gcf, ''tag'', ''formatlist'');' ...
                   'tmpval = get(tmpobj, ''value'');' ...
                   'try, tmpdata{2} = tmpdata{3}{tmpval}; catch, end;' ... %try and catch for custom
                   'set(gcf, ''userdata'', tmpdata);' ...
                   updatefields ...
                   'tmpdata = get(gcf, ''userdata'');' ...
						 'tmpobj1 = findobj(gcf, ''tag'', ''insertcol'');' ... % the lines below
                   'tmpobj2 = findobj(gcf, ''tag'', ''inserttext'');' ... % update the checkbox
                   'try, ' ...                                     % and the edit text box
                   '  if tmpdata{4}(tmpval) == 2,' ...
                   '     set(tmpobj1, ''value'', 1);' ...
                   '  else,' ...
                   '     set(tmpobj1, ''value'', 0);' ...
                   '  end;' ...
                   '  if tmpval == 1,' ... % besa only
                   '     set(tmpobj2, ''string'', ''' int2str(length(chans)) ''');' ...
                   '  else,' ...
                   '     set(tmpobj2, ''string'', '''');' ...
                   '  end;' ...
                   'catch, end;' ... % catch for custom case
                   'tmpobj = findobj(gcf, ''userdata'', ''setfield'');' ...
                   'if tmpval == ' int2str(length(listtype)) ',' ... % disable if non-custom type
                   '   set(tmpobj, ''enable'', ''on'');' ...
                   'else,' ...
                   '   set(tmpobj, ''enable'', ''off'');' ...
                   'end; clear tmpobj tmpobj2 tmpdata tmpval;' ];

   geometry = { [1 1 1] [1 1] [1] [1 1 1] [1 1] [1 1 1] [1 0.3 0.7] [1] [1] };
   listui = { ...
         { 'style' 'text' 'string' 'Filename' } ...
         { 'style' 'edit' 'string' '' 'tag' 'filename' 'horizontalalignment' 'left' } ...
         { 'style' 'pushbutton' 'string' 'Browse' 'callback' ...
           [  '[tmpfile tmppath] = uiputfile(''*'', ''Exporting electrode location file -- pop_writelocs()'');' ...
              'set(findobj(gcbf, ''tag'', ''filename''), ''string'', char([tmppath tmpfile ]));' ...
              'clear tmpfile tmppath;' ] } ...
         { 'style' 'text' 'string' strvcat('Select output file type', ' ', ' ') } ...
         { 'style' 'listbox' 'tag' 'formatlist' 'string' strvcat(listtype) ...
            'value' length(listtype) 'callback' filetypecom } ...
         { 'style' 'text' 'string' 'Select fields to export below' } ...
         { } { 'style' 'pushbutton' 'string' '-> ADD' 'callback' addfieldcom 'userdata' 'setfield' } { } ...
         { 'style' 'listbox' 'tag' 'list1' 'string' strvcat(fieldnames(chans)) 'userdata' 'setfield' } ...
         { 'style' 'listbox' 'tag' 'list2' 'string' '' 'userdata' 'setfield2' } ...
         { } { 'style' 'pushbutton' 'string' 'REMOVE <-' 'callback' rmfieldcom 'userdata' 'setfield' } { } ...
         { 'style' 'text' 'string' 'Insert column names' } ...
         { 'style' 'checkbox' 'tag' 'insertcol' 'value' 1 'userdata' 'setfield' } { } ...
         { 'style' 'text' 'string' 'Enter custom header below' } ...
         { 'style' 'edit' 'userdata' 'setfield' 'tag' 'inserttext' 'horizontalalignment' 'left' 'max' 2 } ...
	};

   inputgui(geometry, listui, 'pophelp(''writelocs'');', ...
      'Exporting electrode location file -- pop_writelocs()', { fieldnames(chans) {} formatinfo formatskip }, 'plot', [1 3 1 1 3 1 1 1 3 ]);
   fig = gcf;

   % set default format
   tmpobj = findobj(fig, 'tag', 'formatlist');
   set(tmpobj, 'value', 6);
   eval(get(tmpobj, 'callback'));

   res = inputgui(geometry, listui, 'pophelp(''writelocs'');', ...
      'Exporting electrode location file -- pop_writelocs()', { listcolformat {} formatinfo formatskip }, fig, [1 3 1 1 3 1 1 1 3 ]);
   if gcf ~= fig, return; end;
   exportfields = get(fig, 'userdata');
   exportfields = exportfields{2};
   close(fig);

   % decode the inputs
   filename = res{1};
   if isempty(filename),
      errordlg2('Error: Empty file name', 'Error');
      return;
   end;
   options = { 'filetype' listtype{res{2}} 'format' exportfields ...
         'header' fastif(res{5}, 'on', 'off') 'customheader' res{6} };
else
	options = varargin;
end;

% generate history
% ----------------
if isempty(inputname(1)) % not a variable name -> probably the structure from pop_chanedit
    writelocs(chans, filename, options{:});
   com = sprintf('pop_writelocs( EEG.chanlocs, ''%s'', %s);', filename, vararg2str(options));
else
    if strcmpi(inputname(1), 'chantmp')
        % do not write file (yet)
        com = sprintf('pop_writelocs( chans, ''%s'', %s);', filename, vararg2str(options));
    else
        writelocs(chans, filename, options{:});
        com = sprintf('pop_writelocs( %s, ''%s'', %s);', inputname(1), filename, vararg2str(options));
    end;
end;


text:
function [fittedCircle, circleNormal] = fitCircle3d(pts)
%FITCIRCLE3D Fit a 3D circle to a set of points.
%
%   [FITTEDCIRCLE, CIRCLENORMAL] = fitCircle3d(PTS)
%
%   Example
%     % points on a 2d circle with noise
%     nop = randi([5 50],1,1);
%     radius = randi([5 25],1,1);
%     points2d = circleToPolygon([0 0 radius], nop);
%     points2d(1,:) = [];
%     points2d = points2d + rand(size(points2d));
%     points2d(:,3)=rand(length(nop),1);
%     % apply random rotation and translation
%     [theta, phi] = randomAngle3d;
%     theta = rad2deg(theta);
%     phi = rad2deg(phi);
%     tfm = eulerAnglesToRotation3d(phi, theta, 0);
%     trans = randi([-250 250],3,1);
%     tfm(1:3,4)=trans;
%     points3d = awgn(transformPoint3d(points2d, tfm),1);
%     % fit 3d circle
%     [fittedCircle, circleNormal] = fitCircle3d(points3d);
%     % plot 3d points and 3d circle
%     figure('Color','w'); hold on; axis equal tight; view(3);
%     xlabel('X');ylabel('Y');zlabel('Z');
%     drawPoint3d(points3d)
%     drawCircle3d(fittedCircle, 'k')
%     drawVector3d(fittedCircle(1:3), circleNormal*fittedCircle(4))
%
%   See also
%   circle3dOrigin, circle3dPosition, circle3dPoint, intersectPlaneSphere
%   drawCircle3d, drawCircleArc3d, drawEllipse3d
%
% ------
% Authors: oqilipo, David Legland
% created: 2017-05-09

% Mean of all points
meanPoint = mean(pts,1);

% Center points by subtracting the meanPoint
centeredPoints = pts - repmat(meanPoint,size(pts,1),1);

% Project 3D data to a plane
[~,~,V]=svd(centeredPoints);
tfmPoints = transformPoint3d(centeredPoints, V');

% Fit a circle to the points in the xy-plane
circleParamter = CircleFitByTaubin(tfmPoints(:,1:2));
center2d = circleParamter(1:2);
radius=circleParamter(3);
center3d = transformPoint3d([center2d, 0], [inv(V'), meanPoint']);
circleNormal = V(:,3)';
[theta, phi, ~] = cart2sph2(circleNormal);
fittedCircle = [center3d radius rad2deg(theta) rad2deg(phi) 0];

end

% Circle Fit (Taubin method)
% version 1.0 (2.24 KB) by Nikolai Chernov
% http://www.mathworks.com/matlabcentral/fileexchange/22678
function Par = CircleFitByTaubin(XY)

%--------------------------------------------------------------------------
%
%     Circle fit by Taubin
%      G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
%                  Space Curves Defined By Implicit Equations, With
%                  Applications To Edge And Range Image Segmentation",
%      IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
%
%     Input:  XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%     Note: this fit does not use built-in matrix functions (except "mean"),
%           so it can be easily programmed in any programming language
%
%--------------------------------------------------------------------------

n = size(XY,1);      % number of data points

centroid = mean(XY);   % the centroid of the data set

%     computing moments (note: all moments will be normed, i.e. divided by n)

Mxx = 0; Myy = 0; Mxy = 0; Mxz = 0; Myz = 0; Mzz = 0;

for i=1:n
    Xi = XY(i,1) - centroid(1);  %  centering data
    Yi = XY(i,2) - centroid(2);  %  centering data
    Zi = Xi*Xi + Yi*Yi;
    Mxy = Mxy + Xi*Yi;
    Mxx = Mxx + Xi*Xi;
    Myy = Myy + Yi*Yi;
    Mxz = Mxz + Xi*Zi;
    Myz = Myz + Yi*Zi;
    Mzz = Mzz + Zi*Zi;
end

Mxx = Mxx/n;
Myy = Myy/n;
Mxy = Mxy/n;
Mxz = Mxz/n;
Myz = Myz/n;
Mzz = Mzz/n;

%    computing the coefficients of the characteristic polynomial

Mz = Mxx + Myy;
Cov_xy = Mxx*Myy - Mxy*Mxy;
A3 = 4*Mz;
A2 = -3*Mz*Mz - Mzz;
A1 = Mzz*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz - Mz*Mz*Mz;
A0 = Mxz*Mxz*Myy + Myz*Myz*Mxx - Mzz*Cov_xy - 2*Mxz*Myz*Mxy + Mz*Mz*Cov_xy;
A22 = A2 + A2;
A33 = A3 + A3 + A3;

xnew = 0;
ynew = 1e+20;
epsilon = 1e-12;
IterMax = 20;

% Newton's method starting at x=0

for iter=1:IterMax
    yold = ynew;
    ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));
    if abs(ynew) > abs(yold)
       disp('Newton-Taubin goes wrong direction: |ynew| > |yold|');
       xnew = 0;
       break;
    end
    Dy = A1 + xnew*(A22 + xnew*A33);
    xold = xnew;
    xnew = xold - ynew/Dy;
    if (abs((xnew-xold)/xnew) < epsilon), break, end
    if (iter >= IterMax)
        disp('Newton-Taubin will not converge');
        xnew = 0;
    end
    if (xnew<0.)
        fprintf(1,'Newton-Taubin negative root:  x=%f\n',xnew);
        xnew = 0;
    end
end

%  computing the circle parameters

DET = xnew*xnew - xnew*Mz + Cov_xy;
Center = [Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy]/DET/2;

Par = [Center+centroid , sqrt(Center*Center'+Mz)];

end    %    CircleFitByTaubin


text:
function handles = initalizeFields(handles)
%INITIALIZEFIELDS    Initialize fonts when CHEBGUI is started.
%
%  Calling sequence:
%       HANDLES = INITIALIZEFIELDS(HANDLES)
%  where
%       HANDLES:    A Matlab handle of the CHEBGUI figure.
%

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Set the default font and font size used:
myFont = 'Monospaced';
mySize = 12;

% Obtain a list of all GUI elements.
names = fieldnames(handles);
% Find what elements are input and button elements:
inputLocs = strfind(names, 'input_');
buttonLocs = strfind(names, 'button_');
% Also want the same font for the iter_list information box.
iterListLoc = strfind(names, 'iter_list');
% And what fields are text fields:
textLocs = strfind(names, 'text_');
% Combine all the locations of elements whose font we wish to change
monospacedLocs = strcat(inputLocs, iterListLoc);
variableWidthLocs = strcat(buttonLocs, textLocs);

% Loop through the elements we want to specify the font of.
for fieldCounter = 1:length(inputLocs)
    if ( ~isempty(monospacedLocs{fieldCounter}) )
        % Access the field values dynamically using the .( ) call.
        set(handles.(names{fieldCounter}), 'FontName', myFont);
        set(handles.(names{fieldCounter}), 'FontSize', mySize + 1);
    end

    if ( ~isempty(variableWidthLocs{fieldCounter}) )
        set(handles.(names{fieldCounter}), 'FontSize', mySize);
    end
end

% Set the string for popup-menu for the choice of plots:
set(handles.popupmenu_bottomFig,'String', ...
    {'Convergence of Newton iteration', ...
    'Coefficients of the solution'}, ...
    'fontsize', handles.fontsizePanels);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function drawAhmPnt(MapFig, Lmk, color, MapOpt)

% DRAWAHMPNT  Draw anchored Homogeneous point landmark in MapFig.

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

global Map

posOffset = [0;0;.2];

% transform to Euclidean
[x,P] = propagateUncertainty(       ...
    Map.x(Lmk.state.r),             ...
    Map.P(Lmk.state.r,Lmk.state.r), ...
    @ahm2euc);

% draw
drawPnt    (MapFig.Lmk(Lmk.lmk).mean,    x,    color.mean)
if MapOpt.showEllip
    drawEllipse(MapFig.Lmk(Lmk.lmk).ellipse, x, P, color.ellip)
end
if MapOpt.showLmkId
    drawLabel  (MapFig.Lmk(Lmk.lmk).label,   x+posOffset, num2str(Lmk.id))
end

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


text:
%2D level 0
%tic
load elements.mat
fprintf('2D triangulation:  %d elements, ',size(elements,1));
fprintf('%d nodes ',max(max(elements)));
[element2edges, edge2nodes]=getEdges(elements);
edge2elements=entryInWhichRows(element2edges);
fprintf('--> %d edges numbered. \n',size(edge2nodes,1));
%toc

%2D level 1
%tic
load elements_1.mat
fprintf('2D triangulation:  %d elements, ',size(elements,1));
fprintf('%d nodes ',max(max(elements)));
[element2edges, edge2nodes]=getEdges(elements);
edge2elements=entryInWhichRows(element2edges);
fprintf('--> %d edges numbered. \n',size(edge2nodes,1));
%toc

%2D level 9
tic
load elements_9.mat
fprintf('2D triangulation:  %d elements, ',size(elements,1));
fprintf('%d nodes ',max(max(elements)));
tic
[element2edges, edge2nodes]=getEdges(elements);
edge2elements=entryInWhichRows(element2edges);
fprintf('--> %d edges numbered. \n',size(edge2nodes,1));
toc

%3D
%tic
elements=[1 2 3 4; 1 2 3 5];
fprintf('3D triangulation:  %d elements, ',size(elements,1));
fprintf('%d nodes ',max(max(elements)));
[element2edges, edge2nodes]=getEdges(elements);
edge2elements=entryInWhichRows(element2edges);
fprintf('--> %d edges numbered. \n',size(edge2nodes,1));
%toc

%3D
tic
load elements3D
fprintf('3D triangulation:  %d elements, ',size(elements,1));
fprintf('%d nodes ',max(max(elements)));
[element2edges, edge2nodes]=getEdges(elements);
edge2elements=entryInWhichRows(element2edges);
fprintf('--> %d edges numbered. \n',size(edge2nodes,1));
toc


text:
function name = p22_name ( )

%*****************************************************************************80
%
%% P22_NAME returns the name of problem 22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'BallMonomial';

  return
end


text:
% rgb.m: translates a colour from multiple formats into matlab colour format
% type 'rgb demo' to get started
%
% [matlabcolor]=rgb(col)
% matlab colors are in the format [R G B]
%
% if 'col' is a string, it is interpreted as
%
% 	[[modifier] descriptor] colour_name
%
% where
%		modifier is one of   (slightly, normal, very, extremely)
%		descriptor is one of (light/pale, normal, dark)
%		colorname is a name of a colour
%			(type 'rgb list' or 'rgb demo' to see them all)
%
% if 'col' is an integer between 0 and &HFFFFFF inclusive,
% it is interpreted as a double word RGB value in the form
% [0][R][G][B]
%
% if 'col' is a negative integer between -1 and -&HFFFFFF
% inclusive, it is interpreted as the complement of a double
% word RGB value in the form [0][B][G][R]
%
% if 'col' is a string of the form 'qbX' or 'qbXX' where X
% is a digit then the number part is interpreted as a qbasic
% color
%
% if 'col' is one of {k,w,r,g,b,y,m,c} a sensible result is
% returned
%
% if 'col' is already in matlab format, it is unchanged

%	VERSION:	06/06/2002
%	AUTHOR:		ben mitch
%	CONTACT:	footballpitch@theconnexion.zzn.com
%	WWW:		www.benmitch.co.uk\matlab (not yet)
%	LOCATION:	figures\colors\

function out=rgb(in)

if isa(in,'char') & length(in)>2 & length(in)<5 & strcmpi('qb',in(1:2))
	out=qbcolor(sscanf(in(3:end),'%i'));
elseif isa(in,'char') & length(in)==1
	out=translatecolorchar(in);
elseif isa(in,'char')
	if strcmp(in,'demo') rgb_demo; return; end
	if strcmp(in,'list') rgb_list; return; end
	out=translatecolorstring(in);
elseif isa(in,'double') & size(in,1)==1 & size(in,2)==1 & abs(in)<16777216
	out=translatecolorRGB(in);
elseif isa(in,'double') & size(in,1)==1 & size(in,2)==3
	out=in;
else
	warning('Unrecognised color format, black assumed');
	out=[0 0 0];
end

function out=translatecolorchar(in)
switch(in)
case 'k', out=[0 0 0];
case 'w', out=[1 1 1];
case 'r', out=[1 0 0];
case 'g', out=[0 1 0];
case 'b', out=[0 0 1];
case 'y', out=[1 1 0];
case 'm', out=[1 0 1];
case 'c', out=[0 1 1];
otherwise
	warning(['Unrecognised colour "' in '", black assumed'])
	out=[0 0 0];
	return;
end

function out=translatecolorstring(in)
args.tokens=rgb_parse(in);
args.N=length(args.tokens);
if args.N>3 warning('Too many words in color description, any more than 3 will be ignored'); end
while(args.N<3)
	args.tokens=[{'normal'};args.tokens];
	args.N=args.N+1;
end

cols=get_cols;
col=[];
for n=1:size(cols,1)
	names=cols{n,1};
	for m=1:length(names)
		if strcmp(args.tokens{3},names{m}) col=cols{n,2}; break; end
	end
	if ~isempty(col) break; end
end

if isempty(col)
	warning(['Unrecognised colour "' args.tokens{3} '", black assumed'])
	out=[0 0 0];
	return;
end

switch args.tokens{1}
case 'slightly', fac=0.75;
case 'normal', fac=0.5;
case 'very', fac=0.25;
case 'extremely', fac=0.125;
otherwise
	warning(['Unrecognised modifier "' args.tokens{1} '", normal assumed'])
	fac=0.5;
end

switch args.tokens{2}
case {'light','pale'}, out=1-(1-col)*fac;
case 'normal', out=col;
case 'dark', out=col*fac;
otherwise
	warning(['Unrecognised descriptor "' args.tokens{2} '", normal assumed'])
	out=col;
end

function out=translatecolorRGB(in)

BGR=0;
if in<0
	in=-in;
	BGR=1;
end

b=bytes4(in);
if BGR out=b(4:-1:2); else out=b(2:4); end

function out=qbcolor(in)

% rgb value from basic colour code
% 0-7 are normal, 8-15 are bright
% 0 - black
% 1 - red,  2 - green,   3 - blue
% 4 - cyan, 5 - magenta, 6 - yellow
% 7 - white

bright=0.5;
if in>7 in=in-8; bright=1; end

switch in
case 0, rgb=[0 0 0];
case 1, rgb=[1 0 0];
case 2, rgb=[0 1 0];
case 3, rgb=[0 0 1];
case 4, rgb=[0 1 1];
case 5, rgb=[1 0 1];
case 6, rgb=[1 1 0];
case 7, rgb=[1 1 1];
otherwise
	warning('Unrecognised QBasic color, black assumed');
	out=[0 0 0];
	return;
end

out=rgb*bright;

function tokens=rgb_parse(str)

% parse string to obtain all tokens
% quoted strings count as single tokens

inquotes=0;
intoken=0;
pos=1;
l=length(str);
st=0;
ed=0;
token='';
tab=char(9);
tokens=cell(0);
while(pos<=l)
	ch=str(pos);
	if inquotes
		if ch=='"'
			inquotes=0;
			tokens={tokens{:} token};
		else
			token=[token ch];
		end
	elseif intoken
		if ch==' ' | ch==tab
			intoken=0;
			tokens={tokens{:} token};
		elseif ch=='"'
			error(['Quote misplace in <' str '>']);
		else
			token=[token ch];
		end
	else
		if ch==' ' | ch==tab
			% do nothing
		elseif ch=='"'
			token='';
			inquotes=1;
		else
			token=ch;
			intoken=1;
		end
	end
	pos=pos+1;
end

if intoken tokens={tokens{:} token}; end
if inquotes error(['Unpaired quotes in <' str '>']); end

tokens=tokens';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEMO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rgb_demo

figure(1)
clf
cols = get_cols;
cols = {cols{:,1}}';
cols = { cols{:}, ...
	'k', ...
	'r', ...
	'g', ...
	'b', ...
	'y', ...
	'm', ...
	'c', ...
	'w', ...
	'', ...
	'extremely dark green', ...
	'very dark green', ...
	'dark green', ...
	'slightly dark green', ...
	'green', ...
	'slightly pale green', ...
	'pale green', ...
	'very pale green', ...
	'extremely pale green', ...
};

height=9;
x=0;
y=0;
for n=1:length(cols)
	rect(x,y,cols{n})
	y=y+1;
	if y==height
		x=x+2;
		y=0;
	end
end
if y==0 x=x-2; end
axis([0 (x+2) 0 height])
title('names on different rows are alternates')

function rect(x,y,col)
if isempty(col) return; end
r=rectangle('position',[x+0.1 y+0.1 1.8 0.8]);
col_=col;
if iscell(col) col=col{1}; end
colrgb=rgb(col);
if strcmp(col(1),'u') & length(col)==2
	t=text(x+1,y+0.5,{'unnamed',['colour (' col(2) ')']});
	set(r,'facecolor',colrgb);
else
	t=text(x+1,y+0.5,col_);
	set(r,'facecolor',colrgb);
	if sum(colrgb)<1.5 set(t,'color',[1 1 1]); end
end
set(t,'horizontalalignment','center')
set(t,'fontsize',10)

function rgb_list
cols=get_cols;
disp(' ')
for n=1:size(cols,1)
	code=cols{n,2};
	str=cols{n,1};
	str_=[];
	for m=1:length(str)
		str_=[str_ str{m} ', '];
	end
	str_=str_(1:end-2);
	if strcmp(str_(1),'u') & length(str_)==2
		str_=['* (' str_(2) ')'];
	end
	disp(['  [' sprintf('%.1f  %.1f  %.1f',code) '] - ' str_])
end
disp([10 '* colours marked thus are not named - if you know their' 10 '  designation, or if you feel sure a colour is mis-named,' 10 '  email me (address via help) or comment at' 10 '  www.mathworks.com/matlabcentral - "rgb demo" to see them' 10])

function cols=get_cols

cols={
	'black', [0 0 0]; ...
	'navy', [0 0 0.5]; ...
	'blue', [0 0 1]; ...
	'u1', [0 0.5 0]; ...
	{'teal','turquoise'}, [0 0.5 0.5]; ...
	'slateblue', [0 0.5 1]; ...
	{'green','lime'}, [0 1 0]; ...
	'springgreen', [0 1 0.5]; ...
	{'cyan','aqua'}, [0 1 1]; ...
	'maroon', [0.5 0 0]; ...
	'purple', [0.5 0 0.5]; ...
	'u2', [0.5 0 1]; ...
	'olive', [0.5 0.5 0]; ...
	{'gray','grey'}, [0.5 0.5 0.5]; ...
	'u3', [0.5 0.5 1]; ...
	{'mediumspringgreen','chartreuse'}, [0.5 1 0]; ...
	'u4', [0.5 1 0.5]; ...
	'sky', [0.5 1 1]; ...
	'red', [1 0 0]; ...
	'u5', [1 0 0.5]; ...
	{'magenta','fuchsia'}, [1 0 1]; ...
	'orange', [1 0.5 0]; ...
	'u6', [1 0.5 0.5]; ...
	'u7', [1 0.5 1]; ...
	'yellow', [1 1 0]; ...
	'u8', [1 1 0.5]; ...
	'white', [1 1 1]; ...
	};

for n=1:size(cols,1)
	if ~iscell(cols{n,1}) cols{n,1}={cols{n,1}}; end
end

% converts a DWORD into a four byte row vector
function out=bytes4(in)

out=[0 0 0 0];
if in<0 | in>(2^32-1)
	warning('DWORD out of range, zero assumed');
	return;
end

N=4;
while(in>0)
	out(N)=mod(in,256);
	in=(in-out(N))/256;
	N=N-1;
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function hmm_out = logisticMarginaliseHMM(hmm_in,iY)
% This function takes a multiple logistic strcuture HMM and returns the
% equivalent HMM with only the marginalised component iY of the
% logisticYdim components.

Xfulldim = hmm_in.train.ndim;
Xdim = Xfulldim - hmm_in.train.logisticYdim;
Xfulldim_new=Xdim+1;
select_vec = [1:Xdim,Xdim+iY];
hmm_out=hmm_in;
hmm_out.train.origlogisticYdim = hmm_in.train.logisticYdim;

%update S and options:
hmm_out.train.S=hmm_out.train.S(1:Xfulldim_new,1:Xfulldim_new);
hmm_out.train.Sind=hmm_out.train.Sind(1:Xfulldim_new,1:Xfulldim_new);
hmm_out.train.logisticYdim = 1;
hmm_out.train.ndim = Xfulldim_new;

%update W:
for st=1:hmm_in.train.K
    hmm_out.state(st).W.Mu_W = hmm_out.state(st).W.Mu_W(select_vec,select_vec);
    hmm_out.state(st).W.S_W = hmm_out.state(st).W.S_W(select_vec,select_vec,select_vec);
    hmm_out.state(st).W.iS_W = hmm_out.state(st).W.iS_W(select_vec,select_vec,select_vec);
end

%update alpha and sigma
for st=1:hmm_out.train.K
    %recall that for logistic setups, alpha has dimension Xdim x logisticYdim
    hmm_out.state(st).alpha.Gam_rate = hmm_out.state(st).alpha.Gam_rate(:,iY);
end

%update psi:
if isfield(hmm_in,'psi')
    if size(hmm_in.psi,2)==hmm_in.train.logisticYdim
        %has psi been calculated specific to this variable:
        hmm_out.psi = hmm_in.psi(:,iY);
    else
        hmm_out = rmfield(hmm_out,'psi');
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [isinside,pt,coord]=linextriangle(p0,p1,plane)
%  [isinside,pt,coord]=linextriangle(p0,p1,plane)
%
%  calculate the intersection of a 3d line (passing two points)
%  with a plane (determined by 3 points)
%
%  author: Qianqian Fang <q.fang at neu.edu>
%  date: 12/12/2008
%
% parameters:
%      p0: a 3d point in form of (x,y,z)
%      p1: another 3d point in form of (x,y,z), p0 and p1 determins the line
%      plane: a 3x3 matrix, each row is a 3d point in form of (x,y,z)
%             this is used to define a plane
% outputs:
%      isinside: a boolean variable, 1 for the intersection is within the
%               3d triangle determined by the 3 points in plane; 0 is outside
%      pt: the coordinates of the intersection pint
%      coord: 1x3 vector, if isinside=1, coord will record the barycentric
%          coordinate for the intersection point within the triangle;
%          otherwise it will be all zeros.
%
% for degenerated lines or triangles, this will stop
%
% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
% this function is part of "metch" toobox, see COPYING for license

[a,b,c,d]=getplanefrom3pt(plane);

if(a*a+b*b+c*c==0.0)
        error('degenerated plane');
end

dl_n=sum([a b c].*(p1-p0));

if(dl_n==0.0)
        error('degenerated line');
end

% solve for the intersection point
t=-(a*p0(1)+b*p0(2)+c*p0(3)+d)/dl_n;
pt=p0+(p1-p0)*t;

dist=sum(abs(diff(plane)));
[md,imax]=sort(dist);
if(md(2)==0.0)
        error('degenerated triangle');
end
goodidx=imax(2:end);

ptproj=pt(goodidx);
mat0=[plane(:,goodidx)',ptproj';1 1 1 1];

isinside=0;
coord=[0 0 0];

det1=det(mat0(:,[4 2 3],:));
det2=det(mat0(:,[1 4 3],:));
if(det1*det2<0)
        return;
end
det3=det(mat0(:,[1 2 4],:));
if(det2*det3<0)
        return;
end
if(det1*det3<0)
        return;
end
isinside=1;
det0=det(mat0(:,1:3));

coord=[det1 det2 det3]/det0;

author: fangq


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the logistic
%  regression exercise. You will need to complete the following functions
%  in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = load('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);

%% ==================== Part 1: Plotting ====================
%  We start the exercise by first plotting the data to understand the
%  the problem we are working with.

fprintf(['Plotting data with + indicating (y = 1) examples and o ' ...
         'indicating (y = 0) examples.\n']);

plotData(X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============ Part 2: Compute Cost and Gradient ============
%  In this part of the exercise, you will implement the cost and gradient
%  for logistic regression. You neeed to complete the code in
%  costFunction.m

%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(X);

% Add intercept term to x and X_test
X = [ones(m, 1) X];

% Initialize fitting parameters
initial_theta = zeros(n + 1, 1);

% Compute and display initial cost and gradient
[cost, grad] = costFunction(initial_theta, X, y);

fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', grad);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 3: Optimizing using fminunc  =============
%  In this exercise, you will use a built-in function (fminunc) to find the
%  optimal parameters theta.

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);

%  Run fminunc to obtain the optimal theta
%  This function will return theta and the cost
[theta, cost] = ...
	fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);

% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('theta: \n');
fprintf(' %f \n', theta);

% Plot Boundary
plotDecisionBoundary(theta, X, y);

% Put some labels
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 4: Predict and Accuracies ==============
%  After learning the parameters, you'll like to use it to predict the outcomes
%  on unseen data. In this part, you will use the logistic regression model
%  to predict the probability that a student with score 45 on exam 1 and
%  score 85 on exam 2 will be admitted.
%
%  Furthermore, you will compute the training and test set accuracies of
%  our model.
%
%  Your task is to complete the code in predict.m

%  Predict probability for a student with score 45 on exam 1
%  and score 85 on exam 2

prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n\n'], prob);

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 3 Assignments/Logistic Regression and Regularization/mlclass-ex2/ex2.m


text:
    function y = x5fun(x1,x2,x3,x4,x6,r1,r4,r5)
             yn = -x3^2-r4^2+r1^2+r5^2-x4^2-2*x2*x6+2*x4*x6;
             yd = 2*(x1-x3);
             y = yn/yd;


text:
function fx = p21_fun ( n, x )

%*****************************************************************************80
%
%% P21_FUN evaluates the integrand for problem 21.
%
%  Interval:
%
%    0 <= X <= 1
%
%  Integrand:
%
%       ( sech (   10.0 * ( x - 0.2 ) ) )^2
%     + ( sech (  100.0 * ( x - 0.4 ) ) )^4
%     + ( sech ( 1000.0 * ( x - 0.6 ) ) )^6
%
%  Exact Integral:
%
%    ( 1 + tanh ( 8 ) * tanh ( 2 ) ) / 10.0 + 2 / 150 + 2 / 1875
%
%  Approximate Integral (20 digits):
%
%    0.21080273631018169851...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner,
%    Comparison of Numerical Quadrature Formulas,
%    in Mathematical Software, edited by John R Rice,
%    Academic Press, 1971.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx = ...
      ( sech (   10.0 * ( x - 0.2 ) ) ).^2 ...
    + ( sech (  100.0 * ( x - 0.4 ) ) ).^4 ...
    + ( sech ( 1000.0 * ( x - 0.6 ) ) ).^6;

  return
end


text:
function D = hmm_kl (hmm_p,hmm_q)
% Computes Kullback-Leibler divergence between two Hidden Markov Model
% distributions, through an approximation (an upper bound) as proposed in
% M. Do (2003). IEEE Signal Processing Letters 10
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

K = length(hmm_p.Pi);
if K~=length(hmm_q.Pi)
    error(['The two HMMs must have the same number of states, ' ...
        'and their order must correspond'])
end
if (hmm_p.train.order ~= hmm_q.train.order) || ...
        (~strcmpi(hmm_p.train.covtype,hmm_q.train.covtype)) || ...
        (length(hmm_p.train.embeddedlags) ~= length(hmm_q.train.embeddedlags)) || ...
        (any(hmm_p.train.embeddedlags ~= hmm_q.train.embeddedlags)) || ...
        (hmm_p.train.zeromean ~= hmm_q.train.zeromean)
   error('The state configuration of the two HMMs must be identical')
end
hmm = hmm_p; setstateoptions;
if isfield(hmm_p.state(1),'W')
    ndim = size(hmm_p.state(1).W.Mu_W,2);
else
    ndim = size(hmm_p.state(1).Omega.Gam_rate,2);
end
S = hmm.train.S==1;
regressed = sum(S,1)>0;

D = 0;
if hmm_p.train.id_mixture
   hmm_p.P = 1/K*ones(K); hmm_q.P = 1/K*ones(K);
end
nu = compute_nu (hmm_p.Pi,hmm_p.P); % weight vector

% Non-state specific stuff
switch train.covtype
    case {'uniquediag','shareddiag'}
        for n = 1:ndim
            if ~regressed(n), continue; end
            D = D + gamma_kl(hmm_p.Omega.Gam_shape,hmm_q.Omega.Gam_shape, ...
                hmm_p.Omega.Gam_rate(n),hmm_q.Omega.Gam_rate(n));
        end
    case {'uniquefull','sharedfull'}
        D = D + wishart_kl(hmm_p.Omega.Gam_rate(regressed,regressed),...
            hmm_q.Omega.Gam_rate(regressed,regressed), ...
            hmm_p.Omega.Gam_shape,hmm_q.Omega.Gam_shape);
    case 'pca'
        D = D + gamma_kl(hmm_p.Omega.Gam_shape,hmm_q.Omega.Gam_shape, ...
                hmm_p.Omega.Gam_rate,hmm_q.Omega.Gam_rate);
end

% State specific stuff
for k = 1:K

    % Trans probabilities
    kk = hmm.train.Pstructure(k,:);
    D = D + nu(k) * dirichlet_kl(hmm_p.Dir2d_alpha(k,kk),hmm_q.prior.Dir2d_alpha(k,kk));

    % State distribution
    hs = hmm_p.state(k);
    hs0 = hmm_q.state(k);

    if ~isempty(hs.W.Mu_W)
        if train.uniqueAR || ndim==1
            if train.uniqueAR || ndim==1
                D = D + nu(k) * gauss_kl(hs.W.Mu_W, hs0.W.Mu_W, hs.W.S_W, hs0.W.S_W);
            else
                D = D + nu(k) * gauss_kl(hs.W.Mu_W, hs0.W.Mu_W, ...
                    permute(hs.W.S_W,[2 3 1]), permute(hs0.W.S_W,[2 3 1]));
            end
        elseif strcmp(train.covtype,'diag') || ...
                strcmp(train.covtype,'uniquediag') || strcmp(train.covtype,'shareddiag') || ...
                strcmp(train.covtype,'pca')
            for n = 1:ndim
                D = D + nu(k) * gauss_kl(hs.W.Mu_W(Sind(:,n),n),hs0.W.Mu_W(Sind(:,n),n), ...
                    permute(hs.W.S_W(n,Sind(:,n),Sind(:,n)),[2 3 1]),...
                    permute(hs0.W.S_W(n,Sind(:,n),Sind(:,n)),[2 3 1]));
            end
        else % full or sharedfull
            mu_w = hs.W.Mu_W';
            mu_w = mu_w(:);
            mu_w0 = hs0.W.Mu_W';
            mu_w0 = mu_w0(:);
            D = D + nu(k) * gauss_kl(mu_w,mu_w0, hs.W.S_W, hs0.W.S_W);
        end
    end

    switch train.covtype
        case 'diag'
            for n=1:ndim
                if ~regressed(n), continue; end
                D = D + nu(k) * gamma_kl(hs.Omega.Gam_shape,hs0.Omega.Gam_shape, ...
                    hs.Omega.Gam_rate(n),hs0.Omega.Gam_rate(n));
            end
        case 'full'
            try
                D = D + nu(k) * wishart_kl(hs.Omega.Gam_rate(regressed,regressed),...
                    hs0.Omega.Gam_rate(regressed,regressed), ...
                    hs.Omega.Gam_shape,hs0.Omega.Gam_shape);
            catch
                error(['Error computing kullback-leibler divergence of the cov matrix - ' ...
                    'Something strange with the data?'])
            end
    end

    if ~isempty(orders) && ~train.uniqueAR && ndim>1
        for n1=1:ndim
            for n2=1:ndim
                if (train.symmetricprior && n2<n1) || S(n1,n2)==0, continue; end
                D = D + nu(k) * gamma_kl(hs.sigma.Gam_shape(n1,n2),hs0.sigma.Gam_shape(n1,n2), ...
                    hs.sigma.Gam_rate(n1,n2),hs0.sigma.Gam_rate(n1,n2));
            end
        end
    end
    if ~isempty(orders)
        for i=1:length(orders)
            D = D + nu(k) * gamma_kl(hs.alpha.Gam_shape,hs0.alpha.Gam_shape, ...
                hs.alpha.Gam_rate(i),hs0.alpha.Gam_rate(i));
        end
    end
end

end

function nu = compute_nu (Pi,P)
eps = 1e-6;
nu = Pi * P;
while true
    nu0 = nu;
    nu = nu * P;
    if mean(nu(:)-nu0(:))<eps, break; end
end
end


text:
function [dimension standard_dev] = fractalvol(data)
%[DIMENSION STANDARD_DEV] = fractalvol(DATA) calculates the fractal
%dimension of the 1 dimensional random walk, DATA. DATA is assumed to be a
%function of its indices.
%   Finds fractal volatility by embedding in the unit square and box
%   counting. Y axis will be rescaled values of DATA, x axis is
%   (1:length(DATA))/length(DATA).
%
%   Uncomment the else construction if you have the parallel computing
%   toolbox and wish to run the code in parallel (it scales linearly in instances).

if size(data,1) < size(data,2)
    data = data';
end

if size(data,2) == 1
    data = [(1:length(data))' data];
end

max1 = max(data(:,1));
min1 = min(data(:,1));
max2 = max(data(:,2));
min2 = min(data(:,2));

%normalize all to unit square
normdata = data;
normdata(:,1) = normdata(:,1)-min1;
normdata(:,1) = normdata(:,1)./(max1-min1);
normdata(:,2) = normdata(:,2)-min2;
normdata(:,2) = normdata(:,2)./(max2-min2);

%make sure nothing falls through the x axis
minwidth = min(diff(normdata(:,1)));
minwidth = log2(minwidth);
minwidth = abs(ceil(minwidth));
minwidth = minwidth-1;

n = zeros(minwidth,1);

% parallcntrl = matlabpool('size');

% if parallcntrl == 0
    for j = 1:minwidth
        width = 2^-j;
        xaxis_pos = 0;
        boxcount = 0;
        while xaxis_pos<1
            indx = (xaxis_pos <= normdata(:,1) &...
                normdata(:,1)<xaxis_pos+width);
            if 1-xaxis_pos == width
                indx(end) = true;
            end

            vertical_column = normdata(indx,2);

            if length(vertical_column) == 1
                boxcount = boxcount + 1;
            else
                rawcount = (max(vertical_column)-min(vertical_column))/width;
                rawcount = rawcount + rem(min(vertical_column),width);
                count = ceil(rawcount);
                boxcount = boxcount + count;
            end
            xaxis_pos = xaxis_pos + width; %advance on x axis
        end
        n(j) = boxcount;
    end
% else
%     parfor j = 1:minwidth
%         width = 2^-j;
%         xaxis_pos = 0;
%         boxcount = 0;
%         while xaxis_pos<1
%             indx = (xaxis_pos <= normdata(:,1) &...
%                 normdata(:,1)<xaxis_pos+width);
%             if 1-xaxis_pos == width
%                 indx(end) = true;
%             end
%
%             vertical_column = normdata(indx,2);
%
%             if length(vertical_column) == 1
%                 boxcount = boxcount + 1;
%             else
%                 rawcount = (max(vertical_column)-min(vertical_column))/width;
%                 rawcount = rawcount + rem(min(vertical_column),width);
%                 count = ceil(rawcount);
%                 boxcount = boxcount + count;
%             end
%             xaxis_pos = xaxis_pos + width; %advance on x axis
%         end
%         n(j) = boxcount;
%     end
% end

r = 2.^-(1:minwidth);
r = r';

s=-gradient(log(n))./gradient(log(r));
IQR = iqr(s);

indx2 = abs(s-median(s)) > IQR/2;

x2 = log(r);
y2 = log(n);

s(indx2)= [];
x2(indx2) = [];
y2(indx2) = [];

%std is the variance of the slope according to OLS theory

X = [ones(size(x2)) x2];
beta = pinv(X)*y2;
C = pinv(X'*X);
e=y2-X*pinv(X)*y2;
sigma = e'*e*C;
sigma = sqrt(sigma);

dimension = -beta(2);
standard_dev = sigma(2,2);