cellml stuff added to handle Beeler-Reuter

Author: shahriariravanian

src/parse.jl

@@ -25,28 +25,28 @@ end
 """
     parse_cn(node)
 
-parse a <cn> node 
+parse a <cn> node
 """
 function parse_cn(node)
     # elements(node) != EzXML.Node[] && error("node cant have elements rn, check if </sep> is in the node")
     # this isnt great might want to check that `sep` actually shows up
     if haskey(node, "type") && elements(node) != EzXML.Node[]
         parse_cn_w_sep(node)
-    else 
+    else
         Meta.parse(node.content)
     end
 end
 
 """
     parse_cn_w_sep(node)
 
-parse a <cn type=".."> node 
+parse a <cn type=".."> node
 
 where type ∈ ["e-notation", "rational", "complex-cartesian", "complex-polar"]
 """
 function parse_cn_w_sep(node)
     # node = clean_attributes(node)
-    txts = findall("//text()", node)
+    txts = findall("text()", node)
     length(txts) != 2 && error("stop, collaborate, and listen!, problem with <cn>")
     x1, x2 = map(x -> Meta.parse(x.content), txts)
     t = node["type"]
@@ -58,14 +58,14 @@ function parse_cn_w_sep(node)
         Complex(x1, x2)
     elseif t == "complex-polar"
         x1 * exp(x2 * im)
-    else 
+    else
         error("$t in parse_cn_w_sep, somethings wrong")
     end
 end
 
 # """
 # remove all attributes neq to `type`
-# this is an issue for undefined namespaces 
+# this is an issue for undefined namespaces
 # https://github.com/JuliaIO/EzXML.jl/issues/156
 # """
 # function clean_attributes(node)
@@ -84,30 +84,61 @@ end
 """
     parse_ci(node)
 
-parse a <ci> node 
+parse a <ci> node
 """
 function parse_ci(node)
     c = Meta.parse(strip(node.content))
     Num(Variable(c))
 end
 
+########## Parse piecewise ###################################################
+
+function parse_piecewise(node)
+    return process_pieces(elements(node))
+end
+
+function process_pieces(pieces)
+    if length(pieces) == 1
+        return process_piece(pieces[1])
+    else
+        return process_piece(pieces[1], process_pieces(pieces[2:end]))
+    end
+end
+
+function process_piece(node)
+    if nodename(node) != "otherwise"
+        error("expect an otherwise")
+    else
+        return parse_node(firstelement(node))
+    end
+end
+
+function process_piece(node, otherwise)
+    if nodename(node) == "otherwise"
+        return parse_node(firstelement(node))
+    elseif nodename(node) == "piece"
+        c = parse_node.(elements(node))
+        return Symbolics.IfElse.ifelse(c[2] > 0.5, c[1], otherwise)
+    end
+end
+
 # """
 #     parse_piecewise(node)
 
-# parse a <piecewise> node 
+# parse a <piecewise> node
 # want to recursively call ifelse on the pieces
 # """
 # function parse_piecewise(node)
 #     es = elements(node)
-#     IfElse.ifelse(a, b, 
-#         IfElse.ifelse(c, d, 
+#     IfElse.ifelse(a, b,
+#         IfElse.ifelse(c, d,
 #             IfElse.ifelse(e, f, otherwise)))
 # end
 
 # """
 #     parse_piece(node)
 
-# parse a <piece> node 
+# parse a <piece> node
 # Each <piece> element contains exactly two children.
 # The conditional is the second child and the return is the first.
 # """
@@ -133,7 +164,7 @@ end
 """
     parse_bvar(node)
 
-parse a <bvar> node 
+parse a <bvar> node
 """
 function parse_bvar(node)
     es = elements(node)
@@ -156,20 +187,20 @@ end
 tagmap = Dict{String,Function}(
     "cn" => parse_cn,
     "ci" => parse_ci,
-    
+
     "degree" => x -> parse_node(x.firstelement), # won't work for all cases
     "bvar" => parse_bvar, # won't work for all cases
-    # "diff" => parse_diff, #inputs are 
+    # "diff" => parse_diff, #inputs are
 
-    # "piecewise" => parse_piecewise, 
-    # "piece" => parse_piece, 
-    # "otherwise" => x-> parse_node(x.firstelement), 
+    "piecewise" => parse_piecewise,
+    # "piece" => parse_piece,
+    # "otherwise" => x-> parse_node(x.firstelement),
 
     "apply" => parse_apply,
     "math" => x -> map(parse_node, elements(x)),
     "vector" => x -> map(parse_node, elements(x)),
 )
