import glob
import copy
import math
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from openalea.mtg.algo import axis
from scipy import optimize, constants
from pathlib import Path
from openalea.hydroroot.display import plot
from openalea.hydroroot.read_file import read_archi_data
from openalea.hydroroot.main import root_builder, hydroroot_flow
from openalea.hydroroot import radius, flux, conductance
from openalea.hydroroot.init_parameter import Parameters
from openalea.hydroroot.conductance import set_conductances, axial
from openalea.hydroroot.water_solute_transport import pressure_calculation, pressure_calculation_no_non_permeating_solutes, \
init_some_MTG_properties, osmotic_p_peg
[docs]
def water_solute_model(parameter, df_archi =None, df_law =None,
df_cnf = None, df_JvP = None, Data_to_Optim = None, Flag_verbose = False,
data_to_use = 'all', output = None, optim_method = 'COBYLA', Flag_debug = False,
Flag_radius = True, Flag_Constraint = True, dK_constraint = -3.0e-2, Flag_w_Lpr = False,
Flag_w_cnf = False):
"""
Perform direct simulations or parameters adjustment to fit data of Jv(P) and/or cut and flow experiments.
Water and solute transport. **Works with constant radial conductivity.**
:param parameter: Parameter - (see :func: Parameters)
:param df_archi: DataFrame (None) - DataFrame with the architecture data (see below structure description)
:param df_law: DataFrame list (None) - DataFrame with the length law data (see below structure description)
:param df_cnf: DataFrame (None) - cut and flow data to fit (see below structure description)
:param df_JvP: DataFrame (None) - Jv(P) data to fit (see below structure description)
:param Data_to_Optim: string list (None) - list of parameters to adjust, if None perform direct simulation, empty list is equivalent to ['K', 'k', 'Js', 'Ps']
:param Flag_verbose: boolean (False) - if True print intermediary results, optimization details, final simulation outputs, etc.
:param data_to_use: string ('all') - data to fit either 'JvP' (Jv(P)), 'cnf' (cut and flow), or 'all' both
:param output: string (None) - if not None output filename
:param optim_method: string ('COBYLA') - solver method used in scipy.optimize.minimize see docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
:param Flag_debug: boolean (False) - if True print intermediary values of parameters adjustment
:param Flag_radius: boolean (True) - if True use diameter recorded in architecture file if present, otherwise use ref_radius
:param Flag_Constraint: boolean (True) - if True apply constraint on axial conductance 1st derivative (with COBYLA constraint may not be respected)
:param dK_constraint: float (-3.0e-2) - lower bound of the axial conductance 1st derivative if Flag_Constraint = True (with COBYLA constraint may not be respected)
:param Flag_w_Lpr: boolean (False) - set weight to 1.0 / len(list_DP) on the JvP objective function
:param Flag_w_cnf: boolean (False) - set weight to len(cut_n_flow_length) on the cut and flow objective function
:return:
- d: DataFrame with results
- g_cut: dictionary with MTG at each cut
**Data_to_Optim list**
- 'K': optimize axial conductance K
- 'k': optimize radial conductivity k
- 'Js': optimize pumping rate Js
- 'Ps': optimize permeability
- 'Klr': optimize axial conductance of laterals Klr if <> than PR
- 'klr': optimize radial conductivity of laterals klr if <> than PR
- 'sigma': optimize the reflection coefficient
['K', 'k', 'Js', 'Ps'] will optimize the four parameters in the list.
But, having too many parameters with this routine using local optimizer from scipy may not lead to any results.
**df_archi column names**
- distance_from_base_(mm), lateral_root_length_(mm), order
**df_law**
- list of 2 dataframe with the length law data: the first for the 1st order laterals on the primary root, the
2nd for the laterals on laterals whatever their order (2nd, 3rd, ...)
- column names: LR_length_mm , relative_distance_to_tip
**df_cnf column names**
- arch: sample name that must be contained in the 'input_file' of the yaml file
- dP_Mpa: column with the working cut and flow pressure (in relative to the base) if constant, may be empty see below
- J0, J1, ..., Jn: columns that start with 'J' containing the flux values, 1st the for the full root, then 1st cut, 2d cut, etc.
- lcut1, ...., lcutn: columns starting with 'lcut' containing the maximum length to the base after each cut, 1st cut, 2d cut, etc. (not the for full root)
- dP0, dP1,.., dPn: column starting with 'dP' containing the working pressure (in relative to the base) of each steps (if not constant): full root, 1st cut, 2d cut, etc.
**df_JvP column names**
- arch: sample name that must be contained in the 'input_file' of the yaml file
- J0, J1, ..., Jn: columns that start with 'J' containing the flux values of each pressure steps
- dP0, dP1,.., dPn: column starting with 'dP' containing the working pressure (in relative to the base) of each steps
**outputfile (csv)**
- column names: 'max_length', 'Jexp cnf (uL/s)', 'Jv cnf (uL/s)', 'surface (m2)', 'length (m)', 'dp', 'Jexp(P)', 'Jv(P)', 'Cbase', 'kpr', 'klr', 'Js', 'Ps', 'F cnf','F Lpr', 'x pr', 'K1st pr', 'x lr', 'K1st lr', 'K pr', 'K lr'
i.e.: max length from the cut to the base, J cnf exp, J cnf sim, root surface, total root length, pressure (Jv(P),
J exp Jv(P), J sim Jv(P), solute concentration at the base (Jv(P)), radial conductivity PR, radial conductivity LR,
pumping rate, permeability, objective fct cnf, objective fct Jv(P), x pr distance to tip for PR, K 1st guess for PR,
the same for laterals if different, then axial conductances for PR and LR.
**Remark**
- radial conductivity: single value or list of 2 values [kpr, klr]: 1st value for PR and 2nd for LR
- axial conductance: it is possible to add K for LR
"""
if Data_to_Optim is None:
Data_to_Optim = [] # direct simulation
elif len(Data_to_Optim) == 0:
Data_to_Optim = ['K', 'k', 'Js', 'Ps'] # if Data_to_Optim = []
Flag_Optim_K = ('K' in Data_to_Optim) # optimize axial conductance K
Flag_Optim_Klr = ('Klr' in Data_to_Optim) # optimize axial conductance of laterals Klr if <> than PR
Flag_Optim_klr = ('klr' in Data_to_Optim) # optimize radial conductivity of laterals klr if <> than PR
Flag_Optim_Js = ('Js' in Data_to_Optim) # optimize pumping rate Js
Flag_Optim_Ps = ('Ps' in Data_to_Optim) # optimize permeability
Flag_Optim_k = ('k' in Data_to_Optim) # optimize radial conductivity k
Flag_Optim_sigma = ('sigma' in Data_to_Optim)
if df_cnf is None:
# force on JvP data
data_to_use = 'JvP'
elif df_JvP is None:
# force on cut and flow data
data_to_use = 'cnf'
# if both are none the simulation will be done using psi_base en ext given parameter and run as one JvP data point
results = {}
g_cut = {}
tip_id = {}
S_g = []
cut_n_flow_length = []
Jexp = []
result_cv = []
def fun_constraint(x):
"""
Calculation of the constraint for the optimize.minimize solver
array of 1 column with non-negative constraints, i.e. every c[i] >= 0
:param x: array of the parameters to optimize
:return: numpy arrays
"""
n = len(axial_data[1])
c = np.ones(n - 1)
l = axial_data[0]
# inequality constraints >= 0 for the axial conductance: line below <=> (x[i+1] - x[i])/(l[i+1] - l[i]) >= dK_constraint
for i in range(n - 1):
c[i] = x[i + 1] - x[i] - dK_constraint * (l[i + 1] - l[i])
# bounds as inequality constraints needed by solver 'COBYLA' but redundant for 'SLSQP' with bounds
for i in range(n, len(x)):
c = np.append(c, (x[i] - bnds[i][0]))
c = np.append(c, (bnds[i][1] - x[i]))
return c
def fun_bound_cobyla(x):
"""
Calculation of the constraint for the optimize.minimize solver COBYLA
bounds expressed as non-negative constraint
array of 1 column with non-negative constraints, i.e. every c[i] >= 0
:param x: array of the parameters to optimize
:return: numpy arrays
"""
c = []
for i in range(len(x)):
c.append(x[i] - bnds[i][0])
c.append(bnds[i][1] - x[i])
return np.array(c)
def fun(x):
"""
Calculation of the objective function (Residual sum of squares) done on Jv(P) and CnF data
:param x: array of the parameters to optimize
:return: F (float), the objective function
"""
# Kx pr, Kx lr, Js, Ps,k lr and k pr may be optimized depends on the Flags_Optim_
_x = x * xini
if iKpr > 0:
axial_data[1] = list(
_x[:iKpr + 1]) # np array : _x[:n] means the n 1st elements index [0;n-1] but _x[n] means _x index n
if not Flag_Optim_Klr:
_axial_lr = None
else:
_axial_lr = axial_lr
_axial_lr[1] = _x[iKpr + 1:iKlr + 1]
if Flag_Optim_Js:
Js = _x[iJs]
else:
Js = J_s
if Flag_Optim_Ps:
Ps = _x[iPs]
else:
Ps = P_s
if Flag_Optim_klr:
klr = _x[iklr]
else:
klr = k[1]
if Flag_Optim_k:
kpr = _x[ikpr]
else:
kpr = k[0]
if Flag_Optim_sigma:
sigma = _x[isig]
else:
sigma = Sigma
# set new K and k in the MTG
g = set_K_and_k(g_cut, axial_data, kpr, axial_lr = _axial_lr, k_lr = klr, nr = Nb_of_roots,
nl = len(cut_n_flow_length))
# run simulation Jv(P) and CnF
JvP, F_JvP, C = Jv_P_calculation(g, sigma, Js, Ps)
JvCnf, F_JvCnF, C = Jv_cnf_calculation(g, sigma, Js, Ps)
F = F_JvP + F_JvCnF
# n = len(axial_data[1])
# c = np.ones(n - 1)
# l = axial_data[0]
# for i in range(n - 1):
# if _x[i+1] - _x[i] - dK_constraint * (l[i+1] - l[i]) < 0.0:
# F += 1.0
if Flag_debug: print('{:0.3e}'.format(F), ' '.join('{:0.3e}'.format(i) for i in _x))
result_cv.append([F, kpr, klr, Ps, Js] + axial_data[1])
return F
def fun_JvP_only(x):
"""
Calculation of the objective function (Residual sum of squares) done on Jv(P) data
:param x: array of the parameters to optimize
:return: F (float), the objective function
"""
# Kx pr, Kx lr, Js, Ps,k lr and k pr may be optimized depends on the Flags_Optim_
_x = x * xini
if iKpr > 0:
axial_data[1] = list(_x[:iKpr + 1]) # np array : _x[:n] means the n 1st elements index [0;n-1] but _x[n] means _x index n
if not Flag_Optim_Klr:
_axial_lr = None
else:
_axial_lr = axial_lr
_axial_lr[1] = _x[iKpr + 1:iKlr + 1]
if Flag_Optim_Js:
Js = _x[iJs]
else:
Js = J_s
if Flag_Optim_Ps:
Ps = _x[iPs]
else:
Ps = P_s
if Flag_Optim_klr:
klr = _x[iklr]
else:
klr = k[1]
if Flag_Optim_k:
kpr = _x[ikpr]
else:
kpr = k[0]
if Flag_Optim_sigma:
sigma = _x[isig]
else:
sigma = Sigma
# set new K and k in the MTG
g = set_K_and_k(g_cut, axial_data, kpr, axial_lr = _axial_lr, k_lr = klr, nr = Nb_of_roots,
nl = len(cut_n_flow_length))
# run simulation Jv(P)
JvP, F, C = Jv_P_calculation(g, sigma, Js, Ps)
if Flag_debug: print('{:0.3e}'.format(F), ' '.join('{:0.3e}'.format(i) for i in _x))
result_cv.append([F, kpr, klr, Ps, Js] + axial_data[1])
return F
def fun_cnf_only(x):
"""
Calculation of the objective function (Residual sum of squares) done on CnF data
:param x: array of the parameters to optimize
:return: F (float), the objective function
"""
# Kx pr, Kx lr, Js, Ps,k lr and k pr may be optimized depends on the Flags_Optim_
_x = x * xini
if iKpr > 0: axial_data[1] = list(
_x[:iKpr + 1]) # np array : _x[:n] means the n 1st elements index [0;n-1] but _x[n] means _x index n
if not Flag_Optim_Klr:
_axial_lr = None
else:
_axial_lr = axial_lr
_axial_lr[1] = _x[iKpr + 1:iKlr + 1]
if Flag_Optim_Js:
Js = _x[iJs]
else:
Js = J_s
if Flag_Optim_Ps:
Ps = _x[iPs]
else:
Ps = P_s
if Flag_Optim_klr:
klr = _x[iklr]
else:
klr = k[1]
if Flag_Optim_k:
kpr = _x[ikpr]
else:
kpr = k[0]
if Flag_Optim_sigma:
sigma = _x[isig]
else:
sigma = Sigma
# set new K and k in the MTG
g = set_K_and_k(g_cut, axial_data, kpr, axial_lr = _axial_lr, k_lr = klr, nr = Nb_of_roots,
nl = len(cut_n_flow_length))
# run simulation Jv(P)
JvCnf, F, C = Jv_cnf_calculation(g, sigma, Js, Ps)
if Flag_debug: print('{:0.3e}'.format(F), ' '.join('{:0.3e}'.format(i) for i in _x))
result_cv.append([F, kpr, klr, Ps, Js] + axial_data[1])
return F
def set_K_and_k(g, axial_pr, k_pr, axial_lr = None, k_lr = None, nr = 1, nl = 0):
"""
set the axial conductance and the radial conductivity of the uncut root and the different cut roots
The vertices of the cut roots from the entire root may have changed therefore the vertices where the cuts are made
must be set with the correct K and k see the code
:param g: (dict) - dictionnary of MTG corresponding to the entire root and the cuts
:param axial_pr: (list) - the axial conductance, list of 2 lists of floats
:param k_pr: (float) - the radial conductivity
:param axial_lr: (list) - if not None the axial conductance of the laterals, list of 2 lists of floats
:param k_lr: (float) - if not None the radial conductivity of the laterals
:param nr: (int) - number of root, because with seminals the measurements may have been done with several roots
:param nl: (int) - number of cuts
:return: g
"""
# if axial_lr is None: axial_lr = axial_pr
# if k_lr is None: k_lr = k_pr
for ig in range(nr):
g[0, ig] = set_conductances(g[0, ig], axial_pr = axial_pr, k0_pr = k_pr, axial_lr = axial_lr, k0_lr = k_lr)
# set the different cut roots
for ic in range(1, nl + 1):
for v in g[ic, ig].vertices_iter(g[0, ig].max_scale()):
vid = g[ic, ig].property('original_vid')[v]
g[ic, ig].property('K')[v] = g[0, ig].property('K')[vid]
g[ic, ig].property('k')[v] = g[0, ig].property('k')[vid]
g[ic, ig].property('K_exp')[v] = g[0, ig].property('K_exp')[
vid] # needed in pressure_calculation if Cpeg because calculation of K
g[ic, ig].property('k0')[v] = g[0, ig].property('k0')[vid]
# difference with the resistance network we do not set k = K the real boundary condition is used
# with the help of label 'cut' see pressure_calculation
return g
def Jv_P_calculation(g, sigma, Js, Ps):
"""
Perform the calculation of the data Jv(P), i.e. for different pressure difference list_DP
Most of the variables are global variables, only the variables that change at each calculation g (MTG) or
that are able to be optimized (sigma, Js, Ps) are passed in arguments
The change of K and k have been taken into account in function set_K_and_k
Return
JvP a dictionnary of outgoing sap flux at each pressure step and for each roots (for the case they are several seminals)
C a dictionnary of the concentration map at each pressure step and for each roots (for the case they are several seminals)
F the objective fonction
:param g: (dict) - dictionnary of MTG corresponding to the entire root and the cuts
:param sigma: (float) - the reflection coefficient
:param Js: (float) - the pumping rate
:param Ps: (float) - the permeability
:return: JvP (dict), C {dict}, F (float)
"""
JvP = {}
C = {}
F = 0.0
data = None
row = None
col = None
for idP in range(len(list_DP)):
Jv = 0.0
for ig in range(Nb_of_roots):
g[0, ig] = flux.flux(g[0, ig], psi_e = psi_base + list_DP[idP], psi_base = psi_base,
invert_model = True)
g[0, ig] = init_some_MTG_properties(g[0, ig], tau = Js, Cini = Cini, t = 1, Ps = Ps)
nb_v = g[0, ig].nb_vertices()
Fdx = 1.0
Fdx_old = 1.
Jv_old = 1.
# Newton-Raphson schemes: in pressure_calculation_no_non_permeating_solutes calculation of dx,
# array with dP and dC variation of the variables between two Newton step. Then the Newton scheme stops when
# Fdx > eps see below
while Fdx > eps:
g[0, ig], dx, data, row, col = pressure_calculation_no_non_permeating_solutes(g[0, ig],
sigma = sigma,
Ce = Ce,
Pe = parameter.exp[
'psi_e'],
Pbase = parameter.exp[
'psi_base'],
Cse = Cse,
dP = list_DP[idP],
C_base = None)
Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
JvP[idP, ig] = g[0, ig].property('J_out')[1]
if abs(JvP[idP, ig] - Jv_old) < 1.0e-4:
break
if abs(Fdx - Fdx_old) < eps:
break
Fdx_old = Fdx
Jv_old = JvP[idP, ig]
Jv += JvP[idP, ig]
C[idP, ig] = copy.deepcopy(g[0, ig].property('C'))
F += w_Lpr * (Jv - list_Jext[idP]) ** 2.0
return JvP, F, C
def Jv_cnf_calculation(g, sigma, Js, Ps):
"""
Perform the calculation of the data CnF, i.e. for different cut length cut_n_flow_length
Most of the variables are global variables, only the variables that change at each calculation g (MTG) or
that are able to be optimizes (sigma, Js, Ps) are passed in arguments
Return
JvCnf a dictionnary of outgoing sap flux at each cut step and for each roots (for the case they are several seminals)
C a dictionnary of the concentration map at each cut step and for each roots (for the case they are several seminals)
F the objective fonction
:param g: (dict) - dictionnary of MTG corresponding to the entire root and the cuts
:param sigma: (float) - the reflection coefficient
:param Js: (float) - the pumping rate
:param Ps: (float) - the permeability
:return: JvCnf (dict), C {dict}, F (float)
"""
ic = 0
JvCnf = {}
C = {}
F = 0.0
Jv = 0.0
data = row = col = None
for ig in range(Nb_of_roots):
g[ic, ig] = flux.flux(g[ic, ig], psi_e = psi_base + DP_cnf[ic], psi_base = psi_base, invert_model = True)
g[ic, ig] = init_some_MTG_properties(g[ic, ig], tau = Js, Cini = Cini, t = 1, Ps = Ps)
nb_v = g[ic, ig].nb_vertices()
Fdx = 1.0
Fdx_old = 1.
Jv_old = 1.
# Newton-Raphson schemes: in pressure_calculation_no_non_permeating_solutes calculation of dx,
# array with dP and dC variation of the variables between two Newton step. Then the Newton scheme stops when
# Fdx > eps see below
while Fdx > eps:
# use pressure_calculation_no_non_permeating_solutes because the root is uncut so no PEG enter the root
g[ic, ig], dx, data, row, col = pressure_calculation_no_non_permeating_solutes(g[ic, ig], sigma = sigma,
Ce = Ce,
Pe = parameter.exp[
'psi_e'],
Pbase = parameter.exp[
'psi_base'],
Cse = Cse,
dP = DP_cnf[ic])
Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
JvCnf[ic, ig] = g[ic, ig].property('J_out')[1]
if abs(JvCnf[ic, ig] - Jv_old) < 1.0e-4:
break
if abs(Fdx - Fdx_old) < eps:
break
Fdx_old = Fdx
Jv_old = JvCnf[ic, ig]
Jv += JvCnf[ic, ig]
C[ic, ig] = copy.deepcopy(g[ic, ig].property('C'))
F += w_cnf * (Jv - Jexp[ic]) ** 2.0
for ic in range(1, len(cut_n_flow_length) + 1):
Jv = 0.0
data = row = col = None
for ig in range(Nb_of_roots):
g[ic, ig] = flux.flux(g[ic, ig], psi_e = psi_base + DP_cnf[ic], psi_base = psi_base,
invert_model = True)
g[ic, ig] = init_some_MTG_properties(g[ic, ig], tau = Js, Cini = Cini, Cpeg_ini = Cpeg_ini, t = 1,
Ps = Ps)
nb_v = g[ic, ig].nb_vertices()
Fdx = 1.0
Fdx_old = 1.
Jv_old = 1.
# Newton-Raphson schemes: in pressure_calculation calculation of dx,
# array with dP, dC and dCpeg (if any) variation of the variables between two Newton step. Then the Newton scheme stops when
# Fdx > eps see below
while Fdx > eps:
g[ic, ig], dx, data, row, col = routine_calculation(g[ic, ig], sigma = sigma,
Ce = Ce, Pe = parameter.exp['psi_e'],
Pbase = parameter.exp['psi_base'],
Cse = Cse, dP = DP_cnf[ic])
Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
JvCnf[ic, ig] = g[ic, ig].property('J_out')[1]
# if Flag_debug: print local_j, Fdx, (Fdx - Fdx_old)
if abs(JvCnf[ic, ig] - Jv_old) < 1.0e-4:
break
if abs(Fdx - Fdx_old) < eps:
break
Fdx_old = Fdx
Jv_old = JvCnf[ic, ig]
Jv += JvCnf[ic, ig]
C[ic, ig] = copy.deepcopy(g[ic, ig].property('C'))
F += w_cnf * (Jv - Jexp[ic]) ** 2.0
return JvCnf, F, C
# architecture file to dataframe
index=''
# if (df_archi is None) and parameter.archi['read_architecture']:
# # architecture with filename in aqua team format
# archi_f = glob.glob(parameter.archi['input_dir'] + parameter.archi['input_file'][0])
# archi_f = archi_f[0]
# df_archi = read_archi_data(archi_f) if parameter.archi['read_architecture'] else None
# index = archi_f.replace(glob.glob(parameter.archi['input_dir'])[0], "")
if parameter.archi['read_architecture']:
# architecture with filename in aqua team format
# archi_f = glob.glob(parameter.archi['input_dir'] + parameter.archi['input_file'][0])
# archi_f = archi_f[0]
fname = str(Path(parameter.archi['input_dir']) / parameter.archi['input_file'][0]) # deals with path in windows
archi_f = glob.glob(fname) # deals with wildcards ? deprecated ???
archi_f = archi_f[0]
if df_archi is None: df_archi = read_archi_data(archi_f)
index = archi_f.replace(glob.glob(parameter.archi['input_dir'])[0], "")
if len(parameter.archi['input_file']) > 0:
index = parameter.archi['input_file'][0]
# length law data: override if necessary
if df_law is not None:
parameter.archi['length_data'] = df_law
# dataframe used to save and export results: cnf and Jv(P)
_col_names = ['max_length', 'Jexp cnf (uL/s)', 'Jv cnf (uL/s)', 'surface (m2)', 'length (m)']
results = {}
for key in _col_names:
results[key] = []
_col_names2 = ['dp', 'Jexp(P)', 'Jv(P)', 'Cbase']
results2 = {}
for key in _col_names2:
results2[key] = []
############################
# get value from yaml file
############################
delta = parameter.archi['branching_delay'][0]
nude_length = parameter.archi['nude_length'][0]
seed = parameter.archi['seed'][0]
axfold = parameter.output['axfold'][0]
radfold = parameter.output['radfold'][0]
# Conductancies: mananging the fact there are or not different values between the primary and laterals
# and the fact there are multiply by axfold and radfold
k = []
if type(parameter.hydro['k0']) != list:
k.append(parameter.hydro['k0'] * radfold)
k.append(None)
else:
if len(parameter.hydro['k0']) > 1:
k.append(parameter.hydro['k0'][0] * radfold)
k.append(parameter.hydro['k0'][1] * radfold)
else:
k.append(parameter.hydro['k0'][0] * radfold)
k.append(None)
exp_axial = parameter.hydro['axial_conductance_data']
axial_data = ([exp_axial[0], exp_axial[1]])
axial_data = list(axial(axial_data, axfold))
if len(exp_axial) == 4:
axial_lr = ([exp_axial[2], exp_axial[3]])
axial_lr = list(axial(axial_lr, axfold))
else:
axial_lr = None #copy.deepcopy(axial_data)
J_s = parameter.solute['J_s']
P_s = parameter.solute['P_s']
Cse = parameter.solute['Cse'] * 1e-9 # mol/m3 -> mol/microL, external permeating solute concentration
Ce = parameter.solute['Ce'] * 1e-9 # mol/m3 -> mol/microL, external non-permeating solute concentration
Cini = Cse # initialization solute concentration into the xylem vessels
Cpeg_ini = Ce # initialization non-permeating solute concentration into the xylem vessels: not 0.0 because more num instability
Sigma = parameter.solute['Sigma'] # reflection coefficient, fixed in this script
Pi_e_peg = osmotic_p_peg(Ce, unit_factor = 8.0e6) # from Ce mol/microL to g/g, external osmotic pressure of non-permeating in MPa
data = None
row = None
col = None
w_cnf = w_Lpr = 1. # weight on cnf cost function
# functions that resolve the matrix system used in the Newton-Raphson scheme
# different function depending on the presence of non-permeating solute, because there is one unknown less Cpeg
routine_calculation = None
if Ce <= 0.:
# no non-permeating solute present
routine_calculation = pressure_calculation_no_non_permeating_solutes
else:
routine_calculation = pressure_calculation
# the objective function calculation to call depending on the data we fit
if data_to_use == "cnf":
fun_objective = fun_cnf_only
elif data_to_use == "JvP":
fun_objective = fun_JvP_only
else:
fun_objective = fun
dK_constraint_max = 6.0e-2 # deprecated
_tol = 5.0e-7 # does not have significant impact !!?? used in some minimize.optimize solver
eps = 1.0e-9 # global: stop criterion for the Newton-Raphson loop in Jv_P_calculation and Jv_cnf_calculation
# Parameter bounds
Kbnds = (1.0e-10, np.inf) # axial conductance
kbnds = (0.01, np.inf) # radial conductivity
Jbnds = (1e-15, np.inf) # Js
Pbnds = (1e-15, np.inf) # Ps
psi_base = parameter.exp['psi_base']
# default value for the pressure difference between the external medium and the base
DP_cnf = []
DP_cnf.append(parameter.exp['psi_e'] - psi_base)
# variables used for the results output see end of script
K = {}
K['x pr'] = axial_data[0]
K['K1st pr'] = axial_data[1]
dK1st = pd.DataFrame(K, columns = ['x pr', 'K1st pr'])
K = {}
if axial_lr:
K['x lr'] = axial_lr[0]
K['K1st lr'] = axial_lr[1]
else:
K['x lr'] = axial_data[0]
K['K1st lr'] = axial_data[1]
dKlr1st = pd.DataFrame(K, columns = ['x lr', 'K1st lr'])
# # architecture file to dataframe
# df_archi = read_archi_data(archi_f) if parameter.archi['read_architecture'] else None
# index = archi_f.replace(glob.glob(parameter.archi['input_dir'])[0], "")
# read the data measurements from data base, cut-n-flow: flux, cut length and pressure difference
if df_cnf is not None:
for key in df_cnf['arch']:
if str(key).lower() in index.lower():
_list = df_cnf[df_cnf.arch == key].filter(regex = '^J').dropna(axis = 1).values.tolist()
Jexp = _list[0] # basal output flux
_list = df_cnf[df_cnf.arch == key].filter(regex = '^lcut').dropna(axis = 1).values.tolist()
cut_n_flow_length = _list[0] # cut lengthes
_list = df_cnf[df_cnf.arch == key].filter(regex = '^dP').dropna(axis = 1).values.tolist()
# the pressure difference is usually constant but sometimes, due to flow meter saturation, it may change
# in that case a list of values is given
if len(_list[0]) != 0:
DP_cnf = _list[0]
else:
DP_cnf = []
DP_cnf.append(parameter.exp['psi_e'] - psi_base) # for compatibility reason with first analysis on arabidopsis
if len(DP_cnf) < len(cut_n_flow_length) + 1: # if constant we create the list with the constant value
for i in range(1, len(cut_n_flow_length) + 1): DP_cnf.append(DP_cnf[0])
parameter.exp['psi_e'] = psi_base + DP_cnf[0]
# read the data measurements from data base Jv(P): flux, pressure
if df_JvP is not None:
for key in df_JvP['arch']:
if str(key).lower() in index.lower():
_list = df_JvP[df_JvP.arch == key].filter(regex = '^J').dropna(axis = 1).values.tolist()
list_Jext = _list[0] # basal output flux
_list = df_JvP[df_JvP.arch == key].filter(regex = '^dP').dropna(axis = 1).values.tolist()
list_DP = _list[0] # delta pressure
## below juste to get data above a minimum dP
# dlpr = pd.DataFrame(list(zip(list_DP, list_Jext)), columns = ['dP', 'Jv'])
# dlpr = dlpr.sort_values('dP')[dlpr['dP']>0.05]
# list_DP = list(dlpr['dP'])
# list_Jext = list(dlpr['Jv'])
else:
list_Jext = [parameter.exp['Jv']]
list_DP = [parameter.exp['psi_e'] - psi_base]
if Flag_w_Lpr: w_Lpr = 1.0 / len(list_DP)
if Flag_w_cnf: w_cnf = len(cut_n_flow_length)
# building the MTG
###################
Nb_of_roots = 2 if "-L" in index else 1 # sometimes thera are 2 roots for a given measurement with seminals
if df_archi is None:
primary_length = parameter.archi['primary_length'][0]
else:
primary_length = 0.
_length = 0
_surface = 0
for ig in range(Nb_of_roots):
if ig == 1:
f2 = archi_f.replace("-L", "-R")
df_archi = read_archi_data(f2) if parameter.archi['read_architecture'] else None
g_cut[0, ig], _p, _l, _s, _seed = root_builder(df = df_archi,
primary_length = parameter.archi['primary_length'][0],
seed = parameter.archi['seed'][0],
delta = parameter.archi['branching_delay'][0],
nude_length = parameter.archi['nude_length'][0],
segment_length = parameter.archi['segment_length'],
length_data = parameter.archi['length_data'],
branching_variability = parameter.archi['branching_variability'],
order_max = parameter.archi['order_max'],
order_decrease_factor = parameter.archi['order_decrease_factor'],
ref_radius = parameter.archi['ref_radius'],
Flag_radius = Flag_radius)
if _p > primary_length: primary_length = _p
_length += _l
_surface += _s
base = {}
for v in g_cut[0, ig]:
base[v] = next(axis(g_cut[0, ig], v))
g_cut[0, ig].properties()['axisbase'] = base
S_g.append(_s)
# case where the primary is shorter than laterals
max_length = primary_length
mylength = g_cut[0, ig].property('mylength')
if max(mylength.values()) > max_length: max_length = max(mylength.values())
# set conductance
g_cut[0, ig] = set_conductances(g_cut[0, ig], axial_pr = axial_data, k0_pr = k[0], axial_lr = axial_lr,
k0_lr = k[1])
# flux calculation without solute transport a way to initialize
g_cut[0, ig] = flux.flux(g_cut[0, ig], psi_e = psi_base + DP_cnf[0], psi_base = psi_base, invert_model = True)
# add properties specific to solute transport
g_cut[0, ig].add_property('C') # permeating solute concentration
g_cut[0, ig].add_property('Cpeg') # non-permeating solute concentration needed if cut-n-flow with them in the medium
g_cut[0, ig].add_property('theta') # see init_some_MTG_properties
g_cut[0, ig].add_property('J_s') # see init_some_MTG_properties, at a certain time I tried varying Js with C
g_cut[0, ig].add_property('P_s') # see init_some_MTG_properties, at a certain time I tried varying Js with C
g_cut[0, ig].add_property('original_vid') # the indices change between the full root and the cut root a way
# to retrieve the original index see set_K_and_k
g_cut[0, ig].add_property('mu') # the viscosity of the sap because could change from the water value when
# non-permeating solute enter the cut roots
# a simple record of the original vertex number in the full architecture
# do this because below when we cut we reindex because equations system is resolved in matrix form on the
# so the vertices need to have proper indices
# MTG
d = {vid: vid for vid in g_cut[0, ig].vertices(g_cut[0, ig].max_scale())}
g_cut[0, ig].properties()['original_vid'] = d
# ############ longitudinal CUTS ####################################
ic = 1
for cut_length in cut_n_flow_length:
# print(cut_length)
tip_id[ic, ig] = flux.segments_at_length(g_cut[0, ig], cut_length, dl = parameter.archi['segment_length'])
g_cut[ic, ig] = flux.cut(g_cut[0, ig], cut_length, parameter.archi['segment_length'])
for i in tip_id[ic, ig]:
v = g_cut[0, ig].parent(i)
g_cut[ic, ig].property('label')[v] = 'cut' # labelling the vertices at cut ends
# Below reindex because the system is resolved in matrix form on the MTG so the vertices need to have proper indices
g_cut[ic, ig].reindex()
i = 0
tip_id[ic, ig] = [] # reinitializing because the cut can be at ramification then one parent for 2 different cut vertices
for vid in g_cut[ic, ig].vertices_iter(g_cut[ic, ig].max_scale()):
if g_cut[ic, ig].label(vid) == 'cut':
tip_id[ic, ig].append(vid)
i += 1
g_cut[ic, ig], surface = radius.compute_surface(g_cut[ic, ig])
S_g.append(surface)
ic += 1
# Optimization
##############
# the parameter are normalized with their inital values to limit scale effect between them, not the best the best would
# be to write the equation in dimensionless form but historicaly hydroroot was not written this way
iKpr = iKlr = iJs = iPs = ikpr = iklr = 0 # indices used to select the correct parameters in the array x, see fun for instance
if Data_to_Optim:
# setting bounds and initial values
ix = -1
bnds = [] # list of tuple for bounds
xini_list = [] # list of initial values of parameters
if Flag_Optim_K:
for var in axial_data[1]:
xini_list.append(var)
ix += len(axial_data[1]) # be careful with axial_data and axial_lr the indices will be use as end of list interval selection => +1
iKpr = int(ix)
if Flag_Optim_Klr:
if axial_lr is None: axial_lr = copy.deepcopy(axial_data)
for var in axial_lr[1]:
xini_list.append(var)
ix += len(axial_lr[1])
iKlr = int(ix)
if Flag_Optim_K or Flag_Optim_Klr:
for i, val in enumerate(xini_list):
bnds.append(Kbnds)
if Flag_Optim_Js:
xini_list.append(J_s)
bnds.append(Jbnds)
ix += 1
iJs = int(ix)
if Flag_Optim_Ps:
xini_list.append(P_s)
bnds.append(Pbnds)
ix += 1
iPs = int(ix)
if Flag_Optim_k:
xini_list.append(k[0])
bnds.append(kbnds)
ix += 1
ikpr = int(ix)
if Flag_Optim_klr:
if k[1] is None: k[1] = copy.deepcopy(k[0])
xini_list.append(k[1])
bnds.append(kbnds)
ix += 1
iklr = int(ix)
if Flag_Optim_sigma:
xini_list.append(Sigma)
if Sigma > 0.0:
b = 1.0 / Sigma
else:
b = 1.0
bnds.append((0.0, b))
ix += 1
isig = int(ix)
xini = np.array(xini_list)
x = np.ones(len(xini)) # the array of parameter that will be optimized, equal unity because we optimize the
# the parameters normalized by their initial value
# array used for constraints see optimize.minimize doc
n = len(x)
n1 = len(axial_data[1])
# linear constraints lb <= A.dot(x) <= ub
A = np.zeros((n, n))
lb = np.full(n, -np.inf)
ub = np.full(n, np.inf)
l = copy.deepcopy(parameter.hydro['axial_conductance_data'][0])
if Flag_Optim_Klr:
l.append(0)
l.append(0)
if Flag_Optim_k: l.append(0)
if Flag_Optim_Klr: l.append(0)
if Flag_Optim_K and Flag_Constraint:
a = dK_constraint # constraint on the 1st derivative
for i in range(n1 - 1): # downward derivative
A[i, i] = -1.
A[i, i + 1] = 1.
lb[i] = a * (l[i + 1] - l[i])
# ub[i] = dK_constraint_max * (l[i + 1] - l[i])
ineq_cons = ({'type': 'ineq', 'fun': fun_constraint}) # !! works for K, k, Ps and Js optimized
else:
# for the COLBYLA solver bounds are not managed as other see fun_bound_cobyla
ineq_cons = {'type': 'ineq', 'fun': fun_bound_cobyla}
constraints = optimize.LinearConstraint(A, lb, ub) if Flag_Constraint else None
if optim_method == 'COBYLA':
res = optimize.minimize(fun_objective, x, method = optim_method, constraints = [ineq_cons])
elif optim_method == 'SLSQP':
res = optimize.minimize(fun_objective, x, bounds = bnds, method = optim_method, constraints = [ineq_cons],
options = {'ftol': 1.0e-9, 'eps': 1e-1})
else:
res = optimize.minimize(fun_objective, x, bounds = bnds, method = optim_method)
# res = optimize.minimize(fun_objective, x, bounds = bnds, method = 'trust-constr', options={'finite_diff_rel_step': 1e-1})
# res = optimize.minimize(fun_objective, x, method='TNC', bounds = bnds)
# res = optimize.minimize(fun_objective, x, bounds = bnds, options={'ftol': _tol, 'eps': 1e-1})
# res = optimize.minimize(fun_objective, x, bounds = bnds, method='nelder-mead', options={'fatol': 1.0e-9})
# optimization results to parameters
n = len(axial_data[1])
_x = res.x * xini
if Flag_Optim_K:
axial_data[1] = list(_x[:iKpr + 1])
if Flag_Optim_Klr:
axial_lr[1] = _x[iKpr + 1:iKlr + 1]
if Flag_Optim_Js: J_s = _x[iJs]
if Flag_Optim_Ps: P_s = _x[iPs]
if Flag_Optim_klr:
k[1] = _x[iklr]
# else:
# k[1] = None
if Flag_Optim_k:
k[0] = _x[ikpr]
if Flag_Optim_sigma:
Sigma = _x[isig]
if Flag_verbose: print(res.x)
# Direct simulation with the optimized values or the values from the yaml file if no optimization asked
g_cut = set_K_and_k(g_cut, axial_data, k[0], axial_lr = axial_lr, k_lr = k[1], nr = Nb_of_roots,
nl = len(cut_n_flow_length))
if data_to_use in ['all', 'JvP']:
JvP, F_JvP, C = Jv_P_calculation(g_cut, Sigma, J_s, P_s)
for idP in range(len(list_DP)):
Jv = 0.0
C_base = 0.0
for ig in range(Nb_of_roots):
# C_base here is in the middle of the 1st MTG element because the boundary condition chosen here is
# dC/dx = 0, so the concentration at the root boundary is the same.
# if there are several roots (as when the experiment is done with 2 seminals), since the MTG elements
# are equals the tital C_base is the average
C_base += C[idP, ig][1] / float(Nb_of_roots)
Jv += JvP[idP, ig]
results2['dp'].append(list_DP[idP])
results2['Jv(P)'].append(Jv)
results2['Jexp(P)'].append(list_Jext[idP])
results2['Cbase'].append(C_base * 1e9)
if data_to_use in ['all', 'cnf']:
JvCnf, F_cnf, C = Jv_cnf_calculation(g_cut, Sigma, J_s, P_s)
for ic in range(len(cut_n_flow_length) + 1):
_surface = 0.
_length = 0.
Jv = 0.0
C_base = 0.0
for ig in range(Nb_of_roots):
g_cut[ic, ig], _s = radius.compute_surface(g_cut[ic, ig])
_l = g_cut[ic, ig].nb_vertices(scale = 1) * parameter.archi['segment_length']
_surface += _s
_length += _l
Jv += JvCnf[ic, ig]
C_base += C[ic, ig][1] / float(Nb_of_roots) # C_base calculation see above JvP case
if ic > 0: max_length = cut_n_flow_length[ic - 1]
results['max_length'].append(max_length)
results['length (m)'].append(_length)
results['surface (m2)'].append(_surface)
results['Jv cnf (uL/s)'].append(Jv)
results['Jexp cnf (uL/s)'].append(Jexp[ic])
## just some parameter calculations for display
for ic in range(int(len(g_cut)/Nb_of_roots)): # FB 250328: added loop over cut
js_tot = 0
for ig in range(Nb_of_roots):
g_cut[ic, ig].add_property('DP')
g_cut[ic, ig].add_property('DC')
g_cut[ic, ig].add_property('jsurf')
g_cut[ic, ig].add_property('js')
g_cut[ic, ig].add_property('DPi')
DP = g_cut[ic, ig].property('DP')
DC = g_cut[ic, ig].property('DC')
DPi = g_cut[ic, ig].property('DPi')
jsurf = g_cut[ic, ig].property('jsurf')
psi_in = g_cut[ic, ig].property('psi_in')
js = g_cut[ic, ig].property('js')
C = g_cut[ic, ig].property('C')
length = g_cut[ic, ig].property('length')
_radius = g_cut[ic, ig].property('radius')
j = g_cut[ic, ig].property('j')
for v in g_cut[ic, ig].vertices_iter(scale = 1):
js[v] = _radius[v] * 2 * np.pi * length[v] * (J_s + P_s * (Cse-C[v]) * 1e9)
js_tot += js[v]
DC[v] = -(Cse - C[v])*1e9
DPi[v] = DC[v] * constants.R * 293 * 1e-6 # 1e-6 because in MPa
DP[v] = parameter.exp['psi_base'] + DP_cnf[0] - psi_in[v]
# psi_in[v] -= psi_base # just to put in relative pressure
jsurf[v] = j[v] / (_radius[v] * 2 * np.pi * length[v])
dr = pd.DataFrame()
dr2 = pd.DataFrame()
F = F2 = 0.0
if Flag_verbose:
pd.set_option('display.max_columns', None)
pd.set_option('display.expand_frame_repr', False)
if data_to_use in ['all', 'JvP']:
dr2 = pd.DataFrame(results2, columns = _col_names2)
# dr2.sort_values(['dp'], inplace=True)
j = np.array(dr2.loc[:, ['Jv(P)', 'Jexp(P)']])
F2 = w_Lpr * np.sum(np.diff(j) ** 2.0)
if Flag_verbose:
print('****** JvP ******')
print(dr2)
if data_to_use in ['all', 'cnf']:
dr = pd.DataFrame(results, columns = _col_names)
j = np.array(dr.loc[:, ['Jv cnf (uL/s)', 'Jexp cnf (uL/s)']])
F = w_cnf * np.sum(np.diff(j) ** 2.0)
if Flag_verbose:
print('****** cut-n-flow ******')
print(dr)
d = pd.concat([dr, dr2], axis = 1).fillna("")
X = {}
X['kpr'] = [k[0]]
if k[1] is None:
X['klr'] = [k[0]]
else:
X['klr'] = [k[1]]
X['Js'] = [J_s]
X['Ps'] = [P_s]
X['F cnf'] = [F]
X['F Lpr'] = [F2]
dX = pd.DataFrame(X, columns = ['kpr', 'klr', 'Js', 'Ps', 'F cnf', 'F Lpr'])
K = {}
# K['x pr'] = axial_data[0]
K['K pr'] = axial_data[1]
dK = pd.DataFrame(K, columns = ['K pr'])
Klr = {}
if axial_lr:
Klr['K lr'] = axial_lr[1]
else:
Klr['K lr'] = axial_data[1]
dKlr = pd.DataFrame(Klr, columns = ['K lr'])
d = pd.concat([d, dX, dK1st, dK, dKlr1st, dKlr], axis = 1).fillna("")
if output is not None: d.to_csv(output, index = False)
if Flag_verbose:
print('****** End ******')
print('objective functions: ', 'F cnf: {:0.2e}'.format(F), 'F JvP: {:0.2e}'.format(F2), 'F tot: {:0.2e}'.format(F+F2))
print(index, ',', 'k: {:0.2f}'.format(k[0]), ',', 'Js: {:0.2e}'.format(J_s), ',', 'Ps: {:0.2e}'.format(P_s), ', K: [',
', '.join('{:0.2e}'.format(i) for i in axial_data[1]), ']')
return d, g_cut
[docs]
def pure_hydraulic_model(parameter = Parameters(), df_archi = None, df_law =None, df_exp = None,
Data_to_Optim = None, output = None, Flag_verbose = False,
Flag_radius = False, Flag_Constraint = True, dK_constraint = -3.0e-2, dk_max = 0.1):
"""
Perform direct simulations or parameters adjustment to fit data of cut and flow experiment.
Water transport only, electrical network analogy
:param parameter: Parameter - (see :func: Parameters)
:param df_archi: DataFrame (None) - DataFrame with the architecture data (see below structure description)
:param df_law: DataFrame list (None) - DataFrame with the length law data (see below structure description)
:param df_exp: DataFrame (None) - data to fit
:param Data_to_Optim: string list (None) - list of parameters to adjust, if None perform direct simulation, ['K', 'k']
:param output: string (None) - if not None output filename
:param Flag_verbose: boolean (False) - if True print intermediary results
:param Flag_radius: boolean (False) - if True use diameter recorded in architecture file if present, otherwise use ref_radius
:param Flag_Constraint: boolean (True) - if True apply constraint on axial conductance 1st derivative
:param dK_constraint: float (-3.0e-2) - lower bound of the axial conductance 1st derivative if Flag_Constraint = True
:param dk_max: float (0.1) - the convergence criteria on
:return:
- df: DataFrame with results
- g_cut: dictionary with MTG at each cut
**df_archi column names**
- distance_from_base_(mm), lateral_root_length_(mm), order
**df_law**
- list of 2 dataframe with the length law data: the first for the 1st order laterals on the primary root, the
2nd for the laterals on laterals whatever their order (2nd, 3rd, ...)
- column names: LR_length_mm , relative_distance_to_tip
**The adjustment is performed as follows**
1. pre-optimization with the adjustment of axfold and radfold, K and k factor, if only k adjustment is asked then step 1 is not performed
2. loop of two successive adjustments: 1st K adjustment then k adjustment. The loop stop when change of k is below dk_max
**Data_to_Optim list of string**
- 'K': optimize axial conductance K only
- 'k': optimize radial conductivity k only
- [ ] empty list or ['K', 'k']: optimize K and k
**df_exp column names**
- arch: sample name that must be contained in the 'input_file' of the yaml file
- J0, ..., Jn: columns containing the flux values of the full root, 1st cut, 2d cut, etc.
- lcut1, ...., lcutn: columns containing the maximum length to the base after each cut, 1st cut, 2d cut, etc. (the primary length of the full root is calculated from the architecture)
**outputfile**
- column names: 'plant', 'cut length (m)', 'max_length', 'k (10-8 m/s/MPa)', 'length (m)', 'surface (m2)', 'Jv (uL/s)', 'Jexp (uL/s)'
- if 'K' in Data_to_Optim add the following: 'x', 'K 1st', 'K optimized' the initial and adjusted K(x)
**Remark**
The routine is designed to work with a single value (float) for parameter.hydro['k0'].
**example**
>>> parameter = Parameters()
>>> filename='parameters_fig-2-B.yml'
>>> parameter.read_file(filename)
>>> fn = 'data/arabido_cnf_data.csv'
>>> df_exp = pd.read_csv(fn, sep = ',', keep_default_na = True)
>>> df = pure_hydraulic_model(parameter,df_exp=df_exp, Flag_verbose=True, Data_to_Optim = ['k', 'K'])
"""
if Data_to_Optim is None:
Data_to_Optim = [] # direct simulation
elif len(Data_to_Optim) == 0:
Data_to_Optim = ['K', 'k'] # if Data_to_Optim = []
Flag_Optim_K = ('K' in Data_to_Optim) # optimize axial conductance K
Flag_Optim_k = ('k' in Data_to_Optim) # optimize radial conductivity k
g_cut = {}
tip_id = {}
cut_n_flow_length = []
_tol = 1.0e-9
def hydro_calculation(g, axfold = 1., radfold = 1., axial_data = None, k_radial = None, psi_base = 0.1, psi_e = 0.1):
if axial_data is None: axial_data = parameter.hydro['axial_conductance_data']
if k_radial is None: k_radial = parameter.hydro['k0']
# compute axial & radial
Kexp_axial_data = conductance.axial(axial_data, axfold)
k_radial_data = conductance.radial(k_radial, axial_data, radfold)
## Step function
# k_radial_data = conductance.radial_step(k_radial,3.0,x_step = 0.02, dx = parameter.archi['segment_length'], scale = radfold)
# compute local jv and psi, global Jv, Keq
g, Keq, Jv = hydroroot_flow(g, segment_length = parameter.archi['segment_length'],
k0 = k_radial,
Jv = _Jv[0],
psi_e = psi_e,
psi_base = psi_base,
axial_conductivity_data = Kexp_axial_data,
radial_conductivity_data = k_radial_data)
return g, Keq, Jv
def fun1(x):
"""
Simulation of the flux at the different cut lengths according to the new parameter value
Implementation 1: only axfold (Kx factor) and radfold (k radial factor) are changed
:param x: the array of adjusted parameters
:return: F the sum((Jv - Jv_exp) ** 2.0)
"""
axfold = x[0]
radfold = x[-1]
g_cut['tot'], Keq, Jv = hydro_calculation(g_cut['tot'], radfold = radfold, axfold = axfold, psi_base = psi_base,
psi_e = psi_base + DP_cnf[0])
F = (Jv - _Jv[0]) ** 2.0
count = 1
for cut_length in cut_n_flow_length:
# _g = g_cut[str(cut_length)].copy() # not necessary and time-consuming
for vid in g_cut[str(cut_length)].vertices_iter(g_cut['tot'].max_scale()):
g_cut[str(cut_length)].property('K')[vid] = g_cut['tot'].property('K')[vid]
g_cut[str(cut_length)].property('k')[vid] = g_cut['tot'].property('k')[vid]
for i in tip_id[str(cut_length)]:
v = g_cut['tot'].parent(i)
g_cut[str(cut_length)].property('k')[v] = g_cut[str(cut_length)].property('K')[v]
g_cut[str(cut_length)] = flux.flux(g_cut[str(cut_length)], Jv = _Jv[count], psi_e = psi_base + DP_cnf[count], psi_base = psi_base,
invert_model = True, cut_and_flow = True)
Jv = g_cut[str(cut_length)].property('J_out')[1]
F += (Jv - _Jv[count]) ** 2.0
count += 1
return F
def fun2(x):
"""
Simulation of the flux at the different cut lengths according to the new parameter value
Implementation 2: only axial_data is changed
:param x: the array of adjusted parameters
:return: F the sum((Jv - Jv_exp) ** 2.0)
"""
# k0 = parameter.hydro['k0']
axial_data = copy.deepcopy(parameter.hydro['axial_conductance_data'])
axial_data[1] = list(x)
g_cut['tot'], Keq, Jv = hydro_calculation(g_cut['tot'], k_radial = k0 ,axial_data = axial_data, psi_base = psi_base,
psi_e = psi_base + DP_cnf[0])
F = (Jv - _Jv[0])**2.0
count = 1
for cut_length in cut_n_flow_length:
# _g = g_cut[str(cut_length)].copy() # not necessary and time-consuming
for vid in g_cut[str(cut_length)].vertices_iter(g_cut['tot'].max_scale()):
g_cut[str(cut_length)].property('K')[vid] = g_cut['tot'].property('K')[vid]
g_cut[str(cut_length)].property('k')[vid] = g_cut['tot'].property('k')[vid]
for i in tip_id[str(cut_length)]:
v = g_cut['tot'].parent(i)
g_cut[str(cut_length)].property('k')[v] = g_cut[str(cut_length)].property('K')[v]
g_cut[str(cut_length)] = flux.flux(g_cut[str(cut_length)], Jv = _Jv[count], psi_e = psi_base + DP_cnf[count], psi_base = psi_base,
invert_model = True, cut_and_flow = True)
Jv = g_cut[str(cut_length)].property('J_out')[1]
F += (Jv - _Jv[count])**2.0
count += 1
return F
def fun3(x):
"""
Simulation of the flux at the different cut lengths according to the new parameter value
Implementation 3: only k is changed
:param x: the array of adjusted parameters
:return: F the sum((Jv - Jv_exp) ** 2.0)
"""
g_cut['tot'] = conductance.compute_k(g_cut['tot'] , k0 = x[0])
g_cut['tot'] = flux.flux(g_cut['tot'], Jv = _Jv[0], psi_e = psi_base + DP_cnf[0], psi_base = psi_base,
invert_model = True)
Jv = g_cut['tot'].property('J_out')[1]
F = (Jv - _Jv[0]) ** 2.0
count = 1
for cut_length in cut_n_flow_length:
# _g = g_cut[str(cut_length)].copy() # not necessary and time-consuming
for vid in g_cut[str(cut_length)].vertices_iter(g_cut['tot'].max_scale()):
g_cut[str(cut_length)].property('K')[vid] = g_cut['tot'].property('K')[vid]
g_cut[str(cut_length)].property('k')[vid] = g_cut['tot'].property('k')[vid]
for i in tip_id[str(cut_length)]:
v = g_cut['tot'].parent(i)
g_cut[str(cut_length)].property('k')[v] = g_cut[str(cut_length)].property('K')[v]
g_cut[str(cut_length)] = flux.flux(g_cut[str(cut_length)], Jv = _Jv[count], psi_e = psi_base + DP_cnf[count], psi_base = psi_base,
invert_model = True, cut_and_flow = True)
Jv = g_cut[str(cut_length)].property('J_out')[1]
F += (Jv - _Jv[count]) ** 2.0
count += 1
return F
# architecture with filename in aqua team format
if len(parameter.archi['input_file']) > 0:
if df_archi is None:
# archi_f = glob.glob(parameter.archi['input_dir'] + parameter.archi['input_file'][0])
# archi_f = archi_f[0]
fname = str(Path(parameter.archi['input_dir']) / parameter.archi['input_file'][0]) # deals with path in windows
archi_f = glob.glob(fname) # deals with wildcards ? deprecated ???
archi_f = archi_f[0]
df_archi = read_archi_data(archi_f) if parameter.archi['read_architecture'] else None
index = archi_f.replace(glob.glob(parameter.archi['input_dir'])[0],"")
else:
index = ''
# length law data: override if necessary
if df_law is not None:
parameter.archi['length_data'] = df_law
psi_e = parameter.exp['psi_e']
psi_base = parameter.exp['psi_base']
columns = ['plant', 'cut length (m)', 'max_length', 'k (10-9 m/s/MPa)', 'length (m)', 'surface (m2)',
'Jv (uL/s)', 'Jexp (uL/s)']
results = {}
for key in columns:
results[key] = []
# read the data measurements from data base
if df_exp is not None:
for key in df_exp['arch']:
if str(key).lower() in index.lower():
_list = df_exp[df_exp.arch == key].filter(regex = '^J').dropna(axis = 1).values.tolist()
# parameter.exp['Jv'] = _list[0][0] # basal output flux full root (uncut)
# _Jv = _list[0][1:] # basal output flux cut root
_Jv = _list[0]
_list = df_exp[df_exp.arch == key].filter(regex = '^lcut').dropna(axis = 1).values.tolist()
cut_n_flow_length = _list[0] # cut lengthes
_list = df_exp[df_exp.arch == key].filter(regex = '^dP').dropna(axis = 1).values.tolist()
# the pressure difference is usually constant but sometimes, due to flow meter saturation, it may change
# in that case a list of values is given
if len(_list[0]) != 0:
DP_cnf = _list[0]
else:
DP_cnf = []
DP_cnf.append(psi_e - psi_base) # for compatibility reason with first analysis on arabidopsis
if len(DP_cnf) < len(cut_n_flow_length)+1: # if constant we create the list with the constant value
for i in range(1, len(cut_n_flow_length) + 1): DP_cnf.append(DP_cnf[0])
else:
_Jv = [parameter.exp['Jv']]
cut_n_flow_length = []
DP_cnf = [psi_e - psi_base]
axfold = parameter.output['axfold'][0]
radfold = parameter.output['radfold'][0]
# g_cut['tot'], primary_length, _length, surface, seed = root_builder(df = df_archi, segment_length = parameter.archi['segment_length'],
# order_decrease_factor = parameter.archi['order_decrease_factor'], ref_radius = parameter.archi['ref_radius'])
g_cut['tot'], primary_length, _length, surface, seed = \
root_builder(df = df_archi,
primary_length = parameter.archi['primary_length'][0],
seed = parameter.archi['seed'][0],
delta = parameter.archi['branching_delay'][0],
nude_length = parameter.archi['nude_length'][0],
segment_length = parameter.archi['segment_length'],
length_data = parameter.archi['length_data'],
branching_variability = parameter.archi['branching_variability'],
order_max = parameter.archi['order_max'],
order_decrease_factor = parameter.archi['order_decrease_factor'],
ref_radius = parameter.archi['ref_radius'],
Flag_radius = Flag_radius)
g_cut['tot'], Keq, Jv = hydro_calculation(g_cut['tot'], psi_base = psi_base, psi_e = psi_base + DP_cnf[0])
###############################################################
#### WARNING : the mtg property 'position' must stay unchanged
#### because the axial conductivity is placed according to it
###############################################################
for cut_length in cut_n_flow_length:
tip_id[str(cut_length)] = \
flux.segments_at_length(g_cut['tot'], cut_length, dl = parameter.archi['segment_length'])
g_cut[str(cut_length)] = \
flux.cut_and_set_conductance(g_cut['tot'], cut_length, parameter.archi['segment_length'])
# g_cut[str(cut_length)], surface = radius.compute_surface(g_cut[str(cut_length)])
axial_data = list(conductance.axial(parameter.hydro['axial_conductance_data'], axfold))
###############################################################################################
## First adjustment: axfold, arfold that are coefficient factor of the radial conductivity k and
## and axial conductance K
###############################################################################################
if Flag_Optim_K:
# pre-optimization needed when there are several parameters so when K is optimized
if Flag_verbose: print("*********** pre_optimization ************************")
x = []
bnds = []
if Flag_Optim_K:
x.append(axfold)
bnds.append((1.0e-20, np.inf))
if Flag_Optim_k:
x.append(radfold)
bnds.append((1.0e-20, np.inf))
res = optimize.minimize(fun1, x, bounds = bnds, options = {'ftol': _tol})
if Flag_Optim_k:
radfold = res.x[-1] # always the last one even if the only one
if Flag_verbose: print('pre-optimization ar: {:0.2e}'.format(res.x[-1]))
if Flag_Optim_K:
axfold = res.x[0]
if Flag_verbose: print('pre-optimization ax: {:0.2e}'.format(res.x[0]))
if Flag_verbose: print("****************************************************************")
## update the conductivities according to the first adjustment
axial_data = list(conductance.axial(parameter.hydro['axial_conductance_data'], axfold))
k0 = parameter.hydro['k0'] *radfold
###############################################################################################
## 2d adjustment:
## -1 axial data adjusted
## -2 radial conductivit adjusted
## - 1 and 2 repeated until the k0 variation is below 0.1
###############################################################################################
x = []
x = axial_data[1]
bnds = []
n = len(x)
for i, val in enumerate(x):
bnds.append((1.0e-20, 1.0))
# linear constraints lb <= A.dot(x) <= ub
A = np.zeros((n, n))
lb = np.full(n, -np.inf)
ub = np.full(n, np.inf)
l = parameter.hydro['axial_conductance_data'][0]
a = dK_constraint # constraint on the 1st derivative
ni = n - 1
for i in range(ni): # downward derivative
A[i, i] = -1.
A[i, i + 1] = 1.
lb[i] = a * (l[i+1]-l[i])
i = ni
A[i, i-1] = -1.
A[i, i] = 1.
lb[i] = a * (l[i]-l[i-1])
k0_old = k0
F_old = (Jv - _Jv[0])**2.0
eps = 1e-9
F = 1.
if not (Flag_Optim_K or Flag_Optim_k):
k0_old2 = k0
else:
k0_old2 = k0 + 10
count2 = 0
while abs(k0-k0_old2) > dk_max:
k0_old2 = k0
# parameter.hydro['k0'] = k0
# #ajout Mistral
# _J_Sim=[]
# lcut = []
# lcut.append(0)
#
# g_cut['tot'] = conductance.compute_k(g_cut['tot'] , k0 = k0)
#
# g_cut['tot'] = flux.flux(g_cut['tot'], Jv = _Jv[0], psi_e = psi_base + DP_cnf[0], psi_base = psi_base,
# invert_model = True)
# _J_Sim.append(g_cut['tot'].property('J_out')[1])
#
# count = 1
# for cut_length in cut_n_flow_length:
# # _g = g_cut[str(cut_length)].copy() # not necessary and time-consuming
# lcut.append(g_cut['tot'].property('position')[1] - cut_length)
# for vid in g_cut[str(cut_length)].vertices_iter(g_cut['tot'].max_scale()):
# g_cut[str(cut_length)].property('K')[vid] = g_cut['tot'].property('K')[vid]
# g_cut[str(cut_length)].property('k')[vid] = g_cut['tot'].property('k')[vid]
#
# for i in tip_id[str(cut_length)]:
# v = g_cut['tot'].parent(i)
# g_cut[str(cut_length)].property('k')[v] = g_cut[str(cut_length)].property('K')[v]
#
# g_cut[str(cut_length)] = flux.flux(g_cut[str(cut_length)], Jv = _Jv[count], psi_e = psi_base + DP_cnf[count], psi_base = psi_base,
# invert_model = True, cut_and_flow = True)
# _J_Sim.append(g_cut[str(cut_length)].property('J_out')[1])
#
#
# count += 1
# plt.figure()
# plt.scatter(lcut,_Jv, c = 'black')
# plt.plot(lcut,_J_Sim, c = 'purple')
# plt.savefig(str(count2) + '.png')
# plt.close()
# count2 += 1
# # end ajout Mistral
## -1 axial data adjusted
#########################
constraints = optimize.LinearConstraint(A, lb, ub)
if Flag_Optim_K:
res = optimize.minimize(fun2, x, bounds = bnds, constraints = constraints, options={'ftol': _tol})
dKx = sum((x-res.x)**2.0)
axial_data[1] = list(res.x)
x = copy.deepcopy(res.x)
if Flag_verbose: print('finished minimize K: [', ', '.join('{:0.2e}'.format(i) for i in res.x), ']')
if Flag_Optim_k:
## -1 radial k adjusted
#######################
resk0 = optimize.minimize(fun3, k0, method = 'Nelder-Mead')
if Flag_verbose: print('finished minimize k0: , {:0.2e}'.format(resk0.x[0])) #,
# 'dk0**2.0 = {:0.2e}'.format((k0-resk0.x[0])**2.), 'dKx**2.0 = {:0.2e}'.format(dKx))
k0 = resk0.x[0]
else:
k0_old2 = k0
# print(resk0)
######################################
## Simulations with Kx and k adjusted
######################################
primary_length = g_cut['tot'].property('position')[1]
g_cut['tot'], Keq, Jv = hydro_calculation(g_cut['tot'], k_radial = k0 ,axial_data = axial_data, psi_base = psi_base,
psi_e = psi_base + DP_cnf[0])
# add some properties for display
g_cut['tot'].add_property('jsurf')
length = g_cut['tot'].property('length')
_radius = g_cut['tot'].property('radius')
j = g_cut['tot'].property('j')
jsurf = g_cut['tot'].property('jsurf')
for v in g_cut['tot'].vertices_iter(scale = 1):
jsurf[v] = j[v] / (_radius[v] * 2 * np.pi * length[v])
results['plant'].append(index)
results['max_length'].append(primary_length)
results['cut length (m)'].append(0.0)
results['k (10-9 m/s/MPa)'].append(k0)
results['length (m)'].append(_length)
results['surface (m2)'].append(surface)
results['Jv (uL/s)'].append(Jv)
results['Jexp (uL/s)'].append(_Jv[0])
count = 1
for cut_length in cut_n_flow_length:
g_cut[str(cut_length)] = g_cut[str(cut_length)].copy()
for vid in g_cut[str(cut_length)].vertices_iter(g_cut['tot'].max_scale()):
g_cut[str(cut_length)].property('K')[vid] = g_cut['tot'].property('K')[vid]
g_cut[str(cut_length)].property('k')[vid] = g_cut['tot'].property('k')[vid]
for i in tip_id[str(cut_length)]:
v = g_cut['tot'].parent(i)
g_cut[str(cut_length)].property('k')[v] = g_cut[str(cut_length)].property('K')[v]
g_cut[str(cut_length)] = flux.flux(g_cut[str(cut_length)], psi_e = psi_base + DP_cnf[count], psi_base = psi_base, invert_model = True)
Jv = g_cut[str(cut_length)].property('J_out')[1]
g_cut[str(cut_length)], surface = radius.compute_surface(g_cut[str(cut_length)])
_length = g_cut[str(cut_length)].nb_vertices(scale = 1) * parameter.archi['segment_length']
# add some properties for display
g_cut[str(cut_length)].add_property('jsurf')
length = g_cut[str(cut_length)].property('length')
_radius = g_cut[str(cut_length)].property('radius')
j = g_cut[str(cut_length)].property('j')
jsurf = g_cut[str(cut_length)].property('jsurf')
for v in g_cut[str(cut_length)].vertices_iter(scale = 1):
jsurf[v] = j[v] / (_radius[v] * 2 * np.pi * length[v])
primary_length = cut_length
results['plant'].append(index)
results['max_length'].append(primary_length)
results['cut length (m)'].append(g_cut['tot'].property('position')[1] - primary_length)
results['k (10-9 m/s/MPa)'].append(k0)
results['length (m)'].append(_length)
results['surface (m2)'].append(surface)
results['Jv (uL/s)'].append(Jv)
results['Jexp (uL/s)'].append(_Jv[count])
count += 1
dresults = pd.DataFrame(results, columns = columns)
optim_results = {}
optim_results['x'] = copy.deepcopy(parameter.hydro['axial_conductance_data'][0])
optim_results['K 1st'] = copy.deepcopy(parameter.hydro['axial_conductance_data'][1])
if Flag_Optim_K:
_x = list(res.x)
optim_results['K optimized'] = copy.deepcopy(_x)
else:
optim_results['K optimized'] = optim_results['K 1st']
doptim = pd.DataFrame(optim_results, columns = ['x', 'K 1st', 'K optimized'])
df = pd.concat([dresults, doptim], axis = 1)
if Flag_verbose:
pd.set_option('display.max_columns', None)
pd.set_option('display.expand_frame_repr', False)
print(dresults)
if Flag_Optim_K: print(doptim)
if output is not None: df.to_csv(output, index = False)
g_cut[0] = g_cut.pop('tot')
icut = 1
for cut_length in cut_n_flow_length:
g_cut[icut] = g_cut.pop(str(cut_length))
icut += 1
return df, g_cut
