Fixing bug in VoltageCalibration

[martin01021207]

Fixed the function for the waveform windows and the physical windows mappings. From a specific case to a general case, now the mapping subranges change depending on the input start_window. And also removed all the hard coding parts in that function.

[fschlueter]

Hey Martin, I still have to digest your code. It might be correct but I feel it is unnecessary complicated. In python you can do something as easy as this:

        samples_idx = (128 * starting_window_channel + numpy.arange(2048)) % 2048
        if starting_window_channel >= 16:
            samples_idx += 2048

Probably with c++ it will be a few more lines (i.e., you have to use a for loop), but with using the modulo operator % it should be possible to make the code much simpler and easier to read.

[fschlueter]

so converted to c++ I think something like this would work:

int start_sample = window_shift * 128;

for (int k = 0; k < 2048; k++) {
  int isamp = (start_sample + k) % 2048
  if (is2ndHalfWindows) {
    isamp += mattak::k::num_radiant_samples;
  }
}

For the old firm where something simular should be possible

[cozzyd]

Yes, for the new firmware, if start window is already in the range [0,31] I think it's just

int isamp = (start_window >=16) * 2048 + ( k + start_window * 128) % 2048;

For the old firmware I think it's

int isamp = (start_window/8) * 1024 + (k + start_window * 128) % 1024;

or something like that...

[martin01021207]

@fschlueter and @cozzyd I just fixed the function the second time following you guys' brilliant ideas, so now the code looks much cleaner!

[cozzyd]

I don't think this is quite right... it loops over iwindow but iwindow is never used?

[martin01021207]

Forgot to change the nested loop to just one loop for all samples... My bad!

[martin01021207]

But now the vectorization part below needs to be modified too since the loop structure is different than before...

[fschlueter]

@martin01021207 while you are working on this can you rename int N to something more explanatory (and avoid capital letter at the beginning)

[fschlueter]

@cozzyd, is there any benefit of using pointers for residTablePerSamp_adc params? We can pass the variables as reference right?

[fschlueter]

@martin01021207 while you are working on this can you rename int N to something more explanatory (and avoid capital letter at the beginning)

and maybe rename isamp to something like lab4_sample or something like that?

[fschlueter]

I also think there is still a bug for old_firmware == true. I do not see how the correct code selects sample in the second buffer, reading it I think it selects samples from the first buffer twice. Is that correct, can you double check? Also in the old firmware is the second starting window always equal the first starting window + 8, or are they independent from each other.

[cozzyd]

@cozzyd, is there any benefit of using pointers for residTablePerSamp_adc params? We can pass the variables as reference right?

well, it's currently an array that has decayed to a pointer, so there is little reason to use a reference there (it would be a reference to a pointer, unless you want to lose the array entirely...). There is a way to rewrite with std::vector that might be simpler, but we can do that later.

[martin01021207]

@martin01021207 while you are working on this can you rename int N to something more explanatory (and avoid capital letter at the beginning)

and maybe rename isamp to something like lab4_sample or something like that?

The N has been there so basically I just followed the convention, as you can see the input argument is N, and there is

if (!out) out = new double[N];

  if (N % mattak::k::radiant_window_size)
  {
    std::cerr << "Not multiple of window size!" << std::endl;
    return 0;
  }

something like that, but I can add something to it.

[martin01021207]

@cozzyd, is there any benefit of using pointers for residTablePerSamp_adc params? We can pass the variables as reference right?

I made residTablePerSamp_adc a pointer also the packed_aveResid_volt following the pointer structure const double *params = packed_fit_params + isamp * (fit_order+1); that Cosmin wrote, but if you were to change it to not using pointers at all, many things would also need to be changed in this script?

[cozzyd]

@martin01021207 while you are working on this can you rename int N to something more explanatory (and avoid capital letter at the beginning)

For some reason, I usually capitalize N even though I don't generally capitalize variables... I don't know why, I guess I feel like it looks better :).

the C++ convention I roughly (but not religiously) follow is:

