diff --git a/ruff.toml b/ruff.toml index 5eca4398..ad3b4d6a 100644 --- a/ruff.toml +++ b/ruff.toml @@ -1,6 +1,12 @@ # ruff.toml line-length = 88 -exclude = ["simopt/gui/*", "notebooks/*"] +exclude = [ + "simopt/gui/*", + "notebooks/*", + "simopt/models/ermexample.py", + "workshop/erm_data_generator.ipynb", + "workshop/workshop.ipynb" +] [lint] select = [ diff --git a/simopt/models/ermexample.py b/simopt/models/ermexample.py new file mode 100644 index 00000000..4c533222 --- /dev/null +++ b/simopt/models/ermexample.py @@ -0,0 +1,203 @@ +"""Example problem of deterministic function with noise. + +Simulate a synthetic problem with a deterministic objective function +evaluated with noise. +""" + +from __future__ import annotations + +from typing import Annotated, ClassVar + +import numpy as np +from pydantic import BaseModel, Field + +from mrg32k3a.mrg32k3a import MRG32k3a +from simopt.base import ( + ConstraintType, + Model, + Objective, + Problem, + RepResult, + VariableType, +) +from simopt.input_models import InputModel + + +class ERMExampleModelConfig(BaseModel): + """Configuration model for ERMExample simulation. + + An empirical risk minimization model for linear regression. + """ + + beta: Annotated[ + tuple[float, ...], + Field( + default=(0.0, 0.0), + description="(intercept, slope) coefficients", + ), + ] + + +class ERMExampleProblemConfig(BaseModel): + """Configuration model for ERMExample Problem. + + Base class to implement simulation-optimization problems. + """ + + initial_solution: Annotated[ + tuple[float, ...], + Field( + default=(0.0, 0.0), + description="initial solution", + ), + ] + budget: Annotated[ + int, + Field( + default=1000, + description="max # of replications for a solver to take", + gt=0, + json_schema_extra={"isDatafarmable": False}, + ), + ] + + +class FileInputModel(InputModel): + def __init__(self, filename): + self.data = np.load(filename) + + def set_rng(self, rng: random.Random) -> None: + self.rng = rng + + def unset_rng(self) -> None: + self.rng = None + + def random(self) -> float: + n_rows = np.shape(self.data)[0] + resample_idx = np.random.choice(n_rows, size=1, replace=True) + resample_x = self.data[resample_idx, 0].item() + resample_y = self.data[resample_idx, 1].item() + return resample_x, resample_y + + +class ERMExampleModel(Model): + """A model that for the empirical risk of a linear regression model.""" + + class_name_abbr: ClassVar[str] = "ERMEXAMPLE" + class_name: ClassVar[str] = "Linear Regression ERM" + config_class: ClassVar[type[BaseModel]] = ERMExampleModelConfig + n_rngs: ClassVar[int] = 1 + n_responses: ClassVar[int] = 1 + + def __init__(self, fixed_factors: dict | None = None) -> None: + """Initialize the model. + + Args: + fixed_factors (dict | None): fixed factors of the model. + If None, use default values. + """ + # Let the base class handle default arguments. + super().__init__(fixed_factors) + self.resample_model = FileInputModel("workshop/erm_data.npy") + + def before_replicate(self, rng_list: list[MRG32k3a]) -> None: # noqa: D102 + self.resample_model.set_rng(rng_list[0]) + + def replicate(self) -> tuple[dict, dict]: + """Evaluate the squared error loss of a single observation. + + Returns: + tuple[dict, dict]: A tuple containing: + - responses (dict): Performance measures of interest, including: + - "sq_error_loss": Squared error loss of a single observation. + - gradients (dict): A dictionary of gradient estimates for + each response. + """ + beta0, beta1 = self.factors["beta"] + x, y = self.resample_model.random() + sq_error_loss = (y - beta0 - beta1 * x) ** 2 + error_loss = y - beta0 - beta1 * x + # gradients wrt beta0 and beta1 + grad_sq_error_loss = (-2 * error_loss, -2 * x * error_loss) + + # Compose responses and gradients. + responses = {"sq_error_loss": sq_error_loss} + gradients = {"sq_error_loss": {"beta": grad_sq_error_loss}} + return responses, gradients + + +class ERMExampleProblem(Problem): + """Base class to implement simulation-optimization problems.""" + + class_name_abbr: ClassVar[str] = "ERM-EXAMPLE-1" + class_name: ClassVar[str] = "Min Empirical Risk" + config_class: ClassVar[type[BaseModel]] = ERMExampleProblemConfig + model_class: ClassVar[type[Model]] = ERMExampleModel + n_objectives: ClassVar[int] = 1 + n_stochastic_constraints: ClassVar[int] = 0 + minmax: ClassVar[tuple[int, ...]] = (-1,) + constraint_type: ClassVar[ConstraintType] = ConstraintType.UNCONSTRAINED + variable_type: ClassVar[VariableType] = VariableType.CONTINUOUS + gradient_available: ClassVar[bool] = True + model_default_factors: ClassVar[dict] = {} + model_decision_factors: ClassVar[set[str]] = {"beta"} + + @property + def optimal_value(self) -> float | None: # noqa: D102 + # Compute optimal beta0 and beta1 + all_data = np.load("workshop/erm_data.npy") + x = all_data[:, 0] + y = all_data[:, 1] + optbeta1, optbeta0 = np.polyfit(x, y, 1) + opttrainingmse = np.mean( + [(yy - optbeta0 - optbeta1 * xx) ** 2 for (xx, yy) in zip(x, y)] + ) + return opttrainingmse + + @property + def optimal_solution(self) -> tuple | None: # noqa: D102 + # Compute optimal beta0 and beta1 + all_data = np.load("workshop/erm_data.npy") + x = all_data[:, 0] + y = all_data[:, 1] + optbeta1, optbeta0 = np.polyfit(x, y, 1) + return (optbeta0, optbeta1) + + @property + def dim(self) -> int: # noqa: D102 + return 2 + + @property + def lower_bounds(self) -> tuple: # noqa: D102 + return (-np.inf,) * self.dim + + @property + def upper_bounds(self) -> tuple: # noqa: D102 + return (np.inf,) * self.dim + + def vector_to_factor_dict(self, vector: tuple) -> dict: # noqa: D102 + return {"beta": vector[:]} + + def factor_dict_to_vector(self, factor_dict: dict) -> tuple: # noqa: D102 + return tuple(factor_dict["beta"]) + + def replicate(self, _x: tuple) -> RepResult: # noqa: D102 + responses, gradients = self.model.replicate() + objectives = [ + Objective( + stochastic=responses["sq_error_loss"], + stochastic_gradients=gradients["sq_error_loss"]["beta"], + ) + ] + return RepResult(objectives=objectives) + + def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: # noqa: D102 + # beta = tuple([rand_sol_rng.uniform(-2, 2) for _ in range(self.dim)]) + beta = tuple( + rand_sol_rng.mvnormalvariate( + mean_vec=[1.0] * self.dim, + cov=np.eye(self.dim).tolist(), + factorized=False, + ) + ) + return beta diff --git a/ty.toml b/ty.toml index b3261959..a3c39d0c 100644 --- a/ty.toml +++ b/ty.toml @@ -19,6 +19,7 @@ exclude = [ "simopt/gui/new_experiment_window.py", "simopt/gui/plot_window.py", "simopt/model.py", + "simopt/models/ermexample.py", "simopt/plot_type.py", "simopt/plots", "simopt/problem.py", diff --git a/workshop/erm_data.npy b/workshop/erm_data.npy new file mode 100644 index 00000000..c117ae3d Binary files /dev/null and b/workshop/erm_data.npy differ diff --git a/workshop/erm_data_generator.ipynb b/workshop/erm_data_generator.ipynb new file mode 100644 index 00000000..fc805e17 --- /dev/null +++ b/workshop/erm_data_generator.ipynb @@ -0,0 +1,43 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "5c2c7b4c", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "n = 10000 # number of observations\n", + "# X ~ Normal(0, 1)\n", + "# Y|X ~ Normal(1 + x, 0.1)\n", + "x = np.random.normal(size=10000)\n", + "y = 1 + x + np.random.normal(loc=0, scale=0.1, size=10000)\n", + "xydata = np.array([[xx, yy] for (xx, yy) in zip(x, y)])\n", + "np.save(\"erm_data.npy\", xydata)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "simopt", + "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.13.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/workshop/workshop.ipynb b/workshop/workshop.ipynb index 2c0527eb..29b52f8b 100644 --- a/workshop/workshop.ipynb +++ b/workshop/workshop.ipynb @@ -24,7 +24,8 @@ "\n", "import sys\n", "\n", - "sys.path.append(\"venv\\\\lib\\\\site-packages\")" + "sys.path.append(\"venv\\\\lib\\\\site-packages\")\n", + "sys.path.append(\"..\")" ] }, { @@ -39,8 +40,8 @@ "import simopt.experiment_base as expbase\n", "from simopt.experiment_base import PlotType\n", "\n", - "# Import Example problem and Random Search and ADAM solvers.\n", - "from simopt.models.example import ExampleProblem\n", + "# Import ERM-Example problem and Random Search and ADAM solvers.\n", + "from simopt.models.ermexample import ERMExampleProblem\n", "from simopt.solvers.adam import ADAM\n", "from simopt.solvers.randomsearch import RandomSearch" ] @@ -51,11 +52,15 @@ "source": [ "In this portion of the workshop, we'll be working with a problem and two solvers.\n", "\n", - "**Problem:** Minimize $||x||^2$ with additive Gaussian noise over $x \\in \\mathbb{R}^2$.\n", + "**Problem:** Find the slope ($\\beta_1$) and intercept ($\\beta_0$) coefficients that minimize the training MSE of a simple linear regression model, i.e., least squares regression:\n", + "\n", + "$$ \\min_{\\beta_0, \\beta_1} \\frac{1}{n}\\sum_{i=1}^n (y_i - (\\beta_0 + \\beta_1 x_i))^2, $$\n", + "\n", + "where $Y(x) = 1 + x + \\epsilon$ and $\\epsilon \\sim N(0, 0.1)$.\n", "\n", "**Solver:** Random Search\n", - "* Randomly samples solutions. For this two-dimensional problem, solutions are sampled from a MVN distribution with mean vector (0, 0) and variance-covariance matrix (1, 0; 0, 1).\n", - "* Takes a fixed number of observations (replications) at each solution.\n", + "* Randomly samples ($\\beta_0$, $\\beta_1$) from a bivariate normal distribution with mean vector (1, 1) and variance-covariance matrix (1, 0; 0, 1).\n", + "* Draws a fixed number of observations from the training dataset (i.e., a minibatch) and computes the squared error loss.\n", "* [Full documentation](https://simopt.readthedocs.io/en/latest/randomsearch.html)\n", "\n", "**Solver:** ADAM\n", @@ -73,11 +78,11 @@ "# CODE CELL [1]\n", "\n", "# Instantiate the problem and the Random Search solver, with specifications.\n", - "my_problem = ExampleProblem(\n", - " fixed_factors={\"initial_solution\": (2.0, 2.0), \"budget\": 200}\n", + "my_problem = ERMExampleProblem(\n", + " fixed_factors={\"initial_solution\": (0.0, 0.0), \"budget\": 2000}\n", ")\n", "my_rand_search_solver = RandomSearch(\n", - " fixed_factors={\"crn_across_solns\": True, \"sample_size\": 10}\n", + " fixed_factors={\"crn_across_solns\": True, \"sample_size\": 100}\n", ")\n", "\n", "# Pair the problem and solver for experimentation.\n", @@ -88,7 +93,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Let's see how Random Search does on this toy problem." + "Let's see how Random Search does on this problem." ] }, { @@ -99,20 +104,16 @@ "source": [ "# CODE CELL [2]\n", "\n", - "# Run 10 macroreplications of Random Search on the Example Problem.\n", + "# Run 10 macroreplications of Random Search on the ERM-Example Problem.\n", "myexperiment.run(n_macroreps=10)\n", "\n", "# Post-process the results.\n", "myexperiment.post_replicate(n_postreps=200)\n", "expbase.post_normalize(experiments=[myexperiment], n_postreps_init_opt=200)\n", "\n", - "# [Results are saved in a file called experiments/outputs/RNDSRCH_on_EXAMPLE-1.pickle.]\n", + "# [Results are saved in a file called experiments//outputs/RNDSRCH_on_ERM-EXAMPLE-1.pickle.]\n", "# [The file is not human-readable, so we'll skip looking at it.]\n", "\n", - "# Record a summary of the results in a human-readable way.\n", - "myexperiment.log_experiment_results()\n", - "# [Go check out the file called experiments/logs/RNDSRCH_on_EXAMPLE-1_experiment_results.txt] # noqa: E501\n", - "\n", "# Plot the (unnormalized) progress curves from the 10 macroreplications.\n", "expbase.plot_progress_curves(\n", " experiments=[myexperiment], plot_type=PlotType.ALL, normalize=False\n", @@ -134,7 +135,7 @@ "\n", "In CODE CELL [1], play around with the arguments when initializing `myproblem` and `mysolver`.\n", "\n", - "Vary factors of the Example problem:\n", + "Vary factors of the ERM-Example problem:\n", "- Change the initial solution.\n", "- Change the budget, i.e., the max number of replications. \n", "\n", @@ -153,14 +154,14 @@ "\n", "### Exercise \\#2\n", "\n", - "1. Open the file simopt/model/example.py in the VS Code editor.\n", - "2. Let's change how random search randomly samples solutions in R^2. For starters, uncomment Line 430\n", + "1. Open the file simopt/model/ermexample.py in the VS Code editor.\n", + "2. Let's change how random search randomly samples solutions in $\\mathbb{R}^2$. For starters, uncomment Line 193\n", "\n", - " `x = tuple([rand_sol_rng.uniform(-2, 2) for _ in range(self.dim)])`\n", + " `beta = tuple([rand_sol_rng.uniform(-2, 2) for _ in range(self.dim)])`\n", "\n", - " and comment out Lines 431-437\n", + " and comment out Lines 194-200\n", " \n", - " `x = tuple(rand_sol_rng.mvnormalvariate(mean_vec=np.zeros(self.dim), cov=np.eye(self.dim), factorized=False))`\n", + " `beta = tuple(rand_sol_rng.mvnormalvariate(mean_vec=[1.0] * self.dim, cov=np.eye(self.dim), factorized=False))`\n", "\n", "3. Restart the kernel using the Restart Button at the top of this notebook. This will ensure the new version of the source code is being imported.\n", "4. Run COMBO CODE CELL [0 + 1 + 2] below (this effectively reruns CODE CELLS [0], [1], and [2]). *How have the plots changed?*\n", @@ -181,17 +182,18 @@ "import sys\n", "\n", "sys.path.append(\"venv\\\\lib\\\\site-packages\")\n", + "sys.path.append(\"..\")\n", "import simopt.experiment_base as expbase\n", "from simopt.experiment_base import PlotType\n", - "from simopt.models.example import ExampleProblem\n", + "from simopt.models.ermexample import ERMExampleProblem\n", "from simopt.solvers.adam import ADAM\n", "from simopt.solvers.randomsearch import RandomSearch\n", "\n", - "my_problem = ExampleProblem(\n", - " fixed_factors={\"initial_solution\": (2.0, 2.0), \"budget\": 200}\n", + "my_problem = ERMExampleProblem(\n", + " fixed_factors={\"initial_solution\": (0.0, 0.0), \"budget\": 2000}\n", ")\n", "my_rand_search_solver = RandomSearch(\n", - " fixed_factors={\"crn_across_solns\": True, \"sample_size\": 10}\n", + " fixed_factors={\"crn_across_solns\": True, \"sample_size\": 100}\n", ")\n", "myexperiment = expbase.ProblemSolver(problem=my_problem, solver=my_rand_search_solver)\n", "\n", @@ -222,8 +224,8 @@ "source": [ "# CODE CELL [3]\n", "\n", - "my_adam_solver = ADAM(fixed_factors={\"crn_across_solns\": True, \"r\": 10})\n", - "# Create a grouping of Example-RandomSearch and Example-ADAM pairs.\n", + "my_adam_solver = ADAM(fixed_factors={\"crn_across_solns\": True, \"r\": 100})\n", + "# Create a grouping of ERM-Example-RandomSearch and ERM-Example-ADAM pairs.\n", "mygroupexperiment = expbase.ProblemsSolvers(\n", " problems=[my_problem], solvers=[my_rand_search_solver, my_adam_solver]\n", ")\n", @@ -236,7 +238,7 @@ "# Record a summary of the results in a human-readable way.\n", "mygroupexperiment.log_group_experiment_results()\n", "# [Go check out the file called\n", - "# experiments/logs/group_RNDSRCH_ADAM_on_EXAMPLE-1_group_experiment_results.txt]\n", + "# experiments//logs/group_RNDSRCH_ADAM_on_ERM-EXAMPLE-1_group_experiment_results.txt]\n", "\n", "# Plot the mean progress curve for each solver from the 10 macroreplications.\n", "expbase.plot_progress_curves(\n", @@ -249,104 +251,11 @@ ")\n", "# [The plot should be displayed in the output produced below.]" ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "#### Your turn.\n", - "\n", - "### Exercise \\#3\n", - "\n", - "1. Open simopt/model/example.py again.\n", - "2. Change the noise in the objective function evaluations to create a slightly different 2D optimization problem. This can be done by changing Line 99: \n", - " \n", - " `fn_eval_at_x = np.linalg.norm(x) ** 2 + noise_rng.normalvariate()`\n", - "\n", - " where `x` is a numpy array of length two. For starters, try passing the argument `sigma=10` into the function call `noise_rng.normalvariate()`. The default value is `sigma=1`, so this has the effect of increasing the common variance of the noise from 1 to 100.\n", - "3. Restart the kernel and run COMBO CODE CELL [0 + 1 + 3] below. *How have the plots changed? Why haven't they changed more?*\n", - "\n", - "4. Next, change the underlying objective function by replacing `np.linalg.norm(x) ** 2` in Line 99 with some other two-dimensional function of `x`, e.g., `1 - np.exp(-np.linalg.norm(x) ** 2)`. (This objective function looks like an upside-down standard bivariate normal pdf, rescaled.)\n", - "5. Depending of your choice of new objective function, you MAY need to change other parts of the code, including:\n", - " * The gradient of `f(x)` in Line 103. For the example given above, this would need to be changed from\n", - " \n", - " `gradients = {\"est_f(x)\": {\"x\": tuple(2 * x)}}`\n", - "\n", - " to\n", - "\n", - " `gradients = {\"est_f(x)\": {\"x\": tuple(2 * x * np.exp(-np.linalg.norm(x) ** 2))}}`\n", - " * If you change the problem to a maxmization problem, you will need to change Line 190 from\n", - " \n", - " `return (-1,)`\n", - " \n", - " to\n", - " \n", - " `return (1,)`.\n", - " * The optimal solution in Line 214. (For the running example, this will not be necessary.)\n", - " * The optimal objective function value in Line 208. (For the running example, this will not be necessary.)\n", - "6. Restart the kernel and run COMBO CODE CELL [0 + 1 + 3] below. *How have the plots changed?*\n", - "\n", - "**Extra for Experts:** Change the dimension of the problem. To do this, you will need to change the dimension of the default initial solution, defined in Line 185." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# COMBO CODE CELL [0 + 1 + 3]\n", - "import os\n", - "\n", - "os.chdir(\"../\")\n", - "import sys\n", - "\n", - "sys.path.append(\"venv\\\\lib\\\\site-packages\")\n", - "import simopt.experiment_base as expbase\n", - "from simopt.experiment_base import PlotType\n", - "from simopt.models.example import ExampleProblem\n", - "from simopt.solvers.adam import ADAM\n", - "from simopt.solvers.randomsearch import RandomSearch\n", - "\n", - "my_problem = ExampleProblem(\n", - " fixed_factors={\"initial_solution\": (2.0, 2.0), \"budget\": 200}\n", - ")\n", - "my_rand_search_solver = RandomSearch(\n", - " fixed_factors={\"crn_across_solns\": True, \"sample_size\": 10}\n", - ")\n", - "my_adam_solver = ADAM(fixed_factors={\"crn_across_solns\": True, \"r\": 10})\n", - "\n", - "mygroupexperiment = expbase.ProblemsSolvers(\n", - " problems=[my_problem], solvers=[my_rand_search_solver, my_adam_solver]\n", - ")\n", - "mygroupexperiment.run(n_macroreps=10)\n", - "mygroupexperiment.post_replicate(n_postreps=200)\n", - "mygroupexperiment.post_normalize(n_postreps_init_opt=200)\n", - "mygroupexperiment.log_group_experiment_results()\n", - "expbase.plot_progress_curves(\n", - " experiments=[\n", - " mygroupexperiment.experiments[0][0],\n", - " mygroupexperiment.experiments[1][0],\n", - " ],\n", - " plot_type=PlotType.MEAN,\n", - " normalize=False,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Other demonstrations, time permitting\n", - "* Walkthrough `replicate()` method of simopt/models/ironore.py to illustrate what the code looks like for a typical stochastic simulation model.\n", - "* Walkthrough `solve()` method of simopt/solvers/ADAM.py to illustrate what the code looks like for a typical solver." - ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "simopt", "language": "python", "name": "python3" }, @@ -360,7 +269,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.13.10" } }, "nbformat": 4,