It turns out that the issue is mathematically much simpler - it is simply use of the ExpectedValue command in a way that it's not meant to be used. For every call to ExpectedValue that you make, the argument is just a constant floating point number. The expected value of a floating point number is always that same floating point number - it's constant! Similarly, the covariance between two constants is always zero - neither of them will ever change.

What you meant to do is take the covariance of the random variables that generated those floating point numbers. I'm looking into how your computation works and will get back to you in a few minutes. (I started writing my previous answer just before you posted yours, and sadly we don't get notified of newer answers or comments when writing answers. But it is still useful information for the case that you described in your initial question, so I'll let the answer stand.)

Erik Postma
Maplesoft.

@serena88: As to your first question / method, ?LinearAlgebra:-DotProduct works only on pairs of Vectors, not on Matrices. It is meant to represent the inner product on real or complex vector spaces. If the product you're trying to compute involves Matrices, use . (the dot), as you found already, or ?LinearAlgebra:-Multiply (which does the same thing, I think).

As to your second question: actually, B . C is well defined: B is a 3x1 Matrix and C a 1x2 Matrix, so you can multiply them to obtain a 3x2 Matrix (since nxk and kxm matrices can be multiplied for any n, k, m to obtain nxm matrices). If you would use, say, C^+ (a 2x1 Matrix) instead of C, then Maple would complain.

Hope this helps,

Erik Postma
Maplesoft.

@serena88: As to your first question / method, ?LinearAlgebra:-DotProduct works only on pairs of Vectors, not on Matrices. It is meant to represent the inner product on real or complex vector spaces. If the product you're trying to compute involves Matrices, use . (the dot), as you found already, or ?LinearAlgebra:-Multiply (which does the same thing, I think).

As to your second question: actually, B . C is well defined: B is a 3x1 Matrix and C a 1x2 Matrix, so you can multiply them to obtain a 3x2 Matrix (since nxk and kxm matrices can be multiplied for any n, k, m to obtain nxm matrices). If you would use, say, C^+ (a 2x1 Matrix) instead of C, then Maple would complain.

Hope this helps,

Erik Postma
Maplesoft.

+1, nice solution. I must admit I didn't try to optimize for speed at all since the demo input [2, 1, 4] was so small - but that's not likely to be the problem that Markiyan wanted to solve :)

It would probably be possible to generate a Compiler:-Compile'able version with some effort, if the highest possible speed is required (iterating over the integers 0 .. 2^n - 1 and using the binary representation to give the signs of the coefficients - with a maximum of 63 or so for n, but that will take way too long anyway), but it's probably better to invest the effort in making a good branch-and-bound strategy. The problem is NP-complete, since Knapsack reduces to it, so solving it perfectly for truly large problem sizes is not going to work.

Erik Postma
Maplesoft.

+1, nice solution. I must admit I didn't try to optimize for speed at all since the demo input [2, 1, 4] was so small - but that's not likely to be the problem that Markiyan wanted to solve :)

It would probably be possible to generate a Compiler:-Compile'able version with some effort, if the highest possible speed is required (iterating over the integers 0 .. 2^n - 1 and using the binary representation to give the signs of the coefficients - with a maximum of 63 or so for n, but that will take way too long anyway), but it's probably better to invest the effort in making a good branch-and-bound strategy. The problem is NP-complete, since Knapsack reduces to it, so solving it perfectly for truly large problem sizes is not going to work.

Erik Postma
Maplesoft.

Indeed the condition m=0 is the problem here: Maple is correctly telling us that there is, in general, no solution to that combination of equations. If we substitute m = 0 in both equations, then we get the solution x(t) = y2(t) = 0 for all t. Otherwise, if we leave m symbolic as is, we can get a solution if we leave the initial condition x(0) = 0 out; it's kind of complicated so I won't include it here. Leaving out the condition y2(0) = 0 still leaves an inconsistent system.

Hope this helps,

Erik Postma
Maplesoft.

Indeed the condition m=0 is the problem here: Maple is correctly telling us that there is, in general, no solution to that combination of equations. If we substitute m = 0 in both equations, then we get the solution x(t) = y2(t) = 0 for all t. Otherwise, if we leave m symbolic as is, we can get a solution if we leave the initial condition x(0) = 0 out; it's kind of complicated so I won't include it here. Leaving out the condition y2(0) = 0 still leaves an inconsistent system.

Hope this helps,

Erik Postma
Maplesoft.

A mini addition to your nice answer: I realized as I saw your code that you can now get some of the speed up of in-kernel code for filling Vectors with a constant, without using ?ArrayTools/Fill, as follows:

(**) bl := Vector(10^7, datatype=float[8]):
(**) kernelopts(bytesused);                 
                                                        82026976

(**) bl[1..5*10^6] := 1:
(**) kernelopts(bytesused);
                                                        82041496

(**) bl[..] := 2:
(**) kernelopts(bytesused);
                                                        82055672

So if you assign a constant to a subrange of an rtable, then no additional allocations are done. (ArrayTools:-Fill is still faster, though, presumably because it is more specialized: you could use the assignment statement for more complicated subsets of the data than ArrayTools:-Fill.)

Hope this helps,

Erik Postma
Maplesoft.

A mini addition to your nice answer: I realized as I saw your code that you can now get some of the speed up of in-kernel code for filling Vectors with a constant, without using ?ArrayTools/Fill, as follows:

(**) bl := Vector(10^7, datatype=float[8]):
(**) kernelopts(bytesused);                 
                                                        82026976

(**) bl[1..5*10^6] := 1:
(**) kernelopts(bytesused);
                                                        82041496

(**) bl[..] := 2:
(**) kernelopts(bytesused);
                                                        82055672

So if you assign a constant to a subrange of an rtable, then no additional allocations are done. (ArrayTools:-Fill is still faster, though, presumably because it is more specialized: you could use the assignment statement for more complicated subsets of the data than ArrayTools:-Fill.)

Hope this helps,

Erik Postma
Maplesoft.

For what it's worth, I've entered this bug into our internal tracking system; that's not a promise that I will have time to do something about it soon, but I'll do my best.

Erik Postma
Maplesoft.

For what it's worth, I've entered this bug into our internal tracking system; that's not a promise that I will have time to do something about it soon, but I'll do my best.

Erik Postma
Maplesoft.

One more little tweak: starting from acer's last version, replace the line

data[y, ix] := data[y, ix] + 1;

by

data[y, ix] := data[y, ix] + iy^(-0.5);

This emphasizes the first couple of iterations over the later ones, and thus shows the 'trails' better as the expression converges towards the attractor. See below. Of course iy^(-0.5) can be replaced by any other function of iy; maybe you'd want to emphasize the first couple of iterations even more, or you'd want to emphasize a few iterations around iy=500, or...

Erik.

@acer Awesome :) I'd like to lighten up the darker areas a bit; so I inserted

map[evalhf,inplace](x -> x^0.8, img_h); 

immediately after the call to worker. Using 50000 iterations per column instead of 500 (using a measly 8 seconds) I got this:

@acer : how is this?

with(ImageTools):
BF8 := proc(ra, rb, xr, yr, n)
local data;
  data := ImageTools:-Create(yr, xr, 1);
  worker(data, HFloat(evalf(ra)), HFloat(evalf(rb)), xr, yr, n);
  return data;
end proc:
worker := proc(
  data :: Array(datatype=float[8], order=C_order),
  ra :: float[8],
  rb :: float[8],
  xr :: posint,
  yr :: posint,
  n :: posint)
local
  ix :: posint,
  iy :: posint,
  r :: float[8],
  rscale :: float[8],
  x :: float[8],
  y :: posint;

  rscale := (rb - ra) / (xr - 1);
  for ix to xr do
    r := ra + rscale * (ix - 1);
    x := HFloat(0.5);
    for iy to 3 do
      # Warmup.
      x := r * x * (1-x);
    end do;
    for iy to n do
      x := r * x * (1-x);
      y := trunc(x * yr) + 1;
      data[y, ix] := data[y, ix] + 1;
    end do;
  end do;
end proc:
worker := Compiler:-Compile(worker):
m := BF8(2.6, 4, 1000, 1000, 250):
map[inplace, evalhf](`/`, m, max(m)):
m3 := m ^~ (1/3):
ImageTools:-Write("abc.tiff", m3):

Leads to this image (converted to png). As you see it's upside down; easy to fix but I didn't bother. If your screen is too narrow (as mine is - I use my monitor in portrait mode), there's a scaled down version below.

Erik.

@Joe Riel : I think I like the idea of a 'let ... in ...' statement and I'd probably use it frequently. It would potentially lead to very high levels of nesting, I imagine, which would be an issue for aesthetical reasons - I would probably request to have your maplev-mode not indent the blocks inside these statements. I wonder whether the memory manager could be made to take advantage of the immutability of these guys. Maybe it would also be nice to have a way to close multiple let... in... scopes at once. (It's interesting to note the difference with the ?use statement, by the way, which creates only a temporary binding (not even a real name) for an expression - it doesn't evaluate the expression at the beginning of the statement, but does so every time it is referenced.)

On the whole I agree with Jacques' comments regarding option remember. I think we now have a slightly better approach with ?option/cache in that it can't grow arbitrarily out of control as easily in terms of memory usage overhead, but it can still be tricky in terms of memory access overhead (and other overhead, and considerations such as mutable data structures and side effects).

Erik.