From ff9d521789da87885edcb2bc18ee802c8b2b5b46 Mon Sep 17 00:00:00 2001 From: Eslick Date: Thu, 28 Mar 2024 18:43:37 -0400 Subject: [PATCH 01/13] Add a function for cubic spline interpolation --- pyomo/contrib/cspline_1d_external/__init__.py | 0 pyomo/contrib/cspline_1d_external/build.py | 32 ++ .../cubic_spline_parameters.py | 366 ++++++++++++++++++ pyomo/contrib/cspline_1d_external/plugins.py | 19 + .../cspline_1d_external/src/CMakeLists.txt | 43 ++ .../cspline_1d_external/src/FindASL.cmake | 106 +++++ .../cspline_1d_external/src/functions.cpp | 144 +++++++ .../cspline_1d_external/test/__init__.py | 0 .../cspline_1d_external/test/function.py | 38 ++ .../cspline_1d_external/test/t1_params.txt | 22 ++ 10 files changed, 770 insertions(+) create mode 100644 pyomo/contrib/cspline_1d_external/__init__.py create mode 100644 pyomo/contrib/cspline_1d_external/build.py create mode 100644 pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py create mode 100644 pyomo/contrib/cspline_1d_external/plugins.py create mode 100644 pyomo/contrib/cspline_1d_external/src/CMakeLists.txt create mode 100644 pyomo/contrib/cspline_1d_external/src/FindASL.cmake create mode 100644 pyomo/contrib/cspline_1d_external/src/functions.cpp create mode 100644 pyomo/contrib/cspline_1d_external/test/__init__.py create mode 100644 pyomo/contrib/cspline_1d_external/test/function.py create mode 100644 pyomo/contrib/cspline_1d_external/test/t1_params.txt diff --git a/pyomo/contrib/cspline_1d_external/__init__.py b/pyomo/contrib/cspline_1d_external/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pyomo/contrib/cspline_1d_external/build.py b/pyomo/contrib/cspline_1d_external/build.py new file mode 100644 index 00000000000..8ca4119dbe1 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/build.py @@ -0,0 +1,32 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import sys +from pyomo.common.cmake_builder import build_cmake_project + + +def build_cspline_1d_external(user_args=[], parallel=None): + return build_cmake_project( + targets=["src"], + package_name="cspline_1d_external", + description="AMPL external interpolation function library", + user_args=["-DBUILD_AMPLASL_IF_NEEDED=ON"] + user_args, + parallel=parallel, + ) + + +class AMPLCspline1DExternalBuilder(object): + def __call__(self, parallel): + return build_cspline_1d_external(parallel=parallel) + + +if __name__ == "__main__": + build_cspline_1d_external(sys.argv[1:]) diff --git a/pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py b/pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py new file mode 100644 index 00000000000..472858b8061 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py @@ -0,0 +1,366 @@ +from pyomo.common.dependencies import attempt_import +import pyomo.environ as pyo + +np, numpy_available = attempt_import("numpy") +plt, plt_available = attempt_import("matplotlib.pyplot") + +import idaes + + +def f(x, alpha, s=None): + """ + Cubic function: + y = a1 + a2 * x + a3 * x^2 + a4 * x^3 + + Optionally if s is provided it is a segment index. + y = a1[s] + a2[s] * x + a3[s] * x^2 + a4[s] * x^3 + + This is used to write constraints more compactly. + + Args: + x: x variable, numeric, numpy array, or Pyomo component + alpha: cubic parameters, numeric or Pyomo component + s: optional segment index + + Returns: + Pyomo expression, numpy array, or float + """ + if s is None: + return alpha[1] + alpha[2] * x + alpha[3] * x**2 + alpha[4] * x**3 + return alpha[s, 1] + alpha[s, 2] * x + alpha[s, 3] * x**2 + alpha[s, 4] * x**3 + + +def fx(x, alpha, s=None): + """ + Cubic function first derivative: + dy/dx = a2 + 2 * a3 * x + 3 * a4 * x^2 + + Optionally if s is provided it is a segment index. + dy/dx = a2[s] + 2 * a3[s] * x + 3 * a4[s] * x^2 + + This is used to write constraints more compactly. + + Args: + x: x variable, numeric, numpy array, or Pyomo component + alpha: cubic parameters, numeric or Pyomo component + s: optional segment index + + Returns: + Pyomo expression, numpy array, or float + """ + if s is None: + return alpha[2] + 2 * alpha[3] * x + 3 * alpha[4] * x**2 + return alpha[s, 2] + 2 * alpha[s, 3] * x + 3 * alpha[s, 4] * x**2 + + +def fxx(x, alpha, s=None): + """ + Cubic function second derivative: + d2y/dx2 = 2 * a3 + 6 * a4 * x + + Optionally if s is provided it is a segment index. + d2y/dx2 = 2 * a3[s] + 6 * a4[s] * x + + This is used to write constraints more compactly. + + Args: + x: x variable, numeric, numpy array, or Pyomo component + alpha: cubic parameters, numeric or Pyomo component + s: optional segment index + + Returns: + Pyomo expression, numpy array, or float + """ + if s is None: + return 2 * alpha[3] + 6 * alpha[4] * x + return 2 * alpha[s, 3] + 6 * alpha[s, 4] * x + + +def cubic_parameters_model( + x_data, + y_data, + x_joints=None, + end_point_constraint=True, + objective_form=False, + name="cubic spline parameters model", +): + """Create a Pyomo model to calculate parameters for a cubic spline. By default + this creates a square linear model, but optionally it can leave off the endpoint + second derivative constraints and add an objective function for fitting data + instead. The purpose of alternative objective form is to allow the spline to + be constrained in other way and/or extend the spline for a particular + extrapolation behavior. + + Args: + x_data: list of x data + y_data: list of y data + x_joints: optional list of joints + end_point_constraint: if true add constraint that second derivative = 0 at + endpoints + objective_form: if true write a least squares objective rather than constraints + to match data + name: optional model name + + + Returns: + Poymo ConcreteModel + """ + n_data = len(x_data) + assert n_data == len(y_data) + if x_joints is None: + n_joints = n_data + n_seg = n_data - 1 + x_joints = x_data + else: + n_joints = len(x_joints) + n_seg = n_joints - 1 + + m = pyo.ConcreteModel(name=name) + # Sets of indexes for joints, segments, and data + m.jnt_idx = pyo.RangeSet(n_joints) + m.seg_idx = pyo.RangeSet(n_joints - 1) + m.dat_idx = pyo.RangeSet(n_data) + + m.x_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(x_data)}) + m.y_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(y_data)}) + m.x = pyo.Param(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_joints)}) + m.alpha = pyo.Var(m.seg_idx, {1, 2, 3, 4}, initialize=1) + + # f_s(x) = f_s+1(x) + @m.Constraint(m.seg_idx) + def y_eqn(blk, s): + if s == m.seg_idx.last(): + return pyo.Constraint.Skip + return f(m.x[s + 1], m.alpha, s) == f(m.x[s + 1], m.alpha, s + 1) + + # f'_s(x) = f'_s+1(x) + @m.Constraint(m.seg_idx) + def yx_eqn(blk, s): + if s == m.seg_idx.last(): + return pyo.Constraint.Skip + return fx(m.x[s + 1], m.alpha, s) == fx(m.x[s + 1], m.alpha, s + 1) + + # f"_s(x) = f"_s+1(x) + @m.Constraint(m.seg_idx) + def yxx_eqn(blk, s): + if s == m.seg_idx.last(): + return pyo.Constraint.Skip + return fxx(m.x[s + 1], m.alpha, s) == fxx(m.x[s + 1], m.alpha, s + 1) + + # Objective function + + # Identify segments used to predict y_data at each x_data. We use search in + # instead of a dict lookup, since we don't want to require the data to be at + # the joints, even though that is almost always the case. + idx = np.searchsorted(x_joints, x_data) + + if end_point_constraint: + add_endpoint_second_derivative_constraints(m) + + # Expression for difference between data and prediction + @m.Expression(m.dat_idx) + def ydiff(blk, d): + s = idx[d - 1] + 1 + if s >= m.seg_idx.last(): + s -= 1 + return m.y_data[d] - f(m.x_data[d], m.alpha, s) + + if objective_form: + # least squares objective + m.obj = pyo.Objective(expr=sum(m.ydiff[d] ** 2 for d in m.dat_idx)) + else: + + @m.Constraint(m.dat_idx) + def match_data(blk, d): + return m.ydiff[d] == 0 + + return m + + +def add_endpoint_second_derivative_constraints(m): + """Usually cubic splines use the endpoint constraints that the second + derivative is zero. This function adds those constraints to a model + """ + + @m.Constraint([m.seg_idx.first(), m.seg_idx.last()]) + def yxx_endpoint_eqn(blk, s): + if s == m.seg_idx.last(): + j = m.jnt_idx.last() + else: + j = s + return fxx(m.x[j], m.alpha, s) == 0 + + +def get_parameters(m, file_name=None): + + joints = [pyo.value(x) for x in m.x.values()] + a1 = [None] * len(m.seg_idx) + a2 = [None] * len(m.seg_idx) + a3 = [None] * len(m.seg_idx) + a4 = [None] * len(m.seg_idx) + for s in m.seg_idx: + a1[s - 1] = pyo.value(m.alpha[s, 1]) + a2[s - 1] = pyo.value(m.alpha[s, 2]) + a3[s - 1] = pyo.value(m.alpha[s, 3]) + a4[s - 1] = pyo.value(m.alpha[s, 4]) + + if file_name is not None: + with open(file_name, "w") as fptr: + fptr.write(f"{len(a1)}\n") + for l in [joints, a1, a2, a3, a4]: + for x in l: + fptr.write(f"{x}\n") + + return joints, a1, a2, a3, a4 + +def _extract_params(m): + """Extract alpha as a plain dict of floats to play nice with vectorized functions""" + alpha = {} + for s in m.seg_idx: + alpha[s] = {} + alpha[s][1] = pyo.value(m.alpha[s, 1]) + alpha[s][2] = pyo.value(m.alpha[s, 2]) + alpha[s][3] = pyo.value(m.alpha[s, 3]) + alpha[s][4] = pyo.value(m.alpha[s, 4]) + return alpha + +def plot_f(m, file_name=None, **kwargs): + """Plot the cspline function. + + Args: + m: Pyomo model with data and parameters + file_name: optional file to save plot to + + Returns: + pyplot object + """ + if not plt_available: + raise ModuleNotFoundError("Matplotlib is not available") + plt.close() + alpha = _extract_params(m) + for s in m.seg_idx: + xvec = np.linspace(m.x[s], m.x[s + 1], 20) + plt.plot(xvec, f(xvec, alpha[s])) + plt.title("f(x)") + x = [] + y = [] + for i in m.dat_idx: + x.append(pyo.value(m.x_data[i])) + y.append(pyo.value(m.y_data[i])) + plt.scatter(x, y) + if file_name is not None: + plt.savefig(file_name, **kwargs) + return plt + +def plot_fx(m, file_name=None, **kwargs): + """Plot the cspline derivative function. + + Args: + m: Pyomo model with data and parameters + file_name: optional file to save plot to + + Returns: + pyplot object + """ + if not plt_available: + raise ModuleNotFoundError("Matplotlib is not available") + plt.close() + alpha = _extract_params(m) + for s in m.seg_idx: + xvec = np.linspace(m.x[s], m.x[s + 1], 20) + plt.plot(xvec, fx(xvec, alpha[s])) + plt.title("f'(x)") + if file_name is not None: + plt.savefig(file_name, **kwargs) + return plt + + +def plot_fxx(m, file_name=None, **kwargs): + """Plot the cspline second derivative function. + + Args: + m: Pyomo model with data and parameters + file_name: optional file to save plot to + + Returns: + pyplot object + """ + if not plt_available: + raise ModuleNotFoundError("Matplotlib is not available") + plt.close() + alpha = _extract_params(m) + for s in m.seg_idx: + xvec = np.linspace(m.x[s], m.x[s + 1], 20) + plt.plot(xvec, fxx(xvec, alpha[s])) + if file_name is not None: + plt.savefig(file_name, **kwargs) + return plt + + +if __name__ == "__main__": + + x_joints = [ + # -2, + # -1, + 1, + 2, + 3, + 4, + 5, + # 7, + # 8, + ] + + x_data = [ + 1, + 2, + 3, + 4, + 5, + ] + + y_data = [ + 2, + 3, + 5, + 2, + 1, + ] + + m = cubic_parameters_model(x_data, y_data, x_joints) + + """ + m.alpha[1, 1].fix(0) + m.alpha[1, 2].fix(0) + m.alpha[1, 3].fix(0) + m.alpha[1, 4].fix(0) + + m.alpha[8, 1].fix(0) + m.alpha[8, 2].fix(0) + m.alpha[8, 3].fix(0) + m.alpha[8, 4].fix(0) + + + @m.Constraint(m.jnt_idx) + def c1(blk, j): + if j > m.seg_idx.last(): + s = j - 1 + else: + s = j + if s >= m.jnt_idx.last() - 2: + return pyo.Constraint.Skip + elif s <= 3: + return pyo.Constraint.Skip + return fxx(m.x[j], m.alpha, s) <= 0 + """ + + + solver_obj = pyo.SolverFactory("clp") + solver_obj.solve(m, tee=True) + + get_parameters(m, "test.txt") + + p = plot_f(m, file_name="f_plot.pdf") + p = plot_fx(m, file_name="fx_plot.pdf") + p = plot_fxx(m, file_name="fxx_plot.pdf") + diff --git a/pyomo/contrib/cspline_1d_external/plugins.py b/pyomo/contrib/cspline_1d_external/plugins.py new file mode 100644 index 00000000000..742fd6afdd4 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/plugins.py @@ -0,0 +1,19 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +from pyomo.common.extensions import ExtensionBuilderFactory +from pyomo.contrib.cspline_1d_external.build import AMPLCspline1DExternalBuilder + + +def load(): + ExtensionBuilderFactory.register("cspline_1d_external")( + AMPLCspline1DExternalBuilder + ) diff --git a/pyomo/contrib/cspline_1d_external/src/CMakeLists.txt b/pyomo/contrib/cspline_1d_external/src/CMakeLists.txt new file mode 100644 index 00000000000..eacf1073372 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/src/CMakeLists.txt @@ -0,0 +1,43 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +cmake_minimum_required(VERSION 3.0) +# CMake 3.0 needed by AMPL/asl build + +PROJECT( ampl_function_demo ) + +INCLUDE( + "${CMAKE_CURRENT_SOURCE_DIR}/../../cspline_1d_external/src/FindASL.cmake") + +# Targets in this project +OPTION(BUILD_EXTERNAL_FCN_LIBRARY + "Build the ASL external function example library" ON) + +IF( BUILD_EXTERNAL_FCN_LIBRARY ) + ADD_LIBRARY( cspline_1d_external SHARED "functions.cpp" ) + TARGET_LINK_LIBRARIES( cspline_1d_external + PUBLIC ${ASL_LIBRARY} ${CMAKE_DL_LIBS}) + TARGET_INCLUDE_DIRECTORIES( cspline_1d_external + PUBLIC ${ASL_INCLUDE_DIR} + INTERFACE . ) + # If you need a CPP directive defined when building the library (e.g., + # for managing __declspec(dllimport) under Windows, uncomment the + # following: + #TARGET_COMPILE_DEFINITIONS( cspline_1d_external PRIVATE BUILDING_ASL_DEMO ) + #SET_TARGET_PROPERTIES( cspline_1d_external PROPERTIES ENABLE_EXPORTS 1 ) + INSTALL( TARGETS cspline_1d_external LIBRARY DESTINATION lib + RUNTIME DESTINATION lib ) + IF( BUILD_AMPLASL ) + # If we are building AMPL/asl (from FindASL), it is possible that we + # are linking against it, so we will add the appropriate dependency + add_dependencies(cspline_1d_external ampl_asl) + ENDIF() +ENDIF() diff --git a/pyomo/contrib/cspline_1d_external/src/FindASL.cmake b/pyomo/contrib/cspline_1d_external/src/FindASL.cmake new file mode 100644 index 00000000000..f413176f1cc --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/src/FindASL.cmake @@ -0,0 +1,106 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +cmake_minimum_required(VERSION 3.0) +# Minimum version inherited from AMPL/asl + +include(ExternalProject) + +# Dependencies that we manage / can install +SET(AMPLASL_TAG "9fb7cb8e4f68ed1c3bc066d191e63698b7d7d1d2" CACHE STRING + "AMPL/asl git tag/branch to checkout and build") +# 9fb7cb8e4f68ed1c3bc066d191e63698b7d7d1d2 corresponds to ASLdate = 20211109 +OPTION(BUILD_AMPLASL + "Download and build AMPL/asl ${AMPLASL_TAG} from GitHub" OFF) + +# Other build / environment options +OPTION(BUILD_AMPLASL_IF_NEEDED + "Automatically enable AMPL/asl build if ASL not found" OFF) +MARK_AS_ADVANCED(BUILD_AMPLASL_IF_NEEDED) + +OPTION(ASL_USE_PKG_CONFIG, + "Use pkgconfig (if present) to attempt to locate the ASL" OFF) +#OPTION(STATIC_LINK "STATIC_LINK" OFF) + +# We need the ASL. We can get it from Ipopt, AMPL/asl, or ASL (netlib) +SET(IPOPT_DIR "" CACHE PATH "Path to compiled Ipopt installation") +SET(AMPLASL_DIR "" CACHE PATH "Path to compiled AMPL/asl installation") +#SET(ASL_NETLIB_DIR "" CACHE PATH "Path to compiled ASL (netlib) installation") + +# Use pkg-config to get the ASL directories from the Ipopt/COIN-OR build +FIND_PACKAGE(PkgConfig) +IF( PKG_CONFIG_FOUND AND ASL_USE_PKG_CONFIG ) + SET(_TMP "$ENV{PKG_CONFIG_PATH}") + SET(ENV{PKG_CONFIG_PATH} "${IPOPT_DIR}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}") + pkg_check_modules(PC_COINASL QUIET coinasl) + SET(ENV{PKG_CONFIG_PATH} "${_TMP}") +ENDIF() + +# cmake does not search LD_LIBRARY_PATH by default. So that libraries +# can be added through mechanisms like 'environment modules', we will explicitly +# add LD_LIBRARY_PATH to the search path +string(REPLACE ":" ";" LD_LIBRARY_DIR_LIST + $ENV{LD_LIBRARY_PATH}:$ENV{DYLD_LIBRARY_PATH} + ) + +# Note: the directory search order is intentional: first the modules we +# are creating, then directories specifically set by the user, and +# finally automatically located installations (e.g., from pkg-config) +FIND_PATH(ASL_INCLUDE_DIR asl_pfgh.h + HINTS "${CMAKE_INSTALL_PREFIX}/include/asl" + "${IPOPT_DIR}/include/coin-or/asl" + "${IPOPT_DIR}/include/coin/ThirdParty" + "${AMPLASL_DIR}/include/asl" + "${PC_COINASL_INCLUDEDIR}" + "${PC_COINASL_INCLUDE_DIRS}" + PATH_SUFFIXES asl +) +FIND_LIBRARY(ASL_LIBRARY NAMES asl coinasl + HINTS "${CMAKE_INSTALL_PREFIX}/lib" + "${IPOPT_DIR}/lib" + "${AMPLASL_DIR}/lib" + "${PC_COINASL_LIBDIR}" + "${PC_COINASL_LIBRARY_DIRS}" + ${LD_LIBRARY_DIR_LIST} +) + +# If BUILD_AMPLASL_IF_NEEDED is set and we couldn't find / weren't +# pointed to an ASL build, then we will forcibly enable the AMPL/asl build +# to provide the ASL. +IF( BUILD_AMPLASL_IF_NEEDED AND (NOT ASL_LIBRARY OR NOT ASL_INCLUDE_DIR) ) + set_property(CACHE BUILD_AMPLASL PROPERTY VALUE ON) +ENDIF() + +IF( BUILD_AMPLASL ) + get_filename_component(ABS_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}" ABSOLUTE) + ExternalProject_Add(ampl_asl + GIT_TAG ${AMPLASL_TAG} + GIT_REPOSITORY https://github.com/ampl/asl.git + CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:STRING=${ABS_INSTALL_PREFIX} + UPDATE_DISCONNECTED TRUE + ) + # Update the ASL paths (if necessary). Since these do not (yet) + # exist, we need to bypass find_path / find_library and explicitly set + # the directories that this build will create. However, we will only + # do this if the paths have not already been set (so users can always + # override what we do here) + IF(NOT ASL_INCLUDE_DIR OR NOT ASL_LIBRARY) + set_property(CACHE ASL_INCLUDE_DIR PROPERTY VALUE + "${ABS_INSTALL_PREFIX}/include/asl") + IF( WIN32 ) + set_property(CACHE ASL_LIBRARY PROPERTY VALUE + "${ABS_INSTALL_PREFIX}/lib/asl.lib") + ELSE() + set_property(CACHE ASL_LIBRARY PROPERTY VALUE + "${ABS_INSTALL_PREFIX}/lib/libasl.a") + ENDIF() + ENDIF() +ENDIF() diff --git a/pyomo/contrib/cspline_1d_external/src/functions.cpp b/pyomo/contrib/cspline_1d_external/src/functions.cpp new file mode 100644 index 00000000000..e65f9b01245 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/src/functions.cpp @@ -0,0 +1,144 @@ +/* ___________________________________________________________________________ + * Pyomo: Python Optimization Modeling Objects + * Copyright (c) 2008-2022 + * National Technology and Engineering Solutions of Sandia, LLC + * Under the terms of Contract DE-NA0003525 with National Technology and + * Engineering Solutions of Sandia, LLC, the U.S. Government retains certain + * rights in this software. + * This software is distributed under the 3-clause BSD License. + * ___________________________________________________________________________ +*/ +#include +#include +#include +#include +#include + +#include "funcadd.h" + +unsigned int n_cspline_1d = 0; +std::unordered_map idx_cspline_1d; +std::vector< std::vector > joints_cspline_1d; +std::vector< std::vector > a1_cspline_1d; +std::vector< std::vector > a2_cspline_1d; +std::vector< std::vector > a3_cspline_1d; +std::vector< std::vector > a4_cspline_1d; +std::vector< unsigned int > n_seg_cspline_1d; + + +int read_parameters_cspline_1d(std::string file_path) { + // Read data table + try{ + // Parameters have been calculated already + return idx_cspline_1d.at(file_path); + } + catch(std::out_of_range const&){ + // Parameters haven't been read yet, so do that. + } + unsigned int i, idx; + unsigned int n; // number of segments + std::ifstream param_file; // parameter input file + + // Set index for cspline param file + idx_cspline_1d[file_path] = n_cspline_1d; + idx = n_cspline_1d; + ++n_cspline_1d; + + // open the parameter file + param_file.open(file_path); + + // get the number of segments + param_file >> n; + n_seg_cspline_1d.resize(n_cspline_1d); + n_seg_cspline_1d[idx] = n; + + joints_cspline_1d.resize(n_cspline_1d); + a1_cspline_1d.resize(n_cspline_1d); + a2_cspline_1d.resize(n_cspline_1d); + a3_cspline_1d.resize(n_cspline_1d); + a4_cspline_1d.resize(n_cspline_1d); + + joints_cspline_1d[idx].resize(n + 1); + a1_cspline_1d[idx].resize(n); + a2_cspline_1d[idx].resize(n); + a3_cspline_1d[idx].resize(n); + a4_cspline_1d[idx].resize(n); + + // get the joints + for(i=0; i < n + 1; ++i){ + param_file >> joints_cspline_1d[idx][i]; + } + + // get the a1 params + for(i=0; i < n; ++i){ + param_file >> a1_cspline_1d[idx][i]; + } + + // get the a2 params + for(i=0; i < n; ++i){ + param_file >> a2_cspline_1d[idx][i]; + } + + // get the a3 params + for(i=0; i < n; ++i){ + param_file >> a3_cspline_1d[idx][i]; + } + + // get the a4 params + for(i=0; i < n; ++i){ + param_file >> a4_cspline_1d[idx][i]; + } + + param_file.close(); + return idx; +} + +extern real cspline_1d(arglist *al) { + const char* data_file = al->sa[-(al->at[1]+1)]; + real x = al->ra[al->at[0]]; + unsigned int idx = read_parameters_cspline_1d(data_file); + size_t seg; + real a1, a2, a3, a4; + + //find segment index + auto lit = std::lower_bound(joints_cspline_1d[idx].begin(), joints_cspline_1d[idx].end(), x); + seg = lit - joints_cspline_1d[idx].begin(); + if(seg > 0) --seg; + if(seg >= n_seg_cspline_1d[idx]) seg = n_seg_cspline_1d[idx] - 1; + + a1 = a1_cspline_1d[idx][seg]; + a2 = a2_cspline_1d[idx][seg]; + a3 = a3_cspline_1d[idx][seg]; + a4 = a4_cspline_1d[idx][seg]; + + // Compute the first derivative, if requested. + if (al->derivs!=NULL) { + al->derivs[0] = a2 + 2*a3*x + 3*a4*x*x; + // Compute the second derivative, if requested. + if (al->hes!=NULL) { + al->hes[0] = 2*a3 + 6*a4*x; + } + } + + return a1 + a2*x + a3*x*x + a4*x*x*x; +} + +// Register external functions defined in this library with the ASL +void funcadd(AmplExports *ae){ + // Arguments for addfunc (this is not fully detailed; see funcadd.h) + // 1) Name of function in AMPL + // 2) Function pointer to C function + // 3) see FUNCADD_TYPE enum in funcadd.h + // 4) Number of arguments + // >=0 indicates a fixed number of arguments + // < 0 indicates a variable length list (requiring at least -(n+1) + // arguments) + // 5) Void pointer to function info + addfunc( + "cspline_1d", + (rfunc)cspline_1d, + FUNCADD_REAL_VALUED|FUNCADD_STRING_ARGS, + 2, + NULL + ); +} diff --git a/pyomo/contrib/cspline_1d_external/test/__init__.py b/pyomo/contrib/cspline_1d_external/test/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/pyomo/contrib/cspline_1d_external/test/function.py b/pyomo/contrib/cspline_1d_external/test/function.py new file mode 100644 index 00000000000..77caa046b47 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/test/function.py @@ -0,0 +1,38 @@ +import pyomo.environ as pyo +from pyomo.common.fileutils import find_library + +if __name__ == "__main__": + lib = find_library("cspline_1d_external") + params = "t1_params.txt" + + m = pyo.ConcreteModel() + m.f = pyo.ExternalFunction(library=lib, function="cspline_1d") + + m.x = pyo.Var(initialize=2) # , bounds=(0.9, 5.1)) + m.y = pyo.Var() + + m.c1 = pyo.Constraint(expr=m.y == m.f(m.x, params)) + + m.o = pyo.Objective(expr=-m.y) + + print(pyo.value(m.f(0.9, params))) + print(pyo.value(m.f(1.01, params))) + print(pyo.value(m.f(1.05, params))) + print(pyo.value(m.f(1.2, params))) + print(pyo.value(m.f(1.3, params))) + print(pyo.value(m.f(1.5, params))) + print(pyo.value(m.f(2, params))) + print(pyo.value(m.f(2.01, params))) + print(pyo.value(m.f(2.1, params))) + print(pyo.value(m.f(3.0, params))) + print(pyo.value(m.f(3.1, params))) + print(pyo.value(m.f(4.9, params))) + print(pyo.value(m.f(5.0, params))) + print(pyo.value(m.f(5.1, params))) + + import idaes + + solver_obj = pyo.SolverFactory("ipopt") + solver_obj.solve(m, tee=True) + + m.display() diff --git a/pyomo/contrib/cspline_1d_external/test/t1_params.txt b/pyomo/contrib/cspline_1d_external/test/t1_params.txt new file mode 100644 index 00000000000..bc6913c6b57 --- /dev/null +++ b/pyomo/contrib/cspline_1d_external/test/t1_params.txt @@ -0,0 +1,22 @@ +4 +1 +2 +3 +4 +5 +0.9999999999999978 +24.714285714286376 +-106.42857142857616 +115.28571428572297 +2.321428571428605 +-33.250000000000966 +97.89285714286154 +-68.39285714286291 +-1.9821428571429038 +15.803571428571882 +-27.91071428571562 +13.660714285715482 +0.6607142857143012 +-2.3035714285714963 +2.55357142857156 +-0.9107142857143654 From d70c1ee48ee2b380337b8802dda726a73306c7ef Mon Sep 17 00:00:00 2001 From: Eslick Date: Mon, 8 Apr 2024 12:38:15 -0400 Subject: [PATCH 02/13] Rename some things, make knots variables. --- .../__init__.py | 0 .../{cspline_1d_external => cspline_external}/build.py | 0 .../cspline_parameters.py} | 9 +++++---- .../{cspline_1d_external => cspline_external}/plugins.py | 0 .../src/CMakeLists.txt | 0 .../src/FindASL.cmake | 0 .../src/functions.cpp | 0 .../test/__init__.py | 0 .../test/function.py | 0 .../test/t1_params.txt | 0 10 files changed, 5 insertions(+), 4 deletions(-) rename pyomo/contrib/{cspline_1d_external => cspline_external}/__init__.py (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/build.py (100%) rename pyomo/contrib/{cspline_1d_external/cubic_spline_parameters.py => cspline_external/cspline_parameters.py} (96%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/plugins.py (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/src/CMakeLists.txt (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/src/FindASL.cmake (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/src/functions.cpp (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/test/__init__.py (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/test/function.py (100%) rename pyomo/contrib/{cspline_1d_external => cspline_external}/test/t1_params.txt (100%) diff --git a/pyomo/contrib/cspline_1d_external/__init__.py b/pyomo/contrib/cspline_external/__init__.py similarity index 100% rename from pyomo/contrib/cspline_1d_external/__init__.py rename to pyomo/contrib/cspline_external/__init__.py diff --git a/pyomo/contrib/cspline_1d_external/build.py b/pyomo/contrib/cspline_external/build.py similarity index 100% rename from pyomo/contrib/cspline_1d_external/build.py rename to pyomo/contrib/cspline_external/build.py diff --git a/pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py similarity index 96% rename from pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py rename to pyomo/contrib/cspline_external/cspline_parameters.py index 472858b8061..671cc94a6d2 100644 --- a/pyomo/contrib/cspline_1d_external/cubic_spline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -123,7 +123,8 @@ def cubic_parameters_model( m.x_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(x_data)}) m.y_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(y_data)}) - m.x = pyo.Param(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_joints)}) + m.x = pyo.Var(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_joints)}) + m.x.fix() m.alpha = pyo.Var(m.seg_idx, {1, 2, 3, 4}, initialize=1) # f_s(x) = f_s+1(x) @@ -239,7 +240,7 @@ def plot_f(m, file_name=None, **kwargs): plt.close() alpha = _extract_params(m) for s in m.seg_idx: - xvec = np.linspace(m.x[s], m.x[s + 1], 20) + xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) plt.plot(xvec, f(xvec, alpha[s])) plt.title("f(x)") x = [] @@ -267,7 +268,7 @@ def plot_fx(m, file_name=None, **kwargs): plt.close() alpha = _extract_params(m) for s in m.seg_idx: - xvec = np.linspace(m.x[s], m.x[s + 1], 20) + xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) plt.plot(xvec, fx(xvec, alpha[s])) plt.title("f'(x)") if file_name is not None: @@ -290,7 +291,7 @@ def plot_fxx(m, file_name=None, **kwargs): plt.close() alpha = _extract_params(m) for s in m.seg_idx: - xvec = np.linspace(m.x[s], m.x[s + 1], 20) + xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) plt.plot(xvec, fxx(xvec, alpha[s])) if file_name is not None: plt.savefig(file_name, **kwargs) diff --git a/pyomo/contrib/cspline_1d_external/plugins.py b/pyomo/contrib/cspline_external/plugins.py similarity index 100% rename from pyomo/contrib/cspline_1d_external/plugins.py rename to pyomo/contrib/cspline_external/plugins.py diff --git a/pyomo/contrib/cspline_1d_external/src/CMakeLists.txt b/pyomo/contrib/cspline_external/src/CMakeLists.txt similarity index 100% rename from pyomo/contrib/cspline_1d_external/src/CMakeLists.txt rename to pyomo/contrib/cspline_external/src/CMakeLists.txt diff --git a/pyomo/contrib/cspline_1d_external/src/FindASL.cmake b/pyomo/contrib/cspline_external/src/FindASL.cmake similarity index 100% rename from pyomo/contrib/cspline_1d_external/src/FindASL.cmake rename to pyomo/contrib/cspline_external/src/FindASL.cmake diff --git a/pyomo/contrib/cspline_1d_external/src/functions.cpp b/pyomo/contrib/cspline_external/src/functions.cpp similarity index 100% rename from pyomo/contrib/cspline_1d_external/src/functions.cpp rename to pyomo/contrib/cspline_external/src/functions.cpp diff --git a/pyomo/contrib/cspline_1d_external/test/__init__.py b/pyomo/contrib/cspline_external/test/__init__.py similarity index 100% rename from pyomo/contrib/cspline_1d_external/test/__init__.py rename to pyomo/contrib/cspline_external/test/__init__.py diff --git a/pyomo/contrib/cspline_1d_external/test/function.py b/pyomo/contrib/cspline_external/test/function.py similarity index 100% rename from pyomo/contrib/cspline_1d_external/test/function.py rename to pyomo/contrib/cspline_external/test/function.py diff --git a/pyomo/contrib/cspline_1d_external/test/t1_params.txt b/pyomo/contrib/cspline_external/test/t1_params.txt similarity index 100% rename from pyomo/contrib/cspline_1d_external/test/t1_params.txt rename to pyomo/contrib/cspline_external/test/t1_params.txt From 7663630ca2c4d1e02eaed8fb54c948ab4964d689 Mon Sep 17 00:00:00 2001 From: Eslick Date: Mon, 8 Apr 2024 12:45:41 -0400 Subject: [PATCH 03/13] Fix some names --- .../cspline_external/cspline_parameters.py | 36 +++++++++---------- .../cspline_external/src/CMakeLists.txt | 16 ++++----- .../cspline_external/src/functions.cpp | 14 ++++---- .../contrib/cspline_external/test/function.py | 2 +- 4 files changed, 34 insertions(+), 34 deletions(-) diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index 671cc94a6d2..dbe3970a530 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -79,7 +79,7 @@ def fxx(x, alpha, s=None): def cubic_parameters_model( x_data, y_data, - x_joints=None, + x_knots=None, end_point_constraint=True, objective_form=False, name="cubic spline parameters model", @@ -94,7 +94,7 @@ def cubic_parameters_model( Args: x_data: list of x data y_data: list of y data - x_joints: optional list of joints + x_knots: optional list of knots end_point_constraint: if true add constraint that second derivative = 0 at endpoints objective_form: if true write a least squares objective rather than constraints @@ -107,23 +107,23 @@ def cubic_parameters_model( """ n_data = len(x_data) assert n_data == len(y_data) - if x_joints is None: - n_joints = n_data + if x_knots is None: + n_knots = n_data n_seg = n_data - 1 - x_joints = x_data + x_knots = x_data else: - n_joints = len(x_joints) - n_seg = n_joints - 1 + n_knots = len(x_knots) + n_seg = n_knots - 1 m = pyo.ConcreteModel(name=name) - # Sets of indexes for joints, segments, and data - m.jnt_idx = pyo.RangeSet(n_joints) - m.seg_idx = pyo.RangeSet(n_joints - 1) + # Sets of indexes for knots, segments, and data + m.jnt_idx = pyo.RangeSet(n_knots) + m.seg_idx = pyo.RangeSet(n_knots - 1) m.dat_idx = pyo.RangeSet(n_data) m.x_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(x_data)}) m.y_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(y_data)}) - m.x = pyo.Var(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_joints)}) + m.x = pyo.Var(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_knots)}) m.x.fix() m.alpha = pyo.Var(m.seg_idx, {1, 2, 3, 4}, initialize=1) @@ -152,8 +152,8 @@ def yxx_eqn(blk, s): # Identify segments used to predict y_data at each x_data. We use search in # instead of a dict lookup, since we don't want to require the data to be at - # the joints, even though that is almost always the case. - idx = np.searchsorted(x_joints, x_data) + # the knots, even though that is almost always the case. + idx = np.searchsorted(x_knots, x_data) if end_point_constraint: add_endpoint_second_derivative_constraints(m) @@ -194,7 +194,7 @@ def yxx_endpoint_eqn(blk, s): def get_parameters(m, file_name=None): - joints = [pyo.value(x) for x in m.x.values()] + knots = [pyo.value(x) for x in m.x.values()] a1 = [None] * len(m.seg_idx) a2 = [None] * len(m.seg_idx) a3 = [None] * len(m.seg_idx) @@ -208,11 +208,11 @@ def get_parameters(m, file_name=None): if file_name is not None: with open(file_name, "w") as fptr: fptr.write(f"{len(a1)}\n") - for l in [joints, a1, a2, a3, a4]: + for l in [knots, a1, a2, a3, a4]: for x in l: fptr.write(f"{x}\n") - return joints, a1, a2, a3, a4 + return knots, a1, a2, a3, a4 def _extract_params(m): """Extract alpha as a plain dict of floats to play nice with vectorized functions""" @@ -300,7 +300,7 @@ def plot_fxx(m, file_name=None, **kwargs): if __name__ == "__main__": - x_joints = [ + x_knots = [ # -2, # -1, 1, @@ -328,7 +328,7 @@ def plot_fxx(m, file_name=None, **kwargs): 1, ] - m = cubic_parameters_model(x_data, y_data, x_joints) + m = cubic_parameters_model(x_data, y_data, x_knots) """ m.alpha[1, 1].fix(0) diff --git a/pyomo/contrib/cspline_external/src/CMakeLists.txt b/pyomo/contrib/cspline_external/src/CMakeLists.txt index eacf1073372..b710bc948b6 100644 --- a/pyomo/contrib/cspline_external/src/CMakeLists.txt +++ b/pyomo/contrib/cspline_external/src/CMakeLists.txt @@ -15,29 +15,29 @@ cmake_minimum_required(VERSION 3.0) PROJECT( ampl_function_demo ) INCLUDE( - "${CMAKE_CURRENT_SOURCE_DIR}/../../cspline_1d_external/src/FindASL.cmake") + "${CMAKE_CURRENT_SOURCE_DIR}/../../cspline_external/src/FindASL.cmake") # Targets in this project OPTION(BUILD_EXTERNAL_FCN_LIBRARY "Build the ASL external function example library" ON) IF( BUILD_EXTERNAL_FCN_LIBRARY ) - ADD_LIBRARY( cspline_1d_external SHARED "functions.cpp" ) - TARGET_LINK_LIBRARIES( cspline_1d_external + ADD_LIBRARY( cspline_external SHARED "functions.cpp" ) + TARGET_LINK_LIBRARIES( cspline_external PUBLIC ${ASL_LIBRARY} ${CMAKE_DL_LIBS}) - TARGET_INCLUDE_DIRECTORIES( cspline_1d_external + TARGET_INCLUDE_DIRECTORIES( cspline_external PUBLIC ${ASL_INCLUDE_DIR} INTERFACE . ) # If you need a CPP directive defined when building the library (e.g., # for managing __declspec(dllimport) under Windows, uncomment the # following: - #TARGET_COMPILE_DEFINITIONS( cspline_1d_external PRIVATE BUILDING_ASL_DEMO ) - #SET_TARGET_PROPERTIES( cspline_1d_external PROPERTIES ENABLE_EXPORTS 1 ) - INSTALL( TARGETS cspline_1d_external LIBRARY DESTINATION lib + #TARGET_COMPILE_DEFINITIONS( cspline_external PRIVATE BUILDING_ASL_DEMO ) + #SET_TARGET_PROPERTIES( cspline_external PROPERTIES ENABLE_EXPORTS 1 ) + INSTALL( TARGETS cspline_external LIBRARY DESTINATION lib RUNTIME DESTINATION lib ) IF( BUILD_AMPLASL ) # If we are building AMPL/asl (from FindASL), it is possible that we # are linking against it, so we will add the appropriate dependency - add_dependencies(cspline_1d_external ampl_asl) + add_dependencies(cspline_external ampl_asl) ENDIF() ENDIF() diff --git a/pyomo/contrib/cspline_external/src/functions.cpp b/pyomo/contrib/cspline_external/src/functions.cpp index e65f9b01245..52e9415f97a 100644 --- a/pyomo/contrib/cspline_external/src/functions.cpp +++ b/pyomo/contrib/cspline_external/src/functions.cpp @@ -18,7 +18,7 @@ unsigned int n_cspline_1d = 0; std::unordered_map idx_cspline_1d; -std::vector< std::vector > joints_cspline_1d; +std::vector< std::vector > knots_cspline_1d; std::vector< std::vector > a1_cspline_1d; std::vector< std::vector > a2_cspline_1d; std::vector< std::vector > a3_cspline_1d; @@ -52,21 +52,21 @@ int read_parameters_cspline_1d(std::string file_path) { n_seg_cspline_1d.resize(n_cspline_1d); n_seg_cspline_1d[idx] = n; - joints_cspline_1d.resize(n_cspline_1d); + knots_cspline_1d.resize(n_cspline_1d); a1_cspline_1d.resize(n_cspline_1d); a2_cspline_1d.resize(n_cspline_1d); a3_cspline_1d.resize(n_cspline_1d); a4_cspline_1d.resize(n_cspline_1d); - joints_cspline_1d[idx].resize(n + 1); + knots_cspline_1d[idx].resize(n + 1); a1_cspline_1d[idx].resize(n); a2_cspline_1d[idx].resize(n); a3_cspline_1d[idx].resize(n); a4_cspline_1d[idx].resize(n); - // get the joints + // get the knots for(i=0; i < n + 1; ++i){ - param_file >> joints_cspline_1d[idx][i]; + param_file >> knots_cspline_1d[idx][i]; } // get the a1 params @@ -101,8 +101,8 @@ extern real cspline_1d(arglist *al) { real a1, a2, a3, a4; //find segment index - auto lit = std::lower_bound(joints_cspline_1d[idx].begin(), joints_cspline_1d[idx].end(), x); - seg = lit - joints_cspline_1d[idx].begin(); + auto lit = std::lower_bound(knots_cspline_1d[idx].begin(), knots_cspline_1d[idx].end(), x); + seg = lit - knots_cspline_1d[idx].begin(); if(seg > 0) --seg; if(seg >= n_seg_cspline_1d[idx]) seg = n_seg_cspline_1d[idx] - 1; diff --git a/pyomo/contrib/cspline_external/test/function.py b/pyomo/contrib/cspline_external/test/function.py index 77caa046b47..462aa0a29ee 100644 --- a/pyomo/contrib/cspline_external/test/function.py +++ b/pyomo/contrib/cspline_external/test/function.py @@ -2,7 +2,7 @@ from pyomo.common.fileutils import find_library if __name__ == "__main__": - lib = find_library("cspline_1d_external") + lib = find_library("cspline_external") params = "t1_params.txt" m = pyo.ConcreteModel() From 16406285ee7df7da051a8afaf87b763f43adb3cb Mon Sep 17 00:00:00 2001 From: Eslick Date: Wed, 10 Apr 2024 16:46:40 -0400 Subject: [PATCH 04/13] Add some other parameters --- .../cspline_external/cspline_parameters.py | 49 ++++++--------- .../cspline_external/src/functions.cpp | 61 +++++++++++++------ .../contrib/cspline_external/test/function.py | 4 +- .../cspline_external/test/t2_params.txt | 42 +++++++++++++ 4 files changed, 108 insertions(+), 48 deletions(-) create mode 100644 pyomo/contrib/cspline_external/test/t2_params.txt diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index dbe3970a530..6faaf75c820 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -117,14 +117,13 @@ def cubic_parameters_model( m = pyo.ConcreteModel(name=name) # Sets of indexes for knots, segments, and data - m.jnt_idx = pyo.RangeSet(n_knots) + m.knt_idx = pyo.RangeSet(n_knots) m.seg_idx = pyo.RangeSet(n_knots - 1) m.dat_idx = pyo.RangeSet(n_data) m.x_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(x_data)}) m.y_data = pyo.Param(m.dat_idx, initialize={i + 1: x for i, x in enumerate(y_data)}) - m.x = pyo.Var(m.jnt_idx, initialize={i + 1: x for i, x in enumerate(x_knots)}) - m.x.fix() + m.x = pyo.Param(m.knt_idx, initialize={i + 1: x for i, x in enumerate(x_knots)}) m.alpha = pyo.Var(m.seg_idx, {1, 2, 3, 4}, initialize=1) # f_s(x) = f_s+1(x) @@ -186,7 +185,7 @@ def add_endpoint_second_derivative_constraints(m): @m.Constraint([m.seg_idx.first(), m.seg_idx.last()]) def yxx_endpoint_eqn(blk, s): if s == m.seg_idx.last(): - j = m.jnt_idx.last() + j = m.knt_idx.last() else: j = s return fxx(m.x[j], m.alpha, s) == 0 @@ -301,15 +300,15 @@ def plot_fxx(m, file_name=None, **kwargs): if __name__ == "__main__": x_knots = [ - # -2, - # -1, + 0, 1, 2, + 2.5, 3, + 3.5, 4, 5, - # 7, - # 8, + 6, ] x_data = [ @@ -328,35 +327,27 @@ def plot_fxx(m, file_name=None, **kwargs): 1, ] - m = cubic_parameters_model(x_data, y_data, x_knots) + m = cubic_parameters_model( + x_data, + y_data, + x_knots, + end_point_constraint=False, + objective_form=True, + ) - """ - m.alpha[1, 1].fix(0) - m.alpha[1, 2].fix(0) + #m.alpha[1, 1].fix(0) + m.alpha[1, 2].fix(0.5) m.alpha[1, 3].fix(0) m.alpha[1, 4].fix(0) - m.alpha[8, 1].fix(0) - m.alpha[8, 2].fix(0) + #m.alpha[8, 1].fix(0) + m.alpha[8, 2].fix(-0.5) m.alpha[8, 3].fix(0) m.alpha[8, 4].fix(0) - @m.Constraint(m.jnt_idx) - def c1(blk, j): - if j > m.seg_idx.last(): - s = j - 1 - else: - s = j - if s >= m.jnt_idx.last() - 2: - return pyo.Constraint.Skip - elif s <= 3: - return pyo.Constraint.Skip - return fxx(m.x[j], m.alpha, s) <= 0 - """ - - - solver_obj = pyo.SolverFactory("clp") + #solver_obj = pyo.SolverFactory("clp") + solver_obj = pyo.SolverFactory("ipopt") solver_obj.solve(m, tee=True) get_parameters(m, "test.txt") diff --git a/pyomo/contrib/cspline_external/src/functions.cpp b/pyomo/contrib/cspline_external/src/functions.cpp index 52e9415f97a..51763b83ad0 100644 --- a/pyomo/contrib/cspline_external/src/functions.cpp +++ b/pyomo/contrib/cspline_external/src/functions.cpp @@ -16,48 +16,61 @@ #include "funcadd.h" +// Don't want to have the overhead of reading parameter files +// each time, so keep params as long as library remains loaded + +// Number of loaded curves unsigned int n_cspline_1d = 0; + +// Parameter file path to unsigned in index map std::unordered_map idx_cspline_1d; + +// Knots for each curve std::vector< std::vector > knots_cspline_1d; + +// Cubic segment parameters for each curve std::vector< std::vector > a1_cspline_1d; std::vector< std::vector > a2_cspline_1d; std::vector< std::vector > a3_cspline_1d; std::vector< std::vector > a4_cspline_1d; -std::vector< unsigned int > n_seg_cspline_1d; +// Number of segments in each curve +std::vector< unsigned int > n_seg_cspline_1d; -int read_parameters_cspline_1d(std::string file_path) { - // Read data table - try{ - // Parameters have been calculated already +// The data file format is a text file with one parameter per line +// Line 1: Number of parameters +// nsegs + 1 lines for knots +// nsegs lines for a1 +// nsegs lines for a2 +// nsegs lines for a3 +// nsegs lines for a4 + +// Function to read parameter file if not read and store params +static int read_parameters_cspline_1d(std::string file_path) { + try{ //if read already return curve index and done. return idx_cspline_1d.at(file_path); } - catch(std::out_of_range const&){ - // Parameters haven't been read yet, so do that. - } - unsigned int i, idx; - unsigned int n; // number of segments + catch(std::out_of_range const&){} // Not read so read it. + unsigned int i, idx, n; //loop counter, curve index, and number of segs std::ifstream param_file; // parameter input file - // Set index for cspline param file + // Set index for cspline param file and increment curve count idx_cspline_1d[file_path] = n_cspline_1d; idx = n_cspline_1d; ++n_cspline_1d; - // open the parameter file + // open the parameter file for input param_file.open(file_path); - // get the number of segments + // get the number of segments and size to vectors param_file >> n; n_seg_cspline_1d.resize(n_cspline_1d); n_seg_cspline_1d[idx] = n; - knots_cspline_1d.resize(n_cspline_1d); a1_cspline_1d.resize(n_cspline_1d); a2_cspline_1d.resize(n_cspline_1d); a3_cspline_1d.resize(n_cspline_1d); a4_cspline_1d.resize(n_cspline_1d); - knots_cspline_1d[idx].resize(n + 1); a1_cspline_1d[idx].resize(n); a2_cspline_1d[idx].resize(n); @@ -90,22 +103,36 @@ int read_parameters_cspline_1d(std::string file_path) { } param_file.close(); + + // Returns curve index return idx; } +// Calculate the cspline 1D function value and derivatives extern real cspline_1d(arglist *al) { + // Get the data file path argument const char* data_file = al->sa[-(al->at[1]+1)]; + // Get the function input real x = al->ra[al->at[0]]; + // Get parameters, read file if needed unsigned int idx = read_parameters_cspline_1d(data_file); - size_t seg; - real a1, a2, a3, a4; + size_t seg; // segment index + real a1, a2, a3, a4; // segment parameters //find segment index auto lit = std::lower_bound(knots_cspline_1d[idx].begin(), knots_cspline_1d[idx].end(), x); seg = lit - knots_cspline_1d[idx].begin(); + + // If off the curve on the low side use first seg (0) otherwise + // subtract 1 to go from knot index to segment index (i.e. seg is + // initially how many knots are below, so 0 is less than first knot, + // 1 is above first knot..., so > 0 we need to decrement to get the + // segment index) if(seg > 0) --seg; + // In off the curve on the high side, use last segment if(seg >= n_seg_cspline_1d[idx]) seg = n_seg_cspline_1d[idx] - 1; + // Get the parameters for the correct segment a1 = a1_cspline_1d[idx][seg]; a2 = a2_cspline_1d[idx][seg]; a3 = a3_cspline_1d[idx][seg]; diff --git a/pyomo/contrib/cspline_external/test/function.py b/pyomo/contrib/cspline_external/test/function.py index 462aa0a29ee..2c974f80e54 100644 --- a/pyomo/contrib/cspline_external/test/function.py +++ b/pyomo/contrib/cspline_external/test/function.py @@ -3,12 +3,12 @@ if __name__ == "__main__": lib = find_library("cspline_external") - params = "t1_params.txt" + params = "t2_params.txt" m = pyo.ConcreteModel() m.f = pyo.ExternalFunction(library=lib, function="cspline_1d") - m.x = pyo.Var(initialize=2) # , bounds=(0.9, 5.1)) + m.x = pyo.Var(initialize=2, bounds=(0.5, 5.5)) m.y = pyo.Var() m.c1 = pyo.Constraint(expr=m.y == m.f(m.x, params)) diff --git a/pyomo/contrib/cspline_external/test/t2_params.txt b/pyomo/contrib/cspline_external/test/t2_params.txt new file mode 100644 index 00000000000..7d7609a3987 --- /dev/null +++ b/pyomo/contrib/cspline_external/test/t2_params.txt @@ -0,0 +1,42 @@ +8 +0 +1 +2 +2.5 +3 +3.5 +4 +5 +6 +1.4999999999908438 +0.9999999999698757 +13.000000002737705 +75.49999998081503 +-230.49999994288345 +55.33333326383919 +66.00000001569855 +3.4999999999416516 +0.5 +2.000000000062904 +-16.00000000408884 +-90.99999997778163 +214.99999994591695 +-29.99999994555961 +-38.00000000945413 +-0.5 +0 +-1.500000000062904 +7.500000002012967 +37.49999999149008 +-64.49999998307611 +5.499999985917197 +7.500000001890827 +0 +0 +0.5000000000209681 +-1.0000000003250102 +-4.999999998921958 +6.333333331585395 +-0.3333333321282528 +-0.5000000001260552 +0 From 8b8628f7d2e0c3ba6310c9d1cef0ab6c2a5f69ca Mon Sep 17 00:00:00 2001 From: Eslick Date: Thu, 11 Apr 2024 17:18:16 -0400 Subject: [PATCH 05/13] Straighten up --- .../cspline_external/cspline_parameters.py | 101 +++++------------- .../cspline_external/src/functions.cpp | 82 +++++++------- ...{function.py => test_external_function.py} | 15 ++- .../test/test_parameter_calculation.py | 51 +++++++++ 4 files changed, 131 insertions(+), 118 deletions(-) rename pyomo/contrib/cspline_external/test/{function.py => test_external_function.py} (61%) create mode 100644 pyomo/contrib/cspline_external/test/test_parameter_calculation.py diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index 6faaf75c820..7bdf6f6ac74 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -1,11 +1,20 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + from pyomo.common.dependencies import attempt_import import pyomo.environ as pyo np, numpy_available = attempt_import("numpy") plt, plt_available = attempt_import("matplotlib.pyplot") -import idaes - def f(x, alpha, s=None): """ @@ -86,22 +95,22 @@ def cubic_parameters_model( ): """Create a Pyomo model to calculate parameters for a cubic spline. By default this creates a square linear model, but optionally it can leave off the endpoint - second derivative constraints and add an objective function for fitting data - instead. The purpose of alternative objective form is to allow the spline to - be constrained in other way and/or extend the spline for a particular - extrapolation behavior. + second derivative constraints and add an objective function for fitting data + instead. The purpose of alternative least squares form is to allow the spline to + be constrained in other that don't require a perfect data match. The knots don't + need to be the same as the x data to allow, for example, additional segments for + extrapolation. Args: x_data: list of x data y_data: list of y data x_knots: optional list of knots - end_point_constraint: if true add constraint that second derivative = 0 at + end_point_constraint: if true add constraint that second derivative = 0 at endpoints objective_form: if true write a least squares objective rather than constraints to match data name: optional model name - Returns: Poymo ConcreteModel """ @@ -147,8 +156,6 @@ def yxx_eqn(blk, s): return pyo.Constraint.Skip return fxx(m.x[s + 1], m.alpha, s) == fxx(m.x[s + 1], m.alpha, s + 1) - # Objective function - # Identify segments used to predict y_data at each x_data. We use search in # instead of a dict lookup, since we don't want to require the data to be at # the knots, even though that is almost always the case. @@ -192,7 +199,9 @@ def yxx_endpoint_eqn(blk, s): def get_parameters(m, file_name=None): - + """Once the model has been solved, this function can be used to extract + the cubic spline parameters. + """ knots = [pyo.value(x) for x in m.x.values()] a1 = [None] * len(m.seg_idx) a2 = [None] * len(m.seg_idx) @@ -213,6 +222,7 @@ def get_parameters(m, file_name=None): return knots, a1, a2, a3, a4 + def _extract_params(m): """Extract alpha as a plain dict of floats to play nice with vectorized functions""" alpha = {} @@ -224,9 +234,10 @@ def _extract_params(m): alpha[s][4] = pyo.value(m.alpha[s, 4]) return alpha + def plot_f(m, file_name=None, **kwargs): """Plot the cspline function. - + Args: m: Pyomo model with data and parameters file_name: optional file to save plot to @@ -252,9 +263,10 @@ def plot_f(m, file_name=None, **kwargs): plt.savefig(file_name, **kwargs) return plt + def plot_fx(m, file_name=None, **kwargs): """Plot the cspline derivative function. - + Args: m: Pyomo model with data and parameters file_name: optional file to save plot to @@ -277,7 +289,7 @@ def plot_fx(m, file_name=None, **kwargs): def plot_fxx(m, file_name=None, **kwargs): """Plot the cspline second derivative function. - + Args: m: Pyomo model with data and parameters file_name: optional file to save plot to @@ -295,64 +307,3 @@ def plot_fxx(m, file_name=None, **kwargs): if file_name is not None: plt.savefig(file_name, **kwargs) return plt - - -if __name__ == "__main__": - - x_knots = [ - 0, - 1, - 2, - 2.5, - 3, - 3.5, - 4, - 5, - 6, - ] - - x_data = [ - 1, - 2, - 3, - 4, - 5, - ] - - y_data = [ - 2, - 3, - 5, - 2, - 1, - ] - - m = cubic_parameters_model( - x_data, - y_data, - x_knots, - end_point_constraint=False, - objective_form=True, - ) - - #m.alpha[1, 1].fix(0) - m.alpha[1, 2].fix(0.5) - m.alpha[1, 3].fix(0) - m.alpha[1, 4].fix(0) - - #m.alpha[8, 1].fix(0) - m.alpha[8, 2].fix(-0.5) - m.alpha[8, 3].fix(0) - m.alpha[8, 4].fix(0) - - - #solver_obj = pyo.SolverFactory("clp") - solver_obj = pyo.SolverFactory("ipopt") - solver_obj.solve(m, tee=True) - - get_parameters(m, "test.txt") - - p = plot_f(m, file_name="f_plot.pdf") - p = plot_fx(m, file_name="fx_plot.pdf") - p = plot_fxx(m, file_name="fxx_plot.pdf") - diff --git a/pyomo/contrib/cspline_external/src/functions.cpp b/pyomo/contrib/cspline_external/src/functions.cpp index 51763b83ad0..78102548a0e 100644 --- a/pyomo/contrib/cspline_external/src/functions.cpp +++ b/pyomo/contrib/cspline_external/src/functions.cpp @@ -20,22 +20,22 @@ // each time, so keep params as long as library remains loaded // Number of loaded curves -unsigned int n_cspline_1d = 0; +unsigned int n_cspline = 0; // Parameter file path to unsigned in index map -std::unordered_map idx_cspline_1d; +std::unordered_map idx_cspline; // Knots for each curve -std::vector< std::vector > knots_cspline_1d; +std::vector< std::vector > knots_cspline; // Cubic segment parameters for each curve -std::vector< std::vector > a1_cspline_1d; -std::vector< std::vector > a2_cspline_1d; -std::vector< std::vector > a3_cspline_1d; -std::vector< std::vector > a4_cspline_1d; +std::vector< std::vector > a1_cspline; +std::vector< std::vector > a2_cspline; +std::vector< std::vector > a3_cspline; +std::vector< std::vector > a4_cspline; // Number of segments in each curve -std::vector< unsigned int > n_seg_cspline_1d; +std::vector< unsigned int > n_seg_cspline; // The data file format is a text file with one parameter per line // Line 1: Number of parameters @@ -46,60 +46,60 @@ std::vector< unsigned int > n_seg_cspline_1d; // nsegs lines for a4 // Function to read parameter file if not read and store params -static int read_parameters_cspline_1d(std::string file_path) { +static int read_parameters_cspline(std::string file_path) { try{ //if read already return curve index and done. - return idx_cspline_1d.at(file_path); + return idx_cspline.at(file_path); } catch(std::out_of_range const&){} // Not read so read it. unsigned int i, idx, n; //loop counter, curve index, and number of segs std::ifstream param_file; // parameter input file // Set index for cspline param file and increment curve count - idx_cspline_1d[file_path] = n_cspline_1d; - idx = n_cspline_1d; - ++n_cspline_1d; + idx_cspline[file_path] = n_cspline; + idx = n_cspline; + ++n_cspline; // open the parameter file for input param_file.open(file_path); // get the number of segments and size to vectors param_file >> n; - n_seg_cspline_1d.resize(n_cspline_1d); - n_seg_cspline_1d[idx] = n; - knots_cspline_1d.resize(n_cspline_1d); - a1_cspline_1d.resize(n_cspline_1d); - a2_cspline_1d.resize(n_cspline_1d); - a3_cspline_1d.resize(n_cspline_1d); - a4_cspline_1d.resize(n_cspline_1d); - knots_cspline_1d[idx].resize(n + 1); - a1_cspline_1d[idx].resize(n); - a2_cspline_1d[idx].resize(n); - a3_cspline_1d[idx].resize(n); - a4_cspline_1d[idx].resize(n); + n_seg_cspline.resize(n_cspline); + n_seg_cspline[idx] = n; + knots_cspline.resize(n_cspline); + a1_cspline.resize(n_cspline); + a2_cspline.resize(n_cspline); + a3_cspline.resize(n_cspline); + a4_cspline.resize(n_cspline); + knots_cspline[idx].resize(n + 1); + a1_cspline[idx].resize(n); + a2_cspline[idx].resize(n); + a3_cspline[idx].resize(n); + a4_cspline[idx].resize(n); // get the knots for(i=0; i < n + 1; ++i){ - param_file >> knots_cspline_1d[idx][i]; + param_file >> knots_cspline[idx][i]; } // get the a1 params for(i=0; i < n; ++i){ - param_file >> a1_cspline_1d[idx][i]; + param_file >> a1_cspline[idx][i]; } // get the a2 params for(i=0; i < n; ++i){ - param_file >> a2_cspline_1d[idx][i]; + param_file >> a2_cspline[idx][i]; } // get the a3 params for(i=0; i < n; ++i){ - param_file >> a3_cspline_1d[idx][i]; + param_file >> a3_cspline[idx][i]; } // get the a4 params for(i=0; i < n; ++i){ - param_file >> a4_cspline_1d[idx][i]; + param_file >> a4_cspline[idx][i]; } param_file.close(); @@ -109,19 +109,19 @@ static int read_parameters_cspline_1d(std::string file_path) { } // Calculate the cspline 1D function value and derivatives -extern real cspline_1d(arglist *al) { +extern real cspline(arglist *al) { // Get the data file path argument const char* data_file = al->sa[-(al->at[1]+1)]; // Get the function input real x = al->ra[al->at[0]]; // Get parameters, read file if needed - unsigned int idx = read_parameters_cspline_1d(data_file); + unsigned int idx = read_parameters_cspline(data_file); size_t seg; // segment index real a1, a2, a3, a4; // segment parameters //find segment index - auto lit = std::lower_bound(knots_cspline_1d[idx].begin(), knots_cspline_1d[idx].end(), x); - seg = lit - knots_cspline_1d[idx].begin(); + auto lit = std::lower_bound(knots_cspline[idx].begin(), knots_cspline[idx].end(), x); + seg = lit - knots_cspline[idx].begin(); // If off the curve on the low side use first seg (0) otherwise // subtract 1 to go from knot index to segment index (i.e. seg is @@ -130,13 +130,13 @@ extern real cspline_1d(arglist *al) { // segment index) if(seg > 0) --seg; // In off the curve on the high side, use last segment - if(seg >= n_seg_cspline_1d[idx]) seg = n_seg_cspline_1d[idx] - 1; + if(seg >= n_seg_cspline[idx]) seg = n_seg_cspline[idx] - 1; // Get the parameters for the correct segment - a1 = a1_cspline_1d[idx][seg]; - a2 = a2_cspline_1d[idx][seg]; - a3 = a3_cspline_1d[idx][seg]; - a4 = a4_cspline_1d[idx][seg]; + a1 = a1_cspline[idx][seg]; + a2 = a2_cspline[idx][seg]; + a3 = a3_cspline[idx][seg]; + a4 = a4_cspline[idx][seg]; // Compute the first derivative, if requested. if (al->derivs!=NULL) { @@ -162,8 +162,8 @@ void funcadd(AmplExports *ae){ // arguments) // 5) Void pointer to function info addfunc( - "cspline_1d", - (rfunc)cspline_1d, + "cspline", + (rfunc)cspline, FUNCADD_REAL_VALUED|FUNCADD_STRING_ARGS, 2, NULL diff --git a/pyomo/contrib/cspline_external/test/function.py b/pyomo/contrib/cspline_external/test/test_external_function.py similarity index 61% rename from pyomo/contrib/cspline_external/test/function.py rename to pyomo/contrib/cspline_external/test/test_external_function.py index 2c974f80e54..3bf60f1eaff 100644 --- a/pyomo/contrib/cspline_external/test/function.py +++ b/pyomo/contrib/cspline_external/test/test_external_function.py @@ -1,12 +1,23 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + import pyomo.environ as pyo from pyomo.common.fileutils import find_library if __name__ == "__main__": lib = find_library("cspline_external") - params = "t2_params.txt" + params = "t1_params.txt" m = pyo.ConcreteModel() - m.f = pyo.ExternalFunction(library=lib, function="cspline_1d") + m.f = pyo.ExternalFunction(library=lib, function="cspline") m.x = pyo.Var(initialize=2, bounds=(0.5, 5.5)) m.y = pyo.Var() diff --git a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py new file mode 100644 index 00000000000..b18005d7af2 --- /dev/null +++ b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py @@ -0,0 +1,51 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +from pyomo.contrib.cspline_external.cspline_parameters import ( + cubic_parameters_model, + get_parameters, +) +import pyomo.environ as pyo + + +if __name__ == "__main__": + + x_data = [ + 1, + 2, + 3, + 4, + 5, + ] + + y_data = [ + 2, + 3, + 5, + 2, + 1, + ] + + m = cubic_parameters_model( + x_data, + y_data, + end_point_constraint=True, + objective_form=False, + ) + + # REMOVE, temporary just to set solver path + import idaes + + # solver_obj = pyo.SolverFactory("clp") + solver_obj = pyo.SolverFactory("ipopt") + solver_obj.solve(m, tee=True) + + get_parameters(m, "test.txt") From 5116d1e0ec6f7097c13c605d8c03a270b7fe9e68 Mon Sep 17 00:00:00 2001 From: Eslick Date: Thu, 11 Apr 2024 17:22:46 -0400 Subject: [PATCH 06/13] Run black --- .../test/test_parameter_calculation.py | 23 ++++--------------- 1 file changed, 4 insertions(+), 19 deletions(-) diff --git a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py index b18005d7af2..ee8ffc28308 100644 --- a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py +++ b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py @@ -18,27 +18,12 @@ if __name__ == "__main__": - x_data = [ - 1, - 2, - 3, - 4, - 5, - ] - - y_data = [ - 2, - 3, - 5, - 2, - 1, - ] + x_data = [1, 2, 3, 4, 5] + + y_data = [2, 3, 5, 2, 1] m = cubic_parameters_model( - x_data, - y_data, - end_point_constraint=True, - objective_form=False, + x_data, y_data, end_point_constraint=True, objective_form=False ) # REMOVE, temporary just to set solver path From 87ee33c5a490e2171622d25aa67f98c0312388a8 Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 12:45:19 -0400 Subject: [PATCH 07/13] A little clean up, erase plotting funcs --- .../cspline_external/cspline_parameters.py | 110 +++--------------- 1 file changed, 18 insertions(+), 92 deletions(-) diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index 7bdf6f6ac74..4d7e45edf22 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -16,7 +16,8 @@ plt, plt_available = attempt_import("matplotlib.pyplot") -def f(x, alpha, s=None): + +def _f_cubic(x, alpha, s=None): """ Cubic function: y = a1 + a2 * x + a3 * x^2 + a4 * x^3 @@ -39,7 +40,7 @@ def f(x, alpha, s=None): return alpha[s, 1] + alpha[s, 2] * x + alpha[s, 3] * x**2 + alpha[s, 4] * x**3 -def fx(x, alpha, s=None): +def _fx_cubic(x, alpha, s=None): """ Cubic function first derivative: dy/dx = a2 + 2 * a3 * x + 3 * a4 * x^2 @@ -62,7 +63,7 @@ def fx(x, alpha, s=None): return alpha[s, 2] + 2 * alpha[s, 3] * x + 3 * alpha[s, 4] * x**2 -def fxx(x, alpha, s=None): +def _fxx_cubic(x, alpha, s=None): """ Cubic function second derivative: d2y/dx2 = 2 * a3 + 6 * a4 * x @@ -97,18 +98,19 @@ def cubic_parameters_model( this creates a square linear model, but optionally it can leave off the endpoint second derivative constraints and add an objective function for fitting data instead. The purpose of alternative least squares form is to allow the spline to - be constrained in other that don't require a perfect data match. The knots don't - need to be the same as the x data to allow, for example, additional segments for - extrapolation. + be constrained in other wats that don't require a perfect data match. The knots + don't need to be the same as the x data to allow, for example, additional segments + for extrapolation. This is not the most computationally efficient way to calculate + parameters, but since it is used to precalculate parameters, speed is not important. Args: x_data: list of x data y_data: list of y data - x_knots: optional list of knots - end_point_constraint: if true add constraint that second derivative = 0 at - endpoints - objective_form: if true write a least squares objective rather than constraints - to match data + x_knots: optional list of knots (default is to use x_data) + end_point_constraint: if True add constraint that second derivative = 0 at + endpoints (default=True) + objective_form: if True write a least squares objective rather than constraints + to match data (default=False) name: optional model name Returns: @@ -118,11 +120,9 @@ def cubic_parameters_model( assert n_data == len(y_data) if x_knots is None: n_knots = n_data - n_seg = n_data - 1 x_knots = x_data else: n_knots = len(x_knots) - n_seg = n_knots - 1 m = pyo.ConcreteModel(name=name) # Sets of indexes for knots, segments, and data @@ -140,21 +140,21 @@ def cubic_parameters_model( def y_eqn(blk, s): if s == m.seg_idx.last(): return pyo.Constraint.Skip - return f(m.x[s + 1], m.alpha, s) == f(m.x[s + 1], m.alpha, s + 1) + return _f_cubic(m.x[s + 1], m.alpha, s) == _f_cubic(m.x[s + 1], m.alpha, s + 1) # f'_s(x) = f'_s+1(x) @m.Constraint(m.seg_idx) def yx_eqn(blk, s): if s == m.seg_idx.last(): return pyo.Constraint.Skip - return fx(m.x[s + 1], m.alpha, s) == fx(m.x[s + 1], m.alpha, s + 1) + return _fx_cubic(m.x[s + 1], m.alpha, s) == _fx_cubic(m.x[s + 1], m.alpha, s + 1) # f"_s(x) = f"_s+1(x) @m.Constraint(m.seg_idx) def yxx_eqn(blk, s): if s == m.seg_idx.last(): return pyo.Constraint.Skip - return fxx(m.x[s + 1], m.alpha, s) == fxx(m.x[s + 1], m.alpha, s + 1) + return _fxx_cubic(m.x[s + 1], m.alpha, s) == _fxx_cubic(m.x[s + 1], m.alpha, s + 1) # Identify segments used to predict y_data at each x_data. We use search in # instead of a dict lookup, since we don't want to require the data to be at @@ -170,7 +170,7 @@ def ydiff(blk, d): s = idx[d - 1] + 1 if s >= m.seg_idx.last(): s -= 1 - return m.y_data[d] - f(m.x_data[d], m.alpha, s) + return m.y_data[d] - _f_cubic(m.x_data[d], m.alpha, s) if objective_form: # least squares objective @@ -195,7 +195,7 @@ def yxx_endpoint_eqn(blk, s): j = m.knt_idx.last() else: j = s - return fxx(m.x[j], m.alpha, s) == 0 + return _fxx_cubic(m.x[j], m.alpha, s) == 0 def get_parameters(m, file_name=None): @@ -233,77 +233,3 @@ def _extract_params(m): alpha[s][3] = pyo.value(m.alpha[s, 3]) alpha[s][4] = pyo.value(m.alpha[s, 4]) return alpha - - -def plot_f(m, file_name=None, **kwargs): - """Plot the cspline function. - - Args: - m: Pyomo model with data and parameters - file_name: optional file to save plot to - - Returns: - pyplot object - """ - if not plt_available: - raise ModuleNotFoundError("Matplotlib is not available") - plt.close() - alpha = _extract_params(m) - for s in m.seg_idx: - xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) - plt.plot(xvec, f(xvec, alpha[s])) - plt.title("f(x)") - x = [] - y = [] - for i in m.dat_idx: - x.append(pyo.value(m.x_data[i])) - y.append(pyo.value(m.y_data[i])) - plt.scatter(x, y) - if file_name is not None: - plt.savefig(file_name, **kwargs) - return plt - - -def plot_fx(m, file_name=None, **kwargs): - """Plot the cspline derivative function. - - Args: - m: Pyomo model with data and parameters - file_name: optional file to save plot to - - Returns: - pyplot object - """ - if not plt_available: - raise ModuleNotFoundError("Matplotlib is not available") - plt.close() - alpha = _extract_params(m) - for s in m.seg_idx: - xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) - plt.plot(xvec, fx(xvec, alpha[s])) - plt.title("f'(x)") - if file_name is not None: - plt.savefig(file_name, **kwargs) - return plt - - -def plot_fxx(m, file_name=None, **kwargs): - """Plot the cspline second derivative function. - - Args: - m: Pyomo model with data and parameters - file_name: optional file to save plot to - - Returns: - pyplot object - """ - if not plt_available: - raise ModuleNotFoundError("Matplotlib is not available") - plt.close() - alpha = _extract_params(m) - for s in m.seg_idx: - xvec = np.linspace(pyo.value(m.x[s]), pyo.value(m.x[s + 1]), 20) - plt.plot(xvec, fxx(xvec, alpha[s])) - if file_name is not None: - plt.savefig(file_name, **kwargs) - return plt From 385a689138c799e8d8b6914f980eb19c38eb9c7d Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 13:08:04 -0400 Subject: [PATCH 08/13] Finish clean up --- .../cspline_external/cspline_parameters.py | 13 ++++++----- .../cspline_external/test/t1_params.txt | 22 ------------------- 2 files changed, 8 insertions(+), 27 deletions(-) delete mode 100644 pyomo/contrib/cspline_external/test/t1_params.txt diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index 4d7e45edf22..c163d156af3 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -13,9 +13,12 @@ import pyomo.environ as pyo np, numpy_available = attempt_import("numpy") -plt, plt_available = attempt_import("matplotlib.pyplot") - +if numpy_available: + searchsorted = np.searchsorted +else: + def searchsorted(x, y): + raise NotImplementedError("cubic_parameters_model() currently relies on Numpy") def _f_cubic(x, alpha, s=None): """ @@ -104,9 +107,9 @@ def cubic_parameters_model( parameters, but since it is used to precalculate parameters, speed is not important. Args: - x_data: list of x data + x_data: sorted list of x data y_data: list of y data - x_knots: optional list of knots (default is to use x_data) + x_knots: optional sorted list of knots (default is to use x_data) end_point_constraint: if True add constraint that second derivative = 0 at endpoints (default=True) objective_form: if True write a least squares objective rather than constraints @@ -159,7 +162,7 @@ def yxx_eqn(blk, s): # Identify segments used to predict y_data at each x_data. We use search in # instead of a dict lookup, since we don't want to require the data to be at # the knots, even though that is almost always the case. - idx = np.searchsorted(x_knots, x_data) + idx = searchsorted(x_knots, x_data) if end_point_constraint: add_endpoint_second_derivative_constraints(m) diff --git a/pyomo/contrib/cspline_external/test/t1_params.txt b/pyomo/contrib/cspline_external/test/t1_params.txt deleted file mode 100644 index bc6913c6b57..00000000000 --- a/pyomo/contrib/cspline_external/test/t1_params.txt +++ /dev/null @@ -1,22 +0,0 @@ -4 -1 -2 -3 -4 -5 -0.9999999999999978 -24.714285714286376 --106.42857142857616 -115.28571428572297 -2.321428571428605 --33.250000000000966 -97.89285714286154 --68.39285714286291 --1.9821428571429038 -15.803571428571882 --27.91071428571562 -13.660714285715482 -0.6607142857143012 --2.3035714285714963 -2.55357142857156 --0.9107142857143654 From 64c4a007df1a3d6f72a257fadf99cbcd75de5b13 Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 13:09:23 -0400 Subject: [PATCH 09/13] Move test params --- .../cspline_external/test/t2_params.txt | 42 ------------------- .../cspline_external/test/test_params.txt | 22 ++++++++++ 2 files changed, 22 insertions(+), 42 deletions(-) delete mode 100644 pyomo/contrib/cspline_external/test/t2_params.txt create mode 100644 pyomo/contrib/cspline_external/test/test_params.txt diff --git a/pyomo/contrib/cspline_external/test/t2_params.txt b/pyomo/contrib/cspline_external/test/t2_params.txt deleted file mode 100644 index 7d7609a3987..00000000000 --- a/pyomo/contrib/cspline_external/test/t2_params.txt +++ /dev/null @@ -1,42 +0,0 @@ -8 -0 -1 -2 -2.5 -3 -3.5 -4 -5 -6 -1.4999999999908438 -0.9999999999698757 -13.000000002737705 -75.49999998081503 --230.49999994288345 -55.33333326383919 -66.00000001569855 -3.4999999999416516 -0.5 -2.000000000062904 --16.00000000408884 --90.99999997778163 -214.99999994591695 --29.99999994555961 --38.00000000945413 --0.5 -0 --1.500000000062904 -7.500000002012967 -37.49999999149008 --64.49999998307611 -5.499999985917197 -7.500000001890827 -0 -0 -0.5000000000209681 --1.0000000003250102 --4.999999998921958 -6.333333331585395 --0.3333333321282528 --0.5000000001260552 -0 diff --git a/pyomo/contrib/cspline_external/test/test_params.txt b/pyomo/contrib/cspline_external/test/test_params.txt new file mode 100644 index 00000000000..36620929582 --- /dev/null +++ b/pyomo/contrib/cspline_external/test/test_params.txt @@ -0,0 +1,22 @@ +4 +1 +2 +3 +4 +5 +1.000000000000001 +24.71428571428584 +-106.4285714285722 +115.28571428571576 +2.3214285714285774 +-33.25000000000018 +97.89285714285785 +-68.39285714285812 +-1.9821428571428679 +15.803571428571509 +-27.9107142857145 +13.660714285714493 +0.6607142857142894 +-2.3035714285714404 +2.55357142857145 +-0.9107142857142996 From bb6f4b28d8d090406be0d2900c212c10abd39fac Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 16:26:32 -0400 Subject: [PATCH 10/13] Add tests for external function --- pyomo/contrib/cspline_external/build.py | 8 +- pyomo/contrib/cspline_external/plugins.py | 6 +- .../cspline_external/src/CMakeLists.txt | 2 +- .../test/test_external_function.py | 137 ++++++++++++++---- .../test/test_params_line.txt | 7 + pyomo/environ/__init__.py | 1 + 6 files changed, 124 insertions(+), 37 deletions(-) create mode 100644 pyomo/contrib/cspline_external/test/test_params_line.txt diff --git a/pyomo/contrib/cspline_external/build.py b/pyomo/contrib/cspline_external/build.py index 8ca4119dbe1..eeccff1dbed 100644 --- a/pyomo/contrib/cspline_external/build.py +++ b/pyomo/contrib/cspline_external/build.py @@ -13,7 +13,7 @@ from pyomo.common.cmake_builder import build_cmake_project -def build_cspline_1d_external(user_args=[], parallel=None): +def build_cspline_external(user_args=[], parallel=None): return build_cmake_project( targets=["src"], package_name="cspline_1d_external", @@ -23,10 +23,10 @@ def build_cspline_1d_external(user_args=[], parallel=None): ) -class AMPLCspline1DExternalBuilder(object): +class AMPLCsplineExternalBuilder(object): def __call__(self, parallel): - return build_cspline_1d_external(parallel=parallel) + return build_cspline_external(parallel=parallel) if __name__ == "__main__": - build_cspline_1d_external(sys.argv[1:]) + build_cspline_external(sys.argv[1:]) diff --git a/pyomo/contrib/cspline_external/plugins.py b/pyomo/contrib/cspline_external/plugins.py index 742fd6afdd4..0231b4b4ca0 100644 --- a/pyomo/contrib/cspline_external/plugins.py +++ b/pyomo/contrib/cspline_external/plugins.py @@ -10,10 +10,10 @@ # ___________________________________________________________________________ from pyomo.common.extensions import ExtensionBuilderFactory -from pyomo.contrib.cspline_1d_external.build import AMPLCspline1DExternalBuilder +from pyomo.contrib.cspline_external.build import AMPLCsplineExternalBuilder def load(): - ExtensionBuilderFactory.register("cspline_1d_external")( - AMPLCspline1DExternalBuilder + ExtensionBuilderFactory.register("cspline_external")( + AMPLCsplineExternalBuilder ) diff --git a/pyomo/contrib/cspline_external/src/CMakeLists.txt b/pyomo/contrib/cspline_external/src/CMakeLists.txt index b710bc948b6..f3a7f9a931a 100644 --- a/pyomo/contrib/cspline_external/src/CMakeLists.txt +++ b/pyomo/contrib/cspline_external/src/CMakeLists.txt @@ -12,7 +12,7 @@ cmake_minimum_required(VERSION 3.0) # CMake 3.0 needed by AMPL/asl build -PROJECT( ampl_function_demo ) +PROJECT( cspline_external ) INCLUDE( "${CMAKE_CURRENT_SOURCE_DIR}/../../cspline_external/src/FindASL.cmake") diff --git a/pyomo/contrib/cspline_external/test/test_external_function.py b/pyomo/contrib/cspline_external/test/test_external_function.py index 3bf60f1eaff..6b1d220c284 100644 --- a/pyomo/contrib/cspline_external/test/test_external_function.py +++ b/pyomo/contrib/cspline_external/test/test_external_function.py @@ -8,42 +8,121 @@ # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ - +import os +import pyomo.common.unittest as unittest import pyomo.environ as pyo -from pyomo.common.fileutils import find_library +from pyomo.common.fileutils import find_library, this_file_dir + +_lib = find_library("cspline_external") +is_pypy = platform.python_implementation().lower().startswith('pypy') -if __name__ == "__main__": - lib = find_library("cspline_external") - params = "t1_params.txt" +@unittest.skipIf(is_pypy, 'Cannot evaluate external functions under pypy') +@unittest.skipIf(not _lib, "cspline library is not available.") +class CsplineExternal1DTest(unittest.TestCase): - m = pyo.ConcreteModel() - m.f = pyo.ExternalFunction(library=lib, function="cspline") + def test_function_call(self): + """Test that the cspline behaves as expected""" + # The parameters for the function call are stored in + # this file. The file is only read once by the external + # function on the first call. + params = os.path.join(this_file_dir(), "test_params.txt") + # this is that data that was used to generate the params + x_data = [1, 2, 3, 4, 5] + y_data = [2, 3, 5, 2, 1] + # Model with the cspline function + m = pyo.ConcreteModel() + m.f = pyo.ExternalFunction(library=_lib, function="cspline") + m.x = pyo.Var(initialize=2) + m.e = pyo.Expression(expr=m.f(m.x, params)) + # We know that spline should go through the data points + # and that f(x) f'(x) and f''(x) for the right and left + # segment at a knot should be the same. + for x, y in zip(x_data, y_data): + delta = 1e-5 + f, fx, fxx = m.f.evaluate_fgh(args=(x, params)) + f_left, fx_left, fxx_left = m.f.evaluate_fgh(args=(x - delta, params)) + f_right, fx_right, fxx_right = m.f.evaluate_fgh(args=(x + delta, params)) + self.assertAlmostEqual(f, y, 8) + # check left/right + self.assertAlmostEqual(f_left, y, 3) + self.assertAlmostEqual(f_right, y, 3) + # check left/right derivatives + self.assertAlmostEqual(fx_left[0], fx[0], 3) + self.assertAlmostEqual(fx_right[0], fx[0], 3) + self.assertAlmostEqual(fxx_left[0], fxx[0], 3) + self.assertAlmostEqual(fxx_right[0], fxx[0], 3) + # we know the endpoint constraints are that f''(0) = 0 + f, fx, fxx = m.f.evaluate_fgh(args=(x_data[0], params)) + self.assertAlmostEqual(fxx[0], 0, 8) + f, fx, fxx = m.f.evaluate_fgh(args=(x_data[-1], params)) + self.assertAlmostEqual(fxx[0], 0, 8) + # check a little more in model context + self.assertAlmostEqual(pyo.value(m.e), 3, 8) - m.x = pyo.Var(initialize=2, bounds=(0.5, 5.5)) - m.y = pyo.Var() + def test_ampl_derivatives(self): + # Make sure the function values and derivatives are right + # based on parameters. This time look at segment mid-points - m.c1 = pyo.Constraint(expr=m.y == m.f(m.x, params)) + # The parameters for the function call are stored in + # this file. The file is only read once by the external + # function on the first call. + params = os.path.join(this_file_dir(), "test_params.txt") + # this is that data that was used to generate the params + x_data = [1, 2, 3, 4, 5] + y_data = [2, 3, 5, 2, 1] + # Model with the cspline function + m = pyo.ConcreteModel() + m.f = pyo.ExternalFunction(library=_lib, function="cspline") + # Read parameters + with open(params, "r") as fptr: + # line 1: number of knots + ns = int(fptr.readline()) + # Make param lists + knots = [None] * (ns + 1) + a = [None] * 4 + for j in range(4): + a[j] = [None] * ns - m.o = pyo.Objective(expr=-m.y) + # Read params + for i in range(ns + 1): + knots[i] = float(fptr.readline()) + for j in range(4): + for i in range(ns): + a[j][i] = float(fptr.readline()) + # Check the value calculated by from parameters to external + # function for each segment. Use the mid. + for i in range(ns): + x = (x_data[i] + x_data[i + 1]) / 2.0 + y = a[0][i] + a[1][i] * x + a[2][i] * x**2 + a[3][i] * x**3 + yx = a[1][i] + 2 * a[2][i] * x + 3 * a[3][i] * x**2 + yxx = 2 * a[2][i] + 6 * a[3][i] * x + f, fx, fxx = m.f.evaluate_fgh(args=(x, params)) + self.assertAlmostEqual(f, y, 8) + self.assertAlmostEqual(fx[0], yx, 8) + self.assertAlmostEqual(fxx[0], yxx, 8) - print(pyo.value(m.f(0.9, params))) - print(pyo.value(m.f(1.01, params))) - print(pyo.value(m.f(1.05, params))) - print(pyo.value(m.f(1.2, params))) - print(pyo.value(m.f(1.3, params))) - print(pyo.value(m.f(1.5, params))) - print(pyo.value(m.f(2, params))) - print(pyo.value(m.f(2.01, params))) - print(pyo.value(m.f(2.1, params))) - print(pyo.value(m.f(3.0, params))) - print(pyo.value(m.f(3.1, params))) - print(pyo.value(m.f(4.9, params))) - print(pyo.value(m.f(5.0, params))) - print(pyo.value(m.f(5.1, params))) + def test_load_multiple_splines(self): + # The last test ensures that you can load multiple splines + # for a model without trouble. You should be able to use + # as many parameter set as you want. The external function + # only reads a file once as long as the library remains loaded. - import idaes + # first spline is the one we used before + params1 = os.path.join(this_file_dir(), "test_params.txt") + # second spline is a line with intercept = 0, slope = 1 + params2 = os.path.join(this_file_dir(), "test_params_line.txt") + # Model with the cspline function + m = pyo.ConcreteModel() + m.f = pyo.ExternalFunction(library=_lib, function="cspline") + m.x = pyo.Var(initialize=1) + m.e1 = pyo.Expression(expr=m.f(m.x, params1)) + m.e2 = pyo.Expression(expr=m.f(m.x, params2)) - solver_obj = pyo.SolverFactory("ipopt") - solver_obj.solve(m, tee=True) + # Make sure both functions return correct value + self.assertAlmostEqual(pyo.value(m.e1), 2, 8) + self.assertAlmostEqual(pyo.value(m.e2), 1, 8) - m.display() + # make sure loading the second set of parameters + # didn't mess up the first + self.assertAlmostEqual(pyo.value(m.e1), 2, 8) + self.assertAlmostEqual(pyo.value(m.e2), 1, 8) diff --git a/pyomo/contrib/cspline_external/test/test_params_line.txt b/pyomo/contrib/cspline_external/test/test_params_line.txt new file mode 100644 index 00000000000..3015890f0c6 --- /dev/null +++ b/pyomo/contrib/cspline_external/test/test_params_line.txt @@ -0,0 +1,7 @@ +1 +0 +1 +0 +1 +0 +0 \ No newline at end of file diff --git a/pyomo/environ/__init__.py b/pyomo/environ/__init__.py index 07b3dfad680..1e426bbb966 100644 --- a/pyomo/environ/__init__.py +++ b/pyomo/environ/__init__.py @@ -39,6 +39,7 @@ def _do_import(pkg_name): 'pyomo.contrib.appsi', 'pyomo.contrib.community_detection', 'pyomo.contrib.cp', + 'pyomo.contrib.cspline_external', 'pyomo.contrib.example', 'pyomo.contrib.fme', 'pyomo.contrib.gdp_bounds', From 002e68a6f479e3ae0304ecb66e749dc0f5345529 Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 17:04:00 -0400 Subject: [PATCH 11/13] Add test for calculate parameters --- .../cspline_external/cspline_parameters.py | 14 +++-- pyomo/contrib/cspline_external/plugins.py | 4 +- .../test/test_external_function.py | 6 +- .../test/test_parameter_calculation.py | 62 +++++++++++++++---- 4 files changed, 65 insertions(+), 21 deletions(-) diff --git a/pyomo/contrib/cspline_external/cspline_parameters.py b/pyomo/contrib/cspline_external/cspline_parameters.py index c163d156af3..64ef92cfdef 100644 --- a/pyomo/contrib/cspline_external/cspline_parameters.py +++ b/pyomo/contrib/cspline_external/cspline_parameters.py @@ -17,9 +17,11 @@ if numpy_available: searchsorted = np.searchsorted else: + def searchsorted(x, y): raise NotImplementedError("cubic_parameters_model() currently relies on Numpy") + def _f_cubic(x, alpha, s=None): """ Cubic function: @@ -101,10 +103,10 @@ def cubic_parameters_model( this creates a square linear model, but optionally it can leave off the endpoint second derivative constraints and add an objective function for fitting data instead. The purpose of alternative least squares form is to allow the spline to - be constrained in other wats that don't require a perfect data match. The knots + be constrained in other wats that don't require a perfect data match. The knots don't need to be the same as the x data to allow, for example, additional segments for extrapolation. This is not the most computationally efficient way to calculate - parameters, but since it is used to precalculate parameters, speed is not important. + parameters, but since it is used to precalculate parameters, speed is not important. Args: x_data: sorted list of x data @@ -150,14 +152,18 @@ def y_eqn(blk, s): def yx_eqn(blk, s): if s == m.seg_idx.last(): return pyo.Constraint.Skip - return _fx_cubic(m.x[s + 1], m.alpha, s) == _fx_cubic(m.x[s + 1], m.alpha, s + 1) + return _fx_cubic(m.x[s + 1], m.alpha, s) == _fx_cubic( + m.x[s + 1], m.alpha, s + 1 + ) # f"_s(x) = f"_s+1(x) @m.Constraint(m.seg_idx) def yxx_eqn(blk, s): if s == m.seg_idx.last(): return pyo.Constraint.Skip - return _fxx_cubic(m.x[s + 1], m.alpha, s) == _fxx_cubic(m.x[s + 1], m.alpha, s + 1) + return _fxx_cubic(m.x[s + 1], m.alpha, s) == _fxx_cubic( + m.x[s + 1], m.alpha, s + 1 + ) # Identify segments used to predict y_data at each x_data. We use search in # instead of a dict lookup, since we don't want to require the data to be at diff --git a/pyomo/contrib/cspline_external/plugins.py b/pyomo/contrib/cspline_external/plugins.py index 0231b4b4ca0..d37c4f92cd9 100644 --- a/pyomo/contrib/cspline_external/plugins.py +++ b/pyomo/contrib/cspline_external/plugins.py @@ -14,6 +14,4 @@ def load(): - ExtensionBuilderFactory.register("cspline_external")( - AMPLCsplineExternalBuilder - ) + ExtensionBuilderFactory.register("cspline_external")(AMPLCsplineExternalBuilder) diff --git a/pyomo/contrib/cspline_external/test/test_external_function.py b/pyomo/contrib/cspline_external/test/test_external_function.py index 6b1d220c284..c48b219decd 100644 --- a/pyomo/contrib/cspline_external/test/test_external_function.py +++ b/pyomo/contrib/cspline_external/test/test_external_function.py @@ -9,14 +9,16 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ import os +import platform import pyomo.common.unittest as unittest import pyomo.environ as pyo from pyomo.common.fileutils import find_library, this_file_dir _lib = find_library("cspline_external") -is_pypy = platform.python_implementation().lower().startswith('pypy') +is_pypy = platform.python_implementation().lower().startswith("pypy") -@unittest.skipIf(is_pypy, 'Cannot evaluate external functions under pypy') + +@unittest.skipIf(is_pypy, "Cannot evaluate external functions under pypy") @unittest.skipIf(not _lib, "cspline library is not available.") class CsplineExternal1DTest(unittest.TestCase): diff --git a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py index ee8ffc28308..a572b8909e1 100644 --- a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py +++ b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py @@ -9,28 +9,66 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ +import os + +from pyomo.common.dependencies import attempt_import from pyomo.contrib.cspline_external.cspline_parameters import ( cubic_parameters_model, get_parameters, ) import pyomo.environ as pyo +from pyomo.common.fileutils import this_file_dir +import pyomo.common.unittest as unittest + +np, numpy_available = attempt_import("numpy") + +# REMOVE, temporary just to set solver path +import idaes -if __name__ == "__main__": +@unittest.skipIf(not numpy_available, "numpy is not available.") +class CsplineExternalParamsTest(unittest.TestCase): + def test_param_gen(self): + x_data = [1, 2, 3, 4, 5] + y_data = [2, 3, 5, 2, 1] - x_data = [1, 2, 3, 4, 5] + m = cubic_parameters_model(x_data, y_data) - y_data = [2, 3, 5, 2, 1] + # This is a linear system of equations by default, so although ipopt + # is not needed, it works perfectly well. + solver_obj = pyo.SolverFactory("ipopt") + solver_obj.solve(m) - m = cubic_parameters_model( - x_data, y_data, end_point_constraint=True, objective_form=False - ) + params = os.path.join(this_file_dir(), "test.txt") + knots, a0, a1, a2, a3 = get_parameters(m, params) + assert len(knots) == 5 + assert len(a0) == 4 + assert len(a1) == 4 + assert len(a2) == 4 + assert len(a3) == 4 - # REMOVE, temporary just to set solver path - import idaes + # Read parameters + with open(params, "r") as fptr: + # line 1: number of knots + ns = int(fptr.readline()) + # Make param lists + knots = [None] * (ns + 1) + a = [None] * 4 + for j in range(4): + a[j] = [None] * ns - # solver_obj = pyo.SolverFactory("clp") - solver_obj = pyo.SolverFactory("ipopt") - solver_obj.solve(m, tee=True) + # Read params + for i in range(ns + 1): + knots[i] = float(fptr.readline()) + for j in range(4): + for i in range(ns): + a[j][i] = float(fptr.readline()) + # Check the value calculated by from parameters + for i, x in enumerate(x_data): + seg = i - 1 + if seg < 0: + seg = 0 + y = a[0][seg] + a[1][seg] * x + a[2][seg] * x**2 + a[3][seg] * x**3 + self.assertAlmostEqual(y, y_data[i], 8) - get_parameters(m, "test.txt") + os.remove(params) From 368ea950428541ab9dbe55ad4c89a694654e7297 Mon Sep 17 00:00:00 2001 From: Eslick Date: Tue, 14 May 2024 17:07:17 -0400 Subject: [PATCH 12/13] Handle no ipopt --- .../cspline_external/test/test_parameter_calculation.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py index a572b8909e1..9974c7fe4c0 100644 --- a/pyomo/contrib/cspline_external/test/test_parameter_calculation.py +++ b/pyomo/contrib/cspline_external/test/test_parameter_calculation.py @@ -16,16 +16,17 @@ cubic_parameters_model, get_parameters, ) +from pyomo.opt import check_available_solvers import pyomo.environ as pyo from pyomo.common.fileutils import this_file_dir import pyomo.common.unittest as unittest np, numpy_available = attempt_import("numpy") -# REMOVE, temporary just to set solver path -import idaes - +@unittest.skipUnless( + check_available_solvers("ipopt"), "The 'ipopt' solver is not available" +) @unittest.skipIf(not numpy_available, "numpy is not available.") class CsplineExternalParamsTest(unittest.TestCase): def test_param_gen(self): From fd9cd4e2730e1f1ebc6b08ccc902ad0698692ec4 Mon Sep 17 00:00:00 2001 From: Eslick Date: Wed, 22 May 2024 08:54:15 -0400 Subject: [PATCH 13/13] Add init file --- pyomo/contrib/cspline_external/src/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 pyomo/contrib/cspline_external/src/__init__.py diff --git a/pyomo/contrib/cspline_external/src/__init__.py b/pyomo/contrib/cspline_external/src/__init__.py new file mode 100644 index 00000000000..e69de29bb2d