+
+
+
+
+from datetime import datetime
+import cantera as ct
+import pandas as pd
+import numpy as np
+from scipy.optimize import minimize
+# noinspection PyPep8Naming
+import xml.etree.ElementTree as ET
+import seaborn as sns
+import matplotlib.pyplot as plt
+import shutil
+import copy
+from inspect import signature
+import os
+import re
+import pkg_resources
+
+sns.set()
+sns.set(font_scale=1.6)
+
+
+[docs]class REEPS:
+ r"""
+ Rare earth elements (REE or RE) Takes in experimental data
+ Returns parameters for GEM
+
+ .. note::
+
+ The order in which the REEs appear in the csv file must be the same
+ order as they appear in the xml, complex_names and
+ rare_earth_ion_names.
+
+ For example, say in exp_csv_filename's csv, RE_1 is Nd RE_2 is Pr,
+ and
+
+ .. code-block:: python
+
+ aq_solvent_name = 'H2O(L)'
+ extractant_name = '(HA)2(org)'
+ diluent_name = 'dodecane'
+
+ Then:
+
+ The csvs column ordering must be:
+
+ [h_i, h_eq, z_i, z_eq, Nd_aq_i, Nd_aq_eq, Nd_d_eq,
+ Pr_aq_i, Pr_aq_eq, Pr_d_eq]
+
+ The aqueous speciesArray must be
+ "H2O(L) H+ OH- Cl- Nd+++ Pr+++"
+
+ The organic speciesArray must be
+ "(HA)2(org) dodecane Nd(H(A)2)3(org) Pr(H(A)2)3(org)"
+
+ .. code-block:: python
+
+ complex_names = ['Nd(H(A)2)3(org)', 'Pr(H(A)2)3(org)']
+ rare_earth_ion_names = ['Nd+++', 'Pr+++']
+
+
+ :param exp_csv_filename: (str) csv file name with experimental data
+
+ In the .csv file, the rows are different experiments and
+ columns are the measured quantities.
+
+ The ordering of the columns needs to be:
+
+ [h_i, h_eq, z_i, z_eq, {RE_1}_aq_i, {RE_1}_aq_eq, {RE_1}_d_eq,
+ {RE_2}_aq_i, {RE_2}_aq_eq, {RE_2}_d_eq,...
+ {RE_N}_aq_i, {RE_N}_aq_eq, {RE_N}_d_eq]
+
+ Naming does not matter, just the order.
+
+ Where {RE_1}-{RE_N} are the rare earth element names of interest
+ i.e. Nd, Pr, La, etc.
+
+ Below is an explanation of the columns.
+
+ +-------+------------+------------------------------------------+
+ | Index | Column | Meaning |
+ +=======+============+==========================================+
+ | 0 | h_i | Initial Concentration of |
+ | | | H+ ions (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 1 | h_eq | Equilibrium concentration of |
+ | | | H+ ions (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 2 | z_i | Initial concentration of |
+ | | | extractant (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 3 | z_eq | Equilibrium concentration of |
+ | | | extractant (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 4 | {RE}_aq_i | Initial concentration of RE ions (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 5 | {RE}_aq_eq | Equilibrium concentration of RE ions |
+ | | | in aqueous phase (mol/L) |
+ +-------+------------+------------------------------------------+
+ | 6 | {RE}_d_eq | Equilibrium Ratio between amount of |
+ | | | RE atoms in organic to aqueous |
+ +-------+------------+------------------------------------------+
+ :param phases_xml_filename: (str) xml file with parameters
+ for equilibrium calc
+
+ Would recommend copying and modifying xmls located in data/xmls
+ or in Cantera's "data" folder
+
+ speciesArray fields need specific ordering.
+
+ In aqueous phase: aq_solvent_name, H+, OH-, Cl-, RE_1, RE_2, ..., RE_N
+
+ For aqueous phase, RE_1-RE_N represent RE ion names i.e. Nd+++, Pr+++
+
+ In organic phase : extractant_name, diluant_name, RE_1, RE_2, ..., RE_N
+
+ For organic phase, RE_1-RE_N represent RE complex names
+ i.e. Nd(H(A)2)3(org), Pr(H(A)2)3(org)
+
+ :param phase_names: (list) names of phases in xml file
+
+ Found in the xml file under <phase ... id={phase_name}>
+
+ :param aq_solvent_name: (str) name of aqueous solvent in xml file
+ :param extractant_name: (str) name of extractant in xml file
+ :param diluant_name: (str) name of diluant in xml file
+ :param complex_names: (list) names of complexes in xml file.
+
+ Ensure the ordering is correct
+ :param rare_earth_ion_names: (list) names of rare earth ions in xml file
+
+ Ensure the ordering is correct
+ :param re_species_list: (list) names of rare earth elements.
+
+ If ``None``, re_species_list will be rare_earth_ion_names without '+'
+ i.e. 'Nd+++'->'Nd'
+
+ Ensure the ordering is correct
+ :param aq_solvent_rho: (float) density of solvent (g/L)
+
+ If ``None``, molar volume/molecular weight is used from xml
+ :param extractant_rho: (float) density of extractant (g/L)
+
+ If ``None``, molar volume/molecular weight is used from xml
+ :param diluant_rho: (float) density of diluant (g/L)
+
+ If ``None``, molar volume/molecular weight is used from xml
+ :param opt_dict: (dict) dictionary containing info about which
+ species parameters are updated to fit model to experimental data
+
+ Should have the format as below
+
+ .. code-block:: python
+
+ opt_dict = {"species1":
+ {"parameter1": "guess1",
+ "parameter2": "guess2",
+ ...
+ "parameterN": "guessN"}
+ "species2":
+ {"parameter1": "guess1",
+ "parameter2": "guess2",
+ ...
+ "parameterN": "guessN"}
+ ...
+ "speciesN":
+ {"parameter1": "guess1",
+ "parameter2": "guess2",
+ ...
+ "parameterN": "guessN"}
+ }
+ :param objective_function: (function or str) function to compute objective
+
+ By default, the objective function is log mean squared error
+ of distribution ratio
+
+ .. code-block:: python
+
+ np.sum((np.log10(d_pred)-np.log10(d_meas))^2)
+
+ Function needs to take inputs:
+
+ .. code-block:: python
+
+ objective_function(predicted_dict, measured_df, kwargs)
+
+ ``kwargs`` is optional
+
+ Function needs to return: (float) value computed by objective function
+
+ Below is the guide for referencing predicted values
+
+ +---------------------------+--------------------------------+
+ | To access | Use |
+ +===========================+================================+
+ | hydrogen ion conc in aq | predicted_dict['h_eq'] |
+ +---------------------------+--------------------------------+
+ | extractant conc in org | predicted_dict['z_eq'] |
+ +---------------------------+--------------------------------+
+ | RE ion eq conc in aq | predicted_dict['{RE}_aq_eq'] |
+ +---------------------------+--------------------------------+
+ | RE complex eq conc in org | predicted_dict['{RE}_org_eq'] |
+ +---------------------------+--------------------------------+
+ | RE distribution ratio | predicted_dict['{RE}_d_eq'] |
+ +---------------------------+--------------------------------+
+
+ Replace "{RE}" with rare earth element i.e. Nd, La, etc.
+
+ For measured values, use the same names, but
+ replace ``predicted_dict`` with ``measured_df``
+ :param optimizer: (function or str) function to perform optimization
+
+ .. note::
+
+ The optimized variables are not directly the species parameters,
+ but instead are first multiplied by the initial guess before
+ sending becoming the species parameters.
+
+ For example, say
+
+ .. code-block:: python
+
+ opt_dict = {'Nd(H(A)2)3(org):'h0':-4.7e6}
+
+ If the bounds on h0 need to be [-4.7e7,-4.7e5], then
+ divide the bounds by the guess and get
+
+ .. code-block:: python
+
+ "bounds": [(1e-1, 1e1)]
+
+ Though fit() returns a structure identical to opt_dict with correct
+ scaled values, in case bounds and constraints are used, you must
+ note the optimized x's are first multiplied by the initial guess
+ before written to the xml.
+
+
+ By default, the optimizer is scipy's optimize function with
+
+ .. code-block:: python
+
+ default_kwargs= {"method": 'SLSQP',
+ "bounds": [(1e-1, 1e1)] * len(x_guess),
+ "constraints": (),
+ "options": {'disp': True,
+ 'maxiter': 1000,
+ 'ftol': 1e-6}}
+
+ Function needs to take inputs:
+ ``optimizer(objective_function, x_guess, kwargs)``
+
+ ``kwargs`` is optional
+
+ Function needs to return: (np.ndarray) Optimized parameters
+
+ :param temp_xml_file_path: (str) path to temporary xml file.
+
+ This xml file is a duplicate of the phases_xml_file name and is
+ modified during the optimization process to avoid changing the original
+ xml file
+
+ default is local temp folder
+ """
+
+ def __init__(self,
+ exp_csv_filename,
+ phases_xml_filename,
+ phase_names,
+ aq_solvent_name,
+ extractant_name,
+ diluant_name,
+ complex_names,
+ rare_earth_ion_names,
+ re_species_list=None,
+ aq_solvent_rho=None,
+ extractant_rho=None,
+ diluant_rho=None,
+ opt_dict=None,
+ objective_function='Log-MSE',
+ optimizer='SLSQP',
+ temp_xml_file_path=None
+ ):
+ self._built_in_obj_list = ['Log-MSE']
+ self._built_in_opt_list = ['SLSQP']
+ self._exp_csv_filename = exp_csv_filename
+ self._phases_xml_filename = phases_xml_filename
+ self._opt_dict = opt_dict
+ self._phase_names = phase_names
+ self._aq_solvent_name = aq_solvent_name
+ self._extractant_name = extractant_name
+ self._diluant_name = diluant_name
+ self._complex_names = complex_names
+ self._rare_earth_ion_names = rare_earth_ion_names
+ self._aq_solvent_rho = aq_solvent_rho
+ self._extractant_rho = extractant_rho
+ self._diluant_rho = diluant_rho
+ self._objective_function = None
+ self.set_objective_function(objective_function)
+ self._optimizer = None
+ self._re_species_list = re_species_list
+ self.set_optimizer(optimizer)
+ if temp_xml_file_path is None:
+ temp_xml_file_path = r'{0}/temp.xml'.format(os.getenv('TEMP'))
+ self._temp_xml_file_path = temp_xml_file_path
+ # Try and except for adding package data to path.
+ # This only works for sdist, not bdist
+ # If bdist is needed, research "manifest.in" python setup files
+ try:
+ shutil.copyfile(self._phases_xml_filename,
+ self._temp_xml_file_path)
+ self._phases = ct.import_phases(self._phases_xml_filename,
+ phase_names)
+ except FileNotFoundError:
+ self._phases_xml_filename = \
+ pkg_resources.resource_filename('reeps',
+ r'..\data\xmls\{0}'.format(
+ phases_xml_filename))
+ shutil.copyfile(self._phases_xml_filename,
+ self._temp_xml_file_path)
+ self._phases = ct.import_phases(self._phases_xml_filename,
+ phase_names)
+ try:
+ self._exp_df = pd.read_csv(self._exp_csv_filename)
+ except FileNotFoundError:
+ self._exp_csv_filename = pkg_resources.resource_filename(
+ 'reeps', r'..\data\csvs\{0}'.format(self._exp_csv_filename))
+ self._exp_df = pd.read_csv(self._exp_csv_filename)
+
+ self._exp_df_columns = ['h_i', 'h_eq', 'z_i', 'z_eq']
+ if self._re_species_list is None:
+ self._re_species_list = []
+ for name in self._rare_earth_ion_names:
+ species = name.replace('+', '')
+ self._re_species_list.append(species)
+ for species in self._re_species_list:
+ self._exp_df_columns.append('{0}_aq_i'.format(species))
+ self._exp_df_columns.append('{0}_aq_eq'.format(species))
+ self._exp_df_columns.append('{0}_d_eq'.format(species))
+ self._exp_df.columns = self._exp_df_columns
+
+ self._in_moles = None
+
+ self._aq_ind = None
+ self._org_ind = None
+ self._re_charges = None
+
+ self.set_in_moles(feed_vol=1)
+ self._predicted_dict = None
+ self.update_predicted_dict()
+
+[docs] @staticmethod
+ def slsqp_optimizer(objective, x_guess):
+ """ The default optimizer for REEPS
+
+ Uses scipy.minimize with options
+
+ .. code-block:: python
+
+ default_kwargs= {"method": 'SLSQP',
+ "bounds": [(1e-1, 1e1)*len(x_guess)],
+ "constraints": (),
+ "options": {'disp': True,
+ 'maxiter': 1000,
+ 'ftol': 1e-6}}
+
+ :param objective: (func) the objective function
+ :param x_guess: (np.ndarray) the initial guess (always 1)
+ :returns: (np.ndarray) Optimized parameters
+ """
+ optimizer_kwargs = {"method": 'SLSQP',
+ "bounds": [(1e-1, 1e1)] * len(x_guess),
+ "constraints": (),
+ "options": {'disp': True,
+ 'maxiter': 1000,
+ 'ftol': 1e-6}}
+ res = minimize(objective, x_guess, **optimizer_kwargs)
+ est_parameters = res.x
+ return est_parameters
+
+[docs] def log_mean_squared_error(self, predicted_dict, meas_df):
+ """Default objective function for REEPS
+
+ Returns the log mean squared error of
+ predicted distribution ratios (d=n_org/n_aq)
+ to measured d.
+
+ np.sum((np.log10(d_pred)-np.log10(d_meas))\**2)
+
+ :param predicted_dict: (dict) contains predicted data
+ :param meas_df: (pd.DataFrame) contains experimental data
+ :return: (float) log mean squared error between predicted and measured
+ """
+ meas = np.concatenate([meas_df['{0}_d_eq'.format(species)].values
+ for species in self._re_species_list])
+ pred = np.concatenate([
+ predicted_dict['{0}_d_eq'.format(species)]
+ for species in self._re_species_list])
+ log_pred = np.log10(pred)
+ log_meas = np.log10(meas)
+ obj = np.sum((log_pred - log_meas) ** 2)
+ return obj
+
+[docs] def get_exp_df(self) -> pd.DataFrame:
+ """Returns the experimental DataFrame
+
+ :return: (pd.DataFrame) Experimental data
+ """
+ return self._exp_df
+
+[docs] def set_exp_df(self, exp_csv_filename):
+ """Changes the experimental DataFrame to input exp_csv_filename data
+ and renames columns to internal REEPS names
+
+
+ h_i, h_eq, z_i, z_eq, {RE}_aq_i, {RE}_aq_eq, {RE}_d
+
+ See class docstring on "exp_csv_filename" for further explanations.
+
+ :param exp_csv_filename: (str) file name/path for experimental data csv
+ """
+ self._exp_csv_filename = exp_csv_filename
+ try:
+ self._exp_df = pd.read_csv(self._exp_csv_filename)
+ except FileNotFoundError:
+ self._exp_csv_filename = pkg_resources.resource_filename(
+ 'reeps', r'..\data\csvs\{0}'.format(self._exp_csv_filename))
+ self._exp_df = pd.read_csv(self._exp_csv_filename)
+ self._exp_df_columns = ['h_i', 'h_eq', 'z_i', 'z_eq']
+ if self._re_species_list is None:
+ self._re_species_list = []
+ for name in self._rare_earth_ion_names:
+ species = name.replace('+', '')
+ self._re_species_list.append(species)
+ for species in self._re_species_list:
+ self._exp_df_columns.append('{0}_aq_i'.format(species))
+ self._exp_df_columns.append('{0}_aq_eq'.format(species))
+ self._exp_df_columns.append('{0}_d_eq'.format(species))
+ self._exp_df.columns = self._exp_df_columns
+ self.set_in_moles(feed_vol=1)
+ self.update_predicted_dict()
+ return None
+
+[docs] def get_phases(self) -> list:
+ """
+ Returns the list of Cantera solutions
+
+ :return: (list) list of Cantera solutions/phases
+ """
+ return self._phases
+
+[docs] def set_phases(self, phases_xml_filename, phase_names):
+ """Change list of Cantera solutions by inputting
+ new xml file name and phase names
+
+ Also runs set_in_moles to set initial molality to 1 g/L
+
+ :param phases_xml_filename: (str) xml file with parameters
+ for equilibrium calc
+ :param phase_names: (list) names of phases in xml file
+ """
+ self._phases_xml_filename = phases_xml_filename
+ self._phase_names = phase_names
+ # Try and except for adding package data to path.
+ # This only works for sdist, not bdist
+ # If bdist is needed, research "manifest.in" python setup files
+ try:
+ shutil.copyfile(self._phases_xml_filename,
+ self._temp_xml_file_path)
+ self._phases = ct.import_phases(self._phases_xml_filename,
+ phase_names)
+ except FileNotFoundError:
+ self._phases_xml_filename = \
+ pkg_resources.resource_filename('reeps',
+ r'..\data\xmls\{0}'.format(
+ phases_xml_filename))
+ shutil.copyfile(self._phases_xml_filename,
+ self._temp_xml_file_path)
+ self._phases = ct.import_phases(self._phases_xml_filename,
+ phase_names)
+ self.set_in_moles(feed_vol=1)
+ self.update_predicted_dict()
+ return None
+
+[docs] def get_opt_dict(self) -> dict:
+ """
+ Returns the dictionary containing optimization information
+
+ :return: (dict) dictionary containing info about which
+ species parameters are updated to fit model to experimental data
+ """
+ return self._opt_dict
+
+[docs] def set_opt_dict(self, opt_dict):
+ """
+ Change the dictionary to input opt_dict.
+
+ opt_dict specifies species parameters to be updated to
+ fit model to data
+
+ See class docstring on "opt_dict" for more information.
+
+ :param opt_dict: (dict) dictionary containing info about which
+ species parameters are updated to fit model to experimental data
+ """
+
+ self._opt_dict = opt_dict
+ return None
+
+[docs] def get_aq_solvent_name(self) -> str:
+ """Returns aq_solvent_name
+
+ :return: aq_solvent_name: (str) name of aqueous solvent in xml file
+ """
+ return self._aq_solvent_name
+
+[docs] def set_aq_solvent_name(self, aq_solvent_name):
+ """ Change aq_solvent_name to input aq_solvent_name
+
+ :param aq_solvent_name: (str) name of aqueous solvent in xml file
+ """
+ self._aq_solvent_name = aq_solvent_name
+ return None
+
+[docs] def get_extractant_name(self) -> str:
+ """Returns extractant name
+
+ :return: extractant_name: (str) name of extractant in xml file
+ """
+ return self._extractant_name
+
+[docs] def set_extractant_name(self, extractant_name):
+ """
+ Change extractant_name to input extractant_name
+ :param extractant_name: (str) name of extractant in xml file
+ """
+ self._extractant_name = extractant_name
+ return None
+
+[docs] def get_diluant_name(self) -> str:
+ """ Returns diluant name
+ :return: diluant_name: (str) name of diluant in xml file
+ """
+ return self._diluant_name
+
+[docs] def set_diluant_name(self, diluant_name):
+ """
+ Change diluant_name to input diluant_name
+
+ :param diluant_name: (str) name of diluant in xml file
+ """
+ self._diluant_name = diluant_name
+ return None
+
+[docs] def get_complex_names(self) -> list:
+ """Returns list of complex names
+
+ :return: complex_names: (list) names of complexes in xml file.
+ """
+ return self._complex_names
+
+[docs] def set_complex_names(self, complex_names):
+ """Change complex names list to input complex_names
+
+ :param complex_names: (list) names of complexes in xml file.
+ """
+ self._complex_names = complex_names
+ return None
+
+[docs] def get_rare_earth_ion_names(self) -> list:
+ """Returns list of rare earth ion names
+
+ :return: rare_earth_ion_names: (list) names of rare earth ions in
+ xml file
+ """
+ return self._rare_earth_ion_names
+
+[docs] def set_rare_earth_ion_names(self, rare_earth_ion_names):
+ """Change list of rare earth ion names to input
+ rare_earth_ion_names
+
+ :param rare_earth_ion_names: (list) names of rare earth ions in
+ xml file
+ """
+ self._rare_earth_ion_names = rare_earth_ion_names
+ return None
+
+[docs] def get_re_species_list(self) -> list:
+ """Returns list of rare earth element names
+
+ :return: re_species_list: (list) names of rare earth elements in
+ xml file
+ """
+ return self._re_species_list
+
+[docs] def set_re_species_list(self, re_species_list):
+ """Change list of rare earth ion names to input
+ rare_earth_ion_names
+
+ :param re_species_list: (list) names of rare earth elements in
+ xml file
+ """
+ self._re_species_list = re_species_list
+ return None
+
+[docs] def get_aq_solvent_rho(self) -> str:
+ """Returns aqueous solvent density (g/L)
+
+ :return: aq_solvent_rho: (float) density of aqueous solvent
+ """
+ return self._aq_solvent_rho
+
+[docs] def set_aq_solvent_rho(self, aq_solvent_rho):
+ """Changes aqueous solvent density (g/L) to input aq_solvent_rho
+
+ :param aq_solvent_rho: (float) density of aqueous solvent
+ """
+ self._aq_solvent_rho = aq_solvent_rho
+ return None
+
+[docs] def get_extractant_rho(self) -> str:
+ """Returns extractant density (g/L)
+
+ :return: extractant_rho: (float) density of extractant
+ """
+ return self._extractant_rho
+
+[docs] def set_extractant_rho(self, extractant_rho):
+ """Changes extractant density (g/L) to input extractant_rho
+
+ :param extractant_rho: (float) density of extractant
+ """
+ self._extractant_rho = extractant_rho
+ return None
+
+[docs] def get_diluant_rho(self) -> str:
+ """Returns diluant density (g/L)
+
+ :return: diluant_rho: (float) density of diluant
+ """
+ return self._diluant_rho
+
+[docs] def set_diluant_rho(self, diluant_rho):
+ """Changes diluant density (g/L) to input diluant_rho
+
+ :param diluant_rho: (float) density of diluant
+ """
+ self._diluant_rho = diluant_rho
+ return None
+
+[docs] def set_in_moles(self, feed_vol):
+ """Function that initializes mole fractions to input feed_vol
+
+ This function is called at initialization
+
+ Sets in_moles to a pd.DataFrame containing initial mole fractions
+
+ Columns for species and rows for different experiments
+
+ This function also calls update_predicted_dict
+
+ :param feed_vol: (float) feed volume of mixture (g/L)
+ """
+ phases_copy = self._phases.copy()
+ exp_df = self._exp_df.copy()
+ solvent_name = self._aq_solvent_name
+ extractant_name = self._extractant_name
+ diluant_name = self._diluant_name
+ solvent_rho = self._aq_solvent_rho
+ extractant_rho = self._extractant_rho
+ diluant_rho = self._diluant_rho
+ re_names = self._rare_earth_ion_names
+ re_species_list = self._re_species_list
+
+ mixed = ct.Mixture(phases_copy)
+ aq_ind = None
+ solvent_ind = None
+ for ind, phase in enumerate(phases_copy):
+ if solvent_name in phase.species_names:
+ aq_ind = ind
+ solvent_ind = phase.species_names.index(solvent_name)
+ if aq_ind is None:
+ raise Exception('Solvent "{0}" not found \
+ in xml file'.format(solvent_name))
+
+ if aq_ind == 0:
+ org_ind = 1
+ else:
+ org_ind = 0
+ self._aq_ind = aq_ind
+ self._org_ind = org_ind
+ extractant_ind = phases_copy[org_ind].species_names.index(
+ extractant_name)
+ diluant_ind = phases_copy[org_ind].species_names.index(diluant_name)
+
+ re_ind_list = [phases_copy[aq_ind].species_names.index(re_name)
+ for re_name in re_names]
+ re_charges = np.array([phases_copy[aq_ind].species(re_ind).charge
+ for re_ind in re_ind_list])
+ self._re_charges = re_charges
+
+ mix_aq = mixed.phase(aq_ind)
+ mix_org = mixed.phase(org_ind)
+ solvent_mw = mix_aq.molecular_weights[solvent_ind] # g/mol
+ extractant_mw = mix_org.molecular_weights[extractant_ind]
+ diluant_mw = mix_org.molecular_weights[diluant_ind]
+ if solvent_rho is None:
+ solvent_rho = mix_aq(aq_ind).partial_molar_volumes[
+ solvent_ind] / solvent_mw * 1e6 # g/L
+ self._aq_solvent_rho = solvent_rho
+ if extractant_rho is None:
+ extractant_rho = mix_org(org_ind).partial_molar_volumes[
+ extractant_ind] / extractant_mw * 1e6
+ self._extractant_rho = extractant_rho
+ if diluant_rho is None:
+ diluant_rho = mix_org(org_ind).partial_molar_volumes[
+ extractant_ind] / extractant_mw * 1e6
+ self._diluant_rho = diluant_rho
+
+ in_moles_data = []
+ aq_phase_solvent_moles = feed_vol * solvent_rho / solvent_mw
+ for index, row in exp_df.iterrows():
+ h_plus_moles = feed_vol * row['h_i']
+ hydroxide_ions = 0
+ rare_earth_moles = np.array([feed_vol * row[
+ '{0}_aq_i'.format(re_species)]
+ for re_species in re_species_list])
+ re_charge_sum = np.sum(re_charges * rare_earth_moles)
+ chlorine_moles = re_charge_sum + h_plus_moles
+ extractant_moles = feed_vol * row['z_i']
+ extractant_vol = extractant_moles * extractant_mw / extractant_rho
+ diluant_vol = feed_vol - extractant_vol
+ diluant_moles = diluant_vol * diluant_rho / diluant_mw
+ complex_moles = np.zeros(len(re_species_list))
+
+ species_moles_aq = [aq_phase_solvent_moles,
+ h_plus_moles,
+ hydroxide_ions,
+ chlorine_moles]
+ species_moles_aq.extend(list(rare_earth_moles))
+ species_moles_org = [extractant_moles, diluant_moles]
+ species_moles_org.extend(list(complex_moles))
+ if aq_ind == 0:
+ species_moles = species_moles_aq + species_moles_org
+ else:
+ species_moles = species_moles_org + species_moles_aq
+ in_moles_data.append(species_moles)
+ self._in_moles = pd.DataFrame(
+ in_moles_data, columns=mixed.species_names)
+ self.update_predicted_dict()
+ return None
+
+[docs] def get_in_moles(self) -> pd.DataFrame:
+ """Returns the in_moles DataFrame which contains the initial mole
+ fractions of each species for each experiment
+
+ :return: in_moles: (pd.DataFrame) DataFrame with initial mole fractions
+ """
+ return self._in_moles
+
+[docs] def set_objective_function(self, objective_function):
+ """Change objective function to input objective_function.
+
+ See class docstring on "objective_function" for instructions
+
+ :param objective_function: (func) Objective function to quantify
+ error between model and experimental data
+ """
+ if not callable(objective_function) \
+ and objective_function not in self._built_in_obj_list:
+ raise Exception(
+ "objective_function must be a function "
+ "or in this strings list: {0}".format(
+ self._built_in_obj_list))
+ if callable(objective_function):
+ if len(signature(objective_function).parameters) < 2:
+ raise Exception(
+ "objective_function must be a function "
+ "with at least 3 arguments:"
+ " f(predicted_dict, experimental_df, kwargs)")
+ if objective_function == 'Log-MSE':
+ objective_function = self.log_mean_squared_error
+ self._objective_function = objective_function
+ return None
+
+[docs] def get_objective_function(self):
+ """Returns objective function
+
+ :return: objective_function: (func) Objective function to quantify
+ error between model and experimental data
+ """
+ return self._objective_function
+
+[docs] def set_optimizer(self, optimizer):
+ """Change optimizer function to input optimizer.
+
+ See class docstring on "optimizer" for instructions
+
+ :param optimizer: (func) Optimizer function to minimize objective
+ function
+ """
+ if not callable(optimizer) \
+ and optimizer not in self._built_in_opt_list:
+ raise Exception(
+ "optimizer must be a function "
+ "or in this strings list: {0}".format(
+ self._built_in_opt_list))
+ if callable(optimizer):
+ if len(signature(optimizer).parameters) < 2:
+ raise Exception(
+ "optimizer must be a function "
+ "with at least 2 arguments: "
+ "f(objective_func,x_guess, kwargs)")
+ if optimizer == 'SLSQP':
+ optimizer = self.slsqp_optimizer
+ self._optimizer = optimizer
+ return None
+
+[docs] def get_optimizer(self):
+ """Returns objective function
+
+ :return: optimizer: (func) Optimizer function to minimize objective
+ function
+ """
+ return self._optimizer
+
+[docs] def get_temp_xml_file_path(self):
+ """Returns path to temporary xml file.
+
+ This xml file is a duplicate of the phases_xml_file name and is
+ modified during the optimization process to avoid changing the original
+ xml file.
+
+ :return: temp_xml_file_path: (str) path to temporary xml file.
+ """
+ return self._temp_xml_file_path
+
+[docs] def set_temp_xml_file_path(self, temp_xml_file_path):
+ """Changes temporary xml file path to input temp_xml_file_path.
+
+ This xml file is a duplicate of the phases_xml_file name and is
+ modified during the optimization process to avoid changing the original
+ xml file.
+
+ :param temp_xml_file_path: (str) path to temporary xml file.
+ """
+ self._temp_xml_file_path = temp_xml_file_path
+ return None
+
+[docs] def update_predicted_dict(self,
+ phases_xml_filename=None,
+ phase_names=None):
+ """Function that computes the predicted equilibrium concentrations
+ the fed phases_xml_filename parameters predicts given the initial
+ mole fractions set by in_moles()
+
+ :param phases_xml_filename: (str)xml file with parameters
+ for equilibrium calc. If ``None``, the
+ current phases_xml_filename is used.
+ :param phase_names: (list) names of phases in xml file.
+ If ``None``, the current phases_names is used.
+ """
+ if phases_xml_filename is None:
+ phases_xml_filename = self._phases_xml_filename
+ if phase_names is None:
+ phase_names = self._phase_names
+ aq_ind = self._aq_ind
+ org_ind = self._org_ind
+ complex_names = self._complex_names
+ extractant_name = self._extractant_name
+ rare_earth_ion_names = self._rare_earth_ion_names
+ in_moles = self._in_moles
+ re_species_list = self._re_species_list
+
+ phases_copy = ct.import_phases(phases_xml_filename, phase_names)
+ mix = ct.Mixture(phases_copy)
+ key_names = ['h_eq', 'z_eq']
+ for re_species in re_species_list:
+ key_names.append('{0}_aq_eq'.format(re_species))
+ key_names.append('{0}_org_eq'.format(re_species))
+ key_names.append('{0}_d_eq'.format(re_species))
+
+ predicted_dict = {'{0}'.format(key_name): []
+ for key_name in key_names}
+
+ for row in in_moles.values:
+ mix.species_moles = row
+ mix.equilibrate('TP', log_level=0)
+ re_org_array = np.array([mix.species_moles[mix.species_index(
+ org_ind, complex_name)] for complex_name in complex_names])
+ re_aq_array = np.array([mix.species_moles[mix.species_index(
+ aq_ind, re_ion_name)] for re_ion_name in rare_earth_ion_names])
+ d_array = re_org_array / re_aq_array
+ hydrogen_ions = mix.species_moles[mix.species_index(aq_ind, 'H+')]
+ extractant = mix.species_moles[mix.species_index(
+ org_ind, extractant_name)]
+ for index, re_species in enumerate(re_species_list):
+ predicted_dict['{0}_aq_eq'.format(
+ re_species)].append(re_aq_array[index])
+ predicted_dict['{0}_org_eq'.format(
+ re_species)].append(re_org_array[index])
+ predicted_dict['{0}_d_eq'.format(
+ re_species)].append(d_array[index])
+ predicted_dict['h_eq'].append(hydrogen_ions)
+ predicted_dict['z_eq'].append(extractant)
+ for key, value in predicted_dict.items():
+ predicted_dict[key] = np.array(value)
+ self._predicted_dict = predicted_dict
+ return None
+
+[docs] def get_predicted_dict(self):
+ """Returns predicted dictionary of species concentrations
+ that xml parameters predicts given current in_moles
+
+ :return: predicted_dict: (dict) dictionary of species concentrations
+ """
+ return self._predicted_dict
+
+ def _internal_objective(self, x, kwargs=None):
+ """
+ Internal objective function. Uses objective function to compute value
+ If the optimizer requires vectorized variables ie pso, this function
+ takes care of it
+
+ :param x: (list) thermo properties varied to minimize objective func
+ :param kwargs: (list) arguments for objective_function
+ """
+ temp_xml_file_path = self._temp_xml_file_path
+ exp_df = self._exp_df
+ objective_function = self._objective_function
+ opt_dict = copy.deepcopy(self._opt_dict)
+ x = np.array(x)
+
+ if len(x.shape) == 1:
+ xs = np.array([x])
+ vectorized_x = False
+ else:
+ vectorized_x = True
+ xs = x
+ objective_values = []
+ for x in xs:
+ i = 0
+ for species_name in opt_dict.keys():
+ for thermo_prop in opt_dict[species_name].keys():
+ opt_dict[species_name][thermo_prop] *= x[i]
+ i += 1
+ self.update_xml(opt_dict, temp_xml_file_path)
+
+ self.update_predicted_dict(temp_xml_file_path)
+ predicted_dict = self.get_predicted_dict()
+ self.update_predicted_dict()
+
+ if kwargs is None:
+ # noinspection PyCallingNonCallable
+ obj = objective_function(predicted_dict, exp_df)
+ else:
+ # noinspection PyCallingNonCallable
+ obj = objective_function(predicted_dict, exp_df, **kwargs)
+ objective_values.append(obj)
+ if vectorized_x:
+ objective_values = np.array(objective_values)
+ else:
+ objective_values = objective_values[0]
+ return objective_values
+
+[docs] def fit(self,
+ objective_function=None,
+ optimizer=None,
+ objective_kwargs=None,
+ optimizer_kwargs=None) -> dict:
+ """Fits experimental to modeled data by minimizing objective function
+ with optimizer. Returns dictionary with opt_dict structure
+
+ :param objective_function: (function) function to compute objective
+ :param optimizer: (function) function to perform optimization
+ :param optimizer_kwargs: (dict) arguments for optimizer
+ :param objective_kwargs: (dict) arguments for objective function
+ :returns opt_dict: (dict) optimized opt_dict. Has identical structure
+ as opt_dict
+ """
+ if objective_function is not None:
+ self.set_objective_function(objective_function)
+ if optimizer is not None:
+ self.set_optimizer(optimizer)
+
+ def objective(x):
+ return self._internal_objective(x, objective_kwargs)
+
+ optimizer = self._optimizer
+ opt_dict = copy.deepcopy(self._opt_dict)
+ i = 0
+ for species_name in opt_dict.keys():
+ for _ in opt_dict[species_name].keys():
+ i += 1
+ x_guess = np.ones(i)
+
+ if optimizer_kwargs is None:
+ # noinspection PyCallingNonCallable
+ est_parameters = optimizer(objective, x_guess)
+ else:
+ # noinspection PyCallingNonCallable
+ est_parameters = optimizer(objective, x_guess, **optimizer_kwargs)
+
+ i = 0
+ for species_name in opt_dict.keys():
+ for thermo_prop in opt_dict[species_name].keys():
+ opt_dict[species_name][thermo_prop] *= est_parameters[i]
+ i += 1
+ return opt_dict
+
+[docs] def update_xml(self,
+ info_dict,
+ phases_xml_filename=None):
+ """updates xml file with info_dict
+
+ :param info_dict: (dict) info in {species_names:{thermo_prop:val}}
+ Requires an identical structure to opt_dict
+ :param phases_xml_filename: (str) xml filename if editing other xml
+ If ``None``, the current xml will be modified and the internal
+ Cantera phases will be refreshed to the new values.
+ """
+ if phases_xml_filename is None:
+ phases_xml_filename = self._phases_xml_filename
+
+ tree = ET.parse(phases_xml_filename)
+ root = tree.getroot()
+ # Update xml file
+ for species_name in info_dict.keys():
+ for thermo_prop in info_dict[species_name].keys():
+ for species in root.iter('species'):
+ if species.attrib['name'] == species_name:
+ for changed_prop in species.iter(thermo_prop):
+ changed_prop.text = str(
+ info_dict[species_name][thermo_prop])
+ now = datetime.now()
+ changed_prop.set('updated',
+ 'Updated at {0}:{1} {2}-{3}-{4}'
+ .format(now.hour, now.minute,
+ now.month, now.day,
+ now.year))
+
+ tree.write(phases_xml_filename)
+ if phases_xml_filename == self._phases_xml_filename:
+ self.set_phases(self._phases_xml_filename, self._phase_names)
+ return None
+
+[docs] def parity_plot(self,
+ compared_value=None,
+ save_path=None,
+ print_r_squared=False):
+ """
+ Parity plot between measured and predicted compared_value.
+ Default compared value is {RE_1}_aq_eq
+
+ :param compared_value: (str) Quantity to compare predicted and
+ experimental data. Can be any column containing "eq" in exp_df i.e.
+ h_eq, z_eq, {RE}_d_eq, etc.
+ :param save_path: (str) save path for parity plot
+ :param print_r_squared: (boolean) To plot or not to plot r-squared
+ value. Prints 2 places past decimal
+ """
+ exp_df = self.get_exp_df()
+ predicted_dict = self.get_predicted_dict()
+ re_species_list = self._re_species_list
+ extractant_name = self.get_extractant_name()
+ re_charges = self._re_charges
+ if compared_value is None:
+ compared_value = '{0}_aq_eq'.format(re_species_list[0])
+ pred = predicted_dict[compared_value]
+ meas = exp_df[compared_value].values
+ min_data = np.min([pred, meas])
+ max_data = np.max([pred, meas])
+ min_max_data = np.array([min_data, max_data])
+ name_breakdown = re.findall('[^_\W]+', compared_value)
+ compared_species = name_breakdown[0]
+ if compared_species == 'h':
+ species_name = '$H^+$'
+ elif compared_species == 'z':
+ species_name = extractant_name
+ else:
+ phase = name_breakdown[1]
+ if phase == 'aq':
+ re_charge = re_charges[re_species_list.index(compared_species)]
+ species_name = '$%s^{%d+}$' % (compared_species, re_charge)
+ elif phase == 'd':
+ species_name = '{0} distribution ratio'.format(
+ compared_species)
+ else:
+ species_name = '{0} complex'.format(compared_species)
+ fig, ax = plt.subplots()
+ p1 = sns.scatterplot(meas, pred, color="r",
+ label="{0} eq. conc. (mol/L)".format(
+ species_name),
+ legend=False)
+ sns.lineplot(min_max_data, min_max_data, color="b", label="")
+ if print_r_squared:
+ p1.text(min_max_data[0],
+ min_max_data[1] * 0.9,
+ '$R^2$={0:.2f}'.format(self.r_squared(compared_value)))
+ plt.legend(loc='lower right')
+ else:
+ plt.legend()
+ ax.set(xlabel='Measured', ylabel='Predicted')
+ plt.show()
+ if save_path is not None:
+ plt.savefig(save_path, bbox_inches='tight')
+ return None
+
+[docs] def r_squared(self, compared_value=None):
+ """r-squared value comparing measured and predicted compared value
+
+ Closer to 1, the better the model's predictions.
+
+ :param compared_value: (str) Quantity to compare predicted and
+ experimental data. Can be any column containing "eq" in exp_df i.e.
+ h_eq, z_eq, {RE}_d_eq, etc.
+ """
+ exp_df = self.get_exp_df()
+ predicted_dict = self.get_predicted_dict()
+ re_species_list = self._re_species_list
+ if compared_value is None:
+ compared_value = '{0}_aq_eq'.format(re_species_list[0])
+ pred = predicted_dict[compared_value]
+ predicted_y = np.array(pred)
+ actual_y = exp_df[compared_value].values
+ num = sum((actual_y - predicted_y) ** 2)
+ den = sum((actual_y - np.mean(actual_y)) ** 2)
+ r_2 = (1 - num / den)
+ return r_2
+