Merge pull request #48 from MartaVanin/master

v0.6.0: Polynomials and GMM

Author: MartaVanin

CHANGELOG.md

@@ -1,5 +1,12 @@
 ## Staged
 
+## v0.6.0
+- drop GaussianMixtures.jl dependency in favour of using MixtureModels from Distributions.jl
+- add support for polynomials in mle (adds Polynomials dependency!)
+- add `detect_ambiguities` at the end of the tests
+- update todo.md
+- extends _PMDSE.rand to host AbstractRNG and, e.g., seed! the random sample generation
+
 ## v0.5.0
 - add support JuMP 0.22, PMD 0.12,0.13, IM 0.7
 - nicer ternary operators

Project.toml

@@ -1,21 +1,21 @@
 name = "PowerModelsDistributionStateEstimation"
 uuid = "d0713e65-ce0c-4a8e-a1da-2ed737d217c5"
 authors = ["Marta Vanin", "Tom Van Acker"]
-version = "0.5.0"
+version = "0.6.0"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GaussianMixtures = "cc18c42c-b769-54ff-9e2a-b28141a64aae"
 InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 PowerModelsDistribution = "d7431456-977f-11e9-2de3-97ff7677985e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -26,12 +26,12 @@ CSV = "~0.6.1"
 DataFrames = "~0.21.0"
 Distributions = "0.22.3 - 0.24.14"
 ForwardDiff = "~0.10.18"
-GaussianMixtures = "0.3.2"
 InfrastructureModels = "~0.6, ~0.7"
 JSON = "~0.18, ~0.19, ~0.20, ~0.21"
 JuMP = "~0.21, ~0.22"
 LoggingExtras = "~0.4"
 Optim = "~1.2.0"
+Polynomials = "~2.0.17"
 PowerModelsDistribution = "~0.11, ~0.12, ~0.13"
 SpecialFunctions = "~0.10.3"
 julia = "^1"

TODO.md

@@ -1,30 +1,43 @@
-## A list of things to do for the next release (0.4.1)
+## A list of things to do for the next release (0.7.0)
 
-- [ ] update basic notebook (add bad data in it), add non-Gaussian notebook and more non-Gaussian docs
+- [ ] MAJOR/BREAKING: upgrade CSV
+- [ ] MAJOR/BREAKING: after CSV, upgrade PMD, Ipopt, JuMP
+      - [ ] move LTS to julia 1.6 and drop support for JuMP < 0.22, and consequently PMD/Ipopt versions
+
+- [ ] import/export pdf, etc. from Distributions.jl? simplify residual constraint
+
+- [ ] update basic notebook (add bad data), add non-Gaussian notebook and more non-Gaussian docs
 
 - [ ] input through array of measurements/dataframe rather than csv
       - allows to create these directly from powerflow results without creating csv files
       - replace mktempdir etc. from the tests and use this once it is up and running
+      - cannot remove CSV dep because of ENWL parsing, but maybe there is a workaround? (like putting ENWL data in another repo except tests, see below)
 
 - [ ] consider deprecating reduced_ac and reduced_ivr after test against @smart_constraint (especially reduced_ac)
 
 - [ ] increase coverage, in particular:
       - test "ls" and "lav" (no weights)
       - test for rand(ExtendedBeta)
-      - test for GMM grad/heslogpdf
+      - test for mles with various distributions
+      - test with more measurement conversions
       - fix test of line 49-50 in pseudo_measurements.jl: in 50, NUMERICAL_ERROR when sol is correct, in 49 EXCEPTION_ACCESS_VIOLATION at 0x2e2075e6 -- mumps_cst_amf_ in windows CI
-      - re-introduce start_values in tests and enable ubuntu tests back in CI (and maybe julia 1.0)
+      - re-introduce start_values in tests and enable ubuntu tests back in CI (julia 1.6 as LTS will move, see below)
       - lines 228, 232-234 in bad_data.jl and move atol to 1.1e-3 in with_errors.jl line 186
 
+- [ ] increase dependency bound on Distributions.jl
+- [ ] nw_id_default: _IM to _PMD. to rm _IM dep, sol_component_value should should be replaced (ask _PMD people to export it?)
+
 ## TODO for future releases/research wishlist
 
