# MethodError: no method matching Float64 when using ForwardDiff.jl despite me pre-allocating arrays to take floating point

``````using Optimization, OptimizationOptimJL, OptimizationBBO, OptimizationMOI, ForwardDiff, ModelingToolkit, Random, OptimizationOptimisers, Statistics, Plots

function MSEg(u, p)
X = p.X
y = p.y
h = u[1]
beta = u[2:end]
N = length(y)
Xbeta = X*beta

krnl = u -> exp.(-(u .* u) ./ 2) ./ sqrt(2*π)

yhat = Array{Float64, 2}(undef, N, 1)
weight_i = Array{Float64, 2}(undef, N, 1)
weight_sum = Array{Float64, 2}(undef, N, 1)

for i = 1:N
u = (Xbeta .- Xbeta[i]) ./ h
Ku = krnl(u)
weight_i[i] = Ku[i]
weight_sum[i] = sum(Ku)
yhat[i] = sum(Ku .* y) ./ sum(Ku)
end

residual = (y .- yhat) ./ (1 .- (weight_i ./ weight_sum))
SBeta = sum(residual.^2)
end

Random.seed!(1)
N = 1000
X = -1 .+ 2 .* rand(N, 2)
beta = [1;3]
Xbeta = X*beta
y = Xbeta.^2 .* (1 .+ sin.(Xbeta)) .+ randn(N, 1)

y_ = y .- X[:,1]
X_ = X[:,2]
beta0 = X_'*X_\X_'*y_
beta0 = [1; beta0]

x_ = X*beta0
hx = median(abs.(x_ .- median(x_)))/0.6745*(4/3/N)^0.2
hy = median(abs.(y .- median(y)))/0.6745*(4/3/N)^0.2
h = sqrt(hy*hx)

param0 = [h;beta0]
u = param0

struct smol
X
y
end

p = smol(X, y)

cons(res, u, p) = (res .= u[2])

lb = [0;-Inf;-Inf]
ub = [Inf;Inf;Inf]

optf = OptimizationFunction(MSEg, Optimization.AutoForwardDiff(), cons = cons)
prob = OptimizationProblem(optf, u, p, lb = lb, ub = ub, lcons = 1, ucons = 1)
``````

I had received this error in the past when I tried to use a for loop to populate elements of arrays and matrices. I had not had trouble in the past with this error in a while but now it is coming back again and I cannot understand why

despite me pre-allocating arrays to take floating point

This is the problem. You need:

``````function MSEg(u::AbstractVector{T}, p) where {T}
# ...
yhat = Array{T, 2}(undef, N, 1)
#...
``````

Your function must not assume that `u` is a `Vector{Float64}`.

I am still confused. I changed the function below and get the same error when running the script

``````function MSEg(u::AbstractVector{T}, p) where {T}
typeof(u)
X = p.X
y = p.y
h = u[1]
beta = u[2:end]
N = length(y)
Xbeta = X*beta

krnl = u -> exp.(-(u .* u) ./ 2) ./ sqrt(2*π)

# yhat = Array{Float64, 2}(undef, N, 1)
# weight_i = Array{Float64, 2}(undef, N, 1)
# weight_sum = Array{Float64, 2}(undef, N, 1)

yhat = Array{T, 2}(undef, N, 1)
weight_i = Array{T, 2}(undef, N, 1)
weight_sum = Array{T, 2}(undef, N, 1)

for i = 1:N
u = (Xbeta .- Xbeta[i]) ./ h
Ku = krnl(u)
weight_i[i] = Ku[i]
weight_sum[i] = sum(Ku)
yhat[i] = sum(Ku .* y) ./ sum(Ku)
end
residual = (y .- yhat) ./ (1 .- (weight_i ./ weight_sum))
SBeta = sum(residual.^2)
end
``````

These won’t be Float64, so the type is incorrect.

You can simplify your `MSEg` code a lot (and there were a bunch of other problems in your code, like the wrong optimizer). How about this:

``````using Optimization
import OptimizationMOI
import ForwardDiff
import Ipopt
import Random
using Statistics

kernel(u) = exp(-(u * u) / 2) / sqrt(2π)

function MSEg(u::AbstractArray{T}, p) where {T}
X, y = p.X, p.y
h, beta = u[1], u[2:end]
Xbeta = X * beta
residuals = zero(T)
for i in 1:length(y)
Ku = kernel.((Xbeta .- Xbeta[i]) ./ h)
yhat = sum(Ku .* y) ./ sum(Ku)
residuals += ((y[i] - yhat) / (1 - (Ku[i] / sum(Ku))))^2
end
return residuals
end

Random.seed!(1)
N = 1000
X = -1 .+ 2 .* rand(N, 2)
beta = [1;3]
Xbeta = X*beta
y = Xbeta.^2 .* (1 .+ sin.(Xbeta)) .+ randn(N, 1)

y_ = y .- X[:,1]
X_ = X[:,2]
beta0 = X_' * X_ \ X_' * y_
beta0 = [1; beta0]

x_ = X * beta0
hx = median(abs.(x_ .- median(x_))) / 0.6745 * (4 / 3 / N)^0.2
hy = median(abs.(y .- median(y))) / 0.6745 * (4 / 3  / N)^0.2
h = sqrt(hy*hx)

struct smol
X
y
end

p = smol(X, y)

cons(res, u, p) = (res .= u[2])

optf = OptimizationFunction(MSEg, Optimization.AutoForwardDiff(), cons = cons)
prob = OptimizationProblem(
optf,
[h; beta0],
p;
lb = [0, -Inf, -Inf],
ub = [Inf, Inf, Inf],
lcons = [1.0],
ucons = [1.0],
)
sol = solve(prob, Ipopt.Optimizer())
``````