text:
function [newL2] = bestMap(L1,L2)
%bestmap: permute labels of L2 to match L1 as good as possible
%   [newL2] = bestMap(L1,L2);
%
%   version 2.0 --May/2007
%   version 1.0 --November/2003
%
%   Written by Deng Cai (dengcai AT gmail.com)

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

L1 = L1(:);
L2 = L2(:);
if size(L1) ~= size(L2)
    error('size(L1) must == size(L2)');
end

Label1 = unique(L1);
nClass1 = length(Label1);
Label2 = unique(L2);
nClass2 = length(Label2);

nClass = max(nClass1,nClass2);
G = zeros(nClass);
for i=1:nClass1
	for j=1:nClass2
		G(i,j) = length(find(L1 == Label1(i) & L2 == Label2(j)));
	end
end

[c,t] = hungarian(-G);
newL2 = zeros(size(L2));
for i=1:nClass2
    newL2(L2 == Label2(i)) = Label1(c(i));
end

author: jwyang


text:
function root_num = p03_root_num ( )

%*****************************************************************************80
%
%% P03_ROOT_NUM returns the number of known roots for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer ROOT_NUM, the number of known roots.
%
  root_num = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function nii_save(img,hdr,FileName)
% nii_save(img,hdr,FileName) save the matrix img in nifti
% img and hdr must have been loaded with nii_load

nii = nii_reset_orient(hdr,img);
nii_tool('save',nii,FileName)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% Gtomo2_wtfmex_test.m
% Test the Gtomo2_wtfmex object

if ~has_aspire
	return
end

% test new vewsion
if 1
	ig = image_geom('nx', 32, 'ny', 30, 'dx', 2);
	sg = sino_geom('par', 'nb', 40, 'na', 20, 'dr', 1.8);
	ig.mask = ig.circ(ig.fov/2) > 0;
	im(ig.mask)
	G = Gtomo2_wtmex(sg, ig, 'chat', 0, 'nthread', 2);
	Gtomo2_test(G, ig.mask) % put it through paces
	clear G
end

if ~isvar('G'),	printm 'setup Gtomo2_wtfmex_test'
	f.dir = test_dir;
	f.dsc = [f.dir 't.dsc'];
	f.wtf = strrep(f.dsc, 'dsc', 'wtf');

	os_run(sprintf('wt -chat 0 dsc 2 >! %s', f.dsc))
	os_run(sprintf('echo y | wt -chat 0 gen %s', f.dsc))

	G = Gtomo2_wtfmex(f.wtf);
	im(331, G.mask, 'mask')
prompt
end

if 1,	printm 'basic tests'
	x = G.mask;
	y = reshape(G * x(:), G.nb, G.na);
	im(332, y)

	y = ones(G.nb, G.na);
	x = reshape(G' * y(:), G.nx, G.ny);
	im(333, x, 'x'), cbar

	x2 = reshape((G.^2)' * y(:), G.nx, G.ny);
	im(337, x2, 'x2'), cbar

	t = reshape(G(:,2718), G.nb, G.na);
	im(334, t)

	t = reshape(G(:,[2718 2710]), [G.nb G.na 2]);
	im(335, t)

	t = reshape(G([202 210],:), [G.nx G.ny 2]);
	im(336, t)
prompt
end

if 1, printm 'adjoint test'

	os_run(sprintf('wt -chat 0 dsc 2 nx 12 ny 14 nb 16 na 10 >! %s', f.dsc))
	os_run(sprintf('echo y | wt -chat 0 gen %s', f.dsc))
	G = Gtomo2_wtfmex(f.wtf);
	test_adjoint(G);
	clear G
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [VMHId] = generateVMHMetAbbr(met, metabolite_structure_rBioNet,metab,rxnDB,customMetAbbrList)
% This function generates VMH ID's based on a metabolites based on
% predefined rules, as we would generally do it manually
%
% INPUT
% met                           Metabolite name
% metabolite_structure_rBioNet  To save time provide rBioNet either as
%                               1) string to mat file to load, e.g.,: /path/to/metabolite_structure_rBioNet.mat
%                               2) structure already in memory: metabolite_structure_rBioNet
%
% metab                         To save time provide rBioNet (as
%                               metab.mat file)
% rxnDB                         To save time provide rBioNet (as
%                               rxn.mat file)
% customMetaboliteList          List of metabolite abbr against which
%                               uniqueness should also be checked
%
% OUTPUT
% VMHId                         New VMH ID
%
% Ines Thiele, 09/2021

% load rBioNet
if 0
    mkdir('data/');
    websave('data/MetaboliteDatabase.txt','https://raw.githubusercontent.com/opencobra/COBRA.papers/master/2021_demeter/input/MetaboliteDatabase.txt');
    websave('data/ReactionDatabase.txt','https://raw.githubusercontent.com/opencobra/COBRA.papers/master/2021_demeter/input/ReactionDatabase.txt');
    createRBioNetDBFromVMHDB('rBioNetDBFolder','data/');
end
if ~exist('metab','var')
    load('data/metab.mat');
    load('data/rxn.mat');
else
    rxn = rxnDB;
end
% load extended rBioNet
if exist('metabolite_structure_rBioNet','var')
    if ~isstruct(metabolite_structure_rBioNet)
        load(metabolite_structure_rBioNet);
    end
else
    load met_strc_rBioNet;
end

% convert input into lower case
VMHId = lower(met);

% remove unneeded parts from the input
[VMHId] = removeJunk(VMHId);
% keep the original but curated input name
[met_VMHId] = VMHId;

% replace metabolite name with known abbreviations from rBioNet
for i = 1 : size(metab,1)
    name = lower(metab(i,2));
    abbr = lower(metab(i,1));
    % clean up a bit names
    [name] = removeJunk(name);
    nameList_rBioNet(i,1) = name;
    % check that abbr does not contain HC#
    if isempty(regexp(abbr{1},'m\d\d\d\d\d')) && isempty(regexp(abbr{1},'hc\d\d\d\d\d'))  && isempty(regexp(abbr{1},'ce\d\d\d\d')) ...
            && isempty(regexp(abbr{1},'cn\d\d\d\d'))  && isempty(regexp(abbr{1},'c\d\d\d\d\d'))&& ~strcmp(abbr{1},'m') && contains(VMHId,name)

        if strcmp(name,'sulfate')
            abbr = 's'; % don't use so4
        else
            new = ['(' abbr{1} ')'];
            VMHId = regexprep(VMHId,name,new);
        end
    end
end

F = fieldnames(metabolite_structure_rBioNet);
for i= 1: length(F)
    name = metabolite_structure_rBioNet.(F{i}).metNames;
    abbr = metabolite_structure_rBioNet.(F{i}).VMHId;
    [name] = removeJunk(name);
    nameList_rBioNet_strc{i,1} = name;
    if isempty(regexp(abbr,'m\d\d\d\d\d')) && isempty(regexp(abbr,'hc\d\d\d\d\d'))  && isempty(regexp(abbr,'ce\d\d\d\d')) ...
            && isempty(regexp(abbr,'cn\d\d\d\d'))  && isempty(regexp(abbr,'c\d\d\d\d\d'))&& ~strcmp(abbr,'m') && contains(VMHId,name)

        if strcmp(name,'sulfate')
            abbr = 's'; % don't use so4
        else

            if ischar(abbr)
                new = ['(' abbr ')'];
            else
                new = ['(' abbr{1} ')'];
            end
            VMHId = regexprep(VMHId,name,new);
        end
    end
end

% these rules are defined as I would have manually implemented them based
% on a metabolite name.
transl = {
    'glucuronide'   '-g-'
    'monophosphate' '-mp-'
    'diphosphate'   '-dp-'
    'triphosphate'   '-tp-'
    'phosphate' '-p-'
    'methyl'    '-m-'
    'uridine'   '-u-'
    'deoxy' '-d-'
    'dehydro'   '-d-'
    'fluoro'    '-f-'
    'uracil'    '-ura-'
    'acetyl'    '-ac-'
    'glutathione'   '-gth-'
    'quinone'   '-q-'
    'thio'  '-t-'
    'pyridil'   'py'
    'amino' '-am-'
    'amine' '-am-'
    'amide' '-ad-'
    'hydroxy'   '-h-'
    'alpha' '-a-'
    'beta'  '-b-'
    'gamma' '-g-'
    'formyl'    '-f-'
    'oxide' '-o-'
    'carboxyl' '-ca-'
    'carboxy' '-ca-'
    'carboxide' '-ca-'
    'oxi'   '-o-'
    'oxy'   '-o-'
    'riboside'  '-rib-'
    'phospho'   '-p-'
    'tauro' '-t-'
    'cheno' '-c-'
    'deaminated'    '-da-'
    'deaminate' '-da-'
    'deamine' '-da-'
    'acid'  '-a-'
    'nicotinic' '-ni-'
    'tadine' 'td'
    'di'    '-d-'
    'carbaldehyde'  '-ca-'
    'carbaldh'  '-ca-'
    'aldehyde'  '-al-'
    'oxo'   '-o-'
    'hydro' '-h-'
    'azole' '-az-'
    'tri'   '-t-'
    'warfarin'  '-warf-'
    'sulphate' '-s-'
    'glucuronate'   '-glc-'
    'benzyl'    '-b-'
    'tetra'    '-t-'
    'enoyl' '-e-'
    'trans' 't'
    'cis' 'c'
    'mono' 'm'
    };
transl = lower(transl);
for i = 1 : size(transl,1)
    VMHId = regexprep(VMHId,transl{i,1},transl{i,2});
end

% no shortening happened
VMHIdO = VMHId;

% choose the first 4 letters
VMHId = regexprep(VMHId,'-([a-z])([a-z])([a-z])([a-z])([a-z]+)-','-$1$2$3$4-');
VMHId = regexprep(VMHId,'-([a-z])([a-z])([a-z])([a-z])([a-z]+)$','-$1$2$3$4-');
VMHId = regexprep(VMHId,'^([a-z])([a-z])([a-z])([a-z])([a-z]+)-','-$1$2$3$4-');
VMHId = regexprep(VMHId,'([a-z])([a-z])([a-z])([a-z])([a-z]+)','-$1$2$3$4-');
% remove numbers between two letters
VMHIdN = VMHId;
% not of the form of DG(15:0/18:2(9Z,12Z)/0:0)
if length(find(strfind(VMHId,'/')))==0
    VMHId = regexprep(VMHId,'([a-z])(\d\d+)([a-z])','$1$3');
else
    VMHId = regexprep(VMHId,'\/','_');
end
VMHId = regexprep(VMHId,'\s','');
VMHId = regexprep(VMHId,'-','');
VMHId = regexprep(VMHId,'\(','');
VMHId = regexprep(VMHId,'\)','');
VMHId = regexprep(VMHId,'\:','');
VMHId = regexprep(VMHId,',','');
VMHId = regexprep(VMHId,'+','_');
VMHId = regexprep(VMHId,'{','');
VMHId = regexprep(VMHId,'}','_');
VMHId = regexprep(VMHId,'__','_');
VMHId = regexprep(VMHId,'^_','');
VMHId = regexprep(VMHId,'!','');
VMHId = regexprep(VMHId,'‐','');

