Commit c9a2d466 authored by hannahm's avatar hannahm
Browse files

gwf version

parent fb0ca9f1
"""
#!/usr/bin/env python
a version of the ARLS filter which will take GW frame files as the
data input to be filtered
It reads in the data in batches as the frame files are opened and
read all at once (no line by line reading)
"""
import numpy as np
#import matplotlib
#matplotlib.use('Agg')
#import matplotlib.pyplot as plt
#from gwpy.timeseries import TimeSeries, TimeSeriesDict
#import gwpy.signal.fft
#from astropy import units as u
# takes in primary.value, reference.value....
def batchARLS(order,timeseries,adap,delayed,P,primary,reference,lambd):
# setup for filter
#lambd = 0.9999
n = len(primary)
#delayed = np.zeros(order)
#adap = np.zeros(order)
#cancelled = np.zeros(n)
I=np.eye(order)
#P=I
print (adap, 'in func')
# run the filter on this batch
for i in range(n):
delayed[0] = reference[i]
cancelled = primary[i] - delayed.dot(adap)
K = np.array( P.dot(delayed)/(lambd+(delayed.dot(P)).dot(delayed)) )[0]
P = (I - np.transpose(np.mat(K)) * np.mat(delayed)).dot(P) / lambd
adap = adap + cancelled*K
delayed[1:order] = delayed[0:order-1]
timeseries.value[i] = cancelled
print (adap, 'in func')
return adap, delayed, P, timeseries
"""
#!/usr/bin/env python
for running the arls filter on gw frame file data - we need to read in
the data in batches and write out as we go
THIS IS NOT WORKING
"""
import os
os.environ.setdefault('PATH', '')
import argparse
import numpy as np
from gwpy.timeseries import TimeSeries, TimeSeriesDict
import arlsGWF
def getCommandLineArguments():
parser = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
#parser.add_argument('--datafile', '-d', dest='dataFileToFilter',
# required=True, type=str,
# help='Path to file containing data to be filtered')
parser.add_argument('--order', '-o', dest='order', required=True, type=int,
help='Order for the filter.')
parser.add_argument('--lambda', '-l', dest='lambd', required=False,
type=float, default=0.9999,
help='Lambda (forgetting factor) for the filter.')
parser.add_argument('--resultsDir', '-r', dest='filtedResultsDir',
required=False,
help='Path to dir where filtered data will be saved')
args = parser.parse_args()
return args
def initialFilterSetUp(order,runDir,time):
"""
setting the initial values of the filter variables
"""
# initial filter values
delayedInit = np.zeros(order)
adapInit = np.zeros(order)
PInit = np.eye(order)
# save the initial values of the filter
np.savetxt('{0}/filter/delayed_{1}.txt'.format(runDir,time),delayedInit)
np.savetxt('{0}/filter/adap_{1}.txt'.format(runDir,time), adapInit)
np.savetxt('{0}/filter/P_{1}.txt'.format(runDir,time), PInit)
#return None ##
return adapInit, delayedInit, PInit
def main():
args = getCommandLineArguments()
order = args.order
lambd = args.lambd
#startTime = 1126621184
startTime = 1166401536
batchDuration = 64#128
# get first 5 characters of time string for locating files
timeDir = str(startTime)[:5]
refDirPath = str('/archive/frames/O1/raw/H1/H-H1_R-{0}/'
.format(timeDir,startTime))
refChannelName = str('H1:PEM-CS_MAINSMON_EBAY_1_DQ')
runDir = "."#/where/to/run/the/filter"
# this will set up the initial values of the filter and save them
#initialFilterSetUp(order,runDir,startTime) ##
adapNow, delayedNow, PNow = initialFilterSetUp(order,args.filtedResultsDir,startTime)
resampleRate = 1024
for i in range(64): # temp -> should calculate batches later
print ('batch number: ', i)
endTime = startTime + batchDuration
# change to define location above when working
# get primary
#primaryFrame = str('/home/hannah.middleton/CW/makeFakeData/scripts/powerTestData/injectionsIntoO1Data/H-H1_HOFT_C02-combined-1126621184-4096.gwf')
#primaryFrame = str('/home/hannah.middleton/CW/makeFakeData/scripts/injTest/H-H1_combinedInj-1126621184-4096.gwf')
primaryFrame = str('/hdfs/frames/O2/hoft_C02/H1/H-H1_HOFT_C02-11664/H-H1_HOFT_C02-1166401536-4096.gwf')
primaryChannelName = str('H1:DCS-CALIB_STRAIN_C02')
primBatch = TimeSeries.read(primaryFrame,primaryChannelName,start=startTime,\
end=endTime,resample=resampleRate)
# get reference
referenceFrame = str('/archive/frames/O2/raw/H1/H-H1_R-11664/H-H1_R-{0}-64.gwf'.format(startTime))
referenceChannelName = str('H1:PEM-CS_MAINSMON_EBAY_1_DQ')
refBatch = TimeSeries.read(referenceFrame,referenceChannelName,start=startTime,\
end=endTime,resample=resampleRate)
maxRef = max(np.atleast_1d(refBatch.value))
scaledRefBatch = refBatch
# temporary to test doing 2 batches at once
#referenceFrameB = str('/archive/frames/O1/raw/H1/H-H1_R-11266/H-H1_R-{0}-64.gwf'.format(startTime+64))
#refBatch128 = TimeSeries.read([referenceFrame,referenceFrameB],referenceChannelName,start=startTime,\
# end=endTime,resample=resampleRate)
#maxRef = max(np.atleast_1d(refBatch128.value))
#scaledRefBatch = refBatch128 / maxRef
# a place for the filtered results to go
filteredBatchEmpty = abs(primBatch * 0.0)
filteredBatchEmpty.name = str('{}_filtered'.format(primaryChannelName))
filteredBatchEmpty.channel = str('{}_filtered'.format(primaryChannelName))
# read in filter status - not used at the moment as we are running continuously but could be handy if we do longer data sets later
#adapNow = np.loadtxt('{0}/filter/adap_{1}.txt'.format(args.filtedResultsDir,startTime)) ##
#delayedNow = np.loadtxt('{0}/filter/delayed_{1}.txt'.format(args.filtedResultsDir,startTime)) ##
#PNow = np.loadtxt('{0}/filter/P_{1}.txt'.format(args.filtedResultsDir,startTime)) ##
#print (adapNow, 'out')
# run the filter
adapNext, delayedNext, PNext, cancelled = arlsGWF.batchARLS(order, \
filteredBatchEmpty, \
adapNow, \
delayedNow, \
PNow, \
primBatch.value, \
scaledRefBatch.value, \
lambd)
#print (adapNext, 'out')
#exit()
cancelled.write('{0}/filtered_{1}-{2}.gwf'.format(args.filtedResultsDir,startTime,endTime))
# save the filter state for the next round.
np.savetxt('{0}/filter/delayed_{1}.txt'.format(args.filtedResultsDir,endTime),delayedNext)
np.savetxt('{0}/filter/adap_{1}.txt'.format(args.filtedResultsDir,endTime), adapNext)
np.savetxt('{0}/filter/P_{1}.txt'.format(args.filtedResultsDir,endTime), PNext)
#print (delayedNow)
adapNow, delayedNow, PNow = adapNext, delayedNext, PNext ##
startTime = startTime + batchDuration
if __name__ == '__main__':
main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment