wf3cte possible inconsistent logic in simulating column readout

[jamienoss]

https://github.com/spacetelescope/hstcal/blob/master/pkg/wfc3/calwf3/wf3cte/wf3cte.c#L1470-L1502

                if ( (ttrap < cte->cte_len) || ( pix_1 >= cte->qlevq_data[w] - 1. ) ){
                    if (pixo[j] >= 0 ){
                        pix_1 = pixo[j] + fcarry; /*shuffle charge in*/
                        fcarry = pix_1 - floor(pix_1); /*carry the charge remainder*/
                        pix_1 = floor(pix_1); /*reset pixel*/
                    }

                    /*HAPPENS AFTER FIRST PASS*/
                    /*SHUFFLE CHARGE IN*/
                    if ( j> 0  ) {
                        if (pixf[j] < pixf[j-1])
                            ftrap *= (pixf[j] /  pixf[j-1]);
                    }

                    /*RELEASE THE CHARGE*/
                    padd_2=0.0;
                    if (ttrap <cte->cte_len){
                        ttrap += 1;
                        padd_2 = Pix(rprof->data,w,ttrap-1) *ftrap;
                    }

                    padd_3 = 0.0;
                    prem_3 = 0.0;
                    if ( pix_1 >= cte->qlevq_data[w]){
                        prem_3 =  cte->dpdew_data[w] / cte->n_par * pixf[j];  /*dpdew is 1 in file */
                        if (ttrap < cte->cte_len)
                            padd_3 = Pix(cprof->data,w,ttrap-1)*ftrap;
                        ttrap=0;
                        ftrap=prem_3;
                    }

                    pixo[j] += padd_2 + padd_3 - prem_3;
                } /*replaces trap continue*/

It appears that the outer conditional exists for efficiency and that the subsequent (inner) conditionals are there to differentiate which of the conditionals in the OR op were true. Since both ttrap and pix_1 are updated yet are subsequently conditionally checked this leads to inconsistent logic.

In addition, should https://github.com/spacetelescope/hstcal/blob/master/pkg/wfc3/calwf3/wf3cte/wf3cte.c#L1493

read

if ( pix_1 >= cte->qlevq_data[w] - 1.)

rather than

if ( pix_1 >= cte->qlevq_data[w] )

as it currently does?

Should the two outer conditions be assigned to temps and only then used? https://github.com/spacetelescope/hstcal/compare/master...jamienoss:issue48-wf3cte-inconsistent-logic-in-sim_colreadout_l

Furthermore, the use of ttrap as a conditional switch looks odd. It is initially assigned the value of cte->cte_len which is exactly what is is compared less than to, however, it (ttrap) can only ever be >= cte->cte_len or 0. The intent of its use would be clearer if the condition read if (ttrap == 0). For example, it makes me wonder if it is being used correctly or if there is room for a bug to exist.

Yet another, ttrap is also used in computing an index reference to cprof->data and rprof->data but is also updated, and thus it is unclear whether this variation was intended.

[jamienoss]

@sosey FYI another possible issue (3 actually) - I'll email Jay Anderson and cc you as per the last one.

[sosey]

@jamienoss good, I may not have time to look at this, this week, so move forward with jay. He'll ultimately have to approve logic changes on this algorithm anyways.

[JayVB]

Ok...

issue#1: "the ttrap conditional switch looks odd"... the (ttrap < cte->cte_len) conditional just means that the trap is still active. This entire loop applies for a particular trap, and we go up through the traps from the highest (electron-level-wise) at w=wmax (or cte->cte_traps-1 for you) to the lowest trap (at w=1). For each trap level, "ttrap" records how long it's been since the trap was filled. cte->cte_len tells us the maximum trail length. All the flux in a trap is pushed out by the time we get to ttrap = cte->cte_len shifts. So the less-than conditional means that the trap has not been fully expunged. It's not true that (as you say) "ttrap can only be >= cte->cte_len or zero". It gets set to zero if it happens that the pixel is above the level of the trap. After ttrap gets set to zero, if the next pixel up the column is below the level of the trap, then ttrap gets incremented by one and electrons get pushed out of the trap and into the pixel.

I don't think you can set the conditionals at the top of the loop. pix_1 is not always the original pixel value. The reason for this has to do with fractional electrons and the fact that the algorithm is "thresholded" at certain charge levels. Since all the comparisons as to whether a pixel will fill a trap are done in terms of integers, and pixels can have fractional electrons (at least the way the routine adds and accounts for them), then the fractional charge above the integer pixel value never gets acted on. This caused some strange thresholding behavior when I was developing it. To avoid this, I keep track of the fractional charge that has not been acted on. I add the remaining fractional charge to the target pixel to determine a new integer pixel value (pix_1), and keep the remainder in "fcarry".

I hope this answers your questions.

[jamienoss]

@JayVB

Furthermore, the use of ttrap as a conditional switch looks odd. It is initially assigned the value of > cte->cte_len which is exactly what is is compared less than to, however, it (ttrap) can only ever be > >= cte->cte_len or 0. The intent of its use would be clearer if the condition read if (ttrap == 0). For > example, it makes me wonder if it is being used correctly or if there is room for a bug to exist.

You are correct that this is not the case, as ttrap is reset outside of the iteration over rows, not within it.

[jamienoss]

@JayVB What would be super useful is if you could rewrite this section of logic making it more explicitly clear of what is intended. Many thanks.

[jamienoss]

Below is @JayVB latest suggestion (fortran)

do j = J1, J2
            if (j.gt.J1) then                          ! if we have an inversion of the the
                if (pixf(j).lt.pixf(j-1))              !    density (ie, a readout-cosmic issue),
     .              ftrap = pixf(j)/pixf(j-1)*ftrap    !    then we don't want to flush too much
                endif                                
                                                       ! set up accounting of pixel value chgs
            prel_1 = 0.                                !    charge to be released                             
            prel_2 = 0.                                !    charge to be flushed out                               
            pgrb_3 = 0.                                !    charge to be trapped

            if (pixo(j).ge.q_w(w)) then                ! filled/refilled trap#W here?
                  if (ttrap.le.(_TDIM_)) then          !    need to flush before filling?    
                      ttrap = ttrap + 1                !       increment time since filled
                      prel_1 = rprof_wt(w,ttrap)*ftrap !       flush out amount for this shift
                      prel_2 = cprof_wt(w,ttrap)*ftrap !       flush out the rest
                  endif 
                  ftrap  = dpde_w(w)/NITs*pixf(j)      !    amount to hold in the trap
                  pgrb_3 = ftrap                       !    subtract that amount held from the pixel
                  ttrap = 0                            !    reset the time-since-filled counter
            else                                    ! trap#W not filled this time
                  if (ttrap.le.(_TDIM_)) then          !    does trap contain charge?
                      ttrap = ttrap + 1                !       release the appropriate number
                      prel_1 = rprof_wt(w,ttrap)*ftrap !       of electrons
                  endif
            endif

            pixo(j) = pixo(j)                        ! make adjustment to the output pixel
     .              + prel_1                         !    add the trail emission 
     .              + prel_2                         !    flush the trap
     .              - pgrb_3        

So, there are some id changes to keep us on our toes. The outer condition has been removed, the fcarry code has also been removed, and the residual conditions combined. Regarding the latter point, if this is the desired logic, then it is in fact, an algorithm change and my original concern with ttrap being updated and later used in a condition seems correct.

