text:
function [ALLEEG cfg] = pop_stat_surrogateStats(ALLEEG,typeproc,varargin)
%
% Return surrogate statistics based on a surrogate distribution (bootstrap, jacknife, etc).
% Each surrogate distribution should approximate the distribution of the
% estimator.
%
%
% Input:
%
%   ALLEEG:         EEGLAB dataset to preprocess. Must contain .CAT.PConn
%                   with surrogate distributions
%   typeproc:       Reserved for future use. Use 0
%
% Optional:
%
%   <'Name',value> pairs as defined in stat_surrogate()
%
% Output:
%
%   ALLEEG:         EEG structure(s) with Stats object stored in
%                   ALLEEG.CAT.Stats
%   cfg:            Argument specification structure.
%
%
% See Also: stat_surrogate()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2010-2011, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

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

if nargin<2
    typeproc = 0;
end

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

% check the dataset
res = hlp_checkeegset(ALLEEG,{'pconn'});
if ~isempty(res)
    error(['SIFT:' fcnName],res{1});
end

if isfield(ALLEEG(1).CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(ALLEEG(1).CAT.configs.(fcnName)) varargin];
end

% reset the defaults if we have multiple EEG datasets
% this ensures that config selections stored in each datset (i.e. Hbase)
% don't conflict with the allowable selection when two datasets are present
if length(ALLEEG)>1 && ~isempty(varargin)
    idx = find(ismember_bc(varargin(1:2:end),'statTest'))*2;
    if ~isempty(idx)
        for k=1:length(idx)
            varargin{idx(k)} = {};
        end
    end
end

if strcmpi(typeproc,'nogui')
    % get the config from function
    cfg = arg_tovals(arg_report('rich',fcnHandle,[{'EEG',ALLEEG},varargin]),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],ALLEEG,varargin{:});

    if isempty(PGh)
        % user chose to cancel
        cfg = [];
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if strcmpi(typeproc,'cfg_only')
    return;
end

% execute the low-level function
[Stats ConnMean] = feval(fcnHandle,'EEG',ALLEEG,cfg);

if isempty(Stats)
    % user canceled
    return;
end

% store statistics and (optionally) the mean of the bootstrap estimator
if length(ALLEEG)>1 && any(strcmp(cfg.statTest.arg_selection,{'Hab'}))
    % create a new dataset containing expected difference between conds
%     EEG2         = ALLEG(2);
%     EEG2.data    = -EEG2.data;  % invert data so average is mean ERP difference
%     EEG2.icaact  = -EEG2.icaact;
%     EEG_new           = pop_mergeset(ALLEEG(1),EEG2,1);

    EEG_new = ALLEEG(1);
    EEG_new.CAT = rmfield(EEG_new.CAT,'PConn');
    % compute ERP condition difference
    for fn={'data','icaact','srcpot'}
        if isfield(EEG_new,fn{1}) && ~isempty(EEG_new.(fn{1}))
            EEG_new.(fn{1}) = mean(ALLEEG(Stats.diffOrder(1)).(fn{1}),3) ...
                            - mean(ALLEEG(Stats.diffOrder(2)).(fn{1}),3);
        end
    end
    EEG_new.trials    = 1;
    EEG_new.epoch     = [];
    EEG_new.setname   = cfg.statTest.datasetOrder;
    EEG_new.condition = cfg.statTest.datasetOrder;
    EEG_new = eeg_checkset(EEG_new);

    EEG_new.CAT.Conn  = ConnMean;  % condition difference
    EEG_new.CAT.configs.(fcnName) = cfg;
    EEG_new.CAT.Stats = Stats;

    EEG_new.CAT.configs.vis_TimeFreqGrid = [];
    EEG_new.CAT.configs.vis_causalBrainMovie3D = [];

    % store the new EEG dataset
    [ALLEEG EEG_new] = eeg_store(ALLEEG,EEG_new,length(ALLEEG)+1);
    ALLEEG = EEG_new;
elseif length(ALLEEG)==1

    if ~isempty(cfg)
        % store the configuration structure
        ALLEEG.CAT.configs.(fcnName) = cfg;
    end

    ALLEEG.CAT.Stats = Stats;

    if ~isempty(ConnMean)
        % replace Conn object with mean of bootstrap distribution
        % insert missing fields into new Conn object
%         extrafields = setdiff_bc(fieldnames(ALLEEG.CAT.Conn),hlp_getConnMethodNames(ALLEEG.CAT.Conn));
%         for i=1:length(extrafields)
%             ConnMean(cnd).(extrafields{i}) = ALLEEG.CAT.Conn.(extrafields{i});
%         end
        ALLEEG.CAT.Conn = ConnMean;

        ALLEEG.CAT.configs.vis_TimeFreqGrid = [];
        ALLEEG.CAT.configs.vis_causalBrainMovie3D = [];
    end
end

repo: REST


