text:
function bst_startup(BrainstormHomeDir, GuiLevel, BrainstormDbDir)
% BST_STARTUP: Start a new Brainstorm Session.
%
% USAGE:  bst_startup(BrainstormHomeDir, GuiLevel=1, BrainstormDbDir=[])
%
% INPUTS:
%    - BrainstormHomeDir : Path to the brainstorm3 folder
%    - GuiLevel          : -1=server, 0=nogui, 1=normal, 2=autopilot
%    - BrainstormDbDir   : Database folder to use by default in this session

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

%% ===== MATLAB CHECK =====
% Parse inputs
if (nargin < 3) || isempty(BrainstormDbDir)
    BrainstormDbDir = [];
end
% If version is too old
MatlabVersion = bst_get('MatlabVersion');
if (MatlabVersion < 701)
    error('Brainstorm needs a version of Matlab >= 7.1');
end
% Is Matlab running (if not it is a compiled version)
isCompiled = bst_iscompiled();
% Compiled version: Force system look and feel
if isCompiled
    try
        javax.swing.UIManager.setLookAndFeel(javax.swing.UIManager.getSystemLookAndFeelClassName());
    catch
        % Whatever....
    end
end

% Startup message
disp(' ');
disp('BST> Starting Brainstorm:');
disp('BST> =================================');

%% ===== BRAINSTORM VERIFICATIONS =====
% Check that no interface is already running
if isappdata(0, 'BrainstormRunning')
    disp('BST> Brainstorm is already running. Restarting...');
    bst_exit();
end
% Initialize shared structure
global GlobalData;
GlobalData = db_template('GlobalData');
GlobalData.Program.GuiLevel = GuiLevel;
GlobalData.DataBase.LastSavedTime = tic();   % Save the current time, to know when to save the database
% Save the software home directory
bst_set('BrainstormHomeDir', BrainstormHomeDir);
% Debugging: show path in compiled application
if isCompiled
    disp(['BST> BrainstormHomeDir = ' BrainstormHomeDir]);
end
% Test for headless mode
if (GuiLevel >= 0) && (java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment.isHeadless() || java.awt.GraphicsEnvironment.getLocalGraphicsEnvironment.isHeadlessInstance())
    disp(' ');
    error(['Cannot create graphic interface.' 10 'If running Brainstorm in headless mode on a distant server, run "brainstorm server".']);
end
% Splash screen
if (GuiLevel == 1)
    bst_splash('show');
end

% === BRAINSTORM VERSION ===
try
    % Get doc folder
    docDir = bst_get('BrainstormDocDir');
    % Open "version.txt"
    fid = fopen(bst_fullfile(docDir, 'version.txt'), 'rt');
    % Get program name
    Name = fgetl(fid);
    Name = Name(3:end); %trim the comment
    % Get version and date
    strVer = fgetl(fid);
    cellVer = textscan(strVer(3:end), 'v. %s %s');
    Version = cellVer{1}{1};
    Release = cellVer{1}{1}(3:end);
    Date = cellVer{2}{1}(2:end-1);
    % Try to get GIT commit
    Commit = fgetl(fid);
    if ischar(Commit) && (length(Commit) >= 40)
        Commit = Commit(10:end);
    else
        Commit = [];
    end
    % Close file
    fclose(fid);
catch
    Name = 'Brainstorm';
    Version = '?';
    Release = '??????';
    Date    = '?';
    Commit  = [];
end
% If the commit is not available from the version.txt file, try to get it from the .git folder (if cloned from github)
if isempty(Commit)
    gitMaster = bst_fullfile(BrainstormHomeDir, '.git', 'refs', 'heads', 'master');
    if exist(gitMaster, 'file')
        fid = fopen(gitMaster, 'rt');
        if (fid >= 0)
            strGit = fgetl(fid);
            if ischar(strGit) && (length(strGit) >= 30)
                Commit = strGit;
            end
            fclose(fid);
        end
    end
end
% Save version in matlab preferences
bstVersion = struct('Name',    Name, ...
                    'Version', Version, ...
                    'Release', Release, ...
                    'Date',    Date, ...
                    'Commit',  Commit);
bst_set('Version', bstVersion);
% Display version number
disp(['BST> Version: ' Date ]);
% Get release date
localRel.year  = str2num(Release(1:2));
localRel.month = str2num(Release(3:4));
localRel.day   = str2num(Release(5:6));

% Check Matlab version
if (MatlabVersion <= 803)
    disp('BST> Warning: For better graphics, use Matlab >= 2014b');
end

%% ===== FORCE COMPILATION OF SOME INTERFACE FILES =====
if (GuiLevel == 1)
    disp('BST> Compiling main interface files...');
    tree_callbacks();
    bst_figures();
    figure_topo();
    figure_3d();
    figure_mri();
    figure_timeseries();
    figure_timefreq();
    bst_colormaps();
    bst_memory();
    bst_navigator();
end

%% ===== EMPTY REPORTS DIRECTORY =====
% Get reports directory
reportsDir = bst_get('UserReportsDir');
% If directory exists
if isdir(reportsDir)
    disp('BST> Deleting old process reports...');
    % List contents of folder
    listDir = dir(bst_fullfile(reportsDir, 'report_*.mat'));
    % If there are files in this folder (older versions of Matlab do not have this datenum field)
    if ~isempty(listDir) && isfield(listDir, 'datenum')
        % Get files that are older than 15 days
        iOldFiles = find(now() - [listDir.datenum] > 15);
        % Delete contents of directory
        for iFile = 1:length(iOldFiles)
            file_delete(bst_fullfile(reportsDir, listDir(iOldFiles(iFile)).name), 1);
        end
    end
end

%% ===== LOAD CONFIG FILE =====
disp('BST> Loading configuration file...');
% Get user database file : brainstorm.mat
dbFile = bst_get('BrainstormDbFile');
% Current DB version
CurrentDbVersion = 5.02;
% Get default colormaps list
sDefColormaps = bst_colormaps('Initialize');
isDbLoaded = 0;
% If file exists: load it
if file_exist(dbFile)
    % Load database file
    try
        bstOptions = load(dbFile);
    catch
        bst_splash('hide');
        java_dialog('msgbox', [...
            'Error: The database file was not saved properly.' 10 10 ...
            'Possible reason: Your hard drive is full or your quota exceeded.' 10 ...
            'Your user options are lost, but your database is probably safe:' 10 ...
            ' - Try to delete files in your home folder' 10 ...
            ' - Change the temporary folder in the Brainstorm preferences' 10 ...
            ' - Import again your database folder (File > Import database).'], 'Database error');
        bstOptions = [];
    end
    % Invalid structure read from dbFile
    if any(~isfield(bstOptions, {'iProtocol', 'ProtocolsListInfo', 'ProtocolsListSubjects', 'ProtocolsListStudies', 'BrainStormDbDir'}))
        disp(['BST> Warning: Ignoring corrupted options file: ' dbFile]);
        bstOptions = [];
    end
else
    bstOptions = [];
end
% Copy saved preferences to current instance
if ~isempty(bstOptions)
    % Add its contents in root app data
    if isfield(bstOptions, 'iProtocol')
        GlobalData.DataBase.iProtocol          = bstOptions.iProtocol;
        GlobalData.DataBase.ProtocolInfo       = bstOptions.ProtocolsListInfo;
        GlobalData.DataBase.ProtocolSubjects   = bstOptions.ProtocolsListSubjects;
        GlobalData.DataBase.ProtocolStudies    = bstOptions.ProtocolsListStudies;
        GlobalData.DataBase.BrainstormDbDir    = bstOptions.BrainStormDbDir;
        GlobalData.DataBase.isProtocolModified = zeros(1, length(bstOptions.ProtocolsListInfo));
        if isfield(bstOptions, 'DbVersion') && ~isempty(bstOptions.DbVersion)
            GlobalData.DataBase.DbVersion = bstOptions.DbVersion;
        end
        if isfield(bstOptions, 'isProtocolLoaded') && ~isempty(bstOptions.isProtocolLoaded)
            GlobalData.DataBase.isProtocolLoaded = bstOptions.isProtocolLoaded;
        else
            GlobalData.DataBase.isProtocolLoaded = ones(1, length(bstOptions.ProtocolsListInfo));
        end
        isDbLoaded = 1;
    end
    % Get saved colormaps
    if isfield(bstOptions, 'Colormaps') && isstruct(bstOptions.Colormaps)
        fNames = fieldnames(sDefColormaps);
        if (length(fieldnames(bstOptions.Colormaps)) ~= length(fNames)) || ~isequal(fieldnames(bstOptions.Colormaps.(fNames{1})), fieldnames(sDefColormaps.(fNames{1})))
            disp('BST> Colormap structure was updated. Fixing...');
        else
            GlobalData.Colormaps = bstOptions.Colormaps;
        end
    end
    % Clone control
    if isfield(bstOptions, 'CloneLock') && ~isempty(bstOptions.CloneLock)
        GlobalData.Program.CloneLock = bstOptions.CloneLock;
    end
    % Get saved preferences
    if isfield(bstOptions, 'Preferences') && isstruct(bstOptions.Preferences)
        GlobalData.Preferences = struct_copy_fields(GlobalData.Preferences, bstOptions.Preferences, 0);
    end
    % Get saved montages
    if isfield(bstOptions, 'ChannelMontages') && isstruct(bstOptions.ChannelMontages) && ...
            all(isfield(bstOptions.ChannelMontages, fieldnames(GlobalData.ChannelMontages))) && ...
            (length(bstOptions.ChannelMontages.Montages) > 20)
        GlobalData.ChannelMontages = bstOptions.ChannelMontages;
        % Reset butterfly plot selection
        panel_montage('SetCurrentMontage', 'MEG', []);
    end
    % Get saved process pipelines
    if isfield(bstOptions, 'Pipelines') && isstruct(bstOptions.Pipelines) && ...
            all(isfield(bstOptions.Pipelines, {'Name', 'Processes'}))
       GlobalData.Processes.Pipelines = bstOptions.Pipelines;
    end
    % Get saved searches
    if isfield(bstOptions, 'Searches') && isstruct(bstOptions.Searches) && ...
            all(isfield(bstOptions.Searches, {'Name', 'Search'}))
       GlobalData.DataBase.Searches.All = bstOptions.Searches;
    end
    % Reset current search filter
    if isfield(GlobalData.Preferences, 'NodelistOptions') && isfield(GlobalData.Preferences.NodelistOptions, 'String') && ~isempty(GlobalData.Preferences.NodelistOptions.String)
        GlobalData.Preferences.NodelistOptions.String = '';
    end
    % Reset previous exploration mode
    if isfield(GlobalData.Preferences, 'Layout') && isfield(GlobalData.Preferences.Layout, 'PreviousExplorationMode')
        GlobalData.Preferences.Layout.PreviousExplorationMode = GlobalData.Preferences.Layout.ExplorationMode;
    end
    % Check database structure for updates
    db_update(CurrentDbVersion);
end
if GlobalData.DataBase.DbVersion == 0
    % Database version is not defined, so it up-to-date
    GlobalData.DataBase.DbVersion = CurrentDbVersion;
