PyMOSO

PyMOSO is software for solving multi-objective simulation optimization (MOSO) problems and for creating, comparing, and testing MOSO algorithms.

Reference

If you use PyMOSO in work leading to publication, please cite the paper which introduces PyMOSO.

Cooper, K., Hunter, S. R. 2020. PyMOSO: Software for Multi-Objective Simulation Optimization with R-PERLE and R-MinRLE. INFORMS Journal on Computing, 32(4): 1101-1108. DOI: http://dx.doi.org/10.1287/ijoc.2019.0902, Preprint: http://web.ics.purdue.edu/~hunter63/PAPERS/pre2019coohun.pdf.

Additional Reading

The initial release of PyMOSO contains solvers that implement four total algorithms, in alphabetical order: R-MinRLE, R-PE, R-PERLE, and R-SPLINE. The algorithms R-MinRLE, R-PE, R-PERLE were introduced in the following paper:

Cooper, K., Hunter, S. R., and Nagaraj, K. 2020. Bi-objective simulation optimization on integer lattices using the epsilon-constraint method in a retrospective approximation framework. INFORMS Journal on Computing, 32(4): 1080-1100. DOI: http://dx.doi.org/10.1287/ijoc.2019.0918, Preprint: http://web.ics.purdue.edu/~hunter63/PAPERS/pre2019coohunnag.pdf.

The algorithm R-SPLINE was introduced in the following paper:

Wang, H., Pasupathy, R., and Schmeiser, B. W. 2013. Integer-ordered simulation optimization using R-SPLINE: Retrospective Search with Piecewise-Linear Interpolation and Neighborhood Enumeration. ACM Transactions on Modeling and Computer Simulation, Vol. 23, No. 3, Article 17 (July 2013), 24 pages. http://dx.doi.org/10.1145/2499913.2499916

We recommend reading these papers to understand the algorithms, what they return, and the algorithm parameter options that we describe in the user manual.

Table of Contents

• PyMOSO
  • Reference
  • Additional Reading
  • Installation
    • Install PyMOSO from the Python Packaging Index using pip
    • Install PyMOSO from the repository using pip
    • Install PyMOSO from Source Code
  • Command Line Interface (CLI)
    • CLI help
    • The listitems command for viewing solvers, testers, and oracles included in PyMOSO
    • The solve command
      • The Example Oracle
      • Table of Algorithm-Specific Parameters
    • The testsolve Command
      • The Example Tester
  • Implementing problems, testers, and algorithms in PyMOSO
    • Implementing PyMOSO Oracles
      • Example Oracle that Wraps a C Simulation
      • Example Wrapper with PyMOSO Random Numbers
    • Implementing PyMOSO Testers
      • Example Metric 1
      • Example Metric 2
    • Implementing PyMOSO Algorithms
      • Template Accelerator
      • Template RA Solver
      • Template MOSO Solver
        • Take Simulation Replications at a Point
        • Find Neighbors and Take Simulation Replications
        • Argsort a Dictionary of Points
        • Select the Minimizer and its Value
        • Use SPLINE to Retrive a Local Minimizer
        • Find the Non-Dominated Points in a Dictionary
        • Randomly Select Points in a Set
    • Using solve and testsolve in Python Programs
      • Minimal solve Example
      • Some solve Examples with Options
      • A testsolve Example
      • Computing a Metric on testsolve Output
  • PyMOSO Object Reference
    • The pymoso.prng.mrg32k3a Module
    • The pymoso.chnbase Module
    • The pymoso.chnutils Module
    • The Oracle Class
    • The MOSOSolver Class
    • The RASolver Class
    • The RLESolver Class

Installation

Since PyMOSO is programmed in Python, every PyMOSO user must first install Python, which can be downloaded from https://www.python.org/downloads/. PyMOSO is now compatible with Python versions 3.10 and higher. Legacy versions (1.0.8 and lower) are compatible with Python 3.6. In the remainder of this section, we assume an appropriate Python version is installed. We discuss three different methods to install PyMOSO: first, from the Python Packaging Index; second, directly from our source code using git; and third, installing PyMOSO from our source code.

Before you begin

PyMOSO requires Python 3.10 or newer, which can be downloaded from https://www.python.org/downloads/. To check which version you have, open a terminal and run python --version. On Windows, use py --version; on some systems, use python3 --version.

We recommend installing PyMOSO into a virtual environment, a self-contained folder that keeps PyMOSO and its setup separate from other Python software on your computer. Many recent operating systems require this and will refuse to install packages without one. To create and activate a virtual environment, open a terminal, navigate to the folder where you want to work, and run:

On macOS and Linux:

python3 -m venv .venv
source .venv/bin/activate

On Windows:

py -m venv .venv
.venv\Scripts\activate

While the environment is active, the commands python and pip refer to the environment's own copy, on every operating system. Run the activation command again whenever you open a new terminal to use PyMOSO.

Install PyMOSO from the Python Packaging Index

We keep stable releases of PyMOSO on the Python Packaging Index (PyPI), and we recommend this method for most users. With your virtual environment active, run:

python -m pip install pymoso

We write python -m pip rather than pip because it guarantees that PyMOSO is installed for the same Python you will run it with.

Install PyMOSO from the repository

Users with git installed can install the most recent development version of PyMOSO directly from our source code:

python -m pip install "git+https://github.com/pymoso/PyMOSO.git"

We consider the latest source to be less stable than the releases we publish to PyPI, and thus we recommend most users install PyMOSO from PyPI.

Install PyMOSO from source code

Users may follow the steps below to install PyMOSO from any version of the source code.

