# How to make this function compatible with ForwardDiff?

**URL:** <https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415>\
**Category:** General Usage\
**Tags:** forwarddiff, autodiff\
**Created:** [August 31, 2021, 9:56am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415 "2021-08-31T09:56:50Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [August 31, 2021, 9:56am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/1 "2021-08-31T09:56:50Z")

</div>

How can I make this function compatible with `ForwardDiff` with respect to `τ`?

```julia
using ForwardDiff

randmat = rand(10, 2)
sto = similar(randmat)

function claytonsample!(sto, τ; randmat=randmat)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

```

right now I get this error:

```julia
julia> ForwardDiff.derivative(τ -> claytonsample!(sto, τ), 0.3)
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#81#82", Float64}, Float64, 1})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::VectorizationBase.Double{T}) where T<:Union{Float16, Float32, Float64, VectorizationBase.Vec{var"#s32", var"#s31"} where {var"#s32", var"#s31"<:Union{Float16, Float32, Float64}}, VectorizationBase.VecUnroll{var"#s30", var"#s29", var"#s28", V} where {var"#s30", var"#s29", var"#s28"<:Union{Float16, Float32, Float64}, V<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit, VectorizationBase.AbstractSIMD{var"#s29", var"#s28"}}}} at /Users/amrods/.julia/packages/VectorizationBase/xmoDw/src/special/double.jl:84
  ...

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [August 31, 2021, 10:49am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/2 "2021-08-31T10:49:57Z")

</div>

Strictly speaking you can’t: this is a function with side effects so it’s not completely clear what it means to differentiate it. Here the problem is that sto has to be of eltype Dual, not Float64. You can eg do `claytonsample(tau; randmat) = claytonsample!(similar(randmat, typeof(tau)), tau; randmat)` [untested] and differentiating this function should work (note the similar command constructs an array just like randmat, but of eltype typeof(tau))

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [August 31, 2021, 11:11am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/3 "2021-08-31T11:11:40Z")

</div>

I see. In this case, is there then a tradeoff between the speed of AD vs the speed gained by avoiding allocations (which is what I tried to do by writing `claytonsample!` that way)?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 31, 2021, 11:17am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/4 "2021-08-31T11:17:19Z")

</div>

You don’t have to choose. You can preallocate in functions you AD. See

> **[GitHub - SciML/PreallocationTools.jl: Tools for building non-allocating...](https://github.com/SciML/PreallocationTools.jl)**
>
> Tools for building non-allocating pre-cached functions in Julia, allowing for GC-free usage of automatic differentiation in complex codes - GitHub - SciML/PreallocationTools.jl: Tools for building ...

which has some helper functions and examples.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [August 31, 2021, 11:31am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/5 "2021-08-31T11:31:03Z")

</div>

If I understand correctly then this should work:

```julia
using ForwarDiff
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)

function claytonsample!(sto, τ; randmat=randmat)
    sto = get_tmp(sto, τ)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

```

yet

```julia
julia> ForwardDiff.derivative(τ -> claytonsample!(sto, τ), 0.3)
ERROR: MethodError: no method matching get_tmp(::Matrix{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#138#139", Float64}, Float64, 1})
Closest candidates are:
  get_tmp(::PreallocationTools.DiffCache, ::T) where T<:ForwardDiff.Dual at /Users/amrods/.julia/packages/PreallocationTools/l1Nsk/src/PreallocationTools.jl:26
  get_tmp(::PreallocationTools.DiffCache, ::Number) at /Users/amrods/.julia/packages/PreallocationTools/l1Nsk/src/PreallocationTools.jl:39

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 31, 2021, 12:36pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/6 "2021-08-31T12:36:04Z")

</div>

You never made a DiffCache.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [August 31, 2021, 1:16pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/7 "2021-08-31T13:16:06Z")

</div>

is this error because you have not defined broadcasting for `DiffCache`?

```julia
using ForwardDiff
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)
stod = dualcache(sto)

function claytonsample!(sto, τ; randmat=randmat)
    sto = get_tmp(sto, τ)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

ForwardDiff.derivative(τ -> claytonsample!(stod, τ), 0.3)

ERROR: LoadError: DimensionMismatch("array could not be broadcast to match destination")
Stacktrace:
  [1] check_broadcast_shape
    @ ./broadcast.jl:520 [inlined]
  [2] check_broadcast_axes
    @ ./broadcast.jl:523 [inlined]
  [3] instantiate
    @ ./broadcast.jl:269 [inlined]
  [4] materialize!
    @ ./broadcast.jl:894 [inlined]
  [5] materialize!
    @ ./broadcast.jl:891 [inlined]
  [6] claytonsample!(sto::PreallocationTools.DiffCache{Matrix{Float64}, Matrix{ForwardDiff.Dual{nothing, Float64, 12}}}, τ::ForwardDiff.Dual{ForwardDiff.Tag{var"#74#75", Float64}, Float64, 1}; randmat::Matrix{Float64})
    @ Main ~/test.jl:1307
  [7] claytonsample!(sto::PreallocationTools.DiffCache{Matrix{Float64}, Matrix{ForwardDiff.Dual{nothing, Float64, 12}}}, τ::ForwardDiff.Dual{ForwardDiff.Tag{var"#74#75", Float64}, Float64, 1})
    @ Main ~/test.jl:1306
  [8] (::var"#74#75")(τ::ForwardDiff.Dual{ForwardDiff.Tag{var"#74#75", Float64}, Float64, 1})
    @ Main ~/test.jl:1319
  [9] derivative(f::var"#74#75", x::Float64)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/5gUap/src/derivative.jl:14
 [10] top-level scope
    @ ~/test.jl:1319
in expression starting at /Users/amrods/test.jl:1319

```

explicit indexing works, but yields `NaN`:

```julia
using ForwardDiff
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)
stod = dualcache(sto)

function claytonsample!(sto, τ; randmat=randmat)
    sto = get_tmp(sto, τ)
    for i in eachindex(randmat)
        sto[i] = randmat[i]
    end
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

julia> ForwardDiff.derivative(τ -> claytonsample!(stod, τ), 0.3)
55×2 Matrix{Float64}:
 0.0 NaN
 0.0 NaN
 0.0 NaN
 0.0 NaN
 0.0 NaN
 ⋮    
 0.0 NaN
 0.0 NaN
 0.0 NaN
 0.0 NaN
 0.0 NaN

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 2, 2021, 12:31pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/9 "2021-09-02T12:31:31Z")

</div>

Well, yeah, you didn’t match the chunk sizes. If you’re differentiation has a single chunk you need a cache of a single chunk as well. Its default is to assume Jacobians w.r.t. a sized variable matching the cache, but if you’re differentiating w.r.t. a scalar you need to have a single chunk. So the following works:

```julia
using ForwardDiff
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)
stod = dualcache(sto, Val{1})

function claytonsample!(sto, τ; randmat=randmat)
    sto = get_tmp(sto, τ)
    @show size(sto), size(randmat)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

ForwardDiff.derivative(τ -> claytonsample!(stod, τ), 0.3)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 2, 2021, 12:39pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/10 "2021-09-02T12:39:59Z")

</div>

Full explanation is now in the docs: [https://github.com/SciML/PreallocationTools.jl/blob/master/README.md#handling-chunk-sizes](https://github.com/SciML/PreallocationTools.jl/blob/master/README.md#handling-chunk-sizes)

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 2, 2021, 11:53pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/11 "2021-09-02T23:53:40Z")

</div>

Notice that `stod = dualcache(sto, Val{1})` here, while in the docs you have `stod = dualcache(sto)`, which results in the broadcast error:

```julia
ERROR: LoadError: DimensionMismatch("array could not be broadcast to match destination")

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 2, 2021, 11:57pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/12 "2021-09-02T23:57:19Z")

</div>

Well yes, because you have to make sure you differentiate with the right chunk size.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 4, 2021, 1:53am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/13 "2021-09-04T01:53:54Z")

</div>

Is there any reason for not defining indexing of `DiffCache`? Look at this example, where I’m getting an error because `getindex` is not defined

```julia
using Optim
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)
stod = dualcache(sto, Val{2})

function claytonsample!(sto, τ; randmat=randmat)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

function claytonsample!(stod::T, τ; randmat=randmat) where {T <: PreallocationTools.DiffCache}
    stod = get_tmp(stod, τ)
    claytonsample!(stod, τ; randmat=randmat)
end

function corr(X)
    x = X[:, 1]
    y = X[:, 2]
    sum(x .* y)/length(x) - sum(x) * sum(y)/length(x)^2
end

function obj(sto, τ; randmat=randmat)
    τ′ = abs(τ) - 1
    claytonsample!(sto, τ′; randmat=randmat)
    corr(sto)
end

opt_ad = optimize(τ -> obj(stod, τ[1]), [5.0, 5.0], BFGS(); autodiff=:forward)
ERROR: LoadError: MethodError: no method matching getindex(::PreallocationTools.DiffCache{Matrix{Float64}, Matrix{ForwardDiff.Dual{nothing, Float64, 2}}}, ::Colon, ::Int64)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 4, 2021, 9:49am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/14 "2021-09-04T09:49:39Z")

</div>

> [@amrods](#):
>
> Is there any reason for not defining indexing of `DiffCache` ? Look at this example, where I’m getting an error because `getindex` is not defined
> 
> ```julia
> 
> ```

No. The documentation shows that you should never act on the DiffCache, always on the tmp that comes out of it. If you’re indexing the DiffCache you’re doing something wrong. No example in the documentation ever shows doing anything remotely like that.

I think you just confused yourself by writing some confusing code. Remember `sto` is a DiffCache, so you need to get its tmp:

```julia
function obj(sto, τ; randmat=randmat)
    τ′ = abs(τ) - 1
    claytonsample!(sto, τ′; randmat=randmat)
    corr(get_tmp(stod, τ))
end

```

Though I would really clean up that code first because it’s very implicit in confusing ways.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 4, 2021, 10:06am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/15 "2021-09-04T10:06:50Z")

</div>

I see. I wanted to have a method for `claytonsample!` when `sto` is a `DiffCache` and another one when `sto` is a `Matrix`. I have to do the same for `obj`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 21, 2021, 9:50pm UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/16 "2021-09-21T21:50:23Z")

</div>

Your simpler code

```julia
using ForwardDiff
using PreallocationTools

randmat = rand(10, 2)
sto = similar(randmat)
stod = dualcache(sto)

function claytonsample!(sto, τ; randmat=randmat)
    sto = get_tmp(sto, τ)
    sto .= randmat
    τ == 0 && return sto

    n = size(sto, 1)
    for i in 1:n
        v = sto[i, 2]
        u = sto[i, 1]
        sto[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
    end
    return sto
end

ForwardDiff.derivative(τ -> claytonsample!(stod, τ), 0.3)

```

now works on the latest version. The chunk size difference is now handled automatically.

---

<div class="post-metadata">

**Author:** ![amrods](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amrods/32/2543_2.png) [@amrods](https://discourse.julialang.org/u/amrods)\
**Post date:** [September 22, 2021, 5:00am UTC](https://discourse.julialang.org/t/how-to-make-this-function-compatible-with-forwarddiff/67415/17 "2021-09-22T05:00:04Z")

</div>

If I can just suggest something, it would be to keep more “transparent” functions available in the examples of the documentation. One of the things I struggled with when I started playing with the package, was understanding what it was doing. Examples with functions using `OrdinaryDiffEq` represent a hurdle for people who would want to use the `PreallocationTools` but not necessarily know or work with differential equations.