{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE DefaultSignatures #-}
--------------------------------------------------------------------------------
-- |
-- Module      :  HGeometry.HyperPlane.Class
-- Copyright   :  (C) Frank Staals
-- License     :  see the LICENSE file
-- Maintainer  :  Frank Staals
--
-- Classes for representing hyperplanes in d-dimensional space.
--
--------------------------------------------------------------------------------
module HGeometry.HyperPlane.Class
  ( HyperPlane_(..)
  , ConstructableHyperPlane_(..)
  , NonVerticalHyperPlane_(..)
  , Plane_, pattern Plane_
  , isParallelTo
  , pointOn

  , HyperPlaneFromPoints(..)
  , showPlaneEquation
  ) where

import Control.Lens hiding (snoc, cons, uncons, unsnoc)
import Data.Default
import Data.Kind (Constraint)
import Data.Type.Ord
import GHC.TypeNats
import HGeometry.Ext
import HGeometry.Point
import HGeometry.Properties
import HGeometry.Vector
import Prelude hiding (head, last)

--------------------------------------------------------------------------------

-- | myLine, myHyperPlane, and myNVHyperPlane all represent the same
-- 2d hyperplane; i.e. a line in R^2.

-- $setup
-- >>> import HGeometry.HyperPlane.NonVertical
-- >>> import HGeometry.HyperPlane
-- >>> import HGeometry.Line.LineEQ
--
-- >>> let myLine      = LineEQ 1 2                          :: LineEQ Double
-- >>> let myLineAgain = HyperPlane (Vector3 2 1 (-1))       :: HyperPlane 2 Double
-- >>> let myLineAsNV  = NonVerticalHyperPlane (Vector2 1 2) :: NonVerticalHyperPlane 2 Double
-- >>> let myVerticalLine = HyperPlane (Vector3 5 (-1) 0)    :: HyperPlane 2 Double
--
-- >>> let myOtherLine   = HyperPlane2 4 3 2                   :: HyperPlane 2 Double
-- >>> let myOtherNVLine = NonVerticalHyperPlane (Vector2 (-3/2 :: Double) (-2))
-- >>> let myPlane       = HyperPlane3 10 2 3 (-1)             :: HyperPlane 3 Double
--
--
-- 'myLine', 'myLineAgain', and 'myLineAsNV' all represent the line y = 1*x + 2
--
-- 'myVerticalLine' represents the vertical line x=5
--
--
-- 'myOtherLine' represents the line 4 + 3*x + 2*y = 0, in other words y = -(3/2)*x - 2
--
-- 'myPlane' represents the plane z = 2*x + 3*y + 10


-- | A class to represent hyperplanes in d-dimensional space.
class ( NumType hyperPlane ~ r
      , Dimension hyperPlane ~ d
      , Has_ Vector_ d r
      , Has_ Vector_ (1+d) r
      -- , NumType (EquationFor hyperPlane) ~ r
      ) => HyperPlane_ hyperPlane d r | hyperPlane -> d
                                      , hyperPlane -> r where
--  {-# MINIMAL hyperPlaneTrough #-}

  -- | Evalute the expression \(a_0 + \sum_i=1^d a_i*p_i \) at the given point.
  --
  --
  -- >>> evalHyperPlaneEquation myLine (Point2 5 10)
  -- -3.0
  evalHyperPlaneEquation     :: (Num r, Point_ point d r) => hyperPlane -> point -> r
  default evalHyperPlaneEquation :: (Num r, Point_ point d r, Has_ Metric_ (d+1) r)
                                 => hyperPlane -> point -> r
  evalHyperPlaneEquation hyperPlane
h point
p = hyperPlane -> Vector (d + 1) r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r) =>
hyperPlane -> Vector (d + 1) r
hyperPlaneEquation hyperPlane
h Vector (d + 1) r -> Vector (d + 1) r -> r
forall vector (d :: Natural) r.
(Metric_ vector d r, Num r) =>
vector -> vector -> r
`dot` (r -> Vector d r -> Vector (d + 1) r
forall vector' vector (d :: Natural) r.
(Vector_ vector d r, Vector_ vector' (d + 1) r) =>
r -> vector -> vector'
cons r
1 (Vector d r -> Vector (d + 1) r) -> Vector d r -> Vector (d + 1) r
forall a b. (a -> b) -> a -> b
$ point
ppoint -> Getting (Vector d r) point (Vector d r) -> Vector d r
forall s a. s -> Getting a s a -> a
^.Getting (Vector d r) point (Vector d r)
forall (d :: Natural) r s.
(Dimension point ~ d, NumType point ~ r, Dimension point ~ d,
 NumType point ~ s) =>
Lens point point (Vector d r) (Vector d s)
forall point point' (d :: Natural) r s.
(HasVector point point', Dimension point ~ d, NumType point ~ r,
 Dimension point' ~ d, NumType point' ~ s) =>
Lens point point' (Vector d r) (Vector d s)
Lens point point (Vector d r) (Vector d r)
vector)
  {-# INLINE evalHyperPlaneEquation #-}


  -- -- | Constructs the hyperplane through d points
  -- hyperPlaneTrough :: Num r => Point_ point d r => Vector d point -> hyperPlane

  -- | A hyperplane \(h) has coefficients \(a_i \in \mathbb{R}\) so that
  -- a \point \( (p_1,..,p_d) \) lies on \(h) iff:
  -- \( a_0  + \sum_i=1^d a_i*p_i = 0 \)
  --
  -- this fuction returns the vector of these coefficients \(\langle a_0,..,a_d \rangle\).
  --
  -- >>> hyperPlaneEquation myLine
  -- Vector3 2.0 1.0 (-1.0)
  -- >>> hyperPlaneEquation myLineAgain
  -- Vector3 2.0 1.0 (-1.0)
  -- >>> hyperPlaneEquation myLineAsNV
  -- Vector3 2.0 1.0 (-1.0)
  --
  --
  -- so for example myLineAsNV actually a line with slope 1 through the point (0,2)
  --
  hyperPlaneEquation :: ( Num r
                        ) => hyperPlane -> Vector (d+1) r
  default hyperPlaneEquation :: ( NonVerticalHyperPlane_ hyperPlane d r
                                , Num r
                                , 1 <= d
                                , Has_ Vector_ (d+1) r, KnownNat (d-1)
                                )
                             => hyperPlane -> Vector (d+1) r
  hyperPlaneEquation hyperPlane
h = r -> Vector d r -> Vector (d + 1) r
forall vector' vector (d :: Natural) r.
(Vector_ vector d r, Vector_ vector' (d + 1) r) =>
r -> vector -> vector'
cons r
a0 Vector d r
a
    where
      a' :: Vector d r
a' = hyperPlane
hhyperPlane
-> Getting (Vector d r) hyperPlane (Vector d r) -> Vector d r
forall s a. s -> Getting a s a -> a
^.Getting (Vector d r) hyperPlane (Vector d r)
forall hyperPlane (d :: Natural) r.
NonVerticalHyperPlane_ hyperPlane d r =>
Lens' hyperPlane (Vector d r)
Lens' hyperPlane (Vector d r)
hyperPlaneCoefficients
      a :: Vector d r
a  = Vector d r
a'Vector d r -> (Vector d r -> Vector d r) -> Vector d r
forall a b. a -> (a -> b) -> b
&(r -> Identity r) -> Vector d r -> Identity (Vector d r)
forall vector (d :: Natural) r.
(Vector_ vector d r, 1 <= d, KnownNat (d - 1)) =>
IndexedLens' Int vector r
IndexedLens' Int (Vector d r) r
last ((r -> Identity r) -> Vector d r -> Identity (Vector d r))
-> r -> Vector d r -> Vector d r
forall s t a b. ASetter s t a b -> b -> s -> t
.~ -r
1
      a0 :: r
a0 = Vector d r
a'Vector d r -> Getting r (Vector d r) r -> r
forall s a. s -> Getting a s a -> a
^.Getting r (Vector d r) r
forall vector (d :: Natural) r.
(Vector_ vector d r, 1 <= d, KnownNat (d - 1)) =>
IndexedLens' Int vector r
IndexedLens' Int (Vector d r) r
last
  {-# INLINE hyperPlaneEquation #-}

  -- | Get the normal vector of the hyperplane. The vector points into
  -- the positive halfspace.
  --
  -- >>> normalVector myVerticalLine
  -- Vector2 1.0 (-0.0)
  -- >>> normalVector myLine
  -- Vector2 (-1.0) 1.0
  normalVector :: (Num r, Eq r, 1 <= d) => hyperPlane -> Vector d r
  default normalVector :: (KnownNat d, Num r, Eq r, 1 <= d, Has_ Vector_ (d+1) r, d <= d+1)
                       => hyperPlane -> Vector d r
  normalVector hyperPlane
h = Vector d r -> Vector d r
forall r vector (d :: Natural).
(Num r, Vector_ vector d r) =>
vector -> vector
negated (Vector d r -> Vector d r)
-> (Vector (d + 1) r -> Vector d r)
-> Vector (d + 1) r
-> Vector d r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector (d + 1) r -> Vector d r
forall (i :: Natural) (d :: Natural) vector vector' r.
(i <= d, Vector_ vector d r, Vector_ vector' i r) =>
vector -> vector'
suffix (Vector (d + 1) r -> Vector d r) -> Vector (d + 1) r -> Vector d r
forall a b. (a -> b) -> a -> b
$ hyperPlane -> Vector (d + 1) r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r) =>
hyperPlane -> Vector (d + 1) r
hyperPlaneEquation hyperPlane
h
  {-# INLINE normalVector #-}
  -- https://en.wikipedia.org/wiki/Normal_(geometry)#Hypersurfaces_in_n-dimensional_space
  -- states that: if the hyperplane is defined as the solution set of a single linear equation
  -- a1*x1 + .. + an*xn = c
  -- then the vector n = [a1,.., an] is a normal.


  -- | Test if a point lies on a hyperplane.
  --
  -- >>> Point2 0 2 `onHyperPlane` myLineAgain
  -- True
  -- >>> Point2 1 3 `onHyperPlane` myLineAgain
  -- True
  -- >>> Point2 1 5 `onHyperPlane` myLineAgain
  -- False
  --
  -- >>> Point2 0 2 `onHyperPlane` myLineAsNV
  -- True
  -- >>> Point2 1 3 `onHyperPlane` myLineAsNV
  -- True
  -- >>> Point2 1 5 `onHyperPlane` myLineAsNV
  -- False
  onHyperPlane     :: (Point_ point d r, Eq r, Num r) => point -> hyperPlane -> Bool
  default onHyperPlane :: ( Point_ point d r, Eq r, Num r
                          ) => point -> hyperPlane -> Bool
  point
q `onHyperPlane` hyperPlane
h = (r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0) (r -> Bool) -> r -> Bool
forall a b. (a -> b) -> a -> b
$ hyperPlane -> point -> r
forall point. (Num r, Point_ point d r) => hyperPlane -> point -> r
forall hyperPlane (d :: Natural) r point.
(HyperPlane_ hyperPlane d r, Num r, Point_ point d r) =>
hyperPlane -> point -> r
evalHyperPlaneEquation hyperPlane
h point
q
  {-# INLINE onHyperPlane #-}

  -- | Test if a point lies on a hyperplane. Returns the sign when evaluating the
  -- hyperplane equation.
  --
  -- >>> Point2 0 2 `onSideTest` myLineAgain
  -- EQ
  -- >>> Point2 1 3 `onSideTest` myLineAgain
  -- EQ
  -- >>> Point2 1 5 `onSideTest` myLineAgain
  -- GT
  -- >>> Point2 4 5 `onSideTest` myLineAgain
  -- LT
  -- >>> Point2 0 0 `onSideTest` HyperPlane2 1 (-1) 0
  -- LT
  --
  -- >>> Point2 1 1 `onSideTest` myVerticalLine
  -- LT
  -- >>> Point2 10 1 `onSideTest` myVerticalLine
  -- GT
  -- >>> Point2 5 20 `onSideTest` myVerticalLine
  -- EQ
  --
  -- >>> Point2 0 1 `onSideTest` myOtherLine
  -- LT
  -- >>> Point2 0 (-2) `onSideTest` myOtherLine
  -- EQ
  -- >>> Point2 1 (-3.5) `onSideTest` myOtherLine
  -- EQ
  -- >>> Point2 1 (-4) `onSideTest` myOtherLine
  -- GT
  onSideTest     :: (Point_ point d r, Ord r, Num r) => point -> hyperPlane -> Ordering
  onSideTest point
q hyperPlane
h = r
0 r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` hyperPlane -> point -> r
forall point. (Num r, Point_ point d r) => hyperPlane -> point -> r
forall hyperPlane (d :: Natural) r point.
(HyperPlane_ hyperPlane d r, Num r, Point_ point d r) =>
hyperPlane -> point -> r
evalHyperPlaneEquation hyperPlane
h point
q
  {-# INLINE onSideTest #-}

-- | Produce a point that lies on the hyperplane. No gurantees are given about which point
--
-- >>> pointOn myLine
-- Point2 (-2.0) 0.0
pointOn   :: forall hyperPlane d r.
             ( HyperPlane_ hyperPlane d r, Eq r, Fractional r, Has_ Additive_ d r
             , FoldableWithIndex Int (Vector d)
             , Has_ Vector_ (d+1) r, d <= d +1, 0 <= (d+1)-1 -- these are silly :(
             )
          => hyperPlane -> Point d r
pointOn :: forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Eq r, Fractional r,
 Has_ Additive_ d r, FoldableWithIndex Int (Vector d),
 Has_ Vector_ (d + 1) r, d <= (d + 1), 0 <= ((d + 1) - 1)) =>
hyperPlane -> Point d r
pointOn hyperPlane
h = case Vector (d + 1) r -> (r, Vector d r)
forall vector' vector (d :: Natural) r.
(Vector_ vector (d + 1) r, Vector_ vector' d r, 0 <= ((d + 1) - 1),
 d <= Dimension vector) =>
vector -> (r, vector')
uncons (Vector (d + 1) r -> (r, Vector d r))
-> Vector (d + 1) r -> (r, Vector d r)
forall a b. (a -> b) -> a -> b
$ hyperPlane -> Vector (d + 1) r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r) =>
hyperPlane -> Vector (d + 1) r
hyperPlaneEquation hyperPlane
h of
              (r
0, Vector d r
_)  -> Point d r
forall point (d :: Natural) r.
(Num r, ConstructablePoint_ point d r) =>
point
origin
              (r
a0, Vector d r
a) -> case (Int -> r -> Bool) -> Vector d r -> Maybe (Int, r)
forall i (f :: * -> *) a.
FoldableWithIndex i f =>
(i -> a -> Bool) -> f a -> Maybe (i, a)
ifind ((r -> Bool) -> Int -> r -> Bool
forall a b. a -> b -> a
const (r -> r -> Bool
forall a. Eq a => a -> a -> Bool
/= r
0)) (Vector d r
a :: Vector d r) of
                           Maybe (Int, r)
Nothing     -> [Char] -> Point d r
forall a. HasCallStack => [Char] -> a
error [Char]
"pointOn: Invalid hyperplane"
                           Just (Int
i,r
ai) ->
                             (Point d r
forall point (d :: Natural) r.
(Num r, ConstructablePoint_ point d r) =>
point
origin :: Point d r)Point d r -> (Point d r -> Point d r) -> Point d r
forall a b. a -> (a -> b) -> b
&(Vector d r -> Identity (Vector d r))
-> Point d r -> Identity (Point d r)
forall (d :: Natural) r s.
(Dimension (Point d r) ~ d, NumType (Point d r) ~ r,
 Dimension (Point d r) ~ d, NumType (Point d r) ~ s) =>
Lens (Point d r) (Point d r) (Vector d r) (Vector d s)
forall point point' (d :: Natural) r s.
(HasVector point point', Dimension point ~ d, NumType point ~ r,
 Dimension point' ~ d, NumType point' ~ s) =>
Lens point point' (Vector d r) (Vector d s)
Lens (Point d r) (Point d r) (Vector d r) (Vector d r)
vector((Vector d r -> Identity (Vector d r))
 -> Point d r -> Identity (Point d r))
-> ((r -> Identity r) -> Vector d r -> Identity (Vector d r))
-> (r -> Identity r)
-> Point d r
-> Identity (Point d r)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Int -> IndexedTraversal' Int (Vector d r) r
forall vector (d :: Natural) r.
Vector_ vector d r =>
Int -> IndexedTraversal' Int vector r
component' Int
i ((r -> Identity r) -> Point d r -> Identity (Point d r))
-> r -> Point d r -> Point d r
forall s t a b. ASetter s t a b -> b -> s -> t
.~ (r -> r
forall a. Num a => a -> a
negate r
a0 r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
ai)
  -- We are trying to find a point p so that a0 + sum_{i=1}^d ai*pi = 0,
  -- in other words, so that
  --
  --    sum_{i=1}^d ai*pi = -a0                          (1)
  --
  -- so if a0 is actually zero, we can simply choose all the pi's in Eq 1 to be zero anyway.
  -- if a0 is not zero, then there must be at least one ai that is non-zero; otherwise
  -- we would have the equation 0 = <notzero>. Hence, we find this ai. We then rewrite Eq (1)
  -- to: ai*pi + sum_{j /= i} aj*pj = -a0. Hence, we pick all the pj's to be zero, and
  -- set pi to -a0/ai.


-- | Class representing hyperplanes with methods to construct the
-- hyperplane.
class HyperPlane_ hyperPlane d r
     => ConstructableHyperPlane_ hyperPlane d r where

  -- | Additional constraints for constructing a hyperplane from its equation.
  type HyperPlaneFromEquationConstraint hyperPlane d r :: Constraint
  type HyperPlaneFromEquationConstraint hyperPlane d r = ()

  -- | Given the coefficients \(a_0,..,a_d\) of the equation, i.e.
  -- so that
  --
  -- \( a_0  + \sum_i=1^d a_i*p_i = 0 \)
  --
  -- construct the hyperplane form it.
  --
  -- >>> hyperPlaneFromEquation (Vector3 10 2 1) :: HyperPlane 2 Int -- the line 2*x + 1*y + 10 = 0
  -- HyperPlane [10,2,1]
  -- >>> hyperPlaneFromEquation $ Vector4 100 5 3 (-1)  :: HyperPlane 3 Int -- the plane 5*x + 3*y + (-1)*z + 100= 0
  -- HyperPlane [100,5,3,-1]
  --
  -- >>> myOtherLine == hyperPlaneFromEquation (Vector3 4 3 2)
  -- True
  hyperPlaneFromEquation :: HyperPlaneFromEquationConstraint hyperPlane d r
                         => Vector (d+1) r -> hyperPlane


  -- | Construct a Hyperplane from a point and a normal. The normal points into the halfplane
  -- for which the side-test is positive.
  --
  -- >>> myVerticalLine == fromPointAndNormal (Point2 5 30) (Vector2 1 0)
  -- True
  fromPointAndNormal     :: ( Point_ point d r, Num r)
                         => point -> Vector d r -> hyperPlane
  default fromPointAndNormal :: ( Point_ point d r, Num r
                                , ConstructableHyperPlane_ hyperPlane d r
                                , HyperPlaneFromEquationConstraint hyperPlane d r
                                , Has_ Metric_ d r
                                , Has_ Vector_ (d+1) r -- this is implied but whatever
                                )
                             => point -> Vector d r -> hyperPlane
  fromPointAndNormal point
q Vector d r
n = Vector (d + 1) r -> hyperPlane
forall hyperPlane (d :: Natural) r.
(ConstructableHyperPlane_ hyperPlane d r,
 HyperPlaneFromEquationConstraint hyperPlane d r) =>
Vector (d + 1) r -> hyperPlane
hyperPlaneFromEquation (Vector (d + 1) r -> hyperPlane) -> Vector (d + 1) r -> hyperPlane
forall a b. (a -> b) -> a -> b
$ r -> Vector d r -> Vector (d + 1) r
forall vector' vector (d :: Natural) r.
(Vector_ vector d r, Vector_ vector' (d + 1) r) =>
r -> vector -> vector'
cons r
a0 (Vector d r -> Vector d r
forall r vector (d :: Natural).
(Num r, Vector_ vector d r) =>
vector -> vector
negated Vector d r
n)
    where
      a0 :: r
a0 = (point
qpoint -> Getting (Vector d r) point (Vector d r) -> Vector d r
forall s a. s -> Getting a s a -> a
^.Getting (Vector d r) point (Vector d r)
forall (d :: Natural) r s.
(Dimension point ~ d, NumType point ~ r, Dimension point ~ d,
 NumType point ~ s) =>
Lens point point (Vector d r) (Vector d s)
forall point point' (d :: Natural) r s.
(HasVector point point', Dimension point ~ d, NumType point ~ r,
 Dimension point' ~ d, NumType point' ~ s) =>
Lens point point' (Vector d r) (Vector d s)
Lens point point (Vector d r) (Vector d r)
vector) Vector d r -> Vector d r -> r
forall vector (d :: Natural) r.
(Metric_ vector d r, Num r) =>
vector -> vector -> r
`dot` Vector d r
n
  {-# INLINE fromPointAndNormal #-}


-- -- | returns a point that lies on the hyperpalne
-- pointOn :: hyperPlane -> Point d r


--------------------------------------------------------------------------------

-- | Non-vertical hyperplanes.
class HyperPlane_ hyperPlane d r => NonVerticalHyperPlane_ hyperPlane d r where
  {-# MINIMAL hyperPlaneCoefficients #-}

  -- | Get the coordinate in dimension \(d\) of the hyperplane at the given position.
  --
  -- >>> evalAt (Point1 1) myLineAsNV
  -- 3.0
  -- >>> evalAt (Point1 10) myLineAsNV
  -- 12.0
  -- >>> evalAt (Point1 5) <$> asNonVerticalHyperPlane myOtherLine
  -- Just (-9.5)
  evalAt     :: ( Num r
                , 1 <= d
                , Point_ point (d-1) r
                ) => point -> hyperPlane -> r
  default evalAt :: ( Num r
                    , 1 <= d
                    , Point_ point (d-1) r
                    , Has_ Metric_ d r
                    , d - 1 <= d, KnownNat (d-1), (d-1) + 1 ~ d
                    ) => point -> hyperPlane -> r
  evalAt point
p hyperPlane
h = (hyperPlane
hhyperPlane
-> Getting (Vector d r) hyperPlane (Vector d r) -> Vector d r
forall s a. s -> Getting a s a -> a
^.Getting (Vector d r) hyperPlane (Vector d r)
forall hyperPlane (d :: Natural) r.
NonVerticalHyperPlane_ hyperPlane d r =>
Lens' hyperPlane (Vector d r)
Lens' hyperPlane (Vector d r)
hyperPlaneCoefficients) Vector d r -> Vector d r -> r
forall vector (d :: Natural) r.
(Metric_ vector d r, Num r) =>
vector -> vector -> r
`dot` Vector (Dimension point) r -> r -> Vector d r
forall vector' vector (d :: Natural) r.
(Vector_ vector d r, Vector_ vector' (d + 1) r) =>
vector -> r -> vector'
snoc (point
ppoint
-> Getting
     (Vector (Dimension point) r) point (Vector (Dimension point) r)
-> Vector (Dimension point) r
forall s a. s -> Getting a s a -> a
^.Getting
  (Vector (Dimension point) r) point (Vector (Dimension point) r)
forall (d :: Natural) r s.
(Dimension point ~ d, NumType point ~ r, Dimension point ~ d,
 NumType point ~ s) =>
Lens point point (Vector d r) (Vector d s)
forall point point' (d :: Natural) r s.
(HasVector point point', Dimension point ~ d, NumType point ~ r,
 Dimension point' ~ d, NumType point' ~ s) =>
Lens point point' (Vector d r) (Vector d s)
Lens
  point
  point
  (Vector (Dimension point) r)
  (Vector (Dimension point) r)
vector) r
1
  {-# INLINE evalAt #-}


  -- | The coefficients \( \langle a_1,..,a_d \rangle \) such that
  -- a point \(p = (p_1,..,p_d) \) lies on the hyperplane the given coefficients iff
  --
  --
  --  \( a_d + \sum_i=1^{d-1} a_i*p_i = p_d \)
  --
  -- >>> view hyperPlaneCoefficients myLineAsNV
  -- Vector2 1.0 2.0
  -- >>> view hyperPlaneCoefficients myLine
  -- Vector2 1.0 2.0
  -- >>> view hyperPlaneCoefficients $ Plane 1 2 3
  -- Vector3 1 2 3
  -- >>> asNonVerticalHyperPlane myOtherLine ^?_Just.hyperPlaneCoefficients
  -- Just (Vector2 (-1.5) (-2.0))
  hyperPlaneCoefficients :: Lens' hyperPlane (Vector d r)

  -- | Test if a point q lies above a non-vertical hyperplane h (i.e. verticalSideTest q h
  -- == GT), on the hyperplane (== EQ), or below (LT).
  --
  -- >>> Point2 0 2 `verticalSideTest` myLineAsNV
  -- EQ
  -- >>> Point2 1 3 `verticalSideTest` myLineAsNV
  -- EQ
  -- >>> Point2 1 5 `verticalSideTest` myLineAsNV
  -- GT
  -- >>> Point2 4 5 `verticalSideTest` myLineAsNV
  -- LT
  --
  -- >>> Point2 0 1 `verticalSideTest` myOtherNVLine
  -- GT
  -- >>> Point2 0 (-2) `verticalSideTest` myOtherNVLine
  -- EQ
  -- >>> Point2 1 (-3.5) `verticalSideTest` myOtherNVLine
  -- EQ
  -- >>> Point2 1 (-4) `verticalSideTest` myOtherNVLine
  -- LT
  verticalSideTest   :: (1 <= d, Point_ point d r, Ord r, Num r)
                     => point -> hyperPlane -> Ordering
  default verticalSideTest :: ( Point_ point d r
                              , Num r, Ord r
                              , 1 <= d
                              , (d-1) + 1 ~ d, d - 1 <= d, KnownNat (d-1) -- silly silly agian :(
                              , Has_ Metric_ d r
                              , Has_ Additive_ (d-1) r
                              ) => point -> hyperPlane -> Ordering
  verticalSideTest point
q hyperPlane
h = (point
qpoint -> Getting r point r -> r
forall s a. s -> Getting a s a -> a
^.Getting r point r
forall point (d :: Natural) r.
(1 <= d, Point_ point d r) =>
IndexedLens' Int point r
IndexedLens' Int point r
dCoord) r -> r -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Point (d - 1) r -> hyperPlane -> r
forall point.
(Num r, 1 <= d, Point_ point (d - 1) r) =>
point -> hyperPlane -> r
forall hyperPlane (d :: Natural) r point.
(NonVerticalHyperPlane_ hyperPlane d r, Num r, 1 <= d,
 Point_ point (d - 1) r) =>
point -> hyperPlane -> r
evalAt (point -> Point (d - 1) r
forall (i :: Natural) point (d :: Natural) r.
(Point_ point d r, i <= d, Has_ Vector_ i r) =>
point -> Point i r
projectPoint point
q :: Point (d-1) r) hyperPlane
h

--------------------------------------------------------------------------------
-- * Helpers for (NonVertical) Planes in R^3.

-- | Shorthand for Non-vertical hyperplanes in R^3
type Plane_ plane = NonVerticalHyperPlane_ plane 3

-- | Destructs a Plane in R^3 into the equation z = ax + by + c
pattern Plane_       :: Plane_ plane r => r -> r -> r -> plane
pattern $mPlane_ :: forall {r} {plane} {r}.
Plane_ plane r =>
plane -> (r -> r -> r -> r) -> ((# #) -> r) -> r
Plane_ a b c <- (view hyperPlaneCoefficients -> (Vector3 a b c))
{-# COMPLETE Plane_ #-}

--------------------------------------------------------------------------------
-- * Functions on Hyperplanes

-- | Test if two hyperplanes are parallel.
isParallelTo       :: ( HyperPlane_ hyperPlane  d r
                      , HyperPlane_ hyperPlane' d r
                      , Has_ Metric_ d r
                      , Num r, Eq r, 1 <= d
                      ) => hyperPlane -> hyperPlane' -> Bool
isParallelTo :: forall hyperPlane (d :: Natural) r hyperPlane'.
(HyperPlane_ hyperPlane d r, HyperPlane_ hyperPlane' d r,
 Has_ Metric_ d r, Num r, Eq r, 1 <= d) =>
hyperPlane -> hyperPlane' -> Bool
isParallelTo hyperPlane
h1 hyperPlane'
h2 = Vector d r -> Vector d r -> Bool
forall r vector (d :: Natural).
(Eq r, Num r, Metric_ vector d r) =>
vector -> vector -> Bool
isScalarMultipleOf (hyperPlane -> Vector d r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r, Eq r, 1 <= d) =>
hyperPlane -> Vector d r
normalVector hyperPlane
h1) (hyperPlane' -> Vector d r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r, Eq r, 1 <= d) =>
hyperPlane -> Vector d r
normalVector hyperPlane'
h2)

-- | Class that tells us we can construct a hyperplane from a vector
-- of points.
class HyperPlaneFromPoints hyperPlane where
  -- | Construct a hyperplane through the given d points.
  hyperPlaneThrough     :: ( Point_ point d r
                           , HyperPlane_ hyperPlane d r
                           , Num r
                           ) => Vector d point -> hyperPlane


--------------------------------------------------------------------------------
-- * Functions on Non-vertical Hyperplanes

-- | Renders the defining equation of a plane in R^3 in some human readable format.
--
-- >>> showPlaneEquation (Plane 5 3 2)
-- "z = 5x + 3y + 2"
showPlaneEquation   :: (Plane_ plane r, Show r) => plane -> String
showPlaneEquation :: forall plane r. (Plane_ plane r, Show r) => plane -> [Char]
showPlaneEquation plane
h = let (Vector3 [Char]
a [Char]
b [Char]
c) = r -> [Char]
forall a. Show a => a -> [Char]
show (r -> [Char]) -> Vector 3 r -> Vector 3 [Char]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> plane
hplane -> Getting (Vector 3 r) plane (Vector 3 r) -> Vector 3 r
forall s a. s -> Getting a s a -> a
^.Getting (Vector 3 r) plane (Vector 3 r)
forall hyperPlane (d :: Natural) r.
NonVerticalHyperPlane_ hyperPlane d r =>
Lens' hyperPlane (Vector d r)
Lens' plane (Vector 3 r)
hyperPlaneCoefficients
                      in [Char]
"z = " [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> [Char]
a [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> [Char]
"x + " [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> [Char]
b [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> [Char]
"y + " [Char] -> [Char] -> [Char]
forall a. Semigroup a => a -> a -> a
<> [Char]
c

--------------------------------------------------------------------------------
-- * Some Instances

instance (HyperPlane_ hyperPlane d r)
         => HyperPlane_ (hyperPlane :+ extra) d r where
  evalHyperPlaneEquation :: forall point.
(Num r, Point_ point d r) =>
(hyperPlane :+ extra) -> point -> r
evalHyperPlaneEquation hyperPlane :+ extra
h = hyperPlane -> point -> r
forall point. (Num r, Point_ point d r) => hyperPlane -> point -> r
forall hyperPlane (d :: Natural) r point.
(HyperPlane_ hyperPlane d r, Num r, Point_ point d r) =>
hyperPlane -> point -> r
evalHyperPlaneEquation (hyperPlane :+ extra
h(hyperPlane :+ extra)
-> Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> hyperPlane
forall s a. s -> Getting a s a -> a
^.Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core)
  {-# INLINE evalHyperPlaneEquation #-}
  hyperPlaneEquation :: Num r => (hyperPlane :+ extra) -> Vector (d + 1) r
hyperPlaneEquation = hyperPlane -> Vector (d + 1) r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r) =>
hyperPlane -> Vector (d + 1) r
hyperPlaneEquation (hyperPlane -> Vector (d + 1) r)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> Vector (d + 1) r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE hyperPlaneEquation #-}
  normalVector :: (Num r, Eq r, 1 <= d) => (hyperPlane :+ extra) -> Vector d r
normalVector = hyperPlane -> Vector d r
forall hyperPlane (d :: Natural) r.
(HyperPlane_ hyperPlane d r, Num r, Eq r, 1 <= d) =>
hyperPlane -> Vector d r
normalVector (hyperPlane -> Vector d r)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> Vector d r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE normalVector #-}
  onHyperPlane :: forall point.
(Point_ point d r, Eq r, Num r) =>
point -> (hyperPlane :+ extra) -> Bool
onHyperPlane point
q = point -> hyperPlane -> Bool
forall point.
(Point_ point d r, Eq r, Num r) =>
point -> hyperPlane -> Bool
forall hyperPlane (d :: Natural) r point.
(HyperPlane_ hyperPlane d r, Point_ point d r, Eq r, Num r) =>
point -> hyperPlane -> Bool
onHyperPlane point
q (hyperPlane -> Bool)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE onHyperPlane #-}
  onSideTest :: forall point.
(Point_ point d r, Ord r, Num r) =>
point -> (hyperPlane :+ extra) -> Ordering
onSideTest point
q = point -> hyperPlane -> Ordering
forall point.
(Point_ point d r, Ord r, Num r) =>
point -> hyperPlane -> Ordering
forall hyperPlane (d :: Natural) r point.
(HyperPlane_ hyperPlane d r, Point_ point d r, Ord r, Num r) =>
point -> hyperPlane -> Ordering
onSideTest point
q (hyperPlane -> Ordering)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> Ordering
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE onSideTest #-}

instance (ConstructableHyperPlane_ hyperPlane d r, Default extra)
          => ConstructableHyperPlane_ (hyperPlane :+ extra) d r where
  type HyperPlaneFromEquationConstraint (hyperPlane :+ extra) d r =
         HyperPlaneFromEquationConstraint hyperPlane d r
  hyperPlaneFromEquation :: HyperPlaneFromEquationConstraint (hyperPlane :+ extra) d r =>
Vector (d + 1) r -> hyperPlane :+ extra
hyperPlaneFromEquation Vector (d + 1) r
v = Vector (d + 1) r -> hyperPlane
forall hyperPlane (d :: Natural) r.
(ConstructableHyperPlane_ hyperPlane d r,
 HyperPlaneFromEquationConstraint hyperPlane d r) =>
Vector (d + 1) r -> hyperPlane
hyperPlaneFromEquation Vector (d + 1) r
v hyperPlane -> extra -> hyperPlane :+ extra
forall core extra. core -> extra -> core :+ extra
:+ extra
forall a. Default a => a
def
  {-# INLINE hyperPlaneFromEquation #-}
  fromPointAndNormal :: forall point.
(Point_ point d r, Num r) =>
point -> Vector d r -> hyperPlane :+ extra
fromPointAndNormal point
p Vector d r
v = point -> Vector d r -> hyperPlane
forall point.
(Point_ point d r, Num r) =>
point -> Vector d r -> hyperPlane
forall hyperPlane (d :: Natural) r point.
(ConstructableHyperPlane_ hyperPlane d r, Point_ point d r,
 Num r) =>
point -> Vector d r -> hyperPlane
fromPointAndNormal point
p Vector d r
v hyperPlane -> extra -> hyperPlane :+ extra
forall core extra. core -> extra -> core :+ extra
:+ extra
forall a. Default a => a
def
  {-# INLINE fromPointAndNormal #-}

instance (NonVerticalHyperPlane_  hyperPlane d r)
         => NonVerticalHyperPlane_ (hyperPlane :+ extra) d r where
  evalAt :: forall point.
(Num r, 1 <= d, Point_ point (d - 1) r) =>
point -> (hyperPlane :+ extra) -> r
evalAt point
p = point -> hyperPlane -> r
forall point.
(Num r, 1 <= d, Point_ point (d - 1) r) =>
point -> hyperPlane -> r
forall hyperPlane (d :: Natural) r point.
(NonVerticalHyperPlane_ hyperPlane d r, Num r, 1 <= d,
 Point_ point (d - 1) r) =>
point -> hyperPlane -> r
evalAt point
p (hyperPlane -> r)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE evalAt #-}
  hyperPlaneCoefficients :: Lens' (hyperPlane :+ extra) (Vector d r)
hyperPlaneCoefficients = (hyperPlane -> f hyperPlane)
-> (hyperPlane :+ extra) -> f (hyperPlane :+ extra)
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core((hyperPlane -> f hyperPlane)
 -> (hyperPlane :+ extra) -> f (hyperPlane :+ extra))
-> ((Vector d r -> f (Vector d r)) -> hyperPlane -> f hyperPlane)
-> (Vector d r -> f (Vector d r))
-> (hyperPlane :+ extra)
-> f (hyperPlane :+ extra)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Vector d r -> f (Vector d r)) -> hyperPlane -> f hyperPlane
forall hyperPlane (d :: Natural) r.
NonVerticalHyperPlane_ hyperPlane d r =>
Lens' hyperPlane (Vector d r)
Lens' hyperPlane (Vector d r)
hyperPlaneCoefficients
  {-# INLINE hyperPlaneCoefficients #-}
  verticalSideTest :: forall point.
(1 <= d, Point_ point d r, Ord r, Num r) =>
point -> (hyperPlane :+ extra) -> Ordering
verticalSideTest point
q = point -> hyperPlane -> Ordering
forall point.
(1 <= d, Point_ point d r, Ord r, Num r) =>
point -> hyperPlane -> Ordering
forall hyperPlane (d :: Natural) r point.
(NonVerticalHyperPlane_ hyperPlane d r, 1 <= d, Point_ point d r,
 Ord r, Num r) =>
point -> hyperPlane -> Ordering
verticalSideTest point
q (hyperPlane -> Ordering)
-> ((hyperPlane :+ extra) -> hyperPlane)
-> (hyperPlane :+ extra)
-> Ordering
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Getting hyperPlane (hyperPlane :+ extra) hyperPlane
-> (hyperPlane :+ extra) -> hyperPlane
forall s (m :: * -> *) a. MonadReader s m => Getting a s a -> m a
view Getting hyperPlane (hyperPlane :+ extra) hyperPlane
forall core extra core' (f :: * -> *).
Functor f =>
(core -> f core') -> (core :+ extra) -> f (core' :+ extra)
core
  {-# INLINE verticalSideTest #-}


-- | Access the last element of a vector
last :: forall vector d r. (Vector_ vector d r, 1 <= d, KnownNat (d-1)
                           ) => IndexedLens' Int vector r
last :: forall vector (d :: Natural) r.
(Vector_ vector d r, 1 <= d, KnownNat (d - 1)) =>
IndexedLens' Int vector r
last = forall (i :: Natural) vector r (d :: Natural).
(i <= (d - 1), KnownNat i, Vector_ vector d r) =>
IndexedLens' Int vector r
component @(d-1)