text:
function hhof = plotHofFeatures(H,hax,maxv,colors,binedges,psize,offset)

if nargin<7,
  offset = [0 0];
end

hhof = [];
for xi = 1:size(H,2),
  cx = (psize/2 + 1 + (xi-1)*psize);
  for yi = 1:size(H,1),
    cy = (psize/2 + 1 + (yi-1)*psize);

    for bini = 1:size(H,3),
      tmp = linspace(binedges(bini),binedges(bini+1),20);
      xcurr = cx + [0,psize/2*cos(tmp),0];
      ycurr = cy + [0,psize/2*sin(tmp),0];
      hhof(yi,xi,bini) = patch(xcurr+offset(1),ycurr+offset(2),colors(bini,:),...
        'LineStyle','none','FaceAlpha',min(1,H(yi,xi,bini)/maxv),'Parent',hax);
    end

  end
end

author: kristinbranson

repo: JAABA


text:
function writeSmoothAbs()
%
% FUNCTION:
%   This function is used to write the code for smoothAbs
%

MaxOrder = 8;  %Should be at least 1
N = MaxOrder + 1;

%% Set up the constraints:

% The value and slope should be 1 at the right transition point, and zero
% for any higher derivatives.
xRight = ones(1,N);
yRight = zeros(1,N);
yRight(1) = 1; yRight(2) = 1;

% The value and slope should be 1 at the right transition point, and zero
% for any higher derivatives.
xLeft = -ones(1,N);
yLeft = zeros(1,N);
yLeft(1) = 1; yLeft(2) = -1;

% All Constraints are expressed across all derivatives
dLevel = 0:N;

%% Begin writing the file:

fid = fopen('smoothAbs.m','w');

fprintf(fid,'function y = smoothAbs(x,alpha,order)\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% y = smoothAbs(x,alpha,order)\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% FUNCTION:\n');
fprintf(fid,'%%   This function is a smooth version of abs(). Smoothing is done using\n');
fprintf(fid,'%%   a piecewise polynomial, which is Nth order continuous at the\n');
fprintf(fid,'%%   transitions, where N is a value between 0 and 8, and input as ''order''\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% INPUTS:\n');
fprintf(fid,'%%   x = a vector of matrix of real numbers\n');
fprintf(fid,'%%   alpha = smoothing parameter: -alpha<x<alpha will be smoothed\n');
fprintf(fid,'%%   order = the order of the smoothing (0=no smoothing)\n');
fprintf(fid,['%%           order must be in the set {0,1,..,' num2str(MaxOrder) '}\n']);
fprintf(fid,'%%\n');
fprintf(fid,'%% OUTPUTS:\n');
fprintf(fid,'%%   y = a smoothed version of abs(x)\n%%\n\n');
fprintf(fid,'switch order \n');
fprintf(fid,'    case 0 %% No Smoothing \n');
fprintf(fid,'        y=abs(x); \n');

for i=2:N
    x = [xLeft(1:i),xRight(1:i)];
    y = [yLeft(1:i),yRight(1:i)];
    d = [dLevel(1:i),dLevel(1:i)];
    [~,pSym] = SolveCoeff(x,y,d);
    str = 'p = [';
    for j=2:length(pSym)  %Leading coefficient is always zero
        if j==length(pSym)
            str = [str char(pSym(j)) ']; \n']; %#ok<*AGROW>
        else
            str = [str char(pSym(j)) ','];
        end
    end
    fprintf(fid,['    case ' num2str(i-1) ' \n']);
    fprintf(fid,'        y=x; \n');
    fprintf(fid,'        c1 = x < -alpha; \n');
    fprintf(fid,'        c2 = x > alpha; \n');
    fprintf(fid,'        y(c1)=-x(c1); \n');
    fprintf(fid,['        ' str]);
    fprintf(fid,'        idx = ~c1&~c2;\n');
    fprintf(fid,'        y(idx)=alpha*polyval(p,x(idx)/alpha); \n');
end

fprintf(fid,'    otherwise \n');
fprintf(fid,'        error(''Order not supported'');\n');
fprintf(fid,'end \n');

fprintf(fid,'end \n');

fclose(fid);
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [beta_draws,sigma_draws,IV_draws,C_draws,D,gamma,storage1,storage2]=irfIV6(betahat,n,m,p,k,T,names,startdate,enddate,X,Y,endo,data_endo,data_exo,const,pref,strctident,IRFt,IRFperiods,It,Bu)

% inputs:  - matrix 'betahat': OLS estimate for beta
%          - matrix 'sigmahat': OLS estimate for sigma
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - cell 'signrestable': table recording the sign restriction input from the user
%          - cell 'signresperiods': table containing the periods corresponding to each restriction
%          - string 'ShockwithInstrument' Name of the shock where the instrument belongs to
% outputs: - cell 'struct_irf_record': record of the draws for the orthogonalised IRFs
%          - matrix 'D_record': record of the draws for the structural matrix D
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - matrix 'hd_record': record of historical decompositions
%          - matrix 'ETA_record': record of structural Shocks
%          - vector 'IVcorrelation': vector of correlation with proxy proxy
%          - vector 'PofIVcorrelation': corresponding P-value

%% Phase 1:% %% IV F-Test
% Load IV and make it comparable with the reduced form errors
[EPSIV,IVcut,EPSt,sigmahatIV,B,EPS,df,sigma_hat,inv_sigma_hat,beginInstrument,EndInstrument,IV,txt,OverlapIVinY,cut1,cut2,cut3,cut4]=bear.loadIV(betahat,k,n,Y,X,T,p,names,startdate,enddate,strctident,pref);

%%  Check strength of instrument
F_stat_test_unc   = NaN*ones(1,n);
signif_test_unc   = NaN*ones(1,n);
Betas_test_unc    = NaN*ones(1,n);
Rsquared_test_unc = NaN*ones(1,n);

for i = 1:n
    yyy      = EPSt(i,:)'; %%get residuals
    xxx_rest = ones(size(EPSt,2),1); %%vector of constants
    xxx      = [ones(size(EPSt,2),1), IVcut]; %%constants plus instrument

    % restricted model
    coeff_rest = (xxx_rest'*xxx_rest)^(-1)*xxx_rest'*yyy;
    resid_rest = yyy - xxx_rest*coeff_rest;
    RSS_rest   = resid_rest'*resid_rest;

    % unrestricted model
    coeff = (xxx'*xxx)^(-1)*xxx'*yyy;
    resid = yyy - xxx*coeff;
    RSS   = resid'*resid;

    % compute statistic
    F_stat   = (RSS_rest-RSS)*(size(yyy,1)-size(xxx,2))/(RSS*(size(xxx,2)-size(xxx_rest,2)));
    F_pvalue = 1-fcdf(F_stat,size(xxx,2)-size(xxx_rest,2),size(yyy,1)-size(xxx,2));

    F_stat_test_unc(i) = F_stat;
    significance = 3*(F_pvalue < 0.01);
    significance = significance + 2*(F_pvalue < 0.05 & F_pvalue > 0.01);
    significance = significance + 1*(F_pvalue < 0.1 & F_pvalue > 0.05);
    signif_test_unc(i) = significance;

    % regression statistics
    coeff  = (xxx'*xxx)^(-1)*xxx'*yyy;
    fitted = xxx*coeff;
    resid  = yyy - fitted;
    yyy_demean           = yyy - mean(yyy);
    Rsquared_test_unc(i) = 1 - resid'*resid/(yyy_demean'*yyy_demean);

    Betas_test_unc(i) = coeff(2);
end

%% Print F Statistic

table_validitytest_rows = char('coefficient', 'significance level', 'number of obs', 'F statistic', 'R squared');
table_validitytest_rows = cellstr(table_validitytest_rows);
table_validitytest_values  = [Betas_test_unc;  signif_test_unc;  length(IVcut)*ones(1,n); F_stat_test_unc; Rsquared_test_unc ];
table_validitytest_columns = cellstr(endo)';

table = array2table(table_validitytest_values);
for kk=1:length(table_validitytest_columns)
    table.Properties.VariableNames(kk) = table_validitytest_columns(1,kk);
end
table_rows = cell2table(table_validitytest_rows);
table=[table_rows,table];

for kk=1
    table.Properties.VariableNames(kk) = cellstr(char('Validitytest'));
end

% print proxy var statistics in the command window and in the results file
filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');
%two empty lines
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
%print Proxy VAR statistics title
fprintf('%s\n','Proxy VAR statistics');
fprintf(fid,'%s\n','Proxy VAR statistics');
fprintf('%25s %15s %15s %15s %15s\n','','coefficient','significance','number of obs','F statistic');
fprintf(fid,'%25s %15s %15s %15s %15s\n','','coefficient','significance','number of obs','F statistic');
% handle the endogenous
   for jj=1:n
      values=[Betas_test_unc(1,jj);  signif_test_unc(1,jj);  length(IVcut); F_stat_test_unc(1,jj)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1}),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1}),values);
   end
%empty line
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fclose(fid);

%% first stage regression

