Hi herclau,

I'm not sure what factor exactly you want to factor out - this one example does not provide enough information to define that uniquely. What would intuitively make sense to me would be to factor out the GCD of the entries of the vector, as follows:

myFactorization1 := proc(v :: rtable, $)
local divisor;
  divisor := foldr(gcd, op(convert(v, list)));
  return divisor, v / divisor;
end proc;

but for the case above, that would return b, <0, 1, 0, 5>. Maybe you would want the LCM of the nonzero entries of the vector? That could be done as follows:

myFactorization2 := proc(v :: rtable, $)
local divisor;
  divisor := foldr(lcm, op(remove(`=`, convert(v, list), 0)));
  return divisor, v / divisor;
end proc;

If you expect to find vectors from a zero-dimensional vector space then you need to adapt both routines, and if you expect all-zero vectors then you need to adapt the last one of these two for that.

There's no really nice solution for getting this as a true factorization that is a product of things, because multiplying a vector by a scalar always evaluates the product. You could instead have the procedures return `%*`(divisor, v / divisor). This is an expression that will turn into the original product once you apply 'value' to it:

a := `%*`(b, c); # this will remain unevaluated through all computations that don't apply value
value(1 + a*d); # this returns 1 + b*c*d

Hope this helps,

Erik Postma
Maplesoft.

Hi Neill,

You've already found that you can do rhs(Eq1) = lhs(Eq1). Just a tiny bit shorter is this:

algsubs((rhs = lhs)(Eq1), Eq2);

This works because if you use something like f = g as a function applied to x, then this evaluates to f(x) = g(x).

As you already suggest in the title, you can also use algsubs. Whenever the thing you are replacing is not just a variable name, algsubs is more reliable than subs. So I'd recommend using that for this case. An alternative would be to roll your own version of algsubs using solve or isolate.

Hope this helps,

Erik Postma
Maplesoft.

Hi afeddersen,

Sorry - I seem to have a habit of answering direct questions with "here is something different you could do that might accomplish your goals in a better manner", and I'm afraid this is going to be such an answer.

The label feature in Maple is great, but it mostly works with direct user interaction: there's no easy access to programmatically retrieve the value of label (x.y) for example, if x and y are integers; instead you have to enter these numbers manually. Loops, of course, are very much programmatic constructs, and thus these two don't mesh very well.

One option would be to just store whatever value is printed out in a ?table or an ?rtable (such as a ?Vector). So if you have a loop that looks like this:

running_total := 0;
for i to 25 do
  running_total := running_total + i:
  i^2 - running_total;
end do;

to print the numbers i*(i-1)/2, you could instead decide to store these in a Vector:

results := Vector(25);
running_total := 0;
 for i to 25 do
   running_total := running_total + i:
   results[i] := i^2 - running_total;
 end do;

Now you can access the ith value as results[i] - with the added advantage over labels that it also works programmatically, not just when you type it directly. Does that help for what you're trying to do?

Erik Postma
Maplesoft.

Hi sund0002,

Glad that the plot worked. For import of text files as matrices, you can use this:

Data := ImportMatrix("/home/epostma/downloads/200-1.1.txt", skiplines=2, datatype=float[8]);
Data := Data[..-2];

The skiplines=2 option is to get rid of the first two lines of the file. Because there's an empty line at the end of the file, you'll need to get rid of the last row of the matrix as well, which happens in the second command.

The datatype=float[8] option is not necessary, but it makes subsequent operations on the matrix faster. This is especially useful if you have larger data sets.

