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

% Last Modified by GUIDE v2.5 15-Jan-2021 12:17:00

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditParallelToFrameCoordSysGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditParallelToFrameCoordSysGUI_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 lvd_EditParallelToFrameCoordSysGUI is made visible.
function lvd_EditParallelToFrameCoordSysGUI_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 lvd_EditParallelToFrameCoordSysGUI (see VARARGIN)

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

    centerUIFigure(hObject);

    coordSys = varargin{1};
    setappdata(hObject, 'coordSys', coordSys);

    lvdData = varargin{2};
    setappdata(hObject, 'lvdData', lvdData);

	populateGUI(handles, coordSys);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, coordSys)
    set(handles.coordSysNameText,'String',coordSys.getName());
    setappdata(handles.lvd_EditFixedInFrameVectorGUI,'frame',coordSys.frame);

    frame = coordSys.frame;
    handles.refFrameTypeCombo.String = ReferenceFrameEnum.getListBoxStr();
    handles.refFrameTypeCombo.Value = ReferenceFrameEnum.getIndForName(frame.typeEnum.name);
    handles.setFrameOptionsButton.TooltipString = sprintf('Current Frame: %s', frame.getNameStr());

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditParallelToFrameCoordSysGUI_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
    if(isempty(handles))
        varargout{1} = false;
    else
        coordSys = getappdata(hObject, 'coordSys');

        coordSys.setName(handles.coordSysNameText.String);

        coordSys.frame = getappdata(hObject,'frame');

        varargout{1} = true;
        close(handles.lvd_EditFixedInFrameVectorGUI);
    end

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

    if(isempty(errMsg))
        uiresume(handles.lvd_EditFixedInFrameVectorGUI);
    else
        msgbox(errMsg,'Invalid Coordinate System Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    if(isempty(strtrim(handles.coordSysNameText.String)))
        errMsg{end+1} = 'Vector name must contain more than white space and must not be empty.';
    end

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

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

% --- Executes during object creation, after setting all properties.
function coordSysNameText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to coordSysNameText 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 key press with focus on lvd_EditParallelToFrameCoordSysGUI or any of its controls.
function lvd_EditFixedInFrameVectorGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditParallelToFrameCoordSysGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditFixedInFrameVectorGUI);
    end

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

% Hints: contents = cellstr(get(hObject,'String')) returns refFrameTypeCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from refFrameTypeCombo
    frameUpdated(handles);

function frameUpdated(handles)
    lvdData = getappdata(handles.lvd_EditFixedInFrameVectorGUI, 'lvdData');

    frame = getappdata(handles.lvd_EditFixedInFrameVectorGUI,'frame'); %curElemSet
    bodyInfo = frame.getOriginBody();

    contents = cellstr(get(handles.refFrameTypeCombo,'String'));
    selFrameType = contents{get(handles.refFrameTypeCombo,'Value')};
    refFrameEnum = ReferenceFrameEnum.getEnumForListboxStr(selFrameType);

    switch refFrameEnum
        case ReferenceFrameEnum.BodyCenteredInertial
            newFrame = bodyInfo.getBodyCenteredInertialFrame();

        case ReferenceFrameEnum.BodyFixedRotating
            newFrame = bodyInfo.getBodyFixedFrame();

        case ReferenceFrameEnum.TwoBodyRotating
            if(frame.typeEnum == ReferenceFrameEnum.TwoBodyRotating)
                newFrame = frame;
            else
                if(not(isempty(bodyInfo.childrenBodyInfo)))
                    primaryBody = bodyInfo;
                    secondaryBody = bodyInfo.childrenBodyInfo(1);
                else
                    primaryBody = bodyInfo.getParBodyInfo();
                    secondaryBody = bodyInfo;
                end

                originPt = TwoBodyRotatingFrameOriginEnum.Primary;

                newFrame = TwoBodyRotatingFrame(primaryBody, secondaryBody, originPt, bodyInfo.celBodyData);
            end

        case ReferenceFrameEnum.UserDefined
            numFrames = lvdData.geometry.refFrames.getNumRefFrames();
            if(numFrames >= 1)
                geometricFrame = lvdData.geometry.refFrames.getRefFrameAtInd(1);
                newFrame = UserDefinedGeometricFrame(geometricFrame, lvdData);
            else
                newFrame = bodyInfo.getBodyCenteredInertialFrame();
                warndlg('There are no geometric frames available.  A body-centered inertial frame will be selected instead.');
            end

        otherwise
            error('Unknown reference frame type: %s', string(refFrameEnum));
    end

    if(not(isempty(newFrame)) && newFrame.typeEnum ~= refFrameEnum)
        refFrameEnum = newFrame.typeEnum;
        handles.refFrameTypeCombo.Value = ReferenceFrameEnum.getIndForName(refFrameEnum.name);
    end

    setappdata(handles.lvd_EditFixedInFrameVectorGUI,'frame',newFrame);
    handles.setFrameOptionsButton.TooltipString = sprintf('Current Frame: %s', newFrame.getNameStr());

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

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

% --- Executes on button press in setFrameOptionsButton.
function setFrameOptionsButton_Callback(hObject, eventdata, handles)
% hObject    handle to setFrameOptionsButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    frame = getappdata(handles.lvd_EditFixedInFrameVectorGUI,'frame');

    newFrame = frame.editFrameDialogUI(EditReferenceFrameContextEnum.ForGeometry);

    setappdata(handles.lvd_EditFixedInFrameVectorGUI,'frame',newFrame);

    frameUpdated(handles);


text:
function c=ref_rdftiii_1(f)
%REF_RDFTIII_1  Reference RDFT by FFT
%   Usage:  c=ref_rdftiii_1(f);
%
%   Compute RDFTII by doing a DFTIII and returning half the coefficients.
%   Only works for real functions.
%
%   The transform is orthonormal

L=size(f,1);
Lhalf=floor(L/2);
Lend=Lhalf*2;

cc=ref_dftiii(f);

c=zeros(size(f));

% Copy the cosine-part of the coefficients.
c(1:2:Lend,:)=sqrt(2)*real(cc(1:Lhalf,:));

% Copy the sine-part of the coefficients.
c(2:2:Lend,:)=-sqrt(2)*imag(cc(1:Lhalf,:));

% If f has an odd length, we must also copy the Niquest-wave
% (it is real)
if mod(L,2)==1
  c(end,:)=real(cc((L+1)/2,:));
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function l = laplacePriorLogProb(prior, x)

% LAPLACEPRIORLOGPROB Log probability of Laplace prior.

% PRIOR

% Compute log prior
l = -prior.precision*sum(abs(x)) - log(2) + log(prior.precision);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function e = p2e (u)
e = u(1:2,:) ./ ([1;1] * u(3,:));
return

author: strawlab

repo: MultiCamSelfCal


text:
function data = set_infer_chanlocs(data,disambiguation_rule)
% Infer the chanlocs subfields (positions and type) from labels.
% Data = set_infer_chanlocs(Data,DisambiguationRule)
%
% In:
%   Data   : some EEGLAB data set with a chanlocs field or the chanlocs field itself
%
%   DisambiguationRule : rule to apply to disambiguate montages when a given set of channel labels
%                        match multiple montage files similarly well. Can be:
%                        * 'coverage' : picks the one with maximum coverage
%                        * 'first' : picks the first one in the list of candidates (out of those that
%                                    have reasonable coverage)
%                        * 'deduce' : uses correlation analysis to deduce the best fit (can be slow)
%                        * filename : uses the montage with the given file name
%                        (default: 'deduce')
%
% Out:
%   Data   : original set with updated chanlocs field, or the updated chanlocs themselves
%
% Notes:
%   Parameters cannot be passed by name to this function.
%   To add a custom cap design to set_infer_chanlocs, see documentation in the
%   bcilab:/resources/caps folder.
%
%
% Examples:
%   % given a data set with labeled channels (but no locations), fill in standard channel locations
%   eeg = set_infer_chanlocs(eeg)
%
%   % given a .chanlocs field of a data set with only labels, fill in standard channel locations
%   chanlocs = set_infer_chanlocs(eeg.chanlocs)
%
%   % given a cell array of channel labels, obtain a chanlocs structure with locations filled in
%   chanlocs = set_infer_chanlocs({'C3','Cz','C4'})
%
%   % given a chanlocs file in a format recognized by EEGLAB, obtain a corresponding chanlocs struct
%   chanlocs = set_infer_chanlocs('mylocs.sfp')
%
% See also:
%   io_loadset
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-05

% path where *.cap files are located
cappath = 'resources:/caps/';

if nargin < 2
    disambiguation_rule = 'deduce'; end

% obtain the chanlocs
signal = [];
if isfield(data,{'head','parts'})
    data = exp_eval(data); end
if isfield(data,'chanlocs')
    locs = data.chanlocs;
    signal = data;
    if ~isfield(data.chaninfo,'labelscheme')
        data.chaninfo.labelscheme = '10-20'; end
else
    locs = data;
end

% channel labels given as a cell array
if iscellstr(locs)
    locs = struct('labels',locs); end;
if ischar(locs)
    locs = readlocs(env_translatepath(locs)); end
if ~isstruct(locs)
    error('The given data is not in an appropriate format (expected either a file name, a cell array of labels, a .chanlocs struct, or an EEGLAB dataset struct).'); end

% --- find best-matching cap if some locs missing ---

