Merge pull request #93 from brianguenter/AddConditionalsBranch

Patch to correct bug introduced when added conditionals in PR #90 (comment)

Project.toml

@@ -1,7 +1,7 @@
 name = "FastDifferentiation"
 uuid = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
 authors = ["BrianGuenter"]
-version = "0.3.17"
+version = "0.4.1"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

README.md

@@ -71,26 +71,44 @@ If you use FD in your work please share the functions you differentiate with me.
 **A**: If you multiply a matrix of **FD** variables times a vector of **FD** variables the matrix vector multiplication loop is effectively unrolled into scalar expressions. Matrix operations on large matrices will generate large executables and long preprocessing time. **FD** functions with up 10⁵ operations should still have reasonable preprocessing/compilation times (approximately 1 minute on a modern laptop) and good run time performance.
 
 **Q**: Does **FD** support conditionals?  
-**A**: **FD** does not yet support conditionals that involve the variables you are differentiating with respect to. You can do this:
+**A**: As of version 0.4.1 **FD** expressions may contain conditionals which involve variables. However, you cannot yet differentiate an expression containing conditionals. A future PR will allow you to differentiate conditional expressions. 
+
+You can use either the builtin `ifelse` function or a new function `if_else`. `ifelse` will evaluate both the true and false branches. By contrast `if_else` has the semantics of `if...else...end`; only one of the true or false branches will be executed. 
+
+This is useful when your conditional is used to prevent exceptions because of illegal input values:
 ```julia
-@variables x y #create FD variables
+julia> f = if_else(x<0,NaN,sqrt(x))
+(if_else  (x < 0) NaN sqrt(x))
+
+julia> g = make_function([f],[x])
 
-julia> f(a,b,c) = a< 1.0 ? cos(b) : sin(c)
-f (generic function with 2 methods)
 
-julia> f(0.0,x,y)
-cos(x)
+julia> g([-1])
+1-element Vector{Float64}:
+ NaN
 
-julia> f(1.0,x,y)
-sin(y)
+julia> g([2.0])
+1-element Vector{Float64}:
+ 1.4142135623730951
+end
 ```
-but you can't do this:
+In this case you wouldn't want to use `ifelse` because it evaluates both the true and false branches and causes a runtime exception:
 ```julia
-julia> f(a,b) = a < b ? cos(a) : sin(b)
-f (generic function with 2 methods)
+julia> f = ifelse(x<0,NaN,sqrt(x))
+(ifelse  (x < 0) NaN sqrt(x))
 
-julia> f(x,y)
-ERROR: MethodError: no method matching isless(::FastDifferentiation.Node{Symbol, 0}, ::FastDifferentiation.Node{Symbol, 0})
+julia> g = make_function([f],[x])
+...
+
+julia> g([-1])
+ERROR: DomainError with -1.0:
+sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
+```
+
+However, you cannot yet compute derivatives of expressions that contain conditionals:
+```julia
+julia> jacobian([f],[x,y])
+ERROR: Your expression contained ifelse. FastDifferentiation does not yet support differentiation through ifelse or any of these conditionals (max, min, copysign, &, |, xor, <, >, <=, >=, !=, ==, signbit, isreal, iszero, isfinite, isnan, isinf, isinteger, !)
 ```
 
 # Release Notes

docs/src/index.md

