MaplePrimes Commons General Technical Discussions

The primary forum for technical discussions.

In the study of the Gödel spacetime model, a tetrad was suggested in the literature [1]. Alas, upon entering the tetrad in question, Maple's Tetrad's package complained that that matrix was not a tetrad! What went wrong? After an exchange with Edgardo S. Cheb-Terrab, Edgardo provided us with awfully useful comments regarding the use of the package and suggested that the problem together with its solution be presented in a post, as others may find it of some use for their work as well.

 

The Gödel spacetime solution to Einsten's equations is as follows.

 

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 858 and is the same as the version installed in this computer, created 2020, October 27, 10:19 hours Pacific Time.`

(1)

with(Physics); with(Tetrads)

_______________________________________________________

 

`Setting `*lowercaselatin_ah*` letters to represent `*tetrad*` indices`

 

((`Defined as tetrad tensors `*`see <a href='http://www.maplesoft.com/support/help/search.aspx?term=Physics,tetrads`*`,' target='_new'>?Physics,tetrads`*`,</a> `*`&efr;`[a, mu]*`, `)*eta[a, b]*`, `*gamma[a, b, c]*`, `)*lambda[a, b, c]

 

((`Defined as spacetime tensors representing the NP null vectors of the tetrad formalism `*`see <a href='http://www.maplesoft.com/support/help/search.aspx?term=Physics,tetrads`*`,' target='_new'>?Physics,tetrads`*`,</a> `*l[mu]*`, `)*n[mu]*`, `*m[mu]*`, `)*conjugate(m[mu])

 

_______________________________________________________

(2)

Working with Cartesian coordinates,

Coordinates(cartesian)

`Systems of spacetime coordinates are:`*{X = (x, y, z, t)}

 

{X}

(3)

the Gödel line element is

 

ds^2 = d_(t)^2-d_(x)^2-d_(y)^2+(1/2)*exp(2*q*y)*d_(z)^2+2*exp(q*y)*d_(z)*d_(t)

ds^2 = Physics:-d_(t)^2-Physics:-d_(x)^2-Physics:-d_(y)^2+(1/2)*exp(2*q*y)*Physics:-d_(z)^2+2*exp(q*y)*Physics:-d_(z)*Physics:-d_(t)

(4)

Setting the metric

Setup(metric = rhs(ds^2 = Physics[d_](t)^2-Physics[d_](x)^2-Physics[d_](y)^2+(1/2)*exp(2*q*y)*Physics[d_](z)^2+2*exp(q*y)*Physics[d_](z)*Physics[d_](t)))

_______________________________________________________

 

`Coordinates: `*[x, y, z, t]*`. Signature: `*`- - - +`

 

_______________________________________________________

 

Physics:-g_[mu, nu] = Matrix(%id = 18446744078354506566)

 

_______________________________________________________

 

`Setting `*lowercaselatin_is*` letters to represent `*space*` indices`

 

[metric = {(1, 1) = -1, (2, 2) = -1, (3, 3) = (1/2)*exp(2*q*y), (3, 4) = exp(q*y), (4, 4) = 1}, spaceindices = lowercaselatin_is]

(5)

The problem appeared upon entering the matrix M below supposedly representing the alleged tetrad.

interface(imaginaryunit = i)

M := Matrix([[1/sqrt(2), 0, 0, 1/sqrt(2)], [-1/sqrt(2), 0, 0, 1/sqrt(2)], [0, 1/sqrt(2), -I*exp(-q*y), I], [0, 1/sqrt(2), I*exp(-q*y), -I]])

Matrix(%id = 18446744078162949534)

(6)

Each of the rows of this matrix is supposed to be one of the null vectors [l, n, m, conjugate(m)]. Before setting this alleged tetrad, Maple was asked to settle the nature of it, and the answer was that M was not a tetrad! With the Physics Updates v.857, a more detailed message was issued:

IsTetrad(M)

`Warning, the given components form a`*null*`tetrad, `*`with a contravariant spacetime index`*`, only if you change the signature from `*`- - - +`*` to `*`+ - - -`*`. 
You can do that by entering (copy and paste): `*Setup(signature = "+ - - -")

 

false

(7)

So there were actually three problems:

1. 

The entered entity was a null tetrad, while the default of the Physics package is an orthonormal tetrad. This can be seen in the form of the tetrad metric, or using the library commands:

eta_[]

Physics:-Tetrads:-eta_[a, b] = Matrix(%id = 18446744078354552462)

(8)

Library:-IsOrthonormalTetradMetric()

true

(9)

Library:-IsNullTetradMetric()

false

(10)
2. 

The matrix M would only be a tetrad if the spacetime index is contravariant. On the other hand, the command IsTetrad will return true only when M represents a tetrad with both indices covariant. For  instance, if the command IsTetrad  is issued about the tetrad automatically computed by Maple, but is passed the matrix corresponding to "`&efr;`[a]^(mu)"  with the spacetime index contravariant,  false is returned:

"e_[a,~mu, matrix]"

Physics:-Tetrads:-e_[a, `~&mu;`] = Matrix(%id = 18446744078297840926)

(11)

"IsTetrad(rhs(?))"

Typesetting[delayDotProduct](`Warning, the given components form a`*orthonormal*`tetrad only if the spacetime index is contravariant. 
You can construct a tetrad with a covariant spacetime index by entering (copy and paste): `, Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = sqrt(2)*exp(-q*y), (3, 4) = -sqrt(2), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}), true).rhs(g[])

 

false

(12)
3. 

The matrix M corresponds to a tetrad with different signature, (+---), instead of Maple's default (---+). Although these two signatures represent the same physics, they differ in the ordering of rows and columns: the timelike component is respectively in positions 1 and 4.

 

The issue, then, became how to correct the matrix M to be a valid tetrad: either change the setup, or change the matrix M. Below the two courses of action are provided.

 

First the simplest: change the settings. According to the message (7), setting the tetrad to be null, changing the signature to be (+---) and indicating that M represents a tetrad with its spacetime index contravariant would suffice:

Setup(tetradmetric = null, signature = "+---")

[signature = `+ - - -`, tetradmetric = {(1, 2) = 1, (3, 4) = -1}]

(13)

The null tetrad metric is now as in the reference used.

eta_[]

Physics:-Tetrads:-eta_[a, b] = Matrix(%id = 18446744078298386174)

(14)

Checking now with the spacetime index contravariant

e_[a, `~&mu;`] = M

Physics:-Tetrads:-e_[a, `~&mu;`] = Matrix(%id = 18446744078162949534)

(15)

At this point, the command IsTetrad  provided with the equation (15), where the left-hand side has the information that the spacetime index is contravariant

"IsTetrad(?)"

`Type of tetrad: `*null

 

true

(16)

Great! one can now set the tetrad M exactly as entered, without changing anything else. In the next line it will only be necessary to indicate that the spacetime index, mu, is contravariant.

Setup(e_[a, `~&mu;`] = M, quiet)

[tetrad = {(1, 1) = -(1/2)*2^(1/2), (1, 3) = (1/2)*2^(1/2)*exp(q*y), (1, 4) = (1/2)*2^(1/2), (2, 1) = (1/2)*2^(1/2), (2, 3) = (1/2)*2^(1/2)*exp(q*y), (2, 4) = (1/2)*2^(1/2), (3, 2) = -(1/2)*2^(1/2), (3, 3) = ((1/2)*I)*exp(q*y), (3, 4) = 0, (4, 2) = -(1/2)*2^(1/2), (4, 3) = -((1/2)*I)*exp(q*y), (4, 4) = 0}]

(17)

 

The tetrad is now the matrix M. In addition to checking this tetrad making use of the IsTetrad command, it is also possible to check the definitions of tetrads and null vectors using TensorArray.

e_[definition]

Physics:-Tetrads:-e_[a, `&mu;`]*Physics:-Tetrads:-e_[b, `~&mu;`] = Physics:-Tetrads:-eta_[a, b]

(18)

TensorArray(Physics:-Tetrads:-e_[a, `&mu;`]*Physics:-Tetrads:-e_[b, `~&mu;`] = Physics:-Tetrads:-eta_[a, b], simplifier = simplify)

Matrix(%id = 18446744078353048270)

(19)

For the null vectors:

l_[definition]

Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-l_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[`~mu`] = 1, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-m_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-mb_[`~mu`] = 0, Physics:-g_[mu, nu] = Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[nu]+Physics:-Tetrads:-l_[nu]*Physics:-Tetrads:-n_[mu]-Physics:-Tetrads:-m_[mu]*Physics:-Tetrads:-mb_[nu]-Physics:-Tetrads:-m_[nu]*Physics:-Tetrads:-mb_[mu]

(20)

TensorArray([Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-l_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[`~mu`] = 1, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-m_[`~mu`] = 0, Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-mb_[`~mu`] = 0, Physics[g_][mu, nu] = Physics:-Tetrads:-l_[mu]*Physics:-Tetrads:-n_[nu]+Physics:-Tetrads:-l_[nu]*Physics:-Tetrads:-n_[mu]-Physics:-Tetrads:-m_[mu]*Physics:-Tetrads:-mb_[nu]-Physics:-Tetrads:-m_[nu]*Physics:-Tetrads:-mb_[mu]], simplifier = simplify)

[0 = 0, 1 = 1, 0 = 0, 0 = 0, Matrix(%id = 18446744078414241910)]

(21)

From its Weyl scalars, this tetrad is already in the canonical form for a spacetime of Petrov type "D": only `&Psi;__2` <> 0

PetrovType()

"D"

(22)

Weyl[scalars]

psi__0 = 0, psi__1 = 0, psi__2 = -(1/6)*q^2, psi__3 = 0, psi__4 = 0

(23)

Attempting to transform it into canonicalform returns the tetrad (17) itself

TransformTetrad(canonicalform)

Matrix(%id = 18446744078396685478)

(24)

Let's now obtain the correct tetrad without changing the signature as done in (13).

Start by changing the signature back to "(- - - +)"

Setup(signature = "---+")

[signature = `- - - +`]

(25)

So again, M is not a tetrad, even if the spacetime index is specified as contravariant.

IsTetrad(e_[a, `~&mu;`] = M)

`Warning, the given components form a`*null*`tetrad, `*`with a contravariant spacetime index`*`, only if you change the signature from `*`- - - +`*` to `*`+ - - -`*`. 
You can do that by entering (copy and paste): `*Setup(signature = "+ - - -")

 

false

(26)

By construction, the tetrad M has its rows formed by the null vectors with the ordering [l, n, m, conjugate(m)]. To understand what needs to be changed in M, define those vectors, independent of the null vectors [l_, n_, m_, mb_] (with underscore) that come with the Tetrads package.

Define(l[mu], n[mu], m[mu], mb[mu], quiet)

and set their components using the matrix M taking into account that its spacetime index is contravariant, and equating the rows of M  using the ordering [l, n, m, conjugate(m)]:

`~`[`=`]([l[`~&mu;`], n[`~&mu;`], m[`~&mu;`], mb[`~&mu;`]], [seq(M[j, 1 .. 4], j = 1 .. 4)])

[l[`~&mu;`] = Vector[row](%id = 18446744078368885086), n[`~&mu;`] = Vector[row](%id = 18446744078368885206), m[`~&mu;`] = Vector[row](%id = 18446744078368885326), mb[`~&mu;`] = Vector[row](%id = 18446744078368885446)]

(27)

"Define(op(?))"

`Defined objects with tensor properties`

 

{Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-Tetrads:-e_[a, mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-gamma_[i, j], Physics:-Tetrads:-gamma_[a, b, c], l[mu], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], m[mu], Physics:-Tetrads:-m_[mu], mb[mu], Physics:-Tetrads:-mb_[mu], n[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(28)

Check the covariant components of these vectors towards comparing them with the lines of the Maple's tetrad `&efr;`[a, mu]

l[], n[], m[], mb[]

l[mu] = Array(%id = 18446744078298368710), n[mu] = Array(%id = 18446744078298365214), m[mu] = Array(%id = 18446744078298359558), mb[mu] = Array(%id = 18446744078298341734)

(29)

This shows the [l_, n_, m_, mb_] null vectors (with underscore) that come with Tetrads package

e_[nullvectors]

Physics:-Tetrads:-l_[mu] = Vector[row](%id = 18446744078354520414), Physics:-Tetrads:-n_[mu] = Vector[row](%id = 18446744078354520534), Physics:-Tetrads:-m_[mu] = Vector[row](%id = 18446744078354520654), Physics:-Tetrads:-mb_[mu] = Vector[row](%id = 18446744078354520774)

(30)

So (29) computed from M is the same as (30) computed from Maple's tetrad.

But, from (30) and the form of Maple's tetrad

e_[]

Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078297844182)

(31)

for the current signature

Setup(signature)

[signature = `- - - +`]

(32)

we see the ordering of the null vectors is [n, m, mb, l], not [l, n, m, mb] used in [1] with the signature (+ - - -). So the adjustment required in  M, resulting in "M^( ')", consists of reordering M's rows to be [n, m, mb, l]

`#msup(mi("M"),mrow(mo("&InvisibleTimes;"),mo("&apos;")))` := simplify(Matrix(4, map(Library:-TensorComponents, [n[mu], m[mu], mb[mu], l[mu]])))

Matrix(%id = 18446744078414243230)

(33)

IsTetrad(`#msup(mi("M"),mrow(mo("&InvisibleTimes;"),mo("&apos;")))`)

`Type of tetrad: `*null

 

true

(34)

Comparing "M^( ')" with the tetrad `&efr;`[a, mu]computed by Maple ((24) and (31), they are actually the same.

References

[1]. Rainer Burghardt, "Constructing the Godel Universe", the arxiv gr-qc/0106070 2001.

[2]. Frank Grave and Michael Buser, "Visiting the Gödel Universe",  IEEE Trans Vis Comput GRAPH, 14(6):1563-70, 2008.


 

Download Godel_universe_and_Tedrads.mw

Hi,
Some people using the Windows platform have had problems installing MapleCloud packages, including the Maplesoft Physics Updates. This problem does not happen in Macintosh or Linux/Unix, also does not happen with all Windows computers but with some of them, and is not a problem of the MapleCloud packages themselves, but a problem of the installer of packages.

I understand that a solution to this problem will be presented within an upcoming Maple dot release.

Meantime, there is a solution by installing a helper library; after that, MapleCloud packages install without problems in all Windows machines. So whoever is having trouble installing MapleCloud packages in Windows and prefers not to wait until that dot release, instead wants to install this helper library, please email me at physics@maplesoft.com

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hello each and everyone,

I am still not able to install/update the Physics package and this since version 713.

I think the phenomenom happens when Maple 2020 was updated to version 2020.1

The 2 ways fail :

1) Physics:-Version(latest)  ==> kernel crashes (lost connection)

2) via the Cloud : it's even worse... "installation" lasts forever (i usually break it after 15 (!) minutes)

I removed totally Maple 2020.1 and reinstalled it, the problem is still there.

I wonder when Maplesoft developpers will fix this problem whitch has been around for a long time now.

I do *not* want to do it by hand as one of this forum contributor suggested me a while ago.

I use the Physics paxkage intensively so i'm hampered in my work.

Kind regards to all.

Jean-Michel

There is an interesting preprint here, on a method for integration of some pseudo-elliptic integrals.

It was mentioned in a (sci.math.symbolic) usenet thread, which can be accessed via Google Groups here.

Inside the paper, the author mentions that the following is possible for some cases -- to first convert to RootOf form.

restart;

kernelopts(version);

`Maple 2020.0, X86 64 LINUX, Mar 4 2020, Build ID 1455132`

ig:=((-3+x^2)*(1-6*x^2+x^4)^(-1/4))/(-1+x^2);

(x^2-3)/((x^4-6*x^2+1)^(1/4)*(x^2-1))

int(ig,x);

int((x^2-3)/((x^4-6*x^2+1)^(1/4)*(x^2-1)), x)

infolevel[int] := 2:

S:=simplify([allvalues(int(convert(ig,RootOf),x))],size):

Stage1: first-stage indefinite integration

Stage2: second-stage indefinite integration

Norman: enter Risch-Norman integrator

Norman: exit Risch-Norman integrator

int/algrisch/int: Risch/Trager's algorithm for algebraic function

