[sage-devel] erf + solve

[Ross Kyprianou]

Question 1:
Is it possible (and reasonable) to have the error function, erf,
return 0 for "erf(0)"?
Currently it returns the expression: erf(0)

Question 2 (related):
The standard normal (or Gaussian) curve has half its (unit) area to
the left (and right) of x==0 as we see here...
sage: gaussian = 1/sqrt(2*pi)*exp( -(1/2)*x^2 )
sage: integrate( gaussian, x, -oo, 0)
1/2

To find the value of t for which the area is 1/2 we might try
sage: solve( integrate(gaussian, x, -oo, t)==1/2, t )

Unfortunately we get the expression
[erf(1/2*sqrt(2)*t) == 0]

which for t==0 reduces to erf(0) which ideally would reduce to 0 hence
Question 1 above ;-)
I suppose Question 2 is:
Although the erf doco suggests this is all done with PARI, is it
possible for [erf(1/2*sqrt(2)*t) == 0] to be made to reduce to
[t==0]?

--
To post to this group, send an email to sage-...@googlegroups.com
To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org

[Ross Kyprianou]

Addendum: I suppose a general query would be how do we incorporate new
knowledge into Sage (narrowing this down to things like (a) closed
form expressions of integrals (b) well known expressions of finite or
infinite sums (c) well known solutions of equations such as the
previous message

Is it a matter of identifying the area (e.g. pynac, pari, maxima) and
proposing a change to that codebase (possibly in C or Lisp)?

[Burcin Erocal]

On Thu, 22 Apr 2010 20:52:53 -0700 (PDT)
Ross Kyprianou <ros...@gmail.com> wrote:

> Addendum: I suppose a general query would be how do we incorporate new
> knowledge into Sage (narrowing this down to things like
>
> (a) closed form expressions of integrals

You should add a new integrator function and register it in the
dictionary sage.symbolic.integration.integral.available_integrators.

At some point we also need to come up with a protocol to allow these
functions to transform the input and pass it on to be processed by the
next one in the queue.

> (b) well known expressions of finite or infinite sums

This is not possible at the moment. We need to setup something like the
framework in sage.symbolic.integration.

> (c) well known solutions of equations such as the previous message

This can be done by adding an _eval_() method to
sage.functions.other.Function_erf. It might look like:

def _eval_(self, x):
if x == 0:
return 0
return BuiltinFunction._eval_default(self, x)

Note that checking if a symbolic expression is equal to 0 might be
expensive. Looking at sage/symbolic/expression.pyx, I see that there
is no way to call the is_zero() method of pynac directly from Python.
Perhaps we should allow this, and call this function if the
argument to _eval_ is a symbolic expression.

> Is it a matter of identifying the area (e.g. pynac, pari, maxima) and
> proposing a change to that codebase (possibly in C or Lisp)?

All of the above can be done in Python/Cython, if not we should fix
that too.


Cheers,
Burcin

[ross kyprianou]

Excellent Burcin - thanks!
Ill try your ideas below.
Ross

[Ross Kyprianou]

> You should add a new integrator function and register it in the
> dictionary sage.symbolic.integration.integral.available_integrators.
>
> At some point we also need to come up with a protocol to allow these
> functions to transform the input and pass it on to be processed by the
> next one in the queue.

So I guess registering means adding (to
sage.symbolic.integration.integral) the last line (shown here)

available_integrators['maxima'] = external.maxima_integrator
available_integrators['sympy'] = external.sympy_integrator
available_integrators['mathematica_free'] =
external.mma_free_integrator
available_integrators['sage'] = external.sage_integrator

and including a corresponding function to
sage.symbolic.integration.external

I think the aim would be that anybody needing a known integral that is
not caught by the other 3 integrators, patching this integral in this
"last chance" "sage" integrator which we have more control over than
the other integrators. (Alternatively, can we patch into sympy as an
alternative and get an upstream update to their integrator instead?)

Ive looked at the 1-2 year old threads in this list on integration and
I guess most of the discussion have been implemented. It seems the
schema above was implemented so we are using existing integrators
rather than opting to recreate another large internal integrator (is
that correct?). If so, it might be that the integrals made available
very quickly with patches to this new sage integrator would ideally
find themselves implemented in maxima and/or sympy eventually and
dropped from this new sage integrator to keep it compact.

Regardless of the responses to the above, is the following what I
should implement?
1) Add to sage.symbolic.integration.integral.available_integrators...
available_integrators['sage'] = external.sage_integrator

2) Include a corresponding function (to return the integral result) in
sage.symbolic.integration.external
def sage_integrator(expression, v, a=None, b=None):
...

