algorithm and example in docs

docs/src/algorithm.md

@@ -5,7 +5,7 @@ Let us consider a classical Ising Hamiltonian
 ```math
 H(\underline{s}_N) =  \sum_{\langle i, j\rangle \in \mathcal{E}} J_{ij} s_i s_j + \sum_{i =1}^N J_{ii} s_i
 ```
-where $\underline{s}_N$ denotes a particular configuration of $N$ binary variables $s_i=\pm 1$. Non-zero couplings $J_{ij} \in \mathbb{R}$ are input parameters of a given problem instance and form a connectivity graph $\mathcal{E}$.
+where $\underline{s}_N$ denotes a particular configuration of $N$ binary variables $s_i=\pm 1$. More generally, we  mark sub-configurations of the first $k$ variables as $\underline{s}_k = (s_1, s_2, \ldots, s_k)$. Non-zero couplings $J_{ij} \in \mathbb{R}$ are input parameters of a given problem instance and form a connectivity graph $\mathcal{E}$.
 
 ## Graphs with large unit cells
 We assume that graph $\mathcal{E}$ forms a quasi-2D lattice. In real life applications such graphs have large unit cells approaching 24 spins. `SpinGlassPEPS.jl` allows for unit cells containing multiple spins. 
@@ -16,18 +16,24 @@ We assume that graph $\mathcal{E}$ forms a quasi-2D lattice. In real life applic
 ```@raw html
 <img src="../images/lattice.pdf" width="200%" class="center"/>
 ```
-Next step in the algorithm is to build clustered Hamiltonian, in which the Ising problem translates to:
+In order to adress this three types of geometries using tensor networks, we represent the problem as a clustered Hamiltonian. To that end we group together sets of variables. In this framework Ising problem translates to:
 ```math
-H(\underline{x}{\bar{N}}) = \sum_{\langle m,n\rangle \in \mathcal{F}} E_{x_m x_n} + \sum_{n=1}^{\bar{N}} E_{x_n}
+H(\underline{x}_{\bar{N}}) = \sum_{\langle m,n\rangle \in \mathcal{F}} E_{x_m x_n} + \sum_{n=1}^{\bar{N}} E_{x_n}
 ```
-$\mathcal{F}$ forms a 2D graph, where we indicate nearest-neighbour interactions with blue lines, and diagonal connections with green lines.
+where $\mathcal{F}$ forms a 2D graph, in which we indicate nearest-neighbour interactions with blue lines, and diagonal connections with green lines in the picture above.
 Each $x_n$ takes $d$ values with  $d=2^4$ for square diagonal, $d=2^{24}$ for Pegasus and $2^{16}$ for Zephyr geometry. 
 $E_{x_n}$ is an intra-node energy of the corresponding binary-variables configuration, and $E_{x_n x_m}$ is inter-node energy.
 
-After creating the clustered Hamiltonian, we can turn it into a PEPS tensor network as shown in the figure below. In all lattices we support, the essential components resembling unit cells are represented by red circles. Spins in adjacent clusters interacted with each other, which is depicted by blue squares. Additionally, we permit diagonal interactions between next nearest neighbors, mediated by green squares.
-
-```@raw html
-<img src="../images/pepstn.pdf" width="200%" class="center"/>
+## Calculating conditional probabilities
+We assume that finding low energy states is equivalent to finding most probable states.
+We represent the probability distribution as a PEPS tensor network.
+```math
+    p(\underline{x}_{\bar{N}}) = \frac{1}{Z} \exp{(-\beta H(\underline{x}_{\bar{N}}))}
+```
+where $Z$ is a partition function and $\beta$ is inverse temperature. Once the PEPS tensor network is constructed, the probability distribution can be obtained by approximately contracting the tensor network, which is described in more details in subsection Tensor network contractions for optimization problems (TODO: link).
+Subsequently, we select only the configurations with the highest marginal probabilities
+```math
+    p(\underline{x}_{n+1}) = p(x_{n+1} | \underline{x}_{n}) \times p(\underline{x}_{n})
 ```
 
 ## Branch and bound search
@@ -37,13 +43,15 @@ By employing branch and bound search strategy iteratively row after row, we addr
 <img src="../images/bb.pdf" width="200%" class="center"/>
 ```
 
-## Calculating conditional probabilities
-
-In order to indentify most probable states we need to calculate the conditional probabilities. Conditional probabilities are obtained by contracting a PEPS tensor network, which, although an NP-hard problem, can be computed approximately. The approach utilized is boundary MPS-MPO, which involves contracting a tensor network row by row and truncating the bond dimension.
+## Tensor network contractions for optimization problems
+Branch and bound search relies on the calculation of conditional probabilities. To that end, we use tensor network techniques. Conditional probabilities are obtained by contracting a PEPS tensor network, which, although an NP-hard problem, can be computed approximately. The approach utilized is boundary MPS-MPO, which involves contracting a tensor network row by row and truncating the bond dimension.
 
 ```@raw html
 <img src="../images/prob.pdf" width="150%" class="center"/>
 ```
+```@raw html
+<img src="../images/explain.pdf" width="75%" class="center"/>
+```
 
 ## References & Related works
 

docs/src/intro.md

@@ -1,19 +1,10 @@
 # Getting started
 Before providing the documentation of the offered functionality, it is good to demonstrate exactly what the package does.
 
-# Basic example
-In this example, we demonstrate how to use the `SpinGlassPEPS.jl` package to obtain a low-energy spectrum search for a spin glass Hamiltonian defined on a square lattice with diagonal interactions on 100 spins.
+## Basic example
+In this example, we demonstrate how to use the `SpinGlassPEPS.jl` package to obtain a low-energy spectrum for a spin glass Hamiltonian defined on a square lattice with diagonal interactions on 100 spins. Let's look at an entire, working code that will do this calculation then discuss the main steps.
 
-The package is used to explore various strategies for solving the problem, and it provides functionalities for performing Hamiltonian clustering, belief propagation, and low-energy spectrum searches using different MPS (Matrix Product State) strategies.
-
-First, we set up the problem by defining the lattice and specifying various parameters such as temperature (β), bond dimension, and search parameters. We also create a clustered Hamiltonian using the specified lattice and perform clustering calculations.
-
-Next, we select the MPS strategy (in this case, Zipper) and other parameters for the network contractor. We create a PEPS (Projected Entangled Pair State) network and initialize the contractor with this network, along with the specified parameters.
-
-Finally, we perform a low-energy spectrum search using the initialized contractor, exploring different branches of the search tree. The example showcases how `SpinGlassPEPS.jl` can be utilized to find the lowest energy configurations for a spin glass system.
-
-
-```@example
+```@julia
 using SpinGlassEngine, SpinGlassTensors, SpinGlassNetworks, SpinGlassPEPS
 using SpinGlassExhaustive
 using Logging
@@ -38,7 +29,7 @@ num_states = 20 # The maximum number of states to be considered during the searc
 
 # MpsParameters 
 bond_dim = 12 # Bond dimension
-max_num_sweeps = 10 # Maximal number of sweeps durion variational compression
+max_num_sweeps = 10 # Maximal number of sweeps during variational compression
 tol_var = 1E-16 # The tolerance for the variational solver used in MPS optimization
 tol_svd = 1E-16 # The tolerance used in singular value decomposition (SVD)
 iters_svd = 2 # The number of iterations to perform in SVD computations
@@ -51,7 +42,12 @@ graduate_truncation = :graduate_truncate # Gradually truncates MPS
 Strategy = Zipper # Strategy to optimize MPS
 transform = rotation(0) # Transformation of the lattice
 Layout = GaugesEnergy # Way of decomposition of the network into MPS
-Sparsity = Sparse # Use sparse mode, when tensors are sparse
+Sparsity = Sparse # Use sparse mode, when tensors are large
+
+# Store parameters in structures
+params = MpsParameters(bond_dim, tol_var, max_num_sweeps, 
+                        tol_svd, iters_svd, iters_var, dtemp_mult, method)
+search_params = SearchParameters(num_states, δp)
 
 # Create Ising graph
 ig = ising_graph(instance)
@@ -61,13 +57,160 @@ cl_h = clustered_hamiltonian(
     spectrum = full_spectrum,
     cluster_assignment_rule=super_square_lattice((m, n, t))
 )
-# Store parameters in structures
+# Build tensor network
+net = PEPSNetwork{SquareCrossSingleNode{Layout}, Sparsity}(m, n, cl_h, transform)
+ctr = MpsContractor{Strategy, NoUpdate}(net, [β], graduate_truncation, params; onGPU=onGPU)
+# Solve using branch and bound search
+sol_peps, s = low_energy_spectrum(ctr, search_params, merge_branches(ctr))
+```
+
+## Steps of the code
+
+The first line
+```@julia
+instance = "$(@__DIR__)/../src/instances/square_diagonal/5x5/diagonal.txt"
+```
+reads instance that is provided in txt format.
+
+Next line defines the problem size
+```@julia
+m, n, t = 5, 5, 4
+```
+In this example, we have number of columns and row, `m` and `n` respectively, equal 5. Parameter `t` tells how many spins are creating a cluster.
+
+`SpinGlassPEPS.jl` enables to perform calculations not only on CPU, but also on GPU. If you want to switch on GPU mode, then type
+```@julia
+onGPU = true
+```
+
+The next part of the code contains parmeters which user should provide before starting the calculations.
+The main parameter is temperature, given as the inverse of temperature.
+```@julia
+β = 1.0
+```
+A higher `β` lets us focus more on low-energy states, but it might make the numerical stability of tensor network contraction a bit shaky. Figuring out the best β depends on the problem, and we might need to try different values in experiments for various instances.
+
+Subsequently, the user can input parameters that will be utilized in exploring the state space, such as the cutoff probability for terminating the search `δp` and the maximum number of states considered during the search (`num_states`).
+```@julia
+# Search parameters
+δp = 0 # The cutoff probability for terminating the search
+num_states = 20 # The maximum number of states to be considered during the search
+```
+
+Another group of parameters describes the method of contracting the network using the boundary MPS-MPO approach.
+```@julia
+bond_dim = 12 # Bond dimension
+max_num_sweeps = 10 # Maximal number of sweeps during variational compression
+tol_var = 1E-16 # The tolerance for the variational solver used in MPS optimization
+tol_svd = 1E-16 # The tolerance used in singular value decomposition (SVD)
+iters_svd = 2 # The number of iterations to perform in SVD computations
+iters_var = 1 # The number of iterations for variational optimization
+dtemp_mult = 2 # A multiplier for the bond dimension
+method = :psvd_sparse # The SVD method to use
+```
+
+We can also choose the arrangement of tensors forming the MPO (`Layout`) and the strategy to optimize boundary MPS (`Strategy`). The user also has the decision-making authority on whether the MPS will be truncated in a gradual manner (`graduate_truncation`). We can initiate our algorithm from various starting points. The parameter responsible for this choice is `transform`, which allows for the rotation and reflection of the tensor network, consequently altering the starting point for the exploration. Last, but not least, the user can also choose whether to use the `Sparse` or `Dense` mode. This choice should depend on the size of the unit cell in the specific problem at hand.
+```@julia
+Layout = GaugesEnergy # Way of decomposition of the network into MPO
+Strategy = Zipper # Strategy to optimize MPS
+graduate_truncation = :graduate_truncate # Gradually truncates MPS
+transform = rotation(0) # Transformation of the lattice
+Sparsity = Sparse # Use sparse mode, when tensors are large
+```
+
+The parameters provided by the user are then stored in data structures `MpsParameters` and `SearchParameters`.
+```@julia
 params = MpsParameters(bond_dim, tol_var, max_num_sweeps, 
                         tol_svd, iters_svd, iters_var, dtemp_mult, method)
 search_params = SearchParameters(num_states, δp)
-# Build tensor network
+```
+
+With this prepared set of parameters, we are ready for the actual computations. The first step is to create the Ising graph. In the Ising graph, nodes are formed by the spin positions, and interactions between them are represented by the edges.
+```@julia
+ig = ising_graph(instance)
+```
+
+Next, we need to translate our problem into a clustered problem, where several spins form an unit cell. This is achieved through the use of the `clustered_hamiltonian` function.
+```@julia
+cl_h = clustered_hamiltonian(
+    ig,
+    spectrum = full_spectrum,
+    cluster_assignment_rule=super_square_lattice((m, n, t))
+)
+```
+
+The next part of the code builds the PEPS tensor network compatible with previously defined clustered Hamiltonian.
+```@julia
 net = PEPSNetwork{SquareCrossSingleNode{Layout}, Sparsity}(m, n, cl_h, transform)
+```
+
+In order to start calculating conditional probabilities we need to define `MpsContractor` structure which stores the elements and information necessary for contracting the tensor network and, consequently, calculating the probability of a given configuration. 
+```@julia
 ctr = MpsContractor{Strategy, NoUpdate}(net, [β], graduate_truncation, params; onGPU=onGPU)
-# Solve using branch and bound search
-# sol_peps, s = low_energy_spectrum(ctr, search_params, merge_branches(ctr))
+```
+
+Finally the call
+```@julia
+sol_peps, s = low_energy_spectrum(ctr, search_params, merge_branches(ctr))
+```
+which runs branch and bound algorithm included in `SpinGlassPEPS.jl` It is actual solver, which iteratively explores the state space in search of the most probable states. The probabilities of a given configuration are calculated approximately through the contractions of the tensor network.
+
+## Expected output
+The function `low_energy_spectrum`, as its output, provides a wealth of information, which we will briefly discuss now.
+```@julia
+sol_peps, schmidt_val = low_energy_spectrum(ctr, search_params, merge_branches(ctr))
+```
+It returns a Solution-type structure (`sol_peps`) and Schmidt values `schmidt_val`. In the `Solution` structure, the following information is recorded:
+* energies - a vector containing the energies of the discovered states
+If you want to display it, you can type:
+```@julia
+println(sol_peps.energies)
+```
+In this case, you should obtain:
+```@julia
+[-215.23679958927175]
+```
+* states - a vector of cluster state configurations corresponding to the energies
+```@julia
+println(sol_peps.states)
+```
+```@julia
+println([[4, 4, 2, 16, 9, 8, 8, 2, 3, 8, 7, 1, 4, 10, 9, 2, 11, 2, 1, 2, 11, 3, 11, 8, 3]])
+```
+* probabilities - the probabilities associated with each discovered state
+```@julia
+println(sol_peps.probabilities)
+```
+```@julia
+[-5.637640487579043]
+```
+* degeneracy
+```@julia
+println(sol_peps.degeneracy)
+```
+```@julia
+[2]
+```
+* largest discarded probability - largest probability below which states are considered to be discarded
+```@julia
+println(sol_peps.largest_discarded_probability)
+```
+```@julia
+-4.70579014404923
+```
+* Schmidt values - a dictionary containing Schmidt spectra for each MPS
+```@julia
+println(schmidt_val)
+```
+which are
+```@julia
+Dict{Any, Any}(5 => Any[1.0, 1.0, 0.0035931343796194565, 0.0015050888555964259, 0.0007184752751868924, 5.2741877514519126e-5, 4.137035816131772e-5, 0.00040017490729592366, 0.00021874495320028077, 6.827849766898342e-5], 4 => Any[1.0, 1.0, 1.704951518711975e-5, 0.00037890798675182353, 0.00011310427642297989, 0.001014257142680146, 0.0012672631937840461, 0.0005487312667512858, 0.0006741839581781018, 0.00012017531445170455], 6 => Any[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 2 => Any[1.0, 1.0, 0.0001405854472578707, 0.0012280075890260514, 0.001177462193268373, 0.0029570115655969827, 0.002997829968910592, 0.0011163442379909382, 0.0010056280784881478, 0.00026431187613365595], 3 => Any[1.0, 1.0, 1.864183962070951e-5, 0.006059161388679921, 0.006793028602573968, 0.012337242616802302, 0.011721497080177857, 0.013791830543357657, 0.020430181282353188, 0.014653186648427675])
+```
+* statistics - a possible warning sign, with values ranging from [0, 2]. A nonzero value signifies that certain conditional probabilities derived from tensor network contraction were negative, suggesting a lack of numerical stability during contraction. Here we display the worst-case scenario.
+```@julia
+println(minimum(values(ctr.statistics)))
+```
+The output should give you
+```@julia
+0.0
 ```

docs/src/sge/api.md

@@ -88,8 +88,6 @@ branch_probability
 exact_conditional_probability
 branch_solution
 gauges_list
-SquareCrossDoubleNode
-SquareSingleNode
 branch_energies
 _equalize
 nodes_search_order_Mps

docs/src/sge/params.md

@@ -31,12 +31,27 @@ The latter, referred to as sparsity, plays a pivotal role in manipulation on lar
 
 # Geometry
 * `SquareSingleNode`
+```@docs
+SquareSingleNode
+```
+
 * `SquareDoubleNode`
+```@docs
+SquareDoubleNode
+```
+
 * `SquareCrossSingleNode`
+```@docs
+SquareCrossSingleNode
+```
+
 * `SquareCrossDoubleNode`
+```@docs
+SquareCrossDoubleNode
+```
 
 # Layout 
-`SpinGlassPEPS.jl` allows for different decompositions of the network into MPS:
+`SpinGlassPEPS.jl` allows for different decompositions of the network into MPOs:
 * `GaugesEnergy`
 * `EnergyGauges`
 * `EngGaugesEng`
@@ -48,7 +63,9 @@ For complex problems, the solution may depend on the choice of decomposition.
 
 # Lattice transformations
 Our package offers users the ability to undergo diverse transformations of PEPS network. Notably, users can apply `rotations`, occurring in multiples of $\frac{\pi}{2}$ radians, and `reflections` along various axes. These transformations include rotations and reflections around the horizontal (x), vertical (y), diagonal, and antidiagonal axes. Transformations are used to contract PEPS and perform search starting from different sites of the lattice. 
-
+```@raw html
+<img src="../images/transform.pdf" width="200%" class="center"/>
+```
 ```@docs
 all_lattice_transformations
 rotation

docs/src/sge/search.md

@@ -1,5 +1,5 @@
 # Branch and bound search
-Here you find fhe main function of the package, which is an actual solver.
+Here you find the main function of the package, which is an actual solver.
 
 ```@docs
 low_energy_spectrum
@@ -10,6 +10,7 @@ Results of the branch and bound search are stored in a Solution structure.
 ```@docs
 Solution
 ```
+
 # Droplet search
 `SpinGlassPEPS.jl` offers the possibility not only finding low lying energy states, but also droplet excitations. In order to search for droplets, one need to choose the option `SingleLayerDroplets` in `merge_branches`.
 ```@docs