int/algrisch/int: entering integrator at time 9.029

int/algrisch/int: function field has degree 4

int/algrisch/int: computation of an integral basis: start time 9.031

int/algrisch/int: computation of an integral basis: end time 9.038

int/algrisch/int: normalization at infinity: start time 9.039

int/algrisch/int: normalization at infinity: end time 9.048

int/algrisch/int: genus of the function field 3

int/algrisch/int: computation of the algebraic part: start time 9.059

int/algrisch/int: computation of the algebraic part: end time 9.060

int/algrisch/int: computation of the transcendental part: start time 9.063

int/algrisch/transcpar: computing a basis for the residues at time 9.068

int/algrisch/residues: computing a splitting field at time 9.068

int/algrisch/transcpar: basis for the residues computed at time 9.103

int/algrisch/transcpar: dimension is 2

int/algrisch/transcpar: building divisors at time 9.300

int/algrisch/transcpar: testing divisors for principality at time 9.605

int/algrisch/goodprime: searching for a good prime at time 9.606

int/algrisch/goodprime: good prime found at time 9.704

int/algrisch/goodprime: searching for a good prime at time 9.704

int/algrisch/goodprime: good prime found at time 9.762

int/algrisch/areprinc: the divisor is principal: time 10.084

int/algrisch/areprinc: the divisor is principal: time 11.833

int/algrisch/transcpar: divisors proven pincipal at time at time 11.833

int/algrisch/transcpar: generators computed at time 11.834

int/algrisch/transcpar: orders are [1 1]

int/algrisch/transcpar: check that the candidate is an actual antiderivative

int/algrisch/transcpar: the antiderivative is elementary

int/algrisch/transcpar: antiderivative is (1/2)*ln((RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^3*_z^3-RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^4+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^2-3*_z^3*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)*_z-5*_z^2+1)/((_z+1)*(_z-1)*_z^2))+(1/2)*RootOf(_Z^2+1)*ln((RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*RootOf(_Z^2+1)*_z^4+RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^3*_z^3-RootOf(_Z^2+1)*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)^2*_z^2+3*_z^3*RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)-5*RootOf(_Z^2+1)*_z^2-RootOf(_Z^4*_z^4-_z^4+6*_z^2-1 index = 1)*_z+RootOf(_Z^2+1))/((_z+1)*(_z-1)*_z^2))

int/algrisch/int: computation of the transcendental part: end time 12.040

int/algrisch/int: exiting integrator for algebraic functions at time 12.041

S[1];

(1/2)*ln(-((x^4-6*x^2+1)^(3/4)*x+(x^4-6*x^2+1)^(1/2)*x^2+(x^4-6*x^2+1)^(1/4)*x^3+x^4-(x^4-6*x^2+1)^(1/2)-3*(x^4-6*x^2+1)^(1/4)*x-5*x^2)/((x+1)*(x-1)))+((1/2)*I)*ln((-(x^4-6*x^2+1)^(3/4)*x+(I*x^2-I)*(x^4-6*x^2+1)^(1/2)+(x^3-3*x)*(x^4-6*x^2+1)^(1/4)-I*x^2*(x^2-5))/((x+1)*(x-1)))

S[2];

(1/2)*ln(-((x^4-6*x^2+1)^(3/4)*x+(x^4-6*x^2+1)^(1/2)*x^2+(x^4-6*x^2+1)^(1/4)*x^3+x^4-(x^4-6*x^2+1)^(1/2)-3*(x^4-6*x^2+1)^(1/4)*x-5*x^2)/((x+1)*(x-1)))-((1/2)*I)*ln((-(x^4-6*x^2+1)^(3/4)*x+(-I*x^2+I)*(x^4-6*x^2+1)^(1/2)+(x^3-3*x)*(x^4-6*x^2+1)^(1/4)+I*x^2*(x^2-5))/((x+1)*(x-1)))

simplify( diff(S[1],x) - ig ), simplify( diff(S[2],x) - ig );

0, 0

 

Download int_pe.mw

I’m extremely pleased to introduce the newest update to the Maple Companion. In this time of wide-spread remote learning, tools like the Maple Companion are more important than ever, and I’m happy that our efforts are helping students (and some of their parents!) with at least one small aspect of their life.  Since we’ve added a lot of useful features since I last posted about this free mobile app, I wanted to share the ones I’m most excited about. 

(If you haven’t heard about the Maple Companion app, you can read more about it here.) 

If you use the app primarily to move math into Maple, you’ll be happy to hear that the automatic camera focus has gotten much better over the last couple of updates, and with this latest update, you can now turn on the flash if you need it. For me, these changes have virtually eliminated the need to fiddle with the camera to bring the math in focus, which sometimes happened in earlier versions.

If you use the app to get answers on your phone, that’s gotten much better, too. You can now see plots instantly as you enter your expression in the editor, and watch how the plot changes as you change the expression. You can also get results to many numerical problems results immediately, without having to switch to the results screen. This “calculator mode” is available even if you aren’t connected to the internet.  Okay, so there aren’t a lot of students doing their homework on the bus right now, but someday!

Speaking of plots, you can also now view plots full-screen, so you can see more of plot at once without zooming and panning, squinting, or buying a bigger phone.

Finally, if English is not you or your students’ first language, note that the app was recently made available in Spanish, French, German, Russian, Danish, Japanese, and Simplified Chinese. 

As always, we’d love you hear your feedback and suggestions. Please leave a comment, or use the feedback forms in the app or our web site.

Visit Maple Companion to learn more, find links to the app stores so you can download the app, and access the feedback form. If you already have it installed, you can get the new release simply by updating the app on your phone.

The ideas here are to allow 3D plotting commands such as plot3d to handle a `size` option similarly to how 2D plotting commands do so, and for the plots:-display command to also handle it for 3D plots.

The size denotes the dimensions of the inlined plotting window, and not the relative lengths of the three axes.

I'd be interested in any new problems introduced with this, eg. export, etc.

restart;

#
# Using ToInert/FromInert
#
# This might go in an initialzation file.
#
try
  __ver:=(u->:-parse(u[7..:-StringTools:-Search(",",u)-1]))(:-sprintf("%s",:-kernelopts(':-version')));
  if __ver>=18.0 and __ver<=2019.2 then
    __KO:=:-kernelopts(':-opaquemodules'=false);
    :-unprotect(:-Plot:-Options:-Verify:-ProcessParameters);
    __KK:=ToInert(eval(:-Plot:-Options:-Verify:-ProcessParameters)):
    __LL:=[:-op([1,2,1,..],__KK)]:
    __NN:=:-nops(:-remove(:-type,[:-op([1,..],__KK)],
                          ':-specfunc(:-_Inert_SET)'))
          +:-select(:-has,[:-seq([__i,__LL[__i]],
                                 __i=1..:-nops(__LL))],
                    "size")[1][1];
    if :-has(:-op([5,2,2,2,1],__KK),:-_Inert_PARAM(__NN)) then
      __KK:=:-subsop([5,2,2,2,1]
                     =:-subs([:-_Inert_PARAM(__NN)=:-NULL],
                              :-op([5,2,2,2,1],__KK)),__KK);
      :-Plot:-Options:-Verify:-ProcessParameters:=
      :-FromInert(:-subsop([5,2,2,3,1]
                  =:-subs([:-_Inert_STRING("size")=:-NULL],
                          :-op([5,2,2,3,1],__KK)),__KK));
      :-print("3D size patch done");
    else
      :-print("3D size patch not appropriate; possibly already done");
    end if;
  else
    :-print(sprintf("3D size patch not appropriate for version %a"),__ver);
  end if;
catch:
  :-print("3D size patch failed");
finally
  :-protect(:-Plot:-Options:-Verify:-ProcessParameters);
  :-kernelopts(':-opaquemodules'=__KO);
end try:

"3D size patch done"

 

P := plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, size=[150,150],
            font=[Times,5], labels=["","",""]):
P;

plots:-display(P, size=[300,300], font=[Times,10]);

#
# inherited from the contourplot3d (the plot3d is unset).
#
plots:-display(
  plots:-contourplot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1,
                       thickness=3, contours=20, size=[800,800]),
  plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, color="Gray",
         transparency=0.1, style=surface)
);

