Custom kets for KPM solvers

[pablosanjose]

With the current design, one can only run the KPM GreenSolver on kets that form a local basis of contacts in an OpenHamiltonian

This means, in particular, that we cannot use a random ket to perform stochastic traces, where the ket has a random amplitude on part or all sites. It also means that we lose the possibility of efficiently computing other things, such as $\langle \psi | G|\psi\rangle$ where $\psi$ is a linear combination of many sites.

We therefore need to design an API to encode a custom ket for use in KPM solvers.

[pablosanjose]

Option 1

Define a new type of SelfEnergy of type Ket that can be attached, e.g. like

oh |> attach(ket(r-> rand()*cis(2π*rand())); sites...)

Then, the KPM apply will construct that ket when it sees a SelfEnergy of this weird type. Pros: little changes to the current, simply extend it with "one more type". Cons: a ket as a SelfEnergy makes no sense whatsoever, conceptually

Option 2

Put the ket into the KPM solver, and if a ket was specified, access it with a special indexing. Something like

g = h |> greenfunction(GS.KPM(order = 1000, ket = ket(r-> rand()*cis(2π*rand()); sites...))

The problem is how to index it. We currently index greenfunctions etc with either an integer (contact index) or a selector. To address a ket this falls short. Unless we equate a ket to a contact somehow, like in Option 1. So Option 2 would imply extending indexing to allow for other options, e.g g(ω)[ket(1)]. This looks ... weird.

EDIT: although, actually, g(ω)[ket(1), ket(2)] as a way to indicate $\langle 1|G(\omega)|2\rangle$ might actually be a good general concept
`

[pablosanjose]

Option 3

Create a Projected{<:GreenFunction} wrapper type that represents p = psi' * g * psi

Here psi is an OrbitalSliceMatrix with the first axis over g's lattice, and the second axis missing. This p::Projected{<:GreenFunction} could have purely integer indexing over the missing axis of psi. By default (fallback), the indexing p(ω)[i,j] becomes (psi' * g(ω)[slice, slice] * psi)[i,j]. For p only Integer indices are accepted, that represent psi columns, not contacts.

One can then also add new custom indexing methods for the g API that accept a portion of psi instead of a CellOrbs to implement specialized methods. Or perhaps replace the basic indexing API to rely on an OrbitalSliceMatrix?

With this design, building p would trigger a recomputation of KPM momenta over psi. Where would they be stored? In the parent g?