From ad49549c44f95b6e0bf79c5a9e2260320de7920a Mon Sep 17 00:00:00 2001 From: Dingcheng Luo <44854707+dc-luo@users.noreply.github.com> Date: Sun, 9 Jul 2023 08:23:24 -0500 Subject: [PATCH] Sphinx dev (#7) * sphinx documentation quickstart * updated autodocs * updated documentation * updated readme * updated readme * updated readme again * more documentation updates * updated documentation * added docs * updated docs * updated opt docs * updated read me and markdown imports * Added installation instructions and new theme * changes to the readme and instructions * updates to docs and added docs for utils --- INSTALL.md | 44 +++ README.md | 21 +- doc/Makefile | 20 ++ doc/conf.py | 50 +++ doc/index.rst | 19 + doc/install.rst | 2 + doc/make.bat | 35 ++ doc/soupy.collectives.rst | 10 + doc/soupy.modeling.rst | 87 +++++ doc/soupy.optimization.rst | 34 ++ doc/soupy.rst | 21 ++ doc/soupy.utils.rst | 36 ++ soupy/collectives/collective.py | 4 +- soupy/modeling/PDEControlProblem.py | 14 + soupy/modeling/augmentedVector.py | 14 +- soupy/modeling/controlCostFunctional.py | 107 +++++- soupy/modeling/controlModel.py | 325 +++++++++--------- soupy/modeling/controlQoI.py | 251 ++++++++++---- soupy/modeling/meanVarRiskMeasure.py | 264 +++++++------- soupy/modeling/meanVarRiskMeasureSAA.py | 116 +++++-- soupy/modeling/penalization.py | 125 ++++++- soupy/modeling/smoothPlusApproximation.py | 9 + soupy/modeling/superquantileRiskMeasureSAA.py | 65 ++-- soupy/modeling/variables.py | 5 +- soupy/optimization/bfgs.py | 62 ++-- soupy/optimization/projectableConstraint.py | 37 +- soupy/optimization/sgd.py | 52 +-- soupy/optimization/steepestDescent.py | 42 ++- soupy/utils/penalizationFiniteDifference.py | 16 +- soupy/utils/qoiFiniteDifference.py | 19 +- soupy/utils/scipyCostWrapper.py | 38 +- soupy/utils/stochasticCostFiniteDifference.py | 6 + 32 files changed, 1423 insertions(+), 527 deletions(-) create mode 100644 INSTALL.md create mode 100644 doc/Makefile create mode 100644 doc/conf.py create mode 100644 doc/index.rst create mode 100644 doc/install.rst create mode 100644 doc/make.bat create mode 100644 doc/soupy.collectives.rst create mode 100644 doc/soupy.modeling.rst create mode 100644 doc/soupy.optimization.rst create mode 100644 doc/soupy.rst create mode 100644 doc/soupy.utils.rst diff --git a/INSTALL.md b/INSTALL.md new file mode 100644 index 0000000..3aa05d8 --- /dev/null +++ b/INSTALL.md @@ -0,0 +1,44 @@ +# Installation + +## Dependencies +`SOUPy` depends on [FEniCS](http://fenicsproject.org/) version 2019.1. and [hIPPYlib](https://hippylib.github.io/) version 3.0.0 or above. + +`FEniCS` needs to be built with the following dependencies enabled: + + - `numpy`, `scipy`, `matplotlib`, `mpi4py` + - `PETSc` and `petsc4py` (version 3.10.0 or above) + - `SLEPc` and `slepc4py` (version 3.10.0 or above) + - PETSc dependencies: `parmetis`, `scotch`, `suitesparse`, `superlu_dist`, `ml`, `hypre` + - (optional): `gmsh`, `mshr`, `jupyter` + +## Recommended installation procedures using conda +1. Install `FEniCS` with appropriate dependencies +``` +conda create -n soupy -c conda-forge fenics==2019.1.0 matplotlib scipy jupyter +``` + +2. Clone the `hIPPYlib` [repository](https://github.com/hippylib/hippylib). +``` +git clone https://github.com/hippylib/hippylib.git +``` + +3. Clone the `SOUPy` [repository](https://github.com/hippylib/soupy/tree/main) +``` +git clone https://github.com/hippylib/soupy.git +``` + +4. Set the path to the `hIPPYlib` and `SOUPy` as environment variables, e.g. +``` +conda activate soupy +conda env config vars set HIPPYLIB_PATH=path/to/hippylib +conda env config vars set SOUPY_PATH=path/to/soupy +``` + +Examples in the `applications` directory can now be run. We refer to the full `FEniCS` [installation instructions](https://hippylib.readthedocs.io/en/3.0.0/installation.html) from `hIPPYlib` for more detail. + +## Build the SOUPy documentation using Sphinx + +Documentation for `SOUPy` can be built using `sphinx`, along with extensions +`myst_parser` and `sphinx_rtd_theme`. These can be installed via `pip`. + +To build simply run `make html` from `doc` folder. diff --git a/README.md b/README.md index 9e4bf01..2e32963 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,19 @@ -# soupy -Stochastic Optimization under Uncertainty in Python. +# Stochastic Optimization under high-dimensional Uncertainty in Python + +**S**tochastic **O**ptimization under high-dimensional **U**ncertainty in **Py**thon—**SOUPy**, +implements scalable algorithms to solve problems of PDE-constrained optimization under uncertainty, with the computational complexity measured in terms of PDE solves independent of the uncertain parameter dimension and optimization variable dimension. + +`SOUPy` is built on the open-source `hIPPYlib` library, which provides state-of-the-art scalable adjoint-based methods for deterministic and Bayesian inverse problems governed by PDEs, which in turn makes use of the `FEniCS` library for high-level formulation, discretization, and scalable solution of PDEs. + +`SOUPy` implements algorithms for the optimization of probability/risk measures such as mean, variance, and superquantile/condition-value-at-risk, subject to PDE constraints. +Sample based approximations of risk measures are supported in `SOUPy`, which leverages MPI for rapid parallel sampling. +Numerical optimization algorithms can be called from `SciPy` or using the built-in optimization algorithms. + +`SOUPy` has been developed and in active development to incorporate advanced approximation algorithms and capabilities, including: +- PDE-constrained operator/tensor/matrix products, +- symbolic differentiation (of appropriate Lagrangians) for the derivation of high order mixed derivatives (via the `FEniCS` interface), +- randomized algorithms for matrix and high order tensor decomposition, +- decomposition of uncertain parameter spaces by mixture models, +- Taylor expansion-based high-dimensional control variates, and +- product convolution approximations, +- common interfaces for random fields, PDE models, probability/risk measures, and control/design/inversion constraints. diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/conf.py b/doc/conf.py new file mode 100644 index 0000000..22681b2 --- /dev/null +++ b/doc/conf.py @@ -0,0 +1,50 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +basedir = os.path.abspath(os.path.join(os.path.dirname(__file__), '..')) +sys.path.insert(0, basedir) + +autodoc_mock_imports = ["dolfin", "matplotlib", "ffc", "ufl", + "petsc4py", "mpi4py", "scipy", "numpy", "hippylib"] +autodoc_default_flags = ['members', 'private-members', 'undoc-members'] +autoclass_content = 'both' + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +project = 'SOUPy' +copyright = "2023, Peng Chen, Dingcheng Luo, Thomas O'Leary-Roseberry, Umberto Villa" +author = "Peng Chen, Dingcheng Luo, Thomas O'Leary-Roseberry, Umberto Villa" +release = '0.1.0' + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = ['sphinx.ext.autodoc', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', + 'sphinx.ext.napoleon', 'myst_parser', 'sphinx.ext.autosummary'] + +templates_path = ['_templates'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] +source_suffix = ['.rst', '.md'] +master_doc = 'index' + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = 'sphinx_rtd_theme' +html_static_path = [] +# html_theme_options = { + # "nosidebar": True +# } +# html_theme_options = { +# "nosidebar": False, "sidebarwidth": "40em" +# } diff --git a/doc/index.rst b/doc/index.rst new file mode 100644 index 0000000..3d70ecb --- /dev/null +++ b/doc/index.rst @@ -0,0 +1,19 @@ +.. SOUPy documentation master file, created by + sphinx-quickstart on Thu Jul 6 18:03:26 2023. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to SOUPy's documentation! +================================= +.. include:: ../README.md + :parser: myst_parser.sphinx_ + +Documentation +--------------------- + +.. toctree:: + :titlesonly: + :maxdepth: 1 + + install + Complete API reference diff --git a/doc/install.rst b/doc/install.rst new file mode 100644 index 0000000..30e6af6 --- /dev/null +++ b/doc/install.rst @@ -0,0 +1,2 @@ +.. include:: ../INSTALL.md + :parser: myst_parser.sphinx_ \ No newline at end of file diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000..32bb245 --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/doc/soupy.collectives.rst b/doc/soupy.collectives.rst new file mode 100644 index 0000000..b0aba79 --- /dev/null +++ b/doc/soupy.collectives.rst @@ -0,0 +1,10 @@ +soupy.collectives +========================== + +soupy.collectives.collective +----------------------------------- + +.. automodule:: soupy.collectives.collective + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/doc/soupy.modeling.rst b/doc/soupy.modeling.rst new file mode 100644 index 0000000..57d7b0b --- /dev/null +++ b/doc/soupy.modeling.rst @@ -0,0 +1,87 @@ +soupy.modeling +========================= + +soupy.modeling.PDEControlProblem +--------------------------------------- + +.. automodule:: soupy.modeling.PDEControlProblem + :members: + :undoc-members: + :show-inheritance: + + +soupy.modeling.controlModel +--------------------------------------- + +.. automodule:: soupy.modeling.controlModel + :members: + :undoc-members: + :show-inheritance: + + +soupy.modeling.controlQoI +--------------------------------------- + +.. automodule:: soupy.modeling.controlQoI + :members: + :undoc-members: + :show-inheritance: + + + +soupy.modeling.controlCostFunctional +--------------------------------------- + +.. automodule:: soupy.modeling.controlCostFunctional + :members: + :undoc-members: + :show-inheritance: + +soupy.modeling.meanVarRiskMeasure +--------------------------------------- + +.. automodule:: soupy.modeling.meanVarRiskMeasure + :members: + :undoc-members: + :show-inheritance: + +soupy.modeling.meanVarRiskMeasureSAA +--------------------------------------- + +.. automodule:: soupy.modeling.meanVarRiskMeasureSAA + :members: + :undoc-members: + :show-inheritance: + +soupy.modeling.variables +------------------------------------------- + +.. automodule:: soupy.modeling.variables + :members: + :undoc-members: + :show-inheritance: + +soupy.modeling.augmentedVector +------------------------------------------- + +.. automodule:: soupy.modeling.augmentedVector + :members: + :undoc-members: + :show-inheritance: + +soupy.modeling.superquantileRiskMeasureSAA +------------------------------------------- + +.. automodule:: soupy.modeling.superquantileRiskMeasureSAA + :members: + :undoc-members: + :show-inheritance: + + +soupy.modeling.penalization +--------------------------------------- + +.. automodule:: soupy.modeling.penalization + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/doc/soupy.optimization.rst b/doc/soupy.optimization.rst new file mode 100644 index 0000000..360122c --- /dev/null +++ b/doc/soupy.optimization.rst @@ -0,0 +1,34 @@ +soupy.optimization +========================== + +soupy.optimization.bfgs +----------------------------------- + +.. automodule:: soupy.optimization.bfgs + :members: + :undoc-members: + :show-inheritance: + +soupy.optimization.projectableConstraint +----------------------------------------- + +.. automodule:: soupy.optimization.projectableConstraint + :members: + :undoc-members: + :show-inheritance: + +soupy.optimization.steepestDescent +----------------------------------- + +.. automodule:: soupy.optimization.steepestDescent + :members: + :undoc-members: + :show-inheritance: + +soupy.optimization.sgd +----------------------------------- + +.. automodule:: soupy.optimization.sgd + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/doc/soupy.rst b/doc/soupy.rst new file mode 100644 index 0000000..0e58ee4 --- /dev/null +++ b/doc/soupy.rst @@ -0,0 +1,21 @@ +soupy package +================ + +Submodules +----------- + +.. toctree:: + + soupy.modeling + soupy.optimization + soupy.collectives + soupy.utils + +Module contents +--------------- + +.. automodule:: soupy + :members: + :undoc-members: + :show-inheritance: + diff --git a/doc/soupy.utils.rst b/doc/soupy.utils.rst new file mode 100644 index 0000000..38791f4 --- /dev/null +++ b/doc/soupy.utils.rst @@ -0,0 +1,36 @@ +soupy.utils module +========================= + +soupy.utils.penalizationFiniteDifference +----------------------------------------------- + +.. automodule:: soupy.utils.penalizationFiniteDifference + :members: + :undoc-members: + :show-inheritance: + + +soupy.utils.qoiFiniteDifference +----------------------------------------------- + +.. automodule:: soupy.utils.qoiFiniteDifference + :members: + :undoc-members: + :show-inheritance: + + +soupy.utils.stochasticCostFiniteDifference +----------------------------------------------- + +.. automodule:: soupy.utils.stochasticCostFiniteDifference + :members: + :undoc-members: + :show-inheritance: + +soupy.utils.scipyCostWrapper module +--------------------------------------- + +.. automodule:: soupy.utils.scipyCostWrapper + :members: + :undoc-members: + :show-inheritance: diff --git a/soupy/collectives/collective.py b/soupy/collectives/collective.py index dd6c28d..554e886 100644 --- a/soupy/collectives/collective.py +++ b/soupy/collectives/collective.py @@ -46,7 +46,9 @@ class MultipleSamePartitioningPDEsCollective: """ def __init__(self, comm, is_serial_check=False): """ - :code:`comm` is :code:`mpi4py.MPI` comm + Constructor: + :param comm: MPI communicator + :type comm : :py:class:`mpi4py.MPI.Comm` """ self.comm = comm self.is_serial_check = is_serial_check diff --git a/soupy/modeling/PDEControlProblem.py b/soupy/modeling/PDEControlProblem.py index 3244689..e5d8f42 100644 --- a/soupy/modeling/PDEControlProblem.py +++ b/soupy/modeling/PDEControlProblem.py @@ -24,6 +24,20 @@ class PDEVariationalControlProblem(hp.PDEVariationalProblem): def __init__(self, Vh, varf_handler, bc, bc0, is_fwd_linear = False, lu_method="mumps"): + """ + Constructor + + :param Vh: List of function spaces the state, parameter, adjoint, and control + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param varf_handler: Variational form handler with :code:`__call__` method + :param bc: List of Dirichlet boundary conditions for the state + :param bc0: List of zeroed Dirichlet boundary conditions + :param is_fwd_linear: Flag indicating whether the forward problem is linear + :type is_fwd_linear: bool + :param lu_method: Method for solving linear systems (default, mumps, etc.) + :type lu_method: str + """ + # assert for class assumptions here assert id(Vh[STATE]) == id(Vh[ADJOINT]), print('Need to have same STATE and ADJOINT spaces') assert len(Vh) == 4 diff --git a/soupy/modeling/augmentedVector.py b/soupy/modeling/augmentedVector.py index 53696af..bcac20d 100644 --- a/soupy/modeling/augmentedVector.py +++ b/soupy/modeling/augmentedVector.py @@ -16,12 +16,18 @@ class AugmentedVector: """ - Class representing an augmented optimization variable (z, t), where t - is a real number. Methods mirror that of dl.Vector() - Assumes the last element in the array is the t variable + Class representing an augmented optimization variable :math:`(z, t)`, \ + where :math:`t` is a real number. \ + Methods mirror that of :code:`dolfin.Vector`. \ + Assumes the last element in the array is the :code:`t` variable + + .. note:: Currently supports only sample parallelism (i.e. serial mesh) """ def __init__(self, v, copy_vector=True): - # Currently worked out only for serial mesh (allow sample parallel) + """ + :param v: :code:`dolfin.Vector` to be augmented + :param copy_vector: If :code:`True`, copy the vector, otherwise use the same memory + """ assert v.mpi_comm().Get_size() == 1 if copy_vector: self.v = v.copy() diff --git a/soupy/modeling/controlCostFunctional.py b/soupy/modeling/controlCostFunctional.py index c122ef3..cac1278 100644 --- a/soupy/modeling/controlCostFunctional.py +++ b/soupy/modeling/controlCostFunctional.py @@ -50,17 +50,22 @@ class DeterministicControlCostFunctional(ControlCostFunctional): """ This class implements a deterministic approximation for the optimal control problem under uncertainty by considering a fixed parameter at the mean of the prior + + .. math:: J(z) := Q(\\bar{m}, z) + P(z) + """ def __init__(self, model, prior, penalization=None): """ - Code implements the deterministic approximation of an optimal control problem - where the random parameter is taken as its mean value - min J(z) := Q(\bar{m}) + P(z) - - :code: `model` an instance of `soupy.ControlModel`, which contains - the PDE problem and qoi - - :code: `prior` a prior for the random parameter. - - :code: `penalization` an optional Penalization object + Constructor + + :param model: control model containing the :code:`soupy.PDEVariationalControlProblem` + and :code:`soupy.ControlQoI` + :type model: :py:class:`soupy.ControlModel` + :param prior: The prior distribution for the random parameter + :type prior: :py:class:`hippylib.Prior` + :param penalization: An optional penalization object for the cost of control + :type penalization: :py:class:`soupy.Penalization` """ self.model = model self.prior = prior @@ -99,11 +104,13 @@ def objective(self, z): def computeComponents(self, z, order=0): """ - Computes the components for the stochastic approximation of the cost - Parameters: - - :code: `z` the control variable - - :code: `order` the order of derivatives needed. - 0 for cost. 1 for grad. 2 for Hessian + Computes the components for the evaluation of the cost functional + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed. + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int """ # Check if a new control variable is used @@ -142,6 +149,17 @@ def computeComponents(self, z, order=0): def cost(self, z, order=0, **kwargs): + """ + Computes the value of the cost functional at given control :math:`z` + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed after evaluation + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + + :return: Value of the cost functional + """ self.computeComponents(z, order=order) objective = self.model.cost(self.x) if self.penalization is None: @@ -152,8 +170,16 @@ def cost(self, z, order=0, **kwargs): def costGrad(self, z, out): """ - Compute the gradient - Assumes self.costValue(z, order=1) has been called + Computes the value of the cost functional at given control :math:`z` + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param out: (Dual of) the gradient of the cost functional + :type out: :py:class:`dolfin.Vector` + + :return: the norm of the gradient in the correct inner product :math:`(g_z,g_z)_{Z}^{1/2}` + + .. note:: Assumes :code:`self.cost` has been called with :code:`order >= 1` """ self.model.evalGradientControl(self.x, out) if self.penalization is not None: @@ -165,10 +191,17 @@ def costGrad(self, z, out): def costHessian(self, z, zhat, out): """ - Apply the the reduced Hessian to the vector :code:`zhat` - evaluated at control variable :code:`z` - Return the result in :code:`out`. - Need to call self.costValue(self, order) with order>=1 before self.costHessian + Apply the the reduced Hessian to the vector :math:`zhat` + evaluated at control variable :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` + :param zhat: The direction of Hessian action + :type zhat: :py:class:`dolfin.Vector` + :param out: The assembled Hessian action + :type out: :py:class:`dolfin.Vector` + + .. note:: Assumes :code:`self.cost` has been called with :code:`order >= 2` """ self.model.setPointForHessianEvaluations(self.x) self.model.applyCz(zhat, self.rhs_fwd) @@ -193,7 +226,22 @@ def costHessian(self, z, zhat, out): class RiskMeasureControlCostFunctional: + """ + This class implements a risk measure cost functional for \ + optimal control problem under uncertainty + + .. math:: J(z) := \\rho[Q(m, z)] + P(z) + + """ def __init__(self, risk_measure, penalization=None): + """ + Constructor + + :param risk_measure: Class implementing the risk measure :math:`\\rho(m,z)` + :type risk_measure: :py:class:`soupy.RiskMeasure` + :param penalization: An optional penalization object for the cost of control + :type penalization: :py:class:`soupy.Penalization` + """ self.risk_measure = risk_measure self.penalization = penalization self.grad_risk = self.risk_measure.generate_vector(CONTROL) @@ -203,6 +251,18 @@ def generate_vector(self, component="ALL"): return self.risk_measure.generate_vector(component) def cost(self, z, order=0, **kwargs): + """ + Computes the value of the cost functional at given control :math:`z` + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed after evaluation + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + :param kwargs: additional arguments, e.g. :code:`rng` for the risk measure computation + + :return: Value of the cost functional + """ self.risk_measure.computeComponents(z, order=order, **kwargs) cost_risk = self.risk_measure.cost() if self.penalization is not None: @@ -214,7 +274,16 @@ def cost(self, z, order=0, **kwargs): def costGrad(self, z, out): """ - First calls cost with order = 1 + Computes the value of the cost functional at given control :math:`z` + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param out: (Dual of) the gradient of the cost functional + :type out: :py:class:`dolfin.Vector` + + :return: the norm of the gradient in the correct inner product :math:`(g_z,g_z)_{Z}^{1/2}` + + .. note:: Assumes :code:`self.cost` has been called with :code:`order >= 1` """ out.zero() self.risk_measure.computeComponents(z, order=1) diff --git a/soupy/modeling/controlModel.py b/soupy/modeling/controlModel.py index 0dce8c0..d8a0d7e 100644 --- a/soupy/modeling/controlModel.py +++ b/soupy/modeling/controlModel.py @@ -26,22 +26,24 @@ class ControlModel: """ This class contains the structure needed to evaluate the control objective - As inputs it takes a :code:`PDEProblem object`, and a :code:`Qoi` object. + As inputs it takes a :py:class:`PDEVariationalControlProblem`, and a :code:`Qoi` object. In the following we will denote with - - :code:`u` the state variable - - :code:`m` the (model) parameter variable - - :code:`p` the adjoint variable - - :code:`z` the control variable - + * :code:`u` the state variable + * :code:`m` the (model) parameter variable + * :code:`p` the adjoint variable + * :code:`z` the control variable """ def __init__(self, problem, qoi): """ - Create a model given: - - problem: the description of the forward/adjoint problem and all the sensitivities - - qoi: the qoi component of the cost functional + Constructor + + :param problem: The PDE problem + :type problem: :py:class:`PDEVariationalControlProblem` + :param qoi: The quantity of interest + :type qoi: :py:class:`ControlQoI` """ self.problem = problem self.qoi = qoi @@ -54,21 +56,15 @@ def __init__(self, problem, qoi): def generate_vector(self, component = "ALL"): """ - By default, return the list :code:`[u,m,p,z]` where: - - - :code:`u` is any object that describes the state variable - - :code:`m` is a :code:`dolfin.Vector` object that describes the parameter variable. \ - (Needs to support linear algebra operations) - - :code:`p` is any object that describes the adjoint variable - - :code:`z` is any object that describes the control variable - - If :code:`component = STATE` return only :code:`u` - - If :code:`component = PARAMETER` return only :code:`m` - - If :code:`component = ADJOINT` return only :code:`p` + :param component: The component of the vector to generate (:code:`soupy.STATE`, \ + :code:`soupy.PARAMETER`, :code:`soupy.ADJOINT`, :code:`soupy.CONTROL`, or :code:`"ALL"`) - If :code:`component = CONTROL` return only :code:`z` + :return: By default, :code:`component == "ALL"` will return the list :code:`[u,m,p,z]` + where each element is a :py:class:`dolfin.Vector` object of the appropriate size + If :code:`component = soupy.STATE` return only :code:`u`. + If :code:`component = soupy.PARAMETER` return only :code:`m`. + If :code:`component = soupy.ADJOINT` return only :code:`p`. + If :code:`component = soupy.CONTROL` return only :code:`z`. """ if component == "ALL": x = [self.problem.generate_state(), @@ -100,31 +96,32 @@ def init_control(self, z): def cost(self, x): """ - Given the list :code:`x = [u,m,p]` which describes the state, parameter, and - adjoint variable compute the cost functional as the sum of - the qoi functional and the regularization functional. - - Return the list [cost functional, regularization functional, qoi functional] - + Evaluate the QoI at the a given point :math:`q(u,m,z)`. + + :param x: The point :code:`x = [u,m,p,z]` at which to evaluate the QoI + :type x: list of :py:class:`dolfin.Vector` objects + + :return: The value of the QoI + .. note:: :code:`p` is not needed to compute the cost functional """ qoi_cost = self.qoi.cost(x) - # reg_cost = self.prior.cost(x[PARAMETER]) - # return [qoi_cost+reg_cost, reg_cost, qoi_cost] return qoi_cost def solveFwd(self, out, x): """ Solve the (possibly non-linear) forward problem. - - Parameters: - - :code:`out`: is the solution of the forward problem (i.e. the state) (Output parameters) - - :code:`x = [u,m,p]` provides - 1) the parameter variable :code:`m` for the solution of the forward problem - 2) the initial guess :code:`u` if the forward problem is non-linear - - .. note:: :code:`p` is not accessed. + + :param out: Solution of the forward problem (state) + :type out: :py:class:`dolfin.Vector` + :param x: The point :code:`x = [u,m,p,z]`. Provides + the parameter and control variables :code:`m` + and :code:`z` for the solution of the forward problem + and the initial guess :code:`u` if the forward problem is non-linear + :type: list of :py:class:`dolfin.Vector` objects + + .. note:: :code:`p` is not accessed. """ self.n_fwd_solve = self.n_fwd_solve + 1 self.problem.solveFwd(out, x) @@ -133,12 +130,16 @@ def solveFwd(self, out, x): def solveAdj(self, out, x): """ Solve the linear adjoint problem. - Parameters: - - :code:`out`: is the solution of the adjoint problem (i.e. the adjoint :code:`p`) (Output parameter) - - :code:`x = [u, m, p]` provides - 1) the parameter variable :code:`m` for assembling the adjoint operator - 2) the state variable :code:`u` for assembling the adjoint right hand side - .. note:: :code:`p` is not accessed + + :param out: Solution of the forward problem (state) + :type out: :py:class:`dolfin.Vector` + :param x: The point :code:`x = [u,m,p,z]`. Provides + the state, parameter and control variables :code:`u`, :code:`m`, + and :code:`z` for the solution assembling the adjoint operator. + Vector :code:`u` is also used to assemble the adjoint right hand side. + :type: list of :py:class:`dolfin.Vector` objects + + .. note:: :code:`p` is not accessed """ self.n_adj_solve = self.n_adj_solve + 1 rhs = self.problem.generate_state() @@ -151,36 +152,38 @@ def solveAdj(self, out, x): def evalGradientParameter(self,x, mg): """ - Evaluate the gradient for the variational parameter equation at the point :code:`x=[u,m,p]`. - Parameters: - - :code:`x = [u,m,p]` the point at which to evaluate the gradient. - - :code:`mg` the variational gradient :math:`(g, mtest)`, mtest being a test function in the parameter space \ - (Output parameter) - - Returns the norm of the gradient in the correct inner product :math:`g_norm = sqrt(g,g)` + Evaluate the :math:`m` gradient action form at the point :math:`(u,m,p,z)` + + :param x: The point :code:`x = [u,m,p,z]`. Provides + the state, parameter, adjoint, and control variables :code:`u`, :code:`m`, :code:`p`, + and :code:`z` for the assembling the gradient action. + Vector :code:`u` is also used to assemble the adjoint right hand side. + :type x: list of :py:class:`dolfin.Vector` objects + :param mg: Dual of the gradient with respect to the parameter i.e. :math:`(g_m, m_{\mathrm{test}})_{M}` + :type mg: :py:class:`dolfin.Vector` + + :return: the norm of the gradient in the correct inner product :math:`(g_m,g_m)_{M}^{1/2}` """ tmp = self.generate_vector(PARAMETER) self.problem.evalGradientParameter(x, mg) self.qoi.grad(PARAMETER,x,tmp) mg.axpy(1., tmp) - # if not qoi_only: - # self.prior.grad(x[PARAMETER], tmp) - # mg.axpy(1., tmp) - - # self.prior.Msolver.solve(tmp, mg) - #self.prior.Rsolver.solve(tmp, mg) return math.sqrt(mg.inner(tmp)) def evalGradientControl(self,x, mg): """ - Evaluate the gradient for the variational parameter equation at the point :code:`x=[u,m,p]`. - Parameters: - - :code:`x = [u,m,p]` the point at which to evaluate the gradient. - - :code:`mg` the variational gradient :math:`(g, mtest)`, mtest being a test function in the parameter space \ - (Output parameter) - - Returns the norm of the gradient in the correct inner product :math:`g_norm = sqrt(g,g)` + Evaluate the :math:`z` gradient action form at the point :math:`(u,m,p,z)` + + :param x: The point :code:`x = [u,m,p,z]`. Provides + the state, parameter, adjoint, and control variables :code:`u`, :code:`m`, :code:`p`, + and :code:`z` for the assembling the gradient action + Vector :code:`u` is also used to assemble the adjoint right hand side. + :type x: list of :py:class:`dolfin.Vector` objects + :param mg: Dual of the gradient with respect to the control i.e. :math:`(g_z, z_{\mathrm{test}})_{Z}` + :type mg: :py:class:`dolfin.Vector` + + :return: the norm of the gradient in the correct inner product :math:`(g_z,g_z)_{Z}^{1/2}` """ tmp = self.generate_vector(CONTROL) self.problem.evalGradientControl(x, mg) @@ -194,15 +197,17 @@ def evalGradientControl(self,x, mg): def setPointForHessianEvaluations(self, x, gauss_newton_approx=False): """ - Specify the point :code:`x = [u,m,p]` at which the Hessian operator (or the Gauss-Newton approximation) - needs to be evaluated. - Parameters: - - :code:`x = [u,m,p]`: the point at which the Hessian or its Gauss-Newton approximation needs to be evaluated. - - :code:`gauss_newton_approx (bool)`: whether to use Gauss-Newton approximation (default: use Newton) + Specify the point :code:`x = [u,m,p,z]` at which the Hessian operator \ + (or the Gauss-Newton approximation) needs to be evaluated. + + :param x: The point :code:`x = [u,m,p,z]` for which the Hessian needs to be evaluated + :type xk: list of :py:class:`dolfin.Vector` objects + :param gauss_newton_approx: whether to use the Gauss-Newton approximation (default: use Newton) + :type gauss_newton_approx: bool .. note:: This routine should either: - - simply store a copy of x and evaluate action of blocks of the Hessian on the fly - - or partially precompute the block of the hessian (if feasible) + 1. simply store a copy of x and evaluate action of blocks of the Hessian on the fly, or + 2. partially precompute the block of the hessian (if feasible) """ self.gauss_newton_approx = gauss_newton_approx self.problem.setLinearizationPoint(x, self.gauss_newton_approx) @@ -214,9 +219,11 @@ def setPointForHessianEvaluations(self, x, gauss_newton_approx=False): def solveFwdIncremental(self, sol, rhs): """ Solve the linearized (incremental) forward problem for a given right-hand side - Parameters: - - :code:`sol` the solution of the linearized forward problem (Output) - - :code:`rhs` the right hand side of the linear system + + :param sol: Solution of the incremental forward problem (state) + :type sol: :py:class:`dolfin.Vector` + :param rhs: Right hand side of the linear system + :type rhs: :py:class:`dolfin.Vector` """ self.n_inc_solve = self.n_inc_solve + 1 self.problem.solveIncremental(sol,rhs, False) @@ -225,9 +232,11 @@ def solveFwdIncremental(self, sol, rhs): def solveAdjIncremental(self, sol, rhs): """ Solve the incremental adjoint problem for a given right-hand side - Parameters: - - :code:`sol` the solution of the incremental adjoint problem (Output) - - :code:`rhs` the right hand side of the linear system + + :param sol: Solution of the incremental adjoint problem (adjoint) + :type sol: :py:class:`dolfin.Vector` + :param rhs: Right hand side of the linear system + :type rhs: :py:class:`dolfin.Vector` """ self.n_inc_solve = self.n_inc_solve + 1 self.problem.solveIncremental(sol,rhs, True) @@ -235,25 +244,27 @@ def solveAdjIncremental(self, sol, rhs): def applyC(self, dm, out): """ - Apply the :math:`C` block of the Hessian to a (incremental) parameter variable, i.e. - :code:`out` = :math:`C dm` + Apply the :math:`C_{m}` block of the Hessian to a (incremental) parameter variable, i.e. + :code:`out` = :math:`C_{m} dm` - Parameters: - - :code:`dm` the (incremental) parameter variable - - :code:`out` the action of the :math:`C` block on :code:`dm` - + :param dm: The (incremental) parameter variable + :type dm: :py:class:`dolfin.Vector` + :param out: The action of the :math:`C_z` block on :code:`dm` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ self.problem.apply_ij(ADJOINT,PARAMETER, dm, out) def applyCz(self, dz, out): """ - Apply the :math:`C` block of the Hessian to a (incremental) control variable, i.e. + Apply the :math:`C_z` block of the Hessian to a (incremental) control variable, i.e. :code:`out` = :math:`C_z dz` - Parameters: - - :code:`dz` the (incremental) control variable - - :code:`out` the action of the :math:`C` block on :code:`dm` + :param dz: The (incremental) control variable + :type dz: :py:class:`dolfin.Vector` + :param out: The action of the :math:`C_z` block on :code:`dz` + :type out: :py:class:`dolfin.Vector` .. note:: This routine assumes that :code:`out` has the correct shape. """ @@ -261,24 +272,28 @@ def applyCz(self, dz, out): def applyCt(self, dp, out): """ - Apply the transpose of the :math:`C` block of the Hessian to a (incremental) adjoint variable. - :code:`out` = :math:`C^t dp` - Parameters: - - :code:`dp` the (incremental) adjoint variable - - :code:`out` the action of the :math:`C^T` block on :code:`dp` - + Apply the transpose of the :math:`C_{m}` block of the Hessian to a (incremental) adjoint variable. + :code:`out` = :math:`C_{m}^T dp` + + :param dp: The (incremental) adjoint variable + :type dp: :py:class:`dolfin.Vector` + :param out: The action of the :math:`C_{m}^T` block on :code:`dp` + :type out: :py:class:`dolfin.Vector` + ..note:: This routine assumes that :code:`out` has the correct shape. """ self.problem.apply_ij(PARAMETER,ADJOINT, dp, out) def applyCzt(self, dp, out): """ - Apply the transpose of the :math:`C_z` block of the Hessian to a (incremental) adjoint variable. - :code:`out` = :math:`C_z^t dp` - Parameters: - - :code:`dp` the (incremental) adjoint variable - - :code:`out` the action of the :math:`C_z^T` block on :code:`dp` - + Apply the transpose of the :math:`C_{z}` block of the Hessian to a (incremental) adjoint variable. + :code:`out` = :math:`C_{z}^T dp` + + :param dp: The (incremental) adjoint variable + :type dp: :py:class:`dolfin.Vector` + :param out: The action of the :math:`C_{z}^T` block on :code:`dp` + :type out: :py:class:`dolfin.Vector` + ..note:: This routine assumes that :code:`out` has the correct shape. """ self.problem.apply_ij(CONTROL, ADJOINT, dp, out) @@ -287,11 +302,12 @@ def applyWuu(self, du, out): """ Apply the :math:`W_{uu}` block of the Hessian to a (incremental) state variable. :code:`out` = :math:`W_{uu} du` - - Parameters: - - :code:`du` the (incremental) state variable - - :code:`out` the action of the :math:`W_{uu}` block on :code:`du` - + + :param du: The (incremental) state variable + :type du: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{uu}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ self.qoi.apply_ij(STATE,STATE, du, out) @@ -307,11 +323,12 @@ def applyWum(self, dm, out): """ Apply the :math:`W_{um}` block of the Hessian to a (incremental) parameter variable. :code:`out` = :math:`W_{um} dm` - - Parameters: - - :code:`dm` the (incremental) parameter variable - - :code:`out` the action of the :math:`W_{um}` block on :code:`du` - + + :param dm: The (incremental) parameter variable + :type dm: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{um}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -328,10 +345,11 @@ def applyWuz(self, dz, out): Apply the :math:`W_{uz}` block of the Hessian to a (incremental) control variable. :code:`out` = :math:`W_{uz} dz` - Parameters: - - :code:`dz` the (incremental) control variable - - :code:`out` the action of the :math:`W_{uz}` block on :code:`du` - + :param dz: The (incremental) control variable + :type dz: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{uz}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -346,11 +364,12 @@ def applyWmu(self, du, out): """ Apply the :math:`W_{mu}` block of the Hessian to a (incremental) state variable. :code:`out` = :math:`W_{mu} du` - - Parameters: - - :code:`du` the (incremental) state variable - - :code:`out` the action of the :math:`W_{mu}` block on :code:`du` - + + :param du: The (incremental) state variable + :type du: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{mu}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -366,10 +385,11 @@ def applyWzu(self, du, out): Apply the :math:`W_{zu}` block of the Hessian to a (incremental) state variable. :code:`out` = :math:`W_{zu} du` - Parameters: - - :code:`du` the (incremental) state variable - - :code:`out` the action of the :math:`W_{zu}` block on :code:`du` - + :param du: The (incremental) state variable + :type du: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{zu}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -381,41 +401,16 @@ def applyWzu(self, du, out): out.axpy(1., tmp) - # def applyR(self, dm, out): - # """ - # Apply the regularization :math:`R` to a (incremental) parameter variable. - # :code:`out` = :math:`R dm` - - # Parameters: - # - :code:`dm` the (incremental) parameter variable - # - :code:`out` the action of :math:`R` on :code:`dm` - - # .. note:: This routine assumes that :code:`out` has the correct shape. - # """ - # self.prior.R.mult(dm, out) - - - # def Rsolver(self): - # """ - # Return an object :code:`Rsovler` that is a suitable solver for the regularization - # operator :math:`R`. - - # The solver object should implement the method :code:`Rsolver.solve(z,r)` such that - # :math:`Rz \approx r`. - # """ - # return self.prior.Rsolver - - def applyWmm(self, dm, out): """ Apply the :math:`W_{mm}` block of the Hessian to a (incremental) parameter variable. :code:`out` = :math:`W_{mm} dm` - Parameters: - - - :code:`dm` the (incremental) parameter variable - - :code:`out` the action of the :math:`W_{mm}` on block :code:`dm` - + :param dm: The (incremental) parameter variable + :type dm: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{mm}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -431,11 +426,11 @@ def applyWzz(self, dz, out): Apply the :math:`W_{zz}` block of the Hessian to a (incremental) control variable. :code:`out` = :math:`W_{zz} dz` - Parameters: - - - :code:`dm` the (incremental) control variable - - :code:`out` the action of the :math:`W_{zz}` on block :code:`dz` - + :param dz: The (incremental) control variable + :type dz: :py:class:`dolfin.Vector` + :param out: The action of the :math:`W_{zz}` block on :code:`du` + :type out: :py:class:`dolfin.Vector` + .. note:: This routine assumes that :code:`out` has the correct shape. """ if self.gauss_newton_approx: @@ -447,6 +442,20 @@ def applyWzz(self, dz, out): out.axpy(1., tmp) def apply_ij(self, i, j, d, out): + """ + Apply the :math:`(i,j)` block of the Hessian to a vector :code:`d` + + :param i: The output variable index :code:`soupy.STATE`, :code:`soupy.PARAMETER`, + :code:`soupy.ADJOINT`, or :code:`soupy.CONTROL` + :type i: int + :param j: The input variable index :code:`soupy.STATE`, :code:`soupy.PARAMETER`, + :code:`soupy.ADJOINT`, or :code:`soupy.CONTROL` + :type j: int + :param d: The vector to which the Hessian is applied + :type d: :py:class:`dolfin.Vector` + :param out: The action of the Hessian on :code:`d` + :type out: :py:class:`dolfin.Vector` + """ if i == STATE and j == STATE: self.applyWuu(d,out) elif i == STATE and j == PARAMETER: diff --git a/soupy/modeling/controlQoI.py b/soupy/modeling/controlQoI.py index 78e292d..029b738 100644 --- a/soupy/modeling/controlQoI.py +++ b/soupy/modeling/controlQoI.py @@ -21,25 +21,29 @@ class ControlQoI(object): """ - Abstract class to model the control quantity of a control problem - In the following :code:`x` will denote the variable :code:`[u, m, p, z]`, denoting respectively - the state :code:`u`, the parameter :code:`m`, the adjoint variable :code:`p`, and the control variable :code:`z` + Abstract class to define the optimization quantity of interest for the \ + optimal control problem under uncertainty \ + In the following :code:`x` will denote the variable :code:`[u, m, p, z]`, \ + denoting respectively the state :code:`u`, the parameter :code:`m`, \ + the adjoint variable :code:`p`, and the control variable :code:`z` - The methods in the class ControlQoI will usually access the state u and possibly the - parameter :code:`m` and control :code: `z`. The adjoint variables will never be accessed. + The methods in the class ControlQoI will usually access the state :code:`u` and possibly the \ + parameter :code:`m` and control :code:`z`. The adjoint variables will never be accessed. """ def cost(self,x): """ - Given x evaluate the cost functional. - Only the state u and (possibly) the parameter m are accessed. """ + Given :code:`x` evaluate the cost functional. Only the state :code:`u` \ + and (possibly) the parameter :code:`m` are accessed. """ raise NotImplementedError("Child class should implement method cost") def grad(self, i, x, out): """ - Given the state and the paramter in :code:`x`, compute the partial gradient of the misfit - functional in with respect to the state (:code:`i == STATE`) or with respect to the parameter (:code:`i == PARAMETER`). + Given the state and the paramter in :code:`x`, compute the partial gradient \ + of the QoI in with respect to the state (:code:`i == soupy.STATE`), \ + parameter (:code:`i == soupy.PARAMETER`), or \ + control (:code:`i == soupy.CONTROL`). """ raise NotImplementedError("Child class should implement method grad") @@ -47,16 +51,14 @@ def grad(self, i, x, out): def setLinearizationPoint(self,x, gauss_newton_approx=False): """ Set the point for linearization. - Inputs: - - :code:`x=[u, m, p]` - linearization point - :code:`gauss_newton_approx (bool)` - whether to use Gauss Newton approximation """ raise NotImplementedError("Child class should implement method setLinearizationPoint") def apply_ij(self,i,j, dir, out): """ - Apply the second variation :math:`\delta_{ij}` (:code:`i,j = STATE,PARAMETER`) of the cost in direction :code:`dir`. + Apply the second variation :math:`\delta_{ij}` (:code:`i,j` = :code:`soupy.STATE`, \ + :code:`soupy.PARAMETER`, :code:`soupy.CONTROL`) \ + of the cost in direction :code:`dir`. """ raise NotImplementedError("Child class should implement method apply_ij") @@ -65,13 +67,23 @@ def apply_ij(self,i,j, dir, out): class L2MisfitVarfHandler: """ - Form handler for the L2 Misfit + Form handler for the L^2 Misfit + + .. math:: \int_{\Omega} \chi (u - u_d)^2 dx + + where :math:`u_d` is the reference state \ + and :math:`\chi` is the characteristic function \ + defining the region of integration """ def __init__(self, ud, chi=None): """ - :code: `ud` is the reference function - :code: `chi` is a characteristic function defining region of integration + Constructor + + :param ud: The reference state + :type ud: :py:class:`dolfin.Function` or :py:class:`dolfin.Expression` + :param chi: The characteristic function defining the region of integration + :type chi: :py:class:`dolfin.Function` or :py:class:`dolfin.Expression` """ self.chi = chi self.ud = ud @@ -86,14 +98,19 @@ def __call__(self, u, m, z): class VariationalControlQoI(ControlQoI): """ - define the quantity of interest and its derivative information + Class for a QoI defined by its variational form """ def __init__(self, mesh, Vh, form_handler): """ - Constructor. - INPUTS: - - mesh: the mesh - - Vh: the finite element space for [state, parameter, adjoint, optimization] variable + Constructor + + :param mesh: The mesh object + :param Vh: List of function spaces for the state, parameter, + adjoint, and optimization variables + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param form_handler: The form handler for the variational form with a + :code:`__call__` method that takes as input the state, parameter, and control variables + as functions and returns the variational form """ self.mesh = mesh self.Vh = Vh @@ -106,9 +123,12 @@ def __init__(self, mesh, Vh, form_handler): def cost(self, x): """ - evaluate the qoi at given x - :param x: [state, parameter, adjoint, optimization] variable - :return: qoi(x) + Evaluate the qoi at given point :math:`q(u,m,z)` + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :return: QoI evaluated at x """ u = hp.vector2Function(x[STATE], self.Vh[STATE]) m = hp.vector2Function(x[PARAMETER], self.Vh[PARAMETER]) @@ -118,16 +138,32 @@ def cost(self, x): def adj_rhs(self, x, rhs): """ - The right hand for the adjoint problem (i.e. the derivative of the Lagrangian funtional - with respect to the state u). - INPUTS: - - x coefficient vector of all variables - - rhs: FEniCS vector to store the rhs for the adjoint problem. + The right hand for the adjoint problem (i.e. the derivative of the Lagrangian functional \ + with respect to the state u). + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :param rhs: The assembled rhs for the adjoint problem. + :type rhs: :py:class:`dolfin.Vector` """ self.grad(STATE, x, rhs) rhs *= -1 def grad(self, i, x, out): + """ + First variation of the QoI with respect to the :code:`i` th variable \ + where :code:`i` is either :code:`soupy.STATE`, :code:`soupy.PARAMETER`, \ + or :code:`soupy.CONTROL`. + + :param i: Index of the variable with respect to which the first variation is taken + :type i: int + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :param out: The assembled first variation + :type out: :py:class:`dolfin.Vector` + """ out.zero() u = hp.vector2Function(x[STATE], self.Vh[STATE]) m = hp.vector2Function(x[PARAMETER], self.Vh[PARAMETER]) @@ -140,12 +176,20 @@ def grad(self, i, x, out): def apply_ij(self,i,j, dir, out): """ - Apply the second variation \delta_ij (i,j = STATE,PARAMETER,CONTROL) of the q.o.i. in direction dir. - INPUTS: - - i,j integer (STATE=0, PARAMETER=1, CONTROL=3) which indicates with respect to which variables differentiate - - dir the direction in which to apply the second variation - - out: FEniCS vector to store the second variation in the direction dir. - NOTE: setLinearizationPoint must be called before calling this method. + Apply the second variation :math:`\\delta_ij` (:code:`i,j` = :code:`soupy.STATE`, \ + :code:`soupy.PARAMETER`, :code:`soupy.CONTROL`) \ + of the QoI in direction :code:`dir`. + + :param i: Index of the output variable + :type i: int + :param j: Index of the input variable + :type j: int + :param dir: The direction in which to apply the second variation + :type dir: :py:class:`dolfin.Vector` + :param out: The assembled second variation + :type out: :py:class:`dolfin.Vector` + + ..note:: :code:`setLinearizationPoint` must be called before calling this method. """ out.zero() @@ -161,14 +205,22 @@ def apply_ij(self,i,j, dir, out): def apply_ijk(self,i,j,k,dir1,dir2, out): """ - Apply the third order variation of the q.o.i. w.r.t. ijk in direction dir1, dir2 for j and k - :param i: STATE or PARAMETER or CONTROL - :param j: - :param k: - :param dir1: - :param dir2: - :param out: - :return: out + Apply the third order variation of the QoI in the \ + :code:`i` th, :code:`j` th, and :code:`k` th variables in directions \ + :code:`dir1` and :code:`dir2`. + + :param i: First variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type i: int + :param j: Second variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type j: int + :param k: Third variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type k: int + :param dir1: Direction for variable :code:`j` + :type dir1: :py:class:`dolfin.Vector` + :param dir2: Direction for variable :code:`k` + :type dir2: :py:class:`dolfin.Vector` + :param out: The assembled third variation + :type out: :py:class:`dolfin.Vector` """ out.zero() @@ -184,8 +236,10 @@ def apply_ijk(self,i,j,k,dir1,dir2, out): def setLinearizationPoint(self, x, gauss_newton_approx=False): """ Specify the linearization point for computation of the second variations in method apply_ij. - INPUTS: - - x = [u,m,p,z] is a list of the state u, parameter m, and adjoint variable p + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` """ for i in range(len(x)): self.x[i].zero() @@ -194,14 +248,22 @@ def setLinearizationPoint(self, x, gauss_newton_approx=False): class L2MisfitControlQoI(ControlQoI): """ - define the quantity of interest and its derivative information + Class for the :math:`L^2(\Omega)` misfit functional, + + .. math:: \\int_\\Omega (u - u_d)^2 dx + + where :math:`u_d` is the reference state. """ def __init__(self, mesh, Vh, ud): """ - Constructor. - INPUTS: - - mesh: the mesh - - Vh: the finite element space for [state, parameter, adjoint, optimization] variable + Constructor + + :param mesh: The mesh object + :param Vh: List of function spaces for the state, parameter, + adjoint, and optimization variables + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param ud: The reference state as a vector + :type ud: :py:class:`dolfin.Vector` """ self.mesh = mesh self.Vh = Vh @@ -220,10 +282,14 @@ def __init__(self, mesh, Vh, ud): def cost(self, x): """ - evaluate the qoi at given x - :param x: [state, parameter, adjoint, optimization] variable - :return: qoi(x) + Evaluate the qoi at given point :math:`q(u,m,z)` + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :return: QoI evaluated at x """ + self.diff.zero() self.diff.axpy(1.0, x[STATE]) self.diff.axpy(-1.0, self.ud) @@ -232,16 +298,33 @@ def cost(self, x): def adj_rhs(self, x, rhs): """ - The right hand for the adjoint problem (i.e. the derivative of the Lagrangian funtional - with respect to the state u). - INPUTS: - - x coefficient vector of all variables - - rhs: FEniCS vector to store the rhs for the adjoint problem. + The right hand for the adjoint problem (i.e. the derivative of the Lagrangian functional \ + with respect to the state u). + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :param rhs: The assembled rhs for the adjoint problem. + :type rhs: :py:class:`dolfin.Vector` """ self.grad(STATE, x, rhs) rhs *= -1 def grad(self, i, x, out): + """ + First variation of the QoI with respect to the :code:`i` th variable \ + where :code:`i` is either :code:`soupy.STATE`, :code:`soupy.PARAMETER`, \ + or :code:`soupy.CONTROL`. + + :param i: Index of the variable with respect to which the first variation is taken + :type i: int + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` + :param out: The assembled first variation + :type out: :py:class:`dolfin.Vector` + """ + out.zero() if i == STATE: self.diff.zero() @@ -254,12 +337,20 @@ def grad(self, i, x, out): def apply_ij(self,i,j, dir, out): """ - Apply the second variation \delta_ij (i,j = STATE,PARAMETER,CONTROL) of the q.o.i. in direction dir. - INPUTS: - - i,j integer (STATE=0, PARAMETER=1, CONTROL=3) which indicates with respect to which variables differentiate - - dir the direction in which to apply the second variation - - out: FEniCS vector to store the second variation in the direction dir. - NOTE: setLinearizationPoint must be called before calling this method. + Apply the second variation :math:`\\delta_ij` (:code:`i,j` = :code:`soupy.STATE`, \ + :code:`soupy.PARAMETER`, :code:`soupy.CONTROL`) \ + of the QoI in direction :code:`dir`. + + :param i: Index of the output variable + :type i: int + :param j: Index of the input variable + :type j: int + :param dir: The direction in which to apply the second variation + :type dir: :py:class:`dolfin.Vector` + :param out: The assembled second variation + :type out: :py:class:`dolfin.Vector` + + ..note:: :code:`setLinearizationPoint` must be called before calling this method. """ out.zero() @@ -272,14 +363,22 @@ def apply_ij(self,i,j, dir, out): def apply_ijk(self,i,j,k,dir1,dir2, out): """ - Apply the third order variation of the q.o.i. w.r.t. ijk in direction dir1, dir2 for j and k - :param i: STATE or PARAMETER or CONTROL - :param j: - :param k: - :param dir1: - :param dir2: - :param out: - :return: out + Apply the third order variation of the QoI in the \ + :code:`i` th, :code:`j` th, and :code:`k` th variables in directions \ + :code:`dir1` and :code:`dir2`. + + :param i: First variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type i: int + :param j: Second variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type j: int + :param k: Third variable index (:code:`soupy.STATE`, :code:`soupy.PARAMETER`, or :code:`soupy.CONTROL`) + :type k: int + :param dir1: Direction for variable :code:`j` + :type dir1: :py:class:`dolfin.Vector` + :param dir2: Direction for variable :code:`k` + :type dir2: :py:class:`dolfin.Vector` + :param out: The assembled third variation + :type out: :py:class:`dolfin.Vector` """ out.zero() @@ -287,8 +386,10 @@ def apply_ijk(self,i,j,k,dir1,dir2, out): def setLinearizationPoint(self, x, gauss_newton_approx=False): """ Specify the linearization point for computation of the second variations in method apply_ij. - INPUTS: - - x = [u,m,p,z] is a list of the state u, parameter m, and adjoint variable p + + :param x: List of vectors :code:`[u, m, p, z]` representing the state, + parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` """ for i in range(len(x)): self.x[i].zero() diff --git a/soupy/modeling/meanVarRiskMeasure.py b/soupy/modeling/meanVarRiskMeasure.py index 44b93d6..f717566 100644 --- a/soupy/modeling/meanVarRiskMeasure.py +++ b/soupy/modeling/meanVarRiskMeasure.py @@ -23,130 +23,146 @@ def meanVarRiskMeasureSettings(data = {}): - # This should be a Parameter - # data['nsamples'] = [100,'Number of Monte Carlo samples'] - data['beta'] = [0,'Weighting factor for variance'] + data['beta'] = [0,'Weighting factor for variance'] - return ParameterList(data) + return ParameterList(data) class MeanVarRiskMeasure(RiskMeasure): - """ - Class for memory efficient evaluation of the Mean + Variance risk measure - E[X] + beta Var[X]. - """ - - def __init__(self, control_model, prior, settings = meanVarRiskMeasureSettings()): - """ - Parameters - - :code: `control_model` control model of problem - - :code: `prior` prior for uncertain parameter - - :code: `settings` additional settings - """ - self.model = control_model - self.prior = prior - self.settings = settings - self.settings.showMe() - # self.n_samples = self.settings['nsamples'] - self.beta = settings['beta'] - - - # Aggregate components for computing cost, grad, hess - self.x = self.model.generate_vector() - self.g = self.model.generate_vector(CONTROL) - self.q_samples = np.zeros(1) - - self.uhat = self.model.generate_vector(STATE) - self.phat = self.model.generate_vector(STATE) - self.zhelp = self.model.generate_vector(CONTROL) - self.rhs_fwd = self.model.generate_vector(STATE) - self.rhs_adj = self.model.generate_vector(STATE) - self.rhs_adj2 = self.model.generate_vector(STATE) - - self.q_bar = 0 - self.g_bar = self.model.generate_vector(CONTROL) - self.qg_bar = self.model.generate_vector(CONTROL) - - # For sampling - self.noise = dl.Vector() - self.prior.init_vector(self.noise, "noise") - - def generate_vector(self, component = "ALL"): - return self.model.generate_vector(component) - - def computeComponents(self, z, order=0, sample_size=100, rng=None): - """ - Computes the components for the stochastic approximation of the cost - Parameters: - - :code: `z` the control variable - - :code: `order` the order of derivatives needed. - 0 for cost. 1 for grad. 2 for Hessian - - :code: `sample_size` number of samples for sample average - - :code: `rng` rng for the sampling (optional) - """ - self.q_samples = np.zeros(sample_size) - if order >= 1: - self.g_bar.zero() - self.qg_bar.zero() - - for i in range(sample_size): - # Assign control - self.x[CONTROL].zero() - self.x[CONTROL].axpy(1.0, z) - - # Sample parameter - if rng is None: - parRandom.normal(1.0, self.noise) - else: - rng.normal(1.0, self.noise) - self.prior.sample(self.noise, self.x[PARAMETER]) - - # Solve state - self.model.solveFwd(self.x[STATE], self.x) - self.q_samples[i] = self.model.cost(self.x) - - if order >= 1: - self.model.solveAdj(self.x[ADJOINT], self.x) - self.model.evalGradientControl(self.x, self.g) - self.g_bar.axpy(1/sample_size, self.g) - self.qg_bar.axpy(self.q_samples[i]/sample_size, self.g) - - # Still need Hessian code - # if i % 10 == 0: - # print(i) - - self.q_bar = np.mean(self.q_samples) - - - def cost(self): - """ - Evaluates the cost given by the risk measure - Assumes :code: `computeComponents` has been called - """ - return self.q_bar + self.beta * np.std(self.q_samples)**2 - - def costGrad(self, g): - """ - Evaluates the gradient by the risk measure - Assumes :code: `computeComponents` has been called with :code: `order>=1` - Parameters - - :code: `g` output vector for the gradient - """ - g.zero() - g.axpy(1.0, self.g_bar) - g.axpy(2*self.beta, self.qg_bar) - g.axpy(-2*self.beta*self.q_bar, self.g_bar) - - def costHessian(self, zhat, Hzhat): - self.model.setPointForHessianEvaluations(self.x) - self.model.applyCz(zhat, self.rhs_fwd) - self.model.solveFwdIncremental(self.uhat, self.rhs_fwd) - self.model.applyWuu(self.uhat, self.rhs_adj) - self.model.applyWuz(zhat, self.rhs_adj2) - self.rhs_adj.axpy(-1., self.rhs_adj2) - self.model.solveAdjIncremental(self.phat, self.rhs_adj) - self.model.applyWzz(zhat, Hzhat) - - self.model.applyCzt(self.phat, self.zhelp) - Hzhat.axpy(1., self.zhelp) - self.model.applyWzu(self.uhat, self.zhelp) - Hzhat.axpy(-1., self.zhelp) \ No newline at end of file + """ + Class for memory efficient evaluation of the Mean + Variance risk measure + + .. math:: \\rho[Q](z) = \mathbb{E}_m[Q(m,z)] + \\beta \mathbb{V}_m[Q(m,z)] + + that allows for stochastic approximations and SGD + """ + + def __init__(self, control_model, prior, settings = meanVarRiskMeasureSettings()): + """ + Constructor: + + :param control_model: control model containing the :code:`soupy.PDEVariationalControlProblem` + and :code:`soupy.ControlQoI` + :type control_model: :py:class:`soupy.ControlModel` + :param prior: The prior distribution for the random parameter + :type prior: :py:class:`hippylib.Prior` + :param settings: additional settings + :type settings: :py:class:`hippylib.ParameterList` + """ + self.model = control_model + self.prior = prior + self.settings = settings + self.settings.showMe() + # self.n_samples = self.settings['nsamples'] + self.beta = settings['beta'] + + + # Aggregate components for computing cost, grad, hess + self.x = self.model.generate_vector() + self.g = self.model.generate_vector(CONTROL) + self.q_samples = np.zeros(1) + + self.uhat = self.model.generate_vector(STATE) + self.phat = self.model.generate_vector(STATE) + self.zhelp = self.model.generate_vector(CONTROL) + self.rhs_fwd = self.model.generate_vector(STATE) + self.rhs_adj = self.model.generate_vector(STATE) + self.rhs_adj2 = self.model.generate_vector(STATE) + + self.q_bar = 0 + self.g_bar = self.model.generate_vector(CONTROL) + self.qg_bar = self.model.generate_vector(CONTROL) + + # For sampling + self.noise = dl.Vector() + self.prior.init_vector(self.noise, "noise") + + def generate_vector(self, component = "ALL"): + return self.model.generate_vector(component) + + def computeComponents(self, z, order=0, sample_size=100, rng=None): + """ + Computes the components for the evaluation of the risk measure + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed. + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + :param rng: The random number generator used for sampling. If :code:`None` the default + uses :py:class:`hippylib.parRandom` + :type rng: :py:class:`hippylib.Random` + """ + + # Check if a new control variable is used + self.q_samples = np.zeros(sample_size) + if order >= 1: + self.g_bar.zero() + self.qg_bar.zero() + + for i in range(sample_size): + # Assign control + self.x[CONTROL].zero() + self.x[CONTROL].axpy(1.0, z) + + # Sample parameter + if rng is None: + parRandom.normal(1.0, self.noise) + else: + rng.normal(1.0, self.noise) + self.prior.sample(self.noise, self.x[PARAMETER]) + + # Solve state + self.model.solveFwd(self.x[STATE], self.x) + self.q_samples[i] = self.model.cost(self.x) + + if order >= 1: + self.model.solveAdj(self.x[ADJOINT], self.x) + self.model.evalGradientControl(self.x, self.g) + self.g_bar.axpy(1/sample_size, self.g) + self.qg_bar.axpy(self.q_samples[i]/sample_size, self.g) + + # Still need Hessian code + # if i % 10 == 0: + # print(i) + + self.q_bar = np.mean(self.q_samples) + + + def cost(self): + """ + Evaluates the value of the risk measure once components have been computed + + :return: Value of the cost functional + + .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` + """ + return self.q_bar + self.beta * np.std(self.q_samples)**2 + + def costGrad(self, g): + """ + Evaluates the value of the risk measure once components have been computed + + :param g: (Dual of) the gradient of the risk measure to store result in + :type g: :py:class:`dolfin.Vector` + + .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` + """ + g.zero() + g.axpy(1.0, self.g_bar) + g.axpy(2*self.beta, self.qg_bar) + g.axpy(-2*self.beta*self.q_bar, self.g_bar) + + def costHessian(self, zhat, Hzhat): + self.model.setPointForHessianEvaluations(self.x) + self.model.applyCz(zhat, self.rhs_fwd) + self.model.solveFwdIncremental(self.uhat, self.rhs_fwd) + self.model.applyWuu(self.uhat, self.rhs_adj) + self.model.applyWuz(zhat, self.rhs_adj2) + self.rhs_adj.axpy(-1., self.rhs_adj2) + self.model.solveAdjIncremental(self.phat, self.rhs_adj) + self.model.applyWzz(zhat, Hzhat) + + self.model.applyCzt(self.phat, self.zhelp) + Hzhat.axpy(1., self.zhelp) + self.model.applyWzu(self.uhat, self.zhelp) + Hzhat.axpy(-1., self.zhelp) \ No newline at end of file diff --git a/soupy/modeling/meanVarRiskMeasureSAA.py b/soupy/modeling/meanVarRiskMeasureSAA.py index de8164f..28b2be5 100644 --- a/soupy/modeling/meanVarRiskMeasureSAA.py +++ b/soupy/modeling/meanVarRiskMeasureSAA.py @@ -26,7 +26,6 @@ def meanVarRiskMeasureSAASettings(data = {}): - # This should be a Parameter data['sample_size'] = [100,'Number of Monte Carlo samples'] data['beta'] = [0,'Weighting factor for variance'] data['seed'] = [1, 'rng seed for sampling'] @@ -52,17 +51,27 @@ def _allocate_sample_sizes(sample_size, comm_sampler): class MeanVarRiskMeasureSAA_MPI(RiskMeasure): """ Class for the mean + variance risk measure using sample average approximation - E[X] + beta Var[X] - with optional sample parallelism using MPI + + .. math:: \\rho[Q](z) = \mathbb{E}_m[Q(m,z)] + \\beta \mathbb{V}_m[Q(m,z)] + + with sample parallelism using MPI + + .. note:: currently does not support simultaneous sample and mesh partition parallelism """ def __init__(self, control_model, prior, settings = meanVarRiskMeasureSAASettings(), comm_sampler=mpi4py.MPI.COMM_WORLD): """ - Parameters - - :code: `control_model` control model of problem - - :code: `prior` prior for uncertain parameter - - :code: `settings` additional settings - - :code: `comm_sampler` MPI communicator for the sampling parallelism + Constructor: + + :param control_model: control model containing the :code:`soupy.PDEVariationalControlProblem` + and :code:`soupy.ControlQoI` + :type control_model: :py:class:`soupy.ControlModel` + :param prior: The prior distribution for the random parameter + :type prior: :py:class:`hippylib.Prior` + :param settings: additional settings + :type settings: :py:class:`hippylib.ParameterList` + :param comm_sampler: MPI communicator for sample parallelism + :type comm_sampler: :py:class:`mpi4py.MPI.Comm` """ self.model = control_model self.prior = prior @@ -134,12 +143,14 @@ def generate_vector(self, component = "ALL"): def computeComponents(self, z, order=0, **kwargs): """ - Computes the components for the stochastic approximation of the cost - Parameters: - - :code: `z` the control variable - - :code: `order` the order of derivatives needed. - 0 for cost. 1 for grad. 2 for Hessian - - :code: `**kwargs` dummy keyword arguments for compatibility + Computes the components for the evaluation of the risk measure + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed. + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + :param kwargs: dummy keyword arguments for compatibility """ # Check if a new control variable is used new_forward_solve = False @@ -219,17 +230,22 @@ def computeComponents(self, z, order=0, **kwargs): def cost(self): """ - Evaluates the cost given by the risk measure - Assumes :code: `computeComponents` has been called + Evaluates the value of the risk measure once components have been computed + + :return: Value of the cost functional + + .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` """ return self.q_bar + self.beta * (self.q2_bar - self.q_bar**2) def costGrad(self, g): """ - Evaluates the gradient by the risk measure - Assumes :code: `computeComponents` has been called with :code: `order>=1` - Parameters - - :code: `g` output vector for the gradient + Evaluates the value of the risk measure once components have been computed + + :param g: (Dual of) the gradient of the risk measure to store result in + :type g: :py:class:`dolfin.Vector` + + .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` """ # print("(proc %d) q_bar = %g" %(self.comm_sampler.Get_rank(), self.q_bar)) g.zero() @@ -256,6 +272,15 @@ def costHessian(self, zhat, Hzhat): # Hzhat.axpy(-1., self.zhelp) def gatherSamples(self, root=0): + """ + Gather the QoI samples on the root process + + :param root: Rank of the process to gather the samples on + :type root: int + + :return: An array of the sample QoI values + :return type: :py:class:`numpy.ndarray` + """ q_all = None if self.comm_sampler.Get_rank() == 0: q_all = np.zeros(self.sample_size) @@ -265,16 +290,22 @@ def gatherSamples(self, root=0): class MeanVarRiskMeasureSAA(RiskMeasure): """ - Class for memory efficient evaluation of the Mean + Variance risk measure - E[X] + beta Var[X]. + Class for the mean + variance risk measure using sample average approximation + + .. math:: \\rho[Q](z) = \mathbb{E}_m[Q(m,z)] + \\beta \mathbb{V}_m[Q(m,z)] """ def __init__(self, control_model, prior, settings = meanVarRiskMeasureSAASettings()): """ - Parameters - - :code: `control_model` control model of problem - - :code: `prior` prior for uncertain parameter - - :code: `settings` additional settings + Constructor: + + :param control_model: control model containing the :code:`soupy.PDEVariationalControlProblem` + and :code:`soupy.ControlQoI` + :type control_model: :py:class:`soupy.ControlModel` + :param prior: The prior distribution for the random parameter + :type prior: :py:class:`hippylib.Prior` + :param settings: additional settings + :type settings: :py:class:`hippylib.ParameterList` """ self.model = control_model self.prior = prior @@ -318,13 +349,17 @@ def generate_vector(self, component = "ALL"): def computeComponents(self, z, order=0, **kwargs): """ - Computes the components for the stochastic approximation of the cost - Parameters: - - :code: `z` the control variable - - :code: `order` the order of derivatives needed. - 0 for cost. 1 for grad. 2 for Hessian - - :code: `**kwargs` dummy keyword arguments for compatibility + + Computes the components for the evaluation of the risk measure + + :param z: the control variable + :type z: :py:class:`dolfin.Vector` + :param order: Order of the derivatives needed. + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + :param kwargs: dummy keyword arguments for compatibility """ + # Check if a new control variable is used new_forward_solve = False self.diff_helper.zero() @@ -387,17 +422,22 @@ def computeComponents(self, z, order=0, **kwargs): def cost(self): """ - Evaluates the cost given by the risk measure - Assumes :code: `computeComponents` has been called + Computes the value of the risk measure once components have been computed + + :return: Value of the cost functional + + .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` """ return self.q_bar + self.beta * np.std(self.q_samples)**2 def costGrad(self, g): """ - Evaluates the gradient by the risk measure - Assumes :code: `computeComponents` has been called with :code: `order>=1` - Parameters - - :code: `g` output vector for the gradient + Computes the value of the risk measure once components have been computed + + :param g: (Dual of) the gradient of the risk measure to store the result in + :type g: :py:class:`dolfin.Vector` + + .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` """ g.zero() g.axpy(1.0, self.g_bar) diff --git a/soupy/modeling/penalization.py b/soupy/modeling/penalization.py index 7bc7fa7..879e020 100644 --- a/soupy/modeling/penalization.py +++ b/soupy/modeling/penalization.py @@ -33,13 +33,20 @@ def hessian(self, z, zhat, out): class MultiPenalization(Penalization): """ Class for a sum of penalization terms + + .. math:: P(z) = \sum_{i=1}^{n} \\alpha_i P_i(z) + """ def __init__(self, Vh, penalization_list, alpha_list=None): """ - - :code: `Vh` function space for STATE, PARAMETER, ADJOINT, CONTROL - - :code: `penalization_list` a list of Penalization objects - - :code: `alpha_list` optional list of weights, assumed to all be 1 - if None is given + Constructor + + :param Vh: List of function spaces the state, parameter, adjoint, and control + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param penalization_list: List of penalization objects + :type penalization_list: list of :py:class:`Penalization` + :param alpha_list: List of weights for each penalization term + :type alpha_list: list of floats """ self.Vh = Vh self.helper = dl.Function(Vh[CONTROL]).vector() @@ -53,18 +60,42 @@ def __init__(self, Vh, penalization_list, alpha_list=None): def cost(self, z): + """ + Compute the penalization at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ cost = 0 for alpha, penalization in zip(self.alpha_list, self.penalization_list): cost += alpha * penalization.cost(z) return cost def grad(self, z, out): + """ + Compute the gradient of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param out: The assembled gradient vector + :type out: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ out.zero() for alpha, penalization in zip(self.alpha_list, self.penalization_list): penalization.grad(z, self.helper) out.axpy(alpha, self.helper) def hessian(self, z, zhat, out): + """ + Compute the Hessian of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param zhat: The direction for Hessian action + :type zhat: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param out: The assembled Hessian action vector + :type out: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ out.zero() for alpha, penalization in zip(self.alpha_list, self.penalization_list): penalization.hessian(z, zhat, self.helper) @@ -73,17 +104,23 @@ def hessian(self, z, zhat, out): class L2Penalization(Penalization): """ - L2 integral over the domain - P(z) = \alpha \int_{\Omega) |z|^2 dx + :math:`L^2(\Omega)` integral over the domain - For finite dimensional controls `z`, this amounts to a little \ell_2 norm - In this case, `Vh[soupy.CONTROL]` needs to be a - `dolfin.VectorFunctionSpace` of reals + .. math:: P(z) = \\alpha \int_{\Omega} |z|^2 dx + + For finite dimensional controls `z`, this amounts to a little \ + :math:`\ell_2` norm \ + In this case, :code:`Vh[soupy.CONTROL]` needs to be a \ + :code:`dolfin.VectorFunctionSpace` of reals """ def __init__(self, Vh, alpha): """ - - :code: `Vh` function space for STATE, PARAMETER, ADJOINT, CONTROL - - :code: `alpha` weighting factor + Constructor + + :param Vh: List of function spaces the state, parameter, adjoint, and control + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param alpha: Weight for the penalization term + :type alpha: float """ self.Vh = Vh self.alpha = alpha @@ -103,6 +140,12 @@ def init_vector(self, v, dim=0): def cost(self, z): + """ + Compute the penalization at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`AugmentedVector` + """ if isinstance(z, AugmentedVector): z_sub = z.get_vector() self.M.mult(z_sub, self.Mz) @@ -112,6 +155,14 @@ def cost(self, z): return self.alpha * z.inner(self.Mz) def grad(self, z, out): + """ + Compute the gradient of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`AugmentedVector` + :param out: The assembled gradient vector + :type out: :py:class:`dolfin.Vector` or :py:class:`AugmentedVector` + """ out.zero() if isinstance(z, AugmentedVector): out_sub = out.get_vector() @@ -123,6 +174,16 @@ def grad(self, z, out): out.axpy(2.0*self.alpha, self.Mz) def hessian(self, z, zhat, out): + """ + Compute the Hessian of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param zhat: The direction for Hessian action + :type zhat: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param out: The assembled Hessian action vector + :type out: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ out.zero() if isinstance(z, AugmentedVector): out_sub = out.get_vector() @@ -138,15 +199,21 @@ def hessian(self, z, zhat, out): class WeightedL2Penalization(Penalization): """ A weighted L2 norm penalization - P(z) = z^T M z - where M is some symmetric positive definite weight matrix + .. math:: P(z) = z^T M z + + where :math:`M` is a symmetric positive definite weight matrix """ def __init__(self, Vh, M, alpha): """ - - :code: `Vh` function space for STATE, PARAMETER, ADJOINT, CONTROL - - :code: `M` weighting matrix with method `mult` - - :code: `alpha` weighting factor - """ + Constructor: + + :param Vh: List of function spaces the state, parameter, adjoint, and control + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param M: Weighting matrix with method :code:`mult` + :type M: :py:class:`dolfin.Matrix` like + :param alpha: Weight for the penalization term + :type alpha: float + """ self.Vh = Vh self.M = M self.alpha = alpha @@ -160,6 +227,12 @@ def init_vector(self, v, dim=0): def cost(self, z): + """ + Compute the penalization at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ if isinstance(z, AugmentedVector): z_sub = z.get_vector() self.M.mult(z_sub, self.Mz) @@ -169,6 +242,14 @@ def cost(self, z): return self.alpha * z.inner(self.Mz) def grad(self, z, out): + """ + Compute the gradient of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param out: The assembled gradient vector + :type out: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ out.zero() if isinstance(z, AugmentedVector): out_sub = out.get_vector() @@ -180,6 +261,16 @@ def grad(self, z, out): out.axpy(2.0*self.alpha, self.Mz) def hessian(self, z, zhat, out): + """ + Compute the Hessian of the penalty at given control :math:`z` + + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param zhat: The direction for Hessian action + :type zhat: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + :param out: The assembled Hessian action vector + :type out: :py:class:`dolfin.Vector` or :py:class:`soupy.AugmentedVector` + """ out.zero() if isinstance(z, AugmentedVector): out_sub = out.get_vector() diff --git a/soupy/modeling/smoothPlusApproximation.py b/soupy/modeling/smoothPlusApproximation.py index 18c2d82..9a3c573 100644 --- a/soupy/modeling/smoothPlusApproximation.py +++ b/soupy/modeling/smoothPlusApproximation.py @@ -14,6 +14,10 @@ import numpy as np class SmoothPlusApproximationSoftplus: + """ + Implements the smooth approximation to the maximum function \ + :math:`\\max(0, t)` using the softplus function + """ def __init__(self, epsilon=1e-2): self.epsilon = epsilon @@ -28,6 +32,11 @@ def hessian(self, t): class SmoothPlusApproximationQuartic: + """ + Implements the smooth approximation to the maximum function \ + :math:`\\max(0, t)` a piecewise quartic function + """ + def __init__(self, epsilon=1e-2): self.epsilon = epsilon diff --git a/soupy/modeling/superquantileRiskMeasureSAA.py b/soupy/modeling/superquantileRiskMeasureSAA.py index 36843ea..c860c9e 100644 --- a/soupy/modeling/superquantileRiskMeasureSAA.py +++ b/soupy/modeling/superquantileRiskMeasureSAA.py @@ -38,7 +38,7 @@ def sampleSuperquantile(samples, beta): def sampleSuperquantileByMinimization(samples, beta, epsilon=1e-2): """ - Evaluate superquantile from samples + Evaluate superquantile from samples by minimization """ quantile = np.percentile(samples, beta * 100) smoothPlus = SmoothPlusApproximationQuartic(epsilon=epsilon) @@ -50,7 +50,6 @@ def sampleSuperquantileByMinimization(samples, beta, epsilon=1e-2): def superquantileRiskMeasureSAASettings(data = {}): - # This should be a Parameter data['sample_size'] = [100,'Number of Monte Carlo samples'] data['beta'] = [0.95, 'Quantile value for superquantile'] data['epsilon'] = [0.01, 'Sharpness of smooth plus approximation'] @@ -77,16 +76,22 @@ def _allocate_sample_sizes(sample_size, comm_sampler): class SuperquantileRiskMeasureSAA_MPI(RiskMeasure): """ Class for the a sample average approximation of the superquantile risk measure (CVaR) - with optional sample parallelism with MPI + with sample parallelism using MPI """ def __init__(self, control_model, prior, settings = superquantileRiskMeasureSAASettings(), comm_sampler=mpi4py.MPI.COMM_WORLD): """ - Parameters - - :code: `control_model` control model of problem - - :code: `prior` prior for uncertain parameter - - :code: `settings` additional settings - - :code: `comm_sampler` MPI communicator for the sampling parallelism + Constructor: + + :param control_model: control model containing the :code:`soupy.PDEVariationalControlProblem` + and :code:`soupy.ControlQoI` + :type control_model: :py:class:`soupy.ControlModel` + :param prior: The prior distribution for the random parameter + :type prior: :py:class:`hippylib.Prior` + :param settings: additional settings given in :code:`superquantileRiskMeasureSAASettings()` + :type settings: :py:class:`hippylib.ParameterList` + :param comm_sampler: MPI communicator for sample parallelism + :type comm_sampler: :py:class:`mpi4py.MPI.Comm` """ self.model = control_model self.prior = prior @@ -164,12 +169,14 @@ def generate_vector(self, component = "ALL"): def computeComponents(self, zt, order=0, **kwargs): """ - Computes the components for the stochastic approximation of the cost - Parameters: - - :code: `z` the control variable - - :code: `order` the order of derivatives needed. - 0 for cost. 1 for grad. 2 for Hessian - - :code: `**kwargs` dummy keyword arguments for compatibility + Computes the components for the evaluation of the risk measure + + :param zt: the control variable with the scalar :code:`t` appended + :type zt: :py:class:`soupy.AugmentedVector` + :param order: Order of the derivatives needed. + :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian + :type order: int + :param kwargs: dummy keyword arguments for compatibility """ z = zt.get_vector() t = zt.get_scalar() @@ -251,18 +258,23 @@ def computeComponents(self, zt, order=0, **kwargs): def cost(self): """ - Evaluates the cost given by the risk measure - Assumes :code: `computeComponents` has been called + Evaluates the value of the risk measure once components have been computed + + :return: Value of the cost functional + + .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` """ t = self.zt.get_scalar() return t + 1/(1-self.beta) * self.s_bar def costGrad(self, gt): """ - Evaluates the gradient by the risk measure - Assumes :code: `computeComponents` has been called with :code: `order>=1` - Parameters - - :code: `g` output vector as an augmented vector + Evaluates the value of the risk measure once components have been computed + + :param g: (Dual of) the gradient of the risk measure to store result in + :type g: :py:class:`dolfin.Vector` + + .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order>=1` """ # print("(proc %d) q_bar = %g" %(self.comm_sampler.Get_rank(), self.q_bar)) dzJ_np = self.sprime_g_bar.get_local()/(1-self.beta) @@ -275,6 +287,15 @@ def costHessian(self, zhat, Hzhat): return def gatherSamples(self, root=0): + """ + Gather the QoI samples on the root process + + :param root: Rank of the process to gather the samples on + :type root: int + + :return: An array of the sample QoI values + :return type: :py:class:`numpy.ndarray` + """ q_all = None if self.comm_sampler.Get_rank() == 0: q_all = np.zeros(self.sample_size) @@ -284,6 +305,10 @@ def gatherSamples(self, root=0): def superquantile(self): """ Evaluate the superquantile using the computed samples + + :return: Value of the superquantile by sampling + + .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` """ q_all = self.gatherSamples(root=0) value = 0.0 diff --git a/soupy/modeling/variables.py b/soupy/modeling/variables.py index adc5f1c..5ec4ef9 100644 --- a/soupy/modeling/variables.py +++ b/soupy/modeling/variables.py @@ -13,7 +13,10 @@ """ Enumerator for the variables of the inverse problem: -- the STATE, PARAMETER, and ADJOINT variables. + - :code:`STATE = 0` for the state variable + - :code:`PARAMETER = 1` for the parameter variable + - :code:`ADJOINT = 2` for the adjoint variable + - :code:`CONTROL = 3` for the control variable """ STATE= 0 PARAMETER = 1 diff --git a/soupy/optimization/bfgs.py b/soupy/optimization/bfgs.py index 43fc114..6901ce6 100644 --- a/soupy/optimization/bfgs.py +++ b/soupy/optimization/bfgs.py @@ -54,6 +54,15 @@ def init_vector(self, x, dim): raise def solve(self, x, b): + """ + Applies the operator :math:`d0 I` + + :param x: solution vector + :type x: :py:class:`dolfin.Vector` + :param b: right-hand side vector + :type b: :py:class:`dolfin.Vector` + """ + # print("\t\t WITHIN IDENTITY") # print("\t\t", self.d0) # print("\t\t", x, x.get_local()) @@ -79,24 +88,24 @@ def __init__(self, parameters=BFGSoperator_ParameterList()): def set_H0inv(self, H0inv): """ Set user-defined operator corresponding to :code:`H0inv` - Input: - :code:`H0inv`: Fenics operator with method :code:`solve()` + + :param H0inv: Fenics operator with method :code:`solve()` """ self.H0inv = H0inv def solve(self, x, b): """ - Solve system: :math:`H_{bfgs} x = b` - where :math:`H_{bfgs}` is the approximation to the Hessian build by BFGS. - That is, we apply - .. math:: - x = (H_{bfgs})^{-1} b = H_k b + Solve system: + :math:`H_{\mathrm{bfgs}} x = b` + where :math:`H_{\mathrm{bfgs}}` is the approximation to the Hessian build by BFGS. + That is, we apply :math:`x = (H_{\mathrm{bfgs}})^{-1} b = H_k b` where :math:`H_k` matrix is BFGS approximation to the inverse of the Hessian. Computation done via double-loop algorithm. - Inputs: - :code:`x = dolfin.Vector` - `[out]` - :code:`b = dolfin.Vector` - `[in]` + :param x: The solution to the system + :type x: :py:class:`dolfin.Vector` + :param b: The right-hand side of the system + :type b: :py:class:`dolfin.Vector` """ A = [] if self.help is None: @@ -126,12 +135,12 @@ def solve(self, x, b): def update(self, s, y): """ Update BFGS operator with most recent gradient update. - To handle potential break from secant condition, update done via damping - - Inputs: - :code:`s = dolfin.Vector` `[in]` - corresponds to update in medium parameters. - :code:`y = dolfin.Vector` `[in]` - corresponds to update in gradient. + + :param s: The update in medium parameters + :type s: :py:class:`dolfin.Vector` + :param y: The update in gradient + :type y: :py:class:`dolfin.Vector` """ damp = self.parameters["BFGS_damping"] memlim = self.parameters["memory_limit"] @@ -197,6 +206,9 @@ def __init__(self, cost_functional, parameters=BFGS_ParameterList()): """ Initialize the BFGS solver. Type :code:`BFGS_ParameterList().showMe()` for default parameters and their description + + :param cost_functional: The cost functional + :param parameters: The parameters for the BFGS solver """ self.cost_functional = cost_functional @@ -213,12 +225,20 @@ def __init__(self, cost_functional, parameters=BFGS_ParameterList()): def solve(self, z, H0inv=RescaledIdentity(), box_bounds=None, constraint_projection=None): """ Solve the constrained optimization problem with initial guess :code:`z` - :code:`z` will be overwritten. - :code:`H0inv`: the initial approximated inverse of the Hessian for the BFGS operator. It has an - optional method :code:`update(x)` that will update the operator - :code:`box_bounds`: Bound constraint (list with two entries: min and max). Can be either a scalar value or a - :code:`constraint_projection` Alternative projectable constraint - Return the solution :code:`z` and :code:`results`, a dictionary of iterates + + :param z: The initial guess + :type z: :py:class:`dolfin.Vector` + :param H0inv: Initial approximation of the inverse of the Hessian. + Has optional method :code:`update(x)` that will update the operator + :param box_bounds: Bound constraint. A list with two entries (min and max). + Can be either a scalar value or a :code:`dolfin.Vector` of the same size as :code:`z` + :type box_bounds: list + :param constraint_projection: Alternative projectable constraint + :type constraint_projection: :py:class:`ProjectableConstraint` + + :return: The optimization solution :code:`z` and a dictionary of information + + .. note:: The input :code:`z` will be overwritten """ if box_bounds is not None: diff --git a/soupy/optimization/projectableConstraint.py b/soupy/optimization/projectableConstraint.py index 11a99c7..c7f7a5a 100644 --- a/soupy/optimization/projectableConstraint.py +++ b/soupy/optimization/projectableConstraint.py @@ -12,21 +12,39 @@ # Software Foundation) version 3.0 dated June 1991. class ProjectableConstraint: + """ + Base class for a constraint for which a projection operator \ + into the feasible set is available + """ def project(self, z): + """ + Projects the vector :code:`z` onto the constraint set + + :param z: + :type z: :py:class:`dolfin.Vector` + """ raise NotImplementedError("Child class should implement project") def cost(self, z): + """ + Returns the amount of violation of the constraint by :code:`z` + + :param z: + :type z: :py:class:`dolfin.Vector` + """ raise NotImplementedError("Child class should implement cost") class InnerProductEqualityConstraint(ProjectableConstraint): """ - Class implements the constraint c^T x - a = 0 + Class implements the constraint :math:`c^T z - a = 0` """ def __init__(self, c, a): """ - - :code: `c` is constraint `dolfin.Vector()` - - :code: `a` value of the constraint + :param c: Constraint vector + :type c: :py:class:`dolfin.Vector` + :param a: Value of the constraint + :type a: float """ self.c = c self.a = a @@ -34,11 +52,22 @@ def __init__(self, c, a): def project(self, z): + """ + Projects the vector :code:`z` onto the hyperplane \ + :math:`\{z : c^T z = a\}` + + :param z: + :type z: :py:class:`dolfin.Vector` + """ factor = (self.c.inner(z) - self.a)/self.c_norm2 z.axpy(-factor, self.c) def cost(self, z): """ - returns the amount of violation of the constraint + Returns constraint violation math:`c^T z - a` for input :code:`z` + :param z: + :type z: :py:class:`dolfin.Vector` + + :return: :math:`c^T z - a` for input :code:`z` """ return self.c.inner(z) - self.a diff --git a/soupy/optimization/sgd.py b/soupy/optimization/sgd.py index c7fd26a..6aa32c9 100644 --- a/soupy/optimization/sgd.py +++ b/soupy/optimization/sgd.py @@ -37,22 +37,18 @@ def SGD_ParameterList(): class SGD: """ - Stochastic Gradient Descent to solve optimization under uncertainty problems - Globalization is performed using the Armijo sufficient reduction condition (backtracking). - The stopping criterion is based on a control on the norm of the gradient. - - The user must provide a model that describes the forward problem, cost functionals, and all the - derivatives for the gradient. - - More specifically the model object should implement following methods: - - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint. - - :code:`cost(x)` -> evaluate the cost functional, report regularization part and misfit separately. - - :code:`solveFwd(out, x)` -> solve the possibly non linear forward problem. - - :code:`solveAdj(out, x)` -> solve the linear adjoint problem. - - :code:`evalGradientParameter(x, out)` -> evaluate the gradient of the parameter and compute its norm. - - :code:`Rsolver()` --> A solver for the regularization term. - - Type :code:`help(Model)` for additional information. + Stochastic Gradient Descent to solve optimization under uncertainty problems \ + The stopping criterion is based on a control on the norm of the gradient.\ + Line search is only possible when the option :code:`stochastic_approximation` \ + is False. + + + The user must provide a cost functional that provides the evaluation and gradient + + More specifically the cost functional object should implement following methods: + - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint, and control. + - :code:`cost(z, order, rng)` -> evaluate the cost functional which allows a given :code:`rng` + - :code:`costGrad(z, g)` -> evaluate the gradient of the cost functional """ termination_reasons = [ "Maximum number of Iteration reached", #0 @@ -63,8 +59,14 @@ class SGD: def __init__(self, cost_functional, parameters = SGD_ParameterList()): """ - Initialize the Stochastic Gradient Descent solver. Type :code:`SGD_ParameterList().showMe()` for list of default parameters - and their descriptions. + Constructor for the SGD solver + + :param cost_functional: The cost functional object + :type cost_functional: :py:class:`soupy.ControlCostFunctional` or similar + :param parameters: The parameters of the solver. + Type :code:`SGD_ParameterList().showMe()` for list of default parameters + and their descriptions. + :type parameters: :py:class:`hippylib.ParameterList`. """ self.cost_functional = cost_functional self.parameters = parameters @@ -78,9 +80,19 @@ def __init__(self, cost_functional, parameters = SGD_ParameterList()): def solve(self, z, box_bounds=None, constraint_projection=None): """ - Solve the constrained optimization problem with initial guess :code:`x = [u,a,p]`. Return the solution :code:`[u,a,p]`. - .. note:: :code:`x` will be overwritten. + Solve the constrained optimization problem using steepest descent with initial guess :code:`z` + + :param z: The initial guess + :type z: :py:class:`dolfin.Vector` + :param box_bounds: Bound constraint. A list with two entries (min and max). + Can be either a scalar value or a :code:`dolfin.Vector` of the same size as :code:`z` + :type box_bounds: list + :param constraint_projection: Alternative projectable constraint + :type constraint_projection: :py:class:`ProjectableConstraint` + + :return: The optimization solution :code:`z` and summary of iterates + .. note:: The input :code:`z` is overwritten """ rel_tol = self.parameters["rel_tolerance"] abs_tol = self.parameters["abs_tolerance"] diff --git a/soupy/optimization/steepestDescent.py b/soupy/optimization/steepestDescent.py index 398a140..d39a13d 100644 --- a/soupy/optimization/steepestDescent.py +++ b/soupy/optimization/steepestDescent.py @@ -41,18 +41,12 @@ class SteepestDescent: Globalization is performed using the Armijo sufficient reduction condition (backtracking). The stopping criterion is based on a control on the norm of the gradient. - The user must provide a model that describes the forward problem, cost functionals, and all the - derivatives for the gradient. - - More specifically the model object should implement following methods: - - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint. - - :code:`cost(x)` -> evaluate the cost functional, report regularization part and misfit separately. - - :code:`solveFwd(out, x)` -> solve the possibly non linear forward problem. - - :code:`solveAdj(out, x)` -> solve the linear adjoint problem. - - :code:`evalGradientParameter(x, out)` -> evaluate the gradient of the parameter and compute its norm. - - :code:`Rsolver()` --> A solver for the regularization term. - - Type :code:`help(Model)` for additional information. + The user must provide a cost functional that provides the evaluation and gradient + + More specifically the cost functional object should implement following methods: + - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint, and control. + - :code:`cost(z)` -> evaluate the cost functional + - :code:`costGrad(z, g)` -> evaluate the gradient of the cost functional """ termination_reasons = [ "Maximum number of Iteration reached", #0 @@ -63,8 +57,14 @@ class SteepestDescent: def __init__(self, cost_functional, parameters = SteepestDescent_ParameterList()): """ - Initialize the Stochastic Gradient Descent solver. Type :code:`SGD_ParameterList().showMe()` for list of default parameters - and their descriptions. + Constructor for the Steepest Descent solver. + + :param cost_functional: The cost functional object + :type cost_functional: :py:class:`soupy.ControlCostFunctional` or similar + :param parameters: The parameters of the solver. + Type :code:`SteepestDescent_ParameterList().showMe()` for list of default parameters + and their descriptions. + :type parameters: :py:class:`hippylib.ParameterList`. """ self.cost_functional = cost_functional self.parameters = parameters @@ -78,9 +78,19 @@ def __init__(self, cost_functional, parameters = SteepestDescent_ParameterList() def solve(self, z, box_bounds=None, constraint_projection=None): """ - Solve the constrained optimization problem with initial guess :code:`x = [u,a,p]`. Return the solution :code:`[u,a,p]`. - .. note:: :code:`x` will be overwritten. + Solve the constrained optimization problem using steepest descent with initial guess :code:`z` + + :param z: The initial guess + :type z: :py:class:`dolfin.Vector` + :param box_bounds: Bound constraint. A list with two entries (min and max). + Can be either a scalar value or a :code:`dolfin.Vector` of the same size as :code:`z` + :type box_bounds: list + :param constraint_projection: Alternative projectable constraint + :type constraint_projection: :py:class:`ProjectableConstraint` + + :return: The optimization solution :code:`z` and a dictionary of information + .. note:: The input :code:`z` is overwritten """ rel_tol = self.parameters["rel_tolerance"] abs_tol = self.parameters["abs_tolerance"] diff --git a/soupy/utils/penalizationFiniteDifference.py b/soupy/utils/penalizationFiniteDifference.py index 31552f0..92a1eb0 100644 --- a/soupy/utils/penalizationFiniteDifference.py +++ b/soupy/utils/penalizationFiniteDifference.py @@ -22,8 +22,20 @@ def penalizationFiniteDifference(Vh, penalization, z, dz, order=1, delta=1e-4, plotting=False): """ - Finite difference checks the gradient, mainly for the state variable - Also computes the gradients for remaining variables --- most of the time will be zero + Finite difference checks the gradient of the penalization + + :param Vh: List of function spaces for the state, parameter, adjoint, and control variables + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param penalization: Penalization term to check + :type penalization: :py:class:`soupy.Penalization` + :param z: The control variable + :type z: :py:class:`dolfin.Vector` or similar + :param dz: The perturbation to the control variable + :type dz: :py:class:`dolfin.Vector` or similar + :param delta: The finite difference step size + :type delta: float + :plotting: If :code:`true`, plots the finite difference Hessian and analytic Hessian + :type plotting: bool """ z1 = dl.Vector(z) diff --git a/soupy/utils/qoiFiniteDifference.py b/soupy/utils/qoiFiniteDifference.py index 9727602..11bf85d 100644 --- a/soupy/utils/qoiFiniteDifference.py +++ b/soupy/utils/qoiFiniteDifference.py @@ -22,8 +22,23 @@ def qoiFiniteDifference(Vh, qoi, x, du, order=1, delta=1e-4, plotting=False): """ - Finite difference checks the gradient, mainly for the state variable - Also computes the gradients for remaining variables --- most of the time will be zero + Finite difference checks the gradient of the quantity of interest with respect to + the state variable, as well as the parameter and control variables + + :param Vh: List of function spaces for the state, parameter, adjoint, and control variables + :type Vh: list of :py:class:`dolfin.FunctionSpace` + :param qoi: Quantity of interest to check + :type qoi: :py:class:`soupy.ControlQoI` + :param x: The list of state, parameter, adjoint, and control variables + :type x: list of :py:class:`dolfin.Vector` or similar + :param du: The perturbation to the state variable + :type du: :py:class:`dolfin.Vector` or similar + :param order: Order of derivative to check. 1 for gradient, 2 for Hessian + :type order: int + :param delta: The finite difference step size + :type delta: float + :plotting: If :code:`true`, plots the finite difference Hessian and analytic Hessian + :type plotting: bool """ x2 = [None, None, None, None] x2[STATE] = dl.Vector(x[STATE]) diff --git a/soupy/utils/scipyCostWrapper.py b/soupy/utils/scipyCostWrapper.py index a820fc9..a87048f 100644 --- a/soupy/utils/scipyCostWrapper.py +++ b/soupy/utils/scipyCostWrapper.py @@ -16,11 +16,17 @@ class ScipyCostWrapper: """ - Class to interface the controlCostFunctional with a - scipy optimizer. Converts inputs to functions taking - and returning numpy arrays + Class to interface the controlCostFunctional with a \ + scipy optimizer. Converts inputs to functions taking \ + and returning :code:`numpy` arrays """ def __init__(self, controlCostFunctional, verbose=False): + """ + Constructor + + :param controlCostFunctional: The cost functional to wrap + :type controlCostFunctional: :py:class:`ControlCostFunctional` + """ self.cost_functional = controlCostFunctional self.z_help = self.cost_functional.generate_vector(CONTROL) self.g_help = self.cost_functional.generate_vector(CONTROL) @@ -29,6 +35,15 @@ def __init__(self, controlCostFunctional, verbose=False): self.n_grad = 0 def cost(self, z_np): + """ + Evaluates the cost functional at given control + + :param z_np: The control as a numpy array + :type z_np: :py:class:`numpy.ndarray` + + :returns: The value of the cost functional + :return type: float + """ self.z_help.set_local(z_np) self.z_help.apply("") cost_value = self.cost_functional.cost(self.z_help, order=0) @@ -40,6 +55,16 @@ def cost(self, z_np): return cost_value def grad(self, z_np): + """ + Evaluates the gradient of the cost functional at given control + + :param z_np: The control as a numpy array + :type z_np: :py:class:`numpy.ndarray` + + :returns: The gradient + :return type: :py:class:`numpy.ndarray` + """ + self.z_help.set_local(z_np) self.z_help.apply("") self.cost_functional.cost(self.z_help, order=1) @@ -54,9 +79,16 @@ def costHessian(self, z_np, zhat_np): pass def function(self): + """ + :returns: A function that evaluates the cost functional at a given control + """ return lambda x : self.cost(x) def jac(self): + """ + :returns: A function that evaluates the gradient of the \ + cost functional at a given control + """ return lambda x : self.grad(x) diff --git a/soupy/utils/stochasticCostFiniteDifference.py b/soupy/utils/stochasticCostFiniteDifference.py index 1649207..bd082db 100644 --- a/soupy/utils/stochasticCostFiniteDifference.py +++ b/soupy/utils/stochasticCostFiniteDifference.py @@ -20,6 +20,9 @@ from ..modeling.augmentedVector import AugmentedVector def stochasticCostFiniteDifference(pde_cost, z, dz, delta=1e-3, sample_size=1): + """ + Finite difference check for a stochastic cost function by fixing the random number generator seed + """ gz = pde_cost.generate_vector(CONTROL) if isinstance(z, AugmentedVector): @@ -49,6 +52,9 @@ def stochasticCostFiniteDifference(pde_cost, z, dz, delta=1e-3, sample_size=1): print("Finite diff gradient: %g" %fd_grad) def SAACostFiniteDifference(pde_cost, z, dz, delta=1e-3): + """ + Finite difference check for a deterministic/SAA cost functional + """ gz = pde_cost.generate_vector(CONTROL) if isinstance(z, AugmentedVector):