text:
function [spike] = ft_read_spike(filename, varargin)

% FT_READ_SPIKE reads spike timestamps and waveforms from various data formats.
%
% Use as
%  [spike] = ft_read_spike(filename, ...)
%
% Additional options should be specified in key-value pairs and can be
%   'spikeformat' = string, described the fileformat (default is automatic)
%
% The following file formats are supported
%   'blackrock_nev'
%   'mclust_t'
%   'neuralynx_ncs'
%   'neuralynx_nse'
%   'neuralynx_nst'
%   'neuralynx_ntt'
%   'neuralynx_nts'
%   'neuroshare'
%   'neurosim_spikes'
%   'nwb'
%   'plexon_ddt'
%   'plexon_nex'
%   'plexon_nex5'
%   'plexon_plx'
%   'wave_clus'
%
% The output spike structure usually contains
%   spike.label     = 1xNchans cell-array, with channel labels
%   spike.waveform  = 1xNchans cell-array, each element contains a matrix (Nleads x Nsamples X Nspikes)
%   spike.waveformdimord = '{chan}_lead_time_spike'
%   spike.timestamp = 1xNchans cell-array, each element contains a vector (1 X Nspikes)
%   spike.unit      = 1xNchans cell-array, each element contains a vector (1 X Nspikes)
% and is described in more detail in FT_DATATYPE_SPIKE
%
% See also FT_DATATYPE_SPIKE, FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT

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

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

if ~exist(filename,'file')
  ft_error('File or directory does not exist')
end

% get the options
spikeformat = ft_getopt(varargin, 'spikeformat');

% optionally get the data from the URL and make a temporary local copy
filename = fetch_url(filename);

if isempty(spikeformat)
  % only do the autodetection if the format was not specified
  spikeformat = ft_filetype(filename);
end

