http://quantlib.414.s1.nabble.com/LexicographicalView-class-help-tp1903p1907.html
----- Original Message -----
Sent: Tuesday, March 19, 2002 9:32
PM
Subject: 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
--