From f19dd557da9a44c001a81f0a200c6259c8572dde Mon Sep 17 00:00:00 2001
From: Daniel Wennberg <daniel.wennberg@gmail.com>
Date: Wed, 28 Feb 2024 11:58:42 -0800
Subject: [PATCH 1/2] Improve finite domain transformation tests

Catch substitution bug
---
 test/inf_integral_tests.jl | 53 ++++++++++++++++++++++++--------------
 1 file changed, 34 insertions(+), 19 deletions(-)

diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl
index 9164fd4..c0763e1 100644
--- a/test/inf_integral_tests.jl
+++ b/test/inf_integral_tests.jl
@@ -50,70 +50,85 @@ problems = (
         domain = (@SVector[-Inf, -Inf], @SVector[Inf, 0]),
         solution = 0.5
     ),
-    (; # 5. single-variable infinite limit: Gaussian
+    (; # 5. multi-variate mixed infinite/finite: Gaussian * quadratic
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x[1]) * x[2]^2,
+        domain = (@SVector[-Inf, 1], @SVector[Inf, 4]),
+        solution = 21.0
+    ),
+    (; # 6. multi-variate mixed semi-infinite lower limit/finite: Gaussian * quadratic
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x[1]) * x[2]^2,
+        domain = (@SVector[0.00, 1], @SVector[Inf, 4]),
+        solution = 10.5
+    ),
+    (; # 7. multi-variate mixed semi-infinite upper limit/finite: Gaussian * quadratic
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x[1]) * x[2]^2,
+        domain = (@SVector[-Inf, 1], @SVector[0.00, 4]),
+        solution = 10.5
+    ),
+    (; # 8. single-variable infinite limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (-Inf, Inf),
         solution = 1.0
     ),
-    (; # 6. single-variable flipped infinite limit: Gaussian
+    (; # 9. single-variable flipped infinite limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (Inf, -Inf),
         solution = -1.0
     ),
-    (; # 7. single-variable semi-infinite upper limit: Gaussian
+    (; # 10. single-variable semi-infinite upper limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (0.00, Inf),
         solution = 0.5
     ),
-    (; # 8. single-variable flipped, semi-infinite upper limit: Gaussian
+    (; # 11. single-variable flipped, semi-infinite upper limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (0.00, -Inf),
         solution = -0.5
     ),
-    (; # 9. single-variable semi-infinite lower limit: Gaussian
+    (; # 12. single-variable semi-infinite lower limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (-Inf, 0.00),
         solution = 0.5
     ),
-    (; # 10. single-variable flipped, semi-infinite lower limit: Gaussian
+    (; # 13. single-variable flipped, semi-infinite lower limit: Gaussian
         f = (x, p) -> pdf(Normal(0.00, 1.00), x),
         domain = (Inf, 0.00),
         solution = -0.5
     ),
-    (; # 11. single-variable infinite limit: Lorentzian
+    (; # 14. single-variable infinite limit: Lorentzian
         f = (x, p) -> 1 / (x^2 + 1),
         domain = (-Inf, Inf),
         solution = pi / 1
     ),
-    (; # 12. single-variable shifted, semi-infinite lower limit: Lorentzian
+    (; # 15. single-variable shifted, semi-infinite lower limit: Lorentzian
         f = (x, p) -> 1 / ((x - 2)^2 + 1),
         domain = (-Inf, 2),
         solution = pi / 2
     ),
-    (; # 13. single-variable shifted, semi-infinite upper limit: Lorentzian
+    (; # 16. single-variable shifted, semi-infinite upper limit: Lorentzian
         f = (x, p) -> 1 / ((x - 2)^2 + 1),
         domain = (2, Inf),
         solution = pi / 2
     ),
-    (; # 14. single-variable flipped, shifted, semi-infinite lower limit: Lorentzian
+    (; # 17. single-variable flipped, shifted, semi-infinite lower limit: Lorentzian
         f = (x, p) -> 1 / ((x - 2)^2 + 1),
         domain = (Inf, 2),
         solution = -pi / 2
     ),
-    (; # 15. single-variable flipped, shifted, semi-infinite upper limit: Lorentzian
+    (; # 18. single-variable flipped, shifted, semi-infinite upper limit: Lorentzian
         f = (x, p) -> 1 / ((x - 2)^2 + 1),
         domain = (2, -Inf),
         solution = -pi / 2
     ),
-    (; # 16. single-variable finite limits: constant
-        f = (x, p) -> 1.0,
-        domain = (1, 3),
-        solution = 2
+    (; # 19. single-variable finite limits: quadratic
+        f = (x, p) -> x^2,
+        domain = (1, 4),
+        solution = 21
     ),
-    (; # 17. single-variable flipped, finite limits: constant
-        f = (x, p) -> 1.0,
-        domain = (3, 1),
-        solution = -2
+    (; # 20. single-variable flipped, finite limits: quadratic
+        f = (x, p) -> x^2,
+        domain = (4, 1),
+        solution = -21
     )
 )
 

From 7a279b5dd96aa19729d41853a19c94e74b83868b Mon Sep 17 00:00:00 2001
From: Daniel Wennberg <daniel.wennberg@gmail.com>
Date: Wed, 28 Feb 2024 12:13:45 -0800
Subject: [PATCH 2/2] Fix incorrect finite domain transformation

---
 src/infinity_handling.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl
index 0cf5c54..25cfa69 100644
--- a/src/infinity_handling.jl
+++ b/src/infinity_handling.jl
@@ -33,7 +33,7 @@ function substitute_t(t::Number, lb::Number, ub::Number)
         lb + t * den * u, den^2 * u
     else
         den = (ub - lb) * oftype(t, 0.5)
-        lb - t * den, den
+        lb + (1 + t) * den, den
     end
 end
 function substitute_t(t::AbstractVector, lb::AbstractVector, ub::AbstractVector)