Page MenuHomec4science

solver.pyx
No OneTemporary

File Metadata

Created
Sat, May 4, 09:07

solver.pyx

#Copyright (c) 2014,2015 Fabien Georget <fabieng@princeton.edu>, Princeton
#University #All rights reserved.
#
#Redistribution and use in source and binary forms, with or without
#modification, are permitted provided that the following conditions are met:
#
#1. Redistributions of source code must retain the above copyright notice, this
#list of conditions and the following disclaimer.
#
#2. Redistributions in binary form must reproduce the above copyright notice,
#this list of conditions and the following disclaimer in the documentation
#and/or other materials provided with the distribution.
#
#3. Neither the name of the copyright holder nor the names of its contributors
#may be used to endorse or promote products derived from this software without
#specific prior written permission.
#
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
#WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
#DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
#CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
#OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
from libcpp.map cimport map as cmap
from libcpp.string cimport string
from libcpp.utility cimport pair
from solver cimport SpecMiCPSolver, AdimensionalSystemSolver
from database cimport DatabaseManager
from constraints cimport SpecMiCPConstraints, AdimensionalSystemConstraints
from solver_options cimport MiCPSolverReturnCode, MiCPPerformance, is_solver_successful
from solution cimport SpecMiCPSolution
cimport units
from cython.operator cimport dereference
cdef class SpecMiCPSolver:
def __cinit__(
self,
DatabaseManager database not None,
SpecMiCPConstraints constraints not None,
SpecMiCPSolution solution
):
"""Initialisation"""
self.database = database
if solution is None:
self.solver = new AdimensionalSystemSolver(
database.get_raw_db(),
dereference(constraints._get())
)
else:
self.solver = new AdimensionalSystemSolver(
database.get_raw_db(),
dereference(constraints._get()),
dereference(solution._get())
)
self.variables = new VectorXd(solution.get_main_variables())
def __dealloc__(self):
if self.solver:
del self.solver
if self.variables:
del self.variables
cdef AdimensionalSystemSolverOptions* get_options(self) nogil:
"""Return the options"""
return &(self.solver.get_options())
cdef AdimensionalSystemOptions* get_system_options(self) nogil:
"""Return the options of the SpecMiCP program"""
return &(self.solver.get_options().system_options)
cdef MiCPSolverOptions* get_micpsolver_options(self) nogil:
"""Return the MiCPSolver options"""
return &(self.solver.get_options().solver_options)
def get_solution(self):
"""Return the solution"""
if not self.solver:
raise RuntimeError("No solver")
solution = SpecMiCPSolution()
solution.set_solution(
self.solver.get_raw_solution(dereference(self.variables)),
self.database.get_raw_db(),
self.solver.get_options().units_set
)
return solution
def initialize_variables(self, water_saturation, log_molalities_component):
"""Initialize the variables
The log_molalities_component can be
- a double : the initial value gor the log_10 of the molalities
- a dictionnary {"label component": log_molalities}
"""
cdef cmap[string, double] map_comp
if not self.variables:
self.variables = new VectorXd()
if hasattr(log_molalities_component, "items"):
for (key, value) in log_molalities_component.items():
map_comp.insert(pair[string, double](key, value))
self.solver.initialise_variables(dereference(self.variables),
float(water_saturation),
map_comp)
else:
self.solver.initialise_variables(dereference(self.variables),
float(water_saturation),
float(log_molalities_component)
)
def solve(self, initialize=True):
"""Solve the problem"""
cdef MiCPPerformance perf
cdef bool init = initialize
if not self.variables:
self.variables = new VectorXd()
with nogil:
perf = self.solver.solve(dereference(self.variables), init)
if not is_solver_successful(perf.return_code):
raise RuntimeError("Error - the system cannot be solved.")
# Options
# =======
# adimensionalsolver options
# --------------------------
def set_length_unit(self, length_unit):
"""Set the length unit"""
if length_unit in ("m", "meter"):
units.set_meter(self.get_options().units_set.length)
elif length_unit in ("dm", "decimeter"):
units.set_decimeter(self.get_options().units_set.length)
elif length_unit in ("cm", "centimeter"):
units.set_centimeter(self.get_options().units_set.length)
else:
raise ValueError("Unit : " + str(length_unit) + "is not recognized")
def enable_restart(self):
"""Allow the solve to be restarted in case of failure"""
self.get_options().allow_restart = True
def disable_restart(self):
"""Failure to solve the problem is critical"""
self.get_options().allow_restart = False
def enable_pcfm(self):
"""Enable the positive continuous fraction method"""
self.get_options().use_pcfm = True
def disable_pcfm(self):
"""Disable the positive continuous fraction method"""
self.get_options().use_pcfm = False
# System options
# ---------------
def enable_non_ideality(self):
"""Enable the nonideality model of the aqueous solution"""
self.get_system_options().non_ideality = True
def disable_non_ideality(self):
"""Disable the nonideality model of the aqueous solution"""
self.get_system_options().non_ideality = False
def set_non_ideality_max_iterations(self, max_iterations):
"""Set the maximum iterations allowed for the non ideality model"""
self.get_system_options().non_ideality_max_iter = int(max_iterations)
def set_non_ideality_tolerance(self, tolerance):
"""Set the tolerance for the non ideality model"""
self.get_system_options().non_ideality_tolerance = float(tolerance)
def set_non_ideality_start_factor(self, start_factor):
"""A factor that governs the non-ideality computation"""
self.get_system_options().start_non_ideality_computation = float(start_factor)
def set_under_relaxation_factor(self, factor):
"""Set the under relaxtion factor """
self.get_system_options().under_relaxation_factor = float(factor)
def set_restart_concentration(self, log_molality):
"""Set the molality used to re-initialize the problem in case of failure"""
self.get_system_options().restart_concentration = float(log_molality)
def set_new_component_concentration(self, log_molality):
"""Set the concentration for a new component"""
self.get_system_options().new_component_concentration = float(log_molality)
def set_cutoff_total_concentrations(self, cutoff):
"""Set the cutoff for the total concentrations"""
self.get_system_options().cutoff_total_concentration = float(cutoff)
# MiCPSolver options
# ------------------
def set_tolerance(self, residuals_tol):
"""Set the residuals tolerance"""
self.get_micpsolver_options().set_tolerance(residuals_tol)
def set_tolerances(self, residuals_tol, step_tol):
""""Set the residuals and the step tolerances"""
self.get_micpsolver_options().set_tolerance(residuals_tol, step_tol)
def set_maximum_iterations(self, max_iterations):
"""Set the maximum number of operations"""
self.get_micpsolver_options().set_maximum_iterations(max_iterations)
def set_maximum_step_length(self, max_step_length, max_iter_at_max_length=None):
"""Set the maximum step length"""
if max_iter_at_max_length is not None:
self.get_micpsolver_options().set_maximum_step_length(
max_step_length, max_iter_at_max_length)
else:
self.get_micpsolver_options().set_maximum_step_length(
max_step_length)
def disable_descent_direction(self):
"""Disable the descent direction condition"""
self.get_micpsolver_options().disable_descent_direction()
def enable_descent_direction(self, factor, power):
"""Enable the descent direction condition"""
self.get_micpsolver_options().enable_descent_direction(factor, power)
def disable_condition_check(self):
"""Disable the descent direction"""
self.get_micpsolver_options().disable_condition_check()
def enable_condition_check(self, threshold):
"""Enable the condition check"""
self.get_micpsolver_options().enable_condition_check(threshold)
def disable_non_monotone_linesearch(self):
"""Disable the non monotone linesearch"""
self.get_micpsolver_options().disable_non_monotone_linesearch()
def enable_non_monotone_linesearch(self):
"""Enable the non monotone linesearch"""
self.get_micpsolver_options().enable_non_monotone_linesearch()
def disable_jacobian_scaling(self):
"""Disable the Jacobian scaling"""
self.get_micpsolver_options().disable_scaling()
def enable_jacobian_scaling(self):
"""Enable the Jacobian scaling"""
self.get_micpsolver_options().enable_scaling()
def disable_crashing(self):
"""Disable the crashing"""
self.get_micpsolver_options().disable_crashing()
def enable_crashing(self):
"""Enable the crashing"""
self.get_micpsolver_options().enable_crashing()

Event Timeline