QuoteFDM is an algorithm, there is no 'mathematical correctness' regarding element indices. Just convention and personal preference. Let the choice of prefered convention up to the user.Mathematicians start at index == 1, C at 0, Fortran at 1. (std::vector always at 0). FDM maths also starts at 1.FDM is an application, consisting of loosely-coupled and configurable components. The application has many algorithms in it.When I deliver V1 those interested can port it to generic and I can provide feedback. But it is not on the critical path, it's just doing the same stuff in different ways.

QuoteOriginally posted by: CuchulainnQuoteFDM is an algorithm, there is no 'mathematical correctness' regarding element indices. Just convention and personal preference. Let the choice of prefered convention up to the user.Mathematicians start at index == 1, C at 0, Fortran at 1. (std::vector always at 0). FDM maths also starts at 1.That does not hold in general and it depends on convention, if anything.When writing an FDM solver for u(t,x) discretized as v(n,m), I was very explicitly treating t=0 as n=0 (the intuitive choice). The problem was interoping with software that couldn't do that (and the resulting mismatch between the notation and the given software -- not elegant). It's quite the same for a barrier at stock price x=0 (and thus m=0). If you have IC or BCs that include 0 in your domain, then it's the most natural choice.If anything, the situation gets different if you want domains with negative numbers (like solving a heat equation on x in [-1, +1]) or complex domains -- but then neither choice helps, IMHO. Then, one has to use an equivalent of std::begin and std::end for each axis -- and the internal representation/indexing of a container doesn't matter anyway.

QuoteOriginally posted by: PolterQuoteOriginally posted by: CuchulainnQuote [@outrun] FDM is an algorithm, there is no 'mathematical correctness' regarding element indices. Just convention and personal preference. Let the choice of prefered convention up to the user.Mathematicians start at index == 1, C at 0, Fortran at 1. (std::vector always at 0). FDM maths also starts at 1.That does not hold in general and it depends on convention, if anything.When writing an FDM solver for u(t,x) discretized as v(n,m), I was very explicitly treating t=0 as n=0 (the intuitive choice). The problem was interoping with software that couldn't do that (and the resulting mismatch between the notation and the given software -- not elegant). It's quite the same for a barrier at stock price x=0 (and thus m=0). If you have IC or BCs that include 0 in your domain, then it's the most natural choice.If anything, the situation gets different if you want domains with negative numbers (like solving a heat equation on x in [-1, +1]) or complex domains -- but then neither choice helps, IMHO. Then, one has to use an equivalent of std::begin and std::end for each axis -- and the internal representation/indexing of a container doesn't matter anyway.Can you give a simple 101 example of what you are describing, especially how begin() pans out in a fd scheme. It's not clear what the advantages would be. What if I want 'ghost points'??I have a Range class in both x and t which then get discretised. So, I distinguish between the PDE domain and the discrete FD domain. No problem. Most numerical analysis texts index from 1, at least those from the 'golden era'. For the V2 feature list. In principle, begin() and MinIndex() are 2 sides of the same coin, the latter being a 'direct' index.

