As a very low-tech solution you could use the ?system (or ?ssystem) command to start a new Pari-GP session with your generated script as an argument - something like

system("gp p");

possibly adorned with paths to gp and the script. This won't work very well if you need to interact repeatedly with the same Pari-GP session, but it sounds like your current solution doesn't do that either -- correct?

Hope this helps,

Erik Postma
Maplesoft.

For your first question (how to see whether an equivalent Matrix is contained in a given list or set), you can use  ?ormap:

AA := <1,2|3,4>;
BB := <5,6|7,8>;
CC := <1,2|3,4>;
DD := <2,3|4,5>;
matrices := [AA, BB];
ormap(verify, matrices, CC, Matrix); # returns true
ormap(verify, matrices, DD, Matrix); # returns false
ormap(LinearAlgebra:-Equal, matrices, CC); # returns true
ormap(LinearAlgebra:-Equal, matrices, DD); # returns false

Hope this helps,

Erik Postma
Maplesoft. 

MyMatrix := (a, b, c) -> Matrix(3,3,(k,l)->MyFunction(a,b,c,k,l));

By the way, I think you wanted MyMatrix(x, y, z) to obtain that matrix - right? What you typed for MySecondFunction doesn't make too much sense - you specify a function of x, y, and z, but neither of those variables occur on the right hand side of the arrow.

Hope this helps,

Erik Postma
Maplesoft.

This is intentional behaviour; it is sometimes used inside the library to select between different versions of the same command, or to select certain options. For example, see ?evalf or ?map or ?Statistics/StandardError. map can even take multiple indices.

The way to access this information is by examining ?procname. That is the name by which the current procedure was invoked. In the example above, procname will be equal to a[b]. To check whether a specific index was given, you could do something like this:

a := proc($)
  if type(procname, 'indexed') and [op(procname)] = ['b'] then
    printf("invoked with index b\n");
  elif not type(procname, 'indexed') then
    printf("invoked without indices\n");
  else
    error "invalid input: expected %1's index to be 'b' or absent, but received %2", op(0, procname), [op(procname)];
  end if;
end proc;

Note that op(procname) is typically enclosed in brackets to make it a list, because you don't know a priori whether someone might try to invoke your procedure with multiple indices.

Hope this helps,

Erik Postma
Maplesoft.

Hi Manohar,

The default method of integration for the stiff case is the Rosenbrock algorithm. The help page for  ?dsolve/rosenbrock lists some details that might be helpful. As a first step, I would set infolevel[dsolve] pretty high, to 5 or even 10, in order to get an idea of where it is that dsolve is spending all its time. The very first clue is, when does it take so much time? Is it during the initial call to dsolve, or is it later when you ask it to evaluate the solution at some point? (I'm assuming you're not using the calling sequence that includes the 'range' option.)

If the time is taken during the initial call, the issue might be computing the Jacobian; with a bit of work and use of the 'jacobian' option, you might be able to get around that. If the time is taken during actual evaluation of the solution, you could try playing with the tolerance options ('abserr', 'relerr', and 'initstep') if you have some flexibility in that.

Finally, the Rosenbrock method is most suitable for the case where you're not looking for extremely high accuracy - if you are already setting the tolerance options in order to get very high accuracy, then you might be better off using method=lsode[x] for some value of x. Or is that what you meant where you were talking about the "islode solver"?

Hope this helps,

Erik Postma
Maplesoft.

Hi alex_01,

One thing I found somewhat surprising (from a mathematical / statistical point of view - I can see how it could seem natural from an experimentalist point of view) is that after taking the subsets of data that are sufficiently far away from the mean, you take just the correlation coefficient of these data sets - which means:

where X and Y are the subsets of data. I would think you might want to consider using the mean of the whole data set together, instead of Mean(X) and Mean(Y). Anyway - maybe it doesn't make much of a difference.

Either way, the condition on which the correlation is conditional, is that and . (If N goes to infinity, then the sample mean approaches the population mean arbitrarily well (for the normal distribution), so we'll ignore that difference.)

Now we can define the bivariate normal PDF (see e.g. Wikipedia), set mu to 0 and sigma to 1 (any different situation can be obtained from that by translation and scaling), restrict its domain to X and Y both greater than w, and start writing down integrals. I did this for a while and got the following horrendous expression for the covariance:

Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*(x-(Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*x/
(Int(Limit(Pi^(1/2)*exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2))-erf((
w-w*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), y = w .. infinity)), x = w .. 
infinity), y = w .. infinity)))*(y-(Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*x/(Int(Limit(Pi^(1/2)*
exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2))-erf((w-w*rho^2+rho^3*y
-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), y = w .. infinity)), x = w .. infinity), y = w .. 
infinity)))/(Int(Limit(Pi^(1/2)*exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/
(2-2*rho^2)^(1/2))-erf((w-w*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), 
y = w .. infinity)), x = w .. infinity), y = w .. infinity)

I'm not very optimistic that we are going to be able to figure out the behaviour of this expression as w grows for various values of rho. It's not unlikely that there are other, better, explanations, though.