-    
+
 function custom_root(x)
     length(x) == 1 ? sqrt(x...) : Base.:^(x[2], x[1])
 end
@@ -180,45 +211,72 @@ function check_ivs(node)
     all(y -> y.content == x[1].content, x)
 end
 
-# need to check the arities 
+H(x) = Symbolics.IfElse.ifelse(x > 0, one(x), zero(x))
+const ϵ = eps(Float64)
+frac(x) = 0.5 - atan(cot(π * x)) / π
+heaviside_or(x) = length(x) == 1 ? x[1] : x[1] + heaviside_or(x[2:end]) - x[1] * heaviside_or(x[2:end])
+
+# need to check the arities
 # units handling??
 applymap = Dict{String,Function}(
     # eq sometimes needs to be ~ and sometimes needs to be =, not sure what the soln is
-    "eq" => x -> Symbolics.:~(x...), # arity 2, 
+    "eq" => x -> Symbolics.:~(x...), # arity 2,
     "times" => Base.prod, # arity 2, but prod fine
-    # "prod" => Base.prod, 
+    # "prod" => Base.prod,
     "divide" => x -> Base.:/(x...),
     "power" => x -> Base.:^(x...),
-    "root" => custom_root, 
+    "root" => custom_root,
     "plus" => x -> Base.:+(x...),
     "minus" => x -> Base.:-(x...),
-    "lt" => x -> Base.foldl(Base.:<, x),
-    "leq" => x -> Base.foldl(Base.:≤, x),
-    "geq" => x -> Base.foldl(Base.:≥, x),
-    "gt" => x -> Base.foldl(Base.:>, x),
+
+    # comparison functions implemented using the Heaviside function
+    "lt" => x -> H(x[2] - x[1] - ϵ),
+    "leq" => x -> H(x[2] - x[1]),
+    "geq" => x -> H(x[1] - x[2]),
+    "gt" => x -> H(x[1] - x[2] - ϵ),
+
+    # "lt" => x -> Base.foldl(Base.:<, x),
+    # "leq" => x -> Base.foldl(Base.:≤, x),
+    # "geq" => x -> Base.foldl(Base.:≥, x),
+    # "gt" => x -> Base.foldl(Base.:>, x),
+
     # "quotient" => x->Base.:div(x...), # broken, RoundingMode
     "factorial" => x -> Base.factorial(x...),
     "max" => x -> Base.max(x...),
     "min" => x -> Base.min(x...),
     "rem" => x -> Base.:rem(x...),
     "gcd" => x -> Base.:gcd(x...),
-    "and" => Base.:&,
-    "or" => Base.:|,
-    "xor" => Base.:⊻,
-    "not" => Base.:!,
+
+    "and" => x -> Base.:*(x...),
+    "or" => heaviside_or,
+    "xor" => x -> x[1]*(one(x[2])-x[2]) + (one(x[1])-x[1])*x[2],
+    "not" => x -> one(x[1]) - x[1],
+
+    # "and" => x -> Base.:&(x...),
+    # "or" => Base.:|,
+    # "xor" => Base.:⊻,
+    # "not" => Base.:!,
+
     "abs" => x -> Base.abs(x...),
     "conjugate" => Base.conj,
     "arg" => Base.angle,
     "real" => Base.real,
     "imaginary" => Base.imag,
     "lcm" =>  x -> Base.lcm(x...),
-    "floor" => x -> Base.floor(x...),
-    "ceiling" => x -> Base.ceil(x...),
+
+    "floor" => x -> x[1] - frac(x[1]),
+    "ceiling" => x -> x[1] - frac(x[1]) + one(x[1]),
+    "round" => x -> x[1] + 0.5 - frac(x[1] + 0.5),
+
+    # "floor" => x -> Base.floor(x...),
+    # "ceiling" => x -> Base.ceil(x...),
+    # "round" => x -> Base.round(x...),
+
     "inverse" => Base.inv,
     "compose" => x -> Base.:∘(x...),
     "ident" => Base.identity,
     "approx" => x -> Base.:≈(x...),
-    
+
     "sin" => x -> Base.sin(x...),
     "cos" => x -> Base.cos(x...),
     "tan" => x -> Base.tan(x...),
@@ -245,7 +303,7 @@ applymap = Dict{String,Function}(
     "arccoth" => x -> Base.acoth(x...),
 
     "exp" => x -> Base.exp(x...),
-    "log" => x -> Base.log10(x...), #  todo handle <logbase> 
+    "log" => x -> Base.log10(x...), #  todo handle <logbase>
     "ln" => x -> Base.log(x...),
 
     "mean" => Statistics.mean,