Hello … it is possible to use optim to solve this ?
R = [0.5 1 2 3] # so for each R I have one S
S = [1 2 3 4 ]
SS(R) = quadgk(a-> R*(a^2) * x[1] + R**a* x[2] ,0,Inf)
for which x[1] and x[2] i have SS = S ??? or [SS - S] → 0 ( minimization )
i have a code using neldermead … but does not work, the coefficients that i found is wrong ! matlab fminsearch find ok, but take days to evaluate …
any ideas or different methods ??
R**a is not a valid Julia expression: what is the exact integrand?
SS = quadgk(a-> R*(a^2) * x[1] + R *a * x[2] ,0,Inf)
it is a simple example … i want only to know the correct code for do that using optim.jl … neldermead.jl is not working … if i know this example, i can apply to my system … and a want to know if you know other better method to do that
What do you mean it’s not working? Is it finding a minimum but the wrong one?
The indefinite integral can be solved analytically and the primitive is
Once you plug in the limits of integration, you only get trivial results for SS, either 0 or \pm\infty according to the values of x, since R is always non-zero. If MATLAB finds any solution, it must be wrong too.
Please test your algorithm on a reasonable integral.
For Optim you can use it like that, it works very well, but keep in mind that Nelder–Mead is a local optimizer, so if you have several minima chances are you won’t find the correct one.
julia> optimize(x->sum(x^2 for x in x),ones(3),NelderMead(),Optim.Options(g_tol=1e-32)).minimizer
3-element Array{Float64,1}:
-1.07756e-16
-6.23653e-17
5.00458e-17
Please don’t double post, or if you do, then at least link to the other discussion. https://stackoverflow.com/questions/46966293/using-neldermead-jl-or-optim-which-is-the-best-for-minimize-julia
my main program is: ( only the important part ) — can you show me the imput to use optim.jl ??
aRe = [1.28; 1.47; 1.65; 1.94; 2.25; 2.39; 2.68; 2.91; 3.13; 3.75; 4.46]
aΩ11e = [2.52; 2.26; 2.07; 1.92; 1.90; 2.08; 2.31; 2.53; 2.68; 2.51; 2.27]function dΩ11(x,g,sigma,T)
R = 0.388/sigma
α = x[3]
ϕ = 0.9*(1-1/(1+(R/x[1])^x[2]))
θ = ϕ*exp(-x[3]*g^2/T)
γ = x[4]
g^5*θ*Qb(g,α)*exp(-g^2/T) + g^5*(1-θ)*γ*Qa(g)*exp(-g^2/T) +
g^5*(1-θ)*(1-γ)*Qc(g)*exp(-g^2/T)
end
Ω11(x,sigma,T) = 1/T^3*quadgk(g-> dΩ11(x,g,sigma,T) ,0 , 10, abstol=1e-3)[1]
function X2(x)
aΩ11 = zeros( lenR )
for i in lenR
aΩ11[i] = afΩ11i
end
sum((aΩ11-aΩ11e).^2)
endfor R in aRe
tic()
sigma=0.388/R
rc0,bc0,gc0 = rc0bc0gc0(sigma)
println(“sigma = $sigma, rc0 = $rc0, bc0 = $bc0, gc0 = $gc0”)
Qa = interpQa(sigma) # retorna a função Qa(g) para este sigma
Qb = interpQb(sigma) # retorna a função Qb(g,α) para este sigma
Qc = interpQc(sigma) # retorna a função Qc(g) para este sigma
push!(afΩ11, x->Ω11(x,sigma,T)) # array de funções
toc()
endxi = ones( lenR )
x, X2min = NelderMead(x->X2(x), [2.50;17.2;0.5;0.5])
#println(x[1]," “,x[2],” “,x[3])
println(x[1],” “,x[2],” “,x[3],” ",x[4])
Triple post: https://github.com/JuliaPy/PyPlot.jl/issues/327
@jonathanBieler’s reply is how
X2min = optimize(X2, [2.50,17.2,0.5,0.5], NelderMead(), Optim.Options(g_tol=1e-32)).minimizer
even though I doubt the result is what you expect since your function X2 does not depend on x.
PS: To properly format code, can you put 3 backticks before and 3 backticks after it?
I know you’re quoting Jonathan above, but allow me to recommend
min = Optim.minimizer(optimizer(f, ..., Method(), Optim.Options())
rather than accessing the fields directly, as this allows for deprecations. Now minimizer might not be that bad, as we’ve corrected the terminology now, but the minimizer field used to be called minimum and the minimum was called f_minimum. This had to be “hard deprecated” because you cannot deprecate field names. So for getting info from a results type instance, it’s advised to use the function based API if you want to avoid sudden problems in your programs should we change something.
Anyway, I don’t get what you’re trying to do, please always provide an MWE. In Optim you can use Nelder Mead using the following syntax
optimize(your_function, initial_x, NelderMead())
you can also set options
optimize(your_function, initial_x, NelderMead(), Optim.Options(kwargs...))
gives this error
MethodError: no method matching NelderMead()
Closest candidates are:
NelderMead(!Matched::Any, !Matched::Any) at C:\Users\Lucas\Desktop\LUCAS\Julia\neldermead.jl:76
NelderMead(!Matched::Any, !Matched::Any, !Matched::Any) at C:\Users\Lucas\Desktop\LUCAS\Julia\neldermead.jl:76
NelderMead(!Matched::Any, !Matched::Any, !Matched::Any, !Matched::Any) at C:\Users\Lucas\Desktop\LUCAS\Julia\neldermead.jl:76
…
include_string(::String, ::String) at loading.jl:515
include_string(::String, ::String, ::Int64) at eval.jl:30
include_string(::Module, ::String, ::String, ::Int64, ::Vararg{Int64,N} where N) at eval.jl:34
(::Atom.##49#53{String,Int64,String})() at eval.jl:50
withpath(::Atom.##49#53{String,Int64,String}, ::String) at utils.jl:30
withpath(::Function, ::String) at eval.jl:38
macro expansion at eval.jl:49 [inlined]
(::Atom.##48#52{Dict{String,Any}})() at task.jl:80
do you know what is it ?
it says that x is not defined … but x is x[1] x[2] x[3] x[4] right ?
not work … tell that x is not defined
We would gladly help you if you provided a minimal example that, except for the optimization part, we can run: the function X2 you provide is incomplete; moreover it does not depend on x so any value of x is a minimizer:
function X2(x)
aΩ11 = zeros( lenR )
for i in lenR # here you probably want for i in 1:lenR
aΩ11[i] = afΩ11i # what is afΩ11i?
end
sum((aΩ11-aΩ11e).^2)
end
As for x not being defined, if you look at the code (corrected after @pkofod’s suggestion), the minimizer is called X2min, so yes, x is not defined.
X2min = Optim.minimizer(optimize(X2, [2.50,17.2,0.5,0.5], NelderMead(), Optim.Options(g_tol=1e-32)))
Your errors point to “at C:\Users\Lucas\Desktop\LUCAS\Julia\neldermead.jl” which has nothing to do with Optim, so I’m not surprised it doesn’t work when you try to use Optim syntax.
all tips you said here not working … i will put my entire code… so, i tried this two method, but none of this works
> using NLsolve, QuadGK, Interpolations, Optim
>
> include("C:\\Users\\Lucas\\Desktop\\LUCAS\\Julia\\brent.jl")
> #include("C:\\Users\\Lucas\\Desktop\\LUCAS\\Julia\\neldermead.jl")
>
>
> # Variáveis globais
>
> aRe = [1.28; 1.47; 1.65; 1.94; 2.25; 2.39; 2.68; 2.91; 3.13; 3.75; 4.46] # Array com os R experimentais
> aΩ11e = [2.52; 2.26; 2.07; 1.92; 1.90; 2.08; 2.31; 2.53; 2.68; 2.51; 2.27] # Array com os Ω11 experimentais
> aΩ11 = zeros(length(aRe)) # Array com os valores calculados com cada função afΩ11[i](sigma)
> lenR = length(aRe)
> afΩ11 = [] # Arrray de funções: uma para cada sigma
> Qa = [] # Variável onde será atribuída a função Qa(g)
> Qb = [] # Variável onde será atribuída a função Qb(g,α)
> Qc = [] # Variável onde será atribuída a função Qc(g)
> T = 0.0987
>
> #
> # Potencial e suas derivadas
> #
>
> Φ(r,sigma) = 2/15*sigma.^9*(1/(r-1)^9-1/(r+1)^9-9/8/r*(1/(r-1)^8-1/(r+1)^8))-
> sigma.^3*(1/(r-1)^3-1/(r+1)^3-3/2/r*(1/(r-1)^2-1/(r+1)^2))
>
> function Φl(r,sigma)
> ((3*sigma.^3)/(2*r^2*(1+r)^2)+(3*sigma.^3)/(r*(1+r)^3)-(3*sigma.^3)/(1+r)^4
> -(3*sigma.^9)/(20*r^2*(1+r)^8)
> -(6*sigma.^9)/(5*r*(1+r)^9)
> +(6*sigma.^9)/(5*(1+r)^10)
> -(3*sigma.^3)/((r-1)^3*r)
> +(6*sigma.^9)/(5*(r-1)^9*r)
> -(3*sigma.^3)/(2*(r-1)^2*r^2)
> +(3*sigma.^9)/(20*(r-1)^8*r^2)
> +(3*sigma.^3)/(r-1)^4-(6*sigma.^9)/(5*(r-1)^10))
> end
>
> function Φll(r,sigma)
> ((-(3*sigma.^3)/(r^3*(1+r)^2))-(6*sigma.^3)/(r^2*(1+r)^3)
> -(9*sigma.^3)/(r*(1+r)^4)+(12*sigma.^3)/(1+r)^5
> +(3*sigma.^9)/(10*r^3*(1+r)^8)
> +(12*sigma.^9)/(5*r^2*(1+r)^9)
> +(54*sigma.^9)/(5*r*(1+r)^10)
> -(12*sigma.^9)/(1+r)^11+(9*sigma.^3)/((r-1)^4*r)
> -(54*sigma.^9)/(5*(r-1)^10*r)
> +(6*sigma.^3)/((r-1)^3*r^2)
> -(12*sigma.^9)/(5*(r-1)^9*r^2)
> +(3*sigma.^3)/((r-1)^2*r^3)
> -(3*sigma.^9)/(10*(r-1)^8*r^3)
> -(12*sigma.^3)/(r-1)^5+(12*sigma.^9)/(r-1)^11)
> end
>
> #
> # Potencial efetivo e suas derivadas
> #
>
> Veff(r,b,g,sigma) = Φ(r,sigma)+b^2*g^2/r^2
> Veffl(r,b,g,sigma) = Φl(r,sigma)-2b^2*g^2/r^3
> Veffll(r,b,g,sigma) = Φll(r,sigma)+6b^2*g^2/r^4
>
> # O rc0, bc0 e gc0 vem das equações
> # Φ(r)+b^2 g^2/r^2 = g^2
> # Φ'(r)-2b^2 g^2/r^3 = 0
> # Φ''(r)+6b^2 g^2/r^4 = 0
> #
> # Ques são equivalentes a
> # Φ(r)+b^2 g^2/r^2 = g^2
> # Φ(r)+rΦ'(r)/2 = g^2
> # 3Φ'(r)+rϕ''(r) = 0
> function rc0bc0gc0(sigma)
> rc0 = brent(r->3Φl(r,sigma)+r*Φll(r,sigma),1.02,5)
> gc0= sqrt(Φ(rc0,sigma)+rc0*Φl(rc0,sigma)/2)
> bc0 = rc0*sqrt(gc0^2-Φ(rc0,sigma))/gc0
> rc0, bc0, gc0
> end
>
> # Função que calcula o r e o b críticos para um dado g
> function frbc(g, sigma)
> rc0,bc0,gc0 = rc0bc0gc0(sigma)
> if g == gc0 return rc0, bc0 end
> if g > gc0 || g == 0 error("g = $g, there is no critical b") end
> rc = rc0
> f(r,sigma) = Φ(r,sigma)+r*Φl(r,sigma)/2-g^2
> while f(rc0,sigma)*f(rc,sigma)>0 rc = rc+1 end
> rc = brent(r->f(r,sigma),rc0,rc)
> bc = sqrt(rc^2*(g^2-Φ(rc,sigma))/g^2)
> rc, bc
> end
>
> # Cálculo do r mínimo
> function Rmin(b,g,sigma)
> if b > 100 # Este caso é quando Φ -> 0 e g^2r^2-b^2g^2+Φ=0
> rmin = b
> return rmin
> end
> rc0, bc0, gc0 = rc0bc0gc0(sigma)
> if g == gc0 return rc0 end # É o caso de b e g críticos
> if g > gc0 # Neste caso não tem extremos
> rinicial = nextfloat(1.0)
> rfinal = rinicial
> while (g^2-Veff(rinicial,b,g,sigma))*(g^2-Veff(rfinal,b,g,sigma)) > 0 rfinal = rfinal+1 end
> rmin = brent(r->g^2-Veff(r,b,g,sigma),rinicial,rfinal)
> return rmin
> end
> # Neste caso Veff tem máximo e mínimo e temos b crítico
> rc, bc = frbc(g, sigma)
> if b == bc return rc end
> if b < bc
> rinicial = nextfloat(1.0)
> rfinal = rinicial
> while (g^2-Veff(rinicial,b,g,sigma))*(g^2-Veff(rfinal,b,g,sigma)) > 0 rfinal = rfinal+1 end
> rmin = brent(r->g^2-Veff(r,b,g,sigma),rinicial,rfinal)
> else
> rfinal = rc
> while (g^2-Veff(rc,b,g,sigma))*(g^2-Veff(rfinal,b,g,sigma)) > 0 rfinal = rfinal+1 end
> rmin = brent(r->g^2-Veff(r,b,g,sigma),rc,rfinal)
> end
> rmin
> end
>
> # Esta função retorna o integrando da função que dá o X (chi)
> # A integração é feita em r (r* nos cálculos)
> # X = pi - 2b*g*integrate(1/(r^2*sqrt(g^2-Veff)),r,rmin,Inf)
> function dX(r,b,g,sigma)
> if g == 0 return 0 end
> rmin = Rmin(b,g,sigma)
> den1 = g^2 - Veff(r,b,g,sigma)
> den2 = -Veffl(rmin,b,g,sigma)*(r-rmin)
> # if den1<0 || den2 < 0
> # warn("g^2-Veff = $den1 -Veffl*(r-rmin) = $den2, r = $r, b = $b, g = $g")
> # end
> den1 = abs(den1)
> den2 = abs(den2)
> # The lines below are
> # 2b*g*/r^2(1/sqrt(den1)-1/sqrt(den2))
> t0 = 2*b*g/r^2
> t1 = 1/sqrt(den1)
> t2 = 1/sqrt(den2)
> if abs(t2/t1) < 0.5
> res = t0*t1*exp(log1p(-t2/t1))
> else
> res = t0*(t1-t2)
> end
> res
> end
>
> function X(b,g,sigma)
> if g > 1e5 return 0 end
> if g == 0 || b == 0 return π end
> rc0, bc0, gc0 = rc0bc0gc0(sigma)
> if b == bc0 && g == gc0 return -Inf end
> if 0 < g <= gc0
> rc, bc = frbc(g, sigma)
> if b == bc return -Inf end
> end
> rmin = Rmin(b,g,sigma)
> if b > 100 return 0 end
> if g < 0.2 && b > bc+0.2
> return -4π*sigma^3/g^2/b^6
> end
> int = 2b*g*pi/(2rmin^(3/2)*sqrt(abs(Veffl(rmin,b,g,sigma))))
> try
> int = int + quadgk(r->dX(r,b,g,sigma), rmin, Inf, abstol=1e-3)[1]
> catch
> println("b=$b, g=$g")
> int = int + quadgk(r->dX(r,b,g,sigma), BigFloat(rmin), Inf, abstol=1e-3)[1]
> end
> pi-int
> end
>
> function dQ(b,g,sigma)
> 2*(1-cos(X(b,g,sigma)))*b
> end
>
> function dQd(b,g,sigma,α)
> 2*(1+1/sqrt(1+α)*sqrt(π*T)/2g*sin(X(b,g,sigma)/2))*b
> end
>
> function dQi(b,g,sigma)
> 2*(1+2/3*sin(X(b,g,sigma)/2))*b
> end
>
> function Q(g,sigma)
> if g > 1e5 return 0 end
> if g == 0 return Inf end
> rc0, bc0, gc0 = rc0bc0gc0(sigma)
> if g > gc0
> int = quadgk(b->dQ(b,g,sigma), 0, Inf, abstol=1e-3)[1]
> else
> rc, bc = frbc(g, sigma)
> int = quadgk(b->dQ(b,g,sigma), 0, bc-1e-3, abstol=1e-3)[1]
> int += quadgk(b->dQ(b,g,sigma), bc+1e-3, Inf, abstol=1e-3)[1]
> end
> end
>
> function Qd(g,sigma,α)
> if g > 1e5 return 0 end
> if g == 0 return Inf end
> rc0, bc0, gc0 = rc0bc0gc0(sigma)
> if g > gc0
> int = quadgk(b->dQ(b,g,sigma), 0, Inf, abstol=1e-3)[1]
> else
> rc, bc = frbc(g, sigma)
> int = quadgk(b->dQd(b,g,sigma,α), 0, bc-1e-3, abstol=1e-3)[1]
> int = int + quadgk(b->dQ(b,g,sigma), bc+1e-3, Inf, abstol=1e-3)[1]
> end
> int
> end
>
> function Qi(g,sigma)
> if g > 1e5 return 0 end
> if g == 0 return Inf end
> rc0, bc0, gc0 = rc0bc0gc0(sigma)
> if g > gc0
> int = quadgk(b->dQ(b,g,sigma), 0, Inf, abstol=1e-3)[1]
> else
> rc, bc = frbc(g, sigma)
> int = quadgk(b->dQi(b,g,sigma), 0, bc-1e-3, abstol=1e-3)[1]
> int = int + quadgk(b->dQ(b,g,sigma), bc+1e-3, Inf, abstol=1e-3)[1]
> end
> int
> end
>
> # Função que interpola valores de Q
> # Isto é necessário porque a função Q é custosa para calcular
> # Então calculamos alguns pontos e interpolamos uma função.
> # Modo de usar:
> #Qa = interpQa(sigma)
> # Depois usamos Qa(g)
> function interpQa(sigma)
> ag = 0.1:1:10 # Array com os valores de g
> mQ = [Q(g,sigma) for g in ag] # Matriz com os valores de Q na grade
> iQ = interpolate(mQ, BSpline(Cubic(Line())), OnGrid()) # Interpola na grade
> sQ = scale(iQ, ag)
> x->sQ[x]
> end
>
> # Modo de usar:
> #Qb = interpQb(sigma)
> # Depois usamos Qb(g,α)
> function interpQb(sigma)
> ag = 0.1:1:10 # Array com os valores de g
> aα = 0:1:5
> mQ = [Qd(g,sigma,α) for g in ag, α in aα] # Matriz com os valores de Q na grade
> iQ = interpolate(mQ, BSpline(Cubic(Line())), OnGrid()) # Interpola na grade
> sQ = scale(iQ, ag, aα)
> (x,y) -> sQ[x,y]
> end
>
> # Modo de usar:
> #Qc = interpQc(sigma)
> # Depois usamos Qc(g)
> function interpQc(sigma)
> ag = 0.1:1:10 # Array com os valores de g
> mQ = [Qi(g,sigma) for g in ag] # Matriz com os valores de Q na grade
> iQ = interpolate(mQ, BSpline(Cubic(Line())), OnGrid()) # Interpola na grade
> sQ = scale(iQ, ag)
> x->sQ[x]
> end
>
> function dΩ11(x,g,sigma,T)
> R = 0.388/sigma
>
> α = 0 # x[1] # x[] são as constantes a serem calculadas pelo metodo neldermead
> ϕ = 0.9*(1-1/(1+(R/x[1])^17.32))
> θ = ϕ*exp(-α*g^2/T)
> #θ(g)=ϕ*exp(-x[1]*g^2/T)
> γ = 1 #x[3]
> g^5*θ*Qb(g,α)*exp(-g^2/T) + g^5*(1-θ)*γ*Qa(g)*exp(-g^2/T) +
> g^5*(1-θ)*(1-γ)*Qc(g)*exp(-g^2/T)
> end
>
> Ω11(x,sigma,T) = 1/T^3*quadgk(g-> dΩ11(x,g,sigma,T) ,0 , 10, abstol=1e-3)[1]
>
> function X2(x)
> aΩ11 = zeros( lenR )
> for i in lenR
> aΩ11[i] = afΩ11[i](x)
> end
> sum((aΩ11-aΩ11e).^2)
> end
>
>
> for R in aRe
> tic()
> sigma=0.388/R
> rc0,bc0,gc0 = rc0bc0gc0(sigma)
> println("sigma = $sigma, rc0 = $rc0, bc0 = $bc0, gc0 = $gc0")
> Qa = interpQa(sigma) # retorna a função Qa(g) para este sigma
> Qb = interpQb(sigma) # retorna a função Qb(g,α) para este sigma
> Qc = interpQc(sigma) # retorna a função Qc(g) para este sigma
> push!(afΩ11, x->Ω11(x,sigma,T)) # array de funções
> toc()
> end
>
> xi = ones( lenR )
> X2min = optimize(X2, 2.50, Optim.Options(g_tol=1e-32)).minimizer
>
> X2min = Optim.minimizer(optimize(X2, 2.50, NelderMead(), Optim.Options(g_tol=1e-32)))
>
>
>
>
> #x, X2min = NelderMead(x->X2(x), [2.50;17.2;0.5;0.5])
> #println(x[1]," ",x[2]," ",x[3])
> println(x[1])
Please do make sure that the code you paste is:
include at the beginning) and some variables are not definedx when the optimal solution is called X2min)The following code works:
using Optim
aΩ11e = [2.52; 2.26; 2.07; 1.92; 1.90; 2.08; 2.31; 2.53; 2.68; 2.51; 2.27]
lenR = length(aΩ11e)
function X2(x)
aΩ11 = zeros(lenR)
#for i in 1:lenR
# aΩ11[i] = afΩ11i # OP does not post values of LHS variable
#end
sum((aΩ11-aΩ11e).^2)
end
X2min = Optim.minimizer(optimize(X2, [2.50, 1.3], NelderMead(), Optim.Options(g_tol=1e-32)))
println(X2min)
o choose only one parameters to fit , x[1] … do you know how to solve this error ?
X2min = Optim.minimizer(optimize(X2, [2.50], NelderMead(), Optim.Options(g_tol=1e-32)))
println(X2min)
and the error is : Use optimize(f, scalar, scalar) for 1D problems
in optimize(::Function, ::Array{Float64,1}, ::Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Void}) at optimize.jl:168
in include_string(::String, ::String) at loading.jl:441
in include_string(::String, ::String, ::Int64) at eval.jl:30
in include_string(::Module, ::String, ::String, ::Int64, ::Vararg{Int64,N}) at eval.jl:34
in (::Atom.##53#56{String,Int64,String})() at eval.jl:50
in withpath(::Atom.##53#56{String,Int64,String}, ::String) at utils.jl:30
in withpath(::Function, ::String) at eval.jl:38
in macro expansion at eval.jl:49 [inlined]
in (::Atom.##52#55{Dict{String,Any}})() at task.jl:60