Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 55 days
Contact mmcdara

MaplePrimes Activity


These are replies submitted by mmcdara

@ecterrab Thank you Edgardo for thi...

mmcdara 7891
February 03 2019
0 0

@ecterrab 
Thank you Edgardo for this clarification.

I was myself impressed by Maple's performance (Mathematica's communication often suggests that it does things better and these 522 examples set things straight).

You wroye '...I haven't seen people complaining about this. People seem to find this design natural.... ': I find nothing wrong with that, especially since, not being a specialist, I cannot criticize your choices in a reasoned way.

Than you again

@Preben Alsholm  I wasn't awar...

mmcdara 7891
February 02 2019
0 0

@Preben Alsholm 

I wasn't aware of that detail.

In fact I was going through the list of PDEs (https://www.12000.org/my_notes/pde_in_CAS/pdse3.htm#x4-30003), looking for the ones that were solved only by one CAS (Mathematica or Maple).
Among them I was intrigued by the test case 19: Although it has been labelled "hand solved = No", it seemed at first sight that particular forms of solutions could easily be obtained.
These are of the form I giae in Problem-19.mw. I was then surprised that Maple only returned one single particular solution of this form.

"pdsolving" the equation without 'build' (the option used in the 12000.org site) delivers effectively more general solutions, even more general than the one I obtained by hand.

I think your remark gives an ending point tothe question
 

@tomleslie @Carl Love @acerFor...

mmcdara 7891
February 02 2019
0 0

@tomleslie 
@Carl Love 
@acer

For information alone, the file "test.mw" produces no error when executed with Maple 2015 (Mac OSX).

test-2_Maple2015.mw

 

 

I don't have Maple 2018 here, only 2...

mmcdara 7891
January 31 2019
0 0

I cn't use Maple 2017 right now, only 2015... but it works well with this latter

> 

restart:
> 

with(CodeGeneration):
> 

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

 (1)
> 

R(exp(x))

cg <- exp(x)

  
> 

 


 

How do you test randomness? Can you give...

mmcdara 7891
January 30 2019
0 0

How do you test randomness?
Can you give an example of such a sequence that would not be random?

@Jalale  You can of course use the...

mmcdara 7891
January 30 2019
0 0

@Jalale 

You can of course use the method proposed by dahr.

But it looks maybe "artificial" because it doesn't uses Histogram (which is no more than a sequence of PLOT(POLYGONS(...))

If You REALLY want to use Histogram you can do this 

restart:
A := [1, 2, 3, 4, 5, 6]:
B := [10, 20, 30, 40, 50, 60]:
C := [seq(seq(n, m=1..B[n]), n in A)]:
Statistics:-Histogram(C)

 

The fastest and the most memory efficien...

mmcdara 7891
January 27 2019
0 0

@Carl Love 

Based on the representation of the number of dice as a sum of powers of 2 in order to reduce the number of convolution products.
The limits in terms of number of dice are extended very far compared to the routine DiceSumDistr, and this for still reasonable calculation times (4.4s & 472 Mb for 65535 dice)

As shown in the final plot, routine Q3 has a run time which is linear in N (so is the memory used)

SumOfDice_3.mw
 

> 

restart:
> 

# Routine used as reference

Q1 := proc(R, M)

  # M is the number of dice

  local N, z, S, s, C, m:

  # initialization

  N := numelems(R):
  z := seq(0, k=1..N):
  S := copy(R):

  if M=1 then
    return [$M..M*N] =~ S /~ 6^M;
  end if:

  # Convolution product

  for m from 2 to M do
    s := [z, S[], z]:
    S := map(n -> add(R[k]*s[k+n], k=1..N), [$1..m*N-m+1]):
  end do:
  return [$M..M*N] =~ S /~ 6^M;
end proc

proc (R, M) local N, z, S, s, C, m; N := numelems(R); z := seq(0, k = 1 .. N); S := copy(R); if M = 1 then return `~`[:-`=`]([`$`(M .. M*N)], ` $`, `~`[:-`/`](S, ` $`, 6^M)) end if; for m from 2 to M do s := [z, S[], z]; S := map(proc (n) options operator, arrow; add(R[k]*s[k+n], k = 1 .. N) end proc, [`$`(1 .. m*N-m+1)]) end do; return `~`[:-`=`]([`$`(M .. M*N)], ` $`, `~`[:-`/`](S, ` $`, 6^M)) end proc

 (1)
> 

# Fastest routine
# (could it be improved?)

Q3 := proc(R, NumberOfDice)

  local M, N, N2, z, S, ConvTable, P, s, C, m:

  uses SignalProcessing:

  M         := floor(log[2](NumberOfDice));
  S         := copy(R):
  K         := 0:  # counter of the number of convolution products
  ConvTable := table([0=S^+]):
  for m from 1 to M do
    S            := Convolution(S, S);
    K            := K+1:
    ConvTable[m] := copy(S):
  end do:


  N2 := convert(NumberOfDice, base, 2);
  P  := [ListTools:-SearchAll(1, %)] -~ 1;
  S  := copy(ConvTable[P[1]]);

  for m in P[2..-1] do
    S := Convolution(S, ConvTable[m]);
    K := K+1:
  end do:

  printf("Number of convolution products = %d\n", K);
  convert~(S, rational) / 6^NumberOfDice;
  return %;
end proc

Warning, `K` is implicitly declared local to procedure `Q3`

  

proc (R, NumberOfDice) local M, N, N2, z, S, ConvTable, P, s, C, m, K; M := floor(log[2](NumberOfDice)); S := copy(R); K := 0; ConvTable := table([0 = S^%T]); for m to M do S := SignalProcessing:-Convolution(S, S); K := K+1; ConvTable[m] := copy(S) end do; N2 := convert(NumberOfDice, base, 2); P := `~`[:-`-`]([ListTools:-SearchAll(1, %)], ` $`, 1); S := copy(ConvTable[P[1]]); for m in P[2 .. -1] do S := SignalProcessing:-Convolution(S, ConvTable[m]); K := K+1 end do; printf("Number of convolution products = %d ", K); `~`[convert](S, rational)/6^NumberOfDice; return % end proc

 (2)
> 

# verification

S := Q3(<seq(1, k=1..6)>, 5):
T := rhs~(Q1([seq(1, k=1..6)], 5)):
<S^+ | Vector[column](T)>;

Number of convolution products = 3

  

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

 (3)
> 

# Carls routine

DiceSumDistr:= (D::list(list(integer)))->
   table((degree=icontent)~([op(expand(mul((add@`~`[`^`])~(_x,D))))])) /~ mul(nops~(D))
;

proc (D::(list(list(integer)))) options operator, arrow; `~`[:-`/`](table(`~`[degree = icontent]([op(expand(mul(`~`[`@`(add, `~`[`^`])](_x, D))))])), ` $`, mul(`~`[nops](D))) end proc

 (4)
> 

# time & memory comparisons

CodeTools:-Usage(Q3(<seq(1, k=1..6)>, 4096)):
CodeTools:-Usage(DiceSumDistr([[$1..6] $ 4096])):

Number of convolution products = 12

memory used=29.79MiB, alloc change=0.88MiB, cpu time=291.00ms, real time=291.00ms, gc time=19.74ms

memory used=138.11MiB, alloc change=108.01MiB, cpu time=2.25s, real time=2.21s, gc time=126.25ms

  
> 

CodeTools:-Usage(Q3(<seq(1, k=1..6)>, 4095)):
CodeTools:-Usage(DiceSumDistr([[$1..6] $ 4095])):

Number of convolution products = 22

memory used=29.71MiB, alloc change=0 bytes, cpu time=273.00ms, real time=273.00ms, gc time=0ns

memory used=137.45MiB, alloc change=32.00MiB, cpu time=2.20s, real time=2.11s, gc time=145.06ms

  
> 

CodeTools:-Usage(Q3(<seq(1, k=1..6)>, 8192)):
CodeTools:-Usage(DiceSumDistr([[$1..6] $ 8192])):

Number of convolution products = 13

memory used=59.49MiB, alloc change=0 bytes, cpu time=510.00ms, real time=510.00ms, gc time=0ns

memory used=461.05MiB, alloc change=119.23MiB, cpu time=9.66s, real time=9.29s, gc time=583.85ms

  
> 

CodeTools:-Usage(Q3(<seq(1, k=1..6)>, 16384)):
CodeTools:-Usage(DiceSumDistr([[$1..6] $ 16384])):

Number of convolution products = 14

memory used=118.87MiB, alloc change=0 bytes, cpu time=1.05s, real time=1.05s, gc time=0ns

memory used=1.63GiB, alloc change=0.64GiB, cpu time=47.79s, real time=47.08s, gc time=1.22s

  
> 

CPU := NULL:
for k from 1 to 20 do
  CPU := CPU, [2^k, CodeTools:-Usage(Q3(<seq(1, k=1..6)>, 2^k), 'output'='cputime')]:
end do:

Number of convolution products = 1
memory used=58.24KiB, alloc change=0 bytes, cpu time=1000.00us, real time=1000.00us, gc time=0ns
Number of convolution products = 2
memory used=76.75KiB, alloc change=0 bytes, cpu time=1000.00us, real time=1000.00us, gc time=0ns
Number of convolution products = 3
memory used=115.42KiB, alloc change=0 bytes, cpu time=2.00ms, real time=2.00ms, gc time=0ns
Number of convolution products = 4
memory used=172.12KiB, alloc change=0 bytes, cpu time=2.00ms, real time=2.00ms, gc time=0ns
Number of convolution products = 5
memory used=340.09KiB, alloc change=0 bytes, cpu time=3.00ms, real time=3.00ms, gc time=0ns
Number of convolution products = 6
memory used=0.63MiB, alloc change=0 bytes, cpu time=6.00ms, real time=6.00ms, gc time=0ns
Number of convolution products = 7
memory used=1.28MiB, alloc change=0 bytes, cpu time=14.00ms, real time=13.00ms, gc time=0ns
Number of convolution products = 8
memory used=2.74MiB, alloc change=0 bytes, cpu time=27.00ms, real time=27.00ms, gc time=0ns
Number of convolution products = 9
memory used=3.79MiB, alloc change=0 bytes, cpu time=34.00ms, real time=33.00ms, gc time=0ns
Number of convolution products = 10
memory used=7.50MiB, alloc change=0 bytes, cpu time=65.00ms, real time=65.00ms, gc time=0ns
Number of convolution products = 11

memory used=14.94MiB, alloc change=0 bytes, cpu time=132.00ms, real time=132.00ms, gc time=0ns
Number of convolution products = 12

memory used=29.77MiB, alloc change=-15.03MiB, cpu time=556.00ms, real time=363.00ms, gc time=276.31ms
Number of convolution products = 13

memory used=59.50MiB, alloc change=0 bytes, cpu time=523.00ms, real time=523.00ms, gc time=0ns
Number of convolution products = 14

memory used=118.87MiB, alloc change=0 bytes, cpu time=1.06s, real time=1.06s, gc time=0ns

Number of convolution products = 15

memory used=237.62MiB, alloc change=6.27MiB, cpu time=2.13s, real time=2.13s, gc time=0ns

Number of convolution products = 16

memory used=475.14MiB, alloc change=8.76MiB, cpu time=4.66s, real time=4.49s, gc time=258.60ms

Number of convolution products = 17

memory used=0.93GiB, alloc change=11.25MiB, cpu time=10.28s, real time=9.44s, gc time=1.34s

Number of convolution products = 18

memory used=1.86GiB, alloc change=27.50MiB, cpu time=19.96s, real time=18.15s, gc time=2.78s

Number of convolution products = 19

memory used=3.71GiB, alloc change=55.00MiB, cpu time=40.04s, real time=36.09s, gc time=5.90s

Number of convolution products = 20

memory used=7.42GiB, alloc change=110.00MiB, cpu time=82.91s, real time=71.75s, gc time=15.73s

  
> 

plots:-loglogplot([CPU], labels=["Number of dice", "cpu time"], labeldirections=[default, vertical])

 
> 

 


 

Download SumOfDice_3.mw

 

An already even faster routine when N=2^...

mmcdara 7891
January 27 2019
0 0

@Carl Love 

This post has been deleted and replaced by the most complete ("The fastest and the most memory efficient for...") below

@Carl Love  Impressive, indeed&nbs...

mmcdara 7891
January 27 2019
0 0

@Carl Love 

Impressive, indeed very fast while being extremely concise.
Versatility is also a strong point.

My only reservation is that you need a solid knowledge of the Maple language to understand what it does (I'm analyzing it step by step).

Comparing your procedure to mine I observed that I can easily reduce the run time of the latter by a factor 4 (by avoiding to manipulate rationals).
In my procedure, if N is the number of dice, I realise N convolution products of increasing length. The total number of operations is of the order of N^2 (up to some multiplicative constant). I had in mind to reduce this number by writting N 
as a sum n1 + n2 + ...np of positive integers (for instance 100 = 64 + 32 + 4) and using the recursive apporach used in the Fast Fourier Transform.
My feeling is that the complexity of the problem could be something like N*log[2](N).
If I have some time I will try to do something around this.

But it seems thay your own procedure has set the bar very high.
Bravo

@Kitonum @JalaleOur previous exchan...

mmcdara 7891
January 27 2019
0 0

@Kitonum @Jalale

Our previous exchange was inspiring. 
I give you what is perhaps the fastest version to calculate the probabilities of the sums of an arbitrary number of dice.

For 100 dice the CPU time is about 0.4 s
(to be compared to the dozen of seconds my "P" procedure and your "Composition" procedure require for 10 dice)

I think this new procedure could still be improved by using mutable structures and, more surely, by using the decomposition of the number of dice in powers of 2.

SumOfDice.mw
 

> 

restart:
> 

# Procedure Q is based on the following result:
#   the pdf of the sum of two RV X and Y is the convolution product
#   of pdf(X) by pdf(Y)
#

Q := proc(M)

  # M is the number of dice

  local R, N, z, S, s, C, m:

  # initialization

  R := [seq(1/6, k=1..6)];  #mass density function of a single dice
  N := numelems(R):
  z := seq(0, k=1..N):
  S := copy(R):

  if M=1 then
    return [$M..M*N] =~ S;
  end if:

  # Convolution product

  for m from 2 to M do
    s := [z, S[], z]:
    S := map(n -> add(R[k]*s[k+n], k=1..N), [$1..m*N-m+1]):
    C := [$m..m*N] =~ S;
  end do:
  return [$M..M*N] =~ S;
end proc

proc (M) local R, N, z, S, s, C, m; R := [seq(1/6, k = 1 .. 6)]; N := numelems(R); z := seq(0, k = 1 .. N); S := copy(R); if M = 1 then return `~`[:-`=`]([`$`(M .. M*N)], ` $`, S) end if; for m from 2 to M do s := [z, S[], z]; S := map(proc (n) options operator, arrow; add(R[k]*s[k+n], k = 1 .. N) end proc, [`$`(1 .. m*N-m+1)]); C := `~`[:-`=`]([`$`(m .. m*N)], ` $`, S) end do; return `~`[:-`=`]([`$`(M .. M*N)], ` $`, S) end proc

 (1)
> 

CodeTools:-Usage(Q(4))

memory used=55.55KiB, alloc change=0 bytes, cpu time=1000.00us, real time=0ns, gc time=0ns

  

[4 = 1/1296, 5 = 1/324, 6 = 5/648, 7 = 5/324, 8 = 35/1296, 9 = 7/162, 10 = 5/81, 11 = 13/162, 12 = 125/1296, 13 = 35/324, 14 = 73/648, 15 = 35/324, 16 = 125/1296, 17 = 13/162, 18 = 5/81, 19 = 7/162, 20 = 35/1296, 21 = 5/324, 22 = 5/648, 23 = 1/324, 24 = 1/1296]

 (2)
> 

CodeTools:-Usage(Q(10)):

memory used=338.62KiB, alloc change=0 bytes, cpu time=3.00ms, real time=2.00ms, gc time=0ns

  
> 

CodeTools:-Usage(Q(100)):

memory used=87.20MiB, alloc change=32.00MiB, cpu time=396.00ms, real time=432.00ms, gc time=10.96ms

  
> 

 


 

Download SumOfDice.mw

 

@Kitonum Agree!I misspoke, I meant ...

mmcdara 7891
January 26 2019
0 0

@Kitonum 

Agree!

I misspoke, I meant "shorter in number ofcommand lines", I should have wrote "more concise".
(The fact that "Composition" was written several years ago does not take away from the fact that "Composition" contains more characters than "P").

I did not look at the computation times but I was aware, when running my procedure "P", that "P" was rather slow. 
So you are undoubtly right by writing that "Composition" is 50 times faster than "P" in the case n=3. 
But is it always the case?

  • For n=8 I got   17s for "Composition" and 27s for "P".
  • For n=9 I got 113s for "Composition" and 45s for "P"

So arguing about cpu time issues is a waste of time.

I keep thinking "P" is simpler to understand to someone who is used to use the package Statistics, and at the same time I have no doubt that "Composition" presents a lot of interesting features on its side (if only to be much faster for small values of n)


My only purpose was to show how the problem you mentioned (which is not Jalale's) could be solved with the tools used by Jalale (Statistics package) ... even if this is not computationally optimal.
 

@Kitonum  Nothing new compared to ...

mmcdara 7891
January 26 2019
0 0

@Kitonum 

Nothing new compared to what you have done. But if all is needed are the  probabilities of the different sums, here is a shorter (maybe not faster) way to proceed

SumOfDice.mw

@vv  This is indeed a solution... ...

mmcdara 7891
January 26 2019
0 0

@vv 

This is indeed a solution...
 

@Carl Love  An astute way to adres...

mmcdara 7891
January 26 2019
0 0

@Carl Love 

An astute way to adress the problem!
Thanks Carl

@acer  I tried to force Maple to re...

mmcdara 7891
January 26 2019
0 0

@acer 
I tried to force Maple to return { b=c if a <> 0, undefined otherwise}.
I hoped that assume will do the job, but all my attempts failed.

I'm curious, is there a way to obtain this result ?

First 134 135 136 137 138 139 140 Last Page 136 of 154