epostma

1364 Reputation

17 Badges

12 years, 305 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

Hi MRI_user,

Glad that this works for you. Just a quick note - with this definition, MySecondFunction and MyMatrix are functionally equivalent. For any values u, v, w, we have that MySecondFunction(u, v, w) and MyMatrix(u, v, w) return identical matrices.

Erik.

Hi MRI_user,

Glad that this works for you. Just a quick note - with this definition, MySecondFunction and MyMatrix are functionally equivalent. For any values u, v, w, we have that MySecondFunction(u, v, w) and MyMatrix(u, v, w) return identical matrices.

Erik.

@Robert Israel : no, the idea would be to plot a formula that happens to render as the program plotting the formula, in the same spirit as Tupper's formula rendering as itself.

The item that makes Tupper's formula a bit of a cheat to me is that n, the number actually encoding the formula, is not part of what's displayed. With this formulation you'd fix that.

I guess one way of constructing the bitmap n in the program would be by using a font and a representation of the program, but I have a, probably unreasonable, hope that it might be possible to somehow construct that number through different means. But as you say, it's probably way too much work given the payoff.

Erik.

@Robert Israel : nice. In TTY maple, interestingly, some of the lines get "swallowed up" by other lines printed over top of them, but you can still read the lines that remain. In GUI it works beautifully of course.

I guess the more difficult challenge is to create a program that does the same without the use of (2D or 3D) TEXT plot structures...

Erik.

Hi alex_01,

I'm not sure I understand what you mean - I have a hunch but let me ask you to confirm it before I start looking at this in detail. Do you want something like a conditional correlation between X and Y given that X = X0, and then compare  that number for values of X0 that are in the tails versus near the median / mean / ... ?

Also, what assumptions can we make on (X, Y)? Can we assume they're in a bivariate normal distribution, or can we assume nothing at all about them, or something in between?

Erik Postma
Maplesoft. 

Does anybody see a way to write a Maple program that, when run, generates a graph that gives its own listing? Self-printing programs or quines (wikipedia) are well-known, but I haven't seen self-graphing ones before.

This could be a skeleton of the program:

n:=9609393799189588849716729621278527547150043396601293066515055192717028023952664246896428421743\
507181212671537827706233559932372808741443078913259639413377234878577357498239266297155171737169951\
652328905382216124032388558661840132355851360488286933379024914542292886670810961844960917051834540\
678277315517054053816273809676025656250169814820834187831638491155902256100036523513703438744618483\
787372381982248498634650331594100549747005931383392264972494617515457283667023697454610146559979337\
98537483143786841806593422227898388722980000748404719:
h:=17:
w:=106:
tupper:=(x,y)->1-floor(frac((floor(y/h)*2^(h*x-(y mod h)))/2)*2):
with(ImageTools):
View(Create(Array(0..h-1,-w..-1,(y,x)->tupper(x,y+n),datatype=float[8],order=C_order))):

Now we'd need to turn this textual representation into an image (using some font), read it in, find the corresponding numbers n, h, w as Robert explains, and find a fixed point of this map. However, it seems unlikely that this will work straightaway - as the number n grows, the needed screen real estate to display n grows as well, and probably more quickly than n itself. So we would need a way to specify n more compactly - using fewer pixels.

Any takers?

Erik Postma
Maplesoft.

Great answer, @longrob. One note: if you need to omit a row, it's even easier to just use

plot(N[2.., 2], N[2.., 3]);

Hope this helps,

Erik Postma
Maplesoft.

Great answer, @longrob. One note: if you need to omit a row, it's even easier to just use

plot(N[2.., 2], N[2.., 3]);

Hope this helps,

Erik Postma
Maplesoft.

You can also just do A[3] to obtain the third row.

Erik Postma
Maplesoft. 

You can also just do A[3] to obtain the third row.

Erik Postma
Maplesoft. 

@longrob : Sorry, I don't really have much windows experience. But I suggest you contact our support team, support@maplesoft.com. They deal with this sort of stuff on a regular basis. (If you're not in Canada or the US you could instead go to http://www.maplesoft.com/contact/international/ to find your local reseller.)

Best of luck,

Erik.

One case in which I could imagine getting a return value of type function would be if fsolve is unsuccessful in finding a root at full precision; it returns unevaluated in that case. (It does this because later manipulation might change the expression inside, and in that case it might be possible to find a solution when you re-evaluate.) A silly example:

fsolve(abs(x) = -1, x);

does some computation, sees that it can't find a root, and returns the original function call. Something similar can happen with roots that are numerically difficult to find with high precision. So, I echo @longrob's suggestion: it might be a good idea if you could show us exactly what fsolve gets as its arguments.

