## 3810 Reputation

6 years, 114 days

## What do you mean by "surrogate optimiza...

@Carl Love

What do you mean by "surrogate optimization built in function"?
A surrogate can be as simple as the linear model Statistics:-LinearFit returns, or as complex as a neural network (Maple 2018 and later) with, in between GPE (gaussian process emulation such as kriging for instance).
Surrogate-based optimization generally uses GPE

As Carl wrote you could use the Kriging interpolation (I think this requires Maple 2018 or later).
But this implementation of the kriging method is more related to geostatistics than to surrogate modeling of an expensive code (mainly: correlation functions not well adpated).

There also exists a (commercial) toolbox dubbed "GlobalOptimization" which  is supposed to implement the EGO (Efficient Global Optimization method based on kriging surrogates).... unfortunately it gives poor results and doesn't work as claimed (in particular it can't be parameterized as described in the help pages; I posted a question about it a few weeks ago).

My personal (and thus questionable) opinion is that Maple is not yet adpated to do Surrogate-based optimization.

## I can't picture out what your proble...

I can't picture out what your problem really is...

1. Were did you find this procedure "print", is it one you wrote somewhere else?
2. Use Matrix instead of matrix with strng type elements.
3. Maybe read the help pages about Matrix and print?

 > restart:
 > problems := Matrix(5, 1, [M1, M2, M3, M4, M5]) (1)
 > # g := "negative effect on the health of the consumer"
 > PA := "probability occurrence of each problems" (2)
 > g := Matrix(1, 3, [1, 2, 3]) (3)
 > Criteria := Matrix(1,3,[ "not serious",  "moderately serious", "very serious"]); #` For each value of g we print the corresponding adjective from the matrix Criteria ` for i from 1 by 1 to 3 do   if  g(i)=i  then        print(Criteria(i));     end if    end do    (4)
 > IN := Matrix(1, 3, [1, 2, 3]) (5)
 > A := Matrix(1, 3, ["without risk", "occasional dange", "very common danger"]) (6)
 > for i to 3 do   if IN(i) = i then     print(A(i))   end if end do   (7)
 >

PS: if my answer is completely out of line, please forget it

## Tally......

An alternative:

For "Tally" see the corresponding help page.
It's ouput is a list of equalities, the lhs being the "name" of an element of Set and the rhs its number of occurences.

 > Set:={{1, 2}, {1, 3}, {1, 4}, {2, 3}}; Statistics:-Tally(op~([Set[]])) (1)
 > # works also for names Set:={{a, b, 2}, {x, 1}, {a, x}, {4, 1, a, 3}}; Statistics:-Tally(op~([Set[]]))  (2)
 > # and even in more general cases Set:={{sin(a), b, 2}, {x, 1}, {a+x, x}, {4, 1, a, 3}}; Statistics:-Tally(op~([Set[]]))  (3)
 >

## I do not think you will be able to obtai...

I do not think you will be able to obtain an analytic expression.
Here is an approximate one:

 > restart:
 > z := sqrt(G*(2-G))+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*sqrt(G*(2-G))*cos(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi)+sqrt(k/(1+k))*(1-G)*k*sin(sqrt((1+k)/k)*arccos(G-1)+(1/5)*Pi) (1)
 > p := plots:-contourplot(z, G=-1..2, k=1..10, contours=, grid=[400, 400]):
 > op(-1, p); c := op(1..-2, op(1, p)): M := convert(op~(1, [c]), Matrix);  (2)
 > plots:-display(p, gridlines=true) > # here is the approximate relation k=f(G) such that z(k=f(G), G) = 0) fit := Statistics:-LinearFit([1, G, G^2, G^3], M, [G]); plots:-display(p, plot(fit, G=min(M[..,1])..max(M[..,1]), color=blue, thickness=7, transparency=0.8, gridlines=true))  > zg := eval(z, k=fit):  # = z(k=fit(G), G): fit_error := evalf(map(u -> eval(zg, G=u), [seq(0..1, 0.01)])) (3)
 >

## Here is a solution where each level curv...

Here is a solution where each level curve is a colred differently

Maybe you would also have a legend for each level curve?
In this case you could ne interested by what acer did here 2D contour plot and legend

PS I removed the ouput for the uploaded display doesn't have colored level curves (but they are)

 > restart:
 > with(plots):
 > f := sin(x*y): r := x = -2 .. 2, y = -2 .. 2: cols := [seq(ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]),i=1..19)]: plots:-display(   plot3d(f, r, axes = framed, style=surface, color=gray),   seq(contourplot3d(f, r, contours=[(j-10)/10], color=cols[j]), j=1..19) );

 >

## Luckily I had already done it, I just wa...

Luckily I had already done it, I just wasted a little time finding the code.
From memory I believe I recoded in Maple the code from R (not completely sure).

Here it is.
I remember having done some tests and it seemed ok by the time , but feel free to inform me if it's not.

 > restart:
 > LSD_Sampling := proc(alpha, N)   local c, V, X, n, u, q:   uses Statistics:   c     := log(1-alpha):   V     := Sample(Uniform(0, 1), N):   X     := Vector[column](N):   for n from 1 to N do      if V[n] > alpha then         X[n] := 1:      else         u := Sample(Uniform(0, 1), 1):         q := 1 - exp(c*u):         if V[n] < q^2 then            X[n] := ceil(1+log(V[n])/log(q)):         elif V[n] < q then            X[n] := 1:         else            X[n] := 2:         end if:      end if:   end do:   return X: end proc:
 > Statistics:-Histogram(LSD_Sampling(0.9, 10000)) >

BTW: I'm curious, what application of the Logarithmic Series distribution are you interested in?

## Even Re(code) doesn(t work. The reason i...

Even Re(code) doesn't work. The reason is: Maple doesn't know if the variables are real or complex !

Assuming all the quanties are real a simple way is

subs(J=I, collect(subs(I=J, code), J));

Other possibility: use assume, or assuming (see help pages)

## I already ask this question a few years ...

I already ask this question a few years ago, you can trace the answers here

The solution I use:

1. I assume you have 2 mw files, let's say File1.mw and File2.mw whoch correspond to different versions of the same code.

2. Open File1.mw and Export it to maple input format (I choose the same name as the mw file)

3. Open File2.mw and Export it to maple input format; two files File1.mpl and File2.mpl are created

4. If not already installed, download NotePad++ (free on Windows ; Notepad++ has contextual coloring for Maple syntax)

5. Open NotePad++ and load File1.mpl and File2.mpl

6. In the menu of NotePad++ choose Run > Compare (this from memory because I do not have NotePad++ on this laptop)
The differences between the two versions are then easily visble

Once done, you have the choice to keep developping with the Java interface of Maple (as usual), or to develop "within" NotePad++.
In this latter case you will have to load an mpl file in a Maple session to verify that its content gives the expected result.
It's up to you, I personally prefer developping "within" a worksheet, saving in a new mw file, and exporting it in a new mpl file.

PS : there exist probably tools apart NotePad++ that can "understand" Maple's syntax and propose contextual coloring.

## I would like to complete acer's answ...

I would like to complete acer's answer. Personally I like to use notations like theta__i instead of theta[i], because i is not only limited to integer values.
But using theta__i instead of theta[i] requires some precautions.
You will find below a few examples.

Be careful, when used in a procedure, variables of the form theta__||n (last example) are not implicitely considered as local but as global: this may cause unwanted behaviours.

 > # differences : a := theta: b := theta__2: whattype(a), whattype(b); printf("%a, %a", theta, theta__2); print(a,b); theta, theta__2 (1)
 > # use of theta__  ; first case: wrong way restart:
 > vars := seq(theta__i, i=1..4);
 > mean := Array(1..4, i -> diff(w(vars), theta__i));
 >  (2)
 > # use of theta__  ; second case: good way restart:
 > vars := seq(theta__||i, i=1..4);
 > mean := Array(1..4, i -> diff(w(vars), theta__||i));
 >  (3)
 > # Why using theta__ instead of theta[..] can be useful? # Because the index n in theta__n can be of non numeric type: restart: MyIndices := [a, b, c, d]; vars := seq(theta__||n, n in MyIndices);
 > mean := Array(1..4, i -> diff(w(vars), vars[i]));   (4)
 > restart:
 > f := proc() local n: for n from 1 to 4 do   theta__||n := n end do: end proc:
 > f(): seq(theta__||i, i=1..4); (5)
 >

## Here are two solutions :restart:vars := ...

Two ways (probably among others)

restart:
vars := seq(theta[i], i=1..4);
mean := Array(1..4, i -> diff(w(vars), theta[i]));

# or
mean := Array(1..4, i -> diff(w(vars), vars[i]));

## In case you would be interested, give a ...

In case you would be interested, give a lookk to this site :

https://oeis.org/?language=english

Next just copy-paste 2,3,4,6,8,9,10,12,14,15,16,18 in the search field and enjoy

## I personally prefer to define f, f2...

I personally prefer to define f, f2 and f3 defore instanciate x and y (even if the result is the same).
The main point is in the definition of f3 which should had been f3 := (x, y) -> f(x, y)+f2(x, y)

 > restart:
 > f  := (x,y) -> x+y; f2 := (x,y) -> x^2+y^3; f3 := (x,y) -> f(x,y)+f2(x,y);   (1)
 > f(x,y) ; f2(x,y) ; f3(x,y)   (2)
 > x := 5; y := 10; f(x, y) ; f2(x, y) ; f3(x, y)     (3)
 > x := 'x': y := 'y': x; y; f(x, y) ; f2(x, y) ; f3(x, y)     (4)
 > restart:
 > x := 5; y := 10;  (5)
 > f  := (x,y) -> x+y; f2 := (x,y) -> x^2+y^3; f3 := (x,y) -> f(x,y)+f2(x,y);   (6)
 > f(x, y) ; f2(x, y) ; f3(x, y)   (7)
 >

## Maybe this?   ...

Maybe this?

 > restart
 > u := sin(x+y):
 > MyPlot := plot3d(u, x = 0 .. 2*Pi, y = 0 .. 2*Pi):
 > plottools:-getdata(MyPlot); (1)
 > M := plottools:-getdata(MyPlot)[-1]; (2)
 > # search Export in the help pages for other formats MyFile := FileTools:-JoinPath(["Desktop/Maple/M.csv"], base=homedir): Export(MyFile, M): (3)
 >

## I believe Christopher has put the finger...

I believe Christopher has put the finger on the problem when saying "...  is it the y portion? "

The structure of the pde is rather simple with the form F . grad(Phi) = 0 (I use Phi instead of w).
For these equations the key is in finding the characteristic curves. The problem seems to be: Maple cannot find a closed form expression for the second component of these curves.