For reference, here are three fortran versions of this function, the original

      subroutine sim_colreadout_l_uvis_w(pixi,pixo,pixf,J1,J2,JDIM,
     .                                   q_w,dpde_w,NITs,
     .                                   tprof_w,rprof_wt,cprof_wt,Ws)
      implicit none

      integer   JDIM
      real      pixi(JDIM)        ! input column array
      real      pixo(JDIM)        ! outout column array
      real      pixf(JDIM)        ! outout column array
      integer   J1, J2

      integer    q_w(_WDIM_)      ! the 5 readout-parameter arrays
      real    dpde_w(_WDIM_)
      integer NITs
      integer tprof_w(_WDIM_)
      real   rprof_wt(_WDIM_,100)
      real   cprof_wt(_WDIM_,100)
      integer Ws

      integer j, t

      real  pix0

      real*8  ftrap
      integer ttrap

      real fpix ! fraction of this pixel involved
      real fopn ! fraction of this trap that is open
      real ffil ! fraction of this trap that gets filled
      real dpix ! the amount of charge that gets transferred

      real*8 padd
      integer w

      real*8 pix_1
      real*8 pix_2
      real*8 pix_3

      real*8 padd_2
      real*8 padd_3
      real*8 prem_3, prem
      real*8 prem_4
      real*8 pmax, prev

      real    fcarry
      real    fcarry0

      fcarry0 = 0.

      if (Ws.gt._WDIM_) stop 'Ws.gt._WDIM_'   ! more traps than we are allowed

c----------------------------------------
c
c figure out which traps we don't need to
c worry about in this column
c
      pmax = 10.
      do j = 0001, JDIM
         pixo(j) = pixi(j)
         if (pixo(j).gt.pmax) pmax = pixo(j)
         enddo

c-----------------------------------------------------------
c
c go thru the traps one at a time (from highest to lowest q)
c and see when they get filled and emptied; adjust the
c pixels accordingly
c
      do W = Ws, 001, -1
         if (q_w(W).gt.pmax) goto 333

         ftrap =   0.0
         ttrap =  _TDIM_

         fcarry = fcarry0
         do j = J1, J2                                  ! go up the column pixel by pixel
            pix_1 = pixo(j)
            if (ttrap.ge.(_TDIM_).and.pix_1.lt.q_w(w)-1) goto 444
            if (pixo(j).ge.0) then
               pix_1 = pixo(j) + fcarry                 ! step 1, shuffle charge in
               fcarry = pix_1 - int(pix_1)
               pix_1  = int(pix_1)
               endif

            if (j.gt.J1) then
               if (pixf(j).lt.pixf(j-1))
     .             ftrap = pixf(j)/pixf(j-1)*ftrap
               endif

            padd_2 = 0.                                 ! step 2, release the charge
            if (ttrap.lt._TDIM_) then
               ttrap = ttrap + 1
               padd_2 = rprof_wt(w,ttrap)*ftrap
               endif

            padd_3 = 0.
            prem_3 = 0.
            if (pix_1.ge.q_w(w)) then
               prem_3 = dpde_w(w)/NITs*pixf(j)
               if (ttrap.lt._TDIM_)
     .             padd_3 = cprof_wt(w,ttrap)*ftrap
               ttrap = 0
               ftrap = prem_3
               endif
            pixo(j) = pixo(j) + padd_2 + padd_3 - prem_3
  444       continue
            enddo!j
  333    continue
         enddo!w

      return
      end

first re-write

      subroutine sim_colreadout_l_uvis_w(pixi,pixo,pixf,J1,J2,JDIM,
     .                                   q_w,dpde_w,NITs,
     .                                   tprof_w,rprof_wt,cprof_wt,Ws)
      implicit none

      integer JDIM
      real    pixi(JDIM)        ! input column array
      real    pixo(JDIM)        ! outout column array
      real    pixf(JDIM)        ! outout column array
      integer J1, J2

      integer    q_w(_WDIM_)    ! the 5 readout-parameter arrays
      real    dpde_w(_WDIM_)
      integer NITs
      integer tprof_w(_WDIM_)
      real    rprof_wt(_WDIM_,100)
      real    cprof_wt(_WDIM_,100)
      integer Ws

      integer j, t

      real  pix0

      real*8  ftrap
      integer ttrap

      integer w

      real pmax

      real    prel_1 ! incidental release of trap
      real    prel_2 ! flush release of trap
      real    pgrb_3 ! fill trap (amount grabbed)

      if (Ws.gt._WDIM_) stop 'Ws.gt._WDIM_'   ! more traps than we are allowed