if (~all(isfield(locs,{'X','Y','Z'})) || all(cellfun('isempty',{locs.X}))) && ~isempty(locs)
    % get a list of all cap files
    capfiles = dir(env_translatepath(cappath));
    capfiles = capfiles(cellfun(@(n)~isempty(strfind(n,'.cap')),{capfiles.name}));
    caps = cellfun(@(n){[cappath n]},{capfiles.name});

    % for each cap file, try to match channel labels
    for c = 1:length(caps)
        try
            cap = caps{c};
            capdata{c} = hlp_microcache('chanlocs',@io_load,cap,'-mat');
            locdb = capdata{c}.CAP;
            fitlocs{c} = locs;
            switch capdata{c}.MATCH
                case {'exact','nocase_exact'}
                    % exact matching (both cases supported)
                    if ~isempty(strfind(capdata{c}.MATCH,'nocase'))
                        [found,idx_in_locdb,idx_in_res] = intersect(lower({locdb.labels}),lower({locs.labels})); %#ok<ASGLU>
                    else
                        [found,idx_in_locdb,idx_in_res] = intersect({locdb.labels},{locs.labels}); %#ok<ASGLU>
                    end
                    for fname = fieldnames(locdb)'
                        [fitlocs{c}(idx_in_res).(fname{1})] = locdb(idx_in_locdb).(fname{1}); end
                    [fitlocs{c}(idx_in_res).type] = deal('EEG');
                case {'substr','nocase_substr'}
                    % sub-string matching (both cases supported)
                    [found,idx_in_res,idx_in_locdb] = hlp_microcache('chanlocs',@do_matching,fitlocs{c},locdb,~isempty(strfind(capdata{c}.MATCH,'nocase'))); %#ok<ASGLU>
                    for fname = fieldnames(locdb)'
                        [fitlocs{c}(idx_in_res).(fname{1})] = locdb(idx_in_locdb).(fname{1}); end %#ok<*AGROW>
                    [fitlocs{c}(idx_in_res).type] = deal('EEG');
            end
            fraction_found(c) = 1-mean(cellfun('isempty',{fitlocs{c}.X}));
        catch e
            fprintf('Not matching channel labels with cap %s due to error: %s',cap,hlp_handleerror(e));
            fraction_found(c) = 0;
        end
    end

    % determine whether we have a clear-cut case
    [sorted_fractions,order] = sort(fraction_found,'descend');
    if max(sorted_fractions) < 0.5
        fprintf('WARNING: only %.0f%% of your channel labels match any known cap design; inferred locations might be wrong. Please consider adding your cap to set_infer_chanlocs.',100*max(sorted_fractions)); end
    if sorted_fractions(1) >= 1.5*sorted_fractions(2)
        % at least twice as many channels matched for this cap as for the next best one: keep it
        bestidx = order(1);
    elseif ~isempty(signal) && size(signal.data,1)==length(locs)
        % not a clear-cut case but we have a signal to work with
        switch disambiguation_rule
            case 'coverage'
                fprintf('Multiple known cap designs match to your data''s channel labels; picking the one with greatest coverage...\n');
                bestidx = order(1);
            case 'first'
                fprintf('Multiple known cap designs match to your data''s channel labels; picking the first one in the list...\n');
                bestidx = 1;
            case 'deduce'
                fprintf('Multiple known cap designs match to your data''s channel labels; determining the best fit...\n');
                quality = hlp_diskcache('montage_quality',@fit_quality,signal,sorted_fractions,order,fitlocs);
                % assess quality of best fit
                if max(quality)<0.3
                    fprintf('WARNING: Did not find a cap design that matches your data. You need to add channel locations for your cap to set_infer_chanlocs to use any function that requires these locations.\n');
                    bestidx = NaN;
                else
                    if max(quality)<0.5
                        fprintf('WARNING: The agreement between your data and the available cap designs is low and the chosen locations might be partially wrong. Either your data is too noisy for a reliable determination or your cap design differs from the known ones.\n');  end
                    bestidx = order(argmax(quality));
                end
            otherwise
                % assume that a montage file name was given
                match = ~cellfun('isempty',strfind(caps,disambiguation_rule));
                if nnz(match) == 1
                    bestidx = find(match);
                elseif nnz(match) > 1
                    fprintf('WARNING: multiple montages match your montage name: picking the first one.\n');
                    bestidx = find(match,1);
                else
                    % assume that this is a file name and try to load it
                    bestidx = NaN;
                    tmp = getfield(hlp_microcache('chanlocs',@io_load,disambiguation_rule,'-mat'),'CAP');
                    if ~isequal({tmp.labels},{locs.labels})
                        error('Currently the given cap file must match the available channel labels exactly.'); end
                end
        end
    elseif sorted_fractions(1)>sorted_fractions(2)
        % not a clear-cut case but minor evidence in favor of best cap
        fprintf('NOTE: multiple cap designs are compatible with your channel labels, using the best-matching one; inferred locations might be wrong.\n');
        bestidx = order(1);
    else
        % no clear winner
        fprintf('WARNING: multiple cap designs match your channel labels equally well, using the first cap in the list; inferred locations may be wrong.\n');
        bestidx = order(1);
    end

    % use the picked locations
    if ~isnan(bestidx)
        fprintf('  using montage %s.\n',caps{bestidx});
        locs = fitlocs{bestidx};
        if isstruct(data) && isfield(data,'data')
            data.chaninfo.labelscheme = capdata{bestidx}.SCHEME;
            data.chaninfo.nosedir = capdata{bestidx}.NOSE;
        end
    end

    % optionally optimize the head center
    if isfield(capdata{bestidx},'RECENTER') && capdata{bestidx}.RECENTER
        try
            if  mean(cellfun('isempty',{locs.X})) < 0.75
                locs = hlp_microcache('chanlocs',@pop_chanedit,locs,'eval','chans = pop_chancenter( chans, [],[]);'); end
        catch e
            fprintf('Failed trying to optimized the head center location for your cap montage due to error: %s\n',e.message);
        end
    end
end

% --- do post-processing ---

if ~isfield(locs,'labels')
    for k=data.nbchan:-1:1
        locs(k).labels = num2str(k); end
end
if ~isfield(locs,'type')
    locs(1).type = []; end
if ~isfield(locs,'theta')
    locs(1).theta = []; end

% reset 'unknown' channel types to empty
for k=1:length(locs)
    if strcmpi(locs(k).type,'unknown')
        locs(k).type = ''; end
end

% if some channel types are missing, try to look them up
empty_type = cellfun('isempty',{locs.type});
if any(empty_type)
    % assign all successfully retrieved EEG channels
    [locs(~cellfun('isempty',{locs.theta}) & empty_type).type] = deal('EEG');

    % assign the most common prefixes
    prefixes = {'eo','EOG','em','EMG','ec','ECG','ex','EXG','gsr','GSR'};
    for p=1:2:length(prefixes)
        [locs(strncmpi(prefixes{p},{locs.labels},length(prefixes{p})) & empty_type).type] = deal(prefixes{p+1}); end

    % also look for occurrence of EO, EM, ECG, EX, and GSR patterns
    patterns = {'EO','EOG','EM','EMG','ECG','ECG','EX','EXG','GSR','GSR','RE','EOG','LE','EOG'};
    for p=1:2:length(patterns)
       [locs(~cellfun('isempty',strfind({locs.labels},patterns{p})) & empty_type).type] = deal(patterns{p+1}); end

   % assign type 'EEG' to all numerically labeled channels (default assumption)
    [locs(cellfun(@(x)~isempty(str2num(x)),{locs.labels})).type] = deal('EEG'); %#ok<ST2NM>

    % assign 'EEG' to all unknown channels that are preceded and followed by EEG channels...
    filled_in = false;
    for k = find(cellfun('isempty',{locs.type}))
        % check if there is an EEG channel before this one (besides other unknowns)
        eeg_before = false;
        for j=k-1:-1:1
            if ~isempty(locs(j).type)
                if strcmp(locs(j).type,'EEG')
                    eeg_before = true; end
                break;
            end
        end
        % check if there is an EEG channel after this one (besides other unknowns)
        eeg_after = false;
        for j=k+1:length(locs)
            if ~isempty(locs(j).type)
                if strcmp(locs(j).type,'EEG')
                    eeg_after = true; end
                break;
            end
        end
        if eeg_before && eeg_after
            locs(k).type = 'EEG';
            filled_in = true;
        end
    end
    if filled_in
        disp('Note: some channels could not clearly be recognized as EEG (not in the 10-5 system), but were preceded and followed by EEG channels; assuming they are EEG.'); end

    % assign 'unknown' to the rest
    [locs(cellfun('isempty',{locs.type})).type] = deal('unknown');
end

if isfield(data,'chanlocs')
    data.chanlocs = locs;
else
    data = locs;
end

% function for sub-string matching between labels and a location struct array
function [found,idx_in_res,idx_in_locdb] = do_matching(fitlocs,locdb,use_lower)
found = {};
idx_in_res = [];
idx_in_locdb = [];
fitlabels = {fitlocs.labels};
dblabels = {locdb.labels};
if use_lower
    fitlabels = lower(fitlabels);
    dblabels = lower(dblabels);
end
for k=1:length(locdb)
    match = ~cellfun('isempty',strfind(fitlabels,dblabels{k}));
    if any(match)
        found{end+1} = locdb(k).labels;
        idx_in_res(end+1) = find(match,1);
        idx_in_locdb(end+1) = k;
    end
end

