RayleighFit crash

[HavardStridBuholdt]

The RayleighFit module is very sensitive to its input and the config values and have been experienced to crash at times, resulting in this error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[40], [line 1](vscode-notebook-cell:?execution_count=40&line=1)
----> [1](vscode-notebook-cell:?execution_count=40&line=1) data_cube.rayleighFit()
      2 # Provides the reference heights

File c:\Users\buholdt\Documents\PicassoPy\tests\..\ppcpy\interface\picassoProc.py:433, in PicassoProc.rayleighFit(self)
    430 print('Start Rayleigh Fit')
    431 logging.warning(f'Potential for differences to matlab code du to numerical issues (subtraction of two small values)')
--> [433](file:///C:/Users/buholdt/Documents/PicassoPy/ppcpy/interface/picassoProc.py:433) self.refH =  rayleighfit.rayleighfit(self)
    434 return self.refH

File c:\Users\buholdt\Documents\PicassoPy\tests\..\ppcpy\calibration\rayleighfit.py:69, in rayleighfit(data_cube)
     56 DPInd = DouglasPeucker(
     57     scaRatio, height, 
     58     config_dict[f'minDecomLogDist{wv}'], np.array(config_dict['heightFullOverlap'])[data_cube.gf(wv, t, tel)], 
     59     config_dict[f'maxDecomHeight{wv}'],
     60     config_dict[f'maxDecomThickness{wv}'], config_dict[f'decomSmoothWin{wv}']
     61     )
     63 #print('DPInd', DPInd)
     64 #DPInd = [133,  332,  333,  533,  675,  875,  1075,  1274,  1275,  1333]
     65 #print('!!!! overwrite DPInd !!!!!!!!!!!!!!!!!!!!')
     66 #print('DPInd', DPInd)
     67 
     68 #p.Results.minRefThickness, p.Results.minRefDeltaExt, p.Results.minRefSNR, flagShowDetail
---> [69](file:///C:/Users/buholdt/Documents/PicassoPy/ppcpy/calibration/rayleighfit.py:69) refHInd = fit_profile(
     70     height, rcs, sig, bg, mSig, DPInd,
     71     config_dict[f'minRefThickness{wv}'], config_dict[f'minRefDeltaExt{wv}'],
     72     config_dict[f'minRefSNR{wv}'], flagShowDetail=False)
     73 print('refHInd', refHInd)
     75 # for debugging reasons
     76 #return rcs, mSig, scaRatio
     77 
     78 # calculate PCR after averaging or average directly PCR -> PCR from average above was cross checked
     79 # with matlab. Above 15km, the relative difference increased to ~ 3%
     80 # continue at https://github.com/PollyNET/Pollynet_Processing_Chain/blob/b3b8ec7726b75d9db6287dcba29459587ca34491/lib/interface/picassoProcV3.m#L788

File c:\Users\buholdt\Documents\PicassoPy\tests\..\ppcpy\calibration\rayleighfit.py:621, in fit_profile(height, sig_aer, pc, bg, sig_mol, dpIndx, layerThickConstrain, slopeConstrain, SNRConstrain, flagShowDetail)
    618 X_val = (np.abs(mean_resid) * np.abs(std_resid) * np.abs(slope_resid) * 
    619          np.abs(msre_resid) * np.abs(Astat) / SNR_ref)
    620 X_val[X_val == 0] = np.nan
--> [621](file:///C:/Users/buholdt/Documents/PicassoPy/ppcpy/calibration/rayleighfit.py:621) indxBest_Int = np.nanargmin(X_val)
    623 # print('hIndB_Test', hIndxB_Test)
    624 # print('hIndT_Test', hIndxT_Test)
    625 # print('mean_resid', mean_resid)
   (...)    630 # print('Astat', Astat)
    631 # print('X_val', X_val)
    633 if flagShowDetail:

File c:\Users\buholdt\AppData\Local\miniconda3\envs\PicassoPy\Lib\site-packages\numpy\lib\_nanfunctions_impl.py:562, in nanargmin(a, axis, out, keepdims)
    560     mask = np.all(mask, axis=axis)
    561     if np.any(mask):
--> [562](file:///C:/Users/buholdt/AppData/Local/miniconda3/envs/PicassoPy/Lib/site-packages/numpy/lib/_nanfunctions_impl.py:562)         raise ValueError("All-NaN slice encountered")
    563 res = np.argmin(a, axis=axis, out=out, keepdims=keepdims)
    564 return res

ValueError: All-NaN slice encountered

The crash has so far only been observed for case: tjk_20250823, where the 1064 total FR channel passes all the tests for one cloud free period but fails to produce a reference height region. Depending on the inputs/config values, this crash may occur for different cloud free periods, usually during daytime and always for the 1064 total FR channel. The tests fails for all of the DP index pairs except one for the issue channel and cloud free period. For the one pair of DP indices it does not fail X_val = 0 (due to msre_resid = 0), and is thus changed to NaN. Which means that we have a full NaN array which we try to take the np.nanmin() and this causes the error shown above. See the following code snippet:

# search the best fit region
X_val = (np.abs(mean_resid) * np.abs(std_resid) * np.abs(slope_resid) * 
             np.abs(msre_resid) * np.abs(Astat) / SNR_ref)
X_val[X_val == 0] = np.nan
indxBest_Int = np.nanargmin(X_val)

Origin: https://github.com/PollyNET/PicassoPy/blob/509f6076e3bd98e5acfbf2c10cf0e621f42c8e7d/ppcpy/calibration/rayleighfit.py#L530C4-L534C39

This crash will only happen in Python as the Matlab equivalent code will return indxBest_Int = 1. This is most likely a bug in the Matlab version. See the following code snippet:

% search the best fit region
X_val = abs(mean_resid) .* abs(std_resid) .* abs(slope_resid) .* ...
        abs(msre_resid) .* abs(Astat) ./ SNR_ref;
% X_val = abs(slope_resid) .* abs(msre_resid) .* abs(Astat);
X_val(X_val == 0) = NaN;
[~, indxBest_Int] = min(X_val);

Origin: https://github.com/PollyNET/Pollynet_Processing_Chain/blob/e6aba147286eaa97480df7b2cfe79eda287ac9d6/lib/calibration/rayleighfit.m#L260C1-L266C1

---

Additionally, the translation of the DouglasPeucker() algorithm to Python looks strange. I fear that parts of this code may be hardcoded to Matlabs 1-based indexing.

[HavardStridBuholdt]

The crash is also highly sensitive to the background.
For example:
Taking the background from range bin [10-240] will lead to the module crashing. However, taking the module form [9-240] will not cause the module to crash. A similar sensitivity is also seen for the first range bin index config values (250 vs 252) in versions of main before merged pull request #43.

The fact that the crash only is observed for the 1064 total FR channel may hint to that the rayleighFit config values for this channel is not optimal, which may cause this unstable behaver.

Another possibility for this unstable behaver is that there seems to be an undetected cloud in the scene at the cloud free period that most commonly crashes --> 06:40 to 07:40 UTC

[martin-rdz]

Starting to collect my thoughts.

The translation of the DP method from matlab to python relied strongly on AI, so I would be not suprised if there are still indexing errrors. I considered using an existing python implemenation of DP, such as this one, but then prefered to have it self-contained.

Also you are right about the missed cloud. Maybe for now we should do the test at a different time period and wait for improvements to the cloud screening later on.

[HavardStridBuholdt]

Following are some further updates to the issue:

A quickfix for the crash was added in Pull request #45. Where now a test is preformed to check if all values of X_val are NaN. If so the reference height indices are set to [np.nan, np.nan]. See the following code snippet:

# search the best fit region
X_val = (np.abs(mean_resid) * np.abs(std_resid) * np.abs(slope_resid) * 
         np.abs(msre_resid) * np.abs(Astat) / SNR_ref)
X_val[X_val == 0] = np.nan

# QuicFix for all NaN arrays:
if np.isnan(X_val).all():
    print("None valid clean region found.")
    return np.nan, np.nan
indxBest_Int = np.nanargmin(X_val)

Origin: https://github.com/PollyNET/PicassoPy/blob/a39d0104551e7a2191fcf1f8d0256a7291c54dec/ppcpy/calibration/rayleighfit.py#L596C5-L605C39

This differs form how it is done in Picasso (Matlab) where the indices connected to the first value of X_val would be used in these scenarios.

Though the implemented fix is successful in addressing the crash issue, it does not tackle the root of the problem, ie. why the rayleigFit procedure is so unstable in this case for the 1064nm channel. A separate branch 48-rayleighfit-crash has been made for the purpose of investigating this further issue further where the original crash is recreated.