acer

32385 Reputation

29 Badges

19 years, 343 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You already knew about assignment statements, and I suspect that you also knew about ways to concatenate names. You've put these actions together.

Note that it is only practical in a modest number, and since the generated names are always global it's far less useful as a programming technique.

Also, each global new name will get checked against the assigned names stored in all .mla and .m files in libname (once at least, up to first evaluation). That seems to scale roughly linearly, IIRC. I once saw an example with a 15sec wait in the computation where the kernel was doing nothing but i/o overhead.

So it's sorta, kinda ok as an off-the-cuff thing at the top-level, but in my opinion doesn't belong in a group of solid programming techniques for serious production line work.

You might consider this (although nothing is fool-proof),

  select(u->is(u,real), evalc~(allvalues~(coef7))) assuming alpha[1,2]::real, alpha[2,6]::real;

or, since it is only the two assumed names as dependent names,

  select(u->is(u,real), evalc~(allvalues~(coef7))) assuming real;

I declared z as local outside the do-loop as, alas, the 2D Input parser doesn't accept your original.

with(LinearAlgebra); T := proc (n::integer) local t0, tn, i, t, t1, z; if 0 < n then t0 := Matrix(1, 1, 1); tn := Matrix(1, 1, 0); for i to n do z := 2^i; t := Matrix([[Add(t0, tn), `~`[`.`](x, t0)], [t0, ZeroMatrix((1/2)*z)]]); t1 := Matrix([[`~`[`.`](y, t0), ZeroMatrix((1/2)*z)], [ZeroMatrix((1/2)*z), ZeroMatrix((1/2)*z)]]); t0 := t; tn := t1 end do else ('T')(n) end if; t0 end proc

proc (n::integer) local t0, tn, i, t, t1, z; if 0 < n then t0 := Matrix(1, 1, 1); tn := Matrix(1, 1, 0); for i to n do z := 2^i; t := Matrix([[Add(t0, tn), `~`[`.`](x, t0)], [t0, LinearAlgebra:-ZeroMatrix((1/2)*z)]]); t1 := Matrix([[`~`[`.`](y, t0), LinearAlgebra:-ZeroMatrix((1/2)*z)], [LinearAlgebra:-ZeroMatrix((1/2)*z), LinearAlgebra:-ZeroMatrix((1/2)*z)]]); t0 := t; tn := t1 end do else ('T')(n) end if; t0 end proc

T(2)

T(2)

NULL

Download G2_ac.mw

You might consider editing your procedures in 1D plaintext Maple Notation (in Execution Groups, or Code Edit Regions).

ps. I am curious, is your task going to go anywhere near the LinearAlgebra[Generic] subpackage? (I know nothing of "Domineering".)

You haven't stated why LinearAlgebra:-LinearSolve is not adequate for your computations. You haven't mentioned whether you are doing purely floating-point computations, and if so whether those are purely real.

I would not be surprised if you needed some large examples (or a great many of them) before the LinearSolve overhead became a bigger bottleneck than anything else in your code -- as a tentative judgement based from general experience with member's code.

But, below is an example that uses the floating-point real tridiagonal linear solver from LAPACK.

This uses two functions, one to do the LU decomposition and another to do the back/forward substitutions for one or more RHSs.

Hopefully I got the integer array's widths correct -- I might have forgetten some detail of the linkage model in use.

As I mention below, this has a little extra overhead (probably avoidable, using other access means).

This was run on Linux. To run on MS-Windows you'd have to changed the name of the external shared library, from "libmkl_rt.so" to the appropriate .dll. And even then it might run amok of MS-Windows awkwardness.

restart;

 

This works on my 64bit Linux, using either Maple 18 or Maple 2019.2.

I don't know whether it works on 64bit MS-Windows.

Mostly I was interested in the possibility of accessing LAPACK's tridiagonal
solver, and wanted to check that it worked in principle. So this may or may

not work on your platform. (If it doesn't, then that's just because Intel's MKL
has a complicated and unhelpful dynamic link model.)

 

The following uses so-called "wrapperless" external calling, which incurs
some overhead that could be avoided with a dedicated external wrapper (bridge).
Also, a dedicated wrapper could extract the diagonals from a packed band

storage Matrix more efficiently.

 