function quality = fit_quality(signal,sorted_fractions,order,fitlocs)
% first calculate bandpass-filtered signal X
attenuation = 80;
frequencies = min(1,2*[0.5 2 45 50]/signal.srate);
w = design_kaiser(frequencies(1),frequencies(2),attenuation,true);
B = design_fir(length(w)-1,[0 frequencies 1],[0 0 1 1 0 0],[],w);
try
    % try to process all data at once
    X = filtfilt_fast(B,1,signal.data')';
catch
    % fall back to channel-by-channel processing (if out of memory)
    for c=signal.nbchan:-1:1
        X(c,:) = filtfilt_fast(B,1,signal.data(c,:)')'; end
end

% for each cap design that matches reasonably well...
for k=1:nnz(sorted_fractions(1) < 2*sorted_fractions)
    o = order(k);
    % get matched locations
    tmp = fitlocs{o};
    [x,y,z] = deal({tmp.X},{tmp.Y},{tmp.Z});
    usable_channels = find(~cellfun('isempty',x) & ~cellfun('isempty',y) & ~cellfun('isempty',z));
    positions = [cell2mat(x(usable_channels));cell2mat(y(usable_channels));cell2mat(z(usable_channels))];
    % calculate channel interpolation matrix
    M = hlp_diskcache('montages',@interpolation_matrix,positions);
    % calculated interpolated signal Y
    Y = M*X(usable_channels,:);
    % calculated median of windowed correlation between X and Y
    window_len = 1*signal.srate;
    wnd = 0:window_len-1;
    offsets = round(1:window_len:size(Y,2)-window_len);
    W = length(offsets);
    corrs = [];
    for o=W:-1:1
        XX = X(usable_channels,offsets(o)+wnd)'; XX = bsxfun(@minus,XX,mean(XX));
        YY = Y(:,offsets(o)+wnd)'; YY = bsxfun(@minus,YY,mean(YY));
        corrs(:,o) = sum(XX.*YY)./(sqrt(sum(XX.^2)).*sqrt(sum(YY.^2)));
    end
    chancorr = median(corrs,2);
    % use average self-correlation across channels as quality index for this cap design
    quality(k) = mean(chancorr);
end

function M = interpolation_matrix(positions)
% interpolation_matrix_version<1.0> -- for the cache
fprintf('  calculating interpolation matrix for your montage; this is a one-time process');
range = 1:length(positions);
M = zeros(length(positions));
for k=range
    M(k,range~=k) = real(sphericalSplineInterpolate(positions(:,range~=k),positions(:,k)))';
    fprintf('.');
end
fprintf('\n');

repo: REST


text:
function [rxnScores, geneScores, hpaScores, arrayScores]=scoreModel(model,hpaData,arrayData,tissue,celltype,noGeneScore,multipleGeneScoring,multipleCellScoring,hpaLevelScores)
% scoreRxns
%   Scores the reactions and genes in a model based on expression data
%   from HPA and/or gene arrays
%
%   Input:
%   model               a model structure
%   hpaData             HPA data structure from parseHPA (opt if arrayData is
%                       supplied, default [])
%   arrayData           gene expression data structure (opt if hpaData is
%                       supplied, default [])
%       genes           cell array with the unique gene names
%       tissues         cell array with the tissue names. The list may not be
%                       unique, as there can be multiple cell types per tissue
%       celltypes       cell array with the cell type names for each tissue
%       levels          GENESxTISSUES array with the expression level for
%                       each gene in each tissue/celltype. NaN should be
%                       used when no measurement was performed
%       threshold       a single value or a vector of gene expression
%                       thresholds, above which genes are considered to be
%                       "expressed". (opt, by default, the mean expression
%                       levels of each gene across all tissues in arrayData
%                       will be used as the threshold values)
%   tissue              tissue to score for. Should exist in either
%                       hpaData.tissues or arrayData.tissues
%   celltype            cell type to score for. Should exist in either
%                       hpaData.celltypes or arrayData.celltypes for this
%                       tissue (opt, default is to use the best values
%                       among all the cell types for the tissue. Use [] if
%                       you want to supply more arguments)
%   noGeneScore         score for reactions without genes (opt, default -2)
%   multipleGeneScoring determines how scores are calculated for reactions
%                       with several genes ('best' or 'average')
%                       (opt, default 'best')
%   multipleCellScoring determines how scores are calculated when several
%                       cell types are used ('best' or 'average')
%                       (opt, default 'best')
%   hpaLevelScores      structure with numerical scores for the expression
%                       level categories from HPA. The structure should have a
%                       "names" and a "scores" field (opt, see code for
%                       default scores)
%
%
%   Output:
%   rxnScores       scores for each of the reactions in model
%   geneScores      scores for each of the genes in model. Genes which are
%                   not in the dataset(s) have -Inf as scores
%   hpaScores       scores for each of the genes in model if only taking hpaData
%                   into account. Genes which are not in the dataset(s)
%                   have -Inf as scores
%   arrayScores     scores for each of the genes in model if only taking arrayData
%                   into account. Genes which are not in the dataset(s)
%                   have -Inf as scores
%
%   Usage: [rxnScores, geneScores, hpaScores, arrayScores]=scoreModel(model,...
%               hpaData,arrayData,tissue,celltype,noGeneScore,multipleGeneScoring,...
%               multipleCellScoring,hpaLevelScores)

if nargin<3
    arrayData=[];
end
if nargin<5
    celltype=[];
else
    celltype=char(celltype);
end
if nargin<6
    noGeneScore=-2;
end
if nargin<7
    multipleGeneScoring='best';
else
    multipleGeneScoring=char(multipleGeneScoring);
end
if nargin<8
    multipleCellScoring='best';
else
    multipleCellScoring=char(multipleCellScoring);
end
if nargin<9
    %The first four are for APE, the other ones for staining
    hpaLevelScores.names={'High' 'Medium' 'Low' 'None' 'Strong' 'Moderate' 'Weak' 'Negative' 'Not detected'};
    hpaLevelScores.scores=[20 15 10 -8 20 15 10 -8 -8];
end

if isempty(hpaData) && isempty(arrayData)
    EM='Must supply hpaData, arrayData or both';
    dispEM(EM);
end
if ~strcmpi(multipleGeneScoring,'best') && ~strcmpi(multipleGeneScoring,'average')
    EM='Valid options for multipleGeneScoring are "best" or "average"';
    dispEM(EM);
end
if ~strcmpi(multipleCellScoring,'best') && ~strcmpi(multipleCellScoring,'average')
    EM='Valid options for multipleCellScoring are "best" or "average"';
    dispEM(EM);
end

%Throw an error if array data for only one tissue is supplied without
%specifying threshold values
if ~isempty(arrayData)
    if numel(unique(arrayData.tissues))<2
        if ~isfield(arrayData,'threshold') || isempty(arrayData.threshold)
            EM='arrayData must contain measurements for at least two celltypes/tissues since the score is calculated based on the expression level compared to the overall average';
            dispEM(EM);
        end
    end
end

%Process arrayData.threshold if necessary
if isfield(arrayData,'threshold') && (numel(arrayData.threshold) == 1)
    % if only a single gene threshold value is provided, then just
    % duplicate this value for all genes.
    arrayData.threshold = arrayData.threshold*ones(size(arrayData.genes));
end

%This is so that the code can ignore which combination of input data that is
%used
if isempty(arrayData)
    arrayData.genes={};
    arrayData.tissues={};
    arrayData.celltypes={};
    arrayData.levels=[];
    arrayData.threshold=[];
end
if isempty(hpaData)
    hpaData.genes={};
    hpaData.tissues={};
    hpaData.celltypes={};
    hpaData.levels={};
    hpaData.types={};
    hpaData.reliabilities={};
    hpaData.gene2Level=[];
    hpaData.gene2Type=[];
    hpaData.gene2Reliability=[];
end

%Check that the tissue exists
if ~ismember(upper(tissue),upper(hpaData.tissues)) && ~ismember(upper(tissue),upper(arrayData.tissues))
    EM='The tissue name does not match';
    dispEM(EM);
end
if any(celltype)
    %Check that both data types has cell type defined if that is to be used
    if ~isfield(hpaData,'celltypes') || ~isfield(arrayData,'celltypes')
        EM='Both hpaData and arrayData must contain cell type information if cell type is to be used';
        dispEM(EM);
    end
    if ~ismember(upper(celltype),upper(hpaData.celltypes)) && ~ismember(upper(celltype),upper(arrayData.celltypes))
        EM='The cell type name does not match';
        dispEM(EM);
    end
end

%Some preprocessing of the structures to speed up a little Remove all
%tissues that are not the correct one
J=~strcmpi(hpaData.tissues,tissue);

%If cell type is supplied, then only keep that cell type
if any(celltype)
    J=J | ~strcmpi(hpaData.celltypes,celltype);
end

hpaData.tissues(J)=[];
if isfield(hpaData,'celltypes')
    hpaData.celltypes(J)=[];
end
if isfield(hpaData,'gene2Level')
    hpaData.gene2Level(:,J)=[];
end
if isfield(hpaData,'gene2Type')
    hpaData.gene2Type(:,J)=[];
end
if isfield(hpaData,'gene2Reliability')
    hpaData.gene2Reliability(:,J)=[];
end

%Remove all genes from the structures that are not in model or that aren't
%measured in the tissue
if ~isempty(hpaData.genes) %This should not be necessary, but the summation is a 0x1 matrix and the other is []
    I=~ismember(hpaData.genes,model.genes) | sum(hpaData.gene2Level,2)==0;
else
    I=[];
end
hpaData.genes(I)=[];
if isfield(hpaData,'gene2Level')
    hpaData.gene2Level(I,:)=[];
end
if isfield(hpaData,'gene2Type')
    hpaData.gene2Type(I,:)=[];
end
if isfield(hpaData,'gene2Reliability')
    hpaData.gene2Reliability(I,:)=[];
end

I=strcmpi(arrayData.tissues,tissue);
%If cell type is supplied, then only keep that cell type
if any(celltype)
    I=I & strcmpi(arrayData.celltypes,celltype);
end

%Remove all genes from the structures that are not in model or that aren't
%measured in the tissue
J=~ismember(arrayData.genes,model.genes) | myAll(isnan(arrayData.levels(:,I)),2);
arrayData.genes(J)=[];
arrayData.levels(J,:)=[];
if isfield(arrayData,'threshold')
    arrayData.threshold(J) = [];
end

%Calculate the scores for the arrayData. These scores are calculated for
%each genes from its fold change between the tissue/celltype(s) in question
%and all other celltypes. This is a lower quality data than protein
%abundance, since a gene that is equally highly expressed in all cell types
%will have a score of 0.0. These scores are therefore only used for genes
%for which there is no HPA data available. The fold changes are transformed
%as min(log(x),10) for x>1 and max(log(x),-5) for x<1 in order to have
%negative scores for lower expressed genes and to scale the scrores to have
%somewhat lower weights than the HPA scores
tempArrayLevels=arrayData.levels;
tempArrayLevels(isnan(tempArrayLevels))=0;
if isfield(arrayData,'threshold') && ~isempty(arrayData.threshold)
    % if provided, the user-supplied expression threshold value(s) will be
    % used as the "average" expression level to which each gene is
    % compared.
    average=arrayData.threshold;
else
    average=sum(tempArrayLevels,2)./sum(~isnan(arrayData.levels),2);
end
if strcmpi(multipleCellScoring,'best')
    current=max(tempArrayLevels(:,I),[],2);
else
    current=sum(tempArrayLevels(:,I),2)./sum(~isnan(arrayData.levels(:,I)),2);
end
if ~isempty(current)
    aScores=5*log(current./average);
else
    aScores=[];
end
aScores(aScores>0)=min(aScores(aScores>0),10);
aScores(aScores<0)=max(aScores(aScores<0),-5);
aScores(isnan(aScores)) = -5;  % NaNs occur when gene expression is zero across all tissues

%Map the HPA levels to scores
[I, J]=ismember(upper(hpaData.levels),upper(hpaLevelScores.names));
if ~all(I)
    EM='There are expression level categories that do not match to hpaLevelScores';
    dispEM(EM);
end
[K, L, M]=find(hpaData.gene2Level);
scores=hpaLevelScores.scores(J);
if strcmpi(multipleCellScoring,'best')
    hScores=max(sparse(K,L,scores(M),numel(hpaData.genes),numel(hpaData.tissues)),[],2);
else
    hScores=mean(sparse(K,L,scores(M),numel(hpaData.genes),numel(hpaData.tissues)),2);
end

%Get the scores for the genes, only use HPA if available
geneScores=inf(numel(model.genes),1)*-1;
hpaScores=geneScores;
arrayScores=geneScores;

[I, J]=ismember(model.genes,hpaData.genes);
hpaScores(I)=hScores(J(I));
geneScores(I)=hScores(J(I));
[I, J]=ismember(model.genes,arrayData.genes);
arrayScores(I)=aScores(J(I));
geneScores(I & myIsInf(geneScores))=aScores(J(I & myIsInf(geneScores)));

%Remove the genes that have no data from the model
I=ismember(model.genes,hpaData.genes) | ismember(model.genes,arrayData.genes);
model.genes(~I)=[];
model.rxnGeneMat(:,~I)=[];

%Map the genes to the HPA/array genes
[hpaExist, hpaMap]=ismember(model.genes,hpaData.genes);
[arrayExist, arrayMap]=ismember(model.genes,arrayData.genes);

%Set the default scores for reactions without genes
rxnScores=ones(numel(model.rxns),1)*noGeneScore;

%Loop through the reactions and calculate the scores
for i=1:numel(model.rxns)
    %Check if it has genes
    I=find(model.rxnGeneMat(i,:));
    if any(I)
        %If any of the genes exist in hpaData, then don't use arrayData
        if any(hpaExist(I))
            %At least one gene was found in HPA
            if strcmpi(multipleGeneScoring,'best')
                rxnScores(i)=max(hScores(hpaMap(I(hpaExist(I)))));
            else
                rxnScores(i)=mean(hScores(hpaMap(I(hpaExist(I)))));
            end
        else
            %Use array data
            if any(arrayExist(I))
                %At least one gene was found in the array data
                if strcmpi(multipleGeneScoring,'best')
                    rxnScores(i)=max(aScores(arrayMap(I(arrayExist(I)))));
                else
                    rxnScores(i)=mean(aScores(arrayMap(I(arrayExist(I)))));
                end
            end
        end
    end
end
end

%This is because isinf and all returns 0x1 for empty set, which gives a
%concatenation error. Do like this instead of having many if statements
function y=myIsInf(x)
y=isinf(x);
if isempty(y)
    y=[];
end
end
function y=myAll(x,dim)
y=all(x,dim);
if isempty(y)
    y=[];
end
end

author: SysBioChalmers

repo: RAVEN


text:
function sts = match(item, spec)

% function sts = match(item, spec)
% This function is an implementation of find to search the cfg tree for
% certain entries.
%
% sts = match(item, spec)
% Spec must be a cell array of struct arrays with one or more fields. Each
% struct must contain two fields - 'name' and 'value'.
% An item matches, if it has a field with the specified field name and the
% contents of this field equals the contents of spec.value. If the field
% name is 'class', an item matches, if its class name is equal to
% spec.value.
% Matches within each struct array are OR-concatenated, while matches
% between struct arrays are AND-concatenated.
% An empty spec always matches.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: match.m 2919 2009-03-23 13:45:41Z volkmar $

rev = '$Rev: 2919 $'; %#ok

% match an empty spec
sts = true;

for k = 1:numel(spec)
    % Assume no match
    sts = false;
    for l = 1:numel(spec{k})
        if strcmp(spec{k}(l).name, 'class')
            sts = strcmp(spec{k}(l).value, class(item));
        elseif any(strcmp(spec{k}(l).name, mysubs_fields))
            sts = isequal(spec{k}(l).value, ...
                          subsref(item, substruct('.', spec{k}(l).name)));
        end;
        if sts
            % OR: success on first match
            break;
        end;
    end;
    if ~sts
        % AND: fail on first non-match
        break;
    end;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function d = p02_d ( m, id, c, w, n, x )

%*****************************************************************************80
%
%% P02_D evaluates a derivative for problem p02.
%
%  Discussion:
%
%    f(x) = 1 / product ( c(1:m)^(-2) + ( x(1:m) - w(1:m) )^2 )
%
%    Default values are:
%
%    c(1:m) = 1
%    w(1:m) = 0.5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    Reidel, 1987, pages 337-340,
%    ISBN: 9027725144,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer ID, the spatial coordinate to be differentiated.
%
%    Input, real C(M,1), W(M,1), the problem parameters.
%
%    Input, integer N, the number of points.
%
%    Input, real X(M,N), the evaluation points.
%
%    Output, real D(N,1), the value of the ID-th derivative component at X.
%
  d = ones ( n, 1 );
  for i = 1 : m
    d = d .* ( c(i) ^ (-2) + ( x(i,:)' - w(i) ) .^ 2 );
  end
  d = 1.0 ./ d .* ( - 2.0 ) .* ( x(id,:)' - w(id) ) ./ ...
    ( c(id) ^ (-2) + ( x(id,:)' - w(id) ) .^ 2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% housekeeping
clearvars
close all
clc

%% choose some options
% scramble the pseudo-random numbers or not
scramble_flag=true;
% optimal polynomial or theoretical ones
optimal_poly=true;
% debug or not
debug_flag=true;
%% read information about the function of interest
[objective,bounds]=ishigami();
% [objective,bounds]=satelli_sobol95();
lb=bounds(:,1);
ub=bounds(:,2);
npar=numel(lb);
param_names=strcat('pp_',num2str((1:npar)'));
%% choose the polynomial order
pol_order=4;
%% choose the expansion_order
expansion_order=3;
%% choose the number of samples
N=2^9;%2048;
%% option 1: generate samples and create the hdmr object

% theta=quasi_monte_carlo.sobol(npar,N,lb,ub,scramble_flag); % number of sub-intervals
theta=quasi_monte_carlo.sobol(lb,ub,N,scramble_flag); % number of sub-intervals
f=objective(theta);
objective_={f,theta};
%% option 2: let the hdmr object generate the samples for you
% objective_={objective,bounds,N};
%% construct the object
obj=hdmr(objective_,param_names,bounds,expansion_order,pol_order,optimal_poly);

%% estimate the object
profile off
profile on
obj=estimate(obj,debug_flag);
profile off
profile viewer
%% plot the fit insample
plot_fit(obj,'insample');
%% plot the fit out of sample
plot_fit(obj,'outofsample');
%% plot the individual effects
for ii=1:size(theta,1)
    figure();
    first_order_effect(obj,'insample',ii);
end

author: jmaih


text:
%% kalman tracking
function k = klmf(k)
for i      = 1 : size(k, 2)                         % for every kalman filter
k(i).s     = kalmanf(k(i).s);                       % update kalman filter
end
end
%% update kalman
function s = kalmanf(s)
s.x        = s.A * s.x;                             % prediction for state vector and covariance
s.P        = s.A * s.P * s.A' + s.Q;                % covariance of the state vector estimate
K          = s.P * s.H' / (s.H * s.P * s.H' + s.R); % compute Kalman gain factor
s.x        = s.x + K * (s.z - s.H * s.x);           % correction based on observation
s.P        = s.P - K * s.H * s.P;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function num_dist = dist_enum ( k, m )

%*****************************************************************************80
%
%% DIST_ENUM returns the number of distributions of indistinguishable objects.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer K, the number of distinguishable "slots".
%
%    Input, integer M, the number of indistinguishable objects.
%
%    Output, integer NUM_DIST, the number of distributions of M
%    indistinguishable objects about K distinguishable slots.
%
  cnk = combin ( m + k - 1, m );

  num_dist = round ( cnk );

  return
end


text:
function costs = surfprofile(problem, x, d1, d2, t1, t2)
% Plot the cost function as a surface over a 2-dimensional subspace.
%
% function surfprofile(problem, x, d1, d2, t1, t2)
% function costs = surfprofile(problem, x, d1, d2, t1, t2)
%
% Evaluates the cost function at points
%
%   gamma(t1, t2) = exponential_x(t1*d1 + t2*d2)
%
% where the exponential map at x is specified by problem.M.exp (retr is
% used instead if needed). d1 and d2 are two tangent vectors to problem.M
% at the point x. The values assigned to t1 and t2 are as specified in the
% two input vectors t1 and t2.
%
% If the function is called with an output, the plot is not drawn and the
% values of the cost are returned in a matrix of size
% length(t1)*length(t2). To plot a surf, call surf(t1, t2, costs.') (notice
% the transpose).

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sep. 1, 2014.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Verify that the problem description is sufficient.
    if ~canGetCost(problem)
        error('It seems no cost was provided.');
    end

    if isfield(problem.M, 'exp')
        expo = problem.M.exp;
        str = 'Exp';
    else
        expo = problem.M.retr;
        str = 'Retr';
    end

    storedb = StoreDB();
    linesearch_fun = @(ta, tb) getCost(problem, ...
                         expo(x, problem.M.lincomb(x, ta, d1, tb, d2)), ...
                         storedb);

    costs = zeros(length(t1), length(t2));
    for i = 1 : length(t1)
        for j = 1 : length(t2)
            costs(i, j) = linesearch_fun(t1(i), t2(j));
        end
    end

    if nargout == 0
        surf(t1, t2, costs.');
        xlabel('t1');
        ylabel('t2');
        zlabel(['f(' str '_x(t1*d1+t2*d2))']);
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function test_suite=SimTest_noddi_test
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;

function TestSetup
setenv('ISDISPLAY','0') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','1')

function test_Sim
disp('===========================================')
disp('Running simulation test for noddi');
disp('testing Simulation Single Voxel Curve...');

Model = str2func('noddi'); Model = Model();
savedModel_fname = fullfile(fileparts(which('qMRLab')),'Test','MoxUnitCompatible','static_savedModelsforRetrocompatibility',['noddi.qmrlab.mat']);
if ~exist(savedModel_fname,'file')
Model.saveObj(savedModel_fname);
else
Model = Model.loadObj(savedModel_fname);
end

disp(class(Model))
try Opt = button2opts(Model.Sim_Single_Voxel_Curve_buttons,1); end
try st = Model.st; catch, try st = mean([Model.lb(:),Model.ub(:)],2); catch, st = ones(length(Model.xnames),1); end; end
if exist('Opt','var') && length(Opt)>1
[Opt(:).SNR] = deal(1000);
else
Opt.SNR=1000;
end
for iopt=1:length(Opt) % Test all simulation options
disp(['Testing ' class(Model) ' simulation option:'])
disp(Opt(iopt))
FitResults = Model.Sim_Single_Voxel_Curve(st,Opt(iopt));
% Compare inputs and outputs
fnm=fieldnames(FitResults);
FitResults = rmfield(FitResults,fnm(~ismember(fnm,Model.xnames))); fnm=fieldnames(FitResults);
[~,FitResults,GroundTruth]=comp_struct(FitResults,mat2struct(st,Model.xnames),[],[],.30);
assertTrue(isempty(FitResults) & isempty(GroundTruth),evalc('FitResults, GroundTruth'))
end
disp ..ok

function TestTeardown
setenv('ISDISPLAY','') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [SetupStruc] = ISM_setup()
%%%%%%   the environment coefficients
SetupStruc.T60 = 0;      %T60, Options: 0, 0.3, 0.6, 0.9
SetupStruc.c = 340;
SetupStruc.Sign_compared = -1;   % choosing the standard signal, non-negtive corresponding to the signal in the select T60, or the original signal
%%%%%%    the microphone coefficients
SetupStruc.Size = 4.35*0.01;
SetupStruc.Channel_Num = 7;      %the number of mics
SetupStruc.Sign_Mid = 1;         %the exist of the center mic
SetupStruc.AddNoiseFlag = 0;
SetupStruc.NoiseSNR = 30;
SetupStruc.Height = 0; %set 0 use the sameheight data, set 1 use the differentheight data

repo: Voice-Separation-and-Enhancement


text:
function ll = hsvargplvmLogLikelihood(model)

F_leaves = hsvargplvmLogLikelihoodLeaves(model.layer{1});
F_nodes = hsvargplvmLogLikelihoodNode(model);
F_entropies = hsvargplvmLogLikelihoodEntropies(model);
% This refers to the KL quantity of the top node. The likelihood part is
% computed in hsvargplvmLogLikelihoodNode.
F_parent = hsvargplvmLogLikelihoodParent(model.layer{model.H});

ll = F_leaves + F_nodes + F_entropies + F_parent;

%--- NEW!!! TEST!
% %{
for h=1:model.H
    % If there's a prior on some parameter, add the (minus) corresponding
    % likelihood.
    for m=1:model.layer{h}.M
        ll = ll + vargplvmParamPriorLogProb(model.layer{h}.comp{m});
    end
end
% %}
%---

end

% The ln p(Y|X) terms
function F_leaves = hsvargplvmLogLikelihoodLeaves(modelLeaves)

F_leaves = 0;
for m=1:modelLeaves.M
    modelLeaves.comp{m}.onlyLikelihood = true;
    F_leaves = F_leaves + vargplvmLogLikelihood(modelLeaves.comp{m});
end
end

% The <ln p(X_h | X_{h-1})>_{q(X_{h-1}} nodes
function F_nodes = hsvargplvmLogLikelihoodNode(model)
F_nodes = 0;
for h=2:model.H
    % It's just like the leaves computation, the only difference is the
    % trace(Y*Y') term which now is replaced by an expectation w.r.t the
    % latent space of the previous layer. However, this replacement is done
    % in hsvargplvmUpdateStats and we dont have to worry here about it
    F_nodes = F_nodes + hsvargplvmLogLikelihoodLeaves(model.layer{h});
end

end

% The H_{q(X_h)} nodes, h ~= H
function F_entropies = hsvargplvmLogLikelihoodEntropies(model)
F_entropies = 0;
for h=1:model.H-1
    vardist = model.layer{h}.vardist;
    F_entropies = F_entropies + 0.5*(vardist.numData*vardist.latentDimension* ...
            (log(2*pi) + 1) + sum(sum(log(vardist.covars))));
end
if isfield(model, 'DEBUG_entropy') && model.DEBUG_entropy
    F_entropies = - F_entropies;
end

end

% The -KL[q(X_H) || p(X_H)]
function F_parent = hsvargplvmLogLikelihoodParent(modelParent)
%if modelParent.M > 1
%    warning('Not implemented multiple models in parent node yet')
%end
% Copied from vargplvmLogLikelihood:
if isfield(modelParent, 'dynamics') & ~isempty(modelParent.dynamics)
        % A dynamics model is being used.
        F_parent = modelVarPriorBound(modelParent);
        F_parent = F_parent + 0.5*modelParent.q*modelParent.N; %%% The constant term!!
else
    varmeans = sum(sum(modelParent.vardist.means.*modelParent.vardist.means));
    varcovs = sum(sum(modelParent.vardist.covars - log(modelParent.vardist.covars)));
    F_parent = -0.5*(varmeans + varcovs) + 0.5*modelParent.q*modelParent.N;
end
end

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function varargout=distitpf(varargin)
%
% For calling details please see v_distitpf.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_distitpf.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: distitpf.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_distitpf(varargin{:});
else
    v_distitpf(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function cadj=ref_spreadadj_1(coef)
%REF_SPREADADJ_1  Symbol of adjoint spreading function.
%   Usage: cadj=ref_spreadadj_1(c,number);
%
%   Development version by FJ for comparison of different implementations
%   cadj=SPREADADJ(c) will compute the symbol cadj of the spreading
%   operator that is the adjoint of the spreading operator with symbol c.
%
%   This is the implementation previously in the toolbox,
%   with the addition of an initialisation of cadj before the loop
%   leading to a huge speed improvement :

L=size(coef,1);

% Matlab cannot handle an FFT of a sparse matrix.
coef=ifft(full(coef));

cadj = zeros(L);
for ii=0:L-1
  for jj=0:L-1
    cadj(ii+1,jj+1)=conj(coef(mod(ii-jj,L)+1,mod(-jj,L)+1));
  end;
end;

cadj=fft(full(cadj));


text:
% Finds the number of k-neighbors (k links away) for every node
%
% INPUTS: adjacency matrix (nxn), start node index, k - number of links
% OUTPUTS: vector of k-neighbors indices
%
% GB: last updated, Oct 7 2012

function kneigh = kneighbors(adj,ind,k)

adjk = adj;
for i=1:k-1; adjk = adjk*adj; end;

kneigh = find(adjk(ind,:)>0);

author: aeolianine

repo: octave-networks-toolbox


text:
function [fileinfo,CHINFO]=cp42sto(stofile,graf);
% CP42STO   Read compressed MGCplus Data files, recorded with CP42
%
%           [fileinfo,CHINFO]=cp42sto([stofile],[graph])
%
%           stofile is the complete filename of the CP42 DAQ file.
%           graph is a string of length 1. If it is different from 'n',
%           then the measurement data of the file is shown.
%           Both input arguments are optional.
%
%           This function reads the compressed data files recorded on a
%           PC card harddisk on CP42. The input argument "stofile" is
%           optional, when it is missing, a user dialog for selecting
%           a file with extension *.st* is shown.
%
%           fileinfo is a struct, containing the FileID, the number of
%           channels "ChannelC", Byte/Line, the data format, the data offset,
%           the time format, the file comment, the comment of the current
%           parameterset, the information about the creation date and the
%           number of used slots.
%
%           CHINFO is a struct, containing information about the channel
%           number, the scaling, offset, unit, signalling bit, subchannel,
%           data, timeaxis, channel name and the channel code.
%
%           With PLOT(CHINFO(k).timeax,CHINFO(k).data) a measurement
%           signal can be displayed.
%           TITLE(CHINFO(k).ChannelName) adds the channel name as
%           graph title
%           YLABEL(CHINFO(k).Unit) puts the used unit on the Y axis.
%
%           There is a lot of information in the mea files about channels
%           and slots. Not all of these information is evaluated.
%           Want more ? Have a look at the word file "CP42MEA-FORMAT.DOC"
%
%           Only CP42 measurement files are supported. These DAQ files may
%           contain time information in following formats:
%           1) No Time Information
%           2) MGCplus Time (64 Bit)
%           3) NTP Time Format
%           The Time Formats "USB Frame Count" and 32 Bit MGC device Time
%           are not supported.
%
%           IMPORTANT:
%           Only Measurement files with format 4 BYTE INT (Intel) are
%           supported.

CHINFO=[];
CHSET=[];
fileinfo=[];

if nargin==0,               % No Input argument
    [fn,pn]=uigetfile('*.st*','Select a compressed CP42 data file');
    if fn==0,
        return
    end
    stofile=[pn,fn];
    graf='n';
end

if nargin==1,               % One Input argument
    if length(stofile)==1,  % If Graphic Option was specified...
        graf=stofile;
        [fn,pn]=uigetfile('*.st*','Select a compressed CP42 data file');
        if fn==0,
            return
        end
        stofile=[pn,fn];
    else                    % CP42 filename was specified
        graf='n';
    end
end

if exist(stofile)~=2,
    disp('Measurement file does not exist!');
    return;
end

fileinfo=struct('FileID','ChannelC','ByteLine','MVLines',...
    'MWformat','MRate','DataOffset','reduction',...
    'cmt','cmtpar','creation','usedslot','filename','filesize');
chaninfo=struct('Channelnumber','Scaling','Offset','Unit',...
    'SigBit','Subchannel','data','timeax',...
    'ChannelName','chncode','SigName','SigTypus');
chsetting=struct('ChannelName','chncode','ampcode','paraflg',...
    'tabquelle','acalflg',...
    'bridget','bridgev','uk','frequenz');

% Open CP42 measurement file (on PC) for reading.
fid = fopen(stofile,'rb');
fseek(fid,0,1);             % Evaluate the filesize in Byte
filesize=ftell(fid);
fseek(fid,0,-1);

% Read Header Block. The header block of all CP42 measurement files (also
% rainflow version CP42S30) is same: 8 long values
a=fread(fid,8,'long');

if a(1)~=6002,
    disp(['Invalid File identifier : ',int2str(a(1))]);
    fclose(fid);
    return;
end

% Only 4 Byte INT Intel format is supported. With newest versions of the
% HBM setup assistant, only this format is upported.
if a(5)~=1253,
    disp('Data not in 4 Byte INT (Intel) Format !');
    fclose(fid);
    return;
end

% The measurement rate is coded with values > 6000. To obtain the sample
% rate in Hz, a switch-case-expression is evaluated.
% Improved: In former version switch-case was used ! (now much shorter)
wrat=[(6300:1:6320), (6321:6341), (6345:6347), (6350:6361);1,2,5,10,25,50,60,75,100,150,200,300,...
        400,600,800,1200,1600,2400,3200,4800,9600,3,4,6,8,15,20,30,48,64,80,96,...
        120,160,192,240,320,384,480,640,960,1920,19200,38400,76800,...
        (1/2),(1/5),(1/10),(1/20),(1/50),(1/100),(1/200),(1/500),(1/1E3),(1/2E3),(1/5E3),(1/1E4)];

k0r=find(wrat(1,:)==a(6));
if isempty(k0r),
    disp(['unknown Measurement Rate : ',int2str(a(6))]);
    fclose(fid);
    return
end
SR=wrat(2,k0r);

noofsamp=floor((filesize-a(8)) ./a(3));
if a(4)~=noofsamp,              % Sometimes some files do not have
    disp('Attention!');         % info about no. of samples
    disp('Number of Samples contained');
    disp('in File seems to be invalid');
    disp(['File : ',int2str(a(4))]);
    disp(['calc : ',int2str(noofsamp)]);
    a(4)=noofsamp;
end                             % (File Repairing)

fileinfo.FileID=a(1);           % File-ID (6001)
fileinfo.ChannelC=a(2);         % Channel Count (incl. Time Channels)
fileinfo.ByteLine=a(3);         % Size MeasVal-Line / Bytes
fileinfo.MVLines=a(4);          % No. of MeasVal-Lines
fileinfo.MWformat=a(5);         % MeasVal-Format
                                % (1253 is required: 4 Byte INT Intel)
fileinfo.MRate=SR;              % Sample Rate in Hz
fileinfo.DataOffset=a(8);       % Data Offset in Byte
fileinfo.reduction=a(7);        % Compression factor
fileinfo.filename=stofile;
fileinfo.filesize=filesize;

% Read Block with channel parameters
for ch=1:a(2),
    chaninfo.Channelnumber=fread(fid,1,'long');     % Channelnumber
    chaninfo.Scaling=fread(fid,1,'float');          % Scale Factor
    chaninfo.Offset=fread(fid,1,'float');           % Offset
    chaninfo.Unit=setstr(fread(fid,[1,4],'char'));  % Unit
    chaninfo.SigBit=fread(fid,1,'short');           % Signal-Bitfeld

    sigbin=flipud(dec2bin(chaninfo.SigBit));        % Convert Signal Bitfield
                                                    % to binary representation
    if length(sigbin)>6,                            % Only the lower 6 Bit
        sigbin=sigbin(1:6);                         % contain sig. types!
    end
    numsigs=sum(sigbin);                            % No. of Signals / Channels
    kitypus=find(sigbin==1);                        % Positions of signals

    sucha=fread(fid,1,'short');                     % Subchannel No. Single
    if sucha==0,                                    % channel amps have got subch.
        chaninfo.Subchannel=1;                      % 0 which is replaced by 1
    else
        chaninfo.Subchannel=sucha;
    end

    for sgtype=1:numsigs,                           % Evaluate Signal-
        typaktuell=kitypus(sgtype);                 % type information
        switch typaktuell
            case {1}
                chaninfo.SigName=' (Gross)';
            case {2}
                chaninfo.SigName=' (Net)';
            case {3}
                chaninfo.SigName=' (PV1)';
            case {4}
                chaninfo.SigName=' (PV2)';
            case {5}
                chaninfo.SigName=' (CPV)';
            case {6}
                chaninfo.SigName=' (DIG)';
        end
        chaninfo.SigTypus=typaktuell;
        CHINFO=[CHINFO,chaninfo];
    end
end

% Read Block with amplifer settings
w0=a(2);
mgctime=0;

for ch=1:w0;
    qtmp=setstr(fread(fid,[1,47],'char'));  % Read ChannelName
    % The channel name contains the name in hyphens and some stupid signs.
    % The channel names are redesigned in the following lines...
    ki34=find(abs(qtmp)==34);               % The channel name is between
    if length(ki34)>=2,                     % two hyphens.
        qtmp=qtmp(ki34(1)+1:ki34(2)-1);
        qtmp=deblank(qtmp);
    end
    ki59=find(abs(qtmp)==59);               % Sometimes, the strings contains
    if ~isempty(ki59),                      % a final semicolon. It is removed
        qtmp=qtmp(1:max(ki59)-1);
    end
    qtmp=deblank(qtmp);                     % Remaining blanks are removed

    chsetting.ChannelName=qtmp;
    chsetting.chncode=fread(fid,1,'char');  % Read Channel Code
    chsetting.ampcode=fread(fid,1,'char');  % Read amplifier code
    chsetting.paraflg=fread(fid,1,'char');  % parflag (1..8: intern, 9: XM001)
    chsetting.tabquelle=fread(fid,1,'char');% Where is the chart from ?
    chsetting.acalflg=fread(fid,1,'char');  % AutoCalibration on / off ?
    fseek(fid,3,0);                         % Ignore loca_, meas_ and stat_flg
    chsetting.bridget=fread(fid,1,'int16'); % Read Bridgetype
    chsetting.bridgev=fread(fid,1,'char');  % Read Bridge supply voltage
    %statu=fseek(fid,167,0);
    statu=fseek(fid,166,0);
    sucha=fread(fid,1,'char');              % Subchannel No.
    if sucha==0,                            % Replace Subchannel 0 by subch. 1
        chsetting.uk=1;
    else
        chsetting.uk=sucha;
    end
    CHSET=[CHSET,chsetting];
end

for qz=1:length(CHSET),                     % Append Chan. Name and Sig. type
    chanakt=CHSET(qz).chncode;
    sucha=CHSET(qz).uk;
    for qz2=1:length(CHINFO),
        chan2=CHINFO(qz2).Channelnumber;
        sucha2=CHINFO(qz2).Subchannel;
        if ((chanakt==chan2) )%& (sucha==sucha2)),
            CHINFO(qz2).ChannelName=[CHSET(qz).ChannelName,CHINFO(qz2).SigName];
        end
    end
end

creation=setstr(fread(fid,[1,26],'char'));  % Date and Time of File Creation
fileinfo.creation=creation;                 % e. g. Wed Sep 03
                                            % 14:43:06 1997

comment=setstr(fread(fid,[1,80],'char'));   % Comment of the measurement
fileinfo.cmt=comment;

commentpar=setstr(fread(fid,[1,80],'char'));% Comment of active parameterset
fileinfo.cmtpar=commentpar;

slotnumber=fread(fid,1,'int');              % No. of used slots for meas.

fseek(fid,a(8),-1);                         % Goto Data Section

nrow=length(CHINFO);                        % If MGC Time Format (a(8)=7000),
                                            % then two time channels are
                                            % contained in CHINFO
nplus=0;

data=fread(fid,[nrow+nplus,a(4)],'int32');  % Read Data of all chan. & all sig.
data=double(data);                          % double for rescaling operation

for p=1:nrow
  data(p,:)=((data(p,:) .*CHINFO(p).Scaling) ./1966080000) - CHINFO(p).Offset;
  CHINFO(p).data=data(p,:);                 % Append MV's to Channel Info
end

                                            % Build "Time" channel (Only
                                            % Index)
timeax=(0:1:a(4)-1) *a(7);

for p=1:length(CHINFO),
    CHINFO(p).timeax=timeax;
end

fclose(fid);                                % Close measurement file

if nargout<2 | graf~='n',
    plotcp42sto(CHINFO);                    % Disp compressed Data
end
% END OF CP42STO

function plotcp42sto(CHINFO);
% PLOTCP42STO
%           plotcp42sto(CHINFO)
%
%           This function displays compressed measurement data of a CP42
%           sto-file in several subplots.

if isempty(CHINFO),                 % No input argument leads outside
    return;
end

cnt=length(CHINFO);
oldinfo=[-1,-1,-1];
graph=0;

% Evaluate Number of needed Subplots
for p=1:cnt,
    kanal=CHINFO(p).Channelnumber;
    unterkanal=CHINFO(p).Subchannel;
    signal=CHINFO(p).SigTypus;
    newinfo=[kanal,unterkanal,signal];

    if ~((newinfo(1)==oldinfo(1)) & (newinfo(2)==oldinfo(2))),
        graph=graph+1;
    end
    oldinfo=newinfo;
end
oldinfo=[-1,-1,-1];

% if more than 6 subplots are needed, the function is left.
if graph>6,
    disp('Too much information for complete graphical presentation');
    return;
end
% Build up the graph
graphakt=0;
set(gcf,'Color',[0.5,1,1]);
for p=1:cnt,
    kanal=CHINFO(p).Channelnumber;
    unterkanal=CHINFO(p).Subchannel;
    signal=CHINFO(p).SigTypus;
    newinfo=[kanal,unterkanal,signal];

    if ~((newinfo(1)==oldinfo(1)) & (newinfo(2)==oldinfo(2))),
        graphakt=graphakt+1;
        hold off;
        subplot(graph,1,graphakt);
        ploha=plot(CHINFO(p).timeax,CHINFO(p).data,'b-');
        set(gca,'FontSize',8,'Color',[1,1,1]);
        yha=ylabel(CHINFO(p).Unit);
        set(yha,'FontSize',8);
        hti=title(CHINFO(p).ChannelName);
        set(hti,'FontSize',8);
        hold on;
    else
        plot(CHINFO(p).timeax,CHINFO(p).data,'m-');
    end
    oldinfo=newinfo;
end
xha=xlabel('Sample Index');
set(xha,'FontSize',8);
set(gcf,'Name','compressed CP42 data file','NumberTitle','Off');

% END OF PLOTCP42STO

function y=dec2bin(x,y);
% DEC2BIN	Conversion from decimal to binary numbers
%
%		Y=DEC2BIN(X,Y)		: single number, recursive
% or	Y=DEC2BIN(X) 		: Row Vector, Single number
%
%		The binary representation is contained in the vector y.
%       The 1st element is the LSB.
%
%		When DEC2BIN is used only with one input argument
%       Y=DEC2BIN(X), the conversion is done iteratively.
%       X may be a row vector. Y is a matrix. Each column
%       represents a binary number.
%       The last row contains the LSB's.

if nargin==2			% Single Value - rec.
	if x>=1
	N=nextpow2(x+1E-4)-1;
	y(1,N+1)=1;
	x=x - (2^N);
		if abs(x)>=0.2
			y=dec2bin(x,y);
		end
	end
else					% Vector - non rec.
	maxx=max(max(x));
	minx=min(min(x));
		if minx<0
			disp('Range Error !');
			return;
		end
	basis=ceil(log(maxx+1E-4)/log(2));
	xk=x;
	y=[];
		for k=(basis-1):(-1):0
			xtmp=floor(xk ./(2^k));
			y=[y;xtmp];
			xk=xk - (2^k)*xtmp;
		end
end                     % Prozedurende "dec2bin"

function K=findgrup(a,modus);
% FINDGRUP
%		K=FINDGRUP(A,[MODUS])
%		Start and Endpositions of 1-Groups of the
%       binary vector A are determined.
%       If the optional argument MODUS (Default = 0)
%       is not used, then leading and finishing
%       1-Groups are neglected.
%
%		Example:
%		(1)	:	>> a=[1,1,1,0,0,0,1,1,0,0,0,1,1,1];
%				>> ki=findgrup(a) % or ki=findgrup(a,0)
%
%				ki=[7;
%				    8]
%
%		(2)	:	>> ki=findgrup(a,1)
%
%				ki=[ 1, 7,12;
%				     3, 8,14]

if nargin==1, modus=0; end

a=a(:).';

if modus~=0, a=[0,a,0]; end

if isempty(find(a)) | isempty(find(~a)),
	K=[];
	return;
end

% Where starts and ends Group of 1 ?
d=diff(a);

kmin=find(d>0) + 1; kmin=kmin(:);
kmax=find(d<0) + 1; kmax=kmax-1; kmax=kmax(:);

if a(1)==1, kmax=kmax(2:length(kmax)); end
if a(length(a))==1, kmin=kmin(1:length(kmin)-1); end
K=[kmin,kmax].';
if modus~=0, K=K-1; end			% END OF FINDGRUP

author: Sable

repo: mcbench-benchmarks


text:
function [M] = ft_connectivity_mim(inputdata, varargin)

% FT_CONNECTIVITY_MIM computes the multivariate interaction measure from a
% data-matrix containing the cross-spectral density. This implements the method
% described in Ewald et al., Estimating true brain connectivity from EEG/MEG data
% invariant to linear and static trasformations in sensor space. Neuroimage, 2012;
% 476:488.
%
% Use as
%   [m] = hcp_connectivity_mim(inputdata, ...)
%
% The input data should be an array organized as
%   Channel x Channel x Frequency
%
% The output m contains the newChannel x newChannel x Frequency connectivity measure,
% with newChannel equal to max(indices).
%
% Additional optional input arguments come as key-value pairs:
%   'indices' = 1xN vector with indices of the groups to which the channels belong,
%               e.g. [1 1 2 2] for a 2-by-2 connectivity between 2 planar MEG channels.
%
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS

% Copyright (C) 2011-2014 by the Human Connectome Project, WU-Minn Consortium (1U54MH091657)
% Copyright (C) 2021 Jan-Mathijs Schoffelen, DCCN
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

indices = ft_getopt(varargin, 'indices');

if isempty(indices) && isequal(size(inputdata(:,:,1)), [2 2])
  % simply assume two channels
  indices = [1 1 2 2];
end

sizein  = size(inputdata);
sizeout = [sizein 1];
sizeout(1:2) = max(indices);

% compute the inverse of the auto terms only once for speed up
invC = cell(sizeout(1),sizeout(3));
for kk = 1:sizeout(3)
  for k = 1:sizeout(1)
    invC{k,kk} = pinv(real(inputdata(indices==k,indices==k,kk)));
  end
end

M = zeros(sizeout);
for kk = 1:sizeout(3)
  for k = 1:sizeout(1)
    for m = 1:sizeout(1)
      indx1 = indices==k;
      indx2 = indices==m;
      %cs_aa_re = real(input(indx1,indx1));
      %cs_bb_re = real(input(indx2,indx2));
      cs_ab_im = imag(inputdata(indx1,indx2,kk));

      %inv_cs_bb_re = pinv(cs_bb_re);
      %inv_cs_aa_re = pinv(cs_aa_re);
      inv_cs_bb_re = invC{m,kk};
      inv_cs_aa_re = invC{k,kk};
      transp_cs_ab_im = transpose(cs_ab_im);
      M(k,m,kk) = trace(inv_cs_aa_re*cs_ab_im*inv_cs_bb_re*transp_cs_ab_im); % try to speed up by dividing calculation in steps
    end
  end
end

% taking the mldivide and mrdivide operators doesn't change the results, but speeds up by a factor of 4 over 1000 iterations (on LM Notebook)
% m = trace(cs_aa_re\cs_ab_im*inv_cs_bb_re*transp_cs_ab_im);

% % % block_a=cs_aa_re\cs_ab_im;
% % % block_b=cs_bb_re\transpose(cs_ab_im);
% % % m = trace(block_a*block_b);

% m = trace(cs_aa_re\cs_ab_im*(cs_bb_re\transpose(cs_ab_im)));


text:
% --------------------------------------------------------
% MDP Tracking
% Copyright (c) 2015 CVGL Stanford
% Licensed under The MIT License [see LICENSE for details]
% Written by Yu Xiang
% --------------------------------------------------------
%
% initialization for testing
function tracker = MDP_initialize_test(tracker, image_width, image_height, dres_det, is_show)

% normalization factor for features
tracker.image_width = image_width;
tracker.image_height = image_height;
tracker.max_width = max(dres_det.w);
tracker.max_height = max(dres_det.h);
tracker.max_score = max(dres_det.r);

tracker.streak_tracked = 0;
tracker.is_show = is_show;

author: yuxng


text:
function [translatedAbundances,normalizedAbundances,unmappedRows]=translateMetagenome2AGORA(MetagenomeAbundancePath,sequencingDepth,reconstructionResource)
% Translates organism identifiers in a published metagenomic or 16S rRNA
% data file with organism abundances (retrieved e.g., from  MetaPhlAn) to
% AGORA pan-model IDs. This will not catch every case since the format of
% input files with abundance data greatly varies. Feel free to modify this
% function and submit a pull request to enable more input files to be
% translated to AGORA. Moreover, slight spelling variations in taxa across
% input files may  result to taxa not being mapped. Check the unmappedRows
% output to identify these cases and modify the function accordingly.
% Pan-models that can be used to create microbiome models in mgPipe can be
% created with the function createPanModels.
% Consider running the function updateTaxonomyInfoAGORA to retrieve the
% most recent taxonomic assignment for the taxa to map. You might need to
% update your input file based on updated taxonomic classifications.
%
% USAGE:
%
%   [translatedAbundances,normalizedAbundances,unmappedRows]=translateMetagenome2AGORA(MetagenomeAbundancePath,sequencingDepth,reconstructionResource)
%
% INPUT:
%   MetagenomeAbundancePath   String containing the path to csv file with
%                             organism abundance data retrieved from
%                             16S rRNA or metagenomic samples (example:
%                             'SRP065497_taxonomy_abundances_v3.0.tsv').
%
% OPTIONAL INPUTS:
%   sequencingDepth           Sequencing depth on the taxonomical level
%                             in the input data (e.g., genus, species).
%                             Allowed inputs are 'Species','Genus',
%                             'Family','Order', 'Class', 'Phylum'.
%                             Default: 'Species'.
%  reconstructionResource     Name of the reconstruction resource to map
%                             the abundances to. Allowed inputs are 'AGORA',
%                             'AGORA2'. Default: 'AGORA'
%
% OUTPUTS:
%   translatedAbundances      Abundances with organism names from the
%                             input file translated to AGORA pan-model IDs
%   normalizedAbundances      Translated abundances normalized so they sum
%                             up to 1 for each sample
%   unmappedRows              Taxa on the selected taxonomical level that
%                             could not be mapped to AGORA pan-models
%
% .. Author: Almut Heinken, 02/2019
%                           10/2019: edited mapping based on an output file
%                           of the NG-Tax pipeline (Ramiro-Garcia et al.,
%                           F1000Res 2016). Made input more intuitive.
%                           11/2021: adapted based on updated taxonomy,
%                           added option to map to either AGORA or AGORA2

if nargin <2
    sequencingDepth = 'Species';
end

if nargin <3
    reconstructionResource = 'AGORA';
end

% read the csv file with the abundance data
metagenome_abundance = readInputTableForPipeline(MetagenomeAbundancePath);
metagenome_abundance{1,1}='';

% list the taxonomical levels
taxLevels={
    'Phylum'
    'Class'
    'Order'
    'Family'
    'Genus'
    'Species'
    };

if strcmp(reconstructionResource,'AGORA')
    infoFile = readInputTableForPipeline('AGORA_infoFile.xlsx');
elseif strcmp(reconstructionResource,'AGORA2')
    infoFile = readInputTableForPipeline('AGORA2_infoFile.xlsx');
else
    error('Incorrect input for reconstruction resource!')
end

% List all taxa in the AGORA resource
agoraTaxa={};
for i=1:length(taxLevels)
    % gets all pan-model IDs
    findTaxCol = find(strcmp(infoFile(1, :), taxLevels{i}));
    allTaxa = unique(infoFile(2:end, findTaxCol));
    allTaxa(:,1)=strrep(allTaxa(:,1),' ','_');
    allTaxa(:,1)=strcat('pan',allTaxa(:,1));
    agoraTaxa=vertcat(agoraTaxa,allTaxa);
end
allTaxa = unique(infoFile(2:end, find(strcmp(infoFile(1, :), 'ModelAGORA'))));
agoraTaxa=vertcat(allTaxa,agoraTaxa);

% format the input file so it can be matched
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),';','|');

% match different phrasing
% This may not catch all cases! Feel free to add additional cases of
% mismatching nomenclature and submit a pull request.
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Candidatus_','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_Family_XI_Incertae_Sedis',' Incertae Sedis XI');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_Family_XIII_Incertae_Sedis',' Incertae Sedis XIII');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'typhimurium','enterica');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_1','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_2','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_3','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_4','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_5','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_6','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_7','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_8','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_9','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),' family','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),' order','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'[','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),']','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),' ','_');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'.','_');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'-','_');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'__','_');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'___','_');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Ruminococcus_gauvreauii_group','Ruminococcus');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Ruminococcus_torques_group','Mediterraneibacter');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Ruminococcus_gnavus_group','Mediterraneibacter');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_coprostanoligenes_group','Eubacterium');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_eligens_group','Lachnospira');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_hallii_group','Anaerobutyricum');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_ruminantium_group','Eubacterium');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_ventriosum_group','Eubacterium');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Eubacterium_xylanophilum_group','Eubacterium');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'Clostridium_innocuum_group','Erysipelatoclostridium');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_sensu_stricto_1','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_sensu_stricto','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_sensu_stricto','');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_cardiffensis','g_Schaalia|s_Actinomyces_cardiffensis');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_georgiae','g_Schaalia|s_Actinomyces_georgiae');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_meyeri','g_Schaalia|s_Actinomyces_meyeri');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_odontolyticus','g_Schaalia|s_Actinomyces_odontolytica');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_odontolyticus','g_Schaalia|s_Actinomyces_odontolytica');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Actinomyces|s_Actinomyces_odontolyticus','g_Schaalia|s_Actinomyces_odontolytica');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Absiella|s_Absiella_dolichum','g_Amedibacillus|s_Absiella_dolichus');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Eubacterium|s_Eubacterium_hallii','g_Anaerobutyricum|s_Eubacterium_hallii');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Clostridium|s_citroniae','g_Lachnoclostridium|s_citroniae');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Clostridium|s_difficile','g_Clostridioides|s_difficile');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Clostridium|s_sordellii','g_Paeniclostridium|s_sordellii');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Clostridium|s_hathewayi','g_Hungatella|s_hathewayi');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Eubacterium|s_biforme','g_Holdemanella|s_biformis');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Eubacterium|s_dolichum','g_Amedibacillus|s_dolichus');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Propionibacterium|s_acnes','g_Cutibacterium|s_acnes');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Propionibacterium|s_granulosum','g_Cutibacterium|s_granulosum');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Ruminococcus|s_gnavus','g_Mediterraneibacter|s_gnavus');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'g_Ruminococcus|s_torques','g_Mediterraneibacter|s_torques');
% remove all rows that are not of the desired sequencing depth and edit the
% taxon descriptions to enable mapping to AGORA
% specific for MetaPhlAn input files
sLevels={
    'Species'   's_'
    'Genus'     'g_'
    'Family'    'f_'
    'Order'     'o_'
    'Class' 	'c_'
    'Phylum'    'p_'
    };