1. Acquire the PyMOSO source code from https://github.com/pymoso/PyMOSO. Users without git can use the green Code button and choose Download ZIP, then unzip the file. Specific versions are available from the repository's Releases page.
2. Open a terminal and navigate into the project directory, the folder that contains the file setup.py.
3. With your virtual environment active, install PyMOSO from the current directory:

python -m pip install .

The trailing . means "this folder."

To confirm any of the three methods worked, run pymoso --version.

Command Line Interface (CLI)

PyMOSO users solving MOSO problems and testing MOSO algorithms may do so using the command line interface. First, we show how to access the included help file. Then, we show how to view the lists of solvers, testers, and oracles installed by default with PyMOSO. Finally, we discuss the solve and testsolve commands.

CLI help

PyMOSO includes a command line help file. The help file shows syntax templates for every PyMOSO command, the available options, and a selection of example invocations. The pymoso --help invocation prints the file to the terminal. The file is also printed when PyMOSO cannot parse an invocation that begins with pymoso. We show the current help file below.

$ pymoso --help
usage: pymoso [-h] [-v] <command> ...

pymoso

positional arguments:
  <command>
    listitems    List the built-in solvers, problems, and testers.
    solve        Solve a MOSO problem.
    testsolve    Test a MOSO algorithm on a MOSO problem.

options:
  -h, --help     show this help message and exit
  -v, --version  show program's version number and exit

Examples:
  pymoso listitems
  pymoso solve ProbTPA RPERLE 4 14
  pymoso solve --budget=100000 --odir=test1  ProbTPB RMINRLE 3 12
  pymoso solve --seed 12345 32123 5322 2 9543 666666666 ProbTPC RPERLE 5 5 5
  pymoso solve --simpar=4 --param betaeps 0.4 ProbTPA RPERLE 30 30
  pymoso solve --param radius 3 ProbTPA RPERLE 45 45
  pymoso testsolve --isp=16 --proc=4 TPATester RPERLE
  pymoso testsolve --isp=20 --proc=10 --metric --crn TPBTester RMINRLE 9 9

Help:
  Use the listitems command to view a list of available solvers, problems, and
  test problems.

$ pymoso solve --help
usage: pymoso solve [-h] [--budget BUDGET] [--odir ODIR] [--crn]
                    [--seed <s> <s> <s> <s> <s> <s>] [--param <param> <val>]
                    [--simpar SIMPAR]
                    <problem> <solver> <x> [<x> ...]

positional arguments:
  <problem>
  <solver>
  <x>

options:
  -h, --help            show this help message and exit
  --budget BUDGET       Set the simulation budget [default: 200]
  --odir ODIR           Set the output file directory name. [default: testrun]
  --crn                 Set if common random numbers are desired.
  --seed <s> <s> <s> <s> <s> <s>
                        Set the random number seed with 6 spaced integers.
  --param <param> <val>
                        Specify a solver-specific parameter <param> <val>.
                        Repeatable.
  --simpar SIMPAR       Set number of parallel processes for simulation
                        replications. [default: 1]

$ pymoso testsolve --help
usage: pymoso testsolve [-h] [--budget BUDGET] [--odir ODIR] [--crn]
                        [--seed <s> <s> <s> <s> <s> <s>]
                        [--param <param> <val>] [--isp ISP] [--proc PROC]
                        [--metric]
                        <tester> <solver> [<x> ...]

positional arguments:
  <tester>
  <solver>
  <x>

options:
  -h, --help            show this help message and exit
  --budget BUDGET       Set the simulation budget [default: 200]
  --odir ODIR           Set the output file directory name. [default: testrun]
  --crn                 Set if common random numbers are desired.
  --seed <s> <s> <s> <s> <s> <s>
                        Set the random number seed with 6 spaced integers.
  --param <param> <val>
                        Specify a solver-specific parameter <param> <val>.
                        Repeatable.
  --isp ISP             Set number of algorithm instances to solve. [default:
                        1]
  --proc PROC           Set number of parallel processes for the algorithm
                        instances. [default: 1]
  --metric              Set if metric computation is desired.

$ pymoso listitems --help
usage: pymoso listitems [-h]

options:
  -h, --help  show this help message and exit

For now, PyMOSO has three commands: listitems, solve, and testsolve, which we explain below.

The listitems command for viewing solvers, testers, and oracles included in PyMOSO

The default installation of PyMOSO includes a selection of solvers, testers, and oracles. Users can view the complete lists of included solvers, testers, and oracles using the pymoso listitems command. We show the current listing below. Test problems A, B, and C refer to those in Cooper et al (2018).

Solver                 Description
*********************  ************************************************************
RMINRLE                A solver using R-MinRLE for integer-ordered MOSO.
RPE                    R-Pe solver for bi-objective simulation optimization.
RPERLE                 R-PERLE solver for bi-objective simulation optimization.
RSPLINE                R-SPLINE solver for single-objective simulation
                       optimization.

Problems               Description                                                   Test Name (if available)
*********************  ************************************************************  ************************
BSProb                 An Oracle that simulates a bus-scheduling problem,            BSTester
                       minimizing total passenger wait time and bus-dispatch cost.
ProbExpensive          An Oracle simulating the same problem as ProbSimpleSO, plus
                       a tunable CPU burn per replication that does not affect
                       g()'s output.
ProbExpensiveHeavy     ProbExpensive at a heavy cost tier: burn_units=60,000.
ProbExpensiveLight     ProbExpensive at a light cost tier: burn_units=2,000.
ProbExpensiveModerate  ProbExpensive at a moderate cost tier: burn_units=20,000.     ExpensiveTester
ProbSimpleSO           An Oracle that simulates the Test Simple SO problem.          SimpleSOTester
ProbTPA                An Oracle that simulates Test Problem A.                      TPATester
ProbTPB                An Oracle that simulates Test Problem B.                      TPBTester
ProbTPC                An Oracle that simulates Test Problem C.                      TPCTester

