text:
function [h] = dtiComputeFiberDensity(h, fiberGroupNum, endptFlag, fgCountFlag)
%
%   [handles] = dtiComputeFiberDensity(handles, [fiberGroupNum=selectFGs], [endPointFlag=false], [fgCountFlag=false])
%
% Calculate how many fibers pass through each of the voxels.  Normalize
% this number to 1.  Return the values as a 3D image attached to the
% handles structure (image is called 'fiber density').
%
%
% RETURNS:
%  handles, with a new background image representing the fiber density.
%
% HISTORY:
% 2003.12.18 RFD (bob@white.stanford.edu) wrote it.
% 2006.06.09 RFD: removed fiber density normalization in here. This allows
% dtiAddBackground to do the normalization and save the original fiber
% counts so that we don't lose these values.

normalise = false; %TODO: expose normalization as an option

if(isempty(h.fiberGroups))
    error('No fiber groups.');
end
if(~exist('fiberGroupNum','var') ||isempty(fiberGroupNum))
    fiberGroupNum = dtiSelectFGs(h, 'Select fiber group(s) for density calculation');
end
if(any(fiberGroupNum>length(h.fiberGroups)))
    error('Invalid fiber group.');
end
if(~exist('endptFlag','var') || isempty(endptFlag))
    endptFlag = false;
end
if(~exist('fgCountFlag','var') || isempty(fgCountFlag))
    fgCountFlag = false;
end

bb = dtiGet(h,'boundingBox');
imSize = diff(bb)+1;
%imSize = size(h.dt6);imSize = imSize(1:3);

%gd = dtiGet(h, 'acpcGrid', h.vec.mat, h.vec.mmPerVoxel, size(h.vec.img(:,:,:,1)));
% We need to get the T1-space image coords from the current anatomy image.
%imgCoords = mrAnatXformCoords(inv(dtiGet(h,'acpcXform')), [gd.X(:) gd.Y(:) gd.Z(:)]);
%xform = inv(h.xformToAcpc);
% TO DO: this should be based on the actual fiber step size, rather than
% assuming that it's 1mm.
mmPerVoxel = [1 1 1];
ac = mrAnatXformCoords(h.xformToAcpc,[0 0 0]);
%xformImgToAcpc = h.xformToAcpc/diag([h.mmPerVoxel 1]);
xformImgToAcpc = diag([mmPerVoxel 1]);
xformImgToAcpc(1:3,4) = ac';

fdImg = dtiComputeFiberDensityNoGUI(h.fiberGroups, xformImgToAcpc, imSize, normalise, fiberGroupNum, endptFlag, fgCountFlag);