mpsDepth=sLevels{find(strcmp(sLevels(:,1),sequencingDepth)),2};
delRows=[];
cnt=1;
for i=2:size(metagenome_abundance,1)
    findSDepth=strsplit(metagenome_abundance{i,1},'|');
    % remove taxon entries without any mapping
    delArray=[];
    delCnt=1;
    for j=1:length(findSDepth)
        if any(strcmp(findSDepth{j},{'s_','g_','f_','o_','c_'}))
            delArray(delCnt)=j;
            delCnt=delCnt+1;
        end
    end
    findSDepth(delArray)=[];
    if length(findSDepth)>1
            % only get the sequencing level you are looking for
            if strncmp(findSDepth{1,end},mpsDepth,2) && length(findSDepth{1,end})>2
                % if the genus information is missing on the species level
                if strncmp(findSDepth{1,end},'s_',2)
                    sname=regexprep(findSDepth{1,end},mpsDepth,'','once');
                    % if species name also contains genus
                    if contains(sname,'_')
                        sname=strsplit(sname,'_');
                        sname=sname{1,2};
                    end
                    gname=regexprep(findSDepth{1,end-1},'g_','','once');
                    if ~strncmp(gname,sname,length(gname))
                        metagenome_abundance{i,1}=strcat(gname,'_',sname);
                    end
                else
                    metagenome_abundance{i,1}=findSDepth{1,end};
                    metagenome_abundance{i,1}=regexprep(metagenome_abundance{i,1},mpsDepth,'','once');
                end
            else
                delRows(cnt,1)=i;
                cnt=cnt+1;
            end
    end