end
% Check that Colormaps are defined
if isempty(GlobalData.Colormaps)
    GlobalData.Colormaps = sDefColormaps;
end
% Check that default montages are loaded
if (length(GlobalData.ChannelMontages.Montages) < 5) || any(~ismember({'CTF LF', 'Bad channels', 'Average reference (L -> R)', 'Scalp current density', 'Scalp current density (L -> R)', 'Head distance'}, {GlobalData.ChannelMontages.Montages.Name}))
    disp('BST> Loading default montages...');
    % Load default selections
    panel_montage('LoadDefaultMontages');
end

%% ===== AUTOMATIC UPDATES =====
% Automatic updates disabled: do not check for internet connection
if ~bst_get('AutoUpdates')
    disp('BST> Warning: Automatic updates are disabled.');
    disp('BST> Warning: Make sure your version of Brainstorm is up to date.');
% Matlab is running: check for updates
elseif ~isCompiled && (GuiLevel == 1)
    % Check internect connection
    fprintf(1, 'BST> Checking internet connectivity... ');
    [GlobalData.Program.isInternet, onlineRel] = bst_check_internet();
    % If no internet connection
    if ~GlobalData.Program.isInternet
        disp('failed');
        disp('BST> Could not check for Brainstorm updates.')
    else
        disp('ok');
        % Determine if release is old (local version > 30 days older than online version)
        daysOnline = onlineRel.year*365 + onlineRel.month*30 + onlineRel.day;
        daysLocal  =  localRel.year*365 +  localRel.month*30 +  localRel.day;
        isOld = ((daysOnline - daysLocal) > 30);
        % Display online version number
        if ((daysOnline - daysLocal) > 1)
            strOnline = datestr([2000+onlineRel.year, onlineRel.month, onlineRel.day, 0, 0, 0]);
            disp(['BST> Update available online: ' strOnline '']);
        end
        % Checking version: download if more then one month old
        if isOld
            disp('BST> Your version of brainstorm is old. Update is required.');
            % Check access rights
            if ~file_attrib(bst_fileparts(BrainstormHomeDir, 1), 'w') || ~file_attrib(BrainstormHomeDir, 'w')
                disp('BST> Brainstorm installation folder is read-only. Cannot update...');
            else
                % Hide splash screen
                bst_splash('hide');
                % Update brainstorm
                isUpdated = bst_update(1);
                % If update successful: Matlab or Brainstorm restart
                if isUpdated
                    return
                end
            end
        end
    end
% Check online connectivity
else
    [GlobalData.Program.isInternet, onlineRel] = bst_check_internet();
end

%% ===== START BRAINSTORM GUI =====
% Get screen configuration
GlobalData.Program.ScreenDef = gui_layout('GetScreenClientArea');
% Create main window (skipped in server mode)
if (GuiLevel >= 0)
    disp('BST> Initializing user interface...');
end
gui_initialize();
% Abort if something went wrong
if isempty(GlobalData.Program.GUI)
    return;
end

%% ===== INITIALIZE DATABASE =====
if ~isDbLoaded
    % Initialize structures
    GlobalData.DataBase.iProtocol          = 0;
    GlobalData.DataBase.ProtocolInfo       = repmat(db_template('ProtocolInfo'), 0);
    GlobalData.DataBase.ProtocolSubjects   = repmat(db_template('ProtocolSubjects'), 0);
    GlobalData.DataBase.ProtocolStudies    = repmat(db_template('ProtocolStudies'), 0);
    GlobalData.DataBase.isProtocolLoaded   = [];
    GlobalData.DataBase.isProtocolModified = [];
end

%% ===== CHECK FOR EEGLAB INSTALL =====
fminPath = lower(which('fminsearch'));
if ~isempty(strfind(fminPath, 'eeglab'))
    strProg = 'EEGLAB';
elseif ~isempty(strfind(fminPath, 'spm'))
    strProg = 'SPM';
elseif ~isempty(strfind(fminPath, 'fieldtrip'))
    strProg = 'FieldTrip';
else
    strProg = [];
end
if ~isempty(strProg)
    if (GuiLevel <= 0)
        disp(['BST> Warning: Some ' strProg ' functions shadow Matlab''s standard functions.']);
        disp(['BST> Warning: Please remove ' strProg ' from your Matlab path.']);
    else
        bst_splash('hide');
        java_dialog('warning', [strProg ' is installed on your system and shadows some standard Matlab functions.' 10 ...
                                'Without access to the function fminsearch, Brainstorm will not run properly.' 10 10 ...
                                'Please remove ' strProg ' from your Matlab path and restart Brainstorm.']);
    end
end

%% ===== START OPENGL =====
if (GuiLevel >= 0)
    fprintf(1, 'BST> Starting OpenGL engine... ');
    [isOpenGL, DisableOpenGL] = panel_options('StartOpenGL');
    % If OpenGL cannot be used: display a warning message
    if ~isOpenGL
        disp('BST> Warning: No OpenGL support available for this computer.');
        disp('BST>          Display will be slow and ugly.');
    % If OpenGL is manually disabled
    elseif isOpenGL && (DisableOpenGL == 1)
        disp('BST> Warning: ');
        disp('BST>    * Using this option causes the display to be slow and ugly.');
        disp('BST>    * Select only if you are experiencing serious display bugs with ');
        disp('BST>    * the full hardware acceleration. To edit this option: ');
        disp('BST>    * Menu: File > Set preferences... > Disable OpenGL rendering.');
    end
end

%% ===== LOAD PLUGINS =====
% Get installed plugins
[InstPlugs, AllPlugs] = bst_plugin('GetInstalled');
% Check installed plugins
if ~isempty(InstPlugs)
    % Display the plugins that are using custom installed path
    iPlugCustom = find([InstPlugs.isLoaded] & ~[InstPlugs.isManaged]);
    for iPlug = iPlugCustom
        disp(['BST> Plugin ' InstPlugs(iPlug).Name ': ' InstPlugs(iPlug).Path]);
        if strcmpi(InstPlugs(iPlug).Name, 'spm12') && isempty(strfind(spm('ver'), 'SPM12'))
            disp(['BST> ** WARNING: Installed version is not SPM12: ' spm('ver') ' **']);
        end
    end
    % Load plugins that should be loaded automatically at startup
    iPlugLoad = find([InstPlugs.AutoLoad] & ~[InstPlugs.isLoaded]);
    if ~isempty(iPlugLoad)
        fprintf('BST> Loading plugins... ');
    end
    for iPlug = iPlugLoad
        bst_plugin('Load', InstPlugs(iPlug), 0);
        fprintf([InstPlugs(iPlug).Name, ' ']);
    end
    if ~isempty(iPlugLoad)
        fprintf('\n');
    end
end

%% ===== PARSE PROCESS FOLDER =====
% Parse process folder
disp('BST> Reading process folder...');
panel_process_select('ParseProcessFolder', 1);

%% ===== INSTALL ANATOMY TEMPLATE =====
% Download ICBM152 template if missing (e.g. when cloning from GitHub)
TemplateDir = fullfile(BrainstormHomeDir, 'defaults', 'anatomy', 'ICBM152');
if ~isCompiled && ~exist(TemplateDir, 'file')
    TemplateName = 'ICBM152_2023b';
    isSkipTemplate = 0;
    % Template file
    ZipFile = bst_fullfile(bst_get('UserDefaultsDir'), 'anatomy', [TemplateName '.zip']);
    % If template is not downloaded yet: download it
    if ~exist(ZipFile, 'file')
        disp('BST> Downloading ICBM152 template...');
        % Download file
        errMsg = gui_brainstorm('DownloadFile', ['http://neuroimage.usc.edu/bst/getupdate.php?t=' TemplateName], ZipFile, 'Download template');
        % Error message
        if ~isempty(errMsg)
            disp(['BST> Error: Could not download template: ' errMsg]);
            isSkipTemplate = 1;
        end
    end
    % If the template is available as a zip file
    if ~isSkipTemplate
        disp('BST> Installing ICBM152 template...');
        % Create folder
        mkdir(TemplateDir);
        % URL: Download zip file
        try
            unzip(ZipFile, TemplateDir);
        catch
            disp(['BST> Error: Could not unzip anatomy template: ' lasterr]);
        end
    end
end

%% ===== LICENSE AGREEMENT =====
if (GuiLevel == 1)
    % Number of days to allow as grace period for renewing license
    GRACE = 15;
    % Get previous agreement date (default: current date)
    if isfield(GlobalData, 'Preferences') && isfield(GlobalData.Preferences, 'DateofAgreement') && ~isempty(GlobalData.Preferences.DateofAgreement)
        DateofAgreement = GlobalData.Preferences.DateofAgreement;
    else
        DateofAgreement = datestr(floor(now) - GRACE - 1);
    end
    % Get number of days since last agreement
    DaysSinceAgree = etime(datevec(now),datevec(DateofAgreement));
    DaysSinceAgree = DaysSinceAgree/(60*60*24); % convert seconds to days

    % If user did not agree to Brainstorm license rencently
    if (DaysSinceAgree >= GRACE)
        % Hide splash screen
        bst_splash('hide');
        % Show license agreement panel
        isOk = bst_license();
        % If user did not agree: exit
        if ~isOk
            clear all
            disp('BST> License agreement unsatisfied. Closing Brainstorm...');
            disp('BST> Type ''brainstorm'' to restart Brainstorm.');
            % Release Brainstorm global mutex
            bst_mutex('release', 'Brainstorm');
            return;
        % Else accept validation for 15 days
        else
            disp('BST> License accepted.');
            GlobalData.Preferences.DateofAgreement = datestr(floor(now));
        end
    end
end

%% ===== TEST MEMORY =====
if ispc && (MatlabVersion >= 706) && (GuiLevel == 1)
    try
        % Get Matlab memory
        usermem = memory();
        % Minimum: 1Gb RAM contiguous
        maxsize = round(usermem.MaxPossibleArrayBytes / 1024 / 1024);
        if (maxsize < 1024) && ~bst_get('IgnoreMemoryWarnings')
            disp(sprintf('BST> Warning: Maximum variable size: %d Mb', maxsize));
            % Hide splash screen
            bst_splash('hide');
            % Display warning
            java_dialog('msgbox', ...
                ['Your system reports to be short on memory.' 10 10 ...
                 'The maximum block of memory that Matlab can allocate is: ' num2str(maxsize) ' Mb.' 10 ...
                 'The recommended minimum for running Brainstorm is 1024 Mb.' 10 ...
                 'Below this limit, you may experiment a lot of "out of memory" errors.' 10 10 ...
                 'Close all the other applications running on this computer to free some memory.' 10 ...
                 'To ignore this message permanently: File > Preferences > Ignore memory warnings.'], ...
                'Memory warning');
        end
    catch
        % Whatever...
    end
end

%% ===== SET DATABASE DIRECTORY =====
isImportDb = 0;
% Get user dir
BrainstormUserDir = bst_get('BrainstormUserDir');
% Get database folder
if isempty(BrainstormDbDir)
    BrainstormDbDir = bst_get('BrainstormDbDir');