QuoteOriginally posted by: CuchulainnQuoteOriginally posted by: PolterQuoteOriginally posted by: CuchulainnQuoteFDM is an algorithm, there is no 'mathematical correctness' regarding element indices. Just convention and personal preference. Let the choice of prefered convention up to the user.Mathematicians start at index == 1, C at 0, Fortran at 1. (std::vector always at 0). FDM maths also starts at 1.That does not hold in general and it depends on convention, if anything.When writing an FDM solver for u(t,x) discretized as v(n,m), I was very explicitly treating t=0 as n=0 (the intuitive choice). The problem was interoping with software that couldn't do that (and the resulting mismatch between the notation and the given software -- not elegant). It's quite the same for a barrier at stock price x=0 (and thus m=0). If you have IC or BCs that include 0 in your domain, then it's the most natural choice.If anything, the situation gets different if you want domains with negative numbers (like solving a heat equation on x in [-1, +1]) or complex domains -- but then neither choice helps, IMHO. Then, one has to use an equivalent of std::begin and std::end for each axis -- and the internal representation/indexing of a container doesn't matter anyway.Can you give a simple 101 example of what you are describing, especially how begin() pans out in a fd scheme. It's not clear what the advantages would be. What if I want 'ghost points'??I have a Range class in both x and t which then get discretised. So, I distinguish between the PDE domain and the discrete FD domain. No problem. Suppose that the notation is as above -- u(t,x) is the solution function, v(n,m) is the discretization.Let's start with a situation where we fix the initial indexes at 0 (for you it'd be 1, doesn't matter here).Suppose the user has specified the domain of "u" by delimiters, for simplicity assume rectangular domain -- t_begin, t_end, x_begin, x_end.At some point, when filling the values given by IC and BCs we want to write at those subsets of the domain.Here, for simplicity assume we have a way to access the values function "u" takes there (IC and BCs are given) with an interface u(t_begin, x), u(t, x_begin), u(t, x_end) (an attempt to access an interior point here not given by IC or BCs is not applicable, might as well assume the behavior is undefined in such a case).At some point, our code will contain something like this:for all interior m: v(0,m) = u(t_begin, x); // ICfor all interior n: v(n,0) = u(t, x_begin); v(n,m_steps) = u(t, x_end); // BCsNow, t_begin in "u" might have been even negative, but it doesn't matter, since we're writing to container "v" anyway.Further, we can abstract from the initial indexes by introducing n_begin, n_end, m_begin, m_end (m_steps above now gives the # of steps from m_begin to m_end):for all interior m: v(n_begin,m) = u(t_begin, x); // ICfor all interior n: v(n,m_begin) = u(t, x_begin); v(n,m_end) = u(t, x_end); // BCsNow, one can probably notice that we have notions of begin/end for each axis -- both (t,x) and (n,m).As you have noticed, you introduce a notion of the range at this point.Now, as a concept, range supports a begin/end access interface to provide delimiting values (modulo dereferencing).So, if you have a discretized range for "m" ((here: first and only, but it doesn't matter) spatial axis of the discrete FD domain) being a discretization of the range for "x" ((here: first and only, but it doesn't matter) spatial axis of the PDE domain), you can use begin(m), end(x), etc... instead of the above.When I was writing my code, I couldn't do that in a standard way (no C++11 then, member functions won't work with everything), but now that we have a generic standard notion of those delimiters, perhaps we can support this (providing suitable overloads for non-C++11 compilers, of course) -- makes sense?

QuoteOriginally posted by: CuchulainnIn principle, begin() and MinIndex() are 2 sides of the same coin, the latter being a 'direct' index. Yes!, modulo dereferencing it's the same thing! The advantage of begin/end is that it's a well-known convention, now also standardized.

Just a remark on notation; the payoff function is not a boundary condition but a terminal (t = T) OR initial condition (t = 0). QuoteF(S,T,K) { return max(S-K,0); }One choice is F(y,0, K) where we march from t = 0 and y = S/(S+a). I think just wait until I write my document and then we can discuss new features.

QuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnJust a remark on notation; the payoff function is not a boundary condition but a terminal (t = T) OR initial condition (t = 0). QuoteF(S,T,K) { return max(S-K,0); }One choice is F(y,0, K) where we march from t = 0 and y = S/(S+a). I think just wait until I write my document and then we can discuss new features.Nice! That could indeed be modelled as a nonuniform grid. There will be many more issues (more complicated cashflows, American, dividends).You're right that we will wait for V1.Yes. There are so many dimensions / issues soon that my head it simply swirls

