Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
917 lines (749 sloc) 38.7 KB
-- Copyright 2012 Justus Sagemüller.
--
-- This file is part of the Cqtx library.
-- This library is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
-- This library is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-- GNU General Public License for more details.
-- You should have received a copy of the GNU General Public License
-- along with this library. If not, see <http://www.gnu.org/licenses/>.
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE ImpredicativeTypes #-}
{-# LANGUAGE MultiParamTypeClasses #-}
-- {-# LANGUAGE FunctionalDependencies#-}
{-# LANGUAGE FlexibleInstances #-}
-- {-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- |
-- Module : CodeGen.CXX.Cqtx.PhqFn
-- Copyright : (c) Justus Sagemüller 2012
-- License : GPL v3
--
-- Maintainer : (@) sagemuej $ smail.uni-koeln.de
-- Stability : experimental
-- Portability : portable
--
-- While the cqtx library provides quite versatile and /in principle/ very
-- low-requirement fitting of general physical functions (namely without the
-- need for any gradients etc.), in practise this is hampered by the cumbersomeness
-- of defining the fittable_phmsqfn subclasses. The macros in @fitfnmacros.h@ can
-- only partially alleviate this: due to the limits of the C preprocessor, most of
-- the calculation load needs to be done in also cumbersome C++ templates or,
-- more likely, at runtime. At the moment, this is implemented rather poorly, in
-- particular the almost bogosort-quality dimensional analysis.
--
-- This Haskell module does not claim to be an optimal solution, but it produces
-- significantly better code than the CPP macros, with very safe and easy invocation.
-- In the future this might become a lot more powerful, since it is also possible
-- to automatically create highly optimised specialised versions of the functions, e.g.
-- squaredistance-calculation on CUDA.
--
module CodeGen.CXX.Cqtx.PhqFn( phqFn, phqFlatMultiIdFn, phqSingleIdMultiFn
-- * General Cqtx code generation
, module CodeGen.CXX.Code
, CqtxCode
, CqtxConfig, withDefaultCqtxConfig
, PhqfnDefining(..)
-- * Length-indexed lists
, IsolenEnd(P), IsolenCons(..)
) where
import CodeGen.CXX.Code
import Control.Monad hiding (forM_)
import Control.Monad.Writer hiding (forM_)
import Control.Monad.Reader hiding (forM_)
import Control.Arrow
import Data.Function
import Data.List (sort, sortBy, intersperse, intercalate, group, find)
import Data.Monoid
import Data.Ratio
import Data.Maybe
import Data.Tuple
import Data.Foldable
import Prelude hiding (foldr, concat, sum, product, any, elem)
-- import Data.Hashable
-- import Data.HashMap
-- | Create a physical function that can be fitted to measured data using the cqtx
-- algorithms such as @evolution_minimizer@. The invocation is similar to the
-- primitive CPP macros, but type-safe and allows full Haskell syntax in the
-- definition – though this can of course be exploited only so much, since phqfns
-- are very limited in their abilities. To make sure these limits are maintained
-- we use universally-quantised arguments (which also suits the implementation very
-- well).
--
-- For example, the standard gaussian peak
-- @&#x1d434;&#x22c5;exp(-(&#x1d465;-&#x1d465;&#x2080;)&#xb2;/(2&#x22c5;&#x1d70e;&#xb2;))@
-- could be defined thus:
--
-- > phqFn "gaussPeak" ("x":."x_0":."\\sigma":."A":.P)
-- > (\ ( x :. x0 :. sigma :. a :.P)
-- > -> let q = (x-x0)/sigma in a * exp(-0.5 * q^2) )
--
-- The use of the type-determined–length lists makes it impossible to accidentally
-- give different numbers of parameter bindings and -labels.
--
-- Avoidance of duplicate calculation, as well as rudimentary optimisation, is taken
-- care for by this preprocessor.
phqFn :: forall paramsList . IsolenList paramsList
=> String -- ^ Name of the resulting phqFn C++ object
-> paramsList String -- ^ Default labels of the parameters
-> (forall x. PhqfnDefining x
=> paramsList x -> x) -- ^ Function definition, as a lambda
-> CqtxCode() -- ^ C++ class and object code for a cqtx fittable physical function corresponding to the given definition.
phqFn fnName sclLabels function
= phqFlatMultiIdFn fnName sclLabels P $ \P -> (P, \scl P -> function scl)
phqSingleIdMultiFn :: forall scalarPrmsList indexedPrmsList
. (IsolenList scalarPrmsList, IsolenList indexedPrmsList)
=> String
-> scalarPrmsList String
-> (String, Maybe Int)
-> ( PhqVarIndexer ->
( indexedPrmsList String
, forall x. PhqfnDefining x
=> scalarPrmsList x -> indexedPrmsList (IdxablePhqDefVar x) -> x) )
-> CqtxCode()
phqSingleIdMultiFn fnName' sclLabels indexer indexedFn
= phqFlatMultiIdFn fnName' sclLabels (indexer:.P) $
\(i:.P) -> first (fmap (,i)) $ indexedFn i
phqFlatMultiIdFn :: forall scalarPrmsList indexerList indexedPrmsList
. (IsolenList scalarPrmsList, IsolenList indexerList, IsolenList indexedPrmsList)
=> String
-> scalarPrmsList String
-> indexerList (String, Maybe Int)
-> ( indexerList PhqVarIndexer ->
( indexedPrmsList (String, PhqVarIndexer)
, forall x. PhqfnDefining x
=> scalarPrmsList x -> indexedPrmsList (IdxablePhqDefVar x) -> x) )
-> CqtxCode()
phqFlatMultiIdFn fnName' sclLabels ixerLabels indexedFn = ReaderT $ codeWith where
codeWith config = do
cxxLine $ " COPYABLE_PDERIVED_CLASS(/*"
cxxLine $ "class*/"++className++",/*: public */fittable_phmsqfn) {"
helpers <- cxxIndent 2 $ do
rangesDecl
offsetMgr
dimFetchDeclares
cxxLine $ " public:"
cxxIndent 2 $ do constructor
paramExamples helpers
functionEval
cxxLine $ "} " ++ (guard (null rangeConstsNeeded) >> fnName) ++ ";"
where
-- function :: forall x . PhqfnDefining x
-- => scalarPrmsList x -> indexedPrmsList (IdxablePhqDefVar x) -> x
ixaParams :: indexedPrmsList (String, PhqVarIndexer)
(ixaParams, function) = indexedFn indexers
fnResultTerm :: PhqFuncTerm
fnResultTerm = function (fmap (\i->PhqFnParameter $ show i) scalarParamIds)
(perfectZipWith (\vni (vnm, PhqVarIndexer ix inm)
-> IdxablePhqDefVar $ q vni vnm ix inm)
offsetConstsNeeded ixaParams)
where q vni vnm ix inm (PhqVarIndexer ix' inm')
| ix'==ix = PhqFnParameter $ vni ++ " + "
++ indizesPrefix ++ makeSafeCXXIdentifier inm
| otherwise = error
$ "In HsMacro-defined phqfunction '"++fnName++"':\n\
\ Using wrong indexer '"++inm'++"' (not adapted range!) \
\for indexing variable '"++vnm++"'.\n\
\ Correct indexer would be '"++inm++"'."
fnDimTrace :: DimTracer
fnDimTrace = function (fmap (VardimVar . PhqIdf) scalarParamIds)
(fmap (\i -> IdxablePhqDefVar
$ \_ -> VardimVar $ PhqIdf i) ixableParamIds)
functionEval :: CXXCode()
functionEval = do
cxxLine $ "auto operator()(const measure& parameters)const -> physquantity {"
cxxIndent 2 $ do
-- cxxLine $ "std::cout << \"evaluate "++fnName++" with parameters\\n\""
-- cxxLine $ " << parameters << std::endl;"
result <- cqtxTermImplementation "parameters" fnResultTerm
cxxSurround "return "result";"
cxxLine $ "}"
paramExamples :: [(Int,CXXExpression)] -> CXXCode()
paramExamples helpers = do
cxxLine $ "auto example_parameterset( const measure& constraints\n\
\ , const physquantity& desiredret)const override -> measure {"
cxxIndent 2 $ do
cxxLine $ "measure example;"
forM_ helpers $ \(n,helper) -> if n < nScalarParams
then do
cxxLine $ "if(!constraints.has(*argdrfs["++show n++"]))"
cxxLine $ " example.let(*argdrfs["++show n++"]) = "
++helper++"(constraints, desiredret);"
else do
let n' = n - nScalarParams
(idxerId, PhqVarIndexer _ i) = ixaParams !!@ n'
idxerRange = ixableParamRanges !!@ n'
offsetQ = offsetConstsNeeded !!@ n'
cxxLine $ "for(unsigned "++i++"=0\
\; "++i++"<"++idxerRange
++"; ++"++i++") {"
cxxIndent 2 $ do
let locatei = "argdrfs["++offsetQ++" + "++i++"]"
cxxLine $ "if(!constraints.has(*"++locatei++")) {"
cxxIndent 2 $ do
cxxLine $ "example.let(*"++locatei++")\
\ = "++helper++"(constraints, desiredret);"
cxxLine $ "}"
cxxLine $ "}"
cxxLine $ "return example;"
cxxLine $ "}"
dimFetchDeclares :: CXXCode [(Int, CXXExpression)]
dimFetchDeclares = forM ( map (False,) (toList scalarParamIds)
++ map (True, ) (toList ixableParamIds) )
$ \(needsIndexer, prmId) -> do
let functionName = "example_"++(guard needsIndexer >> "multi")
++"parameter"++show prmId++"value"
let resultOptions = relevantDimExprsFor (PhqIdf prmId) fnDimTrace
cxxLine $ "auto "++functionName++"( const measure& constraints\n\
\ , const physquantity& desiredret)const -> physquantity {"
cxxIndent 2 $ do
forM_ resultOptions $ \decomp -> do
cxxLine $ "{"
cxxIndent 2 $ do
let refMultiParamVarb n = "reference_for_multiparam_"++show n
cxxLine $ "bool thisdecomp_failed = false;"
forM_ (fixValDecomposition decomp) $ \(e,_) -> case e of
PhqFnParamVal (PhqIdf n) -> if n < nScalarParams
then do
cxxLine $ "if(!constraints.has(*argdrfs["++show n++"]))"
cxxLine $ " thisdecomp_failed = true;"
else do
let n' = n - nScalarParams
(idxerId, PhqVarIndexer _ i) = ixaParams !!@ n'
idxerRange = ixableParamRanges !!@ n'
offsetQ = offsetConstsNeeded !!@ n'
matchFound = "forany_"++i++"_pfound"
cxxLine $ "physquantity "++refMultiParamVarb n'++";"
cxxLine $ "{"
cxxIndent 2 $ do
cxxLine $ "bool "++matchFound++" = false;"
cxxLine $ "for(unsigned "++i++"=0\
\; "++i++"<"++idxerRange
++"; ++"++i++") {"
cxxIndent 2 $ do
let locatei = "argdrfs["++offsetQ++" + "++i++"]"
cxxLine $ "if(constraints.has(*"++locatei++")) {"
cxxIndent 2 $ do
cxxLine $ matchFound++" = true;"
cxxLine $ refMultiParamVarb n'
++" = "++locatei++"(constraints);"
cxxLine $ "break;"
cxxLine $ "}"
cxxLine $ "}"
cxxLine $ "if(!"++matchFound++") thisdecomp_failed = true;"
cxxLine $ "}"
_ -> return ()
cxxLine $ "if(!thisdecomp_failed) {"
cxxIndent 2 $ do
let q (PhqIdf n)
| n < nScalarParams = PhqFnParameter $ show n
| n' <-n-nScalarParams = PhqFnTempRef minBound $ refMultiParamVarb n'
result <- cqtxTermImplementation "constraints" $
calculateDimExpression q decomp
cxxSurround "physquantity result = " result ";"
cxxLine $ "return result.werror(result);"
cxxLine $ "}"
cxxLine $ "}"
cxxLine $ "std::cerr << \"Insufficient constraints given to determine the physical dimension\\n\
\ of parameter \\\"\" << *argdrfs["++show prmId++"] << \"\\\"\
\ in phqfn '"++fnName++"'.\\n Sufficient constraint choices would be:\\n\";"
forM_ resultOptions $ \(DimExpression decomp) -> do
cxxLine $ "std::cerr<<\" \"" ++ concat (intersperse"<<','"$
[ case fixv of
PhqDimlessConst -> ""
PhqFnParamVal (PhqIdf j) -> "<<miniTeX(*argdrfs["++show j++"])"
PhqFnResultVal -> "<<\"<function result>\""
| (fixv,_) <- decomp ]) ++ " << std::endl;"
cxxLine $ "abort();"
cxxLine $ "}"
return (prmId,functionName)
rangesDecl :: CXXCode()
rangesDecl = do
when (not $ null rangeConstsNeeded) .
cxxLine $ "unsigned "
++ intercalate ", " rangeConstsNeeded ++ ";"
when (not . null $ toList offsetConstsNeeded) .
cxxLine $ "unsigned "
++ intercalate ", " (toList offsetConstsNeeded) ++ ";"
rangeConstsNeeded :: [CXXExpression]
rangeConstsNeeded = map ( (indizesPrefix++) . (++indexRangePostfix)
. makeSafeCXXIdentifier . fst )
. filter ( isNothing . snd )
$ toList ixerLabels
offsetConstsNeeded :: indexedPrmsList CXXExpression
offsetConstsNeeded = fmap ( (ixablePPrefix++) . (++ixaOffsetPostfix)
. makeSafeCXXIdentifier . fst )
ixaParams
className
| null rangeConstsNeeded = fnName++"Function"
| otherwise = fnName
fnName = makeSafeCXXIdentifier fnName'
offsetMgr :: CXXCode()
offsetMgr = do
cxxLine $ "void manage_paramarr_offsets() {"
cxxIndent 2 $ do
cxxLine $ "unsigned stackp = "++show(isoLength scalarParamIds)++";"
forM_ (perfectZip offsetConstsNeeded ixableParamRanges)
$ \(osc, rng) -> do
cxxLine $ osc ++ " = stackp;"
cxxLine $ "stackp += " ++ rng ++ ";"
cxxLine $ "argdrfs.resize(stackp);"
cxxLine $ "}"
constructor :: CXXCode()
constructor = do
cxxLine $ className ++"("++args++")"++initialisers++" {"
cxxIndent 2 $ do
cxxLine $ "manage_paramarr_offsets();"
forM_ idxedDefaultLabels $ \(n,label) ->
cxxLine $ "argdrfs["++show n++"] = "++show label++";"
forM_ (perfectZip3 ixaParams ixableParamRanges offsetConstsNeeded)
$ \((ixaLbl, (PhqVarIndexer _ ixn)), rangeq, offsetQ) -> do
cxxLine $ "for (unsigned "++ixn++"=0\
\; "++ixn++" < "++rangeq
++"; ++"++ixn++")"
cxxIndent 2 $ do
cxxLine $ "argdrfs["++offsetQ++" + "++ixn++"] = "
++show ixaLbl++" + LaTeX_subscript("++ixn++");"
cxxLine $ "}"
where cstrArgs = rangeConstsNeeded
args = intercalate ", " $ map ("unsigned "++) cstrArgs
initialisers
| null $ toList ixableParamRanges = ""
| otherwise = "\n : "
++ intercalate ", " ( map (\a -> a++"("++a++")" ) cstrArgs )
indexers :: indexerList PhqVarIndexer
indexers = perfectZipWith PhqVarIndexer (enumFrom' 0) $ fmap fst ixerLabels
scalarParamIds :: scalarPrmsList Int
scalarParamIds = enumFrom' 0
ixableParamIds :: indexedPrmsList Int
ixableParamIds = enumFrom' nScalarParams
ixableParamRanges :: indexedPrmsList CXXExpression
ixableParamRanges = fmap (rngFind . snd) ixaParams
where rngFind (PhqVarIndexer ix _)
= case ixerLabels !!@ ix of
(_, Just n) -> show n
(name, _ ) -> indizesPrefix ++ makeSafeCXXIdentifier name
++ indexRangePostfix
idxedDefaultLabels = perfectZip scalarParamIds sclLabels
nScalarParams = isoLength scalarParamIds
indizesPrefix, ixablePPrefix, indexRangePostfix, ixaOffsetPostfix :: CXXExpression
indizesPrefix = "paramindex_"
ixablePPrefix = "ixaparam_"
indexRangePostfix = "_range"
ixaOffsetPostfix = "_offset"
data PhqIdf = PhqIdf Int deriving (Eq, Ord, Show)
argderefv :: CXXExpression -> CXXExpression -> String
argderefv paramSource n = "argdrfs["++n++"]("++paramSource++")"
data DimTracer = DimlessConstant (Maybe Rational)
| VardimVar PhqIdf
| DimEqualities [(DimTracer,DimTracer)] -- pairs of expressions that should have equal phys-dimension, but not necessarily related to the result
[DimTracer] -- expressions that should have the same phys-dimension as the result
| DimtraceProduct [DimTracer]
| DimtracePower DimTracer Rational
deriving(Show)
unknownDimlessConst :: DimTracer
unknownDimlessConst = DimlessConstant Nothing
beDimLess :: DimTracer -> DimTracer
beDimLess a = DimEqualities [(a, unknownDimlessConst)] [unknownDimlessConst]
biBeDimLess :: DimTracer -> DimTracer -> DimTracer
biBeDimLess a b = DimEqualities (map(,unknownDimlessConst)[a,b]) [unknownDimlessConst]
instance Num DimTracer where
fromInteger = DimlessConstant . Just . fromInteger
a + b = DimEqualities [(a,b)] [a,b]
a - b = DimEqualities [(a,b)] [a,b]
DimtraceProduct l * tr = DimtraceProduct (tr:l)
tr * DimtraceProduct l = DimtraceProduct (tr:l)
a * b = DimtraceProduct [a, b]
negate = id
abs = id
signum a = DimEqualities [(a,a)] [unknownDimlessConst]
instance Fractional DimTracer where
fromRational = DimlessConstant . Just
DimtraceProduct l / tr = DimtraceProduct (recip tr : l)
a / b = DimtraceProduct [a, recip b]
recip a = DimtracePower a (-1)
instance Floating DimTracer where
pi = unknownDimlessConst
exp = beDimLess; log = beDimLess
sqrt a = DimtracePower a (1/2)
a**(DimlessConstant (Just b)) = DimtracePower a b
a**b = biBeDimLess a b
logBase = biBeDimLess
sin = beDimLess; cos = beDimLess; tan = beDimLess
asin = beDimLess; acos = beDimLess; atan = beDimLess
sinh = beDimLess; cosh = beDimLess; tanh = beDimLess
asinh = beDimLess; acosh = beDimLess; atanh = beDimLess
data PhqFixValue = PhqDimlessConst
| PhqFnParamVal PhqIdf
| PhqFnResultVal
deriving(Eq, Ord, Show)
newtype DimExpression = DimExpression { fixValDecomposition :: [(PhqFixValue, Rational)] }
deriving(Eq, Ord, Show)
expExpMap :: (Rational -> Rational) -> DimExpression -> DimExpression
expExpMap f (DimExpression l) = DimExpression[ (a,f r) | (a,r)<-l ]
invDimExp :: DimExpression -> DimExpression
invDimExp = expExpMap negate
primDimExpr :: PhqFixValue -> DimExpression
primDimExpr v = DimExpression[(v,1)]
dimExpNormalForm :: DimExpression -> DimExpression -- DimExpression must always be normalised to this form
dimExpNormalForm (DimExpression l) = DimExpression . reduce . sortf $ l
where sortf = sortBy (compare`on`fst)
reduce ((PhqDimlessConst,_):l') = reduce l'
reduce ((a,r):β@(b,s):l')
| a==b = reduce $ (a,r+s):l'
| otherwise = (a,r) : reduce (β:l')
reduce l' = l'
instance Monoid DimExpression where
mempty = DimExpression[]
mappend (DimExpression a) (DimExpression b) = dimExpNormalForm $ DimExpression(a++b)
mconcat l = dimExpNormalForm . DimExpression $ l >>= fixValDecomposition
dimExprnComplexity :: DimExpression -> Integer
dimExprnComplexity (DimExpression l) = sum $ map complexity l
where complexity (_, fr) = abs(denominator fr) + abs(numerator fr - 1)
compareDimsBasis :: DimExpression -> DimExpression -> Ordering
compareDimsBasis (DimExpression l) (DimExpression r)
= comp (map fst l) (map fst r)
where comp [] [] = EQ
comp _ [] = GT
comp [] _ = LT
comp (l:ls) (r:rs)
| l<r = GT
| l>r = LT
| l==r = comp ls rs
strictlySimplerDimsBasis :: DimExpression -> DimExpression -> Bool
strictlySimplerDimsBasis (DimExpression l) (DimExpression r)
= comp (map fst l) (map fst r)
where comp _ [] = False
comp [] _ = True
comp (l:ls) (r:rs)
| l<r = False
| l>r = comp (l:ls) rs
| l==r = comp ls rs
traceAsValue :: DimTracer -> [DimExpression]
traceAsValue (DimlessConstant _) = return mempty
traceAsValue (VardimVar a) = return (primDimExpr $ PhqFnParamVal a)
traceAsValue (DimEqualities _ v) = v >>= traceAsValue
traceAsValue (DimtraceProduct l) = map mconcat . sequence $ map traceAsValue l
traceAsValue (DimtracePower a q) = map (expExpMap (q*)) $ traceAsValue a
(//-) :: DimExpression -> DimExpression -> DimExpression
a //- b = a <> invDimExp b
extractFixValExp :: PhqFixValue -> DimExpression -> (Rational, DimExpression)
extractFixValExp v = (\(a,b)->(a,DimExpression b)) . go id . fixValDecomposition
where go acc (t@(w,r):l)
| v==w = (r, acc l)
| v>w = go (acc.(t:)) l
go acc e = (0, acc e)
nubDimExprs :: [DimExpression] -> [DimExpression]
nubDimExprs = map head . group . sort
dimExpressionsFor :: PhqIdf -> DimTracer -> [DimExpression]
dimExpressionsFor idf e = go e $ primDimExpr PhqFnResultVal
where go :: DimTracer -> DimExpression -> [DimExpression]
go (DimEqualities mutualEqs resEqs) rev
= nubDimExprs $ ((`go`rev) =<< resEqs) ++ ( do
(ida,idb) <- mutualEqs
(go idb =<< traceAsValue ida)
++ (go ida =<< traceAsValue idb) )
go e@(DimtraceProduct l) rev = nubDimExprs $ do
way <- l
value <- traceAsValue way
invprod <- map invDimExp $ traceAsValue e
go way $ rev<>invprod<>value
go (DimtracePower p q) rev = go p $ expExpMap(/q) rev
go rest rev = do
value <- traceAsValue rest
let (q,p) = extractFixValExp (PhqFnParamVal idf) $ value //- rev
guard (q /= 0)
[expExpMap(/(-q))p]
relevantDimExprsFor :: PhqIdf -> DimTracer -> [DimExpression]
relevantDimExprsFor idf = prune . cplxtySort . dimExpressionsFor idf
where cplxtySort = sortBy (compare `on` negate . dimExprnComplexity)
prune [] = []
prune (e:es)
| any(`strictlySimplerDimsBasis`e)es = prune es
| otherwise = e : prune es
newtype CXXFunc = CXXFunc { wrapCXXFunc :: CXXExpression -> CXXExpression }
newtype CXXInfix = CXXInfix { wrapCXXInfix :: CXXExpression -> CXXExpression -> CXXExpression }
instance Eq CXXFunc where
CXXFunc f==CXXFunc g = f""==g""
instance Eq CXXInfix where
CXXInfix f==CXXInfix g = ""`f`""==""`g`""
cxxFunction :: CXXExpression -> CXXFunc
cxxFunction s = CXXFunc $ \e -> s++"("++e++")"
cxxInfix :: CXXExpression -> CXXInfix
cxxInfix s = CXXInfix $ \e1 e2 -> "("++e1++s++e2++")"
type PhysicalCqtxConst = String
data PhqFuncTerm = PhqFnDimlessConst Double
| PhqFnPhysicalConst PhysicalCqtxConst
| PhqFnTempRef Int CXXExpression
| PhqFnParameter
CXXExpression -- not the parameter itself, but an expression for its /index/ in the argdrfs array.
| PhqFnFuncApply CXXFunc PhqFuncTerm
| PhqFnInfixApply CXXInfix PhqFuncTerm PhqFuncTerm
| PhqFnInfixFoldOverIndexer PhqVarIndexer
CXXInfix
PhqFuncTerm -- init
PhqFuncTerm -- summand
deriving (Eq)
cxxFuncPhqApply :: CXXExpression -> PhqFuncTerm -> PhqFuncTerm
cxxFuncPhqApply s = PhqFnFuncApply $ cxxFunction s
cxxInfixPhqApply :: CXXExpression -> PhqFuncTerm -> PhqFuncTerm -> PhqFuncTerm
cxxInfixPhqApply s = PhqFnInfixApply $ cxxInfix s
isPrimitive :: PhqFuncTerm -> Bool
isPrimitive (PhqFnDimlessConst _) = True
isPrimitive (PhqFnPhysicalConst _) = True
isPrimitive (PhqFnParameter _) = True
isPrimitive _ = False
instance Num PhqFuncTerm where
fromInteger = PhqFnDimlessConst . fromInteger
(+) = cxxInfixPhqApply"+"
(-) = cxxInfixPhqApply"-"
a*(PhqFnDimlessConst 1) = a
(PhqFnDimlessConst 1)*b = b
a*b
| a==b = PhqFnFuncApply (CXXFunc $ \e -> "("++e++").squared()") a
| otherwise = cxxInfixPhqApply"*" a b
negate = cxxFuncPhqApply"-"
abs = cxxFuncPhqApply"abs"
signum = cxxFuncPhqApply"sgn"
instance Fractional PhqFuncTerm where
fromRational = PhqFnDimlessConst . fromRational
(/) = cxxInfixPhqApply"/"
recip = cxxFuncPhqApply"inv"
instance Floating PhqFuncTerm where
pi = PhqFnDimlessConst pi
exp = cxxFuncPhqApply"exp"
log = cxxFuncPhqApply"ln"
sqrt = cxxFuncPhqApply"sqrt"
a**(PhqFnDimlessConst 1) = a
a**(PhqFnDimlessConst x) = PhqFnFuncApply
(CXXFunc $ \base -> "(("++base++").to("++show x++"))") a
a**x = exp(log a*x)
sin = cxxFuncPhqApply"sin"
cos = cxxFuncPhqApply"cos"
tan = cxxFuncPhqApply"tan"
tanh = cxxFuncPhqApply"tanh"
asin = error "asin of physquantities not currently implemented."
acos = error "acos of physquantities not currently implemented."
atan = error "atan of physquantities not currently implemented."
sinh = error "sinh of physquantities not currently implemented."
cosh = error "cosh of physquantities not currently implemented."
asinh = error "asinh of physquantities not currently implemented."
acosh = error "acosh of physquantities not currently implemented."
atanh = error "atanh of physquantities not currently implemented."
-- 'seqPrunePhqFuncTerm' could be implemented much more efficiently by replacing
-- the lists with hash tables. This has low priority, since any function complicated
-- enough for this to take noteworthy time would always take yet a lot more time
-- when used with the cqtx algorithms.
-- instance Hashable PhqFuncTerm where
-- hash (PhqFnDimlessConst a) = hash"DimlessConst"`combine`hash a
-- hash (PhqFnPhysicalConst a) = hash"PhysicalConst"`combine`hash a
-- hash (PhqFnParameter a) = hash"Parameter"`combine`hash a
-- hash (PhqFnFuncApply f a) = hash"FuncApply"`combine`hash f`combine`hash a
-- hash (PhqFnInfixApply f a b) = hash"InfixApply"`combine`hash f`combine`hash a`combine`hash b
-- migrateRs (PhqFnTempRef n) = PhqFnTempRef (n+nLhsRefs)
-- migrateRs (PhqFnFuncApply f a) = PhqFnFuncApply f (migrateRs a)
-- migrateRs (PhqFnInfixApply f a b) = PhqFnInfixApply f (migrateRs a) (migrateRs b)
-- migrateRs c = c
-- | Obtain / eliminate common subexpressions.
seqPrunePhqFuncTerm :: PhqFuncTerm -- The original expression /e/, with possibly duplicate subexpressions.
-> ( PhqFuncTerm -- The \"trimmed\" version of /e/, in which duplicate expressions are replaced by variables, namely
, [(CXXExpression,PhqFuncTerm)] -- The common-subexpr–variables, with the names that were given to them.
)
seqPrunePhqFuncTerm = second ( map . first $ ("tmp"++) . show )
. prune . reverse . go []
where go :: [PhqFuncTerm] -> PhqFuncTerm -> [PhqFuncTerm]
go acc term
| (l, _:r) <- break((==term).expandRefs acc) acc
= refIn acc (length r) : acc
go acc x@(PhqFnFuncApply f a)
| isPrimitive a = x:acc
| otherwise
= let acc' = go acc a
in PhqFnFuncApply f (refIn acc' $ length acc'-1) : acc'
go acc x@(PhqFnInfixApply f a b)
= ifxGo acc x a b $ PhqFnInfixApply f
go acc x@(PhqFnInfixFoldOverIndexer ixd ifx fdInit fdand)
-- = ifxGo acc x fdInit fdand $ PhqFnInfixFoldOverIndexer ixd ifx
| isPrimitive fdInit = x:acc
| otherwise
= let acc' = go acc fdInit
in PhqFnInfixFoldOverIndexer ixd ifx (refIn acc' $ length acc'-1) fdand : acc'
go acc term = term : acc
ifxGo :: [PhqFuncTerm] -> PhqFuncTerm -> PhqFuncTerm -> PhqFuncTerm
-> ( PhqFuncTerm->PhqFuncTerm->PhqFuncTerm )
-> [ PhqFuncTerm ]
ifxGo acc x a b recomb
| isPrimitive a, isPrimitive b = x:acc
| isPrimitive a
= let acc' = go acc b
in recomb a (refIn acc' $ length acc'-1) : acc'
| isPrimitive b
= let acc' = go acc a
in recomb (refIn acc' $ length acc'-1) b : acc'
| otherwise
= let accL = go acc a
accR = go accL b
[nLhsRefs,nRhsRefs] = map length [accL,accR]
in recomb (refIn accR $ nLhsRefs-1)
(refIn accR $ nRhsRefs-1) : accR
expandRefs :: [PhqFuncTerm] -> PhqFuncTerm -> PhqFuncTerm
expandRefs rrl c = expnd c
where reflist = reverse rrl
expnd (PhqFnTempRef n _) = expnd $ reflist!!n
expnd (PhqFnFuncApply f t) = PhqFnFuncApply f $ expnd t
expnd (PhqFnInfixApply f a b) = PhqFnInfixApply f (expnd a) (expnd b)
expnd (PhqFnInfixFoldOverIndexer ixd ifx fdInit fdand)
= PhqFnInfixFoldOverIndexer ixd ifx (expnd fdInit) (expnd fdand)
expnd c = c
refIn rrl n = chase(PhqFnTempRef n $ "tmp"++show n)
where reflist = reverse rrl
chase (PhqFnTempRef n _)
| r'@(PhqFnTempRef _ _)<-reflist!!n = chase r'
chase c = c
prune :: [PhqFuncTerm] -> (PhqFuncTerm, [(Int, PhqFuncTerm)])
prune l = ( inlineIn $ last l, catMaybes $ map zap indexed )
where indexed = zip[0..] l
zap (n, e)
| n`elem`doomed = Nothing
| otherwise = Just (n, inlineIn e)
inlineIn (PhqFnTempRef n' _)
| n'`elem`doomed = inlineIn $ l!!n'
inlineIn (PhqFnFuncApply f e) = PhqFnFuncApply f $ inlineIn e
inlineIn (PhqFnInfixApply f a b) = PhqFnInfixApply f (inlineIn a) (inlineIn b)
inlineIn (PhqFnInfixFoldOverIndexer ixd ifx fdInit fdand)
= PhqFnInfixFoldOverIndexer ixd ifx (inlineIn fdInit) (inlineIn fdand)
inlineIn e = e
doomed = filter ((<=1) . length . referingTo) $ map fst indexed
referingTo n = filter (refersTo n) l
refersTo n (PhqFnTempRef n' _) = n==n'
refersTo n (PhqFnFuncApply _ e) = refersTo n e
refersTo n (PhqFnInfixApply _ a b) = refersTo n a || refersTo n b
refersTo n (PhqFnInfixFoldOverIndexer _ _ a b) = refersTo n a || refersTo n b
refersTo _ _ = False
calculateDimExpression :: (PhqIdf->PhqFuncTerm) -> DimExpression -> PhqFuncTerm
calculateDimExpression prImpl (DimExpression decomp)
= product $ map phqfImplement decomp
where phqfImplement (e, x) = implement e ** fromRational x
implement (PhqFnParamVal pr) = prImpl pr
implement PhqFnResultVal = PhqFnPhysicalConst "desiredret"
implement PhqDimlessConst = PhqFnDimlessConst 1
cqtxTermImplementation :: CXXExpression -> PhqFuncTerm -> CXXCode ( CXXCode() )
cqtxTermImplementation paramSource e = do
forM_ seqChain $ \(rn,se) ->
cxxSurround ("physquantity "++rn++" =") (showE se) (";")
return $ showE result
where (result,seqChain) = seqPrunePhqFuncTerm e
showE (PhqFnDimlessConst x) = cxxLine $ "("++show x++"*real1)"
showE (PhqFnPhysicalConst x) = cxxLine $ "("++x++")"
showE (PhqFnParameter x) = cxxLine $ argderefv paramSource x
showE (PhqFnTempRef _ tmpvn) = cxxLine $ tmpvn
showE (PhqFnFuncApply (CXXFunc f) a) = procCXXCode f $ showE a
showE (PhqFnInfixApply (CXXInfix f) a b) = cxxCombineBlockIndented f (showE a) (showE b)
showE (PhqFnInfixFoldOverIndexer (PhqVarIndexer idxId idxNm)
folder initE summandE ) = do
cxxLine $ "[&]() -> physquantity {"
cxxIndent 2 $ do
let i = indizesPrefix++idxNm
fdAcc = "foldvrb_over_"++idxNm
cxxSurround ("physquantity "++fdAcc++" = ") (showE initE) (";")
cxxLine $ "for(unsigned "++i++"=0\
\; "++i++"<"++i++indexRangePostfix++"\
\; ++"++i++") {"
cxxIndent 2 $ do
foldandResE <- cqtxTermImplementation paramSource summandE
combineAssigner folder fdAcc foldandResE
cxxLine $ "}"
cxxLine $ "return "++fdAcc++";"
cxxLine $ "}()"
-- | Transform a C / C++ infix into the corresponding combined assignment operator,
-- i.e. @+ x y@ maps to @x += y;@. When this is not possible, the explicit form is
-- used, e.g. @x pow@ will yield @x = pow(x, y);@.
combineAssigner :: CXXInfix -> CXXExpression -> CXXCode() -> CXXCode()
combineAssigner (CXXInfix f) acc inc
= case filter(/=' ') $ f"""" of
"+" -> cxxSurround (acc++" += ") inc (";")
"-" -> cxxSurround (acc++" -= ") inc (";")
"*" -> cxxSurround (acc++" *= ") inc (";")
"/" -> cxxSurround (acc++" /= ") inc (";")
_ -> cxxSurround (acc++" = ") (procCXXCode (f acc) inc) (";")
type CqtxConfig = ()
type CqtxCode = ReaderT CqtxConfig CXXCode
withDefaultCqtxConfig :: CqtxCode a -> CXXCode a
withDefaultCqtxConfig = flip runReaderT ()
newtype IdxablePhqDefVar x
= IdxablePhqDefVar {
indexisePhqDefVar :: PhqVarIndexer -> x }
data PhqVarIndexer = PhqVarIndexer
{ phqVarIndexerId :: Int
, phqVarIndexerName :: String
} deriving(Eq)
class (Floating a) => PhqfnDefining a where
sumOverIdx :: PhqVarIndexer
-> ((IdxablePhqDefVar a->a) -> a) -> a
instance PhqfnDefining PhqFuncTerm where
sumOverIdx i summand
= PhqFnInfixFoldOverIndexer i (cxxInfix"+") 0
. forbidDupFold i "sum"
. summand $ \(IdxablePhqDefVar e) -> e i
instance PhqfnDefining DimTracer where
sumOverIdx i summand = summand $
\(IdxablePhqDefVar e) -> e i
forbidDupFold :: PhqVarIndexer -> String -> PhqFuncTerm -> PhqFuncTerm
forbidDupFold i@(PhqVarIndexer _ nmm) fnm = go
where go(PhqFnInfixFoldOverIndexer i' ifx ini tm)
| i'==i = error $ "Duplicate "++fnm++" over index "++nmm
| otherwise = PhqFnInfixFoldOverIndexer i' ifx (go ini) (go tm)
go(PhqFnFuncApply f tm) = PhqFnFuncApply f $ go tm
go(PhqFnInfixApply ifx l r) = PhqFnInfixApply ifx (go l) (go r)
go e = e
class (Functor l, Foldable l) => IsolenList l where
perfectZipWith :: (a->b->c) -> l a -> l b -> l c
buildIsolenList :: (b->b) -> b -> l b
isoLength :: l a -> Int
isoLength = length . toList
enumFrom' :: Enum a => a -> l a
enumFrom' = buildIsolenList succ
(!!@) :: l a -> Int -> a
perfectZip :: IsolenList l => l a -> l b -> l (a,b)
perfectZip = perfectZipWith(,)
perfectZipWith3 :: IsolenList l => (a -> b -> c -> d) -> l a -> l b -> l c -> l d
perfectZipWith3 f la lb = perfectZipWith (uncurry . f) la . perfectZip lb
perfectZip3 :: IsolenList l => l a -> l b -> l c -> l (a,b,c)
perfectZip3 = perfectZipWith3 (,,)
data IsolenEnd a = P deriving (Show)
infixr 5 :.
data IsolenCons l a = a :. l a deriving (Show)
instance Functor IsolenEnd where
fmap _ P = P
instance (Functor l) => Functor (IsolenCons l) where
fmap f (x:.xs) = f x :. fmap f xs
instance Foldable IsolenEnd where
foldr _ = const
instance (Foldable l) => Foldable (IsolenCons l) where
foldr f ini (x:.xs) = x `f` foldr f ini xs
instance IsolenList IsolenEnd where
perfectZipWith _ P P = P
buildIsolenList _ _ = P
isoLength _ = 0
P !!@ (-1) = undefined
instance (IsolenList l) => IsolenList (IsolenCons l) where
perfectZipWith f (x:.xs) (y:.ys) = f x y :. perfectZipWith f xs ys
buildIsolenList f s = s :. buildIsolenList f (f s)
isoLength (_:.xs) = 1 + isoLength xs
(x:._) !!@ 0 = x
(_:.xs) !!@ n = xs !!@ (n-1)
-- instance IsolenList [] where -- obviously unsafe
-- perfectZipWith = zipWith
-- buildIsolenList = iterate
-- enumFrom' = enumFrom
-- singleList = id
-- unfoldToMatch [] _ _ = []
-- unfoldToMatch (_:xs) uff s = let (y,s') = uff s
-- in y : unfoldToMatch xs uff s'
-- buildIsolenList l1 f = unfoldToMatch l1 $ (\b->(b,b)) . f