diff --git a/requirements.txt b/requirements.txt
index 35c4519a8d810178a7ab4194268283be6a056b2a..24a4f4bc3f379b3e16a993590cf807e8e86d9d0b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -24,7 +24,7 @@ rsa~=4.9
pyasn1~=0.4.8
toolz~=0.12.0
utm~=0.7.0
-numpy~=1.22.4
+numpy~=1.23.5
pyzmq~=23.2.1
pexpect~=4.8.0
pymongo~=4.2.0
@@ -49,7 +49,7 @@ HeapDict~=1.0.1
anyio~=3.6.1
sniffio~=1.2.0
Babel~=2.10.3
-setuptools>=65.3.0
+setuptools~=59.6.0
isort~=5.10.1
tomli~=2.0.1
tensorflow~=2.9.1
@@ -64,7 +64,7 @@ ipykernel
nbformat~=5.4.0
cloudpickle~=2.1.0
Cython~=0.29.32
-matplotlib~=3.5.3
+matplotlib~=3.6.2
mpmath
sympy
decorator~=5.1.1
@@ -107,9 +107,9 @@ mistune~=2.0.4
Markdown~=3.4.1
ptyprocess~=0.7.0
autograd~=1.4
-scikit-learn
-tqdm
-Shapely
+scikit-learn~=1.1.2
+tqdm~=4.64.1
+Shapely~=2.0.0
threadpoolctl~=3.1.0
openmdao~=3.16.0
nbclient~=0.6.7
@@ -152,3 +152,4 @@ json5~=0.9.10
click~=8.1.3
windIO~=0.1
fusedwake~=0.1.0
+geopandas~=0.12.2
\ No newline at end of file
diff --git a/topfarm/constraint_components/boundary.py b/topfarm/constraint_components/boundary.py
index 84884910a6ebfe07349a3b1f02a638c509718f44..2b09607ab60cf2e21ea1a526584e7b80db56ce99 100644
--- a/topfarm/constraint_components/boundary.py
+++ b/topfarm/constraint_components/boundary.py
@@ -18,8 +18,8 @@ import matplotlib.pyplot as plt
from matplotlib.pyplot import Circle
# from shapely.ops import cascaded_union
-from shapely.geometry import MultiPoint, Polygon, MultiPolygon
-from shapely.ops import unary_union
+from shapely.geometry import MultiPoint, Polygon, MultiPolygon, LineString
+# from shapely.ops import unary_union
import geopandas as gpd
import topfarm
@@ -33,7 +33,7 @@ class XYBoundaryConstraint(Constraint):
TODO: missing-class-docstring
"""
- def __init__(self, boundary, boundary_type: str = 'convex_hull', units=None, relaxation: bool = False):
+ def __init__(self, boundary, boundary_type: str = 'convex_hull', units=None, relaxation: bool = False, **kwargs):
"""Initialize XYBoundaryConstraint
Parameters
@@ -49,11 +49,13 @@ class XYBoundaryConstraint(Constraint):
- 'polygon': Polygon boundary (may be non-convex). Less suitable for gradient-based optimization \n
- 'rectangle': Smallest axis-aligned rectangle covering the boundary points \n
- 'square': Smallest axis-aligned square covering the boundary points
- - 'multi_polygon': Mulitple polygon boundaries incl. exclusion zones (may be non convex). \n
+ - 'multi_polygon': Multiple polygon boundaries incl. exclusion zones (may be non convex). \n
- 'type_specific': Set of multiple polygon boundaries that depend on the wind turbine type. \n
units :
relaxation :
+ kwargs :
+
"""
self.alpha = None
@@ -109,8 +111,8 @@ class XYBoundaryConstraint(Constraint):
"""
if self.boundary_type in ['multi_polygon', 'turbine_specific']:
- bound_min = np.vstack([(bound).min(0) for bound, _ in self.boundary_comp.boundaries]).min(0)
- bound_max = np.vstack([(bound).max(0) for bound, _ in self.boundary_comp.boundaries]).max(0)
+ bound_min = np.vstack([bound.min(0) for bound, _ in self.boundary_comp.boundaries]).min(0)
+ bound_max = np.vstack([bound.max(0) for bound, _ in self.boundary_comp.boundaries]).max(0)
else:
bound_min = self.boundary_comp.xy_boundary.min(0)
bound_max = self.boundary_comp.xy_boundary.max(0)
@@ -222,6 +224,7 @@ class BoundaryBaseComp(ConstraintComponent, ABC):
self.const_id = const_id
self.units = units
self.relaxation = relaxation
+
if np.any(self.xy_boundary[0] != self.xy_boundary[-1]):
self.xy_boundary = np.r_[self.xy_boundary, self.xy_boundary[:1]]
@@ -260,9 +263,9 @@ class BoundaryBaseComp(ConstraintComponent, ABC):
# calculate distances from each point to each face
args = {x: inputs[x] for x in [topfarm.x_key, topfarm.y_key, topfarm.type_key] if x in inputs}
- boundaryDistances = self.distances(**args)
- outputs['boundaryDistances'] = boundaryDistances
- outputs['penalty_' + self.const_id] = np.sum(np.minimum(boundaryDistances, 0) ** 2)
+ boundary_distances = self.distances(**args)
+ outputs['boundaryDistances'] = boundary_distances
+ outputs['penalty_' + self.const_id] = np.sum(np.minimum(boundary_distances, 0) ** 2)
def compute_partials(self, inputs, partials):
""" return Jacobian dict """
@@ -280,24 +283,31 @@ class BoundaryBaseComp(ConstraintComponent, ABC):
def plot(self, ax):
"""Plot boundary"""
+
if isinstance(self, TurbineSpecificBoundaryComp):
linestyles = ['--', '-']
colors = np.array(['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 10)
for n, t in enumerate(self.types):
- line, = ax.plot(*self.ts_merged_xy_boundaries[n][0][0][0, :], color=colors[t], linewidth=1, label=f'{self.wind_turbines._names[n]} boundary')
+ ax.plot(*self.ts_merged_xy_boundaries[n][0][0][0, :], color=colors[t], linewidth=1,
+ label=f'{self.wind_turbines._names[n]} boundary')
+
for bound, io in self.ts_merged_xy_boundaries[n]:
ax.plot(np.asarray(bound)[:, 0].tolist() + [np.asarray(bound)[0, 0]],
- np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], color=colors[t], linewidth=1, linestyle=linestyles[io])
+ np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], color=colors[t], linewidth=1,
+ linestyle=linestyles[io])
+
elif isinstance(self, MultiPolygonBoundaryComp):
colors = ['--k', 'k']
if self.relaxation != 0:
for bound, io in self.relaxed_polygons():
ax.plot(np.asarray(bound)[:, 0].tolist() + [np.asarray(bound)[0, 0]],
- np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], c='r', linewidth=1, linestyle='--')
+ np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], c='r', linewidth=1,
+ linestyle='--')
ax.plot([], c='r', linewidth=1, linestyle='--', label='Relaxed boundaries')
for bound, io in self.boundaries:
ax.plot(np.asarray(bound)[:, 0].tolist() + [np.asarray(bound)[0, 0]],
np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], colors[io], linewidth=1)
+
else:
ax.plot(self.xy_boundary[:, 0].tolist() + [self.xy_boundary[0, 0]],
self.xy_boundary[:, 1].tolist() + [self.xy_boundary[0, 1]], 'k', linewidth=1)
@@ -399,8 +409,8 @@ class ConvexBoundaryComp(BoundaryBaseComp):
unit_normals = self.unit_normals
# initialize array to hold distances from each point to each face
- dfaceDistance_dx = np.zeros([self.n_wt * self.nVertices, self.n_wt])
- dfaceDistance_dy = np.zeros([self.n_wt * self.nVertices, self.n_wt])
+ dface_distance_dx = np.zeros([self.n_wt * self.nVertices, self.n_wt])
+ dface_distance_dy = np.zeros([self.n_wt * self.nVertices, self.n_wt])
for i in range(0, self.n_wt):
# Determine whether the point is inside or outside each face, and the distances from each face
@@ -411,43 +421,40 @@ class ConvexBoundaryComp(BoundaryBaseComp):
# Determination of the derivatives of the perpendicular distances from a point to the 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]
+ diff_dist_vec_dx = np.vdot(dpa_dx, unit_normals[j]) * unit_normals[j]
+ diff_dist_vec_dy = np.vdot(dpa_dy, unit_normals[j]) * unit_normals[j]
# Calculate derivatives for the sign of perpendicular distances from the point to the 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])
+ dface_distance_dx[i * self.nVertices + j, i] = np.vdot(diff_dist_vec_dx, unit_normals[j])
+ dface_distance_dy[i * self.nVertices + j, i] = np.vdot(diff_dist_vec_dy, unit_normals[j])
# return Jacobian dict
- self.dfaceDistance_dx = dfaceDistance_dx
- self.dfaceDistance_dy = dfaceDistance_dy
+ self.dfaceDistance_dx = dface_distance_dx
+ self.dfaceDistance_dy = dface_distance_dy
- def calculate_distance_to_boundary(self, points):
+ def calculate_distance_to_boundary(self, points) -> NDArray:
"""
:param points: points that you want to calculate the distances from to the faces of the convex hull
- :return face_distace: signed perpendicular distances from each point to each face; + is inside
+ :return face_distance: signed perpendicular distances from each point to each face; + is inside
"""
- nPoints = np.array(points).shape[0]
xy_boundary = self.xy_boundary[:-1]
- nVertices = xy_boundary.shape[0]
vertices = xy_boundary
unit_normals = self.unit_normals
- # initialize array to hold distances from each point to each face
- face_distance = np.zeros([nPoints, nVertices])
# define the vector from the point of interest to the first point of the face
- PA = (vertices[:, na] - points[na])
+ pa = (vertices[:, na] - points[na])
# find perpendicular distances from point to current surface (vector projection)
- dist = np.sum(PA * unit_normals[:, na], 2)
+ dist = np.sum(pa * unit_normals[:, na], 2)
+
# calculate the sign of perpendicular distances from point to current face (+ is inside, - is outside)
d_vec = dist[:, :, na] * unit_normals[:, na]
face_distance = np.sum(d_vec * unit_normals[:, na], 2)
return face_distance.T
- def distances(self, x: NDArray, y: NDArray):
+ def distances(self, x: NDArray, y: NDArray) -> NDArray:
"""
TODO: missing-method-docstring
"""
@@ -467,25 +474,25 @@ class ConvexBoundaryComp(BoundaryBaseComp):
"""
x, y = [np.asarray(state[xyz], dtype=float) for xyz in [topfarm.x_key, topfarm.y_key]]
- dist = self.distances(x, y)
- dist = np.where(dist < 0, np.minimum(dist, -.01), dist)
+ np_dist = self.distances(x, y)
+ np_dist = np.where(np_dist < 0, np.minimum(np_dist, -.01), np_dist)
dx, dy = self.gradients(x, y) # independent of position
dx = dx[:self.nVertices, 0]
dy = dy[:self.nVertices, 0]
- for i in np.where(dist.min(1) < 0)[0]: # loop over turbines that violate edges
+ for i in np.where(np_dist.min(1) < 0)[0]: # loop over turbines that violate edges
# find the smallest movement that where the constraints are satisfied
- d = dist[i]
+ d = np_dist[i]
v = np.linspace(-np.abs(d.min()), np.abs(d.min()), 100)
- X, Y = np.meshgrid(v, v)
- m = np.ones_like(X)
+ np_grid_x, np_grid_y = np.meshgrid(v, v)
+ m = np.ones_like(np_grid_x)
- for dx_, dy_, d in zip(dx, dy, dist.T):
- m = np.logical_and(m, X * dx_ + Y * dy_ >= -d[i])
+ for dx_, dy_, d in zip(dx, dy, np_dist.T):
+ m = np.logical_and(m, np_grid_x * dx_ + np_grid_y * dy_ >= -d[i])
- index = np.argmin(X[m] ** 2 + Y[m] ** 2)
- x[i] += X[m][index]
- y[i] += Y[m][index]
+ index = np.argmin(np_grid_x[m] ** 2 + np_grid_y[m] ** 2)
+ x[i] += np_grid_x[m][index]
+ y[i] += np_grid_y[m][index]
state[topfarm.x_key] = x
state[topfarm.y_key] = y
@@ -517,7 +524,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
@staticmethod
def get_boundary_properties(xy_boundary, inclusion_zone=True):
"""
- TODO: missing-method-docstring
+ TODO: missing-staticmethod-docstring
"""
vertices = np.array(xy_boundary)
@@ -561,6 +568,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
boundary_properties = boundary_properties or self.boundary_properties[1:]
A, B, AB, AB_len, edge_unit_normal, A_normal, B_normal = boundary_properties
+
"""
A: edge start point
B: edge end point
@@ -599,12 +607,12 @@ class PolygonBoundaryComp(BoundaryBaseComp):
# Perpendicular distances to edge (AP dot edge_unit_normal product).
# This is the distance to the edge if not use_A or use_B
- distance = np.sum((AP) * edge_unit_normal, 0)
+ distance = np.sum(AP * edge_unit_normal, 0)
# Update distance for points closer to A
- good_side_of_A = (np.sum((AP * A_normal)[:, use_A], 0) > 0)
- sign_use_A = np.where(good_side_of_A, 1, -1)
- distance[use_A] = (vec_len(AP[:, use_A]) * sign_use_A)
+ good_side_of_a = (np.sum((AP * A_normal)[:, use_A], 0) > 0)
+ sign_use_a = np.where(good_side_of_a, 1, -1)
+ distance[use_A] = (vec_len(AP[:, use_A]) * sign_use_a)
# Update distance for points closer to B
good_side_of_B = np.sum((BP * B_normal)[:, use_B], 0) > 0
@@ -617,14 +625,14 @@ class PolygonBoundaryComp(BoundaryBaseComp):
# Gradient of perpendicular distances to edge.
# This is the gradient if not use_A or use_B
- ddist_dxy = np.tile(edge_unit_normal, (1, len(x), 1))
+ diff_dist_dxy = np.tile(edge_unit_normal, (1, len(x), 1))
# Update gradient for points closer to A or B
- ddist_dxy[:, use_A] = sign_use_A * (AP[:, use_A] / vec_len(AP[:, use_A]))
- ddist_dxy[:, use_B] = sign_use_B * (BP[:, use_B] / vec_len(BP[:, use_B]))
- ddist_dX, ddist_dY = ddist_dxy
+ diff_dist_dxy[:, use_A] = sign_use_a * (AP[:, use_A] / vec_len(AP[:, use_A]))
+ diff_dist_dxy[:, use_B] = sign_use_B * (BP[:, use_B] / vec_len(BP[:, use_B]))
+ diff_dist_dx, diff_dist_dy = diff_dist_dxy
- return distance, ddist_dX, ddist_dY
+ return distance, diff_dist_dx, diff_dist_dy
def calc_distance_and_gradients(self, x: NDArray, y: NDArray):
"""
@@ -633,10 +641,10 @@ class PolygonBoundaryComp(BoundaryBaseComp):
if np.all(np.array([x, y]) == self._cache_input):
return self._cache_output
- distance, ddist_dX, ddist_dY = self._calc_distance_and_gradients(x, y)
+ distance, diff_dist_dx, diff_dist_dy = self._calc_distance_and_gradients(x, y)
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]]
+ self._cache_output = [np.choose(closest_edge_index, v.T) for v in [distance, diff_dist_dx, diff_dist_dy]]
return self._cache_output
@@ -706,14 +714,14 @@ class CircleBoundaryComp(PolygonBoundaryComp):
circle = Circle(self.center, self.radius, color='k', fill=False)
ax.add_artist(circle)
- def distances(self, x, y):
+ def distances(self, x: NDArray, y: NDArray) -> NDArray:
"""
TODO: missing-method-docstring
"""
return self.radius - np.sqrt((x - self.center[0]) ** 2 + (y - self.center[1]) ** 2)
- def gradients(self, x, y):
+ def gradients(self, x: NDArray, y: NDArray) -> (NDArray, NDArray):
"""
TODO: missing-method-docstring
"""
@@ -729,7 +737,23 @@ class CircleBoundaryComp(PolygonBoundaryComp):
class Zone(object):
+ """
+ TODO: missing-class-docstring
+ """
def __init__(self, boundary, dist2wt, geometry_type, incl, name):
+ """
+ TODO: missing-init-docstring
+
+ Parameters
+ ----------
+ boundary:
+ dist2wt:
+ geometry_type:
+ incl:
+ name:
+
+ """
+
self.name = name
self.boundary = boundary
self.dist2wt = dist2wt
@@ -738,11 +762,17 @@ class Zone(object):
class InclusionZone(Zone):
+ """
+ TODO: missing-class-docstring
+ """
def __init__(self, boundary, dist2wt=None, geometry_type='polygon', name=''):
super().__init__(np.asarray(boundary), dist2wt, geometry_type, incl=1, name=name)
class ExclusionZone(Zone):
+ """
+ TODO: missing-class-docstring
+ """
def __init__(self, boundary, dist2wt=None, geometry_type='polygon', name=''):
super().__init__(np.asarray(boundary), dist2wt, geometry_type, incl=0, name=name)
@@ -758,71 +788,59 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
Parameters
----------
- n_wt : TYPE
- DESCRIPTION.
- zones : list
+ n_wt :
+ turbine identification
+ zones :
list of InclusionZone and ExclusionZone objects
const_id : TYPE, optional
- DESCRIPTION. The default is None.
+ The default is None.
units : TYPE, optional
- DESCRIPTION. The default is None.
+ The default is None.
method : {'nearest' or 'smooth_min'}, optional
- 'nearest' calculate the distance to the nearest edge or point'smooth_min'
+ 'nearest' calculate the distance to the nearest edge or point 'smooth_min'
calculates the weighted minimum distance to all edges/points. The default is 'nearest'.
simplify_geometry : bool or dict
- if float, simplification tolerance. if dict, shapely simplify keyword arguments
-
- Returns
- -------
- None.
+ if is a float, simplification tolerance. if dict, shapely simplify keyword arguments
"""
self.zones = zones
- self.bounds_poly, xy_boundaries = self.get_xy_boundaries()
- PolygonBoundaryComp.__init__(self, n_wt, xy_boundary=xy_boundaries[0], const_id=const_id, units=units, relaxation=relaxation)
- # self.bounds_poly = [Polygon(x) for x in xy_boundaries]
- self.incl_excls = [x.incl for x in zones]
- self._setup_boundaries()
self.relaxation = relaxation
self.method = method
- self.xy_multi_boundary = xy_multi_boundary
+ self.bounds_poly, xy_boundaries = self.get_xy_boundaries()
+
+ # inclusion zones only
+ PolygonBoundaryComp.__init__(self, n_wt,
+ xy_boundary=xy_boundaries[0],
+ const_id=const_id,
+ units=units,
+ relaxation=relaxation)
+ self.incl_excls = [x.incl for x in zones]
+ self._setup_boundaries()
+
+ ###################################################################
+ # Geopandas pre processing
list_mpolygon = []
list_zn_type = []
- for poly_i, zn_i in xy_multi_boundary:
+ for poly_i, zn_i in zones:
list_mpolygon.append(Polygon(poly_i))
if zn_i == 'i' or zn_i is True or zn_i == 1:
list_zn_type.append(True)
else:
list_zn_type.append(False)
- geo_mpolygon = gpd.GeoDataFrame(dict(inclusion=list_zn_type, geometry=list_mpolygon))
- geo_mpolygon_inz = gpd.GeoSeries(geo_mpolygon[geo_mpolygon.inclusion == True].unary_union)
-
- shp_cvx_hull = geo_mpolygon_inz.convex_hull
- np_multi_boundary = np.asarray(shp_cvx_hull[0].boundary.xy).T
-
- # inclusion zones only
- PolygonBoundaryComp.__init__(self, n_wt,
- xy_boundary=np_multi_boundary,
- const_id=const_id,
- units=units,
- relaxation=relaxation)
-
- self.bounds_poly = list_mpolygon
+ self.geo_mpolygon = gpd.GeoDataFrame(dict(inclusion=list_zn_type, geometry=list_mpolygon))
self.types_bool = np.array(list_zn_type, dtype=int).tolist()
- self._setup_boundaries()
- self.relaxation = relaxation
- self.method = method
+ ###################################################################
# FIXME : shapely/shapely/geometry/base.py -> line 452, def simplify(self, tolerance, preserve_topology=True)
if simplify_geometry:
self.simplify(simplify_geometry)
- def simplify(self, simplify_geometry):
+ def simplify(self, simplify_geometry: bool = False):
"""
TODO: missing-method-docstring
"""
@@ -831,6 +849,7 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
self.bounds_poly = [rp.simplify(**simplify_geometry) for rp in self.bounds_poly]
else:
self.bounds_poly = [rp.simplify(simplify_geometry) for rp in self.bounds_poly]
+
self._setup_boundaries()
def get_xy_boundaries(self):
@@ -846,10 +865,12 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
buffer = z.dist2wt(**{k: 100 for k in z.dist2wt.__code__.co_varnames})
else:
buffer = 0
+
if z.geometry_type == 'line':
poly = Polygon(LineString(z.boundary).buffer(buffer, join_style=2).exterior)
elif z.geometry_type == 'polygon':
poly = Polygon(z.boundary).buffer(buffer, join_style=2)
+
polygons.append(poly)
bounds.append(np.asarray(poly.exterior.coords))
@@ -871,13 +892,16 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
@staticmethod
def _poly_to_bound(polygons):
+
boundaries = []
for bound in polygons:
x, y = bound.exterior.xy
boundaries.append((np.asarray([x, y]).T[:-1, :], 1))
+
for interior in bound.interiors:
x, y = interior.xy
boundaries.append((np.asarray([x, y]).T[:-1, :], 0))
+
return boundaries
def _calc_resulting_polygons(self, boundary_polygons: list):
@@ -885,7 +909,7 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
Parameters
----------
boundary_polygons :
- list of shapely polygons as specifed or inferred from user input
+ list of shapely polygons as specified or inferred from user input
Returns
-------
list of merged shapely polygons. Resolve issues arising if any are overlapping, touching, or contained in
@@ -950,7 +974,7 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
return domain
@staticmethod
- def sign(dist_ij):
+ def sign(dist_ij: NDArray) -> NDArray:
"""
TODO: docstring
"""
@@ -981,26 +1005,28 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
if np.all(np.array([x, y]) == self._cache_input) & (not self.relaxation):
return self._cache_output
- dist_ij, ddist_d_x, ddist_d_y = self._calc_distance_and_gradients(x, y, self.boundary_properties_list_all)
+ dist_ij, diff_dist_d_x, diff_dist_d_y = self._calc_distance_and_gradients(x, y,
+ self.boundary_properties_list_all)
- dDdk_ijk = np.moveaxis([ddist_d_x, ddist_d_y], 0, -1)
+ diff_dk_ijk = np.moveaxis([diff_dist_d_x, diff_dist_d_y], 0, -1)
sign_i = self.sign(dist_ij)
self._cache_input = np.array([x, y])
- self._cache_output = [dist_ij, dDdk_ijk, sign_i]
+ self._cache_output = [dist_ij, diff_dk_ijk, sign_i]
return self._cache_output
- def calc_relaxation(self, iteration_no=None):
+ def calc_relaxation(self, iteration_no=None) -> float:
"""
- The tupple relaxation contains a first term for the penalty constant
+ The tuple relaxation contains a first term for the penalty constant
and a second term for the n first iterations to apply relaxation.
"""
if iteration_no is None:
iteration_no = self.problem.cost_comp.n_grad_eval + 1
+
return max(0, self.relaxation[0] * (self.relaxation[1] - iteration_no))
- def distances(self, x, y):
+ def distances(self, x: NDArray, y: NDArray) -> NDArray:
"""
TODO: docstring
@@ -1030,9 +1056,9 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
return dist_i
- def gradients(self, x, y):
+ def gradients(self, x: NDArray, y: NDArray) -> NDArray:
"""
- The derivate of the smooth maximum with respect to x and y is calculated with the chain rule:
+ The derivative of the smooth maximum with respect to x and y is calculated with the chain rule:
dS/dk = dS/dD * dD/dk
where S is smooth maximum, D is distance to edge and k is the spacial dimension
@@ -1047,7 +1073,7 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
dist_ij, d_ddk_ijk, _ = self.calc_distance_and_gradients(x, y)
if self.relaxation:
- Dist_ij += self.calc_relaxation()
+ dist_ij += self.calc_relaxation()
if self.method == 'smooth_min':
d_sd_dist_ij = smooth_max_gradient(np.abs(dist_ij), -np.abs(dist_ij).max(), axis=1)
@@ -1063,13 +1089,175 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
return gradients
+class TurbineSpecificBoundaryComp(MultiPolygonBoundaryComp):
+ """
+ TODO: docstring
+ """
+
+ def __init__(self, n_wt: int, wind_turbines, zones: list, const_id: int = None,
+ units=None, relaxation: bool = False, method: str = 'nearest', simplify_geometry: bool = False):
+ """
+ TODO: docstring
+ """
+
+ # self.dependencies = [d or {'type': None, 'multiplier': None, 'ref': None} for d in dependencies]
+ # self.multi_boundary = xy_boundaries = self.get_xy_boundaries(boundaries, geometry_types, incl_excls)
+ self.wind_turbines = wind_turbines
+ self.types = wind_turbines.types()
+
+ self.n_wt = n_wt
+ self.zones = zones
+ self.ts_polygon_boundaries, ts_xy_boundaries = self.get_ts_boundaries()
+ MultiPolygonBoundaryComp.__init__(self, n_wt=n_wt, zones=zones, const_id=const_id, units=units,
+ relaxation=relaxation, method=method, simplify_geometry=simplify_geometry)
+
+ self.ts_merged_polygon_boundaries = self.merge_boundaries()
+ self.ts_merged_xy_boundaries = self.get_ts_xy_boundaries()
+ self.ts_boundary_properties = self.get_ts_boundary_properties()
+ self.ts_item_indices = self.get_ts_item_indices()
+
+ def get_ts_boundaries(self):
+ """
+ TODO: docstring
+ """
+
+ polygons, bounds = [], []
+
+ for t in set(self.types):
+ temp1 = []
+ temp2 = []
+ dist2wt_input = dict(D=self.wind_turbines.diameter(t),
+ H=self.wind_turbines.hub_height(t))
+ for z in self.zones:
+ if hasattr(z.dist2wt, '__code__'):
+ buffer = z.dist2wt(**{k: dist2wt_input[k] for k in z.dist2wt.__code__.co_varnames})
+ else:
+ buffer = 0
+ if z.geometry_type == 'line':
+ poly = Polygon(LineString(z.boundary).buffer(buffer, join_style=2).exterior)
+ elif z.geometry_type == 'polygon':
+ poly = Polygon(z.boundary).buffer(buffer, join_style=2)
+ bound = np.asarray(poly)
+ temp1.append(poly)
+ temp2.append(bound)
+ polygons.append(temp1)
+ bounds.append(temp2)
+
+ return polygons, bounds
+
+ def get_ts_xy_boundaries(self):
+ """
+ TODO: docstring
+ """
+
+ return [self._poly_to_bound(b) for b in self.ts_merged_polygon_boundaries]
+
+ def merge_boundaries(self):
+ """
+ TODO: docstring
+ """
+
+ return [self._calc_resulting_polygons(bounds) for bounds in self.ts_polygon_boundaries]
+
+ def get_ts_boundary_properties(self, ):
+ """
+ TODO: docstring
+ """
+
+ return [[self.get_boundary_properties(bound) for bound, _ in bounds] for bounds in self.ts_merged_xy_boundaries]
+
+ def get_ts_item_indices(self):
+ """
+ TODO: docstring
+ """
+
+ temp = []
+ for bounds in self.ts_merged_xy_boundaries:
+ n_edges = np.asarray([len(bound) for bound, _ in bounds])
+ n_edges_tot = np.sum(n_edges)
+ start_at = np.cumsum(n_edges) - n_edges
+ end_at = start_at + n_edges
+ item_indices = [n_edges_tot, start_at, end_at]
+ temp.append(item_indices)
+
+ return temp
+
+ def calc_distance_and_gradients(self, x: NDArray, y: NDArray, types=None):
+ """
+ TODO: docstring
+ """
+
+ if not isinstance(self._cache_input, type(None)):
+ if np.all(np.array([x, y]) == self._cache_input[0]) & (not self.relaxation) & np.all(
+ np.asarray([types]) == self._cache_input[1]):
+ return self._cache_output
+ if isinstance(types, type(None)):
+ types = np.zeros(self.n_wt)
+
+ dist_i = np.zeros(self.n_wt)
+ sign_i = np.zeros(self.n_wt)
+ diff_dx_i = np.zeros(self.n_wt)
+ diff_dy_i = np.zeros(self.n_wt)
+
+ for t in set(types):
+ t = int(t)
+ idx = (types == t)
+ n_edges_tot, start_at, end_at = self.ts_item_indices[t]
+ dist_ij = np.zeros((sum(idx), n_edges_tot))
+ diff_dk_ijk = np.zeros((sum(idx), n_edges_tot, 2))
+ for n, (bound, bound_type) in enumerate(self.ts_merged_xy_boundaries[t]):
+ sa = start_at[n]
+ ea = end_at[n]
+ distance, diff_dist_dX, diff_dist_dY = self._calc_distance_and_gradients(x[idx], y[idx],
+ self.ts_boundary_properties[t][n][1:])
+ if bound_type == 0:
+ distance *= -1
+ diff_dist_dX *= -1
+ diff_dist_dY *= -1
+ dist_ij[:, sa:ea] = distance
+ diff_dk_ijk[:, sa:ea, 0] = diff_dist_dX
+ diff_dk_ijk[:, sa:ea, 1] = diff_dist_dY
+
+ sign_i[idx] = self.sign(dist_ij)
+ dist_i[idx] = dist_ij[np.arange(sum(idx)), np.argmin(np.abs(dist_ij), axis=1)]
+ diff_dx_i[idx], diff_dy_i[idx] = diff_dk_ijk[np.arange(sum(idx)), np.argmin(np.abs(dist_ij), axis=1), :].T
+
+ self._cache_input = (np.array([x, y]), np.asarray(types))
+ self._cache_output = [dist_i, diff_dx_i, diff_dy_i, sign_i]
+ return self._cache_output
+
+ def distances(self, x: NDArray, y: NDArray, type=None):
+ """
+ TODO: docstring
+ """
+
+ dist_i, _, _, _ = self.calc_distance_and_gradients(x, y, types=type)
+ if self.relaxation:
+ dist_i += self.calc_relaxation()
+ return dist_i
+
+ def gradients(self, x: NDArray, y: NDArray, type=None) -> NDArray:
+ """
+ TODO: docstring
+ """
+
+ dist_i, diff_dx_i, diff_dy_i, _ = self.calc_distance_and_gradients(x, y, types=type)
+ if self.relaxation:
+ dist_i += self.calc_relaxation()
+ diff_dt = -self.relaxation[0]
+ if self.relaxation:
+ gradients = np.diagflat(diff_dx_i), np.diagflat(diff_dy_i), np.ones(self.n_wt) * diff_dt
+ else:
+ gradients = np.diagflat(diff_dx_i), np.diagflat(diff_dy_i)
+
+ return gradients
class LevelSetBoundary:
"""
TODO: docstring
"""
- def __init__(self, xy_multi_boundary: [Polygon, bool], ndiv: int = 100, bounds=[]):
+ def __init__(self, xy_multi_boundary: [Polygon, bool], ndiv: int = 100, bounds: list = [float]):
"""
TODO: docstring
@@ -1094,7 +1282,7 @@ class LevelSetBoundary:
self.np_xgrid: NDArray = np.empty([], dtype=float)
self.np_ygrid: NDArray = np.empty([], dtype=float)
- # Equivalent Geomtry Polygon (Inclusion/ Exclusion) - GeoPandas Polygon/MultiPolygon
+ # Equivalent Geometry Polygon (Inclusion/ Exclusion) - GeoPandas Polygon/MultiPolygon
self.geo_inz = gpd.GeoSeries()
self.geo_exz = gpd.GeoSeries()
self.geo_eqv = gpd.GeoSeries()
@@ -1210,7 +1398,7 @@ class LevelSetBoundary:
TODO: docstring
"""
- # FIXME, shoud be supplied by BoundaryXY Class (in caso of only exclusions zones)
+ # FIXME, should be supplied by BoundaryXY Class (in caso of only exclusions zones)
if len(self.bounds) == 0:
pd_coords = self.geo_eqv.bounds
else:
@@ -1276,7 +1464,7 @@ class LevelSetBoundary:
return x, y
@staticmethod
- def _ks_aggreation(np_gx: NDArray, stab: float = 0., rho: float = 5.) -> float:
+ def _ks_aggregation(np_gx: NDArray, stab: float = 0., rho: float = 5.) -> float:
"""
The Kreisselmeier-Steinhauser (KS) function for aggregate constraints
@@ -1320,9 +1508,9 @@ class LevelSetBoundary:
np_dist_eq = np.zeros(np_dist.shape[1], dtype=float)
for i in range(np_dist.shape[1]):
- np_dist_eq[i] = -self._ks_aggreation(-np_dist[:, i], rho=self.rho)
+ np_dist_eq[i] = -self._ks_aggregation(-np_dist[:, i], rho=self.rho)
- return -self._ks_aggreation(-np_dist.min(axis=0), rho=self.rho)
+ return -self._ks_aggregation(-np_dist.min(axis=0), rho=self.rho)
def ks_dist_diff(self, x: NDArray, y: NDArray, fdtype: str = 'second-order', **kwargs) -> Tuple[NDArray, NDArray]:
"""
@@ -1523,25 +1711,25 @@ def main():
ExclusionZone(e5, name='e5', dist2wt=lambda: 1, geometry_type='line'),
]
- N_points = 50
- xs = np.linspace(-1, 30, N_points)
- ys = np.linspace(-1, 30, N_points)
+ n_points = 50
+ xs = np.linspace(-1, 30, n_points)
+ ys = np.linspace(-1, 30, n_points)
y_grid, x_grid = np.meshgrid(xs, ys)
x = x_grid.ravel()
y = y_grid.ravel()
n_wt = len(x)
- MPBC = MultiPolygonBoundaryComp(n_wt, zones, method='nearest')
- distances = MPBC.distances(x, y)
+ mpbc = MultiPolygonBoundaryComp(n_wt, zones, method='nearest')
+ distances = mpbc.distances(x, y)
delta = 1e-9
- distances2 = MPBC.distances(x + delta, y)
+ distances2 = mpbc.distances(x + delta, y)
dx_fd = (distances2 - distances) / delta
- dx = np.diag(MPBC.gradients(x + delta / 2, y)[0])
+ dx = np.diag(mpbc.gradients(x + delta / 2, y)[0])
plt.figure()
plt.plot(dx_fd, dx, '.')
plt.figure()
- for n, bound in enumerate(MPBC.boundaries):
+ for n, bound in enumerate(mpbc.boundaries):
x_bound, y_bound = bound[0].T
x_bound = np.append(x_bound, x_bound[0])
y_bound = np.append(y_bound, y_bound[0])
@@ -1551,26 +1739,26 @@ def main():
plt.legend()
plt.grid()
plt.axis('square')
- plt.contourf(x_grid, y_grid, distances.reshape(N_points, N_points), np.linspace(-10, 10, 100), cmap='seismic')
+ plt.contourf(x_grid, y_grid, distances.reshape(n_points, n_points), np.linspace(-10, 10, 100), cmap='seismic')
plt.colorbar()
plt.figure()
ax = plt.axes(projection='3d')
ax.contour3D(
x.reshape(
- N_points, N_points), y.reshape(
- N_points, N_points), distances.reshape(
- N_points, N_points), np.linspace(-10, 10, 100), cmap='seismic')
+ n_points, n_points), y.reshape(
+ n_points, n_points), distances.reshape(
+ n_points, n_points), np.linspace(-10, 10, 100), cmap='seismic')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
if 0:
for smpl in [0, 1, 2, 3, 4, 5, 6, 7, 8]:
- MPBC = MultiPolygonBoundaryComp(n_wt, zones, simplify_geometry=smpl)
+ mpbc = MultiPolygonBoundaryComp(n_wt, zones, simplify_geometry=smpl)
plt.figure()
ax = plt.gca()
- MPBC.plot(ax)
+ mpbc.plot(ax)
wind_turbines = WindTurbines(names=['tb1', 'tb2'],
diameters=[80, 120],
@@ -1610,28 +1798,29 @@ def main():
ExclusionZone(b3, geometry_type='line', dist2wt=lambda D: 3 * D, name='river'),
ExclusionZone(b4, geometry_type='line', dist2wt=lambda D, H: max(D * 2, H * 3), name='road'),
]
- N_points = 50
- xs = np.linspace(0, 3000, N_points)
- ys = np.linspace(0, 3000, N_points)
+ n_points = 50
+ xs = np.linspace(0, 3000, n_points)
+ ys = np.linspace(0, 3000, n_points)
y_grid, x_grid = np.meshgrid(xs, ys)
x = x_grid.ravel()
y = y_grid.ravel()
n_wt = len(x)
types = np.zeros(n_wt)
- TSBC = TurbineSpecificBoundaryComp(n_wt, wind_turbines, zones)
- distances = TSBC.distances(x, y, type=types)
+
+ tsbc = TurbineSpecificBoundaryComp(n_wt, wind_turbines, zones)
+ distances = tsbc.distances(x, y, type=types)
delta = 1e-9
- distances2 = TSBC.distances(x + delta, y, type=types)
+ distances2 = tsbc.distances(x + delta, y, type=types)
dx_fd = (distances2 - distances) / delta
- dx = np.diag(TSBC.gradients(x + delta / 2, y, type=types)[0])
+ dx = np.diag(tsbc.gradients(x + delta / 2, y, type=types)[0])
plt.figure()
plt.plot(dx_fd, dx, '.')
plt.figure()
- for ll, t in enumerate(TSBC.types):
- line, = plt.plot(*TSBC.ts_merged_xy_boundaries[ll][0][0][0, :], label=f'type {ll}')
- for n, bound in enumerate(TSBC.ts_merged_xy_boundaries[ll]):
+ for ll, t in enumerate(tsbc.types):
+ line, = plt.plot(*tsbc.ts_merged_xy_boundaries[ll][0][0][0, :], label=f'type {ll}')
+ for n, bound in enumerate(tsbc.ts_merged_xy_boundaries[ll]):
x_bound, y_bound = bound[0].T
x_bound = np.append(x_bound, x_bound[0])
y_bound = np.append(y_bound, y_bound[0])
@@ -1641,9 +1830,9 @@ def main():
plt.grid()
plt.axis('square')
- for ll, t in enumerate(TSBC.types):
+ for ll, t in enumerate(tsbc.types):
plt.figure()
- for n, bound in enumerate(TSBC.ts_merged_xy_boundaries[ll]):
+ for n, bound in enumerate(tsbc.ts_merged_xy_boundaries[ll]):
x_bound, y_bound = bound[0].T
x_bound = np.append(x_bound, x_bound[0])
y_bound = np.append(y_bound, y_bound[0])
@@ -1651,7 +1840,7 @@ def main():
plt.grid()
plt.title(f'type {ll}')
plt.axis('square')
- plt.contourf(x_grid, y_grid, TSBC.distances(x, y, type=t * np.ones(n_wt)).reshape(N_points, N_points), 50)
+ plt.contourf(x_grid, y_grid, tsbc.distances(x, y, type=t * np.ones(n_wt)).reshape(n_points, n_points), 50)
plt.colorbar()
plt.show()
diff --git a/topfarm/drivers/stochastic_gradient_descent_driver.py b/topfarm/drivers/stochastic_gradient_descent_driver.py
index aa5ebbd8683845fc47f4bba7eba359216cceb56b..fab100660b6f107b5a1b4b649dbb9cb4dc931ae7 100644
--- a/topfarm/drivers/stochastic_gradient_descent_driver.py
+++ b/topfarm/drivers/stochastic_gradient_descent_driver.py
@@ -9,16 +9,30 @@ from openmdao.core.driver import Driver, RecordingDebugging
from six import iteritems
import numpy as np
from openmdao.core.analysis_error import AnalysisError
-import topfarm
import time
-from openmdao.utils.concurrent import concurrent_eval
-from openmdao.utils.record_util import create_local_meta
class SGDDriver(Driver):
- def __init__(self, maxiter, **kwargs):
+ def __init__(self, maxiter: int, **kwargs):
+
+ self.obj_list = []
+ self.con_list = []
+ self.obj_and_con_list = []
+
+ self.upper = None
+ self.lower = None
+
+ self.l0 = None
+ self.mid = None
+
+ self.alpha0 = None
+ self._dvlist = None
+ self.abs2prom = None
self.maxiter = maxiter
+
+ self.is_converged = None
+
super().__init__(**kwargs)
# What we support
@@ -73,7 +87,7 @@ class SGDDriver(Driver):
def run(self):
"""
- Excute the stochastic gradient descent algorithm.
+ Execute the stochastic gradient descent algorithm.
Returns
-------
@@ -102,6 +116,7 @@ class SGDDriver(Driver):
lower_bound[i:j] = meta['lower']
upper_bound[i:j] = meta['upper']
x0[i:j] = desvar_vals[name]
+
maxiter = self.maxiter
max_time = self.options['max_time'] or 1e20
assert maxiter < 1e20 or max_time < 1e20, "maxiter or max_time must be set"
@@ -125,7 +140,7 @@ class SGDDriver(Driver):
desvar_info = [(abs2prom[name], *self._desvar_idx[name], lower_bound, upper_bound)
for name, meta in iteritems(desvars)]
desvar_dict = {name: (x0[i:j].copy(), lbound[i:j], ubound[i:j]) for (name, i, j, lbound, ubound) in desvar_info}
- start = time.time()
+ # start = time.time()
self.obj_list = list(self._objs)
self.con_list = []
@@ -137,6 +152,7 @@ class SGDDriver(Driver):
else:
size = meta['global_size'] if meta['distributed'] else meta['size']
self.con_list.append(name)
+
self.obj_and_con_list = self.obj_list + self.con_list
start = time.time()
@@ -146,6 +162,7 @@ class SGDDriver(Driver):
jac_time = time.time() - start
max_iters_from_max_time = int(max_time / jac_time)
self.maxiter = min(self.maxiter, max_iters_from_max_time)
+
j = jac[0]
learning_rate = float(np.copy(self.learning_rate))
self.alpha0 = np.mean(np.abs(j)) / learning_rate
@@ -166,19 +183,22 @@ class SGDDriver(Driver):
self.upper = mid
elif etaM > self.gamma_min:
self.lower = mid
+
self.mid = mid
m = np.zeros(x0.size)
v = np.zeros(x0.size)
self.is_converged = False
+
while (n_iter < self.maxiter) and (not self.is_converged):
obj_value_x1, x1, alpha, learning_rate, m, v, success = self.objective_callback(x1, alpha, learning_rate, m, v, record=True)
n_iter += 1
if disp:
print(n_iter, obj_value_x1)
+
return False
- def objective_callback(self, x, alpha, learning_rate, m, v, record=False):
+ def objective_callback(self, x: np.ndarray, alpha, learning_rate, m, v, record=False):
"""
Evaluate problem objective at the requested point.
@@ -198,6 +218,7 @@ class SGDDriver(Driver):
int
Case number, used for identification when run in parallel.
"""
+
model = self._problem().model
success = 1
diff --git a/topfarm/easy_drivers.py b/topfarm/easy_drivers.py
index a496a779810e3dada2f175602d751b7deca76f91..f1d998e89a84580c795e727c4099080b025edb41 100644
--- a/topfarm/easy_drivers.py
+++ b/topfarm/easy_drivers.py
@@ -19,17 +19,20 @@ import pandas as pd
import openmdao
import scipy
from scipy import optimize
+from scipy.optimize.optimize import OptimizeResult
# from scipy.spatial import distance
from openmdao.drivers import scipy_optimizer as drv
from openmdao.core.driver import Driver, RecordingDebugging
from openmdao.core.constants import INF_BOUND
-from openmdao.drivers.scipy_optimizer import ScipyOptimizeDriver
+from openmdao.drivers import scipy_optimizer
+# from openmdao.drivers.scipy_optimizer import ScipyOptimizeDriver
# from openmdao.drivers.genetic_algorithm_driver import SimpleGADriver # version 2.5.0 has bug see Issue #874
from topfarm.utils import ks_aggreation
from topfarm.drivers.genetic_algorithm_driver import SimpleGADriver
from topfarm.drivers.random_search_driver import RandomSearchDriver
+from topfarm.drivers.stochastic_gradient_descent_driver import SGDDriver
from topfarm.constraint_components.boundary import LevelSetBoundary
try:
@@ -61,6 +64,13 @@ CITATIONS = """
# from watchpoints import watch
+def fmt_option(v):
+ if isinstance(v, str):
+ return v.encode()
+ else:
+ return v
+
+
# noinspection PyCompatibility
class EasyDriverBase(Driver):
"""
@@ -235,7 +245,7 @@ class EasyDriverBase(Driver):
return True
-class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
+class EasyScipyOptimizeDriver(scipy_optimizer.ScipyOptimizeDriver, EasyDriverBase):
"""
TODO: missing-class-docstring
"""
@@ -245,8 +255,8 @@ class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
"""
Parameters
----------
- optimizer : {'COBYLA', 'SLSQP'}
- Gradients are only supported by SLSQP
+ optimizer : {'COBYLA', 'SLSQP', 'SGD'}
+ Gradients are only supported by SLSQP and SGD
maxiter : int
Maximum number of iterations.
tol : float
@@ -254,33 +264,56 @@ class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
disp : bool
Set False to prevent printing of Scipy convergence messages
auto_scale : bool
- Set to true to set ref0 and ref1 to boundaries of the desig variables for the drivers which support it (SLSQP, ).
+ Set to true to set ref0 and ref1 to boundaries of the design variables for the drivers which
+ support it (SLSQP, ).
"""
- ScipyOptimizeDriver.__init__(self)
- self.auto_scale = auto_scale
+ scipy_optimizer.ScipyOptimizeDriver.__init__(self)
+
+ self.auto_scale = auto_scale
self.cite = CITATIONS
- if optimizer == 'SGD':
- from topfarm.drivers.SGD import SGD
- from openmdao.drivers import scipy_optimizer
- from scipy.optimize.optimize import OptimizeResult
+ if optimizer == 'IPOPT':
+ try:
+ from cyipopt.scipy_interface import minimize_ipopt
+ except ImportError:
+ raise ImportError("""Cannot import ipopt wrapper. Please install cyipopt, e.g. via conda
+ Windows: conda install -c pycalphad cyipopt \n Linux/OSX: conda install -c conda-forge cyipopt """)
+
+ ipopt_options = {k: fmt_option(v) for k, v in kwargs.items()}
+
+ def minimize_ipopt_wrapper(*args, maxiter=200, disp=True, **kwargs):
+ ipopt_options.update({'max_iter': self.max_iter or maxiter, 'print_level': int(disp)})
+ return minimize_ipopt(*args, options=ipopt_options, **kwargs)
+
+ kwargs = {}
+ for lst in [scipy_optimizer._optimizers, scipy_optimizer._gradient_optimizers,
+ scipy_optimizer._bounds_optimizers,
+ scipy_optimizer._all_optimizers, scipy_optimizer._constraint_optimizers,
+ scipy_optimizer._constraint_grad_optimizers]:
+ lst.add(minimize_ipopt_wrapper)
+ optimizer = minimize_ipopt_wrapper
+
+ elif optimizer == 'SGD':
sgd_options = {k: fmt_option(v) for k, v in kwargs.items()}
- sgd_opt = SGD(**kwargs)
+ sgd_opt = SGDDriver(**kwargs)
- def minimize_sgd_wrapper(*args, maxiter=200, disp=True, **kwargs):
+ def minimize_sgd_wrapper(*args, maxiter: int =200, disp: bool = True, **kwargs):
sgd_options.update({'max_iter': self.max_iter or maxiter, 'print_level': int(disp)})
s = sgd_opt.run(*args, options=sgd_options, **kwargs)
return OptimizeResult(x=s, fun=args[0](s), jac=kwargs['jac'](s), nit=int(sgd_opt.T),
nfev=None, njev=None, status=1,
message='hello world!', success=1)
kwargs = {}
- for lst in [scipy_optimizer._optimizers, scipy_optimizer._gradient_optimizers, scipy_optimizer._bounds_optimizers,
- scipy_optimizer._all_optimizers, scipy_optimizer._constraint_optimizers, scipy_optimizer._constraint_grad_optimizers]:
+ for lst in [scipy_optimizer._optimizers, scipy_optimizer._gradient_optimizers,
+ scipy_optimizer._bounds_optimizers, scipy_optimizer._all_optimizers,
+ scipy_optimizer._constraint_optimizers, scipy_optimizer._constraint_grad_optimizers]:
lst.add(minimize_sgd_wrapper)
+
optimizer = minimize_sgd_wrapper
self.options.update({'optimizer': optimizer, 'maxiter': self.max_iter or maxiter, 'tol': tol, 'disp': disp})
+
if kwargs:
self.options.update(kwargs)
@@ -295,7 +328,6 @@ class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
if tuple([int(v) for v in scipy.__version__.split(".")]) < (1, 5, 0):
# Upper and lower disturbs SLSQP when running with constraints. Add limits as constraints
model.add_constraint(desvar_name, kwargs.get('lower', None), kwargs.get('upper', None))
- kwargs = {'lower': np.nan, 'upper': np.nan} # Default +/- sys.float_info.max does not work for SLSQP
ref0 = 0
ref1 = 1
@@ -317,6 +349,10 @@ class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
l, u = desvar_values[1], desvar_values[2]
kwargs = {'ref0': ref0, 'ref': ref1, 'lower': l, 'upper': u}
+ else:
+ # Default +/- sys.float_info.max does not work for SLSQP
+ kwargs = {'lower': np.nan, 'upper': np.nan}
+
return kwargs
@property
@@ -367,6 +403,7 @@ def fd_diff(fun: Any, xi_new: NDArray, idx: int = None, perturbation: float = 1.
np_dx = np.abs(dtp * xi_new)
np_zidx = np.where(np_dx < 1.e-9)[0]
+
if len(np_zidx) > 0:
np_dx[np_zidx] = np_dx.mean()
@@ -513,13 +550,13 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
multistart :
Number of MultiStart Looping
lsearch :
- Number of Iteration throung the local optimization search
+ Number of Iteration through the local optimization search
perturbation :
Perturbation value for the multistart
fd:
Finite difference method for approximating the first-order derivatives of the objective function
- **kwargs : Another paramters that could be included
+ **kwargs : Another parameters that could be included
"""
if Problem is not None:
@@ -2075,22 +2112,21 @@ class EasyRandomSearchDriver(RandomSearchDriver, EasyDriverBase):
class EasySGDDriver(SGDDriver, EasyDriverBase):
- def __init__(self, maxiter=100, max_time=600, disp=False, run_parallel=False,
- learning_rate=10, upper=0.1, lower=0, beta1=0.1, beta2=0.2, gamma_min_factor=1e-2,
- speedupSGD=False, sgd_thresh=0.1):
+ def __init__(self, maxiter: int = 100, max_time: int = 600, disp: bool = False, run_parallel: bool = False,
+ learning_rate:float = 10., upper:float = 0.1, lower:float = 0., beta1:float = 0.1, beta2:float = 0.2,
+ gamma_min_factor:float = 1.e-2, speedupSGD:bool = False, sgd_thresh:float = 0.1):
"""Easy initialization of RandomSearchDriver
Parameters
----------
- randomize_func : f(desvar_dict)
- Function to randomize desired variables of desvar_dict
maxiter : int, optional
Maximum iterations
max_time : int, optional
Maximum time in seconds
disp : bool
+
"""
- # self.T = T
+
self.learning_rate = learning_rate
self.upper = upper
self.lower = lower
@@ -2098,5 +2134,6 @@ class EasySGDDriver(SGDDriver, EasyDriverBase):
self.beta2 = beta2
self.gamma_min_factor = gamma_min_factor
self.gamma_min = gamma_min_factor # * learning_rate
+
SGDDriver.__init__(self, maxiter=maxiter, max_time=max_time, disp=disp, run_parallel=run_parallel,
speedupSGD=speedupSGD, sgd_thresh=sgd_thresh)