DGTTRF:=module()
  local ModuleApply, dgttrf_external, longsz;
  longsz:=kernelopts(wordsize)/8;
  dgttrf_external := define_external('dgttrf_',

   '_N'::REF(integer[longsz]),

   '_DL'::ARRAY(1.._N-1,'order'='Fortran_order','datatype'=float[8]),

   '_D'::ARRAY(1.._N,'order'='Fortran_order','datatype'=float[8]),

   '_DU'::ARRAY(1.._N-1,'order'='Fortran_order','datatype'=float[8]),

   '_DU2'::ARRAY(1.._N-2,'order'='Fortran_order','datatype'=float[8]),

   '_IPIV'::ARRAY(1.._N,'datatype'=integer[4]),

   '_INFO'::REF(integer[longsz]),

   ':-LIB'="libmkl_rt.so");

  ModuleApply:=proc(AA::Matrix(datatype=float[8]),{preserve::truefalse:=true})

   local d,dl,du,du2,ipiv,n,A,extinfo,ifail;

   n:=op([1,2],AA);
   if preserve then A:=copy(AA); else A:=AA; end if;

   dl:=Vector(n-1,(i)->A[i+1,i],'datatype'=float[8]);
   d:=Vector(n,(i)->A[i,i],'datatype'=float[8]);
   du:=Vector(n-1,(i)->A[i,i+1],'datatype'=float[8]);
   du2:=Vector(n-2,'datatype'=float[8]);
   ipiv:=Vector(max(1,n),'datatype'=integer[4]);

   extinfo:=0;

   dgttrf_external('n',dl,d,du,du2,ipiv,'extinfo');
   if extinfo=0 then NULL;

   elif extinfo<0 then error "invalid %-n argument",-extinfo;

   elif extinfo>0 then error "error code %1", extinfo; end if;
   return dl,d,du,du2,ipiv;

  end proc:
end module:

DGTTRS:=module()
  local ModuleApply, dgttrs_external, longsz;
  longsz:=kernelopts(wordsize)/8;
  dgttrs_external := define_external('dgttrs_',
   '_ITRANS'::string,

   '_N'::REF(integer[longsz]),
   '_NRHS'::REF(integer[longsz]),

   '_DL'::ARRAY(1.._N-1,'order'='Fortran_order','datatype'='float[8]'),

   '_D'::ARRAY(1.._N,'order'='Fortran_order','datatype'='float[8]'),

   '_DU'::ARRAY(1.._N-1,'order'='Fortran_order','datatype'='float[8]'),

   '_DU2'::ARRAY(1.._N-2,'order'='Fortran_order','datatype'='float[8]'),

   '_IPIV'::ARRAY(1.._N,'datatype'='integer[4]'),
   '_B'::ARRAY(1.._LDB,1.._NRHS,'order'='Fortran_order','datatype'='float[8]'),

   '_LDB'::REF(integer[longsz]),
   '_INFO'::REF('integer[4]'),

   ':-LIB'="libmkl_rt.so");

  ModuleApply := proc(dl,d,du,du2,ipiv,b)

   local itrans,extinfo,ifail,ldb,n,nrhs;
   (itrans,extinfo):="N",0;
   if b::Matrix then
     (n,ldb,nrhs):=op(1,d),op(1,b);
   else
     (n,ldb,nrhs):=op(1,d),op(1,b),1;
   end if;

   dgttrs_external(itrans,'n','nrhs',dl,d,du,du2,ipiv,b,ldb,'extinfo');
   if extinfo=0 then NULL;
   elif extinfo<0 then error "invalid %-n argument",-extinfo;

   elif extinfo>0 then error "error code %1", extinfo; end if;

  end proc:
end module:

 

 

A short example.

 

AT:=LinearAlgebra:-RandomMatrix(3,datatype=float[8],generator=0.0..100.0,shape=band[1,1]);

AT := Matrix(3, 3, {(1, 1) = 81.4723686393179, (1, 2) = 90.5791937075619, (1, 3) = 0., (2, 1) = 12.6986816293506, (2, 2) = 91.3375856139019, (2, 3) = 63.2359246225409, (3, 1) = 0., (3, 2) = 9.75404049994095, (3, 3) = 27.8498218867048})

V:=LinearAlgebra:-RandomMatrix(3,1,generator=-100.0..100.0,datatype=float[8]);

V := Matrix(3, 1, {(1, 1) = 92.9777070398553, (2, 1) = 91.5013670868595, (3, 1) = 9.37630384099677})

