function simState=gravenframes(gpsStart,gpsDur,frConfig,varargin)
% GRAVENFRAMES - generate frame files of GravEn simulations
%
% simState=gravenframes(gpsStart,gpsDur,frConfig,<Graven parameters>)
%   gpsStart - starting GPS time of time-series (seconds)
%   gpsDur - duration of time-series (seconds)
%   frConfig - MATLAB file (or struct) of frame parameters (see below)
%   <Graven parameters> - list of GravEn input parameters
%       inputFileName,seed - Generate signals based on file of sims
%       nSim,inputFileName,seed - nSim random draws from file using seed
%       nSim,inputFileName,external,seed - nSim random draws from file
%           of waveforms which are projected on detectors
%       nSim,simId,ampl,internal,external,seed - nSim signals generated
%           on the fly of type simId chosen from amplitude distribution
%  simState - state of random number generator
%
%   frame parameters are in a 'frparms' structure
%    frparms.instrument - first part of frame file name 
%                               (default  - first letters of detectors)
%    frparms.type       - second part of frame file names (default 'graven')
%    frparms.duration   - frame file duration (default 16 seconds)
%    frparms.bkgd       - background data - 'wn' white-noise 'asq' AS_Q data (default 'none') for unique detector
%    frparms.chanDesc - cell-array of detID:sigID pairs
%                               detID 'H1','H2','L1','GEO','TAMA','VIRGO'
%                               sigID ('AS_Q','DARM_ERR')
%                           - if no sigID, assume raw strain
%                           (default {'H1:AS_Q','H2:AS_Q','L1:AS_Q'} )
%    frparms.chanName - cell-array of channel names in frame file
%                           (default '<detID>:GW-<sigID>')
%    frparms.frameDirName - [OPTIONAL] directory with frame files which
%                               are read for background data
%
% Frame files will be created for frparms.duration periods which cover
%   the time-series.  They will be named
%   <instrument>-<type>-<GPS>-<duration>.gwf'
%
% $Id: gravenframes.m,v 1.10 2005/05/04 13:58:02 thorne Exp $

% Objective: generate frame files of graven sim events
%
% IF # of inputs is wrong
%   CREATE error and exit
% ENDIF
% IF GPS start time or duration are invalid
%   CREATE error and exit
% ENDIF
% SET GPS start, duration to integer values
% GET parameters from configuration file
% GET # of GravEn inputs
if(nargin ~= 5 && nargin ~=6 && nargin ~=7 && nargin ~=9)
    msgId = 'gravenframes:wrongInput';
    error(msgId,'%s: wrong number of inputs',msgId);
end
if(~isnumeric(gpsStart) || gpsStart < 0)
    msgId = 'gravenframes:badInput';
    error(msgId,'%s: gpsStart input is bad',msgId);
end
gpsStart = floor(gpsStart);
if(~isnumeric(gpsDur) || gpsDur < 1)
    msgId = 'gravenframes:badInput';
    error(msgId,'%s: gpsDur input is bad',msgId);
end
gpsDur = floor(gpsDur);
[frparms,detList,sigList,nameList] = checkframeconfig(frConfig);
gravenInput = varargin;
numGrPar = nargin - 3;
numChan = numel(sigList);

fprintf('GRAVENFRAMES: Run GravEn starting at GPS %09d for %d seconds\n',...
    gpsStart,gpsDur);

% IF instrument for frame file is not defined
%   GET instrument from channels
% ENDIF
% PRINT out frame configuration parameters
if(isempty(frparms.instrument))
    newInst = getinstrument(numChan,detList);
    frparms.instrument = newInst;
end
prframeconfig(frparms);

% IF using local frame file directory for background data
%   GET list of frame files in directory
% ENDIF
if(frparms.frameDirFlag == true)
    if(exist(frparms.frameDir,'dir') > 0)
        [gpsTimes,frameFiles,frameDurs]=...
            dir2framelist(frparms.frameDir);
        nFrames = numel(frameFiles);
        if(nFrames > 0)
            fprintf(1,'gravenframes: Use %d frame files in local directory\n',nFrames);
        else
            msgId = 'gravenframes:frameDirEmpty';
            error(msgId,'%s - no frame files in directory %s\n',...
                    msgId,frparms.frameDir);
        end
    else
        msgId = 'gravenframes:noFrameDir';
        error(msgId,'%s - frame directory %s not found\n',...
                msgId,frparms.frameDir);
    end
end

% build list of unique detectors
[numIfo,ifoList,detMatch] = getifolist(numChan,detList);

