Rework BHM-example.jl,

new function `default_starting_vector`

docs/make.jl

@@ -44,7 +44,11 @@ end
 
 makedocs(;
     modules=[Rimu,Rimu.RimuIO],
-    format=Documenter.HTML(prettyurls = false),
+    format=Documenter.HTML(
+        prettyurls = false,
+        size_threshold=500_000, # 500 kB
+        size_threshold_warn=200_000, # 200 kB
+    ),
     pages=[
         "Guide" => "index.md",
         "Examples" => EXAMPLES_PAIRS[sortperm(EXAMPLES_NUMS)],
@@ -69,7 +73,7 @@ makedocs(;
     sitename="Rimu.jl",
     authors="Joachim Brand <j.brand@massey.ac.nz>",
     checkdocs=:exports,
-    doctest=false # Doctests are done while testing.
+    doctest=false, # Doctests are done while testing.
 )
 
 deploydocs(

scripts/BHM-example.jl

@@ -11,103 +11,109 @@
 # `Rimu` for FCIQMC calculation;
 
 using Rimu
-using Random
+using Plots # for plotting
+
+# ## Setting up the model
 
 # Now we define the physical problem:
-# Setting the number of lattice sites `m = 6`;
-# and the number of particles `n = 6`:
-m = n = 6
-# Generating a configuration that particles are evenly distributed:
-aIni = near_uniform(BoseFS{n,m})
-# where `BoseFS` is used to create a bosonic system.
-# The Hamiltonian is defined based on the configuration `aIni`,
-# with additional onsite interaction strength `u = 6.0`
-# and the hopping strength `t = 1.0`:
+# Generating a configuration where 6 particles are evenly distributed in 6 lattice sites:
+aIni = near_uniform(BoseFS{6,6})
+# where `BoseFS` is used to create a bosonic Fock state.
+# The Hamiltonian is defined by specifying the model and the parameters.
+# Here we use the Bose Hubbard model in one dimension and real space:
 Ĥ = HubbardReal1D(aIni; u = 6.0, t = 1.0)
 
+# ## Parameters of the calculation
+
+# Now let's setup the Monte Carlo calculation.
+# We need to decide the number of walkers to use in this Monte Carlo run, which is
+# equivalent to the average one-norm of the coefficient vector:
+targetwalkers = 1_000;
+
+# It is good practice to equilibrate the time series before taking statistics.
+steps_equilibrate = 1_000;
+steps_measure = 2_000;
+
+# The appropriate size of the time step is problem dependent.
+dτ = 0.001;
+
+# ## Defining an observable
+
+# Now let's set up an observable to measure. Here we will measure the projected energy.
+# In additon to the `shift`, the projected energy is a second estimator for the energy.
+
+# We first need to define a projector. Here we use the function `default_starting_vector`
+# to generate a vector with only a single occupied configuration, the same vector that we
+# will use as starting vector for the FCIQMC calculation.
+svec = default_starting_vector(aIni; style=IsStochasticInteger())
+# The choice of the `style` argument already determines the FCIQMC algorithm to use
+# (while it is irrelevant for the projected energy).
 
-# Now let's setup the Monte Carlo settings.
-# The number of walkers to use in this Monte Carlo run:
-targetwalkers = 1_000
-# The number of time steps before doing statistics,
-# i.e. letting the walkers to sample Hilbert and to equilibrate:
-steps_equilibrate = 1_000
-# And the number of time steps used for getting statistics,
-# e.g. time-average of shift, projected energy, walker numbers, etc.:
-steps_measure = 2_000
-
-
-# Set the size of a time step
-dτ = 0.001
-# and we report QMC data every k-th step,
-# set the interval to record QMC data:
-reporting_interval = 1
-
-# Now we prepare initial state and allocate memory.
-# The initial address is defined above as `aIni = near_uniform(Ĥ)`.
-# Putting one of walkers into the initial address `aIni`
-svec = DVec(aIni => 1)
-# Let's plant a seed for the random number generator to get consistent result:
-Random.seed!(17)
-
-# Now let's setup all the FCIQMC strategies.
-
-# Passing dτ and total number of time steps into params:
-params = RunTillLastStep(dτ = dτ, laststep = steps_equilibrate + steps_measure)
-# Strategy for updating the shift:
-s_strat = DoubleLogUpdate(targetwalkers = targetwalkers, ζ = 0.08)
-# Strategy for reporting info:
-r_strat = ReportDFAndInfo(reporting_interval = reporting_interval, info_interval = 100)
-# Strategy for updating dτ:
-τ_strat = ConstantTimeStep()
-# set up the calculation and reporting of the projected energy
-# in this case we are projecting onto the starting vector,
-# which contains a single configuration
-post_step = ProjectedEnergy(Ĥ, copy(svec))
-
-# Print out info about what we are doing:
-println("Finding ground state for:")
-println(Ĥ)
-println("Strategies for run:")
-println(params, s_strat)
-println(τ_strat)
+# Observables are passed into the `lomc!` function with the `post_step` keyword argument.
+post_step = ProjectedEnergy(Ĥ, svec)
+
+# ## Running the calculation
+
+# Seeding the random number generator is sometimes useful in order to get reproducible
+# results
+using Random
+Random.seed!(17);
 
 # Finally, we can start the main FCIQMC loop:
-df, state = lomc!(Ĥ,svec;
-            params,
-            laststep = steps_equilibrate + steps_measure,
-            s_strat,
-            r_strat,
-            τ_strat,
+df, state = lomc!(Ĥ, svec;
+            laststep = steps_equilibrate + steps_measure, # total number of steps
+            dτ,
+            targetwalkers,
             post_step,
 );
 
-# Here is how to save the output data stored in `df` into a `.arrow` file,
-# which can be read in later:
-println("Writing data to disk...")
-save_df("fciqmcdata.arrow", df)
+# `df` is a `DataFrame` containing the time series data.
+
+# ## Analysing the results
 
-# Now let's look at the calculated energy from the shift.
-# Loading the equilibrated data
-qmcdata = last(df,steps_measure);
+# We can plot the norm of the coefficient vector as a function of the number of steps:
+hline([targetwalkers], label="targetwalkers", color=:red, linestyle=:dash)
+plot!(df.steps, df.norm, label="norm", ylabel="norm", xlabel="steps")
+# After some equilibriation steps, the norm fluctuates around the target number of walkers.
 
-# compute the average shift and its standard error
-se = shift_estimator(qmcdata)
+# Now let's look at estimating the energy from the shift.
+# The mean of the shift is a useful estimator of the shift. Calculating the error bars
+# is a bit more involved as correlations have to be removed from the time series.
+# The following code does that:
+se = shift_estimator(df; skip=steps_equilibrate)
 
 # For the projected energy, it a bit more complicated as it's a ratio of two means:
-pe = projected_energy(qmcdata)
+pe = projected_energy(df; skip=steps_equilibrate)
 
 # The result is a ratio distribution. Let's get its median and lower and upper error bars
 # for a 95% confidence interval
 v = val_and_errs(pe; p=0.95)
 
+# Let's visualise these estimators together with the time series of the shift
+plot(df.steps, df.shift, ylabel="energy", xlabel="steps", label="shift")
+
+plot!(x->se.mean, df.steps[steps_equilibrate+1:end], ribbon=se.err, label="shift mean")
+plot!(
+    x -> v.val, df.steps[steps_equilibrate+1:end], ribbon=(v.val_l,v.val_u),
+    label="projected_energy",
+)
+# In this case the projected energy and the shift are close to each other an the error bars
+# are hard to see on this scale.
+
+# The problem was just a toy example, as the dimension of the Hamiltonian is rather small:
+dimension(Ĥ)
+
+# In this case is easy (and more efficient) to calculate the exact ground state energy
+# using standard linear algebra:
+using LinearAlgebra
+exact_energy = eigvals(Matrix(Ĥ))[1]
+
+# Comparing our results for the energy:
 println("Energy from $steps_measure steps with $targetwalkers walkers:
-         Shift: $(se.mean) ± $(se.err);
-         Projected Energy: $(v.val) ± ($(v.val_l), $(v.val_u))")
+         Shift: $(se.mean) ± $(se.err)
+         Projected Energy: $(v.val) ± ($(v.val_l), $(v.val_u))
+         Exact Energy: $exact_energy")
 
-# Finished!
 
 using Test                                      #hide
-@test isfile("fciqmcdata.arrow")                #hide
-@test se.mean ≈ -4.0215 rtol=0.1                #hide
-rm("fciqmcdata.arrow", force=true)              #hide
+@test se.mean ≈ -4.0215 rtol=0.1;               #hide

src/Rimu.jl

@@ -49,6 +49,7 @@ include("StatsTools/StatsTools.jl")
 @reexport using .StatsTools
 
 export lomc!
+export default_starting_vector
 export FciqmcRunStrategy, RunTillLastStep
 export ShiftStrategy, LogUpdate, LogUpdateAfterTargetWalkers
 export DontUpdate, DoubleLogUpdate, DoubleLogUpdateAfterTargetWalkers

src/lomc.jl

@@ -279,6 +279,43 @@ function report_default_metadata!(report::Report, state::QMCState)
     return report
 end
 
+"""
+    default_starting_vector(hamiltonian::AbstractHamiltonian; kwargs...)
+    default_starting_vector(
+        address=starting_address(hamiltonian);
+        style=IsStochasticInteger(),
+        threading=nothing
+    )
+Return a default starting vector for [`lomc!`](@ref). The default choice for the starting
+vector is
+```julia
+v = PDVec(starting_address => 10; style)
+```
+if threading is available or
+
+```julia
+v = DVec(starting_address => 10; style)
+```
+otherwise. See [`PDVec`](@ref), [`DVec`](@ref) and [`StochasticStyle`](@ref).
+"""
+function default_starting_vector(hamiltonian::AbstractHamiltonian; kwargs...)
+    return default_starting_vector(starting_address(hamiltonian); kwargs...)
+end
+function default_starting_vector(address;
+    style=IsStochasticInteger(),
+    threading=nothing,
+)
+    if isnothing(threading)
+        threading = Threads.nthreads() > 1
+    end
+    if threading
+        v = PDVec(address => 10; style)
+    else
+        v = DVec(address => 10; style)
+    end
+    return v
+end
+
 """
     lomc!(ham::AbstractHamiltonian, [v]; kwargs...) -> df, state
     lomc!(state::QMCState, [df]; kwargs...) -> df, state
@@ -363,16 +400,8 @@ julia> metadata(df2, "hamiltonian") # some metadata is automatically added
   be passed for merging with the report `df`.
 
 The default choice for the starting vector is
-```julia
-v = PDVec(starting_address => 10; style)
-```
-if threading is available or
-
-```julia
-v = DVec(starting_address => 10; style)
-```
-otherwise. It triggers the integer walker FCIQMC algorithm. See [`PDVec`](@ref),
-[`DVec`](@ref) and [`StochasticStyle`](@ref).
+`v = default_starting_vector(; address, style, threading)`.
+See [`default_starting_vector`](@ref).
 """
 function lomc!(ham, v; df=DataFrame(), name="lomc!", metadata=nothing, kwargs...)
     state = QMCState(ham, v; kwargs...)
@@ -385,6 +414,7 @@ function lomc!(
     address=starting_address(ham),
     kwargs...
 )
+    v = default_starting_vector(ham; style, threading, address)
     if Threads.nthreads() > 1 && (threading ≠ false) || threading == true
         v = PDVec(address => 10; style)
     else

src/strategies_and_params/poststepstrategy.jl

@@ -24,8 +24,10 @@ abstract type PostStepStrategy end
 """
     post_step(::PostStepStrategy, ::ReplicaState) -> kvpairs
 
-Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs. 
-This function is only called every [`reporting_interval`](@ref) steps, as defined by the `ReportingStrategy`. 
+Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs.
+This function is only called every [`reporting_interval`](@ref) steps, as defined by the
+`ReportingStrategy`.
+
 See also [`PostStepStrategy`](@ref), [`ReportingStrategy`](@ref).
 """
 post_step