text:
function varargout = process_timefreq( varargin )
% PROCESS_TIMEFREQ: Computes the time frequency decomposition of any signal in the database.
%
% USAGE:  sProcess = process_timefreq('GetDescription')
%           sInput = process_timefreq('Run',     sProcess, sInput)
%           TFmask = process_timefreq('GetEdgeEffectMask', Time, Freqs, tfOptions)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Time-frequency (Morlet wavelets)';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Frequency';
    sProcess.Index       = 505;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#Morlet_wavelets';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'results', 'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Options: Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = 'MEG, EEG';
    sProcess.options.sensortypes.InputTypes = {'data'};
    sProcess.options.sensortypes.Group   = 'input';
    % Options: Scouts
    sProcess.options.clusters.Comment = '';
    sProcess.options.clusters.Type    = 'scout_confirm';
    sProcess.options.clusters.Value   = {};
    sProcess.options.clusters.InputTypes = {'results'};
    sProcess.options.clusters.Group   = 'input';
    % Options: Scout function
    sProcess.options.scoutfunc.Comment    = {'Mean', 'Max', 'PCA', 'Std', 'All', 'Scout function:'};
    sProcess.options.scoutfunc.Type       = 'radio_line';
    sProcess.options.scoutfunc.Value      = 1;
    sProcess.options.scoutfunc.InputTypes = {'results'};
    sProcess.options.scoutfunc.Group   = 'input';
    % Options: Time-freq
    sProcess.options.edit.Comment = {'panel_timefreq_options', 'Morlet wavelet options: '};
    sProcess.options.edit.Type    = 'editpref';
    sProcess.options.edit.Value   = [];
    % Options: Normalize
    sProcess.options.normalize2020.Comment = 'Spectral flattening: Multiply output power values by frequency';
    sProcess.options.normalize2020.Type    = 'checkbox';
    sProcess.options.normalize2020.Value   = 0;
    % Old normalize option, for backwards compatibility.
    sProcess.options.normalize.Comment = {'None: Save non-standardized time-frequency maps', '1/f compensation: Multiply output values by frequency'; ...
                                          'none', 'multiply'};
    sProcess.options.normalize.Type    = 'radio_label';
    sProcess.options.normalize.Value   = 'none';
    sProcess.options.normalize.Hidden  = 1;
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Initialize returned values
    OutputFiles = {};
    % Extract method name from the process name
    switch (func2str(sProcess.Function))
        case 'process_timefreq',      strProcess = 'morlet';
        case 'process_hilbert',       strProcess = 'hilbert';
        case 'process_fft',           strProcess = 'fft';
        case 'process_psd',           strProcess = 'psd';
        case 'process_sprint',        strProcess = 'sprint';
        case 'process_ft_mtmconvol',  strProcess = 'mtmconvol';
        otherwise,                    error('Unsupported process.');
    end
    % Get editable options (Edit... button)
    if isfield(sProcess.options, 'edit')
        tfOPTIONS = sProcess.options.edit.Value;
        % If user did not edit values: get default values
        if isempty(tfOPTIONS)
            [bstPanelNew, panelName] = panel_timefreq_options('CreatePanel', sProcess, sInputs);
            jPanel = gui_show(bstPanelNew, 'JavaWindow', panelName, 0, 0, 0);
            drawnow;
            tfOPTIONS = panel_timefreq_options('GetPanelContents');
            gui_hide(panelName);
        end
    % Else: get default options
    else
        tfOPTIONS = bst_timefreq();
        tfOPTIONS.Method = strProcess;
        switch tfOPTIONS.Method
            case 'fft',       tfOPTIONS.Comment = 'FFT';
            case 'psd',       tfOPTIONS.Comment = 'PSD';
            case 'sprint',    tfOPTIONS.Comment = 'SPRiNT';
            case 'morlet',    tfOPTIONS.Comment = 'Wavelet';
            case 'hilbert',   tfOPTIONS.Comment = 'Hilbert';
            case 'mtmconvol', tfOPTIONS.Comment = 'Multitaper';
        end
    end

    % Add other options
    tfOPTIONS.Method = strProcess;
    if isfield(sProcess.options, 'sensortypes')
        tfOPTIONS.SensorTypes = sProcess.options.sensortypes.Value;
    else
        tfOPTIONS.SensorTypes = [];
    end
    if isfield(sProcess.options, 'mirror') && ~isempty(sProcess.options.mirror) && ~isempty(sProcess.options.mirror.Value)
        tfOPTIONS.isMirror = sProcess.options.mirror.Value;
    else
        tfOPTIONS.isMirror = 0;
    end
    if isfield(sProcess.options, 'clusters') && ~isempty(sProcess.options.clusters) && ~isempty(sProcess.options.clusters.Value)
        tfOPTIONS.Clusters = sProcess.options.clusters.Value;
    else
        tfOPTIONS.Clusters = [];
    end
    % Override scouts function
    if isfield(sProcess.options, 'scoutfunc') && isfield(sProcess.options.scoutfunc, 'Value') && ~isempty(sProcess.options.scoutfunc.Value)
        switch lower(sProcess.options.scoutfunc.Value)
            case {1, 'mean'}, tfOPTIONS.ScoutFunc = 'mean';
            case {2, 'max'},  tfOPTIONS.ScoutFunc = 'max';
            case {3, 'pca'},  tfOPTIONS.ScoutFunc = 'pca';
            case {4, 'std'},  tfOPTIONS.ScoutFunc = 'std';
            case {5, 'all'},  tfOPTIONS.ScoutFunc = 'all';
            otherwise,  bst_report('Error', sProcess, [], 'Invalid scout function.');  return;
        end
    else
        tfOPTIONS.ScoutFunc = [];
    end
    % If a time window was specified
    if isfield(sProcess.options, 'timewindow') && ~isempty(sProcess.options.timewindow) && ~isempty(sProcess.options.timewindow.Value) && iscell(sProcess.options.timewindow.Value)
        tfOPTIONS.TimeWindow = sProcess.options.timewindow.Value{1};
    elseif ~isfield(tfOPTIONS, 'TimeWindow')
        tfOPTIONS.TimeWindow = [];
    end
    % If a window length was specified (PSD)
    if isfield(sProcess.options, 'win_length') && ~isempty(sProcess.options.win_length) && ~isempty(sProcess.options.win_length.Value) && iscell(sProcess.options.win_length.Value)
        tfOPTIONS.WinLength  = sProcess.options.win_length.Value{1};
        tfOPTIONS.WinOverlap = sProcess.options.win_overlap.Value{1};
    end
    if isfield(sProcess.options, 'win_std') && ~isempty(sProcess.options.win_std) && ~isempty(sProcess.options.win_std.Value)
        tfOPTIONS.WinStd = sProcess.options.win_std.Value;
        if tfOPTIONS.WinStd
            tfOPTIONS.Comment = [tfOPTIONS.Comment ' std'];
        end
    end
    % If units specified (PSD)
    if isfield(sProcess.options, 'units') && ~isempty(sProcess.options.units) && ~isempty(sProcess.options.units.Value)
        tfOPTIONS.PowerUnits = sProcess.options.units.Value;
    end
    % Multitaper options
    if isfield(sProcess.options, 'mt_taper') && ~isempty(sProcess.options.mt_taper) && ~isempty(sProcess.options.mt_taper.Value)
        if iscell(sProcess.options.mt_taper.Value)
            tfOPTIONS.ft_mtmconvol.taper = sProcess.options.mt_taper.Value{1};
        else
            tfOPTIONS.ft_mtmconvol.taper = sProcess.options.mt_taper.Value;
        end
    end
    if isfield(sProcess.options, 'mt_frequencies') && ~isempty(sProcess.options.mt_frequencies) && ~isempty(sProcess.options.mt_frequencies.Value)
        tfOPTIONS.ft_mtmconvol.frequencies = eval(sProcess.options.mt_frequencies.Value);
        % Add frequencies to comment
        tfOPTIONS.Comment = [tfOPTIONS.Comment, ' ', sProcess.options.mt_frequencies.Value, 'Hz'];
    end
    if isfield(sProcess.options, 'mt_freqmod') && ~isempty(sProcess.options.mt_freqmod) && ~isempty(sProcess.options.mt_freqmod.Value)
        tfOPTIONS.ft_mtmconvol.freqmod = sProcess.options.mt_freqmod.Value{1};
    end
    if isfield(sProcess.options, 'mt_timeres') && ~isempty(sProcess.options.mt_timeres) && ~isempty(sProcess.options.mt_timeres.Value)
        tfOPTIONS.ft_mtmconvol.timeres = sProcess.options.mt_timeres.Value{1};
    end
    if isfield(sProcess.options, 'mt_timestep') && ~isempty(sProcess.options.mt_timestep) && ~isempty(sProcess.options.mt_timestep.Value)
        tfOPTIONS.ft_mtmconvol.timestep = sProcess.options.mt_timestep.Value{1};
    end
    if isfield(sProcess.options, 'measure') && ~isempty(sProcess.options.measure) && ~isempty(sProcess.options.measure.Value)
        tfOPTIONS.Measure = sProcess.options.measure.Value;
        % Add measure to comment
        if strcmpi(tfOPTIONS.Measure, 'none')
            strMeasure = 'complex';
        else
            strMeasure = tfOPTIONS.Measure;
        end
        tfOPTIONS.Comment = [tfOPTIONS.Comment, ' ', strMeasure];
    end
    % if process is SPRiNT
    if isfield(sProcess.options, 'fooof')
       tfOPTIONS.SPRiNTopts = sProcess.options;
    end
    % Output
    if isfield(sProcess.options, 'avgoutput') && ~isempty(sProcess.options.avgoutput) && ~isempty(sProcess.options.avgoutput.Value)
        if sProcess.options.avgoutput.Value
            tfOPTIONS.Output = 'average';
        else
            tfOPTIONS.Output = 'all';
        end
    end
    % Frequency normalization
    if isfield(sProcess.options, 'normalize2020') && ~isempty(sProcess.options.normalize2020)
        if isequal(sProcess.options.normalize2020.Value, 1)
            tfOPTIONS.NormalizeFunc = 'multiply2020';
        elseif ischar(sProcess.options.normalize2020.Value)
            tfOPTIONS.NormalizeFunc = sProcess.options.normalize2020.Value;
        end
    elseif isfield(sProcess.options, 'normalize') && ~isempty(sProcess.options.normalize)
        if isequal(sProcess.options.normalize.Value, 1)
            tfOPTIONS.NormalizeFunc = 'multiply';
        elseif ischar(sProcess.options.normalize.Value)
            tfOPTIONS.NormalizeFunc = sProcess.options.normalize.Value;
        end
    else
        tfOPTIONS.NormalizeFunc = 'none';
    end

    % === EXTRACT CLUSTER/SCOUTS ===
    if ~isempty(tfOPTIONS.Clusters)
        % If cluster function should be applied AFTER time-freq: get all time series
        if strcmpi(tfOPTIONS.ClusterFuncTime, 'after')
            ExtractScoutFunc = 'all';
        else
            ExtractScoutFunc = tfOPTIONS.ScoutFunc;
        end
        AddRowComment = ~isempty(tfOPTIONS.ScoutFunc) && strcmpi(tfOPTIONS.ScoutFunc, 'all');
        % Flip sign only for results
        isflip = strcmpi(sInputs(1).FileType, 'results');
        % Call process
        ClustMat = bst_process('CallProcess', 'process_extract_scout', sInputs, [], ...
            'timewindow',     tfOPTIONS.TimeWindow, ...
            'scouts',         tfOPTIONS.Clusters, ...
            'scoutfunc',      ExtractScoutFunc, ...  % If ScoutFunc is not defined, use the scout function available in each scout
            'isflip',         isflip, ...
            'isnorm',         0, ...
            'concatenate',    0, ...
            'save',           0, ...
            'addrowcomment',  AddRowComment, ...
            'addfilecomment', 0, ...
            'progressbar',    0);
        if isempty(ClustMat)
            bst_report('Error', sProcess, sInputs, 'Cannot access clusters/scouts time series.');
            return;
        end
        % Get data to process
        DataToProcess = {ClustMat.Value};
        tfOPTIONS.TimeVector  = ClustMat(1).Time;
        tfOPTIONS.ListFiles   = {sInputs.FileName};
        tfOPTIONS.nComponents = [ClustMat.nComponents];
        tfOPTIONS.RowNames    = {};
        tfOPTIONS.SurfaceFile = {ClustMat.SurfaceFile};
        for iFile = 1:length(ClustMat)
             tfOPTIONS.RowNames{iFile} = ClustMat(iFile).Description;
        end
        clear ClustMat;
    % === DATA FILES ===
    else
        DataToProcess = {sInputs.FileName};
        tfOPTIONS.TimeVector = in_bst(sInputs(1).FileName, 'Time');
    end

    % === OUTPUT STUDY ===
    if strcmpi(tfOPTIONS.Output, 'average')
        % Get output study
        [sStudy, iStudy, Comment] = bst_process('GetOutputStudy', sProcess, sInputs);
        % If no valid output study can be found
        if isempty(iStudy)
            return;
        end
        % Save all outputs from bst_timefreq in target Study
        tfOPTIONS.iTargetStudy = iStudy;
    else
        tfOPTIONS.iTargetStudy = [];
    end

    % === START COMPUTATION ===
    [OutputFiles, Messages, isError] = bst_timefreq(DataToProcess, tfOPTIONS);
    if ~isempty(Messages)
        if isError
            bst_report('Error', sProcess, sInputs, Messages);
        elseif isempty(OutputFiles)
            bst_report('Warning', sProcess, sInputs, Messages);
        else
            bst_report('Info', sProcess, sInputs, Messages);
            disp(['BST> process_timefreq: ' Messages]);
        end
    end
