TODO

Features:
-make truncate, and orthogonalize that work on copies.
-Support Ncell>1 in InfinitSum{MPO}(InfiniteMPO)  constructor.
-Use an infinite=true flag for iMPOs?
-Add the two site BH model to models lists.

Improvements:
-assign!() is still inefficient.  We need a way to do
  W[r=>1,c=>1]=op where order(op)==2 and order(W)==4
  at the NDTensors level.
-better name for pbc? inf open_chain, boundary_conditions
-iMPO truncate fails with a trivial 2x2 hamiltonian.
Done: -refactor extract_blocks for work with reg_form_Ops
-Gauge fix. Solve s1=A2*A3*A4..AN*A1*s1-A2*A3*A4..AN*b1 - A2*A3*A4..AN-1*b4 .....

-Not Easy?: WARNING: using SparseArrays.permute in module ITensorMPOCompression conflicts with an existing identifier.

-Subtensor
  -Clean up "where {Elt,N}" code.
  -check flux of inserted blocks.
  -support W[il,ir]=op ?
  -Check that spaces are compatable
  -getindex[] returns a view
  -Clean up emptyITensor if statement.  Get poly dispatch working


ITensorInfiniteMPS:
  Done: -Need code for <ψ|iMPO|ψ> sowe can check <ψ|...LWLinv LWLinv LWLinv...|ψ>=<ψ|...WWWWW...|ψ>
  Done: -VUMPS for iMPO
  -Profile vumps iMPO
  Done: -left/right environment for non-tri MPOs (post compression).  More inversions
  Done: -return el/er from vumps.  
  -iMPO analysis for 1) ops on diag 2) unit on diag
  -Do L/R environments one QN link-block at a time?
    Do we need corner I's for every link block?
    Maybe when D id big this becomes irrelevant.
  -Fix L[1] and R[Dw], they are not self consistent right now.
    Just need to subtract e(I| in the test?
  Done: -Fix if N=1 stuff in tdvp_iteration_sequential tdvp_iteration_parallel
  Done: -Or even better do really need two orth sweeps?  Write a demo that compares get_Dw(H) for 1 and two sweeps
   for H=HeisenbergNNN, HubbardNNN and Loic's fqhe
    Yes we need 2 sweeps
  -Support complex InfiniteMPO which is built from InfiniteMPOMatrix, which in turn is built from InfiniteSum{MPO}
  -Why does heisenbergNNNs [1,2,3,4] compress differently than [4,3,2,1]
  -See energy discrepencies [4,1,1,1] and [1,4,1,1] !?
    Does this happen with InfiniteSum{} ?
    What about [4,1] and [1,4]:  Yes also different 

Old:
Done: -Fix Index Tag conventions sections of the docs.
Done: -Fix emprical observations section of the tech notes.
Done: -Put some content in Readme.md
Done: -Test with Julia 1.7 and 1.8

Resolved: How to handle ortho center for N=1 unit cell iMPOs?
  -InfiniteMPO using llim=b-1,rlim=b+1 for b=orthcenter(H) fails when N=1, because 
  left/right/unorth all become indistinguishable.
Done: Rip out direct code for upper?  
  Do Hu->HL by either:
    1) Transpose each W. Might move data in ram ?
    2) Reverse order in lattice -l0-W_1-l1-W_2-l2- ... lN-1-WN-lN -> lN-WN-lN-1-WN-1-......
        This transposes each W without moving any data. 
No: -Do we need to support multiple truncate sweeps for iMPOs?
Done: using reg_form_MPO wrapper internally.
   -Save regular_form upper/lower/irreg/unknown in the MPO struct so we never have to call
    detect_regular_form.
Done: -One last cleanup: impo truncate
Done: -Subtensor cleanup.
Fixed:   iMPO After left orth sweep, 3 iterations Dw reduced from 1 to 4
    -2-body Hamiltonian in trancate fail.
Fixed: -Truncate iMPO transIsing_iMPO NNN=4 has a gauge problem.
  Maybe just need 2 more orth sweeps (Auto iMPO is fine)
Done: -Organize exports and imports.
Done: -Consider default rr_cutoff=1e-14?
Done: -Refactor out all the link fixing and PBC code for MPO making functions.
Done: -Search for TODO items and get them on this list.
Done: -System for hamitonians Hubbard_AutoiMPO => Hubbard_Auto_iMPO
Done: -Clean up "where {Elt,N}" code.
Done: -Eliminate use of tags for finding indices.  Use the network structure instead.  Pass list if indices
Done: -Move Hamiltonians into testing
Done: -Change flip to mirror?
Done: -Fix max(spectrum)  Base.max
Done: -Fix ITensorMPOCompression.transpose
Done: -less use of parse_links.  
  Still need it in translatecell.
