text:
function p03_header ( )

%*****************************************************************************80
%
%% P03_HEADER prints some information about problem 03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    None.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'P03:\n' );
  fprintf ( 1, '  Strang and Persson example #3\n' );
  fprintf ( 1, '  The unit square, with a hole.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The hole is a concentric circle of radius 0.4.\n' );
  fprintf ( 1, '  A uniform mesh density is requested.\n' );
  fprintf ( 1, '  Element sizes tried were 0.4, 0.2, 0.1.\n' );

  fprintf ( 1, '\n' );
  boundary_segment_num = p03_boundary_segment_num ( );
  fprintf ( 1, '  Number of boundary segments = %d\n', boundary_segment_num );
  fixed_num = p03_fixed_num ( );
  fprintf ( 1, '  Number of fixed points = %d\n', fixed_num );
  hole_num = p03_hole_num ( );
  fprintf ( 1, '  Number of holes = %d\n', hole_num );

  return
end


text:
function f = spp_scale_features(f, feat_norm_mean)
% My initial experiments were conducted on features with an average norm
% very close to 20. Using those features, I determined a good range of SVM
% C values to cross-validate over. Features from different layers end up
% have very different norms. We rescale all features to have an average norm
% of 20 (why 20? simply so that I can use the range of C values found in my
% initial experiments), to make the same search range for C reasonable
% regardless of whether these are pool5, fc6, or fc7 features. This strategy
% seems to work well. In practice, the optimal value for C ends up being the
% same across all features.
target_norm = 20;
f = f .* (target_norm / feat_norm_mean);

author: ShaoqingRen


function [ber,rate] = Mrate_method()
% traditional SVD algorithm for rate maximization

global  H Ns V_ropt W_ropt;
[U,~,V] = svd(H);
V_ropt = V(:,1:Ns);
%power constraint
V_ropt = V_ropt / norm(V_ropt,'fro');
W_ropt = U(:,1:Ns);

ber = get_ber(V_ropt,W_ropt);
rate = get_rate(V_ropt,W_ropt);


text:
function test_ft_sourceparcellate

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_sourceparcellate ft_read_atlas ft_sourceinterpolate

% see also https://github.com/fieldtrip/fieldtrip/issues/1753

%%

atlasfilename   = dccnpath('/home/common/matlab/fieldtrip/template/atlas/aal/ROI_MNI_V4.nii');

npos = 38*48*41;
ntime = 20;

source = [];
source.pos = 200*rand(npos,3) - 100;
source.unit = 'mm';
source.time = 1:ntime;
source.pow = randn(npos,ntime);
source.inside = 1:npos;

atlas = ft_read_atlas(atlasfilename);

cfg = [];
cfg.interpmethod = 'nearest';
cfg.parameter = 'tissue';
parcellation = ft_sourceinterpolate(cfg, atlas, source);

cfg = [];
cfg.parcellation = 'tissue';
result = ft_sourceparcellate(cfg, source, parcellation);

%%

x = -5:5;
y = -5:5;
z = -5:5;

[X, Y, Z] = ndgrid(x, y, z);
dim = [length(x) length(y) length(z)];

source = [];
source.dim = dim;
source.pos = [X(:) Y(:) Z(:)];
source.time = linspace(-0.1, 0.5, 15);
source.pow   = randn(prod(dim), 1);
for i=1:prod(dim)
  source.eta{i} = source.pow(i);
end
source.timecourse = randn(prod(dim), 15);
source.moment     = cell(prod(dim), 1);
for i=1:prod(dim)
  source.moment{i} = source.timecourse(i,:);
end
source.powdimord = 'pos';
source.etadimord = '{pos}';
source.timecoursedimord = 'pos_time';
source.momentdimord = '{pos}_ori_time';

atlas = [];
atlas.dim = dim;
atlas.pos = [X(:) Y(:) Z(:)];
atlas.tissue = zeros(prod(dim),1);
roi1 = atlas.pos(:,1)<0;
roi2 = atlas.pos(:,1)>0;
atlas.tissue(roi1) = 1;
atlas.tissue(roi2) = 2;
atlas.tissuelabel{1} = 'left';
atlas.tissuelabel{2} = 'right';

%%

cfg = [];
cfg.method = 'mean';

cfg.parameter = 'pow';
parcel1 = ft_sourceparcellate(cfg, source, atlas);
manual = mean(source.(cfg.parameter)(roi1, :), 1);
assert(isequal(parcel1.pow(1,:), manual));

cfg.parameter = 'eta';
parcel2 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel1.pow(:), parcel2.eta(:), 'abstol', 100*eps));

cfg.parameter = 'timecourse';
parcel3 = ft_sourceparcellate(cfg, source, atlas);
manual = mean(source.(cfg.parameter)(roi1, :), 1);
assert(isequal(parcel3.timecourse(1,:), manual));

cfg.parameter = 'moment';
parcel4 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel3.timecourse(:), parcel4.moment(:), 'abstol', 100*eps));

%%

cfg = [];
cfg.method = 'median';

cfg.parameter = 'pow';
parcel1 = ft_sourceparcellate(cfg, source, atlas);
manual = median(source.(cfg.parameter)(roi1, :), 1);
assert(isequal(parcel1.pow(1,:), manual));

cfg.parameter = 'eta';
parcel2 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel1.pow(:), parcel2.eta(:), 'abstol', 100*eps));

cfg.parameter = 'timecourse';
parcel3 = ft_sourceparcellate(cfg, source, atlas);
manual = median(source.(cfg.parameter)(roi1, :), 1);
assert(isequal(parcel3.timecourse(1,:), manual));

cfg.parameter = 'moment';
parcel4 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel3.timecourse(:), parcel4.moment(:), 'abstol', 100*eps));

%%

cfg = [];
cfg.method = 'min';

cfg.parameter = 'pow';
parcel1 = ft_sourceparcellate(cfg, source, atlas);
manual = min(source.(cfg.parameter)(roi1, :), [], 1);
assert(isequal(parcel1.pow(1,:), manual));

cfg.parameter = 'eta';
parcel2 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel1.pow(:), parcel2.eta(:), 'abstol', 100*eps));

cfg.parameter = 'timecourse';
parcel3 = ft_sourceparcellate(cfg, source, atlas);
manual = min(source.(cfg.parameter)(roi1, :), [], 1);
assert(isequal(parcel3.timecourse(1,:), manual));

cfg.parameter = 'moment';
parcel4 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel3.timecourse(:), parcel4.moment(:), 'abstol', 100*eps));

%%

cfg = [];
cfg.method = 'max';

cfg.parameter = 'pow';
parcel1 = ft_sourceparcellate(cfg, source, atlas);
manual = max(source.(cfg.parameter)(roi1, :), [], 1);
assert(isequal(parcel1.pow(1,:), manual));

cfg.parameter = 'eta';
parcel2 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel1.pow(:), parcel2.eta(:), 'abstol', 100*eps));

cfg.parameter = 'timecourse';
parcel3 = ft_sourceparcellate(cfg, source, atlas);
manual = max(source.(cfg.parameter)(roi1, :), [], 1);
assert(isequal(parcel3.timecourse(1,:), manual));

cfg.parameter = 'moment';
parcel4 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel3.timecourse(:), parcel4.moment(:), 'abstol', 100*eps));

%%

cfg = [];
cfg.method = 'std';

cfg.parameter = 'pow';
parcel1 = ft_sourceparcellate(cfg, source, atlas);
manual = std(source.(cfg.parameter)(roi1, :), 0, 1);
assert(isequal(parcel1.pow(1,:), manual));

cfg.parameter = 'eta';
parcel2 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel1.pow(:), parcel2.eta(:), 'abstol', 100*eps));

cfg.parameter = 'timecourse';
parcel3 = ft_sourceparcellate(cfg, source, atlas);
manual = std(source.(cfg.parameter)(roi1, :), 0, 1);
assert(isequal(parcel3.timecourse(1,:), manual));

cfg.parameter = 'moment';
parcel4 = ft_sourceparcellate(cfg, source, atlas);
assert(isalmostequal(parcel3.timecourse(:), parcel4.moment(:), 'abstol', 100*eps));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s=cosmo_structjoin(varargin)
% joins values in structs or key-value pairs
%
% s=cosmo_structjoin(arg1, arg2, ...)
%
% Inputs:
%   arg{X}      Any of the following:
%               - if a cell, then it should contain structs or key-value
%                 pairs
%               - if a struct, then value=arg{X}.(key) for each key in
%                 fieldnames(arg{X}) is stored as s.(key)=value.
%                 In this case, it is required that the struct is of size
%                 1x1.
%               - if a string, then it is stored as s.(arg{X})=arg{X+1}
%
% Returns:
%   s           Struct with fieldnames and their associated values
%               as key-value pairs. Values in the input are stored from
%               left-to-right.
%
% Example:
%     x=cosmo_structjoin('a',{1;2},'b',{3,4});
%     cosmo_disp(x);
%     %|| .a
%     %||   { [ 1 ]
%     %||     [ 2 ] }
%     %|| .b
%     %||   { [ 3 ]  [ 4 ] }
%     y=cosmo_structjoin(x,'a',66,'x',x,{'c','hello'});
%     cosmo_disp(y);
%     %|| .a
%     %||   [ 66 ]
%     %|| .b
%     %||   { [ 3 ]  [ 4 ] }
%     %|| .x
%     %||   .a
%     %||     { [ 1 ]
%     %||       [ 2 ] }
%     %||   .b
%     %||     { [ 3 ]  [ 4 ] }
%     %|| .c
%     %||   'hello'
%
%     % simulate a function definition function out=f(varargin)
%     % to illustrate overriding default values
%     varargin={'radius',2,'nsamples',12};
%     defaults=struct();
%     defaults.radius=10;
%     defaults.nfeatures=2;
%     params=cosmo_structjoin(defaults,varargin);
%     cosmo_disp(params);
%     %|| .radius
%     %||   [ 2 ]
%     %|| .nfeatures
%     %||   [ 2 ]
%     %|| .nsamples
%     %||   [ 12 ]
%
%     % illustrate overriding values in 'sub-structs' (structs as values in
%     % other structs)
%     v=struct();
%     v.a.foo={1,2};
%     v.b=3;
%     w=struct();
%     w.a.bar=[2,3];
%     w.c=4;
%     j=cosmo_structjoin(v,w);
%     cosmo_disp(j)
%     %|| .a
%     %||   .foo
%     %||     { [ 1 ]  [ 2 ] }
%     %||   .bar
%     %||     [ 2         3 ]
%     %|| .b
%     %||   [ 3 ]
%     %|| .c
%     %||   [ 4 ]
%
%
% Notes:
%  - this function can be used to parse input arguments (including
%    varargin in a "'key1',value1,'key2',value2,..." fashion)
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    % use a wrapper so that recursive calls can be forwarded to the wrapper
    % without using mfilename
    s=structjoin_wrapper(varargin{:});

function s=structjoin_wrapper(varargin)

    s=struct(); % output
    n=numel(varargin);

    k=0;
    while k<n
        % go over all input arguments
        k=k+1;
        v=varargin{k}; %k-th argument

        if iscell(v)
            if isempty(v)
                continue;
            end

            % use recursion
            v=structjoin_wrapper(v{:});
        end

        if isstruct(v)
            if numel(v)~=1
                error(['only singleton structs (of size 1x1) '...
                            'are supported']);
            end

            % overwrite any values in s
            fns=fieldnames(v);

            for j=1:numel(fns);
                fn=fns{j};
                v_fn=v.(fn);

                if isstruct(v_fn) && isfield(s,fn) && isstruct(s.(fn))
                    s.(fn)=structjoin_wrapper(s.(fn),v_fn);
                else
                    s.(fn)=v_fn;
                end
            end
        elseif ischar(v)
            % <key>, <value> pair
            if k+1>n
                % cannot be last argument
                error('Missing argument after key ''%s''', v);
            end

            % move forward to next argument and get value
            k=k+1;
            vv=varargin{k};

            if isstruct(vv) && isfield(s,v) && isstruct(s.(v))
                s.(v)=structjoin_wrapper(s.(v),vv);
            else
                s.(v)=vv;
            end
        else
            error(['Illegal input at position %d: expected cell, struct, ',...
                        'or string'], k);
        end
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% The COBRAToolbox: testConvertCobraToSBML.m
%
% Purpose:
%     - test the convertCobraToSBML function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testConvertCobraToSBML'));
cd(fileDir);

% test variables
sbmlLevel = 0;
sbmlVersion = 0;
compSymbolList = '';
compNameList = '';
model = getDistributedModel('ecoli_core_model.mat');

% function outputs
warning('off', 'all')
  sbmlModel = convertCobraToSBML(model, sbmlLevel, sbmlVersion, compSymbolList, compNameList);
  sbmlModel_2 = convertCobraToSBML(model, sbmlLevel, sbmlVersion);
  assert(length(lastwarn()) > 0)
warning('on', 'all')
% test
assert(isequaln(sbmlModel, sbmlModel_2));

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function [common,system_not_matlab,matlab_not_system] = comparePaths(NAME,printLevel)
%compare the system and matlab paths for a given NAME environment variable
%
% INPUT
% NAME          environment variable, e.g., PATH or LD_LIBRARY_PATH
% printLevel
%
% OUTPUT
% common                cell array
% system_not_matlab
% matlab_not_system
%
% USAGE example
% [common,system_not_matlab,matlab_not_system] = comparePaths('LD_LIBRARY_PATH',1);

if ~exist('printLevel','var')
    printLevel=0;
end

if strcmp(NAME,'PYTHON')
    matlab_split = get_py_path();
    [success,system_response]=system('python -c "import sys; print(sys.path)"');
    system_split = split(system_response,',');
    system_split{1}=strtrim(system_split{1});
    system_split{1}=strrep(system_split{1},'[','');
    system_split{end}=strtrim(system_split{end});
    system_split{end}=strrep(system_split{end},']','');
    system_split{end}=strrep(system_split{end},'''','');
    for i=1:length(system_split)
        system_split{i}=strrep(system_split{i},' ','');
        system_split{i}=strrep(system_split{i},'''','');
    end
else
    [success,system_response]=system(['echo $' NAME]);
    system_split = split(system_response,':');
    for i=1:length(system_split)
        system_split{i} = strtrim(system_split{i});
    end

    matlab_response = getenv(NAME);
    matlab_split = split(matlab_response,':');
    for i=1:length(matlab_split)
        matlab_split{i} = strtrim(matlab_split{i});
    end
end

if printLevel>1
    disp(system_split)
    disp(matlab_split)
end

common = intersect(system_split,matlab_split);
if printLevel>1
    fprintf('%s\n',[NAME ' common to system and matlab'])
    for i=1:length(common)
        fprintf('%s\n',common{i})
    end
end

system_not_matlab = setdiff(system_split,matlab_split);
if printLevel>0
    if isempty(system_not_matlab)
        fprintf('%s\n',['No ' NAME ' in system but not matlab'])
    else
        fprintf('\n%s\n',[NAME ' in system but not matlab'])
        for i=1:length(system_not_matlab)
            fprintf('%s\n',system_not_matlab{i})
        end
    end
end

matlab_not_system = setdiff(matlab_split,system_split);
if printLevel>0
    if isempty(matlab_not_system)
        fprintf('%s\n',['No ' NAME ' in matlab but not system'])
    else
        fprintf('\n%s\n',[NAME ' in matlab but not system'])
        for i=1:length(matlab_not_system)
            fprintf('%s\n',matlab_not_system{i})
        end
    end
end


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

% Last Modified by GUIDE v2.5 06-Aug-2020 15:26:37

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

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

    centerUIFigure(hObject);

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

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

    populateGUI(handles, action, lv);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, action, lv)
    [powerSrcListStr, powerSrcs] = lv.getPowerSrcsListBoxStr();
    set(handles.refPwrCompCombo,'String',powerSrcListStr);

    if(not(isempty(action.src)))
        ind = find(powerSrcs == action.src,1,'first');
    else
        ind = [];
    end

    if(not(isempty(ind)))
        set(handles.refPwrCompCombo,'Value',ind);
    end

    if(action.activeStateToSet)
        set(handles.stateCombo,'Value',1);
    else
        set(handles.stateCombo,'Value',2);
    end

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditActionSetPwrSrcStateGUI_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
        action = getappdata(hObject,'action');
        lv = getappdata(hObject,'lv');

        [~, powerSrcs] = lv.getPowerSrcsListBoxStr();
        ind = get(handles.refPwrCompCombo,'Value');
        action.src = powerSrcs(ind);

        contents = cellstr(get(handles.stateCombo,'String'));
        stateStr = contents{get(handles.stateCombo,'Value')};

        switch stateStr
            case 'Active'
                state = true;
            case 'Inactive'
                state = false;
            otherwise
                error('Invalid state string: %s', stateStr);
        end

        action.activeStateToSet = state;

        varargout{1} = true;
        close(handles.lvd_EditActionSetPwrSrcStateGUI);
    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)
    uiresume(handles.lvd_EditActionSetPwrSrcStateGUI);

% --- 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_EditActionSetPwrSrcStateGUI);

