Skip to content

Conversation

@jhlegarreta
Copy link
Contributor

Ractor PET model: use a base class that contains the essential properties for a PET model and create a derived class that implements the B-Spline correction.

@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch from 2fb555e to 90eb886 Compare August 12, 2025 23:28
@jhlegarreta jhlegarreta requested a review from mnoergaard August 12, 2025 23:38
@jhlegarreta
Copy link
Contributor Author

jhlegarreta commented Aug 12, 2025

Had another look at this. I had introduced a BasePETModel based on the BaseDWIModel. Debugging the thing I currently do not see scenarios where we will need to have the _fit and predict_fit implemented in the BasePETModel, but last time we did discuss the possibility of adding other models. However, as I see things, we wouldn't need to have have _fit and predict_fit implemented in the base class for PET. So, unless you tell me the contrary and you clearly see that it may be beneficial and that you see this immediately @mnoergaard @oesteban, I would leave them as abstract functions and make all the current logic dwell in the BSplinePETModel.

)

# ToDo
# Are timepoints and xlim features that all PET models require ??
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mnoergaard question for you.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would say so, since we are dealing with motion correction, and hence would need temporal information. Currently, the PET model needs both the sampling midpoints and the scan’s total duration to work, so enforcing their presence in the base class keeps the API consistent. When new PET models that do not need temporal information are introduced, we can revisit this requirement; until then, providing dataset.midframe and dataset.total_duration (as done in the unit tests) is the intended usage.

Copy link
Contributor Author

@jhlegarreta jhlegarreta Sep 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I see that this is related to #204 (review). I think the naming should then be made consistent across the PET data class and the model class. Also, I do not see why the model is given the timepoints, if these are hold by the PET data class, i.e. from #204 (review)

(...) midframe is where the dataset stores real-world frame timing; timepoints/_x is the copy of those timings handed to the model

then we should not be providing the model with a copy of them.

If we make the parallel with the DWI class, the gradients are obtained from the dataset, e.g.:

gradient = self._dataset.gradients[:, index]

Also, can the xlim name be made somehow more descriptive or is it a name that is commonly used within the PET domain?

raise NotImplementedError("Fitting with held-out data is not supported")

# ToDo
# Does not make sense to make timepoints be a kwarg if it is provided as a named parameter to __init__
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mnoergaard question for you.


# ToDo
# What is the gtab equivalent of PET ?
model_str = getattr(self, "_model_class", "")
Copy link
Contributor Author

@jhlegarreta jhlegarreta Aug 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This may not be necessary.


# ToDo
# What are the gtab (and S0 if any) equivalent of PET ?
if n_models == 1:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.


# ToDo
# Does not make sense to make timepoints be a kwarg if it is provided as a named parameter to __init__
timepoints = kwargs.get("timepoints", None) or self._x
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

# Does the below apply to PET ? Martin has the return None statement
# if index is None:
# raise RuntimeError(
# f"Model {self.__class__.__name__} does not allow locking.")
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another question for @mnoergaard.

@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch from 90eb886 to afb11da Compare August 13, 2025 12:56
@jhlegarreta
Copy link
Contributor Author

The file test_model_pet.py will eventually need to be integrated into test_model.py as all models are being tested there. Unless we decide to put the dMRI models into a different file well @oesteban ?

@jhlegarreta
Copy link
Contributor Author

jhlegarreta commented Aug 18, 2025

Re #204 (comment) Compared to afb11da, in cbf417d I am discarding the idea of keeping a parallel of the dMRI base model in the PET base model. The PET model tests pass and the notebook runs now. It gives a result that is almost the same as in #203 (comment), with only minor differences in X and Y rotation start/ends. However, the parallelization was lost in the refactoring and the notebook takes almost 2 hours now compared to 30 minutes. So I need to look into that.

@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch from cbf417d to 10af03f Compare August 18, 2025 01:52
@codecov
Copy link

codecov bot commented Aug 18, 2025

Codecov Report

❌ Patch coverage is 70.45455% with 13 lines in your changes missing coverage. Please review.
✅ Project coverage is 75.78%. Comparing base (4a96c3c) to head (7d22d9e).

