take_iv fails with SystemError: <class 'numpy.linalg.LinAlgError'> returned a result with an error set

[jennamo]

filepaths from log:
[ 2023-08-14 14:52:01 ] Writing to file : /data/smurf_data/20230814/ufm_mv9/1692020663/outputs/1692024721.dat
[ 2023-08-14 14:52:01 ] /data/smurf_data/20230814/ufm_mv9/1692020663/outputs/1692024721_mask.txt
full traceback:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
ValueError: On entry to DLASCL parameter number 4 had an illegal value

The above exception was the direct cause of the following exception:

SystemError                               Traceback (most recent call last)
Cell In [48], line 1
----> 1 iva= iv.take_iv(S, cfg)

File /usr/local/src/pysmurf/python/pysmurf/client/util/pub.py:50, in set_action.<locals>.dec.<locals>.wrapper(S, pub_action, *args, **kwargs)
     46             S.pub._action = func.__name__
     48         S.pub._action_ts = S.get_timestamp()
---> 50     rv = func(S, *args, **kwargs)
     52 finally:
     53     if is_top:

File /sodetlib/sodetlib/operations/iv.py:647, in take_iv(S, cfg, bias_groups, overbias_voltage, overbias_wait, high_current_mode, cool_wait, cool_voltage, biases, bias_high, bias_low, bias_step, show_plots, wait_time, run_analysis, run_serially, serial_wait_time, g3_tag, **analysis_kwargs)
    645 _analysis_kwargs = {'save': True, 'update_cfg': True}
    646 _analysis_kwargs.update(analysis_kwargs)
--> 647 analyze_iv(iva, **_analysis_kwargs)
    649 # Save and publish plots
    650 is_interactive = plt.isinteractive()

File /sodetlib/sodetlib/operations/iv.py:375, in analyze_iv(iva, psat_level, save, update_cfg, show_pb)
    371 iva.resp[i] -= norm_fit[1]  # Put resp in real current units
    373 fin_idx = (np.isfinite(iva.i_bias[:sc_idx]) &
    374            np.isfinite(iva.resp[i, :sc_idx]))
--> 375 sc_fit = np.polyfit(iva.i_bias[:sc_idx][fin_idx],
    376                     iva.resp[i, :sc_idx][fin_idx], 1)
    378 # subtract off unphysical y-offset in superconducting branch; this
    379 # is probably due to an undetected phase wrap at the kink between
    380 # the superconducting branch and the transition, so it is
    381 # *probably* legitimate to remove it by hand. 
    382 iva.resp[i, :sc_idx] -= sc_fit[1]

File <__array_function__ internals>:180, in polyfit(*args, **kwargs)

File /usr/local/lib/python3.8/dist-packages/numpy/lib/polynomial.py:668, in polyfit(x, y, deg, rcond, full, w, cov)
    666 scale = NX.sqrt((lhs*lhs).sum(axis=0))
    667 lhs /= scale
--> 668 c, resids, rank, s = lstsq(lhs, rhs, rcond)
    669 c = (c.T/scale).T  # broadcast scale coefficients
    671 # warn on rank reduction, which indicates an ill conditioned matrix

File <__array_function__ internals>:180, in lstsq(*args, **kwargs)

File /usr/local/lib/python3.8/dist-packages/numpy/linalg/linalg.py:2300, in lstsq(a, b, rcond)
   2297 if n_rhs == 0:
   2298     # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
   2299     b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
-> 2300 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
   2301 if m == 0:
   2302     x[...] = 0

File /usr/local/lib/python3.8/dist-packages/numpy/linalg/linalg.py:101, in _raise_linalgerror_lstsq(err, flag)
    100 def _raise_linalgerror_lstsq(err, flag):
--> 101     raise LinAlgError("SVD did not converge in Linear Least Squares")

SystemError: <class 'numpy.linalg.LinAlgError'> returned a result with an error set

[jennamo]

first saw this behavior on crate 1 slot 2 (Mv9), now also seeing on crate 1 slot 7 (Mv15)

[dpdutcher]

I just put in a PR #380

This just skips any channel that fails to polyfit, since that would be junk any way and I'd prefer the whole function to not fail on account of any one channel. I don't know the specifics for why this error occurred though, and this just sweeps the issue under the rug.

[jlashner]

I'm investigating what is causing this, and I think there is a bug in the IV analysis code.

This code is what is selecting the timestream region for a given bias in the IV:

sodetlib/sodetlib/operations/iv.py

Lines 324 to 333 in 5a579b5

	for i in trange(iva.nbiases, disable=(not show_pb)):
	# Start from back because analysis is easier low->high voltages
	for j, bg in enumerate(iva.bias_groups):
	t0, t1 = iva.start_times[bg, -(i+1)], iva.stop_times[bg, -(i+1)]
	chan_mask = iva.bgmap == bg
	m = (t0 < am.timestamps) & (am.timestamps < t1)
	iva.resp[chan_mask, i] = np.nanmean(am.signal[:, m][chan_mask], axis=1)
	if j == 0:
	bias_bits = np.median(am.biases[bg, m])
	iva.v_bias[i] = bias_bits * 2 * iva.meta['rtm_bit_to_volt']

where the start/stop times are set here:

sodetlib/sodetlib/operations/iv.py

Lines 624 to 627 in 5a579b5

	S.set_tes_bias_bipolar_array(_bias_arr)
	start_times[bgs, i] = time.time()
	time.sleep(wait_time)
	stop_times[bgs, i] = time.time()

Based on this, the bias level should never change in a given segment... however if we look at the bias in a selected region, we see that it does change, and it is far from the bias we expect.

The reason being.... if there is a timing system, am.timestamps will be the timestamps determined by the Counters, which are about a second off from the system time:

so we are selecting the wrong region of data for a given bias value.

If we select the region based on am.primary.UnixTime instead of am.timestamps, that is the same "time" units that are being used when we set start/stop times, so we get agreement:

So this issue was happening because for the last few bias points, it was trying to grab a region of the timestream that doesn't have any data, due to the mismatch in "time" units. This becomes a problem and throws this error if the sc_index is close to the end of the timestream, in which case it's trying to polyfit a bunch of nans.

I think this can be solved by just switching am.timestamps to am.primary.UnixTime in the analyze_iv function, which should make IVs more accurate in general (in the presence of a timing system).

[dpdutcher]

What's the behavior in the absence of a timing system?

[jlashner]

Without a timing system, am.timestamps defaults to UnixTime, so this issue doesn't exist

[jlashner]

Just wanted to comment here, but setting up NTP on smurf-servers (timing sever and smurf-srv20) did seem to fix the timestamp difference between Timing System time and UnixTime:

[msilvafe]

That’s good news. @BrianJKoopman found that telegraf can report metrics from ntpq which we can monitor and alarm on in grafana. Let’s make sure Smurf servers get updated when we add this to telegraf (in the vein of daq-discussion #6)