The solve command

The PyMOSO solve command is for solving MOSO problems. Users can solve the built-in problems (use the listitems command to view the built-in problems), however, PyMOSO solve users typically will have their own MOSO problem they wish to solve. Thus, we assume users have implemented a PyMOSO oracle named MyProblem in myproblem.py. In the examples that follow, we assume the MyProblem implementation below, which is a bi-objective oracle with one-dimensional feasible points. See Implementing PyMOSO Oracles for instructions on implementing a MOSO problem as a PyMOSO oracle.

The Example Oracle

# import the Oracle base class
from pymoso.chnbase import Oracle

class MyProblem(Oracle):
    '''Example implementation of a user-defined MOSO problem.'''
    def __init__(self, rng):
        '''Specify the number of objectives and dimensionality of points.'''
        self.num_obj = 2
        self.dim = 1
        super().__init__(rng)

    def g(self, x, rng):
        '''Check feasibility and simulate objective values.'''
        # feasible values for x in this example
        feas_range = range(-100, 101)
        # initialize obj to empty and is_feas to False
        obj = []
        is_feas = False
        # check that dimensions of x match self.dim
        if len(x) == self.dim:
            is_feas = True
            # then check that each component of x is in the range above
            for i in x:
                if not i in feas_range:
                    is_feas = False
        # if x is feasible, simulate the objectives
        if is_feas:
            #use rng to generate random numbers
            z0 = rng.normalvariate(0, 1)
            z1 = rng.normalvariate(0, 1)
            obj1 = x[0]**2 + z0
            obj2 = (x[0] - 2)**2 + z1
            obj = (obj1, obj2)
        return is_feas, obj

The template solve command is pymoso solve oracle solver x0, where oracle is a built-in or user-defined oracle, solver is a built-in or user-defined algorithm, and x0 is a feasible starting point for the solver, with a space between each component. As a first example, we solve the user-defined MyProblem using the built-in R-PERLE starting at the feasible point 97.

pymoso solve myproblem.py RPERLE 97

Similarly, we can solve built-in problems, such as ProbTPA which has two-dimensional feasible points.

pymoso solve ProbTPA RPERLE 40 40

Henceforth, we present solve examples only for solving MyProblem. Since MyProblem is bi-objective, we recommend using the R-PERLE solver. However, for two or more objectives, PyMOSO implements R-MinRLE.

pymoso solve myproblem.py RMINRLE 97

For a single objective problem, PyMOSO implements R-SPLINE. We remark that if given a multi-objective problem, R-SPLINE will simply minimize the first objective. We do not necessarily prohibit such use, but urge that users take care when using R-SPLINE to minimize one objective of a many-objective problem.

pymoso solve myproblem.py RSPLINE 97

Regardless of the chosen solver, PyMOSO creates a new sub-directory of the working directory containing output. There will be a metadata file, indicating the date, time, solver, problem, and any other specified options. In addition, PyMOSO creates a file containing the solver-generated solution. PyMOSO provides additional options for users solving MOSO problems. We present examples of each option below. First, users can specify the name of the output directory.

pymoso solve --odir=OutDirectory myproblem.py RPERLE 45

Users can specify the simulation budget, which is currently set to a default of 200.

pymoso solve --budget=100000 myproblem.py RPERLE 12

Users may specify to take simulation replications in parallel. We only recommend doing so if the user has thought through appropriate pseudo-random number stream control issues (see Implementing PyMOSO Oracles). Furthermore, due to the overhead of parallelization, we only recommend using the parallel simulation replications feature if observations are sufficiently "expensive" to compute, e.g. the simulation takes a half second or more to generate a single observation. We remark that the run-time complexity of the simulation oracle may not perfectly indicate when it is appropriate to use parallelization; other factors include, e.g., the total simulation budget.

pymoso solve --simpar=4 myproblem.py RPERLE 44

Currently, all PyMOSO solvers support using common random numbers. Users may enable the functionality using the --crn option.

pymoso solve --crn myproblem.py RMINRLE 62

We do not recommend this option unless the oracle is implemented to be compatible, that is, the oracle uses PyMOSO's pseudo-random number generator to generate pseudo-random numbers or to provide a seed to an external mrg32k3a generator (see Implementing PyMOSO Oracles).

Users may specify an initial seed to PyMOSO's mrg32k3a pseudo-random number generator. Seeds must be 6 positive integers with spaces. The default is 12345 for each of the 6 components.

pymoso solve --seed 1111 2222 3333 4444 5555 6666 myproblem.py RPERLE 23

Users may specify algorithm-specific parameters (see the papers in which the algorithms were introduced for detailed explanations of the parameters.) All parameters are specified in the form --param name value. For example, the RLE relaxation parameter can be specified and set as betadel to a real number. We refer the reader to the table for the full list of currently available algorithm-specific parameters.

pymoso solve --param betadel 0.2 myproblem.py RPERLE 34

Finally, users may specify any number of options in one invocation. However, all options must be specified in after the solve command and before the myproblem.py argument. Furthermore, any --param options must be the last options. (Note that the \ at the end of the first line continues the command to the second line.)

pymoso solve --crn --simpar=4 --budget=10000 --seed 1 2 3 4 5 6 \
--odir=Exp1 --param mconst 4 --param betadel 0.7 myproblem.py RPERLE 97

Table of Algorithm-Specific Parameters

