Improvements for performance and for synchronous operators

Author: Hilding-Elmqvist

src/BLTandPantelides.jl

@@ -27,18 +27,18 @@ using ..BLTandPantelidesUtilities
 const log = false
 
 """
-    function augmentPath!(G, i, assign, vColour, eColour, vActive)
+    function augmentPath!(G, i, assign, vColour, eColour, vPassive)
 Construction of augmenting path
 
 Reference:
 Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
 of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988). 
 """
-function augmentPath!(G, i, assign, vColour, eColour, vActive)
+function augmentPath!(G, i, assign, vColour, eColour, vPassive)
     # returns pathFound
     # assign: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
     # i: E-node
-    # vActive: set to false has the same effect as deleting V-node and corresponding edges
+    # vPassive: set to != 0 has the same effect as deleting V-node and corresponding edges
     # j: V-node
 
     if log
@@ -50,7 +50,7 @@ function augmentPath!(G, i, assign, vColour, eColour, vActive)
 
     # If a V-node j exists such that edge (i-j) exists and assign[j] == 0
     for j in G[i]
-        if vActive[j] && assign[j] == 0
+        if vPassive[j] == 0 && assign[j] == 0
             pathFound = true
             assign[j] = i
             return pathFound
@@ -59,10 +59,10 @@ function augmentPath!(G, i, assign, vColour, eColour, vActive)
 
     # For every j such that edge (i-j) exists and j is uncoloured
     for j in G[i]
-        if vActive[j] && !vColour[j]
+        if vPassive[j] == 0 && !vColour[j]
             vColour[j] = true
             k = assign[j]
-            pathFound = augmentPath!(G, k, assign, vColour, eColour, vActive)
+            pathFound = augmentPath!(G, k, assign, vColour, eColour, vPassive)
 
             if pathFound 
                 assign[j] = i
@@ -74,11 +74,11 @@ function augmentPath!(G, i, assign, vColour, eColour, vActive)
 end
 
 
-function checkAssign(assign, VSizes, VTypes, ESizes, ETypes, equationsInfix, variableNames, A, vActive=[a == 0 for a in A])
+function checkAssign(assign, VSizes, VTypes, ESizes, ETypes, equationsInfix, variableNames, A, vPassive=A)
     println("Checking assignment")
     assignmentOK = true
     for j in 1:length(assign)
-        if vActive[j]
+        if vPassive[j] == 0
             i = assign[j]
             if i > 0 && VSizes[j] != ESizes[i]
                 assignmentOK = false
@@ -114,10 +114,11 @@ function matching(G, M, vActive=fill(true, M))
     assign::Array{Int,1} = fill(0, M)
     eColour::Array{Bool,1} = fill(false, length(G))
     vColour::Array{Bool,1} = fill(false, M)
+    vPassive::Array{Int,1} = [if va; 0 else 1 end for va in vActive]
     for i in 1:length(G)
         fill!(eColour, false)
         fill!(vColour, false)
-        pathFound = augmentPath!(G, i, assign, vColour, eColour, vActive)
+        pathFound = augmentPath!(G, i, assign, vColour, eColour, vPassive)
     end
     return assign
 end
@@ -142,18 +143,28 @@ of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
 function pantelides!(G, M, A)
     assign::Array{Int,1} = fill(0, M)
     B::Array{Int,1} = fill(0, length(G))
+    eColour::Array{Bool,1} = fill(false, length(G))
+    vColour::Array{Bool,1} = fill(false, M)
     N = length(G)
     N2 = N
     for k in 1:N2
         pathFound = false
         i = k
         while !pathFound
             # Delete all V-nodes with A[.] != 0 and all their incidence edges from the graph
-            vActive::Array{Bool,1} = [a == 0 for a in A]
             # Designate all nodes as "uncoloured"