(dl,d,du,du2,ipiv) := DGTTRF(AT):
X:=copy(V):
DGTTRS(dl,d,du,du2,ipiv,X);
X;

Matrix(3, 1, {(1, 1) = 0.163654386495569e-1, (2, 1) = 1.01175967943733, (3, 1) = -0.176820178759307e-1})

LinearAlgebra:-Norm(AT.X-V);

0.142108547152020037e-13

 

 

A larger example, with many RHSs.

 

N:=2^12;
maxiter:=2^10;
AT:=LinearAlgebra:-RandomMatrix(N,datatype=float[8],generator=0.0..100.0,shape=band[1,1]):
V:=LinearAlgebra:-RandomMatrix(N,maxiter,generator=-100.0..100.0,datatype=float[8]):

4096

1024

XT:=Matrix(N,maxiter,datatype=float[8]):

 

ArrayTools:-Fill(N,0.0,XT):
st:= time[real]():
(dl,d,du,du2,ipiv) := DGTTRF(AT):
for i from 1 to maxiter do
  X:=V[..,i]:
  DGTTRS(dl,d,du,du2,ipiv,X);
  XT[..,i]:=X:
end do:
time[real]()-st;
LinearAlgebra:-Norm(AT.XT-V);

.210

0.974918948348779679e-8

 

Download DGTTRF.mw

There are a few things that can be said.

1) No, it is not necessary to use active int for this. In fact doing so can be unnecessarily very expensive, since it can cause the integral to be attempted symbolically for each new call to the objective procedure (where each might fail, after a while). See first example in the attached.

There can sometimes be a mixup for Mimimize when called with procedure (rather than expression) form for the objective. It happens more rarely then it did long ago. It's caused mostly by some muddle about generating a procedure for evaluating the derivative numerically, via automatic differentiation. My first example does a few things to avoid that.

2) Your second bullet point sounds like it might be an evaluation mix-up. See the last example in the attachment with proc4, which may be what you were looking for here. The point is to keep straight that the K,r,x in globally defined df are themselves global (whereas the K,r in proc4 are procedure parameters). It helps to always test the objective procedure at a single candidate point, as sanity-check. Sorry if I misunderstand you here.

3) There's a battle between the numerical optimizer's wish to establish an optimum that satisfies its optimalitytolerance control option's value versus the numeric integrator's specification to only provide a result that needs to be convergent up to its own accuracy tolerance option epsilon. I know that sounds like a mouthful, but analogous struggles can happen between evalf/Int, fsolve, and Minimize. The "outer" one needs the "inner" one to be accurate enough and stable enough that it can establish its own criteria for convergence. So, in your later example I have forced values for those options. Sometimes it requires fiddling, also because evalf/Int might also require high enough working precision to attain its target accuracy.

Unfortunately, I lost patience fiddling and threw in a scaling factor, to get the integrand to a nice range. It works, but muddies the question of what is really needed.

I have used unapply within the call to evalf(Int(...)) since I've found that all too often it is a useful way to prevent the numeric integrator's zealous approach to checking for discontinuities when presented an integrand in expression form. I didn't check whether it was completly necessary here.

4) I didn't accomodate the case that the numeric integration ever fails to produce a numeric value. The proc3 could (and probably should) test that the numeric integration has returned a numeric value. If it hasn't then it your luck depends on which minimizer you're using, and how it handles ever getting a nonnumeric value for the objective. Some won't handle the discontinuity that a default "large" fallback value can cause. Some won't handle a nonnumeric Float(undefined). DirectSearch handles occasional nonnumeric values from the objective better. 

5) Keep in mind that for the multivariate case Optimization is a local rather than a global optimization package. (The add-on package DirectSearch provides global optimizers.)

The attached worksheet runs pretty quickly. But I had to fiddle to get that.

Practice problem using Optimization:-Minimization

restart

with(Optimization)

proc1 := proc (K, r) local f1, f2, x; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; f1 := exp(.3*x); f2 := K/((K-1)*exp(-r*x)+1.0); evalf(Int((f1-f2)^2, x = 0 .. 3.2)) end proc

Minimize(proc1, 10 .. 100, 0 .. 1)

[0.243305733083052001e-4, Vector[column](%id = 18446883878292498062)]


Actual problem:

f1 := exp(-2.52428056431082+.468950959387562*x-0.403934021717953e-2*x^2); f2 := K*P0/(exp(-r*x)*K-exp(-r*x)*P0+P0); P0 := eval(f1, x = 0)

0.8011593034e-1

An approximation by hand:

plot([f1, eval(f2, [K = 0.25e5, r = .34])], x = 0 .. 38, 0 .. 0.15e5, legend = ["f1", "f2"])

``

proc3 := proc (K, r) local f1, f2, x; global P0; Digits := 15; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; P0 := 0.8011593034e-1; f1 := exp(-2.52428056431082+.468950959387562*x+(-1)*0.403934021717953e-2*x^2); f2 := K*P0/(exp(-r*x)*K-exp(-r*x)*P0+P0); 0.1e-6*evalf(Int(unapply((f1-f2)^2, x), 0 .. 38, epsilon = 0.1e-5, method = _Dexp)) end proc

sol := Minimize(proc3(K, r), K = 0.10e5 .. 0.60e5, r = .2 .. .4, optimalitytolerance = 0.1e-8)

[.184857934786779005, [K = HFloat(17139.92551737672), r = HFloat(0.34703798628407145)]]

plots:-display(
  plot3d(proc3(K, r), K = 1000. .. 60000., r = 0.2 .. 0.4),
  plots:-pointplot3d([eval([K,r,sol[1]],sol[2])],
                     symbolsize=20, symbol=solidsphere, color=red),
  view=0.0 .. 1e1,
  size=[400,400]);

df := (f1-f2)^2

(exp(-2.52428056431082+.468950959387562*x-0.403934021717953e-2*x^2)-0.8011593034e-1*K/(exp(-r*x)*K-0.8011593034e-1*exp(-r*x)+0.8011593034e-1))^2

proc4 := proc (K, r) Digits := 15; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; 0.1e-6*evalf(Int(unapply(eval(df, [:-K = K, :-r = r]), x), 0 .. 38, epsilon = 0.1e-5, method = _Dexp)) end proc

 

proc3(1010, .35), proc4(1010, .35)

35.0925066706428, 35.0925066706428

Minimize(proc4(K, r), K = 0.10e5 .. 0.60e5, r = .2 .. .4, optimalitytolerance = 0.1e-8)

[.184857934786779005, [K = HFloat(17139.92551737672), r = HFloat(0.34703798628407145)]]

 

Download MaplePrimes_Minimization_ac.mw

restart;

expr := y[1](t)*y[2](t)*(y[1](t)+y[2](t))^3:

frontend(diff, [expr, y[1](t)]);

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

