mmcdara

705 Reputation

12 Badges

3 years, 73 days

MaplePrimes Activity


These are answers submitted by mmcdara

I believe that the first thing to do is to simplify your equations...
You could take inspiration from this:
 

NULL

restart:

# do not use with(CodeGeneration) but instead :
#  1/ with(CodeGeneration, Matlab)
#  2: at the end of your worksheet CodeGeneration:-Matlab(...)

# first simplification of dM1

dM1 := subs(sqrt(b^2+z^2)=R,  diff(M1(z), z) = (eval(((2*I)*exp(-sqrt(t^2*v^2+b^2)+(1/2*I)*v^2*t)/sqrt(t^2*v^2+b^2)-(2*I)*exp(-sqrt(t^2*v^2+b^2)-(1/2*I)*v^2*t)/sqrt(t^2*v^2+b^2))*t^2+2*h[1](z)*t+((-(4*I)*t/sqrt(t^2*v^2+b^2)-2)*exp(-sqrt(t^2*v^2+b^2)+(1/2*I)*v^2*t)+(-2+(4*I)*t/sqrt(t^2*v^2+b^2))*exp(-sqrt(t^2*v^2+b^2)-(1/2*I)*v^2*t))*t+h[2](z)+(2*t+(2*I)*t^2/sqrt(t^2*v^2+b^2))*exp(-sqrt(t^2*v^2+b^2)+(1/2*I)*v^2*t)+(2*t-(2*I)*t^2/sqrt(t^2*v^2+b^2))*exp(-sqrt(t^2*v^2+b^2)-(1/2*I)*v^2*t), t = z/v))/v):

print():

# an even more simplified expression

dM1 := diff(M1(z), z) = simplify(expand(rhs(dM1)));

 

diff(M1(z), z) = (2*h[1](z)*z+h[2](z)*v)/v^2

(1)

``


 

Download 1.mw

 

To understand where the abs(1, ...) comes from and to handle correctly the discontinuity of the dervative when 
(.9*sqrt(-c^2+1)-.4358898944*c)/(-.9*sqrt(-c^2+1)-.4358898944*c)) = 0, please look this


 

restart:

diff(log(abs(c)), c)

abs(1, c)/abs(c)

(1)

f := piecewise(c < 0, log(-c), c > 0, log(c))

f := piecewise(c < 0, ln(-c), 0 < c, ln(c))

(2)

diff(f, c)

piecewise(c = 0, undefined, 1/c)

(3)

In your case:

I3 := -sqrt(-c^2+1)*piecewise(F(c) < 0, log(-F(c)), F(c) > 0, log(F(c)), undefined)

I3 := -sqrt(-c^2+1)*piecewise(F(c) < 0, ln(-F(c)), 0 < F(c), ln(F(c)), undefined)

(4)

diff(I3, c)

piecewise(F(c) < 0, ln(-F(c)), 0 < F(c), ln(F(c)), undefined)*c/sqrt(-c^2+1)-sqrt(-c^2+1)*piecewise(F(c) < 0, (diff(F(c), c))/F(c), 0 < F(c), (diff(F(c), c))/F(c), undefined)

(5)

eval(%, F(c) = (.9*sqrt(-c^2+1)-.4358898944*c)/(-.9*sqrt(-c^2+1)-.4358898944*c)):

 

(6)

 


 

Download hint.mw

Unless I'm mistaken, the limit is +infinity if x > 1 and has a finite value if x < 1 (undefined for x=1):

 

restart:

J := Int(Int(exp(-(y-beta[0]-beta[1]*x-b0-b1*x)^2/(2*sigma^2))*exp(-b0)*exp(-b1),b0=0..infinity),b1=0..infinity)

Int(Int(exp(-(1/2)*(-b1*x-x*beta[1]-b0+y-beta[0])^2/sigma^2)*exp(-b0)*exp(-b1), b0 = 0 .. infinity), b1 = 0 .. infinity)

(1)

J0 := op(1,J)

Int(exp(-(1/2)*(-b1*x-x*beta[1]-b0+y-beta[0])^2/sigma^2)*exp(-b0)*exp(-b1), b0 = 0 .. infinity)

(2)

combine(op(1, J0), exp);
T := Student:-Precalculus:-CompleteSquare(%, b0)

exp(-(1/2)*(-b1*x-x*beta[1]-b0+y-beta[0])^2/sigma^2-b0-b1)

 