QuoteYes! One approach for the detailed design is it analyze each concept we use, and try to reduce it to irreducable things that have value in a broader sense than this application (range, grid, value-function, cashlows,..) I expect well have many generic math operators that have a broad value for a wide range of application domains -like the basic set of functions in Matlab- and also commonly known concepts. We'll just have to iterate, and try to come up with something that's scales as ISS gets biggerThis is an approach but I don't think it will scale well to more complex PDE/FDM environments. They are just too complex to handle using C++ alone. A high level system approach is needed and that is what I have done. Is there a recipe for discovering concepts;is it by serendipity (luck, trial and error, iteratve, bottom-up) or by a defined process? And how do we know if the concepts are accurate?In my design I know the design and implementation of 1,2 and 3 factor model reasonably well and whether to use OOP vs GP is for later. I can choose between different implementation options.Once I know the problem domain, life becomes easier, as is always been in application development.I asked this already; let's say I have the PDEs {BS, UVM, FPE}, Factors {1, 2,3 }, fd schemes {CN, ADE, ADI, Soviet Splitting} and mesh {uniform, adaptive} and Domain {truncate, transformation} . How do I build this N:M system? I don'think it can be done and kept understandable at the same time without some standard design techniques. QuoteThese concepts will probably naturally arise when we try to re-use transforms like s/(s+a) on binomial and trinomial trees But more to 2 and 3 factor PDE models.

QuoteOriginally posted by: CuchulainnIs there a recipe for discovering concepts;is it by serendipity (luck, trial and error, iterative, bottom-up) or by a defined process? And how do we know if the concepts are accurate?Realistically, I think it's a combination of both.Theoretically, concept-discovery is very close to the axiom-system-determination in reverse mathematics in mathematical logic: http://en.wikipedia.org/wiki/Reverse_Mathematics.Might be a source of nice analogies (at the end of this post):QuoteIn reverse mathematics, one starts with a framework language and a base theory--a core axiom system--that is too weak to prove most of the theorems one might be interested in, but still powerful enough to develop the definitions necessary to state these theorems. For example, to study the theorem "Every bounded sequence of real numbers has a supremum" it is necessary to use a base system which can speak of real numbers and sequences of real numbers.For each theorem that can be stated in the base system but is not provable in the base system, the goal is to determine the particular axiom system (stronger than the base system) that is necessary to prove that theorem. To show that a system S is required to prove a theorem T, two proofs are required. The first proof shows T is provable from S; this is an ordinary mathematical proof along with a justification that it can be carried out in the system S. The second proof, known as a reversal, shows that T itself implies S; this proof is carried out in the base system. The reversal establishes that no axiom system S' that extends the base system can be weaker than S while still proving T."Generic programming" as understood in the sense meant by Stepanov:QuoteAn algorithm that is defined on all the types in a concept is called a generic algorithm.[...]It is my belief that behind every useful line of code hides a generic algorithm. This belief is based on a personal experience: when I see a piece of code I immediately start asking, "What is the underlying concept that makes this code work?"[...]Now, let us spend some time deriving generic algorithms. The approach is quite simple. First we need to find a useful piece of code. We can use all kinds of sources: existing libraries, Knuth, application code. It should, however, be a useful piece. There should be some evidence that people use it or want to use it. Secondly, we see what makes it work and try to abstract the requirements and identify the concept on which it is really defined."Elements of Programming" by Alexander Stepanov and Paul McJones is the standard reference. I'd recommend the lecture notes (under "Class notes" on http://www.stepanovpapers.com/) as a background reading (the focus is more on the practical C++ part ("how to implement an abstract algorithm?"), while the book's focus is on the theoretical foundations ("what is an abstract algorithm?"). Recall also Alexander Stepanov: STL and Its Design Principleshttp://www.stepanovpapers.com/stepanov-abstrac ... tification and Organization1. Find algorithms and data structures2. Implement them3. Create usable taxonomyQuoteGeneric Programming 1. Take a piece of code2. Write specifications3. Replace actual types with formal types4. Derive requirements for the formal types that imply these specificationsCompare with reverse mathematics above, think concepts and requirements where you see "axiom systems" S, S'..."A base theory--a core axiom system", B is analogous to an initial "piece of code" P with some underlying base concept C_B, not necessarily good enough to make P work correctly in general.Being able "to prove a theorem T" is analogous to making "a piece of code P working."The answer to "the necessary axiom system to prove a theorem", S is analogous to the answer to "what is the underlying concept that makes this code work?", C.Initially, it might be non-trivial to start from anything else than P and C_B.And even then chances are it will take some serious brainstorming and various C', C'', C''', ... until we arrive at desired C.

