Optim -> LeastSquaresOptim + ADerrors update

Manifest.toml

@@ -2,17 +2,23 @@
 
 [[ADerrors]]
 deps = ["BDIO", "FFTW", "FastGaussQuadrature", "ForwardDiff", "LaTeXStrings", "LinearAlgebra", "Nettle", "PGFPlotsX", "Plots", "Printf", "QuadGK", "Roots", "Statistics", "Test", "UnicodePlots"]
-git-tree-sha1 = "f1c9954c17c0e87f185d532163e9e630efd12fe3"
-repo-rev = "9c41ec2e161e4826bfe0752be8bf2a0bc025759c"
+git-tree-sha1 = "2d3b701de1b206f497c632413b5a783fbd65d714"
+repo-rev = "master"
 repo-url = "https://gitlab.ift.uam-csic.es/alberto/aderrors.jl"
 uuid = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 version = "0.1.0"
 
 [[AbstractFFTs]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "8ed9de2f1b1a9b1dee48582ad477c6e67b83eb2c"
+git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "1.0.0"
+version = "1.0.1"
+
+[[Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.2.0"
 
 [[ArgCheck]]
 git-tree-sha1 = "dedbbb2ddb876f899585c4ec4433265e3017215a"
@@ -197,9 +203,9 @@ version = "4.3.1+4"
 
 [[FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
-git-tree-sha1 = "c31e446bf3b12aad2ec8fc500fe19528c148d811"
+git-tree-sha1 = "1b48dbde42f307e48685fa9213d8b9f8c0d87594"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.3.1"
+version = "1.3.2"
 
 [[FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -275,9 +281,9 @@ version = "6.1.2+6"
 
 [[GR]]
 deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
-git-tree-sha1 = "b90b826782cb3ac5b7a7f41b3fd0113180257ed4"
+git-tree-sha1 = "12d971c928b7ecf19b748a2c7df6a365690dbf2c"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.53.0"
+version = "0.55.0"
 
 [[GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt_jll", "Zlib_jll", "libpng_jll"]
@@ -287,9 +293,9 @@ version = "0.53.0+0"
 
 [[GeometryBasics]]
 deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "a28d728c2d825285fe27f38ca322399d35d1a5b9"
+git-tree-sha1 = "4d4f72691933d5b6ee1ff20e27a102c3ae99d123"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-version = "0.3.6"
+version = "0.3.9"
 
 [[Gettext_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
@@ -309,10 +315,10 @@ uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
 version = "1.0.0"
 
 [[HTTP]]
-deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
-git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508"
+deps = ["Base64", "Dates", "IniFile", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
+git-tree-sha1 = "c9f380c76d8aaa1fa7ea9cf97bddbc0d5b15adc2"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "0.8.19"
+version = "0.9.5"
 
 [[IOCapture]]
 deps = ["Logging"]
@@ -386,6 +392,18 @@ git-tree-sha1 = "3a0084cec7bf157edcb45a67fac0647f88fe5eaf"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 version = "0.14.7"
 
+[[LazyArtifacts]]
+deps = ["Pkg"]
+git-tree-sha1 = "4bb5499a1fc437342ea9ab7e319ede5a457c0968"
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+version = "1.3.0"
+
+[[LeastSquaresOptim]]
+deps = ["FiniteDiff", "ForwardDiff", "LinearAlgebra", "Optim", "Printf", "SparseArrays", "Statistics", "SuiteSparse"]
+git-tree-sha1 = "b5a1931bf37616820c4ef71629b6308d67f0b393"
+uuid = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
+version = "0.8.1"
+
 [[LibGit2]]
 deps = ["Printf"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
@@ -467,10 +485,10 @@ uuid = "2fda8390-95c7-5789-9bda-21331edee243"
 version = "0.12.0"
 
 [[MKL_jll]]
-deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
-git-tree-sha1 = "eb540ede3aabb8284cb482aa41d00d6ca850b1f8"
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2020.2.254+0"
+version = "2021.1.1+1"
 
 [[MacroTools]]
 deps = ["Markdown", "Random"]
@@ -531,6 +549,11 @@ git-tree-sha1 = "b3b169cec36b850cfacbefc0b0a5f68a6089eb70"
 uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
 version = "3.4.1+2"
 
+[[NetworkOptions]]
+git-tree-sha1 = "ed3157f48a05543cce9b241e1f2815f7e843d96e"
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
+
 [[Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "a42c0f138b9ebe8b58eba2271c5053773bde52d0"
@@ -614,9 +637,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[PlotThemes]]
 deps = ["PlotUtils", "Requires", "Statistics"]
-git-tree-sha1 = "c6f5ea535551b3b16835134697f0c65d06c94b91"
+git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
 uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "2.0.0"
+version = "2.0.1"
 
 [[PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
@@ -626,9 +649,9 @@ version = "1.0.10"
 
 [[Plots]]
 deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
-git-tree-sha1 = "4797acb266b8d9ff316f4581924e71c6709f152d"
+git-tree-sha1 = "142dd04f5060c04de91cc10ca76dffb291a02426"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "1.10.1"
+version = "1.10.6"
 
 [[PositiveFactorizations]]
 deps = ["LinearAlgebra"]
@@ -679,9 +702,9 @@ version = "1.1.1"
 
 [[RecipesPipeline]]
 deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
-git-tree-sha1 = "9ea2f5bf1b26918b16e9f885bb8e05206bfc2144"
+git-tree-sha1 = "c4d54a78e287de7ec73bbc928ce5eb3c60f80b24"
 uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
-version = "0.2.1"
+version = "0.3.1"
 
 [[Reexport]]
 deps = ["Pkg"]
@@ -709,9 +732,9 @@ version = "0.2.2+1"
 
 [[Roots]]
 deps = ["Printf"]
-git-tree-sha1 = "8f743e4f4368d1d753f3806bf635899dad6b4847"
+git-tree-sha1 = "369e25546984dff5df351bc056fccc30de615080"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "1.0.7"
+version = "1.0.8"
 
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -777,10 +800,10 @@ uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 version = "0.9.6"
 
 [[StructArrays]]
-deps = ["DataAPI", "Tables"]
-git-tree-sha1 = "ad1f5fd155426dcc879ec6ede9f74eb3a2d582df"
+deps = ["Adapt", "DataAPI", "Tables"]
+git-tree-sha1 = "26ea43b4be7e919a2390c3c0f824e7eb4fc19a0a"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.4.2"
+version = "0.5.0"
 
 [[SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -794,14 +817,19 @@ version = "1.0.0"
 
 [[Tables]]
 deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
-git-tree-sha1 = "240d19b8762006ff04b967bdd833269ad642d550"
+git-tree-sha1 = "a716dde43d57fa537a19058d044b495301ba6565"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.2.2"
+version = "1.3.2"
 
 [[Test]]
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[[URIs]]
+git-tree-sha1 = "7855809b88d7b16e9b029afd17880930626f54a2"
+uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
+version = "1.2.0"
+
 [[UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -983,9 +1011,9 @@ version = "1.2.11+18"
 
 [[Zstd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "6f1abcb0c44f184690912aa4b0ba861dd64f11b9"
+git-tree-sha1 = "2c1332c54931e83f8f94d310fa447fd743e8d600"
 uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
-version = "1.4.5+2"
+version = "1.4.8+0"
 
 [[libass_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 ADerrors = "5e92007d-7bf1-471c-8ceb-4591b8b567a9"
 BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, LeastSquaresOptim
 import Statistics: mean
 
 include("juobs_types.jl")

src/juobs_tools.jl

@@ -29,21 +29,6 @@ function apply_rw(data::Vector{<:Array{Float64}}, W::Vector{Matrix{Float64}})
     return data_r
 end
 
-function dobsdp(a::uwreal, p::uwreal) # Compute da / dp
-    
-    if count(p.prop .== true) != 1
-        error("I do not know how to compute this")
-    end
-    
-    for i = 1:min(length(a.der), length(p.der))
-        if (p.prop[i] && a.prop[i])
-            return a.der[i] / p.der[i]
-        end 
-    end 
-    
-    return 0.0
-end
-
 function check_corr_der(obs::Corr, derm::Vector{Corr})   
     g1 = Vector{String}(undef, 0)
     g2 = Vector{String}(undef, 0)
@@ -141,7 +126,7 @@ function corr_obs(cdata::Array{CData, 1}; real::Bool=true, rw::Union{Array{Array
     return Corr(obs, cdata)
 end    
 
-#TODO: VECTORIZE
+#TODO: VECTORIZE, uwreal?
 @doc raw"""
     md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADerrors.wsg)
 
@@ -195,8 +180,8 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
     if !all(id .== id[1])
         error("ids do not match")
     end
-    id = id[1]
-
+    id = ws.id2str[id[1]]
+    
     ivrep = getfield.(ws.fluc[p], :ivrep)
     ivrep1 = fill(ivrep[1], length(ivrep))
     if !all(ivrep .== ivrep1)
@@ -208,6 +193,7 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
         error("Nr obs != Nr md")
     end
 
+    #md_aux as a Matrix + Automatic truncation
     md_aux = md[1][:, 1:ivrep[1]]
     for k = 2:length(md)
         md_aux = cat(md_aux, md[k][:, 1:ivrep[k]], dims=2)
@@ -215,27 +201,21 @@ function md_sea(a::uwreal, md::Vector{Matrix{Float64}}, ws::ADerrors.wspace=ADer
 
     fluc_obs = getfield.(ws.fluc[p], :delta)
     fluc_md = md_aux .- mean(md_aux, dims=2)
-    
+    uwerr(a)
+    fluc_obs = mchist(a, id)
+
     nrw = size(fluc_md, 1)
-    d = a.der[p]
-    nobs = sum(a.prop)
     if nrw == 1
-        der1 = uwreal(0)
-        for k = 1:nobs
-            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
-        end
+        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
         return (der1, der1)
     elseif nrw == 2
-        der1 = uwreal(0)
-        der2 = uwreal(0)
-        for k = 1:nobs
-            der1 = der1 - d[k] * uwreal(fluc_md[1, :] .* fluc_obs[k], id, ivrep)
-            der2 = der2 - d[k] * uwreal(fluc_md[2, :] .* fluc_obs[k], id, ivrep)
-        end
+        der1 = uwreal(-fluc_md[1, :] .* fluc_obs, id, ivrep)
+        der2 = uwreal(-fluc_md[2, :] .* fluc_obs, id, ivrep)
         return (der1, der2)
     else
         return nothing
     end
+    
 end
 
 @doc raw"""
@@ -286,7 +266,7 @@ function md_val(a::uwreal, obs::Corr, derm::Vector{Corr})
 
     corr = getfield(obs, :obs)
 
-    der = [dobsdp(a, corr[k]) for k = 1:length(corr)]
+    der = [derivative(a, corr[k]) for k = 1:length(corr)]
     derm1, derm2 = derm
     return (sum(der .* derm1.obs), sum(der .* derm2.obs))
 end
@@ -442,14 +422,18 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         C = isnothing(wpm) ? [ADerrors.cov(aux[k]) for k = 1:Ndata] : [ADerrors.cov(aux[k], wpm) for k = 1:Ndata]
         chisq_full_cov(p, d) = get_chi2_cov(model, d, C, p, Nalpha)
         min_fun_cov(t) = chisq_full_cov(t, dat)
-        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
+        sol = optimize(min_fun_cov, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
         
-        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, Optim.minimizer(sol), data) : fit_error(chisq_full_cov, Optim.minimizer(sol), data, wpm)
+        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full_cov, sol.minimizer, data) : fit_error(chisq_full_cov, sol.minimizer, data, wpm)
+        println("Chisq / chiexp: ", min_fun_cov(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
     else
         chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
         min_fun(t) = chisq_full(t, dat)
-        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), method=LBFGS())
-        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, Optim.minimizer(sol), data) : fit_error(chisq_full, Optim.minimizer(sol), data, wpm)
+        sol = optimize(min_fun, vcat(fit.param, dat[1:Nalpha*Ndata]), LevenbergMarquardt())
+        
+        (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, sol.minimizer, data) : fit_error(chisq_full, sol.minimizer, data, wpm)
+        println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
+
     end
 
     #### chisq_full, min_fun out of conditional ->
@@ -461,7 +445,6 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         print("\n Fit parameter: ", i, ": ")
         details(upar[i])
     end
-    println("Chisq / chiexp: ", sol.minimum, " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
     return upar
     
 end