Merge pull request #300 from hiddenSymmetries/ml/coil_length_fix

Bug in coil examples and docs: missing "max" argument

docs/source/example_coils.rst

@@ -214,13 +214,19 @@ used and the corresponding terms::
   Jls = [CurveLength(c) for c in base_curves]
 
   # Form the total objective function.
-  objective = Jf + LENGTH_WEIGHT * QuadraticPenalty(sum(Jls), LENGTH_TARGET)
+  objective = Jf + LENGTH_WEIGHT * QuadraticPenalty(sum(Jls), LENGTH_TARGET, "max")
 
 In the last line, we have used the fact that the Optimizable objects
 representing the individual terms in the objective can be scaled by a
 constant and added.  (This feature applies to Optimizable objects that
 have a function ``J()`` returning the objective and, if gradients are
-used, a function ``dJ()`` returning the gradient.)
+used, a function ``dJ()`` returning the gradient.)  Also, the
+``"max"`` option to :obj:`~simsopt.objectives.QuadraticPenalty`
+specifies that the length penalty is active if the coil length is too
+large but not if it is too small. You can instead specify a penalty
+for values that are too small or a regular 2-sided quadratic penalty
+by setting the last argument to ``"min"`` or ``"identity"``
+respectively.
 
 You can check the degrees of freedom that will be varied in the
 optimization by printing the ``dof_names`` property of the objective::
@@ -359,6 +365,7 @@ many of which are illustrated in
 - :obj:`~simsopt.geo.CurveCurveDistance`: Useful for ensuring the minimum coil-to-coil distance is at least a specified target value.
 - :obj:`~simsopt.geo.CurveSurfaceDistance`: Useful for ensuring the minimum coil-to-plasma distance is at least a specified target value.
 - :obj:`~simsopt.geo.ArclengthVariation`: Ensures the curves are parameterized using (approximately) a uniform-arclength parameter.
+- :obj:`~simsopt.geo.LinkingNumber`: Prevents coils from becoming topologically linked to each other.
 
 You can click on any of the links above in this section to see the precise definitions of these objective terms.
 
@@ -449,7 +456,7 @@ the previous examples, we additionally define::
       + LENGTH_WEIGHT * sum(Jls) \
       + DISTANCE_WEIGHT * Jdist \
       + CURVATURE_WEIGHT * sum(Jcs) \
-      + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs) \
+      + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD, "max") for J in Jmscs) \
       + ARCLENGTH_WEIGHT * sum(Jals)
 
 As can be seen here, in the stochastic optimization method,
@@ -516,7 +523,7 @@ surface normal vector computed with the results of the virtual casing principle:
   # fact that Optimizable objects with J() and dJ() functions can be
   # multiplied by scalars and added:
   JF = Jf \
-      + LENGTH_PENALTY * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves)))
+      + LENGTH_PENALTY * sum(QuadraticPenalty(Jls[i], Jls[i].J(), "identity") for i in range(len(base_curves)))
 
 
 The example above uses very minimal coil regularization: only the deviation
@@ -609,5 +616,5 @@ Finally, the objective function takes the form::
   # fact that Optimizable objects with J() and dJ() functions can be
   # multiplied by scalars and added:
   JF = Jf \
-      + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves))) \
+      + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J(), "max") for i in range(len(base_curves))) \
       + DIST_PEN * Jdist

examples/1_Simple/stage_two_optimization_minimal.py

@@ -96,7 +96,7 @@
 # Form the total objective function. To do this, we can exploit the
 # fact that Optimizable objects with J() and dJ() functions can be
 # multiplied by scalars and added:
-JF = Jf + LENGTH_WEIGHT * QuadraticPenalty(sum(Jls), LENGTH_TARGET)
+JF = Jf + LENGTH_WEIGHT * QuadraticPenalty(sum(Jls), LENGTH_TARGET, "max")
 
 B_dot_n = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)
 print('Initial max|B dot n|:', np.max(np.abs(B_dot_n)))
