Merge pull request #25 from yomichi/vectorobservable

update observables

docs/src/lib/public.md

@@ -61,6 +61,25 @@ simple_estimator
 improved_estimator
 ```
 
+## Observables
+```@docs
+MCObservable
+ScalarObservable
+VectorObservable
+MCObservableSet
+makeMCObservable!
+SimpleObservable
+SimpleVectorObservable
+SimpleObservableSet
+SimpleVectorObservableSet
+binning
+Jackknife
+JackknifeVector
+JackknifeSet
+JackknifeVectorSet
+jackknife
+```
+
 ## Snapshots
 ```@docs
 gen_snapshot!

src/lattice/generate.jl

@@ -106,8 +106,8 @@ function generatelattice_std(param)
     bonds = Bond[]
     ib = 0
 
-    use_index_as_sitetype = get(param, "Use Indecies as Site Types", false)
-    use_index_as_bondtype = get(param, "Use Indecies as Bond Types", false)
+    use_index_as_sitetype = get(param, "Use Indicies as Site Types", false)
+    use_index_as_bondtype = get(param, "Use Indicies as Bond Types", false)
 
     for icell in 0:(numcell - 1)
         cellcoord = index2coord(icell, L)

src/observables/MCObservables.jl

@@ -13,8 +13,25 @@ export p_value
 export show, dump_plot
 export merge, merge!
 
+@doc """
+    MCObservable
+
+    abstract type representing observable in Monte Carlo calculation.
+"""
 abstract type MCObservable end
+
+"""
+    ScalarObservable <: MCObservable
+
+    abstract type representing scalar-type observable in Monte Carlo calculation.
+"""
 abstract type ScalarObservable <: MCObservable end
+
+"""
+    VectorObservable <: MCObservable
+
+    abstract type representing vector-type observable in Monte Carlo calculation.
+"""
 abstract type VectorObservable <: MCObservable end
 
 function p_value(X::MCObservable, y::Real; verbose::Bool=false)

src/observables/binning.jl

@@ -13,11 +13,21 @@ mutable struct BinningObservable <: ScalarObservable
     lastbin::Int
     minbinnum::Int
     maxlevel::Int
-end
 
-function BinningObservable(minbinnum::Int=128)
-    return BinningObservable(zeros(0), zeros(0), zeros(1), zeros(1), zeros(Int, 1), 1, 0,
-                             minbinnum, 1)
+    function BinningObservable(minbinnum::Int=128)
+        Base.depwarn("BinningObservable is deprecated. Use obs=SimpleObservable() and binning(obs) instead.",
+                     nothing; force=true)
+
+        raw_ts = zeros(0)
+        bins = zeros(0)
+        sum = zeros(1)
+        sum2 = zeros(1)
+        entries = zeros(Int, 1)
+        binsize = 1
+        lastbin = 0
+        maxlevel = 1
+        return new(raw_ts, bins, sum, sum2, entries, binsize, lastbin, minbinnum, maxlevel)
+    end
 end
 
 function reset!(b::BinningObservable)

src/observables/binningvector.jl

@@ -4,21 +4,30 @@ export BinningVectorObservableSet
 
 mutable struct BinningVectorObservable <: VectorObservable
     ## index of these vectors denotes the level of bins (each bin stores the mean of 2^(i-1) values)
-    raw_ts::Vector{Vector{Float64}}        ## Time series of raw data
-    bins::Vector{Vector{Float64}}          ## Time series of bins
-    sum::Vector{Vector{Float64}}           ## Summation of stored values
-    sum2::Vector{Vector{Float64}}          ## Summation of square of bin's mean
-    entries::Vector{Int}           ## Number of bins
+    raw_ts::Vector{Vector{Float64}} ## Time series of raw data
+    bins::Vector{Vector{Float64}}   ## Time series of bins
+    sum::Vector{Vector{Float64}}    ## Summation of stored values
+    sum2::Vector{Vector{Float64}}   ## Summation of square of bin's mean
+    entries::Vector{Int}            ## Number of bins
     binsize::Int
     lastbin::Int
     minbinnum::Int
     maxlevel::Int
-end
 
-function BinningVectorObservable(minbinnum::Int=128)
-    return BinningVectorObservable(Vector{Float64}[], Vector{Float64}[], Vector{Float64}[],
-                                   Vector{Float64}[],
-                                   zeros(Int, 1), 1, 0, minbinnum, 1)
+    function BinningVectorObservable(minbinnum::Int=128)
+        Base.depwarn("BinningVectorObservable is deprecated. Use obs=SimpleVectorObservable() and binning(obs) instead.",
+                     nothing; force=true)
+
+        raw_ts = Vector{Float64}[]
+        bins = Vector{Float64}[]
+        sum = Vector{Float64}[]
+        sum2 = Vector{Float64}[]
+        entries = zeros(Int, 1)
+        binsize = 1
+        lastbin = 0
+        maxlevel = 1
+        return new(raw_ts, bins, sum, sum2, entries, binsize, lastbin, minbinnum, maxlevel)
+    end
 end
 
 function reset!(b::BinningVectorObservable)

src/observables/jackknife.jl

@@ -1,13 +1,41 @@
 export Jackknife, JackknifeSet
 export jackknife
 