%uniqueFiberInds = unique(fiberInds);
% Calculate how many fibers are associated with each unique fiber index.
% HOW TO DO THIS WITHOUT A LOOP?
%for(ii=uniqueFiberInds')
%    fdImg(ii) = sum(fiberInds==ii);
%end
% USE  [a, b]=hist(uniqueFiberInds)

imName = 'fiber density';
[h,bgNum] = dtiAddBackgroundImage(h, fdImg, imName, mmPerVoxel, xformImgToAcpc);
h = dtiSet(h,'curOverlayNum',bgNum);
return;

repo: vistasoft


text:
function schema
% Schema for the MOVIE package.

%   Copyright 2003 The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2005/01/07 18:00:26 $

schema.package('mmovie');

author: vistalab

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [net, info] = cnn_train_adagrad(net, imdb, getBatch, varargin)
% CNN_TRAIN   Demonstrates training a CNN
%    CNN_TRAIN() is an example learner implementing stochastic gradient
%    descent with momentum to train a CNN for image classification.
%    It can be used with different datasets by providing a suitable
%    getBatch function.

opts.train = [] ;
opts.val = [] ;
opts.numEpochs = 300 ;
opts.batchSize = 256 ;
opts.useGpu = false ;
opts.learningRate = 0.001 ;
opts.continue = false ;
opts.expDir = fullfile('data','exp') ;
opts.conserveMemory = false ;
opts.sync = true ;
opts.prefetch = false ;
opts.weightDecay = 0.0005 ;
opts.errorType = 'multiclass' ;
opts.plotDiagnostics = false ;
opts.delta = 1e-8;
opts.display = 1;
opts.snapshot = 1;
opts.test_interval = 1;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

for i=1:numel(net.layers)
  if ~strcmp(net.layers{i}.type,'conv'), continue; end
  net.layers{i}.filtersMomentum = zeros(size(net.layers{i}.filters), ...
    class(net.layers{i}.filters)) ;
  net.layers{i}.biasesMomentum = zeros(size(net.layers{i}.biases), ...
    class(net.layers{i}.biases)) ; %#ok<*ZEROLIKE>
  if ~isfield(net.layers{i}, 'filtersLearningRate')
    net.layers{i}.filtersLearningRate = 1 ;
  end
  if ~isfield(net.layers{i}, 'biasesLearningRate')
    net.layers{i}.biasesLearningRate = 1 ;
  end
  if ~isfield(net.layers{i}, 'filtersWeightDecay')
    net.layers{i}.filtersWeightDecay = 1 ;
  end
  if ~isfield(net.layers{i}, 'biasesWeightDecay')
    net.layers{i}.biasesWeightDecay = 1 ;
  end
end

if opts.useGpu
  net = vl_simplenn_move(net, 'gpu') ;
  for i=1:numel(net.layers)
    if ~strcmp(net.layers{i}.type,'conv'), continue; end
    net.layers{i}.filtersMomentum = gpuArray(net.layers{i}.filtersMomentum) ;
    net.layers{i}.biasesMomentum = gpuArray(net.layers{i}.biasesMomentum) ;
  end
end

G_f = cell(numel(net.layers), 1);
G_b = cell(numel(net.layers), 1);

for l=1:numel(net.layers)

  if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

  G_f{l} = zeros(size(net.layers{l}.filters), 'single');
  G_b{l} = zeros(size(net.layers{l}.biases), 'single');

end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

rng(0) ;

if opts.useGpu
  one = gpuArray(single(1)) ;
else
  one = single(1) ;
end

info.train.objective = [] ;
info.train.error = [] ;
info.train.topFiveError = [] ;
info.train.speed = [] ;
info.val.objective = [] ;
info.val.error = [] ;
info.val.topFiveError = [] ;
info.val.speed = [] ;

lr = opts.learningRate ;
res = [] ;
for epoch=1:opts.numEpochs

  % fast-forward to where we stopped
  modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
  modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;
  if opts.continue
    if exist(modelPath(epoch),'file')
      if epoch == opts.numEpochs
        load(modelPath(epoch), 'net', 'info') ;
      end
      continue ;
    end
    if epoch > 1
      fprintf('resuming by loading epoch %d\n', epoch-1) ;
      load(modelPath(epoch-1), 'net', 'info') ;
    end
  end

  train = opts.train(randperm(numel(opts.train))) ;
  val = opts.val ;

  info.train.objective(end+1) = 0 ;
  info.train.error(end+1) = 0 ;
  info.train.topFiveError(end+1) = 0 ;
  info.train.speed(end+1) = 0 ;
  info.val.objective(end+1) = 0 ;
  info.val.error(end+1) = 0 ;
  info.val.topFiveError(end+1) = 0 ;
  info.val.speed(end+1) = 0 ;

  for t=1:opts.batchSize:numel(train)
    % get next image batch and labels
    batch = train(t:min(t+opts.batchSize-1, numel(train))) ;
    batch_time = tic ;
    fprintf('training: epoch %02d: processing batch %3d of %3d ...', epoch, ...
            fix(t/opts.batchSize)+1, ceil(numel(train)/opts.batchSize)) ;
    [im, labels] = getBatch(imdb, batch) ;
    if opts.prefetch
      nextBatch = train(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(train))) ;
      getBatch(imdb, nextBatch) ;
    end
    if opts.useGpu
      im = gpuArray(im) ;
    end

    % backprop
    net.layers{end}.class = labels ;
    res = vl_simplenn(net, im, one, res, ...
      'conserveMemory', opts.conserveMemory, ...
      'sync', opts.sync) ;

    % gradient step
    for l=1:numel(net.layers)
      if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

      g_f = (net.layers{l}.filtersLearningRate) * ...
            (opts.weightDecay * net.layers{l}.filtersWeightDecay) * net.layers{l}.filters + ...
            (net.layers{l}.filtersLearningRate) / numel(batch) * res(l).dzdw{1};
      g_b = (net.layers{l}.biasesLearningRate) * ...
            (opts.weightDecay * net.layers{l}.biasesWeightDecay) * net.layers{l}.biases + ...
            (net.layers{l}.biasesLearningRate) / numel(batch) * res(l).dzdw{2};

      G_f{l} = G_f{l} + g_f .^ 2;
      G_b{l} = G_b{l} + g_b .^ 2;

      net.layers{l}.filters = net.layers{l}.filters - lr ./ (opts.delta + sqrt(G_f{l})) .* g_f;
      net.layers{l}.biases  = net.layers{l}.biases - lr ./ (opts.delta + sqrt(G_b{l})) .* g_b;
    end

    % print information
    batch_time = toc(batch_time) ;
    speed = numel(batch)/batch_time ;
    info.train = updateError(opts, info.train, net, res, batch_time) ;

    fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
    n = t + numel(batch) - 1 ;
    switch opts.errorType
      case 'multiclass'
        fprintf(' err %.1f err5 %.1f', ...
          info.train.error(end)/n*100, info.train.topFiveError(end)/n*100) ;
        fprintf('\n') ;
      case 'binary'
        fprintf(' err %.1f', ...
          info.train.error(end)/n*100) ;
        fprintf('\n') ;
      case 'euclideanloss'
        fprintf(' err %.1f', info.train.error(end) / n);
        fprintf('\n') ;
    end

    % debug info
    if opts.plotDiagnostics
      figure(2) ; vl_simplenn_diagnose(net,res) ; drawnow ;
    end
  end % next batch

  % evaluation on validation set
  if epoch == 1 || rem(epoch, opts.test_interval) == 0 || epoch == opts.numEpochs
  for t=1:opts.batchSize:numel(val)
    batch_time = tic ;
    batch = val(t:min(t+opts.batchSize-1, numel(val))) ;
    fprintf('validation: epoch %02d: processing batch %3d of %3d ...', epoch, ...
            fix(t/opts.batchSize)+1, ceil(numel(val)/opts.batchSize)) ;
    [im, labels] = getBatch(imdb, batch) ;
    if opts.prefetch
      nextBatch = val(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(val))) ;
      getBatch(imdb, nextBatch) ;
    end
    if opts.useGpu
      im = gpuArray(im) ;
    end

    net.layers{end}.class = labels ;
    res = vl_simplenn(net, im, [], res, ...
      'disableDropout', true, ...
      'conserveMemory', opts.conserveMemory, ...
      'sync', opts.sync) ;

    % print information
    batch_time = toc(batch_time) ;
    speed = numel(batch)/batch_time ;
    info.val = updateError(opts, info.val, net, res, batch_time) ;

    fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
    n = t + numel(batch) - 1 ;
    switch opts.errorType
      case 'multiclass'
        fprintf(' err %.1f err5 %.1f', ...
          info.val.error(end)/n*100, info.val.topFiveError(end)/n*100) ;
        fprintf('\n') ;
      case 'binary'
        fprintf(' err %.1f', ...
          info.val.error(end)/n*100) ;
        fprintf('\n') ;
      case 'euclideanloss'
        fprintf(' err %.1f', info.val.error(end) / n);
        fprintf('\n') ;
    end
  end
  end

  % save
  info.train.objective(end) = info.train.objective(end) / numel(train) ;
  info.train.error(end) = info.train.error(end) / numel(train)  ;
  info.train.topFiveError(end) = info.train.topFiveError(end) / numel(train) ;
  info.train.speed(end) = numel(train) / info.train.speed(end) ;
  info.val.objective(end) = info.val.objective(end) / numel(val) ;
  info.val.error(end) = info.val.error(end) / numel(val) ;
  info.val.topFiveError(end) = info.val.topFiveError(end) / numel(val) ;
  info.val.speed(end) = numel(val) / info.val.speed(end) ;
  if epoch == 1 || rem(epoch, opts.snapshot) == 0 || epoch == opts.numEpochs
  save(modelPath(epoch), 'net', 'info') ;
  end

  if epoch == 1 || rem(epoch, opts.display) == 0 || epoch == opts.numEpochs
  figure(1) ; clf ;
  subplot(1,2,1) ;
  semilogy(1:epoch, info.train.objective, 'k') ; hold on ;
  semilogy([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.objective([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
  xlabel('training epoch') ; ylabel('energy') ;
  grid on ;
  h=legend('train', 'val') ;
  set(h,'color','none');
  title('objective') ;
  subplot(1,2,2) ;
  switch opts.errorType
    case 'multiclass'
      plot(1:epoch, info.train.error, 'k') ; hold on ;
      plot(1:epoch, info.train.topFiveError, 'k--') ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.topFiveError([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b--') ;
      h=legend('train','train-5','val','val-5') ;
    case 'binary'
      plot(1:epoch, info.train.error, 'k') ; hold on ;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      h=legend('train','val') ;
    case 'euclideanloss'
      plot(1 : epoch, info.train.error, 'k'); hold on;
      plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
      h = legend('train', 'val') ;
  end
  grid on ;
  xlabel('training epoch') ; ylabel('error') ;
  set(h,'color','none') ;
  title('error') ;
  drawnow ;
  print(1, modelFigPath, '-dpdf') ;
  end
end

% -------------------------------------------------------------------------
function info = updateError(opts, info, net, res, speed)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
sz = size(predictions) ;
n = prod(sz(1:2)) ;

labels = net.layers{end}.class ;
info.objective(end) = info.objective(end) + sum(double(gather(res(end).x))) ;
info.speed(end) = info.speed(end) + speed ;
switch opts.errorType
  case 'multiclass'
    [~,predictions] = sort(predictions, 3, 'descend') ;
    error = ~bsxfun(@eq, predictions, reshape(labels, 1, 1, 1, [])) ;
    info.error(end) = info.error(end) +....
      sum(sum(sum(error(:,:,1,:))))/n ;
    info.topFiveError(end) = info.topFiveError(end) + ...
      sum(sum(sum(min(error(:,:,1:5,:),[],3))))/n ;
  case 'binary'
    error = bsxfun(@times, predictions, labels) < 0 ;
    info.error(end) = info.error(end) + sum(error(:))/n ;
  case 'euclideanloss'
    error = euclideanloss(sigmoid(predictions), labels);
    info.error(end) = info.error(end) + error;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef CLIA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation>
% Evolutionary algorithm with cascade clustering and reference point incremental learning

%------------------------------- Reference --------------------------------
% H. Ge, M. Zhao, L. Sun, Z. Wang, G. Tan, Q. Zhang, and C. L. P. Chen, A
% many-objective evolutionary algorithm with two interacting processes:
% Cascade clustering and reference point incremental learning, IEEE
% Transactions on Evolutionary Computation, 2019, 23(4): 572-586.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            global stable_threshold delta crowding_pick_flag;
            [stable_threshold, delta] = Algorithm.ParameterSet([0, 0, 0], 2 * Problem.N);
            [Z, P, A, S, SVM] = initialize(Problem);
            while Algorithm.NotTerminated(S)
                MatingPool = TournamentSelection(2, Problem.N, sum(max(0, P.cons), 2));
                Offspring  = OperatorGA(Problem,P(MatingPool));
                A = update_archive(A, [P, Offspring], Z, ceil(0.33 * Problem.M * Problem.N), Problem);
                % SELECTION OF INDIVIDUALS
                [P, ICA, ICN] = cascade_cluster([P, Offspring], Z, 'PDM', Problem.N, Problem.FE < Problem.maxFE);
                % ADAPTATION OF REFERENCE VECTORS
                [Z, SVM] = incremental_learn(Z, ICA, ICN, A, SVM, Problem);
                if Problem.FE >= Problem.maxFE && crowding_pick_flag
                    S = crowding_pick(update_archive(A, P, Z, [], Problem), Problem.N, 'precise');
                else
                    S = P;
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% waveform data parameters
%clear all
close all
%startup
make_figures = true;

% Geographical coordinates
% coords = ...
%     [28.574013,	-80.57236;
%      28.574182,	-80.57241;
%      28.573894,	-80.572352;
%      28.574004,	-80.572561];
coords = ...
    [28.5740171611,	-80.5723755;
     28.5742216111,	-80.5724168556;
     28.5738811889,	-80.572295775;
     28.5740064361,	-80.57256045];

% SLC40 - SpaceX launch complex
% source.lat = 28.562106;
% source.lon = -80.57718;
source.lon = -80.57719; %-80.57718;
source.lat = 28.56195; %28.562106;

easting = [461.091 473.067 447.133 465.217 465.217 465.217];
northing = [1343.9 1306.23 1320.01 1321.26 1321.26 1321.26];
% Wind tower data - could read this from Excel instead
% get Excel file from
relativeHumidity = 89.5; % percent from NASA weather tower data
temperatureF = 80.65; % 80 Fahrenheit according to weather tower data from NASA
wind_direction_from = 144; % degrees - this is the direction FROM according to NASA, see Lisa Huddleston email of Oct 10th
wind_speed_knots = 5.5; % knots

% now call
rocket_airwave_event_analysis

author: geoscience-community-codes

repo: GISMO


text:
data = dlmread('iris.dat'); % data file
%Dividing it to testing and learning sets to this example.
datalearn=[data(1:25,:);data(51:75,:);data(101:125,:)];
datatest=[data(26:50,:);data(76:100,:);data(126:150,:)];
v=[1:4];
c=5;
measure = 1; % Used measure (see classifier.m)
p = [0.1:0.25:4]; % p parameter range
m = [0.25:0.25:8]; % m parameter range
pl=1; % Do we use plotting w.r.t. parameters p and m or not pl=0 no plotting pl=1 plotting.

[Classification_accuracy,p1,m1,classes]=simclass1(datalearn, datatest,v,c, measure, p, m,pl);
Classification_accuracy

disp('Or if you just want to classify with one parameter values i.e. p=1,m=1')
[Classification_accuracy,p1,m1,classes]=simclass1(datalearn, datatest,v,c, measure, 1, 1,0);
Classification_accuracy

author: Sable

repo: mcbench-benchmarks


text:
function [bodyPlotted] = ma_plotFrame(hFig, mAxes, maData, stateLog, time, prevBody, orbitToPlot, celBodyData, handles)
%ma_plotFrame Summary of this function goes here
%   Detailed explanation goes here
    global GLOBAL_VideoWriter GLOBAL_isRecording;
    persistent usrScMarker hCBodySurf hCBodySurfXform preCbRotateAngle hCbGrdStn hOtherScPts hOtherScOrbits hCbGrdStnText hAnnote hChildBodies hChildBodyOrbits hSunBody hLight prevFrameTime curFrameTime;

    if(~isempty(usrScMarker) && ishandle(usrScMarker))
        delete(usrScMarker);
        usrScMarker = []; %#ok<NASGU>
    end

    [orbitToPlot, stateLog] = ma_getOrbitToPlot(stateLog, time);
    stateLogRow = getStateLogRowAtTime(stateLog, time);

    if(isempty(preCbRotateAngle))
        preCbRotateAngle = 0;
    end

    if(isempty(prevFrameTime) && isempty(curFrameTime))
        prevFrameTime = clock;
        curFrameTime = prevFrameTime;
    else
        curFrameTime = clock;
    end

    bodyPlotted = stateLogRow(8);
    if(bodyPlotted ~= prevBody)
        unfreezeColors(mAxes);
        cla(mAxes,'reset');
        [hCBodySurf] = plotStateLog(stateLog, handles, false, false, false, false, orbitToPlot, [], maData, celBodyData, mAxes, hFig);
        freezeColors(mAxes);
        material(hCBodySurf,'dull');
        set(hCBodySurf,'AmbientStrength',0.1);
        hCBodySurfXform = hgtransform(mAxes);
        set(hCBodySurf,'Parent',hCBodySurfXform);

        preCbRotateAngle = 0;

        if(~isnan(prevBody))
            set(handles.cameraSrcCombo,'Value',1);
            setSrcComboUserData(handles.cameraSrcCombo, handles, false);
            set(handles.cameraTgtCombo,'Value',2);
            setTgtComboUserData(handles.cameraTgtCombo, handles, false);
        end
        setCamSrcTgtComboBoxes(handles, time);

        axisWasReset = true;
    else
        axisWasReset = false;
    end

    bodyInfo = celBodyData.getBodyInfoById(bodyPlotted);
    if(not(isempty(hCBodySurfXform)) && isvalid(hCBodySurfXform))
        M = getBodyXformMatrix(time, bodyInfo, bodyInfo.getBodyCenteredInertialFrame());
        set(hCBodySurfXform,'Matrix',M);
    end

    if(axisWasReset && ~strcmpi(get(handles.cameraTypeCombo,'UserData'),'Inertially Fixed'))
%         axisLimMult = 15;
        axisLimMult = get(handles.axesScaleFactorText,'UserData');
    else
        axisLimMult = 1;
    end

    lims = zeros(0,2);
    lims(1,:) = get(mAxes,'XLim');
    lims(2,:) = get(mAxes,'YLim');
    lims(3,:) = get(mAxes,'ZLim');
    maxMax = axisLimMult * max(max(lims));
    maxMaxArr = [-maxMax maxMax -maxMax maxMax -maxMax maxMax];
    cLim = get(mAxes,'CLim');
    maxMaxArr = [maxMaxArr, cLim];
    axis(maxMaxArr);

    bodyID = stateLogRow(1,8);
    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);

    bgColors = get(handles.backgroundColorCombo,'UserData');
    showGrid = get(handles.showGridCheckbox,'UserData');
    if(strcmpi(showGrid,'on'))
        set(mAxes,'Color',bgColors(1),'XColor',bgColors(2),'YColor',bgColors(3),'ZColor',bgColors(4));
    else
        set(mAxes,'Color',bgColors(1),'XColor',bgColors(1),'YColor',bgColors(1),'ZColor',bgColors(1));
    end
    grid(mAxes,showGrid);

	switch(get(handles.cameraTypeCombo,'UserData'))
        case 'Inertially Fixed'
            azView = get(handles.azOffsetText,'UserData');
            elView = get(handles.elOffsetText,'UserData');

            view(mAxes, [azView,elView]);
        case 'Spacecraft Fixed'
            camSrcStr = get(handles.cameraSrcCombo,'UserData');
            camTgtStr = get(handles.cameraTgtCombo,'UserData');

            [typeStrSrc, ~, idStrSrc] = parseSrcTgtStrings(camSrcStr);
            [typeStrTgt, ~, idStrTgt] = parseSrcTgtStrings(camTgtStr);

            rPosCam = getPositOfObject(maData, time, stateLogRow, typeStrSrc, idStrSrc, bodyInfo, celBodyData);
            rPosTgt = getPositOfObject(maData, time, stateLogRow, typeStrTgt, idStrTgt, bodyInfo, celBodyData);

            camViewVect = rPosTgt - rPosCam;

            [cAz,cEl,cR] = cart2sph(camViewVect(1),camViewVect(2),camViewVect(3));

            azAdd = deg2rad(get(handles.azOffsetText,'UserData'));
            elAdd = deg2rad(get(handles.elOffsetText,'UserData'));
            rngAdd = get(handles.rngOffsetText,'UserData');

            cAz = cAz + azAdd;
            cEl = cEl - elAdd;
            cR = cR + rngAdd;

            [cX,xY,cZ] = sph2cart(cAz,cEl,cR);
            camViewVect = [cX,xY,cZ];

            rPosCam = rPosTgt - camViewVect;

            campos(mAxes,rPosCam); camtarget(mAxes,rPosTgt);
            camva(mAxes,get(handles.fieldOfViewText,'UserData'));
    end

    topLevelBodyInfo = getTopLevelCentralBody(celBodyData);
    rVectSun = getAbsPositBetweenSpacecraftAndBody(time, [0,0,0]', bodyInfo, topLevelBodyInfo, celBodyData);
    if(norm(rVectSun) < 0.01)
        rVectSun = [0 0 0.01];
    end

    rVectSun = normVector(rVectSun);
    if(isempty(hLight) || ~ishghandle(hLight))
        hLight = light('Position',rVectSun,'Parent',mAxes,'style','infinite');
    else
        set(hLight,'Position',rVectSun);
    end

    hold(mAxes,'on');
    markerColor = get(handles.spacecraftMarkerColorCombo,'UserData');
    lineSpec = [get(handles.spacecraftMarkerCombo,'UserData'),markerColor];
    markerSize = get(handles.scMarkerSizeCombo,'UserData');
    usrScMarker = plot3(mAxes, stateLogRow(2),stateLogRow(3),stateLogRow(4),lineSpec,'MarkerSize',markerSize,'MarkerFaceColor',markerColor);
    hold(mAxes,'off');

%     cbSpinAngle = getBodySpinAngle(bodyInfo, time);
%     dAlpha = rad2deg(cbSpinAngle) - preCbRotateAngle;
%     rotate(hCBodySurf,[0,0,1],dAlpha,[0,0,0]);
%     set(hCBodySurf,'NormalMode','auto');
%     preCbRotateAngle = rad2deg(cbSpinAngle);

    if(~isempty(hCbGrdStn(ishghandle(hCbGrdStn))))
        delete(hCbGrdStn(ishghandle(hCbGrdStn))); hCbGrdStn = [];
    end

    if(get(handles.showGroundStationsCheckbox,'UserData')==1)
        stations = maData.spacecraft.stations;
        bdyStnFHdl = @(stnStruct) plotLocOfGrndStation(mAxes, stnStruct, bodyID, time, [0,0,0]', bodyInfo);

        [hCbGrdStn] = cellfun(bdyStnFHdl, stations,'UniformOutput',false);
        hCbGrdStn = [hCbGrdStn{:}];
    end

    delete(hOtherScPts(ishghandle(hOtherScPts)));
    delete(hOtherScOrbits(ishghandle(hOtherScOrbits)));
    if(get(handles.showOtherScCheckbox,'UserData')==1)
        otherSCs = maData.spacecraft.otherSC; %hOtherScOrbits

        otherScFHdl = @(oScStruct) plotLocOfOtherSC(mAxes, oScStruct, bodyID, time, [0,0,0]', bodyInfo, true, true);

        [hOtherScPts, hOtherScOrbits] = cellfun(otherScFHdl, otherSCs,'UniformOutput',false);
        hOtherScPts = [hOtherScPts{:}];
        hOtherScOrbits = [hOtherScOrbits{:}];
    end

    delete(hChildBodies(ishghandle(hChildBodies))); hChildBodies = [];
    showChildBodies = get(handles.showChildBodiesCheckbox,'UserData');
    showChildBodyOrbits = get(handles.showChildBodyOrbitsCheckbox,'UserData');
    if(showChildBodies == 1 || showChildBodyOrbits == 1)
        [children, ~] = getChildrenOfParentInfo(celBodyData, bodyInfo.name);
        for(i=1:length(children)) %#ok<*NO4LP>
            child = children{i};

            if(showChildBodies==1)
                [rVect, ~] = getStateAtTime(child, time, bodyInfo.gm);

                useR = child.radius;
                mColor = colorFromColorMap(child.bodycolor);

                [cX, cY, cZ] = sphere(30);
                hold(mAxes,'on')
                hChildBodies(end+1) = plot3(rVect(1), rVect(2), rVect(3),'Marker','o','MarkerEdgeColor',mColor,'MarkerFaceColor',mColor,'MarkerSize',3); %#ok<AGROW>
                hold(mAxes,'off')
                hChildBodies(end+1) = surface(rVect(1) + useR*cX, rVect(2) + useR*cY, rVect(3) + useR*cZ); %#ok<AGROW>
                colormap(mAxes, child.bodycolor)
                freezeColors(mAxes);
            end

            if(showChildBodyOrbits == 1)
                if(axisWasReset || length(hChildBodyOrbits) < length(children))
                    if(length(hChildBodyOrbits) < length(children))
                        if(i==1)
                            delete(hChildBodyOrbits(ishghandle(hChildBodyOrbits)));
                            hChildBodyOrbits = [];
                        end
                    end
                    hChildBodyOrbits(end+1) = plotBodyOrbit(child, '', bodyInfo.gm, false, mAxes, 0.5, '--'); %#ok<AGROW>
                end
            else
                delete(hChildBodyOrbits(ishghandle(hChildBodyOrbits))); hChildBodyOrbits = [];
            end
        end
    elseif(showChildBodyOrbits ~= 1)
        delete(hChildBodyOrbits(ishghandle(hChildBodyOrbits))); hChildBodyOrbits = [];
    end

    aShowUT = get(handles.annoteShowUTCheckbox,'UserData');
	aShowOrbit = get(handles.annoteShowOrbitCheckbox,'UserData');
    aShowLatLong = get(handles.annoteShowLatLongCheckbox,'UserData');
    aShowMass = get(handles.annoteShowMassCheckbox,'UserData');
    aShowFps = get(handles.showFpsCheckbox,'UserData');
    aShowEvent = get(handles.showEventNameCheckbox,'UserData');

    if(aShowUT==1 || aShowOrbit==1 || aShowLatLong==1 || aShowMass==1 || aShowFps==1 || aShowEvent==1)
        annonStr = cell(0,1);

        if(aShowUT==1)
            [year, day, hour, minute, sec] = convertSec2YearDayHrMnSec(time);
            dateStr = formDateStr(year, day, hour, minute, sec);
            annonStr{end+1} = dateStr;
        end

        if(aShowOrbit==1)
            rVect = stateLogRow(2:4);
            vVect = stateLogRow(5:7);
            gmu = bodyInfo.gm;
            [sma, ecc, inc, raan, arg, tru] = getKeplerFromState(rVect,vVect,gmu);
            oStr = sprintf('Orbit\n   Orbiting=%s\n   SMA=%g km\n   Ecc=%g\n   Inc=%g deg\n   RAAN=%g deg\n   Arg=%g deg\n   Tru=%g deg', bodyInfo.name, sma, ecc, rad2deg(inc), rad2deg(raan), rad2deg(arg), rad2deg(tru));
            annonStr{end+1} = oStr;
        end

        if(aShowEvent==1)
            eventName = maData.script{stateLogRow(13)}.name;
            annonStr{end+1} = sprintf('Event: %s', eventName);
        end

        if(aShowLatLong==1)
            rVect = stateLogRow(2:4)';
            vVect = stateLogRow(5:7)';
            [lat, long, alt] = getLatLongAltFromInertialVect(time, rVect, bodyInfo, vVect);
            ffStr = sprintf('Rotating Frame Position\n   Lat.=%g deg\n   Long.=%g deg\n   Alt.=%g km', rad2deg(lat), rad2deg(long), alt);
            annonStr{end+1} = ffStr;
        end

        if(aShowMass==1)
            totMass = sum(stateLogRow(9:12));
            dryMass = sum(stateLogRow(9));
            fuelOxMass = sum(stateLogRow(10));
            mpMass = sum(stateLogRow(11));
            xeMass = sum(stateLogRow(12));

            propNames = maData.spacecraft.propellant.names;
            if(length(propNames{1}) > 11)
                propNames{1} = propNames{1}(1:11);
            end
            if(length(propNames{2}) > 11)
                propNames{2} = propNames{2}(1:11);
            end
            if(length(propNames{3}) > 11)
                propNames{3} = propNames{3}(1:11);
            end

            massStr = sprintf('Spacecraft Masses\n   Total Mass=%g ton\n   Dry Mass=%g ton\n   %s=%g ton\n   %s=%g ton\n   %s=%g ton', totMass, dryMass, propNames{1}, fuelOxMass, ...
                              propNames{2}, mpMass, propNames{3}, xeMass);
            annonStr{end+1} = massStr;
        end

        if(aShowFps==1)
            fps = 1/etime(curFrameTime,prevFrameTime);
            if(~isfinite(fps))
                fps = 0;
            end
            annonStr{end+1} = sprintf('Frame Rate\n   FPS=%g',fps);
        end

%         annonStr =  strtrim(char(annonStr));
        annonStrN = {};
        for(i=1:length(annonStr))
            ss = strsplit(annonStr{i},char(10))';
            for(j=1:length(ss))
                ss{j} = horzcat(ss{j},'          ');
            end
            annonStrN = vertcat(annonStrN,ss); %#ok<AGROW>
        end

        if(isempty(hAnnote) || ~ishghandle(hAnnote))
            hAnnote = annotation('textbox',[0.77 0.84 0.22 0.1], 'String',annonStrN,'Tag','dataText','VerticalAlignment','top','BackgroundColor','w','FontName','FixedWidth','FontSize',8,'FitBoxToText','on','Margin',5);
        else
            set(hAnnote,'String',annonStrN);
        end
    else
        delete(hAnnote(ishghandle(hAnnote)));
        delete(findall(mAxes,'Tag','dataText'));
    end

    prevFrameTime = curFrameTime;

    if(GLOBAL_isRecording)
        GLOBAL_VideoWriter.writeVideo(getframe(hFig,get(mAxes,'Position')));
    end

%     shading(mAxes,'interp');
    drawnow;
end

function rPosCam = getPositOfObject(maData, time, stateLogRow, typeStr, idStr, bodyInfo, celBodyData)
    switch(typeStr)
        case 'sc'
            rPosCam = [stateLogRow(2),stateLogRow(3),stateLogRow(4)];
        case 'cb'
            rPosCam = [0,0,0];
        case 'osc'
            rPosCam = getRelPositionOfOtherSpacecraft(maData, idStr, stateLogRow(1), bodyInfo, celBodyData);
        case 'grd'
            rPosCam = getRelPositionOfGroundStation(maData, idStr, stateLogRow(1), bodyInfo, celBodyData);
        case 'chBd'
            bodyOther = getBodyInfoByNumber(str2double(idStr), celBodyData);
            [rVect, ~] = getStateAtTime(bodyOther, time, bodyInfo.gm);

            rPosCam = rVect;
        otherwise
            error('Unknown type string when computing position of object for mission animator camera.');
    end

    rPosCam = reshape(rPosCam,1,3);
end

function rPosCam = getRelPositionOfGroundStation(maData, grnID, time, mainScOrbitingBodyInfo, celBodyData)
    lookingFor = str2double(grnID);

    station = [];
    stations = maData.spacecraft.stations;
    for(i=1:length(stations))
        stationTemp = stations{i};
        if(abs(stationTemp.id == lookingFor) < 1E-12)
            station = stationTemp;
            break;
        end
    end

    if(isempty(station))
        error('Did not find ground station with ID %s for use in setting camera src/tgt.', oscID);
    end

    bodyInfoStn = getBodyInfoByNumber(station.parentID, celBodyData);

    lat = station.lat;
    long = station.long;
    alt = station.alt;

    rVectStn = getInertialVectFromLatLongAlt(time, lat, long, alt, bodyInfoStn, [NaN;NaN;NaN]);
    rPosCam = -getAbsPositBetweenSpacecraftAndBody(time, rVectStn, bodyInfoStn, mainScOrbitingBodyInfo, celBodyData); %might need a negative sign here
end

function rPosCam = getRelPositionOfOtherSpacecraft(maData, oscID, time, mainScOrbitingBodyInfo, celBodyData)
    lookingFor = str2double(oscID);

    otherSC = [];
    otherSCs = maData.spacecraft.otherSC;
    for(i=1:length(otherSCs))
        otherSCTemp = otherSCs{i};
        if(abs(otherSCTemp.id - lookingFor) < 1E-12)
            otherSC = otherSCTemp;
            break;
        end
    end

    if(isempty(otherSC))
        error('Did not find other spacecraft with ID %s for use in setting camera src/tgt.', oscID);
    end

    bodyInfoOSC = getBodyInfoByNumber(otherSC.parentID, celBodyData);

    epoch = otherSC.epoch;
    sma = otherSC.sma;
    ecc = otherSC.ecc;
    inc = deg2rad(otherSC.inc);
    raan =  deg2rad(otherSC.raan);
    arg = deg2rad(otherSC.arg);
    mean = deg2rad(otherSC.mean);

    inputOrbit = [sma, ecc, inc, raan, arg, mean, epoch];
    otherScInfo = getBodyInfoStructFromOrbit(inputOrbit);

    [rVectSC, ~] = getStateAtTime(otherScInfo, time, bodyInfoOSC.gm);
    rPosCam = -getAbsPositBetweenSpacecraftAndBody(time, rVectSC, bodyInfoOSC, mainScOrbitingBodyInfo, celBodyData); %might need a negative sign here
end

function [typeStr, nameStr, idStr] = parseSrcTgtStrings(str)
    pattern1 = '([\w]*): ([\w ]*)';
    pattern2 = '(id: ([0-9\.]+)\)*?';
    tokens1 = regexp(str,pattern1,'tokens');
    tokens2 = regexp(str,pattern2,'tokens');

    typeStr = strtrim(tokens1{1}{1});
    nameStr = strtrim(tokens1{1}{2});
    if(~isempty(tokens2))
        idStr = tokens2{1}{1};
    else
        idStr = '';
    end
end

function [hCbGrdStn] = plotLocOfGrndStation(hAxis, stnStruct, validBodyID, ut, posOffset, bodyInfo)
    hCbGrdStn = [];
    if(stnStruct.parentID == validBodyID)
        lat = stnStruct.lat;
        long = stnStruct.long;
        alt = stnStruct.alt;
        color = stnStruct.color;
        markerSym = stnStruct.markerSymbol;
        rVectECI = getInertialVectFromLatLongAlt(ut, lat, long, alt, bodyInfo, [NaN;NaN;NaN]);

        hold on;
        pos = rVectECI+posOffset;
        hCbGrdStn = plot3(hAxis, pos(1), pos(2), pos(3), markerSym, 'Color',color, 'MarkerFaceColor',color, 'MarkerSize', 6);
%         hCbGrdStnText = text(pos(1), pos(2), pos(3),stnStruct.name,'BackgroundColor','w');
        hold off;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
formulae = {'C7H13A2N3O2';'C4H5A2N2O3';'CHA2NO';'C9H17A2N3O2';'C47H68O5';'C41H83NO8P';'C48H93NO8';'C13H24NO10P';...
          'C31H54NO4';'C28H44N2NaO23';'C5H10O3';'C6H11NO4';'C4H9N3O5P';'C55H95AN3O30';'C9H20ANO7P';'C39H44N4O12';'H';'Na';'C19H28O2'};

isotopeAbundance = 0; %use polyisotopic inexact mass i.e. uses all isotopes of each element weighted by natural abundance
generalFormula = 1; %NaN for unknown elements
[getMolecularMassMasses, knownMasses, unknownElements, Ematrix, elements] = getMolecularMass(formulae,isotopeAbundance,generalFormula);

%Jeff Kantor's code:
molweightMasses = molweight(formulae);

%here is the smallest difference obtained by tweaking the options
differenceMW = [0.00221999999999412;0.00127999999997996;0.000340000000001339;0.00281999999995719;0.0140999999999849;0.0121990000000096;0.0144400000000360;0.00379899999990130;0.00933999999989510;0.00847799999996823;0.00149999999999295;0.00183999999998719;0.00117900000000759;0.0166200000001027;0.00259899999997515;0.0118599999999560;0;-1.99999999850320e-06;0.00569999999993343];

%compare = [molweightMasses, getMolecularMassMasses]

%TODO not sure what is the correct code to use
res = molweightMasses - getMolecularMassMasses - differenceMW;
assert(norm(res(isfinite(res)),inf)<1e-11)

author: opencobra

repo: cobratoolbox


text:
function b = cgbsl ( abd, lda, n, ml, mu, ipvt, b, job )

%*****************************************************************************80
%
%% CGBSL solves a complex band system factored by CGBCO or CGBFA.
%
%  Discussion:
%
%    CGBSL can solve A * X = B or hermitan ( A ) * X = B.
%
%    A division by zero will occur if the input factor contains a
%    zero on the diagonal.  Technically this indicates singularity
%    but it is often caused by improper arguments or improper
%    setting of LDA.  It will not occur if the subroutines are
%    called correctly and if CGBCO has set 0.0 < RCOND
%    or CGBFA has set INFO = 0.
%
%    To compute inverse ( A ) * C where C is a matrix with P columns:
%
%      call cgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
%
%      if ( rcond is not too small ) then
%        do j = 1, p
%          call cgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
%        end do
%      end if
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%
%  Parameters:
%
%    Input, complex ABD(LDA,N), the output from CGBCO or CGBFA.
%
%    Input, integer LDA, the leading dimension of ABD.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the number of diagonals below the main diagonal.
%
%    Input, integer MU, the number of diagonals above the main diagonal.
%
%    Input, integer IPVT(N), the pivot vector from CGBCO or CGBFA.
%
%    Input, complex B(N), the right hand side.
%
%    Input, integer JOB.
%    0, to solve A*x = b,
%    nonzero, to solve hermitian(A)*x = b, where hermitian(A) is the
%    conjugate transpose.
%
%    Output, complex B(N), the solution.
%
  m = mu + ml + 1;

  if ( job == 0 )
%
%  JOB = 0, solve A * X = B.
%
%  First solve L * Y = B.
%
    if ( ml ~= 0 )

      for k = 1 : n - 1

        lm = min ( ml, n - k );
        l = ipvt(k);
        t = b(l);

        if ( l ~= k )
          b(l) = b(k);
          b(k) = t;
        end

        b(k+1:k+lm) = b(k+1:k+lm) + t * abd(m+1:m+lm,k);

      end

    end
%
%  Now solve U * X = Y.
%
    for k = n : -1 : 1
      b(k) = b(k) / abd(m,k);
      lm = min ( k, m ) - 1;
      la = m - lm;
      lb = k - lm;
      t = -b(k);
      b(lb:lb+lm-1) = b(lb:lb+lm-1) + t * transpose ( abd(la:la+lm-1,k) );
    end

  else
%
%  JOB = nonzero, solve hermitian(A) * X = B.
%
%  First solve hermitian ( U ) * Y = B.
%
    for k = 1 : n
      lm = min ( k, m ) - 1;
      la = m - lm;
      lb = k - lm;
      t = cdotc ( lm, abd(la,k), 1, b(lb), 1 );
      b(k) = ( b(k) - t ) / conj ( abd(m,k) );
    end
%
%  Now solve hermitian ( L ) * X = Y.
%
    if ( ml ~= 0 )

      for k = n-1 : -1 : 1

        lm = min ( ml, n - k );
        b(k) = b(k) + cdotc ( lm, abd(m+1,k), 1, b(k+1), 1 );
        l = ipvt(k);

        if ( l ~= k )
          t    = b(l);
          b(l) = b(k);
          b(k) = t;
        end

      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% CROSSF - Returns estimates and plots event-related coherence (ERCOH)
%        between two input data time series (X,Y). A lower panel (optionally)
%        shows the coherence phase difference between the processes.
%        In this panel, output by   > crossf(X,Y,...);
%            90 degrees (orange) means X leads Y by a quarter cycle.
%           -90 degrees (blue)   means Y leads X by a quarter cycle.
%        Coherence phase units may be radians, degrees, or msec.
%        Click on any subplot to view separately and zoom in/out.
%
% Function description:
%        Uses EITHER fixed-window, zero-padded FFTs (fastest) OR constant-Q
%        0-padded wavelet DFTs (more even sensitivity across frequencies),
%        both Hanning-tapered.  Output frequency spacing is the lowest
%        frequency ('srate'/'winsize') divided by the 'padratio'.
%
%        If an 'alpha' value is given, then bootstrap statistics are
%        computed (from a distribution of 'naccu' {200} surrogate baseline
%        data epochs) for the baseline epoch, and non-significant features
%        of the output plots are zeroed (and shown in green). The baseline
%        epoch is all windows with center latencies < the given 'baseline'
%        value, or if 'baseboot' is 1, the whole epoch.
% Usage:
%        >> [coh,mcoh,timesout,freqsout,cohboot,cohangles] ...
%                       = crossf(X,Y,frames,tlimits,srate,cycles, ...
%                                        'key1', 'val1', 'key2', val2' ...);
% Required inputs:
%       X       = first single-channel data set (1,frames*nepochs)
%       Y       = second single-channel data set (1,frames*nepochs)
%       frames  = frames per epoch                                 {default: 750}
%       tlimits = [mintime maxtime] (ms) epoch latency limits {def: [-1000 2000]}
%       srate   = data sampling rate (Hz)                          {default: 250}
%       cycles  = 0  -> Use FFTs (with constant window length)
%               = >0 -> Number of cycles in each analysis wavelet
%               = [cycles expfactor] -> if 0 < expfactor < 1,  the number
%                 of wavelet cycles expands with frequency from cycles
%                 If expfactor = 1, no expansion; if = 0, constant
%                 window length (as in FFT)                          {default: 0}
% Optional Coherence Type:
%       'type'  = ['coher'|'phasecoher'] Compute either linear coherence
%                 ('coher') or phase coherence ('phasecoher') also known
%                 as phase coupling factor' {default: 'phasecoher'}.
%       'subitc' = ['on'|'off'] subtract stimulus locked Inter-Trial Coherence
%                 from X and Y. This computes the  'intrinsic' coherence
%                 X and Y not arising from common synchronization to
%                 experimental events. See notes. {default: 'off'}
%       'shuffle' = integer indicating the number of estimates to compute
%                 bootstrap coherence based on shuffled trials. This estimates
%                 the coherence arising only from time locking of X and Y
%                 to experimental events (opposite of 'subitc')      {default: 0}
% Optional Detrend:
%       'detret' = ['on'|'off'], Linearly detrend data within epochs {def: 'off'}
%       'detrep' = ['on'|'off'], Linearly detrend data across trials {def: 'off'}
%
% Optional FFT/DFT:
%       '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  {default: ~frames/8}
%       'timesout' = Number of output latencies (int<frames-winsize)   {def: 200}
%       'padratio' = FFTlength/winsize (2^k)                         {default: 2}
%                    Multiplies the number of output frequencies by
%                    dividing their spacing. When cycles==0, frequency
%                    spacing is (low_frequency/padratio).
%       'maxfreq'  = Maximum frequency (Hz) to plot (& output if cycles>0)
%                    If cycles==0, all FFT frequencies are output  {default: 50}
%       'baseline' = Coherence baseline end latency (ms). NaN -> No baseline
%                      {default:NaN}
%       'powbase'  = Baseline spectrum to log-subtract      {default: from data}
%
% Optional Bootstrap:
%       'alpha'    = If non-0, compute two-tailed bootstrap significance prob.
%                    level. Show non-signif output values as green. {def: 0}
%       'naccu'    = Number of bootstrap replications to compute {def: 200}
%       'boottype' = ['times'|'timestrials'] Bootstrap type: Either shuffle
%                    windows ('times') or windows and trials ('timestrials')
%                    Option 'timestrials' requires more memory  {default: 'times'}
%       'memory'   = ['low'|'high'] 'low' -> decrease memory use {default: 'high'}
%       'baseboot' = Extent of bootstrap shuffling (0=to 'baseline'; 1=whole epoch)
%                    If no baseline is given (NaN), extent of bootstrap shuffling
%                    is the whole epoch                         {default: 0}
%       'rboot'    = Input bootstrap coherence limits (e.g., from CROSSF)
%                    The bootstrap type should be identical to that used
%                    to obtain the input limits. {default: compute from data}
% Optional Scalp Map:
%       'topovec'  = (2,nchans) matrix, plot scalp maps to plot {default: []}
%                    ELSE (c1,c2), plot two cartoons showing channel locations.
%       'elocs'    = Electrode location structure or file for scalp map
%                    {default: none}
%       'chaninfo' = Electrode location additional information (nose position...)
%                    {default: none}
%
% Optional Plot and Compute Features:
%       'compute'   = ['matlab'|'c'] Use C subroutines to speed up the
%                     computation (currently unimplemented) {def: 'matlab'}
%       'savecoher' - [0|1] 1 --> Accumulate the individual trial coherence
%                     vectors; output them as cohangles {default: 0 = off}
%       'plotamp'   = ['on'|'off'], Plot coherence magnitude    {def: 'on'}
%       'maxamp'    = [real] Set the maximum for the amp. scale {def: auto}
%       'plotphase' = ['on'|'off'], Plot coherence phase angle  {def: 'on'}
%       'angleunit' = Phase units: 'ms' -> msec, 'deg' -> degrees,
%                     or 'rad' -> radians                  {default: 'deg'}
%       'title'     = Optional figure title                {default:  none}
%       'vert'      = Latencies to mark with a dotted vertical line
%                                                           {default: none}
%       'linewidth' = Line width for marktimes traces (thick=2, thin=1)
%                                                              {default: 2}
%       'cmax'      = Maximum amplitude for color scale  {def: data limits}
%       'axesfont'  = Axes font size                          {default: 10}
%       'titlefont' = Title font size                          {default: 8}
%
% Outputs:
%       coh         = Matrix (nfreqs,timesout) of coherence magnitudes
%       mcoh        = Vector of mean baseline coherence at each frequency
%       timesout    = Vector of output latencies (window centers) (ms).
%       freqsout    = Vector of frequency bin centers (Hz).
%       cohboot     = Matrix (nfreqs,2) of [lower;upper] coher signif. limits
%                     if 'boottype' is 'trials',  (nfreqs,timesout, 2)
%       cohangle    = (nfreqs,timesout) matrix of coherence angles (in radians)
%       cohangles   = (nfreqs,timesout,trials) matrix of single-trial coherence
%                      angles (in radians), saved and output only if 'savecoher',1
%
% Plot description:
%   Assuming both 'plotamp' and 'plotphase' options are 'on' (=default), the upper panel
%   presents the magnitude of either phase coherence or linear coherence, depending on
%   the 'type' parameter (above). The lower panel presents the coherence phase difference
%   (in degrees). Click on any plot to pop up a new window (using 'axcopy()').
%   -- The upper left marginal panel shows mean coherence during the baseline period
%      (blue), and when significance is set, the significance threshold (dotted black-green).
%   -- The horizontal panel under the coherence magnitude image indicates the maximum
%      (green) and minimum (blue) coherence values across all frequencies. When significance
%      is set (using option 'trials' for 'boottype'), an additional curve indicates the
%      significance threshold (dotted black-green).
%
% Notes: 1) When cycles==0, nfreqs is total number of FFT frequencies.
%        2) As noted above: 'blue' coherence angle -> X leads Y; 'red' -> Y leads X
%        3) The 'boottype' should be ideally 'timesframes', but this creates high
%           memory demands, so the 'times' method must be used in many cases.
%        4) If 'boottype' is 'trials', the average of the complex bootstrap
%           is subtracted from the coherence to compensate for phase differences
%           (the average is also subtracted from the bootstrap distribution).
%           For other bootstraps, this is not necessary since the phase is random.
%        5) If baseline is non-NaN, the baseline is subtracted from
%           the complex coherence. On the left hand side of the coherence
%           amplitude image, the baseline is displayed as a magenta line
%           (if no baseline is selected, this curve represents the average
%           coherence at every given frequency).
%        6) If a out-of-memory error occurs, set the 'memory' option to 'low'
%           (Makes computation time slower; Only the 'times' bootstrap method
%           can be used in this mode).
%
% Authors: Arnaud Delorme, Sigurd Enghoff & Scott Makeig
%          CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
%
% See also: TIMEF

% Copyright (C) 8/1/98  Arnaud Delorme, Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD
%
% 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.

% 11-20-98 defined g.linewidth constant -sm
% 04-01-99 made number of frequencies consistent -se
% 06-29-99 fixed constant-Q freq indexing -se
% 08-13-99 added cohangle plotting -sm
% 08-20-99 made bootstrap more efficient -sm
% 08-24-99 allow nan values introduced by possible EVENTLOCK preproc. -sm
% 03-05-2007 eventlock.m deprecated to eegalign.m. -tf
% 03-16-00 added lead/lag interpretation to help msg - sm & eric visser
% 03-16-00 added AXCOPY feature -sm & tpj
% 04-20-00 fixed Rangle sign for wavelets, added verts array -sm
% 01-22-01 corrected help msg when nargin<2 -sm & arno delorme
% 01-25-02 reformated help & license, added links -ad
% 03-09-02 function restructuration -ad
%  add 'key', val arguments (+ external baseboot, baseline, color axis, angleunit...)
%  add detrending (across time and trials) + 'coher' option for amplitude coherence
%  significance only if alpha is given, plotting options in 'plotamp' and 'plotphase'
% 03-16-02 timeout automatically adjusted if too high -ad
% 04-03-02 added new options for bootstrap -ad

% Note: 3 "objects" (Tf, Coher and Boot) are handled by specific functions under Matlab
%    (Tf) function Tf = tfinit(...) - create object Time Frequency (Tf) associated with some data
%    (Tf) function [Tf, itcvals] = tfitc(...) - compute itc for the selected data
%    (Tf) function [Tf, itcvals] = tfitcpost(Tf, trials) - itc normlisation
%    (Tf) function [Tf, tmpX] = tfcomp(Tf, trials, times) - compute time freq. decomposition
%    (Coher) function Coher = coherinit(...) - initialize coherence object
%    (Coher) function Coher = cohercomp(Coher, tmpX, tmpY, trial, time) - compute coherence
%    (Coher) function Coher = cohercomppost(Coher, trials) - coherence normalization
%    (Boot) function Boot = bootinit(...) - initialize bootstrap object
%    (Boot) function Boot = bootcomp(...) - compute bootstrap
%    (Boot) function [Boot, Rbootout] = bootcomppost(...) - bootstrap normalization
% and by real objects under C++ (C++ code, incomplete)

function [R,mbase,times,freqs,Rbootout,Rangle, trialcoher, Tfx, Tfy] = crossf(X, Y, frame, tlimits, Fs, varwin, varargin)

%varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% ------------------------
% Commandline arg defaults:
% ------------------------
DEFAULT_ANGLEUNIT = 'deg'; % angle plotting units - 'rad', 'ms', or 'deg'
DEFAULT_EPOCH	= 750;			% Frames per epoch
DEFAULT_TIMELIM = [-1000 2000];	% Time range of epochs (ms)
DEFAULT_FS		= 250;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or base on cycles.

% =0: fix window length to nwin
% >0: set window length equal varwin cycles
%     bounded above by winsize, also determines
%     the min. freq. to be computed.

DEFAULT_OVERSMP	= 2;			% Number of times to oversample = vertical resolution
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= 'Event-Related Coherence';			% Figure title
DEFAULT_ALPHA   = NaN;			% Default two-sided significance probability threshold

if (nargin < 2)
   help crossf
   return
end

if ~iscell(X)
	if (min(size(X))~=1 || length(X)<2)
		fprintf('crossf(): X must be a row or column vector.\n');
		return
	elseif (min(size(Y))~=1 || length(Y)<2)
		fprintf('crossf(): Y must be a row or column vector.\n');
		return
	elseif (length(X) ~= length(Y))
		fprintf('crossf(): X and Y must have same length.\n');
		return
	end
end

if (nargin < 3)
   frame = DEFAULT_EPOCH;
elseif (~isnumeric(frame) || length(frame)~=1 || frame~=round(frame))
   fprintf('crossf(): Value of frames must be an integer.\n');
   return
elseif (frame <= 0)
   fprintf('crossf(): Value of frames must be positive.\n');
   return
elseif ~iscell(X) && (rem(length(X),frame) ~= 0)
   fprintf('crossf(): Length of data vectors must be divisible by frames.\n');
   return
end

if (nargin < 4)
   tlimits = DEFAULT_TIMELIM;
elseif (~isnumeric(tlimits) || sum(size(tlimits))~=3)
   error('crossf(): Value of tlimits must be a vector containing two numbers.');
elseif (tlimits(1) >= tlimits(2))
   error('crossf(): tlimits interval must be [min,max].');
end

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

if (nargin < 6)
   varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) || length(varwin)>2)
   error('crossf(): Value of cycles must be a number or a (1,2) vector.');
elseif (varwin < 0)
   error('crossf(): Value of cycles must be either zero or positive.');
end

% consider structure for these arguments
% --------------------------------------
vararginori = varargin;
for index=1:length(varargin)
	if iscell(varargin{index}), varargin{index} = { varargin{index} }; end
end
if ~isempty(varargin)
   try, g = struct(varargin{:});
   catch, error('Argument error in the {''param'', value} sequence'); end;
else
	g = [];
end

try, g.shuffle;    catch, g.shuffle = 0; end
try, g.title;      catch, g.title = DEFAULT_TITLE; end
try, g.winsize;    catch, g.winsize = max(pow2(nextpow2(frame)-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 = ''; end
try, g.alpha;      catch, g.alpha = DEFAULT_ALPHA; end;
try, g.marktimes;  catch, g.marktimes = []; end; % default no vertical lines
try, g.marktimes = g.vert;       catch, g.vert = []; end; % default no vertical lines
try, g.powbase;    catch, g.powbase = nan; end
try, g.rboot;      catch, g.rboot = nan; end
try, g.plotamp;    catch, g.plotamp = 'on'; end
try, g.plotphase;  catch, g.plotphase  = 'on'; end
try, g.plotbootsub;  catch, g.plotbootsub  = 'on'; end
try, g.detrep;     catch, g.detrep = 'off'; end
try, g.detret;     catch, g.detret = 'off'; end
try, g.baseline;   catch, g.baseline = NaN; end
try, g.baseboot;   catch, g.baseboot = 0; end
try, g.linewidth;  catch, g.linewidth = 2; end
try, g.naccu;      catch, g.naccu = 200; end
try, g.angleunit;  catch, g.angleunit = DEFAULT_ANGLEUNIT; end
try, g.cmax;       catch, g.cmax = 0; end; % 0=use data limits
try, g.type;       catch, g.type = 'phasecoher'; end;
try, g.boottype;   catch, g.boottype = 'times'; end;
try, g.subitc;     catch, g.subitc = 'off'; end
try, g.memory;     catch, g.memory = 'high'; end
try, g.compute;    catch, g.compute = 'matlab'; end
try, g.maxamp;     catch, g.maxamp = []; end
try, g.savecoher;  catch, g.savecoher = 0; end
try, g.noinput;    catch, g.noinput = 'no'; end
try, g.chaninfo;   catch, g.chaninfo = []; end

allfields = fieldnames(g);
for index = 1:length(allfields)
	switch allfields{index}
	 case { 'shuffle' 'title' 'winsize' 'pad' 'timesout' 'padratio' 'maxfreq' 'topovec' 'elocs' 'alpha' ...
		  'marktimes' 'vert' 'powbase' 'rboot' 'plotamp' 'plotphase' 'plotbootsub' 'detrep' 'detret' ...
		  'baseline' 'baseboot' 'linewidth' 'naccu' 'angleunit' 'cmax' 'type' 'boottype' 'subitc' ...
		  'memory' 'compute' 'maxamp' 'savecoher' 'noinput' 'chaninfo' };
	  case {'plotersp' 'plotitc' }, disp(['crossf warning: timef option ''' allfields{index} ''' ignored']);
	 otherwise disp(['crossf error: unrecognized option ''' allfields{index} '''']); beep; return;
	end
end

g.tlimits = tlimits;
g.frame   = frame;
g.srate   = Fs;
g.cycles  = varwin(1);
if length(varwin)>1
	g.cyclesfact = varwin(2);
else
	g.cyclesfact = 1;
end
g.type       = lower(g.type);
g.boottype   = lower(g.boottype);
g.detrep     = lower(g.detrep);
g.detret     = lower(g.detret);
g.plotphase  = lower(g.plotphase);
g.plotbootsub = lower(g.plotbootsub);
g.subitc     = lower(g.subitc);
g.plotamp    = lower(g.plotamp);
g.shuffle    = lower(g.shuffle);
g.compute    = lower(g.compute);
g.AXES_FONT  = 10;
g.TITLE_FONT = 14;

% testing arguments consistency
% -----------------------------
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.frame)
   error('Value of winsize must be less than frame length.');
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.frame-g.winsize)
   g.timesout = g.frame-g.winsize;
   disp(['Value of timesout must be <= frame-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 g.maxfreq must be a number.');
elseif (g.maxfreq <= 0)
   error('Value of g.maxfreq must be positive.');
elseif (g.maxfreq > Fs/2)
   fprintf('Warning: input value of g.maxfreq larger that Nyquist frequency %3.4 Hz\n\n',Fs/2);
end

if isempty(g.topovec)
   g.topovec = [];
elseif min(size(g.topovec))==1
   g.topovec = g.topovec(:);
   if size(g.topovec,1)~=2
      error('topovec must be a row or column vector.');
   end
end

if isempty(g.elocs)
   g.elocs = '';
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)
   fprintf('Value of g.alpha is out of the normal range [%g,0.5]\n',2/g.naccu);
   g.naccu = round(2/g.alpha);
   fprintf('  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
      fprintf('Bootstrap analysis will use data in baseline (pre-0) subwindows only.\n')
   else
      fprintf('Bootstrap analysis will use data in all subwindows.\n')
   end
end
switch g.angleunit
   case { 'rad', 'ms', 'deg' },;
   otherwise error('Angleunit must be either ''rad'', ''deg'', or ''ms''');
end;
switch g.type
   case { 'coher', 'phasecoher' 'phasecoher2' },;
   otherwise error('Type must be either ''coher'' or ''phasecoher''');
end;
switch g.boottype
   case { 'times' 'timestrials' 'trials'},;
   otherwise error('Boot type must be either ''times'', ''trials'' or ''timestrials''');
end;
if (~isnumeric(g.shuffle))
   error('Shuffle argument type must be numeric');
end
switch g.memory
   case { 'low', 'high' },;
   otherwise error('memory must be either ''low'' or ''high''');
end
if strcmp(g.memory, 'low') && ~strcmp(g.boottype, 'times')
   error(['Bootstrap type ''' g.boottype ''' cannot be used in low memory mode']);
end

switch g.compute
   case { 'matlab', 'c' },;
   otherwise error('compute must be either ''matlab'' or ''c''');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare 2 conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iscell(X)
	if length(X) ~= 2 || length(Y) ~= 2
		error('crossf: to compare conditions, X and Y input must be 2-elements cell arrays');
	end
	if ~strcmp(g.boottype, 'times')
		disp('crossf warning: The significance bootstrap type is irrelevant when comparing conditions');
	end
	for index = 1:2:length(vararginori)
		if index<=length(vararginori) % needed: if elements are deleted
			%if strcmp(vararginori{index}, 'alpha'), vararginori(index:index+1) = [];
			if strcmp(vararginori{index}, 'title'), vararginori(index:index+1) = [];
			end
		end
	end
	if iscell(g.title)
		if length(g.title) <= 2,
			g.title{3} = 'Condition 2 - condition 1';
		end
	else
		g.title = { 'Condition 1', 'Condition 2', 'Condition 2 - condition 1' };
	end

	fprintf('Running crossf on condition 1 *********************\n');
	fprintf('Note: If an out-of-memory error occurs, try reducing the\n');
	fprintf('      number of time points or number of frequencies\n');
	if ~strcmp(g.type, 'coher')
	   fprintf('Note: Type ''coher'' takes 3 times as much memory as other options!)\n');
        end
	figure;
	subplot(1,3,1); title(g.title{1});
	if ~strcmp(g.type, 'coher')
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1, 'title', ' ',vararginori{:});
	else
		[R1,mbase,times,freqs,Rbootout1,Rangle1, savecoher1, Tfx1, Tfy1] = crossf(X{1}, Y{1}, ...
								frame, tlimits, Fs, varwin, 'savecoher', 1,'title', ' ',vararginori{:});
	end
	R1 = R1.*exp(j*Rangle1); % output Rangle is in radians

	% Asking user for memory limitations
	% if ~strcmp(g.noinput, 'yes')
	%	  tmp = whos('Tfx1');
	%	  fprintf('This function will require an additional %d bytes, do you wish\n', ...
    %     tmp.bytes*6+size(savecoher1,1)*size(savecoher1,2)*g.naccu*8);
	%	  res = input('to continue (y/n) (use the ''noinput'' option to disable this message):', 's');
	%  	if res == 'n', return; end
	% end

	fprintf('\nRunning crossf on condition 2 *********************\n');
	subplot(1,3,2); title(g.title{2});
	if ~strcmp(g.type, 'coher')
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	else
		  [R2,mbase,times,freqs,Rbootout2,Rangle2, savecoher2, Tfx2, Tfy2] = crossf(X{2}, Y{2}, ...
								frame, tlimits, Fs, varwin,'savecoher', 1, 'title', ' ',vararginori{:});
	end
	R2 = R2.*exp(j*Rangle2); % output Rangle is in radians

	subplot(1,3,3); title(g.title{3});
	if isnan(g.alpha)
		plotall(R2-R1, [], [], times, freqs, mbase,  find(freqs <= g.maxfreq), g);
	else
		% accumulate coherence images (all arrays [nb_points * timesout * trials])
		% ---------------------------
		allsavedcoher = zeros(size(savecoher1,1), ...
                          size(savecoher1,2), ...
                          size(savecoher1,3)+size(savecoher2,3));
		allsavedcoher(:,:,1:size(savecoher1,3))     = savecoher1;
		allsavedcoher(:,:,size(savecoher1,3)+1:end) = savecoher2;
		clear savecoher1 savecoher2;

		if strcmp(g.type, 'coher')
			alltfx = zeros(size(Tfx1,1), size(Tfx2,2), size(Tfx1,3)+size(Tfx2,3));
			alltfx(:,:,1:size(Tfx1,3))     = Tfx1;
			alltfx(:,:,size(Tfx1,3)+1:end) = Tfx2;
			clear Tfx1 Tfx2;

			alltfy = zeros(size(Tfy1,1), size(Tfy2,2), size(Tfy1,3)+size(Tfy2,3));
			alltfy(:,:,1:size(Tfy1,3))   = Tfy1;
			alltfy(:,:,size(Tfy1,3)+1:end) = Tfy2;
			clear Tfy1 Tfy2;
		end

		coherimages = zeros(size(allsavedcoher,1), size(allsavedcoher,2), g.naccu);
		cond1trials = length(X{1})/g.frame;
		cond2trials = length(X{2})/g.frame;
		alltrials = [1:cond1trials+cond2trials];
		fprintf('Accumulating bootstrap:');

		% preprocess data
		% ---------------
		switch g.type
		 case 'coher', % take the square of alltfx and alltfy
		  alltfx = alltfx.^2;
		  alltfy = alltfy.^2;
		 case 'phasecoher', % normalize
		  allsavedcoher = allsavedcoher ./ abs(allsavedcoher);
		 case 'phasecoher2', % don't do anything
		end

		if strcmp(g.type, 'coher')
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type, alltfx, alltfy);
		else
			[coherdiff coher1 coher2] = coher2conddiff( allsavedcoher, alltrials, ...
                                                        cond1trials, g.type);
		end
		%figure; g.alpha = NaN; & to check that the new images are the same as the original
		%subplot(1,3,1); plotall(coher1, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%subplot(1,3,2); plotall(coher2, [], [], times, freqs, mbase, find(freqs <= g.maxfreq), g);
		%return;

		for index=1:g.naccu
			if rem(index,10) == 0,  fprintf(' %d',index); end
			if rem(index,120) == 0, fprintf('\n'); end

			if strcmp(g.type, 'coher')
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type, alltfx, alltfy);
			else
				coherimages(:,:,index) = coher2conddiff( allsavedcoher, shuffle(alltrials), ...
                                                        cond1trials, g.type);
			end
		end
		fprintf('\n');

		% create articially a Bootstrap object to compute significance
		Boot = bootinit( [], size(allsavedcoher,1), g.timesout, g.naccu, 0, g.baseboot, ...
														'noboottype', g.alpha, g.rboot);
		Boot.Coherboot.R = coherimages;
		Boot = bootcomppost(Boot, [], [], []);
		g.title = '';
		plotall(coherdiff, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, ...
														find(freqs <= g.maxfreq), g);
	end
	return; % ********************************** END PROCESSING TWO CONDITIONS
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shuffle trials if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if g.shuffle ~= 0
   fprintf('x and y data trials being shuffled %d times\n',g.shuffle);
   XX = reshape(X, 1, frame, length(X)/g.frame);
   YY = Y;
   X = [];
   Y = [];
   for index = 1:g.shuffle
      XX = shuffle(XX,3);
      X = [X XX(:,:)];
      Y = [Y YY];
   end
end

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

% time limits
wintime = 500*g.winsize/g.srate;
times = [g.tlimits(1)+wintime:(g.tlimits(2)-g.tlimits(1)-2*wintime)/(g.timesout-1):g.tlimits(2)-wintime];

%%%%%%%%%%
% baseline
%%%%%%%%%%
if ~isnan(g.baseline)
   baseln = find(times < g.baseline); % subtract means of pre-0 (centered) windows
   if isempty(baseln)
      baseln = 1:length(times); % use all times as baseline
      disp('Bootstrap baseline empty, using the whole epoch.');
   end
   baselength = length(baseln);
else
   baseln = 1:length(times); % use all times as baseline
   baselength = length(times); % used for bootstrap
end

%%%%%%%%%%%%%%%%%%%%
% Initialize objects
%%%%%%%%%%%%%%%%%%%%
tmpsaveall = (~isnan(g.alpha) & isnan(g.rboot) & strcmp(g.memory, 'high')) ...
                 | (strcmp(g.subitc, 'on') & strcmp(g.memory, 'high'));
trials = length(X)/g.frame;
if ~strcmp(g.compute, 'c')
   Tfx = tfinit(X, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, ...
									g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Tfy = tfinit(Y, g.timesout, g.winsize, g.cycles, g.frame, g.padratio, g.detret, ...
									g.srate, g.maxfreq, g.subitc, g.type, g.cyclesfact, tmpsaveall);
   Coher     = coherinit(Tfx.nb_points, trials, g.timesout, g.type);
   Coherboot = coherinit(Tfx.nb_points, trials, g.naccu   , g.type);
   Boot      = bootinit( Coherboot, Tfx.nb_points, g.timesout, g.naccu, baselength, ...
									g.baseboot, g.boottype, g.alpha, g.rboot);
   freqs = Tfx.freqs;
   dispf = find(freqs <= g.maxfreq);
   freqs = freqs(dispf);
else
   freqs = g.srate*g.cycles/g.winsize*[2:2/g.padratio:g.winsize]/2;
end
dispf     = find(Tfx.freqs <= g.maxfreq);

%-------------
% Reserve space
%-------------
% R  = zeros(tfx.nb_points,g.timesout);       % mean coherence
% RR = repmat(nan,tfx.nb_points,g.timesout); % initialize with nans
% Rboot = zeros(tfx.nb_points,g.naccu);  % summed bootstrap coher
% switch g.type
% case 'coher',
%    cumulXY = zeros(tfx.nb_points,g.timesout);
%    cumulXYboot = zeros(tfx.nb_points,g.naccu);
% end;
% if g.bootsub > 0
%    Rboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub); % summed bootstrap coher
%    cumulXYboottrial = zeros(tfx.nb_points, g.timesout, g.bootsub);
% end
% if ~isnan(g.alpha) & isnan(g.rboot)
%    tf.tmpalltimes = repmat(nan,tfx.nb_points,g.timesout);
% end

% --------------------
% Display text to user
% --------------------
fprintf('\nComputing Event-Related ');
switch g.type
    case 'phasecoher',  fprintf('Phase Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'phasecoher2', fprintf('Phase Coherence 2 (ITC) images for %d trials.\n',length(X)/g.frame);
    case 'coher',       fprintf('Linear Coherence (ITC) images for %d trials.\n',length(X)/g.frame);
end
fprintf('The trial latency range is from %4.5g ms before to %4.5g ms after\n     the time-locking event.\n', g.tlimits(1),g.tlimits(2));
fprintf('The frequency range displayed will be %g-%g Hz.\n',min(freqs),g.maxfreq);
if ~isnan(g.baseline)
   if length(baseln) == length(times)
      fprintf('Using the full trial latency range as baseline.\n');
   else
      fprintf('Using trial latencies from %4.5g ms to %4.5g ms as baseline.\n', g.tlimits,g.baseline);
   end
else
   fprintf('No baseline time range was specified.\n');
end
if g.cycles==0
   fprintf('The data window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
   fprintf('The FFT length will be %d samples\n',g.winsize*g.padratio);
else
   fprintf('The window size will be %2.3g cycles.\n',g.cycles);
   fprintf('The maximum window size will be %d samples (%g ms).\n',g.winsize,2*wintime);
end
fprintf('The window will be applied %d times\n',g.timesout);
fprintf('     with an average step size of %2.2g samples (%2.4g ms).\n', Tfx.stp,1000*Tfx.stp/g.srate);
fprintf('Results will be oversampled %d times.\n',g.padratio);
if ~isnan(g.alpha)
   fprintf('Bootstrap confidence limits will be computed based on alpha = %g\n', g.alpha);
else
   fprintf('Bootstrap confidence limits will NOT be computed.\n');
end
switch g.plotphase
case 'on',
    if strcmp(g.angleunit,'deg')
       fprintf(['Coherence angles will be imaged in degrees.\n']);
    elseif strcmp(g.angleunit,'rad')
       fprintf(['Coherence angles will be imaged in radians.\n']);
    elseif strcmp(g.angleunit,'ms')
       fprintf(['Coherence angles will be imaged in ms.\n']);
    end
end
fprintf('\nProcessing trial (of %d): ',trials);

% firstboot = 1;
% Rn=zeros(trials,g.timesout);
% X = X(:)'; % make X and Y column vectors
% Y = Y(:)';
% tfy = tfx;

if strcmp(g.compute, 'c')
   % C PART
   filename = [ 'tmpcrossf' num2str(round(rand(1)*1000)) ];
   f = fopen([ filename '.in'], 'w');
   fwrite(f, tmpsaveall, 'int32');
   fwrite(f, g.detret, 'int32');
   fwrite(f, g.srate, 'int32');
   fwrite(f, g.maxfreq, 'int32');
   fwrite(f, g.padratio, 'int32');
   fwrite(f, g.cycles, 'int32');
   fwrite(f, g.winsize, 'int32');
   fwrite(f, g.timesout, 'int32');
   fwrite(f, g.subitc, 'int32');
   fwrite(f, g.type, 'int32');
   fwrite(f, trials, 'int32');
   fwrite(f, g.naccu, 'int32');
   fwrite(f, length(X), 'int32');
   fwrite(f, X, 'double');
   fwrite(f, Y, 'double');
   fclose(f);

   command = [ '!cppcrosff ' filename '.in ' filename '.out' ];
   eval(command);

   f = fopen([ filename '.out'], 'r');
   size1 = fread(f, 'int32', 1);
   size2 = fread(f, 'int32', 1);
   Rreal = fread(f, 'double', [size1 size2]);
   Rimg  = fread(f, 'double', [size1 size2]);
   Coher.R = Rreal + j*Rimg;
   Boot.Coherboot.R = [];
   Boot.Rsignif = [];
else
   % ------------------------
   % MATLAB PART
   % compute ITC if necessary
   % ------------------------
   if strcmp(g.subitc, 'on')
      for t=1:trials
         if rem(t,10) == 0,  fprintf(' %d',t); end
         if rem(t,120) == 0, fprintf('\n'); end
         Tfx = tfitc( Tfx, t, 1:g.timesout);
         Tfy = tfitc( Tfy, t, 1:g.timesout);
      end;
	  fprintf('\n');
      Tfx = tfitcpost( Tfx, trials);
      Tfy = tfitcpost( Tfy, trials);
   end

   % ---------
   % Main loop
   % ---------
   if g.savecoher,
	   trialcoher = zeros(Tfx.nb_points, g.timesout, trials);
   else
	   trialcoher = [];
   end
   for t=1:trials
      if rem(t,10) == 0,  fprintf(' %d',t); end
      if rem(t,120) == 0, fprintf('\n'); end

      Tfx = tfcomp( Tfx, t, 1:g.timesout);
      Tfy = tfcomp( Tfy, t, 1:g.timesout);
	  if g.savecoher
		  [Coher trialcoher(:,:,t)] = cohercomp( Coher, Tfx.tmpalltimes, ...
                                             Tfy.tmpalltimes, t, 1:g.timesout);
      else
		  Coher = cohercomp( Coher, Tfx.tmpalltimes, Tfy.tmpalltimes, t, 1:g.timesout);
	  end

      Boot = bootcomp( Boot, Coher.Rn(t,:), Tfx.tmpalltimes, Tfy.tmpalltimes);
   end % t = trial
   [Boot Rbootout] = bootcomppost(Boot, Coher.Rn, Tfx.tmpall, Tfy.tmpall);
      % Note that the bootstrap thresholding is actually performed
      %      in the display subfunction plotall()

   Coher  = cohercomppost(Coher, trials);
end

% ----------------------------------
% If coherence, perform the division
% ----------------------------------
% switch g.type
% case 'coher',
%   R = R ./ cumulXY;
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot ./ cumulXYboot;
%   end
%   if g.bootsub > 0
%      Rboottrial = Rboottrial ./ cumulXYboottrial;
%   end
% case 'phasecoher',
%   Rn = sum(Rn, 1);
%   R = R ./ (ones(size(R,1),1)*Rn);               % coherence magnitude
%   if ~isnan(g.alpha) & isnan(g.rboot)
%      Rboot = Rboot / trials;
%   end
%   if g.bootsub > 0
%      Rboottrial = Rboottrial / trials;
%   end
% end

% ----------------
% Compute baseline
% ----------------
mbase = mean(abs(Coher.R(:,baseln)'));     % mean baseline coherence magnitude

% ---------------
% Plot everything
% ---------------
plotall(Coher.R, Boot.Coherboot.R, Boot.Rsignif, times, freqs, mbase, dispf, g);

% --------------------------------------
% Convert output Rangle to degrees or ms - Disabled to keep original default: radians output
% --------------------------------------
% Rangle = angle(Coher.R); % returns radians
% if strcmp(g.angleunit,'ms')  % convert to ms
%    Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));
% elseif strcmp(g.angleunit,'deg')  % convert to deg
%    Rangle = Rangle*180/pi; % convert to degrees
% else % angleunit is 'rad'
%    % Rangle = Rangle;
% end
% Rangle(find(Rraw==0)) = 0; % mask for significance - set angle at non-signif coher points to 0

R = abs(Coher.R);
Rsignif = Boot.Rsignif;
Tfx = permute(Tfx.tmpall, [3 2 1]); % from [trials timesout nb_points]
%                                     to   [nb_points timesout trials]
Tfy = permute(Tfy.tmpall, [3 2 1]);

return; % end crossf() *************************************************

%
% CROSSF plotting functions
% ----------------------------------------------------------------------
function plotall(R, Rboot, Rsignif, times, freqs, mbase, dispf, g)

switch lower(g.plotphase)
case 'on',
   switch lower(g.plotamp),
   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.plotamp),
   case 'on', ordinate1 = 0.1; height = 0.9;  g.plot = 1;
   case 'off', g.plot = 0;
   end;
end;

%
% Compute cross-spectral angles
% -----------------------------
Rangle = angle(R); % returns radians

%
% Optionally convert Rangle to degrees or ms
% ------------------------------------------
if strcmp(g.angleunit,'ms')  % convert to ms
   Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));
   maxangle = max(max(abs(Rangle)));
elseif strcmp(g.angleunit,'deg')  % convert to degrees
   Rangle = Rangle*180/pi; % convert to degrees
   maxangle = 180; % use full-cycle plotting
else
   maxangle = pi;  % radians
end

R = abs(R);

% if ~isnan(g.baseline)
% 	R = R - repmat(mbase',[1 g.timesout]); % remove baseline mean
% end

Rraw = R; % raw coherence (e.g., coherency) magnitude values output

if g.plot
   fprintf('\nNow plotting...\n');
   set(gcf,'DefaultAxesFontSize',g.AXES_FONT)
   colormap(jet(256));

   pos = get(gca,'position'); % plot relative to current axes
   q = [pos(1) pos(2) 0 0];
   s = [pos(3) pos(4) pos(3) pos(4)];
   axis('off')
end

switch lower(g.plotamp)
case 'on'
   %
   % Image the coherence [% perturbations]
   %
   RR = R;
   if ~isnan(g.alpha) % zero out (and 'green out') nonsignif. R values
      RR(find(RR < repmat(Rboot(:),[1 g.timesout]))) = 0;
      Rraw(find(repmat(Rsignif(:),[1,size(Rraw,2)])>=Rraw))=0;
   end

   if g.cmax == 0
      coh_caxis = max(max(R(dispf,:)))*[-1 1];
   else
      coh_caxis = g.cmax*[-1 1];
   end

   h(6) = axes('Units','Normalized', 'Position',[.1 ordinate1 .8 height].*s+q);

   map=hsv(300); % install circular color map - green=0, yellow, orng, red, violet = max
   %                                         cyan, blue, violet = min
   map = flipud([map(251:end,:);map(1:250,:)]);
   map(151,:) = map(151,:)*0.9; % tone down the (0=) green!
   colormap(map);

   imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % plot the coherence image
   if ~isempty(g.maxamp)
	   caxis([-g.maxamp g.maxamp]);
   end
   tmpscale = caxis;

   hold on
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth)
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end
   hold off
   set(h(6),'YTickLabel',[],'YTick',[])
   set(h(6),'XTickLabel',[],'XTick',[])
   %title('Event-Related Coherence')

   h(8) = axes('Position',[.95 ordinate1 .05 height].*s+q);
   cbar(h(8),151:300, [0 tmpscale(2)]); % use only positive colors (gyorv)

   %
   % Plot delta-mean min and max coherence at each time point on bottom of image
   %
   h(10) = axes('Units','Normalized','Position',[.1 ordinate1-0.1 .8 .1].*s+q);
                                                            % plot marginal means below
   Emax = max(R(dispf,:)); % mean coherence at each time point
   Emin = min(R(dispf,:)); % mean coherence at each time point
   plot(times,Emin, times, Emax, 'LineWidth',g.linewidth); hold on;
   plot([times(1) times(length(times))],[0 0],'LineWidth',0.7);
   plot([0 0],[-500 500],'--m','LineWidth',g.linewidth);
   for i=1:length(g.marktimes)
       plot([g.marktimes(i) g.marktimes(i)],[-500 500],'--m','LineWidth',g.linewidth);
   end
   if ~isnan(g.alpha) && strcmp(g.boottype, 'trials')
       % plot bootstrap significance limits (base mean +/-)
      plot(times,mean(Rboot(dispf,:)),'g','LineWidth',g.linewidth); hold on;
      plot(times,mean(Rsignif(dispf,:)),'k:','LineWidth',g.linewidth);
      axis([min(times) max(times) 0 max([Emax(:)' Rsignif(:)'])*1.2])
   else
      axis([min(times) max(times) 0 max(Emax)*1.2])
   end
   tick = get(h(10),'YTick');
   set(h(10),'YTick',[tick(1) ; tick(length(tick))])
   set(h(10),'YAxisLocation','right')
   xlabel('Time (ms)')
   ylabel('coh.')

   %
   % Plot mean baseline coherence at each freq on left side of image
   %
   h(11) = axes('Units','Normalized','Position',[0 ordinate1 .1 height].*s+q);
                                                            % plot mean spectrum
   E = abs(mbase(dispf)); % baseline mean coherence at each frequency
   plot(freqs(dispf),E,'LineWidth',g.linewidth); % plot mbase
   if ~isnan(g.alpha) % plot bootstrap significance limits (base mean +/-)
      hold on
      % plot(freqs(dispf),Rboot(:,dispf)+[E;E],'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rboot  (dispf,:),2),'g','LineWidth',g.linewidth);
      plot(freqs(dispf),mean(Rsignif(dispf,:),2),'k:','LineWidth',g.linewidth);
      axis([freqs(1) freqs(max(dispf)) 0 max([E Rsignif(:)'])*1.2]);
   else             % plot marginal mean coherence only
      if ~isnan(max(E))
         axis([freqs(1) freqs(max(dispf)) 0 max(E)*1.2]);
      end
   end

   tick = get(h(11),'YTick');
   set(h(11),'YTick',[tick(1) ; tick(length(tick))])
   set(h(11),'View',[90 90])
   xlabel('Freq. (Hz)')
   ylabel('coh.')
end

switch lower(g.plotphase)
case 'on'
   %
   % Plot coherence phase lags in bottom panel
   %
   h(13) = axes('Units','Normalized','Position',[.1 ordinate2 .8 height].*s+q);
   Rangle(find(Rraw==0)) = 0; % when plotting, mask for significance
                              % = set angle at non-signif coher points to 0

   imagesc(times,freqs(dispf),Rangle(dispf,:),[-maxangle maxangle]); % plot the
   hold on                                             % coherence phase angles
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth); % zero-time line
   for i=1:length(g.marktimes)
      plot([g.marktimes(i) g.marktimes(i)],[0 freqs(max(dispf))],'--m','LineWidth',g.linewidth);
   end

   ylabel('Freq. (Hz)')
   xlabel('Time (ms)')

   h(14)=axes('Position',[.95 ordinate2 .05 height].*s+q);
   cbar(h(14),0,[-maxangle maxangle]); % two-sided colorbar
end

if g.plot
	try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
    if (length(g.title) > 0) % plot title
        if h(6) ~= 0, axes(h(6)); else axes(h(13)); end
        %h = subplot('Position',[0 0  1 1].*s+q, 'Visible','Off');
        %h(13) = text(-.05,1.01,g.title);
        h(13) = title(g.title);
        %set(h(13),'VerticalAlignment','bottom')
        %set(h(13),'HorizontalAlignment','left')
        set(h(13),'FontSize',g.TITLE_FONT);
    end
   %
   %%%%%%%%%%%%%%% plot topoplot() %%%%%%%%%%%%%%%%%%%%%%%
   %
   if (~isempty(g.topovec))
      h(15) = subplot('Position',[-.1 .43 .2 .14].*s+q);
      if size(g.topovec,2) <= 2
         topoplot(g.topovec(1),g.elocs,'electrodes','off', ...
            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
      else
         topoplot(g.topovec(1,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
      end
      axis('square')

      h(16) = subplot('Position',[.9 .43 .2 .14].*s+q);
      if size(g.topovec,2) <= 2
         topoplot(g.topovec(2),g.elocs,'electrodes','off', ...
            'style', 'blank', 'emarkersize1chan', 10, 'chaninfo', g.chaninfo);
      else
         topoplot(g.topovec(2,:),g.elocs,'electrodes','off', 'chaninfo', g.chaninfo);
      end
      axis('square')
   end

   axcopy(gcf);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  TIME FREQUENCY   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for time freq initialisation
% -------------------------------------
function Tf = tfinit(X, timesout, winsize, ...
   cycles, frame, padratio, detret, srate, maxfreq, subitc, type, cyclesfact, saveall);
Tf.X         = X(:)'; % make X column vectors
Tf.winsize   = winsize;
Tf.cycles    = cycles;
Tf.frame     = frame;
Tf.padratio  = padratio;
Tf.detret    = detret;
Tf.stp       = (frame-winsize)/(timesout-1);
Tf.subitc    = subitc; % for ITC
Tf.type      = type; % for ITC
Tf.saveall   = saveall;
if (Tf.cycles == 0) %%%%%%%%%%%%%% constant window-length FFTs %%%%%%%%%%%%%%%%
   % Tf.freqs = srate/winsize*[1:2/padratio:winsize]/2; % incorrect for padratio > 2
   Tf.freqs = linspace(0, srate/2, length([1:2/padratio:winsize])+1);
   Tf.freqs = Tf.freqs(2:end);
   Tf.win   = hanning(winsize);
   Tf.nb_points = padratio*winsize/2;
else % %%%%%%%%%%%%%%%%%% Constant-Q (wavelet) DFTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   Tf.freqs = srate*cycles/winsize*[2:2/padratio:winsize]/2;
   Tf.win = dftfilt(winsize,maxfreq/srate,cycles,padratio,cyclesfact);
   Tf.nb_points = size(Tf.win,2);
end
Tf.tmpalltimes = zeros(Tf.nb_points, timesout);
trials = length(X)/frame;
if saveall
   Tf.tmpall     = repmat(nan,[trials timesout Tf.nb_points]);
else
   Tf.tmpall = [];
end
Tf.tmpallbool = zeros(trials,timesout);
Tf.ITCdone = 0;
if Tf.subitc
   Tf.ITC  = zeros(Tf.nb_points, timesout);
   switch Tf.type,
	   case { 'coher' 'phasecoher2' }
		Tf.ITCcumul  = zeros(Tf.nb_points, timesout);
   end
end

% function for itc
% ----------------
function [Tf, itcvals] = tfitc(Tf, trials, times);
Tf = tfcomp(Tf, trials, times);
switch Tf.type
   case 'coher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes).^2;
   case 'phasecoher2',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes; % complex coher.
      Tf.ITCcumul(:,times) = Tf.ITCcumul(:,times)+abs(Tf.tmpalltimes);
   case 'phasecoher',
      Tf.ITC(:,times)      = Tf.ITC(:,times) + Tf.tmpalltimes ./ abs(Tf.tmpalltimes);
                                                            % complex coher.
end % ~any(isnan())
return;

function [Tf, itcvals] = tfitcpost(Tf, trials);
switch Tf.type
   case 'coher',       Tf.ITC = Tf.ITC ./ sqrt(trials * Tf.ITCcumul);
   case 'phasecoher2', Tf.ITC = Tf.ITC ./ Tf.ITCcumul;
   case 'phasecoher',  Tf.ITC = Tf.ITC / trials; % complex coher.
end % ~any(isnan())

if Tf.saveall
  Tf.ITC = transpose(Tf.ITC); % do not use ' otherwise conjugate

	%imagesc(abs(Tf.ITC)); colorbar; figure;
	%squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
	%squeeze(Tf.ITC   (1,1,1:Tf.nb_points))
	%Tf.ITC = shiftdim(Tf.ITC, -1);

	Tf.ITC = repmat(shiftdim(Tf.ITC, -1), [trials 1 1]);
	Tf.tmpall = (Tf.tmpall - abs(Tf.tmpall) .* Tf.ITC) ./ abs(Tf.tmpall);

  %	for index = 1:trials
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow; figure;
  %		Tf.tmpall(index,:,:) = (Tf.tmpall(index,:,:) - Tf.tmpall(index,:,:) .* Tf.ITC)./Tf.tmpall(index,:,:);
  %		imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); drawnow;
  %		subplot(10,10, index); imagesc(squeeze(abs(Tf.tmpall(index,:,:)))); caxis([0 1]); drawnow;
  %	end
  %	squeeze(Tf.tmpall(1,1,1:Tf.nb_points))
  %	figure; axcopy;

end
Tf.ITCdone = 1;
return;

% function for time freq decomposition
% ------------------------------------
function [Tf, tmpX] = tfcomp(Tf, trials, times);
% tf is an structure containing all the information about the decomposition
for trial = trials
   for index = times
      if ~Tf.tmpallbool(trial, index) % already computed
         tmpX = Tf.X([1:Tf.winsize]+floor((index-1)*Tf.stp)+(trial-1)*Tf.frame);

         if ~any(isnan(tmpX)) % perform the decomposition
            tmpX = tmpX - mean(tmpX);
            switch Tf.detret, case 'on',
               tmpX = detrend(tmpX);
            end

            if Tf.cycles == 0 % use FFTs
               tmpX = Tf.win .* tmpX(:);
               tmpX = fft(tmpX,Tf.padratio*Tf.winsize);
               tmpX = tmpX(2:Tf.padratio*Tf.winsize/2+1);
            else
               tmpX = transpose(Tf.win) * tmpX(:);
            end
         else
            tmpX = NaN;
         end
         if Tf.ITCdone
            tmpX = (tmpX - abs(tmpX) .* Tf.ITC(:,index)) ./ abs(tmpX);
         end
         Tf.tmpalltimes(:,index) = tmpX;
         if Tf.saveall
            Tf.tmpall(trial, index,:) = tmpX;
            Tf.tmpallbool(trial, index) = 1;
         end
	  else
		  Tf.tmpalltimes(:,index) = Tf.tmpall(trial, index,:);
     end
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    COHERENCE    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for coherence initialisation
% -------------------------------------
function Coher = coherinit(nb_points, trials, timesout, type);
Coher.R  = zeros(nb_points,timesout);       % mean coherence
% Coher.RR = repmat(nan,nb_points,timesout); % initialize with nans
Coher.type = type;
Coher.Rn=zeros(trials,timesout);
switch type
 case 'coher',
  Coher.cumulX = zeros(nb_points,timesout);
  Coher.cumulY = zeros(nb_points,timesout);
 case 'phasecoher2',
  Coher.cumul  = zeros(nb_points,timesout);
end

% function for coherence calculation
% -------------------------------------
% function Coher = cohercomparray(Coher, tmpX, tmpY, trial);
% switch Coher.type
%   case 'coher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY); % complex coher.
%      Coher.cumulXY = Coher.cumulXY + abs(tmpX).*abs(tmpY);
%   case 'phasecoher',
%      Coher.R = Coher.R + tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
%      Coher.Rn(trial,:) = 1;
% end % ~any(ISNAN)

function [Coher,tmptrialcoh] = cohercomp(Coher, tmpX, tmpY, trial, time);
tmptrialcoh = tmpX.*conj(tmpY);
switch Coher.type
   case 'coher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumulX(:,time) = Coher.cumulX(:,time) + abs(tmpX).^2;
      Coher.cumulY(:,time) = Coher.cumulY(:,time) + abs(tmpY).^2;
 case 'phasecoher2',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh; % complex coher.
      Coher.cumul(:,time) = Coher.cumul(:,time) + abs(tmptrialcoh);
   case 'phasecoher',
      Coher.R(:,time) = Coher.R(:,time) + tmptrialcoh ./ abs(tmptrialcoh); % complex coher.
	  %figure; imagesc(abs(tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY))));
      Coher.Rn(trial,time) = Coher.Rn(trial,time)+1;
end % ~any(isnan())

% function for post coherence calculation
% ---------------------------------------
function Coher = cohercomppost(Coher, trials);
switch Coher.type
 case 'coher',
   Coher.R = Coher.R ./ sqrt(Coher.cumulX) ./ sqrt(Coher.cumulY);
 case 'phasecoher2',
   Coher.R = Coher.R ./ Coher.cumul;
 case 'phasecoher',
   Coher.Rn = sum(Coher.Rn, 1);
   Coher.R  = Coher.R ./ (ones(size(Coher.R,1),1)*Coher.Rn); % coherence magnitude
end

% function for 2 conditions coherence calculation
% -----------------------------------------------
function [coherimage, coherimage1, coherimage2] = coher2conddiff( allsavedcoher, alltrials, cond1trials, type, tfx, tfy);
	t1s = alltrials(1:cond1trials);
	t2s = alltrials(cond1trials+1:end);
	switch type
	 case 'coher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sqrt(sum(tfx(:,:,t1s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sqrt(sum(tfx(:,:,t2s),3)) ./ sqrt(sum(tfy(:,:,t1s),3));
	 case 'phasecoher2',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) ./ sum(abs(allsavedcoher(:,:,t1s)),3);
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) ./ sum(abs(allsavedcoher(:,:,t2s)),3);
	 case 'phasecoher',
	  coherimage1 = sum(allsavedcoher(:,:,t1s),3) / cond1trials;
	  coherimage2 = sum(allsavedcoher(:,:,t2s),3) / (size(allsavedcoher,3)-cond1trials);
	end
	coherimage = coherimage2 - coherimage1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BOOTSTRAP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function for bootstrap initialisation
% -------------------------------------
function Boot = bootinit(Coherboot, nb_points, timesout, naccu, baselength, baseboot, boottype, alpha, rboot);
Boot.Rboot       = zeros(nb_points,naccu);  % summed bootstrap coher
Boot.boottype    = boottype;
Boot.baselength  = baselength;
Boot.baseboot    = baseboot;
Boot.Coherboot   = Coherboot;
Boot.naccu       = naccu;
Boot.alpha       = alpha;
Boot.rboot       = rboot;

% function for bootstrap computation
% ----------------------------------
function Boot = bootcomp(Boot, Rn, tmpalltimesx, tmpalltimesy);
if ~isnan(Boot.alpha) && isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'times') % get g.naccu bootstrap estimates for each trial
      goodbasewins = find(Rn==1);
      if Boot.baseboot % use baseline windows only
         goodbasewins = find(goodbasewins<=Boot.baselength);
      end
      ngdbasewins = length(goodbasewins);
      j=1;
      tmpsX = zeros(size(tmpalltimesx,1), Boot.naccu);
      tmpsY = zeros(size(tmpalltimesx,1), Boot.naccu);
      if ngdbasewins > 1
         while j<=Boot.naccu
            s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
            s = goodbasewins(s);
            if ~any(isnan(tmpalltimesx(:,s(1)))) & ~any(isnan(tmpalltimesy(:,s(2))))
               tmpsX(:,j) = tmpalltimesx(:,s(1));
               tmpsY(:,j) = tmpalltimesy(:,s(2));
               j = j+1;
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end
   end
end

% handle other trial bootstrap types
% ----------------------------------
function [Boot, Rbootout] = bootcomppost(Boot, allRn, alltmpsX, alltmpsY);
trials    = size(alltmpsX, 1);
times     = size(alltmpsX, 2);
nb_points = size(alltmpsX, 3);
if ~isnan(Boot.alpha) && isnan(Boot.rboot)
   if strcmp(Boot.boottype, 'trials') % get g.naccu bootstrap estimates for each trial
      fprintf('\nProcessing trial bootstrap (of %d):',times(end));
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      Boot.fullcoherboot = zeros(nb_points, Boot.naccu, times);

      for index=1:times
         if rem(index,10) == 0,  fprintf(' %d',index); end
         if rem(index,120) == 0, fprintf('\n'); end
         for allt=1:trials
            j=1;
            while j<=Boot.naccu
               t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
               if (allRn(t(1),index) == 1) && (allRn(t(2),index) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),index,:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),index,:));
                  j = j+1;
               end
            end
            Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
         end
         Boot.Coherboot = cohercomppost(Boot.Coherboot);  % CHECK IF NECESSARY FOR ALL BOOT TYPE
         Boot.fullcoherboot(:,:,index) = Boot.Coherboot.R;
         Boot.Coherboot = coherinit(nb_points, trials, Boot.naccu, Boot.Coherboot.type);
      end
      Boot.Coherboot.R = Boot.fullcoherboot;
      Boot = rmfield(Boot, 'fullcoherboot');
   elseif strcmp(Boot.boottype, 'timestrials') % handle timestrials bootstrap
      fprintf('\nProcessing time and trial bootstrap (of %d):',trials);
      tmpsX = zeros(size(alltmpsX,3), Boot.naccu);
      tmpsY = zeros(size(alltmpsY,3), Boot.naccu );
      for allt=1:trials
         if rem(allt,10) == 0,  fprintf(' %d',allt); end
         if rem(allt,120) == 0, fprintf('\n'); end
         j=1;
         while j<=Boot.naccu
            t = ceil(rand([1 2])*trials); % random ints [1,g.timesout]
            goodbasewins = find((allRn(t(1),:) & allRn(t(2),:)) ==1);
            if Boot.baseboot % use baseline windows only
               goodbasewins = find(goodbasewins<=baselength);
            end
            ngdbasewins = length(goodbasewins);

            if ngdbasewins>1
               s = ceil(rand([1 2])*ngdbasewins); % random ints [1,g.timesout]
               s=goodbasewins(s);

               if all(allRn(t(1),s(1)) == 1) && all(allRn(t(2),s(2)) == 1)
                  tmpsX(:,j) = squeeze(alltmpsX(t(1),s(1),:));
                  tmpsY(:,j) = squeeze(alltmpsY(t(2),s(2),:));
                  j = j+1;
               end
            end
         end
         Boot.Coherboot = cohercomp(Boot.Coherboot, tmpsX, tmpsY, 1, 1:Boot.naccu);
      end
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   elseif strcmp(Boot.boottype, 'times') % boottype is 'times'
      Boot.Coherboot = cohercomppost(Boot.Coherboot);
   end
