Refactoring and perf improvements (#15)

* Fix array expression integration

* Update manuscript

* fix example

* add benchmark code

* update manuscript

* refactor and improve performance of weingarten calculation

* remove bigint

* almost working refacctored code

* workings tests

* fix trace output

* bump version

* update docs

* fix symbolic traces

* fix docs

* update benchmarks in manuscript

* fix permutation group performance

* remove shows

Author: lpawela

Project.toml

@@ -1,6 +1,6 @@
 name = "IntU"
 uuid = "70c8b9a1-b4c5-4d7a-a56e-8e7e6495d68b"
-version = "0.4.0"
+version = "0.5.0"
 authors = ["Łukasz Pawela", "Zbigniew Puchała"]
 
 [deps]

README.md

@@ -12,12 +12,26 @@ For detailed documentation, please visit [iitis.github.io/IntU.jl](https://iitis
 
 To introduce the main functionality of IntU, consider the problem of averaging
 $|U_{i,j}|^2$ over the unitary group, i.e., computing $\int dU |U_{i,j}|^2 =
-\int dU U_{i,j} U_{i,j}^*$.
+\int dU U_{i,j} U_{k,l}^* dU$.
 
 While numerical approaches (like sampling random matrices) can estimate this,
 they are slow and approximate. IntU provides the **exact** analytic result
 instantly, even for symbolic dimensions.
 
+**New Feature**: You can now integrate matrix-valued expressions directly!
+```julia
+using IntU, Symbolics, LinearAlgebra
+
+d = 3
+@variables U[1:d, 1:d]::Complex
+measure = dU(U, d)
+
+# Integrate the matrix expression U * U'
+# This performs element-wise integration automatically
+res = integrate(collect(U * U'), measure)
+# Output: Identity Matrix (I)
+```
+
 ```julia
 using IntU, Symbolics
 
@@ -161,7 +175,7 @@ U_sym = SymbolicMatrix(:U, false, :U) # unitary
 # Compute ∫ tr(U A U† B) dU
 expr = tr(U_sym * A * U_sym' * B)
 integrate(expr, dU(d))
-# Output: tr(A)*tr(B) / d
+# Output: (tr(A)*tr(B)) / d
 ```
 
 ### Stiefel Manifolds
@@ -192,6 +206,38 @@ A = randomITensor(j, i)
 res = integrate([U, A], dU(2))
 ```
 
+## Running Examples and Benchmarks
+
+IntU.jl comes with a comprehensive set of examples and benchmarks.
+
+### Examples
+
+You can run all examples at once using the provided script:
+
+```bash
+sh examples/runexamples.sh
+```
+
+To run a single example, pass the example number to the script. For instance, to run `01_basics.jl`:
+
+```bash
+sh examples/runexamples.sh 01
+```
+
+### Benchmarks
+
+To run all benchmarks:
+
+```bash
+sh benchmarks/runbenchmarks.sh
+```
+
+To run a single benchmark, pass the number to the script:
+
+```bash
+sh benchmarks/runbenchmarks.sh 01
+```
+
 ## Installation
 
 IntU is tested with Julia 1.11 or later. Installation can be done through the

benchmarks/05_symbolic_trace.jl

@@ -4,9 +4,8 @@ using Symbolics
 
 function run_benchmarks()
     @variables d
-    @variables u_dummy[1:1, 1:1]
-    measure = dU(u_dummy, d)
     U = SymbolicMatrix(:U, false, :U)
+    measure = dU(U, d)
     Ud = U'
 
     # Create constants

benchmarks/09_stress_test_2.jl

@@ -41,8 +41,12 @@ doublefactorial_odd(n::Int) = prod(1:2:n)  # n must be odd
 
 # Symbolic zero test (best-effort; also backs up by numeric spot-checks)
 function sym_is_zero(x)
-    xs = Symbolics.simplify(x)
-    return (xs == 0) || (xs === 0)
+    try
+        xs = Symbolics.simplify(x)
+        return (xs == 0) || (xs === 0)
+    catch
+        return false
+    end
 end
 
 function safe_string(x)
@@ -55,7 +59,12 @@ end
 
 # Evaluate symbolic expression at d = val (if d is a Symbolics variable)
 function eval_at(expr, dvar, val::Int)
-    return Symbolics.simplify(Symbolics.substitute(expr, Dict(dvar => val)))
+    subbed = Symbolics.substitute(expr, Dict(dvar => val))
+    try
+        return Symbolics.simplify(subbed)
+    catch
+        return subbed
+    end
 end
 
 function full_simplify(x)
@@ -64,7 +73,11 @@ function full_simplify(x)
         # Canonicalize by expanding then simplifying
         return Symbolics.simplify(Symbolics.expand(x))
     catch
-        return Symbolics.simplify(x)
+        try
+            return Symbolics.simplify(x)
+        catch
+            return x
+        end
     end
 end
 
@@ -148,6 +161,15 @@ U_cases = [
     ),
 ]
 
+function safe_eq(x, y)
+    try
+        val = (x == y)
+        return val isa Bool ? val : false
+    catch
+        return false
+    end
+end
+
 for (name, expr, expected) in U_cases
     got, bm = bench_integrate(expr, μU; samples = samples)
     diff0 = sym_is_zero(got - expected)
@@ -157,7 +179,7 @@ for (name, expr, expected) in U_cases
         g = eval_at(got, d, dv)
         e = eval_at(expected, d, dv)
         numeric_checks[string(dv)] =
-            Dict("got" => safe_string(g), "expected" => safe_string(e), "ok" => (g == e))
+            Dict("got" => safe_string(g), "expected" => safe_string(e), "ok" => safe_eq(g, e))
     end
 
     println("[U(d)] $name")
@@ -234,7 +256,7 @@ for (name, expr, k) in O_cases
             numeric_checks[string(dv)] = Dict(
                 "got" => safe_string(g),
                 "expected" => safe_string(e),
-                "ok" => (g == e),
+                "ok" => safe_eq(g, e),
             )
         end
     end

benchmarks/16_matrix_integration.jl

@@ -0,0 +1,30 @@
+using IntU
+using Symbolics
+using LinearAlgebra
+using BenchmarkTools
+
+function benchmark_matrix_integration(sizes)
+    println("=== Matrix Integration Benchmark (Haar U * U') ===")
+    @variables d
+
+    for N in sizes
+        println("\nBenchmarking Matrix Size N=$N")
+        # Define N x N symbolic matrix
+        # Note: integration dimension is symbolic 'd', but matrix size is fixed N
+        @variables U[1:N, 1:N]::Complex
+        measure = dU(U, d)
+
+        # We need to collect to ensure we pass a Matrix{Num} to integrate
+        expr = collect(U * U')
+        println("  Expression size: $(size(expr))")
+
+        # Warmup
+        integrate(expr, measure)
+
+        # Benchmark
+        t = @benchmark integrate($expr, $measure)
+        display(t)
+    end
+end
+
+benchmark_matrix_integration([2, 3, 4])

docs/src/api.md

@@ -113,4 +113,7 @@ IntU.character_at_id
 IntU.irrep_dimension
 IntU.get_weingarten_orthogonal_data
 IntU.compute_symplectic_contraction
+IntU.weingarten_orthogonal_val_canonical
+IntU.INTEGRATION_RULES
+IntU.measure_info
 ```

docs/src/symbolic_trace.md

@@ -37,7 +37,7 @@ expr = IntU.tr(U * A * U' * A)
 # 4. Integrate
 res = integrate(expr, measure)
 println(res)
-# Output: (tr(A)^2)/d
+# Output: (tr(A)^2) / d
 ```
 
 ## Implementation Details
@@ -62,14 +62,12 @@ Information Theory without getting bogged down in index hell.
 # Product of traces
 expr = tr(U * A) * tr(U' * B)
 integrate(expr, measure)
-# Output: tr(A B) / d
+# Output: tr(A*B) / d
 
 # Sum of traces
 expr_sum = tr(U * A * U') + tr(B)
 integrate(expr_sum, measure)
-# Output: tr(A)/d * tr(I) + tr(B) = tr(A) + tr(B)
+# Output: (tr(A) / d) * tr(I_d) + tr(B) = tr(A) + tr(B)
 ```
 
 The underlying engine handles the "wiring" of indices across multiple trace cycles automatically.
-
-The underlying engine handles the "wiring" of indices across multiple trace cycles automatically.

docs/src/unitary_integration.md

@@ -121,6 +121,20 @@ println(res3)
 # Output: 1
 ```
 
+### 4. Matrix Integration
+
+New in v0.2: You can integrate matrix-valued expressions directly. The function `integrate` will element-wise integrate any `AbstractArray` passed to it.
+
+```julia
+using LinearAlgebra
+
+# Integrate U * U' (should be identity)
+# We collect the symbolic expression to a Matrix{Num} to ensure it's treated as an array of expressions
+expr_mat = collect(U * U')
+res_mat = integrate(expr_mat, measure)
+# Result is the Identity matrix
+```
+
 ## Potential Pitfalls
 
 -   **Symbolic vs Numeric Dimension**: The dimension $d$ can be symbolic.

examples/01_basics.jl

@@ -1,6 +1,6 @@
-# examples/01_basics.jl
 using IntU
 using Symbolics
+using LinearAlgebra
 
 # 1. Define the dimension and variables
 d_val = 3
@@ -25,6 +25,16 @@ result2 = integrate(expr2, measure)
 println("Result: ", result2)
 println("Expected: ", 1//(d_val^2 - 1))
 
-# 4. Example: General d formula (symbolic d is not fully supported yet, 
+# 5. Example: Matrix Integration
+# You can integrate matrix-valued expressions directly.
+# We use collect() to ensure we pass a standard Julia Matrix of symbolic numbers.
+println("\n5. Example: Matrix Integration")
+println("Integrating U * U' (should be Identity)")
+expr_mat = collect(U * U')
+result_mat = integrate(expr_mat, measure)
+println("Result[1,1]: ", result_mat[1, 1])
+println("Result is Identity? ", result_mat == I)
+
+# 6. Example: General d formula (symbolic d is not fully supported yet,
 # but we can show it matches the theoretical value for a specific d)
 println("\nNote: Theoretical value involves Weingarten function Wg(1^2, d) = 1/(d^2-1)")

examples/09_real_integration.jl

@@ -31,6 +31,17 @@ println("Result (symbolic sum k=1..3): ", res3)
 println("Value at d=3: ", Symbolics.substitute(res3, Dict(d => 3)))
 
 
+println("\n4. Matrix Integration Example")
+println("Integrating O * O^T (should be Identity)")
+# We need to collect symbolic array to standard Matrix{Num} to ensure element-wise integration
+O_mat = collect(O)
+res_mat = integrate(O_mat * O_mat', mO)
+
+# Check first element
+println("Result[1,1]: ", res_mat[1,1])
+println("Expected: 1")
+
+
 # --- Symplectic Group Sp(d) ---
 println("\n--- Symplectic Group Sp(d) ---")
 println("(Note: Sp(d) integration requires d to be even)")