Make struct mutable and other improvements

.github/workflows/CI.yml

@@ -18,8 +18,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
           - '1.7'
+          - '1.8'
+          - '1.9'
           - '1.10'
           - 'nightly'
         os:
@@ -48,3 +49,4 @@ jobs:
       - uses: codecov/codecov-action@v4
         with:
           files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}

.gitignore

@@ -2,3 +2,4 @@
 *.jl.cov
 *.jl.mem
 /Manifest.toml
+benchmark/tune.json

Project.toml

@@ -1,20 +1,18 @@
 name = "UniformIsingModels"
 uuid = "91393fb9-e38b-4a31-9409-aa579b4163ad"
 authors = ["Stefano Crotti <[email protected]> and contributors"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
-OffsetArrays = "1"
-UnPack = "1"
 LogExpFunctions = "0.3"
-julia = "1"
+OffsetArrays = "1"
+julia = "1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -3,21 +3,24 @@
 [![Build Status](https://github.com/stecrotti/UniformIsingModels.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/stecrotti/UniformIsingModels.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/stecrotti/UniformIsingModels.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/stecrotti/UniformIsingModels.jl)
 
-A fully-connected ferromagnetic Ising model with uniform coupling strength, described by a Boltzmann distribution
+A fully-connected ferromagnetic [Ising model](https://en.wikipedia.org/wiki/Ising_model) with uniform coupling strength, described by a Boltzmann distribution
 
->![equation](https://latex.codecogs.com/svg.image?p(\boldsymbol\sigma|J,&space;\boldsymbol{h},&space;\beta)&space;=&space;\frac{1}{Z_{J,&space;\boldsymbol{h},&space;\beta}}\exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j&space;&plus;\sum_{i=1}^Nh_i\sigma_i\right)\right],\quad\boldsymbol\sigma\in\\{-1,1\\}^N)
+$p(\boldsymbol{\sigma}) = \frac{1}{Z} \exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right],\quad \boldsymbol{\sigma}\in\{-1,1\}^N $
 
 is exactly solvable in polynomial time.
 
 
-| Quantity | Cost          |
-| ------------- | ----------- |
-| Normalization, Free energy      |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2)) |
-| Sample a configuration      |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2)) |
-| Average energy, Entropy |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2))  |
-| Distribution of the sum of the N spins | ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2))     |
-| Site magnetizations     | ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^3))     |
-| Pair magnetizations, Correlations |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^5))  |
+| Quantity | Expression | Cost          |
+| ------------- | ----------| ----------- |
+| Normalization | $Z=\sum\limits_{\boldsymbol{\sigma}}\exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right]$ | $\mathcal O (N^2)$ |
+| Free energy | $F = -\frac{1}{\beta}\log Z$ | $\mathcal O (N^2)$ |
+| Sample a configuration | $\boldsymbol{\sigma} \sim p(\boldsymbol{\sigma})$ | $\mathcal O (N^2)$ |
+| Average energy | $U = \sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\left[-\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right]$ | $\mathcal O (N^2)$ |
+| Entropy | $S = -\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\log p(\boldsymbol{\sigma})$ | $\mathcal O (N^2)$ |
+| Distribution of the sum of the N spins | $p_S(s)=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\delta\left(s-\sum_{i=1}^N\sigma_i\right)$ | $\mathcal O (N^2)$ |
+| Site magnetizations     | $m_i=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\sigma_i,\quad\forall i\in\{1,2,\ldots,N\}$ | $\mathcal O (N^3)$ |
+| Correlations     | $r_{ij}=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\sigma_i\sigma_j,\quad\forall j\in\{1,2,\ldots,N\},i<j$ | $\mathcal O (N^5)$ |
+
 
 ## Example
 ```
@@ -37,7 +40,7 @@ x = UniformIsing(N, J, h, β)
 Compute stuff
 ```
 # normalization and free energy
-Z = exp(x.logZ)
+Z = normalization(x)
 F = free_energy(x)
 
 # energy and probability of a configuration
@@ -60,10 +63,9 @@ U = avg_energy(x)
 # entropy
 S = entropy(x)
 
-# pairwise magnetizations <σᵢσⱼ> and correlations
-p = pair_magnetizations(x)
-c = correlations(x)
-
+# correlations <σᵢσⱼ> and covariances <σᵢσⱼ>-<σᵢ><σⱼ>
+p = correlations(x)
+c = covariances(x)
 ```
 
 ## Notes

benchmark/benchmarks.jl

@@ -0,0 +1,20 @@
+using BenchmarkTools
+using UniformIsingModels
+
+SUITE = BenchmarkGroup()
+
+N = 100
+J = 0.5
+h = 1.2 * randn(N)
+β = 2.3
+x = UniformIsing(N, J, h, β)
+
+SUITE["constructor"] = BenchmarkGroup()
+SUITE["constructor"]["constructor"] = @benchmarkable UniformIsing($N, $J, $h, $β)
+
+SUITE["observables"] = BenchmarkGroup()
+SUITE["observables"]["normalization"] = @benchmarkable normalization(x)
+SUITE["observables"]["entropy"] = @benchmarkable entropy(x)
+SUITE["observables"]["site_magnetizations"] = @benchmarkable site_magnetizations(x)
+# SUITE["observables"]["pair_magnetizations"] = @benchmarkable pair_magnetizations(x)
+SUITE["observables"]["sum_distribution"] = @benchmarkable sum_distribution(x)

src/UniformIsingModels.jl

@@ -1,193 +1,20 @@
 module UniformIsingModels
 
-import OffsetArrays: OffsetVector, fill
-import LinearAlgebra: dot
-import Base: show
-import UnPack: @unpack
-import Random: GLOBAL_RNG, AbstractRNG
-import LogExpFunctions: logsumexp, logaddexp, logsubexp
-
-export UniformIsing, energy, normalization, pdf, 
+using OffsetArrays: OffsetVector, fill
+using LinearAlgebra: dot
+using Random: default_rng, AbstractRNG
+using LogExpFunctions: logsumexp, logaddexp
+
+export UniformIsing, nvariables, variables, recompute_partials!,
+        energy, lognormalization, normalization, pdf,
+        avg_energy, entropy, free_energy, 
         site_magnetizations!, site_magnetizations,
-        pair_magnetizations!, pair_magnetizations,
         correlations!, correlations,
+        covariances, covariances,
         sum_distribution!, sum_distribution,
-        sample!, sample,
-        avg_energy, entropy, free_energy
+        sample!, sample
 
 include("accumulate.jl")
-
-struct UniformIsing{T<:Real, U<:OffsetVector}
-    N    :: Int                # number of spins
-    J    :: T                  # uniform coupling strength
-    h    :: Vector{T}          # external fields
-    β    :: T                  # inverse temperature
-    L    :: OffsetVector{U, Vector{U}} # partial sums from the left
-    R    :: OffsetVector{U, Vector{U}} # partial sums from the right
-    dLdB ::  OffsetVector{U, Vector{U}} # Derivative of L wrt β
-    logZ :: T                  # normalization
-    function UniformIsing(N::Int, J::T, h::Vector{T}, β::T=1.0) where T
-        @assert length(h) == N
-        @assert β ≥ 0
-        # L = accumulate_left(h, β)
-        R = accumulate_right(h, β)
-        L, dLdB = accumulate_d_left(h, β)
-        logZ = logsumexp( β*J/2*(s^2/N-1) + L[N][s] for s in -N:N )
-        new{T, eltype(L)}(N, J, h, β, L, R, dLdB, logZ)
-    end
-end
-
-function show(io::IO, x::UniformIsing)
-    @unpack N, J, h, β = x
-    println(io, "UniformIsing with N = $N variables at temperature β = $β")
-end
-
-function energy(x::UniformIsing, σ)
-    s = sum(σ)
-    f = dot(σ, x.h)
-    -( x.J/2*(s^2/x.N-1) + f ) 
-end
-
-pdf(x::UniformIsing, σ) = exp(-x.β*energy(x, σ) - x.logZ)
-
-free_energy(x::UniformIsing) = -1/x.β*x.logZ
-
-function sample_spin(rng::AbstractRNG, p::Real)
-    @assert 0 ≤ p ≤ 1
-    r = rand(rng)
-    r < p ? 1 : -1
-end
-
-# return a sample along with its probability
-function sample!(rng::AbstractRNG, σ, x::UniformIsing)
-    @unpack N, J, h, β, L, R, logZ = x
-    a = 0.0; b = 0
-    f(s) = β*J/2*(s^2/N-1)
-    for i in 1:N
-        tmp_plus = tmp_minus = 0.0
-        for s in -N:N
-            tmp_plus += exp(f(b+1+s) + R[i+1][s])
-            tmp_minus += exp(f(b-1+s) + R[i+1][s])
-        end
-        p_plus = exp(β*h[i]) * tmp_plus
-        p_minus = exp(-β*h[i]) * tmp_minus
-        p_i = p_plus / (p_plus + p_minus)
-        σi = sample_spin(rng, p_i)
-        σ[i] = σi
-        a += h[i]*σi
-        b += σi
-    end
-    p = exp(f(b) + β*a - logZ)
-    @assert a == dot(h, σ); @assert b == sum(σ)
-    σ, p
-end
-sample!(σ, x::UniformIsing) = sample!(GLOBAL_RNG, σ, x)
-sample(rng::AbstractRNG, x::UniformIsing) = sample!(rng, zeros(Int, x.N), x)
-sample(x::UniformIsing) = sample(GLOBAL_RNG, x)
-
-# first store in `p[i]` the quantity log(p(σᵢ=+1)), then transform at the end 
-function site_magnetizations!(p, x::UniformIsing)
-    @unpack N, J, h, β, L, R, logZ = x
-    f(s) = β*J/2*(s^2/N-1)
-    for i in eachindex(p)
-        p[i] = -Inf
-        for sL in -N:N
-            for sR in -N:N
-                s = sL + sR
-                p[i] = logaddexp(p[i], f(s+1) + β*h[i] + L[i-1][sL] + R[i+1][sR])
-            end
-        end
-        # include normalization
-        p[i] = exp(p[i] - logZ)
-        # transform form p(+) to m=2p(+)-1
-        p[i] = 2*p[i] - 1
-    end
-    p
-end
-function site_magnetizations(x::UniformIsing{T,U}) where {T,U} 
-    site_magnetizations!(zeros(T,x.N), x)
-end
-
-# distribution of the sum of all variables as an OffsetVector
-function sum_distribution!(p, x::UniformIsing)
-    @unpack N, J, β, L, logZ = x
-    for s in -N:N
-        p[s] = exp( β*J/2*(s^2/N-1) + L[N][s] - logZ ) 
-    end
-    p
-end
-function sum_distribution(x::UniformIsing{T,U}) where {T,U} 
-    p = fill(zero(T), -x.N:x.N)
-    sum_distribution!(p, x)
-end
-
-function avg_energy(x::UniformIsing{T}) where T
-    @unpack N, J, h, β, L, dLdB, logZ = x
-    Zt = sum( exp( β*J/2*(s^2/N-1) + L[N][s]) * (J/2*(s^2/N-1)+dLdB[N][s]) for s in -N:N)
-    -exp(log(Zt) - logZ)
-end
-
-entropy(x::UniformIsing; kw...) = x.β * (avg_energy(x; kw...) - free_energy(x)) 
-
-function pair_magnetizations!(m, x::UniformIsing{T,U};
-        M = accumulate_middle(x.h, x.β)) where {T,U}
-    @unpack N, J, h, β, L, R, logZ = x
-    Z = exp(logZ)
-    f(s) = β*J/2*(s^2/N-1)
-    for i in 1:N
-        # j = i
-        m[i,i] = 1
-        # j = i+1
-        j = i + 1
-        j > N && break
-        m[i,j] = 0
-        for sL in -N:N
-            for sR in -N:N
-                s = sL + sR
-                m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
-                            exp( f(s-2) - β*(h[i]+h[j]) ) -
-                            exp( f(s)   + β*(h[i]-h[j]) ) - 
-                            exp( f(s)   - β*(h[i]-h[j]) )   ) * 
-                            exp(L[i-1][sL] + R[j+1][sR]) 
-            end
-        end
-        m[i,j] /= Z; m[j,i] = m[i,j]
-        # j > i + 1
-        for j in i+2:N
-            m[i,j] = 0
-            for sM in -N:N
-                for sL in -N:N
-                    for sR in -N:N
-                        s = sL + sM + sR 
-                        m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
-                                  exp( f(s-2) - β*(h[i]+h[j]) ) -
-                                  exp( f(s)   + β*(h[i]-h[j]) ) - 
-                                  exp( f(s)   - β*(h[i]-h[j]) )   ) *
-                                  exp(L[i-1][sL] + M[i+1,j-1][sM] + R[j+1][sR])
-                    end 
-                end
-            end
-            m[i,j] /= Z; m[j,i] = m[i,j]
-        end
-    end
-    m
-end
-function pair_magnetizations(x::UniformIsing{T,U}; kw...) where {T,U} 
-    pair_magnetizations!(zeros(T,x.N,x.N), x; kw...)
-end
-
-function correlations!(c, x::UniformIsing;
-        m = site_magnetizations(x), p = pair_magnetizations(x))
-    N = x.N
-    for i in 1:N
-        for j in 1:N
-            c[i,j] = p[i,j] - m[i]*m[j]
-        end
-    end
-    c
-end
-function correlations(x::UniformIsing{T,U}; kw...) where {T,U} 
-    correlations!(zeros(T,x.N,x.N), x; kw...)
-end
+include("uniform_ising.jl")
 
 end # end module

src/accumulate.jl

@@ -17,30 +17,30 @@ function accumulate_left(h, β)
     accumulate_left!(L, h, β)
 end
 
-function accumulate_d_left!(L, dLdB, h, β)
+function accumulate_d_left!(L, dLdβ, h, β)
     N = length(h)
-    L[0] .= -Inf; L[0][0] = 0.0; dLdB[0] .= 0.0
+    L[0] .= -Inf; L[0][0] = 0.0; dLdβ[0] .= 0.0
     for k in 1:N
         for s in -(k-1):k-1
             L[k][s] = logaddexp( β*h[k] + L[k-1][s-1],
                                 -β*h[k] + L[k-1][s+1] )
-            dLdB[k][s] = (+h[k] + dLdB[k-1][s-1])*exp(+β*h[k]+L[k-1][s-1]) + 
-                       (-h[k] + dLdB[k-1][s+1])*exp(-β*h[k]+L[k-1][s+1])
-            dLdB[k][s] = dLdB[k][s] / exp(L[k][s])
-            isnan(dLdB[k][s]) && (dLdB[k][s] = 0.0)
+            dLdβ[k][s] = (+h[k] + dLdβ[k-1][s-1])*exp(+β*h[k]+L[k-1][s-1]) + 
+                       (-h[k] + dLdβ[k-1][s+1])*exp(-β*h[k]+L[k-1][s+1])
+            dLdβ[k][s] = dLdβ[k][s] / exp(L[k][s])
+            isnan(dLdβ[k][s]) && (dLdβ[k][s] = 0.0)
         end
         L[k][k] = β*h[k] + L[k-1][k-1]
         L[k][-k] = -β*h[k] + L[k-1][-k+1]
-        dLdB[k][k] = h[k] + dLdB[k-1][k-1]
-        dLdB[k][-k] = -h[k] + dLdB[k-1][-k+1]
+        dLdβ[k][k] = h[k] + dLdβ[k-1][k-1]
+        dLdβ[k][-k] = -h[k] + dLdβ[k-1][-k+1]
     end
-    L, dLdB
+    L, dLdβ
 end
 function accumulate_d_left(h, β)
     N = length(h)
     L = OffsetVector([fill(-Inf, -N:N) for i in 0:N], 0:N)
-    dLdB = OffsetVector([fill(0.0, -N:N) for i in 0:N], 0:N)
-    accumulate_d_left!(L, dLdB, h, β)
+    dLdβ = OffsetVector([fill(0.0, -N:N) for i in 0:N], 0:N)
+    accumulate_d_left!(L, dLdβ, h, β)
 end
 
 function accumulate_right!(R, h, β)