Python dev litong

demo/demo_problems_solvers.py

@@ -11,6 +11,13 @@
 
 # Import the ProblemsSolvers class and other useful functions
 from simopt.experiment_base import ProblemsSolvers, plot_solvability_profiles
+from simopt.models.san import SANLongestPath
+from simopt.models.smf_origin import SMF_Max
+from simopt.models.rmitd import RMITDMaxRevenue
+from simopt.models.mm1queue import MM1MinMeanSojournTime
+# from simopt.models.san_1 import SANLongestPath1
+from simopt.experiment_base import ProblemSolver, plot_area_scatterplots, post_normalize, plot_progress_curves, plot_solvability_cdfs, read_experiment_results, plot_solvability_profiles, plot_terminal_scatterplots, plot_terminal_progress
+
 
 # !! When testing a new solver/problem, first go to directory.py.
 # There you should add the import statement and an entry in the respective
@@ -20,28 +27,42 @@
 # Specify the names of the solver and problem to test.
 # These names are strings and should match those input to directory.py.
 # Ex:
-solver_names = ["RNDSRCH", "ASTRODF", "NELDMD"]
-problem_names = ["CNTNEWS-1", "SAN-1"]
+# solver_names = ["ASTRODF", "Boom-PGD", "Boom-FW", "RNDSRCH", "GASSO", "NELDMD"]
+solver_names = ['ACTIVESET', 'Boom-FW', 'Boom-PGD', 'PGD-SS']
+problem_names = ["OPENJACKSON-1", 'SAN-1', 'SMFCVX-1', 'SMF-1', 'CASCADE-1', 'NETWORK-1']
+# problem_names = ["DYNAMNEWS-1", "SSCONT-1", "SAN-1"] "OPENJ-1"
+# problem_names = ["SMF-1", "SAN-1", "RMITD-1", "MM1-1"]
+# problems = [SANLongestPath, SMF_Max, RMITDMaxRevenue, MM1MinMeanSojournTime]
 
 
 # Initialize an instance of the experiment class.
 mymetaexperiment = ProblemsSolvers(solver_names, problem_names)
+# mymetaexperiment = ProblemsSolvers(solver_names=solver_names, problems = problems)
 
-# Write to log file.
-mymetaexperiment.log_group_experiment_results()
+n_solvers = len(mymetaexperiment.experiments)
+n_problems = len(mymetaexperiment.experiments[0])
 
 # Run a fixed number of macroreplications of each solver on each problem.
-mymetaexperiment.run(n_macroreps=3)
+mymetaexperiment.run(n_macroreps=20)
+
 
 print("Post-processing results.")
 # Run a fixed number of postreplications at all recommended solutions.
-mymetaexperiment.post_replicate(n_postreps=50)
+mymetaexperiment.post_replicate(n_postreps=20)
 # Find an optimal solution x* for normalization.
-mymetaexperiment.post_normalize(n_postreps_init_opt=50)
+mymetaexperiment.post_normalize(n_postreps_init_opt=20)
 
 print("Plotting results.")
 # Produce basic plots of the solvers on the problems.
-plot_solvability_profiles(experiments=mymetaexperiment.experiments, plot_type="cdf_solvability")
+plot_solvability_profiles(experiments=mymetaexperiment.experiments, plot_type="cdf_solvability") # cdf_solvability
+# plot_solvability_profiles(experiments=mymetaexperiment.experiments, plot_type="diff_quantile_solvability", ref_solver='RNDSRCH', all_in_one=True, plot_CIs=True, print_max_hw=True) 
+# plot_terminal_scatterplots(experiments=mymetaexperiment.experiments, all_in_one=True)
+
+# Plot the mean progress curves of the solvers on the problems.
+CI_param = True
+for i in range(n_problems):
+    plot_progress_curves([mymetaexperiment.experiments[solver_idx][i] for solver_idx in range(n_solvers)], plot_type="mean", all_in_one=True, plot_CIs=CI_param, print_max_hw=True)
+
 
 # Plots will be saved in the folder experiments/plots.
 print("Finished. Plots can be found in experiments/plots folder.")

docs/gasso.rst

@@ -0,0 +1,61 @@
+Solver: Gradient-Based Adaptive Stochastic Search for Simulation Optimization Over Continuous Space (GASSO)
+================================================================
+
+Description:
+------------
+The solver iteratively generates population of candidate solutions from a sample distribution,
+and uses the performance of sample distribution to update the sampling dsitribution. 
+GASSO has two stages in each iteration: 
+1. Stage I: Generate candidate solutions from some exponential family of distribution, and the 
+2. Stage II: Evaluate candidate solutions, and update the parameter of sampling distribution via 
+direct gradient search. 
+
+Scope:
+------
+* objective_type: single
+
+* constraint_type: box
+
+* variable_type: continuous
+
+Solver Factors:
+---------------
+* crn_across_solns: Use CRN across solutions?
+
+    * Default: True
+
+* N: Number of candidate solutions
+
+    * Default: :math:`50 * \sqrt(dim)`
+
+* M: Number of function evaluations per candidate
+
+    * Default: 10
+
+* K: Number of iterations
+
+    * Default: Budget/(N * M)
+
+* alpha_0: Determines the initial step size
+
+    * Default: 50
+
+* alpha_c: Determines the speed at which the step size decreases
+
+    * Default: 1500
+
+* alpha_p: Determines the rate at which step size gets smaller
+
+    * Default: 0.6
+
+* alpha_k: Step size
+
+    * Default: :math:`\frac{alpha_0}{(k + \alpha_c) ^ {alpha_p}}`
+
+
+References:
+===========
+This solver is adapted from the article Enlu Zhou, Shalabh Bhatnagar (2018).
+Zhou, E., & Bhatnagar, S. (2017). Gradient-based adaptive stochastic search for simulation optimization over continuous space. 
+*INFORMS Journal on Computing, 30(1), 154-167.  
+(https://doi.org/10.1287/ijoc.2017.0771)

docs/openjackson.rst

@@ -0,0 +1,117 @@
+Model: Open Jackson Network
+===============================================
+
+Description:
+------------
+This model represents an Open Jackson Network with Poisson arrival time, exponential service time, and probabilistic routing.
+
+Sources of Randomness:
+----------------------
+There are 3 sources of randomness in this model:
+1. Exponential inter-arrival time of customers at each station.
+2. Exponential service time of customers at each station.
+3. Routing of customers at each station after service.
+
+Model Factors:
+--------------
+* number_queues: The number of queues in the network.
+    * Default: 3
+
+* arrival_alphas: The rate parameter of the exponential distribution for the inter-arrival time of customers at each station.
+    * Default: [1,1,1,1,1]
+
+* service_mus: The rate parameter of exponential distribution for the service time of customers at each station.
+    * Default: [2,2,2,2,2]
+
+* routing_matrix: The routing probabilities for a customer at station i to go to service j after service. 
+    The departure probability from station i is :math: `1 - \sum_{j=1}^{n} (P_{ij})`
+    where n is the number of stations, and P is the routing matrix.
+    * Default: [[0.1, 0.1, 0.2, 0.2, 0],
+                [0.1, 0.1, 0.2, 0.2, 0],
+                [0.1, 0.1, 0, 0.1, 0.3],
+                [0.1, 0.1, 0.1, 0, 0.3],
+                [0.1, 0.1, 0.1, 0.1, 0.2]]
+
+* t_end: The time at which the simulation ends.
+    * Default: 200
+
+* warm_up: The time at which the warm-up period ends. Relevant only when steady_state_initialization is False.
+    * Default: 100
+
+* steady_state_initialization: Whether to initialize with queues sampled from steady state. 
+    If so, we sample geometric distribution with parameter lambdas/service_mus for each queue and initialize the queues with the sample.
+    * Default: True
+
+Below Factors are only relevant when creating random instances of the Model
+
+* density_p: The probability of an edge existing in the graph in the random instance. Higher the value, denser the graph.
+    * Default: 0.5
+
+* random_arrival_parameter: The parameter for the random arrival rate exponential distribution when creating a random instance.
+    * Default: 1
+
+
+Responses:
+----------
+* average_queue_length: The time-average queue length at each station.
+
+References:
+===========
+This model is adapted from Jackson, James R. (1957).
+"Networks of waiting lines". Operations Research. 4 (4): 518–521.
+(doi:10.1287/opre.5.4.518)
+
+Optimization Problem: OpenJacksonMinQueue (OPENJACKSON-1)
+================================================================
+
+Decision Variables:
+-------------------
+* service_mus
+
+Objectives:
+-----------
+Minimize the sum of average queue length at each station.
+
+Constraints:
+------------
+We require that the sum of service_mus at each station to be less than service_rates_budget.
+
+Problem Factors:
+----------------
+* budget: Max # of replications for a solver to take.
+
+  * Default: 1000
+
+* service_rates_budget: Total budget to be allocated to service_mus_budget.
+
+  * Default: 150
+
+Below factors are only relevant when creating random instances of the Problem
+
+* gamma_mean: Scale of the mean of gamma distribution when generating service rates upper bound in random instances.
+
+  * Default: 0.5
+
+* gamma_scale: Shape of gamma distribution when generating service rates upper bound in random instances.
+
+  * Default: 5
+
+Fixed Model Factors:
+--------------------
+* N/A
+
+Starting Solution: 
+------------------
+* initial_solution: lambdas * (service_rates_budget/sum(lambdas))
+
+Random Solutions: 
+-----------------
+Sample a Dirichlet distribution that sum to service_rates_budget - sum(lambdas). Then add lambdas to the sample.
+
+Optimal Solution:
+-----------------
+Unknown
+
+Optimal Objective Function Value:
+---------------------------------
+Unknown

simopt/base.py

@@ -3,11 +3,19 @@
 Summary
 -------
 Provide base classes for solvers, problems, and models.
+This is the modified version to generate and run random model/random problem instance.
 """
 
 import numpy as np
 from copy import deepcopy
-from mrg32k3a.mrg32k3a import MRG32k3a
+import sys
+import os.path as o
+# from mrg32k3a.mrg32k3a import MRG32k3a
+from rng.mrg32k3a import MRG32k3a #when do the multinomial, change to the local
+
+from simopt.auto_diff_util import bi_dict, replicate_wrapper
+
+sys.path.append(o.abspath(o.join(o.dirname(sys.modules[__name__].__file__), "..")))
 
 
 class Solver(object):
@@ -378,7 +386,7 @@ def check_factor_datatype(self, factor_name):
         is_right_type = isinstance(self.factors[factor_name], self.specifications[factor_name]["datatype"])
         return is_right_type
 
-    def attach_rngs(self, rng_list):
+    def attach_rngs(self, random_rng, copy=True):
         """Attach a list of random-number generators to the problem.
 
         Parameters
@@ -387,7 +395,25 @@ def attach_rngs(self, rng_list):
             List of random-number generators used to generate a random initial solution
             or a random problem instance.
         """
-        self.rng_list = rng_list
+        if copy:
+            self.random_rng = [deepcopy(rng) for rng in random_rng]
+        else:
+            self.random_rng = random_rng
+
+    def rebase(self, n_reps):
+        """Rebase the progenitor rngs to start at a later subsubstream index.
+
+        Parameters
+        ----------
+        n_reps : int
+            Substream index to skip to.
+        """
+        new_rngs = []
+        for rng in self.random_rng:
+            stream_index = rng.s_ss_sss_index[0]
+            substream_index = rng.s_ss_sss_index[1]
+            new_rngs.append(MRG32k3a(s_ss_sss_index=[stream_index, substream_index, n_reps]))
+        self.random_rng = new_rngs
 
     def vector_to_factor_dict(self, vector):
         """
@@ -622,8 +648,6 @@ def simulate(self, solution, m=1):
                 # to those of deterministic components of objectives.
                 solution.objectives[solution.n_reps] = [sum(pairs) for pairs in zip(self.response_dict_to_objectives(responses), solution.det_objectives)]
                 if self.gradient_available:
-                    # print(self.response_dict_to_objectives_gradients(vector_gradients))
-                    # print(solution.det_objectives_gradients)
                     solution.objectives_gradients[solution.n_reps] = [[sum(pairs) for pairs in zip(stoch_obj, det_obj)] for stoch_obj, det_obj in zip(self.response_dict_to_objectives_gradients(vector_gradients), solution.det_objectives_gradients)]
                     # solution.objectives_gradients[solution.n_reps] = [[sum(pairs) for pairs in zip(stoch_obj, det_obj)] for stoch_obj, det_obj in zip(self.response_dict_to_objectives(vector_gradients), solution.det_objectives_gradients)]
                 if self.n_stochastic_constraints > 0:
@@ -755,6 +779,21 @@ def check_factor_datatype(self, factor_name):
         """
         is_right_type = isinstance(self.factors[factor_name], self.specifications[factor_name]["datatype"])
         return is_right_type
+
+    def attach_rng(self, random_rng, copy=True):
+        """Attach a list of random-number generators to the problem.
+
+        Parameters
+        ----------
+        rng_list : list [``mrg32k3a.mrg32k3a.MRG32k3a``]
+            List of random-number generators used to generate a random initial solution
+            or a random problem instance.
+        """
+        # self.random_rng = random_rng
+        if copy:
+            self.random_rng = [deepcopy(rng) for rng in random_rng]
+        else:
+            self.random_rng = random_rng 
 
     def replicate(self, rng_list):
         """Simulate a single replication for the current model factors.
@@ -772,6 +811,27 @@ def replicate(self, rng_list):
             Gradient estimate for each response.
         """
         raise NotImplementedError
+
+
+class Auto_Model(Model):
+    """
+    Subclass of Model. 
+    """
+    def __init__(self, fixed_factors):
+        # set factors of the simulation model
+        # fill in missing factors with default values
+        super(Auto_Model, self).__init__(fixed_factors)
+        self.differentiable_factor_names = []
+        for key in self.specifications:
+            if self.specifications[key]["datatype"] == float:
+                self.differentiable_factor_names.append(key)
+        self.bi_dict = bi_dict(self.response_names)
+
+    def innner_replicate(self, rng_list):
+        raise NotImplementedError
+
+    def replicate(self, rng_list, **kwargs):
+        return replicate_wrapper(self, rng_list, **kwargs)
 
 
 class Solution(object):
@@ -826,7 +886,10 @@ class Solution(object):
     def __init__(self, x, problem):
         super().__init__()
         self.x = x
-        self.dim = len(x)
+        if isinstance(x, int) or isinstance(x, float):
+            self.dim = 1
+        else:
+            self.dim = len(x)
         self.decision_factors = problem.vector_to_factor_dict(x)
         self.n_reps = 0
         self.det_objectives, self.det_objectives_gradients = problem.deterministic_objectives_and_gradients(self.x)