text:
% tutorial2_spikehistcoupledGLM.m
%
% This is an interactive tutorial designed to walk you through the steps of
% fitting an autoregressive Poisson GLM (i.e., a spiking GLM with
% spike-history) and a multivariate autoregressive Poisson GLM (i.e., a
% GLM with spike-history AND coupling between neurons).
%
% Data: from Uzzell & Chichilnisky 2004; see README file for details.
%
% Last updated: Mar 10, 2020 (JW Pillow)

% Instructions: Execute each section below separately using cmd-enter.
% For detailed suggestions on how to interact with this tutorial, see
% header material in tutorial1_PoissonGLM.m

%% ====  1. Load the raw data ============

% ------------------------------------------------------------------------
% Be sure to unzip the data file data_RGCs.zip
% (http://pillowlab.princeton.edu/data/data_RGCs.zip) and place it in
% this directory before running the tutorial.
% ------------------------------------------------------------------------
% (Data from Uzzell & Chichilnisky 2004):
datdir = 'data_RGCs/';  % directory where stimulus lives
load([datdir, 'Stim']);    % stimulus (temporal binary white noise)
load([datdir,'stimtimes']); % stim frame times in seconds (if desired)
load([datdir, 'SpTimes']); % load spike times (in units of stim frames)
ncells = length(SpTimes);  % number of neurons (4 for this dataset).
% Neurons #1-2 are OFF, #3-4 are ON.
% -------------------------------------------------------------------------

% Compute some basic statistics on the stimulus
dtStim = (stimtimes(2)-stimtimes(1)); % time bin size for stimulus (s)
nT = size(Stim,1); % number of time bins in stimulus

% See tutorial 1 for some code to visualize the raw data!

%% ==== 2. Bin the spike trains =========================
%
% For now we will assume we want to use the same time bin size as the time
% bins used for the stimulus. Later, though, we'll wish to vary this.
tbins = (.5:nT)*dtStim; % time bin centers for spike train binnning
sps = zeros(nT,ncells);
for jj = 1:ncells
    sps(:,jj) = hist(SpTimes{jj},tbins)';  % binned spike train
end

% Let's just visualize the spike-train auto and cross-correlations
% (Comment out this part if desired!)
clf;
nlags = 30; % number of time-lags to use
for ii = 1:ncells
    for jj = ii:ncells
        % Compute cross-correlation of neuron i with neuron j
        xc = xcorr(sps(:,ii),sps(:,jj),nlags,'unbiased');

        % remove center-bin correlation for auto-correlations (for ease of viz)
        if ii==jj, xc(nlags+1) = 0;
        end

        % Make plot
        subplot(ncells,ncells,(ii-1)*ncells+jj);
        plot((-nlags:nlags)*dtStim,xc,'.-','markersize',20);
        axis tight; drawnow;
        title(sprintf('cells (%d,%d)',ii,jj)); axis tight;
    end
end
xlabel('time shift (s)');

%% ==== 3. Build design matrix: single-neuron GLM with spike-history =========

% Pick the cell to focus on (for now).
cellnum = 3;  % 1-2: OFF, 3-4: ON

% Set the number of time bins of stimulus to use for predicting spikes
ntfilt = 25;  % Try varying this, to see how performance changes!
% Set number of time bins of auto-regressive spike-history to use
nthist = 20;

% Build stimulus design matrix (using 'hankel');
paddedStim = [zeros(ntfilt-1,1); Stim]; % pad early bins of stimulus with zero
Xstim = hankel(paddedStim(1:end-ntfilt+1), Stim(end-ntfilt+1:end));

% Build spike-history design matrix
paddedSps = [zeros(nthist,1); sps(1:end-1,cellnum)];
% SUPER important: note that this doesn't include the spike count for the
% bin we're predicting? The spike train is shifted by one bin (back in
% time) relative to the stimulus design matrix
Xsp = hankel(paddedSps(1:end-nthist+1), paddedSps(end-nthist+1:end));

% Combine these into a single design matrix
Xdsgn = [Xstim,Xsp];

% Let's visualize the design matrix just to see what it looks like
subplot(1,10,1:9);
imagesc(1:(ntfilt+nthist), 1:50, Xdsgn(1:50,:));
xlabel('regressor');
ylabel('time bin of response');
title('design matrix (including stim and spike history)');
subplot(1,10,10);
imagesc(sps(1:50,cellnum));
set(gca,'yticklabel', []);
title('spike count');

% The left part of the design matrix has the stimulus values, the right
% part has the spike-history values.  The image on the right is the spike
% count to be predicted.  Note that the spike-history portion of the design
% matrix had better be shifted so that we aren't allowed to use the spike
% count on this time bin to predict itself!

%% === 4. fit single-neuron GLM with spike-history ==================

% First fit GLM with no spike-history
fprintf('Now fitting basic Poisson GLM...\n');
pGLMwts0 = glmfit(Xstim,sps(:,cellnum),'poisson'); % assumes 'log' link and 'constant'='on'.
pGLMconst0 = pGLMwts0(1);
pGLMfilt0 = pGLMwts0(2:end);

% Then fit GLM with spike history (now use Xdsgn design matrix instead of Xstim)
fprintf('Now fitting Poisson GLM with spike-history...\n');
pGLMwts1 = glmfit(Xdsgn,sps(:,cellnum),'poisson');
pGLMconst1 = pGLMwts1(1);
pGLMfilt1 = pGLMwts1(2:1+ntfilt);
pGLMhistfilt1 = pGLMwts1(ntfilt+2:end);

%%  Make plots comparing filters
ttk = (-ntfilt+1:0)*dtStim; % time bins for stim filter
tth = (-nthist:-1)*dtStim; % time bins for spike-history filter

clf; subplot(221); % Plot stim filters
h = plot(ttk,ttk*0,'k--',ttk,pGLMfilt0, 'o-',ttk,pGLMfilt1,'o-','linewidth',2);
legend(h(2:3), 'GLM', 'sphist-GLM','location','northwest');axis tight;
title('stimulus filters'); ylabel('weight');
xlabel('time before spike (s)');

subplot(222); % Plot spike history filter
colr = get(h(3),'color');
h = plot(tth,tth*0,'k--',tth,pGLMhistfilt1, 'o-');
set(h(2), 'color', colr, 'linewidth', 2);
title('spike history filter');
xlabel('time before spike (s)');
ylabel('weight'); axis tight;

%% Plot predicted rate out of the two models

% Compute predicted spike rate on training data
ratepred0 = exp(pGLMconst0 + Xstim*pGLMfilt0);
ratepred1 = exp(pGLMconst1 + Xdsgn*pGLMwts1(2:end));

% Make plot
iiplot = 1:60; ttplot = iiplot*dtStim;
subplot(212);
stem(ttplot,sps(iiplot,cellnum), 'k'); hold on;
plot(ttplot,ratepred0(iiplot),ttplot,ratepred1(iiplot), 'linewidth', 2);
hold off;  axis tight;
legend('spikes', 'GLM', 'hist-GLM');
xlabel('time (s)');
title('spikes and rate predictions');
ylabel('spike count / bin');

%% === 5. fit coupled GLM for multiple-neuron responses ==================

% First step: build design matrix containing spike history for all neurons

Xspall = zeros(nT,nthist,ncells); % allocate space
% Loop over neurons to build design matrix, exactly as above
for jj = 1:ncells
    paddedSps = [zeros(nthist,1); sps(1:end-1,jj)];
    Xspall(:,:,jj) = hankel(paddedSps(1:end-nthist+1),paddedSps(end-nthist+1:end));
end

% Reshape it to be a single matrix
Xspall = reshape(Xspall,nT,[]);
Xdsgn2 = [Xstim, Xspall]; % full design matrix (with all 4 neuron spike hist)

clf; % Let's visualize 50 time bins of full design matrix
imagesc(1:1:(ntfilt+nthist*ncells), 1:50, Xdsgn2(1:50,:));
title('design matrix (stim and 4 neurons spike history)');
xlabel('regressor');
ylabel('time bin of response');

%% Fit the model (stim filter, sphist filter, coupling filters) for one neuron

fprintf('Now fitting Poisson GLM with spike-history and coupling...\n');

pGLMwts2 = glmfit(Xdsgn2,sps(:,cellnum),'poisson');
pGLMconst2 = pGLMwts2(1);
pGLMfilt2 = pGLMwts2(2:1+ntfilt);
pGLMhistfilts2 = pGLMwts2(ntfilt+2:end);
pGLMhistfilts2 = reshape(pGLMhistfilts2,nthist,ncells);

% So far all we've done is fit incoming stimulus and coupling filters for
% one neuron.  To fit a full population model, redo the above for each cell
% (i.e., to get incoming filters for 'cellnum' = 1, 2, 3, and 4 in turn).

%% Plot the fitted filters and rate prediction

clf; subplot(221); % Plot stim filters
h = plot(ttk,ttk*0,'k--',ttk,pGLMfilt0, 'o-',ttk,pGLMfilt1,...
    ttk,pGLMfilt2,'o-','linewidth',2); axis tight;
legend(h(2:4), 'GLM', 'sphist-GLM','coupled-GLM', 'location','northwest');
title(['stimulus filter: cell ' num2str(cellnum)]); ylabel('weight');
xlabel('time before spike (s)');

subplot(222); % Plot spike history filter
colr = get(h(3),'color');
h = plot(tth,tth*0,'k--',tth,pGLMhistfilts2,'linewidth',2);
legend(h(2:end),'from 1', 'from 2', 'from 3', 'from 4', 'location', 'northwest');
title(['coupling filters: into cell ' num2str(cellnum)]); axis tight;
xlabel('time before spike (s)');
ylabel('weight');

% Compute predicted spike rate on training data
ratepred2 = exp(pGLMconst2 + Xdsgn2*pGLMwts2(2:end));

% Make plot
iiplot = 1:60; ttplot = iiplot*dtStim;
subplot(212);
stem(ttplot,sps(iiplot,cellnum), 'k'); hold on;
plot(ttplot,ratepred0(iiplot),ttplot,ratepred1(iiplot),...
    ttplot,ratepred2(iiplot), 'linewidth', 2);
hold off;  axis tight;
legend('spikes', 'GLM', 'sphist-GLM', 'coupled-GLM', 'location', 'northwest');
xlabel('time (s)');
title('spikes and rate predictions');
ylabel('spike count / bin');

%% 6. Model comparison: log-likelihoood and AIC

% Let's compute loglikelihood (single-spike information) and AIC to see how
% much we gain by adding each of these filter types in turn:

LL_stimGLM = sps(:,cellnum)'*log(ratepred0) - sum(ratepred0);
LL_histGLM = sps(:,cellnum)'*log(ratepred1) - sum(ratepred1);
LL_coupledGLM = sps(:,cellnum)'*log(ratepred2) - sum(ratepred2);

% log-likelihood for homogeneous Poisson model
nsp = sum(sps(:,cellnum));
ratepred_const = nsp/nT;  % mean number of spikes / bin
LL0 = nsp*log(ratepred_const) - nT*sum(ratepred_const);

% Report single-spike information (bits / sp)
SSinfo_stimGLM = (LL_stimGLM - LL0)/nsp/log(2);
SSinfo_histGLM = (LL_histGLM - LL0)/nsp/log(2);
SSinfo_coupledGLM = (LL_coupledGLM - LL0)/nsp/log(2);

fprintf('\n empirical single-spike information:\n ---------------------- \n');
fprintf('stim-GLM: %.2f bits/sp\n',SSinfo_stimGLM);
fprintf('hist-GLM: %.2f bits/sp\n',SSinfo_histGLM);
fprintf('coupled-GLM: %.2f bits/sp\n',SSinfo_coupledGLM);

% Compute AIC
AIC0 = -2*LL_stimGLM + 2*(1+ntfilt);
AIC1 = -2*LL_histGLM + 2*(1+ntfilt+nthist);
AIC2 = -2*LL_coupledGLM + 2*(1+ntfilt+ncells*nthist);
AICmin = min([AIC0,AIC1,AIC2]); % the minimum of these

fprintf('\n AIC comparison (smaller is better):\n ---------------------- \n');
fprintf('stim-GLM: %.1f\n',AIC0-AICmin);
fprintf('hist-GLM: %.1f\n',AIC1-AICmin);
fprintf('coupled-GLM: %.1f\n',AIC2-AICmin);

% These are whopping differencess! Clearly coupling has a big impact in
% terms of log-likelihood, though the jump from stimulus-only to
% own-spike-history is greater than the jump from spike-history to
% full coupling.

