{- ORMOLU_DISABLE -}
-- Implicit CAD. Copyright (C) 2011, Christopher Olah (chris@colah.ca)
-- Copyright (C) 2016, Julia Longtin (julial@turinglace.com)
-- Released under the GNU AGPLV3+, see LICENSE

module Graphics.Implicit.Export.Render.Interpolate (interpolate) where

import Prelude((*), (>), (<), (/=), (+), (-), (/), (==), (&&), abs)

import Graphics.Implicit.Definitions (ℝ, Fastℕ, ℝ2)
import Linear (V2(V2))

default (Fastℕ, ℝ)
-- Consider a function f(x):

{-
   |   \        f(x)
   |    - \
   |_______\________ x
            |
             \
-}

-- The purpose of interpolate is to find the value of x where f(x) crosses zero.
-- This should be accomplished cheaply and accurately.

-- We are given the constraint that x will be between a and b.

-- We are also given the values of f at a and b: aval and bval.

-- Additionaly, we get f (continuous and differentiable almost everywhere),
-- and the resolution of the object (so that we can make decisions about
-- how precise we need to be).

-- While the output will never be used, interpolate will be called
-- in cases where f(x) doesn't cross zero (ie. aval and bval are both
-- positive or negative.

-- Clarification: If f(x) crosses zero, but doesn't necessarily have
-- to do so by intermediate value theorem, it is beyond the scope of this
-- function.

-- If it doesn't cross zero, we don't actually care what answer we give,
-- just that it's cheap.

-- FIXME: accept resolution on multiple axises.

interpolate :: ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ -> ℝ
interpolate :: ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ -> ℝ
interpolate (V2 ℝ
a ℝ
aval) (V2 ℝ
_ ℝ
bval) ℝ -> ℝ
_ ℝ
_ | ℝ
avalforall a. Num a => a -> a -> a
*ℝ
bval forall a. Ord a => a -> a -> Bool
> ℝ
0 = ℝ
a

-- The obvious:
interpolate (V2 ℝ
a  ℝ
0) ℝ2
_ ℝ -> ℝ
_ ℝ
_  = ℝ
a
interpolate ℝ2
_ (V2 ℝ
b  ℝ
0) ℝ -> ℝ
_ ℝ
_  = ℝ
b

-- It may seem, at first, that our task is trivial.
-- Just use linear interpolation!
-- Unfortunately, there's a nasty failure case

{-                   /
                    /
  ________#________/____
  ________________/
-}

-- This is really common for us, for example in cubes,
-- where another variable dominates.

-- It may even be the case that, because we are so close
-- to the side, it looks like we are really close to an
-- answer... And we just give it back.

-- So we need to detect this. And get free goodies while we're
-- at it (shrink domain to guess within fromm (a,b) to (a',b'))
-- :)

interpolate (V2 ℝ
a ℝ
aval) (V2 ℝ
b ℝ
bval) ℝ -> ℝ
f ℝ
_ =
    -- Make sure aval > bval, then pass to interpolateLin
    if ℝ
aval forall a. Ord a => a -> a -> Bool
> ℝ
bval
    then Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateLin Fastℕ
0 (forall a. a -> a -> V2 a
V2 ℝ
a ℝ
aval) (forall a. a -> a -> V2 a
V2 ℝ
b ℝ
bval) ℝ -> ℝ
f
    else Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateLin Fastℕ
0 (forall a. a -> a -> V2 a
V2 ℝ
b ℝ
bval) (forall a. a -> a -> V2 a
V2 ℝ
a ℝ
aval) ℝ -> ℝ
f

-- Yay, linear interpolation!

-- Try the answer linear interpolation gives us...
-- (n is to cut us off if recursion goes too deep)
interpolateLin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateLin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateLin Fastℕ
n (V2 ℝ
a ℝ
aval) (V2 ℝ
b ℝ
bval) ℝ -> ℝ
obj | ℝ
aval forall a. Eq a => a -> a -> Bool
/= ℝ
bval=
    let
        -- Interpolate and evaluate
        mid :: ℝ
        mid :: ℝ
