pairs list is empty for do_hough_transform()

[megan-the-astronomer]

Hello, I'm working through the quickstart guide and getting stuck at every step that calls c.do_hough_transform(). For example,

---

ValueError Traceback (most recent call last)
Cell In[18], line 1
----> 1 c.do_hough_transform()

File /opt/anaconda3/envs/astr230/lib/python3.12/site-packages/rascal/calibrator.py:1545, in Calibrator.do_hough_transform(self, brute_force)
1543 def do_hough_transform(self, brute_force=False):
-> 1545 if self.pairs == []:
1547 logging.warning("pairs list is empty. Try generating now.")
1548 self._generate_pairs()

ValueError: operands could not be broadcast together with shapes (746,2) (0,)

I'm using the development branch downloaded today (2 Feb 2024). Can you please advise?

[cylammarco]

@jveitchmichaelis ping

[jveitchmichaelis]

Hi Megan,

Thanks for trying out Rascal, it'd really help if you could share the data you're trying to calibrate (feel free to email if it's sensitive), along with the script/notebook + config you've been using (eg lamp type, grating/instrument details, etc). Then we can poke around around figure out what's up!

Best
Josh

[megan-the-astronomer]

Thanks for following up! I'd like to use Rascal for my observational astronomy lab class. Our 16" campus telescope has a Lhires III Shelyak spectrograph - full details on the campus observatory webpage: https://www.ruf.rice.edu/~ruco/users.html

For the low resolution spectra, we end up stitching three separate lamp frames together to get good coverage across the whole wavelength range. It's that combined spectrum that we use to determine the wavelength solution. I've put a jupyter notebook showing the stitching and attempted wavelength solution, along with the resulting combined lamp spectrum and our bespoke linelist here: https://github.com/megan-the-astronomer/sandbox

We are using:
astropy 6.0.0
numpy 1.26.3
matplotlib 3.8.2
scipy 1.12.0
rascal 0.3.4

Have I forgotten anything?

Thanks for taking a look!
Best,
Megan

[jveitchmichaelis]

Hey Megan,

Here's a working fit (excuse the zip - I can't upload a notebook here):

spectroscopy_w_rascal.tar.gz

The main adjustments are:

• Filter the line list to your min/max range if you're using user lines, this should probably be fixed to respect min/max wavelengths
• Use the known lines as a prior for the detector range and use that as guesses for the fit min/max (not the same as the atlas), turns out the fit is extremely good and we can use a small range tolerance
• Use a higher threshold for peak finding, if you go too low you'll just try and fit noise. It makes sense to fit based on stronger peaks and then you can always go back and tweak. I used 10.
• Use much smaller values for the hough transform params - the ones in the notebook were very high and caused extremely slow fit initialisation. The settings in the examples should work fine (e.g. nbins = 200).

The moral is basically if you have a good instrument like this, you can afford to give pretty tight fitting parameters and rascal will run a lot better.

# Guesstimate range from known pairs:
x = [103, 165, 990, 1262, 1287]
y=[4921.889, 5015.6343, 6266.495, 6678.276, 6717.043]
x2, x1, x0 = np.polyfit(x, y, deg=2)

pix2wav = lambda p: x0 + x1*p + x2*p*p
print(pix2wav(0))
print(pix2wav(len(spectrum)))

rms = [_y-pix2wav(_x) for (_x, _y) in zip(x, y)]
print(rms)
# No need to give Rascal a very wide range if we know that the min and max wavelength visible are [4900-6720].
# This is not the range of the camera though, so we need to calculate the rough ends.


4765.636480489265
6848.529101379088
[0.11190128288217238, -0.13039686639331194, 0.08590426098180615, -0.46797601360049157, 0.4005673361452864]

atlas = Atlas()

waves = np.array(linelist['wave'])
lines = np.array(linelist['line'])
mask = (waves > 4900) * (waves < 6850)

# Note this doesn't filter by range, so we do it above (it's assumed that if you pass
# a user list, you're sure the lines are good.)
atlas.add_user_atlas(elements=lines[mask],
                     wavelengths=waves[mask],
                     pressure=80000.0,
                     temperature=285.0,)

c.set_atlas(atlas)

# The original values here were extremely high, so the fit was very slow
# no need to set them higher than this (at least for this case)
# Since we know very well what the params of the spectrometer are, we can
# set the range tolerance to be quite small and the min/max range of the
# detector to be based on a prior from the known lines:
c.set_hough_properties(num_slopes=2000,
                       xbins=100,
                       ybins=100,
                       min_wavelength=4700.,
                       max_wavelength=6850.,
                       range_tolerance=100.,
                       linearity_tolerance=50)

c.set_ransac_properties(sample_size=5,
                        top_n_candidate=8)

