gruff-0.1: PeriodScan.hs

module PeriodScan (periodScan, periodNucleus) where

import Control.Parallel.Strategies (parMap, rseq)
import Data.List (genericIndex, nubBy)
import Data.Maybe (listToMaybe)
import Data.Function (on)
import Data.Vec (NearZero, nearZero)
import Text.FShow.Raw (DecimalFormat(nanTest, infTest))

import Fractal.RUFF.Types.Complex (Complex((:+)))

import Utils (safeLast)

periodScan :: (Floating r, Ord r) => r -> Complex r -> Maybe Integer
periodScan r c =
  let cs = [ c + (r:+r), c + (r:+(-r)), c + ((-r):+(-r)), c + ((-r):+r) ]
      zs = iterate (zipWith (\cc z -> z * z + cc) cs) [0,0,0,0]
  in  fmap fst . listToMaybe . dropWhile (not . straddlesOrigin . snd) . zip [0 ..] $ zs

periodNucleus1 :: (Floating r, Ord r) => Integer -> r -> Complex r -> [(r, Complex r)]
periodNucleus1 p r0 c0
  = map fst . filter (straddlesOrigin . snd)
  . parMap rseq (fmap (parMap rseq (\c -> iterate (\z->z * z + c) 0 `genericIndex` p)))
  . fmap (\(r, c) -> ((r, c), [c + (r:+r), c + ((-r):+r), c + ((-r):+(-r)), c + (r:+(-r))]))
  $( [ (r', c0 + d) | i <- [-1,1], j <- [-1,1 ], let d = ((r'*i) :+ (r'*j)) ]
  ++ [ (r', c0 + d) | i <- [  0 ], j <- [  0 ], let d = ((r'*i) :+ (r'*j)) ])
{-
   of
    [] -> []
    xs -> fmap fst . minimumBy (comparing $ offsetFromOrigin . snd) $ xs -}
  where r' = r0 / 2

periodNucleus :: (NearZero r, DecimalFormat r, Floating r, Ord r) => Integer -> r -> Complex r -> Maybe (Complex r)
periodNucleus p r0 c0 =
  let ok (r, (x:+y)) = all ok' [r, x, y]
      ok' x = not (nanTest x || infTest x)
  in fmap snd . safeLast . takeWhile ok . concat . takeWhile (not . null) . iterate (nubBy (approxEq `on` snd) . (uncurry (periodNucleus1 p) =<<)) $ [(r0, c0)]

straddlesOrigin :: (Ord r, Num r) => [Complex r] -> Bool
straddlesOrigin ps = odd . length . filter id . zipWith positiveReal ps $ (drop 1 ps ++ take 1 ps)

positiveReal :: (Ord r, Num r) => Complex r -> Complex r -> Bool
positiveReal (u:+v) (x:+y)
  | v < 0 && y < 0 = False
  | v > 0 && y > 0 = False
  | (u * (y - v) - v * (x - u)) * (y - v) > 0 = True
  | otherwise = False

{-
offsetFromOrigin :: (Floating r) => [Complex r] -> r
offsetFromOrigin = magnitude . sum
-}

approxEq :: (Num c, NearZero c) => c -> c -> Bool
approxEq w z = nearZero (w - z)