Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Symbolic-Numeric Integration

56 views
Skip to first unread message

Peter Luschny

unread,
Feb 7, 2022, 7:04:28 AM2/7/22
to
That might interest some here:

Symbolic-Numeric Integration of Univariate Expressions based on Sparse Regression https://arxiv.org/abs/2201.12468

nob...@nowhere.invalid

unread,
Feb 11, 2022, 5:30:24 AM2/11/22
to

Peter Luschny schrieb:
>
> That might interest some here:
>
> Symbolic-Numeric Integration of Univariate Expressions based on
> Sparse Regression https://arxiv.org/abs/2201.12468
>

The total numbers of indefinite integrals in Table 1 on page 7 do not
agree with the sources; in my files, Hearn's suite comprises 285
integration problems and Timofeev's suite 705 problems. From these
numbers, one computes a modest success rate of 67% for Hearn, and a
disappointingly low 19% for Timofeev. The principle used in selecting
the integration problems must be clearly documented in the paper -
readers will suspect Texas sharpshooting otherwise.

Also, two instances of "Apostle" on page 7 should be corrected to
"Apostol".

Martin.

PS: Are the authors reinventing ancient technology? Their integrator
should perhaps be run against Moses' SIN program.

anti...@math.uni.wroc.pl

unread,
Feb 21, 2022, 4:50:49 PM2/21/22
to
Their results are clearly unimpressive. Concerning "ancient",
their idea is closest to parallel integration, also called
Risch-Norman method. There are two main differences:
- they use floating point linear solver instead of rational one
- they use only heuristics instead of structure theorem (+ heuristic)
to build space of candidate integrals.

IIUC correctly they use floating point because they do not have
rational solver. So maybe their method should be called "very
poor man integrator".

--
Waldek Hebisch

Nasser M. Abbasi

unread,
Sep 3, 2022, 8:56:59 PM9/3/22
to
Fyi,

I was trying to add this integrator to CAS integration tests. But after some testing,
I found big problem and not able to continue. I did submit a ticket on it.

https://github.com/JuliaSymbolics/Symbolics.jl/issues/730

Hopefully when this is resolved, will be able to add Julia integrator.

I get different anti-derivative each time I make a call to integrate so I am not sure
how to handle this and which result to accept. One call can fail to integrate, then
the second call it can integrate it, or give different result and so on.

Examples of this issue are at the above top link.

From

https://symbolicnumericintegration.sciml.ai/dev/

it says "uses a randomized algorithm" so I think this has something to do with it.

But this means my result will not match someone else result if I post them.

If there is a way to fix the random seem to some value, so same anti-derivative is
generated each time, this I can do that, but I do not know how to do this.


--Nasser

Nasser M. Abbasi

unread,
Sep 3, 2022, 9:10:16 PM9/3/22
to
On 9/3/2022 7:56 PM, Nasser M. Abbasi wrote:
> If there is a way to fix the random seem to some value, so same anti-derivative is
> generated each time, this I can do that, but I do not know how to do this.

Typo above:

Meant to fix the random "seed" used before making the call.

This way I post the results, and then say this was generated using this specific seed.

So result shown will be reproducable. Otherwise now, Julia gets
different score each time it is run.

--Nasser


Nasser M. Abbasi

unread,
Sep 3, 2022, 10:00:53 PM9/3/22
to
fyi

I think I found a fix. But doing

Random.seed!(12) #must do this in order to reproduce same result

Then now same antiderivative is generated each time. So now I can add
Julia integrator to CAS integration tests.

--Nasser

Jeff Barnett

unread,
Sep 4, 2022, 12:25:30 AM9/4/22
to
That doesn't seem correct to me since users will get whatever first seed
the system generates (unless they choose their own as you are). The
chances that they will choose 12 as you have is small. I would either
give the system a fail because it is unreliable on that problem or I
would test it for, say, a 100 or so seeds and give the number (a
percent) that it gets correct. The fact that you had to search for a
seed that led to correct behavior voids the system's ability for that
example. It might be that you could find a seed for each problem in the
test suite that leads to correct behavior for the problem. But that
doesn't make the system perfect.

Basically when you are testing a system that uses a "probabilistic"
approach, you need to develop scoring methods that account for that fact.
--
Jeff Barnett