Files with missing lines Patch % Lines
src/nifreeze/model/pet.py 69.04% 9 Missing and 4 partials ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #204      +/-   ##
==========================================
- Coverage   75.95%   75.78%   -0.17%     
==========================================
  Files          24       24              
  Lines        1439     1458      +19     
  Branches      166      168       +2     
==========================================
+ Hits         1093     1105      +12     
- Misses        275      280       +5     
- Partials       71       73       +2     

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

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch 5 times, most recently from 6318c50 to 908faa8 Compare August 18, 2025 23:38
@jhlegarreta
Copy link
Contributor Author

jhlegarreta commented Aug 18, 2025

Re #204 (comment), the throttle issue was due to this line:

n_jobs = n_jobs or 1

which I had borrowed from

n_jobs = n_jobs or 1

Now changed to

n_jobs = n_jobs or min(cpu_count() or 1, 8)

To match more closely

n_jobs = n_jobs or 1

And the speed is back. Results in plots, as I was seeing yesterday (cf. the above comment), almost the same with minor differences:

pet_estimated_motion_pet_model_refactor pet_fd_pet_model_refactor

@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch 3 times, most recently from 4751545 to 19d5956 Compare August 22, 2025 22:39
@mnoergaard mnoergaard marked this pull request as ready for review September 22, 2025 07:52
# Calculate index coordinates in the B-Spline grid
self._n_ctrl = n_ctrl or (len(timepoints) // 4) + 1
def _preproces_data(self) -> np.ndarray:
# ToDo
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the todo here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I recall correctly, the idea there is to preprocess only the data that is required; e.g. if we are using only a subset of the volumes, then we shouldn't preprocess the entire data, e.g.

data, _, gtab = self._dataset[idxmask]
# Select voxels within mask or just unravel 3D if no mask
data = data[brainmask, ...] if brainmask is not None else data.reshape(-1, data.shape[-1])

The data, _, gtab = self._dataset[idxmask] is commented out in here immediately after the ToDo to as a pointer to that idea.


# Convert data into V (voxels) x T (timepoints)
data = data.reshape((-1, data.shape[-1])) if brainmask is None else data[brainmask]
# ToDo
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Allowing _fit to override timepoints through kwargs duplicates information that was already validated and stored on self._x during initialization. Dropping that extra kwarg simplifies the API and avoids inconsistent inputs

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure I understand: maybe the reply is answering a question elsewhere?

data = data.reshape((-1, data.shape[-1])) if brainmask is None else data[brainmask]
# ToDo
# Does not make sense to make timepoints be a kwarg if it is provided as a named parameter to __init__
timepoints = kwargs.get("timepoints", None) or self._x
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should be removed

# ToDo
# Does not make sense to make timepoints be a kwarg if it is provided as a named parameter to __init__
timepoints = kwargs.get("timepoints", None) or self._x
x = np.asarray(timepoints, dtype="float32") / self._xlim * self._n_ctrl
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be replaced with

x = self._x / self._xlim * self._n_ctrl

Copy link
Contributor

@mnoergaard mnoergaard left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @jhlegarreta. Thanks for pushing this forward. I have added some comments. And just to be clear on the notation: midframe is where the dataset stores real-world frame timing; timepoints/_x is the copy of those timings handed to the model; _xlim is the acquisition length used to scale those timings; and x is just the normalized coordinate derived from timepoints when solving or evaluating the B-spline. You agree?

Refactor PET model: use a base class that contains the essential
properties for a PET model and create a derived class that implements
the B-Spline correction.

- Make the `x_lim` private attribute be an array when assigning the
  corresponding value at initialization to enable NumPy broadcasting.
- Remove the unnecessary explicit `float` casting around the number of
  control points: the `dtype="float32"` specifier creates a float array.
  Fixes:
  ```
  Expected type 'str | Buffer | SupportsFloat | SupportsIndex', got 'Literal[0] | None | {__eq__} | int' instead
  ```
  raised by the IDE.
- Use `np.asarray` to avoid extra copies when not necessary.
@jhlegarreta jhlegarreta force-pushed the ref/refactor-pet-model branch from 21026bb to 7d22d9e Compare September 28, 2025 01:20
@oesteban oesteban linked an issue Nov 21, 2025 that may be closed by this pull request
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Finalize PET model merge

2 participants