Done: -Follow Matts kwargs guidelines
Done: -Take out support for upper.
  Done: Convert upper H's to lower durung the conversion to ref_form_MPO
  Strip out lots of code.
  After InfiniteCanonicalMPO implementation
Done: -How to handle eps for is_gauge_fixed functions.
Done: -Split up reg_form_Op, reg_form_MPO and reg_form_iMPO
No:  -Make a separate ITensorsTestHamiltonians module
Done:: -qx.jl  one step:
      Ŵ = ITensor(0.0, ileft, iright, siteinds(Ŵrf))
      Ŵrf⎖ = reg_form_Op(Ŵ, ileft, iright, Ŵrf.ul)
    to
      Ŵrf⎖ = reg_form_Op(ileft, iright, siteinds(Ŵrf)) #default to lower.
Done: -Start using 
    struct InfiniteCanonicalMPO <: AbstractInfiniteMPS
      AL::InfiniteMPO
      C::InfiniteMPO
      AR::InfiniteMPO
    end
  -How to coordintae with reg_form_MPOs?
  -Build reg_form Hamiltonians by default.
  -Follow Matts kwargs guidelines
Done: InfiniteCanonicalMPO Plan:
  -Write a usage case for orthgonalization:
    Hi = InfiniteMPO(model, s; model_kwargs...)
    Ho = orthogonalize!(Hi) returns InfiniteCanonicalMPO
    @test check_ortho(Ho)

    similar for truncation
NO: -Go over all the Problems in VScode.  Are they real?
  Not seeing any now.  Did formatting fix this?
Done: -format

NO: -For iMPOs set llim,rlim=-1,1 so that they are interpreted as non-ortho when N=1.
  Not going to use these
No: -iMPO truncate user HL,HR instead of H,Hm ?

Done: -Remove the bail out block in truncate(MPO)?
Done: -More usage of replacetags(reg_form_OP,...)
Done: -Get rid of the old redim functions
Done: -test_three_body should do better on Dw for inf lattice.
Done: -Get rid of getV1 and setV1.  i.e. make all order(W)==4 by adding dummy indices.
Done: -Move code into compression module.
Done: -DO if inside gauge_fix 
  if !is_gauge_fixed(H,eps)
        gauge_fix!(H)
    end
Done: -Overload ref_form_Op*ITensor can we auto fix ileft,iright?
  Yes: -Will this help clean up truncate(iMPO) code?
Done: -no method matching order(::ITensorMPOCompression.reg_form_Op)
Done: -Fix links in reg_form_iMPO
Done: -make one is_gauge_fixed function for reg_form_MPO and reg_form_iMPO 
Done: -Check in reg_form_Op translator.
Done: -Do ac-orthogonalize! eta = norm(R-I),  do this block by block for block sparse R
Done: -get rid of matrix_state
Done: -check and fix equal_edge_blocks in truncate.  Use redim1 for all orth/truncate code.
Done: -Replace getV with another member in extract_blocks
Done: Make forward(Wrf::ref_form_Op), and backward(Wrf::ref_form_Op)
Done: -Detect_upper_lower should use extract blocks.
Done: -Simplify the whole getM process in MPO truncate.
Done: -Get rid if triangular tests.
Done: -ac_qx option to prime iqx ?
Done: -Use proper 𝐛̂ hatted symbols for OP valued tensors.
Done: -parse_links ... accept n=1 for link indices.
 Done: -Take out the h_mirror feature.
