Carl Love

Carl Love

28050 Reputation

25 Badges

12 years, 335 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@Kitonum But this is not even close to the OP's proposed result.

@sand15athome The following may be already totally obvious to you. But, just in case you don't know...

PLOT isn't a command of the Maple language; it's just an inert data structure that's interpreted by the GUI. (So, perhaps it can be considered a command to the GUI, but to the Maple kernel it's inert.) All commands that produce static (i.e., non-animated) 2-D plots have this data structure as their output. As long as you're aware of all that, I think that it's fine to write your own procedures that produce PLOT output; indeed, I often find it necessary.

@TomM Good work investigating this.

I made your ideas into a procedure. I use randomly selected complex evaluation points selected over a square centered at the origin of user-specified dimension (defaults to 4x4).

restart:

DetUnifloat:= proc(
     M::Matrix(square, polynom(float)),
     {radius::positive:= 2},
     {discardthreshold::positive:= Float(1, 3-Digits)}
)
description
     "A replacement for LinearAlgebra:-Determinant(..., method= unifloat)"
;
local
     R:= rand(evalf(-radius..radius)),
     x:= indets(M, And(name, Not(constant))),
     n:= op([1,1], M),
     maxdeg, mindeg, i, j, E, P, MC;
;
     if nops(x) > 1 then
          error "Only one variable allowed"
     end if;
     x:= `if`(nops(x) > 0, x[], 'x');
     maxdeg:= min(
          add(max(degree~(M[i,..])), i= 1..n),
          add(max(degree~(M[..,j])), j= 1..n)
     );
     mindeg:= max(
          add(min(ldegree~(M[i,..])), i= 1..n),
          add(min(ldegree~(M[..,j])), j= 1..n)
     );
     E:= ['R()' + 'R()'*I $ maxdeg+1];
     P:= CurveFitting:-PolynomialInterpolation(
          E, [seq](LinearAlgebra:-Determinant(eval(M, x= e)), e= E), x
     );
     MC:= max(abs~([coeffs](P,x)));
     select(
          t-> degree(t,x) >= mindeg,
          sort(MC*simplify(fnormal(P/MC, Digits, discardthreshold), 'zero'))
     )
end proc:
     

data_deg_ex := Matrix(6, 6, [[180*sigma2^2, -17980.4676072929*sigma2, 60*sigma2^2, -1792.74593023400*sigma2, -597.004097753022*sigma2, -60*sigma2], [-17980.4676072929*sigma2, 60*sigma2^2, -17980.4676072929*sigma2, -597.004097753022*sigma2, -597.581976744666*sigma2, 5993.48920243096], [60*sigma2^2, -17980.4676072929*sigma2, 180*sigma2^2, -597.581976744666*sigma2, -1791.01229325907*sigma2, -60*sigma2], [-1792.74593023400*sigma2, -597.004097753022*sigma2, -597.581976744666*sigma2, -60*sigma2, 5993.48920243096, 597.581976744666], [-597.004097753022*sigma2, -597.581976744666*sigma2, -1791.01229325907*sigma2, 5993.48920243096, -60*sigma2, 597.004097753022], [-60*sigma2, 5993.48920243096, -60*sigma2, 597.581976744666, 597.004097753022, 60]]);

data_deg_ex := Matrix(6, 6, {(1, 1) = 180*sigma2^2, (1, 2) = -17980.4676072929*sigma2, (1, 3) = 60*sigma2^2, (1, 4) = -1792.74593023400*sigma2, (1, 5) = -597.004097753022*sigma2, (1, 6) = -60*sigma2, (2, 1) = -17980.4676072929*sigma2, (2, 2) = 60*sigma2^2, (2, 3) = -17980.4676072929*sigma2, (2, 4) = -597.004097753022*sigma2, (2, 5) = -597.581976744666*sigma2, (2, 6) = 5993.48920243096, (3, 1) = 60*sigma2^2, (3, 2) = -17980.4676072929*sigma2, (3, 3) = 180*sigma2^2, (3, 4) = -597.581976744666*sigma2, (3, 5) = -1791.01229325907*sigma2, (3, 6) = -60*sigma2, (4, 1) = -1792.74593023400*sigma2, (4, 2) = -597.004097753022*sigma2, (4, 3) = -597.581976744666*sigma2, (4, 4) = -60*sigma2, (4, 5) = 5993.48920243096, (4, 6) = 597.581976744666, (5, 1) = -597.004097753022*sigma2, (5, 2) = -597.581976744666*sigma2, (5, 3) = -1791.01229325907*sigma2, (5, 4) = 5993.48920243096, (5, 5) = -60*sigma2, (5, 6) = 597.004097753022, (6, 1) = -60*sigma2, (6, 2) = 5993.48920243096, (6, 3) = -60*sigma2, (6, 4) = 597.581976744666, (6, 5) = 597.004097753022, (6, 6) = 60})

DetUnifloat(evalf(data_deg_ex));

HFloat(1.8662406249298306e11)*sigma2^8+HFloat(1.4795511563262647e14)*sigma2^7+HFloat(4.0379217068100504e16)*sigma2^6+HFloat(4.393918722504788e18)*sigma2^5+HFloat(1.6722907322306802e20)*sigma2^4

LinearAlgebra:-Determinant(data_deg_ex, method= minor);

186624000000*sigma2^8+147955115634640.*sigma2^7+0.403792170686782e17*sigma2^6+0.439391872250447e19*sigma2^5+0.167229073223072e21*sigma2^4

 


Download unifloat.mw

@tomleslie I think that you may be right---that unlike dsolve, pdsolve's solution procedures don't contain the (correct) partial derivatives. However, Maple says that the requested derivative is 0---it doesn't produce an empty plot or give an error message. It also gives 0 if you use the pds:-value method.

@vv And the value of that Stieltjes integral is not 0?

Can you give an example of a function which can be represented in Maple or any computer-algebra system for which the Stieltjes integral is different from the "normal" integral?

@tomleslie And your formula works for any base if you change the 10 and use a log-based formula for the length. I'm almost sure that this would be faster than the non-base-10 digit-reversal formula that I posted above, regardless of the length; but I don't feel like testing it right now.

In the meantime, here's my any-base digit-reversal:

revNumB:= proc(x::nonnegint, b::And(posint, Not(1)))
local y:= x, r:= irem(y, b, 'y');
     while y <> 0 do
          r:= b*r + irem(y, b, 'y')
     end do
end proc:

@tomleslie The [1,1] in your op[1,1] does nothing. The extreme liberalness of Maple syntax allows you to put a meaningless index on a function that doesn't expect an index. You could change [1,1] to [foo] for example.

@vv You wrote:

P.S. Do you have an example where a local for the variable in seq is mandatory? (probably only for multiple kernels [or threads?] in some circumstances).

The status of the bound variable of seq, add, and mul is the same. Here are two examples of what can go wrong:

restart:

RL:= ['rand(1..1000)()' $ 2^10]:

S1:= Threads:-Map(proc(x)          add(k, k= 1..x) end proc, RL):
S2:= Threads:-Map(proc(x) local k; add(k, k= 1..x) end proc, RL):     
S0:= map(x-> add(k, k= 1..x), RL):     

evalb~([S0=S1, S0=S2]);

[false, true]

k;

926

protect(n);

seq(n, n= 1..9);

Error, attempting to assign to `n` which is protected.  Try declaring `local n`; see ?protect for details.

 

Download local.mw

AFAIK, these are the only examples.

Regarding op(x) vs. x[]:

restart:
L:= [['rand()'] $ 2^23]:
gc(): CodeTools:-Usage([seq(x[], x= L)]):
memory used=64.00MiB, alloc change=64.00MiB, cpu time=1000.00ms, real time=1000.00ms, gc time=0ns
gc(): CodeTools:-Usage([seq(op(x), x= L)]):
memory used=64.00MiB, alloc change=64.00MiB, cpu time=1.72s, real time=1.72s, gc time=0ns
gc(): CodeTools:-Usage(map(x-> x[], L)):
memory used=0.50GiB, alloc change=64.00MiB, cpu time=18.05s, real time=7.82s, gc time=12.14s
gc(): CodeTools:-Usage(map(op, L)):
memory used=0.50GiB, alloc change=64.00MiB, cpu time=14.42s, real time=6.82s, gc time=9.80s
gc(): CodeTools:-Usage(op~(L)):
memory used=384.00MiB, alloc change=64.00MiB, cpu time=12.50s, real time=5.66s, gc time=9.12s

But the most-important conclusion from the above is that seq is much faster than any form of mapping. This is very unfortunate because seq has the most awkward syntax.

@smo4291990 The procs and do loops that you've shown so far are not quite correct. They're almost correct, but computers don't do "almost" (yet). So let's perfect your procs and do loops first. Please write a proc with a do loop that does anything of your choice. Use some reasonable style of indentation; don't write it inline (i.e., as part of a sentence). State what the procedure does.

