diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 88b44c9..1316496 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -9,7 +9,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - python-version: ['3.9', '3.10', '3.11', '3.12'] + python-version: ['3.12'] os: [ ubuntu-latest, macOS-latest, windows-latest ] steps: - uses: actions/checkout@v4 diff --git a/.gitignore b/.gitignore index ee48a4a..b9abd7b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ # Custom todo.txt +amisc_* # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md index 36274e0..b724de0 100644 --- a/docs/CONTRIBUTING.md +++ b/docs/CONTRIBUTING.md @@ -15,8 +15,8 @@ Pull requests are the best way to propose changes to the codebase (bug fixes, ne 1. Fork the repo and create your branch from `main`. 3. If you are adding a feature or making major changes, first create the [issue](https://github.com/eckelsjd/amisc/issues). -4. If you've added code that should be tested, add a [test](https://github.com/eckelsjd/amisc/tests). -5. If you've made major changes, update the [documentation](https://github.com/eckelsjd/amisc/docs). +4. If you've added code that should be tested, add to `/tests`. +5. If you've made major changes, update the `/docs`. 6. Ensure the test suite passes (`pdm run test`). 7. Make sure your code passes lint checks (coming soon). 8. Follow [Conventional commits](https://www.conventionalcommits.org/en/v1.0.0/) guidelines when adding a commit message. @@ -26,27 +26,35 @@ We strongly recommend using [pdm](https://github.com/pdm-project/pdm) to set up ```shell pip install --user pdm + # Fork the repo on Github + git clone https://github.com//amisc.git cd amisc -pdm install # Installs all dev and package dependencies (including amisc itself in -e mode) into local .venv +pdm install git checkout -b + # Make local changes + pdm run test # make sure tests pass git add -A git commit -m "Adding a bugfix or new feature" git push -u origin + # Go to Github and "Compare & Pull Request" on your fork # For your PR to be merged: # squash all your commits on your branch (interactively in an IDE most likely) # rebase to the top of origin/main to include new changes from others - git fetch - git rebase -i main your-branch # for example - # Resolve any conflicts - # Your history now looks something like this: - # o your-branch - # / - # ---o---o---o main + +git fetch +git rebase -i main your-branch # for example + +# Resolve any conflicts +# Your history now looks something like this: +# o your-branch +# / +# ---o---o---o main + # You can delete the branch and fork when your PR has been merged ``` @@ -63,5 +71,5 @@ Start or take part in community [discussions](https://github.com/eckelsjd/amisc/ By contributing, you agree that your contributions will be licensed under its GNU GPLv3 License. ## Releases -The package version is tracked at `amisc.__init__.__version__`. You should not edit this value. Maintainers -will decide when to make a release by bumping this and running `pdm run release`. \ No newline at end of file +The package version is tracked at `amisc.__init__.__version__`. You should not edit this value. The version will be +increased on a case-by-case basis and released depending on the changes being merged. \ No newline at end of file diff --git a/docs/README.md b/docs/README.md index 12ee7c2..3121984 100644 --- a/docs/README.md +++ b/docs/README.md @@ -1,7 +1,7 @@ ![Logo](https://raw.githubusercontent.com/eckelsjd/amisc/main/docs/assets/amisc_logo_text.svg) [![pdm-managed](https://img.shields.io/badge/pdm-managed-blueviolet)](https://pdm-project.org) [![PyPI](https://img.shields.io/pypi/v/amisc?logo=python&logoColor=%23cccccc)](https://pypi.org/project/amisc) -[![Python 3.7](https://img.shields.io/badge/python-3.7+-blue.svg?logo=python&logoColor=cccccc)](https://www.python.org/downloads/) +[![Python 3.12](https://img.shields.io/badge/python-3.12+-blue.svg?logo=python&logoColor=cccccc)](https://www.python.org/downloads/) ![Commits](https://img.shields.io/github/commit-activity/m/eckelsjd/amisc?logo=github) ![build](https://img.shields.io/github/actions/workflow/status/eckelsjd/amisc/deploy.yml?logo=github ) @@ -39,31 +39,13 @@ pdm sync # reads pdm.lock and sets up an identical venv ``` ## Quickstart -```python -from amisc.surrogates import SystemSurrogate -from amisc.utils import UniformRV -import numpy as np - -def fun1(x): - return x ** 2 - -def fun2(y): - return np.sin(y) * np.exp(y) - -x, y, z = UniformRV(0, 1, 'x'), UniformRV(0, 1, 'y'), UniformRV(0, 1, 'z') -model1 = {'name': 'model1', 'model': fun1, 'exo_in': ['x'], 'coupling_out': ['y']} -model2 = {'name': 'model2', 'model': fun2, 'coupling_in': ['y'], 'coupling_out': ['z']} - -system = SystemSurrogate([model1, model2], [x], [y, z]) -system.fit() - -xtest = system.sample_inputs(10) -ytest = system.predict(xtest) +```python title="amisc.examples.tutorial.py" +--8<-- "amisc/examples/tutorial.py:simple" ``` ## Contributing See the [contribution](CONTRIBUTING.md) guidelines. -## Reference +## Citations AMISC paper [[1](https://onlinelibrary.wiley.com/doi/full/10.1002/nme.6958)]. diff --git a/docs/assets/fire_sat.png b/docs/assets/fire_sat.png new file mode 100644 index 0000000..43e5436 Binary files /dev/null and b/docs/assets/fire_sat.png differ diff --git a/docs/css/extra.css b/docs/css/extra.css new file mode 100644 index 0000000..ce970fc --- /dev/null +++ b/docs/css/extra.css @@ -0,0 +1,7 @@ +.md-grid { + max-width: 1440px; +} + +:root > * { + --md-code-hl-keyword-color: #842e21; +} \ No newline at end of file diff --git a/docs/explanation.md b/docs/explanation.md deleted file mode 100644 index fe4d5d3..0000000 --- a/docs/explanation.md +++ /dev/null @@ -1 +0,0 @@ -Coming soon \ No newline at end of file diff --git a/docs/how-to-guides.md b/docs/how-to-guides.md index fe4d5d3..8d665bf 100644 --- a/docs/how-to-guides.md +++ b/docs/how-to-guides.md @@ -1 +1,65 @@ -Coming soon \ No newline at end of file +## Specifying model inputs and outputs +Coming soon. + +## Defining a component +Coming soon. + +## Making a model wrapper function +The examples in the [tutorial](tutorials.md) use the simple function call signatures `ret = func(x)`, where `x` is an `np.ndarray` and +`ret` is a dictionary with the required `y=output` key-value pair. If your model must be executed outside of Python +(such as in a separate `.exe` file), then you can write a Python wrapper function with the same call signature as above +and make any external calls you need inside the function (such as with `os.popen()`). You then pass the wrapper function +to `ComponentSpec` and `SystemSurrogate`. + +!!! Note "Requirements for your wrapper function" + - First argument `x` must be an `np.ndarray` of the model inputs whose **last** dimension is the number of inputs, i.e. `x.shape[-1] = x_dim`. + - You can choose to handle as many other dimensions as you want, i.e. `x.shape[:-1]`. The surrogate will handle the same number + of dimensions you give to your wrapper function (so that `model(x)` and `surrogate(x)` are functionally equivalent). We recommend you handle at least + 1 extra dimension, i.e. `x.shape = (N, x_dim)`. So your wrapper must handle `N` total sets of inputs at a time. The easiest way is to just + write a for loop over `N` and run your model for a single set of inputs at a time. + - Your wrapper function must expect the `x_dim` inputs in a specific order according to how you defined your system. All + system-level exogenous inputs (i.e. those in `system.exo_vars`) must be first and in the order you specified for + `ComponentSpec(exo_in=[first, second, ...])`. All coupling inputs that come from the outputs of other models are next. + Regardless of what order you chose in `ComponentSpec(coupling_in=[one, two, three,...]`, your wrapper **must** expect them + in _sorted_ order according to `system.coupling_vars`. For example, if `system.coupling_vars = [a, b, c]` and + `comp = ComponentSpec(wrapper, coupling_in=[c, a], exo_in=[d, e], coupling_out=[f])`, then `x_dim = 4` and your `wrapper` function + should expect the inputs in `x` to be ordered as `[d, e, a, c]`. + - If you want to pass in model fidelity indices (see $\alpha$ in [theory](theory.md) for details), they must be in the form of a `tuple`, + and your wrapper function should accept the `alpha=...` keyword argument. Specifying `alpha` allows managing a hierarchy of modeling fidelities, if applicable. + - You can pass any number of additional positional arguments. Specify these with `ComponentSpec(model_args=...)`. + - You can pass any number of keyword arguments. Specify these with `ComponentSpec(model_kwargs=...)`. + - If you want to save and keep track of the full output of your model (i.e. if it writes result files to disk), then + you can specify `ComponentSpec(save_output=True)`. When you do this, you must also specify `SystemSurrogate(..., save_dir='path/to/save/dir')`. + You will then get a folder called `save_dir/amisc_timestamp/components/`. This folder will be passed to your + wrapper function as the keyword argument `output_dir=`. Make sure your `wrapper` accepts this keyword (no need to specify it in `ComponentSpec(model_kwargs=...)`; this is done automatically). + You can then have your model write whatever it wants to this folder. You **must** then pass back the names of the files + you created via `ret=dict(files=[your output files, ...])`. The filenames must be in a list and match the order in + which the samples in `x` were executed by the model. + - To assist the adaptive training procedure, you can also optionally have your model compute and return its computational cost via + `ret=dict(cost=cpu_cost)`. The computational cost should be expressed in units of seconds of CPU time (not walltime!) for _one_ model evaluation. + If your model makes use of `n` CPUs in parallel, then the total CPU time would be `n` times the wall clock time. + - The return dictionary of your wrapper can include anything else you want outside of the three fields `(y, files, cost)` discussed here. + Any extra return values will be ignored by the system. + +!!! Example + ```python + def wrapper_func(x, alpha, *args, output_dir=None, **kwargs): + print(x.shape) # (..., x_dim) + + # Your code here, for example: + output = x ** 2 + output_files = ['output_1.json', 'output2.json', ...] + cpu_time = 42 # seconds for one model evaluation + + ret = dict(y=output, files=output_files, cost=cpu_time) + + return ret + ``` + +!!! Warning + Always specify the model at a _global_ scope, i.e. don't use `lambda` or nested functions. When saving to + file, only a symbolic reference to the function signature will be saved, which must be globally defined + when loading back from that save file. + +## Putting it all together +Coming soon. \ No newline at end of file diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js new file mode 100644 index 0000000..32c2132 --- /dev/null +++ b/docs/javascripts/mathjax.js @@ -0,0 +1,16 @@ +window.MathJax = { + tex: { + inlineMath: [['$', '$'], ["\\(", "\\)"]], + displayMath: [["\\[", "\\]"]], + processEscapes: true, + processEnvironments: true + }, + options: { + ignoreHtmlClass: ".*|", + processHtmlClass: "arithmatex" + } +}; + +document$.subscribe(() => { + MathJax.typesetPromise() +}) \ No newline at end of file diff --git a/docs/reference.md b/docs/reference.md deleted file mode 100644 index fe4d5d3..0000000 --- a/docs/reference.md +++ /dev/null @@ -1 +0,0 @@ -Coming soon \ No newline at end of file diff --git a/docs/reference/components.md b/docs/reference/components.md new file mode 100644 index 0000000..6cf06d6 --- /dev/null +++ b/docs/reference/components.md @@ -0,0 +1,3 @@ +::: amisc.component + options: + members_order: source \ No newline at end of file diff --git a/docs/reference/interpolators.md b/docs/reference/interpolators.md new file mode 100644 index 0000000..4a4d028 --- /dev/null +++ b/docs/reference/interpolators.md @@ -0,0 +1,4 @@ +::: amisc.interpolator + options: + filters: [""] + members_order: source \ No newline at end of file diff --git a/docs/reference/overview.md b/docs/reference/overview.md new file mode 100644 index 0000000..1a8dbe6 --- /dev/null +++ b/docs/reference/overview.md @@ -0,0 +1,91 @@ +The `amisc` package takes an object-oriented approach to building a surrogate of a multidisciplinary system. From the +bottom up, you have: + +- **variables** that serve as inputs and outputs for the models, +- **interpolators** that define a specific input → output mathematical relationship to interpolate a function, +- **components** that wrap a model for a single discipline, and a +- **system** that defines the connections between components in a multidisciplinary system. + +The variables, interpolators, and components all have abstract base classes, so that the **system** is ultimately +independent of the specific models, interpolation methods, or underlying variables. As such, the primary top-level object +that users of the `amisc` package will interact with is the `SystemSurrogate`. + +!!! Note + There are already pretty good implementations of the other abstractions that most users will not need to worry about, + but they are provided in this API reference for completeness. The abstractions allow new interpolation + (i.e. function approximation) methods to be implemented if desired, such as neural networks, kriging, etc. + +Here is a class diagram summary of this workflow: + +``` mermaid +classDiagram + namespace Core { + class SystemSurrogate { + +list[BaseRV] exo_vars + +list[BaseRV] coupling_vars + +int refine_level + +fit() + +predict(x) + +sample_inputs(size) + +insert_component(comp) + } + class ComponentSurrogate { + <> + +IndexSet index_set + +IndexSet candidate_set + +list[BaseRV] x_vars + +dict[str: BaseInterpolator] surrogates + +dict[str: float] misc_coeff + +predict(x) + +activate_index(alpha, beta) + +add_surrogate(alpha, beta) + +update_misc_coeff() + } + class BaseInterpolator { + <> + +tuple beta + +list[BaseRV] x_vars + +np.ndarray xi + +np.ndarray yi + +set_yi() + +refine() + +__call__(x) + } + } + class SparseGridSurrogate { + +np.ndarray x_grids + +dict xi_map + +dict yi_map + +get_tensor_grid(alpha, beta) + } + class LagrangeInterpolator { + +np.ndarray x_grids + +np.ndarray weights + +get_grid_sizes() + +leja_1d() + } + class BaseRV { + <> + +tuple bounds + +str units + +float nominal + +pdf(x) + +sample(size) + } + class UniformRV { + +str type + +get_uniform_bounds(nominal) + } + SystemSurrogate o-- "1..n" ComponentSurrogate + ComponentSurrogate o-- "1..n" BaseInterpolator + direction LR + ComponentSurrogate <|-- SparseGridSurrogate + BaseInterpolator <|-- LagrangeInterpolator + SparseGridSurrogate ..> LagrangeInterpolator + BaseRV <|-- UniformRV +``` +Note how the `SystemSurrogate` aggregates the `ComponentSurrogate`, which aggregates the `BaseInterpolator`. In other words, +interpolators can act independently of components, and components can act independently of systems. All three make use +of the random variables (these connections and some RVs are not shown for visual clarity). Currently, the only underlying surrogate +method that is implemented here is Lagrange polynomial interpolation (i.e. the `LagrangeInterpolator`). If one wanted +to use neural networks instead, the only change required is a new implementation of `BaseInterpolator`. \ No newline at end of file diff --git a/docs/reference/system.md b/docs/reference/system.md new file mode 100644 index 0000000..8b1c485 --- /dev/null +++ b/docs/reference/system.md @@ -0,0 +1 @@ +::: amisc.system \ No newline at end of file diff --git a/docs/reference/utilities.md b/docs/reference/utilities.md new file mode 100644 index 0000000..a96117e --- /dev/null +++ b/docs/reference/utilities.md @@ -0,0 +1,3 @@ +# Package utilities + +::: amisc.utils \ No newline at end of file diff --git a/docs/reference/variables.md b/docs/reference/variables.md new file mode 100644 index 0000000..37cfe9b --- /dev/null +++ b/docs/reference/variables.md @@ -0,0 +1 @@ +::: amisc.rv \ No newline at end of file diff --git a/docs/theory.md b/docs/theory.md new file mode 100644 index 0000000..fcc39fe --- /dev/null +++ b/docs/theory.md @@ -0,0 +1 @@ +Coming soon. \ No newline at end of file diff --git a/docs/tutorials.md b/docs/tutorials.md index fe4d5d3..989c12b 100644 --- a/docs/tutorials.md +++ b/docs/tutorials.md @@ -1 +1,26 @@ -Coming soon \ No newline at end of file +## Single component example +Here is an example of interpolating a simple quadratic function. +```python title="amisc.examples.tutorial.py" +--8<-- "amisc/examples/tutorial.py:single" +``` + +## Two component system +Here is a simple example of a two-component multidisciplinary system. +```python title="amisc.examples.tutorial.py" +--8<-- "amisc/examples/tutorial.py:simple" +``` +The first component computes $y=x\sin(\pi x)$. The second component takes the output of the first and computes +$z=1 / (1 + 25y^2)$. The system-level input is $x$ and the system-level outputs are $y$ and $z$. + +!!! Note + Each component always locally returns a dictionary with the output saved as `y=value`. This is not to be confused with the + _system-level_ `y` variable in this example. + +## Fire detection satellite +Here is an example of a three-component fire detection satellite system from [Chauduri (2018)](https://dspace.mit.edu/handle/1721.1/117036). +```python title="amisc.examples.tutorial.py" +--8<-- "amisc/examples/tutorial.py:fire_sat" +``` +We first generate a test set using the ground truth model predictions (and filter any bad values out). Then we train the +surrogate in 10 iterations, and finally plot some results. Here is the output of `plot_slice()`: +![Fire satellite system results](assets/fire_sat.png) \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index f6c6e20..28c2b3b 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -1,22 +1,93 @@ -site_name: Documentation for AMISC +site_name: Adaptive Multi-index Stochastic Collocation +site_url: https://eckelsjd.github.io/amisc/ +site_author: Joshua Eckels +repo_name: eckelsjd/amisc +repo_url: https://github.com/eckelsjd/amisc +copyright: Copyright © 2023 Joshua Eckels theme: - name: "material" + name: material logo: assets/amisc_logo.svg favicon: assets/amisc_logo.svg + palette: + primary: light blue + accent: red + font: + text: Roboto + code: Roboto Mono + features: + - content.code.annotate + - content.code.copy + - navigation.footer + - navigation.instant + - navigation.tabs + - navigation.top + - navigation.tracking + - search.highlight + - search.share + - search.suggest + - toc.follow plugins: + - search: + separator: '[\s\u200b\-_,:!=\[\]()"`/]+|\.(?!\d)|&[lg]t;|(?!\b)(?=[A-Z][a-z])' - mkdocstrings: handlers: python: + paths: [src] options: docstring_style: sphinx - docstring_section_style: list + docstring_section_style: spacy + merge_init_into_class: true + filters: ["!^_"] + show_symbol_type_heading: true # not available to public right now + group_by_category: true nav: - - Home: README.md - - tutorials.md - - How-to Guides: how-to-guides.md - - reference.md - - explanation.md - - Contributing: CONTRIBUTING.md \ No newline at end of file + - Home: + - Getting started: README.md + - Tutorials: tutorials.md + - How-to Guides: how-to-guides.md + - Contributing: CONTRIBUTING.md + - API Reference: + - Overview: reference/overview.md + - System: reference/system.md + - Components: reference/components.md + - Interpolators: reference/interpolators.md + - Variables: reference/variables.md + - Utilities: reference/utilities.md + - Theory: + - Home: theory.md + +extra: + social: + - icon: fontawesome/brands/github + link: https://github.com/eckelsjd + - icon: fontawesome/brands/python + link: https://pypi.org/project/amisc/ + +extra_javascript: + - javascripts/mathjax.js + - https://polyfill.io/v3/polyfill.min.js?features=es6 + - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + +#extra_css: +# - css/extra.css + +markdown_extensions: + - admonition + - pymdownx.superfences: + custom_fences: + - name: mermaid + class: mermaid + format: !!python/name:pymdownx.superfences.fence_code_format + - pymdownx.highlight: + anchor_linenums: true + line_spans: __span + pygments_lang_class: true + - pymdownx.arithmatex: + generic: true + - pymdownx.inlinehilite + - pymdownx.snippets: + dedent_subsections: true + base_path: ['.', './src'] diff --git a/pdm.lock b/pdm.lock index a450178..a972f2d 100644 --- a/pdm.lock +++ b/pdm.lock @@ -5,16 +5,13 @@ groups = ["default", "dev"] strategy = ["cross_platform"] lock_version = "4.4" -content_hash = "sha256:e8c7f0525bdb00302da5e754c8cda9ef49b43edbef1409f37cfce2ef09bd72ec" +content_hash = "sha256:bad2866334854360a4174574ed952149d5018f8759a3082f039fef97abab0d69" [[package]] name = "astroid" version = "3.0.1" requires_python = ">=3.8.0" summary = "An abstract syntax tree for Python with inference support." -dependencies = [ - "typing-extensions>=4.0.0; python_version < \"3.11\"", -] files = [ {file = "astroid-3.0.1-py3-none-any.whl", hash = "sha256:7d5895c9825e18079c5aeac0572bc2e4c83205c95d416e0b4fee8bc361d2d9ca"}, {file = "astroid-3.0.1.tar.gz", hash = "sha256:86b0bb7d7da0be1a7c4aedb7974e391b32d4ed89e33de6ed6902b4b15c97577e"}, @@ -136,6 +133,61 @@ files = [ {file = "colorama-0.4.6.tar.gz", hash = "sha256:08695f5cb7ed6e0531a20572697297273c47b8cae5a63ffc6d6ed5c201be6e44"}, ] +[[package]] +name = "contourpy" +version = "1.2.0" +requires_python = ">=3.9" +summary = "Python library for calculating contours of 2D quadrilateral grids" +dependencies = [ + "numpy<2.0,>=1.20", +] +files = [ + {file = "contourpy-1.2.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:0274c1cb63625972c0c007ab14dd9ba9e199c36ae1a231ce45d725cbcbfd10a8"}, + {file = "contourpy-1.2.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:ab459a1cbbf18e8698399c595a01f6dcc5c138220ca3ea9e7e6126232d102bb4"}, + {file = "contourpy-1.2.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:6fdd887f17c2f4572ce548461e4f96396681212d858cae7bd52ba3310bc6f00f"}, + {file = "contourpy-1.2.0-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:5d16edfc3fc09968e09ddffada434b3bf989bf4911535e04eada58469873e28e"}, + {file = "contourpy-1.2.0-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:1c203f617abc0dde5792beb586f827021069fb6d403d7f4d5c2b543d87edceb9"}, + {file = "contourpy-1.2.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b69303ceb2e4d4f146bf82fda78891ef7bcd80c41bf16bfca3d0d7eb545448aa"}, + {file = "contourpy-1.2.0-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:884c3f9d42d7218304bc74a8a7693d172685c84bd7ab2bab1ee567b769696df9"}, + {file = "contourpy-1.2.0-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:4a1b1208102be6e851f20066bf0e7a96b7d48a07c9b0cfe6d0d4545c2f6cadab"}, + {file = "contourpy-1.2.0-cp310-cp310-win32.whl", hash = "sha256:34b9071c040d6fe45d9826cbbe3727d20d83f1b6110d219b83eb0e2a01d79488"}, + {file = "contourpy-1.2.0-cp310-cp310-win_amd64.whl", hash = "sha256:bd2f1ae63998da104f16a8b788f685e55d65760cd1929518fd94cd682bf03e41"}, + {file = "contourpy-1.2.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:dd10c26b4eadae44783c45ad6655220426f971c61d9b239e6f7b16d5cdaaa727"}, + {file = "contourpy-1.2.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:5c6b28956b7b232ae801406e529ad7b350d3f09a4fde958dfdf3c0520cdde0dd"}, + {file = "contourpy-1.2.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ebeac59e9e1eb4b84940d076d9f9a6cec0064e241818bcb6e32124cc5c3e377a"}, + {file = "contourpy-1.2.0-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:139d8d2e1c1dd52d78682f505e980f592ba53c9f73bd6be102233e358b401063"}, + {file = "contourpy-1.2.0-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:1e9dc350fb4c58adc64df3e0703ab076f60aac06e67d48b3848c23647ae4310e"}, + {file = "contourpy-1.2.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:18fc2b4ed8e4a8fe849d18dce4bd3c7ea637758c6343a1f2bae1e9bd4c9f4686"}, + {file = "contourpy-1.2.0-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:16a7380e943a6d52472096cb7ad5264ecee36ed60888e2a3d3814991a0107286"}, + {file = "contourpy-1.2.0-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:8d8faf05be5ec8e02a4d86f616fc2a0322ff4a4ce26c0f09d9f7fb5330a35c95"}, + {file = "contourpy-1.2.0-cp311-cp311-win32.whl", hash = "sha256:67b7f17679fa62ec82b7e3e611c43a016b887bd64fb933b3ae8638583006c6d6"}, + {file = "contourpy-1.2.0-cp311-cp311-win_amd64.whl", hash = "sha256:99ad97258985328b4f207a5e777c1b44a83bfe7cf1f87b99f9c11d4ee477c4de"}, + {file = "contourpy-1.2.0-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:575bcaf957a25d1194903a10bc9f316c136c19f24e0985a2b9b5608bdf5dbfe0"}, + {file = "contourpy-1.2.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:9e6c93b5b2dbcedad20a2f18ec22cae47da0d705d454308063421a3b290d9ea4"}, + {file = "contourpy-1.2.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:464b423bc2a009088f19bdf1f232299e8b6917963e2b7e1d277da5041f33a779"}, + {file = "contourpy-1.2.0-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:68ce4788b7d93e47f84edd3f1f95acdcd142ae60bc0e5493bfd120683d2d4316"}, + {file = "contourpy-1.2.0-cp312-cp312-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:3d7d1f8871998cdff5d2ff6a087e5e1780139abe2838e85b0b46b7ae6cc25399"}, + {file = "contourpy-1.2.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:6e739530c662a8d6d42c37c2ed52a6f0932c2d4a3e8c1f90692ad0ce1274abe0"}, + {file = "contourpy-1.2.0-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:247b9d16535acaa766d03037d8e8fb20866d054d3c7fbf6fd1f993f11fc60ca0"}, + {file = "contourpy-1.2.0-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:461e3ae84cd90b30f8d533f07d87c00379644205b1d33a5ea03381edc4b69431"}, + {file = "contourpy-1.2.0-cp312-cp312-win32.whl", hash = "sha256:1c2559d6cffc94890b0529ea7eeecc20d6fadc1539273aa27faf503eb4656d8f"}, + {file = "contourpy-1.2.0-cp312-cp312-win_amd64.whl", hash = "sha256:491b1917afdd8638a05b611a56d46587d5a632cabead889a5440f7c638bc6ed9"}, + {file = "contourpy-1.2.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:5fd1810973a375ca0e097dee059c407913ba35723b111df75671a1976efa04bc"}, + {file = "contourpy-1.2.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:999c71939aad2780f003979b25ac5b8f2df651dac7b38fb8ce6c46ba5abe6ae9"}, + {file = "contourpy-1.2.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b7caf9b241464c404613512d5594a6e2ff0cc9cb5615c9475cc1d9b514218ae8"}, + {file = "contourpy-1.2.0-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:266270c6f6608340f6c9836a0fb9b367be61dde0c9a9a18d5ece97774105ff3e"}, + {file = "contourpy-1.2.0-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:dbd50d0a0539ae2e96e537553aff6d02c10ed165ef40c65b0e27e744a0f10af8"}, + {file = "contourpy-1.2.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:11f8d2554e52f459918f7b8e6aa20ec2a3bce35ce95c1f0ef4ba36fbda306df5"}, + {file = "contourpy-1.2.0-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:ce96dd400486e80ac7d195b2d800b03e3e6a787e2a522bfb83755938465a819e"}, + {file = "contourpy-1.2.0-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:6d3364b999c62f539cd403f8123ae426da946e142312a514162adb2addd8d808"}, + {file = "contourpy-1.2.0-cp39-cp39-win32.whl", hash = "sha256:1c88dfb9e0c77612febebb6ac69d44a8d81e3dc60f993215425b62c1161353f4"}, + {file = "contourpy-1.2.0-cp39-cp39-win_amd64.whl", hash = "sha256:78e6ad33cf2e2e80c5dfaaa0beec3d61face0fb650557100ee36db808bfa6843"}, + {file = "contourpy-1.2.0-pp39-pypy39_pp73-macosx_10_9_x86_64.whl", hash = "sha256:be16975d94c320432657ad2402f6760990cb640c161ae6da1363051805fa8108"}, + {file = "contourpy-1.2.0-pp39-pypy39_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b95a225d4948b26a28c08307a60ac00fb8671b14f2047fc5476613252a129776"}, + {file = "contourpy-1.2.0-pp39-pypy39_pp73-win_amd64.whl", hash = "sha256:0d7e03c0f9a4f90dc18d4e77e9ef4ec7b7bbb437f7f675be8e530d65ae6ef956"}, + {file = "contourpy-1.2.0.tar.gz", hash = "sha256:171f311cb758de7da13fc53af221ae47a5877be5a0843a9fe150818c51ed276a"}, +] + [[package]] name = "coverage" version = "7.3.2" @@ -184,7 +236,6 @@ requires_python = ">=3.8" summary = "Code coverage measurement for Python" dependencies = [ "coverage==7.3.2", - "tomli; python_full_version <= \"3.11.0a6\"", ] files = [ {file = "coverage-7.3.2-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:d872145f3a3231a5f20fd48500274d7df222e291d90baa2026cc5152b7ce86bf"}, @@ -221,6 +272,16 @@ files = [ {file = "coverage-7.3.2.tar.gz", hash = "sha256:be32ad29341b0170e795ca590e1c07e81fc061cb5b10c74ce7203491484404ef"}, ] +[[package]] +name = "cycler" +version = "0.12.1" +requires_python = ">=3.8" +summary = "Composable style cycles" +files = [ + {file = "cycler-0.12.1-py3-none-any.whl", hash = "sha256:85cef7cff222d8644161529808465972e51340599459b8ac3ccbac5a854e0d30"}, + {file = "cycler-0.12.1.tar.gz", hash = "sha256:88bb128f02ba341da8ef447245a9e138fae777f6a23943da4540077d3601eb1c"}, +] + [[package]] name = "dill" version = "0.3.7" @@ -232,13 +293,45 @@ files = [ ] [[package]] -name = "exceptiongroup" -version = "1.2.0" -requires_python = ">=3.7" -summary = "Backport of PEP 654 (exception groups)" +name = "fonttools" +version = "4.46.0" +requires_python = ">=3.8" +summary = "Tools to manipulate font files" files = [ - {file = "exceptiongroup-1.2.0-py3-none-any.whl", hash = "sha256:4bfd3996ac73b41e9b9628b04e079f193850720ea5945fc96a08633c66912f14"}, - {file = "exceptiongroup-1.2.0.tar.gz", hash = "sha256:91f5c769735f051a4290d52edd0858999b57e5876e9f85937691bd4c9fa3ed68"}, + {file = "fonttools-4.46.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:d4e69e2c7f93b695d2e6f18f709d501d945f65c1d237dafaabdd23cd935a5276"}, + {file = "fonttools-4.46.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:25852f0c63df0af022f698464a4a80f7d1d5bd974bcd22f995f6b4ad198e32dd"}, + {file = "fonttools-4.46.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:adab73618d0a328b203a0e242b3eba60a2b5662d9cb2bd16ed9c52af8a7d86af"}, + {file = "fonttools-4.46.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2cf923a4a556ab4cc4c52f69a4a2db624cf5a2cf360394368b40c5152fe3321e"}, + {file = "fonttools-4.46.0-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:87c214197712cc14fd2a4621efce2a9c501a77041232b789568149a8a3161517"}, + {file = "fonttools-4.46.0-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:156ae342a1ed1fe38e180de471e98fbf5b2b6ae280fa3323138569c4ca215844"}, + {file = "fonttools-4.46.0-cp310-cp310-win32.whl", hash = "sha256:c506e3d3a9e898caee4dc094f34b49c5566870d5a2d1ca2125f0a9f35ecc2205"}, + {file = "fonttools-4.46.0-cp310-cp310-win_amd64.whl", hash = "sha256:f8bc3973ed58893c4107993e0a7ae34901cb572b5e798249cbef35d30801ffd4"}, + {file = "fonttools-4.46.0-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:982f69855ac258260f51048d9e0c53c5f19881138cc7ca06deb38dc4b97404b6"}, + {file = "fonttools-4.46.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:2c23c59d321d62588620f2255cf951270bf637d88070f38ed8b5e5558775b86c"}, + {file = "fonttools-4.46.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a0e94244ec24a940ecfbe5b31c975c8a575d5ed2d80f9a280ce3b21fa5dc9c34"}, + {file = "fonttools-4.46.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1a9f9cdd7ef63d1b8ac90db335762451452426b3207abd79f60da510cea62da5"}, + {file = "fonttools-4.46.0-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:ca9eceebe70035b057ce549e2054cad73e95cac3fe91a9d827253d1c14618204"}, + {file = "fonttools-4.46.0-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:8be6adfa4e15977075278dd0a0bae74dec59be7b969b5ceed93fb86af52aa5be"}, + {file = "fonttools-4.46.0-cp311-cp311-win32.whl", hash = "sha256:7b5636f5706d49f13b6d610fe54ee662336cdf56b5a6f6683c0b803e23d826d2"}, + {file = "fonttools-4.46.0-cp311-cp311-win_amd64.whl", hash = "sha256:49ea0983e55fd7586a809787cd4644a7ae471e53ab8ddc016f9093b400e32646"}, + {file = "fonttools-4.46.0-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:7b460720ce81773da1a3e7cc964c48e1e11942b280619582a897fa0117b56a62"}, + {file = "fonttools-4.46.0-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:8bee9f4fc8c99824a424ae45c789ee8c67cb84f8e747afa7f83b7d3cef439c3b"}, + {file = "fonttools-4.46.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d3d7b96aba96e05e8c911ce2dfc5acc6a178b8f44f6aa69371ab91aa587563da"}, + {file = "fonttools-4.46.0-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9e6aeb5c340416d11a3209d75c48d13e72deea9e1517837dd1522c1fd1f17c11"}, + {file = "fonttools-4.46.0-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:c779f8701deedf41908f287aeb775b8a6f59875ad1002b98ac6034ae4ddc1b7b"}, + {file = "fonttools-4.46.0-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:ce199227ce7921eaafdd4f96536f16b232d6b580ce74ce337de544bf06cb2752"}, + {file = "fonttools-4.46.0-cp312-cp312-win32.whl", hash = "sha256:1c9937c4dd1061afd22643389445fabda858af5e805860ec3082a4bc07c7a720"}, + {file = "fonttools-4.46.0-cp312-cp312-win_amd64.whl", hash = "sha256:a9fa52ef8fd14d7eb3d813e1451e7ace3e1eebfa9b7237d3f81fee8f3de6a114"}, + {file = "fonttools-4.46.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:c9a0e422ab79e5cb2b47913be6a4b5fd20c4c7ac34a24f3691a4e099e965e0b8"}, + {file = "fonttools-4.46.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:13ac0cba2fc63fa4b232f2a7971f35f35c6eaf10bd1271fa96d4ce6253a8acfd"}, + {file = "fonttools-4.46.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:795150d5edc595e1a2cfb3d65e8f4f3d027704fc2579f8990d381bef6b188eb6"}, + {file = "fonttools-4.46.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d00fc63131dcac6b25f50a5a129758438317e54e3ce5587163f7058de4b0e933"}, + {file = "fonttools-4.46.0-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:3033b55f401a622de2630b3982234d97219d89b058607b87927eccb0f922313c"}, + {file = "fonttools-4.46.0-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:e26e7fb908ae4f622813e7cb32cd2db6c24e3122bb3b98f25e832a2fe0e7e228"}, + {file = "fonttools-4.46.0-cp39-cp39-win32.whl", hash = "sha256:2d0eba685938c603f2f648dfc0aadbf8c6a4fe1c7ca608c2970a6ef39e00f254"}, + {file = "fonttools-4.46.0-cp39-cp39-win_amd64.whl", hash = "sha256:5200b01f463d97cc2b7ff8a1e3584151f4413e98cb8419da5f17d1dbb84cc214"}, + {file = "fonttools-4.46.0-py3-none-any.whl", hash = "sha256:5b627ed142398ea9202bd752c04311592558964d1a765fb2f78dc441a05633f4"}, + {file = "fonttools-4.46.0.tar.gz", hash = "sha256:2ae45716c27a41807d58a9f3f59983bdc8c0a46cb259e4450ab7e196253a9853"}, ] [[package]] @@ -279,6 +372,19 @@ files = [ {file = "GitPython-3.1.40.tar.gz", hash = "sha256:22b126e9ffb671fdd0c129796343a02bf67bf2994b35449ffc9321aa755e18a4"}, ] +[[package]] +name = "griffe" +version = "0.38.1" +requires_python = ">=3.8" +summary = "Signatures for entire Python programs. Extract the structure, the frame, the skeleton of your project, to generate API documentation or find breaking changes in your API." +dependencies = [ + "colorama>=0.4", +] +files = [ + {file = "griffe-0.38.1-py3-none-any.whl", hash = "sha256:334c79d3b5964ade65c05dfcaf53518c576dedd387aaba5c9fd71212f34f1483"}, + {file = "griffe-0.38.1.tar.gz", hash = "sha256:bd68d7da7f3d87bc57eb9962b250db123efd9bbcc06c11c1a91b6e583b2a9361"}, +] + [[package]] name = "idna" version = "3.6" @@ -289,19 +395,6 @@ files = [ {file = "idna-3.6.tar.gz", hash = "sha256:9ecdbbd083b06798ae1e86adcbfe8ab1479cf864e4ee30fe4e46a003d12491ca"}, ] -[[package]] -name = "importlib-metadata" -version = "7.0.0" -requires_python = ">=3.8" -summary = "Read metadata from Python packages" -dependencies = [ - "zipp>=0.5", -] -files = [ - {file = "importlib_metadata-7.0.0-py3-none-any.whl", hash = "sha256:d97503976bb81f40a193d41ee6570868479c69d5068651eb039c40d850c59d67"}, - {file = "importlib_metadata-7.0.0.tar.gz", hash = "sha256:7fc841f8b8332803464e5dc1c63a2e59121f46ca186c0e2e182e80bf8c1319f7"}, -] - [[package]] name = "iniconfig" version = "2.0.0" @@ -335,14 +428,105 @@ files = [ {file = "Jinja2-3.1.2.tar.gz", hash = "sha256:31351a702a408a9e7595a8fc6150fc3f43bb6bf7e319770cbc0db9df9437e852"}, ] +[[package]] +name = "joblib" +version = "1.3.2" +requires_python = ">=3.7" +summary = "Lightweight pipelining with Python functions" +files = [ + {file = "joblib-1.3.2-py3-none-any.whl", hash = "sha256:ef4331c65f239985f3f2220ecc87db222f08fd22097a3dd5698f693875f8cbb9"}, + {file = "joblib-1.3.2.tar.gz", hash = "sha256:92f865e621e17784e7955080b6d042489e3b8e294949cc44c6eac304f59772b1"}, +] + +[[package]] +name = "kiwisolver" +version = "1.4.5" +requires_python = ">=3.7" +summary = "A fast implementation of the Cassowary constraint solver" +files = [ + {file = "kiwisolver-1.4.5-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:05703cf211d585109fcd72207a31bb170a0f22144d68298dc5e61b3c946518af"}, + {file = "kiwisolver-1.4.5-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:146d14bebb7f1dc4d5fbf74f8a6cb15ac42baadee8912eb84ac0b3b2a3dc6ac3"}, + {file = "kiwisolver-1.4.5-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:6ef7afcd2d281494c0a9101d5c571970708ad911d028137cd558f02b851c08b4"}, + {file = "kiwisolver-1.4.5-cp310-cp310-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:9eaa8b117dc8337728e834b9c6e2611f10c79e38f65157c4c38e9400286f5cb1"}, + {file = "kiwisolver-1.4.5-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:ec20916e7b4cbfb1f12380e46486ec4bcbaa91a9c448b97023fde0d5bbf9e4ff"}, + {file = "kiwisolver-1.4.5-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:39b42c68602539407884cf70d6a480a469b93b81b7701378ba5e2328660c847a"}, + {file = "kiwisolver-1.4.5-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:aa12042de0171fad672b6c59df69106d20d5596e4f87b5e8f76df757a7c399aa"}, + {file = "kiwisolver-1.4.5-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:2a40773c71d7ccdd3798f6489aaac9eee213d566850a9533f8d26332d626b82c"}, + {file = "kiwisolver-1.4.5-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:19df6e621f6d8b4b9c4d45f40a66839294ff2bb235e64d2178f7522d9170ac5b"}, + {file = "kiwisolver-1.4.5-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:83d78376d0d4fd884e2c114d0621624b73d2aba4e2788182d286309ebdeed770"}, + {file = "kiwisolver-1.4.5-cp310-cp310-musllinux_1_1_ppc64le.whl", hash = "sha256:e391b1f0a8a5a10ab3b9bb6afcfd74f2175f24f8975fb87ecae700d1503cdee0"}, + {file = "kiwisolver-1.4.5-cp310-cp310-musllinux_1_1_s390x.whl", hash = "sha256:852542f9481f4a62dbb5dd99e8ab7aedfeb8fb6342349a181d4036877410f525"}, + {file = "kiwisolver-1.4.5-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:59edc41b24031bc25108e210c0def6f6c2191210492a972d585a06ff246bb79b"}, + {file = "kiwisolver-1.4.5-cp310-cp310-win32.whl", hash = "sha256:a6aa6315319a052b4ee378aa171959c898a6183f15c1e541821c5c59beaa0238"}, + {file = "kiwisolver-1.4.5-cp310-cp310-win_amd64.whl", hash = "sha256:d0ef46024e6a3d79c01ff13801cb19d0cad7fd859b15037aec74315540acc276"}, + {file = "kiwisolver-1.4.5-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:11863aa14a51fd6ec28688d76f1735f8f69ab1fabf388851a595d0721af042f5"}, + {file = "kiwisolver-1.4.5-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:8ab3919a9997ab7ef2fbbed0cc99bb28d3c13e6d4b1ad36e97e482558a91be90"}, + {file = "kiwisolver-1.4.5-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:fcc700eadbbccbf6bc1bcb9dbe0786b4b1cb91ca0dcda336eef5c2beed37b797"}, + {file = "kiwisolver-1.4.5-cp311-cp311-manylinux_2_12_i686.manylinux2010_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:dfdd7c0b105af050eb3d64997809dc21da247cf44e63dc73ff0fd20b96be55a9"}, + {file = "kiwisolver-1.4.5-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:76c6a5964640638cdeaa0c359382e5703e9293030fe730018ca06bc2010c4437"}, + {file = "kiwisolver-1.4.5-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:bbea0db94288e29afcc4c28afbf3a7ccaf2d7e027489c449cf7e8f83c6346eb9"}, + {file = "kiwisolver-1.4.5-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:ceec1a6bc6cab1d6ff5d06592a91a692f90ec7505d6463a88a52cc0eb58545da"}, + {file = "kiwisolver-1.4.5-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:040c1aebeda72197ef477a906782b5ab0d387642e93bda547336b8957c61022e"}, + {file = "kiwisolver-1.4.5-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:f91de7223d4c7b793867797bacd1ee53bfe7359bd70d27b7b58a04efbb9436c8"}, + {file = "kiwisolver-1.4.5-cp311-cp311-musllinux_1_1_i686.whl", hash = "sha256:faae4860798c31530dd184046a900e652c95513796ef51a12bc086710c2eec4d"}, + {file = "kiwisolver-1.4.5-cp311-cp311-musllinux_1_1_ppc64le.whl", hash = "sha256:b0157420efcb803e71d1b28e2c287518b8808b7cf1ab8af36718fd0a2c453eb0"}, + {file = "kiwisolver-1.4.5-cp311-cp311-musllinux_1_1_s390x.whl", hash = "sha256:06f54715b7737c2fecdbf140d1afb11a33d59508a47bf11bb38ecf21dc9ab79f"}, + {file = "kiwisolver-1.4.5-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:fdb7adb641a0d13bdcd4ef48e062363d8a9ad4a182ac7647ec88f695e719ae9f"}, + {file = "kiwisolver-1.4.5-cp311-cp311-win32.whl", hash = "sha256:bb86433b1cfe686da83ce32a9d3a8dd308e85c76b60896d58f082136f10bffac"}, + {file = "kiwisolver-1.4.5-cp311-cp311-win_amd64.whl", hash = "sha256:6c08e1312a9cf1074d17b17728d3dfce2a5125b2d791527f33ffbe805200a355"}, + {file = "kiwisolver-1.4.5-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:32d5cf40c4f7c7b3ca500f8985eb3fb3a7dfc023215e876f207956b5ea26632a"}, + {file = "kiwisolver-1.4.5-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:f846c260f483d1fd217fe5ed7c173fb109efa6b1fc8381c8b7552c5781756192"}, + {file = "kiwisolver-1.4.5-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:5ff5cf3571589b6d13bfbfd6bcd7a3f659e42f96b5fd1c4830c4cf21d4f5ef45"}, + {file = "kiwisolver-1.4.5-cp312-cp312-manylinux_2_12_i686.manylinux2010_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:7269d9e5f1084a653d575c7ec012ff57f0c042258bf5db0954bf551c158466e7"}, + {file = "kiwisolver-1.4.5-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:da802a19d6e15dffe4b0c24b38b3af68e6c1a68e6e1d8f30148c83864f3881db"}, + {file = "kiwisolver-1.4.5-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:3aba7311af82e335dd1e36ffff68aaca609ca6290c2cb6d821a39aa075d8e3ff"}, + {file = "kiwisolver-1.4.5-cp312-cp312-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:763773d53f07244148ccac5b084da5adb90bfaee39c197554f01b286cf869228"}, + {file = "kiwisolver-1.4.5-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2270953c0d8cdab5d422bee7d2007f043473f9d2999631c86a223c9db56cbd16"}, + {file = "kiwisolver-1.4.5-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:d099e745a512f7e3bbe7249ca835f4d357c586d78d79ae8f1dcd4d8adeb9bda9"}, + {file = "kiwisolver-1.4.5-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:74db36e14a7d1ce0986fa104f7d5637aea5c82ca6326ed0ec5694280942d1162"}, + {file = "kiwisolver-1.4.5-cp312-cp312-musllinux_1_1_ppc64le.whl", hash = "sha256:7e5bab140c309cb3a6ce373a9e71eb7e4873c70c2dda01df6820474f9889d6d4"}, + {file = "kiwisolver-1.4.5-cp312-cp312-musllinux_1_1_s390x.whl", hash = "sha256:0f114aa76dc1b8f636d077979c0ac22e7cd8f3493abbab152f20eb8d3cda71f3"}, + {file = "kiwisolver-1.4.5-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:88a2df29d4724b9237fc0c6eaf2a1adae0cdc0b3e9f4d8e7dc54b16812d2d81a"}, + {file = "kiwisolver-1.4.5-cp312-cp312-win32.whl", hash = "sha256:72d40b33e834371fd330fb1472ca19d9b8327acb79a5821d4008391db8e29f20"}, + {file = "kiwisolver-1.4.5-cp312-cp312-win_amd64.whl", hash = "sha256:2c5674c4e74d939b9d91dda0fae10597ac7521768fec9e399c70a1f27e2ea2d9"}, + {file = "kiwisolver-1.4.5-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:9407b6a5f0d675e8a827ad8742e1d6b49d9c1a1da5d952a67d50ef5f4170b18d"}, + {file = "kiwisolver-1.4.5-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:15568384086b6df3c65353820a4473575dbad192e35010f622c6ce3eebd57af9"}, + {file = "kiwisolver-1.4.5-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:0dc9db8e79f0036e8173c466d21ef18e1befc02de8bf8aa8dc0813a6dc8a7046"}, + {file = "kiwisolver-1.4.5-cp39-cp39-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:cdc8a402aaee9a798b50d8b827d7ecf75edc5fb35ea0f91f213ff927c15f4ff0"}, + {file = "kiwisolver-1.4.5-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:6c3bd3cde54cafb87d74d8db50b909705c62b17c2099b8f2e25b461882e544ff"}, + {file = "kiwisolver-1.4.5-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:955e8513d07a283056b1396e9a57ceddbd272d9252c14f154d450d227606eb54"}, + {file = "kiwisolver-1.4.5-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:346f5343b9e3f00b8db8ba359350eb124b98c99efd0b408728ac6ebf38173958"}, + {file = "kiwisolver-1.4.5-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:b9098e0049e88c6a24ff64545cdfc50807818ba6c1b739cae221bbbcbc58aad3"}, + {file = "kiwisolver-1.4.5-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:00bd361b903dc4bbf4eb165f24d1acbee754fce22ded24c3d56eec268658a5cf"}, + {file = "kiwisolver-1.4.5-cp39-cp39-musllinux_1_1_i686.whl", hash = "sha256:7b8b454bac16428b22560d0a1cf0a09875339cab69df61d7805bf48919415901"}, + {file = "kiwisolver-1.4.5-cp39-cp39-musllinux_1_1_ppc64le.whl", hash = "sha256:f1d072c2eb0ad60d4c183f3fb44ac6f73fb7a8f16a2694a91f988275cbf352f9"}, + {file = "kiwisolver-1.4.5-cp39-cp39-musllinux_1_1_s390x.whl", hash = "sha256:31a82d498054cac9f6d0b53d02bb85811185bcb477d4b60144f915f3b3126342"}, + {file = "kiwisolver-1.4.5-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:6512cb89e334e4700febbffaaa52761b65b4f5a3cf33f960213d5656cea36a77"}, + {file = "kiwisolver-1.4.5-cp39-cp39-win32.whl", hash = "sha256:9db8ea4c388fdb0f780fe91346fd438657ea602d58348753d9fb265ce1bca67f"}, + {file = "kiwisolver-1.4.5-cp39-cp39-win_amd64.whl", hash = "sha256:59415f46a37f7f2efeec758353dd2eae1b07640d8ca0f0c42548ec4125492635"}, + {file = "kiwisolver-1.4.5-pp37-pypy37_pp73-macosx_10_9_x86_64.whl", hash = "sha256:5c7b3b3a728dc6faf3fc372ef24f21d1e3cee2ac3e9596691d746e5a536de920"}, + {file = "kiwisolver-1.4.5-pp37-pypy37_pp73-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:620ced262a86244e2be10a676b646f29c34537d0d9cc8eb26c08f53d98013390"}, + {file = "kiwisolver-1.4.5-pp37-pypy37_pp73-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:378a214a1e3bbf5ac4a8708304318b4f890da88c9e6a07699c4ae7174c09a68d"}, + {file = "kiwisolver-1.4.5-pp37-pypy37_pp73-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:aaf7be1207676ac608a50cd08f102f6742dbfc70e8d60c4db1c6897f62f71523"}, + {file = "kiwisolver-1.4.5-pp37-pypy37_pp73-win_amd64.whl", hash = "sha256:ba55dce0a9b8ff59495ddd050a0225d58bd0983d09f87cfe2b6aec4f2c1234e4"}, + {file = "kiwisolver-1.4.5-pp38-pypy38_pp73-macosx_10_9_x86_64.whl", hash = "sha256:fd32ea360bcbb92d28933fc05ed09bffcb1704ba3fc7942e81db0fd4f81a7892"}, + {file = "kiwisolver-1.4.5-pp38-pypy38_pp73-manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:5e7139af55d1688f8b960ee9ad5adafc4ac17c1c473fe07133ac092310d76544"}, + {file = "kiwisolver-1.4.5-pp38-pypy38_pp73-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:dced8146011d2bc2e883f9bd68618b8247387f4bbec46d7392b3c3b032640126"}, + {file = "kiwisolver-1.4.5-pp38-pypy38_pp73-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:c9bf3325c47b11b2e51bca0824ea217c7cd84491d8ac4eefd1e409705ef092bd"}, + {file = "kiwisolver-1.4.5-pp38-pypy38_pp73-win_amd64.whl", hash = "sha256:5794cf59533bc3f1b1c821f7206a3617999db9fbefc345360aafe2e067514929"}, + {file = "kiwisolver-1.4.5-pp39-pypy39_pp73-macosx_10_9_x86_64.whl", hash = "sha256:e368f200bbc2e4f905b8e71eb38b3c04333bddaa6a2464a6355487b02bb7fb09"}, + {file = "kiwisolver-1.4.5-pp39-pypy39_pp73-manylinux_2_12_i686.manylinux2010_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:e5d706eba36b4c4d5bc6c6377bb6568098765e990cfc21ee16d13963fab7b3e7"}, + {file = "kiwisolver-1.4.5-pp39-pypy39_pp73-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:85267bd1aa8880a9c88a8cb71e18d3d64d2751a790e6ca6c27b8ccc724bcd5ad"}, + {file = "kiwisolver-1.4.5-pp39-pypy39_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:210ef2c3a1f03272649aff1ef992df2e724748918c4bc2d5a90352849eb40bea"}, + {file = "kiwisolver-1.4.5-pp39-pypy39_pp73-win_amd64.whl", hash = "sha256:11d011a7574eb3b82bcc9c1a1d35c1d7075677fdd15de527d91b46bd35e935ee"}, + {file = "kiwisolver-1.4.5.tar.gz", hash = "sha256:e57e563a57fb22a142da34f38acc2fc1a5c864bc29ca1517a88abc963e60d6ec"}, +] + [[package]] name = "markdown" version = "3.5.1" requires_python = ">=3.8" summary = "Python implementation of John Gruber's Markdown." -dependencies = [ - "importlib-metadata>=4.4; python_version < \"3.10\"", -] files = [ {file = "Markdown-3.5.1-py3-none-any.whl", hash = "sha256:5874b47d4ee3f0b14d764324d2c94c03ea66bee56f2d929da9f2508d65e722dc"}, {file = "Markdown-3.5.1.tar.gz", hash = "sha256:b65d7beb248dc22f2e8a31fb706d93798093c308dc1aba295aedeb9d41a813bd"}, @@ -387,6 +571,53 @@ files = [ {file = "MarkupSafe-2.1.3.tar.gz", hash = "sha256:af598ed32d6ae86f1b747b82783958b1a4ab8f617b06fe68795c7f026abbdcad"}, ] +[[package]] +name = "matplotlib" +version = "3.8.2" +requires_python = ">=3.9" +summary = "Python plotting package" +dependencies = [ + "contourpy>=1.0.1", + "cycler>=0.10", + "fonttools>=4.22.0", + "kiwisolver>=1.3.1", + "numpy<2,>=1.21", + "packaging>=20.0", + "pillow>=8", + "pyparsing>=2.3.1", + "python-dateutil>=2.7", +] +files = [ + {file = "matplotlib-3.8.2-cp310-cp310-macosx_10_12_x86_64.whl", hash = "sha256:09796f89fb71a0c0e1e2f4bdaf63fb2cefc84446bb963ecdeb40dfee7dfa98c7"}, + {file = "matplotlib-3.8.2-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:6f9c6976748a25e8b9be51ea028df49b8e561eed7809146da7a47dbecebab367"}, + {file = "matplotlib-3.8.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b78e4f2cedf303869b782071b55fdde5987fda3038e9d09e58c91cc261b5ad18"}, + {file = "matplotlib-3.8.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:4e208f46cf6576a7624195aa047cb344a7f802e113bb1a06cfd4bee431de5e31"}, + {file = "matplotlib-3.8.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:46a569130ff53798ea5f50afce7406e91fdc471ca1e0e26ba976a8c734c9427a"}, + {file = "matplotlib-3.8.2-cp310-cp310-win_amd64.whl", hash = "sha256:830f00640c965c5b7f6bc32f0d4ce0c36dfe0379f7dd65b07a00c801713ec40a"}, + {file = "matplotlib-3.8.2-cp311-cp311-macosx_10_12_x86_64.whl", hash = "sha256:d86593ccf546223eb75a39b44c32788e6f6440d13cfc4750c1c15d0fcb850b63"}, + {file = "matplotlib-3.8.2-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:9a5430836811b7652991939012f43d2808a2db9b64ee240387e8c43e2e5578c8"}, + {file = "matplotlib-3.8.2-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b9576723858a78751d5aacd2497b8aef29ffea6d1c95981505877f7ac28215c6"}, + {file = "matplotlib-3.8.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5ba9cbd8ac6cf422f3102622b20f8552d601bf8837e49a3afed188d560152788"}, + {file = "matplotlib-3.8.2-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:03f9d160a29e0b65c0790bb07f4f45d6a181b1ac33eb1bb0dd225986450148f0"}, + {file = "matplotlib-3.8.2-cp311-cp311-win_amd64.whl", hash = "sha256:3773002da767f0a9323ba1a9b9b5d00d6257dbd2a93107233167cfb581f64717"}, + {file = "matplotlib-3.8.2-cp312-cp312-macosx_10_12_x86_64.whl", hash = "sha256:4c318c1e95e2f5926fba326f68177dee364aa791d6df022ceb91b8221bd0a627"}, + {file = "matplotlib-3.8.2-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:091275d18d942cf1ee9609c830a1bc36610607d8223b1b981c37d5c9fc3e46a4"}, + {file = "matplotlib-3.8.2-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:1b0f3b8ea0e99e233a4bcc44590f01604840d833c280ebb8fe5554fd3e6cfe8d"}, + {file = "matplotlib-3.8.2-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d7b1704a530395aaf73912be741c04d181f82ca78084fbd80bc737be04848331"}, + {file = "matplotlib-3.8.2-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:533b0e3b0c6768eef8cbe4b583731ce25a91ab54a22f830db2b031e83cca9213"}, + {file = "matplotlib-3.8.2-cp312-cp312-win_amd64.whl", hash = "sha256:0f4fc5d72b75e2c18e55eb32292659cf731d9d5b312a6eb036506304f4675630"}, + {file = "matplotlib-3.8.2-cp39-cp39-macosx_10_12_x86_64.whl", hash = "sha256:deaed9ad4da0b1aea77fe0aa0cebb9ef611c70b3177be936a95e5d01fa05094f"}, + {file = "matplotlib-3.8.2-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:172f4d0fbac3383d39164c6caafd3255ce6fa58f08fc392513a0b1d3b89c4f89"}, + {file = "matplotlib-3.8.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:c7d36c2209d9136cd8e02fab1c0ddc185ce79bc914c45054a9f514e44c787917"}, + {file = "matplotlib-3.8.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5864bdd7da445e4e5e011b199bb67168cdad10b501750367c496420f2ad00843"}, + {file = "matplotlib-3.8.2-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:ef8345b48e95cee45ff25192ed1f4857273117917a4dcd48e3905619bcd9c9b8"}, + {file = "matplotlib-3.8.2-cp39-cp39-win_amd64.whl", hash = "sha256:7c48d9e221b637c017232e3760ed30b4e8d5dfd081daf327e829bf2a72c731b4"}, + {file = "matplotlib-3.8.2-pp39-pypy39_pp73-macosx_10_12_x86_64.whl", hash = "sha256:aa11b3c6928a1e496c1a79917d51d4cd5d04f8a2e75f21df4949eeefdf697f4b"}, + {file = "matplotlib-3.8.2-pp39-pypy39_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d1095fecf99eeb7384dabad4bf44b965f929a5f6079654b681193edf7169ec20"}, + {file = "matplotlib-3.8.2-pp39-pypy39_pp73-win_amd64.whl", hash = "sha256:bddfb1db89bfaa855912261c805bd0e10218923cc262b9159a49c29a7a1c1afa"}, + {file = "matplotlib-3.8.2.tar.gz", hash = "sha256:01a978b871b881ee76017152f1f1a0cbf6bd5f7b8ff8c96df0df1bd57d8755a1"}, +] + [[package]] name = "mccabe" version = "0.7.0" @@ -416,7 +647,6 @@ dependencies = [ "click>=7.0", "colorama>=0.4; platform_system == \"Windows\"", "ghp-import>=1.0", - "importlib-metadata>=4.3; python_version < \"3.10\"", "jinja2>=2.11.1", "markdown>=3.2.1", "markupsafe>=2.0.1", @@ -490,18 +720,55 @@ dependencies = [ "Markdown>=3.3", "MarkupSafe>=1.1", "click>=7.0", - "importlib-metadata>=4.6; python_version < \"3.10\"", "mkdocs-autorefs>=0.3.1", "mkdocs>=1.4", "platformdirs>=2.2.0", "pymdown-extensions>=6.3", - "typing-extensions>=4.1; python_version < \"3.10\"", ] files = [ {file = "mkdocstrings-0.24.0-py3-none-any.whl", hash = "sha256:f4908560c10f587326d8f5165d1908817b2e280bbf707607f601c996366a2264"}, {file = "mkdocstrings-0.24.0.tar.gz", hash = "sha256:222b1165be41257b494a9d29b14135d2b7ca43f38161d5b10caae03b87bd4f7e"}, ] +[[package]] +name = "mkdocstrings-python" +version = "1.7.5" +requires_python = ">=3.8" +summary = "A Python handler for mkdocstrings." +dependencies = [ + "griffe>=0.37", + "mkdocstrings>=0.20", +] +files = [ + {file = "mkdocstrings_python-1.7.5-py3-none-any.whl", hash = "sha256:5f6246026353f0c0785135db70c3fe9a5d9318990fc7ceb11d62097b8ffdd704"}, + {file = "mkdocstrings_python-1.7.5.tar.gz", hash = "sha256:c7d143728257dbf1aa550446555a554b760dcd40a763f077189d298502b800be"}, +] + +[[package]] +name = "mkdocstrings" +version = "0.24.0" +extras = ["python"] +requires_python = ">=3.8" +summary = "Automatic documentation from sources, for MkDocs." +dependencies = [ + "mkdocstrings-python>=0.5.2", + "mkdocstrings==0.24.0", +] +files = [ + {file = "mkdocstrings-0.24.0-py3-none-any.whl", hash = "sha256:f4908560c10f587326d8f5165d1908817b2e280bbf707607f601c996366a2264"}, + {file = "mkdocstrings-0.24.0.tar.gz", hash = "sha256:222b1165be41257b494a9d29b14135d2b7ca43f38161d5b10caae03b87bd4f7e"}, +] + +[[package]] +name = "networkx" +version = "3.2.1" +requires_python = ">=3.9" +summary = "Python package for creating and manipulating graphs and networks" +files = [ + {file = "networkx-3.2.1-py3-none-any.whl", hash = "sha256:f18c69adc97877c42332c170849c96cefa91881c99a7cb3e95b7c659ebdc1ec2"}, + {file = "networkx-3.2.1.tar.gz", hash = "sha256:9f1bb5cf3409bf324e0a722c20bdb4c20ee39bf1c30ce8ae499c8502b0b5e0c6"}, +] + [[package]] name = "numpy" version = "1.26.2" @@ -574,6 +841,59 @@ files = [ {file = "pathspec-0.11.2.tar.gz", hash = "sha256:e0d8d0ac2f12da61956eb2306b69f9469b42f4deb0f3cb6ed47b9cce9996ced3"}, ] +[[package]] +name = "pillow" +version = "10.1.0" +requires_python = ">=3.8" +summary = "Python Imaging Library (Fork)" +files = [ + {file = "Pillow-10.1.0-cp310-cp310-macosx_10_10_x86_64.whl", hash = "sha256:1ab05f3db77e98f93964697c8efc49c7954b08dd61cff526b7f2531a22410106"}, + {file = "Pillow-10.1.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:6932a7652464746fcb484f7fc3618e6503d2066d853f68a4bd97193a3996e273"}, + {file = "Pillow-10.1.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a5f63b5a68daedc54c7c3464508d8c12075e56dcfbd42f8c1bf40169061ae666"}, + {file = "Pillow-10.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:c0949b55eb607898e28eaccb525ab104b2d86542a85c74baf3a6dc24002edec2"}, + {file = "Pillow-10.1.0-cp310-cp310-manylinux_2_28_aarch64.whl", hash = "sha256:ae88931f93214777c7a3aa0a8f92a683f83ecde27f65a45f95f22d289a69e593"}, + {file = "Pillow-10.1.0-cp310-cp310-manylinux_2_28_x86_64.whl", hash = "sha256:b0eb01ca85b2361b09480784a7931fc648ed8b7836f01fb9241141b968feb1db"}, + {file = "Pillow-10.1.0-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:d27b5997bdd2eb9fb199982bb7eb6164db0426904020dc38c10203187ae2ff2f"}, + {file = "Pillow-10.1.0-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:7df5608bc38bd37ef585ae9c38c9cd46d7c81498f086915b0f97255ea60c2818"}, + {file = "Pillow-10.1.0-cp310-cp310-win_amd64.whl", hash = "sha256:41f67248d92a5e0a2076d3517d8d4b1e41a97e2df10eb8f93106c89107f38b57"}, + {file = "Pillow-10.1.0-cp311-cp311-macosx_10_10_x86_64.whl", hash = "sha256:1fb29c07478e6c06a46b867e43b0bcdb241b44cc52be9bc25ce5944eed4648e7"}, + {file = "Pillow-10.1.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:2cdc65a46e74514ce742c2013cd4a2d12e8553e3a2563c64879f7c7e4d28bce7"}, + {file = "Pillow-10.1.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:50d08cd0a2ecd2a8657bd3d82c71efd5a58edb04d9308185d66c3a5a5bed9610"}, + {file = "Pillow-10.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:062a1610e3bc258bff2328ec43f34244fcec972ee0717200cb1425214fe5b839"}, + {file = "Pillow-10.1.0-cp311-cp311-manylinux_2_28_aarch64.whl", hash = "sha256:61f1a9d247317fa08a308daaa8ee7b3f760ab1809ca2da14ecc88ae4257d6172"}, + {file = "Pillow-10.1.0-cp311-cp311-manylinux_2_28_x86_64.whl", hash = "sha256:a646e48de237d860c36e0db37ecaecaa3619e6f3e9d5319e527ccbc8151df061"}, + {file = "Pillow-10.1.0-cp311-cp311-musllinux_1_1_aarch64.whl", hash = "sha256:47e5bf85b80abc03be7455c95b6d6e4896a62f6541c1f2ce77a7d2bb832af262"}, + {file = "Pillow-10.1.0-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:a92386125e9ee90381c3369f57a2a50fa9e6aa8b1cf1d9c4b200d41a7dd8e992"}, + {file = "Pillow-10.1.0-cp311-cp311-win_amd64.whl", hash = "sha256:0f7c276c05a9767e877a0b4c5050c8bee6a6d960d7f0c11ebda6b99746068c2a"}, + {file = "Pillow-10.1.0-cp312-cp312-macosx_10_10_x86_64.whl", hash = "sha256:a89b8312d51715b510a4fe9fc13686283f376cfd5abca8cd1c65e4c76e21081b"}, + {file = "Pillow-10.1.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:00f438bb841382b15d7deb9a05cc946ee0f2c352653c7aa659e75e592f6fa17d"}, + {file = "Pillow-10.1.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:3d929a19f5469b3f4df33a3df2983db070ebb2088a1e145e18facbc28cae5b27"}, + {file = "Pillow-10.1.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9a92109192b360634a4489c0c756364c0c3a2992906752165ecb50544c251312"}, + {file = "Pillow-10.1.0-cp312-cp312-manylinux_2_28_aarch64.whl", hash = "sha256:0248f86b3ea061e67817c47ecbe82c23f9dd5d5226200eb9090b3873d3ca32de"}, + {file = "Pillow-10.1.0-cp312-cp312-manylinux_2_28_x86_64.whl", hash = "sha256:9882a7451c680c12f232a422730f986a1fcd808da0fd428f08b671237237d651"}, + {file = "Pillow-10.1.0-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:1c3ac5423c8c1da5928aa12c6e258921956757d976405e9467c5f39d1d577a4b"}, + {file = "Pillow-10.1.0-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:806abdd8249ba3953c33742506fe414880bad78ac25cc9a9b1c6ae97bedd573f"}, + {file = "Pillow-10.1.0-cp312-cp312-win_amd64.whl", hash = "sha256:eaed6977fa73408b7b8a24e8b14e59e1668cfc0f4c40193ea7ced8e210adf996"}, + {file = "Pillow-10.1.0-cp39-cp39-macosx_10_10_x86_64.whl", hash = "sha256:0a026c188be3b443916179f5d04548092e253beb0c3e2ee0a4e2cdad72f66099"}, + {file = "Pillow-10.1.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:04f6f6149f266a100374ca3cc368b67fb27c4af9f1cc8cb6306d849dcdf12616"}, + {file = "Pillow-10.1.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:bb40c011447712d2e19cc261c82655f75f32cb724788df315ed992a4d65696bb"}, + {file = "Pillow-10.1.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1a8413794b4ad9719346cd9306118450b7b00d9a15846451549314a58ac42219"}, + {file = "Pillow-10.1.0-cp39-cp39-manylinux_2_28_aarch64.whl", hash = "sha256:c9aeea7b63edb7884b031a35305629a7593272b54f429a9869a4f63a1bf04c34"}, + {file = "Pillow-10.1.0-cp39-cp39-manylinux_2_28_x86_64.whl", hash = "sha256:b4005fee46ed9be0b8fb42be0c20e79411533d1fd58edabebc0dd24626882cfd"}, + {file = "Pillow-10.1.0-cp39-cp39-musllinux_1_1_aarch64.whl", hash = "sha256:4d0152565c6aa6ebbfb1e5d8624140a440f2b99bf7afaafbdbf6430426497f28"}, + {file = "Pillow-10.1.0-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:d921bc90b1defa55c9917ca6b6b71430e4286fc9e44c55ead78ca1a9f9eba5f2"}, + {file = "Pillow-10.1.0-cp39-cp39-win_amd64.whl", hash = "sha256:cfe96560c6ce2f4c07d6647af2d0f3c54cc33289894ebd88cfbb3bcd5391e256"}, + {file = "Pillow-10.1.0-pp310-pypy310_pp73-macosx_10_10_x86_64.whl", hash = "sha256:937bdc5a7f5343d1c97dc98149a0be7eb9704e937fe3dc7140e229ae4fc572a7"}, + {file = "Pillow-10.1.0-pp310-pypy310_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b1c25762197144e211efb5f4e8ad656f36c8d214d390585d1d21281f46d556ba"}, + {file = "Pillow-10.1.0-pp310-pypy310_pp73-manylinux_2_28_x86_64.whl", hash = "sha256:afc8eef765d948543a4775f00b7b8c079b3321d6b675dde0d02afa2ee23000b4"}, + {file = "Pillow-10.1.0-pp310-pypy310_pp73-win_amd64.whl", hash = "sha256:883f216eac8712b83a63f41b76ddfb7b2afab1b74abbb413c5df6680f071a6b9"}, + {file = "Pillow-10.1.0-pp39-pypy39_pp73-macosx_10_10_x86_64.whl", hash = "sha256:b920e4d028f6442bea9a75b7491c063f0b9a3972520731ed26c83e254302eb1e"}, + {file = "Pillow-10.1.0-pp39-pypy39_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1c41d960babf951e01a49c9746f92c5a7e0d939d1652d7ba30f6b3090f27e412"}, + {file = "Pillow-10.1.0-pp39-pypy39_pp73-manylinux_2_28_x86_64.whl", hash = "sha256:1fafabe50a6977ac70dfe829b2d5735fd54e190ab55259ec8aea4aaea412fa0b"}, + {file = "Pillow-10.1.0-pp39-pypy39_pp73-win_amd64.whl", hash = "sha256:3b834f4b16173e5b92ab6566f0473bfb09f939ba14b23b8da1f54fa63e4b623f"}, + {file = "Pillow-10.1.0.tar.gz", hash = "sha256:e6bf8de6c36ed96c86ea3b6e1d5273c53f46ef518a062464cd7ef5dd2cf92e38"}, +] + [[package]] name = "platformdirs" version = "4.1.0" @@ -612,15 +932,12 @@ summary = "python code static checker" dependencies = [ "astroid<=3.1.0-dev0,>=3.0.1", "colorama>=0.4.5; sys_platform == \"win32\"", - "dill>=0.2; python_version < \"3.11\"", "dill>=0.3.6; python_version >= \"3.11\"", "dill>=0.3.7; python_version >= \"3.12\"", "isort<6,>=4.2.5", "mccabe<0.8,>=0.6", "platformdirs>=2.2.0", - "tomli>=1.1.0; python_version < \"3.11\"", "tomlkit>=0.10.1", - "typing-extensions>=3.10.0; python_version < \"3.10\"", ] files = [ {file = "pylint-3.0.2-py3-none-any.whl", hash = "sha256:60ed5f3a9ff8b61839ff0348b3624ceeb9e6c2a92c514d81c9cc273da3b6bcda"}, @@ -641,6 +958,16 @@ files = [ {file = "pymdown_extensions-10.5.tar.gz", hash = "sha256:1b60f1e462adbec5a1ed79dac91f666c9c0d241fa294de1989f29d20096cfd0b"}, ] +[[package]] +name = "pyparsing" +version = "3.1.1" +requires_python = ">=3.6.8" +summary = "pyparsing module - Classes and methods to define and execute parsing grammars" +files = [ + {file = "pyparsing-3.1.1-py3-none-any.whl", hash = "sha256:32c7c0b711493c72ff18a981d24f28aaf9c1fb7ed5e9667c9e84e3db623bdbfb"}, + {file = "pyparsing-3.1.1.tar.gz", hash = "sha256:ede28a1a32462f5a9705e07aea48001a08f7cf81a021585011deba701581a0db"}, +] + [[package]] name = "pytest" version = "7.4.3" @@ -648,11 +975,9 @@ requires_python = ">=3.7" summary = "pytest: simple powerful testing with Python" dependencies = [ "colorama; sys_platform == \"win32\"", - "exceptiongroup>=1.0.0rc8; python_version < \"3.11\"", "iniconfig", "packaging", "pluggy<2.0,>=0.12", - "tomli>=1.0.0; python_version < \"3.11\"", ] files = [ {file = "pytest-7.4.3-py3-none-any.whl", hash = "sha256:0d009c083ea859a71b76adf7c1d502e4bc170b80a8ef002da5806527b9591fac"}, @@ -813,6 +1138,44 @@ files = [ {file = "requests-2.31.0.tar.gz", hash = "sha256:942c5a758f98d790eaed1a29cb6eefc7ffb0d1cf7af05c3d2791656dbd6ad1e1"}, ] +[[package]] +name = "scikit-learn" +version = "1.3.2" +requires_python = ">=3.8" +summary = "A set of python modules for machine learning and data mining" +dependencies = [ + "joblib>=1.1.1", + "numpy<2.0,>=1.17.3", + "scipy>=1.5.0", + "threadpoolctl>=2.0.0", +] +files = [ + {file = "scikit-learn-1.3.2.tar.gz", hash = "sha256:a2f54c76accc15a34bfb9066e6c7a56c1e7235dda5762b990792330b52ccfb05"}, + {file = "scikit_learn-1.3.2-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:8db94cd8a2e038b37a80a04df8783e09caac77cbe052146432e67800e430c028"}, + {file = "scikit_learn-1.3.2-cp312-cp312-macosx_12_0_arm64.whl", hash = "sha256:61a6efd384258789aa89415a410dcdb39a50e19d3d8410bd29be365bcdd512d5"}, + {file = "scikit_learn-1.3.2-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:cb06f8dce3f5ddc5dee1715a9b9f19f20d295bed8e3cd4fa51e1d050347de525"}, + {file = "scikit_learn-1.3.2-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5b2de18d86f630d68fe1f87af690d451388bb186480afc719e5f770590c2ef6c"}, + {file = "scikit_learn-1.3.2-cp312-cp312-win_amd64.whl", hash = "sha256:0402638c9a7c219ee52c94cbebc8fcb5eb9fe9c773717965c1f4185588ad3107"}, +] + +[[package]] +name = "scipy" +version = "1.11.4" +requires_python = ">=3.9" +summary = "Fundamental algorithms for scientific computing in Python" +dependencies = [ + "numpy<1.28.0,>=1.21.6", +] +files = [ + {file = "scipy-1.11.4-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:028eccd22e654b3ea01ee63705681ee79933652b2d8f873e7949898dda6d11b6"}, + {file = "scipy-1.11.4-cp312-cp312-macosx_12_0_arm64.whl", hash = "sha256:2c6ff6ef9cc27f9b3db93a6f8b38f97387e6e0591600369a297a50a8e96e835d"}, + {file = "scipy-1.11.4-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b030c6674b9230d37c5c60ab456e2cf12f6784596d15ce8da9365e70896effc4"}, + {file = "scipy-1.11.4-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ad669df80528aeca5f557712102538f4f37e503f0c5b9541655016dd0932ca79"}, + {file = "scipy-1.11.4-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:ce7fff2e23ab2cc81ff452a9444c215c28e6305f396b2ba88343a567feec9660"}, + {file = "scipy-1.11.4-cp312-cp312-win_amd64.whl", hash = "sha256:36750b7733d960d7994888f0d148d31ea3017ac15eef664194b4ef68d36a4a97"}, + {file = "scipy-1.11.4.tar.gz", hash = "sha256:90a2b78e7f5733b9de748f589f09225013685f9b218275257f8a8168ededaeaa"}, +] + [[package]] name = "setuptools" version = "69.0.2" @@ -844,13 +1207,13 @@ files = [ ] [[package]] -name = "tomli" -version = "2.0.1" -requires_python = ">=3.7" -summary = "A lil' TOML parser" +name = "threadpoolctl" +version = "3.2.0" +requires_python = ">=3.8" +summary = "threadpoolctl" files = [ - {file = "tomli-2.0.1-py3-none-any.whl", hash = "sha256:939de3e7a6161af0c887ef91b7d41a53e7c5a1ca976325f429cb46ea9bc30ecc"}, - {file = "tomli-2.0.1.tar.gz", hash = "sha256:de526c12914f0c550d15924c62d72abc48d6fe7364aa87328337a31007fe8a4f"}, + {file = "threadpoolctl-3.2.0-py3-none-any.whl", hash = "sha256:2b7818516e423bdaebb97c723f86a7c6b0a83d3f3b0970328d66f4d9104dc032"}, + {file = "threadpoolctl-3.2.0.tar.gz", hash = "sha256:c96a0ba3bdddeaca37dc4cc7344aafad41cdb8c313f74fdfe387a867bba93355"}, ] [[package]] @@ -863,16 +1226,6 @@ files = [ {file = "tomlkit-0.12.3.tar.gz", hash = "sha256:75baf5012d06501f07bee5bf8e801b9f343e7aac5a92581f20f80ce632e6b5a4"}, ] -[[package]] -name = "typing-extensions" -version = "4.8.0" -requires_python = ">=3.8" -summary = "Backported and Experimental Type Hints for Python 3.8+" -files = [ - {file = "typing_extensions-4.8.0-py3-none-any.whl", hash = "sha256:8f92fc8806f9a6b641eaa5318da32b44d401efaac0f6678c9bc448ba3605faa0"}, - {file = "typing_extensions-4.8.0.tar.gz", hash = "sha256:df8e4339e9cb77357558cbdbceca33c303714cf861d1eef15e1070055ae8b7ef"}, -] - [[package]] name = "urllib3" version = "2.1.0" @@ -910,13 +1263,3 @@ files = [ {file = "watchdog-3.0.0-py3-none-win_ia64.whl", hash = "sha256:5d9f3a10e02d7371cd929b5d8f11e87d4bad890212ed3901f9b4d68767bee759"}, {file = "watchdog-3.0.0.tar.gz", hash = "sha256:4d98a320595da7a7c5a18fc48cb633c2e73cda78f93cac2ef42d42bf609a33f9"}, ] - -[[package]] -name = "zipp" -version = "3.17.0" -requires_python = ">=3.8" -summary = "Backport of pathlib-compatible object wrapper for zip files" -files = [ - {file = "zipp-3.17.0-py3-none-any.whl", hash = "sha256:0e923e726174922dce09c53c59ad483ff7bbb8e572e00c7f7c46b88556409f31"}, - {file = "zipp-3.17.0.tar.gz", hash = "sha256:84e64a1c28cf7e91ed2078bb8cc8c259cb19b76942096c8d7b84947690cabaf0"}, -] diff --git a/pyproject.toml b/pyproject.toml index 7a0b2cb..05f1f8d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,8 +7,12 @@ authors = [ ] dependencies = [ "numpy>=1.26.2", + "matplotlib>=3.8.2", + "scipy>=1.11.4", + "scikit-learn>=1.3.2", + "networkx>=3.2.1", ] -requires-python = ">=3.9" +requires-python = ">=3.12" readme = "docs/README.md" license = {file = "LICENSE"} keywords = ["surrogates", "multidisciplinary", "multifidelity", "adaptive", "collocation", "metamodeling"] @@ -40,6 +44,11 @@ missing-class-docstring, missing-function-docstring ''' +[tool.pytest.ini_options] +filterwarnings = [ + "ignore::DeprecationWarning" +] + [tool.pdm.version] source = "file" path = "src/amisc/__init__.py" @@ -47,6 +56,7 @@ path = "src/amisc/__init__.py" [tool.pdm.scripts] release = "python release.py" test = "pytest --cov=amisc tests" +docs = "mkdocs serve" [tool.pdm.dev-dependencies] dev = [ @@ -55,6 +65,6 @@ dev = [ "pytest-cov>=4.1.0", "mkdocs>=1.5.3", "GitPython>=3.1.40", - "mkdocstrings>=0.24.0", "mkdocs-material>=9.5.0", + "mkdocstrings[python]>=0.24.0", ] diff --git a/src/amisc/__init__.py b/src/amisc/__init__.py index 7d228fa..af4d398 100644 --- a/src/amisc/__init__.py +++ b/src/amisc/__init__.py @@ -1,6 +1,17 @@ -"""Adaptive multi-index stochastic collocation (surrogates) for metamodeling of multidisciplinary systems. +"""Adaptive multi-index stochastic collocation for metamodeling/surrogates of multidisciplinary systems. - Author - Joshua Eckels (eckelsjd@umich.edu) - License - GNU GPLv3 """ +import numpy as np + +from amisc.interpolator import BaseInterpolator +from amisc.rv import BaseRV + __version__ = "0.0.8" + +# Custom types that are used frequently +IndexSet = list[tuple[tuple, tuple]] +MiscTree = dict[str: dict[str: float | BaseInterpolator]] +InterpResults = BaseInterpolator | tuple[list[int | tuple | str], np.ndarray, BaseInterpolator] +IndicesRV = list[int | str | BaseRV] | int | str | BaseRV diff --git a/src/amisc/component.py b/src/amisc/component.py new file mode 100644 index 0000000..2272ff8 --- /dev/null +++ b/src/amisc/component.py @@ -0,0 +1,869 @@ +"""`component.py` + +A Component is an `amisc` wrapper around a single discipline model. It manages surrogate construction and optionally +a hierarchy of modeling fidelities that may be available. Concrete component classes all inherit from the base +`ComponentSurrogate` class provided here. Components manage an array of `BaseInterpolator` objects to form a +multifidelity hierarchy. + +Includes +-------- +- `ComponentSurrogate`: the base class that is fundamental to the adaptive multi-index stochastic collocation strategy +- `SparseGridSurrogate`: an AMISC component that manages a hierarchy of `LagrangeInterpolator` objects +- `AnalyticalSurrogate`: a light wrapper around a single discipline model that does not require surrogate approximation +""" +import itertools +import copy +import ast +from pathlib import Path +from abc import ABC, abstractmethod +from concurrent.futures import Executor, ALL_COMPLETED, wait + +import numpy as np +from sklearn.linear_model import Ridge +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import MaxAbsScaler + +from amisc.utils import get_logger +from amisc.rv import BaseRV +from amisc import IndexSet, MiscTree, InterpResults +from amisc.interpolator import BaseInterpolator, LagrangeInterpolator + + +class ComponentSurrogate(ABC): + """The base multi-index stochastic collocation (MISC) surrogate class for a single discipline component model. + + !!! Info "Multi-indices" + A multi-index is a tuple of natural numbers, each specifying a level of fidelity. You will frequently see two + multi-indices: `alpha` and `beta`. The `alpha` (or $\\alpha$) indices specify physical model fidelity and get + passed to the model as an additional argument (e.g. things like discretization level, time step size, etc.). + The `beta` (or $\\beta$) indices specify surrogate refinement level, so typically an indication of the amount of + training data used. Each fidelity index in $\\alpha$ and $\\beta$ increase in refinement from $0$ up to + `max_alpha` and `max_beta`. From the surrogate's perspective, the concatenation of $(\\alpha, \\beta)$ fully + specifies a single fidelity "level". The `ComponentSurrogate` forms an approximation of the model by summing + up over many of these concatenated sets of $(\\alpha, \\beta)$. These lists are stored in a data structure of + `list[ tuple[ tuple, tuple ], ...]`. When $\\alpha$ or $\\beta$ are used as keys in a `dict`, they are cast to + a Python `str` from a `tuple`. + + :ivar index_set: the current active set of multi-indices in the MISC approximation + :ivar candidate_set: all neighboring multi-indices that are candidates for inclusion in `index_set` + :ivar x_vars: list of variables that define the input domain + :ivar ydim: the number of outputs returned by the model + :ivar _model: stores a ref to the model or function that is to be approximated, callable as `ret = model(x)` + :ivar _model_args: additional arguments to supply to the model + :ivar _model_kwargs: additional keyword arguments to supply to the model + :ivar truth_alpha: the model fidelity indices to treat as the "ground truth" model + :ivar max_refine: the maximum level of refinement for each fidelity index in $(\\alpha, \\beta)$ + :ivar surrogates: keeps track of the `BaseInterpolator` associated with each set of $(\\alpha, \\beta)$ + :ivar costs: keeps track of total cost associated with adding a single $(\\alpha, \\beta)$ to the MISC approximation + :ivar misc_coeff: the combination technique coefficients for the MISC approximation + + :vartype index_set: IndexSet + :vartype candidate_set: IndexSet + :vartype x_vars: list[BaseRV] + :vartype ydim: int + :vartype _model: callable[np.ndarray] -> dict + :vartype _model_args: tuple + :vartype _model_kwargs: dict + :vartype truth_alpha: tuple[int, ...] + :vartype max_refine: list[int, ...] + :vartype surrogates: MiscTree + :vartype costs: MiscTree + :vartype misc_coeff: MiscTree + """ + + def __init__(self, x_vars: list[BaseRV] | BaseRV, model: callable, + multi_index: IndexSet = None, + truth_alpha: tuple = (), max_alpha: tuple = (), max_beta: tuple = (), + log_file: str | Path = None, executor: Executor = None, + model_args: tuple = (), model_kwargs: dict = None): + """Construct the MISC surrogate and initialize with any multi-indices passed in. + + !!! Info "Model specification" + The model is a callable function of the form `ret = model(x, *args, **kwargs)`. The return value is a + dictionary of the form `ret = {'y': y, 'files': files, 'cost': cost}`. In the return dictionary, you + specify the raw model output `y` as an `np.ndarray` at a _minimum_. Optionally, you can specify paths to + output files and the average model cost (in units of seconds of cpu time), and anything else you want. + + !!! Warning + If the model has multiple fidelities, then the function signature must be `model(x, alpha, *args, **kwargs)` + ; the first argument after `x` will always be the fidelity indices `alpha`. The rest of `model_args` will + be passed in after (you do not need to include `alpha` in `model_args`, it is done automatically). + + :param x_vars: `[X1, X2, ...]` list of variables specifying bounds/pdfs for each input + :param model: the function to approximate, callable as `ret = model(x, *args, **kwargs)` + :param multi_index: `[((alpha1), (beta1)), ... ]` list of concatenated multi-indices $(\\alpha, \\beta)$ + :param truth_alpha: specifies the highest model fidelity indices necessary for a "ground truth" comparison + :param max_alpha: the maximum model refinement indices to allow, defaults to `(2,...)` if applicable + :param max_beta: the maximum surrogate refinement indices, defaults to `(2,...)` of length `x_dim` + :param log_file: specifies a log file (optional) + :param executor: parallel executor used to add candidate indices in parallel (optional) + :param model_args: optional args to pass when calling the model + :param model_kwargs: optional kwargs to pass when calling the model + """ + self.logger = get_logger(self.__class__.__name__, log_file=log_file) + self.log_file = log_file + self.executor = executor + self.training_flag = None # Keep track of which MISC coeffs are active + # (True=active set, False=active+candidate sets, None=Neither/unknown) + + multi_index = list() if multi_index is None else multi_index + assert self.is_downward_closed(multi_index), 'Must be a downward closed set.' + self.ydim = None + self.index_set = [] # The active index set for the MISC approximation + self.candidate_set = [] # Candidate indices for refinement + self._model = model + self._model_args = model_args + self._model_kwargs = model_kwargs if model_kwargs is not None else {} + self.truth_alpha = truth_alpha + self.x_vars = x_vars if isinstance(x_vars, list) else [x_vars] + max_alpha = (2,)*len(truth_alpha) if max_alpha == () else max_alpha + max_beta = (2,)*len(self.x_vars) if max_beta == () else max_beta + self.max_refine = list(max_alpha + max_beta) # Max refinement indices + + # Initialize important tree-like structures + self.surrogates = dict() # Maps alphas -> betas -> surrogates + self.costs = dict() # Maps alphas -> betas -> wall clock run times + self.misc_coeff = dict() # Maps alphas -> betas -> MISC coefficients + + # Construct vectors of [0,1]^dim(alpha+beta) + Nij = len(self.max_refine) + self.ij = np.zeros((2 ** Nij, Nij), dtype=np.uint8) + for i, ele in enumerate(itertools.product([0, 1], repeat=Nij)): + self.ij[i, :] = ele + + # Initialize any indices that were passed in + multi_index = list() if multi_index is None else multi_index + for alpha, beta in multi_index: + self.activate_index(alpha, beta) + + def activate_index(self, alpha: tuple, beta: tuple): + """Add a multi-index to the active set and all neighbors to the candidate set. + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + """ + # User is responsible for making sure index set is downward-closed + self.add_surrogate(alpha, beta) + ele = (alpha, beta) + if ele in self.index_set: + self.logger.warning(f'Multi-index {ele} is already in the active index set. Ignoring...') + return + + # Add all possible new candidates (distance of one unit vector away) + ind = list(alpha + beta) + new_candidates = [] + for i in range(len(ind)): + ind_new = ind.copy() + ind_new[i] += 1 + + # Don't add if we surpass a refinement limit + if np.any(np.array(ind_new) > np.array(self.max_refine)): + continue + + # Add the new index if it maintains downward-closedness + new_cand = (tuple(ind_new[:len(alpha)]), tuple(ind_new[len(alpha):])) + down_closed = True + for j in range(len(ind)): + ind_check = ind_new.copy() + ind_check[j] -= 1 + if ind_check[j] >= 0: + tup_check = (tuple(ind_check[:len(alpha)]), tuple(ind_check[len(alpha):])) + if tup_check not in self.index_set and tup_check != ele: + down_closed = False + break + if down_closed: + new_candidates.append(new_cand) + + # Build an interpolator for each new candidate + if self.executor is None: # Sequential + for a, b in new_candidates: + self.add_surrogate(a, b) + else: # Parallel + temp_exc = self.executor + self.executor = None + for a, b in new_candidates: + if str(a) not in self.surrogates: + self.surrogates[str(a)] = dict() + self.costs[str(a)] = dict() + self.misc_coeff[str(a)] = dict() + self.parallel_add_candidates(new_candidates, temp_exc) + self.executor = temp_exc + + # Move to the active index set + if ele in self.candidate_set: + self.candidate_set.remove(ele) + self.index_set.append(ele) + new_candidates = [cand for cand in new_candidates if cand not in self.candidate_set] + self.candidate_set.extend(new_candidates) + self.training_flag = None # Makes sure misc coeffs get recomputed next time + + def add_surrogate(self, alpha: tuple, beta: tuple): + """Build a `BaseInterpolator` object for a given $(\\alpha, \\beta)$ + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + """ + # Create a dictionary for each alpha model to store multiple surrogate fidelities (beta) + if str(alpha) not in self.surrogates: + self.surrogates[str(alpha)] = dict() + self.costs[str(alpha)] = dict() + self.misc_coeff[str(alpha)] = dict() + + # Create a new interpolator object for this multi-index (abstract method) + if self.surrogates[str(alpha)].get(str(beta), None) is None: + self.logger.info(f'Building interpolator for index {(alpha, beta)} ...') + x_new_idx, x_new, interp = self.build_interpolator(alpha, beta) + self.surrogates[str(alpha)][str(beta)] = interp + cost = self.update_interpolator(x_new_idx, x_new, interp) # Awkward, but needed to separate the model evals + self.costs[str(alpha)][str(beta)] = cost + if self.ydim is None: + self.ydim = interp.ydim() + + def init_coarse(self): + """Initialize the coarsest interpolation and add to the active index set""" + alpha = (0,) * len(self.truth_alpha) + beta = (0,) * len(self.max_refine[len(self.truth_alpha):]) + self.activate_index(alpha, beta) + + def iterate_candidates(self): + """Iterate candidate indices one by one into the active index set. + + :yields alpha, beta: the multi-indices of the current candidate that has been moved to active set + """ + for alpha, beta in list(self.candidate_set): + # Temporarily add a candidate index to active set + self.index_set.append((alpha, beta)) + yield alpha, beta + del self.index_set[-1] + + def predict(self, x: np.ndarray | float, use_model: str | tuple = None, model_dir: str | Path = None, + training: bool = False, index_set: IndexSet = None) -> np.ndarray: + """Evaluate the MISC approximation at new points `x`. + + !!! Note + By default this will predict the MISC surrogate approximation. However, for convenience you can also specify + `use_model` to call the underlying function instead. + + :param x: `(..., x_dim)` the points to be interpolated, must be within input domain for accuracy + :param use_model: 'best'=high-fidelity, 'worst'=low-fidelity, tuple=a specific `alpha`, None=surrogate (default) + :param model_dir: directory to save output files if `use_model` is specified, ignored otherwise + :param training: if `True`, then only compute with the active index set, otherwise use all candidates as well + :param index_set: a list of concatenated $(\\alpha, \\beta)$ to override `self.index_set` if given, else ignore + :returns y: `(..., y_dim)` the surrogate approximation of the function (or the function itself if `use_model`) + """ + if use_model is not None: + return self._bypass_surrogate(x, use_model, model_dir) + + index_set, misc_coeff = self._combination(index_set, training) # Choose the correct index set and misc_coeff + + y = np.zeros(x.shape[:-1] + (self.ydim,)) + for alpha, beta in index_set: + comb_coeff = misc_coeff[str(alpha)][str(beta)] + if np.abs(comb_coeff) > 0: + func = self.surrogates[str(alpha)][str(beta)] + y += comb_coeff * func(x) + + return y + + def __call__(self, *args, **kwargs): + """Here for convenience so you can also do `ret = surrogate(x)`, just like the `BaseInterpolator`.""" + return self.predict(*args, **kwargs) + + def update_misc_coeffs(self, index_set: IndexSet = None) -> MiscTree: + """Update the combination technique coeffs for MISC using the given index set. + + :param index_set: the index set to consider when computing the MISC coefficients, defaults to the active set + :returns: the MISC coefficients for the given index set ($\\alpha$ -> $\\beta$ -> coeff) + """ + if index_set is None: + index_set = self.index_set + + # Construct a (N_indices, dim(alpha+beta)) refactor of the index_set for arrayed computations + index_mat = np.zeros((len(index_set), len(self.max_refine)), dtype=np.uint8) + for i, (alpha, beta) in enumerate(index_set): + index_mat[i, :] = alpha + beta + index_mat = np.expand_dims(index_mat, axis=0) # (1, Ns, Nij) + + misc_coeff = dict() + for alpha, beta in index_set: + # Add permutations of [0, 1] to (alpha, beta) + alpha_beta = np.array(alpha+beta, dtype=np.uint8)[np.newaxis, :] # (1, Nij) + new_indices = np.expand_dims(alpha_beta + self.ij, axis=1) # (2**Nij, 1, Nij) + + # Find which indices are in the index_set (using np broadcasting comparison) + diff = new_indices - index_mat # (2**Nij, Ns, Nij) + idx = np.count_nonzero(diff, axis=-1) == 0 # (2**Nij, Ns) + idx = np.any(idx, axis=-1) # (2**Nij,) + ij_use = self.ij[idx, :] # (*, Nij) + l1_norm = np.sum(np.abs(ij_use), axis=-1) # (*,) + coeff = np.sum((-1) ** l1_norm) # float + + # Save misc coeff to a dict() tree structure + if misc_coeff.get(str(alpha)) is None: + misc_coeff[str(alpha)] = dict() + misc_coeff[str(alpha)][str(beta)] = coeff + self.misc_coeff[str(alpha)][str(beta)] = coeff + + return misc_coeff + + def get_sub_surrogate(self, alpha: tuple, beta: tuple) -> BaseInterpolator: + """Get the specific sub-surrogate corresponding to the $(\\alpha, \\beta)$ fidelity. + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + :returns: the corresponding `BaseInterpolator` object + """ + return self.surrogates[str(alpha)][str(beta)] + + def get_cost(self, alpha: tuple, beta: tuple) -> float: + """Return the total cost (wall time s) required to add $(\\alpha, \\beta)$ to the MISC approximation. + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + """ + try: + return self.costs[str(alpha)][str(beta)] + except: + return 0.0 + + def update_input_bds(self, idx: int, bds: tuple): + """Update the bounds of the input variable at the given index. + + :param idx: the index of the input variable to update + :param bds: the new bounds + """ + self.x_vars[int(idx)].update_bounds(*bds) + + # Update the bounds in all associated surrogates + for alpha in self.surrogates: + for beta in self.surrogates[alpha]: + self.surrogates[alpha][beta].update_input_bds(idx, bds) + + def save_enabled(self): + """Return whether this model wants to save outputs to file. + + !!! Note + You can specify that a model wants to save outputs to file by providing an `'output_dir'` kwarg. + """ + return self._model_kwargs.get('output_dir') is not None + + def _set_output_dir(self, output_dir: str | Path): + """Update the component model output directory. + + :param output_dir: the new directory for model output files + """ + if output_dir is not None: + output_dir = str(Path(output_dir).resolve()) + self._model_kwargs['output_dir'] = output_dir + for alpha in self.surrogates: + for beta in self.surrogates[alpha]: + self.surrogates[alpha][beta]._model_kwargs['output_dir'] = output_dir + + def __repr__(self): + """Shows all multi-indices in the current approximation and their corresponding MISC coefficients.""" + s = f'Inputs \u2014 {[str(var) for var in self.x_vars]}\n' + if self.training_flag is None: + self.update_misc_coeffs() + self.training_flag = True + + if self.training_flag: + s += '(Training mode)\n' + for alpha, beta in self.index_set: + s += f"[{int(self.misc_coeff[str(alpha)][str(beta)])}] \u2014 {alpha}, {beta}\n" + for alpha, beta in self.candidate_set: + s += f"[-] \u2014 {alpha}, {beta}\n" + else: + s += '(Evaluation mode)\n' + for alpha, beta in self.index_set + self.candidate_set: + s += f"[{int(self.misc_coeff[str(alpha)][str(beta)])}] \u2014 {alpha}, {beta}\n" + return s + + def __str__(self): + """Everyone will view these objects the same way.""" + return self.__repr__() + + def _bypass_surrogate(self, x, use_model, model_dir): + """Bypass surrogate evaluation and use the specified model""" + output_dir = self._model_kwargs.get('output_dir') + if self.save_enabled(): + self._model_kwargs['output_dir'] = model_dir + + alpha_use = {'best': self.truth_alpha, 'worst': (0,) * len(self.truth_alpha)}.get(use_model, use_model) + kwargs = copy.deepcopy(self._model_kwargs) + if len(alpha_use) > 0: + kwargs['alpha'] = alpha_use + ret = self._model(x, *self._model_args, **kwargs) + + if output_dir is not None: + self._model_kwargs['output_dir'] = output_dir + + if not isinstance(ret, dict): + self.logger.warning(f"Function {self._model} did not return a dict of the form {{'y': y}}. Please make sure" + f" you do so to avoid conflicts. Returning the value directly instead...") + + return ret['y'] if isinstance(ret, dict) else ret + + def _combination(self, index_set, training): + """Decide which index set and corresponding misc coefficients to use.""" + misc_coeff = copy.deepcopy(self.misc_coeff) + if index_set is None: + # Use active indices + candidate indices depending on training mode + index_set = self.index_set if training else self.index_set + self.candidate_set + + # Decide when to update misc coefficients + if self.training_flag is None: + misc_coeff = self.update_misc_coeffs(index_set) # On initialization or reset + else: + if (not self.training_flag and training) or (self.training_flag and not training): + misc_coeff = self.update_misc_coeffs(index_set) # Logical XOR cases for training mode + + # Save an indication of what state the MISC coefficients are in (i.e. training or eval mode) + self.training_flag = training + else: + # If we passed in an index set, always recompute misc coeff and toggle for reset on next call + misc_coeff = self.update_misc_coeffs(index_set) + self.training_flag = None + + return index_set, misc_coeff + + @staticmethod + def is_one_level_refinement(beta_old: tuple, beta_new: tuple) -> bool: + """Check if a new `beta` multi-index is a one-level refinement from a previous `beta`. + + !!! Example + Refining from `(0, 1, 2)` to the new multi-index `(1, 1, 2)` is a one-level refinement. But refining to + either `(2, 1, 2)` or `(1, 2, 2)` are not, since more than one refinement occurs at the same time. + + :param beta_old: the starting multi-index + :param beta_new: the new refined multi-index + :returns: whether `beta_new` is a one-level refinement from `beta_old` + """ + level_diff = np.array(beta_new, dtype=int) - np.array(beta_old, dtype=int) + ind = np.nonzero(level_diff)[0] + return ind.shape[0] == 1 and level_diff[ind] == 1 + + @staticmethod + def is_downward_closed(indices: IndexSet) -> bool: + """Return if a list of $(\\alpha, \\beta)$ multi-indices is downward-closed. + + MISC approximations require a downward-closed set in order to use the combination-technique formula for the + coefficients (as implemented here). + + !!! Example + The list `[( (0,), (0,) ), ( (1,), (0,) ), ( (1,), (1,) )]` is downward-closed. You can visualize this as + building a stack of cubes: in order to place a cube, all adjacent cubes must be present (does the logo + make sense now?). + + :param indices: list() of (`alpha`, `beta`) multi-indices + :returns: whether the set of indices is downward-closed + """ + # Iterate over every multi-index + for alpha, beta in indices: + # Every smaller multi-index must also be included in the indices list + sub_sets = [np.arange(tuple(alpha+beta)[i]+1) for i in range(len(alpha) + len(beta))] + for ele in itertools.product(*sub_sets): + tup = (tuple(ele[:len(alpha)]), tuple(ele[len(alpha):])) + if tup not in indices: + return False + return True + + @abstractmethod + def build_interpolator(self, alpha: tuple, beta: tuple) -> InterpResults: + """Return a `BaseInterpolator` object and new refinement points for a given $(\\alpha, \\beta)$ multi-index. + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + :returns: `idx`, `x`, `interp` - list of new grid indices, the new grid points `(N_new, x_dim)`, and the + `BaseInterpolator` object. Similar to `BaseInterpolator.refine()`. + """ + pass + + @abstractmethod + def update_interpolator(self, x_new_idx: list[int | tuple | str], x_new: np.ndarray, interp: BaseInterpolator) -> float: + """Secondary method to actually compute and save model evaluations within the interpolator. + + !!! Note + This distinction with `build_interpolator` was necessary to separately construct the interpolator and be + able to evaluate the model at the new interpolation points. You can see that `parallel_add_candidates` + uses this distinction to compute the model in parallel on MPI workers, for example. + + :param x_new_idx: list of new grid point indices + :param x_new: `(N_new, x_dim)`, the new grid point locations + :param interp: the `BaseInterpolator` object to compute model evaluations with + :returns cost: the cost (in wall time seconds) required to add this `BaseInterpolator` object + """ + pass + + @abstractmethod + def parallel_add_candidates(self, candidates: IndexSet, executor: Executor): + """Defines a function to handle adding candidate indices in parallel. + + !!! Note + While `build_interpolator` can make changes to 'self', these changes will not be saved in the master task + if running in parallel over MPI workers, for example. This method is a workaround so that all required + mutable changes to 'self' are made in the master task, before distributing tasks to parallel workers + using this method. You can pass if you don't plan to add candidates in parallel. + + :param candidates: list of [(alpha, beta),...] multi-indices + :param executor: the executor used to iterate candidates in parallel + """ + pass + + +class SparseGridSurrogate(ComponentSurrogate): + """Concrete MISC surrogate class that maintains a sparse grid composed of smaller tensor-product grids. + + !!! Note + MISC itself can be thought of as an extension to the well-known sparse grid technique, so this class + readily integrates with the MISC implementation in `ComponentSurrogate`. Sparse grids limit the curse + of dimensionality up to about `dim = 10-15` for the input space (which would otherwise be infeasible with a + normal full tensor-product grid of the same size). + + !!! Info "About points in a sparse grid" + A sparse grid approximates a full tensor-product grid $(N_1, N_2, ..., N_d)$, where $N_i$ is the number of grid + points along dimension $i$, for a $d$-dimensional space. Each point is uniquely identified in the sparse grid + by a list of indices $(j_1, j_2, ..., j_d)$, where $j_i = 0 ... N_i$. We refer to this unique identifier as a + "grid coordinate". In the `HashSG` data structure, we use a `str(tuple(coord))` representation to uniquely + identify the coordinate in a hash DS like Python's `dict`. + + :cvar HashSG: a type alias for the hash storage of the sparse grid data (a tree-like DS using dicts) + :ivar curr_max_beta: the current maximum $\\beta$ refinement indices in the sparse grid (for each $\\alpha$) + :ivar x_grids: maps $\\alpha$ indices to a list of 1d grids corresponding to `curr_max_beta` + :ivar xi_map: the sparse grid interpolation points + :ivar yi_map: the function values at all sparse grid points + :ivar yi_nan_map: imputed function values to use when `yi_map` contains `nan` data (sometimes the model fails...) + :ivar yi_files: optional filenames corresponding to the sparse grid `yi_map` data + + :vartype HashSG: dict[str: dict[str: np.ndarray | str]] + :vartype curr_max_beta: dict[str: list[int]] + :vartype x_grids: dict[str: np.ndarray] + :vartype xi_map: HashSG + :vartype yi_map: HashSG + :vartype yi_nan_map: HashSG + :vartype yi_files: HashSG + """ + + HashSG = dict[str: dict[str: np.ndarray | str]] + + def __init__(self, *args, **kwargs): + # Initialize tree-like hash structures for maintaining a sparse grid of smaller tensor-product grids + self.curr_max_beta = dict() # Maps alphas -> current max refinement indices + self.x_grids = dict() # Maps alphas -> list of ndarrays specifying 1d grids corresponding to max_beta + self.xi_map = dict() # Maps alphas -> grid point coords -> interpolation points + self.yi_map = dict() # Maps alphas -> grid point coords -> interpolation function values + self.yi_nan_map = dict() # Maps alphas -> grid point coords -> interpolated yi values when yi=nan + self.yi_files = dict() # Maps alphas -> grid point coords -> model output files (optional) + super().__init__(*args, **kwargs) + + # Override + def predict(self, x, use_model=None, model_dir=None, training=False, index_set=None): + """Need to override `super()` to allow passing in interpolation grids `xi` and `yi`.""" + if use_model is not None: + return self._bypass_surrogate(x, use_model, model_dir) + + index_set, misc_coeff = self._combination(index_set, training) + + y = np.zeros(x.shape[:-1] + (self.ydim,)) + for alpha, beta in index_set: + comb_coeff = misc_coeff[str(alpha)][str(beta)] + if np.abs(comb_coeff) > 0: + # Gather the xi/yi interpolation points/qoi_ind for this sub tensor-product grid + interp = self.surrogates[str(alpha)][str(beta)] + xi, yi = self.get_tensor_grid(alpha, beta) + + # Add this sub tensor-product grid to the MISC approximation + y += comb_coeff * interp(x, xi=xi, yi=yi) + + return y + + def get_tensor_grid(self, alpha: tuple, beta: tuple, update_nan: bool = True) -> tuple[np.ndarray, np.ndarray]: + """Construct the `xi/yi` sub tensor-product grids for a given $(\\alpha, \\beta)$ multi-index. + + :param alpha: model fidelity multi-index + :param beta: surrogate fidelity multi-index + :param update_nan: try to fill `nan` with imputed values, otherwise just return the `nans` in place + :returns: `xi, yi`, of size `(prod(grid_sizes), x_dim)` and `(prod(grid_sizes), y_dim)` respectively, the + interpolation grid points and corresponding function values for this tensor-product grid + """ + interp = self.surrogates[str(alpha)][str(beta)] + grid_sizes = interp.get_grid_sizes(beta) + coords = [np.arange(grid_sizes[n]) for n in range(interp.xdim())] + xi = np.zeros((np.prod(grid_sizes), interp.xdim()), dtype=np.float32) + yi = np.zeros((np.prod(grid_sizes), self.ydim), dtype=np.float32) + for i, coord in enumerate(itertools.product(*coords)): + xi[i, :] = self.xi_map[str(alpha)][str(coord)] + yi_curr = self.yi_map[str(alpha)][str(coord)] + if update_nan and np.any(np.isnan(yi_curr)): + # Try to replace NaN values if they are stored + yi_curr = self.yi_nan_map[str(alpha)].get(str(coord), yi_curr) + yi[i, :] = yi_curr + + return xi, yi + + def get_training_data(self) -> tuple[dict[str: np.ndarray], dict[str: np.ndarray]]: + """Grab all `x,y` training data stored in the sparse grid for each model fidelity level $\\alpha$. + + :returns: `xi`, `yi`, each a `dict` mapping `alpha` indices to `np.ndarrays` + """ + xi, yi = dict(), dict() + for alpha, x_map in self.xi_map.items(): + x = np.zeros((len(x_map), len(self.x_vars))) + y = np.zeros((len(x_map), self.ydim)) + for i, (coord, x_coord) in enumerate(x_map.items()): + x[i, :] = x_coord + y[i, :] = self.yi_nan_map[alpha].get(coord, self.yi_map[alpha][coord]) + + xi[alpha] = x + yi[alpha] = y + + return xi, yi + + def update_yi(self, alpha: tuple, beta: tuple, yi_dict: dict[str: np.ndarray]): + """Helper method to update `yi` values, accounting for possible `nans` by regression imputation. + + :param alpha: the model fidelity indices + :param beta: the surrogate fidelity indices + :param yi_dict: a `dict` mapping `str(coord)` grid coordinates to function values + """ + self.yi_map[str(alpha)].update(yi_dict) + imputer, xdim = None, len(self.x_vars) + for grid_coord, yi in yi_dict.items(): + if np.any(np.isnan(yi)): + if imputer is None: + # Grab all 'good' interpolation points and train a simple linear regression fit + xi_mat, yi_mat = np.zeros((0, xdim)), np.zeros((0, self.ydim)) + for coord, xi in self.xi_map[str(alpha)].items(): + if coord not in self.yi_nan_map[str(alpha)] and coord in self.yi_map[str(alpha)]: + yi_add = self.yi_map[str(alpha)][str(coord)] + xi_mat = np.concatenate((xi_mat, xi.reshape((1, xdim))), axis=0) + yi_mat = np.concatenate((yi_mat, yi_add.reshape((1, self.ydim))), axis=0) + nan_idx = np.any(np.isnan(yi_mat), axis=-1) + xi_mat = xi_mat[~nan_idx, :] + yi_mat = yi_mat[~nan_idx, :] + imputer = Pipeline([('scaler', MaxAbsScaler()), ('model', Ridge(alpha=1))]) + imputer.fit(xi_mat, yi_mat) + x_interp = self.xi_map[str(alpha)][str(grid_coord)].reshape((1, xdim)) + y_interp = np.atleast_1d(np.squeeze(imputer.predict(x_interp))) + nan_idx = np.isnan(yi) + y_interp[~nan_idx] = yi[~nan_idx] # Only keep imputed values where yi is nan + self.yi_nan_map[str(alpha)][str(grid_coord)] = y_interp + + # Go back and try to re-interpolate old nan values as more points get added to the grid + if imputer is not None: + for grid_coord in list(self.yi_nan_map[str(alpha)].keys()): + if grid_coord not in yi_dict: + x_interp = self.xi_map[str(alpha)][str(grid_coord)].reshape((1, xdim)) + y_interp = imputer.predict(x_interp) + self.yi_nan_map[str(alpha)][str(grid_coord)] = np.atleast_1d(np.squeeze(y_interp)) + + # Override + def get_sub_surrogate(self, alpha: tuple, beta: tuple, include_grid: bool = False) -> BaseInterpolator: + """Get the specific sub-surrogate corresponding to the $(\\alpha, \\beta)$ fidelity. + + :param alpha: A multi-index specifying model fidelity + :param beta: A multi-index specifying surrogate fidelity + :param include_grid: whether to add the `xi/yi` interpolation points to the returned `BaseInterpolator` object + :returns: the `BaseInterpolator` object corresponding to $(\\alpha, \\beta)$ + """ + interp = super().get_sub_surrogate(alpha, beta) + if include_grid: + interp.xi, interp.yi = self.get_tensor_grid(alpha, beta) + return interp + + def build_interpolator(self, alpha, beta): + """Abstract method implementation for constructing the tensor-product grid interpolator.""" + # Create a new tensor-product grid interpolator for the base index (0, 0, ...) + if np.sum(beta) == 0: + kwargs = copy.deepcopy(self._model_kwargs) + if len(alpha) > 0: + kwargs['alpha'] = alpha + interp = LagrangeInterpolator(beta, self.x_vars, model=self._model, model_args=self._model_args, + model_kwargs=kwargs, init_grids=True, reduced=True) + x_pt = np.array([float(interp.x_grids[n][beta[n]]) for n in range(interp.xdim())], dtype=np.float32) + self.curr_max_beta[str(alpha)] = list(beta) + self.x_grids[str(alpha)] = copy.deepcopy(interp.x_grids) + self.xi_map[str(alpha)] = {str(beta): x_pt} + self.yi_map[str(alpha)] = dict() + self.yi_nan_map[str(alpha)] = dict() + if self.save_enabled(): + self.yi_files[str(alpha)] = dict() + + return [beta], x_pt.reshape((1, len(self.x_vars))), interp + # Otherwise, all other indices are a refinement of previous grids + + # Look for first multi-index neighbor that is one level of refinement away + refine_tup = None + for beta_old_str in list(self.surrogates[str(alpha)].keys()): + beta_old = ast.literal_eval(beta_old_str) + if self.is_one_level_refinement(beta_old, beta): + idx_refine = int(np.nonzero(np.array(beta, dtype=int) - np.array(beta_old, dtype=int))[0]) + refine_level = beta[idx_refine] + if refine_level > self.curr_max_beta[str(alpha)][idx_refine]: + # Generate next refinement grid and save (refine_tup = tuple(x_new_idx, x_new, interp)) + refine_tup = self.surrogates[str(alpha)][beta_old_str].refine(beta, auto=False) + self.curr_max_beta[str(alpha)][idx_refine] = refine_level + self.x_grids[str(alpha)][idx_refine] = copy.deepcopy(refine_tup[2].x_grids[idx_refine]) + else: + # Access the refinement grid from memory (it is already computed) + num_pts = self.surrogates[str(alpha)][beta_old_str].get_grid_sizes(beta)[idx_refine] + x_refine = self.x_grids[str(alpha)][idx_refine][:num_pts] + refine_tup = self.surrogates[str(alpha)][beta_old_str].refine(beta, x_refine=x_refine, + auto=False) + break # Only need to grab one neighbor + + # Gather new interpolation grid points + x_new_idx, x_new, interp = refine_tup + xn_coord = [] # Save multi-index coordinates of points to compute model at for refinement + xn_pts = np.zeros((0, interp.xdim()), dtype=np.float32) # Save physical x location of new points + for i, multi_idx in enumerate(x_new_idx): + if str(multi_idx) not in self.yi_map[str(alpha)]: + # We have not computed this grid coordinate yet + xn_coord.append(multi_idx) + xn_pts = np.concatenate((xn_pts, x_new[i, np.newaxis, :]), axis=0) # (N_new, xdim) + self.xi_map[str(alpha)][str(multi_idx)] = x_new[i, :] + + return xn_coord, xn_pts, interp + + def update_interpolator(self, x_new_idx, x_new, interp): + """Awkward solution, I know, but actually compute and save the model evaluations here.""" + # Compute and store model output at new refinement points in a hash structure + yi_ret = interp.set_yi(x_new=(x_new_idx, x_new)) + + if self.ydim is None: + for coord_str, yi in yi_ret['y'].items(): + self.ydim = yi.shape[0] + break + + alpha = interp._model_kwargs.get('alpha', ()) + self.update_yi(alpha, interp.beta, yi_ret['y']) + if self.save_enabled(): + self.yi_files[str(alpha)].update(yi_ret['files']) + cost = interp.model_cost * len(x_new_idx) + + return cost + + def parallel_add_candidates(self, candidates: IndexSet, executor: Executor): + """Work-around to make sure mutable instance variable changes are made before/after + splitting tasks using this method over parallel (potentially MPI) workers. You can pass if you are not + interested in such parallel ideas. + + !!! Warning + MPI workers cannot save changes to `self` so this method should only distribute static tasks to the workers. + + :param candidates: list of [(alpha, beta),...] multi-indices + :param executor: the executor used to iterate candidates in parallel + """ + # Do sequential tasks first (i.e. make mutable changes to self), build up parallel task args + task_args = [] + for alpha, beta in candidates: + x_new_idx, x_new, interp = self.build_interpolator(alpha, beta) + task_args.append((alpha, beta, x_new_idx, x_new, interp)) + + def parallel_task(alpha, beta, x_new_idx, x_new, interp): + # Must return anything you want changed in self or interp (mutable changes aren't saved over MPI workers) + logger = get_logger(self.__class__.__name__, log_file=self.log_file, stdout=False) + logger.info(f'Building interpolator for index {(alpha, beta)} ...') + yi_ret = interp.set_yi(x_new=(x_new_idx, x_new)) + model_cost = interp.model_cost if interp.model_cost is not None else 1 + return yi_ret, model_cost + + # Wait for all parallel workers to return + fs = [executor.submit(parallel_task, *args) for args in task_args] + wait(fs, timeout=None, return_when=ALL_COMPLETED) + + # Update self and interp with the results from all workers (and check for errors) + for i, future in enumerate(fs): + try: + a = task_args[i][0] + b = task_args[i][1] + x_new_idx = task_args[i][2] + interp = task_args[i][4] + yi_ret, model_cost = future.result() + interp.model_cost = model_cost + self.surrogates[str(a)][str(b)] = interp + self.update_yi(a, b, yi_ret['y']) + if self.save_enabled(): + self.yi_files[str(a)].update(yi_ret['files']) + self.costs[str(a)][str(b)] = interp.model_cost * len(x_new_idx) + + if self.ydim is None: + for coord_str, yi in self.yi_map[str(a)].items(): + self.ydim = yi.shape[0] + break + except: + self.logger.error(f'An exception occurred in a thread handling build_interpolator{candidates[i]}') + raise + + +class AnalyticalSurrogate(ComponentSurrogate): + """Concrete "surrogate" class that just uses the analytical model (i.e. bypasses surrogate evaluation).""" + + def __init__(self, x_vars, model, *args, **kwargs): + """Initializes a stand-in `ComponentSurrogate` with all unnecessary fields set to empty. + + !!! Warning + This overwrites anything passed in for `truth_alpha`, `max_alpha`, `max_beta`, or `multi_index` since + these are not used for an analytical model. + """ + kwargs['truth_alpha'] = () + kwargs['max_alpha'] = () + kwargs['max_beta'] = () + kwargs['multi_index'] = [] + super().__init__(x_vars, model, *args, **kwargs) + + # Override + def predict(self, x: np.ndarray | float, **kwargs) -> np.ndarray: + """Evaluate the analytical model at points `x`, ignore extra `**kwargs` passed in. + + :param x: `(..., x_dim)` the points to be evaluated + :returns y: `(..., y_dim)` the exact model output at the input points + """ + ret = self._model(x, *self._model_args, **self._model_kwargs) + + if not isinstance(ret, dict): + self.logger.warning(f"Function {self._model} did not return a dict of the form {{'y': y}}. Please make sure" + f" you do so to avoid conflicts. Returning the value directly instead...") + + return ret['y'] if isinstance(ret, dict) else ret + + # Override + def activate_index(self, *args): + """Do nothing""" + pass + + # Override + def add_surrogate(self, *args): + """Do nothing""" + pass + + # Override + def init_coarse(self): + """Do nothing""" + pass + + # Override + def update_misc_coeffs(self, **kwargs): + """Do nothing""" + pass + + # Override + def get_sub_surrogate(self, *args): + """Nothing to return for analytical model""" + return None + + # Override + def get_cost(self, *args): + """Return no cost for analytical model""" + return 0 + + def build_interpolator(self, *args): + """Abstract method implementation, return none for an analytical model""" + return None + + def update_interpolator(self, *args): + """Abstract method implementation, return `cost=0` for an analytical model""" + return 0 + + def parallel_add_candidates(self, *args): + """Abstract method implementation, do nothing""" + pass diff --git a/src/amisc/examples/models.py b/src/amisc/examples/models.py new file mode 100644 index 0000000..3477d1a --- /dev/null +++ b/src/amisc/examples/models.py @@ -0,0 +1,237 @@ +"""Providing some example test functions""" +from pathlib import Path +import pickle +import uuid + +import numpy as np + +from amisc.system import ComponentSpec, SystemSurrogate +from amisc.rv import UniformRV, NormalRV + + +def tanh_func(x, *args, A=2, L=1, frac=4, **kwargs): + """Simple tunable tanh function""" + return {'y': A*np.tanh(2/(L/frac)*(x-L/2)) + A + 1} + + +def ishigami(x, *args, a=7.0, b=0.1, **kwargs): + """For testing Sobol indices: https://doi.org/10.1109/ISUMA.1990.151285""" + return {'y': np.sin(x[..., 0:1]) + a*np.sin(x[..., 1:2])**2 + b*(x[..., 2:3]**4)*np.sin(x[..., 0:1])} + + +def borehole_func(x, *args, **kwargs): + """Model found at https://www.sfu.ca/~ssurjano/borehole.html + :returns vdot: Water flow rate in m^3/yr + """ + rw = x[..., 0] # Radius of borehole (m) + r = x[..., 1] # Radius of influence (m) + Tu = x[..., 2] # Transmissivity (m^2/yr) + Hu = x[..., 3] # Potentiometric head (m) + Tl = x[..., 4] # Transmissivity (m^2/yr) + Hl = x[..., 5] # Potentiometric head (m) + L = x[..., 6] # Length of borehole (m) + Kw = x[..., 7] # Hydraulic conductivity (m/yr) + vdot = 2*np.pi*Tu*(Hu-Hl) / (np.log(r/rw) * (1 + (2*L*Tu/(np.log(r/rw)*Kw*rw**2)) + (Tu/Tl))) + + return {'y': vdot[..., np.newaxis]} + + +def wing_weight_func(x, *args, **kwargs): + """Model found at https://www.sfu.ca/~ssurjano/wingweight.html + :returns Wwing: the weight of the airplane wing (lb) + """ + Sw = x[..., 0] # Wing area (ft^2) + Wfw = x[..., 1] # Weight of fuel (lb) + A = x[..., 2] # Aspect ratio + Lambda = x[..., 3] # Quarter-chord sweep (deg) + q = x[..., 4] # Dynamic pressure (lb/ft^2) + lamb = x[..., 5] # taper ratio + tc = x[..., 6] # Aerofoil thickness to chord ratio + Nz = x[..., 7] # Ultimate load factor + Wdg = x[..., 8] # Design gross weight (lb) + Wp = x[..., 9] # Paint weight (lb/ft^2) + Lambda = Lambda*(np.pi/180) + Wwing = 0.036*(Sw**0.758)*(Wfw**0.0035)*((A/(np.cos(Lambda))**2)**0.6)*(q**0.006)*(lamb**0.04)*\ + (100*tc/np.cos(Lambda))**(-0.3)*((Nz*Wdg)**0.49) + Sw*Wp + + return {'y': Wwing[..., np.newaxis]} + + +def nonlinear_wave(x, *args, env_var=0.1**2, wavelength=0.5, wave_amp=0.1, tanh_amp=0.5, L=1, t=0.25, **kwargs): + """Custom nonlinear model of a traveling Gaussian wave for testing. + + :param x: `(..., x_dim)`, input locations + :param env_var: variance of Gaussian envelope + :param wavelength: sinusoidal perturbation wavelength + :param wave_amp: amplitude of perturbation + :param tanh_amp: amplitude of tanh(x) + :param L: domain length of underlying tanh function + :param t: transition length of tanh function (as fraction of L) + :returns: `(..., y_dim)`, model output + """ + # Traveling sinusoid with moving Gaussian envelope (theta is x2) + env_range = [0.2, 0.6] + mu = env_range[0] + x[..., 1] * (env_range[1] - env_range[0]) + theta_env = 1 / (np.sqrt(2 * np.pi * env_var)) * np.exp(-0.5 * (x[..., 0] - mu) ** 2 / env_var) + ftheta = wave_amp * np.sin((2*np.pi/wavelength) * x[..., 1]) * theta_env + + # Underlying tanh dependence on x1 + fd = tanh_amp * np.tanh(2/(L*t)*(x[..., 0] - L/2)) + tanh_amp + + # Compute model = f(theta, d) + f(d) + y = np.expand_dims(ftheta + fd, axis=-1) # (..., 1) + + return {'y': y} + + +def fire_sat_system(save_dir=None): + """Fire satellite detection system model from Chaudhuri 2018. + + !!! Note "Some modifications" + Orbit will save outputs all the time; Power won't because it is part of an FPI loop. Orbit and Power can + return `np.nan` some of the time (optional, to test robustness of surrogate). Power and Attitude have an `alpha` fidelity + index that controls accuracy of the returned values. `alpha=0,1,2` corresponds to `8,4,0` percent error. + + :param save_dir: where to save model outputs + :returns: a `SystemSurrogate` object for the fire sat MD example system + """ + Re = 6378140 # Radius of Earth (m) + mu = 3.986e14 # Gravitational parameter (m^3 s^-2) + eta = 0.22 # Power efficiency + Id = 0.77 # Inherent degradation of the array + thetai = 0 # Sun incidence angle + LT = 15 # Spacecraft lifetime (years) + eps = 0.0375 # Power production degradation (%/year) + rlw = 3 # Length to width ratio + nsa = 3 # Number of solar arrays + rho = 700 # Mass density of arrays (kg/m^3) + t = 0.005 # Thickness (m) + D = 2 # Distance between panels (m) + IbodyX = 6200 # kg*m^2 + IbodyY = 6200 # kg*m^2 + IbodyZ = 4700 # kg*m^2 + dt_slew = 760 # s + theta = 15 # Deviation of moment axis from vertical (deg) + As = 13.85 # Area reflecting radiation (m^2) + c = 2.9979e8 # Speed of light (m/s) + M = 7.96e15 # Magnetic moment of earth (A*m^2) + Rd = 5 # Residual dipole of spacecraft (A*m^2) + rhoa=5.148e-11 # Atmospheric density (kg/m^3) -- typo in Chaudhuri 2018 has this as 1e11 instead + A = 13.85 # Cross-section in flight (m^2) + Phold = 20 # Holding power (W) + omega = 6000 # Max vel of wheel (rpm) + nrw = 3 # Number of reaction wheels + + def orbit_fun(x, output_dir=None, pct_failure=0): + H = x[..., 0:1] # Altitude (m) + phi = x[..., 1:2] # Target diameter (m) + vel = np.sqrt(mu / (Re + H)) # Satellite velocity (m/s) + dt_orbit = 2*np.pi*(Re + H) / vel # Orbit period (s) + dt_eclipse = (dt_orbit/np.pi)*np.arcsin(Re / (Re + H)) # Eclipse period (s) + theta_slew = np.arctan(np.sin(phi / Re) / (1 - np.cos(phi / Re) + H/Re)) # Max slew angle (rad) + if np.random.rand() < pct_failure: + i = tuple([np.random.randint(0, N) for N in x.shape[:-1]]) + i2 = tuple([np.random.randint(0, N) for N in x.shape[:-1]]) + vel[i + (0,)] = np.nan + theta_slew[i2 + (0,)] = np.nan + y = np.concatenate((vel, dt_orbit, dt_eclipse, theta_slew), axis=-1) + if output_dir is not None: + files = [] + id = str(uuid.uuid4()) + for index in np.ndindex(*x.shape[:-1]): + fname = f'{id}_{index}.pkl' + with open(Path(output_dir) / fname, 'wb') as fd: + pickle.dump({'y': y[index + (slice(None),)]}, fd) + files.append(fname) + return {'y': y, 'files': files} + else: + return {'y': y} + + def power_fun(x, alpha=(0,), *, output_dir=None, pct_failure=0): + pct = 1 - (2 - alpha[0]) * 0.04 if len(alpha) == 1 else 1 # extra pct error term + Po = x[..., 0:1] # Other power sources (W) + Fs = x[..., 1:2] # Solar flux (W/m^2) + dt_orbit = x[..., 2:3] # Orbit period (s) + dt_eclipse = x[..., 3:4] # Eclipse period (s) + Pacs = x[..., 4:5] # Power from attitude control system (W) + Ptot = Po + Pacs + Pe = Ptot + Pd = Ptot + Xe = 0.6 # These are power efficiencies in eclipse and daylight + Xd = 0.8 # See Ch. 11 of Wertz 1999 SMAD + Te = dt_eclipse + Td = dt_orbit - Te + Psa = ((Pe*Te/Xe) + (Pd*Td/Xd)) / Td + Pbol = eta*Fs*Id*np.cos(thetai) + Peol = Pbol * (1 - eps)**LT + Asa = Psa / Peol # Total solar array size (m^2) + L = np.sqrt(Asa*rlw/nsa) + W = np.sqrt(Asa/(rlw*nsa)) + msa = 2*rho*L*W*t # Mass of solar array + Ix = msa*((1/12)*(L**2 + t**2) + (D+L/2)**2) + IbodyX + Iy = (msa/12)*(L**2 + W**2) + IbodyY # typo in Zaman 2013 has this as H**2 instead of L**2 + Iz = msa*((1/12)*(L**2 + W**2) + (D + L/2)**2) + IbodyZ + Itot = np.concatenate((Ix, Iy, Iz), axis=-1) + Imin = np.min(Itot, axis=-1, keepdims=True) + Imax = np.max(Itot, axis=-1, keepdims=True) + if np.random.rand() < pct_failure: + i = tuple([np.random.randint(0, N) for N in x.shape[:-1]]) + i2 = tuple([np.random.randint(0, N) for N in x.shape[:-1]]) + Imin[i2 + (0,)] = np.nan + Asa[i + (0,)] = np.nan + y = np.concatenate((Imin, Imax*pct, Ptot*pct, Asa), axis=-1) + + if output_dir is not None: + files = [] + id = str(uuid.uuid4()) + for index in np.ndindex(*x.shape[:-1]): + fname = f'{id}_{index}.pkl' + with open(Path(output_dir) / fname, 'wb') as fd: + pickle.dump({'y': y[index + (slice(None),)]}, fd) + files.append(fname) + return {'y': y, 'files': files} + else: + return {'y': y} + + def attitude_fun(x, /, alpha=(0,)): + pct = 1 - (2 - alpha[0])*0.04 if len(alpha) == 1 else 1 # extra model fidelity pct error term + H = x[..., 0:1] # Altitude (m) + Fs = x[..., 1:2] # Solar flux + Lsp = x[..., 2:3] # Moment arm for solar radiation pressure + q = x[..., 3:4] # Reflectance factor + La = x[..., 4:5] # Moment arm for aerodynamic drag + Cd = x[..., 5:6] # Drag coefficient + vel = x[..., 6:7] # Satellite velocity + theta_slew = x[..., 7:8] # Max slew angle + Imin = x[..., 8:9] # Minimum moment of inertia + Imax = x[..., 9:10] # Maximum moment of inertia + tau_slew = 4*theta_slew*Imax / dt_slew**2 + tau_g = 3*mu*np.abs(Imax - Imin)*np.sin(2*theta*(np.pi/180)) / (2*(Re+H)**3) + tau_sp = Lsp*Fs*As*(1+q)*np.cos(thetai) / c + tau_m = 2*M*Rd / (Re + H)**3 + tau_a = (1/2)*La*rhoa*Cd*A*vel**2 + tau_dist = np.sqrt(tau_g**2 + tau_sp**2 + tau_m**2 + tau_a**2) + tau_tot = np.max(np.concatenate((tau_slew, tau_dist), axis=-1), axis=-1, keepdims=True) + Pacs = tau_tot*(omega*(2*np.pi/60)) + nrw*Phold + y = np.concatenate((Pacs, tau_tot), axis=-1) * pct + return {'y': y} + + orbit = ComponentSpec(orbit_fun, name='Orbit', exo_in=[0, 1], coupling_out=[0, 1, 2, 3], max_beta=(3, 3), + model_kwargs={'pct_failure': 0}, save_output=True) + power = ComponentSpec(power_fun, name='Power', truth_alpha=(2,), exo_in=[2, 3], max_alpha=(2,), max_beta=(3,)*5, + coupling_in={'Orbit': [1, 2], 'Attitude': [0]}, coupling_out=[4, 5, 6, 7], save_output=True, + model_kwargs={'pct_failure': 0}) + attitude = ComponentSpec(attitude_fun, name='Attitude', truth_alpha=2, max_alpha=2, max_beta=(3,)*10, + exo_in=[0, 3, 4, 5, 6, 7], coupling_in={'Orbit': [0, 3], 'Power': [0, 1]}, + coupling_out=[8, 9]) + + exo_vars = [NormalRV(18e6, 1e6, 'H'), NormalRV(235e3, 10e3, '\u03D5'), NormalRV(1000, 50, 'Po'), + NormalRV(1400, 20, 'Fs'), NormalRV(2, 0.4, 'Lsp'), NormalRV(0.5, 0.1, 'q'), + NormalRV(2, 0.4, 'La'), NormalRV(1, 0.2, 'Cd')] + coupling_vars = [UniformRV(2000, 6000, 'Vsat'), UniformRV(20000, 60000, 'To'), UniformRV(1000, 5000, 'Te'), + UniformRV(0, 4, 'Slew'), UniformRV(0, 12000, 'Imin'), UniformRV(0, 12000, 'Imax'), + UniformRV(0, 10000, 'Ptot'), UniformRV(0, 50, 'Asa'), UniformRV(0, 100, 'Pat'), + UniformRV(0, 5, 'tau_tot')] + surr = SystemSurrogate([orbit, power, attitude], exo_vars, coupling_vars, est_bds=500, save_dir=save_dir) + + return surr diff --git a/src/amisc/examples/tutorial.py b/src/amisc/examples/tutorial.py new file mode 100644 index 0000000..5c3515c --- /dev/null +++ b/src/amisc/examples/tutorial.py @@ -0,0 +1,82 @@ +"""Examples to get started.""" + + +def single_component(): + # --8<-- [start:single] + from amisc.system import SystemSurrogate, ComponentSpec + from amisc.utils import UniformRV + + def fun(x): + return dict(y=x ** 2) + + x = UniformRV(-1, 1) + y = UniformRV(0, 1) + component = ComponentSpec(fun) + system = SystemSurrogate([component], x, y) + + system.fit() + system.predict(0.5) # 0.25 + # --8<-- [end:single] + + +def simple(): + # --8<-- [start:simple] + import numpy as np + + from amisc.system import SystemSurrogate, ComponentSpec + from amisc.utils import UniformRV + + def fun1(x): + return dict(y=x * np.sin(np.pi * x)) + + def fun2(x): + return dict(y=1 / (1 + 25 * x ** 2)) + + x = UniformRV(0, 1, 'x') + y = UniformRV(0, 1, 'y') + z = UniformRV(0, 1, 'z') + model1 = ComponentSpec(fun1, exo_in=x, coupling_out=y) + model2 = ComponentSpec(fun2, coupling_in=y, coupling_out=z) + + inputs = x + outputs = [y, z] + system = SystemSurrogate([model1, model2], inputs, outputs) + system.fit() + + x_test = system.sample_inputs(10) + y_test = system.predict(x_test) + # --8<-- [end:simple] + + +def fire_sat(): + # --8<-- [start:fire_sat] + import numpy as np + import matplotlib.pyplot as plt + + from amisc.examples.models import fire_sat_system + + system = fire_sat_system() + + xtest = system.sample_inputs(100, use_pdf=True) # --> (100, xdim) + ytest = system(xtest, use_model='best') # --> (100, ydim) + use_idx = ~np.any(np.isnan(ytest), axis=-1) + xtest = xtest[use_idx, :] + ytest = ytest[use_idx, :] + test_set = {'xt': xtest, 'yt': ytest} + + system.fit(max_iter=10, test_set=test_set, n_jobs=-1, num_refine=1000) + + print(f'Inputs: {system.exo_vars}') + print(f'Outputs: {system.coupling_vars}') + + # Plots + input_vars = ['H', 'Po'] + output_vars = ['Vsat', 'Asa'] + system.plot_allocation() + system.plot_slice(input_vars, output_vars, show_model=['best', 'worst'], random_walk=True, N=10) + plt.show() + # --8<-- [end:fire_sat] + + +if __name__ == '__main__': + single_component() diff --git a/src/amisc/interpolator.py b/src/amisc/interpolator.py new file mode 100644 index 0000000..6f97a6a --- /dev/null +++ b/src/amisc/interpolator.py @@ -0,0 +1,465 @@ +"""`interpolator.py` + +Provides interpolator classes. Interpolators manage training data and specify how to refine/gather more data. + +Includes +-------- +- `BaseInterpolator`: Abstract class providing basic structure of an interpolator +- `LagrangeInterpolator`: Concrete implementation for tensor-product barycentric Lagrange interpolation +""" +from abc import ABC, abstractmethod +import itertools +import copy + +import numpy as np +from scipy.optimize import direct +from sklearn.linear_model import Ridge +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import MaxAbsScaler + +from amisc.utils import get_logger +from amisc.rv import BaseRV + + +class BaseInterpolator(ABC): + """Base interpolator abstract class. + + !!! Info "Setting the training data" + You can leave the training data `xi`, `yi` empty; they can be iteratively refined later on. + + !!! Info "Model specification" + The model is a callable function of the form `ret = model(x, *args, **kwargs)`. The return value is a dictionary + of the form `ret = {'y': y, 'files': files, 'cost': cost}`. In the return dictionary, you specify the raw model + output `y` as an `np.ndarray` at a _minimum_. Optionally, you can specify paths to output files and the average + model cost (in units of seconds of cpu time), and anything else you want. + + :ivar beta: specifies the refinement level of this interpolator as a set of natural number indices + :ivar x_vars: list of variables that fully determines the input domain of interest for interpolation + :ivar xi: `(Nx, x_dim)`, interpolation points (or knots, training samples, etc.) stored as an array + :ivar yi: `(Nx, y_dim)`, function values at the interpolation points, i.e. the training data + :ivar _model: stores a ref to the model or function that is to be interpolated, callable as `ret = model(x)` + :ivar _model_args: additional arguments to supply to the model + :ivar _model_kwargs: additional keyword arguments to supply to the model + :ivar model_cost: the average total cpu time (in seconds) for a single model evaluation call of one set of inputs + :ivar output_files: tracks model output files corresponding to `yi` training data (for more complex models) + :ivar logger: a logging utility reference + + :vartype beta: tuple[int, ...] + :vartype x_vars: list[BaseRV] + :vartype xi: np.ndarray + :vartype yi: np.ndarray + :vartype _model: callable[np.ndarray] -> dict + :vartype _model_args: tuple + :vartype _model_kwargs: dict + :vartype model_cost: float + :vartype output_files: list[str | Path] + :vartype logger: logging.Logger + """ + + def __init__(self, beta: tuple, x_vars: BaseRV | list[BaseRV], xi=None, yi=None, + model=None, model_args=(), model_kwargs=None): + """Construct the interpolator. + + :param beta: refinement level indices + :param x_vars: list of variables to specify input domain of interpolation + :param xi: `(Nx, xdim)`, interpolation points (optional) + :param yi: `(Nx, ydim)`, the function values at the interpolation points (optional) + :param model: callable as {'y': y} = model(x), with `x = (..., x_dim)`, `y = (..., y_dim)` + :param model_args: optional args for the model + :param model_kwargs: optional kwargs for the model + """ + x_vars = [x_vars] if not isinstance(x_vars, list) else x_vars + self.logger = get_logger(self.__class__.__name__) + self._model = model + self._model_args = model_args + self._model_kwargs = model_kwargs if model_kwargs is not None else {} + self.output_files = [] # Save output files with same indexing as xi, yi + self.xi = xi # Interpolation points + self.yi = yi # Function values at interpolation points + self.beta = beta # Refinement level indices + self.x_vars = x_vars # BaseRV() objects for each input + self.model_cost = None # Total cpu time to evaluate model once (s) + + def update_input_bds(self, idx: int, bds: tuple): + """Update the input bounds at the given index. + + :param idx: the index of the input variable to update + :param bds: the new bounds for the variable + """ + self.x_vars[idx].update_bounds(*bds) + + def xdim(self): + """Get the dimension of the input domain.""" + return len(self.x_vars) + + def ydim(self): + """Get the dimension of the outputs.""" + return self.yi.shape[-1] if self.yi is not None else None + + def save_enabled(self): + """Return whether the underlying model wants to save outputs to file. + + !!! Note + You can specify that a model wants to save outputs to file by providing an `'output_dir'` kwarg. + """ + return self._model_kwargs.get('output_dir') is not None + + def set_yi(self, yi: np.ndarray = None, model: callable = None, + x_new: tuple[list[int | tuple], np.ndarray] = ()) -> dict[str: np.ndarray] | None: + """Set the training data; if `yi=None`, then compute from the model. + + !!! Warning + You would use `x_new` if you wanted to compute the model at these specific locations and store the result. + This will ignore anything passed in for `yi`, and it assumes a model is already specified (or passed in). + + !!! Info + You can pass in integer indices for `x_new` or tuple indices. Integers will index into `self.xi`. Tuples + provide extra flexibility for more complicated indexing, e.g. they might specify indices along different + coordinate directions in an N-dimensional grid. If you pass in a list of tuple indices for `x_new`, the + resulting model outputs will be returned back to you in the form `dict[str: np.ndarray]`. The keys are + string casts of the tuple indices, and the values are the corresponding model outputs. + + :param yi: `(Nx, y_dim)`, training data to set, must match dimension of `self.xi` + :param model: callable function, optionally overrides `self._model` + :param x_new: tuple of `(idx, x)`, where `x` is an `(N_new, x_dim)` array of new interpolation points to + include and `idx` specifies the indices of these points into `self.xi` + :returns: dict[str: np.ndarray] if `idx` contains tuple elements, otherwise `None` + """ + if model is not None: + self._model = model + if self._model is None: + error_msg = 'Model not specified for computing QoIs at interpolation grid points.' + self.logger.error(error_msg) + raise Exception(error_msg) + + # Overrides anything passed in for yi (you would only be using this if yi was set previously) + if x_new: + new_idx = x_new[0] + new_x = x_new[1] + return_y = isinstance(new_idx[0], tuple) # Return y rather than storing it if tuple indices are passed in + ret = dict(y=dict(), files=dict()) + model_ret = self._model(new_x, *self._model_args, **self._model_kwargs) + if not isinstance(model_ret, dict): + self.logger.warning( + f"Function {self._model} did not return a dict of the form {{'y': y}}. Please make sure" + f" you do so to avoid conflicts. Returning the value directly instead...") + model_ret = dict(y=model_ret) + y_new, files_new, cpu_time = model_ret['y'], model_ret.get('files', None), model_ret.get('cost', 1) + + if self.save_enabled(): + for j in range(y_new.shape[0]): + if return_y: + ret['y'][str(new_idx[j])] = y_new[j, :].astype(np.float32) + ret['files'][str(new_idx[j])] = files_new[j] + else: + self.yi[new_idx[j], :] = y_new[j, :].astype(np.float32) + self.output_files[new_idx[j]] = files_new[j] + else: + for j in range(y_new.shape[0]): + if return_y: + ret['y'][str(new_idx[j])] = y_new[j, :].astype(np.float32) + else: + self.yi[new_idx[j], :] = y_new[j, :].astype(np.float32) + + if self.model_cost is None: + self.model_cost = max(1, cpu_time) + + return ret + + # Set yi directly + if yi is not None: + self.yi = yi.astype(np.float32) + return + + # Compute yi + model_ret = self._model(self.xi, *self._model_args, **self._model_kwargs) + if not isinstance(model_ret, dict): + self.logger.warning(f"Function {self._model} did not return a dict of the form {{'y': y}}. Please make sure" + f" you do so to avoid conflicts. Returning the value directly instead...") + model_ret = dict(y=model_ret) + + self.yi, self.output_files, cpu_time = model_ret['y'], model_ret.get('files', list()), model_ret.get('cost', 1) + + if self.model_cost is None: + self.model_cost = max(1, cpu_time) + + @abstractmethod + def refine(self, beta: tuple, auto=True): + """Return a new interpolator with one dimension refined by one level, as specified by `beta`. + + !!! Info "When you want to compute the model manually" + You can set `auto=False`, in which case the newly refined interpolation points `x` will be returned to you + along with their indices, in the form `idx, x, interp = refine(beta, auto=False)`. You might also want to + do this if you did not provide a model when constructing the Interpolator (so `auto=True` won't work). + + :param beta: the new refinement level indices, should only refine one dimension by one level + :param auto: whether to automatically compute and store model at refinement points (default is True) + :returns: `idx` - indices into `xi`, `x` - the new interpolation points, and `interp` - a refined + BaseInterpolator object, just returns `interp` if `auto=True` + """ + pass + + @abstractmethod + def __call__(self, x: np.ndarray | float) -> np.ndarray: + """Evaluate the interpolation at points `x`. + + :param x: `(..., x_dim)`, the points to be interpolated, must be within the input domain for accuracy + :returns y: `(..., y_dim)`, the interpolated function values + """ + pass + + +class LagrangeInterpolator(BaseInterpolator): + """Tensor-product (multivariate) grid interpolator, based on barycentric Lagrange polynomials. + + !!! Info + The refinement level indices `beta` are used in this class to specify anisotropic refinement along each + coordinate direction of the input domain, so `x_dim = len(beta)`. + + :ivar x_grids: univariate Leja sequence points in each 1d dimension + :ivar weights: the barycentric weights corresponding to `x_grids` + :ivar reduced: whether to store `xi` and `yi` training data, can set to `False` to save memory, e.g. if an external + sparse grid data structure manages this data instead + + :vartype x_grids: list[np.ndarray] + :vartype weights: list[np.ndarray] + :vartype reduced: bool + """ + + def __init__(self, beta: tuple, x_vars: BaseRV | list[BaseRV], init_grids=True, reduced=False, **kwargs): + """Initialize a Lagrange tensor-product grid interpolator. + + :param beta: refinement level indices for each input dimension + :param x_vars: list of variables specifying bounds/pdfs for each input x + :param init_grids: whether to compute 1d Leja sequences on initialization + :param reduced: whether to store xi/yi matrices, e.g. set true if storing in external sparse grid structure + :param **kwargs: other optional arguments (see `BaseInterpolator`) + """ + self.weights = [] # Barycentric weights for each dimension + self.x_grids = [] # Univariate nested leja sequences in each dimension + self.reduced = reduced + super().__init__(beta, x_vars, **kwargs) + + if init_grids: + # Construct 1d univariate Leja sequences in each dimension + grid_sizes = self.get_grid_sizes(self.beta) + self.x_grids = [self.leja_1d(grid_sizes[n], self.x_vars[n].bounds(), + wt_fcn=self.x_vars[n].pdf).astype(np.float32) for n in range(self.xdim())] + + for n in range(self.xdim()): + Nx = grid_sizes[n] + bds = self.x_vars[n].bounds() + grid = self.x_grids[n] + C = (bds[1] - bds[0]) / 4.0 # Interval capacity (see Berrut and Trefethen 2004) + xj = grid.reshape((Nx, 1)) + xi = grid.reshape((1, Nx)) + dist = (xj - xi) / C + np.fill_diagonal(dist, 1) # Ignore product when i==j + self.weights.append((1.0 / np.prod(dist, axis=1)).astype(np.float32)) # (Nx,) + + # Cartesian product of univariate grids + if not self.reduced: + self.xi = np.empty((np.prod(grid_sizes), self.xdim()), dtype=np.float32) + for i, ele in enumerate(itertools.product(*self.x_grids)): + self.xi[i, :] = ele + + def refine(self, beta, auto=True, x_refine: np.ndarray = None): + """Return a new interpolator with one dimension refined by one level, specified by `beta`. + + !!! Note + If `self.reduced=True` or `auto=False`, then this function will return tuple indices `idx` corresponding + to the new interpolation points `x`. The tuple indices specify one index along each input dimension. + + :param beta: the new refinement level indices + :param auto: whether to automatically compute model at refinement points + :param x_refine: `(Nx,)` use this array as the refined 1d grid if provided, otherwise compute via `leja_1d` + :returns: `interp` - a `LagrangeInterpolator` with a refined grid (default), otherwise if `auto=False`, + returns `idx, x, interp`, where `idx` and `x` correspond to new interpolation points. + """ + try: + # Initialize a new interpolator with the new refinement levels + interp = LagrangeInterpolator(beta, self.x_vars, model=self._model, model_args=self._model_args, + model_kwargs=self._model_kwargs, init_grids=False, reduced=self.reduced) + + # Find the dimension and number of new points to add + old_grid_sizes = self.get_grid_sizes(self.beta) + new_grid_sizes = interp.get_grid_sizes(beta) + dim_refine = 0 + num_new_pts = 0 + for idx, grid_size in enumerate(new_grid_sizes): + if grid_size != old_grid_sizes[idx]: + dim_refine = idx + num_new_pts = grid_size - old_grid_sizes[idx] + break + + # Add points to leja grid in this dimension + interp.x_grids = copy.deepcopy(self.x_grids) + xi = copy.deepcopy(x_refine) if x_refine is not None else self.leja_1d(num_new_pts, + interp.x_vars[dim_refine].bounds(), + z_pts=interp.x_grids[dim_refine], + wt_fcn=interp.x_vars[dim_refine].pdf) + interp.x_grids[dim_refine] = xi.astype(np.float32) + + # Update barycentric weights in this dimension + interp.weights = copy.deepcopy(self.weights) + Nx_old = old_grid_sizes[dim_refine] + Nx_new = new_grid_sizes[dim_refine] + old_wts = copy.deepcopy(self.weights[dim_refine]) + new_wts = np.zeros(Nx_new, dtype=np.float32) + new_wts[:Nx_old] = old_wts + bds = interp.x_vars[dim_refine].bounds() + C = (bds[1] - bds[0]) / 4.0 # Interval capacity + xi = interp.x_grids[dim_refine] + for j in range(Nx_old, Nx_new): + new_wts[:j] *= (C / (xi[:j] - xi[j])) + new_wts[j] = np.prod(C / (xi[j] - xi[:j])) + interp.weights[dim_refine] = new_wts + + # Copy yi over at existing interpolation points + x_new = np.zeros((0, interp.xdim()), dtype=np.float32) + x_new_idx = [] + tol = 1e-12 # Tolerance for floating point comparison + j = 0 # Use this idx for iterating over existing yi + if not self.reduced: + interp.xi = np.zeros((np.prod(new_grid_sizes), self.xdim()), dtype=np.float32) + interp.yi = np.zeros((np.prod(new_grid_sizes), self.ydim()), dtype=np.float32) + if self.save_enabled(): + interp.output_files = [None] * np.prod(new_grid_sizes) + + old_indices = [np.arange(old_grid_sizes[n]) for n in range(self.xdim())] + old_indices = list(itertools.product(*old_indices)) + new_indices = [np.arange(new_grid_sizes[n]) for n in range(self.xdim())] + new_indices = list(itertools.product(*new_indices)) + for i in range(len(new_indices)): + # Get the new grid coordinate/index and physical x location/point + new_x_idx = new_indices[i] + new_x_pt = np.array([float(interp.x_grids[n][new_x_idx[n]]) for n in range(self.xdim())], + dtype=np.float32) + + if not self.reduced: + # Store the old xi/yi and return new x points + interp.xi[i, :] = new_x_pt + if j < len(old_indices) and np.all(np.abs(np.array(old_indices[j]) - + np.array(new_indices[i])) < tol): + # If we already have this interpolation point + interp.yi[i, :] = self.yi[j, :] + if self.save_enabled(): + interp.output_files[i] = self.output_files[j] + j += 1 + else: + # Otherwise, save new interpolation point and its index + x_new = np.concatenate((x_new, new_x_pt.reshape((1, self.xdim())))) + x_new_idx.append(i) + else: + # Just find the new x indices and return those for the reduced case + if j < len(old_indices) and np.all(np.abs(np.array(old_indices[j]) - + np.array(new_indices[i])) < tol): + j += 1 + else: + x_new = np.concatenate((x_new, new_x_pt.reshape((1, self.xdim())))) + x_new_idx.append(new_x_idx) # Add a tuple() multi-index if not saving xi/yi + + # Evaluate the model at new interpolation points + interp.model_cost = self.model_cost + if self._model is None: + self.logger.warning(f'No model available to evaluate new interpolation points, returning the points ' + f'to you instead...') + return x_new_idx, x_new, interp + elif not auto or self.reduced: + return x_new_idx, x_new, interp + else: + interp.set_yi(x_new=(x_new_idx, x_new)) + return interp + + except Exception as e: + import traceback + tb_str = str(traceback.format_exception(e)) + self.logger.error(tb_str) + raise Exception(f'Original exception in refine(): {tb_str}') + + def __call__(self, x: np.ndarray | float, xi: np.ndarray = None, yi: np.ndarray = None) -> np.ndarray: + """Evaluate the barycentric interpolation at points `x`. + + :param x: `(..., xdim)`, the points to be interpolated, must be within domain of `self.xi` for accuracy + :param xi: `(Ni, xdim)` optional, interpolation grid points to use (e.g. if `self.reduced=True`) + :param yi: `(Ni, ydim)` optional, function values at xi to use (e.g. if `self.reduced=True`) + :returns y: `(..., ydim)`, the interpolated function values + """ + x = np.atleast_1d(x) + if yi is None: + yi = self.yi.copy() + if xi is None: + xi = self.xi.copy() + ydim = yi.shape[-1] + nan_idx = np.any(np.isnan(yi), axis=-1) + if np.any(nan_idx): + # Use a simple linear regression fit to impute missing values (may have resulted from bad model outputs) + imputer = Pipeline([('scaler', MaxAbsScaler()), ('model', Ridge(alpha=1))]) + imputer.fit(xi[~nan_idx, :], yi[~nan_idx, :]) + yi[nan_idx, :] = imputer.predict(xi[nan_idx, :]) + + # Loop over multi-indices and compute tensor-product lagrange polynomials + grid_sizes = self.get_grid_sizes(self.beta) + y = np.zeros(x.shape[:-1] + (ydim,)) # (..., ydim) + indices = [np.arange(grid_sizes[n]) for n in range(self.xdim())] + for i, j in enumerate(itertools.product(*indices)): + L_j = np.empty(x.shape) # (..., xdim) + + # Compute univariate Lagrange polynomials in each dimension + for n in range(self.xdim()): + x_n = x[..., n, np.newaxis] # (..., 1) + x_j = self.x_grids[n] # (Nx,) + w_j = self.weights[n] # (Nx,) + + # Compute the jth Lagrange basis polynomial L_j(x_n) for this x dimension (in barycentric form) + c = x_n - x_j + div_zero_idx = np.abs(c) <= 1e-4 * np.abs(x_j) + 1e-8 # Track where x is at an interpolation pnt x_j + c[div_zero_idx] = 1 # Temporarily set to 1 to avoid divide by zero error + c = w_j / c + L_j[..., n] = c[..., j[n]] / np.sum(c, axis=-1) # (...) same size as original x + + # Set L_j(x==x_j)=1 for the current j and set L_j(x==x_j)=0 for x_j = x_i, i != j + L_j[div_zero_idx[..., j[n]], n] = 1 + L_j[np.any(div_zero_idx[..., [idx for idx in range(grid_sizes[n]) if idx != j[n]]], axis=-1), n] = 0 + + # Add multivariate basis polynomial contribution to interpolation output + L_j = np.prod(L_j, axis=-1, keepdims=True) # (..., 1) + y += L_j * yi[i, :] + + return y + + @staticmethod + def get_grid_sizes(beta: tuple, k: int = 2) -> list[int]: + """Compute number of grid points in each input dimension. + + :param beta: refinement level indices + :param k: level-to-grid-size multiplier (probably just always `k=2`) + :returns: list of grid sizes in each dimension + """ + return [k*beta[i] + 1 for i in range(len(beta))] + + @staticmethod + def leja_1d(N: int, z_bds: tuple, z_pts: np.ndarray = None, wt_fcn: callable = None) -> np.ndarray: + """Find the next `N` points in the Leja sequence of `z_pts`. + + :param N: number of new points to add to the sequence + :param z_bds: bounds on the 1d domain + :param z_pts: current univariate Leja sequence `(Nz,)`, start at middle of `z_bds` if `None` + :param wt_fcn: weighting function, uses a constant weight if `None`, callable as `wt_fcn(z)` + :returns: the Leja sequence `z_pts` augmented by `N` new points + """ + # if wt_fcn is None: + wt_fcn = lambda z: 1 # UPDATE: ignore RV weighting, unbounded pdfs like Gaussian cause problems + if z_pts is None: + z_pts = (z_bds[1] + z_bds[0]) / 2 + N = N - 1 + z_pts = np.atleast_1d(z_pts).astype(np.float32) + + # Construct Leja sequence by maximizing the Leja objective sequentially + for i in range(N): + obj_fun = lambda z: -wt_fcn(np.array(z).astype(np.float32)) * np.prod(np.abs(z - z_pts)) + res = direct(obj_fun, [z_bds]) # Use global DIRECT optimization over 1d domain + z_star = res.x + z_pts = np.concatenate((z_pts, z_star)) + + return z_pts diff --git a/src/amisc/rv.py b/src/amisc/rv.py new file mode 100644 index 0000000..8401a4e --- /dev/null +++ b/src/amisc/rv.py @@ -0,0 +1,340 @@ +"""`rv.py` + +Provides small classes for random variables. + +Includes +-------- +- `BaseRV`: Abstract wrapper class of a random variable +- `UniformRV`: a uniformly-distributed random variable +- `NormalRV`: a normally-distributed random variable +- `ScalarRV`: a stand-in class for a variable with no uncertainty or pdf +- `LogUniformRV`: base 10 log-uniform +- `LogNormalRV`: base 10 log-normal +""" +from abc import ABC, abstractmethod +import random +import string + +import numpy as np + + +class BaseRV(ABC): + """Small wrapper class similar to `scipy.stats` random variables (RVs). + + :ivar id: an identifier for the variable + :ivar bds: the explicit domain bounds of the variable (limits of where you expect to use it) + :ivar nominal: a typical value for this variable (within `bds`) + :ivar tex: latex format for the random variable, i.e. r"$x_i$" + :ivar description: a lengthier description of the variable + :ivar units: assumed units for the variable (if applicable) + :ivar param_type: an additional descriptor for how this rv is used, e.g. calibration, operating, design, etc. + + :vartype id: str + :vartype bds: tuple[float, float] + :vartype nominal: float + :vartype tex: str + :vartype description: str + :vartype units: str + :vartype param_type: str + """ + + def __init__(self, id='', *, tex='', description='Random variable', units='-', + param_type='calibration', nominal=1, domain=(0, 1)): + """Child classes must implement `sample` and `pdf` methods.""" + self.bds = tuple(domain) + self.nominal = nominal + self.id = id if id != '' else 'X_' + ''.join(random.choices(string.ascii_uppercase + string.digits, k=4)) + self.is_custom_id = id != '' # Whether a custom id was assigned to this variable + self.tex = tex + self.description = description + self.units = units + self.param_type = param_type + + def __repr__(self): + return r'{}'.format(f"{self.id} - {self.description} ({self.units})") + + def __str__(self): + return self.id + + def __eq__(self, other): + """Consider two RVs equal if they share the same string id. + + Also returns true when checking if this RV is equal to a string id by itself. + """ + if isinstance(other, BaseRV): + return self.id == other.id + elif isinstance(other, str): + return self.id == other + else: + return NotImplemented + + def __hash__(self): + return hash(self.id) + + def to_tex(self, units=False, symbol=True): + """Return a raw string that is well-formatted for plotting (with tex). + + :param units: whether to include the units in the string + :param symbol: just latex symbol if true, otherwise the full description + """ + s = self.tex if symbol else self.description + if s == '': + s = str(self) + return r'{} [{}]'.format(s, self.units) if units else r'{}'.format(s) + + def bounds(self): + """Return a tuple of the defined domain of this RV.""" + return self.bds + + def update_bounds(self, lb, ub): + """Update the defined domain of this RV to `(lb, ub)`.""" + self.bds = (lb, ub) + + def sample_domain(self, shape: tuple | int) -> np.ndarray: + """Return an array of the given `shape` for random samples over the domain of this RV. + + :param shape: the shape of samples to return + :returns samples: random samples over the domain of the random variable + """ + if isinstance(shape, int): + shape = (shape, ) + return np.random.rand(*shape) * (self.bds[1] - self.bds[0]) + self.bds[0] + + @abstractmethod + def pdf(self, x: np.ndarray) -> np.ndarray: + """Compute the PDF of the RV at the given `x` locations. + + :param x: locations to compute the PDF at + :returns f: the PDF evaluations at `x` + """ + pass + + @abstractmethod + def sample(self, shape: tuple | int, nominal: float = None) -> np.ndarray: + """Draw samples from the PDF. + + :param shape: the shape of the returned samples + :param nominal: a nominal value to use if applicable (i.e. a center for relative Uniform or Normal) + :returns: samples from the PDF of this random variable + """ + pass + + +class ScalarRV(BaseRV): + """A stand-in variable with no uncertainty/pdf, just scalars.""" + + def pdf(self, x): + return np.ones(x.shape) + + def sample(self, shape, nominal=None): + if isinstance(shape, int): + shape = (shape, ) + if nominal is not None: + return np.ones(shape)*nominal + else: + return self.sample_domain(shape) + + +class UniformRV(BaseRV): + """A uniformly-distributed random variable. + + Can be uniformly distributed in one of three ways: between global bounds, relative within a percent, or relative + within a set absolute tolerance. + + :ivar type: specifies the type of uniform distribution, either 'bds', 'pct', or 'tol' as described above + :ivar value: the absolute tolerance or percent uncertainty if type is 'tol' or 'pct' + + :vartype type: str + :vartype value: float + """ + + def __init__(self, arg1: float, arg2: float | str, id='', **kwargs): + """Construct a uniformly-distributed random variable. + + :param arg1: lower bound if specifying U(lb, ub), otherwise a tol or pct if specifying U(+/- tol/pct) + :param arg2: upper bound if specifying U(lb, ub), otherwise a str of either 'tol' or 'pct' + """ + domain = kwargs.get('domain', None) + if isinstance(arg2, str): + self.value = arg1 + self.type = arg2 + else: + self.value = None + self.type = 'bds' + if self.type == 'bds': + domain = (arg1, arg2) if domain is None else tuple(domain) # This means domain overrides (arg1, arg2) + else: + domain = (0, 1) if domain is None else tuple(domain) + kwargs['domain'] = domain + super().__init__(id, **kwargs) + + # Set default nominal value as middle of the domain if not specified + if kwargs.get('nominal', None) is None: + self.nominal = (self.bds[1] + self.bds[0]) / 2 + + def __str__(self): + return self.id if self.is_custom_id else f'U({self.bds[0]}, {self.bds[1]})' + + def get_uniform_bounds(self, nominal: float = None) -> tuple[float, float]: + """Return the correct set of bounds based on type of uniform distribution. + + :param nominal: the center value for relative uniform distributions + :returns: the uniform bounds + """ + match self.type: + case 'bds': + return self.bds + case 'pct': + if nominal is None: + return self.bds + return nominal * (1 - self.value), nominal * (1 + self.value) + case 'tol': + if nominal is None: + return self.bds + return nominal - self.value, nominal + self.value + case other: + raise NotImplementedError(f'self.type = {self.type} not known. Choose from ["pct, "tol", "bds"]') + + def pdf(self, x: np.ndarray, nominal: float = None) -> np.ndarray: + """Compute the pdf for a uniform distribution. + + :param x: locations to compute the pdf at + :param nominal: center location for relative uniform rvs + :returns: the evaluated PDF at `x` + """ + bds = self.get_uniform_bounds(nominal) + den = bds[1] - bds[0] + den = 1 if np.isclose(den, 0) else den + y = np.broadcast_to(1 / den, x.shape).copy() + y[np.where(x > bds[1])] = 0 + y[np.where(x < bds[0])] = 0 + return y + + def sample(self, shape, nominal=None): + if isinstance(shape, int): + shape = (shape, ) + bds = self.get_uniform_bounds(nominal) + return np.random.rand(*shape) * (bds[1] - bds[0]) + bds[0] + + +class LogUniformRV(BaseRV): + """A base 10 log-uniform distributed random variable, only supports absolute bounds.""" + + def __init__(self, log10_a: float, log10_b: float, id='', **kwargs): + """Construct the log-uniform random variable. + + :param log10_a: the lower bound in log10 space + :param log10_b: the upper bound in log10 space + """ + super().__init__(id, **kwargs) + self.bds = (10**log10_a, 10**log10_b) + + def __str__(self): + return self.id if self.is_custom_id else f'LU({np.log10(self.bds[0]): .2f}, {np.log10(self.bds[1]): .2f})' + + def pdf(self, x): + return np.log10(np.e) / (x * (np.log10(self.bds[1]) - np.log10(self.bds[0]))) + + def sample(self, shape, nominal=None): + if isinstance(shape, int): + shape = (shape, ) + lb = np.log10(self.bds[0]) + ub = np.log10(self.bds[1]) + return 10 ** (np.random.rand(*shape) * (ub - lb) + lb) + + +class LogNormalRV(BaseRV): + """A base 10 log-normal distributed random variable. + + :ivar mu: the center of the log-normal distribution + :ivar std: the standard deviation of the log-normal distribution + + :vartype mu: float + :vartype std: float + """ + + def __init__(self, mu: float, std: float, id='', **kwargs): + """Construct the RV with the mean and std of the underlying distribution, + i.e. $\\log_{10}(x) \\sim N(\\mu, \\sigma)$. + + :param mu: the center of the log-normal distribution + :param std: the standard deviation of the log-normal distribution + """ + domain = kwargs.get('domain', None) + if domain is None: + domain = (10 ** (mu - 3*std), 10 ** (mu + 3*std)) # Use a default domain of +- 3std + kwargs['domain'] = domain + super().__init__(id, **kwargs) + self.std = std + self.mu = mu + + def recenter(self, mu: float, std: float = None): + """Move the center of the distribution to `mu` with standard deviation `std` (optional) + + :param mu: the new center of the distribution + :param std: (optional) new standard deviation + """ + self.mu = mu + if std is not None: + self.std = std + + def __str__(self): + return self.id if self.is_custom_id else f'LN_10({self.mu}, {self.std})' + + def pdf(self, x): + return (np.log10(np.e) / (x * self.std * np.sqrt(2 * np.pi))) * \ + np.exp(-0.5 * ((np.log10(x) - self.mu) / self.std) ** 2) + + def sample(self, shape, nominal=None): + if isinstance(shape, int): + shape = (shape, ) + scale = np.log10(np.e) + center = self.mu if nominal is None else nominal + return np.random.lognormal(mean=(1 / scale) * center, sigma=(1 / scale) * self.std, size=shape) + # return 10 ** (np.random.randn(*size)*self.std + center) # Alternatively + + +class NormalRV(BaseRV): + """A normally-distributed random variable. + + :ivar mu: float, the mean of the normal distribution + :ivar std: float, the standard deviation of the normal distribution + + :vartype mu: float + :vartype std: float + """ + + def __init__(self, mu, std, id='', **kwargs): + domain = kwargs.get('domain', None) + if domain is None: + domain = (mu - 2.5*std, mu + 2.5*std) # Use a default domain of +- 2.5std + kwargs['domain'] = domain + super().__init__(id, **kwargs) + self.mu = mu + self.std = std + + # Set default nominal value as the provided mean + if kwargs.get('nominal', None) is None: + self.nominal = mu + + def recenter(self, mu: float, std: float =None): + """Move the center of the distribution to `mu` with standard deviation `std` (optional) + + :param mu: the new center of the distribution + :param std: (optional) new standard deviation + """ + self.mu = mu + if std is not None: + self.std = std + + def __str__(self): + return self.id if self.is_custom_id else f'N({self.mu}, {self.std})' + + def pdf(self, x): + return (1 / (np.sqrt(2 * np.pi) * self.std)) * np.exp(-0.5 * ((x - self.mu) / self.std) ** 2) + + def sample(self, shape, nominal=None): + if isinstance(shape, int): + shape = (shape, ) + center = self.mu if nominal is None else nominal + return np.random.randn(*shape) * self.std + center diff --git a/src/amisc/system.py b/src/amisc/system.py new file mode 100644 index 0000000..64cd5d8 --- /dev/null +++ b/src/amisc/system.py @@ -0,0 +1,1377 @@ +"""`system.py` + +The `SystemSurrogate` is a framework for multidisciplinary models. It manages multiple single discipline component +models and the connections between them. It provides a top-level interface for constructing and evaluating surrogates. + +Features +-------- +- Manages multidisciplinary models in a graph data structure, supports feedforward and feedback connections +- Feedback connections are solved with a fixed-point iteration (FPI) nonlinear solver +- FPI uses Anderson acceleration and surrogate evaluations for speed-up +- Top-level interface for training and using surrogates of each component model +- Adaptive experimental design for choosing training data efficiently +- Convenient testing, plotting, and performance metrics provided to assess quality of surrogates +- Detailed logging and traceback information +- Supports parallel execution with OpenMP and MPI protocols +- Abstract and flexible interfacing with component models + +!!! Info "Model specification" + Models are callable Python wrapper functions of the form `ret = model(x, *args, **kwargs)`, where `x` is an + `np.ndarray` of model inputs (and `*args, **kwargs` allow passing any other required configurations for your model). + The return value is a Python dictionary of the form `ret = {'y': y, 'files': files, 'cost': cost, etc.}`. + In the return dictionary, you specify the raw model output `y` as an `np.ndarray` at a _minimum_. Optionally, you can + specify paths to output files and the average model cost (in seconds of cpu time), and anything else you want. Your + `model()` function can do anything it wants in order to go from `x` → `y`. Python has the flexibility to call + virtually any external codes, or to implement the function natively with `numpy`. + +!!! Info "Component specification" + A component adds some extra configuration around a callable `model`. These configurations are defined in a Python + dictionary, which we give the custom type `ComponentSpec`. At a bare _minimum_, you must specify a callable + `model` and its connections to other models within the multidisciplinary system. The limiting case is a single + component model, for which the configuration is simply `component = ComponentSpec(model)`. +""" +import os +import time +import datetime +import functools +import copy +from datetime import timezone +from pathlib import Path +import random +import string +import pickle +from collections import UserDict +from concurrent.futures import Executor + +import numpy as np +import networkx as nx +import dill +import matplotlib.pyplot as plt +from joblib import Parallel, delayed +from joblib.externals.loky import set_loky_pickler + +from amisc.component import SparseGridSurrogate, ComponentSurrogate, AnalyticalSurrogate +from amisc import IndicesRV, IndexSet +from amisc.utils import ax_default, get_logger +from amisc.rv import BaseRV + + +class ComponentSpec(UserDict): + """Provides a simple extension class of a Python dictionary, used to configure a component model. + + !!! Info "Specifying a list of random variables" + The three fields: `exo_in`, `coupling_in`, and `coupling_out` fully determine how a component fits within a + multidisciplinary system. For each, you must specify a list of variables in the same order as the model uses + them. The model will use all exogenous inputs first, and then all coupling inputs. You can use a variable's + global integer index into the system `exo_vars` or `coupling_vars`, or you can use the `str` id of the variable + or the variable itself. This is summarized in the `amisc.IndicesRV` custom type. + + !!! Example + Let's say you have a model: + ```python + def my_model(x, *args, **kwargs): + print(x.shape) # (3,), so a total of 3 inputs + G = 6.674e-11 + m1 = x[0] # System-level input + m2 = x[1] # System-level input + r = x[2] # Coupling input + F = G*m1*m2 / r**2 + return {'y': F} + ``` + Let's say this model is part of a larger system where `m1` and `m2` are specified by the system, and `r` comes + from a different model that predicts the distance between two objects. You would set the configuration as: + ```python + component = ComponentSpec(my_model, exo_in=['m1', 'm2'], coupling_in=['r'], coupling_out=['F']) + ``` + """ + Options = ['model', 'name', 'exo_in', 'coupling_in', 'coupling_out', 'truth_alpha', 'max_alpha', 'max_beta', + 'surrogate', 'model_args', 'model_kwargs', 'save_output'] + + def __init__(self, model: callable, name: str = '', exo_in: IndicesRV = None, + coupling_in: IndicesRV | dict[str: IndicesRV] = None, coupling_out: IndicesRV = None, + truth_alpha: tuple | int = (), max_alpha: tuple | int = (), max_beta: tuple | int = (), + surrogate: str | ComponentSurrogate = 'lagrange', model_args: tuple = (), model_kwargs: dict = None, + save_output: bool = False): + """Construct the configuration for this component model. + + !!! Warning + Always specify the model at a _global_ scope, i.e. don't use `lambda` or nested functions. When saving to + file, only a symbolic reference to the function signature will be saved, which must be globally defined + when loading back from that save file. + + :param model: the component model, must be defined in a global scope (i.e. in a module or top-level of a script) + :param name: the name used to identify this component model + :param exo_in: specifies the global, system-level (i.e. exogenous/external) inputs to this model + :param coupling_in: specifies the coupling inputs received from other models + :param coupling_out: specifies all outputs of this model (which may couple later to downstream models) + :param truth_alpha: the model fidelity indices to treat as a "ground truth" reference + :param max_alpha: the maximum model fidelity indices to allow for refinement purposes + :param max_beta: the maximum surrogate fidelity indices to allow for refinement purposes + :param surrogate: one of ('lagrange, 'analytical'), or the `ComponentSurrogate` class to use directly + :param model_args: optional arguments to pass to the component model + :param model_kwargs: optional keyword arguments to pass to the component model + :param save_output: whether this model will be saving outputs to file + """ + d = locals() + d2 = {key: value for key, value in d.items() if key in ComponentSpec.Options} + super().__init__(d2) + + def __setitem__(self, key, value): + if key in ComponentSpec.Options: + super().__setitem__(key, value) + else: + raise ValueError(f'"{key}" is not applicable for a ComponentSpec. Try one of {ComponentSpec.Options}.') + + def __delitem__(self, key): + raise TypeError("Not allowed to delete items from a ComponentSpec.") + + +class SystemSurrogate: + """Multidisciplinary (MD) surrogate framework top-level class. + + !!! Note "Accessing individual components" + The `ComponentSurrogate` objects that compose `SystemSurrogate` are internally stored in the `self.graph.nodes` + data structure. You can access them with `get_component(comp_name)`. + + :ivar exo_vars: global list of exogenous/external inputs for the MD system + :ivar coupling_vars: global list of coupling variables for the MD system (including all system-level outputs) + :ivar refine_level: the total number of refinement steps that have been made + :ivar build_metrics: contains data that summarizes surrogate training progress + :ivar root_dir: root directory where all surrogate build products are saved to file + :ivar log_file: log file where all logs are written to by default + :ivar executor: manages parallel execution for the system + :ivar graph: the internal graph data structure of the MD system + + :vartype exo_vars: list[BaseRV] + :vartype coupling_vars: list[BaseRV] + :vartype refine_level: int + :vartype build_metrics: dict + :vartype root_dir: str + :vartype log_file: str + :vartype executor: Executor + :vartype graph: nx.DiGraph + """ + + def __init__(self, components: list[ComponentSpec] | ComponentSpec, exo_vars: list[BaseRV] | BaseRV, + coupling_vars: list[BaseRV] | BaseRV, est_bds: int = 0, save_dir: str | Path = None, + executor: Executor = None, stdout: bool = True, init_surr: bool = True): + """Construct the MD system surrogate. + + !!! Warning + Component models should always use coupling variables in the order they appear in the system-level + `coupling_vars`. + + :param components: list of components in the MD system (using the ComponentSpec class) + :param exo_vars: list of system-level exogenous/external inputs + :param coupling_vars: list of all coupling variables (including all system-level outputs) + :param est_bds: number of samples to estimate coupling variable bounds, do nothing if 0 + :param save_dir: root directory for all build products (.log, .pkl, .json, etc.), won't save if None + :param executor: an instance of a `concurrent.futures.Executor`, use to iterate new candidates in parallel + :param stdout: whether to log to console + :param init_surr: whether to initialize the surrogate immediately when constructing + """ + # Setup root save directory + if save_dir is not None: + timestamp = datetime.datetime.now(tz=timezone.utc).isoformat().split('.')[0].replace(':', '.') + save_dir = Path(save_dir) / ('amisc_' + timestamp) + os.mkdir(save_dir) + self.root_dir = str(save_dir.resolve()) + os.mkdir(Path(self.root_dir) / 'sys') + os.mkdir(Path(self.root_dir) / 'components') + fname = timestamp + 'UTC_sys.log' + self.log_file = str((Path(self.root_dir) / fname).resolve()) + else: + self.root_dir = None + self.log_file = None + self.logger = get_logger(self.__class__.__name__, log_file=self.log_file, stdout=stdout) + self.executor = executor + + # Store system info in a directed graph data structure + self.graph = nx.DiGraph() + self.exo_vars = copy.deepcopy(exo_vars) if isinstance(exo_vars, list) else [exo_vars] + self.coupling_vars = copy.deepcopy(coupling_vars) if isinstance(coupling_vars, list) else [coupling_vars] + self.refine_level = 0 + self.build_metrics = dict() # Save refinement error metrics during training + + # Construct graph nodes + components = [components] if not isinstance(components, list) else components + for k, comp in enumerate(components): + if comp['name'] == '': + comp['name'] = f'Component {k}' + Nk = len(components) + nodes = {comp['name']: comp for comp in components} # work-around since self.graph.nodes is not built yet + for k in range(Nk): + # Add the component as a str() node, with attributes specifying details of the surrogate + comp_dict = components[k] + indices, surr = self._build_component(comp_dict, nodes=nodes) + self.graph.add_node(comp_dict['name'], surrogate=surr, is_computed=False, **indices) + + # Connect all neighbor nodes + for node, node_obj in self.graph.nodes.items(): + for neighbor in node_obj['local_in']: + self.graph.add_edge(neighbor, node) + + # Estimate coupling variable bounds + if est_bds > 0: + self._estimate_coupling_bds(est_bds) + + # Init system with most coarse fidelity indices in each component + if init_surr: + self.init_system() + self.save_to_file('sys_init.pkl') + + def _build_component(self, component: ComponentSpec, nodes=None) -> tuple[dict, ComponentSurrogate]: + """Build and return a `ComponentSurrogate` from a `dict` that describes the component model/connections. + + :param component: specifies details of a component (see `ComponentSpec`) + :param nodes: `dict` of `{node: node_attributes}`, defaults to `self.graph.nodes` + :returns: `connections, surr`: a `dict` of all connection indices and the `ComponentSurrogate` object + """ + nodes = self.graph.nodes if nodes is None else nodes + kwargs = component.get('model_kwargs', {}) + kwargs = {} if kwargs is None else kwargs + + # Set up defaults if this is a trivial one component system + exo_in = component.get('exo_in', None) + coupling_in = component.get('coupling_in', None) + coupling_out = component.get('coupling_out', None) + if len(nodes) == 1: + exo_in = list(np.arange(0, len(self.exo_vars))) + coupling_in = [] + coupling_out = list(np.arange(0, len(self.coupling_vars))) + else: + exo_in = [] if exo_in is None else exo_in + coupling_in = [] if coupling_in is None else coupling_in + coupling_out = [] if coupling_out is None else coupling_out + exo_in = [exo_in] if not isinstance(exo_in, list) else exo_in + coupling_in = [coupling_in] if not isinstance(coupling_in, list | dict) else coupling_in + component['coupling_out'] = [coupling_out] if not isinstance(coupling_out, list) else coupling_out + + # Raise an error if all inputs or all outputs are empty + if len(exo_in) + len(coupling_in) == 0: + raise ValueError(f'Component {component["name"]} has no inputs! Please specify inputs in ' + f'"exo_in" or "coupling_in".') + if len(component['coupling_out']) == 0: + raise ValueError(f'Component {component["name"]} has no outputs! Please specify outputs in ' + f'"coupling_out".') + + # Get exogenous input indices (might already be a list of ints, otherwise convert list of vars to indices) + if len(exo_in) > 0: + if isinstance(exo_in[0], str | BaseRV): + exo_in = [self.exo_vars.index(var) for var in exo_in] + + # Get global coupling output indices for all nodes (convert list of vars to list of indices if necessary) + global_out = {} + for node, node_obj in nodes.items(): + node_use = node_obj if node != component.get('name') else component + coupling_out = node_use.get('coupling_out', None) + coupling_out = [] if coupling_out is None else coupling_out + coupling_out = [coupling_out] if not isinstance(coupling_out, list) else coupling_out + global_out[node] = [self.coupling_vars.index(var) for var in coupling_out] if isinstance( + coupling_out[0], str | BaseRV) else coupling_out + + # Refactor coupling inputs into both local and global index formats + local_in = dict() # e.g. {'Cathode': [0, 1, 2], 'Thruster': [0,], etc...} + global_in = list() # e.g. [0, 2, 4, 5, 6] + if isinstance(coupling_in, dict): + # If already a dict, get local connection indices from each neighbor + for node, connections in coupling_in.items(): + conn_list = [connections] if not isinstance(connections, list) else connections + if isinstance(conn_list[0], str | BaseRV): + global_ind = [self.coupling_vars.index(var) for var in conn_list] + local_in[node] = sorted([global_out[node].index(i) for i in global_ind]) + else: + local_in[node] = sorted(conn_list) + + # Convert to global coupling indices + for node, local_idx in local_in.items(): + global_in.extend([global_out[node][i] for i in local_idx]) + global_in = sorted(global_in) + else: + # Otherwise, convert a list of global indices or vars into a dict of local indices + if len(coupling_in) > 0: + if isinstance(coupling_in[0], str | BaseRV): + coupling_in = [self.coupling_vars.index(var) for var in coupling_in] + global_in = sorted(coupling_in) + for node, node_obj in nodes.items(): + if node != component['name']: + l = list() + for i in global_in: + try: + l.append(global_out[node].index(i)) + except ValueError: + pass + if l: + local_in[node] = sorted(l) + + # Store all connection indices for this component + connections = dict(exo_in=exo_in, local_in=local_in, global_in=global_in, + global_out=global_out.get(component.get('name'))) + + # Set up a component output save directory + if component.get('save_output', False) and self.root_dir is not None: + output_dir = str((Path(self.root_dir) / 'components' / component['name']).resolve()) + if not Path(output_dir).is_dir(): + os.mkdir(output_dir) + kwargs['output_dir'] = output_dir + else: + if kwargs.get('output_dir', None) is not None: + kwargs['output_dir'] = None + + # Initialize a new component surrogate + surr_class = component.get('surrogate', 'lagrange') + if isinstance(surr_class, str): + match surr_class: + case 'lagrange': + surr_class = SparseGridSurrogate + case 'analytical': + surr_class = AnalyticalSurrogate + case other: + raise NotImplementedError(f"Surrogate type '{surr_class}' is not known at this time.") + + # Check for an override of model fidelity indices (to enable just single-fidelity evaluation) + if kwargs.get('hf_override', False): + truth_alpha, max_alpha = (), () + kwargs['hf_override'] = component['truth_alpha'] # Pass in the truth alpha indices as a kwarg to model + else: + truth_alpha, max_alpha = component['truth_alpha'], component['max_alpha'] + max_beta = component.get('max_beta') + truth_alpha = (truth_alpha,) if isinstance(truth_alpha, int) else truth_alpha + max_alpha = (max_alpha,) if isinstance(max_alpha, int) else max_alpha + max_beta = (max_beta,) if isinstance(max_beta, int) else max_beta + + # Assumes input ordering is exogenous vars + sorted coupling vars + x_vars = [self.exo_vars[i] for i in exo_in] + [self.coupling_vars[i] for i in global_in] + surr = surr_class(x_vars, component['model'], truth_alpha=truth_alpha, max_alpha=max_alpha, + max_beta=max_beta, executor=self.executor, log_file=self.log_file, + model_args=component.get('model_args'), model_kwargs=kwargs) + return connections, surr + + def swap_component(self, component: ComponentSpec, exo_add: BaseRV | list[BaseRV] = None, + exo_remove: IndicesRV = None, qoi_add: BaseRV | list[BaseRV] = None, + qoi_remove: IndicesRV = None): + """Swap a new component into the system, updating all connections/inputs. + + !!! Warning "Beta feature, proceed with caution" + If you are swapping a new component in, you cannot remove any inputs that are expected by other components, + including the coupling variables output by the current model. + + :param component: specs of new component model (must replace an existing component with matching `name`) + :param exo_add: variables to add to system exogenous inputs (will be appended to end) + :param exo_remove: indices of system exogenous inputs to delete (can't be shared by other components) + :param qoi_add: system output QoIs to add (will be appended to end of `coupling_vars`) + :param qoi_remove: indices of system `coupling_vars` to delete (can't be shared by other components) + """ + # Delete system exogenous inputs + if exo_remove is None: + exo_remove = [] + exo_remove = [exo_remove] if not isinstance(exo_remove, list) else exo_remove + exo_remove = [self.exo_vars.index(var) for var in exo_remove] if exo_remove and isinstance( + exo_remove[0], str | BaseRV) else exo_remove + + exo_remove = sorted(exo_remove) + for j, exo_var_idx in enumerate(exo_remove): + # Adjust exogenous indices for all components to account for deleted system inputs + for node, node_obj in self.graph.nodes.items(): + if node != component['name']: + for i, idx in enumerate(node_obj['exo_in']): + if idx == exo_var_idx: + error_msg = f"Can't delete system exogenous input at idx {exo_var_idx}, since it is " \ + f"shared by component '{node}'." + self.logger.error(error_msg) + raise ValueError(error_msg) + if idx > exo_var_idx: + node_obj['exo_in'][i] -= 1 + + # Need to update the remaining delete indices by -1 to account for each sequential deletion + del self.exo_vars[exo_var_idx] + for i in range(j+1, len(exo_remove)): + exo_remove[i] -= 1 + + # Append any new exogenous inputs to the end + if exo_add is not None: + exo_add = [exo_add] if not isinstance(exo_add, list) else exo_add + self.exo_vars.extend(exo_add) + + # Delete system qoi outputs (if not shared by other components) + qoi_remove = sorted(self._get_qoi_ind(qoi_remove)) + for j, qoi_idx in enumerate(qoi_remove): + # Adjust coupling indices for all components to account for deleted system outputs + for node, node_obj in self.graph.nodes.items(): + if node != component['name']: + for i, idx in enumerate(node_obj['global_in']): + if idx == qoi_idx: + error_msg = f"Can't delete system QoI at idx {qoi_idx}, since it is an input to " \ + f"component '{node}'." + self.logger.error(error_msg) + raise ValueError(error_msg) + if idx > qoi_idx: + node_obj['global_in'][i] -= 1 + + for i, idx in enumerate(node_obj['global_out']): + if idx > qoi_idx: + node_obj['global_out'][i] -= 1 + + # Need to update the remaining delete indices by -1 to account for each sequential deletion + del self.coupling_vars[qoi_idx] + for i in range(j+1, len(qoi_remove)): + qoi_remove[i] -= 1 + + # Append any new system QoI outputs to the end + if qoi_add is not None: + qoi_add = [qoi_add] if not isinstance(qoi_add, list) else qoi_add + self.coupling_vars.extend(qoi_add) + + # Build and initialize the new component surrogate + indices, surr = self._build_component(component) + surr.init_coarse() + + # Make changes to adj matrix if coupling inputs changed + prev_neighbors = list(self.graph.nodes[component['name']]['local_in'].keys()) + new_neighbors = list(indices['local_in'].keys()) + for neighbor in new_neighbors: + if neighbor not in prev_neighbors: + self.graph.add_edge(neighbor, component['name']) + else: + prev_neighbors.remove(neighbor) + for neighbor in prev_neighbors: + self.graph.remove_edge(neighbor, component['name']) + + self.logger.info(f"Swapped component '{component['name']}'.") + nx.set_node_attributes(self.graph, {component['name']: {'exo_in': indices['exo_in'], 'local_in': + indices['local_in'], 'global_in': indices['global_in'], + 'global_out': indices['global_out'], + 'surrogate': surr, 'is_computed': False}}) + + def insert_component(self, component: ComponentSpec, exo_add: BaseRV | list[BaseRV] = None, + qoi_add: BaseRV | list[BaseRV] = None): + """Insert a new component into the system. + + :param component: specs of new component model + :param exo_add: variables to add to system exogenous inputs (will be appended to end of `exo_vars`) + :param qoi_add: system output QoIs to add (will be appended to end of `coupling_vars`) + """ + if exo_add is not None: + exo_add = [exo_add] if not isinstance(exo_add, list) else exo_add + self.exo_vars.extend(exo_add) + if qoi_add is not None: + qoi_add = [qoi_add] if not isinstance(qoi_add, list) else qoi_add + self.coupling_vars.extend(qoi_add) + + indices, surr = self._build_component(component) + surr.init_coarse() + self.graph.add_node(component['name'], surrogate=surr, is_computed=False, **indices) + + # Add graph edges + neighbors = list(indices['local_in'].keys()) + for neighbor in neighbors: + self.graph.add_edge(neighbor, component['name']) + self.logger.info(f"Inserted component '{component['name']}'.") + + def _save_on_error(func): + """Gracefully exit and save `SystemSurrogate` on any errors.""" + @functools.wraps(func) + def wrap(self, *args, **kwargs): + try: + return func(self, *args, **kwargs) + except: + self.save_to_file('sys_error.pkl') + self.logger.critical(f'An error occurred during execution of {func.__name__}. Saving ' + f'SystemSurrogate object to sys_error.pkl', exc_info=True) + self.logger.info(f'Final system surrogate on exit: \n {self}') + raise + return wrap + _save_on_error = staticmethod(_save_on_error) + + @_save_on_error + def init_system(self): + """Add the coarsest multi-index to each component surrogate.""" + self._print_title_str('Initializing all component surrogates') + for node, node_obj in self.graph.nodes.items(): + node_obj['surrogate'].init_coarse() + # for alpha, beta in list(node_obj['surrogate'].candidate_set): + # # Add one refinement in each input dimension to initialize + # node_obj['surrogate'].activate_index(alpha, beta) + self.logger.info(f"Initialized component '{node}'.") + + @_save_on_error + def fit(self, qoi_ind: IndicesRV = None, num_refine: int = 100, max_iter: int = 20, max_tol: float = 1e-3, + max_runtime: float = 1, save_interval: int = 0, update_bounds: bool = True, test_set: dict = None, + n_jobs: int = 1): + """Train the system surrogate adaptively by iterative refinement until an end condition is met. + + :param qoi_ind: list of system QoI variables to focus refinement on, use all QoI if not specified + :param num_refine: number of samples of exogenous inputs to compute error indicators on + :param max_iter: the maximum number of refinement steps to take + :param max_tol: the max allowable value in relative L2 error to achieve + :param max_runtime: the maximum wall clock time (hr) to run refinement for (will go until all models finish) + :param save_interval: number of refinement steps between each progress save, none if 0 + :param update_bounds: whether to continuously update coupling variable bounds during refinement + :param test_set: `dict(xt=(Nt, x_dim), yt=(Nt, y_dim)` to show convergence of surrogate to the truth model + :param n_jobs: number of cpu workers for computing error indicators (on master MPI task), 1=sequential + """ + qoi_ind = self._get_qoi_ind(qoi_ind) + Nqoi = len(qoi_ind) + max_iter = self.refine_level + max_iter + curr_error = np.inf + t_start = time.time() + test_stats, xt, yt, t_fig, t_ax = None, None, None, None, None + + # Record of (error indicator, component, alpha, beta, num_evals, total added cost (s)) for each iteration + train_record = self.build_metrics.get('train_record', []) + if test_set is not None: + xt, yt = test_set['xt'], test_set['yt'] + xt, yt = self.build_metrics.get('xt', xt), self.build_metrics.get('yt', yt) # Overrides test set param + + # Track convergence progress on a test set and on the max error indicator + err_fig, err_ax = plt.subplots() + if xt is not None and yt is not None: + self.build_metrics['xt'] = xt + self.build_metrics['yt'] = yt + if self.build_metrics.get('test_stats') is not None: + test_stats = self.build_metrics.get('test_stats') + else: + # Get initial perf metrics, (2, Nqoi) + test_stats = np.expand_dims(self.get_test_metrics(xt, yt, qoi_ind=qoi_ind), axis=0) + t_fig, t_ax = plt.subplots(1, Nqoi) if Nqoi > 1 else plt.subplots() + + # Set up a parallel pool of workers, sequential if n_jobs=1 + with Parallel(n_jobs=n_jobs, verbose=0) as ppool: + while True: + # Check all end conditions + if self.refine_level >= max_iter: + self._print_title_str(f'Termination criteria reached: Max iteration {self.refine_level}/{max_iter}') + break + if curr_error == -np.inf: + self._print_title_str(f'Termination criteria reached: No candidates left to refine') + break + if curr_error < max_tol: + self._print_title_str(f'Termination criteria reached: L2 error {curr_error} < tol {max_tol}') + break + if ((time.time() - t_start)/3600.0) >= max_runtime: + actual = datetime.timedelta(seconds=time.time()-t_start) + target = datetime.timedelta(seconds=max_runtime*3600) + self._print_title_str(f'Termination criteria reached: runtime {str(actual)} > {str(target)}') + break + + # Refine surrogate and save progress + refine_res = self.refine(qoi_ind=qoi_ind, num_refine=num_refine, update_bounds=update_bounds, + ppool=ppool) + curr_error = refine_res[0] + if save_interval > 0 and self.refine_level % save_interval == 0: + self.save_to_file(f'sys_iter_{self.refine_level}.pkl') + + # Plot progress of error indicator + train_record.append(refine_res) + error_record = [res[0] for res in train_record] + self.build_metrics['train_record'] = train_record + err_ax.clear(); err_ax.grid(); err_ax.plot(error_record, '-k') + ax_default(err_ax, 'Iteration', r'Relative $L_2$ error indicator', legend=False) + err_ax.set_yscale('log') + if self.root_dir is not None: + err_fig.savefig(str(Path(self.root_dir) / 'error_indicator.png'), dpi=300, format='png') + + # Plot progress on test set + if xt is not None and yt is not None: + stats = self.get_test_metrics(xt, yt, qoi_ind=qoi_ind) + test_stats = np.concatenate((test_stats, stats[np.newaxis, ...]), axis=0) + for i in range(Nqoi): + ax = t_ax if Nqoi == 1 else t_ax[i] + ax.clear(); ax.grid(); ax.set_yscale('log') + ax.plot(test_stats[:, 1, i], '-k') + ax.set_title(self.coupling_vars[qoi_ind[i]].to_tex(units=True)) + ax_default(ax, 'Iteration', r'Relative $L_2$ error', legend=False) + t_fig.set_size_inches(3.5*Nqoi, 3.5) + t_fig.tight_layout() + if self.root_dir is not None: + t_fig.savefig(str(Path(self.root_dir) / 'test_set.png'), dpi=300, format='png') + self.build_metrics['test_stats'] = test_stats + + self.save_to_file(f'sys_final.pkl') + self.logger.info(f'Final system surrogate: \n {self}') + + def get_allocation(self, idx: int = None): + """Get a breakdown of cost allocation up to a certain iteration number during training (starting at 1). + + :param idx: the iteration number to get allocation results for (defaults to last refinement step) + :returns: `cost_alloc, offline_alloc, cost_cum` - the cost alloc per node/fidelity and cumulative training cost + """ + if idx is None: + idx = self.refine_level + if idx > self.refine_level: + raise ValueError(f'Specified index: {idx} is greater than the max training level of {self.refine_level}') + + cost_alloc = dict() # Cost allocation per node and model fidelity + cost_cum = [0.0] # Cumulative cost allocation during training + + # Add initialization costs for each node + for node, node_obj in self.graph.nodes.items(): + surr = node_obj['surrogate'] + base_alpha = (0,) * len(surr.truth_alpha) + base_beta = (0,) * (len(surr.max_refine) - len(surr.truth_alpha)) + base_cost = surr.get_cost(base_alpha, base_beta) + cost_alloc[node] = dict() + if base_cost > 0: + cost_alloc[node][str(base_alpha)] = np.array([1, float(base_cost)]) + cost_cum[0] += float(base_cost) + + # Add cumulative training costs + for i in range(idx): + err_indicator, node, alpha, beta, num_evals, cost = self.build_metrics['train_record'][i] + if cost_alloc[node].get(str(alpha), None) is None: + cost_alloc[node][str(alpha)] = np.zeros(2) # (num model evals, total cpu_time cost) + cost_alloc[node][str(alpha)] += [round(num_evals), float(cost)] + cost_cum.append(float(cost)) + + # Get summary of total offline costs spent building search candidates (i.e. training overhead) + offline_alloc = dict() + for node, node_obj in self.graph.nodes.items(): + surr = node_obj['surrogate'] + offline_alloc[node] = dict() + for alpha, beta in surr.candidate_set: + if offline_alloc[node].get(str(alpha), None) is None: + offline_alloc[node][str(alpha)] = np.zeros(2) # (num model evals, total cpu_time cost) + added_cost = surr.get_cost(alpha, beta) + base_cost = surr.get_sub_surrogate(alpha, beta).model_cost + offline_alloc[node][str(alpha)] += [round(added_cost/base_cost), float(added_cost)] + + return cost_alloc, offline_alloc, np.cumsum(cost_cum) + + def get_test_metrics(self, xt: np.ndarray, yt: np.ndarray, qoi_ind: IndicesRV = None, + training: bool = True) -> np.ndarray: + """Get relative L2 error metric over a test set. + + :param xt: `(Nt, x_dim)` random test set of inputs + :param yt: `(Nt, y_dim)` random test set outputs + :param qoi_ind: list of indices of QoIs to get metrics for + :param training: whether to evaluate the surrogate in training or evaluation mode + :returns: `stats` - `(2, Nqoi)` array → `[num_candidates, rel_L2_error]` for each QoI + """ + qoi_ind = self._get_qoi_ind(qoi_ind) + ysurr = self(xt, training=training) + ysurr = ysurr[:, qoi_ind] + yt = yt[:, qoi_ind] + with np.errstate(divide='ignore', invalid='ignore'): + rel_l2_err = np.sqrt(np.mean((yt - ysurr) ** 2, axis=0)) / np.sqrt(np.mean(yt ** 2, axis=0)) + rel_l2_err = np.nan_to_num(rel_l2_err, posinf=np.nan, neginf=np.nan, nan=np.nan) + num_cands = 0 + for node, node_obj in self.graph.nodes.items(): + num_cands += len(node_obj['surrogate'].index_set) + len(node_obj['surrogate'].candidate_set) + + # Get test stats for each QoI + stats = np.zeros((2, yt.shape[-1])) + self.logger.debug(f'{"QoI idx": >10} {"Iteration": >10} {"len(I_k)": >10} {"Relative L2": >15}') + for i in range(yt.shape[-1]): + stats[:, i] = np.array([num_cands, rel_l2_err[i]]) + self.logger.debug(f'{i: 10d} {self.refine_level: 10d} {num_cands: 10d} {rel_l2_err[i]: 15.5f}') + + return stats + + def _get_qoi_ind(self, qoi_ind: IndicesRV) -> list[int]: + """Small helper to make sure QoI indices are a list of integers.""" + if qoi_ind is None: + qoi_ind = list(np.arange(0, len(self.coupling_vars))) + qoi_ind = [qoi_ind] if not isinstance(qoi_ind, list) else qoi_ind + qoi_ind = [self.coupling_vars.index(var) for var in qoi_ind] if qoi_ind and isinstance( + qoi_ind[0], str | BaseRV) else qoi_ind + + return qoi_ind + + def refine(self, qoi_ind: IndicesRV = None, num_refine: int = 100, update_bounds: bool = True, + ppool: Parallel = None) -> tuple: + """Find and refine the component surrogate with the largest error on system-level QoI. + + :param qoi_ind: indices of system QoI to focus surrogate refinement on, use all QoI if not specified + :param num_refine: number of samples of exogenous inputs to compute error indicators on + :param update_bounds: whether to continuously update coupling variable bounds + :param ppool: a `Parallel` instance from `joblib` to compute error indicators in parallel, None=sequential + :returns refine_res: a tuple of `(error_indicator, component, node_star, alpha_star, beta_star, N, cost)` + indicating the chosen candidate index and incurred cost + """ + self._print_title_str(f'Refining system surrogate: iteration {self.refine_level + 1}') + set_loky_pickler('dill') # Dill can serialize 'self' for parallel workers + temp_exc = self.executor # It can't serialize an executor though, so must save this temporarily + self.set_executor(None) + qoi_ind = self._get_qoi_ind(qoi_ind) + + # Compute entire integrated-surrogate on a random test set for global system QoI error estimation + x_exo = self.sample_inputs((num_refine,)) + y_curr = self(x_exo, training=True) + y_min, y_max = None, None + if update_bounds: + y_min = np.min(y_curr, axis=0, keepdims=True) # (1, ydim) + y_max = np.max(y_curr, axis=0, keepdims=True) # (1, ydim) + + # Find the candidate surrogate with the largest error indicator + error_max, error_indicator = -np.inf, -np.inf + node_star, alpha_star, beta_star, l2_star, cost_star = None, None, None, -np.inf, 0 + for node, node_obj in self.graph.nodes.items(): + self.logger.info(f"Estimating error for component '{node}'...") + candidates = node_obj['surrogate'].candidate_set.copy() + + def compute_error(alpha, beta): + # Helper function for computing error indicators for a given candidate (alpha, beta) + index_set = node_obj['surrogate'].index_set.copy() + index_set.append((alpha, beta)) + y_cand = self(x_exo, training=True, index_set={node: index_set}) + ymin = np.min(y_cand, axis=0, keepdims=True) + ymax = np.max(y_cand, axis=0, keepdims=True) + error = y_cand[:, qoi_ind] - y_curr[:, qoi_ind] + rel_l2 = np.sqrt(np.nanmean(error ** 2, axis=0)) / np.sqrt(np.nanmean(y_cand[:, qoi_ind] ** 2, axis=0)) + rel_l2 = np.nan_to_num(rel_l2, nan=np.nan, posinf=np.nan, neginf=np.nan) + delta_error = np.nanmax(rel_l2) # Max relative L2 error over all system QoIs + delta_work = max(1, node_obj['surrogate'].get_cost(alpha, beta)) # Cpu time (s) + + return ymin, ymax, delta_error, delta_work + + if len(candidates) > 0: + ret = ppool(delayed(compute_error)(alpha, beta) for alpha, beta in candidates) if ppool is not None \ + else [compute_error(alpha, beta) for alpha, beta in candidates] + + for i, (ymin, ymax, d_error, d_work) in enumerate(ret): + if update_bounds: + y_min = np.min(np.concatenate((y_min, ymin), axis=0), axis=0, keepdims=True) + y_max = np.max(np.concatenate((y_max, ymax), axis=0), axis=0, keepdims=True) + alpha, beta = candidates[i] + error_indicator = d_error / d_work + self.logger.info(f"Candidate multi-index: {(alpha, beta)}. Error indicator: " + f"{error_indicator}. L2 error: {d_error}") + + if error_indicator > error_max: + error_max = error_indicator + node_star, alpha_star, beta_star, l2_star, cost_star = node, alpha, beta, d_error, d_work + else: + self.logger.info(f"Component '{node}' has no available candidates left!") + + # Update all coupling variable ranges + if update_bounds: + for i in range(y_curr.shape[-1]): + self._update_coupling_bds(i, (y_min[0, i], y_max[0, i])) + + # Add the chosen multi-index to the chosen component + self.set_executor(temp_exc) + if node_star is not None: + self.logger.info(f"Candidate multi-index {(alpha_star, beta_star)} chosen for component '{node_star}'") + self.graph.nodes[node_star]['surrogate'].activate_index(alpha_star, beta_star) + self.refine_level += 1 + num_evals = round(cost_star / self[node_star].get_sub_surrogate(alpha_star, beta_star).model_cost) + else: + self.logger.info(f"No candidates left for refinement, iteration: {self.refine_level}") + num_evals = 0 + + return l2_star, node_star, alpha_star, beta_star, num_evals, cost_star + + def predict(self, x: np.ndarray | float, max_fpi_iter: int = 100, anderson_mem: int = 10, fpi_tol: float = 1e-10, + use_model: str | tuple | dict = None, model_dir: str | Path = None, verbose: bool = False, + training: bool = False, index_set: dict[str: IndexSet] = None, qoi_ind: IndicesRV = None) -> np.ndarray: + """Evaluate the system surrogate at exogenous inputs `x`. + + !!! Warning + You can use this function to predict outputs for your MD system using the full-order models rather than the + surrogate, by specifying `use_model`. This is convenient because `SystemSurrogate` manages all the + coupled information flow between models automatically. However, it is *highly* recommended to not use + the full model if your system contains feedback loops. The FPI nonlinear solver would be infeasible using + anything more computationally demanding than the surrogate. + + :param x: `(..., x_dim)` the points to get surrogate predictions for + :param max_fpi_iter: the limit on convergence for the fixed-point iteration routine + :param anderson_mem: hyperparameter for tuning the convergence of FPI with anderson acceleration + :param fpi_tol: tolerance limit for convergence of fixed-point iteration + :param use_model: 'best'=highest-fidelity, 'worst'=lowest-fidelity, tuple=specific fidelity, None=surrogate, + specify a `dict` of the above to assign different model fidelities for diff components + :param model_dir: directory to save model outputs if `use_model` is specified + :param verbose: whether to print out iteration progress during execution + :param training: whether to call the system surrogate in training or evaluation mode, ignored if `use_model` + :param index_set: `dict(node=[indices])` to override default index set for a node (only useful for parallel) + :param qoi_ind: list of qoi indices to return, defaults to returning all system `coupling_vars` + :returns y: `(..., y_dim)` the surrogate approximation of the system QoIs + """ + # Allocate space for all system outputs (just save all coupling vars) + x = np.atleast_1d(x) + ydim = len(self.coupling_vars) + y = np.zeros(x.shape[:-1] + (ydim,)) + valid_idx = ~np.isnan(x[..., 0]) # Keep track of valid samples (set to False if FPI fails) + t1 = 0 + output_dir = None + + # Interpret which model fidelities to use for each component (if specified) + if use_model is not None: + if not isinstance(use_model, dict): + use_model = {node: use_model for node in self.graph.nodes} # use same for each component + else: + use_model = {node: None for node in self.graph.nodes} + use_model = {node: use_model.get(node, None) for node in self.graph.nodes} + + # Initialize all components + for node, node_obj in self.graph.nodes.items(): + node_obj['is_computed'] = False + + # Convert system into DAG by grouping strongly-connected-components + dag = nx.condensation(self.graph) + + # Compute component models in topological order + for supernode in nx.topological_sort(dag): + scc = [n for n in dag.nodes[supernode]['members']] + + # Compute single component feedforward output (no FPI needed) + if len(scc) == 1: + if verbose: + self.logger.info(f"Running component '{scc[0]}'...") + t1 = time.time() + + # Gather inputs + node_obj = self.graph.nodes[scc[0]] + exo_inputs = x[..., node_obj['exo_in']] + # for comp_name in node_obj['local_in']: + # assert self.graph.nodes[comp_name]['is_computed'] + coupling_inputs = y[..., node_obj['global_in']] + comp_input = np.concatenate((exo_inputs, coupling_inputs), axis=-1) # (..., xdim) + + # Compute outputs + indices = index_set.get(scc[0], None) if index_set is not None else None + if model_dir is not None: + output_dir = Path(model_dir) / scc[0] + os.mkdir(output_dir) + comp_output = node_obj['surrogate'](comp_input[valid_idx, :], use_model=use_model.get(scc[0]), + model_dir=output_dir, training=training, index_set=indices) + for local_i, global_i in enumerate(node_obj['global_out']): + y[valid_idx, global_i] = comp_output[..., local_i] + node_obj['is_computed'] = True + + if verbose: + self.logger.info(f"Component '{scc[0]}' completed. Runtime: {time.time() - t1} s") + + # Handle FPI for SCCs with more than one component + else: + # Set the initial guess for all coupling vars (middle of domain) + coupling_bds = [rv.bounds() for rv in self.coupling_vars] + x_couple = np.array([(bds[0] + bds[1]) / 2 for bds in coupling_bds]) + x_couple = np.broadcast_to(x_couple, x.shape[:-1] + x_couple.shape).copy() + + adj_nodes = [] + fpi_idx = set() + for node in scc: + for comp_name, local_idx in self.graph.nodes[node]['local_in'].items(): + # Track the global idx of all coupling vars that need FPI + if comp_name in scc: + fpi_idx.update([self.graph.nodes[comp_name]['global_out'][idx] for idx in local_idx]) + + # Override coupling vars from components outside the scc (should already be computed) + if comp_name not in scc and comp_name not in adj_nodes: + # assert self.graph.nodes[comp_name]['is_computed'] + global_idx = self.graph.nodes[comp_name]['global_out'] + x_couple[..., global_idx] = y[..., global_idx] + adj_nodes.append(comp_name) # Only need to do this once for each adj component + x_couple_next = x_couple.copy() + fpi_idx = sorted(fpi_idx) + + # Main FPI loop + if verbose: + self.logger.info(f"Initializing FPI for SCC {scc} ...") + t1 = time.time() + k = 0 + residual_hist = None + x_hist = None + while True: + for node in scc: + # Gather inputs from exogenous and coupling sources + node_obj = self.graph.nodes[node] + exo_inputs = x[..., node_obj['exo_in']] + coupling_inputs = x_couple[..., node_obj['global_in']] + comp_input = np.concatenate((exo_inputs, coupling_inputs), axis=-1) # (..., xdim) + + # Compute component outputs (just don't do this FPI with the real models, please..) + indices = index_set.get(node, None) if index_set is not None else None + comp_output = node_obj['surrogate'](comp_input[valid_idx, :], use_model=use_model.get(node), + model_dir=None, training=training, index_set=indices) + global_idx = node_obj['global_out'] + for local_i, global_i in enumerate(global_idx): + x_couple_next[valid_idx, global_i] = comp_output[..., local_i] + # Can't splice valid_idx with global_idx for some reason, have to loop over global_idx here + + # Compute residual and check end conditions + residual = np.expand_dims(x_couple_next[..., fpi_idx] - x_couple[..., fpi_idx], axis=-1) + max_error = np.max(np.abs(residual[valid_idx, :, :])) + if verbose: + self.logger.info(f'FPI iter: {k}. Max residual: {max_error}. Time: {time.time() - t1} s') + if max_error <= fpi_tol: + if verbose: + self.logger.info(f'FPI converged for SCC {scc} in {k} iterations with {max_error} < tol ' + f'{fpi_tol}. Final time: {time.time() - t1} s') + break + if k >= max_fpi_iter: + self.logger.warning(f'FPI did not converge in {max_fpi_iter} iterations for SCC {scc}: ' + f'{max_error} > tol {fpi_tol}. Some samples will be returned as NaN.') + converged_idx = np.max(np.abs(residual), axis=(-1, -2)) <= fpi_tol + for idx in fpi_idx: + y[~converged_idx, idx] = np.nan + valid_idx = np.logical_and(valid_idx, converged_idx) + break + + # Keep track of residual and x_couple histories + if k == 0: + residual_hist = residual.copy() # (..., xdim, 1) + x_hist = np.expand_dims(x_couple_next[..., fpi_idx], axis=-1) # (..., xdim, 1) + x_couple[:] = x_couple_next[:] + k += 1 + continue # skip anderson accel on first iteration + + # Iterate with anderson acceleration (only iterate on samples that are not yet converged) + converged_idx = np.max(np.abs(residual), axis=(-1, -2)) <= fpi_tol + curr_idx = np.logical_and(valid_idx, ~converged_idx) + residual_hist = np.concatenate((residual_hist, residual), axis=-1) + x_hist = np.concatenate((x_hist, np.expand_dims(x_couple_next[..., fpi_idx], axis=-1)), axis=-1) + mk = min(anderson_mem, k) + Fk = residual_hist[curr_idx, :, k-mk:] # (..., xdim, mk+1) + C = np.ones(Fk.shape[:-2] + (1, mk + 1)) + b = np.zeros(Fk.shape[:-2] + (len(fpi_idx), 1)) + d = np.ones(Fk.shape[:-2] + (1, 1)) + alpha = np.expand_dims(self._constrained_lls(Fk, b, C, d), axis=-3) # (..., 1, mk+1, 1) + x_new = np.squeeze(x_hist[curr_idx, :, np.newaxis, -(mk+1):] @ alpha, axis=(-1, -2)) + for local_i, global_i in enumerate(fpi_idx): + x_couple[curr_idx, global_i] = x_new[..., local_i] + k += 1 + + # Save outputs of each component in SCC after convergence of FPI + for node in scc: + global_idx = self.graph.nodes[node]['global_out'] + for global_i in global_idx: + y[valid_idx, global_i] = x_couple_next[valid_idx, global_i] + self.graph.nodes[node]['is_computed'] = True + + # Return all component outputs (..., Nqoi), samples that didn't converge during FPI are left as np.nan + qoi_ind = self._get_qoi_ind(qoi_ind) + return y[..., qoi_ind] + + def __call__(self, *args, **kwargs): + """Convenience wrapper to allow calling as `ret = SystemSurrogate(x)`.""" + return self.predict(*args, **kwargs) + + def _estimate_coupling_bds(self, num_est: int, anderson_mem: int = 10, fpi_tol: float = 1e-10, + max_fpi_iter: int = 100): + """Estimate and set the coupling variable bounds. + + :param num_est: the number of samples of exogenous inputs to use + :param anderson_mem: FPI hyperparameter (default is usually good) + :param fpi_tol: floating point tolerance for FPI convergence + :param max_fpi_iter: maximum number of FPI iterations + """ + self._print_title_str('Estimating coupling variable bounds') + x = self.sample_inputs((num_est,)) + y = self(x, use_model='best', verbose=True, anderson_mem=anderson_mem, fpi_tol=fpi_tol, + max_fpi_iter=max_fpi_iter) + for i in range(len(self.coupling_vars)): + lb = np.nanmin(y[:, i]) + ub = np.nanmax(y[:, i]) + self._update_coupling_bds(i, (lb, ub), init=True) + + def _update_coupling_bds(self, global_idx: int, bds: tuple, init: bool = False, buffer: float = 0.05): + """Update coupling variable bounds. + + :param global_idx: global index of coupling variable to update + :param bds: new bounds to update the current bounds with + :param init: whether to set new bounds or update existing (default) + :param buffer: fraction of domain length to buffer upper/lower bounds + """ + offset = buffer * (bds[1] - bds[0]) + offset_bds = (bds[0] - offset, bds[1] + offset) + coupling_bds = [rv.bounds() for rv in self.coupling_vars] + new_bds = offset_bds if init else (min(coupling_bds[global_idx][0], offset_bds[0]), + max(coupling_bds[global_idx][1], offset_bds[1])) + self.coupling_vars[global_idx].update_bounds(*new_bds) + + # Iterate over all components and update internal coupling variable bounds + for node_name, node_obj in self.graph.nodes.items(): + if global_idx in node_obj['global_in']: + # Get the local index for this coupling variable within each component's inputs + local_idx = len(node_obj['exo_in']) + node_obj['global_in'].index(global_idx) + node_obj['surrogate'].update_input_bds(local_idx, new_bds) + + def sample_inputs(self, size: tuple | int, comp: str = 'System', use_pdf: bool = False, + nominal: dict[str: float] = None, constants: set[str] = None) -> np.ndarray: + """Return samples of the inputs according to provided options. + + :param size: tuple or integer specifying shape or number of samples to obtain + :param comp: which component to sample inputs for (defaults to full system exogenous inputs) + :param use_pdf: whether to sample from each variable's pdf, defaults to random samples over input domain instead + :param nominal: `dict(var_id=value)` of nominal values for params with relative uncertainty, also can use + to specify constant values for a variable listed in `constants` + :param constants: set of param types to hold constant while sampling (i.e. calibration, design, etc.), + can also put a `var_id` string in here to specify a single variable to hold constant + :returns x: `(*size, x_dim)` samples of the inputs for the given component/system + """ + size = (size, ) if isinstance(size, int) else size + if nominal is None: + nominal = dict() + if constants is None: + constants = set() + x_vars = self.exo_vars if comp == 'System' else self[comp].x_vars + x = np.empty((*size, len(x_vars))) + for i, var in enumerate(x_vars): + # Set a constant value for this variable + if var.param_type in constants or var in constants: + x[..., i] = nominal.get(var, var.nominal) # Defaults to variable's nominal value if not specified + + # Sample from this variable's pdf or randomly within its domain bounds (reject if outside bounds) + else: + lb, ub = var.bounds() + x_sample = var.sample(size, nominal=nominal.get(var, None)) if use_pdf \ + else var.sample_domain(size) + good_idx = (x_sample < ub) & (x_sample > lb) + num_reject = np.sum(~good_idx) + + while num_reject > 0: + new_sample = var.sample((num_reject,), nominal=nominal.get(var, None)) if use_pdf \ + else var.sample_domain((num_reject,)) + x_sample[~good_idx] = new_sample + good_idx = (x_sample < ub) & (x_sample > lb) + num_reject = np.sum(~good_idx) + + x[..., i] = x_sample + + return x + + def plot_slice(self, slice_idx: IndicesRV = None, qoi_idx: IndicesRV = None, show_surr: bool = True, + show_model: list = None, model_dir: str | Path = None, N: int = 50, nominal: dict[str: float] = None, + random_walk: bool = False, from_file: str | Path = None): + """Helper function to plot 1d slices of the surrogate and/or model(s) over the inputs. + + :param slice_idx: list of exogenous input variables or indices to take 1d slices of + :param qoi_idx: list of model output variables or indices to plot 1d slices of + :param show_surr: whether to show the surrogate prediction + :param show_model: also plot model predictions, list() of ['best', 'worst', tuple(alpha), etc.] + :param model_dir: base directory to save model outputs (if specified) + :param N: the number of points to take in the 1d slice + :param nominal: `dict` of `str(var)->nominal` to use as constant value for all non-sliced variables + :param random_walk: whether to slice in a random d-dimensional direction or hold all params const while slicing + :param from_file: path to a .pkl file to load a saved slice from disk + :returns: `fig, ax` with `num_slice` by `num_qoi` subplots + """ + # Manage loading important quantities from file (if provided) + xs, ys_model, ys_surr = None, None, None + if from_file is not None: + with open(Path(from_file), 'rb') as fd: + slice_data = pickle.load(fd) + slice_idx = slice_data['slice_idx'] # Must use same input slices as save file + show_model = slice_data['show_model'] # Must use same model data as save file + qoi_idx = slice_data['qoi_idx'] if qoi_idx is None else qoi_idx + xs = slice_data['xs'] + model_dir = None # Don't run or save any models if loading from file + + # Set default values (take up to the first 3 slices by default) + rand_id = ''.join(random.choices(string.ascii_uppercase + string.digits, k=4)) + if model_dir is not None: + os.mkdir(Path(model_dir) / f'sweep_{rand_id}') + if nominal is None: + nominal = dict() + slice_idx = list(np.arange(0, min(3, len(self.exo_vars)))) if slice_idx is None else slice_idx + qoi_idx = list(np.arange(0, min(3, len(self.coupling_vars)))) if qoi_idx is None else qoi_idx + if isinstance(slice_idx[0], str | BaseRV): + slice_idx = [self.exo_vars.index(var) for var in slice_idx] + if isinstance(qoi_idx[0], str | BaseRV): + qoi_idx = [self.coupling_vars.index(var) for var in qoi_idx] + + exo_bds = [var.bounds() for var in self.exo_vars] + xlabels = [self.exo_vars[idx].to_tex(units=True) for idx in slice_idx] + ylabels = [self.coupling_vars[idx].to_tex(units=True) for idx in qoi_idx] + + # Construct slice model inputs (if not provided) + if xs is None: + xs = np.zeros((N, len(slice_idx), len(self.exo_vars))) + for i in range(len(slice_idx)): + if random_walk: + # Make a random straight-line walk across d-cube + r0 = np.squeeze(self.sample_inputs((1,), use_pdf=False), axis=0) + r0[slice_idx[i]] = exo_bds[slice_idx[i]][0] # Start slice at this lower bound + rf = np.squeeze(self.sample_inputs((1,), use_pdf=False), axis=0) + rf[slice_idx[i]] = exo_bds[slice_idx[i]][1] # Slice up to this upper bound + xs[0, i, :] = r0 + for k in range(1, N): + xs[k, i, :] = xs[k-1, i, :] + (rf-r0)/(N-1) + else: + # Otherwise, only slice one variable + for j in range(len(self.exo_vars)): + if j == slice_idx[i]: + xs[:, i, j] = np.linspace(exo_bds[slice_idx[i]][0], exo_bds[slice_idx[i]][1], N) + else: + xs[:, i, j] = nominal.get(self.exo_vars[j], self.exo_vars[j].nominal) + + # Walk through each model that is requested by show_model + if show_model is not None: + if from_file is not None: + ys_model = slice_data['ys_model'] + else: + ys_model = list() + for model in show_model: + output_dir = None + if model_dir is not None: + output_dir = (Path(model_dir) / f'sweep_{rand_id}' / + str(model).replace('{', '').replace('}', '').replace(':', '=').replace("'", '')) + os.mkdir(output_dir) + ys_model.append(self(xs, use_model=model, model_dir=output_dir)) + if show_surr: + ys_surr = self(xs) if from_file is None else slice_data['ys_surr'] + + # Make len(qoi) by len(inputs) grid of subplots + fig, axs = plt.subplots(len(qoi_idx), len(slice_idx), sharex='col', sharey='row') + for i in range(len(qoi_idx)): + for j in range(len(slice_idx)): + if len(qoi_idx) == 1: + ax = axs if len(slice_idx) == 1 else axs[j] + elif len(slice_idx) == 1: + ax = axs if len(qoi_idx) == 1 else axs[i] + else: + ax = axs[i, j] + x = xs[:, j, slice_idx[j]] + if show_model is not None: + c = np.array([[0, 0, 0, 1], [0.5, 0.5, 0.5, 1]]) if len(show_model) <= 2 else ( + plt.get_cmap('jet')(np.linspace(0, 1, len(show_model)))) + for k in range(len(show_model)): + model_str = str(show_model[k]).replace('{', '').replace('}', '').replace(':', '=').replace("'", '') + model_ret = ys_model[k] + y_model = model_ret[:, j, qoi_idx[i]] + label = {'best': 'High-fidelity' if len(show_model) > 1 else 'Model', + 'worst': 'Low-fidelity'}.get(model_str, model_str) + ax.plot(x, y_model, ls='-', c=c[k, :], label=label) + if show_surr: + y_surr = ys_surr[:, j, qoi_idx[i]] + ax.plot(x, y_surr, '--r', label='Surrogate') + ylabel = ylabels[i] if j == 0 else '' + xlabel = xlabels[j] if i == len(qoi_idx) - 1 else '' + legend = (i == 0 and j == len(slice_idx) - 1) + ax_default(ax, xlabel, ylabel, legend=legend) + fig.set_size_inches(3 * len(slice_idx), 3 * len(qoi_idx)) + fig.tight_layout() + + # Save results (unless we were already loading from a save file) + if from_file is None and self.root_dir is not None: + fname = f's{",".join([str(i) for i in slice_idx])}_q{",".join([str(i) for i in qoi_idx])}' + fname = f'sweep_rand{rand_id}_' + fname if random_walk else f'sweep_nom{rand_id}_' + fname + fdir = Path(self.root_dir) if model_dir is None else Path(model_dir) / f'sweep_{rand_id}' + fig.savefig(fdir / f'{fname}.png', dpi=300, format='png') + save_dict = {'slice_idx': slice_idx, 'qoi_idx': qoi_idx, 'show_model': show_model, 'show_surr': show_surr, + 'nominal': nominal, 'random_walk': random_walk, 'xs': xs, 'ys_model': ys_model, + 'ys_surr': ys_surr} + with open(fdir / f'{fname}.pkl', 'wb') as fd: + pickle.dump(save_dict, fd) + + return fig, axs + + def plot_allocation(self, cmap: str = 'Blues', text_bar_width: float = 0.06, arrow_bar_width: float = 0.02): + """Plot bar charts showing cost allocation during training. + + !!! Warning "Beta feature" + This has pretty good default settings, but it might look terrible for your use. Mostly provided here as + a template for making cost allocation bar charts. Please feel free to copy and edit in your own code. + + :param cmap: the colormap string identifier for `plt` + :param text_bar_width: the minimum total cost fraction above which a bar will print centered model fidelity text + :param arrow_bar_width: the minimum total cost fraction above which a bar will try to print text with an arrow; + below this amount, the bar is too skinny and won't print any text + :returns: `fig, ax`, Figure and Axes objects + """ + # Get total cost (including offline overhead) + train_alloc, offline_alloc, cost_cum = self.get_allocation() + total_cost = cost_cum[-1] + for node, alpha_dict in offline_alloc.items(): + for alpha, cost in alpha_dict.items(): + total_cost += cost[1] + + # Remove nodes with cost=0 from alloc dicts (i.e. analytical models) + remove_nodes = [] + for node, alpha_dict in train_alloc.items(): + if len(alpha_dict) == 0: + remove_nodes.append(node) + for node in remove_nodes: + del train_alloc[node] + del offline_alloc[node] + + # Bar chart showing cost allocation breakdown for MF system at end + fig, axs = plt.subplots(1, 2, sharey='row') + width = 0.7 + x = np.arange(len(train_alloc)) + xlabels = list(train_alloc.keys()) + cmap = plt.get_cmap(cmap) + for k in range(2): + ax = axs[k] + alloc = train_alloc if k == 0 else offline_alloc + ax.set_title('Online training' if k == 0 else 'Overhead') + for j, (node, alpha_dict) in enumerate(alloc.items()): + bottom = 0 + c_intervals = np.linspace(0, 1, len(alpha_dict)) + bars = [(alpha, cost, cost[1] / total_cost) for alpha, cost in alpha_dict.items()] + bars = sorted(bars, key=lambda ele: ele[2], reverse=True) + for i, (alpha, cost, frac) in enumerate(bars): + p = ax.bar(x[j], frac, width, color=cmap(c_intervals[i]), linewidth=1, + edgecolor=[0, 0, 0], bottom=bottom) + bottom += frac + if frac > text_bar_width: + ax.bar_label(p, labels=[f'{alpha}, {round(cost[0])}'], label_type='center') + elif frac > arrow_bar_width: + xy = (x[j] + width / 2, bottom - frac / 2) # Label smaller bars with a text off to the side + ax.annotate(f'{alpha}, {round(cost[0])}', xy, xytext=(xy[0] + 0.2, xy[1]), + arrowprops={'arrowstyle': '->', 'linewidth': 1}) + else: + pass # Don't label really small bars + ax_default(ax, '', "Fraction of total cost" if k == 0 else '', legend=False) + ax.set_xticks(x, xlabels) + ax.set_xlim(left=-1, right=x[-1] + 1) + fig.set_size_inches(2.5*len(x), 4) + fig.tight_layout() + + if self.root_dir is not None: + fig.savefig(Path(self.root_dir) / 'mf_allocation.png', dpi=300, format='png') + + return fig, axs + + def get_component(self, comp_name: str) -> ComponentSurrogate: + """Return the `ComponentSurrogate` object for this component. + + :param comp_name: name of the component to return + :returns: the `ComponentSurrogate` object + """ + return self.graph.nodes[comp_name]['surrogate'] + + def _print_title_str(self, title_str: str): + """Log an important message.""" + self.logger.info('-' * int(len(title_str)/2) + title_str + '-' * int(len(title_str)/2)) + + def save_to_file(self, filename: str, save_dir: str | Path = None): + """Save the SystemSurrogate object to a .pkl file. + + :param filename: filename of the .pkl file to save to + :param save_dir: overrides existing surrogate root directory if provided + """ + if self.root_dir is None and save_dir is None: + # Can't save to file if root_dir is None + return + + save_dir = save_dir if save_dir is not None else str(Path(self.root_dir) / 'sys') + if not Path(save_dir).is_dir(): + save_dir = '.' + + exec_temp = self.executor # Temporarily save executor obj (can't pickle it) + self.set_executor(None) + with open(Path(save_dir) / filename, 'wb') as dill_file: + dill.dump(self, dill_file) + self.set_executor(exec_temp) + self.logger.info(f'SystemSurrogate saved to {(Path(save_dir) / filename).resolve()}') + + def _set_output_dir(self, set_dict: dict[str: str | Path]): + """Set the output directory for each component in `set_dict`. + + :param set_dict: a `dict` of component names (`str`) to their new output directories + """ + for node, node_obj in self.graph.nodes.items(): + if node in set_dict: + node_obj['surrogate']._set_output_dir(set_dict.get(node)) + + def set_root_directory(self, root_dir: str | Path, stdout: bool = True): + """Set the root to a new directory, for example if you move to a new filesystem. + + :param root_dir: new root directory + :param stdout: whether to connect the logger to console (default) + """ + self.root_dir = str(Path(root_dir).resolve()) + log_file = None + if not (Path(self.root_dir) / 'sys').is_dir(): + os.mkdir(Path(self.root_dir) / 'sys') + if not (Path(self.root_dir) / 'components').is_dir(): + os.mkdir(Path(self.root_dir) / 'components') + for f in os.listdir(self.root_dir): + if f.endswith('.log'): + log_file = str((Path(self.root_dir) / f).resolve()) + break + if log_file is None: + fname = datetime.datetime.now(tz=timezone.utc).isoformat().split('.')[0].replace(':', '.') + 'UTC_sys.log' + log_file = str((Path(self.root_dir) / fname).resolve()) + + # Setup the log file + self.log_file = log_file + self.logger = get_logger(self.__class__.__name__, log_file=log_file, stdout=stdout) + + # Update model output directories + for node, node_obj in self.graph.nodes.items(): + surr = node_obj['surrogate'] + surr.logger = get_logger(surr.__class__.__name__, log_file=log_file, stdout=stdout) + surr.log_file = self.log_file + if surr.save_enabled(): + output_dir = str((Path(self.root_dir) / 'components' / node).resolve()) + if not Path(output_dir).is_dir(): + os.mkdir(output_dir) + surr._set_output_dir(output_dir) + + def __getitem__(self, component: str) -> ComponentSurrogate: + """Convenience method to get the `ComponentSurrogate object` from the `SystemSurrogate`. + + :param component: the name of the component to get + :returns: the `ComponentSurrogate` object + """ + return self.get_component(component) + + def __repr__(self): + s = f'----SystemSurrogate----\nAdjacency: \n{nx.to_numpy_array(self.graph, dtype=int)}\n' \ + f'Exogenous inputs: {[str(var) for var in self.exo_vars]}\n' + for node, node_obj in self.graph.nodes.items(): + s += f'Component: {node}\n{node_obj["surrogate"]}' + return s + + def __str__(self): + return self.__repr__() + + def set_executor(self, executor: Executor | None): + """Set a new `concurrent.futures.Executor` object for parallel calls. + + :param executor: the new `Executor` object + """ + self.executor = executor + for node, node_obj in self.graph.nodes.items(): + node_obj['surrogate'].executor = executor + + @staticmethod + def load_from_file(filename: str | Path, root_dir: str | Path = None, executor: Executor = None): + """Load a `SystemSurrogate object` from file. + + :param filename: the .pkl file to load + :param root_dir: folder to use as the root directory, (uses file's second parent directory by default) + :param executor: a `concurrent.futures.Executor` object to set; clears it if None + :returns: the `SystemSurrogate` object + """ + if root_dir is None: + root_dir = Path(filename).parent.parent # Assume root/sys/filename.pkl + + with open(Path(filename), 'rb') as dill_file: + sys_surr = dill.load(dill_file) + sys_surr.set_executor(executor) + sys_surr.set_root_directory(root_dir) + sys_surr.logger.info(f'SystemSurrogate loaded from {Path(filename).resolve()}') + + return sys_surr + + @staticmethod + def _constrained_lls(A: np.ndarray, b: np.ndarray, C: np.ndarray, d: np.ndarray) -> np.ndarray: + """Minimize $||Ax-b||_2$, subject to $Cx=d$, i.e. constrained linear least squares. + + !!! Note + See http://www.seas.ucla.edu/~vandenbe/133A/lectures/cls.pdf for more detail. + + :param A: `(..., M, N)`, vandermonde matrix + :param b: `(..., M, 1)`, data + :param C: `(..., P, N)`, constraint operator + :param d: `(..., P, 1)`, constraint condition + :returns: `(..., N, 1)`, the solution parameter vector `x` + """ + M = A.shape[-2] + dims = len(A.shape[:-2]) + T_axes = tuple(np.arange(0, dims)) + (-1, -2) + Q, R = np.linalg.qr(np.concatenate((A, C), axis=-2)) + Q1 = Q[..., :M, :] + Q2 = Q[..., M:, :] + Q1_T = np.transpose(Q1, axes=T_axes) + Q2_T = np.transpose(Q2, axes=T_axes) + Qtilde, Rtilde = np.linalg.qr(Q2_T) + Qtilde_T = np.transpose(Qtilde, axes=T_axes) + Rtilde_T_inv = np.linalg.pinv(np.transpose(Rtilde, axes=T_axes)) + w = np.linalg.pinv(Rtilde) @ (Qtilde_T @ Q1_T @ b - Rtilde_T_inv @ d) + + return np.linalg.pinv(R) @ (Q1_T @ b - Q2_T @ w) diff --git a/src/amisc/utils.py b/src/amisc/utils.py new file mode 100644 index 0000000..6568925 --- /dev/null +++ b/src/amisc/utils.py @@ -0,0 +1,324 @@ +"""`utils.py` + +Provides some basic utilities for the package. + +Includes +-------- +- `load_variables`: convenience function for loading RVs from a .json config file +- `get_logger`: logging utility with nice formatting +- `ax_default`: plotting utility with nice formatting +- `approx_hess`: finite difference approximation of the Hessian +- `batch_normal_sample`: helper function to sample from arbitrarily-sized Gaussian distribution(s) +- `ndscatter`: plotting utility for n-dimensional data +""" +import json +from pathlib import Path +import logging +import sys + +import matplotlib.pyplot as plt +import matplotlib +from matplotlib.pyplot import cycler +from matplotlib.colors import LinearSegmentedColormap, ListedColormap +from matplotlib.ticker import StrMethodFormatter, AutoLocator, FuncFormatter +import scipy.stats as st +import numpy as np +from numpy.linalg.linalg import LinAlgError + +from amisc.rv import BaseRV, UniformRV, NormalRV, ScalarRV + + +LOG_FORMATTER = logging.Formatter("%(asctime)s \u2014 [%(levelname)s] \u2014 %(name)-20s \u2014 %(message)s") + + +def load_variables(variables: list[str], file: Path | str) -> list[BaseRV]: + """Load a list of BaseRV objects from a variables json `file`. + + :param variables: a list of str ids for variables to find in `file` + :param file: json file to search for variable definitions + :returns rvs: a list of corresponding `BaseRV` objects + """ + with open(Path(file), 'r') as fd: + data = json.load(fd) + + rvs = [] + keys = ['id', 'tex', 'description', 'units', 'param_type', 'nominal', 'domain'] + for str_id in variables: + if str_id in data: + var_info = data.get(str_id) + kwargs = {key: var_info.get(key) for key in keys if var_info.get(key)} + match var_info.get('rv_type', 'none'): + case 'uniform_bds': + bds = var_info.get('rv_params') + rvs.append(UniformRV(bds[0], bds[1], **kwargs)) + case 'uniform_pct': + rvs.append(UniformRV(var_info.get('rv_params'), 'pct', **kwargs)) + case 'uniform_tol': + rvs.append(UniformRV(var_info.get('rv_params'), 'tol', **kwargs)) + case 'normal': + mu, std = var_info.get('rv_params') + rvs.append(NormalRV(mu, std, **kwargs)) + case 'none': + # Make a plain stand-in scalar RV object (no uncertainty) + rvs.append(ScalarRV(**kwargs)) + case other: + raise NotImplementedError(f'RV type "{var_info.get("rv_type")}" is not known.') + else: + raise ValueError(f'You have requested the variable {str_id}, but it was not found in {file}. ' + f'Please add a definition of {str_id} to {file} or construct it on your own.') + + return rvs + + +def get_logger(name: str, stdout=True, log_file: str | Path = None) -> logging.Logger: + """Return a file/stdout logger with the given name. + + :param name: the name of the logger to return + :param stdout: whether to add a stdout handler to the logger + :param log_file: add file logging to this file (optional) + :returns: the logger + """ + logger = logging.getLogger(name) + logger.setLevel(logging.DEBUG) + logger.handlers.clear() + if stdout: + std_handler = logging.StreamHandler(sys.stdout) + std_handler.setFormatter(LOG_FORMATTER) + logger.addHandler(std_handler) + if log_file is not None: + f_handler = logging.FileHandler(log_file, mode='a', encoding='utf-8') + f_handler.setLevel(logging.DEBUG) + f_handler.setFormatter(LOG_FORMATTER) + logger.addHandler(f_handler) + + return logger + + +def approx_hess(func: callable, theta: np.ndarray, pert=0.01) -> np.ndarray: + """Approximate Hessian of `func` at a specified `theta` location. + + :param func: expects to be called as `func(theta) -> (..., y_dim)`, where `y_dim=1` (scalar funcs only) + :param theta: `(..., theta_dim)`, points to linearize model about + :param pert: perturbation for approximate partial derivatives + :returns H: `(..., theta_dim, theta_dim)`, the approximate Hessian `(theta_dim, theta_dim)` at all locations (...) + """ + theta = np.atleast_1d(theta) + shape = theta.shape[:-1] # (*) + theta_dim = theta.shape[-1] # Number of parameters + dtheta = pert * theta + + # Return the Hessians (..., theta_dim, theta_dim) + H = np.empty(shape + (theta_dim, theta_dim)) + + for i in range(theta_dim): + for j in range(i, theta_dim): + # Allocate space at 4 grid points (n1=-1, p1=+1) + theta_n1_n1 = np.copy(theta) + theta_p1_p1 = np.copy(theta) + theta_n1_p1 = np.copy(theta) + theta_p1_n1 = np.copy(theta) + + # Perturbations to theta in each direction + theta_n1_n1[..., i] -= dtheta[..., i] + theta_n1_n1[..., j] -= dtheta[..., j] + f_n1_n1 = func(theta_n1_n1) + + theta_p1_p1[..., i] += dtheta[..., i] + theta_p1_p1[..., j] += dtheta[..., j] + f_p1_p1 = func(theta_p1_p1) + + theta_n1_p1[..., i] -= dtheta[..., i] + theta_n1_p1[..., j] += dtheta[..., j] + f_n1_p1 = func(theta_n1_p1) + + theta_p1_n1[..., i] += dtheta[..., i] + theta_p1_n1[..., j] -= dtheta[..., j] + f_p1_n1 = func(theta_p1_n1) + + res = (f_n1_n1 + f_p1_p1 - f_n1_p1 - f_p1_n1) / np.expand_dims(4 * dtheta[..., i] * dtheta[..., j], + axis=-1) + + # Hessian only computed for scalar functions, y_dim=1 on last axis + H[..., i, j] = np.squeeze(res, axis=-1) + H[..., j, i] = np.squeeze(res, axis=-1) + + return H + + +def batch_normal_sample(mean: np.ndarray | float, cov: np.ndarray | float, size: tuple | int = ()) -> np.ndarray: + """Batch sample multivariate normal distributions (pretty much however you want). + + :param mean: `(..., dim)`, expected values, where dim is the random variable dimension + :param cov: `(..., dim, dim)`, covariance matrices + :param size: shape of additional samples + :returns samples: `(*size, ..., dim)`, samples from multivariate distributions + """ + mean = np.atleast_1d(mean) + cov = np.atleast_2d(cov) + + if isinstance(size, int): + size = (size, ) + shape = size + np.broadcast_shapes(mean.shape, cov.shape[:-1]) + x_normal = np.random.standard_normal((*shape, 1)).astype(np.float32) + samples = np.squeeze(np.linalg.cholesky(cov) @ x_normal, axis=-1) + mean + return samples + + +def ax_default(ax: plt.Axes, xlabel='', ylabel='', legend=True, cmap='tab10'): + """Nice default formatting for plotting X-Y data. + + :param ax: the axes to apply these settings to + :param xlabel: the xlabel to set for `ax` + :param ylabel: the ylabel to set for `ax` + :param legend: whether to show a legend + :param cmap: colormap to use for cycling + """ + plt.rcParams["axes.prop_cycle"] = get_cycle(cmap) + plt.rc('xtick', labelsize='small') + plt.rc('ytick', labelsize='small') + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + ax.tick_params(axis='x', direction='in') + ax.tick_params(axis='y', direction='in') + if legend: + leg = ax.legend(fancybox=True) + frame = leg.get_frame() + frame.set_edgecolor('k') + + +def get_cycle(cmap: str | matplotlib.colors.Colormap, num_colors: int = None): + """Get a color cycler for plotting. + + :param cmap: a string specifier of a matplotlib colormap (or a colormap instance) + :param num_colors: the number of colors to cycle through + """ + use_index = False + if isinstance(cmap, str): + use_index = cmap in ['Pastel1', 'Pastel2', 'Paired', 'Accent', 'Dark2', 'Set1', 'Set2', 'Set3', + 'tab10', 'tab20', 'tab20b', 'tab20c'] + cmap = matplotlib.cm.get_cmap(cmap) + if num_colors is None: + num_colors = cmap.N + if cmap.N > 100: + use_index = False + elif isinstance(cmap, LinearSegmentedColormap): + use_index = False + elif isinstance(cmap, ListedColormap): + use_index = True + if use_index: + ind = np.arange(int(num_colors)) % cmap.N + return cycler("color", cmap(ind)) + else: + colors = cmap(np.linspace(0, 1, num_colors)) + return cycler("color", colors) + + +def ndscatter(samples: np.ndarray, labels: list[str] = None, tick_fmts: list[str] = None, plot='scatter', + cmap='viridis', bins=20, z: np.ndarray = None, cb_label=None, cb_norm='linear', subplot_size=3): + """Triangle scatter plots of n-dimensional samples. + + !!! Warning + Best for `dim < 10`. You can shrink the `subplot_size` to assist graphics loading time. + + :param samples: `(N, dim)` samples to plot + :param labels: list of axis labels of length `dim` + :param tick_fmts: list of str.format() specifiers for ticks, e.g `['{x: ^10.2f}', ...]`, of length `dim` + :param plot: 'hist' for 2d hist plot, 'kde' for kernel density estimation, or 'scatter' (default) + :param cmap: the matplotlib string specifier of a colormap + :param bins: number of bins in each dimension for histogram marginals + :param z: `(N,)` a performance metric corresponding to `samples`, used to color code the scatter plot if provided + :param cb_label: label for color bar (if `z` is provided) + :param cb_norm: `str` or `plt.colors.Normalize`, normalization method for plotting `z` on scatter plot + :param subplot_size: size in inches of a single 2d marginal subplot + :returns fig, axs: the `plt` Figure and Axes objects, (returns an additional `cb_fig, cb_ax` if `z` is specified) + """ + N, dim = samples.shape + x_min = np.min(samples, axis=0) + x_max = np.max(samples, axis=0) + if labels is None: + labels = [f"x{i}" for i in range(dim)] + if z is None: + z = np.zeros(N) + if cb_label is None: + cb_label = 'Performance metric' + + def tick_format_func(value, pos): + if value > 1: + return f'{value:.2f}' + if value > 0.01: + return f'{value:.4f}' + if value < 0.01: + return f'{value:.2E}' + default_ticks = FuncFormatter(tick_format_func) + # if tick_fmts is None: + # tick_fmts = ['{x:.2G}' for i in range(dim)] + + # Set up triangle plot formatting + fig, axs = plt.subplots(dim, dim, sharex='col', sharey='row') + for i in range(dim): + for j in range(dim): + ax = axs[i, j] + if i == j: # 1d marginals on diagonal + # ax.get_shared_y_axes().remove(ax) + ax._shared_axes['y'].remove(ax) + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + ax.spines['left'].set_visible(False) + if i == 0: + ax.get_yaxis().set_ticks([]) + if j > i: # Clear the upper triangle + ax.axis('off') + if i == dim - 1: # Bottom row + ax.set_xlabel(labels[j]) + ax.xaxis.set_major_locator(AutoLocator()) + formatter = StrMethodFormatter(tick_fmts[j]) if tick_fmts is not None else default_ticks + ax.xaxis.set_major_formatter(formatter) + if j == 0 and i > 0: # Left column + ax.set_ylabel(labels[i]) + ax.yaxis.set_major_locator(AutoLocator()) + formatter = StrMethodFormatter(tick_fmts[i]) if tick_fmts is not None else default_ticks + ax.yaxis.set_major_formatter(formatter) + + # Plot marginals + for i in range(dim): + for j in range(dim): + ax = axs[i, j] + if i == j: # 1d marginals (on diagonal) + c = plt.get_cmap(cmap)(0) + if plot == 'kde': + kernel = st.gaussian_kde(samples[:, i]) + x = np.linspace(x_min[i], x_max[i], 1000) + ax.fill_between(x, y1=kernel(x), y2=0, lw=0, alpha=0.3, facecolor=c) + ax.plot(x, kernel(x), ls='-', c=c, lw=1.5) + else: + ax.hist(samples[:, i], edgecolor='black', color=c, density=True, alpha=0.5, + linewidth=1.2, bins='auto') + if j < i: # 2d marginals (lower triangle) + ax.set_xlim([x_min[j], x_max[j]]) + ax.set_ylim([x_min[i], x_max[i]]) + if plot == 'scatter': + sc = ax.scatter(samples[:, j], samples[:, i], s=1.5, c=z, cmap=cmap, norm=cb_norm) + elif plot == 'hist': + ax.hist2d(samples[:, j], samples[:, i], bins=bins, density=True, cmap=cmap) + elif plot == 'kde': + kernel = st.gaussian_kde(samples[:, [j, i]].T) + xg, yg = np.meshgrid(np.linspace(x_min[j], x_max[j], 60), np.linspace(x_min[i], x_max[i], 60)) + x = np.vstack([xg.ravel(), yg.ravel()]) + zg = np.reshape(kernel(x), xg.shape) + ax.contourf(xg, yg, zg, 5, cmap=cmap, alpha=0.9) + ax.contour(xg, yg, zg, 5, colors='k', linewidths=1.5) + else: + raise NotImplementedError('This plot type is not known. plot=["hist", "kde", "scatter"]') + + fig.set_size_inches(subplot_size * dim, subplot_size * dim) + fig.tight_layout() + + # Plot colorbar in standalone figure + if np.max(z) > 0 and plot == 'scatter': + cb_fig, cb_ax = plt.subplots(figsize=(1.5, 6)) + cb_fig.subplots_adjust(right=0.7) + cb_fig.colorbar(sc, cax=cb_ax, orientation='vertical', label=cb_label) + cb_fig.tight_layout() + return fig, axs, cb_fig, cb_ax + + return fig, axs diff --git a/tests/test_component.py b/tests/test_component.py new file mode 100644 index 0000000..597fcc6 --- /dev/null +++ b/tests/test_component.py @@ -0,0 +1,60 @@ +import numpy as np +import matplotlib.pyplot as plt + +from amisc.component import SparseGridSurrogate +from amisc.rv import UniformRV +from amisc.utils import ax_default + + +def test_sparse_grid(plots=False): + """Simple cos test from Jakeman (2022)""" + def model(x, alpha): + alpha = np.atleast_1d(alpha) # (1,) + eps = (1/5) * 2.0**(-alpha[0]) + y = np.cos(np.pi/2 * (x + 4/5 + eps)) + return {'y': y} + + # Construct MISC surrogate from an index set + Ik = [((0,), (0,)), ((0,), (1,)), ((1,), (0,)), ((2,), (0,)), ((1,), (1,)), ((0,), (2,)), ((1,), (2,)), + ((2,), (1,)), ((2,), (2,))] + x_vars = UniformRV(-1, 1) + truth_alpha = (15,) + comp = SparseGridSurrogate(x_vars, model, multi_index=Ik, truth_alpha=truth_alpha) + N = 100 + xg = np.linspace(-1, 1, N).reshape((N, 1)) + yt = comp(xg, use_model=truth_alpha) + y_surr = comp(xg) + l2_error = np.sqrt(np.mean((y_surr - yt)**2)) / np.sqrt(np.mean(yt**2)) + assert l2_error < 0.1 + + # Plot results for each fidelity of the MISC surrogate + if plots: + fig, axs = plt.subplots(3, 3, sharey='row', sharex='col') + for alpha in range(3): + for beta in range(3): + ax = axs[2-alpha, beta] + surr = comp.get_sub_surrogate((alpha,), (beta,), include_grid=True) + s = f'$\hat{{f}}_{{{alpha}, {beta}}}$' + ax.plot(xg, surr(xg), '--k', label=r'{}'.format(s), linewidth=1.5) + s = f'$\hat{{f}}_{alpha}$' + ax.plot(xg, model(xg, alpha)['y'], '--b', label=r'{}'.format(s), linewidth=2) + ax.plot(xg, yt, '-r', label=r'$f$', linewidth=2) + ax.plot(surr.xi, surr.yi, 'or') + xlabel = r'$x$' if alpha == 0 else '' + ylabel = r'$f(x)$' if beta == 0 else '' + ax_default(ax, xlabel, ylabel, legend=True) + + fig.text(0.5, 0.02, r'Increasing surrogate fidelity ($\beta$) $\rightarrow$', ha='center', fontweight='bold') + fig.text(0.02, 0.5, r'Increasing model fidelity ($\alpha$) $\rightarrow$', va='center', fontweight='bold', rotation='vertical') + fig.set_size_inches(3 * 3, 3 * 3) + fig.tight_layout(pad=3, w_pad=1, h_pad=1) + plt.show() + + fig, ax = plt.subplots() + ax.plot(xg, yt, '-r', linewidth=2, label='Model') + ax.plot(xg, y_surr, '--k', linewidth=1.5, label='MISC surrogate') + ax_default(ax, r'$x$', r'$f(x)$', legend=True) + plt.show() + + +# TODO: add tests for testing parallel execution, writing output files, and refining surrogate diff --git a/tests/test_interpolator.py b/tests/test_interpolator.py new file mode 100644 index 0000000..6e13641 --- /dev/null +++ b/tests/test_interpolator.py @@ -0,0 +1,113 @@ +import numpy as np +import matplotlib.pyplot as plt + +from amisc.interpolator import LagrangeInterpolator +from amisc.rv import UniformRV +from amisc.utils import ax_default +from amisc.examples.models import tanh_func, nonlinear_wave + + +def test_tensor_product_1d(plots=False): + beta = (3,) + x_var = UniformRV(0, 1) + x_grid = np.linspace(0, 1, 100).reshape((100, 1)) + y_grid = tanh_func(x_grid)['y'] + interp = LagrangeInterpolator(beta, x_var, model=tanh_func) + interp.set_yi() + y_interp = interp(x_grid) + + # Refine + beta2 = (4,) + interp2 = interp.refine(beta2) + y2_interp = interp2(x_grid) + + # Compute errors + N = 1000 + xtest = np.random.rand(N, 1) + ytest = interp2(xtest) + ytruth = tanh_func(xtest)['y'] + l2_error = np.sqrt(np.mean((ytest - ytruth) ** 2)) / np.sqrt(np.mean(ytruth ** 2)) + assert l2_error < 1e-1 + + # Plot results + if plots: + fig, ax = plt.subplots() + ax.plot(x_grid, y_grid, '-k', label='Model') + ax.plot(interp.xi, interp.yi, 'or', markersize=6, label=r'Training data') + ax.plot(x_grid, y_interp, '-r', label=r'$\beta=3$') + ax.plot(interp2.xi, interp2.yi, 'ob', markersize=4) + ax.plot(x_grid, y2_interp, '-b', label=r'$\beta=4$') + ax_default(ax, r'Input', r'Output', legend=True) + fig.tight_layout() + plt.show() + + +def test_tensor_product_2d(plots=False): + bb_2d_func = lambda x: nonlinear_wave(x, env_var=0.2**2, wave_amp=0.3) + beta = (3, 3) + x_vars = [UniformRV(0, 1), UniformRV(0, 1)] + N = 50 + x_grid = np.linspace(0, 1, N) + xg, yg = np.meshgrid(x_grid, x_grid) + xg = xg.reshape((N, N, 1)) + yg = yg.reshape((N, N, 1)) + x = np.concatenate((xg, yg), axis=-1) + z = bb_2d_func(x)['y'] + + # Set up interpolator + interp = LagrangeInterpolator(beta, x_vars, model=bb_2d_func) + interp.set_yi() + z_interp = interp(x) + error = np.abs(z_interp - z) + + # Refine interpolator + beta2 = (3, 4) + interp2 = interp.refine(beta2) + z2_interp = interp2(x) + error2 = np.abs(z2_interp - z) + vmin = min(np.min(z_interp), np.min(z), np.min(z2_interp)) + vmax = max(np.max(z_interp), np.max(z), np.max(z2_interp)) + emin = min(np.min(error), np.min(error2)) + emax = max(np.max(error), np.max(error2)) + + # Compute errors + N = 1000 + xtest = np.random.rand(N, 2) + ytest = interp2(xtest) + ytruth = bb_2d_func(xtest)['y'] + l2_error = np.sqrt(np.mean((ytest - ytruth)**2)) / np.sqrt(np.mean(ytruth**2)) + assert l2_error < 1e-1 + + if plots: + fig, ax = plt.subplots(2, 3) + c1 = ax[0, 0].contourf(xg.squeeze(), yg.squeeze(), z.squeeze(), 60, cmap='coolwarm', vmin=vmin, vmax=vmax) + plt.colorbar(c1, ax=ax[0, 0]) + ax[0, 0].set_title('True function') + ax_default(ax[0, 0], r'$x_1$', r'$x_2$', legend=False) + c2 = ax[0, 1].contourf(xg.squeeze(), yg.squeeze(), z_interp.squeeze(), 60, cmap='coolwarm', vmin=vmin, vmax=vmax) + ax[0, 1].plot(interp.xi[:, 0], interp.xi[:, 1], 'o', markersize=6, markerfacecolor='green') + plt.colorbar(c2, ax=ax[0, 1]) + ax[0, 1].set_title('Interpolant') + ax_default(ax[0, 1], r'$x_1$', '', legend=False) + c3 = ax[0, 2].contourf(xg.squeeze(), yg.squeeze(), error.squeeze(), 60, cmap='viridis', vmin=emin, vmax=emax) + ax[0, 2].plot(interp.xi[:, 0], interp.xi[:, 1], 'o', markersize=6, markerfacecolor='green') + plt.colorbar(c3, ax=ax[0, 2]) + ax[0, 2].set_title('Absolute error') + ax_default(ax[0, 2], r'$x_1$', '', legend=False) + c1 = ax[1, 0].contourf(xg.squeeze(), yg.squeeze(), z.squeeze(), 60, cmap='coolwarm', vmin=vmin, vmax=vmax) + plt.colorbar(c1, ax=ax[1, 0]) + ax[1, 0].set_title('True function') + ax_default(ax[1, 0], r'$x_1$', r'$x_2$', legend=False) + c2 = ax[1, 1].contourf(xg.squeeze(), yg.squeeze(), z2_interp.squeeze(), 60, cmap='coolwarm', vmin=vmin, vmax=vmax) + ax[1, 1].plot(interp2.xi[:, 0], interp2.xi[:, 1], 'o', markersize=6, markerfacecolor='green') + plt.colorbar(c2, ax=ax[1, 1]) + ax[1, 1].set_title('Refined') + ax_default(ax[1, 1], r'$x_1$', '', legend=False) + c3 = ax[1, 2].contourf(xg.squeeze(), yg.squeeze(), error2.squeeze(), 60, cmap='viridis', vmin=emin, vmax=emax) + ax[1, 2].plot(interp2.xi[:, 0], interp2.xi[:, 1], 'o', markersize=6, markerfacecolor='green') + plt.colorbar(c3, ax=ax[1, 2]) + ax[1, 2].set_title('Absolute error') + ax_default(ax[1, 2], r'$x_1$', '', legend=False) + fig.set_size_inches(15, 11) + fig.tight_layout() + plt.show() diff --git a/tests/test_rv.py b/tests/test_rv.py new file mode 100644 index 0000000..678cace --- /dev/null +++ b/tests/test_rv.py @@ -0,0 +1,20 @@ +from amisc.rv import UniformRV, NormalRV, LogNormalRV, LogUniformRV, ScalarRV + + +# TODO: make sure rvs are sampling correctly + +def test_loading_rvs(): + """Test random variable construction and methods.""" + rvs = [UniformRV(0, 1, 'x'), NormalRV(0, 1, 'y'), LogUniformRV(1, 2, 'z'), LogNormalRV(2, 1, 'a'), + ScalarRV('h')] + samples_pdf = [] + samples_domain = [] + pdf_vals = [] + labels = [] + for v in rvs: + sample = v.sample(5) + samples_pdf.append(sample) + samples_domain.append(v.sample_domain(5)) + pdf_vals.append(v.pdf(sample)) + labels.append(v.to_tex(symbol=False)) + v.update_bounds(*v.bds) diff --git a/tests/test_system.py b/tests/test_system.py new file mode 100644 index 0000000..c53a54a --- /dev/null +++ b/tests/test_system.py @@ -0,0 +1,342 @@ +import time +import warnings +from pathlib import Path +import sys + +import numpy as np +from scipy.linalg import lapack +from scipy.optimize import fsolve +from scipy.stats import gaussian_kde +import matplotlib.pyplot as plt +import pytest + +from amisc.system import SystemSurrogate, ComponentSpec +from amisc.rv import UniformRV +from amisc.utils import ax_default +from amisc.examples.models import fire_sat_system + + +# TODO: Include a swap and insert component test + +@pytest.mark.skipif(not sys.platform.startswith('linux'), reason='not sure why') +def test_fire_sat(plots=True): + """Test the fire satellite coupled system from Chaudhuri (2018)""" + N = 100 + surr = fire_sat_system(save_dir=Path('.')) + xt = surr.sample_inputs(N, use_pdf=True) + yt = surr(xt, use_model='best') + use_idx = ~np.any(np.isnan(yt), axis=-1) + xt = xt[use_idx, :] + yt = yt[use_idx, :] + test_set = {'xt': xt, 'yt': yt} + n_jobs = -1 if sys.platform.startswith('linux') else 1 # parallel on VM runners on Github get stuck for non-linux + surr.fit(max_iter=15, max_tol=1e-2, max_runtime=1/12, test_set=test_set, n_jobs=n_jobs, + num_refine=1000) + + ysurr = surr.predict(xt) + l2_error = np.sqrt(np.nanmean((ysurr-yt)**2, axis=0)) / np.sqrt(np.nanmean(yt**2, axis=0)) + assert np.nanmax(l2_error) < 0.1 + + if plots: + # Plot error histograms + fig, ax = plt.subplots(1, 3) + ax[0].hist(yt[:, 0], color='red', bins=20, edgecolor='black', density=True, linewidth=1.2, label='Truth') + ax[0].hist(ysurr[:, 0], color='blue', bins=20, edgecolor='black', density=True, linewidth=1.2, alpha=0.4, label='Surrogate') + ax[1].hist(yt[:, 7], color='red', bins=20, edgecolor='black', density=True, linewidth=1.2, label='Truth') + ax[1].hist(ysurr[:, 7], color='blue', bins=20, edgecolor='black', density=True, linewidth=1.2, alpha=0.4, label='Surrogate') + ax[2].hist(yt[:, 8], color='red', bins=20, edgecolor='black', density=True, linewidth=1.2, label='Truth') + ax[2].hist(ysurr[:, 8], color='blue', bins=20, edgecolor='black', density=True, linewidth=1.2, alpha=0.4, label='Surrogate') + ax_default(ax[0], 'Satellite velocity ($m/s$)', '', legend=True) + ax_default(ax[1], 'Solar panel area ($m^2$)', '', legend=True) + ax_default(ax[2], 'Attitude control power ($W$)', '', legend=True) + fig.set_size_inches(9, 3) + fig.tight_layout() + + # Plot 1d slices + slice_idx = ['H', 'Po', 'Cd'] + qoi_idx = ['Vsat', 'Asa', 'Pat'] + try: + fig, ax = surr.plot_slice(slice_idx, qoi_idx, show_model=['best', 'worst'], model_dir=surr.root_dir, + random_walk=True, N=10) + except np.linalg.LinAlgError: + print(f'Its alright. Sometimes the random walks are wacky and FPI wont converge.') + + # Plot allocation bar charts + fig, ax = surr.plot_allocation() + plt.show() + + +def test_system_refine(plots=False): + """Test iterative refinement for Figure 5 in Jakeman 2022""" + def coupled_system(): + def f1(x): + return {'y': x * np.sin(np.pi * x)} + def f2(x): + return {'y': 1 / (1 + 25 * x ** 2)} + return f1, f2 + + f1, f2 = coupled_system() + exo_vars = [UniformRV(0, 1)] + coupling_vars = [UniformRV(0, 1), UniformRV(0, 1)] + comp1 = ComponentSpec(f1, name='Model1', exo_in=0, coupling_out=0) + comp2 = ComponentSpec(f2, name='Model2', coupling_in=0, coupling_out=1) + surr = SystemSurrogate([comp1, comp2], exo_vars, coupling_vars, init_surr=True) + + Niter = 4 + x = np.linspace(0, 1, 100).reshape((100, 1)) + y1 = f1(x)['y'] + y2 = f2(x)['y'] + y3 = f2(y1)['y'] + fig, ax = plt.subplots(Niter, 3, sharex='col', sharey='row') + for i in range(Niter): + # Plot actual function values + ax[i, 0].plot(x, y1, '-r', label='$f_1(x)$') + ax[i, 1].plot(x, y2, '-r', label='$f_2(y)$') + ax[i, 2].plot(x, y3, '-r', label='$f(x)$') + + # Plot first component surrogates + comp = surr.get_component('Model1') + ax[i, 0].plot(x, comp(x, training=True), '--k', label='$f_1$ current') + beta_max = 0 + for alpha, beta in comp.index_set: + if beta[0] > beta_max: + beta_max = beta[0] + interp = comp.get_sub_surrogate((), (beta_max,), include_grid=True) + ax[i, 0].plot(interp.xi, interp.yi, 'ok', markersize=8, label='') + for alpha, beta in comp.iterate_candidates(): + comp.update_misc_coeffs() + yJ1 = surr(x, training=True) + ax[i, 0].plot(x, comp(x, training=True), ':b', label='$f_1$ candidate') + interp = comp.get_sub_surrogate(alpha, beta, include_grid=True) + ax[i, 0].plot(interp.xi, interp.yi, 'xb', markersize=8, label='') + comp.update_misc_coeffs() + + # Plot second component surrogates + comp = surr.get_component('Model2') + ax[i, 1].plot(x, comp(x, training=True), '--k', label='$f_2$ current') + beta_max = 0 + for alpha, beta in comp.index_set: + if beta[0] > beta_max: + beta_max = beta[0] + interp = comp.get_sub_surrogate((), (beta_max,), include_grid=True) + ax[i, 1].plot(interp.xi, interp.yi, 'ok', markersize=8, label='') + for alpha, beta in comp.iterate_candidates(): + comp.update_misc_coeffs() + yJ2 = surr(x, training=True) + ax[i, 1].plot(x, comp(x, training=True), '-.g', label='$f_2$ candidate') + interp = comp.get_sub_surrogate(alpha, beta, include_grid=True) + ax[i, 1].plot(interp.xi, interp.yi, 'xg', markersize=8, label='') + comp.update_misc_coeffs() + + # Plot integrated surrogates + ysurr = surr(x, training=True) + ax[i, 2].plot(x, ysurr[:, 1:2], '--k', label='$f_J$') + ax[i, 2].plot(x, yJ1[:, 1:2], ':b', label='$f_{J_1}$') + ax[i, 2].plot(x, yJ2[:, 1:2], '-.g', label='$f_{J_2}$') + ax_default(ax[i, 0], '$x$', '$f_1(x)$', legend=True) + ax_default(ax[i, 1], '$y$', '$f_2(y)$', legend=True) + ax_default(ax[i, 2], '$x$', '$f_2(f_1(x))$', legend=True) + + # Refine the system + surr.refine(qoi_ind=None, num_refine=100) + + ysurr = surr.predict(x) + ytrue = surr.predict(x, use_model='best') + l2_error = np.sqrt(np.mean((ysurr-ytrue)**2, axis=0)) / np.sqrt(np.mean(ytrue**2, axis=0)) + assert np.max(l2_error) < 0.1 + + if plots: + fig.set_size_inches(3.5*3, 3.5*Niter) + fig.tight_layout() + plt.show() + + +def test_feedforward(plots=False): + """Test MD system in Figure 5 in Jakeman 2022""" + def coupled_system(): + def f1(x): + return {'y': x * np.sin(np.pi * x)} + def f2(x): + return {'y': 1 / (1 + 25*x**2)} + def f(x): + return f2(f1(x)['y'])['y'] + return f1, f2, f + + f1, f2, f = coupled_system() + x_in = UniformRV(0, 1) + y1, y2 = UniformRV(0, 1), UniformRV(0, 1) + comp1 = ComponentSpec(f1, exo_in=x_in, coupling_out=y1) + comp2 = ComponentSpec(f2, coupling_in=y1, coupling_out=y2) + surr = SystemSurrogate([comp1, comp2], x_in, [y1, y2], init_surr=False) + + x = np.linspace(0, 1, 100).reshape((100, 1)) + y1 = f1(x)['y'] + y2 = f(x) + y_surr = surr(x, use_model='best') + l2_error = np.sqrt(np.mean((y_surr[:, 1:2] - y2)**2)) / np.sqrt(np.mean(y2**2)) + assert l2_error < 1e-15 + + if plots: + fig, ax = plt.subplots(1, 2) + ax[0].plot(x, y1, '-r', label='$f_1(x)$') + ax[0].plot(np.squeeze(x), y_surr[:, 0], '--k', label='Surrogate') + ax_default(ax[0], '$x$', '$f(x)$', legend=True) + ax[1].plot(x, y2, '-r', label='$f(x)$') + ax[1].plot(np.squeeze(x), y_surr[:, 1], '--k', label='Surrogate') + ax_default(ax[1], '$x$', '$f(x)$', legend=True) + fig.set_size_inches(6, 3) + fig.tight_layout() + plt.show() + + +def test_system_surrogate(plots=False): + """Test the MD system in figure 6 in Jakeman 2022""" + def coupled_system(D1, D2, Q1, Q2): + def f1(x, /, alpha=(0,)): + eps = 10 ** (-float(alpha[0])) + q = np.arange(1, Q1+1).reshape((1,)*len(x.shape[:-1]) + (Q1,)) + return {'y': (x[..., 0, np.newaxis] ** (q-1)) * np.sin(np.sum(x, axis=-1, keepdims=True) + eps)} + + def f2(x, /, alpha=(0,)): + eps = 10 ** (-float(alpha[0])) + q = np.arange(1, Q2+1).reshape((1,)*len(x.shape) + (Q2,)) + prod1 = np.prod(x[..., D2:, np.newaxis] ** (q) - eps, axis=-2) # (..., Q2) + prod2 = np.prod(x[..., :D2], axis=-1, keepdims=True) # (..., 1) + return {'y': prod1 * prod2} + + def f3(x, /, alpha=(0,), D3=D1): + eps = 10 ** (-float(alpha[0])) + prod1 = np.exp(-np.sum((x[..., D3:] - eps) ** 2, axis=-1)) # (...,) + prod2 = 1 + (25/16)*np.sum(x[..., :D3], axis=-1) ** 2 # (...,) + return {'y': np.expand_dims(prod1 / prod2, axis=-1)} # (..., 1) + + def f(x): + # Ground truth (arbitrary high alpha) + alpha = (15,) + x1 = x[..., :D1] + y1 = f1(x1, alpha)['y'] + x2 = np.concatenate((x[..., D1:], y1), axis=-1) + y2 = f2(x2, alpha)['y'] + x3 = np.concatenate((x1, y2), axis=-1) + y3 = f3(x3, alpha)['y'] + return np.concatenate((y1, y2, y3), axis=-1) + + return f1, f2, f3, f + + # Hook up the 'wiring' for this example feedforward system + D1 = 1 + D2 = D1 + Q1 = 1 + Q2 = Q1 + alpha = (15,) + f1, f2, f3, f = coupled_system(D1, D2, Q1, Q2) + comp1 = ComponentSpec(f1, name='Cathode', truth_alpha=alpha, exo_in=list(np.arange(0, D1)), + coupling_out=list(np.arange(0, Q1)), max_alpha=5, max_beta=(3,)*D1) + comp2 = ComponentSpec(f2, name='Thruster', truth_alpha=alpha, exo_in=list(np.arange(D1, D1+D2)), max_alpha=5, + max_beta=(3,)*(D2+Q1), coupling_in={'Cathode': list(np.arange(0, Q1))}, + coupling_out=list(np.arange(Q1, Q1+Q2))) + comp3 = ComponentSpec(f3, name='Plume', truth_alpha=alpha, exo_in=list(np.arange(0, D1)), max_alpha=5, + coupling_in={'Thruster': list(np.arange(0, Q2))}, coupling_out=Q1+Q2, max_beta=(3,)*(D1+Q2)) + exo_vars = [UniformRV(0, 1) for i in range(D1+D2)] + coupling_vars = [UniformRV(0, 1) for i in range(Q1+Q2+1)] + surr = SystemSurrogate([comp1, comp2, comp3], exo_vars, coupling_vars, init_surr=False) + + # Test example + N = 5000 + x = np.random.rand(N, D1+D2) + y = f(x) + y_surr = surr(x, use_model='best') + l2_error = np.sqrt(np.mean((y_surr - y)**2, axis=0)) / np.sqrt(np.mean(y**2, axis=0)) + assert np.max(l2_error) < 1e-15 + + # Show coupling variable pdfs + if plots: + fig, ax = plt.subplots() + ls = ['-r', '--k', ':b'] + pts = np.linspace(0, 1, 100) + for i in range(3): + label_str = f'$\\rho(y_{{{i+1}}})$' + kernel = gaussian_kde(y_surr[:, i]) + # ax[i].hist(y_surr[:, i], density=True, bins=20, color='r', edgecolor='black', linewidth=1.2) + ax.plot(pts, kernel(pts), ls[i], label=label_str) + ax_default(ax, r'$y$', 'PDF', legend=True) + # fig.set_size_inches(9, 3) + fig.tight_layout() + plt.show() + + +def test_fpi(): + """Test fixed point iteration implementation against scipy fsolve.""" + f1 = lambda x: {'y': -x[..., 0:1]**3 + 2 * x[..., 1:2]**2} + f2 = lambda x: {'y': 3*x[..., 0:1]**2 + 4 * x[..., 1:2]**(-2)} + comp1 = ComponentSpec(f1, exo_in=[0], coupling_in={'m2': [0]}, coupling_out=[0], max_beta=(5,), name='m1') + comp2 = ComponentSpec(f2, exo_in=[0], coupling_in={'m1': [0]}, coupling_out=[1], max_beta=(5,), name='m2') + exo_vars = UniformRV(0, 4) + coupling_vars = [UniformRV(1, 10), UniformRV(1, 10)] + surr = SystemSurrogate([comp1, comp2], exo_vars, coupling_vars, init_surr=False) + + # Test on random x against scipy.fsolve + N = 100 + tol = 1e-12 + x0 = np.array([5.5, 5.5]) + exo = surr.sample_inputs(N) + y_surr = surr.predict(exo, use_model='best', anderson_mem=10, max_fpi_iter=200, fpi_tol=tol) # (N, 2) + nan_idx = list(np.any(np.isnan(y_surr), axis=-1).nonzero()[0]) + y_true = np.zeros((N, 2)) + bad_idx = [] + warnings.simplefilter('error') + for i in range(N): + def fun(x): + y1 = x[0] + y2 = x[1] + res1 = -exo[i, 0]**3 + 2*y2**2 - y1 + res2 = 3*exo[i, 0]**2 + 4*y1**(-2) - y2 + return [res1, res2] + + try: + y_true[i, :] = fsolve(fun, x0, xtol=tol) + except Exception as e: + bad_idx.append(i) + + y_surr = np.delete(y_surr, nan_idx + bad_idx, axis=0) + y_true = np.delete(y_true, nan_idx + bad_idx, axis=0) + l2_error = np.sqrt(np.mean((y_surr - y_true)**2, axis=0)) / np.sqrt(np.mean(y_true**2, axis=0)) + assert np.max(l2_error) < tol + + +def test_lls(): + """Test constrained linear least squares routine against scipy lapack.""" + X = 100 + Y = 100 + M = 10 + N = 10 + P = 1 + tol = 1e-8 + + A = np.random.rand(X, Y, M, N) + b = np.random.rand(X, Y, M, 1) + C = np.random.rand(X, Y, P, N) + d = np.random.rand(X, Y, P, 1) + + # custom solver + t1 = time.time() + alpha = np.squeeze(SystemSurrogate._constrained_lls(A, b, C, d), axis=-1) # (*, N) + t2 = time.time() + + # Built in scipy solver + alpha2 = np.zeros((X, Y, N)) + t3 = time.time() + for i in range(X): + for j in range(Y): + Ai = A[i, j, ...] + bi = b[i, j, ...] + Ci = C[i, j, ...] + di = d[i, j, ...] + ret = lapack.dgglse(Ai, Ci, bi, di) + alpha2[i, j, :] = ret[3] + t4 = time.time() + + # Results + diff = alpha - alpha2 + assert np.max(np.abs(diff)) < tol + print(f'Custom CLLS time: {t2-t1} s. Scipy time: {t4-t3} s.') diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 0000000..494d036 --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,57 @@ +import numpy as np +import matplotlib.pyplot as plt + +from amisc.utils import approx_hess, batch_normal_sample, ax_default, get_logger, ndscatter + + +def test_hessian(): + """Test the Hessian of $f(x, y) = 2x^2 + 3xy^3$""" + shape = (11, 12) + dim = 2 + theta0 = np.random.rand(*shape, dim) + fun = lambda theta: 2 * theta[..., 0:1]**2 + 3 * theta[..., 0:1] * theta[..., 1:2]**3 + + H_tilde = approx_hess(fun, theta0) + H_exact = np.empty(shape + (dim, dim)) + H_exact[..., 0, 0] = 4 + H_exact[..., 0, 1] = 9 * theta0[..., 1]**2 + H_exact[..., 1, 0] = 9 * theta0[..., 1]**2 + H_exact[..., 1, 1] = 18 * theta0[..., 0] * theta0[..., 1] + + assert np.allclose(H_tilde, H_exact, rtol=1e-3, atol=1e-3) + + +def test_sample(): + """Test 1d and 2d batch normal sampling""" + dim = 2 + shape = (4, 5) + N = 100000 + mean = np.random.rand(*shape, dim) + cov = np.eye(dim) * 0.01 + samples = batch_normal_sample(mean, cov, N) + assert samples.shape == (N, *shape, dim) + assert np.allclose(np.mean(samples, axis=0), mean, rtol=1e-3, atol=1e-3) + + mean = np.random.rand() + cov = 0.01 + samples = batch_normal_sample(mean, cov, N) + assert np.isclose(mean, np.mean(samples, axis=0), rtol=1e-3, atol=1e-3) + + +def test_logging_and_plotting(): + """Test logging and plotting utils""" + x = np.linspace(0, 1, 100) + fig, ax = plt.subplots() + ax.plot(x, x, label='Hey there:)') + ax_default(ax, 'X label here', 'Y label here', legend=True) + logger = get_logger('tester', stdout=True) + logger.info('Testing logger...') + + mean = np.array([5, -2]) + cov = np.array([[2, 0.4], [0.2, 1]]) + samples = np.random.multivariate_normal(mean, cov.T @ cov, size=100) + yt = samples[:, 0] + samples[:, 1] ** 2 + ysurr = yt + np.random.randn(*yt.shape) + err = np.abs(ysurr - yt) / np.abs(yt) + ndscatter(samples, labels=[r'$\alpha$', r'$\beta$', r'$\gamma$'], plot='scatter', cmap='plasma', + cb_norm='log', z=err)