From d4a11cb367dda4987535c62ae6c26b95df351666 Mon Sep 17 00:00:00 2001 From: caes Date: Thu, 22 Jun 2017 08:11:27 -0400 Subject: [PATCH] more scripts --- scripts/plot.sh | 32 +++++ scripts/psdlag_4bin.py | 45 ++++--- scripts/psdlag_7bin.py | 102 ++++++++++++++++ scripts/psdlag_8bin.py | 102 ++++++++++++++++ scripts/templates/lag_overlay_atlas.gp | 159 +++++++++++++++++++++++++ 5 files changed, 415 insertions(+), 25 deletions(-) create mode 100644 scripts/psdlag_7bin.py create mode 100755 scripts/psdlag_8bin.py create mode 100644 scripts/templates/lag_overlay_atlas.gp diff --git a/scripts/plot.sh b/scripts/plot.sh index ff1e8ff..c83af96 100755 --- a/scripts/plot.sh +++ b/scripts/plot.sh @@ -45,6 +45,38 @@ case $1 in rm $gnuplot_file ;; + "twolags"|"2lags") + gnuplot_file=lag_overlay_atlas.gp + gnuplot_input=$(cat scripts/templates/${gnuplot_file}|perl -pe 's|\n|␤|g') + for tabfile in data/tables/${2}/lag_*.tab; + do + ref_band_extracted=$(basename $tabfile| + sed 's|lag_\([0-9]\{4\}A\)_[0-9]\{4\}A.tab|\1|') + echo_band=$(basename $tabfile| + sed 's|lag_[0-9]\{4\}A_\([0-9]\{4\}A\).tab|\1|') + if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi + gnuplot_input_edit=$(echo "$gnuplot_input"| + sed "s|%FILE%|$tabfile|"| + sed "s|%LABEL%|$echo_band|") + gnuplot_input="${gnuplot_input_edit}" + done + for tabfile in data/tables/${3}/lag_*.tab; + do + ref_band_extracted=$(basename $tabfile| + sed 's|lag_\([0-9]\{4\}A\)_[0-9]\{4\}A.tab|\1|') + echo_band=$(basename $tabfile| + sed 's|lag_[0-9]\{4\}A_\([0-9]\{4\}A\).tab|\1|') + if [[ "$echo_band" == "$ref_band" ]] ; then continue; fi + gnuplot_input_edit=$(echo "$gnuplot_input"| + sed "s|%FILE%|$tabfile|"| + sed "s|%LABEL%|$echo_band|") + gnuplot_input="${gnuplot_input_edit}" + done + echo "$gnuplot_input"|perl -pe 's|␤|\n|g' > ${gnuplot_file} + gnuplot $gnuplot_file + rm $gnuplot_file + ;; + "tophat"|"th") mkdir -p data/tables/ scripts/tophat_fft.pl diff --git a/scripts/psdlag_4bin.py b/scripts/psdlag_4bin.py index b2f44dd..6dfc843 100755 --- a/scripts/psdlag_4bin.py +++ b/scripts/psdlag_4bin.py @@ -22,11 +22,16 @@ t1, l1, l1e = np.loadtxt(ref_file).T # errorbar(t1, l1, yerr=l1e, fmt='o') - -fqL = np.array([0.0049999999, 0.044733049, 0.10747115, - 0.25819945, 0.62032418]) -#fqL = np.array([0.0049999999, 0.018619375, 0.069336227, 0.10747115, 0.62032418]) +#A #fqL = np.logspace(np.log10(0.005),np.log10(0.6),5) + +#B +#fqL = np.array([0.0049999999, 0.018619375, 0.069336227, +# 0.10747115, 0.62032418]) + +#C +fqL = np.array([0.0049999999, 0.044733049, 0.10747115, + 0.22, 0.56]) nfq = len(fqL) - 1 fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.) @@ -40,8 +45,6 @@ p1, p1e = clag.optimize(P1, p1) p1, p1e = clag.errors(P1, p1, p1e) -# xscale('log'); ylim(-4,2) -# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="black") ref_psd = p1 ref_psd_err = p1e @@ -49,8 +52,6 @@ ref_psd_err = p1e t2, l2, l2e = np.loadtxt(echo_file).T -# errorbar(t1, l1, yerr=l1e, fmt='o', color="green") -# errorbar(t2, l2, yerr=l2e, fmt='o', color="black") P2 = clag.clag('psd10r', [t2], [l2], [l2e], dt, fqL) @@ -59,45 +60,39 @@ p2, p2e = clag.optimize(P2, p2) p2, p2e = clag.errors(P2, p2, p2e) -# xscale('log'); ylim(-6,2) -# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="green") -# errorbar(fqd, p2, yerr=p2e, fmt='o', ms=10, color="black") echo_psd = p2 echo_psd_err = p2e -Cx = clag.clag('cxd10r', [[t1,t2]], [[l1,l2]], [[l1e,l2e]], dt, fqL, p1, p2) -p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) # a good starting point generally +Cx = clag.clag('cxd10r', + [[t1,t2]], + [[l1,l2]], + [[l1e,l2e]], + dt, fqL, p1, p2) +# a good starting point generally +p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) p, pe = clag.optimize(Cx, p) phi, phie = p[nfq:], pe[nfq:] lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd) cx, cxe = p[:nfq], pe[:nfq] -cross_spectrm = cx -cross_spectrm_err = cxe +crs_spectrm = cx +crs_spectrm_err = cxe -# xscale('log'); ylim(-2,1) -# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") s, loc, scale = lognorm.fit(lag,loc=.01) -# xscale('log'); ylim(-4,1.5) -# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") -##plot(fqd,norm.pdf(fqd,mu,sigma)) -#plot(fqd,lognorm.pdf(fqd,s,loc,scale)) -# mu,sigma -#plot(ifft(lag)) 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",[cross_spectrm,cross_spectrm_err]) -np.savetxt("timelag.out",[lag,lage]) +np.savetxt("crsspctrm.out",[crs_spectrm,crs_spectrm_err]) +np.savetxt("lag.out",[lag,lage]) diff --git a/scripts/psdlag_7bin.py b/scripts/psdlag_7bin.py new file mode 100644 index 0000000..367621f --- /dev/null +++ b/scripts/psdlag_7bin.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python +import numpy as np +from scipy.stats import norm +from scipy.stats import lognorm +import sys +import getopt +sys.path.insert(1,"/usr/local/science/clag-agn/data/") +import clag +import matplotlib +# %pylab inline + + +#ref_file="data/lc/1367A.lc" +#echo_file="data/lc/2246A.lc" + +ref_file = str(sys.argv[1]) +echo_file = str(sys.argv[2]) + + +dt = 0.01 +t1, l1, l1e = np.loadtxt(ref_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o') + + + +fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029, + 0.25819945, 0.40020915]) +#fqL = np.logspace(np.log10(0.005),np.log10(0.4),8) +nfq = len(fqL) - 1 +fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.) + + + + + +P1 = clag.clag('psd10r', [t1], [l1], [l1e], dt, fqL) +p1 = np.ones(nfq) +p1, p1e = clag.optimize(P1, p1) + +p1, p1e = clag.errors(P1, p1, p1e) + +# xscale('log'); ylim(-4,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="black") + +ref_psd = p1 +ref_psd_err = p1e + + + +t2, l2, l2e = np.loadtxt(echo_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o', color="green") +# errorbar(t2, l2, yerr=l2e, fmt='o', color="black") + + +P2 = clag.clag('psd10r', [t2], [l2], [l2e], dt, fqL) +p2 = np.ones(nfq) +p2, p2e = clag.optimize(P2, p2) + +p2, p2e = clag.errors(P2, p2, p2e) + +# xscale('log'); ylim(-6,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="green") +# errorbar(fqd, p2, yerr=p2e, fmt='o', ms=10, color="black") +echo_psd = p2 +echo_psd_err = p2e + + + + +Cx = clag.clag('cxd10r', [[t1,t2]], [[l1,l2]], [[l1e,l2e]], dt, fqL, p1, p2) +p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) # a good starting point generally +p, pe = clag.optimize(Cx, p) + +phi, phie = p[nfq:], pe[nfq:] +lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd) +cx, cxe = p[:nfq], pe[:nfq] + +cross_spectrm = cx +cross_spectrm_err = cxe + +# xscale('log'); ylim(-2,1) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") + + + + +s, loc, scale = lognorm.fit(lag,loc=.01) +# xscale('log'); ylim(-4,1.5) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") +##plot(fqd,norm.pdf(fqd,mu,sigma)) +#plot(fqd,lognorm.pdf(fqd,s,loc,scale)) +# mu,sigma + + + +#plot(ifft(lag)) + +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",[cross_spectrm,cross_spectrm_err]) +np.savetxt("lag.out",[lag,lage]) diff --git a/scripts/psdlag_8bin.py b/scripts/psdlag_8bin.py new file mode 100755 index 0000000..18cf6f5 --- /dev/null +++ b/scripts/psdlag_8bin.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python +import numpy as np +from scipy.stats import norm +from scipy.stats import lognorm +import sys +import getopt +sys.path.insert(1,"/usr/local/science/clag-agn/data/") +import clag +import matplotlib +# %pylab inline + + +#ref_file="data/lc/1367A.lc" +#echo_file="data/lc/2246A.lc" + +ref_file = str(sys.argv[1]) +echo_file = str(sys.argv[2]) + + +dt = 0.01 +t1, l1, l1e = np.loadtxt(ref_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o') + + + +fqL = np.array([0.0049999999, 0.018619375, 0.044733049, 0.069336227, 0.10747115, 0.16658029, + 0.25819945, 0.40020915, 0.62032418]) +#fqL = np.logspace(np.log10(0.005),np.log10(0.6),9) +nfq = len(fqL) - 1 +fqd = 10**(np.log10( (fqL[:-1]*fqL[1:]) )/2.) + + + + + +P1 = clag.clag('psd10r', [t1], [l1], [l1e], dt, fqL) +p1 = np.ones(nfq) +p1, p1e = clag.optimize(P1, p1) + +p1, p1e = clag.errors(P1, p1, p1e) + +# xscale('log'); ylim(-4,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="black") + +ref_psd = p1 +ref_psd_err = p1e + + + +t2, l2, l2e = np.loadtxt(echo_file).T +# errorbar(t1, l1, yerr=l1e, fmt='o', color="green") +# errorbar(t2, l2, yerr=l2e, fmt='o', color="black") + + +P2 = clag.clag('psd10r', [t2], [l2], [l2e], dt, fqL) +p2 = np.ones(nfq) +p2, p2e = clag.optimize(P2, p2) + +p2, p2e = clag.errors(P2, p2, p2e) + +# xscale('log'); ylim(-6,2) +# errorbar(fqd, p1, yerr=p1e, fmt='o', ms=10, color="green") +# errorbar(fqd, p2, yerr=p2e, fmt='o', ms=10, color="black") +echo_psd = p2 +echo_psd_err = p2e + + + + +Cx = clag.clag('cxd10r', [[t1,t2]], [[l1,l2]], [[l1e,l2e]], dt, fqL, p1, p2) +p = np.concatenate( ((p1+p2)*0.5-0.3,p1*0+0.1) ) # a good starting point generally +p, pe = clag.optimize(Cx, p) + +phi, phie = p[nfq:], pe[nfq:] +lag, lage = phi/(2*np.pi*fqd), phie/(2*np.pi*fqd) +cx, cxe = p[:nfq], pe[:nfq] + +cross_spectrm = cx +cross_spectrm_err = cxe + +# xscale('log'); ylim(-2,1) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") + + + + +s, loc, scale = lognorm.fit(lag,loc=.01) +# xscale('log'); ylim(-4,1.5) +# errorbar(fqd, lag, yerr=lage, fmt='o', ms=10,color="black") +##plot(fqd,norm.pdf(fqd,mu,sigma)) +#plot(fqd,lognorm.pdf(fqd,s,loc,scale)) +# mu,sigma + + + +#plot(ifft(lag)) + +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",[cross_spectrm,cross_spectrm_err]) +np.savetxt("lag.out",[lag,lage]) diff --git a/scripts/templates/lag_overlay_atlas.gp b/scripts/templates/lag_overlay_atlas.gp new file mode 100644 index 0000000..266ac69 --- /dev/null +++ b/scripts/templates/lag_overlay_atlas.gp @@ -0,0 +1,159 @@ +set terminal pdf size 3.5,5 +set output "lag_atlas.pdf" +set termopt enhanced + +set macros + + # Placement of the a,b,c,d labels in the graphs + POS = "at graph 0.61,0.85 font 'Times,12'" + + # x- and ytics for each row resp. column + NOXNUMS = "unset xlabel;\ + set format x ''" + XNUMS = "unset xlabel;\ + set format x '10^{%+3T}'" + XLABEL = "set xlabel 'Temporal Frequency [days^{-1}]' font 'Times' off 0,1;\ + set format x '10^{%+3T}'" + YLABEL = "set ylabel 'Lag [days]' font 'Times' offset 1,2.5;\ + set format y '%.0t'" + NOYNUMS = "set format y ''; unset ylabel" + YNUMS = "set format y '%.0t'" + + + VSET_1 = "set tmargin at screen 0.97; set bmargin at screen 0.825" + VSET_2 = "set tmargin at screen 0.825; set bmargin at screen 0.68" + VSET_3 = "set tmargin at screen 0.68; set bmargin at screen 0.535" + VSET_4 = "set tmargin at screen 0.535; set bmargin at screen 0.39" + VSET_5 = "set tmargin at screen 0.39; set bmargin at screen 0.245" + VSET_6 = "set tmargin at screen 0.245; set bmargin at screen 0.10" + + HSET_1 = "set lmargin at screen 0.15; set rmargin at screen 0.42" + HSET_2 = "set lmargin at screen 0.42; set rmargin at screen 0.69" + HSET_3 = "set lmargin at screen 0.69; set rmargin at screen 0.96" + +unset key +set logscale x +set xtics auto font 'Times,9' offset 0,.5 +set ytics (-1,0,1,2,3,4,5) font 'Times,9' +set ytics add ('' -1.5 1,'' -0.5 1,'' 0.5 1,'' 1.5 1,'' 2.5 1,'' 3.5 1,'' 4.5 1,'' 5.5 1) +set xrange [0.005:0.620]; +set yrange [-2:6] + +# Draw line at origin +set arrow from 0.005,0 to 0.620,0 nohead lt 3 lc rgb 'black' +set pointsize 0 + +set multiplot layout 6,3 rowsfirst + + + # --- GRAPH a + @VSET_1; @HSET_1 + @NOXNUMS; @YNUMS + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH b + @NOXNUMS; @NOYNUMS + @VSET_1; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + #-- GRAPH c + @NOXNUMS; @NOYNUMS + @VSET_1; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH d + @NOXNUMS; @YNUMS + @VSET_2; @HSET_1 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_2; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_2; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @YNUMS + @VSET_3; @HSET_1 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_3; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_3; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @YLABEL + @VSET_4; @HSET_1 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_4; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_4; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @YNUMS + @VSET_5; @HSET_1 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_5; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @NOXNUMS; @NOYNUMS + @VSET_5; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @XNUMS; @YNUMS + @VSET_6; @HSET_1 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @XLABEL; @NOYNUMS + @VSET_6; @HSET_2 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + # --- GRAPH a + @XNUMS; @NOYNUMS + @VSET_6; @HSET_3 + set label 1 '%LABEL%' @POS + plot '%FILE%' using 1:2:($2-$4):($2+$4) with yerrorbars pt 7 ps .3 lt rgb "black" + + +unset multiplot +