text:
function nbrhood=cosmo_spherical_neighborhood(ds, varargin)
% computes neighbors for a spherical searchlight
%
% nbrhood=cosmo_spherical_neighborhood(ds, opt)
%
% Inputs
%   ds                  a dataset struct, either:
%                       - in fmri form (from cosmo_fmri_dataset), when
%                         ds.fa has the fields .i, .j and .k
%                       - in meeg source form (from cosmo_meeg_dataset),
%                         when ds.fa has the field .pos. In this case, the
%                         features must have positions that can be
%                         converted to a grid.
%   'radius', r         } either use a radius of r, or select
%   'count', c          } approximately c voxels per searchlight
%                       Notes:
%                       - These two options are mutually exclusive
%                       - When using this option for an fmri dataset, the
%                         radius r is expressed in voxel units; for an meeg
%                         source dataset, the radius r is in whatever units
%                         the source dataset uses for the positions
%   'progress', p       show progress every p features (default: 1000)
%
% Outputs
%   nbrhood             dataset-like struct without .sa or .samples, with:
%     .a                dataset attributes, from dataset.a
%     .fa               feature attributes with the same fields as fs.fa,
%                       and in addition the fields:
%       .nvoxels        1xP number of voxels in each searchlight
%       .radius         1xP radius in voxel units
%       .center_ids     1xP feature center id
%     .neighbors        Px1 cell so that center2neighbors{k}==nbrs contains
%                       the feature ids of the neighbors of feature k
%                       If the dataset has a field ds.fa.inside, then
%                       features that are not inside are not included as
%                       neighbors in the output
%     .origin           Has fields .a and .fa from input dataset
%
%
% Example:
%     ds=cosmo_synthetic_dataset('type','fmri');
%     radius=1; % radius=3 is typical for 'real-world' searchlights
%     nbrhood=cosmo_spherical_neighborhood(ds,'radius',radius,...
%                                             'progress',false);
%     cosmo_disp(nbrhood)
%     %|| .a
%     %||   .fdim
%     %||     .labels
%     %||       { 'i'  'j'  'k' }
%     %||     .values
%     %||       { [ 1         2         3 ]  [ 1         2 ]  [ 1 ] }
%     %||   .vol
%     %||     .mat
%     %||       [ 2         0         0        -3
%     %||         0         2         0        -3
%     %||         0         0         2        -3
%     %||         0         0         0         1 ]
%     %||     .dim
%     %||       [ 3         2         1 ]
%     %||     .xform
%     %||       'scanner_anat'
%     %|| .fa
%     %||   .nvoxels
%     %||     [ 3         4         3         3         4         3 ]
%     %||   .radius
%     %||     [ 1         1         1         1         1         1 ]
%     %||   .center_ids
%     %||     [ 1         2         3         4         5         6 ]
%     %||   .i
%     %||     [ 1         2         3         1         2         3 ]
%     %||   .j
%     %||     [ 1         1         1         2         2         2 ]
%     %||   .k
%     %||     [ 1         1         1         1         1         1 ]
%     %|| .neighbors
%     %||   { [ 1         4         2 ]
%     %||     [ 2         1         5         3 ]
%     %||     [ 3         2         6 ]
%     %||     [ 4         1         5 ]
%     %||     [ 5         4         2         6 ]
%     %||     [ 6         5         3 ]           }
%     %|| .origin
%     %||   .a
%     %||     .fdim
%     %||       .labels
%     %||         { 'i'
%     %||           'j'
%     %||           'k' }
%     %||       .values
%     %||         { [ 1         2         3 ]
%     %||           [ 1         2 ]
%     %||           [ 1 ]                     }
%     %||     .vol
%     %||       .mat
%     %||         [ 2         0         0        -3
%     %||           0         2         0        -3
%     %||           0         0         2        -3
%     %||           0         0         0         1 ]
%     %||       .dim
%     %||         [ 3         2         1 ]
%     %||       .xform
%     %||         'scanner_anat'
%     %||   .fa
%     %||     .i
%     %||       [ 1         2         3         1         2         3 ]
%     %||     .j
%     %||       [ 1         1         1         2         2         2 ]
%     %||     .k
%     %||       [ 1         1         1         1         1         1 ]
%
%
% Notes:
%   - this function can return neighborhoods with either a fixed number of
%     features, or a fixed radius. When used with a searchlight, the
%     former has the advantage that the number of features is less
%     variable (especially near edges of the brain, in an fmri dataset),
%     which can make it easier to compare result in different regions as
%     the number of features can affect
%     pattern discriminablity. The latter has the advantage that the
%     smoothness of the output maps under the null hypothesis can be more
%     uniformly smooth.
%
% See also: cosmo_fmri_dataset, cosmo_meeg_dataset, cosmo_searchlight
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    check_input(varargin{:});

    defaults=struct();
    defaults.progress=1000;
    opt=cosmo_structjoin(defaults,varargin);

    [use_fixed_radius,radius,voxel_count]=get_selection_params(opt);
    cosmo_check_dataset(ds);

    % ensure not too many features are requested
    feature_mask=get_features_mask(ds);
    nfeatures=sum(feature_mask);
    if nfeatures<voxel_count
        error('Cannot select %d features: only %d are present',...
                    voxel_count, nfeatures);
    end

    % get attributes for output dataset, and the positions and dimension of
    % the grid
    [fdim,fa,pos,grid_dim]=get_spherical_attributes(ds,feature_mask);

    % compute voxel offsets relative to origin
    [sphere_offsets, distances]=get_sphere_offsets(radius);

    % get mapping from linear ids to feature ids
    [lin2feature_ids,lin2feature_mask]=get_lin2feature_ids(grid_dim,pos...
                                        ,feature_mask);

    % number of features associated with each linear id
    feature_id_count=sum(lin2feature_mask,2);

    show_progress=opt.progress>0;

    if show_progress
        clock_start=clock();
        prev_progress_msg='';
    end

    % a position may occur at multiple features; only consider unique
    % positions
    pos(:,~feature_mask)=Inf;
    [center_idxs,unq_pos]=cosmo_index_unique(pos');
    keep_unq_pos=~any(isinf(unq_pos),2);
    center_idxs=center_idxs(keep_unq_pos);
    unq_pos=unq_pos(keep_unq_pos,:);
    nunq_centers=numel(center_idxs);

    % allocate space for output
    ncenters=nunq_centers;
    neighbors=cell(ncenters,1);
    nvoxels=zeros(1,ncenters);
    final_radius=zeros(1,ncenters);
    visited=false(1,ncenters);
    center_ids=zeros(1,ncenters);

    % go over all features
    for k=1:nunq_centers
        variable_radius=NaN;
        if voxel_count==0
            feature_ids=zeros(1,0);
        else
            center_pos=unq_pos(k,:)';

            % - in case of a variable radius, keep growing sphere_offsets
            %   until there are enough voxels selected. This new radius is
            %   kept for every subsequent iteration.
            % - in case of a fixed radius this loop is left after the first
            %   iteration.
            while true
                % add offsets to center
                all_around_pos=bsxfun(@plus, center_pos', sphere_offsets);

                % see which ones are outside the volume
                outside_msk=all_around_pos<=0 | ...
                                bsxfun(@minus,grid_dim,all_around_pos)<0;

                % collapse over 3 dimensions
                feature_outside_msk=any(outside_msk,2);

                % get rid of those outside the volume
                around_pos=all_around_pos(~feature_outside_msk,:);

                % convert to linear indices
                around_lin=fast_sub2ind(grid_dim,around_pos(:,1), ...
                                            around_pos(:,2), ...
                                            around_pos(:,3));

                % convert linear to feature ids
                % (transpose is necessary so that when applying the
                %  mask, the indices remain sorted by distance)
                around_ids_mat=lin2feature_ids(around_lin,:)';
                around_ids_mask=lin2feature_mask(around_lin,:)';
                feature_ids=around_ids_mat(around_ids_mask);

                if use_fixed_radius
                    break; % we're done selecting voxels
                elseif numel(feature_ids)<voxel_count
                    % the current radius is too small.
                    % increase the radius by half a voxel and recompute new
                    % offsets, then try again in the next iteration.
                    radius=radius+.5;
                    [sphere_offsets,distances]=get_sphere_offsets(radius);
                    continue
                end

                % if using variable radius, compute distance of each linear
                % index
                center_distances=distances(~feature_outside_msk);

                % get distance for each feature
                feature_distances=get_distances(center_distances,...
                                            feature_id_count(around_lin));

                % coming here, the radius is variable and enough features
                % were selected. Now decide which voxels to keep,
                % and also compute the metric radius, then leave the while
                % loop.

                nselect=boundary_at_approx(feature_ids,...
                                            feature_distances,voxel_count);
                feature_ids=feature_ids(1:nselect);

                variable_radius=feature_distances(nselect);
                break; % we're done
            end
        end

        % store results
        id=center_idxs{k}(1);

        neighbors{k}=feature_ids(:)';
        nvoxels(k)=numel(feature_ids);
        if use_fixed_radius
            final_radius(k)=radius;
        else
            final_radius(k)=variable_radius;
        end

        visited(k)=true;
        assert(center_ids(k)==0);
        center_ids(k)=id;

        if show_progress && (k==1 || k==nunq_centers || ...
                                        mod(k,opt.progress)==0)
            mean_size=mean(nvoxels(visited));
            msg=sprintf('mean size %.1f', mean_size);
            prev_progress_msg=cosmo_show_progress(clock_start, ...
                                   k/nunq_centers, msg, prev_progress_msg);
        end
    end

    not_visited_ids=find(~visited);
    assert(all(cellfun(@numel,neighbors(not_visited_ids))==0));
    neighbors(not_visited_ids)=repmat({zeros(1,0)},...
                                        1,numel(not_visited_ids));

    % set the dataset and feature attributes
    nbrhood=struct();
    nbrhood.a=ds.a;
    nbrhood.a.fdim=fdim;

    % remove sample dimension if present
    if isfield(nbrhood.a,'sdim')
        nbrhood.a=rmfield(nbrhood.a,'sdim');
    end

    fa_full=cosmo_slice(fa,center_ids,2,'struct');
    nbrhood.fa=cosmo_structjoin('nvoxels',nvoxels,...
                                'radius',final_radius,...
                                'center_ids',center_ids(:)',...
                                fa_full);

    nbrhood.neighbors=neighbors;

    nbrhood=align_nbrhood_to_ds_if_possible(ds,nbrhood);
    origin=struct();
    origin.a=ds.a;
    origin.fa=ds.fa;
    nbrhood.origin=origin;

    cosmo_check_neighborhood(nbrhood,ds);

function nbrhood=align_nbrhood_to_ds_if_possible(ds,nbrhood)
   labels=get_dim_label(ds);

   ds_fa=get_spherical_fa_cell(ds.fa,labels);
   nbrhood_fa=get_spherical_fa_cell(nbrhood.fa,labels);

   [unq_ds,idx_ds]=cosmo_index_unique(ds_fa);
   [unq_nbrhood,idx_nbrhood]=cosmo_index_unique(nbrhood_fa);

   if all(cellfun(@numel,unq_ds)==1) && ...
           isequal(sort(cell2mat(unq_ds)),sort(cell2mat(unq_nbrhood)))
       mp=cosmo_align(nbrhood_fa,ds_fa);

       nbrhood.neighbors=nbrhood.neighbors(mp);
       nbrhood.fa=cosmo_slice(nbrhood.fa,mp,2,'struct');
   end

function feature_distances=get_distances(center_distances,feature_id_count)
    % get distances based on selected features
    n=numel(center_distances);
    assert(n==numel(feature_id_count));

    m=max(feature_id_count);
    if m<=1
        % optimization
        feature_distances=center_distances(feature_id_count==1);
        return
    end

    ds=NaN(m,n);

    for k=1:m
        msk=feature_id_count>=k;
        ds(k,msk)=center_distances(msk);
    end

    feature_distances=ds(~isnan(ds));

function [lin2feature_ids,lin2feature_mask]=get_lin2feature_ids(...
                                            grid_dim,all_pos,center_mask)
    % returns a function that maps linear ids to feature ids
    % the function takes as input linear ids and the distance for each
    % linear id, and returns the feature ids and their corresponding
    % distances

    orig_nvoxels=prod(grid_dim);

    ijk=all_pos(:,center_mask);

    lin_ids=fast_sub2ind(grid_dim, ijk(1,:), ijk(2,:), ijk(3,:));
    [idxs,unq_lin_ids]=cosmo_index_unique(lin_ids');

    mask2full=find(center_mask);
    % lin2feature_ids{k}={i1,...,iN} means that the linear voxel index k
    % corresponds to features i1,...iN
    lin2feature_ids_cell=cell(orig_nvoxels,1);
    for k=1:numel(unq_lin_ids)
        lin_id=unq_lin_ids(k);
        idx=idxs{k}(:)';
        lin2feature_ids_cell{lin_id}=mask2full(idx);
    end

    n_max=max(cellfun(@numel,lin2feature_ids_cell));

    lin2feature_ids=zeros(orig_nvoxels, n_max);
    lin2feature_mask=false(orig_nvoxels,n_max);

    for k=1:numel(unq_lin_ids)
        lin_id=unq_lin_ids(k);
        indices=lin2feature_ids_cell{lin_id};

        cols=1:numel(indices);
        lin2feature_ids(lin_id,cols)=indices;
        lin2feature_mask(lin_id,cols)=true;
    end

function feature_mask=get_features_mask(ds)
    % use .fa.inside if it is present, otherwise an array with only true
    % values
    nfeatures=size(ds.samples,2);

    if cosmo_isfield(ds,'fa.inside')
        inside=ds.fa.inside;

        if size(inside,1)~=1
            error('field .fa.inside must be a row vector');
        end

        if ~islogical(inside)
            error('field .fa.inside must be logical');
        end

        feature_mask=inside;
    else
        feature_mask=true(1,nfeatures);
    end

function lin=fast_sub2ind(sz, i, j, k)
    lin=sz(1)*(sz(2)*(k-1)+(j-1))+i;

function pos=boundary_at_approx(ids, distances, voxel_count)
    % pseudo-random selection of approximatly voxel_count elements
    if voxel_count<=0
        pos=0;
        return
    end

    assert(issorted(distances));

    max_distance=distances(voxel_count);
    first=find(distances<max_distance,1,'last')+1;
    last=find(distances>max_distance,1,'first')-1;

    if isempty(first)
        first=1;
    end

    if isempty(last)
        last=numel(distances);
    end

    delta_first=voxel_count-first;
    delta_last=last-voxel_count;

    if delta_first==delta_last
        % select pseudo-randomly
        if delta_first==0 || mod(sum(ids)+numel(distances),2)==0
            pos=first;
        else
            pos=last;
        end
    elseif delta_first<delta_last
        pos=first;
    else
        pos=last;
    end

    assert(first==1 || distances(first-1)<distances(first));
    assert(last==numel(distances) || distances(last+1)>distances(first));

function [fdim,fa,ijk,orig_dim]=get_spherical_attributes(ds, center_mask)
    % returns fdim, fa, and ijk positions for dataset
    labels=get_dim_label(ds);

    fdim=get_spherical_fdim(ds,labels);
    fa=get_spherical_fa(ds.fa,labels);

    if cosmo_isfield(ds,'fa.inside')
        fa.inside=center_mask;
    end

    small_ds=cosmo_slice(ds,[],1);
    small_ds_vol=cosmo_vol_grid_convert(small_ds, 'tovol');

    ijk=[small_ds_vol.fa.i;small_ds_vol.fa.j;small_ds_vol.fa.k];

    ijk_labels={'i','j','k'};
    [unused,index]=has_fdim_label(small_ds_vol,ijk_labels);

    orig_dim=cellfun(@numel,small_ds_vol.a.fdim.values(index));
    orig_dim=orig_dim(:)';

function [tf,index]=has_fdim_label(ds, label)
    [two,index]=cosmo_dim_find(ds,label,false);
    tf=~isempty(two) && two==2;

function [labels,index]=get_dim_label(ds)
    % get either pos or i, j, and k labels
    possible_labels={{'pos'},{'i';'j';'k'}};
    for j=1:numel(possible_labels)
        labels=possible_labels{j};
        [has_label,index]=has_fdim_label(ds, labels);
        if has_label
            return
        end
    end

    error(['Unable to find dimension labels, either ''pos'' '...
                    'or ''i'', ''j'', and ''k''']);

function fdim=get_spherical_fdim(ds, target_labels)
    first_target_label=target_labels{1};
    [two, index]=cosmo_dim_find(ds,first_target_label,true);

    if two~=2
        error('dimension ''%s'' must be a feature dimension');
    end
    cosmo_isfield(ds,'a.fdim.labels',true);

    dim_labels=ds.a.fdim.labels(:);
    dim_values=ds.a.fdim.values(:);

    nlabels=numel(target_labels);
    idx_labels=(index+(0:(nlabels-1)))';
    if numel(dim_labels)<index+(nlabels-1) || ...
            ~isequal(dim_labels(idx_labels),target_labels)
        error('expected labels %s in .a.fdim.labels(%d:%d)',...
                  cosmo_strjoin(target_labels,', '),...
                idx_labels(1), idx_labels(end));
    end

    fdim=struct();
    fdim.labels=dim_labels(idx_labels);
    fdim.values=dim_values(idx_labels);

    fdim=ensure_row_vector_or_3d_matrix(fdim);

function fdim=ensure_row_vector_or_3d_matrix(fdim)
    labels=fdim.labels;
    nlabels=numel(labels);

    keys={'labels','values'};
    nkeys=numel(keys);
    for k=1:nlabels
        label=labels{k};
        for j=1:nkeys
            key=keys{j};
            value=fdim.(key){k};

            if strcmp(label,'pos') && strcmp(key,'values')
                if size(value,1)~=3
                    error(['''pos'' attribute in .a.fdim.values '...
                                'must be 3xM']);
                end
            else
                if ~isvector(value)
                    error(['''%s'' attribute in .a.fdim.%s must '...
                            'be a vector'],labels,key);
                end
                fdim.(key){k}=value(:)';
            end
        end
    end

function fa=get_spherical_fa(ds_fa, target_labels)
    fa_cell=get_spherical_fa_cell(ds_fa, target_labels);
    fa=cell2struct(fa_cell,target_labels,2);

function fa_cell=get_spherical_fa_cell(ds_fa, target_labels)
    nlabels=numel(target_labels);
    fa_cell=cell(1,nlabels);
    for j=1:nlabels
        label=target_labels{j};
        fa_cell{j}=ds_fa.(label);
    end

function [sphere_offsets, o_distances]=get_sphere_offsets(radius)
    % return offsets and euclidean (and a bit manhattan) distance
    % from origin
    [sphere_offsets, norm2_distances]=cosmo_sphere_offsets(radius);

    % compute manhattan distance
    norm1_distances=sum(abs(sphere_offsets),2);

    % add a tiny bit of manhattan to make distances more varied
    norm12_distances=norm2_distances+1e-5*norm1_distances;

    % ensure distances are sorted
    [o_distances,i]=sort(norm12_distances);
    sphere_offsets=sphere_offsets(i,:);

function check_input(varargin)
    if numel(varargin)<1 || isscalar(varargin{1})
        % change in parameters
        raise_parameter_error();
    end

function [use_fixed_radius,radius,voxel_count]=get_selection_params(opt)
    if isfield(opt,'radius')
        if isfield(opt,'count')
            raise_parameter_error();
        elseif isscalar(opt.radius) && opt.radius>=0
            use_fixed_radius=true;
            radius=opt.radius;
            voxel_count=NaN;
            return
        end
    elseif isfield(opt,'count') && isscalar(opt.count) && ...
                opt.count>=0 && round(opt.count)==opt.count
        use_fixed_radius=false;
        radius=1; % starting point
        voxel_count=opt.count;
        return;
    end

    raise_parameter_error();

function raise_parameter_error()
    name=mfilename();
    error(['Illegal parameters, use one of:\n',...
        '- %s(...,''radius'',r) to use a radius of r voxels\n',...
        '- %s(...,''count'',c) to select c voxels per searchlight\n',...
        '(As of January 2014 the syntax of this function has changed)'],...
            name,name);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [v_res, solution] = MTA_MIQP(OptimizationModel, KOrxn, varargin)
% Returns the OptimizationModel solution of a particular MTA problem and
% an specific model
%
% USAGE:
%
%    [v_res, success, unsuccess] = MTA_MIQP (OptimizationModel, KOrxn, numWorkers, timeLimit, printLevel)
%
% INPUT:
%    OptimizationModel:    Cplex Model struct
%    KOrxn:                perturbation in the model (reactions)
%    numWorkers:           number of threads used by Cplex.
%    FORCE_CPLEX:          1 to force CPLEX solver, 0 (default) for COBRA
%                          solver.
%    printLevel:           1 if the process is wanted to be shown on the
%                          screen, 0 otherwise. Default: 1.
%
% OUTPUTS:
%    Vout:                 Solution flux of MIQP formulation for each case
%    solution:             Cplex solution struct
%
% .. Authors:
%       - Luis V. Valcarcel, 03/06/2015, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Luis V. Valcarcel, 26/10/2018, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Luis V. Valcarcel, 09/03/2021, University of Navarra, CIMA & TECNUN School of Engineering.

p = inputParser; % check the input information
% check requiered arguments
addRequired(p, 'OptimizationModel');
addRequired(p, 'KOrxn');
% Check optional arguments
addParameter(p, 'numWorkers', 0,@(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'timeLimit', inf,@(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'printLevel', 1,@(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'FORCE_CPLEX', 0,@(x)isnumeric(x)&&isscalar(x));
% extract variables from parser
parse(p, OptimizationModel, KOrxn, varargin{:});
numWorkers = p.Results.numWorkers;
timeLimit = p.Results.timeLimit;
printLevel = max(p.Results.printLevel, 0);
FORCE_CPLEX = p.Results.FORCE_CPLEX;

%Indexation of variables
v = OptimizationModel.idx_variables.v;
y_plus_F = OptimizationModel.idx_variables.y_plus_F;
y_minus_F = OptimizationModel.idx_variables.y_minus_F;
y_plus_B = OptimizationModel.idx_variables.y_plus_B;
y_minus_B = OptimizationModel.idx_variables.y_minus_B;
OptimizationModel = rmfield(OptimizationModel,'idx_variables');

% Temporal way: use ibm_cplex if installed until MIQP API for COBRA is
% implemented
global SOLVERS;
global CBT_MIQP_SOLVER
if FORCE_CPLEX || (SOLVERS.ibm_cplex.installed && strcmp(CBT_MIQP_SOLVER,'ibm_cplex'))
    % Generate CPLEX model
    cplex = Cplex('MIQP');
    CplexModel = OptimizationModel;

    b_L(CplexModel.csense == 'E') = CplexModel.b(CplexModel.csense == 'E');
    b_U(CplexModel.csense == 'E') = CplexModel.b(CplexModel.csense == 'E');
    b_L(CplexModel.csense == 'G') = CplexModel.b(CplexModel.csense == 'G');
    b_U(CplexModel.csense == 'G') = inf;
    b_L(CplexModel.csense == 'L') = -inf;
    b_U(CplexModel.csense == 'L') = CplexModel.b(CplexModel.csense == 'L');
    CplexModel.rhs = b_U;
    CplexModel.lhs = b_L;
    CplexModel.Q = CplexModel.F;
    CplexModel.obj = CplexModel.c;
    CplexModel.ctype = CplexModel.vartype;
    CplexModel.sense = 'minimize';

    cplex.Model = CplexModel;
    % include the knock-out reactions
    cplex.Model.lb(KOrxn) = 0;
    cplex.Model.ub(KOrxn) = 0;

    % Cplex Parameter
    if numWorkers>0
        cplex.Param.threads.Cur = numWorkers;
    end
    if printLevel <=1
        cplex.Param.output.clonelog.Cur = 0;
        cplex.DisplayFunc = [];
    elseif printLevel <=2
        cplex.Param.output.clonelog.Cur = 0;
    end
    if timeLimit < 1e75
        cplex.Param.timelimit.Cur = timeLimit;
    end
    %reduce the tolerance
    cplex.Param.mip.tolerances.mipgap.Cur = 1e-5;
    % cplex.Param.mip.tolerances.absmipgap.Cur = 1e-8;
    % cplex.Param.threads.Cur = 16;

    % SOLVE the CPLEX problem if not singular
    try
        cplex.solve();
    catch
        v_res = zeros(length(v),1);
        return
    end

    if cplex.Solution.status ~= 103
        v_res = cplex.Solution.x(v);
        solution = cplex.Solution;
    else
        v_res = zeros(length(v),1);
        solution = nan;
    end

    % clear the cplex object
    delete(cplex)
    clear cplex
else
    % Generate OptimizationModel for this iteration
    MIQPproblem = OptimizationModel;
    % include the knock-out reactions
    MIQPproblem.lb(KOrxn) = 0;
    MIQPproblem.ub(KOrxn) = 0;

    % Solver Parameter
    if timeLimit > 1e75
        timeLimit = 1e75;
    end

    % SOLVE the MIQP problem
    solution = solveCobraMIQP(MIQPproblem, ...
        'timeLimit',timeLimit, 'relMipGapTol',  1e-5, ...
        'printLevel', max(printLevel-1,0), 'logFile', 0,...
        'threads',numWorkers);

    if isnumeric(solution.stat) && solution.stat == 1
        v_res = solution.full(v);
    elseif ischar(solution.stat) && strcmp(solution.stat, 'OPTIMAL')
        v_res = solution.full(v);
    else
        % Use of try for different outputs of COBRA MIQP solver
        try
            v_res = solution.full(v);
        catch
            v_res = zeros(length(v),1);
        end
    end

end

author: opencobra

repo: cobratoolbox


text:
function [alphas,l,rho,d12,ix] = matRad_siddonRayTracer(isocenter, ...
                                    resolution, ...
                                    sourcePoint, ...
                                    targetPoint, ...
                                    cubes)
% siddon ray tracing through 3D cube to calculate the radiological depth
% according to Siddon 1985 Medical Physics
%
% call
%   [alphas,l,rho,d12,vis] = matRad_siddonRayTracer(isocenter, ...
%                               resolution, ...
%                               sourcePoint, ...
%                               targetPoint, ...
%                               cubes)
%
% input
%   isocenter:      isocenter within cube [voxels]
%   resolution:     resolution of the cubes [mm/voxel]
%   sourcePoint:    source point of ray tracing
%   targetPoint:    target point of ray tracing
%   cubes:          cell array of cubes for ray tracing (it is possible to pass
%                   multiple cubes for ray tracing to save computation time)
%
% output (see Siddon 1985 Medical Physics for a detailed description of the
% variales)
%   alphas          relative distance between start and endpoint for the
%                    intersections with the cube
%   l               lengths of intersestions with cubes
%   rho             densities extracted from cubes
%   d12             distance between start and endpoint of ray tracing
%   ix              indices of hit voxels
%
% References
%   [1] http://www.ncbi.nlm.nih.gov/pubmed/4000088
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2015 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Add isocenter to source and target point. Because the algorithm does not
% works with negatives values. This put (resolution.x,resolution.y,resolution.z)
% in the center of first voxel

sourcePoint = sourcePoint + isocenter;
targetPoint = targetPoint + isocenter;

% Save the numbers of planes.
[yNumPlanes, xNumPlanes, zNumPlanes] = size(cubes{1});
xNumPlanes = xNumPlanes + 1;
yNumPlanes = yNumPlanes + 1;
zNumPlanes = zNumPlanes + 1;

% eq 11
% Calculate the distance from source to target point.
d12 = norm(sourcePoint-targetPoint);

% eq 3
% Position of first planes in millimeter. 0.5 because the central position
% of the first voxel is at [resolution.x resolution.y resolution.z]
xPlane_1 = .5*resolution.x;
yPlane_1 = .5*resolution.y;
zPlane_1 = .5*resolution.z;

% Position of last planes in milimiter
xPlane_end = (xNumPlanes - .5)*resolution.x;
yPlane_end = (yNumPlanes - .5)*resolution.y;
zPlane_end = (zNumPlanes - .5)*resolution.z;

% eq 4
% Calculate parametrics values of \alpha_{min} and \alpha_{max} for every
% axis, intersecting the ray with the sides of the CT.
if targetPoint(1) ~= sourcePoint(1)
    aX_1   = (xPlane_1 - sourcePoint(1)) / (targetPoint(1) - sourcePoint(1));
    aX_end = (xPlane_end - sourcePoint(1)) / (targetPoint(1) - sourcePoint(1));
else
    aX_1   = [];
    aX_end = [];
end
if targetPoint(2) ~= sourcePoint(2)
    aY_1   = (yPlane_1 - sourcePoint(2)) / (targetPoint(2) - sourcePoint(2));
    aY_end = (yPlane_end - sourcePoint(2)) / (targetPoint(2) - sourcePoint(2));
else
    aY_1   = [];
    aY_end = [];
end
if targetPoint(3) ~= sourcePoint(3)
    aZ_1   = (zPlane_1 - sourcePoint(3)) / (targetPoint(3) - sourcePoint(3));
    aZ_end = (zPlane_end - sourcePoint(3)) / (targetPoint(3) - sourcePoint(3));
else
    aZ_1   = [];
    aZ_end = [];
end

% eq 5
% Compute the \alpha_{min} and \alpha_{max} in terms of parametric values
% given by equation 4.
alpha_min = max([0 min(aX_1,aX_end) min(aY_1,aY_end) min(aZ_1,aZ_end)]);
alpha_max = min([1 max(aX_1,aX_end) max(aY_1,aY_end) max(aZ_1,aZ_end)]);

% eq 6
% Calculate the range of indeces who gives parametric values for
% intersected planes.
if targetPoint(1) == sourcePoint(1)
    i_min = []; i_max = [];
elseif targetPoint(1) > sourcePoint(1)
    i_min = xNumPlanes - (xPlane_end - alpha_min * (targetPoint(1) - sourcePoint(1)) - sourcePoint(1))/resolution.x;
    i_max = 1          + (sourcePoint(1) + alpha_max * (targetPoint(1) - sourcePoint(1)) - xPlane_1)/resolution.x;
    % rounding
    i_min = ceil(1/1000*(round(1000*i_min)));
    i_max = floor(1/1000*(round(1000*i_max)));
else
    i_min = xNumPlanes - (xPlane_end - alpha_max * (targetPoint(1) - sourcePoint(1)) - sourcePoint(1))/resolution.x;
    i_max = 1          + (sourcePoint(1) + alpha_min * (targetPoint(1) - sourcePoint(1)) - xPlane_1)/resolution.x;
    i_min = ceil(1/1000*(round(1000*i_min)));
    i_max = floor(1/1000*(round(1000*i_max)));
end
if targetPoint(2) == sourcePoint(2)
    j_min = []; j_max = [];
elseif targetPoint(2) > sourcePoint(2)
    j_min = yNumPlanes - (yPlane_end - alpha_min * (targetPoint(2) - sourcePoint(2)) - sourcePoint(2))/resolution.y;
    j_max = 1          + (sourcePoint(2) + alpha_max * (targetPoint(2) - sourcePoint(2)) - yPlane_1)/resolution.y;
    j_min = ceil(1/1000*(round(1000*j_min)));
    j_max = floor(1/1000*(round(1000*j_max)));
else
    j_min = yNumPlanes - (yPlane_end - alpha_max * (targetPoint(2) - sourcePoint(2)) - sourcePoint(2))/resolution.y;
    j_max = 1          + (sourcePoint(2) + alpha_min * (targetPoint(2) - sourcePoint(2)) - yPlane_1)/resolution.y;
    j_min = ceil(1/1000*(round(1000*j_min)));
    j_max = floor(1/1000*(round(1000*j_max)));
end
if targetPoint(3) == sourcePoint(3)
    k_min = []; k_max = [];
elseif targetPoint(3) >= sourcePoint(3)
    k_min = zNumPlanes - (zPlane_end - alpha_min * (targetPoint(3) - sourcePoint(3)) - sourcePoint(3))/resolution.z;
    k_max = 1          + (sourcePoint(3) + alpha_max * (targetPoint(3) - sourcePoint(3)) - zPlane_1)/resolution.z;
    k_min = ceil(1/1000*(round(1000*k_min)));
    k_max = floor(1/1000*(round(1000*k_max)));
else
    k_min = zNumPlanes - (zPlane_end - alpha_max * (targetPoint(3) - sourcePoint(3)) - sourcePoint(3))/resolution.z;
    k_max = 1          + (sourcePoint(3) + alpha_min * (targetPoint(3) - sourcePoint(3)) - zPlane_1)/resolution.z;
    k_min = ceil(1/1000*(round(1000*k_min)));
    k_max = floor(1/1000*(round(1000*k_max)));
end

% eq 7
% For the given range of indices, calculate the paremetrics values who
% represents intersections of the ray with the plane.
if i_min ~= i_max
    if targetPoint(1) > sourcePoint(1)
        alpha_x = (resolution.x*(i_min:1:i_max)-sourcePoint(1)-.5*resolution.x)/(targetPoint(1)-sourcePoint(1));
    else
        alpha_x = (resolution.x*(i_max:-1:i_min)-sourcePoint(1)-.5*resolution.x)/(targetPoint(1)-sourcePoint(1));
    end
else
    alpha_x = [];
end
if j_min ~= j_max
    if targetPoint(2) > sourcePoint(2)
        alpha_y = (resolution.y*(j_min:1:j_max)-sourcePoint(2)-.5*resolution.y)/(targetPoint(2)-sourcePoint(2));
    else
        alpha_y = (resolution.y*(j_max:-1:j_min)-sourcePoint(2)-.5*resolution.y)/(targetPoint(2)-sourcePoint(2));
    end
else
    alpha_y = [];
end
if k_min ~= k_max
    if targetPoint(3) > sourcePoint(3)
        alpha_z = (resolution.z*(k_min:1:k_max)-sourcePoint(3)-.5*resolution.z)/(targetPoint(3)-sourcePoint(3));
    else
        alpha_z = (resolution.z*(k_max:-1:k_min)-sourcePoint(3)-.5*resolution.z)/(targetPoint(3)-sourcePoint(3));
    end
else
    alpha_z = [];
end

% eq 8
% Merge parametrics sets.
alphas = unique([alpha_min alpha_x alpha_y alpha_z alpha_max]);

% eq 10
% Calculate the voxel intersection length.
l = d12*diff(alphas);

% eq 13
% Calculate \alpha_{middle}
alphas_mid = .5*(alphas(1:end-1)+alphas(2:end));

% eq 12
% Calculate the voxel indices: first convert to physical coords
i_mm = sourcePoint(1) + alphas_mid*(targetPoint(1) - sourcePoint(1));
j_mm = sourcePoint(2) + alphas_mid*(targetPoint(2) - sourcePoint(2));
k_mm = sourcePoint(3) + alphas_mid*(targetPoint(3) - sourcePoint(3));
% then convert to voxel index
i = round(i_mm/resolution.x);
j = round(j_mm/resolution.y);
k = round(k_mm/resolution.z);

% Handle numerical instabilities at the borders.
i(i<=0) = 1; j(j<=0) = 1; k(k<=0) = 1;
i(i>xNumPlanes-1) = xNumPlanes-1;
j(j>yNumPlanes-1) = yNumPlanes-1;
k(k>zNumPlanes-1) = zNumPlanes-1;

% Convert to linear indices
ix = j + (i-1)*size(cubes{1},1) + (k-1)*size(cubes{1},1)*size(cubes{1},2);

% obtains the values from cubes
rho = cell(numel(cubes),1);
for i = 1:numel(cubes)
    rho{i} = cubes{i}(ix);
end

repo: matRad


text:
% BUZCODE
%
% Files
%   bz_IsEventDatatype - this function checks that the input variable is of the buzcode event
%   compileBuzcode     - this function recompiles all .c scripts to create the apropriate .mex* files on your local machine

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function sparse_test01 ( )

%*****************************************************************************80
%
%% SPARSE_TEST01 demonstrates a simple use of the SPARSE matrix facility.
%
%  Discussion:
%
%    This unrealistic example sets up a tiny matrix:
%
%    11  12   0   0  15
%     0  22  23   0   0
%    31   0  33  34  35
%
%    Note that we assume here that the entire matrix is available
%    at once, so we only need a single call to SPARSE to define the matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield and John Penny,
%    Numerical Methods Using MATLAB,
%    Prentice Hall, 1999
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST01:\n' );
  fprintf ( 1, '  Demonstrate the use of MATLAB''s SPARSE facility\n' );
  fprintf ( 1, '  to define a tiny sparse matrix.\n' );

  i = [  1,  1,  1,  2,  2,  3,  3,  3,  3 ];
  j = [  1,  2,  5,  2,  3,  1,  3,  4,  5 ];
  s = [ 11, 12, 15, 22, 23, 31, 33, 34, 35 ];
  m = 3;
  n = 5;
  nzmax = 9;

  a = sparse ( i, j, s, m, n, nzmax );

  a

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = line_cvt_lloyd ( n, a, b, it_num, header, x )

%*****************************************************************************80
%
%% LINE_CVT_LLOYD carries out the Lloyd algorithm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%
%    Input, real A, B, the left and right endpoints.
%
%    Input, integer IT_NUM, the number of iterations to take.
%
%    Input, string HEADER, an identifying string.
%
%    Input, real X(N), the initial point locations.
%
%    Output, real X(N), the final point locations.
%
  x = x(:);
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %f\n', k, x(k,1) );
  end
%
%  Initialize the plotting arrays.
%
  step = 1 : it_num;
  e = nan ( it_num, 1 );
  xm = nan ( it_num, 1 );

  for it = 1 : it_num

    x_plot ( 1:n, it ) = x(1:n,1);

    x_new = line_cvt_lloyd_step ( n, a, b, x );

    e(it) = line_cvt_energy ( n, a, b, x );
    e(it) = max ( e(it), eps );
%
%  Display the energy.
%
    figure ( 1 )
    plot ( step, log ( e ), 'm-*' )
    title ( 'Log (Energy)' )
    xlabel ( 'Step' )
    ylabel ( 'Energy' )
    grid
%
%  Compute the generator motion.
%
    xm(it,1) = sum ( ( x_new(:) - x(:) ).^2 ) / n;
%
%  Display the generator motion.
%
    figure ( 2 )
    plot ( step, log ( xm ), 'm-*' )
    title ( 'Log (Average generator motion)' )
    xlabel ( 'Step' )
    ylabel ( 'Motion' )
    grid
%
%  Update the generators.
%
    x(1:n,1) = x_new(1:n,1);

  end

  x_plot(1:n,it_num+1) = x(1:n,1);
%
%  Print the current generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %f\n', k, x(k,1) );
  end
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 3 )

  y = ( 0 : it_num );
  for k = 1 : n
    plot ( x_plot(k,1:it_num+1), y )
    hold on;
  end
  grid on
  hold off;

  title ( 'Generator evolution.' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );
%
%  Save the plots.
%
  figure ( 1 )
  filename = strcat ( header, '_energy.png' );
  print ( '-dpng', filename );
  figure ( 2 )
  filename = strcat ( header, '_motion.png' );
  print ( '-dpng', filename );
  figure ( 3 )
  filename = strcat ( header, '_evolution.png' );
  print ( '-dpng', filename );

  return
end


text:
 function centre=computeCOR(data,geo,angles,slice)
% Code modified from SophiaBeads
% (https://github.com/Sophilyplum/sophiabeads-datasets/blob/master/tools/centre_geom.m)
% Reference:
% T. Liu - "Direct central ray determination in computed microtomography",
% Optical Engineering, April 2009.

% Get central slice
if nargin<4
    slice=floor(size(data,1)/2)+1;
end
data=squeeze(data(slice,:,:));

if size(angles,1)==1
   angles=angles';
end

% Set up coordinate grids for testing the fit to data
[angle_grid, det_grid] = meshgrid(angles',linspace(-geo.sDetector(1)/2+geo.dDetector(1)/2,+geo.sDetector(1)/2-geo.dDetector(1)/2,geo.nDetector(1))');
% Wrap grids in the angular direction to avoid value out of range errors
angle_grid = [(angle_grid - 2*pi) angle_grid (angle_grid + 2*pi)];
det_grid = repmat(det_grid,1,3);
test_data = double(repmat(data,1,3));

% Start search using midpoint at zero
midpoint = 0;
% Vector of precision values to search at
precision = [1 0.1 0.01 0.001 0.0001 0.00001 0.000001];

for i = 1:length(precision)

    COR = (midpoint - 10*precision(i)):precision(i):(midpoint + 10*precision(i)); % values for centre of rotation
    M = zeros(length(COR),1);

    for j = 1:length(COR)

        gamma = atan(linspace(-geo.sDetector(1)/2+geo.dDetector(1)/2,+geo.sDetector(1)/2-geo.dDetector(1)/2,geo.nDetector(1))' / geo.DSD);   % angle of each ray relative to theoretical central ray
        gamma_c = atan(COR(j) / geo.DSD); % angle of assumed centre of rotation to central ray
        gamma_i = gamma - gamma_c;
        beta = 2 * gamma_i + pi;

        s2 = geo.DSD * tan(2 * gamma_c - gamma);
        s2 = repmat(s2, 1, size(angles,2));

        angles_aux = repmat(angles', geo.nDetector(1), 1) + repmat(beta, 1, size(angles,2));
        test = interp2(angle_grid, det_grid, test_data, angles_aux, s2, 'linear', 0);

        nonzero = find(test > 0);
        % We want the number of non-zero values for the average, not the sum of their positions
        M(j) = sum((test(nonzero) - data(nonzero)).^2)*(1/length(nonzero));
    end

    [~, indM] = min(M);   % minimum value and index
    midpoint = COR(indM);   % set midpoint for next search
end

% transform centre to required value
centre =- midpoint * geo.DSO / geo.DSD;

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function excelrecord3fcn(stringdates1, endo, Y, ss_estimates, n, pref)
% script excelrecord3
% records the information contained in the worksheet 'steady state' of the excel spreadsheet 'results.xls'

% create the cell that will be saved on excel
sscell={};

% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates1,1)+3,1);

% loop over variables (horizontal dimension)
for ii=1:n
   % create cell of steady-state record for variable ii
   temp=['steady-state and actual: ' endo{ii,1}];
   ss_i=[temp {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'sample'} {'lw. bound'} {'median'} {'up. bound'};stringdates1 num2cell(Y(:,ii)) num2cell(ss_estimates{ii,1}')];
sscell=[sscell ss_i vertspace];
end

% trim
sscell=sscell(:,1:end-1);

% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),sscell,'steady state','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function imgOut = pyrBlend(img1, img2, weight)
%
%
%        imgOut = pyrBlend(img1, img2, weight)
%
%
%        Input:
%           -img1: an image to be blended
%           -img2: an image to be blended
%           -weight: the weights for img1
%
%        Output:
%           -imgOut: the final blended image
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~isSameImage(img1, img2) || ~isSimilarImage(img1, weight))
   error('pyrBlend: input images are different!');
end

col = size(img1, 3);

p1 = pyrImg3(img1, @pyrLapGen);
p2 = pyrImg3(img2, @pyrLapGen);

g1 = pyrGaussGen(weight);
g2 = pyrGaussGen(1 - weight);

imgOut = zeros(size(img1));

for i=1:col
    tpg1 = pyrLst2OP(p1(i), g1,  @pyrMul);
    tpg2 = pyrLst2OP(p2(i), g2,  @pyrMul);
    tf   = pyrLst2OP(tpg1, tpg2, @pyrAdd);

    imgOut(:,:,i) = pyrVal(tf);
end

end

author: banterle


text:
function rxToggle(rx, nToggles, pauseForSecs);
% For mrRx: quickly toggle between a prescribed slice image and
% a reference slice for comparison.
%
% rxToggle(rx, [nToggles=3]);
%
%
% ras, 01/06.
% ras, 05/09: added nToggles as an input parameter.
if ~exist('rx', 'var') | isempty(rx),
    cfig = findobj('Tag', 'rxControlFig');
    rx = get(cfig, 'UserData');
end

if notDefined('pauseForSecs')
    pauseForSecs = 0.15;
end

if notDefined('nToggles')
	% # of times to alternate b/w interp and ref images
	nToggles = 3;
end

[interp ref] = rxGetComparisonImages(rx);

if ishandle(rx.ui.compareAxes)
    axes(rx.ui.compareAxes)
else
    figure('Name', 'mrRx: rxToggle');
end

rng = [min(interp(:)) max(interp(:))];

% get cur axis limits and keep them, in case user has zoomed image
x = get(gca, 'xlim');
y = get(gca, 'ylim');

for n = 1:nToggles
    imshow(ref, rng);
    xlim(x); ylim(y);
    pause(pauseForSecs);
    imshow(interp, rng);
    xlim(x); ylim(y);
    pause(pauseForSecs);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% PRIME_POOL uses the MATLABPOOL command to run the PRIME code.
%
%  Discussion:
%
%    Output printed by the function appears directly on the screen.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
  n = 10000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_POOL\n' );
  fprintf ( 1, '  Count prime numbers from 1 to %d\n', n );

  matlabpool open local 4

  total = prime_fun ( n );

  matlabpool close

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Total is %d\n', total );


text:
%% HILBERT_BATCH_LOCAL uses the BATCH command to run the HILBERT code locally.
%
%  Discussion:
%
%    The HILBERT code is a function, so first we must write a script
%    called HILBERT_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2012
%
%  Author:
%
%    John Burkardt
%
  clear

  workers = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HILBERT_BATCH_LOCAL\n' );
  fprintf ( 1, '  Run HILBERT_SCRIPT locally.\n' );
%
%  The BATCH command sends the script for execution.
%
  my_job = batch ( 'hilbert_script', ...
    'Configuration', 'local', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'hilbert_fun' }, ...
    'matlabpool', workers );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  The LOAD command makes the script's workspace available.
%  The script created the array "a", so let's ask for its size
%  just to show that it really was created.
%
  load ( my_job );

  size ( a )
%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

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


text:
function prob_test095 ( )

%*****************************************************************************80
%
%% TEST095 tests LAPLACE_CDF, LAPLACE_CDF_INV, LAPLACE_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST095\n' );
  fprintf ( 1, '  For the Laplace PDF:\n' );
  fprintf ( 1, '  LAPLACE_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  LAPLACE_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  LAPLACE_PDF evaluates the PDF;\n' );

  a = 1.0;
  b = 2.0;

  check = laplace_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST095 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = laplace_sample ( a, b, seed );

    pdf = laplace_pdf ( x, a, b );

    cdf = laplace_cdf ( x, a, b );

    x2 = laplace_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function [m,a]=lpcstable(ar)
%LPCSTABLE Test AR coefficients for stability and stabilize if necessary [MA,A]=(AR)
%
% Usage: (1) [m,ar]=lpcstable(ar); % force ar polynolials to be stable
%
%   Input:  ar(:,p+1)  Autoregressive coefficients
% Outputs:  m(:,1)    Mask identifying stable polynomials
%           a(:,p+1)  Stabilized polynomials formed by reflecting unstable
%                       poles in unit circle (with a(:,1)=1)

%      Copyright (C) Mike Brookes 2016
%      Version: $Id: lpcstable.m 8558 2016-09-22 08:22:54Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1] = size(ar);
mm=ar(:,1)~=1;
if any(mm)          % first ensure leading coefficient is always 1
    ar(mm,:)=ar(mm,:)./ar(mm,ones(1,p1));
end
if p1==1
    m=ones(nf,1); % 0'th order filter is always stable
elseif p1==2
    m=abs(ar(:,2))<1; % 1'st order filter
else
    rf = ar;
    k = rf(:,p1); % check final coefficient in range
    m=abs(k)<1;
    if any(m)
        d = (1-k(m).^2).^(-1);
        wj=ones(1,p1-2);
        rf(m,2:p1-1) = (rf(m,2:p1-1)-k(m,wj).*rf(m,p1-1:-1:2)).*d(:,wj);
        for j = p1-2:-1:2
            k(m) = rf(m,j+1);
            m(m)=abs(k)<1;
            if ~any(m), break, end
            d = (1-k(m).^2).^(-1);
            wj=ones(1,j-1);
            rf(m,2:j) = (rf(m,2:j)-k(m,wj).*rf(m,j:-1:2)).*d(:,wj);
        end
    end
end
if nargout>1
    a=ar;
    if ~all(m)
        for i=find(~m)'                 % unstable frames
            z=roots(a(i,:));
            k=abs(z)>1;                 % find any unstable roots
            z(k)=conj(z(k)).^(-1);      % invert them
            a(i,:)=real(poly(z));       % force a real polynomial
        end
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function problem = QUASAR_Problem(v1,v2,barc2,varargin)
% return the cost function and constraint matrices of QUASAR

fprintf('Relax Wahba problem to standard linear SDP...')

params = inputParser;
params.CaseSensitive = false;
params.addParameter('computeConstraintMatrix', 'true', @(x) islogical(x));
params.addParameter('TraceAll','true',@(x) islogical(x));
params.parse(varargin{:});

computeConstraintMatrix = params.Results.computeConstraintMatrix;
TraceAll                = params.Results.TraceAll;

problem.computeConstraintMatrix = computeConstraintMatrix;

N = size(v1,2); % nr vector
% X = [q\tran, q_1\tran, q_2\tran, ..., q_N\tran] has size 1 by Npm, where Npm = a+b+4*N
Npm=4+4*N;

% coefficient matrix that maps vec(qq\tran) to vec(R)
P=[1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1;
   0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0;
   0, 0, 1, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, -1, 0, 0;
   0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 0;
   -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1;
   0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0;
   0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0;
   0, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0, -1, 0, 0, 0;
   -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1];
P=sparse(P);

% build the cost matrix Q_cost
Q_1=zeros(Npm,Npm);
for k=1:N
    idx = 4+blkIndices(k,4);
    P_k = reshape(P'*reshape(v2(:,k)*v1(:,k)',[9,1]),[4,4]);
    ck = 0.5 * ( v1(:,k)'*v1(:,k)+v2(:,k)'*v2(:,k) - barc2 );
    Q_1((1:4),idx) = Q_1((1:4),idx)-0.5*P_k+ck/2*eye(4);
    Q_1(idx,(1:4)) = Q_1(idx,(1:4))-0.5*P_k+ck/2*eye(4);
end

Q_2=zeros(Npm,Npm);
for k=1:N
%     idx = 4+blkIndices(k,4);
    idx = blkIndices(1,4);
    P_k = reshape(P'*reshape(v2(:,k)*v1(:,k)',[9,1]),[4,4]);
    ck = 0.5 * ( v1(:,k)'*v1(:,k)+v2(:,k)'*v2(:,k) + barc2 );
    Q_2(idx,idx) = Q_2(idx,idx) - P_k + ck*eye(4);
end
Q_cost=Q_1+Q_2;
Q_cost=sparse(Q_cost);

Q_cost=Q_cost/barc2;

problem.P = P;
problem.C = {Q_cost};
problem.n = size(Q_cost,1);
n = problem.n;

if computeConstraintMatrix
    % compute all the constraint matrices
    nr_constraints = 3*N*(N+1) + 10*N + 1;
    % nr_licq = 4*N+1;
    nr_licq = 10*N+1;
    nr_redundant = nr_constraints - nr_licq;
    b = sparse([1],[1],[1],nr_constraints,1);
    A = {};
    % trace([Z]_qq) = 1
    if TraceAll
        A_trace  = speye(n,n);
        b(1)     = N + 1;
    else
        A_trace = sparse([1,2,3,4],[1,2,3,4],[1,1,1,1],n,n);
        b(1) = 1;
    end
    A{end+1} = A_trace;

    % [Z]_qiqi - [Z]_qq = 0
    for k = 1:N
        for i = 1:4
            for j = i:4
                row_shift = 4*k;
                col_shift = 4*k;
                tmp = sparse([i,i+row_shift],[j,j+col_shift],[-1,1],n,n);
                tmp = tmp + tmp'; % make it symmetric
                tmp = tmp / norm(tmp,'fro'); % normalize to frobenius norm=1
                A{end+1} = tmp;
            end
        end
    end

    % [Z]_qiqj symmetric
    fprintf('\nBuild symmetric constraints, progress ... ')
    for k1 = 1:N
        fprintf('%d/%d ',k1,N);
        for k2 = k1+1:N+1
            for i = 1:3
                for j = i+1:4
                    row_shift = 4*(k1-1);
                    col_shift = 4*(k2-1);
                    tmp = sparse([i+row_shift,j+row_shift],[j+col_shift,i+col_shift],[-1,1],n,n);
                    tmp = tmp + tmp'; % make it symmetric
                    tmp = tmp / 2; % normalize to frobenius norm=1
                    A{end+1} = tmp;
                end
            end
        end
    end
    fprintf('Done.\n')
    %{
%*************************************************************
%% added by Kim-Chuan Toh
%*************************************************************
    cnt = 1;
    len = 3*N*(N+1); rr = 1/sqrt(2);
    row = zeros(len,1); col=zeros(len,1); val=zeros(len,1);
    fprintf('Build symmetric constraints, progress ... ')
    for k1 = 1:N
        fprintf('%d/%d ',k1,N);
        for k2 = k1+1:N+1
            for i = 1:3
                for j = i+1:4
                    row_shift = 4*(k1-1);
                    col_shift = 4*(k2-1);
                    ii = i+row_shift; ii2 = j+row_shift;
                    jj = j+col_shift; jj2 = i+col_shift;
                    idx = 2*cnt-1:2*cnt;
                    row(idx) = [ii+(jj-1)*jj/2; ii2+(jj2-1)*jj2/2];
                    col(idx) = [cnt; cnt];
                    val(idx) = [-rr; rr];
                    cnt = cnt+1;
                end
            end
        end
    end
    problem.Bt = spconvert([row,col,val;n*(n+1)/2,len,0]);
    fprintf('Done.\n')
%*************************************************************
%*************************************************************
    %}
    assert(length(A) == nr_constraints, 'Incorrect number of constraint matrices in A.');
    problem.Acell = A;
    problem.b = b;
    problem.m_localize = nr_licq;
    problem.m_moment = nr_redundant;
    problem.m = length(A);
    problem.n = size(A{1},1);
    ndelta    = triangle_number(n);
    l         = ndelta - problem.m;
    blk{1,1} = 's';
    blk{1,2} = n;
    problem.blk = blk;

    At       = sparsesvec(blk,problem.Acell);
    problem.At = {At};
    problem.M  = N + 1;

    fprintf('Done.\n')
    fprintf('Linear SDP: n = %d, m = %d, m_localize = %d, m_moment = %d, l = %d.\n',...
        problem.n,problem.m,problem.m_localize,problem.m_moment,l);

%     problem      = rmfield(problem,{'P','computeConstraintMatrix','Acell','m_localize','m_moment'});
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
% GP_SAMPLE_PCG - Draw a sample from a conditional multivariate normal
%                 distribution using (preconditioned) conjugate gradient.
%
% The joint covariance matrix K is assumed to be of the form:
%
%   K = [  K1     K1
%          K1   (K1+K2) ]
%
% A sample is drawn from the conditional distribution
%
%   P(X1|X2) = N( X1 | K1*INV(K1+K2)*X2, K1 - K1*INV(K1+K2)*K1 )
%
% using (preconditioned) conjugate gradient method. This can be efficient
% for sparse covariance matrices.
%
% Usage:
%
% X1 = GP_SAMPLE_PCG(X2, K, K1, L1, L2)
% X1 = GP_SAMPLE_PCG(X2, K, K1, L1, L2, M)
%
% The matrices L1 and L2 are such that K1=L1*L1' and K2=L2*L2'. These can
% be found using, for instance, CHOL or LCHOL.
%
% One can also give function handles such that K(X), K1(X), L1(X) and L2(X)
% return the corresponding matrix-vector products.
%
% M is an optional preconditioner for the conjugate gradient method. To
% speed up the method, INV(M) should approximate INV(K1+K2). It can be
% extremely important to use a good preconditioner. One can also give a
% handle to a function which evaluates and returns M\X.
%
% Z1 = L1*RANDN(N1,1)
% Z2 = L2*RANDN(N2,1)
%
% One can also give a matrix S such that
%
%   K = [  K1       K1*S'
%         S*K1   S*(K1+K2)*S' ]
%
% and then a sample is drawn from the conditional Gaussian distribution
% with mean
%
%   E(X1|X2) = K1 * INV(S*(K1+K2)*S') * S*X2,
%
% and covariance
%
%   COV(X1|X2) = K1 - K1 * S' * INV(S*(K1+K2)*S') * S * K1 )
%
% For instance, S can be an identity matrix from which some rows have been
% removed. This would correspond to having missing values in X2. In any
% case, S should be such that matrix-vector products are fast to evaluate,
% thus, preferably sparse. At least for now, S can not be a function
% handle. Note: the preconditioner might become worse when missing values..
%
% Usage:
%
% X1 = GP_SAMPLE_PCG(X2, K, K1, L1, L2, M, S)

% Last modified 2010-10-29
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function z1 = gaussian_rand_conjgradmv(K1, K2, z1, z2, I, varargin)
%function z1 = gaussian_rand_conjgradmv(x2, K, K1, z1, z2, I, varargin)
% To reduce memory usage: x1=z1
%
% Large variables in memory: x2, z1, z2

options = struct( ...
    'maxiter', 100,  ...
    'tol',     1e-6, ...
    'verbose', false);

% Parse arguments
[options, errmsg] = argparse( options, varargin{:} );
error(errmsg);

if nargin < 6 || isempty(I)
  I = true(size(x2));
end

if isnumeric(K1)
  K1 = @(x) K1*x;
end
if isnumeric(K2)
  K2 = @(x) K2*x;
end

% Assume zero means for now.
%mu1 = 0;
%mu2 = 0;

% Observations are subtracted from Z2!!
%z2 = z2 - x2;
%z2 = z2 - x2 + mu2
%z2 = z2 + z1 - x2 + mu2
z2 = conjgradmv(K2, z2, I, ...
                'tol', options.tol, ...
                'maxiter', options.maxiter, ...
                'verbose', options.verbose);
z1 = z1 - K1(z2);
%z1 = mu1 + z1 - K1(z2);

% x1 = z1;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% Figure 6.19: Polynomial fitting
% Section 6.5.3
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 10/03/05
% (a figure is generated)
%
% Given data u_1,...,u_m and v_1,...,v_m in R, the goal is to fit to the
% data a polynomial of the form
% p(u) = x_1 + x_2*u + ... + x_n*u^{n-1}
% i.e. solve the problem:  minimize ||Ax - v||
% where A is the Vandermonde matrix s.t. Aij = u_i^{j-1}
% Two cases are considered: L2-norm and Linfty-norm

% Input data
n=6;
m=40;
randn('state',0);
% generate 50 ponts ui, vi
u = linspace(-1,1,m);
v = 1./(5+40*u.^2) + 0.1*u.^3 + 0.01*randn(1,m);

% LS fit polynomial x_1 + x_2*u + ... + x_n*u^(n-1) to (ui,vi)
fprintf(1,'Computing optimal polynomial in the case of L2-norm...');

A = vander(u');
A = A(:,m-n+[1:n]);     % last n columns of A
x = A\(v');             % coefficients of the polynomial in the following
                        % order: x = [x_n x_(n-1) ... x_2 x_1]'

fprintf(1,'Done! \n');

% L-infty fit
fprintf(1,'Computing optimal polynomial in the case of Linfty-norm...');

cvx_begin quiet
    variable x1(n)
    minimize (norm(A*x1 - v', inf))
cvx_end

fprintf(1,'Done! \n');

% generates 1000 points in  [-1,1]
u2 = linspace(-1.1,1.1,1000);

% evaluate the interpolating polynomial using Horner's method
vpol = x(1)*ones(1,1000);
vpoll1 = x1(1)*ones(1,1000);
for i = 2:n
  vpol = vpol.*u2 + x(i);
  vpoll1 = vpoll1.*u2 + x1(i);
end;

figure
% plot function and interpolating polynomial
plot(u2, vpol,'-', u, v, 'o', u2, vpoll1,'--');
xlabel('u');
ylabel('p(u)');
title('Fitting of data points with two polynomials of degree 5');
legend('L_2 norm','data points','L_{\infty} norm', 'Location','Best');
% print -deps polapprox.eps

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function vx = VBA_vec(X)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% vx = VBA_vec (X)
% transform array X into a column vector
%
% IN:
%   - X: matrix
% OUT:
%   - vx: vector
%
% /////////////////////////////////////////////////////////////////////////

% avoid weirdly sized empty vector
if isempty(X)
    vx = [];
    return;
end

% return vectorized input
vx = full(X(:));

author: MBB-team

repo: VBA-toolbox


text:
function test_sampled_newton()

    clc;
    clear;
    close all;

    %% Set algorithms
    algorithms = {'L-BFGS-BKT','Newton-CHOLESKY','Newton-INEXACT','Subsamp-Newton-Uniform', 'Subsamp-Newton-RNS', 'Subsamp-Newton-LS'};

    % select problem
    problem_type = 'log_reg';
    %problem_type = 'lin_reg';

    lambda = 0.01;

    % prepare datasets
    if strcmp(problem_type, 'lin_reg')

        % generate synthtic data
        % sample data generating for training: y = w1*x1 + w2*x2 + ... * wd*1
        n = 10000;
        d = 10;
        std = 0.25;
        data = linear_regression_data_generator(n, d, std);

        x_train = data.x_train;
        y_train = data.y_train;
        x_test = data.x_test;
        y_test = data.y_test;

        % solution
        w_opt = pinv(x_train * x_train') * x_train * y_train';
        % for intersect
        d = d + 1;

        % define problem definitions
        problem = linear_regression(x_train, y_train, x_test, y_test, lambda);
        f_opt1 = problem.cost(w_opt);
        fprintf('f_opt: %.16e\n', f_opt1);

        w_opt = problem.calc_solution(100, 'lbfgs');
        f_opt2 = problem.cost(w_opt);
        fprintf('f_opt1: %.16e, f_opt2: %.16e, diff: %.16e\n', f_opt1, f_opt2, f_opt2 - f_opt1);

        if f_opt2 < f_opt1
            f_opt = f_opt2;
        else
            f_opt = f_opt1;
        end

    elseif strcmp(problem_type, 'log_reg')

        % read real-world dataset
        [data_y, data_X] = libsvmread('../data/libsvm/a9a');
        x_train = data_X';
        y_train = data_y';
        d = size(x_train,1);
        n = length(y_train);
        lambda = 0.5;

        if d == 0
            return;
        end

        % define problem definitions
        problem = logistic_regression(x_train, y_train, [], [], lambda);

        % calculate f_opt
        w_opt = problem.calc_solution(100, 'lbfgs');
        f_opt = problem.cost(w_opt);
        fprintf('f_opt: %.24e\n', f_opt);

    else
        return;
    end

    %% initialize
    w_init = zeros(d,1);
    batch_size = 1;
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_optgap = 1e-14;
        options.tol_gnorm = 1e-16;
        options.max_iter = 100;
        options.verbose = 1;
        options.f_opt = f_opt;
        options.store_w = false;
        options.permute_on = 1;
        options.f_opt = f_opt;
        options.batch_size = batch_size;

        switch algorithms{alg_idx}

            % Sub-sampled Newton methods
            case {'Subsamp-Newton-Uniform'}

                options.sub_mode = 'Uniform';
                options.subsamp_hess_size = 100*d;
                [w_list{alg_idx}, info_list{alg_idx}] = subsamp_newton(problem, options);

            case {'Subsamp-Newton-RNS'}

                options.sub_mode = 'RNS';
                options.subsamp_hess_size = 20*d;
                [w_list{alg_idx}, info_list{alg_idx}] = subsamp_newton(problem, options);

            case {'Subsamp-Newton-LS'}

                options.sub_mode = 'LS';
                options.subsamp_hess_size = 20*d;
                options.hess_update_freq = 10;
                [w_list{alg_idx}, info_list{alg_idx}] = subsamp_newton(problem, options);

            % Newton methodss
            case {'Newton-STD'}

                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-DAMP'}

                options.sub_mode = 'DAMPED';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-CHOLESKY-BKT'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-INEXACT-BKT'}

                options.sub_mode = 'INEXACT';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-INEXACT'}

                options.sub_mode = 'INEXACT';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            % BFGS variants
            case {'BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                options.mem_size = 50;
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','optimality_gap', algorithms, w_list, info_list);
    display_graph('time','optimality_gap', algorithms, w_list, info_list);
    %display_graph('time','gnorm', algorithms, w_list, info_list);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
clc,clear all,close all;
[GaitDBDir,probes,numprobe]=testData;
disp(['Processing ',GaitDBDir]);

maxDim=600;
MPCADADim=200;%Dimension of MPCA features b4 LDA
Qval=97;%The Q value, percentage of energy kept
maxK=1;%Maximum number of iterations
strGRDB=[GaitDBDir(12:(end-1))];
%%%%%%%%%%Load gallery gait data
setname='Gal';
load([GaitDBDir,setname]);
Is=size(fea3D);numSpl=Is(4);

fea3D_Train = fea3D;gnd_Train = gnd;clear fea3D
%=================MPCALDA==========================
disp(['Training on gallery']);
[tUs, posIdx, TXmean, Wgt, LDAU] = MPCALDA(fea3D_Train,gnd_Train,MPCADADim,Qval,maxK);
if size(LDAU,2)<maxDim,maxDim=size(LDAU,2);end
testDims=21:maxDim;

newfea=ttm(tensor(fea3D_Train-repmat(TXmean,[ones(1,3), numSpl])),tUs,1:3);
clear fea3D_Train
newfea = reshape(newfea.data,size(newfea,1)*size(newfea,2)*size(newfea,3),numSpl)';
galfea = newfea(:,posIdx);
galfea=galfea*LDAU;
galgnd=gnd;

nDim=length(testDims);
AllSeqR1s=zeros(nDim,numprobe);%Rank 1 recognition rate based on matching gait sequences
AllSeqR5s=zeros(nDim,numprobe);%Rank 5 recognition rate based on matching gait sequences
for iprb=1:numprobe
    %%%%%%%%%%Load probe gait data
    setname=['Prb',probes(iprb)];
    disp(['Testing ',setname]);
    load([GaitDBDir,setname]);Is=size(fea3D);numSpl=Is(4);
    newfea=ttm(tensor(fea3D-repmat(TXmean,[ones(1,3), numSpl])),tUs,1:3);clear fea3D
    newfea = reshape(newfea.data,size(newfea,1)*size(newfea,2)*size(newfea,3),numSpl)';
    prbfea = newfea(:,posIdx);
    prbfea=prbfea*LDAU;
    prbgnd=gnd;
    [SeqR1s,SeqR5s]=MADAll(galfea,galgnd,prbfea,prbgnd,testDims,Wgt);
    AllSeqR1s(:,iprb)=SeqR1s;
    AllSeqR5s(:,iprb)=SeqR5s;
end
MaxSeqR1=max(AllSeqR1s);
MaxSeqR5=max(AllSeqR5s);
disp('Rank 1 identification rate for Probe A to G')
disp(MaxSeqR1)
disp('Rank 5 identification rate for Probe A to G')
disp(MaxSeqR5)
disp(['Mean rank 1 identification rate=',num2str(mean(MaxSeqR1))])
disp(['Mean rank 5 identification rate=',num2str(mean(MaxSeqR5))])

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Start by maximum global manipulability at a pose
% Start by maximizing manipulability locally at each time step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [qq, manips, index1]=path_planning_max_manip_global(direction)
close all;
global robot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%initial and end poses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q0 = [-1.3 -0.7 0.0 1.8 0.3 1.0 0.3]';
T0=directkinematic(robot, q0);
p0 = T0(1:3,4);

%movement in meters from start pose
%this is point 1
deltaV = [0.2 0.85 0.1]';
ds=0.02;

drawrobot3d(robot, q0)
draw_axes(T0, 'Xpiece', 'Ypiece', 'Zpiece', 1.2);

%The inverse kinematics is solved for the following reason
%if using the moore penrose
% this allows the joint positions to be continuous
%can be solved with moore_penrose or with q3_0
pp0 = [p0 p0 + deltaV];
plot3(pp0(1,:),pp0(2,:),pp0(3,:),'k')

mov = 0.0:ds:norm(deltaV);
Ti=T0;
q = q0;
%initial pose and manipulability
qq = [];
pp = [];
pph = [];
j=1;
%select direction of movement
if strcmp(direction, 'forth')
	K = 1:length(mov);
else
    K = length(mov):-1:1;
end
for i=K%length(mov)
    fprintf('Move %d out of %d\n', i, length(mov))
    %update to next point in trajectory
    Ti(1:3,4) = p0 + mov(i)*deltaV;
    %draw_axes(Ti, 'Xpiece', 'Ypiece', 'Zpiece', 1.2);
    %solve inverse kinematics from starting q
    %q = pinv(J)*X
    q = inverse_kinematics_sawyer(robot, Ti, q, 'moore_penrose');
    if j==1
        %q = optimize_manip_global_mcl(robot, q, Ti, 'max_manip');
        q = optimize_manip_global(robot, q, 'max_manip');
        %this is just a refinement on the previous estimate
        q = optimize_manip_local(robot, q, 'max_manip');
    else
        %max manipulability along the null space at that particular pose
        q = optimize_manip_local(robot, q, 'max_manip');
    end
    j=j+1;
    Th=directkinematic(robot, q);
    %is q continuous?
    %drawrobot3d(robot, q)
    %draw_axes(Ti, 'Xpiece', 'Ypiece', 'Zpiece', 1.2);
    %plot3(pp0(1,:),pp0(2,:),pp0(3,:),'k')
    qq = [qq q];
    pp = [pp Ti(1:3,4)];
    pph = [pph Th(1:3,4)];
end

manips = compute_manip(robot, qq);
%total index along trajectory
index1=mean(manips)
index2=prod(manips)

for i=1:5:size(qq,2)
     drawrobot3d(robot, qq(:,i))
     draw_axes(Ti, 'Xpiece', 'Ypiece', 'Zpiece', 1.2);
end
plot3(pp(1,:),pp(2,:),pp(3,:))
plot3(pph(1,:),pph(2,:),pph(3,:),'r')
figure,
plot(manips)
title('manipulability index at each movement')

figure,
plot(qq')
title('joint positions')
legend('q_1', 'q_2','q_3','q_4','q_5','q_6','q_7')

'max manipulability'
max_manip = max(manips)

'min manipulability'
min_manip = min(manips)

'mean manipulability'
mean_manip = mean(manips)

'semi sum manipulability'
semi_manip = (max_manip-min_manip)/2

save test_no_max_manip.mat

function animate_local(robot, q)
global configuration

h=figure(configuration.figure.robot);, hold on,

%get adjusted view
[az,el] = view;
for j=1:size(q, 2),
    clf(h);
    qj=q(:,j);
    view(az,el);

    drawrobot3d(robot, qj);

    %pause to get a nice view
    pause(0.1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check whether orientation has been reached
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function reach = reached_orientation(Qf, Qi)
Q = Qf-Qi;
reach = sqrt(Q(1)^2 + Q(2)^2 + Q(3)^2 + Q(4)^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check whether orientation has been reached
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function reach = reached_position(Pf, Pi)
P = Pf-Pi;
reach = sqrt(P(1)^2 + P(2)^2 + P(3)^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalize vector if possible.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = normalize(v)
d = sqrt(v(1)^2+v(2)^2+v(3)^2);
if d>0
    v = v/d;
end
v=v(:);

repo: ARTE


text:
function S = spm_cfg_eeg_contrast
% configuration file for computing contrast over epochs
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_cfg_eeg_contrast.m 5377 2013-04-02 17:07:57Z vladimir $

D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the EEG mat file.'};

c = cfg_entry;
c.tag = 'c';
c.name = 'Contrast coefficients';
c.strtype = 'r';
c.num = [1 inf];
c.help = {'Enter the contrast vector.'};

label = cfg_entry;
label.tag = 'label';
label.name = 'New condition label';
label.strtype = 's';
label.help = {'Enter the label for the condition derived by applying the contrast.'};

contrast = cfg_branch;
contrast.tag = 'contrast';
contrast.name = 'Contrast';
contrast.val = {c label};

contrasts = cfg_repeat;
contrasts.tag = 'contrasts';
contrasts.name = 'Contrasts';
contrasts.help = {'Each contrast defines a new condition in the output file.'};
contrasts.values  = {contrast};
contrasts.num     = [1 Inf];

weighted = cfg_menu;
weighted.tag = 'weighted';
weighted.name = 'Weight by replications';
weighted.labels = {'Yes', 'No'};
weighted.values = {1 , 0};
weighted.val = {1};
weighted.help = {'Weight the contrast by the numner of replications from which the average was computed.'};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the output dataset. Default prefix is ''w''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'w'};

S = cfg_exbranch;
S.tag = 'contrast';
S.name = 'Contrast over epochs';
S.val = {D contrasts weighted prefix};
S.help = {'Computes contrasts over EEG/MEG epochs.'};
S.prog = @eeg_contrast;
S.vout = @vout_eeg_contrast;
S.modality = {'EEG'};

function out = eeg_contrast(job)
% construct the S struct
S.D = job.D{1};
S.c = cat(1, job.contrast(:).c);
S.label = {job.contrast.label};

S.weighted = job.weighted;
S.prefix   = job.prefix;

out.D = spm_eeg_contrast(S);

out.Dfname = {fullfile(out.D)};

function dep = vout_eeg_contrast(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Contrast of M/EEG epochs';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Contrast Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
close all;
clearvars;
clc;
rng default;
rng(100);
degrees = 1:9;

n = 10;
sigma = 0.25;
[x_train, t_train] = spx.data.synthetic.func.sinusoid('n', n, 'sigma', sigma);

n2 = 100;
[x_test, t_test] = spx.data.synthetic.func.sinusoid('n', n2, 'sigma', 0);

training_errors  = [];
test_errors = [];

for iter=1:numel(degrees)
    degree = degrees(iter);
    features_train = spx.ml.features.polynomial(x_train, degree);
    model = spx.ml.models.linear.LinearRegression;
    model.fit(features_train, t_train);
    y_train = model.predict(features_train);
    features_test = spx.ml.features.polynomial(x_test, degree);
    y_test = model.predict(features_test);
    rmse_train = rms(y_train - t_train);
    rmse_test = rms(y_test - t_test);
    training_errors(end+1) = rmse_train;
    test_errors(end+1) = rmse_test;
end

hold on;
plot(degrees, training_errors, 'o-');
plot(degrees, test_errors, 'o-');
legend({'Training', 'Test'});
xlabel('Degree');
ylabel('RMSE');

author: indigits

repo: sparse-plex


text:
function varargout = process_granger1n( varargin )
% PROCESS_GRANGER1N: Compute the Granger causality between all the pairs of signals, in one file.

% @=============================================================================
% 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, 2012-2020

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Bivariate Granger causality NxN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 661;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 1);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === GRANGER ORDER
    sProcess.options.grangerorder.Comment = 'Maximum Granger model order (default=10):';
    sProcess.options.grangerorder.Type    = 'value';
    sProcess.options.grangerorder.Value   = {10, '', 0};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'granger';
    OPTIONS.RemoveEvoked = sProcess.options.removeevoked.Value;
    OPTIONS.GrangerOrder = sProcess.options.grangerorder.Value{1};
    OPTIONS.pThresh      = 0.05;
    OPTIONS.Freqs        = 0;

    % Compute metric
    OutputFiles = bst_connectivity({sInputA.FileName}, [], OPTIONS);
end

%% ===== TEST FUNCTION =====
function Test() %#ok<DEFNU>
    % Start a new report
    bst_report('Start');
    % Get test datasets
    sFile = process_simulate_ar('Test');
    % Coherence process
    sTmp = bst_process('CallProcess', 'process_granger1n', sFile, [], ...
        'timewindow',   [], ...    % All the time in input
        'grangerorder', 10, ...
        'outputmode',   1);        % Save individual results (one file per input file)
    % Snapshot: spectrum
    bst_process('CallProcess', 'process_snapshot', sTmp, [], ...
        'target',       11, ...  % Connectivity matrix (image)
        'modality',     1, 'orient', 1, 'time', 0, 'contact_time', [-40, 110], 'contact_nimage', 16, ...
        'Comment',      [sFile.Comment, ': ' sTmp.Comment]);
    % Save and display report
    ReportFile = bst_report('Save', sFile);
    bst_report('Open', ReportFile);
end

author: brainstorm-tools


text:
function [ o, c, e ] = polynomial_mul ( m, o1, c1, e1, o2, c2, e2 )

%*****************************************************************************80
%
%% POLYNOMIAL_MUL multiplies two polynomials.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer O1, the "order" of polynomial 1.
%
%    Input, real C1(O1), the coefficients of polynomial 1.
%
%    Input, integer E1(O1), the indices of the exponents of polynomial 1.
%
%    Input, integer O2, the "order" of polynomial 2.
%
%    Input, real C2(O2), the coefficients of polynomial 2.
%
%    Input, integer E2(O2), the indices of the exponents of polynomial 2.
%
%    Output, integer O, the "order" of the polynomial product.
%
%    Output, real C(O), the coefficients of the polynomial product.
%
%    Output, integer E(O), the indices of the exponents of the
%    polynomial product.
%
  o = o1 * o2;
  e = zeros ( o, 1 );
  c = zeros ( o, 1 );

  o = 0;
  for j = 1 : o2
    for i = 1 : o1
      o = o + 1;
      c(o) = c1(i) * c2(j);
      f1(1:m) = mono_unrank_grlex ( m, e1(i) );
      f2(1:m) = mono_unrank_grlex ( m, e2(j) );
      f(1:m) = f1(1:m) + f2(1:m);
      e(o) = mono_rank_grlex ( m, f );
    end
  end

  [ c, e ] = polynomial_sort ( o, c, e );
  [ o, c, e ] = polynomial_compress ( o, c, e );

  return
end


text:
clear all; close all; clc;

% Define path to a urdf file
path_to_urdf = 'ur10e.urdf';

% Generate functions for dynamics based on Lagrange method
% Note that it might take some time
% generate_rb_dynamics(path_to_urdf);
generate_friction_eq();

% Generate regressors for inverse dynamics of the robot, friction and load
% Note that it might take some time
% generate_rb_regressor(path_to_urdf);
% generate_load_regressor(path_to_urdf);

% Run tests
test_rb_inverse_dynamics()
test_base_params()

% Perform QR decompostion in order to get base parameters of the robot
include_motor_dynamics = 1;
[pi_lgr_base, baseQR] = base_params_qr(include_motor_dynamics);

% Estimate drive gains
drive_gains = estimate_drive_gains(baseQR, 'PC-OLS');
% Or use those found in the paper by De Luca
% drive_gains = [14.87; 13.26; 11.13; 10.62; 11.03; 11.47];

% Estimate dynamic parameters
path_to_est_data = 'ur-20_02_10-30sec_12harm.csv';      idxs = [635, 3510];
% path_to_data = 'ur-20_02_12-40sec_12harm.csv';    idxs = [500, 4460];
% path_to_data = 'ur-20_02_05-20sec_8harm.csv';     idxs = [320, 2310];
% path_to_data = 'ur-20_02_12-50sec_12harm.csv';    idxs = [355, 5090];
sol = estimate_dynamic_params(path_to_est_data, idxs, ...
                              drive_gains, baseQR, 'PC-OLS');

% Validate estimated parameters
path_to_val_data = 'ur-20_01_17-ptp_10_points.csv';     idxs = [700, 4200];

rre = validate_dynamic_params(path_to_val_data, idxs, ...
                              drive_gains, baseQR, sol.pi_b, sol.pi_fr)

author: shamilmamedov


text:
function c = tapas_hgf_jget_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Hierarchical Gaussian Filter (HGF) model for the jumping
% Gaussian estimation task (JGET).
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_hgf_jget_plotTraj(est)
%
% where est is the stucture returned by fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mux_0      row vector of initial values of mu_x (in ascending order of levels)
%         est.p_prc.sax_0      row vector of initial values of sigma_x (in ascending order of levels)
%         est.p_prc.mua_0      row vector of initial values of mu_alpha (in ascending order of levels)
%         est.p_prc.saa_0      row vector of initial values of sigma_alpha (in ascending order of levels)
%         est.p_prc.kau        kappa_u
%         est.p_prc.omu        omega_u
%         est.p_prc.kax        row vector of kappa_x (in ascending order of levels)
%         est.p_prc.omx        row vector of omega_x (in ascending order of levels)
%         est.p_prc.kaa        row vector of kappa_alpha (in ascending order of levels)
%         est.p_prc.oma        row vector of omega_alpha (in ascending order of levels)
%
%         est.traj.mux         mux (rows: trials, columns: levels)
%         est.traj.sax         sigma_x (rows: trials, columns: levels)
%         est.traj.muxhat      prediction of mu_x (rows: trials, columns: levels)
%         est.traj.saxhat      precisions of predictions of x (rows: trials, columns: levels)
%         est.traj.wx          weighting factors for x (rows: trials, columns: levels)
%         est.traj.dax         volatility prediction errors in x (rows: trials, columns: levels)
%         est.traj.mua         mu_alpha (rows: trials, columns: levels)
%         est.traj.saa         sigma_alpha (rows: trials, columns: levels)
%         est.traj.muahat      prediction of mu_alpha (rows: trials, columns: levels)
%         est.traj.saahat      precisions of predictions of alpha (rows: trials, columns: levels)
%         est.traj.wa          weighting factors for alpha (rows: trials, columns: levels)
%         est.traj.daa         volatility prediction errors in alpha (rows: trials, columns: levels)
%         est.traj.dau         input prediction error
%
%
%
% Tips:
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013-2014 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 = 'hgf_jget';

% Number of levels (minimum: 2)
c.n_levels = 2;

% Sufficient statistics of Gaussian parameter priors

% PLACEHOLDER VALUES
% It is often convenient to set some priors to values
% derived from the inputs. This can be achieved by
% using placeholder values. The available placeholders
% are:
%
% 99991   Value of the first input
%         Usually a good choice for mux_0mu(1)
% 99992   Variance of the first 20 inputs
%         Usually a good choice for mux_0sa(1)
% 99993   Log-variance of the first 20 inputs
%         Usually a good choice for logsax_0mu(1) and mua_0mu(1)
% 99994   Log-variance of the first 20 inputs minus two
%         Usually a good choice for omxmu(1)

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first level, this is usually best
% kept fixed to 1 (determines origin on x_i-scale).
c.mux_0mu = [99991, 1];
c.mux_0sa = [    0, 0];

c.logsax_0mu = [log(3), log(0.1)];
c.logsax_0sa = [      0,       0];

c.mua_0mu = [log(1), 1];
c.mua_0sa = [   0.1, 0];

c.logsaa_0mu = [log(3), log(0.1)];
c.logsaa_0sa = [      0,       0];

% Kappas
% Format: row vector of length n_levels-1 (except kappa_u: scalar)
% This should be fixed (preferably to 1) if the observation model
% does not use mu_i+1 (kappa then determines the scaling of x_i+1).
c.logkaumu = log(1);
c.logkausa = 0;

c.logkaxmu = [log(1)];
c.logkaxsa = [     0];

c.logkaamu = [log(1)];
c.logkaasa = [     0];

% Omegas
% Format: row vector of length n_levels (except omega_u: scalar)
c.omumu = 0;
c.omusa = 0;

c.omxmu = [  0,  -7];
c.omxsa = [5^2,   1];

c.omamu = [  0,  -7];
c.omasa = [5^2,   1];

% Gather prior settings in vectors
c.priormus = [
    c.mux_0mu,...
    c.logsax_0mu,...
    c.mua_0mu,...
    c.logsaa_0mu,...
    c.logkaumu,...
    c.logkaxmu,...
    c.logkaamu,...
    c.omumu,...
    c.omxmu,...
    c.omamu,...
         ];

c.priorsas = [
    c.mux_0sa,...
    c.logsax_0sa,...
    c.mua_0sa,...
    c.logsaa_0sa,...
    c.logkausa,...
    c.logkaxsa,...
    c.logkaasa,...
    c.omusa,...
    c.omxsa,...
    c.omasa,...
         ];

% Check whether we have the right number of priors
expectedLength = 8*c.n_levels;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength;
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_hgf_jget;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_hgf_jget_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
% test_trifacet_signed_area.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% Triangle with positive area
x = [
    0 0
    1 0
    0 1
    ];

tri = [1 2 3];

a = trifacet_signed_area(tri, x)

% Triangle with negative area
tri = [1 3 2];

a = trifacet_signed_area(tri, x)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus