Using Numeric.Integration.TanhSinh for N-dimensional integration

Question

I am using Numeric.Integration.TanhSinh for numerical integration in Haskell. This defines a function

parTrap :: (Double -> Double) -> Double -> Double -> [Result]

where the first argument is the 1-dimensional function to be integrated, then upper and lower bounds. I have a wrapper which transforms this function as

ttrap f xmin xmax = (ans, err)
   where
      res = absolute 1e-6 $ parTrap f xmin xmax
      ans = result res
      err = errorEstimate res

To integrate a 2-dimensional function, I can use

 ttrap2 f y1 y2 g1 g2 = ttrap h y1 y2 -- f ylower yupper (fn for x lower) (fn for x upper) 
       where
           h y = fst $ ttrap (flip f y) (g1 y) (g2 y)

ttrap2_fixed f y1 y2 x1 x2 = ttrap2 f y1 y2 (const x1) (const x2)

The idea of ttrap2_fixed is that I can now do a double integral where the function is (Double -> Double -> Double), and bounds are y1 y2 x1 x2.

Using this pattern, I can define higher order integration functions

ttrap3_fixed :: (Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap3_fixed f z1 z2 y1 y2 x1 x2 = fst $ ttrap h z1 z2
  where
    h z  = fst $ ttrap2_fixed (f z)  x1 x2 y1 y2


ttrap4_fixed :: (Double -> Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap4_fixed f w1 w2 z1 z2 y1 y2 x1 x2 = fst $ ttrap h w1 w2
  where
    h w  = ttrap3_fixed (f w) z1 z2 x1 x2 y1 y2


ttrap5_fixed :: (Double -> Double -> Double -> Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double
ttrap5_fixed f u1 u2 w1 w2 z1 z2 y1 y2 x1 x2 = fst $ ttrap h u1 u2
  where
    h u  = ttrap4_fixed (f u) w1 w2 z1 z2 x1 x2 y1 y2

However, I would like to integrate a function of type

f :: [Double] -> Double

with the idea being that the dimensionality of the function can vary within the program. Ideally, I would like a function with type

int_listfn :: ([Double] -> Double) -> [(Double, Double)] -> Double

where I can integrate the multidimensional function over a list of tuples of bounds. As part of this, it seems like I need to use something like continuation passing style to construct the integrator function, but this is where I am getting stuck. Thanks in advance.

---

An example of a f that I would like to integrate would be something like

> let f [x,y,z] = x**3.0 * sin(y) + (1.0/z)

Consider bounds

> let bounds = [(1.0,2.0),(2.0,4.0), (0.0,3.0)]
> int_listfn f bounds

This should be equivalent to calculating

---

EDIT: Adding another example function

f1 :: Double -> Double
f1 x = 1.0 * x


fn_maker :: [Double -> Double] -> ([Double] -> Double)
fn_maker inlist = myfn
  where
    myfn xlist = product $ zipWith (\f x -> f x) inlist xlist


m = 4

f_list = fn_maker (replicate f1 m)

f_list has type [Double] -> Double and is equivalent to f x y z w = x * y * z * w. I thought that the type [Double] -> Double would be appropriate because the dimensionality of the function, m, is a parameter. Perhaps I need to change the design.

---

List function to curried functions, @fizruk I have been using this in my code and seems to work, although I need to keep track of which function to call based on the size of the bounds list.

static_1 :: ([Double] -> Double) -> (Double -> Double)
static_1 f = f'
  where
    f' x = f [x]


static_2 :: ([Double] -> Double) -> (Double -> Double -> Double)
static_2 f = f'
  where
    f' x y = f [x,y]


static_3 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double)
static_3 f = f'
  where
    f' x y z = f [x,y,z]


static_4 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double -> Double)
static_4 f = f'
  where
    f' x y z w = f [x,y,z,w]


static_5 :: ([Double] -> Double) -> (Double -> Double -> Double -> Double -> Double -> Double)
static_5 f = f'
  where
    f' x y z w u = f [x,y,z,w,u]


int_listfn :: ([Double] -> Double) -> [(Double, Double)] -> Double
int_listfn f bounds  
  | f_dim == 1 = fst $ ttrap (static_1 f) (fst (bounds !! 0)) (snd (bounds !! 0))
  | f_dim == 2 = fst $ ttrap2_fixed  (static_2 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1))
  | f_dim == 3 = ttrap3_fixed  (static_3 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2))
  | f_dim == 4 = ttrap4_fixed  (static_4 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2)) (fst (bounds !! 3)) (snd (bounds !! 3))
  | f_dim == 5 = ttrap5_fixed  (static_5 f) (fst (bounds !! 0)) (snd (bounds !! 0)) (fst (bounds !! 1)) (snd (bounds !! 1)) (fst (bounds !! 2)) (snd (bounds !! 2)) (fst (bounds !! 3)) (snd (bounds !! 3)) (fst (bounds !! 3)) (snd (bounds !! 3))
  | otherwise = error "Unsupported integral size"
  where
    f_dim = length bounds

Answer 1

The idea

For convenience I introduce these type aliases:

type Res   = (Double, Double)
type Bound = (Double, Double)

Let's have a closer look at the types of ttrap_fixedN:

ttrap :: (Double -> Double) -> Double -> Double -> Res
ttrap_fixed2 :: (Double -> Double -> Double) -> Double -> Double -> Double -> Double -> Res

Obviously, we can pair bounds to get a shorter and cleaner version:

ttrap :: (Double -> Double) -> Bound -> Res
ttrap_fixed2 :: (Double -> Double -> Double) -> Bound -> Bounds -> Res

Furthermore, we can collect Bounds together for N > 1:

ttrap_fixed2 :: (Double -> Double -> Double) -> (Bound, Bound) -> Res
ttrap_fixed3 :: (Double -> Double -> Double -> Double) -> (Bound, Bound, Bound) -> Res

Notice how we made all ttrap_fixedN functions to take precisely 2 arguments. Also notice that the type of the second argument (the arity of tuple with Bounds) depends on the first (the arity of function to integrate).

Now it should be clear that a general ttrap_fixed function would depend on the arity of given function and we need type class to implement that kind of polymorphism. Besided integrate method (this is the general ttrap_fixed) that class would need an associated type synonym for the second argument:

class Integrable r where
  type Bounds r :: *
  integrate :: r -> Bounds r -> Res

Solution sketch

We are going to need following extensions:

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}

import Numeric.Integration.TanhSinh

Here's the only function from the question I'll use:

ttrap :: (Double -> Double) -> Double -> Double -> Res
ttrap f xmin xmax = (ans, err)
   where
      res = absolute 1e-6 $ parTrap f xmin xmax
      ans = result res
      err = errorEstimate res

Define Integrable type class:

class Integrable r where
  type Bounds r :: *
  integrate :: r -> Bounds r -> Res

And it's instances

instance Integrable Double where
  type Bounds Double = ()
  integrate x _ = (x, 0)

instance Integrable r => Integrable (Double -> r) where
  type Bounds (Double -> r) = (Bound, Bounds r)
  integrate f ((xmin, xmax), args) = ttrap g xmin xmax
    where
      g x = fst $ integrate (f x) args

Go test it:

main :: IO ()
main = do
  let f :: Double -> Double -> Double -> Double
      f x y z = x ** 3.0 * sin(y) + (1.0/z)
      zbounds = (1.0, 2.0)
      ybounds = (2.0, 4.0)
      xbounds = (0.0, 3.0)
      res = integrate f (xbounds, (ybounds, (zbounds, ())))
  print res -- prints (8.9681929657648,4.732074732061164e-10)

Currying Bounds

As you might have noticed, we now have not very convenient nested tuples for Bounds. It would be nice if integrate f would return a curried function, such as:

integrate (*) :: Bound -> Bound -> Res

instead of

integrate (*) :: (Bound, (Bound, ())) -> Res

Unfortunately, I failed to find any simple way to refactor Integrable to allow that. However, we can solve this issue with another typeclass hackery. The idea is to introduce a polymorhic curryBounds:

class CurryBounds bs where
  type Curried bs a :: *
  curryBounds :: (bs -> a) -> Curried bs a

The instances are straightforward:

instance CurryBounds () where
  type Curried () a = a
  curryBounds f = f ()

instance CurryBounds bs => CurryBounds (b, bs) where
  type Curried (b, bs) a = b -> Curried bs a
  curryBounds f = \x -> curryBounds (\xs -> f (x, xs))

Now we can define a nicer version of integrate:

integrate' :: (Integrable r, CurryBounds (Bounds r)) => r -> Curried (Bounds r) Res
integrate' = curryBounds . integrate

Example

>>> let f x y z = x**3.0 * sin(y) + (1.0/z) :: Double
>>> integrate' f (0, 3) (2, 4) (1, 2)
(8.9681929657648,4.732074732061164e-10)

Answer 2

Fizruk's answer seems to be correct (although I haven't fact checked the details).... But I thought I should add a complementary set of comments that explain why what you are trying to do won't work.

The function type to integrate has type

[Double]->Double

This doesn't quite do what you want it to do. You want to f to have a preset domain dimensionality, but in the world of Haskell the type signature you gave would actually take in any number of values.... The following would be a valid consistent definition of f

f [] = 1.0
f [x] = x+1
f [x, y] = x+y

If you just passed this value into in integration function, it would have no way to know which line in the definition you meant to specify (you could add this dimensionality in as an extra parameter, but it is kind of overkill to use a mechanism that can take an array of arbitrary length then use a param to fix the length.... Although considering the alternatives, this might be the best solution).

You could use tuples instead. This would solve the dimensionality question, but writing code that works across arbitrary length tuples is hard, and generally needs some GHC extension.

I like Fizruk's approach, it is clean and very sleek, and you get to use fully curried functions as they were intended to be used.... The only downside is that you will have to learn slightly more advanced Haskell to understand how it works (although perhaps that is an advantage for you).

---

@SteveJB pointed out that the dimensionality is given as the length of the bounds param (an obvious fact that I had missed).... But I still think that the array approach has some problems.

First, I would naturally think to use recursion here (ie- to integrate an in N dimensions, sum the N-1 dimensions at spaced slices), but that is really hard to do with an array function and sized bounds. Once you peeled off the outermost bounds, the size would change and the wrong "subfunction" in f would be used. You could solve this by pulling out the size of the bounds on a wrapper function, but the internal function would need a dimensionality param, anyway :). Also, it isn't so easy to "curry" a function of the type [Double]->Double in the way you are using it. Ultimately you might have to pass in all the slice positions.

Again, I want to stress that you could make this work, but it would probably be messier than what Fizruk did.