[math] Greetings from a newcomer (but not to numerics)

[Al Chou]

Greetings to the commons-math community! I emailed Robert Donkin privately a
few days ago asking, as I take it others have, about the reasoning behind
creating a numerics library from scratch rather than incorporating an existing
one. I think I understand that reasoning now, but despite not wanting to
dampen anyone's enthusiasm, I do want to caution those who have not done a lot
of numerical computing that it can very easily be done wrong. A big reason why
there's a lot of legacy code in scientific computing is that it's hard to get
numerical algorithms right, so once you do, there's great inertia towards
changing anything (there are of course other big reasons, such as the fact that
many computational scientists are not actually as computer-savvy as they are
science-savvy, so they're not very facile at creating new code).

As an example, the high school quadratic formula
r = ( -b +- sqrt( b^2 - 4 * a * c ) ) / ( 2 * a )
is extremely inaccurate, because of finite-precision arithmetic, when
4 * a * c << b^2 ,
because of the subtraction of nearly-equal values. But in high school they
never teach you that fact.

On a more positive note, let me recommend _Numerical Recipes_, _Numerical
Methods that [Usually] Work_ (which incidentally presents an alternate
quadratic formula for use in the regime where the traditional one fails), and
_Real Computing Made REAL_ as easy-to-read and down-to-earth references that
show how easy it is for the naive implementation to be woefully inadequate as
well as teach how to do numerical math correctly.

I'd like to participate in commons-math in an advisory capacity of some sort,
as I can't in good faith commit to being able to contribute code to the
project. Robert indicated that such a role would be useful to the project, so
I hope you all feel the same!


Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Phil Steitz]

Al Chou wrote:
> Greetings to the commons-math community! I emailed Robert Donkin privately a
> few days ago asking, as I take it others have, about the reasoning behind
> creating a numerics library from scratch rather than incorporating an existing
> one. I think I understand that reasoning now, but despite not wanting to
> dampen anyone's enthusiasm, I do want to caution those who have not done a lot
> of numerical computing that it can very easily be done wrong.

No question about that. As I have stated a couple of times, however, I
do not personally see commons-math as a numerics library. There will
certainly be numerical considerations to worry about, though, and your
point is well taken.

A big reason why
> there's a lot of legacy code in scientific computing is that it's hard to get
> numerical algorithms right, so once you do, there's great inertia towards
> changing anything (there are of course other big reasons, such as the fact that
> many computational scientists are not actually as computer-savvy as they are
> science-savvy, so they're not very facile at creating new code).
>
Whence the wonderful proliferation of Fortran code in Java ;-)

> As an example, the high school quadratic formula
> r = ( -b +- sqrt( b^2 - 4 * a * c ) ) / ( 2 * a )
> is extremely inaccurate, because of finite-precision arithmetic, when
> 4 * a * c << b^2 ,
> because of the subtraction of nearly-equal values. But in high school they
> never teach you that fact.
>
> On a more positive note, let me recommend _Numerical Recipes_, _Numerical
> Methods that [Usually] Work_ (which incidentally presents an alternate
> quadratic formula for use in the regime where the traditional one fails), and
> _Real Computing Made REAL_ as easy-to-read and down-to-earth references that
> show how easy it is for the naive implementation to be woefully inadequate as
> well as teach how to do numerical math correctly.

No question that the first of these at least and Knuth are classics to
refer to. I am not familiar with the second two -- thanks for the tip. I
also refer to Burden and Faires and Atkinson's _Numerical Analysis_
texts. Do you know of any decent web numerical analysis references? I
have not been able to find much. It would be nice to have some of these
as well so we could refer to them in the docs and discussion. In any
case, it is probably a good idea to add a bibliography section to the
web site.

>
> I'd like to participate in commons-math in an advisory capacity of some sort,
> as I can't in good faith commit to being able to contribute code to the
> project. Robert indicated that such a role would be useful to the project, so
> I hope you all feel the same!

Please, please, please! Whatever time you have to review
implementations, comment on strategies, and patiently point out our
numerical blunders will be greatly appreciated.

Your frank assessment of project direction and the whole question of
whether or not we should implement the things on the task list
(http://jakarta.apache.org/commons/sandbox/math/tasks.html) would also
be appreciated.

Phil
>
> Al
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Search - Faster. Easier. Bingo.
> http://search.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[J.Pietschmann]

Phil Steitz wrote:
> Do you know of any decent web numerical analysis references?

The only web reference on mathematics which is advanced enough
and reasonably usable (not a bunch of grad papers) seems to be
wolfram research. It's not too heavy on numerics though.

If you have questions, just ask (me: math MS, 5y in CAD/CAM/CNC,
ODE/PDE stuff up to 1997).

J.Pietschmann

[Phil Steitz]

J.Pietschmann wrote:
> Phil Steitz wrote:
>
>> Do you know of any decent web numerical analysis references?
>
>
> The only web reference on mathematics which is advanced enough
> and reasonably usable (not a bunch of grad papers) seems to be
> wolfram research. It's not too heavy on numerics though.

Agreed. Too bad.

>
> If you have questions, just ask (me: math MS, 5y in CAD/CAM/CNC,
> ODE/PDE stuff up to 1997).

Will do, but we are going to have to have some way to refer to standard
references. Guess we will have to use books :-(

One question. Probably the most numerically challenging (and
duplicative of existing stuff) thing on the task list is RealMatrix
solve() and inverse(). I was planning to use Gaussian elimination with
scaled column pivoting (cf. Burden and Faires, _Numerical Analysis_ alg.
6.3) to implement solve() directly and use this (inefficiently) to
implemeent inverse(). Is this reasonable? If not, can you refer me to a
well-documented algorithm elsewhere that is at least as good as this
numerically and maybe more performant? If the answer is
LU-decomposition using the Crout-Doolittle algorithm (what Jama
http://math.nist.gov/javanumerics/jama/doc/ and some other
implementations do) can you explain why this is better?

Phil

>
> J.Pietschmann
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Mark R. Diggory]

Phil Steitz wrote:
> J.Pietschmann wrote:
>
>> Phil Steitz wrote:
>>
>>> Do you know of any decent web numerical analysis references?
>>
>>
>>
>> The only web reference on mathematics which is advanced enough
>> and reasonably usable (not a bunch of grad papers) seems to be
>> wolfram research. It's not too heavy on numerics though.
>
>
> Agreed. Too bad.
>

Yea, they wouldn't want to share those trade secrets... ;-)

>>
>> If you have questions, just ask (me: math MS, 5y in CAD/CAM/CNC,
>> ODE/PDE stuff up to 1997).
>
>
> Will do, but we are going to have to have some way to refer to standard
> references. Guess we will have to use books :-(
>

I found this resource:
Numerical Recipies online (C, Fortran)
http://www.nr.com/nronline_switcher.html



> One question. Probably the most numerically challenging (and
> duplicative of existing stuff) thing on the task list is RealMatrix
> solve() and inverse(). I was planning to use Gaussian elimination with
> scaled column pivoting (cf. Burden and Faires, _Numerical Analysis_ alg.
> 6.3) to implement solve() directly and use this (inefficiently) to
> implemeent inverse(). Is this reasonable? If not, can you refer me to a
> well-documented algorithm elsewhere that is at least as good as this
> numerically and maybe more performant? If the answer is
> LU-decomposition using the Crout-Doolittle algorithm (what Jama
> http://math.nist.gov/javanumerics/jama/doc/ and some other
> implementations do) can you explain why this is better?
>
> Phil
>
>>
>> J.Pietschmann
>>
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
>> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>>
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Brent Worden]

> -----Original Message-----
> From: Mark R. Diggory [mailto:***@latte.harvard.edu]
> Sent: Sunday, May 25, 2003 12:05 PM
>
> I found this resource:
> Numerical Recipies online (C, Fortran)
> http://www.nr.com/nronline_switcher.html
>

Other online resources that I've found useful are:

Guide to Available Mathematical Software
http://gams.nist.gov/

StatLib
http://lib.stat.cmu.edu/

NetLib
http://www.netlib.org/

The WWW Virtual Library: Mathematics
http://web.math.fsu.edu/Science/math.html

The WWW Virtual Library: RANDOM NUMBERS and MONTE CARLO METHODS
http://random.mat.sbg.ac.at/links/

The WWW Virtual Library: Statistics
http://www.stat.ufl.edu/vlib/statistics.html

Brent Worden
http://www.brent.worden.org

[J.Pietschmann]

Phil Steitz wrote:
> One question. Probably the most numerically challenging (and
> duplicative of existing stuff) thing on the task list is RealMatrix
> solve() and inverse(). I was planning to use Gaussian elimination with
> scaled column pivoting (cf. Burden and Faires, _Numerical Analysis_ alg.
> 6.3) to implement solve() directly and use this (inefficiently) to
> implemeent inverse(). Is this reasonable?

Uh. There is no one-size-fits-all in solving linear equations,
otherwise there wouldn't be so much papers about this.

> If the answer is
> LU-decomposition using the Crout-Doolittle algorithm (what Jama
> http://math.nist.gov/javanumerics/jama/doc/ and some other
> implementations do) can you explain why this is better?

Ouch again. Of course if the matrix is well behaved the answer
is LU decomposition, although I never met Crout-Doolittle before,
so I can't explain it. All I used myself and read about were
simple derivations from the original Gauss algorithm itself. The
problem is to inexpensively detect whether the matrix is well
behaved.
If the matrix turns out to be nearly singular, use Householder or
singular value decomposition, which are both more expensive, both
because of those pesky transcendent functions and harder to
localize memory access.
Also, all methods mentioned so far are for dense, small matrices
which can be kept in memory and wont cause too much cache
thrashing. I'm not sure whether cache thrashing is an issue
with Java programs though. Sparse matrices require usually
special treatment (look up "fill-in").

HTH a bit.

J.Pietschmann

[Brent Worden]

> Phil Steitz wrote:
> > One question. Probably the most numerically challenging (and
> > duplicative of existing stuff) thing on the task list is RealMatrix
> > solve() and inverse(). I was planning to use Gaussian elimination with
> > scaled column pivoting (cf. Burden and Faires, _Numerical
> Analysis_ alg.
> > 6.3) to implement solve() directly and use this (inefficiently) to
> > implemeent inverse(). Is this reasonable?
>
> Uh. There is no one-size-fits-all in solving linear equations,
> otherwise there wouldn't be so much papers about this.
>
> > If the answer is
> > LU-decomposition using the Crout-Doolittle algorithm (what Jama
> > http://math.nist.gov/javanumerics/jama/doc/ and some other
> > implementations do) can you explain why this is better?
>
> Ouch again. Of course if the matrix is well behaved the answer
> is LU decomposition, although I never met Crout-Doolittle before,
> so I can't explain it. All I used myself and read about were
> simple derivations from the original Gauss algorithm itself. The
> problem is to inexpensively detect whether the matrix is well
> behaved.
> If the matrix turns out to be nearly singular, use Householder or
> singular value decomposition, which are both more expensive, both
> because of those pesky transcendent functions and harder to
> localize memory access.
> Also, all methods mentioned so far are for dense, small matrices
> which can be kept in memory and wont cause too much cache
> thrashing. I'm not sure whether cache thrashing is an issue
> with Java programs though. Sparse matrices require usually
> special treatment (look up "fill-in").
>
> HTH a bit.
>
> J.Pietschmann

It should be noted for least square estimates (the primary reason we're
looking for a linear equation solution is needed), the matrix that needs to
be solved/inverted is symmetric and positive-definite. Correct me if I'm
wrong, but I think Cholesky decomposition is a good candidate for this class
of matrices.

Brent Worden
http://www.brent.worden.org

[J.Pietschmann]

Brent Worden wrote:
> It should be noted for least square estimates (the primary reason we're
> looking for a linear equation solution is needed), the matrix that needs to
> be solved/inverted is symmetric and positive-definite. Correct me if I'm
> wrong, but I think Cholesky decomposition is a good candidate for this class
> of matrices.
>
Yes, definitely ok.

J.Pietschmann

[Phil Steitz]

J.Pietschmann wrote:
> Brent Worden wrote:
>
>> It should be noted for least square estimates (the primary reason we're
>> looking for a linear equation solution is needed), the matrix that
>> needs to
>> be solved/inverted is symmetric and positive-definite. Correct me if I'm
>> wrong, but I think Cholesky decomposition is a good candidate for this
>> class
>> of matrices.
>>
> Yes, definitely ok.
>

I would prefer to implement something that is reasonable for a large
class of matrices (not just pos definite). The RealMatrixImpl class
requires in-memory storage, so we are not going to be working on *large*
matrices in that implementation.

Do any of you have implementations that you would like to contribute?

If not, I will imitate the Jama implementation, using the Doolittle
method to do LU-decomposition (described, for example here:
http://www.ecs.fullerton.edu/~mathews/numerical/linear/dol/dol.html),
unless there are strong objections to that.

One more thing that Al pointed out is that inverse() is really a lot
less important than solve(). Even for the regression, the problem is
really to solve the normal equations.

Phil

> J.Pietschmann
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Al Chou]

--- Phil Steitz <***@steitz.com> wrote:
> J.Pietschmann wrote:
> > Brent Worden wrote:
> >
> >> It should be noted for least square estimates (the primary reason we're
> >> looking for a linear equation solution is needed), the matrix that
> >> needs to
> >> be solved/inverted is symmetric and positive-definite. Correct me if I'm
> >> wrong, but I think Cholesky decomposition is a good candidate for this
> >> class
> >> of matrices.
> >>
> > Yes, definitely ok.
> >
>
> I would prefer to implement something that is reasonable for a large
> class of matrices (not just pos definite). The RealMatrixImpl class
> requires in-memory storage, so we are not going to be working on *large*
> matrices in that implementation.
>
> Do any of you have implementations that you would like to contribute?
>
> If not, I will imitate the Jama implementation, using the Doolittle
> method to do LU-decomposition (described, for example here:
> http://www.ecs.fullerton.edu/~mathews/numerical/linear/dol/dol.html),
> unless there are strong objections to that.

I would agree that something fundamental like LU decomposition plus backsolve
(I think that's what Crout and/or Doolittle describe) should be part of the
package first. But we need to keep in mind for the near future that, as J.
Pietschmann points out, there is no single algorithm that's best for all linear
algebra problems. One of the hardest things about numerical math is
categorizing the problem at hand and choosing the best algorithm you can find.


> One more thing that Al pointed out is that inverse() is really a lot
> less important than solve(). Even for the regression, the problem is
> really to solve the normal equations.
>
> Phil

Yes, I would definitely concentrate on solving the system of equations and
leave matrix inversion for later.



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Brent Worden]

> I would prefer to implement something that is reasonable for a large
> class of matrices (not just pos definite). The RealMatrixImpl class
> requires in-memory storage, so we are not going to be working on *large*
> matrices in that implementation.
I have no qualms about that at all.

> Do any of you have implementations that you would like to contribute?

Here is a Cholesky Decomposition class usable by commons-math which can be
used or modified as you see fit:

package org.apache.commons.math;

/**
* Solve linear systems of equations using Cholesky decomposition.
*
* See
http://ikpe1101.ikp.kfa-juelich.de/briefbook_data_analysis/node33.html
* for further information.
*
* @author Brent Worden
*/
public class CholeskyDecomposition {
/**
* Default constructor. Does nothing as this is a utility class.
*/
public CholeskyDecomposition(){
super();
}

/**
* Decompose a symmetric, positive definite matrix into using Cholesky
* factorization. Using this technique, mat is decomposed into LL',
where
* L is a lower triangular matrix.
* @param mat the matrix to be decomposed.
* @return the lower triangular matrix of the Cholesky factorization.
* @throws IllegalArgumentException if mat is not square or is singular.
*/
public RealMatrix decompose(RealMatrix mat){
if((mat == null) ||
(mat.getColumnDimension() <= 0) ||
(mat.getColumnDimension() != mat.getRowDimension()))
{
throw new IllegalArgumentException("matrix must be square.");
}

int n = mat.getColumnDimension(); // size of square matrix
double[][] a = mat.getData(); // original matrix

// decomposition
RealMatrix decomposition = new RealMatrixImpl(n, n);
double[][] data = decomposition.getData();

data[0][0] = Math.sqrt(a[0][0]);
for (int i = 1; i < n; i++) {
data[i][0] = a[i][0] / data[0][0];
}
for (int i = 0; i < n; i++) {
double sum = a[i][i];
for (int j = 0; j < i; j++) {
sum = sum - (data[i][j] * data[i][j]);
}
if(sum <= 0.0){
// TODO: should this be a specialized exception?
throw new IllegalArgumentException("can not decompose
singular matrix.");
}
data[i][i] = Math.sqrt(sum);
for (int j = n - 1; j > i; j--) {
sum = a[j][i];
for (int k = 0; k < i; k++) {
sum = sum - (data[j][k] * data[i][k]);
}
data[j][i] = sum / data[i][i];
}
}

return decomposition;
}

/**
* Take the output of {@link
CholeskyDecomposition#decompose(RealMatrix)}
* and solves a linear set of equations, LL'x = b.
* @param mat LL' matrix generated from Cholesky factorization.
* @param b ??? what is the common name for this vector ???
* @return the solution vector, x
* @throws IllegalArgumentException if mat is not square or is not the
same
* size as b.
*/
public double[] solveUsingDecomposition(RealMatrix mat, double[] b){
if((mat == null) ||
(mat.getColumnDimension() <= 0) ||
(mat.getColumnDimension() != mat.getRowDimension()) ||
(b == null) ||
(b.length != mat.getColumnDimension()))
{
throw new IllegalArgumentException(
"matrix must be square and the same dimension as b
vector.");
}

int n = b.length;
double[] x = new double[n];
double[][] data = mat.getData();

for (int i = 0; i < n; i++) {
double sum = b[i];
for (int j = i - 1; j >= 0; j--) {
sum = sum - (data[i][j] * x[j]);
}
x[i] = sum / data[i][i];
}
for (int i = n - 1; i >= 0; i--) {
double sum = x[i];
for (int j = i + 1; j < n; j++) {
sum = sum - (data[j][i] * x[j]);
}
x[i] = sum / data[i][i];
}

return x;
}

/**
* Solves a linear set of equations, Ax = b. A must be of a symmetric,
* positive definite matrix.
* @param a the symmetric, positive definite matrix
* @param b ??? what is the common name for this vector ???
* @return the solution vector, x
*/
public double[] solve(RealMatrix a, double[] b){
RealMatrix decomposition = decompose(a);
return solveUsingDecomposition(decomposition, b);
}
}

And here is a simple unit test:

package org.apache.commons.math;

import junit.framework.TestCase;

/**
* @author Brent Worden
*/
public class CholeskyDecompositionTest extends TestCase {

/**
* Constructor for CholeskyDecompositionTest.
* @param name
*/
public CholeskyDecompositionTest(String name) {
super(name);
}

public void testSolve(){
double[][] data = {
{ 1, -1, 1},
{-1, 10, -10},
{ 1, -10, 14}
};
double[] b = {2, -2, 6};
double[] expected = {2, 1, 1};

CholeskyDecomposition cd = new CholeskyDecomposition();
double[] x = cd.solve(new RealMatrixImpl(data), b);

for (int i = 0; i < x.length; i++) {
assertEquals(expected[i], x[i], 1e-2);
}
}
}

Brent Worden
http://www.brent.worden.org

[Al Chou]

--- Phil Steitz <***@steitz.com> wrote:
> Al Chou wrote:
> > Greetings to the commons-math community! I emailed Robert Donkin privately
> a
> > few days ago asking, as I take it others have, about the reasoning behind
> > creating a numerics library from scratch rather than incorporating an
> existing
> > one. I think I understand that reasoning now, but despite not wanting to
> > dampen anyone's enthusiasm, I do want to caution those who have not done a
> lot
> > of numerical computing that it can very easily be done wrong.
>
> No question about that. As I have stated a couple of times, however, I
> do not personally see commons-math as a numerics library. There will
> certainly be numerical considerations to worry about, though, and your
> point is well taken.

I'll have to catch up on your posts, Phil, as I don't immediately understand
how a math library can be not a numerics library.


> A big reason why
> > there's a lot of legacy code in scientific computing is that it's hard to
> get
> > numerical algorithms right, so once you do, there's great inertia towards
> > changing anything (there are of course other big reasons, such as the fact
> that
> > many computational scientists are not actually as computer-savvy as they
> are
> > science-savvy, so they're not very facile at creating new code).
>
> Whence the wonderful proliferation of Fortran code in Java ;-)

Remind me sometime to tell you the story my former supervisor told me about an
"expert" he once worked with whose personal conception of OOP was completely
orthogonal to the facilities provided for it in Java, which was the
implementation language they were working in.


> > On a more positive note, let me recommend _Numerical Recipes_, _Numerical
> > Methods that [Usually] Work_ (which incidentally presents an alternate
> > quadratic formula for use in the regime where the traditional one fails),
> and
> > _Real Computing Made REAL_ as easy-to-read and down-to-earth references
> that
> > show how easy it is for the naive implementation to be woefully inadequate
> as
> > well as teach how to do numerical math correctly.
>
> No question that the first of these at least and Knuth are classics to
> refer to. I am not familiar with the second two -- thanks for the tip. I

In the preface to the first editions of _Numerical Recipes_, the authors
acknowledge the influence of Forman Acton's _Numerical Methods that [Usually]
Work_ on their presentation style. I discovered the first edition of it by
accident when I was looking in the library for a copy of _NR_, which they were
all out of at the time, much to my benefit. It was reissued in 1990, and he
followed up with _Real Computing..._, essentially a practitioner's coursebook
in avoiding errors before they occur, in the late '90's. Amazon has excerpts
at
http://www.amazon.com/exec/obidos/tg/detail/-/0691036632/104-7924268-4227944?vi=glance


> also refer to Burden and Faires and Atkinson's _Numerical Analysis_
> texts. Do you know of any decent web numerical analysis references? I
> have not been able to find much. It would be nice to have some of these
> as well so we could refer to them in the docs and discussion. In any
> case, it is probably a good idea to add a bibliography section to the
> web site.

Not general references (I confess I haven't read a lot of them -- they're often
hard reading and not very helpful in terms of actually implementing anything,
e.g., Bulirsch and Stoer); we should make a point of looking for some.
_Numerical Recipes_ is available online at http://www.nr.com/ . Herbert Wilf
has made several of his books available for download, though they're
sufficiently advanced and specialized (though quite readable, from my skimming)
that one probably wouldn't need to dive into them right away. The
sci.math.num-analysis FAQ is at
http://www.mathcom.com/corpdir/techinfo.mdir/index.html and lists books (not
necessarily online) at
http://www.mathcom.com/corpdir/techinfo.mdir/scifaq/q165.html#q165, including
this review of Acton's first book:

[Daniel Pick] This book is almost worth its price just for the cathartic
interlude in the middle of the book on what not to compute. You should require
your students to read it, learn it, live it. You may find just giving them the
railroad problem found at the beginning of the book a worthwhile exercise.
[Bill Frensley] Amen, brother! The only complaint that I have about Acton's
interlude is that after demolishing the notion of "fitting exponential data,"
he fails to point out that this is the inverse Laplace transform problem.
Perhaps if everyone read this and made the connection, we would be spared the
monthly "is there any good algorithm for the inverse Laplace transform?"

For the curious, the railroad rail problem is paraphrased by
http://krietz.org/AC/Spring03/MAT182/Calc1.pdf as follows:

"This is a famous problem in numerical analysis adapted from Numerical Methods
That . . . Work, by Forman Acton. You don't need to know any numerical analysis
to solve it. I give you all required information.

"A straight railroad rail 1 mile = 5280 feet long is firmly fixed at both ends
along a flat piece of ground. During the next day, the rail heats up and
expands by one additional foot, to 5281 feet long. The ends of the rail are
still fixed at 5280 feet apart (since the ground is not going to stretch the
same way), so this causes the rail to bow up into an arc of a circle.

"What is the maximum height of the rail above its former position?"


> > I'd like to participate in commons-math in an advisory capacity of some
> sort,
> > as I can't in good faith commit to being able to contribute code to the
> > project. Robert indicated that such a role would be useful to the project,
> so
> > I hope you all feel the same!
>
> Please, please, please! Whatever time you have to review
> implementations, comment on strategies, and patiently point out our
> numerical blunders will be greatly appreciated.

Thanks! Having been a computational physicist and not a numerical
mathematician, I've been a user rather than a producer of numerical algorithms
(though computational science makes one a producer of programs that always seem
to require work to make existing code play together correctly), so I'm heavy on
practice and light on theory, though I audited every numerical computing class
offered, both undergrad and grad, when I was at UCLA. So, I won't be of much
use proving things, but I hope my experience in evaluating algorithms'
appropriateness for specific problems based on careful reading and thought will
help out the project.


> Your frank assessment of project direction and the whole question of
> whether or not we should implement the things on the task list
> (http://jakarta.apache.org/commons/sandbox/math/tasks.html) would also
> be appreciated.
>
> Phil
> >
> > Al

I'll briefly comment on these now, prefacing with my initials [ADC]. I'm not
at all strong in probability and statistics, so I won't have much to say about
those components unless I read up fast.

* Add quantiles (1,5,10,25,50,75,90,95,99) to all Univaraite implementations
and bootstrap confidence intervals to StoredUnivariate implementations.

* Add higher order moments to StoredUnivariate (UnivariateImpl if possible).

* t-test statistic needs to be added and we should probably add the capability
of actually performing t- and chi-square tests at fixed significance levels
(.1, .05, .01, .001).

* numerical approximation of the t- and chi-square distributions to enable
user-supplied significance levels.

* The RealMatrixImpl class is missing some key method implementations. The
critical thing is inversion. We need to implement a numerically sound
inversion algorithm. This will enable solve() and also support general linear
regression.
[ADC] Solution of linear systems is essentially never done via matrix
inversion. Commonly an LU (lower-upper triangular) decomposition is performed,
followed by an iteration through the variables, backsolving for each one in
terms of the preceding ones (the first one being trivially solved as x_i = v_i
via the LU decomposition, where v_i is a known value). However, for those
(rare, if I remember correctly) occasions when a matrix inverse itself is
actually needed, there are good and bad ways -- the least efficient being the
recursive method taught in undergraduate linear algebra.

* ComplexNumber interface and implementation. The only tricky thing here is
making division numerically sound and what extended value topology to adopt.
[ADC] I should have something to say here, having had some exposure to special
function evaluation, but sadly will have to brush up a lot before having
anything to contribute.

* Bivariate Regression, corellation. This could be done with simple formulas
manipulating arrays and this is probably what we should aim for in an initial
release. Down the road, we should use the RealMatrixImpl solve() to support
general linear regression.

* Binomial coefficients An "exact" implementation that is limited to what can
be stored in a long exists. This should be extended to use BigIntegers and
potentially to support logarithmic representations.

* Newton's method for finding roots
[ADC] It would also be good to provide the contextual infrastructure
(bracketing a root) as well as methods that don't rely on being able to compute
the derivative of the function. Someday I would also like to re-create a
complex equation solver based on inverse quadratic interpolation (sadly, I
can't remember the original author's name) I once had in my toolkit.

* Exponential growth and decay
[ADC] Not sure what is supposed to be provided here.

* Polynomial Interpolation (curve fitting)
[ADC] Rational function interpolation and cubic spline, too.

* Sampling from Collections
[ADC] Probably brings up the thorny topic of (good) random number generation,
about which Hamming said, "The generation of random numbers is too important to
be left to chance."

* Addition of a Arithmetic Geometric Mean



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Brent Worden]

> * Newton's method for finding roots
> [ADC] It would also be good to provide the contextual infrastructure
> (bracketing a root) as well as methods that don't rely on being
> able to compute
> the derivative of the function. Someday I would also like to re-create a
> complex equation solver based on inverse quadratic interpolation (sadly, I
> can't remember the original author's name) I once had in my toolkit.
Good points. I would prefer to implement the secant or false-position
methods prior to newton's because of its reliance on the derivative. The
other methods have a wider area of application because of the have fewer
constraints on the input function. Plus, they can be used internally by
other areas of the project, such as inverting cumulative distribution
functions.

> * Sampling from Collections
> [ADC] Probably brings up the thorny topic of (good) random number
> generation,
> about which Hamming said, "The generation of random numbers is
> too important to
> be left to chance."
Knuth (the sorting and searching volume I think) provides some insight and
algorithms for simple random samples. One algorithm that is quite
impressive, can generate a simple random sample from a collection of unknown
size.

Brent Worden
http://www.brent.worden.org

[Al Chou]

--- Brent Worden <***@worden.org> wrote:
> > * Newton's method for finding roots
> > [ADC] It would also be good to provide the contextual infrastructure
> > (bracketing a root) as well as methods that don't rely on being
> > able to compute
> > the derivative of the function. Someday I would also like to re-create a
> > complex equation solver based on inverse quadratic interpolation (sadly, I
> > can't remember the original author's name) I once had in my toolkit.
> Good points. I would prefer to implement the secant or false-position
> methods prior to newton's because of its reliance on the derivative. The
> other methods have a wider area of application because of the have fewer
> constraints on the input function. Plus, they can be used internally by

The best derivative-less algorithm I know of is the van Wijngaarden - Dekker -
Brent method described in _NR_. I have been pondering lately the issue of
implementing algorithms that I know of from _NR_ without violating their
license, which explicitly restricts redistribution of source code based on
their published code. But I don't know whether a clean-room implementation is
worth the time or likely to be as good as a published implementation. Any
thoughts from out there?


> > * Sampling from Collections
> > [ADC] Probably brings up the thorny topic of (good) random number
> > generation,
> > about which Hamming said, "The generation of random numbers is
> > too important to
> > be left to chance."
> Knuth (the sorting and searching volume I think) provides some insight and
> algorithms for simple random samples. One algorithm that is quite
> impressive, can generate a simple random sample from a collection of unknown
> size.
>
> Brent Worden
> http://www.brent.worden.org

Yes, Knuth has some good implementations, as does _NR_, which quotes him. As
someone says later in this thread, we may not be able to do sufficiently better
than what's already in the JDK to warrant trying, though we should check to
make sure what the JDK has before we decide.



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[J.Pietschmann]

Al Chou wrote:
> The best derivative-less algorithm I know of is the van Wijngaarden - Dekker -
> Brent method described in _NR_. I have been pondering lately the issue of
> implementing algorithms that I know of from _NR_ without violating their
> license, which explicitly restricts redistribution of source code based on
> their published code. But I don't know whether a clean-room implementation is
> worth the time or likely to be as good as a published implementation. Any
> thoughts from out there?


I can supply an implementation based on H.M.Antia: "Numerical
Methods for Scientists and Engineers" which does not have such
a restriction (the sample code is FORTRAN, anyway).

Then there is a method which uses quadratic interpolation (instead
of inverse quadratic interpolation as Brent), which requires
solving a quadratic but is still performant if square roots
are cheap (and they are in modern hardware, of the same order as
a FP division). The method is still derivative-less but should
converge quadratically (same as Newton) as long as the function
has a smooth second order derivative. I've never seen this analyzed
in a textbook though.

J.Pietschmann

[Al Chou]

--- "J.Pietschmann" <***@yahoo.de> wrote:
> Al Chou wrote:
> > The best derivative-less algorithm I know of is the van Wijngaarden -
> Dekker -
> > Brent method described in _NR_. I have been pondering lately the issue of
> > implementing algorithms that I know of from _NR_ without violating their
> > license, which explicitly restricts redistribution of source code based on
> > their published code. But I don't know whether a clean-room implementation
> is
> > worth the time or likely to be as good as a published implementation. Any
> > thoughts from out there?
>
>
> I can supply an implementation based on H.M.Antia: "Numerical
> Methods for Scientists and Engineers" which does not have such
> a restriction (the sample code is FORTRAN, anyway).

Cool. FYI, _NR_ was first published in FORTRAN 77 and was ported to all other
languages it supports, so it's not the implementation language that matters
when considering _NR_'s license -- you're not even allowed to port their code
to other languages and distribute the ported code without permission.


> Then there is a method which uses quadratic interpolation (instead
> of inverse quadratic interpolation as Brent), which requires
> solving a quadratic but is still performant if square roots
> are cheap (and they are in modern hardware, of the same order as
> a FP division). The method is still derivative-less but should
> converge quadratically (same as Newton) as long as the function
> has a smooth second order derivative. I've never seen this analyzed
> in a textbook though.

Interesting. Where's that method from?



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Phil Steitz]

Al Chou wrote:
> --- Phil Steitz <***@steitz.com> wrote:
>
>>Al Chou wrote:
>>
>>>Greetings to the commons-math community! I emailed Robert Donkin privately
>>
>>a
>>
>>>few days ago asking, as I take it others have, about the reasoning behind
>>>creating a numerics library from scratch rather than incorporating an
>>
>>existing
>>
>>>one. I think I understand that reasoning now, but despite not wanting to
>>>dampen anyone's enthusiasm, I do want to caution those who have not done a
>>
>>lot
>>
>>>of numerical computing that it can very easily be done wrong.
>>
>>No question about that. As I have stated a couple of times, however, I
>>do not personally see commons-math as a numerics library. There will
>>certainly be numerical considerations to worry about, though, and your
>>point is well taken.
>
>
> I'll have to catch up on your posts, Phil, as I don't immediately understand
> how a math library can be not a numerics library.
>

Of course, what commons-math becomes will be what the contributors
choose to make it; but I don't think that it *has* to evolve into a
primarily numerics package or to aspire to a broad numerics scope. The
random data stuff, for example, is not very numerically intensive, nor
is most of the basic stats material, nor any combinatorics or other
discrete math or logic that we might decide to implement. There is no
question, however, that *some* of the things that we are already doing
involve numerical methods and we need to be careful about this stuff.

>
>
>> A big reason why
>>
>>>there's a lot of legacy code in scientific computing is that it's hard to
>>
>>get
>>
>>>numerical algorithms right, so once you do, there's great inertia towards
>>>changing anything (there are of course other big reasons, such as the fact
>>
>>that
>>
>>>many computational scientists are not actually as computer-savvy as they
>>
>>are
>>
>>>science-savvy, so they're not very facile at creating new code).
>>
>>Whence the wonderful proliferation of Fortran code in Java ;-)
>
>
> Remind me sometime to tell you the story my former supervisor told me about an
> "expert" he once worked with whose personal conception of OOP was completely
> orthogonal to the facilities provided for it in Java, which was the
> implementation language they were working in.
>
>
>
>>>On a more positive note, let me recommend _Numerical Recipes_, _Numerical
>>>Methods that [Usually] Work_ (which incidentally presents an alternate
>>>quadratic formula for use in the regime where the traditional one fails),
>>
>>and
>>
>>>_Real Computing Made REAL_ as easy-to-read and down-to-earth references
>>
>>that
>>
>>>show how easy it is for the naive implementation to be woefully inadequate
>>
>>as
>>
>>>well as teach how to do numerical math correctly.
>>
>>No question that the first of these at least and Knuth are classics to
>>refer to. I am not familiar with the second two -- thanks for the tip. I
>
>
> In the preface to the first editions of _Numerical Recipes_, the authors
> acknowledge the influence of Forman Acton's _Numerical Methods that [Usually]
> Work_ on their presentation style. I discovered the first edition of it by
> accident when I was looking in the library for a copy of _NR_, which they were
> all out of at the time, much to my benefit. It was reissued in 1990, and he
> followed up with _Real Computing..._, essentially a practitioner's coursebook
> in avoiding errors before they occur, in the late '90's. Amazon has excerpts
> at
> http://www.amazon.com/exec/obidos/tg/detail/-/0691036632/104-7924268-4227944?vi=glance
>
>
>
>>also refer to Burden and Faires and Atkinson's _Numerical Analysis_
>>texts. Do you know of any decent web numerical analysis references? I
>>have not been able to find much. It would be nice to have some of these
>>as well so we could refer to them in the docs and discussion. In any
>>case, it is probably a good idea to add a bibliography section to the
>>web site.
>
>
> Not general references (I confess I haven't read a lot of them -- they're often
> hard reading and not very helpful in terms of actually implementing anything,
> e.g., Bulirsch and Stoer); we should make a point of looking for some.
> _Numerical Recipes_ is available online at http://www.nr.com/ . Herbert Wilf
> has made several of his books available for download, though they're
> sufficiently advanced and specialized (though quite readable, from my skimming)
> that one probably wouldn't need to dive into them right away. The
> sci.math.num-analysis FAQ is at
> http://www.mathcom.com/corpdir/techinfo.mdir/index.html and lists books (not
> necessarily online) at
> http://www.mathcom.com/corpdir/techinfo.mdir/scifaq/q165.html#q165, including
> this review of Acton's first book:
>
> [Daniel Pick] This book is almost worth its price just for the cathartic
> interlude in the middle of the book on what not to compute. You should require
> your students to read it, learn it, live it. You may find just giving them the
> railroad problem found at the beginning of the book a worthwhile exercise.
> [Bill Frensley] Amen, brother! The only complaint that I have about Acton's
> interlude is that after demolishing the notion of "fitting exponential data,"
> he fails to point out that this is the inverse Laplace transform problem.
> Perhaps if everyone read this and made the connection, we would be spared the
> monthly "is there any good algorithm for the inverse Laplace transform?"
>
> For the curious, the railroad rail problem is paraphrased by
> http://krietz.org/AC/Spring03/MAT182/Calc1.pdf as follows:
>
> "This is a famous problem in numerical analysis adapted from Numerical Methods
> That . . . Work, by Forman Acton. You don't need to know any numerical analysis
> to solve it. I give you all required information.
>
> "A straight railroad rail 1 mile = 5280 feet long is firmly fixed at both ends
> along a flat piece of ground. During the next day, the rail heats up and
> expands by one additional foot, to 5281 feet long. The ends of the rail are
> still fixed at 5280 feet apart (since the ground is not going to stretch the
> same way), so this causes the rail to bow up into an arc of a circle.
>
> "What is the maximum height of the rail above its former position?"
>
>
>
>>>I'd like to participate in commons-math in an advisory capacity of some
>>
>>sort,
>>
>>>as I can't in good faith commit to being able to contribute code to the
>>>project. Robert indicated that such a role would be useful to the project,
>>
>>so
>>
>>>I hope you all feel the same!
>>
>>Please, please, please! Whatever time you have to review
>>implementations, comment on strategies, and patiently point out our
>>numerical blunders will be greatly appreciated.
>
>
> Thanks! Having been a computational physicist and not a numerical
> mathematician, I've been a user rather than a producer of numerical algorithms
> (though computational science makes one a producer of programs that always seem
> to require work to make existing code play together correctly), so I'm heavy on
> practice and light on theory, though I audited every numerical computing class
> offered, both undergrad and grad, when I was at UCLA. So, I won't be of much
> use proving things, but I hope my experience in evaluating algorithms'
> appropriateness for specific problems based on careful reading and thought will
> help out the project.
>

I am certain that it will. I also have a fair amount of applied
experience, but most of it in other languages (probably obvious from my
code -- he he) and a few years back (like my mathematical career). This
will be interesting.

>
>
>>Your frank assessment of project direction and the whole question of
>>whether or not we should implement the things on the task list
>>(http://jakarta.apache.org/commons/sandbox/math/tasks.html) would also
>>be appreciated.
>>
>>Phil
>>
>>>Al
>>
>
> I'll briefly comment on these now, prefacing with my initials [ADC]. I'm not
> at all strong in probability and statistics, so I won't have much to say about
> those components unless I read up fast.
>
> * Add quantiles (1,5,10,25,50,75,90,95,99) to all Univaraite implementations
> and bootstrap confidence intervals to StoredUnivariate implementations.
>
> * Add higher order moments to StoredUnivariate (UnivariateImpl if possible).
>
> * t-test statistic needs to be added and we should probably add the capability
> of actually performing t- and chi-square tests at fixed significance levels
> (.1, .05, .01, .001).
>
> * numerical approximation of the t- and chi-square distributions to enable
> user-supplied significance levels.
>
> * The RealMatrixImpl class is missing some key method implementations. The
> critical thing is inversion. We need to implement a numerically sound
> inversion algorithm. This will enable solve() and also support general linear
> regression.
> [ADC] Solution of linear systems is essentially never done via matrix
> inversion. Commonly an LU (lower-upper triangular) decomposition is performed,
> followed by an iteration through the variables, backsolving for each one in
> terms of the preceding ones (the first one being trivially solved as x_i = v_i
> via the LU decomposition, where v_i is a known value). However, for those
> (rare, if I remember correctly) occasions when a matrix inverse itself is
> actually needed, there are good and bad ways -- the least efficient being the
> recursive method taught in undergraduate linear algebra.
>
> * ComplexNumber interface and implementation. The only tricky thing here is
> making division numerically sound and what extended value topology to adopt.
> [ADC] I should have something to say here, having had some exposure to special
> function evaluation, but sadly will have to brush up a lot before having
> anything to contribute.
>
> * Bivariate Regression, corellation. This could be done with simple formulas
> manipulating arrays and this is probably what we should aim for in an initial
> release. Down the road, we should use the RealMatrixImpl solve() to support
> general linear regression.
>
> * Binomial coefficients An "exact" implementation that is limited to what can
> be stored in a long exists. This should be extended to use BigIntegers and
> potentially to support logarithmic representations.
>
> * Newton's method for finding roots
> [ADC] It would also be good to provide the contextual infrastructure
> (bracketing a root) as well as methods that don't rely on being able to compute
> the derivative of the function. Someday I would also like to re-create a
> complex equation solver based on inverse quadratic interpolation (sadly, I
> can't remember the original author's name) I once had in my toolkit.

Yes. One thing that I am struggling with here is how to do something
simple and useful without function pointers. How should we think about
the interface for rootfinding (however we do it) in Java?

>
> * Exponential growth and decay
> [ADC] Not sure what is supposed to be provided here.

Just simple computations to support mostly financial applications.

>
> * Polynomial Interpolation (curve fitting)
> [ADC] Rational function interpolation and cubic spline, too.

Yes. Here again, step 0 is how to think about the interface/domain in Java.

>
> * Sampling from Collections
> [ADC] Probably brings up the thorny topic of (good) random number generation,
> about which Hamming said, "The generation of random numbers is too important to
> be left to chance."

This is a place where we could spend *lots* of time, recruit some number
theorists and end up with something marginally better than what ships
with the JDK. If someone wants to do this, great, but I think we can
probably assemble some useful stuff just leveraging the JDK generators
(a la RandomData). Just writing the relatively trivial code to generate
a random sub-collection from a collection using the JDK PRNG would
probably be useful to a lot of applications. It could be, however, that
this really belongs in commons-collections.

>
> * Addition of a Arithmetic Geometric Mean
>
>
>
> Al
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Search - Faster. Easier. Bingo.
> http://search.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Mark R. Diggory]

Phil Steitz wrote:

>> * Sampling from Collections
>> [ADC] Probably brings up the thorny topic of (good) random number
>> generation,
>> about which Hamming said, "The generation of random numbers is too
>> important to
>> be left to chance."
>
>
> This is a place where we could spend *lots* of time, recruit some number
> theorists and end up with something marginally better than what ships
> with the JDK. If someone wants to do this, great, but I think we can
> probably assemble some useful stuff just leveraging the JDK generators
> (a la RandomData). Just writing the relatively trivial code to generate
> a random sub-collection from a collection using the JDK PRNG would
> probably be useful to a lot of applications. It could be, however, that
> this really belongs in commons-collections.
>

Hmm, but again, even in Java there are tested Random Number generation
packages that are already available to the community in the public
domain. Much has already been explored in COLT/RngPack. providing a
number of sound strategies that are by far better than the default
generator in Java. For both random data generation and shuffling.

http://hoschek.home.cern.ch/hoschek/colt/V1.0.3/doc/cern/jet/random/engine/package-summary.html
http://hoschek.home.cern.ch/hoschek/colt/V1.0.3/doc/edu/cornell/lassp/houle/RngPack/package-summary.html

Again I wonder on the position of developers of such packages and thier
possible interest in contributing to Apache.

It does make my head spin when considering the issues that seem to arise
when it comes to using code from the public domain vs building it in
house. Is it more responsible to develop your own code or leverage and
acknowledge someone elses packages and efforts.

On top of this, with publicly published source, duplicating
functionality found in other packages and what is definable as plagarism
is very grey in nature. Say, if I read Mersenne Twister and write a
class that does the same thing, is that plagarism if I don't include the
Colt LGPL license? I certainly used the code as an example to build my
own class.

Is it really out of the realm of possibility to have LGPL and GPL code
mixed in with Apache licensed code? I'm sure this is a hot debate to be
had. I've heard Apache developers saying that Apache is compatable with
GPL and LGPL. If its true, can such licenses "coexist" in the Apache
source tree? Then using other's Random Number generators is just a
matter of keeping thier license/authorship in place on the java source.

-Mark

[Mark R. Diggory]

My apologies for not prefixing the package properly. I'll try not to do
it again...

-Mark

Mark R. Diggory wrote:
> Phil Steitz wrote:
>
>>> * Sampling from Collections
>>> [ADC] Probably brings up the thorny topic of (good) random number
>>> generation,
>>> about which Hamming said, "The generation of random numbers is too
>>> important to
>>> be left to chance."
>>
>>
>>
>> This is a place where we could spend *lots* of time, recruit some
>> number theorists and end up with something marginally better than what
>> ships with the JDK. If someone wants to do this, great, but I think we
>> can probably assemble some useful stuff just leveraging the JDK
>> generators (a la RandomData). Just writing the relatively trivial
>> code to generate a random sub-collection from a collection using the
>> JDK PRNG would probably be useful to a lot of applications. It could
>> be, however, that this really belongs in commons-collections.
>>
>
> Hmm, but again, even in Java there are tested Random Number generation
> packages that are already available to the community in the public
> domain. Much has already been explored in COLT/RngPack. providing a
> number of sound strategies that are by far better than the default
> generator in Java. For both random data generation and shuffling.
>
> http://hoschek.home.cern.ch/hoschek/colt/V1.0.3/doc/cern/jet/random/engine/package-summary.html
>
> http://hoschek.home.cern.ch/hoschek/colt/V1.0.3/doc/edu/cornell/lassp/houle/RngPack/package-summary.html
>
>
> Again I wonder on the position of developers of such packages and thier
> possible interest in contributing to Apache.
>
> It does make my head spin when considering the issues that seem to arise
> when it comes to using code from the public domain vs building it in
> house. Is it more responsible to develop your own code or leverage and
> acknowledge someone elses packages and efforts.
>
> On top of this, with publicly published source, duplicating
> functionality found in other packages and what is definable as plagarism
> is very grey in nature. Say, if I read Mersenne Twister and write a
> class that does the same thing, is that plagarism if I don't include the
> Colt LGPL license? I certainly used the code as an example to build my
> own class.
>
> Is it really out of the realm of possibility to have LGPL and GPL code
> mixed in with Apache licensed code? I'm sure this is a hot debate to be
> had. I've heard Apache developers saying that Apache is compatable with
> GPL and LGPL. If its true, can such licenses "coexist" in the Apache
> source tree? Then using other's Random Number generators is just a
> matter of keeping thier license/authorship in place on the java source.
>
> -Mark
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[J.Pietschmann]

Mark R. Diggory wrote:
> Is it really out of the realm of possibility to have LGPL and GPL code
> mixed in with Apache licensed code? I'm sure this is a hot debate to be
> had.
Yes.

> I've heard Apache developers saying that Apache is compatable with
> GPL
Wrong, definitely.

> and LGPL.
Undecided, but declared "no".

> If its true, can such licenses "coexist" in the Apache
> source tree?
No.

GPL is designed to infect everything which comes near GPL'd
code, so if you want to make sure your code can be used in
closed source software and don't want to risk a FSF law suit,
stay at a safe distance. The more the difference between
executable code, source and data is blurred, the farther
a "safe distance" is. E.g. don't even host GPL'd source on
hardware for which your organization is responsible, lest
someone get inconvenient ideas.

The LGPL was designed at times where dynamic linking was in
its infancy and static machine code linking was the rule.
Because the term "linking" is not well defined for Java, and
the FSF deliberately avoids taking a clear position on what
"linking" in the sense of the LGPL means for Java programs,
keeping a safe distance was declared to be a good idea too.

However, copyright law allows for rerelease by the copyright
owner under another license, or dual licensing. Often just
asking the author (copyright holder) for "changing the license"
of a copy to APL or BSDL is enough to get it granted. Public
domain (unlicensed) software is more dangerous in this regard,
because the copyright owner can revoke the public domain status
at any time and even invalidate works which were derived in the
past.

For further questions, ask the PMC or the board. IIRC there
was also a page about the Apache licensing strategy. Further
discussion should take place at the community or the general
list.

J.Pietschmann

[Stefan Bodewig]

On Mon, 26 May 2003, J. Pietschmann <***@yahoo.de> wrote:

>> I've heard Apache developers saying that Apache is compatable with
>> GPL
> Wrong, definitely.

Uhm, no, not that hard. The ASF says that we think the Apache
Software License should be compatible with the (L)GPL when it comes to
a (L)GPLed distribution. The FSF says no here.

Bundling (L)GPLed code with any other code requires the bundle to be
under (L)GPL as well.

The remaining issue is compile-time or even just runtime dependency on
(L)GPLed code without redistributing the code in question. The GPL is
quite clear that this requires the dependening code to be GPLed as
well, the LGPL is less clear, especially in the Java context, but ...

The Apache Software License as such may be acceptable for the author
of a LGPLed library, but it (the Apache Software License) has
specifically been designed so that anybody else can redistribute our
code under almost any conditions he/she chooses.

The later certainly doesn't satisfy the LGPL requirements which means
that our users get restricted in the coice of their license in a way
that extends the restrictions of the Apache Software License, which is
not acceptable to the ASF.

Stefan

[Al Chou]

--- Phil Steitz <***@steitz.com> wrote:
> Al Chou wrote:
> > ... I'm heavy on
> > practice and light on theory, though I audited every numerical computing
> class
> > offered, both undergrad and grad, when I was at UCLA. So, I won't be of
> much
> > use proving things, but I hope my experience in evaluating algorithms'
> > appropriateness for specific problems based on careful reading and thought
> will
> > help out the project.
> >
>
> I am certain that it will. I also have a fair amount of applied
> experience, but most of it in other languages (probably obvious from my
> code -- he he) and a few years back (like my mathematical career). This
> will be interesting.

Cool. All my numerical code in grad school was in Pascal(!) and FORTRAN 77,
but since then I've written in neither of those. Ruby is my language of choice
these days, though I started with sed, AWK, and PERL in grad school and have
learned enough TCL, VBScript, and Python to be dangerous to myself. I've coded
some in Java and was heavily into trying to be well-read in Java and C++ a
couple of years ago but have given that up for the time being. I also spent
many hours on the train learning Scheme from _The Scheme Programming Language_
(http://www.scheme.com/tspl2d/index.html) on my Handspring but never quite got
the hang of doing things recursively as required by the exercises.


> > * Newton's method for finding roots
> > [ADC] It would also be good to provide the contextual infrastructure
> > (bracketing a root) as well as methods that don't rely on being able to
> compute
> > the derivative of the function. Someday I would also like to re-create a
> > complex equation solver based on inverse quadratic interpolation (sadly, I
> > can't remember the original author's name) I once had in my toolkit.
>
> Yes. One thing that I am struggling with here is how to do something
> simple and useful without function pointers. How should we think about
> the interface for rootfinding (however we do it) in Java?

I'm accustomed to the equivalent of a function pointer being used in FORTRAN.
Is passing as a parameter an object that has methods to provide the required
function values (and derivative[s] if necessary) not a good model?


> > * Exponential growth and decay
> > [ADC] Not sure what is supposed to be provided here.
>
> Just simple computations to support mostly financial applications.

Beyond computing a * exp( b * x ) I can't picture what is needed here. Can you
please give an example?


> > * Polynomial Interpolation (curve fitting)
> > [ADC] Rational function interpolation and cubic spline, too.
>
> Yes. Here again, step 0 is how to think about the interface/domain in Java.

I also want us to keep in mind that there's a distinction between interpolation
and curve fitting. A fitted curve doesn't have to pass through all the input
data points, but an interpolating curve must.


> > * Sampling from Collections
> > [ADC] Probably brings up the thorny topic of (good) random number
> generation,
> > about which Hamming said, "The generation of random numbers is too
> important to
> > be left to chance."
>
> This is a place where we could spend *lots* of time, recruit some number
> theorists and end up with something marginally better than what ships
> with the JDK. If someone wants to do this, great, but I think we can
> probably assemble some useful stuff just leveraging the JDK generators
> (a la RandomData). Just writing the relatively trivial code to generate
> a random sub-collection from a collection using the JDK PRNG would
> probably be useful to a lot of applications. It could be, however, that
> this really belongs in commons-collections.

As I said in an earlier post, we should briefly assess the quality of the JDK's
PRNG before deciding whether it needs a replacement.



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Mark R. Diggory]

Hello Al,

I think any participation concerning quality and soundness of Numerical
Algorithms will always be needed and welcomed. The math commons
component is different than most other components because its
implementation is so involved with using the correct algorithm for a
numerical calculation. This is seldom a question of a opinion or
personal taste and (like you say) more deeply rooted in finite-precision
arithmetic.

I think we're starting out with very rudimentary implementations right
now and working towards both clearly defining numerical accuracy and
optimizing the design.

I may be a bit overzealous to put things into [math]. My last email
tended to suggest adding allot more capabilities than are currently
planed. Phil reclarified these and it help me to get perspective again.

> The Math project shall create and maintain a library of lightweight,
> self-contained mathematics and statistics components addressing the
> most common practical problems not immediately available in the Java
> programming language or commons-lang. The guiding principles for
> commons-math will be:
>
> 1. Real-world application use cases determine priority
> 2. Emphasis on small, easily integrated components rather than large
> libraries with complex dependencies
> 3. All algorithms are fully documented and follow generally accepted
> best practices
> 4. In situations where multiple standard algorithms exist, use the
> Strategy pattern to support multiple implementations
> 5. No external dependencies beyond Commons components and the JDK


I think your observations about legacy code are warranted and (as Phil
pointed out) its clear now that we initially do now want to replicate
other entire legacy projects out there and fall into the trap of
implementing "all kinds" of math functionality when others have
accomplished that already.

In terms of precision, I can already see simple examples of what you are
talking about in our code where we calculate moments in Univar.

sum/n

where if sum << n there will be loss in precision. Say your averaging a
very large number of very small double values. There will be loss of
precision that may be problematic.

-Mark

Al Chou wrote:
> Greetings to the commons-math community! I emailed Robert Donkin privately a
> few days ago asking, as I take it others have, about the reasoning behind
> creating a numerics library from scratch rather than incorporating an existing
> one. I think I understand that reasoning now, but despite not wanting to
> dampen anyone's enthusiasm, I do want to caution those who have not done a lot
> of numerical computing that it can very easily be done wrong. A big reason why
> there's a lot of legacy code in scientific computing is that it's hard to get
> numerical algorithms right, so once you do, there's great inertia towards
> changing anything (there are of course other big reasons, such as the fact that
> many computational scientists are not actually as computer-savvy as they are
> science-savvy, so they're not very facile at creating new code).
>
> As an example, the high school quadratic formula
> r = ( -b +- sqrt( b^2 - 4 * a * c ) ) / ( 2 * a )
> is extremely inaccurate, because of finite-precision arithmetic, when
> 4 * a * c << b^2 ,
> because of the subtraction of nearly-equal values. But in high school they
> never teach you that fact.
>
> On a more positive note, let me recommend _Numerical Recipes_, _Numerical
> Methods that [Usually] Work_ (which incidentally presents an alternate
> quadratic formula for use in the regime where the traditional one fails), and
> _Real Computing Made REAL_ as easy-to-read and down-to-earth references that
> show how easy it is for the naive implementation to be woefully inadequate as
> well as teach how to do numerical math correctly.
>
> I'd like to participate in commons-math in an advisory capacity of some sort,
> as I can't in good faith commit to being able to contribute code to the
> project. Robert indicated that such a role would be useful to the project, so
> I hope you all feel the same!
>
>
> Al
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Search - Faster. Easier. Bingo.
> http://search.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Mark R. Diggory]

In fact the higher the moment you calculate (variance, skew, kurtosis)
on the below set of numbers, the greater the loss of precision, this is
because for values less than one

sumquad << sumcube << sumsq << sum

-Mark

Mark R. Diggory wrote:

> In terms of precision, I can already see simple examples of what you are
> talking about in our code where we calculate moments in Univar.
>
> sum/n
>
> where if sum << n there will be loss in precision. Say your averaging a
> very large number of very small double values. There will be loss of
> precision that may be problematic.
>
> -Mark
>

[Brent Worden]

> In fact the higher the moment you calculate (variance, skew, kurtosis)
> on the below set of numbers, the greater the loss of precision, this is
> because for values less than one
>
> sumquad << sumcube << sumsq << sum
>
> -Mark
>
> Mark R. Diggory wrote:

And when the values are greater than one, we run the additional risk of
overflow. Also, because sumsq and n*xbar^2 are relatively large and
relatively equal, subtracting the two, as done in computing variance,
results in loss of precision as well.

One possible way to limit these problems it by using central moments in lieu
of raw moments. Since central moments are expected values, they tend to
converge to a finite value as the sample size increases. They only time
they wouldn't converge is when the data is drawn from a distribution where
those higher moments don't exist.

There are easy formulas for skewness and kurtosis based on the central
moments which could be used for the stored, univariate implementations:
http://mathworld.wolfram.com/Skewness.html
http://mathworld.wolfram.com/Kurtosis.html

As for the rolling implementations, there might be some more research
involved before using this method because of their memoryless property. But
for starters, the sum and sumsq can easily be replaced with there central
moment counterparts, mean and variance. There are formulas that update those
metrics when a new value is added. Weisberg's "Applied Linear Regression"
outlines two such updating formulas for mean and sum of squares which are
numerically superior to direct computation and the raw moment methods.

mean[0] = 0
mean[m + 1] = mean[m] + ((1 / (m + 1)) * (x[m + 1] - mean[m]))

ss[0] = 0
ss[m + 1] = ss[m] + ((m / (m + 1)) * (x[m + 1] - mean[m])^2)

where mean[m] is the mean for the first m values
x[m] is the m-th value
and ss[m] is the sum of squares for the first m values

The sum of squares formula could then be used to derive a similar formula
for variance:

var[0] = 0.0
var[m + 1] = ((m - 1) / m) * var[m] + ((1 / (m + 1)) * (x[m + 1] -
mean[m])^2)

where var[m] is the sample variance for the first m values


I'd be surprised if similar updating formulas didn't exist for the third and
forth central moments. I'll look into it further.

Brent Worden
http://www.brent.worden.org

[Phil Steitz]

Brent Worden wrote:
>>In fact the higher the moment you calculate (variance, skew, kurtosis)
>>on the below set of numbers, the greater the loss of precision, this is
>>because for values less than one
>>
>>sumquad << sumcube << sumsq << sum
>>
>>-Mark
>>
>>Mark R. Diggory wrote:
>
>
> And when the values are greater than one, we run the additional risk of
> overflow. Also, because sumsq and n*xbar^2 are relatively large and
> relatively equal, subtracting the two, as done in computing variance,
> results in loss of precision as well.
>
> One possible way to limit these problems it by using central moments in lieu
> of raw moments. Since central moments are expected values, they tend to
> converge to a finite value as the sample size increases. They only time
> they wouldn't converge is when the data is drawn from a distribution where
> those higher moments don't exist.
>
> There are easy formulas for skewness and kurtosis based on the central
> moments which could be used for the stored, univariate implementations:
> http://mathworld.wolfram.com/Skewness.html
> http://mathworld.wolfram.com/Kurtosis.html
>
> As for the rolling implementations, there might be some more research
> involved before using this method because of their memoryless property. But
> for starters, the sum and sumsq can easily be replaced with there central
> moment counterparts, mean and variance. There are formulas that update those
> metrics when a new value is added. Weisberg's "Applied Linear Regression"
> outlines two such updating formulas for mean and sum of squares which are
> numerically superior to direct computation and the raw moment methods.
>
Why exactly, are these numerically superior? For what class of
examples? Looks like lots more operations to me, especially in the
UnivariateImpl case where the mean, variance are computed only when
demanded -- i.e., there is no requirement to generate mean[0],
mean[1],...etc. I understand that adding (or better, swapping)
operations can sometimes add precision, but I am having a hard time
seeing exactly where the benefit is in this case, especially given the
amount of additional computation required.

> mean[0] = 0
> mean[m + 1] = mean[m] + ((1 / (m + 1)) * (x[m + 1] - mean[m]))
>
> ss[0] = 0
> ss[m + 1] = ss[m] + ((m / (m + 1)) * (x[m + 1] - mean[m])^2)
>
> where mean[m] is the mean for the first m values
> x[m] is the m-th value
> and ss[m] is the sum of squares for the first m values
>
> The sum of squares formula could then be used to derive a similar formula
> for variance:
>
> var[0] = 0.0
> var[m + 1] = ((m - 1) / m) * var[m] + ((1 / (m + 1)) * (x[m + 1] -
> mean[m])^2)
>
> where var[m] is the sample variance for the first m values
>
>
> I'd be surprised if similar updating formulas didn't exist for the third and
> forth central moments. I'll look into it further.
>
> Brent Worden
> http://www.brent.worden.org
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Al Chou]

[Phil Steitz]

Al Chou wrote:
> --- Phil Steitz <***@steitz.com> wrote:
>
>>Brent Worden wrote:
>>
>>>>In fact the higher the moment you calculate (variance, skew, kurtosis)
>>>>on the below set of numbers, the greater the loss of precision, this is
>>>>because for values less than one
>>>>
>>>>sumquad << sumcube << sumsq << sum
>>>>
>>>>-Mark
>>>>
>>>>Mark R. Diggory wrote:
>>>
>>>
>>>And when the values are greater than one, we run the additional risk of
>>>overflow. Also, because sumsq and n*xbar^2 are relatively large and
>>>relatively equal, subtracting the two, as done in computing variance,
>>>results in loss of precision as well.
>>>
>>>One possible way to limit these problems it by using central moments in
>>
>>lieu
>>
>>>of raw moments. Since central moments are expected values, they tend to
>>>converge to a finite value as the sample size increases. They only time
>>>they wouldn't converge is when the data is drawn from a distribution where
>>>those higher moments don't exist.
>>>
>>>There are easy formulas for skewness and kurtosis based on the central
>>>moments which could be used for the stored, univariate implementations:
>>>http://mathworld.wolfram.com/Skewness.html
>>>http://mathworld.wolfram.com/Kurtosis.html
>>>
>>>As for the rolling implementations, there might be some more research
>>>involved before using this method because of their memoryless property.
>>
>>But
>>
>>>for starters, the sum and sumsq can easily be replaced with there central
>>>moment counterparts, mean and variance. There are formulas that update
>>
>>those
>>
>>>metrics when a new value is added. Weisberg's "Applied Linear Regression"
>>>outlines two such updating formulas for mean and sum of squares which are
>>>numerically superior to direct computation and the raw moment methods.
>>>
>>
>>Why exactly, are these numerically superior? For what class of
>>examples? Looks like lots more operations to me, especially in the
>>UnivariateImpl case where the mean, variance are computed only when
>>demanded -- i.e., there is no requirement to generate mean[0],
>>mean[1],...etc. I understand that adding (or better, swapping)
>>operations can sometimes add precision, but I am having a hard time
>>seeing exactly where the benefit is in this case, especially given the
>>amount of additional computation required.
>
>
> _Numerical Recipes in C_, 2nd ed. p. 613
> (http://lib-www.lanl.gov/numerical/bookcpdf/c14-1.pdf) explains:
> "Many textbooks use the binomial theorem to expand out the definitions [of
> statistical quantities] into sums of various powers of the data, ...[,] but
> this can magnify the roundoff error by a large factor... . A clever way to
> minimize roundoff error, especially for large samples, is to use the corrected
> two-pass algorithm [1]: First calculate x[bar, the mean of x], then calculate
> [formula for variance in terms of x - xbar.] ..."
>
> [1] "Algorithms for Computing the Sample Variance: Analysis and
> Recommendations", Chan, T.F., Golub, G.H., and LeVeque, R.J. 1983, American
> Statistician, vol. 37, pp. 242?247.
>
> Unfortunately, Google could not find the above article online, though it is
> cited by articles in CiteSeer, which led to some more generally interesting
> resources listed below. The corrected two-pass algorithm is only given for a
> non-rolling dataset in _NR_, though Brent may have derived a rolling version in
> his last post in this thread (I'm not sure at a glance, I'd have to get out
> pencil and paper to check).
>
> * "Numerical Analysis", Gordon K. Smyth, May 1997
> http://citeseer.nj.nec.com/91928.html
> Abstract: This article discusses the basic ideas of accuracy and describes
> briefly key topics in numerical analysis which are treated at more depth in
> separate articles. Pointers to available software are given at the end of the
> article.
>
> * Bibliography of Accuracy and Stability of Numerical Algorithms, ...
> http://www.ma.man.ac.uk/~higham/asna/asna2.pdf
>

Thany you, Al!

I am convinced by this for that at least for the variance and higher
order moments, whenever we actually store the values (that would include
UnivariateImpl with finite window), we should use the "corrected one
pass" formula (14.1.8) from
http://lib-www.lanl.gov/numerical/bookcpdf/c14-1.pdf. It is a clever
idea to explicitly correct for error in the mean.

This could easily be patched into AbstractStoreUnivariate. We will have
to think a little bit about how to incorporate it into UnivariateImpl
for just the finite window case.

One point on using _NR_. I am OK with grabbing formula (14.1.8), since
there is a reference attributing the source of the formula to an AmStat
article, which means that the _NR_ cannot claim ownership of the
formula. As long as we include attribution to the AmStat authors, I see
no problem including it. We need to be careful about lifting anything
directly from _NR_ (or anywhere else for that matter) where we do not
have a "generally available, no explicit license" expression of the
algorithm/formula in question. My philosophy on this sort of thing is
that 99% of the algorithms of real practical value are in the public
domain and as long as we concentrate on implementing *algorithms*
instead of imitating/translating/stealing code, we will be OK. For
example, Brent's Cholesky decomp is his own implementation, based on his
understanding of the algorithm so this is OK. Cutting and pasting from
Jama would not be OK (even though Jama is in fact in the public domain).
Similarly for my use of inversion, etc to turn uniform in to
exponential, poission deviates, etc.or anything else that I have
submitted. Fortunately, almost all of the really useful things in
mathematics were developed before all of this "IP license-mania" stuff
took root. I don't worry too much about the descendents of Gauss,
Lagrange, et al coming after us, despite the fact that their claims
would have *much* more merit than most of what people are extracting
licensing fees for today.

I am still unsure about how the "inductive" formula below would improve
things for the no-memory, infinite window version of univariate. It
still looks like a lot of extra computation to me for questionable gain.
Note that to emulate the formula above, you might want to (after some
predetermined initial number of iterates) subtract the (changing)
correction factor from (14.1.8). This would certainly have an
interesting impact on the mean estimate :-). Might be a good idea to
dig up the AmStat article and see what impact this would have on the
arguments there. In any case, I would appreciate any ideas that you
have on how exactly the inductive sums below would improve accuracy in
the no-memory case.

Thanks for looking into this.

Phil

>
>
>>>mean[0] = 0
>>>mean[m + 1] = mean[m] + ((1 / (m + 1)) * (x[m + 1] - mean[m]))
>>>
>>>ss[0] = 0
>>>ss[m + 1] = ss[m] + ((m / (m + 1)) * (x[m + 1] - mean[m])^2)
>>>
>>>where mean[m] is the mean for the first m values
>>> x[m] is the m-th value
>>> and ss[m] is the sum of squares for the first m values
>>>
>>>The sum of squares formula could then be used to derive a similar formula
>>>for variance:
>>>
>>>var[0] = 0.0
>>>var[m + 1] = ((m - 1) / m) * var[m] + ((1 / (m + 1)) * (x[m + 1] -
>>>mean[m])^2)
>>>
>>>where var[m] is the sample variance for the first m values
>>
>
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Search - Faster. Easier. Bingo.
> http://search.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Al Chou]

--- Phil Steitz <***@steitz.com> wrote:
> Al Chou wrote:
> > --- Phil Steitz <***@steitz.com> wrote:
> >
> >>Brent Worden wrote:
> >>
> >>>>In fact the higher the moment you calculate (variance, skew, kurtosis)
> >>>>on the below set of numbers, the greater the loss of precision, this is
> >>>>because for values less than one
> >>>>
> >>>>sumquad << sumcube << sumsq << sum
> >>>>
> >>>>-Mark
> >>>>
> >>>>Mark R. Diggory wrote:
> >>>
> >>>
> >>>And when the values are greater than one, we run the additional risk of
> >>>overflow. Also, because sumsq and n*xbar^2 are relatively large and
> >>>relatively equal, subtracting the two, as done in computing variance,
> >>>results in loss of precision as well.
> >>>
> >>>One possible way to limit these problems it by using central moments in
> >>
> >>lieu
> >>
> >>>of raw moments. Since central moments are expected values, they tend to
> >>>converge to a finite value as the sample size increases. They only time
> >>>they wouldn't converge is when the data is drawn from a distribution where
> >>>those higher moments don't exist.
> >>>
> >>>There are easy formulas for skewness and kurtosis based on the central
> >>>moments which could be used for the stored, univariate implementations:
> >>>http://mathworld.wolfram.com/Skewness.html
> >>>http://mathworld.wolfram.com/Kurtosis.html
> >>>
> >>>As for the rolling implementations, there might be some more research
> >>>involved before using this method because of their memoryless property.
> >>
> >>But
> >>
> >>>for starters, the sum and sumsq can easily be replaced with there central
> >>>moment counterparts, mean and variance. There are formulas that update
> >>
> >>those
> >>
> >>>metrics when a new value is added. Weisberg's "Applied Linear Regression"
> >>>outlines two such updating formulas for mean and sum of squares which are
> >>>numerically superior to direct computation and the raw moment methods.
> >>>
> >>
> >>Why exactly, are these numerically superior? For what class of
> >>examples? Looks like lots more operations to me, especially in the
> >>UnivariateImpl case where the mean, variance are computed only when
> >>demanded -- i.e., there is no requirement to generate mean[0],
> >>mean[1],...etc. I understand that adding (or better, swapping)
> >>operations can sometimes add precision, but I am having a hard time
> >>seeing exactly where the benefit is in this case, especially given the
> >>amount of additional computation required.
> >
> >
> > _Numerical Recipes in C_, 2nd ed. p. 613
> > (http://lib-www.lanl.gov/numerical/bookcpdf/c14-1.pdf) explains:
> > "Many textbooks use the binomial theorem to expand out the definitions [of
> > statistical quantities] into sums of various powers of the data, ...[,] but
> > this can magnify the roundoff error by a large factor... . A clever way to
> > minimize roundoff error, especially for large samples, is to use the
> corrected
> > two-pass algorithm [1]: First calculate x[bar, the mean of x], then
> calculate
> > [formula for variance in terms of x - xbar.] ..."
> >
> > [1] "Algorithms for Computing the Sample Variance: Analysis and
> > Recommendations", Chan, T.F., Golub, G.H., and LeVeque, R.J. 1983, American
> > Statistician, vol. 37, pp. 242?247.
> >
> > Unfortunately, Google could not find the above article online, though it is
> > cited by articles in CiteSeer, which led to some more generally interesting
> > resources listed below. The corrected two-pass algorithm is only given for
> a
> > non-rolling dataset in _NR_, though Brent may have derived a rolling
> version in
> > his last post in this thread (I'm not sure at a glance, I'd have to get out
> > pencil and paper to check).
> >
> > * "Numerical Analysis", Gordon K. Smyth, May 1997
> > http://citeseer.nj.nec.com/91928.html
> > Abstract: This article discusses the basic ideas of accuracy and describes
> > briefly key topics in numerical analysis which are treated at more depth in
> > separate articles. Pointers to available software are given at the end of
> the
> > article.
> >
> > * Bibliography of Accuracy and Stability of Numerical Algorithms, ...
> > http://www.ma.man.ac.uk/~higham/asna/asna2.pdf
> >
>
> Thany you, Al!
>
> I am convinced by this for that at least for the variance and higher
> order moments, whenever we actually store the values (that would include
> UnivariateImpl with finite window), we should use the "corrected one
> pass" formula (14.1.8) from
> http://lib-www.lanl.gov/numerical/bookcpdf/c14-1.pdf. It is a clever
> idea to explicitly correct for error in the mean.

Glad to help.


> One point on using _NR_. I am OK with grabbing formula (14.1.8), since
> there is a reference attributing the source of the formula to an AmStat
> article, which means that the _NR_ cannot claim ownership of the
> formula. As long as we include attribution to the AmStat authors, I see
> no problem including it. We need to be careful about lifting anything
> directly from _NR_ (or anywhere else for that matter) where we do not
> have a "generally available, no explicit license" expression of the
> algorithm/formula in question. My philosophy on this sort of thing is
> that 99% of the algorithms of real practical value are in the public
> domain and as long as we concentrate on implementing *algorithms*
> instead of imitating/translating/stealing code, we will be OK. For
> example, Brent's Cholesky decomp is his own implementation, based on his
> understanding of the algorithm so this is OK. Cutting and pasting from
> Jama would not be OK (even though Jama is in fact in the public domain).
> Similarly for my use of inversion, etc to turn uniform in to
> exponential, poission deviates, etc.or anything else that I have
> submitted. Fortunately, almost all of the really useful things in
> mathematics were developed before all of this "IP license-mania" stuff
> took root. I don't worry too much about the descendents of Gauss,
> Lagrange, et al coming after us, despite the fact that their claims
> would have *much* more merit than most of what people are extracting
> licensing fees for today.

I agree we need to be scrupulously careful about how we use _NR_ as a resource,
and in an earlier post I had asked about how we should. I completely agree
with your assessment of the legitimacy of using a published formula that is
simply quoted by the authors of _NR_. My problem if I personally were
contributing code to commons-math is that I've lived so long with _NR_ that it
would be difficult to distance myself from the actual code presented in the
books. Guess it's a good thing I don't really have time to code for the
project. If I do end up coding, I may have to take pains to write new code in
TDD fashion or some such, to avoid any licensing issues that might arise if I
attempt to write code using a reference program.

What I'm not clear on, and I'm not sure I agree with J. Pietschmann's
assessment, is the allowability of incorporating code whose authors have
explicitly placed it in the public domain or licensed in a manner compatible
with Apache. But my opinions are heavily influenced by my career in
computational physics, where implementing an algorithm from scratch was the
last thing we ever wanted to do, given that almost any algorithm one would want
to use (not counting new ones or ones not widely used) has already been coded
and tested better already by someone else whose primary interest at the time
was exactly that. The buy vs. build decision in computational physics was
heavily biased towards "buy", because our interest was physics, not necessarily
numerical mathematics for its own sake. Had I been a numericist primarily, I'd
probably be more biased towards "build". Still, I can't help wanting to avoid
duplication when the new duplicate is likely to be of lower quality than code
that already exists and is free to use.


> I am still unsure about how the "inductive" formula below would improve
> things for the no-memory, infinite window version of univariate. It
> still looks like a lot of extra computation to me for questionable gain.
> Note that to emulate the formula above, you might want to (after some
> predetermined initial number of iterates) subtract the (changing)
> correction factor from (14.1.8). This would certainly have an
> interesting impact on the mean estimate :-). Might be a good idea to
> dig up the AmStat article and see what impact this would have on the
> arguments there. In any case, I would appreciate any ideas that you
> have on how exactly the inductive sums below would improve accuracy in
> the no-memory case.
>
> Thanks for looking into this.
>
> Phil
>
> >
> >
> >>>mean[0] = 0
> >>>mean[m + 1] = mean[m] + ((1 / (m + 1)) * (x[m + 1] - mean[m]))
> >>>
> >>>ss[0] = 0
> >>>ss[m + 1] = ss[m] + ((m / (m + 1)) * (x[m + 1] - mean[m])^2)
> >>>
> >>>where mean[m] is the mean for the first m values
> >>> x[m] is the m-th value
> >>> and ss[m] is the sum of squares for the first m values
> >>>
> >>>The sum of squares formula could then be used to derive a similar formula
> >>>for variance:
> >>>
> >>>var[0] = 0.0
> >>>var[m + 1] = ((m - 1) / m) * var[m] + ((1 / (m + 1)) * (x[m + 1] -
> >>>mean[m])^2)
> >>>
> >>>where var[m] is the sample variance for the first m values

I'll have to do some more research and/or symbolic calculations by hand to see
what I think about the inductive formulas, which I'd never seen before, in the
infinite window case.



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Mark R. Diggory]

Al Chou wrote:
> I agree we need to be scrupulously careful about how we use _NR_ as a resource,
> and in an earlier post I had asked about how we should. I completely agree
> with your assessment of the legitimacy of using a published formula that is
> simply quoted by the authors of _NR_. My problem if I personally were
> contributing code to commons-math is that I've lived so long with _NR_ that it
> would be difficult to distance myself from the actual code presented in the
> books. Guess it's a good thing I don't really have time to code for the
> project. If I do end up coding, I may have to take pains to write new code in
> TDD fashion or some such, to avoid any licensing issues that might arise if I
> attempt to write code using a reference program.
>
> What I'm not clear on, and I'm not sure I agree with J. Pietschmann's
> assessment, is the allowability of incorporating code whose authors have
> explicitly placed it in the public domain or licensed in a manner compatible
> with Apache. But my opinions are heavily influenced by my career in
> computational physics, where implementing an algorithm from scratch was the
> last thing we ever wanted to do, given that almost any algorithm one would want
> to use (not counting new ones or ones not widely used) has already been coded
> and tested better already by someone else whose primary interest at the time
> was exactly that. The buy vs. build decision in computational physics was
> heavily biased towards "buy", because our interest was physics, not necessarily
> numerical mathematics for its own sake. Had I been a numericist primarily, I'd
> probably be more biased towards "build". Still, I can't help wanting to avoid
> duplication when the new duplicate is likely to be of lower quality than code
> that already exists and is free to use.

I tend to agree, the goal is again, an Open Source library of really
solid implementations. I believe if authorship and acknowledgment is
maintained, the authors of public domain code may have an interest in
seeing the inclusion and improvement of their implementations. I think
this should be done in a case-by-case fashion. I've currently taken a
more proactive approach and have begun to approach some of the authors
at CERN and RndPack about their possible interest in becoming involved
and/or providing code from their current projects. If it is available,
if we can get permission, and if we can surmount any license
transitions. Then it would be more beneficial to utilize that which is
available. Its important to note that everything that is currently in
the package was "donated" by the current development team and
contributors. Community is the Apache Spirit.

-Mark

[Al Chou]

--- "Mark R. Diggory" <***@latte.harvard.edu> wrote:
> Al Chou wrote:
> > I agree we need to be scrupulously careful about how we use _NR_ as a
> resource,
> > and in an earlier post I had asked about how we should. I completely agree
> > with your assessment of the legitimacy of using a published formula that is
> > simply quoted by the authors of _NR_. My problem if I personally were
> > contributing code to commons-math is that I've lived so long with _NR_ that
> it
> > would be difficult to distance myself from the actual code presented in the
> > books. Guess it's a good thing I don't really have time to code for the
> > project. If I do end up coding, I may have to take pains to write new code
> in
> > TDD fashion or some such, to avoid any licensing issues that might arise if
> I
> > attempt to write code using a reference program.
> >
> > What I'm not clear on, and I'm not sure I agree with J. Pietschmann's
> > assessment, is the allowability of incorporating code whose authors have
> > explicitly placed it in the public domain or licensed in a manner
> compatible
> > with Apache. But my opinions are heavily influenced by my career in
> > computational physics, where implementing an algorithm from scratch was the
> > last thing we ever wanted to do, given that almost any algorithm one would
> want
> > to use (not counting new ones or ones not widely used) has already been
> coded
> > and tested better already by someone else whose primary interest at the
> time
> > was exactly that. The buy vs. build decision in computational physics was
> > heavily biased towards "buy", because our interest was physics, not
> necessarily
> > numerical mathematics for its own sake. Had I been a numericist primarily,
> I'd
> > probably be more biased towards "build". Still, I can't help wanting to
> avoid
> > duplication when the new duplicate is likely to be of lower quality than
> code
> > that already exists and is free to use.
>
> I tend to agree, the goal is again, an Open Source library of really
> solid implementations. I believe if authorship and acknowledgment is
> maintained, the authors of public domain code may have an interest in
> seeing the inclusion and improvement of their implementations. I think
> this should be done in a case-by-case fashion. I've currently taken a
> more proactive approach and have begun to approach some of the authors
> at CERN and RndPack about their possible interest in becoming involved
> and/or providing code from their current projects. If it is available,
> if we can get permission, and if we can surmount any license
> transitions. Then it would be more beneficial to utilize that which is
> available. Its important to note that everything that is currently in
> the package was "donated" by the current development team and
> contributors. Community is the Apache Spirit.
>
> -Mark

I hope your efforts are fruitful, Mark. You raise a good point, in that the
authors whose work comprises COLT probably never imagined that the license they
agreed to distribute their code under would end up hindering others trying to
use it. Maybe it's a matter of talking to them and getting them to release
under multiple licenses.


Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Phil Steitz]

Al Chou wrote:
> --- "Mark R. Diggory" <***@latte.harvard.edu> wrote:
>
>>Al Chou wrote:
>>
>>>I agree we need to be scrupulously careful about how we use _NR_ as a
>>
>>resource,
>>
>>>and in an earlier post I had asked about how we should. I completely agree
>>>with your assessment of the legitimacy of using a published formula that is
>>>simply quoted by the authors of _NR_. My problem if I personally were
>>>contributing code to commons-math is that I've lived so long with _NR_ that
>>
>>it
>>
>>>would be difficult to distance myself from the actual code presented in the
>>>books. Guess it's a good thing I don't really have time to code for the
>>>project. If I do end up coding, I may have to take pains to write new code
>>
>>in
>>
>>>TDD fashion or some such, to avoid any licensing issues that might arise if
>>
>>I
>>
>>>attempt to write code using a reference program.
>>>
>>>What I'm not clear on, and I'm not sure I agree with J. Pietschmann's
>>>assessment, is the allowability of incorporating code whose authors have
>>>explicitly placed it in the public domain or licensed in a manner
>>
>>compatible
>>
>>>with Apache. But my opinions are heavily influenced by my career in
>>>computational physics, where implementing an algorithm from scratch was the
>>>last thing we ever wanted to do, given that almost any algorithm one would
>>
>>want
>>
>>>to use (not counting new ones or ones not widely used) has already been
>>
>>coded
>>
>>>and tested better already by someone else whose primary interest at the
>>
>>time
>>
>>>was exactly that. The buy vs. build decision in computational physics was
>>>heavily biased towards "buy", because our interest was physics, not
>>
>>necessarily
>>
>>>numerical mathematics for its own sake. Had I been a numericist primarily,
>>
>>I'd
>>
>>>probably be more biased towards "build". Still, I can't help wanting to
>>
>>avoid
>>
>>>duplication when the new duplicate is likely to be of lower quality than
>>
>>code
>>
>>>that already exists and is free to use.
>>
>>I tend to agree, the goal is again, an Open Source library of really
>>solid implementations. I believe if authorship and acknowledgment is
>>maintained, the authors of public domain code may have an interest in
>>seeing the inclusion and improvement of their implementations. I think
>>this should be done in a case-by-case fashion. I've currently taken a
>>more proactive approach and have begun to approach some of the authors
>>at CERN and RndPack about their possible interest in becoming involved
>>and/or providing code from their current projects. If it is available,
>>if we can get permission, and if we can surmount any license
>>transitions. Then it would be more beneficial to utilize that which is
>>available. Its important to note that everything that is currently in
>>the package was "donated" by the current development team and
>>contributors. Community is the Apache Spirit.
>>
>>-Mark
>
>
> I hope your efforts are fruitful, Mark. You raise a good point, in that the
> authors whose work comprises COLT probably never imagined that the license they
> agreed to distribute their code under would end up hindering others trying to
> use it. Maybe it's a matter of talking to them and getting them to release
> under multiple licenses.
>
From my perspective, "relicensing" by itself will not do us much good.
We also have to think about refactoring, versioning and support. If
people want to contribute and can get past whatever legal hurdles they
need to donate their stuff, great. But I think it will be best to have
contributors join the community -- and be willing to put some effort
into refactoring/extending as necessary. We need to keep focused on
assembling an easy-to-use, high-value toolkit. If we get tangled up in
dependencies on large libraries or bring in too much diverse and/or
specialized material, this will make commons-math both harder to use and
harder to maintain.

I think that we should keep focused on getting an intitial release
together that does not include substantially more than what we have on
the task list. There is not really that much left to implement to get
us there. Once we have gotten past the huge hurdle of a Jakarta
release, we can listen to the feedback of the developers who use
commons-math and use this feedback to set priorities for what else to add.

Just my not-too-apache-experienced humble opinion.

Phil

> Al
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Search - Faster. Easier. Bingo.
> http://search.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Al Chou]

--- Phil Steitz <***@steitz.com> wrote:
> Al Chou wrote:
> > --- "Mark R. Diggory" <***@latte.harvard.edu> wrote:
> >
> >>Al Chou wrote:
> >>
> >>>I agree we need to be scrupulously careful about how we use _NR_ as a
> >>
> >>resource,
> >>
> >>>and in an earlier post I had asked about how we should. I completely
> agree
> >>>with your assessment of the legitimacy of using a published formula that
> is
> >>>simply quoted by the authors of _NR_. My problem if I personally were
> >>>contributing code to commons-math is that I've lived so long with _NR_
> that
> >>
> >>it
> >>
> >>>would be difficult to distance myself from the actual code presented in
> the
> >>>books. Guess it's a good thing I don't really have time to code for the
> >>>project. If I do end up coding, I may have to take pains to write new
> code
> >>
> >>in
> >>
> >>>TDD fashion or some such, to avoid any licensing issues that might arise
> if
> >>
> >>I
> >>
> >>>attempt to write code using a reference program.
> >>>
> >>>What I'm not clear on, and I'm not sure I agree with J. Pietschmann's
> >>>assessment, is the allowability of incorporating code whose authors have
> >>>explicitly placed it in the public domain or licensed in a manner
> >>
> >>compatible
> >>
> >>>with Apache. But my opinions are heavily influenced by my career in
> >>>computational physics, where implementing an algorithm from scratch was
> the
> >>>last thing we ever wanted to do, given that almost any algorithm one would
> >>
> >>want
> >>
> >>>to use (not counting new ones or ones not widely used) has already been
> >>
> >>coded
> >>
> >>>and tested better already by someone else whose primary interest at the
> >>
> >>time
> >>
> >>>was exactly that. The buy vs. build decision in computational physics was
> >>>heavily biased towards "buy", because our interest was physics, not
> >>
> >>necessarily
> >>
> >>>numerical mathematics for its own sake. Had I been a numericist
> primarily,
> >>
> >>I'd
> >>
> >>>probably be more biased towards "build". Still, I can't help wanting to
> >>
> >>avoid
> >>
> >>>duplication when the new duplicate is likely to be of lower quality than
> >>
> >>code
> >>
> >>>that already exists and is free to use.
> >>
> >>I tend to agree, the goal is again, an Open Source library of really
> >>solid implementations. I believe if authorship and acknowledgment is
> >>maintained, the authors of public domain code may have an interest in
> >>seeing the inclusion and improvement of their implementations. I think
> >>this should be done in a case-by-case fashion. I've currently taken a
> >>more proactive approach and have begun to approach some of the authors
> >>at CERN and RndPack about their possible interest in becoming involved
> >>and/or providing code from their current projects. If it is available,
> >>if we can get permission, and if we can surmount any license
> >>transitions. Then it would be more beneficial to utilize that which is
> >>available. Its important to note that everything that is currently in
> >>the package was "donated" by the current development team and
> >>contributors. Community is the Apache Spirit.
> >>
> >>-Mark
> >
> >
> > I hope your efforts are fruitful, Mark. You raise a good point, in that
> the
> > authors whose work comprises COLT probably never imagined that the license
> they
> > agreed to distribute their code under would end up hindering others trying
> to
> > use it. Maybe it's a matter of talking to them and getting them to release
> > under multiple licenses.
> >
> From my perspective, "relicensing" by itself will not do us much good.
> We also have to think about refactoring, versioning and support. If
> people want to contribute and can get past whatever legal hurdles they
> need to donate their stuff, great. But I think it will be best to have
> contributors join the community -- and be willing to put some effort
> into refactoring/extending as necessary. We need to keep focused on
> assembling an easy-to-use, high-value toolkit. If we get tangled up in
> dependencies on large libraries or bring in too much diverse and/or
> specialized material, this will make commons-math both harder to use and
> harder to maintain.
>
> I think that we should keep focused on getting an intitial release
> together that does not include substantially more than what we have on
> the task list. There is not really that much left to implement to get
> us there. Once we have gotten past the huge hurdle of a Jakarta
> release, we can listen to the feedback of the developers who use
> commons-math and use this feedback to set priorities for what else to add.
>
> Just my not-too-apache-experienced humble opinion.

I think that even within the focused scope of the initial release there are
portions that could be boosted by not having to write original code (for
instance, a Newton's method or other root finder). Many entries found in the
Netlib repository are from ACM publications, which are subject to the following
policy:

Submittal of an algorithm for publication in one of the ACM
Transactions implies that unrestricted use of the algorithm within a
computer is permissible. General permission to copy and distribute
the algorithm without fee is granted provided that the copies are not
made or distributed for direct commercial advantage. The ACM
copyright notice and the title of the publication and its date appear,
and notice is given that copying is by permission of the Association
for Computing Machinery. To copy otherwise, or to republish, requires
a fee and/or specific permission.

Krogh, F. Algorithms Policy. ACM Tran. Math. Softw. 13(1987),
183-186.

Thus it seems to me that high quality implementations with Apache-compatible
licensing could be had by using ACM-published implementations as references.
The above is quoted from a package that transliterated the Dekker root finder
(which was published in ACM Transactions on Mathematical Software) from Algol
to Fortran and C (though I just looked at both transliterations, and they're
horrible to read, so it might take as much work to refactor them into something
comprehensible as it would to implement the algorithm from scratch).


Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[D***@nascopgh.com]

> --- Phil Steitz <***@steitz.com> wrote:
>
> Many entries found in
> the
> Netlib repository are from ACM publications, which
> are subject to the
> following
> policy:
>
> Submittal of an algorithm for publication
> in one of the ACM
> Transactions implies that unrestricted use of
> the algorithm within a
> computer is permissible. General permission
> to copy and distribute
> the algorithm without fee is granted provided
> that the copies are
> not
> made or distributed for direct commercial
> advantage. The
> ACM
> copyright notice and the title of the
> publication and its date
> appear,
> and notice is given that copying is by
> permission of the
> Association
> for Computing Machinery. To copy otherwise, or
> to republish,
> requires
> a fee and/or specific permission.
>

I don't know if I would call that Apache-compatible licensing terms. The
language seems to specifically preclude distribution for profit, which is
a restriction over and above the Apache license.

[Al Chou]

--- ***@nascopgh.com wrote:
> > --- Phil Steitz <***@steitz.com> wrote:
> >
> > Many entries found in
> > the
> > Netlib repository are from ACM publications, which
> > are subject to the
> > following
> > policy:
> >
> > Submittal of an algorithm for publication
> > in one of the ACM
> > Transactions implies that unrestricted use of
> > the algorithm within a
> > computer is permissible. General permission
> > to copy and distribute
> > the algorithm without fee is granted provided
> > that the copies are
> > not
> > made or distributed for direct commercial
> > advantage. The
> > ACM
> > copyright notice and the title of the
> > publication and its date
> > appear,
> > and notice is given that copying is by
> > permission of the
> > Association
> > for Computing Machinery. To copy otherwise, or
> > to republish,
> > requires
> > a fee and/or specific permission.
>
> I don't know if I would call that Apache-compatible licensing terms. The
> language seems to specifically preclude distribution for profit, which is
> a restriction over and above the Apache license.

Darn, I should have read the Apache license first. I didn't realize until this
morning how un-restrictive it is. Oh, well.


Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Mark R. Diggory]

I caution at assuming this is 'restrictive'. I would like to point out
an important clause in this license that you should consider, I've
received emails from ACM members who are following our project (I am
actually an inactive ACM member myself). ACM is a member based
organization just like Apache.

>To copy otherwise, or to republish, requires
>a fee and/or specific permission.
>
Basically what this is saying is "talk to us". ACM is suggesting
involvement and acknowledgment of their efforts in organizing and
archiving these algorithms. I think often these license clauses (while
legally protecting the the license') are also grounds for establishing
'legal' avenues of involvement and partnership. This clause basically
makes it the responsibility of ACM to decide if the terms and conditions
of the license are to be applied. It gives them jurisdiction to alter
the conditions for the license on a case by case basis. As such, if we
have an interest in using ACM material, we should contact ACM and get an
official position on the usage of such material for an Open Source
Apache project and the legal bindings they would want in such a
relationship..

We also have to consider here, what *copies are not made or distributed
for direct commercial advantage* means in this case as well. This really
appears similar (yet vague) to the Apache clause that suggests that you
cannot use the Apache name as a label to promote your project. IE you
can say 'Foobar, Powered by Apache', but you can't say 'Apache FooBar ',
which suggests that Apache endorses your FooBar as an Apache product. I
feel this is suggesting that you can't use the ACM algorithms as a basis
to gain commercial advantage over another product, not that you are
excluded from using them in a commercial product or to use them here?

Perhaps, as ACM has such a large repository of public source. Apache and
ACM should get together and establish an avenue of opinion when it comes
to this case and other cases where ACM algorithms may be applied in
Apache Source code. Remember, the core necessity of Open Source
licensing is about protecting the authors rights, not about restricting
the reuse and development of Open Source code.

-Mark

Al Chou wrote:

>--- ***@nascopgh.com wrote:
>
>
>>>--- Phil Steitz <***@steitz.com> wrote:
>>>
>>>Many entries found in
>>>the
>>>Netlib repository are from ACM publications, which
>>>are subject to the
>>>following
>>>policy:
>>>
>>> Submittal of an algorithm for publication
>>>in one of the ACM
>>> Transactions implies that unrestricted use of
>>>the algorithm within a
>>> computer is permissible. General permission
>>>to copy and distribute
>>> the algorithm without fee is granted provided
>>>that the copies are
>>>not
>>> made or distributed for direct commercial
>>> advantage. The
>>>ACM
>>> copyright notice and the title of the
>>>publication and its date
>>>appear,
>>> and notice is given that copying is by
>>>permission of the
>>>Association
>>> for Computing Machinery. To copy otherwise, or
>>>to republish,
>>>requires
>>> a fee and/or specific permission.
>>>
>>>
>>I don't know if I would call that Apache-compatible licensing terms. The
>>language seems to specifically preclude distribution for profit, which is
>>a restriction over and above the Apache license.
>>
>>
>
>Darn, I should have read the Apache license first. I didn't realize until this
>morning how un-restrictive it is. Oh, well.
>
>
>Al
>
>=====
>Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
>__________________________________
>Do you Yahoo!?
>The New Yahoo! Search - Faster. Easier. Bingo.
>http://search.yahoo.com
>
>---------------------------------------------------------------------
>To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
>For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>
>
>

[Tim O'Brien]

On Tue, 2003-05-27 at 11:00, Mark R. Diggory wrote:
> Basically what this is saying is "talk to us". ACM is suggesting
> involvement and acknowledgment of their efforts in organizing and
> archiving these algorithms...<snip>...Open Source
> Apache project and the legal bindings they would want in such a
> relationship..

-1

I welcome greater participation especially from the ACM, but I also
encourage people not to "jump the gun". Commons-math is a sandbox
component, nothing guarantees promotion to Commons proper. Commons-math
is not attempting to be a repository for all numerical algorithms - at
least not in the current proposal.

The scope as defined in the proposal is very limited to encourage
developers to submit patches focused on a small set of realistic goals.
aimed at a real-world requirements.

I'd like to see us explore an idea of providing an interface to other
implementations of algorithms (i.e. commons-logging), but I think that
is something for *the future*. Right now, we need focus, unit tests,
pages of detailed documentation, and attention to detail.

[Mark R. Diggory]

Tim O'Brien wrote:

>On Tue, 2003-05-27 at 11:00, Mark R. Diggory wrote:
>
>
>>Basically what this is saying is "talk to us". ACM is suggesting
>>involvement and acknowledgment of their efforts in organizing and
>>archiving these algorithms...<snip>...Open Source
>>Apache project and the legal bindings they would want in such a
>>relationship..
>>
>>
>
>-1
>
>I welcome greater participation especially from the ACM, but I also
>encourage people not to "jump the gun". Commons-math is a sandbox
>component, nothing guarantees promotion to Commons proper. Commons-math
>is not attempting to be a repository for all numerical algorithms - at
>least not in the current proposal.
>
>

+1 Yes I do agree that we should not focus on being a repository for all
numerical algoithms.

>The scope as defined in the proposal is very limited to encourage
>developers to submit patches focused on a small set of realistic goals.
>aimed at a real-world requirements.
>
>
IMHO, I don't consider this discussion to be "jumping the gun" or out of
scope of the issues realted to this project. My primary point is that if
your interested in using/developing code based on an algorithm published
by the ACM, that instead of looking at the license and say, nope, can't
do that. You contact the responsible parties at ACM and ask them
directly concerning if you can apply it in this case. Seems overly
"assuming" to do otherwise. Finally, I wanted to state is that there was
also a probible future of activity between ACM and Apache Math
developers considering such cases.

>I'd like to see us explore an idea of providing an interface to other
>implementations of algorithms (i.e. commons-logging), but I think that
>is something for *the future*. Right now, we need focus, unit tests,
>pages of detailed documentation, and attention to detail.
>
>
Yes, I think that would be intersting, I suspect in the future we will
see cases wher ethis can be applied.

Cheers,
-Mark

[Phil Steitz]

--- Tim O'Brien <***@discursive.com> wrote:
> On Tue, 2003-05-27 at 11:00, Mark R. Diggory wrote:
> > Basically what this is saying is "talk to us". ACM is suggesting
> > involvement and acknowledgment of their efforts in organizing and
> > archiving these algorithms...<snip>...Open Source
> > Apache project and the legal bindings they would want in such a
> > relationship..
>
> -1
>
> I welcome greater participation especially from the ACM, but I also
> encourage people not to "jump the gun". Commons-math is a sandbox
> component, nothing guarantees promotion to Commons proper. Commons-math
> is not attempting to be a repository for all numerical algorithms - at
> least not in the current proposal.
>
> The scope as defined in the proposal is very limited to encourage
> developers to submit patches focused on a small set of realistic goals.
> aimed at a real-world requirements.
>
> I'd like to see us explore an idea of providing an interface to other
> implementations of algorithms (i.e. commons-logging), but I think that
> is something for *the future*. Right now, we need focus, unit tests,
> pages of detailed documentation, and attention to detail.
>
I agree. I at least plan to keep plugging away providing simple standard-
algorithm implementations for stuff on the task list and filling out the docs
and unit tests. While things like rootfinding for non-differentiable functions
may eventually have a place and may benefit from algorithms that someone can
claim copyright ownership of, if no one else does it before I get to it, I will
translate my simple newton's method implementation (which is trivial) and
submit it. I would appreciate input on what a nice Java interface would look
like for rootfinding, initally assuming that the function has a derivative, but
ideally extensible to support strategies that do not require differentiability.

Phil

>
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>


__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Brent Worden]

> and unit tests. While things like rootfinding for
> non-differentiable functions
> may eventually have a place and may benefit from algorithms that
> someone can
> claim copyright ownership of, if no one else does it before I get
> to it, I will
> translate my simple newton's method implementation (which is trivial) and
> submit it. I would appreciate input on what a nice Java interface
> would look
> like for rootfinding, initally assuming that the function has a
> derivative, but
> ideally extensible to support strategies that do not require
> differentiability.

In the chi-square patch, I created a root finding utility class. I used the
bisection method to provide a default mechanism for computing inverse CDFs.
It's driven by a simple Function interface. Check it out and see if it's
something you can use or improve.

The relevant types are org.apache.jakarta.commons.math.RootFinding and
org.apache.jakarta.commons.math.Function and there it's utilized in
org.apache.jakarta.commons.math.stat.distribution.AbstractContinuousDistribu
tion.

Let me know what you think.

Brent Worden
http://www.brent.worden.org

[Phil Steitz]

Brent Worden wrote:
>>and unit tests. While things like rootfinding for
>>non-differentiable functions
>>may eventually have a place and may benefit from algorithms that
>>someone can
>>claim copyright ownership of, if no one else does it before I get
>>to it, I will
>>translate my simple newton's method implementation (which is trivial) and
>>submit it. I would appreciate input on what a nice Java interface
>>would look
>>like for rootfinding, initally assuming that the function has a
>>derivative, but
>>ideally extensible to support strategies that do not require
>>differentiability.
>
>
> In the chi-square patch, I created a root finding utility class. I used the
> bisection method to provide a default mechanism for computing inverse CDFs.
> It's driven by a simple Function interface. Check it out and see if it's
> something you can use or improve.
>
> The relevant types are org.apache.jakarta.commons.math.RootFinding and
> org.apache.jakarta.commons.math.Function and there it's utilized in
> org.apache.jakarta.commons.math.stat.distribution.AbstractContinuousDistribu
> tion.
>
> Let me know what you think.
>

Looks fine to me,at least. I was looking for some magical way to avoid
the "Function" interface requirement; but after thinking about this some
more and looking at your implementation, I think that is about as
convenient as we can make it for users. Newton's method or other
algorithms could be added as alternative RootFinding strategies. I
didn't think of bisection, which is a great idea for the CDF inversion
application, since these functions tend to be well-behaved. Thanks for
setting this up.

Phil

> Brent Worden
> http://www.brent.worden.org
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Al Chou]

--- Brent Worden <***@worden.org> wrote:
> > and unit tests. While things like rootfinding for
> > non-differentiable functions
> > may eventually have a place and may benefit from algorithms that
> > someone can
> > claim copyright ownership of, if no one else does it before I get
> > to it, I will
> > translate my simple newton's method implementation (which is trivial) and
> > submit it. I would appreciate input on what a nice Java interface
> > would look
> > like for rootfinding, initally assuming that the function has a
> > derivative, but
> > ideally extensible to support strategies that do not require
> > differentiability.
>
> In the chi-square patch, I created a root finding utility class. I used the
> bisection method to provide a default mechanism for computing inverse CDFs.
> It's driven by a simple Function interface. Check it out and see if it's
> something you can use or improve.
>
> The relevant types are org.apache.jakarta.commons.math.RootFinding and
> org.apache.jakarta.commons.math.Function and there it's utilized in
> org.apache.jakarta.commons.math.stat.distribution.AbstractContinuousDistribu
> tion.
>
> Let me know what you think.
>
> Brent Worden
> http://www.brent.worden.org


Brent,

I think I have a few improvements that came to mind while I was coding a
Ridders' method implementation. Should I send you a diff file privately?


Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
http://calendar.yahoo.com

[J.Pietschmann]

Brent Worden wrote:
> In the chi-square patch, I created a root finding utility class. I used the
> bisection method to provide a default mechanism for computing inverse CDFs.
> It's driven by a simple Function interface. Check it out and see if it's
> something you can use or improve.

Here's my take on root finding. Note that BrentSolver uses only
inverse quadratic interpolation rather than the full Brent algorithm.

THere are a few issues with accuracy, there should really be a
relative accuracy too, as well as some plausiblity checks.