By the way, note that you now have a Matrix, not an Array. This should make the last step of our previous conversation unnecessary (but it won't hurt either): the conversion from Array to Vector.

Hope this helps,

Erik Postma
Maplesoft.

Hi Morten,

The easiest way of doing what you ask is to run map(assign, A, 0). However, whenever variables have a mathematical meaning (as they probably do if they occur in a matrix), there is often a better solution than just assigning a value to them. For example, if you want to know the value of an expression where the variables have that value, you can use a set of equations that provide values for the variables:
    values := map(`=`, indets(A, name), 0);
and then evaluate the expression at that set of values:
    result := eval((a + b + 2)^4, values);
This makes it much easier to also evaluate the expression at other values of your variables: if you had values assigned to a and b, then result := (a + b + 2)^4 just assigns the value 16 to result, with no way of retrieving the original expression in terms of a and b.

Hope this helps,

Erik Postma
Maplesoft.

The best solution for this sort of thing is generally ?plots/surfdata , except that requires that you already have data for a uniform grid of points, and your data is not in a uniform grid. So you'll need to fit your data to a surface first. And frankly, the arrangement of data that you show here is not going to be enough to get a very useful fit in between the data points that are given. But maybe you have more data that would just be too voluminous to post? Let's have a look at the options at any rate.

The usual suspects for constructing such a fit are to be found in the ?CurveFitting package, but here most procedures either again require data for a uniform grid of points, or they accept only one-dimensional input data whereas you have 2D input. Another option would be the ?Statistics/Regression subpackage. That provides two functions that could be useful here: ?Statistics/LinearFit and ?Statistics/NonlinearFit. Both can fit a nonlinear expression to data; the first command requires that that expression is linear in the parameters that are to be fitted but provides globally optimal results, the second command is more widely applicable but provides only locally optimal results. The ?Statistics/Fit command provides a way to let Maple choose which of the two commands is appropriate.

A final option, one that could be considered nonparametric, is to use a kernel density fit. Again, Maple has a facility for that (?Statistics/KernelDensity), but it only allows for 1D input / 1D output data. But this would maybe be useful to construct ourselves. Essentially what want to do then, is: given x- and y-coordinates, take a weighted average of the c-values corresponding to our input points, where the weight for c_i depends on the vector (x - a_i, y - b_i). We'll do that as follows:

a := <.19375, .19494, .19581, .19632, .19646, .19624, .19565, .19472, .19349, .19198, .19024,  .18833>;
b := <-0.1517e-1, -0.1515e-1, -0.1512e-1, -0.1509e-1, -0.1506e-1, -0.1503e-1, -0.15e-1, -0.1497e-1, -0.1495e-1, -0.1492e-1, -0.1491e-1, -0.149e-1>;
c := <0.112346e-10, 0.111963e-10, 0.111544e-10, 0.111101e-10, 0.110649e-10, 0.110201e-10, 0.10977e-10, 0.10937e-10, 0.109011e-10, 0.108705e-10, 0.108458e-10, 0.108277e-10>;

kdFit := proc(xdata, ydata, zdata)
local xradius, yradius;
  xradius := (max - min)(xdata);
  yradius := (max - min)(ydata);
  return proc(x, y)
  local d, values;
    values := zip((xx, yy) -> exp(- 10 * (xx^2 + yy^2)), (xdata -~ x) / xradius, (ydata -~ y) / yradius);
    return (values . zdata) / add(d, d in values);
  end proc;
end proc;
ourFit := kdFit(a, b, c);
p1 := plots[contourplot3d](ourFit, min(a) .. max(a), min(b) .. max(b), axes=boxed, style=surfacecontour);
p2 := plots[pointplot3d](a, b, c, color=black, symbolsize=25);
plots[display](p1, p2);

Instead of the 3D contour plot that I plotted here, you can also create a regular, 'flat', contour plot by using plots[contourplot] - in which case you're probably less interested in the points, so you just need p1.

Hope this helps,

Erik Postma
Maplesoft.

Edit: had a typo in the code - fixed now.

It's not quite what you're asking for since it's also runtime behaviour, but I'd say it is somewhat in the spirit of what you're asking since it might allow you to catch a mistake before it actually leads to a wrong result. There are some underused safety-checking features that you can turn on with the command kernelopts(assertlevel = 2) (or by starting Maple with the -A 2 command line option): type assertions on local variables, assignments, and procedure return values. This makes Maple check (at runtime) the types of certain values, if you provide the type that they need to have. For example:

p := proc();
local a :: posint;
  a := -3/2;
  return a;
end proc;

q := proc()
local a;
  a :: posint := -3/2;
  return a;
end proc;

r := proc() :: posint;
local a;
  a := -3/2;
  return a;
end proc;

Each of these will throw an error when called; which admittedly is not as good as a compile time check (or rather, in Maple that would be an evaluation-(or simplification-)of-the-procedure type check), but it is at least earlier than when three procedures down the road you use the returned value to index a list.

Of course, parameter type checking (including delimiting the argument sequence by $ if you don't expect additional arguments) falls into the same category.

Hope that this slightly off-topic response still helps,

Erik Postma
Maplesoft.

Here's what I get if I copy and paste the 2D output from Maple:

This is an image! Note I have extended typesetting turned off in Maple, so that's presumably why BesselJ shows up as BesselJ instead of as J_0. I'm not sure exactly how copy and paste works in general, but I think the "providing" program (Maple, in this case) can offer a number of different formats for the copied content, and apparently the editing widget selects the image format. It's quite a wide image though; about double the width of the expression itself. So less suitable for in running text () because of the abundant amount of whitespace.

And here is what I get if I type or paste into the popup window:

It doesn't hang for me, but it didn't show anything useful in the little preview thingy until I also typed something into the box (a space, which I subsequently removed again). (That's probably because I pasted using the middle mouse button; when I retried using Ctrl-V it did update immediately.)

Note that the first image doesn't seem to have any useful "alt" text; the second one does. I concur with Alec and Acer that the second image looks horrible; the first is quite a bit better.

By the way, this is on Firefox 3.5.9 in Ubuntu 9.10.

Erik Postma
Maplesoft.

As in many other questions (and posts?), the first thing I am inclined to ask is: why do you want this, and could you do without it?

Maybe efficiency is of no concern, or the number of lists is always known beforehand to be very small. But if that is not the case, then Maple's lists are inefficient for building up incrementally; moreover, there is really no way to get around that if you need to keep all the intermediate results c[1], c[2], ... . But I find it hard to imagine a scenario where that would be absolutely necessary. For example, if you need access to only some of these partial lists after you have run the loop, then it might be better to just construct those partial lists on the fly: wherever in your program you had c[i], use seq(op(a[j]), j = 1 .. i) locally there. This is still not very efficient, but if you need only a few of these lists, it might be good enough.

If you need to do something to many (or all) partial lists, then the best solution might be to just store your data differently. For a few releases now, we have had automatically growing Array data structures (see the very last paragraph of ?rtable_indexing), which are very convenient to use; if you store all the entries of a[1], a[2], ..., in one big Array and use a separate Array (with datatype = integer[4] for bonus points) to store the last index of the elements belonging conceptually to a[1], a[2], etc., you can get much better performance. So the lists a[1] = [1, 2, 3], a[2] = [4, 5], a[3] = [], a[4] = [9, 4] would be represented with A1 := Array([1, 2, 3, 4, 5, 9, 4]) and A2 := Array([3, 5, 5, 7], datatype = integer[4]). You would then have access to the equivalent of c[i] by taking A1[1 .. A2[i]], or even better, ArrayTools[Alias](A1, [1 .. A2[i]]) (which would not involve any copying at all), or just A1 and knowing that you need to look only up to index A2[i]. Of course, all of this would result in an Array rather than a list, so you'd need to adapt the algorithm that operates on c[i].

Hope this helps,

Erik Postma
Maplesoft.

The easiest way to do that is by using a different command than Histogram: the  ?Statistics/TallyInto  command, which Histogram uses internally anyway:

a := [1, 3, 4, 6, 7, 8, 9, 4, 5, 9, 8];
minmax := min(a) - 1/2 .. max(a) + 1/2;
frequencies := map(rhs, Statistics[TallyInto](a, minmax, bins = (rhs - lhs)(minmax)));

This is much more reliable than obtaining the values from the Histogram itself, because the exact plot data structure that is used could be different in the next Maple release. However, if you need to do it, it can be done as follows (no guarantees about the future, of course):

valuesFromHistogram := proc(plotStructure :: specfunc(anything, PLOT), $)
local polygons, quadrangles;
  polygons := indets(plotStructure, 'specfunc'('anything', 'POLYGONS'));
  if nops(polygons) <> 1 then
    error "expected 1 POLYGONS structure, found %1", nops(polygons);
  end if;
  quadrangles := select(type, [op(polygons[1])], 'list'(['numeric', 'numeric']));
  return map(max, map[3](map2, op, 2, quadrangles));
end proc;

valuesFromHistogram(b);

Hope this helps,

Erik Postma
Maplesoft.

Hi Nick,

One issue with what you posted; this may be just a typo in copying things to the site. You write step = 0.1, which is just an equation, not an assignment; so Maple evaluates that equation but doesn't do anything with it. That means that in the UpdatePosition procedure, Maple calculates expressions involving step as a symbolic variable (since it doesn't have a value assigned to it). That's relatively expensive.

If you replace the equation with the assignment step := 0.1; you already get a speedup; for me, the loop takes 0.00135 seconds (measured by running it 1000 times). But using the Compiler[Compile] command to compile a (slightly rewritten) version into C and run that, I got it down to 0.000039 seconds (measured running it 10^4 times). Here's what I did:

pos := Matrix(3, 241, datatype=float[8]); # Need datatype=float[8] for the compiler to be able to deal with it.
vel := Matrix(3, 241, datatype=float[8]);

step := 0.1;

pos(1, 1) := 2; pos(2, 1) := 0; pos(3, 1) := 5;  #arbitrary values
vel(1, 1) := 13; vel(2, 1) := 0; vel(3, 1) := 9;

UpdateAllPositions := proc (maxi :: posint, pos :: Matrix, vel :: Matrix, step :: float[8])
local i :: posint;
  for i from 2 to maxi do
    pos[1, i] := pos[1, i-1]+vel[1, i]*step; #x position
    pos[2, i] := pos[2, i-1]+vel[2, i]*step; #y position
    pos[3, i] := pos[3, i-1]+vel[3, i]*step  #z position
  end do:
end proc;

UpdateAllPositionsC := Compiler:-Compile(UpdateAllPositions):

ttt := time():
to 10000 do
  UpdateAllPositionsC(241, pos, vel, step)
end do:
(time() - ttt)/10000;

By the way, since the velocities for i > 0 are all zero, there is an even faster solution - just copy the first column to all other columns (probably quickest if you use ArrayTools). But I guess that would be considered cheating, right?

A third option would be to use LinearAlgebra to do these operations; anyone want to try that?

Hope this helps,

Erik Postma
Maplesoft.

Edit: I meant to add a fourth option; if what you are modelling is essentially a discretized differential equation (as suggested by the position and velocity names for the matrices), possibly with discrete switching, it might be worth investigating the use of ?dsolve/numeric. It has had a very sophisticated facility for dealing with such discrete events, explained in ?dsolve/numeric/Events, and provides a very natural, accurate, and efficient way to deal with such stuff.

I tried a few random (positive) values for the other parameters, and for a couple of them, there was a unique real solution for i, but not with these values:

y = 2
a = 9
b = 1/10
c = 1/10
d = 5
g = 3/4
h = 1/8

Here, there are still a number of solutions (Maple finds 51 of them), but they're all complex, and I guess that would not be valid in the context of neural networks. Here's a plot of the right hand side of the equation minus y:

And indeed, we can see that if y > g + h, then substituting anything positive for the exponentials is not going to lead to a solution. So this seems to indicate that the function is not invertible -- unless this violates some constraint on the values that the parameters should have?

Erik.

There are two degrees of freedom involved here, expressible through three quantities: the number of bars, the width of each bar, and their product: the width of the range of data that is covered. If you use the ?Histogram command from the Statistics package, then you can either set these explicitly or select a rule of thumb that Maple should use to select these parameters; the default is one of the rules. In more detail:

• The range option is set to deduce by default, which includes all given data points. You can select a smaller (or indeed larger) range as well.
• The bincount option selects the number of "bins" the data is collected into - in other words, the number of bars. If this option is given, then the next option is ignored:
• The binwidth option selects either the absolute width of a bar (specified as a positive number), or the bin width calculation method as one of the keywords Sturges, Scott, or FreedmanDiaconis. The default is Sturges.

So if you want to have a fixed bar width, just specify the binwidth option and don't specify the bincount option.

Hope this helps,

Erik Postma
Maplesoft.

If I do , then I get both labels:

Hope this helps,

Erik Postma
Maplesoft.

Hi Chewyraver,

I haven't tried to manipulate your equation for very long, but I think it's not unlikely that there is no closed form solution for i. So that means that you will have to make do with the RootOf expression.

As you can see on the ?RootOf help page, RootOf is a general notation for "the solution of this equation". It may sound a bit like cheating to return such a general notation from solve ("What is the solution to this equation?" "The solution to this equation is the solution to this equation."), but in many cases there just doesn't exist a combination of functions that have a name and have been studied in the mathematical literature that is an inverse to the expression you give solve. This may well be a case like that. And in fact, even if there is a function with that name, that only means that some of the properties of that function are known. Take for example the equation ; the solution for z in terms of y is given by Lambert's W function, as . From a purely algebraic point of view, you could say that this is the definition of Lambert's W function, and the solution would be nothing more than a tautology. In some sense, the same is true for the square root, which could be defined as (one branch of) the inverse function of the square. So knowing to call the solution LambertW(y) is only useful because there might be other algorithms that can take advantage of that fact - for example, we know the derivative of LambertW, and we can evaluate it numerically.

However, most Maple commands can handle these RootOf expressions fairly well; in some sense, they are just another specially named function. You can differentiate them, and generally you can evaluate them numerically as well - provided you substitute values for the other parameters, as would be true for the LambertW function as well. You write that you don't know how to handle them; the natural question is then: what do you want to do with them? If you let us know, we might be able to help you.

Hope this helps,

Erik Postma
Maplesoft.

P.S.: you know that to enter an exponential, you need to type exp(...) or use the button on the expression palette, right? If you just type "e ^ x", then you get the variable e raised to the x'th power, but e is just another variable, not the base of the natural logarithm.