Nasser M. Abbasi

unread,
Sep 4, 2022, 12:41:54 AM9/4/22
to
On 9/3/2022 11:25 PM, Jeff Barnett wrote:

>>
>> Random.seed!(12) #must do this in order to reproduce same result
>>
>> Then now same antiderivative is generated each time. So now I can add
>> Julia integrator to CAS integration tests.
>

> That doesn't seem correct to me since users will get whatever first seed
> the system generates (unless they choose their own as you are). The
> chances that they will choose 12 as you have is small. I would either
> give the system a fail because it is unreliable on that problem or I
> would test it for, say, a 100 or so seeds and give the number (a
> percent) that it gets correct. The fact that you had to search for a
> seed that led to correct behavior voids the system's ability for that
> example. It might be that you could find a seed for each problem in the
> test suite that leads to correct behavior for the problem. But that
> doesn't make the system perfect.
>
> Basically when you are testing a system that uses a "probabilistic"
> approach, you need to develop scoring methods that account for that fact.

But I have documented that this result is based on using seed `12`.
Any other value could very well generate new result.

So if someone wants to obtain the same result as I have obtained,
they have to use same seed.

For an example

-----------------------
using Symbolics
using SymbolicNumericIntegration
using Dates
import Random

Random.seed!(12) #must do this in order to reproduce same result

@syms x

println("before calling integrate")

time1 = now();
anti,n_unsolved,residual_error =integrate(1 / (((1 - 2x)^3)*((2 + 3x)^3)*((3 + 5x)^2)),x)
time2 = now();
time_used_in_sec = (Dates.DateTime(time2)-Dates.DateTime(time1))/ Millisecond(1) * (1 / 1000)

println("time used =",time_used_in_sec)

println("anti = ",anti)
println("n_unsolved = ",n_unsolved)
println("residual_error = ",residual_error)
--------------------------------

Now calling this many times, same result is obtained:

---------------
>julia bug_2.jl
before calling integrate
time used =70.639
anti = 0
n_unsolved = 1.6602780101345385e13((x - (1//2))^-1) + 2.8289042577780735e15(((3//5) + x)^-1) + 8.301390050672688e12(((1//4) + x^2 - x)^-1) - 1.2230893145047396e6(((4//9) + (4//3)*x + x^2)^-1) - 1.697342553676207e15(((9//25) + (6//5)*x + x^2)^-1) - 1.7402011587196218e6x*(((4//9) + (4//3)*x + x^2)^-1) - 1.6602780101345375e13x*(((1//4) + x^2 - x)^-1) - 2.82890425603774e15x*(((9//25) + (6//5)*x + x^2)^-1)
residual_error = Inf

-----------------
>julia bug_2.jl
before calling integrate
time used =71.749
anti = 0
n_unsolved = 1.6602780101345385e13((x - (1//2))^-1) + 2.8289042577780735e15(((3//5) + x)^-1) + 8.301390050672688e12(((1//4) + x^2 - x)^-1) - 1.2230893145047396e6(((4//9) + (4//3)*x + x^2)^-1) - 1.697342553676207e15(((9//25) + (6//5)*x + x^2)^-1) - 1.7402011587196218e6x*(((4//9) + (4//3)*x + x^2)^-1) - 1.6602780101345375e13x*(((1//4) + x^2 - x)^-1) - 2.82890425603774e15x*(((9//25) + (6//5)*x + x^2)^-1)
residual_error = Inf

---------------
>julia bug_2.jl
before calling integrate
time used =69.586
anti = 0
n_unsolved = 1.6602780101345385e13((x - (1//2))^-1) + 2.8289042577780735e15(((3//5) + x)^-1) + 8.301390050672688e12(((1//4) + x^2 - x)^-1) - 1.2230893145047396e6(((4//9) + (4//3)*x + x^2)^-1) - 1.697342553676207e15(((9//25) + (6//5)*x + x^2)^-1) - 1.7402011587196218e6x*(((4//9) + (4//3)*x + x^2)^-1) - 1.6602780101345375e13x*(((1//4) + x^2 - x)^-1) - 2.82890425603774e15x*(((9//25) + (6//5)*x + x^2)^-1)
residual_error = Inf


Without fixing the seed, different result would have obtained.

As long as same result can be obtained, I think this is valid approach?

> Basically when you are testing a system that uses a "probabilistic"
> approach, you need to develop scoring methods that account for that fact.

I think this is beyond the scope of what I am doing now and I do not know
how to design such a test myself at this time.

May be someone else could do a test based on this approach you are
suggesting. It sounds interesting.

Regards
--Nasser

Nasser M. Abbasi

unread,
Sep 4, 2022, 12:43:29 AM9/4/22
to
On 9/3/2022 11:25 PM, Jeff Barnett wrote:

>>
>> Random.seed!(12) #must do this in order to reproduce same result
>>
>> Then now same antiderivative is generated each time. So now I can add
>> Julia integrator to CAS integration tests.
>


> That doesn't seem correct to me since users will get whatever first seed
> the system generates (unless they choose their own as you are). The
> chances that they will choose 12 as you have is small. I would either
> give the system a fail because it is unreliable on that problem or I
> would test it for, say, a 100 or so seeds and give the number (a
> percent) that it gets correct. The fact that you had to search for a
> seed that led to correct behavior voids the system's ability for that
> example. It might be that you could find a seed for each problem in the
> test suite that leads to correct behavior for the problem. But that
> doesn't make the system perfect.
>
> Basically when you are testing a system that uses a "probabilistic"
> approach, you need to develop scoring methods that account for that fact.


But I have documented that this result is based on using seed `12`.
Any other value could very well generate new result.

So if someone wants to obtain the same result as I have obtained,
they have to use same seed.

For an example

-----------------------
using Symbolics
using SymbolicNumericIntegration
using Dates
import Random

Random.seed!(12) #must do this in order to reproduce same result

> Basically when you are testing a system that uses a "probabilistic"
> approach, you need to develop scoring methods that account for that fact.

Nasser M. Abbasi

unread,
Sep 4, 2022, 1:24:26 AM9/4/22
to
On 9/3/2022 11:25 PM, Jeff Barnett wrote:

Hello again:

> I would either
> give the system a fail because it is unreliable on that problem

Sorry, I do not understand the above. Suppose it returns some
antiderivative on one call. How is the program supposed to be it
failed or not?

The program has no verification on the anti-derivative that it satisfies
the integral. It just checkes the result and looks to see if the number of
terms that were not integrated is zero or not to decide if it was
able to integrate it fully or not.

> or I
> would test it for, say, a 100 or so seeds and give the number (a
> percent) that it gets correct.


So you are saying to run the same integral 100 times, each with different seed
(say 1,2,3,4....,100 ?) saving all the results, then do histogram and
pick the anti-derivative that shows up the most frequent?

And if they are all different, what to do? Say it failed or pick one by random?

Sure, I could do that, but the time this will take will be prohibitive.

I have 12,500 integrals (picked out of 85,500 integrals, since Julia can
only do univariate ones with constant coefficients).

Lets say each call to integrate takes 25 seconds on average (I do
not have the average for Julia yet as the tests are still running, it
could be more or it could be much more) but looking at the terminal
now, some are more and some less. So have to wait to find
out what the average is.

Using 25 seconds per on call, this will result in 31,250,000 seconds or
about one year to finish just the Julia test. So it is not practical
for me to do such a setup. Doing parallel processing is beyond what
I can do on my PC.

So I think using one seed, as long as it is documented what it is, and
the result is reproducible can be considered as a compromise.

This is until Julia obtains an standard symbolic integrator. I read
somewhere that they are working on doing that for the future.

The fact that you had to search for a
> seed that led to correct behavior voids the system's ability for that
> example. It might be that you could find a seed for each problem in the
> test suite that leads to correct behavior for the problem. But that
> doesn't make the system perfect.
>
> Basically when you are testing a system that uses a "probabilistic"
> approach, you need to develop scoring methods that account for that fact.

--Nasser

Jeff Barnett

unread,
Sep 4, 2022, 1:32:46 AM9/4/22
to
On 9/3/2022 10:43 PM, Nasser M. Abbasi wrote:
> On 9/3/2022 11:25 PM, Jeff Barnett wrote:
>
>>>
>>> Random.seed!(12) #must do this in order to reproduce same result
>>>
>>> Then now same antiderivative is generated each time. So now I can add
>>> Julia integrator to CAS integration tests.
>>
>
>
>> That doesn't seem correct to me since users will get whatever first seed
>> the system generates (unless they choose their own as you are). The
>> chances that they will choose 12 as you have is small. I would either
>> give the system a fail because it is unreliable on that problem or I
>> would test it for, say, a 100 or so seeds and give the number (a
>> percent) that it gets correct. The fact that you had to search for a
>> seed that led to correct behavior voids the system's ability for that
>> example. It might be that you could find a seed for each problem in the
>> test suite that leads to correct behavior for the problem. But that
>> doesn't make the system perfect.
>>
>> Basically when you are testing a system that uses a "probabilistic"
>> approach, you need to develop scoring methods that account for that fact.
>
>
> But I have documented that this result is based on using seed `12`.
> Any other value could very well generate new result.
>
> So if someone wants to obtain the same result as I have obtained,
> they have to use same seed.

I've snipped everything below since I wasn't clear and I think there is
a misunderstanding. I assumed, tacitly, that people look at the results
of running test suites "to understand how well they can expect them to
work" they do not look at the results to learn how to solve exactly the
problems in that suite. What about a similar problem which is properly
solved with initial seed 32457689912 but fails with 12? They look at
your work to learn which systems are trustworthy and worth the price.

A zillion years ago I wrote a program to form fault trees from systems
where some of the components had circular dependencies. That fault tree
then was reduced to minimal fault sets. This gibberish meant find sets
of components such that if all failed, the system failed and that no
subset of fault set led to system failure. This is yet another
processing of and or expressions to get in a particular form and, of
course, this is one of those BP complete problems. The algorithm I
developed order all the variables and did some magic that allowed the
output tree to share, multiple times, most of its nodes dramatically
saving space. The exact amount of space but more importantly, the
runtime depended on the variable order chosen. There was no way to
guarantee a polynomial run time since it was an NP complete problem. For
engineering work I would try an ordering and after a while I'd start it
with another and I usually found solutions in very good time. The
problem was that one day it solved a problem in few seconds and the next
day it might run for hours and not terminate. Most of the work I did was
trying to develop heuristics to find good orders that led to good
performance. However, those heuristics often depended on what order they
examined parameters.

The problem with a test suite for the above is the users didn't care
about all these internal machinations; they just wanted to know if the
system would work. I could have published test cases and a good ordering
heuristic for each one but that would be useless information to an
engineer that was trying to do a reliability forecast. The work was for
the design of ultra performance aircraft and was to evaluate the
reliability of the plane and its safety and effectiveness. The joke in
all of this was that their current methods for dealing with circular
dependency were not correct; what I mean is you have in effect several
boolean equations in many variables and you substitute one of the
equations into another and an inconsistency pops out.

Let me finish by repeating a line from above: "What about a similar
problem which is properly solved with initial seed 32457689912 but fails
with 12?" Test suites must tell someone how well things work; They are
not to teach hacking.
--
Jeff Barnett

Jeff Barnett

unread,
Sep 4, 2022, 1:36:14 AM9/4/22
to
Hi again. I wrote a lengthy reply above. I'll refer it to you now.
--
Jeff Barnett

Nasser M. Abbasi

unread,
Sep 4, 2022, 3:21:25 AM9/4/22
to
On 9/3/2022 11:25 PM, Jeff Barnett wrote:
>I would either
> give the system a fail because it is unreliable on that problem or I
> would test it for, say, a 100 or so seeds and give the number (a
> percent) that it gets correct.

I've tried your idea. I've run this integral 100 times, using
100 different seeds and printed the anti-derivative each time

-------------------------
using Symbolics
using SymbolicNumericIntegration
using Dates
import Random
@syms x
integrand = 1 / (((1 - 2*x)^3)*((2 + 3*x)^3)*((3 + 5*x)^2))
for i=1:100
Random.seed!(i)
anti,n_unsolved,residual_error =integrate(integrand,x)
println("seed = ",i,", anti=",anti)
end
----------------------------

This is what I get

---------------------------
seed = 1, anti=0
seed = 2, anti=0

seed = 3, anti=0.0032245661147886633(((1//4) + x^2 - x)^-1) + 1.3932951722611113(((4//9) + (4//3)*x + x^2)^-1) + 4.2331946787976005log((3//5) + x) + 2.3517748215530303(x^3)*(((9//25) + (6//5)*x + x^2)^-1) + 3.0398064936094906x*(((4//9) + (4//3)*x + x^2)^-1) - 0.21458141730725133x - 0.9634428771843939(((9//25) + (6//5)*x + x^2)^-1) - 0.010517325162243041log(x - (1//2)) - 4.274386808487491log((2//3) + x) - 0.00644913222957772x*(((1//4) + x^2 - x)^-1) - 2.1371934042401866(x^3)*(((4//9) + (4//3)*x + x^2)^-1) - 2.452377064399745x*(((9//25) + (6//5)*x + x^2)^-1)

seed = 4, anti=-0.016330430066172893log(x - (1//2))
seed = 5, anti=0
Error from sparse_fitLinearAlgebra.SingularException(4)

seed = 6, anti=0
seed = 7, anti=0
seed = 8, anti=0
seed = 9, anti=0

seed = 10, anti=0.059870841012718994x - 0.014879657188586044(((9//25) + (6//5)*x + x^2)^-1) - 0.013484342743641575log(x - (1//2)) - 0.1077675138225977log((3//5) + x) - 0.0032459258832005468x*(((9//25) + (6//5)*x + x^2)^-1) - 0.059870841012343586(x^3)*(((9//25) + (6//5)*x + x^2)^-1)

seed = 11, anti=0
seed = 12, anti=0
seed = 13, anti=0
seed = 14, anti=0
seed = 15, anti=0
seed = 16, anti=0
seed = 17, anti=0
seed = 18, anti=-0.0026325376637571993log((3//5) + x)
seed = 19, anti=0
seed = 20, anti=0
seed = 21, anti=0
seed = 22, anti=0
seed = 23, anti=0
seed = 24, anti=0
seed = 25, anti=0

seed = 26, anti=1.719315621107581x + 2.349316120556791log((3//5) + x) + 2.578973431661372log(x - (1//2)) + 2.406019727054061x*(((1//4) + x^2 - x)^-1) - 0.9880954108885848(((1//4) + x^2 - x)^-1) - 1.7193156211075782(x^3)*(((1//4) + x^2 - x)^-1)

seed = 27, anti=0
seed = 28, anti=-67.51266553547356log((3//5) + x)
seed = 29, anti=0
seed = 30, anti=0
seed = 31, anti=0
seed = 32, anti=0
seed = 33, anti=0
Error from sparse_fitLinearAlgebra.SingularException(4)

seed = 34, anti=0
seed = 35, anti=0
seed = 36, anti=0
seed = 37, anti=0
seed = 38, anti=0
seed = 39, anti=0
seed = 40, anti=0.06317647317345293log((3//5) + x)
seed = 41, anti=0
seed = 42, anti=0
seed = 43, anti=0
seed = 44, anti=0

seed = 45, anti=(0.014378665475956662 - 0.0645731363698366im)*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1) + (-0.0995313498923154 + 0.07922575111068406im)*x*((72 + 132x + 2262(x^4) + 6743(x^5) - 754(x^2) - 234(x^6) - 1641(x^3) - 9180(x^7) - 5400(x^8))^-1)
Error from sparse_fitLinearAlgebra.SingularException(4)
Error from sparse_fitLinearAlgebra.SingularException(3)

seed = 46, anti=0
seed = 47, anti=0
seed = 48, anti=0.06276392360387156log((3//5) + x)
seed = 49, anti=-133.00981698474402log((3//5) + x)
seed = 50, anti=0
Error from sparse_fitLinearAlgebra.SingularException(4)

seed = 51, anti=0

seed = 52, anti=0.10587351232085687(((3//5) + x)^-1) + 1.3864505792612496log((3//5) + x) - 0.002356341847725731log(x - (1//2)) - 1.345612601895796log((2//3) + x)

seed = 53, anti=0
seed = 54, anti=24.16357143891812log((3//5) + x)
seed = 55, anti=0
Error from sparse_fitLinearAlgebra.SingularException(4)
seed = 56, anti=0
seed = 57, anti=0

seed = 58, anti=0.20123985741641454x + 0.003878647339941809x*(((4//9) + (4//3)*x + x^2)^-1) - 0.05704085952566808(((4//9) + (4//3)*x + x^2)^-1) - 0.40247971483032596log((2//3) + x) - 0.20123985741402833(x^3)*(((4//9) + (4//3)*x + x^2)^-1)

seed = 59, anti=0.006561881810512138x - 0.007225855430077368(((9//25) + (6//5)*x + x^2)^-1) - 0.011811387258923282log((3//5) + x) - 0.006561881810514191(x^3)*(((9//25) + (6//5)*x + x^2)^-1) - 0.009680814931677178x*(((9//25) + (6//5)*x + x^2)^-1)

seed = 60, anti=0
seed = 61, anti=0
seed = 62, anti=0
seed = 63, anti=0
seed = 64, anti=0
seed = 65, anti=0
seed = 66, anti=0
seed = 67, anti=0
seed = 68, anti=0
seed = 69, anti=0
seed = 70, anti=-771.818391946339log(x - (1//2))
seed = 71, anti=0
seed = 72, anti=2.683289125488832log((3//5) + x)
seed = 73, anti=0
seed = 74, anti=0
seed = 75, anti=0
seed = 76, anti=-0.009088132845265414log(x - (1//2))
seed = 77, anti=0
Error from sparse_fitLinearAlgebra.SingularException(3)
seed = 78, anti=0
seed = 79, anti=0
Error from sparse_fitLinearAlgebra.SingularException(4)
seed = 80, anti=0
seed = 81, anti=-0.0021911954170486065((x - (1//2))^-1) - 0.006597505528434938log(x - (1//2))
seed = 82, anti=0

seed = 83, anti=0.00853042316622113(((1//4) + x^2 - x)^-1) + 0.1222253941760599x + 0.12923006955074817(((9//25) + (6//5)*x + x^2)^-1) + 0.2713198455831079x*(((9//25) + (6//5)*x + x^2)^-1) + 0.0331534845235633(x^3)*(((1//4) + x^2 - x)^-1) - 0.049730226785344944log(x - (1//2)) - 0.2796819816593179log((3//5) + x) - 0.02534921746333307x*(((1//4) + x^2 - x)^-1) - 0.15537887869961337(x^3)*(((9//25) + (6//5)*x + x^2)^-1)

seed = 84, anti=0
seed = 85, anti=-0.006923734706044644log((3//5) + x)
seed = 86, anti=0.009191532695211821log((3//5) + x)
seed = 87, anti=0.04867187748837544log((3//5) + x)
seed = 88, anti=0
seed = 89, anti=0
seed = 90, anti=0

seed = 91, anti=0.005796247487652043x - 0.0055676726126699(((9//25) + (6//5)*x + x^2)^-1) - 0.01043324547777366log((3//5) + x) - 0.005796247487652001(x^3)*(((9//25) + (6//5)*x + x^2)^-1) - 0.0071928052588951085x*(((9//25) + (6//5)*x + x^2)^-1)
seed = 92, anti=0.025581076356940465(((9//25) + (6//5)*x + x^2)^-1) + 0.2484241701152666log((3//5) + x) + 0.13801342784197507(x^3)*(((9//25) + (6//5)*x + x^2)^-1) - 0.1380134278416556x - 0.007049706761546733x*(((9//25) + (6//5)*x + x^2)^-1)
seed = 93, anti=0

seed = 94, anti=0.03466331592299737x - 0.05509270729327145(((4//9) + (4//3)*x + x^2)^-1) - 0.06932663184599135log((2//3) + x) - 0.03466331592299913(x^3)*(((4//9) + (4//3)*x + x^2)^-1) - 0.06723314275190788x*(((4//9) + (4//3)*x + x^2)^-1)
seed = 95, anti=0
seed = 96, anti=-105277.6007622171log(x - (1//2))
seed = 97, anti=0
seed = 98, anti=0
seed = 99, anti=0

seed = 100, anti=0.15059932742674897x + 0.22589899114011922log(x - (1//2)) + 0.12969704840228988x*(((1//4) + x^2 - x)^-1) - 0.04602360827280188(((1//4) + x^2 - x)^-1) - 0.1505993274267452(x^3)*(((1//4) + x^2 - x)^-1)
-------------------------------

Looking at the above, each result which is not zero is different.

anti=0 shows up most often. This implies that we should pick anti=0.

Which implies Julia failed to integrate it, since 0 can't be antiderivative
unless the integrand is zero.

Does this sound a correct approach and what you had in mind? I am
still not sure I could do the above for all 12,500 integrals due to long
time it will take. May be I could do it for say 20 different seeds?

--Nasser







Nasser M. Abbasi

unread,
Sep 4, 2022, 7:47:31 AM9/4/22
to
Ok, I've decided to use the same settings as Julia's own SymbolicNumericIntegration.jl
test script.

I've just increased the number of trials from 10 to 20.

https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/docs/src/index.md

The following is the prototype I will be using on the 12,500 integrals. No seed
setting is used at all in this. Just like with Julia's own test script at

https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/test/runtests.jl

All the other options are the same as Julia's test script.

----------

using Symbolics
using SymbolicNumericIntegration
@syms x
integrand1 = 1 / (((1 - 2*x)^3)*((2 + 3*x)^3)*((3 + 5*x)^2))

anti,n_unsolved,residual_error =integrate( integrand1, x ;
symbolic = true,
verbose = false,
homotopy = true,
num_steps = 2,
num_trials = 20)
--------------------

This runs fast also, taking 70 seconds on my PC for 20 trials. It
failes to integrate it, giving anti=0.

The description of options are from the above link:
=========================================
integrate has the form integrate(y; kw...) or integrate(y, x; kw...),
where y is the integrand and the optional x is the variable of integration.
The keyword parameters are:

abstol (default 1e-6): the error tolerance to accept a solution.

symbolic (default true): if true, pure symbolic integration is attempted first.

bypass (default false): if true, the whole expression is considered
at once and not per term.

num_steps (default 2): one plus the number of expanded basis to check
(if num_steps is 1, only the main basis is checked).

num_trials (default 5): the number of attempts to solve the integration
numerically for each basis set.

show_basis (default false): print the basis set, useful for debugging.
Only works if verbose is also set.

homotopy (default: true as of version 0.7.0): uses the continuous
Homotopy operators to generate the integration candidates.

verbose (default false): if true, prints extra (and voluminous!)
debugging information.

radius (default 1.0): the starting radius to generate random test points.

opt (default STLSQ(exp.(-10:1:0))): the optimizer passed to sparse_regression!.

max_basis (default 110): the maximum number of expression in the basis.

complex_plane (default true): random test points are generated on the
complex plane (only over the real axis if complex_plane is
false).
-------------------------

This integrator is different than the ones I've used before, so it took sometime
getting used to how it works and how to test it.

So the the above is what I will be using, unless someone has objection or
suggestion to change something.

Final result should be ready in 1-2 weeks. All other 8 CAS systems have
already finished the univariate 12,500 integrals. Only one to do
now is Julia using the above setup and then work on generating
the reports and statistics and the web pages.


--Nasser

nob...@nowhere.invalid

unread,
Sep 4, 2022, 10:22:40 AM9/4/22
to

"Nasser M. Abbasi" schrieb:
> [...]
>

I think you are wasting your time.

This rational integral is solved instantaneously on Derive 6.10:

INT(1/((1 - 2*x)^3*(2 + 3*x)^3*(3 + 5*x)^2), x) =

- 290625*LN(5*x + 3)/14641 + 333639*LN(3*x + 2)/16807
- 274224*LN(2*x - 1)/246071287
- (1523948040*x^4 + 458007084*x^3 - 957482214*x^2 - 147486147*x
+ 160532983)
/(6391462*(2*x - 1)^2*(3*x + 2)^2*(5*x + 3))

The fourth antiderivative differentiates back to:

816521503/(25000000000*(1 - 2*x))

which obviously has very little to do with the integrand.

The first result differentiates back to:

(349038822330000*x^8 + 579297885968250*x^7 - 53181422066825*x^6
- 508545473526155*x^5 - 208696958441257*x^4 + 54883553843119*x^3
+ 45917333548640*x^2 + 8977272033166*x + 1048954387812)
/(250000000000*(1 - 2*x)^3*(3*x + 2)^3*(5*x + 3)^3) - 3/10000000000

So the integrator here is claiming that, on the scale of the integrand,
3/10000000000 vanishes, which can be accepted, and that:

(349038822330000*x^8 + 579297885968250*x^7 - 53181422066825*x^6
- 508545473526155*x^5 - 208696958441257*x^4 + 54883553843119*x^3
+ 45917333548640*x^2 + 8977272033166*x + 1048954387812)
/(250000000000*(5*x + 3))

equals unity over some interval of interest; over the range -1/2 < x <
+1/2, it actually wobbles between 0 and +6, and it already exceeds +100
at x = +1. Ugh.

Your statistics on the results of this integrator would therefore be a
case of garbage in - garbage out, unless you simply count all solution
attempts as failures, which would turn your's into an easy task.

Martin.

PS: As a comment on a recent fricas-devel post, Derive 6.10 solves
INT(1/PRODUCT(x - i*a, i, 0, 10), x) rapidly for both the factored and
the expanded integrand:

INT(1/(x*(x - a)*(x - 2*a)*(x - 3*a)*(x - 4*a)*(x - 5*a)*(x - 6*a)
*(x - 7*a)*(x - 8*a)*(x - 9*a)*(x - 10*a)), x)

INT(1/(3628800*a^10*(x - 10*a)) - 1/(362880*a^10*(x - 9*a))
+ 1/(80640*a^10*(x - 8*a)) - 1/(30240*a^10*(x - 7*a))
+ 1/(17280*a^10*(x - 6*a)) - 1/(14400*a^10*(x - 5*a))
+ 1/(17280*a^10*(x - 4*a)) - 1/(30240*a^10*(x - 3*a))
+ 1/(80640*a^10*(x - 2*a)) - 1/(362880*a^10*(x - a))
+ 1/(3628800*a^10*x), x)

Jeff Barnett

unread,
Sep 4, 2022, 1:24:56 PM9/4/22
to
What I'm trying to say and so far haven't been clear is that the
purposes of a test suite applied to multiple products include

1. Helping users and potential users decide which product is most
reliable, i.e., which ones give trustworthy answers and "I CAN'T DO IT"
when appropriate.

2. ------ decide which product is most robust, i.e., which one solves
the highest percentage of problems.

3. ------ decide which product is the most automatic, i.e., which one
needs the least user wizardry and divine guidance.

With the situation that you have uncovered, you need to think through
how to test such products. To be fair to others that don't make you
guess anything, perhaps the right thing to do is run things once with
system responsible for seeds and all such matters. Its fair to answer
questions such as "is x positive, negative, or zero?" but not "guess a
seed".

What I have been doing is help you to appreciate that the kind of
testing you do depends, radically, on what you want the test results to
influence and how. The hints of suggestions in prior messages were
basically to help you think this problem through more carefully. They
really were hints and not suggestions.

By the way I think its terrific that you are willing to spend your time
do this as "community service". Perhaps some of that community reading
this thread can offer suggestions about what results they find most
helpful and how the test construction would prove most helpful to them
and others. When you are testing systems whose behavior is
probabilistic, things get incredible difficult -- particularly in
domains where probabilities make no sense as they don't here.

I think I'll drop out of this thread now because I'm not a serious user
of the systems you are testing and I'm only guessing what is important.
Good luck with whatever you decide to do.
--
Jeff Barnett

Richard Fateman

unread,
Sep 13, 2022, 7:02:42 PM9/13/22
to
Can someone explain how this is even plausible as a method to use for symbolic integration?

Some probabilistic algorithms are used in computer algebra systems.. They work like this:
The system will pick a "random" number (perhaps a random prime number less than half the word size of
the computer you are using). If it picks a lucky number, the algorithm finishes fast.
If it pick an unlucky number, the algorithm notices that this is an unlucky number, and picks another one.
If it runs out of picks, the algorithm gives up with "failure of algorithm".

At no time is it acceptable for the algorithm to give a wrong answer, or for the answer to depend on the
particular first choice of a number.

It seems that this program is just wrong if it depends on the random number generator's seed being 12.

RJF
0 new messages