Commit fb0ca9f1 authored by hannahm's avatar hannahm
Browse files

starting to collect scripts

parent aff2f93d
"""
#!/usr/bin/env python
Sofia's ARLS code (Sat Oct 13 08:32:12 2018) with some modifications to
prevent reading in the whole file at once.
This will take a file of format XXXXXXX
and create filtered version of the data.
It will write the results as it goes.
"""
import numpy as np
import scipy.fftpack as ff
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
def runningARLS(fileName,order,filteredResultsDir='./',scaleReference=1.,lambd=0.9999):
resultsFile=open(filteredResultsDir+'/filteredData.dat','w')
resultsFile.write('time\tprim\tref\tfilt\n')
with open('{0}'.format(fileName)) as fileobject:
# take first line as these are the headers
firstLine = fileobject.readline().split()
#n = len(primary)
delayed = np.zeros(order)
adap = np.zeros(order)
#cancelled = np.zeros(n)
I=np.eye(order)
P=I
for ki,line in enumerate(fileobject):
columns = line.split()
time, primary, reference = float(columns[0]),\
float(columns[1]),\
float(columns[2])
delayed[0] = reference/scaleReference
cancelled = primary-delayed.dot(adap)
#delayed[0] = reference[ki]
#cancelled[ki] = primary[ki]-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
#print K
delayed[1:order] = delayed[0:order-1]
wLength=10240
if ki/wLength == ki//wLength and ki!=0:
np.savetxt('{0}/weights/adap_{1}.txt'.format(filteredResultsDir,ki), adap)
#plt.clf()
#plt.bar(np.arange(0,len(adap),1), adap)
#plt.savefig('{0}/weights/weights_checking-{1}.png'.format(filteredResultsDir,ki))
else: pass
"""
columns = line.split()
time, primary, reference = float(columns[0]),\
float(columns[1]),\
float(columns[2])
delayed[0] = reference/scaleReference
cancelled = primary-delayed.dot(adap)
#delayed[0] = reference[ki]
#cancelled[ki] = primary[ki]-delayed.dot(adap)
K =np.array( P.dot(delayed)/(lambd+(delayed.dot(P)).dot(delayed)))[0]
adap = adap + cancelled*K
P = (I-np.transpose(np.mat(K))*np.mat(delayed)).dot(P)/lambd
#print K
delayed[1:order] = delayed[0:order-1]
wLength=10240
if ki/wLength == ki//wLength and ki!=0:
np.savetxt('{0}/weights/adap_{1}.txt'.format(filteredResultsDir,ki), adap)
#plt.clf()
#plt.bar(np.arange(0,len(adap),1), adap)
#plt.savefig('{0}/weights/weights_checking-{1}.png'.format(filteredResultsDir,ki))
else: pass
print (ki)
columns = line.split()
time, primary, reference = float(columns[0]),\
float(columns[1]),\
float(columns[2])
delayed[0] = reference/scaleReference
pi = P.dot(delayed )
#print (pi)
K =np.array( P.dot(delayed)/(lambd+(delayed.dot(P)).dot(delayed)))[0]
#print (adap)
#print(delayed)
cancelled = primary - delayed.dot(adap)
print (cancelled)
print(K)
adap = adap + cancelled * K
P = (1./lambd) * (P - np.transpose(np.mat(K)).dot(delayed).dot(P))
#P = (1./lambd) * (I - np.transpose(np.mat(K))*np.mat(delayed)).dot(P)
delayed[1:order]=delayed[0:order-1]
"""
resultsFile.write('{0}\t{1}\t{2}\t{3}\n'.format(time,primary,
reference,cancelled))
resultsFile.close()
return [cancelled, adap]
"""
#!/usr/bin/env python
To run the filter on the data - this is specific for NOT ligo
data and only reads in ascii files
Need to adapt to read in frame files for real LIGO data
"""
import argparse
import arlsASCII
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', dest='lambd',required=False,
type=float, default=0.9999,
help='set the forgetting factor for the filter')
parser.add_argument('--resultsDir', '-r', dest='filtedResultsDir',
required=False,
help='Path to dir where filtered data will be saved')
parser.add_argument('--refScale', dest='refScalingFactor',required=False,
type=float, default=1.0,
help='add a scaling factor if reference is too large')
args = parser.parse_args()
return args
def main():
# command line arguments
# setup for the filter
args = getCommandLineArguments()
unfilteredData = args.dataFileToFilter
order = args.order
resultsLocation = args.filtedResultsDir
print("running the ARLS...")
arlsASCII.runningARLS(unfilteredData, order, filteredResultsDir=resultsLocation, scaleReference=args.refScalingFactor,lambd=args.lambd)
print("done")
if __name__ == '__main__':
main()
# using a different way to make the fake data after Sofia's idea that the noise is making a mess of things for us with the simple sinusoid.
import numpy as np
import matplotlib.pyplot as plt
#import cPickle as pickle
#import plotFFT
# the primary measured signal containing the GW, clutter line and noise
def getPrimarySignal(refSignal,gwSignal,clutterAmp=2.,noiseAmp=2.,psiRefOne=0.235,refSignalTwo=False,psiRefTwo=0.343):
## parameters ##
# psi: this is to shift the reference signal before adding to the primary
# so that is it NOT the same phase as the measured reference signal
# we give to the ANC code
# noiseAmp: the amplitude of the random noise added to the primary signal
# shifting the clutter line in phase before making the primary signal
clutterLine = clutterAmp * (refSignal * np.exp(1.j*psiRefOne)).real
# second reference channel
if len(refSignalTwo)!=1:
clutterLine2 = clutterAmp * (refSignalTwo * np.exp(1.j*psiRefTwo)).real
primarySignal = clutterLine + clutterLine2 + gwSignal + noiseAmp*np.random.randn(len(refSignal))
return primarySignal
# the primary signal = clutter + gw + noise
primarySignal = clutterLine + gwSignal + noiseAmp*np.random.randn(len(refSignal))
return primarySignal # clutter line for debugging only, not used in ANC
# this is the small gravitational wave signal at freqGW
def getGWSignal(time,freqGW=4.0, h=0.05, phi=0.34):
## parameters ##
# freqGW: the gravitatioanl wave line's frequency
# h: gw amplitude
# phi: phase of gw signal#
g = h * np.sin((2. * np.pi * freqGW * time) + phi)
return g
# the reference signal
def getReferenceSignal(time,freqC=5.2,N=5,wn=np.array([0.41,0.04,0.22,0.57,0.35]),an = np.array([-0.03,0.10,-0.16,0.08,-0.07]),phin=np.array([0.0,0.1,0.2,0.3,0.4]),refAmp=0.1,refNoiseAmp=0.1,twoChannels=False):
## parameters
# N: the number of weighted sinusoids
# wn: weight of sinusoids
# an: wandering frequency of clutter line
# refNoiseAmp: amplitude of noise added to reference signal
# need to loop over t here so that we can do the sum over wn and an
complexRef = np.zeros(len(time))
for i,t in enumerate(time):
complexRef[i] = sum( wn * np.exp(2.j * np.pi * (freqC + an)*t + phin) )
realRef = refAmp*complexRef.real + refNoiseAmp*np.random.randn(len(complexRef.real)) # added some noise
# Second Ref Channel - phase shifted and different noise
if twoChannels==True:
phaseShift = np.random.rand()
realRef2 = refAmp*(complexRef* np.exp(1.j*phaseShift)).real + refNoiseAmp*np.random.randn(len(complexRef.real))
return realRef,realRef2
#plt.plot(time,realRef)
#plt.show()
return realRef # complex ref not needed
def writeLines(times,primary,ref,inj,name='tmp.dat',ref2=False):
dataFile = open(name,'w')
if len(ref2)!=1:
dataFile.write('time\tprimary\treference1\treference2\tinjected\n')
for t,p,r1,r2,i in zip(times,primary,ref,ref2,inj):
dataFile.write('{0}\t{1}\t{2}\t{3}\t{4}\n'.format(t,p,r1,r2,i))
else:
dataFile.write('time\tprimary\treference\tinjected\n')
for t,p,r,i in zip(times,primary,ref,inj):
dataFile.write('{0}\t{1}\t{2}\t{3}\n'.format(t,p,r,i))
dataFile.close()
return
def saveParams(dataParams, gWParams, refSigParams, primarySigParams, name='tmp.dat',directory='./'):
paramFile = open('{0}{1}'.format(directory,name),'w')
for key, value in dataParams.iteritems():
paramFile.write("{0}\t{1}\n".format(key, value))
for key, value in gWParams.iteritems():
paramFile.write("{0}\t{1}\n".format(key, value))
for key, value in refSigParams.iteritems():
paramFile.write("{0}\t{1}\n".format(key, value))
for key, value in primarySigParams.iteritems():
paramFile.write("{0}\t{1}\n".format(key, value))
paramFile.close()
return
#testing here
import configparser
import argparse
parser = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('--inifile', '-i', required=True, type=str,
help='Path to ini file')
args = parser.parse_args()
iniFile = str(args.inifile)
config = configparser.ConfigParser()
config.read(iniFile)
###########################
##### data parameters #####
###########################
startTime = float(config['data']['start_time'])
endTime = float(config['data']['end_time'])
srate = float(config['data']['srate'])
timeStep = 1./srate
times = np.arange(startTime,endTime,timeStep)
seed = int(config['data']['seed'])
np.random.seed(seed)
####################################
###### ref/clutter parameters ######
####################################
fC = float(config['clutter']['centralFreq'])
weightSines = np.atleast_1d([float(str_val) for str_val in config['clutter']['weightSines'].split(',')])
freqsSines = np.atleast_1d([float(str_val) for str_val in config['clutter']['freqSines'].split(',')])
phaseSines = np.atleast_1d([float(str_val) for str_val in config['clutter']['phaseSines'].split(',')])
noSinusoids = int(len(weightSines)) # N
noiseAmpRef = float(config['clutter']['noiseAmp']) #0.01
referenceAmp = float(config['clutter']['refAmp']) #0.1 #0.1#1.E-19
#realRefSig, realRefSigTwo = getReferenceSignal(times, freqC=fC, N=noSinusoids, wn=weightSines, an=freqsSines, refAmp=referenceAmp, refNoiseAmp=noiseAmpRef, twoChannels=True)
################################
##### Injection parameters #####
################################
fGW = float(config['injection']['freq'] )#59.98
hAmp = float(config['injection']['hAmp']) #5e-19 # THIS IS THE h AMPLITUDE IN THE PRIMARY
phiGW = float(config['injection']['phi']) #0.0
#gravWaveSignal = getGWSignal(times,freqGW=fGW,h=hAmp,phi=phiGW)
###############################
##### Primary parameters ######
###############################
primaryNoiseAmp = float(config['primary']['primNoiseAmp']) #1.E-18 #0.00000001#2.0#1.E-19 #2.0#0.01# 2.0
clutterAmpInPrimary = float(config['primary']['clutterAmpInPrim']) #THIS IS THE REF AMPLITUDE IN THE PRIMARY
# save parameters used to generate fake data
dirName = '.'.format(fC,fGW)
# trying to avoid loading everything into memory at once and save each time step as we go
#dataSaveFileName='{0}/fC{1}_fGW{2}_hAmp{3}.dat'.format(dirName,fC,fGW,hAmp)
dataSaveFileName='data.dat'.format(fC,fGW,hAmp)
rollingSaveFile = open(dirName+'/'+dataSaveFileName, 'w')
rollingSaveFile.write('time\tprimary\treference1\tinjected\n')
batchSize=10240
batchNo=int((srate*endTime)/batchSize)
print ('no. batches: ', batchNo)
tStart=0.0
for b in range(batchNo):
tEnd = tStart+(batchSize/srate)
times = np.arange(tStart,tEnd,timeStep)
realRefRoll, realRefTwoRoll = getReferenceSignal(times, freqC=fC, N=noSinusoids, \
wn=weightSines, an=freqsSines, \
phin=phaseSines,\
refAmp=referenceAmp, \
refNoiseAmp=noiseAmpRef, \
twoChannels=True)
gWRoll = getGWSignal(times,freqGW=fGW,h=hAmp,phi=phiGW)
primRoll = getPrimarySignal(realRefRoll, gWRoll, clutterAmp=clutterAmpInPrimary, \
noiseAmp=primaryNoiseAmp,refSignalTwo=realRefTwoRoll, \
psiRefOne=0.111)
for t,p,r,i in zip(times,primRoll,realRefRoll,gWRoll):
rollingSaveFile.write('{0}\t{1}\t{2}\t{3}\n'.format(t,p,r,i))
tStart=tEnd
rollingSaveFile.close()
"""
# seed, dataParams, refSigParams, gWSigParams, primarySigParams
dataParams = { "randomSeed": seed,
"startTime": startTime,
"endTime": endTime,
"timeStep": timeStep,
"samplingRate": int(srate) }
gWParams = { "gwFrequency": fGW,
"gwAmplitude": hAmp,
"gwPhase": phiGW}
refSigParams = { "clutterFrequency": fC,
"numberSinusoids": noSinusoids,
"sinusoidWeights": weightSines,
"sinusoidFrequencies": freqsSines,
"referenceSigNoiseAmp": noiseAmpRef}
primarySigParams = { "primarySigNoiseAmp": primaryNoiseAmp,
"clutterAmpInPrimary": clutterAmpInPrimary}
dataDuration = int(endTime - startTime)
runDetails = 'srate{0}_dur{1}_cAmpInPrim{2}_hAmp{3}'.format(int(dataParams["samplingRate"]),dataDuration,clutterAmpInPrimary,hAmp)
paramFileName = ('params_{1}.dat'.format(dirName,runDetails))
dataFileName = ('data_{1}.dat'.format(dirName,runDetails))
#print 'writing fake data to file'
#writeLines(times,primarySignal,realRefSig,gravWaveSignal,name=dataFileName,ref2=realRefSigTwo)
#print 'done'
print ('saving data parameters')
saveParams(dataParams, gWParams, refSigParams, primarySigParams,name=paramFileName,directory=dirName)
print ('done')
"""
#plt.plot(times,primarySignal,alpha=0.7,color='r')
#plt.plot(times,clutter, alpha=0.7, color='c')
#plt.plot(times,realRefSig,alpha=0.7,color='g')
#plt.plot(times,gravWaveSignal,alpha=0.7,color='k')
#plt.show()
exit()
def makeFFT(time,value):
n = len(time)
print(n/2)
k = np.arange(n)
T = max(time)
frq = k/T
frq = frq[range(int(n/2))]
fft = np.fft.fft(value) / n
fft = fft[range(int(n/2))]
return frq, fft
#colours
orange='#E69F00'
blue='#0072B2'
green='#009E73'
ltblue='#56B4E6'
vermillion='#D55E00'
yellow='#F0E442'
pink='#CC79A7'
dataIdentifier='fC60.0_fGW59.98_hAmp1e-19'
dataToPlot = 'ARLS/data/testSet/{}.dat'.format(dataIdentifier)
dataToPlot = 'test.dat'
data = np.genfromtxt(dataToPlot,names=True)
time = data['time']
prim = data['primary']
ref = data['reference1']
inj = data['injected']
frqP,fftPrim = makeFFT(time,prim)
frqR,fftRef = makeFFT(time,ref)
frqI,fftInj = makeFFT(time,inj)
plt.clf()
plt.clf()
fig, ax = plt.subplots(3, 1)
#plt.vlines(60.01,1E-24,1E-18,color='g',alpha=0.3)
ax[0].plot(frqP,abs(fftPrim),alpha=0.6,color='k',label='primary',lw=2)
#ax[0].plot(frqF,abs(fftFilt),alpha=0.7,color=orange,label='filtered',lw=2)
#ax[0].plot(frqF,abs(fftFilt),alpha=0.7,color=orange,label='filtered',lw=2)
ax[0].legend()
#ax[0].set_xlim(59.5,60.5)
#plt.xlim(1,130)
#ax[0].set_ylim(1E-22,3E-18)
ax[0].set_xlabel('frequency / Hz')
ax[0].set_ylabel('|asd| (fft)')
ax[0].set_yscale('log')
ax[0].set_xlim(58, 62)
ax[2].plot(frqI,abs(fftInj),alpha=0.6,color='b',label='injection',lw=2)
ax[2].legend()
ax[2].set_xlim(58, 62)
print(fftRef)
print (frqR)
ax[1].plot(frqR,abs(fftRef),alpha=0.6,color=vermillion,label='reference',lw=2)
#ax[1].set_xlim(59.5,60.5)
ax[1].set_ylim(min(abs(fftRef)),max(abs(fftRef)))
ax[1].legend()
#ax[1].set_ylim(1E-5,3E-2)
ax[1].set_ylabel('|asd| (fft)')
ax[1].set_xlabel('frequency / Hz')
ax[1].set_yscale('log')
ax[1].set_xlim(58,62)
plt.tight_layout()
plt.savefig('test2.pdf')
plt.show()
[data]
start_time = 0.0
end_time = 1800.
srate = 1024
seed = 428
[clutter]
centralFreq = 60.0
freqSines = -0.495, -0.496, -0.497, -0.498, -0.499, -0.5, -0.501, -0.502, -0.503, -0.504, -0.505, -0.295, -0.296, -0.297, -0.298, -0.299, -0.3, -0.301, -0.302, -0.303, -0.304, -0.305, -0.005, -0.004, -0.003, -0.002, -0.001, 0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.295, 0.296, 0.297, 0.298, 0.299, 0.3, 0.301, 0.302, 0.303, 0.304, 0.305, 0.495, 0.496, 0.497, 0.498, 0.499, 0.5, 0.501, 0.502, 0.503, 0.504, 0.505
weightSines = 0.01, 0.01, 0.01, 0.01, 0.01, 0.1, 0.01, 0.01, 0.01, 0.01, 0.01, 0.05, 0.05, 0.05, 0.05, 0.05, 0.2, 0.05, 0.05, 0.05, 0.05, 0.05, 0.1, 0.1, 0.1, 0.1, 0.15, 2.0, 0.15, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.2, 0.05, 0.05, 0.05, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01, 0.01, 0.1, 0.01, 0.01, 0.01, 0.01, 0.01
phaseSines = 1.88464592, 1.16476356, 0.59260117, 0.4459515 , 0.98339819, 0.35539645, 1.80759105, 0.94270568, 1.35353234, 0.55839407, 0.12532757, 0.80392952, 0.46307127, 0.40150183, 1.77987356, 2.45427921, 0.77851773, 0.46634746, 1.47223271, 2.54901271, 1.3066802 , 2.79936517, 1.56310523, 0.58387607, 2.10591965, 1.0308672 , 1.44644829, 0.79843965, 0.68986151, 1.23491135, 0.12295101, 3.0966532 , 1.87212376, 0.26505531, 1.86960676, 2.2326574 , 2.55250239, 2.01754675, 1.04635494, 1.04327534, 1.85965675, 1.41435024, 2.4473795 , 2.22478731, 2.36732106, 1.55086739, 1.74647516, 2.02080924, 0.16702836, 2.3573375 , 2.12666702, 0.37478988, 2.86760093, 2.64067115, 2.70584321
noiseAmp = 0.5
refAmp = 0.5
[injection]
freq = 59.98
hAmp = 0.0
phi= 0.0
[primary]
primNoiseAmp = 2E-17
clutterAmpInPrim = 7E-18
#hAmp = 1e-19
#weightSines = 0.2, 0.6,0.1,0.1,0.1,0.1,0.1,0.1,0.1, 0.8, 0.6, 0.2
#freqSines = -0.5, -0.3,-0.25,-0.2,-0.19,-0.18,-0.17,-0.16,-0.15,0.00, 0.3, 0.5
#phaseSines = 0.1, 0.4,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.9, 0.4, 0.3
# 1.0 #, 0.09, 0.2, 0.04, 0.35
# 0.00,-0.0002,0.002,0.0001
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