|
| 1 | +using ITensors.ITensorMPS: ITensorMPS, AbstractProjMPO, MPO, MPS |
| 2 | +using ITensors.ITensorMPS: linkinds, replaceinds |
| 3 | +using ITensors: ITensors, OneITensor |
| 4 | +import ITensorTDVP: alternating_update, rproj, lproj |
| 5 | + |
| 6 | +""" |
| 7 | +A ReducedFitProblem represents the projection |
| 8 | +of an MPS `input_state` onto the basis of a different MPS `state`. |
| 9 | +`state` may be an approximation of `input_state`. |
| 10 | +``` |
| 11 | + *--*--*- -*--*--*--*--*--* <state| |
| 12 | + | | | | | | | | | | | |
| 13 | + o--o--o- -o--o--o--o--o--o |input_state> |
| 14 | +``` |
| 15 | +""" |
| 16 | +mutable struct ReducedFitProblem <: AbstractProjMPO |
| 17 | + lpos::Int |
| 18 | + rpos::Int |
| 19 | + nsite::Int |
| 20 | + input_state::MPS |
| 21 | + environments::Vector{ITensor} |
| 22 | +end |
| 23 | + |
| 24 | +function ReducedFitProblem(input_state::MPS) |
| 25 | + lpos = 0 |
| 26 | + rpos = length(input_state) + 1 |
| 27 | + nsite = 2 |
| 28 | + environments = Vector{ITensor}(undef, length(input_state)) |
| 29 | + return ReducedFitProblem(lpos, rpos, nsite, input_state, environments) |
| 30 | +end |
| 31 | + |
| 32 | +function lproj(P::ReducedFitProblem)::Union{ITensor,OneITensor} |
| 33 | + (P.lpos <= 0) && return OneITensor() |
| 34 | + return P.environments[P.lpos] |
| 35 | +end |
| 36 | + |
| 37 | +function rproj(P::ReducedFitProblem)::Union{ITensor,OneITensor} |
| 38 | + (P.rpos >= length(P) + 1) && return OneITensor() |
| 39 | + return P.environments[P.rpos] |
| 40 | +end |
| 41 | + |
| 42 | + |
| 43 | +function Base.copy(reduced_operator::ReducedFitProblem) |
| 44 | + return ReducedFitProblem( |
| 45 | + reduced_operator.lpos, |
| 46 | + reduced_operator.rpos, |
| 47 | + reduced_operator.nsite, |
| 48 | + copy(reduced_operator.input_state), |
| 49 | + copy(reduced_operator.environments), |
| 50 | + ) |
| 51 | +end |
| 52 | + |
| 53 | +Base.length(reduced_operator::ReducedFitProblem) = length(reduced_operator.input_state) |
| 54 | + |
| 55 | +function ITensorMPS.set_nsite!(reduced_operator::ReducedFitProblem, nsite) |
| 56 | + reduced_operator.nsite = nsite |
| 57 | + return reduced_operator |
| 58 | +end |
| 59 | + |
| 60 | +function ITensorMPS.makeL!(reduced_operator::ReducedFitProblem, state::MPS, k::Int) |
| 61 | + # Save the last `L` that is made to help with caching |
| 62 | + # for DiskProjMPO |
| 63 | + ll = reduced_operator.lpos |
| 64 | + if ll ≥ k |
| 65 | + # Special case when nothing has to be done. |
| 66 | + # Still need to change the position if lproj is |
| 67 | + # being moved backward. |
| 68 | + reduced_operator.lpos = k |
| 69 | + return nothing |
| 70 | + end |
| 71 | + # Make sure ll is at least 0 for the generic logic below |
| 72 | + ll = max(ll, 0) |
| 73 | + L = lproj(reduced_operator) |
| 74 | + while ll < k |
| 75 | + L = L * reduced_operator.input_state[ll+1] * dag(state[ll+1]) |
| 76 | + reduced_operator.environments[ll+1] = L |
| 77 | + ll += 1 |
| 78 | + end |
| 79 | + # Needed when moving lproj backward. |
| 80 | + reduced_operator.lpos = k |
| 81 | + return reduced_operator |
| 82 | +end |
| 83 | + |
| 84 | +function ITensorMPS.makeR!(reduced_operator::ReducedFitProblem, state::MPS, k::Int) |
| 85 | + # Save the last `R` that is made to help with caching |
| 86 | + # for DiskProjMPO |
| 87 | + rl = reduced_operator.rpos |
| 88 | + if rl ≤ k |
| 89 | + # Special case when nothing has to be done. |
| 90 | + # Still need to change the position if rproj is |
| 91 | + # being moved backward. |
| 92 | + reduced_operator.rpos = k |
| 93 | + return nothing |
| 94 | + end |
| 95 | + N = length(state) |
| 96 | + # Make sure rl is no bigger than `N + 1` for the generic logic below |
| 97 | + rl = min(rl, N + 1) |
| 98 | + R = rproj(reduced_operator) |
| 99 | + while rl > k |
| 100 | + R = R * reduced_operator.input_state[rl-1] * dag(state[rl-1]) |
| 101 | + reduced_operator.environments[rl-1] = R |
| 102 | + rl -= 1 |
| 103 | + end |
| 104 | + reduced_operator.rpos = k |
| 105 | + return reduced_operator |
| 106 | +end |
| 107 | + |
| 108 | + |
| 109 | +struct ReducedFitMPSsProblem <: AbstractProjMPO |
| 110 | + problems::Vector{ReducedFitProblem} |
| 111 | + coeffs::Vector{<:Number} |
| 112 | +end |
| 113 | + |
| 114 | +function ReducedFitMPSsProblem( |
| 115 | + input_states::AbstractVector{MPS}, |
| 116 | + coeffs::AbstractVector{<:Number}, |
| 117 | +) |
| 118 | + ReducedFitMPSsProblem(ReducedFitProblem.(input_states), coeffs) |
| 119 | +end |
| 120 | + |
| 121 | +function Base.copy(reduced_operator::ReducedFitMPSsProblem) |
| 122 | + return ReducedFitMPSsProblem(reduced_operator.problems, reduced_operator.coeffs) |
| 123 | +end |
| 124 | + |
| 125 | +function Base.getproperty(reduced_operator::ReducedFitMPSsProblem, sym::Symbol) |
| 126 | + if sym === :nsite |
| 127 | + return getfield(reduced_operator, :problems)[1].nsite |
| 128 | + end |
| 129 | + return getfield(reduced_operator, sym) |
| 130 | +end |
| 131 | + |
| 132 | + |
| 133 | +Base.length(reduced_operator::ReducedFitMPSsProblem) = length(reduced_operator.problems[1]) |
| 134 | + |
| 135 | +function ITensorMPS.set_nsite!(reduced_operator::ReducedFitMPSsProblem, nsite) |
| 136 | + for p in reduced_operator.problems |
| 137 | + ITensorMPS.set_nsite!(p, nsite) |
| 138 | + end |
| 139 | + return reduced_operator |
| 140 | +end |
| 141 | + |
| 142 | +function ITensorMPS.makeL!(reduced_operator::ReducedFitMPSsProblem, state::MPS, k::Int) |
| 143 | + for p in reduced_operator.problems |
| 144 | + ITensorMPS.makeL!(p, state, k) |
| 145 | + end |
| 146 | + return reduced_operator |
| 147 | +end |
| 148 | + |
| 149 | + |
| 150 | +function ITensorMPS.makeR!(reduced_operator::ReducedFitMPSsProblem, state::MPS, k::Int) |
| 151 | + for p in reduced_operator.problems |
| 152 | + ITensorMPS.makeR!(p, state, k) |
| 153 | + end |
| 154 | + return reduced_operator |
| 155 | +end |
| 156 | + |
| 157 | + |
| 158 | + |
| 159 | +function _contract(P::ReducedFitProblem, v::ITensor)::ITensor |
| 160 | + itensor_map = Union{ITensor,OneITensor}[lproj(P)] |
| 161 | + push!(itensor_map, rproj(P)) |
| 162 | + |
| 163 | + # Reverse the contraction order of the map if |
| 164 | + # the first tensor is a scalar (for example we |
| 165 | + # are at the left edge of the system) |
| 166 | + if dim(first(itensor_map)) == 1 |
| 167 | + reverse!(itensor_map) |
| 168 | + end |
| 169 | + |
| 170 | + # Apply the map |
| 171 | + Hv = v |
| 172 | + for it in itensor_map |
| 173 | + Hv *= it |
| 174 | + end |
| 175 | + return Hv |
| 176 | +end |
| 177 | + |
| 178 | +function contract_operator_state_updater(operator::ReducedFitProblem, init; internal_kwargs) |
| 179 | + state = ITensor(true) |
| 180 | + for j = (operator.lpos+1):(operator.rpos-1) |
| 181 | + state *= operator.input_state[j] |
| 182 | + end |
| 183 | + state = _contract(operator, state) |
| 184 | + return state, (;) |
| 185 | +end |
| 186 | + |
| 187 | +function contract_operator_state_updater( |
| 188 | + operator::ReducedFitMPSsProblem, |
| 189 | + init; |
| 190 | + internal_kwargs, |
| 191 | +) |
| 192 | + states = ITensor[] |
| 193 | + for (p, coeff) in zip(operator.problems, operator.coeffs) |
| 194 | + res = contract_operator_state_updater(p, init; internal_kwargs) |
| 195 | + push!(states, coeff * res[1]) |
| 196 | + end |
| 197 | + return sum(states), (;) |
| 198 | +end |
| 199 | + |
| 200 | + |
| 201 | +function contract_fit(input_state::MPS, init::MPS; coeff::Number = 1, kwargs...) |
| 202 | + links = ITensors.sim.(linkinds(init)) |
| 203 | + init = replaceinds(linkinds, init, links) |
| 204 | + reduced_operator = ReducedFitProblem(input_state) |
| 205 | + return alternating_update( |
| 206 | + reduced_operator, |
| 207 | + init; |
| 208 | + updater = contract_operator_state_updater, |
| 209 | + kwargs..., |
| 210 | + ) |
| 211 | +end |
| 212 | + |
| 213 | + |
| 214 | +function fit( |
| 215 | + input_states::AbstractVector{MPS}, |
| 216 | + init::MPS; |
| 217 | + coeffs::AbstractVector{<:Number} = ones(Int, length(input_states)), |
| 218 | + kwargs..., |
| 219 | +) |
| 220 | + links = ITensors.sim.(linkinds(init)) |
| 221 | + init = replaceinds(linkinds, init, links) |
| 222 | + reduced_operator = ReducedFitMPSsProblem(input_states, coeffs) |
| 223 | + return alternating_update( |
| 224 | + reduced_operator, |
| 225 | + init; |
| 226 | + updater = contract_operator_state_updater, |
| 227 | + kwargs..., |
| 228 | + ) |
| 229 | +end |
| 230 | + |
| 231 | +function fit( |
| 232 | + input_states::AbstractVector{MPO}, |
| 233 | + init::MPO; |
| 234 | + coeffs::AbstractVector{<:Number} = ones(Int, length(input_states)), |
| 235 | + kwargs..., |
| 236 | +) |
| 237 | + :MPO |
| 238 | + to_mps(Ψ::MPO) = MPS([x for x in Ψ]) |
| 239 | + |
| 240 | + res = fit(to_mps.(input_states), to_mps(init); coeffs = coeffs, kwargs...) |
| 241 | + return MPO([x for x in res]) |
| 242 | +end |
0 commit comments