mid = ℝ
a forall a. Num a => a -> a -> a
+ (ℝ
bforall a. Num a => a -> a -> a
-ℝ
a)forall a. Num a => a -> a -> a
*ℝ
avalforall a. Fractional a => a -> a -> a
/(ℝ
avalforall a. Num a => a -> a -> a
-ℝ
bval)
        midval :: ℝ
midval = ℝ -> ℝ
obj ℝ
mid
    -- Are we done?
    in if ℝ
midval forall a. Eq a => a -> a -> Bool
== ℝ
0
    then ℝ
mid
    --
    else let
        (ℝ
a', ℝ
a'val, ℝ
b', ℝ
b'val, ℝ
improveRatio) =
            if ℝ
midval forall a. Ord a => a -> a -> Bool
> ℝ
0
                then (ℝ
mid, ℝ
midval, ℝ
b, ℝ
bval, ℝ
midvalforall a. Fractional a => a -> a -> a
/ℝ
aval)
                else (ℝ
a, ℝ
aval, ℝ
mid, ℝ
midval, ℝ
midvalforall a. Fractional a => a -> a -> a
/ℝ
bval)

    -- some times linear interpolate doesn't work,
    -- because one side is very close to zero and flat
    -- we catch it because the interval won't shrink when
    -- this is the case. To test this, we look at whether
    -- the replaced point evaluates to substantially closer
    -- to zero than the previous one.
    in if ℝ
improveRatio forall a. Ord a => a -> a -> Bool
< ℝ
0.3 Bool -> Bool -> Bool
&& Fastℕ
n forall a. Ord a => a -> a -> Bool
< Fastℕ
4
    -- And we continue on.
    then Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateLin (Fastℕ
nforall a. Num a => a -> a -> a
+Fastℕ
1) (forall a. a -> a -> V2 a
V2 ℝ
a' ℝ
a'val) (forall a. a -> a -> V2 a
V2 ℝ
b' ℝ
b'val) ℝ -> ℝ
obj
    -- But if not, we switch to binary interpolate, which is
    -- immune to this problem
    else Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateBin (Fastℕ
nforall a. Num a => a -> a -> a
+Fastℕ
1) (forall a. a -> a -> V2 a
V2 ℝ
a' ℝ
a'val) (forall a. a -> a -> V2 a
V2 ℝ
b' ℝ
b'val) ℝ -> ℝ
obj

-- And a fallback:
interpolateLin Fastℕ
_ (V2 ℝ
a  ℝ
_) ℝ2
_ ℝ -> ℝ
_ = ℝ
a

-- Now for binary searching!
interpolateBin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ

-- The termination case:

interpolateBin :: Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateBin Fastℕ
5 (V2 ℝ
a ℝ
aval) (V2 ℝ
b ℝ
bval) ℝ -> ℝ
_ =
    if forall a. Num a => a -> a
abs ℝ
aval forall a. Ord a => a -> a -> Bool
< forall a. Num a => a -> a
abs ℝ
bval
    then ℝ
a
    else ℝ
b

-- Otherwise, have fun with mid!

interpolateBin Fastℕ
n (V2 ℝ
a ℝ
aval) (V2 ℝ
b ℝ
bval) ℝ -> ℝ
f =
    let
        mid :: ℝ
        mid :: ℝ
mid = (ℝ
aforall a. Num a => a -> a -> a
+ℝ
b)forall a. Fractional a => a -> a -> a
/ℝ
2
        midval :: ℝ
midval = ℝ -> ℝ
f ℝ
mid
    in if ℝ
midval forall a. Ord a => a -> a -> Bool
> ℝ
0
    then Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateBin (Fastℕ
nforall a. Num a => a -> a -> a
+Fastℕ
1) (forall a. a -> a -> V2 a
V2 ℝ
mid ℝ
midval) (forall a. a -> a -> V2 a
V2 ℝ
b ℝ
bval) ℝ -> ℝ
f
    else Fastℕ -> ℝ2 -> ℝ2 -> (ℝ -> ℝ) -> ℝ
interpolateBin (Fastℕ
nforall a. Num a => a -> a -> a
+Fastℕ
1) (forall a. a -> a -> V2 a
V2 ℝ
a ℝ
aval) (forall a. a -> a -> V2 a
V2 ℝ
mid ℝ
midval) ℝ -> ℝ
f