How to do multiple structures mesh-valued regression?#

NB: an alternative way to using a for loop on How to do mesh-valued regression?.

[1]:

import pyvista as pv
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

import polpo.preprocessing.dict as ppdict
import polpo.preprocessing.pd as ppd
from polpo.models import DictMeshes2Comps, ObjectRegressor
from polpo.preprocessing.learning import NestedDictsToXY
from polpo.preprocessing.load.pregnancy.jacobs import (
    MeshLoader,
    TabularDataLoader,
)
from polpo.preprocessing.mesh.registration import RigidAlignment

[KeOps] Warning : cuda was detected, but driver API could not be initialized. Switching to cpu only.

[2]:

STATIC_VIZ = True

if STATIC_VIZ:
    pv.set_jupyter_backend("static")

Loading meshes#

[3]:

subject_id = "01"

[4]:

mesh_loader = (
    MeshLoader(
        subject_subset=[subject_id],
        as_mesh=True,
    )
    + ppdict.ExtractUniqueKey()  # session, struct
    + ppdict.NestedDictSwapper()  # struct, session
)

prep_pipe = ppdict.DictMap(RigidAlignment(max_iterations=500))


pipe = mesh_loader + prep_pipe

meshes = pipe()

Loading tabular data#

[5]:

pipe = TabularDataLoader(subject_id=subject_id)

df = pipe()

INFO: Data has already been downloaded... using cached file ('/home/luisfpereira/.herbrain/data/maternal/maternal_brain_project_pilot/rawdata/28Baby_Hormones.csv').

Here, we filter the tabular data.

[6]:

session_selector = ppd.DfIsInFilter("stage", ["post"], negate=True)

predictor_selector = (
    session_selector + ppd.ColumnsSelector("gestWeek") + ppd.SeriesToDict()
)

[7]:

x_dict = predictor_selector(df)

Merge data#

We get the data in the proper format for fitting.

[8]:

dict_pipe = NestedDictsToXY()

# meshes_ : dict[list]
X, meshes_ = dict_pipe([x_dict, meshes])

Create and fit regressor#

[9]:

pca = PCA(n_components=4)

objs2y = DictMeshes2Comps(n_pipes=len(meshes), dim_reduction=pca)

[10]:

model = ObjectRegressor(LinearRegression(fit_intercept=True), objs2y=objs2y)

[11]:

model.fit(X, meshes_)

[11]:

ObjectRegressor(objs2y=AdapterPipeline(steps=[('step_0', BiDictToValuesList()),
                                              ('step_1',
                                               MapTransformer(par_steps=[AdapterPipeline(steps=[('step_0',
                                                                                                 BiMeshesToVertices()),
                                                                                                ('step_1',
                                                                                                 FunctionTransformer(func=<function stack at 0x7828b0bd0870>)),
                                                                                                ('step_2',
                                                                                                 BiFlattenButFirst()),
                                                                                                ('step_3',
                                                                                                 StandardScaler(with_std=False)),
                                                                                                ('step_4',
                                                                                                 PCA(n_components=4))]),
                                                                         Adapte...
                                                                                                 BiFlattenButFirst()),
                                                                                                ('step_3',
                                                                                                 StandardScaler(with_std=False)),
                                                                                                ('step_4',
                                                                                                 PCA(n_components=4))]),
                                                                         AdapterPipeline(steps=[('step_0',
                                                                                                 BiMeshesToVertices()),
                                                                                                ('step_1',
                                                                                                 FunctionTransformer(func=<function stack at 0x7828b0bd0870>)),
                                                                                                ('step_2',
                                                                                                 BiFlattenButFirst()),
                                                                                                ('step_3',
                                                                                                 StandardScaler(with_std=False)),
                                                                                                ('step_4',
                                                                                                 PCA(n_components=4))])])),
                                              ('step_2', BiHstack())]))
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Evaluate fit#

model.predict outputs meshes, but we know LinearRegression sees PCA components. We can evaluate r2_score by applying transform.

NB: these are values on the training data.

[12]:

meshes_pred = model.predict(X)

y_true = objs2y.transform(meshes_)
y_pred = objs2y.transform(meshes_pred)

scores = r2_score(y_true, y_pred, multioutput="raw_values")

dict(zip(meshes.keys(), scores.reshape(len(meshes), -1)))

[12]:

{'L_Caud': array([0.19577571, 0.43087607, 0.03419797, 0.01492182]),
 'L_Pall': array([0.35915275, 0.04855135, 0.24886381, 0.05712201]),
 'L_Amyg': array([0.00954852, 0.25933387, 0.3971549 , 0.00524618]),
 'L_Accu': array([3.66642485e-01, 9.99477298e-03, 5.68769931e-02, 1.89650868e-04]),
 'L_Hipp': array([0.04356304, 0.00644152, 0.24583455, 0.00145247]),
 'L_Puta': array([0.66647764, 0.01291934, 0.03550199, 0.02858547]),
 'L_Thal': array([0.01415544, 0.53246226, 0.01803246, 0.19446471]),
 'R_Caud': array([0.28438218, 0.25672153, 0.06444699, 0.03906277]),
 'R_Pall': array([0.52617965, 0.02088388, 0.06184909, 0.07608363]),
 'R_Amyg': array([0.31313143, 0.10727699, 0.01111522, 0.00072599]),
 'R_Accu': array([0.26589627, 0.01162288, 0.01511771, 0.20374866]),
 'R_Hipp': array([0.03846713, 0.0608323 , 0.25813036, 0.0021322 ]),
 'R_Puta': array([0.31350877, 0.13438504, 0.12970855, 0.13376527]),
 'R_Thal': array([0.29409485, 0.19884905, 0.00710386, 0.04048809]),
 'BrStem': array([0.43173304, 0.00093647, 0.12750402, 0.02389235])}