Support MPO-(MPS with dangling site indices) contraction

I made a fresh start after #44 was merged into the main branch.

For QTT, we may want to contract an MPO (with two site indices at each vertex) and an MPS with dangling site indices. I found several issues, but let me sout them out one by one. The first issue is the following.

A minimum example for naive algoritm works:

nbit = 3
sites = siteinds("Qubit", nbit)

H = randomMPO(sites) + randomMPO(sites)
H = replaceprime(H, 1=>2, 0=>1)

psi = randomMPO(sites) + randomMPO(sites)

# MPO-MPO
Hpsi_ref = contract(H, psi; alg="naive")

# MPO-MPS
Hpsi_ref2 = contract(H, MPS([A for A in psi]); alg="naive")

@show H
@show psi
@show Hpsi_ref
@show Hpsi_ref2
@assert MPS([A for A in Hpsi_ref]) ≈ Hpsi_ref2

H = MPO
[1] ((dim=2|id=58|"Qubit,Site,n=1")'', (dim=2|id=58|"Qubit,Site,n=1")', (dim=2|id=68|"Link,l=1")')
[2] ((dim=2|id=656|"Qubit,Site,n=2")'', (dim=2|id=656|"Qubit,Site,n=2")', (dim=2|id=121|"Link,l=2")', (dim=2|id=68|"Link,l=1")')
[3] ((dim=2|id=270|"Qubit,Site,n=3")'', (dim=2|id=270|"Qubit,Site,n=3")', (dim=2|id=121|"Link,l=2")')

psi = MPO
[1] ((dim=2|id=58|"Qubit,Site,n=1")', (dim=2|id=58|"Qubit,Site,n=1"), (dim=2|id=787|"Link,l=1"))
[2] ((dim=2|id=656|"Qubit,Site,n=2")', (dim=2|id=656|"Qubit,Site,n=2"), (dim=2|id=46|"Link,l=2"), (dim=2|id=787|"Link,l=1"))
[3] ((dim=2|id=270|"Qubit,Site,n=3")', (dim=2|id=270|"Qubit,Site,n=3"), (dim=2|id=46|"Link,l=2"))

Hpsi_ref = MPO
[1] ((dim=2|id=58|"Qubit,Site,n=1")'', (dim=2|id=58|"Qubit,Site,n=1"), (dim=4|id=716|"CMB,Link"))
[2] ((dim=2|id=656|"Qubit,Site,n=2")'', (dim=2|id=656|"Qubit,Site,n=2"), (dim=4|id=206|"CMB,Link"), (dim=4|id=716|"CMB,Link"))
[3] ((dim=2|id=270|"Qubit,Site,n=3")'', (dim=2|id=270|"Qubit,Site,n=3"), (dim=4|id=206|"CMB,Link"))

Hpsi_ref2 = MPS
[1] ((dim=2|id=58|"Qubit,Site,n=1")'', (dim=2|id=58|"Qubit,Site,n=1"), (dim=4|id=596|"CMB,Link"))
[2] ((dim=2|id=656|"Qubit,Site,n=2")'', (dim=2|id=656|"Qubit,Site,n=2"), (dim=4|id=816|"CMB,Link"), (dim=4|id=596|"CMB,Link"))
[3] ((dim=2|id=270|"Qubit,Site,n=3")'', (dim=2|id=270|"Qubit,Site,n=3"), (dim=4|id=816|"CMB,Link"))

fitting algorithm does not work:

# MPO-MPS
contract(H, MPS([A for A in psi]); alg="fit")

BoundsError: attempt to access 1-element Vector{Index{Int64}} at index [2]

Stacktrace:
  [1] getindex
    @ ./array.jl:924 [inlined]
  [2] replaceinds(is::Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, inds1::Vector{Index{Int64}}, inds2::Vector{Index{Int64}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/indexset.jl:570
  [3] replaceinds(::NDTensors.DenseTensor{Float64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, ::Vector{Index{Int64}}, ::Vararg{Vector{Index{Int64}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/itensor.jl:1480
  [4] replaceinds(::NDTensors.DenseTensor{Float64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, ::Vector{Index{Int64}}, ::Vector{Index{Int64}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/itensor.jl:1479
  [5] replaceinds!(::ITensor, ::Vector{Index{Int64}}, ::Vararg{Vector{Index{Int64}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/itensor.jl:1484
  [6] replaceinds!(::ITensor, ::Vector{Index{Int64}}, ::Vector{Index{Int64}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/itensor.jl:1483
  [7] contract(::ITensors.Algorithm{:fit}, A::MPO, psi0::MPS; init_state::MPS, nsweeps::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensorNetworks /sharehome/shinaoka/git/ITensorNetworks.jl/src/treetensornetworks/solvers/contract_operator_state.jl:39
  [8] contract
    @ /sharehome/shinaoka/git/ITensorNetworks.jl/src/treetensornetworks/solvers/contract_operator_state.jl:13 [inlined]
  [9] contract(A::MPO, ψ::MPS; alg::String, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors /sharehome/shinaoka/.julia/packages/ITensors/fpBnt/src/mps/mpo.jl:567
 [10] top-level scope
    @ In[14]:2

The replacement of site indices does not work with an MPS with dangling site indices. Shall we just remove this part?
https://github.com/mtfishman/ITensorNetworks.jl/blob/main/src/treetensornetworks/solvers/contract_operator_state.jl#L38-L42

If we move those lines, it will be user's responsibility to provide an appropriate init_state with the same site indices as H * psi0. The following default value assumes that init_state lives in the same space as psi0, which is the case if H is a time-evoluation operator. For a more general MPO, this default value is not appropriate. What do you think about removing these lines for the replacement and adding a check on the site indces on init_state?

https://github.com/mtfishman/ITensorNetworks.jl/blob/main/src/treetensornetworks/solvers/contract_operator_state.jl#L17

Start in src/treetensornetworks/solvers/contract_operator_state.jl, especially the site-index replacement around lines 38–42 and the default init_state handling near line 17. Reproduce the provided MPO-MPS example with alg="fit" and use the BoundsError as the initial failure; done means the fitting contraction handles dangling site indices or clearly validates init_state compatibility.