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)

p := plots:-contourplot(z, G=-1..2, k=1..10, contours=[0], grid=[400, 400]):

op(-1, p);
c := op(1..-2, op(1, p)):
M := convert(op~(1, [c]), Matrix);

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[4](map(u -> eval(zg, G=u), [seq(0..1, 0.01)]))

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)
);  

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)[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))

# differences :

a := theta[2]:
b := theta__2:

whattype(a), whattype(b);

printf("%a, %a", theta[2], theta__2);
print(a,b);

theta[2], theta__2

# 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));

# 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));

# 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]));

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);

restart:

f  := (x,y) -> x+y;
f2 := (x,y) -> x^2+y^3;
f3 := (x,y) -> f(x,y)+f2(x,y);

f(x,y) ; f2(x,y) ; f3(x,y)

x := 5;
y := 10;
f(x, y) ; f2(x, y) ; f3(x, y)

x := 'x':
y := 'y':
x;
y;
f(x, y) ; f2(x, y) ; f3(x, y)

restart:

x := 5;
y := 10;

f  := (x,y) -> x+y;
f2 := (x,y) -> x^2+y^3;
f3 := (x,y) -> f(x,y)+f2(x,y);

f(x, y) ; f2(x, y) ; f3(x, y)

restart

u := sin(x+y):

MyPlot := plot3d(u, x = 0 .. 2*Pi, y = 0 .. 2*Pi):

plottools:-getdata(MyPlot);

M := plottools:-getdata(MyPlot)[-1];

# search Export in the help pages for other formats

MyFile := FileTools:-JoinPath(["Desktop/Maple/M.csv"], base=homedir):
Export(MyFile, M):

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;

# 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)));

# 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);

# Characteristic curve, component 3


eq3 := subs({x(s)=U(s), z(s)=W(s)}, diff(W(s), s) = C(s));

dsolve(eq3, W(s))

# 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[1]*(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[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)+K[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)*alpha-K[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)*(_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[1]+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[1]+y(x)*_C1+x*K[1]-K[1]*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[1]-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.

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

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)

a := {x6,x4}:  
b := {x5,x8}:
c := {x12,x13}:

p(a, b, c)

Error, (in p) It's not specified what to do when A inter B = {}

a := {x6,x4}:  
b := {x4, x5,x8}:
c := {x4, x12,x13}:

p(a, b, c)

restart:

with(Statistics):
with(combinat):

# How is the game managed (I do not know about it)?
#
# Strategy 1:
# One decides to form the red team first by randomly picking 5 out of 15 members
# Next one decides to form the green team by randomly picking 5 out of the 10 remaining members
# The 5 remaining ones form the last (orange) team
#
#
# Strategy 2
# One randomly form 3 new teams of 5 members and randomly decide to color thes teams as red, green, orange
 

ProbabilityToPutFiveBluesInTheRedTeam := mul((6-n)/(15-n), n=0..4);
ProbabilityToPutTheLastBlueInTheOrangeTeam := 1/2;

ProbabilityStrategy1 := ProbabilityToPutFiveBluesInTheRedTeam * ProbabilityToPutTheLastBlueInTheOrangeTeam ;

# simulation strategy 1

OriginalBlueTeam   := [B$6]:
OriginalYellowTeam := [Y$9]:
candidates         := [OriginalBlueTeam[], OriginalYellowTeam[]]:


N := 10^6:


NewTeams := Matrix(N, 3):
for n from 1 to N do
   x := randperm(candidates):
   NewTeams[n, 1] := numelems(select(has, x[ 1.. 5], B)):
   NewTeams[n, 2] := numelems(select(has, x[ 6..10], B)):
   NewTeams[n, 3] := numelems(select(has, x[11..15], B)):
end do:


DesiredNumberOfBluesInNewTeams := [5, 0, 1]; # teams ordered red, green, orange

FrequencyStrategy1 := add( map(n -> if [entries(NewTeams[n,..], nolist)] = [5, 0, 1] then 1 end if, [$1..N]) ) / N

evalf([ProbabilityStrategy1, FrequencyStrategy1]);

ProbabilityToPutFiveBluesInTheSameTeam   := mul((6-n)/(15-n), n=0..4) * 3;
ProbabilityToPutTheLastBlueInAnotherTeam := 1;

ProbabilityStrategy2 := ProbabilityToPutFiveBluesInTheSameTeam * ProbabilityToPutTheLastBlueInAnotherTeam ;

# simulation strategy 2

OriginalBlueTeam   := [B$6]:
OriginalYellowTeam := [Y$9]:
candidates         := [OriginalBlueTeam[], OriginalYellowTeam[]]:


N := 10^6:


NewTeams := Matrix(N, 3):
for n from 1 to N do
   x := randperm(candidates):
   NewTeams[n, 1] := numelems(select(has, x[ 1.. 5], B)):
   NewTeams[n, 2] := numelems(select(has, x[ 6..10], B)):
   NewTeams[n, 3] := numelems(select(has, x[11..15], B)):
end do:


DesiredPartitions := { permute([5, 0, 1])[] };

FrequencyStrategy2 := add( map(n -> if member([entries(NewTeams[n,..], nolist)], DesiredPartitions) then 1 end if, [$1..N]) ) / N;

evalf([ProbabilityStrategy2, FrequencyStrategy2]);