Physics:-diff(expr, y[1](t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

subs(y1t=y[1](t),diff(subs(y[1](t)=y1t,expr),y1t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

thaw(diff(subs(y[1](t)=freeze(y[1](t)),expr),freeze(y[1](t))));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

By "group" do you mean sort, or separate?

If you mean sort then you could use the sort command. If you mean separate then you could use the select (or selectremove command).

In any case, it would help if you explained how you wish to compare (for sorting) or discriminate (for selecting) between the eigenvalues.

Are they all purely real?

Here is an example of sorting a list of lists by the real value in each inner list's second entry.

L := [ ["SA",3], ["HF",2], ["JU",7], ["AK",1] ]:

sort( L, (a,b) -> a[2] < b[2] );

    [["AK", 1], ["HF", 2], ["SA", 3], ["JU", 7]]

sort( L, (a,b) -> a[2] > b[2] );

    [["JU", 7], ["SA", 3], ["HF", 2], ["AK", 1]]

If the second entries were complex-valued then you could also sort them by modulus or argument (or some combination, using and, etc).

The sort command has some powerful functionality, and I've only given a simple example. There are also matters like efficiency (eg. see the key option), ease-of-use, etc, which may or may not be important to your case.

sort( L, key = (a -> a[2]) );

    [["AK", 1], ["HF", 2], ["SA", 3], ["JU", 7]]

You might consider looking at the Help page for the sort command.

Here is a very simple example of separation according to a test on the second entry of each inner list. Again, you could make the selection procedure more complicated.

(L1,L2) := selectremove(a->type(a[2],odd), L):

L1;

      [["SA", 3], ["JU", 7], ["AK", 1]]

L2;

                [["HF", 2]]

If by "group" you meant something more complicated then you should tell us what you mean in detail. It might be that you're after functionality attainable using a command from the ListTools package (eg, its Classify, or Group commands, etc).

Check your syntax. It looks like some of your bracketed terms might accidentally be entered as part of function calls rather than as parts of any multiplication.

For example, it might be that you have the function call,
    UA_C(T_c-T_C)
while what you intended was the multiplication,
    UA_C*(T_c-T_C)

Otherwise, you won't be able to eliminate T_c, because it's part of an abstract function call.

This syntax issue might be a problem in more than one place (eg, also on UA_H(...) at least...).

We can't tell for sure because you have only provided us with an image of your code rather than actual code. You can upload and attach a Document using the green up-arrow in the Mapleprimes editor's menubar.

Your condition that all the variables are positive integers violates the constraint. It is violated even if the values are all as small as possible, ie. they have a value of 1.

restart

with(Optimization):

Maximize(5*x1+2*x2+7*x3+6*x4+x5+6*x6+8*x7+6*x8, {x1 >= 1, x2 >= 1, x3 >= 1, x4 >= 1, x5 >= 1, x6 >= 1, x7 >= 1, x8 >= 1, 10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8 <= 22}, assume = nonnegint)

Error, (in Optimization:-LPSolve) no feasible point found for LP subproblem

eval(10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8, [x1 = 1, x2 = 1, x3 = 1, x4 = 1, x5 = 1, x6 = 1, x7 = 1, x8 = 1]);

61

``

Positive_integer_ac.mw

Perhaps you meant this:

Maximize(5*x1+2*x2+7*x3+6*x4+x5+6*x6+8*x7+6*x8,
         {10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8<=22},
         assume = nonnegint);

    [56, [x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0,
          x6 = 0, x7 = 7, x8 = 0] ]

If you state which solution you want (and what ranges for x and m) then you could adapt the rootfinding call appropriately.

The `plot` below could also be made more clear (filling in the gap). But you could also figure out how much accuracy you want.

The code has two ways to do the rootfinding (Calculus1:-Roots which utilizes fsolve's avoid option, and fsolve with its maxsols option which utilizes RootFinding:-NextZero).

restart;

N := q -> (exp(q)/sqrt(q))*sqrt(3/2)*Dw(sqrt(3*q/2));

proc (q) options operator, arrow; exp(q)*sqrt(3/2)*Dw(sqrt((3/2)*q))/sqrt(q) end proc

Dw := x -> exp(-x^2)*Int(exp(t^2), t=0..x);

proc (x) options operator, arrow; exp(-x^2)*(Int(exp(t^2), t = 0 .. x)) end proc

this := Dw(sqrt(3*q/2));

exp(-(3/2)*q)*(Int(exp(t^2), t = 0 .. (1/2)*6^(1/2)*q^(1/2)))

eq3a := q=(24/(n-3))*(x*k*m)/t;
eq3b := eval(eq3a, [t=6/10,n=6,k=13/10]);

q = 24*x*k*m/((n-3)*t)

q = (52/3)*x*m

eq4a := m = 1/2*((exp(q)/(q*N(q)))-1-1/q);

m = (1/6)*6^(1/2)/(q^(1/2)*exp(-(3/2)*q)*(Int(exp(t^2), t = 0 .. (1/2)*6^(1/2)*q^(1/2))))-1/2-(1/2)/q

eq4b := combine(value(eval(eq4a, eq3b))) assuming real;

m = (1/26)*26^(1/2)*exp(26*x*m)/(erfi((x*m)^(1/2)*26^(1/2))*(x*m*Pi)^(1/2))-(3/104)/(x*m)-1/2

M := proc(X)
  local res;
  option remember;
  res:=Student:-Calculus1:-Roots(eval(eq4b,x=X),m=-0.6 .. 0.6,numeric);
  #res:=[fsolve(eval(eq4b,x=X),m=-0.6 .. 0.6,maxsols=4)];
  res := [res[], undefined$(4-nops(res))];
end proc:
M1 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[1];
end proc:
M2 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[2];
end proc:
M3 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[3];
end proc:
M4 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[4];
end proc:

M(1.1);
M1(1.1), M2(1.1), M3(1.1), M4(1.1);

[-.4404638037, -0.4037552807e-1, .1382804516, .3786155189]

-.4404638037, -0.4037552807e-1, .1382804516, .3786155189

plot([M1(x), M2(x), M3(x), M4(x)], x = -3 .. 3, adaptive=false, numpoints=60);

plots:-implicitplot(rhs(eq4b)-lhs(eq4b), x=-3 .. 3, m=-0.5 .. 0.5, gridrefine=3);

 

Download dwson_m_x.mw

You can also do a similar thing for finding x as function of m. (It's actually easier and slightly quicker, avoids the gap, and the plot's axes could even be transposed).  dwson.mw

You could use an Empirical distribution or choose for this from RandomTools (or ProbabilityTable from the Statistics package).

I don't know whether you want individual (scalar) results, or lists, or Vectors. I don't know whether you want integers or floats. I don't know whether you'd prefer using Statistics or RandomTools, and I don't know whether you want easy ways to get just a few values or efficient ways to get many values.

So here are various ways.

restart;

with(RandomTools):

G := distribution(Empirical(<1,2,3,4,5>,
                            probabilities=[0.2,0.5,0.1,0.1,0.1])):

Generate(G);

HFloat(4.0)

seq(trunc(Generate(G)), i=1..10);

5, 1, 5, 2, 1, 2, 2, 5, 5, 1

L1 := Generate(list(choose([1$2,2$5,3,4,5]),10));

[2, 3, 2, 4, 5, 1, 3, 5, 4, 1]

# This is not fast
N := 10^5:
L2 := Generate(list(choose([1$2,2$5,3,4,5]),N)):

map(s->lhs(s)=evalf[2](rhs(s)/N), Statistics:-Tally(L2));

[1 = .20, 2 = .50, 3 = 0.99e-1, 4 = .10, 5 = .10]

with(Statistics):

H := RandomVariable(Empirical(<1,2,3,4,5>,
                              probabilities=[0.2,0.5,0.1,0.1,0.1])):

Sample(H, 1);

Vector(1, {(1) = 2.0})

trunc(Sample(H, 1)[1]);

3

V1 := trunc~(Vector(Sample(H, 10),datatype=anything));

Vector[row](10, {(1) = 1, (2) = 2, (3) = 2, (4) = 2, (5) = 3, (6) = 1, (7) = 3, (8) = 2, (9) = 5, (10) = 4})

L3 := trunc~(convert(Sample(H, 10),list));

[4, 4, 1, 2, 2, 2, 5, 3, 2, 1]

# This is not terribly slow
N := 10^5:
V2 := trunc~(Vector(Sample(H, N),datatype=anything)):

map(s->lhs(s)=evalf[2](rhs(s)/N), Statistics:-Tally(V2));

[1 = .20, 2 = .50, 3 = .10, 4 = .10, 5 = .10]

CodeTools:-Usage( Generate(list(choose([1$2,2$5,3,4,5]), 2*10^5)) ):

memory used=331.13MiB, alloc change=1.53MiB, cpu time=2.88s, real time=2.67s, gc time=439.41ms

CodeTools:-Usage( trunc~(Vector(Sample(H, 2*10^5),datatype=anything)) ):

memory used=19.85MiB, alloc change=6.11MiB, cpu time=68.00ms, real time=68.00ms, gc time=0ns

 

Download empirical_ac1.mw

You might also consider throwing in additional arrows normal to the line whose slope is the linear objective c1*x1+2*x2 .

One nice thing to illustrate is that the maximum is attained by moving [x1,x2] as far as possible (within the feasible region) in the direction of the normal to the objective line.

Below, I translate the line (so that it passes through the optimial point, give the current c1 value) since that gives a visual cue.

There are all kinds of variations on doing this.

If you really want to get fancy you could animate by the rotation angle (and compute the effective c1 value on the fly from that), which would produce a more even visual experience.

Another finesse might be to have a chunk of frames for the special case that the line is parallel to a side of the feasible region. In that case the whole edge of the region is optimal, and could be briefly all colored in green. Having a chunk of constant frames would provide a visual pause.

restart;

cnsts := [ 4*x1 + x2 <= 12,
             x1 - x2 >= 2,
                  x1 >= 0,
                  x2 >= 0 ]:

feasibleRegion := plots:-inequal(cnsts, x1 = 1 .. 4, x2 = -0.5 .. 1.5, nolines):

L := c1*x1+2*x2;

c1*x1+2*x2

m := solve(L,x2)/x1;

-(1/2)*c1

b := x2 - m*x1;

(1/2)*c1*x1+x2

Lform := Typesetting:-Typeset(Typesetting:-EV(L)):

 

All := proc(C1)
  local Nor,opt,P1,P2;
  uses plots, Optimization;
  opt := Maximize(eval(L,c1=C1), cnsts, x1=1..4, x2=-0.5..1.5);
  P1 := plot(eval(m*x1 + eval(b,opt[2]), c1=C1),
             x1=1..4, x2=-0.5..1.5, color=red,
             title=typeset('max'(L)=evalf[5](opt[1]),"\n",
                           "occurs at ",evalf[5](fnormal(opt[2])),"\n",
                           "when ", 'c1'=evalf[5](C1)));
  Nor := arrow(eval([x1,x2],opt[2]),
               <[C1,2]>, length=0.3, head_width=0.075, color=blue);
  P2 := pointplot([rhs~(opt[2])], color=green,
                  symbolsize=15, symbol=solidcircle);
  display(P1, P2, Nor, feasibleRegion,
          scaling=constrained);
end proc:

All(-0.4);

plots:-animate(All, [c1], c1=-20..20, paraminfo=false);

 

 

Download LP_anim2.mw

The basics of showing the feasible region, objective line, and optimal point are not complicated. There are several straightforward approaches to that, depending on how clear you want it to be for a reader who does not already know the underlying math.

Naturally, the code gets longer as you add in normal lines, a complicated title, etc.

Eventually you might discover that it's easier to learn and use a little Maple language syntax than to try and use the context-menus (Clickable Math, right-panel) for everything.

But here are some ways, including a little more use of the context-menus. (There are, of course, many other ways, including doing it fully by language commands.)

restart

with(LinearAlgebra)

tau := Matrix(3, 3, {(1, 1) = 200, (1, 2) = 0, (1, 3) = 0, (2, 1) = -50, (2, 2) = -800, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 200})

 

tau-lambda*(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1})) = Matrix(%id = 18446883947032617798)"(->)"-(-200+lambda)^2*(800+lambda)"(->)"-(-200+lambda)^2*(800+lambda) = 0"(->)"[[lambda = -800], [lambda = 200], [lambda = 200]]"(->)"results``NULL

sigma := map[2](eval, lambda, results)

[-800, 200, 200]

sigma[1]

-800

sigma[2]

200

sigma[3]

200

eval(lambda, results[1])

-800

tau-lambda*(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1})) = Matrix(%id = 18446883947011398582)"(->)"-(-200+lambda)^2*(800+lambda)"(->)"-(-200+lambda)^2*(800+lambda) = 0"(->)"eqn

R := [solve(eqn, lambda)]

[-800, 200, 200]

R[1]

-800

``

Download maple_ac.mw

You could try changing the syntax to, say,

   evalf(Int((r,a)->P(r)*a, [0..1, 0..1])):

or,

   evalf(Int('P'(r)*a, [r=0..1, a=0..1]));

Do you see that when you changed from single-argument to two-argument you also changed the calling sequence? You changed it to an invalid mashup of operator and expression form.

Alternatively, you could use P(r)*a as the expression, but you'd first have to modify the defn of P so that it returned unevaluated if passed symbol r rather than a number. One of your attempts this way encountered what's known as premature evaluation, where P received nonnumeric r and was not set up to deal with that. For example,

P := proc(r)
   if not r::numeric then return 'procname'(r); end if;
   ...proceed as usual...   
end proc

I prefer this last approach.

But you cannot sensibly use P*a as the integrand, where P is a procedure and a is an expression in one of more of the variables of integration.

If you upload your actual data file as an attachment in a Comment on the Question then I can easily show you all three methods at work.

As a side note, procedure P could be far more efficient if you wrote it to not have to re-form the data every time it gets called. I'd prefer to have the data, to check, before re-writing it all. You could substitute the actual data in lieu of a dummy name in the proc, or you could do it like, say,
   unapply('CurveFitting:-ArrayInterpolation'(...),numeric)
where use is again made of uneval quotes.

It's difficult to be precise when you haven't shown us the details of the actual data.

But I'm curious, did you want to escape a backslash in "VAGRANTASP\Adam" ?

Your final output seems like it may contain "VAGRANTASP\\Adam" with an escaped backslash. But your 2D Input contains a piecewise that compares against "VAGRANTASP\Adam" instead. Is that mismatch an oversight?

First 130 131 132 133 134 135 136 Last Page 132 of 336