```{-# LANGUAGE CPP #-}

{- |

Maintainer : Joachim Breitner <mail@joachim-breitner.de>
Stability  : stable
Portability: portable
-}
module Optimisation.CirclePacking (
packCircles
-- * Example
-- \$example
) where

import Data.List (sortBy, find)
import Data.Ord (comparing)
import Data.Maybe (fromMaybe)
#else
import Language.Fay.Prelude

-- Not provided in Fay set, it seems
comparing :: (b -> Double) -> b -> b -> Ordering
comparing p x y = compare (p x) (p y)

fromMaybe     :: a -> Maybe a -> a
fromMaybe d x = case x of {Nothing -> d;Just v  -> v}
#endif

type Circle a = (Double, a)
type Coordinate = (Double, Double)
type PlacedCircle a = ((Double, a), Coordinate)
type TouchingCircles a = [(PlacedCircle a, PlacedCircle a)]

{- | 'packCircles' takes a list of circles and a function that yields the

It returns a list of all circles, in unspecified order, together with
coordinates such that they do not overlap but sit as tight as possible,
filling a large circle.

Finding the optimal solution to this is NP hard, so only
heuristics are feasible. This particular
implementation is neither very good nor very fast,
compared to the state of the art in research. Nevertheless
it is simple to use and gives visually acceptable results.

The heuristics begins by placing the largest circle first, and the
next-to-largest next to it. From then on it adds circles by considering all
points where the circle to be added would touch two circles but overlap with
none, and picks the one that is closest to the center of mass of the current
placements.

-}
packCircles :: (a -> Double) -> [a] -> [(a, (Double, Double))]
map (\t -> case t of (x,p) -> (snd x, p)) .
-- Forget the pairs of touching circles
fst .
-- Run the main algorithm
go .
-- Look at large circles last
map (\x -> (radiusFunction x, x))

-- Just for a nicer name

-- Place the tail, then try to place the head
go :: [Circle a] -> ([PlacedCircle a], TouchingCircles a)
go  [] = ([],[])
go  (c:cs) = case go cs of
(placed, pairs) -> case place c placed pairs of
(cp, newPairs) -> (cp : placed, newPairs ++ pairs)

place :: Circle a -> [PlacedCircle a] -> TouchingCircles a ->
(PlacedCircle a, TouchingCircles a)
place c [] _ = ((c, (0,0)), [])

place c [cp'@(c',(x,y))] _ =
let cp = (c, (x + radius c + radius c', y))
in (cp, [(cp, cp')])

place c placed pairs = (cp, newPairs)
where
newPairs = [ (cp, cp') | cp' <- placed, touching cp cp' ]
cp = (c, p)
p = fromMaybe (error "packCircles: The end of the real plane has been reached?") \$
find (\p' -> all (valid p') placed) \$
sortBy (comparing centerDistance)
[ p' | (c1, c2) <- pairs,
touching c1 c2,
p' <- near c1 c2
]
centerDistance (x,y) = sqrt ((centerx - x)**2  + (centery - y)**2)

centerx = sum [x * (radius c')**2 | (c',(x,_)) <- placed] / area
centery = sum [y * (radius c')**2 | (c',(_,y)) <- placed] / area
area = sum [(radius c')**2 | (c',_) <- placed]

valid (x1,y1) (c2,(x2,y2)) =
sqrt ((x2 - x1)**2  + (y2-y1)**2) >= (radius c + radius c2) * (1-eps)

touching (c1,(x1,y1)) (c2, (x2, y2)) =
sqrt ((x2 - x1)**2  + (y2-y1)**2) <= (radius c1 + radius c2) * (1+eps)

near (c1,(x1,y1)) (c2, (x2, y2)) = [(c1x,c1y), (c2x,c2y)]
where
base = sqrt ((x2 - x1)**2  + (y2 - y1)**2)
--From http://stackoverflow.com/a/11356687/946226
ad_length = (base**2 + lat1**2 - lat2**2)/(2 * base)
h = sqrt (abs (lat1**2 - ad_length**2))
dx = x1 + ad_length * (x2 - x1)/base
dy = y1 + ad_length * (y2 - y1)/base

c1x = dx + h * (y2 - y1) / base
c1y = dy - h * (x2 - x1) / base

c2x = dx - h * (y2 - y1) / base
c2y = dy + h * (x2 - x1) / base

eps :: Double
eps = 0.00001

{-\$example

The following code demonstrates how one can use 'packCircles' together with
the diagrams library:

>import Diagrams.Prelude
>import Diagrams.Backend.SVG.CmdLine
>
>import Optimisation.CirclePacking
>
>colorize = zipWith fc \$
>    cycle [red,blue,yellow,magenta,cyan,bisque,firebrick,indigo]
>
>objects = colorize \$
>    [ circle r  | r <- [0.1,0.2..1.6] ] ++
>    [ hexagon r | r <- [0.1,0.2..0.7] ] ++
>    [ decagon r | r <- [0.1,0.2..0.7] ]
>
>-- Just a approximation, diagram objects do not have an exact radius
>radiusApproximation o = maximum [ radius (e (CircleFrac alpha)) o | alpha <- [0,0.1..1.0]]
>
>main = defaultMain \$
>    position \$ map (\(o,(x,y)) -> (p2 (x,y),o)) \$