Simplify lomc! arguments (#223)

* avoid `params`

* fix tests, kwarg address

* doc changes

* more doc changes

* Rework BHM-example.jl,
new function `default_starting_vector`

* add tests

* fix tests

* improve docs

* simplify code

* remove params from G2-example.jl

---------

Co-authored-by: Joachim Brand <joachim.brand@gmail.com>

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,112 @@
 # `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 with blocking transformations:
+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)
+# For the projected energy, it a bit more complicated as it's a ratio of fluctuationg
+# quantities:
+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]
+# Read more about `Rimu.jl`s capabilities for exact diagonalisation in the example
+# "Exact diagonalisation".
+
+# 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

scripts/G2-example.jl

@@ -53,21 +53,14 @@ targetwalkers = 100;
 
 # Other FCIQMC parameters and strategies are the same as before
 dτ = 0.001
-reporting_interval = 1
 svec = DVec(aIni => 1)
-Random.seed!(17)
-params = RunTillLastStep(dτ = dτ, laststep = steps_equilibrate + steps_measure)
-s_strat = DoubleLogUpdate(targetwalkers = targetwalkers, ζ = 0.08)
-r_strat = ReportDFAndInfo(reporting_interval = reporting_interval, info_interval = 100)
-τ_strat = ConstantTimeStep();
+Random.seed!(17);
 
 # Now we run the main FCIQMC loop:
 df, state = lomc!(H, svec;
-            params,
+            dτ,
             laststep = steps_equilibrate + steps_measure,
-            s_strat,
-            r_strat,
-            τ_strat,
+            targetwalkers,
             replica,
 );
 

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