diff --git a/Data/LinearMap/HerMetric.hs b/Data/LinearMap/HerMetric.hs
--- a/Data/LinearMap/HerMetric.hs
+++ b/Data/LinearMap/HerMetric.hs
@@ -8,6 +8,7 @@
 {-# LANGUAGE StandaloneDeriving         #-}
 {-# LANGUAGE ConstraintKinds            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE UnicodeSyntax              #-}
 
 
 
@@ -24,6 +25,10 @@
   , spanHilbertSubspace
   , spanSubHilbertSpace
   , IsFreeSpace
+  -- * One-dimensional axes and product spaces
+  , factoriseMetric, factoriseMetric'
+  , productMetric, productMetric'
+  , metricAsLength, metric'AsLength
   -- * Utility for metrics
   , transformMetric, transformMetric'
   , dualiseMetric, dualiseMetric'
@@ -40,6 +45,8 @@
   -- * Fundamental requirements
   , MetricScalar
   , FiniteDimensional(..)
+  -- * Misc
+  , Stiefel1(..)
   ) where
     
 
@@ -366,6 +373,26 @@
   (v,w)<.>^(v',w') = v<.>^v' + w<.>^w'
   functional f = (functional $ f . (,zeroV), functional $ f . (zeroV,))
   doubleDual = id; doubleDual'= id
+instance (SmoothScalar s, Ord s, KnownNat n) => HasMetric' (s^n) where
+  type DualSpace (s^n) = s^n
+  (<.>^) = (<.>)
+  functional = fnal
+   where fnal :: ∀ s n . (SmoothScalar s, KnownNat n) => (s^n -> s) -> s^n
+         fnal f =     FreeVect . Arr.generate n $
+            \i -> f . FreeVect . Arr.generate n $ \j -> if i==j then 1 else 0
+          where Tagged n = theNatN :: Tagged n Int
+  doubleDual = id; doubleDual'= id
+instance (HasMetric v, s~Scalar v) => HasMetric' (FinVecArrRep t v s) where
+  type DualSpace (FinVecArrRep t v s) = FinVecArrRep t (DualSpace v) s
+  FinVecArrRep v <.>^ FinVecArrRep w = HMat.dot v w
+  functional = fnal
+   where fnal :: ∀ v . HasMetric v =>
+                 (FinVecArrRep t v (Scalar v) -> Scalar v)
+                       -> FinVecArrRep t (DualSpace v) (Scalar v)
+         fnal f = FinVecArrRep . (n HMat.|>)
+                     $ (f . FinVecArrRep) <$> HMat.toRows (HMat.ident n)
+         Tagged n = dimension :: Tagged v Int
+  doubleDual = id; doubleDual'= id
 
 
 
@@ -448,7 +475,51 @@
                      
 
 
-spanHilbertSubspace :: forall s v w
+-- | Project a metric on each of the factors of a product space. This works by
+--   projecting the eigenvectors into both subspaces.
+factoriseMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
+               => HerMetric (v,w) -> (HerMetric v, HerMetric w)
+factoriseMetric (HerMetric Nothing) = (HerMetric Nothing, HerMetric Nothing)
+factoriseMetric met = (sumV *** sumV) . unzip
+                   $ (projector.fst &&& projector.snd) <$> eigenSpan' met
+
+factoriseMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
+               => HerMetric' (v,w) -> (HerMetric' v, HerMetric' w)
+factoriseMetric' met = (sumV *** sumV) . unzip
+                   $ (projector'.fst &&& projector'.snd) <$> eigenSpan met
+
+productMetric :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
+               => HerMetric v -> HerMetric w -> HerMetric (v,w)
+productMetric (HerMetric Nothing) (HerMetric Nothing) = HerMetric Nothing
+productMetric (HerMetric (Just mv)) (HerMetric (Just mw))
+        = HerMetric . Just $ HMat.diagBlock [mv, mw]
+productMetric (HerMetric Nothing) (HerMetric (Just mw))
+        = HerMetric . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]
+ where (Tagged dv) = dimension :: Tagged v Int
+productMetric (HerMetric (Just mv)) (HerMetric Nothing)
+        = HerMetric . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]
+ where (Tagged dw) = dimension :: Tagged w Int
+
+productMetric' :: ∀ v w . (HasMetric v, HasMetric w, Scalar v ~ ℝ, Scalar w ~ ℝ)
+               => HerMetric' v -> HerMetric' w -> HerMetric' (v,w)
+productMetric' (HerMetric' Nothing) (HerMetric' Nothing) = HerMetric' Nothing
+productMetric' (HerMetric' (Just mv)) (HerMetric' (Just mw))
+        = HerMetric' . Just $ HMat.diagBlock [mv, mw]
+productMetric' (HerMetric' Nothing) (HerMetric' (Just mw))
+        = HerMetric' . Just $ HMat.diagBlock [HMat.konst 0 (dv,dv), mw]
+ where (Tagged dv) = dimension :: Tagged v Int
+productMetric' (HerMetric' (Just mv)) (HerMetric' Nothing)
+        = HerMetric' . Just $ HMat.diagBlock [mv, HMat.konst 0 (dw,dw)]
+ where (Tagged dw) = dimension :: Tagged w Int
+
+metricAsLength :: HerMetric ℝ -> ℝ
+metricAsLength = recip . (`metric`1)
+
+metric'AsLength :: HerMetric' ℝ -> ℝ
+metric'AsLength = recip . (`metric'`1)
+
+
+spanHilbertSubspace :: ∀ s v w
    . (HasMetric v, Scalar v ~ s, IsFreeSpace w, Scalar w ~ s)
       => HerMetric v   -- ^ Metric to induce the inner product on the Hilbert space.
           -> [v]       -- ^ @n@ linearly independent vectors, to span the subspace @w@.
@@ -477,3 +548,12 @@
           -> Option (Embedding (Linear s) w v)
 spanSubHilbertSpace = spanHilbertSubspace euclideanMetric'
 
+
+-- | The /n/-th Stiefel manifold is the space of all possible configurations of
+--   /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised
+--   vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking
+--   requires the containing space to be at least metric (if not Hilbert); we would
+--   however like to be able to use this concept also in spaces with no inner product,
+--   therefore we define this space not as normalised vectors, but rather as all
+--   vectors modulo scaling by positive factors.
+newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v }
diff --git a/Data/Manifold/Cone.hs b/Data/Manifold/Cone.hs
new file mode 100644
--- /dev/null
+++ b/Data/Manifold/Cone.hs
@@ -0,0 +1,332 @@
+-- |
+-- Module      : Data.Manifold.Cone
+-- Copyright   : (c) Justus Sagemüller 2015
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+
+{-# LANGUAGE FlexibleInstances        #-}
+{-# LANGUAGE UndecidableInstances     #-}
+{-# LANGUAGE TypeFamilies             #-}
+{-# LANGUAGE FunctionalDependencies   #-}
+{-# LANGUAGE FlexibleContexts         #-}
+{-# LANGUAGE LiberalTypeSynonyms      #-}
+{-# LANGUAGE GADTs                    #-}
+{-# LANGUAGE RankNTypes               #-}
+{-# LANGUAGE TupleSections            #-}
+{-# LANGUAGE ConstraintKinds          #-}
+{-# LANGUAGE PatternGuards            #-}
+{-# LANGUAGE TypeOperators            #-}
+{-# LANGUAGE UnicodeSyntax            #-}
+{-# LANGUAGE MultiWayIf               #-}
+{-# LANGUAGE ScopedTypeVariables      #-}
+{-# LANGUAGE RecordWildCards          #-}
+{-# LANGUAGE CPP                      #-}
+
+
+module Data.Manifold.Cone where
+    
+
+
+import Data.List
+import qualified Data.Vector.Generic as Arr
+import qualified Data.Vector
+import Data.Maybe
+import Data.Semigroup
+import Data.Function (on)
+import Data.Fixed
+
+import Data.VectorSpace
+import Data.LinearMap
+import Data.LinearMap.HerMetric
+import Data.MemoTrie (HasTrie(..))
+import Data.AffineSpace
+import Data.Basis
+import Data.Complex hiding (magnitude)
+import Data.Void
+import Data.Tagged
+import Data.Manifold.Types.Primitive
+
+import Data.CoNat
+import Data.VectorSpace.FiniteDimensional
+
+import qualified Numeric.LinearAlgebra.HMatrix as HMat
+
+import qualified Prelude
+import qualified Control.Applicative as Hask
+
+import Control.Category.Constrained.Prelude hiding ((^))
+import Control.Arrow.Constrained
+import Control.Monad.Constrained
+import Data.Foldable.Constrained
+
+import Data.Manifold.PseudoAffine
+import Data.Embedding
+
+
+
+type ConeVecArr m = FinVecArrRep Cℝay (CℝayInterior m) (Scalar (Needle m))
+type ConeNeedle m = Needle (ConeVecArr m)
+type SConn'dConeVecArr m = FinVecArrRep Cℝay (ℝ, Interior m) ℝ
+
+
+class ( Semimanifold m, Semimanifold (Interior (Interior m))
+      , Semimanifold (ConeVecArr m)
+      , Interior (ConeVecArr m) ~ ConeVecArr m )
+           => ConeSemimfd m where
+  {-# MINIMAL (fromCℝayInterior | fromCD¹Interior)
+            , (toCℝayInterior | toCD¹Interior) #-}
+  type CℝayInterior m :: *
+  
+  fromCℝayInterior :: ConeVecArr m -> Cℝay m
+  fromCℝayInterior = projCD¹ToCℝay . fromCD¹Interior
+  fromCD¹Interior :: ConeVecArr m -> CD¹ m
+  fromCD¹Interior = embCℝayToCD¹ . fromCℝayInterior
+  
+  toCℝayInterior :: Cℝay m -> Option (ConeVecArr m)
+  toCℝayInterior = toCD¹Interior . embCℝayToCD¹
+  toCD¹Interior :: CD¹ m -> Option (ConeVecArr m)
+  toCD¹Interior = toCℝayInterior . projCD¹ToCℝay
+
+  
+
+
+
+instance (ConeSemimfd m) => Semimanifold (Cℝay m) where
+  type Needle (Cℝay m) = ConeNeedle m
+  type Interior (Cℝay m) = ConeVecArr m
+  fromInterior = fromCℝayInterior
+  toInterior = toCℝayInterior
+  translateP = ctp
+   where ctp :: Tagged (Cℝay m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)
+         ctp = Tagged ctp'
+          where Tagged ctp' = translateP
+                  :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)
+  
+instance (ConeSemimfd m) => Semimanifold (CD¹ m) where
+  type Needle (CD¹ m) = ConeNeedle m
+  type Interior (CD¹ m) = ConeVecArr m
+  fromInterior = fromCD¹Interior
+  toInterior = toCD¹Interior
+  translateP = ctp
+   where ctp :: Tagged (CD¹ m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)
+         ctp = Tagged ctp'
+          where Tagged ctp' = translateP
+                  :: Tagged (ConeVecArr m) (ConeVecArr m -> ConeNeedle m -> ConeVecArr m)
+
+instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (Cℝay m) where
+  p.-~.i = (.-~.i) =<< toInterior p
+instance (ConeSemimfd m, SmoothScalar (Scalar (Needle m))) => PseudoAffine (CD¹ m) where
+  p.-~.i = (.-~.i) =<< toInterior p
+
+
+instance ConeSemimfd (ZeroDim ℝ) where
+  type CℝayInterior (ZeroDim ℝ) = ℝ
+  fromCℝayInterior (FinVecArrRep qb) | HMat.size qb == 0  = Cℝay 1 Origin
+                                     | x <- qb HMat.! 0   = Cℝay (bijectℝtoℝplus x) Origin 
+  toCℝayInterior (Cℝay 0 Origin) = Hask.empty
+  toCℝayInterior (Cℝay y Origin) = pure . FinVecArrRep $ 1 HMat.|>[bijectℝplustoℝ y]
+instance ConeSemimfd ℝ where
+  type CℝayInterior ℝ = ℝ²
+  fromCℝayInterior (FinVecArrRep qb) = Cℝay (q'+b') (q'-b')
+   where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb
+  toCℝayInterior (Cℝay 0 _) = Hask.empty
+  toCℝayInterior (Cℝay h x) = pure . FinVecArrRep 
+                              . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h+x, h-x]
+  fromCD¹Interior (FinVecArrRep qb) = CD¹ (bijectℝplustoIntv $ q'+b') (q'-b')
+   where [q', b'] = HMat.toList $ HMat.cmap ((/2) . bijectℝtoℝplus) qb
+  toCD¹Interior (CD¹ h x) = pure . FinVecArrRep
+                              . HMat.cmap bijectℝplustoℝ $ HMat.fromList [h'+x, h'-x]
+   where h' = bijectIntvtoℝplus h
+
+instance ConeSemimfd S⁰ where
+  type CℝayInterior S⁰ = ℝ
+  fromCℝayInterior xa | x>0        = Cℝay x PositiveHalfSphere
+                      | otherwise  = Cℝay (-x) NegativeHalfSphere
+   where x = getFinVecArrRep xa HMat.! 0
+  toCℝayInterior (Cℝay x PositiveHalfSphere) = return . FinVecArrRep $ HMat.scalar x
+  toCℝayInterior (Cℝay x NegativeHalfSphere) = return . FinVecArrRep . HMat.scalar $ -x
+  fromCD¹Interior xa | x>0        = CD¹ (bijectℝtoIntv x) PositiveHalfSphere
+                     | otherwise  = CD¹ (-bijectℝtoIntv x) NegativeHalfSphere
+   where x = getFinVecArrRep xa HMat.! 0
+  toCD¹Interior (CD¹ 1 _) = Hask.empty
+  toCD¹Interior (CD¹ x PositiveHalfSphere)
+        = return . FinVecArrRep . HMat.scalar $ bijectIntvtoℝ x
+  toCD¹Interior (CD¹ x NegativeHalfSphere)
+        = return . FinVecArrRep . HMat.scalar $ -bijectℝtoIntv x
+
+
+instance ConeSemimfd S¹ where
+  type CℝayInterior S¹ = ℝ²
+  fromCℝayInterior (FinVecArrRep xy) = Cℝay r (S¹ $ atan2 y x)
+   where r = HMat.norm_2 xy
+         [x,y] = HMat.toList xy
+  toCℝayInterior (Cℝay r (S¹ φ)) = return . FinVecArrRep
+                    . HMat.scale r $ HMat.fromList [cos φ, sin φ]
+  fromCD¹Interior (FinVecArrRep xy) = CD¹ (bijectℝtoIntv r) (S¹ $ atan2 y x)
+   where r = HMat.norm_2 xy
+         [x,y] = HMat.toList xy
+  toCD¹Interior (CD¹ 1 _) = Hask.empty
+  toCD¹Interior (CD¹ r (S¹ φ)) = return . FinVecArrRep
+                    . HMat.scale r' $ HMat.fromList [cos φ, sin φ]
+   where r' = bijectIntvtoℝ r
+
+
+instance ConeSemimfd S² where
+  type CℝayInterior S² = ℝ³
+  fromCℝayInterior (FinVecArrRep xyz) = Cℝay r (S² (acos $ z/r) (atan2 y x))
+   where r = HMat.norm_2 xyz
+         [x,y,z] = HMat.toList xyz
+  toCℝayInterior (Cℝay r (S² ϑ φ)) = return . FinVecArrRep
+                    . HMat.scale r $ HMat.fromList [w*x₀, w*y₀, z₀]
+   where x₀ = cos φ; y₀ = sin φ; z₀ = cos ϑ; w = sin ϑ
+
+                                      
+
+
+-- | Products of simply connected spaces.
+instance ( PseudoAffine x, PseudoAffine y
+         , WithField ℝ HilbertSpace (Interior x), WithField ℝ HilbertSpace (Interior y)
+         , LinearManifold (FinVecArrRep Cℝay (ℝ, (Interior x, Interior y)) ℝ)
+         ) => ConeSemimfd (x,y) where
+  type CℝayInterior (x,y) = (ℝ, (Interior x, Interior y))
+  fromCℝayInterior = simplyCncted_fromCℝayInterior
+  toCℝayInterior = simplyCncted_toCℝayInterior
+
+instance ( KnownNat n ) => ConeSemimfd (ℝ^n) where
+  type CℝayInterior (ℝ^n) = (ℝ, ℝ^n)
+  fromCℝayInterior = simplyCncted_fromCℝayInterior
+  toCℝayInterior = simplyCncted_toCℝayInterior
+
+instance ( HilbertSpace (FinVecArrRep t v ℝ) ) => ConeSemimfd (FinVecArrRep t v ℝ) where
+  type CℝayInterior (FinVecArrRep t v ℝ) = (ℝ, FinVecArrRep t v ℝ)
+  fromCℝayInterior = simplyCncted_fromCℝayInterior
+  toCℝayInterior = simplyCncted_toCℝayInterior
+
+
+  
+instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x)
+         , HilbertSpace (CℝayInterior x)
+         , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ)
+         ) => ConeSemimfd (CD¹ x) where
+  type CℝayInterior (CD¹ x) = (ℝ, ConeVecArr x)
+  fromCℝayInterior i = Cℝay h (embCℝayToCD¹ o)
+   where (Cℝay h o) = simplyCncted_fromCℝayInterior i
+  toCℝayInterior (Cℝay _ (CD¹ 1 _)) = Hask.empty
+  toCℝayInterior (Cℝay h p) = simplyCncted_toCℝayInterior $ Cℝay h (projCD¹ToCℝay p)
+  
+  
+instance ( WithField ℝ ConeSemimfd x, PseudoAffine (Cℝay x)
+         , HilbertSpace (CℝayInterior x)
+         , HilbertSpace (FinVecArrRep Cℝay (CℝayInterior x) ℝ)
+         ) => ConeSemimfd (Cℝay x) where
+  type CℝayInterior (Cℝay x) = (ℝ, ConeVecArr x)
+  fromCℝayInterior = simplyCncted_fromCℝayInterior
+  toCℝayInterior = simplyCncted_toCℝayInterior
+  
+  
+simplyCncted_fromCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x))
+        => SConn'dConeVecArr x -> Cℝay x
+simplyCncted_fromCℝayInterior (FinVecArrRep ri) = Cℝay h . fromInterior . fromPackedVector
+                         $ subtract (h/n) `Arr.map` Arr.tail cmps
+   where h = Arr.sum cmps
+         cmps = bijectℝtoℝplus `HMat.cmap` ri
+         n = fromIntegral $ Arr.length cmps
+  
+simplyCncted_toCℝayInterior :: (PseudoAffine x, WithField ℝ HilbertSpace (Interior x))
+        => Cℝay x -> Option (SConn'dConeVecArr x)
+simplyCncted_toCℝayInterior (Cℝay h v) | h/=0, Option (Just vi) <- toInterior v 
+   = let cmps'' = asPackedVector vi
+         cmps' = (+ h/n) `HMat.cmap` cmps''
+         cmps = (h - Arr.sum cmps') `Arr.cons` cmps
+         n = fromIntegral $ Arr.length cmps
+     in return $ FinVecArrRep (bijectℝplustoℝ `Arr.map` cmps)
+simplyCncted_toCℝayInterior (Cℝay _ _) = Hask.empty
+
+
+-- Some essential homeomorphisms
+bijectℝtoℝplus      , bijectℝplustoℝ
+ , bijectIntvtoℝplus, bijectℝplustoIntv
+ ,     bijectIntvtoℝ, bijectℝtoIntv
+               :: ℝ -> ℝ
+
+bijectℝplustoℝ x = x - 1/x
+bijectℝtoℝplus y = y/2 + sqrt(y^2/4 + 1)
+
+-- [0, 1[ ↔ ℝ⁺
+bijectℝplustoIntv y = 1 - recip (y+1)
+bijectIntvtoℝplus x = recip(1-x) - 1
+
+-- ]-1, 1[ ↔ ℝ  (Similar to 'tanh', but converges less quickly towards ±1.)
+bijectℝtoIntv y | y>0        = -1/(2*y) + sqrt(1/(4*y^2) + 1)
+                | y<0        = -1/(2*y) - sqrt(1/(4*y^2) + 1)
+                | otherwise  = 0
+                 -- 0 = x² + x/y - 1
+                 -- x = -1/2y ± sqrt(1/4y² + 1)
+bijectIntvtoℝ x = x / (1-x^2)
+
+embCℝayToCD¹ :: Cℝay m -> CD¹ m
+embCℝayToCD¹ (Cℝay h m) = CD¹ (bijectℝplustoIntv h) m
+
+projCD¹ToCℝay :: CD¹ m -> Cℝay m
+projCD¹ToCℝay (CD¹ h m) = Cℝay (bijectIntvtoℝplus h) m
+
+-- instance (WithScalar ℝ PseudoAffine m) => Semimanifold (Cℝay m) where
+--   type Needle (Cℝay m) = (Needle m, ℝ)
+--   type Interior (Cℝay m) = (Interior m, ℝ)
+-- 
+--   fromInterior (im, d)
+--      | d>38       = Cℝay m d  -- from 38 on, the +1 is numerically
+--                               -- insignificant against the exponential.
+--      | otherwise  = cℝay m (log $ exp d + 1)
+--                -- note that (for the same reason we can shortcut above 38)
+--                -- such negative arguments will actually yield the value zero.
+--                -- This means we're actually reaching the “infinitely far”
+--                -- rim rather quickly. This might be a problem, but normally
+--                -- shouldn't really matter much.
+--                -- It would perhaps be better to have homeomorphism that
+--                -- approaches -1/x in the negative limit, but such a
+--                -- function doesn't seem as easy to come by.
+--    where m = fromInterior im
+--   toInterior (Cℝay m q)
+--      | q>38       = fmap (,q) im
+--      | q>0        = fmap (, log $ exp d - 1) im
+--      | otherwise  = Hask.empty
+--    where im = toInterior m
+
+stiefel1Project :: LinearManifold v =>
+             DualSpace v       -- ^ Must be nonzero.
+                 -> Stiefel1 v
+stiefel1Project = Stiefel1
+
+stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v
+stiefel1Embed (Stiefel1 n) = normalized n
+  
+
+class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualSpace v))
+          => HasUnitSphere v where
+  type UnitSphere v :: *
+  stiefel :: UnitSphere v -> Stiefel1 v
+  stiefel = Stiefel1 . embed
+  unstiefel :: Stiefel1 v -> UnitSphere v
+  unstiefel = coEmbed . getStiefel1N
+
+instance HasUnitSphere ℝ  where type UnitSphere ℝ  = S⁰
+instance HasUnitSphere (FinVecArrRep t ℝ ℝ) where type UnitSphere (FinVecArrRep t ℝ ℝ)   = S⁰
+
+instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹
+instance HasUnitSphere (FinVecArrRep t ℝ² ℝ) where type UnitSphere (FinVecArrRep t ℝ² ℝ) = S¹
+
+instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S²
+instance HasUnitSphere (FinVecArrRep t ℝ³ ℝ) where type UnitSphere (FinVecArrRep t ℝ³ ℝ) = S²
+
+-- instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where
+--   embed = embed . unstiefel
+--   coEmbed = stiefel . coEmbed
+
+
+
+
diff --git a/Data/Manifold/Griddable.hs b/Data/Manifold/Griddable.hs
new file mode 100644
--- /dev/null
+++ b/Data/Manifold/Griddable.hs
@@ -0,0 +1,182 @@
+-- |
+-- Module      : Data.Manifold.Griddable
+-- Copyright   : (c) Justus Sagemüller 2015
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE DeriveFunctor              #-}
+{-# LANGUAGE DeriveFoldable             #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE FunctionalDependencies     #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE GADTs                      #-}
+{-# LANGUAGE RankNTypes                 #-}
+{-# LANGUAGE TupleSections              #-}
+{-# LANGUAGE ParallelListComp           #-}
+{-# LANGUAGE UnicodeSyntax              #-}
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE PatternGuards              #-}
+{-# LANGUAGE LambdaCase                 #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE LiberalTypeSynonyms        #-}
+{-# LANGUAGE RecordWildCards            #-}
+{-# LANGUAGE DataKinds                  #-}
+
+
+module Data.Manifold.Griddable (GridAxis(..), Griddable(..)) where
+
+
+import Data.List hiding (filter, all, elem, sum)
+import Data.Maybe
+import qualified Data.Map as Map
+import qualified Data.Vector as Arr
+import Data.List.NonEmpty (NonEmpty(..))
+import Data.List.FastNub
+import qualified Data.List.NonEmpty as NE
+import Data.Semigroup
+
+import Data.VectorSpace
+import Data.LinearMap
+import Data.LinearMap.HerMetric
+import Data.LinearMap.Category
+import Data.AffineSpace
+import Data.Basis
+import Data.Complex hiding (magnitude)
+import Data.Void
+import Data.Tagged
+import Data.Proxy
+
+import Data.SimplicialComplex
+import Data.Manifold.Types
+import Data.Manifold.Types.Primitive ((^), (^.))
+import Data.Manifold.PseudoAffine
+import Data.Manifold.TreeCover (Shade(..), fullShade, shadeCtr, shadeExpanse)
+    
+import Data.Embedding
+import Data.CoNat
+
+import qualified Prelude as Hask hiding(foldl, sum, sequence)
+import qualified Control.Applicative as Hask
+import qualified Control.Monad       as Hask hiding(forM_, sequence)
+import Data.Functor.Identity
+import Control.Monad.Trans.State
+import Control.Monad.Trans.Writer
+import Control.Monad.Trans.Class
+import qualified Data.Foldable       as Hask
+import Data.Foldable (all, elem, toList, sum)
+import qualified Data.Traversable as Hask
+import Data.Traversable (forM)
+
+import qualified Numeric.LinearAlgebra.HMatrix as HMat
+
+import Control.Category.Constrained.Prelude hiding
+     ((^), all, elem, sum, forM, Foldable(..), Traversable)
+import Control.Arrow.Constrained
+import Control.Monad.Constrained hiding (forM)
+import Data.Foldable.Constrained
+
+import Text.Printf
+import GHC.Generics (Generic)
+
+
+data GridAxis m g = GridAxInterval (Shade m)
+                  | GridAxCons (Shade m) g (GridAxis m g)
+                  | GridAxisClosed g (GridAxis m g)
+             deriving (Hask.Functor)
+
+gshmap :: (Shade m -> Shade n) -> GridAxis m g -> GridAxis n g
+gshmap f (GridAxInterval i) = GridAxInterval $ f i
+gshmap f (GridAxCons i g ax) = GridAxCons (f i) g $ gshmap f ax
+gshmap f (GridAxisClosed g ax) = GridAxisClosed g $ gshmap f ax
+
+axisEnumFromStepTo :: (ℝ->a) -> ℝ -> ℝ -> ℝ -> GridAxis ℝ a
+axisEnumFromStepTo f l st r
+    | l' > r   = GridAxInterval $ intvl2Shade (Interval l l')
+    | otherwise  = GridAxCons (intvl2Shade $ Interval l l')
+                              (f l') $ axisEnumFromStepTo f l' st r
+ where l' = l+st
+
+axisGrLength :: GridAxis m a -> Int
+axisGrLength (GridAxInterval _) = 0
+axisGrLength (GridAxCons _ _ ax) = 1 + axisGrLength ax
+axisGrLength (GridAxisClosed _ ax) = axisGrLength ax
+
+class (WithField ℝ Manifold m) => Griddable m g where
+  data GriddingParameters m g :: *
+  mkGridding :: GriddingParameters m g -> Int -> Shade m -> [GridAxis m g]
+
+
+instance Griddable ℝ String where
+  data GriddingParameters ℝ String = ℝGridParam
+  mkGridding ℝGridParam n (Shade c expa') = [ax]
+   where l = c - expa
+         r = c + expa
+         
+         expa = metric'AsLength expa'
+         
+         (Just ax) = find ((>=n) . axisGrLength)
+                $ [ let qe = 10^^lqe' * nb
+                    in axisEnumFromStepTo (prettyFloatShow lqe')
+                         ( qe * fromIntegral (floor $ l / qe) ) qe r
+                  | lqe' <- [lqe - 1, lqe - 2 ..], nb <- [5, 2, 1] ]
+         
+         lqe = lqef expa :: Int
+         lqef n | n > 0      = floor $ lg   n
+                | n < 0      = floor $ lg (-n)
+
+
+instance (Griddable m a, Griddable n a) => Griddable (m,n) a where
+  data GriddingParameters (m,n) a = PairGriddingParameters {
+               fstGriddingParams :: GriddingParameters m a
+             , sndGriddingParams :: GriddingParameters n a }
+  mkGridding (PairGriddingParameters p₁ p₂) n (Shade (c₁,c₂) e₁e₂)
+          = gshmap ( uncurry fullShade . (                  (,c₂).(^.shadeCtr)
+                                         &&& (`productMetric'`e₂).(^.shadeExpanse)) )
+              <$> g₁s
+         ++ gshmap ( uncurry fullShade . (                  (c₁,).(^.shadeCtr)
+                                         &&& ( productMetric' e₁).(^.shadeExpanse)) )
+              <$> g₂s
+   where g₁s = mkGridding p₁ n $ fullShade c₁ e₁
+         g₂s = mkGridding p₂ n $ fullShade c₂ e₂
+         (e₁,e₂) = factoriseMetric' e₁e₂ 
+
+prettyFloatShow :: Int -> Double -> String
+prettyFloatShow _ 0 = "0"
+prettyFloatShow preci x
+    | preci >= 0, preci < 4  = show $ round x
+    | preci < 0, preci > -2  = printf "%.1f" x
+    | otherwise   = case ceiling (0.01 + lg (abs x/10^^(preci+1))) + preci of
+                        0    | preci < 0  -> printf ("%."++show(-preci)++"f") x
+                        expn | expn>preci -> printf ("%."++show(expn-preci)++"f*10^%i")
+                                                      (x/10^^expn)                 expn
+                             | otherwise  -> printf ("%i*10^%i")
+                                                      (round $ x/10^^expn :: Int)  expn
+
+
+
+data Interval = Interval { ivLBound, ivRBound :: ℝ }
+
+shade2Intvl :: Shade ℝ -> Interval
+shade2Intvl sh = Interval l r
+ where c = sh ^. shadeCtr
+       expa = metric'AsLength $ sh ^. shadeExpanse
+       l = c - expa; r = c + expa
+
+intvl2Shade :: Interval -> Shade ℝ
+intvl2Shade (Interval l r) = fullShade c (projector' expa)
+ where c = (l+r) / 2
+       expa = (r-l) / 2
+       
+
+lg :: Floating a => a -> a
+lg = logBase 10
+
diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs
--- a/Data/Manifold/PseudoAffine.hs
+++ b/Data/Manifold/PseudoAffine.hs
@@ -17,6 +17,14 @@
 -- diffeomorphic. At the moment, we mainly focus on /region-wise differentiable functions/,
 -- which are a promising compromise between flexibility of definition and provability of
 -- analytic properties. In particular, they are well-suited for visualisation purposes.
+-- 
+-- The classes in this module are mostly aimed at manifolds /without boundary/.
+-- Manifolds with boundary (which we call @MWBound@, never /manifold/!)
+-- are more or less treated as a disjoint sum of the interior and the boundary.
+-- To understand how this module works, best first forget about boundaries – in this case,
+-- @'Interior' x ~ x@, 'fromInterior' and 'toInterior' are trivial, and
+-- '.+~|', '|-~.' and 'betweenBounds' are irrelevant.
+-- The manifold structure of the boundary itself is not considered at all here.
 
 {-# LANGUAGE FlexibleInstances        #-}
 {-# LANGUAGE UndecidableInstances     #-}
@@ -84,10 +92,12 @@
 import Data.Manifold.Types.Primitive
 
 import Data.CoNat
+import Data.VectorSpace.FiniteDimensional
 
 import qualified Numeric.LinearAlgebra.HMatrix as HMat
 
 import qualified Prelude
+import qualified Control.Applicative as Hask
 
 import Control.Category.Constrained.Prelude hiding ((^))
 import Control.Arrow.Constrained
@@ -100,11 +110,13 @@
 infix 6 .-~.
 infixl 6 .+~^, .-~^
 
-class (AdditiveGroup (Needle x)) => Semimanifold x where
+class ( AdditiveGroup (Needle x), Interior (Interior x) ~ Interior x )
+          => Semimanifold x where
+  {-# MINIMAL ((.+~^) | fromInterior), toInterior, translateP #-}
   -- | The space of &#x201c;natural&#x201d; ways starting from some reference point
   --   and going to some particular target point. Hence,
   --   the name: like a compass needle, but also with an actual length.
-  --   For affine space, 'Needle' is simply the space of
+  --   For affine spaces, 'Needle' is simply the space of
   --   line segments (aka vectors) between two points, i.e. the same as 'Diff'.
   --   The 'AffineManifold' constraint makes that requirement explicit.
   -- 
@@ -112,9 +124,37 @@
   --   used somewhat synonymously).
   type Needle x :: *
   
-  -- | Generalised translation operation.
-  (.+~^) :: x -> Needle x -> x
+  -- | Manifolds with boundary are a bit tricky. We support such manifolds,
+  --   but carry out most calculations only in “the fleshy part” – the
+  --   interior, which is an “infinite space”, so you can arbitrarily scale paths.
+  -- 
+  --   The default implementation is @'Interior' x = x@, which corresponds
+  --   to a manifold that has no boundary to begin with.
+  type Interior x :: *
+  type Interior x = x
   
+  -- | Generalised translation operation. Note that the result will always also
+  --   be in the interior; scaling up the needle can only get you ever /closer/
+  --   to a boundary.
+  (.+~^) :: Interior x -> Needle x -> x
+  (.+~^) = addvp
+   where addvp :: ∀ x . Semimanifold x => Interior x -> Needle x -> x
+         addvp p = fromInterior . tp p
+          where (Tagged tp) = translateP :: Tagged x (Interior x -> Needle x -> Interior x)
+    
+  -- | 'id' sans boundary.
+  fromInterior :: Interior x -> x
+  fromInterior p = p .+~^ zeroV 
+  
+  toInterior :: x -> Option (Interior x)
+  
+  -- | The signature of '.+~^' should really be @'Interior' x -> 'Needle' x -> 'Interior' x@,
+  --   only, this is not possible because it only consists of non-injective type families.
+  --   The solution is this tagged signature, which is of course rather unwieldy. That's
+  --   why '.+~^' has the stronger, but easier usable signature. Without boundary, these
+  --   functions should be equivalent, i.e. @translateP = Tagged (.+~^)@.
+  translateP :: Tagged x (Interior x -> Needle x -> Interior x)
+  
   -- | Shorthand for @\\p v -> p .+~^ 'negateV' v@, which should obey the /asymptotic/ law
   --   
   -- @
@@ -126,17 +166,19 @@
   --   as /O/ (/&#x3b7;/&#xb2;). For large vectors, it will however behave differently,
   --   except in flat spaces (where all this should be equivalent to the 'AffineSpace'
   --   instance).
-  (.-~^) :: x -> Needle x -> x
+  (.-~^) :: Interior x -> Needle x -> x
   p .-~^ v = p .+~^ negateV v
 
--- | This is the class underlying manifolds. ('Manifold' only adds an extra constraint that
---   would be circular if it was in a single class. You can always just use 'Manifold'
---   as a constraint in your signatures, but you must /define/ only 'PseudoAffine' for
---   manifold types &#x2013; the 'Manifold' instance follows universally from this.)
+  
+-- | This is the class underlying manifolds. ('Manifold' only precludes boundaries
+--   and adds an extra constraint that would be circular if it was in a single
+--   class. You can always just use 'Manifold' as a constraint in your signatures,
+--   but you must /define/ only 'PseudoAffine' for manifold types &#x2013;
+--   the 'Manifold' instance follows universally from this, if @'Interior x ~ x@.)
 --   
---   The interface is almost identical to the better-known 'AffineSpace' class, but unlike
---   in the mathematical definition of affine spaces we don't require associativity 
---   of '.+~^' with '^+^' &#x2013; except in an asymptotic sense for small vectors.
+--   The interface is (boundaries aside) almost identical to the better-known
+--   'AffineSpace' class, but we don't require associativity of '.+~^' with '^+^'
+--   &#x2013; except in an /asymptotic sense/ for small vectors.
 --   
 --   That innocent-looking change makes the class applicable to vastly more general types:
 --   while an affine space is basically nothing but a vector space without particularly
@@ -146,22 +188,38 @@
 --   manifolds in their usual maths definition (with an atlas of charts: a family of
 --   overlapping regions of the topological space, each homeomorphic to the 'Needle'
 --   vector space or some simply-connected subset thereof).
-class Semimanifold x => PseudoAffine x where
+class ( Semimanifold x, Semimanifold (Interior x)
+      , Needle (Interior x) ~ Needle x, Interior (Interior x) ~ Interior x)
+        => PseudoAffine x where
   -- | The path reaching from one point to another.
-  --   Should only yield 'Nothing' if the points are on disjoint segments of a
-  --   non&#x2013;path-connected manifold. Otherwise, the identity
+  --   Should only yield 'Nothing' if
+  -- 
+  --   * The points are on disjoint segments of a non&#x2013;path-connected space.
+  -- 
+  --   * Either of the points is on the boundary. Use '|-~.' to deal with this.
+  -- 
+  --   On manifolds, the identity
   --   
   -- @
   -- p .+~^ (q.-~.p) &#x2261; q
   -- @
   --   
   --   should hold, at least save for floating-point precision limits etc..
-  (.-~.) :: x -> x -> Option (Needle x)
+  -- 
+  --   '.-~.' and '.+~^' only really work in manifolds without boundary. If you consider
+  --   the path between two points, one of which lies on the boundary, it can't really
+  --   be possible to scale this path any longer – it would have to reach “out of the
+  --   manifold”. To adress this problem, these functions basically consider only the
+  --   /interior/ of the space.
+  (.-~.) :: x -> Interior x -> Option (Needle x)
+
   
+  
+  
 
 -- | See 'Semimanifold' and 'PseudoAffine' for the methods.
-class (PseudoAffine m, LinearManifold (Needle m)) => Manifold m
-instance (PseudoAffine m, LinearManifold (Needle m)) => Manifold m
+class (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m
+instance (PseudoAffine m, LinearManifold (Needle m), Interior m ~ m) => Manifold m
 
 type LocallyScalable s x = ( PseudoAffine x, (Needle x) ~ Needle x
                            , HasMetric (Needle x)
@@ -173,7 +231,7 @@
 --   
 --   (Actually, 'LinearManifold' is stronger than 'VectorSpace' at the moment, since
 --   'HasMetric' requires 'FiniteDimensional'. This might be lifted in the future.)
-type LinearManifold x = ( PseudoAffine x, Needle x ~ x, HasMetric x )
+type LinearManifold x = ( PseudoAffine x, Interior x ~ x, Needle x ~ x, HasMetric x )
 
 -- | Require some constraint on a manifold, and also fix the type of the manifold's
 --   underlying field. For example, @WithField &#x211d; 'HilbertSpace' v@ constrains
@@ -185,12 +243,12 @@
 type WithField s c x = ( c x, s ~ Scalar (Needle x) )
 
 -- | The 'RealFloat' class plus manifold constraints.
-type RealDimension r = ( PseudoAffine r, Needle r ~ r
+type RealDimension r = ( PseudoAffine r, Interior r ~ r, Needle r ~ r
                        , HasMetric r, DualSpace r ~ r, Scalar r ~ r
                        , RealFloat r )
 
 -- | The 'AffineSpace' class plus manifold constraints.
-type AffineManifold m = ( PseudoAffine m, AffineSpace m
+type AffineManifold m = ( PseudoAffine m, Interior m ~ m, AffineSpace m
                         , Needle m ~ Diff m, LinearManifold (Diff m) )
 
 -- | A Hilbert space is a /complete/ inner product space. Being a vector space, it is
@@ -200,7 +258,8 @@
 --   but since 'Manifold's are at the moment confined to finite dimension, they are in
 --   fact (trivially) complete.)
 type HilbertSpace x = ( LinearManifold x, InnerSpace x
-                      , Needle x ~ x, DualSpace x ~ x, Floating (Scalar x) )
+                      , Interior x ~ x, Needle x ~ x, DualSpace x ~ x
+                      , Floating (Scalar x) )
 
 -- | An euclidean space is a real affine space whose tangent space is a Hilbert space.
 type EuclidSpace x = ( AffineManifold x, InnerSpace (Diff x)
@@ -215,46 +274,98 @@
 type Metric' x = HerMetric' (Needle x)
 
 
--- | Interpolate between points, approximately linearly.
-palerp :: (PseudoAffine x, VectorSpace (Needle x))
-    => x -> x -> Option (Scalar (Needle x) -> x)
-palerp p1 p2 = fmap (\v t -> p1 .+~^ t *^ v) $ p2 .-~. p1
+-- | Interpolate between points, approximately linearly. For
+--   points that aren't close neighbours (i.e. lie in an almost
+--   flat region), the pathway is basically undefined – save for
+--   its end points.
+-- 
+--   A proper, really well-defined (on global scales) interpolation
+--   only makes sense on a Riemannian manifold, as geodesics.
+--   This is a task to be tackled in the future.
+palerp :: ∀ x. Manifold x
+    => Interior x -> Interior x -> Option (Scalar (Needle x) -> x)
+palerp p1 p2 = case (fromInterior p2 :: x) .-~. p1 of
+  Option (Just v) -> return $ \t -> p1 .+~^ t *^ v
+  _ -> Hask.empty
 
 
 
 #define deriveAffine(t)          \
 instance Semimanifold (t) where { \
   type Needle (t) = Diff (t);      \
-  (.+~^) = (.+^) };                 \
-instance PseudoAffine (t) where {    \
+  fromInterior = id;                \
+  toInterior = pure;                 \
+  translateP = Tagged (.+^);          \
+  (.+~^) = (.+^) };                    \
+instance PseudoAffine (t) where {       \
   a.-~.b = pure (a.-.b);      }
 
 deriveAffine(Double)
 deriveAffine(Rational)
 
+instance SmoothScalar s => Semimanifold (FinVecArrRep t b s) where
+  type Needle (FinVecArrRep t b s) = FinVecArrRep t b s
+  type Interior (FinVecArrRep t b s) = FinVecArrRep t b s
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+^)
+  (.+~^) = (.+^)
+instance SmoothScalar s => PseudoAffine (FinVecArrRep t b s) where
+  a.-~.b = pure (a.-.b)
+  
+
 instance Semimanifold (ZeroDim k) where
   type Needle (ZeroDim k) = ZeroDim k
+  fromInterior = id
+  toInterior = pure
   Origin .+~^ Origin = Origin
   Origin .-~^ Origin = Origin
+  translateP = Tagged (.+~^)
 instance PseudoAffine (ZeroDim k) where
   Origin .-~. Origin = pure Origin
 
 instance (Semimanifold a, Semimanifold b) => Semimanifold (a,b) where
   type Needle (a,b) = (Needle a, Needle b)
+  type Interior (a,b) = (Interior a, Interior b)
   (a,b).+~^(v,w) = (a.+~^v, b.+~^w)
   (a,b).-~^(v,w) = (a.-~^v, b.-~^w)
+  fromInterior (i,j) = (fromInterior i, fromInterior j)
+  toInterior (a,b) = fzip (toInterior a, toInterior b)
+  translateP = tp
+   where tp :: ∀ a b . (Semimanifold a, Semimanifold b)
+                     => Tagged (a,b) ( (Interior a, Interior b) 
+                                    -> (Needle a, Needle b)
+                                    -> (Interior a, Interior b) )
+         tp = Tagged $ \(a,b) (v,w) -> (ta a v, tb b w)
+          where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a)
+                Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b)
 instance (PseudoAffine a, PseudoAffine b) => PseudoAffine (a,b) where
   (a,b).-~.(c,d) = liftA2 (,) (a.-~.c) (b.-~.d)
 
 instance (Semimanifold a, Semimanifold b, Semimanifold c) => Semimanifold (a,b,c) where
   type Needle (a,b,c) = (Needle a, Needle b, Needle c)
+  type Interior (a,b,c) = (Interior a, Interior b, Interior c)
   (a,b,c).+~^(v,w,x) = (a.+~^v, b.+~^w, c.+~^x)
   (a,b,c).-~^(v,w,x) = (a.-~^v, b.-~^w, c.-~^x)
+  fromInterior (i,j,k) = (fromInterior i, fromInterior j, fromInterior k)
+  toInterior (a,b,c) = liftA3 (,,) (toInterior a) (toInterior b) (toInterior c)
+  translateP = tp
+   where tp :: ∀ a b v . (Semimanifold a, Semimanifold b, Semimanifold c)
+                     => Tagged (a,b,c) ( (Interior a, Interior b, Interior c) 
+                                      -> (Needle a, Needle b, Needle c)
+                                      -> (Interior a, Interior b, Interior c) )
+         tp = Tagged $ \(a,b,c) (v,w,x) -> (ta a v, tb b w, tc c x)
+          where Tagged ta = translateP :: Tagged a (Interior a -> Needle a -> Interior a)
+                Tagged tb = translateP :: Tagged b (Interior b -> Needle b -> Interior b)
+                Tagged tc = translateP :: Tagged c (Interior c -> Needle c -> Interior c)
 instance (PseudoAffine a, PseudoAffine b, PseudoAffine c) => PseudoAffine (a,b,c) where
   (a,b,c).-~.(d,e,f) = liftA3 (,,) (a.-~.d) (b.-~.e) (c.-~.f)
 
 instance (MetricScalar a, KnownNat n) => Semimanifold (FreeVect n a) where
   type Needle (FreeVect n a) = FreeVect n a
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   (.+~^) = (.+^)
 instance (MetricScalar a, KnownNat n) => PseudoAffine (FreeVect n a) where
   a.-~.b = pure (a.-.b)
@@ -262,6 +373,9 @@
 
 instance Semimanifold S⁰ where
   type Needle S⁰ = ℝ⁰
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   p .+~^ Origin = p
   p .-~^ Origin = p
 instance PseudoAffine S⁰ where
@@ -271,6 +385,9 @@
 
 instance Semimanifold S¹ where
   type Needle S¹ = ℝ
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   S¹ φ₀ .+~^ δφ
      | φ' < 0     = S¹ $ φ' + tau
      | otherwise  = S¹ $ φ'
@@ -282,8 +399,25 @@
      | otherwise   = pure δφ
    where δφ = φ₁ - φ₀
 
+instance Semimanifold D¹ where
+  type Needle D¹ = ℝ
+  type Interior D¹ = ℝ
+  fromInterior = D¹ . tanh
+  toInterior (D¹ x) | abs x < 1  = return $ atanh x
+                    | otherwise  = Hask.empty
+  translateP = Tagged (+)
+instance PseudoAffine D¹ where
+  D¹ 1 .-~. _ = Hask.empty
+  D¹ (-1) .-~. _ = Hask.empty
+  D¹ x .-~. y
+    | abs x < 1  = return $ atanh x - y
+    | otherwise  = Hask.empty
+
 instance Semimanifold S² where
   type Needle S² = ℝ²
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   S² ϑ₀ φ₀ .+~^ δv
      | ϑ₀ < pi/2  = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv
      | otherwise  = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv
@@ -304,6 +438,9 @@
 
 instance Semimanifold ℝP² where
   type Needle ℝP² = ℝ²
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   ℝP² r₀ φ₀ .+~^ (δr, δφ)
    | r₀ > 1/2   = case r₀ + δr of
                    r₁ | r₁ > 1     -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi)
@@ -323,16 +460,16 @@
    | otherwise  = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) )
 
 
-instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)
-             => Semimanifold (CD¹ m) where
-  type Needle (CD¹ m) = (Needle m, ℝ)
-  CD¹ h₀ m₀ .+~^ (h₁δm, δh)
-      = let h₁ = min 1 . max 1e-300 $ h₀+δh; δm = h₁δm^/h₁
-        in CD¹ h₁ (m₀.+~^δm)
-instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)
-             => PseudoAffine (CD¹ m) where
-  CD¹ h₁ m₁ .-~. CD¹ h₀ m₀
-     = fmap ( \δm -> (h₁*^δm, h₁-h₀) ) $ m₁.-~.m₀
+-- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)
+--              => Semimanifold (CD¹ m) where
+--   type Needle (CD¹ m) = (Needle m, ℝ)
+--   CD¹ h₀ m₀ .+~^ (h₁δm, δh)
+--       = let h₁ = min 1 . max 1e-300 $ h₀+δh; δm = h₁δm^/h₁
+--         in CD¹ h₁ (m₀.+~^δm)
+-- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)
+--              => PseudoAffine (CD¹ m) where
+--   CD¹ h₁ m₁ .-~. CD¹ h₀ m₀
+--      = fmap ( \δm -> (h₁*^δm, h₁-h₀) ) $ m₁.-~.m₀
                                
 
 
@@ -403,6 +540,9 @@
            in (z, f'*.*g', devfg)
 
 
+instance (RealDimension s) => EnhancedCat (->) (Differentiable s) where
+  arr (Differentiable f) x = let (y,_,_) = f x in y
+
 instance (MetricScalar s) => Cartesian (Differentiable s) where
   type UnitObject (Differentiable s) = ZeroDim s
   swap = Differentiable $ \(x,y) -> ((y,x), lSwap, const zeroV)
@@ -694,6 +834,11 @@
 
 instance (RealDimension s) => EnhancedCat (PWDiffable s) (Differentiable s) where
   arr = globalDiffable
+instance (RealDimension s) => EnhancedCat (->) (PWDiffable s) where
+  arr (PWDiffable g) x = let (_,Differentiable f) = g x
+                             (y,_,_) = f x 
+                         in y
+
                 
 instance (RealDimension s) => Cartesian (PWDiffable s) where
   type UnitObject (PWDiffable s) = ZeroDim s
diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs
new file mode 100644
--- /dev/null
+++ b/Data/Manifold/Riemannian.hs
@@ -0,0 +1,244 @@
+-- |
+-- Module      : Data.Manifold.Riemannian
+-- Copyright   : (c) Justus Sagemüller 2015
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+-- Riemannian manifolds are manifolds equipped with a 'Metric' at each point.
+-- That means, these manifolds aren't merely topological objects anymore, but
+-- have a geometry as well. This gives, in particular, a notion of distance
+-- and shortest paths (geodesics) along which you can interpolate.
+-- 
+-- Keep in mind that the types in this library are
+-- generally defined in an abstract-mathematical spirit, which may not always
+-- match the intuition if you think about manifolds as embedded in ℝ³.
+-- (For instance, the torus inherits its geometry from the decomposition as
+-- @'S¹' × 'S¹'@, not from the “doughnut” embedding; the cone over @S¹@ is
+-- simply treated as the unit disk, etc..)
+
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE DeriveFunctor              #-}
+{-# LANGUAGE DeriveFoldable             #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE FunctionalDependencies     #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE GADTs                      #-}
+{-# LANGUAGE RankNTypes                 #-}
+{-# LANGUAGE TupleSections              #-}
+{-# LANGUAGE ParallelListComp           #-}
+{-# LANGUAGE UnicodeSyntax              #-}
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE PatternGuards              #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE LiberalTypeSynonyms        #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
+
+
+module Data.Manifold.Riemannian  where
+
+
+import Data.List hiding (filter, all, elem, sum)
+import Data.Maybe
+import qualified Data.Map as Map
+import qualified Data.Vector as Arr
+import Data.List.NonEmpty (NonEmpty(..))
+import Data.List.FastNub
+import qualified Data.List.NonEmpty as NE
+import Data.Semigroup
+import Data.Ord (comparing)
+import Control.DeepSeq
+
+import Data.VectorSpace
+import Data.LinearMap
+import Data.LinearMap.HerMetric
+import Data.LinearMap.Category
+import Data.AffineSpace
+import Data.Basis
+import Data.Complex hiding (magnitude)
+import Data.Void
+import Data.Tagged
+import Data.Proxy
+
+import Data.Manifold.Types
+import Data.Manifold.Types.Primitive ((^), embed, coEmbed)
+import Data.Manifold.PseudoAffine
+import Data.VectorSpace.FiniteDimensional
+    
+import Data.Embedding
+import Data.CoNat
+
+import qualified Prelude as Hask hiding(foldl, sum, sequence)
+import qualified Control.Applicative as Hask
+import qualified Control.Monad       as Hask hiding(forM_, sequence)
+import Data.Functor.Identity
+import Control.Monad.Trans.State
+import Control.Monad.Trans.Writer
+import Control.Monad.Trans.Class
+import qualified Data.Foldable       as Hask
+import Data.Foldable (all, elem, toList, sum)
+import qualified Data.Traversable as Hask
+import Data.Traversable (forM)
+
+import qualified Numeric.LinearAlgebra.HMatrix as HMat
+
+import Control.Category.Constrained.Prelude hiding
+     ((^), all, elem, sum, forM, Foldable(..), Traversable)
+import Control.Arrow.Constrained
+import Control.Monad.Constrained hiding (forM)
+import Data.Foldable.Constrained
+
+import GHC.Generics (Generic)
+
+
+class PseudoAffine x => Geodesic x where
+  geodesicBetween ::
+          x -- ^ Starting point; the interpolation will yield this at -1.
+       -> x -- ^ End point, for +1.
+            -- 
+            --   If the two points are actually connected by a path...
+       -> Option (D¹ -> x) -- ^ ...then this is the interpolation function. Attention: 
+                           --   the type will change to 'Differentiable' in the future.
+
+interpolate :: (Geodesic x, IntervalLike i) => x -> x -> Option (i -> x)
+interpolate a b = (. toClosedInterval) <$> geodesicBetween a b
+
+
+
+
+#define deriveAffineGD(x)                                         \
+instance Geodesic x where {                                        \
+  geodesicBetween a b = return $ alerp a b . (/2) . (+1) . xParamD¹ \
+ }
+
+deriveAffineGD (ℝ)
+
+instance Geodesic (ZeroDim ℝ) where
+  geodesicBetween Origin Origin = return $ \_ -> Origin
+
+instance (Geodesic a, Geodesic b) => Geodesic (a,b) where
+  geodesicBetween (a,b) (α,β) = liftA2 (&&&) (geodesicBetween a α) (geodesicBetween b β)
+
+instance (Geodesic a, Geodesic b, Geodesic c) => Geodesic (a,b,c) where
+  geodesicBetween (a,b,c) (α,β,γ)
+      = liftA3 (\ia ib ic t -> (ia t, ib t, ic t))
+           (geodesicBetween a α) (geodesicBetween b β) (geodesicBetween c γ)
+
+instance (KnownNat n) => Geodesic (FreeVect n ℝ) where
+  geodesicBetween (FreeVect v) (FreeVect w)
+      = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2
+                            in FreeVect $ Arr.zipWith (\vi wi -> μv*vi + μw*wi) v w
+
+instance (PseudoAffine v) => Geodesic (FinVecArrRep t v ℝ) where
+  geodesicBetween (FinVecArrRep v) (FinVecArrRep w)
+   | HMat.size v>0 && HMat.size w>0
+      = return $ \(D¹ t) -> let μv = (1-t)/2; μw = (t+1)/2
+                            in FinVecArrRep $ HMat.scale μv v + HMat.scale μw w
+
+instance (Geodesic v, WithField ℝ HilbertSpace v)
+             => Geodesic (Stiefel1 v) where
+  geodesicBetween (Stiefel1 p') (Stiefel1 q')
+      = (\f -> \(D¹ t) -> Stiefel1 . f . D¹ $ g * tan (ϑ*t))
+            <$> geodesicBetween p q
+   where p = normalized p'; q = normalized q'
+         l = magnitude $ p^-^q
+         ϑ = asin $ l/2
+         g = sqrt $ 4/l^2 - 1
+
+
+instance Geodesic S⁰ where
+  geodesicBetween PositiveHalfSphere PositiveHalfSphere = return $ const PositiveHalfSphere
+  geodesicBetween NegativeHalfSphere NegativeHalfSphere = return $ const NegativeHalfSphere
+  geodesicBetween _ _ = Hask.empty
+
+instance Geodesic S¹ where
+  geodesicBetween (S¹ φ) (S¹ ϕ)
+    | abs (φ-ϕ) < pi  = (>>> S¹) <$> geodesicBetween φ ϕ
+    | φ > 0           = (>>> S¹ . \ψ -> signum ψ*pi - ψ)
+                        <$> geodesicBetween (pi-φ) (-ϕ-pi)
+    | otherwise       = (>>> S¹ . \ψ -> signum ψ*pi - ψ)
+                        <$> geodesicBetween (-pi-φ) (pi-ϕ)
+
+
+instance Geodesic (Cℝay S⁰) where
+  geodesicBetween p q = (>>> fromℝ) <$> geodesicBetween (toℝ p) (toℝ q)
+   where toℝ (Cℝay h PositiveHalfSphere) = h
+         toℝ (Cℝay h NegativeHalfSphere) = -h
+         fromℝ x | x>0        = Cℝay x PositiveHalfSphere
+                 | otherwise  = Cℝay (-x) NegativeHalfSphere
+
+instance Geodesic (CD¹ S⁰) where
+  geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)
+   where toI (CD¹ h PositiveHalfSphere) = h
+         toI (CD¹ h NegativeHalfSphere) = -h
+         fromI x | x>0        = CD¹ x PositiveHalfSphere
+                 | otherwise  = CD¹ (-x) NegativeHalfSphere
+
+instance Geodesic (Cℝay S¹) where
+  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)
+   where fromP = fromInterior
+         toP w = case toInterior w of {Option (Just i) -> i}
+
+instance Geodesic (CD¹ S¹) where
+  geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q)
+   where toI (CD¹ h (S¹ φ)) = (h*cos φ, h*sin φ)
+         fromI (x,y) = CD¹ (sqrt $ x^2+y^2) (S¹ $ atan2 y x)
+
+instance Geodesic (Cℝay S²) where
+  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)
+   where fromP = fromInterior
+         toP w = case toInterior w of {Option (Just i) -> i}
+
+instance Geodesic (CD¹ S²) where
+  geodesicBetween p q = (>>> fromI) <$> geodesicBetween (toI p) (toI q :: ℝ³)
+   where toI (CD¹ h sph) = h *^ embed sph
+         fromI v = CD¹ (magnitude v) (coEmbed v)
+
+#define geoVSpCone(c,t)                                               \
+instance (c) => Geodesic (Cℝay (t)) where {                            \
+  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q) \
+   where { fromP (x,0) = Cℝay 0 x                                        \
+         ; fromP (x,h) = Cℝay h (x^/h)                                    \
+         ; toP (Cℝay h w) = ( h*^w, h ) } } ;                              \
+instance (c) => Geodesic (CD¹ (t)) where {                                  \
+  geodesicBetween p q = (>>> fromP) <$> geodesicBetween (toP p) (toP q)      \
+   where { fromP (x,0) = CD¹ 0 x                                              \
+         ; fromP (x,h) = CD¹ h (x^/h)                                          \
+         ; toP (CD¹ h w) = ( h*^w, h ) } }
+
+geoVSpCone ((), ℝ)
+geoVSpCone ((), ℝ⁰)
+geoVSpCone ((WithField ℝ HilbertSpace a, WithField ℝ HilbertSpace b, Geodesic (a,b)), (a,b))
+geoVSpCone (KnownNat n, FreeVect n ℝ)
+geoVSpCone ((Geodesic v, WithField ℝ HilbertSpace v), FinVecArrRep t v ℝ)
+
+
+
+
+-- | One-dimensional manifolds, whose closure is homeomorpic to the unit interval.
+class WithField ℝ PseudoAffine i => IntervalLike i where
+  toClosedInterval :: i -> D¹ -- Differentiable ℝ i D¹
+
+instance IntervalLike D¹ where
+  toClosedInterval = id
+instance IntervalLike (CD¹ S⁰) where
+  toClosedInterval (CD¹ h PositiveHalfSphere) = D¹ h
+  toClosedInterval (CD¹ h NegativeHalfSphere) = D¹ (-h)
+instance IntervalLike (Cℝay S⁰) where
+  toClosedInterval (Cℝay h PositiveHalfSphere) = D¹ $ tanh h
+  toClosedInterval (Cℝay h NegativeHalfSphere) = D¹ $ -tanh h
+instance IntervalLike (CD¹ ℝ⁰) where
+  toClosedInterval (CD¹ h Origin) = D¹ $ h*2 - 1
+instance IntervalLike (Cℝay ℝ⁰) where
+  toClosedInterval (Cℝay h Origin) = D¹ $ 1 - 2/(h+1)
+instance IntervalLike ℝ where
+  toClosedInterval x = D¹ $ tanh x
+
diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs
--- a/Data/Manifold/TreeCover.hs
+++ b/Data/Manifold/TreeCover.hs
@@ -34,7 +34,9 @@
 
 module Data.Manifold.TreeCover (
        -- * Shades 
-         Shade, shadeCtr, shadeExpanse, fullShade, pointsShades
+         Shade(..), Shade'(..)
+       -- ** Lenses and constructors
+       , shadeCtr, shadeExpanse, shadeNarrowness, fullShade, fullShade', pointsShades
        -- * Shade trees
        , ShadeTree(..), fromLeafPoints
        -- * Simple view helpers
@@ -42,14 +44,14 @@
        -- ** Auxiliary types
        , SimpleTree, Trees, NonEmptyTree, GenericTree(..)
        -- * Misc
-       , sShSaw, chainsaw, HasFlatView(..)
+       , sShSaw, chainsaw, HasFlatView(..), shadesMerge, smoothInterpolate
        -- ** Triangulation-builders
        , TriangBuild, doTriangBuild, singleFullSimplex, autoglueTriangulation
        , AutoTriang, elementaryTriang, breakdownAutoTriang
     ) where
 
 
-import Data.List hiding (filter, all, elem, sum)
+import Data.List hiding (filter, all, elem, sum, foldr1)
 import Data.Maybe
 import qualified Data.Map as Map
 import qualified Data.Vector as Arr
@@ -87,14 +89,14 @@
 import Control.Monad.Trans.Writer
 import Control.Monad.Trans.Class
 import qualified Data.Foldable       as Hask
-import Data.Foldable (all, elem, toList, sum)
+import Data.Foldable (all, elem, toList, sum, foldr1)
 import qualified Data.Traversable as Hask
 import Data.Traversable (forM)
 
 import qualified Numeric.LinearAlgebra.HMatrix as HMat
 
 import Control.Category.Constrained.Prelude hiding
-     ((^), all, elem, sum, forM, Foldable(..), Traversable)
+     ((^), all, elem, sum, forM, Foldable(..), foldr1, Traversable)
 import Control.Arrow.Constrained
 import Control.Monad.Constrained hiding (forM)
 import Data.Foldable.Constrained
@@ -116,17 +118,48 @@
 -- 
 --   For a /precise/ description of an arbitrarily-shaped connected subset of a manifold,
 --   there is 'Region', whose implementation is vastly more complex.
-data Shade x = Shade { shadeCtr :: !x
-                     , shadeExpanse :: !(Metric' x) }
+data Shade x = Shade { _shadeCtr :: !(Interior x)
+                     , _shadeExpanse :: !(Metric' x) }
 
+-- | A &#x201c;co-shade&#x201d; can describe ellipsoid regions as well, but unlike
+--   'Shade' it can be unlimited / infinitely wide in some directions.
+--   It does OTOH need to have nonzero thickness, which 'Shade' needs not.
+data Shade' x = Shade' { _shade'Ctr :: !(Interior x)
+                       , _shade'Narrowness :: !(Metric x) }
+
+class IsShade shade where
+--  type (*) shade :: *->*
+  -- | Access the center of a 'Shade' or a 'Shade''.
+  shadeCtr :: Functor f (->) (->) => (Interior x->f (Interior x)) -> shade x -> f (shade x)
+--  -- | Convert between 'Shade' and 'Shade' (which must be neither singular nor infinite).
+--  unsafeDualShade :: WithField ℝ Manifold x => shade x -> shade* x
+
+instance IsShade Shade where
+  shadeCtr f (Shade c e) = fmap (`Shade`e) $ f c
+
+shadeExpanse :: Functor f (->) (->) => (Metric' x -> f (Metric' x)) -> Shade x -> f (Shade x)
+shadeExpanse f (Shade c e) = fmap (Shade c) $ f e
+
+instance IsShade Shade' where
+  shadeCtr f (Shade' c e) = fmap (`Shade'`e) $ f c
+
+shadeNarrowness :: Functor f (->) (->) => (Metric x -> f (Metric x)) -> Shade' x -> f (Shade' x)
+shadeNarrowness f (Shade' c e) = fmap (Shade' c) $ f e
+
 instance (AffineManifold x) => Semimanifold (Shade x) where
   type Needle (Shade x) = Diff x
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   Shade c e .+~^ v = Shade (c.+^v) e
   Shade c e .-~^ v = Shade (c.-^v) e
 
 fullShade :: WithField ℝ Manifold x => x -> Metric' x -> Shade x
 fullShade ctr expa = Shade ctr expa
 
+fullShade' :: WithField ℝ Manifold x => x -> Metric x -> Shade' x
+fullShade' ctr expa = Shade' ctr expa
+
 subshadeId' :: WithField ℝ Manifold x
                    => x -> NonEmpty (DualSpace (Needle x)) -> x -> (Int, HourglassBulb)
 subshadeId' c expvs x = case x .-~. c of
@@ -151,6 +184,7 @@
 pointsShades :: WithField ℝ Manifold x => [x] -> [Shade x]
 pointsShades = map snd . pointsShades' zeroV
 
+
 pseudoECM :: WithField ℝ Manifold x => NonEmpty x -> (x, ([x],[x]))
 pseudoECM (p₀ NE.:| psr) = foldl' ( \(acc, (rb,nr)) (i,p)
                                   -> case p.-~.acc of 
@@ -171,22 +205,52 @@
               <$> mapM (.-~.ctr) ps
        
 
-minusLogOcclusion :: (PseudoAffine x, HasMetric (Needle x)
-             , s ~ (Scalar (Needle x)), RealDimension s )
+-- | Attempt to reduce the number of shades to fewer (ideally, a single one).
+--   In the simplest cases these should guaranteed cover the same area;
+--   for non-flat manifolds it only works in a heuristic sense.
+shadesMerge :: WithField ℝ Manifold x
+                 => ℝ -- ^ How near (inverse normalised distance, relative to shade expanse)
+                      --   two shades must be to be merged. If this is zero, any shades
+                      --   in the same connected region of a manifold are merged.
+                 -> [Shade x] -- ^ A list of /n/ shades.
+                 -> [Shade x] -- ^ /m/ &#x2264; /n/ shades which cover at least the same area.
+shadesMerge fuzz (sh₁@(Shade c₁ e₁) : shs) = case extractJust tryMerge shs of
+          (Just mg₁, shs') -> shadesMerge fuzz
+                                $ shs'++[mg₁] -- Append to end to prevent undue weighting
+                                              -- of first shade and its mergers.
+          (_, shs') -> sh₁ : shadesMerge fuzz shs' 
+ where tryMerge (Shade c₂ e₂)
+           | Option (Just v) <- c₁.-~.c₂
+           , Option (Just v') <- c₂.-~.c₁
+           , [e₁',e₂'] <- recipMetric<$>[e₁, e₂] 
+           , b₁ <- metric e₂' v
+           , b₂ <- metric e₁' v
+           , fuzz*b₁*b₂ <= b₁ + b₂
+                  = Just $ let cc = c₂ .+~^ v ^/ 2
+                               Option (Just cv₁) = c₁.-~.cc
+                               Option (Just cv₂) = c₂.-~.cc
+                           in Shade cc . sumV $ [e₁, e₂] ++ projector'<$>[cv₁, cv₂] 
+           | otherwise  = Nothing
+shadesMerge _ shs = shs
+
+minusLogOcclusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )
                 => Shade x -> x -> s
 minusLogOcclusion (Shade p₀ δ) = occ
  where occ p = case p .-~. p₀ of
-         Option(Just vd) -> metricSq δinv vd
+         Option(Just vd) | mSq <- metricSq δinv vd
+                         , mSq == mSq  -- avoid NaN
+                         -> mSq
          _               -> 1/0
        δinv = recipMetric δ
   
 -- | Check the statistical likelyhood of a point being within a shade.
-occlusion :: (PseudoAffine x, HasMetric (Needle x)
-             , s ~ (Scalar (Needle x)), RealDimension s )
+occlusion :: ( Manifold x, s ~ (Scalar (Needle x)), RealDimension s )
                 => Shade x -> x -> s
 occlusion (Shade p₀ δ) = occ
  where occ p = case p .-~. p₀ of
-         Option(Just vd) -> exp . negate $ metricSq δinv vd
+         Option(Just vd) | mSq <- metricSq δinv vd
+                         , mSq == mSq  -- avoid NaN
+                         -> exp (negate mSq)
          _               -> zeroV
        δinv = recipMetric δ
 
@@ -215,11 +279,6 @@
 flipHour :: Hourglass s -> Hourglass s
 flipHour (Hourglass u l) = Hourglass l u
 
-newtype Hourglasses s = Hourglasses {
-             getHourglasses :: NonEmpty (Hourglass s) }
-    deriving (Generic, Hask.Functor, Hask.Foldable)
-instance (NFData s) => NFData (Hourglasses s)
-
 data HourglassBulb = UpperBulb | LowerBulb
 oneBulb :: HourglassBulb -> (a->a) -> Hourglass a->Hourglass a
 oneBulb UpperBulb f (Hourglass u l) = Hourglass (f u) l
@@ -256,6 +315,9 @@
 -- | Experimental. There might be a more powerful instance possible.
 instance (AffineManifold x) => Semimanifold (ShadeTree x) where
   type Needle (ShadeTree x) = Diff x
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   PlainLeaves xs .+~^ v = PlainLeaves $ (.+^v)<$>xs 
   OverlappingBranches n sh br .+~^ v
         = OverlappingBranches n (sh.+~^v)
@@ -279,25 +341,23 @@
    where ne (PlainLeaves []) = False; ne _ = True
 
 
--- | Build a really quite nicely balanced tree from a cloud of points, on
---   any real manifold.
+-- | Build a quite nicely balanced tree from a cloud of points, on any real manifold.
 -- 
 --   Example:
 -- 
 -- @
--- > :m +Graphics.Dynamic.Plot.R2 Data.Manifold.TreeCover Data.VectorSpace Data.AffineSpace
--- > import Diagrams.Prelude ((^&), P2, R2, circle, fc, (&), moveTo, green)
---  
--- > let testPts0 = [0^&0, 0^&1, 1^&1, 1^&2, 2^&2] :: [P2]  -- Generate sort-of&#x2013;random point cloud
--- > let testPts1 = [p .+^ v^/3 | p<-testPts0, v <- [0^&0, (-1)^&1, 1^&2]]
--- > let testPts2 = [p .+^ v^/4 | p<-testPts1, v <- [0^&0, (-1)^&1, 1^&2]]
--- > let testPts3 = [p .+^ v^/5 | p<-testPts2, v <- [0^&0, (-2)^&1, 1^&2]]
--- > let testPts4 = [p .+^ v^/7 | p<-testPts3, v <- [0^&1, (-2)^&1, 1^&2]]
--- > length testPts4
---     405
+-- > :m +Graphics.Dynamic.Plot.R2 Data.Manifold.TreeCover Data.VectorSpace Data.AffineSpace 
+-- > import Diagrams.Prelude ((^&), p2, r2, P2, circle, fc, (&), moveTo, opacity)
 -- 
--- > plotWindow [ plot . onlyNodes $ fromLeafPoints testPts4
--- >            , plot [circle 0.06 & moveTo p & fc green :: PlainGraphics | p <- testPts4] ]
+-- >   -- Generate sort-of&#x2013;random cloud of lots of points
+-- > let testPts0 = p2 \<$\> [(0,0), (0,1), (1,1), (1,2), (2,2)] :: [P2 Double]
+-- > let testPts1 = [p .+^ v^/3 | p\<-testPts0, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]
+-- > let testPts2 = [p .+^ v^/4 | p\<-testPts1, v \<- r2\<$\>[(0,0), (-1,1), (1,2)]]
+-- > let testPts3 = [p .+^ v^/5 | p\<-testPts2, v \<- r2\<$\>[(0,0), (-2,1), (1,2)]]
+-- > let testPts4 = [p .+^ v^/7 | p\<-testPts3, v \<- r2\<$\>[(0,1), (-1,1), (1,2)]]
+-- 
+-- > plotWindow [ plot [ shapePlot $ circle 0.06 & moveTo p & opacity 0.3 | p <- testPts4 ]
+-- >            , plot . onlyNodes $ 'fromLeafPoints' testPts4 ]
 -- @
 -- 
 -- <<images/examples/simple-2d-ShadeTree.png>>
@@ -311,7 +371,7 @@
                                          Just redBrchs
                                            -> OverlappingBranches
                                                   (length xs) rShade
-                                                  (branchProc (shadeExpanse rShade) redBrchs)
+                                                  (branchProc (_shadeExpanse rShade) redBrchs)
                                          _ -> PlainLeaves xs
                      partitions -> DisjointBranches (length xs)
                                    . NE.fromList
@@ -406,6 +466,9 @@
   BaryCoords v .+^ w = BaryCoords $ v ^+^ w
 instance (KnownNat n) => Semimanifold (BaryCoords n) where
   type Needle (BaryCoords n) = FreeVect n ℝ
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   (.+~^) = (.+^)
 instance (KnownNat n) => PseudoAffine (BaryCoords n) where
   (.-~.) = pure .: (.-.)
@@ -874,6 +937,109 @@
 sShSaw _ _ = error "`sShSaw` is not supposed to cut anything else but `OverlappingBranches`"
 
 
+
+-- | Essentially the same as @(x,y)@, but not considered as a product topology.
+--   The 'Semimanifold' etc. instances just copy the topology of @x@, ignoring @y@.
+data x`WithAny`y
+      = WithAny { _untopological :: y
+                , _topological :: !x  }
+ deriving (Hask.Functor)
+
+instance (Semimanifold x) => Semimanifold (x`WithAny`y) where
+  type Needle (WithAny x y) = Needle x
+  type Interior (WithAny x y) = Interior x `WithAny` y
+  WithAny y x .+~^ δx = WithAny y $ x.+~^δx
+  fromInterior (WithAny y x) = WithAny y $ fromInterior x
+  toInterior (WithAny y x) = fmap (WithAny y) $ toInterior x
+  translateP = tpWD
+   where tpWD :: ∀ x y . Semimanifold x => Tagged (WithAny x y)
+                            (Interior x`WithAny`y -> Needle x -> Interior x`WithAny`y)
+         tpWD = Tagged `id` \(WithAny y x) δx -> WithAny y $ tpx x δx
+          where Tagged tpx = translateP :: Tagged x (Interior x -> Needle x -> Interior x)
+            
+instance (PseudoAffine x) => PseudoAffine (x`WithAny`y) where
+  WithAny _ x .-~. WithAny _ ξ = x.-~.ξ
+
+instance (AffineSpace x) => AffineSpace (x`WithAny`y) where
+  type Diff (WithAny x y) = Diff x
+  WithAny _ x .-. WithAny _ ξ = x.-.ξ
+  WithAny y x .+^ δx = WithAny y $ x.+^δx 
+
+instance (VectorSpace x, Monoid y) => VectorSpace (x`WithAny`y) where
+  type Scalar (WithAny x y) = Scalar x
+  μ *^ WithAny y x = WithAny y $ μ*^x 
+
+instance (AdditiveGroup x, Monoid y) => AdditiveGroup (x`WithAny`y) where
+  zeroV = WithAny mempty zeroV
+  negateV (WithAny y x) = WithAny y $ negateV x
+  WithAny y x ^+^ WithAny υ ξ = WithAny (mappend y υ) (x^+^ξ)
+
+instance (AdditiveGroup x) => Hask.Applicative (WithAny x) where
+  pure x = WithAny x zeroV
+  WithAny f x <*> WithAny t ξ = WithAny (f t) (x^+^ξ)
+  
+instance (AdditiveGroup x) => Hask.Monad (WithAny x) where
+  return x = WithAny x zeroV
+  WithAny y x >>= f = WithAny r $ x^+^q
+   where WithAny r q = f y
+
+shadeWithoutAnything :: Shade (x`WithAny`y) -> Shade x
+shadeWithoutAnything (Shade (WithAny _ b) e) = Shade b e
+
+-- | This is to 'ShadeTree' as 'Data.Map.Map' is to 'Data.Set.Set'.
+type x`Shaded`y = ShadeTree (x`WithAny`y)
+
+stiWithDensity :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y)
+         => x`Shaded`y -> x -> Cℝay y
+stiWithDensity (PlainLeaves lvs)
+  | [locShape@(Shade baryc expa)] <- pointsShades $ _topological <$> lvs
+       = let nlvs = fromIntegral $ length lvs :: ℝ
+             indiShapes = [(Shade p expa, y) | WithAny y p <- lvs]
+         in \x -> let lcCoeffs = [ occlusion psh x | (psh, _) <- indiShapes ]
+                      dens = sum lcCoeffs
+                  in mkCone dens . linearCombo . zip (snd<$>indiShapes)
+                       $ (/dens)<$>lcCoeffs
+stiWithDensity (DisjointBranches _ lvs)
+           = \x -> foldr1 qGather $ (`stiWithDensity`x)<$>lvs
+ where qGather (Cℝay 0 _) o = o
+       qGather o _ = o
+stiWithDensity (OverlappingBranches n (Shade (WithAny _ bc) extend) brs) = ovbSWD
+ where ovbSWD x = case x .-~. bc of
+           Option (Just v)
+             | dist² <- metricSq ε v
+             , dist² < 9
+             , att <- exp(1/(dist²-9)+1/9)
+               -> qGather att $ fmap ($x) downPrepared
+           _ -> coneTip
+       ε = recipMetric extend
+       downPrepared = dp =<< brs
+        where dp (DBranch _ (Hourglass up dn))
+                 = fmap stiWithDensity $ up:|[dn]
+       qGather att contribs = mkCone (att*dens)
+                 $ linearCombo [(v, d/dens) | Cℝay d v <- NE.toList contribs]
+        where dens = sum (hParamCℝay <$> contribs)
+
+
+smoothInterpolate :: (WithField ℝ Manifold x, WithField ℝ LinearManifold y)
+             => NonEmpty (x,y) -> x -> y
+smoothInterpolate l = \x ->
+             case ltr x of
+               Cℝay 0 _ -> defy
+               Cℝay _ y -> y
+ where defy = linearCombo [(y, 1/n) | WithAny y _ <- l']
+       n = fromIntegral $ length l'
+       l' = (uncurry WithAny . swap) <$> NE.toList l
+       ltr = stiWithDensity $ fromLeafPoints l'
+
+
+coneTip :: (AdditiveGroup v) => Cℝay v
+coneTip = Cℝay 0 zeroV
+
+mkCone :: AdditiveGroup v => ℝ -> v -> Cℝay v
+mkCone 0 _ = coneTip
+mkCone h v = Cℝay h v
+
+
 foci :: [a] -> [(a,[a])]
 foci [] = []
 foci (x:xs) = (x,xs) : fmap (second (x:)) (foci xs)
@@ -909,4 +1075,14 @@
         ]
   superFlatView = foldMap go . flatView
    where go (t,ds) = t : ds
+
+
+
+
+
+
+extractJust :: (a->Maybe b) -> [a] -> (Maybe b, [a])
+extractJust f [] = (Nothing,[])
+extractJust f (x:xs) | Just r <- f x  = (Just r, xs)
+                     | otherwise      = second (x:) $ extractJust f xs
 
diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs
--- a/Data/Manifold/Types.hs
+++ b/Data/Manifold/Types.hs
@@ -72,6 +72,7 @@
 
 import Data.Manifold.Types.Primitive
 import Data.Manifold.PseudoAffine
+import Data.Manifold.Cone
 import Data.LinearMap.HerMetric
 import Data.VectorSpace.FiniteDimensional
 
@@ -85,20 +86,14 @@
 #define deriveAffine(c,t)                \
 instance (c) => Semimanifold (t) where {  \
   type Needle (t) = Diff (t);              \
-  (.+~^) = (.+^) };                         \
-instance (c) => PseudoAffine (t) where {     \
+  fromInterior = id;                        \
+  toInterior = pure;                         \
+  translateP = Tagged (.+~^);                 \
+  (.+~^) = (.+^) };                            \
+instance (c) => PseudoAffine (t) where {        \
   a.-~.b = pure (a.-.b);      }
 
 
--- | The /n/-th Stiefel manifold is the space of all possible configurations of
---   /n/ orthonormal vectors. In the case /n/ = 1, simply the subspace of normalised
---   vectors, i.e. equivalent to the 'UnitSphere'. Even so, it strictly speaking
---   requires the containing space to be at least metric (if not Hilbert); we would
---   however like to be able to use this concept also in spaces with no inner product,
---   therefore we define this space not as normalised vectors, but rather as all
---   vectors modulo scaling by positive factors.
-newtype Stiefel1 v = Stiefel1 { getStiefel1N :: DualSpace v }
-
 newtype Stiefel1Needle v = Stiefel1Needle { getStiefel1Tangent :: HMat.Vector (Scalar v) }
 newtype Stiefel1Basis v = Stiefel1Basis { getStiefel1Basis :: Int }
 s1bTrie :: forall v b. FiniteDimensional v => (Stiefel1Basis v->b) -> Stiefel1Basis v:->:b
@@ -163,19 +158,12 @@
 
 instance (WithField k LinearManifold v, Real k) => Semimanifold (Stiefel1 v) where 
   type Needle (Stiefel1 v) = Stiefel1Needle v
+  fromInterior = id
+  toInterior = pure
+  translateP = Tagged (.+~^)
   Stiefel1 s .+~^ Stiefel1Needle n = Stiefel1 . fromPackedVector . HMat.scale (signum s'i)
    $ if| ν==0      -> s' -- ν'≡0 is a special case of this, so we can otherwise assume ν'>0.
--- --  | ν<=1      -> let -- κ = (-1 − 1/(ν−1)) / ν'
---                        -- m ∝         spro +         κ · n
---                        --   ∝ (1−ν) · spro + (1−ν) · κ · n
---                        --   = (1−ν) · spro + (-(1−ν) − -1)/ν' · n
---                        m = HMat.scale (1-ν) spro + HMat.scale (ν/ν') n
---                    in insi (1-ν) m
-       | ν<=2      -> let -- κ = (1/(ν−1) − 1) / ν'
-                          -- m ∝       - spro +         κ · n
-                          --   ∝ (1−ν) · spro + (ν−1) · κ · n
-                          --   = (1−ν) · spro + (1 − (ν−1))/ν' · n
-                          m = HMat.scale ιmν spro + HMat.scale ((1-abs ιmν)/ν') n
+       | ν<=2      -> let m = HMat.scale ιmν spro + HMat.scale ((1-abs ιmν)/ν') n
                           ιmν = 1-ν 
                       in insi ιmν m
        | otherwise -> let m = HMat.scale ιmν spro + HMat.scale ((abs ιmν-1)/ν') n
@@ -199,8 +187,7 @@
             s'i | v <- HMat.scale (recip s'i) delis - tpro
                 , absv <- l2norm v
                 , absv > 0
-                       -> let μ -- = (1 − recip (|v| + 1)) / |v| for sgn sᵢ = sgn tᵢ
-                                   = (signum (t'i/s'i) - recip(absv + 1)) / absv
+                       -> let μ = (signum (t'i/s'i) - recip(absv + 1)) / absv
                           in HMat.scale μ v
                 | t'i/s'i > 0  -> samePoint
                 | otherwise    -> antipode
@@ -214,34 +201,21 @@
          samePoint = (d-1) HMat.|> repeat 0
          antipode = (d-1) HMat.|> (2 : repeat 0)
 
-l2norm :: MetricScalar s => HMat.Vector s -> s
-l2norm = realToFrac . HMat.norm_2
 
+instance ( WithField ℝ HilbertSpace x ) => ConeSemimfd (Stiefel1 x) where
+  type CℝayInterior (Stiefel1 x) = x
+  fromCℝayInterior (FinVecArrRep v) = case HMat.size v of
+      0 -> Cℝay 0 $ Stiefel1 zeroV
+      _ -> Cℝay (HMat.norm_2 v) $ Stiefel1 (fromPackedVector v)
+  toCℝayInterior (Cℝay 0 _) = pure zeroV
+  toCℝayInterior (Cℝay l (Stiefel1 v))
+        = pure.FinVecArrRep $ HMat.scale (l/HMat.norm_2 v') v'
+   where v' = asPackedVector v
 
-stiefel1Project :: LinearManifold v =>
-             DualSpace v       -- ^ Must be nonzero.
-                 -> Stiefel1 v
-stiefel1Project = Stiefel1
 
-stiefel1Embed :: HilbertSpace v => Stiefel1 v -> v
-stiefel1Embed (Stiefel1 n) = normalized n
-  
-
-class (PseudoAffine v, InnerSpace v, NaturallyEmbedded (UnitSphere v) (DualSpace v))
-          => HasUnitSphere v where
-  type UnitSphere v :: *
-  stiefel :: UnitSphere v -> Stiefel1 v
-  stiefel = Stiefel1 . embed
-  unstiefel :: Stiefel1 v -> UnitSphere v
-  unstiefel = coEmbed . getStiefel1N
-
-instance HasUnitSphere ℝ  where type UnitSphere ℝ  = S⁰
-instance HasUnitSphere ℝ² where type UnitSphere ℝ² = S¹
-instance HasUnitSphere ℝ³ where type UnitSphere ℝ³ = S²
+l2norm :: MetricScalar s => HMat.Vector s -> s
+l2norm = realToFrac . HMat.norm_2
 
-instance (HasUnitSphere v, v ~ DualSpace v) => NaturallyEmbedded (Stiefel1 v) v where
-  embed = embed . unstiefel
-  coEmbed = stiefel . coEmbed
 
 
 
diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs
--- a/Data/Manifold/Types/Primitive.hs
+++ b/Data/Manifold/Types/Primitive.hs
@@ -37,7 +37,7 @@
         , Projective1, Projective2
         , Disk1, Disk2, Cone, OpenCone
         -- * Linear manifolds
-        , ZeroDim(..)
+        , ZeroDim(..), isoAttachZeroDim
         , ℝ⁰, ℝ, ℝ², ℝ³
         -- * Hyperspheres
         , S⁰(..), S¹(..), S²(..)
@@ -49,7 +49,7 @@
         , CD¹(..), Cℝay(..)
         -- * Utility (deprecated)
         , NaturallyEmbedded(..)
-        , GraphWindowSpec(..), Endomorphism, (^), EqFloating
+        , GraphWindowSpec(..), Endomorphism, (^), (^.), EqFloating
    ) where
 
 
@@ -60,6 +60,8 @@
 import Data.Void
 import Data.Monoid
 
+import Control.Applicative (Const(..))
+
 import qualified Prelude
 
 import Control.Category.Constrained.Prelude hiding ((^))
@@ -67,11 +69,13 @@
 import Control.Monad.Constrained
 import Data.Foldable.Constrained
 
+import Data.Embedding
 
 
 
 
 
+
 type EqFloating f = (Eq f, Ord f, Floating f)
 
 
@@ -101,6 +105,13 @@
   decompose Origin = []
   decompose' Origin = absurd
 
+{-# INLINE isoAttachZeroDim #-}
+isoAttachZeroDim :: ( WellPointed c, UnitObject c ~ (), ObjectPair c a ()
+                    , Object c (ZeroDim k), ObjectPair c a (ZeroDim k)
+                    , PointObject c (ZeroDim k) )
+                       => Isomorphism c a (a, ZeroDim k)
+isoAttachZeroDim = second (Isomorphism (const Origin) terminal) . attachUnit
+
 -- | The zero-dimensional sphere is actually just two points. Implementation might
 --   therefore change to @ℝ⁰ 'Control.Category.Constrained.+' ℝ⁰@: the disjoint sum of two
 --   single-point spaces.
@@ -247,7 +258,14 @@
   () <.> () = 0
 
 
+infixr 8 ^
 
 (^) :: Num a => a -> Int -> a
 (^) = (Prelude.^)
+
+
+infixl 8 ^.
+{-# INLINE (^.) #-}
+(^.) :: s -> (forall f . Prelude.Functor f => (a->f a) -> s->f s) -> a
+o ^. g = getConst (g Const o)
 
diff --git a/Data/VectorSpace/FiniteDimensional.hs b/Data/VectorSpace/FiniteDimensional.hs
--- a/Data/VectorSpace/FiniteDimensional.hs
+++ b/Data/VectorSpace/FiniteDimensional.hs
@@ -10,13 +10,16 @@
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE MultiParamTypeClasses      #-}
 {-# LANGUAGE TypeOperators              #-}
 {-# LANGUAGE TupleSections              #-}
 {-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE PolyKinds                  #-}
 {-# LANGUAGE UndecidableInstances       #-}
 {-# LANGUAGE StandaloneDeriving         #-}
 {-# LANGUAGE ConstraintKinds            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE UnicodeSyntax              #-}
 
 
 
@@ -24,6 +27,7 @@
 module Data.VectorSpace.FiniteDimensional (
     FiniteDimensional(..)
   , SmoothScalar 
+  , FinVecArrRep(..), concreteArrRep, (⊗), splitArrRep
   ) where
     
 
@@ -31,6 +35,7 @@
 
 import Prelude hiding ((^))
 
+import Data.AffineSpace
 import Data.VectorSpace
 import Data.LinearMap
 import Data.Basis
@@ -42,7 +47,10 @@
     
 import Data.Manifold.Types.Primitive
 import Data.CoNat
+import Data.Embedding
 
+import Control.Arrow
+
 import qualified Data.Vector as Arr
 import qualified Numeric.LinearAlgebra.HMatrix as HMat
 
@@ -55,7 +63,6 @@
                       , Num(HMat.Vector s), HMat.Indexable(HMat.Vector s)s
                       , HMat.Normed(HMat.Vector s) )
 
-
 -- | Many linear algebra operations are best implemented via packed, dense 'HMat.Matrix'es.
 --   For one thing, that makes common general vector operations quite efficient,
 --   in particular on high-dimensional spaces.
@@ -160,4 +167,100 @@
   fromPackedVector arr = FreeVect (Arr.convert arr)
   -- asPackedMatrix = _ -- could be done quite efficiently here!
                                                           
+
+
+-- | Semantically the same as @'Tagged' tag refvs@, but directly uses the
+--   packed-vector array representation.
+-- 
+--   The tag should really be kind-polymorphic, but at least GHC-7.8 doesn't quite
+--   handle the associated types of the manifold classes then.
+newtype FinVecArrRep (tag :: * -> *) refvs scalar
+      = FinVecArrRep { getFinVecArrRep :: HMat.Vector scalar }
+
+instance (SmoothScalar s) => AffineSpace (FinVecArrRep t b s) where
+  type Diff (FinVecArrRep t b s) = FinVecArrRep t b s
+  (.-.) = (^-^)
+  (.+^) = (^+^)
+  
+instance (SmoothScalar s) => AdditiveGroup (FinVecArrRep t b s) where
+  zeroV = FinVecArrRep $ 0 HMat.|> []
+  negateV (FinVecArrRep v) = FinVecArrRep $ negate v
+  FinVecArrRep v ^+^ FinVecArrRep w
+   | HMat.size v == 0  = FinVecArrRep w
+   | HMat.size w == 0  = FinVecArrRep w
+   | otherwise         = FinVecArrRep $ v + w
+
+instance (SmoothScalar s) => VectorSpace (FinVecArrRep t b s) where
+  type Scalar (FinVecArrRep t b s) = s
+  μ *^ FinVecArrRep v = FinVecArrRep $ HMat.scale μ v
+
+instance (SmoothScalar s) => InnerSpace (FinVecArrRep t b s) where
+  FinVecArrRep v <.> FinVecArrRep w
+   | HMat.size v == 0  = 0
+   | HMat.size w == 0  = 0
+   | otherwise         = v`HMat.dot`w
+
+concreteArrRep :: (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)
+           => Isomorphism (->) r (FinVecArrRep t r s)
+concreteArrRep = Isomorphism (FinVecArrRep     . asPackedVector)
+                             (fromPackedVector . getFinVecArrRep)
+
+(⊗) :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w
+                   , Scalar v ~ s, Scalar w ~ s )
+          => FinVecArrRep t v s -> FinVecArrRep t w s -> FinVecArrRep t (v,w) s
+FinVecArrRep v ⊗ FinVecArrRep w
+  | HMat.size v + HMat.size w == 0  = FinVecArrRep v
+  | HMat.size v == 0                = FinVecArrRep $ HMat.vjoin [HMat.konst 0 nv, w]
+  | HMat.size w == 0                = FinVecArrRep $ HMat.vjoin [v, HMat.konst 0 nw]
+  | otherwise                       = FinVecArrRep $ HMat.vjoin [v,w]
+ where Tagged nv = dimension :: Tagged v Int
+       Tagged nw = dimension :: Tagged w Int
+
+splitArrRep :: ∀ t s v w . ( SmoothScalar s, FiniteDimensional v, FiniteDimensional w
+                   , Scalar v ~ s, Scalar w ~ s )
+          => FinVecArrRep t (v,w) s -> (FinVecArrRep t v s, FinVecArrRep t w s)
+splitArrRep (FinVecArrRep vw)
+  | HMat.size vw == 0   = (FinVecArrRep vw, FinVecArrRep vw)
+  | otherwise           = ( FinVecArrRep $ HMat.subVector 0 nv vw
+                          , FinVecArrRep $ HMat.subVector nv nw vw )
+ where Tagged nv = dimension :: Tagged v Int
+       Tagged nw = dimension :: Tagged w Int
+                  
+
+instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)
+                 => HasBasis (FinVecArrRep t r s) where
+  type Basis (FinVecArrRep t r s) = Basis r
+  basisValue = (concreteArrRep$->$) . basisValue
+  decompose = decompose . (concreteArrRep$<-$)
+  decompose' = decompose' . (concreteArrRep$<-$)
+
+instance (SmoothScalar s, FiniteDimensional r, Scalar r ~ s)
+                 => FiniteDimensional (FinVecArrRep t r s) where
+  dimension = d
+   where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) Int
+         d = Tagged n
+          where Tagged n = dimension :: Tagged r Int
+  indexBasis = d
+   where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Int -> Basis r)
+         d = Tagged n
+          where Tagged n = indexBasis :: Tagged r (Int -> Basis r)
+  basisIndex = d
+   where d :: ∀ t r s . FiniteDimensional r => Tagged (FinVecArrRep t r s) (Basis r -> Int)
+         d = Tagged n
+          where Tagged n = basisIndex :: Tagged r (Basis r -> Int)
+  asPackedVector = apv
+   where apv :: ∀ t r s . (FiniteDimensional r, SmoothScalar s)
+                     => FinVecArrRep t r s -> HMat.Vector s
+         apv (FinVecArrRep v)
+             | HMat.size v == 0  = HMat.konst 0 n
+             | otherwise         = v
+          where Tagged n = dimension :: Tagged r Int
+  fromPackedVector = FinVecArrRep
+
+
+instance (NaturallyEmbedded m r, FiniteDimensional r, s ~ Scalar r)
+                 => NaturallyEmbedded m (FinVecArrRep t r s) where
+  embed = (concreteArrRep$<-$) . embed
+  coEmbed = coEmbed . (concreteArrRep$->$)
+                     
 
diff --git a/images/examples/simple-2d-ShadeTree.png b/images/examples/simple-2d-ShadeTree.png
Binary files a/images/examples/simple-2d-ShadeTree.png and b/images/examples/simple-2d-ShadeTree.png differ
diff --git a/manifolds.cabal b/manifolds.cabal
--- a/manifolds.cabal
+++ b/manifolds.cabal
@@ -1,8 +1,8 @@
 Name:                manifolds
-Version:             0.1.3.1
+Version:             0.1.5.0
 Category:            Math
 Synopsis:            Working with manifolds in a direct, embedding-free way.
-Description:         Manifolds, a generalisation of the notion of \"smooth curves\" or sufaces,
+Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,
                      are topological spaces /locally homeomorphic to a vector space/. This gives
                      rise to what is actually the most natural / mathematically elegant way of dealing
                      with them: calculations can be carried out locally, in connection with Riemannian
@@ -41,7 +41,6 @@
                      , vector-space>=0.8
                      , MemoTrie
                      , vector
-                     , vector-algorithms
                      , hmatrix >= 0.16 && < 0.18
                      , containers
                      , comonad
@@ -68,8 +67,11 @@
                      Data.LinearMap.HerMetric
                      -- Data.Manifold.Visualisation.R3.GLUT
                      Data.Manifold.Types
+                     Data.Manifold.Griddable
+                     Data.Manifold.Riemannian
   Other-modules:   Data.List.FastNub
                    Data.Manifold.Types.Primitive
+                   Data.Manifold.Cone
                    Data.CoNat
                    Data.Embedding
                    Data.LinearMap.Category