+"""
+    Jackknife <: ScalarObservable
+
+    Jackknife resampling observable.
+"""
 mutable struct Jackknife <: ScalarObservable
     xs::Vector{Float64}
 end
 
-Jackknife(jk::Jackknife, f::Function) = Jackknife(map(f, jk.xs))
+"""
+    Jackknife(f::Function, jks::Jackknife...)
+
+    Construct a Jackknife observable by applying `f` to `jks`.
+    For example, `Jackknife(mean, jk1, jk2, jk3)` returns a Jackknife observable of the means of `jk1`, `jk2`, and `jk3`.
+"""
+function Jackknife(f::Function, jks::Jackknife...)
+    if isempty(jks)
+        return Jackknife(zeros(0))
+    end
+    nbins = count(jks[1])
+    xs = zeros(nbins)
+    for i in 1:nbins
+        arg = [jk.xs[i] for jk in jks]
+        xs[i] = f(arg...)
+    end
+    return Jackknife(xs)
+end
+
+function Jackknife(jk::Jackknife, f::Function)
+    Base.depwarn("Jackknife(jk::Jackknife, f::Function) is deprecated. Use Jackknife(f, jk) instead.",
+                 :Jackknife)
+    return Jackknife(f, jk)
+end
 
-function jk_helper(xs)
+function jk_helper(xs::Vector{Float64})
     s = sum(xs)
     n = length(xs) - 1
     ret = similar(xs)
@@ -33,30 +61,32 @@ end
 count(jk::Jackknife) = length(jk.xs)
 
 mean(jk::Jackknife) = mean(jk.xs)
-function stderror(jk::Jackknife)
-    n = count(jk)
-    if n < 2
-        return NaN
-    else
-        m2 = sum(abs2, jk.xs)
-        m2 /= n
-        m = mean(jk)
-        sigma2 = m2 - m * m
-        sigma2 *= n - 1
-        sigma2 = maxzero(sigma2)
-        return sqrt(sigma2)
-    end
-end
+# function stderror(jk::Jackknife)
+#     n = count(jk)
+#     if n < 2
+#         return NaN
+#     else
+#         m2 = sum(abs2, jk.xs)
+#         m2 /= n
+#         m = mean(jk)
+#         sigma2 = m2 - m*m
+#         sigma2 *= n-1
+#         sigma2 = maxzero(sigma2)
+#         return sqrt(sigma2)
+#     end
+# end
 function var(jk::Jackknife)
     n = count(jk)
     if n < 2
         return NaN
     else
         m = mean(jk)
-        return sum(abs2, jk.xs .- m) * (n - 1)
+
+        return mean(abs2, jk.xs .- m) * (n - 1)
     end
 end
 stddev(jk::Jackknife) = sqrt(var(jk))
+stderror(jk::Jackknife) = stddev(jk)
 
 function confidence_interval(jk::Jackknife, confidence_rate::Real)
     q = 0.5 + 0.5 * confidence_rate
@@ -90,16 +120,15 @@ unary_functions = (:-,
                    :sqrt, :cbrt)
 
 for op in unary_functions
-    @eval Base.$op(jk::Jackknife) = Jackknife(jk, $op)
+    @eval Base.$op(jk::Jackknife) = Jackknife($op, jk)
 end
 
 binary_functions = (:+, :-, :*, :/, :\)
 
 for op in binary_functions
-    @eval Base.$op(jk::Jackknife, rhs::Real) = Jackknife(jk, lhs -> ($op)(lhs, rhs))
-    @eval Base.$op(lhs::Real, jk::Jackknife) = Jackknife(jk, rhs -> ($op)(lhs, rhs))
-    @eval Base.$op(lhs::Jackknife, rhs::Jackknife) = Jackknife(broadcast($op, lhs.xs,
-                                                                         rhs.xs))
+    @eval Base.$op(jk::Jackknife, rhs::Real) = Jackknife(lhs -> ($op)(lhs, rhs), jk)
+    @eval Base.$op(lhs::Real, jk::Jackknife) = Jackknife(rhs -> ($op)(lhs, rhs), jk)
+    @eval Base.$op(lhs::Jackknife, rhs::Jackknife) = Jackknife($op, lhs, rhs)
 end
 
 import Base.^
@@ -109,9 +138,25 @@ import Base.^
 ^(lhs::Integer, rhs::Jackknife) = Jackknife(lhs .^ (rhs.xs))
 ^(lhs::Jackknife, rhs::Jackknife) = Jackknife((lhs.xs) .^ (rhs.xs))
 
+"""
+    JackknifeSet
+
+    Alias of `MCObservableSet{Jackknife}`.
+"""
 const JackknifeSet = MCObservableSet{Jackknife}
 
+@doc doc"""
+    jackknife(obs::ScalarObservable)
+
+    Construct a Jackknife observable from a scalar observable
+"""
 jackknife(obs::ScalarObservable) = Jackknife(obs)
+
+@doc doc"""
+    jackknife(obsset::MCObservableSet)
+
+    Construct a JackknifeSet from a MCObservableSet
+"""
 function jackknife(obsset::MCObservableSet)
     JK = JackknifeSet()
     for (k, v) in obsset