end

%% ===== GET EDGE EFFECT MASK =====
function TFmask = GetEdgeEffectMask(Time, Freqs, tfOptions) %#ok<DEFNU>
    % Get time vector
    if ~iscell(Time)
        t = Time;
    else
        t = mean(cell2mat(Time(:,2:3)), 2)';
    end
    % Get frequency vector
    if ~iscell(Freqs)
        f = Freqs;
    else
        FreqBands = process_tf_bands('GetBounds', Freqs);
        f = mean(FreqBands, 2)';
    end

    % Morlet wavelets
    if isfield(tfOptions, 'Method') && strcmpi(tfOptions.Method, 'morlet') && isfield(tfOptions, 'MorletFc') && isfield(tfOptions, 'MorletFwhmTc')
        FWHM_t = tfOptions.MorletFwhmTc * tfOptions.MorletFc ./ f;
        TF_timeres = repmat(FWHM_t' ./ 2, [1,length(t)]);
        TFmask = (bst_bsxfun(@minus, t - t(1), TF_timeres) > 0);
        TFmask = (TFmask & bst_flip(TFmask, 2));

    % Hilbert transform
    elseif isfield(tfOptions, 'Method') && strcmpi(tfOptions.Method, 'hilbert')
        if iscell(Time)
            disp('BST> Edge effects map for Hilbert process is not available yet when time bands are selected.');
            TFmask = [];
        else
            % Sampling frequency
            sfreq = 1 ./ (t(2) - t(1));
            % Compute the transients for each bandpass filter
            TFmask = zeros(size(Freqs,1), length(t));
            for i = 1:size(Freqs,1)
                % Compute the filter specifications
                [tmp, FiltSpec] = process_bandpass('Compute', [], sfreq, FreqBands(i,1), FreqBands(i,2), 'bst-hfilter-2019');
                % Only the values outside of the transients are valid
                TFmask(i,(t - t(1) > FiltSpec.transient) & (t(end) - t > FiltSpec.transient)) = 1;
            end
        end
    else
        TFmask = [];
    end
end

author: brainstorm-tools


text:
% gigascience -> EMG signal figure

clear all; clc; close all;

dd='C:\Users\Doyeunlee\Desktop\Analysis\01_Raw data\';
% filelist={'sub01','sub02','sub03','sub04','sub05','sub06','sub07','sub08','sub09','sub10'};
filelist={'session1_sub4_multigrasp_MI_EMG'};
ival=[0 4500];

% selected_class = [11 12 13 15 16 17 18 20 21 22 23 43 44 45 47 48 49 51 52 53];
% selected_class = {'F3','F1','Fz','F2','F4','FC3','FC1','FCz','FC2','FC4','C3','C1', 'Cz', 'C2', 'C4','CP3','CP1','CPz','CP2','CP4','P3','P1','Pz','P2','P4'};
% selected_class = {'C3','C1', 'Cz', 'C2', 'C4'};
% selected_class = {'Cz'};
selected_class ={'EMG_1','EMG_2','EMG_3','EMG_4','EMG_5','EMG_6','EMG_ref'};
% selected_class = [1 2 3 4 5 6 7];

fold = 5;

for sub = 1:length(filelist)
    [cnt,mrk,mnt]=eegfile_loadMatlab([dd '\' filelist{sub}]);

    cnt=proc_selectChannels(cnt,selected_class);
    cnt=proc_filtButter(cnt,5,[10 500]);
    epo=cntToEpo(cnt,mrk,ival);

    classes=size(epo.className,2);
    trial=size(epo.x,3)/2/(classes-1);

    eachClassFold_no=trial/fold;

end

fs = 250;
% plot(psd(spectrum.periodogram, cntReach.x, 'Fs', fs));
% psdest = psd(spectrum.periodogram, cntReach.x, 'Fs',fs);
% plot(psdest.Frequencies,psdest.Data);
% xlabel('Hz'); grid on;
% ylim([0 60]);
% xlim([0 60]);
% hold on;
%
% fs = 250;
% figure();
% plot(psd(spectrum.periodogram, cntReach_ME.x, 'Fs', fs));
% psdest = psd(spectrum.periodogram, cntReach_ME.x, 'Fs',fs);
% plot(psdest.Frequencies,psdest.Data);
% xlabel('Hz'); grid on;
% ylim([0 60]);
% xlim([0 60]);
%
% xdft = fft(cntReach_ME.x);
% xdft = xdft(1:length(cntReach_ME.x)/2+1);
% freq = 0:fs/length(cntReach_ME.x):fs/2;
% figure();
% plot(freq,abs(xdft));
% xlabel('Hz');

% T = 1/fs;
% tx = [0:length(cnt.x)-1]/fs;
% figure();
% subplot(211);
% plot(tx,cnt.x); xlabel('Time(s)'); ylabel('Amplitude(uV)');

% mean_EEGsig = mean(cnt.x);
% max_value = max(cnt.x);
% mean_value = mean(cnt.x);
% threshold=(max_value - mean_value)/2;
% subplot(212), plot(tx, mean_value);
% xlabel('Time (s)'); ylabel('Amplitude(uV)');

% T = 1/fs;
% tx = [0:length(cntReach.x)-1]/fs;
% figure();
% % subplot(211);
% plot(tx,cntReach.x); xlabel('Time(s)'); ylabel('Amplitude(uV)');

epo_1 = cntToEpo(cnt, mrk, ival);
fv_1 = proc_rectifyChannels(epo_1);
fv_1 = proc_movingAverage(fv_1, 100, 'centered');
fv_1 = proc_baseline(fv_1, [0 500]);

result = mean(fv_1.x,3);
a = result;

X = abs(a);
X_1 = movmean(X, 30);

plot(X_1);
% plot(X_1,'LineWidth',1.5);
% ylim([0 200]);
xticks([0 1250 3750 6250 8750 11250]);
xlim([0 11251]);
% ylim([0 35]); % reaching
ylim([0 35]); % grasp
% ylim([0 10]); % twist

author: PatternRecognition

repo: OpenBMI


text:
function prctmus = gp_predprctmu(gp, x, y, varargin)
%GP_PREPRCTMU  Percentiles of the distribution of the location parameter
%
%  Description
%    PRCTMU = GP_PREDPRCTMU(GP, X, Y, XT, OPTIONS)
%    takes a GP structure together with matrix X of training inputs
%    and vector Y of training targets, and evaluates the
%    percentiles of the distribution of the location parameter at
%    test inputs XT.
%
%    PRCTMU = GP_PREDPRCTMU(GP, X, Y, OPTIONS)
%    evaluates the percentiles of the distribution of the location
%    parameter at training inputs X.
%
%    OPTIONS is optional parameter-value pair
%      prct   - percentiles to be computed (default = [5 50 95])
%      nsamp  - determines the number of samples used by GP_RND in case of
%               MCMC or IA (default = 5000).
%      predcf - index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn)
%      tstind - a vector defining, which rows of X belong to which
%               training block in *IC type sparse models. Default is [].
%               See also GP_PRED.
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.
%
%  See also
%    GP_PRED, GP_PAK, GP_UNPAK
%
% Copyright (c) 2011 Ville Tolvanen
% Copyright (c) 2011-2012 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'GP_PREDPRCTMU';
  ip.addRequired('gp',@(x) isstruct(x) || iscell(x));
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('prct', [5 50 95], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('nsamp', 5000, @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                   isvector(x) && isreal(x) && all(isfinite(x)&x>0))
  ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                   (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
  if numel(varargin)==0 || isnumeric(varargin{1})
    % inputParser should handle this, but it doesn't
    ip.parse(gp, x, y, varargin{:});
  else
    ip.parse(gp, x, y, [], varargin{:});
  end
  xt=ip.Results.xt;
  z = ip.Results.z;
  zt = ip.Results.zt;
  prct = ip.Results.prct;
  nsamp = ip.Results.nsamp;
  predcf=ip.Results.predcf;
  tstind=ip.Results.tstind;
  if isempty(xt)
    xt=x;
    if isempty(tstind)
      if iscell(gp)
        gptype=gp{1}.type;
      else
        gptype=gp.type;
      end
      switch gptype
        case {'FULL' 'VAR' 'DTC' 'SOR'}
          tstind = [];
        case {'FIC' 'CS+FIC'}
          tstind = 1:size(x,1);
        case 'PIC'
          if iscell(gp)
            tstind = gp{1}.tr_index;
          else
            tstind = gp.tr_index;
          end
      end
    end
    if isempty(zt)
      zt=z;
    end
  end

  % pass these forward
  options=struct();
  if ~isempty(z);options.z=z;end
  if ~isempty(zt);options.zt=zt;end
  if ~isempty(predcf);options.predcf=predcf;end
  if ~isempty(tstind);options.tstind=tstind;end

  [tn, nin] = size(x);

  if iscell(gp) || numel(gp.jitterSigma2)>1
    % gp_array or MCMC samples
    % the combined latent posterior is not Gaussian
    sampft = gp_rnd(gp, x, y, xt, 'nsamp', nsamp, options);
    if iscell(gp)
      % in the next step we need just one gp from the array
      gp=gp{1};
    end
    if isfield(gp, 'monotonic') && gp.monotonic
        [gp,x,y,z,xt,zt] = gp.fh.setUpDataForMonotonic(gp,x,y,z,xt,zt);
    end
    if isfield(gp.lik.fh,'trcov')
      % Gaussian likelihood
      prctmus = prctile(sampft', prct)';
    else
      prctmus = prctile(gp.lik.fh.invlink(gp.lik, sampft, zt)', prct)';
    end
  else
    % single GP
    if isfield(gp, 'monotonic') && gp.monotonic
        [gp,x,y,z,xt,zt] = gp.fh.setUpDataForMonotonic(gp,x,y,z,xt,zt);
    end
    % the latent posterior is Gaussian
    [Eft, Varft] = gp_pred(gp, x, y, xt, 'tstind', tstind, options);
    prct = prct./100;
    prct = norminv(prct, 0, 1);
    if isfield(gp.lik.fh,'trcov')
      % Gaussian likelihood
      prctmus = bsxfun(@plus, Eft, bsxfun(@times, sqrt(Varft), prct));
    else
      % Non-Gaussian likelihood
      np = length(prct);
      prctmus = zeros(size(Eft,1),np);
      for i=1:np
        prctmus(:,i) = gp.lik.fh.invlink(gp.lik, Eft+prct(i).*sqrt(Varft), zt);
      end
    end
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [Res] = rma2(P, Vec, Mat, Pth)

% rma2 replaces all descriptors from Mat except the ones in the given descriptor vector in the given position Pth and returns
%      and returns the vector with lower S_err, using multiple linear regression analysis.
%      (The initial descriptor is changed even if it has the lower S_err).
%
%	   Input:
%             P             Property vector
%             Vec          Descriptor vector
%             Mat           Descriptor matrix
%             Pth           Path to follow
%
%     Returns:
%
%             Res           Vector containing in the first place the
%                           S_err, and afterwards the corresponding descriptor vector
%
%
%
% Andrew G. Mercader
% INIFTA, La Plata, Argentina
% Created: 30 Jan 2007

if (nargin < 4)
   error('the function requires at least 4 input variables. Type ''help rma2''.');
end

[k, n_m] = size(Mat);

Num=[1:n_m];
Num(Vec)=[];

[k,n_n]=size(Num);

Smin=10000;   % A very big number compared to a normal S necessary just to start the program
for j=1:n_n;
Vec(Pth)=Num(j);
Ser=rms(P,Vec,Mat);
    if (Ser<Smin)
         Smin=Ser;
        desc=Num(j);
    end
end
Vec(Pth)=desc;
Res=[Smin, Vec];
%End of rma2

author: Sable

repo: mcbench-benchmarks


text:
classdef MimWSDataset < MimModel
    methods (Access = protected)
        function value = run(obj)
            seriesUid = obj.Parameters.seriesUid;
            value = obj.Callback.getMim().CreateDatasetFromUid(seriesUid);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function g = p01_g ( n, x )

%*****************************************************************************80
%
%% P01_G evaluates the gradient for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  th = p01_th ( x );

  r = sqrt ( x(1) * x(1) + x(2) * x(2) );
  t = x(3) - 10.0 * th;
  s1 = 5.0 * t / ( pi * r * r );

  g(1) = 200.0 * ( x(1) - x(1) / r + x(2) * s1 );
  g(2) = 200.0 * ( x(2) - x(2) / r - x(1) * s1 );
  g(3) = 2.0 * ( 100.0 * t + x(3) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
vectors_file_name = add_data_dir_base('GoogleNews_vectors_norm.mat');

word2vec_sqrt = false;
if word2vec_sqrt
	vectors_file_name = fname_concat(vectors_file_name, '_sqrt');
end

is_sampled = true;
% without normalizing we had numeric problems in the LMM
normalize  = true;

if is_sampled
    vectors_file_name = fname_concat(vectors_file_name, '_sampled');
end

output(1, 'reading vectors file %s\n', vectors_file_name);
load(vectors_file_name);
output(1, 'done\n');

TargetDim = 300;

%desiredTrnSizeForICA = 512000;

output(1, 'Use ICA to reduce to dimension %d\n', TargetDim);

[icasig, A, trans_matrix] = fastica (vectors_, 'numOfIC', TargetDim);
vectors_ = trans_matrix * vectors_;

[dim, num_vectors] = size(vectors_)

if normalize
    % normalize
    output(1, 'normalizing...\n');
    for i = 1:num_vectors
      vectors_(:,i) = norma(vectors_(:,i));
    end
end

output(1, 'saving to file...\n');

post_fix = sprintf('_ica_%d', TargetDim);
ica_vectors_file_name = fname_concat(vectors_file_name, post_fix);
trans_matrix_file_name = add_data_dir_base(strcat('trans_matrix', post_fix));

save(ica_vectors_file_name, 'vectors_', '-v7.3');
output(1, 'saved to file: %s\n', ica_vectors_file_name);

save(trans_matrix_file_name, 'trans_matrix', '-v7.3');
output(1, 'saved to file: %s\n', trans_matrix_file_name);

clear vectors_;
clear trans_matrix;

author: layumi

repo: Image-Text-Embedding


text:
function [ p, t ] = cvt_square_uniform ( n, sample_num, delaunay_display )

%*****************************************************************************80
%
%% CVT_SQUARE_UNIFORM demonstrates how a CVT can be computed and displayed in MATLAB.
%
%  Discussion:
%
%    This simple example carries out an iterative CVT calculation in the
%    unit square, with a uniform density.  The initial placement of the
%    generators is random.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review,
%    Volume 41, 1999, pages 637-676.
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, integer SAMPLE_NUM, the number of sample points.
%
%    Input, logical DELAUNAY_DISPLAY, is TRUE (nonzero) if the Delaunay
%    triangulation is to be displayed.
%
%    Output, real P(N,2), the location of the generators.
%
%    Output, integer T(NT,3), information defining the Delaunay
%    triangulation of the generators.  NT is the number of triangles,
%    which varies depending on the arrangement of the generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_SQUARE_UNIFORM:\n' );
  fprintf ( 1, '  A simple demonstration of a CVT computation\n' );
  fprintf ( 1, '  (Centroidal Voronoi Tessellation)\n' );
  fprintf ( 1, '  in a square, with a uniform density.\n' );

  if ( nargin < 1 )
    n = 100;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_SQUARE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of N was supplied.\n' );
    fprintf ( 1, '  N is the number of generators.\n' );
    fprintf ( 1, '  A default value N = %d will be used.\n', n );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of generators = %d\n',  n );
  end

  if ( nargin < 2 )
    sample_num = 1000 * n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_SQUARE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of SAMPLE_NUM was supplied.\n' );
    fprintf ( 1, '  SAMPLE_NUM is the number of sample points.\n' );
    fprintf ( 1, '  A default value SAMPLE_NUM = %d will be used.\n', ...
      sample_num );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified number of sample points = %d\n', ...
      sample_num );
  end

  if ( nargin < 3 )
    delaunay_display = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_SQUARE_UNIFORM - Note:\n' );
    fprintf ( 1, '  No value of DELAUNAY_DISPLAY was supplied.\n' );
    fprintf ( 1, '  DELAUNAY_DISPLAY is TRUE (nonzero) if the\n' );
    fprintf ( 1, '  Delaunay triangulation is also to be displayed.\n' );
    fprintf ( 1, '  A default value DELAUNAY_DISPLAY = %d will be used.\n', ...
      delaunay_display );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  User specified DELAUNAY_DISPLAY = %d\n', ...
      delaunay_display );
  end
%
%  This switch is set to 1 (TRUE) if the ACCUMARRAY command is available.
%  That speeds up the calculation a lot.  If you don't have the ACCUMARRAY
%  command, just set this to 0.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_SQUARE_UNIFORM:\n' );
  fprintf ( 1, '  MATLAB''s ACCUMARRAY command can be used for faster\n' );
  fprintf ( 1, '  computation.  This command is not available in\n' );
  fprintf ( 1, '  some versions of MATLAB.  If ACCUMARRAY is available,\n' );
  fprintf ( 1, '  simply make sure that the ACCUMARAY_AVAILABLE variable\n' );
  fprintf ( 1, '  is set to 1!\n' );

  accumarray_available = 1;

  if ( accumarray_available )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The ACCUMARRAY command will be used.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The ACCUMARRAY command will NOT be used.\n' );
  end
%
%  Clear the figure screen, if already open.
%
  clf
%
%  Randomize the initial locations of the generators in the unit square.
%  If another region is used, then this initialization should be changed.
%
  p = rand ( n, 2 );
  plot ( p(:,1), p(:,2), 'b.' );
  axis ( [ 0.0, 1.0, 0.0, 1.0 ] )
  line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
    'Color', 'r', 'LineWidth', 2 );
  title_string = sprintf ( 'Initial Generators' );
  title ( title_string );
  axis equal
  drawnow

  string = input ( 'RETURN, or Q to quit: ', 's' );

  if ( string == 'q' | string == 'Q' )
    return
  end

  it = 0;

  while ( 1 )
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = delaunay ( p(:,1), p(:,2) );
%
%  Display the Delaunay triangulation, if requested.
%
    if ( delaunay_display )
      subplot ( 1, 2, 2 )
      trimesh ( t, p(:,1), p(:,2), zeros(n,1), 'EdgeColor', 'b' )
      hold on
      plot ( p(:,1), p(:,2), 'b.' );
      axis ( [ -0.0, 1.0, -0.0, 1.0 ] )
      line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
        'Color', 'r', 'LineWidth', 2 );
      title_string = sprintf ( 'Delaunay, step %d', it );
      title ( title_string );
      axis equal
      view ( 2 )
      hold off
    end
%
%  Display the CVT generators, and the associated Voronoi diagram.
%
    if ( delaunay_display )
      subplot ( 1, 2, 1 )
    end

    voronoi ( p(:,1), p(:,2), t );
    axis ( [ -0.0, 1.0, -0.0, 1.0 ] )
    line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ], ...
      'Color', 'r', 'LineWidth', 2 );
    title_string = sprintf ( 'Voronoi, step %d', it );
    title ( title_string );
    axis equal
    axis tight
    drawnow
%
%  Generate sample points.
%
%  These sample points implicitly define the geometry of the region.
%  If the region is not a unit square, then the range of the sample
%  data must be changed.
%
%  The data is sampled uniformly.  If a nonuniform density is desired,
%  then the sampling must be done in a biased way.
%
    xs = rand ( sample_num, 1 );
    ys = rand ( sample_num, 1 );
%
%  For each sample point, find K, the index of the nearest generator.
%  We do this efficiently by using the Delaunay information with
%  Matlab's DSEARCH command, rather than a brute force nearest neighbor
%  computation.
%
    k(1:sample_num,1) = dsearch ( p(:,1), p(:,2), t, xs, ys );
%
%  The centroid of the Voronoi region associated with each generator
%  is approximated by the average of the sample points it was closest to.
%
    if ( accumarray_available )

      count(1:n) = accumarray ( k, ones(sample_num,1) );
      centroid(1,1:n) = accumarray ( k, xs );
      centroid(2,1:n) = accumarray ( k, ys );

    else

      count(1:n) = 0;
      centroid(1,1:n) = 0.0;
      centroid(2,1:n) = 0.0;

      for i = 1 : sample_num
        j = k(i);
        count(j) = count(j) + 1;
        centroid(1,j) = centroid(1,j) + xs(i);
        centroid(2,j) = centroid(2,j) + ys(i);
      end

    end
%
%  Replace the generators by the centroids.
%
    p(1:n,1) = ( centroid(1,1:n) ./ count(1:n) )';
    p(1:n,2) = ( centroid(2,1:n) ./ count(1:n) )';

    string = input ( 'RETURN, or Q to quit: ', 's' );

    if ( string == 'q' | string == 'Q' )
      break
    end

    it = it + 1;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_SQUARE_UNIFORM:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function p = square_uniform ( n )

%*****************************************************************************80
%
%% SQUARE_UNIFORM returns sample points from the unit square.
%
%  Discussion:
%
%    This routine returns N points sampled uniformly at random
%    from within the unit square.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points to generate.
%
%    Output, real P(N,2), the sample points.
%
  p(1:n,1:2) = rand(n,2);

  return
end


text:
% Test trajectory primitives

mdl_puma560

[q,qd,qdd] = jtraj(qz, qr, 20)
q
qd
qdd

[q,qd,qdd] = jtraj(qz, qr, 20, 0.1*mat(ones(1,6)), -0.1*mat(ones(1,6)) )
q
qd
qdd

[q,qd,qdd] = jtraj(qz, qr, [0:0.2:10])
q

t1 = trotx(0.1) * transl(0.2, 0.3, 0.4)
t1
t2 = troty(-0.3) * transl(-0.2, -0.3, 0.6)
t2
ctraj(t1, t2, 5)
ctraj(t1, t2, [0:0.1:1])

author: petercorke

repo: robotics-toolbox-matlab


text:
classdef nnspnorm < nntest
  methods (Test)
   function basic(test)
      h = 13 ;
      w = 17 ;
      d = 4 ;
      n = 5 ;
      param = [3, 3, 0.1, 0.75] ;
      x = test.randn(h,w,d,n,'single') ;
      y = vl_nnspnorm(x, param) ;
      dzdy = test.rand(h, w, d, n) ;
      dzdx = vl_nnspnorm(x, param, dzdy) ;
      test.der(@(x) vl_nnspnorm(x,param), x, dzdy, dzdx, test.range * 1e-3) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [ chg mag geometry ] = import_chgcar( filename )
%IMPORT_CHGCAR Import a VASP CHGCAR file.
%   [chg,mag,geometry] = import_chgcar(filename)
%   Import a VASP CHGCAR file. If no filename is specified, data is read
%   from CHGCAR. chg and mag are three dimensional arrays containing the
%   charge magnetization densities in Bohr magneton per cubic Angstrom.
%   Note that these are not the same units as the CHGCAR file. geometry is
%   a struct describing the cell geometry; see IMPORT_POSCAR for a detailed
%   description.
%
%   See also IMPORT_POSCAR, IMPORT_LOCPOT.

% todo:
% check compatibility with non-spin-polarized files
% extract chemical symbols
% what about AECAR and ELFCAR files?

  if nargin == 0
      filename='CHGCAR';
  end

  fid = fopen(filename);
  if fid==-1
    error(['File ' filename ' not found']);
  end

    geometry = import_poscar(fid);

    vol = abs(dot(geometry.lattice(1,:),cross(geometry.lattice(2,:),geometry.lattice(3,:))));
    natoms = sum(geometry.atomcount);

    fgetl(fid); % blank line

    gridsize = fscanf(fid, '%d %d %d', [3 1])';

    chg = fscanf(fid, '%f', [prod(gridsize,2) 1])';
    chg = reshape(chg,gridsize);
    chg = chg/vol;

    fgetl(fid); % empty string (or padding)

    pos = ftell(fid);
    line = fgetl(fid);
    fseek(fid,pos,'bof');
    if line(1)=='a'
        for i = 1:natoms
            line = fgetl(fid);
            nentries = sscanf(line,['augmentation occupancies ' num2str(i)...
                ' %d']); % number of occupancy entries
            fscanf(fid, '%f', [nentries 1]);
            fgetl(fid); % empty string
        end
        fscanf(fid, '%f', [natoms 1]); % don't know what these are
        fgetl(fid); % empty string
    end

    pos = ftell(fid);
    fgetl(fid);
    if ~feof(fid)
        fseek(fid,pos,'bof');
        gridsize = fscanf(fid, '%d %d %d', [3 1])';

        mag = fscanf(fid, '%f', [prod(gridsize,2) 1])';
        mag = reshape(mag,gridsize);
        mag = mag/vol;
    else
        mag = [];
    end

    fclose(fid);

end


text:
clear('pi');
y           = load('data/M3_21.dat')';
airline     = ssm_airline;
airline     = estimate(y, airline, 0.1);
theta       = -airline.param(1);
Theta12     = (-airline.param(2))^(1/12);
zetavar     = airline.param(3);

param0  = [theta+Theta12 -theta*Theta12 Theta12 Theta12 zetavar];
result4 = zeros(31, 6);
n       = 1;
% 4-5-1
for i = 1 : 6
    genair          = ssm_genair(4, 5, i);
    [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
    result4(n, :)   = [roots([1 -genair.param(1:2)])' genair.param(3:4) reallog(genair.param(5))/2 -2*logL-13*reallog(2*pi)+2*4];
    n               = n + 1;
end
% 4-4-2
for i = 1 : 5
    for j = i+1 : 6
        genair          = ssm_genair(4, 4, [i j]);
        [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
        result4(n, :)   = [roots([1 -genair.param(1:2)])' genair.param(3:4) reallog(genair.param(5))/2 -2*logL-13*reallog(2*pi)+2*4];
        n               = n + 1;
    end
end
% 4-3-3
for i = 2 : 5
    for j = i+1 : 6
        genair          = ssm_genair(4, 3, [1 i j]);
        [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
        result4(n, :)   = [roots([1 -genair.param(1:2)])' genair.param(3:4) reallog(genair.param(5))/2 -2*logL-13*reallog(2*pi)+2*4];
        n               = n + 1;
    end
end
fline   = '%.5f\t\t%.5f\t\t%.5f\t\t%.5f\t\t%.5f\t\t%.5f\n';
fprintf(1, fline, result4');

param0  = [theta Theta12 Theta12 zetavar];
result3 = zeros(41, 5);
n       = 1;
% 3-5-1
for i = 1 : 6
    genair          = ssm_genair(3, 5, i);
    [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
    result3(n, :)   = [genair.param(1:3) reallog(genair.param(4))/2 -2*logL-13*reallog(2*pi)+2*3];
    n               = n + 1;
end
% 3-4-2
for i = 1 : 5
    for j = i+1 : 6
        genair          = ssm_genair(3, 4, [i j]);
        [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
        result3(n, :)   = [genair.param(1:3) reallog(genair.param(4))/2 -2*logL-13*reallog(2*pi)+2*3];
        n               = n + 1;
    end
end
% 3-3-3
for i = 1 : 4
    for j = i+1 : 5
        for k = j+1 : 6
            genair          = ssm_genair(3, 3, [i j k]);
            [genair logL]   = estimate(y, genair, param0, [], 'fmin', 'simplex', 'disp', 'off');
            result3(n, :)   = [genair.param(1:3) reallog(genair.param(4))/2 -2*logL-13*reallog(2*pi)+2*3];
            n               = n + 1;
        end
    end
end
fline   = '%.5f\t\t%.5f\t\t%.5f\t\t%.5f\t\t%.5f\n';
fprintf(1, fline, result3');

repo: REST


text:
function a = biw ( n )

%*****************************************************************************80
%
%% BIW returns the BIW matrix.
%
%  Discussion:
%
%    BIW is a bidiagonal matrix of Wilkinson.   Originally, this matrix
%    was considered for N = 100.
%
%  Formula:
%
%    if ( I == J )
%      A(I,J) = 0.5 + I / ( 10 * N )
%    else if ( J == I+1 )
%      A(I,J) = -1.0
%    else
%      A(I,J) = 0
%
%  Example:
%
%    N = 5
%
%    0.52 -1.00  0.00  0.00  0.00
%    0.00  0.54 -1.00  0.00  0.00
%    0.00  0.00  0.56 -1.00  0.00
%    0.00  0.00  0.00  0.58 -1.00
%    0.00  0.00  0.00  0.00  0.60
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    a(i,i) = 0.5 + i / ( 10 * n );
  end

  for i = 1 : n - 1
    a(i,i+1) = - 1.0;
  end

  return
end


text:
function X = QPbyUseSFC(waypts, ts, decomp)

%% condition
traj.n_order = 7;
traj.n_poly = size(waypts, 1) - 1;
traj.p0 = waypts(1,:);
traj.pe = waypts(end,:);
traj.v0 = [0,0,0];
traj.ve = [0,0,0];
traj.a0 = [0,0,0];
traj.ae = [0,0,0];

%% trajectory plan
[minSnapValue, px, py, pz] = minimum_snap_three_axis_SFC(traj, ts, decomp);

disp(['minSnapValue is : ',num2str(minSnapValue)]);

X = [px py pz];
end

% v0 = [1*3]
function [minValue, px, py, pz] = minimum_snap_three_axis_SFC(traj, ts, decomp)
    % Dimension
    dim = 3;
    n_order = traj.n_order;
	p0 = traj.p0;
	pe = traj.pe;
    v0 = traj.v0;
    ve = traj.ve;
    a0 = traj.a0;
    ae = traj.ae;

	n_poly = traj.n_poly;
	n_coef = n_order+1;

	%% compute Q,  Q is a Symmetric matrix
    % Q_1 =   [0    0   0]
    %         [0    f   f]
    %         [0    f   f]          % f =  integral of ts(i) ~ ts(i+1)
    % Q_x =   [Q_1  0     0    0]
    %         [0    Q_2   0    0]
    %         [0    0   Q_...  0]
    %         [0    0     0  Q_n]   % n = n_poly
    % Q_all = [Q_x  0     0]
    %         [0    Q_y   0]
    %         [0    0   Q_z]        % Q_x == Q_y == Q_z
    Q_x = [];
	for i=1:n_poly
	    Q_x = blkdiag(Q_x,calc_Q(n_order,4,ts(i),ts(i+1)));
    end
    zeroM = zeros(size(Q_x));
    Q_all = [Q_x zeroM zeroM; zeroM Q_x zeroM; zeroM zeroM Q_x];
    xlen = size(Q_all,1);
	f = zeros(xlen,1);      %% min (1/2p^TQp + f^Tp)

    %% compute Aeq x = beq
    % beacuse Aeq_x == Aeq_y == Aeq_z
    % Aeq = [Aeq_x  0      0]
    %       [0    Aeq_y    0]
    %       [0      0  Aeq_z]
	neq = 8;  %% (8 equations)
    eqNum_x = (7*(n_poly-1)+neq);
	Aeq_x = zeros(eqNum_x, n_coef*n_poly);
	beq = zeros(dim*eqNum_x, 1);

	% start/terminal pva constraints  (8 equations)
	Aeq_x(1:4,1:n_coef) = [calc_dc(ts(1),n_order,0);
	                     calc_dc(ts(1),n_order,1);
	                     calc_dc(ts(1),n_order,2);
	                     calc_dc(ts(1),n_order,3)];
	Aeq_x(5:8,n_coef*(n_poly-1)+1:n_coef*n_poly) = ...
	                    [calc_dc(ts(end),n_order,0);
	                     calc_dc(ts(end),n_order,1);
	                     calc_dc(ts(end),n_order,2);
	                     calc_dc(ts(end),n_order,3)];
	beq(1:8,1)                              = [p0(1),v0(1),a0(1),0,pe(1),ve(1),ae(1),0]';
    beq((eqNum_x   + 1):(eqNum_x   + 8),1)  = [p0(2),v0(2),a0(2),0,pe(2),ve(2),ae(2),0]';
    beq((eqNum_x*2 + 1):(eqNum_x*2 + 8),1)  = [p0(3),v0(3),a0(3),0,pe(3),ve(3),ae(3),0]';

	% continuous constraints  ((n_poly-1)*7 equations)
	for i=1:n_poly-1
		t_derc_p = calc_dc(ts(i+1),n_order,0);
        t_derc_v = calc_dc(ts(i+1),n_order,1);
        t_derc_a = calc_dc(ts(i+1),n_order,2);
        t_derc_j = calc_dc(ts(i+1),n_order,3);  % jerk
        t_derc_4 = calc_dc(ts(i+1),n_order,4);
        t_derc_5 = calc_dc(ts(i+1),n_order,5);
        t_derc_6 = calc_dc(ts(i+1),n_order,6);
        Aeq_x(neq+1,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_p,-t_derc_p];
        Aeq_x(neq+2,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_v,-t_derc_v];
        Aeq_x(neq+3,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_a,-t_derc_a];
        Aeq_x(neq+4,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_j,-t_derc_j];
        Aeq_x(neq+5,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_4,-t_derc_4];
        Aeq_x(neq+6,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_5,-t_derc_5];
        Aeq_x(neq+7,n_coef*(i-1)+1:n_coef*(i+1))=[t_derc_6,-t_derc_6];
        neq = neq + 7;
    end

    zeroM = zeros(size(Aeq_x));
    Aeq = [Aeq_x zeroM zeroM; zeroM Aeq_x zeroM; zeroM zeroM Aeq_x];

    %% compute Ax <= b
    A = [];
    b = [];
    ieq = 1;
    xll = xlen / dim;  % 1/3 of xlen,  = n_coef * (n_poly-1)
    for i = 1 : n_poly
        planesCell = decomp.lines_{i}.polyhedron_.polys_;
        [~, lenj] = size(planesCell);
        for j = 1 : lenj
            n = planesCell{j}.n_;
            p = planesCell{j}.p_;
            nx = zeros(1, xlen);
            % add (Denominator + 1) points on the poly's Ax<b
            Denominator = 2;  % Denominator > 0
            for k = 0 : Denominator
                tvec_mid = calc_dc(ts(i) +  k*(ts(i+1) - ts(i))/Denominator, n_order, 0);   % [1*8]vector
                nx(1, 8*(i-1) + 1 : 8*(i-1) + 8)                = tvec_mid.*n(1);
                nx(1, 8*(i-1) + 1 + xll : 8*(i-1) + 8 + xll)    = tvec_mid.*n(2);
                nx(1, 8*(i-1) + 1 + 2*xll : 8*(i-1) + 8 + 2*xll)= tvec_mid.*n(3);
                A(ieq,:) = nx;
                b(ieq,:) = dot(n, p);
                % dot(plane.n, p1 - plane.p) < 0
                % => dot(plane.n, p1) - dot(plane.n, plane.p) < 0
                % => dot(plane.n, p1) < dot(plane.n, plane.p)
                %% Ax < b
                ieq = ieq + 1;
            end
        end
    end

%     % Modify the number of iterations
    options = optimoptions('fmincon');
    options.MaxIterations = 2000;   % Defaults = 1000
    lb = []; ub = []; x0 = [];
    p = quadprog(Q_all,f,A,b,Aeq,beq, lb,ub,x0, options);

% 	p = quadprog(Q_all,f,A,b,Aeq,beq);

    minValue = p'*Q_all*p;

    px = p(1:xll,1);
    py = p(xll + 1:2*xll,1);
    pz = p(xll*2 + 1:3*xll,1);
end

author: LenaShengzhen

repo: AerialRobotics


text:
%% Copyright (C) 2015 Carnë Draug
%% Copyright (C) 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deffn Constant eulergamma ()
%% Return Euler--Mascheroni constant.
%%
%% @example
%% @group
%% eulergamma
%%   @result{} (sym) γ
%%
%% vpa (eulergamma ())
%%   @result{} (sym) 0.57721566490153286060651209008240
%% @end group
%% @end example
%%
%% @seealso{catalan}
%% @end deffn

%% Author: Carnë Draug
%% Keywords: symbolic, constants

function g = eulergamma ()

  if (nargin ~= 0)
    print_usage ();
  end

  g = pycall_sympy__ ('return sympy.S.EulerGamma,');
end

%!error catalan (sym(1))
%!assert (double (eulergamma ()) > 0.577215664901)
%!assert (double (eulergamma ()) < 0.577215664902)

repo: octsympy


text:
function [X, present] = block(X, past, weights, hyperParameters)
% block   Transformer block for GPT-2
%
%   [X, present] = block(X, past, weights, hyperParameters) computes a
%   GPT-2 style transformer block on the input X as described in [1] (see
%   Section 2.3). One difference between this style of transformer block
%   and others is that this block uses layer normalization at the
%   beginning.
%
%   Inputs:
%       X               - A (numFeatures*numHeads)-by-numInputSubwords
%                         input array.
%       past            - A numFeatures-by-numPastSubwords-by-numHeads-by-2
%                         array. This contains the 'keys' and 'values' for
%                         past subwords. These are needed to predict future
%                         outputs in an autoregressive manner. 'keys' are
%                         stored in past(:,:,:,1) and 'values' are stored
%                         in past(:,:,:,2).
%       weights         - The weights for the transformer block stored in a
%                         struct. In this block we have:
%                           - ln_1_g_0: Weight vector for the first layer
%                             normalization.
%                           - ln_1_b_0: Bias vector for the first layer
%                             normalization.
%                           - ln_2_g_0: Weight vector for the second layer
%                             normalization.
%                           - ln_2_b_0: Bias vector for the second layer
%                             normalization.
%                         In the attention sub-block:
%                           - attn_c_attn_w_0: A weight matrix for the
%                             first fully connected layer.
%                           - attn_c_attn_b_0: A bias vector for the first
%                             fully connected layer.
%                           - attn_c_proj_w_0: A weight matrix for the
%                             final fully connected layer.
%                           - attn_c_proj_b_0: A bias vector for the final
%                             fully connected layer.
%                         In the multi-layer perceptron block:
%                           - mlp_c_fc_w_0: A weight matrix for the first
%                             fully connected layer.
%                           - mlp_c_fc_b_0: A bias vector for the first
%                             fully connected layer.
%                           - mlp_c_proj_w_0: A weight matrix for the
%                             second fully connected layer.
%                           - mlp_c_proj_b_0: A bias vector for the second
%                             fully connected layer.
%       numHeads        - The number of attention heads. This is a
%                         hyper-parameter.
%
%   Outputs:
%       Z               - A (numFeatures*numHeads)-by-numInputSubwords
%                         output array.
%       present         - A numFeatures-by-numAllSubwords-by-numHeads-by-2
%                         array. This contains the 'keys' and 'values' that
%                         are created from inputs. These need to passed
%                         back in as the 'past' input if we want to predict
%                         future outputs in an autoregressive manner. 'keys'
%                         are stored in present(:,:,:,1) and 'values' are
%                         stored in present(:,:,:,2).
%
%   References:
%
%   [1] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei,
%       Ilya Sutskever, "Language Models are Unsupervised Multitask
%       Learners",
%       https://cdn.openai.com/better-language-models/language_models_are_unsupervised_multitask_learners.pdf

XNorm1 = transformer.layer.normalization(X, ...
    weights.ln_1_g_0, weights.ln_1_b_0);

[A, present] = transformer.layer.attention(XNorm1, past, weights, hyperParameters);

X = X + A;

XNorm2 = transformer.layer.normalization(X, ...
    weights.ln_2_g_0, weights.ln_2_b_0);

M = transformer.layer.multiLayerPerceptron(XNorm2, weights);

X = X + M;

end

author: matlab-deep-learning

repo: transformer-models


text:
clear; close all; clc;

confmat = magic(3); % sample data
% plotting
plotConfMat(confmat, {'Dog', 'Cat', 'Horse'});

author: opencobra

repo: cobratoolbox


text:
function tapas_sem_display_posterior_summary(results, summary)
%% Display the posterior summary.
%
% Input
%       results     -- Results structure
%       summary     -- Summary structure of the posterior.
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2019
%

CONG = 0;
INCONG = 1;

data = results.data;

% Offset of the data
if ~isfield(data(1).y, 'offset')
    offset = 0;
else
    offset = data(1).y.offset;
end

% Assume data in milliseconds and a scaling of 100ms
if ~isfield(data(1).y, 'scale')
    scale = 100;
else
    scale = data(1).y.scale;
end

% Collapse all the data
t = arrayfun(@(x) x.y.t, data, 'UniformOutput', false);
t = vertcat(t{:});

t = scale * t + offset;

% (a)ctions
a = arrayfun(@(x) x.y.a, data, 'UniformOutput', false);
a = vertcat(a{:});

% (t)rial (t)ype
tt = arrayfun(@(x) x.u.tt, data, 'UniformOutput', false);
tt = vertcat(tt{:});

% Plot only correct trials
ttype = mod(tt, 2);
correct = a == ttype;

%% Normalized fits

fits = summary.fits;
nf = numel(fits);

% Normalized fits
nfits = cell(nf, 1);

for i = 1:numel(fits)
    nfits{i} = tapas_sem_normalized_fits(data(i), fits{i});
end

%% Display a table with the summaries
tapas_display_table(summary.summaries, 'Posterior summary')

%% Plot all the responses combined
fig = figure('name', 'Group responses and fits');
fig.Color = 'w';

all_y = struct('t', t, 'a', a);
all_u = struct('tt', tt);
edges = tapas_sem_plot_responses(all_y, all_u);
dt = edges(2) - edges(1);

plot_predicted_normalized_fits(nfits, dt, offset, scale);

%% Delta plots
fig = figure('name', 'Group delta plots');
fig.Color = 'w';

hold on
% Plot all together
tapas_sem_empirical_delta_plot(...
    t(correct & (ttype == CONG)), ...
    t(correct & (ttype == INCONG)));

xlabel('time')
ylabel('Delta incong. RT - cong. RT')
% Make delta plots
time = scale * fits{1}(1).t + offset;

nt = numel(time);

cong = zeros(nt, 1);
incong = zeros(nt, 1);

for i = 1:numel(nfits)
    for j = 1:numel(nfits{i})
        cong = cong + nfits{i}(j).pro;
        incong = incong + nfits{i}(j).anti;
    end
end

time = reshape(time, numel(time), 1);

tapas_sem_predicted_delta_plot(time, cong, incong);
legend({'Empirical delta', 'Predicted delta'})
end

function plot_predicted_normalized_fits(nfits, dt, offset, scale)

fits = nfits{1};

for i = 2:numel(nfits)
    for j = 1:numel(fits)
        fits(j).pro = fits(j).pro + nfits{i}(j).pro;
        fits(j).anti = fits(j).anti + nfits{i}(j).anti;
    end
end

for j = 1:numel(fits)
    fits(j).t = fits(j).t * scale + offset;
end

tapas_sem_plot_fits(fits, dt/scale)

end

author: translationalneuromodeling

repo: tapas


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcSharpen(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

choice = questdlg('Apply to all slices?','All Slices','No','Yes','No');
if isempty(choice)
    warndlg('Image sharpening is cancelled.');
    return;
end
% Handle response
H =  fspecial('unsharp');
switch choice
    case 'No'
        handles.TMatrix(:,:,handles.V.Slice) = imfilter(handles.BMatrix,H,'replicate');
    case 'Yes'
        if length(handles.V.DimSize)>2
            for i= 1: handles.V.DimSize(3)
                handles.TMatrix(:,:,i) = imfilter(handles.TMatrix(:,:,i),H,'replicate');
                MU_update_waitbar(handles.Progress_axes,i,handles.V.DimSize(3));
            end
        else
            handles.TMatrix = imfilter(handles.BMatrix,H,'replicate');
        end
end

MergeM=get(handles.Matrix_name_edit,'String');
set(handles.Matrix_name_edit,'String',[MergeM '_shp']);

% update current display matrix
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function CPD = set_params(CPD, varargin)
% SET_PARAMS Set the parameters (fields) for a softmax_CPD object
% CPD = set_params(CPD, name/value pairs)
%
% The following optional arguments can be specified in the form of name/value pairs:
% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y), Q1=ns(dps(1)), Q2=ns(dps(2)), ...
%   where dps are the discrete parents; if there are no discrete parents, we set Q1=1.)
%
% weights - (W(:,j,a,b,...) - W(:,j',a,b,...)) is ppn to dec. boundary
%           between j,j' given Q1=a,Q2=b,... [ randn(X,Y,Q1,Q2,...) ]
% offset  - (offset(j,a,b,...) - offset(j',a,b,...)) is the offset to dec. boundary
%           between j,j' given Q1=a,Q2=b,... [ randn(Y,Q1,Q2,...) ]
% clamped     - 'yes' means don't adjust params during learning ['no']
% max_iter    - the maximum number of steps to take [10]
% verbose     - 'yes' means print the LL at each step of IRLS ['no']
% wthresh     - convergence threshold for weights [1e-2]
% llthresh    - convergence threshold for log likelihood [1e-2]
% approx_hess - 'yes' means approximate the Hessian for speed ['no']
%
% e.g., CPD = set_params(CPD,'offset', zeros(ns(i),1));

args = varargin;
nargs = length(args);
glimsz = prod(CPD.sizes(CPD.dpndx));
for i=1:2:nargs
  switch args{i},
   case 'discrete',     str='nothing to do';
   case 'clamped',      CPD = set_clamped(CPD, strcmp(args{i+1}, 'yes'));
   case 'max_iter',     CPD.max_iter = args{i+1};
   case 'verbose',      CPD.verbose = strcmp(args{i+1}, 'yes');
   case 'max_iter',     CPD.max_iter = args{i+1};
   case 'wthresh',      CPD.wthresh = args{i+1};
   case 'llthresh',     CPD.llthresh = args{i+1};
   case 'approx_hess',  CPD.approx_hess = strcmp(args{i+1}, 'yes');
   case 'weights',      for q=1:glimsz, CPD.glim{q}.w1 = args{i+1}(:,:,q); end;
   case 'offset',
    if glimsz == 1
      CPD.glim{1}.b1 = args{i+1};
    else
      for q=1:glimsz, CPD.glim{q}.b1 = args{i+1}(:,q); end;
    end
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [A]=rician_LMMSE_filter(M,S,k)

%Based on Aja-Fernandez et al. 2008. Restoration of DWI Data Using a Rician
%LMMSE Estimator

%%
% Sg=1;
pad_step=(k-1)/2;
% hg=gaussiankernel(2,k,Sg); %Equivalent to: hg = fspecial('gaussian',[k k], S);

% hg=gauss_kernel(k,nd,f,m);
% hg=gauss_kernel(k,2,Sg,'sigma');
hg=gauss_kernel(k,2,1.5,'width');

Mij2=M.^2; Mij4=M.^4;

%% Gaussian filtering M.^2

[Mij2_gauss]=padrep(Mij2,pad_step.*ones(1,3)); %Padding array
Mij2_gauss=convn(Mij2_gauss,hg,'same');
Mij2_gauss=Mij2_gauss(pad_step+1:(size(Mij2_gauss,1)-pad_step),pad_step+1:(size(Mij2_gauss,2)-pad_step),pad_step+1:(size(Mij2_gauss,3)-pad_step)); %Trimming back to normal size

%% Gaussian filtering M.^4

[Mij4_gauss]=padrep(Mij4,pad_step.*ones(1,3)); %Padding array
Mij4_gauss=convn(Mij4_gauss,hg,'same');
Mij4_gauss=Mij4_gauss(pad_step+1:(size(Mij4_gauss,1)-pad_step),pad_step+1:(size(Mij4_gauss,2)-pad_step),pad_step+1:(size(Mij4_gauss,3)-pad_step)); %Trimming back to normal size

%% Calculating K

Kij_nom=(4.*(S.^2)).*(Mij2_gauss-(S.^2));
Kij_denom=Mij4_gauss-(Mij2_gauss.^2);
L=Kij_denom~=0; %Logic index of non zero values
Kij=zeros(size(Kij_denom));
Kij(L)=1-(Kij_nom(L)./Kij_denom(L));
Kij(Kij<0)=0;

%% Calculating signal A

A2= Mij2_gauss-2.*(S.^2) + Kij.*(Mij2-Mij2_gauss) ;
A2(A2<0)=0; %Smaller then zero set to zero
A=sqrt(A2);

%% END

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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 3 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [V,F,TV,TF,obj] = unbake_normal_map(V,F,TV,TF,nim,varargin)
  % UNBAKE_NORMAL_MAP Given a uv-mapped (TV,TF) 3D surface (V,F) mesh and a
  % corresponding normal map image (nim), recover a high-resolution 3D geometry
  % that "explains" the normal map. Resolution (and runtime) is controlled by
  % the resolution of the normal map image (hint: start with a downsampled
  % image).
  %
  % [V,F,TV,TF,obj] = unbake_normal_map(V,F,TV,TF,nim,varargin)
  %
  % Inputs:
  %   V  #V by 3 list of 3D vertex positions
  %   F  #F by 3 list of triangle indices into V
  %   TV  #TV by 3 list of 2D uv vertex positions
  %   TF  #F by 3 list of triangle indices into TV
  %   nim  #nim by #nim by 3 image of *object-space* normals (in range [0,1])
  % Outupts:
  %   V  #V by 3 list of 3D vertex positions
  %   F  #F by 3 list of triangle indices into V
  %   TV  #TV by 3 list of 2D uv vertex positions
  %   TF  #F by 3 list of triangle indices into TV
  %
  % See also: normal_to_displacement_map, tangent_to_object_normal_map
  %

  function [sel,V,F,TV,TF] = upsample_helper(V,F,TV,TF)
    mean_area = mean(doublearea(TV,TF));
    pixel_area = 1./(size(nim,1)*size(nim,2));
    dblA = doublearea(TV,TF);
    sel = find(dblA>pixel_area);
    if isempty(sel)
      return;
    end
    % Selected edges on 3D mesh (V,F)
    [uE,~,EMAP] = unique(sort(reshape(F(:,[2 3 1 3 1 2]),[],2),2),'rows');
    uEM = sparse(EMAP,1,repmat(sparse(sel,1,1,size(F,1),1),3,1),size(uE,1),1)>0;
    MF = full(reshape(uEM(EMAP),[],3));
    % Selected edges on texture mesh (TV,TF)
    [uE,~,EMAP] = unique(sort(reshape(TF(:,[2 3 1 3 1 2]),[],2),2),'rows');
    uEM = sparse(EMAP,1,repmat(sparse(sel,1,1,size(TF,1),1),3,1),size(uE,1),1)>0;
    MTF = full(reshape(uEM(EMAP),[],3));
    M = MF | MTF;
    [TV,TF] = upsample(TV,TF,'OnlySelected',M);
    mprev = size(F,1);
    [V,F] = upsample(V,F,'OnlySelected',M);
    if ~quiet
      fprintf('#F: %g -> %g\n',mprev,size(F,1));
    end
  end

  max_iters = 1000;
  quiet = true;
  quiet = false;
  max_step_size_seen = -inf;
  method = 'lbfgs';
  progressive_upsampling = true;
  fix_symmetry = false;

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'FixSymmetry','Method','Quiet','Progressive'}, ...
    {'fix_symmetry','method','quiet','progressive_upsampling'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % Progressive Upsampling is better: similar amount of time for a far better
  % solution.
  if ~progressive_upsampling
    while true
      [sel,V,F,TV,TF] = upsample_helper(V,F,TV,TF);
      if isempty(sel)
        break;
      end
    end
  end

  cross2 = @(a,b,c) ...
    [a(:,2).*b(:,3)-a(:,3).*b(:,2), ...
     a(:,3).*b(:,1)-a(:,1).*b(:,3), ...
     a(:,1).*b(:,2)-a(:,2).*b(:,1)];
  nrmlev = @(p1,p2,p3) cross2(p2-p1,p3-p1);
  unrml = @(V,F) normalizerow(nrmlev(V(F(:,1),:),V(F(:,2),:),V(F(:,3),:)));

  while true
    % Compute normal averaged over all texture-space pixels landing in each
    % triangle: ~/uniform quadrature in texture space
    [X,Y] = meshgrid(linspace(0,1,size(nim,2)),linspace(1,0,size(nim,1)));
    I = in_element_aabb(TV(:,1:2),TF,[X(:) Y(:)]);
    Nnim = reshape(nim,[],3);
    TNsum = sparse( ...
      repmat(I(I>0),1,3),repmat(1:3,sum(I>0),1),Nnim(I>0,:),size(TF,1),3);
    TNcount =  sparse(I(I>0),1,1,size(TF,1),1);
    TN = normalizerow(full(TNsum./TNcount)*2-1);
    % Some triangles might not have received any samples
    N = unrml(V,F);
    TN(any(isnan(TN),2),:) = N(any(isnan(TN),2),:);
    A = diag(sparse(doublearea(V,F)));
    if fix_symmetry
      % flip if aggregate normal points in opposite direction from initial
      % normal
      D = A*sum(TN.*N,2);
      [C,CF] = connected_components(TF);
      CD = accumarray(CF',D)./accumarray(CF',diag(A));
      %tsurf(TF,[TV(:,1:2) C'],'CData',1*(CD(CF)<0.8));
      %T = normals(TV(:,1:2)*eye(2,3),TF);
      %tsurf(TF,[TV(:,1:2) C'],'CData',1*(T(:,3)>0));
      flipped = (CD(CF)<0.8);
      tsurf(TF(flipped,:),[TV(:,1:2) C'],'FaceVertexCData',N(flipped,:)*0.5+0.5);
      hold on;
      surf(X,Y,0*X,falpha(1,0),'FaceColor','texturemap','CData',nim);
      hold off;
      %TN(D<0,:) = -TN(D<0,:);
      error
    end
    ATN = A*TN;

    vec = @(X) X(:);
    Asqr = @(X) sum(sum(X.*(A*X)));
    obj_fun = @(V) Asqr(unrml(V,F)-TN);
    grad_fun = @(V) normal_gradient(V,F,vec(A*unrml(V,F))-vec(ATN));

    r3 = @(X) reshape(X,[],3);
    % This is awful...
    %options = optimoptions('fminunc','SpecifyObjectiveGradient',true,'Display','iter');
    %[V,obj1] = fminunc(@(V) deal(obj_fun(r3(V)),grad_fun(r3(V))),vec(V),options);

    switch method
    case 'lbfgs'
      % This thing actually seems to work
      [V,obj1] = fminlbfgs( ...
        @(V) fun_and_grad(V,@(V) obj_fun(r3(V)),@(V) grad_fun(r3(V))),vec(V), ...
        struct('Display','final','GradObj','on','TolFun',1e-6,'TolX',1e-5));
    case 'gd'
      %% My gradient descent is faster
      %[V,obj1] = fminlbfgs( ...
      %  @(V) fun_and_grad(V,@(V) obj_fun(r3(V)),@(V) grad_fun(r3(V))),vec(V), ...
      %  struct('HessUpdate','steepdesc','Display','final','GradObj','on','TolX',5e-5));
      [V,obj1] = fmingd( ...
        @(V) fun_and_grad(V,@(V) obj_fun(r3(V)),@(V) grad_fun(r3(V))),vec(V));
    end
    V = r3(V);

    %obj1 = obj_fun(V);
    %prev_step = 0;
    %step_size = 1;
    %prev_V = V;
    %for iter = 1:max_iters
    %  %method = 'gradient-descent';
    %  %method = 'stochastic-gradient-descent';
    %  %method = 'momentum';
    %  method = 'nesterov';
    %  %method = 'gauss-newton';
    %  %method = 'levenberg-marquardt';
    %  switch method
    %  case {'gradient-descent','gauss-newton','levenberg-marquardt','momentum'}
    %    % Recompute current normals and gradients
    %    dNdVT = normal_gradient(V,F);
    %    grad = dNdVT*(vec(A*unrml(V,F))-vec(ATN));
    %    switch method
    %    case {'gauss-newton','levenberg-marquardt'}
    %      max_step_size = 100;
    %      % This works really well for small meshes but takes forever/diverges for
    %      % big meshes
    %      K = (dNdVT*A*dNdVT');
    %      switch method
    %      case 'levenberg-marquardt'
    %        K = (dNdVT*A*dNdVT');
    %        T = diag(diag(K));
    %        lambda = 1;
    %        K = K + lambda*T;
    %      end
    %      % This works really well for small meshes but takes forever/diverges for
    %      % big meshes
    %      step = reshape(K\(-grad),[],3);
    %    case 'gradient-descent'
    %      % Gradient descent
    %      step = -reshape(grad,[],3);
    %      max_step_size = 1;
    %    case 'momentum'
    %      % Fairly sensitive to mu
    %      mu = 0.5;
    %      step = mu*prev_step + (1-mu)*(-reshape(grad,[],3));
    %      % works better if I use the non line search step here...
    %      prev_step = step;
    %      max_step_size = 1;
    %    end
    %  case 'nesterov'
    %    % following AQP paper
    %    eta = 10;
    %    theta = (1-sqrt(1/eta))/(1+sqrt(1/eta));
    %    V = (1+theta)*V - theta*prev_V;
    %    dNdVT = normal_gradient(V,F);
    %    grad = dNdVT*(vec(A*unrml(V,F))-vec(ATN));
    %    step = -reshape(grad,[],3);
    %    max_step_size = 1;
    %  case 'stochastic-gradient-descent'
    %    % not impressed with this... maybe it'll be worthwhile if to keep things
    %    % in memory
    %    %
    %    % Especially doesn't seem to play well with the "stalling" criteria.
    %    batch_size = min(size(F,1),10000);
    %    P = randperm(size(F,1));
    %    P = P(1:batch_size);
    %    dNdVT = normal_gradient(V,F(P,:));
    %    grad = dNdVT*(vec(A(P,P)*unrml(V,F(P,:)))-vec(ATN(P,:)));
    %    step = -reshape(grad,[],3);
    %    max_step_size = 1;
    %  end
    %  prev_V = V;
    %  [step_size,V,obj,step] = line_search(obj_fun,@(V) V,V,step,max_step_size);
    %  max_step_size_seen = max(max_step_size_seen,step_size);
    %  V = reshape(V,[],3);
    %  if step_size == 0
    %    warning('step size == 0');
    %    break;
    %  end
    %
    %  obj0 = obj1;
    %  obj1 = obj_fun(V);
    %  if ~quiet
    %    fprintf('  iter: % 5d, obj: %g\n',iter,obj1);
    %  end
    %  if obj0 - obj1 < 1e-5
    %    break
    %  end
    %end
    %if ~quiet
    %  fprintf('converged after %d iterations to %g\n',iter,obj1);
    %end

    if ~quiet
      tsurf(F,V,fsoft,'FaceVertexCData',repmat(orange,size(F,1),1),falpha(1,0.1));
      view(2);axis equal;camlight;
      %apply_ambient_occlusion([],'AddLights',false,'Factor',1);
      drawnow;
    end

    if progressive_upsampling
      [sel,V,F,TV,TF] = upsample_helper(V,F,TV,TF);
      if isempty(sel)
        break;
      end
    else
      % only one iteration needed: we upsampled everything at the beginning
      break;
    end
  end

  obj = obj_fun(V);
  if ~quiet
    fprintf('Final objective: %g\n',obj);
  end

  if ~quiet
    fprintf('max_step_size_seen: %g\n',max_step_size_seen);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [ names, sizes ] = cvx_dispvar( v, name, isdual )
[ names, sizes ] = dispvar_( v, name, isdual );
if nargout <= 1 && ~isempty( names ),
    mlen = max( cellfun( @numel, names ) ) + 3;
    spc = ' '; spc(1,2:mlen+3) = ' ';
    names = cellfun( @(x)[x,spc(1:mlen-numel(x))], names, 'UniformOutput', false );
    names = strcat( names, sizes );
end
if nargout == 0,
    names = sprintf( '   %s\n', names{:} );
    fprintf( 1, '\n%s\n', names );
    clear names
end

function [ names, sizes ] = dispvar_( v, name, isdual )
names = {};
sizes = {};
switch class( v ),
    case 'struct',
        if numel( v ) > 1,
            sv = size( v );
            sstr = struct( 'type', '()', 'subs', {cell(1,length(sv))} );
            isvec = all(sv(2:end)==1);
            for k = 1 : prod(sv),
                [ sstr.subs{:} ] = ind2sub( sv, k );
                if isvec,
                    label = sprintf( '%s(%d)', name, k );
                else
                    label = sprintf( ',%d', sstr.subs{:} );
                    label = sprintf( '%s(%s)', name, label(2:end) );
                end
                [ name2, size2 ] = dispvar_( v(k), sstr, label, isdual );
                names = [ names, name2 ]; %#ok
                sizes = [ sizes, size2 ]; %#ok
                if k == 1, name(:) = ' '; end
            end
        else
            fn = fieldnames( v );
            if ~isempty( name ), name( end + 1 ) = '.'; end
            names = {}; sizes = {};
            for k = 1 : length( fn ),
                [ name2, size2 ] = dispvar_( v.(fn{k}), [ name, fn{k} ], isdual );
                names = [ names, name2 ]; %#ok
                sizes = [ sizes, size2 ]; %#ok
                if k == 1, name(1:end-1) = ' '; end
            end
        end
    case 'cell',
        sv = size( v );
        sstr = cell( 1, length(sv) );
        isvec = all(sv(2:end)==1);
        for k = 1 : prod(sv),
            [ sstr{:} ] = ind2sub( sv, k );
            if isvec,
                label = sprintf( '%s{%d}', name, k );
            else
                label = sprintf( ',%d', sstr{:} );
                label = sprintf( '%s{%s}', name, label(2:end) );
            end
            [ name2, size2 ] = dispvar_( v{k}, label, isdual );
            names = [ names, name2 ]; %#ok
            sizes = [ sizes, size2 ]; %#ok
            if k == 1, name( 1 : end ) = ' '; end
        end
    case 'cvx';
        names{1} = name;
        sizes{1} = type(v);
    case 'double',
        names{1} = name;
        if isdual,
            sizes{1} = 'unassigned';
        else
            sizes{1} = [ type(cvx(v)), ' constant' ];
        end
    otherwise,
        names{1} = name;
        sizes{1} = 'unknown';
end

author: yu-jiang

repo: radpbook

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


text:
function model = gpReversibleDynamicsCreate(q, d, latentVals, options)

% GPREVERSIBLEDYNAMICSCREATE Create a reversible dynamics model.
% FORMAT
% DESC creates a Gaussian process model for dealing with reversible dynamics
% in the latent space of a GP-LVM.
% ARG q : the latent space dimension.
% ARG q : the latent space dimension.
% ARG X : the latent variables.
% ARG options : options structure as defined by gpOptions.m.
% RETURN model : model structure containing the Gaussian process.
%
% SEEALSO : gpCreate, gpReversibleDynamicsOptions,
% gpReversibleDynamicsLatentGradients, gpReversibleDynamicsSetLatentValues,
% gpReversibleDynamicsLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2005

% FGPLVM

if nargin < 4
  options = gpReversibleDynamicsOptions('ftc');
end

diffs = latentVals(2:end, :) - latentVals(1:end-1, :);
X = [latentVals(2:end-1, :) diffs(1:end-1, :)];
y = diffs(2:end, :);
model = gpCreate(2*q, d, X, y, options);
model.learn = 0;
model.type = 'gpReversibleDynamics';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function imageplot(M,str, a,b,c)

% imageplot - diplay an image and a title
%
% Example of usages:
%   imageplot(M);
%   imageplot(M,title);
%   imageplot(M,title,1,2,1);   % to make subplot(1,2,1);
%
%   If you want to display several images:
%       imageplot({M1 M2}, {'title1', 'title2'});
%
%   Copyright (c) 2007 Gabriel Peyre

if nargin<2
    str = [];
end

nbdims = 2;
if size(M,1)==1 || size(M,2)==1
    nbdims = 1;
end

if iscell(M)
    q = length(M);
    if nargin<5
        c = 1;
        a = ceil(q/4);
        b = ceil(q/a);
    end
    if (c-1+q)>(a*b)
        warning('a and c parameters not large enough');
        a = ceil((c-1+q)/4);
        b = ceil((c-1+q)/a);
    end
    for i=1:q
        if iscell(str)
            str1 = str{i};
        else
            str1 = str;
        end
        imageplot(M{i},str1, a,b,c-1+i);
    end
    global axlist;
    if not(isempty(axlist))
        linkaxes(axlist, 'xy');
    end
    return;
end

if nargin==5
    global axlist;
    global imageplot_size;
    if c==1 || isempty(imageplot_size) || imageplot_size~=size(M,1)
        clear axlist;
        global axlist;
        axlist = [];
        imageplot_size = size(M,1);
    end
    axlist(end+1) = subplot(a,b,c);
end

if nbdims==1
    plot(M); axis tight;
else
    if size(M,3)==2
        M = cat(3,M, zeros(size(M,1),size(M,2)));
    end
    if size(M,3)==1
        colormap gray(256);
    else
        colormap jet(256);
    end
    imagesc(rescale(M)); axis image; axis off;
end
if not(isempty(str))
    title(str);
end

if nargin==5 && c==a*b
    linkaxes(axlist, 'xy');
end

author: gpeyre

repo: matlab-toolboxes


text:
function ex9mbvp
%EX9MBVP  Example 9 of the BVP tutorial, solved as a multi-point BVP
%   This boundary value problem is the subject of Chapter 8 of
%   C.C. Lin and L.A. Segel, Mathematics Applied to Deterministic
%   Problems in the Natural Sciences, SIAM, Philadelphia, 1988.
%   The ODEs
%
%      v' = (C - 1)/n
%      C' = (vC - min(x,1))/eta
%
%   are solved on the interval [0, lambda]. The boundary conditions
%   are v(0) = 0, C(lambda) = 1, and continuity of v(x) and C(x) at
%   x = 1.  Example EX9BVP shows how this three-point BVP is
%   reformulated for solution with BVP4C present in MATLAB 6.0.
%   Starting with MATLAB 7.0, BVP4C solves multi-point BVPs directly.
%
%   The quantity of most interest is the emergent osmolarity Os =
%   1/v(lambda).  The parameters are related to another parameter
%   kappa by eta = lambda^2/(n*kappa^2).  Lin and Segel develop an
%   approximate solution for Os valid for "small" n.  Here the BVP is
%   solved for a range of kappa when lambda = 2 and n = 0.005.  The
%   computed Os is compared to the approximation of Lin and Segel.

% Copyright 2004, The MathWorks, Inc.

  % Known parameters, visible in nested functions.
  n = 5e-2;
  lambda = 2;

  % Initial mesh - duplicate the interface point x = 1.
  xinit = [0, 0.25, 0.5, 0.75, 1, 1, 1.25, 1.5, 1.75, 2];   lambda = 2;

  sol = bvpinit(xinit,[1 1]);

  fprintf(' kappa    computed Os  approximate Os \n')
  for kappa = 2:5
    eta = lambda^2/(n*kappa^2);
    % After creating function handles, the new value of eta
    % will be used in nested functions.
    sol = bvp4c(@ex9mode,@ex9mbc,sol);

    K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
    approx = 1/(1 - K2);
    computed = 1/sol.y(1,end);
    fprintf('  %2i    %10.3f    %10.3f \n',kappa,computed,approx);
  end

  figure
  plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--')
  legend('v(x)','C(x)')
  title('A three-point BVP.')
  xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.'])
  ylabel('v and C')

  % --------------------------------------------------------------------------
  % Nested functions
  %

  function dydx = ex9mode(x,y,region)
  %EX9MODE  ODE function for Example 9 of the BVP tutorial.
  % Here the problem is solved directly, as a three-point BVP.
    dydx = zeros(2,1);
    dydx(1) = (y(2) - 1)/n;
    % The definition of C'(x) depends on the region.
    switch region
      case 1    % x in [0 1]
        dydx(2) = (y(1)*y(2) - x)/eta;
      case 2    % x in [1 lambda]
        dydx(2) = (y(1)*y(2) - 1)/eta;
    end
  end % ex9mode

  % --------------------------------------------------------------------------

  function res = ex9mbc(YL,YR)
  %EX9MBC  Boundary conditions for Example 9 of the BVP tutorial.
  % Here the problem is solved directly, as a three-point BVP.
    res = [ YL(1,1)             % v(0) = 0
            YR(1,1) - YL(1,2)   % continuity of v(x) at x = 1
            YR(2,1) - YL(2,2)   % continuity of C(x) at x = 1
            YR(2,end) - 1    ]; % C(lambda) = 1
  end % ex9mbc

  % --------------------------------------------------------------------------

end  % ex9mbvp

author: Sable

repo: mcbench-benchmarks


text:
function pcut = addEvalVariableCuts(p)

pcut = p;
if ~isempty(p.evalMap)
    pcut = emptyNumericalModel;
    for i = 1:length(p.evalMap)
        y = p.evalVariables(i);
        x = p.evalMap{i}.variableIndex;
        xL = p.lb(x);
        xU = p.ub(x);

        % Generate a convex hull polytope
        if xL<xU
            if ~isempty(p.evalMap{i}.properties.convexhull)
                % A convex hull generator function is available!
                % Might be able to reuse hull from last run node
                if isfield(p.evalMap{i},'oldhull') && isequal(p.evalMap{i}.oldhull.xL,xL) && isequal(p.evalMap{i}.oldhull.xU,xU)
                    [Ax,Ay,b,K] = getOldHull(p,i);
                else
                    [Ax,Ay,b,K,p] = updateHull(xL,xU,p,i);
                    if isempty(Ax)
                        % Operator bounder does not cover this interval so
                        % use the sample-based instead
                        [Ax,Ay,b,K,p] = convexhullSampled(xL,xU,p,i);
                    end
                end
            else
               [Ax,Ay,b,K] = convexhullSampled(xL,xU,p,i);
            end
            if ~isempty(b)
                if isempty(K)
                    % Compatibility with old code
                    K.f = 0;
                    K.l = length(b);
                end
                F_structemp = zeros(size(b,1),length(p.c)+1);
                F_structemp(:,1+y) = -Ay;
                F_structemp(:,1+x) = -Ax;
                F_structemp(:,1) = b;
                localModel = createNumericalModel(F_structemp,K);
                pcut = mergeNumericalModels(pcut,localModel);
            end
        end
    end

    pcut = mergeNumericalModels(p,pcut);
end

function [Ax,Ay,b,K] = getOldHull(p,i);

Ax = p.evalMap{i}.oldhull.Ax;
Ay = p.evalMap{i}.oldhull.Ay;
b = p.evalMap{i}.oldhull.b;
K = p.evalMap{i}.oldhull.K;

function [Ax,Ay,b,K,p] = updateHull(xL,xU,p,i);

try
    [Ax,Ay,b,K]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
catch
    [Ax,Ay,b]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
    if ~isempty(Ax)
        problem = find(any(isinf([Ax Ay b]),2) | any(isnan([Ax Ay b]),2));
        Ax(problem,:) = [];
        Ay(problem,:) = [];
        b(problem) = [];
    end
    K = [];
end
p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i);

function p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i)
p.evalMap{i}.oldhull.xL = xL;
p.evalMap{i}.oldhull.xU = xU;
p.evalMap{i}.oldhull.Ax = Ax;
p.evalMap{i}.oldhull.Ay = Ay;
p.evalMap{i}.oldhull.b = b;
p.evalMap{i}.oldhull.K = K;

function [Ax,Ay,b,K,p] = convexhullSampled(xL,xU,p,i)

if length(xL)>1
    Ax = [];
    Ay = [];
    b = [];
    K = [];
    return
end
% sample function
z = linspace(xL,xU,100);

if isequal(p.evalMap{i}.fcn,'power_internal2')
    % Special code for automatically converting sigmonial
    % terms to be solvable with bmibnb
    fz = feval(p.evalMap{i}.fcn,z,p.evalMap{i}.arg{2});

else
    arg = p.evalMap{i}.arg;
    arg{1} = z;
    fz = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
    % end
    [minval,minpos] = min(fz);
    [maxval,maxpos] = max(fz);
    xtestmin = linspace(z(max([1 minpos-5])),z(min([100 minpos+5])),100);
    xtestmax = linspace(z(max([1 maxpos-5])),z(min([100 maxpos+5])),100);
    arg{1} = xtestmin;
    fz1 = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
    arg{1} = xtestmax;
    fz2 = real(feval(p.evalMap{i}.fcn,arg{1:end-1}));
    z = [z(:);xtestmin(:);xtestmax(:)];
    fz = [fz(:);fz1(:);fz2(:)];
    [z,sorter] = sort(z);
    fz = fz(sorter);
    [z,ii,jj]=unique(z);
    fz = fz(ii);
end

[Ax,Ay,b] = convexhullFromSampled(z,fz,xL,xU);
K = [];

p = saveOldHull(xL,xU,Ax,Ay,b,K,p,i);

author: shamilmamedov


text:
function varargout = ma_MissionNotesGUI(varargin)
% MA_MISSIONNOTESGUI MATLAB code for ma_MissionNotesGUI.fig
%      MA_MISSIONNOTESGUI, by itself, creates a new MA_MISSIONNOTESGUI or raises the existing
%      singleton*.
%
%      H = MA_MISSIONNOTESGUI returns the handle to a new MA_MISSIONNOTESGUI or the handle to
%      the existing singleton*.
%
%      MA_MISSIONNOTESGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MA_MISSIONNOTESGUI.M with the given input arguments.
%
%      MA_MISSIONNOTESGUI('Property','Value',...) creates a new MA_MISSIONNOTESGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ma_MissionNotesGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to ma_MissionNotesGUI_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 ma_MissionNotesGUI

% Last Modified by GUIDE v2.5 18-Jul-2015 20:44:51

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

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

% --- Executes just before ma_MissionNotesGUI is made visible.
function ma_MissionNotesGUI_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 ma_MissionNotesGUI (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);
handles.ma_MainGUI = varargin{1};
guidata(hObject, handles);

setupGui(handles);

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

function setupGui(handles)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    set(handles.missionNotesText,'String',maData.notes);

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

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

function missionNotesText_Callback(hObject, eventdata, handles)
% hObject    handle to missionNotesText (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 missionNotesText as text
%        str2double(get(hObject,'String')) returns contents of missionNotesText as a double

% --- Executes during object creation, after setting all properties.
function missionNotesText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to missionNotesText (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in saveCloseButton.
function saveCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    maData = getappdata(handles.ma_MainGUI,'ma_data');
    maData.notes = get(handles.missionNotesText,'String');
    setappdata(handles.ma_MainGUI,'ma_data',maData);

    close(handles.ma_MissionNotesGUI);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function i4vec_print ( n, a, title )

%*****************************************************************************80
%
%% I4VEC_PRINT prints an I4VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, integer A(N), the vector to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '%6d: %6d\n', i, a(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function kappa = surfaceCurvature(kappa1, kappa2, theta)
%SURFACECURVATURE Curvature on a surface from angle and principal curvatures.
%
%   usage:
%   KAPPA = surfaceCurvature(KAPPA1, KAPPA2, THETA)
%   return the curvature KAPPA of surface with respect to direction THETA.

%   KAPPA1 and KAPPA2 are the principal curvatures of the surface at the
%   considered point. THETA is angle of direction relative to angle of
%   first principal curvature KAPPA1.
%
%   Examples:
%   K = surfaceCurvature(KAPPA1, KAPPA2, 0) returns KAPPA1.
%   K = surfaceCurvature(KAPPA1, KAPPA2, pi/2) returns KAPPA2.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 07/04/2004.
%

%   HISTORY
%   20/04/2004 change name and add doc.
%   14/06/2004 correct creation date

kappa = kappa1 * cos(theta).^2 + kappa2 * sin(theta).^2;


text:
function out = spm_run_tissue_volumes(cmd, job)
% SPM job execution function for Tissue Volumes
%
% See also: spm_cfg_tissue_volumes, spm_summarise
%__________________________________________________________________________
% Copyright (C) 2013-2018 Wellcome Trust Centre for Neuroimaging

% Ged Ridgway
% $Id: spm_run_tissue_volumes.m 7460 2018-10-29 15:55:12Z john $

switch lower(cmd)
    %----------------------------------------------------------------------
    case 'exec'

        %%
        mat = job.matfiles;
        T   = job.tmax;
        msk = char(job.mask);
        outf = job.outf;
        if isempty(msk)
            msk = 'all';
        end

        %%
        N = numel(mat);
        vol = nan(N, T);
        for n = 1:N
            res = load(mat{n});

            % check for previously computed volumes (to save time if there)
            if isfield(res, 'volumes') && strcmp(res.volumes.mask, msk) ...
                    && numel(res.volumes.litres) >= T
                vol(n, :) = res.volumes.litres(1:T);
                continue
            end

            % determine number of tissue classes Kb
            if isfield(res, 'mg'),
                Kb  = max(res.lkp);
            else
                Kb  = size(res.intensity(1).lik, 2);
            end
            tc  = false(Kb, 4);

            % look for existing mwc files
            mwc = cell(T, 1);
            for t = 1:T
                mwc{t} = spm_file(res.image(1).fname, 'prefix', ['mwc' num2str(t)],'ext','nii');
                if ~exist(mwc{t}, 'file')
                    tc(t, 4) = true; % i.e. need to produce this mwc
                end
            end

            % produce mwc files if required
            if any(tc(:, 4))
                for f = 1:numel(res.image)
                    if ~exist(res.image(f).fname, 'file')
                        error('Original image no longer found at:\n%s\n', res.image(f).fname)
                    end
                end
                spm_preproc_write8(res, tc);
            end

            % compute tissue volumes
            vol(n, :) = spm_summarise(mwc, msk, 'litres');

            % add to mat file for future reuse
            volumes.litres  = vol(n, :);
            volumes.mask    = msk;
            save(mat{n}, 'volumes', '-append')

            % if mwc newly created above, delete now
            for t = 1:T
                if tc(t, 4), spm_unlink(mwc{t}), end
            end
        end

        %% Put into output structure for use with dependencies
        for t = 1:T
            out.(sprintf('vol%d', t)) = vol(:, t);
        end
        out.vol_sum = sum(vol, 2); % (total intracranial volume if T=1:3)

        %% Optionally save in CSV format
        if ~isempty(outf)
            [pth, nam, ext] = spm_fileparts(outf);
            if isempty(ext), ext = '.csv'; end
            fnm = fullfile(pth, [nam ext]);
            fid = fopen(fnm, 'wt');
            if fid < 0, error('Failed to open %s\n', fnm); end
            delim = ',';
            fprintf(fid, 'File');
            fprintf(fid, [delim 'Volume%d'], 1:T);
            for n = 1:N
                fprintf(fid, '\n''%s''', mat{n});
                fprintf(fid, [delim '%d'], vol(n, :));
            end
            fprintf(fid, '\n');
            fclose(fid);
        end

        %% Display in command window
        fprintf('\nSegmentation files:\n');
        fprintf('\t%s\n', mat{:});
        fprintf('\nVolumes (litres):\n');
        disp(vol);
        %------------------------------------------------------------------
    case 'vout'
        try
            T = job.tmax;
        catch
            T = 3;
        end
        out(T+1) = cfg_dep;
        for t = 1:T
            out(t).sname        = num2str(t);
            out(t).src_output   = substruct('.', sprintf('vol%d', t));
            out(t).tgt_spec     = cfg_findspec({
                {'strtype','e', 'strtype','r'}
                });
        end
        out(T+1).sname      = 'Sum';
        out(T+1).src_output = substruct('.', 'vol_sum');
        out(T+1).tgt_spec   = cfg_findspec({
            {'strtype','e', 'strtype','r'}
            });
end

author: spm