Erik.

@vicaut : You can use Kamel's idea as a component to building a solution that includes all the roots. In fact, I would suggest using the ?RootFinding/NextZero command; it has somewhat better numerical behaviour and its behaviour is a bit better specified than fsolve for this type of question (where you want to "avoid" many roots). Define the following (Maple 15) procedure:

allRoots := proc(f :: procedure, rng :: range, $)
local result, root;
  result := Vector(0, 'datatype = float');
  root := RootFinding:-NextZero(f, lhs(rng), 'maxdistance' = rhs(rng) - lhs(rng));
  while root <> FAIL and root <= rhs(rng) do
    result(numelems(result) + 1) := root;
    root := RootFinding:-NextZero(f, root, 'maxdistance' = rhs(rng) - root);
  end do;
  return result;
end proc;

Now allRoots(rhs(ODE_4), 10 .. 20) should give you a Vector with all the (100 or so) roots that you were talking about.

Hope this helps,

Erik Postma
Maplesoft.

Here is an alternative approach that is fully in-memory. It won't give you a very beautiful plot, but the idea is there. My example (in Maple 15) uses the ODE with parameter a.

Actually, in production code, I'd use option cache rather than remember because it is less likely to blow up with memory use. However, here option remember will use less memory than the plot structures involved, so it will do for this example.

sys := [a*diff(f(x), x, x) + diff(f(x),x) + f(x) = 0, f(0)=0, D(f)(0) = 1];
solver := dsolve(sys, numeric, parameters=[a]):
solvedf := x -> eval(f(:-x), solver(x));

zerolist := proc(aa)
global solver, solvedf, parameters, a;
option remember;
local x0, result;
  solver(parameters=[a=aa]);
  x0 := RootFinding:-NextZero(solvedf, 0);
  result := Vector(0);
  while x0 < 20 do
    result(numelems(result) + 1) := x0;
    x0 := RootFinding:-NextZero(solvedf, x0);
  end do;
  return result;
end proc;

dispatcher := i -> proc(aa)
local result;
  result := zerolist(aa);
  if numelems(result) >= i then
    return result[i];
  else
    return NULL;
  end if;
end proc;

plot(dispatcher(1), 1/5 .. 5):
# This does adaptive plotting, populating the remember table of zerolist in an appropriate way.
rememberedXValues := sort([indices(op(4, eval(zerolist)), 'nolist')]);
maxlength := max(map(numelems, [entries(op(4, eval(zerolist)), 'nolist')]));

plot([seq(dispatcher(i), i = 1 .. maxlength)], 1/5 .. 5, 'sample' = rememberedXValues, 
     'adaptive = false', 'view = [DEFAULT, 0..15]');

This results in the following graph:

Hope this helps,

Erik Postma
Maplesoft.

Hi zbay,

The easiest solution to achieve what you ask for is to simply run the Statistics:-Mean and Statistics:-StandardDeviation procedures. So you load the Statistics package (by typing with(Statistics): ), enter your values into a  ?Vector V, then type Mean(V), StandardDeviation(V); to get the two values.

If you want to use the procedure you have written already, then I'd adapt it as follows:

myStandardDeviation := proc(S)
local x, mean;
  mean := Statistics:-Mean(S);
  # or if you prefer, mean := add(x, x in S) / numelems(S);
  return sqrt(add((x - mean)^2, x in S) / (numelems(S) - 1));
end proc;

Note the following things:

1. In cases such as this, where you're adding a concrete set of numbers instead of constructing a symbolically described sum, ?add is more appropriate than  ?sum .
2. You forgot to specify which the running variable is in the sum.
3. I use Maple 15's command  ?numelems , which makes this automatically work for lists, Vectors, sets, and essentially all data structures you can think of.

Hope this helps,

Erik Postma
Maplesoft.
  

Hi mutmurat,

It's highly unlikely that the number of variables is going to be the limiting factor in your computation. You can easily, say, differentiate simple expressions in thousands of variables. However, if the computation requires a lot of computations of large, complicated expressions, which are subsequently used in building ever larger expressions, there are circumstances where this will quickly grind to a halt. (A notorious, innocent-seeming example is inverting a symbolic matrix. Some of these matrices look quite innocent but symbolic inverses cannot fit into memory of a current machine.)

I'm aware that this does not really answer your underlying question. Your best bet is probably simply trying your computation. Be sure to ask here again if you get stuck, and provide lots of details about what you're doing!

Hope this helps,

Erik Postma
Maplesoft. 

Try the RealDomain package.

RealDomain:-solve(x^4=1); # returns -1, 1

It's not perfect, but it's pretty good.

Hope this helps,

Erik Postma
Maplesoft.

One way that you can sometimes attack these problems is by reducing the amount of normalization that Maple does. For matrix inversion, for example, every potential pivot has to be checked for being equal to 0, which is by default done by converting to a factored normal form. The method of normalization used in LinearAlgebra is configurable, however, as is the method of zero testing, and that can sometimes lead to a solution in a case like this. However, the result, since it is not normalized as much as you would like, is generally not very useful.