all the best

[Burcin Erocal]

Hi Ross,

On Sun, 2 May 2010 19:43:17 -0700 (PDT)
Ross Kyprianou <ros...@gmail.com> wrote:

> > You should add a new integrator function and register it in the
> > dictionary sage.symbolic.integration.integral.available_integrators.
> >
> > At some point we also need to come up with a protocol to allow these
> > functions to transform the input and pass it on to be processed by
> > the next one in the queue.
>
> So I guess registering means adding (to
> sage.symbolic.integration.integral) the last line (shown here)
>
> available_integrators['maxima'] = external.maxima_integrator
> available_integrators['sympy'] = external.sympy_integrator
> available_integrators['mathematica_free'] =
> external.mma_free_integrator
> available_integrators['sage'] = external.sage_integrator
>
> and including a corresponding function to
> sage.symbolic.integration.external

You are on the right track. Sorry for not being more explicit before.

Note that most of this is new and we are just developing the necessary
framework as we go along. Your example is a good test case, so please
keep on trying, sending emails, and poking people (me) to work on this.

I suggest giving a more specific name to your function. Assuming that
you're working on the erf function, erf_integrator() might be a good
name.

You should definitely add your function to the list of available
integrators. This will let people use only your function if they want
to, by using the algorithm keyword argument of the top level integrate
function.

The new function you implement doesn't need to go in the file
sage/symbolic/integration/external.py. That file contains the
routines which call external systems for integration. I suggest putting
a new file in sage/symbolic/integration, or adding it somewhere close
to the definition of the erf() function.

Later we can also consider letting functions automatically register
custom integration routines, similar to the way they register custom
numeric evaluation or differentiation methods. I don't really like this
idea, since the order we call these functions might be important.


Can you post some example code (your integrator function) so I have
something to experiment with?

> I think the aim would be that anybody needing a known integral that is
> not caught by the other 3 integrators, patching this integral in this
> "last chance" "sage" integrator which we have more control over than
> the other integrators. (Alternatively, can we patch into sympy as an
> alternative and get an upstream update to their integrator instead?)
>
> Ive looked at the 1-2 year old threads in this list on integration and
> I guess most of the discussion have been implemented. It seems the
> schema above was implemented so we are using existing integrators
> rather than opting to recreate another large internal integrator (is
> that correct?). If so, it might be that the integrals made available
> very quickly with patches to this new sage integrator would ideally
> find themselves implemented in maxima and/or sympy eventually and
> dropped from this new sage integrator to keep it compact.

We are definitely aiming to implement symbolic integration natively in
Sage. There are some existing implementations even. It's just a matter
of getting things cleaned up and submitted to Sage. Since the work
is done for research, it can be hard to get them ready for public
consumption. :)

> Regardless of the responses to the above, is the following what I
> should implement?
> 1) Add to sage.symbolic.integration.integral.available_integrators...
> available_integrators['sage'] = external.sage_integrator

You should also add it to the lists on line 60 and 146 of
sage/symbolic/integration/integral.py.

> 2) Include a corresponding function (to return the integral result) in
> sage.symbolic.integration.external
> def sage_integrator(expression, v, a=None, b=None):
> ...

Thanks.

Burcin

[ross kyprianou]

Burcin

> Your example is a good test case, so please
> keep on trying, sending emails, and poking people (me) to work on this.

> Can you post some example code (your integrator function) so I have
> something to experiment with?

Id like to do as much as possible. This might be a good example for me
to start in development.
Ill can start with specs (the integral formulas Ive listed) and try to
express these in code and send you that first.
(And to keep this reply short, Ill just state I acknowledge all your
points above about the implementation details)

> the order we call these functions might be important.

Im sure this is important

> We are definitely aiming to implement symbolic integration natively in Sage.

Wow. Sounds like a huge project. Same or bigger scope than Maxima integration?

>There are some existing implementations even. It's just a matter
> of getting things cleaned up and submitted to Sage. Since the work
> is done for research, it can be hard to get them ready for public
> consumption. :)

Integration has been discussed on this list but havent seen any detail
of this new integration system
Im very interested because the two areas I need to get most into (for
my thesis) is
1) integration (particularly of functions associated with common
probability distributions and their multivariate counterparts) and
2) signal processing (which Im looking out for open source libraries -
I know about scipy but Im sure there may be more out there to collate
and bring together consistently)
Id like to think Ill be doing sage development in integration and
signal processing before long

all the best
Ross

[Burcin Erocal]

Hi Ross,

On Wed, 5 May 2010 01:02:38 +0930
ross kyprianou <ros...@gmail.com> wrote:

> > Your example is a good test case, so please
> > keep on trying, sending emails, and poking people (me) to work on
> > this. Can you post some example code (your integrator function) so
> > I have something to experiment with?
>
> Id like to do as much as possible. This might be a good example for me
> to start in development.
> Ill can start with specs (the integral formulas Ive listed) and try to
> express these in code and send you that first.

A first step might be to put your examples on the wiki. The wiki
supports jsmath, so you can just write the integrals and the expected
answers in a table. Here is a possible address:

http://wiki.sagemath.org/symbolics/integration/examples

You can just click on the link and start editing. :)

> > We are definitely aiming to implement symbolic integration natively
> > in Sage.
> Wow. Sounds like a huge project. Same or bigger scope than Maxima
> integration?

There is no real "project." The purpose of switching to pynac for the
symbolics backend was to enable interested people to implement basic
symbolics functionality in Sage. IMHO, ruling out the possibility of
having native integration code in Sage is not a good long term strategy.

> >There are some existing implementations even. It's just a matter
> > of getting things cleaned up and submitted to Sage. Since the work
> > is done for research, it can be hard to get them ready for public
> > consumption. :)
> Integration has been discussed on this list but havent seen any detail
> of this new integration system
> Im very interested because the two areas I need to get most into (for
> my thesis) is
> 1) integration (particularly of functions associated with common
> probability distributions and their multivariate counterparts) and

It is most likely that the examples you're interested in will not be
covered by the implementations that come out of mathematical research.
AFAICT, your best bet is to prepare tables of integrals you want
handled, or point to already existing tables in books or elsewhere, and
work towards a solution based on pattern matching.

Once you put up a few examples, perhaps we can see if there is any
other algorithmic solution to your problems.


Thanks.

Burcin

[ross kyprianou]

Burcin

> BTW, you could also try to make a patch for the erf() function, to
> add the _eval_() method I mentioned in this thread on sage-devel.
> Don't worry about details I wrote in
> that email, just put an _eval_ method that returns 0 when the input is 0

Heres the patch. I hope I created it properly. I used Mercurial Queues.
Assuming its ok, as my first patch, it was an ideal for learning about
creating patches.
Well done to the people that created [1]. It was very clear and meant
I did this without the usual newbie questions.
(So that documentation was an excellent investment in time Guys ! :-)

(And let me know if the patch can be improved (even if you can make
the improvement yourself quickly)
Everything Im doing at this stage is worthwhile experience
e.g. If you think its appropriate for me to create a ticket, to go
through that exercise, Id be happy to)

I know youre busy so Ill work on the integrals until I here from you
cheers
Ross

[1] Walking Through the Development Process
http://www.sagemath.org/doc/developer/walk_through.html#chapter-walk-through

[Burcin Erocal]

Hi Ross,

On Thu, 6 May 2010 22:12:55 +0930
ross kyprianou <ros...@gmail.com> wrote:

> > BTW, you could also try to make a patch for the erf() function, to
> > add the _eval_() method I mentioned in this thread on sage-devel.
> > Don't worry about details I wrote in
> > that email, just put an _eval_ method that returns 0 when the input
> > is 0
>
> Heres the patch. I hope I created it properly. I used Mercurial
> Queues. Assuming its ok, as my first patch, it was an ideal for
> learning about creating patches.

This looks great! Now we have a brave new developer stepping into the
world of symbolics. :)

> Well done to the people that created [1]. It was very clear and meant
> I did this without the usual newbie questions.
> (So that documentation was an excellent investment in time Guys ! :-)

I agree. Thanks especially to Minh and Rob Beezer for the walk through.
Apologies if I'm forgetting anyone who contributed.

> (And let me know if the patch can be improved (even if you can make
> the improvement yourself quickly)
> Everything Im doing at this stage is worthwhile experience
> e.g. If you think its appropriate for me to create a ticket, to go
> through that exercise, Id be happy to)

Yes, please create a ticket in the symbolics component and attach your
patch there.

A few short comments:

* The first line of the class Function_erf (line 27 in
sage/functions/other.py) is

_eval_ = BuiltinFunction._eval_default

Since you're redefining the function _eval_(), you can remove this
line.

* You should add doctests to check if your function returns the right
type of result for different representations of 0. For example, what
happens when you evaluate erf(0.), erf(CC(0)), etc.

* Can you time the new erf function on float, complex, CC, RR, inputs
to see if it slows down? I.e.,

%timeit t = erf(CC.1)
%timeit t = erf(RR(1))


Thanks a lot for your efforts. Keep up the good work!

Burcin