% DEFINE arrays to hold the time-series data
% SET offset from segment start to begin simulation
%       Typically about 1 second
% SET random seed = last GravEn input
% LOOP over unique detectors
%   SET sample rate for unique detector
%   SET random state = saved random state
permSeed = gravenInput{numGrPar};
SIM_OFFSET_SEC = 1;
timeData = cell(1,numChan);
for kIfo = 1:numIfo
    detId = char(ifoList{kIfo});
    sampRate = getdetrate(detId);
    gravenGps = gpsStart;

%   IF using white noise background
%       GET white noise data
%   ELSEIF using AS_Q background
%       SET raw data channel to default for unique detector
%       GET unique detector data using data start, duration from LDR
%       IF unique detector data is NOT ok
%           CREATE error message
%       ENDIF
%   ELSE
%       CREATE empty background
%   ENDIF
    switch frparms.bkgd
        case 'wn'
            fprintf(1,' Get white-noise data vector\n');
            dataError = 0;
            randnState = floor(abs(permSeed(1,1)*100));
            [hVector,tmpState] = wnseries(gpsDur,sampRate,1,...
                randnState);
        case 'asq'
            if(frparms.altChannelFlag)
                sigChanStr = frparms.altChannel;
                sigChan = chanstruct(sigChanStr);
            else
                sigChanStr = sprintf('%s:LSC-AS_Q',detId);
                sigChan = chanstruct(sigChanStr);
            end
            fprintf(1,' Get raw data vector\n');
            if(frparms.frameDirFlag)
                [hVector,sampRate,dataError] =...
                    chanvector(sigChan,gpsStart,gpsDur,...
                    gpsTimes,frameFiles,frameDurs);
            else
                [hVector,sampRate,dataError] = ...
                    chanvector(sigChan,gpsStart,gpsDur);
            end
            if(dataError)
                msgId = 'gravenframes:badData';
                error(msgId,'%s: error %d in reading background data',msgId,dataError);
            end
        case 'none'
            hVector = zeros(sampRate*gpsDur,1);
        otherwise
            hVector = zeros(sampRate*gpsDur,1);
    end

%   IF length of unique detector data does not match data duration
%       CREATE error and exit
%   ENDIF
%   CREATE GravEn log file for this unique detector
%   SET starting limits of simulation inside of data vector
%       to allow for spill-over during projection
%   IF first unique detector
%       SET GravEn parameters to generate simulated injections
%       SAVE log file as input to subsequent GravEn calls
%   ELSE
%       SET GravEn parameters to use first unique detector's log file
%           as input (to get same injections)
%   ENDIF
%   CALL GravEn to get time-series with injections
    dataDurSamples = length(hVector);
    if (gpsDur*sampRate ~= dataDurSamples)
        msgId = 'gravenframes:badData';
        fmt = '%s:Raw data at GPS %d Dur %d sample rate %d has wrong # of samples %d - should be %d';
        error(msgId,fmt,msgId,gpsStart, gpsDur, sampRate,...
                 dataDurSamples, gpsDur*sampRate);
    end
    logFileName=sprintf('Sims-%s-%09d-%d-%s.txt',frparms.type,...
            gpsStart,gpsDur,detId);
    if(exist(logFileName,'file') > 0)
        delete(logFileName);
    end
    
    simOffsetSamples = SIM_OFFSET_SEC * sampRate;
    simSampleRange = [simOffsetSamples,(dataDurSamples-2*simOffsetSamples)];
    if(kIfo == 1)
        [grVars,numVars] = mkgrvars(logFileName,gravenGps,...
            detId,sampRate,simSampleRange,gravenInput);
        rerunInput = cell(1,2);
        rerunInput{1,1} = logFileName;
        rerunInput{1,2} = 0;
    else
        [grVars,numVars] = mkgrvars(logFileName,gravenGps,...
            detId,sampRate,simSampleRange,rerunInput);
    end
    fprintf(1,'GRAVENFRAMES: run GravEn for detector %s\n',detId);
    fprintf(1,'   GravEn logfile: %s\n',logFileName);
    [sims, simState] = graven(grVars{:});

%  GravEn returns each sim as a short-time series starting at a given
%  sample number.  We now need to place these sims onto the background time
%   series
%   LOOP through the sims
%       CALCULATE where to add sim to full time series
%       ADD sim to time series
%   ENDLOOP
%   DEFINE frame file data for this unique detector and channel
%   STORE time-series with injections
% ENDLOOP over unique detectors
    numSims = numel(sims);
    asqVector = hVector;
    darmVector = hVector;
    for ksim = 1:numSims
        simIndexVector = (sims(ksim).startSamp - 1) + (1:length(sims(ksim).h));