# Some options should still act as 2D-plot-specific.
#
try plot3d(sin(x)*y^2, x=-Pi..Pi, y=-1..1, legend="Q");
    print("Not OK");
catch:
if StringTools:-FormatMessage(lastexception[2..-1])
   ="the legend option is not available for 3-D plots"
then print("OK"); else print("Not OK"); error; end if; end try;

"OK"

 

Download 3Dsize_hotedit.mw

If this works fine then it might be a candidate for inclusion in an initialization file, so that it's
automatically available.

Hi, 

This is more of an open discussion than a real question. Maybe it would gain to be displaced in the post section?

Working with discrete random variables I found several inconsistencies or errors.
In no particular order: 

  • The support of a discrete RV is not defined correctly (a real range instead of a countable set)
  • The plot of the probability function (which, in my opinion, would gain to be renamed "Probability Mass Function, see https://en.wikipedia.org/wiki/Probability_mass_function) is not correct.
  • The  ProbabiliytFunction of a discrte rv of EmpiricalDistribution can be computed at any point, but its formal expression doesn't exist (or at least is not accessible).
  • Defining the discrete rv "toss of a fair dice"  with EmpiricalDistribution and DiscreteUniform gives different results.


The details are given in the attached file and I do hope that the companion text is clear enough to point the issues.
I believe there is no major issues here, but that Maple suffers of some lack of consistencies in the treatment of discrete (at least some) rvs. Nothing that could easily be fixed.


As I said above, if some think this question has no place here and ought to me moved to the post section, please feel free to do it.

Thanks for your attention.


 

restart:

with(Statistics):


Two alternate ways to define a discrete random variable on a finite set
of equally likely outcomes.

Universe    := [$1..6]:
toss_1_dice := RandomVariable(EmpiricalDistribution(Universe));
TOSS_1_DICE := RandomVariable(DiscreteUniform(1, 6));

_R

 

_R0

(1)


Let's look to the ProbabilityFunction of each RV

ProbabilityFunction(toss_1_dice, x);
ProbabilityFunction(TOSS_1_DICE, x);

"_ProbabilityFunction[Typesetting:-mi("x",italic = "true",mathvariant = "italic")]"

 

piecewise(x < 1, 0, x <= 6, 1/6, 6 < x, 0)

(2)


It looks like the procedure ProbabilityFunction is not an attribute of RV with EmpiticalDistribution.
Let's verify

law := [attributes(toss_1_dice)][3]:
lprint(exports(law))

Conditions, ParentName, Parameters, CDF, DiscreteValueMap, Mean, Median, Mode, ProbabilityFunction, Quantile, Specialize, Support, RandomSample, RandomVariate

 


Clearly ProbabilityFunction is an attribute of toss_1_dice.

In fact it appears the explanation of the difference of behaviours relies upon different definitions
of the set of outcomes of toss_1_dice and TOSS_1_DICE

LAW := [attributes(TOSS_1_DICE)][3]:
exports(LAW):

law:-Conditions;
LAW:-Conditions;

[(Vector(6, {(1) = 1, (2) = 2, (3) = 3, (4) = 4, (5) = 5, (6) = 6}))::rtable]

 

[1 < 6]

(3)


From :-Conditions one can see that toss_1_dice is realy a discrete RV defined on a countable set of outcomes,
but that nothing is said about the set over which TOSS_1_DICE is defined.

The truly discrete definition of toss_1_dice is confirmed here :
(the second result is correct

ProbabilityFinction(toss_1_dice, x) = {0 if x < 1, 0 if x > 6, 1/6 if x::integer, 0 otherwise

ProbabilityFunction~(toss_1_dice, Universe);
ProbabilityFunction~(toss_1_dice, [seq(0..7, 1/2)]);

[1/6, 1/6, 1/6, 1/6, 1/6, 1/6]

 

[0, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0, 1/6, 0, 0]

(4)


One can also see that the Support of both of these RVs are wrong

(see for instance https://en.wikipedia.org/wiki/Discrete_uniform_distribution)

There should be {1, 2, 3, 4, 5, 6}, not a RealRange.

Support(toss_1_dice);
Support(TOSS_1_DICE);

RealRange(1, 6)

 

RealRange(1, 6)

(5)

 

0

 

{1, 2, 3, 4, 5, 6}

 

 


Now this is the surprising ProbabilityFunction of TOSS_1_DICE.
This obviously wrong result probably linked to the weak definition of the conditions for this RB.

# plot(ProbabilityFunction(TOSS_1_DICE, x), x=0..7);
plot(ProbabilityFunction(TOSS_1_DICE, x), x=0..7, discont=true)

 


These differences of treatments raise a lot of questions :
    -  Why is a DiscreteUniform RV not defined on a countable set?
    -  Why does the ProbabilityFunction of an EmpiricalDistribution return no result
        if its second parameter is not set to one  its outcomes.

 All this without even mentioning the wrong plot shown above.
 

I believe something which would work like the module below would be much better than what is done

right now

 

EmpiricalRV := module()
export MassDensityFunction, PlotMassDensityFunction, Support:

MassDensityFunction := proc(rv, x)
  local u, v, N:
  u := [attributes(rv)][3]:
  if u:-ParentName = EmpiricalDistribution then
    v := op([1, 1], u:-Conditions);
    N := numelems(v):
    return piecewise(op(op~([seq([x=v[n], 1/N], n=1..N)])), 0)
  else
    error "The random variable does not have an EmpiricalDistribution"
  end if
end proc:

PlotMassDensityFunction := proc(rv, x1, x2)
  local u, v, a, b:
  u := [attributes(rv)][3]:
  if u:-ParentName = EmpiricalDistribution then
    v := op([1, 1], u:-Conditions);
    a := select[flatten](`>=`, v, x1);
    b := select[flatten](`<=`, a, x2);
    PLOT(seq(CURVES([[n, 0], [n, 1/numelems(v)]], COLOR(RGB, 0, 0, 1), THICKNESS(3)), n in b), VIEW(x1..x2, default))
  else
    error "The random variable does not have an EmpiricalDistribution"
  end if
end proc:

Support := proc(rv, x1, x2)
  local u, v, a, b:
  u := [attributes(rv)][3]:
  if u:-ParentName = EmpiricalDistribution then
    v := op([1, 1], u:-Conditions);
    return {entries(v, nolist)}
  else
    error "The random variable does not have an EmpiricalDistribution"
  end if
end proc:

end module:
 

EmpiricalRV:-MassDensityFunction(toss_1_dice, x);
 

piecewise(x = 1, 1/6, x = 2, 1/6, x = 3, 1/6, x = 4, 1/6, x = 5, 1/6, x = 6, 1/6, 0)

(6)

f := unapply(EmpiricalRV:-MassDensityFunction(toss_1_dice, x), x):
f(2);
f(5/2);
 

1/6

 

0

(7)

EmpiricalRV:-PlotMassDensityFunction(toss_1_dice, 0, 7);

 

 


 

Download Discrete_RV.mw

 

 

I'm particularly interested in data analysis and more specifically in statistical analysis of computer code outputs.

One of the main activity of this very broad field is named Uncertainty Propagation. In a few words it consists in perturbing the inputs of a computational code in order to understand (and quantify) how these perturbations propagates through the outputs of this code.

At the core of uncertainty propagation is the ability to generate large numbers of "random" variations of the inputs. Knowing that these entries can be counted in tens, one sees that the first problem consists in generating "random" points in a space of potentially very large dimension.

Even among my mathematician colleagues an impressive number of them is completely ignorant of the way "random" numbers are generated. I guess that a lot of Mapleprimes' users are too. My purpose is not to give a course on this topic and the affording litterature is vast enough for everyone interested might find informations of any level of complexity.
Among those who have some knowledge about Pseudo Random Numbers Generators (PRNG), only a few of them know that a PRNG has to pass severe tests ("tests of randomness") before the streams of number it generates might  be qualified as "reasonably random" and therefore this PRNG might be released.

One of most famous example of a bad PRNG is given by "randu" (IBM 1966, and probably used in Fortran libraries during more than 30 years), this same PRNG that Knuth qualified himself as the "infamous generator".

These tests of randomness are generally gathered in dedicated libraries and Diehard is probably tone of the most known of them.
Diehard has originally been developed by George Marsaglia more than twenty years ago and it's still widely ued today.

I recently decided, not because I have doubts about the quality of the work done by Maplesoft, to test the Maple's PRNG named "Mersenne Twister". First, because it can do no harm to publish quantitative information that allows everyone to know that it is using a proven PRNG; second, because the (very simple) approach used here can fill the gaps I have mentioned above.

Mersenne Twister (often dubbed mt19937) is considered as a very good PRNG; it is used in a lot of applications (including finance where it is not so rare to sample input spaces of dimensions larger than 1000... ok I know, mt19937 is often considered as a poor candidate for cryptography applications, but it's not my concern here).

I have thus decided to spend some time to run the Diehard suite of tests on a sequence of integers numbers generated by RandomTools[MersenneTwister].


 

restart:


DIEHARD tests suite for Pseudo Random Numbers Generators (PRNG)

Reference: http://webhome.phy.duke.edu/~rgb/General/dieharder.php

The installation procedure (Mac OSX) can be found here
    https://gist.github.com/blixt/9abfafdd0ada0f4f6f26
or here
    http://macappstore.org/dieharder/

For other operating systems, please search on the web pages.


dieharder [-h]   # for inline help
dieharder -l      # to get the lists all the avaliable tests




A description of the many tests can be found here:
    https://en.wikipedia.org/wiki/Diehard_tests
    https://sites.google.com/site/astudyofentropy/background-information/the-tests/dieharder-test-descriptions
    https://www.stata.com/support/cert/diehard/randnumb_mt64.out

General theory about PRNG testing can be found here (a reference among many):
    http://liu.diva-portal.org/smash/get/diva2:740158/FULLTEXT01.pdf

or here (more oriented to the NIST test suite)
    https://www.random.org/analysis/Analysis2005.pdf
    https://nvlpubs.nist.gov/nistpubs/legacy/sp/nistspecialpublication800-22r1a.pdf



In a terminal window execute the following commands for an exhaustive testing ("-a" option).
The "-g 202" option means that the generator is replaced by a text format input file
(use dieharder -h for more details).

cd //..../Desktop/DIEHARD

dieharder -g 202 -f SomeAsciiFile -a > //..../Desktop/DIEHARD/TheResultFile.txt

Be carefull, the complete testing takes several hours (about 5 on my computer)



__________________________________________________________________________________
 


Maple's Mersenne Twister Generator

Maple help page : RandomTools[MersenneTwister][GenerateInteger]
(see rincluded references to the Mersenne Twister PRNG).

Note: in the sequel this generator will be dubbed mt19937


The Mersenne Twister is implemented in many softwares.
It is higly likely that this PRNG (and the others these softwares propose) have been intensively
tested with one of the existing PRNG testing libraries.
Unfortunately only a few editors have made public the results of these tests (probably because
the implementation in itself is rarely questioned... but a code typo is always a possibility).

One exception is ths software STATA.
A summary of the results can be found here
   https://www.stata.com/support/cert/diehard/.
A complete description of the results of the tests passed is given here
   https://www.stata.com/support/cert/diehard/randnumb_mt64.out

The classical pattern of the performances of mt19937 can be found here

   http://www2.ic.uff.br/~celso/artigos/pjo6.ps.

and the table below comes from it (P means "Passed", F means "Failed"):


____________________________________________________________________________


In the Maple code below, a sequence of N UnsignedInt32 numbers is generated from the
Maple's Mersenne Twister and the result is exported in an ASCII file.
The Seed is set to 1 (SetState(state=1)) to compare, with a small value of N (let's say N=10)
the sequence produced by Maple's mt19937 with the the sequence of the same length generated
by Diehard's mt19937.
To generate this later sequence and save it in file Diehard_mt19937, just run in a terminan window
the command (-S 1 means "seed = 1", -t 10 means "a sequence of length 10"):
   dieharder -S 1 -B -o -t 10 > Diehard_mt19937

About the value of N:

In http://webhome.phy.duke.edu/~rgb/General/dieharder.php it's recommend that N be at least
equal to 2.5 million; STATA used N=3 million.
Other web sources say this value is too small.
For N=10 million the Maple's mt19937 doesn't pass the tests successfully.
I used here N=50 million (the resulting ASCII file has size 537 Mo).



Name of the input file.

The file generated by Maple is named Maple_mt19937_N=5e7.txt



One important thing is the preamble of a licit input file.

This preamble must have 6 lines (the value 10 right to count must be set to the value of N).
A licit preamble is of the form.

#==================================================================

# some text indicating the generator used

#==================================================================

type: d

count: 10

numbit: 32

As Maple_mt19937_N=5e7.txt is generated from an ExportMatrix command, this preamble is added
by hand.
 


Running multiple Diehard tests

To run the same tests used to qualify STATA's Mersenne Twister, open a terminal window,
go to the directory that contains input file Maple_mt19937_N=5e7.txt and run this script:

 for i in {0,1,2,3,4,8,9,10,11,12,13,14,15,16}; do

    dieharder -g 202 -f Maple_mt19937_N=5e7.txt -d $i >> Diehard___Maple_mt19937_N=5e7

 done ;

The results are then forked in the ASCII file Diehard___Maple_mt19937_N=5e7

 

with(RandomTools[MersenneTwister]):

dir := cat("/", currentdir(), "Desktop/DIEHARD/"):
InputFile := cat(dir, "Maple_mt19937_N=5e7.txt"):

SetState(state=1);

N := 5*10^7:

st := time():
S := convert([seq(GenerateUnsignedInt32(), i=1..N)], Matrix)^+;
time()-st;

S := Vector(4, {(1) = ` 50000000 x 1 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

84.526

(1)

st := time():
ExportMatrix(InputFile, S, format=rectangular, mode=ascii);
time()-st;

537066525

 

61.926

(2)


Diehard's results


Full test suite (about 5 hours of computational time)

Command :
dieharder -g 202 -f Maple_mt19937_N=5e7.txt -a > Diehard___ALL___Maple_mt19937_N=5e7


The results are compared to those obtained for Diehard's mt19937.
Two ways are used :

  - 1 - In a first stage one generates a stream of PRN and store it in an ASCII file (just as we did with Maple).
         The whole suite of tests is then run on this file.
         Commands (-g 013 codes for mt19937):

         dieharder -S 1 -g 013 -o -t 50000000 > Diehard_mt19937_N=5e7.txt
         dieharder -g 202 -f Diehard_mt19937_N=5e7.txt -a > Diehard___ALL___Diehard_mt19937_N=5e7



  - 2 - The whole suite is run by invoking directectly mt19937 "online"
         Commands :
         dieharder -S 1 -g 013 -t 50000000 -a > Diehard___ALL___Online


A UNIX diff command has been used to verify that the two files Maple_mt19937_N=5e7.txt and
 Diehard_mt19937_N=5e7.txt were identical (thet were).

Note that the Diehard doens't responds identically depending on the stream of random numbers comes from a file
or is generated online (this last [- 2 -] situation seems to give better results).-

Résumé (114 tests):
   - * - Maple's  and Diehard's  mt19937 respond exactly the same way when the stream of random
          numbers is read from an ASCII file (8 tests failed (******) and 6 weak (**)).
   - * - Diehard's  mt19937 fails 0 test and is weak on 4 tests when the stream is generated online
 

 

restart:

dir := currentdir():
FromMapleFile     := cat(dir, "Diehard___ALL___Maple_mt19937_N=5e7"):
FromDiehardFile   := cat(dir, "Diehard___ALL___diehard_mt19937_N=5e7"):
FromDiehardNoFile := cat(dir, "Diehard___ALL___Online"):


printf("                           ======================|======================|======================|\n"):
printf("                          |   From Maple's file  | From Diehard's File  | Diehard online test  |\n"):
printf("==========================|======================|======================|======================|\n"):
printf("          test       ntup | p.value   Assessment | p.value   Assessment | p.value   Assessment |\n"):
printf("==========================|======================|======================|======================|\n"):


for k from 1 to 9 do
  LMF  := readline(FromMapleFile):
  LDF  := readline(FromDiehardFile):
  LDNF := readline(FromDiehardNoFile):
end do:


while LMF <> 0 do
  if StringTools:-Search("|", LMF) > 0 then
    res := StringTools:-StringSplit(LMF, "|")[[1, 2, 5, 6]];
    printf("%-20s  %3d | %1.7f ", res[1], parse(res[2]), parse(res[3]));
      if StringTools:-Search("WEAK"  , res[4]) > 0 then printf("    **     |")
    elif StringTools:-Search("FAILED", res[4]) > 0 then printf("  ******   |")
    else printf("  PASSED   |")
    end if:
  end if:
  LMF  := readline(FromMapleFile):

  if StringTools:-Search("|", LDF) > 0 then
    res := StringTools:-StringSplit(LDF, "|")[[5, 6]];
    printf(" %1.7f ", parse(res[1]));
      if StringTools:-Search("  WEAK"  , res[2]) > 0 then printf("     **    |")
    elif StringTools:-Search("  FAILED", res[2]) > 0 then printf("   ******  |")
    else printf("   PASSED  |")
    end if:
  end if:
  LDF  := readline(FromDiehardFile):
                     
  if StringTools:-Search("|", LDNF) > 0 then
    res := StringTools:-StringSplit(LDNF, "|")[[5, 6]];
    printf(" %1.7f ", parse(res[1]));
      if StringTools:-Search("WEAK"  , res[2]) > 0 then printf("     **    |")
    elif StringTools:-Search("FAILED", res[2]) > 0 then printf("   ******    |")
    else printf("   PASSED  |")
    end if:
    printf("\n"):
  end if:
  LDNF := readline(FromDiehardNoFile):


end do:

                           ======================|======================|======================|
                          |   From Maple's file  | From Diehard's File  | Diehard online test  |
==========================|======================|======================|======================|
          test       ntup | p.value   Assessment | p.value   Assessment | p.value   Assessment |
==========================|======================|======================|======================|
   diehard_birthdays    0 | 0.9912651   PASSED   | 0.9912651    PASSED  | 0.8284550    PASSED  |
      diehard_operm5    0 | 0.1802226   PASSED   | 0.1802226    PASSED  | 0.5550587    PASSED  |
  diehard_rank_32x32    0 | 0.3099035   PASSED   | 0.3099035    PASSED  | 0.9575440    PASSED  |
    diehard_rank_6x8    0 | 0.2577249   PASSED   | 0.2577249    PASSED  | 0.3915666    PASSED  |
   diehard_bitstream    0 | 0.5519218   PASSED   | 0.5519218    PASSED  | 0.9999462      **    |
        diehard_opso    0 | 0.1456442   PASSED   | 0.1456442    PASSED  | 0.7906533    PASSED  |
        diehard_oqso    0 | 0.4882425   PASSED   | 0.4882425    PASSED  | 0.9574014    PASSED  |
         diehard_dna    0 | 0.0102880   PASSED   | 0.0102880    PASSED  | 0.5149193    PASSED  |
diehard_count_1s_str    0 | 0.1471956   PASSED   | 0.1471956    PASSED  | 0.9517290    PASSED  |
diehard_count_1s_byt    0 | 0.1158707   PASSED   | 0.1158707    PASSED  | 0.1568255    PASSED  |
 diehard_parking_lot    0 | 0.1148982   PASSED   | 0.1148982    PASSED  | 0.1611173    PASSED  |
    diehard_2dsphere    2 | 0.9122204   PASSED   | 0.9122204    PASSED  | 0.2056657    PASSED  |
    diehard_3dsphere    3 | 0.9385972   PASSED   | 0.9385972    PASSED  | 0.3620517    PASSED  |
     diehard_squeeze    0 | 0.2686977   PASSED   | 0.2686977    PASSED  | 0.8611266    PASSED  |
        diehard_sums    0 | 0.1602355   PASSED   | 0.1602355    PASSED  | 0.5103248    PASSED  |
        diehard_runs    0 | 0.1235328   PASSED   | 0.1235328    PASSED  | 0.9402086    PASSED  |
        diehard_runs    0 | 0.6341956   PASSED   | 0.6341956    PASSED  | 0.3274267    PASSED  |
       diehard_craps    0 | 0.0243605   PASSED   | 0.0243605    PASSED  | 0.1844482    PASSED  |
       diehard_craps    0 | 0.2952043   PASSED   | 0.2952043    PASSED  | 0.1407422    PASSED  |
 marsaglia_tsang_gcd    0 | 0.0000000   ******   | 0.0000000    ******  | 0.5840531    PASSED  |
 marsaglia_tsang_gcd    0 | 0.0000000   ******   | 0.0000000    ******  | 0.8055035    PASSED  |
         sts_monobit    1 | 0.9397218   PASSED   | 0.9397218    PASSED  | 0.9018886    PASSED  |
            sts_runs    2 | 0.8092469   PASSED   | 0.8092469    PASSED  | 0.2247600    PASSED  |
          sts_serial    1 | 0.2902851   PASSED   | 0.2902851    PASSED  | 0.9223063    PASSED  |
          sts_serial    2 | 0.9541680   PASSED   | 0.9541680    PASSED  | 0.6140772    PASSED  |
          sts_serial    3 | 0.4090798   PASSED   | 0.4090798    PASSED  | 0.2334754    PASSED  |
          sts_serial    3 | 0.5474851   PASSED   | 0.5474851    PASSED  | 0.7370361    PASSED  |
          sts_serial    4 | 0.7282286   PASSED   | 0.7282286    PASSED  | 0.2518826    PASSED  |
          sts_serial    4 | 0.9905724   PASSED   | 0.9905724    PASSED  | 0.6876253    PASSED  |
          sts_serial    5 | 0.8297711   PASSED   | 0.8297711    PASSED  | 0.2123014    PASSED  |
          sts_serial    5 | 0.9092172   PASSED   | 0.9092172    PASSED  | 0.3532615    PASSED  |
          sts_serial    6 | 0.4976615   PASSED   | 0.4976615    PASSED  | 0.9967160      **    |
          sts_serial    6 | 0.9853355   PASSED   | 0.9853355    PASSED  | 0.5537414    PASSED  |
          sts_serial    7 | 0.9675717   PASSED   | 0.9675717    PASSED  | 0.3804243    PASSED  |
          sts_serial    7 | 0.4446567   PASSED   | 0.4446567    PASSED  | 0.0923678    PASSED  |
          sts_serial    8 | 0.7254384   PASSED   | 0.7254384    PASSED  | 0.4544030    PASSED  |
          sts_serial    8 | 0.8984816   PASSED   | 0.8984816    PASSED  | 0.7501155    PASSED  |
          sts_serial    9 | 0.8255134   PASSED   | 0.8255134    PASSED  | 0.4260288    PASSED  |
          sts_serial    9 | 0.6609663   PASSED   | 0.6609663    PASSED  | 0.5622308    PASSED  |
          sts_serial   10 | 0.9984397     **     | 0.9984397      **    | 0.5789212    PASSED  |
          sts_serial   10 | 0.7987434   PASSED   | 0.7987434    PASSED  | 0.8599317    PASSED  |
          sts_serial   11 | 0.5552886   PASSED   | 0.5552886    PASSED  | 0.3546752    PASSED  |
          sts_serial   11 | 0.4417852   PASSED   | 0.4417852    PASSED  | 0.5042245    PASSED  |
          sts_serial   12 | 0.3843880   PASSED   | 0.3843880    PASSED  | 0.6723639    PASSED  |
          sts_serial   12 | 0.1514682   PASSED   | 0.1514682    PASSED  | 0.9428701    PASSED  |
          sts_serial   13 | 0.5396454   PASSED   | 0.5396454    PASSED  | 0.5793677    PASSED  |
          sts_serial   13 | 0.9497671   PASSED   | 0.9497671    PASSED  | 0.3370774    PASSED  |
          sts_serial   14 | 0.3616613   PASSED   | 0.3616613    PASSED  | 0.4372343    PASSED  |
          sts_serial   14 | 0.3996251   PASSED   | 0.3996251    PASSED  | 0.5185021    PASSED  |
          sts_serial   15 | 0.3847188   PASSED   | 0.3847188    PASSED  | 0.3188851    PASSED  |
          sts_serial   15 | 0.1012968   PASSED   | 0.1012968    PASSED  | 0.1631942    PASSED  |
          sts_serial   16 | 0.9974802     **     | 0.9974802      **    | 0.6645914    PASSED  |
          sts_serial   16 | 0.1157822   PASSED   | 0.1157822    PASSED  | 0.3465564    PASSED  |
         rgb_bitdist    1 | 0.4705599   PASSED   | 0.4705599    PASSED  | 0.8627740    PASSED  |
         rgb_bitdist    2 | 0.7578920   PASSED   | 0.7578920    PASSED  | 0.3296790    PASSED  |
         rgb_bitdist    3 | 0.9934502   PASSED   | 0.9934502    PASSED  | 0.5558012    PASSED  |
         rgb_bitdist    4 | 0.3674201   PASSED   | 0.3674201    PASSED  | 0.1607977    PASSED  |
         rgb_bitdist    5 | 0.7930273   PASSED   | 0.7930273    PASSED  | 0.9999802      **    |
         rgb_bitdist    6 | 0.8491477   PASSED   | 0.8491477    PASSED  | 0.3774760    PASSED  |
         rgb_bitdist    7 | 0.1537432   PASSED   | 0.1537432    PASSED  | 0.4715169    PASSED  |
         rgb_bitdist    8 | 0.9454030   PASSED   | 0.9454030    PASSED  | 0.9890644    PASSED  |
         rgb_bitdist    9 | 0.2017856   PASSED   | 0.2017856    PASSED  | 0.0571014    PASSED  |
         rgb_bitdist   10 | 0.9989305     **     | 0.9989305      **    | 0.4575834    PASSED  |
         rgb_bitdist   11 | 0.4441883   PASSED   | 0.4441883    PASSED  | 0.4960057    PASSED  |
         rgb_bitdist   12 | 0.7074388   PASSED   | 0.7074388    PASSED  | 0.6808850    PASSED  |
rgb_minimum_distance    2 | 0.9604056   PASSED   | 0.9604056    PASSED  | 0.8859729    PASSED  |
rgb_minimum_distance    3 | 0.5143592   PASSED   | 0.5143592    PASSED  | 0.3266204    PASSED  |
rgb_minimum_distance    4 | 0.3779106   PASSED   | 0.3779106    PASSED  | 0.3537417    PASSED  |
rgb_minimum_distance    5 | 0.4861264   PASSED   | 0.4861264    PASSED  | 0.9032057    PASSED  |
    rgb_permutations    2 | 0.9206310   PASSED   | 0.9206310    PASSED  | 0.8052940    PASSED  |
    rgb_permutations    3 | 0.9299743   PASSED   | 0.9299743    PASSED  | 0.2209750    PASSED  |
    rgb_permutations    4 | 0.8330345   PASSED   | 0.8330345    PASSED  | 0.5819945    PASSED  |
    rgb_permutations    5 | 0.2708879   PASSED   | 0.2708879    PASSED  | 0.9276941    PASSED  |
      rgb_lagged_sum    0 | 0.0794660   PASSED   | 0.0794660    PASSED  | 0.9918681    PASSED  |
      rgb_lagged_sum    1 | 0.5279555   PASSED   | 0.5279555    PASSED  | 0.1304600    PASSED  |
      rgb_lagged_sum    2 | 0.0433872   PASSED   | 0.0433872    PASSED  | 0.1149961    PASSED  |
      rgb_lagged_sum    3 | 0.0028004     **     | 0.0028004      **    | 0.2731577    PASSED  |
      rgb_lagged_sum    4 | 0.0000074     **     | 0.0000074      **    | 0.8978870    PASSED  |
      rgb_lagged_sum    5 | 0.1332411   PASSED   | 0.1332411    PASSED  | 0.2065880    PASSED  |
      rgb_lagged_sum    6 | 0.0412128   PASSED   | 0.0412128    PASSED  | 0.7611867    PASSED  |
      rgb_lagged_sum    7 | 0.0225446   PASSED   | 0.0225446    PASSED  | 0.4810145    PASSED  |
      rgb_lagged_sum    8 | 0.0087433   PASSED   | 0.0087433    PASSED  | 0.3120378    PASSED  |
      rgb_lagged_sum    9 | 0.0000000   ******   | 0.0000000    ******  | 0.1334315    PASSED  |
      rgb_lagged_sum   10 | 0.4147842   PASSED   | 0.4147842    PASSED  | 0.2334790    PASSED  |
      rgb_lagged_sum   11 | 0.0206564   PASSED   | 0.0206564    PASSED  | 0.6491578    PASSED  |
      rgb_lagged_sum   12 | 0.0755835   PASSED   | 0.0755835    PASSED  | 0.5332069    PASSED  |
      rgb_lagged_sum   13 | 0.3112028   PASSED   | 0.3112028    PASSED  | 0.4194447    PASSED  |
      rgb_lagged_sum   14 | 0.0000000   ******   | 0.0000000    ******  | 0.2584573    PASSED  |
      rgb_lagged_sum   15 | 0.0890059   PASSED   | 0.0890059    PASSED  | 0.0007064      **    |
      rgb_lagged_sum   16 | 0.2962076   PASSED   | 0.2962076    PASSED  | 0.1344984    PASSED  |
      rgb_lagged_sum   17 | 0.2696070   PASSED   | 0.2696070    PASSED  | 0.2242021    PASSED  |
      rgb_lagged_sum   18 | 0.0826388   PASSED   | 0.0826388    PASSED  | 0.0450341    PASSED  |
      rgb_lagged_sum   19 | 0.0000000   ******   | 0.0000000    ******  | 0.5508302    PASSED  |
      rgb_lagged_sum   20 | 0.0101437   PASSED   | 0.0101437    PASSED  | 0.4290150    PASSED  |
      rgb_lagged_sum   21 | 0.1417859   PASSED   | 0.1417859    PASSED  | 0.1624411    PASSED  |
      rgb_lagged_sum   22 | 0.0160264   PASSED   | 0.0160264    PASSED  | 0.5204838    PASSED  |
      rgb_lagged_sum   23 | 0.0535167   PASSED   | 0.0535167    PASSED  | 0.6571892    PASSED  |
      rgb_lagged_sum   24 | 0.0000000   ******   | 0.0000000    ******  | 0.8578906    PASSED  |
      rgb_lagged_sum   25 | 0.8453426   PASSED   | 0.8453426    PASSED  | 0.3568988    PASSED  |
      rgb_lagged_sum   26 | 0.2113484   PASSED   | 0.2113484    PASSED  | 0.9755715    PASSED  |
      rgb_lagged_sum   27 | 0.1903762   PASSED   | 0.1903762    PASSED  | 0.4356739    PASSED  |
      rgb_lagged_sum   28 | 0.0733066   PASSED   | 0.0733066    PASSED  | 0.8354990    PASSED  |
      rgb_lagged_sum   29 | 0.0000000   ******   | 0.0000000    ******  | 0.1716599    PASSED  |
      rgb_lagged_sum   30 | 0.0932124   PASSED   | 0.0932124    PASSED  | 0.0732090    PASSED  |
      rgb_lagged_sum   31 | 0.0000000   ******   | 0.0000000    ******  | 0.3497910    PASSED  |
      rgb_lagged_sum   32 | 0.0843455   PASSED   | 0.0843455    PASSED  | 0.5441949    PASSED  |
     rgb_kstest_test    0 | 0.4399862   PASSED   | 0.4399862    PASSED  | 0.9766581    PASSED  |
     dab_bytedistrib    0 | 0.0748312   PASSED   | 0.0748312    PASSED  | 0.7035800    PASSED  |
             dab_dct  256 | 0.0919474   PASSED   | 0.0919474    PASSED  | 0.3985889    PASSED  |
        dab_filltree   32 | 0.1227533   PASSED   | 0.1227533    PASSED  | 0.7390925    PASSED  |
        dab_filltree   32 | 0.6819630   PASSED   | 0.6819630    PASSED  | 0.1773611    PASSED  |
       dab_filltree2    0 | 0.1774773   PASSED   | 0.1774773    PASSED  | 0.2088828    PASSED  |
       dab_filltree2    1 | 0.1718216   PASSED   | 0.1718216    PASSED  | 0.2257006    PASSED  |
        dab_monobit2   12 | 0.9999881     **     | 0.9999881      **    | 0.8084149    PASSED  |

 

 


 

Download DIEHARD_test_of_MAPLE_MersenneTwister.mw

A lot of supplementary details are given in the attached file.
I let the readers discover by themselves if Maple's implementation of the Mersenne Twister PRNG is correct or not.
Beyond this exercise, I hope this work will be useful to people who could be tempted to test their own generator.

 

 

After updating to Maple2019.2, the user initialization files maple.ini (Windows) and .mapleinit (Linux), which are placed into the user's home dierectories, are no longer read in upon startup. The behavior of Maple2019.1 was correct, the files were read.

We have just released updates to Maple and MapleSim.

Maple 2019.2 includes corrections and improvements to a variety of areas in the product, including a new “Go to page ____” option in print preview (that am personally quite pleased about), sections are expanded by default when printing or exporting, a fix to a problem using non-executable math with text in document mode that sometimes made it impossible to advance to a new line using Enter, improvements to VectorCalculus, select, abs and other math functions, support for macOS Catalina, and more.  We recommend that all Maple 2019 users install these updates.

This update is available through Tools>Check for Updates in Maple, and is also available from our website on the Maple 2019.2 download page, where you can also find more details.

For MapleSim users, the MapleSim 2019.2 family of products includes enhancements in the areas of model development and toolchain connectivity, including substantial enhancements to the MapleSim CAD toolbox.   For more details and download instructions, visit the MapleSim 2019.2 download page.

Although the graph of a parametrized surface can be viewed and manipulated on the computer screen as a surface in 3D, it is not quite suitable for printing on a 3D printer since such a surface has zero thickness, and thus it does not correspond to physical object.

To produce a 3D printout of a surface, it needs to be endowed with some "thickness".  To do that, we move every point from the surface in the direction of that point's nomral vector by the amount ±T/2, where T is the desired thickness.  The locus of the points thus obtained forms a thin shell of thickness T around the original surface, thus making it into a proper solid. The result then may be saved into a file in the STL format and be sent to a 3D printner for reproduction.

The worksheet attached to this post provides a facility for translating a parametrized surface into an STL file.  It also provides a command for viewing the thickened object on the screen.  The details are documented within that worksheet.

Here are a few samples.  Each sample is shown twice—one as it appears within Maple, and another as viewed by loading the STL file into MeshLab which is a free mesh viewing/manipulation software.

 

Here is the worksheet that produced these:  thicken.mw

 

 

I experienced a significant obstacle while trying to generate independent random samples with Statistics:-Sample on different nodes of a Grid multi-processing environment. After many hours of trial-and-error, I discovered an astonishing workaround, and I achieved excellent time and memory performance. Since this seems like a generally useful computation, I thought that it was worthy of a Post.

This Post is also worth reading to learn how to use Grid when you need to initialize a substantial environment on each node before using Grid:-Map or Grid:-Seq.

All remaining details are in the following worksheet.
 

How to use Statistics:-Sample in the `Grid` environment

Author: Carl Love <carl.j.love@gmail.com> 1 August 2019

 

I experienced a significant obstacle while trying to generate indenpendent random samples with Statistics:-Sample on the nodes of a multi-processor Grid (on a single computer). After several hours of trial-and-error, I discovered that two things are necessary to do this:

1. 

The random number generator needs to be seeded differently in each node. (The reason for this is easy to understand.)

2. 

The random variables generated by Statistics:-RandomVariable need to have different names in each node. This one is mind-boggling to me. Afterall, each node has its own kernel and so its own memory It's as if the names of random variables are stored in a disk file which all kernels access. And also the generator has been seeded differently in each node.

 

Once these things were done, the time and memory performance of the computation were excellent.

restart
:

Digits:= 15
:

#Specify the size of the computation:
(n1,n2,n3):= (100, 100, 1000):
# n1 = size of each random sample;
# n2 = number of samples in a batch;
# n3 = number of batches.

#
#Procedure to initialize needed globals on each node:
Init:= proc(n::posint)
local node:= Grid:-MyNode();
   #This is wrapped in parse so that it'll act globally. Otherwise, an environment
   #variable would be reset when this procedure ends.
   parse("Digits:= 15;", 'statement');

   randomize(randomize()+node); #Initialize independent RNG for this node.
   #If repeatability of results is desired, remove the inner randomize().

   (:-X,:-Y):= Array(1..n, 'datatype'= 'hfloat') $ 2;

   #Perhaps due to some oversight in the design of Statistics, it seems necessary that
   #r.v.s in different nodes **need different names** in order to be independent:
   N||node:= Statistics:-RandomVariable('Normal'(0,1));
   :-TRS:= (X::rtable)-> Statistics:-Sample(N||node, X);
   #To verify that different names are needed, change N||node to N in both lines.
   #Doing so, each node will generate identical samples!

   #Perform some computation. For the pedagogical purpose of this worksheet, all that
   #matters is that it's some numeric computation on some Arrays of random Samples.
   :-GG:= (X::Array, Y::Array)->
      evalhf(
         proc(X::Array, Y::Array, n::posint)
         local s, k, S:= 0, p:= 2*Pi;
            for k to n do
               s:= sin(p*X[k]);  
               S:= S + X[k]^2*cos(p*Y[k])/sqrt(2-sin(s)) + Y[k]^2*s
            od
         end proc
         (X, Y, n)
      )      
   ;
   #Perform a batch of the above computations, and somehow numerically consolidate the
   #results. Once again, pedagogically it doesn't matter how they're consolidated.  
   :-TRX1:= (n::posint)-> add(GG(TRS(X), TRS(Y)), 1..n);
   
   #It doesn't matter much what's returned. Returning `node` lets us verify that we're
   #actually running this on a grid.
   return node
end proc
:

The procedure Init above uses the :- syntax to set variables globally for each node. The variables set are X, Y, N||node, TRS, GG, and TRX1. Names constructed by concatenation, such as N||node, are always global, so :- isn't needed for those.

#
#Time the initialization:
st:= time[real]():
   #Send Init to each node, but don't run it yet:
   Grid:-Set(Init)
   ;
   #Run Init on each node:
   Nodes:= Grid:-Run(Init, [n1], 'wait');
time__init_Grid:= time[real]() - st;

Array(%id = 18446745861500764518)

1.109

The only purpose of array Nodes is that it lets us count the nodes, and it lets us verify that Grid:-MyNode() returned a different value on each node.

num_nodes:= numelems(Nodes);

8

#Time the actual execution:
st:= time[real]():
   R1:= [Grid:-Seq['tasksize'= iquo(n3, num_nodes)](TRX1(k), k= [n2 $ n3])]:
time__run_Grid:= time[real]() - st

4.440

#Just for comparison, run it sequentially:
st:= time[real]():
   Init(n1):
time__init_noGrid:= time[real]() - st;

st:= time[real]():
   R2:= [seq(TRX1(k), k= [n2 $ n3])]:
time__run_noGrid:= time[real]() - st;

0.16e-1

24.483

R1 and R2 will be different because different random numbers were used, but they should have similar histograms.

plots:-display(
   Statistics:-Histogram~(
      <R1 | R2>, #side-by-side plots
      'title'=~ <<"With Grid\n"> | <"Without Grid\n">>,
      'gridlines'= false
   )
);

(Plot output deleted because MaplePrimes cannot handle side-by-side plots!)

They look similar enough to me!

 

Let's try to quantify the benefit of using Grid:

speedup_factor:= time__run_noGrid / time__run_Grid;

5.36319824753560

Express that as a fraction of the theoretical maximum speedup:

efficiency:= speedup_factor / num_nodes;

.670399780941950

I think that that's really good!

 

The memory usage of this code is insignificant, which can be verified from an external memory monitor such as Winodws Task Manager. It's just a little bit more than that needed to start a kernel on each node. It's also possible to measure the memory usage programmatically. Doing so for a Grid:-Seq computation is a little bit beyond the scope of this worksheet.

 


 

Download GridRandSample.mw

Here are the histograms:

Hare in the forest

The rocket flies

  

Быльнов_raketa_letit.mws

1 2 3 4 5 6 7 Last Page 1 of 76