-- [ ] try rescale ENWL database per MVA                    (?)
-- [ ] test example notebooks                               (?)
+- [ ] facilitate change ENWL database power_base           (?)
+- [ ] add tests on example notebooks                       (?)
 - [ ] intuitive/automatic inclusion of load/transfo models (?)
       - or MV/LV notebook?
 - [ ] additional (convex?) formulations                    (?)
 - [ ] advanced bad data functionalities                    (?)
 - [ ] consider other robust estimators, e.g. schweppe huber(?)
       - with MOI/JuMP complementarity constraints ?
       - or ConditionalJuMP.jl ? Or other?
-- [ ] find/build 4-wire test case and test                  
+- [ ] find/build 4-wire test case and test          
+- [ ] de-localize ENWL dataset to other repo (except feeders used in tests)       
+- [ ] impl. own Gauss-Newton solver and matrix functions   (?) 

src/PowerModelsDistributionStateEstimation.jl

@@ -9,18 +9,17 @@ module PowerModelsDistributionStateEstimation
 
 # import pkgs
 import Base: maximum, minimum, rand
-import CSV
-import DataFrames
+import CSV, DataFrames
 import Distributions
 import Distributions: logpdf, gradlogpdf
 import ForwardDiff
-import GaussianMixtures
 import InfrastructureModels
 import JuMP
 import LinearAlgebra
 import LinearAlgebra: diag, diagm, I
 import Logging, LoggingExtras
 import Optim
+import Polynomials as _Poly
 import PowerModelsDistribution
 import PowerModelsDistribution: _has_nl_expression #need this to use @smart_constraint
 import Random
@@ -35,7 +34,6 @@ export optimizer_with_attributes
 const _CSV = CSV
 const _DFS = DataFrames
 const _DST = Distributions
-const _GMM = GaussianMixtures
 const _IM  = InfrastructureModels
 const _PMD = PowerModelsDistribution
 const _PMDSE = PowerModelsDistributionStateEstimation
@@ -79,6 +77,7 @@ include("form/reduced_ivr.jl")
 include("io/distributions.jl")
 include("io/measurement_parser.jl")
 include("io/network_parser.jl")
+include("io/polynomials.jl")
 include("io/postprocessing.jl")
 
 include("prob/se.jl")

src/core/constraint.jl

