fixed sampling base

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "fe3590824bb3d254aeac25ec225a0e78cf93c68c"
+project_hash = "8855e078cc96f94368ad9ff0a8d624d7133c88f7"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245"

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 EmbeddedGraphs = "9c1af47c-29f5-11e9-0b47-9f5334384f20"
 HaltonSequences = "13907d55-377f-55d6-a9d6-25ac19e11b95"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"

example/pypsaeur/main.jl

@@ -69,7 +69,7 @@ println("\tWind\t", avg_wind / avg_load * 100,  "%")
 println("\tSolar\t", avg_solar / avg_load * 100, "%")
 println("\tRenewables\t", avg_renewables / avg_load * 100, "%")
 println("\tAll sources\t", avg_prod / avg_load * 100, "%")
-# [Q]: why is this number so low? 
+# [Q]: why is this number so low? how much is going into storage?
 
 
 ## Calculate power injection and transmission distribution factor
@@ -91,3 +91,11 @@ K = incidence_matrix(g, oriented=true)
 B = diagm( ones(nv(g)) ) #???
 H = Ω * K' * pinv(B)
 # [Q]: what is the nodal succeptance matrix B?
+
+# power flow
+f = Matrix( hcat([ H * [p[i, 2:end]...]  for i=1:size(p,1)]... )')
+
+## Time averages injection and flow
+avg_p = mean.(eachcol(p[!, 2:end]))
+avg_f = mean(f, dims=1)[1,:]
+

src/MultilevelChainSampler.jl

-# This file is part of the MultilevelChainSampler.jl package, which is licensed under the GPLv3 license.
 
 module MultilevelChainSampler
 
+include("base.jl")
+export value, propose, evaluate
+export Sample
+export StaticProposal, RandomWalk, CyclicWalk
+export LogDensity
 
 end
\ No newline at end of file

src/___debug__.jl

+# Developent file, run this to load the components of the package into the current environment 
+
+## Copy from src/MultilevelChainSampler.jl
+
+# commented out to avoid redefinition error
+#----------------------------------------
+#module MultilevelChainSampler 
+
+include("base.jl")
+export State
+export LogDensity, logdensity
+
+include("propose.jl")
+export initialize, propose
+export StaticProposal
+export RandomWalkProposal 
+export CyclicWalkProposal
+
+#end (module) 
+#----------------------------------------
+

src/base.jl

 
-# State  
-abstract type AbstractState end
-struct State{X} <: AbstractState
-    state::X
+
+using Random
+using Distributions
+using AbstractMCMC
+
+## === Sample 
+abstract type AbstractSample end
+
+# single level sample
+abstract type SimpleSample <: AbstractSample end
+
+""" simplest sample, just wraps a value
+   (this is used for conventional single-level MCMC)"""
+struct Sample{X} <: SimpleSample
+    sample::X
 end
-Base.eltype(s::State{X}) where {X} = X 
-Base.convert(T::Type, s::State) = convert(T, s.state)
+Base.eltype(s::Sample{X}) where {X} = X
+value(s::Sample) = s.sample
+
+
+## === Proposal 
+abstract type AbstractProposal{S<:AbstractSample,issymmetric} end
+SymmetricProposal = AbstractProposal{S,true} where {S<:AbstractSample}
+AsymmetricProposal = AbstractProposal{S,false} where {S<:AbstractSample}
+    
+"""
+    To define a proposal implement the sampling methods 
+        - initialization x₀ ~ p₀(⋅) by `x₀ = propose(p)`  
+        - and transition y ~ q(⋅|x) by `y  = propose(p, x)`
 
+    For asymmetric proposals also implement method
+        for calculating  q(y|x)  by `logpdf(p,x,y)` 
+"""
+function propose end
+propose(p::AbstractProposal; kwargs...) = propose(Random.default_rng(), p; kwargs...)
+propose(p::AbstractProposal, x::AbstractSample; kwargs...) = propose(Random.default_rng(), p, x; kwargs...)
 
-# log probability density function 
-abstract type AbstractLogDensity end
-struct LogDensity{F <: Function, S <: AbstractState} <: AbstractLogDensity
-    logdensity::F
+logpdfratio(p::SymmetricProposal, x::AbstractSample, y::AbstractSample) = 0
+logpdfratio(p::AsymmetricProposal, x::AbstractSample, y::AbstractSample) = logpdf(p,x,y) - logpdf(p,y,x)
+
+## Single-level proposals
+SimpleProposal = AbstractProposal{S,issymmetric} where {S<:SimpleSample,issymmetric}
+SimpleSymmetricProposal = SymmetricProposal{S} where {S<:SimpleSample}
+
+
+# Static, sample independently of previous state 
+struct StaticProposal{D<:Distribution,S} <: SimpleSymmetricProposal{S}
+    distribution::D
+end
+function StaticProposal(d::Distribution)
+    S = Sample{typeof(rand(d))}
+    return StaticProposal{typeof(d),S}(d)
 end
+propose(rng::AbstractRNG, p::StaticProposal) = Sample(rand(rng, p.distribution))
+propose(rng::AbstractRNG, p::StaticProposal, x) = propose(rng, p)
+#logpdf(p::StaticProposal, x::Sample) = logpdf(p.distribution, x.sample) # won't be used
+
+
+# Random Walk 
+# (Caveat: only if step has zero mean is this actually symmetric! ) 
+struct RandomWalk{D_init <: Distribution, D_step <: Distribution, F<:Function, S} <: SimpleSymmetricProposal{S}
+    init::D_init
+    step::D_step
+    agg::F
+end 
+function RandomWalk(init::Distribution, step::Distribution, agg::Function=(+))
+    Tinit = typeof(rand(init))
+    Tstep = typeof(rand(step))
+
+    # check that aggretation method yields propper type 
+    if ! hasmethod(agg, (Tinit, Tstep)) 
+        throw( ArgumentError("Aggregator $agg does not support ($Tinit, $Tstep) arguments.") )
+    end
 
-# assert that the function can be called with the state
-function LogDensity(S::Type{<:AbstractState}, f::Function) 
-    if hasmethod(f, (S,))
-        return LogDensity{typeof(f),S}(f)
-    else 
-        throw(ArgumentError("The function $f does not accept a $S argument"))
+    # check that type is stable over chain 
+    T = typeof(agg(rand(init), rand(step)))
+    if T != Tinit 
+        throw( ArgumentError("Aggreted type $T differes from initial type $Tinit.") )
     end
+
+    return RandomWalk{typeof(init), typeof(step), typeof(agg), Sample{T}}(init, step, agg)
 end
-LogDensity(T::Type, f::Function) = LogDensity(State{T}, x -> f(x.state))
+propose(rng::AbstractRNG, p::RandomWalk) = Sample(rand(rng, p.init))
+propose(rng::AbstractRNG, p::RandomWalk, x::Sample) = Sample( p.agg(x.sample, rand(rng, p.step)) )
 
+CyclicWalk(a::Real, b::Real, s::Real=(b-a)/2) = RandomWalk(Uniform(a,b), Uniform(0,s), (x,y)->mod(x+y-a,b-a)+a)
 
-# call logdensity with a state or a value, ignore arguments
-(d::AbstractLogDensity)(x; kwargs...) = logdensity(d, x; kwargs...)
-logdensity(d::LogDensity{F,S}, x::S) where {F <: Function, S <: AbstractState} = d.logdensity(x)
 
-## for distribution just wrap the logpdf function
-using Distributions
-LogDensity(d::UnivariateDistribution)     = LogDensity(Float64, x::Float64 -> logpdf(d, x))
-LogDensity(d::MultivariateDistribution)   = LogDensity(Vector{Float64}, x::Vector{Float64} -> logpdf(d, x))
+
+## === LogDensity
+abstract type AbstractLogDensity{S <: AbstractSample} <: AbstractMCMC.AbstractModel end
+Base.eltype(d::AbstractLogDensity{S}) where {S} = S 
+
+""" Define an unnormalized logarithmic probability density 
+    function `f(x)` by implementing  method `evaluate(f, x)`   
+    where `x` is a sample that matches the sample type `eltype(f)`.
+""" 
+function evaluate end
+(d::AbstractLogDensity{S})(x::S; kwargs...) where {S} = evaluate(d, x; kwargs...)
+
+
+## Single-level log density 
+abstract type SimpleLogDensity{S <: SimpleSample} <: AbstractLogDensity{S} end
 
 
-## =============== tests 
-d = LogDensity(Float64, x::Float64 -> x^2)
-d(State(2.0))
+# wrap a distribution, where logpdf is known
+struct DistrLogDensity{D <: Distribution, S <: SimpleSample} <: SimpleLogDensity{S}
+    d::D
+end
+function LogDensity(d::Distribution)
+    S = Sample{typeof(rand(d))}
+    return DistrLogDensity{typeof(d), S}(d)
+end
+evaluate(f::DistrLogDensity{D,S}, x::S) where {D <: Distribution,S <: Sample} = logpdf(f.d, x.sample) 
 
-d = LogDensity(Normal(0,1))
-d(State(2.0)) 
 
-d = LogDensity(MvNormal([0.0, 0.0], [1.0 0.5; 0.5 1.0]))
-d(State([2.0, 2.0]))
\ No newline at end of file
+# wrap a specific method of a function 
+struct FuncLogDensity{F <: Function, S <: SimpleSample} <: SimpleLogDensity{S}
+    f :: F
+end
+function LogDensity(f::Function, T::Type)
+    if ! isa(T, SimpleSample) T = State{T} end 
+
+    if hasmethod(f, (T,)) 
+        return FuncLogDensity{typeof(f), T}(f) 
+
+    elseif hasmethod(f, (eltype(T),))
+        g = x -> f(x.state)
+        return FuncLogDensity{typeof(g), T}(g)
+    else
+        throw(ArgumentError("Function $f does not support $T argument."))
+    end
+end
+evaluate(d::FuncLogDensity{F,S}, x::S) where {F <: Function, S} = d.f(x)

src/hastings.jl

+
+include("base.jl")
+#include("propose.jl")
+
+using AbstractMCMC
+using Distributions
+using Random
+
+struct MetropolisHastings{P <: AbstractProposal} <: AbstractMCMC.AbstractSampler
+    proposal :: P
+end
+
+function AbstractMCMC.step(rng::AbstractRNG, model::AbstractLogDensity, s::MetropolisHastings; kwargs...)
+    sample = initialize(rng, s.proposal)
+    logp = logdensity(model, sample) 
+    state = (sample, logp, true)
+    return sample, state
+end
+
+function AbstractMCMC.step(rng::AbstractRNG, model::AbstractLogDensity, s::MetropolisHastings, state; kwargs...)
+    sample, logp, accept  = state
+    new_sample = propose(s.proposal, sample)
+    new_logp   = logdensity(model, new_sample)
+    q_logratio = logdensityratio(s.proposal, new_sample, sample)
+
+    if log(rand(rng)) < new_logp - logp + q_logratio
+        return new_sample, (new_sample, new_logp, true)
+    else
+        return sample, (sample, logp, false)
+    end
+end
+
+
+## =============== tests ===========================
+#using Test
+#p = StaticProposal(Uniform(-1, 1))
+#s = MetropolisHastings(p)
+#m = LogDensity(Normal(0, 1))
+#@time c = sample(m, s, 100000)
\ No newline at end of file

src/hierarchy.jl

-abstract type Multilevel ... end
\ No newline at end of file
+## Multi-level state 
+abstract type HierarchicalState <: AbstractState end
+
+#=
+# by default, the multilevel state wraps a tuple of states
+struct MultilevelState{T <: Tuple} <: HierarchicalState
+    levels :: T
+end  
+MultilevelState(levels...; kwargs...) = MultilevelState(levels)
+
+course(::MultilevelState) = MultilevelState(s.levels[1:end-1])
+fine(s::MultilevelState)   = s.levels[end]
+combine(s::MultilevelState, fine_state) = MultilevelState(s.levels..., fine_state)
+level(s::MultilevelState) = length(s.levels)
+
+# Multi-level log density
+abstract type HierarchicalLogDensity{S <: HierarchicalState} <: AbstractLogDensity{S} end
+
+struct MultilevelLogDensity{F, S} <: HierarchicalLogDensity{S}
+    logdensities :: F
+end
+
+logdensity(d::HierarchicalLogDensity, state::HierarchicalState; level=length(state)) = logdensity(d, state)
+=#
+#level(s::MultilevelState) = length(s.states)
+#Base.eltype(states::MultilevelState) = eltype(states.states)
+
+## construct from single level states
+#MultilevelState(states::State ... ) = MultilevelState(getindex.([states ... ], :state))
+
+## recursive construction
+#course(states::MultilevelState) = MultilevelState(states[1:end-1])
+#fine(states::MultilevelState)   = states[end]
+#combine(states::MultilevelState, state) = MultilevelState(vcat(states.states ... , state)) 
+
+
+
+#=
+## Multi-level log density
+abstract type HierarchicalLogDensity{S} <: AbstractLogDensity{S} end
+struct MultilevelLogDensity{F, S} <: HierarchicalLogDensity{S}
+    logdensities :: F
+end
+logdensity(d::MultilevelLogDensity, state::MultilevelState) = logdensity(d.logdensities, state.states
+
+
+#logdensity(d::HierarchicalLogDensity, state::HierarchicalState, level::Int)
+
+function logdensity(d::MultilevelLogDensity, state::HierarchicalState, level::Int)
+    logdensity(d.logdensities[level], state.states[level]
+    
+end
+#MultilevelLogDensity(logdensities::AbstractVector{<:AbstractLogDensity}) = MultilevelLogDensity(logdensities
+
+=#
+
+#abstract type AbstractMultilevelLogDensity <: LogDensity{S} end
+
+#abstract type AbstractMultilevelLogDensity <: LogDensity{S} end
+
+#struct MultilevelLogDensity <: AbstractMultilevelLogDensity
+#end
+

test/sanity.jl

+using Test
+using Distributions
+using Statistics
+
+using Revise
+using MultilevelChainSampler
+
+
+function sample_chain(p, n=100)
+    x = [propose(p)] 
+    for i in 1:n
+        push!(x, propose(p, x[end]))
+    end
+    x = value.(x)
+    return x
+end
+
+# Test src/base.jl
+@testset "base" begin
+
+s = Sample(1.0)
+@test value(s) == 1.0
+
+p = StaticProposal(Normal(0,1)) 
+c = sample_chain(p, 10000)
+@test ≈( mean(c), 0.0 , atol=0.1 )
+@test ≈( std(c),  1.0 , atol=0.1 )
+
+p = CyclicWalk(0,1, .1)
+c = sample_chain(p, 10000)
+@test ≈( mean(c) , 0.5 , atol=0.1 ) 
+@test ≈( var(c) , 0.08333333333333333 , atol=0.1 )
+
+
+d = LogDensity(Uniform(0,1))
+@test evaluate(d, Sample(0.0)) == d(Sample(0.0))
+
+end