{-# LANGUAGE OverloadedStrings #-}
module Exercise10 where

import System.IO
import Data.Char
import qualified Data.Text as T
import Data.List (find)
import Data.Maybe (fromMaybe)
import Data.Map as Map (lookup)

import Turtle ( black, move, sit, turn, Ln, Turtle )

{-Wettbewerb imports-}
import Data.Complex
import Data.Array
import Data.Colour.RGBSpace.HSV
import Data.Colour.RGBSpace
import Codec.Picture
import Data.Word
import System.Random

-- data type for L-System
data Rule = Char :->: String                -- context-free and deterministic! 
  deriving (Eq,Show)

data LSystem = LSystem {
                  start :: String, 
                  rules :: [Rule]           -- constraint: unique left sides
                  }
  deriving Eq

instance Show LSystem where
  show (LSystem s r) = unlines $ ["Start: " ++ show s, "Rules: "] ++ map show r


apply :: Char -> Turtle -> Turtle
apply 'F' t = move 10 t
apply 'G' t = move 10 t
apply '+' t = turn 30 t
apply '-' t = turn (-30) t
apply '*' t = turn 15 t
apply '~' t = turn (-15) t
apply _ t = sit t

-- Use apply to convert movements of turtle to GL lines.
-- Try changing the color to red! 
execute :: LSystem -> Integer -> [Ln]
execute ls n = let (pen, ang, pnt, ln, bps) = lines (black, 0, 0, [],[]) $ expandLSystem ls n in ln
  where
    lines t [] = t 
    lines t (x:xs) = lines (apply x t) xs

-- sample LSystems for (manual) testing
emptySystem :: LSystem
emptySystem = LSystem [] []

dragoncurve :: LSystem
dragoncurve = LSystem "FX" ['X' :->: "X+++YF+++", 'Y' :->: "---FX---Y"]

kochcurve :: LSystem
kochcurve = LSystem "F" ['F' :->: "F+++F---F---F+++F"]

sierpinski :: LSystem
sierpinski = LSystem "F++++G++++G" ['F' :->: "F++++G----F----G++++F", 'G' :->: "GG"]

-- finds the first occurrence of a fitting rule
findRule :: [Rule] -> Char -> Rule 
findRule [] c = c :->: [c]
findRule (r@(rc :->: _):rs) c
    | rc == c   = r
    | otherwise = findRule rs c

-- expands the L-System n times
expandLSystem :: LSystem -> Integer -> String
expandLSystem l = expandStrRec (start l) (rules l)
 
expandStrRec :: String -> [Rule] -> Integer -> String
expandStrRec s _ 0  = s
expandStrRec s rs depth = expandStrRec expandedString rs (depth - 1)
    where expandedString = expandStr s rs


expandStr :: String -> [Rule] -> String
expandStr [] rs     = []
expandStr (c:ss) rs = ns ++ expandStr ss rs
    where (cc :->: ns) = findRule rs c

-- updating LSystem via command
update :: LSystem -> IO LSystem
update ls@(LSystem st rls) 
    = do    isReady <- hReady stdin
            if isReady then
                do  input <- getLine
                    let (comType, tempCommand) = seperateCommand input
                        command = T.drop 1 tempCommand    -- removes the Whitespace in front of the command
                    case T.toLower comType of
                        "start" -> update $ LSystem (T.unpack command) rls
                        "rule"  -> do   let newRule = parseRule command
                                        case newRule of
                                            Nothing -> do   putStrLn "Error parsing rule"
                                                            update ls
                                            Just r  -> do   let newRules = updateRules r rls
                                                            update $ LSystem st newRules
                        "clear" -> update emptySystem
                        "print" -> do   putStrLn $ show ls
                                        update ls
                        _       -> do   putStrLn "Error parsing command"
                                        update ls
            else
                return ls

seperateCommand :: String -> (T.Text, T.Text)
seperateCommand = T.breakOn " " . T.pack

test = T.count " " "asd asd "

parseRule :: T.Text -> Maybe Rule
parseRule t
    -- If there is no " -> " present or the rhs is "", the command is wrong
    | T.null s          = Nothing
    -- According to the Zulip Post of Lukas Stevens, no Spaces are allowed
    | T.count " " s > 0 = Nothing 
    -- If the first argument is not a char, the command is wrong
    | T.length tC /= 1  = Nothing
    | otherwise         = Just $ c :->: T.unpack s
    where   (tC, tempS) = T.breakOn " -> " t
            s = T.drop 4 tempS  -- removes " -> "
            strC = T.unpack tC
            (c:_) = strC

updateRules :: Rule -> [Rule] -> [Rule]
updateRules r@(c :->: s) rs = r:removedRs
    where removedRs = removeSameCharRules c rs

-- Removes all rules of the List, which have the strucutre c :->: _
removeSameCharRules :: Char -> [Rule] -> [Rule]
removeSameCharRules _ [] = []
removeSameCharRules c (r@(ruleC :->: _):rs)
    | c == ruleC    = removeSameCharRules c rs
    | otherwise     = r : removeSameCharRules c rs
-- add the WETT ... TTEW tags if you want to participate in the wettbewerb!

{-MCCOMMENT
Wettbewerb File: ./wett/haskell.png (haskell_quadratic.png for lower resolution)

The created pictures are a possible representation of Complex functions.
Since Complex Functions can be seen as R^2 -> R^2 we cannot represent those functions w/o information loss in a < 4 dimenesional setting.
However complex analysis is only interested in certain behaviors of complex funtions in which mulitplication with a real Value c does not change anything.
Hence E. Wegert suggest this way of plotting complex functions in https://www.springer.com/de/book/9783034801799#otherversion=9783034801805.
The color only depends on the argument phi of z = r * e ^ (i * phi) and is directly mapped to the corresponding hsv value (and some shades are added to add contrast).
To my own surprise I actually started to like this representation aswell.
E. Wegert however seems to have fallen in love and even sells calendars with such plots: https://tu-freiberg.de/en/fakult1/ana/institute/institute-of-applied-analysis/organisation/complex-beauties

Some of the Information which can be gathered from a plot:
Zeros of the function (negative rotation)
Poles of the function (positive rotation)
The degree of the Pole / Zero (Number of rotations)
Saddlepoints of the function
Even some coefficients of the power series can be calculated with little additional information.

The final image is the plot of a function f : C -> C, f(z) = (z - z1) * ... * (z - z31) * (1 / (z - w1)) * ... * (1 / (z - w31))
The Wettbewerb file used a strategy, where when poles and zeroes are close to each other, they surrounding will stay very smooth in comparison to functions,
where the pole and zero are far from each other. This technique is discussed in the calender entry of february 2011. (The calender can be found in ./wett/etc/cal11.pdf)
The final function has 31 zeros and the same amount of poles, where each pole is located very close to a zero value.
To create the Wettbewerb file, the command "saveComplexBeauty "FILENAME" lampda" has to be run (it will however overwrite files with the same name).
I created the file with high resolution (resX = 8192) but this took quite a while due to the high number of zero and pole values => high "degree" of the function.
Do test if the file was actually created with my code, you might want to reduce the resX to 400 in the saveComplexBeauty function.

How caotic (and narcotic) plots of functions where poles and zero values are further from each other can be seen in various plots in "./wett/etc".
While those plots are more intersting from a complex analysis point of few, they might not be considered as hübsch as the Abgabe file.
-}

{-WETT-}
type ComplexFunction = Complex Double -> Complex Double

hueFromComplex :: Complex Double -> Double
hueFromComplex z = fixDeg degPhi
    where   phi = pi + phase z  -- (0, 2pi]
            normPhi = phi / 2.0 / pi    -- (0, 1]
            degPhi = normPhi * 360.0 + 180.0    -- reverses the +pi in the first line
            fixDeg :: Double -> Double
            fixDeg d
                | d >= 360   = d - 360
                | otherwise = d

-- Calculates the rgb value from the hue and adds the shaded partitions
rgbFromHue :: Int -> Double -> RGB Double
rgbFromHue partitions h = hsv h 1 value                         -- Example partitions = 20
    where   partitionSize = 360 / fromIntegral partitions       -- 360 / 20             -> 18
            fh = floor h                                        -- 312.12               -> 312
            numPart = floor $ h / partitionSize                 -- 312.12 -> 312 / 18   = 17
            modValue = h - fromIntegral numPart * partitionSize -- 312.12 -> 17 * 18    = 6.12
            l = 0.7 -- value in [0.7, 1] as suggested by E. Wegert
            u = 1
            normValue = modValue / partitionSize
            value = l + (u - l) * normValue

-- Converts between the RGB Formats
pixelRGB8FromRGB :: RGB Double -> PixelRGB8
pixelRGB8FromRGB c = PixelRGB8 (doubleToWord8 r) (doubleToWord8 g) (doubleToWord8 b)
    where   r = channelRed c
            g = channelGreen c
            b = channelBlue c

-- (0, 1] -> {0, ..., 255}
doubleToWord8 :: Double -> Word8
doubleToWord8 d = round $ d * 255

-- Creates a Complex Beatuy with min realValue minA, min imValue minB, ...
-- resX specifies the pixel on the real axis, resY ...
-- f is the function which we plot
-- partitions specifies into how many partitions the circle should be divided for the "Wave" - effect
-- I have NO idea where, but somewhere I conjugate the value, so I had to add another conjugate in the function                                                v
createComplexBeauty :: (Double, Double) -> (Double, Double) -> Int -> Int -> ComplexFunction -> Int -> Image PixelRGB8
createComplexBeauty (minA, maxA) (minB, maxB) resX resY f partitions = generateImage (\x -> pixelRGB8FromRGB . rgbFromHue partitions . hueFromComplex . f . conjugate . complexFromPos x) resX resY
    where   diffA = maxA - minA
            stepA = diffA / fromIntegral resX
            diffB = maxB - minB
            stepB = diffB / fromIntegral resY
            -- indicies are 0 based
            complexFromPos :: Int -> Int -> Complex Double
            complexFromPos i j = (minA + fromIntegral i * stepA) :+ (minB + fromIntegral j * stepB)

saveComplexBeauty :: String -> ComplexFunction -> IO()
saveComplexBeauty functionName f = do
        g <- newStdGen
        let resX = 500 -- I used 8192
            (minA, maxA, minB, maxB) = (-1.55, 1.55, -1.15, 1.15)
            resY = floor $ fromIntegral resX * (maxB - minB) / (maxA - minA)
            img = ImageRGB8 $ createComplexBeauty (minA, maxA) (minB, maxB) resX resY f 20
            path = "./wett/" ++ functionName ++ ".png"
        savePngImage path img

-- Used for the Blaschke Function and other Plots in ./wett/etc
createOtherExamples :: IO()
createOtherExamples = do
        g <- newStdGen
        let rs = randoms g :: [Double]
            as = take 50 rs
            bs = take 50 $ drop 50 rs
            randomVals = zipWith (\a b -> ((-1) + 2*a) :+ ((-1) + 2*b)) as bs
            f = blaschkeFunction randomVals -- This is the Blaschke Function
            img = ImageRGB8 $ createComplexBeauty (-1.55, 1.55) (-1.55, 1.55) 1000 1000 (\z -> (z+0.5) * (z-(0.5 :+ 0.5))^2/(z-(0.5 :+ (-0.5)))^2) 20
            path = "./wett/etc/(z + 0.5)(z - (0.5 + 0.5i))^2(z - (0.5 - 0.5i))^-2.png"
        savePngImage path img
{-TTEW-}

{-Different usable functions-}
square :: ComplexFunction
square z = z^2

lampda :: ComplexFunction
lampda z =  createSmallOffsetEffect straight offsetStraight z *     -- Long bar from the Lampda
            createSmallOffsetEffect small offsetSmall z *           -- Short bar from the Lampda
            createSmallOffsetEffect (map (+ offSetToLeft) small) (map (+ offSetToLeft) offsetSmall) z * -- Lower part of >
            product [z - conjugate(w + offSetToLeft) | w <- drop 1 small] * product [1/(z - conjugate(w + offSetToLeft)) | w <- drop 1 offsetSmall] *   -- Upper part of >
            createSmallOffsetEffect upper offsetUpper z *           -- Upper equal line
            createSmallOffsetEffect lower offsetLower z             -- Lower equal line
            
    where   -- number of points on long bar
            nStr = 10
            (straight, offsetStraight) = line nStr offset ((-0.80) :+ 1) (0.80 :+ (-1))

            -- number of points on short bars (lampda + >)
            nSmall = 5
            (small, offsetSmall) = line nSmall offset 0 ((-0.80) :+ (-1))
            -- Use this in lower > to get the matching orientation effect (replace small with smallReverse in "createSmallOffsetEffect (map (+ offSetToLeft) small) (...) z")
            (smallReverse, offsetSmallReversed) = line nSmall (-offset) 0 ((-0.80) :+ (-1))
            -- Offset to the left for the >
            offSetToLeft = (-0.5) :+ 0

            -- number of points of upper / lower bar in equal sign
            nUpper = 4
            nLower = 3
            (upper, offsetUpper) = line nUpper offset (0 :+ 0.3) (1.4 :+ 0.3)
            (lower, offsetLower) = line nLower offset (0.4 :+ (-0.1)) (1.4 :+ (-0.1))   -- For the different orienation: (-offset) (0.4 :+ (-0.1)) (1.4 :+ (-0.1))
            
            -- Offset between Zeroes and Poles
            offset = 0.2

--      numberOfDots -> offSet -> Start -> End -> (normalLine, offsetLine)
line :: Int -> Double -> Complex Double -> Complex Double -> ([Complex Double], [Complex Double])
line n offSet start end = (normalLine, offsetLine)
    where   way = end - start
            step = way / (fromIntegral (n-1) :+ 0)
            normalLine = [start + step * (fromIntegral t :+ 0) | t <- [0..n-1]]
            offsetLine = [start + step * ((fromIntegral t + offSet) :+ 0) | t <- [0..n-1]]

-- Creates a complex function with zero values zeroes and poles poles (and no other zeroes / poles)
createSmallOffsetEffect :: [Complex Double] -> [Complex Double] -> Complex Double -> Complex Double
createSmallOffsetEffect zeroes poles z = product [z - w | w <- zeroes] * product [1/(z - w) | w <- poles]

spiral :: ComplexFunction
spiral z = product [1/(z - w) | w <- offsetAlongCurve] * product [z - w | w <- alongCurve]

alongCurve :: [Complex Double]
alongCurve = map complCurve [0..20]

offsetAlongCurve :: [Complex Double]
offsetAlongCurve = map (complCurve . (+0.15)) [0..19]

-- [0, inf) -> C
-- Once circ = 40
complCurve :: Double -> Complex Double
complCurve t = mkPolar (3 + (sin (2 * nt) + cos (3 * nt)) / 2) nt
    where   nt = t * pi * 2 / 19

blaschkeFunction :: [Complex Double] -> ComplexFunction
blaschkeFunction vals z = product [(z - w) / (1 - conjugate w * z) | w <- vals]