From 5360d3f82a7ed4c1b756c2311e7e8145e4ac6584 Mon Sep 17 00:00:00 2001 From: lauraporta <29216006+lauraporta@users.noreply.github.com> Date: Thu, 25 Apr 2024 18:28:44 +0100 Subject: [PATCH] Add methods to fit the models and scripts to find optimal weights given the fold changes --- pyproject.toml | 7 +- translocator_models/fit_the_model.py | 193 ++++++++++++++++++++++ translocator_models/fitting_methods.py | 214 +++++++++++++++++++++++++ 3 files changed, 413 insertions(+), 1 deletion(-) create mode 100644 translocator_models/fit_the_model.py create mode 100644 translocator_models/fitting_methods.py diff --git a/pyproject.toml b/pyproject.toml index 38324a4..10b20c0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,12 @@ readme = "README.md" requires-python = ">=3.9.0" dynamic = ["version"] -dependencies = [] +dependencies = [ + "pandas", + "numpy", + "scipy", + "seaborn" +] license = {text = "BSD-3-Clause"} diff --git a/translocator_models/fit_the_model.py b/translocator_models/fit_the_model.py new file mode 100644 index 0000000..50eed88 --- /dev/null +++ b/translocator_models/fit_the_model.py @@ -0,0 +1,193 @@ +import pandas as pd +import seaborn as sns +from fitting_methods import ( + fit_fold_changes_to_data, + get_predicted_fold_changes_matched_dataset, + get_predicted_fold_changes_passive_same_luminance, + get_predicted_fold_changes_visual_flow, +) +from models import arithmetic_sum_model + +# =============================================== +# Mean fold changes across datasets +# visual flow, 39 clusters +visual_flow = {} +visual_flow["RV"] = 2.52 +visual_flow["RVT"] = 2.42 +visual_flow["V"] = 1.88 +visual_flow["VT"] = 2.13 + +# passive same luminance, 151 clusters +passive_same_luminance = {} +passive_same_luminance["T"] = 1.44 +passive_same_luminance["V"] = 1.16 +passive_same_luminance["VT"] = 1.76 + +# Re-balanced fold changes, together in the same dataset +matched_dataset = {} +matched_dataset["T"] = 1.74 +matched_dataset["V"] = 1.88 +matched_dataset["VT"] = 2.13 +matched_dataset["RV"] = 2.52 +matched_dataset["RVT"] = 2.42 +matched_dataset["RV_slip"] = 2.61 +matched_dataset["RVT_slip"] = 3.045 + +# Weights of the arithmetic sum model, fitted on the passive same luminance +# dataset +b0 = 0.08 +b1 = 0.74 +b2 = 0.53 + +# =============================================== +# Fitting parameters to visual flow dataset +print("===============================================") +print("Visual flow dataset") + +best_result = fit_fold_changes_to_data( + data=[ + visual_flow["V"], + visual_flow["VT"], + visual_flow["RV"], + visual_flow["RVT"], + ] +) + +k, c, w1, w2, w3 = best_result.x + +print(f"fc(v) = {k:.2f} * v + {c:.2f}") +print(f"v = {w1:.2f} (VF > 0) + {w2:.2f} (T - R > 0)(T - R) + {w3:.2f} R") + +predicted_fc = get_predicted_fold_changes_visual_flow(best_result.x) + + +print(f"V: {predicted_fc[0]:.2f}, mean from data: {visual_flow['V']}") +print(f"VT: {predicted_fc[1]:.2f}, mean from data: {visual_flow['VT']}") +print(f"RV: {predicted_fc[2]:.2f}, mean from data: {visual_flow['RV']}") +print(f"RVT: {predicted_fc[3]:.2f}, mean from data: {visual_flow['RVT']}") + + +# =============================================== +# Fitting parameters to passive_same_luminance dataset +print("===============================================") +print("Passive same luminance dataset") + + +best_result = fit_fold_changes_to_data( + data=[ + passive_same_luminance["V"], + passive_same_luminance["VT"], + passive_same_luminance["T"], + ], + dataset="passive_same_luminance", +) + + +k, c, w1, w2, w3 = best_result.x + +print(f"fc(v) = {k:.2f} * v + {c:.2f}") +print(f"v = {w1:.2f} (VF > 0) + {w2:.2f} (T - R > 0)(T - R) + {w3:.2f} R") + +predicted_fc = get_predicted_fold_changes_passive_same_luminance(best_result.x) + + +print( + f"V: {predicted_fc[0]:.2f}, " + + "mean from data: {passive_same_luminance['V']}" +) +print( + f"VT: {predicted_fc[1]:.2f}, " + + "mean from data: {passive_same_luminance['VT']}" +) +print( + f"T: {predicted_fc[2]:.2f}, " + + "mean from data: {passive_same_luminance['T']}" +) + + +arithmetic_sum = {} +arithmetic_sum["VT"] = arithmetic_sum_model( + fc_T_VStatic=predicted_fc[2], + fc_VF=predicted_fc[0], + b0=b0, + b1=b1, + b2=b2, +) + +print( + "VT as predicted with the arithmetic sum + rate model: " + + f"{arithmetic_sum['VT']:.2f}" +) + + +# =============================================== +# matched dataset +print("===============================================") +print("Visual flow + mismatch matched dataset + T") + + +best_result = fit_fold_changes_to_data( + data=[ + matched_dataset["T"], + matched_dataset["V"], + matched_dataset["VT"], + matched_dataset["RV"], + matched_dataset["RVT"], + matched_dataset["RV_slip"], + matched_dataset["RVT_slip"], + ], + dataset="matched", +) + + +k, c, w1, w2, w3 = best_result.x + +print(f"fc(v) = {k:.2f} * v + {c:.2f}") +print(f"v = {w1:.2f} (VF > 0) + {w2:.2f} (T - R > 0)(T - R) + {w3:.2f} R") + +predicted_fc = get_predicted_fold_changes_matched_dataset(best_result.x) + +print(f"T: {predicted_fc[0]:.2f}, mean from data: {matched_dataset['T']}") +print(f"V: {predicted_fc[1]:.2f}, mean from data: {matched_dataset['V']}") +print(f"VT: {predicted_fc[2]:.2f}, mean from data: {matched_dataset['VT']}") +print(f"RV: {predicted_fc[3]:.2f}, mean from data: {matched_dataset['RV']}") +print(f"RVT: {predicted_fc[4]:.2f}, mean from data: {matched_dataset['RVT']}") +print( + f"RV_slip: {predicted_fc[5]:.2f}, " + + "mean from data: {matched_dataset['RV_slip']}" +) +print( + f"RVT_slip: {predicted_fc[6]:.2f}, " + + "mean from data: {matched_dataset['RVT_slip']}" +) + + +arithmetic_sum = {} +arithmetic_sum["VT"] = arithmetic_sum_model( + fc_T_VStatic=predicted_fc[0], + fc_VF=predicted_fc[1], + b0=b0, + b1=b1, + b2=b2, +) + +print( + "VT as predicted with the arithmetic sum + rate model: " + + f"{arithmetic_sum['VT']:.2f}" +) + +df = pd.DataFrame( + { + "trial_type": matched_dataset.keys(), + "data": matched_dataset.values(), + "predictions": predicted_fc, + }, +) +df = pd.melt( + df, + id_vars="trial_type", + value_vars=["data", "predictions"], + var_name="group", + value_name="fold_change", +) +sns.barplot(df, x="trial_type", y="fold_change", hue="group") diff --git a/translocator_models/fitting_methods.py b/translocator_models/fitting_methods.py new file mode 100644 index 0000000..1e4ee76 --- /dev/null +++ b/translocator_models/fitting_methods.py @@ -0,0 +1,214 @@ +import numpy as np +from models import rate_based_model +from scipy.optimize import minimize + + +def single_fit(params, data, dataset="visual_flow"): + def objective_function(params): + if dataset == "visual_flow": + predicted_fc = get_predicted_fold_changes_visual_flow(params) + elif dataset == "passive_same_luminance": + predicted_fc = get_predicted_fold_changes_passive_same_luminance( + params + ) + elif dataset == "matched": + predicted_fc = get_predicted_fold_changes_matched_dataset(params) + residuals = data - predicted_fc + return np.sum(np.abs(residuals.ravel())) + + bounds = ( + (0, 3), # k + (1, 1), # c = 1 in our definition + (1, 1), # w1 = 1 in our definition + (0, 3), # w2 + (0, 3), # w3 + ) + result = minimize( + objective_function, + params, + bounds=bounds, + method="Nelder-Mead", + ) + return result + + +def get_predicted_fold_changes_visual_flow(params): + k, c, w1, w2, w3 = params + + V = rate_based_model( + VF=1, + T=0, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + VT = rate_based_model( + VF=1, + T=1, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RV = rate_based_model( + VF=1, + T=0, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RVT = rate_based_model( + VF=1, + T=1, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + return np.asarray([V, VT, RV, RVT]) + + +def get_predicted_fold_changes_passive_same_luminance(params): + k, c, w1, w2, w3 = params + + V = rate_based_model( + VF=1, + T=0, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + VT = rate_based_model( + VF=1, + T=1, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + T = rate_based_model( + VF=0, + T=1, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + return np.asarray([V, VT, T]) + + +def get_predicted_fold_changes_matched_dataset(params): + k, c, w1, w2, w3 = params + + V = rate_based_model( + VF=1, + T=0, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + VT = rate_based_model( + VF=1, + T=1, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RV = rate_based_model( + VF=1, + T=0, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RVT = rate_based_model( + VF=1, + T=1, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RVT_slip = rate_based_model( + VF=2, + T=2, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + RV_slip = rate_based_model( + VF=2, + T=0, + R=1, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + T = rate_based_model( + VF=0, + T=1, + R=0, + k=k, + c=c, + w1=w1, + w2=w2, + w3=w3, + ) + + return np.asarray([T, V, VT, RV, RVT, RV_slip, RVT_slip]) + + +def fit_fold_changes_to_data( + data, + dataset="visual_flow", + initial_params=np.asarray([0.85, 1, 1, 0.6, 1]), +): + result = single_fit(initial_params, data, dataset) + + return result