In this case, we define C as above and then do:

Normalizer := x -> x; # No normalization at all.
Testzero := testeq; # Fast heuristic zero testing - still extremely reliable.
inv := LinearAlgebra:-MatrixInverse(C):

The computation is now actually pretty fast. But note that I put a colon after that last statement, not a semicolon. That's very much intentional, since Maple can't possibly print any of the entries of the matrix - they are much too big!  For example,

length(inv[1, 1]);

returns 1611432722. That is, writing out the [1, 1] entry would require around 10^10 words of memory! The reason Maple can deal with such expressions at all is that it uses a DAG structure for representing expressions, in which repeating subexpressions use memory only once.

So can you use inv at all? Yes, you can. You can, for example, evaluate inv at a point in (x1, ..., y4)-space by substituting numbers in. However, it's not certain that this is actually faster than substituting the values into C and then inverting the matrix numerically. I only did a preliminary test, as follows:

variables := convert(indets(C), list):
mypoint := map(v -> v = RandomTools:-Generate(float(range=0..1, method=uniform)), variables);
CodeTools:-Usage(eval(inv, mypoint), iterations=1000):
# memory used=0.54MiB, alloc change=0 bytes, cpu time=3.88ms, real time=3.89ms
CodeTools:-Usage(eval(C, mypoint)^(-1), iterations=1000):
# memory used=63.66KiB, alloc change=0 bytes, cpu time=890.00us, real time=892.00us

This test is unfair, since the matrix inversion happens in a BLAS (in optimized external code) whereas evaluating inv happens in interpreted Maple code. So maybe the evaluation of inv could be improved by setting things up so that you can use evalhf or even Compiler:-Compile, but that would take some work.

All in all, I think in this case you're better off substituting the values first.

Hope this helps,

Erik Postma
Maplesoft.

Hi Andreas,

I entered the data in your example; this looks like data that would be good to fit with an exponential function. (The Excel notation is a little puzzling, but it seems that that's what happens there too.) An exponential function is of course nonlinear, but the fitting process itself is linear: you take the logarithms of all the observations, do a linear fit, and then apply exp to the result. Maple supports this using the ?Statistics/ExponentialFit command:

with(Statistics):
data := <0.07, 2,0.16, 0.24, 0.36, 0.5, 0.8, 1.2, 1.8, 8, 4, 6, 9, 15, 20, 30, 43, 66.9, 100, 150>;
PointPlot(data, axis[2]=[mode=log]); # to show that an exponential fit will work
result := ExponentialFit(<seq(1 .. 20)>, data, x);

This assigns result := 0.0866208931860853*exp(0.365097831283158*x). In order to find the R^2 value for this (linear!) fit, we need the residual sum of squares (easy to obtain with ExponentialFit) and the so-called total corrected sum of squares. For this total corrected sum of squares, we can't use ExponentialFit and we'll need to do the transformation by hand. We can obtain the total sum of squares by seeing that the second central moment is the average of the squared deviations from the mean, whereas the total sum of squares is the sum of the squared deviations from the mean - thus they differ by the factor given by the number of data points. Putting this together, we obtain:

ssRes := Statistics:-ExponentialFit(<seq(1..20)>, data, x, 'output = residualsumofsquares');
logdata := log~(data);
ssTotal := numelems(logdata) * CentralMoment(logdata, 2);
Rsquared := 1 - ssRes / ssTotal;

This results in a value of 0.9167199079, which is what Excel also seems to get.

If you have a more esoteric fitting function and you need to use ?Statistics/NonlinearFit or, better, ?Statistics/Fit, you can use the same approach to find the pseudo-R^2, except you would take the central moment of data - since there is not necessarily an equivalent of logdata that makes sense in that case.

Hope this helps,

Erik Postma
Maplesoft.

Hi hthomas10,

Do you mean for a matrix θ? If so, you'll want to use ?LinearAlgebra/NullSpace.

Hope this helps,

Erik Postma
Maplesoft.

Hi @jimmyinhmb,

A very interesting question indeed, and leading to a fascinating discussion. I noticed that one of your questions has not been answered at all, and one only implicitly. First your question 1. To do an in-place sort of an rtable, you can call sort['inplace']. Example:

restart():
v := Statistics:-Sample(Normal(0, 1), 10^7):
CodeTools:-Usage(sort['inplace'](v));

tells me:

memory used=0.52KiB, alloc change=0 bytes, cpu time=2.05s, real time=2.05s

Question 3 about ArrayTools:-Alias was answered implicitly already, by acer's use of it: that warning is there because if you really do your best, you can crash Maple's kernel with it. But in order to do that, you need to:

1. manually clear the attributes of the returned rtable, so that it doesn't reference the original rtable anymore;
2. make sure there is no reference to the original rtable anywhere else;
3. cause a garbage collection to occur.

If you're not actually trying to do that, then use of ArrayTools:-Alias is perfectly safe.

Hope this helps,

Erik Postma
Maplesoft.