class names: CapitalizedCamelCase method names: uncapitalizedCamelCase variable names: snake_case (but sometimes I don't do this for single-character names, like sometimes I use x/X for time-domain/fourier transform, etc.).

[martin01021207]

I also think there is still a bug for old_firmware == true. I do not see how the correct code selects sample in the second buffer, reading it I think it selects samples from the first buffer twice. Is that correct, can you double check? Also in the old firmware is the second starting window always equal the first starting window + 8, or are they independent from each other.

I've tested Cosmin's formula and it does what we want it to do. The start_window is not changing because it's directly from the input argument.

[fschlueter]

@martin01021207 I am still confused. I tried to replicated the algorithm in python. And I do not get the expected result.

def func(start_window, new=True):
    nSamplesPerGroup = 2048
    nWindowsPerGroup = 16

    window_order = []

    if not new:
        nSamplesPerGroup /= 2
        nWindowsPerGroup /= 2
        isamp_A = int(start_window / nWindowsPerGroup) * nSamplesPerGroup
    else:
        isamp_A = (start_window >= nWindowsPerGroup) * nSamplesPerGroup

    print("isamp_A", isamp_A)
    for i in range(2048):
        isamp_B = (i + start_window * 128) % nSamplesPerGroup

        isamp_lab4 = isamp_A + isamp_B

        window = int(isamp_lab4 / 128)
        if not len(window_order) or window_order[-1] != window:
            window_order.append(window)

    print(window_order)

if __name__ == "__main__":
    func(3, False)
    func(3, True)

Maybe I made a mistake, please appologize if I did. I also assume that this line in c++ isamp_A = (start_window / nWindowsPerGroup) * nSamplesPerGroup makes an integer division.

[fschlueter]

@cozzyd, is there any benefit of using pointers for residTablePerSamp_adc params? We can pass the variables as reference right?

well, it's currently an array that has decayed to a pointer, so there is little reason to use a reference there (it would be a reference to a pointer, unless you want to lose the array entirely...). There is a way to rewrite with std::vector that might be simpler, but we can do that later.

Ah I understand, this way we do not copy at all anything which is much faster, right? I thought we would create an array of the parameter per sample (and thus already copying it) and then pass is as reference (to not copy a second time - but I see now that not copying at all is probably desired here)

[martin01021207]

@martin01021207 I am still confused. I tried to replicated the algorithm in python. And I do not get the expected result.

def func(start_window, new=True):
    nSamplesPerGroup = 2048
    nWindowsPerGroup = 16

    window_order = []

    if not new:
        nSamplesPerGroup /= 2
        nWindowsPerGroup /= 2
        isamp_A = int(start_window / nWindowsPerGroup) * nSamplesPerGroup
    else:
        isamp_A = (start_window >= nWindowsPerGroup) * nSamplesPerGroup

    print("isamp_A", isamp_A)
    for i in range(2048):
        isamp_B = (i + start_window * 128) % nSamplesPerGroup

        isamp_lab4 = isamp_A + isamp_B

        window = int(isamp_lab4 / 128)
        if not len(window_order) or window_order[-1] != window:
            window_order.append(window)

    print(window_order)

if __name__ == "__main__":
    func(3, False)
    func(3, True)

Maybe I made a mistake, please appologize if I did. I also assume that this line in c++ isamp_A = (start_window / nWindowsPerGroup) * nSamplesPerGroup makes an integer division.

@fschlueter You were right about it, Felix. Your testing was more thorough than mine so that's why you found the issue but I didn't realize it, thank you for that! I've revised the formula in the function from your example script, so it should look like this, and it should now work the way we want:

def func(start_window, new=True):
    nSamplesPerGroup = 2048
    nWindowsPerGroup = 16

    window_order = []

    isamp_A = (start_window >= nWindowsPerGroup) * nSamplesPerGroup

    if not new:
        nSamplesPerGroup /= 2
        nWindowsPerGroup /= 2

    print("isamp_A", isamp_A)
    for i in range(2048):
        isamp_B = (i + start_window * 128) % nSamplesPerGroup

        if not new:
            isamp_B += (i >= nSamplesPerGroup) * nSamplesPerGroup

        isamp_lab4 = isamp_A + isamp_B

        window = int(isamp_lab4 / 128)

        if not len(window_order) or window_order[-1] != window:
            window_order.append(window)

    print(window_order)

[fschlueter]

@cozzyd The code for the old firmware just works if the second start window is first start window + 8. Is that always given?

[fschlueter]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

[fschlueter]

in the function adcTablePerSample is there a reason why only the pointer for the volt table but not the adc table is locally copied?

[fschlueter]

I was also wondering if we can optimise the calibration. I feel like there has to be a more efficient way. One ansatz I was thinking about, might it be possible that we not always create the residTablePerSamp_adc for the entire bias-scan range?

[cozzyd]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

[fschlueter]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

[fschlueter]

I was also wondering if we can optimise the calibration. I feel like there has to be a more efficient way. One ansatz I was thinking about, might it be possible that we not always create the residTablePerSamp_adc for the entire bias-scan range?

Another why would be to cache the adcTables for each channel and sample. What do you think @cozzyd, in python that make around 150 MB (getsizeof(np.random.normal(0, 1, (24, 4096, 196))) / 1e6)

[cozzyd]

I was also wondering if we can optimise the calibration. I feel like there has to be a more efficient way. One ansatz I was thinking about, might it be possible that we not always create the residTablePerSamp_adc for the entire bias-scan range?

Another why would be to cache the adcTables for each channel and sample. What do you think @cozzyd, in python that make around 150 MB (getsizeof(np.random.normal(0, 1, (24, 4096, 196))) / 1e6)

Yes, profiling suggests adc_table_per_sample is the dominant contributor to run time.

Caching the adc_table_per_sample is probably the right thing to do, or at least limiting the range computed to the relevant region (which would probably be around a factor of 5-10). Vectorizing evalPars would probably also get a factor of 3-4 or so in speed in calculating adc_table_per_sample (and help elsewhere too). We can use float math instead of double math to bring down the size to ~77 MB.

[cozzyd]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

voltTable[i] I believe is an ADC count? Right @martin01021207

[fschlueter]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

voltTable[i] I believe is an ADC count? Right @martin01021207

and resid_adc as well?

[martin01021207]

in the function adcTablePerSample is there a reason why only the pointer for the volt table but not the adc table is locally copied?

@fschlueter The biggest reason we calculate the average residual in the fitting way volt -> adc is that volt points are equidistant, so we want them to be in the x-axis. It means that the voltTable (x-axis) doesn't change, and the adcTable (y-axis) is y + Δy, i.e. y from the polynomial and Δy from the residual.

[martin01021207]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

@fschlueter This function is not wrong, because evalPars() can do either adc -> volt or volt -> adc depending on how you do the fitting, when there's no residuals you don't flip the axes so your fits are directly adc -> volt, then the fit parameters become coming from f(adc), this is the original calculation that Cosmin wrote; on the other hand when you do calculate the residuals, you flip the axes so your fits become volt -> adc in order to easily get the average residuals, so then your fit parameters become coming from f(volt).

[martin01021207]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

voltTable[i] I believe is an ADC count? Right @martin01021207

@cozzyd Here the voltTable[i] is just the volt points from the bias scan, and evalPars() is used to find the corresponding adc from the input volt, since the fits are done in the way volt -> adc in the isUsingResid == true case, all the fit parameters are from f(volt) not f(adc).

[martin01021207]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

Unless it has changed, it should be adc->volt

Than this function is wrong, or do I completely missunderstand things here:

static double* adcTablePerSample(int order, int npoints, const double * par, const double * resid_volt, const double * resid_adc)
{
  const double *voltTable = resid_volt;
  double *adcTable = new double[npoints];

  for (int i = 0; i < npoints; i++)
  {
    adcTable[i] = evalPars(voltTable[i], order, par) + resid_adc[i];
  }

  return adcTable;
}

voltTable[i] I believe is an ADC count? Right @martin01021207

and resid_adc as well?

@fschlueter The resid_adc here are literally just the average residuals (Δy) that we calculated after the fitting.

[martin01021207]

@martin01021207 and @cozzyd this part of the code (only used when isUsingResid == false) is wrong, right?

    else
    {
      out[i] = evalPars(adc, fit_order, params);
    }

evalPars is the function Volt -> ADC and not the other way around?

@fschlueter That's the original voltage calibration when there's no residuals, so the fits are done directly adc -> volt, adc would be the x-axis in this case. Only when we want the average residual we flip the axes, so volt becomes the x-axis.

evalPars() doesn't have to be volt -> adc, it can be used as adc -> volt, depends on what parameters go into the function.

[fschlueter]

in the function adcTablePerSample is there a reason why only the pointer for the volt table but not the adc table is locally copied?

@fschlueter The biggest reason we calculate the average residual in the fitting way volt -> adc is that volt points are equidistant, so we want them to be in the x-axis. It means that the voltTable (x-axis) doesn't change, and the adcTable (y-axis) is y + Δy, i.e. y from the polynomial and Δy from the residual.

My question was purely a c++ one, not a physics one.

[fschlueter]

Okay, I will go on and merge this PR, the remaining c++ questions can be solved in a different PR.

[martin01021207]

in the function adcTablePerSample is there a reason why only the pointer for the volt table but not the adc table is locally copied?

@fschlueter The biggest reason we calculate the average residual in the fitting way volt -> adc is that volt points are equidistant, so we want them to be in the x-axis. It means that the voltTable (x-axis) doesn't change, and the adcTable (y-axis) is y + Δy, i.e. y from the polynomial and Δy from the residual.

My question was purely a c++ one, not a physics one.

Those local copies are not necessary now, same answer as my last reply. And I've removed them.