% import .xdf file--separate this process since importing takes time and annoying
%
% 09/07/2017 Makoto. Copied to /data/mobi. Process it again without using HJR.
% 05/29/2017 Makoto. Used. 801 was renamed as 831 for make the file ID incremental.
% 04/20/2017 Makoto. Used.
% 07/19/2016 Makoto. Created for busy John.

cd /data/mobi/Audiomaze/makoto/p00_originalData
allFolders = dir;
allFolders = allFolders(3:end);
for subjIdx = 9; %1:length(allFolders)
    
    % Move to the next subject.
    cd /data/mobi/Audiomaze/makoto/p00_originalData
    currentSubjId = allFolders(subjIdx).name;
    cd(currentSubjId)
    
    allXdfs = dir('*.xdf');
    for n = 1:length(allXdfs)
        
        % Determine current .xdf file
        tmpXdf = allXdfs(n).name;
        currentTaskId = tmpXdf(end-5:end-4);
        
        % Determine save path and name
        tmpDestination = ['/data/mobi/Audiomaze/makoto/p01_imported'];
        
        % Import the datasets without 'HandleJitterRemoval'.
        dataWithoutLinearlizedTimestamps = load_xdf(tmpXdf, 'HandleJitterRemoval', false);

        % Identify the indices
        streamNames = cell(length(dataWithoutLinearlizedTimestamps),1);
        for m = 1:length(dataWithoutLinearlizedTimestamps)
            streamNames{m,1} = dataWithoutLinearlizedTimestamps{m}.info.name;
        end
        fileplayIdx        = find(strcmp(streamNames, 'fileplay'));
        MAX_markersIdx     = find(strcmp(streamNames, 'MAX_markers'));
        PhaseSpaceIdx      = find(strcmp(streamNames, 'PhaseSpace'));
        trialInfoIdx       = find(strcmp(streamNames, 'trialInfo'));
        overheadIdx        = find(strcmp(streamNames, 'overhead'));
        buoysIdx           = find(strcmp(streamNames, 'buoys'));
        flourishIdx        = find(strcmp(streamNames, 'flourish'));
        noisepitchIdx      = find(strcmp(streamNames, 'noisepitch'));
        trialResultsIdx    = find(strcmp(streamNames, 'trialResults'));
        headwallIdx        = find(strcmp(streamNames, 'headwall'));
        handproximityIdx   = find(strcmp(streamNames, 'handproximity'));
        behaviorDataIdx    = find(strcmp(streamNames, 'behaviorData'));
        HEDtagsIdx         = find(strcmp(streamNames, 'HEDtags'));
        BioSemiIdx         = find(strcmp(streamNames, 'BioSemi'));
        AudioCaptureWinIdx = find(strcmp(streamNames, 'AudioCaptureWin'));
        
        % Skip this data if PhaseSpaceIdx is empty
        if isempty(PhaseSpaceIdx)
            continue
        end
        
        % Correct timestamps to biosemi stream
        constantValue = dataWithoutLinearlizedTimestamps{BioSemiIdx}.time_stamps(1);
        dataWithoutLinearlizedTimestamps{fileplayIdx}.time_stamps        = dataWithoutLinearlizedTimestamps{fileplayIdx}.time_stamps        - constantValue;
        if ~isempty(MAX_markersIdx) % Subj755 does not have this.
            dataWithoutLinearlizedTimestamps{MAX_markersIdx}.time_stamps     = dataWithoutLinearlizedTimestamps{MAX_markersIdx}.time_stamps     - constantValue;
        end
        dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_stamps      = dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_stamps      - constantValue;
        dataWithoutLinearlizedTimestamps{trialInfoIdx}.time_stamps       = dataWithoutLinearlizedTimestamps{trialInfoIdx}.time_stamps       - constantValue;
        dataWithoutLinearlizedTimestamps{overheadIdx }.time_stamps       = dataWithoutLinearlizedTimestamps{overheadIdx }.time_stamps       - constantValue;
        dataWithoutLinearlizedTimestamps{buoysIdx}.time_stamps           = dataWithoutLinearlizedTimestamps{buoysIdx}.time_stamps           - constantValue;
        dataWithoutLinearlizedTimestamps{flourishIdx}.time_stamps        = dataWithoutLinearlizedTimestamps{flourishIdx}.time_stamps        - constantValue;
        dataWithoutLinearlizedTimestamps{noisepitchIdx}.time_stamps      = dataWithoutLinearlizedTimestamps{noisepitchIdx}.time_stamps      - constantValue;
        dataWithoutLinearlizedTimestamps{trialResultsIdx}.time_stamps    = dataWithoutLinearlizedTimestamps{trialResultsIdx}.time_stamps    - constantValue;
        dataWithoutLinearlizedTimestamps{headwallIdx}.time_stamps        = dataWithoutLinearlizedTimestamps{headwallIdx}.time_stamps        - constantValue;
        dataWithoutLinearlizedTimestamps{handproximityIdx}.time_stamps   = dataWithoutLinearlizedTimestamps{handproximityIdx}.time_stamps   - constantValue;
        dataWithoutLinearlizedTimestamps{behaviorDataIdx}.time_stamps    = dataWithoutLinearlizedTimestamps{handproximityIdx}.time_stamps   - constantValue;
        dataWithoutLinearlizedTimestamps{HEDtagsIdx}.time_stamps         = dataWithoutLinearlizedTimestamps{HEDtagsIdx}.time_stamps         - constantValue;
        dataWithoutLinearlizedTimestamps{BioSemiIdx}.time_stamps         = dataWithoutLinearlizedTimestamps{BioSemiIdx}.time_stamps         - constantValue;
        dataWithoutLinearlizedTimestamps{AudioCaptureWinIdx}.time_stamps = dataWithoutLinearlizedTimestamps{AudioCaptureWinIdx}.time_stamps - constantValue;
        
        % Export to EEG
        EEG = eeg_emptyset;
        EEG.setname  = ['Subj' currentSubjId];
        EEG.filename = ['Subj' currentSubjId];
        EEG.data     = dataWithoutLinearlizedTimestamps{BioSemiIdx}.time_series;
        EEG.data     = EEG.data(2:131,:); % Ch1 is flat. Excluding Ch1, there are 128ch + 2ch (L and R mastoid, respectively). 
        EEG.times    = dataWithoutLinearlizedTimestamps{BioSemiIdx}.time_stamps;
        EEG.srate    = str2num(dataWithoutLinearlizedTimestamps{BioSemiIdx}.info.nominal_srate);
        EEG.xmin     = 0;
        EEG.xmax     = EEG.times(end);
        EEG.pnts     = length(EEG.times);
        EEG.nbchan   = size(EEG.data,1);
        
        % Import channel locations
        channelLocationFile = dir('*.sfp');
        EEG = pop_chanedit(EEG, 'load',{[pwd filesep channelLocationFile.name] 'filetype' 'autodetect'},'delete', 135, 'delete', 134, 'eval', 'chans = pop_chancenter( chans, [],[]);');
        
        % Create empty EEG.event structure
        EEG.event = struct('type', {}, 'latency',{});
        
        % Import Hand Proximity: No touching == 999, touching == 0
        handProximity            = cellfun(@str2num, dataWithoutLinearlizedTimestamps{handproximityIdx}.time_series(1,:)');
        diffHandProximity        = diff(handProximity);
        wallTouchIdx             = find(diffHandProximity < -990);
        wallTouchLatencyInSecond = dataWithoutLinearlizedTimestamps{handproximityIdx}.time_stamps(wallTouchIdx);
        withinBiosemiIdx         = wallTouchLatencyInSecond < EEG.xmax;
        wallTouchLatencyInSecond = wallTouchLatencyInSecond(withinBiosemiIdx);
        convertedLatency = zeros(length(wallTouchLatencyInSecond),1);
        for m = 1:length(wallTouchLatencyInSecond)
            convertedIdx = find(EEG.times - wallTouchLatencyInSecond(m)>0, 1); % From LSL to EEG.times
            convertedLatency(m,1) = EEG.times(convertedIdx)*EEG.srate;
        end
        convertedLatency = num2cell(convertedLatency);
        [EEG.event(1,1:length(convertedLatency)).latency] = convertedLatency{:};
        [EEG.event(1,1:length(convertedLatency)).type]    = deal('wallTouch');
        
        % Import MaxMarkers (to EEG.etc.maxMarker)
        if ~isempty(MAX_markersIdx)
            EEG.etc.audiomaze.maxMarker.events  = dataWithoutLinearlizedTimestamps{MAX_markersIdx}.time_series;
            EEG.etc.audiomaze.maxMarker.latency = dataWithoutLinearlizedTimestamps{MAX_markersIdx}.time_stamps;
        end
        
        % Import flourish.
        flourish                 = cellfun(@str2num, dataWithoutLinearlizedTimestamps{flourishIdx}.time_series(1,:)');
        flourishLatencyInSecond  = dataWithoutLinearlizedTimestamps{flourishIdx}.time_stamps;
        withinBiosemiIdx         = flourishLatencyInSecond < EEG.xmax;
        flourish                 = flourish(withinBiosemiIdx);
        flourishLatencyInSecond  = flourishLatencyInSecond(withinBiosemiIdx);
        convertedLatency = zeros(length(flourishLatencyInSecond),1);
        for m = 1:length(flourishLatencyInSecond)
            convertedIdx = find(EEG.times - flourishLatencyInSecond(m)>0, 1); % From LSL to EEG.times
            convertedLatency(m,1) = EEG.times(convertedIdx)*EEG.srate;
        end
        convertedLatency = num2cell(convertedLatency);
        currentEventLength = length(EEG.event);
        [EEG.event(1, currentEventLength+1:currentEventLength+length(convertedLatency)).latency] = convertedLatency{:};
        for m = 1:length(flourish)
            EEG.event(1, currentEventLength+m).type = ['flourish' num2str(flourish(m))];
        end
        
        % Sort events
        EEG = eeg_checkset(EEG, 'eventconsistency');
        
        % Make urevents
        EEG = eeg_checkset(EEG, 'makeur');
        
        % Import PhaseSpace data: resample PhaseSpace data using interp1 so that it
        % shares the time stamps of EEG.times (i.e. same length and same
        % interval.)
        %         phaseSpaceData = dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_series;
        %         phaseSpaceTime = dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_stamps;
        %         samplingUnitInSec = 1/dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.info.effective_srate;
        %         phaseSpaceSameSizeTimeStamp = 0:samplingUnitInSec:phaseSpaceTime(end);
        %         onsetIdx = find(phaseSpaceSameSizeTimeStamp>=dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_stamps(1),1);
        %         phaseSpaceSameSizeTimeSeries = zeros(size(dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_series,1), length(phaseSpaceSameSizeTimeStamp));
        %         phaseSpaceSameSizeTimeSeries(:,onsetIdx:onsetIdx+size(phaseSpaceData,2)-1) = phaseSpaceData;
        %         withinBiosemiIdx = phaseSpaceSameSizeTimeStamp >= EEG.xmin & phaseSpaceSameSizeTimeStamp <= EEG.xmax;
        %         phaseSpaceSameSizeTimeStamp  = phaseSpaceSameSizeTimeStamp(withinBiosemiIdx);
        %         phaseSpaceSameSizeTimeSeries = phaseSpaceSameSizeTimeSeries(:,withinBiosemiIdx);
        %         eegNormalizedPhaseSpace = (interp1(1:length(phaseSpaceSameSizeTimeStamp), phaseSpaceSameSizeTimeSeries', 1:EEG.pnts))';
        %         EEG.etc.phaseSpace = eegNormalizedPhaseSpace;
        EEG.etc.audiomaze.phaseSpaceData = dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_series;
        EEG.etc.audiomaze.phaseSpaceTime = dataWithoutLinearlizedTimestamps{PhaseSpaceIdx}.time_stamps;
        
        
        % Keep all junks under EEG.etc.audiomaze for later analysis.
        EEG.etc.audiomaze.fileplay        = dataWithoutLinearlizedTimestamps{fileplayIdx};
        EEG.etc.audiomaze.trialInfo       = dataWithoutLinearlizedTimestamps{trialInfoIdx};
        EEG.etc.audiomaze.noisepitch      = dataWithoutLinearlizedTimestamps{noisepitchIdx};
        EEG.etc.audiomaze.trialResults    = dataWithoutLinearlizedTimestamps{trialResultsIdx};
        EEG.etc.audiomaze.headwall        = dataWithoutLinearlizedTimestamps{headwallIdx};
        EEG.etc.audiomaze.handproximity   = dataWithoutLinearlizedTimestamps{handproximityIdx};
        EEG.etc.audiomaze.behaviorData    = dataWithoutLinearlizedTimestamps{behaviorDataIdx};
        EEG.etc.audiomaze.HEDtags         = dataWithoutLinearlizedTimestamps{HEDtagsIdx};
        EEG.etc.audiomaze.AudioCaptureWin = dataWithoutLinearlizedTimestamps{AudioCaptureWinIdx};
        
        % Save the data.
        EEG.subject = currentSubjId;
        pop_saveset(EEG, 'filename', ['subj' currentSubjId '_' currentTaskId], 'filepath', '/data/mobi/Audiomaze/makoto/p01_imported');
    end
end

disp('Done.')
cd /data/mobi/Audiomaze/makoto/p01_imported