More explicitly:

test_3_solver (tests.test_hubbardFermiMatrix.TestHubbardFermiMatrix)
Test Ax=b solvers. ... [2020-08-07 14:25] Testing solveM on C20
libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: Incompatible buffer format: mismatched elemental data type
Abort trap: 6

and

> python3 fullHMCEvolution.py
libc++abi.dylib: terminating with uncaught exception of type pybind11::stop_iteration:
Abort trap: 6

There is some room for defining the chemical potential. It might be better to change it. See #17

It'd be nice to be able to say for s in isle.Species instead of for s in (isle.Species.PARTICLE, isle.Species.HOLE), for example.

Fermion matrices for Species.PARTICLE are actually hole matrices and vice versa in either exp or dia discretization. This is due to a sign error in the hopping parameter kappa.

Need to find out which one is correct (compare to other code) and fix the broken one.

Is this visible on slack?

Just writing this down someplace safe. Discussing with Alexei Bazavov today he mentioned to me the Wang and Landau algorithm (PRL) which is a method for constructing the density of states.

Quoting Wikipedia,

The major consequence is that this sampling distribution leads to a simulation where the energy barriers are invisible. This means that the algorithm visits all the accessible states (favorable and less favorable) much faster than a Metropolis algorithm.

but neither is it just random sampling. I have to say I do not yet understand it---I'd have to code up an example myself to really grok it. But it may be of interest!

There are too many vectors being copied around the place! Most importantly:

• hmc (initial state and integrator)
• solver (copy vectors in for multiple RHSs, maybe more on Python side)

It can be very cheap to calculate the energy inside of a proposer. It might already be calculated by an integrator anyway. So just return it along with field and momentum.

Currently, the convention for basis transformations is inconsistent.

In onePointFunctions it's vector --> U vector but in singleParticleCorrelator and spinSpinCorrelator it's matrix --> U† matrix U. This causes lots of very annoying problem.

We should harmonize, unifying on matrix --> U matrix U†.

If computing trajectories takes long, it is useful to have regular updates in the log file to keep track of progress.

Talking with Alexei Bazavov he pointed me to two possible techniques of interest. The first is parallel tempering (first paper)

The idea is that you start two (or more) separate ensembles with different actions (for example, introducing a staggered mass, or with a different temperature, or a different U...) that I'll call A and B.

You separately evolve ensemble A according to action A and ensemble B according to action B. Then, at some point, you proposing exchanging the field configurations.

You imagine that you start with action S_A(phi_A) + S_B(phi_B) and Metropolis accept-reject the exchange according to the new action S_A(phi_B) + S_B(phi_A). You can do this with multiple ensembles, rather than just 2.

This could be implemented via MPI: the two (or more) streams are trivially parallelizable except for a once-in-a-while swap step. (Alexei points out that it's easier to MPI communicate action parameters than field configurations).

Swapping with a staggered mass ensemble (for example) would formally remove ergodicity problems. Changing T could help overcome tunneling difficulties for first-order phase transitions. However, there can be a bad overlap problem---if the two actions like field configurations that are too different you'll hardly accept.

For large lattices, not all measurement results can be kept in RAM at the same time. Especially correlators run out of memory quickly and crash Isle.

Possible solution: keep only one (or a few) results in memory and write to HDF5 early. This could be done with a class which can be given numpy arrays one after the other and writes them to disk as necessary.

This would require:

• Changing the way measurements are saved. They would need to do it themselves instead of waiting for a call to save. This means the measurements need to know file name / HDF5 group at initialisation.
• Computing the number of configurations to measure on ahead of time to avoid resizing datasets and fragmenting files.

This might allow for some other improvements:

• Handle configurations in the same way i.e. collect several and save them together. Maybe difficult because configs are saved into separate datasets.
• Consolidate measurement output to reduce number of file accesses.

Do you have ideas / suggestions on this?

In most measurements that require propagators, like SpinSpinCorrelator and DeterminantCorrelators one has to pass a particle AND a hole propagator. It would be much simpler if the propagator.AllToAll returned just a dict from species to propagators.

Right now SingleParticleCorrelator takes just one species. But usually, you call it twice (see measure.py), once for particles, once for holes. So we could clean that up similarly, allowing just one call in measure.py` by passing the dict.

The only concern is speed---requiring an inverse for both species means the user always pays for both. We could add another measurement propagator.OneSpecies (or something like it).

This would be a breaking change, but not difficult to implement.

Can we implement a Schur complement solver, as in http://arxiv.org/abs/1803.05478 ? We have encountered some numerical instability when Ubeta becomes too large---presumably the Schur solver does not suffer?

(Unfortunately it is only for exp discretization.)

For example, when µ=0 on bipartite lattices, phi --> –phi? Let's say I want to measure on both, which I know is allowed because of https://arxiv.org/abs/1812.09268 IV.B

The single particle correlator measurement does not work for C20. Running it with Nt=16 gives

Traceback (most recent call last):
  File "docs/examples/measure.py", line 72, in <module>
    main()
  File "docs/examples/measure.py", line 69, in main
    measState(measurements, configRange=args.configs)
  File "/home/jl/Uni/PhD/cns/isle/src/isle/drivers/meas.py", line 29, in __call__
    self.mapOverConfigs(measurements, configRange)
  File "/home/jl/Uni/PhD/cns/isle/src/isle/drivers/meas.py", line 55, in mapOverConfigs
    measurement(phi, action, i)
  File "/home/jl/Uni/PhD/cns/isle/src/isle/meas/singleParticleCorrelator.py", line 42, in __call__
    propagators = res.reshape([len(self.irreps), nt, nt, len(self.irreps)])
ValueError: cannot reshape array of size 81920 into shape (16,16,16,16)

Don't know whether this is an old bug and no one noticed or whether it was introduced in 7bf5ad3.
I couln't reproduce this with any other lattice.

Implement a way to concatenate actions directly instead of collecting them in Hamiltonian.

Reasons:

• nicer to use
• doesn't force adding kinetic term => lift dependency on HMC

Needs:

• a new action for HMC kinetic term
• a new 'meta action' to store individual concrete actions
  (Can this work with the old hybrid pointer?)

For numeric reproducibility it would be good to store the einsum paths that different observables use. (It depends on the lattice size, nt, etc.)

Annoyingly the output of np.einsum_path can't just be dumped to h5.

import h5py as h5
import numpy as np

M = np.arange(1024).reshape(4,4,8,8)

path, summary = np.einsum_path('abde,acdf,bcef->',M,M,M,optimize="optimal")

# path = ['einsum_path', (0, 1), (0, 1)]
# print(summary) yields
#   Complete contraction:  abde,acdf,bcef->
#         Naive scaling:  6
#     Optimized scaling:  6
#      Naive FLOP count:  9.830e+04
#  Optimized FLOP count:  6.758e+04
#   Theoretical speedup:  1.455
#  Largest intermediate:  1.024e+03 elements
#--------------------------------------------------------------------------
#scaling                  current                                remaining
#--------------------------------------------------------------------------
#   6             acdf,abde->bcef                              bcef,bcef->
#   4                 bcef,bcef->                                       ->

# Another example
path, summary = np.einsum_path('abde,acdf,bcef,gghh->',M,M,M,M,optimize="optimal")

but because it's a mixed array path can't easily be written to an h5 path.

I tried pickling the whole shebang, but I couldn't figure out how to send the pickle.dump to an h5 field. There is an hdf5 pickle but I'd be afraid to use it and h5py at the same time.

How should this be accomplished? Sometimes many different contractions share an optimization (as in spin-spin correlators). Could be a string attribute (str(path) works).