%
%  NOTE: Use an intermediate array to hold the (original data + sim),
%          then set new data = intermediate array
%       The difference is that d(n) = d(n) + sims requires a copy
%       and can fragment memory because the because the lvalue is
%       also an rvalue. The alternative using an intermediate
%       variable (should) avoid the need for a large copy and the
%       concommittant memory fragmentation.
%
        simAndH = hVector(simIndexVector) + sims(ksim).h;
        hVector(simIndexVector) = simAndH;
        
        simAndAsq = asqVector(simIndexVector) + sims(ksim).AS_Q;
        asqVector(simIndexVector) = simAndAsq;

        simAndDarm = darmVector(simIndexVector) + sims(ksim).DARM_ERR;
        darmVector(simIndexVector) = simAndDarm;
    end
    clear simAndH;
    clear simAndAsq;
    clear simAndDarm;
    clear sims;
    for iChan = 1:numChan
        if(detMatch(iChan) == kIfo)
            sigId = char(sigList{iChan});
            switch sigId
                case 'AS_Q'
                    timeData{iChan} = asqVector;
                case 'DARM_ERR'
                    timeData{iChan} = darmVector;
                otherwise
                    timeData{iChan} = hVector;
            end
            chName = char(nameList{1,iChan});
            frData(iChan) = struct('channel',chName,...
                'data',[],...
                'type','d',...
                'mode','a');
        end
    end
    clear hVector;
    clear asqVector;
    clear darmVector;
end

% SET # of frames = data duration / frame length
% LOOP over frames
%   CLEAR new frame flag
%   SET frame file name
%   SET GPS start of frame
%   SET GPS end of frame
% ENDLOOP
numFrames = ceil(gpsDur / frparms.duration);
fprintf('GRAVENFRAMES: create %d frame files\n',numFrames);
newFrameFlag = zeros(1,numFrames);
frameBegGps = zeros(1,numFrames);
frameEndGps = zeros(1,numFrames);
frameFile = cell(1,numFrames);
for iFr = 1:numFrames
    newFrameFlag(1,iFr) = true;
    frameBegGps(1,iFr) = gpsStart + (iFr-1)*frparms.duration;
    frameEndGps(1,iFr) = min((gpsStart+gpsDur),...
        (frameBegGps(1,iFr) + frparms.duration));
    frameDur(1,iFr) = frameEndGps(1,iFr) - frameBegGps(1,iFr);
    frStr = sprintf('%s-%s-%09d-%d.gwf',...
        frparms.instrument,frparms.type,...
        frameBegGps(1,iFr),frameDur(1,iFr));
    frameFile{1,iFr} = frStr;
end

% LOOP over frames
%   SET beginning sample of frame
%   SET end sample of frame
%   LOOP over unique detectors
%       SET frame data = time-series from beg sample to end sample
%   ENDLOOP
%   CREATE frame file with frame data from each unique detector
% ENDLOOP
for iFr = 1:numFrames
    myFrame = frameFile{1,iFr};
    for iChan = 1:numChan
        detId = char(detList{iChan});
        sampRate = getdetrate(detId);
%        
        begSamp = (frameBegGps(1,iFr) - gpsStart)*sampRate + 1;
        endSamp = (begSamp - 1) + frameDur(1,iFr)*sampRate;
        begSamp = max(begSamp,1);
        endSamp = min(endSamp,numel(timeData{iChan}));
%        
        frData(iChan).data = timeData{iChan}(begSamp:endSamp);
        if(iFr == numFrames)
            timeData{iChan} = [];
        end
    end
    if(iFr == numFrames)
        clear timeData;
    end
    if(exist(myFrame,'file') > 0)
        delete(myFrame);
    end
    fprintf(1,'GRAVENFRAMES: Create frame file %s\n',myFrame);
    mkframe(myFrame,frData,'n',frameDur(1,iFr),...
            frameBegGps(1,iFr));
end
return

function [grVars,numVars] = mkgrvars(logFileName,gps,detId,sampFreq,startSamp,grInput)
% MKGRVARS - create input parameters for call to GravEn
%
% grVars = mkgrvars(logFileName,gps,detId,sampFreq,startSamp,grInput)
%
% logFileName       string with name of log file (must include file 
%                   extention)
% gps               GPS start time to allow selection of calibration
% detId             detector identifier string; type 'help getifo' for list 
%                   of currently supported detectors 
% sampFreq          sampling frequency of detId
% startSamp         1x2 vector defining range of allowed start times in
%                   samples wrt the beginning of data set
% grInput          special GravEn inputs 

% IF # of special inputs = 2
%   CREATE 6-element GravEn input
% ELSEIF # of special inputs = 3
%   CREATE 7-element GravEn input
% ELSEIF # of special inputs = 4
%   CREATE 9-element GravEn input
% ELSEIF # of special inputs = 6
%   CREATE 11-element GravEn input
% ELSE
%   CREATE error
% ENDIF
numSpec = numel(grInput);
if(numSpec == 2)
    grVars = cell(1,6);
    grVars{1} = logFileName;
    grVars{2} = gps;
    grVars{3} = detId;
    grVars{4} = grInput{1};    %inputFileName
    grVars{5} = sampFreq;
    grVars{6} = grInput{2};    %seed
elseif(numSpec == 3)
    grVars = cell(1,7);
    grVars{1} = grInput{1};    %nSim
    grVars{2} = logFileName;
    grVars{3} = gps;
    grVars{4} = detId;
    grVars{5} = grInput{2};    %inputFileName
    grVars{6} = sampFreq;
    grVars{7} = grInput{3};    %seed
elseif(numSpec == 4)
    grVars = cell(1,9);
    grVars{1} = grInput{1};    %nSim
    grVars{2} = logFileName;
    grVars{3} = gps;
    grVars{4} = detId;
    grVars{5} = grInput{2};    %inputFileName
    grVars{6} = startSamp;
    grVars{7} = grInput{3};    %external
    grVars{8} = sampFreq;
    grVars{9} = grInput{4};    %seed
elseif(numSpec == 6)
    grVars = cell(1,11);
    grVars{1} = grInput{1};    %nSim
    grVars{2} = logFileName;
    grVars{3} = gps;
    grVars{4} = detId;
    grVars{5} = grInput{2};    %simId
    grVars{6} = grInput{3};    %ampl
    grVars{7} = startSamp;
    grVars{8} = grInput{4};    %internal
    grVars{9} = grInput{5};    %external
    grVars{10} = sampFreq;
    grVars{11} = grInput{6};    %seed
else
    error('gravenframes - bad GravEn input');
end
numVars = numel(grVars);
return

function instName = getinstrument(numChan,detList)
% GETINSTRUMENT - builds frame file instrument from detector list
%   instName = getinstrument(numChan,detList)
%
%       numChan - number of channels
%       detList - list of detectors for each channel
%     instName - instrument name (sorted list of first characters
%                       of detectors)

% set instrument = first letters of detectors
% LOOP over channels
%   GET first letter of detector for channel
%   IF instrument name is empty
%       START instrument name with letter
%   ELSE
%       IF letter is not in instrument name
%           ADD letter to instrument name
%       ENDIF
%   ENDIF    
% ENDLOOP
% SORT characters alphabetically
newInst = [];
for iChan = 1:numChan
    detId = char(detList{iChan});
    firstChar = upper(detId(1:1));
    if(isempty(newInst))
        newInst(1:1) = firstChar;
        lenInst = numel(newInst);
    else
        if( isempty(strfind(newInst,firstChar)))
            newInst(lenInst+1:lenInst+1) = firstChar;
            lenInst = numel(newInst);
        end
    end
end
newInst = char(newInst);
instName = sort(newInst);
return

function [numIfo,ifoList,detMatch] = getifolist(numChan,detList)
% GETIFOLIST - build list of unique IFOs from detectors
%   [numIfo,ifoList,detMatch] = getifolist(numChan,detList)
%
%       numChan - number of channels
%       detList - list of detectors for each channel
%     numIfo - # of unique IFOs
%     ifoList - list of unique IFOs
%     detMatch - pointers to unique IFO list
       
% LOOP over signal list
%   GET unique detector for signal
%   IF no unique detector defined
%       DEFINE first unique detector
%   ELSE
%       CLEAR detector match flag
%       LOOP over unique detectors
%           IF signal's detector matches one in list
%               SET detector match flag
%           ENDIF
%       ENDLOOP
%       IF detector match flag not set
%           ADD signal's detector to list of unique detectors
%       ENDIF
%   ENDIF
% ENDLOOP
numIfo = 0;
detMatch = zeros(1,numChan);
for iChan = 1:numChan
    detId = detList{iChan};
    if(numIfo < 1)
        ifoList = cell(1,1);
        numIfo = 1;
        ifoList{numIfo} = detId;
        detMatch(iChan) = numIfo;
    else
        matchIfo = false;
        for iChk = 1:numIfo
            if(strcmp(detId,ifoList{iChk}) == true)
                detMatch(iChan) = iChk;
                matchIfo = true;
                break
            end
        end
        if(matchIfo == false)
            numIfo = numIfo + 1;
            ifoList{numIfo} = detId;
            detMatch(iChan) = numIfo;
        end
    end
end
return