-            eColour::Array{Bool,1} = fill(false, length(G))
-            vColour::Array{Bool,1} = fill(false, M)
-            pathFound = augmentPath!(G, i, assign, vColour, eColour, vActive)
+            if length(eColour) == length(G)
+                fill!(eColour, false)
+            else
+                eColour = fill(false, length(G))
+            end
+            if length(vColour) == length(M)
+                fill!(vColour, false)
+            else
+                vColour = fill(false, M)
+            end
+
+            pathFound = augmentPath!(G, i, assign, vColour, eColour, A)
             if !pathFound
                 if log
                     println("\nDifferentiate:")

src/Symbolic.jl

@@ -20,7 +20,6 @@ using Unitful
 using Measurements
 using MonteCarloMeasurements
 
-
 """
     e = removeBlock(ex)
 
@@ -40,6 +39,20 @@ function removeBlock(ex::Expr)
     end
 end
 
+prependPar(ex, prefix, parameters=[], inputs=[]) = ex
+
+function prependPar(ex::Symbol, prefix, parameters=[], inputs=[]) 
+    if prefix == nothing; ex elseif ex in [:time, :instantiatedModel, :_leq_mode, :_x]; ex else Expr(:ref, prefix, QuoteNode(ex)) end
+end
+
+function prependPar(ex::Expr, prefix, parameters=[], inputs=[])
+    if isexpr(ex, :.)
+        Expr(:ref, prependPar(ex.args[1], prefix, parameters, inputs), QuoteNode(ex.args[2].value))
+    else
+        Expr(ex.head, [prependPar(arg, prefix, parameters, inputs) for arg in ex.args]...)
+    end
+end
+
 """
     e = makeDerVar(ex)
 
@@ -48,20 +61,29 @@ Recursively converts der(x) to Symbol(:(der(x))) in expression `ex`
 * `ex`: Expression or array of expressions
 * `return `e`: ex with der(x) converted 
 """
-makeDerVar(ex, parameters, inputs=[]) = if typeof(ex) in [Symbol, Expr] && ex in parameters; prepend(ex, :(_p)) 
-    elseif typeof(ex) in [Symbol, Expr] && ex in inputs; prepend(ex, :(_p)) else ex end
+makeDerVar(ex, parameters, inputs, evaluateParameters=false) = if typeof(ex) in [Symbol, Expr] && 
+    ( ex in keys(parameters) || ex in keys(inputs) ); prependPar(ex, :(_p), parameters, inputs) 
+    else ex end
 
-function makeDerVar(ex::Expr, parameters=[], inputs=[])
+function makeDerVar(ex::Expr, parameters, inputs, evaluateParameters=false)
     if ex.head == :call && ex.args[1] == :der
         Symbol(ex)
-	elseif isexpr(ex, :.) && ex in parameters
-		prepend(ex, :(_p))
-	elseif isexpr(ex, :.) && ex in inputs
-		prepend(ex, :(_p))
+	elseif isexpr(ex, :.) && ex in keys(parameters)
+        if evaluateParameters
+            parameters[ex]
+        else
+            prependPar(ex, :(_p), parameters, inputs)
+        end
+	elseif isexpr(ex, :.) && ex in keys(inputs)
+        if evaluateParameters
+            inputs[ex]
+        else
+            prependPar(ex, :(_p), parameters, inputs)
+        end
     elseif ex.head == :.
         Symbol(ex)
     else
-        Expr(ex.head, [makeDerVar(arg, parameters, inputs) for arg in ex.args]...)
+        Expr(ex.head, [makeDerVar(arg, parameters, inputs, evaluateParameters) for arg in ex.args]...)
     end
 end
 
@@ -106,6 +128,7 @@ nClocks = 0
 nSamples = 0
 previousVars = []
 preVars = []
+holdVars = []
 
 function resetEventCounters()
     global nCrossingFunctions
@@ -114,22 +137,25 @@ function resetEventCounters()
     global nSamples 
 	global previousVars
     global preVars
+    global holdVars
     nCrossingFunctions = 0
     nAfter = 0
     nClocks = 0
     nSamples = 0
 	previousVars = []
     preVars = []
+    holdVars = []
 end
 
 function getEventCounters()
     global nCrossingFunctions
     global nAfter
     global nClocks
-    global nSamples 
+    global nSamples
 	global previousVars
     global preVars
-    return (nCrossingFunctions, nAfter, nClocks, nSamples, previousVars, preVars)
+    global holdVars
+    return (nCrossingFunctions, nAfter, nClocks, nSamples, previousVars, preVars, holdVars)
 end
 
 substituteForEvents(ex) = ex
@@ -141,6 +167,7 @@ function substituteForEvents(ex::Expr)
     global nSamples 
 	global previousVar
     global preVars
+    global holdVars
     if ex.head in [:call, :kw]
         if ex.head == :call && ex.args[1] == :positive
             nCrossingFunctions += 1
@@ -172,6 +199,15 @@ function substituteForEvents(ex::Expr)
             else
                 error("The previous function presently takes two arguments: $ex")
             end
+        elseif ex.head == :call && ex.args[1] == :hold
+            push!(holdVars, ex.args[2])
+            nHold = length(holdVars)
+            if length(ex.args) == 3
+                :(hold($(substituteForEvents(ex.args[2])), $(substituteForEvents(ex.args[3])), instantiatedModel, $nHold))
+            else
+#                error("The hold function takes two or three arguments, hold(v, clock) or hold(expr, start, clock) : $ex")
+                error("The hold function takes two arguments, hold(v, clock): $ex")
+            end
         elseif ex.head == :call && ex.args[1] in [:initial, :terminal]
             if length(ex.args) == 1
                 :($(ex.args[1])(instantiatedModel))
@@ -245,8 +281,8 @@ Finds the linear `factor` and `rest` if `ex` is `linear` with regards to `x` (ex
 * `return (rest, factor, linear)`:  
 """
 linearFactor(ex, x) = (ex, 0, true)
-linearFactor(ex::Symbol, x) =  if ex == x; (0, 1, true) else (ex, 0, true) end
-function linearFactor(ex::Expr, x)
+linearFactor(ex::Symbol, x::Incidence) =  if ex == x; (0, 1, true) else (ex, 0, true) end
+function linearFactor(ex::Expr, x::Incidence)
     if ex.head == :block
         linearFactor(ex.args[1], x)    
     elseif ex.head == :macrocall && ex.args[1] == Symbol("@u_str")
@@ -374,7 +410,7 @@ Tests if `ex` is `linear` with regards to `x` (ex == factor*x + rest) and checks
 * `ex`: Expression
 * `return (linear, constant)` 
 """
-function isLinear(equ::Expr, x)
+function isLinear(equ::Expr, x::Incidence)
     (rest, factor, linear) = linearFactor(equ, x)
     (linear, typeof(factor) != Expr)
 end
@@ -389,7 +425,7 @@ If `ex` is `linear` with regards to all `incidence` (Symbols and der(...)), the
 * `return (incidence, coefficients, rest, linear)` 
 """
 
-function getCoefficients(ex)
+function getCoefficients(ex::Expr)
     incidence = Incidence[]
     findIncidence!(ex, incidence)
     coefficients = []
@@ -400,10 +436,12 @@ function getCoefficients(ex)
         findIncidence!(coefficients, crossIncidence)
         if x in crossIncidence
             linear = false
+            break
         end
         (rest, factor, linearRest) = linearFactor(rest, x)
         if !linearRest
             linear = false
+            break
         end
         push!(coefficients, factor)
     end
@@ -412,6 +450,10 @@ end
 
 substitute(substitutions, ex) = begin nex = get(substitutions, ex, ex); if nex != ex; nex = substitute(substitutions, nex) else nex end; nex end 
 
+substitute(substitutions, ex::Vector{Symbol}) = [substitute(substitutions, e) for e in ex]
+
+substitute(substitutions, ex::Vector{Expr}) = [substitute(substitutions, e) for e in ex]
+
 substitute(substitutions, ex::MonteCarloMeasurements.StaticParticles) = ex
 
 substitute(substitutions, ex::Measurements.Measurement{Float64}) = ex