J.Pietschmann

[Phil Steitz]

J.Pietschmann wrote:
> Brent Worden wrote:
>
>> In the chi-square patch, I created a root finding utility class. I
>> used the
>> bisection method to provide a default mechanism for computing inverse
>> CDFs.
>> It's driven by a simple Function interface. Check it out and see if it's
>> something you can use or improve.
>
>
> Here's my take on root finding. Note that BrentSolver uses only
> inverse quadratic interpolation rather than the full Brent algorithm.

Can you provide a math reference desribing this? I have been referring
to Atkinson and I am having a hard time following the implementation.
What exactly do you mean by "only inverse quadratic interpolation"? Am
I correct in assuming that what you mean by this is that you *never* use
Secant or Bisection? Looks to me like this is the case, but I am having
a little difficulty following the code. If it is the case, what impact
does this have on convergence/stability?

>
> THere are a few issues with accuracy, there should really be a
> relative accuracy too, as well as some plausiblity checks.

What do you have in mind re: plausibility?

I guess I like your rootfinding framework better than Brent's (Here I
mean Brent Worden, the famous commons-math contributor, not the
numerical analyst). I suggest that we agree to use your interfaces and
UnivariateRealSolverImpl base,include Brent's bisection implementation
strategy and modify his CDF inversion to use the new rootfinding framework.

I do have a few small questions/comments on the framework:

1. Having solve() *require* brackets makes it impossible (or at least,
un-natural) to add Newton's method or any other method that just starts
with one initial value. Brent W. includes a method that will try to
find a bracket from initial values. I suppose that we could add a
solve() with only one initial value and if necessary push it out to a
bracket. Personally, given the multiple good implementations available
now, I am OK with dropping Newton altogether, so this could be a moot
point at least for the initial release.

2. I am curious as to why the "result" property is there. How exactly
do we expect clients to make use of this?

3. What were you thinking about putting into the base solve()
implementation as diagnostics? Do you mean things like checking to make
sure the bracket is good?

4. Thinking of the developers who will use this stuff, I feel like we
need a way to insulate them from having to think about rootfinding
strategy choice. As Al has pointed out, we are not (at least yet ;-))
in the AI business, so we can't automagically make the best choice for
them; but it might be a good idea to somehow specify a default strategy.
Unfortunately, the "safe" choice would likely be bisection. We
could obviously annoint bisection as the "default" by naming it
something like "Solver", but there might be simpler/better ways to do this.
<flame-invitation> I bet that in most real world applications the
difference in convergence speed is not a big deal and guaranteed
convergence/domain of application is more important than expected speed
of convergence. Consider, for example, our own use in CDF inversion to
get critical values for stat tests. </flame-invitation>.


Phil

>
> J.Pietschmann
>
>
> ------------------------------------------------------------------------
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org

[J.Pietschmann]

Phil Steitz wrote:
> Can you provide a math reference desribing this?

H.M.Antia: Num. Methods for Scientists and Engineers.

> I have been referring
> to Atkinson and I am having a hard time following the implementation.
> What exactly do you mean by "only inverse quadratic interpolation"?

Brent's method requires a bracketed root as a start. The primary
method for shrinking the bracket is inverse quadratic interpolation.
This means you get three points
x0,y0 x1,y1 x2,y2 x0<x2<x1
and interpolate a parabola
x=a*y^2+b*y+c
and because your goal is to find the x for y=0, your next estimate
for the root is c (set y=0). The convergence is of the order 1.8,
which is better than the secant method (~1.4 if non bracketing, 1..1.2
on average if bracketing).
The full Brent algorithm measures how well the interval shrinks, and
resorts to bisection if progress is too slow. I didn't implement this
because I forgot to take the book with me and had only a hazy memory
of the control parameters. All in all the algorithm combines
near-guaranteed convergence (occasionally problems might falsely
detected as ill-conditioned) with near Newton-speed.

> What do you have in mind re: plausibility?
If the interval to search is (x0,x1), then absAccuracy should be of
the order of max(abs(x0),abs(x1))*relAccuracy. Achievable relative
accuracy depends on underlying hardware, although nowadays basically
everyone uses IEEE 754 (means, uh, 52 bit for double? Damn, have to
look it up!). Both relative and absolute accuracy of function
evaluation are also important, which is the reason to let the user
override defaults.
This means if relAcc is given then reject absAcc if
max(abs(x0),abs(x1))*relAcc+c*absAcc == max(abs(x0),abs(x1))*relAcc
for some predermined constant c in 0.2..1.

> I guess I like your rootfinding framework better than Brent's (Here I
> mean Brent Worden, the famous commons-math contributor, not the
> numerical analyst). I suggest that we agree to use your interfaces and
> UnivariateRealSolverImpl base,include Brent's bisection implementation
> strategy and modify his CDF inversion to use the new rootfinding framework.

No argument against :-) I took special care for the JavaDocs, which
seems to pay off...

> I do have a few small questions/comments on the framework:
>
> 1. Having solve() *require* brackets makes it impossible (or at least,
> un-natural) to add Newton's method or any other method that just starts
> with one initial value.

Why? Start with -4000000,[phone] or so. The algorithm is not
*required* to start with a bracket, just with an interval. Individual
solver implementations should document whether they require a bracket.
Apart from this, there is always the possiblity to add another method.

> 2. I am curious as to why the "result" property is there. How exactly
> do we expect clients to make use of this?

The client can ask for the last result as long as no further attempt
to solve for another root was made. I found this handy. I don't insist
on keeping it.

> 3. What were you thinking about putting into the base solve()
> implementation as diagnostics? Do you mean things like checking to make
> sure the bracket is good?\

No, just adding a message indicating that an unimplemented method
was called. Some frameworks don't display the type of the exception
to the user and rely on the message.

> 4. Thinking of the developers who will use this stuff, I feel like we
> need a way to insulate them from having to think about rootfinding
> strategy choice. As Al has pointed out, we are not (at least yet ;-))
> in the AI business, so we can't automagically make the best choice for
> them; but it might be a good idea to somehow specify a default strategy.
> Unfortunately, the "safe" choice would likely be bisection.

Brent's algorithm is quite safe in this respect.
I'll see whether I can complete the implementation near term. Unfortunately
I'll go on vacation on saturday and will be offline until 2003-06-09.


J.Pietschmann

[Brent Worden]

"J.Pietschmann" <***@yahoo.de> wrote in message
news:***@yahoo.de...
> Phil Steitz wrote:
> > Can you provide a math reference desribing this?
>
> H.M.Antia: Num. Methods for Scientists and Engineers.
>
> > I have been referring
> > to Atkinson and I am having a hard time following the implementation.
> > What exactly do you mean by "only inverse quadratic interpolation"?
>
> Brent's method requires a bracketed root as a start. The primary
> method for shrinking the bracket is inverse quadratic interpolation.
> This means you get three points
> x0,y0 x1,y1 x2,y2 x0<x2<x1
> and interpolate a parabola
> x=a*y^2+b*y+c
> and because your goal is to find the x for y=0, your next estimate
> for the root is c (set y=0). The convergence is of the order 1.8,
> which is better than the secant method (~1.4 if non bracketing, 1..1.2
> on average if bracketing).
> The full Brent algorithm measures how well the interval shrinks, and
> resorts to bisection if progress is too slow. I didn't implement this
> because I forgot to take the book with me and had only a hazy memory
> of the control parameters. All in all the algorithm combines
> near-guaranteed convergence (occasionally problems might falsely
> detected as ill-conditioned) with near Newton-speed.
>
> > What do you have in mind re: plausibility?
> If the interval to search is (x0,x1), then absAccuracy should be of
> the order of max(abs(x0),abs(x1))*relAccuracy. Achievable relative
> accuracy depends on underlying hardware, although nowadays basically
> everyone uses IEEE 754 (means, uh, 52 bit for double? Damn, have to
> look it up!). Both relative and absolute accuracy of function
> evaluation are also important, which is the reason to let the user
> override defaults.
> This means if relAcc is given then reject absAcc if
> max(abs(x0),abs(x1))*relAcc+c*absAcc == max(abs(x0),abs(x1))*relAcc
> for some predermined constant c in 0.2..1.
>
> > I guess I like your rootfinding framework better than Brent's (Here I
> > mean Brent Worden, the famous commons-math contributor, not the
> > numerical analyst). I suggest that we agree to use your interfaces and
> > UnivariateRealSolverImpl base,include Brent's bisection implementation
> > strategy and modify his CDF inversion to use the new rootfinding
framework.
>
> No argument against :-) I took special care for the JavaDocs, which
> seems to pay off...

I agree. The this looks like a very solid framework. One suggestion I
would like to make, is instead of a both a firstDirevative and
secondDerivate method to evaluate the derivates. Create a single
getDerivate() method that returns a UnivariateRealFunction. That way if a
user needs the n-th derivate, they just call the getDerivate() method n
times, once on the original function and once on each of the returned
functions. That way, common-math supports creating whatever degree derivate
a
method might need without ever having to change the framework. We provide
the maximum amout of derivate creation support to the user while providing
us with the minimual amount of maintenance.

>
> > I do have a few small questions/comments on the framework:
> >
> > 1. Having solve() *require* brackets makes it impossible (or at least,
> > un-natural) to add Newton's method or any other method that just starts
> > with one initial value.
>
> Why? Start with -4000000,[phone] or so. The algorithm is not
> *required* to start with a bracket, just with an interval. Individual
> solver implementations should document whether they require a bracket.
> Apart from this, there is always the possiblity to add another method.

Why not allow the user to supply any number of initial values and design the
solvers to compensate as best they can when not enough values are provided.
Each algorithm has criteria about the initial values that must be met before
root finding can be attempted. If the user provided initial values do not
satisfy the criteria, the solver should first attempt to morph the initial
values into a set of values that do satisfy the criteria. If this can not
be accomplish, then the solver would raise an exception. This would take
away some of the user's responsibility of knowing how these algorithms
actually work and place it on us to create more robust components.

>
> > 2. I am curious as to why the "result" property is there. How exactly
> > do we expect clients to make use of this?
>
> The client can ask for the last result as long as no further attempt
> to solve for another root was made. I found this handy. I don't insist
> on keeping it.
>
> > 3. What were you thinking about putting into the base solve()
> > implementation as diagnostics? Do you mean things like checking to make
> > sure the bracket is good?\
>
> No, just adding a message indicating that an unimplemented method
> was called. Some frameworks don't display the type of the exception
> to the user and rely on the message.
>
> > 4. Thinking of the developers who will use this stuff, I feel like we
> > need a way to insulate them from having to think about rootfinding
> > strategy choice. As Al has pointed out, we are not (at least yet ;-))
> > in the AI business, so we can't automagically make the best choice for
> > them; but it might be a good idea to somehow specify a default strategy.

A simple and flexible way to make the "best" choice for the user, is to
control the creation of solver instances.
This could be done via a factory or factory method. The creation could be
based on the class of function the user want to evaluate (real, polynomial,
rational, etc.) Also, we could create a solver chain that provides the
ability to link various solvers together so if a solver fails to find a
root, the next solver in the chain is given a chance to solve it. Either,
eventually one of the solvers finds a root or all solvers fail and an
exception is raised. The complexity of creating the "best chain" would be
contained in the factory or factory method, totally hidden from the user.
This again would relieve users of having some numerical method knowledge
about root finding methods.

> > Unfortunately, the "safe" choice would likely be bisection.
>
> Brent's algorithm is quite safe in this respect.

I'm by no means beholden to bisection. I just used it because I needed
something to use for the distribution work. Brent's algorithm is
guarranteed to converge to a root, provided one is bracketed by the initial
values. So, if a bracket can be provided or generated, it is just as good
as bisection as far as safety and it should almost always beat it in speed
of convergence.

> I'll see whether I can complete the implementation near term.
Unfortunately
> I'll go on vacation on saturday and will be offline until 2003-06-09.
>
>
> J.Pietschmann

Brent Worden
http://www.brent.worden.org

[Phil Steitz]

Brent Worden wrote:
> "J.Pietschmann" <***@yahoo.de> wrote in message
> news:***@yahoo.de...
>
>>Phil Steitz wrote:
>>
>>>Can you provide a math reference desribing this?
>>
>>H.M.Antia: Num. Methods for Scientists and Engineers.
>>
>>
>>> I have been referring
>>>to Atkinson and I am having a hard time following the implementation.
>>>What exactly do you mean by "only inverse quadratic interpolation"?
>>
>>Brent's method requires a bracketed root as a start. The primary
>>method for shrinking the bracket is inverse quadratic interpolation.
>>This means you get three points
>> x0,y0 x1,y1 x2,y2 x0<x2<x1
>>and interpolate a parabola
>> x=a*y^2+b*y+c
>>and because your goal is to find the x for y=0, your next estimate
>>for the root is c (set y=0). The convergence is of the order 1.8,
>>which is better than the secant method (~1.4 if non bracketing, 1..1.2
>>on average if bracketing).
>>The full Brent algorithm measures how well the interval shrinks, and
>>resorts to bisection if progress is too slow. I didn't implement this
>>because I forgot to take the book with me and had only a hazy memory
>>of the control parameters. All in all the algorithm combines
>>near-guaranteed convergence (occasionally problems might falsely
>>detected as ill-conditioned) with near Newton-speed.
>>
>>
>>>What do you have in mind re: plausibility?
>>
>>If the interval to search is (x0,x1), then absAccuracy should be of
>>the order of max(abs(x0),abs(x1))*relAccuracy. Achievable relative
>>accuracy depends on underlying hardware, although nowadays basically
>>everyone uses IEEE 754 (means, uh, 52 bit for double? Damn, have to
>>look it up!). Both relative and absolute accuracy of function
>>evaluation are also important, which is the reason to let the user
>>override defaults.
>>This means if relAcc is given then reject absAcc if
>> max(abs(x0),abs(x1))*relAcc+c*absAcc == max(abs(x0),abs(x1))*relAcc
>>for some predermined constant c in 0.2..1.
>>
>>
>>>I guess I like your rootfinding framework better than Brent's (Here I
>>>mean Brent Worden, the famous commons-math contributor, not the
>>>numerical analyst). I suggest that we agree to use your interfaces and
>>>UnivariateRealSolverImpl base,include Brent's bisection implementation
>>>strategy and modify his CDF inversion to use the new rootfinding
>>
> framework.
>
>>No argument against :-) I took special care for the JavaDocs, which
>>seems to pay off...
>
>
> I agree. The this looks like a very solid framework. One suggestion I
> would like to make, is instead of a both a firstDirevative and
> secondDerivate method to evaluate the derivates. Create a single
> getDerivate() method that returns a UnivariateRealFunction. That way if a
> user needs the n-th derivate, they just call the getDerivate() method n
> times, once on the original function and once on each of the returned
> functions. That way, common-math supports creating whatever degree derivate
> a
> method might need without ever having to change the framework. We provide
> the maximum amout of derivate creation support to the user while providing
> us with the minimual amount of maintenance.

Mathematically, that is certainly a natural suggestion. I am not sure,
however, how easy it will be to implement in UnivariateRealFunction
implementations (contract will be tricky regarding existence) or how
often it will actually be used. Given that the rootfinding solutions
that we seem to be converging on (pun intended) for the first release do
not require derivative computations,
one thing we might consider would be to remove derivatives entirely from
the UnivariateRealFunction interface for now.

It may be more natural for an R->R function to mark itself (infinitely?)
differentiable by implementing a DifferentiableUnivariateRealFunction
(getting a bit long, but you get the idea) interface that extends
UnivariateRealFunction than to not throw exceptions or return null or
"NaF" (not a function, of course) when you ask it for its derivative.
>
>
>>>I do have a few small questions/comments on the framework:
>>>
>>>1. Having solve() *require* brackets makes it impossible (or at least,
>>>un-natural) to add Newton's method or any other method that just starts
>>>with one initial value.
>>
>>Why? Start with -4000000,[phone] or so. The algorithm is not
>>*required* to start with a bracket, just with an interval. Individual
>>solver implementations should document whether they require a bracket.
>>Apart from this, there is always the possiblity to add another method.
>
>
> Why not allow the user to supply any number of initial values and design the
> solvers to compensate as best they can when not enough values are provided.
> Each algorithm has criteria about the initial values that must be met before
> root finding can be attempted. If the user provided initial values do not
> satisfy the criteria, the solver should first attempt to morph the initial
> values into a set of values that do satisfy the criteria. If this can not
> be accomplish, then the solver would raise an exception. This would take
> away some of the user's responsibility of knowing how these algorithms
> actually work and place it on us to create more robust components.
>
>
>>>2. I am curious as to why the "result" property is there. How exactly
>>>do we expect clients to make use of this?
>>
>>The client can ask for the last result as long as no further attempt
>>to solve for another root was made. I found this handy. I don't insist
>>on keeping it.
>>
>>
>>>3. What were you thinking about putting into the base solve()
>>>implementation as diagnostics? Do you mean things like checking to make
>>>sure the bracket is good?\
>>
>>No, just adding a message indicating that an unimplemented method
>>was called. Some frameworks don't display the type of the exception
>>to the user and rely on the message.
>>
>>
>>>4. Thinking of the developers who will use this stuff, I feel like we
>>>need a way to insulate them from having to think about rootfinding
>>>strategy choice. As Al has pointed out, we are not (at least yet ;-))
>>>in the AI business, so we can't automagically make the best choice for
>>>them; but it might be a good idea to somehow specify a default strategy.
>>
>
> A simple and flexible way to make the "best" choice for the user, is to
> control the creation of solver instances.
> This could be done via a factory or factory method. The creation could be
> based on the class of function the user want to evaluate (real, polynomial,
> rational, etc.) Also, we could create a solver chain that provides the
> ability to link various solvers together so if a solver fails to find a
> root, the next solver in the chain is given a chance to solve it. Either,
> eventually one of the solvers finds a root or all solvers fail and an
> exception is raised. The complexity of creating the "best chain" would be
> contained in the factory or factory method, totally hidden from the user.
> This again would relieve users of having some numerical method knowledge
> about root finding methods.

Nice idea; but do we really need to do this? See below.

>
>
>>> Unfortunately, the "safe" choice would likely be bisection.
>>
>>Brent's algorithm is quite safe in this respect.
>
>
> I'm by no means beholden to bisection. I just used it because I needed
> something to use for the distribution work. Brent's algorithm is
> guarranteed to converge to a root, provided one is bracketed by the initial
> values. So, if a bracket can be provided or generated, it is just as good
> as bisection as far as safety and it should almost always beat it in speed
> of convergence.


I agree. I suggest that we make the necessary adjustments, test the heck
out of it, beg Al, J. and whatever other numerical analysts we can find
to review it, and annoint Richard Brent's algorithm as our Solver and,
while leaving the extensibility in the framework, we do not add any
other strategies to the initial release.

>
>
>>I'll see whether I can complete the implementation near term.
>
> Unfortunately
>
>>I'll go on vacation on saturday and will be offline until 2003-06-09.
>>
>>
>>J.Pietschmann
>
>
> Brent Worden
> http://www.brent.worden.org
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>
>

[J.Pietschmann]

Brent Worden wrote:
> I agree. The this looks like a very solid framework. One suggestion I
> would like to make, is instead of a both a firstDirevative and
> secondDerivate method to evaluate the derivates. Create a single
> getDerivate() method that returns a UnivariateRealFunction. That way if a
> user needs the n-th derivate, they just call the getDerivate() method n
> times, once on the original function and once on each of the returned
> functions. That way, common-math supports creating whatever degree derivate a
> method might need without ever having to change the framework. We provide
> the maximum amout of derivate creation support to the user while providing
> us with the minimual amount of maintenance.

Given that numerical algorithms try to avoid derivatives
for a number of reasons, a generic method to calculate
potentially arbitrary derivatives seems to be overkill.
In fact, only very few methods even use a second derivative.

Furthermore, the implementor of the function should know best
how the derivatives are calculated. Note that the contract
allows for numerical derivatives (although this is discouraged).
If a function object is returned, implementation is decoupled,
and higher derivatives might be implemented by someone who
does not match the precision and performance perceptions of
the implementor of the base function.

I'd like to keep it all in one object for this reasons.
I know this may seem to be a bit inconvenient for people who
want to find a zero of a derivative of a function they just
see lying around, but you definitely don't want to solve a
numerical derivative for a root, you'll just get garbage.

> Why not allow the user to supply any number of initial values and design the
> solvers to compensate as best they can when not enough values are provided.
> Each algorithm has criteria about the initial values that must be met before
> root finding can be attempted. If the user provided initial values do not
> satisfy the criteria, the solver should first attempt to morph the initial
> values into a set of values that do satisfy the criteria. If this can not
> be accomplish, then the solver would raise an exception. This would take
> away some of the user's responsibility of knowing how these algorithms
> actually work and place it on us to create more robust components.

