text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% write_gipl_volume.m
% Function for writing Guys Image Processing Lab (Gipl) files with given
% volumn, file name of the output, and voxel size.
%
% for example,
% I=uint8(rand(64,64,64)*256);
% 1: gipl_write_volume(I);
% 2: gipl_write_volume(I,'random.gipl',[2 2 2];

% Wan Jing
% 11/06/2008 - create
%============================================

function gipl_write_volume(I,fname,scales,origin)

% define supported image types and bit length.
image_types=struct('double',65,'float',65,'single',64,'logical',1,'int8',7,'uint8',8,'int16',15,'uint16',16,'int32',32,'uint32',31);
bit_lengths=struct('double',64,'float',64,'single',32,'logical',1,'int8',8,'uint8',8,'int16',16,'uint16',16,'int32',32,'uint32',32);

% Sizes
    sizes=size(I);
    while(length(sizes)<4), sizes=[sizes 1]; end
% Scales
    if(exist('scales','var')==0), scales=ones(1,length(sizes)); end;
    while(length(scales)<4), scales=[scales 0]; end
% Offset
    offset=256;
% Image Type
    image_type=getfield(image_types,class(I));
% File size
    fsize=offset+prod(sizes)*getfield(bit_lengths,class(I))/8;
% Patient
    patient='Generated by Matlab';
    while(length(patient)<80), patient=[patient ' ']; end
% Matrix
    matrix=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% orientation
    orientation=0;
% voxel_min
    voxmin=min(I(:));
% voxel_max
    voxmax=max(I(:));
% origin
 %   origin=[0 0 0 0];
% pixval_offset
    pixval_offset=0;
% pixval_cal
    pixval_cal=0;
% interslicegap
    interslicegap=0;
% user_def2
    user_def2=0;
% par2
    par2=0;
% magic_number
    magic_number=4026526128;

trans_type{1}='binary'; trans_type{7}='char'; trans_type{8}='uchar'; trans_type{15}='short';
trans_type{16}='ushort'; trans_type{31}='uint'; trans_type{32}='int'; trans_type{64}='float';
trans_type{65}='double'; trans_type{144}='C_short'; trans_type{160}='C_int'; trans_type{192}='C_float';
trans_type{193}='C_double'; trans_type{200}='surface'; trans_type{201}='polygon';

trans_orien{0+1}='UNDEFINED'; trans_orien{1+1}='UNDEFINED_PROJECTION'; trans_orien{2+1}='AP_PROJECTION';
trans_orien{3+1}='LATERAL_PROJECTION'; trans_orien{4+1}='OBLIQUE_PROJECTION';  trans_orien{8+1}='UNDEFINED_TOMO';
trans_orien{9+1}='AXIAL'; trans_orien{10+1}='CORONAL'; trans_orien{11+1}='SAGITTAL'; trans_orien{12+1}='OBLIQUE_TOMO';

% display information about the volumn
disp(['filename : ' num2str(fname)]);
disp(['filesize : ' num2str(fsize)]);
disp(['sizes : ' num2str(sizes)]);
disp(['scales : ' num2str(scales)]);
disp(['image_type : ' num2str(image_type) ' - ' trans_type{image_type}]);
disp(['patient : ' patient]);
disp(['matrix : ' num2str(matrix)]);
disp(['orientation : ' num2str(orientation) ' - ' trans_orien{orientation+1}]);
disp(['voxel min : ' num2str(voxmin)]);
disp(['voxel max : ' num2str(voxmax)]);
disp(['origin : ' num2str(origin)]);
disp(['pixval_offset : ' num2str(pixval_offset)]);
disp(['pixval_cal : ' num2str(pixval_cal)]);
disp(['interslicegap : ' num2str(interslicegap)]);
disp(['user_def2 : ' num2str(user_def2)]);
disp(['par2 : ' num2str(par2)]);
disp(['offset : ' num2str(offset)]);

% save volumn to *.gipl file
fout=fopen(fname,'wb','ieee-be');
fwrite(fout,uint16(sizes),'ushort'); % 4
fwrite(fout,uint16(image_type),'ushort'); % 1
fwrite(fout,single(scales),'float'); % 4
fwrite(fout,patient,'char'); % 80
fwrite(fout,single(matrix),'float'); % 20
fwrite(fout,uint8(orientation),'uint8'); % 1
fwrite(fout,uint8(par2),'uint8'); % 1
fwrite(fout,double(voxmin),'double'); % 1
fwrite(fout,double(voxmax),'double'); % 1
fwrite(fout,double(origin),'double'); % 4
fwrite(fout,single(pixval_offset),'float'); % 1
fwrite(fout,single(pixval_cal),'float'); % 1
fwrite(fout,single(interslicegap),'float'); % 1
fwrite(fout,single(user_def2),'float'); % 1
fwrite(fout,uint32(magic_number),'uint'); % 1
fwrite(fout,I, trans_type{image_type});
fclose('all');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%
% This demo shows the designed of filters and some figures in
%
% Ha T. Nguyen and Minh N. Do, Hybrid Filter Banks with Fractional Delays:
% Minimax Design and Application to Multichannel Sampling, vol. 56, no. 7,
% pp. 3180-3190, July 2008.

initialization;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Design of (IIR for FIR) synthesis filters F_i(z)

% Design IIR synthesis filters F_i(z)
[num, den, gamma] = designIIR(phi, D, m0, h, M);

% Design FIR synthesis filters F_i(z)
% [num, den, gamma] = designFIR(phi, D, m0, h, M, n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demo of the approximation of the high-resolution signal
% Filtering using the synthesis filters F_i(z)
y0hat = 0*y0;

for i = 1:N
    z = filter(num{i}, den{i}, upsample(x{i},M));
    y0hat = y0hat + z(1:L);
end

% Approximation error
e = y0 - y0hat;

% Compute the errors
max(e)
mean(e(1:200).^2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results,

% Plot Fig. 9
figure(9);
bode(numphi, denphi);
title('Bode diagram of \Phi(s)')

% Plot Fig. 10
plot_equi_filters;

% Plot of Fig. 11
% (plot only the first 2 filters if more than 2 are available)
figure(11);
freqz(num{1}, den{1});
hold on;
freqz(num{2}, den{2});

% figure;
% plot(e(1:200));
% xlabel('sample');
% ylabel('error');
% title('Approximation error')

% Plot Fig. 12
figure(12);
plot(e(1:50));
hold on;
plot(y0(1:50), 'r--')
xlabel('sample');
ylabel('error');
legend('error', 'signal')
title('Approximation error vs the high resolution signal')

author: Sable

repo: mcbench-benchmarks


text:
% PURPOSE: detects ICA process in a dataset
%
% FORMAT
%
% value = isica(EEG);
%
% if value = 1 then the EEG dataset is ICA-ed. Otherwise, value = 0
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% July 7, 2011

function value = isica(EEG)
value     = 0; % no ica-ed file by default
if nargin<1
        help isica
        return
end
icafields = {'icaact','icawinv','icasphere','icaweights','icachansind','specicaact','icasplinefile'};
nif       = length(icafields);
for k=1:nif
        if isfield(EEG, icafields{k}) && ~isempty(EEG.(icafields{k}))
                value = 1; % ica-ed file detected
                break
        end
end
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% Example for the algorithms described in
%
%       D. Kressner, M. Steinlechner, A. Uschmajew:
%		Low-rank tensor methods with subspace correction for symmetric eigenvalue problems
%		SIAM J. Sci. Comput., 36(5):A2346-A2368, 2014.
%
% Code to produce Figure 4.4: Three Eigenvalues for Henon-Heiles with
% with n=28, d=10

%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt
% =========================================================================
clear all
close all

n = 28;
d = 10;
A = TTeMPS_op_NN_hermite(n, d);
p = 3;
r = 1;

% Run block eigenvalue procedure:
% =========================================================================

if ~exist('hh_3_blk_hermite.mat','file')
    opts = struct( 'maxiter', 3, ...
                   'maxrank', 40, ...
                   'tol', 1e-8, ...
                   'tolOP', 1e-3, ...
                   'tolLOBPCG', 1e-6, ...
                   'maxiterLOBPCG', 500, ...
                   'verbose', true , ...
                   'precInner', true);

    rng(11)
    rr = [1, 1 * ones(1, d-1), 1];
    [X_blk_hermite, C_blk_hermite, evalue_blk_hermite, residuums_blk_hermite, micro_res_blk_hermite, objective_blk_hermite, t_blk_hermite] = block_eigenvalue( A, p, rr, opts);
    save('hh_3_blk_hermite', 'X_blk_hermite', 'C_blk_hermite', 'evalue_blk_hermite', 'residuums_blk_hermite', 'micro_res_blk_hermite', 'objective_blk_hermite','t_blk_hermite');
else
    load('hh_3_blk_hermite.mat')
end

% Run EVAMEn:
% =========================================================================

if ~exist('hh_3_evamen_hermite.mat','file')
    opts = struct( 'maxiter', 3, ...
                   'maxrank', 40, ...
                   'maxrankRes', 2, ...
                   'tol', 1e-8, ...
                   'tolOP', 1e-3, ...
                   'tolLOBPCG', 1e-6, ...
                   'maxiterLOBPCG', 500, ...
                   'verbose', true , ...
                   'precInner', true);
    rng(11)
    rr = [1, 1 * ones(1, d-1), 1];
    [X_evamen_hermite, C_evamen_hermite, evalue_evamen_hermite, residuums_evamen_hermite, micro_res_evamen_hermite, objective_evamen_hermite, t_evamen_hermite] = amen_eigenvalue( A, 1, p, rr, opts);
    save('hh_3_evamen_hermite', 'X_evamen_hermite', 'C_evamen_hermite', 'evalue_evamen_hermite', 'residuums_evamen_hermite', 'micro_res_evamen_hermite', 'objective_evamen_hermite','t_evamen_hermite');
else
    load('hh_3_evamen_hermite.mat')
end

% Prepare data for plotting:
% =========================================================================

evalue_end = repmat(evalue_blk_hermite(:,end), [1,size(evalue_blk_hermite,2)-1]);
ev_blk_hermite = abs(evalue_blk_hermite(:,1:end-1) - evalue_end);
ev_evamen_hermite = abs(evalue_evamen_hermite(:,1:end-1) - evalue_end);

% Plot vs. Iterations
% =========================================================================
f = figure
set(0,'defaultlinelinewidth',2)
subplot(1,2,1)

semilogy( sqrt(sum(micro_res_blk_hermite.^2, 1)), '-b' )
hold on
semilogy( sqrt(sum(micro_res_evamen_hermite.^2, 1)), '-k' )

semilogy( sum(ev_blk_hermite,1), '--b' )
semilogy( sum(ev_evamen_hermite,1), '--k' )

res_blk_hermite         = sqrt(sum(micro_res_blk_hermite.^2, 1))
res_evamen_hermite  = sqrt(sum(micro_res_evamen_hermite.^2, 1))

semilogy((d-1):(d-1):length(micro_res_blk_hermite),res_blk_hermite(:,(d-1):(d-1):end),'ob')
semilogy((d-1):(d-1):length(micro_res_evamen_hermite),res_evamen_hermite(:,(d-1):(d-1):end),'ok')

semilogy((d-1):(d-1):length(ev_blk_hermite),sum(ev_blk_hermite(:,(d-1):(d-1):end),1),'ob')
semilogy((d-1):(d-1):length(ev_evamen_hermite),sum(ev_evamen_hermite(:,(d-1):(d-1):end),1),'ok')

set(gca,'fontsize',20)
xlabel('Microiterations')
ylabel('Residual and eigenvalue error')

% Plot vs. Time
% =========================================================================

subplot(1,2,2)
semilogy( t_blk_hermite, sqrt(sum(micro_res_blk_hermite.^2, 1)), '-b' )
hold on
semilogy( t_evamen_hermite, sqrt(sum(micro_res_evamen_hermite.^2, 1)), '-k' )

semilogy( t_blk_hermite, sum(ev_blk_hermite,1), '--b' )
semilogy( t_evamen_hermite, sum(ev_evamen_hermite,1), '--k' )

semilogy(t_blk_hermite((d-1):(d-1):end),res_blk_hermite(:,(d-1):(d-1):end),'ob')
semilogy(t_evamen_hermite((d-1):(d-1):end),res_evamen_hermite(:,(d-1):(d-1):end),'ok')

semilogy(t_blk_hermite((d-1):(d-1):end),sum(ev_blk_hermite(:,(d-1):(d-1):end),1),'ob')
semilogy(t_evamen_hermite((d-1):(d-1):end),sum(ev_evamen_hermite(:,(d-1):(d-1):end),1),'ok')

semilogy(t_blk_hermite((d-1):(d-1):end),        res_blk_hermite(:,(d-1):(d-1):end),'ob')
semilogy(t_evamen_hermite((d-1):(d-1):end), res_evamen_hermite(:,(d-1):(d-1):end),'ok')

semilogy(t_blk_hermite((d-1):(d-1):end),        sum(ev_blk_hermite(:,(d-1):(d-1):end),1),'ob')
semilogy(t_evamen_hermite((d-1):(d-1):end), sum(ev_evamen_hermite(:,(d-1):(d-1):end),1),'ok')

h_leg = legend('Res. err., Block-ALS',...
       'Res. err. EVAMEn, prec.',...
       'EV. err., Block-ALS',...
       'EV. err. EVAMEn, prec.')
set(gca,'fontsize',20)
set(h_leg, 'fontsize',16)
xlabel('Time [s]')
ylabel('Residual and eigenvalue error')

set(f, 'Position', [0 0 1200 700])

author: NicolasBoumal

repo: manopt


text:
clear
clc

fitnessfcn = @my_first_multi;   % Function handle to the fitness function
nvars = 2;                      % Number of decision variables
lb = [-5,-5];                   % Lower bound
ub = [5,5];                     % Upper bound
A = []; b = [];                 % No linear inequality constraints
Aeq = []; beq = [];             % No linear equality constraints
options = gaoptimset('ParetoFraction',0.3,'PopulationSize',100,'Generations',200,'StallGenLimit',200,'TolFun',1e-100,'PlotFcns',@gaplotpareto);

[x,fval] = gamultiobj(fitnessfcn,nvars, A,b,Aeq,beq,lb,ub,options);

author: HuangCongQing


text:
function [gKern, gVarmeans, gVarcovars, gInd] = rbfardjitVardistPsi1Gradient(rbfard2Kern, vardist, Z, covGrad, learnInducing)

% RBFARDJITVARDISTPSI1GRADIENT description.

% VARGPLVM

if nargin < 5
    learnInducing = 1;
end

 [gKern, gVarmeans, gVarcovars, gInd] = rbfard2VardistPsi1Gradient(rbfard2Kern, vardist, Z, covGrad, learnInducing);

author: SheffieldML

repo: vargplvm


text:
function [tests,pass,perf]=test_sumann(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
test_string={'report=sumann(''mitdb/100'',''atr'');'};
clean_up={};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function [ c, seed ] = ch_uniform ( clo, chi, seed )

%*****************************************************************************80
%
%% CH_UNIFORM returns a random character in a given range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CLO, CHI, the minimum and maximum acceptable characters.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, character C, the randomly chosen character.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ d, seed ] = r8_uniform_01 ( seed );

  c = round ( ( 1.0 - d ) * clo + d * chi );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, x, fx ] = psi_values ( n_data )

%*****************************************************************************80
%
%% PSI_VALUES returns some values of the Psi or Digamma function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      PolyGamma[x]
%
%    or
%
%      PolyGamma[0,x]
%
%    PSI(X) = d ln ( Gamma ( X ) ) / d X = Gamma'(X) / Gamma(X)
%
%    PSI(1) = -Euler's constant.
%
%    PSI(X+1) = PSI(X) + 1 / X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 11;

  fx_vec = [ ...
     -0.5772156649015329E+00, ...
     -0.4237549404110768E+00, ...
     -0.2890398965921883E+00, ...
     -0.1691908888667997E+00, ...
     -0.6138454458511615E-01, ...
      0.3648997397857652E-01, ...
      0.1260474527734763E+00, ...
      0.2085478748734940E+00, ...
      0.2849914332938615E+00, ...
      0.3561841611640597E+00, ...
      0.4227843350984671E+00 ];

  x_vec = [ ...
     1.0E+00, ...
     1.1E+00, ...
     1.2E+00, ...
     1.3E+00, ...
     1.4E+00, ...
     1.5E+00, ...
     1.6E+00, ...
     1.7E+00, ...
     1.8E+00, ...
     1.9E+00, ...
     2.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
% Prepares DSS calculation environment. Currently only compiles the
% testkeypress MEX function.

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

disp('Initializing DSS environment');
disp('Compiling testkeypress executable');

mex src/testkeypress.c;

disp('Environment initialized successfully');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function show3DLD(X, adj, color, varargin)

% function show3DLD(X, adj, color, varargin)
%
% This file is created by Tianfan Xue (tfxue@mit.edu).

para.showNumber = false;
para.showAxis = true;
para.fontsize = 20;
para = propval(varargin, para);

if (nargin == 2) || isempty(color)
    color = [1,0,0];
end

X1 = X(1,:); X2 = X(2,:); X3 = X(3,:);
plot3(X1,X2,X3,'.','Color',color,'MarkerSize',10);
hold on;
line(X1(adj)', X2(adj)', X3(adj)', 'Color', color,'LineWidth',3);
if para.showNumber
    for i=1:size(X,2)
        text(X1(i),X2(i),X3(i),sprintf('%d',i),'FontSize',para.fontsize);
    end
end
if para.showAxis
    xlabel('x');
    ylabel('y');
    zlabel('z');
end
axis equal

end

author: jiajunwu


text:
clc
clear all
close all

% To run this script we recommend running example2 first to create all the
% necessary outpouts from GLMsingle that are going to be reused here.

% This script shows how to save GLMsingle outputs as NIFTI files. NIFTI is
% neuroimaging format often used to store neuroimaging data. This script
% requires a NiFTi toolbox
% https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
% or alternatively one can use matlab's native support for nifti files
% (explaiend at the end of the script).

% load GLMsingle results
load('./example2outputs/GLMsingle/TYPED_FITHRF_GLMDENOISE_RR.mat')

% we will save following estiamtes (modelmd corresponds to single trial betas)
val2save = {'modelmd';'R2';'HRFindex';'FRACvalue'};
niftidir = './example2outputs/GLMsingle/GLMDENOISE_RR_nifti_files';
mkdir(niftidir);

voxelsize = [1.8 1.8 1.8];
origin = [0 0 0];

for v = 1 : length(val2save)

    nii = make_nii(eval(val2save{v}),voxelsize,origin);
    save_nii(nii,sprintf('%s/%s.nii',niftidir,val2save{v}))

end

% if you don't want to change the reference frame and position of the data
% (done by setting a new origin) you can reuse a NIFTI file (e.g T1 image)
% as a template for saving your results. This is usefull if your functional
% data is already aligned with the anatomy.

% template = load_nii('T1.nii')
% template is a matlab structure where "vol" field stores the data matrix.
% template.vol = eval(val2save{v}); swap the data to your GLMsingle output
% save_nii(template,sprintf('%s/%s.nii',niftidir,val2save{v})); save nifti

% additionally you can use matlab's native support for reading and saving
% nifti files (niftiread and niftiwrite). This skips making the nifti file
% and the header information can be directly added to niftiwrite function.
% info.Description = 'Modified using MATLAB R2017b';
% niftiwrite(nii,sprintf('%s/%s.nii',niftidir,val2save{v})), info);
% where info is a header.

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function clOut = get_cluster_volume(clIn)
% :Usage:
% ::
%
%    clOut = get_cluster_volume(clIn)
%
% Adapted from:
%
% FORMAT clusters = ihb_getClusters
%
% Get cluster information (use [SPM,VOL,xX,xCon,xSDM] = spm_getSPM;)
%
% ..
%    01.03.01    Sergey Pakhomov
%    01.08.01    last modified
%
%    Only change is to eliminate stuff in beginning - just update the cluster.
%    Tor Wager, 10/18/01
% ..

clOut = clIn;

%----------------------------------------------------------------------------------
% Define default voxel size for canonical image
% This is OK to leave at 1's, because the head plotting utility adjusts
% to millimeters.
%----------------------------------------------------------------------------------
ihbdfl_tal_x_vox = 1;
ihbdfl_tal_y_vox = 1;
ihbdfl_tal_z_vox = 1;

%----------------------------------------------------------------------------------
% Prepare MNI volume data xMNI, yMNI, zMNI, and volMNI
%----------------------------------------------------------------------------------
rMin = min(clOut.XYZmm');
rMax = max(clOut.XYZmm');
voxMNI = clOut.voxSize;

% make all voxMNI positive: Correct bug with negative x voxel size
voxMNI = abs(voxMNI);

[xMNI, yMNI, zMNI] = meshgrid(...
    rMin(1) - voxMNI(1) : voxMNI(1) : rMax(1) + voxMNI(1),...
    rMin(2) - voxMNI(2) : voxMNI(2) : rMax(2) + voxMNI(2),...
    rMin(3) - voxMNI(3) : voxMNI(3) : rMax(3) + voxMNI(3));

volMNI = zeros(size(xMNI));

[x y z] = size(xMNI);

for curVox = 1:clOut.numVox
    tmp = clOut.XYZmm(:, curVox);
    ind1 = (tmp(1) - rMin(1))/voxMNI(1) + 1;
    ind2 = (tmp(2) - rMin(2))/voxMNI(2) + 1;
    ind3 = (tmp(3) - rMin(3))/voxMNI(3) + 1;
    index = round(x*y*ind3 + x*ind1 + ind2 + 1);  % tor changed to avoid rounding error sometimes - 11/13
    volMNI(index) = clOut.Z(curVox);

    % waitbar stuff too much of a pain.
    %nTotalProcessed = nTotalProcessed + 1;
    %waitbar(nTotalProcessed/nTotalVox);
end

sizeMNI = size(volMNI);
sizeMNIprod = prod(sizeMNI);
xMNIreshaped = reshape(xMNI, [1 sizeMNIprod]);
yMNIreshaped = reshape(yMNI, [1 sizeMNIprod]);
zMNIreshaped = reshape(zMNI, [1 sizeMNIprod]);
xyzReshaped = [xMNIreshaped; yMNIreshaped; zMNIreshaped];

% tor eliminated this transformation line to keep in MNI space
%xyzTal = ihb_Mni2Tal(xyzReshaped);

xyzTal = xyzReshaped;

xx = reshape(xyzTal(1, :), sizeMNI);
yy = reshape(xyzTal(2, :), sizeMNI);
zz = reshape(xyzTal(3, :), sizeMNI);

clOut.xMin = ihbdfl_tal_x_vox*floor(min(xyzTal(1, :))/ihbdfl_tal_x_vox);
clOut.yMin = ihbdfl_tal_y_vox*floor(min(xyzTal(2, :))/ihbdfl_tal_y_vox);
clOut.zMin = ihbdfl_tal_z_vox*floor(min(xyzTal(3, :))/ihbdfl_tal_z_vox);

clOut.xMax = ihbdfl_tal_x_vox*ceil(max(xyzTal(1, :))/ihbdfl_tal_x_vox);
clOut.yMax = ihbdfl_tal_y_vox*ceil(max(xyzTal(2, :))/ihbdfl_tal_y_vox);
clOut.zMax = ihbdfl_tal_z_vox*ceil(max(xyzTal(3, :))/ihbdfl_tal_z_vox);

[clOut.xTal, clOut.yTal, clOut.zTal] = meshgrid(clOut.xMin:ihbdfl_tal_x_vox:clOut.xMax,...
                                                clOut.yMin:ihbdfl_tal_y_vox:clOut.yMax,...
                                                clOut.zMin:ihbdfl_tal_z_vox:clOut.zMax);
volTal = interp3(xx, yy, zz, volMNI, clOut.xTal, clOut.yTal, clOut.zTal, 'linear');
index = find(isnan(volTal(:,:,:)));
volTal(index) = 0;
clOut.vTal = volTal;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%SETPOSTPROC (Re)set POSTPROC field of a datafile
%
%   A = SETPOSTPROC(A,MAPPING)
%   A = SETPOSTPROC(A)
%
% INPUT
%   A        - Datafile
%   POSTPROC - Postprocessing mapping command
%
% OUTPUT
%   A       - Datafile
%
% DESCRIPTION
% Sets the mappings stored in A.POSTPROC. The size of the datafile
% A is set to the output size of PRMAPPING.
% A call without MAPPING clears A.POSTPROC. The size of the datafile
% A is reset to undefined (0).
%
% The mappings in A.POSTPROC may be extended by ADDPOSTPROC.
%
% Mappings in A.POSTPROC are stored only and executed just
% after A is converted from a DATAFILE into a PRDATASET.
%
% SEE ALSO (PRTools Guide)
% DATAFILES, SETPREPROC, ADDPOSTPROC.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [x,P,K] = KalmanFilterX_UpdatePDA(xPred,PPred,Y,W,yPred,S,Pxy)
% KALMANFILTERX_UPDATEPDA Perform the discrete-time Probabilistic Data
% Association (PDA) KF update step, under the assumption of additive process
% noise, for multiple measurements (as a Gaussian Mixture)
%
% Parameters
% ----------
% xPred: column vector
%   The (xDim x 1) predicted state estimate.
% PPred: matrix
%   The (xDim x xDim) predicted state covariance matrix.
% Y: matrix
%   The (yDim x nY) measurement vector.
% W: row vector
%   The (1 x nY+1) measurement association/mixture weights
%   vector. (dummy measurement assumed at index 1)
% yPred: column vector
%   The (yDim x 1) predicted measurement estimate.
% S: matrix
%   The (yDim x yDim) innovation covariance matrix.
% Pxy: matrix
%   The (xDim x yDim) cross-covariance matrix.
%
% Returns
% -------
% x: column vector
%   The (xDim x 1) state estimate at the current time-step.
% P: matrix
%   The (xDim x xDim) state covariance matrix at the current
%   time-step.
% K: matrix
%   The (xDim x yDim) Kalman gain matrix at the current
%   time-step.
%
%October 2017 Lyudmil Vladimirov, University of Liverpool.

    % Get size of observation vector
    [yDim,nY] = size(Y);

    % Compute Kalman gain
    K = Pxy/S;

    % Compute innovation mean and (cross) covariance
    innov_err       = Y - yPred(:,ones(1,nY));
    tot_innov_err   = innov_err*W(2:end)';
    Pc              = PPred - K*S*K';
    Pgag            = K*((innov_err.*W(ones(yDim,1),2:end))*innov_err' - tot_innov_err*tot_innov_err')*K';

    % Compute filtered estimates
    x    = xPred + K*tot_innov_err;
    P    = W(1)*PPred + (1-W(1))*Pc + Pgag;
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
%MDL_PLANAR3 Create model of a simple planar 3-link mechanism
%
% MDL_PLANAR2 is a script that creates the workspace variable p3 which
% describes the kinematic characteristics of a simple redundant planar
% 3-link mechanism.
%
% Also defines the vector:
%   qz   corresponds to the zero joint angle configuration.
%
%
% Notes::
% - Moves in the XY plane.
% - No dynamics in this model.
%
% See also SerialLink, mdl_twolink, mdl_planar1, mdl_planar2.

% MODEL: generic, planar, 3DOF, standard_DH

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

a1 = 1;
a2 = 1;
a3 = 1;

p3 = SerialLink([
    Revolute('d', 0, 'a', a1, 'alpha', 0, 'standard')
    Revolute('d', 0, 'a', a2, 'alpha', 0, 'standard')
    Revolute('d', 0, 'a', a3, 'alpha', 0, 'standard')
    ], ...
    'name', 'three link');
qz = [0 0 0];

author: Allopart

repo: rbpf-gmapping


text:
function [ x, w ] = square_symq ( degree, n )

%*****************************************************************************80
%
%% SQUARE_SYMQ returns a symmetric quadrature rule for the square.
%
%  Discussion:
%
%    This procedure returns a quadrature rule for smooth functions
%    on the unit square [-1,1]^2.
%
%    All quadratures are accurate to 15 digits
%    All weights are positive and inside the square
%
%    The nodes are symmetrically arranged.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the degree of the quadrature rule.
%    1 <= DEGREE <= 20.
%
%    Input, integer N, the number of nodes.
%    This can be determined by a call to RULE_FULL_SIZE(DEGREE).
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  if ( degree == 1 )
    [ x, w ] = rule01 ( n );
  elseif ( degree == 2 )
    [ x, w ] = rule02 ( n );
  elseif ( degree == 3 )
    [ x, w ] = rule03 ( n );
  elseif ( degree == 4 )
    [ x, w ] = rule04 ( n );
  elseif ( degree == 5 )
    [ x, w ] = rule05 ( n );
  elseif ( degree == 6 )
    [ x, w ] = rule06 ( n );
  elseif ( degree == 7 )
    [ x, w ] = rule07 ( n );
  elseif ( degree == 8 )
    [ x, w ] = rule08 ( n );
  elseif ( degree == 9 )
    [ x, w ] = rule09 ( n );
  elseif ( degree == 10 )
    [ x, w ] = rule10 ( n );
  elseif ( degree == 11 )
    [ x, w ] = rule11 ( n );
  elseif ( degree == 12 )
    [ x, w ] = rule12 ( n );
  elseif ( degree == 13 )
    [ x, w ] = rule13 ( n );
  elseif ( degree == 14 )
    [ x, w ] = rule14 ( n );
  elseif ( degree == 15 )
    [ x, w ] = rule15 ( n );
  elseif ( degree == 16 )
    [ x, w ] = rule16 ( n );
  elseif ( degree == 17 )
    [ x, w ] = rule17 ( n );
  elseif ( degree == 18 )
    [ x, w ] = rule18 ( n );
  elseif ( degree == 19 )
    [ x, w ] = rule19 ( n );
  elseif ( degree == 20 )
    [ x, w ] = rule20 ( n );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SQUARESYMVQ - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of DEGREE.\n' );
    error ( 'SQUARE_SYMQ - Fatal error!' );
  end

  w_sum = sum ( w(1:n) );

  w(1:n) = 4.0 * w(1:n) / w_sum;

  return
end


text:
classdef RWMOP25 < PROBLEM
% <multi> <real> <constrained>
% Process synthesis problem

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 2;
            obj.lower    = [0,-0.49];
            obj.upper    = [1.6,1.49];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1); x2 = round(x(:,2));
            % Objective function
            f(:,1) = x2 + 2*x1;
            f(:,2) = -x1.^2 -x2;
            % Constraints
            g(:,1) = f(:,2) +1.25;
            g(:,2) = x1 + x2 -1.6;
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [3.2000000e+00  -1.2500000e+00];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function title = p00_title ( problem )

%*****************************************************************************80
%
%% P00_TITLE returns a title for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the number of the problem.
%
%    Output, string TITLE, a title for the problem.
%
  if ( problem == 1 )
    title = p01_title ( );
  elseif ( problem == 2 )
    title = p02_title ( );
  elseif ( problem == 3 )
    title = p03_title ( );
  elseif ( problem == 4 )
    title = p04_title ( );
  elseif ( problem == 5 )
    title = p05_title ( );
  elseif ( problem == 6 )
    title = p06_title ( );
  elseif ( problem == 7 )
    title = p07_title ( );
  elseif ( problem == 8 )
    title = p08_title ( );
  elseif ( problem == 9 )
    title = p09_title ( );
  elseif ( problem == 10 )
    title = p10_title ( );
  elseif ( problem == 11 )
    title = p11_title ( );
  elseif ( problem == 12 )
    title = p12_title ( );
  elseif ( problem == 13 )
    title = p13_title ( );
  elseif ( problem == 14 )
    title = p14_title ( );
  elseif ( problem == 15 )
    title = p15_title ( );
  elseif ( problem == 16 )
    title = p16_title ( );
  elseif ( problem == 17 )
    title = p17_title ( );
  elseif ( problem == 18 )
    title = p18_title ( );
  elseif ( problem == 19 )
    title = p19_title ( );
  elseif ( problem == 20 )
    title = p20_title ( );
  elseif ( problem == 21 )
    title = p21_title ( );
  elseif ( problem == 22 )
    title = p22_title ( );
  elseif ( problem == 23 )
    title = p23_title ( );
  elseif ( problem == 24 )
    title = p24_title ( );
  elseif ( problem == 25 )
    title = p25_title ( );
  elseif ( problem == 26 )
    title = p26_title ( );
  elseif ( problem == 27 )
    title = p27_title ( );
  elseif ( problem == 28 )
    title = p28_title ( );
  elseif ( problem == 29 )
    title = p29_title ( );
  elseif ( problem == 30 )
    title = p30_title ( );
  elseif ( problem == 31 )
    title = p31_title ( );
  elseif ( problem == 32 )
    title = p32_title ( );
  elseif ( problem == 33 )
    title = p33_title ( );
  elseif ( problem == 34 )
    title = p34_title ( );
  elseif ( problem == 35 )
    title = p35_title ( );
  elseif ( problem == 36 )
    title = p36_title ( );
  elseif ( problem == 37 )
    title = p37_title ( );
  elseif ( problem == 38 )
    title = p38_title ( );
  elseif ( problem == 39 )
    title = p39_title ( );
  elseif ( problem == 40 )
    title = p40_title ( );
  elseif ( problem == 41 )
    title = p41_title ( );
  elseif ( problem == 42 )
    title = p42_title ( );
  elseif ( problem == 43 )
    title = p43_title ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_TITLE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of PROBLEM = %d', problem );
    error ( 'P00_TITLE - Fatal error!' );
  end

  return
end


text:
function point = cart2geod(src, curve)
%CART2GEOD Convert cartesian coordinates to geodesic coord.
%
%   PT2 = cart2geod(PT1, CURVE)
%   PT1 is the point to transform, in Cartesian coordinates (same system
%   used for the curve).
%   CURVE is a N-by-2 array which represents coordinates of curve vertices.
%
%   The function first compute the projection of PT1 on the curve. Then,
%   the first geodesic coordinate is the length of the curve to the
%   projected point, and the second geodesic coordinate is the
%   distance between PT1 and it projection.
%
%
%   TODO : add processing of points not projected on the curve.
%   -> use the closest end
%
%   See also
%   polylines2d, geod2cart, curveLength
%

% ---------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% created the 08/04/2004.
% Copyright 2010 INRA - Cepia Software Platform.

%   HISTORY
%   15/02/2007 replace minDistance by minDistancePoints

% parametrization approximation
t = parametrize(curve);

% compute distance between each src point and the curve
[dist, ind] = minDistancePoints(src, curve);

% convert to 'geodesic' coordinate
point = [t(ind) dist];

% Old version:
% for i=1:size(pt1, 1)
%     [dist, ind] = minDistance(src(i,:), curve);
%     point(i,:) = [t(ind) dist];
% end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function dateStringUtc = local_time_to_utc( dateVectorLocal, format )
%
% local_time_to_utc -- convert local time/date information to UTC.
%
% dateStringUtc = local_time_to_utc( dateVectorLocal ) converts a vector of local
%    date/time information to UTC dates and times.  The format of dateVectorLocal can be
%    either Matlab datenumbers (N x 1 or 1 x N) or else a vector of Matlab date vectors
%    (ie, N x 6, for example [2008 10 22 13 41 00 ; 2008 10 22 13 41 12]).  The time/date
%    values will be converted to UTC and returned in 'dd-mmm-yyyy HH:MM:SS Z' format.
%
% dateStringUtc = local_time_to_utc( dateVectorLocal, format ) returns the UTC time / date
%    values in a format specified by numeric argument format.  The format values match
%    those used by the Matlab datestr function.  Supported formats are as follows:
%
%    Number           String                     Example
%    ===========================================================================
%       0             'dd-mmm-yyyy HH:MM:SS Z'   01-Mar-2000 15:45:17 Z
%      13             'HH:MM:SS Z'               15:45:17 Z
%      15             'HH:MM Z'                  15:45 Z
%      21             'mmm.dd,yyyy HH:MM:SS Z'   Mar.01,2000 15:45:17 Z
%      30 (ISO 8601)  'yyyymmddTHHMMSSZ'         20000301T154517Z
%      31             'yyyy-mm-dd HH:MM:SS Z'    2000-03-01 15:45:17 Z
%
% Version date:  2008-November-10.
%

% Modification History:
%
%    2008-November-12, PT:
%        switch to use of Java methods and eliminate check for unix vs PC (java methods
%        work on all platforms).
%    2008-November-10, PT:
%        added a bunch of disp statements which should help diagnose mysterious smoke test
%        failures in PMD and PA.
%
%=========================================================================================

% arguments -- if the argument is a vector of datenums, convert to a vector of time/date
% vectors

  if ( size(dateVectorLocal,2) ~= 6 || max(dateVectorLocal(:,2)) > 12 )
      dateVectorLocal = datevec(dateVectorLocal) ;
  end

% check that the format number is supported

  if ( nargin == 1)
      format = 0 ;
  end
  supportedFormats = [0 13 15 21 30 31] ;
  if ~ismember(format,supportedFormats)
      error('common:localTimeToUtc:unsupportedFormat' , ...
          'local_time_to_utc:  unsupported format requested' ) ;
  end

% call the converter vector

  dateVectorUtc = local_date_vector_to_utc( dateVectorLocal ) ;

% convert the date vectors to the desired format

  dateStringUtc = datestr(dateVectorUtc,format) ;
  stringLength = size(dateStringUtc,2) ;
  nStrings = size(dateStringUtc,1) ;

% append the 'Z'

  dateStringUtc = append_z_to_date_string( dateStringUtc, format ) ;

return

% and that's it!

%
%
%

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

% subfunction which converts a Matlab date vector in local time to one in UTC

function dateVectorUtc = local_date_vector_to_utc( dateVectorLocal )

  nDates = size(dateVectorLocal,1) ;
  dateVectorUtc = zeros(nDates,6) ;

% import the Java classes needed for this process

  import java.text.SimpleDateFormat ;
  import java.util.Date ;
  import java.util.TimeZone ;

% instantiate a SimpleDateFormat object with a fixed time/date format and UTC time zone

  utcFormatObject = SimpleDateFormat('yyyy-MM-dd HH:mm:ss') ;
  utcFormatObject.setTimeZone(TimeZone.getTimeZone('UTC')) ;

% loop over date strings

  for iDate = 1:nDates

      dateVec = dateVectorLocal(iDate,:) ;

%     instantiate a Java Date class object with the local time.  Note that Java year is
%     year since 1900, and Java month is zero-based

      localDateObject = Date(dateVec(1)-1900, dateVec(2)-1, dateVec(3), ...
                             dateVec(4), dateVec(5), dateVec(6)) ;

%     convert the date object to a string in the correct format and in UTC

      dateStringUtc = char(utcFormatObject.format(localDateObject)) ;

%     pick through the resulting string and extract the data we want, converting to
%     numbers as we go

      dateVectorUtc(iDate,1) = str2num(dateStringUtc(1:4)) ;
      dateVectorUtc(iDate,2) = str2num(dateStringUtc(6:7)) ;
      dateVectorUtc(iDate,3) = str2num(dateStringUtc(9:10)) ;
      dateVectorUtc(iDate,4) = str2num(dateStringUtc(12:13)) ;
      dateVectorUtc(iDate,5) = str2num(dateStringUtc(15:16)) ;
      dateVectorUtc(iDate,6) = str2num(dateStringUtc(18:19)) ;

  end % loop over dates

return

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

% append a 'Z', signifying UTC, to a Matlab date string.

function datestrZ = append_z_to_date_string( datestr, format )

  nStrings = size(datestr,1) ;
  stringLength = size(datestr,2) ;
  if (format ~= 30)
      addString = repmat(' Z',nStrings,1) ;
      stringSizeIncrease = 2 ;
  else
      addString = repmat('Z',nStrings,1) ;
      stringSizeIncrease = 1 ;
  end
  datestrZ = [datestr(:) ; addString(:)] ;
  datestrZ = reshape(datestrZ,nStrings,stringLength+stringSizeIncrease) ;

return


text:
function test_suite = testAssertExceptionThrown
%testAssertExceptionThrown Unit tests for assertExceptionThrown

%   Steven L. Eddins
%   Copyright 2008 The MathWorks, Inc.

initTestSuite;

function test_happyCase
assertExceptionThrown(...
    @() error('MyProd:MyFun:MyId', 'my message'), 'MyProd:MyFun:MyId');

function test_wrongException
assertExceptionThrown(@() assertExceptionThrown(...
    @() error('MyProd:MyFun:MyId', 'my message'), ...
    'MyProd:MyFun:DifferentId'), 'assertExceptionThrown:wrongException');

function test_noException
assertExceptionThrown(@() assertExceptionThrown(@() sin(pi), 'foobar'), ...
    'assertExceptionThrown:noException');

author: ricedsp


text:
function [region_obj, results_table, varargout] = table(w, varargin)
% Print a table of all regions in a thresholded image. Return a Matlab table object
% and labeled region object. Used with fmri_data, statistic_image, and
% image_vector objects.
%
% This function has two "modes":
% - Without an optional atlas object input, it uses region.table to generate a table based
% on contiguous clusters (separated by those with positive and negative peak values)
%
% - With an optional atlas object input, it uses image_vector.subdivide_by_atlas to
% generate a table based on atlas regions covered by the input image. The table has
% one row per atlas region. The region object is also divided by the atlas,
% so one large contiguous blob in the image may cover multiple rows with
% different labeled regions.
%
% :Usage:
% ::
%
%    [poscl, negcl] = table(cl, [optional inputs])
%
% - By default, region.table() separates clusters into subregions with positive
%   and negative values. Thus, the number of rows may not match the original number of regions.
%   To turn this feature off, use 'nosep'.
% - By default, region.table() re-sorts the regions to group the table rows by macro-scale brain
%   structures (cortex, basal ganglia, etc.). So the ordering in the table may not match
%   the original region object. To turn this feature off, use 'nosort'.
%
% :Optional inputs:
%
%   **atlas_obj**
%       Any atlas-class object, e.g., loaded by load_atlas().
%
%   **k:**
%        Print only regions with k or more contiguous voxels
%
%   **nosep:**
%        do not separate cl with pos and neg effects based on peak in .val
%
%   **names:**
%        name clusters manually before printing to table and output; saves in .shorttitle field
%
%   **forcenames:**
%        force manual naming of cl by removing existing names in .shorttitle field
%
%   **nosort:**
%        Do not sort rows by network/brain lobe [default is to sort]
%
%   **legacy:**
%        force manual naming of cl by removing existing names in .shorttitle field
%
%   **nolegend:**
%        omit table legend
%
% :Outputs:
%
%   Returns region objects for cl with pos and neg effects
%   - autolabeled if Neuroimaging_Pattern_Masks and atlas tools are available on Matlab path
%   - limited by size if entered
%   - manually named if requested (see optional inputs)
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2011  Tor Wager
%
%     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/>.
% ..
% :Examples:
% -------------------------------------------------------------------------
% ::
% Example 1:
% % Complete group analysis of a standard dataset
% % Do analysis and prep results region object:
%
%   img_obj = load_image_set('emotionreg');         % Load a dataset
%   t = ttest(img_obj, .005, 'unc');                % Do a group t-test
%   t = threshold(t, .005, 'unc', 'k', 10);         % Re-threshold with extent threshold of 10 contiguous voxels
%   r = region(t);                                  % Turn t-map into a region object with one element per contig region
%
%   Label regions and print a table:
%   [r, region_table, table_legend_text] = autolabel_regions_using_atlas(r);
%                                                   % Label regions. Can be skipped because 'table' below attempts to do this automatically
%   table(r);                                       % Print a table of results using new region names
%
%   [rpos, rneg] = table(r);                        % Print and table and return region object separated into regions with positive vs. negative statistic values (from .Z field)

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..
%    July 2018:  Autolabel update and "new 2018 version", Tor Wager. Also added legend text.

n_cols = 140;                       % 140 good for HTML reports
sep_str = repmat('_', 1, n_cols);   % see textwrap

k = 0;
dosep = true;
donames = false;        % name clusters before printing to table and output; saves in .shorttitle field (legacy only)
forcenames = false;     % force naming of cl by removing existing names in .shorttitle field (legacy only)
dolegacy = false;
dosortrows = true;          % sort rows by area
dolegend = true;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            % functional commands
            case {'k', 'maxsize'}, k = varargin{i+1};
            case 'nosep', dosep = 0;
            case {'names', 'name', 'donames'}, donames = true;
            case 'forcenames', forcenames = true;

            case {'nosort', 'nosortrows'}, dosortrows = false;

            case {'legacy', 'dolegacy'}, dolegacy = true;

            case 'nolegend', dolegend = false;

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

[poscl, negcl, results_table] = table(cl, varargin)

[w_atlas, w_region_obj] = subdivide_by_atlas(w, atl);

% Make a table of labeled regions with significant voxels in the weight map

[region_obj, region_table] = autolabel_regions_using_atlas(w_region_obj, atl);
region_table = region_table(:, 1:3);


text:
% dftfilt3() - discrete complex wavelet filters
%
% Usage:
%   >> [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin)
%
% Inputs:
%   freqs    - vector of frequencies of interest.
%   cycles   - cycles array. If cycles=0, then the Hanning tapered Short-term FFT is used.
%              If one value is given and cycles>0, all wavelets have
%              the same number of cycles. If two values are given, the
%              two values are used for the number of cycles at the lowest
%              frequency and at the highest frequency, with linear or
%              log-linear interpolation between these values for intermediate
%              frequencies
%   srate    - sampling rate (in Hz)
%
% Optional Inputs: Input these as 'key/value pairs.
%   'cycleinc' - ['linear'|'log'] increase mode if [min max] cycles is
%              provided in 'cycle' parameter. {default: 'linear'}
%   'winsize'  Use this option for Hanning tapered FFT or if you prefer to set the length of the
%              wavelets to be equal for all of them (e.g., to set the
%              length to 256 samples input: 'winsize',256). {default: [])
%              Note: the output 'wavelet' will be a matrix and it may be
%              incompatible with current versions of timefreq and newtimef.
%   'timesupport' The number of temporal standard deviation used for wavelet lengths {default: 7)
%
% Output:
%   wavelet - cell array or matrix of wavelet filters
%   timeresol - temporal resolution of Morlet wavelets.
%   freqresol - frequency resolution of Morlet wavelets.
%
% Note: The length of the window is always made odd.
%
% Authors: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 3/28/2003
%          Rey Ramirez, SCCN/INC/UCSD, La Jolla, 9/26/2006

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

%
% Revision 1.12 2006/09/25  rey r
% Almost complete rewriting of dftfilt2.m, changing both Morlet and Hanning
% DFT to be more in line with conventional implementations.
%
% Revision 1.11  2006/09/07 19:05:34  scott
% further clarified the Morlet/Hanning distinction -sm
%
% Revision 1.10  2006/09/07 18:55:15  scott
% clarified window types in help msg -sm
%
% Revision 1.9  2006/05/05 16:17:36  arno
% implementing cycle array
%
% Revision 1.8  2004/03/04 19:31:03  arno
% email
%
% Revision 1.7  2004/02/25 01:45:55  arno
% sinus test
%
% Revision 1.6  2004/02/15 22:23:08  arno
% implementing morlet wavelet
%
% Revision 1.5  2003/05/09 20:55:10  arno
% adding hanning function
%
% Revision 1.4  2003/04/29 16:02:54  arno
% header typos
%
% Revision 1.3  2003/04/29 01:09:16  arno
% debug imaginary part
%
% Revision 1.2  2003/04/28 23:01:13  arno
% *** empty log message ***
%
% Revision 1.1  2003/04/28 22:46:49  arno
% Initial revision
%

function [wavelet,cycles,freqresol,timeresol] = dftfilt3( freqs, cycles, srate, varargin);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rey fixed all input parameter sorting.
if nargin < 3
    error(' A minimum of 3 arguments is required');
end;
numargin=length(varargin);
if rem(numargin,2)
    error('There is an uneven number key/value inputs. You are probably missing a keyword or its value.')
end
varargin(1:2:end)=lower(varargin(1:2:end));

% Setting default parameter values.
cycleinc='linear';
winsize=[];
timesupport=7;  % Setting default of 7 temporal standard deviations for wavelet's length.

for n=1:2:numargin
    keyword=varargin{n};
    if strcmpi('cycleinc',keyword)
        cycleinc=varargin{n+1};
    elseif strcmpi('winsize',keyword)
        winsize=varargin{n+1};
        if ~mod(winsize,2)
            winsize=winsize+1; % Always set to odd length wavelets and hanning windows;
        end
    elseif strcmpi('timesupport',keyword)
        timesupport=varargin{n+1};
    else
        error(['What is ' keyword '? The only legal keywords are: type, cycleinc, winsize, or timesupport.'])
    end
end
if isempty(winsize) & cycles==0
    error('If you are using a Hanning tapered FFT, please supply the winsize input-pair.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute number of cycles at each frequency
% ------------------------------------------
type='morlet';
if length(cycles) == 1 & cycles(1)~=0
    cycles = cycles*ones(size(freqs));
elseif length(cycles) == 2
    if strcmpi(cycleinc, 'log') % cycleinc
         cycles = linspace(log(cycles(1)), log(cycles(2)), length(freqs));
         cycles = exp(cycles);
         %cycles=logspace(log10(cycles(1)),log10(cycles(2)),length(freqs)); %rey
    else
        cycles = linspace(cycles(1), cycles(2), length(freqs));
    end;
end;
if cycles==0
    type='sinus';
end

sp=1/srate; % Rey added this line (i.e., sampling period).
% compute wavelet
for index = 1:length(freqs)
    fk=freqs(index);
    if strcmpi(type, 'morlet') % Morlet.
        sigf=fk/cycles(index); % Computing time and frequency standard deviations, resolutions, and normalization constant.
        sigt=1./(2*pi*sigf);
        A=1./sqrt(sigt*sqrt(pi));
        timeresol(index)=2*sigt;
        freqresol(index)=2*sigf;
        if isempty(winsize) % bases will be a cell array.
            tneg=[-sp:-sp:-sigt*timesupport/2];
            tpos=[0:sp:sigt*timesupport/2];
            t=[fliplr(tneg) tpos];
            psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t));
            wavelet{index}=psi;  % These are the wavelets with variable number of samples based on temporal standard deviations (sigt).
        else % bases will be a matrix.
            tneg=[-sp:-sp:-sp*winsize/2];
            tpos=[0:sp:sp*winsize/2];
            t=[fliplr(tneg) tpos];
            psi=A.*(exp(-(t.^2)./(2*(sigt^2))).*exp(2*i*pi*fk*t));
            wavelet(index,:)=psi; % These are the wavelets with the same length.
            % This is useful for doing time-frequency analysis as a matrix vector or matrix matrix multiplication.
        end
    elseif strcmpi(type, 'sinus') % Hanning
        tneg=[-sp:-sp:-sp*winsize/2];
        tpos=[0:sp:sp*winsize/2];
        t=[fliplr(tneg) tpos];
        win = exp(2*i*pi*fk*t);
        wavelet(index,:) = win .* hanning(winsize)';
        %wavelet{index} = win .* hanning(winsize)';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end;
end;

% symmetric hanning function
function w = hanning(n)
if ~rem(n,2)
    w = .5*(1 - cos(2*pi*(1:n/2)'/(n+1)));
    w = [w; w(end:-1:1)];
else
    w = .5*(1 - cos(2*pi*(1:(n+1)/2)'/(n+1)));
    w = [w; w(end-1:-1:1)];
end

author: PatternRecognition

repo: OpenBMI


text:
function [L, W, A, BE] = calc_tps(refx, refy, holx, holy, calc_BE) %coord in column vector style
%calculate 2D tps warp
n = length(refx);
P = [ ones(n,1) refx refy];
V = [holx'; holy'];
K = zeros(n,n);
BE = 0 ;
for i=1:n,
    for j=1:n,
        r_2= (refx(i) - refx(j))^2 + (refy(i) -refy(j))^2 ;
        if r_2 > 0, K(i,j) = r_2 * log(r_2); end;
    end;
end;
for i=1:n, K(i,i)=0; end;
L = [K P;P' zeros(3,3)];
W=zeros(2,n);
A=zeros(2,3);
Y=[V zeros(2,3)]';
temp = ( inv(L)*Y )';
W = temp(1:2,1:n);
A = temp(1:2,n+1:n+3);
if calc_BE == 0, return; end;

temp = zeros(n,n); % server as Ln-1
temp1 = inv(L);
Ln_inv = temp1(1:n,1:n);
[bas_u, s, bas_v] = svd(Ln_inv);%col of bas_v basis vectors
B = zeros(2,n-3); %storage of coefs of decomp of V
for i=1:n-3,
    B(1,i) = V(1,:)*bas_v(:,i);
    B(2,i) = V(2,:)*bas_v(:,i);
end;
BE = 0;
for i=1:n-3,
    BE = BE + s(i,i)*( B(1,i)^2 + B(2,i)^2 );
end;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function circles2d(varargin)
%CIRCLES2D Description of functions operating on circles
%
%   Circles are represented by their center and their radius:
%   C = [xc yc r];
%   One sometimes considers orientation of circle, by adding an extra
%   boolean value in 4-th position, with value TRUE for direct (i.e.
%   turning Counter-clockwise) circles.
%
%   Circle arcs are represented by their center, their radius, the starting
%   angle and the angle extent, both in degrees:
%   CA = [xc yc r theta0 dtheta];
%
%   Ellipses are represented by their center, their 2 semi-axis length, and
%   their angle (in degrees) with Ox direction.
%   E = [xc yc A B theta];
%
%   See also:
%   ellipses2d, createCircle, createDirectedCircle,
%   isPointInCircle, isPointOnCircle, enclosingCircle, circumCircle
%   intersectLineCircle, intersectCircles, radicalAxis
%   circleToPolygon, circleArcToPolyline
%   drawCircle, drawCircleArc
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2010 INRA - Cepia Software Platform.

help('circles2d');


text:
function [ value, ifault ] = gammds ( x, p )

%*****************************************************************************80
%
%% GAMMDS computes the incomplete Gamma integral.
%
%  Discussion:
%
%    The parameters must be positive.  An infinite series is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Chi Leung Lau.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Chi Leung Lau,
%    Algorithm AS 147:
%    A Simple Series for the Incomplete Gamma Integral,
%    Applied Statistics,
%    Volume 29, Number 1, 1980, pages 113-114.
%
%  Parameters:
%
%    Input, real X, P, the arguments of the incomplete
%    Gamma integral.  X and P must be greater than 0.
%
%    Output, real VALUE, the value of the incomplete
%    Gamma integral.
%
%    Output, integer IFAULT, error flag.
%    0, no errors.
%    1, X <= 0 or P <= 0.
%    2, underflow during the computation.
%
  e = 1.0E-09;
  uflo = 1.0E-37;
%
%  Check the input.
%
  if ( x <= 0.0 )
    ifault = 1;
    value = 0.0;
    return
  end

  if ( p <= 0.0 )
    ifault = 1;
    value = 0.0;
    return
  end
%
%  ALNGAM is the natural logarithm of the gamma function.
%
  arg = p * log ( x ) - alngam ( p + 1.0 ) - x;

  if ( arg < log ( uflo ) )
    value = 0.0;
    ifault = 2;
    return
  end

  f = exp ( arg );

  if ( f == 0.0 )
    value = 0.0;
    ifault = 2;
    return
  end

  ifault = 0;
%
%  Series begins.
%
  c = 1.0;
  value = 1.0;
  a = p;

  while ( 1 )

    a = a + 1.0;
    c = c * x / a;
    value = value + c;

    if ( c <= e * value )
      break;
    end

  end

  value = value * f;

  return
end


text:
% Test file for chebtech1 constructor.
% Here, we check populate().  (This function is not user-facing.)

function pass = test_constructor(pref)

% Get preferences:
if ( nargin < 1 )
    pref = chebtech.techPref();
end

% Initialize with default data:
data = chebtech.parseDataInputs(struct());

%%
% Test on a scalar-valued function:
pref.refinementFunction = 'nested';
f = @(x) sin(x);
g = populate(chebtech1, f, data, pref);
x = chebtech1.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(1) = norm(f(x) - values, inf) < 10*vscale(g)*eps;

% Test on an array-valued function:
pref.refinementFunction = 'nested';
f = @(x) [sin(x) cos(x) exp(x)];
g = populate(chebtech1, f, data, pref);
x = chebtech1.chebpts(length(g.coeffs));
values = g.coeffs2vals(g.coeffs);
pass(2) = norm(f(x) - values, inf) < 10*max(vscale(g)*eps);

%%
% Test on a scalar-valued function:
pref.refinementFunction = 'resampling';
f = @(x) sin(x);
[g, values] = populate(chebtech1, f, data, pref);
x = chebtech1.chebpts(length(values));
pass(3) = norm(f(x) - values, inf) < 10*vscale(g)*eps;

% Test on an array-valued function:
pref.refinementFunction = 'resampling';
f = @(x) [sin(x) cos(x) exp(x)];
[g, values] = populate(chebtech1, f, data, pref);
x = chebtech1.chebpts(length(values));
pass(4) = norm(f(x) - values, inf) < 10*max(vscale(g)*eps);

%%
% Some other tests:

% This should fail with an error:
try
    f = @(x) x + NaN;
    populate(chebtech1, f, data, pref);
    pass(5) = false;
catch ME
    pass(5) = strcmp(ME.message, 'Too many NaNs/Infs to handle.');
end

% As should this:
try
    f = @(x) x + Inf;
    populate(chebtech1, f, data, pref);
    pass(6) = false;
catch ME
    pass(6) = strcmp(ME.message, 'Too many NaNs/Infs to handle.');
end

% Check that things don't crash if pref.minSamples and pref.maxLength are equal.
try
    pref.minSamples = 8;
    pref.maxLength = 8;
    populate(chebtech1, @sin, data, pref);
    pass(7) = true;
catch
    pass(7) = false;
end

% Test logical-valued functions:
f = chebtech1(@(x) x > -2);
g = chebtech1(1);
pass(8) = normest(f - g) < eps;

f = chebtech1(@(x) x < -2);
g = chebtech1(0);
pass(9) = normest(f - g) < eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [headers, fid] = readPLXHeaders(filename, fullread)
% readPLXHeaders - Read file headers and channel information from PLX file.
%
% [headers, fid] = readPLXHeaders(filename)
%
% readPLXHeaders will read the header data from a PLX file and return a
% structure with fields corresponding to each piece of header data.
%
% Input:
%   filename - Filename of a PLX file. A valid file-ID can be provided
%              instead of a filename. In this case, the file-ID will not be
%              closed after reading.
%   fullread - Verify the contents of the file header by reading the entire
%              file. This effects the values of headers.tscounts,
%              headers.wfcounts, headers.evcounts, and headers.slowcounts.
%              It will also verify that headers.chans(n).nunits is
%              accurate, and that headers.numPointsWave is accurate.
%              (default = false)
%
% Output:
%   headers - PLX file header information.
%   fid     - File-ID of the
%
% WARNING: Plexon didn't have much foresight for the number of
% units/channel, and only left enough space in the headers for 4
% units/channel. Therefore, the data stored in the headers (and output as
% headers.tscounts and headers.wfcounts) can't account for more than 4
% units/channel. The only way to get an accurate count of the number of
% timestamps, waveforms, and events is to read the entire data file. If you
% call this function with the 'fullread' argument set false (default) it
% will only read the values provided in the header. If you set 'fullread'
% to true, it will read the entire file to verify the information in the
% header. This is similar to the behavior of the Plexon function called
% plx_info.
%
% Note: This function behaves slightly different than plx_info. In
% particular:
%     tscounts   = 27 x m array. The entry at tscounts(c,u) refers to
%                  channel c and unit u-1. Channel numbers start at 1.
%                  Unit 0 (in row 1) is unsorted waves.
%     wfcounts   = 27 x m array. The entry at wfcounts(c,u) refers to
%                  channel c and unit u-1. Channel numbers start at 1.
%                  Unit 0 (in row 1) is unsorted waves.
%     evcounts   = 1 x m array. The entry at evcounts(c) refers to
%                  channel c. Channel numbers start at 1. m is the event
%                  with the largest channel number (usually 259).
%     slowcounts = 1 x m+1 array. The entry at slowcounts(c) refers to
%                  channel c+1. Channel numbers start at 0. m is one
%                  greater than the largest slow channel number.
%
% See also writePLXFile, checkPLXData, plx_info, plx_information.
%
% AUTHOR: Benjamin Kraus (bkraus@bu.edu, ben@benkraus.com)
% Copyright (c) 2011, Benjamin Kraus
% $Id$

DEBUG = true;

if (nargin == 0 || isempty(filename))
    FilterSpec = {'*.plx', 'Plexon PLX File (*.plx)';
                  '*', 'All Files'};
    [fname, pathname] = uigetfile(FilterSpec, 'Select a Plexon PLX file');
    if(fname == 0); headers = struct([]); return; end
    filename = strcat(pathname, fname);
end

if(~ischar(filename));
    % If the filename input is not a character string, check to see if it
    % is a valid file ID instead. We can do this by calling 'frewind'. If
    % 'frewind' succeeds, then we know we have a valid file ID.
    fid = filename;
    frewind(fid);
else
    fid = fopen(filename, 'r');
    if(fid ~= -1 && nargout < 2); c = onCleanup(@()fclose(fid)); end
end

if(fid == -1); error('readPLXHeaders:FileError','Error opening file'); end

if(nargin < 2); fullread = false; shortread = false;
elseif(islogical(fullread)); shortread = false;
elseif(ischar(fullread))
    switch lower(fullread)
        case 'true';   fullread = true;  shortread = false;
        case 'false';  fullread = false; shortread = false;
        case 'full';   fullread = true;  shortread = false;
        case 'short';  fullread = false; shortread = true;
        otherwise;     fullread = false; shortread = false;
    end
end

% Determine file size:
fseek(fid,0,'eof');
fsize = ftell(fid);
frewind(fid);

% File header is 256 bytes

magic = fread(fid, 4, '*char')';
if(~strcmp(magic,'PLEX'))
    error('readPLXHeaders:InvalidPLXFile','This file is not a valid .plx file');
end

headers.version = fread(fid, 1, 'int32');
headers.comment = deblank(fread(fid, 128, '*char')');
headers.ADFrequency = fread(fid, 1, 'int32');
headers.numDSPChannels = fread(fid, 1, 'int32');
headers.numEventChannels = fread(fid, 1, 'int32');
headers.numSlowChannels = fread(fid, 1, 'int32');
headers.numPointsWave = fread(fid, 1, 'int32');
headers.numPointsPreThr = fread(fid, 1, 'int32');

YR = fread(fid, 1, 'int32');
MO = fread(fid, 1, 'int32');
DA = fread(fid, 1, 'int32');
HR = fread(fid, 1, 'int32');
MI = fread(fid, 1, 'int32');
SC = fread(fid, 1, 'int32');
headers.date = datenum([YR, MO, DA, HR, MI, SC]);

headers.fastread = fread(fid, 1, 'int32');
headers.waveformfreq = fread(fid, 1, 'int32');
headers.lasttimestamp = fread(fid, 1, 'double');

if(headers.version >= 103)
    headers.trodalness = fread(fid, 1, 'uint8');
    headers.datatrodalness = fread(fid, 1, 'uint8');
    headers.bitsperspikesample = fread(fid, 1, 'uint8');
    headers.bitsperslowsample = fread(fid, 1, 'uint8');
    headers.spikeMaxMagnitudeMV = fread(fid, 1, 'ushort');
    headers.slowMaxMagnitudeMV = fread(fid, 1, 'ushort');
else
    headers.trodalness = 1;
    headers.datatrodalness = 1;
    headers.bitsperspikesample = 12;
    headers.bitsperslowsample = 12;
    headers.spikeMaxMagnitudeMV = 3000;
    headers.slowMaxMagnitudeMV = 5000;
    fread(fid, 8, 'uint8');
end

if(headers.version >= 105)
    headers.spikePreAmpGain = fread(fid, 1, 'ushort');
else
    headers.spikePreAmyGain = 1000;
    fread(fid, 1, 'ushort');
end

if(headers.version >= 106)
    headers.acquiringsoftware = deblank(fread(fid, 18, '*char')');
    headers.processingsoftware = deblank(fread(fid, 18, '*char')');
    fread(fid, 10, 'uint8');
else
    headers.acquiringsoftware = '';
    headers.processingsoftware = '';
    fread(fid, 46, 'char');
end

% Current position = 256

% Timestamps, waveforms, event counts total 7248 bytes.

tscounts = fread(fid, [5, 130], 'int32');
wfcounts = fread(fid, [5, 130], 'int32');
evcounts = fread(fid, [1, 512], 'int32');

headers.tscounts = tscounts(:,2:min(headers.numDSPChannels+1,130));
headers.wfcounts = wfcounts(:,2:min(headers.numDSPChannels+1,130));
headers.evcounts = evcounts(2:300);
headers.slowcounts = evcounts(301:min(300+headers.numSlowChannels,512));

% Current position = 7504

% If instructed to do a 'short' read, stop here.
if(shortread); return; end

% Each channel header is 1020 bytes.
headers.chans = struct([]);
for ii = 1:headers.numDSPChannels;
    headers.chans(ii,1).name = deblank(fread(fid, 32, '*char')');
    headers.chans(ii,1).signame = deblank(fread(fid, 32, '*char')');
    headers.chans(ii,1).channel = fread(fid, 1, 'int32');
    headers.chans(ii,1).wfrate = fread(fid, 1, 'int32');
    headers.chans(ii,1).sig = fread(fid, 1, 'int32');
    headers.chans(ii,1).ref = fread(fid, 1, 'int32');
    headers.chans(ii,1).gain = fread(fid, 1, 'int32');
    headers.chans(ii,1).filter = fread(fid, 1, 'int32');
    headers.chans(ii,1).threshold = fread(fid, 1, 'int32');
    headers.chans(ii,1).method = fread(fid, 1, 'int32');
    headers.chans(ii,1).nunits = fread(fid, 1, 'int32');
    headers.chans(ii,1).template = fread(fid, [5, 64], 'short');
    headers.chans(ii,1).fit = fread(fid, 5, 'int32');
    headers.chans(ii,1).sortwidth = fread(fid, 1, 'int32');
    headers.chans(ii,1).boxes = zeros(5, 2, 4);
    headers.chans(ii,1).boxes(:) = fread(fid, 5*2*4, 'short');
    headers.chans(ii,1).sortbeg = fread(fid, 1, 'int32');
    if(headers.version>=105)
        headers.chans(ii,1).comment = deblank(fread(fid, 128, '*char')');
    else
        headers.chans(ii,1).comment = '';
        fread(fid, 128, 'char');
    end

    if(headers.version>=106)
        headers.chans(ii,1).srcid = fread(fid, 1, 'uchar');
        fread(fid, 1, 'uchar');
        headers.chans(ii,1).chanid = fread(fid, 1, 'ushort');
        fread(fid, 10, 'int32');
    else fread(fid, 11, 'int32');
    end
end

% Current position = 7504 + 1020*nchans

% Each event channel header is 296 bytes
headers.evchans = struct([]);
for ii = 1:headers.numEventChannels;
    headers.evchans(ii,1).name = deblank(fread(fid, 32, '*char')');
    headers.evchans(ii,1).num = fread(fid, 1, 'int32');

    if(headers.version>=105)
        headers.evchans(ii,1).comment = deblank(fread(fid, 128, '*char')');
    else
        headers.evchans(ii,1).comment = '';
        fread(fid, 128, 'char');
    end

    if(headers.version>=106)
        headers.evchans(ii,1).srcid = fread(fid, 1, 'uchar');
        fread(fid, 1, 'uchar');
        headers.evchans(ii,1).chanid = fread(fid, 1, 'ushort');
        fread(fid, 32, 'int32');
    else fread(fid, 33, 'int32');
    end
end

% Current position = 7504 + 1020*nchans + 296*nevchans

% Each slow channel header is 296 bytes
headers.slowchans = struct([]);
for ii = 1:headers.numSlowChannels;
    headers.slowchans(ii,1).name = deblank(fread(fid, 32, '*char')');
    headers.slowchans(ii,1).channel = fread(fid, 1, 'int32');
    headers.slowchans(ii,1).ADFrequency = fread(fid, 1, 'int32');
    headers.slowchans(ii,1).adgain = fread(fid, 1, 'int32');
    headers.slowchans(ii,1).enabled = fread(fid, 1, 'int32');
    headers.slowchans(ii,1).preampgain = fread(fid, 1, 'int32');

    if(headers.version>=104)
        headers.slowchans(ii,1).spikechannel = fread(fid, 1, 'int32');
    else
        headers.slowchans(ii,1).spikechannel = 0;
        fread(fid, 1, 'int32');
    end

    if(headers.version>=105)
        headers.slowchans(ii,1).comment = deblank(fread(fid, 128, '*char')');
    else
        headers.slowchans(ii,1).comment = '';
        fread(fid, 128, 'char');
    end

    if(headers.version>=106)
        headers.slowchans(ii,1).srcid = fread(fid, 1, 'uchar');
        fread(fid, 1, 'uchar');
        headers.slowchans(ii,1).chanid = fread(fid, 1, 'ushort');
        fread(fid, 27, 'int32');
    else fread(fid, 28, 'int32');
    end
end

% Current position = 7504 + 1020*nchans + 296*nevchans + 296*nslowchans

% Keep track of the current position.
headers.datastart = ftell(fid);

% Unless instructed to read the entire file, stop here.
if(~fullread); return; end
if(DEBUG); fprintf(1,'Starting full file read.\n'); end

% If we are still going, then we were instructed to read the entire file to
% get accurate counts for spikes, events, and waveforms.

% Lets start by looking up the channel numbers.
if(headers.numDSPChannels>0)
    channums = vertcat(headers.chans.channel);
    maxchan = max(channums);
else
    maxchan = 0;
end

if(headers.numEventChannels>0)
    evnums = vertcat(headers.evchans.num);
    maxev = max(evnums);
else
    maxev = 0;
end

if(headers.numSlowChannels>0)
    slownums = vertcat(headers.slowchans.channel);
    maxslow = max(slownums);
else
    maxslow = 0;
end

tscounts   = zeros(27, maxchan);
wfcounts   = zeros(27, maxchan);
evcounts   = zeros(1,  maxev);
slowcounts = zeros(1,  maxslow+1);
slowfrags  = zeros(1,  maxslow+1);

% Read data headers and count number of each type of channel.
fseek(fid,headers.datastart,'bof');

numpointswave = headers.numPointsWave;

    function fillbuffer()
        dbufloc(1) = dbufloc-ndbuf;
        ndbuf(1) = min(floor((fsize-ftell(fid))/2),dbufsize);
        dbuf(1:ndbuf) = fread(fid,ndbuf,'*short');
    end

datablocks = 0;

dbufsize = 1024*1024*2;
dbuf = int16(zeros(dbufsize,1));
ndbuf = 0;
dbufloc = 1;

fillbuffer();

% type    = db(1);   db(end-7);
% upperts = db(2);   db(end-6);
% lowerts = db(3:4); db(end-5:end-4);
% chan    = db(5);   db(end-3);
% unit    = db(6);   db(end-2);
% nwaves  = db(7);   db(end-1);
% nwords  = db(8);   db(end);
dbh = int16(zeros(8,1));
needvals = 0;
havevals = 0;

while(ndbuf > 0)
    datablocks = datablocks + 1;
    if(DEBUG && mod(datablocks,100000)==0);
        fprintf('Block: %.0f\n',datablocks);
    end

    if(ndbuf-dbufloc>=7)
        dbh(1:8) = dbuf(dbufloc:dbufloc+7);
        dbufloc = dbufloc+8;
    else
        havevals(1) = max(0,ndbuf-dbufloc+1);
        dbh(1:havevals) = dbuf(dbufloc:ndbuf);
        needvals(1) = 8-havevals;
        dbufloc = dbufloc+havevals;
        fillbuffer();
        if(ndbuf-dbufloc+1<needvals); break; end
        dbh(havevals+1:8) = dbuf(dbufloc:dbufloc+needvals-1);
        dbufloc = dbufloc+needvals;
    end
    dbufloc = dbufloc+double(dbh(7)*dbh(8));

    if(dbh(1)==1)
        tscounts(dbh(6)+1,dbh(5)) = tscounts(dbh(6)+1,dbh(5)) + 1;
        wfcounts(dbh(6)+1,dbh(5)) = wfcounts(dbh(6)+1,dbh(5)) + dbh(7);
        if(dbh(8) ~= numpointswave(end))
            numpointswave(end+1) = dbh(8); %#ok<AGROW>
        end
    elseif(dbh(1) == 4)
        evcounts(dbh(5)) = evcounts(dbh(5)) + 1;
    elseif(dbh(1) == 5)
        slowfrags(dbh(5)+1) = slowfrags(dbh(5)+1) + 1;
        slowcounts(dbh(5)+1) = slowcounts(dbh(5)+1) + double(dbh(7)*dbh(8));
    end
end

if(~isscalar(numpointswave))
    headers.numPointsWave = max(numpointswave);

    warning('readPLXHeaders:NumPointsWaveMismatch',...
        'Number of points per wave is inconsistent or doesn''t match header information.');
end

headers.tscounts = tscounts;
headers.wfcounts = wfcounts;
headers.evcounts = evcounts;
headers.slowfrags = slowfrags;
headers.slowcounts = slowcounts;

for ii = 1:headers.numDSPChannels
    headers.chans(ii,1).nunits = find(tscounts(:,ii)>0,1,'last')-1;
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Load original BFM
load('01_MorphableModel.mat');

% Load modified FW expression model
load('3DDFA_Release/Matlab/Model_Expression.mat');

% Load FW to BFM mapping
load('util/map_tddfa_to_basel.mat');
% Fix zero-based indexing
map_tddfa_to_basel=map_tddfa_to_basel+1;

model.shapePC = zeros(3*53490,10);
for i=1:5
model.shapePC(:,i) = shapePC(:,i).*shapeEV(i);
end

%% Extrapolate expressions to mouth interior via Laplace-Beltrami
vertices = double(reshape(shapeMU,3,53490)');
[L,~]=LaplaceBeltrami(vertices',tl');
% Selection matrix for BFM vertices to FW
S = sparse(1:length(map_tddfa_to_basel),double(map_tddfa_to_basel),ones(length(map_tddfa_to_basel),1),length(map_tddfa_to_basel),53490);
for i=1:5
    expression = ([L; S]\[L*vertices; S*vertices+reshape(w_exp(:,i),3,53215)'])-vertices;
    expression = expression';
    model.shapePC(:,i+5) = expression(:);
end
%% Load UV coordinates and resample model
load('util/BFM_UV.mat');
[ newmodel ] = resampleModel( shapeMU,model.shapePC,[],tl,UV,112 );
newmodel.shapeMU=newmodel.shapeMU./1000;
newmodel.shapePC=newmodel.shapePC./1000;
save model.mat -struct newmodel
%% Generate random mesh
clear FV
FV.faces = newmodel.faces;
alpha = randn(10,1);
FV.vertices = reshape(newmodel.shapeMU+newmodel.shapePC*alpha,3,112^2)';
figure; patch(FV, 'FaceColor', [1 1 1], 'EdgeColor', 'none', 'FaceLighting', 'phong'); axis equal; light; axis tight

author: anilbas


text:
function res = bf_write_spmeeg(BF, S)
% Writes out beamformer results as M/EEG dataset
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_write_spmeeg.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    mode         = cfg_menu;
    mode.tag     = 'mode';
    mode.name    = 'Writing mode';
    mode.help    = {'How to generate the output'};
    mode.labels  = {
        'Write new'
        'Online montage original data'
        'Copy + online montage'
        }';
    mode.values  = {
        'write'
        'online'
        'onlinecopy'
        }';
    mode.val = {'write'};

    modality         = cfg_menu;
    modality.tag     = 'modality';
    modality.name    = 'Modality';
    modality.help    = {'What modality to output'};
    modality.labels  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.values  = {
        'MEG'
        'MEGPLANAR'
        'EEG'
        }';
    modality.val = {'MEG'};

    none = cfg_const;
    none.tag = 'none';
    none.name = 'None';
    none.val  = {0};

    addchannels      = cfg_choice;
    addchannels.tag  = 'addchannels';
    addchannels.name = 'Extra channels to add';
    addchannels.values  = {none, spm_cfg_eeg_channel_selector};
    addchannels.val  = {none};

    prefix         = cfg_entry;
    prefix.tag     = 'prefix';
    prefix.name    = 'Filename Prefix';
    prefix.help    = {'Specify the string to be prepended to the output (if relevant).'};
    prefix.strtype = 's';
    prefix.num     = [1 Inf];
    prefix.val     = {'B'};

    spmeeg      = cfg_branch;
    spmeeg.tag  = 'spmeeg';
    spmeeg.name = 'SPM M/EEG dataset';
    spmeeg.val  = {mode, modality, addchannels, prefix};

    res = spmeeg;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

D = BF.data.D;

if isfield(S.addchannels, 'channels')
    addchannels = D.chanlabels(D.selectchannels(spm_cfg_eeg_channel_selector(S.addchannels.channels)));
else
    addchannels = {};
end

usemontage = 0;

if isfield(BF.output, 'montage')

    usemontage = 1;

    montage = BF.output.montage.(S.modality);

    [sel1, sel2] = spm_match_str(montage.labelorg, addchannels);
    for i = 1:length(sel1)
        montage.labelnew(end+1) = addchannels(sel2(i));
        montage.tra(end+1, sel1(i)) = 1;

        montage.chantypenew(end+1) = chantype(D, D.indchannel(addchannels(sel2(i))))';
        montage.chanunitnew(end+1) = units(D, D.indchannel(addchannels(sel2(i))))';

    end
    addchannels(sel2) = [];

    if ~isempty(addchannels)
        montage.labelorg = [montage.labelorg(:); addchannels(:)];
        montage.labelnew = [montage.labelnew(:); addchannels(:)];

        montage.chantypenew = [montage.chantypenew(:); chantype(D, D.indchannel(addchannels))'];
        montage.chanunitnew = [montage.chanunitnew(:); units(D, D.indchannel(addchannels))'];

        na = numel(addchannels);
        montage.tra((end+1):(end+na), (end+1):(end+na)) = eye(na);
    end

    montage.chantypeold = chantype(D, D.indchannel(montage.labelorg))';
    montage.chanunitold = units(D, D.indchannel(montage.labelorg))';
elseif isfield(BF.output, 'sourcedata')
    ftdata = BF.output.sourcedata.(S.modality).ftdata;

    if ~isempty(addchannels)
        addind = D.indchannel(addchannels);
        ftdata.label = [ftdata.label;D.chanlabels(addind)'];
        if numel(ftdata.trial) == length(BF.features.trials)
            for i = 1:numel(ftdata.trial)
                ftdata.trial{i} = [ftdata.trial{i};D(addind, D.indsample(ftdata.time{i}(1)):D.indsample(ftdata.time{i}(end)), BF.features.trials(i))];
            end
        else
            error('Cannot match trials for adding extra channels.')
        end
    end

    Ds = spm_eeg_ft2spm(ftdata, [S.prefix D.fname]);
    Ds = trialonset(Ds, ':', D.trialonset);

    if isfield(BF.output.sourcedata.(S.modality), 'events')
        for i = 1:Ds.ntrials
            evold = events(D, BF.features.trials(i));
            if iscell(evold)
                evold = evold{1};
            end

            evnew = BF.output.sourcedata.(S.modality).events;
            if iscell(evnew)
                evnew = evnew{i};
            end

            ev = spm_cat_struct(evold, evnew);

            Ds = events(Ds, i, ev);
        end
    else
        if ~isempty(D.events)
            Ds = events(Ds, ':', D.events);
        end
    end

    if isfield(BF.sources, 'voi')
        Ds = chantype(Ds, ':', 'LFP');
    else
        Ds = chantype(Ds, ':', 'SRC');
    end

    if ~isempty(addchannels)
         Ds = chantype(Ds, (Ds.nchannels-length(addind)+1):Ds.nchannels, D.chantype(addind));
    end

    save(Ds);

    D = Ds;
else
    error('Unsupported option');
end

if usemontage
    if D.montage('getindex')
        vmontage = rmfield(D.montage('getmontage'), 'channels');
        montage  =  ft_apply_montage(vmontage, montage);
    end

    S1 = [];
    S1.montage = montage;
    S1.prefix  = S.prefix; % ignored for online
    S1.keepsensors = false;
    S1.keepothers  = false;

    switch S.mode
        case 'write'
            S1.mode    = 'write';
        case 'online'
            S1.mode = 'switch';
        case 'onlinecopy'
            S1.mode    = 'switch';
            D = copy(D, [S.prefix D.fname]);
    end

    S1.D = D;
    D = spm_eeg_montage(S1);
end

res.files = {fullfile(D)};

author: spm


text:
function []=panel5cfdisp(N,n,T,Units,endo,Ymat,stringdates2,decimaldates2,Fstartlocation,Fendlocation,cforecast_estimates,pref)

% preliminary task: reshape Ymat
Ymat=reshape(Ymat,T,n,N);
% preliminary task: gather in a cell the values to be plotted
% initiate the cell
plotdata={};
% because forecasts have to be computed for each unit, loop over units
for ii=1:N
% each cell entry is a matrix of actual and forecast values
% this matrix comprises 4 rows: the first row is actual data, while the three other rows are the estimates (point estimates and confidence bands) for the forecasts
% also, this matrix has a number of rows equal to the dimension of decimaldates2, which comprises the total period sample+forecasts
   % loop over variables
   for jj=1:n
      plotdata{jj,1,ii}=nan(4,size(decimaldates2,1));
      % record actual sample values
      plotdata{jj,1,ii}(1,1:T)=Ymat(:,jj,ii)';
      % copy the last point of the actual sample for the forecast part of the matrix (required to have a clean plot)
      plotdata{jj,1,ii}(:,Fstartlocation-1)=repmat(Ymat(Fstartlocation-1,jj,ii),4,1);
      % record forecast, lower bound
      plotdata{jj,1,ii}(2,Fstartlocation:Fendlocation)=cforecast_estimates{jj,1,ii}(1,:);
      % record forecast, point estimate
      plotdata{jj,1,ii}(3,Fstartlocation:Fendlocation)=cforecast_estimates{jj,1,ii}(2,:);
      % record forecast, upper bound
      plotdata{jj,1,ii}(4,Fstartlocation:Fendlocation)=cforecast_estimates{jj,1,ii}(3,:);
   end
end

% then plot the figure
if pref.plot
forecast=figure('Tag','BEARresults');
set(forecast,'Color',[0.9 0.9 0.9]);
set(forecast,'name','conditional forecasts');
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % increment count
   count=count+1;
   % then plot
   subplot(N,n,count)
   hold on
   Xpatch=[decimaldates2(Fstartlocation-1:Fendlocation,1)' fliplr((decimaldates2(Fstartlocation-1:Fendlocation,1))')];
   Ypatch=[plotdata{jj,1,ii}(2,Fstartlocation-1:Fendlocation) fliplr(plotdata{jj,1,ii}(4,Fstartlocation-1:Fendlocation))];
   Fpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(Fpatch,'facealpha',0.6);
   set(Fpatch,'edgecolor','none');
   plot(decimaldates2,plotdata{jj,1,ii}(3,:),'Color',[0.4 0.4 1],'LineWidth',2);
   plot(decimaldates2,plotdata{jj,1,ii}(1,:),'Color',[0 0 0],'LineWidth',2);
   hold off
   set(gca,'XLim',[decimaldates2(1,1) decimaldates2(end,1)],'FontName','Times New Roman');
   set(gca,'XGrid','on');
   set(gca,'YGrid','on');
      % top labels
      if count<=n
      title(endo{count,1},'FontWeight','normal');
      end
      % side labels
      if jj==1
      ylabel(Units{ii,1},'FontWeight','normal');
      end
   end
end
end

% save on Excel
% create the cell that will be saved on excel
cforecastcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},size(stringdates2,1)+3,1);
horzspace=repmat({''},3,6*n);
% loop over units
for ii=1:N
% initiate the cell of results
unitcell={};
   % loop over endogenous variables (horizontal dimension)
   for jj=1:n
   % create a header
   header=[{[Units{ii,1} ': ' endo{jj,1}]} {''} {''} {''} {''};{''} {''} {''} {''} {''};{''} {'actual'} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   endocell=[[header;stringdates2 num2cell((plotdata{jj,1,ii})')] vertspace];
   % concatenate to the previous parts of unitcell
   unitcell=[unitcell endocell];
   end
% concatenate to the previous parts of afcell
cforecastcell=[cforecastcell;horzspace;unitcell];
end
% trim
cforecastcell=cforecastcell(4:end,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),cforecastcell,'cond forecasts','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
% read_rdf() - read RDF-formatted EEG files.
%
% Usage:
%   >> [eeg,ev,header] = read_rdf(filename);
%
% Inputs:
%   filename - EEG data file in RDF format
%
% Outputs:
%   eeg - eeg data (array in size of [chan_no timepoint];
%   ev  - event structure
%      ev.sample_offset[] - event offsets in samples
%                           from the first sample (0)
%      ev.event_code[]    - event codes (integers)
%   header - data structure for header information
%      header.ch_no     - number of channels
%      header.sample_no - number of samples

% Notes:
%
% Authors: Jeng-Ren Duann, CNL/Salk & INC/UCSD, 2002-12-12
%          with help from Andrey Vankov, creator of the RDF file format.

function [eeg,ev,header] = read_rdf(filename)

    if nargin < 1
        help read_rdf;
        return;
    end;

  eeg = [];
  ev = [];
  header = [];

  fp = fopen(filename,'rb','ieee-le');
  if fp == -1,
    disp('read_RDF(): Cannot open data file...!');
    return;
  end

  fseek(fp,6,-1);
  header.ch_no = fread(fp,1,'uint16');

  cnt = 0;
  ev_cnt = 0;
  while(~feof(fp)),
    tag = fread(fp,1,'uint32');
    if length(tag) == 0,
      break;
    end
    if tag == hex2dec('f0aa55'),
      cnt = cnt + 1;
      disp(['block ' num2str(cnt) ' found']);
      % read ch_no and block length
      fseek(fp,2,0);
      ch_no = fread(fp,1,'uint16');
      block_size = power(2,fread(fp,1,'uint16'));
      % read events
      fseek(fp,62,0);
      for i=1:110,
	samp_off = fread(fp,1,'char');
	cond_code = fread(fp,1,'char');
	ev_code = fread(fp,1,'uint16');
	if samp_off ~= 0,
	  ev_cnt = ev_cnt + 1;
	  ev(ev_cnt).sample_offset = samp_off + (cnt-1)*128;
	  ev(ev_cnt).event_code = ev_code;
	end
      end
      data = fread(fp,ch_no*block_size,'int16');
      data = reshape(data,ch_no,block_size);
      eeg = [eeg data];
    end
  end

  fclose(fp);
  header.sample_no = size(eeg,2);


text:
function batchImport(source,destn,zipFlag,mergeScansFlag)
%---------------------------------------------
%INPUTS
%source = 'path/to/source';
%destn = 'path/to/destn';
%zipFlag = 'No';
%mergeScansFlag = 'No';

%NOTE:
%1. To import selected list of folders:
%source: dirsToImportC (cell array of paths to directories to be imported)
%2.To open a pool of parallel processors:
%pool = parpool(4);

if iscell(source)

    dirsToImportC = source;

else

    % Define file and directory filter patterns
    fileFilterC = {'dose','scout','cor','sag','report','_nac','mip',...
        'coronal','sagittal','screen_save'};
    fileFilterC = {};
    dirFilterC = {'ot'};

    % Get all directories and files
    tic,
    [filesInCurDir,dirsInCurDir] = rdir(source);
    toc

    % Convert to lower case
    dirNameC = {dirsInCurDir.name};
    dirFullPathC = {dirsInCurDir.fullpath};

    % determine back or forward slash
    if ispc
        slashType = '\';
    else
        slashType = '/';
    end

    % Filter out files
    indV = zeros(1,length(dirNameC));
    for i = 1:length(fileFilterC)
        indC = strfind(dirNameC,fileFilterC{i});
        indV1 = ~cellfun(@isempty, indC);
        indV = indV | indV1;
        sum(indV)
    end

    % Filter out directories
    for i = 1:length(dirFilterC)
        indC = strfind(dirFullPathC,[slashType,dirFilterC{i},slashType]);
        indV1 = ~cellfun(@isempty, indC);
        indV = indV | indV1;
    end

    % list of directories after filtering NAC's, etc.
    dirsToImportC = dirFullPathC(~indV);

    % filter directories containing no files
    indV1 = false(length(dirsToImportC),1);
    tic,
    for dirNum = 1:length(dirsToImportC)
        dirS = dir(dirFullPathC{dirNum});
        if ~all([dirS.isdir])
            indV1(dirNum) = 1;
        end
    end
    toc
    dirsToImportC = dirsToImportC(indV1);
end

%% Import DICOM to CERR
tic
hWaitbar = NaN;
% Import all the dirs
for dirNum = 1:length(dirsToImportC)
    try
        init_ML_DICOM
        %hWaitbar = waitbar(0,'Scanning Directory Please wait...');
        sourceDir = dirsToImportC{dirNum};
        selected = 'all';
        if strcmpi(selected,'all')
            combinedDcmdirS = getCombinedDir(sourceDir);
            % Pass the java dicom structures to function to create CERR plan
            try
                optName = fullfile(getCERRPath,'CERROptions.json');
                optS    = opts4Exe(optName);
                planC = dcmdir2planC(combinedDcmdirS,mergeScansFlag,optS);
            end
        else
            patient = scandir_mldcm(sourceDir, hWaitbar, 1);
            %close(hWaitbar);
            dcmdirS = struct(['patient_' num2str(1)],patient.PATIENT(1));
            for j = 2:length(patient.PATIENT)
                dcmdirS.(['patient_' num2str(j)]) = patient.PATIENT(j);
            end
            patNameC = fieldnames(dcmdirS);
            planC = dcmdir2planC(patient.PATIENT,mergeScansFlag);
        end

        indexS = planC{end};

%         % build the filename for storing planC
%         mrn = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.PatientID;
%         studyDscr = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.StudyDescription;
%         seriesDscr = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.SeriesDescription;
%         modality = planC{indexS.scan}(1).scanInfo(1).DICOMHeaders.Modality;

%         outFileName = [mrn,'~',studyDscr,'~','_FSPost_',seriesDscr,'~',modality];
%         outFileName = strrep(outFileName,'\','-');
%         outFileName = strrep(outFileName,'/','-');
%         outFileName = strrep(outFileName,':','-');
        %outFileName = mrn;   % store file names as MRNs
        %------------TMEP FIX FOR FILES WITH dots in name-------------
        [~,outFileName,endbit] = fileparts(sourceDir);  % store file names as DICOM directory names
        outFileName = [outFileName,endbit];
        %-----------------------------------------------------
        fullOutFileName = fullfile(destn,outFileName);

        %Check for duplicate name of fullOutFileName
        dirOut = dir(destn);
        allOutNames = {dirOut.name};
        if any(strcmpi([outFileName,'.mat'],allOutNames))
            fullOutFileName = [outFileName,'_duplicate_',num2str(rand(1))];
        end
        if strcmpi(zipFlag,'Yes')
            saved_fullFileName = fullfile(destn,[outFileName,'.mat.bz2']);
        else
            saved_fullFileName = fullfile(destn,[outFileName,'.mat']);
        end
        if ~exist(fileparts(saved_fullFileName),'dir')
            mkdir(fileparts(saved_fullFileName))
        end
        save_planC(planC,[], 'passed', saved_fullFileName);

    catch

        disp(['Cannot convert ',sourceDir])

    end
end

toc
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%                                            Lden = symbcholden(L,dense,DAt)
% SYMBCHOLDEN  Creates Lden.{LAD, perm,dz, sign, first}
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, dpr1fact

function Lden = symbcholden(L,dense,DAt)
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

% ------------------------------------------------------------
% Symbolic forward Cholesky of dense columns, in order
% [LP, Q-blk, Q-norm, Q-tr]
% ------------------------------------------------------------
i1 = dense.l + 1;
i2 = i1 + length(dense.q);
LAD = [symbfwblk(L,dense.A(:,1:i1-1)), symbfwblk(L,DAt.denq),...
    symbfwblk(L,dense.A(:,i2:end)),symbfwblk(L,dense.A(:,i1:i2-1))];
% ------------------------------------------------------------
% Incremental ordering heuristic, excluding the Lorentz-trace cols
% ------------------------------------------------------------
[perm, dz] = incorder(LAD(:,1:length(dense.cols)));
% ------------------------------------------------------------
% Insert the trace cols with a "-1"-factor just after corresponding
% Lorentz-block columns
% ------------------------------------------------------------
Lden = finsymbden(LAD,perm,dz,i1);

repo: REST


text:
function imagenet_subset(dataset, subset)
% IMAGENET_SUBSET Modifies a ImageNet dataset to use a subset of classes

  rng(0);
  if isnumeric(subset)
    if isscalar(subset)
      % pick random subset with given size
      idx = randperm(1000, subset);
    else  % passed it in directly
      idx = subset;
    end
  else
    % named subsets
    switch subset
    case 'fruits'
      idx = [323, 320, 324, 229, 322, 326, 321, 318, 331, 319];
    case 'animals'
      idx = [1347, 1123, 1161, 1276, 1215];
    otherwise
      error('Unknown subset.');
    end
  end

  % index to synset mapping
  synsets = [];
  load([dataset.dataDir '/ILSVRC2014_devkit/data/meta_clsloc.mat'], 'synsets');

  if ischar(idx)
    % find classes by name
    subset_names = strtrim(strsplit(lower(idx), ','));
    synset_names = {synsets.words};

    valid = false(numel(synset_names), 1);
    for i = 1:numel(synset_names)
      % each synset name contains multiple synonims, separated by commas
      synonyms = strtrim(strsplit(lower(synset_names{i}), ','));
      valid(i) = any(ismember(subset_names, synonyms));
    end

    idx = find(valid);
    assert(numel(idx) == numel(subset_names), 'At least one class was not found.');
  end

  % synset IDs (e.g. 'n01542...', a string) of the selected classes
  synset_ids = {synsets(idx).WNID};

  % binary vector of which synsets/classes to keep
  valid_classes = ismember(dataset.classes.name, synset_ids);

  assert(nnz(valid_classes) == numel(synset_ids), 'At least one class was not found.');

  % remove all other classes, and their samples
  invalid_samples = ~ismember(dataset.labels, find(valid_classes));
  dataset.labels(invalid_samples) = [];
  dataset.sets(invalid_samples) = [];
  dataset.filenames(invalid_samples) = [];
  dataset.augmentImage(invalid_samples) = [];

  dataset.trainSet = find(dataset.sets == 1);
  dataset.valSet = find(dataset.sets == 2);

  % remap labels to 1:N
  map = cumsum(double(valid_classes));
  dataset.labels = map(dataset.labels);
end

repo: curveball


text:
function Q = spm_dcm_csd_Q(csd)
% Precision of cross spectral density
% FORMAT Q = spm_dcm_csd_Q(csd)
%
% csd{i}   - [cell] Array of complex cross spectra
% Q        - normalised precision
%--------------------------------------------------------------------------
% This routine returns the precision of complex cross spectra based upon
% the asymptotic results described in Camba-Mendez & Kapetanios (2005):
% In particular, the scaled difference between the sample spectral density
% (g) and the predicted density (G);
%
% e = vec(g - G)
%
% is asymptotically complex normal, where the covariance between e(i,j) and
% e(u,v) is given by Q/h and:
%
% Q = G(i,u)*G(j,u):  h = 2*m + 1
%
% Here m represent the number of averages from a very long time series.The
% inverse of the covariance is thus a scaled precision, where the
% hyperparameter (h) plays the role of the degrees of freedom (e.g., the
% number of averages comprising the estimate). In this routine, we use the
% sample spectral density to create a frequency specific precision matrix
% for the vectorised spectral densities - under the assumption that the
% former of this sample spectral density resembles the predicted spectral
% density (which will become increasingly plausible with convergence).
%
% Camba-Mendez, G., & Kapetanios, G. (2005). Estimating the Rank of the
% Spectral Density Matrix. Journal of Time Series Analysis, 26(1), 37-48.
% doi: 10.1111/j.1467-9892.2005.00389.x
%__________________________________________________________________________
% Copyright (C) 2018 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_csd_Q.m 7751 2019-12-06 11:59:09Z peter $

%-Check for cell arrays
%--------------------------------------------------------------------------
if iscell(csd)
    CSD   = spm_zeros(csd{1});
    n     = numel(csd);
    for i = 1:n
        CSD = CSD + csd{i};
    end
    Q  = spm_dcm_csd_Q(CSD/n);
    Q  = kron(eye(n,n),Q);
    return
end

%-Get precision
%--------------------------------------------------------------------------
SIZ     = size(csd);
Qn      = spm_length(csd);
Q       = spalloc(Qn,Qn,SIZ(1)*prod(SIZ(2:end))^2);
[w,i,j] = ind2sub(SIZ,1:Qn);
for Qi  = 1:Qn
    for Qj = 1:Qn
        if w(Qi) == w(Qj)
            Q(Qi,Qj) = csd(w(Qi),i(Qi),i(Qj))*csd(w(Qi),j(Qi),j(Qj));
        end
    end
end
Q       = inv(Q + norm(Q,1)*speye(size(Q))/32);