(As an aside - the opening bracket and closing parenthesis don't match. That must have been a copy-and-paste error.)

Hope this helps,

Erik Postma
Maplesoft.

One case in which I could imagine getting a return value of type function would be if fsolve is unsuccessful in finding a root at full precision; it returns unevaluated in that case. (It does this because later manipulation might change the expression inside, and in that case it might be possible to find a solution when you re-evaluate.) A silly example:

fsolve(abs(x) = -1, x);

does some computation, sees that it can't find a root, and returns the original function call. Something similar can happen with roots that are numerically difficult to find with high precision. So, I echo @longrob's suggestion: it might be a good idea if you could show us exactly what fsolve gets as its arguments.

(As an aside - the opening bracket and closing parenthesis don't match. That must have been a copy-and-paste error.)

Hope this helps,

Erik Postma
Maplesoft.

@longrob : Do any other commands involving Compiler:-Compile work? (See e.g. the ?Compiler/Compile help page.) If not, please look at the instructions at http://www.maplesoft.com/support/faqs/detail.aspx?sid=99396. On OS X, you'd have to install XCode; on 64-bit windows, you'll want to follow these instructions. The 32-bit windows installation comes with a working compiler and on Linux a compiler is usually already installed in the system and picked up automatically.

Erik.

@longrob : I'll answer your questions in random order.

  • Determining the time spent in different parts of a program is called profiling, and in Maple this functionality is provided by the  CodeTools:-Profiling  package - in particular, the CodeTools:-Profiling:-Profile and showstat commands. I think the help pages probably have enough detail to get you started (now that you know where to look! :) ).
  • The difference between versions 2 and 3 are twofold: first of all, the RandomVariables are created only once (creating such a random variable is not a big computation, but the administrative overhead is substantial) and then remembered. This is accomplished in the variableStorage Vector and the getVariables procedure. Secondly, using the profiling mentioned in the previous bullet point, I saw that the line defining expr also took substantial time. This could really only come from the sqrt function call. By enclosing sqrt in forward quotes, I prevented Maple from doing any actual computation there; when the expression enters Sample, and via that enters the external code, all computation happens there anyway.
  • I think leaving out uses S = Statistics; was simply a copy and paste error - I'm pretty sure I had it in my original version 3 code as well. But it doesn't matter - thanks for catching it!

Finally, I realized late last night that part of the evaluation in the external C code still happens inside Maple - the external code calls back into the interpreted Maple environment to evaluate part of the expressions. So I came up with the following:

ncubem4 := module()
local
    ModuleApply, normalSample, sampleStorage, processor;

    normalSample := Statistics:-Sample('Uniform(0, 1)');
    sampleStorage := Vector(0, 'datatype = float');

    processor := proc(s :: Vector(datatype = float[8]),
                      n :: integer,
                      N :: integer)
    local
        i :: integer,
        intermediate :: float[8],
        j :: integer,
        result :: float[8];

        result := 0.;
        for i to 2*n*N by 2*n do
            intermediate := 0.;
            for j from i to i + 2*n - 1 by 2 do
                intermediate := intermediate + (s[j] - s[j+1])^2;
            end do;
            result := result + sqrt(intermediate);
        end do;

        return result / N;
    end proc;

    processor := Compiler:-Compile(processor);

    ModuleApply := proc(n :: posint, N :: posint, m)
        if numelems(sampleStorage) < 2*n*N then
            # Sample storage is too small to contain sufficiently many
            # sampled numbers: allocate a new Vector.
            sampleStorage := normalSample(2*n*N);
        elif numelems(sampleStorage) > 2*n*N then
            # Sample storage is bigger than needed - fill only an
            # initial segment.
            normalSample(ArrayTools:-Alias(sampleStorage, [2*n*N]));
        else
            # Sample storage is exactly the right size - fill it up
            # with fresh samples.
            normalSample(sampleStorage);
        end if;

        return m * processor(sampleStorage, n, N);
    end proc;
end module;
CodeTools:-Usage(ncubem4(15, 1000, 1), iterations=1000):
memory used=1.20KiB, alloc change=0.75MiB, cpu time=2.30ms, real time=2.31ms

One caveat is that this will only work in hardware float mode; if you set Digits higher than evalhf(Digits) (which is currently 15 on all platforms) or you set UseHardwareFloats to false, the code will die miserably. However, this is stuff you shouldn't do with high precision anyway: the Monte Carlo scheme destroys the high precision anyway.

The numbers here are a bit skewed since this is my work machine instead of home - the numbers I get here for the other versions are: 200.64ms for ncubem, 31.55ms for ncubem2, 7.10ms for ncubem3. So here we have another factor of 3 over ncubem3.

Erik.

5 6 7 8 9 10 11 Last Page 7 of 20