c----------------------------------------
c
c figure out which traps we don't need to
c worry about in this column
c
      pmax = 10.
      do j = 0001, JDIM
         pixo(j) = pixi(j)
         if (pixo(j).gt.pmax) pmax = pixo(j)
         enddo

c-----------------------------------------------------------
c
c go thru the traps one at a time (from highest to lowest q)
c and see when they get filled and emptied; adjust the
c pixels accordingly
c
      do W = 001, Ws                                 ! I changed the order of this do loop

         if (q_w(W).gt.pmax) goto 777                ! if trap too high to matter, break

         ftrap =   0.0
         ttrap =  _TDIM_ + 1

         do j = J1, J2                               ! go up the column pixel by pixel

                                                         ! set up accounting of pixel value chgs
            prel_1 = 0.                              !    charge to be released                             
            prel_2 = 0.                              !    charge to be flushed out                               
            pgrb_3 = 0.                              !    charge to be trapped

            if (ttrap.le.(_TDIM_)) then              ! if the trap contains charge, let it
                ttrap = ttrap + 1                    !    release the appropriate number
                prel_1 = rprof_wt(w,ttrap)*ftrap     !    of electrons
                endif

            if (j.gt.J1) then                        ! if we have an inversion of the
                if (pixf(j).lt.pixf(j-1))            !    trap density (ie, a readout cosmic),
     .              ftrap = pixf(j)/pixf(j-1)*ftrap  !    then we don't want to flush too much
                endif                                

            if (pixo(j).ge.q_w(w)) then              ! filled/refilled trap here?
                if (ttrap.lt._TDIM_) then            !    need to flush?    
                    prel_2 = cprof_wt(w,ttrap)*ftrap !       amount to flush out of the trap
                    endif 
                ftrap  = dpde_w(w)/NITs*pixf(j)      !    amount to hold in the trap
                pgrb_3 = ftrap                       !    subtract that amount held from the pixel
                ttrap = 0                            !    reset the time-since-filled counter
                endif

            pixo(j) = pixo(j)                        ! make adjustment to the output pixel
     .              + prel_1                         !    add the trail emission 
     .              + prel_2                         !    flush the trap
     .              - pgrb_3                         !    fill the trap
            enddo!j
         enddo!w

  777 continue                                       ! come here if we break out of loop

      return
      end

and the latest.

c------------------------------------------------------------
c
c This is the workhorse subroutine; it simulates the readout
c of one column of pixels, given in array pixi(), and outputs 
c this to another array, pixo(), using a single iteration of 
c the CTE model at (1/NITs) intensity.  
c
c This routine can be called successively to do the transfer 
c in NITs steps.
c
      subroutine sim_colreadout_l_uvis_w(pixi,pixo,pixf,J1,J2,JDIM,
     .                                   q_w,dpde_w,NITs,
     .                                   tprof_w,rprof_wt,cprof_wt,Ws)
      implicit none

      integer JDIM                 ! number of pixels in column
      real    pixi(JDIM)           ! input column array
      real    pixo(JDIM)           ! outout column array
      real    pixf(JDIM)           ! scaling of model for each pixel
      integer J1, J2               ! bottom and top pixel in column

                                   ! the five readout-parameter arrays 
      integer    q_w(_WDIM_)       !    the charge level for trap#W
      real    dpde_w(_WDIM_)       !    the amt of charge this trap grabs
      integer NITs                 !    the num of iterations (dpde-->dpde/NITs)
      integer tprof_w(_WDIM_)      !    <apparently not used; supposed to be length of trail>
      real    rprof_wt(_WDIM_,100) !    the trap emission for t=T
      real    cprof_wt(_WDIM_,100) !    the amount left in trap after t=T emission
      integer Ws

      integer j                    ! pixel location up the column

      real*8  ftrap                ! total number of electrons in the trap
      integer ttrap                ! shifts since the trap was last filled

      integer w                    ! trap counter
      real    pmax                 ! max pixel value in the column; tells us the highest relevant trap numbrer
      integer Wf                   ! highest relevant trap number

      real    prel_1               ! amount in incidental release from trap
      real    prel_2               ! amount in flush release of trap
      real    pgrb_3               ! amount grabbed by filling trap 