@@ -23,7 +23,7 @@ function constraint_mc_residual(pm::_PMD.AbstractUnbalancedPowerModel, i::Int; n
 
     for (idx, c) in enumerate(conns)
         if (occursin("ls", crit) || occursin("lav", crit)) && isa(dst[idx], _DST.Normal)
-            μ, σ = occursin("w", crit) ? (_DST.mean(dst[idx]), _DST.std(dst[idx])) :  (_DST.mean(dst[idx]), 1.0)
+            μ, σ = occursin("w", crit) ? (_DST.mean(dst[idx]), _DST.std(dst[idx])) : (_DST.mean(dst[idx]), 1.0)
         end
         if isa(dst[idx], Float64)
             JuMP.@constraint(pm.model, var[c] == dst[idx])
@@ -48,11 +48,12 @@ function constraint_mc_residual(pm::_PMD.AbstractUnbalancedPowerModel, i::Int; n
                 res[idx] * rsc * σ >= - (var[c] - μ)
             )
         elseif crit == "mle"
-            pkg_id = ( isa(dst[idx], ExtendedBeta{Float64}) || isa(dst[idx], _GMM.GMM) ) ? _PMDSE : _DST
-            lb = ( !isa(dst[idx], _GMM.GMM) && !isinf(pkg_id.minimum(dst[idx])) ) ? pkg_id.minimum(dst[idx]) : -10
-            ub = ( !isa(dst[idx], _GMM.GMM) && !isinf(pkg_id.maximum(dst[idx])) ) ? pkg_id.maximum(dst[idx]) : 10
-            if isa(dst[idx], _GMM.GMM) lb = _PMD.ref(pm, nw, :meas, i, "min") end
-            if isa(dst[idx], _GMM.GMM) ub = _PMD.ref(pm, nw, :meas, i, "max") end
+            #TODO: enforce min and max in the meas dictionary and just with haskey make it optional for extendedbeta
+            pkg_id = any([ dst[idx] isa d for d in [ExtendedBeta{Float64}, _Poly.Polynomial]]) ? _PMDSE : _DST
+            lb = ( !isa(dst[idx], _DST.MixtureModel) && !isinf(pkg_id.minimum(dst[idx])) ) ? pkg_id.minimum(dst[idx]) : -10
+            ub = ( !isa(dst[idx], _DST.MixtureModel) && !isinf(pkg_id.maximum(dst[idx])) ) ? pkg_id.maximum(dst[idx]) : 10
+            if any([ dst[idx] isa d for d in [_DST.MixtureModel, _Poly.Polynomial]]) lb = _PMD.ref(pm, nw, :meas, i, "min") end
+            if any([ dst[idx] isa d for d in [_DST.MixtureModel, _Poly.Polynomial]]) ub = _PMD.ref(pm, nw, :meas, i, "max") end
 
             shf = abs(Optim.optimize(x -> -pkg_id.logpdf(dst[idx],x),lb,ub).minimum)
             f = Symbol("df_",i,"_",c)

src/io/distributions.jl

@@ -136,8 +136,10 @@ end
 rand(dst::ExtendedBeta, N::Int) =
     (dst.max - dst.min) .* rand(_DST.Beta(dst.α, dst.β), N) .+ dst.min
 
-## Gamma
-# functions
+rand(r::_RAN.AbstractRNG, dst::ExtendedBeta, N::Int) =
+    (dst.max - dst.min) .* rand(r, _DST.Beta(dst.α, dst.β), N) .+ dst.min
+
+## Gamma heslogpdf
 function heslogpdf(d::_DST.Gamma{T}, x::Real) where T<:Real
     if _DST.insupport(_DST.Gamma, x)
         α, θ = _DST.params(d)
@@ -146,30 +148,24 @@ function heslogpdf(d::_DST.Gamma{T}, x::Real) where T<:Real
         zero(T)
     end
 end
-
-## logpdf, gradlogpdf and heslogpdf for Gaussian mixture
-#  The GMM type is from package GaussianMixtures
-function logpdf(d::_GMM.GMM{Float64,Array{Float64,2}}, x::Real)
-    σ = sqrt.(d.Σ)
-    μ = d.μ
-    g = _DST.Normal.(μ, σ)
-    log(sum([d.w[i]*_DST.pdf(g[i],x) for i in 1:d.n]))
-end
-
-function gradlogpdf(d::_GMM.GMM{Float64,Array{Float64,2}}, x::Real)
-    σ = sqrt.(d.Σ)
-    μ = d.μ
-    γ = d.w./(sqrt(2*π)*σ)
-    sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:d.n])/sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:d.n])
+##
+# gradlogpdf and heslogpdf for Gaussian mixtures done directly 
+# with Distributions.jl
+function gradlogpdf(d::_DST.MixtureModel, x::Real)
+    @assert all([c isa _DST.Normal for c in d.components]) "Only Gaussian mixture models are supported!"
+    σ = [c.σ for c in d.components]
+    μ = [c.μ for c in d.components]
+    γ = d.prior.p./(sqrt(2*π)*σ)
+    sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])/sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)])
 end
 
-function heslogpdf(d::_GMM.GMM{Float64,Array{Float64,2}}, x::Real)
-    σ = sqrt.(d.Σ)
-    μ = d.μ
-    γ = d.w./(sqrt(2*π)*σ)
-    p1 = sum([γ[i]*(x-μ[i])^2*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-4) - γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:d.n])    
-    p2 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:d.n]) 
-    p3 = sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:d.n])^2 
-    p4 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:d.n])^2
+function heslogpdf(d::_DST.MixtureModel, x::Real)
+    σ = [c.σ for c in d.components]
+    μ = [c.μ for c in d.components]
+    γ = d.prior.p./(sqrt(2*π)*σ)
+    p1 = sum([γ[i]*(x-μ[i])^2*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-4) - γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])    
+    p2 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)]) 
+    p3 = sum([-γ[i]*(x-μ[i])*exp(-0.5*(x-μ[i])^2/σ[i]^2)*σ[i]^(-2) for i in 1:length(d.components)])^2 
+    p4 = sum([γ[i]*exp(-0.5*(x-μ[i])^2/σ[i]^2) for i in 1:length(d.components)])^2
     return p1/p2-p3/p4
 end

src/io/polynomials.jl