The user should have some rough ideas of the interval of the root.
It could be quite a large interval, but then the solver may complain.
There is no "best" algorithm, which can take arbitrary functions
and start values and produce a root with a reasonable amount of work.
Also, many functions have multiple roots. An exhaustive search for
all roots is often impossible. Therefore, the interval is mostly to
prevent the solver from finding a root outside of the domain of
interest for the user. There is no way to get the user's domain of
interest without explicit input.


J.Pietschmann

[D***@nascopgh.com]

> "Mark R. Diggory" <***@latte.harvard.edu> wrote:
>
> >To copy otherwise, or to republish, requires
> >a fee and/or specific permission.
> >
> Basically what this is saying is "talk to us". ACM
> is suggesting
> involvement and acknowledgment of their efforts in
> organizing and
> archiving these algorithms. I think often these
> license clauses (while
> legally protecting the the license') are also
> grounds for establishing
> 'legal' avenues of involvement and partnership.

Absolutely. I'm an ACM member as well, and it's a great organization. It's
certainly possible they *might* choose to donate code to the Apache
project. I merely meant to point out that the license has restrictions
that would prohibit simply incorportating the software into an Apache
Group project w/out such permission from ACM.

> As such, if we
> have an interest in using ACM material, we should
> contact ACM and get an
> official position on the usage of such material for
> an Open Source
> Apache project and the legal bindings they would
> want in such a
> relationship..

Well, again, my understanding of the Apache project and its mission leads
me to believe that those "legal bindings" would have to be "you may
release this under the Apache license" w/ no additional restrictions
placed by the ACM.

>
> We also have to consider here, what *copies are not
> made or distributed
> for direct commercial advantage* means in this case
> as well.

IANAL, but the plain language seems quite clear. "Direct commercial
advantage" means you are selling the software as a product (with or w/out
source, alone or in combination w/ other software), not e.g. teaching a
course (for which you and some institution get paid) or using the code for
some other purpose which happens to generate revenue but where the primary
activity is not distribution of the code (like leasing computer time on a
supercomputer which happens to have the software installed as a library --
no distribution there).

>
> Remember, the core necessity of
> Open Source
> licensing is about protecting the authors rights,
> not about restricting
> the reuse and development of Open Source code.
>

Careful there pilgrim! That kind of talk starts license wars (as very
reasonable people can quite strongly disagree about "the core neccessity
of OSS" -- it's the tension between authors' rights and users' rights that
causes the split between GPL proponents and Apache license proponents).

To end on a non-flamewar-inspiring note, the correct thing to do WRT ACM
stuff is obviously to ask the ACM on a case-by-case basis if they'd be
interested in donating it.

Dave

[Mark R. Diggory]

***@nascopgh.com wrote:

>>"Mark R. Diggory" <***@latte.harvard.edu> wrote:
>>
>>
>>
>>>To copy otherwise, or to republish, requires
>>>a fee and/or specific permission.
>>>
>>>
>>>
>>Basically what this is saying is "talk to us". ACM
>>is suggesting
>>involvement and acknowledgment of their efforts in
>>organizing and
>>archiving these algorithms. I think often these
>>license clauses (while
>>legally protecting the the license') are also
>>grounds for establishing
>>'legal' avenues of involvement and partnership.
>>
>>
>
>Absolutely. I'm an ACM member as well, and it's a great organization. It's
>certainly possible they *might* choose to donate code to the Apache
>project. I merely meant to point out that the license has restrictions
>that would prohibit simply incorportating the software into an Apache
>Group project w/out such permission from ACM.
>
>
True

>
>
>>As such, if we
>>have an interest in using ACM material, we should
>>contact ACM and get an
>>official position on the usage of such material for
>>an Open Source
>>Apache project and the legal bindings they would
>>want in such a
>>relationship..
>>
>>
>
>Well, again, my understanding of the Apache project and its mission leads
>me to believe that those "legal bindings" would have to be "you may
>release this under the Apache license" w/ no additional restrictions
>placed by the ACM.
>
>
Yes, that would be quite direct and obvious.

>
>
>>We also have to consider here, what *copies are not
>>made or distributed
>>for direct commercial advantage* means in this case
>>as well.
>>
>>
>
>IANAL, but the plain language seems quite clear. "Direct commercial
>advantage" means you are selling the software as a product (with or w/out
>source, alone or in combination w/ other software), not e.g. teaching a
>course (for which you and some institution get paid) or using the code for
>some other purpose which happens to generate revenue but where the primary
>activity is not distribution of the code (like leasing computer time on a
>supercomputer which happens to have the software installed as a library --
>no distribution there).
>
>
I still am not convinced that direct commercial advantage means --> sell
a tool/source based on it for a profit. This is the danger of poorly
worded licenses. What does commercial *advantage* mean? (Hope I'm not
sounding too "Clintonesque").

>
>
>>Remember, the core necessity of
>>Open Source
>>licensing is about protecting the authors rights,
>>not about restricting
>>the reuse and development of Open Source code.
>>
>>
>>
>
>Careful there pilgrim! That kind of talk starts license wars (as very
>reasonable people can quite strongly disagree about "the core neccessity
>of OSS" -- it's the tension between authors' rights and users' rights that
>causes the split between GPL proponents and Apache license proponents).
>
>
"Well, thems a fightn' words there cowboy!" I can't promise, but I'll
try not to make so many "generalized" statements in the future ;-)

>To end on a non-flamewar-inspiring note, the correct thing to do WRT ACM
>stuff is obviously to ask the ACM on a case-by-case basis if they'd be
>interested in donating it.
>
>Dave
>
>
See, my response to priorities. I do agree.

Cheers,
-Mark

[Stefan Bodewig]

On Tue, 27 May 2003, Mark R. Diggory <***@latte.harvard.edu>
wrote:

>>To copy otherwise, or to republish, requires
>>a fee and/or specific permission.
>>
> Basically what this is saying is "talk to us".

Yep, but it would also say that to people republishing the Apache code
that republishes or copies the licensed code - this is not a
restruction for the ASF (I'm sure they'd we'd get a positive response)
but on the things you could do with the ASF's software.

This may be subtle point, but I want to stress it again. It's not
enough that something can be distributed under the Apache Software
License, there must not be any additional requirements on top of the
Apache Software License if it is going to be distributed by the ASF.

Stefan

[Al Chou]

Hi, Mark,

Thanks for the welcome. Robert was right about Jakarta being focused on
community!

Thanks also for the level-setting. I didn't know exactly how long commons-math
had been in existence (I learned about it in the latest Jakarta newsletter)
until I skimmed last night through the subject lines in the commons-dev list
archive containing the word "math". I'm glad to have started my participation
fairly early in the project.

I think I understand a lot better now the intent of commons-math. In essence,
it's supposed to be a small set of sharp tools for solving the most commonly
encountered math situations in real-world programming rather than an attempt to
build a comprehensive math library.

I'll have to look at Univar to see what you're calculating before I'll have any
comments aside from the general one that generally the higher the moment, the
less often used it is. As for loss of precision in averaging a bunch of small
values, Jon Bentley suggested in _Programming Pearls_ (or was it _More
Programming Pearls_?) summing them in order of ascending magnitude to avoid, as
much as possible, loss of precision. I've never seen that technique mentioned
in numerical analysis circles, where sorting seems to be an extremely rare
operation (in my experience, at least). Unfortunately, if you have some values
of opposite sign and nearly equal magnitude, this technique will not prevent
loss of precision when summing those values. I guess it would be best to sum
values of the same sign only, but that just passes the precision buck to the
final addition of the sum of the positive values to the sum of the negative
values. Without further knowledge of the specific problem, it's difficult to
find a way to eliminate the possibility of subtracting two nearly equal values
(the best way is to analytically determine what the small remainder is when the
big portions that are equal and of opposite sign have been canceled out -- the
key to the railroad rail problem, incidentally).


Al


--- "Mark R. Diggory" <***@latte.harvard.edu> wrote:
> Hello Al,
>
> I think any participation concerning quality and soundness of Numerical
> Algorithms will always be needed and welcomed. The math commons
> component is different than most other components because its
> implementation is so involved with using the correct algorithm for a
> numerical calculation. This is seldom a question of a opinion or
> personal taste and (like you say) more deeply rooted in finite-precision
> arithmetic.
>
> I think we're starting out with very rudimentary implementations right
> now and working towards both clearly defining numerical accuracy and
> optimizing the design.
>
> I may be a bit overzealous to put things into [math]. My last email
> tended to suggest adding allot more capabilities than are currently
> planed. Phil reclarified these and it help me to get perspective again.
>
> > The Math project shall create and maintain a library of lightweight,
> > self-contained mathematics and statistics components addressing the
> > most common practical problems not immediately available in the Java
> > programming language or commons-lang. The guiding principles for
> > commons-math will be:
> >
> > 1. Real-world application use cases determine priority
> > 2. Emphasis on small, easily integrated components rather than large
> > libraries with complex dependencies
> > 3. All algorithms are fully documented and follow generally accepted
> > best practices
> > 4. In situations where multiple standard algorithms exist, use the
> > Strategy pattern to support multiple implementations
> > 5. No external dependencies beyond Commons components and the JDK
>
>
> I think your observations about legacy code are warranted and (as Phil
> pointed out) its clear now that we initially do now want to replicate
> other entire legacy projects out there and fall into the trap of
> implementing "all kinds" of math functionality when others have
> accomplished that already.
>
> In terms of precision, I can already see simple examples of what you are
> talking about in our code where we calculate moments in Univar.
>
> sum/n
>
> where if sum << n there will be loss in precision. Say your averaging a
> very large number of very small double values. There will be loss of
> precision that may be problematic.
>
> -Mark

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
The New Yahoo! Search - Faster. Easier. Bingo.
http://search.yahoo.com

[Shapira, Yoav]

Howdy,

>4. Thinking of the developers who will use this stuff, I feel like we
>need a way to insulate them from having to think about rootfinding
>strategy choice. As Al has pointed out, we are not (at least yet ;-))
>in the AI business, so we can't automagically make the best choice for
>them; but it might be a good idea to somehow specify a default
strategy.
> Unfortunately, the "safe" choice would likely be bisection. We
>could obviously annoint bisection as the "default" by naming it
>something like "Solver", but there might be simpler/better ways to do
this.
><flame-invitation> I bet that in most real world applications the
>difference in convergence speed is not a big deal and guaranteed
>convergence/domain of application is more important than expected speed
>of convergence. Consider, for example, our own use in CDF inversion to
>get critical values for stat tests. </flame-invitation>.

As someone who plans to use commons-math but does not plan to contribute
code to it (though I do plan on testing, commenting in this forum, and
using in production), I definitely need guaranteed convergence more than
convergence speed.

Yoav Shapira



This e-mail, including any attachments, is a confidential business communication, and may contain information that is confidential, proprietary and/or privileged. This e-mail is intended only for the individual(s) to whom it is addressed, and may not be saved, copied, printed, disclosed or used by anyone else. If you are not the(an) intended recipient, please immediately delete this e-mail from your computer system and notify the sender. Thank you.

[Al Chou]

--- "Shapira, Yoav" <***@mpi.com> wrote:
> Howdy,
>
> >4. Thinking of the developers who will use this stuff, I feel like we
> >need a way to insulate them from having to think about rootfinding
> >strategy choice. As Al has pointed out, we are not (at least yet ;-))
> >in the AI business, so we can't automagically make the best choice for
> >them; but it might be a good idea to somehow specify a default
> strategy.
> > Unfortunately, the "safe" choice would likely be bisection. We
> >could obviously annoint bisection as the "default" by naming it
> >something like "Solver", but there might be simpler/better ways to do
> this.
> ><flame-invitation> I bet that in most real world applications the
> >difference in convergence speed is not a big deal and guaranteed
> >convergence/domain of application is more important than expected speed
> >of convergence. Consider, for example, our own use in CDF inversion to
> >get critical values for stat tests. </flame-invitation>.
>
> As someone who plans to use commons-math but does not plan to contribute
> code to it (though I do plan on testing, commenting in this forum, and
> using in production), I definitely need guaranteed convergence more than
> convergence speed.
>
> Yoav Shapira
>
>
>
> This e-mail, including any attachments, is a confidential business
> communication, and may contain information that is confidential, proprietary
> and/or privileged. This e-mail is intended only for the individual(s) to
> whom it is addressed, and may not be saved, copied, printed, disclosed or
> used by anyone else. If you are not the(an) intended recipient, please
> immediately delete this e-mail from your computer system and notify the
> sender. Thank you.


Phil,

Brent never replied about my patches to his bisection code, and now I have an
(as-yet untested) implementation of Ridders' method, which is a simple but
substantial improvement on bisection and not as complex as Brent or its
descendants. Should I just create attachments in Bugzilla so folks can examine
my code?


Al


=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
http://calendar.yahoo.com

[Phil Steitz]

Al Chou wrote:
> --- "Shapira, Yoav" <***@mpi.com> wrote:
>
>>Howdy,
>>
>>
>>>4. Thinking of the developers who will use this stuff, I feel like we
>>>need a way to insulate them from having to think about rootfinding
>>>strategy choice. As Al has pointed out, we are not (at least yet ;-))
>>>in the AI business, so we can't automagically make the best choice for
>>>them; but it might be a good idea to somehow specify a default
>>
>>strategy.
>>
>>> Unfortunately, the "safe" choice would likely be bisection. We
>>>could obviously annoint bisection as the "default" by naming it
>>>something like "Solver", but there might be simpler/better ways to do
>>
>>this.
>>
>>><flame-invitation> I bet that in most real world applications the
>>>difference in convergence speed is not a big deal and guaranteed
>>>convergence/domain of application is more important than expected speed
>>>of convergence. Consider, for example, our own use in CDF inversion to
>>>get critical values for stat tests. </flame-invitation>.
>>
>>As someone who plans to use commons-math but does not plan to contribute
>>code to it (though I do plan on testing, commenting in this forum, and
>>using in production), I definitely need guaranteed convergence more than
>>convergence speed.
>>
>>Yoav Shapira
>>
>>
>>
>>This e-mail, including any attachments, is a confidential business
>>communication, and may contain information that is confidential, proprietary
>>and/or privileged. This e-mail is intended only for the individual(s) to
>>whom it is addressed, and may not be saved, copied, printed, disclosed or
>>used by anyone else. If you are not the(an) intended recipient, please
>>immediately delete this e-mail from your computer system and notify the
>>sender. Thank you.
>
>
>
> Phil,
>
> Brent never replied about my patches to his bisection code, and now I have an
> (as-yet untested) implementation of Ridders' method, which is a simple but
> substantial improvement on bisection and not as complex as Brent or its
> descendants. Should I just create attachments in Bugzilla so folks can examine
> my code?
>
>
> Al
>

Yes. We are getting a bit backed up here. I think the best thing to do
is to wait until one of the committers commits Brent's stuff and then
submit a Bugzilla patch against the cvs source.

I would suggest the following order of events:

1. commit 20279 (Brent's CDF, chi-square stuff, including bisection)

2. submit your patch to Brent's bisection

3. commit J.Pietschmann's submission

4. Brent's stuff refactored to fit into J.Pietschmann's framework

5. Declare rootfinding victory :-) (at least for the first release,
modulo test/doc/validation)


Phil


>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
> http://calendar.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[O'brien, Tim]

Proposal: move all Univariate and Bivariate interfaces and
implementations into the package named "stat".

IMO, this qualifies as a "Long Term Plan", which is subject to lazy
consensus. Are there any binding -1s to moving anything
"**/*Univariate*.java" into stat. (UnivariateFunction.java from Brent's
patch is not involved in this change).

[Phil Steitz]

O'brien, Tim wrote:
> Proposal: move all Univariate and Bivariate interfaces and
> implementations into the package named "stat".
>
> IMO, this qualifies as a "Long Term Plan", which is subject to lazy
> consensus. Are there any binding -1s to moving anything
> "**/*Univariate*.java" into stat. (UnivariateFunction.java from Brent's
> patch is not involved in this change).
>
>
I agree with the proposal. We will probably want to move some more
stuff there as well, but this is a good start.

Phil

>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Mark R. Diggory]

+1 I think its a good starting decision for package structure.

While we're at it, maybe matrix related Classes should get moved into a
similar package for consistency?

-M.

Phil Steitz wrote:

> O'brien, Tim wrote:
>
>> Proposal: move all Univariate and Bivariate interfaces and
>> implementations into the package named "stat".
>> IMO, this qualifies as a "Long Term Plan", which is subject to lazy
>> consensus. Are there any binding -1s to moving anything
>> "**/*Univariate*.java" into stat. (UnivariateFunction.java from Brent's
>> patch is not involved in this change).
>>
>>
> I agree with the proposal. We will probably want to move some more
> stuff there as well, but this is a good start.
>
> Phil
>
>>
>>
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
>> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>>
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Brent Worden]

> Brent,
>
> I think I have a few improvements that came to mind while I was coding a
> Ridders' method implementation. Should I send you a diff file privately?
>
>
> Al
>
> =====
> Albert Davidson Chou

Al,

Sorry I didn't get back to you. I got bogged down doing other things. I
would agree with what Phil suggested and let the committer get caught up
with all the patches, create the patch for your stuff and commit it then.

thanks,

Brent Worden
http://www.brent.worden.org