(Nota: Maple find the characteristic curves in closed form if alpha=0 ... and successes too in solving the pde)

 > restart:
 > pde :=  diff(Phi(x,y,z),x)         +         (y^2- a*exp(alpha*x)*(x*y-1))*diff(Phi(x,y,z),y)         +         (c*exp(beta*x)*z^2+b*exp(-beta*x))*diff(Phi(x,y,z),z)         =         0; (1)
 > # Write the pde as F . grad(Phi) = 0 where F is the vector (A, B, C) # Next parameterize F as a function of s : F(s) = (A(s), B(s), C(s)) A := s -> 1; B := s -> (y(s)^2- a*exp(alpha*x(s))*(x(s)*y(s)-1)); C := s -> (c*exp(beta*x(s))*z(s)^2+b*exp(-beta*x(s)));   (2)
 > # Define the characteristic curve associated to the pde as the curve defined by (U(s), V(s), W(s)) # where : #   diff(U(s), s) = A(s) #   diff(V(s), s) = B(s) #   diff(W(s), s) = C(s) # # Then Phi(s) is constant along each characteristic curve # # Characteristic curve, component 1 eq1 := diff(U(s), s) = A(s); dsolve(eq1, U(s)): U := unapply(rhs(%), s);  (3)
 > # Characteristic curve, component 3 eq3 := subs({x(s)=U(s), z(s)=W(s)}, diff(W(s), s) = C(s)); dsolve(eq3, W(s))  (4)
 > # Characteristic curve, component 2 eq2 := subs({x(s)=U(s), y(s)=V(s)}, diff(V(s), s) = expand(B(s))); infolevel[dsolve] := 4: dsolve(eq2, V(s)) Methods for first order ODEs: --- Trying classification methods --- trying a quadrature trying 1st order linear trying Bernoulli trying separable trying inverse linear trying homogeneous types: trying Chini differential order: 1; looking for linear symmetries trying exact Looking for potential symmetries trying Riccati trying inverse_Riccati trying 1st order ODE linearizable_by_differentiation --- Trying Lie symmetry methods, 1st order ---  -> Computing symmetries using: way = 4  -> Computing symmetries using: way = 2  -> Computing symmetries using: way = 6 trying symmetry patterns for 1st order ODEs -> trying a symmetry pattern of the form [F(x)*G(y), 0] -> trying a symmetry pattern of the form [0, F(x)*G(y)] -> trying symmetry patterns of the forms [F(x),G(y)] and [G(y),F(x)]  -> Computing symmetries using: way = HINT    -> Calling odsolve with the ODE diff(y(x) x) = 2*y(x)/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x) = y(x)/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+y(x)*alpha y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+exp(alpha*(x+_C1))*a*K*(x+_C1) y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       <- quadrature successful    -> Calling odsolve with the ODE diff(y(x) x)+(y(x)*_C1*alpha+y(x)*alpha*x+y(x)+2*K)/(x+_C1) y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+K y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+y(x)*alpha-K y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+y(x)*(_C1*alpha+alpha*x+1)/(x+_C1) y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x) = 0 y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x) = -y(x)*alpha y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)-2*y(x)/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)-(x*alpha*K+y(x))/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)-(x*_C1*alpha*K+y(x)*_C1+x*K-K*alpha)/(_C1*x-1) y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)-y(x)/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)-y(x)*_C1/(_C1*x-1) y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful    -> Calling odsolve with the ODE diff(y(x) x)+(exp(alpha*_C1)*a*K-2*y(x))/x y(x)       *** Sublevel 2 ***       Methods for first order ODEs:       --- Trying classification methods ---       trying a quadrature       trying 1st order linear       <- 1st order linear successful -> trying a symmetry pattern of the form [F(x),G(x)] -> trying a symmetry pattern of the form [F(y),G(y)] -> trying a symmetry pattern of the form [F(x)+G(y), 0] -> trying a symmetry pattern of the form [0, F(x)+G(y)] -> trying a symmetry pattern of the form [F(x),G(x)*y+H(x)] -> trying a symmetry pattern of conformal type
 > # Note there is no explicit solution found: I suppose this is the reason why Maple cannot solve the pde
 > eq2 := subs({x(s)=U(s), y(s)=V(s), alpha=0}, diff(V(s), s) = expand(B(s))); infolevel[dsolve] := 4: dsolve(eq2, V(s)) Methods for first order ODEs: --- Trying classification methods --- trying a quadrature trying 1st order linear trying Bernoulli trying separable trying inverse linear trying homogeneous types: trying Chini Chini's absolute invariant is: a*(s+_C1)^2 differential order: 1; looking for linear symmetries trying exact Looking for potential symmetries Found potential symmetries: [0 1] [exp(-V)/a exp(-V)*(s+_C1)] Proceeding with integration step. (5)
 > pdsolve(subs(alpha=0, pde), Phi(x,y,z))
 Methods for first order ODEs: --- Trying classification methods --- trying a quadrature trying 1st order linear trying Bernoulli trying separable trying inverse linear trying homogeneous types: trying Chini Chini's absolute invariant is: a*x^2 differential order: 1; looking for linear symmetries trying exact Looking for potential symmetries Found potential symmetries: [0 1] [exp(-y) exp(-y)*a*x] Proceeding with integration step. Methods for first order ODEs: --- Trying classification methods --- trying a quadrature trying 1st order linear trying Bernoulli trying separable trying inverse linear trying homogeneous types: trying Chini Chini's absolute invariant is: beta^2/(c*b) <- Chini successful Methods for first order ODEs: --- Trying classification methods --- trying a quadrature trying 1st order linear <- 1st order linear successful (6)
 >

## Maybe this ?Watchout: you did not say wh...

Maybe this ?

Watchout: you did not say what to do if A inter B is the empty set

 > restart:
 > p := proc(A::set, B::set, C::set)    local AB, ABC:    AB := A intersect B:    if AB <> {} then       ABC := AB intersect C:       if ABC <> {} then          return ABC       else          return `union`(seq(AB *~ C[n], n=1..nops(c)))       end if;    else       error  "It's not specified what to do when A inter B = {} "    end if: end proc:
 > a := {x6,x4,x2,x7,x8,x9,x10}:   b := {x2,x3,x5,x8}: c := {x4,x9,x11,x12,x13}:
 > p(a, b, c) (1)
 > a := {x6,x4}:   b := {x5,x8}: c := {x12,x13}: p(a, b, c)
 > a := {x6,x4}:   b := {x4, x5,x8}: c := {x4, x12,x13}: p(a, b, c) (2)
 >