text:
function [N,Mmax,catalog,IDAll,HistoryAll] = ETESProject3(cat,vEtas,vAfter,StartDateProj,EndDate,FaultParam)

    %This program is like ETESProject.m except that background earthquakes are
    %taken from an input grid, in the same manner as in TotalAftSim.8, instead
    %of smoothed directly from the catalog.  See the companion MakeRateGrid.m
    %program also in this directory.
    %
    %This program starts with A real earthquake catalog and then
    %projects what seismicity will follow further in time.  This is A Monte
    %Carlo modeler -- ideally it should be run multiple times to get average
    %results and error.  It also has A background component.  The background
    %component is set by BckgrndR -- this gives the number of earthquakes that
    %you want each year as
    %background earthquakes.  The background earthquakes will have an r^1.37
    %distribution in one dimension distribution away from hypocenters in the
    %input catalog which are identified as being background earthquakes
    %themselves (e.g. not within one fault length and 10 years of A bigger
    %earthquake).
    %
    %Earthquakes that are smaller than magb are treated as point sources;
    %earthquakes that are larger than magb are modeled as planes.  Where the
    %middle of the fault plane is not well constrained we use the median of the
    %first two days of the aftershock sequence.  For simplicity, large
    %simulated aftershocks are assigned the same focal mechanisms as their
    %mainshocks.  Where the mainshocks are too small for their focal mechanism
    %to be listed we use 75% 332 degree strike aqnd 25% 242 degree strike,
    %based on approximate assesement of trends in the Southern California
    %catalog, with A 90 degree dip.
    %
    %Written by Karen Felzer, 2007.
    %
    %Input:
    %
    % cat: Real eartquake catalog.  Ten column format: year, month, day, hour,
    % minute, second, lat, lon, depth, magnitude.
    %
    %RateGrid: Background rate, spatially varying.  minlon,maxlon, minlat,
    %maxlat, earthquake rate of M>=4.
    %
    %CX, P, A: Omori's law parameters, c, p, and A.  They should be entered
    %as direct triggering parameters, with A recommended p = 1.37.  The value
    %of A should be appropriate for the rate of A mainshock of magnitude M
    %triggering aftershocks of magnitude >=M.  For Mmin = 2.5, we recommend CX = 0.095, P = 1.34, and
    %A = 0.008. (from Felzer and Kilb, in preparation, 2007).
    %
    % Mmin: Minimum earthquake magnitude to report and simulate
    % MminR: Minimum earthquake magnitude for reporting in the output catalog
    % maxm: Maximum earthquake magnitude
    % DMax: Maximum distance for aftershocks
    % EndDate: Last date for the simulated catalog, to be entered as float
    %StartDate; First date for the projected catalog
    % (e.g. 1975.54  'January 1, 2007').  The simulated catalog will begin immediately
    % after the last date in the input catalog cat.
    %
    %magb: Cutoff magnitude for the point source vs. plane source
    %representation of the mainshock
    %
    %FaultParam: Fault parameters for the M>magb earthquakes, listed in order
    %of their occurrence in the catalog.  Columns give fault length, width,
    %strike, dip, year, month, day.
    %
    %
    %OUTPUT:
    %
    %N: total number of aftershocks in the simulated catalog.  Comparing N with
    %the total length of the catalog produced gives the percentage of the total
    %that is aftershocks.
    %
    %MMax: Maximum magnitude in the simulated catalog
    %
    %catalog: The simulated catalog
    %columns 1-10: years, month, day, hour, minute, second, lat, lon, depth,
    %magnitude.  Column 11: ID number for this earthquake.  Column12: If this
    %earthquake is an aftershock, column 12 gives the ID number of its
    %mainshock.  Column 13: If the earthquake is an aftershock, the distance
    %between its epicenter and the nearest point on the fault plane of the
    %mainshock that triggered it.
    %
    %IDAll: The ID numbers for each earthquake; same output as is listed in
    %Column 11 of the catalog
    %
    %HistoryAll: Mainshock of each earthquake; same output as is listed in
    %Column 13 of the catalog.

    %%%%Starting the program
    %set up
    CX=vEtas(1);
    P=vEtas(2);
    A=vEtas(3);

    Mmin=vAfter(1);
    MminR=vAfter(2);
    maxm=vAfter(3);
    DMax=vAfter(4);
    magb=vAfter(5);

    N = zeros(1,1);
    Mmax = zeros(1,1);
    tdistAll = [];

    [fYr1, nMn1, nDay1, nHr1, nMin1, nSec1]=decyear2mat(StartDateProj);
    StartDateProj=datenum(floor(fYr1),nMn1,nDay1,nHr1, nMin1,nSec1);
    [fYr2, nMn2, nDay2, nHr2, nMin2, nSec2]=decyear2mat(EndDate);
    EndDate=datenum(floor(fYr2),nMn2,nDay2,nHr2, nMin2,nSec2);

    T = EndDate - StartDateProj;  %Gives longest time that aftershocks may be generated over; for simplicity all aftershocks will
    %be generated over this time and then be trimmed as needed.

    T3 = EndDate;   %Gives the end time of the simulation

    T2 = StartDateProj;  %The starting date of the projected catalog

    %And limiting the big earthquakes in FaultParam to before the beginning of
    %the simulation

    % tp = datenum(FaultParam(:,5),FaultParam(:,6),FaultParam(:,7));
    % xp = find(tp>T2 | tp<datenum(StartDate));
    % changes vst
    FaultParam = [NaN NaN NaN NaN NaN NaN NaN];
    % FaultParam(xp,:) = [];

    %Finding the average number of aftershocks expected in time T for M - Maft
    %= 0.

    if(P ~= 0)
        NM = (A/(1-P))*((T+CX)^(1-P) - CX^(1-P));
    else
        NM = A*(log(T+CX) - log(CX));
    end

    %Setting up rotation matrices for later so that simulated M>magb mainshocks will have
    %either 332 degree strike or 242 degree strike.  Based on A quick eye
    %assessment of the So Cal fault map, 75% of mainshocks will be assigned the
    %332 strike and 25% the 242 strike.

    z1 = deg2rad(303);
    lambda1 = [cos(z1) sin(z1); -sin(z1) cos(z1)];

    z2 = deg2rad(213);
    lambda2 = [cos(z2) sin(z2); -sin(z2) cos(z2)];

    %And converting catalog times to Julian days

    tdist = datenum(cat(:,1),cat(:,2),cat(:,3),cat(:,4),cat(:,5),cat(:,6));

    m = cat(:,10);

    %eliminating catalog earthquakes that are past the start of the projection

    A=[];
    cat(A,:) = [];
    tdist(A) = [];
    m(A) = [];

    %And changing the lat, lon to xy
    catxy = cat;
    locCen(1) = mean(cat(:,7));
    locCen(2) = mean(cat(:,8));
    catxy(:,7:8) = latlon2xy2(catxy(:,7:8),locCen);
    Locs = [catxy(:,7) catxy(:,8) catxy(:,9)];
    NBack=size(Locs,1);

    %And initiating ID, history, and fault distance tracking matrices, and
    %putting magnitudes in A vector
    ID = [1:1:length(tdist)];
    History = zeros(size(ID));
    FaultDist = -1*ones(size(ID));

    aa2 = find(m>=MminR & tdist>=T2 & tdist<T3);  %we only report earthquakes in the final catalog for m>=MminR

    if(~isempty(aa2))
        tdistAll(1:length(aa2)) = tdist(aa2)';
        LocsAll(1:length(aa2),:) = Locs(aa2,:);
        mAll(1:length(aa2),:) = m(aa2);
        IDAll(1:+length(aa2),:) = ID(aa2)';
        HistoryAll(1:length(aa2),:) = 0;
        FaultDistAll(1:+length(aa2),:) = -1;
    end

    % And doing some final setting up for the ETAS run
    N = 0;
    Mmax = 0;
    ID = ID';
    History = History';
    FaultDist = FaultDist';

    %Now running the ETAS simulation***********%%%%%%%%%%%%

    %The loop gets all of the aftershocks of each earthquake occurring within
    %the specified duration of the catalog and over the minimum magnitude for
    %the simulation

    xm = find(m>=MminR);
    maxID = max(ID);

    countAll = length(tdistAll)+1;

    A = find(m>=Mmin & tdist<=T3);

    NIttr = 1;

    while(~isempty(A))

        IDold = ID(A);

        %Running the small point source and larger extended source earthquakes
        %separtely, through different versions of the aftershock program.
        %First, the small point source earthquakes

        a2 = find(m(A)<magb);

        a2 = A(a2);

        tdistSave = tdist;
        IDSave = ID;
        mSave = m;
        LocsSave = Locs;

        try
            [tdist,m,Locs,ID,History,FaultDist] = GAft3(tdist(a2),maxID,ID(a2),m(a2),Mmin,maxm,NM,T,Locs(a2,:),DMax,P,CX);
        catch
            disp('No point source events');
        end
        %And then doing the larger, plane source earthquakes

        a3 = find(mSave(A)>=magb);

        a3 = A(a3);

        if(~isempty(a3))

            if(NIttr>0)  %These are simulated large earthquakes which need FaultParam assigned.  Doing widths and lengths from Wells & Coppersmith

                FaultParam = BuildFaultParam(mSave(a3),LocsSave(a3,:),lambda1,lambda2);
            end

            %And getting the aftershocks for the big mainshocks

            try
                [tdist2,m2,Locs2,ID2,History2,FaultDist2] = GAft4(tdistSave(a3),maxID,IDSave(a3),mSave(a3),Mmin,maxm,NM,T,LocsSave(a3,:),DMax,FaultParam,P,CX);
            catch
                disp('No plane source events');
            end

        else
            tdist2 = [];
            m2 = [];
            Locs2 = [];
            ID2 = [];
            History2 = [];
            FaultDist2 = [];
        end

        %And putting the results from the two runs together

        tdist = [tdist; tdist2];

        m = [m; m2];
        Locs = [Locs; Locs2];
        ID = [ID ID2];
        History = [History'; History2];
        FaultDist = [FaultDist; FaultDist2];

        maxID = max([maxID max(ID)]);

        A = find(tdist>T2 & tdist<=T3);

        aa2 = find(m(A)>=MminR);  %we only report earthquakes in the final catalog for m>=MminR

        aa2 = A(aa2);

        if(~isempty(aa2))

            tdistAll(countAll:countAll+length(aa2)-1) = tdist(aa2)';

            LocsAll(countAll:countAll+length(aa2)-1,:) = Locs(aa2,:);

            mAll(countAll:countAll+length(aa2)-1,:) = m(aa2);

            IDAll(countAll:countAll+length(aa2)-1,:) = ID(aa2)';

            HistoryAll(countAll:countAll+length(aa2)-1,:) = History(aa2);

            FaultDistAll(countAll:countAll+length(aa2)-1,:) = FaultDist(aa2);

        end

        countAll = countAll+length(aa2);

        N = N + length(aa2);

        if(max(m(A))>Mmax)
            Mmax = max(m(A));
        end

        N;
        Mmax;

        NIttr = NIttr + 1;

    end

    %And translating the results into A traditional earthquake catalog

    % tdays = datenum(StartDate) + tdistAll;

    tvect = datevec(tdistAll);

    latlon = xy2latlon(LocsAll(:,1:2),locCen);

    catalog = [tvect latlon LocsAll(:,3) mAll IDAll HistoryAll FaultDistAll];

    catalog = sortrows(catalog,[1 2 3 4 5 6]);

    %cat = gquakes(catalog,StartDate,EndDate);

    %At this point solving for the background rate, RateGrid, if it was
    %initially input as zero.
end

    %--------------------------------------%%%

    %Getting aftershock times, locations, and mags for point source mainshocks

function [tdist,m,Locs,ID,History,FaultDist] = GAft3(tdist,maxID,IDold,m,Mmin,maxm,NM,T,Locs,DMax,P,CX)

    %First getting the number of aftershocks produced by each mainshock; stored
    %in the vector A.

    A = NM.*10.^(m - Mmin);

    try
        A = poissinv(rand(size(A,1),1),A);
    catch
        poisstest.NM=NM;
        poisstest.A=A;
        poisstest.m=m;
        save poisstest.mat poisstest -mat
        A = poissinv(rand(size(A,1),1),A);
    end

    counter = 1;

    x = find(A>0);

    tnew2 = [];

    %Making A vector with all of the mainshock times and locations

    if(~isempty(x))

        Axx = A(x);
        Axx = cumsum(Axx);
        Sp = [1; Axx(1:end-1)+1];
        Ep = Axx;
        tdist = tdist(x);
        Locs = Locs(x,:);
        xl = Locs(:,1);
        yl = Locs(:,2);
        zl = Locs(:,3);

        %The loop below is the time limiting factor for the whole calculation in
        %matlab

        AddOn = zeros(Axx(end),1);

        for j=1:length(x)

            AddOn(Sp(j):Ep(j)) = tdist(j);
            lx(Sp(j):Ep(j)) = xl(j);
            ly(Sp(j):Ep(j)) = yl(j);
            lz(Sp(j):Ep(j)) = zl(j);
            History(Sp(j):Ep(j)) = IDold(x(j));

        end

        AddDist = [lx' ly' lz'];

        %IN THIS VERSION: For each aftershock that will be produced, listing the
        %location of the mainshock that will produce it.

        %NOTE: COMMENTING THIS OUT RIGHT NOW TO SAVE TIME

        %for j = 1:length(tdist(x))

        %   idx = A(x(j));

        %   LocsT = [Locs(x(j),1)*ones(idx,1) Locs(x(j),2)*ones(idx,1) Locs(x(j),3)*ones(idx,1)];

        %  AddDist(counter:counter+A(x(j))-1,1:3) = LocsT;
        %  History(counter:counter+A(x(j))-1,1) = IDold(x(j));
        %  counter = counter + A(x(j));

        %end

        %Summing the total number of aftershocks

        AX = sum(A);

        %Now getting distances of all of the aftershocks from the faults
        %Allowing A closest approach of 1 m.

        dnew = GPowDistR(1.37,AX,0.001,DMax);

        FaultDist = dnew;

        %Converting the distances into x, y using A random theta.
        %Although not completely accutate, to save computation time
        %just keeping z the same as the generating point on the fault
        %so that we don't need to worry about staying within seismogenic
        %depth.  NOTE: there is code in the program for Joan that can be used
        %to fix this!  GET IT FIXED!!

        theta = rand(length(dnew),1).*2*pi;

        xnew = dnew.*cos(theta);
        ynew = dnew.*sin(theta);
        znew = zeros(length(ynew),1);

        LocsNew = AddDist + [xnew ynew znew];

    end

    %Now generating the aftershock times

    if(~isempty(x))

        %Then calculating times for all of the aftershocks, in terms of time since
        %the mainshock

        %tnew = GPowDistR(1.37,AX,0.0104,T);

        tnew = GPowDistRc(P,CX,AX,0,T);

        %And finally adding each aftershock time to the time of its mainshock to
        %get the absolute time of the aftershock in the catalog

        tnew2 = tnew + AddOn;

    end

    %And assigning the aftershock magnitudes

    if(~isempty(tnew2))

        m = GetMag(length(tnew2),Mmin,maxm);

        tdist = tnew2;
        Locs = LocsNew;

    else
        tdist = [];
        m = [];
        Locs = [];
        ID = [];
        History = [];
        FaultDist = [];
    end

    %And generating ID numbers for all of the aftershocks

    if(~isempty(x))

        ID = [maxID+1:1:maxID+AX+1];
    end
end
    %---------------------------------------

    %Getting aftershock times, locations, and mags for fault plane mainshocks

function [tdist,m,Locs,ID,History,FaultDist] = GAft4(tdist,maxID,IDold,m,Mmin,maxm,NM,T,Locs,DMax,FaultParam,P,CX)

    %First getting the number of aftershocks produced by each mainshock; stored
    %in the vector A.

    A = NM.*10.^(m - Mmin);

    A = poissinv(rand(length(A),1),A);

    counter = 1;

    x = find(A>0);

    tnew2 = [];

    %And setting up the fault plane

    %And for each aftershock that will be produced, picking A random reference
    %point on the mainshock fault that the distance of the aftershock from the
    %fault will be measured from.

    %IN THIS VERSION: For each aftershock that will be produced, listing the
    %location of the mainshock that will produce it.

    for j = 1:length(tdist(x))

        FaultY = FaultParam(j,4);
        FaultZ = FaultParam(j,5);

        YLow = FaultParam(j,2) - FaultY/2;
        ZLow = FaultParam(j,3) - FaultZ/2;

        %And getting points on the fault that will generate the earthquakes

        rpy = rand(A(x(j)),1).*FaultY;
        rpz = rand(A(x(j)),1)*FaultZ;
        rpx = zeros(size(rpy,1),1);

        rpy = rpy - FaultY/2;

        %Then rotating to the proper strike and dip

        %matrix for rotation around strike

        z1 = deg2rad(FaultParam(x(j),6));
        lambda1 = [cos(z1) sin(z1); -sin(z1) cos(z1)];

        z2 = deg2rad(FaultParam(x(j),7) - 90);  %rotation around dip

        if(FaultParam(x(j),6)>180)
            z2 = deg2rad(90-FaultParam(x(j),7));
        end

        lambda2 = [cos(z2) sin(z2); -sin(z2) cos(z2)];

        %And rotating.  Fist going to A 0 degree strike, and then to A 90 degree
        %dip.
        rot = lambda1*[rpx rpy]';

        rpx = rot(1,:)';
        rpy = rot(2,:)';

        rot = lambda2*[rpz rpy]';

        rpz = rot(1,:)';
        rpy = rot(2,:)';

        Lx = rpx + FaultParam(x(j),1);
        Ly = rpy + FaultParam(x(j),2);
        Lz = rpz + FaultParam(x(j),3);

        %And assigning the points on the fault from which aftershocks will be generated
        %And placing ID values in the history matrix

        AddDist(counter:counter+A(x(j))-1,1:3) = [Lx Ly Lz];
        History(counter:counter+A(x(j))-1,1) = IDold(x(j));
        counter = counter + A(x(j));

    end

    if(~isempty(x))

        %Summing the total number of aftershocks

        AX = sum(A);

        %Now getting distances of all of the aftershocks from the faults
        %Allowing A closest approach of 1 m.

        dnew = GPowDistR(1.37,AX,0.001,DMax);

        FaultDist = dnew;

        %Converting the distances into x, y using A random theta.
        %Although not completely accurate, to save computation time
        %just keeping z the same as the generating point on the fault
        %so that we don't need to worry about staying within seismogenic
        %depth.  NOTE: there is code in the program for Joan that can be used
        %to fix this!  GET IT FIXED!!

        theta = rand(length(dnew),1).*2*pi;

        xnew = dnew.*cos(theta);
        ynew = dnew.*sin(theta);
        znew = zeros(length(ynew),1);

        LocsNew = AddDist + [xnew ynew znew];

    end

    %Now generating the aftershock times

    %Making A vector with all of the mainshock times

    counter = 1;

    for j = 1:length(tdist(x))
        AddOn(counter:counter+A(x(j))-1) = tdist(x(j));
        counter = counter + A(x(j));
    end

    if(~isempty(x))

        %Then calculating times for all of the aftershocks, in terms of time since
        %the mainshock

        %tnew = GPowDistR(1.37,AX,0.0104,T);

        tnew = GPowDistRc(P,CX,AX,0,T);

        %And finally adding each aftershock time to the time of its mainshock to
        %get the absolute time of the aftershock in the catalog

        tnew2 = tnew + AddOn';

    end

    %And getting aftershock magnitudes

    if(~isempty(tnew2))

        m = GetMag(length(tnew2),Mmin,maxm);

        tdist = tnew2;
        Locs = LocsNew;

    else
        tdist = [];
        m = [];
        Locs = [];
        ID = [];
        History = [];
        FaultDist = [];
    end

    %And generating ID numbers for all of the aftershocks

    if(~isempty(x))

        ID = [maxID+1:1:maxID+AX+1];
    end
end

    %------------------------------------

    %Getting aftershock magnitudes

function m = GetMag(L,Mmin,maxm)

    m = Mmin - log10(rand(L,1));

    x = find(m>maxm);

    for i=1:length(x)

        while(m(x(i))>maxm)
            m(x(i)) = Mmin - log10(rand(1));

        end
    end
end

author: CelsoReyes


text:
function roi = dtiRoiFromImg(roiImg, imgXform, bb)
%
% roi = dtiRoiFromImg(roiImg, imgXform, bb)
%
% If not the identity matrix, imgXform is typically set to xformToAcpc. This
% has the effect of forcing the roiImg to have the same voxel size as the
% image that xformToAcpc is based on.
%
% bb is the bounding box, defined in the xformed spaced. Defaults to
% [min(coords)-10; max(coords)+10].
%
% [roiImg, imgXform, bb] = dtiRoiToImg(roi);
% % Do some processing on the ROI
% perimImg = bwperim(roiImg);
% perimRoi = dtiRoiFromImg(roiImg, imgXform, bb);
%
% HISTORY:
% 2010.03.11 RFD wrote it.

[coords(:,1), coords(:,2), coords(:,3)] = ind2sub(size(roiImg), find(roiImg));
coords(:,1) = coords(:,1) + bb(1,1) - 1;
coords(:,2) = coords(:,2) + bb(1,2) - 1;
coords(:,3) = coords(:,3) + bb(1,3) - 1;
coords = mrAnatXformCoords(inv(imgXform), coords);
roi = dtiNewRoi('roiFromImg','g', coords);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout=mod(varargin)
%MOD (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        if ~isa(varargin{2},'double')
            error('MOD is currently only supported for DOUBLE second argument');
        end

        x = varargin{1};
        y = varargin{2};
        % Some boring code for scalarization
        if prod(size(x))==1 & prod(size(y))>1
            x = x*ones(size(y));
        elseif prod(size(y))==1 & prod(size(x))>1
            y = y*ones(size(x));
        end
        if ~all(size(x) == size(y))
            error('Matrix dimensions must agree.');
        end
        dim = size(x);
        x = reshape(x,prod(dim),1);
        y = reshape(y,prod(dim),1);
        z = [];
        % Create one variable for each element
        for i = 1:length(x)
            xi = extsubsref(x,i);
            yi = extsubsref(y,i);
            inarg = {xi,yi};
            z = [z;yalmip('define',mfilename,inarg{:})];
        end
        z = reshape(z,dim);
        varargout{1} = z;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph

        % Description using epigraphs
        t = varargin{2};
        x = varargin{3};
        y = varargin{4};

        % t = mod(x,y), i.e. t = x - n*y, n = floor(x/y)
        n = intvar(1,1);
        [M,m] = derivebounds(x);
        F_floor = [(x/y) - 1 <= n <= (x/y)];
        if ~isinf(m) && isa(y,'double')
            F_floor = [F_floor, floor(m/y) <= n];
        end
        if ~isinf(M) && isa(y,'double')
            F_floor = [F_floor, n <= floor(M/y)];
        end
        F = [t == x - y*n, F_floor];

        varargout{1} = F;
        varargout{2} = struct('convexity','none','monotonicity','none','definiteness','none','model','integer');
        varargout{3} = [x(:);y(:)];

    otherwise
        error('Strange type on first argument in SDPVAR/MOD');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function C = tolUnique(A, tol)
%TOLUNIQUE   UNIQUE with a tolerance for checking floating-point equality.
%   C = TOLUNIQUE(A, TOL) for a real vector A is the same as UNIQUE(A) except
%   that two values are deemed "equal" if they are within an absolute
%   difference TOL of one another.  If this happens, the mean of the two values
%   is placed in C.
%
%   C = TOLUNIQUE(A) uses a default tolerance of 100*EPS*MAX(NORM(A, Inf)).

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

if ( nargin < 3 )
    tol = 100*eps*max(norm(A, Inf));
end

% TODO:  This might not do the right thing if there are multiple values that
% are all close to one another, say, 2, 2 + 2*eps, 2 + 4*eps, and 2 + 6*eps,
% and tol = 3*eps.  This code will merge all of those into one value:  2 + eps
% (the average of 2 and 2 + 2*eps), and one can argue that behavior is wrong.
% We can resolve this if this ever becomes an issue in practice.
C = unique(A);
ind = find(diff(C) < tol);
C(ind) = (C(ind) + C(ind+1))/2;
C(ind+1) = [];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [pos, tri] = remove_double_vertices(pos, tri)

% REMOVE_DOUBLE_VERTICES removes double vertices from a triangular, tetrahedral or
% hexahedral mesh, renumbering the vertex-indices for the elements.
%
% Use as
%   [pos, tri] = remove_double_vertices(pos, tri)
%   [pos, tet] = remove_double_vertices(pos, tet)
%   [pos, hex] = remove_double_vertices(pos, hex)
%
% See also REMOVE_VERTICES

% Copyright (C) 2004-2022, Robert Oostenveld and Jan-Mathijs Schoffelen
%
% 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$

npos = size(pos, 1);
[dum, keeppos, i2] = unique(pos, 'rows');
clear dum

numb    = zeros(1,npos);
numb(keeppos) = 1:length(keeppos);

% re-index the indices in tri
tri = keeppos(i2(tri));
tri = numb(tri);

% remove the vertices and triangles
pos = pos(keeppos, :);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
seed = 1;
rand('state', seed);
randn('state', seed);

obs_model = 'unique';  % each cell has a unique label (essentially fully observable)
%obs_model = 'four'; % each cell generates 4 observations, NESW

% Generate the true network, and a randomization of it
realnet = mk_map_hhmm('p', 0.9, 'obs_model', obs_model);
rndnet = mk_rnd_map_hhmm('obs_model', obs_model);
eclass = realnet.equiv_class;
U = 1; A = 2; C = 3; F = 4; onodes = 5;

ss = realnet.nnodes_per_slice;
T = 100;
evidence = sample_dbn(realnet, 'length', T);
ev = cell(ss,T);
ev(onodes,:) = evidence(onodes,:);

infeng = jtree_dbn_inf_engine(rndnet);

if 0
% suppose we do not observe the final finish node, but only know
% it is more likely to be on that off
ev2 = ev;
infeng = enter_evidence(infeng, ev2, 'soft_evidence_nodes', [F T], 'soft_evidence',  {[0.3 0.7]'});
end

learnednet = learn_params_dbn_em(infeng, {evidence}, 'max_iter', 5);

disp('real model')
disp_map_hhmm(realnet)

disp('learned model')
disp_map_hhmm(learnednet)

disp('rnd model')
disp_map_hhmm(rndnet)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function grads = imLocalGradient(img, pos, sigma)
% Compute gradient for chosen locations within image.
%
%   The aim of this function is to compute gradient within images at
%   specified positions without requiring to compute the whole gradient
%   image. This greatly saves memory when image is large and number of
%   position is small.
%
%   GRAD = imLocalGradient(IMG, POS)
%   Evaluates the gradient within the image IMG at the specified position
%   POS. The result GRAD is a 1-by-2 row vector containing gradient along
%   the x and y diections.
%   If POS is a N-by-2 array of coordinates, the result is a N-by-2 array
%   of vector components.
%
%   GRAD = imLocalGradient(IMG, POS, SIGMA)
%   Also specifies the sigma parameter for computing gradient.
%
%   Example
%     imLocalGradient
%
%   See also
%     imGradientFilter, imLaplacian, gradientKernels, gradientKernels3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-06-01,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE.

%% Input arguments

% process optional inputs
if nargin == 2
    sigma = 2;
end

% compute output dimensions from input arguments
nPos = size(pos, 1);
nd = ndims(img);

% allocate memory for result
grads = zeros(nPos, nd);

%% Main processing

if nd == 2
    % create gradient kernels
    [kx, ky] = gradientKernels(sigma);

    % computation will be performed using point-wise multiplication, so we need
    % to use the symetric kernels
    kx = kx(end:-1:1, end:-1:1);
    ky = ky(end:-1:1, end:-1:1);

    % size of the kernels along each dimension
    ks1 = size(kx, 1);
    ks2 = size(kx, 2);
    kr1 = floor((ks1 - 1) / 2);
    kr2 = floor((ks2 - 1) / 2);

    % iterate over position
    for iPos = 1:nPos
        % position in array indexing (position is x,y indexing)
        pos1 = round(pos(iPos, 2));
        pos2 = round(pos(iPos, 1));

        % clamp indexing of sub-image
        inds1 = max(min(pos1-kr1:pos1+kr1, size(img, 1)), 1);
        inds2 = max(min(pos2-kr2:pos2+kr2, size(img, 2)), 1);

        % extract region to process (using replication of borders if necessary)
        sub = double(img(inds1, inds2));

        % evaluate gradient by multiplying with gradient kernels
        gx = sub .* kx;
        gy = sub .* ky;
        grads(iPos, :) = [sum(gx(:)) sum(gy(:))];
    end

else
    % create gradient kernels
    [kx, ky, kz] = gradientKernels3d(sigma);

    % computation will be performed using point-wise multiplication, so we need
    % to use the symetric kernels
    kx = kx(end:-1:1, end:-1:1, end:-1:1);
    ky = ky(end:-1:1, end:-1:1, end:-1:1);
    kz = kz(end:-1:1, end:-1:1, end:-1:1);

    % size of the kernels along each dimension
    ks1 = size(kx, 1);
    ks2 = size(kx, 2);
    ks3 = size(kx, 3);
    kr1 = floor((ks1 - 1) / 2);
    kr2 = floor((ks2 - 1) / 2);
    kr3 = floor((ks3 - 1) / 2);

    % iterate over position
    for iPos = 1:nPos
        % position in array indexing (position is x,y indexing)
        pos1 = round(pos(iPos, 2));
        pos2 = round(pos(iPos, 1));
        pos3 = round(pos(iPos, 3));

        % clamp indexing of sub-image
        inds1 = max(min(pos1-kr1:pos1+kr1, size(img, 1)), 1);
        inds2 = max(min(pos2-kr2:pos2+kr2, size(img, 2)), 1);
        inds3 = max(min(pos3-kr3:pos3+kr3, size(img, 3)), 1);

        % extract region to process (using replication of borders if necessary)
        sub = double(img(inds1, inds2, inds3));

        % evaluate gradient by multiplying with gradient kernels
        gx = sub .* kx;
        gy = sub .* ky;
        gz = sub .* kz;
        grads(iPos, :) = [sum(gx(:)) sum(gy(:)) sum(gz(:))];
    end
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% function [Ps,Pb] = pb_qam_koh_ray (EbN0db,M,L,map);
% ------------------------------------------------------------------------------
% EINGABE:
%   EbN0db: Signal/Noise pro Infobit [db]
%           (Skalar, Spalten- oder Zeilenvektor)
%        M: Stufigkeit der Modulation (default=4)
%           (Skalar)
%        L: Anzahl der Pfade (default=1)
%           (Skalar)
%      map: Art des Mappings
%           'gray' (default)
%             Gray-Codierung
%           'nat'
%             natuerliches Mapping
%           (string)
%
% AUSGABE:
%   Ps: Symbolfehlerwahrscheinlichkeit (Spaltenvektor)
%   Pb: Bitfehlerwahrscheinlichkeit (Spaltenvektor)
%
% ANMERKUNGEN:
%  - benoetigt pb_pam_ray_co.m
%
% AUTOR: Juergen Rinas,  08.12.1998
% Erweitert um Symbolfehlerraten und Mapping von Volker Kuehn, 04.01.01
% ------------------------------------------------------------------------------

function [Ps,Pb] = pb_qam_koh_ray (EbN0db,M,L,map)

if (sscanf(version,'%d')<5)
  disp(['# Achtung: ',mfilename,' wurde unter Matlab 5.1 entwickelt.']);
end;
if (nargin<2) M=4; end;        % default: M=4
if (nargin<3) L=1; end;        % default: L=1
if (nargin<4) map='gray'; end; % default: map='gray'
K=log2(M);
if (K~=fix(K))
  disp(['# Achtung(',mfilename,') M ist keine Potenz von 2.']);
end;
sqrtM=sqrt(M);
if (sqrtM~=fix(sqrtM))
  disp(['# Achtung(',mfilename,') sqrt(M) ist keine ganze Zahl.']);
end;

[tmp,Pb]=pb_pam_ray(EbN0db,sqrt(M),L,map);
Ps=2*tmp-tmp.^2;

% EOF --------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function out=real(x)

precision=x(1).precision;
out_rval=cell(size(x));
out_ival=cell(size(x));

for ii=1:numel(x)
 imag=false;
 [out_rval{ii},xival]=getVals(x,ii);
end % for ii=1:max(ex,
out=class(struct('rval',out_rval,...
                 'ival',out_ival,...
                 'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
function test_suite = test_distancePointEdge3d
%TEST_DISTANCEPOINTEDGE3D  Test case for the file distancePointEdge3d
%
%   Test case for the file distancePointEdge3d

%   Example
%   test_distancePointEdge3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-04-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Simple(testCase) %#ok<*DEFNU>
% Test call of function without argument

edge = [10 10 10 30 10 10];

p0 = [10 10 10];
exp0 = 0;
testCase.assertEqual(exp0, distancePointEdge3d(p0, edge), 'AbsTol', .01);

p1 = [30 10 10];
exp1 = 0;
testCase.assertEqual(exp1, distancePointEdge3d(p1, edge), 'AbsTol', .01);

p2 = [0 0 0];
exp2 = 10 * sqrt(3);
testCase.assertEqual(exp2, distancePointEdge3d(p2, edge), 'AbsTol', 1e-14);

p3 = [40 20 20];
exp3 = 10 * sqrt(3);
testCase.assertEqual(exp3, distancePointEdge3d(p3, edge), 'AbsTol', 1e-14);

function testArray(testCase)

edge = [10 10 10 30 10 10];
pts = [10 10 10;30 10 10;0 0 0;40 20 20];
exp = [0 0 10*sqrt(3) 10*sqrt(3)]';
testCase.assertEqual(exp, distancePointEdge3d(pts, edge), 'AbsTol', 1e-14);

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function dz = singlePendulumRhs(~,z,g,l)
% This function is used inside of ode45 for running a simulation of a
% single pendulum. The first argument (time) is not used.

th = z(1,:);
w = z(2,:);

dth = w;
dw = singlePendulumDynamics(th,g,l);

dz = [dth;dw];

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function data = rrgen(tmax, Pe, Pn, seed, filename);
% data = rrgen(tmax, Pe, Pn, seed, 'filename');
% generates artificial rr time series over a time tmax seconds
% and dumps the output to 'filename'. If filename is not specified,
% the output is returned into the data array. If it is specified, then the
% data array is zero-valued and of length(tmax). Pe and Pn are the
% probability of ectopy and noise.
% Defaults : tmax=300; NoEctopic=1; NoNoise=1; seed=1;
% P(ectopy), Pe = 0.0003;  ...  Probability of ectopy ~ 1 per hr
% P(noise),  Pn = 0.0048;  ...  Probability of noise ~ 16 per hr
% P(noise) in sleep is Pn/12, P(ectopy) in sleep is unchanged.
%
% Requires rrgen binary - compilation of rrgen.c on your system:
% gcc -Wall rrgen.c -lm -o rrgen
%
% (C) P.E. McSharry & G.D. Clifford 2002 under the
% GNU Public License. Contact P. McSharry (patrick@mcsharry.net)
% or G. Clifford (gari@mit.edu)
% Also see - McSharry P.E., Clifford G.D., Tarassenko L.,
% Smith L.: Method for generating an artificial RR tachogram of a typical
% healthy human over 24-hours, Computers in Cardiology, 29:225-228, IEEE
% Computer Society Press, September 2002.

if nargin < 5
  filename = 'tmp.dat';
end
if nargin < 4
  seed=1;
end
if nargin < 3
    Pn=0.0048;
end
if nargin < 2
    Pe=0.0003;
end
if nargin < 1
  tmax=300;
end

if tmax<115
    fprintf('Warning, tmax minimum value is 115, using this as default\n')
end

tmax=round(tmax);
eval(['!rrgen ' int2str(seed) ' ' int2str(tmax) ' ' int2str(Pe) ' ' int2str(Pn) ' > ' filename]);
fid = fopen(filename,'r');
data = zeros(1,tmax);
for(i=1:tmax)
 fscanf(fid,'%f\n',data(i));
end
fclose(fid);

% if not file name read in from tmp.dat and then delete it.
if nargin < 5
 load tmp.dat;
 data = tmp;
 if(ispc==1)
  system('del tmp.dat');
 else
  system('rm tmp.dat');
 end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function calcAllSeparateTextures_batchHN(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,scale_mat,algo_cell,Ng_mat,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% function calcAllSeparateTextures_batchHN(pathData,pathText,namePT,nameCT,nameROI,outcomes,featType,scale_mat,Ng_mat,nBatch,matlabPATH)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes SEPARATE texture features for all patients, for
% all different combinations of the following texture extraction parameters:
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Quantization algorith: Type of quantization algorithm used.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for SEPARATE scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx>, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathData: Full path to the HN sData files directory.
%              --> Ex: '/myProject/WORKSPACE/DATA'
% 2. pathText: Full path to the HN non texture features directory.
%              --> Ex: '/myProject/WORKSPACE/FEATURES/TEXTURES'
% 3. namePT: Cell of strings of all PET sData files to read
%            --> Ex: {'HGJ_001_PT.PTscan.mat';'HGJ_022_PT.PTscan.mat'}
% 4. namePT: Cell of strings of all CT sData files to read
%            --> Ex: {'HGJ_001_CT.CTscan.mat';'HGJ_022_CT.CTscan.mat'}
% 5. nameROI: Cell of strings specifying the ROI names to analyze for the
%             patients defined by "namePT" and "nameCT"
%             --> Ex: {'GTV';'GTV-P'}
% 6. outcomes: Structure specifying the status (1 or 0) for different
%              outcomes in HN cancer. Contains: outcomes.Failure,
%              outcomes.Locoregional, outcomes.Distant. See ref.[1] for
%              more details.
% 7. featType: Either 'GTVp' for primary GTV, or 'GTVtot' for primaty GTV +
%              nodal GTVs
%              --> Ex: 'GTVp'
% 8. scale_mat: Array vector specifying the different 'Scale' values to test.
%               --> Ex: [1,2,3,4,5]
% 9. Ng_mat: Array vector specifying the different 'Ng' values to test.
%            --> Ex: [8,16,32,64]
% 10. nBatch: Number of parallel batch.
%             --> Ex: 8
% 11.  matlabPATH: Full path to the MATLAB excutable on the system.
%      --> Ex: 'matlab' (symbolic link)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: March 2016
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
ind = strfind(namePT{1},'_'); cohortID = namePT{1}(1:ind-1);
cd(pathText), mkdir(['batchLog_',cohortID,'_',featType,'_SepText']), cd(['batchLog_',cohortID,'_',featType,'_SepText']), pathBatch = pwd;
time = 60; % Number of seconds to wait before checking if parallel computations are done
nameOutcomes = fieldnames(outcomes); nOutcomes = numel(nameOutcomes);
scans = {'PT','CT'}; nScans = numel(scans);

% PRODUCE BATCH COMPUTATIONS
nPatient = numel(namePT); valid = ones(nPatient,1);
for i = 1:nPatient
    if isempty(nameROI{i})
        valid(i) = 0;
    end
end
indValid = find(valid); nPatient = numel(indValid);
if nPatient < nBatch
    nBatch = nPatient;
end
[patients] = batchPatients(nPatient,nBatch);
save('workspace','pathData','pathText','namePT','nameCT','nameROI','patients','indValid','featType','scale_mat','algo_cell','Ng_mat'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    fprintf(fid,['calcAllSeparateTextures_HN(pathData,pathText,namePT(indValid(patients{',num2str(i),'})),nameCT(indValid(patients{',num2str(i),'})),nameROI(indValid(patients{',num2str(i),'})),featType,scale_mat,algo_cell,Ng_mat)\n']);
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

% GROUPING RESULTS FROM ALL BATCH
nPatient = numel(namePT);
names = {namePT,nameCT};
for scan = 1:nScans
    cd(pathText)
    if exist(['HGJ_001_',scans{scan},'_',featType,'_text.mat'],'file')
        temp = load(['HGJ_001_',scans{scan},'_',featType,'_text']); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    else
        temp = load(['HGJ_001_',scans{scan},'_','GTVp','_text']); temp = struct2cell(temp); temp = temp{1}; % In order to get the necessary 'nameType' and 'nameFeature' fields
    end
    nameType = fieldnames(temp.Experiment1); nameType(end) = []; nType = numel(nameType); % All texture types are the same
    text = cell(numel(scale_mat),numel(algo_cell),numel(Ng_mat));
    tempText = cell(1,nPatient); % Cell used to load patient textures only once
    for p = 1:nPatient
        ind = strfind(names{scan}{p},'.'); namePatientScan = names{scan}{p}(1:ind(1)-1);
        if exist([namePatientScan,'_',featType,'_text.mat'],'file')
            load([namePatientScan,'_',featType,'_text']) % Variable 'textures' is now in MATLAB workspace
        else
            load([namePatientScan,'_','GTVp','_text']) % Variable 'textures' is now in MATLAB workspace
        end
        tempText{p} = textures;
    end
    experiment = 0;
    for s = 1:numel(scale_mat)
        for a = 1:numel(algo_cell)
            for n = 1:numel(Ng_mat)
                text{s,a,n} = struct;
                experiment = experiment + 1;
                strExperiment = ['Experiment',num2str(experiment)];
                for t = 1:nType
                    nameFeature = fieldnames(temp.(strExperiment).(nameType{t})); nFeature = numel(nameFeature);
                    for f = 1:nFeature
                        data = zeros(nPatient,1);
                        for p = 1:nPatient
                            data(p,1) = tempText{p}.(strExperiment).(nameType{t}).(nameFeature{f});
                        end
                        text{s,a,n}.(nameType{t}).(nameFeature{f}).Data = data;
                        for o = 1:nOutcomes
                            [text{s,a,n}.(nameType{t}).(nameFeature{f}).Spearman.(nameOutcomes{o}).rs,text{s,a,n}.(nameType{t}).(nameFeature{f}).Spearman.(nameOutcomes{o}).p] = corr(data,outcomes.(nameOutcomes{o}),'type','Spearman');
                        end
                    end
                end
            end
        end
    end
    cd .., save(['text_',cohortID,'_',scans{scan},'_',featType],'text')
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function fignum = xsec_fig()
    % XSEC_FIG returns the handle to the cross-section figure
    fignum=findobj('Type','Figure','-and','Name','Cross -Section');
    if isempty(fignum)
        fignum=figure('Name','Cross -Section');
    end


text:
n = 500;

% You need to return the following variables correctly.
x = zeros(n, 1);

word_indices=[100,212,312,423,345,46,37];

for i=1:length(word_indices)
	if (word_indices(i))
		x(i)=1
	else
		x(i)=0;
	end
end

author: yhyap

repo: machine-learning-coursera


text:
% mff_exportevents - export MFF EEG event from EEGLAB structure to
%                    'Events_exported_from_EEGLAB.xml' file
% Usage:
%   mff_exportevents(EEG, mffFile);
%
% Inputs:
%  EEG     - EEGLAB structure
%  mffFile - filename/foldername for the MFF file (MFF file/folder must
%            already exist)

% This file is part of mffmatlabio.
%
% mffmatlabio 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.
%
% mffmatlabio 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 mffmatlabio.  If not, see <https://www.gnu.org/licenses/>.

function mff_exportevents(EEG, mffFile)

p = fileparts(which('mff_importsignal.m'));
warning('off', 'MATLAB:Java:DuplicateClass');
javaaddpath(fullfile(p, 'MFF-1.2.2-jar-with-dependencies.jar'));
warning('on', 'MATLAB:Java:DuplicateClass');

if isfield(EEG.etc, 'recordingtime')
    begTime = EEG.etc.recordingtime;
else
    begTime = now;
end
if isfield(EEG.etc, 'timezone')
    timeZone = EEG.etc.timezone;
else
    timeZone = '00:00';
end
srate    = EEG.srate;

% create a factory
mfffactorydelegate = javaObject('com.egi.services.mff.api.LocalMFFFactoryDelegate');
mfffactory = javaObject('com.egi.services.mff.api.MFFFactory', mfffactorydelegate);

% determine the number of files
events   = EEG.event;
if isfield(EEG.event, 'name')
    allNames = { events.name };
    allNames = cellfun(@num2str, allNames, 'UniformOutput', false);
    uniqueNames = unique( allNames );
else
    uniqueNames = { '' };
end

tracktype = '';
for iFile = 1:length(uniqueNames)

    % Create Signal object and read in event track file.
    fileName = 'Events_exported_from_EEGLAB.xml';
    if length(uniqueNames) > 1
        eventtrackfilename = fullfile(mffFile, [ fileName(1:end-4) '_' num2str(iFile) '.xml' ]);
        indEvents = strmatch( uniqueNames{iFile}, allNames, 'exact');
        if isfield(events, 'tracktype')
            tracktype = events(indEvents(1)).tracktype;
        end
    else
        eventtrackfilename = fullfile(mffFile, fileName);
        if isfield(events, 'tracktype')
            tracktype = events(1).tracktype;
        end
        indEvents = 1:length(events);
    end

    if ~all(cellfun(@(x)isequal(x, 'boundary'), {events(indEvents).type}))
        eventtracktype = javaObject('com.egi.services.mff.api.MFFResourceType', javaMethod('valueOf', 'com.egi.services.mff.api.MFFResourceType$MFFResourceTypes', 'kMFF_RT_EventTrack'));
        if mfffactory.createResourceAtURI(eventtrackfilename, eventtracktype)
            disp('Success at creating the event file');
        else
            disp('Could not create event file');
        end

        fprintf('Exporting %d events...\n', length(events));
        eventtrackObj = mfffactory.openResourceAtURI(eventtrackfilename, eventtracktype);
        if isfield(events, 'name') && ~isempty(uniqueNames{iFile})
            eventtrackObj.setName(events(1).name);
        end
        if isfield(events, 'tracktype') && ~isempty(tracktype)
            eventtrackObj.setTrackType(tracktype);
        end

        jList = javaObject('java.util.ArrayList');
        addLatency = 0;
        multiplier = 24*60*60*srate;
        keyWarning = true;
        for iEvent = 1:length(events)
            if strcmpi(events(iEvent).type, 'boundary') % do not export boundary events
                addLatency = addLatency + events(iEvent).duration;
            else
                if ~isfield(events, 'name') || strcmpi(events(iEvent).name, uniqueNames{iFile})
                    eventObj = javaObject('com.egi.services.mff.api.Event');

                    % Get keys for event and display key codes
                    if isfield(events(iEvent), 'code')
                        eventObj.setCode(        events(iEvent).code );
                    else
                        eventObj.setCode(num2str(events(iEvent).type));
                    end
                    if isfield(events(iEvent), 'description'),  eventObj.setDescription( events(iEvent).description ); end
                    if isfield(events(iEvent), 'duration'   ),  eventObj.setDuration(    events(iEvent).duration ); end
                    if isfield(events(iEvent), 'label'),        eventObj.setLabel(       events(iEvent).label ); end
                    if isfield(events(iEvent), 'sourcedevice'), eventObj.setSourceDevice(events(iEvent).sourcedevice ); end

                    if isfield(events, 'mffkeysbackup') && ~isempty(events(iEvent).mffkeysbackup)
                        jListKeys = javaObject('java.util.ArrayList');
                        tmpmffkeysbackup = events(iEvent).mffkeysbackup;
                        tmpmffkeysbackup(tmpmffkeysbackup == 10) = ' '; % remove carriage return
                        tmpKeys = eval(tmpmffkeysbackup);
                        for iKey = 1:length(tmpKeys);
                            keyObj = javaObject('com.egi.services.mff.api.Key');
                            keyObj.setCode(tmpKeys(iKey).code);
                            keyObj.setData(tmpKeys(iKey).data);
                            keyObj.setDataType(tmpKeys(iKey).datatype);
                            keyObj.setDescription(tmpKeys(iKey).description);
                            jListKeys.add(keyObj);
                        end
                        eventObj.setKeys(jListKeys);
                    elseif keyWarning
                        disp('Warning: no MFF backup key structure found - MFF keys not exported');
                        disp('         import keys when loading the data to be able to export them');
                        keyWarning = false;
                    end

                    realLatency = (events(iEvent).latency + addLatency)/multiplier+begTime;
                    tmp = mff_encodetime(realLatency, timeZone);
                    if isfield(events(iEvent), 'begintime')
                        if ~isequal(events(iEvent).begintime, tmp) && ~isempty(events(iEvent).begintime)
                            fprintf('Note: exported event time %d differ from original one %s vs %s\n', iEvent, events(iEvent).begintime, tmp);
                        end
                    end
                    eventObj.setBeginTime(mff_encodetime(realLatency, timeZone));

                    jList.add(eventObj);
                end
            end
        end
        eventtrackObj.setEvents(jList);
        eventtrackObj.saveResource();
    end
end


text:
% TIMEF - Returns estimates and plots of mean event-related spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) changes
%           across event-related trials (epochs) of a single input time series.
%        * Uses either fixed-window, zero-padded FFTs (fastest), wavelet
%           0-padded DFTs (both Hanning-tapered), OR multitaper spectra ('mtaper').
%        * For the wavelet and FFT methods, output frequency spacing
%           is the lowest frequency ('srate'/'winsize') divided by 'padratio'.
%           NaN input values (such as returned by EVENTLOCK) are ignored.
%        * If 'alpha' is given, then bootstrap statistics are computed
%           (from a distribution of 'naccu' surrogate data trials) and
%           non-significant features of the output plots are zeroed out
%           (i.e., plotted in green).
%        * Given a 'topovec' scalp map weights vector and an 'elocs' electrode
%           location file or structure, the figure also shows a TOPOPLOT
%           image of the specified scalp map.
%
%        * Note: Left-click on subplots to view and zoom in separate windows.
% Usage:
%        >> [ersp,itc,powbase,times,freqs,erspboot,itcboot,itcphase] = ...
%                timef(data,frames,tlimits,srate,cycles,...
%                                 'key1',value1,'key2',value2, ... );
% NOTE:
%        * For more detailed information about TIMEF, >> timef details
%        * Default values may differ when called from POP_TIMEF
%
% Required inputs:
%       data        = Single-channel data vector (1,frames*ntrials) (required)
%       frames      = Frames per trial                     {def|[]: datalength}
%       tlimits     = [mintime maxtime] (ms) Epoch time limits
%                      {def|[]: from frames,srate}
%       srate       = data sampling rate (Hz)                  {def:250}
%       cycles      = If 0 -> Use FFTs (with constant window length) {0 = FFT}
%                     If >0 -> Number of cycles in each analysis wavelet
%                     If [wavecycles factor] -> wavelet cycles increase with
%                     frequency  beginning at wavecyles (0<factor<1; factor=1
%                     -> no increase, standard wavelets; factor=0 -> fixed epoch
%                     length, as in FFT.  Else, 'mtaper' -> multitaper decomp.
%
%    Optional Inter-Irial Coherence (ITC) type:
%       'type'      = ['coher'|'phasecoher'] Compute either linear coherence
%                      ('coher') or phase coherence ('phasecoher') also known
%                      as the phase coupling factor           {'phasecoher'}.
%    Optional detrending:
%       'detret'    = ['on'|'off'], Detrend data in time.               {'off'}
%       'detrep'    = ['on'|'off'], Detrend data across trials          {'off'}
%
%    Optional FFT/DFT parameters:
%       'winsize'   = If cycles==0: data subwindow length (fastest, 2^n<frames);
%                     If cycles >0: *longest* window length to use. This
%                      determines the lowest output frequency       {~frames/8}
%       'timesout'  = Number of output times (int<frames-winframes)       {200}
%       'padratio'  = FFT-length/winframes (2^k)                            {2}
%                      Multiplies the number of output frequencies by
%                      dividing their spacing. When cycles==0, frequency
%                      spacing is (low_freq/padratio).
%       'maxfreq'   = Maximum frequency (Hz) to plot (& to output if cycles>0)
%                      If cycles==0, all FFT frequencies are output.      {50}
%       'baseline'  = Spectral baseline window center end-time (in ms).    {0}
%       'powbase'   = Baseline spectrum (power, not dB) to normalize the data.
%                      {def|NaN->from data}
%
%    Optional multitaper parameters:
%       'mtaper'    = If [N W], performs multitaper decomposition.
%                      (N is the time resolution and W the frequency resolution;
%                      maximum taper number is 2NW-1). Overwrites 'winsize' and
%                      'padratio'.
%                     If [N W K], uses K Slepian tapers (if possible).
%                      Phase is calculated using standard methods.
%                      The use of mutitaper with wavelets (cycles>0) is not
%                      recommended (as multiwavelets are not implemented).
%                      Uses Matlab functions DPSS, PMTM.      {no multitaper}
%
%    Optional bootstrap parameters:
%       'alpha'     = If non-0, compute two-tailed bootstrap significance prob.
%                     level. Show non-signif. output values in green       {0}
%       'naccu'     = Number of bootstrap replications to accumulate       {200}
%       'baseboot'  = Bootstrap baseline to subtract (1 -> use 'baseline'(above)
%                                                     0 -> use whole trial) {1}
%    Optional scalp map:
%       'topovec'   = Scalp topography (map) to plot                     {none}
%       'elocs'     = Electrode location file for scalp map
%                     File should be ascii in format of  >> topoplot example
%                     May also be an EEG.chanlocs struct.
%                     {default: file named in icadefs.m}
%    Optional plotting parameters:
%       'hzdir'     = ['up'|'down'] Direction of the frequency axes; reads default
%                     from icadefs.m                                     {'up'}
%       'plotersp'  = ['on'|'off'] Plot power spectral perturbations     {'on'}
%       'plotitc'   = ['on'|'off'] Plot inter trial coherence            {'on'}
%       'plotphase' = ['on'|'off'] Plot sign of the phase in the ITC panel, i.e.
%                     green->red, pos.-phase ITC, green->blue, neg.-phase ITC {'on'}
%       'erspmax'   = [real dB] set the ERSP max. for the scale (min= -max){auto}
%       'itcmax'    = [real<=1] set the ITC maximum for the scale          {auto}
%       'title'     = Optional figure title                                {none}
%       'marktimes' = Non-0 times to mark with a dotted vertical line (ms) {none}
%       'linewidth' = Line width for 'marktimes' traces (thick=2, thin=1)  {2}
%       'pboot'     = Bootstrap power limits (e.g., from TIMEF)    {from data}
%       'rboot'     = Bootstrap ITC limits (e.g., from TIMEF)      {from data}
%       'axesfont'  = Axes text font size                                  {10}
%       'titlefont' = Title text font size                                 {8}
%       'vert'      = [times_vector] -> plot vertical dashed lines at given ms.
%       'verbose'   = ['on'|'off'] print text                              {'on'}
%
%    Outputs:
%            ersp   = Matrix (nfreqs,timesout) of log spectral diffs. from
%                     baseline (in dB).  NB: Not masked for significance.
%                     Must do this using erspboot
%            itc    = Matrix of inter-trial coherencies (nfreqs,timesout)
%                     (range: [0 1]) NB: Not masked for significance.
%                     Must do this using itcboot
%          powbase  = Baseline power spectrum (NOT in dB, used to norm. the ERSP)
%            times  = Vector of output times (sub-window centers) (in ms)
%            freqs  = Vector of frequency bin centers (in Hz)
%         erspboot  = Matrix (2,nfreqs) of [lower;upper] ERSP significance diffs
%          itcboot  = Matrix (2,nfreqs) of [lower;upper] ITC thresholds (not diffs)
%          itcphase = Matrix (nfreqs,timesout) of ITC phase (in radians)
%
% Plot description:
%   Assuming both 'plotersp' and 'plotitc' options are 'on' (= default).
%   The upper panel presents the data ERSP (Event-Related Spectral Perturbation)
%   in dB, with mean baseline spectral activity (in dB) subtracted. Use
%   "'baseline', NaN" to prevent TIMEF from removing the baseline.
%   The lower panel presents the data ITC (Inter-Trial Coherence).
%   Click on any plot axes to pop up a new window (using 'AXCOPY')
%   -- Upper left marginal panel presents the mean spectrum during the baseline
%      period (blue), and when significance is set, the significance threshold
%      at each frequency (dotted green-black trace).
%   -- The marginal panel under the ERSP image shows the maximum (green) and
%      minimum (blue) ERSP values relative to baseline power at each frequency.
%   -- The lower left marginal panel shows mean ITC across the imaged time range
%      (blue), and when significance is set, the significance threshold (dotted
%      green-black).
%   -- The marginal panel under the ITC image shows the ERP (which is produced by
%      ITC across the data spectral pass band).
%
% Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
% Known problems:
%   Significance masking currently fails for linear coherence.
%
% See also: CROSSF

% Copyright (C) 1998 Sigurd Enghoff, Scott Makeig, Arnaud Delorme,
% CNL / Salk Institute 8/1/98-8/28/01
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 10-19-98 avoided division by zero (using MIN_ABS) -sm
% 10-19-98 improved usage message and commandline info printing -sm
% 10-19-98 made valid [] values for tvec and g.elocs -sm
% 04-01-99 added missing freq in freqs and plots, fixed log scaling bug -se & -tpj
% 06-29-99 fixed frequency indexing for constant-Q -se
% 08-24-99 reworked to handle NaN input values -sm
% 12-07-99 adjusted ERPtimes to plot ERP under ITC -sm
% 12-22-99 debugged ERPtimes, added BASE_BOOT -sm
% 01-10-00 debugged BASE_BOOT=0 -sm
% 02-28-00 added NOTE on formula derivation below -sm
% 03-16-00 added AXCOPY feature -sm & tpj
% 04-16-00 added multiple marktimes loop -sm
% 04-20-00 fixed ITC cbar limits when specified in input -sm
% 07-29-00 changed frequencies displayed msg -sm
% 10-12-00 fixed bug in freqs when cycles>0 -sm
% 02-07-01 fixed inconsistency in BASE_BOOT use -sm
% 08-28-01 matlab 'key' value arguments -ad
% 08-28-01 multitaper decomposition -ad
% 01-25-02 reformated help & license -ad
% 03-08-02 debug & compare to old timef function -ad
% 03-16-02 timeout automatically adjusted if too high -ad
% 04-02-02 added 'coher' option -ad

function [P,R,mbase,times,freqs,Pboot,Rboot,Rphase,PA] = timef(X,frames,tlimits,Fs,varwin,varargin);

% Note: undocumented arg PA is output of 'phsamp','on'

%varwin,winsize,g.timesout,g.padratio,g.maxfreq,g.topovec,g.elocs,g.alpha,g.marktimes,g.powbase,g.pboot,g.rboot)

% ITC:   Normally, R = |Sum(Pxy)| / (Sum(|Pxx|)*Sum(|Pyy|)) is linear coherence.
%        But here, we consider:  Phase(Pyy) = 0 and |Pyy| = 1 -> Pxy = Pxx
%        Giving, R = |Sum(Pxx)|/Sum(|Pxx|), the inter-trial coherence (ITC)
%        Also called 'phase-locking factor' by Tallon-Baudry et al. (1996),
%        the ITC is the phase coherence between the data time series and the
%        time-locking event time series.

% Read system-wide / dir-wide constants:
icadefs

% Constants set here:
ERSP_CAXIS_LIMIT = 0;           % 0 -> use data limits; else positive value
                                % giving symmetric +/- caxis limits.
ITC_CAXIS_LIMIT  = 0;           % 0 -> use data limits; else positive value
                                % giving symmetric +/- caxis limits.

% Commandline arg defaults:
DEFAULT_EPOCH	= NaN;		% Frames per trial
DEFAULT_TIMLIM  = NaN;	                % Time range of g.frames (ms)
DEFAULT_FS	= 250;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or fixed number of cycles.
				% =0: fix window length to that determined by nwin
				% >0: set window length equal to varwin cycles
				%     Bounded above by winsize, which determines
				%     the min. freq. to be computed.
DEFAULT_OVERSMP	= 2;			% Number of times to oversample frequencies
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= '';			% Figure title
DEFAULT_ELOC    = 'chan.locs';	% Channel location file
DEFAULT_ALPHA   = NaN;			% Percentile of bins to keep
DEFAULT_MARKTIME= NaN;

% Font sizes:
AXES_FONT       = 10;           % axes text FontSize
TITLE_FONT      = 8;

if nargout>7
   Rphase = []; % initialize in case Rphase asked for, but ITC not computed
end

if (nargin < 1)
	help timef
	return
end

if ischar(X) && strcmp(X,'details')
   more on
   help timefdetails
   more off
   return
end

if (min(size(X))~=1 || length(X)<2)
	error('Data must be a row or column vector.');
end

if nargin < 2 || isempty(frames) || isnan(frames)
	frames = DEFAULT_EPOCH;
elseif (~isnumeric(frames) || length(frames)~=1 || frames~=round(frames))
	error('Value of frames must be an integer.');
elseif (frames <= 0)
	error('Value of frames must be positive.');
elseif (rem(length(X),frames) ~= 0)
	error('Length of data vector must be divisible by frames.');
end
if isnan(frames) || isempty(frames)
    frames = length(X);
end

if nargin < 3 || isempty(tlimits) || isnan(tlimits(1))
	tlimits = DEFAULT_TIMLIM;
elseif (~isnumeric(tlimits) || sum(size(tlimits))~=3)
	error('Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
	error('tlimits interval must be ascending.');
end

if (nargin < 4)
	Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) || length(Fs)~=1)
	error('Value of srate must be a number.');
elseif (Fs <= 0)
	error('Value of srate must be positive.');
end

if isempty(tlimits) || isnan(tlimits(1))
   hlim = 1000*frames/Fs;  % fit default tlimits to srate and frames
   tlimits = [0 hlim];
end

framesdiff = frames - Fs*(tlimits(2)-tlimits(1))/1000;
if abs(framesdiff) > 1
        error('Given time limits, frames and sampling rate are incompatible');
elseif framesdiff ~= 0
   	tlimits(1) = tlimits(1) - 0.5*framesdiff*1000/Fs;
   	tlimits(2) = tlimits(2) + 0.5*framesdiff*1000/Fs;
    	fprintf('Adjusted time limits slightly, to [%.1f,%.1f] ms, to match frames and srate.\n',tlimits(1),tlimits(2));
end

if (nargin < 5)
	varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) || length(varwin)>2)
	error('Value of cycles must be a number.');
elseif (varwin < 0)
	error('Value of cycles must be zero or positive.');
end

% consider structure for these arguments
% --------------------------------------
if ~isempty(varargin)
    try, g = struct(varargin{:});
    catch, error('Argument error in the {''param'', value} sequence'); end;
end
g.tlimits = tlimits;
g.frames   = frames;
g.srate   = Fs;
g.cycles  = varwin(1);
if length(varwin)>1
	g.cyclesfact = varwin(2);
else
	g.cyclesfact = 1;
end

try, g.title;      catch, g.title = DEFAULT_TITLE; end
try, g.winsize;    catch, g.winsize = max(pow2(nextpow2(g.frames)-3),4); end
try, g.pad;        catch, g.pad = max(pow2(nextpow2(g.winsize)),4); end
try, g.timesout;   catch, g.timesout = DEFAULT_NWIN; end
try, g.padratio;   catch, g.padratio = DEFAULT_OVERSMP; end
try, g.maxfreq;    catch, g.maxfreq = DEFAULT_MAXFREQ; end
try, g.topovec;    catch, g.topovec = []; end
try, g.elocs;      catch, g.elocs = DEFAULT_ELOC; end
try, g.alpha;      catch, g.alpha = DEFAULT_ALPHA; end;
try, g.marktimes;  catch, g.marktimes = DEFAULT_MARKTIME; end
try, g.powbase;    catch, g.powbase = NaN; end
try, g.pboot;      catch, g.pboot = NaN; end
try, g.rboot;      catch, g.rboot = NaN; end
try, g.plotersp;   catch, g.plotersp = 'on'; end
try, g.plotitc;    catch, g.plotitc  = 'on'; end
try, g.detrep;     catch, g.detrep = 'off'; end
try, g.detret;     catch, g.detret = 'off'; end
try, g.baseline;   catch, g.baseline = 0; end
try, g.baseboot;   catch, g.baseboot = 1; end
try, g.linewidth;  catch, g.linewidth = 2; end
try, g.naccu;      catch, g.naccu = 200; end
try, g.mtaper;     catch, g.mtaper = []; end
try, g.vert;       catch, g.vert = []; end
try, g.type;       catch, g.type = 'phasecoher'; end
try, g.phsamp;     catch, g.phsamp = 'off'; end
try, g.plotphase;  catch, g.plotphase = 'on'; end
try, g.itcmax;     catch, g.itcmax = []; end
try, g.erspmax;    catch, g.erspmax = []; end
try, g.verbose;    catch, g.verbose = 'on'; end
try, g.chaninfo;   catch, g.chaninfo = []; end
try, g.hzdir;      catch, g.hzdir = HZDIR; end; % default from icadefs
lasterr('');

% testing arguments consistency
% -----------------------------
if strcmp(g.hzdir,'up')
    g.hzdir = 'normal';
elseif strcmp(g.hzdir,'down')
    g.hzdir = 'reverse';
else
    error('unknown ''hzdir'' value - not ''up'' or ''down''');
end

switch lower(g.verbose)
    case { 'on', 'off' }, ;
    otherwise error('verbose must be either on or off');
end
if (~ischar(g.title))
	error('Title must be a string.');
end

if (~isnumeric(g.winsize) || length(g.winsize)~=1 || g.winsize~=round(g.winsize))
	error('Value of winsize must be an integer number.');
elseif (g.winsize <= 0)
	error('Value of winsize must be positive.');
elseif (g.cycles == 0 && pow2(nextpow2(g.winsize)) ~= g.winsize)
	error('Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (g.winsize > g.frames)
	error('Value of winsize must be less than frames per epoch.');
end

if (~isnumeric(g.timesout) || length(g.timesout)~=1 || g.timesout~=round(g.timesout))
	error('Value of timesout must be an integer number.');
elseif (g.timesout <= 0)
	error('Value of timesout must be positive.');
end
if (g.timesout > g.frames-g.winsize)
	g.timesout = g.frames-g.winsize;
	disp(['Value of timesout must be <= frames-winsize, timeout adjusted to ' int2str(g.timesout) ]);
end

if (~isnumeric(g.padratio) || length(g.padratio)~=1 || g.padratio~=round(g.padratio))
	error('Value of padratio must be an integer.');
elseif (g.padratio <= 0)
	error('Value of padratio must be positive.');
elseif (pow2(nextpow2(g.padratio)) ~= g.padratio)
	error('Value of padratio must be an integer power of two [1,2,4,8,16,...]');
end

if (~isnumeric(g.maxfreq) || length(g.maxfreq)~=1)
	error('Value of maxfreq must be a number.');
elseif (g.maxfreq <= 0)
	error('Value of maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
	myprintf(g.verbose,['Warning: value of maxfreq reduced to Nyquist rate' ...
		 ' (%3.2f)\n\n'], Fs/2);
	g.maxfreq = Fs/2;
end

if isempty(g.topovec)
	g.topovec = [];
	if isempty(g.elocs)
		error('Channel location file must be specified.');
	end
end
if isempty(g.elocs)
	g.elocs = DEFAULT_ELOC;
elseif (~ischar(g.elocs)) && ~isstruct(g.elocs)
	error('Channel location file must be a valid text file.');
end

if (~isnumeric(g.alpha) || length(g.alpha)~=1)
	error('timef(): Value of g.alpha must be a number.\n');
elseif (round(g.naccu*g.alpha) < 2)
	myprintf(g.verbose,'Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
    g.naccu = round(2/g.alpha);
	myprintf(g.verbose,'  Increasing the number of bootstrap iterations to %d\n',g.naccu);
end
if g.alpha>0.5 || g.alpha<=0
    error('Value of g.alpha is out of the allowed range (0.00,0.5).');
end
if ~isnan(g.alpha)
   if g.baseboot > 0
     myprintf(g.verbose,'Bootstrap analysis will use data in baseline (pre-0 centered) subwindows only.\n')
   else
     myprintf(g.verbose,'Bootstrap analysis will use data in all subwindows.\n')
   end
end
if ~isnumeric(g.vert)
    error('vertical line(s) option must be a vector');
else
	if ~isempty(g.vert)
        if min(g.vert(:)) < g.tlimits(1) || max(g.vert(:)) > g.tlimits(2)
            error('vertical line(s) time out-of-bound');
        end
	end
end

if ~isnan (g.rboot)
  if size(g.rboot) == [1,1]
    if g.cycles == 0
        g.rboot = g.rboot*ones(g.winsize*g.padratio/2);
    end
  end
end

if ~isempty(g.mtaper) % mutitaper, inspired from Bijan Pesaran matlab function
  if length(g.mtaper) < 3
        %error('mtaper argument must be [N W] or [N W K]');

    if g.mtaper(1) * g.mtaper(2) < 1
        error('mtaper 2 first arguments'' product must be higher than 1');
    end
    if length(g.mtaper) == 2
        g.mtaper(3) = floor( 2*g.mtaper(2)*g.mtaper(1) - 1);
    end
    if length(g.mtaper) == 3
        if g.mtaper(3) > 2 * g.mtaper(1) * g.mtaper(2) -1
            error('mtaper number too high (maximum (2*N*W-1))');
        end
    end
    disp(['Using ' num2str(g.mtaper(3)) ' tapers.']);
    NW = g.mtaper(1)*g.mtaper(2);   % product NW
    N  = g.mtaper(1)*g.srate;
    [e,v] = dpss(N, NW, 'calc');
    e=e(:,1:g.mtaper(3));
    g.alltapers = e;
  else
    g.alltapers = g.mtaper;
    disp('mtaper argument not [N W] or [N W K]; considering raw taper matrix');
  end

  g.winsize = size(g.alltapers, 1);
  g.pad = max(pow2(nextpow2(g.winsize)),256); % pad*nextpow

  %nfk = floor([0 g.maxfreq]./g.srate.*g.pad); % not used any more
  %g.padratio = 2*nfk(2)/g.winsize;

  g.padratio = g.pad/g.winsize;

  %compute number of frequencies
  %nf = max(256, g.pad*2^nextpow2(g.winsize+1));
  %nfk = floor([0 g.maxfreq]./g.srate.*nf);
  %freqs = linspace( 0, g.maxfreq, diff(nfk)); % this also work in the case of a FFT

end;

switch lower(g.plotphase)
    case { 'on', 'off' }, ;
    otherwise error('plotphase must be either on or off');
end
switch lower(g.plotersp)
    case { 'on', 'off' }, ;
    otherwise error('plotersp must be either on or off');
end
switch lower(g.plotitc)
    case { 'on', 'off' }, ;
    otherwise error('plotitc must be either on or off');
end
switch lower(g.detrep)
    case { 'on', 'off' }, ;
    otherwise error('detrep must be either on or off');
end
switch lower(g.detret)
    case { 'on', 'off' }, ;
    otherwise error('detret must be either on or off');
end
switch lower(g.phsamp)
    case { 'on', 'off' }, ;
    otherwise error('phsamp must be either on or off');
end
if ~isnumeric(g.linewidth)
    error('linewidth must be numeric');
end
if ~isnumeric(g.naccu)
    error('naccu must be numeric');
end
if ~isnumeric(g.baseline)
    error('baseline must be numeric');
end
switch g.baseboot
    case {0,1}, ;
    otherwise, error('baseboot must be 0 or 1');
end
switch g.type
    case { 'coher', 'phasecoher', 'phasecoher2' },;
    otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;
if isnan(g.baseline)
    g.unitpower = 'uV/Hz';
else
    g.unitpower = 'dB';
end

if (g.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
    freqs = linspace(0, g.srate/2, g.padratio*g.winsize/2+1);
    freqs = freqs(2:end);
    win = hanning(g.winsize);

    P  = zeros(g.padratio*g.winsize/2,g.timesout); % summed power
    PP = zeros(g.padratio*g.winsize/2,g.timesout); % power
    R  = zeros(g.padratio*g.winsize/2,g.timesout); % mean coherence
    RR = zeros(g.padratio*g.winsize/2,g.timesout); % (coherence)
    Pboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap power
    Rboot = zeros(g.padratio*g.winsize/2,g.naccu); % summed bootstrap coher
    Rn = zeros(1,g.timesout);
    Rbn = 0;

	switch g.type
	    case { 'coher' 'phasecoher2' },
           cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
           cumulXboot = zeros(g.padratio*g.winsize/2,g.naccu);
        case 'phasecoher'
           switch g.phsamp
             case 'on'
               cumulX = zeros(g.padratio*g.winsize/2,g.timesout);
           end
    end;

else % %%%%%%%%%%%%%%%%%% cycles>0, Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%

    freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
    dispf = find(freqs <= g.maxfreq);
    freqs = freqs(dispf);

    win = dftfilt(g.winsize,g.maxfreq/g.srate,g.cycles,g.padratio,g.cyclesfact);
    P = zeros(size(win,2),g.timesout);       % summed power
    R = zeros(size(win,2),g.timesout);       % mean coherence
    PP = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
    RR = repmat(NaN,size(win,2),g.timesout); % initialize with NaN
    Pboot = zeros(size(win,2),g.naccu);  % summed bootstrap power
    Rboot = zeros(size(win,2),g.naccu);  % summed bootstrap coher
    Rn = zeros(1,g.timesout);
    Rbn = 0;

	switch g.type
	  case { 'coher' 'phasecoher2' },
           cumulX = zeros(size(win,2),g.timesout);
           cumulXboot = zeros(size(win,2),g.naccu);
      case 'phasecoher'
           switch g.phsamp
             case 'on'
               cumulX = zeros(size(win,2),g.timesout);
           end
   end;
end

switch g.phsamp
  case 'on'
    PA = zeros(size(P,1),size(P,1),g.timesout); % NB: (freqs,freqs,times)
end                                             %       phs   amp

wintime = 1000/g.srate*(g.winsize/2); % (1000/g.srate)*(g.winsize/2);
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];
ERPtimes = [g.tlimits(1):(g.tlimits(2)-g.tlimits(1))/(g.frames-1):g.tlimits(2)+0.000001];
ERPindices = [];
for ti=times
 [tmp indx] = min(abs(ERPtimes-ti));
 ERPindices  = [ERPindices indx];
end
ERPtimes = ERPtimes(ERPindices); % subset of ERP frames on t/f window centers

if ~isempty(find(times < g.baseline))
   baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
else
   baseln = 1:length(times); % use all times as baseline
end
if ~isnan(g.alpha) && length(baseln)==0
  myprintf(g.verbose,'timef(): no window centers in baseline (times<%g) - shorten (max) window length.\n', g.baseline)
  return
elseif ~isnan(g.alpha) && g.baseboot
  myprintf(g.verbose,'   %d bootstrap windows in baseline (center times < %g).\n',...
          length(baseln), g.baseline)
end
dispf = find(freqs <= g.maxfreq);
stp = (g.frames-g.winsize)/(g.timesout-1);

myprintf(g.verbose,'Computing Event-Related Spectral Perturbation (ERSP) and\n');
switch g.type
    case 'phasecoher',  myprintf(g.verbose,'  Inter-Trial Phase Coherence (ITC) images based on %d trials\n',length(X)/g.frames);
    case 'phasecoher2', myprintf(g.verbose,'  Inter-Trial Phase Coherence 2 (ITC) images based on %d trials\n',length(X)/g.frames);
    case 'coher',       myprintf(g.verbose,'  Linear Inter-Trial Coherence (ITC) images based on %d trials\n',length(X)/g.frames);
end
myprintf(g.verbose,'  of %d frames sampled at %g Hz.\n',g.frames,g.srate);
myprintf(g.verbose,'Each trial contains samples from %d ms before to\n',g.tlimits(1));
myprintf(g.verbose,'  %.0f ms after the timelocking event.\n',g.tlimits(2));
myprintf(g.verbose,'The window size used is %d samples (%g ms) wide.\n',g.winsize,2*wintime);
myprintf(g.verbose,'The window is applied %d times at an average step\n',g.timesout);
myprintf(g.verbose,'  size of %g samples (%g ms).\n',stp,1000*stp/g.srate);
myprintf(g.verbose,'Results are oversampled %d times; the %d frequencies\n',g.padratio,length(dispf));
myprintf(g.verbose,'  displayed are from %2.1f Hz to %3.1f Hz.\n',freqs(dispf(1)),freqs(dispf(end)));
if ~isnan(g.alpha)
  myprintf(g.verbose,'Only significant values (bootstrap p<%g) will be colored;\n',g.alpha)
  myprintf(g.verbose,'  non-significant values will be plotted in green\n');
end

trials = length(X)/g.frames;
baselength = length(baseln);
myprintf(g.verbose,'\nOf %d trials total, processing trial:',trials);

% detrend over epochs (trials) if requested
% -----------------------------------------
switch g.detrep
    case 'on'
        X = reshape(X, g.frames, length(X)/g.frames);
        X = X - mean(X,2)*ones(1, length(X(:))/g.frames);
        X = X(:)';
end;

for i=1:trials
    if (rem(i,100)==0)
        myprintf(g.verbose,'\n');
    end
	if (rem(i,10) == 0)
		myprintf(g.verbose,'%d',i);
	elseif (rem(i,2) == 0)
		myprintf(g.verbose,'.');
	end

    ERP = blockave(X,g.frames); % compute the ERP trial average

    Wn = zeros(1,g.timesout);
	for j=1:g.timesout,
		tmpX = X([1:g.winsize]+floor((j-1)*stp)+(i-1)*g.frames);
                                                      % pull out data g.frames
		tmpX = tmpX - mean(tmpX); % remove the mean for that window
        switch g.detret, case 'on', tmpX = detrend(tmpX); end
		if ~any(isnan(tmpX))
		  if (g.cycles == 0) % FFT
            if ~isempty(g.mtaper)   % apply multitaper (no hanning window)
                tmpXMT = fft(g.alltapers .* ...
                             (tmpX(:) * ones(1,size(g.alltapers,2))), g.pad);
			    %tmpXMT = tmpXMT(nfk(1)+1:nfk(2),:);
			    tmpXMT = tmpXMT(2:g.padratio*g.winsize/2+1,:);
                PP(:,j) = mean(abs(tmpXMT).^2, 2);
                  % power; can also ponderate multitaper by their eigenvalues v
		        tmpX = win .* tmpX(:);
               	tmpX = fft(tmpX, g.pad);
    		    tmpX = tmpX(2:g.padratio*g.winsize/2+1);
            else
               % TF and MC (12/2006): Calculation changes made so that
                % power can be correctly calculated from ERSP.
		        tmpX = win .* tmpX(:);
               	tmpX = fft(tmpX,g.padratio*g.winsize);
                tmpX = tmpX / g.winsize;    % TF and MC (12/11/2006): normalization, divide by g.winsize
    		    tmpX = tmpX(2:g.padratio*g.winsize/2+1);
                PP(:,j) = 2/0.375*abs(tmpX).^2; % power
                % TF and MC (12/14/2006): multiply by 2 account for negative frequencies,
                % Counteract the reduction by a factor 0.375
                % that occurs as a result of cosine (Hann) tapering. Refer to Bug 446
            end;
          else % wavelet
            if ~isempty(g.mtaper)  % apply multitaper
			    tmpXMT = g.alltapers .* (tmpX(:) * ones(1,size(g.alltapers,2)));
			    tmpXMT = transpose(win) * tmpXMT;
                PP(:,j) = mean(abs(tmpXMT).^2, 2); % power
		        tmpX = transpose(win) * tmpX(:);
            else
		        tmpX = transpose(win) * tmpX(:);
                PP(:,j) = abs(tmpX).^2; % power
            end
          end

          if abs(tmpX) < eps    % If less than smallest possible machine value
                                % (i.e. if it's zero) then call it 0.
		        RR(:,j) = zeros(size(RR(:,j)));
          else
		      switch g.type
		        case { 'coher' },
		          RR(:,j) = tmpX;
                  cumulX(:,j) = cumulX(:,j)+abs(tmpX).^2;
		        case { 'phasecoher2' },
		          RR(:,j) = tmpX;
                  cumulX(:,j) = cumulX(:,j)+abs(tmpX);
		        case 'phasecoher',
		          RR(:,j) = tmpX ./ abs(tmpX); % normalized cross-spectral vector
         		  switch g.phsamp
             		  case 'on'
                	    cumulX(:,j) = cumulX(:,j)+abs(tmpX); % accumulate for PA
          		  end
              end
         end
          Wn(j) = 1;
        end

        switch g.phsamp
         case 'on' % PA (freq x freq x time)
          PA(:,:,j) = PA(:,:,j)  + (tmpX ./ abs(tmpX)) * ((PP(:,j)))';
                                           % cross-product: unit phase (column)
                                           %                times amplitude (row)
        end
	end % window

	if ~isnan(g.alpha) % save surrogate data for bootstrap analysis
        j = 1;
        goodbasewins = find(Wn==1);
        if g.baseboot % use baseline windows only
          goodbasewins = find(goodbasewins<=baselength);
        end
        ngdbasewins = length(goodbasewins);
        if ngdbasewins>1
		  while j <= g.naccu
            i=ceil(rand*ngdbasewins);
            i=goodbasewins(i);
			Pboot(:,j) = Pboot(:,j) + PP(:,i);
		    Rboot(:,j) = Rboot(:,j) + RR(:,i);
		    switch g.type
		        case 'coher',       cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX).^2;
		        case 'phasecoher2', cumulXboot(:,j) = cumulXboot(:,j)+abs(tmpX);
            end
            j = j+1;
          end
          Rbn = Rbn + 1;
	    end
	end % bootstrap

    Wn = find(Wn>0);
    if length(Wn)>0
	  P(:,Wn) = P(:,Wn) + PP(:,Wn); % add non-NaN windows
	  R(:,Wn) = R(:,Wn) + RR(:,Wn);
	  Rn(Wn) = Rn(Wn) + ones(1,length(Wn)); % count number of addends
    end
end % trial

% if coherence, perform the division
% ----------------------------------
switch g.type
 case 'coher',
  R = R ./ ( sqrt( trials*cumulX ) );
  if ~isnan(g.alpha)
	  Rboot = Rboot ./ ( sqrt( trials*cumulXboot ) );
  end
 case 'phasecoher2',
  R = R ./ ( cumulX );
  if ~isnan(g.alpha)
	  Rboot = Rboot ./ cumulXboot;
  end;
 case 'phasecoher',
  R = R ./ (ones(size(R,1),1)*Rn);
end;

switch g.phsamp
 case 'on'
  tmpcx(1,:,:) = cumulX; % allow ./ below
  for j=1:g.timesout
    PA(:,:,j) = PA(:,:,j) ./ repmat(PP(:,j)', [size(PP,1) 1]);
  end
end

if min(Rn) < 1
  myprintf(g.verbose,'timef(): No valid timef estimates for windows %s of %d.\n',...
                         int2str(find(Rn==0)),length(Rn));
  Rn(find(Rn<1))==1;
  return
end
P = P ./ (ones(size(P,1),1) * Rn);

if isnan(g.powbase)
  myprintf(g.verbose,'\nComputing the mean baseline spectrum\n');
  mbase = mean(P(:,baseln),2)';
else
  myprintf(g.verbose,'Using the input baseline spectrum\n');
  mbase = g.powbase;
end
if ~isnan( g.baseline(1) ) && ~isnan( mbase(1) )
    P = 10 * (log10(P) - repmat(log10(mbase(1:size(P,1)))',[1 g.timesout])); % convert to (10log10) dB
else
    P = 10 * log10(P);
end

Rsign = sign(imag(R));
if nargout > 7
   for lp = 1:size(R,1)
       Rphase(lp,:) = rem(angle(R(lp,:)),2*pi); % replaced obsolete phase() -sm 2/1/6
   end
   Rphase(find(Rphase>pi))  = 2*pi-Rphase(find(Rphase>pi));
   Rphase(find(Rphase<-pi)) = -2*pi-Rphase(find(Rphase<-pi));
end

R = abs(R); % convert coherence vector to magnitude

if ~isnan(g.alpha) % if bootstrap analysis included . . .
    if Rbn>0
	   i = round(g.naccu*g.alpha);
       if isnan(g.pboot)
            Pboot = Pboot / Rbn; % normalize
            if ~isnan( g.baseline )
	            Pboot = 10 * (log10(Pboot) - repmat(log10(mbase)',[1 g.naccu]));
            else
                Pboot = 10 * log10(Pboot);
            end;
            Pboot = sort(Pboot');
	        Pboot = [mean(Pboot(1:i,:)) ; mean(Pboot(g.naccu-i+1:g.naccu,:))];
       else
            Pboot = g.pboot;
       end

      if isnan(g.rboot)
	       Rboot = abs(Rboot) / Rbn;
	       Rboot = sort(Rboot');
	       Rboot = mean(Rboot(g.naccu-i+1:g.naccu,:));
      else
           Rboot = g.rboot;
      end
    else
      myprintf(g.verbose,'No valid bootstrap trials...!\n');
    end
end

switch lower(g.plotitc)
   case 'on',
       switch lower(g.plotersp),
          case 'on', ordinate1 = 0.67; ordinate2 = 0.1; height = 0.33; g.plot = 1;
          case 'off', ordinate2 = 0.1; height = 0.9; g.plot = 1;
       end;
   case 'off', ordinate1 = 0.1; height = 0.9;
       switch lower(g.plotersp),
          case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
          case 'off', g.plot = 0;
       end;
end;

if g.plot
    myprintf(g.verbose,'\nNow plotting...\n');
    set(gcf,'DefaultAxesFontSize',AXES_FONT)
    colormap(jet(256));
    pos = get(gca,'position');
    q = [pos(1) pos(2) 0 0];
    s = [pos(3) pos(4) pos(3) pos(4)];
end

switch lower(g.plotersp)
 case 'on'
    %
    %%%%%%% image the ERSP %%%%%%%%%%%%%%%%%%%%%%%%%%
    %
	h(1) = subplot('Position',[.1 ordinate1 .9 height].*s+q);

	PP = P;            % PP will be ERSP power after
	if ~isnan(g.alpha) % zero out nonsignif. power differences
		PP(find((PP > repmat(Pboot(1,:)',[1 g.timesout])) ...
                    & (PP < repmat(Pboot(2,:)',[1 g.timesout])))) = 0;
	end

    if ERSP_CAXIS_LIMIT == 0
	   ersp_caxis = [-1 1]*1.1*max(max(abs(P(dispf,:))));
    else
       ersp_caxis = ERSP_CAXIS_LIMIT*[-1 1];
    end

    if ~isnan( g.baseline )
        imagesc(times,freqs(dispf),PP(dispf,:),ersp_caxis);
    else
        imagesc(times,freqs(dispf),PP(dispf,:));
    end
    set(gca,'ydir',g.hzdir);  % make frequency ascend or descend
	if ~isempty(g.erspmax)
		caxis([-g.erspmax g.erspmax]);
	end

	hold on
	plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth); % plot time 0
    if ~isnan(g.marktimes) % plot marked time
     for mt = g.marktimes(:)'
	   plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
     end
    end
	hold off
	set(h(1),'YTickLabel',[],'YTick',[])
	set(h(1),'XTickLabel',[],'XTick',[])
	if ~isempty(g.vert)
		for index = 1:length(g.vert)
			line([g.vert(index), g.vert(index)], [min(freqs(dispf)) max(freqs(dispf))], 'linewidth', 1, 'color', 'm');
		end
	end

	h(2) = gca;
	h(3) = cbar('vert'); % ERSP colorbar axes
	set(h(2),'Position',[.1 ordinate1 .8 height].*s+q)
	set(h(3),'Position',[.95 ordinate1 .05 height].*s+q)
	title([ 'ERSP(' g.unitpower ')' ])

	E = [min(P(dispf,:));max(P(dispf,:))];
	h(4) = subplot('Position',[.1 ordinate1-0.1 .8 .1].*s+q); % plot marginal ERSP means
                                                    % below the ERSP image
	plot(times,E,[0 0],...
	     [min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3], ...
             '--m','LineWidth',g.linewidth)
	axis([min(times) max(times) ...
             min(E(1,:))-max(max(abs(E)))/3 max(E(2,:))+max(max(abs(E)))/3])

	tick = get(h(4),'YTick');
	set(h(4),'YTick',[tick(1) ; tick(end)])
	set(h(4),'YAxisLocation','right')
    set(h(4),'TickLength',[0.020 0.025]);
	xlabel('Time (ms)')
	ylabel( g.unitpower )

	E = 10 * log10(mbase(dispf));
	h(5) = subplot('Position',[0 ordinate1 .1 height].*s+q); % plot mean spectrum
                                                    % to left of ERSP image
    plot(freqs(dispf),E,'LineWidth',g.linewidth)
	if ~isnan(g.alpha)
        hold on;
		plot(freqs(dispf),Pboot(:,dispf)+[E;E],'g', 'LineWidth',g.linewidth);
		plot(freqs(dispf),Pboot(:,dispf)+[E;E],'k:','LineWidth',g.linewidth)
	end

	axis([freqs(1) freqs(max(dispf)) min(E)-max(abs(E))/3 max(E)+max(abs(E))/3])
	tick = get(h(5),'YTick');
    if (length(tick)>1)
	   set(h(5),'YTick',[tick(1) ; tick(end-1)])
    end
    set(h(5),'TickLength',[0.020 0.025]);
	set(h(5),'View',[90 90])
	xlabel('Frequency (Hz)')
	ylabel( g.unitpower )
    if strcmp(g.hzdir,'normal')
        freqdir = 'reverse';
    else
        freqdir = 'normal';
    end
    set(h(5),'xdir',freqdir);  % make frequency ascend or descend
end

switch lower(g.plotitc)
  case 'on'
    %
    %%%%%%%%%%%% Image the ITC %%%%%%%%%%%%%%%%%%
    %
	h(6) = subplot('Position',[.1 ordinate2 .9 height].*s+q); % ITC image

	RR = R; % RR is the masked ITC (R)
	if ~isnan(g.alpha)
		RR(find(RR < repmat(Rboot(1,:)',[1 g.timesout]))) = 0;
	end

    if ITC_CAXIS_LIMIT == 0
	   coh_caxis = min(max(max(R(dispf,:))),1)*[-1 1]; % 1 WAS 0.4 !
    else
       coh_caxis = ITC_CAXIS_LIMIT*[-1 1];
    end

	if exist('Rsign') && strcmp(g.plotphase, 'on')
		imagesc(times,freqs(dispf),Rsign(dispf,:).*RR(dispf,:),coh_caxis); % <---
	else
		imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % <---
	end
	if ~isempty(g.itcmax)
		caxis([-g.itcmax g.itcmax]);
	end
	tmpcaxis = caxis;
    set(gca,'ydir',g.hzdir);  % make frequency ascend or descend

	hold on
	plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
    if ~isnan(g.marktimes)
     for mt = g.marktimes(:)'
	   plot([mt mt],[0 freqs(max(dispf))],'--k','LineWidth',g.linewidth);
     end
    end
	hold off
	set(h(6),'YTickLabel',[],'YTick',[])
	set(h(6),'XTickLabel',[],'XTick',[])
	if ~isempty(g.vert)
		for index = 1:length(g.vert)
			line([g.vert(index), g.vert(index)], ...
                  [min(freqs(dispf)) max(freqs(dispf))], ...
                  'linewidth', 1, 'color', 'm');
		end
	end

	h(7) = gca;
	h(8) = cbar('vert');
	%h(9) = get(h(8),'Children');
	set(h(7),'Position',[.1 ordinate2 .8 height].*s+q)
	set(h(8),'Position',[.95 ordinate2 .05 height].*s+q)
	set(h(8),'YLim',[0 tmpcaxis(2)]);
	title('ITC')

    %
    %%%%% plot the ERP below the ITC image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
	% E = mean(R(dispf,:));

    ERPmax = max(ERP);
    ERPmin = min(ERP);
    ERPmax = ERPmax + 0.1*(ERPmax-ERPmin);
    ERPmin = ERPmin - 0.1*(ERPmax-ERPmin);
	h(10) = subplot('Position',[.1 ordinate2-0.1 .8 .1].*s+q); % ERP

    plot(ERPtimes,ERP(ERPindices),...
         [0 0],[ERPmin ERPmax],'--m','LineWidth',g.linewidth);
    hold on; plot([times(1) times(length(times))],[0 0], 'k');
    axis([min(ERPtimes) max(ERPtimes) ERPmin ERPmax]);

	tick = get(h(10),'YTick');
	set(h(10),'YTick',[tick(1) ; tick(end)])
    set(h(10),'TickLength',[0.02 0.025]);
	set(h(10),'YAxisLocation','right')
	xlabel('Time (ms)')
    ylabel('\muV')
    if (~isempty(g.topovec))
      if length(g.topovec) ~= 1, ylabel(''); end; % ICA component
    end

	E = mean(R(dispf,:)');
	h(11) = subplot('Position',[0 ordinate2 .1 height].*s+q); % plot the marginal mean
                                                    % ITC left of the ITC image
	if ~isnan(g.alpha)
		plot(freqs(dispf),E,'LineWidth',g.linewidth); hold on;
        plot(freqs(dispf),Rboot(dispf),'g', 'LineWidth',g.linewidth);
        plot(freqs(dispf),Rboot(dispf),'k:','LineWidth',g.linewidth);
        axis([freqs(1) freqs(max(dispf)) 0 max([E Rboot(dispf)])+max(E)/3])
	else
		plot(freqs(dispf),E,'LineWidth',g.linewidth)
		axis([freqs(1) freqs(max(dispf)) min(E)-max(E)/3 max(E)+max(E)/3])
	end

	tick = get(h(11),'YTick');
	set(h(11),'YTick',[tick(1) ; tick(length(tick))])
	set(h(11),'View',[90 90])
        set(h(11),'TickLength',[0.020 0.025]);
	xlabel('Frequency (Hz)')
	ylabel('ERP')
    if strcmp(g.hzdir,'normal')
        freqdir = 'reverse';
    else
        freqdir = 'normal';
    end
    set(gca,'xdir',freqdir);  % make frequency ascend or descend

    %
    %%%%%%%%%%%%%%% plot a topoplot() %%%%%%%%%%%%%%%%%%%%%%%
    %
	if (~isempty(g.topovec))
		h(12) = subplot('Position',[-.1 .43 .2 .14].*s+q);
		if length(g.topovec) == 1
		  topoplot(g.topovec,g.elocs,'electrodes','off', ...
			   'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
		else
		  topoplot(g.topovec,g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
		end
		    axis('square')
	end
end; % switch

if g.plot
	try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
    if (length(g.title) > 0)
	    axes('Position',pos,'Visible','Off');
	    h(13) = text(-.05,1.01,g.title);
	    set(h(13),'VerticalAlignment','bottom')
	    set(h(13),'HorizontalAlignment','left')
        set(h(13),'FontSize',TITLE_FONT);
    end

    axcopy(gcf);
end

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

function myprintf(verbose, varargin)
    if strcmpi(verbose, 'on')
        fprintf(varargin{:});
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [v] = nanvar(x, normalizeflag, dim, flag)

% [V] = NANVAR(X, NORMALIZEFLAG, DIM, FLAG) computes the variance,
% across all cells in x along the dimension dim, accounting for NaNs.
% Normalizeflag = 1 normalizes by N, normalizeflag = [] or 0 normalizes by N-1.
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells. If flag==1, the mean will
% be subtracted first (default behavior, but to save time on already demeaned data, it
% can be set to 0).

if nargin<2 || isempty(normalizeflag)
  normalizeflag = 0;
end

D = checkinput(x);
if nargin<3 || isempty(dim),
  dim = find(D,1,'last');
else
  % check whether the requested dim can be used
  if ~D(dim), error('data can not be concatenated across dimension %d',dim); end
end

if nargin<4,
  flag = 1;
end

if flag,
  m    = nanmean(x, dim);
  x    = cellvecadd(x, -m);
end

nx   = max(size(x));
nsmp = cellfun(@sum,        isfinite(x), repmat({dim},1,nx), 'UniformOutput', 0);
ssmp = cellfun(@nansumsq,   x,           repmat({dim},1,nx), 'UniformOutput', 0);

if normalizeflag
  N = sum(nsmp);
else
  N = sum(nsmp)-1;
end

v = sum(cell2mat(ssmp), dim)./N;

function [s] = nansumsq(x, dim)

s = nansum(x.^2, dim);


text:
function test_suite = test_regression_hier

%   Run specific demo, save values and compare the results to the expected.
%   Works for both xUnit Test Framework package by Steve Eddins and for
%   the built-in Unit Testing Framework (as of Matlab version 2013b).
%
%   See also
%     TEST_ALL, DEMO_REGRESSION_HIER
%
% Copyright (c) 2011-2012 Ville Tolvanen
% Copyright (c) 2014 Tuomas Sivula

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  % Check which package to use
  if exist('initTestSuite', 'file')
    initTestSuite_path = which('initTestSuite');
  else
    initTestSuite_path = '';
  end
  if ~isempty(initTestSuite_path) && ...
     exist([initTestSuite_path(1:end-15) 'runtests'], 'file')
    % xUnit package
    initTestSuite;
  else
    % Built-in package
    % Use all functions except the @setup
    tests = localfunctions;
    tests = tests(~cellfun(@(x)strcmp(func2str(x),'setup'),tests));
    test_suite = functiontests(tests);
  end
end

% -------------
%     Tests
% -------------

function testRunDemo(testCase)
  % Run the correspondin demo and save the values. Note this test has to
  % be run at lest once before the other test may succeed.
  rundemo(getName())
end

function testPredictionMissingData(testCase)
  verifyVarsEqual(testCase, getName(), {'Eff'}, {@mean}, ...
    'RelTolElement', 0.01, 'RelTolRange', 0.005)
end

% ------------------------
%     Helper functions
% ------------------------

function testCase = setup
  % Helper function to suply empty array into variable testCase as an
  % argument for each test function, if using xUnit package. Not to be
  % used with built-in test framework.
  testCase = [];
end

function name = getName
  % Helperfunction that returns the name of the demo, e.g. 'binomial1'.
  name = mfilename;
  name = name(6:end);
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function essential_distMinAnglePair_test
resetRands(3)
flagDegenerateCase=true;
k=2;

e3=[0;0;1];
Q1=rot_randn([],[],2);
if flagDegenerateCase
    Q1b=[Q1(:,:,1);Q1(:,:,2)];
    Q2b=essential_randomVerticalMotion(Q1b);
    Q2=cat(3,Q2b(1:3,:),Q2b(4:6,:));
else
    Q2=rot_randn([],[],2);
end
Rzt=@(t) rot(t*e3);

Q21tk=@(t,k) Rzt(t)*essential_flipAmbiguity_R1(Q2(:,:,1),k);
Q22tk=@(t,k) Rzt(t)*essential_flipAmbiguity_R2(Q2(:,:,2),k);

figure(1)
[tMin,fMin,tBreak1,tBreak2,Q2Flip]=essential_distMinAnglePair([Q1(:,:,1);Q1(:,:,2)],[Q2(:,:,1);Q2(:,:,2)],k);
tMin=modAngle(tMin);
ft=@(t) (rot_dist(Q1(:,:,1),Q21tk(t,k))^2+rot_dist(Q1(:,:,2),Q22tk(t,k))^2);
dft=@(t) 2*e3'*(Q1(:,:,1)*logrot(Q1(:,:,1)'*Q21tk(t,k))+Q1(:,:,2)*logrot(Q1(:,:,2)'*Q22tk(t,k)));
check_der(ft,dft,'angle')
hold on
plot(tBreak1,ft(tBreak1),'r+')
plot(tBreak2,ft(tBreak2),'g+')

plot(tMin,fMin,'kx','MarkerSize',20)

hold off

author: NicolasBoumal

repo: manopt