ROEChargeInjection does not work via wrapper

[Jammy2211]

If I pass ROEChargeInjection to the wrapper, it does not calls it set_express_matrix_from_rows_and_express and set_store_trap_states_matrix methods, instead reverting to the methods of the default ROE class.

I have put a print(roe->type) before the following line in add_cti.cpp`:

roe->set_express_matrix_from_rows_and_express(n_rows, express, offset);

And the type is indeed type 1, indicating that the Python wrapper appears to be passing the correct type.

My best guess is that the inheritance structure is failing to overwrite the methods correctly?

[Jammy2211]

I have got it working via the following hacks, which are on the branch https://github.com/jkeger/arctic/tree/hack_good:

class ROEChargeInjection : public ROE {
   public:
    ROEChargeInjection(
        std::valarray<double>& dwell_times = dwell_times_default,
        bool empty_traps_between_columns = true,
        bool force_release_away_from_readout = true,
        bool use_integer_express_matrix = false);
    ~ROEChargeInjection(){};

    void set_express_matrix_from_rows_and_express(
        int n_rows, int express = 0, int offset = 0);
    void set_store_trap_states_matrix();
};

Is now:

class ROEChargeInjection : public ROE {
   public:
    ROEChargeInjection(
        std::valarray<double>& dwell_times = dwell_times_default,
        bool empty_traps_between_columns = true,
        bool force_release_away_from_readout = true,
        bool use_integer_express_matrix = false);
    ~ROEChargeInjection(){};

    void set_express_matrix_from_rows_and_express(
        int n_rows, int express = 0, int offset = 0);
    void set_express_matrix_from_rows_and_express_ci(
        int n_rows, int express = 0, int offset = 0);
    void set_store_trap_states_matrix();
    void set_store_trap_states_matrix_ci();
};

Anddddddddddddddddddd

// ========
// ROEChargeInjection::
// ========
/*
    Class ROEChargeInjection.
    Modified ROE for charge injection modes.
    Instead of charge starting in each pixel and moving different distances to
    the readout register, for charge injection the electrons are directly
    created at the far end of the CCD, then are all transferred the same number
    of times through the full image of pixels to the readout register.
    Parameters
    ----------
    Same as ROE, but empty_traps_for_first_transfers is automatically false,
    since only the single leading pixel of the charge injection in each column
    will see the untouched traps.
*/
ROEChargeInjection::ROEChargeInjection(
    std::valarray<double>& dwell_times, bool empty_traps_between_columns,
    bool force_release_away_from_readout, bool use_integer_express_matrix)
    : ROE(dwell_times, empty_traps_between_columns, false,
          force_release_away_from_readout, use_integer_express_matrix) {

    type = roe_type_charge_injection;
}

/*
    See ROE::set_express_matrix_from_rows_and_express().
    For charge injection, all charges are clocked the same number of times
    through all the pixels to the readout register.
*/
void ROEChargeInjection::set_express_matrix_from_rows_and_express(
    int n_rows, int express, int offset) {

    int n_transfers = n_rows + offset;

    // Set default express to all transfers, and check no larger
    if (express == 0)
        express = n_transfers;
    else
        express = std::min(express, n_transfers);

    // Compute the multiplier factors
    double max_multiplier = (double)n_transfers / express;
    if (use_integer_express_matrix) max_multiplier = ceil(max_multiplier);

    // Initialise an array with enough pixels to contain the supposed image,
    // including offset
    std::valarray<double> tmp_express_matrix(max_multiplier, express * n_transfers);

    // Adjust integer multipliers to correct the total number of transfers
    if ((use_integer_express_matrix) && (n_transfers % express != 0)) {
        double current_n_transfers;
        double reduced_multiplier;

        for (int express_index = express - 1; express_index >= 0; express_index--) {
            // Count the current number of transfers for this pixel
            current_n_transfers = 0.0;
            for (int i = 0; i <= express_index; i++) {
                current_n_transfers += tmp_express_matrix[i * n_transfers];
            }

            // Reduce the multipliers until no longer have too many transfers
            if (current_n_transfers <= n_transfers) break;

            reduced_multiplier =
                std::max(0.0, max_multiplier + n_transfers - current_n_transfers);

            tmp_express_matrix[std::slice(
                express_index * n_transfers, n_transfers, 1)] = reduced_multiplier;
        }
    }

    n_express_passes = express;

    // Finished if no offset
    if (offset == 0) {
        express_matrix = tmp_express_matrix;
        return;
    }

    // Remove the offset (which is not represented in the image pixels)
    std::valarray<double> express_matrix_trim(0.0, n_express_passes * n_rows);

    // Copy the post-offset slices of each row
    for (int express_index = 0; express_index < n_express_passes; express_index++) {
        express_matrix_trim[std::slice(express_index * n_rows, n_rows, 1)] =
            tmp_express_matrix[std::slice(express_index * n_transfers, n_rows, 1)];
    }

    express_matrix = express_matrix_trim;
}

/*
    See ROE::set_store_trap_states_matrix().
    For charge injection, the first charge cloud in each column will always
    encounter empty traps in every new pixel. So no need to store trap states
    between transfers.
*/
void ROEChargeInjection::set_store_trap_states_matrix() {
    store_trap_states_matrix = std::valarray<bool>(false, express_matrix.size());
}

// ========
// ROETrapPumping::
// ========
/*
    Class ROETrapPumping.
    Modified ROE for trap pumping (AKA pocket pumping) modes.
    Instead of clocking the charge in all pixels towards the readout register,
    trap pumping shifts the charge back and forth, to end up in the same place
    they began. If one pixel (and phase) contains traps and its neighbours do
    not, then a charge dipole can be created as charge is captured and released
    asymmetrically in the active and adjacent pixel. This allows direct study of
    the trap species in a CCD.
    The number of phases is assumed to be half the number of steps in the clock
    sequence, which must be even.
    Currently, this algorithm assumes that only a single pixel is active and
    contains traps.
    Parameters
    ----------
    Same as ROE, but empty_traps_between_columns is automatically true, and
    force_release_away_from_readout is automatically false.
    n_pumps : int
        The number of times the charge is pumped back and forth.
    The diagram below illustrates the steps in the clocking sequence produced
    by ROE::set_clock_sequence() in this mode, for three phases. The first three
    steps are the same as the standard case. However, now there are three
    additional steps in which the charge cloud is shifted one step further into
    the next pixel, but then back to where it began in the original pixel,
    instead of continuing on towards the readout register.
    This means that, unlike in the standard case, the traps in this pixel can
    capture charge that originated in a different pixel, as shown in step 3,
    where the high phase in pixel p contains the charge cloud that started in
    pixel p+1.
    Three phases
    ============
                #     Pixel p-1      #       Pixel p      #     Pixel p+1      #
    Step         Phase2 Phase1 Phase0 Phase2 Phase1 Phase0 Phase2 Phase1 Phase0
    0           +             +------+             +------+             +------+
    Capture from|             |      |             |   p  |             |      |
    Release to  |             |      |  p-1     p  |   p  |             |      |
                +-------------+      +-------------+      +-------------+      +
    1                  +------+             +------+             +------+
    Capture from       |      |             |   p  |             |      |
    Release to         |      |          p  |   p  |   p         |      |
                -------+      +-------------+      +-------------+      +-------
    2           +------+             +------+             +------+             +
    Capture from|      |             |   p  |             |      |             |
    Release to  |      |             |   p  |   p     p+1 |      |             |
                +      +-------------+      +-------------+      +-------------+
    3           +             +------+             +------+             +------+
    Capture from|             |      |             |  p+1 |             |      |
    Release to  |             |      |   p     p+1 |  p+1 |             |      |
                +-------------+      +-------------+      +-------------+      |
    4           +------+             +------+             +------+             +
    Capture from|      |             |   p  |             |      |             |
    Release to  |      |             |   p  |   p     p+1 |      |             |
                +      +-------------+      +-------------+      +-------------+
    5                  +------+             +------+             +------+
    Capture from       |      |             |   p  |             |      |
    Release to         |      |          p  |   p  |   p         |      |
                -------+      +-------------+      +-------------+      +-------
    Below are corresponding illustrations for two and four phases. As in the
    standard case, an even number of phases leads to released charge being split
    in one phase each step between twp pixels.
    Two phases
    ==========
                  #  Pixel p-1  #   Pixel p   #  Pixel p+1  #
    Step           Phase1 Phase0 Phase1 Phase0 Phase1 Phase0
    0             +      +------+      +------+      +------+
    Capture from  |      |      |      |   p  |      |      |
    Release to    |      |      | p-1&p|   p  |      |      |
                  +------+      +------+      +------+      +
    1             +------+      +------+      +------+      +
    Capture from  |      |      |   p  |      |      |      |
    Release to    |      |      |   p  | p&p+1|      |      |
                  +      +------+      +------+      +------+
    2             +      +------+      +------+      +------+
    Capture from  |      |      |      |  p+1 |      |      |
    Release to    |      |      | p&p+1|  p+1 |      |      |
                  +------+      +------+      +------+      +
    3             +------+      +------+      +------+      +
    Capture from  |      |      |   p  |      |      |      |
    Release to    |      |      |   p  | p&p+1|      |      |
                  +      +------+      +------+      +------+
    Four phases
    ===========
               Pixel p-1      #          Pixel p          #         Pixel p+1
    Step         Phase1 Phase0 Phase3 Phase2 Phase1 Phase0 Phase3 Phase2 Phase1
    0                  +------+                    +------+                    +
    Capture from       |      |                    |   p  |                    |
    Release to         |      |  p-1   p-1&p    p  |   p  |                    |
                -------+      +--------------------+      +--------------------+
    1           +------+                    +------+                    +------+
    Capture from|      |                    |   p  |                    |      |
    Release to  |      |        p-1&p    p  |   p  |   p                |      |
                +      +--------------------+      +--------------------+      +
    2           +                    +------+                    +------+
    Capture from|                    |   p  |                    |      |
    Release to  |                p   |   p  |   p    p&p+1       |      |
                +--------------------+      +--------------------+      +-------
    3                         +------+                    +------+
    Capture from              |   p  |                    |      |
    Release to                |   p  |   p    p&p+1  p+1  |      |
                --------------+      +--------------------+      +--------------
    4                  +------+                    +------+                    +
    Capture from       |      |                    |  p+1 |                    |
    Release to         |      |   p    p&p+1   p+1 |  p+1 |                    |
                -------+      +--------------------+      +--------------------+
    5                         +------+                    +------+
    Capture from              |   p  |                    |      |
    Release to                |   p  |   p    p&p+1  p+1  |      |
                --------------+      +--------------------+      +--------------
    6           +                    +------+                    +------+
    Capture from|                    |   p  |                    |      |
    Release to  |                p   |   p  |   p    p&p+1       |      |
                +--------------------+      +--------------------+      +-------
    7           +------+                    +------+                    +------+
    Capture from|      |                    |   p  |                    |      |
    Release to  |      |        p-1&p    p  |   p  |   p                |      |
                +      +--------------------+      +--------------------+      +
*/
ROETrapPumping::ROETrapPumping(
    std::valarray<double>& dwell_times, int n_pumps,
    bool empty_traps_for_first_transfers, bool use_integer_express_matrix)
    : ROE(dwell_times, true, empty_traps_for_first_transfers, false,
          use_integer_express_matrix) {

    type = roe_type_trap_pumping;
    this->n_pumps = n_pumps;

    if (n_steps % 2 != 0)
        error("The number of steps for trap pumping (%d) must be even", n_steps);
    n_phases = n_steps / 2;
}

/*
    See ROE::set_express_matrix_from_rows_and_express().
    For trap pumping, instead of charge being transferred from pixel to pixel
    until readout, only the back-and-forth pumping clock sequence is repeated a
    number of times for the active pixel(s) containing charge.
    Currently, this algorithm assumes that only a single pixel is active and
    contains traps. So, rather than the number of pixels and offset (which must
    be 1 and 0) the number of express passes is controlled by the number of
    pumps back and forth.
    Conveniently, the required express multipliers are the same as the ones for
    the pixel furthest from readout with standard clocking.
*/
void ROETrapPumping::set_express_matrix_from_rows_and_express(
    int n_rows, int express, int offset) {

    if (offset != 0) error("Trap pumping requires the offset (%d) to be 0", offset);

    // Set default express to all transfers, and check no larger
    if (express == 0)
        express = n_pumps;
    else
        express = std::min(express, n_pumps);

    // Start with the standard express matrix for n_transfers = n_pumps
    ROE::set_express_matrix_from_rows_and_express(n_pumps, express, offset);

    // Extract the relevant express multipliers of the final pixel
    std::valarray<double> tmp_col(0.0, n_express_passes);
    for (int express_index = 0; express_index < n_express_passes; express_index++) {
        tmp_col[express_index] = express_matrix[n_pumps - 1 + express_index * n_pumps];
    }

    // Extract the non-zero elements if doing first transfers separately
    if ((empty_traps_for_first_transfers) && (express < n_pumps)) {
        std::valarray<double> tmp_col_2 = tmp_col[tmp_col != 0.0];

        n_express_passes = express + 1;
        tmp_col.resize(n_express_passes);
        // Set the non-zero elements, which won't include the final entry for
        // mismatched integer multipliers
        if ((use_integer_express_matrix) && (n_pumps % express != 0)) {
            // (Using slice here doesn't compile on mac for some odd reason...)
            for (int i = 0; i < n_express_passes - 1; i++) tmp_col[i] = tmp_col_2[i];

            // Put back the final zero
            tmp_col[n_express_passes - 1] = 0.0;
        } else
            tmp_col = tmp_col_2;
    }

    express_matrix.resize(n_rows * n_express_passes);
    // Set multipliers for all rows, even though only one row will be active and
    // actually used
    for (int row_index = 0; row_index < n_rows; row_index++) {
        express_matrix[std::slice(row_index, n_express_passes, n_rows)] = tmp_col;
    }
}

/*
    See ROE::set_store_trap_states_matrix().
    For trap pumping, trap states must be stored after every pump sequence of
    the same trap, until the final pump.
*/

Is now:

// ========
// ROEChargeInjection::
// ========
/*
    Class ROEChargeInjection.

    Modified ROE for charge injection modes.

    Instead of charge starting in each pixel and moving different distances to
    the readout register, for charge injection the electrons are directly
    created at the far end of the CCD, then are all transferred the same number
    of times through the full image of pixels to the readout register.

    Parameters
    ----------
    Same as ROE, but empty_traps_for_first_transfers is automatically false,
    since only the single leading pixel of the charge injection in each column
    will see the untouched traps.
*/
ROEChargeInjection::ROEChargeInjection(
    std::valarray<double>& dwell_times, bool empty_traps_between_columns,
    bool force_release_away_from_readout, bool use_integer_express_matrix)
    : ROE(dwell_times, empty_traps_between_columns, false,
          force_release_away_from_readout, use_integer_express_matrix) {

    type = roe_type_charge_injection;
}

/*
    See ROE::set_express_matrix_from_rows_and_express().

    For charge injection, all charges are clocked the same number of times
    through all the pixels to the readout register.
*/
void ROEChargeInjection::set_express_matrix_from_rows_and_express(
    int n_rows, int express, int offset) {

    int n_transfers = n_rows + offset;

    // Set default express to all transfers, and check no larger
    if (express == 0)
        express = n_transfers;
    else
        express = std::min(express, n_transfers);

    // Compute the multiplier factors
    double max_multiplier = (double)n_transfers / express;
    if (use_integer_express_matrix) max_multiplier = ceil(max_multiplier);

    // Initialise an array with enough pixels to contain the supposed image,
    // including offset
    std::valarray<double> tmp_express_matrix(max_multiplier, express * n_transfers);

    // Adjust integer multipliers to correct the total number of transfers
    if ((use_integer_express_matrix) && (n_transfers % express != 0)) {
        double current_n_transfers;
        double reduced_multiplier;

        for (int express_index = express - 1; express_index >= 0; express_index--) {
            // Count the current number of transfers for this pixel
            current_n_transfers = 0.0;
            for (int i = 0; i <= express_index; i++) {
                current_n_transfers += tmp_express_matrix[i * n_transfers];
            }

            // Reduce the multipliers until no longer have too many transfers
            if (current_n_transfers <= n_transfers) break;

            reduced_multiplier =
                std::max(0.0, max_multiplier + n_transfers - current_n_transfers);

            tmp_express_matrix[std::slice(
                express_index * n_transfers, n_transfers, 1)] = reduced_multiplier;
        }
    }

    n_express_passes = express;

    // Finished if no offset
    if (offset == 0) {
        express_matrix = tmp_express_matrix;
        return;
    }

    // Remove the offset (which is not represented in the image pixels)
    std::valarray<double> express_matrix_trim(0.0, n_express_passes * n_rows);

    // Copy the post-offset slices of each row
    for (int express_index = 0; express_index < n_express_passes; express_index++) {
        express_matrix_trim[std::slice(express_index * n_rows, n_rows, 1)] =
            tmp_express_matrix[std::slice(express_index * n_transfers, n_rows, 1)];
    }

    express_matrix = express_matrix_trim;
}

void ROEChargeInjection::set_express_matrix_from_rows_and_express_ci(
    int n_rows, int express, int offset) {

    int n_transfers = n_rows + offset;

    // Set default express to all transfers, and check no larger
    if (express == 0)
        express = n_transfers;
    else
        express = std::min(express, n_transfers);

    // Compute the multiplier factors
    double max_multiplier = (double)n_transfers / express;
    if (use_integer_express_matrix) max_multiplier = ceil(max_multiplier);

    // Initialise an array with enough pixels to contain the supposed image,
    // including offset
    std::valarray<double> tmp_express_matrix(max_multiplier, express * n_transfers);

    // Adjust integer multipliers to correct the total number of transfers
    if ((use_integer_express_matrix) && (n_transfers % express != 0)) {
        double current_n_transfers;
        double reduced_multiplier;

        for (int express_index = express - 1; express_index >= 0; express_index--) {
            // Count the current number of transfers for this pixel
            current_n_transfers = 0.0;
            for (int i = 0; i <= express_index; i++) {
                current_n_transfers += tmp_express_matrix[i * n_transfers];
            }

            // Reduce the multipliers until no longer have too many transfers
            if (current_n_transfers <= n_transfers) break;

            reduced_multiplier =
                std::max(0.0, max_multiplier + n_transfers - current_n_transfers);

            tmp_express_matrix[std::slice(
                express_index * n_transfers, n_transfers, 1)] = reduced_multiplier;
        }
    }

    n_express_passes = express;

    // Finished if no offset
    if (offset == 0) {
        express_matrix = tmp_express_matrix;
        return;
    }

    // Remove the offset (which is not represented in the image pixels)
    std::valarray<double> express_matrix_trim(0.0, n_express_passes * n_rows);

    // Copy the post-offset slices of each row
    for (int express_index = 0; express_index < n_express_passes; express_index++) {
        express_matrix_trim[std::slice(express_index * n_rows, n_rows, 1)] =
            tmp_express_matrix[std::slice(express_index * n_transfers, n_rows, 1)];
    }

    express_matrix = express_matrix_trim;
}

/*
    See ROE::set_store_trap_states_matrix().

    For charge injection, the first charge cloud in each column will always
    encounter empty traps in every new pixel. So no need to store trap states
    between transfers.
*/
void ROEChargeInjection::set_store_trap_states_matrix() {
    store_trap_states_matrix = std::valarray<bool>(false, express_matrix.size());
}
void ROEChargeInjection::set_store_trap_states_matrix_ci() {
    store_trap_states_matrix = std::valarray<bool>(false, express_matrix.size());
}

[rjmassey]

Great, thanks. If this works, please merge it.
We've got all the mac compilation issues fixed on /dev/ branch, and are very soon going to pull everything into master. Your timing is perfect.

[Jammy2211]

The inheritance structure of the original code was broken, my solution involved a lot of copy and pasting doing bad thing...

Happy to merge it but its not good code and I would assume other parts of the code might be broken?

[rjmassey]

Oh, OK I see - sorry. This might be one to ask Jacob about. Could be really easy for him to fix.

[jkeger]

I should be able to look at this during next week if that's okay.

[jkeger]

I don't think the problem is happening for me? Sorry it took me a while to check it out. And perhaps I'm just missing something. But I'm on the latest dev branch (I made a quick new commit but that's only to avoid an irrelevant compiler warning), and as a simple example, in test/test_arctic.py I swapped line 795 cti.ROE( for cti.ROEChargeInjection(. I then added in src/roe.cpp:
printf("### ROE::set_express_matrix_from_rows_and_express() \n");
printf("### ROEChargeInjection::set_express_matrix_from_rows_and_express() \n");
printf("### ROETrapPumping::set_express_matrix_from_rows_and_express() \n");
at the start of their respective functions, on lines 203, 568, and 803. The output prints the desired ROEChargeInjection:: version.

So can you check that on yours, making sure you're not on an older branch, and if it is broken there or breaks with a different test case then send me a MWE and I'll take another look :)

[Jammy2211]

Definitely not fixed on the master branch.

I am running this script, and have checked via my hack that the output likelihood should give -561207.1121553452 for ROEChargeInjection. When I use the master branch it gives -655871.4574066445, which is the value I got using the ROE.

https://github.com/Jammy2211/autocti_workspace_test/blob/main/imaging_ci/profiling/likelihood/parallel_x3.py

The methods of void ROEChargeInjection are never called or used irrespective of the Python input.

[jkeger]

I'm confused, why are you using master, the latest version which is where I'm saying it seems to work is dev?

[Jammy2211]

master and dev look in sync (I think Richard maybe did a merge?) to me, but I will try dev in a bit.

[jkeger]

Oh okay, perhaps he did. Could you try doing something more explicit like my check above? I don't know what those numbers you printed mean. And if there is still a problem, I have a conference this coming week but if you could make a more minimal example then I can try to have a look.

[rjmassey]

I merged them today, because everything else was working. Having said that, I also just found a bug, in that charge is not conserved in an image - either the FPR is too small or the EPER too big (which might explain James’s amazing recent results). I’m investigating that separately/at the moment. On 2 Mar 2022, at 15:30, Jacob Kegerreis ***@***.******@***.***>> wrote: [EXTERNAL EMAIL] Oh okay, perhaps he did. Could you try doing something more explicit like my check above? I don't know what those numbers you printed mean. And if there is still a problem, I have a conference this coming week but if you could make a more minimal example then I can try to have a look. — Reply to this email directly, view it on GitHub<#19 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AE5BZNS5PRMQBQSFXQGPC4TU56CRRANCNFSM5K4AF7XQ>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you commented.Message ID: ***@***.***>

[Jammy2211]

For me, the following script gives 0.0 at the bottom, meaning that the addition of CTI is identical for the two ROE's:

import numpy as np
from arcticpy.src import cti
import autocti as ac


image_pre_cti = np.array([[1.0, 1.0, 1.0],
                         [0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0]])

parallel_trap_0 = ac.TrapInstantCapture(density=0.07275, release_timescale=0.8)

parallel_trap_list = [parallel_trap_0]

parallel_ccd = ac.CCD(
    well_fill_power=0.58, well_notch_depth=0.0, full_well_depth=200000.0
)


image_post_cti_0 = cti.add_cti(
    image=image_pre_cti,
    parallel_traps=parallel_trap_list,
    parallel_ccd=parallel_ccd,
    parallel_roe=ac.ROE()
)


image_post_cti_1 = cti.add_cti(
    image=image_pre_cti,
    parallel_traps=parallel_trap_list,
    parallel_ccd=parallel_ccd,
    parallel_roe=ac.ROEChargeInjection()
)


print(np.max(abs(image_post_cti_0 - image_post_cti_1)))

Do you get the same behaviour?

[jkeger]

Okay thanks that looks clear, I'll have a look when I can, hopefully before the conference starts, but let me know if it's urgent.

[Jammy2211]

Got my hacked branch so not urgent :)

[rjmassey]

This currently works perfectly on dev branch. But only because no classes are inherited from anywhere else (I think). For example, ROEChargeInjection() is its own class, and has nothing to do with ROE(). roe.cpp certainly looks very different from the snippet pasted above. Jacob, did you fix that? It works, but is it ideal? I can try to be more clever if that would be worthwhile.