Parameter Name	Default Value	Affected Solvers	Description
mconst	2	RPERLE, RMINRLE, RPE, RSPLINE	Initialize the sample size and subsequent schedule of sample sizes.
bconst	8	RPERLE, RMINRLE, RPE, RSPLINE	Initialize the search sampling limit and subsequent schedule of limits.
radius	1	RPERLE, RMINRLE, RPE, RSPLINE	Sets radius that determines a point's neighborhood.
betadel	0.5	RPERLE, RMINRLE	Roughly, affects how likely it is for RLE to keep its given solution.
betaeps	0.5	RPERLE, RPE	Roughly, affects how likely PE will perform a search from a point.

The testsolve Command

The PyMOSO testsolve command tests algorithms on problems using a PyMOSO tester. Users can test built-in or user-defined solvers with built-in or user-defined testers. In the examples that follow, we assume users have implemented MyProblem as in The Example Oracle and the corresponding tester named MyTester in mytester.py, shown in The Example Tester. See Implementing PyMOSO Testers for instructions on implementing a user-defined tester, including a metric for comparing algorithms, in PyMOSO.

The Example Tester

import sys, os
sys.path.insert(0,os.path.dirname(__file__))
# use hausdorff distance (dh) as an example metric
from pymoso.chnutils import dh
# import the MyProblem oracle
from myproblem import MyProblem

# optionally, define a function to randomly choose a MyProblem feasible x0
def get_ranx0(rng):
    val = rng.choice(range(-100, 101))
    x0 = (val, )
    return x0

# compute the true values of x, for computing the metric
def true_g(x):
    '''Compute the objective values.'''
    obj1 = x[0]**2
    obj2 = (x[0] - 2)**2
    return obj1, obj2

# define an answer as appropriate for the metric
myanswer = {(0, 4), (4, 0), (1, 1)}

class MyTester(object):
    '''Example tester implementation for MyProblem.'''
    def __init__(self):
        self.ranorc = MyProblem
        self.answer = myanswer
        self.true_g = true_g
        self.get_ranx0 = get_ranx0

    def metric(self, eles):
        '''Metric to be computed per retrospective iteration.'''
        epareto = [self.true_g(point) for point in eles]
        haus = dh(epareto, self.answer)
        return haus

The template testsolve command is pymoso testsolve tester solver where tester is a built-in or user-defined tester, and solver is a built-in or user-defined solver. Users may also specify an x0, as in the solve command, if the tester does not implement the function to generate feasible points. As a first example, we test RPERLE on MyProblem using MyTester. Since some options are compatible with both solve and testsolve, we include those options in this example.

pymoso testsolve --budget=999 --odir=exp1 \
--crn --seed 1 2 3 4 5 6 mytester.py RPERLE

Users may want to compute some metric on the algorithm-generated solutions. If a metric is defined as part of the tester, such as in MyTester, the testsolve command can compute the metric on every algorithm iteration using the --metric option.

pymoso testsolve --metric mytester.py RPERLE

The testsolve command cannot perform simulation replications in parallel. However, testers can apply the solvers to independent sample paths of the problems. For example, to test RPERLE on 100 independent sample paths of MyProblem, compute the metrics for each sample path, and use common random numbers in each sample path, use the following command.

pymoso testsolve --crn --metric --isp=100 mytester.py RPERLE

PyMOSO can perform independent algorithm runs in parallel. Use the --proc option to specify the number of processes available to PyMOSO.

pymoso testsolve --crn --metric --isp=100 --proc=20 mytester.py RPERLE

We remark here that, to ensure the algorithm runs remain independent using PyMOSO's pseudo-random number generator (see Implementing PyMOSO Oracles), researchers should set the total simulation budget so that the included algorithms do not surpass 200 retrospective approximation (RA) iterations. For reference, using the default settings, the sample size at every point in the 200th RA iteration is almost 380 million.

The testsolve command creates a results file for each independent sample path. The file contains the solutions generated at every algorithm iteration, such that the solution of iteration 2 is on line 2, iteration 10 on line 10, and so forth. If --metric is specified, PyMOSO generates a second file for each independent sample path containing the collection of triples (iteration number, simulations used at end of iteration, metric).

Implementing problems, testers, and algorithms in PyMOSO

To use PyMOSO, users solving MOSO problems must implement a PyMOSO oracle, and users testing MOSO algorithms should implement, at least, a PyMOSO oracle and tester. In this section, we provide template Python code to help users quickly implement oracles, testers, and perhaps solvers in PyMOSO.

Implementing PyMOSO Oracles

Usually, implementing a PyMOSO oracle implies implementing a Monte Carlo simulation oracle as a black box function while following the PyMOSO rules put forth in this section. For reference, we discuss the example PyMOSO oracle MyProblem in The Example Oracle. Users may copy the code in The Example Oracle and re-implement the function g as needed. We now list the basic requirements of every g implementation.

1. The function g must be an instance method of an Oracle sub-class, and thus take self as its first parameter.
2. The function g must take an arbitrarily-named second parameter which is a tuple of length self.dim and represents a point. Stylistically, PyMOSO consistently names this parameter x.
3. The function g must take an arbitrarily-named third parameter which is a modified Python random.Random object. Stylistically, PyMOSO consistently names this parameter rng.
4. The function g must return a boolean first and a tuple of length self.num_obj second.
   • The boolean is True if x is feasible, and False otherwise.
   • If x is feasible, the tuple contains a single observation of every objective. If x is not feasible, each element in the tuple is None.

If users already have an implemented simulation oracle, they may find it convenient to implement g as wrapper which calls that simulation from Python. As an example, suppose a user has implemented a simulation in C which is compiled to a C library called mysim.so and placed in the working directory. Suppose further that the simulation function takes the following as parameters: an array of integers representing a point and an unsigned integer representing the number of observations to take at x. The function output is defined as struct Simout with members feas set to 0 or 1, obj a double array set to the mean of the observed objective values, and var a double array set to the sample variance of the observed objective values. Then users can modify the template to wrap the C function struct Simout c_func(int x, int n) as in the example.