end
metagenome_abundance(delRows,:)=[];

% Remove unclassified organisms
agoraTaxa(strncmp('unclassified', agoraTaxa, 12)) = [];
metagenome_abundance(~cellfun(@isempty, strfind(metagenome_abundance(:,1),'_uncl')),:) = [];

% replace unknown species with pan-genus model
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),'_cf','');

% summarize duplicate entries
[uniqueA,i,j] = unique(metagenome_abundance(:,1));
n  = accumarray(j(:),1);
Dupes=uniqueA(find(n>1));
delArray=[];
cnt=1;
for i=1:length(Dupes)
    indexToDupes = find(strcmp(metagenome_abundance(:,1),Dupes{i}));
    for j=2:length(indexToDupes)
        for k=2:size(metagenome_abundance,2)
            if contains(version,'(R202') % for Matlab R2020a and newer
                metagenome_abundance{indexToDupes(1),k}=metagenome_abundance{indexToDupes(1),k}+metagenome_abundance{indexToDupes(j),k};
            else
                metagenome_abundance{indexToDupes(1),k}=num2str(str2double(metagenome_abundance{indexToDupes(1),k})+str2double(metagenome_abundance{indexToDupes(j),k}));
            end
        end
        delArray(cnt,1)=indexToDupes(j);
        cnt=cnt+1;
    end
