fixing up python script for clag

scripts/analyze_lightcurve.py

@@ -16,122 +16,97 @@ except getopt.GetoptError:
     print 'analyze_lightcure.py <reference curve> <compared curve>'
     sys.exit(2)
 
-## load the first light curve
-lc1_table = np.loadtxt(args[0],skiprows=1)
-
-# works if first two entries represent minimum spacing, from example
-# dt = lc1_table[1,0] - lc1_table[0, 0]
-
 # Time resolution determined from inspection and testing. This script
 # does not expect evenly spaced data in time.
-dt = 0.1
+dt = 0.01
-
-
-# _ = plot(lc1_table[:,0], lc1_table[:,1])
-# _ = plot(lc1_table[:,0], lc1_table[:,3])
-
-
-# Split the light curve into segments #
-#seg_length = 256
-#index = np.arange(len(data)).reshape((-1, seg_length))
-
-# For now, instead of splitting up the curves, the program will assume 
-# that the data list is shorter than 256 elemements. so,
-
-index = np.arange(len(lc1_table)).reshape(-1,len(lc1_table))
-
-lc1_time = [lc1_table[i, 0] for i in index]
-lc1_strength  = [lc1_table[i, 1] for i in index]
-lc1_strength_err = [lc1_table[i, 2] for i in index]
 
 #### Get the psd for the first light curve ####
 
 # These bin values determined summer 2016 for STORM III optical/UV lightcurves
-# fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029, 0.25819945, 0.40020915, 0.62032418])
+fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029, 0.25819945, 0.40020915, 0.62032418])
 
-# using utilities to set up frequency bins #
+#A general rules for fqL is as follows:
-seg_length = 256;
+#
-fqL = np.logspace(np.log10(1.1/seg_length),np.log10(.5/dt),7)
+#    define f1, f2 as the two extreme frequencies allowed by the data. i.e.
-fqL = np.concatenate(([0.5/seg_length], fqL))
+#    f1=1/T with T being the length of observation in time units, and
+#    f2=0.5/Δt
+#
+#	The frequency limits where the psd/lag can be constrained are about
+#	~0.9f1−0.5f2. The 0.9 factor doesn't depend on the data much, but
+#	values in the range ~[0.9-1.1] are ok. The factor in front of f2
+#	depends on the quality of the data, for low qualily data, use ~0.1--
+#	0.2, and for high quality data, increase it up to 0.9−−1.
+#
+#	Always include two dummy bins at the low and high frequencies and
+#	ignore them. The first and last bins are generally biased, So I suggest
+#	using the first bin as 0.5f1−0.9f1 (or whatever value you use instead
+#	of 0.9f1, see second point above), and the last bin should be
+#	0.5f2−2∗f2 (or whatever value instead of 0.5f2, see second point
+#	above). So the frequency bins should be something like:
+#	[0.5f1,0.9f1,...,0.5f2,2f2], the bins in between can be devided as
+#	desired.
+#
+#fqd is the bin center
+#
+# If lightcurves need to be split:
+# seg_length = 256;
+# fqL = np.logspace(np.log10(1.1/seg_length),np.log10(.5/dt),7)
+# fqL = np.concatenate(([0.5/seg_length], fqL))
+#
 
 nfq = len(fqL) - 1
+fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.)
+
+
+## load the first light curve
+lc1_time, lc1_strength , lc1_strength_err = np.loadtxt(args[0],skiprows=1)
+
+# for pylab: errorbar(t1,l1,yerr=l1e,fmt='o')
+
+# Used throughout
+initial_args = np.ones(nfq)
 
 
 ## initialize the psd class for multiple light curves ##
-P1  = clag.clag('psd', lc1_time, lc1_strength, lc1_strength_err, dt, fqL)
+P1  = clag.clag('psd10r', [lc1_time], [lc1_strength], [lc1_strength_err], dt, fqL)
-
+ref_psd, ref_psd_err = clag.optimize(P1, initial_args)
-
+ref_psd, ref_psd_err = clag.errors(P1, ref_psd, ref_psd_err)
-
-
-
-## initial parameters, start with ones
-inpars = np.ones(nfq)
-
-## print the loglikelihood for the input values ##
-P1.logLikelihood(inpars)
-
-
-
-## Now do the fitting and find the best fit psd values at the given frequency bins ##
-ref_psd, ref_psd_err = clag.optimize(P1, inpars)
-
-
-
 
 ## plot ##
-fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.)
+#xscale('log'); ylim(-4,2)
-
+#errorbar(fqd, ref_psd, yerr=ref_psd_err, fmt='o', ms=10)
-#loglog(fqd, 0.1*fqd**(-1.5), label='input psd')
-#errorbar(fqd[1:-1], ref_psd[1:-1], yerr=ref_psd_err[1:-1], fmt='o', ms=10, label='fit')
-
-# load second lightcurve
-
-# This would work if both curves are in same file
-# lc2_strength  = [lc1_table[i, 3] for i in index]
-# lc2_strength_err = [lc1_table[i, 4] for i in index]
-
-# But, they aren't, so,
 
 # Load second light curve
-lc2_table = np.loadtxt(args[1],skiprows=1)
+lc2_time, lc2_strength, lc2_strength_err = np.loadtxt(args[1],skiprows=1)
-
+P2  = clag.clag('psd10r', [lc2_time], [lc2_strength], [lc2_strength_err], dt, fqL)
-index = np.arange(len(lc2_table)).reshape(-1,len(lc2_table))
+echo_psd, echo_psd_err = clag.optimize(P2, initial_args)
-
+echo_psd, echo_psd_err = clag.errors(P2, echo_psd, echo_psd_err)
-lc2_time = [lc2_table[i, 0] for i in index]
-lc2_strength  = [lc2_table[i, 1] for i in index]
-lc2_strength_err = [lc2_table[i, 2] for i in index]
-
-
-
-## Now do the second light curve
-
-P2  = clag.clag('psd', lc2_time, lc2_strength, lc2_strength_err, dt, fqL)
-
-echo_psd, echo_psd_err = clag.optimize(P2, inpars)
-
-
-
 
 
 
 ### Now the cross spectrum ###
 ### We also give it the calculated psd values as input ###
-Cx = clag.clag('cxd', 
+Cx = clag.clag('cxd10r',
-               [list(i) for i in zip(lc1_time,lc1_time)], 
+				[[lc1_time,lc1_time]], 
-               [list(i) for i in zip(lc1_strength,lc2_strength)],
+            	[[lc1_strength,lc2_strengt]],
-               [list(i) for i in zip(lc1_strength_err,lc2_strength_err)], 
+            	[[lc1_strength_err,lc2_strength_err]], 
                dt, fqL, ref_psd, echo_psd)
 
-inpars = np.concatenate( (0.3*(ref_psd*echo_psd)**0.5, ref_psd*0+1.) )
+#Cx_vals = np.concatenate( (0.3*(ref_psd*echo_psd)**0.5, ref_psd*0+1.) )
-p, pe = clag.optimize(Cx, inpars)
+Cx_vals = np.concatenate( ((ref_psd+echo_psd)*0.5-0.3,ref_psd*0+0.1) )
-phi, phie = p[nfq:], pe[nfq:]
+Cx_vals, Cx_err = clag.optimize(Cx, Cx_vals)
+
+#?????? %autoreload
+Cx_vals, Cx_err = clas.errors(Cx,Cx_vals,Cx_err)
+
+phi, phie = Cx_vals[nfq:], Cx_err[nfq:]
 lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd)
-cx, cxe   = p[:nfq], pe[:nfq]
+cross_spectrm, cross_spectrm_err = Cx_vals[:nfq], Cx_err[:nfq]
 
 np.savetxt("freq.out",fqL.reshape((-1,len(fqL))))
 np.savetxt("ref_psd.out",[ref_psd,ref_psd_err])
 np.savetxt("echo_psd.out",[echo_psd,echo_psd_err])
-np.savetxt("crsspctrm.out",[cx,cxe])
+np.savetxt("crsspctrm.out",[cross_spectrm,cross_spectrm_err])
 np.savetxt("timelag.out",[lag,lage])