switch spikeformat
  case {'neurosim_spikes' 'neurosim_ds'}
    spike = read_neurosim_spikes(filename);

  case {'neuralynx_ncs' 'plexon_ddt'}
    % these files only contain continuous data
    ft_error('file does not contain spike timestamps or waveforms');

  case 'matlab'
    % plain MATLAB file with a single variable in it
    load(filename, 'spike');

  case 'mclust_t'
    fp = fopen(filename, 'rb', 'ieee-le');
    H = ReadHeader(fp);
    fclose(fp);
    % read only from one file
    S = read_mclust_t({filename});
    spike.hdr = H(:);
    [p, f, x] = fileparts(filename);
    spike.label     = {f};  % use the filename as label for the spike channel
    spike.timestamp = S;
    spike.hdr       = H;

  case 'wave_clus'
    load(filename, 'cluster_class', 'spikes', 'par'); % load the mat file
    clusters = sort(unique(cluster_class(:,1))); % detected clusters
    clusters(clusters==0) = []; % remove rejected cluster (indexed by zeros)
    nclust = numel(clusters);
    [p, f, x] = fileparts(filename);
    t = tokenize(f, '_'); % extract channel name
    spike.label     = cell(1,nclust);
    spike.unit      = cell(1,nclust);
    spike.waveform  = cell(1,nclust);
    spike.timestamp = cell(1,nclust);
    spike.hdr       = par;
    for cl = 1:nclust
      unit_idx                  = cluster_class(:,1)==cl;
      spike.label{cl}           = [t{2} '-' num2str(cl)];
      spike.timestamp{cl}       = cluster_class(unit_idx,2)';
      spike.waveform{cl}(1,:,:) = spikes(unit_idx,:)';
      spike.unit{cl}            = cluster_class(unit_idx,1)';
    end
    fprintf('note that wave_clus timestamps are typically expressed in millisec and not in samples\n')

  case 'neuralynx_nse'
    % single channel file, read all records
    nse = read_neuralynx_nse(filename);
    if isfield(nse.hdr, 'NLX_Base_Class_Name')
      spike.label   = {nse.hdr.NLX_Base_Class_Name};
    else
      spike.label   = {nse.hdr.AcqEntName};
    end
    spike.timestamp = {nse.TimeStamp};
    spike.waveform  = {nse.dat};
    spike.unit      = {nse.CellNumber};
    spike.hdr       = nse.hdr;

  case 'neuralynx_nst'
    % single channel stereotrode file, read all records
    nst = read_neuralynx_nst(filename, 1, inf);
    if isfield(nst.hdr, 'NLX_Base_Class_Name')
      spike.label   = {nst.hdr.NLX_Base_Class_Name};
    else
      spike.label   = {nst.hdr.AcqEntName};
    end
    spike.timestamp = {nst.TimeStamp};
    spike.waveform  = {nst.dat};
    spike.unit      = {nst.CellNumber};
    spike.hdr       = nst.hdr;

  case 'neuralynx_ntt'
    % single channel stereotrode file, read all records
    ntt = read_neuralynx_ntt(filename);
    if isfield(ntt.hdr, 'NLX_Base_Class_Name')
      spike.label   = {ntt.hdr.NLX_Base_Class_Name};
    else
      spike.label   = {ntt.hdr.AcqEntName};
    end
    spike.timestamp = {ntt.TimeStamp};
    spike.waveform  = {ntt.dat};
    spike.unit      = {ntt.CellNumber};
    spike.hdr       = ntt.hdr;

  case 'neuralynx_nts'
    % single channel file, read all records
    nts = read_neuralynx_nts(filename);
    if isfield(nte.hdr, 'NLX_Base_Class_Name')
      spike.label   = {nts.hdr.NLX_Base_Class_Name};
    else
      spike.label   = {nts.hdr.AcqEntName};
    end
    spike.timestamp = {nts.TimeStamp(:)'};
    spike.waveform  = {zeros(0,length(nts.TimeStamp))};  % does not contain waveforms
    spike.unit      = {zeros(0,length(nts.TimeStamp))};  % does not contain units
    spike.hdr       = nts.hdr;

  case 'plexon_nex'
    % a single file can contain multiple channels of different types
    hdr  = read_plexon_nex(filename);
    typ  = [hdr.VarHeader.Type];
    chan = 0;

    spike.label     = {};
    spike.timestamp = {};
    spike.waveform  = {};
    spike.unit      = {};

    for i=1:length(typ)
      if typ(i)==0
        % neurons, only timestamps
        nex = read_plexon_nex(filename, 'channel', i);
        nspike = length(nex.ts);
        chan = chan + 1;
        spike.label{chan}     = deblank(hdr.VarHeader(i).Name);
        spike.waveform{chan}  = zeros(0, nspike);
        spike.unit{chan}      = nan(1,nspike);
        spike.timestamp{chan} = nex.ts;
      elseif typ(i)==3
        % neurons, timestamps and waveforms
        nex = read_plexon_nex(filename, 'channel', i);
        chan = chan + 1;
        nspike = length(nex.ts);
        spike.label{chan}     = deblank(hdr.VarHeader(i).Name);
        spike.waveform{chan}  = permute(nex.dat,[3 1 2]);
        spike.unit{chan}      = nan(1,nspike);
        spike.timestamp{chan} = nex.ts;
      end
    end
    spike.hdr = hdr;

  case 'plexon_nex5'
    % a single file can contain multiple channels of different types
    hdr  = read_nex5(filename);
    typ  = [hdr.VarHeader.Type];
    chan = 0;

    spike.label     = {};
    spike.timestamp = {};
    spike.waveform  = {};
    spike.unit      = {};

    for i=1:length(typ)
      if typ(i)==0
        % neurons, only timestamps
        nex = read_nex5(filename, 'channel', i);
        nspike = length(nex.ts);
        chan = chan + 1;
        spike.label{chan}     = deblank(hdr.VarHeader(i).Name);
        spike.waveform{chan}  = zeros(0, nspike);
        spike.unit{chan}      = nan(1,nspike);
        % spike.timestamp{chan} are the raw timestamps as recorded by the hardware system
        spike.timestamp{chan} = nex.ts;
      elseif typ(i)==3
        % waveform variables: timestamps and waveforms
        nex = read_nex5(filename, 'channel', i);
        chan = chan + 1;
        nspike = length(nex.ts);
        spike.label{chan}     = deblank(hdr.VarHeader(i).Name);
        spike.waveform{chan}  = permute(nex.dat,[3 1 2]);
        spike.unit{chan}      = nan(1,nspike);
        % spike.timestamp{chan} are the raw timestamps as recorded by the hardware system
        spike.timestamp{chan} = nex.ts;
      end
    end
    spike.hdr = hdr;

  case 'plexon_plx'
    % read the header information
    hdr   = read_plexon_plx(filename);
    nchan = length(hdr.ChannelHeader);
    typ   = [hdr.DataBlockHeader.Type];
    unit  = [hdr.DataBlockHeader.Unit];
    chan  = [hdr.DataBlockHeader.Channel];

    for i=1:nchan
      % select the data blocks that contain spike waveforms and that belong to this channel
      sel = (typ==1 & chan==hdr.ChannelHeader(i).Channel);

      if any(sel)
        % get the timestamps that correspond with this spike channel
        tsl = [hdr.DataBlockHeader(sel).TimeStamp];
        tsh = [hdr.DataBlockHeader(sel).UpperByteOf5ByteTimestamp];
        % convert the 16 bit high timestamp into a 32 bit integer
        ts = timestamp_plexon(tsl, tsh);
        spike.timestamp{i} = ts;
        spike.unit{i}      = unit(sel);
      else
        % this spike channel is empty
        spike.timestamp{i} = [];
        spike.unit{i}      = [];
      end
    end
    for i=1:nchan
      spike.label{i}    = deblank(hdr.ChannelHeader(i).Name);
      spike.waveform{i} = permute(read_plexon_plx(filename, 'ChannelIndex', i, 'header', hdr),[3 1 2]);
    end
    spike.hdr = hdr;

  case 'neuroshare' % NOTE: still under development
    % check that the required neuroshare toolbox is available
    ft_hastoolbox('neuroshare', 1);

    tmp = read_neuroshare(filename, 'readspike', 'yes');
    spike.label = {tmp.hdr.entityinfo(tmp.list.segment).EntityLabel};
    for i=1:length(spike.label)
      spike.waveform{i}  = tmp.spikew.data(:,:,i);
      spike.timestamp{i} = tmp.spikew.timestamp(:,i)';
      spike.unit{i}      = tmp.spikew.unitID(:,i)';
    end

  case 'neuroscope'
    % the information about the spikes is represented in:
    % x.clu.y or x.res.y (containing the timing +cluster info)
    % x.spk.y (containing the waveform info)
    % x.fet.y (containing features: do we need this?)

    if isfolder(filename)
      tmp = dir(filename);
      filenames = {tmp.name}';
    end

    % read the header
    filename_hdr = filenames{~cellfun('isempty',strfind(filenames,'.xml'))};
    hdr          = ft_read_header(fullfile(filename,filename_hdr), 'headerformat', 'neuroscope_xml');
    spikegroups  = hdr.orig.spikeGroups;
    fsample      = hdr.orig.rates.wideband;

    filename_clu = filenames(~cellfun('isempty',strfind(filenames,'.clu')));
    filename_spk = filenames(~cellfun('isempty',strfind(filenames,'.spk')));

    % FIXME should we do a sanity check on whether the clu and spk actually
    % belong together?

    c = cell(numel(filename_clu),1);
    w = cell(numel(filename_spk),1);
    for k = 1:numel(filename_clu)
      c{k} = LoadSpikeTimes(fullfile(filename,filename_clu{k}), fsample);
    end
    for k = 1:numel(filename_spk)
      w{k} = LoadSpikeWaveforms(fullfile(filename,filename_spk{k}),numel(spikegroups.groups{k}),spikegroups.nSamples(k));
    end

    spike = [];
    spike.label = cell(hdr.orig.spikeGroups.nGroups,1);
    spike.hdr   = hdr;
    spike.unit      = cell(1,numel(spike.label));
    spike.waveform  = cell(1,numel(spike.label));
    spike.timestamp = cell(1,numel(spike.label));

    for k = 1:numel(spike.label)
      sel = find(c{k}(:,3)>1); % values >1 corresponds to individual units, 0 = noise, 1 = MUA

      % the times are defined in s, convert to original time stamps
      timestamps = c{k}(sel,1) * hdr.orig.rates.wideband;
      if any(abs(timestamps-round(timestamps))>1e-5)
        ft_error('there seems to be a mismatch between the spike times and the expected integer-valued timestamps');
      end

      spike.timestamp{k} = round(timestamps(:))';
      spike.waveform{k}  = permute(w{k}(sel,:,:), [2 3 1]);
      spike.unit{k}      = c{k}(sel,3)';
      spike.label{k}     = sprintf('spikegroup%03d',k);
    end

  case 'nwb'
    ft_hastoolbox('MatNWB', 1);
    spike = read_nwb_spike(filename);

  case {'blackrock_nev'}
    % use the NPMK toolbox for the file reading
    ft_hastoolbox('NPMK', 1);

    % ensure that the filename contains a full path specification,
    % otherwise the low-level function fails
    [p, f, x] = fileparts(filename);
    if ~isempty(p)
      % this is OK
    elseif isempty(p)
      filename = which(filename);
    end

    % 'nosave' prevents the automatic conversion of the .nev file as a .mat file
    nev = openNEV(filename, 'nosave');

    nchan = length(nev.ElectrodesInfo);
    for i=1:nchan
      spike.label{i} = deblank(nev.ElectrodesInfo(i).ElectrodeLabel(:)');
      if isfield(nev.Data, 'Spikes')
        % select the spikes that were detected on this electrode
        sel = nev.Data.Spikes.Electrode;
        spike.timestamp{i} = nev.Data.Spikes.TimeStamp(sel);
        spike.unit{i}      = nev.Data.Spikes.Unit(sel);
        spike.waveform{i}  = nev.Data.Spikes.Waveform(sel);
      end
    end

  otherwise
    ft_error(['unsupported data format (' spikeformat ')']);
end

% add the waveform
if isfield(spike,'waveform')
  spike.dimord = '{chan}_lead_time_spike';
end


text:
function [data_train, labels_train, data_test, labels_test, raw_test, PC, means_norm, stds_norm, vid_ids_test, success_test] = ...
    Prepare_HOG_AU_data_generic_dynamic(train_users, test_users, au_train, rest_aus, root, features_dir)

%% This should be a separate function?

input_train_label_files = cell(numel(train_users),1);
input_test_label_files = cell(numel(test_users),1);

% This is for loading the labels
for i=1:numel(train_users)
    input_train_label_files{i} = [root, '/ActionUnit_Labels/', train_users{i}, '/', train_users{i}];
end

% This is for loading the labels
for i=1:numel(test_users)
    input_test_label_files{i} = [root, '/ActionUnit_Labels/', test_users{i}, '/', test_users{i}];
end

% First extracting the labels
[train_geom_data] = Read_geom_files_dynamic(train_users, features_dir);
[test_geom_data] = Read_geom_files_dynamic(test_users, features_dir);

% Reading in the HOG data
[train_data, tracked_inds_hog, vid_ids_train] = Read_HOG_files_dynamic(train_users, features_dir);
[test_data, success_test, vid_ids_test] = Read_HOG_files_dynamic(test_users, features_dir);

train_data = cat(2, train_data, train_geom_data);
raw_test = cat(2, test_data, test_geom_data);

% Extracting the labels
labels_train = extract_au_labels(input_train_label_files, au_train);
labels_test = extract_au_labels(input_test_label_files, au_train);

labels_other = zeros(size(labels_train,1), numel(rest_aus));

% This is used to pick up activity of other AUs for a more 'interesting'
% data split and not only neutral expressions for negative samples
if(numel(input_train_label_files) > 0)
    for i=1:numel(rest_aus)
        labels_other(:,i) = extract_au_labels(input_train_label_files, rest_aus(i));
    end

    % can now extract the needed training labels (do not rebalance validation
    % data)

    % make sure the same number of positive and negative samples is taken
    reduced_inds = false(size(labels_train,1),1);
    reduced_inds(labels_train > 0) = true;

    % make sure the same number of positive and negative samples is taken
    pos_count = sum(labels_train > 0);
    neg_count = sum(labels_train == 0);

    % pos_count = pos_count * 8;

    num_other = floor(pos_count / (size(labels_other, 2)));

    inds_all = 1:size(labels_train,1);

    for i=1:size(labels_other, 2)+1

        if(i > size(labels_other, 2))
            % fill the rest with a proportion of neutral
            inds_other = inds_all(sum(labels_other,2)==0 & ~labels_train);
            num_other_i = min(numel(inds_other), pos_count - sum(labels_train(reduced_inds,:)==0));
        else
            % take a proportion of each other AU
            inds_other = inds_all(labels_other(:, i) & ~labels_train);
            num_other_i = min(numel(inds_other), num_other);
        end
        inds_other_to_keep = inds_other(round(linspace(1, numel(inds_other), num_other_i)));
        reduced_inds(inds_other_to_keep) = true;

    end

    % Remove invalid ids based on CLM failing or AU not being labelled
    reduced_inds(~tracked_inds_hog) = false;

    labels_train = labels_train(reduced_inds);
    train_data = train_data(reduced_inds,:);

end

geom_size = max(size(train_geom_data,2), size(test_geom_data,2));

% Peforming zone specific masking
if(au_train < 8 || au_train == 43 || au_train == 45) % upper face AUs ignore bottom face
    % normalise the data
    pca_file = '../../pca_generation/generic_face_upper.mat';
    load(pca_file);
elseif(au_train > 9) % lower face AUs ignore upper face and the sides
    % normalise the data
    pca_file = '../../pca_generation/generic_face_lower.mat';
    load(pca_file);
elseif(au_train == 9) % Central face model
    % normalise the data
    pca_file = '../../pca_generation/generic_face_rigid.mat';
    load(pca_file);
end

PC_n = zeros(size(PC)+geom_size);
PC_n(1:size(PC,1), 1:size(PC,2)) = PC;
PC_n(size(PC,1)+1:end, size(PC,2)+1:end) = eye(geom_size);
PC = PC_n;

means_norm = cat(2, means_norm, zeros(1, geom_size));
stds_norm = cat(2, stds_norm, ones(1, geom_size));

data_test = bsxfun(@times, bsxfun(@plus, raw_test, -means_norm), 1./stds_norm);
data_test = data_test * PC;

if(numel(train_data > 0))
    data_train = bsxfun(@times, bsxfun(@plus, train_data, -means_norm), 1./stds_norm);
    data_train = data_train * PC;
else
   data_train = [];
end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function out = isinf(f)
%ISINF   Test if a TRIGTECH is unbounded.
%   ISINF(F) returns TRUE if F has any infinite values and FALSE otherwise.
%
% See also ISFINITE, ISNAN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check if any values are infinite:
out = any(isinf(f.values(:)));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

for c = 1:p
	temp = X .^ c;
	X_poly(:, c) = temp;

% =========================================================================

end

author: sfvsfv

repo: Mathematical-modeling


text:
function check_for_cd_arcs(onodes, cnodes, dag)
% CHECK_FOR_CD_ARCS Raise an error if there are any C->D links where the C node is hidden.
% check_for_cd_arcs(onodes, cnodes, dag)
%
% We cannot convert the logistic/softmax function (C->D CPD) to a Gaussian potential
% unless we use the variational approximation discussed in
% "A variational approximation for Bayesian networks with discrete and continuous latent
% variables", K. Murphy, UAI 1999.

n = length(dag);
hnodes = mysetdiff(1:n, onodes);
chid = myintersect(cnodes, hnodes);
dnodes = mysetdiff(1:n, cnodes);
for i=chid(:)'
  dcs = myintersect(children(dag, i), dnodes);
  if ~isempty(dcs)
    error(['hidden cts node ' num2str(i) ' has a discrete child']);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [hd_record,hd_estimates]=panel3hd(beta_gibbs,D_record,strshocks_record,It,Bu,Yi,Xi,N,n,m,p,k,T,HDband)

% initiate the cells storing the results
hd_record={};
hd_estimates={};
% because historical decomposition has to be computed for ach unit, loop over units
for ii=1:N
% run the Gibbs sampler for historical decomposition
hd_record(:,:,ii)=bear.hdecomp(beta_gibbs(:,:,ii),D_record(:,:,ii),strshocks_record(:,:,ii),It,Bu,Yi(:,:,ii),Xi(:,:,ii),n,m,p,k,T);
% then obtain point esimates and credibility intervals
hd_estimates(:,:,ii)=bear.hdestimates(hd_record(:,:,ii),n,T,HDband);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function results = vl_test_pegasos(varargin)
% VL_TEST_KDTREE
vl_test_init ;

function s = setup()
randn('state',0) ;

s.biasMultiplier = 10 ;
s.lambda = 0.01 ;

Np = 10 ;
Nn = 10 ;
Xp = diag([1 3])*randn(2, Np) ;
Xn = diag([1 3])*randn(2, Nn) ;
Xp(1,:) = Xp(1,:) + 2 + 1 ;
Xn(1,:) = Xn(1,:) - 2 + 1 ;

s.X = [Xp Xn] ;
s.y = [ones(1,Np) -ones(1,Nn)] ;
%s.w = exact_solver(s.X, s.y, s.lambda, s.biasMultiplier)
s.w = [1.181106685845652 ;
       0.098478251033487 ;
       -0.154057992404545 ] ;

function test_problem_1(s)
for conv = {@single,@double}
  vl_twister('state',0) ;
  conv = conv{1} ;
  [w b info] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'MaxIterations', 100000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'BiasLearningRate', .1) ;

  % test input
  vl_assert_equal(info.biasMultiplier,s.biasMultiplier);
  vl_assert_almost_equal(info.biasLearningRate,.1,1e-3);
  vl_assert_almost_equal(conv([w; b]), conv(s.w), 0.1) ;
end

function test_continue_training(s)
for conv = {@single,@double}
  conv = conv{1} ;

  vl_twister('state',0) ;
  [w b] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'MaxIterations', 3000, ...
                 'BiasMultiplier', s.biasMultiplier) ;

  vl_twister('state',0) ;
  [w1 b1] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'StartingIteration', 1, ...
                 'MaxIterations', 1500, ...
                  'BiasMultiplier', s.biasMultiplier) ;
  [w2 b2] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                  'StartingIteration', 1501, ...
                  'StartingModel', w1, ...
                  'StartingBias', b1, ...
                  'MaxIterations', 3000, ...
                  'BiasMultiplier', s.biasMultiplier) ;
  vl_assert_almost_equal([w; b],[w2; b2],1e-7) ;
end

function test_continue_training_with_perm(s)
perm = uint32(randperm(size(s.X,2))) ;
for conv = {@single,@double}
  conv = conv{1} ;

  vl_twister('state',0) ;
  [w b] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'MaxIterations', 3000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'Permutation', perm) ;

  vl_twister('state',0) ;
  [w1 b1] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'StartingIteration', 1, ...
                 'MaxIterations', 1500, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                  'Permutation', perm) ;
  [w2 b2] = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                  'StartingIteration', 1501, ...
                  'StartingModel', w1, ...
                  'StartingBias', b1, ...
                  'MaxIterations', 3000, ...
                  'BiasMultiplier', s.biasMultiplier, ...
                  'Permutation', perm) ;

  vl_assert_almost_equal([w; b],[w2; b2],1e-7) ;
end

function test_homkermap(s)
for conv = {@single,@double}
  vl_twister('state',0) ;
  conv = conv{1} ;
  sxe = vl_homkermap(conv(s.X), 1, 'kchi2', 'gamma', .5) ;
  [we be] = vl_pegasos(sxe, int8(s.y), s.lambda, ...
                 'MaxIterations', 100000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'BiasLearningRate', .1) ;
  vl_twister('state',0) ;
  [w b] = vl_pegasos(s.X, int8(s.y), s.lambda, ...
                 'MaxIterations', 100000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'BiasLearningRate', .1,...
                 'homkermap',1,...
                 'gamma',.5,...
                 'kchi2') ;

  vl_assert_almost_equal([w; b],[we; be], 1e-7) ;
end

function test_diagnostic(s)
for conv = {@single,@double}
  vl_twister('state',0) ;
  conv = conv{1} ;

  x = 0;
  dhandle = @(x,stat) (assert(stat.elapsedTime == 0 || stat.elapsedTime ~= 0)) ;

  [w b] = vl_pegasos(s.X, int8(s.y), s.lambda, ...
                        'MaxIterations', 100000, ...
                        'BiasMultiplier', s.biasMultiplier, ...
                        'BiasLearningRate', .1) ;
  vl_twister('state',0) ;
  [wd bd] = vl_pegasos(s.X, int8(s.y), s.lambda, ...
                          'MaxIterations', 100000, ...
                          'BiasMultiplier', s.biasMultiplier, ...
                          'BiasLearningRate', .1,...
                          'DiagnosticFunction',dhandle,...
                          'DiagnosticCallRef',x) ;

  vl_assert_almost_equal([w; b], [wd; bd], 1e-7) ;
end

function test_epsilon(s)
for conv = {@single,@double}
  vl_twister('state',0) ;
  conv = conv{1} ;

  [w b info] = vl_pegasos(s.X, int8(s.y), s.lambda, ...
                        'MaxIterations', 1000000, ...
                        'BiasMultiplier', s.biasMultiplier, ...
                        'BiasLearningRate', .1) ;
  vl_twister('state',0) ;
  [we be infoe] = vl_pegasos(s.X, int8(s.y), s.lambda, ...
                          'MaxIterations', 1000000, ...
                          'Epsilon',1e-7,...
                          'BiasMultiplier', s.biasMultiplier, ...
                          'BiasLearningRate', .1) ;

  vl_assert_almost_equal([w; b], [we; be], 1e-2) ;
  assert(info.iterations > infoe.iterations);
end

function w = exact_solver(X, y, lambda, biasMultiplier)
N = size(X,2) ;
model = svmtrain(y', [(1:N)' X'*X], sprintf(' -c %f -t 4 ', 1/(lambda*N))) ;
w = X(:,model.SVs) * model.sv_coef ;
w(3) = - model.rho / biasMultiplier ;
format long ;
disp('model w:')
disp(w)

author: Cloud-CV

repo: object-proposals


text:
function varargout = quad2d(varargin)
%QUAD2D  Complete definite integral of SPHEREFUN.
%   I = QUAD2D( F ), returns the definite integral of a SPHEREFUN integrated
%   over its domain of definition.
%
%   I = QUAD2D(F, a, b, c, d), returns the definite integral of a SPHEREFUN.
%   Integrated over the domangle [a b] x [c d].
%
% See also SPHEREFUN/INTEGRAL2, SPHEREFUN/SUM2, SPHEREFUN/INTEGRAL.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = quad2d@separableApprox(varargin{:});

end


text:
function   r = assortativity_bin(CIJ,flag)
% ASSORTATIVITY_BIN      Assortativity coefficient
%
%   r = assortativity(CIJ,flag);
%
%   The assortativity coefficient is a correlation coefficient between the
%   degrees of all nodes on two opposite ends of a link. A positive
%   assortativity coefficient indicates that nodes tend to link to other
%   nodes with the same or similar degree.
%
%   Inputs:     CIJ,    binary directed/undirected connection matrix
%               flag,   0, undirected graph: degree/degree correlation
%                       1, directed graph: out-degree/in-degree correlation
%                       2, directed graph: in-degree/out-degree correlation
%                       3, directed graph: out-degree/out-degree correlation
%                       4, directed graph: in-degree/in-degree correlation
%
%   Outputs:    r,      assortativity coefficient
%
%   Notes: The function accepts weighted networks, but all connection
%   weights are ignored. The main diagonal should be empty. For flag 1
%   the function computes the directed assortativity described in Rubinov
%   and Sporns (2010) NeuroImage.
%
%   Reference:  Newman (2002) Phys Rev Lett 89:208701
%               Foster et al. (2010) PNAS 107:10815�10820
%
%   Olaf Sporns, Indiana University, 2007/2008
%   Vassilis Tsiaras, University of Crete, 2009
%   Murray Shanahan, Imperial College London, 2012
%   Mika Rubinov, University of Cambridge, 2012

if (flag==0)                        % undirected version
    deg = degrees_und(CIJ);
    [i,j] = find(triu(CIJ,1)>0);
    K = length(i);
    degi = deg(i);
    degj = deg(j);

else                                % directed versions
    [id,od] = degrees_dir(CIJ);
    [i,j] = find(CIJ>0);
    K = length(i);

    switch flag
        case 1
            degi = od(i);
            degj = id(j);
        case 2
            degi = id(i);
            degj = od(j);
        case 3
            degi = od(i);
            degj = od(j);
        case 4
            degi = id(i);
            degj = id(j);
    end
end

% compute assortativity
r = ( sum(degi.*degj)/K - (sum(0.5*(degi+degj))/K)^2 ) / ...
    ( sum(0.5*(degi.^2+degj.^2))/K - (sum(0.5*(degi+degj))/K)^2 );

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function J = jacobian( f, g )
%JACOBIAN   Jacobian determinant of two CHEBFUN2.
%   J = JACOBIAN(F,G) returns the Jacobian determinant of the Jacobian matrix.
%
%   Note we return the determinant of the Jacobian matrix and not the Jacobian
%   matrix itself.
%
% See also CHEBFUN2V/JACOBIAN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( f ) || isempty( g ) )
    J = [];
    return
end

% Call CHEBFUN2V/JACOBIAN():
J = jacobian( chebfun2v( { f, g } ) );

end


text:
function F = uniaxial(d,rate)
% deformationGradientTensor representing uniaxial compression / tension
%
% Syntax
%
%   F = velocityGradientTensor.uniaxial(d,rate)
%
% Input
%  d - @vector3d expansion direction
%  rate - deformation rate
%
% Output
%  F - @deformationGradientTensor
%

if nargin < 2, rate = 2; end

rot = orientation.map(xvector,d);

F = deformationGradientTensor(diag([rate,1./sqrt(rate),1./sqrt(rate)]));

F = rotate(F,inv(rot));

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function adcDwellTime = getAdcDwellTime(MRSIStruct)
    adcDwellTime = MRSIStruct.adcDwellTime;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function movePTPEllipse_YZ_1(iiwa,c,ratio,theta,velocity,accel,TefTool)
%% This funciton is used to
% Move end-effector of the robot on an ellipse the plane of the ellipse is
% parallel to the XY plane of the robot base.

%% Syntax:
% movePTPEllipse_YZ(iiwa,c,ratio,theta,velocity,accel,TefTool)

%% Arreguments:
% iiwa: KST obejct.
% c: 2x1 vector, is a vector defining the displacment of the the center
% point of the ellipse in the YZ plane in relation to the starting point of the ellipse,
% position units in (mm).
% ratio: the radious ratio (a/b) of the ellipse.
% theta: is the angle of the ellipe part , for drawing a complete
% ellipse this angle is equal to 2*pi.
% TefTool: 4X4 Transform matrix of the end-efector in the reference frame
% of the flange of the KUKA iiwa, position units in (mm).
% accel: is the acceleration (mm/sec2).
% velocity: is the motion velocity (mm/sec).

% Copy right: Mohammad SAFEEA
% 30-April-2017

%% Convert inputs to a column vectors
c=colVec(c);
ratio=colVec(ratio);
theta=colVec(theta);
velocity=colVec(velocity);
accel=colVec(accel);

%% Check input variables
if(size(c,1)~=2)
fprintf('Error: the vector (c) shall be a 2x1 vector \n');
return;
end

if(size(ratio,1)~=1)
fprintf('Error: the ratio shall be a scalar  \n');
return;
end

if(size(theta,1)~=1)
fprintf('Error: the variable (theta) shall be a scalar  \n');
return;
end

if(size(TefTool,1)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(TefTool,2)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(velocity,1)~=1)
fprintf('Error: the velocity on the path shall be a scalar  \n');
return;
end

if(size(accel,1)~=1)
fprintf('Error: the acceleration on the path shall be a scalar  \n');
return;
end

dir=[0;1;0];
c=[0;
    c(1);
    c(2)];

% perform the elliptical motion
movePTPEllipse(iiwa,c,dir,ratio,theta,velocity,accel,TefTool);
end

function y=colVec(x)
if size(x,2)==1
    y=x;
else
    y=x';
end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [annualRet,annualCov,annualStd] = ComputeHistoricalStats(prices)

returns = tick2ret(prices,[],'continuous');
numReturns = size(returns,1);

[annualRet,annualCov] = geom2arith(mean(returns),cov(returns),numReturns);

annualRet = annualRet';
annualStd = sqrt(diag(annualCov));

author: Sable

repo: mcbench-benchmarks


text:
function rfs = getVarReceptiveFields(obj, var)
%GETVARRECEPTIVEFIELDS Get the receptive field of a variable
%   RFS = GETVARRECEPTIVEFIELDS(OBJ, VAR) gets the receptivie fields RFS of
%   all the variables of the DagNN OBJ into variable VAR. VAR is a variable
%   name or index.
%
%   RFS has one entry for each variable in the DagNN following the same
%   format as has DAGNN.GETRECEPTIVEFIELDS(). For example, RFS(i) is the
%   receptive field of the i-th variable in the DagNN into variable VAR. If
%   the i-th variable is not a descendent of VAR in the DAG, then there is
%   no receptive field, indicated by `rfs(i).size == []`. If the receptive
%   field cannot be computed (e.g. because it depends on the values of
%   variables and not just on the network topology, or if it cannot be
%   expressed as a sliding window), then `rfs(i).size = [NaN NaN]`.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi. All rights reserved.
%
% This file is part of the VLFeat library and is made available under the
% terms of the BSD license (see the COPYING file).

if ~isnumeric(var)
  var_n = obj.getVarIndex(var) ;
  if isnan(var_n)
    error('Variable %s not found.', var);
  end
  var = var_n;
end
nv = numel(obj.vars) ;
nw = numel(var) ;
rfs = struct('size', cell(nw, nv), 'stride', cell(nw, nv), 'offset', cell(nw,nv)) ;