% --- Executes on selection change in refPwrCompCombo.
function refPwrCompCombo_Callback(hObject, eventdata, handles)
% hObject    handle to refPwrCompCombo (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 refPwrCompCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from refPwrCompCombo

% --- Executes during object creation, after setting all properties.
function refPwrCompCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to refPwrCompCombo (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 selection change in stateCombo.
function stateCombo_Callback(hObject, eventdata, handles)
% hObject    handle to stateCombo (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 stateCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from stateCombo

% --- Executes during object creation, after setting all properties.
function stateCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stateCombo (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 key press with focus on lvd_EditActionSetPwrSrcStateGUI or any of its controls.
function lvd_EditActionSetPwrSrcStateGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditActionSetPwrSrcStateGUI (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_EditActionSetPwrSrcStateGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [Ain, Cin, bin, fin, center] = initialize_components(Y, K, tau, options, P)

% Initalize components using a greedy approach followed by hierarchical
% alternative least squares (HALS) NMF. Optional use of spatio-temporal
% downsampling to boost speed.

%Input:
%Y          d1 x d2 x T movie, raw data
%K          number of neurons to extract (default value: 30)
%tau        standard deviation of neuron size (default value: 5)

%options    fine-tuning parameters (optional)
%           options.init_method: method of initialization ('greedy','sparse_NMF','HALS')
%           options.
%           options.nIter: number of iterations for shape tuning (default 5)
%           options.gSiz: size of kernel (default 2*tau + 1)
%           options.ssub: spatial downsampling factor (default 1)
%           options.tsub: temporal downsampling factor (default 1)
%           options.nb: rank of background component (default 1)
%           options.save_memory: flag for processing data in chunks to save memory (default 0)
%           options.windowSiz: size of spatial window when computing the median (default 32 x 32)
%           options.chunkSiz: number of timesteps to be processed simultaneously if on save_memory mode (default: 100)
%           options.med_app: number of timesteps to be interleaved for fast (approximate) median calculation (default: 1, no approximation)
%           options.rem_prct: percentile to be removed before initialization (default: 20)

% P         parameter struct used for normalization by noise and user feed component centroids (optional)

%
%Output:
%Ain        (d1*d2) x K matrix, location of each neuron
%Cin        T x K matrix, calcium activity of each neuron
%center     K x 2 matrix, inferred center of each neuron
%bin        (d1*d2) X nb matrix, initialization of spatial background
%fin        nb X T matrix, initalization of temporal background
%res        d1 x d2 x T movie, residual
%
%Authors: Eftychios A. Pnevmatikakis and Pengchen Zhou, with inputs from Weijian Yang

defoptions = CNMFSetParms;
if nargin < 4 || isempty(options); options = defoptions; end
if nargin < 2 || isempty(K)
    K = 30;
    fprintf('Number of components to be detected not specified. Using the default value 30. \n');
end
if nargin < 3 || isempty(tau)
    options.gSig = 5;
    fprintf('Standard deviation for neuron not specified. Using the default value 5. \n');
else options.gSig = tau;
end

if ~isfield(options,'init_method'); options.init_method = 'greedy'; end
if ~isfield(options,'rem_prct') || isempty(options.rem_prct); options.rem_prct = defoptions.rem_prct; end
% downsample the data

if ~isfield(options, 'ssub'); options.ssub = 1; end; ssub = options.ssub;
if ssub == 1; fprintf('No spatial downsampling is performed. Consider spatial downsampling if the field of view is very large. \n'); end
if ~isfield(options, 'tsub'), options.tsub = 1; end; tsub = options.tsub;
if tsub == 1; fprintf('No temporal downsampling is performed. Consider temporal downsampling if the recording is very long. \n'); end

if ~isfield(options,'noise_norm') || isempty(options.noise_norm)
    options.noise_norm = defoptions.noise_norm; % normalization by noise (true if P is present)
end

if nargin < 5
    if options.noise_norm
        warning('Normalization by noise value is not performed since noise values are not provided. \n');
    end
    options.noise_norm = false;
end

ndimsY = ndims(Y)-1;
sY = size(Y);
d = sY(1:ndimsY);
T = sY(end);

if options.noise_norm
    mY = mean(Y,ndims(Y));
    norm_image = mY + median(mY(:)) + 1e-4;
    min_noise = norm_image;
    %min_noise = prctile(P.sn(P.sn>0),options.noise_norm_prctile);
    %Y = bsxfun(@times,Y,reshape(1./max(P.sn,min_noise),d));
    Y = bsxfun(@times,Y,reshape(1./double(min_noise),d));
end

ds = d;
ds(1:2) = ceil(d(1:2)/ssub); % do not subsample along z axis
%d1s = ceil(d1/ssub);        %size of downsampled image
%d2s = ceil(d2/ssub);
Ts = floor(T/tsub);         %reduced number of frames
% spatial downsampling
fprintf('starting resampling \n')
if ssub~=1;
    if ndimsY == 2; Y_ds = imresize(Y, [ds(1), ds(2)], 'box'); end
    if ndimsY == 3;
        Y_ds = zeros([ds(1:2),T,ds(end)]);
        for z = 1:ds(3)
            Y_ds(:,:,:,z) = imresize(squeeze(Y(:,:,z,:)), [ds(1), ds(2)], 'box');
        end
        Y_ds = permute(Y_ds,[1,2,4,3]);
    end
else
    Y_ds = Y;
end
% temporal downsampling
if tsub~=1
    if ndimsY == 2; Y_ds = squeeze(mean(reshape(Y_ds(:, :, 1:(Ts*tsub)),ds(1), ds(2), tsub, Ts), 3)); end
    if ndimsY == 3; Y_ds = squeeze(mean(reshape(Y_ds(:, :, :, 1:(Ts*tsub)),ds(1), ds(2), ds(3), tsub, Ts), 4)); end
end

options_ds = options;
options_ds.d1 = ds(1);
options_ds.d2 = ds(2);

if strcmpi(options.init_method,'greedy')
    % run greedy method
    if nargin < 5 || ~isfield(P,'ROI_list')
        ROI_list = [];
    else
        ROI_list = round(P.ROI_list/ssub);
        K = size(ROI_list,1);
    end
    fprintf('Initializing components with greedy method \n');
    [Ain, Cin, bin, fin] = greedyROI(Y_ds, K, options, ROI_list);
elseif strcmpi(options.init_method, 'greedy_corr')
    fprintf('Initializing components with greedy_corr method \n');
    [Ain, Cin, bin, fin] = greedyROI_corr(Y_ds, K, options);
elseif strcmpi(options.init_method,'sparse_NMF')
    % run sparse_NMF method
    fprintf('Initializing components with sparse NMF \n');
    [Ain,Cin,bin,fin] = sparse_NMF_initialization(Y_ds,K,options_ds);
elseif strcmpi(options.init_method,'HALS')
    fprintf('Initializing components with HALS \n');
    [Ain,Cin,bin,fin] = HALS_initialization(Y_ds,K,options_ds);
else
    error('Unknown initialization method')
end

% refine with HALS
fprintf('Refining initial estimates with HALS...');
[Ain, Cin, bin, fin] = HALS(Y_ds, full(Ain), Cin, bin, fin, options_ds);
fprintf('  done \n');
%% upsample Ain, Cin, bin, fin
if nargout == 5
    if ndimsY == 2; center = ssub*com(Ain,ds(1),ds(2)); else center = ssub*com(Ain,ds(1),ds(2),ds(3)); end
end

Ain = imresize(reshape(full(Ain), [ds(1),ds(2), size(Ain,2)*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]); %,prod(d)/d(1)/d(2)*sum(K)]);
Ain = max(sparse(reshape(Ain, prod(d), [])),0);

bin = imresize(reshape(bin,[ds(1),ds(2), options.nb*prod(ds)/ds(1)/ds(2)]),[d(1),d(2)]);
bin = max(double(reshape(bin,prod(d),[])),0);

if options.noise_norm
    %Ain = bsxfun(@times,Ain,double(max(P.sn(:),min_noise)));
    %bin = bsxfun(@times,bin,max(P.sn(:),min_noise));
    Ain = bsxfun(@times,Ain,double(min_noise(:)));
    bin = bsxfun(@times,bin,double(min_noise(:)));
end
Cin = max(imresize(Cin, [size(Cin, 1), Ts*tsub]),0);
fin = max(imresize(fin, [options.nb, Ts*tsub]),0);
if T ~= Ts*tsub
    Cin = padarray(Cin, [0, T-Ts*tsub], 'post');
    fin = padarray(fin, [0, T-Ts*tsub], fin(end), 'post');
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [A,xy]= cut_grid_data

% Generate a cut-grid graph for the ICM 2006 talk example
%
% The graph has 64 nodes and 95 edges
% A is an n x m incidence matrix (n is number of nodes, m is number of edges)
% xy is the location data
%
% Original code by Arpita Ghosh, modified for ICM06 talk by Almir Mutapcic

n  = 8;
r1 = 2;
y  = ones(n,1) * (1:n);
x  = y';
x  = x(:);
y  = y(:);
dx = x * ones(1,n^2);
dy = y * ones(1,n^2);
xy = [ x, y ];

% Find the adjacency matrix, manually deleting edges to get down to size
Adj1 = tril( ( dx - dx' ) .^ 2 + ( dy - dy' ) .^2 < r1, -1 );
Adj1(49,41) = 0;
Adj1(50,42) = 0;
Adj1(16,8)  = 0;
Adj1(24,16) = 0;
Adj1(15,7)  = 0;
Adj1(23,15) = 0;
Adj1(10,1)  = 0;
Adj1(21,13) = 0;
Adj1(13,5)  = 0;
Adj1(22,14) = 0;
Adj1(14,6)  = 0;
Adj1(51,43) = 0;
Adj1(52,44) = 0;
Adj1(53,45) = 0;
Adj1(54,46) = 0;
Adj1(42,41) = 0;
Adj1(34,33) = 0;
Adj1(26,25) = 0;

% Build the incidence matrix
[i,j] = find(Adj1);
m = length(i);
A = sparse( [i;j], [1:m,1:m]', [ones(m,1);-ones(m,1)] );

author: yu-jiang

repo: radpbook

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


text:
function c = tapas_gaussian_obs_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Gaussian noise observation model for continuous responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The Gaussian noise observation model assumes that responses have a Gaussian distribution around
% the inferred mean of the relevant state. The only parameter of the model is the noise variance
% (NOT standard deviation) zeta.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'tapas_gaussian_obs';

% Sufficient statistics of Gaussian parameter priors
%
% Zeta
c.logzemu = log(0.005);
c.logzesa = 0.1;

% Gather prior settings in vectors
c.priormus = [
    c.logzemu,...
         ];

c.priorsas = [
    c.logzesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_gaussian_obs;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_gaussian_obs_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
function [T, M, S, DISTR, df] = dtiDirTestStat(g1, g2, Y, mask)

% Computes voxel-wise test statistics for two groups from a data array of
% DTI directions.
%
%   [T, M, S, DISTR, df] = dtiDirTestStat(g1, g2, DIR_ARRAY, [MASK])
%
% Input:
%   g1, g2      List of indices that correspond to each group out of 1:N
%                   E.g: g1 = 1:7, g2 = 8:14, N = 14
%   DIR_ARRAY   Data array of size XxYxZx3xN (or nx3xN), where X, Y, Z are the volume
%                   dimensions and N is the number of subjects.
%                   The directions are contained in the dimension of size 3.
%   MASK        Optional XxYxZ binary array. Values of M and S are computed
%                   where mask = 1; in other voxels, M and S are set to NaN.
%                   Default is entire volume.
%
% Output:
%   T           XxYxZx1 (or nx1) array of test statistics (NaN where mask = 0)
%   M           XxYxZx3x2 (or nx3x2) array of mean directions for both groups (NaN where mask = 0)
%   S           XxYxZx1 (or nx1) array of dispersions (NaN where mask = 0)
%   DISTR       The string 'f'
%   df          The degrees of freedom of the f distribution
%
% Utilities:    ndfun.m, dtiEig.m, dti33to6.m
%
% WARNING: If using Pentium 4, eliminate NaN's from array before running
% (processor bug).
%
% Copyright by Armin Schwartzman, 2004

% HISTORY:
%   2004.06.23 ASH (armins@stanford.edu) wrote it.
%   2006.07.18 ASH (armins@stanford.edu) added indexed format capability.
%

% Check inputs
if (ndims(Y)==2 | ndims(Y)==3),
    Ind = 1;    % Data in indexed nx3xN format
    Y = shiftdim(Y, -2);
else
    Ind = 0;    % Data in XxYxZx3xN format
end
if (ndims(Y)<4 | ndims(Y)>5),
    error('Wrong input format');
end
if (~exist('mask')),
    mask = ones([size(Y,1) size(Y,2) size(Y,3)]);
end

% Computations
N1 = length(g1);
N2 = length(g2);
N  = N1 + N2;
p = size(Y, 4); % should be 3
DISTR = 'f';
df = [p-1, (p-1)*(N-2)];

Y = permute(Y, [4 6 5 1 2 3]); % permutation required by ndfun
Tbar1 = shiftdim(sum(ndfun('mult', Y(:,1,g1,:,:,:), permute(Y(:,1,g1,:,:,:), [2 1 3:6])), 3), 3) / N1;
Tbar2 = shiftdim(sum(ndfun('mult', Y(:,1,g2,:,:,:), permute(Y(:,1,g2,:,:,:), [2 1 3:6])), 3), 3) / N2;
Tbar = (N1 * Tbar1 + N2 * Tbar2)/(N1 + N2);
[vec1, val1] = dtiEig(dti33to6(Tbar1)); % ndfun('eig', Tbar1);
[vec2, val2] = dtiEig(dti33to6(Tbar2)); % ndfun('eig', Tbar2);
[vec, val] = dtiEig(dti33to6(Tbar)); % ndfun('eig', Tbar);
M = cat(5, vec1(:,:,:,:,1), vec2(:,:,:,:,1));
S = (N - N1*val1(:,:,:,1) - N2*val2(:,:,:,1)) / (N-2);;
T = (N1*val1(:,:,:,1) + N2*val2(:,:,:,1) - N*val(:,:,:,1)) ./ S;
M(~repmat(mask,[1 1 1 3 2])) = NaN;
S(~mask) = NaN;
T(~mask) = NaN;

% Adjust output
if Ind,
    T = shiftdim(T, 2);
    M = shiftdim(M, 2);
    S = shiftdim(S, 2);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = lnDiffCumGaussian(u, uprime)

% LNDIFFCUMGAUSSIAN Log of the difference between two cumulative Gaussians.
% FORMAT
% DESC computes the logarithm of the difference between two
% cumulative Gaussian distributions.
% ARG u1 : the argument of the first (positive) cumulative
% Gaussian.
% ARG u2 : the argument of the second (negative) cumulative
% Gaussian.
% RETURN f : where f = log(cumGaussian(u1) - cumGaussian(u2)).
%
% SEEALSO : cumGaussian, gaussOverDiffCumGaussian, lnCumGaussian
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% NDLUTIL

f = log(gaussOverDiffCumGaussian(u, uprime, 1)+1e-300) ...
    + .5*u.*u + .5*log(2*pi);
f=-f;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [lf] = eeg_halfspace_dipole(dippos, elc, vol)

% EEG_HALFSPACE_DIPOLE calculate the leadfield on electrode positions elc
% for a dipole at position dippos. The halfspace solution requires a plane dividing a
% conductive zone (cond > 0), from a non-coductive zone (cond = 0).
%
% Use as
%   [lf] = eeg_halfspace_dipole(dippos, elc, vol)
%
% See also EEG_INFINITE_DIPOLE, EEG_INFINITE_MONOPOLE, EEG_HALFSPACE_MONOPOLE

% Copyright (C) 2011, Cristiano Micheli and Robert Oostenveld
%
% 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$

if isfield(vol, 'pos')
  % this is for forward/backward compatibility
  vol.pnt = vol.pos;
  vol = rmfield(vol, 'pos');
end

siz = size(dippos);
if any(siz==1)
  % positions are specified as a single vector
  Ndipoles = prod(siz)/3;
  dippos = dippos(:)'; % ensure that it is a row vector
elseif siz(2)==3
  % positions are specified as a Nx3 matrix -> reformat to a single vector
  Ndipoles = siz(1);
  dippos = dippos';
  dippos = dippos(:)'; % ensure that it is a row vector
else
  ft_error('incorrect specification of dipole locations');
end

Nelc     = size(elc,1);
lf       = zeros(Nelc,3*Ndipoles);

for i=1:Ndipoles
  % this is the position of dipole "i"
  dip1 = dippos((1:3) + 3*(i-1));

  % find the position of a mirror dipole symmetric to the plane
  dip2 = get_mirror_pos(dip1, vol);

  % compute the potential of the original and the mirror dipole
  lf1 = eeg_infinite_dipole(dip1, elc, vol);
  lf2 = eeg_infinite_dipole(dip2, elc, vol);

  % the z-direction of the mirror dipole should be swapped
  lf2(:,3) = -lf2(:,3);

  % take the sum of the two dipoles
  lf = lf1 + lf2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function P2 = get_mirror_pos(P1,vol)
% calculates the position of a point symmetric to pnt with respect to a plane

% define the plane
pnt = vol.pnt;
ori = vol.ori; % already normalized

if abs(dot(P1-pnt,ori))<eps
  ft_warning(sprintf ('point %f %f %f lies in the symmetry plane',P1(1),P1(2),P1(3)))
  P2 = P1;
else
  % define the plane in parametric form
  % define a non colinear vector vc with respect to the plane normal
  vc = [1 0 0];
  if abs(cross(ori, vc, 2))<eps
    vc = [0 1 0];
  end
  % define plane's direction vectors
  v1 = cross(ori, vc, 2);  v1 = v1/norm(v1);
  v2 = cross(pnt, ori, 2); v2 = v2/norm(v2);
  plane = [pnt v1 v2];

  % distance plane-point P1
  d = abs(dot(ori, plane(:,1:3)-P1(:,1:3), 2));

  % symmetric point
  P2 = P1 + 2*d*ori;
end


text:
function X = uplus(X)
%UPLUS (overloaded)

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef math_model_opf_acpi_legacy < mp.math_model_opf_acpi & mp.mm_shared_opf_legacy
%MP.MATH_MODEL_OPF_ACPI_LEGACY  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACPI_LEGACY ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        %% constructor
        function obj = math_model_opf_acpi_legacy()
            obj@mp.math_model_opf_acpi();
            if nargin > 0 && isstruct(mpc)
                obj.mpc = mpc;
            end

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if om.var is empty) in ADD_VAR(), DISPLAY() and
            %%              INIT_INDEXED_NAME(), after object construction,
            %%              but before object use.
        end

        function obj = add_named_set(obj, varargin)
            % call parent method (also checks for valid type for named set)
            add_named_set@mp.math_model_opf_acpi(obj, varargin{:});
            obj.add_named_set_legacy(varargin{:});
        end

        function obj = def_set_types(obj)
            obj.def_set_types_legacy();
        end

        function obj = init_set_types(obj)
            init_set_types@mp.math_model_opf_acpi(obj);
            obj.init_set_types_legacy();
        end

        function obj = build(obj, nm, dm, mpopt)
            obj.mpc = dm.source;
            build@mp.math_model_opf_acpi(obj, nm, dm, mpopt);
            obj.build_legacy(nm, dm, mpopt);
        end

        function obj = add_vars(obj, nm, dm, mpopt)
            add_vars@mp.math_model_opf_acpi(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined variables
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_vars(nm, dm, mpopt);
            end
        end

        function add_system_costs(obj, nm, dm, mpopt)
            add_system_costs@mp.math_model_opf_acpi(obj, nm, dm, mpopt);    %% call parent

            %% legacy user-defined costs
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_costs(nm, dm, 0);
            end
        end

        function obj = add_system_constraints(obj, nm, dm, mpopt)
            %% call parent
            add_system_constraints@mp.math_model_opf_acpi(obj, nm, dm, mpopt);

            %% legacy user-defined constraints
            if isfield(dm.userdata, 'legacy_opf_user_mods')
                obj.add_legacy_user_constraints_ac(nm, dm, mpopt);
            end
        end

        function names = legacy_user_var_names(obj)
            names = {'Va', 'Vm', 'Pg', 'Qg'};
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% UPDATE generates html help files for all classes and updates funct.html

d = '/Users/marcelge/Code/DMLT/html';
cd(d);

% read implemented methods

D = dir('../+dml/*.m');

% write method help

s = which('help2html');

f = strfind(s,'/');
if isempty(f)
  cd(d);
  error('path not found');
end

ds = s(1:(f(end)-1));
cd(ds);

for i=1:length(D)
  nme = D(i).name(1:(end-2));
  s = help2html(['dml.' nme]);
  fname = fullfile(d,[nme '.html']);
  bf = exist(fname);
  fid = fopen(fname,'w');
  fprintf(fid,'%s',s);
  fclose(fid);
  if ~bf
    fprintf('added %s.html\n',nme);
    cd(d);
    system(sprintf('git add %s.html',nme));
    cd(ds);
  end
end

% write funct.m

fid = fopen(fullfile(d,'funct.m'),'w');
fprintf(fid,'%%%% DML Toolbox Functions\n%%\n');
for i=1:length(D)
  nme = D(i).name(1:(end-2));
  x = help(sprintf('dml.%s',nme));
  b = strfind(lower(x),nme)+length(nme)+1; b = b(1);
  e = min(b+50,(strfind(x,'.')-1)); e = e(1);
  x = x(b:e);
  fprintf(fid,['%% * <' nme '.html |dml.' nme '|> - ' x '\n']);
end
fclose(fid);

cd(d);

publish('funct.m','outputDir',d);


text:
function [VertexLabelIds, LabelMap, AtlasName] = in_label_bs(FileName)
% IN_LABEL_SVREG: Import an atlas from an SVReg-labelled surface
%
% USAGE: in_label_bs(FileName, Verbosity=1) : Load labeled vertices from given file
%
% INPUT
%       - FileName : full path of SVReg surface file to read
%
% OUTPUT
%       - LabelIds :  vector of vertex label ids
%       - LabelMap : java.util.HashMap of label names and colors for corresponding vertex index.
%           - LabelMap.get(id): cell containing label name and color of id i

% @=============================================================================
% 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: Andrew Krause, 2013

% ===== Read XML Label Description File =====
fPath = bst_fileparts(FileName);

st=strfind(FileName,'.svreg.');
ed=strfind(FileName,'.dfs');
AtlasName = FileName(st+7:ed-1);

XmlFile = file_find(fPath, ['brainsuite_labeldescription*',AtlasName,'.xml']);
if isempty(XmlFile)
    fprintf(1, 'BST> For Atlas %s could not find XML label description file brainsuite_labeldescription.xml\n', AtlasName);
    fprintf(1, 'BST> Only label Ids will be used for this atlas without label names\n');
    LabelMap = [];
else
    LabelMap = generate_label_map(XmlFile);
end

% ===== Get Vertex Labels from DFS FILE =====
VertexLabelIds = get_vertex_label_ids(FileName);

end

%% ======================================================================================
%  ===== HELPER FUNCTIONS ===============================================================
%  ======================================================================================
%% Generate mapping between label id and label name and color
function labelMap = generate_label_map(xmlFileName)
    import javax.xml.xpath.*

    try
        docNode = xmlread(xmlFileName);
    catch
        error('Could not open XML Label Description: %s', xmlFileName);
    end

    factory = XPathFactory.newInstance;
    xPath = factory.newXPath;

    idsExpression    = xPath.compile('//label/@id');
    tagsExpression   = xPath.compile('//label/@fullname');
    colorsExpression = xPath.compile('//label/@color');
    idsList    = idsExpression.evaluate(docNode, XPathConstants.NODESET);
    tagsList   = tagsExpression.evaluate(docNode, XPathConstants.NODESET);
    colorsList = colorsExpression.evaluate(docNode, XPathConstants.NODESET);

    numLabels = idsList.getLength;
    labelMap = java.util.HashMap;
    for i = 1:numLabels
        idNode = idsList.item(i-1);
        tagNode = tagsList.item(i-1);
        colorNode = colorsList.item(i-1);

        % Convert color from hexadecimal to RGB [0, 1]
        colorHex = char(colorNode.getNodeValue);
        while length(colorHex) < 8, colorHex = [colorHex '0']; end
        color = [hex2dec(colorHex(3:4)) hex2dec(colorHex(5:6)) hex2dec(colorHex(7:8))] ./ 255;

        % Add ids, names, and colors to hashmap
        id = char(idNode.getNodeValue);
        labelName = char(tagNode.getNodeValue);
        if strncmp(labelName, 'L.', 2) || strncmp(labelName, 'R.', 2)
            labelName = labelName(4:end);
        end
        labelMap.put(id, {labelName, color});
    end
end

%% Return the label ids of each vertex
function Labels = get_vertex_label_ids(s)
    fid=fopen(s,'r');
    if (fid<0)
        error('Unable to read surface: %s',s);
    end

    %% ===== Read number of vertices in surface and labels =====
    fseek(fid, 28, -1); % Skip to nVertices
    nVertices = fread(fid,1,'int32');
    fseek(fid, 20, 0); % Skip to labelOffset
    LabelOffset = fread(fid,1,'int32');

    if (LabelOffset > 0)
        fseek(fid, LabelOffset, -1);
        Labels = fread(fid, nVertices,'uint16');
    else
        error('Surface does not have label attributes');
    end
    fclose(fid);
end

author: brainstorm-tools


text:
function [all] = VBA_classification(X,y,k,verbose,options,sparse)
% performs binary classification using VBA
% function [pv,stat,df,all] = VBA_classification(X,y,k,verbose,options)
% In brief, this function fits the following logistic regresion model:
%   p(y=1) = s(X*beta)
% where y is the binary data, s is the standard sigmoid mapping, X is the
% design matrix and beta are unknown weight parameters.
% Statistical testing is performed using cross-validation scheme, such that
% the p-value computes the probability of the test score under the null H0.
% In cross-validation, the test score is the binary accuracy outcome
% sampled across test sets. The null distribution accounts for potential
% imbalance in the data in that the first-order moment of the corresponding
% binomial distribution is specified in terms of the sample mean of the
% data. VBA_classification uses a k-fold cross-validation scheme, which
% involves partitioning the data into k subsamples of equal size, each of
% which is used as the test set in turn.
% Note: VBA_classification also performs a fully Bayesian model inversion
% on the entire dataset (full-data inversion).
% IN:
%   - X: nXp design matrix
%   - y: nX1 binary data matrix
%   - k: the number of folds for the k-fold cross-validation scheme. If k
%   is set to n (default), then VBA_classification uses a leave-one-out
%   scheme. If k is set to 0, then the cross-validation scheme is entirely
%   skipped (only the full-data inversion is performed).
%   - verbose: flag for displaying results (default is 0)
%   - options: structure containing VBA's options (can be used for passing,
%   e.g., priors on classification weight parameters)
%   - sparse: when sparse=1, VBA_classification uses sparsifying priors
%   (default=0).
% OUT:
%   - all: structure array with fields:
%       .stat: structure of summary statistics:
%           .pv: classical p-value on classifier accuracy
%           .success: nx1 vector of successful classifications
%           .pa: cross-validation prediction accuracy
%           .bpa: balanced prediction accuracy
%           .pBayes: Bayesian exceedance prob. on prediction accuracy
%       .P: pxk matrix of estimated classification weights
%       .r: data imbalance (r=0.5 means balanced data)
%       .in: structure storing the inputs to VBA_classification
%       .date: date stamp
%       .dt: computing time
%       .handles: a structure containing the handles of the graphical
%       objects (filled-in by VBA_classication_display.m)
% Note: cross-validation results can be displayed using the following
% command line: VBA_classification_display(all).

% fill in default I/O
tStart = tic;
all = [];
[n,m] = size(y);
[n0,p] = size(X);

try,k;catch,k=n;end
try,verbose;catch,verbose=0;end
try % check whether VBA's options structure has been sepecified
    options;
catch
    options = [];
end
try % use sparsify tansform?
    sparse = ~~sparse(1);
catch
    sparse = 0;
end

% check basic numerical requirements
try
    if VBA_isWeird (y)
        disp('Error: data contains weird values!')
        return
    end
end
if ~ VBA_isBinary (y)
    disp('Error: data should be binary!')
    return
end
if ~isequal(n,n0)
    disp('Error: design matrix has to have as many rows as the data matrix!')
    return
end
if ~isequal(m,1)
    disp('Error: data should be a nx1 vector!')
    return
end
if k>n || k<0
    disp('Warning: number of folds reduced to data length!')
    k = n;
end

% randomly re-order the data (for k-fold partitioning)
if k ~= n
    io = randperm(n);
    y = y(io);
    X = X(io,:);
end

% specify options and priors for VBA inversion
dim.p = n;
dim.n_t = 1;
dim.n_phi = p;
dim.n_theta = 0;
dim.n = 0;
g_fname = @g_classif0;
options.sources = struct('type',1,'out',1);
options.DisplayWin = 0;
options.verbose = 0;
options.inG.X = X';
options.inG.sparse = sparse;
options.n0 = 0; % number of dummy counts

if ~isequal(k,0) % performing cross-validation scheme
    if verbose
        disp(' ')
        fprintf(1,'Performing k-folds cross-validation scheme...')
        fprintf(1,'%6.2f %%',0)
        et0 = clock; % get time
    end
    sizeFolds = floor(n./k);
    acc = zeros(n,1); % test accuracy
    acc0 = zeros(n,1); % training accuracy
    Eg = zeros(n,1); % out-of-sample prediction E[y] (on test data)
    Vg = zeros(n,1); % out-of-sample prediction V[y] (on test data)
    P = zeros(p,k);
    for i=1:k
        if i<k
            itest = (i-1)*sizeFolds+1:i*sizeFolds;
        else % last test fold contains all remaining data
            itest = (i-1)*sizeFolds+1:n;
        end
        options.isYout = zeros(n,1);
        options.isYout(itest) = 1;
        [posterior,out] = VBA_NLStateSpaceModel(y,[],[],g_fname,dim,options);
        Eg(itest) = out.suffStat.gx(itest);
        Vg(itest) = out.suffStat.vy(itest);
        ytest = Eg(itest)>=0.5;
        acc(itest) = [y(itest)==ytest];
        acc0(itest) = out.fit.acc;
        if sparse
            P(:,i) = VBA_sparsifyPrior (posterior.muPhi);
        else
            P(:,i) = posterior.muPhi;
        end
        if verbose
            fprintf(1,repmat('\b',1,8))
            fprintf(1,'%6.2f %%',floor(100*i/k))
        end
    end
    if verbose
        fprintf(1,repmat('\b',1,8))
        fprintf(1,[' OK (took ',num2str(etime(clock,et0)),' seconds).'])
        fprintf(1,'\n')
    end
end

% performing full-data inversion
if verbose
    fprintf(1,'Performing whole-data inversion...')
    et0 = clock; % get time
end
options.isYout = zeros(n,1);
[posterior,out] = VBA_NLStateSpaceModel(y,[],[],g_fname,dim,options);
if verbose
    fprintf(1,[' OK (took ',num2str(etime(clock,et0)),' seconds).'])
    fprintf(1,'\n')
end

% wrap-up
all.in.X = X;
all.in.y = y;
all.in.k = k;
all.in.sparse = sparse;
all.date = clock;
all.dt = toc(tStart);
all.posterior = posterior;
all.out = out;
if ~isequal(k,0)
    r = mean(y);
    if r<0.5
        r = 1-r;
    end
    nok = sum(acc);
    [pdf0,cdf0] = VBA_binomial(nok,n,r);
    all.stat.pv = 1 - cdf0; % classical p-value on classifier accuracy
    all.stat.success = acc; % nx1 vector of successful classifications
    all.stat.pa = nok./n; % cross-validation prediction accuracy
    all.stat.bpa = 0.5*(sum(acc.*y)./sum(y) + sum(acc.*(1-y))./sum(1-y)); % balanced class. acc.
    all.stat.pBayes = VBA_PPM(nok+options.n0,n-nok+options.n0,r,'beta',0); % Bayesian exceedance prob.
    all.Eg = Eg;
    all.Vg = Vg;
    all.P = P; % set of estimated weights (for each k-fold)
    all.r = r; % data imbalance (r=0.5 means balanced data)
    if verbose
        [all] = VBA_classification_display(all);
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function P = putcolor00(A00, sizeP)
%------------------------------------------------------------------------------
%
% Gridfunction A00 is upsampled.
%
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 6, 1999.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[n, m] = size(A00);
if nargin == 2
  nP = sizeP(1);
  mP = sizeP(2);
  if nP < 2*n-1
    error(' putcolor00 - 1st dimension of P too small ')
  end
  if mP < 2*m-1
    error(' putcolor00 - 2nd dimension of P too small ')
  end
elseif nargin == 1
  nP = 2*n-1;
  mP = 2*m-1;
else
  error(' putcolor00 - wrong number of arguments ')
end
P=reshape(linspace(0,0,nP*mP),nP,mP);
P(1:2:nP, 1:2:mP)=A00;
%------------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function tests = test_spx_checks
  tests = functiontests(localfunctions);
end

function setupOnce(testCase)
    import spx.matrix;
end

function test_is_square(testCase)
    import spx.matrix;
    a = randn(3, 4);
    verifyFalse(testCase, matrix.is_square(a));
    verifyTrue(testCase, matrix.is_square(zeros(3, 3)));
end

function test_is_symmetric(testCase)
    import spx.matrix;
    a  = randn(3, 4);
    verifyFalse(testCase, matrix.is_symmetric(a));
    a = [1 2; 3 4];
    verifyFalse(testCase, matrix.is_symmetric(a));
    a = [1 2; 2 4];
    verifyTrue(testCase, matrix.is_symmetric(a));
    a = [1 2i; 2i 4];
    verifyTrue(testCase, matrix.is_symmetric(a));
    a = [1 2i; -2i 4];
    verifyFalse(testCase, matrix.is_symmetric(a));
end

function test_is_hermitian(testCase)
    import spx.matrix;
    a  = randn(3, 4);
    verifyFalse(testCase, matrix.is_hermitian(a));
    a = [1 2; 3 4];
    verifyFalse(testCase, matrix.is_hermitian(a));
    a = [1 2; 2 4];
    verifyTrue(testCase, matrix.is_hermitian(a));
    a = [1 2i; 2i 4];
    verifyFalse(testCase, matrix.is_hermitian(a));
    a = [1 2i; -2i 4];
    verifyTrue(testCase, matrix.is_hermitian(a));
end

function test_is_pd(testCase)
    import spx.matrix;
    a = randn(3, 3);
    verifyFalse(testCase, matrix.is_positive_definite(a));
    a  = [
    2 -1 0;
    -1 2 -1;
    0 -1 2];
    verifyTrue(testCase, matrix.is_positive_definite(a));
end

author: indigits

repo: sparse-plex


text:
function title = p13_title ( )

%*****************************************************************************80
%
%% P13_TITLE returns a title for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Trigonometric Function.';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function events = in_events_video(sFile, ChannelMat, EventFile, format)
% IN_EVENTS_VIDEO: Read video events information from a text file
%
% USAGE:  events = in_events_video(sFile, ChannelMat, EventFile)
%
%   EventFile must be text file in the form hh:mm:ss:ff

% @=============================================================================
% 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, Elizabeth Bock, 2012-2014

%% ===== READ FILE =====
% Read text file
fid = fopen(EventFile);

nEvents = 0;
while(1)
    tline = fgetl(fid);
    if tline<0
        break;
    end
    % Convert string hh:mm:ss:ff -> hhmmssff
    spl = str_split(tline,':');
    nEvents = nEvents + 1;
    EventsMat(nEvents) = str2double([spl{1:4}]);
end
fclose(fid)

%% ===== CONVERT VIDEO TIME TO FILE TIME =====
% read the video channel from the recording
iVideo = find(strcmpi({ChannelMat.Channel(:).Type}, 'Video'));
if isempty(iVideo)
    error('No video time channel in this file');
end
[F, TimeVector] = in_fread(sFile, ChannelMat, [], round(sFile.prop.times .* sFile.prop.sfreq), iVideo);

% Read channel data for each event
iEvents = bst_closest(EventsMat, F);
eveTimes = TimeVector(iEvents);

%% ===== CONVERT TO BRAINSTORM STRUCTURE =====
% Initialize list of events
events = db_template('event');
% Ask for a label
res = java_dialog('input', 'Please enter a label for this event:', 'Event Label');
if isempty(res)
    events.label = '1';
else
    events.label = res;
end
events.times      = eveTimes;
events.epochs     = ones(1, length(eveTimes)); % Epoch: set as 1 for all the occurrences
events.color      = [];
events.reactTimes = [];
events.select     = 1;
events.channels   = [];
events.notes      = [];


text:
function test_ft_math

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_math

%% create some test data
raw1 = [];
raw1.time = {[1:10], [1:10], [1:10]};
raw1.trial = {ones(2,10), ones(2,10), ones(2,10)};
raw1.label = {'chan01';'chan02'};
raw1.trialinfo = rand(3,4);

timelock1.label  = {'chan1'; 'chan2'};
timelock1.time   = 1:5;
timelock1.dimord = 'chan_time';
timelock1.avg    = ones(2,5);
timelock1.cfg    = struct([]);

timelock2 = timelock1;
timelock2.avg  = ones(2,5)*2;

timelock3 = timelock1;
timelock3.var = ones(2,5);

source1 = [];
source1.pos = randn(10,3);
source1.pow = randn(10,1);
source1.powdimord = 'pos';
source1.mom = cell(10,1);
for i=1:10
  source1.mom{i} = ones(3, 20)*1;
end
source1.momdimord = '{pos}_ori_time';
source1.time = 1:20;

source2 = source1;
for i=1:10
  source2.mom{i} = ones(3, 20)*2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with one raw input

cfg=[];
cfg.showcallinfo = 'no';
cfg.parameter = 'trial';

cfg.operation = 'log10';
tmp = ft_math(cfg, raw1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');

cfg.scalar = pi;

cfg.operation = 'add';
tmp = ft_math(cfg, raw1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');

cfg.operation = 'subtract';
tmp = ft_math(cfg, raw1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, raw1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');

cfg.operation = 'divide';
tmp = ft_math(cfg, raw1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with one timelock input

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'avg';

cfg.operation = 'log10';
tmp = ft_math(cfg, timelock1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.scalar = pi;

cfg.operation = 'add';
tmp = ft_math(cfg, timelock1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'subtract';
tmp = ft_math(cfg, timelock1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, timelock1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'divide';
tmp = ft_math(cfg, timelock1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with one timelock input, multiple parameters

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = {'avg', 'var'};

cfg.operation = 'log10';
tmp = ft_math(cfg, timelock3);
assert(isfield(tmp, 'avg') && isfield(tmp, 'var'), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');
assert(isequal(tmp.var, log10(timelock3.var)));
assert(isequal(tmp.avg, log10(timelock3.avg)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with one source input

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'mom'; % note that this is {pos}_ori_time

cfg.operation = 'log10';
tmp = ft_math(cfg, source1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.scalar = pi;

cfg.operation = 'add';
tmp = ft_math(cfg, source1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'subtract';
tmp = ft_math(cfg, source1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, source1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'divide';
tmp = ft_math(cfg, source1);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with two timelock inputs

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'avg';

cfg.operation = 'add';
tmp = ft_math(cfg, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'subtract';
tmp = ft_math(cfg, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'divide';
tmp = ft_math(cfg, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with two source inputs

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'mom';

cfg.operation = 'add';
tmp = ft_math(cfg, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'subtract';
tmp = ft_math(cfg, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'divide';
tmp = ft_math(cfg, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with more than two timelock inputs

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'avg';

cfg.operation = 'add';
tmp = ft_math(cfg, timelock1, timelock2, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, timelock1, timelock2, timelock1, timelock2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, 'dimord'), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with more than two source inputs

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'mom';

cfg.operation = 'add';
tmp = ft_math(cfg, source1, source2, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

cfg.operation = 'multiply';
tmp = ft_math(cfg, source1, source2, source1, source2);
assert(isfield(tmp, cfg.parameter), 'the output parameter is missing');
assert(isfield(tmp, [cfg.parameter 'dimord']), 'the output dimord is missing');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check the numerical output of the operation

cfg = [];
cfg.parameter = 'trial';

cfg.operation = 'log10';
tmp = ft_math(cfg, raw1);
assert(tmp.trial{1}(1)==0);

cfg.operation = 'multiply';
cfg.scalar = -1;
tmp = ft_math(cfg, raw1);
assert(tmp.trial{1}(1)==-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% idem for a timelock input

cfg = [];
cfg.parameter = 'avg';
cfg.operation = 'add';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==3);

cfg.operation = 'subtract';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==-1);

cfg.operation = 'divide';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==1/2);

cfg.operation = 'multiply';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==2);

cfg.operation = 'log10';
tmp = ft_math(cfg, timelock1);
assert(tmp.avg(1)==0);

cfg.operation = 'log(x1)';
tmp = ft_math(cfg, timelock1);
assert(tmp.avg(1)==0);

cfg.operation = '(x1-x2)/(x1+x2)';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==-1/3);

cfg.scalar = 2;
cfg.operation = '(x1+x2)^s';
tmp = ft_math(cfg, timelock1, timelock2);
assert(tmp.avg(1)==9);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% idem for a source structure with a cell-array

cfg = [];
cfg.parameter = 'mom';
cfg.operation = 'add';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==3);

cfg.operation = 'subtract';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==-1);

cfg.operation = 'divide';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==1/2);

cfg.operation = 'multiply';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==2);

cfg.operation = 'log10';
tmp = ft_math(cfg, source1);
assert(tmp.mom{1}(1)==0);

cfg.operation = 'log(x1)';
tmp = ft_math(cfg, source1);
assert(tmp.mom{1}(1)==0);

cfg.operation = '(x1-x2)/(x1+x2)';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==-1/3);

cfg.scalar = 2;
cfg.operation = '(x1+x2)^s';
tmp = ft_math(cfg, source1, source2);
assert(tmp.mom{1}(1)==9);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% do operation with a matrix

cfg = [];
cfg.showcallinfo = 'no';
cfg.parameter    = 'avg';

cfg.matrix = zeros(size(timelock1.avg));

cfg.operation = 'add';
tmp = ft_math(cfg, timelock1);
assert(isequal(timelock1.avg, tmp.avg));
cfg.operation = 'subtract';
tmp = ft_math(cfg, timelock1);
assert(isequal(timelock1.avg, tmp.avg));

cfg.matrix = ones(size(timelock1.avg));

cfg.operation = 'multiply';
tmp = ft_math(cfg, timelock1);
assert(isequal(timelock1.avg, tmp.avg));
cfg.operation = 'divide';
tmp = ft_math(cfg, timelock1);
assert(isequal(timelock1.avg, tmp.avg));

try
  failed = false;
  cfg = [];
  cfg.showcallinfo = 'no';
  cfg.parameter    = 'avg';
  cfg.operation    = 'log10';
  cfg.matrix       = randn(3);
  tmp = ft_math(cfg, timelock1);
catch
  failed = true;
end
assert(failed, 'this should have failed');


text:
function [x cost info] = trustregions(problem, x, options)
% Riemannian trust-regions solver for optimization on manifolds.
%
% function [x cost info] = trustregions(problem)
% function [x cost info] = trustregions(problem, x0)
% function [x cost info] = trustregions(problem, x0, options)
% function [x cost info] = trustregions(problem, [], options)
%
% This is the Riemannian Trust-Region solver (with tCG inner solve), named
% RTR. This solver will attempt to minimize the cost function described in
% the problem structure. It requires the availability of the cost function
% and of its gradient. It will issue calls for the Hessian. If no Hessian
% nor approximate Hessian is provided, a standard approximation of the
% Hessian based on the gradient will be computed. If a preconditioner for
% the Hessian is provided, it will be used.
%
% For a description of the algorithm and theorems offering convergence
% guarantees, see the references below. Documentation for this solver is
% available online at:
%
% http://www.manopt.org/solver_documentation_trustregions.html
%
%
% The initial iterate is x0 if it is provided. Otherwise, a random point on
% the manifold is picked. To specify options whilst not specifying an
% initial iterate, give x0 as [] (the empty matrix).
%
% The two outputs 'x' and 'cost' are the last reached point on the manifold
% and its cost. Notice that x is not necessarily the best reached point,
% because this solver is not forced to be a descent method. In particular,
% very close to convergence, it is sometimes preferable to accept very
% slight increases in the cost value (on the order of the machine epsilon)
% in the process of reaching fine convergence. In practice, this is not a
% limiting factor, as normally one does not need fine enough convergence
% that this becomes an issue.
%
% The output 'info' is a struct-array which contains information about the
% iterations:
%   iter (integer)
%       The (outer) iteration number, or number of steps considered
%       (whether accepted or rejected). The initial guess is 0.
%	cost (double)
%       The corresponding cost value.
%	gradnorm (double)
%       The (Riemannian) norm of the gradient.
%	numinner (integer)
%       The number of inner iterations executed to compute this iterate.
%       Inner iterations are truncated-CG steps. Each one requires a
%       Hessian (or approximate Hessian) evaluation.
%	time (double)
%       The total elapsed time in seconds to reach the corresponding cost.
%	rho (double)
%       The performance ratio for the iterate.
%	rhonum, rhoden (double)
%       Regularized numerator and denominator of the performance ratio:
%       rho = rhonum/rhoden. See options.rho_regularization.
%	accepted (boolean)
%       Whether the proposed iterate was accepted or not.
%	stepsize (double)
%       The (Riemannian) norm of the vector returned by the inner solver
%       tCG and which is retracted to obtain the proposed next iterate. If
%       accepted = true for the corresponding iterate, this is the size of
%       the step from the previous to the new iterate. If accepted is
%       false, the step was not executed and this is the size of the
%       rejected step.
%	Delta (double)
%       The trust-region radius at the outer iteration.
%	cauchy (boolean)
%       Whether the Cauchy point was used or not (if useRand is true).
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached at each (outer) iteration.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below
%       this. For well-scaled problems, a rule of thumb is that you can
%       expect to reduce the gradient norm by 8 orders of magnitude
%       (sqrt(eps)) compared to the gradient norm at a "typical" point (a
%       rough initial iterate for example). Further decrease is sometimes
%       possible, but inexact floating point arithmetic will eventually
%       limit the final accuracy. If tolgradnorm is set too low, the
%       algorithm may end up iterating forever (or at least until another
%       stopping criterion triggers).
%   maxiter (1000)
%       The algorithm terminates if maxiter (outer) iterations were executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%	miniter (3)
%       Minimum number of outer iterations (used only if useRand is true).
%	mininner (1)
%       Minimum number of inner iterations (for tCG).
%	maxinner (problem.M.dim() : the manifold's dimension)
%       Maximum number of inner iterations (for tCG).
%	Delta_bar (problem.M.typicaldist() or sqrt(problem.M.dim()))
%       Maximum trust-region radius. If you specify this parameter but not
%       Delta0, then Delta0 will be set to 1/8 times this parameter.
%   Delta0 (Delta_bar/8)
%       Initial trust-region radius. If you observe a long plateau at the
%       beginning of the convergence plot (gradient norm VS iteration), it
%       may pay off to try to tune this parameter to shorten the plateau.
%       You should not set this parameter without setting Delta_bar.
%	useRand (false)
%       Set to true if the trust-region solve is to be initiated with a
%       random tangent vector. If set to true, no preconditioner will be
%       used. This option is set to true in some scenarios to escape saddle
%       points, but is otherwise seldom activated.
%	kappa (0.1)
%       Inner kappa convergence tolerance.
%	theta (1.0)
%       Inner theta convergence tolerance.
%	rho_prime (0.1)
%       Accept/reject ratio : if rho is at least rho_prime, the outer
%       iteration is accepted. Otherwise, it is rejected. In case it is
%       rejected, the trust-region radius will have been decreased.
%       To ensure this, rho_prime must be strictly smaller than 1/4.
%   rho_regularization (1e3)
%       Close to convergence, evaluating the performance ratio rho is
%       numerically challenging. Meanwhile, close to convergence, the
%       quadratic model should be a good fit and the steps should be
%       accepted. Regularization lets rho go to 1 as the model decrease and
%       the actual decrease go to zero. Set this option to zero to disable
%       regularization (not recommended). See in-code for the specifics.
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information. statsfun is
%       called with the point x that was reached last, after the
%       accept/reject decision. See comment below.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (2)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent. 3 and above includes a
%       display of the options structure at the beginning of the execution.
%   debug (false)
%       Set to true to allow the algorithm to perform additional
%       computations for debugging purposes. If a debugging test fails, you
%       will be informed of it, usually via the command window. Be aware
%       that these additional computations appear in the algorithm timings
%       too.
%   storedepth (20)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb. If the
%       caching features of Manopt are not used, this is irrelevant. If
%       memory usage is an issue, you may try to lower this number.
%       Profiling may then help to investigate if a performance hit was
%       incured as a result.
%
% Notice that statsfun is called with the point x that was reached last,
% after the accept/reject decision. Hence: if the step was accepted, we get
% that new x, with a store which only saw the call for the cost and for the
% gradient. If the step was rejected, we get the same x as previously, with
% the store structure containing everything that was computed at that point
% (possibly including previous rejects at that same point). Hence, statsfun
% should not be used in conjunction with the store to count operations for
% example. Instead, you could use a global variable and increment that
% variable directly from the cost related functions. It is however possible
% to use statsfun with the store to compute, for example, alternate merit
% functions on the point x.
%
% See also: steepestdescent conjugategradient manopt/examples

% This file is part of Manopt: www.manopt.org.
% This code is an adaptation to Manopt of the original GenRTR code:
% RTR - Riemannian Trust-Region
% (c) 2004-2007, P.-A. Absil, C. G. Baker, K. A. Gallivan
% Florida State University
% School of Computational Science
% (http://www.math.fsu.edu/~cbaker/GenRTR/?page=download)
% See accompanying license file.
% The adaptation was executed by Nicolas Boumal.
%
% Change log:
%
%   NB April 3, 2013:
%       tCG now returns the Hessian along the returned direction eta, so
%       that we do not compute that Hessian redundantly: some savings at
%       each iteration. Similarly, if the useRand flag is on, we spare an
%       extra Hessian computation at each outer iteration too, owing to
%       some modifications in the Cauchy point section of the code specific
%       to useRand = true.
%
%   NB Aug. 22, 2013:
%       This function is now Octave compatible. The transition called for
%       two changes which would otherwise not be advisable. (1) tic/toc is
%       now used as is, as opposed to the safer way:
%       t = tic(); elapsed = toc(t);
%       And (2), the (formerly inner) function savestats was moved outside
%       the main function to not be nested anymore. This is arguably less
%       elegant, but Octave does not (and likely will not) support nested
%       functions.
%
%   NB Dec. 2, 2013:
%       The in-code documentation was largely revised and expanded.
%
%   NB Dec. 2, 2013:
%       The former heuristic which triggered when rhonum was very small and
%       forced rho = 1 has been replaced by a smoother heuristic which
%       consists in regularizing rhonum and rhoden before computing their
%       ratio. It is tunable via options.rho_regularization. Furthermore,
%       the solver now detects if tCG did not obtain a model decrease
%       (which is theoretically impossible but may happen because of
%       numerical errors and/or because of a nonlinear/nonsymmetric Hessian
%       operator, which is the case for finite difference approximations).
%       When such an anomaly is detected, the step is rejected and the
%       trust region radius is decreased.
%
%   NB Dec. 3, 2013:
%       The stepsize is now registered at each iteration, at a small
%       additional cost. The defaults for Delta_bar and Delta0 are better
%       defined. Setting Delta_bar in the options will automatically set
%       Delta0 accordingly. In Manopt 1.0.4, the defaults for these options
%       were not treated appropriately because of an incorrect use of the
%       isfield() built-in function.

% Verify that the problem description is sufficient for the solver.
if ~canGetCost(problem)
    warning('manopt:getCost', ...
            'No cost provided. The algorithm will likely abort.');
end
if ~canGetGradient(problem)
    warning('manopt:getGradient', ...
            'No gradient provided. The algorithm will likely abort.');
end
if ~canGetHessian(problem)
    warning('manopt:getHessian:approx', ...
            'No Hessian provided. Using an approximation instead.');
end

% Define some strings for display
tcg_stop_reason = {'negative curvature',...
                   'exceeded trust region',...
                   'reached target residual-kappa',...
                   'reached target residual-theta',...
                   'dimension exceeded',...
                   'model increased'};

% Set local defaults here
localdefaults.verbosity = 2;
localdefaults.maxtime = inf;
localdefaults.miniter = 3;
localdefaults.maxiter = 1000;
localdefaults.mininner = 1;
localdefaults.maxinner = problem.M.dim();
localdefaults.tolgradnorm = 1e-6;
localdefaults.kappa = 0.1;
localdefaults.theta = 1.0;
localdefaults.rho_prime = 0.1;
localdefaults.useRand = false;
localdefaults.rho_regularization = 1e3;

% Merge global and local defaults, then merge w/ user options, if any.
localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

% Set default Delta_bar and Delta0 separately to deal with additional
% logic: if Delta_bar is provided but not Delta0, let Delta0 automatically
% be some fraction of the provided Delta_bar.
if ~isfield(options, 'Delta_bar')
    if isfield(problem.M, 'typicaldist')
        options.Delta_bar = problem.M.typicaldist();
    else
        options.Delta_bar = sqrt(problem.M.dim());
    end
end
if ~isfield(options,'Delta0')
    options.Delta0 = options.Delta_bar / 8;
end

% Check some option values
assert(options.rho_prime < 1/4, ...
        'options.rho_prime must be strictly smaller than 1/4.');
assert(options.Delta_bar > 0, ...
        'options.Delta_bar must be positive.');
assert(options.Delta0 > 0 && options.Delta0 < options.Delta_bar, ...
        'options.Delta0 must be positive and smaller than Delta_bar.');

% It is sometimes useful to check what the actual option values are.
if options.verbosity >= 3
    disp(options);
end

% Create a store database
storedb = struct();

tic();

% If no initial point x is given by the user, generate one at random.
if ~exist('x', 'var') || isempty(x)
    x = problem.M.rand();
end

%% Initializations

% k counts the outer (TR) iterations. The semantic is that k counts the
% number of iterations fully executed so far.
k = 0;

% initialize solution and companion measures: f(x), fgrad(x)
[fx fgradx storedb] = getCostGrad(problem, x, storedb);
norm_grad = problem.M.norm(x, fgradx);

% initialize trust-region radius
Delta = options.Delta0;

% Save stats in a struct array info, and preallocate
% (see http://people.csail.mit.edu/jskelly/blog/?x=entry:entry091030-033941)
if ~exist('used_cauchy', 'var')
    used_cauchy = [];
end
stats = savestats(problem, x, storedb, options, k, fx, norm_grad, Delta);
info(1) = stats;
info(min(10000, options.maxiter+1)).iter = [];

% ** Display:
if options.verbosity == 2
   fprintf(['%3s %3s      %5s                %5s     ',...
            'f: %e   |grad|: %e\n'],...
           '   ','   ','     ','     ', fx, norm_grad);
elseif options.verbosity > 2
   fprintf('************************************************************************\n');
   fprintf('%3s %3s    k: %5s     num_inner: %5s     %s\n',...
           '','','______','______','');
   fprintf('       f(x) : %e       |grad| : %e\n', fx, norm_grad);
   fprintf('      Delta : %f\n', Delta);
end

% **********************
% ** Start of TR loop **
% **********************
while true

	% Start clock for this outer iteration
    tic();

    % Run standard stopping criterion checks
    [stop reason] = stoppingcriterion(problem, x, options, info, k+1);

    % If the stopping criterion that triggered is the tolerance on the
    % gradient norm but we are using randomization, make sure we make at
    % least miniter iterations to give randomization a chance at escaping
    % saddle points.
    if stop == 2 && options.useRand && k < options.miniter
        stop = 0;
    end

    if stop
        if options.verbosity >= 1
            fprintf([reason '\n']);
        end
        break;
    end

    if options.verbosity > 2 || options.debug > 0
        fprintf('************************************************************************\n');
    end

    % *************************
    % ** Begin TR Subproblem **
    % *************************

    % Determine eta0
    if ~options.useRand
        % Pick the zero vector
        eta = problem.M.zerovec(x);
    else
        % Random vector in T_x M (this has to be very small)
        eta = problem.M.lincomb(x, 1e-6, problem.M.randvec(x));
        % Must be inside trust-region
        while problem.M.norm(x, eta) > Delta
            eta = problem.M.lincomb(x, sqrt(sqrt(eps)), eta);
        end
    end

    % solve TR subproblem
    [eta Heta numit stop_inner storedb] = ...
                     tCG(problem, x, fgradx, eta, Delta, options, storedb);
    srstr = tcg_stop_reason{stop_inner};

    % This is only computed for logging purposes, because it may be useful
    % for some user-defined stopping criteria. If this is not cheap for
    % specific application (compared to evaluating the cost), we should
    % reconsider this.
    norm_eta = problem.M.norm(x, eta);

    if options.debug > 0
        testangle = problem.M.inner(x, eta, fgradx) / (norm_eta*norm_grad);
    end

    % If using randomized approach, compare result with the Cauchy point.
    % Convergence proofs assume that we achieve at least the reduction of
    % the Cauchy point. After this if-block, either all eta-related
    % quantities have been changed consistently, or none of them have
    % changed.
    if options.useRand
        used_cauchy = false;
        % Check the curvature,
        [Hg storedb] = getHessian(problem, x, fgradx, storedb);
        g_Hg = problem.M.inner(x, fgradx, Hg);
        if g_Hg <= 0
            tau_c = 1;
        else
            tau_c = min( norm_grad^3/(Delta*g_Hg) , 1);
        end
        % and generate the Cauchy point.
        eta_c  = problem.M.lincomb(x, -tau_c * Delta / norm_grad, fgradx);
        Heta_c = problem.M.lincomb(x, -tau_c * Delta / norm_grad, Hg);

        % Now that we have computed the Cauchy point in addition to the
        % returned eta, we might as well keep the best of them.
        mdle  = fx + problem.M.inner(x, fgradx, eta) ...
                   + .5*problem.M.inner(x, Heta,   eta);
        mdlec = fx + problem.M.inner(x, fgradx, eta_c) ...
                   + .5*problem.M.inner(x, Heta_c, eta_c);
        if mdle > mdlec
            eta = eta_c;
            Heta = Heta_c; % added April 11, 2012
            used_cauchy = true;
        end
    end

	% Compute the retraction of the proposal
	x_prop  = problem.M.retr(x, eta);

	% Compute the function value of the proposal
	[fx_prop storedb] = getCost(problem, x_prop, storedb);

	% Will we accept the proposed solution or not?
    % Check the performance of the quadratic model against the actual cost.
    rhonum = fx - fx_prop;
    rhoden = -problem.M.inner(x, fgradx, eta) ...
             -.5*problem.M.inner(x, eta, Heta);

    % Heuristic -- added Dec. 2, 2013 (NB) to replace the former heuristic.
    % This heuristic is documented in the book by Conn Gould and Toint on
    % trust-region methods, section 17.4.2.
    % rhonum measures the difference between two numbers. Close to
    % convergence, these two numbers are very close to each other, so
    % that computing their difference is numerically challenging: there may
    % be a significant loss in accuracy. Since the acceptance or rejection
    % of the step is conditioned on the ratio between rhonum and rhoden,
    % large errors in rhonum result in a large error in rho, hence in
    % erratic acceptance / rejection. Meanwhile, close to convergence,
    % steps are usually trustworthy and we should transition to a Newton-
    % like method, with rho=1 consistently. The heuristic thus shifts both
    % rhonum and rhoden by a small amount such that far from convergence,
    % the shift is irrelevant and close to convergence, the ratio rho goes
    % to 1, effectively promoting acceptance of the step.
    % The rationale is that close to convergence, both rhonum and rhoden
    % are quadratic in the distance between x and x_prop. Thus, when this
    % distance is on the order of sqrt(eps), the value of rhonum and rhoden
    % is on the order of eps, which is indistinguishable from the numerical
    % error, resulting in badly estimated rho's.
    % For abs(fx) < 1, this heuristic is invariant under offsets of f but
    % not under scaling of f. For abs(fx) > 1, the opposite holds. This
    % should not alarm us, as this heuristic only triggers at the very last
    % iterations if very fine convergence is demanded.
    rho_reg = max(1, abs(fx)) * eps * options.rho_regularization;
    rhonum = rhonum + rho_reg;
    rhoden = rhoden + rho_reg;

    if options.debug > 0
        fprintf('DBG:     rhonum : %e\n', rhonum);
        fprintf('DBG:     rhoden : %e\n', rhoden);
    end

    % This is always true if a linear, symmetric operator is used for the
    % Hessian (approximation) and if we had infinite numerical precision.
    % In practice, nonlinear approximations of the Hessian such as the
    % built-in finite difference approximation and finite numerical
    % accuracy can cause the model to increase. In such scenarios, we
    % decide to force a rejection of the step and a reduction of the
    % trust-region radius. We test the sign of the regularized rhoden since
    % the regularization is supposed to capture the accuracy to which
    % rhoden is computed: if rhoden were negative before regularization but
    % not after, that should not be (and is not) detected as a failure.
    model_decreased = (rhoden >= 0);

    if ~model_decreased
        srstr = [srstr ', model did not decrease']; %#ok<AGROW>
    end

    rho = rhonum / rhoden;

    if options.debug > 0
        m = @(x, eta) ...
          getCost(problem, x, storedb) + ...
          getDirectionalDerivative(problem, x, eta, storedb) + ...
          .5*problem.M.inner(x, getHessian(problem, x, eta, storedb), eta);
        zerovec = problem.M.zerovec(x);
        actrho = (fx - fx_prop) / (m(x, zerovec) - m(x, eta));
        fprintf('DBG:   new f(x) : %e\n', fx_prop);
        fprintf('DBG: actual rho : %e\n', actrho);
        fprintf('DBG:   used rho : %e\n', rho);
    end

    % Choose the new TR radius based on the model performance
    trstr = '   ';
    % If the actual decrease is smaller than 1/4 of the predicted decrease,
    % then reduce the TR radius.
    if rho < 1/4 || ~model_decreased
        trstr = 'TR-';
        Delta = Delta/4;
    % If the actual decrease is at least 3/4 of the precicted decrease and
    % the tCG (inner solve) hit the TR boundary, increase the TR radius.
    elseif rho > 3/4 && (stop_inner == 1 || stop_inner == 2)
        trstr = 'TR+';
        Delta = min(2*Delta, options.Delta_bar);
    end
    % Otherwise, keep the TR radius constant.

    % Choose to accept or reject the proposed step based on the model
    % performance.
    if model_decreased && rho > options.rho_prime
        accept = true;
        accstr = 'acc';
        x = x_prop;
        fx = fx_prop;
        [fgradx storedb] = getGradient(problem, x, storedb);
        norm_grad = problem.M.norm(x, fgradx);
    else
        accept = false;
        accstr = 'REJ';
    end

    % Make sure we don't use too much memory for the store database
    storedb = purgeStoredb(storedb, options.storedepth);

    % k is the number of iterations we have accomplished.
    k = k + 1;

    % Log statistics for freshly executed iteration.
    % Everything after this in the loop is not accounted for in the timing.
    stats = savestats(problem, x, storedb, options, k, fx, norm_grad, ...
                      Delta, info, rho, rhonum, rhoden, accept, numit, ...
                      norm_eta, used_cauchy);
    info(k+1) = stats; %#ok<AGROW>

    % ** Display:
    if options.verbosity == 2,
        fprintf(['%3s %3s   k: %5d     num_inner: %5d     ', ...
        'f: %e   |grad|: %e   %s\n'], ...
        accstr,trstr,k,numit,fx,norm_grad,srstr);
    elseif options.verbosity > 2,
        if options.useRand && used_cauchy,
            fprintf('USED CAUCHY POINT\n');
        end
		fprintf('%3s %3s    k: %5d     num_inner: %5d     %s\n', ...
				accstr, trstr, k, numit, srstr);
		fprintf('       f(x) : %e     |grad| : %e\n',fx,norm_grad);
		if options.debug > 0
			fprintf('      Delta : %f          |eta| : %e\n',Delta,norm_eta);
		end
		fprintf('        rho : %e\n',rho);
    end
    if options.debug > 0,
        fprintf('DBG: cos ang(eta,gradf): %d\n',testangle);
        if rho == 0
            fprintf('DBG: rho = 0, this will likely hinder further convergence.\n');
        end
    end

end  % of TR loop (counter: k)

% Restrict info struct-array to useful part
info = info(1:k+1);

if (options.verbosity > 2) || (options.debug > 0),
   fprintf('************************************************************************\n');
end
if (options.verbosity > 0) || (options.debug > 0)
    fprintf('Total time is %f [s] (excludes statsfun)\n', info(end).time);
end

% Return the best cost reached
cost = fx;

end

% Routine in charge of collecting the current iteration stats
function stats = savestats(problem, x, storedb, options, k, fx, ...
                           norm_grad, Delta, info, rho, rhonum, ...
                           rhoden, accept, numit, norm_eta, used_cauchy)
    stats.iter = k;
    stats.cost = fx;
    stats.gradnorm = norm_grad;
    stats.Delta = Delta;
    if k == 0
        stats.time = toc();
        stats.rho = inf;
        stats.rhonum = NaN;
        stats.rhoden = NaN;
        stats.accepted = true;
        stats.numinner = NaN;
        stats.stepsize = NaN;
        if options.useRand
            stats.cauchy = false;
        end
    else
        stats.time = info(k).time + toc();
        stats.rho = rho;
        stats.rhonum = rhonum;
        stats.rhoden = rhoden;
        stats.accepted = accept;
        stats.numinner = numit;
        stats.stepsize = norm_eta;
        if options.useRand,
          stats.cauchy = used_cauchy;
        end
    end

    % See comment about statsfun above: the x and store passed to statsfun
    % are that of the most recently accepted point after the iteration
    % fully executed.
    stats = applyStatsfun(problem, x, storedb, options, stats);

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function ans = is_row (v);
% ans = is_row (v);
%returns 1 if v is a row vector
%returns 0 if v is a column vector
%returns -1 if v is a matrix
%
% Ziad Saad July 20 97

if (size (v,1) ~= 1 & size (v,2) ~= 1),
	ans = -1;
end;

if (size (v,1) == 1),
	ans = 1;
	elseif (size (v,2) == 1),
		ans = 0;
end;

return;

author: fieldtrip