end
metagenome_abundance(delArray,:)=[];

% match and find overlapping IDs
% metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),mpsDepth,'');
metagenome_abundance(:,1)=strrep(metagenome_abundance(:,1),' ','_');
metagenome_abundance(2:end,1)=strcat('pan',metagenome_abundance(2:end,1));

% remove and report taxa in abundance data table that could not be found in
% AGORA and/or are not of the desired sequencing depth

findTaxCol = find(strcmp(infoFile(1, :), sequencingDepth));
allTaxa = unique(infoFile(2:end, findTaxCol));
allTaxa(:,1)=strrep(allTaxa(:,1),' ','_');
allTaxa(:,1)=strcat('pan',allTaxa(:,1));
[unmappedRows,IA] = setdiff(metagenome_abundance(:,1),allTaxa);
IA(1,:)=[];
metagenome_abundance(IA,:)=[];
translatedAbundances=metagenome_abundance;

% normalize the abundances so that sum for each individual is 1
normalizedAbundances=translatedAbundances;
for i=2:size(translatedAbundances,2)
    if contains(version,'(R202') % for Matlab R2020a and newer
        if sum(cell2mat(translatedAbundances(2:end,i)))>0
            for j=2:size(translatedAbundances,1)
                normalizedAbundances{j,i}=cell2mat(translatedAbundances(j,i))/sum(cell2mat(translatedAbundances(2:end,i)));
            end
        else
            for j=2:size(translatedAbundances,1)
                normalizedAbundances{j,i}=translatedAbundances{j,i};
            end
        end
    else
    if sum(str2double(translatedAbundances(2:end,i)))>0
        for j=2:size(translatedAbundances,1)
            normalizedAbundances{j,i}=num2str(str2double(translatedAbundances(j,i))/sum(str2double(translatedAbundances(2:end,i))));
        end
    else
        for j=2:size(translatedAbundances,1)
            normalizedAbundances{j,i}=translatedAbundances{j,i};
        end
    end
    end
end

% export the translated abundances
writetable(cell2table(translatedAbundances),'translatedAbundances','FileType','text','WriteVariableNames',false,'Delimiter','tab');

end

author: opencobra

repo: cobratoolbox