ENH: Add GP error analysis experiment script

[jhlegarreta]

Add GP error analysis experiment script.

[jhlegarreta]

A few notes:

• Pardon my not very elaborated commit message.
• Started placing the methods in new modules, but I do not think that is useful: even if e.g. a crossing fiber simulation script would benefit from some methods, these methods are not required for the framework to work/not necessary for the package, they would require adding tests, etc. So unnecessary work. So I would sacrifice that for having duplicated code across script in this case. At least for now.
• The script can be launched as e.g.

pyhon /path/to/eddymotion/scripts/dwi_estimation_error_analysis.py \ 
  60 1000 100 \
  --evals1 0.0015 0.0003 0.0003 \
  --snr 20 \
  --repeats 3 \
  --kfold 1 3 5 15 30

The result of the above command should be:

Unclear why the error is not monotonically increasing.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 66.18%. Comparing base (a61c5e1) to head (edc0331).
Report is 18 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #227   +/-   ##
=======================================
  Coverage   66.18%   66.18%           
=======================================
  Files          19       19           
  Lines         905      905           
  Branches      112      112           
=======================================
  Hits          599      599           
  Misses        264      264           
  Partials       42       42           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jhlegarreta]

Using the std deviation value of scikit-learn's GP regressor class (committed in
ea93c4a for the records, but most likely will drop when we get this sorted out) we get the below results, which is again not what we would expect

[jhlegarreta]

@oesteban Adopted scikit learn's KFold to create folds in acf683b. One would call the script as e.g.

pyhon /path/to/eddymotion/scripts/dwi_estimation_error_analysis.py \
  60 1000 100 --evals1 0.0015 0.0003 0.0003 --repeats 3 --snr 20 --kfold 3 5 7 10 20 120

• Kept the repeats parameter but not using them now.
• The RMSE increases as SNR decreases, which seems OK (SNR None means no noise, so the error is smaller than with SNR 20, whose error is smaller than with SNR 10).
• The trend of each output is not what I would expect: error and std dev do not decrease as we add more volumes to the training set and predict less volumes each time (i.e. at increasing values of N; at 120 we are training with 119 volumes and predicting 1).

[oesteban]

@jhlegarreta I've moved things to use scikit-learn's cross-validation to avoid errors on our side implementing the CV framework:

1. Updated the model to take arrays instead of gradient tables (ff08c26). We probably will want to add a checkpoint in fit and check if y is a gradient table (and extract bvecs accordingly). But ideally, we want the model to self adapt to incoming ndarrays or DIPY's gradient tables.
2. Then, I use the default k-fold CV (b963b96).

Running this:

python scripts/dwi_estimation_error_analysis.py 60 1000 100 --evals1 0.0015 0.0003 0.0003 --repeats 3 --snr 20 --kfold 5 10 60 --repeats 200

I get

{5: (-6.5451513798923475, 0.7829974075194627), 10: (-6.704245824324263, 1.451348158720791), 60: (-6.339094677419443, 2.883843582961101)}

which I believe is consistent with the expectations:

• RMSE is negative because scikit-learn negates it.
• Error is constant around 6. something, however the spread grows with larger folds. This makes sense to me because the more folds the more fits and the bigger the probability of getting RMSEs further from the average. Also the average is with only two values, so it's also easier to get larger errors. With repeats=200 then k=5 means 1000 fits while k=60 means 12000 fits. I'll try to compensate for this later.

Regarding the number of directions in training, the interpretation of the k-fold is also different from what I think we were discussing:

k = 5 -> 5 folds, 120/5=24 bvecs per test fold -> 96 bvecs for training.
k = 60 -> 60 folds, 120/60 = 2 bvecs per test fold -> 118 bvecs for training.

I keep receiving the following warnings though:

/data/home/oesteban/workspace/scikit-learn/sklearn/gaussian_process/kernels.py:442: ConvergenceWarning: The optimal value found for dimension 0 of parameter a is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.
/data/home/oesteban/workspace/scikit-learn/sklearn/gaussian_process/kernels.py:452: ConvergenceWarning: The optimal value found for dimension 0 of parameter lambda_s is close to the specified upper bound 10000.0. Increasing the bound and calling fit again may find a better value.
/data/home/oesteban/workspace/scikit-learn/sklearn/gaussian_process/kernels.py:442: ConvergenceWarning: The optimal value found for dimension 0 of parameter sigma_sq is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.

and:

/data/home/oesteban/workspace/scikit-learn/sklearn/gaussian_process/_gpr.py:659: ConvergenceWarning: lbfgs failed to converge (status=2):
ABNORMAL_TERMINATION_IN_LNSRCH.

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
  _check_optimize_result("lbfgs", opt_res)

WDYT?

[oesteban]

Decreasing the SNR did not help with the warnings. My hypothesis being that noisier data would regularize the problem.

[oesteban]

Okay, decreasing the SNR had the expected effect on the cross-validation:

python scripts/dwi_estimation_error_analysis.py 60 1000 100 --evals1 0.0015 0.0003 0.0003 --repeats 3 --snr 10 --kfold 5 10 60 --repeats 200

Results in

{5: (-13.634124763730938, 1.0844600539135052), 10: (-13.739925211912865, 1.6677184190441654), 60: (-13.892240109267137, 4.902059389908273)}

[jhlegarreta]

@jhlegarreta I've moved things to use scikit-learn's cross-validation to avoid errors on our side implementing the CV framework:
(...)
which I believe is consistent with the expectations:

Looks like it does at least concerning the folds: the larger the value of k, the smaller (although marginally) the mean of the error.

Regarding the number of directions in training, the interpretation of the k-fold is also different from what I think we were discussing:

Not sure 🤔. My expectation is that the more bvecs for training (the larger k), the smaller the error, and that is what we see with your numbers.

What I got wrong at the beginning was the intuition with the SNR, but I think the trend with the SNR shown in the plots is correct.

TBH, in the plots that I pasted, I have the impression that if we have a very close look at the SNR 10 case (and may be the SNR 20 as well), there is a very slight decrease at N (k) 120 wrt to 20 (have not checked quantitatively). So unless I am missing something, this looks consistent with your results. But I'd expect the error to be notably, and not marginally, smaller.

With repeats=200 then k=5 means 1000 fits while k=60 means 12000 fits. I'll try to compensate for this later.

Since my last commit/the adoption of proper k-folds, we are not using the repeats, so not sure how to assess this. Please ignore my comment if you believe that yours makes sense.

I keep receiving the following warnings though:

Yes, I am also getting the abnormal termination notice. Increasing the tolerance bounds did not result in any benefit for me, so there is something else that we are missing, either in a previous step or in the controls of scikit-learn's GP.

In DIPY we were getting warnings when using a given CVXPY optimizer and they went away when switching to another one. In that case the warning itself was suggesting to switch to another one; not sure if in this case we have that control.

Okay, decreasing the SNR had the expected effect on the cross-validation:

So the smaller the SNR, the larger the error, which is expected, and matches what we saw in the plots if I read well ❔.

[jhlegarreta]

31e5972 adds the capability to generate a carpet plot with the ground truth and the estimated signal (see an example in nipreps/nireports#137) and a correlation plot, which looks like this for the noiseless and k=5 case:

The calls are commented since after pulling the most recent commits I've seen that we no longer have the data available; had done this prior to your latest commits.

[oesteban]

@jhlegarreta I've found several issues - see 10a576d

I'll keep digging the kernel.theta setting. Keep you posted.

[oesteban]

@jhlegarreta Progress scikit-learn/scikit-learn#30098

EDITED: it was a false alarm.

[oesteban]

(the optimizer issue is still standing though)

[jhlegarreta]

Thanks for all this investigation.

Let me know if you are able to push the latest changes so that I can continue investigating on my end.

Two notes:

• In Andersson 2015 they reparameterize the model to estimate lambda and sigma in order to avoid negative values (p.168).
• If the gradient vector covariance lead turns out to be the good one towards fixing things with the GP, we will need to add one step to get the signal predictions: eq. 7 in Andersson 2015 (p. 167).

---

Unrelated note: eventually, after we manage to fix the current issues and if it makes sense to bring the implementation closer to DIPY's reconstruction models, we may want to use DIPY's cross validation:
https://github.com/dipy/dipy/blob/1961a91e2130aa9e2891a5610a79e2c9aca24d88/dipy/reconst/cross_validation.py#L58
https://github.com/dipy/dipy/blob/1961a91e2130aa9e2891a5610a79e2c9aca24d88/doc/examples/kfold_xval.py

[oesteban]

Currently, I have:

python scripts/dwi_estimation_error_analysis.py 60 1000 100 --evals1 0.0015 0.0003 0.0003 --snr 20 --kfold 3 5 10 15 20 25 30 40 50 60 --repeats 1000
{3: (-5.323663473768337, 0.5078300442425429), 5: (-5.592537291455177, 0.6344872918456032), 10: (-6.457823460261486, 1.1560694555678308), 15: (-6.0008898214350985, 1.2415790700142275), 20: (-4.785889809484692, 1.2121841573793941), 25: (-5.02350508651579, 1.7419826969482066), 30: (-4.902673405605454, 1.6647001726929596), 40: (-4.886979705096015, 1.9020697112634937), 50: (-5.276554893079895, 2.636081796110276), 60: (-5.393546086268717, 2.5409755115661166)}

[oesteban]

Running a comparison right now - let's keep moving :)