# Turing TypeErrror Issue

**URL:** https://discourse.julialang.org/t/turing-typeerrror-issue/76032
**Category:** Probabilistic Programming
**Created:** [February 8, 2022, 3:25pm UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032 "2022-02-08T15:25:10Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [February 8, 2022, 3:25pm UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032/1 "2022-02-08T15:25:10Z")

</div>

I’m experimenting with using Turing to handle some Gaussian random fields problems, and I’m struggling a bit with the interface. The following is supposed to do a Bayesian posterior sampling of the problem

y= u(x\_i) + \eta\_i

where u has a prior mean zero Gaussian with covariance (-d^2/dx^2)^{-1}. I’m trying to use some old code I wrote for FFT based sampling, which might be the problem. Here’s what I’m running:

```julia
using Turing
using Random
using Distributions
using FFTW

N = 32;
γ = 0.01;
x = LinRange(0,1,N+1)[2:end-1];
uᵗ(x) = x*(1-x);

Random.seed!(500); # set a seed for reproducibility
y_data = @. uᵗ(x) + γ * randn()

"""
`build_field` - Build a mean zero Gaussian random field with the (-Δ)^{-α} covariance operator in dimension one

### Fields
`ξ` - Vector of N(0,1) values  
`α=1` - Smoothness parameter
"""
function build_field(ξ; α=1)
    N = length(ξ)
    
    uhat = zeros(ComplexF64,2*N); # preallocate space

    # construct the eigenvalues
    k = 1:N;
    λ = @. 1/(π*k)^(2*α);

    # fill in the nonzero entries
    # NOTE we need to multiply by 2 *N for FFT scaling
    @. uhat[2:N+1] = 2 * N * sqrt(λ) * sqrt(2) * ξ;

    # invert and get the relevant imaginary part
    u = imag.(ifft(uhat))[N+2:end];
    return u
end

@model function mean_recovery(y_data)
    ξ ~ MvNormal(zeros(N), 1.)
    u = build_field(ξ);
    
    for i in 1:length(y_data)
       y_data[i]~Normal(u[i], γ^2)
    end
    
end

model=mean_recovery(y_data)

chain = sample(model, HMC(0.1, 10), 10^4)

```

The code then seems to run for a few moments before spitting out:

```julia
TypeError: in new, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 11}

Stacktrace:... # omitted for brevity, but I can put in if requested.

```

---

<div class="post-metadata">

### Author: ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [February 8, 2022, 8:13pm UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032/2 "2022-02-08T20:13:50Z")

</div>

> [@gideonsimpson](#):
>
> ` uhat = zeros(ComplexF64,2*N); # preallocate space`

You’ve constructed this to have a `ComplexF64` eltype, but this is incompatible with operator-overloading automatic differentiation (AD) packages like ForwardDiff, ReverseDiff, and Tracker, which use a custom numeric type to compute gradients. Instead, try

```julia
uhat = zeros(complex(eltype(ξ)),2*N);

```

The next issue you’ll have is that not all AD packages support complex numbers. Zygote definitely does, but it doesn’t support mutation. But if you want to differentiate an FFT, I think Zygote is your only option. So here’s a non-mutating `build_field`:

```julia
function build_field(ξ; α=one(eltype(ξ)))
    N = length(ξ)
    
    # construct the eigenvalues
    πk = π * (1:N);
    # NOTE we need to multiply by 2 *N for FFT scaling
    c = 2N * sqrt(2)
    umid = @. c * ξ / πk^α;
    uhat = [0; umid; zeros(N - 1)]

    # invert and get the relevant imaginary part
    u = @views imag.(ifft(uhat)[N+2:end]);
    return u
end

```

See [https://turing.ml/dev/docs/using-turing/autodiff](https://turing.ml/dev/docs/using-turing/autodiff) for details.

Also, from your simulated data, it seems `γ` might be intended to be a standard deviation? Note that `Normal` takes a standard deviation as a parameter, so perhaps you mean to use `Normal(u[i], γ) ` as your likelihood?

---

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [February 9, 2022, 1:39am UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032/3 "2022-02-09T01:39:50Z")

</div>

I tried your suggested version of `build_field`, but now this generates the error:

```julia
type DataType has no field mutable

Stacktrace:
  [1] getproperty
    @ ./Base.jl:37 [inlined]
  [2] adjoint
    @ ~/.julia/packages/Zygote/zowrf/src/lib/lib.jl:281 [inlined]
  [3] _pullback
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined]
# keeps going for a bit

```

But the switch to `Zygote` does appear to have addressed the other issue.

Also, are any of the interpolation packages compatible with `Zygote` for automatic differentiation?

---

<div class="post-metadata">

### Author: ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)
#### Post date: [February 9, 2022, 8:53am UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032/4 "2022-02-09T08:53:35Z")

</div>

> [@gideonsimpson](#):
>
> I tried your suggested version of `build_field` , but now this generates the error:

> **Odd, this runs fine for me:**
>
> ```julia
> using Turing, Random, Distributions, FFTW, Zygote
> Turing.setadbackend(:zygote)
> 
> N = 32;
> γ = 0.01;
> x = LinRange(0,1,N+1)[2:end-1];
> uᵗ(x) = x*(1-x);
> 
> Random.seed!(500); # set a seed for reproducibility
> y_data = @. uᵗ(x) + γ * randn()
> 
> """
> `build_field` - Build a mean zero Gaussian random field with the (-Δ)^{-α} covariance operator in dimension one
> 
> ### Fields
> `ξ` - Vector of N(0,1) values  
> `α=1` - Smoothness parameter
> """
> function build_field(ξ; α=one(eltype(ξ)))
> N = length(ξ)
>     
> # construct the eigenvalues
> πk = π * (1:N);
> # NOTE we need to multiply by 2 *N for FFT scaling
> c = 2N * sqrt(2)
> umid = @. c * ξ / πk^α;
> uhat = [0; umid; zeros(N - 1)]
> 
> # invert and get the relevant imaginary part
> u = @views imag.(ifft(uhat)[N+2:end]);
> return u
> end
> 
> @model function mean_recovery(y_data)
> ξ ~ filldist(Normal(), N)
> u = build_field(ξ);
> y_data .~ Normal.(u, γ)    
> end
> 
> model=mean_recovery(y_data)
> 
> chain = sample(model, NUTS(), MCMCThreads(), 10^3, 4)
> 
> ```
> 
> Note that I switched to NUTS and dropped the number of requested samples. This model samples so well with NUTS that the 4000 draws for most parameters have effective sample sizes of over 10,000.

> [@gideonsimpson](#):
>
> Also, are any of the interpolation packages compatible with `Zygote` for automatic differentiation?

I’m not familiar with the interpolation packages. I’d say just try and see.

---

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [February 9, 2022, 7:26pm UTC](https://discourse.julialang.org/t/turing-typeerrror-issue/76032/5 "2022-02-09T19:26:29Z")

</div>

Never mind. I think one of my packages was out of date. `Interpolations` appears to work.