Example Oracle that Wraps a C Simulation

from ctypes import CDLL, c_double, c_uint, c_int, Structure
import os.path
libname = 'mysim.so'
libabspath = os.path.dirname(os.path.abspath(__file__)) + os.path.sep + dll_name
libobj = CDLL(libabspath)

class Simout(Structure):
    _fields_ = [("feas", c_int), ("obj", c_double*2), ("var", c_double*2)]
csimout = libobj.c_func
csimout.restype = Simout

from pymoso.chnbase import Oracle

class MyProblem(Oracle):
    '''Example implementation of a user-defined MOSO problem.'''
    def __init__(self, rng):
        '''Specify the number of objectives and dimensionality of points.'''
        self.num_obj = 2
        self.dim = 1
        super().__init__(rng)

    def g(self, x, rng):
        '''Check feasibility and simulate objective values.'''
        is_feasible = True
        objective_values = (None, None)
        # g takes only one observation so set the c_func parameter to 1
        c_n = c_uint(1)
        # c_func requires is an integer so convert it -- this is a 1D example
        c_x = c_int(x[0])
        # call the C function
        mysimout = csimout(c_x, c_n)
        if not mysimout.feas:
            is_feasible = False
        else:
            is_feasible = True
        if is_feasible:
            objective_values = tuple(mysimout.obj)
        return is_feasible, objective_values

The C wrapper example is a valid PyMOSO oracle which wraps a C function. However, PyMOSO algorithms cannot enable common random numbers on this oracle. Furthermore, PyMOSO cannot guarantee that observations are independent when taken in parallel. To enable these properties, the external simulation must use mrg32k3a as the generator and must accept a user-specified seed.

Suppose the library mysim.so also implements the function set_simseed which accepts a long array representing an mrg32k3a seed. We modify the wrapper for compatibility with common random numbers and to guarantee independence of parallel observations. The new wrapper demonstrates using rng.get_seed() to return the current mrg32k3a seed.

Example Wrapper with PyMOSO Random Numbers

from ctypes import CDLL, c_double, c_uint, c_int, Structure, c_long
import os.path
libname = 'mysim.so'
libabspath = os.path.dirname(os.path.abspath(__file__)) + os.path.sep + dll_name
libobj = CDLL(libabspath)

class Simout(Structure):
    _fields_ = [("feas", c_int), ("obj", c_double*2), ("var", c_double*2)]
csimout = libobj.c_func
csetseed = libobj.set_simseed
csimout.restype = Simout

from pymoso.chnbase import Oracle

class MyProblem(Oracle):
    '''Example implementation of a user-defined MOSO problem.'''
    def __init__(self, rng):
        '''Specify the number of objectives and dimensionality of points.'''
        self.num_obj = 2
        self.dim = 1
        super().__init__(rng)

    def g(self, x, rng):
        '''Check feasibility and simulate objective values.'''
        is_feasible = True
        objective_values = (None, None)
        # get the PyMOSO seed from rng
        seed = rng.get_seed()
        # convert the seed to c_long array
        c_longarr = c_long*6
        c_seed = c_longarr(seed[0], seed[1], seed[2], seed[3], seed[4], seed[5])
        # use the library function to set the sim seed
        csetseed(c_seed)
        # g takes only one observation so set the c_func parameter to 1
        c_n = c_uint(1)
        # c_func requires is an integer so convert it -- this is a 1D example
        c_x = c_int(x[0])
        # call the C function
        mysimout = csimout(c_x, c_n)
        if not mysimout.feas:
            is_feasible = False
        else:
            is_feasible = True
        if is_feasible:
            objective_values = tuple(mysimout.obj)
        return is_feasible, objective_values

Alternatively, if the number of required pseudo-random numbers is known, users can use rng.random() to generate pseudo-random numbers and then pass them to an external simulation if such functionality is supported.

The rng object is implemented as a sub-class of Python's random.Random class, thus the official Python documentation for random applies to rng and is found at https://docs.python.org/3/library/random.html. In addition to rng using mrg32k3a as its generator, we also implement rng.normalvariate such that it uses the Beasley-Springer-Moro algorithm (Law 2015, p. 458) to approximate the inverse of the standard normal cumulative distribution function.

When using rng, to ensure independent sampling of observations, PyMOSO "jumps" forward in the pseudo-random number stream after obtaining every simulation replication. Each jump is of fixed size 2^76 pseudo-random numbers. Thus, we require that every simulation replication use fewer than 2^76 pseudo-random numbers. We ensure independence among parallel replications by "giving" each processor a stream (an rng), each of which is 2^127 pseudo-random numbers apart. When using the current PyMOSO algorithms that rely on RA, each RA iteration begins the next available independent stream 2^127, where PyMOSO accounts for the possibility of parallel computation within an RA iteration. Thus, in a given RA iteration, a user may simulate 100 million points at a sample size of 1 million, without common random numbers, and easily not reach the limit.

Implementing PyMOSO Testers

Consider again the example tester. As a minimal valid PyMOSO tester, users may do nothing but assign the MyTester member self.ranorc to a PyMOSO oracle, such as MyProblem, in Line 27. However, we expect most users to leverage PyMOSO features by implementing metrics and feasible point generators. The function get_ranx0 allows the tester to generate feasible points to MyProblem and metric allows the tester to compute a metric on sets returned by a solver. Researchers may implement any number of additional supporting functions, including members and methods of the tester class. The true_g function is an example of such a supporting function, which is used to compute the example metric.

First, we list the rules for implementing a feasible point generator.

1. The function is arbitrarily named but must be set to the self.get_ranx0 member of a tester.
2. The function must take a single parameter, an arbitrarily named random.Random object we suggest naming rng.
3. The function must return a tuple with length corresponding to the self.dim member of the self.ranorc member of the tester.

Since a researcher's desired metric depends on the algorithm capabilities and problem complexity, PyMOSO allows researchers to implement any metric they choose. We provide three example metrics, but first, we list the implementation rules of the metric function.

1. The \inline{metric} function must be an instance method of a tester, and thus take \inline{self} as its first parameter.
2. The second parameter of \inline{metric} is arbitrarily named and is a Python set of tuples.
3. PyMOSO does not enforce the return value of \inline{metric}, but we recommend a scalar real number.

The metric implemented in the example tester is the Hausdorff distance from (a) the true image of an estimated solution returned by an algorithm, to (b) the true solution hard-coded as myanswer.

For an example of a different metric, consider a MOSO problem that has more than one local efficient set (LES) and such that each LES contains no members of another LES. Since an algorithm that converges to a LES is may find only one LES, we may define the metric to compute the Hausdorff distance between the true image of the estimated solution and the "closest" true LES, as follows. Let self.answer be implemented as a list of sets, and assume a self.true_g implementation. Then the following example implements the described metric.

Example Metric 1

def metric(self, eles):
    # use the distance to the closest set.
    epareto = [self.true_g(point) for point in eles]
    # self.soln is a list of sets
    dist_list = [dh(epareto, les) for les in self.answer]
    return min(dist_list)

For single-objective problems with one correct solution x*, a simple metric that takes an estimated solution X is |g(X) - g(x*)|, which we implement in the second example metric assuming an appropriate implementation of self.answer and self.true_g.

Example Metric 2

def metric(self, singleton_set):
    # single objective algorithms still return a set
    point, = singleton_set
    # let self.soln be a real number
    dist = abs(self.true_g(point) - self.answer)
    return dist

Implementing PyMOSO Algorithms

Researchers can implement simulation optimization algorithms in the PyMOSO framework. PyMOSO provides support for algorithms in three categories:

1. PyMOSO provides strong support for implementing new MOSO algorithms that rely on RLE in an RA framework.
2. PyMOSO provides strong support for implementing general RA algorithms.
3. PyMOSO provides basic support, such as pseudo-random number control, for implementing other simulation optimization algorithms. We provide templates of algorithms implemented in each of these three categories, along with example code snippets.

In the first category, programmers can use PyMOSO to create new RA algorithms that use RLE for convergence. The novel part of these algorithms, created by the user, will be the accel function which should collect points to send to RLE for certification. Here, we list the rules for accel.

1. The accel function must be an instance method of an RLESolver object, and thus its first parameter must be self.
2. The second parameter is arbitrarily named and is a set of tuples. We recommend naming the parameter warm_start, as it represents the sample-path solution of the previous RA iteration.
3. The return value must be a set of tuples representing feasible points; we do not recommend any particular name. In every RA iteration, PyMOSO will first call accel(self, warm_start) and send the returned set to rle(self, candidate_les). The return value must be a set of tuples. The implementer does not need to implement or call RLE, as in the template accelerator.

Template Accelerator

from pymoso.chnbase import RLESolver

# create a subclass of RLESolver
class MyAccel(RLESolver):
    '''Example implementation of an RLE accelerator.'''

    def accel(self, warm_start):
        '''Return a collection of points to send to RLE.'''
        # implement algorithm logic here and return a set
        return warm_start

In the second category, algorithm designers can quickly implement any RA algorithm by sub-classing RASolver and implementing the spsolve function, as shown in the RA solver template. The algorithm can be a single-objective algorithm. PyMOSO cannot guarantee the convergence of such algorithms. The template is technically valid in PyMOSO but is probably not effective.

Template RA Solver

As written below, this template cannot reach a normal budget-exhausted stop when run via solve(), at any budget: spsolve never calls self.estimate, so self.num_calls never advances past 0. It does not run forever, though -- the retrospective-approximation loop's own iteration-dependent sample-size formula (calc_b) eventually overflows Python's float range and raises OverflowError (thousands of fast, simulation-free iterations in, not immediately, and naming neither this template nor the real cause). A real implementation must call self.estimate (directly, or via self.upsample/self.spline) to make progress against the budget; see KNOWN_ISSUES.md for the exact failure this produces as written.

from pymoso.chnbase import RASolver

class MyRAAlg(RASolver):
    '''Template implementation of an RA solver.'''

    def spsolve(self, warm_start):
        '''Return the sample path solution.'''
        # implement algorithm logic here and return a set
        return warm_start

Though analogous to those of an RLESolver.accel method, for completeness, we list the requirements for an RASolver.spsolve method.

1. The spsolve function must be an instance method of an RASolver object, and thus its first parameter must be self.
2. The second parameter is arbitrarily named and is a set of tuples. We recommend naming the parameter warm_start as it represents the sample-path solution of the previous RA iteration.
3. The return value must be a set of tuples representing feasible points; we do not recommend any particular name.

In the third category, PyMOSO can accommodate any simulation optimization algorithm by implementing the solve function of a MOSOSolver sub-class as shown. It does not have to be a multi-objective algorithm. PyMOSO will require users to send an initial feasible point x0 whether or not the algorithm needs it. The initial feasible point x0 is accessed through self.x0 which is a tuple. We now list the rules for implementing any MOSOSolver.solve function.

1. The solve function must be an instance method of MOSOSolver, and thus take self as its first parameter.
2. The second parameter is the simulation budget, a natural number.
3. The solve function must return a dictionary (we name it results in our example) with at least 3 keys: 'itersoln', 'simcalls', 'endseed'. Researchers may track additional data and add it to results as desired.
   • The 'itersoln' key itself corresponds to a dictionary with a key for each algorithm iteration labeled {0, 1, ...}. The value at each iteration is a set containing the estimated solution at the end of the iteration.
   • The 'simcalls' key itself corresponds to a dictionary with a key for each algorithm iteration labeled {0, 1, ...}. The value at each iteration is a natural number containing the cumulative number of simulation replications taken at the end of the iteration.
   • The 'endseed' key corresponds to a tuple of length 6, representing an mrg32k3a seed. The algorithm programmer should ensure the stream generated by results['endseed'] is independent of all streams used by the algorithm.

Researchers may use the MOSO template to implement new simulation optimization algorithms.

Template MOSO Solver

This is a structural sketch, not a runnable template: it shows the required shape of a MOSOSolver.solve() (loop against the budget, return a results dict) rather than a working algorithm. Do not call solve() on this class as written.

from pymoso.chnbase import MOSOSolver

class MyMOSOAlg(MOSOSolver):
    '''Structural sketch of a MOSO solver, illustrative only -- not meant to be run as written.'''

    def solve(self, budget):
        while self.num_calls <= budget:
            pass  # implement algorithm logic and return the results
        return results

For convenience, in the mini-sections below, we also provide some example code snippets that we find useful when implementing algorithms in PyMOSO. They work without modification when using the templates above that inherit \inline{RLESolver} or \inline{RASolver}, but some functions may require implementation or modification for use in a MOSOSolver. For reference, this section contains a list of most objects accessible to PyMOSO programmers.

Take Simulation Replications at a Point

from pymoso.chnutils import get_nbors
# pretend x has not yet been visited in this RA iteration and is feasible
x = (1, 1, 1)

# self.m is the sample size of the current RA iteration
m = self.m
# self.num_calls is the cumulative number of simulations used till now
start_num_calls = self.num_calls
# use estimate to sample x and put results in self.gbar and self.sehat
isfeas, fx, se = self.estimate(x)
calls_used = self.num_calls - start_num_calls
print(m == calls_used) # True
print(fx == self.gbar[x]) # True
print(se == self.sehat[x]) # True

# estimate will not simulate again in subsequent visits to a point
start_num_calls = self.num_calls
isfeas, fx, se = self.estimate(x)
calls_used = self.num_calls - start_num_calls
print(calls_used == 0) # True

Find Neighbors and Take Simulation Replications

# neighborhood radiuss
r = self.nbor_rad
nbors = get_nbors(x, r)
self.upsample(nbors)
for n in nbors:
  print(n in self.gbar) # True if n feasible else False

# upsample also returns the feasible subset
nbors = self.upsample(nbors)

Argsort a Dictionary of Points

# 0 index for first objective
sorted_feas = sorted(nbors | {x}, key=lambda t: self.gbar[t][0])

Select the Minimizer and its Value

xmin = sorted_feas[0]
fxmin = self.gbar[xmin]

Use SPLINE to Retrive a Local Minimizer

# no constraints and minimize the 2nd objective
x0 = (2, 2, 2)
isfeas, fx, sex = self.estimate(x0)
# the suppressed value is the set visited along SPLINE's trajectory
_, xmin, fxmin, sexmin = self.spline(x0, float('inf'), 1, 0)
print(self.gbar[xmin] == fxmin) # True

Find the Non-Dominated Points in a Dictionary

from pymoso.chnutils import get_nondom
nondom = get_nondom(self.gbar)

Randomly Select Points in a Set

solver_rng = self.sprn
# pick up to 5 points -- returns a list, not a set.
ran_pts = solver_rng.sample(list(nondom), min(5, len(nondom)))
one_in_five = solver_rng.choice(ran_pts)

Using solve and testsolve in Python Programs

Users may the solve and testsolve functions within a Python program.

Minimal solve Example

# import the solve function
from pymoso.chnutils import solve
# import the module containing the RPERLE implementation
import pymoso.solvers.rperle as rp
# import MyProblem - myproblem.py should usually be in the script directory
import myproblem as mp

# specify an x0. In MyProblem, it is a tuple of length 1
x0 = (97,)
soln = solve(mp.MyProblem, rp.RPERLE, x0)
print(soln)

Some solve Examples with Options

# example for specifying budget and seed
budget=10000
seed = (111, 222, 333, 444, 555, 666)
soln1 = solve(mp.MyProblem, rp.RPERLE, x0, budget=budget, seed=seed)

# specify crn and simpar
soln2 = solve(mp.MyProblem, rp.RPERLE, x0, crn=True, simpar=4)

# specify algorithm specific parameters
soln3 = solve(mp.MyProblem, rp.RPERLE, x0, radius=2, betaeps=0.3, betadel=0.4)

# mix them
soln4 = solve(mp.MyProblem, rp.RPERLE, x0, crn=True, seed=seed, radius=5)

A testsolve Example

# import the testsolve functions
from pymoso.chnutils import testsolve
# import the module containing RPERLE
import pymoso.solvers.rperle as rp
# import the MyTester class
from mytester import MyTester

# testsolve always requires an x0 positionally, even though it is only
# actually used when ranx0=True is passed (otherwise every independent
# sample path uses this x0 unchanged); MyTester can generate its own
# per sample path via get_ranx0 if you pass ranx0=True below.
x0 = (1, )
run_data = testsolve(MyTester, rp.RPERLE, x0, isp=100, crn=True, radius=2)

Computing a Metric on testsolve Output

Programmers must compute their metric. Here, run_data is a dictionary of the form described here and we compute the metric on the 5th iteration of of the 12th independent algorithm instance.

iter5_soln = run_data[0][11]['itersoln'][4]
isp12_iter5_metric = MyTester().metric(iter5_soln)

