19 views

Skip to first unread message

Apr 22, 2010, 10:27:55 PM4/22/10

to sage-devel

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

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

Apr 22, 2010, 11:52:53 PM4/22/10

to sage-devel

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)?

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)?

Apr 24, 2010, 6:12:02 AM4/24/10

to sage-...@googlegroups.com

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
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

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

framework in sage.symbolic.integration.

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

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)?

that too.

Cheers,

Burcin

Apr 26, 2010, 3:05:33 AM4/26/10

to sage-...@googlegroups.com

Excellent Burcin - thanks!

Ill try your ideas below.

Ross

Ill try your ideas below.

Ross

May 2, 2010, 10:43:17 PM5/2/10

to sage-devel

> 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
> 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.

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

May 3, 2010, 2:02:48 PM5/3/10

to sage-...@googlegroups.com

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

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

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.

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

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):

> ...

Burcin

May 4, 2010, 11:32:38 AM5/4/10

to sage-...@googlegroups.com

Burcin

> Your example is a good test case, so please

> keep on trying, sending emails, and poking people (me) to work on this.

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

> 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
> something to experiment with?

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.

> 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. :)

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

May 4, 2010, 2:45:19 PM5/4/10

to sage-...@googlegroups.com

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

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.

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?

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

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

May 6, 2010, 8:42:55 AM5/6/10

to sage-...@googlegroups.com

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

> 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

May 6, 2010, 11:56:16 AM5/6/10

to sage-...@googlegroups.com

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

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.

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 ! :-)

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)

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

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu