Allow non-vectorized model functions?

[kzapfe]

It seems it cannot. When I use a model that multiplies a parameter within itself, a no method error appears.

model(x,p)=(1/(sqrt(2*pi)*p[1])) * exp(-( (x-p[2])/(sqrt(2)*p[1])) ^2)
xdata=collect(cuac[1][2:end]) #my bins
ydata=cuac[2] ;    #my histograms

And I get this beautiful error

LoadError: MethodError: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})
hile loading In[104], in expression starting on line 1
 in power_by_squaring at intfuncs.jl:78
 in model at In[103]:1
 in f at /home/karel/.julia/v0.4/LsqFit/src/curve_fit.jl:38
 in levenberg_marquardt at /home/karel/.julia/v0.4/Optim/src/levenberg_marquardt.jl:38
 in lmfit at /home/karel/.julia/v0.4/LsqFit/src/curve_fit.jl:29
 in curve_fit at /home/karel/.julia/v0.4/LsqFit/src/curve_fit.jl:39

I have already tried to write the function in different ways, but it seems that somewhere, inside the fit, the routine tries to multiply the parameters without knowing what are they.

My Julia is 0.4.3 and LsqFit is the last one avaible.

[mlubin]

Here's a reduced example that might reveal what's going on:

julia> x = ones(2)
2-element Array{Float64,1}:
 1.0
 1.0

julia> x^2
ERROR: MethodError: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})
...

[kzapfe]

Yes, exactly that is it, but it is going internally on the curve_fit routine and I seem to not be able to make it go away. I actually have a simpler example (adjusting a degree 2 polynomial) and it seems to trigger the same problem. Maybe I am not understanding the syntaxis right...

 model(x,p)=p[3]+p[2]*x+p[1]*x*x
    ydata=map(x->log(x),yys)
    rara=curve_fit(model,xxs,ydata, [0,0,0])

returns the same error.

[davidssmith]

@mlubin was just being nice. He was trying to tell you the bug is on your end. When you define the function, you need to use elemental operations like .* because curve_fit internally is evaluating the function for vector values.

Specifically, you need to write

model(x,p)=(1/(sqrt(2*pi)*p[1])) * exp(-( (x-p[2])/(sqrt(2)*p[1])) .^2)

Note the extra period before the power.

[kzapfe]

Oh, I know he was nice, anyone who responds is doing good. Was I rude? in that case my apologies, probably is because I write a bit rough. Going back to the issue, yes, it seems to be that I don't understand the point syntaxis in Julia well. How is that p[1] is NOT an element of an array? I put the point in every operator now and it works, but I still don't understand what it is supposed to mean (well I have another error now...)

[davidssmith]

No, you were not rude. Sorry if I implied that.

p[1] is an element of an array, but in the context of the evaluation of the model function, it is a scalar. The argument x of the function appears to be a scalar to you, but it is a vector inside curve_fit.

Here is a simple example.

f(x,p) = x*p[1] + p[2]
f(1,[2,3]) = 1*2 + 3 = 5
f([1,2,3],[2,3]) = [5,7,9]

You are thinking of the function like the second line, while inside curve_fit it is being called like in the third line.

[CaBeOli]

I had this problem, but I didn't understand how to use the point yet. I put in diferente places but the ERROR continuing coming. Someone can help or say where I can find a tutorial about that?
My code:

`model(t, p) = ((p[1]*p[2])/(sqrt((1/sqrt(p[3]p[4]))^2-(p[2]/(2p[3]))^2)*p[3]))*sin(sqrt((1/sqrt(p[3]p[4]))^2-(p[2]/(2p[3]))^2)*t)exp(-(p[2]/(2p[3]))*t)

fit = curve_fit(model, time, Volts, [Vm, R, L, C])`

[blakejohnson]

@CaBeOli, when your model is evaluated, the t is a Vector{Float64}. So, anywhere you have a t on both sides of a binary operator, you should use the element-wise syntax. In your case, this should be sufficient:

model(t, p) = ((p[1]*p[2])/(sqrt((1/sqrt(p[3]p[4]))^2 -
              (p[2]/(2p[3]))^2)*p[3])) * sin(sqrt((1/sqrt(p[3]p[4]))^2-(p[2]/(2p[3]))^2)*t)
                                      .* exp(-(p[2]/(2p[3]))*t)

[blakejohnson]

Maybe the more appropriate question if whether we should modify LsqFit.jl to not require the user to write their model in a vectorized form?

[IlyaOrson]

Maybe with julia v0.6 it is enough to exemplify the use of the new @. macro in the README example?

[pkofod]

Maybe with julia v0.6 it is enough to exemplify the use of the new @. macro in the README example?

I agree.