gravenframes

 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 $

0001 function simState=gravenframes(gpsStart,gpsDur,frConfig,varargin)
0002 % GRAVENFRAMES - generate frame files of GravEn simulations
0003 %
0004 % simState=gravenframes(gpsStart,gpsDur,frConfig,<Graven parameters>)
0005 %   gpsStart - starting GPS time of time-series (seconds)
0006 %   gpsDur - duration of time-series (seconds)
0007 %   frConfig - MATLAB file (or struct) of frame parameters (see below)
0008 %   <Graven parameters> - list of GravEn input parameters
0009 %       inputFileName,seed - Generate signals based on file of sims
0010 %       nSim,inputFileName,seed - nSim random draws from file using seed
0011 %       nSim,inputFileName,external,seed - nSim random draws from file
0012 %           of waveforms which are projected on detectors
0013 %       nSim,simId,ampl,internal,external,seed - nSim signals generated
0014 %           on the fly of type simId chosen from amplitude distribution
0015 %  simState - state of random number generator
0016 %
0017 %   frame parameters are in a 'frparms' structure
0018 %    frparms.instrument - first part of frame file name
0019 %                               (default  - first letters of detectors)
0020 %    frparms.type       - second part of frame file names (default 'graven')
0021 %    frparms.duration   - frame file duration (default 16 seconds)
0022 %    frparms.bkgd       - background data - 'wn' white-noise 'asq' AS_Q data (default 'none') for unique detector
0023 %    frparms.chanDesc - cell-array of detID:sigID pairs
0024 %                               detID 'H1','H2','L1','GEO','TAMA','VIRGO'
0025 %                               sigID ('AS_Q','DARM_ERR')
0026 %                           - if no sigID, assume raw strain
0027 %                           (default {'H1:AS_Q','H2:AS_Q','L1:AS_Q'} )
0028 %    frparms.chanName - cell-array of channel names in frame file
0029 %                           (default '<detID>:GW-<sigID>')
0030 %    frparms.frameDirName - [OPTIONAL] directory with frame files which
0031 %                               are read for background data
0032 %
0033 % Frame files will be created for frparms.duration periods which cover
0034 %   the time-series.  They will be named
0035 %   <instrument>-<type>-<GPS>-<duration>.gwf'
0036 %
0037 % $Id: gravenframes.m,v 1.10 2005/05/04 13:58:02 thorne Exp $
0038 
0039 % Objective: generate frame files of graven sim events
0040 %
0041 % IF # of inputs is wrong
0042 %   CREATE error and exit
0043 % ENDIF
0044 % IF GPS start time or duration are invalid
0045 %   CREATE error and exit
0046 % ENDIF
0047 % SET GPS start, duration to integer values
0048 % GET parameters from configuration file
0049 % GET # of GravEn inputs
0050 if(nargin ~= 5 && nargin ~=6 && nargin ~=7 && nargin ~=9)
0051     msgId = 'gravenframes:wrongInput';
0052     error(msgId,'%s: wrong number of inputs',msgId);
0053 end
0054 if(~isnumeric(gpsStart) || gpsStart < 0)
0055     msgId = 'gravenframes:badInput';
0056     error(msgId,'%s: gpsStart input is bad',msgId);
0057 end
0058 gpsStart = floor(gpsStart);
0059 if(~isnumeric(gpsDur) || gpsDur < 1)
0060     msgId = 'gravenframes:badInput';
0061     error(msgId,'%s: gpsDur input is bad',msgId);
0062 end
0063 gpsDur = floor(gpsDur);
0064 [frparms,detList,sigList,nameList] = checkframeconfig(frConfig);
0065 gravenInput = varargin;
0066 numGrPar = nargin - 3;
0067 numChan = numel(sigList);
0068 
0069 fprintf('GRAVENFRAMES: Run GravEn starting at GPS %09d for %d seconds\n',...
0070     gpsStart,gpsDur);
0071 
0072 % IF instrument for frame file is not defined
0073 %   GET instrument from channels
0074 % ENDIF
0075 % PRINT out frame configuration parameters
0076 if(isempty(frparms.instrument))
0077     newInst = getinstrument(numChan,detList);
0078     frparms.instrument = newInst;
0079 end
0080 prframeconfig(frparms);
0081 
0082 % IF using local frame file directory for background data
0083 %   GET list of frame files in directory
0084 % ENDIF
0085 if(frparms.frameDirFlag == true)
0086     if(exist(frparms.frameDir,'dir') > 0)
0087         [gpsTimes,frameFiles,frameDurs]=...
0088             dir2framelist(frparms.frameDir);
0089         nFrames = numel(frameFiles);
0090         if(nFrames > 0)
0091             fprintf(1,'gravenframes: Use %d frame files in local directory\n',nFrames);
0092         else
0093             msgId = 'gravenframes:frameDirEmpty';
0094             error(msgId,'%s - no frame files in directory %s\n',...
0095                     msgId,frparms.frameDir);
0096         end
0097     else
0098         msgId = 'gravenframes:noFrameDir';
0099         error(msgId,'%s - frame directory %s not found\n',...
0100                 msgId,frparms.frameDir);
0101     end
0102 end
0103 
0104 % build list of unique detectors
0105 [numIfo,ifoList,detMatch] = getifolist(numChan,detList);
0106 
0107 % DEFINE arrays to hold the time-series data
0108 % SET offset from segment start to begin simulation
0109 %       Typically about 1 second
0110 % SET random seed = last GravEn input
0111 % LOOP over unique detectors
0112 %   SET sample rate for unique detector
0113 %   SET random state = saved random state
0114 permSeed = gravenInput{numGrPar};
0115 SIM_OFFSET_SEC = 1;
0116 timeData = cell(1,numChan);
0117 for kIfo = 1:numIfo
0118     detId = char(ifoList{kIfo});
0119     sampRate = getdetrate(detId);
0120     gravenGps = gpsStart;
0121 
0122 %   IF using white noise background
0123 %       GET white noise data
0124 %   ELSEIF using AS_Q background
0125 %       SET raw data channel to default for unique detector
0126 %       GET unique detector data using data start, duration from LDR
0127 %       IF unique detector data is NOT ok
0128 %           CREATE error message
0129 %       ENDIF
0130 %   ELSE
0131 %       CREATE empty background
0132 %   ENDIF
0133     switch frparms.bkgd
0134         case 'wn'
0135             fprintf(1,' Get white-noise data vector\n');
0136             dataError = 0;
0137             randnState = floor(abs(permSeed(1,1)*100));
0138             [hVector,tmpState] = wnseries(gpsDur,sampRate,1,...
0139                 randnState);
0140         case 'asq'
0141             if(frparms.altChannelFlag)
0142                 sigChanStr = frparms.altChannel;
0143                 sigChan = chanstruct(sigChanStr);
0144             else
0145                 sigChanStr = sprintf('%s:LSC-AS_Q',detId);
0146                 sigChan = chanstruct(sigChanStr);
0147             end
0148             fprintf(1,' Get raw data vector\n');
0149             if(frparms.frameDirFlag)
0150                 [hVector,sampRate,dataError] =...
0151                     chanvector(sigChan,gpsStart,gpsDur,...
0152                     gpsTimes,frameFiles,frameDurs);
0153             else
0154                 [hVector,sampRate,dataError] = ...
0155                     chanvector(sigChan,gpsStart,gpsDur);
0156             end
0157             if(dataError)
0158                 msgId = 'gravenframes:badData';
0159                 error(msgId,'%s: error %d in reading background data',msgId,dataError);
0160             end
0161         case 'none'
0162             hVector = zeros(sampRate*gpsDur,1);
0163         otherwise
0164             hVector = zeros(sampRate*gpsDur,1);
0165     end
0166 
0167 %   IF length of unique detector data does not match data duration
0168 %       CREATE error and exit
0169 %   ENDIF
0170 %   CREATE GravEn log file for this unique detector
0171 %   SET starting limits of simulation inside of data vector
0172 %       to allow for spill-over during projection
0173 %   IF first unique detector
0174 %       SET GravEn parameters to generate simulated injections
0175 %       SAVE log file as input to subsequent GravEn calls
0176 %   ELSE
0177 %       SET GravEn parameters to use first unique detector's log file
0178 %           as input (to get same injections)
0179 %   ENDIF
0180 %   CALL GravEn to get time-series with injections
0181     dataDurSamples = length(hVector);
0182     if (gpsDur*sampRate ~= dataDurSamples)
0183         msgId = 'gravenframes:badData';
0184         fmt = '%s:Raw data at GPS %d Dur %d sample rate %d has wrong # of samples %d - should be %d';
0185         error(msgId,fmt,msgId,gpsStart, gpsDur, sampRate,...
0186                  dataDurSamples, gpsDur*sampRate);
0187     end
0188     logFileName=sprintf('Sims-%s-%09d-%d-%s.txt',frparms.type,...
0189             gpsStart,gpsDur,detId);
0190     if(exist(logFileName,'file') > 0)
0191         delete(logFileName);
0192     end
0193     
0194     simOffsetSamples = SIM_OFFSET_SEC * sampRate;
0195     simSampleRange = [simOffsetSamples,(dataDurSamples-2*simOffsetSamples)];
0196     if(kIfo == 1)
0197         [grVars,numVars] = mkgrvars(logFileName,gravenGps,...
0198             detId,sampRate,simSampleRange,gravenInput);
0199         rerunInput = cell(1,2);
0200         rerunInput{1,1} = logFileName;
0201         rerunInput{1,2} = 0;
0202     else
0203         [grVars,numVars] = mkgrvars(logFileName,gravenGps,...
0204             detId,sampRate,simSampleRange,rerunInput);
0205     end
0206     fprintf(1,'GRAVENFRAMES: run GravEn for detector %s\n',detId);
0207     fprintf(1,'   GravEn logfile: %s\n',logFileName);
0208     [sims, simState] = graven(grVars{:});
0209 
0210 %  GravEn returns each sim as a short-time series starting at a given
0211 %  sample number.  We now need to place these sims onto the background time
0212 %   series
0213 %   LOOP through the sims
0214 %       CALCULATE where to add sim to full time series
0215 %       ADD sim to time series
0216 %   ENDLOOP
0217 %   DEFINE frame file data for this unique detector and channel
0218 %   STORE time-series with injections
0219 % ENDLOOP over unique detectors
0220     numSims = numel(sims);
0221     asqVector = hVector;
0222     darmVector = hVector;
0223     for ksim = 1:numSims
0224         simIndexVector = (sims(ksim).startSamp - 1) + (1:length(sims(ksim).h));
0225 %
0226 %  NOTE: Use an intermediate array to hold the (original data + sim),
0227 %          then set new data = intermediate array
0228 %       The difference is that d(n) = d(n) + sims requires a copy
0229 %       and can fragment memory because the because the lvalue is
0230 %       also an rvalue. The alternative using an intermediate
0231 %       variable (should) avoid the need for a large copy and the
0232 %       concommittant memory fragmentation.
0233 %
0234         simAndH = hVector(simIndexVector) + sims(ksim).h;
0235         hVector(simIndexVector) = simAndH;
0236         
0237         simAndAsq = asqVector(simIndexVector) + sims(ksim).AS_Q;
0238         asqVector(simIndexVector) = simAndAsq;
0239 
0240         simAndDarm = darmVector(simIndexVector) + sims(ksim).DARM_ERR;
0241         darmVector(simIndexVector) = simAndDarm;
0242     end
0243     clear simAndH;
0244     clear simAndAsq;
0245     clear simAndDarm;
0246     clear sims;
0247     for iChan = 1:numChan
0248         if(detMatch(iChan) == kIfo)
0249             sigId = char(sigList{iChan});
0250             switch sigId
0251                 case 'AS_Q'
0252                     timeData{iChan} = asqVector;
0253                 case 'DARM_ERR'
0254                     timeData{iChan} = darmVector;
0255                 otherwise
0256                     timeData{iChan} = hVector;
0257             end
0258             chName = char(nameList{1,iChan});
0259             frData(iChan) = struct('channel',chName,...
0260                 'data',[],...
0261                 'type','d',...
0262                 'mode','a');
0263         end
0264     end
0265     clear hVector;
0266     clear asqVector;
0267     clear darmVector;
0268 end
0269 
0270 % SET # of frames = data duration / frame length
0271 % LOOP over frames
0272 %   CLEAR new frame flag
0273 %   SET frame file name
0274 %   SET GPS start of frame
0275 %   SET GPS end of frame
0276 % ENDLOOP
0277 numFrames = ceil(gpsDur / frparms.duration);
0278 fprintf('GRAVENFRAMES: create %d frame files\n',numFrames);
0279 newFrameFlag = zeros(1,numFrames);
0280 frameBegGps = zeros(1,numFrames);
0281 frameEndGps = zeros(1,numFrames);
0282 frameFile = cell(1,numFrames);
0283 for iFr = 1:numFrames
0284     newFrameFlag(1,iFr) = true;
0285     frameBegGps(1,iFr) = gpsStart + (iFr-1)*frparms.duration;
0286     frameEndGps(1,iFr) = min((gpsStart+gpsDur),...
0287         (frameBegGps(1,iFr) + frparms.duration));
0288     frameDur(1,iFr) = frameEndGps(1,iFr) - frameBegGps(1,iFr);
0289     frStr = sprintf('%s-%s-%09d-%d.gwf',...
0290         frparms.instrument,frparms.type,...
0291         frameBegGps(1,iFr),frameDur(1,iFr));
0292     frameFile{1,iFr} = frStr;
0293 end
0294 
0295 % LOOP over frames
0296 %   SET beginning sample of frame
0297 %   SET end sample of frame
0298 %   LOOP over unique detectors
0299 %       SET frame data = time-series from beg sample to end sample
0300 %   ENDLOOP
0301 %   CREATE frame file with frame data from each unique detector
0302 % ENDLOOP
0303 for iFr = 1:numFrames
0304     myFrame = frameFile{1,iFr};
0305     for iChan = 1:numChan
0306         detId = char(detList{iChan});
0307         sampRate = getdetrate(detId);
0308 %
0309         begSamp = (frameBegGps(1,iFr) - gpsStart)*sampRate + 1;
0310         endSamp = (begSamp - 1) + frameDur(1,iFr)*sampRate;
0311         begSamp = max(begSamp,1);
0312         endSamp = min(endSamp,numel(timeData{iChan}));
0313 %
0314         frData(iChan).data = timeData{iChan}(begSamp:endSamp);
0315         if(iFr == numFrames)
0316             timeData{iChan} = [];
0317         end
0318     end
0319     if(iFr == numFrames)
0320         clear timeData;
0321     end
0322     if(exist(myFrame,'file') > 0)
0323         delete(myFrame);
0324     end
0325     fprintf(1,'GRAVENFRAMES: Create frame file %s\n',myFrame);
0326     mkframe(myFrame,frData,'n',frameDur(1,iFr),...
0327             frameBegGps(1,iFr));
0328 end
0329 return
0330 
0331 function [grVars,numVars] = mkgrvars(logFileName,gps,detId,sampFreq,startSamp,grInput)
0332 % MKGRVARS - create input parameters for call to GravEn
0333 %
0334 % grVars = mkgrvars(logFileName,gps,detId,sampFreq,startSamp,grInput)
0335 %
0336 % logFileName       string with name of log file (must include file
0337 %                   extention)
0338 % gps               GPS start time to allow selection of calibration
0339 % detId             detector identifier string; type 'help getifo' for list
0340 %                   of currently supported detectors
0341 % sampFreq          sampling frequency of detId
0342 % startSamp         1x2 vector defining range of allowed start times in
0343 %                   samples wrt the beginning of data set
0344 % grInput          special GravEn inputs
0345 
0346 % IF # of special inputs = 2
0347 %   CREATE 6-element GravEn input
0348 % ELSEIF # of special inputs = 3
0349 %   CREATE 7-element GravEn input
0350 % ELSEIF # of special inputs = 4
0351 %   CREATE 9-element GravEn input
0352 % ELSEIF # of special inputs = 6
0353 %   CREATE 11-element GravEn input
0354 % ELSE
0355 %   CREATE error
0356 % ENDIF
0357 numSpec = numel(grInput);
0358 if(numSpec == 2)
0359     grVars = cell(1,6);
0360     grVars{1} = logFileName;
0361     grVars{2} = gps;
0362     grVars{3} = detId;
0363     grVars{4} = grInput{1};    %inputFileName
0364     grVars{5} = sampFreq;
0365     grVars{6} = grInput{2};    %seed
0366 elseif(numSpec == 3)
0367     grVars = cell(1,7);
0368     grVars{1} = grInput{1};    %nSim
0369     grVars{2} = logFileName;
0370     grVars{3} = gps;
0371     grVars{4} = detId;
0372     grVars{5} = grInput{2};    %inputFileName
0373     grVars{6} = sampFreq;
0374     grVars{7} = grInput{3};    %seed
0375 elseif(numSpec == 4)
0376     grVars = cell(1,9);
0377     grVars{1} = grInput{1};    %nSim
0378     grVars{2} = logFileName;
0379     grVars{3} = gps;
0380     grVars{4} = detId;
0381     grVars{5} = grInput{2};    %inputFileName
0382     grVars{6} = startSamp;
0383     grVars{7} = grInput{3};    %external
0384     grVars{8} = sampFreq;
0385     grVars{9} = grInput{4};    %seed
0386 elseif(numSpec == 6)
0387     grVars = cell(1,11);
0388     grVars{1} = grInput{1};    %nSim
0389     grVars{2} = logFileName;
0390     grVars{3} = gps;
0391     grVars{4} = detId;
0392     grVars{5} = grInput{2};    %simId
0393     grVars{6} = grInput{3};    %ampl
0394     grVars{7} = startSamp;
0395     grVars{8} = grInput{4};    %internal
0396     grVars{9} = grInput{5};    %external
0397     grVars{10} = sampFreq;
0398     grVars{11} = grInput{6};    %seed
0399 else
0400     error('gravenframes - bad GravEn input');
0401 end
0402 numVars = numel(grVars);
0403 return
0404 
0405 function instName = getinstrument(numChan,detList)
0406 % GETINSTRUMENT - builds frame file instrument from detector list
0407 %   instName = getinstrument(numChan,detList)
0408 %
0409 %       numChan - number of channels
0410 %       detList - list of detectors for each channel
0411 %     instName - instrument name (sorted list of first characters
0412 %                       of detectors)
0413 
0414 % set instrument = first letters of detectors
0415 % LOOP over channels
0416 %   GET first letter of detector for channel
0417 %   IF instrument name is empty
0418 %       START instrument name with letter
0419 %   ELSE
0420 %       IF letter is not in instrument name
0421 %           ADD letter to instrument name
0422 %       ENDIF
0423 %   ENDIF
0424 % ENDLOOP
0425 % SORT characters alphabetically
0426 newInst = [];
0427 for iChan = 1:numChan
0428     detId = char(detList{iChan});
0429     firstChar = upper(detId(1:1));
0430     if(isempty(newInst))
0431         newInst(1:1) = firstChar;
0432         lenInst = numel(newInst);
0433     else
0434         if( isempty(strfind(newInst,firstChar)))
0435             newInst(lenInst+1:lenInst+1) = firstChar;
0436             lenInst = numel(newInst);
0437         end
0438     end
0439 end
0440 newInst = char(newInst);
0441 instName = sort(newInst);
0442 return
0443 
0444 function [numIfo,ifoList,detMatch] = getifolist(numChan,detList)
0445 % GETIFOLIST - build list of unique IFOs from detectors
0446 %   [numIfo,ifoList,detMatch] = getifolist(numChan,detList)
0447 %
0448 %       numChan - number of channels
0449 %       detList - list of detectors for each channel
0450 %     numIfo - # of unique IFOs
0451 %     ifoList - list of unique IFOs
0452 %     detMatch - pointers to unique IFO list
0453        
0454 % LOOP over signal list
0455 %   GET unique detector for signal
0456 %   IF no unique detector defined
0457 %       DEFINE first unique detector
0458 %   ELSE
0459 %       CLEAR detector match flag
0460 %       LOOP over unique detectors
0461 %           IF signal's detector matches one in list
0462 %               SET detector match flag
0463 %           ENDIF
0464 %       ENDLOOP
0465 %       IF detector match flag not set
0466 %           ADD signal's detector to list of unique detectors
0467 %       ENDIF
0468 %   ENDIF
0469 % ENDLOOP
0470 numIfo = 0;
0471 detMatch = zeros(1,numChan);
0472 for iChan = 1:numChan
0473     detId = detList{iChan};
0474     if(numIfo < 1)
0475         ifoList = cell(1,1);
0476         numIfo = 1;
0477         ifoList{numIfo} = detId;
0478         detMatch(iChan) = numIfo;
0479     else
0480         matchIfo = false;
0481         for iChk = 1:numIfo
0482             if(strcmp(detId,ifoList{iChk}) == true)
0483                 detMatch(iChan) = iChk;
0484                 matchIfo = true;
0485                 break
0486             end
0487         end
0488         if(matchIfo == false)
0489             numIfo = numIfo + 1;
0490             ifoList{numIfo} = detId;
0491             detMatch(iChan) = numIfo;
0492         end
0493     end
0494 end
0495 return