Done: -Test with pivoting
Ac blocking code:
  Done 2) -add_dummy_links! return irs,ils for left and right links.  Less dag calls.  harmonize ortho/gauge_fix!
    loops.  No more N+1 indexes.  
  Done: 1) -make remove_dummy_links!, remove_dummy_links, has_dummy_links.  name: edge_links?
  Done: 3) -use vector(t) .. make a short function to get the Dw=1 index first.
  -extract_blocks onlt swap b,c blocks if user requests.  Handle ilb,ilf based on ms
  Done: 4) -ac_qx check that dh is near interger.
  Done: 7) -set_𝐀̂𝐜̂_block  use space(ilf) to get first and last QN blocks.  Use a simplified redim function
  Done: 8)   grow(il,R,ir)  should do the same thing.
  Done: 5) -set_𝐀̂𝐜̂_block  lower-left==upper-right for preserving corners and one row/col
  Done: 6)  same in gauge_fix!(W...)
  Done: 9) -gauge_fix! use a set_c function
  Done: 10) -Fix all noprime(A,tags="Link") calls.
  Done: 11) -Refactor set_𝐀̂𝐜̂_block into
        Done: a) newQ() which makes Wp, preserves columns and Is
        Done:b) call set_𝐀̂𝐜̂_block or set_𝐀̂b_block
        Done:c) set_𝐀̂𝐜̂_block calls set_𝐜̂_block is nr==1
        Done:d) set_𝐀̂b_block calls set_b_block is nr==1
  Done: -Make a struct MPO_with_links,  has ils, irs, ul, d0,dN
  Done:   Same Op_with_links il,ir,ul
    -Ask Matt if he is already doing this?
  Done:   -Do QR pivot PR first.
Fixed: -Incomplete block transfer: truncate (transIsing_MPO,upper,"S=1/2"),
Done: -Replace RL with X in bock_qx.
Done: -refactor gauge fix in Block QX
Done: -Clean up getM()

Done: -MPO truncate Heisenberg AutoMPO
Done: -MPO truncate ... don't switch to dense for svd(M) ... need to preserve QN structure
Done: -Get Heisenberg AutoMPO iMPOs working (Need to handle case where NNN>N)
Done: -Upper reg-form from autompo
Done: -iMPO from AutoMPO
Done: -Test alternating S=1/2 and S=1 sites for MPO.
Done: -Test alternating S=1/2 and S=1 sites for iMPO.
Done: -Test Heisenberg model from autoMPO with non-trivial QN space structure on Link indices.
Done: -Test Fermions/Hubberd.

Done: -If user doesn;t set rr_cutoff, set to 1e-15
Done: -Enforce expected type for all get(kwargs) calls
Done: - use maxlinkdim for maxlinkdim
Done: -Review high level algos
  MPO truncate should iterate until no bailouts.
-Done: convergence in truncate for finite MPO?
  -Keep sweeping until there are no bailouts.
 
No: -rr_cutoff=0.0 not allowed in truncate?
Not a bug: -Handle 3 site unit cell for iMPO with QNs.
  This crashed when creating sites inside ITensorInfiniteMPS with error:
  Dividing ITensors.QNVal(SzParity, 1, 2) by 3, the resulting QN value is not an integer
Done -iMPO, if Ncell=1 the left/right tags of H[1] are the same.  We need to prime one of them
  prior to truncation, maybe during or prior to orthogonalization.
- Check for unit ops along the diagonal: For rectangluar matrices we seem to be looking a the  wrong diagonal.
Done: -Can we use linkind(MPO,bond::Integer) to get link indices? Less use or parsle_links?
Done: -User llim and rlim in MPO/iMPO structs to save ortho status.
  -unortho llim=0, rlim=N+1
  -Ortho on b llim=b-1, rlim=b+1
  -Useful functions in abtractmps.jl: leftlim rightlim setleftlim! setrightlim! ortho_lims
    set_ortho_lims! reset_ortho_lims! isortho orthocenter 
    linkind linkinds defaultlinktags hasdefaultlinktags eachlinkinds eachsiteinds
  -Encountered ambiguity problems with iMPO when N=1.
Done: -Clean up index matching in SetV
Done: -Pare down unit test for shorter run time.
  -Very hard, lots of CPU time on compiling inside ITensors for GS energy tests
  -Using a sysimage for ITensors does not seem to help!
Done: -Profile and code coverage.
   Spends most of the time in detect_regular_form.  We should keep this info in the MPO
   struct.
Done: -Gs/Ss returned by truncate should be Diag not Dense tensors.
Done: -Use new is_upper functions using combiners to get rank-2 tensors.
   Not sure the combiners do what we want here for an op-valued matrix.
   Using norm and no slice to test for zeros
Done: -Add some sort of Results section at the end ot tech docs.
Done:  -rules for smart_sweep
Done:  -orthog inside truncate
Done: -Switch from epsrr to rr_cutoff
Done: -Clean up SetV code.  Use dummy Dw=1 index trick instead of SetV1 code.
Done: -clean up kwargs=Dict(kwargs) is not needed.
Done: -qx_iterate! clean up GS init loop

Done: -Verbose mode show all Dws and steps orth-right etc.
  Done: -subtensor function
  	-Must handle block sparse
  	-I only need it for link indices
  	-Do the work in NDTensors

Done: -Does block_qx really need to know n?
  Instead of r,c index we should really think in terms of f,r = forward, reverse
  R gets transferred in f/forward direction
Done: -Implement @show for bond_spectrum 
Done: -Clean up the whole parse links mess to handle iMPOs
Done: -Support transIsing_iMPO
Done: -For the make hamiltonian functions put ul,hx,J into kwargs
  This will tie in with using the Model("model_name") framework Matt is using.
Done: -Review all examples.  Add iMPO example.
Can't do (for now: -Reproduce 30 site example from Parker paper.
  I estimate Dw>>1000 for N=30.  They must be building H_hand more efficiently.
  Anyway for N=9 we get Dw=206 which truncates nicely to Dw=9, takes about a minute just to build up H though!
Done:  Need to build H-three_body from scratch.
Done:  Need code for W_1 + W_2
Done: -Try using cutoff instead of truncate=false for autoMPO
Done: -Clean up iMPO truncate code.
Done: -Move i_orthog code into orthogonalize.jl files
Done: -Review generated docs, upgrade for iMPO support.
Done: -QN support for InfiniteSum code.
Done: -Use Itensor spectrum instead of bond_spectrum.
Done: -Use sweep(H,lr) instead of the if blocks.
Done: -Refactor the H[n]*U dag(U)*H[n+1] code for iMPO truncate.

Done: 1) Not seeing any rank reduction
Done: 2) Upgrade tests to use ~ and atol.
Done: 3) Upgrade tests to use for loops.
Done: 4) ITensorMPOCompression/test/hamiltonians.jl
Done: Deal with: Calling `inner(x::MPS, A::MPO, y::MPS)` where the site indices of the `MPS`
│ `x` and the `MPS` resulting from contracting `MPO` `A` with `MPS` `y` don't
│ match is deprecated as of ITensors v0.3 and will result in an error in ITensors
│ v0.4. The most common cause of this is something like the following:
Done: 6) orthogonalize! -> ITensors.orthogonalize! for proper overloading
Done: 7) truncate! -> ITensors.truncate! for proper overloading
Done: 8) Use iq returned from qr()
Yes:  9) Do we still need disable_debug_checks before calling dmrg() ?


Done: -make block_qx return link name qx instead of ql,lq,qr,rq
Done: -calculate getV,setV offsets from matrix_state
Done: -We need pure upper/lower enums without diagonal and full options
Done: -further simplify the compress H loop
Doen: -Need to save SV spectrum somewhere
Done: -Look into the MPS compession code for inspiration.
Done:    rename compress! to trunctate! (supprts the svd args)
        cutoff  Using a cutoff allows the SVD algorithm to truncate as many states as possible while still ensuring a certain accuracy.

        maxdim –- integer M. If the number of singular values exceeds M, only the largest M will be retained.

        mindim –- integer m. At least m singular values will be retained, even if some fall below the cutoff
    rename canonical! to orthogonalize! and orthogonalize (support left/right rankreveal)

Done: -Figure out how to turn off asserts for performance.
-consistent variable names for link indexs: ir,ic,iq,iu,iv ?
Fixed: #grow(i1,W,i2) relies on index ordering inside W ... is that bad?? (hint:yes!)
Done: -grow(i1,W,i2) --> grow(W,i1,i2)
Done: -Fix known fail for compressing when H has 8 NNN  Do tidy up items first.
Done -Get RQ decomp going. in order to ....
Done -Support upper-tri MPOs
Done: -Support rank releavling QR,RQ,QL,LQ
Done: -Clean up canonical and compress code.
Done: -reduce two block_qx functions to one
Done: -Do we need to copy W in block_qx? Yes
-Comment indices in orthogonalize!(...)
Done: -use index returned by ql,qr,lq,rq
Done -We need rank reaving QR/QL/RQ/LQ throughout, otherwise we get zero col/rows away from the edges which causes no end of trouble.
-is_regular is looking along the wrong diagonal for upper rectangular.  Not sure any more
Fixed: -for lower,right SetV seems to require the last column(above I)  should be preserved ... 
  but why?       It should be zeros.  Index ordering?
Done: -Make rank reaviling an option for block_qx
Done: -For obc MPO try compressing 2..N-1, what do we end up with for W[N].  Then try N..2, what does W[1] look like??  Can we get around PBC issue this way?
Done: -Test canonical with rr_cutoff=0.0 and compare.  need an maxlinkdim(H::MPO) function.  
-Use  # Exact energy for transverse field Ising model J=-1, hx=-1
    # with open boundary conditions at criticality
    energy_exact = 1.0 - 1.0 / sin(π / (4 * N + 2))
    @test abs((energy - energy_exact) / energy_exact) < 1e-4
    Can we get fast convergence for small lattices?