c----------------------------------------
c
c bounds checking
c
      if (Ws.gt._WDIM_) stop 'Ws.gt._WDIM_' 

c----------------------------------------
c
c figure out which traps we don't need to
c worry about in this column
c
      pmax = 10.
      do j = 0001, JDIM
         pixo(j) = pixi(j)
         if (pixo(j).gt.pmax) pmax = pixo(j)
         enddo

c----------------------------------------
c
c figure the highest trap number we need to
c concern ourselves with
c
      Wf = 1
      do W = 1, Ws
         if (pmax.ge.q_w(W)) Wf = W
         enddo

c-----------------------------------------------------------
c
c go thru the traps one at a time (from highest to lowest q)
c and see when they get filled and emptied; adjust the
c pixel values accordingly
c
      do w = Wf, 1, -1                                 ! it's better to go backwards, now that it's easy

         ftrap =   0.0                                 ! initialize the flux in the trap to zero
         ttrap =  _TDIM_ + 1                           ! initialize the time-since-fluxh to the max

         do j = J1, J2                                 ! go up the column pixel by pixel

            if (j.gt.J1) then                          ! if we have an inversion of the the
                if (pixf(j).lt.pixf(j-1))              !    density (ie, a readout-cosmic issue),
     .              ftrap = pixf(j)/pixf(j-1)*ftrap    !    then we don't want to flush too much
                endif                                
                                                       ! set up accounting of pixel value chgs
            prel_1 = 0.                                !    charge to be released                             
            prel_2 = 0.                                !    charge to be flushed out                               
            pgrb_3 = 0.                                !    charge to be trapped

            if (pixo(j).ge.q_w(w)) then                ! filled/refilled trap#W here?
                  if (ttrap.le.(_TDIM_)) then          !    need to flush before filling?    
                      ttrap = ttrap + 1                !       increment time since filled
                      prel_1 = rprof_wt(w,ttrap)*ftrap !       flush out amount for this shift
                      prel_2 = cprof_wt(w,ttrap)*ftrap !       flush out the rest
                      endif 
                  ftrap  = dpde_w(w)/NITs*pixf(j)      !    amount to hold in the trap
                  pgrb_3 = ftrap                       !    subtract that amount held from the pixel
                  ttrap = 0                            !    reset the time-since-filled counter
               else                                    ! trap#W not filled this time
                  if (ttrap.le.(_TDIM_)) then          !    does trap contain charge?
                      ttrap = ttrap + 1                !       release the appropriate number
                      prel_1 = rprof_wt(w,ttrap)*ftrap !       of electrons
                      endif
               endif

            pixo(j) = pixo(j)                        ! make adjustment to the output pixel
     .              + prel_1                         !    add the trail emission 
     .              + prel_2                         !    flush the trap
     .              - pgrb_3                         !    fill the trap

            enddo!j

         enddo!w

      return
      end

[jamienoss]

@JayVB Hi, I've just noticed something that needs clarifying, please.

In your original Fortran and the current C implementation you compare ttrap < trailLength however, in your new version you compare ttrap <= trailLength. Was this intended?

[JayVB]

Good catch. I think in both of the instances, you are right that it should be "less-than" not "less-than-equal-to".

In practice, it isn't hurting us at all since the arrays have 100 elements along the trail, which is more room than the formal trail length, and there is formally no flux in the trail beyond the formal trail length. But you are right that this is unnecessarily sloppy.

Please do change it to less-than in both cases. I'll do the same in my code. Again, it should have null effect since our trail lengths are smaller than the maximum of 100 that I allow in the parameter files, but we may as well keep the code consistent.

[jamienoss]

Cool, will do. Cheers :)