@@ -7,12 +7,16 @@ CurrentModule = FastDifferentiation
 [![Build Status](https://github.com/brianguenter/FastDifferentiation.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/brianguenter/FastDifferentiation.jl/actions/workflows/CI.yml?query=branch%3Amain)
 
 
-FastDifferentiation (**FD**) is a package for generating efficient executables to evaluate derivatives of Julia functions. It can also generate efficient true symbolic derivatives for symbolic analysis. Unlike forward and reverse mode automatic differentiation **FD** automatically generates efficient derivatives for arbitrary function types: ℝ¹->ℝ¹, ℝ¹->ℝᵐ, ℝⁿ->ℝ¹, and ℝⁿ->ℝᵐ, m≠1,n≠1. 
+FastDifferentiation (**FD**) is a package for generating efficient executables to evaluate derivatives of Julia functions. It can also generate efficient true symbolic derivatives for symbolic analysis. Unlike forward and reverse mode automatic differentiation **FD** automatically generates efficient derivatives for arbitrary function types: ℝ¹->ℝ¹, ℝ¹->ℝᵐ, ℝⁿ->ℝ¹, and ℝⁿ->ℝᵐ. 
+
+For f:ℝⁿ->ℝᵐ with n,m large **FD** may have better performance than conventional AD algorithms because the **FD** algorithm finds expressions shared between partials and computes them only once. In some cases **FD** derivatives can be as efficient as manually coded derivatives (see the Lagrangian dynamics example in the [D*](https://www.microsoft.com/en-us/research/publication/the-d-symbolic-differentiation-algorithm/) paper or the [Benchmarks](@ref) section of the documentation for another example).
+
+
+ **FD** may take much less time to compute symbolic derivatives than Symbolics.jl even in the ℝ¹->ℝ¹ case. The executables generated by **FD** may also be much faster (see [Symbolic Processing](@ref)). 
+
 
-For f:ℝⁿ->ℝᵐ with n,m large FD may have better performance than conventional AD algorithms because the **FD** algorithm finds expressions shared between partials and computes them only once. In some cases **FD** derivatives can be as efficient as manually coded derivatives (see the Lagrangian dynamics example in the [D*](https://www.microsoft.com/en-us/research/publication/the-d-symbolic-differentiation-algorithm/) paper or the [Benchmarks](@ref) section of the documentation for another example).
 
 
- **FD** may take much less time to compute symbolic derivatives than Symbolics.jl even in the ℝ¹->ℝ¹ case. The executables generated by **FD** may also be much faster (see [Symbolic Processing](@ref)[^1]. 
 
 You should consider using FastDifferentiation when you need: 
 * a fast executable for evaluating the derivative of a function and the overhead of the preprocessing/compilation time is swamped by evaluation time.
@@ -36,8 +40,48 @@ This is **beta** software being modified on a daily basis. Expect bugs and frequ
 ## Notes about special derivatives
 The derivative of `|u|` is `u/|u|` which is NaN when `u==0`. This is not a bug. The derivative of the absolute value function is undefined at 0 and the way **FD** signals this is by returning NaN.
 
+## Conditionals
+
+As of version 0.4.1 **FD** allows you to create expressions with conditionals using either the builtin `ifelse` function or a new function `if_else`. `ifelse` will evaluate both inputs. By contrast `if_else` has the semantics of `if...else...end`; only the true or false branch will be executed. This is useful when your conditional is used to prevent exceptions because of illegal input values:
+```julia
+julia> f = if_else(x<0,NaN,sqrt(x))
+(if_else  (x < 0) NaN sqrt(x))
+
+julia> g = make_function([f],[x])
+
+
+julia> g([-1])
+1-element Vector{Float64}:
+ NaN
+
+julia> g([2.0])
+1-element Vector{Float64}:
+ 1.4142135623730951
+end
+```
+In this case you wouldn't want to use `ifelse` because it evaluates both the true and false branches and causes a runtime exception:
+```julia
+julia> f = ifelse(x<0,NaN,sqrt(x))
+(ifelse  (x < 0) NaN sqrt(x))
+
+julia> g = make_function([f],[x])
+...
+
+julia> g([-1])
+ERROR: DomainError with -1.0:
+sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
+```
+
+However, you cannot yet compute derivatives of expressions that contain conditionals:
+```julia
+julia> jacobian([f],[x,y])
+ERROR: Your expression contained ifelse. FastDifferentiation does not yet support differentiation through ifelse or any of these conditionals (max, min, copysign, &, |, xor, <, >, <=, >=, !=, ==, signbit, isreal, iszero, isfinite, isnan, isinf, isinteger, !)
+```
+This may be a breaking change for some users. In previous versions the expression `x==y` returned a `Bool`. Some data structures, such as `Dict` use `==` by default to determine if two entries are the same. This will no longer work since `x==y` will now return an expression graph. Use an `IdDict` instead since this uses `===`.
+
+A future PR will add support for differentiating through conditionals.
+
 
-[^1]: I am working with the SciML team to see if it is possible to integrate **FD** differentiation directly into Symbolics.jl.
 
 
 

docs/src/limitations.md

@@ -1,72 +1,5 @@
 # Limitations
-**FD** does not support expressions with conditionals on **FD** variables. For example, you can do this:
-```julia
-julia> f(a,b,c) = a< 1.0 ? cos(b) : sin(c)
-f (generic function with 2 methods)
-
-julia> f(0.0,x,y)
-cos(x)
-
-julia> f(1.0,x,y)
-sin(y)
-```
-but you can't do this:
-```julia
-julia> f(a,b) = a < b ? cos(a) : sin(b)
-f (generic function with 2 methods)
-
-julia> f(x,y)
-ERROR: MethodError: no method matching isless(::FastDifferentiation.Node{Symbol, 0}, ::FastDifferentiation.Node{Symbol, 0})
-
-Closest candidates are:
-  isless(::Any, ::DataValues.DataValue{Union{}})
-   @ DataValues ~/.julia/packages/DataValues/N7oeL/src/scalar/core.jl:291
-  isless(::S, ::DataValues.DataValue{T}) where {S, T}
-   @ DataValues ~/.julia/packages/DataValues/N7oeL/src/scalar/core.jl:285
-  isless(::DataValues.DataValue{Union{}}, ::Any)
-   @ DataValues ~/.julia/packages/DataValues/N7oeL/src/scalar/core.jl:293
-  ...
-```
-This is because the call `f(x,y)` creates an expression graph. At graph creation time the **FD** variables `x,y` are unevaluated variables with no specific value so they cannot be compared with any other value.
-
-The algorithm can be extended to work with conditionals applied to **FD** variables but the processing time and graph size may grow exponentially with conditional nesting depth. A future version may allow for limited conditional nesting. See [Future Work](@ref) for a potential long term solution to this problem.
 
 **FD** does not support looping internally. All operations with loops, such as matrix vector multiplication, are unrolled into scalar operations. The corresponding executable functions generated by `make_function` have size proportional to the number of operations. 
 
-Expressions with ≈10⁵ scalar operations have reasonable symbolic preprocessing and compilation times.  Beyond this size LLVM compilation time can become extremely long and eventually the executables become so large that their caching behavior is not good and performance declines. 
-
-A possible solution to this problem is to do what is called rerolling: detecting repreating indexing patterns in the **FD** expressions and automatically generating loops to replace inlined code. This rerolling step would be performed on the **FD** expressions graphs before function compliation.
-
-It is not necessary to completely undo the unrolling back to the original expresssion, just to reduce code size enough to get reasonable compilation times and better caching behavior.
-
-For example, in this matrix vector multiplication
-```julia
-
-julia> a = make_variables(:a,3,3)
-3×3 Matrix{FastDifferentiation.Node}:
- a1_1  a1_2  a1_3
- a2_1  a2_2  a2_3
- a3_1  a3_2  a3_3
-
-julia> b = make_variables(:b,3)
-3-element Vector{FastDifferentiation.Node}:
- b1
- b2
- b3
-
-julia> a*b
-3-element Vector{Any}:
- (((a1_1 * b1) + (a1_2 * b2)) + (a1_3 * b3))
- (((a2_1 * b1) + (a2_2 * b2)) + (a2_3 * b3))
- (((a3_1 * b1) + (a3_2 * b2)) + (a3_3 * b3))
-```
-the goal is to replace  concrete index numbers with symbolic index variables that represent offsets rather than absolute indices
-```julia
-
-[
-  a[i,j]*b[j] + a[i, j+1]*b[j+1] + a[i,j+2]*b[j+2]
-  a[i+1,j]*b[j] + a[i+1, j+1]*b[j+1] + a[i+1,j+2]*b[j+2]
-  a[i+2,j]*b[j] + a[i+2, j+1]*b[j+1] + a[i+2,j+2]*b[j+2]
-]
-```
-and then to extract looping structure from these indices.
+Expressions with ≈10⁵ scalar operations have reasonable symbolic preprocessing and compilation times.  Beyond this size LLVM compilation time can become extremely long and eventually the executables become so large that their caching behavior is not good and performance declines.

src/CodeGeneration.jl

@@ -8,12 +8,68 @@ then `sparsity = (nelts-nzeros)/nelts`.
 
 Frequently used in combination with a call to `make_function` to determine whether to set keyword argument `init_with_zeros` to false."""
 function sparsity(sym_func::AbstractArray{<:Node})
-    zeros = mapreduce(x -> is_zero(x) ? 1 : 0, +, sym_func)
+    zeros = mapreduce(x -> is_identically_zero(x) ? 1 : 0, +, sym_func)
     tot = prod(size(sym_func))
     return zeros == 0 ? 1.0 : (tot - zeros) / tot
 end
 export sparsity
 
+
+function _dag_to_function!(node, local_body, variable_to_index, node_to_var)
+
+    tmp = get(node_to_var, node, nothing)
+
+    if tmp === nothing #if node not in node_to_var then it hasn't been visited. Otherwise it has so don't recurse.
+        node_to_var[node] = node_symbol(node, variable_to_index)
+
+        if is_tree(node)
+            if value(node) === if_else #special case code generation for if...else. Need to generate nested code so only the statements in the true or false branch will be executed.
+                true_body = Expr(:block)
+                false_body = Expr(:block)
+                if_cond_var = _dag_to_function!(children(node)[1], local_body, variable_to_index, node_to_var)
+
+                true_node = children(node)[2]
+                false_node = children(node)[3]
+
+                if is_leaf(true_node) #handle leaf nodes properly 
+                    if is_constant(true_node)
+                        temp_val = value(true_node)
+                    else
+                        temp_val = node_to_var[true_node]
+                    end
+
+                    push!(true_body.args, :($(gensym(:s)) = $(temp_val))) #seems roundabout to use an assignment when really just want the value of the node but couldn't figure out how to make this work with Expr
+                else
+                    _dag_to_function!(children(node)[2], true_body, variable_to_index, node_to_var)
+                end
+
+                if is_leaf(false_node)
+                    if is_constant(false_node)
+                        temp_val = value(false_node)
+                    else
+                        temp_val = node_to_var[false_node]
+                    end
+                    push!(false_body.args, :($(gensym(:s)) = $(temp_val))) #seems roundabout to use an assignment when really just want the value of the node but couldn't figure out how to make this work with Expr
+                else
+                    _dag_to_function!(children(node)[3], false_body, variable_to_index, node_to_var)
+                end
+
+                statement = :($(node_to_var[node]) = if $(if_cond_var)
+                    $(true_body)
+                else
+                    $(false_body)
+                end)
+            else
+                args = _dag_to_function!.(children(node), Ref(local_body), Ref(variable_to_index), Ref(node_to_var))
+                statement = :($(node_to_var[node]) = $(Symbol(value(node)))($(args...)))
+            end
+            push!(local_body.args, statement)
+        end
+    end
+
+    return node_to_var[node]
+end
+
 """
     function_body!(
         dag::Node,
@@ -37,31 +93,12 @@ end
 ```
 and the second return value will be the constant value.
 """
-function function_body!(dag::Node, variable_to_index::IdDict{Node,Int64}, node_to_var::Union{Nothing,IdDict{Node,Union{Symbol,Real,Expr}}}=nothing)
+function function_body!(dag::Node, variable_to_index::IdDict{Node,Int64}, node_to_var::Union{Nothing,IdDict{Node,Union{Symbol,Real,Expr}}}=nothing, body::Expr=Expr(:block))
     if node_to_var === nothing
         node_to_var = IdDict{Node,Union{Symbol,Real,Expr}}()
     end
 
-    body = Expr(:block)
-
-    function _dag_to_function(node)
-
-        tmp = get(node_to_var, node, nothing)
-
-        if tmp === nothing #if node not in node_to_var then it hasn't been visited. Otherwise it has so don't recurse.
-            node_to_var[node] = node_symbol(node, variable_to_index)
-
-            if is_tree(node)
-                args = _dag_to_function.(children(node))
-                statement = :($(node_to_var[node]) = $(Symbol(value(node)))($(args...)))
-                push!(body.args, statement)
-            end
-        end
-
-        return node_to_var[node]
-    end
-
-    return body, _dag_to_function(dag)
+    return body, _dag_to_function!(dag, body, variable_to_index, node_to_var)
 end
 
 function zero_array_declaration(array::StaticArray{S,<:Any,N}) where {S,N}
@@ -124,8 +161,8 @@ function make_Expr(func_array::AbstractArray{T}, input_variables::AbstractVector
     node_to_var = IdDict{Node,Union{Symbol,Real,Expr}}()
     body = Expr(:block)
 
-    num_zeros = count(is_zero, (func_array))
-    num_const = count((x) -> is_constant(x) && !is_zero(x), func_array)
+    num_zeros = count(is_identically_zero, (func_array))
+    num_const = count((x) -> is_constant(x) && !is_identically_zero(x), func_array)
 
 
     zero_threshold = 0.5
@@ -166,7 +203,7 @@ function make_Expr(func_array::AbstractArray{T}, input_variables::AbstractVector
     for (i, node) in pairs(func_array)
         # skip all terms that we have computed above during construction
         if is_constant(node) && initialization_strategy === :const || # already initialized as constant above
-           is_zero(node) && (initialization_strategy === :zero || !init_with_zeros) # was already initialized as zero above or we don't want to initialize with zeros
+           is_identically_zero(node) && (initialization_strategy === :zero || !init_with_zeros) # was already initialized as zero above or we don't want to initialize with zeros
             continue
         end
         node_body, variable = function_body!(node, node_to_index, node_to_var)
@@ -185,11 +222,11 @@ function make_Expr(func_array::AbstractArray{T}, input_variables::AbstractVector
 
     # wrap in function body
     if in_place
-        return :((result, input_variables) -> @inbounds begin
+        return :((result, input_variables::AbstractArray) -> @inbounds begin
             $body
         end)
     else
-        return :((input_variables) -> @inbounds begin
+        return :((input_variables::AbstractArray) -> @inbounds begin
             $body
         end)
     end
@@ -248,9 +285,9 @@ function make_Expr(A::SparseMatrixCSC{T,Ti}, input_variables::AbstractVector{S},
         push!(body.args, :(return result))
 
         if in_place
-            return :((result, input_variables) -> $body)
+            return :((result, input_variables::AbstractArray) -> $body)
         else
-            return :((input_variables) -> $body)
+            return :((input_variables::AbstractArray) -> $body)
         end
     end
 end
@@ -316,7 +353,20 @@ function make_function(func_array::AbstractArray{T}, input_variables::AbstractVe
     vars = variables(func_array) #all unique variables in func_array
     all_input_vars = vcat(input_variables...)
 
-    @assert vars ⊆ all_input_vars "Some of the variables in your function (the func_array argument) were not in the input_variables argument. Every variable that is used in your function must have a corresponding entry in the input_variables argument."
+    #Because FD defines == operator for Node, which does not return a boolean, many builtin Julia functions will not work as expected. For example:
+    #  vars ⊆ all_input_vars errors because internally issubset tests for equality between the node values using ==, not ===. == returns a Node value but the issubset function expects a Bool.
+
+    temp = Vector{eltype(vars)}(undef, 0)
+
+    input_dict = IdDict(zip(all_input_vars, all_input_vars))
+    for one_var in vars
+        value = get(input_dict, one_var, nothing)
+        if value === nothing
+            push!(temp, one_var)
+        end
+    end
+
+    @assert length(temp) == 0 "The variables $temp were not in the input_variables argument to make_function. Every variable that is used in your function must have a corresponding entry in the input_variables argument."
 
     @RuntimeGeneratedFunction(make_Expr(func_array, all_input_vars, in_place, init_with_zeros))
 end