WIP

Project.toml

@@ -6,6 +6,7 @@ version = "0.3.6"
 [deps]
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FastMPOContractions = "f6e391d2-8ffa-4d7a-98cd-7e70024481ca"
+ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d"

notebook/sum.ipynb

@@ -0,0 +1,89 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using Revise\n",
+    "using Quantics\n",
+    "using ITensors"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "MethodError",
+     "evalue": "MethodError: no method matching tags(::Nothing)\n\nClosest candidates are:\n  tags(!Matched::Index)\n   @ ITensors ~/.julia/packages/ITensors/nUzgj/src/index.jl:204\n",
+     "output_type": "error",
+     "traceback": [
+      "MethodError: no method matching tags(::Nothing)\n",
+      "\n",
+      "Closest candidates are:\n",
+      "  tags(!Matched::Index)\n",
+      "   @ ITensors ~/.julia/packages/ITensors/nUzgj/src/index.jl:204\n",
+      "\n",
+      "\n",
+      "Stacktrace:\n",
+      "  [1] replacebond!(M::MPS, b::Int64, phi::ITensor; normalize::Bool, swapsites::Nothing, ortho::String, which_decomp::Nothing, mindim::Int64, maxdim::Int64, cutoff::Float64, eigen_perturbation::Nothing, svd_alg::Nothing, use_absolute_cutoff::Nothing, use_relative_cutoff::Nothing, min_blockdim::Nothing)\n",
+      "    @ ITensors.ITensorMPS ~/.julia/packages/ITensors/nUzgj/src/lib/ITensorMPS/src/mps.jl:563\n",
+      "  [2] replacebond!\n",
+      "    @ ~/.julia/packages/ITensors/nUzgj/src/lib/ITensorMPS/src/mps.jl:534 [inlined]\n",
+      "  [3] region_update!(nsite_val::Val{2}, reverse_step_val::Val{false}, reduce_operator::Quantics.ReducedSumProblem, state::MPS, b::Int64; updater::typeof(Quantics.add_operator_state_updater), updater_kwargs::@NamedTuple{}, current_time::Nothing, time_step::Nothing, outputlevel::Int64, normalize::Bool, direction::Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}, noise::Int64, which_decomp::Nothing, svd_alg::Nothing, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)\n",
+      "    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/4kskX/src/sweep_update.jl:372\n",
+      "  [4] region_update!(reduce_operator::Quantics.ReducedSumProblem, state::MPS, b::Int64; updater::Function, updater_kwargs::@NamedTuple{}, nsite::Int64, reverse_step::Bool, current_time::Nothing, outputlevel::Int64, time_step::Nothing, normalize::Bool, direction::Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}, noise::Int64, which_decomp::Nothing, svd_alg::Nothing, cutoff::Float64, maxdim::Int64, mindim::Int64, maxtruncerr::Float64)\n",
+      "    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/4kskX/src/sweep_update.jl:196\n",
+      "  [5] sub_sweep_update(direction::Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}, reduce_operator::Quantics.ReducedSumProblem, state::MPS; updater::Function, updater_kwargs::@NamedTuple{}, which_decomp::Nothing, svd_alg::Nothing, sweep::Int64, current_time::Nothing, time_step::Nothing, nsite::Int64, reverse_step::Bool, normalize::Bool, observer!::ITensorTDVP.EmptyObserver, outputlevel::Int64, maxdim::Int64, mindim::Int64, cutoff::Float64, noise::Int64)\n",
+      "    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/4kskX/src/sweep_update.jl:117\n",
+      "  [6] sweep_update(order::ITensorTDVP.TDVPOrder{2, Base.Order.ForwardOrdering()}, reduce_operator::Quantics.ReducedSumProblem, state::MPS; current_time::Nothing, time_step::Nothing, kwargs::@Kwargs{updater::typeof(Quantics.add_operator_state_updater), updater_kwargs::@NamedTuple{}, nsite::Int64, reverse_step::Bool, sweep::Int64, observer!::ITensorTDVP.EmptyObserver, normalize::Bool, outputlevel::Int64, maxdim::Int64, mindim::Int64, cutoff::Float64, noise::Int64})\n",
+      "    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/4kskX/src/sweep_update.jl:35\n",
+      "  [7] macro expansion\n",
+      "    @ ~/.julia/packages/ITensorTDVP/4kskX/src/alternating_update.jl:71 [inlined]\n",
+      "  [8] macro expansion\n",
+      "    @ ./timing.jl:395 [inlined]\n",
+      "  [9] alternating_update(operator::Quantics.ReducedSumProblem, init::MPS; updater::Function, updater_kwargs::@NamedTuple{}, nsweeps::Int64, checkdone::Returns{Bool}, write_when_maxdim_exceeds::Nothing, nsite::Int64, reverse_step::Bool, time_start::Nothing, time_step::Nothing, order::Int64, observer!::ITensorTDVP.EmptyObserver, sweep_observer!::ITensorTDVP.EmptyObserver, outputlevel::Int64, normalize::Bool, maxdim::Int64, mindim::Int64, cutoff::Float64, noise::Int64)\n",
+      "    @ ITensorTDVP ~/.julia/packages/ITensorTDVP/4kskX/src/alternating_update.jl:70\n",
+      " [10] alternating_update\n",
+      "    @ ~/.julia/packages/ITensorTDVP/4kskX/src/alternating_update.jl:25 [inlined]\n",
+      " [11] #variational_add#35\n",
+      "    @ ~/git/Quantics.jl/src/sum.jl:18 [inlined]\n",
+      " [12] top-level scope\n",
+      "    @ ~/git/Quantics.jl/notebook/sum.ipynb:4"
+     ]
+    }
+   ],
+   "source": [
+    "sites = siteinds(\"S=1/2\", 10)\n",
+    "Ψ = random_mps(sites)\n",
+    "\n",
+    "Quantics.variational_add(Ψ; cutoff=1e-10, maxdim=2, nsweeps=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.10.4",
+   "language": "julia",
+   "name": "julia-1.10"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.10.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

src/Quantics.jl

@@ -17,13 +17,15 @@ import LinearAlgebra: I
 using StaticArrays
 
 import FastMPOContractions
+import ITensorTDVP
 
 using EllipsisNotation
 
 function __init__()
 end
 
 include("util.jl")
+include("sum.jl")
 include("tag.jl")
 include("binaryop.jl")
 include("mul.jl")

src/sum.jl

@@ -0,0 +1,131 @@
+using ITensors: ITensor
+using ITensors.ITensorMPS: ITensorMPS, AbstractProjMPO, MPO, MPS
+
+function add_operator_state_updater(operator, init; internal_kwargs)
+  state = ITensor(true)
+  for j in (operator.lpos + 1):(operator.rpos - 1)
+    state *= operator.input_state[j]
+  end
+  return state, (;)
+end
+
+function variational_add(
+  input_state::MPS;
+  init=input_state,
+  kwargs...,
+)
+  reduced_operator = ReducedSumProblem(input_state)
+  return ITensorTDVP.alternating_update(
+    reduced_operator, init; updater=add_operator_state_updater, kwargs...
+  )
+end
+
+"""
+A ReducedSumProblem represents fitting to MPO `operator`
+onto an MPS `input_state` but "projected" by
+the basis of a different MPS `state` (which
+could be an approximation to operator|state>).
+
+As an implementation of the AbstractProjMPO
+type, it supports multiple `nsite` values for
+one- and two-site algorithms.
+```
+     *--*--*-      -*--*--*--*--*--* <state|
+     |  |  |  |  |  |  |  |  |  |  |
+     o--o--o-      -o--o--o--o--o--o |input_state>
+```
+"""
+mutable struct ReducedSumProblem <: AbstractProjMPO
+  lpos::Int
+  rpos::Int
+  nsite::Int
+  input_state::MPS
+  environments::Vector{ITensor}
+end
+
+function ReducedSumProblem(input_state::MPS)
+  lpos = 0
+  rpos = length(input_state) + 1
+  nsite = 2
+  environments = Vector{ITensor}(undef, length(input_state))
+  return ReducedSumProblem(lpos, rpos, nsite, input_state, environments)
+end
+
+function Base.getproperty(reduced_operator::ReducedSumProblem, sym::Symbol)
+  # This is for compatibility with `AbstractProjMPO`.
+  if sym == :LR
+    return getfield(reduced_operator, :environments)
+  end
+  return getfield(reduced_operator, sym)
+end
+
+function Base.copy(reduced_operator::ReducedSumProblem)
+  return ReducedSumProblem(
+    reduced_operator.lpos,
+    reduced_operator.rpos,
+    reduced_operator.nsite,
+    copy(reduced_operator.input_state),
+    copy(reduced_operator.environments),
+  )
+end
+
+Base.length(reduced_operator::ReducedSumProblem) = length(reduced_operator.input_state)
+
+function ITensorMPS.set_nsite!(reduced_operator::ReducedSumProblem, nsite)
+  reduced_operator.nsite = nsite
+  return reduced_operator
+end
+
+function ITensorMPS.makeL!(reduced_operator::ReducedSumProblem, state::MPS, k::Int)
+  # Save the last `L` that is made to help with caching
+  # for DiskProjMPO
+  ll = reduced_operator.lpos
+  if ll ≥ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if lproj is
+    # being moved backward.
+    reduced_operator.lpos = k
+    return nothing
+  end
+  # Make sure ll is at least 0 for the generic logic below
+  ll = max(ll, 0)
+  L = lproj(reduced_operator)
+  while ll < k
+    L =
+      L *
+      reduced_operator.input_state[ll + 1] *
+      dag(state[ll + 1])
+    reduced_operator.environments[ll + 1] = L
+    ll += 1
+  end
+  # Needed when moving lproj backward.
+  reduced_operator.lpos = k
+  return reduced_operator
+end
+
+function ITensorMPS.makeR!(reduced_operator::ReducedSumProblem, state::MPS, k::Int)
+  # Save the last `R` that is made to help with caching
+  # for DiskProjMPO
+  rl = reduced_operator.rpos
+  if rl ≤ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if rproj is
+    # being moved backward.
+    reduced_operator.rpos = k
+    return nothing
+  end
+  N = length(state)
+  # Make sure rl is no bigger than `N + 1` for the generic logic below
+  rl = min(rl, N + 1)
+  R = rproj(reduced_operator)
+  while rl > k
+    R =
+      R *
+      reduced_operator.input_state[rl - 1] *
+      dag(state[rl - 1])
+    reduced_operator.environments[rl - 1] = R
+    rl -= 1
+  end
+  reduced_operator.rpos = k
+  return reduced_operator
+end