% check that this abbr does not exist yet
[VMH_existance,rBioNet_existance] = checkAbbrExists({VMHId},metab,rxn,metabolite_structure_rBioNet);
% if the abbr already exists, try the version with the internal numbers
if ~isempty(find(contains(VMH_existance(:,3),'1'))) ||  ~isempty(find(contains(rBioNet_existance(:,3),'1')))
    % check whether it is the same metabolite based on name match (not
    % more)
    if isempty(find(ismember(nameList_rBioNet_strc,met_VMHId))) &&  isempty(find(ismember(nameList_rBioNet,met_VMHId)))
        % choose a 5 letter abbr
        VMHId = VMHIdO;
        % add 1 to the end
        VMHId = [VMHId '1'];
        [VMH_existance,rBioNet_existance] = checkAbbrExists({VMHId},metab,rxn,metabolite_structure_rBioNet);
        if ~isempty(find(contains(VMH_existance(:,3),'1'))) ||  ~isempty(find(contains(rBioNet_existance(:,3),'1')))
            VMHId = VMHIdO;
            % add 1 to the end
            VMHId = [VMHId '2'];
            [VMH_existance,rBioNet_existance] = checkAbbrExists({VMHId},metab,rxn,metabolite_structure_rBioNet);
        end

        VMHId = regexprep(VMHId,'\s','');
        VMHId = regexprep(VMHId,'-','');
        VMHId = regexprep(VMHId,'\(','');
        VMHId = regexprep(VMHId,'\)','');
        VMHId = regexprep(VMHId,'\:','');
        VMHId = regexprep(VMHId,',','');
        VMHId = regexprep(VMHId,'+','_');
        VMHId = regexprep(VMHId,'{','');
        VMHId = regexprep(VMHId,'}','_');
VMHId = regexprep(VMHId,'‐','');

    end
end

% now check also against the costum metabolite abbr list
if exist('customMetAbbrList','var')
    if ~isempty(customMetAbbrList) &&  length(strmatch(VMHId,customMetAbbrList,'exact'))>0
        % choose a 5 letter abbr
        VMHId = VMHIdO;
        % add 1 to the end
        VMHId = [VMHId '1'];
        num = 2;
        [VMH_existance,rBioNet_existance] = checkAbbrExists({VMHId},metab,rxn,metabolite_structure_rBioNet);
        while ~isempty(find(contains(VMH_existance(:,3),'1'))) ||  ~isempty(find(contains(rBioNet_existance(:,3),'1')))
            VMHId = VMHIdO;
            % add 1 to the end
            VMHId = [VMHId num2str(num)];
            num = num+1;
            [VMH_existance,rBioNet_existance] = checkAbbrExists({VMHId},metab,rxn,metabolite_structure_rBioNet);
        end

        VMHId = regexprep(VMHId,'\s','');
        VMHId = regexprep(VMHId,'-','');
        VMHId = regexprep(VMHId,'\(','');
        VMHId = regexprep(VMHId,'\)','');
        VMHId = regexprep(VMHId,'\:','');
        VMHId = regexprep(VMHId,',','');
        VMHId = regexprep(VMHId,'+','_');
        VMHId = regexprep(VMHId,'{','');
        VMHId = regexprep(VMHId,'}','_');
        VMHId = regexprep(VMHId,'/','_');
        VMHId = regexprep(VMHId,'‐','');

        % Check here again if it is in customMetAbbrList (needs first
        % converted with regexprep)
         while ~isempty(find(ismember(customMetAbbrList,VMHId))>0)
            VMHId = [VMHId '1'];
         end
    end
end

function [name] = removeJunk(name)
% remove parts of the metabolite name (in input as well as in rBioNet) that
% is not informative but makes mapping more challenging

