carstenbauer / montecarlo.jl Goto Github PK

Currently DQMC (and MC) is the root type for a simulation, meaning it holds a bunch stuff together while also representing the algorithm used. This leads to some problems of availability, i.e. in a couple of places we want to dispatch on the algorithm, but it's not available (we dispatch on Type{<: MonteCarloFlavor} instead) and in some others we need parameters that have not yet been generated.

I've thought about revamping things a bit in #107 and I wanted to collect some of those thoughts here.

Currently existing structures:

Lattice

The lattice is self contained and can be created with just a lattice size and lattice type. All good.

Model

The model currently holds the lattice but only uses it create some temporary storage. Maybe this could be moved into an init!() function and become self contained? Difficulty here is eltype.

Measurements

Measurements are kind of a mess right now. I want to avoid re-running time displaced greens iterators and duplicating lattice iterators which has made things quite complicated. Both however are constructed when the simulation is ran, so they should be fairly detached. What isn't detached is the recorder, i.e. LogBinner or FullBinner from BinningAnalysis, as those are created with a zero-element. Maybe we could add some functionality to BinningAnalysis to resize those at a later point, but that leaves an eltype problem.

LatticeIterator

Requires the lattice and sometimes the model to know the number of flavors.

Updates

Are fully detached at this point.

Update scheduler

Only requires updates on construction.

Configuration recorder

Requires type information of the compressed configuration.

(Unequal time) Stack

Requires a bunch of type information (for greens, hopping and interaction matrices) from model + algorithm. Needs number of sites, flavors, time slices and safe_mult on initialization. Conceptually part of the algorithm.

New or revamped structures

Algorithm / MonteCarloFlavor / Engine

I think it would be good to bundle up things that are specifically needed by the algorithm. For DQMC this would be the stack, unequal time stack, some parameters and some analysis variables.

(Monte Carlo) Simulation

This would be the outer most wrapper which brings everything together. It would contain:

• Some (or all?) parameters.
• The configuration as both MC and DQMC have it
• A configuration recorder
• Measurements
• the algorithm / Monte Carlo flavor / engine or whatever we may call it
• the model if we can detach it from the algorithm
• the update scheduler
• the lattice
• perhaps some more analysis tools

Parameters

Parameters are a bit weird, conceptually they should be attached to the types that need them but it might also be good to have them all in one place. Sorting the parameters we have to structs they belong to:

• Lattice: linear system size L (, number of sites N)
• Model: model parameters U, t, µ etc
• DQMC: performed checks (sign problem, propagation error), safe_mult, delta_tau, number of time slices
• Simulation: number of thermalization & measurement sweeps, measure rate, print rate, saving options (checkpoint, early save, filename, overwrite, rename)
• Recorder: recording rate (this should probably be detached from measure rate)
• Measurements (measure rate?)
• Scheduler: (types and order of updates, adaptive configuration)

It might also be interesting to make "Simulation recipes" where you can define parameters (+ some types like model, lattice, etc) in a Dict-like fashion and auto-complete missing ones. I could see something like this working well with what I messed around with in https://github.com/ffreyer/ParameterFiles.jl to generate a set of parallel simulations.

A useful concept perhaps is a deferred creation of structs. One could generate a Specification(constructor, args, kwargs) and later call constructor(extra_args, args...; kwargs...) to create the object.

Currently one can pass measurements as measurements = :default, which calls default_measurements(mc, model) or as a Dict{Symbol, AbstractMeasurements}. All the default measurements are however constructed as Measurement(mc, model) meaning they can't be constructed before calling MC or DQMC. It would make sense to allow passing a dictionary of types, which then internally call their constructors.

It would be good to test lattice iterators and LatPhysLattice, but I'm not sure how with LatticePhysics being an unregistered add-on.

I've noticed that currently, there are only a few points left where MC and DQMC explicitly interface with the lattice. Other than in build_checkerboard(lattice) only a few calls to length(lattice) (or m.l.sites) remain. We should remove these in favor of nsites(model).

Maybe we should also rethink build_checkerboard(), since it is the only thing giving requirements to the lattice when writing a new model.

Since T is a static matrix we should not be repeating this computation. Same goes for the inverse.

https://github.com/crstnbr/MonteCarlo.jl/blob/ab564454400ea9d921153117875a37710279beac/src/flavors/DQMC/slice_matrices.jl#L17-L21