QuoteOriginally posted by: PolterQuoteOriginally posted by: CuchulainnIs there a recipe for discovering concepts;is it by serendipity (luck, trial and error, iterative, bottom-up) or by a defined process? And how do we know if the concepts are accurate?Realistically, I think it's a combination of both.Theoretically, concept-discovery is very close to the axiom-system-determination in reverse mathematics in mathematical logic: http://en.wikipedia.org/wiki/Reverse_Mathematics.Might be a source of nice analogies (at the end of this post):QuoteIn reverse mathematics, one starts with a framework language and a base theory--a core axiom system--that is too weak to prove most of the theorems one might be interested in, but still powerful enough to develop the definitions necessary to state these theorems. For example, to study the theorem "Every bounded sequence of real numbers has a supremum" it is necessary to use a base system which can speak of real numbers and sequences of real numbers.For each theorem that can be stated in the base system but is not provable in the base system, the goal is to determine the particular axiom system (stronger than the base system) that is necessary to prove that theorem. To show that a system S is required to prove a theorem T, two proofs are required. The first proof shows T is provable from S; this is an ordinary mathematical proof along with a justification that it can be carried out in the system S. The second proof, known as a reversal, shows that T itself implies S; this proof is carried out in the base system. The reversal establishes that no axiom system S' that extends the base system can be weaker than S while still proving T."Generic programming" as understood in the sense meant by Stepanov:QuoteAn algorithm that is defined on all the types in a concept is called a generic algorithm.[...]It is my belief that behind every useful line of code hides a generic algorithm. This belief is based on a personal experience: when I see a piece of code I immediately start asking, "What is the underlying concept that makes this code work?"[...]Now, let us spend some time deriving generic algorithms. The approach is quite simple. First we need to find a useful piece of code. We can use all kinds of sources: existing libraries, Knuth, application code. It should, however, be a useful piece. There should be some evidence that people use it or want to use it. Secondly, we see what makes it work and try to abstract the requirements and identify the concept on which it is really defined."Elements of Programming" by Alexander Stepanov and Paul McJones is the standard reference. I'd recommend the lecture notes (under "Class notes" on http://www.stepanovpapers.com/) as a background reading (the focus is more on the practical C++ part ("how to implement an abstract algorithm?"), while the book's focus is on the theoretical foundations ("what is an abstract algorithm?"). Recall also Alexander Stepanov: STL and Its Design Principleshttp://www.stepanovpapers.com/stepanov-abstrac ... tification and Organization1. Find algorithms and data structures2. Implement them3. Create usable taxonomyQuoteGeneric Programming 1. Take a piece of code2. Write specifications3. Replace actual types with formal types4. Derive requirements for the formal types that imply these specificationsCompare with reverse mathematics above, think concepts and requirements where you see "axiom systems" S, S'..."A base theory--a core axiom system", B is analogous to an initial "piece of code" P with some underlying base concept C_B, not necessarily good enough to make P work correctly in general.Being able "to prove a theorem T" is analogous to making "a piece of code P working."The answer to "the necessary axiom system to prove a theorem", S is analogous to the answer to "what is the underlying concept that makes this code work?", C.Initially, it might be non-trivial to start from anything else than P and C_B.And even then chances are it will take some serious brainstorming and various C', C'', C''', ... until we arrive at desired C.That's very clear and excellent perspective. It's the way the humn brain chunks and classifies information and is not specific to mathematics. Software industry as well.Two approaches1. Bottom-up, concept discovery and 'inctemental programing' 2. Top-down decomposition and concepts known already (because there is already an implementation)Both appraches should converge to the same result in C++.In the case of 2. some patterns already exist and OO implementations, which we can formulate as concepts. Outrun's fdm function argument list is essentially the same.Known of the problem's requirements is a pre-requisite for choices 1 and 2.