@@ -142,5 +142,9 @@ def fun(dofs):
 B_dot_n = np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)
 print('Final max|B dot n|:', np.max(np.abs(B_dot_n)))
 
+total_curve_length = sum(func.J() for func in Jls)
+print("Sum of lengths of base curves:", total_curve_length)
+assert total_curve_length < 1.1 * LENGTH_TARGET
+
 # Save the optimized coil shapes and currents so they can be loaded into other scripts for analysis:
 bs.save(OUT_DIR + "biot_savart_opt.json")

examples/2_Intermediate/stage_two_optimization.py

@@ -124,7 +124,7 @@
     + CC_WEIGHT * Jccdist \
     + CS_WEIGHT * Jcsdist \
     + CURVATURE_WEIGHT * sum(Jcs) \
-    + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs)
+    + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD, "max") for J in Jmscs)
 
 # We don't have a general interface in SIMSOPT for optimisation problems that
 # are not in least-squares form, so we write a little wrapper function that we

examples/2_Intermediate/stage_two_optimization_finite_beta.py

@@ -121,7 +121,7 @@
 # fact that Optimizable objects with J() and dJ() functions can be
 # multiplied by scalars and added:
 JF = Jf \
-    + LENGTH_PENALTY * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves)))
+    + LENGTH_PENALTY * sum(QuadraticPenalty(Jls[i], Jls[i].J(), "identity") for i in range(len(base_curves)))
 
 # We don't have a general interface in SIMSOPT for optimisation problems that
 # are not in least-squares form, so we write a little wrapper function that we

examples/2_Intermediate/stage_two_optimization_stochastic.py

@@ -165,7 +165,7 @@ def pprint(*args, **kwargs):
     + LENGTH_WEIGHT * sum(Jls) \
     + DISTANCE_WEIGHT * Jdist \
     + CURVATURE_WEIGHT * sum(Jcs) \
-    + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs) \
+    + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD, "max") for J in Jmscs) \
     + ARCLENGTH_WEIGHT * sum(Jals)
 
 # We don't have a general interface in SIMSOPT for optimisation problems that

examples/3_Advanced/stage_two_optimization_finitebuild.py

@@ -132,7 +132,7 @@
 # fact that Optimizable objects with J() and dJ() functions can be
 # multiplied by scalars and added:
 JF = Jf \
-    + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves))) \
+    + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J(), "max") for i in range(len(base_curves))) \
     + DIST_PEN * Jdist
 
 # We don't have a general interface in SIMSOPT for optimisation problems that

src/simsopt/objectives/utilities.py

@@ -93,16 +93,16 @@ def dJ(self):
 
 class QuadraticPenalty(Optimizable):
 
-    def __init__(self, obj, cons=0., f="min"):
+    def __init__(self, obj, cons=0., f="identity"):
         r"""
         A quadratic penalty function of the form :math:`0.5f(\text{obj}.J() - \text{cons})^2` for an underlying objective ``obj``
-        and wrapping function `f`. This can be used to implement a barrier penalty function for (in)equality
-        constrained optimization problem. The wrapping function defaults to 'min'.
+        and wrapping function ``f``. This can be used to implement a barrier penalty function for (in)equality
+        constrained optimization problem. The wrapping function defaults to ``"identity"``.
 
         Args:
             obj: the underlying objective. It should provide a ``.J()`` and ``.dJ()`` function.
             cons: constant
-            f: the function that wraps the difference :math:`obj-\text{cons}`.  The options are 'min', 'max', or 'identity'
+            f: the function that wraps the difference :math:`obj-\text{cons}`.  The options are ``"min"``, ``"max"``, or ``"identity"``.
                which respectively return :math:`\min(\text{obj}-\text{cons}, 0)`, :math:`\max(\text{obj}-\text{cons}, 0)`, and :math:`\text{obj}-\text{cons}`.
         """
         Optimizable.__init__(self, x0=np.asarray([]), depends_on=[obj])