@@ -0,0 +1,6 @@
+# maximum(p::_Poly.Polynomial) = Base.maximum(p)
+# minimum(p::_Poly.Polynomial) = Base.minimum(p)
+pdf(p::_Poly.Polynomial) = error("`pdf` of a polynomial is not definied because it is expected that the polynomial used is an approximation of the `logpdf` instead!")
+logpdf(p::_Poly.Polynomial, x::Real) = p(x)
+gradlogpdf(p::_Poly.Polynomial, x::Real) = _Poly.derivative(p, 1)(x)
+heslogpdf(p::_Poly.Polynomial, x::Real) = _Poly.derivative(p, 2)(x)

test/distributions.jl

@@ -26,14 +26,22 @@ dst = ExtendedBeta(2.0,4.0,10.0,100.0)
 @test isapprox(_PMDSE.heslogpdf(dst,50.0),                -0.001825,      atol=1e-8)
 
 # Testing logpdf/gradlogpdf/heslogpdf of gaussian mixture models
-beta = _PMDSE.ExtendedBeta(1.6339, 20.9022, -1.0e-6, 8.268e-5)
-gmm = _GMM.GMM(1, rand(beta, 100000))
-g1 = _DST.Normal(gmm.μ[1], sqrt(gmm.Σ[1]))
+beta = _PMDSE.ExtendedBeta(1.6339, 20.9022, 0.001, 2.0)
+gmm = _DST.MixtureModel([_DST.Normal(0.14585265, sqrt(0.01144241))], [1.0])
+g1 = _DST.Normal(gmm.components[1].μ, gmm.components[1].σ)
 @test isapprox(_DST.logpdf(g1, beta.min),      _PMDSE.logpdf(gmm, beta.min),     atol = 1e-8)
 @test isapprox(_DST.logpdf(g1, beta.max),      _PMDSE.logpdf(gmm, beta.max),     atol = 1e-8)
 @test isapprox(_DST.gradlogpdf(g1, beta.min),  _PMDSE.gradlogpdf(gmm, beta.min), atol = 1e-8)    
 @test isapprox(_DST.gradlogpdf(g1, beta.max),  _PMDSE.gradlogpdf(gmm, beta.max), atol = 1e-8)
 @test isapprox(_PMDSE.heslogpdf(g1, beta.min), _PMDSE.heslogpdf(gmm, beta.min),  atol = 1e-2)
 @test isapprox(_PMDSE.heslogpdf(g1, beta.max), _PMDSE.heslogpdf(gmm, beta.max),  atol = 1e-2)
 
+p = _Poly.Polynomial([0,1,1])
+@test isapprox(_PMDSE.logpdf(p, 0.0), 0.0, atol=1e-8)
+@test isapprox(_PMDSE.logpdf(p, 1.0), 2.0, atol=1e-8)
+@test isapprox(_PMDSE.gradlogpdf(p, 1.0), 3.0, atol=1e-8)
+@test isapprox(_PMDSE.gradlogpdf(p, 0.0), 1.0, atol=1e-8)
+@test isapprox(_PMDSE.heslogpdf(p, 0.0), 2.0, atol=1e-8)
+@test isapprox(_PMDSE.heslogpdf(p, 1.0), 2.0, atol=1e-8)
+
 end

test/runtests.jl

@@ -7,19 +7,20 @@
 ################################################################################
 
 # using pkgs
-using Distributions, GaussianMixtures
+using Distributions
 using HDF5
 using Ipopt
+using Polynomials
 using PowerModels, PowerModelsDistribution
 using PowerModelsDistributionStateEstimation
 #using SCS #removed while SDP tests are not active
 using Test
 
 # pkg const
 const _DST = Distributions
-const _GMM = GaussianMixtures
 const _PMD = PowerModelsDistribution
 const _PMDSE = PowerModelsDistributionStateEstimation
+const _Poly = Polynomials
 
 #network and feeder from ENWL for tests
 ntw, fdr = 4, 2
@@ -54,3 +55,7 @@ ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer,"max_cpu_time"=>300.0,
     include("with_errors.jl")
 
 end
+ambiguities = Test.detect_ambiguities(PowerModelsDistributionStateEstimation);
+if !isempty(ambiguities)
+    println("ambiguities detected: $ambiguities")
+end