SuperFunction() does not work properly when more than two functions are used

[StephenLesage]

SuperFunction() does not work properly when more than two functions are used. Here is example code to reproduce the issue:

import matplotlib.pyplot as plt

from gdt.missions.fermi.gbm.phaii import GbmPhaii
from gdt.missions.fermi.gbm.response import GbmRsp2
from gdt.missions.fermi.gbm.collection import GbmDetectorCollection
from gdt.core.background.fitter import BackgroundFitter
from gdt.core.background.binned import Polynomial
from gdt.core.spectra.fitting import SpectralFitterPstat
from gdt.core.spectra.functions import SmoothlyBrokenPowerLaw, PowerLaw, GaussLine
from gdt.core.plot.model import ModelFit

if __name__ == "__main__":

    b0_phaii = GbmPhaii.open('./data/glg_cspec_b0_bn221009553_v00.pha')
    b1_phaii = GbmPhaii.open('./data/glg_cspec_b1_bn221009553_v00.pha')
    n4_phaii = GbmPhaii.open('./data/glg_cspec_n4_bn221009553_v00.pha')
    n6_phaii = GbmPhaii.open('./data/glg_cspec_n6_bn221009553_v00.pha')
    n8_phaii = GbmPhaii.open('./data/glg_cspec_n8_bn221009553_v00.pha')
    cspecs = GbmDetectorCollection.from_list([b0_phaii, b1_phaii, n4_phaii, n6_phaii, n8_phaii])

    b0_rsp2 = GbmRsp2.open('./data/glg_cspec_b0_bn221009553_v13.rsp2')
    b1_rsp2 = GbmRsp2.open('./data/glg_cspec_b1_bn221009553_v13.rsp2')
    n4_rsp2 = GbmRsp2.open('./data/glg_cspec_n4_bn221009553_v13.rsp2')
    n6_rsp2 = GbmRsp2.open('./data/glg_cspec_n6_bn221009553_v13.rsp2')
    n8_rsp2 = GbmRsp2.open('./data/glg_cspec_n8_bn221009553_v13.rsp2')
    rsps = GbmDetectorCollection.from_list([b0_rsp2, b1_rsp2, n4_rsp2, n6_rsp2, n8_rsp2])

    bkgd_times = [(-150.0, -2.0), (1620.0, 2370.0)]
    poly_orders = [4, 4, 3, 4, 4]
    src_time = (310, 320)
    erange_nai = (45.0, 900.0)
    erange_bgo = (400.0, 39_000.0)

    backfitters = [BackgroundFitter.from_phaii(cspec, Polynomial, time_ranges=bkgd_times) for cspec in cspecs]
    backfitters = GbmDetectorCollection.from_list(backfitters, dets=cspecs.detector())

    [backfitter.fit(order=poly_order) for backfitter, poly_order in zip(backfitters, poly_orders)]

    bkgds = backfitters.interpolate_bins(cspecs.data()[0].tstart, cspecs.data()[0].tstop)
    bkgds = GbmDetectorCollection.from_list(bkgds, dets=cspecs.detector())

    bkgd_specs = bkgds.integrate_time(*src_time)

    phas = cspecs.to_pha(time_ranges=src_time, nai_kwargs={'energy_range':erange_nai}, bgo_kwargs={'energy_range':erange_bgo})
    rsps_interp = [rsp.interpolate(pha.tcent) for rsp, pha in zip(rsps, phas)]

    print('\n----- SBPL + PL Fit -----')

    SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw()
    SBPL_PL.default_values = [0.06, 12.0, 9.0, 0.5, -4.4, 100.0, 0.06, -1.7, 87.0]

    specfitter = SpectralFitterPstat(phas, bkgds.to_list(), rsps_interp, method='L-BFGS-B')
    specfitter.fit(SBPL_PL, options={'maxiter': 10_000})

    modelplot = ModelFit(fitter=specfitter)
    plt.savefig('./SBPL_PL_cnts.png')

    modelplot = ModelFit(fitter=specfitter, view='nufnu')
    modelplot.nufnu_spectrum(num_samples=500)
    modelplot.ax.grid(which='both')
    plt.savefig('./SBPL_PL_nuFnu.png')

    print('\n----- SBPL + PL + Gline Fit -----')

    # SBPL_PL_Gline = SmoothlyBrokenPowerLaw()+PowerLaw()+GaussLine()
    # SBPL_PL_Gline = (SmoothlyBrokenPowerLaw()+PowerLaw())+GaussLine()
    SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw()
    SBPL_PL_Gline = SBPL_PL + GaussLine()

    SBPL_PL_Gline.default_values = [0.06, 12.0, 9.0, 0.5, -4.4, 100.0, 0.06, -1.7, 87.0, 0.11, 8_975.0, 2_649.0]

    specfitter = SpectralFitterPstat(phas, bkgds.to_list(), rsps_interp, method='L-BFGS-B')
    specfitter.fit(SBPL_PL_Gline, options={'maxiter': 10_000})

    modelplot = ModelFit(fitter=specfitter)
    plt.savefig('./SBPL_PL_Gline_cnts.png')

    modelplot = ModelFit(fitter=specfitter, view='nufnu')
    modelplot.nufnu_spectrum(num_samples=500)
    modelplot.ax.grid(which='both')
    plt.savefig('./SBPL_PL_Gline_nuFnu.png')

Here is the failure trace:

Traceback (most recent call last):
  File "/Users/stephenlesage/Fermi-GBM/projects/gaussian_line/python_sims/test.py", line 87, in <module>
    modelplot = ModelFit(fitter=specfitter, view='nufnu')
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 78, in __init__
    self.set_fit(fitter, resid=resid)
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 215, in set_fit
    self.nufnu_spectrum()
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 173, in nufnu_spectrum
    self._plot_spectral_model(**kwargs)
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 288, in _plot_spectral_model
    label=comps[i], color=self.colors[i+1],
IndexError: list index out of range

I tried the following syntax suggestions and got the same error message: 1) SBPL_PL_Gline = SmoothlyBrokenPowerLaw()+PowerLaw()+GaussLine() 2) SBPL_PL_Gline = (SmoothlyBrokenPowerLaw()+PowerLaw())+GaussLine() 3) SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw() SBPL_PL_Gline = SBPL_PL + GaussLine()

It appears to be an issue in model.py, specifically that the number of components "num_comp" does not reflect the number of function components "comps" being fed to SuperFunction(). Implementing print statements reveals:

num_comp = 3
comps = ['SmoothlyBrokenPowerLaw+PowerLaw', 'GaussLine']

I can provide the data files to run the example code if necessary.

[AdamGoldstein-USRA]

This looks like a problem in updating the names attribute in SuperFunction. This is really a bug in gdt-core, so I will open up an issue there and reference this one.

[AdamGoldstein-USRA]

@StephenLesage the fix was implemented in PR USRA-STI/gdt-core#35 and merged into main. You should be able to pip install the main branch from the repo and verify this solves the issue.