exp(-(1/2)*(b0-(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)*sigma^2)^2/sigma^2-(1/2)*(-b1*x-x*beta[1]+y-beta[0])^2/sigma^2-b1+(1/2)*(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)^2*sigma^2)

(3)

T_e := op(T);

-(1/2)*(b0-(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)*sigma^2)^2/sigma^2-(1/2)*(-b1*x-x*beta[1]+y-beta[0])^2/sigma^2-b1+(1/2)*(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)^2*sigma^2

(4)

Remainder := add(op(2..-1, T_e)):
Gauss := op(1, T_e):

J0 := exp(Remainder)*Int(exp(Gauss), b0=0..infinity);

exp(-(1/2)*(-b1*x-x*beta[1]+y-beta[0])^2/sigma^2-b1+(1/2)*(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)^2*sigma^2)*(Int(exp(-(1/2)*(b0-(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)*sigma^2)^2/sigma^2), b0 = 0 .. infinity))

(5)

# the integral is obviously equal to 1/(2*sigma*sqrt(2*Pi))

J0_sol := exp(Remainder)/(2*sigma*sqrt(2*Pi))

(1/4)*exp(-(1/2)*(-b1*x-x*beta[1]+y-beta[0])^2/sigma^2-b1+(1/2)*(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)^2*sigma^2)*2^(1/2)/(sigma*Pi^(1/2))

(6)

J1 := Int(J0_sol, b1=0..+infinity)

Int((1/4)*exp(-(1/2)*(-b1*x-x*beta[1]+y-beta[0])^2/sigma^2-b1+(1/2)*(-(1/2)*(2*b1*x+2*x*beta[1]-2*y+2*beta[0])/sigma^2-1)^2*sigma^2)*2^(1/2)/(sigma*Pi^(1/2)), b1 = 0 .. infinity)

(7)

J1 := expand(combine(J1, exp));

(1/4)*2^(1/2)*exp((1/2)*sigma^2)*exp(beta[0])*exp(beta[1]*x)*(Int(exp(b1*x)/exp(b1), b1 = 0 .. infinity))/(sigma*Pi^(1/2)*exp(y))

(8)

J := value(%);

(1/4)*2^(1/2)*exp((1/2)*sigma^2)*exp(beta[0])*exp(beta[1]*x)*(limit((exp(b1*x-b1)-1)/(-1+x), b1 = infinity))/(sigma*Pi^(1/2)*exp(y))

(9)

eval(J) assuming x < 1

-(1/4)*2^(1/2)*exp((1/2)*sigma^2)*exp(beta[0])*exp(beta[1]*x)/(sigma*Pi^(1/2)*exp(y)*(-1+x))

(10)

eval(J) assuming x > 1

2^(1/2)*exp((1/2)*sigma^2)*exp(beta[0])*exp(beta[1]*x)*infinity/(sigma*Pi^(1/2)*exp(y))

(11)

 


 

Download error_function.mw

Replace ; by : at the end of the command

  1. https://www.maplesoft.com/support/help/Maple/view.aspx?path=ProgrammingGuide/Chapter15
    next choose Multithread programming and/or Parallel programming (it depends on what type of "parallelization" you are interested in).
  2.  I'm sure there is a little two parts course on the Posts section of Mapleprimes (it was pretty well done... author dohashi ??? ).
    Unfortunately iI couldn't retrieve it.
  3. on Mapleprimes site :
    Parallel Programming Contents Page
    The Current Limitations of Parallel Programming In Maple

     

@radaar 
"Will Grid:-Map use a lot of memory? " : It depends on the way you have coded.
From experience, reducing the computation time with Grid is very easy, but avoiding used memory inflation can be a hard task.

The first thing to do is:

  1. Run your code on a single proc and catch the memory used from the task monitor (if you use Wondows).
    Look also to the runing process dubbed "mserver" and rnote the memory it uses ; let's notoe it S.
     
  2. Do the same by using Grid on P processors. You will have P+1 mserver processes: look to the memory size they use (depending on your programming you could have P processors with size S and one "the parent one" with size P*S: this means you will have to search in your code why it consumes so many memory when run on P processors).
     
  3. I don't have Windows right now to give you the correct command I use to use, but you can catch programmatically the total memory used (the command can be run with ssystem("command"))


More generally, if you are interested to use Grid, here are a few advices (from experience):

  1. Avoid using convert(blablabla, X) where X is a set, list, array, matrix, ...; from experience this can dramatically inflate the memory. So think first to the good structure blablabla must have.
     
  2. Avoid solve and prefer fsolve. I do not understand why solve tends to increase memory (maybe a point related to the remember table solve uses?) , but it seems to be the fact.
     
  3. grid:-Run and Grid:-Launch seem to do the same thing, but it's not true and a code optimized for one of them can give poor results when run with the other one.
     
  4. Forget garbage collector gc: it doesn't work on multiple processors
     