% If database folder was passed in input
else
    if isequal(BrainstormDbDir, 'local')
        BrainstormDbDir = bst_fullfile(BrainstormUserDir, 'local_db');
    end
    % Save brainstorm directory
    bst_set('BrainstormDbDir', BrainstormDbDir);
    disp(['BST> Database folder: ' BrainstormDbDir]);
end
% If folder is not defined yet: ask user to set it
if isempty(BrainstormDbDir)
    % Hide splash screen
    bst_splash('hide');
    % Display message: first startup
    java_dialog('msgbox', ['It is the first time you run Brainstorm on this Matlab installation.' 10 10 ...
                           'First of all, you need to create a new directory to store the Brainstorm database,' 10 ...
                           'called for instance "brainstorm_db".' 10 10 ...
                           'IMPORTANT NOTES: ' 10 ...
                           '- Do not create this database directory in the Brainstorm program directory' 10 ...
                           '- The database directory must contain only files created by Brainstorm.' 10 ...
                           '- Do not put your original data files and personal results in the database directory.' 10 ...
                           '- Do not put any file in the Brainstorm program directory.'], 'Brainstorm setup');
    % Set database folder
    BrainstormDbDir = gui_brainstorm('SetDatabaseFolder');
    % If no directory selected : exit
    if isempty(BrainstormDbDir)
        % Release Brainstorm global mutex
        bst_mutex('release', 'Brainstorm');
        return
    end
    % Check if there are protocols in this folder
    if ~isempty(file_find(BrainstormDbDir, 'brainstormsubject*.mat', 4))
        % Ask if user wants to import all the database
        isImportDb = java_dialog('confirm', ['This folder already contains Brainstorm protocols.' 10 10 ...
                                             'Load all these protocols now ?' 10 10]);
    end
end

%% ===== INITIALIZATION DONE =====
disp('BST> Loading current protocol...');
% Get handle to the main window
jFrame = bst_get('BstFrame');
% If the GUI is requested by the user
if (GuiLevel == 1)
    % Show main Brainstorm window
    jFrame.setVisible(1);
    % Weird thing with some macs: by default, window is "always on top"
    if strncmp(computer,'MAC',3)
        jFrame.setVisible(0);
        jFrame.setAlwaysOnTop(0);
        jFrame.setVisible(1);
    end
end
% Headless mode: load protocol data
if (GuiLevel == -1)
    if ~bst_get('isProtocolLoaded')
        db_load_protocol(GlobalData.DataBase.iProtocol);
    end
% Regular GUI: Prepare protocol explorer panel
else
    % Display Brainstorm version
    jFrame.setTitle(['Brainstorm ' Date]);
    % Read the protocols list in UserDataBase
    gui_brainstorm('UpdateProtocolsList');
    % Load the selected protocol
    gui_brainstorm('SetCurrentProtocol', GlobalData.DataBase.iProtocol);
    % Update permanent panels (to disable them)
    panel_surface('UpdatePanel');
    panel_scout('UpdatePanel');
    panel_cluster('UpdatePanel');
end
% Get decoration size
GlobalData.Program.DecorationSize = gui_layout('GetDecorationSize', jFrame);
disp('BST> =================================');
disp(' ');
% Set a flag to mark that brainstorm is now running
setappdata(0, 'BrainstormRunning', 1);
% Hide spash screen
bst_splash('hide');
% Make sure that figure named "Brainstorm" is hidden
hMutex = bst_mutex('get', 'Brainstorm');
set(hMutex, 'Visible', 'off');

%% ===== DELETE OLD PLUGINS =====
% Start with GUI
if (GuiLevel == 1)
    % Add bst_duneuro (now called 'duneuro')
    AllPlugs(end+1).Name = 'bst_duneuro';
    AllPlugs(end+1).Name = 'nirstorm';
    iOldInstall = find(cellfun(@(c)exist(fullfile(BrainstormUserDir,c),'file'), {AllPlugs.Name}));
    % Some old plugins were detected: ask user
    if ~isempty(iOldInstall)
        OldPlugPath = cellfun(@(c)fullfile(BrainstormUserDir,c), {AllPlugs(iOldInstall).Name}, 'UniformOutput', 0);
        isConfirm = java_dialog('confirm', ['The plugin system was updated.' 10 10 ...
            'All the Brainstorm plugins are now managed from the menu "Plugins".' 10 ...
            'The old plugins cannot be used anymore and need to be deleted.' 10 ...
            'Next time you need them, they will be downloaded again automatically.' 10 10 ...
            'Delete the old plugins listed below?' 10 ...
            sprintf(' - %s\n', OldPlugPath{:})], 'Plugin manager');
        % Delete plugin
        if isConfirm
            % Delete files in $HOME/.brainstorm/plugname
            file_delete(OldPlugPath, 1, 3);
            % NIRSTORM: Call uninstall function (it previously installed functions in various $HOME/.brainstorm subfolders)
            if any(strcmpi({AllPlugs(iOldInstall).Name},'nirstorm')) && exist('uninstall_nirstorm')
                cur_dir=pwd;
                cd(bst_get('UserProcessDir'));
                uninstall_nirstorm();
                file_delete(fullfile(bst_get('UserProcessDir'),'uninstall_nirstorm.m'),1);
                file_delete(fullfile(bst_get('UserProcessDir'),{'dg_voronoi.mexa64','dg_voronoi.mexglx'}),1);
                cd(cur_dir);
            end
        end
    end
end

%% ===== RELOAD ALL DATABASE =====
% Create file to indicate that Brainstorm was started
StartFile = bst_fullfile(BrainstormUserDir, 'is_started.txt');
% If import database was mandatory
if isImportDb
    disp('BST> Reloading database...');
    db_import(BrainstormDbDir);
% Check if the program was closed unexpectedly
elseif file_exist(StartFile) && (GuiLevel == 1)
    % Delete this file
    if (file_delete(StartFile, 1) == 1) && ~file_exist(StartFile)
        isReloadCurrent = java_dialog('confirm', ...
            ['Brainstorm was not closed properly, some modifications to the current' 10 ...
             'protocol might not have been saved properly in the dabase.' 10 10 ...
             'Reload the current protocol ?' 10 10]);
        if isReloadCurrent
            db_reload_database(GlobalData.DataBase.iProtocol);
        end
    else
        disp('BST> Warning: Brainstorm is already running from a different Matlab session.');
    end
end
% Open a new file to track if Brainstorm is opened
fid = fopen(StartFile, 'w');
fwrite(fid, ['Brainstorm started: ' datestr(now) 10]);
fwrite(fid, ['User: ' char(java.lang.System.getProperty('user.name')) 10]);
try
    jLocalHost = java.net.InetAddress.getLocalHost();
    fwrite(fid, ['Host: ' char(jLocalHost.getHostName()) ' / ' char(jLocalHost.getHostAddress()) ' (' char(java.lang.System.getProperty('os.name')) ' ' char(java.lang.System.getProperty('os.version')) ')' 10]);
catch
    fwrite(fid, ['Host: Unknown (' char(java.lang.System.getProperty('os.name')) ' ' char(java.lang.System.getProperty('os.version')) ')' 10]);
end
fclose(fid);

%% ===== SET TEMPORARY FOLDER =====
% Get temp folder
TmpDir = bst_get('BrainstormTmpDir');
% If folder is not defined yet: ask user to set it
if ~isempty(strfind(TmpDir, '/home/bic/'))
    % Hide splash screen
    bst_splash('hide');
    % Display message: first startup
    java_dialog('msgbox', [...
        'Warning: You should change the temporary directory.' 10 10 ...
        'The temporary folder used by Brainstorm is set by default to:' 10 '    ' TmpDir 10 ...
        'You only have 1Gb of available space in the folder /home/bic/, ' 10 ...
        'this is might not be enough for Brainstorm to work properly.' 10 10 ...
        'Please create a folder "brainstorm_tmp" on your local hard drive:' 10 ...
        '    /export01/data/username/brainstorm_tmp or ' 10 ...
        '    /export02/username/brainstorm_tmp', 10 10 ...
        'Then change the temporary directory in the next window.' 10 ], ...
       'BIC workstation: Incorrect temporary folder');
    % Edit preferences
    gui_show('panel_options', 'JavaWindow', 'Brainstorm preferences', [], 1, 0, 0);
end

% Empty temporary folder with confirmation (if nogui/server: display warning)
gui_brainstorm('EmptyTempFolder', 1);

%% ===== PREPARE BUG REPORTING =====
% % Get current configuration
% BugReportOptions = bst_get('BugReportOptions');
% % If incomplete: ask user to complete it
% if BugReportOptions.isEnabled && (isempty(BugReportOptions.SmtpServer) || isempty(BugReportOptions.UserEmail))
%     gui_show_dialog('Bug reporting', @panel_bug);
% end

%% ===== COMPILED MODE: WAIT =====
if isCompiled
    %     % Loop to wait for the end
    %     while brainstorm('status')
    %         pause(2);
    %     end
    %     % Exit Matlab
    %     quit force;

    % Wait until the main Brainstorm mutex gets deleted
    waitfor(bst_mutex('get', 'Brainstorm'));
end

author: brainstorm-tools


text:
function res = CompMat(a);
%COMPMAT      Ostrowski's comparison matrix
%
%   res = CompMat(a)
%

% written   8/12/94     S.M. Rump
% modified 09/22/02     S.M. Rump  check square matrix and interval matrices
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  [m n] = size(a);
  if m~=n
    error('Comparison matrix only for square matrices')
  end

  res = -mag(a);
  res(1:n+1:n*n) = mig(diag(a));

repo: OpenMAS


text:
function varargout = GUIEditTextualFileHeader(varargin)
% GUIEditTextualFileHeader : GUI for editing the SGY header
%
% Call :
%   SH=GUIEditTextualFileHeader(SH);
%
%
% Example :
%   [Data,STH,SH]=ReadSegy('841_m.sgy');
%   SH=GUIEditTextualFileHeader(SH);
%
%

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

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

% Last Modified by GUIDE v2.5 05-Jan-2009 17:27:47

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GUIEditTextualFileHeader_OpeningFcn, ...
                   'gui_OutputFcn',  @GUIEditTextualFileHeader_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);

if nargin & isstr(varargin{1})
  gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

% End initialization code - DO NOT EDIT

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

% Choose default command line output for GUIEditTextualFileHeader
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% Load TextualFileHeader if given
if length(varargin)>0
  if nargin & isstruct(varargin{1})
    if isfield(varargin{1},'TextualFileHeader')
    data=guidata(gcf);
    data.SegyHeader=varargin{1};
    data.SegyHeaderOrg=data.SegyHeader;

    if isempty(find(data.SegyHeader.TextualFileHeader>200))
        HeaderType=2; % ASCII
    else
        HeaderType=1; % EBCDIC
    end
    set(handles.popType,'Value',HeaderType);

    guidata(gcf, data);
    UpdateText(hObject, handles)
    end
  end
end

% UIWAIT makes GUIEditTextualFileHeader wait for user response (see UIRESUME)
uiwait(handles.figTextualHeader);

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