for w = 1:numel(var)
  rfs(w,var(w)).size = [1 1] ;
  rfs(w,var(w)).stride = [1 1] ;
  rfs(w,var(w)).offset = [1 1] ;
end

for l = obj.executionOrder
  % visit all blocks and get their receptive fields
  in = obj.layers(l).inputIndexes ;
  out = obj.layers(l).outputIndexes ;
  blockRfs = obj.layers(l).block.getReceptiveFields() ;

  for w = 1:numel(var)
    % find the receptive fields in each of the inputs of the block
    for i = 1:numel(in)
      for j = 1:numel(out)
        rf = composeReceptiveFields(rfs(w, in(i)), blockRfs(i,j)) ;
        rfs(w, out(j)) = resolveReceptiveFields([rfs(w, out(j)), rf]) ;
      end
    end
  end
end
end

% -------------------------------------------------------------------------
function rf = composeReceptiveFields(rf1, rf2)
% -------------------------------------------------------------------------
if isempty(rf1.size) || isempty(rf2.size)
  rf.size = [] ;
  rf.stride = [] ;
  rf.offset = [] ;
  return ;
end

rf.size = rf1.stride .* (rf2.size - 1) + rf1.size ;
rf.stride = rf1.stride .* rf2.stride ;
rf.offset = rf1.stride .* (rf2.offset - 1) + rf1.offset ;
end

% -------------------------------------------------------------------------
function rf = resolveReceptiveFields(rfs)
% -------------------------------------------------------------------------

rf.size = [] ;
rf.stride = [] ;
rf.offset = [] ;

for i = 1:numel(rfs)
  if isempty(rfs(i).size), continue ; end
  if isnan(rfs(i).size)
    rf.size = [NaN NaN] ;
    rf.stride = [NaN NaN] ;
    rf.offset = [NaN NaN] ;
    break ;
  end
  if isempty(rf.size)
    rf = rfs(i) ;
  else
    if ~isequal(rf.stride,rfs(i).stride)
      % incompatible geometry; this cannot be represented by a sliding
      % window RF field and may denotes an error in the network structure
      rf.size = [NaN NaN] ;
      rf.stride = [NaN NaN] ;
      rf.offset = [NaN NaN] ;
      break;
    else
      % the two RFs have the same stride, so they can be recombined
      % the new RF is just large enough to contain both of them
      a = rf.offset - (rf.size-1)/2 ;
      b = rf.offset + (rf.size-1)/2 ;
      c = rfs(i).offset - (rfs(i).size-1)/2 ;
      d = rfs(i).offset + (rfs(i).size-1)/2 ;
      e = min(a,c) ;
      f = max(b,d) ;
      rf.offset = (e+f)/2 ;
      rf.size = f-e+1 ;
    end
  end
end
end

repo: LapSRN


text:
function options = dsdpset(varargin)
%DSDPSET  Create or alter the options for Optimization with DSDP
%
% options = dsdpset('param1',value1,'param2',value2,...) creates an DSDP
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the DSDP
% default.
%
% options = dsdpset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = dsdpset() creates an options structure with all fields set to
% dsdpset defaults.
%
% dsdpset() prints a list of all possible fields and their function.
%
% See supplied DSDP Documentation for further details of these options.

%   Copyright (C) 2013 Jonathan Currie (IPL)

% Print out possible values of properties.
if (nargin == 0) && (nargout == 0)
    printfields();
    return
end
%Names and Defaults
Names = {'drho';'rpos';'r0';'penalty';'rho';'dbound';'gaptol';'rtol';'mu0';'maxtrust';...
         'steptol';'ptol';'pnormtol';'reuse';'zbar';'dlbound';'ybound';'fixed'};
Defaults = {1;0;-1;1e8;5;1e20;1e-7;1e-6;-1;1e10;0.05;1e-4;1e30;4;[];[];[];[]};

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case {'r0','penalty','rho','mu0','dbound','pnormtol','zbar','dlbound','ybound'}
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar non negative double
    case {'rtol','ptol','maxtrust','gaptol','steptol'}
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar 0/1
    case {'rpos','drho'}
        err = opticheckval.checkScalar01(value,field);
    %Scalar non-negative integer
    case 'reuse'
        err = opticheckval.checkScalarIntNonNeg(value,field);
    %Double column vector
    case 'fixed'
        err = opticheckval.checkDualCol(value,field);
    otherwise
        err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf('          rpos: [ Use Penalty Parameter to enforce Feasibility {0}, 1 ] \n');