%step 2: Regress the first reduced form shock on the instrument
Shock = EPSIV(:,1);
[nobs,~] = size(IVcut);
XX = [ones(nobs,1) IVcut];
[~,nvar] = size(XX);
%get OLS estimate
XpXi = (XX'*XX)\eye(nvar);
betaIV=XpXi*(XX'*Shock);
%get predicted value
IVpred = XX*betaIV;

%% second stage regression
%step 3: Regress the other reduced form shocks on the predicted value
ImpactIRFIV = zeros(n,1);
ImpactIRFIV(1,1) = 1;

for hh=2:n
Shock = EPSIV(:,hh);
[nobs,~] = size(IVpred);
IVpredtemp = [ones(nobs,1) IVpred];
[~,nvar] = size(IVpredtemp);
IVpIVi = (IVpredtemp'*IVpredtemp)\eye(nvar);
betaIV2=IVpIVi*(IVpredtemp'*Shock);
ImpactIRFIV(hh,1) = betaIV2(2,1); %should be equal to Mu from 2:end
end

% step 5: Create the structural matrix and only fill the first column as
% this is the only one identified
%Dols=eye(n,n);
%%another way to retrieve b11 (the scalar that scales the IRF to be a 1sdt Shock) is simply
C=chol(bear.nspd(sigmahatIV),'lower');

if IRFt==5 % D & gamma output
    % step 5: Create the structural matrix and only fill the first column as
% this is the only one identified
D=eye(n,n);
% Step 6: Replace the first Column in the Cholesky Decomposition by
%the structural impact matrix computed above
% D(1:end,1) = ImpactIRFIV*b11;

gamma=eye(n);
%%another way to retrieve b11 (the scalar that scales the IRF to be a 1sdt Shock) is simply
b=ImpactIRFIV;
%%Recover the vector q that maps the first column of C into b such that Cq=b;
q = C\b;
%%b11 is the euclidian length of q
b11q = 1/norm(q);
if strctident.stdeviation ==1
D(1:end,1) = ImpactIRFIV*b11q;
else
D(1:end,1) = ImpactIRFIV;
end
end

%% Preliminiaries for sign restrictions
Acc=It-Bu;  %Minimum accepted draws for the rotation matrix

%% Storage cells and other preliminiaries for the bootstrap part
endo_artificial = zeros(T+p,n); %initialize artifical data matrix
%preliminaries for IV block bootstrap
BlockSize = p; %this can be made frequency specific %only used for bootstraptype2
nBlock = ceil((T)/BlockSize); %only used for bootstraptype2
startofcommonsample= find(OverlapIVinY==1, 1, 'first');
endofcommonsample= find(OverlapIVinY==1, 1, 'last');
%EPShatcut = EPS(OverlapIVinY,:);
beta_draws = nan(k*n,Acc*10);
sigma_draws=nan(n^2,Acc*10);
IV_draws=nan(n,Acc*10);
C_draws=nan(n^2,Acc*10);
beta = betahat;
B=reshape(beta,k,n);

%%%%while AA<=BB with BB=1000;
%%%%%%%%%%%%%%%%%%%%Bootstrap to approximate distribution of reduced form and IV coefficients%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hbar=bear.parfor_progressbar(Acc*10,'IV progress');
for AA=1:Acc*10 %parfor loop possible here? %%%% are 100000 draws necessary? parfor AA=1:Acc*100
    stationary=0;
while stationary==0
        %%Get reduced form coefficients and IV column
        %% STEP 1: generate the artificial data
        if strctident.bootstraptype==1
            % Initialize wild bootstrap by drawing residuals and flipping the sign at
            % random get reduced form residuals
            EPS  = Y-X*B;
            rotationvector = 1-2*(rand(T,1)>0.5); %the rotation vector is of the same size as EPS. We need to determine where EPS starts in terms of IV and where it ends
            EPSrotate = EPS.*(rotationvector*ones(1,n));
            %create IV specific rotation vector that is 1 if the samples dont
            %overlapp and the EPS rotation vector otherwhise
            %         IVrotationvector = rotationvector(OverlapIVinY);
            %         IVrotate = IV.*IVrotationvector;

            IVrotationvector = ones(length(IV),1);
            for kk=startofcommonsample:(endofcommonsample-startofcommonsample)
                if OverlapIVinY(kk,1)==1
                    IVrotationvector(kk,1) = rotationvector(kk,1);
                else
                    IVrotationvector(kk,1) = 1;
                end
            end
            IVrotate = IV.*IVrotationvector;

        else
            %%Moving block bootstrap bootstrap
            EPS  = Y-X*B;
            IVblocks = zeros(BlockSize,size(IV,2),length(IV)-BlockSize+1); %create moving blocks for the instrument
            for yy = 1:length(IV)-BlockSize+1
                IVblocks(:,:,yy) = IV(yy:BlockSize+yy-1,:);
            end

            EPSblocks = zeros(BlockSize,n,T-BlockSize+1); %create moving blocks for the reduced form errors
            for yy = 1:T-BlockSize+1
                EPSblocks(:,:,yy) = EPS(yy:BlockSize+yy-1,:);
            end

            %the first block that can be accomodated by the instrument (i.e. where reduced form errors and instrument overlap
            %is given by EPSblocks(:,:,startofcommonsaple);

            %the last block of reduced form errors that can be accomodated by
            %the instrument is EPSblocks(:,:,endofcommonsample-BlockSize);

            %preliminaries for centering of the bootstrapped VAR errors
            centering = zeros(BlockSize,n);
            for yy = 1:BlockSize
                centering(yy,:) = mean(EPS(yy:T-BlockSize+yy,:,1),1); %running mean
            end

            centering = repmat(centering,[nBlock,1]);
            centering = centering(1:T,:);

            %center the bootstrapped proxy variables
            IVcentering = zeros(BlockSize,size(IV,2));
            for yy = 1:BlockSize %moving average of the instrument
                subIV = IV(yy:size(IV,1)-BlockSize+yy,:);
                IVcentering(yy,:) = mean(subIV((subIV(:,1) ~= 0),1),1); %running mean
            end
            IVcentering = repmat(IVcentering,[nBlock,1]);
            IVcentering = IVcentering(1:size(IV,1),:);

            pickblock = randi([startofcommonsample endofcommonsample-BlockSize+1],nBlock,1);
            pickblockIV = pickblock-startofcommonsample+1;
            %now check if we can replace some blocks in the beginning of the sample
            %with other residuals, that are not accomodated by the instrument
            missing_beginning = startofcommonsample;
            %missing_end = T-endofcommonsample;

            blocks_beginning = floor(missing_beginning/BlockSize); %number of blocks that we can possibly replace in the beginning
            blocks_end = floor((T-endofcommonsample)/BlockSize); %number of blocks that we can possibly replace in the end

            if blocks_beginning > 0 %if the instrument and the sample dont coincide for at least the length of the first block
                pickblock_beginning = randi([1 startofcommonsample],blocks_beginning,1);
                pickblock(1:blocks_beginning,:)=pickblock_beginning(1:blocks_beginning);
            end

            if blocks_end > 0 %if we can resample a block from the blocks that are not accomodated by the instrument do so
                pickblock_end = randi([endofcommonsample-BlockSize+2 T-BlockSize+1],blocks_end,1);
                howmanyblock = length((nBlock-blocks_end) : nBlock) ;
                if howmanyblock > length(pickblock_end)
                pickblock((nBlock-blocks_end)+1 : nBlock)=pickblock_end(1:blocks_end);
                else
                pickblock((nBlock-blocks_end) : nBlock)=pickblock_end(1:blocks_end);
                end
            end

            EPSrotate = zeros(nBlock*BlockSize,n);
            IVrotate = zeros(nBlock*BlockSize,size(IV,2));
            for yy = 1:nBlock
                EPSrotate(1+BlockSize*(yy-1):BlockSize*yy,:) = EPSblocks(:,:,pickblock(yy,1));
                IVrotate(1+BlockSize*(yy-1):BlockSize*yy,:) = IVblocks(:,:,pickblockIV(yy,1));
            end

            EPSrotate = EPSrotate(1:T,:); %drop last pseudo observations
            IVrotate = IVrotate(1:T,:); %drop last pseudo observations
%             IVrotate = IVrotate(1:length(IV),:);
            IVrotate = IVrotate(OverlapIVinY,:); %drop observations where instrument is not available
            %EPScutrotate = EPSrotate(OverlapIVinY,:);
            %center the bootstrapped residuals and proxies by subtracting
            %the mean
            EPSrotate = EPSrotate - centering;
            for yy = 1:size(IVrotate,2)
                IVrotate((IVrotate(:,yy)~=0),yy) = IVrotate((IVrotate(:,yy)~=0),yy) - IVcentering((IVrotate(:,yy)~=0),yy);
            end

        end

        %% STEP 2: Generate artificial data
        [endo_artificial_gen]=bear.gen_fake_data(endo_artificial,data_endo,B,EPSrotate,const,p,n,T);

        %% STEP 3: estimate reduced form VAR on artificial data.
        [~,betadraw,sigmadraw,Xdraw,~,Ydraw]=bear.olsvar(endo_artificial_gen,data_exo,const,p);
        %%%%% do we need this test for irft6 only? or is it also useful for
        %%%%% IRFt5
        if IRFt==5
        stationary=1;
        elseif IRFt==6
        [stationary]=bear.checkstable(betadraw,n,p,k); %switches stationary to 0, if the draw is not stationary
        end
end
        if IRFt==5
            %% STEP 4: identify the model
    [Ddraw]=bear.irfiv_ols_for_bootstrap_GK(txt,beginInstrument,EndInstrument,names,IVrotate,betadraw,sigmadraw,m,n,Xdraw,Ydraw,k,p,enddate,startdate,endo,IRFperiods,pref,cut1,cut2,cut3,cut4,strctident);
    %% Step 5: Calculate impulse responses and store them
    [~,ortirfmatrix]=bear.irfsim(betadraw,Ddraw,n,m,p,k,IRFperiods);
       for jj=1:IRFperiods
       storage1{AA,1}(:,:,jj)=ortirfmatrix(:,:,jj);
       end
       storage2{AA,1}=Ddraw;
%        beta_gibbs(:,AA)=betadraw;
%        sigma_gibbs(:,AA)=bear.vec(sigmadraw);
        elseif IRFt==6
        %% Step 4: Identifity the proxy VAR equation
        [~,~,q,C]=irfiv_ols_for_ivres(txt,names,IVrotate, betadraw,sigmadraw,m,n,Xdraw,Ydraw,k,p,enddate,startdate, endo, cut1, cut2,cut3,cut4);
beta_draws(:,AA)=betadraw;
sigma_draws(:,AA)=bear.vec(sigmadraw);
%D_draws(:,AA)=bear.vec(D);
IV_draws(:,AA)=q/norm(q);
C_draws(:,AA)=bear.vec(C);
%IVrotate_draws(:,AA)=IVrotate;
        end
%%%AA = AA+1 ;
hbar.iterate(1); % update progress by one iteration
end
close(hbar);   %close progress bar

%%%%statistics here?

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = plotcoeffs(f, varargin)
%PLOTCOEFFS   Display coefficients graphically.
%   PLOTCOEFFS(F) calls PLOTCOEFFS(F.SMOOTHPART).
%
%   This method is mainly aimed at Chebfun developers.
%
% See also PLOT.

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

% Deal with an empty input:
if ( isempty(f) )
    if ( nargout == 1 )
        [varargout{1:nargout}] = plot([]);
    end
    return
end

% Call PLOTCOEFFS() on the SMOOTHPART of F.
[varargout{1:nargout}] = plotcoeffs(f.smoothPart, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% input: uv_c1 N x 2 point array in c1 frame,
% output: uv_c2 N by 2 matrix in C2 image coordinates
% R_c2c1 the rotation from C1 frame to C2 frame, t_c2c1 the translation
% from c1 frame to c2 frame, it is the position of the camera c1 origin in
% c2 frame, n is the average normal direction of the feature point in the
% C1 frame, d is the distance from the feature to the c1 optical center
% along the minus n direction
% for normalize() input is 2 x N matrix, output is 2 by N matrix,
% project_points2() input 3 by N matrix, output 2 by N matrix
function uv_c2=rotate_with_dist_fc_c2c1(cam, uv_c1, R_c2c1, t_c2c1, n, d)

xn=normalize(uv_c1',cam.fc,cam.cc,cam.kc,cam.alpha);
xn2=((R_c2c1-(t_c2c1*n'/d)))*[xn;ones(1,size(xn,2))];
uv_c2=project_points2(xn2,zeros(3,1),zeros(3,1),cam.fc,cam.cc,cam.kc,cam.alpha);
uv_c2=uv_c2';

repo: ekfmonoslam


text:
%  qtkp_shift_solve -- Solves the problem (alpha@T1@T2@...@Tp-lambda*I)*y=c
%  where @ stands for the kronecker product for quasi-triangular matrices Ti
%
%  ::
%
% 	y = qtkp_shift_solve(T,n,c,lambda)
% 	y = qtkp_shift_solve(T,n,c,lambda,alpha)
% 	y = qtkp_shift_solve(T,[],...)
%
%  Args:
%
%     - **T** [cell array]: each element Ti, i=1,2,...,p is a
%       quasi-triangular (square) matrix
%
%     - **n** [vector|{[]}]: number of rows (or columns) for Ti,
%       i=1,2,...,p.
%
%     - **c** [vector]: right-hand side of the equality
%
%     - **lambda** [scalar]: See formula above
%
%     - **alpha** [scalar|2 x 2 matrix|{1}]: First element of the kronecker
%       product
%
%  Returns:
%     :
%
%     - **y** [N x 1 vector]: where N=n1 x n2 x ... x np
%
%  See also : utils.kronecker.kp_shift_solve
%

author: jmaih


text:
function sqi = stdsqi(signal)
%stdSQI Standard Deviation SQI
%
% Returns the normalized standard deviation of a signal on the upper
% 5% of the signals amplitude
%
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%
% Output:
%   sqi:            resulting sSQI for segment
%
% --
% fecgsyn toolbox, version 1.2, March 2017
% Released under the GNU General Public License
%
% Copyright (C) 2017  Joachim Behar & Fernando Andreotti
% Department of Engineering Science, University of Oxford
% joachim.behar@oxfordalumni.org, fernando.andreotti@eng.ox.ac.uk
%
%
% For more information visit: http://www.fecgsyn.com
%
% Referencing this work
% (SQI indices)
% Andreotti, F., Gräßer, F., Malberg, H., & Zaunseder, S. (2017). Non-Invasive Fetal ECG Signal Quality Assessment for
% Multichannel Heart Rate Estimation. IEEE Trans. Biomed. Eng., (in press).
%
% (FECGSYN Toolbox)
% Behar, J., Andreotti, F., Zaunseder, S., Li, Q., Oster, J., & Clifford, G. D. (2014). An ECG Model for Simulating
% Maternal-Foetal Activity Mixtures on Abdominal ECG Recordings. Physiol. Meas., 35(8), 1537–1550.
%
%
% Last updated : 15-03-2017
%
% 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/>.

y = sort(signal,'descend');
nfac = median(y(1:round(0.05*length(y))));
sqi=std(signal)/nfac;
sqi(sqi>1) = 1; % limit to one, in the exceptional case it is larger than unit

end

author: fernandoandreotti


text:
function obj = removeContact(obj, contact)
    % Remove contact constraints defined for the system
    %
    % Parameters:
    % contact: the name of the contact @type cellstr

    % validate the input argument
    if ischar(contact)
        contact_name = {contact};
    elseif isa(contact, 'ContactFrame')
        contact_name = {contact.Name};
    elseif iscell(contact)
        contact_name = cell(numel(contact),1);
        for i=1:numel(contact)
            c = contact{i};
            if ischar(c)
                contact_name{i} = c;
            elseif isa(c, 'ContactFrame')
                contact_name{i} = c.Name;
            else
                error('You must provide the name of the contact or the ContactFrame object.')
            end
        end
    else
        error('You must provide the name of the contact or the ContactFrame object.')
    end

    for i=1:length(contact_name)
        name = contact_name{i};

        % remove the associated holonomic constraints
        obj.removeHolonomicConstraint(name);

        % remove the associated unilateral constraints
        % 1. friction cone
        fc_cstr = ['fc', name];
        if isfield(obj.UnilateralConstraints, fc_cstr)
            obj.removeUnilateralConstraint(fc_cstr);
        end
        % 2. ZMP constraints
        zmp_cstr = ['zmp', name];
        if isfield(obj.UnilateralConstraints, zmp_cstr)
            obj.removeUnilateralConstraint(zmp_cstr);
        end
    end
end

author: ayonga

repo: frost-dev

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


text:
 function [xs, info] = pl_iot(x, Ab, data, R, varargin)
%function [xs, info] = pl_iot(x, Ab, data, R, [options])
%|
%| Generic penalized-likelihood minimization,
%| for arbitrary negative log-likelihood with convex non-quadratic penalty,
%| via incremental optimization transfer using separable quadratic surrogates.
%|
%| cost(x) = sum_i h(data_i; [Ax]_i) + R(x),
%|
%| in
%|	x	[np 1]		initial estimate
%|	Ab	[nd np]		system matrix (see Gblock)
%|	data	{cell}		whatever data is needed for the likelihood
%|				data{1} must have size "[nb,na]", so for 3D
%|				one must use reshaper(yi, '2d').
%|	R			penalty object (see Robject.m)
%|
%| options
%|	dercurv	{function_handle} function returning derivatives and curvatures
%|				of negative log-likeihood via the call:
%|		[deriv curv] = dercurv(data, Ab{m}*x, curvtype, iblock, nblock)
%|				or choose 'trl' or 'wls' (default)
%|	niter	#		# total iterations (default: 1+0, pure OS!)
%|	os	#		how many "warmup" OS iterations (default: 1)
%|	riter	#		# of penalty subiterations (default: 3).
%|	curvtype ''		curvature type:
%|					'pc' : precomputed (default)
%|					'oc' : optimal (ensures monotonic)
%|	pixmax	(value)		upper bound on pixel values (default: infinity)
%|	pixmin	(value)		lower bound on pixel values (default: 0)
%|	isave	[]		list of iterations to archive
%|				(default: [] 'last)
%|	gi	[nd]		precomputed Ab * 1 factors.
%|	chat	(value)		verbosity (default: 0)
%|	userfun			user defined function handle (see default below)
%|
%| out
%|	xs	[np nsave]	estimates each (saved) iteration
%|	info	[niter+1 ?]	userfun output.  default is cpu times
%|
%| Copyright 2005-3-8, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test')
	run_mfile_local('pwls_example')
return
end
if nargin < 4, ir_usage, end

cpu etic

Ab = block_ob(Ab, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Ab, 'n');
starts = subset_start(nblock);

% defaults
arg.dercurv = 'wls';
arg.curvtype = 'pc';
arg.riter = 3;
arg.os = 1;
arg.niter = arg.os + 0; % default is pure OS, no IOT
arg.isave = [];
arg.pixmin = 0;
arg.pixmax = inf;
arg.gi = [];
arg.update_even_if_denom_0 = true;
arg.chat = false;
arg.userfun = @userfun_default;

% options
arg = vararg_pair(arg, varargin);
arg.isave = iter_saver(arg.isave, arg.niter);

if streq(arg.dercurv, 'wls')
	arg.dercurv = @wls_dercurv;
elseif streq(arg.dercurv, 'trl')
	arg.dercurv = @trl_dercurv;
end

if ~isa(arg.dercurv, 'function_handle')
	error 'dercurv not function handle?'
end

ticker(mfilename, 1, arg.niter)

[nb na] = size(data{1});

% g_i = sum_j g_ij.  caution: this requires real g_ij and g_ij >= 0
if isempty(arg.gi)
	arg.gi = sum(Ab');
	arg.gi = max(reale(arg.gi), 0); % trick: just in case...
	arg.gi = reshape(arg.gi, [nb na]);
end

% precompute likelihood-term denominator if needed
if streq(arg.curvtype, 'pc')
	[dummy curvi] = feval(arg.dercurv, data, 0, 'pc');
	ldenom = Ab' * col(arg.gi .* curvi); % one denominator shared by all subsets
else
	ldenoms = zeros(np, nblock);
end

np = length(x);
x = min(x, arg.pixmax);
x = max(x, arg.pixmin);
xs = zeros(np, length(arg.isave));
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = x;
end

%info = zeros(arg.niter, ?); % trick: do not initialize since size may change

%
% precompute gradient-related state vectors, usually by OS-SPS warmup iterations
%
vvm = zeros(np,nblock);

if arg.os == 0
	for iblock=1:nblock
%		vvm(:,iblock) = Ab{ii}
		error 'not done'
	end
end

%
% SPS-OS iterations to "warm up" (often 1 or 2 iterations suffices)
%
for iter = 1:arg.os
	for iset = 1:nblock
		ticker([mfilename ': os'], [iter iset], [arg.os nblock])

		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;	% l=A*x "line integrals"
		li = reshape(li, nb, length(ia));
		[dhi curvi] = feval(arg.dercurv, data, li, ...
			arg.curvtype, iblock, nblock);

		if ~streq(arg.curvtype, 'pc') % on-the-fly curvatures
			ldenoms(:,iblock) = ...
				Ab{iblock}' * col(arg.gi(:,ia) .* curvi);
			ldenom = nblock * ldenoms(:,iblock);
		end

		lnum = ldenom .* x - nblock * (Ab{iblock}' * dhi(:));

		if iter == arg.os % save last gradient-related state vectors
			vvm(:,iblock) = lnum / nblock;
		end

		x = inner_update(x, lnum, ldenom, R, arg);
	end
	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
	info(1+iter,:) = feval(arg.userfun);
	if arg.chat, printf('Range %g %g', min(x), max(x)), end
end

%
% At this point we have initialized vv_m and ldenoms
% and initialized x based on the last subset.
% It is often logical to update x using all subsets now.
% This is "almost free" since no new likelihood gradients are used.
%
vv = sum(vvm,2);
if streq(arg.curvtype, 'pc')
	ldenoms = ldenom / nblock;
else
	ldenom = sum(ldenoms,2);
end
if 1
	x = inner_update(x, vv, ldenom, R, arg);
	if any(arg.isave == -1)
		xs(:, arg.isave == -1) = x;
	end
end

%
% IOT iterations
%
for iter = (1+arg.os):arg.niter
	for iset = 1:nblock
		ticker([mfilename ': iot'], [iter iset], [arg.niter nblock])

		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;	% l=A*x "line integrals"
		li = reshape(li, nb, length(ia));
		[dhi curvi] = feval(arg.dercurv, data, li, ...
			arg.curvtype, iblock, nblock);

		if ~streq(arg.curvtype, 'pc')
			ldenom = ldenom - ldenoms(:,iblock);
			ldenoms(:,iblock) = ...
				Ab{iblock}' * col(arg.gi(:,ia) .* curvi);
			ldenom = ldenom + ldenoms(:,iblock);
		end

		vv = vv - vvm(:,iblock);
		if streq(arg.curvtype, 'pc')
			vvm(:,iblock) = ldenoms .* x - (Ab{iblock}' * dhi(:));
		else
			vvm(:,iblock) = ldenoms(:,iblock) .* x - (Ab{iblock}' * dhi(:));
		end
		vv = vv + vvm(:,iblock);

		x = inner_update(x, vv, ldenom, R, arg);
	end

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
	info(1+iter,:) = feval(arg.userfun);
	if arg.chat, printf('Range %g %g', min(x), max(x)), end
end

%
% inner_update()
% given SPS for likelihood, update x using regularization subiterations
%
function x = inner_update(x, lnum, ldenom, R, arg)

for ii=1:arg.riter
	if isempty(R)
		num = lnum;
		rdenom = 0;
	else
		rdenom = R.denom(R, x);
		num = lnum - R.cgrad(R, x) + rdenom .* x;
	end

	% update
	if arg.update_even_if_denom_0
		x = num ./ (ldenom + rdenom);
	else
		old = x(ldenom == 0);
		x = num ./ (ldenom + rdenom);
		x(ldenom == 0) = old;
	end
	x = max(x, arg.pixmin);	% lower bound
	x = min(x, arg.pixmax);	% upper bound
end

% default user function.
% using this evalin('caller', ...) trick, one can compute anything of interest
function out = userfun_default
%x = evalin('caller', 'x');
out = [cpu('etoc')];

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% Filename and Format string
% This code loads in the contents of the data file with textscan a block at a
% time and converts the cellarray from textscan into a double array

filename='waferdata.csv';
FormatString=['%*f%*f%*f%*f' repmat('%f',1,9) repmat('%*f',1,7)];

%% Parameters
Headers=3;
NumLines=1e6; % Total number of lines to read
BlockLines=10000; % Size of block
NumBlocks=NumLines/BlockLines; % Number of blocks

%% Open file
fid = fopen(filename); % Open file

%% Preallocate data
data=zeros(NumLines,9); % Pre-allocate space for data

%% Get Headers
cellchunk=textscan(fid,'%s',3,'delimiter','\n'); % Get first 3 lines

%% Read in blocks
for Block=1:NumBlocks
    data((Block-1)*BlockLines+(1:BlockLines),:) = ...
	cell2mat(textscan(fid,FormatString,BlockLines,'delimiter',','));
    disp(['Block ' num2str(Block) ]); % Display current block to show progress
end

%% Close file
fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
%ROTZ Rotation about Z axis
%
% R = ROTZ(THETA) is an SO(3) rotation matrix (3x3) representing a rotation of THETA
% radians about the z-axis.
%
% R = ROTZ(THETA, 'deg') as above but THETA is in degrees.
%
% See also ROTX, ROTY, ANGVEC2R, ROT2.

% 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

function R = rotz(t, deg)
    if nargin > 1 && strcmp(deg, 'deg')
        t = t *pi/180;
    end

    ct = cos(t);
    st = sin(t);
    R = [
        ct  -st  0
        st   ct  0
        0    0   1
        ];

author: Allopart

repo: rbpf-gmapping


text:
function x = p19_start ( i, x )

%*****************************************************************************80
%
%% P19_START returns a starting point for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the starting point.
%
%    Output, real X, the starting point.
%
  if ( i == 1 )
    x = 0.0;
  elseif ( i == 2 )
    x = 1.0;
  elseif ( i == 3 )
    x = 0.5;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P19_START - Fatal error!\n' );
    fprintf ( 1, '  Illegal starting point index = %d\n', i );
    error ( 'P19_START - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = log(f)
%LOG   Natural logarithm of a CHEBFUN3.
%   LOG(F) is the natural logarithm of F. This function gives an error
%   message if the function passes through or becomes numerically close to
%   zero.
%
% See also CHEBFUN3/COMPOSE.

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

% Empty check:
if ( isempty(f) )
    return
end

if ( isreal(f) )
    % Positive/negative test.
    ss = singleSignTest(f);  % Returns TRUE if there is no sign change.
    if ( ~ss )
        error('CHEBFUN:CHEBFUN3:log:notSmooth', ...
            'Sign change detected. Unable to represent the result.');
    end
end

f = compose(f, @log);

end


text:
%% Constrained Total Variation Projection with ADDM
% Test for ADMM algorithm convergence on a randomized example.

addpath('../');
addpath('../toolbox/');

%%
% We want to solve a problem of the form

%%
%   |min_{x in C, norm(K*x,1)<=1} -<x,w> + epsilon/2*norm(x)^2

%%
% Here |K| is some sort of "gradient" operator (here we use a random
% matrix) and C

%%
% This can be written as

%%
%   |min_x F(K*x) + G(x)|

%%
% where |F(u) = i_{norm(u,1)<=1}| and
%  |G(x) = i_{C}(x) - <f,w> + epsilon/2*norm(f)^2|.

%%
% Inner product shortcut.

dotp = @(u,v)sum(u(:).*v(:));

%%
% Dimension of the problem.

n = 300;

%%
% Number of computed "gradient".

p = 500;

%%
% Number of affine constraint

r = 10;

%%
% Regularization.

epsilon = .1;

%%
% Gradient operator.

K = randn(p,n);

%%
% Linear function to optimize.

w = randn(n,1);

%%
% Constraint operator.

A = randn(r,n);
y = randn(r,1)*0;

%%
% Projector on A*f=y.

pA = A'*(A*A')^(-1);
ProjC = @(f)f + pA*(y-A*f);

%%
% Projection on L1 ball.

ProxF = @(u,rho)perform_l1ball_projection(u,1);
ProxFS = compute_dual_prox(ProxF);

%%
% Proximal operator of G.

ProxG = @(x,tau)ProjC( (x+tau*w)/(1+tau*epsilon) );

%%
% Callback to record information during the iterations.

F = @(x)-dotp(x,w) + epsilon/2*norm(x(:))^2;
Constr = @(x)norm(K*x,1);
options.report = @(x)struct('F', F(x), 'Constr', Constr(x));

%%
% Run the algorihtm.

options.niter = 5000;
[f,R] = perform_admm(zeros(n,1), K, K', ProxFS, ProxG, options);

%%
% Retrieve the F and constraint function values.

f = s2v(R,'F');
constr = s2v(R,'Constr');

%%
% Display.

clf;
subplot(2,1,1);
plot(f);
title('Energy');
subplot(2,1,2);
plot(log10(abs(constr-1)));
title('Constraint');

author: gpeyre

repo: matlab-toolboxes


text:
function VAR = MoveCluster(VAR,i_fish,sourceIDs,destIDs)
nClus = length(VAR(i_fish).ClusGroup{sourceIDs{1}});
clus_list = intersect(1:nClus,sourceIDs{2});

Source = VAR(i_fish).ClusGroup{sourceIDs{1}}(clus_list);

dest_clus_list = destIDs{2}(1:length(clus_list));
VAR(i_fish).ClusGroup{destIDs{1}}(dest_clus_list) = Source;

VAR(i_fish).ClusGroup{sourceIDs{1}}(clus_list) = [];
end

author: xiuyechen

repo: FishExplorer


text:
function [MSE, Entropy, SNRdB, sdV] = QuantSNR (Yq, Xq, FPDF, ScaleF)
% Calculate the SNR (dB) for a quantizer defined by a table for a
% signal with a given probability density function.
%
% Yq - Nlev quantizer output levels (in ascending order)
% Xq - Nlev-1 quantizer decision levels (in ascending order)
% FPDF - Cell array of function handles {Farea, Fmean, Fvar}
% Farea - Pointer to a function to calculate the integral of a probability
%   density function in an interval.
%                 b
%   Farea(a,b) = Int p(x) dx,
%                 a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% Fmean - Pointer to a function to calculate the mean of a probability
%   density function in an interval.
%                 b
%   Fmean(a,b) = Int x p(x) dx,
%                 a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% Fvar - Pointer to a function to calculate the second moment of a
%   probability density function in an interval.
%                b
%   Fvar(a,b) = Int x^2 p(x) dx,
%                a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% ScaleF - Scaling factor for the signal (optional, default 1)

% MSE - Resulting mean square error
% Entropy - Resulting entropy (bits)
% SNRdB - Resulting signal-to-noise ratio dB
% sdV - Peak to average ratio for the quantizer. The peak level (overload)
%   V of the quantizer is defined here as the the distance from the mean to
%   the largest output level plus a half step size. Then
%     sdV = sd / V,
%   where sd is the standard deviation of the signal.

if (nargin <= 3)
  ScaleF = 1;
end

% The scale factor represents a scaling of the signal
% For convenience we scale the quantizer by 1/ScaleF and later multipy
% the computed MSE by ScaleF^2 achieve the same effect as scaling the
% signal

% Quantizers are really a combination of two operations
%  - determine an interval k, such that Xq(k-1) <= x < Xq(k)
%  - output a value Yq(k)
% Consider a scaled signal g*x, then k can be determined as the value
% which satisfies
%   Xq(k-1)/g <= x < Xq(k)/g
% The quantization error for the scaled signal is
%   e = g*x - Yq(k)
%     = g * (x - Yq(k)/g)
% The bracketed term is the error for a signal x (not scaled) applied
% to a quantizer with break points Xq(.)/g and output levels Yq(.)/g.

Nlev = length(Yq);
if (ScaleF ~=1)
  if (Nlev > 1)
    Xq = Xq / ScaleF;
  end
  Yq = Yq / ScaleF;
end

% Calculate the mean square error
MSE = ScaleF^2 * QuantMSE (Yq, Xq, FPDF);

% Calculate the entropy
if (nargout >= 2)
  Entropy = QuantEntropy (Xq, FPDF{1});
end

% Calculate the SNR in dB
if (nargout >= 3)
  Fmean = FPDF{2};
  Fvar = FPDF{3};
  Xmean = feval(Fmean, -Inf, Inf);
  Svar = feval(Fvar, -Inf, Inf) - Xmean^2;
  SNRdB = 10 * log10(ScaleF^2 * Svar / MSE);
end

% Calculate the peak to average ratio for the quantizer
if (nargout >= 4)
  if (Nlev > 1)
    VL = (Yq(1)-Xq(1)) + Yq(1);
    VH = (Yq(end)-Xq(end)) + Yq(end);
    V = max(abs(VL-Xmean), abs(VH-Xmean));
  else
    V = Inf;
  end
  sdV = sqrt(Svar) / V;
end

return


text:
function D = driving_function_mono_wfs_fs(x0,nx0,xs,f,conf)
%DRIVING_FUNCTION_MONO_WFS_FS driving signal for a focused source in WFS
%
%   Usage: D = driving_function_mono_wfs_fs(x0,nx0,xs,f,[conf])
%
%   Input parameters:
%       x0          - position of the secondary sources / m [nx3]
%       nx0         - directions of the secondary sources / m [nx3]
%       xs          - position of virtual focused source / m [nx3]
%       f           - frequency of the monochromatic source / Hz
%       conf        - optional configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   DRIVING_FUNCTION_MONO_WFS_FS(x0,xs,f,src,conf) returns WFS driving signals
%   for the given secondary sources, the virtual focused source position and the
%   frequency f. For 3D and 2.5D the default behavior is to use a focused point
%   source as source model, for 2D a focused line source is used instead.
%
%   See also: driving_function_mono_wfs, driving_function_imp_wfs_ps
%
%   References:
%       Spors, Wierstorf, Geier, Ahrens (2009) - "Physical and Perceptual
%       Properties of Focused Sources in Wave Field Synthesis", 127th Convention
%       of the Audio Engineering Society, Paper 7914,
%       http://www.aes.org/e-lib/browse.cfm?elib=15109
%
%       Verheijen (1997) - "Sound Reproduction by Wave Field Synthesis", PhD
%       thesis, TU Delft,
%       http://resolver.tudelft.nl/uuid:9a35b281-f19d-4f08-bec7-64f6920a3821
%
%       Wierstorf (2014) - "Perceptual Assessment of Sound Field Synthesis",
%       PhD thesis, TU Berlin, https://doi.org/10.14279/depositonce-4310

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargmatrix(x0,nx0,xs);
isargpositivescalar(f);
isargstruct(conf);
if size(xs,2)~=3
    error('%s: size(xs,2) has to be 3 and not %i.',upper(mfilename),size(xs,2));
end

%% ===== Configuration ===================================================
xref = conf.xref;
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation =====================================================
% Calculate the driving function in time-frequency domain

% Frequency
omega = 2*pi*f;

if strcmp('2D',dimension) || strcmp('3D',dimension)

    % === 2- or 3-Dimensional ============================================

    % For 2D the default focussed source should be a line sink
    if strcmp('2D',dimension) && strcmp('default',driving_functions)
        driving_functions = 'line_sink';
    end

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % D using a point source and the approximation w/c|x0-xs|>>1
        %
        %            1  i w (x0-xs) nx0
        % D(x0,w) = --- --- ----------- e^(i w/c |x0-xs|)
        %           2pi  c   |x0-xs|^2
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-fd-wfs-focused-3d
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = 1./(2.*pi) .* (1i.*omega)./c ...
            .* vector_product(x0-xs,nx0,2) ./ r.^2 ...
            .* exp(+1i.*omega./c.*r);
        %
    case 'point_sink'
        % D using a point sink as source model
        %
        % D(x0,w) =
        %
        %  1  /i w      1    \  (x0-xs) nx0
        % --- |--- + ------- |  ----------- e^(i w/c |x0-xs|)
        % 2pi \ c    |x0-xs| /   |x0-xs|^2
        %
        % See Wierstorf (2014), eq.(2.71)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = 1./(2.*pi) .* ( -1i.*omega./c + 1./r ) ...
            .* vector_product(x0-xs,nx0,2) ./ r.^2 ...
            .* exp(+1i.*omega./c.*r);
        %
    case 'line_sink'
        % D using a line sink as source model
        %
        %              iw (x0-xs)nk  (1)/ w         \
        % D(x0,w) =  - -- --------- H1  | - |x0-xs| |
        %              2c  |x0-xs|      \ c         /
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = -1i.*omega./(2.*c) ...
            .* vector_product(x0-xs,nx0,2) ./ r ...
            .* besselh(1,1,omega./c.*r);
        %
    case 'legacy'
        % --- Old SFS Toolbox default ------------------------------------
        % D using an approximated point sink as source model
        %
        %            1  i w (x0-xs) nx0
        % D(x0,w) = --- --- ------------- e^(i w/c |x0-xs|)
        %           2pi  c  |x0-xs|^(3/2)
        %
        % See Wierstorf (2014), eq.(2.73)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = 1./(2.*pi) .* 1i.*omega./c ...
            .* vector_product(x0-xs,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i.*omega./c.*r);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a focused source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    % Reference point
    xref = repmat(xref,[size(x0,1) 1]);

    switch driving_functions
    case {'default', 'reference_circle'}
        % Driving function with two stationary phase approximations,
        % reference to circle around the focused source with radius |xref-xs|
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        %
        % 2.5D correction factor
        %         _____________
        %        |        r
        % g0 = _ |1 + ---------
        %       \|    |xref-xs|
        %
        g0 = sqrt( 1 + r./vector_norm(xref-xs,2) );
        %                       ___     ___
        %                      | 1     |-iw  (xs-x0) nx0
        % D_2.5D(x0,w) = g0  _ |---  _ |--- ------------- e^(i w/c |x0-xs|)
        %                     \|2pi   \| c  |x0-xs|^(3/2)
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-fd-wfs-focused-25d
        %
        % Driving signal
        D = 1./sqrt(2.*pi) .* sqrt(-1i.*omega./c) .* g0 ...
            .* vector_product(xs-x0,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i.*omega./c.*r);
        %
    case 'reference_point'
        % Driving function with only one stationary phase approximation,
        % reference to one point in field
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % 2.5D correction factor
        %         _____________________
        %        |      |xref-x0|
        % g0 = _ |---------------------
        %       \|||xref-x0| - |xs-x0||
        %
        % Verheijen (1997), eq. (A.14)
        %
        g0 = sqrt( vector_norm(xref-x0,2) ./ abs(vector_norm(x0-xref,2) - r) );
        %                       ___     ___
        %                      | 1     |-iw  (xs-x0) nx0
        % D_2.5D(x0,w) = g0  _ |---  _ |--- ------------- e^(i w/c |x0-xs|)
        %                     \|2pi   \| c  |x0-xs|^(3/2)
        %
        % Driving signal
        D = 1./sqrt(2.*pi) .* sqrt(-1i.*omega./c) .* g0 ...
            .* vector_product(xs-x0,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i.*omega./c.*r);
        %
    case {'reference_line', 'delft1988'}
        % Driving function with two stationary phase approximations,
        % reference to a line parallel to a LINEAR secondary source distribution
        %
        % Distance ref-line to linear ssd
        dref = abs( vector_product(xref-x0,nx0,2) );
        % Distance source and linear ssd
        ds = abs( vector_product(xs-x0,nx0,2) );
        %
        % 2.5D correction factor
        %        ______________________
        % g0 = \| d_ref / (d_ref - d_s)
        %
        % See Start (1997), eq. (3.16)
        %
        g0 = sqrt( dref ./ (dref - ds) );
        %                       ___     ___
        %                      | 1     |-iw  (xs-x0) nx0
        % D_2.5D(x0,w) = g0  _ |---  _ |--- ------------- e^(i w/c |x0-xs|)
        %                     \|2pi   \| c  |x0-xs|^(3/2)
        %
        % See Verheijen (1997), eq. (2.29b)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = 1./sqrt(2.*pi) .* sqrt(-1i.*omega./c) .* g0 ...
            .* vector_product(xs-x0,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i.*omega./c.*r);
        %
    case 'legacy'
        % --- Old SFS Toolbox default ------------------------------------
        % 2.5D correction factor
        %        ______________
        % g0 = \| 2pi |xref-x0|
        %
        g0 = sqrt(2*pi*vector_norm(xref-x0,2));
        %
        % D_2.5D using an approximated point sink as source model
        %                        ___
        %                 g0    |iw   (x0-xs) nx0
        % D_2.5D(x0,w) = ---  _ |--- ------------- e^(i w/c |x0-xs|)
        %                2pi   \| c  |x0-xs|^(3/2)
        %
        % See Wierstorf (2014), eq.(2.74)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = 1./(2.*pi) .* sqrt(1i.*omega./c) .* g0 ...
            .* vector_product(x0-xs,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i*omega./c.*r);
        %
    case 'spors2009eq7'
        % --- SFS Toolbox ------------------------------------------------
        % 2.5D correction factor
        %        ______________
        % g0 = \| 2pi |xref-x0|
        %
        g0 = sqrt(2*pi*vector_norm(xref-x0,2));
        %
        % D_2.5D using a line sink with point amplitude characteristic as source
        % and the large argument approximation of the driving function above.
        % This results in the "traditional" driving function, derived in
        % Verheijen (1997), see Spors (2009) eq. 7.
        %                       ______
        %                      |   w  |  (x0-xs) nx0
        % D_2.5D(x0,w) = -g0 _ |------  ------------- e^(i w/c |x0-xs|)
        %                     \|2pi ic  |x0-xs|^(3/2)
        %
        % See Spors et al. (2009), eq.(7)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = -g0 .* sqrt( omega./(2.*pi.*c.*1i) ) ...
            .* vector_product(x0-xs,nx0,2) ./ r.^(3./2) ...
            .* exp(+1i*omega./c.*r);
        %
    case 'spors2009eq6'
        % --- Spors 2009 -------------------------------------------------
        % 2.5D correction factor
        %        ______________
        % g0 = \| 2pi |xref-x0|
        %
        g0 = sqrt(2*pi*vector_norm(xref-x0,2));
        %
        % D_2.5D using a line sink with point source amplitude characteristics as
        % source
        %
        %                    iw (x0-xs) nx0   (1)/ w         \
        % D_2.5D(x0,w) = -g0 -- -----------  H1  | - |x0-xs| |
        %                    2c   |x0-xs|        \ c         /
        %
        % See Spors et al. (2009), eq.(6)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = -g0 .* 1i.*omega./(2.*c) ...
            .* vector_product(x0-xs,nx0,2) ./ r ...
            .* besselh(1,1,omega./c.*r);
        %
    case 'point_sink'
        % --- Point Sink -------------------------------------------------
        % 2.5D correction factor
        %        ______________
        % g0 = \| 2pi |xref-x0|
        %
        g0 = sqrt(2*pi*vector_norm(xref-x0,2));
        %
        % D_2.5D using a point sink as source model
        %
        % D_2.5D(x0,w) =
        %             ___       ___
        %    g0  /   |i w|     | c |    1    \  (x0-xs) nx0
        %   ---  | _ |---  + _ |---  ------- |  ----------- e^(i w/c |x0-xs|)
        %   2pi  \  \| c      \|i w  |x0-xs| /   |x0-xs|^2
        %
        % See ...
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Driving signal
        D = g0./(2.*pi) .* (sqrt(1i.*omega./c) + sqrt(c./(1i.*omega)) ./ r) ...
            .* vector_product(x0-xs,nx0,2) ./ r.^2 ...
            .* exp(+1i*omega./c.*r);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D focused source.'],upper(mfilename),driving_functions);
    end
else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function linpack_c_test10 ( )

%*****************************************************************************80
%
%% TEST10 tests CGTSL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  For a complex tridiagonal matrix:\n' );
  fprintf ( 1, '  CGTSL solves a linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  seed = 123456789;

  c(1) = 0.0;
  [ c(2:n), seed ] = c4vec_uniform_01 ( n-1, seed );

  [ e(1:n-1), seed ] = c4vec_uniform_01 ( n-1, seed );
  e(n) = 0.0;

  d(1:n) = 0.0;

  d(1:n-1) = d(1:n-1) - 2.0 * e(1:n-1);
  d(2:n)   = d(2:n) - 2.0 * c(2:n);
%
%  Set the desired solution.
%
  for i = 1 : n
    x(i) = complex ( i, 10 * i );
  end
%
%  Compute the corresponding right hand side.
%
  b(1) = d(1) * x(1) + e(1) * x(2);
  for i = 2 : n-1
    b(i) = c(i) * x(i-1) + d(i) * x(i) + e(i) * x(i+1);
  end
  b(n) = c(n) * x(n-1) + d(n) * x(n);
%
%  Solve the tridiagonal system.
%
  [ b, info ] = cgtsl ( n, c, d, e, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed                     Exact\n' );
  fprintf ( 1, '  Solution                     Solution\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  (%8f  %8f)  (%8f  %8f)\n', ...
      real ( b(i) ), imag ( b(i) ), real ( x(i) ), imag ( x(i) ) );
  end

  return
end


text:
function t_data = parameterize_arc_length ( m, data_num, p_data )

%*****************************************************************************80
%
%% PARAMETERIZE_ARC_LENGTH parameterizes data by pseudo-arclength.
%
%  Discussion:
%
%    A parameterization is required for the interpolation.
%
%    This routine provides a parameterization by computing the
%    pseudo-arclength of the data, that is, the Euclidean distance
%    between successive points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Input, real P_DATA(M,DATA_NUM), the data values.
%
%    Output, real T_DATA(DATA_NUM), parameter values
%    assigned to the data.
%
  t_data(1) = 0.0;
  for j = 2 : data_num
    t_data(j) = t_data(j-1) ...
      + sqrt ( sum ( ( p_data(1:m,j) - p_data(1:m,j-1) ).^2 ) );
  end

  return
end


text:
%% Copyright (C) 2016-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defun sinhint (@var{x})
%% Numerical sinhint function.
%%
%% Example:
%% @example
%% @group
%% sinhint (1.1)
%%   @result{} ans = 1.1767
%% @end group
%% @end example
%%
%% @strong{Note} this function may be slow for large numbers of inputs.
%% This is because it is not a native double-precision implementation
%% but rather the numerical evaluation of the Python @code{mpmath} function
%% @code{shi}.
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @seealso{@@sym/sinhint}
%% @end defun

function y = sinhint (x)
  if (nargin ~= 1)
    print_usage ();
  end
  cmd = { 'L = _ins[0]'
          'A = [complex(mpmath.shi(x)) for x in L]'
          'return A,' };
  c = pycall_sympy__ (cmd, num2cell (x(:)));
  y = reshape (cell2mat (c), size (x));
end

%!error sinhint (1, 2)

%!test
%! x = 1.1;
%! y = sym(11)/10;
%! A = sinhint (x);
%! B = double (sinhint (y));
%! assert (A, B, -4*eps);

%!test
%! y = [2 3 sym(pi); exp(sym(1)) 5 6];
%! x = double (y);
%! A = sinhint (x);
%! B = double (sinhint (y));
%! assert (A, B, -4*eps);

%!test
%! % maple:
%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
%! % > for a in A do evalf(Shi(a)) end do;
%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
%! A = [  0.45769171128668800567 + 1.8332099215048436239*1i
%!        0.60888490416819746440 + 1.2482232175376056201*1i
%!        0.13577763724269399110e42
%!        1.6583475942188740493*1i
%!       -0.00010000000005553888891 + 0.10000000016666111119e-5*1i
%!       -0.74701205140887967022e7 - 0.10381447401236722090e8*1i ];
%! B = sinhint (x);
%! assert (A, B, -eps)

%!assert (sinhint (inf), inf)
%!assert (sinhint (-inf), -inf)

repo: octsympy


text:
function ft_write_data(filename, dat, varargin)

% FT_WRITE_DATA exports electrophysiological data such as EEG to a file.
%
% Use as
%   ft_write_data(filename, dat, ...)
%
% The specified filename can contain the filename extension. If it has no filename
% extension not, it will be added automatically.
%
% Additional options should be specified in key-value pairs and can be
%   'header'       = header structure that describes the data, see FT_READ_HEADER
%   'event'        = event structure that corresponds to the data, see FT_READ_EVENT
%   'chanindx'     = 1xN array, for selecting a subset of channels from header and data
%   'dataformat'   = string, see below
%   'append'       = boolean, not supported for all formats
%
% The supported dataformats for writing are
%   edf
%   gdf
%   anywave_ades
%   brainvision_eeg
%   neuralynx_ncs
%   neuralynx_sdma
%   plexon_nex
%   fcdc_matbin
%   fcdc_mysql
%   fcdc_buffer
%   flac, m4a, mp4, oga, ogg, wav (audio formats)
%   matlab
%   homer_nirs
%   snirf
%
% For EEG data, the input data is assumed to be scaled in microvolt.
% For NIRS data, the input data is assumed to represent optical densities.
%
% See also FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, FT_WRITE_EVENT

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

global data_queue    % for fcdc_global
global header_queue  % for fcdc_global
global db_blob       % for fcdc_mysql
if isempty(db_blob)
  db_blob = 0;
end

% get the options
append        = ft_getopt(varargin, 'append', false);
nbits         = ft_getopt(varargin, 'nbits', 16); % for audio
chanindx      = ft_getopt(varargin, 'chanindx');
hdr           = ft_getopt(varargin, 'header');
evt           = ft_getopt(varargin, 'event');
dataformat    = ft_getopt(varargin, 'dataformat');

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

if startsWith(dataformat, 'audio_')
  % support for  audio formats is implemented in a generic fashion
  dataformat = dataformat(7:end);
end

% convert 'yes' or 'no' string into boolean
append = istrue(append);

% ensure that the directory exists if we want to write to a file
if ~ismember(dataformat, {'empty', 'fcdc_global', 'fcdc_buffer', 'fcdc_mysql'})
  isdir_or_mkdir(fileparts(filename));
end

% determine the data size
[nchans, nsamples] = size(dat);

% ensure that the header is (reasonably) complete
if ~isfield(hdr, 'nChans')
  if isfield(hdr, 'label')
    hdr.nChans = length(hdr.label);
  else
    hdr.nChans = nchans;
  end
end

if ~isfield(hdr, 'label')
  hdr.label = arrayfun(@num2str, 1:hdr.nChans, 'UniformOutput', false)';
end

if ~isfield(hdr, 'chantype')
  % use a helper function which has some built in intelligence
  hdr.chantype = ft_chantype(hdr);
end

if ~isfield(hdr, 'chanunit')
  % use a helper function which has some built in intelligence
  hdr.chanunit = ft_chanunit(hdr);
end

if ~isempty(chanindx)
  % the header (and possibly the data) correspond to the original multichannel file
  hdr.label    = hdr.label(chanindx);
  hdr.chantype = hdr.chantype(chanindx);
  hdr.chanunit = hdr.chanunit(chanindx);
  hdr.nChans   = length(chanindx);
  if length(chanindx)==nchans
    % assume that the data already represents the desired subset of channels
  else
    % assume that the data corresponds to the original multichannel file
    dat = dat(chanindx,:);
  end
end

switch dataformat

  case 'empty'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % just pretend that we are writing the data, this is only for debugging
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [numC, numS] = size(dat);
    ft_info('Pretending to write %i samples from %i channels...\n',numS,numC);
    % Insert a small delay to make this more realitic for testing purposes
    % The time for writing to an actual location will differ and depend on
    % the amount of data
    pause(0.001);

  case 'fcdc_global'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % store it in a global variable, this is only for debugging
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(evt)
      ft_error('writing events is not supported here, please see FT_WRITE_EVENT');
    end

    if ~isempty(hdr)
      header_queue = hdr;
    end

    if isempty(data_queue) || ~append
      data_queue = dat;
    else
      data_queue = cat(2, data_queue, dat);
    end

  case 'fcdc_buffer'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % write to a network transparent buffer for realtime analysis
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(evt)
      ft_error('writing events is not supported here, please see FT_WRITE_EVENT');
    end

    [host, port] = filetype_check_uri(filename);

    type = {
      'char'
      'uint8'
      'uint16'
      'uint32'
      'uint64'
      'int8'
      'int16'
      'int32'
      'int64'
      'single'
      'double'
      };

    wordsize = {
      1 % 'char'
      1 % 'uint8'
      2 % 'uint16'
      4 % 'uint32'
      8 % 'uint64'
      1 % 'int8'
      2 % 'int16'
      4 % 'int32'
      8 % 'int64'
      4 % 'single'
      8 % 'double'
      };

    % this should only be done the first time
    if ~append && ~isempty(hdr)
      % reformat the header into a buffer-compatible format
      packet.fsample   = hdr.Fs;
      packet.nchans    = hdr.nChans;
      packet.nsamples  = 0;
      packet.nevents   = 0;
      packet.data_type = find(strcmp(type, class(dat))) - 1; % zero-offset
      if isfield(hdr,'label') && iscell(hdr.label)
        packet.channel_names = hdr.label;
      end
      if isfield(hdr,'siemensap')
        if isa(hdr.siemensap, 'uint8')
          packet.siemensap = hdr.siemensap;
        else
          %          try
          %            packet.siemensap = matlab2sap(hdr.siemensap);
          %          catch
          warning 'Ignoring field "siemensap"';
          %          end
        end
      end
      if isfield(hdr,'nifti_1')
        if isa(hdr.nifti_1, 'uint8')
          packet.nifti_1 = hdr.nifti_1;
        else
          try
            packet.nifti_1 = encode_nifti1(hdr.nifti_1);
          catch
            warning 'Ignoring field "nifti_1"';
          end
        end
      end
      if isfield(hdr,'ctf_res4')
        if isa(hdr.ctf_res4, 'uint8')
          packet.ctf_res4 = hdr.ctf_res4;
        else
          warning 'Ignoring non-uint8 field "ctf_res4"';
        end
      end

      % try to put_hdr and initialize if necessary
      try
        % try writing the packet
        buffer('put_hdr', packet, host, port);

      catch
        if contains(lasterr, 'Buffer size N must be an integer-valued scalar double.')
          % this happens if the MATLAB75/toolbox/signal/signal/buffer
          % function is used instead of the FieldTrip buffer
          ft_error('the FieldTrip buffer mex file was not found on your path, it should be in fieldtrip/fileio/private');

        elseif contains(lasterr, 'failed to create socket') && (strcmp(host, 'localhost') || strcmp(host, '127.0.0.1'))

          % start a local instance of the TCP server
          ft_warning('starting FieldTrip buffer on %s:%d', host, port);
          buffer('tcpserver', 'init', host, port);
          pause(1);

          % rewrite the packet until success
          success = false;
          while ~success
            try
              % it may take some time before the TCP server is fully initialized
              % try writing the packet again
              buffer('put_hdr', packet, host, port);
              success = true;
            catch
              success = false;
            end
          end
        end % if strfind...

      end % try

    end % writing header

    if ~isempty(dat)
      MAXNUMSAMPLE = 600000; % see buffer.h
      if size(dat,2)>MAXNUMSAMPLE
        ft_error('number of samples exceeds the size of the ring buffer');
      end
      % reformat the data into a buffer-compatible format
      packet.nchans    = size(dat,1);
      packet.nsamples  = size(dat,2);
      packet.data_type = find(strcmp(type, class(dat))) - 1; % zero-offset
      packet.bufsize   = numel(dat) * wordsize{strcmp(type, class(dat))};
      packet.buf       = dat;
      buffer('put_dat', packet, host, port);
    end

  case {'brainvision_eeg', 'brainvision_vhdr'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % combination of *.eeg and *.vhdr file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    % the header should at least contain the following fields
    %   hdr.label
    %   hdr.nChans
    %   hdr.Fs
    write_brainvision_eeg(filename, hdr, dat, evt);

  case 'fcdc_matbin'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiplexed data in a *.bin file (ieee-le, 64 bit floating point values),
    % accompanied by a MATLAB V6 file containing the header
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [path, file, ext] = fileparts(filename);
    headerfile = fullfile(path, [file '.mat']);
    datafile   = fullfile(path, [file '.bin']);

    if append && exist(headerfile, 'file') && exist(datafile, 'file')
      % read the existing header and perform a sanity check
      old = load(headerfile);
      assert(old.hdr.nChans==size(dat,1));

      % update the existing header
      hdr          = old.hdr;
      hdr.nSamples = hdr.nSamples + nsamples;

      if isfield(old, 'event')
        event = old.event;
      else
        event = [];
      end
      % append the existing and the new events
      event = appendstruct(event, evt);

      save(headerfile, 'hdr', 'event', '-v6');

      % update the data file
      fid = fopen_or_error(datafile,'ab','ieee-le');
      fwrite(fid, dat, hdr.precision);
      fclose(fid);

    else
      hdr.nSamples = nsamples;
      hdr.nTrials  = 1;
      if ~isfield(hdr, 'precision')
        hdr.precision = 'double';
      end

      % rename the variable name for the new events
      event = evt;

      % write the header and events to the file
      save(headerfile, 'hdr', 'event', '-v6');

      % write the data file
      fid = fopen_or_error(datafile,'wb','ieee-le');
      fwrite(fid, dat, hdr.precision);
      fclose(fid);
    end

  case 'fcdc_mysql'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % write to a MySQL server listening somewhere else on the network
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(evt)
      ft_error('writing events is not supported here, please see FT_WRITE_EVENT');
    end

    % check that the required low-level toolbox is available
    ft_hastoolbox('mysql', 1);
    db_open(filename);

    if ~isempty(hdr) && isempty(dat)
      % insert the header information into the database
      if db_blob
        % insert the structure into the database table as a binary blob
        db_insert_blob('fieldtrip.header', 'msg', hdr);
      else
        % make a structure with the same elements as the fields in the database table
        s             = struct;
        s.Fs          = hdr.Fs;           % sampling frequency
        s.nChans      = hdr.nChans;       % number of channels
        s.nSamples    = hdr.nSamples;     % number of samples per trial
        s.nSamplesPre = hdr.nSamplesPre;  % number of pre-trigger samples in each trial
        s.nTrials     = hdr.nTrials;      % number of trials
        s.label       = mxSerialize(hdr.label);
        try
          s.msg = mxSerialize(hdr);
        catch
          ft_warning(lasterr);
        end
        db_insert('fieldtrip.header', s);
      end

    elseif isempty(hdr) && ~isempty(dat)
      dim = size(dat);
      if numel(dim)==2
        % ensure that the data dimensions correspond to ntrials X nchans X samples
        dim = [1 dim];
        dat = reshape(dat, dim);
      end
      ntrials = dim(1);
      for i=1:ntrials
        if db_blob
          % insert the data into the database table as a binary blob
          db_insert_blob('fieldtrip.data', 'msg', reshape(dat(i,:,:), dim(2:end)));
        else
          % create a structure with the same fields as the database table
          s = struct;
          s.nChans   = dim(2);
          s.nSamples = dim(3);
          try
            s.data = mxSerialize(reshape(dat(i,:,:), dim(2:end)));
          catch
            ft_warning(lasterr);
          end
          % insert the structure into the database
          db_insert('fieldtrip.data', s);
        end
      end

    else
      ft_error('you should specify either the header or the data when writing to a MySQL database');
    end

  case 'matlab'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plain MATLAB file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file '.mat']);

    if  append &&  exist(filename, 'file')
      % read the previous header and data from MATLAB file
      prev = load(filename);
      % do a sanity chjeck to ensure that the file content is consistent with the new data
      if ~isempty(hdr) && ~isequal(hdr, prev.hdr)
        ft_error('inconsistent header');
      end

    elseif append && ~exist(filename, 'file')
      % file does not yet exist, which is not a problem
      prev = [];

    elseif ~append &&  exist(filename, 'file')
      % file already exists, delete it and make a new one further down
      ft_warning('deleting existing file ''%s''', filename);
      delete(filename);
      prev = [];

    elseif ~append && ~exist(filename, 'file')
      % file does not yet exist, which is not a problem
      prev = [];
    end

    if isfield(prev, 'dat')
      % append the new data to the previous data from from the MATLAB file
      dat = cat(2, prev.dat, dat);
    end
    if isfield(prev, 'event')
      % append the new events to the previous events from from the MATLAB file
      event = cat(2, prev.event, evt);
    else
      % rename the variable name for the new events
      event = evt;
    end

    % write the data, header and events to the file
    save(filename, 'dat', 'hdr', 'event');

  case 'mff'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % MFF files using Phillips plugin
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ft_hastoolbox('mffmatlabio', 1);
    mff_fileio_write(filename, hdr, dat, evt);

  case 'neuralynx_sdma'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The first version of this file format contained in the first 8 bytes the
    % channel label as string. Subsequently it contained 32 bit integer values.
    %
    % The second version of this file format starts with 8 bytes describing (as
    % a space-padded string) the data type. The channel label is contained in
    % the filename as dataset.chanlabel.bin.
    %
    % The third version of this file format starts with 7 bytes describing (as
    % a zero-padded string) the data type, followed by the 8th byte which
    % describes the downscaling for the 8 and 16 bit integer representations.
    % The downscaling itself is represented as uint8 and should be interpreted as
    % the number of bits to shift. The channel label is contained in the
    % filename as dataset.chanlabel.bin.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end

    statuschannel = {
      'stx'
      'pid'
      'siz'
      'tsh'
      'tsl'
      'cpu'
      'ttl'
      'x01'
      'x02'
      'x03'
      'x04'
      'x05'
      'x06'
      'x07'
      'x08'
      'x09'
      'x10'
      'crc'
      };

    dirname = filename;
    clear filename
    [path, file] = fileparts(dirname);
    for i=1:hdr.nChans
      downscale(i) = 0;
      if ~isempty(strmatch(hdr.label{i}, statuschannel))
        format{i} = 'uint32';
      else
        format{i} = 'int32';
      end
      filename{i} = fullfile(dirname, [file '.' hdr.label{i} '.bin']);
    end

    if ~isfolder(dirname)
      mkdir(dirname);
    end

    % open and write to the output files, one for each selected channel
    fid = zeros(hdr.nChans,1);
    for j=1:hdr.nChans

      if append==false
        fid(j) = fopen_or_error(filename{j}, 'wb', 'ieee-le'); % open the file
        magic = format{j};                               % this used to be the channel name
        magic((end+1):8) = 0;                            % pad with zeros
        magic(8) = downscale(j);                         % number of bits to shift
        fwrite(fid(j), magic(1:8));                      % write the 8-byte file header
      else
        fid(j) = fopen_or_error(filename{j}, 'ab', 'ieee-le');    % open the file for appending
      end % if append

      % convert the data into the correct class
      buf = dat(j,:);
      if ~strcmp(class(buf), format{j})
        switch format{j}
          case 'int16'
            buf = int16(buf);
          case 'int32'
            buf = int32(buf);
          case 'single'
            buf = single(buf);
          case 'double'
            buf = double(buf);
          case 'uint32'
            buf = uint32(buf);
          otherwise
            ft_error('unsupported format conversion');
        end
      end

      % apply the scaling, this corresponds to bit shifting
      buf = buf ./ (2^downscale(j));

      % write the segment of data to the output file
      fwrite(fid(j), buf, format{j}, 'ieee-le');

      fclose(fid(j));
    end % for each channel

  case {'flac' 'm4a' 'mp4' 'oga' 'ogg' 'wav'}
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This writes data Y to a Windows WAVE file specified by the file name
    % WAVEFILE, with a sample rate of FS Hz and with NBITS number of bits.
    % NBITS must be 8, 16, 24, or 32.  For NBITS < 32, amplitude values
    % outside the range [-1,+1] are clipped
    %
    % Supported extensions for AUDIOWRITE are:
    %   .flac
    % 	.m4a
    % 	.mp4
    % 	.oga
    % 	.ogg
    % 	.wav
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file '.' dataformat]);

    switch dataformat
      case {'m4a' 'mp4' 'oga' 'ogg'}
        options = {};
      otherwise
        options = {'BitsPerSample', nbits};
    end % switch

    audiowrite(filename, dat', hdr.Fs, options{:});

  case 'plexon_nex'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % single or mulitple channel Plexon NEX file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end
    if nchans~=1
      ft_error('only supported for single-channel data');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file, '.nex']);

    % construct a NEX structure with  the required parts of the header
    nex.hdr.VarHeader.Type       = 5; % continuous
    nex.hdr.VarHeader.Name       = hdr.label{1};
    nex.hdr.VarHeader.WFrequency = hdr.Fs;
    if isfield(hdr, 'FirstTimeStamp')
      nex.hdr.FileHeader.Frequency = hdr.Fs * hdr.TimeStampPerSample;
      nex.var.ts = hdr.FirstTimeStamp;
    else
      ft_warning('no timestamp information available');
      nex.hdr.FileHeader.Frequency  = nan;
      nex.var.ts = nan;
    end
    nex.var.indx = 0;
    nex.var.dat  = dat;

    write_plexon_nex(filename, nex);

    if 0
      % the following code snippet can be used for testing
      [nex2.var, nex2.hdr] = read_plexon_nex(filename, 'channel', 1);
    end

  case 'neuralynx_ncs'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % single channel Neuralynx NCS file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end
    if nchans>1
      ft_error('only supported for single-channel data');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file, '.ncs']);

    LABEL      = hdr.label{1};  % single channel
    ADCHANNEL  = -1;            % unknown
    FSAMPLE    = hdr.Fs;
    RECORDNSMP = 512;
    RECORDSIZE = 1044;

    % cut the downsampled LFP data into record-size pieces
    nrecords = ceil(nsamples/RECORDNSMP);
    fprintf('construct ncs with %d records\n', nrecords);

    % construct a ncs structure with all header details and the data in it
    ncs                = [];
    ncs.NumValidSamp   = ones(1,nrecords) * RECORDNSMP;   % except for the last block
    ncs.ChanNumber     = ones(1,nrecords) * ADCHANNEL;
    ncs.SampFreq       = ones(1,nrecords) * FSAMPLE;
    ncs.TimeStamp      = zeros(1,nrecords,'uint64');

    if rem(nsamples, RECORDNSMP)>0
      % the data length is not an integer number of records, pad the last record with zeros
      dat = cat(2, dat, zeros(nchans, nrecords*RECORDNSMP-nsamples));
      ncs.NumValidSamp(end) = rem(nsamples, RECORDNSMP);
    end

    ncs.dat = reshape(dat, RECORDNSMP, nrecords);

    for i=1:nrecords
      % timestamps should be 64 bit unsigned integers
      ncs.TimeStamp(i) = uint64(hdr.FirstTimeStamp) + uint64((i-1)*RECORDNSMP*hdr.TimeStampPerSample);
    end

    % add the elements that will go into the ascii header
    ncs.hdr.CheetahRev            = '4.23.0';
    ncs.hdr.NLX_Base_Class_Type   = 'CscAcqEnt';
    ncs.hdr.NLX_Base_Class_Name   = LABEL;
    ncs.hdr.RecordSize            = RECORDSIZE;
    ncs.hdr.ADChannel             = ADCHANNEL;
    ncs.hdr.SamplingFrequency     = FSAMPLE;

    % write it to a file
    fprintf('writing to %s\n', filename);
    write_neuralynx_ncs(filename, ncs);

    if false
      % the following code snippet can be used for testing
      ncs2 = read_neuralynx_ncs(filename, 1, inf);
    end

  case 'gdf'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiple channel GDF file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file, '.gdf']);

    write_gdf(filename, hdr, dat);

  case 'edf'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiple channel European Data Format file
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end

    [path, file, ext] = fileparts(filename);
    filename = fullfile(path, [file, '.edf']);

    write_edf(filename, hdr, dat);

  case 'anywave_ades'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % see http://meg.univ-amu.fr/wiki/AnyWave:ADES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end
    if ~isempty(evt)
      ft_error('writing events is not supported');
    end

    dattype = unique(hdr.chantype);
    datunit = cell(size(dattype));
    for i=1:numel(dattype)
      unit = hdr.chanunit(strcmp(hdr.chantype, dattype{i}));
      if ~all(strcmp(unit, unit{1}))
        ft_error('channels of the same type with different units are not supported');
      end
      datunit{i} = unit{1};
    end

    % only change these after checking channel types and units
    chantype = adestype(hdr.chantype);
    dattype  = adestype(dattype);

    % ensure that all channels have the right scaling
    for i=1:size(dat,1)
      switch chantype{i}
        case 'MEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T');
        case 'Reference'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T');
        case 'GRAD'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'T/m');
        case 'EEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'SEEG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'EMG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        case 'ECG'
          dat(i,:) = dat(i,:) * ft_scalingfactor(hdr.chanunit{i}, 'V');
        otherwise
          % FIXME I am not sure what scaling to apply
      end
    end

    [p, f, x] = fileparts(filename);
    filename = fullfile(p, f); % without extension
    mat2ades(dat, filename, hdr.Fs, hdr.label, chantype, dattype, datunit);

  case 'homer_nirs'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % https://www.nitrc.org/plugins/mwiki/index.php/homer2:Homer_Input_Files#NIRS_data_file_format
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    % convert the input arguments into a FieldTrip raw data structure
    data = [];
    data.hdr = hdr;
    data.trial{1} = dat;
    data.time{1} = ((1:hdr.nSamples*hdr.nTrials)-1)/hdr.Fs;
    data.label = hdr.label;
    data.sampleinfo = [1 size(dat,2)];

    % convert the raw data structure to Homer format
    nirs = fieldtrip2homer(data, 'event', evt);

    % Homer files are MATLAB files in disguise
    % see https://www.nitrc.org/plugins/mwiki/index.php/homer2:Homer_Input_Files#NIRS_data_file_format
    save(filename, '-struct', 'nirs'); % save the fields as individual variables in the file

  case 'snirf'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % https://github.com/fNIRS/snirf
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if append
      ft_error('appending data is not yet supported for this data format');
    end

    % this uses the SNIRF reading functions from the Homer3 toolbox
    ft_hastoolbox('homer3', 1);

    % construct a time axis that matches the data, it starts at 0 seconds
    time = ((1:hdr.nSamples*hdr.nTrials)-1)/hdr.Fs;

    % divide data in nirs channels, stimulus channels and auxillary channels
    seldat  = startsWith(hdr.chantype, 'nirs');
    selstim = strcmp(hdr.chantype, 'stimulus');
    selaux  = ~seldat & ~selstim;

    % create empty SnirfClass
    snirf = SnirfClass();

    % collect information for creation of snirf file
    source_idx   = find(contains(hdr.opto.optotype, {'transmitter', 'source'}));
    detector_idx = find(contains(hdr.opto.optotype, {'receiver', 'detector'}));
    tra = hdr.opto.tra;
    tra_t = tra'; % transpose tra matrix to get indices of wavelength by row (thus by channel)
    wl_idx = find(tra_t>0);
    all_wavelengths = hdr.opto.wavelength(tra_t(wl_idx));
    split = nanmedian(all_wavelengths);
    WL1.values = all_wavelengths(all_wavelengths<split);
    WL2.values = all_wavelengths(all_wavelengths>split);
    WL1.nominal = round(median(WL1.values),-1);
    WL2.nominal = round(median(WL2.values),-1);
    ft_warning('assuming that the nominal wavelengths are %d and %d nm', WL1.nominal, WL2.nominal)

    % metaDataTags
    snirf.metaDataTags(1).tags.LengthUnit = hdr.opto.unit;
    snirf.metaDataTags(1).tags.TimeUnit = 's';
    snirf.metaDataTags(1).tags.FrequencyUnit = 'Hz';

    % data
    snirf.data(1).dataTimeSeries = dat(seldat,:)'; % <number of time points> x <number of channels>
    snirf.data(1).time = time';                    % <number of time points x 1> (can also be  represented as <start time x sample time spacing>

    % measurementList
    for i=1:size(tra,1)
      source = find(tra(i,:)>0);
      detector = find(tra(i,:)<0);
      snirf.data.measurementList(i).sourceIndex   = find(source_idx==source);
      snirf.data.measurementList(i).detectorIndex = find(detector_idx==detector);
      %     snirf.data.measurementList(i).wavelengthActual = all_wavelengths(i); % this is not yet supported by the snirf toolbox
      if any(all_wavelengths(i)==WL1.values)
        snirf.data.measurementList(i).wavelengthIndex = 1;
      else
        snirf.data.measurementList(i).wavelengthIndex = 2;
      end
      snirf.data.measurementList(i).dataType = 99999;
      snirf.data.measurementList(i).dataTypeLabel = 'dOD';
    end
    ft_warning('assuming that the input data represents (change in) optical densities')

    % sort the channels according to wavelengths, because this is the way that homer handles data
    [dum, idx] = sort(([snirf.data.measurementList(:).wavelengthIndex]));
    % update the data accordingly
    snirf.data.measurementList = snirf.data.measurementList(idx);
    snirf.data.dataTimeSeries=snirf.data.dataTimeSeries(:, idx);

    % stim
    if ~isempty(evt)
      % select events with a string value, the type will be a string
      sel = cellfun(@ischar, {evt.value});
      if any(sel)
        evt_string = evt(sel);
        evt_names = unique({evt_string(:).value}); % if the values are strings, this propably contains the event names
        for i=1:length(evt_names)
          snirf.stim(i).name = evt_names{i};
          evt_idx = find(strcmp({evt_string(:).value}, evt_names{i}));
          starttime = ([evt_string(evt_idx).sample]-1)/hdr.Fs;
          duration = [evt_string(evt_idx).duration]/hdr.Fs;
          if isempty(duration)
            duration = zeros(1, length(starttime));
          end
          value = ones(1,length(evt_idx));
          snirf.stim(i).data = [starttime' duration' value'];
        end
      end
      % select events with a numeric value, the type will be a string
      sel = cellfun(@isnumeric, {evt.value});
      if any(sel)
        evt_numeric = evt(sel);
        evt_names = unique({evt_numeric(:).type});
        for i=1:length(evt_names)
          snirf.stim(i).name = evt_names{i};
          evt_idx = find(strcmp({evt_numeric(:).type}, evt_names{i}));
          starttime = ([evt_numeric(evt_idx).sample]-1)/hdr.Fs;
          duration = [evt_numeric(evt_idx).duration]/hdr.Fs;
          if isempty(duration)
            duration = zeros(1, length(starttime));
          end
          value = [evt_numeric(evt_idx).value];
          if isempty(value)
            value = ones(1, length(starttime));
          end
          snirf.stim(i).data = [starttime' duration' value'];
        end
      end
    end

    % probe
    snirf.probe(1).wavelengths = [WL1.nominal WL2.nominal];
    if all(hdr.opto.optopos(:,3)==0)
      snirf.probe(1).sourcePos2D    = hdr.opto.optopos(source_idx, 1:2);
      snirf.probe(1).detectorPos2D  = hdr.opto.optopos(detector_idx, 1:2);
    else
      snirf.probe(1).sourcePos3D    = hdr.opto.optopos(source_idx, 1:3);
      snirf.probe(1).detectorPos3D  = hdr.opto.optopos(detector_idx, 1:3);
      layoutpos = getorthoviewpos(hdr.opto.optopos, 'ras', 'superior');
      snirf.probe(1).sourcePos2D    = layoutpos(source_idx, 1:2);
      snirf.probe(1).detectorPos2D  = layoutpos(detector_idx, 1:2);
    end
    snirf.probe(1).sourceLabels     = hdr.opto.optolabel(source_idx);
    snirf.probe(1).detectorLabels   = hdr.opto.optolabel(detector_idx);

    % aux
    if sum(selaux)~=0
      auxdata = dat(selaux,:);
      auxlabels = hdr.label(selaux);
      for i=1:sum(selaux)
        snirf.aux(i).name = auxlabels{i}; % check if correct format
        snirf.aux(i).dataTimeSeries = auxdata(i,:)';
        snirf.aux(i).time = time';
      end
    end

    % save .snirf file
    snirf.Save(filename)

  otherwise
    ft_error('unsupported data format');
end % switch dataformat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function type = adestype(type)
for i=1:numel(type)
  switch lower(type{i})
    case 'meggrad'
      type{i} = 'MEG'; % this is for CTF and BTi/4D
    case 'megmag'
      type{i} = 'MEG'; % this is for Neuromag and BTi/4D
    case 'megplanar'
      type{i} = 'GRAD'; % this is for Neuromag
    case {'refmag' 'refgrad'}
      type{i} = 'Reference';
    case 'eeg'
      type{i} = 'EEG';
    case {'seeg' 'ecog' 'ieeg'}
      type{i} = 'SEEG'; % all intracranial channels
    case 'ecg'
      type{i} = 'ECG';
    case 'emg'
      type{i} = 'EMG';
    case 'trigger'
      type{i} = 'Trigger';
    case 'source'
      type{i} = 'Source'; % virtual channel
    otherwise
      type{i} = 'Other';
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Population,RankSolution] = EnvironmentalSelection_noCon(Population,N,alpha,gamma,para)

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Parameter
    popSize   = length(Population);
    NumSeq    = (1:popSize)';
    RankConvg = zeros(popSize,1);
    RankDivs  = zeros(popSize,1);

    %% Modify the infeasible solutions
    PopObj = Population.objs;
    PopCon = Population.cons;
    z      = min(PopObj,[],1);
    n      = max(PopObj,[],1);

    Infeasible_all = any(PopCon>0,2);
    phi_max = max(sum(max(0,PopCon(Infeasible_all,:)),2));

    M          = length(z);
    [W,~]      = UniformPoint(N,M);
    [~,Region] = min(pdist2(PopObj-z,W,'cosine'),[],2);
    PopObj_2   = PopObj;
    for i = 1:size(W,1)
        index = find(Region==i);
        if (~isempty(index))
            Objs_temp  = PopObj_2(index,:);
            Cons_temp  = PopCon(index,:);
            Infeasible = any(Cons_temp>0,2);
            if (sum(Infeasible)~=0)
                F_max = max(Objs_temp,[],1);
                PopObj_2(index(Infeasible),:) = Objs_temp(Infeasible,:)+(sum(max(0,Cons_temp(Infeasible,:)),2)/phi_max).^(exp(para)/max(gamma,0.000001)).*(F_max-Objs_temp(Infeasible,:));
            end
        end
    end

    %% Non-dominated sorting
    Dominate = false(popSize);
    for i = 1 : popSize-1
        for j = i+1 : popSize
            k = any(PopObj_2(i,:)<PopObj_2(j,:)) - any(PopObj_2(i,:)>PopObj_2(j,:));
            if k == 1
                Dominate(i,j) = true;
            elseif k == -1
                Dominate(j,i) = true;
            end
        end
    end

    %% Calculate S(i)
    S = sum(Dominate,2);

    %% Calculate R(i)
    R = zeros(1,popSize);
    for i = 1 : popSize
        R(i) = sum(S(Dominate(:,i)));
    end
    FrontNo = R + 1;

    %% Calculate the crowding distance of each solution
    PopObj   = Population.objs;
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Distance = sort(Distance,2);
    CrowdDis = Distance(:,floor(sqrt(popSize)));

    %% Add a middle column
    MiddleLevel = zeros(popSize,1);

    %% Environmental selection
    Next = FrontNo == 1;
    if sum(Next) <= N
        [~,indx_Convg] = sortrows([FrontNo',-CrowdDis]);
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
        MiddleLevel(Temp(Del)) = 1;
        [~,indx_Convg] = sortrows([FrontNo',MiddleLevel,-CrowdDis]);
    end
    RankConvg(indx_Convg) = NumSeq;

    %% Environmental selection -- diversity
    FrontNo_D = ones(popSize,1);
    for i = 1:size(W,1)
        index = find(Region==i);
        if (~isempty(index))
            Objs_temp = PopObj_2(index,:);
            g_temp = sum((Objs_temp-z).*W(i,:),2);
            [~,index_FrontNo_D] = sort(g_temp);
            FrontNo_D(index(index_FrontNo_D)) = (1:length(g_temp))';
        end
    end
    [~,indx_divs] = sortrows([FrontNo_D,-CrowdDis]);
    RankDivs(indx_divs) = NumSeq;

    %% Population for next generation
    RankSolution = alpha*RankConvg+(1-alpha)*RankDivs;
    [~,Rank]     = sort(RankSolution);
    Population   = Population(Rank(1:N));
    RankSolution = 1 : N;
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO