mmcdara

how to manage the informations saved dur...

Asked: mmcdara 7891 Product: Maple

November 26 2020

1 3

Hi,

When I do large computations in a worksheet and that I manually save it, I get a message in a pop-up window which asks me if I want to save these calculations too:

What is it that is really saved if I answer "YES"
As I always answer "NO" how can I force MAPLE to ApplyGlobally this answer (I wasn't capable to fin it out in the Preferences>Interface menu).
is this preference accounted for during AutoSaving?

Always in the goal to save only the necessary informations (in a mw format, export to mpl format could of course be a solution), how can I remove all the outpouts (and more specifically all the plots) in the current worksheet?
Presently I'm doing this and, before my last saving, I turn TF to false and rexecute all _display commands ... not very practical, I agree

restart

_display := proc(tf, graph)
if tf then print(graph) end if:
end proc;

TF := true: # or false if no plot is required

# for each plot:
p := display( something ):
_display(TF, p)

Thanks in advance

Can Maple represent 2D density plots?...

Asked: mmcdara 7891 Product: Maple

November 24 2020

1 10

Suppose you do this

restart:
with(plots):
with(Statistics):
A := Sample(RandomVariable(Normal(0, 1)), 10^4):
B := Sample(RandomVariable(Normal(0, 1)), 10^4):
ScatterPlot(A, B, symbol=point);

The plot shows a greater density of dots around (0, 0) and they become increasingly sparse as one moves away from (0, 0).
A pretty way to represent such clouds is to do "2D density plots".
You will find some example from R here 2d-density-plot-with-ggplot2.html

For the moment I do something like this (which is to be adpated for almost any particular case).
My main question is: do recent versions of MAPLE possess something similar to R's features?

>

restart:

>	interface(version);

(1)

>	with(Statistics):

>

UseHardwareFloats := false:

N := 10^5:
A := Sample(RandomVariable(BetaDistribution(2, 6)), N):
B := Sample(RandomVariable(BetaDistribution(6, 2)), N):

nb := 50:

mA := min(A):
hA := Range(A)/(nb-1):
mB := min(B):
hB := Range(B)/(nb-1):

cA := floor~((A-~mA)/~hA)+~1:
cB := floor~((B-~mB)/~hB)+~1:

NP := Matrix(nb$2, 0):
for n from 1 to N do
NP[cA[n], cB[n]] := NP[cA[n], cB[n]]+1
end do:

NP:

>	plots:-matrixplot( NP, heights=histogram, colorscheme=["Blue", "White", "Red"], orientation=[270, 0, 0], style=patchnogrid, lightmodel=none, labels=["A", "B", ""], axis[1]=[tickmarks=[nb+1=max(B), 1=mB]], axis[2]=[tickmarks=[nb+1=max(A), 1=mA]] )

>

Download 2D_density_plot.mw

How to manage the number of colors color...

Asked: mmcdara 7891 Product: Maple 2015

November 13 2020

3 1

Hi,

I often need to draw several curves on the same figure and the question of which colors to choose to distinguish them is crucial.
In particular when these curves depend on a parameter I would like to color them by using a "regular" change of color.
The trick I use is this one

# here I want to plot N=10 different curves

N := 10:

# step 1: generate a list of colours using colorscheme
# (I was hoping for 10 levels of colors but I got 19)

plot(x, x=0..1, colorscheme=["Gold", "Blue"], numpoints=N):
MyColors := op([1,2, 2], %);
M := numelems(MyColors[..,1]);

# step 2: do the plots of interest, for instance

plots:-display(
  seq(
    plot(m*x^2, x=0..1, color=ColorTools:-Color([entries(MyColors[round(m/N*M)], nolist)])), 
    m=1..N
  )
);

My question is: How can we make that MyColors be a matrix whose the number of rows equal N (this to avoid the round(m/N*M) operation)?
(I will accept any other strategy)

Thanks in advance

How to plot the difference between 2 ds...

Asked: mmcdara 7891 Product: Maple 2015

November 05 2020

1 10

Hi,

I have an ODE system with a free parameter P, that I can only solve numerically.
To simplify let's suppose this parameter has only 2 values P1 and P2.
Let SOL1 the solution of this system when P=P1 and SOL2 its solution for P=P2:

Again for a sake of simplification let's suppose the independent variable is t and the dependent one is x(t).

How can I plot the difference x1(t)-x2(t) where x1(t) is the solution corresponding to SOL1 and x2(t) the solution corresponding to SOL2
(sorry I screwed something up and changed the config of my keyboard making the interrogation point unaccessible)

Here is a detailed example

Thanks in advance

>

restart:

>	with(plots):

>

# This is a 1D problem.
# A body with unit mass is submitted to an acceleration and to a solid-solid friction resistance
# The friction model is a simple Coulomb's

sys := {
         diff(x(t), t) = v(t),
         diff(v(t), t) = t - piecewise(v(t) > 0, f, v(t) < 0, -f, 0),
         x(0) = 0,
         v(0) = 0
       };

sol__0 := dsolve(sys, numeric, parameters=[f], maxfun=10^6)

sys := {diff(v(t), t) = t-piecewise(0 < v(t), f, v(t) < 0, -f, 0), diff(x(t), t) = v(t), v(0) = 0, x(0) = 0}

(1)

>	# This system doesn't seem to be solved formally sol__exact := dsolve(sys, {v(t), x(t)})

>

# Attempt of solving it numerically.
#
# Let's observe that the "v" solution must be an increasing function of t
# so the friction model could be replaced by piecewise(v(t) > 0, f, 0) for v(t) >=0.
# But this is a particular case (think of a sine acceleration for instance) and
# I keep using the initial Coulomb's model.
#
# The main problem, as it appears below is the discontinuity of the friction effect at v(t)=0
# (thus at the initial time).

sol := dsolve(sys, numeric, parameters=[f], maxfun=10^6):
sol(parameters=[1]):
sol(0.1);

Error, (in sol) cannot evaluate the solution further right of 0.15837836e-1, maxfun limit exceeded (see ?dsolve,maxfun for details)

>

# Clearly increasing maxfum will help to simulate on a longer time but it's not the
# most efficient way to handle this problem.
#
# The classical way is to use a regularized Coulomb's model ; for instance

sys := {
             diff(x(t), t) = v(t),
             diff(v(t), t) = t - f *tanh(v(t)/c),
             x(0) = 0,
             v(0) = 0
           };

sol__1 := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8)

(2)

>

# The smaller c the closer the solution to the one invoking the Coulomb's model.
#
# The question is: how can we measure as close to sol__exact we are?
# The idea is to choose a "small" value of c, lets say c__0, and to declare that the solution
# obtained with c__0 is "the exact solution".
#
# Studying the convergence of the solution as c__0 becomes smaller and smaller will help
# choosing c__0.
# Of course the computational time is an increasing function of 1/c__0 : the smaller c__0
# the higher this time.

c__0 := 1e-8: # arbitrary choice

sol__ref := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8):
sol__ref(parameters=[1, c__0]):

tmax := 10:

CodeTools:-Usage(sol__ref(tmax)):
numfun__0 := sol__ref(numfun);

odeplot(sol__ref, [t, v(t)], t=0..tmax)

memory used=33.20KiB, alloc change=0 bytes, cpu time=34.76s, real time=36.19s, gc time=0ns

>

# Assuming c__0 = 1e-8 is a value small enough to consider the solution plotted above
# is "the solution for the true Coulomb's model".
#
# Here are a few solutions for larger values of c:

cs   := 10^~[$-6..-1]:
tmax := 10:

printf("c = 1e-08 number of calls = %d\n", numfun__0):

for i from 1 to nops(cs) do
  sol__1(parameters=[1,cs[i]]):
  sol__1(tmax):
  couleur := ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):
  Xplot__||i := odeplot(
                         sol__1, [t, x(t)], t=0..tmax,
                         legend=cs[i],
                         color=couleur
                       ):
  printf("c = %1.0e number of calls = %d\n", cs[i], sol__1(numfun))
end do:

LY := plottools:-transform((x, y) -> [x, log(y)]):
LXY := plottools:-transform((x, y) -> [log(x), log(y)]):

display(
         seq(Xplot__||i , i=1..nops(cs)),
         view=[0..3, 0..2], labels=[t, x], title="x(t) (lin-lin)"
       );

c = 1e-08 number of calls = 50831166

c = 1e-06 number of calls = 1308504

c = 1e-05 number of calls = 126848
c = 1e-04 number of calls = 12825
c = 1e-03 number of calls = 1513
c = 1e-02 number of calls = 456
c = 1e-01 number of calls = 246

>

# From this plot it seems that values of c less than 1e-4 give roughly the same curves.
#
# Remembering that the cimputational time is an increasing function of 1/c, we face the
# practical dilema of getting a solution "very close" to the ideal one (Coulomb's model)
# in a reasonable amount of time.
# Thus we need a small value of c nut not too small.
#
# The natural idea is to plot the difference between x(t) obtained with some value of c
# and x(t) obtained with our "reference value" c = c__0.
#
# How can we do this efficiently?
# (unfortunately piecewise and array/Array outputs are not allowed for parameterized solutions)
#
# For instance this is very slow

gap := proc(c, s1, s2, n)
  local s, h, k:
  sol__1(parameters=[1, c]):
  sol__1(s2):
  h := Matrix(n, 2):
  k := 1:
  for s from s1 to s2 by (s2-s1)/(n-1) do
    h[k, 1] := s:
    h[k, 2] := eval(x(t), sol__1(s)) - eval(x(t), sol__ref(s)):
    k       := k+1:
  end do:
  return h
end proc:

GAP := CodeTools:-Usage( gap(1e-5, 0, 1, 11) ):

Statistics:-ScatterPlot(GAP[..,1], GAP[..,2], symbol=solidcircle, color=blue)

memory used=399.35KiB, alloc change=0 bytes, cpu time=34.34s, real time=34.81s, gc time=0ns

>

Download odeplot_differences.mw

How to fasten this simulation?...

Asked: mmcdara 7891 Product: Maple 2015

October 30 2020

2 3

Hi,

I want to simulate a decision rule (let's say a production control strategy to fix the ideas) based on the outcomes of a binomial random variable.

One stage of this simulation involves operations which correspond to the NaiveProc described below. This procedure is a straightforward and rather naive translation in Maple of these operations.
It is given here to help you to understand what is to be done.

This coding being extremely slow, I implemented the required operations differently.
Given the fact that I'm only concerned by binomial samplings with outcomes 0 and 1, I basically replaced these samplings by three operations U:=... , Got01:=... , Got0:= ... .
Each new implementation is about 200 times faster than the naive one.

My question is: Can we do still better?
For information the same simulation realized with R runs in less than 15s when REP = 10^8 (more than 100 times faster than what I'm capable to do with Maple)

Thanks for your help.

>

restart:

>	with(Statistics):

>

# This is the naive coding of the problem
#
# It has the advantage to explain clearly (at least I hope so) what is to be done.

NaiveProc := proc(REP, a, b, N)
  local roll, rep, p, K, N0, N1:
  roll := rand(a .. b);
  N0   := 0:
  N1   := 0:
  for rep from 1 to REP do
    p := roll();
    K := Sample(Binomial(N, p), 1)[1]:
    if K=0 then
      N0 := N0+1:
    elif K=1 then
      N1 := N1+1:
    end if:
  end do:
  N0, N1:
end proc:

CodeTools:-Usage(NaiveProc(10^3, 0, 5e-4, 25000));

memory used=78.92MiB, alloc change=68.00MiB, cpu time=1.21s, real time=6.30s, gc time=24.80ms

(1)

>

# As the previous code is very slow I implemented several variants based on the fact
# that I care only on values of K equal to 0 and 1 and that it is then useless to sample
# a Binomial random variable to do this.
#
# Each of them is much faster than the naive coding... but can we even faster?
# For comparison, the same simulation realized in R takes less than 15s to run with REP=10^8

pi := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k;

f1 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Select(t -> is(t <= 0), U-P01):
  Got0 := Select(t -> is(t <= 0), U-P0 );
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f1(10^5, 0, 5e-4, 25000));

print():

f2 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := select[flatten](`<=`, U-P01, 0):
  Got0 := select[flatten](`<=`, U-P0 , 0):
  N01   := numelems(Got01):
  N0    := numelems(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f2(10^5, 0, 5e-4, 25000));

print():

f3 := proc(REP, a, b, N)
  local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:
  Q     := Sample(Uniform(a, b), REP)^+:
  P0    := pi~(N, Q, 0):
  P1    := pi~(N, Q, 1):
  P01   := P0 + P1:
  U     := Sample(Uniform(0, 1), REP)^+:
  Got01 := Vector(REP, i -> `if`(U[i] <= P01[i], 1, 0)):
  Got0 := Vector(REP, i -> `if`(U[i] <= P0 [i], 1, 0)):
  N01   := add(Got01):
  N0    := add(Got0):
  N1    := N01-N0:
  N0, N1;
end proc:

CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));

memory used=376.83MiB, alloc change=296.77MiB, cpu time=4.12s, real time=4.03s, gc time=430.38ms

memory used=173.26MiB, alloc change=0 bytes, cpu time=2.30s, real time=1.98s, gc time=461.22ms

memory used=154.19MiB, alloc change=0 bytes, cpu time=2.28s, real time=2.01s, gc time=432.40ms

(2)

>

Download MakeItFaster.mw

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

7891 Reputation

22 Badges

MaplePrimes Activity

These are questions asked by mmcdara

how to manage the informations saved dur...

Can Maple represent 2D density plots?...

How to manage the number of colors color...

How to plot the difference between 2 ds...

How to fasten this simulation?...

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