Merge branch '7-implement-polygon-boundary' into 'master'

Resolve "implement polygon boundary"

Closes #7

See merge request !11

topfarm/constraint_components/boundary_component.py

@@ -6,15 +6,14 @@ from scipy.spatial import ConvexHull
 class BoundaryComp(ExplicitComponent):
 
     def __init__(self, vertices, nTurbines, boundary_type='convex_hull'):
-
         super(BoundaryComp, self).__init__()
-
-        self.calculate_boundary_and_normals(vertices, boundary_type)
         self.nTurbines = nTurbines
+        self.vertices = np.array(vertices)
         self.nVertices = self.vertices.shape[0]
+        self.calculate_boundary_and_normals(vertices, boundary_type)
+        self.calculate_gradients()
 
     def calculate_boundary_and_normals(self, vertices, boundary_type):
-
         if boundary_type == 'convex_hull':
             # find the points that actually comprise a convex hull
             hull = ConvexHull(list(vertices))
@@ -59,6 +58,33 @@ class BoundaryComp(ExplicitComponent):
 
         self.vertices, self.unit_normals = vertices, unit_normals
 
+    def calculate_gradients(self):
+        unit_normals = self.unit_normals
+
+        # initialize array to hold distances from each point to each face
+        dfaceDistance_dx = np.zeros([self.nTurbines * self.nVertices, self.nTurbines])
+        dfaceDistance_dy = np.zeros([self.nTurbines * self.nVertices, self.nTurbines])
+
+        for i in range(0, self.nTurbines):
+            # determine if point is inside or outside of each face, and distance from each face
+            for j in range(0, self.nVertices):
+
+                # define the derivative vectors from the point of interest to the first point of the face
+                dpa_dx = np.array([-1.0, 0.0])
+                dpa_dy = np.array([0.0, -1.0])
+
+                # find perpendicular distance derivatives from point to current surface (vector projection)
+                ddistanceVec_dx = np.vdot(dpa_dx, unit_normals[j]) * unit_normals[j]
+                ddistanceVec_dy = np.vdot(dpa_dy, unit_normals[j]) * unit_normals[j]
+
+                # calculate derivatives for the sign of perpendicular distance from point to current face
+                dfaceDistance_dx[i * self.nVertices + j, i] = np.vdot(ddistanceVec_dx, unit_normals[j])
+                dfaceDistance_dy[i * self.nVertices + j, i] = np.vdot(ddistanceVec_dy, unit_normals[j])
+
+        # return Jacobian dict
+        self.dfaceDistance_dx = dfaceDistance_dx
+        self.dfaceDistance_dy = dfaceDistance_dy
+
     def calculate_distance_to_boundary(self, points):
         """
         :param points: points that you want to calculate the distance from to the faces of the convex hull
@@ -121,29 +147,151 @@ class BoundaryComp(ExplicitComponent):
         outputs['boundaryDistances'] = self.calculate_distance_to_boundary(locations)
 
     def compute_partials(self, inputs, partials):
+        # return Jacobian dict
+        partials['boundaryDistances', 'turbineX'] = self.dfaceDistance_dx
+        partials['boundaryDistances', 'turbineY'] = self.dfaceDistance_dy
 
-        unit_normals = self.unit_normals
 
-        # initialize array to hold distances from each point to each face
-        dfaceDistance_dx = np.zeros([self.nTurbines * self.nVertices, self.nTurbines])
-        dfaceDistance_dy = np.zeros([self.nTurbines * self.nVertices, self.nTurbines])
+class PolygonBoundaryComp(BoundaryComp):
 
-        for i in range(0, self.nTurbines):
-            # determine if point is inside or outside of each face, and distance from each face
-            for j in range(0, self.nVertices):
+    def __init__(self, vertices, nTurbines):
 
-                # define the derivative vectors from the point of interest to the first point of the face
-                dpa_dx = np.array([-1.0, 0.0])
-                dpa_dy = np.array([0.0, -1.0])
+        super(BoundaryComp, self).__init__()
 
-                # find perpendicular distance derivatives from point to current surface (vector projection)
-                ddistanceVec_dx = np.vdot(dpa_dx, unit_normals[j]) * unit_normals[j]
-                ddistanceVec_dy = np.vdot(dpa_dy, unit_normals[j]) * unit_normals[j]
+        self.nTurbines = nTurbines
+        vertices = np.array(vertices)
+        self.nVertices = vertices.shape[0]
+
+        def edges_counter_clockwise(vertices):
+            if np.any(vertices[0] != vertices[-1]):
+                vertices = np.r_[vertices, vertices[:1]]
+
+            x1, y1 = vertices[:-1].T
+            x2, y2 = vertices[1:].T
+            double_area = np.sum((x1 - x2) * (y1 + y2))  # 2 x Area (+: counterclockwise
+            assert double_area != 0, "Area must be non-zero"
+            if double_area < 0:  #
+                return edges_counter_clockwise(vertices[::-1])
+            else:
+                return vertices[:-1], x1, y1, x2, y2
+
+        self.vertices, self.x1, self.y1, self.x2, self.y2 = edges_counter_clockwise(vertices)
+        self.min_x, self.min_y = np.min([self.x1, self.x2], 0), np.min([self.y1, self.y2], )
+        self.max_x, self.max_y = np.max([self.x1, self.x2], 1), np.max([self.y1, self.y2], 0)
+        self.dx = self.x2 - self.x1
+        self.dy = self.y2 - self.y1
+        self.x2y1 = self.x2 * self.y1
+        self.y2x1 = self.y2 * self.x1
+        self.length = ((self.y2 - self.y1)**2 + (self.x2 - self.x1)**2)**0.5
+        self.edge_unit_vec = (np.array([self.dy, -self.dx]) / self.length)
+        v = np.hstack((self.edge_unit_vec, self.edge_unit_vec[:, :1]))
+        self.xy2_vec = v[:, :-1] + v[:, 1:]
+        self.xy1_vec = np.hstack((self.xy2_vec[:, -1:], self.xy2_vec[:, 1:]))
+
+        self.dEdgeDist_dx = -self.dy / self.length
+        self.dEdgeDist_dy = self.dx / self.length
+        self._cache_input = None
+        self._cache_output = None
 
-                # calculate derivatives for the sign of perpendicular distance from point to current face
-                dfaceDistance_dx[i * self.nVertices + j, i] = np.vdot(ddistanceVec_dx, unit_normals[j])
-                dfaceDistance_dy[i * self.nVertices + j, i] = np.vdot(ddistanceVec_dy, unit_normals[j])
+    def setup(self):
 
-        # return Jacobian dict
-        partials['boundaryDistances', 'turbineX'] = dfaceDistance_dx
-        partials['boundaryDistances', 'turbineY'] = dfaceDistance_dy
+        # Explicitly size input arrays
+        self.add_input('turbineX', np.zeros(self.nTurbines), units='m',
+                       desc='x coordinates of turbines in global ref. frame')
+        self.add_input('turbineY', np.zeros(self.nTurbines), units='m',
+                       desc='y coordinates of turbines in global ref. frame')
+
+        # Explicitly size output array
+        # (vector with positive elements if turbines outside of hull)
+        self.add_output('boundaryDistances', np.zeros(self.nTurbines),
+                        desc="signed perpendicular distance from each turbine to each face CCW; + is inside")
+
+        self.declare_partials('boundaryDistances', ['turbineX', 'turbineY'])
+        #self.declare_partials('boundaryDistances', ['boundaryVertices', 'boundaryNormals'], method='fd')
+
+    def calc_distance_and_gradients(self, x, y):
+        """
+        distance point(x,y) to edge((x1,y1)->(x2,y2))
+        +/-: inside/outside 
+        case (x,y) closest to edge: 
+            distance = edge_unit_vec dot (x1-x,y1-y)
+            ddist_dx = -(y2-y2)/|edge|
+            ddist_dy = (x2-x2)/|edge|
+        case (x,y) closest to (x1,y1) (and (x2,y2)):
+            sign = sign of distance to nearest edge
+            distance = sign * (x1-x^2 + y1-y)^2)^.5
+            ddist_dx = sign * 2*x-2*x1 / (2 * distance^.5)
+            ddist_dy = sign * 2*y-2*y1 / (2 * distance^.5)
+        """
+        if np.all(np.array([x, y]) == self._cache_input):
+            return self._cache_output
+        
+        X, Y = [np.tile(xy, (len(self.x1), 1)).T for xy in [x, y]]  # dim = (ntb, nEdges)
+        X1, Y1, X2, Y2, ddist_dX, ddist_dY = [np.tile(xy, (len(x), 1))
+                                              for xy in [self.x1, self.y1, self.x2, self.y2, self.dEdgeDist_dx, self.dEdgeDist_dy]]
+
+        # perpendicular distance to edge (dot product)
+        d12 = (self.x1 - X) * self.edge_unit_vec[0] + (self.y1 - Y) * self.edge_unit_vec[1]
+
+        # nearest point on edge
+        px = X + d12 * self.edge_unit_vec[0]
+        py = Y + d12 * self.edge_unit_vec[1]
+
+        # distance to start and end points
+        d1 = np.sqrt((self.x1 - X)**2 + (self.y1 - Y)**2)
+        d2 = np.sqrt((self.x2 - X)**2 + (self.y2 - Y)**2)
+
+        # use start or end point if nearest point is outside edge
+        use_xy1 = (((self.dx != 0) & (px < self.x1) & (self.x1 < self.x2)) |
+                   ((self.dx != 0) & (px > self.x1) & (self.x1 > self.x2)) |
+                   ((self.dx == 0) & (py < self.y1) & (self.y1 < self.y2)) |
+                   ((self.dx == 0) & (py > self.y1) & (self.y1 > self.y2)))
+        use_xy2 = (((self.dx != 0) & (px > self.x2) & (self.x2 > self.x1)) |
+                   ((self.dx != 0) & (px < self.x2) & (self.x2 < self.x1)) |
+                   ((self.dx == 0) & (py > self.y2) & (self.y2 > self.y1)) |
+                   ((self.dx == 0) & (py < self.y2) & (self.y2 < self.y1)))
+
+        px[use_xy1] = X1[use_xy1]
+        py[use_xy1] = Y1[use_xy1]
+        px[use_xy2] = X2[use_xy2]
+        py[use_xy2] = Y2[use_xy2]
+
+        distance = d12.copy()
+        v = (px[use_xy1] - X[use_xy1]) * self.xy1_vec[0, np.where(use_xy1)[1]] + (py[use_xy1] - Y[use_xy1]) * self.xy1_vec[1, np.where(use_xy1)[1]]
+        sign_use_xy1 = np.choose(v >= 0, [-1, 1])
+        v = (px[use_xy2] - X[use_xy2]) * self.xy2_vec[0, np.where(use_xy2)[1]] + (py[use_xy2] - Y[use_xy2]) * self.xy2_vec[1, np.where(use_xy2)[1]]
+        sign_use_xy2 = np.choose(v >= 0, [-1, 1])
+
+        d12[use_xy2]
+        d12[:, 1:][use_xy2[:, :-1]]
+
+        distance[use_xy1] = sign_use_xy1 * d1[use_xy1]
+        distance[use_xy2] = sign_use_xy2 * d2[use_xy2]
+
+        length = np.sqrt((X1[use_xy1] - X[use_xy1])**2 + (Y1[use_xy1] - Y[use_xy1])**2)
+        ddist_dX[use_xy1] = sign_use_xy1 * (2 * X[use_xy1] - 2 * X1[use_xy1]) / (2 * length)
+        ddist_dY[use_xy1] = sign_use_xy1 * (2 * Y[use_xy1] - 2 * Y1[use_xy1]) / (2 * length)
+
+        length = np.sqrt((X2[use_xy2] - X[use_xy2])**2 + (Y2[use_xy2] - Y[use_xy2])**2)
+        ddist_dX[use_xy2] = sign_use_xy2 * (2 * X[use_xy2] - 2 * X2[use_xy2]) / (2 * length)
+        ddist_dY[use_xy2] = sign_use_xy2 * (2 * Y[use_xy2] - 2 * Y2[use_xy2]) / (2 * length)
+
+        closest_edge_index = np.argmin(np.abs(distance), 1)
+
+        self._cache_input = np.array([x,y])
+        self._cache_output = [np.choose(closest_edge_index, v.T) for v in [distance, ddist_dX, ddist_dY]]
+        return self._cache_output
+
+    def compute(self, inputs, outputs):
+        turbineX = inputs['turbineX']
+        turbineY = inputs['turbineY']
+
+        outputs['boundaryDistances'] = self.calc_distance_and_gradients(turbineX, turbineY)[0]
+
+    def compute_partials(self, inputs, partials):
+        turbineX = inputs['turbineX']
+        turbineY = inputs['turbineY']
+
+        _, dx, dy = self.calc_distance_and_gradients(turbineX, turbineY)
+        partials['boundaryDistances', 'turbineX'] = np.diagflat(dx)
+        partials['boundaryDistances', 'turbineY'] = np.diagflat(dy)

topfarm/cost_models/dummy.py

@@ -45,9 +45,10 @@ class DummyCost(ExplicitComponent):
 
 
 class DummyCostPlotComp(PlotComp):
-    def __init__(self, optimal, memory=10):
-        super().__init__(memory)
+    def __init__(self, optimal, memory=10, delay=0.001):
+        super().__init__(memory, delay)
         self.optimal = optimal
+        
 
     def init_plot(self, boundary):
         PlotComp.init_plot(self, boundary)

Colonel @ 63d550ac

-Subproject commit b3d528b5bddb5d62f4a6d9f476d7989f5d0f2ad7
+Subproject commit 63d550ac397e2058fc278ecf1488323440a786fc

topfarm/plotting.py

@@ -18,17 +18,15 @@ def mypause(interval):
                 canvas.draw()
             canvas.start_event_loop(interval)
             return
-        
-        
+
+
 class PlotComp(ExplicitComponent):
-    """
-    Evaluates the equation f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3.
-    """
     colors = ['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 100
 
-    def __init__(self, memory=10):
+    def __init__(self, memory=10, delay=0.001):
         ExplicitComponent.__init__(self)
         self.memory = memory
+        self.delay = delay 
         self.history = []
         self.counter = 0
 
@@ -55,16 +53,11 @@ class PlotComp(ExplicitComponent):
         plt.ylim(ylim)
 
     def compute(self, inputs, outputs):
-        """
-        f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3
-
-        Optimal solution (minimum): x = 6.6667; y = -7.3333
-        """
         x = inputs['turbineX']
         y = inputs['turbineY']
         cost = inputs['cost']
         self.history = [(x.copy(), y.copy())] + self.history[:self.memory]
-        
+
         boundary = inputs['boundary']
         self.init_plot(boundary)
         plt.title(cost)
@@ -74,9 +67,26 @@ class PlotComp(ExplicitComponent):
             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)
-        
-        if self.counter ==0:
+
+        if self.counter == 0:
             plt.pause(.01)
-        mypause(0.01)
-        
+        mypause(self.delay)
+
         self.counter += 1
+        
+class NoPlot(PlotComp):
+    def __init__(self, *args, **kwargs):
+        ExplicitComponent.__init__(self)
+        
+    def show(self):
+        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
\ No newline at end of file

topfarm/topfarm.py

@@ -3,7 +3,8 @@ import time
 from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp
 
 import numpy as np
-from topfarm.constraint_components.boundary_component import BoundaryComp
+from topfarm.constraint_components.boundary_component import BoundaryComp,\
+    PolygonBoundaryComp
 from topfarm.constraint_components.spacing_component import SpacingComp
 import warnings
 
@@ -20,9 +21,12 @@ class TopFarm(object):
                  driver_options={'optimizer': 'SLSQP'}):
 
         self.initial_positions = turbines = np.array(turbines)
-        
+
         n_wt = turbines.shape[0]
-        self.boundardy_comp = BoundaryComp(boundary, n_wt, boundary_type)
+        if boundary_type == 'polygon':
+            self.boundardy_comp = PolygonBoundaryComp(boundary, n_wt)
+        else:
+            self.boundardy_comp = BoundaryComp(boundary, n_wt, boundary_type)
         self.problem = prob = Problem()
         indeps = prob.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
         indeps.add_output('turbineX', turbines[:, 0], units='m')
@@ -45,6 +49,7 @@ class TopFarm(object):
             plot_comp.n_wt = n_wt
             plot_comp.n_vertices = self.boundardy_comp.vertices.shape[0]
             prob.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
+
         self.plot_comp = plot_comp
         prob.model.add_constraint('wtSeparationSquared', lower=np.zeros(int(((n_wt - 1.) * n_wt / 2.))) + (min_spacing)**2)
         prob.model.add_constraint('boundaryDistances', lower=np.zeros(self.boundardy_comp.nVertices * n_wt))
@@ -71,13 +76,13 @@ class TopFarm(object):
         with warnings.catch_warnings():  # suppress OpenMDAO/SLSQP warnings
             warnings.filterwarnings('ignore', "Inefficient choice of derivative mode.  You chose 'rev' for a problem with")
             self.problem.run_model()
-        print ("Evaluated in\t%.3fs" % (time.time() - t))
+        print("Evaluated in\t%.3fs" % (time.time() - t))
         return self.get_cost(), self.turbine_positions
 
     def optimize(self):
         t = time.time()
         self.problem.run_driver()
-        print ("Optimized in\t%.3fs" % (time.time() - t))
+        print("Optimized in\t%.3fs" % (time.time() - t))
         return np.array([self.problem['turbineX'], self.problem['turbineY']]).T
 
     def get_cost(self):
@@ -86,7 +91,7 @@ class TopFarm(object):
     @property
     def boundary(self):
         b = self.boundardy_comp.vertices
-        return np.r_[b,b[:1]]
+        return np.r_[b, b[:1]]
 
     @property
     def turbine_positions(self):
@@ -96,7 +101,7 @@ class TopFarm(object):
 if __name__ == '__main__':
     from topfarm.cost_models.dummy import DummyCostPlotComp, DummyCost
     from topfarm.plotting import PlotComp
-    
+
     n_wt = 4
     random_offset = 5
     optimal = [(3, -3), (7, -7), (4, -3), (3, -7), (-3, -3), (-7, -7), (-4, -3), (-3, -7)][:n_wt]