https://github.com/crstnbr/MonteCarlo.jl/blob/ab564454400ea9d921153117875a37710279beac/src/flavors/DQMC/slice_matrices.jl#L31-L39

Quick benchmark using my own model w/ 16 sites and 20 slices:

                                                       Time                   Allocations      
                                               ──────────────────────   ───────────────────────
               Tot / % measured:                  535495s / 0.00%           5.17GiB / 26.9%    

Section                                ncalls     time   %tot     avg     alloc   %tot      avg
───────────────────────────────────────────────────────────────────────────────────────────────
run!                                        1    3.65s   100%   3.65s   1.39GiB  100%   1.39GiB
...
 slice_matrix!   (precomputed)           23.6k    268ms  8.04%  11.3μs    377KiB  0.03%    16.3B
 slice_matrix!   (not precomputed)       23.6k    513ms  14.1%  21.7μs    377KiB  0.03%    16.3B
───────────────────────────────────────────────────────────────────────────────────────────────

I'm worried that greens() is returning I - G^T instead of G ...

#6 currently includes a function save_measurements!(mc, filename) which calls a function save!(::AbstractMeasurement) on each measurement, saving their observables to some JLD file filename.

It would be good to have a more general save!(mc), which saves MC/DQMC related parameters, model related parameters and measurements. Saving models could be done similarly to measurements, i.e. we define a default save!(::Model) and allow users to to define their own method.

Currently something like DQMC(model, beta=2.2, delta_tau=0.1, safe_mult=10) fails, because we require that N * safe_mult * delta_tau == beta. We can simplify this requirement to N * delta_tau == beta if we either

• perform an irregular stabilization at the first (last?) time propagation or
• disconnect time steps from stabilization

Currently we have

https://github.com/crstnbr/MonteCarlo.jl/blob/1c931b57588f44a2ddfe3b547e01aae3332fc221/src/flavors/MC/MC.jl#L224

https://github.com/crstnbr/MonteCarlo.jl/blob/1c931b57588f44a2ddfe3b547e01aae3332fc221/src/flavors/MC/MC.jl#L228

i.e. ΔE, Δsite swaps to Δsite, ΔE when given to accept_local. We should keep the same order.

We should allow users to specify a particle density instead of chemical potential in DQMC.

Use 4 spaces for 1 tab consistently.

Also, make sure that (almost) no line exceeds the 100 characters ruler.

I want lattice and greens iterators to be created just once since they can be rather large. In the case of CombinedGreensIterator it's also very worth it to group measurements that use it.

The current implementation achieves this by saving types instead of objects in measurements and later resolving those to types. This works but is kind of awkward and complicated. With the addition of wrapper types for lattice iterators this has become even more awkward. Some of those wrappers require data from user input, leading to half-constructed types.

The alternative would be to use a global lookup array or perhaps dictionary. Conceptually we would want something like lookup = Dict(hash(dqmc, lattice_iterator_type) => lattice_iterator), though hash(dqmc) would change as measurements are added.

Afterwards, install CompatHelper.

Specifically @timeit_debug which, by default, is optimized out.

It keeps happening ... inconsistently ... Why? :(

See:
https://github.com/crstnbr/MonteCarlo.jl/runs/1446792885
...
https://github.com/crstnbr/MonteCarlo.jl/runs/1486481645
...
https://github.com/crstnbr/MonteCarlo.jl/runs/1813601955

10_000 + 20_000 sweeps seem to be enough with the current error tolerance

Can we use udt! instead of udt here?
https://github.com/crstnbr/MonteCarlo.jl/blob/310391725eb0607d8d40c4cc4e10217a0d19a5a2/src/flavors/DQMC/stack.jl#L242

Since we have multiple interfaces it would be a good idea to define one consistent InterfaceError <: Exception.

See https://github.com/JuliaLang/julia/blob/4111609c2b49372d5139c0eb9119b738bf4d73d5/base/boot.jl#L293 and https://github.com/JuliaLang/julia/blob/c6da87ff4bc7a855e217856757ad3413cf6d1f79/base/errorshow.jl#L172 for MethodError

At some point we really have to update the docs!

Same as #15 but for models.

Currently the observables are :conf, :Greens and :BosonEnergy. We should document what these are and what the actual data they contain means. In particular, I'm not sure how to interpret the matrix returned by e.g.

julia> meas[:Greens].obs |> mean
8×8 Array{Float64,2}:
  0.486632   -0.307531    -0.0152945   0.0458478   0.0299359    0.0384947  -0.0139994  -0.25898  
 -0.310227    0.538009    -0.304965    0.0196503   0.049056    -0.0135113   0.0716403   0.0193743
 -0.0263509  -0.269518     0.518684   -0.237224   -0.0447043    0.0557874   0.0241328   0.0390698
  0.0654728   0.00375744  -0.317623    0.409047   -0.331587     0.0385827   0.042445   -0.0055775
 -0.018052    0.0400457    0.0192787  -0.237427    0.605781    -0.309867    0.0151516   0.0356679
  0.032978    0.00207081   0.0471798   0.0185239  -0.247373     0.357371   -0.261951    0.0109985
  0.0188517   0.0305705   -0.0226769   0.041055    0.00292663  -0.342347    0.561436   -0.235713 
 -0.335868    0.0088788    0.0530074  -0.0274625   0.0661984    0.0356364  -0.361961    0.438607 

I suppose that one dimension is the imaginary time and another is the spatial coordinate, but I'm not sure which is which.

Use Int64 to be safe from overflows.

We can implement an optimal version (using Bits) later.

A slice matrix should be given by B_l = exp(- 0.5 * delta_tau * T) * exp(- delta_tau * V_l). If I'm not mistaken exp(- delta_tau * V_l) transforms to exp(- lambda * conf).

With that, the sign of power should either change in multiply_slice_matrix..., slice_matrix or interaction_matrix_exponential!, and in ΔE_boson = -2. * lambda * conf[i, slice] should lose the - in propose_local.

As far as I can tell though, this error just sends conf -> -conf without causing other errors.

Create input xml and parse it.

Could use https://github.com/crstnbr/dqmc/blob/master/lattices/gen_square_lattice_xml.jl.

Update to Documenter 0.24.

Since we dropped Travis we have to deploy docs via a github action.

Although <:AbstractLattice is allowed as a type parameter it doesn't currently work for arbitrary lattices.

Among others, rand is assuming a cubic lattice.

julia> m = IsingModel(L=3, dims=2, l=l, neighs=zeros(Int, 0,0))
2D-Ising model, L=3 (36 sites)

julia> mc = MC(m, T=0.1234)
Monte Carlo simulation
Model: 2D-Ising model, L=3 (36 sites)
Beta: 8.1 (T ≈ 0.123)

Check where we can drop mutable.

Currently calling greens(dqmc) always performs calculations. Since most measurements require the the greens function, the calculation is repeated frequently. It would be good to save the result and only recalculate it when necessary.

It might be good to also add a function for I - G and give it the same behavior.

Create/update MC and DQMC constructors for Dicts and NamedTuples.

Because of all the LoopVectorization stuff I added...

You can currently create am attractive Hubbard model with any lattice, but the docstring and the constructors do nott reflect that.

After JuliaRegistries/General#2640 is merged, the installation instructions simplify to

] add https://github.com/crstnbr/MonteCarlo.jl

(both StableDQMC and MonteCarloObservable are registered then)

I don't think we use any functionality from MonteCarloObservable that isn't in BinningAnalysis, so we should switch to just requiring BinningAnalysis.

Should allow us to speed up mul!(C, A, B), for example.

StructArrays.jl might also be useful for complex matrices (see https://youtu.be/qz2kJdVDWi0?t=900)

We use the symmetric Trotter approximation to factorize propagators B ~ exp(T + V) = exp(T/2) exp(V) exp(T/2). An update in DQMC changes exp(V) -> exp(V') = (I + delta) exp(V). In order to write the determinant ratio in terms of the Greens function, we need to be able to isolate (1 + delta), i.e we need to have (I+delta) B or B (I + delta). For this, we introduce an effective Greens function with B ~ exp(T/2) exp(T/2) exp(V), which can be mapped back to the normal Greens function by applying G_normal = exp(-T/2) G exp(T/2). This allows us to pull (I + delta) to the right if both exp(V) and (I + delta) are diagonal. Otherwise we cannot permute them and we're stuck with exp(T/2) exp(T/2) (I + delta) exp(V).

If we were to change the effective Greens function to use B = exp(V) exp(T/2) exp(T/2) this wouldn't be a problem. The normal Greens function would now be G_normal = exp(T/2) G exp(-T/2) and MC updates would need to use G = G - G (I + delta (I - G))^-1 delta (I-G)) instead of G = G - (I -G) (I + delta (I - G))^-1 delta G).

This is because the udt products are not finalized in calculate_slice_matrix_chain and calculate_slice_matrix_chain_dagger. To fix this, repeat

 U *= spdiagm(0 => D)
 U, D, Tnew = decompose_udt(U)
 T =  Tnew * T
 svs[:,svc] = log.(D)
 svc += 1

after the loop.

Hello Carsten,
I am about to start writing MC simulations to find the shape of a surface that minimizes a given energy ( such as bending energy or free energy of contained liquid ).
I wonder how much this fits into the scope of this package and what exactly has to be done to make this work.
As far as I see, I would need to implement at least a non-static lattice and setup a corresponding model type. Would that be all?

This keeps getting in the way when making changes to existing structs. We could just as easily save :DQMC and dispatch on Val(:DQMC) for file loading. (Or use if-statements if dispatch is slow)

With that change we can probably also clean up the DQMC constructors to some degree. If we keep the complete() / make_concrete() functions, we can probably also make a macro to generate them with the relevant constructors.

Except for maybe the configurations, we do not need to store the full time series.

I.e. this file: https://github.com/crstnbr/MonteCarlo.jl/blob/master/src/inplace_udt.jl

Simulate attractive Hubbard model as in https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.161120 and show s-wave SC phase by (3D) plotting equal-time pairing correlations.

Came up in #47.

See CI log here:
https://github.com/crstnbr/MonteCarlo.jl/runs/316807579

Excerpt:

[ Info: Running DQMC β=1.0, 10k + 20k sweeps, ≈1min
[ Info: Running ED
Attractive Hubbard Model (ED): Test Failed at D:\a\MonteCarlo.jl\MonteCarlo.jl\test\ED\ED_tests.jl:93
  Expression: isapprox(G_DQMC[i, j], G_ED[i, j], atol=0.025, rtol=0.1)
   Evaluated: isapprox(0.1197478321658655, 0.09341371846249304; atol=0.025, rtol=0.1)
Stacktrace:
 [1] top-level scope at D:\a\MonteCarlo.jl\MonteCarlo.jl\test\ED\ED_tests.jl:93
 [2] top-level scope at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\Test\src\Test.jl:1113
 [3] top-level scope at D:\a\MonteCarlo.jl\MonteCarlo.jl\test\ED\ED_tests.jl:74

@ffreyer Can you check/fix this?

https://github.com/crstnbr/MonteCarlo.jl/blob/1c931b57588f44a2ddfe3b547e01aae3332fc221/src/flavors/MC/MC.jl#L229

Instead of this it would be better to have acc_rate += 1 and later divide by length(conf). Because

• we don't sweep through sites, we sweep through the configuration
• we can remove nsites() from the model requirements
• it's (slightly) faster

What to compare to? Maybe https://juliaphysics.github.io/PhysicsTutorials.jl/tutorials/machine_learning/ml_ising/ml_ising.html?

In the process of getting a feel for using the package, I tried running the following script:

using MonteCarlo, Distributions, DataFrames

Tdist = Normal(MonteCarlo.IsingTc, .64)
n_Ts = 2^5 # 2^8
Ts = sort!(rand(Tdist, n_Ts))
Ts = Ts[Ts.>=1.2]
Ts = Ts[Ts.<=3.8]
therm = 5*10^1  # 10^4
sweeps = 5*10^1 # 10^3

df = DataFrame(L=Int[], T=Float64[], BE=Float64[])

for L in 2 .^ [3, 4, 5] 
    println("L = ", L)
    for (i, T) in enumerate(Ts)
    	println("\t T = ", T)
	beta = 1/T

        model = HubbardModelAttractive(dims=2, L=L)
        dqmc = DQMC(model, beta=beta, delta_tau = (beta/100), slices=100, checkerboard=true)
        run!(dqmc, sweeps=sweeps, thermalization=therm, verbose=false)

	meas = MonteCarlo.measurements(dqmc)[:ME]
	push!(df, [L, T, mean(meas[:BosonEnergy].obs)])    
    end
end

and found that for L = 2^3 and L = 2^4, I was getting 98% CPU utilization almost constantly across 8 threads, but for L = 2^5, CPU utilization dropped down to about 10% per thread, very occasionally spiking up to 50%.

This is not a very advanced analysis and I'll soon try and dig into what's going on to see if this is an actual problem or just that for large lattices, its unavoidable for the CPUs to need to do more waiting, but in the meantime here's an issue to track in case anyone knows what may be going on.

It's not necessary for local updates and for global ones, it can be implemented and called on the model site if necessary.