fprintf('          drho: [ Use Dynamic Strategy to choose Rho: 0, {1} ] \n');
fprintf('            r0: [ Set initial value for variable r in Dual {-1} ] \n');
fprintf('       penalty: [ Penalty Parameter Gamma {1e8} \n');
fprintf('           rho: [ Potential Parmater: {5} ] \n');
fprintf('           mu0: [ Barrier Parameter: {-1} ] \n');
fprintf('          rtol: [ Classify Dual as feasible if r is less than this tolerance: {1e-6} ] \n');
fprintf('          ptol: [ Classify Primal as feasible if infeasibility is less than this tolerance: {1e-4} ] \n');
fprintf('      maxtrust: [ Maximum Trust Radius on Step Direction: {1e10} ] \n');
fprintf('        gaptol: [ Convergence Gap Tolerance: {1e-7} ] \n');
fprintf('       steptol: [ Terminate solver if step length in Dual is below this tolerance: {0.05} ] \n');
fprintf('        dbound: [ Terminate solver if Dual Objective greater than this Value {1e20} ] \n');
fprintf('      pnormtol: [ Terminate the solver when duality gap is sufficiently small and pnorm is less than this quantity: {1e30} ] \n');
fprintf('         reuse: [ Number of times the Hessian of the Barrier Function will be reused {4} ] \n');
fprintf('          zbar: [ Upper bound on objective value at the solution {[]} ] \n');
fprintf('       dlbound: [ Lower bound on dual value {[]} ] \n');
fprintf('        ybound: [ Bound on the dual variables y {[]} ] \n');
fprintf('         fixed: [ Matrix of fixed variables, column 1 variable indices, column 2, fixed values {[]} ] \n');
fprintf('\n');

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef ContourFitting < handle
    %CONTOURFITTING  Contour Fitting algorithm using Fourier descriptors
    %
    % Contour fitting matches two contours `z_a` and `z_b` minimizing distance
    % `d(z_a, z_b) = sum_n (a_n - s * b_n * exp(j * (n * alpha + phi)))^2`
    % where `a_n` and `b_n` are Fourier descriptors of `z_a` and `z_b` and `s`
    % is a scaling factor and `phi` is angle rotation and `alpha` is starting
    % point factor adjustement.
    %
    % ## References:
    % [PersoonFu1977]:
    % > E Persoon and King-Sun Fu. "Shape discrimination using fourier
    % > descriptors". IEEE Transactions on Pattern Analysis and Machine
    % > Intelligence, 7(3):170-179, 1977.
    %
    % [BergerRaghunathan1998]:
    % > L Berger, V A Raghunathan, C Launay, D Ausserre, and Y Gallot.
    % > "Coalescence in 2 dimensions: experiments on thin copolymer films and
    % > numerical simulations". The European Physical Journal B - Condensed
    % > Matter and Complex Systems, 2(1):93-99, 1998.
    %
    % See also: cv.ContourFitting.ContourFitting, cv.matchShapes,
    %  cv.ShapeContextDistanceExtractor
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % number of Fourier descriptors used in
        % cv.ContourFitting.estimateTransformation equal to number of contour
        % points after resampling.
        CtrSize
        % number of Fourier descriptors used for optimal curve matching in
        % cv.ContourFitting.estimateTransformation when using vector of points
        FDSize
    end

    %% ContourFitting
    methods
        function this = ContourFitting(varargin)
            %CONTOURFITTING  Create ContourFitting object
            %
            %     obj = cv.ContourFitting()
            %     obj = cv.ContourFitting('OptionName',optionValue, ...)
            %
            % ## Options
            % * __CtrSize__ number of contour points after resampling.
            %   default 1024
            % * __FDSize__ number of Fourier descriptors. default 16
            %
            % See also: cv.ContourFitting.estimateTransformation
            %
            this.id = ContourFitting_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.ContourFitting
            %
            if isempty(this.id), return; end
            ContourFitting_(this.id, 'delete');
        end

        function [alphaPhiST, d] = estimateTransformation(this, src, ref, varargin)
            %ESTIMATETRANSFORMATION  Fits two closed curves using Fourier descriptors
            %
            %     [alphaPhiST, d] = obj.estimateTransformation(src, ref)
            %     [...] = obj.estimateTransformation(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __src__ Contour defining first shape (source), or Fourier
            %   descriptors if `FD` is true.
            % * __ref__ Contour defining second shape (target), or Fourier
            %   descriptors if `FD` is true.
            %
            % ## Output
            % * __alphaPhiST__ transformation as a 5-elements vector
            %   `[alpha, phi, s, Tx, Ty]`, where:
            %   * __alpha__ starting point factor adjustement
            %   * __phi__ angle rotation in radian
            %   * __s__ scaling factor
            %   * __Tx__, __Ty__ the translation
            % * __d__ distance between `src` and `ref` after matching.
            %
            % ## Options
            % * __FD__ If false then `src` and `ref` are contours, and
            %   if true `src` and `ref` are Fourier descriptors. default false
            %
            % When `FD` is false, it applies cv.ContourFitting.contourSampling
            % and cv.ContourFitting.fourierDescriptor to compute Fourier
            % descriptors.
            %
            % More details in [PersoonFu1977] and [BergerRaghunathan1998].
            %
            % See also: cv.ContourFitting.transformFD,
            %  cv.ContourFitting.contourSampling,
            %  cv.ContourFitting.fourierDescriptor
            %
            [alphaPhiST, d] = ContourFitting_(this.id, 'estimateTransformation', src, ref, varargin{:});
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.ContourFitting.empty, cv.ContourFitting.load
            %
            ContourFitting_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if algorithm object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the algorithm object is empty
            %   (e.g. in the very beginning or after unsuccessful read).
            %
            % See also: cv.ContourFitting.clear, cv.ContourFitting.load
            %
            b = ContourFitting_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.ContourFitting.load
            %
            ContourFitting_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.ContourFitting.save
            %
            ContourFitting_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.ContourFitting.save, cv.ContourFitting.load
            %
            name = ContourFitting_(this.id, 'getDefaultName');
        end
    end

    %% Getters/Setters
    methods
        function value = get.CtrSize(this)
            value = ContourFitting_(this.id, 'get', 'CtrSize');
        end
        function set.CtrSize(this, value)
            ContourFitting_(this.id, 'set', 'CtrSize', value);
        end

        function value = get.FDSize(this)
            value = ContourFitting_(this.id, 'get', 'FDSize');
        end
        function set.FDSize(this, value)
            ContourFitting_(this.id, 'set', 'FDSize', value);
        end
    end

    %% Static functions
    methods (Static)
        function out = contourSampling(src, numElt)
            %CONTOURSAMPLING  Contour sampling
            %
            %     out = cv.ContourFitting.contourSampling(src, numElt)
            %
            % ## Input
            % * __src__ input contour, vector of 2D points stored in numeric
            %   array Nx2/Nx1x2/1xNx2 or cell array of 2-element vectors
            %   `{[x,y], ...}`.
            % * __NumElt__ number of points in `out` contour.
            %
            % ## Output
            % * __out__ output contour with `numElt` points.
            %
            % See also: cv.findContours, cv.approxPolyDP
            %
            out = ContourFitting_(0, 'contourSampling', src, numElt);
        end

        function dst = fourierDescriptor(src, varargin)
            %FOURIERDESCRIPTOR  Fourier descriptors for planed closed curves
            %
            %     dst = cv.ContourFitting.fourierDescriptor(src)
            %     dst = cv.ContourFitting.fourierDescriptor(src, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __src__ input contour, vector of 2D points stored in numeric
            %   array Nx2/Nx1x2/1xNx2 or cell array of 2-element vectors
            %   `{[x,y], ...}`.
            %
            % ## Output
            % * __dst__ 2-channel array of type `single` and length `NumElt`.
            %
            % ## Options
            % * __NumElt__ number of rows in `dst` or cv.getOptimalDFTSize
            %   rows if `NumElt=-1`. default -1
            % * __NumFD__ number of FD to return in `dst`,
            %   `dst = [FD(1...NumFD/2) FD(NumFD/2-NumElt+1...:NumElt)]`.
            %   default -1 (return all of FD as is).
            %
            % For more details about this implementation, please see
            % [PersoonFu1977].
            %
            % See also: cv.dft
            %
            dst = ContourFitting_(0, 'fourierDescriptor', src, varargin{:});
        end

        function dst = transformFD(src, t, varargin)
            %TRANSFORMFD  Transform a contour
            %
            %     dst = cv.ContourFitting.transformFD(src, t)
            %     dst = cv.ContourFitting.transformFD(src, t, 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __src__ contour, or Fourier descriptors if `FD` is true.
            % * __t__ 1x5 transform matrix given by
            %   cv.ContourFitting.estimateTransformation method.
            %
            % ## Output
            % * __dst__ 2-channel matrix of type `double` and `NumElt` rows.
            %
            % ## Options
            % * __FD__ if true `src` are Fourier descriptors, if false `src`
            %   is a contour. default true
            %
            % See also: cv.ContourFitting.estimateTransformation,
            %  cv.ContourFitting.contourSampling,
            %  cv.ContourFitting.fourierDescriptor
            %
            dst = ContourFitting_(0, 'transformFD', src, t, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function blockface_correct_perspective(indir, files55, tform, outdir)
% blockface_correct_perspective  Correct perspective of Bewster angle
% blockface images.
%
% blockface_correct_perspective(INDIR, FILES55, TFORM, OUTDIR)
%
%   INDIR is a string with the directory where the input files are kept.
%
%   FILES55 is the result of a dir() command with a list of the files to
%   correct.
%
%   TFORM is a projective2d object with the transform correction for the
%   frames.
%
%   OUTDIR is a string with the directory for output images. INDIR and
%   OUTDIR must be different to avoid overwriting the input files.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.2.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/>.

% check arguments
narginchk(4, 4);
nargoutchk(0, 0);

% safety to prevent us from overwriting the source files
if (strcmp(indir, outdir))
    error('Input and output directory are the same, and files would be overwritten')
end

% if destination directory doesn't exist, create it
if (isempty(dir(outdir)))
    ok = mkdir(outdir);
    if (~ok)
        error(['Cannot create output directory: ' outdir])
    end
end

% correct perspective
for I = 1:length(files55)

    % load image
    im55 = imread([indir filesep files55(I).name]);

    % make output image the same size and coordinates as input image
    ra = imref2d(size(im55));

    % correct perspective
    im55to90 = imwarp(im55, tform, 'cubic', ...
        'OutputView', ra, 'FillValues', double(min(im55(:))));

%     % plot results
%     subplot(2, 1, 1)
%     imagesc(im55)
%     subplot(2, 1, 2)
%     imagesc(im55to90)

    % save result
    imwrite(im55to90, [outdir filesep files55(I).name]);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


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

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function view=initDataValsSlider(view,curDVIndex)
%
%       view=initDataValsSlider(view,[curDVIndex])
%
% Installs the slider and text box that control the current datavals .
% The slider is placed on the mid-left, under the datavals control slider.
%
% The first time through, the callback and related textbox are installed
%
% It is also called repeatedly as dataTypes or FFTCheckBox are switched to re-display
%   the slider that represents datavals information
%
% ARW Wrote it 030606:
% Based on BW's initScanSlider

% The sliderHandle always exists because it is created in openXXWindow
% Set the step size for clicking on the arrows and the trough
sliderHandle = view.ui.datavals.sliderHandle;

nDataVals=numDataVals(view); % This looks at dataTYPES and also FFTCheckBox to see if we're currently displaying FFT or Time-domain data
if nDataVals < 2
    %   Can't have a slider value of 1.  Sorry.
    nDataVals=2;
    set(sliderHandle,'Enable','off');
else
    set(sliderHandle,'Enable','on');
end

% By setting the max to a little more than nDataValss, we should
% never quite reach it by clicking the arrow or the trough
set(sliderHandle,'Min',1,'Max',nDataVals);

% Clicking an arrow moves by an amount (max - min)*sliderstep(1)
% Cllicking in the trough moves by an amount (max - min)*sliderstep(2)
% setSlider traps any values == to the max value and queries the user
%    this is a cheap way of setting a value by hand.
sliderStep(1) = 1/(nDataVals - 1);
sliderStep(2) = 2*sliderStep(1);
set(sliderHandle,'sliderStep',sliderStep);

% If we are first opening the window, create the text box.
if ~isfield(view.ui.datavals,'textHandle')
    % Can this be inside the if below so that it only happens
	% when the window is first opened?
	% djh, 9/5/01
	callBackStr = ...
        ['val = get(',view.name,'.ui.datavals.sliderHandle,''value'');'...
            'val = round(val);'...
            'setSlider(',view.name,',',view.name,'.ui.datavals,val);'...
            view.name '=setCurDataValIndex(',view.name,',val);'...
            view.name,'=refreshScreen(',view.name,');'];
	set(sliderHandle,'CallBack',callBackStr);

   % Position the text box to the right of the slider
   pos = get(sliderHandle,'Position');
   l2 = pos(1) + pos(3) + .01;
   w2 = pos(3)/4; b = pos(2); h = pos(4);
   position = [l2 b w2 h];

   view.ui.datavals.textHandle = ...
      uicontrol('Style','text',...
      'Units','normalized',...
      'Position',position);
end

% This text tells the user how many datavals indices there are
% in the current dataType and domain.  Here, we set the value
%
str = sprintf('%.0f',nDataVals);
set(view.ui.datavals.textHandle,'String',str);

% Synchronize the curScan with the slider state
%

if (~isfield(view,'dataValIndex'))
    view = setCurDataValIndex(view,1);
end

if ~exist('curDataValIndex','var')
   curDataValIndex = getCurDataValIndex(view);
end

if curDataValIndex > numDataVals(view)
    curDataValIndex = 1;
end
setSlider(view,view.ui.datavals,curDataValIndex);
view = setCurDataValIndex(view,curDataValIndex);

% ras 08/05: make the scan slider update the curScan field as well:
cb = 'val=str2num(get(gcbo,''String'')); ';
cb = [cb sprintf('%s=setCurDataValIndex(%s,val); ',view.name,view.name)];
cb = [cb sprintf('%s=refreshScreen(%s);',view.name,view.name)];
set(view.ui.datavals.labelHandle,'Callback',cb);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%returns the first component of the frequency representation of the signal
classdef SpectralSpread < Algorithm

    methods (Access = public)

        function obj = SpectralSpread()
            obj.name = 'SpectralSpread';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        %receives a fourier transform
        function result = compute(~,Y)
            N = length(Y);
            Y = abs(Y) / N;
            Y = Y / max(Y);

            m = ((1/(2*N))*(1:N))';

            % compute the spectral spread
            C = sum(m .* Y) / (sum(Y)+eps);
            result = sqrt(sum(((m-C).^2) .* Y)/ (sum(Y)+eps));
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 11 * n;
            memory = n;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
classdef TestGlob
    %TestGlob

    methods (Static)
        function test_1
            pattern = fullfile(mexopencv.root(),'test','left*.jpg');
            files = cv.glob(pattern, 'Recursive',false);
            validateattributes(files, {'cell'}, {'vector'});
            assert(iscellstr(files));
        end

        function test_error_argnum
            try
                cv.glob();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
function [diet] = makeDietFormatConsistant(model,diet)
% This function takes a diet input and checks it for formatting consisteny
% for input into setFoodConstraints
%
% USAGE:
%
%   [diet] = makeDietFormatConsistant(model,diet)
%
% INPUTS
%   model:             A COBRA model
%
%   diet:            A nx2 cell array containing n dietary components and
%                    the corresponding flux
% OUTPUT
%   diet:           A diet of consistent format for setFoodConstraints
%
% Authors:
%   Bronson R. Weston 2022

load fdTable.mat
load fdCategoriesTable.mat
foods=[fdTable.Properties.VariableNames(2:end),fdCategoriesTable.Properties.VariableNames(2:end)];
for i=1:length(diet(:,1))
    if ~isempty(find(strcmp(foods,diet{i,1})))
        diet{i,1}=['Food_EX_',diet{i,1},'[d]'];
    elseif ~isempty(find(strcmp(model.rxns,['Diet_EX_' diet{i,1} '[d]'])))
                disp('C')
        diet{i,1}=['Diet_EX_',diet{i,1},'[d]'];
    end
end
end

author: opencobra

repo: cobratoolbox


text:
function test_failed=test_spread
%TEST_SPREAD  Test spreading function.
%
%  This script runs a throrough test of the SPREADOP routine,
%  testing it on a range of input parameters.
global LTFAT_TEST_TYPE;
disp(' ===============  TEST_SPREAD ================');

disp('--- Used subroutines ---');

which comp_col2diag

Lr=[12, 13 ];

test_failed=0;

condNoLim.('single') = 100;
condNoLim.('double') = 100;

for ii=1:length(Lr);

L=Lr(ii);

spfraction=.5;

for rtype=1:2

  if rtype==1
    rname='REAL ';
  else
    rname='CMPLX';
  end;

  for sptype=1:2

    if sptype==1
      spname='FULL  ';
    else
      spname='SPARSE';
      if ~strcmpi(LTFAT_TEST_TYPE,'double')
          disp(sprintf('Skipping. Cannot work with sparse matrices of type %s.',LTFAT_TEST_TYPE));
          break;
      end
    end;

    condNo = 1e10;
    while condNo>condNoLim.(LTFAT_TEST_TYPE)
    if rtype==1
      if sptype==1
        coef=tester_rand(L,L);
        T=tester_rand(L,L);
        coef2=tester_rand(L,L);
      else
        coef=tester_sprand(L,L,spfraction);
        T=tester_sprand(L,L,spfraction);
        coef2=tester_sprand(L,L,spfraction);
      end;
    else
      if sptype==1
        coef=tester_crand(L,L);
        T=tester_crand(L,L);
        coef2=tester_crand(L,L);
      else
        coef=tester_sprand(L,L,spfraction);
        T=tester_sprand(L,L,spfraction);
        coef2=tester_sprand(L,L,spfraction);
      end;
    end;

      coeftmp=ifft(full(coef))*L;
      % The following matrix is inverted in spreadinv. We want a nicer cond
      % number in order not to fail
      Ttmp=comp_col2diag(coeftmp);
      condNo = cond(Ttmp);
    end

    for W=1:3

      if rtype==1
        f=tester_rand(L,W);
      else
        f=tester_crand(L,W);
      end;

      % ---------- Reference testing ------------

      fs=spreadop(f,coef);
      fs2=ref_spreadop(f,coef,1);

      fsdiff=fs-fs2;
      res=norm(fsdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf('REF %s %s L:%3i W:%2i %0.5g %s',rname,spname,L,W,res,fail);
      disp(s)

      % -------- Inversion testing -------------

      r=spreadinv(fs,coef);

      rdiff=f-r;
      res=norm(rdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf('INV %s %s L:%3i W:%2i %0.5g %s',rname,spname,L,W,res,fail);
      disp(s)

      % -------- Twisted convolution -------------------

      f1=spreadop(spreadop(f,coef2),coef);
      f2=spreadop(f,tconv(coef,coef2));

      rdiff=f1-f2;
      res=norm(rdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf('TWI %s %s L:%3i W:%2i %0.5g %s',rname,spname,L,W,res,fail);
      disp(s)

      % -------- Spreading function ---------------------

      coef=spreadfun(T);

      rdiff=T*f-spreadop(f,coef);

      res=norm(rdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf('FUN %s %s L:%3i W:%2i %0.5g %s',rname,spname,L,W,res,fail);
      disp(s)

      % -------- Adjoint operator -----------------------

      cadj=spreadadj(coef);

      rdiff=T'*f-spreadop(f,cadj);

      res=norm(rdiff(:));
      [test_failed,fail]=ltfatdiditfail(res,test_failed);
      s=sprintf('ADJ %s %s L:%3i W:%2i %0.5g %s',rname,spname,L,W,res,fail);
      disp(s)

    end;

  end;

end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [afixed,sqnorm,ierr] = lsearch (afloat,L,D,Chi2,ncands)
%LSEARCH: Integer ambiguity resolution, search
%
% This routine finds the integer vector which is closest to a given
% float vector, in a least squares sence. This is the search-step in
% integer ambiguity resolution. It is best to perform this search only
% on ambiguities which have been decorrelated using LAMBDA.
%
% Input arguments:
%    afloat : Float ambiguities (\hat{a})
%    L      : LtDL-decomposition of the decorrelated
%    D        variance-covariance matrix of the ambiguities
%    Chi2   : Size of the search ellipsoid
%    ncands : Number of requested candidates
%
% Output arguments:
%    afixed : Estimated integers (matrix)
%    sqnorm : Corresponding squared norms (vector, sorted)
%    ierr   : Error code: 0: No errors found
%                         1: Not enough candidates found

% ----------------------------------------------------------------------
% File.....: lsearch.m
% Date.....: 19-MAY-1999
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

% -------------------------------
% --- Initializing statements ---
% -------------------------------

Linv      = inv(L);
Dinv      = 1./D;

True      = 1;
False     = 0;

n         = max(size(afloat));

right     = [zeros(n,1) ; Chi2];
left      = [zeros(n+1,1)];
dq        = [Dinv(2:n)./Dinv(1:n-1) 1/Dinv(n)];

cand_n    = False;
c_stop    = False;
endsearch = False;

ncan      = 0;

i         = n + 1;
iold      = i;
ierr      = 0;

afixed = zeros(n,ncands);
sqnorm = zeros(1,ncands);

dist = [];

% ----------------------------------
% --- Start the main search-loop ---
% ----------------------------------

while ~ (endsearch);

   i = i - 1;

   if iold <= i
      lef(i) = lef(i) + Linv(i+1,i);
   else
      lef(i) = 0;
      for j = i+1:n;
         lef(i) = lef(i) + Linv(j,i)*dist(j,1);
      end;
   end;
   iold = i;

   right(i)  = (right(i+1) - left(i+1)) * dq(i);
   reach     = sqrt(right(i));
   delta     = afloat(i) - reach - lef(i);
   dist(i,1) = ceil(delta) - afloat(i);

   if dist(i,1) > reach - lef(i)

%     ----------------------------------------------------
%     --- There is nothing at this level, so backtrack ---
%     ----------------------------------------------------

      cand_n = False;
      c_stop = False;

      while (~ c_stop) & (i < n);

         i = i + 1;
         if dist(i) < endd(i);
            dist(i) = dist(i) + 1;
            left(i) = (dist(i) + lef(i)) ^ 2;
            c_stop = True;
            if i == n; cand_n = True; end;
         end;

      end;

      if (i == n) & (~ cand_n); endsearch = True; end;

   else

%     ----------------------------
%     --- Set the right border ---
%     ----------------------------

      endd(i) = reach - lef(i) - 1;
      left(i) = (dist(i,1) + lef(i)) ^ 2;

   end

   if i == 1;

%     -------------------------------------------------------------------
%     --- Collect the integer vectors and corresponding               ---
%     --- squared distances, add to vectors "afixed" and "sqnorm" if: ---                             ---
%     --- * Less then "ncands" candidates found so far                ---
%     --- * The squared norm is smaller than one of the previous ones ---
%     -------------------------------------------------------------------

      t       = Chi2 - (right(1)-left(1)) * Dinv(1);
      endd(1) = endd(1) + 1;

      while dist(1) <= endd(1);

         if ncan < ncands;

            ncan             = ncan + 1;
            afixed(1:n,ncan) = dist + afloat;
            sqnorm(ncan)     = t;

         else

            [maxnorm,ipos] = max(sqnorm);
            if t < maxnorm;
               afixed(1:n,ipos) = dist + afloat;
               sqnorm(ipos)     = t;
            end;

         end;

         t       = t + (2 * (dist(1) + lef(1)) + 1) * Dinv(1);
         dist(1) = dist(1) + 1;

      end;

%     -------------------------
%     --- And backtrack ... ---
%     -------------------------

      cand_n = False;
      c_stop = False;

      while (~ c_stop) & (i < n);

         i = i + 1;

         if dist(i) < endd(i);
            dist(i) = dist(i) + 1;
            left(i) = (dist(i) + lef(i)) ^ 2;
            c_stop = True;
            if i == n; cand_n = True; end;
         end;

      end;

      if (i == n) & (~ cand_n); endsearch = True; end;

   end;

end;

% ----------------------------------------------------------------------
% --- Sort the resulting candidates, according to the norm
% ----------------------------------------------------------------------

tmp    = sortrows ([sqnorm' afixed']);
sqnorm = tmp(:,1)';
afixed = tmp(:,2:n+1)';

% ------------------------
% --- Check for errors ---
% ------------------------

if ncan < ncands; ierr = 1; end;

% ----------------------------------------------------------------------
% End of routine: lsearch
% ----------------------------------------------------------------------

author: yandld


text:
function table = distableAll(board)

n = numel(board);
table = NaN * zeros(n);

for i = 1:n
    table(i,i) = 0;
    table(i, board(i).Neighbors) = 1;
    d = 1;
    N = -1;
    while nnz(~isnan(table(i,:))) ~= N;
        N = nnz(~isnan(table(i,:)));
        for j = find(table(i,:) == d)
            for k = board(j).Neighbors
                if i ~= k && isnan(table(i, k))
                    table(i, k) = d + 1;
                end
            end
        end
        d = d + 1;
    end
end


text:
function y = daub16_transform ( n, x )

%*****************************************************************************80
%
%% DAUB16_TRANSFORM computes the DAUB16 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2 and at least 4.
%
%    Input, real X(N), the vector to be transformed.
%
%    Output, real Y(N), the transformed vector.
%
  c = [ ...
     5.441584224310400E-02; ...
     3.128715909142999E-01; ...
     6.756307362972898E-01; ...
     5.853546836542067E-01; ...
    -1.582910525634930E-02; ...
    -2.840155429615469E-01; ...
     4.724845739132827E-04; ...
     1.287474266204784E-01; ...
    -1.736930100180754E-02; ...
    -4.408825393079475E-02; ...
     1.398102791739828E-02; ...
     8.746094047405776E-03; ...
    -4.870352993451574E-03; ...
    -3.917403733769470E-04; ...
     6.754494064505693E-04; ...
    -1.174767841247695E-04 ];
  p = 15;
  y(1:n,1) = x(1:n);
  m = n;
  q = floor ( ( p - 1 ) / 2 );

  while ( 4 <= m )

    i = 1;
    z(1:m,1) = 0.0;

    for j = 1 : 2 : m - 1

      mh = floor ( m / 2 );
      for k = 0 : 2 : p - 1
        j0 = i4_wrap ( j + k,     1, m );
        j1 = i4_wrap ( j + k + 1, 1, m );
        z(i,1)    = z(i,1)    + c(  k+1) * y(j0) + c(  k+2) * y(j1);
        z(i+mh,1) = z(i+mh,1) + c(p-k+1) * y(j0) - c(p-k  ) * y(j1);
      end

      i = i + 1;

    end

    y(1:m,1) = z(1:m);

    m = floor ( m / 2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lattice_rule_test10 ( )

%*****************************************************************************80
%
%% LATTICE_RULE_TEST10 tests LATTICE_NP1;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994, page 78-80, pages 32-40, 145-147.
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATTICE_RULE_TEST10\n' );
  fprintf ( 1, '  LATTICE_NP1 applies a lattice rule to a\n' );
  fprintf ( 1, '  nonperiodic function using a nonlinear transformation,\n' );
  fprintf ( 1, '  to integrate a function over the unit square.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The spatial dimension DIM_NUM = %d\n', dim_num );

  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  exact = e_01_2d ( dim_num, a, b );

  z(1:dim_num) = [ 1, 2 ];
  i4vec_print ( dim_num, z, '  The lattice generator vector:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '         M    EXACT       ESTIMATE    ERROR\n' );
  fprintf ( 1, '\n' );

  for k = 3 : 18

    m = fibonacci ( k );

    quad = lattice_np1 ( dim_num, m, z, @f_01_2d );

    error = abs ( exact - quad );

    fprintf ( 1, '  %8d  %10.6f  %10.6f  %10.6e\n', m, exact, quad, error );

  end

  return
end


text:
function [pX,gX,pY,gY,X,Y,U] = get_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb,N,np,lx,ly)

% legacy code
s = warning ('on');
warning ('*** The function `get_MCMC_predictiveDensity_fb` is now deprecated. Please use `VBA_MCMC_predictiveDensity_fb` instead (same syntax).')
warning (s);

% fallback
switch nargin
    case 11
        [pX,gX,pY,gY,X,Y] = VBA_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb,N,np,lx,ly);
    case 10
        [pX,gX,pY,gY,X,Y] = VBA_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb,N,np,lx);
    case 9
        [pX,gX,pY,gY,X,Y] = VBA_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb,N,np);
    case 8
        [pX,gX,pY,gY,X,Y] = VBA_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb,N);
    case 7
        [pX,gX,pY,gY,X,Y] = VBA_MCMC_predictiveDensity_fb(f_fname,g_fname,u,n_t,options,dim,fb);
    otherwise
        error('VBA_MCMC_predictiveDensity_fb: wrong number of arguments');
end

author: MBB-team

repo: VBA-toolbox


text:
function hrf = boyntonHIRF(t, n, tau, delay)
%
%    hrf = boyntonHIRF(t, [n=3], [tau=1.08], [delay=2.05]])
%
%Purpose:
%   Compute the Boynton et al. HIRF function.  The return values include
% both the hrf and the values of time and the parameters used in the
% computation.
%
% Equation:  (Eq 3 from Boynton & Heeger, J Neurosci 1996)
%   h(t) = [(t/tau) ^ (n-1) * exp(-t/tau)] / [tau(n-1)!]
%
% Inputs:
%   t: t should be in units of SECONDS, and reflect the time window in
%   which to estimate the HRF.
%
%   n: exponent in eq. above. (corresponds to z in standard gamma
%   functions).
%
%   tau: time constant in eq.
%
%   delay: additional delay before onset of gamma function. This is an
%   added heuristic, which may vary from subject to subject.
%
%
% Outputs:
%   hrf: estimate of the HRF sampled at t seconds.
%   IMPORTANT: If you are going to convolve
%   this for fMRI analysis, and the TR of your data is not 1, you will
%   need to resample both the input t and hrf to match the MR frames.
%
% written 2005 by wandell.
% ras, 01/2007: heavily modified: no parms struct, each arg can be
% specified separately; clarifies difference between seconds and frames,
% doesn't modify t.
if notDefined('n'),     n = 3;          end
if notDefined('tau'),   tau = 1.08;     end
if notDefined('delay'), delay = 2.05;   end

% ras 01/07: this is sometimes nice, but let's just turn it off
% altogether for now
% verbose = prefsVerboseCheck;
% if verbose,
%     disp('Boynton HIRF')
% end

% initialize the HRF to be zeros, the same size as t:
hrf = zeros( size(t) );

% The HRF is not specified for t < 0 secs. In addition, any
% values before the delay should also be zero.
% So, we only sample the HRF below, for time points after [0+delay].
% We call this sampling vector x to distinguish it from t.
x = t - delay;
x = x(x>0);  % not defined for x<0

% main computation (per Boynton & Heeger, 1996)
tmpHrf = (x/tau).^(n-1) .* exp(-(x/tau)) / (tau*(factorial(n-1)));

% paste in tmpHRF into appropriate indices in hrf, corresponding
% to the (shifted after delay) time points:
hrf(end-length(tmpHrf)+1:end) = tmpHrf;

return;


text:
function [x_pred, y_pred]= convert_trans_and_scale(pred_point_raw, mid_point, scale)
   % 1. The input heatmap is 248 by 248, which is rescaled from a crop, introducing
   %    a scale
   % 2. The crop comes from the original image, and introduce a translation
   % Summary: This script is to do the backwards mapping, back to the original position.
   % (1) First, find the position in the cropped image
   % (2) Second, find the position in the original image

   % Use (scale) for the scale. [180, 280] is current value.
   % Consult [gen_cropped_test_images.m]. Make sure use same scale.
   % Officially, scale is w.r.t. 200 px height
   human_bbox_wid = 180 * scale;
   human_bbox_ht = 280 * scale;

   % Prepare [x_min, y_min] for translation
   x_min = mid_point(1) - human_bbox_wid/2.0;
   y_min = mid_point(2) - human_bbox_ht/2.0;

   % Based on translation, scale, and output, find input
   [x_pred, y_pred] = find_pos_for_org_image(pred_point_raw, human_bbox_wid, human_bbox_ht, x_min, y_min);
end

function [x0, y0] = find_pos_for_org_image(pos_248, human_bbox_wid, human_bbox_ht, x_min, y_min)
    [x1, y1] = find_pos_for_cropped_image(pos_248, human_bbox_wid, human_bbox_ht);

    x0 = x1 + x_min;
    y0 = y1 + y_min;

    pos_org = [x0, y0];
end

function [x1, y1] = find_pos_for_cropped_image(pos_248, human_bbox_wid, human_bbox_ht)
   x2 = pos_248(1);
   y2 = pos_248(2);

   x1 = int64((x2 * human_bbox_wid)/248.0);
   y1 = int64((y2 * human_bbox_ht)/248.0);

   pos_cropped = [x1, y1];
end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
% ------------------------------------------------------------------------
% Copyright (C) 2013 University of Southern California, SAIL, U.S.
% Author: Maarten Van Segbroeck
% Mail: maarten@sipi.usc.edu
% Date: 2013-28-10
% ------------------------------------------------------------------------
% ltsv = FE_LTSV(sam,fs,R,M,ltsvThr,ltsvSlope);
%      Generates a Long-Term Spectral Variability (LTSV) stream of speech compute on the
%      Gammatone filtered time-frequency representation of speech (different from original
%      LTSV stream in [1]). The LTSV is also processed by a sigmoid function to assure
%      a probability value between 0 and 1 on the current signal
%
%      --IN--
%      sam: sample file
%      fs: sampling frequency
%      R: context window parameter [default:50
%      M: smoothing parameter [default:M]
%      ltsvThr: threshold of sigmoid function [default:0.5]
%      ltsvSlope: slope of sigmoid function [default:0.2]
%
%      --OUT--
%      ltsv: LTSV stream of the signal

function ltsv=FE_LTSV(sam,fs,R,M,gt,ltsvThr,ltsvSlope);

FrameLen=0.032*fs;
FrameShift=0.01*fs;
Nfft=2^ceil(log(FrameLen)/log(2));
Nfft=128;
S=FE_GT(sam,fs,Nfft,FrameLen);

if ~exist('R', 'var')
 R=50;
end
if ~exist('M', 'var')
 M=10;
end
if ~exist('ltsvThr', 'var')
 ltsvThr=0.5;
end
if ~exist('ltsvSlope', 'var')
 ltsvSlope=0.2;
end

append_ndx=size(S,2):-1:max(size(S,2)-20,1);
S=[S S(:,append_ndx)];

% frequency smoothing
S_M=[];
for k=1:size(S,1)
	S_M(k,:)=(smooth(S(k,:),M,'moving'));
end
% normalized spectrogram
S_R=[];
for t=1:size(S,2)
	ndx=[ones(1,R/2-t) max(1,t-R/2+1):min(size(S,2),t+R/2) ones(1,R/2+t-size(S,2))];
%	ndx=[t:min(size(S,2),t+R) ones(1,R+t-size(S,2))];
	S_R(:,t)=S_M(:,t)./sum(S_M(:,ndx),2);
end

% entropy measure of normalized spectrogram over R consecutive frames, ending at current frame
E_R=-100*S_R.*log(100*S_R);
L=var(E_R);

L=[];
N=[1];
skip=0;
for n=1:length(N)
  for i=1:N(n)
        fndx=round([Nfft*(i-1)/(2*N(n))+1:Nfft*i/(2*N(n))]);
	fndx=unique(min(max(fndx,1),Nfft));
	L(i+skip,:)=var(E_R(fndx,:));
  end
  skip=skip+N(n);
end

S=S(:,1:end-length(append_ndx));
L=L(:,1:end-length(append_ndx));
if exist('ltsvThr', 'var') && exist('ltsvSlope', 'var')
 L=sigmoid(2*L-ltsvThr,ltsvSlope);
end
ltsv=L;

author: jtkim-kaist

repo: VAD

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


text:
%COPY_GEOMETRY_OBJECT Copy geometry object(s) with transformation.
%
%   [ GEOM, NEW_TAGS ] = COPY_GEOMETRY_OBJECT( TAGS, GEOM, TRSFARGS,
%   CONNECT ) Copies the geometry object(s) in GEOM with corresponding
%   TAGS. Optional TRSFARGS are passed to transform the objects, and
%   CONNECTing boundaries for 3D CAD geometries (default false).
%
%   See also GEOM_APPLY_TRANSFORMATION

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function q = project2FR_ref(q,qCS,q_ref)
% projects quaternions to a fundamental region
%
% Syntax
%   project2FundamentalRegion(q,CS,q_ref) % to FR around reference rotation
%
% Input
%  q        - @quaternion
%  CS1, CS2 - crystal @symmetry
%  q_ref    - reference @quaternion single or size(q) == size(q_ref)
%
% Output
%  q     - @quaternion
%  omega - rotational angle to reference quaternion
%

qCS = qCS.rot;

s = size(q);
q.a = q.a(:); q.b = q.b(:); q.c = q.c(:); q.d = q.d(:);
try q.i = q.i(:); end %#ok<TRYNC>

% compute distance to reference orientation
if ~isempty(q_ref)
  q_ref.a = q_ref.a(:); q_ref.b = q_ref.b(:);
  q_ref.c = q_ref.c(:); q_ref.d = q_ref.d(:);

  co2 = abs(quat_dot(q,q_ref));

else

  co2 = abs(q.a);

end

% may be we can skip something
minAngle = reshape(abs(qCS.angle),[],1);
minAngle = min([inf;minAngle(minAngle > 1e-3)]);
notInside = co2 < cos(minAngle/4);

% maybe we can skip everything
if any(notInside) && length(qCS) > 1

  % restrict to quaternion which are not yet it FR
  if length(q) == numel(notInside)
    q_sub = q.subSet(notInside);
  else
    q_sub = q;
  end

  % compute all distances to the fundamental regions
  if ~isempty(q_ref)

    % if q_ref was a list of reference rotations
    if length(q_ref) == numel(notInside)
      omegaSym  = abs(quat_dot_outer(inv(q_ref.subSet(notInside)).*q_sub,qCS));
    else
      omegaSym  = abs(quat_dot_outer(inv(q_ref).*q_sub,qCS));
    end

  else
    omegaSym  = abs(quat_dot_outer(q_sub,qCS));
  end

  % find symmetry elements with minimum distance
  [~,nx] = max(omegaSym,[],2);

  % project to fundamental region
  qn = times(q_sub, reshape(inv(qCS.subSet(nx)),size(q_sub)),0);

  % replace projected quaternions
  q.a(notInside) = qn.a;
  q.b(notInside) = qn.b;
  q.c(notInside) = qn.c;
  q.d(notInside) = qn.d;
end

% ensure correct sign
if isempty(q_ref)
  changeSign = q.a < 0;
else
  changeSign = q.a .* q_ref.a + q.b .* q_ref.b + q.c .* q_ref.c + q.d .* q_ref.d < 0;
end
q.a(changeSign) = -q.a(changeSign);
q.b(changeSign) = -q.b(changeSign);
q.c(changeSign) = -q.c(changeSign);
q.d(changeSign) = -q.d(changeSign);

q = reshape(q,s);

% some testing code
% cs = crystalSymmetry('432')
% q = quaternion.rand(100,1);
% q_ref = quaternion.rand(100,1);
% q_proj = project2FR_ref(q,quaternion(cs),q_ref);
% hist(angle(q_proj,q_ref)/degree)

author: mtex-toolbox

repo: mtex


text:
% spline 举例
clear; clc;

xi=[27.7, 28, 29, 30]; % 插值节点
yi=[4.1, 4.3, 4.1, 3.0]; % 节点处的函数值
df0=3.0; dfn=-4.0; % 边界条件
pp=spline(xi,[df0, yi, dfn]);

xh=27.7:0.1:30; % 需要插值的点
yh=ppval(pp,xh); % 通过插值求得的近似值
plot(xi,yi,'r+',xh,yh,'o-b','LineWidth',1.5,'MarkerSize',12);

repo: NumericalAnalysis


text:
classdef Estimator < handle
    %ESTIMATOR  Rotation estimator base class
    %
    % It takes features of all images, pairwise matches between all images
    % and estimates rotations of all cameras.
    %
    % Note: The coordinate system origin is implementation-dependent, but you
    % can always normalize the rotations in respect to the first camera, for
    % instance.
    %
    % See also: cv.Stitcher, cv.BundleAdjuster
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = Estimator(estimatorType, varargin)
            %ESTIMATOR  Constructor
            %
            %     obj = cv.Estimator(estimatorType)
            %     obj = cv.Estimator(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __estimatorType__ Estimator type. One of:
            %   * __HomographyBasedEstimator__ Homography based rotation
            %     estimator.
            %   * __AffineBasedEstimator__ Affine transformation based
            %     estimator. This estimator uses pairwise transformations
            %     estimated by matcher to estimate final transformation for
            %     each camera.
            %
            % The following are options for the various algorithms:
            %
            % ### `HomographyBasedEstimator`
            % * __IsFocalsEstimated__ default false
            %
            % See also: cv.Estimator.estimate
            %
            this.id = Estimator_(0, 'new', estimatorType, varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.Estimator
            %
            if isempty(this.id), return; end
            Estimator_(this.id, 'delete');
        end

        function typename = typeid(this)
            %TYPEID  Name of the C++ type (RTTI)
            %
            %     typename = obj.typeid()
            %
            % ## Output
            % * __typename__ Name of C++ type
            %
            typename = Estimator_(this.id, 'typeid');
        end
    end

    %% Estimator
    methods
        function [cameras,success] = estimate(this, features, pairwise_matches)
            %ESTIMATE  Estimates camera parameters
            %
            %     cameras = obj.estimate(features, pairwise_matches)
            %     [cameras, success] = obj.estimate(...)
            %
            % ## Input
            % * __features__ Features of images. See cv.FeaturesFinder.
            % * **pairwise_matches** Pairwise matches of images. See
            %   cv.FeaturesMatcher.
            %
            % ## Output
            % * __cameras__ Estimated camera parameters. Structure that
            %   describes camera parameters with the following fields:
            %   * __aspect__ Aspect ratio.
            %   * __focal__ Focal length.
            %   * __ppx__ Principal point X.
            %   * __ppy__ Principal point Y.
            %   * __R__ 3x3 camera rotation matrix.
            %   * __t__ 3x1 camera translation vector.
            %   * __K__ 3x3 camera intrinsic parameters.
            % * __success__ True in case of success, false otherwise.
            %
            % See also: cv.Estimator.Estimator
            %
            [cameras,success] = Estimator_(this.id, 'estimate', features, pairwise_matches);
        end
    end

    %% Auto-calibration methods
    methods (Static)
        function [K,success] = calibrateRotatingCamera(Hs)
            %CALIBRATEROTATINGCAMERA  Calibrate rotating camera
            %
            %     K = cv.Estimator.calibrateRotatingCamera(Hs)
            %     [K,success] = cv.Estimator.calibrateRotatingCamera(Hs)
            %
            % ## Input
            % * __Hs__ Cell-array of 3x3 double matrices.
            %
            % ## Output
            % * __K__ 3x3 double matrix.
            % * __success__ True in case of success, false otherwise.
            %
            [K,success] = Estimator_(0, 'calibrateRotatingCamera', Hs);
        end

        function focals = estimateFocal(features, pairwise_matches)
            %ESTIMATEFOCAL  Estimates focal lengths for each given camera
            %
            %     focals = cv.Estimator.estimateFocal(features, pairwise_matches)
            %
            % ## Input
            % * __features__ Features of images.
            % * **pairwise_matches** Matches between all image pairs.
            %
            % ## Output
            % * __focals__ Estimated focal lengths for each camera, vector of
            %   doubles.
            %
            focals = Estimator_(0, 'estimateFocal', features, pairwise_matches);
        end

        function [f0, f1, f0_ok, f1_ok] = focalsFromHomography(H)
            %FOCALSFROMHOMOGRAPHY  Tries to estimate focal lengths from the given homography under the assumption that the camera undergoes rotations around its centre only
            %
            %     [f0, f1] = cv.Estimator.focalsFromHomography(H)
            %     [f0, f1, f0_ok, f1_ok] = cv.Estimator.focalsFromHomography(H)
            %
            % ## Input
            % * __H__ Homography, 3x3 double matrix.
            %
            % ## Output
            % * __f0__ Estimated focal length along X axis.
            % * __f1__ Estimated focal length along Y axis.
            % * **f0_ok** True, if `f0` was estimated successfully, false
            %   otherwise.
            % * **f1_ok** True, if `f1` was estimated successfully, false
            %   otherwise.
            %
            % ## References
            % > Heung-Yeung Shum and Richard Szeliski. "Construction of
            % > of Panoramic Image Mosaics with Global and Local Alignment".
            %
            [f0, f1, f0_ok, f1_ok] = Estimator_(0, 'focalsFromHomography', H);
        end
    end

end


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function ising_mom = ising_model_moments(Q)

n_vars = size(Q,1);

% Generate all binary vectors
bin_vals = dec2bin(0:2^n_vars-1)-'0';
bin_vals(bin_vals == 0) = -1;
n_vectors = size(bin_vals,1);

% Compute values of PDF
pdf_vals = zeros(n_vectors,1);
for i=1:n_vectors
	pdf_vals(i) = exp(bin_vals(i,:)*Q*bin_vals(i,:)');
end

% Compute normalizing constant
norm_const = sum(pdf_vals);

% Generate matrix to store moments
ising_mom = zeros(n_vars,n_vars);

% Compute second moment for each pair of values
for i=1:n_vars
	for j=1:n_vars
		bin_pair = bin_vals(:,i).*bin_vals(:,j);
		ising_mom(i,j) = sum(bin_pair.*pdf_vals)/norm_const;
	end
end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
function [randCat] = syn_randomize_catalog(mCatalog, bLon, bLat, bDepth, bTimes, nMagnitudes, fBValue, fMc, fInc)
    %  Randomizes a given catalog.
    %
    %  [randCat] = syn_randomize_catalog(mCatalog, bLon, bLat, bDepth,
    %                                                   bTimes, nMagnitudes, fBValue, fMc, fInc)
    %
    %
    % Input parameters:
    %   mCatalog          Catalog for randomizing
    %   bLon              Perturb longitudes (true) or leave longitudes unchanged (false)
    %   bLat              Perturb latitudes (true) or leave latitudes unchanged (false)
    %   bDepth            Perturb depths (true) or leave depths unchanged (false)
    %   bTimes            Perturb focal times (true) or leave focal times unchanged (false)
    %   nMagnitudes       Magnitude switch
    %                     1: Leave magnitudes unchanged
    %                     2: Generate new magnitudes according to parameters fBValue, fMc, fInc
    %                     3: Perturb magnitudes
    %   fBValue           b-value for new magnitudes
    %   fMc               magnitude of completeness for new magnitudes
    %   fInc              magnitude increment for new magnitudes
    %
    % Output parameters:
    %   randCat     Randomized catalog
    %
    % Danijel Schorlemmer
    % April 29, 2002

    randCat = copy(mCatalog);
    if isnumeric(nMagnitudes)
        if nMagnitudes == 3
            nMagnitudes = 'perturb';
        else
            nMagnitudes = 'create';
        end
    end
    % Permute longitudes
    if bLon
        randCat.Longitude = randCat.Longitude(randperm(randCat.Count));
    end

    % Permute latitudes
    if bLat
        randCat.Latitude = randCat.Latitude(randperm(randCat.Count));
    end

    % Permute depths
    if bDepth
        randCat.Depth = randCat.Depth(randperm(randCat.Count));
    end

    % Permute times
    if bTimes
        randCat.Date = randCat.Date(randperm(randCat.Count));
    end

    if nMagnitudes == "perturb"  %perturb magnitudes
        randCat.Magnitude = randCat.Magnitude(randperm(randCat.Count));
    elseif nMagnitudes == "create" % create new magnitudes
        randCat.Magnitude= syn_create_magnitudes(randCat.Count, fBValue, fMc, fInc);
    end
end

author: CelsoReyes


text:
function x = z1(y)
% linear normalize between 0 and 1
x = (y-min(y(:)))/(max(y(:))-min(y(:)))+eps;

% for multidimensional stuff, this normalizes each column to between 0 and
% 1 independent of other columns

% T=length(y);
% y=y';
% miny=min(y);
% x = (y-repmat(miny,T,1))./repmat(max(y)-min(y),T,1);

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function value = p20_exact ( dim_num )

%*****************************************************************************80
%
%% P20_EXACT returns the exact integral for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  aval = 0.0;
  aval = p20_r8 ( 'G', 'A', aval );

  bval = 0.0;
  bval = p20_r8 ( 'G', 'B', bval );

  p = 0.0;
  p = p20_r8 ( 'G', 'P', p );

  value = 0.0;
  exponent = dim_num + p;

  minus_one = -1.0;
  for i = 0 : dim_num
    minus_one = - minus_one;
    value = value + minus_one * r8_choose ( dim_num, i ) ...
      * ( ( dim_num - i ) * bval + i * aval )^exponent;
  end

  for i = 1 : dim_num
    value = value / ( p + i );
  end

  return
end


text:
function value = r8_pi

%*****************************************************************************80
%
%% R8_PI returns the value of pi.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the value of pi.
%
  value = 3.14159265358979323846264338327950288419716939937510;

  return
end


text:
function mne_write_inverse_sol_w(stem,inv,sol)
%
% function mne_write_inverse_sol_w(stem,inv,sol)
%
% Save static inverse solution data into stc files
%
% stem      - Stem for the w files
% inv       - The inverse operator structure (can be the forward operator as well)
% sol       - The solution matrix (number of locations)
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.2  2009/03/04 22:04:31  msh
%   Fixed typos in comments
%
%   Revision 1.1  2009/03/04 03:50:15  msh
%   Added mne_write_inverse_sol_w
%
%   Revision 1.2  2006/09/14 22:12:48  msh
%   Added output of the files written.
%
%   Revision 1.1  2006/05/05 03:50:40  msh
%   Added routines to compute L2-norm inverse solutions.
%   Added mne_write_inverse_sol_stc to write them in stc files
%   Several bug fixes in other files
%
%

me='MNE:mne_write_inverse_sol_w';
global FIFF;
if isempty(FIFF)
   FIFF = fiff_define_constants();
end

if nargin ~= 3
    error(me,'Incorrect number of arguments');
end

if length(sol) ~= inv.nsource
    error(me,'The solution data cannot correspond to these source spaces');
end

off = 0;
for k = 1:length(inv.src)
    off = off + inv.src(k).nuse;
    if (inv.src(k).id ~= FIFF.FIFFV_MNE_SURF_LEFT_HEMI && ...
            inv.src(k).id ~= FIFF.FIFFV_MNE_SURF_RIGHT_HEMI)
        error(me,'Source space hemispheres not properly assigned.');
    end
end
if off ~= inv.nsource
    error(me,'The source spaces are inconsistent with other inverse/forward operator data');
end
%
%   Write a separate stc file for each source space
%
off = 0;
for k = 1:length(inv.src)
    if (inv.src(k).id == FIFF.FIFFV_MNE_SURF_LEFT_HEMI)
        outname = sprintf('%s-lh.w',stem);
    elseif (inv.src(k).id == FIFF.FIFFV_MNE_SURF_RIGHT_HEMI)
        outname = sprintf('%s-rh.w',stem);
    end
    w.vertices = double(inv.src(k).vertno - 1);
    w.data     = sol(off+1:off+inv.src(k).nuse);
    mne_write_w_file(outname,w);
    off = off + inv.src(k).nuse;
    fprintf(1,'Wrote %s\n',outname);
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip