Dual{::ForwardDiff...} type error when JuMP tries autodiff User-defined Function

Specific Domains Optimization (Mathematical)
1

Hi, I am trying to optimize User-defined function with and I got this error:

ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##84#86")){typeof(min_f)},Float64},Float64,11})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:741
  Float64(::Int8) at float.jl:60
  ...
Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##84#86")){typeof(min_f)},Float64},Float64,11}) at ./number.jl:7
 [2] setindex! at ./array.jl:769 [inlined]
 [3] MC36SM_Mindaugo_SS(::Float64, ::NTuple{32,ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##84#86")){typeof(min_f)},Float64},Float64,11}}, ::Array{Float64,1}, ::Float64, ::Float64) at /home/kest/Projects/gjJulia/testMC36SM.jl:90

I figured out that error is because of trying to do (in User-defined function for JuMP):

g = zeros(36, 4);
.....
g[i,j] = Go[1]

g type is 2D Array of Float64, but Go[1] (optimization variable) passed by JuMP autodiff type is:

Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##84#86")){typeof(min_f)},Float64}}(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0)

Maybe somebody have suggestion how to solve such type mismatch?

2

Please read: Please read: make it easier to help you, particularly about including a minimum working example.

For your error, see ForwardDiff: Using on functions that are defined on floats64

Basically, you need things typed as Real instead of Float64.

3 Likes
3

Thanks @odow, error is solved. Sorry not including MWE, I just taught it is simple thing and would be shorter to read without MWE and I didn’t find that post explaining ForwardDiff uses on Float64. Now I see that is mentioned in ForwardDiff documentation, I am just starting to use Julia and didn’t understand that after reading JuMP User-defined function’s doc. Thanks again.

4

I have made a pull request to improve JuMP’s documentation: https://github.com/JuliaOpt/JuMP.jl/pull/2000

5

I noticed that if I use Real instead of Float calculation is ~400 fold slower. Also there are large increase in allocations. But only difference in codes is that I use Real instead of Float (by default) in few Arrays.

  33.443 ms (2354 allocations: 2.15 MiB)
  13.397 s (616406763 allocations: 9.19 GiB)

also

julia> @btime z = zeros(Real,36,36);
  7.131 ÎĽs (1 allocation: 10.25 KiB)

julia> @btime z = zeros(36,36);
  1.613 ÎĽs (1 allocation: 10.25 KiB)

MWE of our’s model demonstrating Real vs Float benchmark:

using LinearAlgebra, BenchmarkTools

function MC36SM_Mindaugo_SS_Real(vj::Float64, par_tuple, p::Array{Real}, pc1c2::Float64, pc2c1::Float64)
    PPP = zeros(Real, 36, 36);
    p1 = 0;
    p2 = 0;
    p3 = 0;
    p4 = 0;

    pr = collect(par_tuple);
    par = reshape(pr, 4, 8);

    A =  par[:,1];
    v0 = par[:,2];
    Go = par[:,3];
    Gc = par[:,4];
    Ro = par[:,5];
    Rc = par[:,6];
    Pt = par[:,7];
    P = par[:,8];

    K = Pt[1];
    g = zeros(Real, 36, 4);
    gg = zeros(Real, 36); #buffer for the sums
    v = zeros(Real, 36, 4);
    k = zeros(Real, 36, 4);
    R = zeros(Real, 36, 4);

   # //states conductances
    for i1 = 1:2
        for i2 = 1:3
            for i3 = 1:3
                for i4 = 1:2
                    if (i1 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Go[1];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Ro[1];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Gc[1];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Rc[1];
                    end
                    if (i2 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Go[2];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Ro[2];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Gc[2];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Rc[2];
                    end
                    if (i3 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Go[3];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Ro[3];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Gc[3];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Rc[3];
                    end
                    if (i4 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Go[4];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Ro[4];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Gc[4];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Rc[4];
                    end
                end
            end
        end
    end

    for a = 1:4 #// voltages and rectification
    #    gg = 1 ./ sum( (1 ./ g), dims = 2 ); #sums rows
        fill!(gg, 0.0)
       	row_sum!(gg, g) #happens twice so took it out into function
    #    v=vj*[gg gg gg gg]./g;
        for j = 1:4, i = 1:36 #loop so no allocations
            @inbounds v[i, j] = vj * gg[i] / g[i, j];
        end
        g = g .* exp.(v ./ R);
    end

    for i = 1:36
        for j = 1:4
            @inbounds k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
        end
    end

   # //------------Matrica P-------------------------------------------------
    for i1 = 1:2 for i2 = 1:3 for i3 = 1:3 for i4 = 1:2
                    for j1 = 1:2 for j2 = 1:3 for j3 = 1:3 for j4 = 1:2
                                    i = (i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1;
                                    j = (j1 - 1) * 18 + (j2 - 1) * 6 + (j3 - 1) * 2 + (j4 - 1) * 1 + 1;


      # %-------be Pt---------------------
      # % p1
                                    if (i1 == 1)
                                        if (j1 == 1)
                                            p1 = (1 - K * k[i,1] / (1 + k[i,1]));
                                        else
                                            p1 = K * k[i,1] / (1 + k[i,1]);
                                        end
                                    else
                                        if (j1 == 2)
                                            p1 = (1 - K / (1 + k[i,1]));
                                        else
                                            p1 = K / (1 + k[i,1]);
                                        end
                                    end
   # %------------------------------
   #   % p2
                                    if (i2 == 1)
                                        if (j2 == 1)
                                            p2 = (1 - K * k[i,2] / (1 + k[i,2]));
                                        end
                                        if (j2 == 2)
                                            p2 = K * k[i,2] / (1 + k[i,2]);
                                        end
                                        if (j2 == 3)
                                            p2 = 0;
                                        end
                                    end

                                    if (i2 == 2)
                                        if (j2 == 1)
                                            p2 = (K / (1 + k[i,2])) * (1 - pc1c2);
                                        end
                                        if (j2 == 2)
                                            p2 = (1 - K / (1 + k[i,2])) * (1 - pc1c2);
                                        end
                                        if (j2 == 3)
                                            p2 = pc1c2;
                                        end
                                    end

                                    if (i2 == 3)
                                        if (j2 == 1)
                                            p2 = 0;
                                        end
                                        if (j2 == 2)
                                            p2 = pc2c1;
                                        end
                                        if (j2 == 3)
                                            p2 = 1 - pc2c1;
                                        end
                                    end
     # %----------------------------------------
     # % p3
                                    if (i3 == 1)
                                        if (j3 == 1)
                                            p3 = (1 - K * k[i,3] / (1 + k[i,3]));
                                        end
                                        if (j3 == 2)
                                            p3 = K * k[i,3] / (1 + k[i,3]);
                                        end
                                        if (j3 == 3)
                                            p3 = 0;
                                        end
                                    end

                                    if (i3 == 2)
                                        if (j3 == 1)
                                            p3 = (K / (1 + k[i,3])) * (1 - pc1c2);
                                        end
                                        if (j3 == 2)
                                            p3 = (1 - K / (1 + k[i,3])) * (1 - pc1c2);
                                        end
                                        if (j3 == 3)
                                            p3 = pc1c2;
                                        end
                                    end

                                    if (i3 == 3)
                                        if (j3 == 1)
                                            p3 = 0;
                                        end
                                        if (j3 == 2)
                                            p3 = pc2c1;
                                        end
                                        if (j3 == 3)
                                            p3 = 1 - pc2c1;
                                        end
                                    end
     # %----------------------------
     # %p4
                                    if (i4 == 1)
                                        if (j4 == 1)
                                            p4 = (1 - K * k[i,4] / (1 + k[i,4]));
                                        else
                                            p4 = K * k[i,4] / (1 + k[i,4]);
                                        end
                                    else
                                        if (j4 == 2)
                                            p4 = (1 - K / (1 + k[i,4]));
                                        else
                                            p4 = K / (1 + k[i,4]);
                                        end
                                    end

                                    PPP[j,i] = p1 * p2 * p3 * p4;

                                end
                            end
                        end
                    end
                end
            end
        end
    end

    d_g = 100000;
    g_old = 100000;
    g_final = 0;
    i = 0;
    tmp = similar(p);
    # ggg = zeros(MVector{36,Float64});
#    ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
    fill!(gg, 0.0);
    row_sum!(gg, g)
    while (d_g / (g_old + g_final) > .001e-5)
        i = i + 1;
    #    p=p*PPP;
        mul!(tmp, PPP, p)
        p, tmp = tmp, p;
        g_final = dot(gg, p); # (p*ggg)[1];
        d_g = abs(g_old - g_final);
        g_old = g_final;
    #    print(" ",i,", ", d_g / (g_old + g_final), " ")
    end
    # print(i, " ");
    # println(i);
    return g_final, i;
end

function MC36SM_Mindaugo_SS(vj::Float64, par_tuple, p::Array{Float64}, pc1c2::Float64, pc2c1::Float64)
    PPP = zeros(36, 36);
    p1 = 0;
    p2 = 0;
    p3 = 0;
    p4 = 0;

    pr = collect(par_tuple);
    par = reshape(pr, 4, 8); #[]... - Tuple to array, and reshape into dimensions

    A =  par[:,1];
    v0 = par[:,2];
    Go = par[:,3];
    Gc = par[:,4];
    Ro = par[:,5];
    Rc = par[:,6];
    Pt = par[:,7];
    P = par[:,8];

    K = Pt[1];
    g = zeros(36, 4);
    gg = zeros(36); #buffer for the sums
    v = zeros(36, 4);
    k = zeros(36, 4);
    R = zeros(36, 4);

   # //states conductances
    for i1 = 1:2
        for i2 = 1:3
            for i3 = 1:3
                for i4 = 1:2
                    if (i1 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Go[1];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Ro[1];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Gc[1];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Rc[1];
                    end
                    if (i2 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Go[2];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Ro[2];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Gc[2];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Rc[2];
                    end
                    if (i3 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Go[3];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Ro[3];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Gc[3];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Rc[3];
                    end
                    if (i4 == 1)
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Go[4];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Ro[4];
                    else
                        g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Gc[4];
                        R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Rc[4];
                    end
                end
            end
        end
    end

    for a = 1:4 #// voltages and rectification
    #    gg = 1 ./ sum( (1 ./ g), dims = 2 ); #sums rows
        fill!(gg, 0.0)
       	row_sum!(gg, g) #happens twice so took it out into function
    #    v=vj*[gg gg gg gg]./g;
        for j = 1:4, i = 1:36 #loop so no allocations
            @inbounds v[i, j] = vj * gg[i] / g[i, j];
        end
        g = g .* exp.(v ./ R);
    end

    for i = 1:36
        for j = 1:4
            @inbounds k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
        end
    end

   # //------------Matrica P-------------------------------------------------
    for i1 = 1:2 for i2 = 1:3 for i3 = 1:3 for i4 = 1:2
                    for j1 = 1:2 for j2 = 1:3 for j3 = 1:3 for j4 = 1:2
                                    i = (i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1;
                                    j = (j1 - 1) * 18 + (j2 - 1) * 6 + (j3 - 1) * 2 + (j4 - 1) * 1 + 1;


      # %-------be Pt---------------------
      # % p1
                                    if (i1 == 1)
                                        if (j1 == 1)
                                            p1 = (1 - K * k[i,1] / (1 + k[i,1]));
                                        else
                                            p1 = K * k[i,1] / (1 + k[i,1]);
                                        end
                                    else
                                        if (j1 == 2)
                                            p1 = (1 - K / (1 + k[i,1]));
                                        else
                                            p1 = K / (1 + k[i,1]);
                                        end
                                    end
   # %------------------------------
   #   % p2
                                    if (i2 == 1)
                                        if (j2 == 1)
                                            p2 = (1 - K * k[i,2] / (1 + k[i,2]));
                                        end
                                        if (j2 == 2)
                                            p2 = K * k[i,2] / (1 + k[i,2]);
                                        end
                                        if (j2 == 3)
                                            p2 = 0;
                                        end
                                    end

                                    if (i2 == 2)
                                        if (j2 == 1)
                                            p2 = (K / (1 + k[i,2])) * (1 - pc1c2);
                                        end
                                        if (j2 == 2)
                                            p2 = (1 - K / (1 + k[i,2])) * (1 - pc1c2);
                                        end
                                        if (j2 == 3)
                                            p2 = pc1c2;
                                        end
                                    end

                                    if (i2 == 3)
                                        if (j2 == 1)
                                            p2 = 0;
                                        end
                                        if (j2 == 2)
                                            p2 = pc2c1;
                                        end
                                        if (j2 == 3)
                                            p2 = 1 - pc2c1;
                                        end
                                    end
     # %----------------------------------------
     # % p3
                                    if (i3 == 1)
                                        if (j3 == 1)
                                            p3 = (1 - K * k[i,3] / (1 + k[i,3]));
                                        end
                                        if (j3 == 2)
                                            p3 = K * k[i,3] / (1 + k[i,3]);
                                        end
                                        if (j3 == 3)
                                            p3 = 0;
                                        end
                                    end

                                    if (i3 == 2)
                                        if (j3 == 1)
                                            p3 = (K / (1 + k[i,3])) * (1 - pc1c2);
                                        end
                                        if (j3 == 2)
                                            p3 = (1 - K / (1 + k[i,3])) * (1 - pc1c2);
                                        end
                                        if (j3 == 3)
                                            p3 = pc1c2;
                                        end
                                    end

                                    if (i3 == 3)
                                        if (j3 == 1)
                                            p3 = 0;
                                        end
                                        if (j3 == 2)
                                            p3 = pc2c1;
                                        end
                                        if (j3 == 3)
                                            p3 = 1 - pc2c1;
                                        end
                                    end
     # %----------------------------
     # %p4
                                    if (i4 == 1)
                                        if (j4 == 1)
                                            p4 = (1 - K * k[i,4] / (1 + k[i,4]));
                                        else
                                            p4 = K * k[i,4] / (1 + k[i,4]);
                                        end
                                    else
                                        if (j4 == 2)
                                            p4 = (1 - K / (1 + k[i,4]));
                                        else
                                            p4 = K / (1 + k[i,4]);
                                        end
                                    end

                                    PPP[j,i] = p1 * p2 * p3 * p4;

                                end
                            end
                        end
                    end
                end
            end
        end
    end

    d_g = 100000;
    g_old = 100000;
    g_final = 0;
    i = 0;
    tmp = similar(p);
    # ggg = zeros(MVector{36,Float64});
#    ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
    fill!(gg, 0.0);
    row_sum!(gg, g)
    while (d_g / (g_old + g_final) > .001e-5)
        i = i + 1;
    #    p=p*PPP;
        mul!(tmp, PPP, p)
        p, tmp = tmp, p;
        g_final = dot(gg, p); # (p*ggg)[1];
        d_g = abs(g_old - g_final);
        g_old = g_final;
    end
    return g_final, i;
end

@inline function row_sum!(buffer, g) #fill the buffer in a loop -> no allocations
    for j = 1:size(g, 2), i = 1:size(g, 1)
        @inbounds buffer[i] += 1 / g[i, j];
    end
    buffer .= 1 ./ buffer;
end

# end

function getGj_Real(N, vj, par)
    gj = zeros(Real, N);
    pc1c2 = 0;
    pc2c1 = 0;
    pc1c2 = 0.01;
    pc2c1 = 0.001;
    ppp = zeros(Real, 36);
    ppp[1] = 1;
    iter = 0;               
    for i = 1:N
        gj[i], .. = MC36SM_Mindaugo_SS_Real(vj[i], par, ppp, pc1c2, pc2c1);
    end
    return gj
end

function getGj_Float(N, vj, par)
    gj = zeros(N);
    pc1c2 = 0;
    pc2c1 = 0;
    pc1c2 = 0.01;
    pc2c1 = 0.001;
    ppp = zeros(36);
    ppp[1] = 1;
    iter = 0;               
    for i = 1:N
        gj[i], .. = MC36SM_Mindaugo_SS(vj[i], par, ppp, pc1c2, pc2c1);
    end
    return gj
end

parAllFast = [ .15 20 20e-9 2e-9  1e10  1e10  .5 -1 
               .15 20 20e-9 2e-9  1e10  1e10  .5  1 ];
parAllSlow = [ .1  30 40e-9 1e-15 1e10  1e10  .5 -1
               .1  30 40e-9 1e-15 1e10  1e10  .5  1 ];
par = [parAllFast; parAllSlow];
N = 100;
vj_e = collect(LinRange(0, 120, N));

@btime getGj_Float(N, vj_e, par);
@btime getGj_Real(N, vj_e, par);

Maybe someone could comment about such problem?

6

The problem is that typing the elements of the vector as Real is not a concrete type.

Instead of

function foo(x)
    y = Vector{Real}(undef, 2)
    y[1] = x
    y[2] = 2 * x
    return sum(y)
end

Type the elements of the vector based on the input arguments

function foo(x::T) where {T <: Real}
    y = Vector{T}(undef, 2)
    y[1] = x
    y[2] = 2 * x
    return sum(y)
end
2 Likes
7

Thanks, that indeed solved performance problem.

8

Hi, I have a ion channel model (of cardiac physiology). I try to use JuMP for curve fitting and solving model parameters from experimental data. But I got similar error as mentioned before in this thread:

ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("#min_f#18")){StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1},Scene,getfield(Main, Symbol("#my_model#17")){Array{Float64,1}}},Float64},Float64,1})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:718
  Float64(::Int8) at float.jl:60

I_ext = zeros(25000) is global array in main function (not shown there), that is difference from previous posts in this thread.
Error appears on line u, v, w, I_ion_all[i] = fenton(u, v, w, I_ext[i];... in following file (function to calculate all data points):

function runFenton(I_ext::Array{T,1};
                    tv1m = 19.6,
                    tv2m = 1250.,
                    twp = 870.0,
                    twm = 41.0, #original
                    td = 0.25,
                    to = 12.5,
                    tr = 33.33,
                    tsi = 30.0,
                    xk = 10.0,
                    ucsi = 0.85,
                    uc = 0.13,
                    uv = 0.04
                    ) where {T <: Real}
u_all = zeros(T,length(I_ext))
I_ion_all = zeros(T,length(I_ext))
u::T = 0.; v::T = 0.; w::T = 0.;
for i = 1:length(I_ext)
u, v, w, I_ion_all[i] = fenton(u, v, w, I_ext[i];
                                tv1m = tv1m,
                                tv2m = tv2m,
                                twp = twp,
                                twm = twm,
                                td = td,
                                to = to,
                                tr = tr,
                                tsi = tsi,
                                xk = xk,
                                ucsi = ucsi,
                                uc = uc,
                                uv = uv)
u_all[i] = u
# print(u, " ")
end
return u_all, I_ion_all
end

Function to calculate one iteration:

function fenton(u, v, w, I_ext;
    tvp = 3.33,
    tv1m = 19.6,
    tv2m = 1250.,
    twp = 870.0,
    twm = 41.0, #original
    td = 0.25,
    to = 12.5,
    tr = 33.33,
    tsi = 30.0,
    xk = 10.0,
    ucsi = 0.85,
    uc = 0.13,
    uv = 0.04
    )

    Vfi = 15e-3#nernst potential of fast gate
    V0 = -85e-3#resting membrane potential

# numerical parameters
    timestep = 0.02 # size of time step in ms

#currents
    jfi = 0
    jso = 0
    jsi = 0

#calculation
    p = 0; # heaviside functions
    q = 0;

    if(u >= uc)
        p = 1;
    end

    if(u >= uv)
        q = 1.;
    end

    dv = (1 - p) * (1 - v) / ((1 - q) * tv1m + tv2m * q) - p * v / tvp;
    dw = (1 - p) * (1 - w) / twm - p * w / twp;

    jfi = -v * p * (u - uc) * (1 - u) / td;
    jso = u * (1 - p) / to + p / tr;
    jsi = -w * (1 + tanh(xk * (u - ucsi))) / (2. * tsi);

    u = u - timestep * (jfi + jso + jsi + I_ext)
    v = v + timestep * dv; # solving/updating v and w
    w = w + timestep * dw;

    return u, v, w, jfi + jso + jsi
end

I’ve tried many parameter combinations with where {T <: Real} but no success. It works only if I define I_ext inside runFenton function and initialize with zeros(T,some_integer), is it not possible to define it outside function and access with global scope or pass as parameter?
I_ext is constant during all calculation. I am optimizing parameters which I choose using runFenton() Keyword Arguments. I did not post all code for not making post very long.
Maybe somebody have any insights why that error appears?

9

Nice start on getting a minimum working example. But there are a few things you can do to make it better. (You may want to take a read of https://discourse.julialang.org/t/psa-make-it-easier-to-help-you again.)

  1. Simplify even more! You have lots of parameters. Try replacing them with 1.0, and if still errors, remove the parameter.
  2. Make the code fully reproducible. That means including all the using Package lines. You mention JuMP, but also ForwardDiff. It isn’t obvious how you are using them.
  3. Put everything in one script. I should be able to copy-and-paste your code into a REPL and have it trigger the error. I shouldn’t need to grab lots of different functions.
  4. Post the full error, even if you don’t think it’s helpful. Often it has some useful lines in there for people who know where to look.

Once you’ve simplified the problem to address the above points, we can figure out what’s going on.

10

Perhaps a starting point would be to examine the types of elements you are wanting to differentiate over. Julia won’t allow you to take a derivative over something you assert as a Float64 because it has a special number type to deal with derivatives. I would suggest printing out the types of each element within your code before each call to ForwardDiff is made. Perhaps, you are accidentally asserting the type of one of elements you are differentiating over as a Float64. Below is an example of what I mean.

julia> function f(x::Array{T}) where {T<:Real}
       println(eltype(x))
       y =  x[1]^2 + x[2]^2
       end
f (generic function with 1 method)

julia> ForwardDiff.gradient!(Array{Float64}(undef, 1, 2), x -> f(x), [3.0, 4.0])

#this is the type of x
ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##21#22")),Float64},Float64,2}

#this is the result 
1Ă—2 Array{Float64,2}:
 6.0  8.0

julia> function f_bad(x::Array{Float64})
       println(eltype(x))
       y =  x[1]^2 + x[2]^2
       end
f_bad (generic function with 3 methods)

julia> ForwardDiff.gradient!(Array{Float64}(undef, 1, 2), x -> f_bad(x), [3.0, 4.0])

#this is the type of x
ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##23#24")),Float64},Float64,2}

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##23#24")),Float64},Float64,2})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:718
  Float64(::Int8) at float.jl:60

The reason for this failure is that

julia> ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##23#24")),Float64},Float64,2} <: Float64
false

julia> ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##23#24")),Float64},Float64,2} <: Real
true

(Stacktrace omitted)

Note that even though you are passing [3.0, 4.0] in both cases, the element type of x is internally being converted to a ForwardDiff.Dual type because you are calling ForwardDiff over x.

1 Like
11

@bashonubuntu, thanks for valuable insight, situation is getting more clear now, yes I am using Float64 array I_ext[] which is defined outside objective function, if I define it in objective function using I_ext=zeros(T,some_number) and where {T<:Real} in function declaration, then it works without errors, but I want to pass I_ext[] (which is not changing during calculation) from the outside of the objective function, should I change type of that array? What type it should be then?
Objective function is using variable u for calculation of fit value. I_ext is used for calculation of u:

Actually there are few functions nested in objective function, but I think that doesn’t change the whole principle.

12

It worked. Actually I just added ::T before td (which is optimization parameter, = 0.25 is only for convenience in case I choose another parameter to optimize):

function runFenton(I_ext::Array{Float64,1};...,td::T = 0.25,...) where {T <: Real}
...
u_all = zeros(T,length(I_ext)) 
...
end

I think key line is u_all = zeros(T,length(I_ext)) since u_all defines return values.

But if I want to use td argument definition without :T, I got error LoadError: UndefVarError: T not defined at line with …zeros(…):

function runFenton(I_ext::Array{Float64,1};...,td = 0.25,...) where {T <: Real}
...
u_all = zeros(T,length(I_ext)) 
...
end

My question, is it possible to omit ::T before td, but still have T <: Real declaration somehow to be possible to use T in u_all = zeros(T,length(I_ext))?

13

Why define the parameter outside the function though? If you want to keep it fixed in the gradient calculation, you can do so and ignore the value of the gradient for that element (which will be 0.0). Please see an example below.

julia> function f(x::Array{T}) where {T<:Real}
       println(eltype(x))
       y =  x[1]^2 + x[2]^2
       end
f (generic function with 1 method)

julia> ForwardDiff.gradient!(Array{Float64}(undef, 1, 3), x -> f(x), [3.0, 4.0, 4.5])
ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##11#12")),Float64},Float64,3}
1Ă—3 Array{Float64,2}:
 6.0  8.0  0.0
14

Sorry, in previous post I forgot end in function definition, I corrected it. So u_all = zeros(T,length(I_ext)) is inside the function.

15

I am still a bit confused as er what to do. I have code in which I initialize a few matrices as, say,

A = zeros(10,10)

Is that where my problem is? Am I just to convert it to Reals? But why would

A = zeros(Real, 10, 10);

work?

16

Is that where my problem is?

Probably.

But why would A = zeros(Real, 10, 10); work?

Because now it can store a Float64 and a Dual number.

17

But do I just need the inputs of the objective function to be real, or every single element used in the function itself. For example, what follows still gives me the same error:

using JuMP, Ipopt;
function objfunc(x1::Real,x2::Real,x3::Real)
    return sum(SSresid(mod_pars,x1,x2,x3).^2);
end
model = Model(with_optimizer(Ipopt.Optimizer, max_iter=100, print_level=0));
#adding variables
@variable(model, 0.8 <= x1 <= 1.2);
@variable(model, 95 <= x2 <= 100);
@variable(model, 1e-6 <= x3 <= 0.01);
#register function
register(model, :objfunc, 3, objfunc, autodiff=true);
@NLobjective(model, Min, objfunc(x1,x2,x3));
@time xsol_jump = JuMP.optimize!(model);
18

SSresid also needs to be able to handle Dual.

I suggest you formulate a minimum working example (I don’t know what SSresid or mod_pars are), and post as a new question, rather than repurposing this old issue.

19

@danicaratelli Have you found solution to this problem? If yes, please share for our learning. Thanks.

20

@Humphrey_Lee I have switched to a different solver, in part because I could not figure this out quickly enough. I now use NLsolve.