% Get default command line output from handles structure
data=guidata(gcf);
varargout{1} =     data.SegyHeader;
close(handles.figTextualHeader);

% --- Executes during object creation, after setting all properties.
function eTextHeader_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eTextHeader (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

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

% Hints: get(hObject,'String') returns contents of eTextHeader as text
%        str2double(get(hObject,'String')) returns contents of
%        eTextHeader as a double
%eTextHeader_Callback(handles.eTextHeader, eventdata, handles)
  HS=get(hObject,'String');
  VAL=get(hObject,'Value');

  data=guidata(hObject);
  data.HS=HS;
  data.VAL=VAL;
  %disp(sprintf('LINE=''%s''',HS(VAL,:)))
  guidata(hObject, data);

  % UPDATE TEXT STRING
  set(handles.eEditLine,'String',HS(VAL,:))

% --- Executes during object creation, after setting all properties.
function popType_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popType (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in popType.
function popType_Callback(hObject, eventdata, handles)
% hObject    handle to popType (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popType contents as cell array
%        contents{get(hObject,'Value')} returns selected item from
%        popType
  UpdateText(hObject,handles)

function UpdateText(hObject, handles)
  data=guidata(hObject);

  HeaderType=get(handles.popType,'Value');

  if HeaderType==1,
    txt=char(ebcdic2ascii(data.SegyHeader.TextualFileHeader));
  else
    txt=char(data.SegyHeader.TextualFileHeader');
  end
  %disp('--')
  %reshape(txt,80,40)'
  %disp('--')
  %disp(size(txt))
  try
    s=[];
    for i=1:40;
      s=[s,txt((i-1)*80+1:(i-1)*80+80)];
      try
          if i~=40, s=[s,'|']; end
      catch
          keyboard
      end
    end
  catch
    if (HeaderType==1)
      warndlg(['The Textual Header is most probably NOT EBCDIC formatted'],'Textual Header Info')
    else
      warndlg(['The Textual Header is most probably NOT ASCII formatted'],'Textual Header Info')
    end
  end

  try
    set(handles.eTextHeader,'String',s)
  catch
    disp('WHUUPS')
  end

% --- Executes during object creation, after setting all properties.
function eEditLine_CreateFcn(hObject, eventdata, handles)
% hObject    handle to eEditLine (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

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

% Hints: get(hObject,'String') returns contents of eEditLine as text
%        str2double(get(hObject,'String')) returns contents of eEditLine as a double

data=guidata(hObject);
if ~isfield(data,'HS')
    eTextHeader_Callback(handles.eTextHeader, eventdata, handles)
    data=guidata(hObject);
end

try

    NEWSTR=get(hObject,'String');
    n=length(NEWSTR);
    if n<80
        NEWSTR(n+1:80)=' ';
    end

    n=length(NEWSTR);
    if n>=80
        data.HS(data.VAL,:)=NEWSTR(1:80);
    else
        data.HS(data.VAL,1:n)=NEWSTR(1:n);
    end

    HStransp=data.HS';

    if get(handles.popType,'Value')==1;
        % EBCDIC
        data.SegyHeader.TextualFileHeader=double(ascii2ebcdic(HStransp(:)));
    else
        % ASCII
        data.SegyHeader.TextualFileHeader=double(HStransp(:));
    end
    guidata(gcf, data);

    UpdateText(hObject, handles)

    eTextHeader_Callback(handles.eTextHeader, eventdata, handles)
catch
    keyboard

end

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

% --- Executes on button press in pbReset.
function pbReset_Callback(hObject, eventdata, handles)
% hObject    handle to pbReset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
data=guidata(gcf);
data.SegyHeader=data.SegyHeaderOrg;
guidata(gcf, data);
UpdateText(hObject, handles)

% --- Executes on button press in pbFont.
function pbFont_Callback(hObject, eventdata, handles)
% hObject    handle to pbFont (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
d = uisetfont(handles.eTextHeader);
% Apply those settings to c2
set(handles.eEditLine,d);

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function value = p29_autonomous ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


function outPutFeatures = yolov3v4Predict(cfg_file,weight_file,image)
% 功能：yolov3快速输出检测特征，同darknet官网输出结果方式保持一致
% 输入：
%      cfg_file, 指定的cfg后缀的模型描述文件
%      weight_file, 指定的.weights后缀的二进制文件
%      image ：输入网络的图像数据,单张图像(H*W*C)或者批量图像(H*W*C*bs)
% 输出：
%     outPutFeatures： M*(5+nc)或者bs*M*(5+nc)大小
%     ,为bs*M*[x,y,w,h,Pobj,p1,p2,...,pn]大小的形式矩阵，如果是单张图像检测，则输出大小为M*(5+nc)，否则是bs*M*(5+nc),
%     其中，M为检测框的数量，bs为图片数量，nc为训练网络dlnet类别数量，x,y,w,h分别是输入图片上对应的x,y,width,height，Pobj
%     为物体概率，p1,p2,...,pn分别为对应coco.names类别概率
%
% author: cuixingxing
% email:cuixingxing150@gmail.com
% 2020.4.22创建
% 2020.5.2 修改
% 2020.5.13 minor fix
%

persistent dlnet yolovLayerArray netInputSize
if isempty(dlnet)

    %% import network and predict
    [layerGraphYolo,hyperParams] = importDarknetWeights(cfg_file,weight_file);
    dlnet = dlnetwork(layerGraphYolo); % 65秒左右
    % analyzeNetwork(layerGraphYolo)% visualize network
    % exportDarkNetNetwork(dlnet,hyperParams,'temp.cfg','temp.weights')

     %% get yolo index
    yoloIdx = [];
    for i = 1:length(dlnet.Layers)
        currentLayerType = class(dlnet.Layers(i));
        if strcmpi(currentLayerType,'yoloV3Layer')
            yoloIdx = [yoloIdx;i];
        end
    end
    assert(~isempty(yoloIdx),'输入网络非yolov3/4网络！')

    yolovLayerArray = dlnet.Layers(yoloIdx);
    netInputSize = dlnet.Layers(1).InputSize(1:2);
end

inputSize = [size(image,1),size(image,2)];
scale = inputSize./netInputSize;% [heightScale,widthScale]

img = imresize(im2single(image),netInputSize);% [0,1]数据，保持与训练时候同大小和类型
dlX = dlarray(img,'SSCB');% 大小为h*w*c*bs,注意是否归一化要看与训练时候图像一致
if(canUseGPU())
    dlX = gpuArray(dlX);% 推送到GPU上
end

numsYOLO = length(yolovLayerArray);
outFeatureMaps = cell(numsYOLO,1);
[outFeatureMaps{:}] = predict(dlnet,dlX,'Outputs',dlnet.OutputNames);% h*w*c*bs,matlab输出方式排列
outPutFeatures = [];
for i = 1:numsYOLO
    currentYOLOLayer = yolovLayerArray(i);
    currentFeatureMap = outFeatureMaps{i};

    % 由于yolov3Layer类里面predict函数未改变类属性，故重新给属性赋值
    currentYOLOLayer.numY = size(currentFeatureMap,1);
    currentYOLOLayer.numX = size(currentFeatureMap,2);
    currentYOLOLayer.stride = max(currentYOLOLayer.imageSize)./max(currentYOLOLayer.numX,...
        currentYOLOLayer.numY);

    % reshape currentFeatureMap到有意义的维度，h*w*c*bs --> h*w*(5+nc)*na*bs
    % --> bs*na*h*w*(5+nc),最终的维度方式与darknet官网兼容
    bs = size(currentFeatureMap,4);
    h = currentYOLOLayer.numY;
    w = currentYOLOLayer.numX;
    na = currentYOLOLayer.nAnchors;
    nc = currentYOLOLayer.classes;
    currentFeatureMap = reshape(currentFeatureMap,h,w,5+nc,na,bs);% h*w*(5+nc)*na*bs
    currentFeatureMap = permute(currentFeatureMap,[5,4,1,2,3]);% bs*na*h*w*(5+nc)

    [~,~,yv,xv] = ndgrid(1:bs,1:na,0:h-1,0:w-1);% yv,xv大小都为bs*na*h*w，注意顺序，后面做加法维度标签要对应
    gridXY = cat(5,xv,yv);% 第5维上扩展，大小为bs*na*h*w*2, x,y从1开始的索引
    currentFeatureMap(:,:,:,:,1:2) = sigmoid(currentFeatureMap(:,:,:,:,1:2)) + gridXY; % 大小为bs*na*h*w*2,预测对应xy
    anchor_grid = currentYOLOLayer.anchorsUse/currentYOLOLayer.stride; % 此处anchor_grid大小为na*2
    anchor_grid = reshape(anchor_grid,1,currentYOLOLayer.nAnchors,1,1,2);% 此处anchor_grid大小为1*na*1*1*2，方便下面相乘
    currentFeatureMap(:,:,:,:,3:4) = exp(currentFeatureMap(:,:,:,:,3:4)).*anchor_grid;% 大小为bs*na*h*w*2
    currentFeatureMap(:,:,:,:,1:4) = currentFeatureMap(:,:,:,:,1:4)*currentYOLOLayer.stride; % 预测的bboxes
    currentFeatureMap(:,:,:,:,5:end) = sigmoid(currentFeatureMap(:,:,:,:,5:end)); % 预测的scores

    if currentYOLOLayer.classes == 1
        currentFeatureMap(:,:,:,:,6) = 1;
    end
    currentFeatureMap = reshape(currentFeatureMap,bs,[],5+currentYOLOLayer.classes);% bs*N*(5+nc)

    if isempty(outPutFeatures)
        outPutFeatures = currentFeatureMap;
    else
        outPutFeatures = cat(2,outPutFeatures,currentFeatureMap);% bs*M*(5+nc)
    end
end

%% 坐标转换到原始图像上
outPutFeatures = extractdata(outPutFeatures);% bs*M*(5+nc) ,为[x_center,y_center,w,h,Pobj,p1,p2,...,pn]
outPutFeatures(:,:,[1,3]) = outPutFeatures(:,:,[1,3])*scale(2);% x_center,width
outPutFeatures(:,:,[2,4]) = outPutFeatures(:,:,[2,4])*scale(1);% y_center,height
outPutFeatures(:,:,1) = outPutFeatures(:,:,1) -outPutFeatures(:,:,3)/2;%  x
outPutFeatures(:,:,2) = outPutFeatures(:,:,2) -outPutFeatures(:,:,4)/2; % y
outPutFeatures = squeeze(outPutFeatures); % 如果是单张图像检测，则输出大小为M*(5+nc)，否则是bs*M*(5+nc)
if(canUseGPU())
    outPutFeatures = gather(outPutFeatures); % 推送到CPU上
end
end


text:
function [ ns, xyz ] = sphere_cubed_grid_points_face ( n, i1, j1, k1, i2, ...
  j2, k2, ns, xyz )

%*****************************************************************************80
%
%% SPHERE_CUBED_GRID_POINTS_FACE: points on a face of a cubed sphere grid.
%
%  Discussion:
%
%    This function generates points on a face of a cubed sphere grid, and
%    appneds them to a growing list.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sections into which each face of
%    the cube is to be divided.
%
%    Input, integer I1, J1, K1, I2, J2, K2, the logical indices, between 0
%    and N, of two corners of the face grid.  It is guaranteed that I1 <= I2,
%    J1 <= J2, and K1 <= K2.
%
%    Input, integer NS, the initial number of points.
%
%    Input, real XYZ(NS,3), distinct points on the unit sphere
%    generated by a cubed sphere grid.
%
%    Output, integer NS, the final number of points.
%
%    Output, real XYZ(NS,3), distinct points on the unit sphere
%    generated by a cubed sphere grid, including the points just genereated.
%
  for i = i1 : i2

    if ( i1 < i2 )
      xc = tan ( ( 2 * i - n ) * 0.25 * pi / n );
    elseif ( i1 == 0 )
      xc = -1.0;
    elseif ( i1 == n )
      xc = +1.0;
    else
      xc = 0.0;
    end

    for j = j1 : j2

      if ( j1 < j2 )
        yc = tan ( ( 2 * j - n ) * 0.25 * pi / n );
      elseif ( j1 == 0 )
        yc = -1.0;
      elseif ( j1 == n )
        yc = +1.0;
      else
        yc = 0.0;
      end

      for k = k1 : k2

        if ( k1 < k2 )
          zc = tan ( ( 2 * k - n ) * 0.25 * pi / n );
        elseif ( k1 == 0 )
          zc = -1.0;
        elseif ( k1 == n )
          zc = +1.0;
        else
          zc = 0.0;
        end

        xyzn = sqrt ( xc^2 + yc^2 + zc^2 );

        ns = ns + 1;
        xyz(ns,1) = xc / xyzn;
        xyz(ns,2) = yc / xyzn;
        xyz(ns,3) = zc / xyzn;

      end
    end
  end

  return
end


text:
%
%   FILE NAME:    integrateRxnAssocViaIDsAndMets.m
%
%   PURPOSE: Integrate rxn associations from mets and external rxn ids
%

load('MNXRxns.mat');          % load MNX reactions
load('mergedModel.mat');      % load mergedModel
load('mapRxnResults.mat');    % load automatic rxn association via mets

% Treat the MNX assoc overlapped by rxn-ids and via mets as confirmed
mergedModel.confirmedMNXID=cell(numel(mergedModel.rxns),1);
mergedModel.confirmedMNXID(:)={''};
mergedModel.confirmedFilteredMNXID=cell(numel(mergedModel.rxns),1);
mergedModel.confirmedFilteredMNXID(:)={''};
mergedModel.unResolved=zeros(numel(mergedModel.rxns),1);

numUnresloved=0;
numConfirmed=0;
numToCheck=0;
for i=1:length(mergedModel.rxns)
		if ~isempty(mergedModel.rxnAssocMNXID{i}) || ~isempty(results.rxnMNXID{i})
				overlap=intersect(mergedModel.rxnAssocMNXID{i},results.rxnMNXID{i});

				if ~isempty(overlap)
						numConfirmed=numConfirmed+1;
						mergedModel.confirmedMNXID{i}=overlap;
						mergedModel.confirmedFilteredMNXID{i}=filterBalancedRxns(overlap,MNXRxns);
				else
						numToCheck=numToCheck+1;
				end

		else  % These rxns have no association both by rxn-ids and mets
				numUnresloved=numUnresloved+1;
				mergedModel.unResolved(i)=1;
		end
end
%numConfirmed=2506;
%numUnresloved=726;
%numToCheck=739;

% Treat the cases need to be checked
filledRxnAssocViaMets.rxnHMRID=cell(1000,1);
filledRxnAssocViaMets.rxnMNXID=cell(1000,1);
filledRxnAssocViaMets.rxnFilteredMNXID=cell(1000,1);
filledRxnAssocViaMets.notes=cell(1000,1);
needToCheck.rxnHMRID=cell(1000,1);
needToCheck.rxnAssocMNXID=cell(1000,1);
needToCheck.rxnViaMetsMNXID=cell(1000,1);
needToCheck.rxnViaMetsFilteredMNXID=cell(1000,1);
needToCheck.notes=cell(1000,1);

numToCheck=0;
numFilledViaMets=0;
for i=1:numel(mergedModel.rxns)
		if ~isempty(mergedModel.rxnAssocMNXID{i}) || ~isempty(results.rxnMNXID{i})
				overlap=intersect(mergedModel.rxnAssocMNXID{i},results.rxnMNXID{i});
				if isempty(overlap)
						if ~isempty(mergedModel.rxnAssocMNXID{i})
								numToCheck=numToCheck+1;
								needToCheck.rxnHMRID{numToCheck}=mergedModel.rxns{i};
								needToCheck.rxnAssocMNXID{numToCheck}=mergedModel.rxnAssocMNXID{i};
								needToCheck.notes{numToCheck}=results.notes{i};

								% Both have assoc, but non-overlap (202) could be caused by wrong assoc by rxn-ids
								if ~isempty(results.rxnMNXID{i})
										needToCheck.rxnViaMetsMNXID{numToCheck}=results.rxnMNXID{i};
										needToCheck.rxnViaMetsFilteredMNXID{numToCheck}=filterBalancedRxns(results.rxnMNXID{i},MNXRxns);

								% No assoc via mets, keep origianl ones (248) Check the manual curation
								else
										needToCheck.rxnViaMetsMNXID{numToCheck}='';
										needToCheck.rxnViaMetsFilteredMNXID{numToCheck}='';
								end

						% No rxn assoc, fill by assoc via mets (prefer to true type) (296)
						else
								numFilledViaMets=numFilledViaMets+1;
								filledRxnAssocViaMets.rxnHMRID{numFilledViaMets}=mergedModel.rxns{i};
								filledRxnAssocViaMets.rxnMNXID{numFilledViaMets}=results.rxnMNXID{i};
								filledRxnAssocViaMets.rxnFilteredMNXID{numFilledViaMets}=filterBalancedRxns(results.rxnMNXID{i},MNXRxns);
								filledRxnAssocViaMets.notes{numFilledViaMets}=results.notes{i};
						end
				end
		end
end
%numToCheck=444;
%numFilledViaMets=295;

needToCheck.rxnHMRID=needToCheck.rxnHMRID(1:numToCheck);
needToCheck.rxnAssocMNXID=needToCheck.rxnAssocMNXID(1:numToCheck);
needToCheck.rxnViaMetsMNXID=needToCheck.rxnViaMetsMNXID(1:numToCheck);
needToCheck.rxnViaMetsFilteredMNXID=needToCheck.rxnViaMetsFilteredMNXID(1:numToCheck);
needToCheck.notes=needToCheck.notes(1:numToCheck);

filledRxnAssocViaMets.rxnHMRID=filledRxnAssocViaMets.rxnHMRID(1:numFilledViaMets);
filledRxnAssocViaMets.rxnMNXID=filledRxnAssocViaMets.rxnMNXID(1:numFilledViaMets);
filledRxnAssocViaMets.rxnFilteredMNXID=filledRxnAssocViaMets.rxnFilteredMNXID(1:numFilledViaMets);
filledRxnAssocViaMets.notes=filledRxnAssocViaMets.notes(1:numFilledViaMets);

%save('mergedModel.mat','mergedModel');
%save('filledRxnAssocViaMets_20180522.mat','filledRxnAssocViaMets');
%save('needToCheck_20180522.mat','needToCheck');

%% sub functions
%Filter MNX rxns according to their balance status
function filteredMNXrxns = filterBalancedRxns(MNXrxns,mnx)

		countNum=numel(MNXrxns);
		if countNum==1  % there is only one MNX id, skip balance check
				filteredMNXrxns=MNXrxns;
		elseif countNum>1   % there are several MNX ids
				[a, b]=ismember(MNXrxns,mnx.MNX_ID);
				% The priority order of Balance status: true > NA > ambiguous
				if find(strcmp('true',mnx.Balance(b)));
						trueHits=find(strcmp('true',mnx.Balance(b)));
						filteredMNXrxns=MNXrxns(trueHits);
				elseif find(strcmp('NA',mnx.Balance(b)))
						NAHits=find(strcmp('NA',mnx.Balance(b)));
						filteredMNXrxns=MNXrxns(NAHits);
				elseif find(strcmp('ambiguous',mnx.Balance(b)))
						amHits=find(strcmp('ambiguous',mnx.Balance(b)));
						filteredMNXrxns=MNXrxns(amHits);
				else
						filteredMNXrxns=MNXrxns;
				end
		end
end

author: SysBioChalmers

repo: Human-GEM


text:
function [HMDBId_new] = convertOld2NewHMDB(HMDBId)
% This function converts the old style HMDB ids to the new style
% old style 'HMDB06525'
% new style 'HMDB0006525' -- 7 digits -- fill up old ID to new ID with 0
%
% INPUT
% HMDBId    HMDB id
%
% OUTPUT
% HMDBId_new    new style HMDB id
%
% Ines Thiele 03/2022
HMDBId_new = '';
if length(HMDBId) < 12 % old style
    % remove HMDB part
    id = regexprep(HMDBId,'HMDB','');
    % add the missing
    z = '';
    for i = 1 : (7-length(id))
        z =[z '0'];
    end
    HMDBId_new = ['HMDB' z id];
end

author: opencobra

repo: cobratoolbox


text:
%compute speed in tsmhsme tailsm-headsm direction
function [data,units]=compute_veltsmhsm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
veltsmhsm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    veltsmhsm{1,i}=trx(larva).velmag.*(cos(trx(larva).velang).*cos(trx(larva).tailsmheadsmang(1,1:end-1))+sin(trx(larva).velang).*sin(trx(larva).tailsmheadsmang(1,1:end-1)));
end

units=parseunits('mm/s');
data=veltsmhsm;

author: kristinbranson

repo: JAABA


text:
function exact = p04_exact ( )

%*****************************************************************************80
%
%% P04_EXACT returns the estimated integral for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 0.47942822668880166736;

  return
end


text:
function [ h_p ] = rl_2D_plot_policy(X,Y,Q,actions,dx,dy,h_p)
%RL_2D_PLOT_POLICY
%
%   input -----------------------------------------------------------------
%
%       o X       : meshgrid X values.
%
%       o Y       : meshgrid Y values
%
%       o Q       : (num_states x num_actions)
%
%       o actions : (M x 2)
%
%       o h_p     : plot function handle
%

[~,idx] = max(Q,[],2);
uv      = 0.5 .* actions(idx,:);

if isempty(h_p)
    h_p = quiver(X,Y,uv(:,1),uv(:,2),0,'Color',[0 0 0],'Linewidth',2);
else

    set(h_p,'UData',uv(:,1),'VData',uv(:,2));

end

end

author: epfl-lasa


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the second part
%  of the exercise which covers regularization with logistic regression.
%
%  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 X values and the third column
%  contains the label (y).

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

plotData(X, y);

% Put some labels
hold on;

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

% Specified in plot order
legend('y = 1', 'y = 0')
hold off;

%% =========== Part 1: Regularized Logistic Regression ============
%  In this part, you are given a dataset with data points that are not
%  linearly separable. However, you would still like to use logistic
%  regression to classify the data points.
%
%  To do so, you introduce more features to use -- in particular, you add
%  polynomial features to our data matrix (similar to polynomial
%  regression).
%

% Add Polynomial Features

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1
lambda = 100;

% Compute and display initial cost and gradient for regularized logistic
% regression
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);

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

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

%% ============= Part 2: Regularization and Accuracies =============
%  Optional Exercise:
%  In this part, you will get to try different values of lambda and
%  see how regularization affects the decision coundart
%
%  Try the following values of lambda (0, 1, 10, 100).
%
%  How does the decision boundary change when you vary lambda? How does
%  the training set accuracy vary?
%

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1 (you should vary this)
lambda = 100;

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

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

% Plot Boundary
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

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

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

author: vugsus

repo: coursera-machine-learning


text:
function beta = l1plusTV(H, Ht, A, At, y, sizeImage, TVlambda, p, normfac, err, insweep, decfac, tol)

% Algorithm for solving problems of the form:
% min ||x||_p + TVlambda*TV(A'x)s.t. ||y-Hx||_2 < err

% Inputs
% H - Forward Measurement Operator
% Ht - Backward Measurement Operator
% A - Sparsifying Forward Transform
% At - Sparsifying Backward Transform
% y - collected data
% sizeImage - size of Image matrix
% TVlambda - multiplication factor for TV norm
% p - non-convex norm (default 1)
% normfac - highest singular value of A (default 1)
% err - two norm mismatch (default 1e-6)
% insweep - number of sweeps for internal loop (default 50)
% decfac - lambda decrease factor for outside loop (default 0.5)
% tol - tolerance for convergence (default 1e-4)

% copyright (c) Angshul Majumdar 2010

N = sizeImage;
c = 4;

if nargin < 9
    p = 1;
end
if nargin < 10
    normfac = 1;
end
if nargin < 11
    err = 1e-6;
end
if nargin < 12
    insweep = 50;
end
if nargin < 13
    decfac = 0.5;
end
if nargin < 14
    tol = 1e-4;
end

alpha = 1.1*normfac;
x_initial = Ht(y);
z_initial = zeros(length(TV(x_initial)),1);

x = x_initial; z = z_initial;

lambdaInit = decfac*max(abs(Ht(y))); lambda = lambdaInit;

f_current = norm(y-H(x_initial)) + lambda*norm(x,p) + TVlambda*norm(TV(At(x)),1);

while lambda > lambdaInit*tol
    lambda
    for i = 1:insweep
        f_previous = f_current;

        Wtilde = abs(x).^(p-2);
        W = 1 + (lambda/2)*Wtilde;

        b = x + (1/alpha)*(Ht(y-H(x)));
        btilde = (1./sqrt(W)).*b;
        z = (c*z + weightedTV(b-weightedTVtrans(z)))./((2*alpha/TVlambda)*(abs(weightedTV(x))).^(2-p)+c);
        x = 1./sqrt(W).*btilde - (1./sqrt(W)).*weightedTVtrans(z);

        f_current = norm(y-H(x_initial)) + lambda*norm(x,p) + TVlambda*norm(TV(x),1);

        if norm(f_current-f_previous)/norm(f_current + f_previous)<tol
            break;
        end
    end
    if norm(y-H(x))<err
        break;
    end
    lambda = decfac*lambda;
end
beta = x;

%%%%%% Functions for computing TV and TVtranspose %%%%%%%

    function wFvec = weightedTV(Xvec)
        wFvec = TV((1./sqrt(W)).*Xvec);
    end

    function wXvec = weightedTVtrans(Fvec)
        wXvec = (1./sqrt(W)).*TVtrans(Fvec);
    end

    function Fvec = TV(Xvec)
        Ivec = At(Xvec); % image pixels
        X = reshape(Ivec,N(1),N(2)); % sparse transform coefficients

        FX = zeros(N);
        FY = zeros(N);

        FX(1:end-1,:) = diff(X);
        FY(:,1:end-1) = diff(X')';

        Fvec = [FX(:)' FY(:)']';
    end

    function Xvec = TVtrans(Fvec)
        FX = reshape(Fvec(1:end/2),N(1),N(2));
        FY = reshape(Fvec(1+end/2:end),N(1),N(2));

        DX = zeros(N(1)+1,N(2));
        DY = zeros(N(1),N(2)+1);
        DX(2:end-1,:) = FX(1:end-1,:);
        DY(:,2:end-1) = FY(:,1:end-1);

        XT = -diff(DX);
        YT = -diff(DY')';

        X = XT + YT; % image
        Xvec = A(X(:));
    end
end

author: Sable

repo: mcbench-benchmarks


text:
classdef SIBEA < ALGORITHM
% <multi> <real/integer/label/binary/permutation>
% Simple indicator-based evolutionary algorithm

%------------------------------- Reference --------------------------------
% E. Zitzler, D. Brockhoff, and L. Thiele, The hypervolume indicator
% revisited: On the design of Pareto-compliant indicators via weighted
% integration, Proceedings of the International Conference on Evolutionary
% Multi-Criterion Optimization, 2007, 862-876.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Liangli Zhen

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

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = randi(length(Population),1,Problem.N);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [Ip] = go_prepare_image_Aplus(I, x1,y1,x2,y2,s,K1, K2)

Ic = imcrop(I,[x1 y1 x2 y2]);

I(y1,x1:x1+x2+1,:) = 255;
I(y1+1,x1:x1+x2+1,:) = 255;

I(y1+y2,x1:x1+x2+1,:) = 255;
I(y1+y2+1,x1:x1+x2+1,:) = 255;

I(y1:y1+y2,x1,:) = 255;
I(y1:y1+y2,x1+1,:) = 255;
I(y1:y1+y2,x1+x2,:) = 255;
I(y1:y1+y2,x1+x2+1,:) = 255;

Is = imresize(Ic,s,'nearest');

Ip = uint8(255*ones(size(I,1)+K1, size(I,2)+K2, 3));

Ip(1:size(I,1),1:size(I,2),:) = I;

Is(1:end,1,:) = 255;
Is(1:end,end,:) = 255;
Is(1,1:end,:) = 255;
Is(end,1:end,:) = 255;

Ip(size(Ip,1)-size(Is,1)+1:size(Ip,1),size(Ip,2)-size(Is,2)+1:size(Ip,2),:)=Is;

author: thomas-koehler

repo: SupER


text:
function [trl, conditionlabels, S] = spm_eeg_definetrial(S)
% Definition of trials based on events
% FORMAT [trl, conditionlabels, S] = spm_eeg_definetrial(S)
% S                 - input structure (optional)
% (optional) fields of S:
%   S.D             - MEEG object or filename of M/EEG mat-file
%   S.timewin       - time window {in PST ms}
%   S.trialdef      - structure array for trial definition with fields (optional)
%       S.trialdef.conditionlabel - string label for the condition
%       S.trialdef.eventtype      - string
%       S.trialdef.eventvalue     - string, numeric or empty
%       S.trialdef.trlshift       - shift the triggers by a fixed amount {ms}
%                                   (e.g. projector delay).
%   S.reviewtrials  - review individual trials after selection [yes/no: 1/0]
%   S.save          - save trial definition [yes/no: 1/0]
%
% OUTPUT:
%   trl             - Nx3 matrix [start end offset]
%   conditionlabels - Nx1 cell array of strings, label for each trial
%   S               - modified configuration structure (for history)
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak, Robert Oostenveld
% $Id: spm_eeg_definetrial.m 7132 2017-07-10 16:22:58Z guillaume $

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm_figure('FindWin','Interactive'); spm_clf('Interactive');

%-Get MEEG object
%--------------------------------------------------------------------------
try
    D = S.D;
catch
    [D, sts] = spm_select(1, 'mat', 'Select M/EEG mat file');
    if ~sts, trl = []; conditionlabels = {}; S = []; return; end
    S.D = D;
end

D = spm_eeg_load(D);

%-Get input parameters
%--------------------------------------------------------------------------
if ~isequal(D.type, 'continuous')
    error('Trial definition requires continuous dataset as input.');
end

event    = events(D, 1, 'samples');
fsample  = D.fsample;

if isempty(event)
    error('No event information was found in the input.');
end

if ~isfield(S, 'timewin')
    S.timewin = spm_input('Time window (ms)', '+1', 'r', [], 2);
end

pretrig  = S.timewin(1);
posttrig = S.timewin(2);

if ~isfield(S, 'trialdef')
    S.trialdef = [];
    ncond = spm_input('How many conditions?', '+1', 'n', '1');
    for i = 1:ncond
        OK = false;
        pos = '+1';
        while ~OK
            conditionlabel = spm_input(['Label of condition ' num2str(i)], pos, 's');
            selected = spm_eeg_select_event_ui(event);
            if isempty(conditionlabel) || isempty(selected)
                pos = '-1';
            else
                shift = spm_input('Shift triggers (ms)', pos, 'r', '0');
                for j = 1:size(selected, 1)
                    S.trialdef = [S.trialdef ...
                        struct('conditionlabel', conditionlabel, ...
                        'eventtype', selected{j, 1}, ...
                        'eventvalue', selected{j, 2}, ...
                        'trlshift', shift)];
                    OK = true;
                end
            end
        end
    end
end

for i = 1:length(S.trialdef)
    if ~isfield(S.trialdef(i),'trlshift')
        trlshift(i) = 0;
    else
        trlshift(i) = round(S.trialdef(i).trlshift * fsample/1000); % assume passed as ms
    end
end

%-Build trl based on selected events
%--------------------------------------------------------------------------
trl = [];
conditionlabels = {};
for i=1:numel(S.trialdef)

    if ischar(S.trialdef(i).eventvalue)
        % convert single string into cell-array, otherwise the intersection does not work as intended
        S.trialdef(i).eventvalue = {S.trialdef(i).eventvalue};
    end

    sel = [];
    % select all events of the specified type and with the specified value
    for j=find(strcmp(S.trialdef(i).eventtype, {event.type}))
        if isempty(S.trialdef(i).eventvalue)
            sel = [sel j];
        elseif ~isempty(intersect(event(j).value, S.trialdef(i).eventvalue))
            sel = [sel j];
        end
    end

    for j=1:length(sel)
        % override the offset of the event
        trloff = round(0.001*pretrig*fsample);
        % also shift the begin sample with the specified amount
        if ismember(event(sel(j)).type, {'trial', 'average'})
            % In case of trial events treat the 0 time point as time of the
            % event rather than the beginning of the trial
            trlbeg = event(sel(j)).sample - event(sel(j)).offset + trloff;
        else
            trlbeg = event(sel(j)).sample + trloff;
        end
        trldur = round(0.001*(-pretrig+posttrig)*fsample);
        trlend = trlbeg + trldur;

        % Added by Rik in case wish to shift triggers (e.g, due to a delay
        % between trigger and visual/auditory stimulus reaching subject).
        trlbeg = trlbeg + trlshift(i);
        trlend = trlend + trlshift(i);

        % add the beginsample, endsample and offset of this trial to the list
        trl = [trl; trlbeg trlend trloff];
        conditionlabels{end+1} = S.trialdef(i).conditionlabel;
    end
end

%-Sort the trl in right temporal order
%--------------------------------------------------------------------------
if isempty(trl)
    warning('No trials found.');
else
    [junk, sortind] = sort(trl(:,1));
    trl             = trl(sortind, :);
    conditionlabels = conditionlabels(sortind);
end

%-Review selected trials
%--------------------------------------------------------------------------
if ~isfield(S, 'reviewtrials')
    S.reviewtrials = spm_input('Review individual trials?','+1','yes|no',[1 0], 0);
end

if S.reviewtrials && ~isempty(trl)
    eventstrings = cell(size(trl,1),1);
    for i=1:size(trl,1)
        eventstrings{i} = [num2str(i) ' Label: ' conditionlabels{i} ' Time (sec): ' num2str((trl(i, 1)- trl(i, 3))./fsample)];
    end

    selected = find(trl(:,1)>0);

    [indx,ok] = listdlg('ListString', eventstrings, 'SelectionMode', 'multiple', 'InitialValue', ...
        selected, 'Name', 'Select events', 'ListSize', [300 300]);

    if ok
        trl = trl(indx, :);
        conditionlabels = conditionlabels(indx);
    end
end

%-Create trial definition file
%--------------------------------------------------------------------------
if ~isfield(S, 'save')
    S.save = spm_input('Save trial definition?','+1','yes|no',[1 0], 0);
end

if S.save
    [trlfilename, trlpathname] = uiputfile( ...
        {'*.mat', 'MATLAB File (*.mat)'}, 'Save trial definition as');

    if ~isequal(trlfilename,0) && ~isequal(trlpathname,0)
        trialdef = S.trialdef;
        timewin  = S.timewin;
        source   = D.fname;

        save(fullfile(trlpathname, trlfilename),...
            'trl', 'conditionlabels', 'trialdef', 'source', 'timewin', ...
            spm_get_defaults('mat.format'));
    end
end

%-Cleanup
%--------------------------------------------------------------------------
spm_clf('Interactive');

author: spm


text:
function test_gppde_2d

% $$$ x_g = -5:0.5:5;
% $$$ [x1_g, x2_g] = meshgrid(
T = 1;
L = 1;
%[x1_g, x2_g] = meshgrid(linspace(-0*T,3*T,60), linspace(-1*L,2*L,40));
min_L = 0.0*L;
max_L = 1.0*L;
min_T = 0.0*T;
max_T = 1.0*T;
[x1_g, x2_g] = meshgrid(linspace(min_T,max_T,40), linspace(min_L,max_L,40));
X_g = [x1_g(:), x2_g(:)];
[x1_f, x2_f] = meshgrid(linspace(0,T,45), linspace(0,L,45));
X_f = [x1_f(:), x2_f(:)];

% Define the partial differential equation
N_g = size(X_g, 1);
switch 1
 case 1
  %
  % HEAT EQUATION (d_t - d_x^2 = 0)
  %
  theta = [1 1/10]; % [scale, length scale]
  s2_y = theta(1)^2 * 1e-5;

  D = [1 0;     % differentiation
       0 2];
  alpha = [-1;   % coefficients
           0.02]; % (<- heat conduction)

  % Spatial boundary fixed to zero
  N_y = 30;
  bound1 = [linspace(min_T,max_T,N_y)', min_L*ones(N_y,1)];
  bound2 = [linspace(min_T,max_T,N_y)', max_L*ones(N_y,1)];
  X_y = [bound1; bound2];
  y = zeros(size(X_y,1),1);
  % Initial function at time T=min_T
  T_obs = min_T + 0.1*(max_T-min_T);
  bound_init = [T_obs*ones(N_y,1), linspace(min_L,max_L,N_y)'];
  X_y = [X_y; bound_init];
  y_init = 2*exp( -0.5*(bound_init(:,2)-0.5*(min_L+max_L)).^2 / 0.1^2);
  y = [y; y_init] + sqrt(s2_y)*randn(size(X_y,1),1);

  equation = 'Heat equation';

 case 2
  %
  % WAVE EQUATION (d_t^2 - d_x^2 = 0)
  %
  theta = [1 1/30]; % [scale, length scale]
  s2_y = theta(1)^2 * 1e-4;

  D = [2 0;
       0 2];
  alpha = [-1;
           16]; % speed ^ 2
  % Spatial boundary fixed to zero
  N_y = 100;
  bound1 = [linspace(min_T,max_T,N_y)', min_L*ones(N_y,1)];
  bound2 = [linspace(min_T,max_T,N_y)', max_L*ones(N_y,1)];
  X_y = [bound1; bound2];
  y = zeros(size(X_y,1),1);
% $$$   % Initial function at time T=min_T
% $$$   bound_init = [min_T*ones(N_y,1), linspace(min_L,max_L,N_y)'];
% $$$   X_y = [X_y; bound_init];
% $$$   y_init = 2*exp( -0.5*(bound_init(:,2)-0.5*(min_L+max_L)).^2 / 0.03^2);
% $$$   y = [y; y_init];
  y = y + sqrt(s2_y)*randn(size(X_y,1),1);

  equation = 'Wave equation';

 case 3
  %
  % LAPLACE EQUATION (d_x^2 + d_y^2 = 0)
  %
  theta = [1 1/20]; % [scale, length scale]
  s2_y = theta(1)^2 * 1e-5;

  D = [2 0;
       0 2];
  alpha = [1;
           1];
  X_y = zeros(0,2);
  y = zeros(0,1);

end
g = zeros(N_g,1);

% Inference: Get posterior predictive distribution
[m_f, V_f] = gppde(X_y, y, s2_y, theta, X_g, g, alpha, D, X_f);
% $$$ figure
% $$$ imagesc(V_f);

% Draw posterior predictive samples
samples = 1;
I_f = speye(length(V_f));
s2_f = 1e-8; % numerical reasons
L_f = chol(V_f + s2_f*I_f, 'lower');
f = bsxfun(@plus, L_f * randn(length(L_f),samples), m_f);
figure
for i = 1:samples
  subplot(ceil(sqrt(samples)), ceil(sqrt(samples)), i);
%  pcolor(reshape(f(:,i), size(x1_f)));
  contourf(x1_f, x2_f, reshape(f(:,i), size(x1_f)), 100);
  shading('flat')
  map_colormap();
  cl = max( abs(f(:,i)) );
  set(gca, 'clim', [-cl, cl])
  title(equation)
  xlabel('time');
  xlabel('time');
  ylabel('location');
  set(gca, 'xtick', [], 'ytick', []);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function copyLossConvertToClipboard(lossConvert)
%copyLossConvertToClipboard Summary of this function goes here
%   Detailed explanation goes here
    global GLOBAL_LossConvertClipboard;

    GLOBAL_LossConvertClipboard = lossConvert;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [vertexIndices,distSq] = GetPointsWithinSphere(point,radius,xyzData)

% CHECK INPUTS
assert(numel(point) == 3,   'Point must be specified as a 3D point [3x1]');
assert(size(xyzData,2) == 3,'Vertex data must be specified a list of 3D vertices [nx3]');
assert(numel(radius) == 1,  'The radius is a scalar comparitor.');

% DEFAULT VALUES
% internalPoints = xyzData;               % Return the full set
vertexIndices = 1:1:size(xyzData,1);      % Return a substitute index vector
distSq = inf(1,3);

% DEFAULT RESPONSE
if isinf(radius)
    return
end
% VECTORISED SQUARED DISTANCE CALCULATION
distSq = (xyzData(:,1) - point(1)).^2 + (xyzData(:,2) - point(2)).^2 + (xyzData(:,3) - point(3)).^2;
% SQUARED DISTANCE COMPARITOR
ind = distSq < radius^2;
% THE REDUCED POINT CLOUD INSIDE THE RADIUS
[vertexIndices,~] = find(ind);
end


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%

function variable

% small number of numeric value comparison
mx = false;
eps = 1e-5;

v1 = [5;2.8;2];
v2 = [1,5.01,3;6,5,4];

%%% constructor
a1 = ocl.casadi.CasadiVariable.Matrix([3,1],mx);
a2 = ocl.casadi.CasadiVariable.Matrix([2,3],mx);
s1 = a1.value;
s2 = a2.value;

f = casadi.Function('f',{s1},{s1});
assert(isequal(full(f(v1)),v1))
f = casadi.Function('f',{s2},{s2});
assert(isequal(full(f(v2)),v2))

%%% horzcat, vertcat, subsref
aTest = [a1,a1,a1;a2];
aTest = aTest(3:5,2);
vTest = [v1,v1,v1;v2];
vTest = vTest(3:5,2);
f = casadi.Function('f',{s1,s2},{aTest(2:3).value});
assert(isequal(full(f(v1,v2)),vTest(2:3)))

%%% subsasgn
aTest(2:3) = 11;
vTest(2:3) = 11;
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% norm
aTest = norm(aTest);
vTest = norm(vTest);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% uplus, uminus
%aTest = -+(+aTest);
%vTest = -+(+vTest);
%f = casadi.Function('f',{s1,s2},{aTest.value});
%assert(isequal(full(f(v1,v2)),vTest))

%%% sum
aTest = [a1,a1,a1;a2];
aTest = sum(aTest);
vTest = [v1,v1,v1;v2];
vTest = sum(vTest);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% mtimes
aTest = [a1,a1,a1;a2] * 4 * s1 * [2,3,4];
vTest = [v1,v1,v1;v2] * 4 * v1 * [2,3,4];
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% transpose
aTest = aTest.';
vTest = vTest.';
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% reshape
aTest = reshape(aTest,15,1);
vTest = reshape(vTest,15,1);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))
aTest = reshape(aTest,[3,5]);
vTest = reshape(vTest,[3,5]);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% triu
aTest = triu(aTest);
vTest = triu(vTest);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% repmat
aTest = repmat(aTest,2,3);
vTest = repmat(vTest,2,3);
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% mpower
% mpower is supported since casadi 3.3
aTest = aTest(:,1:6)^2;
vTest = vTest(:,1:6)^2;
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(all(all(vTest-full(f(v1,v2))< eps)))

%%% mldivide (solve in casadi)
A = [0.2625    0.9289    0.5785;
     0.8010    0.7303    0.2373;
     0.0292    0.4886    0.4588];
b = [0.9631,0.5468,0.5211]';

aA = ocl.casadi.CasadiVariable.Matrix([3,3],mx);
ab = ocl.casadi.CasadiVariable.Matrix([3,1],mx);

sA = aA.value;
sb = ab.value;

aTest = aA\ab;
vTest = A\b;
f = casadi.Function('f',{sA,sb},{aTest.value});
assert(all( abs(full(f(A,b))-vTest) <= eps))

%%% cross
aTest = cross(a1,a2(2,:).');
vTest = cross(v1,v2(2,:).');
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% dot
aTest = dot(a1,a2(2,:).');
vTest = dot(v1,v2(2,:).');
f = casadi.Function('f',{s1,s2},{aTest.value});
assert(isequal(full(f(v1,v2)),vTest))

%%% inv
aTest = inv(aA);
vTest = inv(A);
f = casadi.Function('f',{sA},{aTest.value});
assert( all(all( abs(full(f(A))-vTest) <= eps)))

%%% det
aTest = det(aA);
vTest = det(A);
f = casadi.Function('f',{sA},{aTest.value});
assert( abs(full(f(A))-vTest) <= eps)

%%% trace
aTest = trace(aA);
vTest = trace(A);
f = casadi.Function('f',{sA},{aTest.value});
assert(isequal(full(f(A)),vTest))

%%% diag
aTest = diag(aA);
vTest = diag(A);
f = casadi.Function('f',{sA},{aTest.value});
assert(isequal(full(f(A)),vTest))

%%% polyval
aTest = polyval([2;5;4],a1);
vTest = polyval([2;5;4],v1);
f = casadi.Function('f',{s1},{aTest.value});
assert(isequal(full(f(v1)),vTest))

ab.set([2,5,4].');
aTest = polyval(ab,a1);
vTest = polyval([2,5,4].',v1);
f = casadi.Function('f',{s1},{aTest.value});
assert(isequal(full(f(v1)),vTest))

%%% jacobian
% test against finite diff jacobian
aTest = jacobian(testJacobianFun(a1),a1);
vTest = finiteDiffJac(@testJacobianFun,v1);
f = casadi.Function('f',{s1},{aTest.value});
assert(all(all( abs(full(f(v1))-vTest) <= eps)))

%%% jtimes
aTest = jtimes(testJacobianFun(a1),a1,a1);
vTest = vTest * v1;
f = casadi.Function('f',{s1},{aTest.value});
assert(all(all( abs(full(f(v1))-vTest) <= 10*eps)))

%%% plus,minus,times,power,rdivide
aTest = a1(1).^2*a1(2)+a1(3)-a1(2)./a1(3).\a1(1);
vTest = v1(1).^2*v1(2)+v1(3)-v1(2)./v1(3).\v1(1);
f = casadi.Function('f',{s1},{aTest.value});
assert(isequal(full(f(v1)),vTest))

%%% abs,sqrt,sin,cos,tan,atan,asin,acos,atanh,asinh,acosh,exp,log,tanh,cosh,sinh
aTest = tanh(acosh(atan(a1(1)).^2));
aTest = sinh(aTest * sqrt(a1(2))+asinh(abs(a1(3))-log(sin(a1(2)))));
aTest = cosh(acos(asin(aTest * exp(atanh(cos(a1(3)))).\tan(a1(1))+1)));

vTest = tanh(acosh(atan(v1(1)).^2));
vTest = sinh(vTest * sqrt(v1(2))+asinh(abs(v1(3))-log(sin(v1(2)))));
vTest = cosh(acos(asin(vTest * exp(atanh(cos(v1(3)))).\tan(v1(1))+1)));

f = casadi.Function('f',{s1},{aTest.value});
assert(abs(full(f(v1)) - vTest) < eps)

%%% atan2, times
aTest = atan2(a1(1).^2.*a1(2)+a1(3),atan(a1(1)).^2);
vTest = atan2(v1(1).^2.*v1(2)+v1(3),atan(v1(1)).^2);
f = casadi.Function('f',{s1},{aTest.value});
assert(isequal(full(f(v1)),vTest))

%%% (:)
%aTest = [a1,a1,a1;a2];
%vTest = [v1,v1,v1;v2];
%aTest(:) = 2;
%vTest(:) = 2;
%f = casadi.Function('f',{s1,s2},{aTest.value});
%assert(isequal(full(f(v1,v2)),vTest))

% keep type
assert(isa(aTest,'ocl.casadi.CasadiVariable'));

end

function v = testJacobianFun(x)
  v = x*x(1)+cross([x(1);x(3)^2;x(2)],x);
end

function J = finiteDiffJac(functionHandle,x)
  FDeps = 1e-6;
  fx = functionHandle(x);
  Ncols = numel(x);
  Nrows = numel(fx);
  J = zeros(Nrows,Ncols);
  for k=1:Ncols
    dx = zeros(Ncols,1);
    dx(k) = FDeps;
    J(:,k) = (functionHandle(x+dx) - fx) / FDeps;
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text: %PLUS Mapping overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% batch
totaltime = tic; %#ok<NASGU>
isFullData = 1;
data_masterdir = GetCurrentDataDir();

M_stimrange = GetStimRange();

range_fish =  8:18; % range_fish = GetFishRange();

%% custom params here:
% numK1 = 20;
masterthres = 0.7;

%%
M_regthres = {0.7,0.5};
M_place = {1,2,3,1,4,5,6,7};
M_stimname = {'4x4','PT','OMR','defS','Spt','DF','Lm','Dot'};
%%
for i_count = 1,%%%%%%%%%
    masterthres = M_regthres{i_count};
%     clusParams = struct('merge',masterthres,'cap',masterthres,'reg1',masterthres,...
%         'reg2',masterthres,'minSize',10,'k1',numK1);

    for i_stimrange = 1:8,
        if i_stimrange == 1,
            M_stimrange = GetStimRange('5');
        elseif i_stimrange == 2,
            M_stimrange = GetStimRange('P');
        elseif i_stimrange == 3,
            M_stimrange = GetStimRange('O');
        elseif i_stimrange == 4,
            M_stimrange = GetStimRange('M');
        elseif i_stimrange == 5,
            M_stimrange = GetStimRange('S');
        elseif i_stimrange == 6,
            M_stimrange = GetStimRange('D');
        elseif i_stimrange == 7,
            M_stimrange = GetStimRange('L');
        elseif i_stimrange == 8,
            M_stimrange = GetStimRange('Y');
        end

        for i = 1:length(range_fish),
            i_fish = range_fish(i);
            disp(i_fish);

            % check this loop
            stimrange = M_stimrange{i_fish};
            if isempty(stimrange),
                continue;
            end

            % Load fish
            LoadFullFish(hfig,i_fish,isFullData);

            %% 1.
            % setup
            absIX = getappdata(hfig,'absIX');

            i_ClusGroup = 2;
            i_Cluster = 1;

            % Load cluster data
            [cIX_load,gIX] = LoadCluster_Direct(i_fish,i_ClusGroup,i_Cluster,absIX);

            %% partitions for CV
            timelists = getappdata(hfig,'timelists');
            timelists_names = getappdata(hfig,'timelists_names');
            periods = getappdata(hfig,'periods');

            M_stim = M_stimrange{i_fish};

            timelistsCV_raw = cell(length(M_stim),2);
            timelistsCV = cell(1,2);

            for k_stim = 1:length(M_stim), % :3
                i_stim = M_stim(k_stim);
                TL = timelists{i_stim};
                period = periods(i_stim);
                nrep = size(TL,2)/periods(i_stim); % integer
                n = floor(nrep/2);
                if n>0,
                    timelistsCV_raw{k_stim,1} = TL(1:n*period);
                    timelistsCV_raw{k_stim,2} = TL(1+n*period:2*n*period);% before 12/5/16: TL(1+n*period):TL(2*n*period);
                else % for spont, only one period
                    halfperiod = floor(period/2);
                    timelistsCV_raw{k_stim,1} = TL(1:halfperiod);
                    timelistsCV_raw{k_stim,2} = TL(1+halfperiod:2*halfperiod);
                end
            end
            timelistsCV{1} = horzcat(timelistsCV_raw{:,1});
            timelistsCV{2} = horzcat(timelistsCV_raw{:,2});
            assert(length(timelistsCV{1})==length(timelistsCV{2}));

            %%
            for k = 1:2,% CV halves
                tIX = timelistsCV{k};
                M = GetTimeIndexedData_Default_Direct(hfig,cIX_load,tIX);
                M_0 = GetTimeIndexedData_Default_Direct(hfig,[],tIX,'isAllCells');

                % ------custom code here---------
                isWkmeans = true;
                isMakeFoxels = true;

                [cIX,gIX] = AutoClustering(cIX_load,gIX,M_0,cIX_load,isWkmeans,[],...
                    isMakeFoxels,masterthres);

                % save cluster
                name = ['Auto_',M_stimname{i_stimrange},'_M',num2str(masterthres),'_CV',num2str(k)];
                clusgroupID = 3+k; % 4 and 5
                clusIDoverride = M_place{i_stimrange};
                SaveCluster_Direct(cIX,gIX,absIX,i_fish,name,clusgroupID,clusIDoverride);
            end
        end
    end
end
SaveVARwithBackup();
totaltime = toc;
disp(totaltime);

author: xiuyechen

repo: FishExplorer


text:
seed = 0;
rand('state', seed);
randn('state', seed);

chars = ['a', 'c', 'g', 't'];
motif = 'accca';
motif_length = length(motif);
motif_code = zeros(1, motif_length);
for i=1:motif_length
  motif_code(i) = find(chars == motif(i));
end

[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_length', length(motif));
%[bnet_init, Qnodes, Fnodes, Onode] = mk_motif_hhmm('motif_pattern', motif);
ss = bnet_init.nnodes_per_slice;

% We generate a training set by creating uniform sequences,
% and inserting a single motif at a random location.
ntrain = 100;
T = 20;
cases = cell(1, ntrain);

if 1
  % uniform background
  background_dist = normalise(ones(1, length(chars)));
end
if 0
  % use a constant background
  background_dist = zeros(1, length(chars));
  m = find(chars=='t');
  background_dist(m) = 1.0;
end
if 0
  % use a background skewed away from the motif
  p = 0.01; q = (1-(2*p))/2;
  background_dist = [p p q q];
end

unif_pos = normalise(ones(1, T-length(motif)));
cases = cell(1, ntrain);
data = zeros(1,T);
for i=1:ntrain
  data = sample_discrete(background_dist, 1, T);
  L = sample_discrete(unif_pos, 1, 1);
  data(L:L+length(motif)-1) = motif_code;
  cases{i} = cell(ss, T);
  cases{i}(Onode,:) = num2cell(data);
end
disp('sample training cases')
for i=1:5
  chars(cell2num(cases{i}(Onode,:)))
end

engine_init = hmm_inf_engine(bnet_init);

[bnet_learned, LL, engine_learned] = ...
    learn_params_dbn_em(engine_init, cases, 'max_iter', 100, 'thresh', 1e-2);
%			'anneal', 1, 'anneal_rate', 0.7);

% extract the learned motif profile
eclass = bnet_learned.equiv_class;
CPDO=struct(bnet_learned.CPD{eclass(Onode,1)});
fprintf('columns = chars, rows = states\n');
profile_learned = squeeze(CPDO.CPT(2,:,:))
[m,ndx] = max(profile_learned, [], 2);
map_motif_learned = chars(ndx)
back_learned = squeeze(CPDO.CPT(1,1,:))'
%map_back_learned = chars(argmax(back_learned))

CPDO_init = struct(bnet_init.CPD{eclass(Onode,1)});
profile_init = squeeze(CPDO_init.CPT(2,:,:));
back_init = squeeze(CPDO_init.CPT(1,1,:))';

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt