Merge branch 'fix_disp_and_tests' into 'master'

Part 1 of an update to include example_3_turbine_type_optimization.py

See merge request !61

topfarm/_topfarm.py

@@ -16,17 +16,31 @@ from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp
 
 
 class TopFarmProblem(Problem):
-    def __init__(self, cost_comp, record_id, expected_cost):
+    def __init__(self, cost_comp, driver, plot_comp, record_id, expected_cost):
         Problem.__init__(self)
         if isinstance(cost_comp, TopFarmProblem):
             cost_comp = cost_comp.as_component()
         cost_comp.parent = self
-        self.record_id = record_id
         self.cost_comp = cost_comp
+        
+        if isinstance(driver, list):
+            driver = DOEDriver(ListGenerator(driver))
+        elif isinstance(driver, DOEGenerator):
+            driver = DOEDriver(generator=driver)
+        self.driver = driver
+
+        self.plot_comp = plot_comp
+
+        self.record_id = record_id
+
         self.indeps = self.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
         self.model.add_subsystem('cost_comp', cost_comp, promotes=['*'])
         self.model.add_objective('cost', scaler=1 / abs(expected_cost))
 
+        if plot_comp:
+            self.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
+            plot_comp.problem = self
+
     @property
     def cost(self):
         return self['cost'][0]
@@ -55,10 +69,10 @@ class TopFarmProblem(Problem):
     def evaluate(self, state={}):
         t = time.time()
         self.update_state(state)
-        rec = ListRecorder()
-        self.driver.add_recorder(rec)
+        self.recorder = ListRecorder()
+        self.driver.add_recorder(self.recorder)
         self.run_model()
-        self.driver._rec_mgr._recorders.remove(rec)
+        self.driver._rec_mgr._recorders.remove(self.recorder)
         print("Evaluated in\t%.3fs" % (time.time() - t))
         return self.cost, self.state
 
@@ -81,8 +95,8 @@ class TopFarmProblem(Problem):
             try:
                 self.update_state({k: self.recorder.get(k)[-1] for k in self.state.keys() if k not in state})
             except ValueError:
-                pass # loaded state does not fit into dimension of current state
-            
+                pass  # loaded state does not fit into dimension of current state
+
         self.driver.add_recorder(self.recorder)
         self.run_driver()
         self.cleanup()
@@ -127,8 +141,8 @@ class TopFarmProblem(Problem):
 
 
 class TurbineTypeOptimizationProblem(TopFarmProblem):
-    def __init__(self, cost_comp, turbineTypes, lower, upper, driver, record_id=None, expected_cost=1):
-        TopFarmProblem.__init__(self, cost_comp, record_id, expected_cost)
+    def __init__(self, cost_comp, turbineTypes, lower, upper, driver, plot_comp=None, record_id=None, expected_cost=1):
+        TopFarmProblem.__init__(self, cost_comp, driver, plot_comp, record_id, expected_cost)
         self.turbineTypes = turbineTypes
         self.lower = lower
         self.upper = upper
@@ -142,13 +156,15 @@ class TurbineTypeOptimizationProblem(TopFarmProblem):
 
         self.model.add_constraint('turbineType', lower=lim[:, 0], upper=lim[:, 1])
         self.indeps.add_output('turbineType', np.array(turbineTypes).astype(np.int))
-        self.driver = driver
+
         self.model.add_design_var('turbineType', lower=lim[:, 0], upper=lim[:, 1])
+        if self.plot_comp:
+            plot_comp.n_wt = n_wt
         self.setup(check=True, mode='fwd')
 
     @property
     def state(self):
-        state = {'turbineType': np.array(self['turbineType']).astype(np.int)}
+        state = {'turbineType': np.round(self['turbineType']).astype(np.int)}
         state.update(TopFarmProblem.state.fget(self))
         return state
 
@@ -156,11 +172,16 @@ class TurbineTypeOptimizationProblem(TopFarmProblem):
 class TurbineXYZOptimizationProblem(TopFarmProblem):
     def __init__(self, cost_comp, turbineXYZ, boundary_comp, min_spacing=None,
                  driver=ScipyOptimizeDriver(), plot_comp=None, record_id=None, expected_cost=1):
-        self.turbineXYZ = turbineXYZ
-        TopFarmProblem.__init__(self, cost_comp, record_id, expected_cost)
-
         turbineXYZ = np.array(turbineXYZ)
+        self.turbineXYZ = turbineXYZ
         self.n_wt = n_wt = turbineXYZ.shape[0]
+
+        if plot_comp:
+            if plot_comp == "default":
+                plot_comp = PlotComp()
+
+        TopFarmProblem.__init__(self, cost_comp, driver, plot_comp, record_id, expected_cost)
+
         turbineXYZ = np.hstack((turbineXYZ, np.zeros((n_wt, 4 - turbineXYZ.shape[1]))))
         self.min_spacing = min_spacing
 
@@ -170,7 +191,7 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
         self.boundary_comp = boundary_comp
         boundary_comp.setup_as_constraints(self)
 
-        do = driver.options
+        do = self.driver.options
         if len(boundary_comp.xy_boundary) > 0:
 
             ref0_x, ref0_y = self.boundary_comp.xy_boundary.min(0)
@@ -203,20 +224,13 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
         self.indeps.add_output('turbineY', turbineXYZ[:, 1], units='m')
         self.indeps.add_output('turbineZ', turbineXYZ[:, 2], units='m')
 
-        self.driver = driver
-
         if plot_comp:
-            if plot_comp == "default":
-                plot_comp = PlotComp()
-
             plot_comp.n_wt = n_wt
             if self.boundary_comp:
                 plot_comp.n_vertices = len(self.boundary_comp.xy_boundary)
             else:
                 plot_comp.n_vertices = 0
-            self.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
 
-        self.plot_comp = plot_comp
         self.setup(check=True, mode='fwd')
 
     @property
@@ -252,10 +266,6 @@ class InitialXYZOptimizationProblem(TurbineXYZOptimizationProblem):
                  driver=None, plot_comp=None):
         #          if driver is None:
         #              driver = DOEDriver(shuffle_generator(self, 10))
-        if isinstance(driver, list):
-            driver = DOEDriver(ListGenerator(driver))
-        elif isinstance(driver, DOEGenerator):
-            driver = DOEDriver(generator=driver)
         TurbineXYZOptimizationProblem.__init__(self, cost_comp, turbineXYZ,
                                                boundary_comp, min_spacing,
                                                driver=driver, plot_comp=plot_comp)

topfarm/constraint_components/constrained_generator.py

@@ -2,6 +2,30 @@ from openmdao.drivers.doe_generators import DOEGenerator, ListGenerator
 import numpy as np
 
 
+
+class ConstrainedDiscardXYZGenerator(DOEGenerator):
+    def __init__(self, generator):
+        DOEGenerator.__init__(self)
+        if isinstance(generator, list):
+            generator = ListGenerator(generator)
+        self.generator = generator
+        
+    def __call__(self, design_vars, model=None):
+        xy_boundary_comp = model._get_subsystem('xy_bound_comp')
+        spacing_comp = model._get_subsystem('spacing_comp')
+        for xyz_tuple_lst in self.generator(design_vars, model=model):
+            x, y, z = ([np.array(t[1]) for t in xyz_tuple_lst] + [0])[:3]
+            if spacing_comp:
+                dist = spacing_comp._compute(x, y)
+                if np.any(dist < spacing_comp.min_spacing**2):
+                    continue
+                 
+            if xy_boundary_comp:
+                dist_to_boundary =xy_boundary_comp.distances(x,y)
+                if np.any(dist_to_boundary<0):
+                    continue
+            yield xyz_tuple_lst
+
 class ConstrainedXYZGenerator(DOEGenerator):
     def __init__(self, generator, n_iter=1000, step_size=0.1, offset=0.5, verbose=False):
         DOEGenerator.__init__(self)
@@ -17,16 +41,15 @@ class ConstrainedXYZGenerator(DOEGenerator):
         xy_boundary_comp = model._get_subsystem('xy_bound_comp')
         spacing_comp = model._get_subsystem('spacing_comp')
         for xyz_tuple_lst in self.generator(design_vars, model=model):
-            x, y, z = [np.array(t[1]) for t in xyz_tuple_lst]
+            x, y, z = ([np.array(t[1]).astype(np.float) for t in xyz_tuple_lst] + [0])[:3]
             if spacing_comp:
                 for j in range(self.n_iter):
                     dist = spacing_comp._compute(x, y)
                     dx, dy = spacing_comp._compute_partials(x, y)
-                    dx, dy = dx[0], dy[0]
-                    index = np.argmin(dist)
+                    index = int(np.argmin(dist))
                     if dist[index] < spacing_comp.min_spacing**2 or j == 0:
-                        x[index] += dx[index] * self.step_size
-                        y[index] += dy[index] * self.step_size
+                        x += dx[index] * self.step_size
+                        y += dy[index] * self.step_size
                         if xy_boundary_comp:
                             x, y, z = xy_boundary_comp.move_inside(x, y, z)
                     else:

topfarm/plotting.py

@@ -21,15 +21,19 @@ def mypause(interval):
 class PlotComp(ExplicitComponent):
     colors = ['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 100
 
-    def __init__(self, memory=10, delay=0.001, plot_initial=True,
-                 animate=False):
+    def __init__(self, memory=10, delay=0.001, plot_initial=True, ax=None):
         ExplicitComponent.__init__(self)
+        self.ax_ = ax
         self.memory = memory
         self.delay = delay
         self.plot_initial = plot_initial
         self.history = []
         self.counter = 0
-        self.animate = animate
+        self.by_pass = False
+
+    @property
+    def ax(self):
+        return self.ax_ or plt.gca()
 
     def show(self):
         plt.show()
@@ -38,66 +42,92 @@ class PlotComp(ExplicitComponent):
         self.add_input('turbineX', np.zeros(self.n_wt), units='m')
         self.add_input('turbineY', np.zeros(self.n_wt), units='m')
         self.add_input('cost', 0.)
-        self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
-
-    def init_plot(self, boundary):
-        plt.cla()
-        plt.axis('equal')
-        b = np.r_[boundary[:], boundary[:1]]
-        plt.plot(b[:, 0], b[:, 1], 'k')
-        mi = b.min(0)
-        ma = b.max(0)
+        if hasattr(self, 'n_vertices'):
+            self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
+
+    def init_plot(self, limits):
+        self.ax.cla()
+        self.ax.axis('equal')
+        mi = limits.min(0)
+        ma = limits.max(0)
         ra = ma - mi
         ext = .1
         xlim, ylim = np.array([mi - ext * ra, ma + ext * ra]).T
-        plt.xlim(xlim)
-        plt.ylim(ylim)
+        self.ax.set_xlim(xlim)
+        self.ax.set_ylim(ylim)
+
+    def plot_boundary(self, boundary):
+        b = np.r_[boundary[:], boundary[:1]]
+        plt.plot(b[:, 0], b[:, 1], 'k')
+
+    def plot_history(self, x, y):
+        rec = self.problem.recorder
+        if rec.num_cases > 0:
+            x = np.r_[rec['turbineX'][-self.memory:], [x]]
+            y = np.r_[rec['turbineY'][-self.memory:], [y]]
+            for c, x_, y_ in zip(self.colors, x.T, y.T):
+                self.ax.plot(x_, y_, '--', color=c)
+
+    def plot_initial2current(self, x0, y0, x, y):
+        rec = self.problem.recorder
+        if rec.num_cases > 0:
+            x0 = rec['turbineX'][0]
+            y0 = rec['turbineY'][0]
+            for c, x0_, y0_, x_, y_ in zip(self.colors, x0, y0, x, y):
+                self.ax.plot((x0_, x_), (y0_, y_), '-', color=c)
+
+    def plot_current_position(self, x, y):
+        for c, x_, y_ in zip(self.colors, x, y):
+            self.ax.plot(x_, y_, 'o', color=c, ms=5)
+            self.ax.plot(x_, y_, 'xk', ms=4)
+
+    def set_title(self, cost0, cost):
+        rec = self.problem.recorder
+        if cost0 != 0:
+            self.ax.set_title("%d: %f (%.2f%%)" % (rec.num_cases, cost, (cost0 - cost) / cost0 * 100))
+        else:
+            self.ax.set_title("%d: %f" % (rec.num_cases, cost))
+
+    def get_initial(self):
+        rec = self.problem.recorder
+        if rec.num_cases > 0:
+            x0 = rec['turbineX'][0]
+            y0 = rec['turbineY'][0]
+            cost0 = rec['cost'][0]
+            return x0, y0, cost0
 
     def compute(self, inputs, outputs):
-        x = inputs['turbineX']
-        y = inputs['turbineY']
-        cost = inputs['cost'][0]
-        if not hasattr(self, "initial"):
-            self.initial = np.array([x, y]).T, cost
-        self.history = [(x.copy(), y.copy())] + self.history[:self.memory]
-
-        boundary = inputs['boundary']
-        self.init_plot(boundary)
-        plt.title("%d: %f (%.2f%%)" % (self.counter, cost,
-                  (self.initial[1]-cost)/self.initial[1]*100))
-
-        history_arr = np.array(self.history)
-        for i, c, x_, y_ in zip(range(len(x)), self.colors, x, y):
-            if self.plot_initial:
-                plt.plot([self.initial[0][i, 0], x_],
-                         [self.initial[0][i, 1], y_], '-', color=c, lw=1)
-            plt.plot(history_arr[:, 0, i], history_arr[:, 1, i], '.--',
-                     color=c, lw=1)
-            plt.plot(x_, y_, 'o', color=c, ms=5)
-            plt.plot(x_, y_, 'x' + 'k', ms=4)
-            self.temp = '../__animations__'
-            if not os.path.exists(self.temp):
-                os.makedirs(self.temp)
-            if self.animate:
-                path = os.path.join(self.temp,
-                                    'plot_{:05d}.png'.format(self.counter))
-                plt.savefig(path)
-
-        if self.counter == 0:
-            plt.pause(.01)
-        mypause(self.delay)
-
-        self.counter += 1
-
-    def run_animate(self, anim_time=10, verbose=False):
-        N = anim_time/self.counter
-        string = 'ffmpeg -f image2 -r 1/'
-        string += '{} -i {}//plot_%05d.png'.format(N, self.temp)
-        string += ' -vcodec mpeg4 -y {}//animation.mp4'.format(self.temp)
-        if verbose:
-            print('\nCreating animation:')
-            print(string)
-        os.system(string)
+        if self.by_pass is False:
+            x = inputs['turbineX']
+            y = inputs['turbineY']
+            cost = inputs['cost'][0]
+
+            if 'boundary' in inputs:
+                boundary = inputs['boundary']
+                self.init_plot(boundary)
+                self.plot_boundary(boundary)
+            else:
+                self.init_plot(np.array([x, y]).T)
+
+            initial = self.get_initial()
+
+            if initial is not None:
+                x0, y0, cost0 = initial
+                if self.plot_initial:
+                    self.plot_initial2current(x0, y0, x, y)
+                if self.memory > 0:
+                    self.plot_history(x, y)
+            else:
+                cost0 = cost
+
+            self.plot_current_position(x, y)
+            self.set_title(cost0, cost)
+
+            if self.counter == 0:
+                plt.pause(.01)
+            mypause(self.delay)
+
+            self.counter += 1
 
 
 class NoPlot(PlotComp):
@@ -108,10 +138,35 @@ class NoPlot(PlotComp):
         pass
 
     def setup(self):
-        self.add_input('turbineX', np.zeros(self.n_wt), units='m')
-        self.add_input('turbineY', np.zeros(self.n_wt), units='m')
         self.add_input('cost', 0.)
-        self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
 
     def compute(self, inputs, outputs):
         pass
+
+
+class TurbineTypePlotComponent(PlotComp):
+    colors = ['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 100
+    markers = np.array(list(".ov^<>12348spP*hH+xXDd|_"))
+
+    def __init__(self, turbine_type_names, memory=10, delay=0.001, plot_initial=True):
+        self.turbine_type_names = turbine_type_names
+        PlotComp.__init__(self, memory=memory, delay=delay, plot_initial=plot_initial)
+
+    def setup(self):
+        PlotComp.setup(self)
+        self.add_input('turbineType', np.zeros(self.n_wt, dtype=np.int))
+
+    def compute(self, inputs, outputs):
+        self.types = np.asarray(inputs['turbineType'], dtype=np.int)
+        PlotComp.compute(self, inputs, outputs)
+
+    def init_plot(self, limits):
+        PlotComp.init_plot(self, limits)
+        for m, n in zip(self.markers, self.turbine_type_names):
+            self.ax.plot([], [], m + 'k', label=n)
+        self.ax.legend()
+
+    def plot_current_position(self, x, y):
+        for m, c, x_, y_ in zip(self.markers[self.types], self.colors, x, y):
+            #self.ax.plot(x_, y_, 'o', color=c, ms=5)
+            self.ax.plot(x_, y_,  m + 'k', markeredgecolor=c, markeredgewidth=1, ms=8)

topfarm/recorders.py

@@ -84,6 +84,10 @@ class ListRecorder(BaseRecorder):
             self._abs2meta[name]['explicit'] = True
             if name in states:
                 self._abs2meta[name]['explicit'] = False
+                
+    @property
+    def num_cases(self):
+        return len(self.driver_iteration_lst)
 
     @property
     def driver_cases(self):

topfarm/tests/test_topfarm_utils/test_parellel_runner.py

 from openmdao.drivers.doe_driver import DOEDriver
-from openmdao.drivers.doe_generators import FullFactorialGenerator,\
-    ListGenerator, UniformGenerator
+from openmdao.drivers.doe_generators import UniformGenerator
 import pytest
 from topfarm._topfarm import InitialXYZOptimizationProblem
 import numpy as np
 from topfarm.cost_models.dummy import DummyCost
-from topfarm.tests import npt, uta
-from topfarm.constraint_components.constrained_generator import ConstrainedXYZGenerator
+from topfarm.tests import npt
 from topfarm.constraint_components.boundary_component import BoundaryComp
 from topfarm.parallel_runner import ParallelRunner
 

topfarm/tests/test_with_dummy.py

@@ -6,52 +6,49 @@ import unittest
 import pytest
 import numpy as np
 from topfarm.cost_models.dummy import DummyCost, DummyCostPlotComp
-
-
-class Test(unittest.TestCase):  # unittest version
-
-    def test_optimize_4tb(self):
-        """Optimize 4-turbine layout and check final positions
-
-            The farm boundaries and min spacing are chosen such that the desired
-            turbine positions are not within the boundaries or constraints.
-        """
-
-        # test options
-        dec_prec = 4  # decimal precision for comparison
-
-        # given
-        boundary = [(0, 0), (6, 0), (6, -10), (0, -10)]  # turbine boundaries
-        initial = [[6, 0], [6, -8], [1, 1], [-1, -8]]  # initial turbine layouts
-        desired = [[3, -3], [7, -7], [4, -3], [3, -7]]  # desired turbine layouts
-        optimal = np.array([[2.5, -3], [6, -7], [4.5, -3], [3, -7]])  # optimal turbine layout
-        min_spacing = 2  # min distance between turbines
-
-        # when
-        tf = TopFarm(initial, DummyCost(desired, ['turbineX','turbineY']), min_spacing,
-                     boundary=boundary, record_id=None)
-        tf.optimize()
-        tb_pos = tf.turbine_positions
-
-        # then
-        tol = 1e-6
-        self.assertGreater(sum((tb_pos[2] - tb_pos[0])**2), 2**2 - tol)  # check min spacing
-        self.assertLess(tb_pos[1][0], 6 + tol)  # check within border
-        np.testing.assert_array_almost_equal(tb_pos[:,:2], optimal, dec_prec)
-
-    def testDummyCostPlotComp(self):
-        if os.name == 'posix' and "DISPLAY" not in os.environ:
-            pytest.xfail("No display")
-        
-        desired = [[3, -3], [7, -7], [4, -3], [3, -7]]
-        tf = TopFarm(turbines=[[6, 0], [6, -8], [1, 1], [-1, -8]],
-                     cost_comp=DummyCost(desired, ['turbineX','turbineY']),
-                     min_spacing=2,
-                     boundary=[(0, 0), (6, 0), (6, -10), (0, -10)],
-                     plot_comp = DummyCostPlotComp(desired),
-                     record_id=None)
-        tf.evaluate()
-
-
-if __name__ == "__main__":
-    unittest.main()
+from topfarm.tests import uta
+
+
+def test_optimize_4tb():
+    """Optimize 4-turbine layout and check final positions
+
+        The farm boundaries and min spacing are chosen such that the desired
+        turbine positions are not within the boundaries or constraints.
+    """
+
+    # test options
+    dec_prec = 4  # decimal precision for comparison
+
+    # given
+    boundary = [(0, 0), (6, 0), (6, -10), (0, -10)]  # turbine boundaries
+    initial = [[6, 0], [6, -8], [1, 1], [-1, -8]]  # initial turbine layouts
+    desired = [[3, -3], [7, -7], [4, -3], [3, -7]]  # desired turbine layouts
+    optimal = np.array([[2.5, -3], [6, -7], [4.5, -3], [3, -7]])  # optimal turbine layout
+    min_spacing = 2  # min distance between turbines
+
+    # when
+    tf = TopFarm(initial, DummyCost(desired, ['turbineX', 'turbineY']), min_spacing,
+                 boundary=boundary, record_id=None)
+    tf.optimize()
+    tb_pos = tf.turbine_positions
+
+    # then
+    tol = 1e-6
+    uta.assertGreater(sum((tb_pos[2] - tb_pos[0])**2), 2**2 - tol)  # check min spacing
+    uta.assertLess(tb_pos[1][0], 6 + tol)  # check within border
+    np.testing.assert_array_almost_equal(tb_pos[:, :2], optimal, dec_prec)
+
+
+def testDummyCostPlotComp():
+    if os.name == 'posix' and "DISPLAY" not in os.environ:
+        pytest.xfail("No display")
+
+    desired = [[3, -3], [7, -7], [4, -3], [3, -7]]
+    tf = TopFarm(turbines=[[6, 0], [6, -8], [1, 1], [-1, -8]],
+                 cost_comp=DummyCost(desired, ['turbineX', 'turbineY']),
+                 min_spacing=2,
+                 boundary=[(0, 0), (6, 0), (6, -10), (0, -10)],
+                 plot_comp=DummyCostPlotComp(desired),
+                 record_id=None)
+    tf.evaluate()
+    tf.optimize()