From 76788cc7c6f7342398911c7a9e0d7129bc6e1595 Mon Sep 17 00:00:00 2001
From: Joachim Brand <joachim.brand@gmail.com>
Date: Fri, 12 Jan 2024 22:32:37 +1300
Subject: [PATCH] Add support for initiator rule in default starting vector

---
 src/lomc.jl  | 31 +++++++++++++++++++++----------
 test/lomc.jl |  3 +++
 2 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/src/lomc.jl b/src/lomc.jl
index c7e0e6197..d27d8fa97 100644
--- a/src/lomc.jl
+++ b/src/lomc.jl
@@ -284,19 +284,24 @@ end
     default_starting_vector(
         address=starting_address(hamiltonian);
         style=IsStochasticInteger(),
+        initiator=NonInitiator(),
         threading=nothing
     )
 Return a default starting vector for [`lomc!`](@ref). The default choice for the starting
 vector is
 ```julia
-v = PDVec(address => 10; style)
+v = PDVec(address => 10; style, initiator)
 ```
-if threading is available or
-
+if threading is available, or otherwise
 ```julia
 v = DVec(address => 10; style)
 ```
-otherwise. See [`PDVec`](@ref), [`DVec`](@ref) and [`StochasticStyle`](@ref).
+if `initiator == NonInitiator()`, and
+```julia
+v = InitiatorDVec(address => 10; style, initiator)
+```
+if not. See [`PDVec`](@ref), [`DVec`](@ref), [`InitiatorDVec`](@ref),
+[`StochasticStyle`](@ref), and [`InitiatorRule`].
 """
 function default_starting_vector(
     hamiltonian::AbstractHamiltonian;
@@ -307,14 +312,17 @@ end
 function default_starting_vector(address;
     style=IsStochasticInteger(),
     threading=nothing,
+    initiator=NonInitiator(),
 )
     if isnothing(threading)
         threading = Threads.nthreads() > 1
     end
     if threading
-        v = PDVec(address => 10; style)
-    else
+        v = PDVec(address => 10; style, initiator)
+    elseif initiator isa NonInitiator
         v = DVec(address => 10; style)
+    else
+        v = InitiatorDVec(address => 10; style, initiator)
     end
     return v
 end
@@ -336,6 +344,8 @@ Alternatively, a `QMCState` can be passed in to continue a previous simulation.
 * `address = starting_address(ham)` - set starting address for default `v` and `shift`.
 * `style = IsStochasticInteger()` - set [`StochasticStyle`](@ref) for default `v`; unused
   if `v` is specified.
+* `initiator = NonInitiator()` - set [`InitiatorRule`](@ref) for default `v`; unused if `v`
+  is specified.
 * `threading` - default is to use multithreading and
   [MPI](https://juliaparallel.org/MPI.jl/latest/) if multiple threads are available. Set to
   `true` to force [`PDVec`](@ref) for the starting vector, `false` for serial computation;
@@ -403,9 +413,9 @@ 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
-`v = default_starting_vector(; address, style, threading)`.
-See [`default_starting_vector`](@ref), [`PDVec`](@ref), [`DVec`](@ref), and
-[`StochasticStyle`](@ref).
+`v = default_starting_vector(; address, style, threading, initiator)`.
+See [`default_starting_vector`](@ref), [`PDVec`](@ref), [`DVec`](@ref),
+[`StochasticStyle`](@ref), and [`InitiatorRule`](@ref).
 """
 function lomc!(ham, v; df=DataFrame(), name="lomc!", metadata=nothing, kwargs...)
     state = QMCState(ham, v; kwargs...)
@@ -416,9 +426,10 @@ function lomc!(
     style=IsStochasticInteger(),
     threading=nothing,
     address=starting_address(ham),
+    initiator=NonInitiator(),
     kwargs...
 )
-    v = default_starting_vector(address; style, threading)
+    v = default_starting_vector(address; style, threading, initiator)
     return lomc!(ham, v; address, kwargs...) # pass address for setting the default shift
 end
 # For continuation, you can pass a QMCState and a DataFrame
diff --git a/test/lomc.jl b/test/lomc.jl
index 6e8b5735b..48eae7eb7 100644
--- a/test/lomc.jl
+++ b/test/lomc.jl
@@ -220,6 +220,9 @@ using DataFrames
 
         df, state = lomc!(H; laststep=100, style = IsDeterministic())
         @test StochasticStyle(state.replicas[1].v) isa IsDeterministic
+
+        df, state = lomc!(H; laststep=1, threading=false, initiator=Initiator())
+        @test state.replicas[1].v isa InitiatorDVec
     end
 
     @testset "ShiftStrategy" begin