c.do_hough_transform()

res = c.fit()

{'fit_coeff': array([ 4.76868364e+03,  1.47624662e+00,  1.03868046e-04, -9.63339994e-08,
         2.97181764e-11]),
 'matched_peaks': [103.22040451624993,
  164.88369873051113,
  182.81148021279787,
  206.796932007569,
  231.4945354929042,
  280.6978560168899,
  376.69140210736606,
  418.44204522817375,
  525.0062835313796,
  587.7932569466659,
  610.2958267051915,
  629.3283765419278,
  658.7847550006928,
  696.2197128830677,
  716.496563150403,
  735.9203008572999,
  777.43588151281,
  797.5007498988053,
  833.5482900049134,
  862.909824383491,
  877.3178886179553,
  908.2883778170061,
  921.7795682128838,
  957.1513742167108,
  989.6916861294919,
  1014.9468041874037,
  1034.5028502825865,
  1066.574478690382,
  1079.2561712002414,
  1148.1889907224484,
  1165.5460003238143,
  1209.1631087449775,
  1261.6928041302358,
  1287.1993876547626],
 'matched_atlas': [4921.931,
  5015.6779,
  5037.751,
  5080.383,
  5116.503,
  5188.612,
  5330.7775,
  5400.5617,
  5562.766,
  5656.659,
  5689.8163,
  5719.225,
  5764.4188,
  5820.1558,
  5852.4878,
  5881.895,
  5944.8342,
  5975.534,
  6029.997,
  6074.3377,
  6096.1631,
  6143.0626,
  6163.5939,
  6217.2812,
  6266.495,
  6304.789,
  6334.4278,
  6382.9917,
  6402.246,
  6506.5281,
  6532.8822,
  6598.9529,
  6678.276,
  6717.043],
 'rms': 1.2715412885062707,
 'residual': array([ 0.13549909, -1.17130016,  3.72337615, -2.77177064, -1.61982266,
         0.68871456,  4.18325099, -2.11472924, -2.09685253, -0.37737011,
         0.72619036,  0.29039485, -0.07334182,  1.13927862, -0.35810561,
        -0.23697891, -0.09480144, -0.32274432, -0.06822963,  0.13451326,
         0.15904267,  0.20924158,  0.12288241,  0.02054029,  0.08145666,
         0.01930509, -0.02021648, -0.04823012, -0.11176909, -0.07140436,
        -0.1498031 , -0.16358807,  0.14694353,  0.09042893]),
 'peak_utilisation': 0.7391304347826086,
 'atlas_utilisation': 0.8095238095238095,
 'success': True}

Fit looks good, the unmatched peaks look either non-existent or are in the noise.

Hope that helps!

[jveitchmichaelis]

@cylammarco still a bug in the residual plot?

[cylammarco]

I think I fixed it in the new version not the current one.

[jveitchmichaelis]

That's true, @megan-the-astronomer you should try the notebook with the new_ransac branch. I think that's what I testd with.

[sugar-joh]

Hi, I have the same issue here and I am using the 0.3.8 version. I only know the wavelength is between ~5000 and 9500. Can you please advise?
data.txt

[jveitchmichaelis]

Hi @sugar-joh,

If you could provide more information that'd be a lot of help:

• What is your calibration source/line list? Can you provide that as well as your spectrum?
• Have you tried with the suggestions I made above?
  • Filter the potential lines to within the range that you're interested in (i.e. no lines outside 5000-9500)
• Please provide a code example with how you're setting up the solver.

I'd suggest trying out the new_ransac branch and modfying the example I provided above to your system.

We can support more easily if we know what you've tried!

Best
Josh

[sugar-joh]

Hi, I have try to use the new branch but it is not working. I attached the calibration line list and the spectrum as well.
Here is the code, I barely change anything from the tutorial except the min/max wavelength.

peaks, _ = find_peaks(spectrum.flux, prominence=200)
peaks_refined = refine_peaks(spectrum=spectrum.flux, peaks=peaks)
c = Calibrator(peaks_refined, spectrum.flux)
c.set_calibrator_properties(num_pix=len(spectrum.flux),
plotting_library='matplotlib',
log_level='info')
c.plot_arc()
atlas = Atlas()
atlas.add_user_atlas(elements=element,
wavelengths=arc_wav,
pressure=61700.,
temperature=276.55)
c.set_atlas(atlas)
c.set_hough_properties(num_slopes=2000,
xbins=50,
ybins=50,
min_wavelength=5500.,
max_wavelength=9500.,
range_tolerance=100.,
linearity_tolerance=50)
c.set_ransac_properties(sample_size=5,
top_n_candidate=8)
c.do_hough_transform()

arc_wav.txt
data.txt