If you follow my instructions, I promise that I'll work through this exercise with you.

@vv

Vote up.

I was able to improve your procedure, making it

  • take less time and memory,
  • handle the degenerate cases k > n and k = 0, and
  • more robust by doing type checking and by using a local variable.

Here's the code:

#VV's method, modified by Carl
rcomb2:= proc(n::nonnegint, k::nonnegint)
local i;
     if k > n then ()
     elif k=0 then []
     else
          proc(n,k,p:= [])
               if k=n then [seq(1..n), p[]]
               elif k=1 then seq([i,p[]], i= 1..n)
               else thisproc(n-1, k, p), thisproc(n-1, k-1, [n,p[]])
               fi
          end proc(n,k)
     fi
end proc:

With these modifications, the above is a tiny bit faster than my seq method, but uses more memory. I modified my Answer above to include both your original version and my modification in the overall comparison.

 

@vv Your rcomb seems to produce the correct number of tuples, but most of them are shorter than they should be:

rcomb(5,3);

     [1, 2, 3], [1, 2], [1, 3], [2, 3], [1, 2], [1, 3], [2, 3], [1, 4], [2, 4], [3, 4]

@Joe Riel Since I thought about it so long and hard, please let me know if you know of any way that uses the OP's simple recursion, but using tables or rtables rather than sets or lists.

@sand15athome 

So, I see that you figured out how to use ||, which allows you to combine some features of both indices and subscripts. In general, there are a few drawbacks to using ||, but they don't apply to your situation. The two drawbacks that I can think of to using || are

  • variables created with || are always global, and
  • if you use hundreds or thousands of such variables, you will notice a slowdown that won't happen if you use indexed variables.

Here is your worksheet with my ideas for improvement:

restart:

with(Statistics):
N := 2:
for n in ['n', $1..N] do
#     ^^^^^^^^^^^^^^^
   X__E__||n := RandomVariable(Normal(mu__E__||n, sigma__E__||n));
end do:

PDF(X__E__1, x__E__1);  # good, I wanted this

(1/2)*2^(1/2)*exp(-(1/2)*(x__E__1-mu__E__1)^2/sigma__E__1^2)/(Pi^(1/2)*sigma__E__1)

n := 'n':
PDF(X__E__||n, x__E__||n); # I thought to obtain the same output as befor with subscript n instead of 1

(1/2)*2^(1/2)*exp(-(1/2)*(x__E__n-mu__E__n)^2/sigma__E__n^2)/(Pi^(1/2)*sigma__E__n)

Product(PDF(X__E__||n, x__E__||n), n=1..N); # which does not work due to the previous "problem"

Product((1/2)*2^(1/2)*exp(-(1/2)*(x__E__n-mu__E__n)^2/sigma__E__n^2)/(Pi^(1/2)*sigma__E__n), n = 1 .. 2)

mul(PDF(X__E__||n, x__E__||n), n=1..N);   # mul instead of `*` as suggested by Carl
#I suggested that you use `mul` instead of `product`, not instead of `*`.

(1/2)*exp(-(1/2)*(x__E__1-mu__E__1)^2/sigma__E__1^2)*exp(-(1/2)*(x__E__2-mu__E__2)^2/sigma__E__2^2)/(Pi*sigma__E__1*sigma__E__2)

# but what this does not give the same result ?

product('PDF(X__E__||n, x__E__||n)', n=1..N);
#       ^                         ^
#Quotes needed to prevent premature evaluation; but, better to use `mul` instead.

(1/2)*exp(-(1/2)*(x__E__1-mu__E__1)^2/sigma__E__1^2)*exp(-(1/2)*(x__E__2-mu__E__2)^2/sigma__E__2^2)/(Pi*sigma__E__1*sigma__E__2)

 


Download WhatHappensHere3.mw

 

@vv

With Maple, modifying the code for experimentation is easy. Maple's arbitrary-length integer arithmetic is quite fast; it uses GMP. Maple's memory management isn't the greatest, but at least you get automatic garbage collection on parallel processors. I suspect that things like reversing the digits of an arbitary-length integer in an arbitrary base are very tedious in C.

If I was going to write a program intended to run for weeks or months operating on numbers with millions of digits, I'd use C. But I'd definitely write that program first in Maple, perfect it, and then hand translate it to C.

First 379 380 381 382 383 384 385 Last Page 381 of 709