diff --git a/.idea/dictionaries/Titus.xml b/.idea/dictionaries/Titus.xml index c7fc181..ff0cf04 100644 --- a/.idea/dictionaries/Titus.xml +++ b/.idea/dictionaries/Titus.xml @@ -12,6 +12,7 @@ extractant ftol kmol + llepe lmse matplotlib maxiter @@ -25,6 +26,7 @@ seaborn slsqp thermo + vals viridis xmls diff --git a/.idea/workspace.xml b/.idea/workspace.xml index e28ee0d..d5984e8 100644 --- a/.idea/workspace.xml +++ b/.idea/workspace.xml @@ -1,11 +1,37 @@ - + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - + - + - + - + - - - - - + + + + + @@ -248,7 +275,14 @@ @@ -269,8 +303,8 @@ - - @@ -290,43 +324,48 @@ - + - - + + - + + + + + + - - - + + + - + - - - + + + - + - - - + + + - + - - - + + + - + - + @@ -335,31 +374,39 @@ - + + + + + - - + + - - - + + + - + + + + + - - - + + + - - + + diff --git a/data/xmls/twophase.xml b/data/xmls/twophase.xml index c10d135..ea2a037 100644 --- a/data/xmls/twophase.xml +++ b/data/xmls/twophase.xml @@ -50,7 +50,7 @@ 298.14999999999998 - -4704703.645715787 + -4704699.156668724 1117.965 0.0 diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle index 021b817..176b86d 100644 Binary files a/docs/build/doctrees/environment.pickle and b/docs/build/doctrees/environment.pickle differ diff --git a/docs/build/doctrees/index.doctree b/docs/build/doctrees/index.doctree index 2378fff..9a0fc6a 100644 Binary files a/docs/build/doctrees/index.doctree and b/docs/build/doctrees/index.doctree differ diff --git a/docs/build/doctrees/modules/reeps.doctree b/docs/build/doctrees/modules/reeps.doctree index fc7c101..55a0dd7 100644 Binary files a/docs/build/doctrees/modules/reeps.doctree and b/docs/build/doctrees/modules/reeps.doctree differ diff --git a/docs/build/html/.buildinfo b/docs/build/html/.buildinfo index 73e9205..d0738d2 100644 --- a/docs/build/html/.buildinfo +++ b/docs/build/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: b8f412ebe5809ccc0d7953f47cdeb07e +config: 09c8b5742a1f96cda53a1848d8ac9ec1 tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/docs/build/html/_modules/index.html b/docs/build/html/_modules/index.html index 142a213..ca853b9 100644 --- a/docs/build/html/_modules/index.html +++ b/docs/build/html/_modules/index.html @@ -8,7 +8,7 @@ - Overview: module code — reeps 1.0.0 documentation + Overview: module code — LLEPE 1.0.0 documentation @@ -48,7 +48,7 @@ - reeps + LLEPE @@ -88,7 +88,7 @@

Searchers

@@ -103,7 +103,7 @@ @@ -150,7 +150,8 @@

All modules for which code is available

-
diff --git a/docs/build/html/_modules/reeps/reeps.html b/docs/build/html/_modules/reeps/reeps.html index 810545e..bb9a06e 100644 --- a/docs/build/html/_modules/reeps/reeps.html +++ b/docs/build/html/_modules/reeps/reeps.html @@ -8,7 +8,7 @@ - reeps.reeps — reeps 1.0.0 documentation + reeps.reeps — LLEPE 1.0.0 documentation @@ -48,7 +48,7 @@ - reeps + LLEPE @@ -88,7 +88,7 @@

Searchers

@@ -103,7 +103,7 @@ @@ -167,6 +167,7 @@ import os import re import pkg_resources +from .utils import set_size sns.set() sns.set(font_scale=1.6) @@ -211,14 +212,15 @@ rare_earth_ion_names = ['Nd+++', 'Pr+++'] - :param exp_csv_filename: (str) csv file name with experimental data + :param exp_data: (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, + [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] @@ -272,7 +274,7 @@ :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 @@ -280,7 +282,7 @@ 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. @@ -300,27 +302,26 @@ :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 + Should have the format as below. Dictionary keys under user defined + parameter name must be named as shown below ('upper_element_name', + 'upper_attrib_name', etc.). 'attrib_name's and 'attrib_value's can + be None. {} denotes areas for user to fill in. .. 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"} - } + opt_dict = {"{user_defined_name_for_parameter_1}": + {'upper_element_name': {param_upper_element}, + 'upper_attrib_name': {param_upper_attrib_name}, + 'upper_attrib_value': {param_upper_attrib_value}, + 'lower_element_name': {param_lower_element}, + 'lower_attrib_name': {param_lower_attrib_name}, + 'lower_attrib_value': {param_lower_attrib_value}, + 'input_format': {str format to input input_value} + 'input_value': {guess_value}}, + "{user_defined_name_for_parameter_2}": + ... + ... + } :param objective_function: (function or str) function to compute objective By default, the objective function is log mean squared error @@ -337,7 +338,7 @@ 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 @@ -381,14 +382,8 @@ "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', @@ -403,7 +398,8 @@ ``kwargs`` is optional - Function needs to return: (np.ndarray) Optimized parameters + Function needs to return: ((np.ndarray, float)) Optimized parameters, + objective_function value :param temp_xml_file_path: (str) path to temporary xml file. @@ -412,10 +408,13 @@ xml file default is local temp folder + + :param dependant_params_dict: (dict) dictionary containing information + about parameters dependant on opt_dict """ def __init__(self, - exp_csv_filename, + exp_data, phases_xml_filename, phase_names, aq_solvent_name, @@ -430,11 +429,12 @@ opt_dict=None, objective_function='Log-MSE', optimizer='SLSQP', - temp_xml_file_path=None + temp_xml_file_path=None, + dependant_params_dict=None, ): self._built_in_obj_list = ['Log-MSE'] self._built_in_opt_list = ['SLSQP'] - self._exp_csv_filename = exp_csv_filename + self._exp_data = exp_data self._phases_xml_filename = phases_xml_filename self._opt_dict = opt_dict self._phase_names = phase_names @@ -454,6 +454,7 @@ 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 + self._dependant_params_dict = dependant_params_dict # 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 @@ -471,12 +472,15 @@ 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) + if isinstance(self._exp_data, str): + try: + self._exp_df = pd.read_csv(self._exp_data) + except FileNotFoundError: + self._exp_data = pkg_resources.resource_filename( + 'reeps', r'..\data\csvs\{0}'.format(self._exp_data)) + self._exp_df = pd.read_csv(self._exp_data) + else: + self._exp_df = self._exp_data.copy() self._exp_df_columns = ['h_i', 'h_eq', 'z_i', 'z_eq'] if self._re_species_list is None: @@ -488,7 +492,12 @@ 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 + for species in self._re_species_list: + self._exp_df['{0}_org_eq'.format(species)] = \ + self._exp_df['{0}_aq_eq'.format(species)] \ + * self._exp_df['{0}_d_eq'.format(species)] self._in_moles = None @@ -500,16 +509,18 @@ self._predicted_dict = None self.update_predicted_dict() -
[docs] @staticmethod - def slsqp_optimizer(objective, x_guess): +
[docs] @staticmethod + def scipy_minimize(objective, x_guess, optimizer_kwargs=None): """ The default optimizer for REEPS - Uses scipy.minimize with options + Uses scipy.minimize + + By default, options are .. code-block:: python default_kwargs= {"method": 'SLSQP', - "bounds": [(1e-1, 1e1)*len(x_guess)], + "bounds": [(1e-1, 1e1)]*len(x_guess), "constraints": (), "options": {'disp': True, 'maxiter': 1000, @@ -517,17 +528,20 @@ :param objective: (func) the objective function :param x_guess: (np.ndarray) the initial guess (always 1) - :returns: (np.ndarray) Optimized parameters + :param optimizer_kwargs: (dict) dictionary of options for minimize + :returns: ((np.ndarray, float)) Optimized parameters, + objective_function value """ - optimizer_kwargs = {"method": 'SLSQP', - "bounds": [(1e-1, 1e1)] * len(x_guess), - "constraints": (), - "options": {'disp': True, - 'maxiter': 1000, - 'ftol': 1e-6}} + if optimizer_kwargs is None: + 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
+ return est_parameters, res.fun
[docs] def log_mean_squared_error(self, predicted_dict, meas_df): """Default objective function for REEPS @@ -549,7 +563,8 @@ 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) + log_diff = (log_pred - log_meas) ** 2 + obj = np.sum(log_diff) return obj
[docs] def get_exp_df(self) -> pd.DataFrame: @@ -559,7 +574,7 @@ """ return self._exp_df
-
[docs] def set_exp_df(self, exp_csv_filename): +
[docs] def set_exp_df(self, exp_data): """Changes the experimental DataFrame to input exp_csv_filename data and renames columns to internal REEPS names @@ -568,15 +583,19 @@ See class docstring on "exp_csv_filename" for further explanations. - :param exp_csv_filename: (str) file name/path for experimental data csv + :param exp_data: (str or pd.DataFrame) + file name/path or DataFrame 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_data = exp_data + if isinstance(self._exp_data, str): + try: + self._exp_df = pd.read_csv(self._exp_data) + except FileNotFoundError: + self._exp_data = pkg_resources.resource_filename( + 'reeps', r'..\data\csvs\{0}'.format(self._exp_data)) + self._exp_df = pd.read_csv(self._exp_data) + else: + self._exp_df = exp_data.copy() self._exp_df_columns = ['h_i', 'h_eq', 'z_i', 'z_eq'] if self._re_species_list is None: self._re_species_list = [] @@ -588,6 +607,10 @@ 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 + for species in self._re_species_list: + self._exp_df['{0}_org_eq'.format(species)] = \ + self._exp_df['{0}_aq_eq'.format(species)] \ + * self._exp_df['{0}_d_eq'.format(species)] self.set_in_moles(feed_vol=1) self.update_predicted_dict() return None
@@ -604,7 +627,7 @@ """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 + Also runs set_in_moles to set feed volume to 1 L :param phases_xml_filename: (str) xml file with parameters for equilibrium calc @@ -720,8 +743,8 @@
[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 + + :return: rare_earth_ion_names: (list) names of rare earth ions in xml file """ return self._rare_earth_ion_names
@@ -730,7 +753,7 @@ """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 + :param rare_earth_ion_names: (list) names of rare earth ions in xml file """ self._rare_earth_ion_names = rare_earth_ion_names @@ -756,15 +779,15 @@
[docs] def get_aq_solvent_rho(self) -> str: """Returns aqueous solvent density (g/L) - - :return: aq_solvent_rho: (float) density of aqueous solvent + + :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 + :param aq_solvent_rho: (float) density of aqueous solvent """ self._aq_solvent_rho = aq_solvent_rho return None
@@ -772,14 +795,14 @@
[docs] def get_extractant_rho(self) -> str: """Returns extractant density (g/L) - :return: extractant_rho: (float) density of extractant + :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 + :param extractant_rho: (float) density of extractant """ self._extractant_rho = extractant_rho return None
@@ -787,14 +810,14 @@
[docs] def get_diluant_rho(self) -> str: """Returns diluant density (g/L) - :return: diluant_rho: (float) density of diluant + :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 + :param diluant_rho: (float) density of diluant """ self._diluant_rho = diluant_rho return None
@@ -810,7 +833,7 @@ This function also calls update_predicted_dict - :param feed_vol: (float) feed volume of mixture (g/L) + :param feed_vol: (float) feed volume of mixture (L) """ phases_copy = self._phases.copy() exp_df = self._exp_df.copy() @@ -963,7 +986,7 @@ "with at least 2 arguments: " "f(objective_func,x_guess, kwargs)") if optimizer == 'SLSQP': - optimizer = self.slsqp_optimizer + optimizer = self.scipy_minimize self._optimizer = optimizer return None
@@ -998,6 +1021,23 @@ self._temp_xml_file_path = temp_xml_file_path return None +
[docs] def get_dependant_params_dict(self): + """ + Returns the dependant_params_dict + :return: dependant_params_dict: (dict) dictionary containing + information about parameters dependant on opt_dict + """ + return self._dependant_params_dict
+ +
[docs] def set_dependant_params_dict(self, dependant_params_dict): + """ + Sets the dependant_params_dict + :param dependant_params_dict: (dict) dictionary containing information + about parameters dependant on opt_dict + """ + self._dependant_params_dict = dependant_params_dict + return None
+
[docs] def update_predicted_dict(self, phases_xml_filename=None, phase_names=None): @@ -1080,6 +1120,7 @@ exp_df = self._exp_df objective_function = self._objective_function opt_dict = copy.deepcopy(self._opt_dict) + dep_dict = copy.deepcopy(self._dependant_params_dict) x = np.array(x) if len(x.shape) == 1: @@ -1093,9 +1134,14 @@ 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] + if not np.isnan( + x[i]): # if nan, do not update xml with nan + opt_dict[species_name][thermo_prop] *= x[i] i += 1 - self.update_xml(opt_dict, temp_xml_file_path) + + self.update_xml(opt_dict, + temp_xml_file_path, + dependant_params_dict=dep_dict) self.update_predicted_dict(temp_xml_file_path) predicted_dict = self.get_predicted_dict() @@ -1118,16 +1164,19 @@ objective_function=None, optimizer=None, objective_kwargs=None, - optimizer_kwargs=None) -> dict: + optimizer_kwargs=None) -> tuple: """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 + If 'None', last set objective or default function is used :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 'None', last set optimizer or default is used + :param optimizer_kwargs: (dict) optional arguments for optimizer + :param objective_kwargs: (dict) optional arguments + for objective function + :returns tuple: (opt_dict (dict), opt_value (float)) + optimized opt_dict: Has identical structure as opt_dict """ if objective_function is not None: self.set_objective_function(objective_function) @@ -1147,21 +1196,24 @@ if optimizer_kwargs is None: # noinspection PyCallingNonCallable - est_parameters = optimizer(objective, x_guess) + est_parameters, obj_value = optimizer(objective, x_guess) else: # noinspection PyCallingNonCallable - est_parameters = optimizer(objective, x_guess, **optimizer_kwargs) + est_parameters, obj_value = 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
+ return opt_dict, obj_value
[docs] def update_xml(self, info_dict, - phases_xml_filename=None): + phases_xml_filename=None, + dependant_params_dict=None): """updates xml file with info_dict :param info_dict: (dict) info in {species_names:{thermo_prop:val}} @@ -1169,20 +1221,42 @@ :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. + :param dependant_params_dict: (dict) dictionary containing information + about parameters dependant on info_dict """ if phases_xml_filename is None: phases_xml_filename = self._phases_xml_filename + new_dict = copy.deepcopy(info_dict) + dep_dict = dependant_params_dict + if dep_dict is not None: + for species_name in dep_dict.keys(): + for thermo_prop in dep_dict[species_name]: + mod_func = \ + dep_dict[species_name][thermo_prop]['function'] + mod_kwargs = \ + dep_dict[species_name][thermo_prop]['kwargs'] + ind_vars = \ + dep_dict[species_name][thermo_prop]['ind_vars'] + ind_vals = [new_dict[ind_var[0]][ind_var[1]] + for ind_var in ind_vars] + + new_dict[species_name] = {} + new_dict[species_name][thermo_prop] = {} + new_dict[species_name][thermo_prop] = \ + mod_func(ind_vals, **mod_kwargs) + # print(mod_func(ind_vals, **mod_kwargs)) + # print(new_dict) 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_name in new_dict.keys(): + for thermo_prop in new_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]) + new_dict[species_name][thermo_prop]) now = datetime.now() changed_prop.set('updated', 'Updated at {0}:{1} {2}-{3}-{4}' @@ -1195,10 +1269,213 @@ self.set_phases(self._phases_xml_filename, self._phase_names) return None
+ def _internal_objective_ver2(self, x, kwargs=None): + """ + ver2 generalizes to handle accessing parameters. ver1 assumes species + parameter is modified. ver2 assumes parameter is accessed by going + through two levels: upper and lower + 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) + dep_dict = copy.deepcopy(self._dependant_params_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: + for ind, param_name in enumerate(opt_dict.keys()): + if not np.isnan( + x[ind]): # if nan, do not update xml with nan + opt_dict[param_name]['input_value'] *= x[ind] + + self.update_xml_ver2(opt_dict, + temp_xml_file_path, + dependant_params_dict=dep_dict) + + 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_ver2(self, + objective_function=None, + optimizer=None, + objective_kwargs=None, + optimizer_kwargs=None) -> tuple: + """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 + If 'None', last set objective or default function is used + :param optimizer: (function) function to perform optimization + If 'None', last set optimizer or default is used + :param optimizer_kwargs: (dict) optional arguments for optimizer + :param objective_kwargs: (dict) optional arguments + for objective function + :returns tuple: (opt_dict (dict), opt_value (float)) + 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_ver2(x, objective_kwargs) + + optimizer = self._optimizer + opt_dict = copy.deepcopy(self._opt_dict) + x_guess = np.ones(len(list(opt_dict.keys()))) + + if optimizer_kwargs is None: + # noinspection PyCallingNonCallable + est_parameters, obj_value = optimizer(objective, x_guess) + else: + # noinspection PyCallingNonCallable + est_parameters, obj_value = optimizer(objective, + x_guess, + optimizer_kwargs) + for ind, param_name in enumerate(opt_dict.keys()): + opt_dict[param_name]['input_value'] *= est_parameters[ind] + + return opt_dict, obj_value
+ +
[docs] def update_xml_ver2(self, + info_dict, + phases_xml_filename=None, + dependant_params_dict=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. + :param dependant_params_dict: (dict) dictionary containing information + about parameters dependant on info_dict + """ + if phases_xml_filename is None: + phases_xml_filename = self._phases_xml_filename + new_dict = copy.deepcopy(info_dict) + dep_dict = dependant_params_dict + + if dep_dict is not None: + new_dict.update(dep_dict) + for param_name in dep_dict.keys(): + mod_func = \ + dep_dict[param_name]['function'] + mod_kwargs = \ + dep_dict[param_name]['kwargs'] + if isinstance(dep_dict[param_name]['independent_params'], str): + ind_param_names = [dep_dict[ + param_name]['independent_params']] + else: + ind_param_names = \ + dep_dict[param_name]['independent_params'] + ind_vals = [new_dict[ind_param_name]['input_value'] + for ind_param_name in ind_param_names] + if mod_kwargs is None: + new_dict[param_name]['input_value'] = mod_func(ind_vals) + else: + new_dict[param_name]['input_value'] = \ + mod_func(ind_vals, + **mod_kwargs) + tree = ET.parse(phases_xml_filename) + root = tree.getroot() + # Update xml file + for key in list(new_dict.keys()): + d = new_dict[key] + now = datetime.now() + if (d['upper_attrib_name'] is not None + and d['lower_attrib_name'] is not None): + for child1 in root.iter(d['upper_element_name']): + if (child1.attrib[d['upper_attrib_name']] + == d['upper_attrib_value']): + for child2 in child1.iter(d['lower_element_name']): + if (child1.attrib[d['lower_attrib_name']] + == d['lower_attrib_value']): + child2.text = d['input_format'].format( + d['input_value']) + child2.set('updated', + 'Updated at {0}:{1} {2}-{3}-{4}' + .format(now.hour, now.minute, + now.month, now.day, + now.year)) + elif (d['upper_attrib_name'] is None + and d['lower_attrib_name'] is not None): + for child1 in root.iter(d['upper_element_name']): + for child2 in child1.iter(d['lower_element_name']): + if (child1.attrib[d['lower_attrib_name']] + == d['lower_attrib_value']): + child2.text = d['input_format'].format( + d['input_value']) + child2.set('updated', + 'Updated at {0}:{1} {2}-{3}-{4}' + .format(now.hour, now.minute, + now.month, now.day, + now.year)) + elif (d['upper_attrib_name'] is not None + and d['lower_attrib_name'] is None): + for child1 in root.iter(d['upper_element_name']): + if (child1.attrib[d['upper_attrib_name']] + == d['upper_attrib_value']): + for child2 in child1.iter(d['lower_element_name']): + child2.text = d['input_format'].format( + d['input_value']) + child2.set('updated', + 'Updated at {0}:{1} {2}-{3}-{4}' + .format(now.hour, now.minute, + now.month, now.day, + now.year)) + else: + for child1 in root.iter(d['upper_element_name']): + for child2 in child1.iter(d['lower_element_name']): + child2.text = d['input_format'].format( + d['input_value']) + child2.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, + c_data=None, + c_label=None, + plot_title=None, save_path=None, - print_r_squared=False): + print_r_squared=False, + data_labels=None, + legend=True): """ Parity plot between measured and predicted compared_value. Default compared value is {RE_1}_aq_eq @@ -1206,9 +1483,25 @@ :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 plot_title: (str or boolean) + + If None (default): Plot title will be generated from compared_value + Recommend to just explore. If h_eq, plot_title is + "H^+ eq conc". + + If str: Plot title will be plot_title string + + If "False": No plot title + :param c_data: (list or np.ndarray) data for color axis + :param c_label: (str) label for color axis :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 + :param data_labels: labels for the data such as paper's name where + experiment is pulled from. + :param legend: whether to display legend for data_labels. Has no + effect if data_labels is None + :return fig, ax: returns the figure and axes objects """ exp_df = self.get_exp_df() predicted_dict = self.get_predicted_dict() @@ -1217,45 +1510,98 @@ 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 + pred = pd.DataFrame(predicted_dict)[compared_value].fillna(0).values + meas = exp_df[compared_value].fillna(0).values + name_breakdown = re.findall('[^_\W]+', compared_value) + compared_species = name_breakdown[0] + if compared_species == 'h': + feed_molarity = exp_df['h_i'].fillna(0).values + elif compared_species == 'z': + feed_molarity = exp_df['z_i'].fillna(0).values + else: + feed_molarity = exp_df[ + '{0}_aq_i'.format(compared_species)].fillna(0).values + if isinstance(data_labels, list): + combined_df = pd.DataFrame({'pred': pred, + 'meas': meas, + 'label': data_labels, + 'feed_molarity': feed_molarity}) + elif isinstance(c_data, str): + combined_df = pd.DataFrame({'pred': pred, + 'meas': meas, + c_data: exp_df[c_data].values, + 'feed_molarity': feed_molarity}) + else: + combined_df = pd.DataFrame({'pred': pred, + 'meas': meas, + 'feed_molarity': feed_molarity}) + + combined_df = combined_df[(combined_df['feed_molarity'] != 0)] + meas = combined_df['meas'].values + pred = combined_df['pred'].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^+$' + default_title = '$H^+$ eq. conc. (mol/L)' elif compared_species == 'z': - species_name = extractant_name + default_title = '{0} eq. conc. (mol/L)'.format(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) + default_title = '$%s^{%d+}$ eq. conc. (mol/L)' \ + % (compared_species, re_charge) elif phase == 'd': - species_name = '{0} distribution ratio'.format( + default_title = '{0} distribution ratio'.format( compared_species) else: - species_name = '{0} complex'.format(compared_species) + default_title = '{0} complex eq. conc. (mol/L)'.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 isinstance(data_labels, list): + unique_labels = list(set(data_labels)) + for label in unique_labels: + filtered_data = combined_df[combined_df['label'] == label] + filtered_meas = filtered_data['meas'] + filtered_pred = filtered_data['pred'] + ax.scatter(filtered_meas, filtered_pred, label=label) + if legend: + ax.legend(loc='best') + + elif c_data is not None: + if isinstance(c_data, str): + c_data = combined_df[c_data].values + p1 = ax.scatter(meas, pred, c=c_data, alpha=1, cmap='viridis') + c_bar = fig.colorbar(p1, format='%.2f') + if c_label is not None: + c_bar.set_label(c_label, rotation=270, labelpad=20) + else: + sns.scatterplot(meas, pred, color="r", + legend=False) + ax.plot(min_max_data, min_max_data, color="b", label="") + if print_r_squared: - p1.text(min_max_data[0], + ax.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() + # plt.legend(loc='lower right') + # else: + # plt.legend() + ax.set(xlabel='Measured', ylabel='Predicted') + if plot_title is None: + ax.set_title(default_title) + elif isinstance(plot_title, str): + ax.set_title(plot_title) + set_size(8, 6) + plt.tight_layout() plt.show() if save_path is not None: plt.savefig(save_path, bbox_inches='tight') - return None
+ return fig, ax
[docs] def r_squared(self, compared_value=None): """r-squared value comparing measured and predicted compared value @@ -1264,20 +1610,87 @@ :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. + h_eq, z_eq, {RE}_d_eq, etc. default is {RE}_aq_eq """ 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] + pred = pd.DataFrame(predicted_dict)[compared_value].fillna(0).values predicted_y = np.array(pred) - actual_y = exp_df[compared_value].values + actual_y = exp_df[compared_value].fillna(0).values + name_breakdown = re.findall('[^_\W]+', compared_value) + compared_species = name_breakdown[0] + if compared_species == 'h': + feed_molarity = exp_df['h_i'].fillna(0).values + elif compared_species == 'z': + feed_molarity = exp_df['z_i'].fillna(0).values + else: + feed_molarity = exp_df[ + '{0}_aq_i'.format(compared_species)].fillna(0).values + combined_df = pd.DataFrame({'pred': predicted_y, + 'meas': actual_y, + 'in_moles': feed_molarity}) + combined_df = combined_df[(combined_df['in_moles'] != 0)] + actual_y = combined_df['meas'].values + predicted_y = combined_df['pred'].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
+ if den == 0: + r_2 = 0 + else: + r_2 = (1 - num / den) + return r_2 + +
[docs] @staticmethod + def plot_3d_data(x_data, + y_data, + z_data, + c_data=None, + x_label=None, + y_label=None, + z_label=None, + c_label=None): + """ + + :param x_data: (list) list of data for x axis + :param y_data: (list) list of data for y axis + :param z_data: (list) list of data for z axis + :param c_data: (list) list of data for color axis + :param x_label: (str) label for x axis + :param y_label: (str) label for y axis + :param z_label: (str) label for z axis + :param c_label: (str) label for color axis + :return: + """ + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + if c_data is None: + ax.plot(x_data, y_data, z_data, 'o') + else: + p1 = ax.scatter(x_data, + y_data, + z_data, 'o', c=c_data, + cmap='viridis', alpha=1) + c_bar = fig.colorbar(p1) + if c_label is not None: + c_bar.set_label(c_label, rotation=270, labelpad=20) + if x_label is None: + ax.set_xlabel('x', labelpad=15) + else: + ax.set_xlabel(x_label, labelpad=15) + if y_label is None: + ax.set_ylabel('y', labelpad=15) + else: + ax.set_ylabel(y_label, labelpad=15) + if z_label is None: + ax.set_zlabel('z', labelpad=15) + else: + ax.set_zlabel(z_label, labelpad=15) + plt.show() + return fig, ax
diff --git a/docs/build/html/_sources/index.rst.txt b/docs/build/html/_sources/index.rst.txt index 0cef716..9a973a9 100644 --- a/docs/build/html/_sources/index.rst.txt +++ b/docs/build/html/_sources/index.rst.txt @@ -3,12 +3,11 @@ You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. -Welcome to reeps's docs! - the Rare earth element parameter searcher -==================================================================== -REEPS is a package for thermodynamic parameter estimation specifically -for rare earth element extraction modeling. +Welcome to llepe's docs! - the Liquid-Liquid Extraction Parameter Estimator +=========================================================================== +LLEPE is a package for thermodynamic parameter estimation for liquid-liquid extraction modeling -REEPS takes experimental data in a csv and system data in a xml. +LLEPE takes experimental data in a csv and system data in a xml. The package then uses Cantera, another python package, to simulate equilibrium. diff --git a/docs/build/html/_sources/modules/reeps.rst.txt b/docs/build/html/_sources/modules/reeps.rst.txt index 9be4450..bba5c13 100644 --- a/docs/build/html/_sources/modules/reeps.rst.txt +++ b/docs/build/html/_sources/modules/reeps.rst.txt @@ -1,8 +1,8 @@ .. _reeps: -.. automodule:: reeps +.. automodule:: llepe -REEPS +LLEPE ===== @@ -10,6 +10,6 @@ REEPS Parameters ---------- -.. autoclass:: REEPS +.. autoclass:: LLEPE :members: :inherited-members: \ No newline at end of file diff --git a/docs/build/html/genindex.html b/docs/build/html/genindex.html index 195d0e4..d1e1f7b 100644 --- a/docs/build/html/genindex.html +++ b/docs/build/html/genindex.html @@ -9,7 +9,7 @@ - Index — reeps 1.0.0 documentation + Index — LLEPE 1.0.0 documentation @@ -49,7 +49,7 @@ - reeps + LLEPE @@ -89,7 +89,7 @@

Searchers

@@ -104,7 +104,7 @@ @@ -169,7 +169,7 @@

F

@@ -177,41 +177,43 @@

G

@@ -219,7 +221,18 @@

L

+
@@ -231,7 +244,7 @@ module @@ -240,7 +253,11 @@

P

+
@@ -248,18 +265,7 @@

R

-
@@ -267,41 +273,43 @@

S

@@ -309,11 +317,11 @@

U

diff --git a/docs/build/html/guide/install.html b/docs/build/html/guide/install.html index 36945cc..0eed037 100644 --- a/docs/build/html/guide/install.html +++ b/docs/build/html/guide/install.html @@ -8,7 +8,7 @@ - Installation — reeps 1.0.0 documentation + Installation — LLEPE 1.0.0 documentation @@ -36,7 +36,7 @@ - + @@ -50,7 +50,7 @@ - reeps + LLEPE @@ -97,7 +97,7 @@

Searchers

@@ -112,7 +112,7 @@ @@ -204,7 +204,7 @@ pip install -e .[docs,tests] - + diff --git a/docs/build/html/guide/quickstart.html b/docs/build/html/guide/quickstart.html index 69b4264..ed14293 100644 --- a/docs/build/html/guide/quickstart.html +++ b/docs/build/html/guide/quickstart.html @@ -8,7 +8,7 @@ - Getting Started — reeps 1.0.0 documentation + Getting Started — LLEPE 1.0.0 documentation @@ -35,7 +35,7 @@ - + @@ -50,7 +50,7 @@ - reeps + LLEPE @@ -90,7 +90,7 @@

Searchers

@@ -105,7 +105,7 @@ @@ -193,7 +193,7 @@ aqueous phase.