let f m x be a computationally expensive function that takes two arguments and is computed using a numerical algorithm. f is a differentiable function that has the property that the derivate of f m x is equal to f m+1 x.
In order to lower the computational cost of calling this function thousand of times, we can expand it in a 6 terms Taylor series around xk,
f(m,x) = f(m,xk + dx) = f(m,xk) -f(m+1,xk)*dx + 0.5*f(m+2,xk)(dx)^2 - 1/6 * f(m+3,xk)(dx)^3 + 1/24 * f(m+4,xk)(dx)^4 + 1/120 * f(m+5,xk)(dx)^5......
Where f(m+1,xk) is the first derivate of f m.
Using the Taylor series we can generate a grid of xk points [0.1, 0.2, ..xk,.. xmax] in mmax + 6 dimesions, where xmax and mmax are the maximum value of x and m, respectively. Basically, the idea is to populate the nodes of the grid calling the expensive function f and the value of the function between nodes are calculated using the above Taylor series expansion.
This grid can be generated using a lazy map,
import Control.Applicative ((<$>),(<*>))
import Data.Map.Lazy as M
data Fun = Fun Double Double deriving (Ord,Eq) -- data representing the point f m x
generateGrid :: (Double -> Double -> Double)-> Int -> Int -> Double -> Map Fun Double
generateGrid f mMax xmax delta = M.fromList $ zip keys vals
where keys = Fun <$> [fromIntegral x | x <- [0..mMax]] <*> [delta*fromIntegral(i) | i <- [0..xmax]]
vals = fmap (\(Fun m x) -> f m x) keys
funTaylor :: Map Fun Double -> Double -> Double -> Double
funTaylor grid m x = (fun m xk) - (fun (m+1) xk)*dx + 0.5*(fun (m+2) xk)*dx^2 -
(1/6)*(fun (m+3) xk)*dx^3 + (1/24)*(fun (m+4) xk)*dx^4 - (1/120)*(fun (m+5) xk)*dx^5
where delta = 0.1
dx = x-xk
xk = (calcClosestPoint x )* delta
fun = calcutateFun grid
calcClosestPoint r = fromIntegral $ floor $ x /delta
calculateFun :: Map Fun Double -> Double -> Double -> Double
calculateFun = ....
We have a grid where each node is a thunk of the expensive function f. Now every time that we need to calculate f we call funTaylor instead of f.
Notice that's this is a virtual grid in the sense that every time that we call funTaylor m x, only the points involved in the Taylor series are evaluated while the rest of the grid remains unevaluated.
Now suppose that I don't know a priori what are the number xmax and mmax, this virtual space is huge therefore guessing the limits and generating millions of thunks is not a good option.
The idea then is to be able to construct the grid "on the fly", every time that the funTaylor function is invocated new entries are added to the Map. The problem is that the funTaylor function is called from several parts in my code, from several functions that are evaluated in parallel. Then, I need to share the grid between all the functions that call funTaylor.
One alternative is to create a Channel or Queue that would work like a server between the grid and the functions that call funTaylor, every time a new value funTaylor m x is required it's sent to a Queue that update the grid avoiding the problem of sharing the State-grid between several functions.
The question is then: What is the best way to generate the grid "on the fly?"