Done: -add Printf for 
Done: -simplify MPOpbc code
Done: -Get rid of MPOpbs struct.
Done: -Prevent Documter.jl from showing ITensor docs.
Done: -Seeing flakey GS tests.  probably need to seed the randomMPS for consistency.
Done: -Add a big lattice test and look at SVs
Plan:
  90% Document what I have so far
  Done: Test with fixed autoMPO Hamiltonians (start the heisenberg with arbitrary neighbours)
  Done: Test with DMRG GS calculations (should just work, but who knows right?)
  Test with Fermion/Hubbard  Hamiltonians
  Test with the 3 body hamiltonians used in the Parker paper (If I can figure what they were doing)
  Test with QN's (new area for me)
  50% Refactor the code.  A lot of things can be simplified once there are lots of unit tests in place.
  Done: Figure how to disable 100's of asserts for production mode usage.  
  Done: Incoroprate Ham functions into module
  Done: Organize exports
  90% Add examples to docs (needs Ham functions)
  

Matts emails:
Nice meeting you and thanks again for offering to help us out with the development of ITensor!

Here are the Julia development tools that were brought up, and some other ones I thought of:

Julia debuggers (the primary ones appear to be the one built into the Julia VSCode extension and Debugger.jl):
https://opensourc.es/blog/basics-debugging/
https://www.julia-vscode.org/docs/stable/userguide/debugging/
https://github.com/JuliaDebug/Debugger.jl

Julia vim plugin:
https://github.com/JuliaEditorSupport/julia-vim

OhMyREPL for nice syntax highlighting in the Julia REPL:
https://github.com/KristofferC/OhMyREPL.jl

Revise is absolutely critical for developing Julia packages, since it allows you to edit code when you are developing Julia packages and have the changes immediately reflected outside the package, without restarting the Julia session:
https://github.com/timholy/Revise.jl

I would also recommend using the JuliaMono font to help with viewing unicode (many fonts have limited unicode support):
https://juliamono.netlify.app/

Here is the profile visualization tool I use:
https://github.com/timholy/ProfileView.jl

For a graph visualization of the profile you can use:
https://github.com/JuliaPerf/PProf.jl
https://opensourc.es/blog/constraint-solver-profiling/

And here is the code I wrote for making an operator valued matrix from an MPO:
https://github.com/ITensor/ITensorInfiniteMPS.jl/blob/main/src/vumps_mpo.jl#L36-L64

We actually have an immediate application of the MPO compression in Julia (there is an intern working on a project where he is renormalizing a local Hamiltonian by a unitary circuit, so it would be great to automatically compress the MPO instead of computing the new Hamiltonian more manually). I also have some quantum chemistry applications where that would be useful. I think that would be a top priority, though of course feel free to work on whatever peaks your interest most.

I don't think you need any particular Hamiltonian. A useful test would be to start with any local Hamiltonian (say the Heisenberg model) and time evolve it by any local local circuit, it could even be totally random gates. I'm attaching an example you could start from.

Unfortunately there is a slightly circular logic in this example, since the lines:

UHU = apply(U, H; apply_dag=true, cutoff=1e-15)

that performs the circuit evolution also does an MPO truncation, however it uses the naive form of truncation that treats the MPO as a state. It would also take the MPO out of the upper/lower triangular form. So I guess we would have to think about how the basis rotation should actually be applied in this case. Ultimately the goal would be for the apply have an option where you can specify the MPO should be treated as a local Hamiltonian and truncated using the Parker et al. algorithm.

Let me know if you think that example makes sense and if there is a way you could make use of it. Perhaps it is best to think about it in terms of just applying one of the gates of the circuit at a time, and truncating just the bond where it is applied. In that way, we could engineer the gate so it is applied at first without truncation and then truncated using the Parker et al. algorithm.

Also, maybe it is helpful to know that you can turn the circuit into an MPO like this:

Id = MPO(OpSum() + ("Id", 1), s)
Umpo = apply(U, Id; cutoff=1e-15)

That form could make it easier to apply the circuit (or just a gate) to the Hamiltonian MPO to keep it in a form where it can be used with the Parker et al. algorithm.


julia --sysimage /home/janr/.julia/sysimages/sys_itensors.so

open("/home/janr/prof1.txt", "w") do s
           Profile.print(IOContext(s, :displaysize => (24, 500)),format=:flat)
       end