PyMOSO Object Reference

The pymoso.prng.mrg32k3a Module

The pymoso.prng.mrg32k3a module exposes the pseudo-random number generator and functions to manipulate it.

Object	Description
MRG32k3a	Sub-class of random.Random, defines all rng objects.
get_next_prnstream(seed)	Return an rng object seeded 2^127 steps from the input seed.
jump_substream(rng)	Seed the input rng object 2^76 steps forward.

The pymoso.chnbase Module

The pymoso.chnbase module implements the base classes for oracles and solvers. Programmers should sub-class these when creating new PyMOSO implementations.

Class	Description
Oracle	Base class for implementing oracles.
RLESolver	Base class for implementing solvers using RLE.
RASolver	Base class for implementing RA solvers.
MOSOSolver	Base class for all solvers.

The pymoso.chnutils Module

The pymoso.chnutils contains useful functions for implementing and testing algorithms.

Function	Description
solve(oracle, solver, x0, **kwargs)	See here for instructions.
testsolve(tester, solver, x0, **kwargs)	See here for instructions.
does_weak_dominate(g, h, relg, relh)	All inputs are tuples of equal length. Returns True if g weakly dominates h with the relaxations.
does_dominate(g, h, relg, relh)	Returns True if g dominates h with the relaxations.
does_strict_dominate(g, h, relg, relh)	Returns True if g strictly dominates h with the relaxations.
get_nondom(obj_dict)	Input: a dictionary with tuples for keys and values. The keys are feasible points; the values are their objective values. Return: a set of tuples representing non-dominated points.
get_nbors(x, r)	Input: a tuple x, a positive real scalar r indicating the neighborhood radius. Return: Set of tuples which are the neighbors.
get_setnbors(S, r)	Input: a set of tuples, and the neighborhood radius. Return: the union of get_nbors(s, r) for every s in S.
dh(A, B)	Returns the Hausdorff distance between set A and set B.
edist(x1, x2)	Return the Euclidean distance from x1 to x2.
get_metric(results, tester)	Input: results is a dictionary, the output of each sample path of testsolve. tester must implement metric. Returns: The set of triples (iteration, simulation count, metric) for an algorithm run.

The Oracle Class

When implementing RASolver algorithms, programmers may not need to access Oracle objects directly at all. When implementing MOSOSolver algorithms, programmers will use (or wrap) hit and crn_advance().

Member/Method	Description
num_obj	A positive integer, the number of objectives.
dim	A positive integer, the dimensionality of feasible points.
rng	An instance of MRG32k3a.
hit(x, n)	Take n observations of x. Return: True, and a tuple containing the mean of the observations for each objective, and a tuple containing the standard error for each objective if x is feasible. The function handles CRN internally.
set_crnflag(bool)	Turn CRN on (True) or off.
set_crnold(state)	Save the rng state as the CRN baseline, e.g. for an algorithm iteration. Get the state using rng.getstate().
crn_reset()	Back the oracle rng to the CRN baseline.
crn_advance()	If CRN is on, reset, and then jump to the next independent pseudo-random stream and save the new baseline, e.g. before starting a new algorithm iteration.
crn_setobs()	Set an intermediate CRN for individual oracle observations.
crn_nextobs()	Jump the rng forward, e.g. after taking an observation, and crn_setobs the seed.
crn_check()	f CRN is on, return to the baseline. Otherwise, use crn_nextobs before taking the next observation.

The MOSOSolver Class

The class provides a basic structure for implementing new MOSO algorithms in PyMOSO.

Member/Method	Description
orc	The oracle object for the solver to solve.
dim	Number of dimensions of points in the feasible set. Should match self.orc
num_obj	Similarly, the number of objectives in self.orc.
num_calls	A running count of the number of simulation replications taken of self.orc.
x0	A feasible starting point. This point is additionally supplied to algorithms that don't need one.

The RASolver Class

The class implements a common structure for all RA algorithms, including: caching of simulation replications, scheduling and updating of sample sizes and limits, and a wrapper to Oracle.hit.

Member/Method	Description
sprn	An instance of MRG32k3a for the solver to use.
nbor_rad	The neighborhood radius used by solvers seeking local optimality.
gbar	A dictionary where every key and value is a tuple. The keys are feasible points, values are their objective values. gbar is "wiped" every retrospective iteration.
sehat	Exactly like gbar except the values are standard errors.
m	The sample size of the current iteration.
calc_m(nu)	Compute the sample size of the current iteration. RA algorithms automatically do this every iteration and assign the value to `m'.
b	The searching sample limit of the current iteration.
calc_b(nu)	Exactly as calc_m but for the searching sample limit.
estimate(x, c, obj)	The estimate function is essentially a smart wrapper for self.orc.hit. Inputs: tuple x to sample, c a feasibility constraint, obj the objective to constrain. Return: same as Oracle.hit. Retrieves or saves the results from/to gbar and sehat as appropriate. Returns not feasible if the otherwise feasible result is not less than the constraint.
upsample(S)	A version of estimate for sets. Returns the feasible subset of S.
spline(x, c, obmin, obcon)	Return a sample path local minimizer. Input: a feasible start, constraint, objective to minimize, objective to constrain. Return: a set of tuples of the trajectory, the minimizer tuple, the minimum tuple, the standard error tuple.

The RLESolver Class

The sub-class of \inlinetwo{RASolver} adds RLE and its relaxation.

Member/Method	Description
betadel	Affects the relaxation values computed in RLE.
calc_delta(se)	Computes the RLE relaxation given a standard error, using self.m and self.betadel.
rle(candidate_les)	Input: set of tuples, Returns: set of tuples. Finds the LES at sample size self.m.