Skip to content

Instantly share code, notes, and snippets.

@mniip
Last active November 28, 2020 00:23
Show Gist options
  • Save mniip/39c96123b24e771ec3087dd214106e78 to your computer and use it in GitHub Desktop.
Save mniip/39c96123b24e771ec3087dd214106e78 to your computer and use it in GitHub Desktop.
{-
This program generates a poly-line approximation to a De Rham curve. The curve
is specified using a collection of affine contractions. Under the assumption
that all fixed points are in the unit circle, the program is able to guarantee
an upper bound on the deviation of the actual curve from the poly-line
approximation.
The output is a series of lines with three columns each: the value of the
parameter in [0, 1], the real part, and the imaginary part of the curve.
Any real number can be specified as a fraction of two real numbers, e.g. 1/2,
and any complex number can be specified in either of the forms:
1,2
1+2i
1+i2
Usage: DeRham [--q | --qout] [--cts] --tolerance TOL [--epsilon EPS]
(--cesaro z | --koch-peano z | --generic [a:b:]c:d:e:f |
--arrowhead a:b:c | --affine a:b:c:d:e:f |
--triangle a:b:c:d:e:f)
Available options:
-h,--help Show this help text
--q Do intermediate arithmetic in Q. This is slow but
allows epsilon below 1e-14
--qout Also output the curve points as elements of Q, in a/b
format
--cts Don't insert empty lines in places where continuity
could not be ensured
--tolerance TOL Upper bound on how much the curve is allowed to
deviate from the poly-line
--epsilon EPS Precision with which points on the curve are
calculated. Also continuity will not be ensured for
parameterization steps smaller than EPS. (Default:
1e-14)
--cesaro z Use Cesaro affine maps corresponding to z, a complex
number, to generate Cesaro-Faber curves. Uses maps:
w <- zw
w <- z + (1-z)w
This requires |z| < 1, |1-z| < 1.
--koch-peano z Use Koch-Peano affine maps corresponding to z, a
complex number, to generate Peano and Koch curves.
Uses maps:
w <- zw*
w <- z + (1-z)w*
This requires |z| < 1, |1-z| < 1.
--generic [a:b:]c:d:e:f Use a generic pair of affine maps that have fixpoints
at 0 and 1:
[1] [ 1 0 0 ] [1]
[x] <- [ 0 a d ] [x]
[y] [ 0 b c ] [y]
[1] [ 1 0 0 ] [1]
[x] <- [ a 1-a e ] [x]
[y] [ b -b f ] [y]
where a, b, c, d, e, f are real numbers. This puts
the midpoint of the curve at a,b. The values of a and
b are 1/2,1 unless otherwise specified.
--arrowhead a:b:c Use three maps that generates the arrowhead curve by
mapping a triangle to three corners of itself. The
complex numbers a, b, and c choose how the sides of
the triangle are partitioned. Setting a=b=c=1/2
generates the standard Sierpinski triangle.
Incompatible with --q.
--affine a:b:c:d:e:f Specify an arbitrary affine map using real numbers a,
b, c, d, e, f:
[1] [ 1 0 0 ] [1]
[x] <- [ e a b ] [x]
[y] [ f c d ] [y]
Repeat this option multiple times to specify multiple
maps in this way.
--triangle a:b:c:d:e:f Specify an arbitrary affine map that would map the
triangle a,b,c to the triangle d,e,f; where a, b, c,
d, e, f are complex numbers. Repeat this option
multiple times to specify multiple maps in this way.
-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveFunctor #-}
import Control.Monad
import Data.List
import Data.Ratio
import System.IO
import Options.Applicative
import qualified Options.Applicative.Help as H
import Text.ParserCombinators.ReadP hiding (option, many)
import Numeric
-- | ''Data.Complex'' requires 'RealFrac' for 'Num'
data Complex a = a :+ a deriving (Eq, Show, Functor)
infix 6 :+
instance Num a => Num (Complex a) where
(a :+ b) + (c :+ d) = (a + c) :+ (b + d)
(a :+ b) - (c :+ d) = (a - c) :+ (b - d)
(a :+ b) * (c :+ d) = (a * c - b * d) :+ (a * d + b * c)
negate (a :+ b) = (negate a :+ negate b)
abs = error "abs"
signum = error "signum"
fromInteger n = (fromInteger n :+ 0)
cis :: Floating a => a -> Complex a
cis phi = cos phi :+ sin phi
magnitudeSq :: Num a => Complex a -> a
magnitudeSq (x :+ y) = x * x + y * y
-- | Affine map of the plane
data Affine a = Affine
{ dxdx :: !a
, dxdy :: !a
, dydx :: !a
, dydy :: !a
, x0 :: !a
, y0 :: !a
} deriving (Eq, Ord, Show, Functor)
applyAffine :: Num a => Affine a -> Complex a -> Complex a
applyAffine (Affine dxdp dxdq dydp dydq x0 y0) (dp :+ dq)
= (x0 + dxdp * dp + dxdq * dq) :+ (y0 + dydp * dp + dydq * dq)
-- | >>> affineCenter f = applyAffine f (0 :+ 0)
affineCenter :: Affine a -> Complex a
affineCenter (Affine _ _ _ _ x0 y0) = x0 :+ y0
-- | Frobenius norm squared
affineNormSq :: Num a => Affine a -> a
affineNormSq (Affine dxdx dxdy dydx dydy _ _)
= dxdx^2 + dxdy^2 + dydx^2 + dydy^2
-- | >>> applyAffine (f <> g) v = applyAffine f (applyAffine g v)
instance Num a => Semigroup (Affine a) where
(Affine dxdp dxdq dydp dydq x0 y0) <> (Affine dpdx dpdy dqdx dqdy p0 q0)
= Affine
(dxdp * dpdx + dxdq * dqdx)
(dxdp * dpdy + dxdq * dqdy)
(dydp * dpdx + dydq * dqdx)
(dydp * dpdy + dydq * dqdy)
(x0 + dxdp * p0 + dxdq * q0)
(y0 + dydp * p0 + dydq * q0)
instance Num a => Monoid (Affine a) where
mempty = Affine 1 0 0 1 0 0
-- | Assuming all fixpoints are in the unit circle around 0+0i, returns an affine map that turns the unit circle into a neighborhood of the curve whose radius is smaller than given tolerance
deRhamCurveNbhd :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> Affine a
deRhamCurveNbhd maps !tolerance = go mempty
where
n = length maps
go !a !t
| affineNormSq a < tolerance^2 = a
| i <- min (n - 1) $ floor (fromIntegral n * t)
= go (a <> maps !! i) (fromIntegral n * t - fromIntegral i)
-- | Assuming all fixpoints are in the unit circle around 0+0i, returns an affine map that turns the unit circle into a neighborhood corresponding to the given uncertainty in the argument
deRhamCurvePrecision :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> Affine a
deRhamCurvePrecision maps = go mempty
where
n = length maps
go !a !dt !t
| dt >= 1 = a
| i <- min (n - 1) $ floor (fromIntegral n * t)
= go (a <> maps !! i) (fromIntegral n * dt) (fromIntegral n * t - fromIntegral i)
tabulateCurve :: (Ord a, RealFrac a) => [Affine a] -> a -> a -> (a, a) -> [Maybe (a, Complex a)]
tabulateCurve maps !machineEps !tolerance (a, b) = let ea = eval a; eb = eval b in Just (a, ea) : go ea eb a b
where
eval !t = affineCenter $ deRhamCurveNbhd maps machineEps t
go !ea !eb !a !b
| magnitudeSq (ea - eb) < tolerance^2 && affineNormSq precision < tolerance^2 = [Just (mid, em), Just (b, eb)]
| abs (a - mid) < machineEps || abs (b - mid) < machineEps =
if magnitudeSq (ea - eb) < tolerance^2
then [Just (b, eb)]
else [Nothing, Just (b, eb)]
| otherwise = go ea em a mid ++ go em eb mid b
where
dt = (b - a) / 2
mid = (a + b) / 2
precision = deRhamCurvePrecision maps dt mid
em = eval mid
-- | >>> applyAffine (affineFromImages (input1, output1) (input2, output2) (input3, output3)) input_i = output_i
affineFromImages :: Fractional a => (Complex a, Complex a) -> (Complex a, Complex a) -> (Complex a, Complex a) -> Affine a
affineFromImages (x1 :+ y1, z1 :+ w1) (x2 :+ y2, z2 :+ w2) (x3 :+ y3, z3 :+ w3) =
Affine (a z1 z2 z3) (b z1 z2 z3) (a w1 w2 w3) (b w1 w2 w3) (e z1 z2 z3) (e w1 w2 w3)
where
det = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
a p1 p2 p3 = (p1 * (y2 - y3) + p2 * (y3 - y1) + p3 * (y1 - y2)) / det
b p1 p2 p3 = (p1 * (x3 - x2) + p2 * (x1 - x3) + p3 * (x2 - x1)) / det
e p1 p2 p3 = (p1 * (x2 * y3 - x3 * y2) + p2 * (x3 * y1 - x1 * y3) + p3 * (x1 * y2 - x2 * y1)) / det
cesaro :: Num a => Complex a -> [Affine a]
cesaro (a :+ b) =
[ Affine a (-b) b a 0 0
, Affine (1 - a) b (-b) (1 - a) a b
]
kochPeano :: Num a => Complex a -> [Affine a]
kochPeano (a :+ b) =
[ Affine a b b (-a) 0 0
, Affine (1 - a) (-b) (-b) (a - 1) a b
]
generic :: Num a => a -> a -> a -> a -> a -> a -> [Affine a]
generic alpha beta delta epsilon zeta eta =
[ Affine alpha delta beta epsilon 0 0
, Affine (1 - alpha) zeta (-beta) eta alpha beta
]
arrowhead :: Floating a => Complex a -> Complex a -> Complex a -> [Affine a]
arrowhead alpha beta gamma =
[ affineFromImages (a, a) (b, q) (c, r)
, affineFromImages (a, r) (b, b) (c, p)
, affineFromImages (a, p) (b, q) (c, c)
]
where
a = 0 :+ 0
b = cis (pi / 3)
c = 1 :+ 0
p = (1 - alpha) * b + alpha * c
q = (1 - beta) * c + beta * a
r = (1 - gamma) * a + gamma * b
data OptCurves
= Cesaro (Complex Rational)
| KochPeano (Complex Rational)
| Generic Rational Rational Rational Rational Rational Rational
| Arrowhead (Complex Rational) (Complex Rational) (Complex Rational)
| Maps [OptMap]
deriving Show
data OptMap
= AffineMap Rational Rational Rational Rational Rational Rational
| TriangleMap (Complex Rational, Complex Rational, Complex Rational) (Complex Rational, Complex Rational, Complex Rational)
deriving Show
parseCurve :: OptCurves -> Either [Affine Double] [Affine Rational]
parseCurve (Cesaro z) = Right $ cesaro z
parseCurve (KochPeano z) = Right $ kochPeano z
parseCurve (Generic a b c d e f) = Right $ generic a b c d e f
parseCurve (Arrowhead a b c) = Left $ arrowhead (fmap fromRational a) (fmap fromRational b) (fmap fromRational c)
parseCurve (Maps maps) = Right $ map parseMap maps
where
parseMap (AffineMap a b c d e f) = Affine a b c d e f
parseMap (TriangleMap (x1 :+ y1, x2 :+ y2, x3 :+ y3) (z1 :+ w1, z2 :+ w2, z3 :+ w3))
= affineFromImages (x1 :+ y1, z1 :+ w1) (x2 :+ y2, z2 :+ w2) (x3 :+ y3, z3 :+ w3)
isContraction :: (Ord a, Fractional a) => a -> Affine a -> Bool
isContraction q (Affine dxdx dxdy dydx dydy _ _) = if discr > 0
then abs (dxdx + dydy) < 2 * q && (2 * q - abs (dxdx + dydy))^2 > discr
else 4 * q^2 > discr + (dxdx + dydy)^2
where discr = (dxdx - dydy)^2 + 4 * dxdy * dydx
readP :: ReadP a -> ReadM a
readP p = eitherReader $ \str -> case readP_to_S p str of
(x, ""):_ -> return x
_ -> Left $ "no parse: " ++ show str
fraction :: ReadP Rational
fraction = liftA2 (/) (float <* skipSpaces <* char '/') (skipSpaces *> float) <|> float
where
float = readS_to_P readFloat
signed :: Num a => ReadP a -> ReadP a
signed c = (negate <$> (char '-' *> c)) <|> c
point :: Num a => ReadP a -> ReadP (Complex a)
point c = liftA2 (:+) (signed c <* skipSpaces <* char ',') (skipSpaces *> signed c) <|> expr
where
expr =
liftA2 (+) (sterm <* skipSpaces <* char '+') (skipSpaces *> term) <|>
liftA2 (-) (term <* skipSpaces <* char '-') (skipSpaces *> term) <|>
sterm
sterm = (negate <$> (char '-' *> skipSpaces *> term)) <|> term
term =
((0 :+) <$> (char 'i' *> skipSpaces *> c)) <|>
((0 :+ 1) <$ char 'i') <|>
((0 :+) <$> (c <* skipSpaces <* char 'i')) <|>
((:+ 0) <$> c)
(<:*:>) :: ReadP (a -> b) -> ReadP a -> ReadP b
c <:*:> d = (c <* skipSpaces <* char ':') <*> (skipSpaces *> d)
infixl 4 <:*:>
data OptComp = DoubleComp | QComp | QOut deriving (Eq, Show)
data Opts = Opts
{ optComp :: !OptComp
, optContinuous :: !Bool
, optTolerance :: !Rational
, optEpsilon :: !Rational
, optCurves :: OptCurves
} deriving Show
optsParser :: Parser Opts
optsParser = pure Opts
<*> (flag' QComp (long "q" <> help qHelp)
<|> flag DoubleComp QOut (long "qout" <> help qOutHelp))
<*> switch (long "cts" <> help ctsHelp)
<*> option fractionP (long "tolerance" <> metavar "TOL" <> help toleranceHelp)
<*> (option fractionP (long "epsilon" <> metavar "EPS" <> help epsilonHelp) <|> pure 1e-14)
<*> curvesParser
where
curvesParser =
(Cesaro <$> option pointP
(long "cesaro" <> metavar "z" <> helpDoc (Just cesaroHelp)))
<|> (KochPeano <$> option pointP
(long "koch-peano" <> metavar "z" <> helpDoc (Just kochPeanoHelp)))
<|> ((\(a, b, c, d, e, f) -> Generic a b c d e f) <$> option genericP
(long "generic" <> metavar "[a:b:]c:d:e:f" <> helpDoc (Just genericHelp)))
<|> ((\(a, b, c) -> Arrowhead a b c) <$> option arrowheadP
(long "arrowhead" <> metavar "a:b:c" <> helpDoc (Just arrowheadHelp)))
<|> (Maps <$> some mapParser)
fractionP = readP fraction
pointP = readP $ point fraction
genericP = readP $ g4 <|> g6
where
g4 = ((,,,,,) 0.5 1.0) <$> c <:*:> c <:*:> c <:*:> c
g6 = (,,,,,) <$> c <:*:> c <:*:> c <:*:> c <:*:> c <:*:> c
c = signed fraction
arrowheadP = readP $ (,,) <$> c <:*:> c <:*:> c
where c = point fraction
mapParser =
((\(a, b, c, d, e, f) -> AffineMap a b c d e f) <$> option affineP
(long "affine" <> metavar "a:b:c:d:e:f" <> helpDoc (Just affineHelp)))
<|> ((\(a, b) -> TriangleMap a b) <$> option trianglesP
(long "triangle" <> metavar "a:b:c:d:e:f" <> helpDoc (Just triangleHelp)))
affineP = readP $ (,,,,,) <$> c <:*:> c <:*:> c <:*:> c <:*:> c <:*:> c
where c = signed fraction
trianglesP = readP $ (,) <$> ((,,) <$> c <:*:> c <:*:> c) <:*:> ((,,) <$> c <:*:> c <:*:> c)
where c = point fraction
qHelp = "Do intermediate arithmetic in Q. This is slow but allows epsilon below 1e-14"
qOutHelp = "Also output the curve points as elements of Q, in a/b format"
ctsHelp = "Don't insert empty lines in places where continuity could not be ensured"
toleranceHelp = "Upper bound on how much the curve is allowed to deviate from the poly-line"
epsilonHelp = "\
\Precision with which points on the curve are calculated. Also continuity \
\will not be ensured for parameterization steps smaller than EPS. \
\(Default: 1e-14)"
cesaroHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Use Cesaro affine maps corresponding to z, a complex number, to \
\generate Cesaro-Faber curves. Uses maps:"
, H.indent 4 (H.vsep [H.text "w <- zw", H.text "w <- z + (1-z)w"])
, H.extractChunk $ H.paragraph "\
\This requires |z| < 1, |1-z| < 1."
]
kochPeanoHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Use Koch-Peano affine maps corresponding to z, a complex number, to \
\generate Peano and Koch curves. Uses maps:"
, H.indent 4 (H.vsep [H.text "w <- zw*", H.text "w <- z + (1-z)w*"])
, H.extractChunk $ H.paragraph "\
\This requires |z| < 1, |1-z| < 1."
]
matrix xss = H.hsep $ map (H.vsep . map H.text)
([replicate height "["] ++ xss ++ [replicate height "]"])
where
height = case xss of
[] -> 1
(xs:_) -> length xs
genericHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Use a generic pair of affine maps that have fixpoints at 0 and 1:"
, H.indent 4 $ H.vsep
[ H.text "[1] [ 1 0 0 ] [1]"
, H.text "[x] <- [ 0 a d ] [x]"
, H.text "[y] [ 0 b c ] [y]"
, H.empty
, H.text "[1] [ 1 0 0 ] [1]"
, H.text "[x] <- [ a 1-a e ] [x]"
, H.text "[y] [ b -b f ] [y]"
]
, H.extractChunk $ H.paragraph "\
\where a, b, c, d, e, f are real numbers. This puts the midpoint of the \
\curve at a,b. The values of a and b are 1/2,1 unless otherwise \
\specified."
]
arrowheadHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Use three maps that generates the arrowhead curve by mapping a triangle \
\to three corners of itself. The complex numbers a, b, and c choose how \
\the sides of the triangle are partitioned. Setting a=b=c=1/2 generates \
\the standard Sierpinski triangle. Incompatible with --q."
]
affineHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Specify an arbitrary affine map using real numbers a, b, c, d, e, f:"
, H.indent 4 $ H.vsep
[ H.text "[1] [ 1 0 0 ] [1]"
, H.text "[x] <- [ e a b ] [x]"
, H.text "[y] [ f c d ] [y]"
]
, H.extractChunk $ H.paragraph "\
\Repeat this option multiple times to specify multiple maps in this way."
]
triangleHelp = H.vsep
[ H.extractChunk $ H.paragraph "\
\Specify an arbitrary affine map that would map the triangle a,b,c to \
\the triangle d,e,f; where a, b, c, d, e, f are complex numbers. Repeat \
\this option multiple times to specify multiple maps in this way."
]
optsInfo :: ParserInfo Opts
optsInfo = info (helper <*> optsParser) $
headerDoc (Just $ H.vsep
[ H.extractChunk $ H.paragraph "\
\This program generates a poly-line approximation to a De Rham curve. The \
\curve is specified using a collection of affine contractions. Under the \
\assumption that all fixed points are in the unit circle, the program is \
\able to guarantee an upper bound on the deviation of the actual curve \
\from the poly-line approximation."
, H.empty
, H.extractChunk $ H.paragraph "\
\The output is a series of lines with three columns each: the value of the \
\parameter in [0, 1], the real part, and the imaginary part of the curve."
, H.empty
, H.extractChunk (H.paragraph "\
\Any real number can be specified as a fraction of two real numbers, e.g. \
\1/2, and any complex number can be specified in either of the forms:")
, H.indent 4 (H.vsep [H.text "1,2", H.text "1+2i", H.text "1+i2"])
])
main :: IO ()
main = do
opts <- execParser optsInfo
let maps = parseCurve (optCurves opts)
case optComp opts of
DoubleComp -> go opts id $ either id (map $ fmap fromRational) maps
_ -> case maps of
Left maps -> do
hPutStrLn stderr "Not rational maps: "
hPutStrLn stderr $ show maps
hPutStrLn stderr "Try:"
hPutStrLn stderr $ formatMaps maps
Right maps -> go opts fromRational maps
where
go :: (RealFrac a, Show a) => Opts -> (a -> Double) -> [Affine a] -> IO ()
go opts fmt maps = do
forM_ (filter (not . isContraction 0.999) maps) $ \a -> do
hPutStrLn stderr $ "Not a contraction (q > 0.999): " ++ show a
forM_ (tabulateCurve maps (fromRational $ optEpsilon opts)
(fromRational $ optTolerance opts) (0, 1)) $ \m -> case m of
Just (t, x :+ y) -> do
if optComp opts == QOut
then putStrLn $ formatQ t ++ " " ++ formatQ x ++ " " ++ formatQ y
else putStrLn $ show (fmt t) ++ " " ++ show (fmt x) ++ " " ++ show (fmt y)
Nothing -> unless (optContinuous opts) $ do
putStrLn ""
formatQ x = show (numerator q) ++ "/" ++ show (denominator q)
where q = toRational x
formatMaps maps = intercalate " " $ map ("--affine " ++) $ map mkMap maps
mkMap (Affine a b c d e f) = intercalate ":" $ map show [a,b,c,d,e,f]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment