acer

The solution seems simple, but you may not like to hear it. The repeated nesting of unapply, function application, and evalf@simplify is unnecessary and is causing your code to be thousands of times less efficient than it ought to be.

Rewrite the entire worksheet and completely change your methodology. Forbid yourself from creating operators or procedures, altogether. Forbid yourself from using `unapply`. Forbid yourself from calling `evalf` on a symbolic expression which you then subsequently pass to `simplify`.

Use expressions instead of operators, function calls, and unapply. Call simplifiy sparingly. Never call evalf on a symbolic expression which is not expected to produce a purely numeric nonsymbolic result. Never call `simplify` on a symbolic expression with floating-point coefficients within it (and hopefully never create such expressions).

If you are bent on using some procedure, then let it be just a one procedure which would take as its arguments value for parameters ymin, yman, amin, amax, eps, etc. This would wrap around almost the whole sheet. Possibly, just possibly, consider having another inner procedure which might take as its arguments numeric values for y, a, d, etc. But that would be for moving the whole computation entirely to pure numerics -- and I don't know whether you want that or not, so investigate whether you want the symbolic approach first, done efficiently.

The define_external call assigned to `init` is unusual. It passes the 'MAPLE' option (as if it were calling a custom wrapper which did all its own argument processing), but also passes all the typed parameters for the external function (as if it were utilizing wrapperless external calling). Usually a define_external call has either one of those or the other, but not both. I could be wrong, but my guess is that you don't want/need the 'MAPLE' option for `init`, while you do for `FindZerosOfFirstAndSecondDerivatives` which is more obviously a custom wrapper.

acer

correct...

Answered: acer 32333

November 07 2011

4 0

In your first example , you have not copied the Vector assigned to `a`, when you assign that also to `b`. All that you've accomplished is to make the value of `b` be the same object.

This behaviour, in the first example, should not seem mysterious. The most natural way to get `b` to be assigned a copy of the Vector assigned to `a` is to use the copy command (a top-level command, not part of any package, which knows how to copy various sorts of object).

Note that there is a flip-side to this copying business, which is related to uniqueness. A Vector is a so-called mutable data structure, which means that its entries can be changed in-place. This is related, in some respects, to the fact there can be more than one distinct length-2 Vector with the entries x and y. Changes to one such Vector does not change the other.

>  restart:

> V1:=Vector([x,y]):

> V2:=copy(V1):

> {V1,V2}; # This surprises some people, but it shouldn't.

                                 /[x]  [x]\ 
                                { [ ], [ ] }
                                 \[y]  [y]/ 

> evalb(V1=V2); # This surprises some people, but it shouldn't.

                                    false

> V1[2]:=17:

> V2[2];
                                      y

So, some people are surprised that mere assignment does not make distinct copies. And some people are surprised that separate instances are recognized as being distinct.

You can contrast all that you've seen here for rtables (of which a Vector is one kind, the other two being Matrix and Array), with the behaviour of lists and tables (the latter of which also form the data structure used for lowercase array, matrix, and vector). How do these differ, you may ask? Well, for lists mutability is faked, which is much discusses as the reason why doing list-entry assignment is highly inefficient and therefore quite dubious. When you assign into just an entry of a list then the maple kernel actually goes and replaces the whole structure with a new object behind the scenes. And tables have last-name-eval, which makes them behave differently especially when passed as arguments to a procedure.

Which brings us to your second example. The ArrayTools:-Copy behaviour is a little unusual. But the broader patterns of behaviour are motivated in part to efficiency needs and to argument passing for procedures. I'll try and explain:

Suppose you create a length-2 Vector V with entries x and y, where x and y are not yet assigned. And then you subsequently assign a value to x, let's say a value of 5.

> restart:

> V := Vector([x,y]):

> x := 5:

The Vector V still has entries of x and y. If you access V[1], then you get the fully evaluated result of 5. But the entry of V is still just x. When you cause V to be printed, then really you're just seeing how the `print/rtable` command sees V when in turn it has to access the entries. You can still see the stored entries of V using `lprint`.

> V[1];

                                  5

> V;

                                     [5]
                                     [ ]
                                     [y]

> lprint(V);
Vector[column](2, {1 = x, 2 = y}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

Now, what would happen if we passed V as a argument to a procedure. Would we want all entries of V to be evaluated fully (as happens for some other structures...), or not?

Suppose we have an enormous Vector V, and the procedure will only access a small number of the entries a small number of times. In this scenario, we would very much want the entries of V to not all be fully evaluated, as far as what is seen inside the procedure. We'd want behaviour like at the top-level, outside the procedure, where evaluation only occurs upon individual entry access.

This is, in fact, the default behaviour to how rtables (eg. Vectors) are passed.

> f:=proc(S)
>    lprint(S);
>    S[1];
> end proc:

> f(V);
Vector[column](2, {1 = x, 2 = y}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

                                      5

Now suppose that all entries of V will be accessed by the algorithm within the procedure. and that this will happen many times repeatedly. In this scenario, we'd very much prefer it if the entries of V have already been fully evaluated already, inside whatever is assigned to that formal parameter of the procedure. Because then each subsequent access would be fast and not require the full evaluation and resolution. There is a (rarely mentioned) command to handle this (somewhat uncommon) situation, called rtable_eval. Since these situations relate to efficiency, it's important to be able to resolve the values of the rtable, inplace, without necessitating extra memory allocation that is entailed with fully copying and replacement of the structure.

> restart:

> V := Vector([x,y]):

> y := 13:

> lprint(V);
Vector[column](2, {1 = x, 2 = y}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

> f := proc(S)
>    rtable_eval(S,inplace);
>    lprint(S);
> end proc:

> f(V);
Vector[column](2, {1 = x, 2 = 13}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

> lprint(V); # different!
Vector[column](2, {1 = x, 2 = 13}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

All the same fundamental issues were also present for the older (now deprecated) table-based array, matrix, and vector objects. But since these had last-name-eval they faced the problems slightly differently. You can still see the workarounds, though. For example, several `linalg` package commands would accept a matrix argument A and immediately make an assignment to the local B of the form B:=eval(A). Sometimes that was done so that an inplace algorithm could act on local copy B, but sometimes it was done because the entries of last-name-eval matrix A were not all fully evaluated like gets done for other arguments to the procedure.

Now, the ArrayTools:-Copy command is interesting. It is not a kernel builtin, but it is a compiled routine. It appears to have some unusual behaviour, and I'd label your example as an anomaly. I'd stick with the top-level `copy` command instead, if I were you and I wanted to easily create a distinct copy of a Vector/Matrix/Array. Or assign the entries of the new Vector by some other means, if desired.

> restart:
> a := Vector([x,y]):
> b := copy(a):
> b;

                                     [x]
                                     [ ]
                                     [y]

> y:=13:
> a;

                                    [x ]
                                    [  ]
                                    [13]

> b;

                                    [x ]
                                    [  ]
                                    [13]

> restart:
> a := Vector([x,y]):
> b := Vector([0,0]):
> b[1],b[2]:=a[1],a[2]:
> b;

                                     [x]
                                     [ ]
                                     [y]

> y:=13:
> a;

                                    [x ]
                                    [  ]
                                    [13]

> b;

                                    [x ]
                                    [  ]
                                    [13]

> restart:
> a := Vector([x,y]):
> b := Vector([0,0]):
> b[1..2]:=a[1..2]:
> b;

                                     [x]
                                     [ ]
                                     [y]

> y:=13:
> a;

                                    [x ]
                                    [  ]
                                    [13]

> b;

                                    [x ]
                                    [  ]
                                    [13]

acer

add...

Answered: acer 32333

November 04 2011

1 0

  add(points[2*k,2],k=1..45);

acer

possibility...

Answered: acer 32333

November 03 2011

2 2

I'm not sure that I understand why you wrap MyTargetMultArgsFunctionExt inside G.

You might try these two (distinct and separate) approaches:

1) Get rid of G altogether, and make MyTargetMultArgsFunctionExt itself the object, and call Minimize(MyTargetMultArgsFunctionExt), which I have not tried on a bivariate external-call example.

2) Use the so-called Matrix form. You'd change the definition of MyTargetMultArgsFunctionExt below to be the define_external, but the rest would stay as I have it.

restart:

MyTargetMultArgsFunctionExt:=proc(x,y) (x-2)^2+(y-3)^4; end proc:
#MyTargetMultArgsFunctionExt := define_external(...);

G:=proc(V::Vector) MyTargetMultArgsFunctionExt(V[1],V[2]); end proc:

dG := proc(xy::Vector,grad::Vector)
  grad[1] := fdiff( MyTargetMultArgsFunctionExt, [1], [xy[1],xy[2]] );
  grad[2] := fdiff( MyTargetMultArgsFunctionExt, [2], [xy[1],xy[2]] );
  NULL;
end proc:

Optimization:-NLPSolve(2,G,objectivegradient=dG);

               [                      -19  [2.00000000000000]]
               [1.70100225269204398 10   , [                ]]
               [                           [2.99997969157599]]

acer

sparse storage format for float[8] Matri...

Answered: acer 32333

October 30 2011

1 0

The internal format of a storage=sparse, datatype=float[8] Matrix is a triplet of float[8] dense storage Vectors all of equal length. The first Vector represents the row coordinate "i" of each stored entry, the second Vector represents the column coordinate "j" of each stored entry, and the third Vector represents the value of the entry stored. In other words, each stored entry is given by a tuple. This is sometimes called coordinate list (COO) format.

Some computational routines (NAG, etc) may need the data to be in so-called compressed row storage (or CRS) format instead. It's not too hard to convert from one format to the other in a C wrapper, although it's a cost that you'd only want to do once, up front, per matrix.

There are functions in the Maple (OpenMaple, I guess) external API which allow one to access the triplet of Vectors with the coordinate list (COO) storage. I don't recall that there already is any pre-built, direct way to access these Vectors, from within Maple itself. But a way could be constructed, with a customized externally-called C function utility. I'm being a little loose with terminology. The data is in three contiguous arrays in memory, and the external API rtable constructor could be used to create three Maple Vectors which point at the respective addresses of the start of those three arrays (and such Vectors could be tagged as foreign, to prevent inadventant collection of the data which is shared with the original sparse Maple Matrix.)

I don't know if there is a direct way to construct a sparse Matlab matrix from such Vectors (or from such a triple of pointers to contiguous memory arrays, in either format described above), using Matlab commands. But if there were, then aliases (or, maybe safer, copies) of those three Vectors might be passed via the Matlab link.

If that's the case, then a C wrapper could be written which, when used via Maple's external calling mechanism, could take a sparse, float[8] Matrix and return three Vectors which pointed at aliases or copies of the (row, column, and value) data portions.

In summary, yes, you can find pointers to the three data portions by which a sparse, float[8] Matrix is stored in Maple. And you can also turn those into three distinct Maple Vectors. But, how to convert those into a sparse Matlab matrix in a running session, is another question to which I do not have an answer.

acer

please try.....

Answered: acer 32333

October 30 2011

1 1

This should not normally occur, and there is not something extra and special that one needs to do to enable this functionality.

Note that Fedora 14 is not listed as an officially supported 64bit Linux platform for Maple 15. But perhaps we might have some luck with figuring it out.

I see a few possibilities about what may be the cause: the external shared library actually is not being found, or the `SmallFactors` symbol is not being resolved within it, or some other dependency is not being resolved when that `SmallFactors` symbol is dlopen'd, or there is some mismatch in dependent libraries expected by that shared lib and whatever is on your OS.

Let's try and find out what precisely is the error when the define_external is made. In the routine from which that error is issued, there is a try..catch in the source code. Any error from the define_external call is caught and rethrown in the form you quoted. So let's try and see what the original (caught) error message is.

Could you please issue this command in your Maple, and report back with whatever that produces,

define_external('SmallFactors','MAPLE','LIB'="libmodLA.so");

Can you successfully run any other examples from help-pages of any commands in the LinearAlgebra[Modular] package? Is there a `processor` (or `processor32`) executable in any bin.XXXX subdirectory of your Maple installation, and if so could you run it from a terminal shell and report whatever it says? If there's no such executable, then could you say what CPU your machine has, with the exact model number?

acer

a few more...

Answered: acer 32333

October 27 2011

0 0

The original list mentioned `maptype` and `typematch` which, while important, don't come up as often as several of the other listed commands. If those are in, then maybe `foldl`, `foldr`, and `curry` could be in too.

acer

parameters vs arguments...

Answered: acer 32333

October 25 2011

4 0

I feel that (in Maple, at least) the significant difference between parameters (inside a proc) and arguments (to that proc) lies more in how you expect to use them than it does in some special coding practice. There are always exceptions, of course, but my preference is to manage parameters as formal arguments.

There are a few minor sins which I think this helps avoid. (The severity of the sin varies with circumstance and context, of course.)

One of those is declaing globals in a proc, when it's not necessary.
[We've seen the worksheets here, that have procs each with 2 formal parameters and a dozen declared globals. And at worst the results are confusion as to which `x` is which, debugging nightmares, slow code, and certainly no evalhf'ability.]
Using a procedure without knowing why.
[We've seen this too on mapleprimes: it starts off simple, and then somehow quotes get added because they seem to help, and then before you know it you're creating a procedure which you only ever call once. At this point, the central threads of purpose may have been lost. NB. This case in point is not even close to the extreme, which is something like,
```
fexpr:=unapply(expr,x);
plot(fexpr(x),x=a..b);  # and fexpr is never again used
```
Having to do a full restart,or even just have Maple do more "work", just to make changes to an example. Or put another way: the less re-calculation required in order to re-run with different parameter values, the better. And re-creating a new instance of the whole procedure, just to use different parameter values, is certainly more re-calculation. (Yes, I realize that there may be other quite different memory management reasons why one has to do a full restart. But it still be be decent practice to not force the requirement...)

So I might try something like this:

restart:

U := proc(c, s) if s=1 then log(c); else (c^(1-s)-1)/(1-s); end if; end proc:

Parameters1 := [ s = 2 ]:
U(c, subs(Parameters1,s));

                              1    
                            - - + 1
                              c    

# And now, without having to re-create any procedure,

Parameters1 := [ s = 1 ]:
U(c, subs(Parameters1,s));

                             ln(c)

As a very minor quibble -- if I may offer more opinion on coding style -- I suggest trying to use different names for formal parameters than are expected to be used as supplied arguments. Fr example, to help keep thing straight in our heads,

restart:

U := proc(C, S) if S=1 then log(C); else (C^(1-S)-1)/(1-S); end if; end proc:

Parameters1 := [ s = 2 ]:
U(c, subs(Parameters1,s));
                              1    
                            - - + 1
                              c    
Parameters1 := [ s = 1 ]:
U(c, subs(Parameters1,s));
                             ln(c)

The formal parameters of procedure U are now named differently, and may help avoids some confusion, but the way that U gets called and the final form results are the same.

acer

testeq...

Answered: acer 32333

October 23 2011

0 0

See testeq

acer

print, or inert form...

Answered: acer 32333

October 23 2011

1 0

Oftentimes, one can simply `print` the procedure, and copy & paste to your favourite editor.

Or you could use `ToInert` and `FromInert` which can potentially be demanding. (Eg, figure out the need for subsop([5,1,1,2,1]=... in this.)

The first example below works ok. But the second example is a problem, because you'd need to know or figure out that `i` is the first declared local, otherwise it's hard to match only the right thing.

restart:

p := proc() local i,y; y:=(i=30); end proc;

                p := proc() local i, y; y := i = 30 end proc;

FromInert(subsindets(ToInert(eval(p)),
                     specfunc(Or(specfunc(anything,_Inert_LOCAL),
                                 identical(_Inert_INTPOS(30))),
                              _Inert_EQUATION),
                     z->subs(30=50,z)));

                   proc() local i, y; y := i = 50 end proc;

p := proc() local i,j,y,x; y:=(i=30); x:=(j=30); end proc;

      p := proc() local i, j, y, x; y := i = 30; x := j = 30; end proc;

FromInert(subsindets(ToInert(eval(p)),
                     specfunc(Or(specfunc(anything,_Inert_LOCAL),
                                 identical(_Inert_INTPOS(30))),
                              _Inert_EQUATION),
                     z->subs(30=50,z)));

        proc() local i, j, y, x; y := i = 50; x := j = 50; end proc;

As a last resort, if you are using the commandline interface and have quit the session without having savelib'd your procedure `p`, you could try the session history file.

$ "C:/Program Files (x86)/Maple 15/bin.win/cmaple.exe"
    |\^/|     Maple 15 (IBM INTEL NT)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2011
 \  MAPLE  /  All rights reserved. Maple is a trademark of
   Waterloo Maple Inc.
      |       Type ? for help.
> interface(historyfile);
                        "C:\Users\Rubrum/.maple_history"

I suspect that in Maple 12 the result of the above command might be like ".maple_history" and be relative (implicitly) to currentdir().

acer

fnormal...

Answered: acer 32333

October 21 2011

3 1

See fnormal

Eg,

> restart:

> z := 2e-15*x - 3e-15*x*y + 5*x^2*y:

> fnormal(z,10,1e-15);

                    -15         -15          2  
                2 10    x - 3 10    x y + 5 x  y

> fnormal(z,10,1e-14);

                                2  
                             5 x  y

> fnormal(z); # session defaults

                                2  
                             5 x  y

acer

initialvalues (or use a global solver...

Answered: acer 32333

October 20 2011

1 1

The first thing to note is that success should not necessarily be measured by virtue of *all* the computed optimizing parameter values being very close to your original generating function. That's because the effects of some terms may swamp out the effects of others, the latter of which are hence somewhat free to vary. A reasonable gauge of success is to compare plots, instead.

Statistics:-NonlinearFit uses a "local optimization" method, and so the result may be locally optimal but not globally optimal. In other words, there may be other local optima which are better than the local optimum that is found.

There are two basic ways deal with this. 1) get and install a global optimization package such as DirectSearch or GlobalOptimization, or 2) supply a judicious/lucky set of initial values for the parameters.

Sometimes you can use ExponentialFit or LogarithmicFit to help guide you in choosing initial values for some of the parameters. For example,

fitting101.mw

acer

curve fitting...

Answered: acer 32333

October 20 2011

2 4

You mentioned other "forms", like an ellipse.

restart:

L := [[.241, 0.2e-2], [.241, 0.3e-2], [.242, 0.], [.243, 0.5e-2],
      [.246, -0.3e-2], [.246, 0.6e-2], [.247, 0.6e-2], [.248, -0.4e-2],
      [.248, 0.6e-2], [.25, -0.5e-2], [.252, -0.6e-2], [.253, 0.6e-2],
      [.254, 0.6e-2], [.255, 0.6e-2], [.261, -0.9e-2], [.261, 0.5e-2],
      [.265, -0.1e-1], [.266, 0.4e-2], [.269, 0.3e-2], [.271, -0.11e-1],
      [.272, -0.11e-1], [.272, 0.2e-2], [.278, -0.11e-1], [.279, -0.11e-1],
      [.28, -0.11e-1], [.28, -0.2e-2], [.283, -0.1e-1], [.284, -0.9e-2],
      [.284, -0.6e-2]]:

P1:=plot(L, style=point, color=black):

K:=Statistics:-Fit(a*(x-x0)^2-b*(x-x0)*(y-y0)+c*(y-y0)^2-1,
                   Matrix(L), Vector(nops(L)),[x,y],
                   initialvalues=[a=0,b=0,c=0,x0=0.26,y0=0]);

                                                       2                             
          K := 3163.59634684427 (x - 0.262991041730829) 

               + 9227.98709575456 (x - 0.262991041730829) (y + 0.00246191128098791)
                                                                  
                                                           2                             
               + 19957.6913974776 (y + 0.00246191128098791)  - 1              


ellK:=plots:-implicitplot(K,x=0.24..0.3,y=-0.02..0.02):

plots:-display(P1,ellK);

One may need to experiment with the initial values (x0 and y0 being easiest to guesstimate), though a global solver might do well automatically, for other data sets.

There is more than just one way to fit an ellipse to a hull. Above, it fits the points themselves in usual least-squares distance kind of way. But for your original, smaller, data set you can get quite different results if instead you fit the midpoints of the line segments of the hull, or if you fit the original points and the midpoints, or if you fit by minimizing the sum of areas (sum of integrals integral) between line segments of the hull and the ellipse curve (there too: a choice between minimal distance to ellipse from each point, or distance taken orthogonally from each point on hull's line segment).

acer

circle...

Answered: acer 32333

October 18 2011

2 4

noisycircle.mw

It's not clear that using the cross (orthogonal x and y coordinates) of two univariate normal distributions for the noise is best. It may be ok for the parametrized circle as the curve, on account of symmetry, but I suspect it might surprise if used for some other curves.

Here is some code that is reasonably fast and memory careful. The 2D plot data-structure for 30000 pts can be generated here in less than 1/10th of a second on an Intel i7.

(It could likely be done even faster, by using the Compiler on a procedure which populated the final `XY` Matrix. But that looks like overkill, especially when one considers how long `plot` takes to produce the `PLOT` structure, or how long the GUI takes to insert and display the graphic.)

But throwing up the plot onto the Worksheet takes much longer, or at least it does on Maple 12 for Linux 64bit. In Maple 15.01 it is quite fast on Windows 64bit, to insert the plot into a pre-made Plot Component. (Coincidentally, I had recently been preparing a blog post on plot insertion and rendering speeds in the Worksheet. This is a nice example to add there.)

In Maple 12, this looks really good at symbolsize=4 in the Worksheet. Also, it seemed to render more quickly for symbol=solidcircle which is interesting. But in Maple 12 the symbolsize=4 plot exported with right-click content-menu to a blank plot. In Maple 15 that same symbolsize (for the same solidcircle) appears much larger(!) and not nearly as nice. And in Maple 15, smaller symbol sizes produce a blank plot for me(!). I intend to submit bug reports.

Here's a screenshot of what it looks like in Maple 12.02,

And here's a screenshot of what it looks like in Maple 15.01. Once again: both have symbolsize=4 and symbol=solidcircle, this M12 plot exported as empty with right-click context-menu gif/jpeg driver, and in M15 no smaller symbol shows as nonempty even in the Worksheet.

Personally, I think that M12 offered the more sophisticated/professional appearance here, and it certainly had the more flexible symbol size functionality. All we need now is the best of both worlds: proper non-empty export using all drivers, small symbol size in the Worksheet, and the same good speed for inline plots as are obtained for Plot Component updating.

Now, what are some alternatives for adding noise a point on a parametrized curve in 2D? For any point, one could take an angle theta from a uniform distribution and a radius from a normal (or other) distribution. Or, should the normal of the curve at that point make an effect on the direction of the noise component? What would one expect, for an ellispse, say?

acer

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

32333 Reputation

29 Badges

Social Networks and Content at Maplesoft.com

MaplePrimes Activity

These are answers submitted by acer

base64 of dotm...

revision...

correct...

add...

possibility...

sparse storage format for float[8] Matri...

please try.....

a few more...

parameters vs arguments...

testeq...

print, or inert form...

fnormal...

initialvalues (or use a global solver...

curve fitting...

circle...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again

From:
To:

Custom Message (optional):