%% Advanced exercises:
% --------------------
% 1. Write code to simulate spike trains from the fitted spike-history GLM.
% Simulate a raster of repeated responses from the stim-only GLM and
% compare to raster from the spike-history GLM

% 2. Write code to simulate the 4-neuron population-coupled GLM. There are
% now 16 spike-coupling filters (including self-coupling), since each
% neuron has 4 incoming coupling filters (its own spike history coupling
% filter plus coupling from three other neurons.  How does a raster of
% responses from this model compare to the two single-neuron models?

% 3. Compute a non-parametric estimate of the spiking nonlinearity for each
% neuron. How close does it look to exponential now that we have added
% spike history? Rerun your simulations using different non-parametric
% nonlinearity for each neuron. How much improvement do you see in terms of
% log-likelihood, AIC, or PSTH % variance accounted for (R^2) when you
% simulate repeated responses?

author: pillowlab

repo: GLMspiketraintutorial


text:
function jed = ss_to_jed_unix ( s )

%*****************************************************************************80
%
%% SS_TO_JED_UNIX converts a UNIX SS date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S, the UNIX date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%
  jed_epoch = epoch_to_jed_unix ( );

  d = s / ( 24.0 * 60.0 * 60.0 );

  jed = jed_epoch + d;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f=ref_idwilt(c,g,a,M)
%REF_DWILT  Reference Inverse Discrete Wilson Transform
%   Usage:  f=ref_idwilt(c,g,a,M);
%

% Setup transformation matrix.

L=size(g,1);
W=size(c,2);
N=L/a;

F=zeros(L,M*N);

% Weight coefficients.

l=(0:L-1).';

pif=0;

if 1
  % This version uses sines and cosines to express the basis functions.

  for n=0:N/2-1
    % Do the unmodulated coefficient.
    F(:,2*M*n+1)=circshift(g,2*a*n);

    % Setting this to -n*a should produce a time-invariant transform.
    timeinv=0; %-n*a;

    % m odd case
    for m=1:2:M-1
      F(:,m+2*M*n+1)   = sqrt(2)*sin(pi*m/M*(l+timeinv)+pif).*circshift(g,2*n*a);
      F(:,m+2*M*n+M+1) = sqrt(2)*cos(pi*m/M*(l+timeinv)+pif).*circshift(g,(2*n+1)*a);
    end;

    % m even case
    for m=2:2:M-1
      F(:,m+2*M*n+1)     = sqrt(2)*cos(pi*m/M*(l+timeinv)+pif).*circshift(g,2*n*a);
      F(:,m+2*M*n+M+1)   = sqrt(2)*sin(pi*m/M*(l+timeinv)+pif).*circshift(g,(2*n+1)*a);
    end;

    % Most modulated coefficient, Nyquest frequency.
    if mod(M,2)==0
      F(:,M+2*M*n+1)=(-1).^(l+timeinv).*circshift(g,2*n*a);
    else
      F(:,M+2*M*n+1)=(-1).^(l+timeinv).*circshift(g,(2*n+1)*a);
    end;
  end;

else

  % This version uses a cosine,

  for n=0:N/2-1
    % Do the unmodulated coefficient.
    F(:,2*M*n+1)=circshift(g,2*a*n);

    timeinv=-n*a;

    % m odd case
    for m=1:2:M-1
      F(:,m+2*M*n+1)   = sqrt(2)*cos(pi*m/M*(l+timeinv-M/2)+pif).*circshift(g,2*n*a);
      F(:,m+2*M*n+M+1) = sqrt(2)*sin(pi*m/M*(l+timeinv-M/2-a)+pif).*circshift(g,(2*n+1)*a);
    end;

    % m even case
    for m=2:2:M-1
      F(:,m+2*M*n+1)     = sqrt(2)*cos(pi*m/M*(l+timeinv-M/2)+pif).*circshift(g,2*n*a);
      F(:,m+2*M*n+M+1)   = sqrt(2)*sin(pi*m/M*(l+timeinv-M/2-a)+pif).*circshift(g,(2*n+1)*a);
    end;

    % Most modulated coefficient, Nyquest frequency.
    if mod(M,2)==0
      F(:,M+2*M*n+1)=(-1).^(l+timeinv).*circshift(g,2*n*a);
    else
      F(:,M+2*M*n+1)=(-1).^(l+timeinv-a).*circshift(g,(2*n+1)*a);
    end;
  end;

end;

f=F*c;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [res, plotind] = coor2D(this, ind, val, mindist)
% returns x and y coordinates of channels in 2D plane
% FORMAT coor2D(this)
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak, Laurence Hunt
% $Id: coor2D.m 7445 2018-10-12 13:24:48Z vladimir $

megind = indchantype(this, {'MEG', 'MEGPLANAR', 'MEGCOMB'});
eegind = indchantype(this, {'EEG'});
otherind = setdiff(1:nchannels(this), [megind eegind]);

if nargin==1 || isempty(ind)
    if nargin<3 || (size(val, 2)<nchannels(this))
        if ~isempty(megind)
            ind = megind;
        elseif ~isempty(eegind)
            ind = eegind;
        else
            ind = 1:nchannels(this);
        end
    else
        ind = 1:nchannels(this);
    end
elseif ischar(ind)
    switch upper(ind)
        case 'MEG'
            ind = megind;
        case 'EEG'
            ind = eegind;
        case ':'
            ind = 1:nchannels(this);
        otherwise
            ind = otherind;
    end
end

if nargin < 3 || isempty(val)
    if ~isempty(intersect(ind, megind))

        if this.montage.Mind==0
            if ~any(cellfun('isempty', {this.channels(megind).X_plot2D this.channels(megind).Y_plot2D}))
                meg_xy = [this.channels(megind).X_plot2D; this.channels(megind).Y_plot2D];
            elseif all(cellfun('isempty', {this.channels(megind).X_plot2D this.channels(megind).Y_plot2D}))
                meg_xy = grid(length(megind));
            else
                error('Either all or none of MEG channels should have 2D coordinates defined.');
            end
        else
            if ~any(cellfun('isempty', {this.montage.M(this.montage.Mind).channels(megind).X_plot2D ...
                    this.montage.M(this.montage.Mind).channels(megind).Y_plot2D}))
                meg_xy = [this.montage.M(this.montage.Mind).channels(megind).X_plot2D; ...
                    this.montage.M(this.montage.Mind).channels(megind).Y_plot2D];
            elseif all(cellfun('isempty', {this.montage.M(this.montage.Mind).channels(megind).X_plot2D...
                    this.montage.M(this.montage.Mind).channels(megind).Y_plot2D}))
                meg_xy = grid(length(megind));
            else
                error('Either all or none of EEG channels should have 2D coordinates defined.');
            end
        end
    end

    if ~isempty(intersect(ind, eegind))
        if this.montage.Mind==0
            if ~any(cellfun('isempty', {this.channels(eegind).X_plot2D this.channels(eegind).Y_plot2D}))
                eeg_xy = [this.channels(eegind).X_plot2D; this.channels(eegind).Y_plot2D];
            elseif all(cellfun('isempty', {this.channels(eegind).X_plot2D this.channels(eegind).Y_plot2D}))
                eeg_xy = grid(length(eegind));
            else
                error('Either all or none of EEG channels should have 2D coordinates defined.');
            end
        else
            if ~any(cellfun('isempty', {this.montage.M(this.montage.Mind).channels(eegind).X_plot2D ...
                    this.montage.M(this.montage.Mind).channels(eegind).Y_plot2D}))
                eeg_xy = [this.montage.M(this.montage.Mind).channels(eegind).X_plot2D; ...
                    this.montage.M(this.montage.Mind).channels(eegind).Y_plot2D];
            elseif all(cellfun('isempty', {this.montage.M(this.montage.Mind).channels(eegind).X_plot2D...
                    this.montage.M(this.montage.Mind).channels(eegind).Y_plot2D}))
                eeg_xy = grid(length(eegind));
            else
                error('Either all or none of EEG channels should have 2D coordinates defined.');
            end
        end
    end

    if ~isempty(intersect(ind, otherind))
        other_xy = grid(length(otherind));
    end

    xy = zeros(2, length(ind));
    plotind = zeros(1, length(ind));
    for i = 1:length(ind)
        [found, loc] = ismember(ind(i), megind);
        if found
            xy(:, i) = meg_xy(:, loc);
            plotind(i) = 1;
        else
            [found, loc] = ismember(ind(i), eegind);
            if found
                xy(:, i) = eeg_xy(:, loc);
                plotind(i) = 2;
            else
                [found, loc] = ismember(ind(i), otherind);
                if found
                    xy(:, i) = other_xy(:, loc);
                    plotind(i) = 3;
                end
            end
        end
    end

    if nargin > 3 && ~isempty(mindist)
        xy = shiftxy(xy,mindist);
    end

    res = xy;
else
    if this.montage.Mind==0
        this = getset(this, 'channels', 'X_plot2D', ind, val(1, :));
        this = getset(this, 'channels', 'Y_plot2D', ind, val(2, :));
    else
        this.montage.M(this.montage.Mind) = getset(this.montage.M(this.montage.Mind), 'channels', 'X_plot2D', ind, val(1, :));
        this.montage.M(this.montage.Mind) = getset(this.montage.M(this.montage.Mind), 'channels', 'Y_plot2D', ind, val(2, :));
    end
    res = this;
end

function xy = grid(n)

ncol = ceil(sqrt(n));
x = 0:(1/(ncol+1)):1;
x = 0.9*x+0.05;
x = x(2:(end-1));
y = fliplr(x);
[X, Y] = meshgrid(x, y);
xy = [X(1:n); Y(1:n)];

function xy = shiftxy(xy,mindist)

x = xy(1,:);
y = xy(2,:);

l=1;
i=1; %filler
mindist = mindist/0.999; % limits the number of loops
while (~isempty(i) && l<50)
    xdiff = repmat(x,length(x),1) - repmat(x',1,length(x));
    ydiff = repmat(y,length(y),1) - repmat(y',1,length(y));
    xydist= sqrt(xdiff.^2 + ydiff.^2); %euclidean distance between all sensor pairs

    [i,j] = find(xydist<mindist*0.999);
    rm=(i<=j); i(rm)=[]; j(rm)=[]; %only look at i>j

    for m = 1:length(i);
        if (xydist(i(m),j(m)) == 0)
            diffvec = [mindist./sqrt(2) mindist./sqrt(2)];
        else
            xydiff = [xdiff(i(m),j(m)) ydiff(i(m),j(m))];
            diffvec = xydiff.*mindist./xydist(i(m),j(m)) - xydiff;
        end
        x(i(m)) = x(i(m)) - diffvec(1)/2;
        y(i(m)) = y(i(m)) - diffvec(2)/2;
        x(j(m)) = x(j(m)) + diffvec(1)/2;
        y(j(m)) = y(j(m)) + diffvec(2)/2;
    end
    l = l+1;
end

xy = [x; y];

author: spm


text:
function hermite_polynomial_test06 ( )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST06 tests H_QUADRATURE_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST06:\n' );
  fprintf ( 1, '  H_QUADRATURE_RULE computes the quadrature rule\n' );
  fprintf ( 1, '  associated with H(n,x);\n' );

  n = 7;
  [ x, w ] = h_quadrature_rule ( n );

  r8vec2_print ( n, x, w, '      X            W' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use the quadrature rule to estimate:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Q = Integral ( -oo < X < +00 ) X^E exp(-X^2) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   E       Q_Estimate      Q_Exact\n' );
  fprintf ( 1, '\n' );

  for e = 0 : 2 * n - 1
    if ( e == 0 )
      f = ones ( n, 1 );
    else
      f(1:n) = x(1:n).^e;
    end
    q = w' * f;
    q_exact = h_integral ( e );
    fprintf ( 1, '  %2d  %14g  %14g\n', e, q, q_exact );
  end

  return
end


text:
function [lambda] = spm_mci_adjoint_int (U,P,M,V,djdx,tol)
% Integrate adjoint equation
% FORMAT [lambda] = spm_mci_adjoint_int (U,P,M,V,djdx,tol)
%
% U         Inputs
% P         Parameters
% M         Model structure
% V         states
% djdx      derivative of log likelihood wrt states
% tol       tolerances
%
% lambda    adjoint parameters, at times M.t
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: spm_mci_adjoint_int.m 6548 2015-09-11 12:39:47Z will $

init_t=M.t(1);
final_t=M.T;

tol.abs=1e-3;
tol.rel=1e-5;

% parameters for the intergrator
options = odeset('AbsTol',tol.abs,'RelTol',tol.rel);

% integrate adjoint equation backwards
lam_init=zeros(M.n,1);
[TT, ld] = ode15s(@(t,lam) integrate_adjoints(t,lam,djdx,V,U,P,M),[final_t init_t], lam_init, options);

% interpolate to M.t
% (Note, interp1q is not accurate enough here
% as TT typically has less resolution than M.t)
lambda = interp1(TT,ld,M.t,'spline');

end

%--------------------------------------------------------------------------
function DlDt = integrate_adjoints(t,lam,djdx,V,U,P,M)

Nx=M.n;
DlDt = zeros(Nx,1);

% Interpolate state and dj/dx to current time point
for d=1:Nx,
    local_djdx(d)=interp1q(M.t,djdx(:,d),t);
    v(d)=interp1q(M.t,V(:,d),t);
end

% Find nearest time point for which we have pre-computed input
if isempty(U)
    ut=[];
else
    [tmp,ind]=min(abs(t-M.t));
    ut=U(:,ind);
end

% Evaluate state Jacobian
if isfield(M,'dfdx')
    Fx = feval(M.dfdx,v(:),ut,P,M);
else
    Fx = spm_diff(M.f,v(:),ut,P,M,1);
end

DlDt=local_djdx-lam'*Fx;
DlDt=DlDt';

end


text:
function [lc,dflag,dattype]=loadcell(fname,delim,exclusions,options);
%function [lc,dflag,numdata]=loadcell(fname,delim,exclusions);
%
%  loadcell loads a cell array with character delimited
%  data, which can have variable length lines and content.
%  Numeric values are converted from string to double
%  unless options is a string containing 'string'.
%
%  loadcell is for use with small datasets. It is not optimised
%  for large datasets.
%
%  fname is the filename to be loaded
%
%  delim is/are the relevant delimiter(s). If char(10) is included
%  newlines are simply treated as delimiters and a 1-d array is created.
%
%  exclusions  are the set of characters to be treated as paired
%    braces: line ends or delimiters within braces are ignored.
%    braces are single characters and any brace can pair with
%    any other brace: no type pair checking is done.
%
%  options can be omitted or can contain 'string' if no numeric
%    conversion is required, 'single' if multiple adjacent seperators
%    should not be treated as one, 'free' if all linefeeds should be stripped
%    first and 'empty2num' if empty fields are to be treated as numeric
%    zeros rather than an empty character set. Combine options using
%    concatenation.
%
%  lc is a cell array containing the loaded data.
%
%  dflag is a set of flags denoting the (i,j) values where data was entered
%  dflag(i,j)=1 implies lc(i,j) was loaded from the data, and not just set
%  to empty, say, by default.
%
%  numdata is an array  numdata(i,j)=NaN implies
%    lc(i,j) is a string, otherwise it stores the number at i,j.
%    This will occur regardless of whether the 'string' option is set.
%
%  lc will return -1 if the file is not found or could not be
%  opened.
%
%  Hint: numdata+(1/dflag-1) provides a concise descriptor for the numeric data
%  Inf=not loaded
%  NaN=was string or empty set.
%  otherwise numeric
%
%  EXAMPLE
%
%[a,b]=loadcell('resultsfile',[',' char(9)],'"','single-string');
%   will load file 'resultsfile' into variable a, treating any of tab or
%   comma as delimiters. Delimiters or carriage returns lying
%   between two double inverted commas will be ignored. Two adjacent delimiters
%   will count twice, and all data will be kept as a string.
%
%   Note: in space-separated data 'single' would generally be omitted,
%   wheras in comma-separated data it would be included.
%
%   Note the exclusion characters will remain in the final data, and any data
%   contained within or containing exclusion characters will not be
%   converted to numerics.
%
%   (c) Amos Storkey 2002
%    v b160702

% MATLAB is incapable of loading variable length lines or variable type values
% with a whole file command under the standard library sets. This mfile
% fills that gap.
if (nargin<4)
    options=' ';
end;
dflag = [];
%Open file
fid=fopen(fname,'rt');
%Cannot open: return -1
if (fid<0)
  lc=-1;
else
  fullfile=fread(fid,'uchar=>char')';
  %Strip LF if free is set
  if ~isempty(findstr(options,'free'))
      fullfile=strrep(fullfile,char(10),'');
  end;
  %Find all delimiters
  delimpos=[];
  for s=1:length(delim)
    delimpos=[delimpos find(fullfile==delim(s))];
  end
  %Find all eol
  endpos=find(fullfile==char(10));
  endpos=setdiff(endpos,delimpos);
  %find all exclusions
  xclpos=[];
  for s=1:length(exclusions);
    xclpos=[xclpos find(fullfile==exclusions(s))];
  end
  sort(xclpos);
  xclpos=[xclpos(1:2:end-1);xclpos(2:2:end)];
  %Combine eol and delimiters
  jointpos=union(delimpos,endpos);
  t=1;
  %Remove delim/eol within exclusion pairs
  removedelim=[];
  for s=1:length(jointpos)
    if any((jointpos(s)>xclpos(1,:)) & (jointpos(s)<xclpos(2,:)))
      removedelim(t)=jointpos(s);
      t=t+1;
    end;

  end
  %and add start point
  jointpos=[0 setdiff(jointpos,removedelim)];
  i=1;
  j=1;
  posind=1;
  multflag=isempty(findstr(options,'single'));
  stringflag=~isempty(findstr(options,'string'));
  emptnum=~isempty(findstr(options,'empty2num'));
  %Run through
  while (posind<(length(jointpos)))
    %Get current field
    tempstr=fullfile(jointpos(posind)+1:jointpos(posind+1)-1);
    %If empty only continue if adjacent delim count.
    if ~(isempty(tempstr) & multflag);
      %This ij is set
      dflag(i,j)=1;
      %Convert to num
      tempno=str2double([tempstr]);
      %If emptystring convert to zero if emptnum set
      if (isempty(tempstr) & emptnum)
          tempno=0;
      end;
      %Set dattype to no (or NaN if not a num
      dattype(i,j)=tempno;
      %If NaN set lc to string else to num if stringflag not set
      if (isnan(tempno) |  stringflag)
        lc{i,j}=tempstr;
      else
        lc{i,j}=tempno;
      end;
      %Next j
      j=j+1;
    end;
    %If eol inc i and reset j
    if ismember(jointpos(posind+1),endpos)
        i=i+1;
        j=1;
    end;
    %Inc to next delim
    posind=posind+1;
  end;
end;
%Logicalise dflag
dflag=logical(dflag);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [P orderstruct] = mmtimes(varargin)
% P = mmtimes(M1, M2, ... Mn)
%   return a chain matrix product P = M1*M2* ... *Mn
%
% {Mi} are matrices with compatible dimension: size(Mi,2) = size(Mi+1,1)
%
% Because the matrix multiplication is associative; the chain product can
% be carried out with different order, leading to the same result (up to
% round-off error). MMTIMES uses "optimal" order of binary product to
% reduce the computational effort (probably the accuracy is also improved).
%
% The function assumes the cost of the product of (m x n) with (n x p)
% matrices is (m*n*p). This assumption is typically true for full matrix.
%
% Notes:
%   Scalar matrix are groupped together, and the rest will be
%   multiplied with optimal order.
%
%   To get the the structure that stores the best order, call with the
%   second outputs:
%   >> [P orderstruct] = mmtimes(M1, M2, ... Mn);
%   % This structure can be used later if the input matrices have the
%   % same sizes as those in the first call (but with different contents)
%   >> P = mmtimes(M1, M2, ... Mn, orderstruct);
%
% See also: mtimes
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Orginal: 19-Jun-2010
%          20-Jun-2010: quicker top-down algorithm
%          23-Jun-2010: treat the case of scalars
%          16-Aug-2010: passing optimal order as output/input argument

Matrices = varargin;

buildexpr = false;
if ~isempty(Matrices) && isstruct(Matrices{end})
    orderstruct = Matrices{end};
    Matrices(end) = [];
else
    % Detect scalars
    iscst = cellfun('length',Matrices) == 1;
    if any(iscst)
        % scalars are multiplied apart
        cst = prod([Matrices{iscst}]);
        Matrices = Matrices(~iscst);
    else
        cst = 1;
    end
    % Size of matrices
    szmats = [cellfun('size',Matrices,1) size(Matrices{end},2)];
    s = MatrixChainOrder(szmats);

    orderstruct = struct('cst', cst, ...
                         's', s, ...
                         'szmats', szmats);

    if nargout>=2
        % Prepare to build the string expression
        vnames = arrayfun(@inputname, 1:nargin, 'UniformOutput', false);
        % Default names, e.g., M1, M2, ..., for inputs that is not single variable
        noname = cellfun('isempty', vnames);
        vnames(noname) = arrayfun(@(i) sprintf('M%d', i), find(noname), 'UniformOutput', false);
        if any(iscst)
            % String '(M1*M2*...)' for constants
            cstexpr = strcat(vnames(iscst),'*');
            cstexpr = strcat(cstexpr{:});
            cstexpr = ['(' cstexpr(1:end-1) ')'];
        else
            cstexpr = '';
        end
        vnames = vnames(~iscst);
        buildexpr = true;
    end
end

if ~isempty(Matrices)
    P = ProdEngine(1,length(Matrices),orderstruct.s,Matrices);
    if orderstruct.cst~=1
        P = orderstruct.cst*P;
    end
    if buildexpr
        expr = Prodexpr(1,length(Matrices),orderstruct.s,vnames);
        if ~isempty(cstexpr)
            % Concatenate the constant expression in front
            expr = [cstexpr '*' expr];
        end
        orderstruct.expr = expr;
    end
else
    P = orderstruct.cst;
    if nargout>=2
        orderstruct.expr = cstexpr;
    end
end

end % mmtimes

%%
function [s qmin] = MatrixChainOrder(szmats)
% Find the best ordered chain-product, the best splitting index
% of M(i)*...*M(j) is stored in s(j,i) of the array s (only the lower
% part is filled)
% Top-down dynamic programming, complexity O(n^3)

n = length(szmats)-1;
s = zeros(n);

pk = szmats(2:n);
ij = (0:n-1)*(n+1)+1;
left = zeros(1,n-1);
right = zeros(1,n-1);
L = 1;
while true % off-diagonal offset
    q = zeros(size(pk));
    for j=1:n-L % this is faster and BSXFUN or product with DIAGONAL matrix
        q(:,j) = (szmats(j)*szmats(j+L+1))*pk(:,j);
    end
    q = q + left + right;
    [qmin loc] = min(q, [], 1);
    s(ij(1:end-L)+L) = (1:n-L)+loc;

    if L<n-1
        pk = [pk(:,1:end-1);
              pk(end,2:end)];
        left = [left(:,1:end-1);
                qmin(1:end-1)];
        right = [qmin(2:end);
                 right(:,2:end)];
        L = L+1;
    else
        break
    end % if
end % while-loop

end % MatrixChainOrder

%%
function P = ProdEngine(i,j,s,Matrices)
% Perform matrix product from the optimal order, recursive engine
if i==j
    P = Matrices{i};
else
    k = s(j,i);
    P = ProdEngine(i,k-1,s,Matrices)*ProdEngine(k,j,s,Matrices);
end

end

%%
function expr = Prodexpr(i,j,s,vnames)
% Return the string expression of the optimal order
if i==j
    expr = vnames{i};
else
    k = s(j,i);
    expr = ['(' Prodexpr(i,k-1,s,vnames) '*' Prodexpr(k,j,s,vnames) ')'];
end

end % Prodexpr

author: Sable

repo: mcbench-benchmarks


text:
function airway_tree = PTKCloseBranchesInTree(airway_tree, closing_size_mm, image_size, reporting)
    % PTKCloseBranchesInTree. Takes a segmented airway tree and
    %     performs a morphological closing on each segment, and between each segment
    %     and its child segments.
    %
    % This function is used by PTKAirwayRegionGrowingWithExplosionControl.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if ~exist('reporting', 'var')
        reporting = [];
    end

    number_of_segments = airway_tree.CountBranches;

    reporting.ShowProgress('Closing branches in the airway tree');
    segments_to_do = airway_tree;
    segment_index = 0;

    while ~isempty(segments_to_do)
        segment = segments_to_do(end);
        segments_to_do(end) = [];
        segments_to_do = [segments_to_do, segment.Children];

        if ~isempty(reporting)
            if reporting.HasBeenCancelled
                reporting.Error(CoreReporting.CancelErrorId, 'User cancelled');
            end

            progress_value = round(100*(segment_index/number_of_segments));
            reporting.UpdateProgressValue(progress_value);
        end
        segment_index = segment_index + 1;

        CloseSegment(segment, closing_size_mm, image_size);
    end
end

function CloseSegment(segment, closing_size_mm, image_size)
    voxel_indices = segment.GetAcceptedVoxels;

    if ~isempty(voxel_indices)
        all_points = GetClosedIndicesForSegmentAndChildren(segment, closing_size_mm, image_size);
        new_points = setdiff(all_points, voxel_indices);
        segment.AddClosedPoints(unique(new_points));
    end
end

function result_points = GetClosedIndicesForSegmentAndChildren(current_segment, closing_size_mm, image_size)
    segment_indices = current_segment.GetAcceptedVoxels;
    if isempty(current_segment.Children)
        result_points = GetClosedIndices(segment_indices, closing_size_mm, image_size);
    else
        result_points = [];
        for child_segment  = current_segment.Children
            child_indices = child_segment.GetAcceptedVoxels;
            all_indices = [segment_indices; child_indices];
            new_points_all = GetClosedIndices(all_indices, closing_size_mm, image_size);
            new_points_children = GetClosedIndices(child_indices, closing_size_mm, image_size);

            new_points = setdiff(new_points_all, new_points_children);
            result_points = [result_points; new_points];
        end
    end
    result_points = unique(result_points);
end

function new_points = GetClosedIndices(voxel_indices, closing_size_mm, image_size)
    [offset, segment_image, ~] = MimImageCoordinateUtilities.GetMinimalImageForIndices(voxel_indices', image_size);
    border_size = 3;
    bordered_segment_image = PTKImage(segment_image);
    bordered_segment_image.AddBorder(border_size);
    bordered_segment_image.BinaryMorph(@imclose, closing_size_mm);

    border_offset = border_size;
    bordered_image_size = bordered_segment_image.ImageSize;
    all_points = find(bordered_segment_image.RawImage(:));
    new_points = MimImageCoordinateUtilities.OffsetIndices(int32(all_points), -border_offset + int32(offset), int32(bordered_image_size), int32(image_size));
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [R_sq]=R_squared(y,yf)

% function [R_squared]=R_squared(y,yf)
% ------------------------------------------------------------------------
% This function calculates the coefficient of determination (R-Squared) for
% the (e.g. measurement) data |y| and the data |yf| (e.g. a model fit).
% NaN entries in either inputs are ignored.
%
%
% Change log:
% 2008/10/05
% 2019/06/26 Updated description
% 2019/06/26 Updated handling of NaN entries in the data
% ------------------------------------------------------------------------

%%

L=~isnan(y) & ~isnan(yf); %Logic for non-NaN entries

y_bar=mean(y(L)); %Mean of y
SS_err=sum((y(L)-yf(L)).^2); %Sum of differences Of y with yf
SS_tot=sum((y(L)-y_bar).^2); %Sum of differences of y with mean
R_sq=1-(SS_err./SS_tot); %R-squared value

%%
% _*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:
%% cart2im
% Below is a demonstration of the features of the |cart2im| function

%%
clear; close all; clc;

%% Syntax
% |[I,J,K]=cart2im(X,Y,Z,v);|

%% Description
% This function converts the cartesian coordinates X,Y,Z to image
% coordinates I,J,K using the voxel dimension v.
%
% X,Y,Z can be scalars, vectors or matrices.
% v is a vector of length 3 where v(1), v(2) and v(3) correspond to the
% voxel dimensions in the x,y and z direction respectively.

%% Examples

%%
% Plot settings
cMap=gjet(250);
faceAlpha1=0.5;
edgeColor1='none';
edgeColor2='none';
fontSize=15;

%% Example: Using |im2cart| and |cart2im| to convert from image to real world coordinates

% Get a 3D image
load mri;
M=squeeze(D); %example image data set
v=[2 2 5]; %example voxel size, not voxels are ellongated in slice direction

%%
% The voxels to display can be specified as a list (vector) of voxels
% numbers (linear indices) or using a mask (logic array).

%Defining row, column and slice indicices for slice patching
sliceIndexI=round(size(M,1)/2); %(close to) middle row
sliceIndexJ=round(size(M,2)/2); %(close to) middle column
sliceIndexK=round(size(M,3)/2); %(close to) middle slice

%Defining "masks" i.e. logic arrays with ones for voxels of interest
logicSliceI=false(size(M));
logicSliceI(sliceIndexI,:,:)=1;
logicSliceI=logicSliceI & M>0;

logicSliceJ=false(size(M));
logicSliceJ(:,sliceIndexJ,:)=1;
logicSliceJ=logicSliceJ & M>0;

logicSliceK=false(size(M));
logicSliceK(:,:,sliceIndexK)=1;
logicSliceK=logicSliceK & M>0;

% Creating patch data
[Fx,Vx,Cx]=ind2patch(logicSliceJ,M,'sj');
[Fy,Vy,Cy]=ind2patch(logicSliceI,M,'si');
[Fz,Vz,Cz]=ind2patch(logicSliceK,M,'sk');

% Convert image coordinates to cartesian coordinates
[Vx(:,1),Vx(:,2),Vx(:,3)]=im2cart(Vx(:,2),Vx(:,1),Vx(:,3),v);
[Vy(:,1),Vy(:,2),Vy(:,3)]=im2cart(Vy(:,2),Vy(:,1),Vy(:,3),v);
[Vz(:,1),Vz(:,2),Vz(:,3)]=im2cart(Vz(:,2),Vz(:,1),Vz(:,3),v);

h8=cFigure;
title('MRI visualisation, slices and voxels in cartesian coordinates with aid of voxel size');
xlabel('X (mm)');ylabel('Y (mm)'); zlabel('Z (mm)'); hold on;
hp2= patch('Faces',Fx,'Vertices',Vx,'FaceColor','flat','CData',Cx,'EdgeColor',edgeColor2,'FaceAlpha',faceAlpha1);
hp3= patch('Faces',Fy,'Vertices',Vy,'FaceColor','flat','CData',Cy,'EdgeColor',edgeColor2,'FaceAlpha',faceAlpha1);
hp4= patch('Faces',Fz,'Vertices',Vz,'FaceColor','flat','CData',Cz,'EdgeColor',edgeColor2,'FaceAlpha',faceAlpha1);
axis equal; view(3); axis tight; axis vis3d; grid on;
colormap(gray(250)); colorbar;
camlight headlight;
set(gca,'fontSize',fontSize);
drawnow;

%%
% Get example cartesian coordinates to map to image coordinates

vMid=mean(Vz,1);
plotV(vMid,'b.','MarkerSize',50);

%%
% Map to image coordinates using |cart2im|
[i,j,k]=im2cart(vMid(:,2),vMid(:,1),vMid(:,3),v);

%Image coordinates for the point are
disp(num2str([i j k]));

%Which means the indices for the voxel containing the point are
disp(num2str(round([i j k])));
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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/>.


text:
function dif_print ( ntab, xtab, diftab, title )

%*****************************************************************************80
%
%% DIF_PRINT prints the polynomial represented by a divided difference table.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NTAB, the dimension of the arrays DIFTAB and XTAB.
%
%    Input, real XTAB(NTAB), the X values for the polynomial.
%
%    Input, real DIFTAB(NTAB), the divided difference table
%    for the polynomial.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  p(x) =                           %14f\n', diftab(1) );

  for i = 2 : ntab
    fprintf ( 1, '       + ( x - %14f) * ( %14f\n', xtab(i-1), diftab(i) );
  end

  fprintf ( 1, '  ' );
  for i = 1 : ntab-1
    fprintf ( 1, ')' );
  end
  fprintf ( 1, '\n' );

  return
end


text:
function r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_PRINT_SOME prints out a portion of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%7d  ', i );

      for j = j2lo : j2hi
        fprintf ( 1, '%12g  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
%LAPLACIAN_PERRINX   Compute surface Laplacian of EEG data.
% [surf_lap,G,H] = laplacian_perrinX(data,x,y,z[,leg_order,smoothing]);
%
% INPUTS    :
%      data : EEG data (can be N-D, but first dimension must be electrodes)
%     x,y,z : x,y,z coordinates of electrode positions (e.g., [EEG.chanlocs.X])
%
% (optional inputs)
% leg_order : order of Legendre polynomial (default is 20 [40 for >100 electrodes])
% smoothing : G smoothing parameter (lambda), set to 1e-5 by default
%
%
% OUTPUTS   :
%  surf_lap : the surface Laplacian (second spatial derivative)
% (optional outputs)
%       G,H : G and H matrices
%
%  This is an implementation of algorithms described by
%    Perrin, Pernier, Bertrand, and Echallier (1989). PubMed #2464490

% mikexcohen@gmail.com

function [surf_lap,G,H] = laplacian_perrinX(data,x,y,z,varargin) % vararg order: leg_order,smoothing

numelectrodes = numel(x);

if nargin<4
    help laplacian_perrinX
    error('Read help file!')
end

%% compute G and H matrices

% initialize
G=zeros(numelectrodes);
H=zeros(numelectrodes);
cosdist=zeros(numelectrodes);

% default parameters for +/- 100 electrodes
if numelectrodes>100
    m=3; leg_order=40;
else
    m=4; leg_order=20;
end

if numel(varargin)>0 && ~isempty(varargin{1})
    leg_order=varargin{1};
end

% scale XYZ coordinates to unit sphere
[junk,junk,spherical_radii] = cart2sph(x,y,z);
maxrad = max(spherical_radii);
x = x./maxrad;
y = y./maxrad;
z = z./maxrad;

for i=1:numelectrodes
    for j=i+1:numelectrodes
        cosdist(i,j) = 1 - (( (x(i)-x(j))^2 + (y(i)-y(j))^2 + (z(i)-z(j))^2 ) / 2 );
    end
end
cosdist = cosdist+cosdist' + eye(numelectrodes);

% compute Legendre polynomial
legpoly = zeros(leg_order,numelectrodes,numelectrodes);
for ni=1:leg_order
    temp = legendre(ni,cosdist);
    legpoly(ni,:,:) = temp(1,:,:);
end

% precompute electrode-independent variables
twoN1  = 2*(1:leg_order)+1;
gdenom = ((1:leg_order).*((1:leg_order)+1)).^m;
hdenom = ((1:leg_order).*((1:leg_order)+1)).^(m-1);

for i=1:numelectrodes
    for j=i:numelectrodes

        g=0; h=0;

        for ni=1:leg_order
            % compute G and H terms
            g = g + (twoN1(ni)*legpoly(ni,i,j)) / gdenom(ni);
            h = h - (twoN1(ni)*legpoly(ni,i,j)) / hdenom(ni);
        end
        G(i,j) =  g/(4*pi);
        H(i,j) = -h/(4*pi);
    end
end

% mirror matrix
G=G+G'; H=H+H';

% correct for diagonal-double
G = G-eye(numelectrodes)*G(1)/2;
H = H-eye(numelectrodes)*H(1)/2;

%% compute laplacian

% reshape data to electrodes X time/trials
orig_data_size = squeeze(size(data));
if any(orig_data_size==1)
    data=data(:);
else
    data = reshape(data,orig_data_size(1),prod(orig_data_size(2:end)));
end

% smoothing constant
if numel(varargin)==2
    smoothing=varargin{2};
else
    smoothing=1e-5;
end

% add smoothing constant to diagonal
% (change G so output is unadulterated)
Gs = G + eye(numelectrodes)*smoothing;

% compute C matrix
GsinvS = sum(inv(Gs));
dataGs = data'/Gs;
C      = dataGs - (sum(dataGs,2)/sum(GsinvS))*GsinvS;

% compute surface Laplacian (and reshape to original data size)
surf_lap = reshape((C*H')',orig_data_size);

%% end

author: mikexcohen

repo: AnalyzingNeuralTimeSeries


text:
function MountainCarPlot( x,a,steps )
subplot(2,1,2);
set(gco,'BackingStore','off')  % for realtime inverse kinematics
set(gco,'Units','data')
xplot =-1.6:0.05:0.6;
yplot =sin(3*xplot);
%Mountain
h = area(xplot,yplot,-1.1);
set(h,'FaceColor',[.1 .7 .1])
hold on
% Car  [1 .7 .1]
plot([x(1)-0.075 x(1)+0.075] ,[sin(3*(x(1)-0.075))+0.2  sin(3*(x(1)+0.075))+0.2 ],'-','LineWidth',10,'Color',[1 .7 .1]);
% wheels
plot(x(1)-0.05,sin(3*(x(1)-0.05))+0.06,'ok','markersize',12,'MarkerFaceColor',[.5 .5 .5]);
plot(x(1)+0.05,sin(3*(x(1)+0.05))+0.06,'ok','markersize',12,'MarkerFaceColor',[.5 .5 .5]);

%Goal
plot(0.45,sin(3*0.5)+0.1,'-pk','markersize',15,'MarkerFaceColor',[1 .7 .1]);
% direction of the force
if (a<0)
      plot(x(1)-0.08-0.05,sin(3*(x(1)-0.05))+0.2,'<k','MarkerFaceColor','g','markersize',10);
elseif (a>0)
      plot(x(1)+0.08+0.05,sin(3*(x(1)+0.05))+0.2,'>k','MarkerFaceColor','g','markersize',10);
end

%ctitle(strcat ('Step: ',int2str(steps)));
%-----------------------
axis([-1.6 0.6 -1.1 1.5]);
drawnow
hold off

author: epfl-lasa


text:
function ctranspose(f) %#ok<*INUSD>
%CTRANSPOSE  CHEBTECH objects are not transposable, so this method will throw an
%   error.

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

error('CHEBFUN:CHEBTECH:ctranspose:notPossible', ...
    'CHEBTECH objects are not transposable.')

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [StopFlag, Status] = StopCritDuplo(nfxk,Niter, Nmap, T, MaxNumIter, MaxNumMapEval, TimeLimit, epsilon, Stopping_Crit)
% Function checking that one of the stopping criteria
% holds to terminate LLM and GLM. It perepares the status determining why
% the algorithm is stopped.
%
% USAGE:
%
%    [StopFlag, Status] = StopCritDuplo(nfxk,Niter, Nmap, T, MaxNumIter, MaxNumMapEval, TimeLimit, epsilon, Stopping_Crit)
%
% INPUTS:
%    nhxk:             the norm 2 of `h(xk)`
%    Niter:            the number of iterations
%    Nmap:             the number of mapping calls
%    T:                the running time
%    MaxNumIter:       maximum number of iterations
%    MaxNumMapEval:    maximum number of function evaluations
%    TimeLimit:        maximum running time
%    epsilon:          accuracy parameter
%    Stopping_Crit:    stopping criterion:
%
%                        1. stop if :math:`||nfxk|| \leq \epsilon`
%                        2. stop if `MaxNumIter` is reached
%                        3. stop if `MaxNumMapEval` is reached
%                        4. stop if `TimeLimit` is reached
%                        5. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
% OUTPUTS:
%    StopFlag:         1: if one of the stopping criteria holds, 0: if none of the stopping criteria holds
%    Status:           the reason of the scheme termination

switch Stopping_Crit

  case 1
    if nhxk <= epsilon
      StopFlag = 1;
      Status   = 'A solution of nonlinear system is found.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 2
    if Niter >= MaxNumIter
      StopFlag = 1;
      Status   = 'Maximum number of iterations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 3
    if Nmap >= MaxNumMapEval
      StopFlag = 1;
      Status   = 'Maximum number of mapping evaluations is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 4
    if T >= TimeLimit
      StopFlag = 1;
      Status   = 'Time limit is reached.';
    else
      StopFlag = 0;
      Status   = [];
    end

  case 5
    if (nfxk <= epsilon || Niter >= MaxNumIter)
      StopFlag = 1;
      if Niter < MaxNumIter
          Status = 'a solution is found.';
      else
          Status = 'Maximum number of iterations is reached.';
      end
    else
      StopFlag = 0;
      Status   = [];
    end

end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% End of StopCritDuplo.m %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
%EEGLAB graphic interface functions (popfunc folder):
%  eeg_addnewevents     - EEG_ADDNEWEVENTS Add new events to EEG structure. Both EEG.event and...
%  eeg_amplitudearea    - Resamples an ERP average using spline interpolation...
%  eeg_chaninds         - Look up channel indices in a EEG structure...
%  eeg_context          - Returns (in output 'delays') a matrix giving, for each event of specified...
%  eeg_countepochs      - EEG_COUNTEPOCHS Count how many epochs there are of each type...
%  eeg_decodechan       - Given an input EEG dataset structure, output a new EEG data structure...
%  eeg_dipselect        - Select component dipoles from an EEG dataset with...
%  eeg_eegrej           - Reject porition of continuous data in an EEGLAB...
%  eeg_emptyset         - Initialize an EEG dataset structure with default values.
%  eeg_epoch2continuous - Convert epoched dataset to continuous dataset...
%  eeg_epochformat      - Convert the epoch information of a dataset from struct...
%  eeg_eventformat      - Convert the event information of a dataset from struct...
%  eeg_eventhist        - Return or plot histogram of event or urevent field values.
%  eeg_eventtable       - Returns all events contained in the EEG structure (and...
%  eeg_eventtypes       - Return a list of event or urevent types in a dataset and...
%  eeg_getepochevent    - Return dataset event field values for all events...
%  eeg_getica           - Get ICA component activation. Recompute if necessary.
%  eeg_insertbound      - Insert boundary event in an EEG event structure.
%  eeg_insertboundold   - Insert boundary event in an EEG event structure.
%  eeg_interp           - Interpolate data channels...
%  eeg_laplac           - Gives the laplacian for the data contained in EEG...
%  eeg_lat2point        - Convert latencies in time units relative to the...
%  eeg_latencyur        - Transform latency of sample point in the continuous...
%  eeg_matchchans       - Find closest channels in a larger EEGLAB chanlocs structure...
%  eeg_mergechan        - Merge channel structure while preserving channel...
%  eeg_mergelocs        - Merge channel structure while preserving channel...
%  eeg_mergelocs_diffstruct - Merge channel structure while preserving channel...
%  eeg_multieegplot     - Produce an EEGPLOT of a the average of an epoched dataset...
%  eeg_oldica           - Report, return or add to oldicaweights and oldicasphere...
%  eeg_point2lat        - Convert latency in data points to latency in ms relative...
%  eeg_pv               - Compute EEG.data 'percent variance ' (pv) of the whole EEG data versus the projections...
%  eeg_pvaf             - Compute EEG.data 'percent variance accounted for' (pvaf) by specified components.
%  eeg_rejmacro         - Internal EEGLAB macro for all pop_ functions that...
%  eeg_rejsuperpose     - Superpose rejections of a EEG dataset.
%  eeg_timeinterp       - Perform spline interpolation of a portion...
%  eeg_topoplot         - Plot scalp map...
%  eeg_urlatency        - Find the original (ur) latency of a time point in...
%  getchanlist          - Obtain indices of specified channel types.
%  importevent          - Import experimental events from data file or Matlab...
%  pop_adjustevents     - Adjust event offset of all (or specified) events...
%  pop_autorej          - Perform automatic artifact epoch detection and rejection...
%  pop_averef           - Convert an EEG dataset to average reference.
%  pop_biosig           - Import data files into EEGLAB using BIOSIG toolbox...
%  pop_biosig16         - Import data files into EEGLAB using BIOSIG toolbox...
%  pop_biosig16ying     - Import data files into EEGLAB using BIOSIG toolbox...
%  pop_chancenter       - Recenter cartesian X,Y,Z channel coordinates...
%  pop_chancoresp       - Define correspondences between two channel locations structures...
%  pop_chanedit         - Edit the channel locations structure of an EEGLAB dataset,...
%  pop_chanevent        - Import event latencies from the rising and/or falling 'edge'...
%  pop_chansel          - Pop up a graphic interface to select channels...
%  pop_comments         - Edit comments...
%  pop_compareerps      - Compare the (ERP) averages of two datasets.
%  pop_comperp          - Compute the grand average ERP waveforms of multiple datasets...
%  pop_copyset          - Copy the current EEG dataset into another dataset.
%  pop_crossf           - Return estimates and plots of event-related spectral coherence...
%  pop_editeventfield   - Add/remove/rename/modify a field in the event structure...
%  pop_editeventvals    - Edit events contained in an EEG dataset structure.
%  pop_editset          - Edit EEG dataset structure fields.
%  pop_eegfilt          - Interactively filter EEG dataset data using EEGFILT...
%  pop_eegplot          - Visually inspect EEG data using a scrolling display.
%  pop_eegthresh        - Reject artifacts by detecting outlier values.  This has...
%  pop_envtopo          - Plot envelope of an averaged EEG epoch, plus scalp maps...
%  pop_epoch            - Convert a continuous EEG dataset to epoched data by extracting...
%  pop_erpimage         - Draw an ERP-image plot of a given EEG channel or independent...
%  pop_eventstat        - Computes and plots statistical characteristics of an EEG event,...
%  pop_expevents        - Export events to CSV file...
%  pop_expica           - Export ICA weights or inverse matrix...
%  pop_export           - Export EEG dataset...
%  pop_fileio           - Import data files into EEGLAB using FileIO...
%  pop_fileiodir        - Import directory into EEGLAB using FileIO...
%  pop_findmatchingcomps - Find components with the highest...
%  pop_headplot         - Plot one or more spherically-splined EEG field maps...
%  pop_icathresh        - Main menu for choosing threshold for component...
%  pop_importdata       - Import data from a Matlab variable or disk file by calling...
%  pop_importegimat     - Import EGI Matlab segmented file...
%  pop_importepoch      - Export epoch and/or epoch event information to the event...
%  pop_importerplab     - Import ERPLAB event list file and bin file into...
%  pop_importevent      - Import events into an EEG dataset. If the EEG dataset...
%  pop_importpres       - Append Presentation event file information into an EEGLAB dataset...
%  pop_interp           - Interpolate data channels...
%  pop_jointprob        - Reject artifacts in an EEG dataset using joint...
%  pop_loadbci          - Import a BCI2000 ascii file into EEGLAB...
%  pop_loadset          - Load an EEG dataset. If no arguments, pop up an input window.
%  pop_mergeset         - Merge two or more datasets. If only one argument is given,...
%  pop_newcrossf        - Return estimates and plots of event-related spectral coherence...
%  pop_newset           - Edit/save EEG dataset structure information.
%  pop_newtimef         - Returns estimates and plots of event-related (log) spectral...
%  pop_plotdata         - Plot average of EEG channels or independent components in...
%  pop_plottopo         - Plot one or more concatenated multichannel data epochs...
%  pop_prop             - Plot the properties of a channel or of an independent...
%  pop_readegi          - Load a EGI EEG file (pop out window if no arguments).
%  pop_readlocs         - Load a EGI-format EEG file (pop up an interactive window if no arguments).
%  pop_readsegegi       - Load a segmented EGI EEG file. Pop up query...
%  pop_rejchan          - Reject artifacts channels in an EEG dataset using joint...
%  pop_rejchanspec      - Reject artifacts channels in an EEG dataset using...
%  pop_rejcont          - Reject continuous portions of data based on spectrum...
%  pop_rejepoch         - Reject pre-labeled trials in a EEG dataset.
%  pop_rejkurt          - Rejection of artifact in a dataset using kurtosis...
%  pop_rejspec          - Rejection of artifact in a dataset using...
%  pop_rejtrend         - Measure linear trends in EEG data; reject data epochs...
%  pop_reref            - Convert an EEG dataset to average reference or to a...
%  pop_resample         - Resample dataset (pop up window).
%  pop_rmbase           - Remove channel baseline means from an epoched or...
%  pop_rmdat            - Remove continuous data around specific events...
%  pop_runica           - Run an ICA decomposition of an EEG dataset using RUNICA,...
%  pop_runscript        - Run Matlab script...
%  pop_saveh            - Save the EEGLAB session command history stored in ALLCOM...
%  pop_saveset          - Save one or more EEG dataset structures...
%  pop_select           - Given an input EEG dataset structure, output a new EEG data structure...
%  pop_selectcomps      - Display components with button to visualize their...
%  pop_selectevent      - Find events in an EEG dataset. If the dataset...
%  pop_signalstat       - Computes and plots statistical characteristics of a signal,...
%  pop_snapread         - Load an EEG SnapMaster file (pop out window if no arguments).
%  pop_spectopo         - Plot spectra of specified data channels or components.
%  pop_subcomp          - Remove specified components from an EEG dataset.
%  pop_timef            - Returns estimates and plots of event-related (log) spectral...
%  pop_timtopo          - Call the TIMTOPO function for epoched EEG datasets.
%  pop_topochansel      - Pop up a topographic interface to select channels...
%  pop_topoplot         - Plot scalp map(s) in a figure window. If number of input...
%  pop_writeeeg         - Write EEGLAB dataset to disk in EDF/GDF or BDF format...
%  pop_writelocs        - Load a EGI EEG file (pop out window if no arguments).

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function par_accept_results(tag, payload)
% Accept a result from the scheduler and store it in a result table.
% par_accept_result(Tag, Payload)
%
% This function stores the results produced by the scheduler in a global result table for later
% pickup by par_endschedule. The reason is that the scheduler, being a Java class, cannot hold very
% large (or very many) results since the MATLAB Java memory is limited to a maximum of 64GB.
% Therefore we need to ferry the results out of Java as they are being produced.
%
% The function also attempts to check if the results contain errors, and prints them as they are
% being received.
%
% In:
%   Tag: tag string that uniquely identifies the result; should begin with 'tag__'.
%
%   Payload: payload uint8 array of the same format as returned by the Scheduler

global tracking;
% atomically store payload in global results table
tracking.parallel.results.(tag) = payload;

if length(payload) < 16384
    try
        % for short results we check if we got an error, and if so, we print it out
        decoded = hlp_deserialize(typecast(payload,'uint8'));
        if all(isfield(decoded{2},{'message','identifier','stack'}))
            fprintf('Got an error during parallel execution: %s\n',hlp_handleerror(decoded{2})); end
    catch e %#ok<NASGU>
    end
end

repo: REST


text:
function example_sel2html(fn)
%Convert Dick Benson's example_sel text files to HTML documents

if nargin < 1
    fn = 'example_sel.txt';
end
f = textread(fn,'%s','delimiter','\n','whitespace','');

ln = [];
for n = 1:length(f)
    [st1,fn1,tk1] = regexp(f{n},'^([^\$]*)');
    [st2,fn2,tk2] = regexp(f{n},'\$(.*)$');
    % Blank line.  Also, ignore initialization lines (#)
    if isempty(tk1) | (strncmp(f{n},'#',1) & length(f{n})>1)
        ln(n).text = '';
        ln(n).code = '';
        ln(n).style = '';
    else
        ln(n).text = deblank(f{n}(tk1{1}(1):tk1{1}(2)));
        ln(n).text = strrep(ln(n).text,'>','>');
        ln(n).text = strrep(ln(n).text,'<','<');

        % Determine the indent depending on how many spaces there are
        [st,fn] = regexp(ln(n).text,'^([\s]*)');
        if isempty(fn)
            ln(n).style = 'level1';
        elseif fn<=4
            ln(n).style = 'level2"';
        else
            ln(n).style = 'level3"';
        end

        if ~isempty(tk2)
            ln(n).code = f{n}(tk2{1}(1):tk2{1}(2));
        else
            ln(n).code = '';
        end

    end

    %Look for initialization commands:  Line starts with #
    if strncmp(f{n},'#',1) & length(f{n})>1
        evalin('base',f{n}(2:end),'disp([''fail to execute string: '',f{n}(2:end)])');
        disp('Evaluating Initialization commands');
    end;

end

s{1} = '';
s{end+1} = '<title>Model Based Design and the Path to FPGAs</title>';
s{end+1} = '';
s{end+1} = '';
s{end+1} = '';

for n = 1:length(ln)
    if isempty(ln(n).code)
        if isempty(ln(n).text)
            % No text, no code
            s{end+1} = '
';
        else
            % Text, no code
            s{end+1} = sprintf('
%s
', ...
                ln(n).style,ln(n).text);
        end
    else
        % Text and code
        s{end+1} = sprintf('
%s
', ...
            ln(n).style,ln(n).code,ln(n).text);
    end
end
s{end+1} = '';

web(['text://' s{:}]);

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/9565-model-based-design-and-fpga-implementation-with-simulink/Model Based Design and the FPGA Implementation in Simulink/edge_detection_menu.m


text:
function [n, nm, nl, ts, names, m] = nex_marker(filename, varname)
% nex_marker(filename, varname): Read a marker variable from a .nex file
%
% [n, nm, nl, ts, names, m] = nex_marker(filename, varname)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   varname - variable name
%
%           continuous (a/d) data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector d.
% OUTPUT:
%   n - number of markers
%   nm - number of fields in each marker
%   nl - number of characters in each marker field
%   ts - array of marker timestamps (in seconds)
%   names - names of marker fields ([nm 64] character array)
%   m - character array of marker values [n nl nm]

n = 0;
nm = 0;
nl = 0;
ts = 0;
m = 0;
names = 0;

if(nargin ~= 2)
   disp('2 input arguments are required')
   return
end

if(ischar(filename) == 0)
   disp('input arguments should be character arrays')
   return
end

if(ischar(varname) == 0)
   disp('input arguments should be character arrays')
   return
end

if(isempty(filename))
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
	filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r');
if(fid == -1)
	disp('cannot open file');
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
name = zeros(1, 64);
found = 0;
for i=1:nvar
	type = fread(fid, 1, 'int32');
	var_version = fread(fid, 1, 'int32');
	name = fread(fid, [1 64], 'char');
	offset = fread(fid, 1, 'int32');
	n = fread(fid, 1, 'int32');
	dummy = fread(fid, 32, 'char');
	adfreq = fread(fid, 1, 'double');
	adtomv = fread(fid, 1, 'double');
	npw = fread(fid, 1, 'int32');
	nm = fread(fid, 1, 'int32');
	nl = fread(fid, 1, 'int32');
	dummy = fread(fid, 68, 'char');
	name = char(name);
	name = deblank(name);
	k = strcmp(name, deblank(varname));
	if(k == 1)
		if type ~= 6
			disp(sprintf('%s is not a marker variable', deblank(varname)));
			return;
		end
		found = 1;
		fseek(fid, offset, 'bof');
		ts = fread(fid, [1 n], 'int32');
		names = zeros(1,64);
		m = zeros(n, nl, nm);
		for j=1:nm
			names(j, :) = fread(fid, [1 64], 'char');
			for p = 1:n
				m(p, :, j) = fread(fid, [1 nl], 'char');
			end
		end
		break
	end
end

fclose(fid);

if found == 0
	disp('did not find variable in the file');
else
	names = char(names);
	m = char(m);
	ts = ts/freq;
	disp(strcat('number of markers = ', num2str(n)));
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function z = ge(x,y)
%GE Greater than or equal for sptensors.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Observations for sparse matrix case.
% The result of a >= 5 is sparse.
% The result of a >= 0 is sparse.
% The result of a >= full(a) is sparse.

%% Case 1: One argument is a scalar
if isscalar(y)
    subs1 = x.subs((x.vals >= y),:);
    if y <= 0
        subs2 = setdiff(allsubs(x),x.subs,'rows');
    else
        subs2 = [];
    end
    z = sptensor([subs1;subs2],true,size(x));
    return;
end

% Call back with the arguments reversed.
if isscalar(x)
    z = le(y,x);
    return;
end

%% Case 2: Both x and y are tensors of some sort
% Check that the sizes match
if ~isequal(x.size,y.size)
    error('Size mismatch');
end

% Case 2a: Two sparse tensors
if isa(x,'sptensor') && isa(y,'sptensor')
    z = le(y,x);
    return;
end

% Case 2b: One dense tensor
if isa(y,'tensor')

    % x zero
    subs1 = find(y <= 0);
    subs1 = setdiff(subs1,x.subs,'rows');

    % x nonzero
    subs2 = x.subs(x.vals >= y(x.subs,'extract'),:);

    % assemble
    z = sptensor([subs1;subs2],true,size(x));

    return;

end

%% Otherwise
error('The arguments must be two sptensors or an sptensor and a scalar.');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ac = compressLogicalArray(a)
ac = [];
if isempty(a)
    return;
end
if all(a==1)
    ac = [1,length(a), 1];
    return
end
if all(a==0)
    ac = [1,length(a), 0];
    return
end
k = find(a==0);
kk = 1;

% Generate all zero intervals
for ii = 1:length(k)
    if ii == 1
        iS = k(ii);
        iE = k(ii);
    end
    if ii>1 && k(ii-1)+1 == k(ii)
        iE = iE+1;
    elseif ii>1 && k(ii-1)+1 < k(ii)
        ac(kk,:) = [iS, iE, 0]; %#ok<*AGROW>
        kk = kk+1;
        iS = k(ii);
        iE = k(ii);
    end
    if ii==length(k)
        ac(kk,:) = [iS, iE, 0];
    end
end

% Generate all one intervals
kk = 1;
for ii = 1:size(ac,1)
    if ii == 1
        if ac(ii,1) == 1
            ac2(kk,:) = ac(ii,:);
        else
            ac2(kk,:) = [1, ac(ii,1)-1, 1];
            kk = kk+1;
            ac2(kk,:) = ac(ii,:);
        end
    else
        ac2(kk,:) = [ac(ii-1,2)+1, ac(ii,1)-1, 1];
        kk = kk+1;
        ac2(kk,:) = ac(ii,:);
    end
    kk = kk+1;
    if ii==size(ac,1) && ac(ii,2)<length(a)
        ac2(kk,:) = [ac(ii,2)+1, length(a), 1];
        kk = kk+1;
    end
end

ac = ac2;

author: BUNPC


text:
function dy = lorenz(t,y)

% Lorenz's ODE

dy = [ -10*(y(1)-y(2)); ...
       28*y(1)-y(2)-y(1)*y(3); ...
       y(1)*y(2)-2.6666666666*y(3)];

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function point_num = xyz_example_size ( )

%*****************************************************************************80
%
%% XYZ_EXAMPLE_SIZE sizes an example XYZ dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer POINT_NUM, the number of points.
%
  point_num = 101;

  return
end


text:
Nneurons = [10,10]; % the number of neurons in each sequence
Dt = [10,10]; % the number of time steps between each neuron in the sequence
Pseq = [1,1]; % the probability of the sequence occuring
NeuronNoise = 0.001; % the noise in a neurons firing rate
SeqNoiseTime = [0,0]; % the noise in the sequence aka jitter (p of each neuron jittered 1 dt)
SeqNoiseNeuron = [1,1]; % the probability that a neuron participates in a given seq
T = 10000;
shared  = 1;
%%
V = data.sequences(T,Nneurons,Dt,NeuronNoise,SeqNoiseTime,SeqNoiseNeuron,shared);

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function CPD = root_CPD(bnet, self, val)
% ROOT_CPD Make a conditional prob. distrib. which has no parameters.
% CPD = ROOT_CPD(BNET, NODE_NUM, VAL)
%
% The node must not have any parents and is assumed to always be observed.
% It is a way of modelling exogenous inputs to a model.
% VAL is the value to which the root is clamped (default: [])

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = init_fields;
  CPD = class(CPD, 'root_CPD', generic_CPD(1));
  return;
elseif isa(bnet, 'root_CPD')
  % This might occur if we are copying an object.
  CPD = bnet;
  return;
end
CPD = init_fields;

if nargin < 3, val = []; end

ns = bnet.node_sizes;
ps = parents(bnet.dag, self);
if ~isempty(ps)
  error('root CPDs should have no parents')
end

CPD.self = self;
CPD.val = val;
CPD.sizes = ns(self);

clamped = 1;
CPD = class(CPD, 'root_CPD', generic_CPD(clamped));

%%%%%%%%%%%

function CPD = init_fields()
% This ensures we define the fields in the same order
% no matter whether we load an object from a file,
% or create it from scratch. (Matlab requires this.)

CPD.self = [];
CPD.val = [];
CPD.sizes = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Demo 2: With missing data

load data/faceSource.mat
load data/faceTargetMissing.mat

% Specify that surface normals are available and can be used.
Options.useNormals = 1;

% Specify that the source deformations should be plotted.
Options.plot = 1;

% Perform non-rigid ICP
[pointsTransformed, X] = nricp(Source, TargetMissing, Options);

author: charlienash

repo: nricp

save_path: github-repos/MATLAB/charlienash-nricp


text:
function [multiLineCols] = PTlinecmap(nColors)
%% [multiLineCols] = PTlinecmap(nColors)
% input number of color lines, output colormap

% ----------------------------------------------------------------------------------
% "THE BEER-WARE LICENSE" (Revision 42):
% <brian.white@queensu.ca> wrote this file. As long as you retain this notice you
% can do whatever you want with this stuff. If we meet some day, and you think
% this stuff is worth it, you can buy me a beer in return. -Brian White
% ----------------------------------------------------------------------------------

    cmap=flipud(colormap(jet));
    multiLineCols=(downsample(cmap,ceil(length(cmap)/nColors)));
    for i = find(multiLineCols(:,1) > .5 & multiLineCols(:,2) > .7 & multiLineCols(:,3) < .3)
        multiLineCols(i,:) = multiLineCols(i,:) * .78;
    end
    round(100/nColors)
    multiLineCols = repmat(multiLineCols, round(100/nColors),1);% repeats colormap to be 100 rows long

end

repo: PIDtoolbox


text:
function out = binsumcol(in, numinbin)
% Rebins image by summing numinbin bins together along columns of the
% 2D data.
% out = binsumcol(in, numinbin)

di=0;
if strcmp(class(in),'dip_image')
    di=1;
    in=double(in);
end
tmp = in;

for ii=2:numinbin
    tmp = tmp + circshift(in,[-1*(ii-1) 0]);
end

nb = floor(size(in,1)/numinbin);
ix=numinbin*(1:nb)-(numinbin-1);
out = tmp(ix,:);
if di
    out=dip_image(out);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% Fun with remap
% Maps a rectangular image into a circle using cv.remap
%
% To build the mapping we do the following simple geometry calculation:
%
% * shift the origin (0,0) coordinate to the center of the image
% * convert the cartesian x-y coordinates to polar r-theta form
% * determine the length of the line with the same theta that touches the
%   border of the rectangle
% * retain the theta, but scale the r value based on how much the line has to
%   shrink to fit into the circle.
% * convert the modified r-theta values back into x-y coordinates
%
% Sources:
%
% * <http://sidekick.windforwings.com/2012/12/opencv-fun-with-remap.html>
%

%%
% source image
fname = which('coins.png');
if isempty(fname)
    fname = fullfile(mexopencv.root(), 'test', 'apple.jpg');
end
src = imread(fname);

%%
% coordinates space, shifted so that origin is at the center of the image
[h,w,~] = size(src);
c = [w,h] / 2;
[X,Y] = meshgrid(single(1:w)-c(1), single(1:h)-c(2));

%%
% convert cartesian coordinates to polar form (r,theta)
R = hypot(Y,X);
T = atan(Y./X);  % NOTE: atan not atan2

% handle NaN case when atan(0/0) in the center
T(isnan(T)) = 0;

%%
% scale R
radius = min(h,w) / 2;
D = min(abs(c(1)./cos(T)), abs(c(2)./sin(T)));
R = R ./ (radius ./ D);

%%
% remap points
map_X = c(1) + sign(X) .* R .* cos(abs(T));
map_Y = c(2) + sign(Y) .* R .* sin(abs(T));

%%
% destination image
dst = cv.remap(src, map_X, map_Y);

%%
% show result
subplot(121), imshow(src), title('image')
subplot(122), imshow(dst), title('remapped')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function MimSaveAsMatlab(image_data, path, filename, reporting)
    % MimSaveAsMatlab. Saves an image as a Matlab matrix
    %
    %     Syntax
    %     ------
    %
    %         MimSaveAsMatlab(image_data, path, filename, reporting)
    %
    %             image_data      is a PTKImage (or PTKDicomImage) class containing the image
    %                             to be saved
    %             path, filename  specify the location to save the files.
    %             reporting       an object implementing CoreReportingInterface
    %                             for reporting progress and warnings
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if ~isa(image_data, 'PTKImage')
        reporting.Error('MimSaveAsMatlab:InputMustBePTKImage', 'Requires a PTKImage as input');
    end

    full_filename = fullfile(path, filename);
    value = [];
    value.image = image_data.RawImage;
    MimDiskUtilities.Save(full_filename, value);
end


text:
function result = RMS_amplitude(data)

if size(data,1) < size(data,2)
    data = data';
end
result = sqrt((data'*data)/(length(data)));

author: alexandrebarachant


text:
function [yd] = dm2yd(dm)
% Convert length from decimeters to yards.
% Chad A. Greene 2012
yd = dm*0.1093613298338;

author: Sable

repo: mcbench-benchmarks


text:
function [ cc_est, r_min, r_max ] = cube_random_centralize ( m, n, oc, cc, cr )

%*****************************************************************************80
%
%% CUBE_RANDOM_CENTRALIZE centralizes random data from a cube surface.
%
%  Discussion:
%
%    We generate N points uniformly at random on the surface of the
%    M dimensional cube, of center CC and radius R.
%
%    The cube is aligned with the coordinate axes.
%
%    We seek to estimate the value of CC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of random points to generate.
%
%    Input, real OC(M,1), the observation point.
%
%    Input, real CC(M,1), the center of the cube.
%
%    Input, real CR, the radius of the cube.
%
%    Output, real CC_EST(M,1), the estimated cube center.
%
%    Output, real R_MIN, R_MAX, the minimum and maximum distances
%    between the observation point and sample points on the cube surface.
%

%
%  Destroy all row vectors!
%
  cc = cc(:);
  oc = oc(:);
%
%  Compute sample points on the cube surface.
%
  x = cube_surface_sample ( m, n, cc, cr );
%
%  The center estimate is simply the average.
%
  cc_est(1:m,1) = sum ( x, 2 ) / n;
%
%  For each point on the cube surface, the "radius" ROC is the norm of X - OC.
%
  roc = x - repmat ( oc, 1, n );
  roc = roc.^2;
  roc = sum ( roc, 1 );
  roc = sqrt ( roc );
%
%  Return the minimum and maximum radiuses.
%
  r_min = min ( roc );
  r_max = max ( roc );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function truncated_svd_bfgs_test(A, p)

    % Generate some random data to test the function if none is given.
    if ~exist('A', 'var') || isempty(A)
        A = randn(402, 500);
    end
    if ~exist('p', 'var') || isempty(p)
        p = 5;
    end

    % Retrieve the size of the problem and make sure the requested
    % approximation rank is at most the maximum possible rank.
    [m, n] = size(A);
    assert(p <= min(m, n), 'p must be smaller than the smallest dimension of A.');

    % Define the cost and its derivatives on the Grassmann manifold

    tuple.U = grassmannfactory(m, p);
    tuple.V = grassmannfactory(n, p);

    M = productmanifold(tuple);

    problem.M = M;
    problem.cost  = @cost;
    problem.egrad = @egrad;
    problem.ehess = @ehess;

    % Cost function
    function f = cost(X)
        U = X.U;
        V = X.V;
        f = -.5*norm(U'*A*V, 'fro')^2;
    end

    % Euclidean gradient of the cost function
    function g = egrad(X)
        U = X.U;
        V = X.V;
        AV = A*V;
        AtU = A'*U;
        g.U = -AV*(AV'*U);
        g.V = -AtU*(AtU'*V);
    end

    % Euclidean Hessian of the cost function
    function h = ehess(X, H)
        U = X.U;
        V = X.V;
        Udot = H.U;
        Vdot = H.V;
        AV = A*V;
        AtU = A'*U;
        AVdot = A*Vdot;
        AtUdot = A'*Udot;
        h.U = -(AVdot*AV'*U + AV*AVdot'*U + AV*AV'*Udot);
        h.V = -(AtUdot*AtU'*V + AtU*AtUdot'*V + AtU*AtU'*Vdot);
    end

    conjugategradient(problem);
    %     trustregions(problem);

    problem.precon = preconBFGS(problem);
    problem.linesearch = @(x, xdot, storedb, key) 1;
    options.beta_type ='steep';
    conjugategradient(problem, [],options);
end

author: NicolasBoumal

repo: manopt


text:
function [tree,asgn] = vl_hikmeans(data,K,nleaves)
% VL_HIKMEANS  Hierachical integer K-means
%   [TREE,ASGN] = VL_HIKMEANS(DATA,K,NLEAVES) applies recursive integer
%   K-menas to cluster the data DATA, returing a structure TREE
%   representing the clusters and a vector ASGN with the data to
%   cluster assignments. The depth of the recursive partition is
%   computed so that at least NLEAVES are generated.
%
%   VL_HIKMEANS() is built on top of VL_IKMEANS() and requires the data to
%   be of class UINT8.
%
%   TREE is a structure representing the hierarchical clusters.  Each
%   node of the tree is also a structure with fields:
%
%   DEPTH::
%     Depth of the tree (only at the root node)
%
%   CENTERS::
%     K cluster centers
%
%   SUB::
%     Array of K node structures representing subtrees
%     (this field is missing at leaves).
%
%   ASGN is a matrix with one column per datum and height equal to the
%   depth of the tree. Each column encodes the branch of the tree that
%   correspond to each datum.
%
%   Example::
%     ASGN(:,7) = [1 5 3] means that the tree as depth equal to 3 and
%     that the datum X(:,7) corresponds to the branch
%     ROOT->SUB(1)->SUB(5)->SUB(3).
%
%   See also: VL_HIKMEANSPUSH(), VL_HIKMEANSHIST(), VL_IKMEANS(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: Cloud-CV

repo: object-proposals


text:
function [dat] = ft_fetch_data(data, varargin)

% FT_FETCH_DATA mimics the behavior of FT_READ_DATA, but for a FieldTrip
% raw data structure instead of a file on disk.
%
% Use as
%   [dat] = ft_fetch_data(data, ...)
%
% See also FT_READ_DATA, FT_FETCH_HEADER, FT_FETCH_EVENT

% Copyright (C) 2009-2013, Jan-Mathijs Schoffelen, Robert Oostenveld
% Copyright (C) 2008, Esther Meeuwissen
%
% 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$

% check whether input is data
skipcheckdata = ft_getopt(varargin, 'skipcheckdata');
if isempty(skipcheckdata) || skipcheckdata ~= 1
  data = ft_checkdata(data, 'datatype', 'raw', 'hassampleinfo', 'yes');
end

% get the options
if true
  p = inputParser;
  p.KeepUnmatched = true;
  addOptional(p, 'header', []);
  addOptional(p, 'begsample', []);
  addOptional(p, 'endsample', []);
  addOptional(p, 'chanindx', []);
  addOptional(p, 'allowoverlap', false);
  parse(p,varargin{:});
  hdr           = p.Results.header;
  begsample     = p.Results.begsample;
  endsample     = p.Results.endsample;
  chanindx      = p.Results.chanindx;
  allowoverlap  = p.Results.allowoverlap;
else
  hdr          = ft_getopt(varargin, 'header');
  begsample    = ft_getopt(varargin, 'begsample');
  endsample    = ft_getopt(varargin, 'endsample');
  chanindx     = ft_getopt(varargin, 'chanindx');
  allowoverlap = ft_getopt(varargin, 'allowoverlap', false);
end

allowoverlap = istrue(allowoverlap);

if isempty(hdr)
  hdr = ft_fetch_header(data);
end

if isempty(begsample)
  begsample = min(data.sampleinfo(:,1));
end

if isempty(endsample)
  endsample = max(data.sampleinfo(:,2));
end

if isempty(chanindx)
  chanindx = 1:hdr.nChans;
end

% get trial definition according to original data file
if isfield(data, 'sampleinfo')
  trl = data.sampleinfo;
else
  ft_error('data does not contain a consistent trial definition, fetching data is not possible');
end
trlnum = length(data.trial);

% start with the output data being all NaN
dat = nan(numel(chanindx), endsample-begsample+1);

if trlnum>1
  % original implementation, used when the input data has multiple trials

  trllen = trl(:,2) - trl(:,1) + 1;

  % check whether data.trial is consistent with trl
  if size(trl,1)~=length(data.trial)
    ft_error('trial definition is not internally consistent')
  elseif any(trllen~=(trl(:,2)-trl(:,1)+1))
    ft_error('trial definition is not internally consistent')
  end

  minchan = min(chanindx);
  maxchan = max(chanindx);
  if minchan<1 || maxchan>hdr.nChans
    ft_error('selected channels are not present in the data')
  end

  % determine for each sample in the data where it originates from
  if begsample<1
    buflen  = endsample;
  else
    buflen  = endsample - begsample + 1;
  end
  buffer    = zeros(1, buflen, 'int32');

  count     = buffer;
  trialnum  = buffer;
  samplenum = buffer;

  for trllop = 1:trlnum
    trlbeg = trl(trllop,1);
    trlend = trl(trllop,2);
    if trlbeg>endsample || trlend<begsample
      % skip this piece, it is not interesting because the requested range falls completely outside
    else
      trlbeg = trlbeg - begsample + 1;
      trlend = trlend - begsample + 1;
      offset = 0;
      overlap = 0;
      if trlbeg < 1
        offset = 1 - trlbeg;
        trlbeg = 1;
      end
      if trlend > buflen
        overlap = trlend - buflen;
        trlend = buflen;
      end
      % make vector with 0= no sample of old data, 1= one sample of old data, 2= two samples of old data, etc
      count(trlbeg:trlend) = count(trlbeg:trlend) + 1;
      % make vector with 1's if samples belong to trial 1, 2's if samples belong to trial 2 etc. overlap/ no data --> Nan
      trialnum(trlbeg:trlend) = trllop;
      % make samplenum vector with samplenrs for each sample in the old trials
      samplenum(trlbeg:trlend) = (1 + offset):(trllen(trllop) - overlap);
    end
  end

  % overlap --> NaN
  % trialnum(count>1)  = NaN;
  % samplenum(count>1) = NaN;

  % check if all samples are present and are not present twice or more
  if any(count>1)
    if ~allowoverlap
      % ft_error('some of the requested samples occur twice in the data');
      % this  can be considered OK if the overlap has exactly identical values
      sel = find(count>1); % must be row vector
      for smplop=sel
        % find in which trials the sample occurs
        seltrl = find(smplop>=trl(:,1) + 1 - begsample & ...
                      smplop<=trl(:,2) + 1 - begsample);  % which trials, requires the adjustment with begsample, if different from 1, JM 20180116
        selsmp = smplop - trl(seltrl,1) + begsample; % which sample in each of the trials, requires the adjustment with begsample, rather than 1
        for i=2:length(seltrl)
          % compare all occurences to the first one
          % consider also mutual occurring NaNs as equal values
          if ~all(isequaln(data.trial{seltrl(i)}(:,selsmp(i)), data.trial{seltrl(1)}(:,selsmp(1))))
            ft_error('some of the requested samples occur twice in the data and have conflicting values');
          end
        end
      end
    else
      ft_warning('samples present in multiple trials, using only the last occurence of each sample')
    end
  end

  %   if any(count==0)
  %     ft_warning('not all requested samples are present in the data, filling with NaNs');
  %   end

  % construct the output data array
  % dat = nan(length(chanindx), length(samplenum));
  % for smplop=1:length(samplenum)
  %   if samplenum(smplop)==0
  %    dat(:, smplop) = nan;
  %   else
  %    dat(:, smplop) = data.trial{trialnum(smplop)}(chanindx,samplenum(smplop));
  %   end
  % end

  % the following piece of code achieves the same as the commented code above,
  % but much smaller. rather than looping over samples it loops over the blocks
  % of samples defined by the original trials

  utrl = unique(trialnum);
  utrl(~isfinite(utrl)) = 0;
  utrl(utrl==0) = [];
  if length(utrl)==1
    ok   = trialnum==utrl;
    smps = samplenum(ok);
    dat(:,ok) = data.trial{utrl}(chanindx,smps);
  else
    for xlop=1:length(utrl)
      ok   = trialnum==utrl(xlop);
      smps = samplenum(ok);
      dat(:,ok) = data.trial{utrl(xlop)}(chanindx,smps);
    end
  end

else
  % only one trial is present in the input data, so it's quite simple and can be done much faster

  % get the indices
  begindx  = begsample - trl(1) + 1;
  endindx  = endsample - trl(1) + 1;
  if endindx<=0 || begindx>size(data.trial{1},2)
    % requested data samples outside the data present
    return;
  end

  tmptrl = trl([1 2]) - [trl(1) trl(1)]+1; % ignore offset in case it's present

  datbegindx = max(1,                     trl(1)-begsample+1);
  datendindx = min(endsample-begsample+1, trl(2)-begsample+1);

  %   if begsample<trl(1) || endsample>trl(2)
  %     ft_warning('not all requested samples are present in the data, filling with NaNs');
  %   end

  if begsample >= trl(1) && begsample <= trl(2)
    % the begin sample can be found in the available data
    if endsample >= trl(1) && endsample <= trl(2)
      dat(:,datbegindx:datendindx) = data.trial{1}(chanindx,begindx:endindx);
    else
      dat(:, datbegindx:datendindx) = data.trial{1}(chanindx,begindx:tmptrl(2));
    end
  elseif endsample >= trl(1) && endsample <= trl(2)
    % the end sample can be found in the available data
    dat(:, datbegindx:datendindx) = data.trial{1}(chanindx,tmptrl(1):endindx);
  else
    % neither the begin, nor the end sample are in the available data
    dat(:, datbegindx:datendindx) = data.trial{1}(chanindx,tmptrl(1):tmptrl(2));
  end

end % if trlnum is multiple or one

author: spm


text:
function [hdr,be] = read_hdr_raw(hname)
% Read a NIFTI header
% FORMAT [hdr,be] = read_hdr_raw(hname)
% hname - filename of image's header
% hdr   - a structure containing header info
% be    - whether big-endian or not
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: read_hdr_raw.m 6995 2017-01-27 15:57:40Z guillaume $

hdr = [];
be  = [];
sts = true;

%-Open header file
%--------------------------------------------------------------------------
fp  = fopen(hname,'r','native');
if fp==-1
    return;
end

%-Detect endianness and file format
%--------------------------------------------------------------------------
[unused,unused,mach] = fopen(fp);
if strncmp(mach,'ieee-be',7)
    be = true;
else
    be = false;
end

sb  = false;
fseek(fp,0,'bof');
d   = fread(fp,1,'*int32');
if isempty(d), fclose(fp); return; end
if d == 348
    fmt = 'nifti1';
elseif swapbytes(d) == 348
    fmt = 'nifti1';
    sb = true;
elseif d == 540
    fmt = 'nifti2';
elseif swapbytes(d) == 540
    fmt = 'nifti2';
    sb = true;
else
    %fclose(fp); return;
    warning('Cannot recognise format. Trying Analyze.');
    fmt = 'analyze';
end

%-Swap bytes if necessary
%--------------------------------------------------------------------------
if sb
    be = ~be;
    if be, mach = 'ieee-be';
    else   mach = 'ieee-le';
    end
    fclose(fp);
    fp = fopen(hname,'r',mach);
    if fp==-1
        return;
    end
end

%-Read magic string
%--------------------------------------------------------------------------
fseek(fp,0,'bof');
switch fmt
    case {'analyze','nifti1'}
        fseek(fp,344,'bof');
        mgc = fread(fp,4,'uint8')';
        if numel(mgc)~=4, fclose(fp); return; end
        switch char(mgc(1:3))
            case {'ni1','n+1'}
                org = niftistruc('nifti1');
            otherwise
                org = mayostruc;
        end
    case 'nifti2'
        fseek(fp,4,'bof');
        mgc = fread(fp,8,'uint8')';
        if numel(mgc)~=8 || ...
           ~ismember(char(mgc(1:3)),{'ni2','n+2'}) || ...
           ~isequal(mgc(5:8),[13 10 26 10]) % uint8(sprintf('\r\n\032\n'))
            fclose(fp);
            return;
        end
        org = niftistruc('nifti2');
end

%-Read header fields
%--------------------------------------------------------------------------
fseek(fp,0,'bof');
for i=1:length(org)
    field = fread(fp,org(i).len,['*' org(i).dtype.prec])';
    if numel(field) ~= org(i).len
        disp([length(field) org(i).len]);
        field = org(i).def;
        sts = false;
    end
    hdr.(org(i).label) = feval(org(i).dtype.conv,field);
end

%-Read header extensions
%--------------------------------------------------------------------------
try
    ext = read_extensions(fp,hdr);
    if ~isempty(ext)
        hdr.ext = ext;
    end
end

%-Close header file
%--------------------------------------------------------------------------
try, fclose(fp); catch, sts = false; end

%-Report problems
%--------------------------------------------------------------------------
if ~sts
    fprintf('There was a problem reading the header of\n');
    fprintf('"%s".\n', hname);
end

%==========================================================================
function ext = read_extensions(fp,hdr)

%-Read first header extension (if any)
%--------------------------------------------------------------------------
ext = []; n = 0;
exts = fread(fp,4,'*uint8');
if numel(exts) ~= 4, exts = uint8([0 0 0 0]); end
if exts(1) ~= 0
    n = n + 1;
    ext(n).esize = fread(fp,1,'*int32');
    ext(n).ecode = fread(fp,1,'*int32');
    ext(n).edata = fread(fp,ext(n).esize-8,'*uint8');
    if numel(ext(n).edata) ~= ext(n).esize-8
        fprintf('Cannot read header extension.');
        ext(n)   = [];
        n        = n - 1;
    else
        while ext(n).edata(end) == 0, ext(n).edata(end) = []; end
    end
end


text:
function [M, Mx, My, Mz] = convectionTerm3D(u)
% This function uses the central difference scheme to discretize a 2D
% convection term in the form \grad (u \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%   [M, Mx, My, Mz] = convectionTerm3D(u)
%
% PARAMETERS:
%   u   - FaceVariable
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nx = u.domain.dims(1);
Ny = u.domain.dims(2);
Nz = u.domain.dims(3);
G=reshape((1:(Nx+2)*(Ny+2)*(Nz+2)), Nx+2, Ny+2, Nz+2);
DXe = repmat(u.domain.cellsize.x(3:end), 1, Ny, Nz);
DXw = repmat(u.domain.cellsize.x(1:end-2), 1, Ny, Nz);
DXp = repmat(u.domain.cellsize.x(2:end-1), 1, Ny, Nz);
DYn = repmat(u.domain.cellsize.y(3:end)', Nx, 1, Nz);
DYs = repmat(u.domain.cellsize.y(1:end-2)', Nx, 1, Nz);
DYp = repmat(u.domain.cellsize.y(2:end-1)', Nx, 1, Nz);
DZ = zeros(1,1,Nz+2);
DZ(1,1,:) = u.domain.cellsize.z;
DZf=repmat(DZ(1,1,3:end), Nx, Ny, 1);
DZb=repmat(DZ(1,1,1:end-2), Nx, Ny, 1);
DZp=repmat(DZ(1,1,2:end-1), Nx, Ny, 1);

% define the vectors to stores the sparse matrix data
iix = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
jjx = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
sx = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
iiy = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
jjy = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
sy = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
iiz = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
jjz = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
sz = zeros(3*(Nx+2)*(Ny+2)*(Nz+2),1);
mnx = Nx*Ny*Nz;	mny = Nx*Ny*Nz;   mnz = Nx*Ny*Nz;

% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = u.xvalue(2:Nx+1,:,:)./(DXp+DXe);
uw = u.xvalue(1:Nx,:,:)./(DXp+DXw);
vn = u.yvalue(:,2:Ny+1,:)./(DYp+DYn);
vs = u.yvalue(:,1:Ny,:)./(DYp+DYs);
wf = u.zvalue(:,:,2:Nz+1)./(DZp+DZf);
wb = u.zvalue(:,:,1:Nz)./(DZp+DZb);

% calculate the coefficients for the internal cells
AE = reshape(ue,mnx,1);
AW = reshape(-uw,mnx,1);
AN = reshape(vn,mny,1);
AS = reshape(-vs,mny,1);
AF = reshape(wf,mnz,1);
AB = reshape(-wb,mnz,1);
APx = reshape((ue.*DXe-uw.*DXw)./DXp,mnx,1);
APy = reshape((vn.*DYn-vs.*DYs)./DYp,mny,1);
APz = reshape((wf.*DZf-wb.*DZb)./DZp,mnz,1);

% build the sparse matrix based on the numbering system
rowx_index = reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mnx,1); % main diagonal x
iix(1:3*mnx) = repmat(rowx_index,3,1);
rowy_index = reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mny,1); % main diagonal y
iiy(1:3*mny) = repmat(rowy_index,3,1);
rowz_index = reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mnz,1); % main diagonal z
iiz(1:3*mnz) = repmat(rowz_index,3,1);
jjx(1:3*mnx) = [reshape(G(1:Nx,2:Ny+1,2:Nz+1),mnx,1); reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mnx,1); reshape(G(3:Nx+2,2:Ny+1,2:Nz+1),mnx,1)];
jjy(1:3*mny) = [reshape(G(2:Nx+1,1:Ny,2:Nz+1),mny,1); reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mny,1); reshape(G(2:Nx+1,3:Ny+2,2:Nz+1),mny,1)];
jjz(1:3*mnz) = [reshape(G(2:Nx+1,2:Ny+1,1:Nz),mnz,1); reshape(G(2:Nx+1,2:Ny+1,2:Nz+1),mnz,1); reshape(G(2:Nx+1,2:Ny+1,3:Nz+2),mnz,1)];
sx(1:3*mnx) = [AW; APx; AE];
sy(1:3*mny) = [AS; APy; AN];
sz(1:3*mnz) = [AB; APz; AF];

% build the sparse matrix
kx = 3*mnx;
ky = 3*mny;
kz = 3*mnz;
Mx = sparse(iix(1:kx), jjx(1:kx), sx(1:kx), (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2));
My = sparse(iiy(1:ky), jjy(1:ky), sy(1:ky), (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2));
Mz = sparse(iiz(1:kz), jjz(1:kz), sz(1:kz), (Nx+2)*(Ny+2)*(Nz+2), (Nx+2)*(Ny+2)*(Nz+2));
M = Mx + My + Mz;

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool