diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 3e6786a..3b870c5 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -16,7 +16,7 @@ jobs: - name: Install dependencies run: | - pip install sphinx sphinx_rtd_theme myst_parser sphinx-autoapi ipython + pip install sphinx sphinx_rtd_theme myst_parser sphinx-autoapi ipython sphinx-gallery pip install .[dev] - name: List files @@ -45,3 +45,9 @@ jobs: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: docs/_build/html force_orphan: true + + - name: Upload documentation artifact + uses: actions/upload-artifact@v4 + with: + name: sphinx-html + path: docs/_build/html diff --git a/docs/Makefile b/docs/Makefile index c51c47d..c76bca6 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -34,6 +34,6 @@ latexpdf-noplot: clean: rm -rf $(BUILDDIR)/* - rm -rf source/auto_examples/ - rm -rf source/auto_tutorials/ - rm -rf source/api/ \ No newline at end of file +# rm -rf auto_examples/ + rm -rf auto_tutorials/ + rm -rf api/ \ No newline at end of file diff --git a/docs/auto_examples/auto_examples_jupyter.zip b/docs/auto_examples/auto_examples_jupyter.zip new file mode 100644 index 0000000..2daefde Binary files /dev/null and b/docs/auto_examples/auto_examples_jupyter.zip differ diff --git a/docs/auto_examples/auto_examples_python.zip b/docs/auto_examples/auto_examples_python.zip new file mode 100644 index 0000000..8098402 Binary files /dev/null and b/docs/auto_examples/auto_examples_python.zip differ diff --git a/docs/auto_examples/ex_genz_bcs.codeobj.json b/docs/auto_examples/ex_genz_bcs.codeobj.json new file mode 100644 index 0000000..98eb677 --- /dev/null +++ b/docs/auto_examples/ex_genz_bcs.codeobj.json @@ -0,0 +1,1578 @@ +{ + "GenzOscillatory": [ + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func.genz", + "module_short": "pytuq.func.genz", + "name": "GenzOscillatory" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func", + "module_short": "pytuq.func", + "name": "GenzOscillatory" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "GenzOscillatory" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func.genz", + "module_short": "pytuq.func.genz", + "name": "GenzBase" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func", + "module_short": "pytuq.func", + "name": "GenzBase" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "GenzBase" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func.func", + "module_short": "pytuq.func.func", + "name": "Function" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.func", + "module_short": "pytuq.func", + "name": "Function" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "Function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.func.genz", + "module_short": "pytuq.func.genz", + "name": "GenzOscillatory" + } + ], + "PCE": [ + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq.surrogates", + "module_short": "pytuq.surrogates", + "name": "PCE" + }, + { + "is_class": true, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCE" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE" + } + ], + "ax1": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes" + } + ], + "ax1.plot": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.plot" + } + ], + "ax1.set_xlabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_xlabel" + } + ], + "ax1.set_ylabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_ylabel" + } + ], + "ax2": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes" + } + ], + "ax2.plot": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.plot" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.plot" + } + ], + "ax2.set_xlabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_xlabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_xlabel" + } + ], + "ax2.set_ylabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Axes.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes._base", + "module_short": "matplotlib.axes._base", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.axes", + "module_short": "matplotlib.axes", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "_AxesBase.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.set_ylabel" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.set_ylabel" + } + ], + "copy.deepcopy": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "copy", + "module_short": "copy", + "name": "deepcopy" + } + ], + "eta_opt": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "float64" + } + ], + "etas": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "fig1": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "Figure" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Figure" + } + ], + "fig1.add_axes": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "Figure.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Figure.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "FigureBase.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "FigureBase.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.add_axes" + } + ], + "fig2": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "Figure" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Figure" + } + ], + "fig2.add_axes": [ + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "Figure.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Figure.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.figure", + "module_short": "matplotlib.figure", + "name": "FigureBase.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "FigureBase.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.artist", + "module_short": "matplotlib.artist", + "name": "Artist.add_axes" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib", + "module_short": "matplotlib", + "name": "Artist.add_axes" + } + ], + "func": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.func.genz", + "module_short": "pytuq.func.genz", + "name": "GenzOscillatory" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.func", + "module_short": "pytuq.func", + "name": "GenzOscillatory" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "GenzOscillatory" + } + ], + "func_dim": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "func_weights": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "kfold_cv": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + } + ], + "kfold_split": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + } + ], + "math.sqrt": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "math", + "module_short": "math", + "name": "sqrt" + } + ], + "n_trn": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "n_tst": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "noise_std": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "float" + } + ], + "np.argmin": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "argmin" + } + ], + "np.array": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "array" + } + ], + "np.array_split": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "array_split" + } + ], + "np.atleast_2d": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "atleast_2d" + } + ], + "np.concatenate": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "concatenate" + } + ], + "np.min": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "min" + } + ], + "np.power": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ufunc" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "power" + } + ], + "np.random.RandomState": [ + { + "is_class": true, + "is_explicit": false, + "module": "numpy.random.mtrand", + "module_short": "numpy.random", + "name": "RandomState" + }, + { + "is_class": true, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "RandomState" + }, + { + "is_class": true, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "RandomState" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "RandomState" + } + ], + "np.random.normal": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "normal" + } + ], + "np.random.seed": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "builtin_function_or_method" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "seed" + } + ], + "np.square": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ufunc" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "square" + } + ], + "np.squeeze": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "squeeze" + } + ], + "np.std": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "_ArrayFunctionDispatcher" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "std" + } + ], + "np.subtract": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ufunc" + }, + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "subtract" + } + ], + "optimize_eta": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + } + ], + "order": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "pce_surr": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates", + "module_short": "pytuq.surrogates", + "name": "PCE" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCE" + } + ], + "pce_surr.build": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE.build" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates", + "module_short": "pytuq.surrogates", + "name": "PCE.build" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCE.build" + } + ], + "pce_surr.evaluate": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE.evaluate" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates", + "module_short": "pytuq.surrogates", + "name": "PCE.evaluate" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCE.evaluate" + } + ], + "pce_surr.pcrv.mindices": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "pce_surr.pcrv.printInfo": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.rv.pcrv", + "module_short": "pytuq.rv.pcrv", + "name": "PCRV.printInfo" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.rv", + "module_short": "pytuq.rv", + "name": "PCRV.printInfo" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCRV.printInfo" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.rv.mrv", + "module_short": "pytuq.rv.mrv", + "name": "MRV.printInfo" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.rv", + "module_short": "pytuq.rv", + "name": "MRV.printInfo" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "MRV.printInfo" + } + ], + "pce_surr.set_training_data": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE.set_training_data" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates", + "module_short": "pytuq.surrogates", + "name": "PCE.set_training_data" + }, + { + "is_class": false, + "is_explicit": false, + "module": "pytuq", + "module_short": "pytuq", + "name": "PCE.set_training_data" + } + ], + "plt.errorbar": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "errorbar" + } + ], + "plt.figure": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "figure" + } + ], + "plt.legend": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "legend" + } + ], + "plt.plot": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "plot" + } + ], + "plt.savefig": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "savefig" + } + ], + "plt.subplots": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "subplots" + } + ], + "plt.tick_params": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "tick_params" + } + ], + "plt.xlabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "xlabel" + } + ], + "plt.xscale": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "xscale" + } + ], + "plt.ylabel": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "ylabel" + } + ], + "plt.yscale": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "matplotlib.pyplot", + "module_short": "matplotlib.pyplot", + "name": "yscale" + } + ], + "qmc.LatinHypercube": [ + { + "is_class": true, + "is_explicit": false, + "module": "scipy.stats._qmc", + "module_short": "scipy.stats._qmc", + "name": "LatinHypercube" + }, + { + "is_class": true, + "is_explicit": false, + "module": "scipy.stats", + "module_short": "scipy.stats", + "name": "LatinHypercube" + }, + { + "is_class": true, + "is_explicit": false, + "module": "scipy", + "module_short": "scipy", + "name": "LatinHypercube" + }, + { + "is_class": true, + "is_explicit": false, + "module": "scipy.stats._qmc", + "module_short": "scipy.stats._qmc", + "name": "QMCEngine" + }, + { + "is_class": true, + "is_explicit": false, + "module": "scipy.stats", + "module_short": "scipy.stats", + "name": "QMCEngine" + }, + { + "is_class": true, + "is_explicit": false, + "module": "scipy", + "module_short": "scipy", + "name": "QMCEngine" + }, + { + "is_class": true, + "is_explicit": false, + "module": "abc", + "module_short": "abc", + "name": "ABC" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats.qmc", + "module_short": "scipy.stats.qmc", + "name": "LatinHypercube" + } + ], + "rmse_trn": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "float64" + } + ], + "rmse_tst": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "float64" + } + ], + "rng": [ + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats._qmc", + "module_short": "scipy.stats._qmc", + "name": "LatinHypercube" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats", + "module_short": "scipy.stats", + "name": "LatinHypercube" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy", + "module_short": "scipy", + "name": "LatinHypercube" + } + ], + "rng.random": [ + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats._qmc", + "module_short": "scipy.stats._qmc", + "name": "LatinHypercube.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats", + "module_short": "scipy.stats", + "name": "LatinHypercube.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy", + "module_short": "scipy", + "name": "LatinHypercube.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats._qmc", + "module_short": "scipy.stats._qmc", + "name": "QMCEngine.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy.stats", + "module_short": "scipy.stats", + "name": "QMCEngine.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "scipy", + "module_short": "scipy", + "name": "QMCEngine.random" + }, + { + "is_class": false, + "is_explicit": false, + "module": "abc", + "module_short": "abc", + "name": "ABC.random" + } + ], + "rng_seed": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "int" + } + ], + "root_mean_squared_error": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "function" + }, + { + "is_class": false, + "is_explicit": false, + "module": "sklearn.metrics", + "module_short": "sklearn.metrics", + "name": "root_mean_squared_error" + } + ], + "x_trn": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "x_tst": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "y_trn": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "y_trn_approx": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "y_trn_mM": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ], + "y_tst": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ndarray" + } + ], + "y_tst_approx": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "dict" + } + ], + "y_tst_mM": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "list" + } + ] +} \ No newline at end of file diff --git a/docs/auto_examples/ex_genz_bcs.ipynb b/docs/auto_examples/ex_genz_bcs.ipynb new file mode 100644 index 0000000..3772d88 --- /dev/null +++ b/docs/auto_examples/ex_genz_bcs.ipynb @@ -0,0 +1,267 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n# Function Approximation with Sparse Regression\n\nIn this tutorial, we demonstrate how to approximate a function with sparse regression by constructing a Polynomial Chaos (PC) surrogate with \nBayesian compressive sensing (BCS). The function we will approximate here is the Genz Oscillatory function, defined as:\n\n\\begin{align}f(x) = \\cos\\left(2 \\pi s + \\sum_{i=1}^d w_i x_i\\right)\\end{align}\n\nWhere:\n\n- $s$: The shift parameter (``self.shift`` in the class).\n- $w_i$: The weights for each dimension (``self.weights`` in the class).\n- $x_i$: The input variables.\n- $d$: The dimensionality of the input $x$ (number of components in $x$).\n\nThrough 2 different build processes, we will construct two different PC surrogates to demonstrate effects of the BCS eta hyperparameter on model results. \nThe first build process will demonstrate most simply the construct-build-evaluate process when using BCS for our PC surrogate, along with a given a eta of 1e-10.\nThe second build process will select the most optimal eta for BCS through cross-validation algorithms (exposed here), which will soon be implemented in PyTUQ under-the-hood.\nAfterwards, we will evaluate both models on testing and training data, returning parity plots and a Root Mean Square Error for each evaluation.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import os\nimport sys\n\nimport numpy as np\nimport copy\nimport math\nimport pytuq.utils.funcbank as fcb\nfrom matplotlib import pyplot as plt\nfrom sklearn.metrics import root_mean_squared_error\n\nfrom pytuq.surrogates.pce import PCE\nfrom pytuq.utils.maps import scale01ToDom\nfrom pytuq.func.genz import GenzOscillatory" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Setting a random number generator seed:\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Random number generator\nfrom scipy.stats import qmc\nrng_seed = 43" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Constructing PC surrogate and generating data\nAfter importing GenzOscillatory from ``pytuq.func.genz``, we generate the Genz function below, along with training data and testing data with output noise. \nThis data and the corresponding Genz function will be used to create the same PC surrogate fitted in both examples: \n(1) The first surrogate will be fitted using BCS with a given eta, and (2) the second surrogate will be fitted using BCS with the most optimal eta.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Use Genz Oscillatory function in multiple dimensions for the true model\nfunc_dim = 4\nfunc_weights = [1.0/(i+1)**2 for i in range(func_dim)]\nfunc = GenzOscillatory(shift=0.25, weights=func_weights)\nnoise_std = 0.1\nrng = qmc.LatinHypercube(d=func_dim, seed=rng_seed)\n\n# Training data\nnp.random.seed(42)\nn_trn = 70\nx_trn = rng.random(n=n_trn) # random numbers in [0,1]^d\ny_trn = func(x_trn) + np.random.normal(0, noise_std, size = (n_trn,1))\n\n# Testing data\nn_tst = 10000\nx_tst = rng.random(n=n_tst) # random numbers in [0,1]^d\ny_tst = func(x_tst) + np.random.normal(0, noise_std, size = (n_tst,1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With a stochastic dimensionality of 4 (defined above) and a polynomial order of 4, we construct the PC surrogate that will be used in both builds. \nYou have the option of printing the PC surrogate's full basis, before BCS selects and retains the most significant PC coefficients to reduce the basis.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# (1) Construct a PC surrogate\norder = 4\npce_surr = PCE(func_dim, order, 'LU', verbose = 1)\n\n# Optional verbosity output:\nprint(\"Full Basis and Current Coefficients:\")\npce_surr.pcrv.printInfo()\nprint(\"Number of Basis Terms:\", len(pce_surr.pcrv.mindices[0]))\n\n# (1.5) Set training data\npce_surr.set_training_data(x_trn, y_trn[:,0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## BCS with default settings (default eta)\nHere, we call the PCE class method of ``build()`` to build the linear regression model used to fit the surrogate. \nWith the flag ``regression='bcs'``, we choose the BCS method for the fitting. A user-defined eta of 1e-10 is also passed in.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# (2) Build the linear regression object for fitting\npce_surr.build(regression='bcs', eta=1.e-10)\n\n# Optional verbosity output:\nprint(\"Retained Basis and Coefficients:\")\npce_surr.pcrv.printInfo()\nprint(\"Number of retained basis terms:\", len(pce_surr.pcrv.mindices[0]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After fitting, we evaluate the PCE using our training and testing data. To analyze the model's goodness of fit, \nwe calculate the root mean square error between the surrogate results and the training and testing data.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# (3) Evaluate the PC model\ny_trn_approx = pce_surr.evaluate(x_trn)\ny_tst_approx = pce_surr.evaluate(x_tst)\n\n# Evaluate goodness of fit with RMSE\nrmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx[\"Y_eval\"])\nprint(\"The training RMSE error in the PCE BCS approximation is %.2e\"%rmse_trn)\n\nrmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx[\"Y_eval\"])\nprint(\"The testing RMSE error in the PCE BCS approximation is %.2e\"%rmse_tst)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice above how the training RMSE error is almost half that of the testing RMSE error. This shows that our current model is overfitting, \nlearning the training data with noise too well. To address this issue, we can explore selecting a better eta for the BCS fitting.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Plot the surrogate model's output vs. the training data output\n\ny_trn_mM = [y_trn[:,0].min(),y_trn[:,0].max()]\n\n\nfig1 = plt.figure(figsize=(8,6))\nax1 = fig1.add_axes([0.15, 0.15, 0.80, 0.75])\n\n\nax1.plot(y_trn[:,0],y_trn_approx[\"Y_eval\"],\".\")\nax1.plot(y_trn_mM,y_trn_mM) # Diagonal line\n\nax1.set_xlabel(\"Train Data y\", size=16)\nax1.set_ylabel(\"Predicted y\", size=16);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Plot the surrogate model's output vs. the testing data output\n\ny_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()]\n\n\nfig2 = plt.figure(figsize=(8,6))\nax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75])\n\n\nax2.plot(y_tst[:,0],y_tst_approx[\"Y_eval\"],\".\")\nax2.plot(y_tst_mM,y_tst_mM) # Diagonal line\n\nax2.set_xlabel(\"Test Data y\", size=16)\nax2.set_ylabel(\"Predicted y\", size=16);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "BCS with optimal eta (found through cross-validation) \n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n In this section, we use the same PC surrogate, ``pce_surr``, for the second build. We call the PCE class method of ``build()`` \n to build the linear regression model used to fit the surrogate. With the flag ``regression='bcs'``, we choose the BCS method for the fitting. \n\n#############################################################################\n Instead of using a default eta, we call the cross-validation algorithm, ``optimize_eta()``, to choose the most optimal eta below. \n\n - With the flag ``plot=True``, the CV algorithm produces a graph of the training and testing (validation) data's RMSE values for each eta. The eta with the smallest RMSE for the validation data is the one chosen as the optimal eta.\n\n#############################################################################\n Before that, we expose the cross-validation algorithm ``optimize_eta`` and its two helper functions, ``kfold_split`` and ``kfold_cv`` that will be called in this section.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def kfold_split(nsamples,nfolds,seed=13):\n \"\"\"\n Return dictionary of training and testing pairs using k-fold cross-validation.\n\n Args:\n nsamples (int): Total number of training samples.\n nfolds (int): Number of folds to use for k-fold cross-validation.\n seed (int, optional): Random seed for reproducibility. Defaults to 13.\n\n Returns:\n dict: A dictionary where each key is the fold number (0 to nfolds-1) \n and each value is a dictionary with:\n - \"train index\" (np.ndarray): Indices of training samples.\n - \"val index\" (np.ndarray): Indices of validation samples.\n \"\"\"\n # Returns split data where each data is one fold left out\n KK=nfolds\n rn = np.random.RandomState(seed)\n\n # Creating a random permutation of the samples indices list\n indp=rn.permutation(nsamples)\n\n # Split the permuted indices into KK (or # folds) equal-sized chunks \n split_index=np.array_split(indp,KK)\n\n # Dictionary to hold the indices of the training and validation samples\n cvindices = {}\n\n # create testing and training folds\n for j in range(KK):\n # Iterating through the number of folds\n fold = j\n # Iterate through # folds, if i != fold number, \n newindex = [split_index[i] for i in range(len(split_index)) if i != (fold)]\n train_ind = np.array([],dtype='int64')\n for i in range(len(newindex)): train_ind = np.concatenate((train_ind,newindex[i]))\n test_ind = split_index[fold]\n cvindices[j] = {'train index': train_ind, 'val index': test_ind}\n\n return cvindices" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def kfold_cv(x,y,nfolds=3,seed=13):\n \"\"\"\n Splits data into training/testing pairs for kfold cross-val\n x is a data matrix of size n x d1, d1 is dim of input\n y is a data matrix of size n x d2, d2 is dim of output\n\n Args:\n x (np.ndarray): \n Input matrix with shape (n, d1) or 1D array with shape (n,).\n Each row is a sample; columns are input features.\n y (np.ndarray):\n Target array with shape (n,) for single-output, or (n, d2) for\n multi-output. If 1D, it is internally reshaped to (n, 1) before\n slicing; outputs are `np.squeeze`d per fold.\n nfolds (int, optional): Number of folds for cross-validation. Defaults to 3.\n seed (int, optional):\n Random seed for reproducible shuffling in `kfold_split`. Defaults to 13.\n \"\"\"\n if len(x.shape)>1:\n n,d1 = x.shape\n else:\n n=x.shape\n ynew = np.atleast_2d(y)\n if len(ynew) == 1: ynew = ynew.T # change to shape (n,1)\n _,d2 = ynew.shape\n cv_idx = kfold_split(n,nfolds,seed)\n\n kfold_data = {}\n for k in cv_idx.keys():\n kfold_data[k] = {\n 'xtrain': x[cv_idx[k]['train index']],\n 'xval': x[cv_idx[k]['val index']],\n 'ytrain': np.squeeze(ynew[cv_idx[k]['train index']]),\n 'yval': np.squeeze(ynew[cv_idx[k]['val index']])\n } # use squeeze to return 1d array\n\n # set train and test to the same if 1 fold\n if nfolds == 1:\n kfold_data[k]['xtrain'] = kfold_data[k]['xval']\n kfold_data[k]['ytrain'] = kfold_data[k]['yval']\n\n return kfold_data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def optimize_eta(pce, etas, nfolds, verbose=0, plot=False):\n \"\"\"\n Choose the optimum eta for Bayesian compressive sensing. Calculates the RMSE\n for each eta for a specified number of folds. Selects the eta with the lowest\n RMSE after averaging the RMSEs over the folds.\n\n Input:\n y: 1D numpy array (vector) with function, evaluated at the\n sample points [#samples,]\n x: N-dimensional NumPy array with sample points [#samples,\n #dimensions]\n etas: NumPy array or list with the threshold for stopping the\n algorithm. Smaller values retain more nonzero\n coefficients\n verbose: Flag for print statements\n plot: Flag for whether to generate a plot for eta optimization\n\n Output:\n eta_opt: Optimum eta\n\n \"\"\"\n # Split data in k folds -> Get dictionary of data split in training + testing folds\n kfold_data = kfold_cv(pce._x_train, pce._y_train, nfolds)\n\n # Each value has data for 1 fold. Each value is a list of the RMSEs for each possible eta in the fold. \n RMSE_list_per_fold_tr = [] \n\n # Same but for testing data\n RMSE_list_per_fold_test = []\n\n # Make a copy of the PCE object to run the cross-validation algorithm on\n pce_copy = copy.deepcopy(pce)\n pce_copy.verbose = 0\n\n # Loop through each fold\n for i in range(nfolds):\n\n # Get the training and validation data\n x_tr = kfold_data[i]['xtrain']\n y_tr = kfold_data[i]['ytrain']\n x_test = kfold_data[i]['xval']\n y_test = kfold_data[i]['yval']\n \n # As we conduct BCS for this fold with each separate eta, the RMSEs will be added to these lists\n RMSE_per_eta_tr = [] \n RMSE_per_eta_test = [] \n\n # Set the x and y training data for the copied PCE object\n pce_copy.set_training_data(x_tr, y_tr)\n\n # Loop through each eta\n for eta in etas:\n\n # Conduct the BCS fitting. The object is automatically updated with new multiindex and coefficients received from the fitting.\n cfs = pce_copy.build(regression = 'bcs', eta=eta)\n\n # Evaluate the PCE object at the training and validation points \n y_tr_eval = (pce_copy.evaluate(x_tr))['Y_eval']\n y_test_eval = (pce_copy.evaluate(x_test))['Y_eval']\n\n # Print statement for verbose flag\n if verbose > 1:\n print(\"Fold \" + str(i + 1) + \", eta \" + str(eta) + \", \" + str(len(cfs)) + \" terms retained out of a full basis of size \" + str(len(pce.pcrv.mindices[0])))\n \n # Calculate the RMSEs for the training and validation points.\n # Append the values into the list of etas per fold.\n MSE = np.square(np.subtract(y_tr, y_tr_eval)).mean()\n RMSE = math.sqrt(MSE)\n RMSE_per_eta_tr.append(RMSE)\n\n MSE = np.square(np.subtract(y_test, y_test_eval)).mean()\n RMSE = math.sqrt(MSE)\n RMSE_per_eta_test.append(RMSE)\n\n # Now, append the fold's list of RMSEs for each eta into the list carrying the lists for all folds \n RMSE_list_per_fold_tr.append(RMSE_per_eta_tr)\n RMSE_list_per_fold_test.append(RMSE_per_eta_test)\n\n # After compiling the RMSE data for each eta from all the folds, we find the eta with the lowest validation RMSE to be our optimal eta.\n # Compute the average and standard deviation of the training and testing RMSEs over the folds\n avg_RMSE_tr = np.array(RMSE_list_per_fold_tr).mean(axis=0)\n avg_RMSE_test = np.array(RMSE_list_per_fold_test).mean(axis=0)\n\n std_RMSE_tr = np.std(np.array(RMSE_list_per_fold_tr), axis=0)\n std_RMSE_test = np.std(np.array(RMSE_list_per_fold_test), axis=0)\n\n # Choose the eta with lowest RMSE across all folds' testing data\n eta_opt = etas[np.argmin(avg_RMSE_test)]\n\n # Plot RMSE vs. eta for training and testing RMSE\n if plot:\n\n fig, ax = plt.subplots(figsize=(10,10))\n\n plt.errorbar(etas, avg_RMSE_tr, xerr=None, yerr=std_RMSE_tr, linewidth=2, markersize=8, capsize=8, label=('Training'))\n plt.errorbar(etas, avg_RMSE_test, xerr=None, yerr=std_RMSE_test, linewidth=2, markersize=8, capsize=8, label=('Validation'))\n\n plt.plot(eta_opt, np.min(avg_RMSE_test), marker=\"o\", markersize=15, color='black', label=(\"Optimum\"))\n\n plt.xlabel(\"Eta\",fontsize=20)\n plt.ylabel(\"RMSE\",fontsize=20)\n\n # Change size of tick labels\n plt.tick_params(axis='both', labelsize=16)\n\n plt.xscale('log')\n plt.yscale('log')\n\n # Create legend\n plt.legend(loc='upper left')\n\n # Save\n plt.savefig('eta_opt.pdf', format='pdf', dpi=1200)\n\n return eta_opt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# We first create a list of possible etas to pass in: [1e-16, 1e-15, ... , 1e-2, 1e-1, 1]\netas = 1/np.power(10,[i for i in range(0,16)])\n\n# Then, we call the function to choose the optimal eta:\neta_opt = optimize_eta(pce_surr, etas, 10, plot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, with the optimal eta obtained, we run the fitting again. Then, we evaluate the PCE and produce a parity plot for the predicted output vs. the testing data. \nNotice that the larger eta, 10e-2, retained fewer basis terms (6) compared to the smaller user-defined eta of 10e-10 (which retained 20 basis terms).\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Build the linear regression object for fitting\npce_surr.build(regression='bcs', eta=eta_opt)\n\n# Optional verbosity output:\nprint(\"Retained Basis and Coefficients:\")\npce_surr.pcrv.printInfo()\nprint(\"Number of retained basis terms:\", len(pce_surr.pcrv.mindices[0]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Evaluate the PC model with training and testing data\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "y_trn_approx = pce_surr.evaluate(x_trn)\ny_tst_approx = pce_surr.evaluate(x_tst)\n\n# Evaluate goodness of fit with RMSE\nrmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx[\"Y_eval\"])\nprint(\"The training RMSE error in the PCE BCS approximation is %.2e\"%rmse_trn)\n\nrmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx[\"Y_eval\"])\nprint(\"The testing RMSE error in the PCE BCS approximation is %.2e\"%rmse_tst)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "While the training RMSE error was almost half that of the testing RMSE error for the first fitting, the RMSE errors here are much closer to each other in value. \nThis suggests that the model has more effectively generalized to the unseen data; a better eta has improved performance.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Plot the surrogate model's output vs. the testing data output\ny_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()]\n\n\nfig2 = plt.figure(figsize=(8,6))\nax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75])\n\n\nax2.plot(y_tst[:,0],y_tst_approx[\"Y_eval\"],\".\")\nax2.plot(y_tst_mM,y_tst_mM) # Diagonal line\n\nax2.set_xlabel(\"Test Data y\", size=16)\nax2.set_ylabel(\"Predicted y\", size=16); \n\n# sphinx_gallery_thumbnail_number = 2" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/docs/auto_examples/ex_genz_bcs.py b/docs/auto_examples/ex_genz_bcs.py new file mode 100644 index 0000000..1dd7b14 --- /dev/null +++ b/docs/auto_examples/ex_genz_bcs.py @@ -0,0 +1,420 @@ +r""" +Function Approximation with Sparse Regression +=============================================== + +In this tutorial, we demonstrate how to approximate a function with sparse regression by constructing a Polynomial Chaos (PC) surrogate with +Bayesian compressive sensing (BCS). The function we will approximate here is the Genz Oscillatory function, defined as: + +.. math:: + + f(x) = \cos\left(2 \pi s + \sum_{i=1}^d w_i x_i\right) + +Where: + +- :math:`s`: The shift parameter (``self.shift`` in the class). +- :math:`w_i`: The weights for each dimension (``self.weights`` in the class). +- :math:`x_i`: The input variables. +- :math:`d`: The dimensionality of the input :math:`x` (number of components in :math:`x`). + +Through 2 different build processes, we will construct two different PC surrogates to demonstrate effects of the BCS eta hyperparameter on model results. +The first build process will demonstrate most simply the construct-build-evaluate process when using BCS for our PC surrogate, along with a given a eta of 1e-10. +The second build process will select the most optimal eta for BCS through cross-validation algorithms (exposed here), which will soon be implemented in PyTUQ under-the-hood. +Afterwards, we will evaluate both models on testing and training data, returning parity plots and a Root Mean Square Error for each evaluation. +""" +# %% + +import os +import sys + +import numpy as np +import copy +import math +import pytuq.utils.funcbank as fcb +from matplotlib import pyplot as plt +from sklearn.metrics import root_mean_squared_error + +from pytuq.surrogates.pce import PCE +from pytuq.utils.maps import scale01ToDom +from pytuq.func.genz import GenzOscillatory + +# %% +# Setting a random number generator seed: + +# Random number generator +from scipy.stats import qmc +rng_seed = 43 + +# %% +# Constructing PC surrogate and generating data +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# After importing GenzOscillatory from ``pytuq.func.genz``, we generate the Genz function below, along with training data and testing data with output noise. +# This data and the corresponding Genz function will be used to create the same PC surrogate fitted in both examples: +# (1) The first surrogate will be fitted using BCS with a given eta, and (2) the second surrogate will be fitted using BCS with the most optimal eta. + +# Use Genz Oscillatory function in multiple dimensions for the true model +func_dim = 4 +func_weights = [1.0/(i+1)**2 for i in range(func_dim)] +func = GenzOscillatory(shift=0.25, weights=func_weights) +noise_std = 0.1 +rng = qmc.LatinHypercube(d=func_dim, seed=rng_seed) + +# Training data +np.random.seed(42) +n_trn = 70 +x_trn = rng.random(n=n_trn) # random numbers in [0,1]^d +y_trn = func(x_trn) + np.random.normal(0, noise_std, size = (n_trn,1)) + +# Testing data +n_tst = 10000 +x_tst = rng.random(n=n_tst) # random numbers in [0,1]^d +y_tst = func(x_tst) + np.random.normal(0, noise_std, size = (n_tst,1)) + +# %% +# With a stochastic dimensionality of 4 (defined above) and a polynomial order of 4, we construct the PC surrogate that will be used in both builds. +# You have the option of printing the PC surrogate's full basis, before BCS selects and retains the most significant PC coefficients to reduce the basis. + +# (1) Construct a PC surrogate +order = 4 +pce_surr = PCE(func_dim, order, 'LU', verbose = 1) + +# Optional verbosity output: +print("Full Basis and Current Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of Basis Terms:", len(pce_surr.pcrv.mindices[0])) + +# (1.5) Set training data +pce_surr.set_training_data(x_trn, y_trn[:,0]) + +# %% +# BCS with default settings (default eta) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Here, we call the PCE class method of ``build()`` to build the linear regression model used to fit the surrogate. +# With the flag ``regression='bcs'``, we choose the BCS method for the fitting. A user-defined eta of 1e-10 is also passed in. + +# (2) Build the linear regression object for fitting +pce_surr.build(regression='bcs', eta=1.e-10) + +# Optional verbosity output: +print("Retained Basis and Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + +# %% +# After fitting, we evaluate the PCE using our training and testing data. To analyze the model's goodness of fit, +# we calculate the root mean square error between the surrogate results and the training and testing data. + +# (3) Evaluate the PC model +y_trn_approx = pce_surr.evaluate(x_trn) +y_tst_approx = pce_surr.evaluate(x_tst) + +# Evaluate goodness of fit with RMSE +rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) +print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + +rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) +print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + +# %% +# Notice above how the training RMSE error is almost half that of the testing RMSE error. This shows that our current model is overfitting, +# learning the training data with noise too well. To address this issue, we can explore selecting a better eta for the BCS fitting. + +# Plot the surrogate model's output vs. the training data output + +y_trn_mM = [y_trn[:,0].min(),y_trn[:,0].max()] + + +fig1 = plt.figure(figsize=(8,6)) +ax1 = fig1.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax1.plot(y_trn[:,0],y_trn_approx["Y_eval"],".") +ax1.plot(y_trn_mM,y_trn_mM) # Diagonal line + +ax1.set_xlabel("Train Data y", size=16) +ax1.set_ylabel("Predicted y", size=16); + +# %% + +# Plot the surrogate model's output vs. the testing data output + +y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + +fig2 = plt.figure(figsize=(8,6)) +ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") +ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + +ax2.set_xlabel("Test Data y", size=16) +ax2.set_ylabel("Predicted y", size=16); + +# %% +# BCS with optimal eta (found through cross-validation) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# In this section, we use the same PC surrogate, ``pce_surr``, for the second build. We call the PCE class method of ``build()`` +# to build the linear regression model used to fit the surrogate. With the flag ``regression='bcs'``, we choose the BCS method for the fitting. +# +############################################################################## +# Instead of using a default eta, we call the cross-validation algorithm, ``optimize_eta()``, to choose the most optimal eta below. +# +# - With the flag ``plot=True``, the CV algorithm produces a graph of the training and testing (validation) data's RMSE values for each eta. The eta with the smallest RMSE for the validation data is the one chosen as the optimal eta. +# +############################################################################## +# Before that, we expose the cross-validation algorithm ``optimize_eta`` and its two helper functions, ``kfold_split`` and ``kfold_cv`` that will be called in this section. + +def kfold_split(nsamples,nfolds,seed=13): + """ + Return dictionary of training and testing pairs using k-fold cross-validation. + + Args: + nsamples (int): Total number of training samples. + nfolds (int): Number of folds to use for k-fold cross-validation. + seed (int, optional): Random seed for reproducibility. Defaults to 13. + + Returns: + dict: A dictionary where each key is the fold number (0 to nfolds-1) + and each value is a dictionary with: + - "train index" (np.ndarray): Indices of training samples. + - "val index" (np.ndarray): Indices of validation samples. + """ + # Returns split data where each data is one fold left out + KK=nfolds + rn = np.random.RandomState(seed) + + # Creating a random permutation of the samples indices list + indp=rn.permutation(nsamples) + + # Split the permuted indices into KK (or # folds) equal-sized chunks + split_index=np.array_split(indp,KK) + + # Dictionary to hold the indices of the training and validation samples + cvindices = {} + + # create testing and training folds + for j in range(KK): + # Iterating through the number of folds + fold = j + # Iterate through # folds, if i != fold number, + newindex = [split_index[i] for i in range(len(split_index)) if i != (fold)] + train_ind = np.array([],dtype='int64') + for i in range(len(newindex)): train_ind = np.concatenate((train_ind,newindex[i])) + test_ind = split_index[fold] + cvindices[j] = {'train index': train_ind, 'val index': test_ind} + + return cvindices + +# %% + +def kfold_cv(x,y,nfolds=3,seed=13): + """ + Splits data into training/testing pairs for kfold cross-val + x is a data matrix of size n x d1, d1 is dim of input + y is a data matrix of size n x d2, d2 is dim of output + + Args: + x (np.ndarray): + Input matrix with shape (n, d1) or 1D array with shape (n,). + Each row is a sample; columns are input features. + y (np.ndarray): + Target array with shape (n,) for single-output, or (n, d2) for + multi-output. If 1D, it is internally reshaped to (n, 1) before + slicing; outputs are `np.squeeze`d per fold. + nfolds (int, optional): Number of folds for cross-validation. Defaults to 3. + seed (int, optional): + Random seed for reproducible shuffling in `kfold_split`. Defaults to 13. + """ + if len(x.shape)>1: + n,d1 = x.shape + else: + n=x.shape + ynew = np.atleast_2d(y) + if len(ynew) == 1: ynew = ynew.T # change to shape (n,1) + _,d2 = ynew.shape + cv_idx = kfold_split(n,nfolds,seed) + + kfold_data = {} + for k in cv_idx.keys(): + kfold_data[k] = { + 'xtrain': x[cv_idx[k]['train index']], + 'xval': x[cv_idx[k]['val index']], + 'ytrain': np.squeeze(ynew[cv_idx[k]['train index']]), + 'yval': np.squeeze(ynew[cv_idx[k]['val index']]) + } # use squeeze to return 1d array + + # set train and test to the same if 1 fold + if nfolds == 1: + kfold_data[k]['xtrain'] = kfold_data[k]['xval'] + kfold_data[k]['ytrain'] = kfold_data[k]['yval'] + + return kfold_data + +# %% +def optimize_eta(pce, etas, nfolds, verbose=0, plot=False): + """ + Choose the optimum eta for Bayesian compressive sensing. Calculates the RMSE + for each eta for a specified number of folds. Selects the eta with the lowest + RMSE after averaging the RMSEs over the folds. + + Input: + y: 1D numpy array (vector) with function, evaluated at the + sample points [#samples,] + x: N-dimensional NumPy array with sample points [#samples, + #dimensions] + etas: NumPy array or list with the threshold for stopping the + algorithm. Smaller values retain more nonzero + coefficients + verbose: Flag for print statements + plot: Flag for whether to generate a plot for eta optimization + + Output: + eta_opt: Optimum eta + + """ + # Split data in k folds -> Get dictionary of data split in training + testing folds + kfold_data = kfold_cv(pce._x_train, pce._y_train, nfolds) + + # Each value has data for 1 fold. Each value is a list of the RMSEs for each possible eta in the fold. + RMSE_list_per_fold_tr = [] + + # Same but for testing data + RMSE_list_per_fold_test = [] + + # Make a copy of the PCE object to run the cross-validation algorithm on + pce_copy = copy.deepcopy(pce) + pce_copy.verbose = 0 + + # Loop through each fold + for i in range(nfolds): + + # Get the training and validation data + x_tr = kfold_data[i]['xtrain'] + y_tr = kfold_data[i]['ytrain'] + x_test = kfold_data[i]['xval'] + y_test = kfold_data[i]['yval'] + + # As we conduct BCS for this fold with each separate eta, the RMSEs will be added to these lists + RMSE_per_eta_tr = [] + RMSE_per_eta_test = [] + + # Set the x and y training data for the copied PCE object + pce_copy.set_training_data(x_tr, y_tr) + + # Loop through each eta + for eta in etas: + + # Conduct the BCS fitting. The object is automatically updated with new multiindex and coefficients received from the fitting. + cfs = pce_copy.build(regression = 'bcs', eta=eta) + + # Evaluate the PCE object at the training and validation points + y_tr_eval = (pce_copy.evaluate(x_tr))['Y_eval'] + y_test_eval = (pce_copy.evaluate(x_test))['Y_eval'] + + # Print statement for verbose flag + if verbose > 1: + print("Fold " + str(i + 1) + ", eta " + str(eta) + ", " + str(len(cfs)) + " terms retained out of a full basis of size " + str(len(pce.pcrv.mindices[0]))) + + # Calculate the RMSEs for the training and validation points. + # Append the values into the list of etas per fold. + MSE = np.square(np.subtract(y_tr, y_tr_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_tr.append(RMSE) + + MSE = np.square(np.subtract(y_test, y_test_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_test.append(RMSE) + + # Now, append the fold's list of RMSEs for each eta into the list carrying the lists for all folds + RMSE_list_per_fold_tr.append(RMSE_per_eta_tr) + RMSE_list_per_fold_test.append(RMSE_per_eta_test) + + # After compiling the RMSE data for each eta from all the folds, we find the eta with the lowest validation RMSE to be our optimal eta. + # Compute the average and standard deviation of the training and testing RMSEs over the folds + avg_RMSE_tr = np.array(RMSE_list_per_fold_tr).mean(axis=0) + avg_RMSE_test = np.array(RMSE_list_per_fold_test).mean(axis=0) + + std_RMSE_tr = np.std(np.array(RMSE_list_per_fold_tr), axis=0) + std_RMSE_test = np.std(np.array(RMSE_list_per_fold_test), axis=0) + + # Choose the eta with lowest RMSE across all folds' testing data + eta_opt = etas[np.argmin(avg_RMSE_test)] + + # Plot RMSE vs. eta for training and testing RMSE + if plot: + + fig, ax = plt.subplots(figsize=(10,10)) + + plt.errorbar(etas, avg_RMSE_tr, xerr=None, yerr=std_RMSE_tr, linewidth=2, markersize=8, capsize=8, label=('Training')) + plt.errorbar(etas, avg_RMSE_test, xerr=None, yerr=std_RMSE_test, linewidth=2, markersize=8, capsize=8, label=('Validation')) + + plt.plot(eta_opt, np.min(avg_RMSE_test), marker="o", markersize=15, color='black', label=("Optimum")) + + plt.xlabel("Eta",fontsize=20) + plt.ylabel("RMSE",fontsize=20) + + # Change size of tick labels + plt.tick_params(axis='both', labelsize=16) + + plt.xscale('log') + plt.yscale('log') + + # Create legend + plt.legend(loc='upper left') + + # Save + plt.savefig('eta_opt.pdf', format='pdf', dpi=1200) + + return eta_opt + +# %% + +# We first create a list of possible etas to pass in: [1e-16, 1e-15, ... , 1e-2, 1e-1, 1] +etas = 1/np.power(10,[i for i in range(0,16)]) + +# Then, we call the function to choose the optimal eta: +eta_opt = optimize_eta(pce_surr, etas, 10, plot=True) + +# %% +# Now, with the optimal eta obtained, we run the fitting again. Then, we evaluate the PCE and produce a parity plot for the predicted output vs. the testing data. +# Notice that the larger eta, 10e-2, retained fewer basis terms (6) compared to the smaller user-defined eta of 10e-10 (which retained 20 basis terms). + +# Build the linear regression object for fitting +pce_surr.build(regression='bcs', eta=eta_opt) + +# Optional verbosity output: +print("Retained Basis and Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + +# %% +# Evaluate the PC model with training and testing data +y_trn_approx = pce_surr.evaluate(x_trn) +y_tst_approx = pce_surr.evaluate(x_tst) + +# Evaluate goodness of fit with RMSE +rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) +print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + +rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) +print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + +# %% +# While the training RMSE error was almost half that of the testing RMSE error for the first fitting, the RMSE errors here are much closer to each other in value. +# This suggests that the model has more effectively generalized to the unseen data; a better eta has improved performance. + +# Plot the surrogate model's output vs. the testing data output +y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + +fig2 = plt.figure(figsize=(8,6)) +ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") +ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + +ax2.set_xlabel("Test Data y", size=16) +ax2.set_ylabel("Predicted y", size=16); + +# sphinx_gallery_thumbnail_number = 2 \ No newline at end of file diff --git a/docs/auto_examples/ex_genz_bcs.py.md5 b/docs/auto_examples/ex_genz_bcs.py.md5 new file mode 100644 index 0000000..1ae0028 --- /dev/null +++ b/docs/auto_examples/ex_genz_bcs.py.md5 @@ -0,0 +1 @@ +f4d64dfb6a6e4bd8172e2db42886d735 \ No newline at end of file diff --git a/docs/auto_examples/ex_genz_bcs.rst b/docs/auto_examples/ex_genz_bcs.rst new file mode 100644 index 0000000..d32031f --- /dev/null +++ b/docs/auto_examples/ex_genz_bcs.rst @@ -0,0 +1,844 @@ + +.. DO NOT EDIT. +.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. +.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: +.. "auto_examples/ex_genz_bcs.py" +.. LINE NUMBERS ARE GIVEN BELOW. + +.. only:: html + + .. note:: + :class: sphx-glr-download-link-note + + :ref:`Go to the end ` + to download the full example code. + +.. rst-class:: sphx-glr-example-title + +.. _sphx_glr_auto_examples_ex_genz_bcs.py: + + +Function Approximation with Sparse Regression +=============================================== + +In this tutorial, we demonstrate how to approximate a function with sparse regression by constructing a Polynomial Chaos (PC) surrogate with +Bayesian compressive sensing (BCS). The function we will approximate here is the Genz Oscillatory function, defined as: + +.. math:: + + f(x) = \cos\left(2 \pi s + \sum_{i=1}^d w_i x_i\right) + +Where: + +- :math:`s`: The shift parameter (``self.shift`` in the class). +- :math:`w_i`: The weights for each dimension (``self.weights`` in the class). +- :math:`x_i`: The input variables. +- :math:`d`: The dimensionality of the input :math:`x` (number of components in :math:`x`). + +Through 2 different build processes, we will construct two different PC surrogates to demonstrate effects of the BCS eta hyperparameter on model results. +The first build process will demonstrate most simply the construct-build-evaluate process when using BCS for our PC surrogate, along with a given a eta of 1e-10. +The second build process will select the most optimal eta for BCS through cross-validation algorithms (exposed here), which will soon be implemented in PyTUQ under-the-hood. +Afterwards, we will evaluate both models on testing and training data, returning parity plots and a Root Mean Square Error for each evaluation. + +.. GENERATED FROM PYTHON SOURCE LINES 25-40 + +.. code-block:: Python + + + import os + import sys + + import numpy as np + import copy + import math + import pytuq.utils.funcbank as fcb + from matplotlib import pyplot as plt + from sklearn.metrics import root_mean_squared_error + + from pytuq.surrogates.pce import PCE + from pytuq.utils.maps import scale01ToDom + from pytuq.func.genz import GenzOscillatory + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 41-42 + +Setting a random number generator seed: + +.. GENERATED FROM PYTHON SOURCE LINES 42-47 + +.. code-block:: Python + + + # Random number generator + from scipy.stats import qmc + rng_seed = 43 + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 48-53 + +Constructing PC surrogate and generating data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +After importing GenzOscillatory from ``pytuq.func.genz``, we generate the Genz function below, along with training data and testing data with output noise. +This data and the corresponding Genz function will be used to create the same PC surrogate fitted in both examples: +(1) The first surrogate will be fitted using BCS with a given eta, and (2) the second surrogate will be fitted using BCS with the most optimal eta. + +.. GENERATED FROM PYTHON SOURCE LINES 53-72 + +.. code-block:: Python + + + # Use Genz Oscillatory function in multiple dimensions for the true model + func_dim = 4 + func_weights = [1.0/(i+1)**2 for i in range(func_dim)] + func = GenzOscillatory(shift=0.25, weights=func_weights) + noise_std = 0.1 + rng = qmc.LatinHypercube(d=func_dim, seed=rng_seed) + + # Training data + np.random.seed(42) + n_trn = 70 + x_trn = rng.random(n=n_trn) # random numbers in [0,1]^d + y_trn = func(x_trn) + np.random.normal(0, noise_std, size = (n_trn,1)) + + # Testing data + n_tst = 10000 + x_tst = rng.random(n=n_tst) # random numbers in [0,1]^d + y_tst = func(x_tst) + np.random.normal(0, noise_std, size = (n_tst,1)) + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 73-75 + +With a stochastic dimensionality of 4 (defined above) and a polynomial order of 4, we construct the PC surrogate that will be used in both builds. +You have the option of printing the PC surrogate's full basis, before BCS selects and retains the most significant PC coefficients to reduce the basis. + +.. GENERATED FROM PYTHON SOURCE LINES 75-88 + +.. code-block:: Python + + + # (1) Construct a PC surrogate + order = 4 + pce_surr = PCE(func_dim, order, 'LU', verbose = 1) + + # Optional verbosity output: + print("Full Basis and Current Coefficients:") + pce_surr.pcrv.printInfo() + print("Number of Basis Terms:", len(pce_surr.pcrv.mindices[0])) + + # (1.5) Set training data + pce_surr.set_training_data(x_trn, y_trn[:,0]) + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + Constructed PC Surrogate with the following attributes: + ['LU', 'LU', 'LU', 'LU'] PC Random Variable(pdim=1, sdim=4) + + Full Basis and Current Coefficients: + [[0 0 0 0] + [1 0 0 0] + [0 1 0 0] + [0 0 1 0] + [0 0 0 1] + [2 0 0 0] + [1 1 0 0] + [1 0 1 0] + [1 0 0 1] + [0 2 0 0] + [0 1 1 0] + [0 1 0 1] + [0 0 2 0] + [0 0 1 1] + [0 0 0 2] + [3 0 0 0] + [2 1 0 0] + [2 0 1 0] + [2 0 0 1] + [1 2 0 0] + [1 1 1 0] + [1 1 0 1] + [1 0 2 0] + [1 0 1 1] + [1 0 0 2] + [0 3 0 0] + [0 2 1 0] + [0 2 0 1] + [0 1 2 0] + [0 1 1 1] + [0 1 0 2] + [0 0 3 0] + [0 0 2 1] + [0 0 1 2] + [0 0 0 3] + [4 0 0 0] + [3 1 0 0] + [3 0 1 0] + [3 0 0 1] + [2 2 0 0] + [2 1 1 0] + [2 1 0 1] + [2 0 2 0] + [2 0 1 1] + [2 0 0 2] + [1 3 0 0] + [1 2 1 0] + [1 2 0 1] + [1 1 2 0] + [1 1 1 1] + [1 1 0 2] + [1 0 3 0] + [1 0 2 1] + [1 0 1 2] + [1 0 0 3] + [0 4 0 0] + [0 3 1 0] + [0 3 0 1] + [0 2 2 0] + [0 2 1 1] + [0 2 0 2] + [0 1 3 0] + [0 1 2 1] + [0 1 1 2] + [0 1 0 3] + [0 0 4 0] + [0 0 3 1] + [0 0 2 2] + [0 0 1 3] + [0 0 0 4]] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. + 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. + 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] + Number of Basis Terms: 70 + + + + +.. GENERATED FROM PYTHON SOURCE LINES 89-93 + +BCS with default settings (default eta) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +Here, we call the PCE class method of ``build()`` to build the linear regression model used to fit the surrogate. +With the flag ``regression='bcs'``, we choose the BCS method for the fitting. A user-defined eta of 1e-10 is also passed in. + +.. GENERATED FROM PYTHON SOURCE LINES 93-102 + +.. code-block:: Python + + + # (2) Build the linear regression object for fitting + pce_surr.build(regression='bcs', eta=1.e-10) + + # Optional verbosity output: + print("Retained Basis and Coefficients:") + pce_surr.pcrv.printInfo() + print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + Regression method: bcs + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Retained Basis and Coefficients: + [[2 0 0 0] + [1 0 3 0] + [0 2 0 2] + [0 0 1 2] + [4 0 0 0] + [3 1 0 0] + [1 1 0 1] + [2 0 0 2] + [2 0 1 1] + [0 0 0 4] + [1 0 0 3] + [0 0 0 1] + [1 1 0 0] + [1 0 0 1] + [0 3 1 0] + [0 1 0 2] + [0 1 1 2] + [2 0 2 0] + [0 1 0 1] + [0 0 0 0]] [ 0.11844696 -0.1286659 -0.23890894 0.20561299 -0.65563148 1.41773821 + -0.20218904 -0.30625401 0.40503892 -0.27250245 0.99544673 -0.90463019 + -2.75925351 -0.66386396 0.13068211 -0.96868699 -0.65718534 -0.27253516 + 3.07988105 0.12550357] + Number of retained basis terms: 20 + + + + +.. GENERATED FROM PYTHON SOURCE LINES 103-105 + +After fitting, we evaluate the PCE using our training and testing data. To analyze the model's goodness of fit, +we calculate the root mean square error between the surrogate results and the training and testing data. + +.. GENERATED FROM PYTHON SOURCE LINES 105-117 + +.. code-block:: Python + + + # (3) Evaluate the PC model + y_trn_approx = pce_surr.evaluate(x_trn) + y_tst_approx = pce_surr.evaluate(x_tst) + + # Evaluate goodness of fit with RMSE + rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) + print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + + rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) + print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + The training RMSE error in the PCE BCS approximation is 1.72e-01 + The testing RMSE error in the PCE BCS approximation is 2.73e-01 + + + + +.. GENERATED FROM PYTHON SOURCE LINES 118-120 + +Notice above how the training RMSE error is almost half that of the testing RMSE error. This shows that our current model is overfitting, +learning the training data with noise too well. To address this issue, we can explore selecting a better eta for the BCS fitting. + +.. GENERATED FROM PYTHON SOURCE LINES 120-136 + +.. code-block:: Python + + + # Plot the surrogate model's output vs. the training data output + + y_trn_mM = [y_trn[:,0].min(),y_trn[:,0].max()] + + + fig1 = plt.figure(figsize=(8,6)) + ax1 = fig1.add_axes([0.15, 0.15, 0.80, 0.75]) + + + ax1.plot(y_trn[:,0],y_trn_approx["Y_eval"],".") + ax1.plot(y_trn_mM,y_trn_mM) # Diagonal line + + ax1.set_xlabel("Train Data y", size=16) + ax1.set_ylabel("Predicted y", size=16); + + + + +.. image-sg:: /auto_examples/images/sphx_glr_ex_genz_bcs_001.png + :alt: ex genz bcs + :srcset: /auto_examples/images/sphx_glr_ex_genz_bcs_001.png + :class: sphx-glr-single-img + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + + Text(71.09722222222221, 0.5, 'Predicted y') + + + +.. GENERATED FROM PYTHON SOURCE LINES 137-153 + +.. code-block:: Python + + + # Plot the surrogate model's output vs. the testing data output + + y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + + fig2 = plt.figure(figsize=(8,6)) + ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + + ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") + ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + + ax2.set_xlabel("Test Data y", size=16) + ax2.set_ylabel("Predicted y", size=16); + + + + +.. image-sg:: /auto_examples/images/sphx_glr_ex_genz_bcs_002.png + :alt: ex genz bcs + :srcset: /auto_examples/images/sphx_glr_ex_genz_bcs_002.png + :class: sphx-glr-single-img + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + + Text(71.09722222222221, 0.5, 'Predicted y') + + + +.. GENERATED FROM PYTHON SOURCE LINES 154-166 + +BCS with optimal eta (found through cross-validation) + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + In this section, we use the same PC surrogate, ``pce_surr``, for the second build. We call the PCE class method of ``build()`` + to build the linear regression model used to fit the surrogate. With the flag ``regression='bcs'``, we choose the BCS method for the fitting. + +############################################################################# + Instead of using a default eta, we call the cross-validation algorithm, ``optimize_eta()``, to choose the most optimal eta below. + + - With the flag ``plot=True``, the CV algorithm produces a graph of the training and testing (validation) data's RMSE values for each eta. The eta with the smallest RMSE for the validation data is the one chosen as the optimal eta. + +############################################################################# + Before that, we expose the cross-validation algorithm ``optimize_eta`` and its two helper functions, ``kfold_split`` and ``kfold_cv`` that will be called in this section. + +.. GENERATED FROM PYTHON SOURCE LINES 166-208 + +.. code-block:: Python + + + def kfold_split(nsamples,nfolds,seed=13): + """ + Return dictionary of training and testing pairs using k-fold cross-validation. + + Args: + nsamples (int): Total number of training samples. + nfolds (int): Number of folds to use for k-fold cross-validation. + seed (int, optional): Random seed for reproducibility. Defaults to 13. + + Returns: + dict: A dictionary where each key is the fold number (0 to nfolds-1) + and each value is a dictionary with: + - "train index" (np.ndarray): Indices of training samples. + - "val index" (np.ndarray): Indices of validation samples. + """ + # Returns split data where each data is one fold left out + KK=nfolds + rn = np.random.RandomState(seed) + + # Creating a random permutation of the samples indices list + indp=rn.permutation(nsamples) + + # Split the permuted indices into KK (or # folds) equal-sized chunks + split_index=np.array_split(indp,KK) + + # Dictionary to hold the indices of the training and validation samples + cvindices = {} + + # create testing and training folds + for j in range(KK): + # Iterating through the number of folds + fold = j + # Iterate through # folds, if i != fold number, + newindex = [split_index[i] for i in range(len(split_index)) if i != (fold)] + train_ind = np.array([],dtype='int64') + for i in range(len(newindex)): train_ind = np.concatenate((train_ind,newindex[i])) + test_ind = split_index[fold] + cvindices[j] = {'train index': train_ind, 'val index': test_ind} + + return cvindices + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 209-253 + +.. code-block:: Python + + + def kfold_cv(x,y,nfolds=3,seed=13): + """ + Splits data into training/testing pairs for kfold cross-val + x is a data matrix of size n x d1, d1 is dim of input + y is a data matrix of size n x d2, d2 is dim of output + + Args: + x (np.ndarray): + Input matrix with shape (n, d1) or 1D array with shape (n,). + Each row is a sample; columns are input features. + y (np.ndarray): + Target array with shape (n,) for single-output, or (n, d2) for + multi-output. If 1D, it is internally reshaped to (n, 1) before + slicing; outputs are `np.squeeze`d per fold. + nfolds (int, optional): Number of folds for cross-validation. Defaults to 3. + seed (int, optional): + Random seed for reproducible shuffling in `kfold_split`. Defaults to 13. + """ + if len(x.shape)>1: + n,d1 = x.shape + else: + n=x.shape + ynew = np.atleast_2d(y) + if len(ynew) == 1: ynew = ynew.T # change to shape (n,1) + _,d2 = ynew.shape + cv_idx = kfold_split(n,nfolds,seed) + + kfold_data = {} + for k in cv_idx.keys(): + kfold_data[k] = { + 'xtrain': x[cv_idx[k]['train index']], + 'xval': x[cv_idx[k]['val index']], + 'ytrain': np.squeeze(ynew[cv_idx[k]['train index']]), + 'yval': np.squeeze(ynew[cv_idx[k]['val index']]) + } # use squeeze to return 1d array + + # set train and test to the same if 1 fold + if nfolds == 1: + kfold_data[k]['xtrain'] = kfold_data[k]['xval'] + kfold_data[k]['ytrain'] = kfold_data[k]['yval'] + + return kfold_data + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 254-370 + +.. code-block:: Python + + def optimize_eta(pce, etas, nfolds, verbose=0, plot=False): + """ + Choose the optimum eta for Bayesian compressive sensing. Calculates the RMSE + for each eta for a specified number of folds. Selects the eta with the lowest + RMSE after averaging the RMSEs over the folds. + + Input: + y: 1D numpy array (vector) with function, evaluated at the + sample points [#samples,] + x: N-dimensional NumPy array with sample points [#samples, + #dimensions] + etas: NumPy array or list with the threshold for stopping the + algorithm. Smaller values retain more nonzero + coefficients + verbose: Flag for print statements + plot: Flag for whether to generate a plot for eta optimization + + Output: + eta_opt: Optimum eta + + """ + # Split data in k folds -> Get dictionary of data split in training + testing folds + kfold_data = kfold_cv(pce._x_train, pce._y_train, nfolds) + + # Each value has data for 1 fold. Each value is a list of the RMSEs for each possible eta in the fold. + RMSE_list_per_fold_tr = [] + + # Same but for testing data + RMSE_list_per_fold_test = [] + + # Make a copy of the PCE object to run the cross-validation algorithm on + pce_copy = copy.deepcopy(pce) + pce_copy.verbose = 0 + + # Loop through each fold + for i in range(nfolds): + + # Get the training and validation data + x_tr = kfold_data[i]['xtrain'] + y_tr = kfold_data[i]['ytrain'] + x_test = kfold_data[i]['xval'] + y_test = kfold_data[i]['yval'] + + # As we conduct BCS for this fold with each separate eta, the RMSEs will be added to these lists + RMSE_per_eta_tr = [] + RMSE_per_eta_test = [] + + # Set the x and y training data for the copied PCE object + pce_copy.set_training_data(x_tr, y_tr) + + # Loop through each eta + for eta in etas: + + # Conduct the BCS fitting. The object is automatically updated with new multiindex and coefficients received from the fitting. + cfs = pce_copy.build(regression = 'bcs', eta=eta) + + # Evaluate the PCE object at the training and validation points + y_tr_eval = (pce_copy.evaluate(x_tr))['Y_eval'] + y_test_eval = (pce_copy.evaluate(x_test))['Y_eval'] + + # Print statement for verbose flag + if verbose > 1: + print("Fold " + str(i + 1) + ", eta " + str(eta) + ", " + str(len(cfs)) + " terms retained out of a full basis of size " + str(len(pce.pcrv.mindices[0]))) + + # Calculate the RMSEs for the training and validation points. + # Append the values into the list of etas per fold. + MSE = np.square(np.subtract(y_tr, y_tr_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_tr.append(RMSE) + + MSE = np.square(np.subtract(y_test, y_test_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_test.append(RMSE) + + # Now, append the fold's list of RMSEs for each eta into the list carrying the lists for all folds + RMSE_list_per_fold_tr.append(RMSE_per_eta_tr) + RMSE_list_per_fold_test.append(RMSE_per_eta_test) + + # After compiling the RMSE data for each eta from all the folds, we find the eta with the lowest validation RMSE to be our optimal eta. + # Compute the average and standard deviation of the training and testing RMSEs over the folds + avg_RMSE_tr = np.array(RMSE_list_per_fold_tr).mean(axis=0) + avg_RMSE_test = np.array(RMSE_list_per_fold_test).mean(axis=0) + + std_RMSE_tr = np.std(np.array(RMSE_list_per_fold_tr), axis=0) + std_RMSE_test = np.std(np.array(RMSE_list_per_fold_test), axis=0) + + # Choose the eta with lowest RMSE across all folds' testing data + eta_opt = etas[np.argmin(avg_RMSE_test)] + + # Plot RMSE vs. eta for training and testing RMSE + if plot: + + fig, ax = plt.subplots(figsize=(10,10)) + + plt.errorbar(etas, avg_RMSE_tr, xerr=None, yerr=std_RMSE_tr, linewidth=2, markersize=8, capsize=8, label=('Training')) + plt.errorbar(etas, avg_RMSE_test, xerr=None, yerr=std_RMSE_test, linewidth=2, markersize=8, capsize=8, label=('Validation')) + + plt.plot(eta_opt, np.min(avg_RMSE_test), marker="o", markersize=15, color='black', label=("Optimum")) + + plt.xlabel("Eta",fontsize=20) + plt.ylabel("RMSE",fontsize=20) + + # Change size of tick labels + plt.tick_params(axis='both', labelsize=16) + + plt.xscale('log') + plt.yscale('log') + + # Create legend + plt.legend(loc='upper left') + + # Save + plt.savefig('eta_opt.pdf', format='pdf', dpi=1200) + + return eta_opt + + + + + + + + +.. GENERATED FROM PYTHON SOURCE LINES 371-378 + +.. code-block:: Python + + + # We first create a list of possible etas to pass in: [1e-16, 1e-15, ... , 1e-2, 1e-1, 1] + etas = 1/np.power(10,[i for i in range(0,16)]) + + # Then, we call the function to choose the optimal eta: + eta_opt = optimize_eta(pce_surr, etas, 10, plot=True) + + + + +.. image-sg:: /auto_examples/images/sphx_glr_ex_genz_bcs_003.png + :alt: ex genz bcs + :srcset: /auto_examples/images/sphx_glr_ex_genz_bcs_003.png + :class: sphx-glr-single-img + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + Warning: Sigma matrix has a negative diagonal element. Setting them to zero, but this may lead to inaccuracies. + + + + +.. GENERATED FROM PYTHON SOURCE LINES 379-381 + +Now, with the optimal eta obtained, we run the fitting again. Then, we evaluate the PCE and produce a parity plot for the predicted output vs. the testing data. +Notice that the larger eta, 10e-2, retained fewer basis terms (6) compared to the smaller user-defined eta of 10e-10 (which retained 20 basis terms). + +.. GENERATED FROM PYTHON SOURCE LINES 381-390 + +.. code-block:: Python + + + # Build the linear regression object for fitting + pce_surr.build(regression='bcs', eta=eta_opt) + + # Optional verbosity output: + print("Retained Basis and Coefficients:") + pce_surr.pcrv.printInfo() + print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + Regression method: bcs + Retained Basis and Coefficients: + [[1 0 0 0] + [2 0 0 0] + [0 1 0 1] + [1 0 3 0] + [0 2 0 2] + [0 0 1 2]] [-1.15692903 0.24292952 -0.29524043 -0.15503765 0.13276008 0.06118484] + Number of retained basis terms: 6 + + + + +.. GENERATED FROM PYTHON SOURCE LINES 391-392 + +Evaluate the PC model with training and testing data + +.. GENERATED FROM PYTHON SOURCE LINES 392-402 + +.. code-block:: Python + + y_trn_approx = pce_surr.evaluate(x_trn) + y_tst_approx = pce_surr.evaluate(x_tst) + + # Evaluate goodness of fit with RMSE + rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) + print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + + rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) + print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + + + + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + The training RMSE error in the PCE BCS approximation is 8.11e-02 + The testing RMSE error in the PCE BCS approximation is 1.10e-01 + + + + +.. GENERATED FROM PYTHON SOURCE LINES 403-405 + +While the training RMSE error was almost half that of the testing RMSE error for the first fitting, the RMSE errors here are much closer to each other in value. +This suggests that the model has more effectively generalized to the unseen data; a better eta has improved performance. + +.. GENERATED FROM PYTHON SOURCE LINES 405-420 + +.. code-block:: Python + + + # Plot the surrogate model's output vs. the testing data output + y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + + fig2 = plt.figure(figsize=(8,6)) + ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + + ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") + ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + + ax2.set_xlabel("Test Data y", size=16) + ax2.set_ylabel("Predicted y", size=16); + + # sphinx_gallery_thumbnail_number = 2 + + +.. image-sg:: /auto_examples/images/sphx_glr_ex_genz_bcs_004.png + :alt: ex genz bcs + :srcset: /auto_examples/images/sphx_glr_ex_genz_bcs_004.png + :class: sphx-glr-single-img + + +.. rst-class:: sphx-glr-script-out + + .. code-block:: none + + + Text(71.09722222222221, 0.5, 'Predicted y') + + + + +.. rst-class:: sphx-glr-timing + + **Total running time of the script:** (0 minutes 5.760 seconds) + + +.. _sphx_glr_download_auto_examples_ex_genz_bcs.py: + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-example + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download Jupyter notebook: ex_genz_bcs.ipynb ` + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download Python source code: ex_genz_bcs.py ` + + .. container:: sphx-glr-download sphx-glr-download-zip + + :download:`Download zipped: ex_genz_bcs.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/docs/auto_examples/ex_genz_bcs.zip b/docs/auto_examples/ex_genz_bcs.zip new file mode 100644 index 0000000..b37d268 Binary files /dev/null and b/docs/auto_examples/ex_genz_bcs.zip differ diff --git a/docs/auto_examples/ex_nn.codeobj.json b/docs/auto_examples/ex_nn.codeobj.json new file mode 100644 index 0000000..8de0fd2 --- /dev/null +++ b/docs/auto_examples/ex_nn.codeobj.json @@ -0,0 +1,137 @@ +{ + "NN": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.nn", + "module_short": "pytuq.surrogates.nn", + "name": "NN" + } + ], + "Poly": [ + { + "is_class": false, + "is_explicit": false, + "module": "quinn.nns.rnet", + "module_short": "quinn.nns.rnet", + "name": "Poly" + } + ], + "Sine": [ + { + "is_class": false, + "is_explicit": false, + "module": "quinn.func.funcs", + "module_short": "quinn.func.funcs", + "name": "Sine" + } + ], + "__name__": [ + { + "is_class": false, + "is_explicit": false, + "module": "builtins", + "module_short": "builtins", + "name": "str" + } + ], + "myrc": [ + { + "is_class": false, + "is_explicit": false, + "module": "quinn.utils.plotting", + "module_short": "quinn.utils.plotting", + "name": "myrc" + } + ], + "np.array": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "array" + } + ], + "np.pi": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "pi" + } + ], + "np.random.rand": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "rand" + } + ], + "np.random.seed": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "seed" + } + ], + "np.tile": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "tile" + } + ], + "scale01ToDom": [ + { + "is_class": false, + "is_explicit": false, + "module": "quinn.utils.maps", + "module_short": "quinn.utils.maps", + "name": "scale01ToDom" + } + ], + "torch.cuda.is_available": [ + { + "is_class": false, + "is_explicit": false, + "module": "torch.cuda", + "module_short": "torch.cuda", + "name": "is_available" + } + ], + "torch.device": [ + { + "is_class": false, + "is_explicit": false, + "module": "torch", + "module_short": "torch", + "name": "device" + } + ], + "torch.double": [ + { + "is_class": false, + "is_explicit": false, + "module": "torch", + "module_short": "torch", + "name": "double" + } + ], + "torch.set_default_dtype": [ + { + "is_class": false, + "is_explicit": false, + "module": "torch", + "module_short": "torch", + "name": "set_default_dtype" + } + ] +} \ No newline at end of file diff --git a/docs/auto_examples/ex_nn.ipynb b/docs/auto_examples/ex_nn.ipynb new file mode 100644 index 0000000..070395d --- /dev/null +++ b/docs/auto_examples/ex_nn.ipynb @@ -0,0 +1,43 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n# Residual Neural Network Construction\n\nThis example demonstrates how to use the Neural Network wrapper class, ``pytuq.surrogates.nn`` with scalar valued functions.\nThe constructor of the NN class accepts in an optional dictionary, ``net_options``, to specify additional hyperparameters. \nThe ``build()`` and ``evaluate()`` functions similarly accept dictionaries and explicit keyword arguments during their respective function calls.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import sys\nimport torch\nimport numpy as np\n\nfrom pytuq.surrogates.nn import NN\nfrom quinn.solvers.nn_vi import NN_VI\nfrom quinn.nns.rnet import RNet, Poly\nfrom quinn.utils.plotting import myrc\nfrom quinn.utils.maps import scale01ToDom\nfrom quinn.func.funcs import Sine, Sine10, blundell\n\ndef main():\n \"\"\"Main function.\"\"\"\n torch.set_default_dtype(torch.double)\n myrc()\n\n #################################################################################\n #################################################################################\n\n # defaults to cuda:0 if available\n device_id='cuda:0'\n device = torch.device(device_id if torch.cuda.is_available() else 'cpu')\n print(\"Using device\",device)\n\n nall = 15 # total number of points\n trn_factor = 0.9 # which fraction of nall goes to training\n ntst = 13 # separate test set\n ndim = 1 # input dimensionality\n datanoise = 0.02 # Noise in the generated data\n true_model, nout = Sine, 1 # Scalar valued output example\n\n #################################################################################\n #################################################################################\n\n # Domain: defining range of input variable\n domain = np.tile(np.array([-np.pi, np.pi]), (ndim, 1))\n\n np.random.seed(111)\n\n # Generating x-y training, validation, and testing data\n xall = scale01ToDom(np.random.rand(nall, ndim), domain)\n if true_model is not None:\n yall = true_model(xall, datanoise=datanoise)\n\n if ntst > 0:\n np.random.seed(100)\n xtst = scale01ToDom(np.random.rand(ntst, ndim), domain)\n if true_model is not None:\n ytst = true_model(xtst, datanoise=datanoise)\n\n # (1) Initialize neural network with optional parameters for object instantiation\n net_options = {'wp_function': Poly(0),\n 'indim': ndim, \n 'outdim': nout,\n 'layer_pre': True,\n 'layer_post': True,\n 'biasorno': True,\n 'nonlin': True,\n 'mlp': False, \n 'final_layer': None,\n 'device': device,\n }\n\n # Pass in unpacked net_options dict to constructor through kwargs\n nnet = NN('RNet', 3, 3, **net_options)\n \n # (1.5) Split data into training and validation, assign through member functions\n ntrn = int(trn_factor * nall)\n indperm = range(nall) # np.random.permutation(range(nall))\n indtrn = indperm[:ntrn]\n indval = indperm[ntrn:]\n xtrn, xval = xall[indtrn, :], xall[indval, :]\n ytrn, yval = yall[indtrn, :], yall[indval, :]\n\n nnet.set_validation_data(xval, yval) # optional\n nnet.set_training_data(xtrn, ytrn)\n\n # (2) 1st build: Call build function (defaults to UQ method of variational inference) with optional parameters for fitting\n nnet.build(datanoise=datanoise, lrate=0.01, batch_size=None, nsam=1, nepochs=5000, verbose=False)\n\n results = nnet.evaluate(xtst, nsam = 100, msc = 2) # Return samples of predictions with variance + covariance\n print(\"Y_eval:\", results['Y_eval'], end=\"\\n\\n\") # Printing only samples of predictions\n\n # (3) 2nd build: Example of throwing an exception when passing in a network option through build()\n # fit_options = {'datanoise': 0.05,\n # 'outdim': 3,\n # 'lrate': 0.01,\n # 'batch_size': None,\n # 'nsam': 1,\n # 'nepochs': 300,\n # 'verbose': False\n # }\n\n # nnet.build(**fit_options)\n # results = nnet.evaluate(xtst, nsam = 100, msc = 2)\n # print(results)\n\n\nif __name__ == '__main__':\n main()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/docs/auto_examples/ex_nn.py b/docs/auto_examples/ex_nn.py new file mode 100644 index 0000000..db32364 --- /dev/null +++ b/docs/auto_examples/ex_nn.py @@ -0,0 +1,109 @@ +r""" +Residual Neural Network Construction +===================================== + +This example demonstrates how to use the Neural Network wrapper class, ``pytuq.surrogates.nn`` with scalar valued functions. +The constructor of the NN class accepts in an optional dictionary, ``net_options``, to specify additional hyperparameters. +The ``build()`` and ``evaluate()`` functions similarly accept dictionaries and explicit keyword arguments during their respective function calls. +""" + +import sys +import torch +import numpy as np + +from pytuq.surrogates.nn import NN +from quinn.solvers.nn_vi import NN_VI +from quinn.nns.rnet import RNet, Poly +from quinn.utils.plotting import myrc +from quinn.utils.maps import scale01ToDom +from quinn.func.funcs import Sine, Sine10, blundell + +def main(): + """Main function.""" + torch.set_default_dtype(torch.double) + myrc() + + ################################################################################# + ################################################################################# + + # defaults to cuda:0 if available + device_id='cuda:0' + device = torch.device(device_id if torch.cuda.is_available() else 'cpu') + print("Using device",device) + + nall = 15 # total number of points + trn_factor = 0.9 # which fraction of nall goes to training + ntst = 13 # separate test set + ndim = 1 # input dimensionality + datanoise = 0.02 # Noise in the generated data + true_model, nout = Sine, 1 # Scalar valued output example + + ################################################################################# + ################################################################################# + + # Domain: defining range of input variable + domain = np.tile(np.array([-np.pi, np.pi]), (ndim, 1)) + + np.random.seed(111) + + # Generating x-y training, validation, and testing data + xall = scale01ToDom(np.random.rand(nall, ndim), domain) + if true_model is not None: + yall = true_model(xall, datanoise=datanoise) + + if ntst > 0: + np.random.seed(100) + xtst = scale01ToDom(np.random.rand(ntst, ndim), domain) + if true_model is not None: + ytst = true_model(xtst, datanoise=datanoise) + + # (1) Initialize neural network with optional parameters for object instantiation + net_options = {'wp_function': Poly(0), + 'indim': ndim, + 'outdim': nout, + 'layer_pre': True, + 'layer_post': True, + 'biasorno': True, + 'nonlin': True, + 'mlp': False, + 'final_layer': None, + 'device': device, + } + + # Pass in unpacked net_options dict to constructor through kwargs + nnet = NN('RNet', 3, 3, **net_options) + + # (1.5) Split data into training and validation, assign through member functions + ntrn = int(trn_factor * nall) + indperm = range(nall) # np.random.permutation(range(nall)) + indtrn = indperm[:ntrn] + indval = indperm[ntrn:] + xtrn, xval = xall[indtrn, :], xall[indval, :] + ytrn, yval = yall[indtrn, :], yall[indval, :] + + nnet.set_validation_data(xval, yval) # optional + nnet.set_training_data(xtrn, ytrn) + + # (2) 1st build: Call build function (defaults to UQ method of variational inference) with optional parameters for fitting + nnet.build(datanoise=datanoise, lrate=0.01, batch_size=None, nsam=1, nepochs=5000, verbose=False) + + results = nnet.evaluate(xtst, nsam = 100, msc = 2) # Return samples of predictions with variance + covariance + print("Y_eval:", results['Y_eval'], end="\n\n") # Printing only samples of predictions + + # (3) 2nd build: Example of throwing an exception when passing in a network option through build() + # fit_options = {'datanoise': 0.05, + # 'outdim': 3, + # 'lrate': 0.01, + # 'batch_size': None, + # 'nsam': 1, + # 'nepochs': 300, + # 'verbose': False + # } + + # nnet.build(**fit_options) + # results = nnet.evaluate(xtst, nsam = 100, msc = 2) + # print(results) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/auto_examples/ex_nn.rst b/docs/auto_examples/ex_nn.rst new file mode 100644 index 0000000..2804ae2 --- /dev/null +++ b/docs/auto_examples/ex_nn.rst @@ -0,0 +1,157 @@ + +.. DO NOT EDIT. +.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. +.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: +.. "auto_examples/ex_nn.py" +.. LINE NUMBERS ARE GIVEN BELOW. + +.. only:: html + + .. note:: + :class: sphx-glr-download-link-note + + :ref:`Go to the end ` + to download the full example code. + +.. rst-class:: sphx-glr-example-title + +.. _sphx_glr_auto_examples_ex_nn.py: + + +Residual Neural Network Construction +===================================== + +This example demonstrates how to use the Neural Network wrapper class, ``pytuq.surrogates.nn`` with scalar valued functions. +The constructor of the NN class accepts in an optional dictionary, ``net_options``, to specify additional hyperparameters. +The ``build()`` and ``evaluate()`` functions similarly accept dictionaries and explicit keyword arguments during their respective function calls. + +.. GENERATED FROM PYTHON SOURCE LINES 9-109 + +.. code-block:: Python + + + import sys + import torch + import numpy as np + + from pytuq.surrogates.nn import NN + from quinn.solvers.nn_vi import NN_VI + from quinn.nns.rnet import RNet, Poly + from quinn.utils.plotting import myrc + from quinn.utils.maps import scale01ToDom + from quinn.func.funcs import Sine, Sine10, blundell + + def main(): + """Main function.""" + torch.set_default_dtype(torch.double) + myrc() + + ################################################################################# + ################################################################################# + + # defaults to cuda:0 if available + device_id='cuda:0' + device = torch.device(device_id if torch.cuda.is_available() else 'cpu') + print("Using device",device) + + nall = 15 # total number of points + trn_factor = 0.9 # which fraction of nall goes to training + ntst = 13 # separate test set + ndim = 1 # input dimensionality + datanoise = 0.02 # Noise in the generated data + true_model, nout = Sine, 1 # Scalar valued output example + + ################################################################################# + ################################################################################# + + # Domain: defining range of input variable + domain = np.tile(np.array([-np.pi, np.pi]), (ndim, 1)) + + np.random.seed(111) + + # Generating x-y training, validation, and testing data + xall = scale01ToDom(np.random.rand(nall, ndim), domain) + if true_model is not None: + yall = true_model(xall, datanoise=datanoise) + + if ntst > 0: + np.random.seed(100) + xtst = scale01ToDom(np.random.rand(ntst, ndim), domain) + if true_model is not None: + ytst = true_model(xtst, datanoise=datanoise) + + # (1) Initialize neural network with optional parameters for object instantiation + net_options = {'wp_function': Poly(0), + 'indim': ndim, + 'outdim': nout, + 'layer_pre': True, + 'layer_post': True, + 'biasorno': True, + 'nonlin': True, + 'mlp': False, + 'final_layer': None, + 'device': device, + } + + # Pass in unpacked net_options dict to constructor through kwargs + nnet = NN('RNet', 3, 3, **net_options) + + # (1.5) Split data into training and validation, assign through member functions + ntrn = int(trn_factor * nall) + indperm = range(nall) # np.random.permutation(range(nall)) + indtrn = indperm[:ntrn] + indval = indperm[ntrn:] + xtrn, xval = xall[indtrn, :], xall[indval, :] + ytrn, yval = yall[indtrn, :], yall[indval, :] + + nnet.set_validation_data(xval, yval) # optional + nnet.set_training_data(xtrn, ytrn) + + # (2) 1st build: Call build function (defaults to UQ method of variational inference) with optional parameters for fitting + nnet.build(datanoise=datanoise, lrate=0.01, batch_size=None, nsam=1, nepochs=5000, verbose=False) + + results = nnet.evaluate(xtst, nsam = 100, msc = 2) # Return samples of predictions with variance + covariance + print("Y_eval:", results['Y_eval'], end="\n\n") # Printing only samples of predictions + + # (3) 2nd build: Example of throwing an exception when passing in a network option through build() + # fit_options = {'datanoise': 0.05, + # 'outdim': 3, + # 'lrate': 0.01, + # 'batch_size': None, + # 'nsam': 1, + # 'nepochs': 300, + # 'verbose': False + # } + + # nnet.build(**fit_options) + # results = nnet.evaluate(xtst, nsam = 100, msc = 2) + # print(results) + + + if __name__ == '__main__': + main() + +.. _sphx_glr_download_auto_examples_ex_nn.py: + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-example + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download Jupyter notebook: ex_nn.ipynb ` + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download Python source code: ex_nn.py ` + + .. container:: sphx-glr-download sphx-glr-download-zip + + :download:`Download zipped: ex_nn.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/docs/auto_examples/ex_nn.zip b/docs/auto_examples/ex_nn.zip new file mode 100644 index 0000000..d3840b3 Binary files /dev/null and b/docs/auto_examples/ex_nn.zip differ diff --git a/docs/auto_examples/ex_pce.codeobj.json b/docs/auto_examples/ex_pce.codeobj.json new file mode 100644 index 0000000..dfe5350 --- /dev/null +++ b/docs/auto_examples/ex_pce.codeobj.json @@ -0,0 +1,65 @@ +{ + "PCE": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.surrogates.pce", + "module_short": "pytuq.surrogates.pce", + "name": "PCE" + } + ], + "fcb.sin4": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.utils.funcbank", + "module_short": "pytuq.utils.funcbank", + "name": "sin4" + } + ], + "np.array": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "array" + } + ], + "np.ones": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy", + "module_short": "numpy", + "name": "ones" + } + ], + "np.random.rand": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "rand" + } + ], + "np.random.seed": [ + { + "is_class": false, + "is_explicit": false, + "module": "numpy.random", + "module_short": "numpy.random", + "name": "seed" + } + ], + "scale01ToDom": [ + { + "is_class": false, + "is_explicit": false, + "module": "pytuq.utils.maps", + "module_short": "pytuq.utils.maps", + "name": "scale01ToDom" + } + ] +} \ No newline at end of file diff --git a/docs/auto_examples/ex_pce.ipynb b/docs/auto_examples/ex_pce.ipynb new file mode 100644 index 0000000..014222d --- /dev/null +++ b/docs/auto_examples/ex_pce.ipynb @@ -0,0 +1,79 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n# Polynomial Chaos Expansion Construction\n\nThis tutorial demonstrates how to create a Polynomial Chaos Expansion (PCE) surrogate model for scalar valued\nfunctions. We will use the ``pytuq.surrogates.pce`` wrapper class to approximate the function $\\sin^4(x)$.\n\nPyTUQ provides a number of utilities for your workflows, including modules for mapping, test functions, and metric comparisons --\nthese, along with others, can be found in the ``utils`` directory. You can also explore the ``func`` directory for additional sample functions to use.\n\nThis example below outlines how to:\n\n- Define basic parameters for the surrogate model\n- Use a sample function of $\\sin^4(x)$ from the ``utils`` directory\n- Set up a PCE surrogate model\n- And evaluate the performance of your model\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np\nimport pytuq.utils.funcbank as fcb\n\n\nfrom pytuq.surrogates.pce import PCE\nfrom pytuq.utils.maps import scale01ToDom\nfrom pytuq.lreg.anl import anl\nfrom pytuq.lreg.lreg import lsq\nfrom pytuq.utils.mindex import get_mi" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## #######################################################\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "N = 14 # Number of data points to generate\norder = 4 # Polynomial order\ntrue_model = fcb.sin4 # Function to approximate \ndim = 3 # Dimensionality of input data\n# dim = 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## #######################################################\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "np.random.seed(42) \n\n# Domain: defining range of input variable\ndomain = np.ones((dim, 2))*np.array([-1.,1.])\n\n# Generating x-y data\nx = scale01ToDom(np.random.rand(N,dim), domain)\ny = true_model(x)\n\n# Testing PCE class:\n\n# (1) Construct polynomial chaos expansion\npce = PCE(dim, order, 'LU')\n# pce = PCE(dim, order, 'HG')\n# pce = PCE(dim, order, ['LU', 'HG']) # dim = 2\n# pce = PCE(dim, order, ['HG', 'LU', 'HG', 'LU']) # dim = 4\n\npce.set_training_data(x, y)\n\n# (2) Pick method for linear regression object, defaulting to least squares regression\n# print(pce.build(regression = 'anl'))\n# print(pce.build(regression = 'anl', method = 'vi'))\nprint(pce.build())\n# print(pce.build(regression = 'lsq'))\n# print(pce.build(regression = 'bcs'))\n\n# (3) Make predictions for data points and print results:\nresults = pce.evaluate(x)\n\n# np.random.seed(45) # Generate a single random data point within the domain:\n# single_point = scale01ToDom(np.random.rand(1, dim), domain)\n# results = pce.evaluate(single_point)\n\n# Generate random input for evaluation:\n# x_eval = scale01ToDom(np.random.rand(10,dim), domain)\n# results = pce.evaluate(x_eval)\n\nprint(results)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/docs/auto_examples/ex_pce.py b/docs/auto_examples/ex_pce.py new file mode 100644 index 0000000..68672f8 --- /dev/null +++ b/docs/auto_examples/ex_pce.py @@ -0,0 +1,82 @@ +r""" +Polynomial Chaos Expansion Construction +======================================== + +This tutorial demonstrates how to create a Polynomial Chaos Expansion (PCE) surrogate model for scalar valued +functions. We will use the ``pytuq.surrogates.pce`` wrapper class to approximate the function :math:`\sin^4(x)`. + +PyTUQ provides a number of utilities for your workflows, including modules for mapping, test functions, and metric comparisons -- +these, along with others, can be found in the ``utils`` directory. You can also explore the ``func`` directory for additional sample functions to use. + +This example below outlines how to: + +- Define basic parameters for the surrogate model +- Use a sample function of :math:`\sin^4(x)` from the ``utils`` directory +- Set up a PCE surrogate model +- And evaluate the performance of your model + +""" + +import numpy as np +import pytuq.utils.funcbank as fcb + + +from pytuq.surrogates.pce import PCE +from pytuq.utils.maps import scale01ToDom +from pytuq.lreg.anl import anl +from pytuq.lreg.lreg import lsq +from pytuq.utils.mindex import get_mi + + +######################################################## +######################################################## +######################################################## + +N = 14 # Number of data points to generate +order = 4 # Polynomial order +true_model = fcb.sin4 # Function to approximate +dim = 3 # Dimensionality of input data +# dim = 1 + +######################################################## +######################################################## +######################################################## + +np.random.seed(42) + +# Domain: defining range of input variable +domain = np.ones((dim, 2))*np.array([-1.,1.]) + +# Generating x-y data +x = scale01ToDom(np.random.rand(N,dim), domain) +y = true_model(x) + +# Testing PCE class: + +# (1) Construct polynomial chaos expansion +pce = PCE(dim, order, 'LU') +# pce = PCE(dim, order, 'HG') +# pce = PCE(dim, order, ['LU', 'HG']) # dim = 2 +# pce = PCE(dim, order, ['HG', 'LU', 'HG', 'LU']) # dim = 4 + +pce.set_training_data(x, y) + +# (2) Pick method for linear regression object, defaulting to least squares regression +# print(pce.build(regression = 'anl')) +# print(pce.build(regression = 'anl', method = 'vi')) +print(pce.build()) +# print(pce.build(regression = 'lsq')) +# print(pce.build(regression = 'bcs')) + +# (3) Make predictions for data points and print results: +results = pce.evaluate(x) + +# np.random.seed(45) # Generate a single random data point within the domain: +# single_point = scale01ToDom(np.random.rand(1, dim), domain) +# results = pce.evaluate(single_point) + +# Generate random input for evaluation: +# x_eval = scale01ToDom(np.random.rand(10,dim), domain) +# results = pce.evaluate(x_eval) + +print(results) \ No newline at end of file diff --git a/docs/auto_examples/ex_pce.rst b/docs/auto_examples/ex_pce.rst new file mode 100644 index 0000000..02bf7b9 --- /dev/null +++ b/docs/auto_examples/ex_pce.rst @@ -0,0 +1,143 @@ + +.. DO NOT EDIT. +.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. +.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: +.. "auto_examples/ex_pce.py" +.. LINE NUMBERS ARE GIVEN BELOW. + +.. only:: html + + .. note:: + :class: sphx-glr-download-link-note + + :ref:`Go to the end ` + to download the full example code. + +.. rst-class:: sphx-glr-example-title + +.. _sphx_glr_auto_examples_ex_pce.py: + + +Polynomial Chaos Expansion Construction +======================================== + +This tutorial demonstrates how to create a Polynomial Chaos Expansion (PCE) surrogate model for scalar valued +functions. We will use the ``pytuq.surrogates.pce`` wrapper class to approximate the function :math:`\sin^4(x)`. + +PyTUQ provides a number of utilities for your workflows, including modules for mapping, test functions, and metric comparisons -- +these, along with others, can be found in the ``utils`` directory. You can also explore the ``func`` directory for additional sample functions to use. + +This example below outlines how to: + +- Define basic parameters for the surrogate model +- Use a sample function of :math:`\sin^4(x)` from the ``utils`` directory +- Set up a PCE surrogate model +- And evaluate the performance of your model + +.. GENERATED FROM PYTHON SOURCE LINES 19-31 + +.. code-block:: Python + + + import numpy as np + import pytuq.utils.funcbank as fcb + + + from pytuq.surrogates.pce import PCE + from pytuq.utils.maps import scale01ToDom + from pytuq.lreg.anl import anl + from pytuq.lreg.lreg import lsq + from pytuq.utils.mindex import get_mi + + + +.. GENERATED FROM PYTHON SOURCE LINES 32-34 + +####################################################### +####################################################### + +.. GENERATED FROM PYTHON SOURCE LINES 34-41 + +.. code-block:: Python + + + N = 14 # Number of data points to generate + order = 4 # Polynomial order + true_model = fcb.sin4 # Function to approximate + dim = 3 # Dimensionality of input data + # dim = 1 + + +.. GENERATED FROM PYTHON SOURCE LINES 42-44 + +####################################################### +####################################################### + +.. GENERATED FROM PYTHON SOURCE LINES 44-82 + +.. code-block:: Python + + + np.random.seed(42) + + # Domain: defining range of input variable + domain = np.ones((dim, 2))*np.array([-1.,1.]) + + # Generating x-y data + x = scale01ToDom(np.random.rand(N,dim), domain) + y = true_model(x) + + # Testing PCE class: + + # (1) Construct polynomial chaos expansion + pce = PCE(dim, order, 'LU') + # pce = PCE(dim, order, 'HG') + # pce = PCE(dim, order, ['LU', 'HG']) # dim = 2 + # pce = PCE(dim, order, ['HG', 'LU', 'HG', 'LU']) # dim = 4 + + pce.set_training_data(x, y) + + # (2) Pick method for linear regression object, defaulting to least squares regression + # print(pce.build(regression = 'anl')) + # print(pce.build(regression = 'anl', method = 'vi')) + print(pce.build()) + # print(pce.build(regression = 'lsq')) + # print(pce.build(regression = 'bcs')) + + # (3) Make predictions for data points and print results: + results = pce.evaluate(x) + + # np.random.seed(45) # Generate a single random data point within the domain: + # single_point = scale01ToDom(np.random.rand(1, dim), domain) + # results = pce.evaluate(single_point) + + # Generate random input for evaluation: + # x_eval = scale01ToDom(np.random.rand(10,dim), domain) + # results = pce.evaluate(x_eval) + + print(results) + +.. _sphx_glr_download_auto_examples_ex_pce.py: + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-example + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download Jupyter notebook: ex_pce.ipynb ` + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download Python source code: ex_pce.py ` + + .. container:: sphx-glr-download sphx-glr-download-zip + + :download:`Download zipped: ex_pce.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/docs/auto_examples/ex_pce.zip b/docs/auto_examples/ex_pce.zip new file mode 100644 index 0000000..1f6c5da Binary files /dev/null and b/docs/auto_examples/ex_pce.zip differ diff --git a/docs/auto_examples/images/sphx_glr_ex_genz_bcs_001.png b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_001.png new file mode 100644 index 0000000..99f7409 Binary files /dev/null and b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_001.png differ diff --git a/docs/auto_examples/images/sphx_glr_ex_genz_bcs_002.png b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_002.png new file mode 100644 index 0000000..a590523 Binary files /dev/null and b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_002.png differ diff --git a/docs/auto_examples/images/sphx_glr_ex_genz_bcs_003.png b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_003.png new file mode 100644 index 0000000..ce0031b Binary files /dev/null and b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_003.png differ diff --git a/docs/auto_examples/images/sphx_glr_ex_genz_bcs_004.png b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_004.png new file mode 100644 index 0000000..796bf33 Binary files /dev/null and b/docs/auto_examples/images/sphx_glr_ex_genz_bcs_004.png differ diff --git a/docs/auto_examples/images/thumb/sphx_glr_ex_genz_bcs_thumb.png b/docs/auto_examples/images/thumb/sphx_glr_ex_genz_bcs_thumb.png new file mode 100644 index 0000000..8f2d01b Binary files /dev/null and b/docs/auto_examples/images/thumb/sphx_glr_ex_genz_bcs_thumb.png differ diff --git a/docs/auto_examples/images/thumb/sphx_glr_ex_nn_thumb.png b/docs/auto_examples/images/thumb/sphx_glr_ex_nn_thumb.png new file mode 100644 index 0000000..8a5fed5 Binary files /dev/null and b/docs/auto_examples/images/thumb/sphx_glr_ex_nn_thumb.png differ diff --git a/docs/auto_examples/images/thumb/sphx_glr_ex_pce_thumb.png b/docs/auto_examples/images/thumb/sphx_glr_ex_pce_thumb.png new file mode 100644 index 0000000..8a5fed5 Binary files /dev/null and b/docs/auto_examples/images/thumb/sphx_glr_ex_pce_thumb.png differ diff --git a/docs/auto_examples/index.rst b/docs/auto_examples/index.rst new file mode 100644 index 0000000..92ded99 --- /dev/null +++ b/docs/auto_examples/index.rst @@ -0,0 +1,105 @@ +:orphan: + +Software Tutorials +=================== + +Welcome to the PyTUQ Tutorials page! Here, you'll find a series of comprehensive tutorials designed to +showcase the functionality of PyTUQ and guide you through its core features. + +Surrogates +~~~~~~~~~~ + +PyTUQ supports the construction of various surrogate models, including Polynomial Chaos Expansions for multivariate +random variables and Neural Networks, which can be accessed through the Quantification of Uncertainties in Neural Networks (QUiNN) library. + + +.. raw:: html + +
+ +.. thumbnail-parent-div-open + +.. raw:: html + +
+ +.. only:: html + + .. image:: /auto_examples/images/thumb/sphx_glr_ex_pce_thumb.png + :alt: + + :ref:`sphx_glr_auto_examples_ex_pce.py` + +.. raw:: html + +
Polynomial Chaos Expansion Construction
+
+ + +.. raw:: html + +
+ +.. only:: html + + .. image:: /auto_examples/images/thumb/sphx_glr_ex_nn_thumb.png + :alt: + + :ref:`sphx_glr_auto_examples_ex_nn.py` + +.. raw:: html + +
Residual Neural Network Construction
+
+ + +.. raw:: html + +
+ +.. only:: html + + .. image:: /auto_examples/images/thumb/sphx_glr_ex_genz_bcs_thumb.png + :alt: + + :ref:`sphx_glr_auto_examples_ex_genz_bcs.py` + +.. raw:: html + +
Function Approximation with Sparse Regression
+
+ + +.. thumbnail-parent-div-close + +.. raw:: html + +
+ + +.. toctree:: + :hidden: + + /auto_examples/ex_pce + /auto_examples/ex_nn + /auto_examples/ex_genz_bcs + + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-gallery + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download all examples in Python source code: auto_examples_python.zip ` + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download all examples in Jupyter notebooks: auto_examples_jupyter.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/docs/auto_examples/sg_execution_times.rst b/docs/auto_examples/sg_execution_times.rst new file mode 100644 index 0000000..d48caef --- /dev/null +++ b/docs/auto_examples/sg_execution_times.rst @@ -0,0 +1,43 @@ + +:orphan: + +.. _sphx_glr_auto_examples_sg_execution_times: + + +Computation times +================= +**00:05.760** total execution time for 3 files **from auto_examples**: + +.. container:: + + .. raw:: html + + + + + + + + .. list-table:: + :header-rows: 1 + :class: table table-striped sg-datatable + + * - Example + - Time + - Mem (MB) + * - :ref:`sphx_glr_auto_examples_ex_genz_bcs.py` (``ex_genz_bcs.py``) + - 00:05.760 + - 0.0 + * - :ref:`sphx_glr_auto_examples_ex_nn.py` (``ex_nn.py``) + - 00:00.000 + - 0.0 + * - :ref:`sphx_glr_auto_examples_ex_pce.py` (``ex_pce.py``) + - 00:00.000 + - 0.0 diff --git a/docs/conf.py b/docs/conf.py index 7d4651f..7132b2e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -55,7 +55,7 @@ 'IPython.sphinxext.ipython_console_highlighting', 'IPython.sphinxext.ipython_directive'] -# extensions += ['sphinx_gallery.gen_gallery'] +extensions += ['sphinx_gallery.gen_gallery'] numpydoc_show_class_members = False @@ -68,7 +68,7 @@ autoapi_type = "python" autoapi_template_dir = "_templates/autoapi" # Templates for AutoAPI documentation # autoapi_add_toctree_entry = False # Adding the generateed documentation into the TOC Tree -suppress_warnings = ["autoapi"] +suppress_warnings = ["autoapi", 'ref.python', 'toc.excluded'] autoapi_options = [ "members", "undoc-members", @@ -81,6 +81,15 @@ autoapi_keep_files = False # Keep the AutoAPI generated files on the filesystem # ----------------------------------------------------------------------------- +# Sphinx Gallery +# ----------------------------------------------------------------------------- +sphinx_gallery_conf = { + 'examples_dirs': '../examples/surrogates', # path to your example scripts. WIP: currently generating just surrogate examples + 'gallery_dirs': 'auto_examples', # path to where to save gallery generated output + 'filename_pattern': r'ex_genz_bcs\.py', # prefix for files that should be executed. Currently only 1 file being executed + 'ignore_pattern': r'ex_nn_json\.py', # exclude this specific file + # 'run_stale_examples': False, # default behavior is to only rebuild changed examples. Make True to force rebuild of examples (e.g. if underlying code changed) +} # Add any paths that contain templates here, relative to this directory. # templates_path = ['_templates'] @@ -96,6 +105,7 @@ # -- Options for HTML output ------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output +# https://sphinx-rtd-theme.readthedocs.io/en/stable/configuring.html# # html_theme = 'furo' html_theme = 'sphinx_rtd_theme' @@ -103,13 +113,15 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -# html_static_path = ['_static'] +html_static_path = ['_static'] +html_logo = 'pytuq logo.png' # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'sphinx' html_theme_options = { 'prev_next_buttons_location': 'bottom', + 'logo_only': True, 'style_external_links': False, 'vcs_pageview_mode': '', 'flyout_display': 'hidden', diff --git a/docs/index.rst b/docs/index.rst index 2c512f4..b83dc87 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -48,6 +48,12 @@ Check out the :ref:`Getting Started ` section for further information, in misc/installation misc/about +.. toctree:: + :maxdepth: 4 + :caption: Tutorials + + auto_examples/index + .. toctree:: :maxdepth: 4 :caption: API Reference diff --git a/docs/pytuq logo.png b/docs/pytuq logo.png new file mode 100644 index 0000000..e59864d Binary files /dev/null and b/docs/pytuq logo.png differ diff --git a/examples/surrogates/README.rst b/examples/surrogates/README.rst new file mode 100644 index 0000000..c638b71 --- /dev/null +++ b/examples/surrogates/README.rst @@ -0,0 +1,11 @@ +Software Tutorials +=================== + +Welcome to the PyTUQ Tutorials page! Here, you'll find a series of comprehensive tutorials designed to +showcase the functionality of PyTUQ and guide you through its core features. + +Surrogates +~~~~~~~~~~~ + +PyTUQ supports the construction of various surrogate models, including Polynomial Chaos Expansions for multivariate +random variables and Neural Networks, which can be accessed through the Quantification of Uncertainties in Neural Networks (QUiNN) library. \ No newline at end of file diff --git a/examples/surrogates/ex_genz_bcs.py b/examples/surrogates/ex_genz_bcs.py new file mode 100644 index 0000000..1dd7b14 --- /dev/null +++ b/examples/surrogates/ex_genz_bcs.py @@ -0,0 +1,420 @@ +r""" +Function Approximation with Sparse Regression +=============================================== + +In this tutorial, we demonstrate how to approximate a function with sparse regression by constructing a Polynomial Chaos (PC) surrogate with +Bayesian compressive sensing (BCS). The function we will approximate here is the Genz Oscillatory function, defined as: + +.. math:: + + f(x) = \cos\left(2 \pi s + \sum_{i=1}^d w_i x_i\right) + +Where: + +- :math:`s`: The shift parameter (``self.shift`` in the class). +- :math:`w_i`: The weights for each dimension (``self.weights`` in the class). +- :math:`x_i`: The input variables. +- :math:`d`: The dimensionality of the input :math:`x` (number of components in :math:`x`). + +Through 2 different build processes, we will construct two different PC surrogates to demonstrate effects of the BCS eta hyperparameter on model results. +The first build process will demonstrate most simply the construct-build-evaluate process when using BCS for our PC surrogate, along with a given a eta of 1e-10. +The second build process will select the most optimal eta for BCS through cross-validation algorithms (exposed here), which will soon be implemented in PyTUQ under-the-hood. +Afterwards, we will evaluate both models on testing and training data, returning parity plots and a Root Mean Square Error for each evaluation. +""" +# %% + +import os +import sys + +import numpy as np +import copy +import math +import pytuq.utils.funcbank as fcb +from matplotlib import pyplot as plt +from sklearn.metrics import root_mean_squared_error + +from pytuq.surrogates.pce import PCE +from pytuq.utils.maps import scale01ToDom +from pytuq.func.genz import GenzOscillatory + +# %% +# Setting a random number generator seed: + +# Random number generator +from scipy.stats import qmc +rng_seed = 43 + +# %% +# Constructing PC surrogate and generating data +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# After importing GenzOscillatory from ``pytuq.func.genz``, we generate the Genz function below, along with training data and testing data with output noise. +# This data and the corresponding Genz function will be used to create the same PC surrogate fitted in both examples: +# (1) The first surrogate will be fitted using BCS with a given eta, and (2) the second surrogate will be fitted using BCS with the most optimal eta. + +# Use Genz Oscillatory function in multiple dimensions for the true model +func_dim = 4 +func_weights = [1.0/(i+1)**2 for i in range(func_dim)] +func = GenzOscillatory(shift=0.25, weights=func_weights) +noise_std = 0.1 +rng = qmc.LatinHypercube(d=func_dim, seed=rng_seed) + +# Training data +np.random.seed(42) +n_trn = 70 +x_trn = rng.random(n=n_trn) # random numbers in [0,1]^d +y_trn = func(x_trn) + np.random.normal(0, noise_std, size = (n_trn,1)) + +# Testing data +n_tst = 10000 +x_tst = rng.random(n=n_tst) # random numbers in [0,1]^d +y_tst = func(x_tst) + np.random.normal(0, noise_std, size = (n_tst,1)) + +# %% +# With a stochastic dimensionality of 4 (defined above) and a polynomial order of 4, we construct the PC surrogate that will be used in both builds. +# You have the option of printing the PC surrogate's full basis, before BCS selects and retains the most significant PC coefficients to reduce the basis. + +# (1) Construct a PC surrogate +order = 4 +pce_surr = PCE(func_dim, order, 'LU', verbose = 1) + +# Optional verbosity output: +print("Full Basis and Current Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of Basis Terms:", len(pce_surr.pcrv.mindices[0])) + +# (1.5) Set training data +pce_surr.set_training_data(x_trn, y_trn[:,0]) + +# %% +# BCS with default settings (default eta) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# Here, we call the PCE class method of ``build()`` to build the linear regression model used to fit the surrogate. +# With the flag ``regression='bcs'``, we choose the BCS method for the fitting. A user-defined eta of 1e-10 is also passed in. + +# (2) Build the linear regression object for fitting +pce_surr.build(regression='bcs', eta=1.e-10) + +# Optional verbosity output: +print("Retained Basis and Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + +# %% +# After fitting, we evaluate the PCE using our training and testing data. To analyze the model's goodness of fit, +# we calculate the root mean square error between the surrogate results and the training and testing data. + +# (3) Evaluate the PC model +y_trn_approx = pce_surr.evaluate(x_trn) +y_tst_approx = pce_surr.evaluate(x_tst) + +# Evaluate goodness of fit with RMSE +rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) +print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + +rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) +print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + +# %% +# Notice above how the training RMSE error is almost half that of the testing RMSE error. This shows that our current model is overfitting, +# learning the training data with noise too well. To address this issue, we can explore selecting a better eta for the BCS fitting. + +# Plot the surrogate model's output vs. the training data output + +y_trn_mM = [y_trn[:,0].min(),y_trn[:,0].max()] + + +fig1 = plt.figure(figsize=(8,6)) +ax1 = fig1.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax1.plot(y_trn[:,0],y_trn_approx["Y_eval"],".") +ax1.plot(y_trn_mM,y_trn_mM) # Diagonal line + +ax1.set_xlabel("Train Data y", size=16) +ax1.set_ylabel("Predicted y", size=16); + +# %% + +# Plot the surrogate model's output vs. the testing data output + +y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + +fig2 = plt.figure(figsize=(8,6)) +ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") +ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + +ax2.set_xlabel("Test Data y", size=16) +ax2.set_ylabel("Predicted y", size=16); + +# %% +# BCS with optimal eta (found through cross-validation) +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# In this section, we use the same PC surrogate, ``pce_surr``, for the second build. We call the PCE class method of ``build()`` +# to build the linear regression model used to fit the surrogate. With the flag ``regression='bcs'``, we choose the BCS method for the fitting. +# +############################################################################## +# Instead of using a default eta, we call the cross-validation algorithm, ``optimize_eta()``, to choose the most optimal eta below. +# +# - With the flag ``plot=True``, the CV algorithm produces a graph of the training and testing (validation) data's RMSE values for each eta. The eta with the smallest RMSE for the validation data is the one chosen as the optimal eta. +# +############################################################################## +# Before that, we expose the cross-validation algorithm ``optimize_eta`` and its two helper functions, ``kfold_split`` and ``kfold_cv`` that will be called in this section. + +def kfold_split(nsamples,nfolds,seed=13): + """ + Return dictionary of training and testing pairs using k-fold cross-validation. + + Args: + nsamples (int): Total number of training samples. + nfolds (int): Number of folds to use for k-fold cross-validation. + seed (int, optional): Random seed for reproducibility. Defaults to 13. + + Returns: + dict: A dictionary where each key is the fold number (0 to nfolds-1) + and each value is a dictionary with: + - "train index" (np.ndarray): Indices of training samples. + - "val index" (np.ndarray): Indices of validation samples. + """ + # Returns split data where each data is one fold left out + KK=nfolds + rn = np.random.RandomState(seed) + + # Creating a random permutation of the samples indices list + indp=rn.permutation(nsamples) + + # Split the permuted indices into KK (or # folds) equal-sized chunks + split_index=np.array_split(indp,KK) + + # Dictionary to hold the indices of the training and validation samples + cvindices = {} + + # create testing and training folds + for j in range(KK): + # Iterating through the number of folds + fold = j + # Iterate through # folds, if i != fold number, + newindex = [split_index[i] for i in range(len(split_index)) if i != (fold)] + train_ind = np.array([],dtype='int64') + for i in range(len(newindex)): train_ind = np.concatenate((train_ind,newindex[i])) + test_ind = split_index[fold] + cvindices[j] = {'train index': train_ind, 'val index': test_ind} + + return cvindices + +# %% + +def kfold_cv(x,y,nfolds=3,seed=13): + """ + Splits data into training/testing pairs for kfold cross-val + x is a data matrix of size n x d1, d1 is dim of input + y is a data matrix of size n x d2, d2 is dim of output + + Args: + x (np.ndarray): + Input matrix with shape (n, d1) or 1D array with shape (n,). + Each row is a sample; columns are input features. + y (np.ndarray): + Target array with shape (n,) for single-output, or (n, d2) for + multi-output. If 1D, it is internally reshaped to (n, 1) before + slicing; outputs are `np.squeeze`d per fold. + nfolds (int, optional): Number of folds for cross-validation. Defaults to 3. + seed (int, optional): + Random seed for reproducible shuffling in `kfold_split`. Defaults to 13. + """ + if len(x.shape)>1: + n,d1 = x.shape + else: + n=x.shape + ynew = np.atleast_2d(y) + if len(ynew) == 1: ynew = ynew.T # change to shape (n,1) + _,d2 = ynew.shape + cv_idx = kfold_split(n,nfolds,seed) + + kfold_data = {} + for k in cv_idx.keys(): + kfold_data[k] = { + 'xtrain': x[cv_idx[k]['train index']], + 'xval': x[cv_idx[k]['val index']], + 'ytrain': np.squeeze(ynew[cv_idx[k]['train index']]), + 'yval': np.squeeze(ynew[cv_idx[k]['val index']]) + } # use squeeze to return 1d array + + # set train and test to the same if 1 fold + if nfolds == 1: + kfold_data[k]['xtrain'] = kfold_data[k]['xval'] + kfold_data[k]['ytrain'] = kfold_data[k]['yval'] + + return kfold_data + +# %% +def optimize_eta(pce, etas, nfolds, verbose=0, plot=False): + """ + Choose the optimum eta for Bayesian compressive sensing. Calculates the RMSE + for each eta for a specified number of folds. Selects the eta with the lowest + RMSE after averaging the RMSEs over the folds. + + Input: + y: 1D numpy array (vector) with function, evaluated at the + sample points [#samples,] + x: N-dimensional NumPy array with sample points [#samples, + #dimensions] + etas: NumPy array or list with the threshold for stopping the + algorithm. Smaller values retain more nonzero + coefficients + verbose: Flag for print statements + plot: Flag for whether to generate a plot for eta optimization + + Output: + eta_opt: Optimum eta + + """ + # Split data in k folds -> Get dictionary of data split in training + testing folds + kfold_data = kfold_cv(pce._x_train, pce._y_train, nfolds) + + # Each value has data for 1 fold. Each value is a list of the RMSEs for each possible eta in the fold. + RMSE_list_per_fold_tr = [] + + # Same but for testing data + RMSE_list_per_fold_test = [] + + # Make a copy of the PCE object to run the cross-validation algorithm on + pce_copy = copy.deepcopy(pce) + pce_copy.verbose = 0 + + # Loop through each fold + for i in range(nfolds): + + # Get the training and validation data + x_tr = kfold_data[i]['xtrain'] + y_tr = kfold_data[i]['ytrain'] + x_test = kfold_data[i]['xval'] + y_test = kfold_data[i]['yval'] + + # As we conduct BCS for this fold with each separate eta, the RMSEs will be added to these lists + RMSE_per_eta_tr = [] + RMSE_per_eta_test = [] + + # Set the x and y training data for the copied PCE object + pce_copy.set_training_data(x_tr, y_tr) + + # Loop through each eta + for eta in etas: + + # Conduct the BCS fitting. The object is automatically updated with new multiindex and coefficients received from the fitting. + cfs = pce_copy.build(regression = 'bcs', eta=eta) + + # Evaluate the PCE object at the training and validation points + y_tr_eval = (pce_copy.evaluate(x_tr))['Y_eval'] + y_test_eval = (pce_copy.evaluate(x_test))['Y_eval'] + + # Print statement for verbose flag + if verbose > 1: + print("Fold " + str(i + 1) + ", eta " + str(eta) + ", " + str(len(cfs)) + " terms retained out of a full basis of size " + str(len(pce.pcrv.mindices[0]))) + + # Calculate the RMSEs for the training and validation points. + # Append the values into the list of etas per fold. + MSE = np.square(np.subtract(y_tr, y_tr_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_tr.append(RMSE) + + MSE = np.square(np.subtract(y_test, y_test_eval)).mean() + RMSE = math.sqrt(MSE) + RMSE_per_eta_test.append(RMSE) + + # Now, append the fold's list of RMSEs for each eta into the list carrying the lists for all folds + RMSE_list_per_fold_tr.append(RMSE_per_eta_tr) + RMSE_list_per_fold_test.append(RMSE_per_eta_test) + + # After compiling the RMSE data for each eta from all the folds, we find the eta with the lowest validation RMSE to be our optimal eta. + # Compute the average and standard deviation of the training and testing RMSEs over the folds + avg_RMSE_tr = np.array(RMSE_list_per_fold_tr).mean(axis=0) + avg_RMSE_test = np.array(RMSE_list_per_fold_test).mean(axis=0) + + std_RMSE_tr = np.std(np.array(RMSE_list_per_fold_tr), axis=0) + std_RMSE_test = np.std(np.array(RMSE_list_per_fold_test), axis=0) + + # Choose the eta with lowest RMSE across all folds' testing data + eta_opt = etas[np.argmin(avg_RMSE_test)] + + # Plot RMSE vs. eta for training and testing RMSE + if plot: + + fig, ax = plt.subplots(figsize=(10,10)) + + plt.errorbar(etas, avg_RMSE_tr, xerr=None, yerr=std_RMSE_tr, linewidth=2, markersize=8, capsize=8, label=('Training')) + plt.errorbar(etas, avg_RMSE_test, xerr=None, yerr=std_RMSE_test, linewidth=2, markersize=8, capsize=8, label=('Validation')) + + plt.plot(eta_opt, np.min(avg_RMSE_test), marker="o", markersize=15, color='black', label=("Optimum")) + + plt.xlabel("Eta",fontsize=20) + plt.ylabel("RMSE",fontsize=20) + + # Change size of tick labels + plt.tick_params(axis='both', labelsize=16) + + plt.xscale('log') + plt.yscale('log') + + # Create legend + plt.legend(loc='upper left') + + # Save + plt.savefig('eta_opt.pdf', format='pdf', dpi=1200) + + return eta_opt + +# %% + +# We first create a list of possible etas to pass in: [1e-16, 1e-15, ... , 1e-2, 1e-1, 1] +etas = 1/np.power(10,[i for i in range(0,16)]) + +# Then, we call the function to choose the optimal eta: +eta_opt = optimize_eta(pce_surr, etas, 10, plot=True) + +# %% +# Now, with the optimal eta obtained, we run the fitting again. Then, we evaluate the PCE and produce a parity plot for the predicted output vs. the testing data. +# Notice that the larger eta, 10e-2, retained fewer basis terms (6) compared to the smaller user-defined eta of 10e-10 (which retained 20 basis terms). + +# Build the linear regression object for fitting +pce_surr.build(regression='bcs', eta=eta_opt) + +# Optional verbosity output: +print("Retained Basis and Coefficients:") +pce_surr.pcrv.printInfo() +print("Number of retained basis terms:", len(pce_surr.pcrv.mindices[0])) + +# %% +# Evaluate the PC model with training and testing data +y_trn_approx = pce_surr.evaluate(x_trn) +y_tst_approx = pce_surr.evaluate(x_tst) + +# Evaluate goodness of fit with RMSE +rmse_trn = root_mean_squared_error(y_trn[:,0],y_trn_approx["Y_eval"]) +print("The training RMSE error in the PCE BCS approximation is %.2e"%rmse_trn) + +rmse_tst = root_mean_squared_error(y_tst[:,0],y_tst_approx["Y_eval"]) +print("The testing RMSE error in the PCE BCS approximation is %.2e"%rmse_tst) + +# %% +# While the training RMSE error was almost half that of the testing RMSE error for the first fitting, the RMSE errors here are much closer to each other in value. +# This suggests that the model has more effectively generalized to the unseen data; a better eta has improved performance. + +# Plot the surrogate model's output vs. the testing data output +y_tst_mM = [y_tst[:,0].min(),y_tst[:,0].max()] + + +fig2 = plt.figure(figsize=(8,6)) +ax2 = fig2.add_axes([0.15, 0.15, 0.80, 0.75]) + + +ax2.plot(y_tst[:,0],y_tst_approx["Y_eval"],".") +ax2.plot(y_tst_mM,y_tst_mM) # Diagonal line + +ax2.set_xlabel("Test Data y", size=16) +ax2.set_ylabel("Predicted y", size=16); + +# sphinx_gallery_thumbnail_number = 2 \ No newline at end of file diff --git a/examples/surrogates/ex_nn.py b/examples/surrogates/ex_nn.py index 9d3a64e..db32364 100644 --- a/examples/surrogates/ex_nn.py +++ b/examples/surrogates/ex_nn.py @@ -1,8 +1,11 @@ -#!/usr/bin/env python - -"""This file is for testing the NN wrapper class with scalar valued functions. - Calling the NN constructor passes in an example net_options dictionary, while the build and evaluate - sections pass in both dictionaries and explicit keyword arguments during their respective function calls.""" +r""" +Residual Neural Network Construction +===================================== + +This example demonstrates how to use the Neural Network wrapper class, ``pytuq.surrogates.nn`` with scalar valued functions. +The constructor of the NN class accepts in an optional dictionary, ``net_options``, to specify additional hyperparameters. +The ``build()`` and ``evaluate()`` functions similarly accept dictionaries and explicit keyword arguments during their respective function calls. +""" import sys import torch diff --git a/examples/surrogates/ex_nn_json.py b/examples/surrogates/ex_nn_json.py index e4502d9..3fb2ca3 100644 --- a/examples/surrogates/ex_nn_json.py +++ b/examples/surrogates/ex_nn_json.py @@ -1,10 +1,13 @@ -#!/usr/bin/env python - -"""This file is for testing the NN wrapper class with scalar valued functions. - Calls the NN constructor with no optional keyword arguments before passing in an example json file with user-defined - fitting and building options. The build and evaluate sections pass in explicit keyword arguments - during their respective function calls to demonstrate updates to the neural network fitting/evaluating options. - When requested with a provided filename, the updated options are printed out to a json file.""" +r""" +Residual Neural Network Construction (JSON) +============================================= + +This file is for testing the NN wrapper class with scalar valued functions. +Calls the NN constructor with no optional keyword arguments before passing in an example json file with user-defined +fitting and building options. The build and evaluate sections pass in explicit keyword arguments +during their respective function calls to demonstrate updates to the neural network fitting/evaluating options. +When requested with a provided filename, the updated options are printed out to a json file. +""" import sys import torch diff --git a/examples/surrogates/ex_pce.py b/examples/surrogates/ex_pce.py index 72771bc..68672f8 100644 --- a/examples/surrogates/ex_pce.py +++ b/examples/surrogates/ex_pce.py @@ -1,6 +1,21 @@ -#!/usr/bin/env python +r""" +Polynomial Chaos Expansion Construction +======================================== -"""This file is for testing the PCE wrapper class with scalar valued functions.""" +This tutorial demonstrates how to create a Polynomial Chaos Expansion (PCE) surrogate model for scalar valued +functions. We will use the ``pytuq.surrogates.pce`` wrapper class to approximate the function :math:`\sin^4(x)`. + +PyTUQ provides a number of utilities for your workflows, including modules for mapping, test functions, and metric comparisons -- +these, along with others, can be found in the ``utils`` directory. You can also explore the ``func`` directory for additional sample functions to use. + +This example below outlines how to: + +- Define basic parameters for the surrogate model +- Use a sample function of :math:`\sin^4(x)` from the ``utils`` directory +- Set up a PCE surrogate model +- And evaluate the performance of your model + +""" import numpy as np import pytuq.utils.funcbank as fcb