text:
%This Matlab script can be used to reproduce Figure 7.18 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%% Define parameters
M_max = 15;          % number of antenna columns/rows
d_H = 1/2;          % horizontal antenna spacing
h_BS = 25;          % height of the antenna array
h_UE = 1.5;         % height of the UE
CellRadius = 500;   % cell radius
K = 500;            % number of UEs
fc = 2.6e9;         % center frequency
bw = 20e6;          % bandwidth
N = 1024;           % number of subcarriers

% Compute derived parameters
h = h_BS - h_UE;    % effective height
d_V = 1/2;          % vertical antenna spacing
d = d_V;

%Simulate new channel reponses or load existing results (the latter is the
%preset)
simulate = false;

%% Simulate channels
if (simulate)
    % Draw random positions
    A = rand(1,K)*CellRadius^2;
    B = rand(1,K)*pi*2/3 - pi/3;
    POS = repmat(sqrt(A)',[1,2]).* [cos(B); sin(B)]'; %UE positions uniformly distributed on a disc of radius R

    % Create a new layout
    s = simulation_parameters; % Create object with general parameters
    s.center_frequency = fc;
    s.sample_density = 1;
    s.use_absolute_delays = 1;

    l = layout(s);  % Create new layout from general parameters
    l.no_rx = K;  % Generate KUEs
    l.rx_array.generate('omni'); % use omnidirectional antennas at the UEs

    l.randomize_rx_positions(CellRadius, h_UE, h_UE, 1); % randomly distribute UEs within a disc of radius R and height 1.5m
    for k=1:K
        l.rx_position(:,k) = [POS(k,1); POS(k,2);  h_UE]; % use the same positions as for the other simulations
    end

    for i=1:l.no_rx % for each receiver
        l.track(i).generate('linear',0,0) % define a linear track consiting of only one position
        l.track(i).scenario = '3GPP_3D_UMa_LOS'; % select the Urban Macrocell NLOS scenario for the link
    end

    % Simulate channel for planar array
    l.tx_array.generate('3gpp-3d', 1, M_max, M_max, fc, 1, 0, d);
    [h_channel, h_parset, h_cb] = l.get_channels();
    h_freq = h_channel.fr(bw,N);
    H_plan = zeros(K,M_max^2, N);
    for k=1:K
        H_plan(k,:,:) = squeeze(h_freq{k});
    end

    % Simulate channel for horizontal array
    l.tx_array.generate('3gpp-3d', 1, 1, M_max^2, fc, 1, 0, d);
    [h_channel, h_parset, h_cb] = l.get_channels();
    h_freq = h_channel.fr(bw,N);
    H_hori = zeros(K,M_max^2, N);
    for k=1:K
        H_hori(k,:,:) = squeeze(h_freq{k});
    end

    % Simulate channel for vertical array
    l.tx_array.generate('3gpp-3d', 1, M_max^2, 1, fc, 1, 0, d);
    [h_channel, h_parset, h_cb] = l.get_channels();
    h_freq = h_channel.fr(bw,N);
    H_vert = zeros(K,M_max^2, N);
    for k=1:K
        H_vert(k,:,:) = squeeze(h_freq{k});
    end
end

%% Compute channel orthogonality
if (~simulate)
    load('section7_orthogonality_3GPP.mat')
else

    IT = 100000;
    COR_3GPP_PLAN = zeros(1,M_max);
    COR_3GPP_HORI = zeros(1,M_max);
    COR_3GPP_VERT = zeros(1,M_max);

    for it=1:IT

        %Output simulation progress
        disp([num2str(it) ' iterations out of ' num2str(IT)]);

        % Draw random position indices
        posInd = randperm(K,2);

        % For all antenna size
        for m = 1:M_max
            M_H = m;
            M_V = m;
            M = M_H*M_V;
            indices = repmat((0:M_H-1)*M_max + 1, [M_V,1]) + repmat((0:M_V-1)', [1,M_H]);
            indices = indices(:);
            for n=1:N
                % Planar array
                h1 = H_plan(posInd(1),indices,n);
                h2 = H_plan(posInd(2),indices,n);
                COR_3GPP_PLAN(m) = abs(h1*h2')^2 / (norm(h1)^2*norm(h2)^2)/IT/N + COR_3GPP_PLAN(m);

                % Horizontal array
                h1 = H_hori(posInd(1),1:m^2,n);
                h2 = H_hori(posInd(2),1:m^2,n);
                COR_3GPP_HORI(m) = abs(h1*h2')^2 / (norm(h1)^2*norm(h2)^2)/IT/N + COR_3GPP_HORI(m);

                % Vertical array
                h1 = H_vert(posInd(1),1:m^2,n);
                h2 = H_vert(posInd(2),1:m^2,n);
                COR_3GPP_VERT(m) = abs(h1*h2')^2 / (norm(h1)^2*norm(h2)^2)/IT/N + COR_3GPP_VERT(m);
            end
        end
    end
end

%% Load channel measurements
hor_out=load('section7_measurement_horizontal.mat');
vert_out=load('section7_measurement_vertical.mat');
planar_square = load('section7_measurement_planar.mat');

x_measured=(1:8).^2;

%Compute average UE correlation
COR_MEAS_HORI = mean(hor_out.crossco).^2;
COR_MEAS_PLAN = mean(planar_square.crossco).^2;
COR_MEAS_VERT = mean(vert_out.crossco).^2;

%% Plot the simulation results
figure;
hold on; box on;

xvals = (1:M_max).^2;

plot(xvals, COR_3GPP_VERT, 'LineWidth', 1, 'Color', 'k', 'LineStyle', '-', 'Marker', 'x');
plot(x_measured, COR_MEAS_VERT(x_measured), 'LineWidth', 1, 'Color', 'k', 'LineStyle', '--', 'Marker', 'x');

plot(xvals, COR_3GPP_PLAN, 'LineWidth', 1, 'Color', 'b', 'LineStyle', '-', 'Marker', 'o');
plot(x_measured, COR_MEAS_PLAN, 'LineWidth', 1, 'Color', 'b', 'LineStyle', '--', 'Marker', 'o');

plot(x_measured, COR_MEAS_HORI(x_measured), 'LineWidth', 1, 'Color', 'r', 'LineStyle', '--','Marker', 'square');
plot(xvals, COR_3GPP_HORI, 'LineWidth', 1, 'Color', 'r', 'LineStyle', '-','Marker', 'square');

COR_IID = 1./xvals;
plot(xvals, COR_IID, 'LineWidth', 1, 'Color', 'k', 'LineStyle', '-.');

legend({'3GPP Vertical', 'Meas. Vertical', '3GPP Planar', 'Meas. Planar', 'Meas. Horizontal', '3GPP Horizontal', 'i.i.d.'}, 'location', 'southwest');
set(gca, 'XScale', 'log', 'YScale', 'log');
xlabel('Number of antennas (M)');
ylabel('Average UE correlation');
set(gca,'xtick',[1,4,9,16,25,36,49,64,121,225]);
xlim([1,M_max^2]);
ylim([COR_IID(end), 1]);

author: emilbjornson

repo: massivemimobook


text:
function G=GradSymbolic(p)
%create gradient profile

global VMgd

% Initialize parameters
GradLine=p.GradLine;
GradXEqu=p.GradXEqu;
GradYEqu=p.GradYEqu;
GradZEqu=p.GradZEqu;

% Initialize display grid
X=VMgd.xgrid;
Y=VMgd.ygrid;
Z=VMgd.zgrid;

try
    eval(['G(:,:,:,1)=' GradXEqu ';']);
    eval(['G(:,:,:,2)=' GradYEqu ';']);
    eval(['G(:,:,:,3)=' GradZEqu ';']);
catch me
    G = zeros(size(VMgd.xgrid));
    error_msg{1,1}='ERROR!!! Creating gradient profile using symbolic equation fails!';
    error_msg{2,1}=me.message;
    errordlg(error_msg);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function cMap=gjet(varargin)

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

cMap=[51  105 232;... %Blue
      0   153 37;... %Green
      238 178 17;... %Yelow
      213 15  37;... %red
    ];
cMap=cMap./255;

[cMap]=resampleColormap(cMap,n);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% script for computing window features

% set up path
if ispc,
  JCtrax_path = 'E:\Code\JCtrax';
  FlyBowlAnalysis_path = 'E:\Code\FlyBowlAnalysis';
  rootdatadir = 'E:\Code\Jdetect\larva\fly_data\TrainingData';
  settingsdir = 'E:\Code\FlyBowlAnalysis\settings';
else
  JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
  FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
  settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
end
addpath(fullfile(JCtrax_path,'misc'));
addpath(fullfile(JCtrax_path,'filehandling'));
addpath(FlyBowlAnalysis_path);

%% load test data

labelmatnames = {
  'perframe_labeledsharpturns_movie01_pBDPGAL4U_TrpA_Rig1Plate10BowlA_20110202T105734_fly1_11.mat'
  'perframe_labeledsharpturns_movie01_pBDPGAL4U_TrpA_Rig1Plate10BowlA_20110202T105734_fly2_08.mat'
  'perframe_labeledsharpturns_movie02_pBDPGAL4U_TrpA_Rig2Plate14BowlA_20110202T110111_fly1_01.mat'
  'perframe_labeledsharpturns_movie02_pBDPGAL4U_TrpA_Rig2Plate14BowlA_20110202T110111_fly2_15.mat'
  'perframe_labeledsharpturns_movie03_pBDPGAL4U_TrpA_Rig2Plate14BowlB_20110202T110116_fly1_12.mat'
  'perframe_labeledsharpturns_movie03_pBDPGAL4U_TrpA_Rig2Plate14BowlB_20110202T110116_fly2_03.mat'
  'perframe_labeledsharpturns_movie07_GMR_14G03_AE_01_TrpA_Rig2Plate14BowlA_20110202T102516_fly2_17.mat'
  'perframe_labeledsharpturns_movie08_GMR_14H07_AE_01_TrpA_Rig1Plate10BowlB_20110202T084231_fly1_09.mat'
  'perframe_labeledsharpturns_movie09_GMR_16E02_AE_01_TrpA_Rig1Plate10BowlC_20110202T140143_fly1_09.mat'
  };
nexps = numel(labelmatnames);

paramsfilename = 'SharpturnPerFrameParams_AllParams.xml';

%% set parameters
[params,cellparams] = ReadPerFrameParams(paramsfilename);
fns = fieldnames(cellparams);

%% do it

X = cell(1,nexps);
Y = cell(1,nexps);
feature_names = {};

for expi = 1:nexps,

  labelfilename = fullfile(rootdatadir,labelmatnames{expi});
  fprintf('Computing window features for file %s...\n',labelmatnames{expi});
  if ~exist(labelfilename,'file'),
    warning('File %s does not exist, skipping',labelfilename);
    continue;
  end
  labeldata = load(labelfilename);
  for i = 1:numel(fns),
    fn = fns{i};
    fprintf('Per-frame feature %s...\n',fn);
    j = find(strcmp(labeldata.perframefns,fn));
    if isempty(j),
      error('Unknown perframe field %s',fn);
    end
    [x_curr,feature_names_curr] = ...
      ComputeWindowFeatures(labeldata.perframedata{j},cellparams.(fn){:});
    if expi == 1,
      feature_names = [feature_names,cellfun(@(s) [{fn},s],feature_names_curr,'UniformOutput',false)];
    end

    % crop out labeled part of video
    x_curr = x_curr(:,labeldata.trk_labelstart:labeldata.trk_labelend-1);

    % store
    X{expi} = [X{expi};x_curr];
  end

  % store labels
  Y{expi} = labeldata.label;

end
nfeatures = numel(feature_names);

%% display names for stats

statnames = cell(1,nfeatures);
for i = 1:nfeatures,
  fn = '';
  for j = 1:numel(feature_names{i}),
    if j == 1,
    elseif mod(j,2) == 0,
      fn = [fn,', '];
    else
      fn = [fn,'='];
    end
    if ischar(feature_names{i}{j}),
      fn = [fn,feature_names{i}{j}];
    else
      fn = [fn,num2str(feature_names{i}{j})];
    end
  end
  statnames{i} = fn;
end

%% combine stats to plot everything

X_combine = X{1};
Y_combine = Y{1};
expi_combine = ones(1,size(X{1},2));
for expi = 2:nexps,
  if isempty(X{expi}),
    continue;
  end
  X_combine = [X_combine,nan([nfeatures,1]),X{expi}];
  Y_combine = [Y_combine,nan,Y{expi}];
  expi_combine = [expi_combine,nan,expi+zeros(1,size(X{expi},2))];
end
N = numel(Y_combine);

%% interactive figure in which we can plot some of the window features

handles = struct;
handles.fig = 1;
axespos = [.05,.1,.7,.85];
border = .01;
listpos = [axespos(1)+axespos(3)+border,...
  axespos(2),...
  1-(axespos(1)+axespos(3)+2*border),...
  axespos(4)];
figure(handles.fig);
clf;
handles.ax = axes('Parent',handles.fig,'Units','Normalized','Position',axespos);
handles.listbox = uicontrol(handles.fig,'Style','listbox',...
  'Units','normalized','Position',listpos,...
  'String',statnames,...
  'Max',nfeatures,'Min',0,'Value',[]);
xlim = [0,N+1];
ylim = [min(X_combine(:)),max(X_combine(:))];
set(handles.fig,'ToolBar','figure');
set(handles.ax,'XLim',xlim,'YLim',ylim);
hold(handles.ax,'on');
[i0,i1] = get_interval_ends(Y_combine==1);
i1 = i1 - 1;
for i = 1:numel(i0),
  patch([i0(i)-.5,i1(i)+.5,i1(i)+.5,i0(i)-.5,i0(i)-.5],ylim([1,1,2,2,1]),...
    [1,.7,.7],'EdgeColor','none');
end
[i0,i1] = get_interval_ends(isnan(Y_combine));
i1 = i1 - 1;
for i = 1:numel(i0),
  patch([i0(i)-.5,i1(i)+.5,i1(i)+.5,i0(i)-.5,i0(i)-.5],ylim([1,1,2,2,1]),...
    [.5,.5,.5],'EdgeColor','none');
end
handles.title = title('');
handles.data = [];
setappdata(handles.listbox,'handles',handles);
setappdata(handles.listbox,'x',X_combine);
set(handles.listbox,'Callback',@DebugComputeWindowFeatures_ListBoxCallback);

%% interactive figure in which we can histogram some of the window features

nbins = 100;
edges = nan(nfeatures,nbins+1);
centers = nan(nfeatures,nbins);
frac = nan(nfeatures,nbins,2);
for i = 1:nfeatures,
  edges(i,:) = linspace(min(X_combine(i,:)),max(X_combine(i,:)),nbins+1);
  centers(i,:) = (edges(i,1:end-1)+edges(i,2:end))/2;
  edges(i,end) = inf;
  for j = 1:2,
    counts = histc(X_combine(i,Y_combine==j-1),edges(i,:));
    frac(i,:,j) = counts(1:end-1) / sum(counts(1:end-1));
  end
end

handles = struct;
handles.fig = 2;
axespos = [.05,.1,.7,.85];
border = .01;
listpos = [axespos(1)+axespos(3)+border,...
  axespos(2),...
  1-(axespos(1)+axespos(3)+2*border),...
  axespos(4)];
figure(handles.fig);
clf;
handles.ax = axes('Parent',handles.fig,'Units','Normalized','Position',axespos);
handles.listbox = uicontrol(handles.fig,'Style','listbox',...
  'Units','normalized','Position',listpos,...
  'String',statnames,...
  'Max',1,'Min',0,'Value',1);
set(handles.fig,'ToolBar','figure');
hold(handles.ax,'on');
handles.data = [];
setappdata(handles.listbox,'handles',handles);
setappdata(handles.listbox,'frac',frac);
setappdata(handles.listbox,'centers',centers);
setappdata(handles.listbox,'legends',{'none','sharpturn'});
DebugComputeWindowFeatures_ListBoxCallback_Hist(handles.listbox);
set(handles.listbox,'Callback',@DebugComputeWindowFeatures_ListBoxCallback_Hist);

author: kristinbranson

repo: JAABA


text:
%% test modulate for all oasis functions.
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]}; % colors
plot_cvx = false;

%% example 7:  constrained-foopsi, AR1
g = 0.95;         % AR coefficient
noise = .3;
T = 3000;
framerate = 30;
firerate = 0.5;
b = 0;              % baseline
N = 1;              % number of trials
seed = 13;          % seed for genrating random variables
[y, true_c, true_s] = gen_data(g, noise, T, framerate, firerate, b, N, seed);

% cvx solution
[c_cvx, s_cvx] = constrained_foopsi(y, g, noise);
% case 1: all parameters are known
[c_oasis, s_oasis] = deconvolveCa(y, 'ar1', g, 'constrained', 'sn', noise);  %#ok<*ASGLU>

figure('name', 'constrained-FOOPSI, AR1, known: g, sn', 'papersize', [15, 4]);
plot_cvx = true;
show_results;
plot_cvx = false;

% case 2: nothing is known, estimate g with auto-correlation method
[c_oasis, s_oasis,options] = deconvolveCa(y, 'ar1', 'constrained');

fprintf('true gamma:        %.3f\n', g);
fprintf('estimated gamma:   %.3f\n', options.pars);

figure('name', 'FOOPSI, AR1, estimated: g, sn', 'papersize', [15, 4]);
show_results;

% case 3: nothing is know, estimate g with auto-correlation method first
% and then update it to minimize the RSS
[c_oasis, s_oasis, options] = deconvolveCa(y, 'ar1', 'constrained', ...
    'optimize_pars');

fprintf('true gamma:        %.3f\n', g);
fprintf('estimated gamma:   %.3f\n', options.pars);

figure('name', 'FOOPSI, AR1, estimated: g, sn, update:g', 'papersize', [15, 4]);
show_results;

% case 4: nothing is know, estimate g with auto-correlation method first
% and then update it to minimize the RSS, the baseline is also unknown
true_b = 0.5;
[c_oasis, s_oasis, options] = deconvolveCa(y+true_b, 'ar1', g,...
    'constrained','optimize_b', 'sn', noise);
fprintf('true gamma:        %.3f\n', g);
fprintf('estimated gamma:   %.3f\n', options.pars);
fprintf('true b:       %.3f\n', true_b);
fprintf('estimated b:       %.3f\n', options.b);
fprintf('tuning parameter:  %.3f\n', options.lambda);

figure('name', 'FOOPSI, AR1, estimated: g, sn, lambda', 'papersize', [15, 4]);
show_results;

% case 5: nothing is know, estimate g with auto-correlation method first
% and then update it to minimize the RSS, the baseline is also unknown
true_b = 0.5;
[c_oasis, s_oasis, options] = deconvolveCa(y+true_b, 'ar1',...
    'constrained','optimize_b', 'optimize_pars');
fprintf('true gamma:        %.3f\n', g);
fprintf('estimated gamma:   %.3f\n', options.pars);
fprintf('estimated b:       %.3f\n', options.b);
fprintf('tuning parameter:  %.3f\n', options.lambda);

figure('name', 'FOOPSI, AR1, estimated: g, sn, lambda, update:g', 'papersize', [15, 4]);
show_results;

%%

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [A,b] = pwamodel(F,x)

[dummy,aux,s,model] = export(F,[]);

A = -model.F_struc(:,2:end);
b = model.F_struc(:,1);

keep_these = find(ismember(aux.used_variables,getvariables(x)));

A = [A(:,keep_these) A(:,setdiff(1:size(A,2),keep_these))];
[A,b] = fourier_motzkin(A,b,length(keep_these));

function [A,b] = fourier_motzkin(A,b,m)

while size(A,2)>m
    [A,b] = fourier_motzkin_1(A,b,m);
    [aux,i] = unique([A b],'rows');
    A = A(i,:);
    b = b(i);
end

function [Aout,bout] = fourier_motzkin_1(A,b,m)

for i = m+1:size(A,2)
    less = find(A(:,i)>0);
    larger = find(A(:,i)<0);
    t(i-m) =length(less)*length(larger);
end

[minn,remove] = min(t);
remove = remove+m;
keep = setdiff(1:size(A,2),remove);

less = find(A(:,remove)>0);
larger = find(A(:,remove)<0);
notinvolved = find(A(:,remove)==0);

Aout = A(notinvolved,keep);
bout = b(notinvolved);

for i = 1:length(less)
    for j = 1:length(larger)
        Aout = [Aout;A(less(i),keep)/abs(A(less(i),remove)) + A(larger(j),keep)/abs(A(larger(j),remove))];
        bout = [bout;b(less(i))+b(larger(j))];
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function r8mat_solve2_test ( )

%*****************************************************************************80
%
%% R8MAT_SOLVE2_TEST tests R8MAT_SOLVE2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 4;

  a1 = [ ...
    1.0, 3.0; ...
    2.0, 4.0 ]';
  a2 = [ ...
    2.0, 1.0, 1.0; ...
    1.0, 1.0, 0.0; ...
    1.0, 0.0, 1.0 ]';
  a3 = [ ...
    1.0, 2.0, 1.0, 3.0; ...
    0.0, 1.0, 2.0, 1.0; ...
    0.0, 0.0, 3.0, 2.0; ...
    1.0, 3.0, 0.0, 1.0 ]';
  a4 = [ ...
    2.0, 1.0, 3.0; ...
    4.0, 2.0, 6.0; ...
    1.0, 4.0, 5.0 ]';
  b1 = [ 5.0, 11.0 ];
  b2 = [ 4.0, 2.0, 2.0 ];
  b3 = [ 5.0, 11.0, 16.0, 15.0 ];
  b4 = [ 13.0, 17.0, 20.0 ];
  n_test = [ 2, 3, 4, 3 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_SOLVE2_TEST\n' );
  fprintf ( 1, '  R8MAT_SOLVE2 is a linear solver.\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    n = n_test ( test );

    if ( test == 1 )
      a(1:n,1:n) = a1(1:n,1:n);
      b(1:n) = b1(1:n);
    elseif ( test == 2 )
      a(1:n,1:n) = a2(1:n,1:n);
      b(1:n) = b2(1:n);
    elseif ( test == 3 )
      a(1:n,1:n) = a3(1:n,1:n);
      b(1:n) = b3(1:n);
    elseif ( test == 4 )
      a(1:n,1:n) = a4(1:n,1:n);
      b(1:n) = b4(1:n);
    end

    r8vec_print ( n, b, '  Right hand side:' );

    [ x, ierror ] = r8mat_solve2 ( n, a, b );

    fprintf ( 1, '\n' );
    if ( ierror == 0 )
      fprintf ( 1, '  The system is nonsingular.\n' );
    elseif ( ierror == 1 )
      fprintf ( 1, '  The system is singular, but consistent.\n' );
    elseif ( ierror == 2 )
      fprintf ( 1, '  The system is singular and inconsistent.\n' );
    end

    r8vec_print ( n, x, '  Computed solution:' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function K = sdlfmvXsdrbfKernCompute(sdlfmvKern, sdrbfKern, t1, t2)

% SDLFMVXSDRBFKERNCOMPUTE Cross kernel between a SDLFMV and a SDRBF kernels.
% FORMAT
% DESC computes cross kernel terms between the velocity switching dynamical
% LFM and the switching dynamicl RBF.
% ARG sdlfmvKern : the kernel structure associated with the velocity of the
% first SDLFM kernel.
% ARG sdrbfKern : the kernel structure associated with the SDRBF kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between the velocity switching dynamical
% LFM and the switching dynamicl RBF.
% ARG sdlfmvKern : the kernel structure associated with the velocity of the
% first SDLFM kernel.
% ARG sdrbfKern : the kernel structure associated with the SDRBF kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : sdlfmvKernParamInit, sdlfmvKernCompute, sdlfmKernParamInit
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 4
    t2 = t1;
end

K = sdlfmXsdrbfKernCompute(sdlfmvKern, sdrbfKern, t1, t2, 'Vel');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f = beale ( x )

%*****************************************************************************80
%
%% BEALE computes the Beale function.
%
%  Discussion:
%
%    This function has a global minimizer:
%
%      X = ( 3.0, 0.5 )
%
%    for which
%
%      F(X) = 0.
%
%    For a relatively easy computation, start the search at
%
%      X = ( 1.0, 1.0 )
%
%    A harder computation starts at
%
%      X = ( 1.0, 4.0 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Evelyn Beale,
%    On an Iterative Method for Finding a Local Minimum of a Function
%    of More than One Variable,
%    Technical Report 25,
%    Statistical Techniques Research Group,
%    Princeton University, 1958.
%
%    Richard Brent,
%    Algorithms for Minimization with Derivatives,
%    Dover, 2002,
%    ISBN: 0-486-41998-3,
%    LC: QA402.5.B74.
%
%  Parameters:
%
%    Input, real X(N,2), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f1 = 1.5   - x(:,1) .* ( 1.0 - x(:,2)    );
  f2 = 2.25  - x(:,1) .* ( 1.0 - x(:,2).^2 );
  f3 = 2.625 - x(:,1) .* ( 1.0 - x(:,2).^3 );

  f = f1.^2 + f2.^2 + f3.^2;

  return
end


text:
%EXPORT_XLS Export grid and data in Excel format.
%
%   EXPORT_XLS( FILENAME, FEA, EXPR, SOLNUM, FID_LOG )
%
%   Export of grid points and corresponding solution/expression data
%   in Excel format. The output file will consist of a header line
%   with the space coordinate names, dependent variable names, and
%   given expressions (columns). This is followed by the actual data
%   rows where one line corresponds to the data for the specific grid
%   point.
%
%   FILENAME gives the output file name. FEA must be a valid FEA data
%   struct, EXPR is optionally a cell array with expressions to
%   process (the grid coordinates and dependent variables in the FEA
%   struct will be prepended to EXPR). SOLNUM is optionally the
%   solution number of which to evaluate the expressions. FID_LOG is
%   an optional log file handle for message output (negative for gui
%   output or empty for no output).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
load('menpo_mtcnn.mat')

% Find the ground truth bboxes
bboxes_gt = bboxes;

bboxes_det = dets;

non_detected = bboxes_det(:,3) == 0;

% Removing the outliers
widths_gt = bboxes_gt(:,3) - bboxes_gt(:,1);
widths_det = bboxes_det(:,3) - bboxes_det(:,1);

heights_gt = bboxes_gt(:,4) - bboxes_gt(:,2);
heights_det = bboxes_det(:,4) - bboxes_det(:,2);

tx_gt =  bboxes_gt(:,1);
ty_gt =  bboxes_gt(:,2);

tx_det = bboxes_det(:,1);
ty_det = bboxes_det(:,2);

bad_det_1 = abs(1 - widths_gt ./ widths_det) > 0.5;
bad_det_2 = abs(1 - heights_gt ./ heights_det) > 0.5;

bad_det_3 = abs((tx_gt - tx_det) ./ widths_det) > 0.4;
bad_det_4 = abs((ty_gt - ty_det) ./ heights_det) > 0.5;

non_detected = non_detected | bad_det_1 | bad_det_2 | bad_det_3 | bad_det_4;

% if the width is quite different from detection then it failed

bboxes_gt = bboxes_gt(~non_detected,:);
bboxes_det = bboxes_det(~non_detected,:);

%% some visualisations
% a = 1;
% plot(gt_labels(a,:,1), gt_labels(a,:,2), '.r');
% hold on;
% bbox = bboxes_detector(a,:);
% % bbox(2) = -bbox(2);
% rectangle('Position', bbox);
% hold off;
% axis equal;

% Want to find out what scaling and translation would lead to the smallest
% RMSE error between initialised landmarks and gt landmarks TODO

% Find the width and height mappings
widths_gt = bboxes_gt(:,3) - bboxes_gt(:,1);
widths_det = bboxes_det(:,3) - bboxes_det(:,1);

heights_gt = bboxes_gt(:,4) - bboxes_gt(:,2);
heights_det = bboxes_det(:,4) - bboxes_det(:,2);

s_width = widths_det \ widths_gt;
s_height = heights_det \ heights_gt;

tx_gt =  bboxes_gt(:,1);
ty_gt =  bboxes_gt(:,2);

tx_det = bboxes_det(:,1);
ty_det = bboxes_det(:,2);

s_tx = median((tx_gt - tx_det) ./ widths_det);
s_ty = median((ty_gt - ty_det) ./ heights_det);

%%
new_widths = widths_det * s_width;
new_heights = heights_det * s_height;
new_tx = widths_det * s_tx + tx_det;
new_ty = heights_det * s_ty + ty_det;

overlaps = zeros(numel(widths_det), 1);
new_overlaps = zeros(numel(widths_det), 1);

for i=1:numel(widths_det)
    bbox_gt = bboxes_gt(i,:);
    bbox_old = bboxes_det(i,:);
    overlaps(i) = overlap(bbox_gt, bbox_old);
    bbox_new = [new_tx(i), new_ty(i), new_tx(i) + new_widths(i), new_ty(i) + new_heights(i)];
    new_overlaps(i) = overlap(bbox_gt, bbox_new);
end

fprintf('Orig - %.3f, now - %.3f\n', mean(overlaps), mean(new_overlaps));

author: TadasBaltrusaitis

repo: OpenFace


text:
function data=syms2v(symbols,numberOfZeros)
%% syms2v : Symbols (nRows,N) to Vector of size (1,L)
% data=syms2v(symbols,numberOfZeros)
%% INPUTS: symbols, and numberOfZeros; as follows
% where L = number of the elements in the vector [data]
%% Outputs:
% data: is the original vector without the originally attached zeros
[r c]=size(symbols);
d=reshape(symbols,1,r*c);
data=d(1:end-numberOfZeros);

author: Sable

repo: mcbench-benchmarks


text:
%Metabolite
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 1/5/2011
classdef Metabolite < edu.stanford.covert.cell.sim.MoleculeCountState
    %Constants
    properties (Constant)
        optionNames             = {   %names of properties that are options
            'verbosity'
            'seed'
            };
        fixedConstantNames      = {   %names of process properties that are considered fixed constants
            'molecularWeights'
            'baseCounts'
            'lengths'
            'halfLives'
            'compartments'
            'meanNTPConcentration'
            'meanNDPConcentration'
            'meanNMPConcentration'
            };
        fittedConstantNames     = {   %names of process properties that are considered fitted constants, and should be stored with the simulation as such
            'biomassComposition'
            'biomassProduction'
            'byproducts'
            };
        stateNames              = {   %names of properties which are part of the simulation's state
            'counts'
            'processRequirements'
            'processAllocations'
            'processUsages'
            };
        dependentStateNames     = {}; %names of properties which can be calculated from the simulation's state
    end

    %indices
    properties
        hydrophobicIndexs    %indices of hydrophobic counts
        aminoAcidIndexs      %amino acids indices within counts
        dnmpIndexs           %dAMP, dCMP, dGMP, dTMP indices within counts
        dndpIndexs           %dADP, dCDP, dGDP, dTDP indices within counts
        dntpIndexs           %dATP, dCTP, dGTP, dTTP indices within counts
        nmpIndexs            %AMP, CMP, GMP, UMP indices within counts
        ndpIndexs            %ADP, CDP, GDP, UDP indices within counts
        ntpIndexs            %ATP, CTP, GTP, UTP indices within counts
        ampIndexs            %AMP index within counts
        adpIndexs            %ADP index within counts
        atpIndexs            %ATP index within counts
        phosphateIndexs      %Pi index within counts
        diphosphateIndexs    %PPi index within counts
        waterIndexs          %H2O index within counts
        hydrogenIndexs       %H index within counts
        dr5pIndexs           %DR5P index within counts
        unmodifiedBaseIndexs %AD, CSN, GN, THY indices within counts
        m6ADIndexs           %m6AD index within counts
        lipidIndexs          %lipid indices within counts
        polyamineIndexs      %polyamine indices within counts
        carbohydrateIndexs   %carbohydrate indices within counts
        vitaminIndexs        %vitamin indices within counts
        ionIndexs            %ion indices within counts
        antibioticIndexs     %antibiotic indices within counts
    end

    %constants
    properties
        setCounts                                     %experimentally set metabolic conditions (concentration, time when concentration is applied)

        experimentalBiomassComposition                %experimentally observed biomass composition (molecules / cell)
        experimentalBiomassCompositionWeightFractions %experimentally observed biomass composition (weight fraction)
        experimentalBiomassCompositionMolFractions    %experimentally observed biomass composition (mol fraction)
        biomassComposition                            %metabolite composition of biomass (mollecules / cell) (metabolites X compartments)
        biomassProduction                             %biomass production
        byproducts                                    %expected energy and other byproducts of networks (ADP, AMP, PPi, Pi, H)
        processBiomassProduction                      %expected metabolic demand of each process
        processByproduct                              %expected metabolic byproducts produced by each process

        meanNTPConcentration                          %mean ATP, GTP concentration (mM)
        meanNDPConcentration                          %mean ADP, GDP concentration (mM)
        meanNMPConcentration                          %mean AMP, CMP, GMP, UMP concentration (mM)

        processWholeCellModelIDs                      %IDs of processes
    end

    %constants
    properties (Dependent)
        metabolismProduction       %biomass production

        experimentalNMPComposition %experimentally determined mol fractions of NMPs
        experimentalAAComposition  %experimentally determined mol fractions of AAs

        nmpComposition             %mol fractions of NMPs
        aaComposition              %mol fractions of AAs

        wetWeight                  %water weight (g)
    end

    %state
    properties
        processRequirements        %instantaneous metabolic requirements of each process
        processAllocations         %instantaneous metabolites allocated to each process
        processUsages              %instantaneous metabolic usages of each process
    end

    %constructor
    methods
        function this = Metabolite(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.MoleculeCountState(wholeCellModelID, name);
        end
    end

    methods
        function initializeConstants(this, knowledgeBase, simulation)
            import edu.stanford.covert.cell.sim.constant.Condition;
            import edu.stanford.covert.util.ConstantUtil;

            this.initializeConstants@edu.stanford.covert.cell.sim.MoleculeCountState(knowledgeBase, simulation);

            %ids, names, weights
            this.wholeCellModelIDs = {knowledgeBase.metabolites.wholeCellModelID}';
            this.names             = {knowledgeBase.metabolites.name}';
            this.molecularWeights  = [knowledgeBase.metabolites.molecularWeight]';

            %indices
            this.hydrophobicIndexs    = find([knowledgeBase.metabolites.hydrophobic]);
            this.aminoAcidIndexs      = this.getIndexs({'ALA'; 'ARG'; 'ASN'; 'ASP'; 'CYS'; 'GLN'; 'GLU'; 'GLY'; 'HIS'; 'ILE'; 'LEU'; 'LYS'; 'MET'; 'PHE'; 'PRO'; 'SER'; 'THR'; 'TRP'; 'TYR'; 'VAL'; 'FMET'});
            this.dnmpIndexs           = this.getIndexs({'DAMP'; 'DCMP'; 'DGMP'; 'DTMP'});
            this.dndpIndexs           = this.getIndexs({'DADP'; 'DCDP'; 'DGDP'; 'DTDP'});
            this.dntpIndexs           = this.getIndexs({'DATP'; 'DCTP'; 'DGTP'; 'DTTP'});
            this.nmpIndexs            = this.getIndexs({'AMP';  'CMP';  'GMP';  'UMP'});
            this.ndpIndexs            = this.getIndexs({'ADP';  'CDP';  'GDP';  'UDP'});
            this.ntpIndexs            = this.getIndexs({'ATP';  'CTP';  'GTP';  'UTP'});
            this.ampIndexs            = this.getIndexs({'AMP'});
            this.adpIndexs            = this.getIndexs({'ADP'});
            this.atpIndexs            = this.getIndexs({'ATP'});
            this.phosphateIndexs      = this.getIndexs({'PI'});
            this.diphosphateIndexs    = this.getIndexs({'PPI'});
            this.waterIndexs          = this.getIndexs({'H2O'});
            this.hydrogenIndexs       = this.getIndexs({'H'});
            this.dr5pIndexs           = this.getIndexs({'DR5P'});
            this.unmodifiedBaseIndexs = this.getIndexs({'AD'; 'CSN'; 'GN'; 'THY'});
            this.m6ADIndexs           = this.getIndexs({'m6AD'});
            this.lipidIndexs          = find(strcmp({knowledgeBase.metabolites.category}', 'lipid'));
            this.polyamineIndexs      = find(strcmp({knowledgeBase.metabolites.category}', 'polyamine'));
            this.carbohydrateIndexs   = find(strcmp({knowledgeBase.metabolites.category}', 'carbohydrate'));
            this.vitaminIndexs        = find(strcmp({knowledgeBase.metabolites.category}', 'vitamin'));
            this.ionIndexs            = find(strcmp({knowledgeBase.metabolites.category}', 'ion'));
            this.antibioticIndexs     = find(strcmp({knowledgeBase.metabolites.category}', 'antibiotic'));

            %media
            this.setCounts = zeros(0, 6);
            for i = 1:knowledgeBase.numMetabolites
                if isempty(knowledgeBase.metabolites(i).mediaCompartments); continue; end;

                value = zeros(length(knowledgeBase.metabolites(i).mediaCompartments), 6);
                value(:, Condition.objectIndexs)      = knowledgeBase.metabolites(i).idx;
                value(:, Condition.compartmentIndexs) = [knowledgeBase.metabolites(i).mediaCompartments.idx];
                value(:, Condition.valueIndexs)       = knowledgeBase.metabolites(i).mediaConcentrations;
                value(:, Condition.initialTimeIndexs) = knowledgeBase.metabolites(i).mediaInitialTimes;
                value(:, Condition.finalTimeIndexs)   = knowledgeBase.metabolites(i).mediaFinalTimes;
                value(:, Condition.objectCompartmentIndexs) = sub2ind([knowledgeBase.numMetabolites knowledgeBase.numCompartments],...
                    value(:, Condition.objectIndexs),...
                    value(:, Condition.compartmentIndexs));

                this.setCounts = [this.setCounts; value];
            end

            this.setCounts(:, Condition.valueIndexs) = round(...
                this.setCounts(:, Condition.valueIndexs) * ...
                edu.stanford.covert.util.ConstantUtil.nAvogadro / 1000 * simulation.state('Geometry').chamberVolume);

            %biomass
            this.experimentalBiomassComposition = -knowledgeBase.biomassComposition;
            this.experimentalBiomassCompositionMolFractions = sum(this.experimentalBiomassComposition, 2) / sum(this.experimentalBiomassComposition(:));
            this.experimentalBiomassCompositionWeightFractions = ...
                sum(this.experimentalBiomassCompositionMolFractions, 2) .* this.molecularWeights / ...
                (sum(this.experimentalBiomassCompositionMolFractions, 2)' * this.molecularWeights);

            %processes
            this.processWholeCellModelIDs = {knowledgeBase.processes.wholeCellModelID}';
        end

        function allocateMemory(this, numTimePoints)
            import edu.stanford.covert.util.SparseMat;

            this.allocateMemory@edu.stanford.covert.cell.sim.MoleculeCountState(numTimePoints);

            nProcesses = numel(this.processWholeCellModelIDs);
            this.processRequirements = SparseMat([], [], [numel(this.wholeCellModelIDs) * this.compartment.count  nProcesses  numTimePoints]);
            this.processAllocations = SparseMat([], [], [numel(this.wholeCellModelIDs) * this.compartment.count  nProcesses  numTimePoints]);
            this.processUsages = SparseMat([], [], [numel(this.wholeCellModelIDs) * this.compartment.count  nProcesses  numTimePoints]);
        end
    end

    %getters
    methods
        function value = get.metabolismProduction(this)
            value = this.biomassProduction - this.byproducts;
        end

        %experimental composition
        function value = get.experimentalNMPComposition(this)
            value = this.experimentalBiomassCompositionMolFractions(this.nmpIndexs);
            value = value / sum(value);
        end

        function value = get.experimentalAAComposition(this)
            value = this.experimentalBiomassCompositionMolFractions(this.aminoAcidIndexs);
            value = value / sum(value);
        end

        %composition
        function value = get.nmpComposition(this)
            value = this.biomassComposition(this.nmpIndexs, :);
            value = value / sum(value(:));
        end

        function value = get.aaComposition(this)
            value = sum(this.biomassComposition(this.aminoAcidIndexs, :), 2);
            value = value / sum(value(:));
        end

        %weight
        function value = calcDryWeight(this)
            if size(this.counts, 3) == 1
                value = this.molecularWeights' * this.counts;
            else
                value = multiprod(this.molecularWeights', this.counts, [1 2], [1 2]);
            end
            value = ...
                + value / edu.stanford.covert.util.ConstantUtil.nAvogadro ...
                - this.wetWeight;
        end

        function value = get.wetWeight(this)
            value = this.counts(this.waterIndexs, :, :) * this.molecularWeights(this.waterIndexs) / edu.stanford.covert.util.ConstantUtil.nAvogadro;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% eeg_checkchanlocs() - Check the consistency of the channel locations structure
%                  of an EEGLAB dataset.
%
% Usage:
%  >> EEG = eeg_checkchanlocs( EEG, 'key1', value1, 'key2', value2, ... );
%  >> [chanlocs chaninfo] = eeg_checkchanlocs( chanlocs, chaninfo, 'key1', value1, 'key2', value2, ... );
%
% Inputs:
%   EEG      - EEG dataset
%   chanlocs - EEG.chanlocs structure
%   chaninfo - EEG.chaninfo structure
%
% Outputs:
%   EEG        - new EEGLAB dataset with updated channel location structures
%                EEG.chanlocs, EEG.urchanlocs, EEG.chaninfo
%   chanlocs   - updated channel location structure
%   chaninfo   - updated chaninfo structure
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, March 2, 2011

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

% Hey Arno -- this is a quick fix to make an analysis work for Makoto
% I think the old version had a bug...

function [chans, chaninfo, chanedit]= eeg_checkchanlocs(chans, chaninfo);

if nargin < 1
    help eeg_checkchanlocs;
    return;
end;

if nargin < 2
    chaninfo = [];
end;

processingEEGstruct = 0;
if isfield(chans, 'data')
    processingEEGstruct = 1;
    tmpEEG = chans;
    chans = tmpEEG.chanlocs;
    chaninfo = tmpEEG.chaninfo;
end;

if ~isfield(chans, 'datachan')
    [chanedit,dummy,complicated] = insertchans(chans, chaninfo);
else
    chanedit = chans;
    complicated = true;
end;

nosevals       = { '+X' '-X' '+Y' '-Y' };
if ~isfield(chaninfo, 'plotrad'), chaninfo.plotrad = []; end;
if ~isfield(chaninfo, 'shrink'),  chaninfo.shrink = [];  end;
if ~isfield(chaninfo, 'nosedir'), chaninfo.nosedir = nosevals{1}; end;

% handles deprecated fields
% -------------------------
plotrad  = [];
if isfield(chanedit, 'plotrad'),
    plotrad = chanedit(1).plotrad;
    chanedit = rmfield(chanedit, 'plotrad');
    if isstr(plotrad) & ~isempty(str2num(plotrad)), plotrad = str2num(plotrad); end;
    chaninfo.plotrad = plotrad;
end;
if isfield(chanedit, 'shrink') && ~isempty(chanedit(1).shrink)
    shrinkorskirt = 1;
    if ~isstr(chanedit(1).shrink)
        plotrad = 0.5/(1-chanedit(1).shrink); % convert old values
    end;
    chanedit = rmfield(chanedit, 'shrink');
    chaninfo.plotrad = plotrad;
end;

% set non-existent fields to []
% -----------------------------
fields    = { 'labels' 'theta' 'radius' 'X'   'Y'   'Z'   'sph_theta' 'sph_phi' 'sph_radius' 'type' 'ref' 'urchan' };
fieldtype = { 'str'    'num'   'num'    'num' 'num' 'num' 'num'       'num'     'num'        'str'  'str' 'num'    };
check_newfields = true; %length(fieldnames(chanedit)) < length(fields);
if ~isempty(chanedit)
    for index = 1:length(fields)
        if check_newfields && ~isfield(chanedit, fields{index})
            % new field
            % ---------
            if strcmpi(fieldtype{index}, 'num')
                chanedit = setfield(chanedit, {1}, fields{index}, []);
            else
                for indchan = 1:length(chanedit)
                    chanedit = setfield(chanedit, {indchan}, fields{index}, '');
                end;
            end;
        else
            % existing fields
            % ---------------
            allvals = {chanedit.(fields{index})};
            if strcmpi(fieldtype{index}, 'num')
                if ~all(cellfun('isclass',allvals,'double'))
                    numok = cellfun(@isnumeric, allvals);
                    if any(numok == 0)
                        for indConvert = find(numok == 0)
                            chanedit = setfield(chanedit, {indConvert}, fields{index}, []);
                        end;
                    end;
                end
            else
                strok = cellfun('isclass', allvals,'char');
                if strcmpi(fields{index}, 'labels'), prefix = 'E'; else prefix = ''; end;
                if any(strok == 0)
                    for indConvert = find(strok == 0)
                        try
                              strval   = [ prefix num2str(getfield(chanedit, {indConvert}, fields{index})) ];
                              chanedit = setfield(chanedit, {indConvert}, fields{index}, strval);
                        catch
                            chanedit = setfield(chanedit, {indConvert}, fields{index}, '');
                        end;
                    end;
                end;
            end;
        end;
    end;
end;
if ~isequal(fieldnames(chanedit)',fields)
    try
        chanedit = orderfields(chanedit, fields);
    catch, end;
end

% check if duplicate channel label
% --------------------------------
if isfield(chanedit, 'labels')
    tmp = sort({chanedit.labels});
    if any(strcmp(tmp(1:end-1),tmp(2:end)))
        disp('Warning: some channels have the same label');
    end
end;

% check for empty channel label
% -----------------------------
if isfield(chanedit, 'labels')
    indEmpty = find(cellfun(@isempty, {chanedit.labels}));
    if ~isempty(indEmpty)
        tmpWarning = warning('backtrace');
        warning backtrace off;
        warning('channel labels should not be empty, creating unique labels');
        warning(tmpWarning);
        for index = indEmpty
            chanedit(index).labels = sprintf('E%d', index);
        end;
    end;
end;

% remove fields
% -------------
if isfield(chanedit, 'sph_phi_besa'  ), chanedit = rmfield(chanedit, 'sph_phi_besa'); end;
if isfield(chanedit, 'sph_theta_besa'), chanedit = rmfield(chanedit, 'sph_theta_besa'); end;

% reconstruct the chans structure
% -------------------------------
if complicated
    [chans chaninfo.nodatchans] = getnodatchan( chanedit );
    if ~isfield(chaninfo, 'nodatchans'), chaninfo.nodatchans = []; end;
    if isempty(chanedit)
        for iField = 1:length(fields)
            chanedit = setfield(chanedit, fields{iField}, []);
        end;
    end;
else
    chans = rmfield(chanedit,'datachan');
    chaninfo.nodatchans = [];
end

if processingEEGstruct
    tmpEEG.chanlocs = chans;
    tmpEEG.chaninfo = chaninfo;
    chans = tmpEEG;
end;

% ---------------------------------------------
% separate data channels from non-data channels
% ---------------------------------------------
function [chans, fidsval] = getnodatchan(chans)
if isfield(chans,'datachan')
    [chans(cellfun('isempty',{chans.datachan})).datachan] = deal(0);
    fids = [chans.datachan] == 0;
    fidsval = chans(fids);
    chans = rmfield(chans(~fids),'datachan');
else
    fids = [];
end;

% ----------------------------------------
% fuse data channels and non-data channels
% ----------------------------------------
function [chans, chaninfo,complicated] = insertchans(chans, chaninfo, nchans)
if nargin < 3, nchans = length(chans); end;
[chans.datachan] = deal(1);
complicated = false;        % whether we need complicated treatment of datachans & co further down the road.....

if isfield(chans,'type')
    mask = strcmpi({chans.type},'FID') | strcmpi({chans.type},'IGNORE');
    if any(mask)
        [chans(mask).datachan] = deal(0);
        complicated = true;
    end
end
if length(chans) > nchans & nchans ~= 0 % reference at the end of the structure
    chans(end).datachan = 0;
    complicated = true;
end;
if isfield(chaninfo, 'nodatchans')
    if ~isempty(chaninfo.nodatchans) && isstruct(chaninfo.nodatchans)
        chanlen = length(chans);
        for index = 1:length(chaninfo.nodatchans)
            fields = fieldnames( chaninfo.nodatchans );
            ind = chanlen+index;
            for f = 1:length( fields )
                chans = setfield(chans, { ind }, fields{f}, getfield( chaninfo.nodatchans, { index },  fields{f}));
            end;
            chans(ind).datachan = 0;
            complicated = true;
        end;
        chaninfo = rmfield(chaninfo, 'nodatchans');

        % put these channels first
        % ------------------------
        % tmp = chans(chanlen+1:end);
        % chans(length(tmp)+1:end) = chans(1:end-length(tmp));
        % chans(1:length(tmp)) = tmp;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ basis, filt_sig, param ] = gsp_eigenspace_estimation( G, k, param )
%GSP_EIGENSPACE_ESTIMATION Estimation of first eigenvectors of any graph Laplacian
%   Usage:  basis = gsp_eigenspace_estimation(G,k);
%           [basis, approx_U, param] = gsp_eigenspace_estimation(G,k,param);
%
%   Input parameters :
%         G          : Graph structure.
%         k          : Dimension of the subspace.
%         param      : Optional parameters
%   Output parameters:
%         basis      : Approximated basis of k first eigenvectors
%         filt_sig   : Filtered random signal
%         param      : Optional parameters (with new entries)
%
%   'gsp_eigenspace_estimation(G,k)' computes an estimation of the first
%   $k$ eigenvectors of the Laplacian of G using Gaussian random signal
%   filtering, following the FEARS method described in paratte2017fast.
%
%
%   Example:::
%
%         G = gsp_sensor(256);
%         G = gsp_estimate_lmax(G);
%         k = 8;
%         param.order = 100;
%         Uk_est = gsp_eigenspace_estimation(G, k, param);
%         G = gsp_compute_fourier_basis(G);
%         proj_energy = norm(Uk_est' * G.U(:, 1:k), 'fro');
%
%
%   Additional parameters
%   ---------------------
%
%   * *param.filter*  : Select the filter to be used for the computation.
%     * 'lp-ch'   : Chebyshev polynomial approximation
%     * 'lp-jch'  : Jackson-Chebyshev polynomial approximation
%     * 'expwin'  : Exponentially decreasing polynomial approximation Default: 'lp-jch'
%   * *param.order* : Degree of the polynomial approximation (default=50).
%   * *param.lk_est_method* : Select the version of lk estimation.
%     * 'fast'      : Accelerated method using local uniformity assumption
%     * 'std'  : Usual method using dichotomy all the time Default: 'fast'
%   * *param.R* : Random matrix to use (of size N > d, d >= k)
%     (default: Gaussian(0, 1/k) of size Nxk)
%   * *param.pcoef* : Polynomial coefficients if already known.
%   * *param.lk* : Estimated value of lambda_k if already known.
%   * *param.verbose* : Verbosity level (0 no log - 1 display warnings) (default 1).
%
%
%   References: paratte2016fast
%

% Author: Johan Paratte, Lionel Martin
% Date: 3 November 2016

if nargin < 3, param = struct; end
if ~isfield(param, 'filter'), param.filter = 'lp-jch'; end
if ~isfield(param, 'order'), param.order = 50; end
if ~isfield(param, 'lk_est_method'), param.lk_est_method = 'fast'; end
if ~isfield(param, 'R'), param.R = randn(G.N, k)/sqrt(k); end
if ~isfield(param, 'verbose'), param.verbose = 1; end

assert(size(param.R, 1) == G.N && size(param.R, 2) >= k, 'The optional parameter R has wrong size.');

if ~isfield(param, 'pcoefs')
    if param.verbose, disp('Polynomial filtering required. Computing polynomial coefficients...'); end

    if ~isfield(G, 'lmax')
        G = gsp_estimate_lmax(G);
        warning(['GSP_EIGENSPACE_ESTIMATION: The variable lmax is not available.', ...
            'The function will compute it for you. However, if you apply ', ...
            'many time this function, you should precompute it using the ', ...
            'function: gsp_estimate_lmax.']);
    end

    if ~isfield(param, 'lk')
        if param.verbose, fprintf('Estimation of lambda_k');
        end

        tic;
        switch param.lk_est_method
            case 'fast'
                [param.lk, info] = gsp_fast_estimate_lk(G, k, param);
                if param.verbose, disp('using our accelerated method.'); end

            case 'std'
                [~, param.lk, ~, ~, nb_iter_lk, k_est_lk] = gsp_estimate_lk(G, k, param);
                if param.verbose, disp('using the standard method.'); end

            otherwise
                error('Unknown method for lk_est_method.');
        end
        t = toc;

        if param.verbose
            fprintf(['* Estimated lk: %d\n', ...
            '* Time to estimate lk: %f sec\n', ...
            '* in %d iterations with k_est=%d (target=%d)\n'], ...
            param.lk, t, mean(info.calls), mean(info.k_est), k);
        end
    else
        warning('lambda_k was provided to the method from param.');
    end

    tic;
    switch param.filter
        case 'lp-ch'
            [param.pcoefs, ~] = jackson_cheby_poly_coefficients(0, param.lk, [0, G.lmax], param.order);

        case 'lp-jch'
            [~, param.pcoefs] = jackson_cheby_poly_coefficients(0, param.lk, [0, G.lmax], param.order);

        case 'expwin'
            ew = gsp_design_expwin(G, param.lk/G.lmax);
            param.pcoefs = gsp_cheby_coeff(G, ew, param.order);

        otherwise
            error('Unknown filter type!');
    end

    t = toc;
    if param.verbose, fprintf('* Time to compute polynomial coefficients: %f sec.\n', t); end

else
    if param.verbose, warning('pcoef was provided to the method from param.'); end
end

if param.verbose, disp('Filtering random signals...'); end
tic;
filt_sig = gsp_cheby_op(G, param.pcoefs, param.R);
t = toc;
if param.verbose, fprintf('* Time to filter random signals: %f sec.\n', t); end

if param.verbose, disp('Computing the SVD for eigenspace recovery...'); end
tic;
[basis, ~, ~] = svd(filt_sig, 'econ');
t = toc;
if param.verbose, fprintf('* Time to compute SVD: %f sec.\n', t); end

end

repo: gspbox


text:
function [DVH,volume_VOI_first,volume_VOI_check] = dicomrt_dvhcal(VOI,study,dose_xmesh,dose_ymesh,dose_zmesh,dvhselect)
% dicomrt_dvhcal(VOI,study,dose_xmesh,dose_ymesh,dose_zmesh,dvhselect)
%
% Calculate DOSE-VOLUME-HISTOGRAMS for VOIs and cell_case_study.
% This function calculate DVHs for all the available VOIs or for a selected number of them.
%
% study is the RTPLAN or the MC generated dataset containing the dose matrix
% ct is the CT dataset
% dose_xmesh,dose_ymesh,dose_zmesh are x-y-z coordinates of the center of the dose-voxels
% VOIs is a cell array and contains all the Volumes of Interest
% dvhselect is an OPTIONAL vector which contains the number of VOI to calculate DVHs for.
%
% NOTE: if dvhselect is omitted DVHs are calculated for all the VOIs
%
% DVHs are stored in a cell array with the following structure:
%
%  -----------------------------
%  | [DVH 1] | [3D dose mask]  |
%  |         -------------------
%  |         | [dvh data]      |
%  |         -------------------
%  |         | VOI volume      |
%  |         -------------------
%  |         | Voxel volume    |
%  -----------------------------
%  |   ...   |     ...         |
%  -----------------------------
%  | [DVH n] | [3D dose mask]  |
%  |         -------------------
%  |         | [dvh data]      |
%  |         -------------------
%  |         | VOI volume      |
%  |         -------------------
%  |         | Voxel volume    |
%  -----------------------------
%
% [dvh data] is a 2 columns vector with following structure:
%
% -----------------
% | dose | volume |
% | (Gy) |  (cc)  |
% -----------------
% |      |        |
% |      |        |
% |      |        |
% |      |        |
% -----------------
%
%
% Example:
%
% dvh=dicomrt_dvhcal(VOI,dose,dose_xmesh,dose_ymesh,dose_zmesh) returns
% the dvhs for the DVH for all the VOIs in VOI and the dose distribution in the
% 3D matrix 'dose'.
%
% See also dicomrt_loaddose, roifilt2, roipoly, dicomrt_dvhplot
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check case and set-up some parameters and variables
[dose_temp,type_dose,doselabel,PatientPosition]=dicomrt_checkinput(study);
dose=dicomrt_varfilter(dose_temp);
[VOI_temp]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

% Define cell array
DVH=cell(size(VOI,1),2);

% Write VOIs labels
for k=1:size(VOI,1)
    DVH{k,1}=VOI{k,1};
end

% Retrieve dmax
dmax=max(max(max(dose)));

% Build dose grid
dosegrid=zeros(double(int32(dmax)),1);
%for i=1:double(int32(dmax))                    % use with hist
%    x(i)=dmax/double(int32(dmax))*(0.5+i-1);   % hist uses centers
%end                                            %
for i=1:double(int32(dmax))                     % use with histc
    dosegrid(i)=i*dmax/double(int32(dmax));     % histc uses edges
end                                             %

% Calculate voxel volume
[voxelvolume, status]=dicomrt_voxelvolumecal(study,dose_xmesh,dose_ymesh,dose_zmesh);

volume_VOI_check=[];
volume_VOI_first=[];

for k=1:length(dvhselect) % loop over VOIs
    % Print header
    disp(['DVH calculation for VOI: ',VOI{dvhselect(k),1}]);
    % Retrieve info necessary to get volume_VOI
    [mask_VOI,volume_VOI,mask4VOI]=dicomrt_mask(VOI_temp,dose_temp,dose_xmesh,dose_ymesh,dose_zmesh,...
        dvhselect(k),'nan','y'); % dose matrix masked in dicomrt_Vlevel
    % Build frequency DVH data
    % case 1: uniform slice thickness
    volume_VOI_first=[volume_VOI_first,volume_VOI];
    if status==0
        vgrid=histc(reshape(mask_VOI{2,1},1,size(mask_VOI{2,1},1)*size(mask_VOI{2,1},2)*size(mask_VOI{2,1},3)),dosegrid);
        % get back voxels with dose<min(dosegrid) cause histc does not
        % account for it
        temp_min=find(mask_VOI{2,1}<dosegrid(1) & mask_VOI{2,1}~=0);
        %temp_min2=find(mask_VOI{2,1}<dosegrid(1));
        if isempty(temp_min)~=1
            vgrid(1)=vgrid(1)+length(temp_min);
        end
        vgrid=vgrid.*voxelvolume(1);
        temp2=cumsum(vgrid);
        volume_VOI_check=[volume_VOI_check,temp2(end)];
    else % case 2: non-uniform slice thickness
        vgrid=zeros(length(dosegrid),1);
        for i=1:size(mask_VOI{2,1},3)
            temp_vgrid=histc(reshape(mask_VOI{2,1}(:,:,i),1,size(mask_VOI{2,1}(:,:,i),1)*size(mask_VOI{2,1}(:,:,i),2)),dosegrid);
            % get back voxels with dose<min(dosegrid) cause histc does not
            % account for it
            temp_min=find(mask_VOI{2,1}(:,:,i)<dosegrid(1) & mask_VOI{2,1}~=0);
            if isempty(temp_min)~=1
                temp_vgrid(1)=temp_vgrid(1)+length(temp_min);
            end
            vgrid=vgrid+temp_vgrid'.*voxelvolume(i);
        end
        temp2=cumsum(vgrid);
        volume_VOI_check=[volume_VOI_check,temp2(end)];
        %for i=1:length(dosegrid)
        %    if i==1
        %        vgrid(i)=dicomrt_Vlevel(mask_VOI,0,dose_xmesh,dose_ymesh,dose_zmesh,...
        %            [[0 dosegrid(i)] 0 0 0],VOI_temp,dvhselect(k),1,volume_VOI);
        %    else
        %        vgrid(i)=dicomrt_Vlevel(mask_VOI,0,dose_xmesh,dose_ymesh,dose_zmesh,...
        %            [[dosegrid(i-1) dosegrid(i)] 0 0 0],VOI_temp,dvhselect(k),1,volume_VOI);
        %    end
        %end
    end
    vgrid=dicomrt_makevertical(vgrid);
    % Store DVH data
    DVH{dvhselect(k),2}{1,1}=mask_VOI;
    DVH{dvhselect(k),2}{2,1}=[dosegrid,vgrid];
    DVH{dvhselect(k),2}{3,1}=volume_VOI;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function H = ent_g(x, biascorrect)
% ENT_G Entropy of a Gaussian variable in bits
%   H = ent_g(x) returns the entropy of a (possibly
%   multidimensional) Gaussian variable x with bias correction.
%   Rows of x correspond to samples, columns to dimensions/variables.
%   (Samples first axis)
%
%   biascorrect : true / false option (default true) which specifies
%   whether bias correction should be applied to the estimated entropy.

if isvector(x)
    x = x(:);
end
if ndims(x)~=2
    error('ent_g: input arrays should be 2d')
end
[Ntrl, Nvarx] = size(x);

if nargin < 2
    % default is to apply bias correction
    biascorrect = true;
end

% demean data
gx.m = sum(x,1)/Ntrl;
x = bsxfun(@minus,x,gx.m);

% covariance
C = (x'*x) / (Ntrl - 1);
chC = chol(C);

% entropy in nats
HX = sum(log(diag(chC))) + 0.5*Nvarx*(log(2*pi)+1);

ln2 = log(2);
if biascorrect
    psiterms = psi((Ntrl - (1:Nvarx))/2) / 2;
    dterm = (ln2 - log(Ntrl-1)) / 2;
    HX = (HX - Nvarx*dterm - sum(psiterms));
end

% convert to bits
H = HX / ln2;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Generates test fixtures for whole cell simulation.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 5/31/2011
function generateTestFixtures()
import edu.stanford.covert.cell.sim.runners.SimulationRunner;

%create runner which helps do the tasks below
runner = SimulationRunner(...
    'useCachedKb', false, ...
    'useCachedSim', false, ...
    'cacheKb', true, ...
    'cacheSimulation', true);

%construction KB and simulation
sim = runner.constructKbAndSimulation();

% verify initial growth rate distribution
if ~runner.isSimulationFitted(sim)
    throw(MException('Simulation:notFitted', 'Simulation not properly fitted'));
end

%save simulation parameters in JSON format
runner.cacheSimulationParametersJson(sim);

%generate simulation test fixtures
runner.cacheSimulationTestFixtures(sim);


text:
function tree = copy(tree,subuid,uid)
% XMLTREE/COPY Copy Method (copy a subtree in another branch)
% FORMAT tree = copy(tree,subuid,uid)
%
% tree      - XMLTree object
% subuid    - UID of the subtree to copy
% uid       - UID of the element where the subtree must be duplicated
%__________________________________________________________________________
%
% Copy a subtree to another branch.
% The tree parameter must be in input AND in output.
%__________________________________________________________________________
% Copyright (C) 2002-2015  http://www.artefact.tk/

% Guillaume Flandin
% $Id: copy.m 6480 2015-06-13 01:08:30Z guillaume $

%error(nargchk(2,3,nargin));

if nargin == 2
    uid = parent(tree,subuid);
end

l = length(tree);
tree = sub_copy(tree,subuid,uid);
tree.tree{uid}.contents = [tree.tree{uid}.contents l+1];

% to have the copy next to the original and not at the end?
%  contents = get(tree,parent,'contents');
%  i = find(contents==uid);
%  tree = set(tree,parent,'contents',[contents(1:i) l+1 contents(i+1:end)]);

%==========================================================================
function tree = sub_copy(tree,uid,p)

    l = length(tree);
    tree.tree{l+1} = tree.tree{uid};
    tree.tree{l+1}.uid = l+1;
    tree.tree{l+1}.parent = p;
    tree.tree{l+1}.contents = [];
    if isfield(tree.tree{uid},'contents')
        contents = get(tree,uid,'contents');
        m = length(tree);
        for i=1:length(contents)
            tree.tree{l+1}.contents = [tree.tree{l+1}.contents m+1];
            tree = sub_copy(tree,contents(i),l+1);
            m = length(tree);
        end
    end


text:
function errors = test_wiener()

errors = 0;
gsp_reset_seed
N = 20;
M = 3;
sigma = 0.3;

G = gsp_sensor(N);
G = gsp_compute_fourier_basis(G);
g = @(x) sin(x);
psd = @(x) g(x).^2;

Mask = rand(N,1)>0.5;

x = gsp_filter_analysis(G,g,randn(N,M));

Mop =@(x) bsxfun(@times,Mask,x);

y1 = Mop(x);

y2 = Mop(x+sigma*randn(N,M));
clear paramopt
paramopt.maxit = 1000;
paramopt.tol = 1e-12;
sol11 = gsp_wiener_inpainting(G,y1,Mask,psd,0,paramopt);
sol12 = gsp_wiener_inpainting_exact(G,y1,Mask,psd,0);

if norm(sol11-sol12,'fro')/norm(sol11,'fro') > 1e-10
    errors = errors + 1;
    norm(sol11-sol12,'fro')/norm(sol11,'fro')
    warning('Test Wiener opt 1 - ERROR');
else
    disp('Test Wiener opt 1 - OK');
end

sol21 = gsp_wiener_inpainting(G,y2,Mask,psd,sigma.^2,paramopt);
sol22 = gsp_wiener_inpainting_exact(G,y2,Mask,psd,sigma.^2);

if norm(sol21-sol22,'fro')/norm(sol22,'fro')> 1e-10
    errors = errors + 1;
    norm(sol21-sol22,'fro')/norm(sol22,'fro')
    warning('Test Wiener opt 2 - ERROR');
else
    disp('Test Wiener opt 2 - OK');
end

%%

g = gsp_design_expwin(G,0.2);
psd = @(x) g(x).^2;

Mask = rand(N,1)>0.5;
x = gsp_filter_analysis(G,g,randn(N,M));

Mop =@(x) bsxfun(@times,Mask,x);

y1 = Mop(x);

y2 = Mop(x+sigma*randn(N,M));

paramopt.tol = 1e-12;

%%
paramopt.maxit = 1000;
paramopt.gamma = 0.01;
sol11 = gsp_wiener_inpainting(G,y1,Mask,psd,0,paramopt);
sol12 = gsp_wiener_inpainting_exact(G,y1,Mask,psd,0);

if norm(sol11-sol12,'fro')/norm(sol11,'fro') > 1e-3
    errors = errors + 1;
    norm(sol11-sol12,'fro')/norm(sol11,'fro')
    warning('Test Wiener opt 3 - ERROR');
else
    disp('Test Wiener opt 3 - OK');
end

if norm(x-sol12,'fro')/norm(x,'fro') > 1e-10
    errors = errors + 1;
    norm(x-sol12,'fro')/norm(x,'fro')
    warning('Test Wiener opt 4 - ERROR');
else
    disp('Test Wiener opt 4 - OK');
end

%%
% figure(1)
% subplot(121)
% gsp_plot_signal_spectral(G,gsp_gft(G,x(:,1)))
% subplot(122)
% gsp_plot_signal_spectral(G,gsp_gft(G,sol12(:,1)))

%%
paramopt.gamma = 0.1;

paramopt.maxit = 1000;

sol21 = gsp_wiener_inpainting(G,y2,Mask,psd,sigma.^2,paramopt);
sol22 = gsp_wiener_inpainting_exact(G,y2,Mask,psd,sigma.^2);

if norm(sol21-sol22,'fro')/norm(sol22,'fro')> 1e-10
    errors = errors + 1;
    norm(sol21-sol22,'fro')/norm(sol22,'fro')
    warning('Test Wiener opt 5 - ERROR');
else
    disp('Test Wiener opt 5 - OK');
end

end


text:
function [varargout] = ft_plot_patch(hdat, vdat, varargin)

% FT_PLOT_PATCH plot a colored shape, similar to the MATLAB patch() function. It is
% similar in usage as ft_plot_vector, and they can be combined, for example,
% to plot an area equivalent to a SEM or STD-DEV around a line.
%
% Use as
%   ft_plot_patch(X, Y, ...)
% where X and Y are similar as the input to the MATLAB patch() function.
%
% Optional arguments should come in key-value pairs and can include
%   'axis'            = draw the local axis,  can be 'yes', 'no', 'xy', 'x' or 'y'
%   'box'             = draw a box around the local axes, can be 'yes' or 'no'
%   'tag'             = string, the name assigned to the object. All tags with the same name can be deleted in a figure, without deleting other parts of the figure.
%   'facecolor'       = see MATLAB standard patch properties
%   'facealpha'       = see MATLAB standard patch properties (note, approx. transparency can be achieved using 'facecolor')
%   'edgecolor'       = see MATLAB standard patch properties (default is 'none') (equivalent to 'linecolor' in PLOT)
%   'linestyle'       = see MATLAB standard patch properties
%   'linewidth'       = see MATLAB standard patch properties
%
% The color of the patchand the edges (i.e. border lines) can be specified in a variety of ways
%   - as a string with one character per line that you want to plot. Supported colors are the same as in PATCH, i.e. 'bgrcmykw'.
%   - as an 'RGB triplet', a 1x3 vector with values between 0 and 1
%   - as 'none' if you do not want the face of the patch to be filled (useful when you want to plot an empty box).
%
% It is possible to plot the object in a local pseudo-axis (c.f. subplot), which is specfied as follows
%   'hpos'            = horizontal position of the center of the local axes
%   'vpos'            = vertical position of the center of the local axes
%   'width'           = width of the local axes
%   'height'          = height of the local axes
%   'hlim'            = horizontal scaling limits within the local axes
%   'vlim'            = vertical scaling limits within the local axes
%
% Example
%   hdat = [1:10 10:-1:1];
%   vdat = rand(1,10);
%   vdat = [vdat vdat(end:-1:1)+1];
%   ft_plot_patch(hdat, vdat)
%
% See also FT_PLOT_VECTOR, PATCH, PLOT

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

% get the optional input arguments
hpos            = ft_getopt(varargin, 'hpos');
vpos            = ft_getopt(varargin, 'vpos');
width           = ft_getopt(varargin, 'width');
height          = ft_getopt(varargin, 'height');
hlim            = ft_getopt(varargin, 'hlim', 'maxmin');
vlim            = ft_getopt(varargin, 'vlim', 'maxmin');
axis            = ft_getopt(varargin, 'axis', false);
box             = ft_getopt(varargin, 'box', false);
tag             = ft_getopt(varargin, 'tag', '');
parent          = ft_getopt(varargin, 'parent', []);
facecolor       = ft_getopt(varargin, 'facecolor', 'b');
facealpha       = ft_getopt(varargin, 'facealpha', 1);
edgecolor       = ft_getopt(varargin, 'edgecolor', 'none');
linestyle       = ft_getopt(varargin, 'linestyle', 'none');
linewidth       = ft_getopt(varargin, 'linewidth', .5);

% convert the yes/no strings into boolean values
box = istrue(box);

% color management
if ischar(facecolor) && exist([facecolor '.m'], 'file')
  facecolor = feval(facecolor);
end
if ischar(edgecolor) && exist([edgecolor '.m'], 'file')
  edgecolor = feval(edgecolor);
end

% this should be a string, because valid options include yes, no, xy, x, y
if isequal(axis, true)
  axis = 'yes';
elseif isequal(axis, false)
  axis = 'no';
end

% everything is added to the current figure
holdflag = ishold;
if ~holdflag
  hold on
end

if ischar(hlim)
  switch hlim
    case 'maxmin'
      hlim = [min(hdat) max(hdat)];
    case 'maxabs'
      hlim = max(abs(hdat));
      hlim = [-hlim hlim];
    otherwise
      ft_error('unsupported option for hlim')
  end % switch
end % if ischar

if ischar(vlim)
  switch vlim
    case 'maxmin'
      vlim = [min(vdat(:)) max(vdat(:))];
    case 'maxabs'
      vlim = max(abs(vdat(:)));
      vlim = [-vlim vlim];
    otherwise
      ft_error('unsupported option for vlim')
  end % switch
end % if ischar

if vlim(1)==vlim(2)
  % vertical scaling cannot be determined, behave consistent to the plot() function
  vlim = [-1 1];
end

% these must be floating point values and not integers, otherwise the scaling fails
hdat = double(hdat);
vdat = double(vdat);
hlim = double(hlim);
vlim = double(vlim);

if isempty(hpos) && ~isempty(hlim)
  hpos = (hlim(1)+hlim(2))/2;
end
if isempty(vpos) && ~isempty(vlim)
  vpos = (vlim(1)+vlim(2))/2;
end

if isempty(width) && ~isempty(hlim)
  width = hlim(2)-hlim(1);
end

if isempty(height) && ~isempty(vlim)
  height = vlim(2)-vlim(1);
end

% first shift the horizontal axis to zero
if any(hlim) ~= 0
  hdat = hdat - (hlim(1)+hlim(2))/2;
  % then scale to length 1
  if hlim(2)-hlim(1)~=0
    hdat = hdat ./ (hlim(2)-hlim(1));
  else
    hdat = hdat /hlim(1);
  end
  % then scale to the new width
  hdat = hdat .* width;
end
% then shift to the new horizontal position
hdat = hdat + hpos;

if any(vlim) ~= 0
  % first shift the vertical axis to zero
  vdat = vdat - (vlim(1)+vlim(2))/2;
  % then scale to length 1
  vdat = vdat / (vlim(2)-vlim(1));
  % then scale to the new width
  vdat = vdat .* height;
end
% then shift to the new vertical position
vdat = vdat + vpos;

% plot the patch
h = patch(hdat, vdat, facecolor, 'facealpha', facealpha, 'edgecolor', edgecolor, 'linestyle', linestyle, 'linewidth', linewidth);

if box
  % this plots a box around the original hpos/vpos with appropriate width/height
  x1 = hpos - width/2;
  x2 = hpos + width/2;
  y1 = vpos - height/2;
  y2 = vpos + height/2;

  X = [x1 x2 x2 x1 x1];
  Y = [y1 y1 y2 y2 y1];
  h = line(X, Y);
  set(h, 'color', 'k');

  % this plots a box around the original hpos/vpos with appropriate width/height
  % boxposition = zeros(1,4);
  % boxposition(1) = hpos - width/2;
  % boxposition(2) = hpos + width/2;
  % boxposition(3) = vpos - height/2;
  % boxposition(4) = vpos + height/2;
  % ft_plot_box(boxposition, 'facecolor', 'none', 'edgecolor', 'k');

  % this plots a box around the complete data
  % boxposition = zeros(1,4);
  % boxposition(1) = hlim(1);
  % boxposition(2) = hlim(2);
  % boxposition(3) = vlim(1);
  % boxposition(4) = vlim(2);
  % ft_plot_box(boxposition, 'hpos', hpos, 'vpos', vpos, 'width', width, 'height', height, 'hlim', hlim, 'vlim', vlim);
end

if ~isempty(axis) && ~strcmp(axis, 'no')
  switch axis
    case {'yes' 'xy'}
      xaxis = true;
      yaxis = true;
    case {'x'}
      xaxis = true;
      yaxis = false;
    case {'y'}
      xaxis = false;
      yaxis = true;
    otherwise
      ft_error('invalid specification of the "axis" option')
  end

  if xaxis
    % x-axis should touch 0,0
    xrange = hlim;
    if sign(xrange(1))==sign(xrange(2))
      [dum, minind] = min(abs(hlim));
      xrange(minind) = 0;
    end
    ft_plot_line(xrange, [0 0], 'hpos', hpos, 'vpos', vpos, 'hlim', hlim, 'vlim', vlim, 'width', width, 'height', height);
  end
  if yaxis
    % y-axis should touch 0,0
    yrange = vlim;
    if sign(yrange(1))==sign(yrange(2))
      [dum, minind] = min(abs(vlim));
      yrange(minind) = 0;
    end
    ft_plot_line([0 0], yrange, 'hpos', hpos, 'vpos', vpos, 'hlim', hlim, 'vlim', vlim, 'width', width, 'height', height);
  end
end

set(h, 'tag', tag);

if ~isempty(parent)
  set(h, 'Parent', parent);
end

% the (optional) output is the handle
if nargout == 1
  varargout{1} = h;
end

if ~holdflag
  hold off
end


text:
function [ly] = ft2ly(ft)
% Convert length from feet to light years.
% Chad A. Greene 2012
ly = ft*3.221738542107e-17;


text:
function y = perform_dct2_transform(x,dir,dim)

% perform_dct2_transform - perform the DCT2 transform
%
%   y = perform_dct2_transform(x);
%
%   Copyright (c) 2007

if nargin<2
    dir=1;
end
if nargin<3
    dim = 1;
    if size(x,1)>1 & size(x,2)>1
        dim = 2;
    end
end

n = size(x,1);
if dim==1
    % perform 1D DCT
    p = prod(size(x))/n;
    if 0
        %% implementation with size x 4 FFT
        y = zeros(4*n,p);
        y(2:2:2*n,:) = x(:,:);
        y(4*n:-2:2*n+1,:) = x(:,:);
        y = fft(y,[],1);
        y = real( y(1:n,:) ) / sqrt(2*n);
        y(1,:) = y(1,:)/sqrt(2);
    else
        %% implementation with size x 2 FFT
        w = exp( -1i*pi/(2*n)*(0:n-1)');
        y = cat(1, x(:,:),zeros(n,p) );
        y = fft(y, [], 1);
        y = y(1:n,:);
        y = sqrt(2/n) * real( y.*repmat( w, [1 p] ) );
        y(1,:) = y(1,:)/sqrt(2);
    end
    y = reshape(y,size(x));
    return;
end

%% apply on 1st dimension
s = size(x); x = x(:,:);
y = perform_dct2_transform(x,dir,1);
y = permute(y,[2 1 3:10]);
y = perform_dct2_transform(y,dir,1);
y = permute(y,[2 1 3:10]);

return;

if size(x,2)==1 || dim==1
    % 1D transform
    y = zeros(4*n,size(x,2),size(x,3),size(x,4));
    y(2:2:2*n,:,:,:) = x; y(4*n:-2:2*n+1,:,:,:) = x;
    y = fft(y,[],1);
    y = real( y(1:n,:,:,:) ) / sqrt(2*n);
    y(:,1,:,:) = y(:,1,:,:)/sqrt(2);
else
    if dir==1
        if  0
        y = x;
        y(1:end/2,:,:,:) = x(1:2:end,:,:,:);
        y(end:-1:end/2+1,:,:,:) = x(2:2:end,:,:,:);
        y = fft(y,[],1);
        y = real( y.*repmat( exp(-1i*pi/(2*n)*(0:n-1)'), [1,size(x,2),size(x,3),size(x,4)]  ) );
        end

        % transform on x
        y = zeros(4*n,size(x,2),size(x,3),size(x,4));
        y(2:2:2*n,:,:,:) = x; y(4*n:-2:2*n+1,:,:,:) = x;
        y = fft(y,[],1);
        y = real( y(1:n,:,:,:) );
        % transform on y
        z = zeros(size(x,1),4*n,size(x,3),size(x,4));
        z(:,2:2:2*n,:,:) = y; z(:,4*n:-2:2*n+1,:,:) = y;
        z = fft(z,[],2);
        y = real( z(:,1:n,:,:) ) / (2*n);
        % y(:,1,:,:) = y(:,1,:,:)/sqrt(2);
        % y(1,:,:,:) = y(1,:,:,:)/sqrt(2);
    else
        error('Not yet implemented (DCT3)');

    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function area = triangle_area_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA_2D computes the area of a triangle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real AREA, the absolute area of the triangle.
%
  area = 0.5 * abs ( ...
      t(1,1) * ( t(2,2) - t(2,3) ) ...
    + t(1,2) * ( t(2,3) - t(2,1) ) ...
    + t(1,3) * ( t(2,1) - t(2,2) ) );

  return
end


text:
function index = end(X,position,numindices)
%END (overloaded)

if nargin == 2
    index = X.dim(position);
else
    if numindices == 1
        index = prod(X.dim);
    elseif numindices == length(X.dim)
        index = X.dim(position);
    else
        error('Indexing logic not supported yet');
    end
end


text:
% pop_rejchanspec() - reject artifacts channels in an EEG dataset using
%                  channel spectrum. The average spectrum for all selected
%                  is computed and a threshold is applied.
%
% Usage:
%   >> pop_rejchanspec( INEEG ) % pop-up interative window mode
%   >> [OUTEEG, indelec] = pop_rejchanspec( INEEG, 'key', 'val');
%
% Inputs:
%   INEEG       - input EEGLAB dataset
%
% Optional inputs:
%   'freqlims'  - [min max] frequency limits. May also be an array where
%                 each row defines a different set of limits. Default is
%                 35 to the Niquist frequency of the data.
%   'stdthresh' - [max] positive threshold in terms of standard deviation.
%                 Default is 5.
%   'absthresh' - [max] positive threshold in terms of spectrum units
%                 (overides the option above).
%   'averef'    - ['on'|'off'] 'on' computes average reference before
%                 applying threshold. Default is 'off'.
%   'plothist'  - ['on'|'off'] 'on' plot the histogram of values along
%                 with the threshold.
%   'plotchans'  - ['on'|'off'] 'on' plot the channels scrollplot with
%                 selected channels for rejection in red. Allow selected
%                 channels rejection with the 'REJECT' button.
%   'elec'      - [integer array] only include specific channels. Default
%                 is to use all channels.
%   'specdata'  - [fload array] use this array containing the precomputed
%                 spectrum instead of computing the spectrum. Default is
%                 empty.
%   'specfreqs' - [fload array] frequency array for precomputed spectrum
%                 above.
%   'verbose'   - ['on'|'off'] display information. Default is 'off'.
%
% Outputs:
%   OUTEEG    - output dataset with updated joint probability array
%   indelec   - indices of rejected electrodes
%   specdata  - data spectrum for the selected channels
%   specfreqs - frequency array for spectrum above
%
% Author: Arnaud Delorme, CERCO, UPS/CNRS, 2008-

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

function [EEG allrmchan specdata specfreqs com] = pop_rejchanspec(EEG, varargin)

if nargin < 1
    help pop_rejchanspec;
    return;
end;
allrmchan = [];
specdata  = [];
specfreqs = [];
com       = '';
if nargin < 2
    uilist = { { 'style' 'text' 'string' 'Electrode (number(s); Ex: 2 4 5)' } ...
               { 'style' 'edit' 'string' ['1:' int2str(EEG.nbchan)] } ...
               { 'style' 'text' 'string' 'Frequency limits [min max]' } ...
               { 'style' 'edit' 'string' [ '35 ' int2str(floor(EEG.srate/2)) ] } ...
               { 'style' 'text' 'string' 'Standard dev. threshold limits [max]' } ...
               { 'style' 'edit' 'string' '5' } ...
               { 'style' 'text' 'string' 'OR absolute threshold limit [min max]' } ...
               { 'style' 'edit' 'string' '' } ...
               { 'style' 'text' 'string' 'Compute average reference first (check=on)' } ...
               { 'style' 'checkbox' 'string' '' 'value' 0 } { } ...
               { 'style' 'text' 'string' 'Plot histogram of power values (check=on)' } ...
               { 'style' 'checkbox' 'string' '' 'value' 0 } { } ...
               { 'style' 'text' 'string' 'Plot channels scrollplot (check=on)' } ...
               { 'style' 'checkbox' 'string' '' 'value' 0 } { } ...
             };

    geom = { [2 1] [2 1] [2 1] [2 1] [2 0.3 0.7] [2 0.3 0.7] [2 0.3 0.7] };
    result = inputgui( 'uilist', uilist, 'geometry', geom, 'title', 'Reject channel using spectrum -- pop_rejchanspec()', ...
        'helpcom', 'pophelp(''pop_rejchan'')');
    if isempty(result), return; end;

    options = { 'elec' eval( [ '[' result{1} ']' ] ) 'stdthresh' str2num(result{3}) 'freqlims' str2num(result{2}) };
    if ~isempty(result{4})
        options = { options{:} 'absthresh' str2num(result{4}) };
    end;
    if result{5},
         options = { options{:} 'averef', 'on' };
    end;
    if result{6},
         options = { options{:} 'plothist', 'on' };
    end;
    % Begin: Added by Romain on 22 July 2010
    if result{7},
         options = { options{:} 'plotchans', 'on' };
    end;
    % End: Added by Romain on 22 July 2010

else
    options = varargin;
end;

% decode options
% --------------
opt = finputcheck( options, { 'averef'    'string'    { 'on';'off' }       'off';
                              'plothist'  'string'    { 'on';'off' }       'off';
                              'plotchans' 'string'    { 'on';'off' }       'off';
                              'verbose'   'string'    { 'on';'off' }       'off';
                              'elec'      'integer'   []                   [1:EEG.nbchan];
                              'freqlims'  'real'   []                      [35 EEG.srate/2];
                              'specdata'  'real'   []                      [];
                              'specfreqs' 'real'   []                      [];
                              'absthresh' 'real'   []                      [];
                              'stdthresh' 'real'   []                      5 }, 'pop_rejchanspec');
if isstr(opt), error(opt); end;

% compute average referecne if necessary
if strcmpi(opt.averef, 'on')
     NEWEEG = pop_reref(EEG, [], 'exclude', setdiff([1:EEG.nbchan], opt.elec));
else NEWEEG = EEG;
end;
if isempty(opt.specdata)
    [tmpspecdata specfreqs] = pop_spectopo(NEWEEG, 1, [], 'EEG' , 'percent', 100, 'freqrange',[0 EEG.srate/2], 'plot', 'off');
    % add back 0 channels
    devStd = std(EEG.data(:,:), [], 2);
    if any(devStd == 0)
        goodchan  = find(devStd ~= 0);
        specdata  = zeros(length(opt.elec), size(tmpspecdata,2));
        specdata(goodchan,:) = tmpspecdata;
    else
        specdata = tmpspecdata;
    end;
else
    specdata  = opt.specdata;
    specfreqs = opt.specfreqs;
end;

if size(opt.stdthresh,1) == 1 && size(opt.freqlims,1) > 1
    opt.stdthresh = ones(length(opt.stdthresh), size(opt.freqlims,1))*opt.stdthresh;
end;

allrmchan = [];
for index = 1:size(opt.freqlims,1)
    % select frequencies, compute median and std then reject channels
    % ---------------------------------------------------------------
    [tmp fbeg] = min(abs(specfreqs - opt.freqlims(index,1)));
    [tmp fend] = min(abs(specfreqs - opt.freqlims(index,2)));
    selectedspec = mean(specdata(opt.elec, fbeg:fend), 2);
    if ~isempty(opt.absthresh)
        rmchan = find(selectedspec <= opt.absthresh(1) | selectedspec >= opt.absthresh(2));
    else
        m = median(selectedspec);
        s = std( selectedspec);
        nbTresh = size(opt.stdthresh);
        if length(opt.stdthresh) > 1
            rmchan = find(selectedspec <= m+s*opt.stdthresh(index,1) | selectedspec >= m+s*opt.stdthresh(index,2));
        else
            rmchan = find(selectedspec > m+s*opt.stdthresh(index));
        end
    end;

    % print out results
    % -----------------
    if isempty(rmchan)
         textout = sprintf('Range %2.1f-%2.1f Hz: no channel removed\n',  opt.freqlims(index,1), opt.freqlims(index,2));
    else textout = sprintf('Range %2.1f-%2.1f Hz: channels %s removed\n', opt.freqlims(index,1), opt.freqlims(index,2), int2str(opt.elec(rmchan')));
    end;
    fprintf(textout);
    if strcmpi(opt.verbose, 'on')
        for inde = 1:length(opt.elec)
            if ismember(inde, rmchan)
                 fprintf('Elec %s power: %1.2f *\n', EEG.chanlocs(opt.elec(inde)).labels, selectedspec(inde));
            else fprintf('Elec %s power: %1.2f\n', EEG.chanlocs(opt.elec(inde)).labels  , selectedspec(inde));
            end;
        end;
    end;
    allrmchan = [ allrmchan rmchan' ];

    % plot histogram
    % --------------
    if strcmpi(opt.plothist, 'on')
        figure; hist(selectedspec);
        hold on; yl = ylim;
        if ~isempty(opt.absthresh)
            plot([opt.absthresh(1) opt.absthresh(1)], yl, 'r');
            plot([opt.absthresh(2) opt.absthresh(2)], yl, 'r');
        else
            if length(opt.stdthresh) > 1
                threshold1 =  m+s*opt.stdthresh(index,1);
                threshold2 =  m+s*opt.stdthresh(index,2);
                plot([m m], yl, 'g');
                plot([threshold1 threshold1], yl, 'r');
                plot([threshold2 threshold2], yl, 'r');
            else
                threshold =  m+s*opt.stdthresh(index,1);
                plot([threshold threshold], yl, 'r');
            end
        end;
        title(textout);
    end;

end;
allrmchan = unique_bc(allrmchan);

com = sprintf('EEG = pop_rejchan(EEG, %s);', vararg2str(options));
if strcmpi(opt.plotchans, 'on')
    tmpcom = [ 'EEGTMP = pop_select(EEG, ''nochannel'', [' num2str(opt.elec(allrmchan)) ']);' ];
    tmpcom = [ tmpcom ...
            'LASTCOM = ' vararg2str(com) ';' ...
            '[ALLEEG EEG CURRENTSET tmpcom] = pop_newset(ALLEEG, EEGTMP, CURRENTSET);' ...
            '   if ~isempty(tmpcom),' ...
            '     EEG = eegh(LASTCOM, EEG);' ...
            '     eegh(tmpcom);' ...
            '     eeglab(''redraw'');' ...
            '  end; clear EEGTMP tmpcom;' ];

    colors = cell(1,length(opt.elec)); colors(:) = { 'k' };
    colors(allrmchan) = { 'r' }; colors = colors(end:-1:1);
    fprintf('%d electrodes labeled for rejection\n', length(find(allrmchan)));
    tmpchanlocs = EEG.chanlocs;
    if ~isempty(EEG.chanlocs), tmplocs = EEG.chanlocs(opt.elec); tmpelec = { tmpchanlocs(opt.elec).labels }';
    else                       tmplocs = []; tmpelec = mattocell([1:EEG.nbchan]');
    end;
    eegplot(EEG.data(opt.elec,:,:), 'srate', EEG.srate, 'title', 'Scroll component activities -- eegplot()', ...
        'limits', [EEG.xmin EEG.xmax]*1000, 'color', colors, 'eloc_file', tmplocs, 'command', tmpcom);
else
    EEG = pop_select(EEG, 'nochannel', opt.elec(allrmchan));
end;

if nargin < 2
    allrmchan = sprintf('EEG = pop_rejchanspec(EEG, %s);', vararg2str(options));
end;


text:
function [irf_record,D_record,gamma_record,struct_irf_record,irf_estimates,D_estimates,gamma_estimates,strshocks_record,strshocks_estimates]=...
    panel5irf(Y,Xdot,theta_gibbs,sigma_gibbs,Xi,It,Bu,IRFperiods,IRFband,N,n,m,p,k,T,IRFt,favar)

% create first the cells storing the results
irf_record={};
irf_estimates={};

% deal with shocks in turn (there are N*n shocks with this model)
for ii=1:N*n

   % start iterating
   for jj=1:It-Bu

   % draw theta from its posterior distribution
   theta=theta_gibbs(:,jj);

   % create a matrix of zeros of dimension p*(N*n)
   Ysim=zeros(p,N*n);
   %  step 2: set the value of the last row, column i, equal to 1
   Ysim(p,ii)=1;

      % step 4: for each iteration kk, repeat the algorithm for periods T+1 to T+h
      for kk=1:IRFperiods-1

      % use the function lagx to obtain a matrix temp, containing the endogenous regressors
      temp=bear.lagx(Ysim,p-1);

      % define the vector X
      Xsim=[temp(end,:) zeros(1,m)];

      % obtain Xbar et Xtilde
      Xbar=kron(speye(N*n),Xsim);
      Xtilde=Xbar*Xi;

      % obtain the predicted value for T+kk
      yp=Xtilde*theta;

      % concatenate yp at the top of Y
      Ysim=[Ysim;yp'];

      % repeat until values are obtained for T+h
      end

      % record the results from current iteration in cell irf_record
      % loop over variables
      for kk=1:N*n
      % consider column kk of matrix Ysim and trim the (p-1) initial periods: what remains is the series of IRFs for period T to period T+h-1, for variable kk
      temp=Ysim(p:end,kk);
      % record these values in the corresponding matrix of irf_record
      irf_record{kk,ii}(jj,:)=temp';
      end

   % then go for next iteration
   end

% conduct the same process with shocks in other variables
end

% then apply the structural decomposition (if applicable)

% if IRFs have been set to an unrestricted VAR (IRFt=1):
if IRFt==1
% run a pseudo Gibbs sampler to obtain records for D and gamma (for the trivial SVAR)
[D_record,gamma_record]=bear.irfunres(N*n,It,Bu,sigma_gibbs);
struct_irf_record=[];
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(irf_record,N*n,IRFperiods,IRFband,IRFt,[],[],favar);
% if IRFs have been set to an SVAR with Choleski identification (IRFt=2):
elseif IRFt==2
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record,D_record,gamma_record]=bear.irfchol(sigma_gibbs,irf_record,It,Bu,IRFperiods,N*n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,N*n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);
elseif IRFt==3
% run the Gibbs sampler to transform unrestricted draws into orthogonalised draws
[struct_irf_record,D_record,gamma_record]=bear.irftrig(sigma_gibbs,irf_record,It,Bu,IRFperiods,N*n,favar);
% compute posterior estimates
[irf_estimates,D_estimates,gamma_estimates]=bear.irfestimates(struct_irf_record,N*n,IRFperiods,IRFband,IRFt,D_record,gamma_record,favar);
end

% also, if a s structural identification was implemented, compute structural shocks
strshocks_record={};
strshocks_estimates={};
if IRFt~=1
% run the Gibbs sampler
[strshocks_record]=bear.strshockspan5(theta_gibbs,Xi,D_record,Y,Xdot,N*n,k,T,It,Bu);
% obtain point estimates and credibility intervals
[strshocks_estimates]=bear.strsestimates(strshocks_record,N*n,T,IRFband);
% reshape
strshocks_estimates=reshape(strshocks_estimates,[n 1 N]);
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [ a, seed ] = r8sp_random ( m, n, nz_num, row, col, seed )

%*****************************************************************************80
%
%% R8SP_RANDOM sets a random R8SP matrix.
%
%  Discussion:
%
%    The R8SP storage format stores the row, column and value of each nonzero
%    entry of a sparse matrix.
%
%    It is possible that a pair of indices (I,J) may occur more than
%    once.  Presumably, in this case, the intent is that the actual value
%    of A(I,J) is the sum of all such entries.  This is not a good thing
%    to do, but I seem to have come across this in MATLAB.
%
%    The R8SP format is used by CSPARSE ("sparse triplet"), DLAP/SLAP
%    ("nonsymmetric SLAP triad"), by MATLAB, and by SPARSEKIT ("COO" format).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer NZ_NUM, the number of nonzero elements in the matrix.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
%    of the nonzero elements.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(NZ_NUM), the nonzero elements of the matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ a, seed ] = r8vec_uniform_01 ( nz_num, seed );

  return
end


text:
% function x = synsq_cwt_iw(Tx, fs, opt, Cs, freqband)
%
% Inverse Synchrosqueezing transform of Tx with associated
% frequencies in fs and curve bands in time-frequency plane
% specified by Cs and freqband.  This implements Eq. 5 of [1].
%
% 1. G. Thakur, E. Brevdo, N.-S. Fučkar, and H.-T. Wu,
% "The Synchrosqueezing algorithm for time-varying spectral analysis: robustness
%  properties and new paleoclimate applications," Signal Processing, 93:1079-1094, 2013.
%
% Input:
%   Tx, fs: See help synsq_cwt_fw
%   opt: options structure (see help synsq_cwt_fw)
%      opt.type: type of wavelet used in synsq_cwt_fw
%
%      other wavelet options (opt.mu, opt.s) should also match
%      those used in synsq_cwt_fw
%	Cs: (optional) curve centerpoints
%	freqs: (optional) curve bands
%
% Output:
%   x: components of reconstructed signal, and residual error
%
% Example:
%   [Tx,fs] = synsq_cwt_fw(t, x, 32); % Synchrosqueezing
%   Txf = synsq_filter_pass(Tx, fs, -Inf, 1); % Pass band filter
%   xf = synsq_cwt_iw(Txf, fs);  % Filtered signal reconstruction
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo, Gaurav Thakur
%---------------------------------------------------------------------------------
function x = synsq_cwt_iw(Tx, fs, opt, Cs, freqband)
    if nargin<3, opt = struct(); end
    if ~isfield(opt, 'type'), opt.type = 'morlet'; end
	if nargin<4, Cs = ones(size(Tx,2),1); end
	if nargin<5, freqband = size(Tx,1); end

    % Find the admissibility coefficient Cpsi
    Css = synsq_adm(opt.type, opt);

	% Invert Tx around curve masks in the time-frequency plane to recover
	% individual components; last one is the remaining signal
    % Integration over all frequencies recovers original signal
	% factor of 2 is because real parts contain half the energy
	x = zeros(size(Cs,1),size(Cs,2)+1);
	TxMask = zeros(size(Tx));
	TxRemainder = Tx;
	for n=[1:size(Cs,2)]
		TxMask = zeros(size(Tx));
		UpperCs=min(max(Cs(:,n)+freqband(:,n),1),length(fs));
		LowerCs=min(max(Cs(:,n)-freqband(:,n),1),length(fs));
		%Cs==0 corresponds to no curve at that time, so this removes such points from the inversion
		UpperCs(find(Cs(:,n)<1))=1;
		LowerCs(find(Cs(:,n)<1))=2;
		for m=[1:size(Tx,2)]
			TxMask(LowerCs(m):UpperCs(m),m) = Tx(LowerCs(m):UpperCs(m), m);
			TxRemainder(LowerCs(m):UpperCs(m),m) = 0;
		end
		% Due to linear discretization of integral in log(fs), this becomes a simple normalized sum.
		x(:,n) = 1/Css*sum(real(TxMask),1).';
	end
	x(:,n+1) = 1/Css*sum(real(TxRemainder),1).';
	x = x.';
end

author: ebrevdo

repo: synchrosqueezing


text:
function z = isequal( x, y )

% This only works if we treat CVX objects as immutable.
z = x.id_ == y.id_;

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
%ISLIMIT Test if HAL joint angles are within their limits
%
% Returns a matrix the same size as the matrix of joint angles input
% (where each row is a joint angle set), whose elements are false if
% within joint limits or true if not. The formulae by Lenarcic & Umek
% (1994) are used to calculate the shoulder range of motion, with the
% qlim values representing the limits when other joint angles are zero.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% This file modifies file(s) from The Robotics Toolbox for MATLAB (RTB)
% by Peter Corke (www.petercorke.com), see file for statement
%
% Upper arm range of motion:
%  Lenarcic, J. & Umek, A. 1994, 'Simple Model of Human Arm Reachable
%   Workspace', IEEE Transactions on systems, man, and cybernetics,
%   vol. 24, no. 8, pp. 1239-46.
%
% Syntax:
%  (1) I = hal.islimit(q)
%
% Outputs:
%  I : logical matrix, elements which are true are angles outside limit
%
% Inputs:
%  q : Joint row vector, or trajectory matrix of joint row vectors
%
% See also HAL.reachable

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 1993-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license,
% Modified 16/6/2014 (HAL is a subclass of SerialLink)

function I = islimit(hal, q)

I = false(size(q));

for i = [1 4 5 6 7]
    I(:,i) = hal.qlim(i,2) <= q(:,i)|q(:,i) <= hal.qlim(i,1);
end

toplim(:,1) = hal.qlim(2,2) - q(:,1)/6;
toplim(:,2) = hal.qlim(3,2) + ...
    4*q(:,1)/9 - 5*q(:,2)/9 + 10*q(:,1).*q(:,2)/(9*pi);

lowlim(:,1) = hal.qlim(2,1) + q(:,1)/3;
lowlim(:,2) = hal.qlim(3,1) + ...
    7*q(:,1)/9 - q(:,2)/9 + 4*q(:,1).*q(:,2)/(9*pi);

I(:,[2 3]) =  toplim <= q(:,[2 3])|q(:,[2 3]) <= lowlim;

end

author: Allopart

repo: rbpf-gmapping