So, no real answers here, just be very careful: a code that requires a small amount of memory on a single processor can be very greedy on multiple processors (i already crashed my 64 Gb machine with "raw" coding when a good coding used less than 1Gb)

Here are two ways to do this. 
 


 

restart:

with(Statistics):

C := RandomVariable(Cauchy(0, 1)):

# The easiest way

N := 10^4:
C_sample := Sample(C, N, method=[envelope, range=-1..1]):
Histogram(C_sample);

[Mean, Variance](C_sample);

CodeTools:-Usage( Sample(C, N, method=[envelope, range=-1..1]) ):

 

[HFloat(-3.747279823038963e-4), HFloat(0.28029046512049793)]

 

memory used=507.56KiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns

 

# A lengthy way

p := Probability(C < -1, numeric);
q := 1-Probability(C > 1, numeric);

U_sample := Sample(Uniform(p, q), N):

C_cdf  := unapply(CDF(C, z), z):
C_icdf := unapply(solve(C_cdf(u)=z, u), z);

C_sample := C_icdf~(U_sample):
Histogram(C_sample);

[Mean, Variance](C_sample);

CodeTools:-Usage( C_icdf~(U_sample) ):

HFloat(0.25)

 

HFloat(0.75)

 

proc (z) options operator, arrow; -cot(z*Pi) end proc

 

 

[HFloat(-0.008152438124121577), HFloat(0.2737286335656307)]

 

memory used=2.21MiB, alloc change=0 bytes, cpu time=30.00ms, real time=30.00ms, gc time=0ns

 

 


 

Download CauchySampling.mw

Completing Carl Love's reply: if you want to pick "without replacement" and "without order constraints" (that is the nth picked number is not necessarily greater than the previous one) P numbers from a list L you must use combinat:-randperm
 

restart

L := [$1..100]:

P:= 3:
combinat:-randperm(L)[1..P]

[13, 71, 67]

(1)

# compare this

for k from 1 to 10 do
  combinat:-randcomb([$1..5], 5)
end do;

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

 

[1, 2, 3, 4, 5]

(2)

# to this

for k from 1 to 10 do
  combinat:-randperm([$1..5], 5)
end do;

[2, 4, 5, 1, 3]

 

[5, 2, 1, 4, 3]

 

[5, 1, 2, 3, 4]

 

[5, 1, 3, 4, 2]

 

[2, 4, 1, 3, 5]

 

[2, 4, 5, 1, 3]

 

[5, 4, 2, 3, 1]

 

[1, 5, 2, 3, 4]

 

[1, 2, 3, 4, 5]

 

[1, 3, 2, 5, 4]

(3)

 


 

Download pick.mw

If you really want to use colorscheme you can use this (the key point is style=point)

N  := numelems(x1):
xy := convert(Matrix(2, N, [x1, y1])^+,listlist):

plot(xy, style=point, colorscheme=["zgradient", ["Orange", "Red", "NavyBlue"]]);  # for instance


Another way to blavk/Blue coloring:

with(plots):
display(seq(pointplot([x1[k], y1[k]], color=piecewise(z1[k]<0.5, "Black", "RoyalBlue")), k=1..N))

 


 

restart

fe := Int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)-5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x = 0 .. al)

Int((1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)-5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2)^(1/2), x = 0 .. al)

(1)

# this to have a broad idea of a search interval

plot(evalf(subs(al=a, fe)), a=0..1, gridlines=true)

 

g := proc(a)
  evalf(subs(al=a, fe))-0.5
end proc;

proc (a) evalf(subs(al = a, fe))-.5 end proc

(2)

fsolve(g(al), al=0..1);

.1550088024

(3)

evalf(subs(al=%, fe))

.5000000002

(4)

 


 

Download int.mw

@radaar   
@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 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])

problems := Matrix(5, 1, {(1, 1) = M1, (2, 1) = M2, (3, 1) = M3, (4, 1) = M4, (5, 1) = M5})

(1)

# g := "negative effect on the health of the consumer"

PA := "probability occurrence of each problems"

"probability occurrence of each problems"

(2)

g := Matrix(1, 3, [1, 2, 3])

