Implement the DDS algorithm for spotpy

.gitignore

@@ -32,6 +32,9 @@ wheels/
 .installed.cfg
 *.egg
 
+# Python Unittest tool `py.test`
+.pytest_cache/*
+
 # PyInstaller
 #  Usually these files are written by a python script from a template
 #  before PyInstaller builds the exe, so as to inject date/other infos into it.

spotpy/algorithms/__init__.py

@@ -30,4 +30,5 @@
 from .abc import abc         # Artificial Bee Colony
 from .fscabc import fscabc   # Fitness Scaling Artificial Bee Colony
 from .dream import dream     # DiffeRential Evolution Adaptive Metropolis
-from .list_sampler import list_sampler # Samples from given spotpy database
+from .list_sampler import list_sampler  # Samples from  given spotpy database
+from .dds import dds         # Dynamically Dimensioned Search algorithm by Bryan Tolson.

spotpy/examples/dds/benchmark_dds.py

@@ -0,0 +1,117 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+from pprint import pprint
+import numpy as np
+import matplotlib.pylab as plt
+import json
+
+import time
+
+try:
+    import spotpy
+except ImportError:
+    import sys
+
+    sys.path.append(".")
+    import spotpy
+
+from spotpy.examples.spot_setup_hymod_python import spot_setup
+
+
+
+spot_setup = spot_setup()
+
+# Create samplers for every algorithm:
+results = []
+
+benchmarks_dict = []
+benchmarks_duration = {"dds":[], "sceua":[], "dds_like":[],"sceua_like":[]}
+reps = [300, 1000, 3000, 4000, 5000, 10000]
+
+
+for rep in reps:
+
+    timeout = 10  # Given in Seconds
+
+    parallel = "seq"
+    dbformat = "csv"
+
+    start = time.time()
+    dds_sampler = spotpy.algorithms.DDS(spot_setup, parallel=parallel, dbname='DDS', dbformat=dbformat, sim_timeout=timeout)
+    dds_sampler.sample(rep, trials=1)
+    results.append(dds_sampler.getdata())
+    dds_elapsed = time.time() - start
+
+    start = time.time()
+    sceua_sampler = spotpy.algorithms.sceua(spot_setup, parallel=parallel, dbname='SCEUA', dbformat=dbformat,
+                                            sim_timeout=timeout, alt_objfun=None)
+    sceua_sampler.sample(rep)
+    results.append(sceua_sampler.getdata())
+    sceua_elapsed = time.time() - start
+
+
+    print("#########################################")
+
+    #print(dds_elapsed, dds_sampler.status.params)
+
+    print(sceua_elapsed, sceua_sampler.status.params)
+
+    benchmarks_dict.append({
+        "rep": rep,
+        "dds_time": dds_elapsed,
+        "sceua_time": sceua_elapsed,
+        "dds_like": dds_sampler.status.objectivefunction,
+        "sceua_like": sceua_sampler.status.objectivefunction,
+        "dds_param": list(dds_sampler.status.params),
+        "sceua_param": list(sceua_sampler.status.params)
+    })
+    benchmarks_duration["dds"].append(dds_elapsed)
+    benchmarks_duration["sceua"].append(sceua_elapsed)
+    benchmarks_duration["sceua_like"].append(sceua_sampler.status.objectivefunction)
+    benchmarks_duration["dds_like"].append(dds_sampler.status.objectivefunction)
+
+print(json.dumps(benchmarks_dict))
+
+
+def autolabel(rects):
+    """
+    Attach a text label above each bar displaying its height
+    """
+    for rect in rects:
+        height = rect.get_height()
+        ax.text(rect.get_x() + rect.get_width()/2., 1.05*height,
+                '%d' % int(height),
+                ha='center', va='bottom')
+
+
+
+fig = plt.figure(figsize=(10, 6))
+ax = plt.subplot(111)
+
+rep_labels = [str(j) for j in reps]
+x_pos = [i for i, _ in enumerate(rep_labels)]
+
+
+X = np.arange(len(benchmarks_duration["dds"]))
+dds_plot = ax.bar(x_pos, benchmarks_duration["dds_like"], color = 'b', width = 0.45)
+sceua_plot = ax.bar([j+0.45 for j in x_pos], benchmarks_duration["sceua_like"], color = 'g', width = 0.45)
+
+#dds_plot = ax.bar(x_pos, benchmarks_duration["dds"], color = 'b', width = 0.45)
+#sceua_plot = ax.bar([j+0.45 for j in x_pos], benchmarks_duration["sceua"], color = 'g', width = 0.45)
+
+
+
+plt.xticks(x_pos, rep_labels)
+plt.legend(("DDS", "SCEUA"))
+plt.xlabel("Repetitions")
+plt.ylabel("Best Objective Function Value")
+
+autolabel(dds_plot)
+autolabel(sceua_plot)
+
+plt.show()
+plt.savefig("MPI_TEST")
+#
+

spotpy/examples/dds/dds_parallel.py

@@ -0,0 +1,46 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+import numpy as np
+import sys
+import os
+import matplotlib.pylab as plt
+import json
+import time
+
+try:
+    import spotpy
+except ImportError:
+    import sys
+
+    sys.path.append(".")
+    import spotpy
+
+from spotpy.examples.spot_setup_hymod_python import spot_setup
+
+spot_setup = spot_setup()
+
+
+path = os.path.abspath(os.path.dirname(__file__))
+json_path = path + "/dds_parallel_data.json"
+benchmarks_duration = json.load(open(json_path))
+
+
+rep = int(sys.argv[1])
+timeout = 10  # Given in Seconds
+parallel = "mpi"
+dbformat = "csv"
+start = time.time()
+dds_sampler = spotpy.algorithms.dds(spot_setup, parallel=parallel, dbname='DDS', dbformat=dbformat, sim_timeout=timeout)
+dds_sampler.sample(rep, trials=1)
+dds_elapsed = time.time() - start
+print(dds_elapsed)
+
+benchmarks_duration["dds_duration"].append(dds_elapsed)
+benchmarks_duration["dds_like"].append(dds_sampler.status.objectivefunction)
+benchmarks_duration["rep"].append(rep)
+
+print(benchmarks_duration)
+
+json.dump(benchmarks_duration, open(json_path,"w"))

spotpy/examples/dds/dds_parallel_plot.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+import matplotlib.pylab as plt
+import json
+import matplotlib as mp
+
+data_normalizer = mp.colors.Normalize()
+color_map = mp.colors.LinearSegmentedColormap(
+    "my_map",
+    {
+        "red": [(0, 1.0, 1.0),
+                (1.0, .5, .5)],
+        "green": [(0, 0.5, 0.5),
+                  (1.0, 0, 0)],
+        "blue": [(0, 0.50, 0.5),
+                 (1.0, 0, 0)]
+    }
+)
+
+
+def autolabel(ax, rects):
+    """
+    Attach a text label above each bar displaying its height
+    """
+    for rect in rects:
+        height = rect.get_height()
+        ax.text(rect.get_x() + rect.get_width() / 2., 1.05 * height,
+                '%f' % height,
+                ha='center', va='bottom')
+
+
+def subplot(data, name, ylabel):
+    fig = plt.figure(figsize=(20, 6))
+    ax = plt.subplot(111)
+    rep_labels = [str(j) for j in reps]
+    x_pos = [i for i, _ in enumerate(rep_labels)]
+    X = np.arange(len(data))
+    ax_plot = ax.bar(x_pos, data, color=color_map(data_normalizer(data)), width=0.45)
+
+    plt.xticks(x_pos, rep_labels)
+    plt.xlabel("Repetitions")
+    plt.ylabel(ylabel)
+
+    autolabel(ax, ax_plot)
+    plt.savefig(name + ".png")
+
+
+parallel_data = json.loads('{"dds_duration": [1.1293659210205078, 3.254117250442505, 4.888171672821045, 18.719818592071533, 34.56907820701599, 169.47716689109802, 337.86882615089417, 1644.955144405365, 3348.948029756546], "rep": [30, 500, 1000, 5000, 10000, 50000, 100000, 500000, 1000000], "dds_like": [-8384.884435178812, -8269.480874403698, -8268.453892284442, -8268.51195094138, -8269.65509041187, -8268.1421690868, -8267.791798085422, -8267.79178644684, -8268.141980514703]}')
+reps = parallel_data["rep"]
+subplot(parallel_data["dds_duration"], "DDS_PARALLEL_DURATION_all", "Duration of Run in Seconds")
+subplot(parallel_data["dds_like"], "DDS_PARALLEL_OBJECTIVEFUNCTION_all", "Best Objective Function Value")

spotpy/examples/dds/dds_parallel_run.py

@@ -0,0 +1,8 @@
+import subprocess
+import os
+path = os.path.abspath(os.path.dirname(__file__))
+for r in [500,1000,5000,10000,50000,100000,500000]:
+    args = ["mpirun", "-c 6", "python", path + "/dds_parallel.py", str(r)]
+    print(args)
+    subprocess.run(args)
+    exit(8)

spotpy/examples/spot_setup_dds.py

@@ -0,0 +1,83 @@
+import spotpy
+from spotpy.parameter import Uniform
+from spotpy.objectivefunctions import rmse
+import numpy as np
+
+
+def ackley10(vector):
+    length = len(vector)
+    sum1 = 0
+    sum2 = 0
+    for i in range(length):
+        sum1 = sum1 + vector[i] ** 2
+        sum2 = sum2 + np.cos(2 * np.pi * vector[i])
+    return -1*(-20 * np.exp(-0.2 * (sum1 / length) ** 0.5) - np.exp(sum2 / length))
+
+
+def griewank10(vector):
+    sum1 = 0
+    term2 = 1
+    term3 = 1
+
+    for i in range(len(vector)):
+        sum1 = sum1 + (vector[i] ** 2) / 4000
+        term2 = term2 * np.cos(vector[i] / (i + 1) ** 0.5)
+
+    return -1*(sum1 - term2 + term3)
+
+
+class spot_setup(object):
+    """
+        Setup for a simple example to run DDS Algorithm
+    """
+
+    def __init__(self):
+        self.params = None
+        self.objfunc = None
+
+    def _objfunc_switcher(self, name):
+        """
+        Set new parameter and objective function while setup is instanced in a test case
+        :param name: function name which overwrites initial objective function
+        :return:
+        """
+
+        if name == "ackley":
+            self.objfunc = ackley10
+            self.params = [Uniform(str(j), -2, 2, 1.5, 3.0, -2, 2, doc=str(j) + ' value of Rosenbrock function')
+                           for j in range(10)]
+        elif name == "griewank":
+            self.objfunc = griewank10
+            self.params = [Uniform('d' + str(j), -500, 700, 1.5, 3.0, -500, 700,
+                                   doc=str(j) + 'distinc parameter within a boundary', as_int=True)
+                           for j in range(2)] + [Uniform('c' + str(j), -500, 700, 1.5, 3.0, -500, 700,
+                                                         doc=str(j) + 'continuous parameter within a boundary')
+                                                 for j in range(8)]
+
+    def parameters(self):
+        if self.params is None:
+            self.params = [
+                Uniform("0", -10, 10, 1.5, 3.0, -10, 10, doc='x value of Rosenbrock function'),
+                Uniform("1", -10, 10, 1.5, 3.0, -10, 10, doc='y value of Rosenbrock function'),
+                Uniform("z", -10, 10, 1.5, 3.0, -10, 10, doc='z value of Rosenbrock function')]
+        return spotpy.parameter.generate(self.params)
+
+    def simulation(self, vector):
+        x = np.array(vector)
+        # simulations = [sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)]
+        simulations = x * np.random.rand(len(vector))
+        #simulations = x * np.sum(vector)
+        return simulations
+
+    def evaluation(self):
+        # observations = [0]
+        observations = [2, 3, 4]
+        return observations
+
+    def objectivefunction(self, simulation, evaluation, params):
+
+        if self.objfunc is None:
+            return -1*rmse(evaluation, simulation)
+        else:
+            pars, names = params
+            return self.objfunc(pars)