[O'brien, Tim]

In general, code will be applied if it is submitted as a patch and we
can all agree on an approach. Sending a tarball to the list with some
proposed code is not helpful - I'm not scolding anyone for doing that, I
just want to make it clear that patches are always preferred.

Also, just a reminder - don't use tabs, and always use the LICENSE.txt
as a prefix to all classes.

Thanks



On Fri, 2003-05-30 at 10:22, Brent Worden wrote:
> > Brent,
> >
> > I think I have a few improvements that came to mind while I was coding a
> > Ridders' method implementation. Should I send you a diff file privately?
> >
> >
> > Al
> >
> > =====
> > Albert Davidson Chou
>
> Al,
>
> Sorry I didn't get back to you. I got bogged down doing other things. I
> would agree with what Phil suggested and let the committer get caught up
> with all the patches, create the patch for your stuff and commit it then.
>
> thanks,
>
> Brent Worden
> http://www.brent.worden.org
>
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[O'brien, Tim]

Sorry, some clarification.

I guess what I'm trying to say here is. Please submit patches. Even if
your changes rely on patches which have yet to be committed - you are
allowed to submit a patch which depends on another patch being applied.

Someone recently mailed me a large JAR file in which all Javadoc errors
had been fixed. These changes were a few days out of date, and
attempting to manually reconcile the changes with the current code-base
--- possible, but we all have a limited amount of time.

You'll help Mark and I. MArk, I see you are in the avail file and have
the necessary karma. Update that project.xml file.

Tim

[Mark R. Diggory]

O'brien, Tim wrote:

>Sorry, some clarification.
>
>I guess what I'm trying to say here is. Please submit patches. Even if
>your changes rely on patches which have yet to be committed - you are
>allowed to submit a patch which depends on another patch being applied.
>
>
patches are much easier to apply, even when out of date...

>You'll help Mark and I. MArk, I see you are in the avail file and have
>the necessary karma. Update that project.xml file.
>
>

Hmm, it appears I'm having a little difficulty commiting. I'm going back
through the commiter info to see if I'm missing something. I keep
getting the following error on a freshly checked out version of math.
I've simply edited the project.xml file and am trying to commit it,
updating it doesn't cause any conflicts, so I know its not a conflict
causing this.

The server reported an error while performing the "cvs commit" command.
math: cvs server: Pre-commit check failed
math: cvs [server aborted]: correct above errors first!

Any tips?
-Mark

[Mark R. Diggory]

My checkout was still using /home/cvspublic. Once I switched to
/home/cvs everything worked fine.

-Mark

Mark R. Diggory wrote:
>
> The server reported an error while performing the "cvs commit" command.
> math: cvs server: Pre-commit check failed
> math: cvs [server aborted]: correct above errors first!
>
> Any tips?
> -Mark
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[Tim O'Brien]

Eeek. I'll respond on the list only in that it might be of reference to
someone in the future. 95% of my problems are solved via
groups.google.com...

1. Make sure you are not hitting the public repository - you should
either hit CVS using pserver and portforwarding, or ext with CVS_RSH=ssh

2. Look at the avail file in CVSROOT - this is the font from which karma
flows - make sure your name in on the same line as
jakarta-commons-sandbox

Tim

On Fri, 2003-05-30 at 16:39, Mark R. Diggory wrote:
>
>
> O'brien, Tim wrote:
>
> >Sorry, some clarification.
> >
> >I guess what I'm trying to say here is. Please submit patches. Even if
> >your changes rely on patches which have yet to be committed - you are
> >allowed to submit a patch which depends on another patch being applied.
> >
> >
> patches are much easier to apply, even when out of date...
>
> >You'll help Mark and I. MArk, I see you are in the avail file and have
> >the necessary karma. Update that project.xml file.
> >
> >
>
> Hmm, it appears I'm having a little difficulty commiting. I'm going back
> through the commiter info to see if I'm missing something. I keep
> getting the following error on a freshly checked out version of math.
> I've simply edited the project.xml file and am trying to commit it,
> updating it doesn't cause any conflicts, so I know its not a conflict
> causing this.
>
> The server reported an error while performing the "cvs commit" command.
> math: cvs server: Pre-commit check failed
> math: cvs [server aborted]: correct above errors first!
>
> Any tips?
> -Mark
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>
>

[Mark R. Diggory]

Tim O'Brien wrote:
> Eeek. I'll respond on the list only in that it might be of reference to
> someone in the future. 95% of my problems are solved via
> groups.google.com...
>
> 1. Make sure you are not hitting the public repository - you should
> either hit CVS using pserver and portforwarding, or ext with CVS_RSH=ssh
>

Yes, I was dragging around /home/cvspublic directories from my old
anonymous checkout. This is resolved in the following line of perl on
the commiters page, since I didn't use it to convert the anonymous
checkout on my crappy Windowz laptop. I ran into problems

% cd your-project-working-directory
% find . -type d -name CVS \
> | perl -pe 's,(.*CVS$),$1/Root,' \
> | xargs perl -pi -e \
> 's,:pserver:anoncvs(.+:/home/cvs)public,:ext:username$1,'
% cvs -q update -dP

thanks again for the tip. -M.

> 2. Look at the avail file in CVSROOT - this is the font from which karma
> flows - make sure your name in on the same line as
> jakarta-commons-sandbox
>
> Tim
>
> On Fri, 2003-05-30 at 16:39, Mark R. Diggory wrote:
>
>>
>>O'brien, Tim wrote:
>>
>>
>>>Sorry, some clarification.
>>>
>>>I guess what I'm trying to say here is. Please submit patches. Even if
>>>your changes rely on patches which have yet to be committed - you are
>>>allowed to submit a patch which depends on another patch being applied.
>>>
>>>
>>
>>patches are much easier to apply, even when out of date...
>>
>>
>>>You'll help Mark and I. MArk, I see you are in the avail file and have
>>>the necessary karma. Update that project.xml file.
>>>
>>>
>>
>>Hmm, it appears I'm having a little difficulty commiting. I'm going back
>>through the commiter info to see if I'm missing something. I keep
>>getting the following error on a freshly checked out version of math.
>>I've simply edited the project.xml file and am trying to commit it,
>>updating it doesn't cause any conflicts, so I know its not a conflict
>>causing this.
>>
>>The server reported an error while performing the "cvs commit" command.
>> math: cvs server: Pre-commit check failed
>> math: cvs [server aborted]: correct above errors first!
>>
>>Any tips?
>>-Mark
>>
>>
>>
>>---------------------------------------------------------------------
>>To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
>>For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>>
>>
>
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[b***@worden.org]

On Fri, 30 May 2003 20:14:23 +0200, "J.Pietschmann" wrote:

>
> Brent Worden wrote:
> > I agree. The this looks like a very solid framework. One
> suggestion I
> > would like to make, is instead of a both a firstDirevative and
> > secondDerivate method to evaluate the derivates. Create a single
> > getDerivate() method that returns a UnivariateRealFunction. That
> way if a
> > user needs the n-th derivate, they just call the getDerivate()
> method n
> > times, once on the original function and once on each of the
returned
> > functions. That way, common-math supports creating whatever degree
> derivate a
> > method might need without ever having to change the framework. We
> provide
> > the maximum amout of derivate creation support to the user while
> providing
> > us with the minimual amount of maintenance.
>
> Given that numerical algorithms try to avoid derivatives
> for a number of reasons, a generic method to calculate
> potentially arbitrary derivatives seems to be overkill.
> In fact, only very few methods even use a second derivative.

Just because we the creators of the framework see no need for higher
order derivates, the end users, whose needs no one can fully
enumerate, may have a quite reasonable use for them. Your approach
limits them to the 2nd derivate with no easy means a creating anything
higher. My approach doesn't limit the end user and gives them the
flexibility to create any derivate they feel they need.

>
> Furthermore, the implementor of the function should know best
> how the derivatives are calculated. Note that the contract
> allows for numerical derivatives (although this is discouraged).
> If a function object is returned, implementation is decoupled,

How are they decoupled. If you create a function, you control the
implementation for creating the derivative function object and is
directly coupled to your object. The implementor, you, is in control
of all details of the function and its derivate.

> and higher derivatives might be implemented by someone who
> does not match the precision and performance perceptions of
> the implementor of the base function.

Since the same implementor is in control of both the function and its
derivate that is fault of that one implementor. Furthermore, why
should we care? If the user wants to create a function with less or
more precision or performance, whose to stop them? Your approach
can't prevent it and it might be all the precision/performance they
need.

>
> I'd like to keep it all in one object for this reasons.
> I know this may seem to be a bit inconvenient for people who
> want to find a zero of a derivative of a function they just
> see lying around, but you definitely don't want to solve a
> numerical derivative for a root, you'll just get garbage.

Again, you're pigeonholing the needs of the user. The user might have
needs we can't possibly anticipate. That's the exactly reason we
should make creating derivatives convenient and easy for the user. We
shouldn't take the arrogant position that we know best and prohibit
users from doing the things they want to do.

> > Why not allow the user to supply any number of initial values and
> design the
> > solvers to compensate as best they can when not enough values are
> provided.
> > Each algorithm has criteria about the initial values that must be
> met before
> > root finding can be attempted. If the user provided initial values
> do not
> > satisfy the criteria, the solver should first attempt to morph the
> initial
> > values into a set of values that do satisfy the criteria. If this
> can not
> > be accomplish, then the solver would raise an exception. This
would
> take
> > away some of the user's responsibility of knowing how these
> algorithms
> > actually work and place it on us to create more robust components.
>
> The user should have some rough ideas of the interval of the root.
> It could be quite a large interval, but then the solver may complain.
> There is no "best" algorithm, which can take arbitrary functions
> and start values and produce a root with a reasonable amount of
> work.

I never said anything about a best algorithm. What i did say is every
solver should take an arbitrary function and start value(s) and make
its best effort to return a root. If that means generating starting
values, the solver should have the intelligence to do such a thing.

> Also, many functions have multiple roots. An exhaustive search for
> all roots is often impossible.

When did I introduce finding all roots into this discussion?

> Therefore, the interval is mostly to
> prevent the solver from finding a root outside of the domain of
> interest for the user. There is no way to get the user's domain of
> interest without explicit input.

Why can't the domain of interest and explicit input be a single value?
That's exactly how I did the CDF stuff using a root finding method
which needs a bracket and it worked fine. And it's a generic process
that could be incorporated into all the solvers.

Brent Worden
http://www.brent.worden.org/

[Al Chou]

--- ***@worden.org wrote:
> On Fri, 30 May 2003 20:14:23 +0200, "J.Pietschmann" wrote:
> > Brent Worden wrote:
> > > I agree. The this looks like a very solid framework. One
> > suggestion I
> > > would like to make, is instead of a both a firstDirevative and
> > > secondDerivate method to evaluate the derivates. Create a single
> > > getDerivate() method that returns a UnivariateRealFunction. That
> > way if a
> > > user needs the n-th derivate, they just call the getDerivate()
> > method n
> > > times, once on the original function and once on each of the
> returned
> > > functions. That way, common-math supports creating whatever degree
> > derivate a
> > > method might need without ever having to change the framework. We
> > provide
> > > the maximum amout of derivate creation support to the user while
> > providing
> > > us with the minimual amount of maintenance.
> >
> > Given that numerical algorithms try to avoid derivatives
> > for a number of reasons, a generic method to calculate
> > potentially arbitrary derivatives seems to be overkill.
> > In fact, only very few methods even use a second derivative.
>
> Just because we the creators of the framework see no need for higher
> order derivates, the end users, whose needs no one can fully
> enumerate, may have a quite reasonable use for them. Your approach
> limits them to the 2nd derivate with no easy means a creating anything
> higher. My approach doesn't limit the end user and gives them the
> flexibility to create any derivate they feel they need.
>
> >
> > Furthermore, the implementor of the function should know best
> > how the derivatives are calculated. Note that the contract
> > allows for numerical derivatives (although this is discouraged).
> > If a function object is returned, implementation is decoupled,
>
> How are they decoupled. If you create a function, you control the
> implementation for creating the derivative function object and is
> directly coupled to your object. The implementor, you, is in control
> of all details of the function and its derivate.
>
> > and higher derivatives might be implemented by someone who
> > does not match the precision and performance perceptions of
> > the implementor of the base function.
>
> Since the same implementor is in control of both the function and its
> derivate that is fault of that one implementor. Furthermore, why
> should we care? If the user wants to create a function with less or
> more precision or performance, whose to stop them? Your approach
> can't prevent it and it might be all the precision/performance they
> need.
>
> >
> > I'd like to keep it all in one object for this reasons.
> > I know this may seem to be a bit inconvenient for people who
> > want to find a zero of a derivative of a function they just
> > see lying around, but you definitely don't want to solve a
> > numerical derivative for a root, you'll just get garbage.
>
> Again, you're pigeonholing the needs of the user. The user might have
> needs we can't possibly anticipate. That's the exactly reason we
> should make creating derivatives convenient and easy for the user. We
> shouldn't take the arrogant position that we know best and prohibit
> users from doing the things they want to do.
>
> > > Why not allow the user to supply any number of initial values and
> > design the
> > > solvers to compensate as best they can when not enough values are
> > provided.
> > > Each algorithm has criteria about the initial values that must be
> > met before
> > > root finding can be attempted. If the user provided initial values
> > do not
> > > satisfy the criteria, the solver should first attempt to morph the
> > initial
> > > values into a set of values that do satisfy the criteria. If this
> > can not
> > > be accomplish, then the solver would raise an exception. This
> would
> > take
> > > away some of the user's responsibility of knowing how these
> > algorithms
> > > actually work and place it on us to create more robust components.
> >
> > The user should have some rough ideas of the interval of the root.
> > It could be quite a large interval, but then the solver may complain.
> > There is no "best" algorithm, which can take arbitrary functions
> > and start values and produce a root with a reasonable amount of
> > work.
>
> I never said anything about a best algorithm. What i did say is every
> solver should take an arbitrary function and start value(s) and make
> its best effort to return a root. If that means generating starting
> values, the solver should have the intelligence to do such a thing.
>
> > Also, many functions have multiple roots. An exhaustive search for
> > all roots is often impossible.
>
> When did I introduce finding all roots into this discussion?
>
> > Therefore, the interval is mostly to
> > prevent the solver from finding a root outside of the domain of
> > interest for the user. There is no way to get the user's domain of
> > interest without explicit input.
>
> Why can't the domain of interest and explicit input be a single value?
> That's exactly how I did the CDF stuff using a root finding method
> which needs a bracket and it worked fine. And it's a generic process
> that could be incorporated into all the solvers.
>
> Brent Worden
> http://www.brent.worden.org/
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org


It seems to me that typical usage of the package will be (in pseudocode):

findRoot( f, x0, x1 ) ;

I personally want to be able to specify just a function and a region of
interest (no guarantee it brackets a root) and ask for a root, with an
indication if it turns out no root is found. I think we would unnecessarily
complicate our first version of the interface by introducing the concept of
derivatives, which I hazard to guess most of the users of this package would
not be able to specify. Of the root-finding algorithms that I am aware of that
have been coded for commons-math so far (viz., bisection, Brent's method,
Ridders' method), none require evaluation of a function's derivatives, so why
not leave it at that for now, as Phil said? We will hear from those more
advanced users who want more explicit control -- if any such exist -- soon
enough after release of the package, and we can address their needs then.


To follow the derivative thread a step anyway, I would prefer to add a second
parameter to the derivative method, which is the order of the derivative, e.g.,

getDerivative( f, n ) ;

That way only one method call would be required on the user's part to get
whatever order derivative they want, rather than having to do n calls to
getDerivative( f ) explicitly. Also (and I don't see any immediate use for
it), such an interface would allow for non-integer orders of derivatives.



Al

=====
Albert Davidson Chou

Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
http://calendar.yahoo.com

[Phil Steitz]

Al Chou wrote:
> --- ***@worden.org wrote:
>
>>On Fri, 30 May 2003 20:14:23 +0200, "J.Pietschmann" wrote:
>>
>>>Brent Worden wrote:
>>>
>>>>I agree. The this looks like a very solid framework. One
>>>
>>>suggestion I
>>>
>>>>would like to make, is instead of a both a firstDirevative and
>>>>secondDerivate method to evaluate the derivates. Create a single
>>>>getDerivate() method that returns a UnivariateRealFunction. That
>>>
>>>way if a
>>>
>>>>user needs the n-th derivate, they just call the getDerivate()
>>>
>>>method n
>>>
>>>>times, once on the original function and once on each of the
>>>
>>returned
>>
>>>>functions. That way, common-math supports creating whatever degree
>>>
>>>derivate a
>>>
>>>>method might need without ever having to change the framework. We
>>>
>>>provide
>>>
>>>>the maximum amout of derivate creation support to the user while
>>>
>>>providing
>>>
>>>>us with the minimual amount of maintenance.
>>>
>>>Given that numerical algorithms try to avoid derivatives
>>>for a number of reasons, a generic method to calculate
>>>potentially arbitrary derivatives seems to be overkill.
>>>In fact, only very few methods even use a second derivative.
>>
>>Just because we the creators of the framework see no need for higher
>>order derivates, the end users, whose needs no one can fully
>>enumerate, may have a quite reasonable use for them. Your approach
>>limits them to the 2nd derivate with no easy means a creating anything
>>higher. My approach doesn't limit the end user and gives them the
>>flexibility to create any derivate they feel they need.
>>
>>
>>>Furthermore, the implementor of the function should know best
>>>how the derivatives are calculated. Note that the contract
>>>allows for numerical derivatives (although this is discouraged).
>>>If a function object is returned, implementation is decoupled,
>>
>>How are they decoupled. If you create a function, you control the
>>implementation for creating the derivative function object and is
>>directly coupled to your object. The implementor, you, is in control
>>of all details of the function and its derivate.
>>
>>
>>>and higher derivatives might be implemented by someone who
>>>does not match the precision and performance perceptions of
>>>the implementor of the base function.
>>
>>Since the same implementor is in control of both the function and its
>>derivate that is fault of that one implementor. Furthermore, why
>>should we care? If the user wants to create a function with less or
>>more precision or performance, whose to stop them? Your approach
>>can't prevent it and it might be all the precision/performance they
>>need.
>>
>>
>>>I'd like to keep it all in one object for this reasons.
>>>I know this may seem to be a bit inconvenient for people who
>>>want to find a zero of a derivative of a function they just
>>>see lying around, but you definitely don't want to solve a
>>>numerical derivative for a root, you'll just get garbage.
>>
>>Again, you're pigeonholing the needs of the user. The user might have
>>needs we can't possibly anticipate. That's the exactly reason we
>>should make creating derivatives convenient and easy for the user. We
>>shouldn't take the arrogant position that we know best and prohibit
>>users from doing the things they want to do.
>>
>>
>>>>Why not allow the user to supply any number of initial values and
>>>
>>>design the
>>>
>>>>solvers to compensate as best they can when not enough values are
>>>
>>>provided.
>>>
>>>>Each algorithm has criteria about the initial values that must be
>>>
>>>met before
>>>
>>>>root finding can be attempted. If the user provided initial values
>>>
>>>do not
>>>
>>>>satisfy the criteria, the solver should first attempt to morph the
>>>
>>>initial
>>>
>>>>values into a set of values that do satisfy the criteria. If this
>>>
>>>can not
>>>
>>>>be accomplish, then the solver would raise an exception. This
>>>
>>would
>>
>>>take
>>>
>>>>away some of the user's responsibility of knowing how these
>>>
>>>algorithms
>>>
>>>>actually work and place it on us to create more robust components.
>>>
>>>The user should have some rough ideas of the interval of the root.
>>>It could be quite a large interval, but then the solver may complain.
>>>There is no "best" algorithm, which can take arbitrary functions
>>>and start values and produce a root with a reasonable amount of
>>>work.
>>
>>I never said anything about a best algorithm. What i did say is every
>>solver should take an arbitrary function and start value(s) and make
>>its best effort to return a root. If that means generating starting
>>values, the solver should have the intelligence to do such a thing.
>>
>>
>>>Also, many functions have multiple roots. An exhaustive search for
>>>all roots is often impossible.
>>
>>When did I introduce finding all roots into this discussion?
>>
>>
>>>Therefore, the interval is mostly to
>>>prevent the solver from finding a root outside of the domain of
>>>interest for the user. There is no way to get the user's domain of
>>>interest without explicit input.
>>
>>Why can't the domain of interest and explicit input be a single value?
>> That's exactly how I did the CDF stuff using a root finding method
>>which needs a bracket and it worked fine. And it's a generic process
>>that could be incorporated into all the solvers.
>>
>>Brent Worden
>>http://www.brent.worden.org/
>>
>>---------------------------------------------------------------------
>>To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
>>For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>
>
>
> It seems to me that typical usage of the package will be (in pseudocode):
>
> findRoot( f, x0, x1 ) ;
>
> I personally want to be able to specify just a function and a region of
> interest (no guarantee it brackets a root) and ask for a root, with an
> indication if it turns out no root is found. I think we would unnecessarily
> complicate our first version of the interface by introducing the concept of
> derivatives, which I hazard to guess most of the users of this package would
> not be able to specify. Of the root-finding algorithms that I am aware of that
> have been coded for commons-math so far (viz., bisection, Brent's method,
> Ridders' method), none require evaluation of a function's derivatives, so why
> not leave it at that for now, as Phil said? We will hear from those more
> advanced users who want more explicit control -- if any such exist -- soon
> enough after release of the package, and we can address their needs then.
>

+1

We are really talking about two things here -- rootfinding and
representing R->R functions. The first requries a bit of the second,
but not much -- not even first derivatives if we skip Newton. I would
prefer to focus on getting the rootfinding framework in place and
postpone (or continue asynchronously ;-) what promises to be a lively
discussion of how we should be thinking about functions. Note that R->R
is just the beginning ...

FWIW, I think that it will be best for both us and our users if we try
to stay as close as possible to the applications that we are working
with and bring in mathematical abstractions only as we really need them.

>
> To follow the derivative thread a step anyway, I would prefer to add a second
> parameter to the derivative method, which is the order of the derivative, e.g.,
>
> getDerivative( f, n ) ;
>
> That way only one method call would be required on the user's part to get
> whatever order derivative they want, rather than having to do n calls to
> getDerivative( f ) explicitly. Also (and I don't see any immediate use for
> it), such an interface would allow for non-integer orders of derivatives.
>
>
>
> Al
>
> =====
> Albert Davidson Chou
>
> Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
>
> __________________________________
> Do you Yahoo!?
> Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
> http://calendar.yahoo.com
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-***@jakarta.apache.org
> For additional commands, e-mail: commons-dev-***@jakarta.apache.org
>

[J.Pietschmann]

***@worden.org wrote:
> Just because we the creators of the framework see no need for higher
> order derivates, the end users, whose needs no one can fully
> enumerate, may have a quite reasonable use for them.

Well, I designed the interface as a helper for a variety of
numerical algorithms, not for general use. It's the way the
user supplies his function to solvers. Solvers rarely need
higher derivatives. If a solver has to use a third or higher
derivative, we can still extend the interface. If someone else
wants to use an algorithm which needs higher derivatives, he
can use his own mechanism or extend what is available on his
own.

> If you create a function, you control the
> implementation for creating the derivative function object and is
> directly coupled to your object. The implementor, you,

Not quite. I expect a variety of prebuild function objects
and perhaps some mechanisms for more generic definition of
functions, but if push comes to shove, it's the user of
j-c-math who has to code a class for supplying his function
to the solver. I want to think him of the calculations to
be closely related.

[rearranged]

> Why can't the domain of interest and explicit input be a single value?
> That's exactly how I did the CDF stuff using a root finding method
> which needs a bracket and it worked fine. And it's a generic process
> that could be incorporated into all the solvers.

Nice, build a modular system: an interval(or bracket)-from-start-value
builder and solvers which need an interval (or more likely a bracket).
You can provide convenience objects which combines everything.

> When did I introduce finding all roots into this discussion?

You didn't.
A solver might not find the root closest to the start value,
with no indication that there is a root closer to the start
value than the reported root. I mean something like start
value 0, expected root near 0.5, reported root 1E15. I wager
a guess this will irritate a few people. OTOH, letting the
solver do heavy checking every time just in case may be seen
as an unnecessary overhead to more advanced users, and because
this *will* be equivalent to the problem of finding *all*
roots in general, it is often impossible.
Note that the most difficult part in root finding is not finding
the root, but detecting when you stray out of the domain of
applicability of the algorithm and detecting ill conditioned
problems.

Having the user to supply an interval is no guarantee that he
knows what he is doing, but if the interval size is much larger
than the required accuracy suggests, implementations could just
bail out.

J.Pietschmann