name = regexprep(name,'conjugate','');
name = regexprep(name,'acceptor','');
name = regexprep(name,'donor','');
name = regexprep(name,'intermediate','');
name = regexprep(name,'\(1\-\)','');
name = regexprep(name,'\(1\+\)','');
name = regexprep(name,'\(2\-\)','');
name = regexprep(name,'\(2\+\)','');
name = regexprep(name,'\(3\-\)','');
name = regexprep(name,'\(3\+\)','');
name = regexprep(name,'\(4\-\)','');
name = regexprep(name,'\(4\+\)','');
name = regexprep(name,'\(1\)','');
name = regexprep(name,'\(2\)','');
name = regexprep(name,'^ ','');
name = regexprep(name,' $','');
name = regexprep(name,'*','');
name = regexprep(name,'''','');
name = regexprep(name,'β','b');
name = regexprep(name,'α','a');
name = regexprep(name,'Δ','d');
name = regexprep(name,'δ','d');
name = regexprep(name,'γ','g');
name = regexprep(name,'ω','o');
name = regexprep(name,'\.','_');
% these edits are specific to the metabolon input file
name = regexprep(name,'\[1\]','');
name = regexprep(name,'\[2\]','');
name = regexprep(name,'′','');
name = regexprep(name,'\[cis or trans\]','');
% remove any square brackets but retain text
name = regexprep(name,'\[','');
name = regexprep(name,'\]','');
name = regexprep(name,';','_');
name = regexprep(name,'!','');
name = regexprep(name,'‐','');
name = regexprep(name,'â','a');

author: opencobra

repo: cobratoolbox


text:
function [R,t,E] = psth(data,sig,plt,T,err,t)
%  function to plot trial averaged rate smoothed by
%  a Gaussian kernel - visual check on stationarity
%  Usage: [R,t,E] = psth(data,sig,plt,T,err,t)
%
%  Inputs:
% Note that all times have to be consistent. If data
% is in seconds, so must be sig and t. If data is in
% samples, so must sig and t. The default is seconds.
% data            structural array of spike times
% sig             std dev of Gaussian (default 50ms)
%                 (minus indicates adaptive with
%                  approx width equal to mod sig)
% plt = 'n'|'r' etc      (default 'r')
% T is the time interval (default all)
% err - 0 = none
%       1 = Poisson
%       2 = Bootstrap over trials (default)
% (both are based on 2* std err rather than 95%)
% t   = times to evaluate psth at
%
% The adaptive estimate works by first estimating
% the psth with a fixed kernel width (-sig) it
% then alters the kernel width so that the number
% of spikes falling under the kernel is the same on
% average but is time dependent.  Reagions rich
% in data therefore have their kernel width reduced
%
% Outputs:
%
% R = rate
% t = times
% E = errors (standard error)

if nargin <1 ; error('I need data!');end
[data]=padNaN(data); % create a zero padded data matrix from input structural array
sz=size(data);
% transposes data so that the longer dimension is assumed to be time
% Procedure used to bring input data into form compatible with that required
% for Murray's routine
if sz(1)>sz(2); data=data';end;
if nargin < 2; sig = 0.05;end
if nargin < 3; plt = 'r';end
if nargin < 5; err = 2; end

if isempty(sig); sig = 0.05;end
if isempty(plt); plt = 'r';end
if isempty(err); err = 2; end

adapt = 0;
if sig < 0; adapt = 1; sig = -sig; end

%  to avoid end effects increase time interval by the width
%  of the kernel (otherwise rate is too low near edges)

if nargin < 4;
  T(1) = min(data(:,1));
  T(2) = max(max(data));
else
  T(1) = T(1)-4*sig;
  T(2) = T(2)+4*sig;
end

% calculate NT and ND and filter out times of interest

NT = length(data(:,1));
if NT < 4 && err == 2
  disp('Using Poisson errorbars as number of trials is too small for bootstrap')
  err = 1;
end

m = 1;
D = zeros(size(data));
ND=zeros(1,NT);
for n=1:NT
  indx = find(~isnan(data(n,:)) & data(n,:)>=T(1) & data(n,:)<=T(2));
  ND(n) = length(indx);
  D(n,1:ND(n)) = data(n,indx);
  m = m + ND(n);
end
N_tot = m;
N_max = max(ND);
D = D(:,1:N_max);

% if the kernel density is such that there are on average
% one or less spikes under the kernel then the units are probably wrong

L = N_tot/(NT*(T(2)-T(1)));
if 2*L*NT*sig < 1 || L < 0.1
  disp('Spikes very low density: are the units right? is the kernel width sensible?')
  disp(['Total events: ' num2str(fix(100*N_tot)/100) ' sig: ' ...
        num2str(fix(1000*sig)) 'ms T: ' num2str(fix(100*T)/100) ' events/sig: ' ...
        num2str(fix(100*N_tot*sig/(T(2)-T(1)))/100)])
end

%    Smear each spike out
%    std dev is sqrt(rate*(integral over kernal^2)/trials)
%    for Gaussian integral over Kernal^2 is 1/(2*sig*srqt(pi))

if nargin < 6
  N_pts =  fix(5*(T(2)-T(1))/sig);
  t = linspace(T(1),T(2),N_pts);
else
  N_pts = length(t);
end

RR = zeros(NT,N_pts);
f = 1/(2*sig^2);
for n=1:NT
  for m=1:ND(n)
    RR(n,:) = RR(n,:) + exp(-f*(t-D(n,m)).^2);
  end
end
RR = RR*(1/sqrt(2*pi*sig^2));
if NT > 1; R = mean(RR); else R = RR;end

if err == 1
  E = sqrt(R/(2*NT*sig*sqrt(pi)));
elseif err == 2
  Nboot = 10;
  mE = 0;
  sE = 0;
  for b=1:Nboot
    indx = floor(NT*rand(1,NT)) + 1;
    mtmp = mean(RR(indx,:));
    mE = mE + mtmp;
    sE = sE + mtmp.^2;
  end
  E = sqrt((sE/Nboot - mE.^2/Nboot^2));
end

% if adaptive warp sig so that on average the number of spikes
% under the kernel is the same but regions where there is
% more data have a smaller kernel

if adapt
  sigt = mean(R)*sig./R;
  RR = zeros(NT,N_pts);
  f = 1./(2*sigt.^2);
  for n=1:NT
    for m=1:ND(n)
      RR(n,:) = RR(n,:) + exp(-f.*(t-D(n,m)).^2);
    end
    RR(n,:) = RR(n,:).*(1./sqrt(2*pi*sigt.^2));
  end
  if NT > 1; R = mean(RR); else R = RR;end

  if err == 1
    E = sqrt(R./(2*NT*sigt*sqrt(pi)));
  elseif err == 2
    Nboot = 10;
    mE = 0;
    sE = 0;
    for b=1:Nboot
      indx = floor(NT*rand(1,NT)) + 1;
      mtmp = mean(RR(indx,:));
      mE = mE + mtmp;
      sE = sE + mtmp.^2;
    end
    E = sqrt((sE/Nboot - mE.^2/Nboot^2));
  end
end

if plt == 'n';return;end
plot(t,R,plt)
hold on
if err > 0
  plot(t,R+2*E,'g')
  plot(t,R-2*E,'g')
end
%axis([T(1)+(4*sig) T(2)-(4*sig) 0 1.25*max(R)])
axis([T(1)+(4*sig) T(2)-(4*sig) 0 max(R+2*E)+10])
xlabel('time (s)')
ylabel('rate (Hz)')
title(['Trial averaged rate : Gaussian Kernel :'  ...
	    ' sigma = ' num2str(1000*sig) 'ms'])
hold off

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Tracker performance evaluation tool for LaSOT
% 07/19/2018 by Heng Fan

clc; clear; close all;

addpath('./utils/');
addpath('./sequence_evaluation_config/');

tmp_mat_path  = './tmp_mat/';          % path to save temporary results
path_anno     = './annos/';            % path to annotations
path_att      = './annos/att/';        % path to attribute
rp_all        = './tracking_results/'; % path to tracking results
save_fig_path = './res_fig/';          % path to result figures
save_fig_suf  = 'eps';                 % suffix of figures, 'png' or 'eps'

att_name      = {'Illumination Variation', 'Partial Occlusion', 'Deformation', ...
                 'Motion Blur', 'Camera Motion', 'Rotation', 'Background Clutter', ...
                 'Viewpoint Change', 'Scale Variation', 'Full Occlusion', 'Fast Motion', ...
                 'Out-of-View', 'Low Resolution', 'Aspect Ration Change'};
att_fig_name  = {'IV', 'POC', 'DEF', 'MB', 'CM', 'ROT', 'BC', ...
                 'VC', 'SV', 'FOC', 'FM', 'OV', 'LR', 'ARC'};

% 'all' --- evaluation with the whole benchmark
% 'test_set' --- evaluation with training subset
evaluation_dataset_type = 'all';

% use normalization or not
norm_dst = false;

trackers   = config_tracker();
sequences  = config_sequence(evaluation_dataset_type);
plot_style = config_plot_style();

num_seq = numel(sequences);
num_tracker = numel(trackers);

% load tracker info
name_tracker_all = cell(num_tracker, 1);
for i = 1:num_tracker
    name_tracker_all{i} = trackers{i}.name;
end

% load sequence info
name_seq_all = cell(num_seq, 1);
for i = 1:num_seq
    name_seq_all{i} = sequences{i};
    seq_att         = dlmread(fullfile(path_att, [sequences{i} '.txt']));
    if i == 1
        att_all = zeros(num_seq, numel(seq_att));
    end
    att_all(i, :) = seq_att;
end

% parameters for evaluation
metric_type_set = {'error', 'overlap'};
eval_type       = 'OPE';
% ranking_type    = 'AUC';
ranking_type    = 'threshold';   % change it to 'AUC' for success plots
rank_num        = 35;

threshold_set_error   = 0:50;
if norm_dst
    threshold_set_error = threshold_set_error / 100;
end
threshold_set_overlap = 0:0.05:1;

for i = 1:numel(metric_type_set)
    % error (for distance plots) or overlap (for success plots)
    metric_type = metric_type_set{i};
    switch metric_type
        case 'error'
            threshold_set = threshold_set_error;
            rank_idx      = 21;
            x_label_name  = 'Location error threshold';
            y_label_name  = 'Precision';
        case 'overlap'
            threshold_set = threshold_set_overlap;
            rank_idx      = 11;
            x_label_name  = 'Overlap threshold';
            y_label_name  = 'Success rate';
    end

%     if strcmp(metric_type, 'error') && strcmp(ranking_type, 'AUC') % for ranking_type = 'AUC'
    if strcmp(metric_type, 'overlap') && strcmp(ranking_type, 'threshold')  % for ranking_type = 'threshold'
        continue;
    end

    t_num = numel(threshold_set);

    % we only use OPE for evaluation
    plot_type = [metric_type '_' eval_type];

    switch metric_type
        case 'error'
            title_name = ['Precision plots of ' eval_type];
            if norm_dst
                title_name = ['Normalized ' title_name];
            end

            if strcmp(evaluation_dataset_type, 'all')
                title_name = [title_name ' on LaSOT'];
            else
                title_name = [title_name ' on LaSOT Testing Set'];
            end
        case 'overlap'
            title_name = ['Success plots of ' eval_type];

            if strcmp(evaluation_dataset_type, 'all')
                title_name = [title_name ' on LaSOT'];
            else
                title_name = [title_name ' on LaSOT Testing Set'];
            end
    end

    dataName = [tmp_mat_path 'aveSuccessRatePlot_' num2str(num_tracker) ...
                'alg_'  plot_type '.mat'];

    % evaluate tracker performance
    if ~exist(dataName, 'file')
        eval_tracker(sequences, trackers, eval_type, name_tracker_all, ...
                    tmp_mat_path, path_anno, rp_all, norm_dst);
    end

    % plot performance
    load(dataName);
    num_tracker = size(ave_success_rate_plot, 1);

    if rank_num > num_tracker || rank_num <0
        rank_num = num_tracker;
    end

    fig_name= [plot_type '_' ranking_type];
    idx_seq_set = 1:numel(sequences);

    % draw and save the overall performance plot
    plot_draw_save(num_tracker, plot_style, ave_success_rate_plot, ...
                   idx_seq_set, rank_num, ranking_type, rank_idx, ...
                   name_tracker_all, threshold_set, title_name, ...
                   x_label_name, y_label_name, fig_name, save_fig_path, ...
                   save_fig_suf);

    % draw and save the per-attribute performance plot
    att_trld = 0;
    att_num  = size(att_all, 2);
    for att_idx = 1:att_num    % for each attribute
        idx_seq_set = find(att_all(:, att_idx) > att_trld);
        if length(idx_seq_set) < 2
            continue;
        end
        disp([att_name{att_idx} ' ' num2str(length(idx_seq_set))]);

        fig_name   = [att_fig_name{att_idx} '_'  plot_type '_' ranking_type];
        title_name = ['Plots of ' eval_type ': ' att_name{att_idx} ' (' num2str(length(idx_seq_set)) ')'];

        switch metric_type
            case 'overlap'
                title_name = ['Success plots of ' eval_type ' - ' att_name{att_idx} ' (' num2str(length(idx_seq_set)) ')'];
            case 'error'
                title_name = ['Precision plots of ' eval_type ' - ' att_name{att_idx} ' (' num2str(length(idx_seq_set)) ')'];
                if norm_dst
                    title_name = ['Normalized ' title_name];
                end
        end

        plot_draw_save(num_tracker, plot_style, ave_success_rate_plot, ...
                   idx_seq_set, rank_num, ranking_type, rank_idx, ...
                   name_tracker_all, threshold_set, title_name, ...
                   x_label_name, y_label_name, fig_name, save_fig_path, ...
                   save_fig_suf);
    end
end

author: HengLan


text:
function rx = rxObliqueRx(rx);
%
% rx = rxObliqueRx(rx);
%
% Add an Rx positioned (for most alignment situations in KGS / Wandell
% Lab land) roughly perpindicular to the calcarine sulcus,
% for aligning oblique inplanes to a volume anatomy using mrRx.
%
% ras 01/06.
if ieNotDefined('rx')
    cfig = findobj('Tag','rxControlFig');
    rx = get(cfig,'UserData');
end

% % store the existing settings
% rx = rxStore(rx, 'Last Setting');

% ras 01/06:
% I find these rotations / flips cover the most common case,
% where the first slice is the most anterior, pretty well:
voxRatio = rx.rxVoxelSize ./ rx.volVoxelSize;
trans = rx.volDims/2 - rx.rxDims/2 + [60 0 0];
rot = deg2rad([90 63 90]);
scale = [1 -1 1];
newXform = affineBuild(trans, rot, scale, [0 0 0]);

% if there's an Rx Figure open, set it to sagittal view, which
% is easiest to prescribe off of:
if ishandle(rx.ui.volOri(1))
    selectButton(rx.ui.volOri, 3);
end

rx = rxSetXform(rx, newXform, 1);
rx = rxStore(rx, 'Oblique Rx');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
clc;

close all;
clear all;

disp('Running all the experiments together')

%featuriziation
% comment out when done with featurizer
promptresponse = input('would you like to featurize? [Y/N enter for N]', 's');
if (promptresponse=='y' || promptresponse=='Y')

    [num,txt,raw] = xlsread('data\final104.xls');
    descriptions = raw(2:size(raw,1),2);
    style_ratings = num(1:size(num,1),1);
    comfort_ratings = num(1:size(num,1),4);
    overal_ratings = num(1:size(num,1),5);
    [featurs, headers] = featurize_bigram(descriptions,20,0,1);
    csvwrite('data\forWeka_featuresonly.csv', featurs);
    csvwrite('data\headers.csv', headers);
end

ccaregression_generalized
ccaregression_generalized_adjusted
linearRegression_generalized
naivbayes_generalized
regresstree_generalized
%svm_generalized

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [g1, g2] = rbfinfwhiteXwhiteKernGradient(rbfKern, whiteKern, t1, varargin)

% RBFINFWHITEXWHITEKERNGRADIENT Compute gradient between the RBF-WHITE kernel
% (with integration limits between minus infinity and infinity) and the
% WHITE kernel.
% FORMAT
% DESC computes the gradient of an objective function with respect to cross
% kernel terms between RBF-WHITE and WHITE kernels for the multiple output
% kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE
% kernel.
% ARG whiteKern : the kernel structure associated with the WHITE
% kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of RBF-WHITE kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of WHITE kernel.
%
% FORMAT
% DESC computes the gradient of an objective function with respect
% to cross kernel terms between RBF-WHITE and WHITE kernels for
% the multiple output kernel.
% ARG rbfKern : the kernel structure associated with the RBF-WHITE
% kernel.
% ARG whiteKern : the kernel structure associated with the WHITE
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of objective function with respect to kernel
% parameters of RBF-WHITE kernel.
% RETURN g2 : gradient of objective function with respect to kernel
% parameters of WHITE kernel.
%
% SEEALSO : multiKernParamInit, multiKernCompute, rbfinfwhiteKernParamInit,
% whiteKernParamInit
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 5
    t2 = t1;
else
    t2 = varargin{1};
end
covGrad = varargin{end};

if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if rbfKern.variance ~= whiteKern.variance
  error('Kernels cannot be cross combined if they have different variances.')
end

g1 = zeros(1, 2);
g2 = 0; % The only parameter of the WHITE kernel (its variance) is already
        % accounted for in g1

% Parameters required for further computations
variance = rbfKern.variance;
invWidth = rbfKern.inverseWidth;

T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
deltaT = T1 - T2;

K = exp(-0.5*invWidth*(deltaT.^2));

% Gradient w.r.t. the inverse width
g1(1) = (0.5*variance/sqrt(2*pi)) ...
    * sum(sum((1/sqrt(invWidth)-sqrt(invWidth)*(deltaT.^2)) .* K .* covGrad));

% Gradient w.r.t. sigma_r^2
g1(2) = sqrt(invWidth/(2*pi)) * sum(sum(K .* covGrad));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = getpid(varargin)
% VL_GETPID Get MATLAB process ID
%   VL_GETPID() returns MATLAB process ID.
%
%   Algorithm::
%     VL_GETPID() returns the result of the getpid() system call on
%     UNIX-like operating systems and of GetCurrentProcessID() on
%     Windows.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_getpid(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [res,errorL2]=FBP(proj,geo,angles,varargin)
%TODO docs FBP
%
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%                     and
%                     https://www.mathworks.com/matlabcentral/fileexchange/view_license?file_info_id=35548
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Kyungsang Kim, modified by Ander Biguri
%--------------------------------------------------------------------------

% Assertion exploiting lazy evaluation
if isfield(geo,'mode') && ~strcmpi(geo.mode,'parallel')
     assert(false,'Only use FBP for parallel beam CT')
end
geo=checkGeo(geo,angles);

[filter,parker]=parse_inputs(proj,geo,angles,varargin);
geo.filter=filter;

%Input is data,geosize,angles

if size(geo.offDetector,2)==1
    offset=repmat(geo.offDetector,[1 size(angles,2)]);
else
    offset=geo.offDetector;
end

%% Weight
%proj=data
proj=permute(proj,[2 1 3]);

%% filter
proj_filt = filtering(proj,geo,angles,parker); % Not sure if offsets are good in here
%RMFIELD Remove fields from a structure array.
geo=rmfield(geo,'filter');
%% backproject

res=Atb((proj_filt),geo,angles)*geo.DSO(1)/geo.DSD(1);

if nargout>1
     error=proj-Ax(res,geo,angles);
     errorL2=norm(error(:));
end

end

function [filter, parker]=parse_inputs(proj,geo,alpha,argin)
opts=     {'filter','parker'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:FBP:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option isnot default, then extranc value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
         if isempty(ind)
            error('CBCT:FDK:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end

    switch opt
        % % % % % % % Verbose
        case 'parker'
            if default
                parker=0;
            else
                parker=val;
            end

        case 'filter'
            if default
                filter='ram-lak';
            else
                if  ~ischar( val)
                    error('CBCT:FDK:InvalidInput','Invalid filter')
                end
                filter=val;
            end

        otherwise
            error('CBCT:FDK:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in FAK()']);
    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function outfilenames = SaveFigLotsOfWays(hfig,basename,formats)

if exist('savefig1','var'),
  savefigfun = @savefig1;
else
  savefigfun = @savefig;
end

if nargin < 3,
  formats = {'pdf','png','svg','fig'};
end
[path,basename] = myfileparts(basename);

outfilenames = {};
if ismember('svg',formats),
  filename = [basename,'.svg'];
  plot2svg(fullfile(path,filename),hfig,'png');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'svg'});

if ismember('ai',formats),
  filename = [basename,'.ai'];
  saveas(hfig,fullfile(path,filename),'ai');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'ai'});

if ismember('fig',formats),
  filename = [basename,'.fig'];
  saveas(hfig,fullfile(path,filename),'fig');
  outfilenames{end+1} = fullfile(path,filename);
end
formats = setdiff(formats,{'fig'});

tmpbasename = regexprep(basename,'[^a-zA-Z_0-9]','_');
for i = 1:numel(formats),
  filename = [basename,'.',formats{i}];
  tmpfilename = [tmpbasename,'.',formats{i}];
  try
    savefig_pa(tmpfilename,hfig,formats{i});
    if ~isempty(path),
      movefile(tmpfilename,fullfile(path,filename));
    end
  catch ME,
    warning('Error while creating %s: %s',filename,getReport(ME));
    continue;
  end
  outfilenames{end+1} = fullfile(path,filename);
end

author: kristinbranson

repo: JAABA


text:
%% Radix-4 FFT Test Script
% This file runs three versions of a Radix-4 FFT written in MATLAB:
%
%   radix4FFT1_Float.m computes a radix-4 FFT for floating point data types
%
%   radix4FFT2_FixPt.m computes a fixed point radix-4 FFT (requires Fixed
%   Point Toolbox)
%
%   radix4FFT3_FixPtEML.m is an Embedded MATLAB version of the radix-4 FFT
%   that can be used in Simulink.  You can also generate C code for this
%   code (using Real Time Workshop).  This version can also be compiled
%   into a MEX'd executable that runs significantly faster than the fixed
%   point code.
%
% For a description of the radix-4 FFT algorithm see the following link to
% DSPDesignLine.com:
%
%       http://www.dspdesignline.com/showArticle.jhtml;jsessionid=5DBROAJQ3
%       SIWCQSNDLOSKHSCJUNN2JVN?articleID=204400676&queryText=freescale+fft
%
%% Set up a signal
% The test signal is the sum of 2 sinusoids plus some noise.  The test
% signal is a complex signal with real and imaginary components.
a=256;   % FFT Length
t=0:1/a:(a-1)/a;
s=(sin(2*pi*12*t)+.8*j*sin(2*pi*4.25*t)+.01*randn(size(t)))/2;
s = floor(s*16384)/16384;   %Quantized sum of 2 sines plus noise
subplot(211)
plot(t,real(s))
xlabel('time (sec)');ylabel('Amplitude');title('Test Signal - Real Component')
subplot(212)
plot(t,imag(s))
xlabel('time (sec)');ylabel('Amplitude');title('Test Signal - Imag Component')

%% Test Original Floating Point Code
% Run the Radix-4 algorithm as a floating point implementation.
%
% The function radix4FFT1_Float.m accepts real or complex floating point
% inputs and calculates the FFT.
S = radix4FFT1_Float(s);
S = bitrevorder(S);

% Calculate FFT using MATLAB function
Y = fft(s);

% Compare accuracy of Radix-4 FFT to MATLAB's FFT
errs = double(S) - Y;
Sig = sum(abs(Y).^2)/a;
Noise = sum(abs(errs).^2)/a;
SNR = 10*log10(Sig/Noise);
sprintf('SNR for 2 floating point FFT methods is: %6.2f dB', SNR)

figure
plotPYYf(s,a)

% The radix-4 FFT is nearly identical in accuracy to MATLAB's built in FFT
% computation.
%% Set Fixed Point Parameters
% The Fixed Point Toolbox can convert a test signal to a true fixed point
% data type with specified word lengths and fractional scaling. In this
% case the test signal is a 16 bit integer with 15 fractional bits.
wl = 16;
sfi=fi(s,1,wl,wl-1);    % Data is Q16,15
sfi.RoundMode = 'nearest';        % Fixed Point Rounding, etc.
sfi.OverflowMode = 'wrap';
sfi.ProductMode = 'KeepMSB';
sfi.ProductWordLength = wl*2;
sfi.SumMode = 'KeepMSB';
sfi.SumWordLength = wl*2;
plot(t,real(sfi))
xlabel('time (sec)');ylabel('Amplitude');title('Test Signal - Real Component')

%% Fixed Point Test
% Run the Radix-4 algorithm as a fixed point implementation.
%
% The function radix4FFT2_FixPt.m accepts fixed or floating point inputs
% and calculates the FFT
SFI = radix4FFT2_FixPt(sfi);
SFI = bitrevorder(SFI);

% Calculate FFT using MATLAB function
y = fft(s);

% Compare results
errs = double(SFI) - y;
Sig = sum(abs(y).^2)/a;
Noise = sum(abs(errs).^2)/a;
SNR = 10*log10(Sig/Noise);
sprintf('SNR for fixed vs floating point methods is: %6.2f dB', SNR)

plotPYYf(double(sfi),a)

% As expected, the accuracy of the fixed point FFT is reduced compared to
% the floating point computation.
%% Set Fixed Point Parameters
% The accuracy of the FFT for a shorter word length can be tested by
% changing the input signal to a 14 bit integer with 13 fractional bits.
wl = 14;
sfi=fi(s,1,wl,wl-1);    % Data is Q16,15
sfi.RoundMode = 'nearest';        % Fixed Point Rounding, etc.
sfi.OverflowMode = 'wrap';
sfi.ProductMode = 'KeepMSB';
sfi.ProductWordLength = wl*2;
sfi.SumMode = 'KeepMSB';
sfi.SumWordLength = wl*2;
plot(t,real(sfi))
xlabel('time (sec)');ylabel('Amplitude');title('Test Signal - Real Component')

%% Fixed Point Test
% Run the Radix-4 FFT algorithm with a lower precision data type.
%
SFI = radix4FFT2_FixPt(sfi);
SFI = bitrevorder(SFI);

% Calculate FFT using MATLAB function
y = fft(s);

% Compare results
errs = double(SFI) - y;
Sig = sum(abs(y).^2)/a;
Noise = sum(abs(errs).^2)/a;
SNR = 10*log10(Sig/Noise);
sprintf('SNR for fixed vs floating point methods is: %6.2f dB', SNR)

plotPYYf(double(sfi),a)

% In this case there is an approximate loss of 12 dB of accuracy compared
% to the 16 bit FFT computation

%% Use emlmex to compile code into executable
% The fixed point FFT code can be accelerated by changing the algorithm to
% an EML-compliant algorithm and compiling it with the emlmex command. The
% emlmex command will produce a mex'd version of the MATLAB algorithm.
emlmex -o radix4FFT3_MX -eg {sfi}  radix4FFT3_FixPtEML

%% Show speed of non-compiled code
tic;SFI = radix4FFT3_FixPtEML(sfi);toc
tic;SFI = radix4FFT3_FixPtEML(sfi);toc
tic;SFI = radix4FFT3_FixPtEML(sfi);toc

%% Show speed of compiled code
tic;SFI = radix4FFT3_MX(sfi);toc
tic;SFI = radix4FFT3_MX(sfi);toc
tic;SFI = radix4FFT3_MX(sfi);toc

% The MEX'd version of the FFT code runs over 600 times faster than the
% fixed point MATLAB algorithm.
%
%% Generate C source code for radix-4 FFT
% EML compliant algorithms can be used to generate C source code using Real
% Time Workshop.
rtwcfg = emlcoder.RTWConfig

emlc -v -s rtwcfg -eg {sfi} -o radix4FFT_C radix4FFT3_FixPtEML

% The source code can be inspected in the emcprj directory

%% Use the EML compliant FFT code in Simulink
% The attached Simulink model hiperman_4_2007b_FixPt.mdl uses the radix-4
% FFT algorithm in the OFDM_RX block (requires Simulink and Communications
% Blockset).

author: Sable

repo: mcbench-benchmarks


text:
function nProtocols = db_import(BrainstormDbDir)
% DB_IMPORT: Import a full database from a folder on the hard-drive
%
% USAGE:  db_import(BrainstormDbDir)
%         db_import()                : Ask for the database directory

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

% Authors: Francois Tadel, 2010

% No directory specified
nProtocols = 0;
if (nargin < 1) || isempty(BrainstormDbDir)
    % Default directory
    if ispc
        defDir = 'C:\';
    else
        defDir = bst_get('UserDir');
    end
    % Open 'Select directory' dialog
    BrainstormDbDir = bst_uigetdir(defDir, 'Please select database directory.');
    if isempty(BrainstormDbDir) || ~ischar(BrainstormDbDir)
        return
    end
end

% Get all the subfolders
files = dir(BrainstormDbDir);
dirList = files([files.isdir] == 1);
if isempty(dirList)
    return
end
iProtocols = [];

% Loop on each subfolder
for i = 1:length(dirList)
    % Ignore folders that start with a "."
    if (dirList(i).name(1) == '.')
        continue;
    end
    % Protocol folder
    ProtocolDir = bst_fullfile(BrainstormDbDir, dirList(i).name);
    % Try to get a subject folder and a study folder
    [subjectFile, iDepthSubject] = file_find(ProtocolDir, 'brainstormsubject*.mat', 3);
    [studyFile, iDepthStudy]     = file_find(ProtocolDir, 'brainstormstudy*.mat',   3);
    % If not both files are found, it's not a brainstorm protocol
    if isempty(subjectFile) || isempty(studyFile) || (iDepthSubject ~= 3) || (iDepthSubject ~= 3)
        continue;
    end

    % Initialize protocol structure
    sProtocol = db_template('ProtocolInfo');
    [tmp__,sProtocol.Comment] = bst_fileparts(ProtocolDir, 1);
    sProtocol.SUBJECTS = bst_fileparts(bst_fileparts(subjectFile), 1);
    sProtocol.STUDIES  = bst_fileparts(bst_fileparts(studyFile), 1);
    % Load protocol
    tmp = db_edit_protocol('load', sProtocol);
    if ~isempty(tmp) && ~isequal(tmp, -1)
        iProtocols(end+1) = tmp;
    end
end

% Select the last protocol
if ~isempty(iProtocols)
    gui_brainstorm('SetCurrentProtocol', iProtocols(1));
end

author: brainstorm-tools


text:
function phiFaceAverage = tvdMean2D(phi, u, FL)
% This function gets the value of the field variable phi defined
% over the MeshStructure and calculates the TVD average on
% the cell faces, based on the direction of the velocity vector for a uniform mesh.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Written by Ali A. Eftekhari
% See the license file

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;
uy = u.yvalue;

% check the size of the variable and the mesh dimension
Nx = u.domain.dims(1);
Ny = u.domain.dims(2);
dx=repmat(0.5*(u.domain.cellsize.x(1:end-1)+u.domain.cellsize.x(2:end)), 1, Ny);
dy=repmat(0.5*(u.domain.cellsize.y(1:end-1)+u.domain.cellsize.y(2:end))', Nx, 1);

%
phiX_p = zeros(Nx+1,Ny);
phiX_m = zeros(Nx+1,Ny);
phiY_p = zeros(Nx,Ny+1);
phiY_m = zeros(Nx,Ny+1);

% calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
% x direction
dphiX_p = (phi.value(2:Nx+2, 2:Ny+1)-phi.value(1:Nx+1, 2:Ny+1))./dx;
rX_p = dphiX_p(1:end-1,:)./fsign(dphiX_p(2:end,:));
phiX_p(2:Nx+1,:) = phi.value(2:Nx+1, 2:Ny+1)+0.5*FL(rX_p).* ...
    (phi.value(3:Nx+2,2:Ny+1)-phi.value(2:Nx+1, 2:Ny+1));
phiX_p(1, :) = (phi.value(1, 2:Ny+1)+phi.value(2, 2:Ny+1))/2; % left boundary
% y direction
dphiY_p = (phi.value(2:Nx+1, 2:Ny+2)-phi.value(2:Nx+1, 1:Ny+1))./dy;
rY_p = dphiY_p(:,1:end-1)./fsign(dphiY_p(:,2:end));
phiY_p(:,2:Ny+1) = phi.value(2:Nx+1, 2:Ny+1)+0.5*FL(rY_p).* ...
    (phi.value(2:Nx+1,3:Ny+2)-phi.value(2:Nx+1, 2:Ny+1));
phiY_p(:,1) = (phi.value(2:Nx+1,1)+phi.value(2:Nx+1,2))/2; % Bottom boundary

% calculate the upstream to downstream gradient ratios for u<0 (- ratio)
% x direction
rX_m = dphiX_p(2:end,:)./fsign(dphiX_p(1:end-1,:));
phiX_m(1:Nx,:) = phi.value(2:Nx+1, 2:Ny+1)+0.5*FL(rX_m).* ...
    (phi.value(1:Nx, 2:Ny+1)-phi.value(2:Nx+1, 2:Ny+1));
phiX_m(Nx+1,:) = (phi.value(end, 2:Ny+1)+phi.value(end-1, 2:Ny+1))/2; % right boundary
% y direction
rY_m = dphiY_p(:,2:end)./fsign(dphiY_p(:,1:end-1));
phiY_m(:,1:Ny) = phi.value(2:Nx+1, 2:Ny+1)+0.5*FL(rY_m).* ...
    (phi.value(2:Nx+1, 1:Ny)-phi.value(2:Nx+1, 2:Ny+1));
phiY_m(:, Ny+1) = (phi.value(2:Nx+1, end)+phi.value(2:Nx+1, end-1))/2; % top boundary

% calculate the average value
xvalue = (ux>0).*phiX_p+ ...
                        (ux<0).*phiX_m+ ...
                        0.5*(ux==0).*(phi.value(1:Nx+1,2:Ny+1)+phi.value(2:Nx+2,2:Ny+1));
yvalue = (uy>0).*phiY_p+ ...
                        (uy<0).*phiY_m+ ...
                        0.5*(uy==0).*(phi.value(2:Nx+1,1:Ny+1)+phi.value(2:Nx+1,2:Ny+2));
phiFaceAverage=FaceVariable(phi.domain, xvalue, yvalue, []);
end

function phi_out = fsign(phi_in)
% This function checks the value of phi_in and assigns an eps value to the
% elements that are less than or equal to zero, while keeping the signs of
% the nonzero elements
    phi_out = (abs(phi_in)>=eps).*phi_in+eps*(phi_in==0)+eps*(abs(phi_in)<eps).*sign(phi_in);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function Agal=step_mg_diff_setup(x,y)
%mg_diff_setup_step   GMG diffusion problem on step domain
%   Agal=mg_diff_setup_step(x,y)
%   input
%          x       x coordinate vector for coarse grid
%          y       y coordinate vector for coarse grid
%   output
%          Agal    discrete diffusion operator
%
%   IFISS function: AR; 20 January, 2003.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
n=length(y)-1; np=n/2; nq=n/4;
nvtx= np*(np*1) + (5*np+1)*(n+1);
% negative x-values
xneg=x(1:n/2);yneg=y(1:n/2);
xpos=x(n/2+1:end);ypos=y(n/2+1:end);
[Xneg,Ypos]=meshgrid(xneg,ypos);
xx=reshape(Xneg',np*(np+1),1);
yy=reshape(Ypos',np*(np+1),1);
xyleft=[xx(:),yy(:)];
% assembly process
kx = 1;
ky = 1;
mel=0;
% loop over 2x2 macroelements
for j=1:nq
   for i=1:nq
      mref=np*(ky-1)+kx;
      pref=nq*(j-1)+i;
      mel=mel+1;
      nvv(1) = mref;
      nvv(2) = mref+2;
      nvv(3) = mref+2*np+2;
      nvv(4) = mref+2*np;
      nvv(5) = mref+1;
      nvv(6) = mref+np+2;
      nvv(7) = mref+2*np+1;
      nvv(8)=  mref+np;
      nvv(9)=  mref+np+1;
      npp(1) = pref;
      npp(2) = pref+1;
      npp(3) = pref+nq+1;
      npp(4) = pref+nq;
      mv(mel,1:9)=nvv(1:9);
      mp(mel,1:4)=npp(1:4);
      kx = kx + 2;
   end
   ky = ky + 2;
   kx = 1;
end
%
% correction along the internal boundary
mref=2*np*(3*np+1)+1;
pref=2*nq*(3*nq+1)+1;
for mel=nq:nq:nq*nq;
   nvv=mv(mel,:);
   npp=mp(mel,:);
   nvv(2) = mref;
   nvv(3) = mref+10*np+2;
   nvv(6) = mref+5*np+1;
   npp(2) = pref;
   npp(3) = pref+5*nq+1;
   mv(mel,1:9)=nvv(1:9);
   mp(mel,1:4)=npp(1:4);
   mref=mref+10*np+2;
   pref=pref+5*nq+1;
end
%
% positive x_values
[Xpos,Y]=meshgrid(xpos,y);
xx=reshape(Xpos',(5*np+1)*(n+1),1);
yy=reshape(Y',(5*np+1)*(n+1),1);
xyright=[xx(:),yy(:)];
xy=[xyleft;xyright];
%
kx = 1;
ky = 1;
mel=nq*nq;
for j=1:np
   for i=1:5*nq
      mref = (5*np+1)*(ky-1)+kx + np*(np+1);
      pref = (5*nq+1)*(j-1)+i +  nq*(nq+1);
      mel=mel+1;
      nvv(1) = mref;
      nvv(2) = mref+2;
      nvv(3) = mref+10*np+4;
      nvv(4) = mref+10*np+2;
      nvv(5) = mref+1;
      nvv(6) = mref+5*np+3;
      nvv(7) = mref+10*np+3;
      nvv(8)=  mref+5*np+1;
      nvv(9)=  mref+5*np+2;
      npp(1) = pref;
      npp(2) = pref+1;
      npp(3) = pref+5*nq+2;
      npp(4) = pref+5*nq+1;
      mv(mel,1:9)=nvv(1:9);
      mp(mel,1:4)=npp(1:4);
      kx = kx + 2;
   end
   ky = ky + 2;
   kx = 1;
end
%
% compute boundary vertices
% six boundary edges
k1=find( xy(:,1) <0  & xy(:,2)==0 );
e1=[]; for k=1:mel, if any(mv(k,5)==k1), e1=[e1,k]; end, end
ef1=ones(size(e1));
%
k2=find( xy(:,1)==0  & xy(:,2)<=0 );
e2=[]; for k=1:mel, if any(mv(k,8)==k2), e2=[e2,k]; end, end
ef2=4*ones(size(e2));
%
k3=find( xy(:,1) >0  & xy(:,2)==-1);
e3=[]; for k=1:mel, if any(mv(k,5)==k3), e3=[e3,k]; end, end
ef3=ones(size(e3));
%
k5=find( xy(:,2)==1 );
e5=[]; for k=1:mel, if any(mv(k,7)==k5), e5=[e5,k]; end, end
ef5=3*ones(size(e5));
%
k6=find( xy(:,1)==-1 & xy(:,2)<1   & xy(:,2) >0 );
e6=[]; for k=1:mel, if any(mv(k,8)==k6), e6=[e6,k]; end, end
ef6=4*ones(size(e6));
%
bound=sort([k1;k2;k3;k5;k6]);
mbound=[e1',ef1';e2',ef2';e3',ef3';;e5',ef5';e6',ef6'];
%
%% set up matrices for Q1 approximation
[ev,ebound]=mg_q1grid(x,y,xy,mv,bound,mbound);
[A,M,fdummy] = mg_q1diff(xy,ev);
%
% impose zero boundary conditions
[Agal] = mg_zerobc(A,xy,bound);
return


text:
function params = RightStancePreProcess(sys, t, x, controller, params) %#ok<INUSL>

    y = struct2array(sys.VirtualConstraints);

    for i=1:numel(y)
        if strcmp(y(i).PhaseType,'TimeBased')
            params.(y(i).PhaseParamName) = t + params.(y(i).PhaseParamName);
        end
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function bc = specific_bc(xbd,ybd)
%forwardstep_bc  Forward facing step flow boundary condition
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies streamfunction associated with Forward step flow
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bc=0*xbd;
k=find(xbd==5); bc(k)=(ybd(k).^3 -3*ybd(k) -2)/6;
k=find(ybd==1); bc(k)=-2/3;
return


text:
function [U,Usteps] = conformalized_mean_curvature_flow(V,F,varargin)
  % CONFORMALIZED_MEAN_CURVATURE_FLOW Flow a surface according to "Can mean
  % curvature flow be made non-singular?" [Kazhdan et al. 2012]
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices into V
  %   Optional:
  %     'MaxIter' followed by maximum number of iterations {100}
  %     'MinDiff' followed by minimum difference between iterations {1e-13}
  %     'Conformalize' followed by whether to rebuild _just_ the mass
  %       matrix each step {true}
  %     'delta' followed by delta value, should roughly be in range
  %       [1e-13,1e13] {1}
  %     'LaplacianType' followed by 'cotangent' of 'uniform'.
  %     'V0' followed by #V by 3 mesh positions to treat as initial mesh (to
  %       build laplacian from)
  %     'RescaleOutput' followed by whether to scale output to match input
  %       (otherwise scaled to unit surface area and moved to origin for
  %       numerical robustness) {false}
  % Outputs:
  %   U  #V by dim list of new vertex positions
  %   Usteps  #V by dim by iterations list of vertex positions during flow
  %
  function L = laplacian(V,F)
    switch laplacian_type
    case 'cotangent'
      L = cotmatrix(V,F);
    case 'uniform'
      A = adjacency_matrix(F);
      L = A - diag(sparse(sum(A,2)));
    end
  end

  % default values
  delta = 1;
  max_iter = 100;
  laplacian_type = 'cotangent';
  until_self_intersection_free = false;
  V0 = V;
  rescale_output = false;
  min_diff = 1e-13;
  conformalize = true;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'MaxIter','Conformalize','delta','LaplacianType','V0', ...
      'RescaleOutput', 'UntilSelfIntersectionFree','MinDiff'}, ...
    {'max_iter','conformalize','delta','laplacian_type','V0', ...
      'rescale_output', 'until_self_intersection_free','min_diff'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  L = laplacian(V,F);

  switch size(F,2)
  case 3
    SF = F;
  case 4
    SF = boundary_faces(F);
  end

  if nargout > 1
    Usteps = zeros([size(V) max_iter]);
  end
  U = V;
  iter = 1;
  while true
    if nargout > 1
      Usteps(:,:,iter) = U;
    end
    if until_self_intersection_free
      [~,~,IF] = selfintersect(U,SF,'DetectOnly',true,'FirstOnly',true);
      if isempty(IF)
        break;
      end
    end
    U_prev = U;
    % 'full' seems slight more stable than 'barycentric' which is more stable
    % than 'voronoi'
    M = massmatrix(U,F,'barycentric');
    if ~conformalize
      L = laplacian(V,F);
    end
    U = (M-delta(min(iter,end))*L)\(M*U);
    area = sum(doublearea(U,SF)*0.5);
    c = sum(bsxfun(@times,0.5*doublearea(U,SF)/area,barycenter(U,SF)));
    U = bsxfun(@minus,U,c);
    U = U/sqrt(area);
    % Use difference from previous as stopping criterion
    % Better would be to look for convergence while factoring out M??bius
    % transformation.
    % Q: Stop when no change in angles?
    d = trace(((U-U_prev)'*M*(U-U_prev)).^2);
    if d < min_diff
      warning('converged...');
      break;
    end
    % Volume of unit area sphere: pi^-0.5/6
    %[c,vol] = centroid(U,F);
    %tsurf(F,U);
    %%[vol pi^-0.5/6]
    %title(sprintf('%g',sum(doublearea(U,F)*0.5)));
    %drawnow;
    if iter >= max_iter
      warning('Max iterations (%d) exceeded without convergence',max_iter);
      break;
    end
    iter = iter + 1;
  end

  if nargout > 1
    Usteps = Usteps(:,:,1:iter);
  end

  if rescale_output
    area = sum(doublearea(V,SF)*0.5);
    c = sum(bsxfun(@times,0.5*doublearea(V,SF)/area,barycenter(V,SF)));
    U = U*sqrt(area);
    U = bsxfun(@plus,U,c);
    if nargout > 1
      for iter = 2:size(Usteps,3)
        Usteps(:,:,iter) = Usteps(:,:,iter)*sqrt(area);
        Usteps(:,:,iter) = bsxfun(@plus,Usteps(:,:,iter),c);
      end
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
function MDP = DEM_demo_MDP_rule
% Demo of active inference for visual salience
%__________________________________________________________________________
%
% This routine simulates a crude form of consciousness using active
% inference and structure learning to vitiate ignorance or nescience. This
% entails learning the hyper parameters of causal structure generating
% outcomes and then using Bayesian model reduction (during sleep) to
% minimise complexity.
%
% We first set up an abstract problem in which an agent has to respond
% according to rules (identify the correct colour depending upon one of
% three rules that are specified by the colour of a cue in the centre of
% vision). If the rule is centre, the colour is always green; however, if
% the colour of the Centre cue is red, the correct colour is on the left
% (and on the right if the queue is blue). Simulations are provided when
% the agent knows the rules. This is then repeated in the absence
% (nescience) of any knowledge about the rules to see if the agent can
% learn causal structure through Bayesian belief updating of the likelihood
% array (A).
%
% We then consider the improvement in performance (in terms of variational
% free energy, its constituent parts and performance) following Bayesian
% model reduction of the likelihood model (heuristically, like slow wave
% sleep), followed by a restitution of posterior beliefs during fictive
% active inference (as in REM sleep). Finally, we address the communication
% of the implicit  structure learning to a conspecific or child to
% demonstrate the improvement under instruction.
%
% see also: DEM_demo_MDP_habits.m and spm_MPD_VB_X.m
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_MDP_rule.m 7679 2019-10-24 15:54:07Z spm $

% set up and preliminaries
%==========================================================================

% specify the (rule-based) generative process (environment)
%==========================================================================
rng('default')

% prior beliefs about initial states (in terms of counts_: D
%--------------------------------------------------------------------------
D{1} = [1 1 1]';        % rule:   {'left','centre','right'}
D{2} = [1 1 1]';        % what:   {'red','green','blue'}
D{3} = [0 0 0 1]';      % where:  {'left','centre','right','null'}
D{4} = [0 0 0 1]';      % report: {'red','green','blue','undecided'}

% probabilistic mapping from hidden states to outcomes: A
%--------------------------------------------------------------------------
Nf    = numel(D);
Ns    = zeros(1,Nf);
for f = 1:Nf
    Ns(f) = numel(D{f});
end
for f1 = 1:Ns(1)                     % rule
    for f2 = 1:Ns(2)                 % correct colour
        for f3 = 1:Ns(3)             % location of fixation
            for f4 = 1:Ns(3)         % decision

                % A{1} what: {'red','green','blue','null'}
                %==========================================================
                if f3 == 4, A{1}(4,f1,f2,f3,f4)  = 1;  end
                if f3 == 2, A{1}(f1,f1,f2,f3,f4) = 1;  end
                if f3 == 1
                    if f1 == 1
                        A{1}(f2,f1,f2,f3,f4) = 1;
                    else
                        A{1}(1:3,f1,f2,f3,f4) = 1/3;
                    end
                end
                if f3 == 3
                    if f1 == 3
                        A{1}(f2,f1,f2,f3,f4) = 1;
                    else
                        A{1}(1:3,f1,f2,f3,f4) = 1/3;
                    end
                end

                % A{2} where: {'left','centre','right','null'}
                %----------------------------------------------------------
                A{2}(f3,f1,f2,f3,f4) = 1;

                % A{3} feedback: {'null','right','wrong'}
                %----------------------------------------------------------
                if f4 == 4,
                    A{3}(1,f1,f2,f3,f4) = 1;                             % undecided
                else
                    if f1 == 2 && f4 == 2,  A{3}(2,f1,f2,f3,f4) = 1; end % right
                    if f1 == 2 && f4 ~= 2,  A{3}(3,f1,f2,f3,f4) = 1; end % wrong
                    if f1 == 1 && f4 == f2, A{3}(2,f1,f2,f3,f4) = 1; end % right
                    if f1 == 1 && f4 ~= f2, A{3}(3,f1,f2,f3,f4) = 1; end % wrong
                    if f1 == 3 && f4 == f2, A{3}(2,f1,f2,f3,f4) = 1; end % right
                    if f1 == 3 && f4 ~= f2, A{3}(3,f1,f2,f3,f4) = 1; end % wrong
                end
            end
        end
    end
end
Ng    = numel(A);
for g = 1:Ng
    No(g) = size(A{g},1);
    A{g}  = double(A{g});
end

% controlled transitions: B{f} for each factor
%--------------------------------------------------------------------------
for f = 1:Nf
    B{f} = eye(Ns(f));
end

% control states B(3): where {'left','centre','right','null'}
%--------------------------------------------------------------------------
for k = 1:Ns(3)
    B{3}(:,:,k) = 0;
    B{3}(k,:,k) = 1;
end

% control states B(4): report {'red','green','blue','undecided'}
%--------------------------------------------------------------------------
for k = 1:Ns(4)
    B{4}(:,:,k) = 0;
    B{4}(k,:,k) = 1;
end

% allowable policies (specified as the next action) U
%--------------------------------------------------------------------------
U(1,1,:)  = [1 1 1 4]';         % sample left
U(1,2,:)  = [1 1 2 4]';         % sample centre
U(1,3,:)  = [1 1 3 4]';         % sample right
U(1,4,:)  = [1 1 4 1]';         % return to centre and report red
U(1,5,:)  = [1 1 4 2]';         % return to centre and report green
U(1,6,:)  = [1 1 4 3]';         % return to centre and report blue

% priors: (utility) C; the agent expects to avoid mistakes
%--------------------------------------------------------------------------
for g = 1:Ng
    C{g}  = zeros(No(g),1);
end
% and expects itself to make a decision after the fifth observation
%--------------------------------------------------------------------------
C{3} = [ 0  0  0  0 -8 -8;
         0  0  0  0  0  0;
        -4 -4 -4 -4 -4 -4];

% MDP Structure
%--------------------------------------------------------------------------
mdp.T = size(C{3},2);           % number of moves
mdp.U = U;                      % allowable policies
mdp.A = A;                      % observation model
mdp.B = B;                      % transition probabilities
mdp.C = C;                      % preferred outcomes
mdp.D = D;                      % prior over initial states

mdp.Aname = {'what','where','feedback'};
mdp.Bname = {'rule','colour','where','decision'};

% illustrate a single trial
%==========================================================================
mdp  = spm_MDP_check(mdp);
MDP  = spm_MDP_VB_X(mdp);

% show belief updates (and behaviour)
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1'); clf
spm_MDP_VB_trial(MDP);

% illustrate phase-precession and responses
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf
spm_MDP_VB_LFP(MDP,[],1);

% illustrate behaviour
%--------------------------------------------------------------------------
subplot(3,2,3)
spm_MDP_rule_plot(MDP)

% illustrate a sequence of trials
%==========================================================================
clear MDP

% create structure array
%--------------------------------------------------------------------------
N     = 8;
for i = 1:N
    MDP(i) = mdp;
end

% Solve an example sequence under different initial states
%==========================================================================
MDP   = spm_MDP_VB_X(MDP);

% illustrate behavioural responses and neuronal correlates
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf
spm_MDP_VB_game(MDP);

[F,Fu] = spm_MDP_F(MDP);

subplot(6,1,4), plot(1:N,F),  xlabel('trial'), spm_axis tight, title('Free energy','Fontsize',16)
subplot(6,1,5), plot(1:N,Fu), xlabel('trial'), spm_axis tight, title('Confidence','Fontsize',16)
subplot(6,1,6), spm_MDP_plot_moves(MDP),       spm_axis tight
title('Action ','Fontsize',16), legend({'saccades','hits'})

% illustrate phase-amplitude (theta-gamma) coupling
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf
spm_MDP_VB_LFP(MDP);

% illustrate behaviour in more detail
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 5'); clf;
for i = 1:min(N,15)
    subplot(5,3,i), spm_MDP_rule_plot(MDP(i));
    axis square
end

% repeat with rule learning
%==========================================================================

% retain knowledge about the cue but remove knowledge about the rules
%--------------------------------------------------------------------------
for f1 = 1:Ns(1)                     % rule
    for f2 = 1:Ns(2)                 % correct colour
        for f3 = 1:Ns(3)             % location of fixation
            for f4 = 1:Ns(3)         % decision

                % A{1} what: {'red','green','blue','null'}
                %==========================================================
                if f3 == 4, a{1}(4,f1,f2,f3,f4)  = 128;  end
                if f3 == 2, a{1}(f1,f1,f2,f3,f4) = 128;  end
                if f3 == 1
                    a{1}(1:3,f1,f2,f3,f4) = 1;
                end
                if f3 == 3
                    a{1}(1:3,f1,f2,f3,f4) = 1;
                end
            end
        end
    end
end

a{2}   = A{2}*128;
a{3}   = A{3}*128;
mda    = mdp;
mda.a  = a;
mda.a0 = a;

% create structure array
%--------------------------------------------------------------------------
clear MDP
N     = 32;
for i = 1:N
    MDP(i) = mda;
end

% Solve - an example sequence
%==========================================================================
rng('default')
MDP  = spm_MDP_VB_X(MDP);

% show responses to a difficult (right) trial before and after learning
%--------------------------------------------------------------------------
for i = 1:N
    s(:,i) = MDP(i).s(:,1);
end
i = find(ismember(s(1:2,:)',s(1:2,1)','rows'));

spm_figure('GetWin','Figure 6 - before');
spm_MDP_VB_LFP(MDP(i(1)),[],2);
subplot(3,2,2), set(gca,'YLim',[-.1 1])
subplot(3,2,4), set(gca,'YLim',[-.1 .2])
subplot(3,2,3), spm_MDP_rule_plot(MDP(i(1)))
title(['trial ' num2str(i(1))],'FontSize',16)

spm_figure('GetWin','Figure 6 - after' );
spm_MDP_VB_LFP(MDP(i(end)),[],2);
subplot(3,2,2), set(gca,'YLim',[-.1 1]);
subplot(3,2,4), set(gca,'YLim',[-.1 .2])
subplot(3,2,3), spm_MDP_rule_plot(MDP(i(end)))
title(['trial ' num2str(i(end))],'FontSize',16)

% show learning effects in terms of underlying uncertainty
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 7 - before'); clf; spm_MDP_VB_trial(MDP(i(1)));
spm_figure('GetWin','Figure 7 - after' ); clf; spm_MDP_VB_trial(MDP(i(end)));

% illustrate phase-amplitude (theta-gamma) coupling
%--------------------------------------------------------------------------
n  = min(N,4);
spm_figure('GetWin','Figure 4a'); clf
spm_MDP_VB_LFP(MDP(1:n),[],2);
for i = 1:n
    subplot(4,n,i + 3*n), spm_MDP_rule_plot(MDP(i));
    axis square
end

% show trial-by-trial measures in terms of free energy terms
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 8'); clf;
spm_MDP_VB_game(MDP);

[F,Fu] = spm_MDP_F(MDP);

subplot(6,1,4), plot(1:N,F),  xlabel('trial'), spm_axis tight, title('Free energy','Fontsize',16)
subplot(6,1,5), plot(1:N,Fu), xlabel('trial'), spm_axis tight, title('Confidence','Fontsize',16)
subplot(6,1,6), spm_MDP_plot_moves(MDP),       spm_axis tight
title('Action ','Fontsize',16), legend({'saccades','hits'})

% illustrate Bayesian model reduction
%--------------------------------------------------------------------------
OPTIONS.g = 1;
OPTIONS.f = 2;
OPTIONS.T = 3;
OPTIONS.m = @(i,i1,i2,i3,i4) i == i2;

for n = 1:N
    sdp{n} = spm_MDP_VB_sleep(MDP(n),OPTIONS);
    cor(n) = corr(spm_vec(MDP(n).A{1} > 0),spm_vec(sdp{n}.a{1}) > 0);
end
[m,n] = max(cor);

spm_figure('GetWin','Figure 9'); clf; str = sprintf('Sleep (trial %d)',n);
for n = n
    subplot(3,2,1), spm_MDP_A_plot(MDP(n).A);
    subplot(3,2,2), spm_MDP_A_plot(MDP(n).a0), title('Before','Fontsize',16)
    subplot(3,2,3), spm_MDP_A_plot(MDP(n).a),  title('After', 'Fontsize',16)
    subplot(3,2,4), spm_MDP_A_plot(sdp{n}.a),  title(str,     'Fontsize',16)
end

% return % here for short demo

% Bayesian model reduction with dreaming
%--------------------------------------------------------------------------
OPTIONS.o = {MDP.o};
rdp       = spm_MDP_VB_sleep(MDP(n),OPTIONS);
sdp       = sdp{n};

% illustrate benefits of sleep (with and without dreaming)
%==========================================================================
RDP   = MDP(n:end);
RDP   = rmfield(RDP,{'o','u'});
for i = 1:length(RDP)
    RDP(i).s = RDP(i).s(:,1);
end
SDP    = RDP;
SDP(1).a = sdp.a;
RDP(1).a = rdp.a;
SDP    = spm_MDP_VB_X(SDP);
RDP    = spm_MDP_VB_X(RDP);

% plot performance and without sleep
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 10'); clf;

[F,Fu] = spm_MDP_F(MDP);

subplot(4,1,1),plot(1:N,F, ':'),  xlabel('trial'), spm_axis tight, hold on
subplot(4,1,2),plot(1:N,Fu,':'),  xlabel('trial'), spm_axis tight, hold on

% and after SWS sleep
%--------------------------------------------------------------------------
[F,Fu] = spm_MDP_F(SDP);

subplot(4,1,1),plot(n:N,F ,'-.'), xlabel('trial'), spm_axis tight, hold on
subplot(4,1,2),plot(n:N,Fu,'-.'), xlabel('trial'), spm_axis tight, hold on

% and after REM sleep
%--------------------------------------------------------------------------
[F,Fu] = spm_MDP_F(RDP);

subplot(4,1,1),plot(n:N,F ,'-'),  xlabel('trial'), spm_axis tight, hold on
title('Free energy (with sleep)','Fontsize',16)
subplot(4,1,2),plot(n:N,Fu,'-'),  xlabel('trial'), spm_axis tight, hold on
title('Confidence (with sleep)','Fontsize',16)

% and correct responses
%--------------------------------------------------------------------------
hit = @(MDP) any(MDP.o(3,:) == 2) & ~any(MDP.o(3,:) == 3);
k   = min(Fu);
hold on
for i = 1:numel(RDP)
    j = i + n - 1;
    if hit(SDP(i)), plot(j,k,'hr','MarkerSize',8); end
end
for i = 1:(n - 1)
    if hit(MDP(i)), plot(i,k,'hb','MarkerSize',8); end
end
hold off

% illustrate the benefits of instruction (communication of model)
%==========================================================================

% create instructed agent
%--------------------------------------------------------------------------
n     = 1:16;
NDP   = MDP(n);
NDP   = rmfield(NDP,{'o','u'});
for i = 1:length(NDP)
    NDP(i).s = NDP(i).s(:,1);
    NDP(i).a = sdp.a0;
end
NDP   = spm_MDP_VB_X(NDP);

% plot performance and without received wisdom
%--------------------------------------------------------------------------
[F,Fu] = spm_MDP_F(MDP(n));

subplot(4,2,5),plot(n,F ,':'), xlabel('trial'), spm_axis tight, hold on
subplot(4,2,6),plot(n,Fu,':'), xlabel('trial'), spm_axis tight, hold on

% and after instruction
%--------------------------------------------------------------------------
[F,Fu] = spm_MDP_F(NDP(n));
k      = min(Fu);

subplot(4,2,5),plot(n,F),       xlabel('trial'), spm_axis tight
title('Free energy (instruction)','Fontsize',16)
subplot(4,2,6),plot(n,Fu),      xlabel('trial'), spm_axis tight
title('Confidence (instruction)','Fontsize',16)

% and correct responses
%--------------------------------------------------------------------------
for i = n
    hold on
    if hit(NDP(i)); plot(i,k,'hr','MarkerSize',8), end
    if hit(MDP(i)); plot(i,k,'hb','MarkerSize',8), end
    hold off
end

% run multiple subjects
%==========================================================================
BMR.g   = 1;
BMR.f   = 2;
BMR.T   = 3;
BMR.m   = @(i,i1,i2,i3,i4) i == i2;
OPT.BMR = BMR;

N     = 32;
Ns    = 64;
for m = 1:Ns

    % create structure array and solve
    %----------------------------------------------------------------------
    clear MDP OPTIONS
    for i = 1:N
        MDP(i) = mda;
    end
    rng(m)
    MDP  = spm_MDP_VB_X(MDP);

    % free energy and confidence
    %----------------------------------------------------------------------
    [F,Fu]   = spm_MDP_F(MDP);
    Fm(:,m)  = F(:);
    Fum(:,m) = Fu(:);

    % find run of correct responses
    %----------------------------------------------------------------------
    for i = 1:N
        if hit(MDP(i)); h(i,m) = 1; else, h(i,m) = 0;  end
    end

    % repeat with BMR
    %----------------------------------------------------------------------
    RDP = MDP;
    RDP = rmfield(RDP,{'u','o'});
    for i = 1:N
        RDP(i).a  = mda.a;
        RDP(i).a0 = mda.a0;
        RDP(i).s  = RDP(i).s(:,1);
    end
    rng(m)
    RDP      = spm_MDP_VB_X(RDP,OPT);
    [F,Fu]   = spm_MDP_F(RDP);
    Rm(:,m)  = F(:);
    Rum(:,m) = Fu(:);

    % look for instances of BMR
    %----------------------------------------------------------------------
    vA    = spm_vec(A{1});
    for i = 1:N
       c(i,1) = corr(vA,(spm_vec(RDP(i).a{1}) > 0));
    end
    bmr(:,m) = diff(c);

    % preferred locations
    %----------------------------------------------------------------------
    for i = 1:N
        if hit(RDP(i)); r(i,m) = 1; else, r(i,m) = 0;  end
    end

    % find run of correct responses
    %----------------------------------------------------------------------
    for i = 1:N
        if hit(RDP(i)); r(i,m) = 1; else, r(i,m) = 0;  end
    end

    % find non-learners
    %----------------------------------------------------------------------
    R = spm_conv(r,2,0);
    H = spm_conv(h,2,0);

    % show results
    %----------------------------------------------------------------------
    spm_figure('GetWin','Figure 11');clf
    subplot(4,1,1)
    b = bar(mean(R(:,:),2)); set(b,'EdgeColor','w','FaceColor',[1 1 1]*.0),hold on
    b = bar(mean(H(:,:),2)); set(b,'EdgeColor','w','FaceColor',[1 1 1]*.8),hold off
    xlabel('trial'), ylabel('probability of correct'), axis([1/2 (N + 1/2) 1/3 1]);
    title('Average performance','Fontsize',16)

    subplot(4,1,2)
    spm_plot_ci(mean(Rm'),var(Rm')), hold on
    plot(mean(Fm'),'r'),              hold off
    xlabel('trial'), ylabel('free energy'); set(gca,'XLim',[1 N])
    title('Average free energy','Fontsize',16)

    subplot(4,1,3)
    spm_plot_ci(mean(Rum'),var(Rum')), hold on
    plot(mean(Fum'),'r'),              hold off
    xlabel('trial'), ylabel('confidence');  set(gca,'XLim',[1 N])
    title('Average confidence','Fontsize',16)

    % show individual performance
    %----------------------------------------------------------------------
    subplot(4,1,4), image(32*(r'))
    xlabel('trial'), ylabel('subject')
    title('Aha moments','Fontsize',16)

    % plot model updates
    %----------------------------------------------------------------------
    hold on
    for i = 1:m
        j = find(bmr(:,i) > 0) + 1;
        try, plot(j(1),i,'.m','MarkerSize',32), end
        try, plot(j(2),i,'.r','MarkerSize',32), end
        j = find(bmr(:,i) < 0);
        plot(j, (j - j + i),'.b','MarkerSize',32)
    end
    hold off, drawnow

    save paper

end

return

% confidence - negatively over policies
%--------------------------------------------------------------------------
for i = 1:numel(MDP)
    p     = MDP(i).R;
    Fu(i) = sum(sum(p.*log(p)));
end

return

function spm_MDP_plot_moves(MDP)
for i = 1:numel(MDP)
    m(i) = sum(~~diff(MDP(i).u(3,:)));
end
h = bar(m); set(h,'EdgeColor','w','FaceColor',[1 1 1]*.9), hold on

% and correct responses
%--------------------------------------------------------------------------
hit   = @(MDP) any(MDP.o(3,:) == 2) & ~any(MDP.o(3,:) == 3);
for i = 1:numel(MDP)
    if hit(MDP(i))
        h = plot(i,0,'hr','MarkerSize',8);
        set(h,'MarkerFaceColor',[1 1/2 0])
    end
end
return

function spm_MDP_A_plot(A)
% assemble key parts of the likelihood array
%==========================================================================
for i = 1:3
    for j = 1:3;
        a{i,j} = squeeze(A{1}(:,i,:,j,4));
        a{i,j} = a{i,j}*diag(1./sum(a{i,j}));
    end
end
a = spm_cat(a);
imagesc(a);
title( 'Sample: left - center - right', 'FontSize',16)
ylabel('Rule: left - center - right','FontSize',14)
xlabel('Correct color', 'FontSize',14)
set(gca,'XTick',1:9)
set(gca,'YTick',1:12)
set(gca,'XTicklabel',repmat(['r','g','b'],[1 3])')
set(gca,'YTicklabel',repmat(['r','g','b',' '],[1 3])')
axis image

return

function spm_MDP_rule_plot(MDP)
% illustrates visual search graphically
%==========================================================================

% locations
%--------------------------------------------------------------------------
x{1} = [-1  0; 0  1; 1  0; 0  0];
x{2} = [-1 -1; 0 -1; 1 -1; 0 -2]/2;
col  = {'r','g','b','c'};

% plot cues
%--------------------------------------------------------------------------
if strcmp('replace',get(gca,'Nextplot'))

    % plot cues
    %----------------------------------------------------------------------
    s     = MDP.s;hold off
    for i = 1:length(MDP.D{3})
        a = MDP.A{1}(:,s(1),s(2),i,1);
        j = find(rand < cumsum(a),1);
        plot(x{1}(i,1),x{1}(i,2),['.',col{j}],'MarkerSize',32), hold on
    end

    % plot choices
    %----------------------------------------------------------------------
    for i = 1:length(MDP.D{4})
        a = find(MDP.A{3}(:,s(1),s(2),4,i));
        if a == 2
            plot(x{2}(i,1),x{2}(i,2),['.m'],'MarkerSize',32), hold on
        end
        plot(x{2}(i,1),x{2}(i,2),['.',col{i}],'MarkerSize',16), hold on

    end
    axis([-2 2 -2 2]);

end

% Extract and plot eye movements and choice
%--------------------------------------------------------------------------
for i = 1:numel(MDP.o(2,:))
    X(i,:) = x{1}(MDP.o(2,i),:);
end
plot(X(:,1),X(:,2),'k')
for i = 1:numel(MDP.s(4,:))
    X(i,:) = x{2}(MDP.s(4,i),:);
end
plot(X(:,1),X(:,2),'k')
axis([-2 2 -2 2]);

return

author: spm


text:
classdef FCP3 < PROBLEM
% <multi> <real> <constrained>
% Benchmark constrained MOP proposed by Jiawei Yuan

%------------------------------- Reference --------------------------------
% J. Yuan, H. Liu, Y. Ong, and Z. He, Indicator-based evolutionary
% algorithm for solving constrained multi-objective optimization problems,
% IEEE Transactions on Evolutionary Computation, 2022, 26(2): 379-391.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Jiawei Yuan

    methods
        %% Initialization
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D)
                obj.D = 30;
            end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            g = 1 + 9*mean(PopDec(:,2:end),2);
            t = mod(floor(100*g),2);
            g = g + t.*(g-9).^2;
            PopObj(:,1) = cos(0.5*pi*PopDec(:,1)).*g;
            PopObj(:,2) = sin(0.5*pi*PopDec(:,1)).*g;
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            g = 1 + 9*mean(PopDec(:,2:end),2);
            t = mod(floor(100*g),2);
            g = g + t.*(g-9).^2;
            Dis    = abs(9-g);
            %%%%% Type-II constraints
            y1     = Dis.^2-0.25;
            y2     = 1./(Dis+1e-6).*(1.2+sin(Dis*pi));
            PopCon = min([y1,y2],[],2);
        end
        %% Sample reference points on Pareto front
        function P = GetOptimum(obj,N)
            t = 0.5*pi*(0:1/N:1)';
            P=8.5*[cos(t),sin(t)];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [coef]=comp_dwiltiv(coef2,a)
%COMP_DWILTIV  Compute Discrete Wilson transform type IV.
%

M=size(coef2,1)/2;
N=size(coef2,2);
W=size(coef2,3);
L=N*a;

coef=zeros(M,N,W,assert_classname(coef2));

% --- m is even ---------
coef(1:2:M,1:2:N,:)= 1/sqrt(2)*(exp(-i*pi/4)*coef2(1:2:M,1:2:N,:)+exp(-i*pi*3/4)*coef2(2*M:-2:M+1,1:2:N,:));
coef(1:2:M,2:2:N,:)= 1/sqrt(2)*(exp(i*pi/4)*coef2(1:2:M,2:2:N,:)+exp(i*pi*3/4)*coef2(2*M:-2:M+1,2:2:N,:));

% --- m is odd ----------
coef(2:2:M,1:2:N,:)= 1/sqrt(2)*(exp(i*pi/4)*coef2(2:2:M,1:2:N,:)+exp(i*pi*3/4)*coef2(2*M-1:-2:M+1,1:2:N,:));
coef(2:2:M,2:2:N,:)= 1/sqrt(2)*(exp(-i*pi/4)*coef2(2:2:M,2:2:N,:)+exp(-i*pi*3/4)*coef2(2*M-1:-2:M+1,2:2:N,:));

coef=reshape(coef,M*N,W);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% ***************************************************************************
% Linear Interpolation
% ***************************************************************************
% Author: Chaobin
% Email:  chaubyZou@163.com
% Date: October 2020
% ***************************************************************************
% Language: Matlab
% Also available in: Python
% Required library: None
% ***************************************************************************

classdef LinearInterpolation < handle
    properties
        name = 'linear interpolation';
        q_via = [];
        t_via = [];
    end

    methods
        % crate the objective
        % name: string
        % q_via: N x 3 array
        % t_via: N x 1 array
        function obj = LinearInterpolation(name, q_via, t_via)
            obj.name = name;
            obj.q_via = q_via;
            obj.t_via = t_via;
            if size(q_via(:,1)) ~= length(t_via)
                error('The q_via and t_via must have a same length');
            end
        end

        % linar interpolation with two data points
        % q0: the first data point
        % q1: the second data point
        % t0: the time of the first data point
        % t1: the time of the second data point
        % a0, a1: parameters
        function [a0, a1] = linear(obj, q0, q1, t0, t1)
            if abs(t0 - t1) < 1e-6
                error('t0 and t1 must be different');
            end
            a0 = q0;
            a1 = (q1 - q0)/(t1 - t0);
        end

        % linar interpolation for all data points
        % t: specified time
        % q: 1 x 3 array, output of the interpolation at time t
        function q = getPosition(obj, t)
            if (t < obj.t_via(1)) || (t > obj.t_via(end))
                error('The specific time error, time ranges error');
            end

            j = find(obj.t_via >= t, 1, 'first'); % find the index of t1
            if j == 1
                i = 1;
                j = 2;
            else
                i = j-1;
            end

            % position
            q0 = obj.q_via(i,1);
            t0 = obj.t_via(i);
            q1 = obj.q_via(j,1);
            t1 = obj.t_via(j);
            [a0, a1] = obj.linear(q0, q1, t0, t1);
            q(1, 1) = a0 + a1*(t - t0);

            % velocity
            q(1, 2) = a1;

            % acceleration
            q(1, 3) = 0; % for linear model, the acceleration is infinite, here we set to zero
        end
    end
end

author: chauby

repo: PolynomialInterpolation


text:
clear;

% Sample Code to test the OptSpace program :

fprintf(1,'Generating Matrix...\n');
n = 1000;
m = 1000;
r = 3;
tol = 1e-8 ;

eps = 10*r*log10(n);
% Generate the matrix with the given parameters
U = randn(n,r);
V = randn(m,r);
Sig = eye(r) ;
M0 = U*Sig*V';

% Select the entries independently with prob. eps/sqrt(mn)
E = 1 - ceil( rand(n,m) - eps/sqrt(m*n)  ) ;
M_E = sparse(M0.*E) ;

fprintf(1,'Testing without noise...\n');
% Call the OptSpace function
[X S Y dist] = OptSpace(sparse(M_E),[],[],tol);

% Compute the Frobenius norm
err_nonoise = norm(X*S*Y' - M0,'fro')/sqrt(m*n);

% Add noise
sig = .01 ;
M = M0 + sig*randn(n,m);
M_E = sparse(M.*E) ;

fprintf(1,'Testing with noise...\n');
% Call the OptSpace function
[X S Y dist] = OptSpace(sparse(M_E),[],20,tol);

% Compute the Frobenius norm
err_noise = norm(X*S*Y' - M0,'fro')/sqrt(m*n);

fprintf(1,'RMSE (without noise)           : %e\n',err_nonoise);
fprintf(1,'RMSE (with noise var %f) : %e\n',sig,err_noise);

%%
numr = size(M,1);
numc = size(M,2);
I = randi([0 1],numr,numc);
MI = M.*I;
show_2dvideo(MI,m,n);
tol = 1e-8;
[X,S,Y,dist] = OptSpace(sparse(MI),[],20,tol);
L = X*S*Y';
show_2dvideo(L,m,n);

author: andrewssobral

repo: lrslibrary


text:
function varargout = lap(varargin)
%LAP   Laplacian of a CHEBFUN2.
%   L = LAP(F) returns a CHEBFUN2 representing the Laplacian of F.
%
%   This is shorthand for LAPLACIAN(F).
%
% See also CHEBFUN2/LAPLACIAN.

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

% Call SEPARABLEAPPROX/LAP:
[varargout{1:nargout}] = lap@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [windows,window_radii,windowi2radiusi,nradii] = SetWindowParameters(...
  windows,window_offsets,...
  window_radii,min_window_radius,...
  max_window_radius,nwindow_radii)

if isempty(windows),

  % use radii and offsets if windows not input
  if isempty(window_radii),
    % use min, max, n if radii not input
    window_radii = unique(round(logspace(log10(min_window_radius+1),log10(max_window_radius+1),nwindow_radii)))-1;
  end
  if isempty(window_radii),
    error('window_radii is empty.');
  end
  if isempty(window_offsets),
    error('window_offsets is empty.');
  end

  % take all combinations of radii and offsets
  [all_radii,all_offsets] = meshgrid(window_radii,window_offsets);
  all_radii = all_radii(:);
  all_offsets = all_offsets(:);

  % offsets are fractions of radii
  all_offsets = round(all_offsets.*(max(all_radii,1)));

  % window is the pair of radius, offset
  windows = [all_radii,all_offsets];

  % make sure these are unique -- rounding to nearest frame might make them
  % not unique
  windows = unique(windows,'rows');

  % window i for frame t is
  % [t-windows(i,1)+windows(i,2), t+windows(i,1)+windows(i,2)]

end

if isempty(windows),
  error('windows is empty.');
end
if size(windows,2) ~= 2,
  error('windows must be nwindows x 2.');
end

[window_radii,~,windowi2radiusi] = unique(windows(:,1));
nradii = numel(window_radii);


text:
function [seg, SEG_l] = idpLinSegment(l,t)

% IDPLINSEGMENT  IDP line endpoints.
%   IDPLINSEGMENT(L,T) returns the 3D segment corresponding to the IDP line
%   L at abscissas T = [t1;t2].
%
%   [s,S_l] = IDPLINSEGMENT(...) returns the Jacobian wrt L.
%
%   See also IDPLINENDPOINTS, IDPLIN2SEG, IDPLIN2IDPPNTS.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% abscissas
t1 = t(1);
t2 = t(2);

% support 3d segment
[s, S_l] = idpLin2seg(l);

% support points
p1 = s(1:3);
p2 = s(4:6);

P1_l = S_l(1:3,:);
P2_l = S_l(4:6,:);

% 3d endpoints
e1 = (1-t1)*p1 + t1*p2;
e2 = (1-t2)*p1 + t2*p2;

E1_p1 = 1-t1;
E1_p2 = t1;
E2_p1 = 1-t2;
E2_p2 = t2;

E1_l = E1_p1*P1_l + E1_p2*P2_l;
E2_l = E2_p1*P1_l + E2_p2*P2_l;

% 3d segment
seg   = [e1;e2];
SEG_l = [E1_l;E2_l];

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
classdef testPrintingComparator < handle

    properties (Access = protected, Abstract)
        fileOutputName
        iter
    end

    methods (Access = protected)

        function hasPassed = hasPassed(obj)
            hasMshChanged = obj.compareFile('.msh');
            hasResChanged = obj.compareFile('.res');
            filesHaveChanged = hasMshChanged || hasResChanged;
            hasPassed = ~filesHaveChanged;
        end

        function selectComputedVar(obj)
        end

    end

    methods (Access = private)

       function sF = obtainSavedPrintedFile(obj,ext)
           out   = obj.fileOutputName;
           fullOutSavedName = [out,'.flavia',ext];
           sF = fullfile('tests','PrintingTests',out,fullOutSavedName);
        end

        function oF = obtainOutPutFile(obj,ext)
           out   = obj.fileOutputName;
           fullOutName = [out,num2str(obj.iter),'.flavia',ext];
           oF = fullfile('Output',out,fullOutName);
        end

        function hasChanged = compareFile(obj,ext)
            sF = obj.obtainSavedPrintedFile(ext);
            oF = obj.obtainOutPutFile(ext);
            hasChanged = FileComparator().areFilesDifferent(sF,oF);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function value = mono_upto_enum ( m, n )

%*****************************************************************************80
%
%% MONO_UPTO_ENUM enumerates monomials in M dimensions of degree up to N.
%
%  Discussion:
%
%    For M = 2, we have the following values:
%
%    N  VALUE
%
%    0    1
%    1    3
%    2    6
%    3   10
%    4   15
%    5   21
%
%    In particular, VALUE(2,3) = 10 because we have the 10 monomials:
%
%      1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the maximum degree.
%
%    Output, integer VALUE, the number of monomials in D variables,
%    of total degree N or less.
%
  value = nchoosek ( n + m, n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_laplacian( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e6*pref.techPrefs.chebfuneps;

% Test with different parity of m,p
% Example 1:
f = ballfun(@(x,y,z)x.^2+y.^2+z.^2);
g = laplacian(f);
exact = ballfun(@(x,y,z)6);
pass(1) = norm( g - exact ) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end


text:
function [ y, ierror ] = y_check_republican ( y )

%*****************************************************************************80
%
%% Y_CHECK_REPUBLICAN checks a Republican year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, the year, which must not be 0.
%
%    Output, integer IERROR, is 0 if Y is legal, and 1 otherwise.
%
  if ( 0 < y )
    ierror = 0;
  else
    ierror = 1;
  end

  return
end