Skip to content

Commit

Permalink
Fix bad initial value for Yamamoto method
Browse files Browse the repository at this point in the history
Change initial value from the automatic middle of the boundary space (500 & 500) to (1, 1), as an initial value of 500 can lead to numerical errors.

Also adapted yamamoto.py to PEP8
  • Loading branch information
r.jaepel committed Dec 1, 2023
1 parent bdb1fe0 commit 705aef8
Showing 1 changed file with 9 additions and 11 deletions.
20 changes: 9 additions & 11 deletions CADETProcess/tools/yamamoto.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
from CADETProcess import CADETProcessError
from CADETProcess.processModel import TubularReactorBase, StericMassAction


__all__ = ['GradientExperiment', 'plot_experiments', 'YamamotoResults', 'fit_parameters']


Expand Down Expand Up @@ -66,7 +65,7 @@ def calculate_normalized_gradient_slope(self, column_volume, total_porosity):
"""
slope = (self.c_salt_end - self.c_salt_start) / self.gradient_volume
vol_factor = column_volume - total_porosity*column_volume
vol_factor = column_volume - total_porosity * column_volume
return slope * vol_factor

def plot(self, fig=None, ax=None, sec_ax=None):
Expand Down Expand Up @@ -118,9 +117,9 @@ def yamamoto_equation(log_c_salt_at_max_M, lambda_, nu, k_eq):
DESCRIPTION.
"""
lambda_M = lambda_/1000
lambda_M = lambda_ / 1000

return np.multiply((nu + 1), log_c_salt_at_max_M) - np.log10(k_eq * lambda_M**nu * (nu + 1))
return np.multiply((nu + 1), log_c_salt_at_max_M) - np.log10(k_eq * lambda_M ** nu * (nu + 1))


class YamamotoResults:
Expand Down Expand Up @@ -154,8 +153,8 @@ def plot(self, fig=None, ax=None):
nu = self.characteristic_charge[i_p]

x = [
min(self.log_c_salt_at_max_M[:, i_p])*1.05,
max(self.log_c_salt_at_max_M[:, i_p])*0.95
min(self.log_c_salt_at_max_M[:, i_p]) * 1.05,
max(self.log_c_salt_at_max_M[:, i_p]) * 0.95
]

y = yamamoto_equation(
Expand Down Expand Up @@ -212,7 +211,7 @@ def fit_parameters(experiments, column):

for i_p in range(experiments[0].n_proteins):
c_salt_at_max = [exp.c_salt_at_max[i_p] for exp in experiments]
log_c_salt_at_max_M[:, i_p] = np.log10(np.array(c_salt_at_max)/1000)
log_c_salt_at_max_M[:, i_p] = np.log10(np.array(c_salt_at_max) / 1000)

nu[i_p], k_eq[i_p] = _fit_yamamoto(
log_c_salt_at_max_M[:, i_p], log_gradient_slope, column.binding_model.capacity
Expand All @@ -224,7 +223,7 @@ def fit_parameters(experiments, column):
yamamoto_results = YamamotoResults(
column, experiments,
log_gradient_slope, log_c_salt_at_max_M
)
)

return yamamoto_results

Expand Down Expand Up @@ -254,8 +253,7 @@ def _fit_yamamoto(log_c_salt_at_max_M, log_gradient_slope, lambda_):
def yamamoto_wrapper(c_s, nu, k_eq):
return yamamoto_equation(c_s, lambda_, nu, k_eq)

results, pcov = curve_fit(
yamamoto_wrapper, log_c_salt_at_max_M, log_gradient_slope, bounds=bounds
)
results, pcov = curve_fit(f=yamamoto_wrapper, xdata=log_c_salt_at_max_M, ydata=log_gradient_slope, bounds=bounds,
p0=(1, 1))

return results

0 comments on commit 705aef8

Please sign in to comment.