From c24d7da5940bef92b408103e36ad32188d6e3d58 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Sun, 19 May 2019 21:07:10 -0400 Subject: [PATCH 01/14] Initial work on ACAImage support for flickering pixels --- chandra_aca/aca_image.py | 105 ++++++++++++++++++++++++++- chandra_aca/data/flicker_cdf.fits.gz | Bin 0 -> 6370 bytes 2 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 chandra_aca/data/flicker_cdf.fits.gz diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index c9a5609..005e8d2 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -3,6 +3,7 @@ from math import floor from itertools import count, chain from copy import deepcopy +from pathlib import Path import six from six.moves import zip @@ -309,6 +310,108 @@ def col0(self): def col0(self, value): self.meta['IMGCOL0'] = np.int64(value) + @classmethod + def _read_flicker_cdfs(cls): + """Read flickering pixel model cumulative distribution functions + and associated metadata. Set up class variables accordingly. + + """ + from astropy.io import fits + + filename = Path(__file__).parent / 'data' / 'flicker_cdf.fits.gz' + with fits.open(filename) as hdus: + hdu = hdus[0] + hdr = hdu.header + + # Get the main data, which is an n_cdf * n_cdf_x array. Each row corresponds + # to the CDF for a particular bin range in e-/sec, e.g. 200 to 300 e-/sec. + # CDF will go from 0.0 to 1.0 + cls.flicker_cdfs = hdu.data + + # CDF_x is the x-value of the distribution, namely the log-amplitude change + # in pixel value due to a flicker event. + cls.flicker_cdf_x = np.linspace(hdr['cdf_x0'], hdr['cdf_x1'], hdr['n_cdf_x']) + + # CDF bin range (e-/sec) for each for in flicker_cdfs. + cdf_bins = [] + for ii in range(hdr['n_bin']): + cdf_bins.append(hdr[f'cdf_bin{ii}']) + cls.flicker_cdf_bins = np.array(cdf_bins) + + def flicker_init(self, mean_flicker_time=10000, flicker_scale=1.0): + """Initialize instance variables to allow for flickering pixel updates. + + :param mean_flicker_time: mean flickering time (sec, default=10000) + :param flicker_scale: multiplicative factor beyond model default for + flickering amplitude (default=1.0) + """ + if not hasattr(self, 'flicker_cdf_bins'): + self._read_flicker_cdfs() + + self.mean_flicker_time = mean_flicker_time + self.flicker_scale = flicker_scale + + # Make a flattened view of the image for easier update processing. + # Also store the initial pixel values, since flicker updates are + # relative to the initial value, not the current value (which would + # end up random-walking). + self.flicker_vals = self.view(np.ndarray).ravel() + self.flicker_vals0 = self.flicker_vals.copy() + self.flicker_n_vals = len(self.flicker_vals) + + # Get the index to the CDFs which is appropriate for each pixel + # based on its initial value. + self.flicker_cdf_idxs = np.searchsorted(self.flicker_cdf_bins, + self.flicker_vals0) - 1 + + # Create an array of time (secs) until next flicker for each pixel + # This is drawing from an exponential distribution. For the initial + # time assume the flickering is randomly phased within that interval. + phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) + rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) + t_flicker = -np.log(1.0 - rand_unifs) * mean_flicker_time + self.flicker_times = t_flicker * phase + + # Pixels where self.flicker_cdf_idxs == 0 have val < 50 (no CDF) and are + # modeled as not flickering. Make a mask to indicate which ones flicker. + self.flicker_mask = self.flicker_cdf_idxs >= 0 + + def flicker_update(self, dt): + """ + Propagate the image forward by ``dt`` seconds and update any pixels + that have flickered during that interval. + + TO DO: use numba for this once numba with np.interp is available in Ska3. + (E.g. 0.43 has it). + + :param dt: time (secs) to propagate image + """ + if not hasattr(self, 'flicker_times'): + self.flicker_init() + + self.flicker_times[self.flicker_mask] -= dt + + # When flicker_times < 0 that means a flicker occurs + idxs = np.where(self.flicker_times < 0)[0] + + # Random uniform used for (1) distribution of flickering amplitude + # via the CDFs and (2) distribution of time to next flicker. + rand_ampls = np.random.uniform(0.0, 1.0, size=len(idxs)) + rand_times = np.random.uniform(0.0, 1.0, size=len(idxs)) + + for idx, rand_time, rand_ampl in zip(idxs, rand_times, rand_ampls): + # Determine the new value after flickering and set in array view. + # First get the right CDF from the list of CDFs based on the pixel value. + cdf_idx = self.flicker_cdf_idxs[idx] + y = np.interp(fp=self.flicker_cdf_x, + xp=self.flicker_cdfs[cdf_idx], + x=rand_ampl) + val = 10 ** (np.log10(self.flicker_vals0[idx]) + y) + self.flicker_vals[idx] = val + + # Get the new time before next flicker + t_flicker = -np.log(1.0 - rand_time) * self.mean_flicker_time + self.flicker_times[idx] = t_flicker def _prep_6x6(img, bgd=None): @@ -410,7 +513,7 @@ def __init__(self, filename=None): if filename is None: filename = os.path.join(os.path.dirname(__file__), 'data', 'aca_psf_lib.dat') dat = Table.read(filename, format='ascii.basic', guess=False) - self.dat = dat + self.dat = dat.as_array() # Sub-pixel grid spacing in pixels. This assumes the sub-pixels are # all the same size and square, which is indeed the case. diff --git a/chandra_aca/data/flicker_cdf.fits.gz b/chandra_aca/data/flicker_cdf.fits.gz new file mode 100644 index 0000000000000000000000000000000000000000..cc22a2174902a51214734c0d80444775f9c30e6d GIT binary patch literal 6370 zcmV<87#-&yiwFoe=;2%f|7L7yV{2t{Ut?ruE@o+Ta{vGU0RR8&)Ikb@KorLDtGo+X zj;YncMWw<(NuUdNGNpnGXJWKPw@-vty&;@&6Zk*CkB@m`N-x)Q1P}mH25YzrUJKcF z&=y$Y~0Sea_(VXM3w8i25Tp zZNnZN5 z<0ZwBCW`kXn)?RP#yn5oBXizk@EpW*47?5EZ45jY@mvGX3$SzOdp_d%2HqC&b_U)q zz&;m!Zy$VPTHg5cTUuIL`h(s800960>{#hPm0i>}Wv+zE^rVD}6p7}uOd(Q+L=sV` zM9EYV3MG=1hzuziQc;F9nF*CK8In1}!MSCY>0Re^e|g{MA9z1J_m}Ip_g;JL&F{MQ zy4GPLKT4Bm=1b=^g3`~xVaoUvMd;}cDrLbszM1dgAaMPi2SM9wUq8-|f*8nZs7cXa zZOqg4I#$jgT_Bn()ZheioSgkxMGmm#ylJ}F_*PIhe?DNz`5GwmiX4Z(ZUlAq`9a0h zL!dcow9YlR2KMugHvgjB1_SeTgE6yD!C2dKT|1uyn3V>-_o6p+H-J4Q_`&%H`NR0e3Y{DGof#hx>Q;1U~Yegvbl513|}PAm(9EnsVw& zh<`7XR`YlTBvkx<(0>U74ks6l{ z5Du#y?_bvh!mlOH`rH#BA`+Uo2{|B=4sV#8ybMHUZIJJ?SRnEfRxB;!2jcT1B{7jS zAnGsQ4~!52;wLAyYw0f_y3B3**?55HyB6kFQ3AxE{8~rTH$V`A-}0S#f%tvSA^mGN z5Mx@Q>+dW8F`*!5-&_U6q{XFf&fP#v%@%QKrUUVZ&1S^$9AbXf$G|!uX2RF(F@6Tb z>_w3X_xC`|Z9QWhCXP6#&B^frh>;__Zwj zi7?{m1*zmx7-srkAEm-U!n*01DW#c2Dv6^tO%-dQkTNOmwZhAFG3C$Sd0_HWx>CmN z4&3B;@w?Krt!+<$UzX#^_){?uA?H_VM=sPT~_X&L8oK4EN6@mc1_+brqfZL~=GFS^z;hvc6mG?;p;eoifQS`Djh%Msg zSaw_oo}B)|pUaNAiQb~u>h~F5j)f-$=A40y&?b)`jggR3!EyiU!~w{iQ**0Z_Z^DX z7x}Cz;)YMI$4eo{9lp%LYy>S7YJL0CqeVACV?&{~rNK#PZy&!tIb#pq-`D0Cb-sY! z02%XF3(C;H``XcqoA&_K#rNl!=2{r=JoBdUYz+(=ABl|Zd;~*FUo2nt!vuz&7Sk+s zqk$%V)>zi;3DCT|-Y9HQ1lpHFHtSbc1HrMkx;gqF5IgG&E_kv6VU@@;nwkZK_lQh8 zlNu1=HeKPf@<1f3KB+6$1|qxgrXm;er&w+*Ww#a(l~q5znVW(59{fFYG7fpKZ2Z2> z9eFSG_(TT{h~9un8P6*~47}MTDn0}RE!pban)S&0<+dxgt_5PG+_CbWEb{*Z_lb}r zK#aez2(_L@{AG-^npFbC!p;-NzMu{^ zh7a~DA|BfL=*lI;a9uZ-NJPd3|D}%^^@;MEV!f2vM-xWT#rc%Ex$&RrXTvEAX5)6n zg6Dwg=%0%Nt&+gZ?XhS4H6JkN{9#$}p9L1XOJ!U~`hj)%bg+P2J+KyEGW>@*AJ{HT zuavHT0n5mEez}7R>;ij1KwnaFg^LFWlItthZTHQy^MaK>+=or3bwD`ejrnkp5QzT0 z4@i=68>f_&*SyRIh3?2P)jw0PJu4zHK}`WD&s`tCQ}F|}+_%qetG9s`v*}Kr^b??G zqq$70ZUY#Q@rRG)I~iWe26Hlx6*+IWt$kf=VGKBud7K|K9dEu>3YR*3RrRFLfsc#d z`FI!s{|u97J?#Z>>*zUqk>-PNPqy4Y^OO-paWG{*R-1(QLcfL5rD2dX)A-bG|38rK z&R zbZk<0YE@>3U+S?NlNM8nMfxAF_DbYJ2piFPGbAi9N|96mk=f?6rN;jlgOUDTE* zcM9@s3aS%D8{xjOZRIevy92O6@MgE7n);D3DnNeCCd$$Mi@BfdOm2uKd&OzLprFH@DFc z(8dJyw!j^2NZo`e)xIwS>(_$&^G$q)h}|E!l->WOkMGJKUCx>llpYd?QI`C}BE*1WZ9404dj)X0EIIQ*#szrD z_$7j$w;O(}0Y1jQd{|~y%rcf74*ZP$XF(uU;&k4~UJzvD5djfSX*S2%b`anF$3Dk8 z2_&}np7!b929o6ZI>tFAuz_)oJIIpfY$eyX>$Yt>>>mre$oV}!8L@E#Cqa#oM;Wxq z_5F|UY!^r?0)z7U?46~yVB+w;x?#oRtN?-bvVG z=q-i5d#-Vv+4Vpz$=l`nq67vHZ{@AKtN}yUk{mS?sX(LbuTPBL2ei_Fpup%KK&%g! zt@K6TW2cqk9f3aSZl~Dsm~0?k^#+e*aA6$?xv|s({Y`Cqm8>u7=10jp_IUJ7oiq!< zMATJZ`B%z?Heh~v`_ ztCt|}nFap@Y9R8;AJOnYT$Qok;~S!aGQ7e(^XhN6^xi~#UC|?#gjgA3xPbbYryAu| zWBpvXL)-th6S2%)+Whyw^ii~^A@vxu8>NoKmVlzGJl7;BU5v3KrT3Fw_*T^~ltIS) z2FmbLv8>-#>XcE&`UT3k?cGOZ@$r;NGVafM&Y-xCKa^QAe!kz!&0PY0+24J@T<;;D zq;LdS|LzAiWhH6;nj5f)ktYt=!-qw^oU(y~k?$*TPNwovm;D0nzj=Yj)Rx`&wE*z` z%@0!KJ_SZzGuY0E?*@v&!M9>N^+A=K*Cyxnd55OLFKvTE*Y^A6?+pf1a{f3ueoZUUOsL z-tL+#h`vG6(e!>ICA@yrjowd)N4nsB<~x(k$Ab|U6jr;@eUj+guE_+%Ex{Rg1rhgI z`EK+@Jg#I|I!dn>ToP7z{z7oe@0nV}jJ|loTZm0#`YD?K(noFR#=|nNQh4m1^B zJ*V}8GQt>trc5fi+qb;Wr_3w-zK1lzCZ6>cTq4{$h|3lRjct}Lsw6Vl;MXQIRolPHMvzT&K+stE};bK_r`7r{$KXZBqWYayp_ zfOE9128!7P=_^pR{g(gRrToxPQ{=gFl^S#qeGad)U7+WdZpKN!8tAWCZokr12L@MJ zOlnFs!;s5*li_f0pc!ovzbfz$Xgy=W9Y65>)L&J{@(|x`A6Y-3Qh>(ZaApr4x4`oL^1f#`U~d1)`= zeWC6m^n-KnKOJ^=LhN)clih~8;ga0ZI*YiuHuWz0vH2&hOfRucEv)P}E=OOvV57N1 zegN_H6R`(w|I)|s#RFBZ#G8PLM9xgNZ8y{20v}^s2ZG_&zA<*?AjFu@1F_nh1Fn-< zu$r;{9VBj;rX;D?fMlb3Zp%h(Shr!*w`GFGu)#yPfA1UxWIwabd;EF_n?#2N(ki!t zys|^sg`+!Pu_?yU6`&jJOKWBG`<}86zjB{AQ zY4o3NGp-D93maQf)SLt!jPtI*wLQDEm=X)&Hks#sXt0uJ<_J7G>bbWhr462~ZsFM7 z`WUjd_IYeP9|r~Dmo{*f9EUG024rNsfQ&?}YY?6mhZP{UGw zTW#fr!SsW@FVR=fxI^Xd6xaZ*gjMW)H~Jg-tm+j^=;vKIdpRPJHxZ>T_GVxooPQv1 zdkyxxjKt`Z+U`jUd|iO6?-h0nS8E{}*@@hQZ2cT{bZ=NS62G5b^ZC-Ge# ztK8{WAA@{X;MK~&dN6*kb#*8ev1)yNH^xnHI|nHXA@0pp3BbO2;wtBOGxqNjAB<8r zV*Q!yE>|kWcX>+d;gSUGzo(v^?e)ceXL>Ou-wNaYC}rx*yhq+!Y^V=M|1*0%?zsFD zuPJSDA+>o+iFrz)lokHO62}+ zOO!LWaI(R6MjlJpQ5@kV{N)|&V&q#6yFI3~nM2RR9x^W|99h2c0(+7Z;Zq zKp)d|U|EPR4D7nEzE(LKhGdFU#BUD)4Wi;={VRahtKMqy8F{xq0h(A0fN&2G!J~nQ z^)F9m$9FZ0`o+~1@on{F#Orw=z8qbA`#G;=PjrE~<+;-iR>KGU&7HjNEzf-5)VK&A9>J*)ct@NNpd)(I2bk zus;}&+LrKo2kJv5OvY|GV$X28DAujXU_GU3tOrwyDewI;|5P)(>dqUe6CdkR;sokM zI&Xi13gTDM+9UMeshlsGF2=uWR`}b`LcGVU>9&v_tb4QJQEJ@SSImB!9ovNbpA)c5 zG(^4485Ekjjv(HV;*P|=X|D3dhmQ~FeF)XD;|#rD*ITm)>*M^*w7L!SKD~_Je?Qiv z1$NhM=jnZiMniuX*1ZMaoey)4A{NBxwaEQneSpA3j=o_^1gs@dZbb0n<#{7eWQ;># z`x_Os3x#4hT^kQxOd;pwfoaYYb$-GW>eXqZ@@q!a0|01|>;{@O6syT3G zoHGv2jQb3M8=2p|rB9!w@-l$UAIR7j1ff9<^AD1^;V$@HvJYW~2Y13M>a&yKk)${0NW&_4TF}`RxS|R& z$UMc9`Ss>C%~03IYV53a2D-&c<5E9K19iW2)RGhWF!*K9(wL7gfVM;Nw1_138{f|y zx|>!Fgz#sU_|mUH9Ft8?|5*-%bG!=A&N3iEG8K1kMt`5Up)mXw`k7ZJ1J&R8q0fpl z8184ox^gVcdkV+Jq76$q(Pw>X7LHwsb?2-56_=M-C+Y+{bm>07^s>>A>&;=w*x(F#1YIBRg9Dar k`oH@4@B8oj@B8oj@B9Dt{SN>D0RR630Jqw4#{ew=0IF=)IsgCw literal 0 HcmV?d00001 From b9e2b2d65173a01479000b7497e28720c405f1db Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Mon, 20 May 2019 08:24:18 -0400 Subject: [PATCH 02/14] Support the flicker amplitude scale factor --- chandra_aca/aca_image.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 005e8d2..10ec9eb 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -341,6 +341,12 @@ def _read_flicker_cdfs(cls): def flicker_init(self, mean_flicker_time=10000, flicker_scale=1.0): """Initialize instance variables to allow for flickering pixel updates. + The ``flicker_scale`` can be interpreted as follows: if the pixel + was going to flicker by a multiplicative factor of (1 + x), now + make it flicker by (1 + x * flicker_scale). This applies for flickers + that increase the amplitude. For flickers that make the value smaller, + then it would be 1 / (1 + x) => 1 / (1 + x * flicker_scale). + :param mean_flicker_time: mean flickering time (sec, default=10000) :param flicker_scale: multiplicative factor beyond model default for flickering amplitude (default=1.0) @@ -382,7 +388,7 @@ def flicker_update(self, dt): that have flickered during that interval. TO DO: use numba for this once numba with np.interp is available in Ska3. - (E.g. 0.43 has it). + (E.g. 0.43 has it). This will probably be substantially faster. :param dt: time (secs) to propagate image """ @@ -406,6 +412,16 @@ def flicker_update(self, dt): y = np.interp(fp=self.flicker_cdf_x, xp=self.flicker_cdfs[cdf_idx], x=rand_ampl) + + if self.flicker_scale != 1.0: + # Express the multiplicative change as (1 + x) and change + # it to be (1 + x * scale). This makes sense for positive y, + # so use abs(y) and then flip the sign back at the end. For + # negative y this is the same as doing this trick expressing the + # change as 1 / (1 + x). + dy = (10 ** np.abs(y) - 1.0) * self.flicker_scale + 1.0 + y = np.log10(dy) * np.sign(y) + val = 10 ** (np.log10(self.flicker_vals0[idx]) + y) self.flicker_vals[idx] = val From d7aaadbb944241a8f8feb9daec1484e2aeb1cd00 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 07:23:29 -0400 Subject: [PATCH 03/14] Flicking CDFs that are clipped to a factor of 2 --- chandra_aca/data/flicker_cdf.fits.gz | Bin 6370 -> 6662 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/chandra_aca/data/flicker_cdf.fits.gz b/chandra_aca/data/flicker_cdf.fits.gz index cc22a2174902a51214734c0d80444775f9c30e6d..9646b15d0b6ee1bdb9d40d597004997c9c3b243d 100644 GIT binary patch literal 6662 zcmV+h8u{fPiwFpaD&<@P|7L7yV{2t{Ut?ruE@o+Ta{vGU0RR8&)VmIXFc<~kSNUH+ zZMg{pgBMJiC^3P7ja0$Jfz}ipe0vFrlYe7qSxh(&@a43pX~s_HV*(HZk{niWmEtB< z(?TcV!14^*j*E(`3eap>er5&gEdsCixKjKD74JB7uSX3(l5Pr?fu^6+d*az#Y6bH4 zk?&0{+Pf63-lM2T34zfWJV(o<61G1>hM8bSIk$kO!< zJ#=3)FD3!Hktq=o zNrMKWNwY@MJWGSpL?z7{rP54E8a2=JIp^k~B-Fb*zw_mNpMT)}@a!+wZ|}9&z1Mx8 zb?tSn-MNBit^vaNPc>ZH|OUqRr#IC+FCBw)=+~ttr@4O|5k8WJ6 zus8!Rn_pe)y0!|hHZS~eSpG2*qbGibc-Lp;$Ss~Iv9=qD;Jl5}F^#1%eC#MM(hJkfUy1i0+9Wc{)KDc@qOwiH z7S*2|^1LEbQ8U%b&ON;awV&1Za&RS~zEjai%uWwY{8B^wjRj~f6x(~NPXev|)JR%K z5jrFyg2H5V(WT_#tof1;-63*u$6gnrr;GQzPI4IfcrW}qCG-gWqE|RAeN{0a*s;I( zhc^b7+FxIjGKIm`byI%JH)6=y`^lcCQyAW{P)vex2*VjJmGebDV`MwGpr*wujCg;T z{2(-lkzwv_Iih?RmFBqK&KZnR`#CvvZ)!05&S{pD#u!F_%PDDCD}uIte>&d*E6~&` z+yiI%fack6xJ1nrv@E6j!mkHFYaa~vWH$wUuIs5Lhbqv8L#yy874#EBiz!28(DfY` zhE%zN?z+nC&3jeQV|!N^XkP^V9mlOmAzRSfdge9juK@k$oU2{39AGTr4Yct448|sz z%V+AUz&NmdM8{Vj45e9hv);6Wp(7qJr+X0?HnCCHe@uYko1B*!^&X6vL91mx>R@Du z@BF#-F{tV9CTqS(5*(A^YwsuclVe$L0Ktj+CBt$AC!cXLz7zbT&;BEd;MATuljrbt z*8dt%lh1ikc0L?;SZJ0}d2o4tIc%BpWH{357fAm&31{Yd0o-hhjOebr;Kj7Z!bgtx zr$>Pde0?5i98ss>e_#n?K*<6DEt(B0T%!^A$Zy;^%^X37+vko)QV2!KCLyIyXy&Je<^SMxjYW> z)$5UWI{*ozK7-$nIv~k9dt6w^3@NsbeI5^ckUF?yz$@!J(#@Airnw(S#_AG#`*9Iu zQk%c>6v`rdq#&_OGy-o{(RrUY{=!>@O9q0XRVc7``FLA)GYT0$57Pb+ejIGkZ{`|D zaYBEMrRiIguD;TGTK5vlBCl0Gt-Xf|S#`~)8*K1lLvif9?ggmYb60=*?j)*>V~VrB zT||vf7ys-;PSkGf>Dsz59CgFTy~NipMB{F@_n5ODU+T^32M%vQYpKHNhtU_&@hQC` zQtcVKep+qYt-yvJ3oo<$FkSS1nH6_kVJrH%?t2#RAH=}=XI2|deZ;`D#K5uPml!N+ z{d42{5QcOZ{FLWC41L3sv`tqqY<_O7ZJ`E+2TRR5r8F^eT0v6JehEgBC43B%%P_il zXzP{>su;bww)y+?9*pMsxyN{)1#JNi8y}Mc?X>0L$$L_u`4yUr@Dd)>tMQ17_Jh9k za5MWqB%hDZ6YUV^0Nut~neVqY=&vOF?1w8rZ>hJHb?OFVPL!DYNER4dHdk~Q3xRQ} zy!w9IG#Iy}g!v{M!El?+z>c+GypWG_i(U*yfys2;oikuGmW;iYRt00!A)2R_4U8#` zVGey`P&{#YUTp55Hga;tXLW(vlNb1Er9UV+@7JT7Ux7MTAeW`o0qR=G_T3HlNxuJ6 zdD-h2!5>MdiunkR`{L*D3_nbFiRI4la3wct^q8Wc5B$IR!QGHW(M1B1&1`F z_ru<~L_~GH3miFDexF~k7!T%PE3T=)dByWAYfgT+oZ?T}@39YVfpdm!RKCLVNxP0x z^;-DIMk}4SZi3&+oOKJ1eL%p$&E4}lq!D2Olu5y}0yv z?j10koxcQ!J_jR7P&trmJs35-_Qf6Kz8~#aY5I33sJT6x^{Q@yTDxdhl9K=^328~U zK4DPDE@yHIll;D%W3ux+x!3iMYy^)EDD!5&V;RYyoJzBX7MX(b$_sIlPXZORI(xbD zEKrdFyuYlD~SK|jo?tsaqpJ| zX)c57H3@z-c$mRI!yBcI@(nY5kt}NQn&HWs8wT?K;>VtT`9AZqM>xtt#am;em+QYl znK}1|s#(?f1!KYw2 zRG)BtK^shRPa6ljU4of=*}i(t&#>f}ZdBKghfP>!LhImD*h!_Zo%)jihaNvsJ;Ne6 zH7*vlJ>&srTK}ORuAXqUn9BJs=K#+md)W4hb-^e5=;rG8Q}AEZxqU95Gad&`#LfO} ziJ;_~x0b7|@O04vAyH8_JbUf+PH9INqQtv(Ik3*kuJ5|H<77d>vhD$2stJn%8%c;gtr3k)!hH)pwR2N zpLup7ik>G)8Tryty0m_7SGqmQ+gCOlD%B94%zv89Wru2afj73v?)WTe*k_w*jXJH~ z)0VE5XsCF(t#af!zSNZHEl6uYYpX?A2yFv8U3cxcT{nzwkFU;6FU-*^nJHKGj~x1P z*tOdiPNIKojr5451qRuttetvoF~|fC4E65Jsqt?iI*K!bXAjX)28!if)1NTPlM^vl zndqd*z*PE|C7=n*UHg9h5NI(kLU?(%g3jYJ>&{dk=my$P5)xEE&zkt+r%dX>g#WJj z9|6W5I&Y1|YB2O3`Ud9?fDx*}U7vOqjM9C=VV*?a^z&W#A$bPWf|y^;C?oC&g{!W&IS6)z(_8S%W2H7|~Fsy<9U&G^1*Y7}} z?SO^UU=BjGa*K||??w1MrT8^2-w~PEadyy-ff#>&`ZQ}G4lJQ)uCOlFq#tdVQ51eXh&~cy==e*T7}AO3z1mR5B{SX^oR5b zft#wog^GYaQZcqKYZe&dwuvS;yTC9goTp+E3`T^@S*I*}Fe=}UjSdj~GZtoN&>IJ8 z*$|i2zGXz;D0Yp_N&|IrWA>z=5U3k%F>>mpUYI0pdQAHO%JGJ#75jEjer@Fg8jYZy zDK|uj691T_6MLz`Aoyp(O@KEj7$ zW&4%?fT}tljP3(rai~6 zQKSw{hac+sXhbl{T_<5T!TbiA%Q(SLf<`xU2zKg^AF(D#WmvD5`4>Od3+xsg^V|su z7D|g(?=(E$0a@l;8pmSUL?pht;3RXs9%nN{_a$ZtL5XRXgmMArl|s42xUAL_;5N~Z zt5Z4E%1L=pKYRSaQV9jzkdjey>@0*fYoD$et!(AVyU=6pzrCV(x|i)3ZnON1So=*i zS3h4a(*}##-#R*eyoQZi{vf{!8|;|puY=>|odJG({qgYEg6gY`TW}c}OgtC074B_i z9>szjC9`kQ)jiEkooMLRe>rWa*KE-UivwpprtvIirR~! znc0$hlp-j1 zq>0R3M|6&E88J!IN6(U~Y9boY`<03}|7e5$Tgh|17iD3f{8_H_QqpJ4i!W3-m4{*e zw>Gbi-oc1yd_a{XJ4QLp^Z&GJlK%Qh{__6Epy};ZKUzY8HZ)l@$*=_doVa+lG3gWH zwqD&p>LPtw*wxGQJ{X5I6T&m(!7vdrEU(;4bP}V6P4ycX9Y4|wR<(dyX#Bx6L>-iP zszl~Pl1C?#X%AXO2p{hJ>VvKC*A=@!*(t}`Xub#Kwr872G09K=Bjx*LQbC10 zxR+Q@^vuh;{GEMgL8blHFrn#!db44E#t#Qj1xy2|VuhDK|Ac@lueF%9?=PrIJFszn%uhKFlQeHm7vXCxRVqJr0Ba;)jIO z)Sn8jEs$m5$YH zNJ8;7&#yfP>QOH8>`jhb9x6XPvl^9`M@`r`-&MC9)Uo`VUZl`Ivq}e+k1Z`Yp>!WRnpc`-PeHZHvde5!3;Tvy&ail`!Am=tP z$ge3E*a${iHqSna^ug_}XM65PfttT?ae~QVP{Kx>{hx_0IvBVxvho(FGvD4Pm67}T zu1fjvCZbnn-fI17P@Z>>$&G9Q^`!Y>+UgY27wkDwutEk@#@D%vI?D+Ua^&80oFe+@ zx7>bNlCK|mXIC^gll;z{?ZVkk?tk&3OGekp{U7Tn+}H!EWn%XHyM!-of``WsXn|@! z?l7=`=%^0EcMT0B?|Zr5J!U6*Y{1Py>oy(KP%GP!c??jaq7ECvzk;Hjf4A%wsTWgC zg-dtYk-VSHbBUekwP`V}u+|uY7c={kXL#YjA?if67<6r!MNZ~i9 zC-OMU!b{f`GuEd%<2rM$0!^!xf&a)nh7NPx7PpSL)Scg`gS$+-B8(o{JZhP@1}4mX zelQh06*_q?2c~8Ak+ZeBaPRi(7Yi-eVb1iohUGb0$^MaVxG%M_=hn4W+_$}6$i`2@ z{nC)Iz}7#oUDRML7P=U=7nZ6NU0Da)(36rXtHfd3uCSM`=(G69+=z= zw^cm_=bN(=0|VZ{MI}5hsag_lC-+yDvnRmgU}H!^eCpL`0R*x3 zKh3@0@%*eUp0V~tR5{MxLHXi&!h+3oyFGa6>UH`hR|{UT_NB1w8B3%z_s|a`$7sbv zu{$ovA2qle;Uk1%{?H3+)2mSaTgb9S&JZ6jo{7HNZH(H>jcR)jUO_{siM;7seY9wV z@Y`&?hAv*_heA)X>y8?o{pe@E{m(BCRSd*j+-@2A6+`NJVx!ynF&v$%_PS3EBdmBv z9ov*8MFT)%y(g_X_E}&ZFX&2ow`5YLK`%c#RiS1M2LHrob(c$E7$-zeq)CC1GQFua zhxFS`-!E~5^@5sRRj8a+2Ws03k;rn=PaK_Rb#a9B#}`7be0oLnzt(o?<9CQ|zuVCp zHBIu?JpO#Ax&rB=RsQVz77Z$}Zh1{_H>nrdr>Nx@K)uizJMlb`=qL$o=lP@$PvZ-V zks>-j<4;48jwR7kJ&U`Z5Iyzol`RPkRH<*Kw8#-q@68yOFJ}?{99DarOYZs4gMJ}J zgm<-`WzRm4xaxLI#?%tMR@WjP-$LqCy}d?K52-f|!hz|dO9@{hk`>xXU-V~ndG{wz zg7F8R{w8t%l}c9;CjG|WGLdCU?1VoT!mIC-K4_|UxobSpbJIp{N9rB_#g9`Q{N71v zD!9r*?fDJ2is#khHgm2Eqlc-@AGh1#9&`NzEFMk#&2Hm`4b!d;+s4}PML+dn=XOXa z#w7{%>t+@vU~l0Wpw|!#``RNH_k{0<<353FT7?C0dN!%$-zx{FU#Y^^6Z-Jryv~rQ zsWToVu8aS9a5Ek<&#Az}ah^3&jpfkcefH!sD(>$g-GBcd$^|0ZC+Tt z4Q^3~UF^Ku;K7Q=^X|{U_WFMKGUK*^Z>Mf((8N0UvCj8DCKpw)IvN42a{~=hjQ`nk z8o|;x9TW{NA%yk3@D19l@7>&j2-bTZE4Vc<~)WPEe21DMq=dY3x18da~L(s-tyv84QRS|-ae7?0DVo{ zVeyJ;(4R=TzxeYC3=+(2lW8!HazqP#)&gVZ*aML>{|oS{dx(S4!(n#J&=!mdbKx9~ z7ofQMn#3*;eIbyw#?g)Rz0!qydG-l_Qrb53Q1mz`^49q=S3zmc5^Gf_{qxM{Tx&@B zj9R0efkdx7h{=y$Y~0Sea_(VXM3w8i25Tp zZNnZN5 z<0ZwBCW`kXn)?RP#yn5oBXizk@EpW*47?5EZ45jY@mvGX3$SzOdp_d%2HqC&b_U)q zz&;m!Zy$VPTHg5cTUuIL`h(s800960>{#hPm0i>}Wv+zE^rVD}6p7}uOd(Q+L=sV` zM9EYV3MG=1hzuziQc;F9nF*CK8In1}!MSCY>0Re^e|g{MA9z1J_m}Ip_g;JL&F{MQ zy4GPLKT4Bm=1b=^g3`~xVaoUvMd;}cDrLbszM1dgAaMPi2SM9wUq8-|f*8nZs7cXa zZOqg4I#$jgT_Bn()ZheioSgkxMGmm#ylJ}F_*PIhe?DNz`5GwmiX4Z(ZUlAq`9a0h zL!dcow9YlR2KMugHvgjB1_SeTgE6yD!C2dKT|1uyn3V>-_o6p+H-J4Q_`&%H`NR0e3Y{DGof#hx>Q;1U~Yegvbl513|}PAm(9EnsVw& zh<`7XR`YlTBvkx<(0>U74ks6l{ z5Du#y?_bvh!mlOH`rH#BA`+Uo2{|B=4sV#8ybMHUZIJJ?SRnEfRxB;!2jcT1B{7jS zAnGsQ4~!52;wLAyYw0f_y3B3**?55HyB6kFQ3AxE{8~rTH$V`A-}0S#f%tvSA^mGN z5Mx@Q>+dW8F`*!5-&_U6q{XFf&fP#v%@%QKrUUVZ&1S^$9AbXf$G|!uX2RF(F@6Tb z>_w3X_xC`|Z9QWhCXP6#&B^frh>;__Zwj zi7?{m1*zmx7-srkAEm-U!n*01DW#c2Dv6^tO%-dQkTNOmwZhAFG3C$Sd0_HWx>CmN z4&3B;@w?Krt!+<$UzX#^_){?uA?H_VM=sPT~_X&L8oK4EN6@mc1_+brqfZL~=GFS^z;hvc6mG?;p;eoifQS`Djh%Msg zSaw_oo}B)|pUaNAiQb~u>h~F5j)f-$=A40y&?b)`jggR3!EyiU!~w{iQ**0Z_Z^DX z7x}Cz;)YMI$4eo{9lp%LYy>S7YJL0CqeVACV?&{~rNK#PZy&!tIb#pq-`D0Cb-sY! z02%XF3(C;H``XcqoA&_K#rNl!=2{r=JoBdUYz+(=ABl|Zd;~*FUo2nt!vuz&7Sk+s zqk$%V)>zi;3DCT|-Y9HQ1lpHFHtSbc1HrMkx;gqF5IgG&E_kv6VU@@;nwkZK_lQh8 zlNu1=HeKPf@<1f3KB+6$1|qxgrXm;er&w+*Ww#a(l~q5znVW(59{fFYG7fpKZ2Z2> z9eFSG_(TT{h~9un8P6*~47}MTDn0}RE!pban)S&0<+dxgt_5PG+_CbWEb{*Z_lb}r zK#aez2(_L@{AG-^npFbC!p;-NzMu{^ zh7a~DA|BfL=*lI;a9uZ-NJPd3|D}%^^@;MEV!f2vM-xWT#rc%Ex$&RrXTvEAX5)6n zg6Dwg=%0%Nt&+gZ?XhS4H6JkN{9#$}p9L1XOJ!U~`hj)%bg+P2J+KyEGW>@*AJ{HT zuavHT0n5mEez}7R>;ij1KwnaFg^LFWlItthZTHQy^MaK>+=or3bwD`ejrnkp5QzT0 z4@i=68>f_&*SyRIh3?2P)jw0PJu4zHK}`WD&s`tCQ}F|}+_%qetG9s`v*}Kr^b??G zqq$70ZUY#Q@rRG)I~iWe26Hlx6*+IWt$kf=VGKBud7K|K9dEu>3YR*3RrRFLfsc#d z`FI!s{|u97J?#Z>>*zUqk>-PNPqy4Y^OO-paWG{*R-1(QLcfL5rD2dX)A-bG|38rK z&R zbZk<0YE@>3U+S?NlNM8nMfxAF_DbYJ2piFPGbAi9N|96mk=f?6rN;jlgOUDTE* zcM9@s3aS%D8{xjOZRIevy92O6@MgE7n);D3DnNeCCd$$Mi@BfdOm2uKd&OzLprFH@DFc z(8dJyw!j^2NZo`e)xIwS>(_$&^G$q)h}|E!l->WOkMGJKUCx>llpYd?QI`C}BE*1WZ9404dj)X0EIIQ*#szrD z_$7j$w;O(}0Y1jQd{|~y%rcf74*ZP$XF(uU;&k4~UJzvD5djfSX*S2%b`anF$3Dk8 z2_&}np7!b929o6ZI>tFAuz_)oJIIpfY$eyX>$Yt>>>mre$oV}!8L@E#Cqa#oM;Wxq z_5F|UY!^r?0)z7U?46~yVB+w;x?#oRtN?-bvVG z=q-i5d#-Vv+4Vpz$=l`nq67vHZ{@AKtN}yUk{mS?sX(LbuTPBL2ei_Fpup%KK&%g! zt@K6TW2cqk9f3aSZl~Dsm~0?k^#+e*aA6$?xv|s({Y`Cqm8>u7=10jp_IUJ7oiq!< zMATJZ`B%z?Heh~v`_ ztCt|}nFap@Y9R8;AJOnYT$Qok;~S!aGQ7e(^XhN6^xi~#UC|?#gjgA3xPbbYryAu| zWBpvXL)-th6S2%)+Whyw^ii~^A@vxu8>NoKmVlzGJl7;BU5v3KrT3Fw_*T^~ltIS) z2FmbLv8>-#>XcE&`UT3k?cGOZ@$r;NGVafM&Y-xCKa^QAe!kz!&0PY0+24J@T<;;D zq;LdS|LzAiWhH6;nj5f)ktYt=!-qw^oU(y~k?$*TPNwovm;D0nzj=Yj)Rx`&wE*z` z%@0!KJ_SZzGuY0E?*@v&!M9>N^+A=K*Cyxnd55OLFKvTE*Y^A6?+pf1a{f3ueoZUUOsL z-tL+#h`vG6(e!>ICA@yrjowd)N4nsB<~x(k$Ab|U6jr;@eUj+guE_+%Ex{Rg1rhgI z`EK+@Jg#I|I!dn>ToP7z{z7oe@0nV}jJ|loTZm0#`YD?K(noFR#=|nNQh4m1^B zJ*V}8GQt>trc5fi+qb;Wr_3w-zK1lzCZ6>cTq4{$h|3lRjct}Lsw6Vl;MXQIRolPHMvzT&K+stE};bK_r`7r{$KXZBqWYayp_ zfOE9128!7P=_^pR{g(gRrToxPQ{=gFl^S#qeGad)U7+WdZpKN!8tAWCZokr12L@MJ zOlnFs!;s5*li_f0pc!ovzbfz$Xgy=W9Y65>)L&J{@(|x`A6Y-3Qh>(ZaApr4x4`oL^1f#`U~d1)`= zeWC6m^n-KnKOJ^=LhN)clih~8;ga0ZI*YiuHuWz0vH2&hOfRucEv)P}E=OOvV57N1 zegN_H6R`(w|I)|s#RFBZ#G8PLM9xgNZ8y{20v}^s2ZG_&zA<*?AjFu@1F_nh1Fn-< zu$r;{9VBj;rX;D?fMlb3Zp%h(Shr!*w`GFGu)#yPfA1UxWIwabd;EF_n?#2N(ki!t zys|^sg`+!Pu_?yU6`&jJOKWBG`<}86zjB{AQ zY4o3NGp-D93maQf)SLt!jPtI*wLQDEm=X)&Hks#sXt0uJ<_J7G>bbWhr462~ZsFM7 z`WUjd_IYeP9|r~Dmo{*f9EUG024rNsfQ&?}YY?6mhZP{UGw zTW#fr!SsW@FVR=fxI^Xd6xaZ*gjMW)H~Jg-tm+j^=;vKIdpRPJHxZ>T_GVxooPQv1 zdkyxxjKt`Z+U`jUd|iO6?-h0nS8E{}*@@hQZ2cT{bZ=NS62G5b^ZC-Ge# ztK8{WAA@{X;MK~&dN6*kb#*8ev1)yNH^xnHI|nHXA@0pp3BbO2;wtBOGxqNjAB<8r zV*Q!yE>|kWcX>+d;gSUGzo(v^?e)ceXL>Ou-wNaYC}rx*yhq+!Y^V=M|1*0%?zsFD zuPJSDA+>o+iFrz)lokHO62}+ zOO!LWaI(R6MjlJpQ5@kV{N)|&V&q#6yFI3~nM2RR9x^W|99h2c0(+7Z;Zq zKp)d|U|EPR4D7nEzE(LKhGdFU#BUD)4Wi;={VRahtKMqy8F{xq0h(A0fN&2G!J~nQ z^)F9m$9FZ0`o+~1@on{F#Orw=z8qbA`#G;=PjrE~<+;-iR>KGU&7HjNEzf-5)VK&A9>J*)ct@NNpd)(I2bk zus;}&+LrKo2kJv5OvY|GV$X28DAujXU_GU3tOrwyDewI;|5P)(>dqUe6CdkR;sokM zI&Xi13gTDM+9UMeshlsGF2=uWR`}b`LcGVU>9&v_tb4QJQEJ@SSImB!9ovNbpA)c5 zG(^4485Ekjjv(HV;*P|=X|D3dhmQ~FeF)XD;|#rD*ITm)>*M^*w7L!SKD~_Je?Qiv z1$NhM=jnZiMniuX*1ZMaoey)4A{NBxwaEQneSpA3j=o_^1gs@dZbb0n<#{7eWQ;># z`x_Os3x#4hT^kQxOd;pwfoaYYb$-GW>eXqZ@@q!a0|01|>;{@O6syT3G zoHGv2jQb3M8=2p|rB9!w@-l$UAIR7j1ff9<^AD1^;V$@HvJYW~2Y13M>a&yKk)${0NW&_4TF}`RxS|R& z$UMc9`Ss>C%~03IYV53a2D-&c<5E9K19iW2)RGhWF!*K9(wL7gfVM;Nw1_138{f|y zx|>!Fgz#sU_|mUH9Ft8?|5*-%bG!=A&N3iEG8K1kMt`5Up)mXw`k7ZJ1J&R8q0fpl z8184ox^gVcdkV+Jq76$q(Pw>X7LHwsb?2-56_=M-C+Y+{bm>07^s>>A>&;=w*x(F#1YIBRg9Dar k`oH@4@B8oj@B8oj@B9Dt{SN>D0RR630Jqw4#{ew=0IF=)IsgCw From b5f65b058a4fe9dd299115633fc5fa38e6090823 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 07:23:52 -0400 Subject: [PATCH 04/14] Rename mean_flicker_time to flicker_mean_time --- chandra_aca/aca_image.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 10ec9eb..e2025d2 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -338,7 +338,7 @@ def _read_flicker_cdfs(cls): cdf_bins.append(hdr[f'cdf_bin{ii}']) cls.flicker_cdf_bins = np.array(cdf_bins) - def flicker_init(self, mean_flicker_time=10000, flicker_scale=1.0): + def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): """Initialize instance variables to allow for flickering pixel updates. The ``flicker_scale`` can be interpreted as follows: if the pixel @@ -347,14 +347,14 @@ def flicker_init(self, mean_flicker_time=10000, flicker_scale=1.0): that increase the amplitude. For flickers that make the value smaller, then it would be 1 / (1 + x) => 1 / (1 + x * flicker_scale). - :param mean_flicker_time: mean flickering time (sec, default=10000) + :param flicker_mean_time: mean flickering time (sec, default=10000) :param flicker_scale: multiplicative factor beyond model default for flickering amplitude (default=1.0) """ if not hasattr(self, 'flicker_cdf_bins'): self._read_flicker_cdfs() - self.mean_flicker_time = mean_flicker_time + self.flicker_mean_time = flicker_mean_time self.flicker_scale = flicker_scale # Make a flattened view of the image for easier update processing. @@ -375,7 +375,7 @@ def flicker_init(self, mean_flicker_time=10000, flicker_scale=1.0): # time assume the flickering is randomly phased within that interval. phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) - t_flicker = -np.log(1.0 - rand_unifs) * mean_flicker_time + t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time self.flicker_times = t_flicker * phase # Pixels where self.flicker_cdf_idxs == 0 have val < 50 (no CDF) and are @@ -426,7 +426,7 @@ def flicker_update(self, dt): self.flicker_vals[idx] = val # Get the new time before next flicker - t_flicker = -np.log(1.0 - rand_time) * self.mean_flicker_time + t_flicker = -np.log(1.0 - rand_time) * self.flicker_mean_time self.flicker_times[idx] = t_flicker From eeb8e31d276be303f919e075fb4f70bc35a08f10 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 09:34:57 -0400 Subject: [PATCH 05/14] Add numba implementation for about 6 times speedup --- chandra_aca/aca_image.py | 87 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index e2025d2..ae7a18b 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -8,6 +8,7 @@ import six from six.moves import zip +import numba import numpy as np from astropy.utils.compat.misc import override__dir__ @@ -326,7 +327,7 @@ def _read_flicker_cdfs(cls): # Get the main data, which is an n_cdf * n_cdf_x array. Each row corresponds # to the CDF for a particular bin range in e-/sec, e.g. 200 to 300 e-/sec. # CDF will go from 0.0 to 1.0 - cls.flicker_cdfs = hdu.data + cls.flicker_cdfs = hdu.data.astype(np.float64) # CDF_x is the x-value of the distribution, namely the log-amplitude change # in pixel value due to a flicker event. @@ -429,6 +430,90 @@ def flicker_update(self, dt): t_flicker = -np.log(1.0 - rand_time) * self.flicker_mean_time self.flicker_times[idx] = t_flicker + def flicker_update_numba(self, dt): + _flicker_update_numba(dt, len(self.flicker_vals), + self.flicker_vals0, + self.flicker_vals, + self.flicker_mask, + self.flicker_times, + self.flicker_cdf_idxs, + self.flicker_cdf_x, + self.flicker_cdfs, + self.flicker_scale, + self.flicker_mean_time) + + +@numba.jit(nopython=True) +def _flicker_update_numba(dt, nvals, + flicker_vals0, + flicker_vals, + flicker_mask, + flicker_times, + flicker_cdf_idxs, + flicker_cdf_x, + flicker_cdfs, + flicker_scale, + flicker_mean_time): + """ + Propagate the image forward by ``dt`` seconds and update any pixels + that have flickered during that interval. + """ + for ii in range(nvals): # nvals + if not flicker_mask[ii]: + continue + + flicker_times[ii] -= dt + + if flicker_times[ii] > 0: + continue + + # Random uniform used for (1) distribution of flickering amplitude + # via the CDFs and (2) distribution of time to next flicker. + rand_ampl = np.random.uniform(0.0, 1.0) + rand_time = np.random.uniform(0.0, 1.0) + + # Determine the new value after flickering and set in array view. + # First get the right CDF from the list of CDFs based on the pixel value. + cdf_idx = flicker_cdf_idxs[ii] + y = np_interp(yin=flicker_cdf_x, + xin=flicker_cdfs[cdf_idx], + xout=rand_ampl) + + if flicker_scale != 1.0: + # Express the multiplicative change as (1 + x) and change + # it to be (1 + x * scale). This makes sense for positive y, + # so use abs(y) and then flip the sign back at the end. For + # negative y this is the same as doing this trick expressing the + # change as 1 / (1 + x). + dy = (10 ** np.abs(y) - 1.0) * flicker_scale + 1.0 + y = np.log10(dy) * np.sign(y) + + val = 10 ** (np.log10(flicker_vals0[ii]) + y) + flicker_vals[ii] = val + + # Get the new time before next flicker + t_flicker = -np.log(1.0 - rand_time) * flicker_mean_time + flicker_times[ii] = t_flicker + + +@numba.jit(nopython=True) +def np_interp(yin, xin, xout): + idx = np.searchsorted(xin, xout) + + if idx == 0: + return yin[0] + + if idx == len(xin): + return yin[-1] + + x0 = xin[idx - 1] + x1 = xin[idx] + y0 = yin[idx - 1] + y1 = yin[idx] + yout = (xout - x0) / (x1 - x0) * (y1 - y0) + y0 + + return yout + def _prep_6x6(img, bgd=None): """ From a80da64dc22d0071da7af9c0f86b0a0d85804a6e Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 13:43:28 -0400 Subject: [PATCH 06/14] Make a single flicker_update() interface and default to numba --- chandra_aca/aca_image.py | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index ae7a18b..912dbfa 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -383,19 +383,37 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): # modeled as not flickering. Make a mask to indicate which ones flicker. self.flicker_mask = self.flicker_cdf_idxs >= 0 - def flicker_update(self, dt): + def flicker_update(self, dt, use_numba=True): """ Propagate the image forward by ``dt`` seconds and update any pixels that have flickered during that interval. - TO DO: use numba for this once numba with np.interp is available in Ska3. - (E.g. 0.43 has it). This will probably be substantially faster. + This has the option to use one of two implementations. The default is + to use the numba-based version which is about 6 times faster. The + vectorized version is left in for reference. :param dt: time (secs) to propagate image + :param use_numba: use the numba version of updating (default=True) """ if not hasattr(self, 'flicker_times'): self.flicker_init() + if use_numba: + _flicker_update_numba(dt, len(self.flicker_vals), + self.flicker_vals0, + self.flicker_vals, + self.flicker_mask, + self.flicker_times, + self.flicker_cdf_idxs, + self.flicker_cdf_x, + self.flicker_cdfs, + self.flicker_scale, + self.flicker_mean_time) + else: + self._flicker_update_vectorized(dt) + + def _flicker_update_vectorized(self, dt): + self.flicker_times[self.flicker_mask] -= dt # When flicker_times < 0 that means a flicker occurs @@ -430,18 +448,6 @@ def flicker_update(self, dt): t_flicker = -np.log(1.0 - rand_time) * self.flicker_mean_time self.flicker_times[idx] = t_flicker - def flicker_update_numba(self, dt): - _flicker_update_numba(dt, len(self.flicker_vals), - self.flicker_vals0, - self.flicker_vals, - self.flicker_mask, - self.flicker_times, - self.flicker_cdf_idxs, - self.flicker_cdf_x, - self.flicker_cdfs, - self.flicker_scale, - self.flicker_mean_time) - @numba.jit(nopython=True) def _flicker_update_numba(dt, nvals, From d9395070f9ef2fd39910199bb658c956ce4a410b Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 13:54:02 -0400 Subject: [PATCH 07/14] Add comments about model review and provenance --- chandra_aca/aca_image.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 912dbfa..0eeb247 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -316,6 +316,7 @@ def _read_flicker_cdfs(cls): """Read flickering pixel model cumulative distribution functions and associated metadata. Set up class variables accordingly. + """ from astropy.io import fits @@ -348,6 +349,14 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): that increase the amplitude. For flickers that make the value smaller, then it would be 1 / (1 + x) => 1 / (1 + x * flicker_scale). + The flicker_cdf file here was created using: + /proj/sot/ska/www/ASPECT/ipynb/chandra_aca/flickering-pixel-model.ipynb + + Examples and performance details at: + /proj/sot/ska/www/ASPECT/ipynb/chandra_aca/flickering-implementation.ipynb + + The model was reviewed and approved at SS&AWG on 2019-05-22. + :param flicker_mean_time: mean flickering time (sec, default=10000) :param flicker_scale: multiplicative factor beyond model default for flickering amplitude (default=1.0) From adbfd79f29c7318fcb1147dacb63d00d82957dee Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 13:55:30 -0400 Subject: [PATCH 08/14] Update version to 4.26 --- chandra_aca/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chandra_aca/__init__.py b/chandra_aca/__init__.py index 8490495..cff429d 100644 --- a/chandra_aca/__init__.py +++ b/chandra_aca/__init__.py @@ -1,7 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from .transform import * -__version__ = '4.25.2' +__version__ = '4.26' def test(*args, **kwargs): From c78571058837bce175597a63f6bfe8c13048d49d Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 14:29:54 -0400 Subject: [PATCH 09/14] Add tests and ability to set random seed --- chandra_aca/aca_image.py | 12 +++++++++- chandra_aca/tests/test_aca_image.py | 37 +++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 0eeb247..f1b91e7 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -340,7 +340,7 @@ def _read_flicker_cdfs(cls): cdf_bins.append(hdr[f'cdf_bin{ii}']) cls.flicker_cdf_bins = np.array(cdf_bins) - def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): + def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None): """Initialize instance variables to allow for flickering pixel updates. The ``flicker_scale`` can be interpreted as follows: if the pixel @@ -360,6 +360,7 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): :param flicker_mean_time: mean flickering time (sec, default=10000) :param flicker_scale: multiplicative factor beyond model default for flickering amplitude (default=1.0) + :param seed: random seed for reproducibility (default=None => no seed) """ if not hasattr(self, 'flicker_cdf_bins'): self._read_flicker_cdfs() @@ -367,6 +368,10 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0): self.flicker_mean_time = flicker_mean_time self.flicker_scale = flicker_scale + if seed is not None: + np.random.seed(seed) + _numba_random_seed(seed) + # Make a flattened view of the image for easier update processing. # Also store the initial pixel values, since flicker updates are # relative to the initial value, not the current value (which would @@ -458,6 +463,11 @@ def _flicker_update_vectorized(self, dt): self.flicker_times[idx] = t_flicker +@numba.jit(nopython=True) +def _numba_random_seed(seed): + np.random.seed(seed) + + @numba.jit(nopython=True) def _flicker_update_numba(dt, nvals, flicker_vals0, diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 0b43c92..7c7f31e 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -325,3 +325,40 @@ def test_aca_image_operators(): [6, -93, -192, 9], [10, -289, -388, 13], [14, 15, 16, 17]]) + + +def test_flicker_numba(): + a = ACAImage(np.linspace(0, 800, 9).reshape(3, 3)) + a.flicker_init(flicker_mean_time=1000, flicker_scale=1.5, seed=10) + for ii in range(10): + a.flicker_update(100.0, use_numba=True) + + assert np.all(np.round(a) == [[0, 81, 200], + [326, 176, 609], + [659, 720, 1043]]) + + +def test_flicker_vectorized(): + a = ACAImage(np.linspace(0, 800, 9).reshape(3, 3)) + a.flicker_init(flicker_mean_time=1000, flicker_scale=1.5, seed=10) + for ii in range(10): + a.flicker_update(100.0, use_numba=False) + + assert np.all(np.round(a) == [[0, 111, 200], + [219, 436, 531], + [470, 829, 822]]) + + +def test_flicker_no_seed(): + """Make sure results vary when seed is not supplied""" + a = ACAImage(np.linspace(0, 800, 9).reshape(3, 3)) + a.flicker_init(flicker_mean_time=300) + for ii in range(10): + a.flicker_update(100.0) + + b = ACAImage(np.linspace(0, 800, 9).reshape(3, 3)) + b.flicker_init(flicker_mean_time=300) + for ii in range(10): + b.flicker_update(100.0) + + assert np.any(a != b) From 7b15422618425231a25aa4b8c5e41641e4c76139 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 14:32:33 -0400 Subject: [PATCH 10/14] Add file creation comment (and fix extra whitespace) --- chandra_aca/aca_image.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index f1b91e7..530d672 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -316,6 +316,8 @@ def _read_flicker_cdfs(cls): """Read flickering pixel model cumulative distribution functions and associated metadata. Set up class variables accordingly. + The flicker_cdf file here was created using: + /proj/sot/ska/www/ASPECT/ipynb/chandra_aca/flickering-pixel-model.ipynb """ from astropy.io import fits From 48f8b0f18f3ab0a05faac82311047ae57f7f7fb4 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 22 May 2019 14:40:17 -0400 Subject: [PATCH 11/14] Revert unrelated change in PSF library (probably leftover) --- chandra_aca/aca_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 530d672..04be562 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -641,7 +641,7 @@ def __init__(self, filename=None): if filename is None: filename = os.path.join(os.path.dirname(__file__), 'data', 'aca_psf_lib.dat') dat = Table.read(filename, format='ascii.basic', guess=False) - self.dat = dat.as_array() + self.dat = dat # Sub-pixel grid spacing in pixels. This assumes the sub-pixels are # all the same size and square, which is indeed the case. From ad4362708b3c300364a234f7f286248d95319dc9 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Sat, 25 May 2019 17:25:04 -0400 Subject: [PATCH 12/14] Make the flicker mask be ACAImage and be independent of cdf_idxs --- chandra_aca/aca_image.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 04be562..ce48130 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -382,6 +382,13 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None): self.flicker_vals0 = self.flicker_vals.copy() self.flicker_n_vals = len(self.flicker_vals) + # Make a bool ACAImage like self to allow convenient mask/unmask of + # pixels to flicker. This is used in annie. Also make the corresponding + # 1-d ravelled version. + self.flicker_mask = ACAImage(np.ones(self.shape, dtype=np.bool), + row0=self.row0, col0=self.col0) + self.flicker_mask_vals = self.flicker_mask.view(np.ndarray).ravel() + # Get the index to the CDFs which is appropriate for each pixel # based on its initial value. self.flicker_cdf_idxs = np.searchsorted(self.flicker_cdf_bins, @@ -395,10 +402,6 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None): t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time self.flicker_times = t_flicker * phase - # Pixels where self.flicker_cdf_idxs == 0 have val < 50 (no CDF) and are - # modeled as not flickering. Make a mask to indicate which ones flicker. - self.flicker_mask = self.flicker_cdf_idxs >= 0 - def flicker_update(self, dt, use_numba=True): """ Propagate the image forward by ``dt`` seconds and update any pixels @@ -418,7 +421,7 @@ def flicker_update(self, dt, use_numba=True): _flicker_update_numba(dt, len(self.flicker_vals), self.flicker_vals0, self.flicker_vals, - self.flicker_mask, + self.flicker_mask_vals, self.flicker_times, self.flicker_cdf_idxs, self.flicker_cdf_x, @@ -430,10 +433,11 @@ def flicker_update(self, dt, use_numba=True): def _flicker_update_vectorized(self, dt): - self.flicker_times[self.flicker_mask] -= dt + self.flicker_times[self.flicker_mask_vals] -= dt # When flicker_times < 0 that means a flicker occurs - idxs = np.where(self.flicker_times < 0)[0] + ok = (self.flicker_times < 0) & (self.flicker_cdf_idxs >= 0) + idxs = np.where(ok)[0] # Random uniform used for (1) distribution of flickering amplitude # via the CDFs and (2) distribution of time to next flicker. @@ -474,7 +478,7 @@ def _numba_random_seed(seed): def _flicker_update_numba(dt, nvals, flicker_vals0, flicker_vals, - flicker_mask, + flicker_mask_vals, flicker_times, flicker_cdf_idxs, flicker_cdf_x, @@ -486,7 +490,10 @@ def _flicker_update_numba(dt, nvals, that have flickered during that interval. """ for ii in range(nvals): # nvals - if not flicker_mask[ii]: + # cdf_idx is -1 for pixels with dark current in range that does not flicker. + # Don't flicker those ones or pixels that are masked out. + cdf_idx = flicker_cdf_idxs[ii] + if cdf_idx < 0 or not flicker_mask_vals[ii]: continue flicker_times[ii] -= dt @@ -501,7 +508,6 @@ def _flicker_update_numba(dt, nvals, # Determine the new value after flickering and set in array view. # First get the right CDF from the list of CDFs based on the pixel value. - cdf_idx = flicker_cdf_idxs[ii] y = np_interp(yin=flicker_cdf_x, xin=flicker_cdfs[cdf_idx], xout=rand_ampl) From 09e7cf3c502afa2477bd37643abd03ba777ee076 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 29 May 2019 16:08:02 -0400 Subject: [PATCH 13/14] Use simpler expression suggested by Larry --- chandra_aca/aca_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index ce48130..e094ce5 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -461,7 +461,7 @@ def _flicker_update_vectorized(self, dt): dy = (10 ** np.abs(y) - 1.0) * self.flicker_scale + 1.0 y = np.log10(dy) * np.sign(y) - val = 10 ** (np.log10(self.flicker_vals0[idx]) + y) + val = self.flicker_vals0[idx] * 10 ** y self.flicker_vals[idx] = val # Get the new time before next flicker From 00be9b70e207ad577d0247620284c310eb5acc20 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 29 May 2019 17:19:41 -0400 Subject: [PATCH 14/14] Add test mode for comparison to C implementation --- chandra_aca/aca_image.py | 44 ++++++++++++++++++++--------- chandra_aca/tests/test_aca_image.py | 36 +++++++++++++++++++++++ 2 files changed, 66 insertions(+), 14 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index e094ce5..462060b 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -369,8 +369,9 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None): self.flicker_mean_time = flicker_mean_time self.flicker_scale = flicker_scale + self.test_idx = 1 if seed == -1 else 0 - if seed is not None: + if seed is not None and seed != -1: np.random.seed(seed) _numba_random_seed(seed) @@ -395,16 +396,21 @@ def flicker_init(self, flicker_mean_time=10000, flicker_scale=1.0, seed=None): self.flicker_vals0) - 1 # Create an array of time (secs) until next flicker for each pixel - # This is drawing from an exponential distribution. For the initial - # time assume the flickering is randomly phased within that interval. - phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) - rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) - t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time + if seed == -1: + # Special case of testing, use a constant flicker_mean_time initially + t_flicker = np.ones(shape=(self.flicker_n_vals,)) * flicker_mean_time + phase = 1.0 + else: + # This is drawing from an exponential distribution. For the initial + # time assume the flickering is randomly phased within that interval. + phase = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) + rand_unifs = np.random.uniform(0.0, 1.0, size=self.flicker_n_vals) + t_flicker = -np.log(1.0 - rand_unifs) * flicker_mean_time + self.flicker_times = t_flicker * phase def flicker_update(self, dt, use_numba=True): - """ - Propagate the image forward by ``dt`` seconds and update any pixels + """Propagate the image forward by ``dt`` seconds and update any pixels that have flickered during that interval. This has the option to use one of two implementations. The default is @@ -418,7 +424,9 @@ def flicker_update(self, dt, use_numba=True): self.flicker_init() if use_numba: - _flicker_update_numba(dt, len(self.flicker_vals), + _flicker_update_numba(dt, + len(self.flicker_vals), + self.test_idx, self.flicker_vals0, self.flicker_vals, self.flicker_mask_vals, @@ -428,6 +436,8 @@ def flicker_update(self, dt, use_numba=True): self.flicker_cdfs, self.flicker_scale, self.flicker_mean_time) + if self.test_idx > 0: + self.test_idx += 1 else: self._flicker_update_vectorized(dt) @@ -475,7 +485,7 @@ def _numba_random_seed(seed): @numba.jit(nopython=True) -def _flicker_update_numba(dt, nvals, +def _flicker_update_numba(dt, nvals, test_idx, flicker_vals0, flicker_vals, flicker_mask_vals, @@ -501,10 +511,16 @@ def _flicker_update_numba(dt, nvals, if flicker_times[ii] > 0: continue - # Random uniform used for (1) distribution of flickering amplitude - # via the CDFs and (2) distribution of time to next flicker. - rand_ampl = np.random.uniform(0.0, 1.0) - rand_time = np.random.uniform(0.0, 1.0) + if test_idx > 0: + # Deterministic and reproducible but bouncy sequence that is reproducible in C + # (which has a different random number generator). + rand_ampl = np.abs(np.sin(float(ii + test_idx))) + rand_time = np.abs(np.cos(float(ii + test_idx))) + else: + # Random uniform used for (1) distribution of flickering amplitude + # via the CDFs and (2) distribution of time to next flicker. + rand_ampl = np.random.uniform(0.0, 1.0) + rand_time = np.random.uniform(0.0, 1.0) # Determine the new value after flickering and set in array view. # First get the right CDF from the list of CDFs based on the pixel value. diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 7c7f31e..c55ebd1 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -362,3 +362,39 @@ def test_flicker_no_seed(): b.flicker_update(100.0) assert np.any(a != b) + + +def test_flicker_test_sequence(): + """Test the deterministic testing sequence that allows for + cross-implementation comparison. This uses the seed=-1 pathway. + """ + a = ACAImage(np.linspace(0, 800, 9).reshape(3, 3)) + a.flicker_init(seed=-1, flicker_mean_time=15) + dt = 10 + assert np.all(np.round(a) == [[0, 100, 200], + [300, 400, 500], + [600, 700, 800]]) + + a.flicker_update(dt) + assert np.all(np.round(a) == [[0, 100, 200], + [300, 400, 500], + [600, 700, 800]]) + + a.flicker_update(dt) + assert np.all(np.round(a) == [[0, 80, 229], + [447, 374, 531], + [952, 693, 815]]) + a.flicker_update(dt) + assert np.all(np.round(a) == [[0, 80, 229], + [278, 374, 531], + [592, 693, 815]]) + + a.flicker_update(dt) + assert np.all(np.round(a) == [[0, 80, 185], + [278, 374, 531], + [592, 693, 815]]) + + assert np.allclose(a.flicker_times, + np.array([15., 49.06630191, 48.34713118, 38.34713118, + 28.34713118, 1.03039723, 26.30875397, 16.30875397, + 7.40192939]))