g := Matrix(1, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 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    
 

Criteria := Matrix(1, 3, {(1, 1) = "not serious", (1, 2) = "moderately serious", (1, 3) = "very serious"})

 

"not serious"

 

"moderately serious"

 

"very serious"

(4)

IN := Matrix(1, 3, [1, 2, 3])

IN := Matrix(1, 3, {(1, 1) = 1, (1, 2) = 2, (1, 3) = 3})

(5)

A := Matrix(1, 3, ["without risk", "occasional dange", "very common danger"])

A := Matrix(1, 3, {(1, 1) = "without risk", (1, 2) = "occasional dange", (1, 3) = "very common danger"})

(6)

for i to 3 do
  if IN(i) = i then
    print(A(i))
  end if
end do

"without risk"

 

"occasional dange"

 

"very common danger"

(7)

 


 

Download DontUnderstandYourProblem.mw

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

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 = 3, 2 = 2, 3 = 2, 4 = 1]

(1)

# works also for names

Set:={{a, b, 2}, {x, 1}, {a, x}, {4, 1, a, 3}};

Statistics:-Tally(op~([Set[]]))

{{1, x}, {a, x}, {2, a, b}, {1, 3, 4, a}}

 

[1 = 2, 2 = 1, 3 = 1, 4 = 1, a = 3, x = 2, b = 1]

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

{{1, x}, {x, a+x}, {2, b, sin(a)}, {1, 3, 4, a}}

 

[1 = 2, 2 = 1, 3 = 1, 4 = 1, a = 1, sin(a) = 1, x = 2, b = 1, a+x = 1]

(3)

 


 

Download tally.mw

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)

(G*(2-G))^(1/2)+(1-G)*(arccos(G-1)+(1/5)*Pi)+k*(G*(2-G))^(1/2)*cos(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)+(k/(1+k))^(1/2)*(1-G)*k*sin(((1+k)/k)^(1/2)*arccos(G-1)+(1/5)*Pi)

(1)

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

"AXESLABELS[Typesetting:-mi("G",italic = "true",mathvariant = "italic"), Typesetting:-mi("k",italic = "true",mathvariant = "italic")]"

 

Matrix(%id = 18446744078394417014)

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

HFloat(5.167125933526261)-HFloat(6.356106757782967)*G+HFloat(3.0043651079986984)*G^2-HFloat(0.767811803757935)*G^3

 

 

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

[0.15e-1, 0.3e-2, 0.5e-2, 0.2e-2, 0.2e-2, 0., 0.2e-2, -0.4e-2, -0.3e-2, -0.3e-2, -0.2e-2, 0.1e-2, -0.4e-2, -0.4e-2, -0.1e-2, -0.4e-2, -0.2e-2, 0.5e-3, -0.27e-2, -0.16e-2, -0.35e-2, -0.6e-3, 0., 0.2e-3, -0.16e-2, 0.8e-3, 0.5e-3, -0.22e-2, -0.6e-3, -0.20e-2, 0.8e-3, 0.10e-2, -0.5e-3, 0.16e-2, -0.6e-3, 0.42e-2, 0.42e-2, 0.30e-2, 0.95e-3, 0.309e-2, 0.246e-2, 0.280e-2, 0.157e-2, 0.433e-2, 0.3612e-2, 0.191e-2, 0.247e-2, 0.188e-2, 0.351e-2, 0.24e-3, 0.256e-2, 0.295e-2, 0.275e-2, 0.299e-2, 0.65e-3, 0.271e-2, 0.229e-2, 0.58e-3, 0.40e-3, 0.24e-2, 0.12e-2, 0.19e-2, 0.5e-3, 0.6e-3, 0.10e-2, -0.2e-3, -0.6e-3, -0.10e-3, -0.92e-3, -0.98e-3, -0.19e-3, -0.199e-2, -0.288e-2, -0.349e-2, -0.258e-2, -0.166e-2, -0.225e-2, -0.358e-2, -0.265e-2, -0.176e-2, -0.248e-2, -0.318e-2, -0.339e-2, -0.304e-2, -0.336e-2, -0.198e-2, -0.170e-2, -0.224e-2, -0.183e-2, -0.236e-2, -0.79e-3, -0.107e-2, -0.133e-2, -0.57e-3, 0.30e-3, 0.29e-3, 0.140e-2, 0.1598e-2, 0.2962e-2, 0.3399e-2, 0.4400e-2]

(3)

 


 

Download LevelCurve2.mw

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

 

 


 

Download contourplot.mw

1 2 3 4 5 6 7 Last Page 1 of 10