end

% test if precomputed
if ~isnan(Boot.alpha) && isnan(Boot.rboot) % if bootstrap analysis included . . .
   % 'boottype'='times' or 'timestrials', size(R)=nb_points*naccu
   % 'boottype'='trials',                 size(R)=nb_points*naccu*times
   Boot.Coherboot.R = abs (Boot.Coherboot.R);
   Boot.Coherboot.R = sort(Boot.Coherboot.R,2);

   % compute bootstrap significance level
   i = round(Boot.naccu*Boot.alpha);
   Boot.Rsignif = mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2); % significance levels for Rraw
   Boot.Coherboot.R = squeeze(mean(Boot.Coherboot.R(:,Boot.naccu-i+1:Boot.naccu),2));
   if size(Boot.Coherboot.R, 2) == 1
	   Rbootout(:,2) = Boot.Coherboot.R;
   else
	   Rbootout(:,:,2) = Boot.Coherboot.R;
   end
   % BEFORE
   %Rboot = [mean(Rboot(1:i,:)) ; mean(Rboot(g.naccu-i+1:g.naccu,:))];
elseif ~isnan(Boot.rboot)
	Boot.Coherboot.R = Boot.rboot;
	Boot.Rsignif     = Boot.rboot;
	Rbootout         = Boot.rboot;
else
	Boot.Coherboot.R = [];
	Boot.Rsignif     = [];
	Rbootout         = [];
end % NOTE: above, mean ?????

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

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab