Posted by
Luigi Ballabio-4 on
Mar 19, 2002; 1:33pm
URL: http://quantlib.414.s1.nabble.com/LexicographicalView-class-help-tp1903p1904.html
Re: [Quantlib-users] LexicographicalView class
help!!
Hi Toyin,
sorry
for the delay---I pushed your mail in the stack and it was a while
before I could pop it...
At 7:42 PM +0000 3/14/02, Toyin Akin wrote:
Your doc
refers to the possibility of implementing 2D Finite differences using
the Lexicographical
class. It
also states that no 2D operator is currently implemented apart from
the identity operator.
(Page 20 of
the pdf, version 0.3.0a5)
Could you
point me to the location of where this 2D identity operator is
implemented?
Ouch, that was a leftover from code that isn't there anymore.
There's no 2D identity operator.
This would
give me ideas of how the other 2D operators can be implemented with
respect to the
Lexicographical class.
Hmm, I probably wrote that in a misleading way. Let's see if I
can write this straight. Apologies if I'm stating the
obvious.
Representing a discretized 1-D function f(x) is
straightforward---one just chooses a series of x_i, i in 0..N-1, and
the function is then the array
[f_0,...,f_{N-1}] where f_i = f(x_i). 1-D differential operators
are discretized just as easily. Let's take the x_i are evenly spaced
with x_{i+1} - x_i = h.
The discrete derivative in x_i is f'_i = (x_{i+1}-x_{i-1})/2h
which translates in the differential operator
|?
? |
|-1 0
1
|
1 | -1 0
1 |
- |
..... |
2h | -1
0 1 |
|
-1 0 1|
| ? ?|
where the ? depend on which boundary conditions are chosen.
2-D functions and operators are a different breed of cats.
Given a function f(x,y), one way to start is to discretize the x
and y domain as [x_0,...,x_{N-1}] and [y_0,...,y_{M-1}]. We end up
with a collection of points (x_i,y_j) and with the corresponding
f(x_i,y_j). Hic sunt leones. We want to put these points in an array
[f_0,...,f_{N*M-1}] (note: not a matrix!) so that we can use linear
algebra to apply a differential operator D to the function. However,
there is no obvious ordering of the points (or rather, there are more
than one), or put differently, there is no single mapping of the (i,j)
into the indexes of the array f.
Let's say we have [x_0,x_1,x_2] and [y_0,...,y_3]. The
corresponding points are
f_00 f_10 f_20
f_01 f_11 f_21
f_02 f_12 f_22
f_03 f_13 f_23
Lexicographical mappings give the array
[f_00 f_10 f_20 f_01 f_11 f_21 f_02 f_12 f_22 f_03 f_13
f_23];
reverse lexicographical mapping gives
[f_00 f_01 f_02 f_03 f_10 f_11 f_12 f_13 f_20 f_21 f_22
f_23];
red-black mapping gives
[f_00 f_20 f_11 f_02 f_22 f_13 f_10 f_01 f_21 f_12 f_03
f_23]
(consider the points above as on a chessboard; order
lexicographically all the red points first, then the black).
There are many other which have different applications. For
instance, lexicographical orderings seem most natural. However,
red-black ordering has the advantage that in a 5-star differential
formula (that is, a derivative in one point depends on the points
immediately above, below and besides it but not the ones touching its
corners only) all points in the formula except the central one have
the same color. Thus, the resulting differential operator has tighter
bands than the one resulting from lexicographical ordering.
Now, as for your question: I'm afraid I have only half an
answer.
Mapping classes such as LexicographicalView have the purpose of
making it easier to manage discretized 2-D functions. Its use is
straightforward enough: given [x_0,x_{N-1}] and [y_0,y_{M-1}], you can
write:
Array f(N*M); // discretized function
LexicographicalView view(f.begin(),f.end(),N); // put
a "view" over f
for (Size i=0; i<N; i++)
for (Size j=0; j<M; j++)
view[i][j] =
f(x[i],y[j]); // view will put this value at the
// right position in f
had you used RedBlackView (which doesn't exist yet, alas) you
would have written exactly the same for loop; however, the discrete
values would have ended up at different positions in the f array
corresponding to a red-black ordering. Thus, these classes save you
the hassle of calculating such indexes yourself.
As for implementing 2-D differential operators, it is clear that,
say, the first derivative with respect to x will have a different
matricial form depending on the mapping used for storing the function
into its array (I would give some examples, but they would be too
painful for me to write and for you read in plain ASCII).
Therefore, there is some coupling between constructing the array
and applying the operator: they both must be aware of the mapping used
(note that this is not a problem introduced by the LexicographicalView
class; the problem would be the same had we build the array by
hand).
I admit I don't have a design clear in my mind right now.
Of course, one possibility would be to let the user shoot his own
foot (which is actually kind of appealing since it would promote
natural selection :) but I will discard it for the time being. Other
possibilities all involve carrying view information with the
array.
We could pass the view to the operator and let it operate on the
array through the view, but this would prevent composing operators
unless we jump through some major hoops.
Another possibility I like better is to declare the operators as
templates and specialize them for the view classes. This would allow
us to write code such as:
template <class View>
result_type doSomething(whatever) {
Array f(N*M);
View view(f.begin(),f.end(),N);
... // fill the array
Operator L = Dx<View>() +
2.0*Dy<View>() // simple matricial algebra
L.applyTo(f);
... // etc.
}
I will have to think about this and put together a QuEP.
Time to go to bed now.
Bye,
Luigi
--