epostma

1534 Reputation

19 Badges

16 years, 255 days
Maplesoft

Social Networks and Content at Maplesoft.com

Maple Application Center
I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are replies submitted by epostma

Thanks Robert, I enjoyed reading this account.

Two notes, one mathematical and one Maple-related. First the mathematical one: I think you would expect 3 degrees of freedom, right? There are 10 equations, but there is a dependency between them: all four row-sums must needs sum to the same as all four column-sums. So there are really only 9 independent equations.

The Maple point is that the combinatorial question you raise is solved by the Iterator package, in particular the BoundedComposition command. You supply it with a k-tuple of bounds (b1, ..., bk), and a sum s, all nonnegative integers. It then finds all k-tuples (a1, ..., ak) of nonnegative integers with ai <= bi that sum to s. This raises one problem: we want integers from 1-9, and  BoundedComposition will also give us zeroes -- but that's easily dealt with: we just use integers from 0-8 instead, adjust the requested sum (by subtracting k=4), and presto. For example:

with(Iterator):
k := 4;
B := BoundedComposition([8 $ k], 25-k):
Number(B); # result: 324
Print(B, 'showrank');

tells us that there are 324 such tuples with sum 25, and will print all of them. To do a computation with them, you would use:

for b in BoundedComposition([8 $ k], 25 - k) do
  tup := b +~ 1;
  # Now use the tuple, tup, which will contain suitable integers 1-9. (It's an Array.)
end do;

Of course, if one number is given already, you could generate 3-tuples.

@C_R If this is what you want, the best option is probably to use convert/system to convert every unit to one consistent system, and then use Carl's solution below:

expr := sin(t * Unit(1/min)) * Unit(mm);
to_SI := evalindets(expr, specfunc(Units:-Unit), e -> convert(e, 'system', 'SI'));
result := eval(to_SI, Units:-Unit = 1);
# result is now sin(t/60) / 1000

@C_R For what it's worth, here is another teaser screenshot...

Hi everyone,

I think Tom and Carl have already answered the actual question, so I won't add an answer... but I will add a little teaser image of what the current in-development version of Maple looks like when applied to this worksheet:

If everything goes as expected, you can expect this to work without any tinkering being necessary, in the next major release of Maple.

@mike7944 You're welcome!

Hi Mike,

I should have mentioned explicitly - this was written to work in Maple 2021. Over the last couple of years we've made some additions to the Maple syntax. One of them is allowing full if statements (including arbitrary substatements) inside expressions; it looks like that's what you're running into. If you need to use this code in an older version, you can replace this fragment:

 local result := add(nonconst) + (
    if is(const = 0) then
      0;
    else
      const := evalf(const);
      if type(const, ':-float') then
        frem(const, 2.*Pi);
      else
        frem(Re(const), 2.*Pi) + I*Im(const);
      end if;
    end if);

with the following:

 local result := add(nonconst);
    if is(const <> 0) then
      const := evalf(const);
      if type(const, ':-float') then
        result := result + frem(const, 2.*Pi);
      else
        result := result + frem(Re(const), 2.*Pi) + I*Im(const);
      end if;
    end if;

You may also need to move all local declarations to the top of the procedure and module; for example, the module would start with a line saying local ModuleApply, process_trig; and the word local would be removed from the lines assigning to ModuleApply and process_trig.

Please let me know if you're still having trouble -- I don't have Maple 2017 handy so I can't easily check right now, but I'll dig it out if needed.

@Carl Love Yeah, the change should be in 2020.0 and later.

Best,
Erik.

@vs140580 Another option would be to use the command Iterator:-CartesianProduct. If you want to deal with the results one by one, you can use it like this:

m := 4; # or more
S := {0, 1, 2}; # or something else
for v in Iterator:-CartesianProduct(S $ m) do
  print(v); # v will be an Array with one result
end do;

Note that the Array in v will be reused between iterations; only its contents will change.

@Carl Love Thanks again, Carl, for all your time on this.

@epostma 

One more thing that I meant to add: I also looked up who introduced this bug. As is so often the case with these sorts of things, I found it was introduced in a change made more than a few years ago that rewrote some working code unnecessarily... submitted by one Erik Postma :(

@mmcdara Strictly speaking, the Gaussian family also has a specific sampler - it just happens to use the ziggurat algorithm, which could have been implemented as a generic algorithm, but isn't really. So I would say that all pre-defined distributions have a specific sampler implemented.

@Carl Love Thanks Carl - once I get around to fixing this, I'll be sure to re-run that worksheet. I expect they're all the same underlying bug; many of the distributions have their samples generated by sampling the normal distribution as one of the steps.

@mmcdara Ziggurat is used for the exponential and normal distributions.

For beta, it samples two gamma random variates.

For gamma, it looks a little tricky - I can't easily summarize it here. The general case involves generating a truncated normally distributed random variate and a uniform random variate and combining them in various ways... I'm sure if I looked into it for longer I'd recognize a commonly recommended algorithm but I haven't done that yet.

For Cauchy, it's interesting: if a hardware float value is requested (which is the common case, it usually happens when Digits <= 15), then there is a small procedure that uses 8 predefined constants (again - I haven't looked into it in detail but when I need to I'm sure I could recognize the algorithm). If a software float value is requested, those constants can't be pre-defined at arbitrarily high precision; it uses the inverse transformation method in that case.

Hi mmcdara,

Thanks for filing this report. I'll start investigating this shortly.

You ask about the Box-Mueller method. I can't say why it wasn't chosen (the person who implemented this code left Maplesoft long ago), but I'm not eager to switch over to it, if we can manage to fix this bug. The current method is a ziggurat algorithm (see e.g. https://en.wikipedia.org/wiki/Ziggurat_algorithm) which I think is a well-respected choice. The fallback algorithm for the tail appears to be the one that Marsaglia suggests, at least according to that wikipedia page.

 

Erik Postma
Manager of the mathematical software group
Maintainer of the Statistics package

I'd just like to add that we're looking into this issue; we've entered it into our internal systems. Thanks everyone for reporting and investigating it!

The previous post that _Maxim_ mentioned was removed exactly because it dealt with the same issue as this post, and this post had some comments and answers at that time whereas the other one didn't.

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