Cusp finding: implement bound consistency improvement

brush-strokes/src/lib/Math/Bezier/Stroke.hs

@@ -1636,7 +1636,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                    | otherwise
                    = ( [ d, u ], ( "s", s_mid ) )
              in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) )
-                  ( makeBox1Consistent =<< bisGuesses )
+                  ( doStrategy =<< bisGuesses )
       where
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
@@ -1655,6 +1655,65 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
                 !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
             in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
+    -- Attempting to implement Algorithm 6 "Heuristic to apply bound-consistency"
+    -- from the paper
+    -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+    doStrategy box@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+      | eps_box1 > eps_box2
+      = if | w_t > ebs_box1 || w_s > ebs_box1
+           -> makeBox1Consistent box
+           | w_t > ebs_box2 || w_s > ebs_box2
+           -> [ makeBox2Consistent box ]
+           | otherwise
+           -> [ box ]
+      | otherwise
+      = if | w_t > ebs_box2 || w_s > ebs_box2
+           -> [ makeBox2Consistent box ]
+           | w_t > ebs_box1 || w_s > ebs_box1
+           -> makeBox1Consistent box
+           | otherwise
+           -> [ box ]
+      where
+        eps_box1 = 0.4
+        eps_box2 = 0.1
+        w_t = t_hi - t_lo
+        w_s = s_hi - s_lo
+
+    -- An implementation of "bound-consistency" from the paper
+    -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+    makeBox2Consistent :: Box -> Box
+    makeBox2Consistent x = ( `State.evalState` False ) $ doLoop 0.25 x
+      where
+        doLoop :: Double -> Box -> State Bool Box
+        doLoop lambda x = do
+          x'  <- boundConsistency get_t set_t lambda x
+          x'' <- boundConsistency get_s set_s lambda x'
+          modified <- State.get
+          let lambda' = if modified then lambda else 0.5 * lambda
+          if lambda' < 0.001
+          then return x''
+          else do { State.put False ; doLoop lambda' x'' }
+
+        boundConsistency :: ( Box -> 𝕀 Double )
+                         -> ( 𝕀 Double -> Box -> Box )
+                         -> Double -> Box -> State Bool Box
+        boundConsistency getter setter lambda box = do
+          let x@( 𝕀 x_inf x_sup ) = getter box
+              c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
+              c2 = lambda * x_inf + ( 1 - lambda ) * x_sup
+              x'_inf =
+                case makeBox1Consistent ( setter ( 𝕀 x_inf c1 ) box ) of
+                  []  -> c1
+                  x's -> minimum $ map ( inf . getter ) x's
+              x'_sup =
+                case makeBox1Consistent ( setter ( 𝕀 c2 x_sup ) box ) of
+                  []  -> c2
+                  x's -> maximum $ map ( sup . getter ) x's
+              x' = 𝕀 x'_inf x'_sup
+          when ( width x - width x' >= eps_eq ) $
+            State.put True
+          return $ setter x' box
+
     -- An implementation of "bc_enforce" from the paper
     -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
     --
@@ -1863,10 +1922,6 @@ allNarrowingOperators eps_eq eps_bis eqs =
       [ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ]
       ]
   where
-    get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup
-    set_t ( 𝕀 t_inf t_sup ) ( _, i, s ) = ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), i, s )
-    get_s ( _ , _, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) = 𝕀 s_inf s_sup
-    set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )
     ff'_t (t, i, s) =
       let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup )
               , _D22_dx = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) )
@@ -1897,3 +1952,8 @@ allNarrowingOperators eps_eq eps_bis eqs =
               , _D12_dy = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) )
               } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
       in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )
+
+get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup
+set_t ( 𝕀 t_inf t_sup ) ( _, i, s ) = ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), i, s )
+get_s ( _ , _, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) = 𝕀 s_inf s_sup
+set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )