tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

@acer 

My (rather elderly) machine has 4 cores/8 threads

I know that kernelopts(numcpus) can be set. Like I said in my explanation, if I leave numcpus on 4 for my machine and then run the example, Windows Task Manager shows all eight threads running at nearly 100%. This is the behaviour I have always experienced when using parallelization (that is when I bothered to check). Since all eight threads seem to be running flat out anyway. I always assumed that there would be no point/advantage in trying to set kernelopts(numcpus) to 8.

The most puzzling issue (for me) is that the OP reports strange timing behaviour when changing the "size" of integers in the problem. NB this doesn't change the number of calculations in the problem, arrays are the same size etc.

In my test example

  1. When filling a 4X10,000,000 Array using the large "seed" value (374894756873546859847556), one of the rows in the array calculates integers whose lengths range from 23..93 digits.
  2. When filling a 4X10,000,000 Array using the small "seed" value (700), the same row in the array calculates integer whose lengths range from 2..72 digits.
  3. Timings for either sequential or parallel code didn't noticeably change. For both of the above cases, I get about 55secs for sequential code and 23secs with paralleization. So parallelization gives  ~2.5X speedup

The OP's test case used a 10X  bigger array, ie 4x100,000,000, so the crude "number of calculations" increases by 10X. I would have expected a corresponding 10X increase in execution times. It is also true that increasing the size of the array in the way also increases the length of the "problem" integers. For this test case, OP reports

  1. Using the large "seed" (ie 374894756873546859847556), the increased array size means that returned integer lengths now ought to be 23..103 (compare 23..93 above)
  2. Sequential code runs in 1570secs and parallel code takes 570secs, so about a 2.75X speed-up which makes sense
  3. Using the "small" seed (ie 700), and modifying the calculation somewhat, so that returned integer lengths are smaller, is where OP reports anomalous behaviour. Sequential code runs in  392secs and parallel code is 705secs - ie sequential is nearly 2X faster than parallel. I have trouble believing this
  4. I can't replicate either of these test cases because my RAM usage hits 16GB (which is all I have on this machine), and everything just "hangs". The mserver.exe is at 14.8GB and I have to kill its process tree to regain control.
  5. Further diagnosis of this anomaly is going to need someone with a bigger machine!

@Ronan 

I am not an expert on parallel programming in Maple or anywhere else, so treat all of the following observations with considerable care - some of them may be incorrect!!!

 

Impact of the value of 'n'

Your experience with this (really!!) surprised me.

In my previous answer, I had n=374894756873546859847556 and nelems=10000000. Sequential code ran in ~55secs and parallel code took ~24secs.

I have just repeated this with n=700 and nelems=10000000. Sequential code ran in ~51secs and parallel code took ~24secs.

In other words the execution time for either serial or parallel code does not *seem* to be very dependent on the value of the parameter 'n'. This is actually what I would expect. Maple's integer arithmetic is very quick, and provided one is only doing "simple" arithmetic operations, the size of the integers doesn't matter much. Why you are seeing such a huge diference in execution times between the cases n=374894756873546859847556 and n=700 is a bit of a puzzle for me!

I can only suggest that you take the worksheet I posted previously, record the execution times, then change the two assignments to 'n' to 700, re-execute and record the execution times again - just to see if the discrepancy still exists. I would expect your execution times to be generally  a bit quicker than mine, because I have quite an "elderly" machine

 

Parallelization Overheads

As a general rule one can view parallelization as a three stage process

  1. Set up and intialise the parallel tasks. There is actually a fair amount of "housekeeping" to do here - think stuff like memory control/allocation
  2. Run the parallel tasks
  3. Combine the results from the parallel tasks, and do more "housekeeping" - garbage collection etc

Items (1) and (3) are "overhead" and cost time which is not required in a straightforward sequential caclulation. The "hope" is that item (2) above gives sufficient speed-up to offset the time spent in items (1) and (3). Depending on the calcuation being performed, sometimes combining the results can get a bit "painful".

A side effect of the above is that it is actually rather difficult to wrok out how much time is being taken at each stage of the process - it is possible but gets a bit 'icky'

 

Number of Tasks

This is one which can get really confusing. I set uo the parallelization to use four parallel tasks, mainly becuase I haev a four core/eight thread machine (standard Intel hyperthreading). So far as I am aware, Maple will only parallelize up to the number of physical cores (ie in my case four). You can check how many cores  Maple  'thinks' you have available by executing the command kernelopts(numcpus).

Although I only use four parallel Tasks, if I check the Windows Task manager while the parallel code is executing, all eight threads are close to maxing out. My interpreation of this is that Maple works with "cores" and Windows/Intel enable the (hyper)threading

Althpugh I manually set up the problem to use four tasks, this was just to save me some extra coding. It is possible to use the return value from kernelopts(numcpus) to determine the "optimal" number of tasks and then subdivide the problem accordingly. This just takes a few extra lines of code which I was too lazy to write. It is also not difficult to handle the case where the problem "size" is not a multiple of the number of physical cores, so not all of the tasks end up being exactly the same length. This is just more housekeeping.

@Teep 

the attached, where a new graph is created with vertices appropriately(?) labelled and colored

restart; with(GraphTheory)

I wish to represent the flow across a simple network from a set of sources to a number of  destinations / sinks.

In the example herein, there are 4 sources and 10 possible sinks.

 

The matrix, A, represents the flow configuration as a binary output; 1 for active and 0, otherwise.

For instance, sink numbers 2, 4, 6, 7, 8 and 9 are fed from source 1 (row #1), sink number 10 is supplied from source number 2 (row #2), source number 3 (row #3)is redundant and source number 4 (row #4) supplies sink locations 3 and 5.

 

I would like to assign labels to each source and sink; e.g. Sources labeled A,B,C, ... and sinks 1, 2, 3, ...

This will be very convenient especially for large matrices.

The directed graph could also illustrate the source nodes (lettered) and the corresponding sink nodes (numbered).

 

Does anyone know of a simple routine that can do this?

 

 

A := Matrix(4, 10, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 1, (1, 5) = 0, (1, 6) = 1, (1, 7) = 1, (1, 8) = 1, (1, 9) = 1, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1, (4, 4) = 0, (4, 5) = 1, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0})

The square adjacency matrix, B

B := Matrix(10, A)

Matrix(%id = 18446744074375031254)

(1)

G := Graph(B); DrawGraph(G, style = spring)

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744074375014390), `GRAPHLN/table/1`, 0)

 

 

 

 

Draw the graph of each source and its corresponding sinks.

hSel := proc (spV) options operator, arrow; Graph(select(proc (i) options operator, arrow; i[1] = spV end proc, Edges(G))) end proc

X := DrawGraph(hSel(1), style = spring); DrawGraph(hSel(2), style = spring); DrawGraph(hSel(4), style = spring)

 

 

 

sinks := {seq(`if`(`and`(InDegree(G, j) > 0, OutDegree(G, j) = 0), j, NULL), `in`(j, Vertices(G)))}; sources := {seq(`if`(`and`(OutDegree(G, j) > 0, InDegree(G, j) = 0), j, NULL), `in`(j, Vertices(G)))}; neither := `minus`(convert(Vertices(G), set), `union`(sinks, sources))
nSource := 1nSink := 1
nNeither := 1

for j in Vertices(G) do if member(j, sinks) then HighlightVertex(G, j, stylesheet = [color = red]); C[j] := cat("sink", nSink); nSink := nSink+1 elif member(j, sources) then HighlightVertex(G, j, stylesheet = [color = green]); C[j] := cat("source", nSource); nSource := nSource+1 else HighlightVertex(G, j, stylesheet = [color = yellow]); C[j] := cat("neither", nNeither); nNeither := nNeither+1 end if end do

G2 := RelabelVertices(G, convert(C, list))
DrawGraph(G2, style = spring, size = [800, 800])

 

interface(rtablesize = 10)

 

Download sinkSource2.mw

trying to achieve here.

The best I can get based on the info you provide is shown in the attached, which for reasons I don't understand, won't display on this site :-(

Download anim2.mw

 

In the terminology of Graph Theory

  1. a source
    1. has 1 or more outward-directed edges - ie the OutDegree for a source is >=1
    2. exactly zero inward directed edges - ie the InDegree for a source is exactly 0
  2. a sink
    1. has 1 or more inward-directed edges - ie the InDegree for a sink is >=1
    2. exactly zero outward directed edges - ie the OutDegree for a sink is exactly 0
  3. any vertex which has both an inward-directed edge (InDegree>=1) and an outward-directed edge (OutDegree>=1) is nether a sink nor a soucre

It is relatively trivial to use the InDegree() and OutDegree() commands from the GraphTheory() package with the above logical conditions in order to identify which vertices ares sink, sources or neither for any particular graph.

For your example, I have added an execution group in the attached whihc performs this "classification" of vertices. This clearly demonstrates that vertices 3, 5, 6, 7, 8, 9, 10 are sinks, vertex 1 is a source, and vertices 2, 4 are neither sink nor source.

If we can agree on the definition of "sink", source", and "neither", then we can (probably?) come up with a systematic way to label them. However the simple idea of "letters" for sinks and "numerals" for sources is incomplete, becuase it does not address the issue of vertices which are neither sinks nor sources :-   How should this last category of vertices be labelled??

In the attached, ignoire the the interface(rtablesize=10) command. It only exists because there is a bug on this site which means Maple 2019 worksheets won't display correctly without it: very annoying!!

restart; with(GraphTheory)

I wish to represent the flow across a simple network from a set of sources to a number of  destinations / sinks.

In the example herein, there are 4 sources and 10 possible sinks.

 

The matrix, A, represents the flow configuration as a binary output; 1 for active and 0, otherwise.

For instance, sink numbers 2, 4, 6, 7, 8 and 9 are fed from source 1 (row #1), sink number 10 is supplied from source number 2 (row #2), source number 3 (row #3)is redundant and source number 4 (row #4) supplies sink locations 3 and 5.

 

I would like to assign labels to each source and sink; e.g. Sources labeled A,B,C, ... and sinks 1, 2, 3, ...

This will be very convenient especially for large matrices.

The directed graph could also illustrate the source nodes (lettered) and the corresponding sink nodes (numbered).

 

Does anyone know of a simple routine that can do this?

 

 

A := Matrix(4, 10, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 1, (1, 5) = 0, (1, 6) = 1, (1, 7) = 1, (1, 8) = 1, (1, 9) = 1, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 1, (4, 4) = 0, (4, 5) = 1, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0})

The square adjacency matrix, B

B := Matrix(10, A)

Matrix(%id = 18446744074375031254)

(1)

G := Graph(B); DrawGraph(G, style = spring)

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744074375014390), `GRAPHLN/table/1`, 0)

 

 

 

 

Draw the graph of each source and its corresponding sinks.

hSel := proc (spV) options operator, arrow; Graph(select(proc (i) options operator, arrow; i[1] = spV end proc, Edges(G))) end proc

X := DrawGraph(hSel(1), style = spring); DrawGraph(hSel(2), style = spring); DrawGraph(hSel(4), style = spring)

 

 

 

sinks := seq(`if`(`and`(InDegree(G, j) > 0, OutDegree(G, j) = 0), j, NULL), `in`(j, Vertices(G))); sources := seq(`if`(`and`(OutDegree(G, j) > 0, InDegree(G, j) = 0), j, NULL), `in`(j, Vertices(G))); neither := (`minus`(convert(Vertices(G), set), `union`({sinks}, {sources})))[]

3, 5, 6, 7, 8, 9, 10

 

1

 

2, 4

(2)

interface(rtablesize = 10)

 

Download sinkSource.mw

 

 

 

@JAMET 

Code will need substantial "restructuring" in order to "possibly" work, mainly because objects in the geometry package don't "like" being defined with variables Such "restructuring" may be possible, by using a procedure for the calculation which accepts the unknown 't0' as an argument.

I haven't even attempted this (yet) because a more fundamental problem is that you use the point 'N' in the definition of the triangle 'T', and the point 'N' is nowhere defined in your worksheet

@nm 

I would describe mine as pretty vanilla for 64-bit Windows 7. I pretty much just let the installer "do its thing" Maple 2019 (and all sub-directories) is installed at C:\Program Files\Maple2019.

From memory the only things I did manually with this update installation were

  1. Copy the "DirectSearch"related files from C:\Program Files\Maple 2018\lib to C:\Program Files\Maple 2019\lib
  2. Run the file C:\Program Files\Maple 2019\MapleToolbox2019.0WindowsX64Installer.exe to establish Maple2019 as the symbolic toolbox for my Matlab installation

Now that Physics updates are available for Maple 2019, I have just downloaded/installed the latest update using the "cloud" icon (upper right corner of Maple GUI). This has correctly downloaded/installed the v.333 Physics update and (by default) has placed the associated files in my "user" directory at C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates.

The fact that these updates end up in my "user" directory rather than being installed somewhere under the default MAple installation at C:\Program Files\Maple 2019 has never bothered me much. I assume that it makes life easier if (for any reason) I wanted to revert to whatever Physics version actually shipped with Maple 2019 - just delete/move/rename the stuff in the "user" directory

The problem with superfluous output in your original post (which I could replicate before) has now "gone away".

Only problem I still have, is that in order to get files from Maple 2019, like the attached, to display correctly on this site, I have to include the completely "superfluous" command interfacer(rtablesize=10), somewhere in the worksheet. If this particular  "bug" persists on this site, I can see it becoming really, *really* annoying!
 

pde := a*diff(w(x,y),x) +  b*arctan(lambda*y)*diff(w(x,y),y) =  a*arctan(mu*x)^m+arctan(beta*y)^k:
sol := pdsolve(pde,w(x,y)):

#another example

pde := a*diff(w(x,y),x) +  b*arccot(lambda*y)*diff(w(x,y),y) =  a*arccot(mu*x)^m+arccot(beta*y)^k:
sol:=pdsolve(pde,w(x,y)):

kernelopts(version);

`Maple 2019.0, X86 64 WINDOWS, Mar 9 2019, Build ID 1384062`

(1)

Physics:-Version();

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, March 18, 10:15 hours, version in the MapleCloud: 333, version installed in this computer: 333.`

(2)

interface(rtablesize=10);

[10, 10]

(3)

 


 

Download inttest.mw

 

Re-executing the worksheet ( using !!! with no restart command), the superflous output only appears on the first execution. On second and subsequent executions, it has disappeared!!

If you include a 'restart;' the superfluous output appears every time

@garth 

I posted a response to the other question which you asked here

https://www.mapleprimes.com/questions/226632-How-To-Use-Eulers-Method-With-A-Step

I thought you would be able to use that response to figure out how to do the same thing for the ODE which you post in this query.

Nothing in my original response is "clever". The differential equation posted there (and the one posted here) are pretty trivial and both can be solved exactly. Since both can be solved exactly it is not "obvious" why you would want to solve them numerically. Such numerical solutions will only ever be an approximation to the exact ones.

I freely admit that I had to check Wikipedia to remind myself whick (exactly) is Euler's method for solving ODEs. At my age, remembering the differences between Euler's method, Gear's method, Runge-Kutta methods etc etc gets harder. Luckily for me the Wikipedia page at

https://en.wikipedia.org/wiki/Euler_method

gives a pretty clear description. So all I did was code it (with no real attempt at efficiency/generality etc etc - I just wrote something which "worked")

There are various "pitfalls" associated with numerical solutions, and the obvious one in this case is that you have requested Δt=1.0 which I ( and Carl Love in his response to this thread) think is way too big to get a numerical solution which is even close to the exact solution.


In the attached I have expanded my response to your other post, so that it now deals with the ODE in this question. For the current ODE the attached has a procedure which will return the "numerical" solution using Euler's method for any given stop time and step size. The attached evaluates solutions for step sizes of 0.01, 0.1 and 1.0. I would characterise these solutions as "good", "marginal" and "awful", respectively.

  restart;
  with(plots):
################################################
# Solve the simple RC problem exactly, simply
# because it can be done!
#
  V:= z-> 2*cos(3*z):
  R:= 4:
  C:= 0.5:
  ode:= diff( V__c(t),t)= (V(t)-V__c(t))/(R*C):
  sol:= dsolve( {ode, V__c(0)=2}):
#
# Plot this exact solution
#
  p1:= plot( rhs(sol),
             t=0..10,
             color=red,
             title="Exact solution",
             titlefont=[times, bold, 24],
             thickness=5,
             labelfont= [times, bold, 20],
             size=[1000, 500]
           );
####################################################
# Solve the system numerically using Euler's method
# with a step size of 0.1
#
  V__c:=Array(0..100):
  V__c[0]:=2:
  t__step:=0.1:
  for i from 1 by 1 to 100 do
      V__c[i]:=V__c[i-1]+t__step*(V((i-1)*t__step)-V__c[i-1])/(R*C)
  od:
#
# Plot this numerical solution
#
  p2:=pointplot( [seq( [ t__step*i, V__c[i]], i=0..100)],
                 symbol=solidcircle,
                 symbolsize=10,
                 color=blue,
                 title="Numerical solution",
                 titlefont=[times, bold, 20],
                 size=[1000, 500]
               );
#####################################################
# Display the "exact" and the "numerical" solution
# on the same graph - pretty good agreement!
#
  display( [p1, p2],
           title="Exact and Numerical Solutions",
           titlefont=[times, bold, 20],
           size=[1000, 500]
         )

 

 

 

###############################################################
#
# Op's other differential equation
#
# Solve exactly and plot the result, just becasue it is trivial
# to do and usefule for comparisons with numerical methods
#
  ode:=diff(y(t),t)=(3-y(t))*(y(t)+1):
  sol:=dsolve( [ode, y(0)=4] ):
  p1:= plot( rhs(sol),
             t=0..5,
             color=red,
             title="Exact solution",
             titlefont=[times, bold, 24],
             thickness=5,
             labelfont= [times, bold, 20],
             size=[1000, 500]
           );
#
# Now define a *quick and dirty" procedure which will numerically
# solve this ode for a given step size and stop time
#
  getODEsol:= proc( tstop, tstep )
                    uses plots:
                    local yn:= Array(0..round(tstop/tstep)),
                          i:
                    yn[0]:= 4:
                    for i from 1 by 1 to round(tstop/tstep) do
                        yn[i]:=yn[i-1]+tstep*(3-yn[i-1])*(yn[i-1]+1);
                    od:
                    return pointplot
                           ( [ seq
                               ( [ tstep*i, yn[i]],
                                   i=0..round(tstop/tstep)
                               )
                             ],
                             symbol=solidcircle,
                             symbolsize=10,
                             color=blue
                           );
              end proc:
#
# Plot the numerical solutions from the above
# procedure for step sizes of 0.01, 0.1 and 1.0.
#
# Include the exact solution obtained earlier on
# each plot just for comparison purposes
#
  stopTime:= 5.0:
  stepSize:= [0.01, 0.1, 1.0]:
  for k from 1 by 1 to numelems(stepSize) do
      display
      ( [ p1, getODEsol(stopTime, stepSize[k])],
        title= typeset
              ( "Exact and Numerical Solutions\n",
                " ( Step Size=", stepSize[k], " )"),
        titlefont=[times, bold, 20],
        size=[1000, 500]
      );
  end do;

 

 

 

 

 

Download odeEuler.mw

 

@Josolumoh 

I have better things to do than fix basic syntax because you cannot be bothered to read the manual. A (non-exhaustive) list of problems you should address contains

  1. Notice that in your definition of the function 'B', the Maple output contains the names 'simplify' and 'unapply'. In other words, these have not been interpreted, nor used, as commands. Did you ever look at the output of this function definition?? Why does this happen?? Because you have decided to use the combination of Maple's 2-D input and code indentation. I wouldn't recommend this combination unless you really really know what you are doing (and you don't!). Maple's 2D input parser can place all sorts of interestiing interpretations on the "whitespace" which you have  inserted purely for indentation purposes. In particular never put any "whitespace" (or anything can be interpreted as "whitespace") between a Maple command and its associated opening parenthesis: so for example, use 'simplify(', and not 'simplify ('
  2. Your final for loop contains references to the function 'SUBS', which in turn refers to the name 'val_b[n]'. The only other reference to the latter in this worksheet is in a line which has been commented out. So 'val_b[n]' will never evaluate to anything "useful", hence 'SUBS' will never evaluate to anything "useful"...and so on through the whole of this for loop

@Josolumoh 

returns a sequence - the clue is in the name.

If you want the output returned as a list, then enclose it in '[]' like

[ seq( X__B(X__0, SX__1[k], SX__2[k]), k=1..N) ];

If you want the output returned as a vector, then enclose it in '<>' like

<seq( X__B(X__0, SX__1[k], SX__2[k]), k=1..N) >;

@Josolumoh 

to my previous post where it says

# The 'local k;' construct in the function definition
# is necessary to suppress a (superfluous) warning
# message in Maple 2019. If using Maple 2018 or earlier
# the function definition should be written as
#
# X__B:= (x, y, z)-> exp(add(beta[k]*[x,y,z][k+1], k = 0 .. 2));
#
  X__B:= (x, y, z)-> local k;
                     exp
                     ( add
                       ( beta[k]*[x,y,z][k+1],
                         k = 0 .. 2
                       )
                     );

Now explain clearly which part you don't understand

*Seems* to be doing exactl what I woulld expect, it is difficult for me to work out what you think is going wrong. As a *guess* from your description, you seem to be confusing two fundamentally different operations, which can be summarised as.

FIRST SCENARIO

  1. Define a couple of random variables.
  2. Sample these random variables;
  3. Define a function of these random variables;
  4. Sample this function. The values obtained here  will bear no relation to the samples obtained at step (2) above, since "new" samples of the random variables defined at step(1) above will be generated

SECOND SCENARIO

  1. Define a couple of random variables.
  2. Sample these random variables;
  3. Define a function
  4. Apply this function to the samples obtained at step(2) above. In this operation, the function defined at step (3) above is not *technically* a random variable and cannot be meaningfull "sampled". All it does is use the sample values which were generated at step (2) above

It is possible to define the function (ie step (3) in either of the above "scenarios"), so that either of the above "scenarios" can be handled. This is done in the attached.

(Be aware the attachment was developed in Maple 2019, which is very new and has a few quirks. For example, ignore the rtablesize command at the end: this is just to overcome a problem about displaying worksheet contents on this site)

  restart;
  with(Statistics):
  #randomize();
  B:= unapply
      ( simplify
        ( binomial(x + r - 1, x)*(b/2/(1 + d*x + b/2))^r*((d*x + 1)/(1 + d*x + b/2))^x/(d*x + 1) ),
        [r, b, d]
      );
  beta:= Array( 0..2,
                [0.25, 0.6, -0.2]
              );
  N:= 5;

proc (r, b, d) options operator, arrow; binomial(x+r-1, x)*(b/(2*d*x+b+2))^r*2^x*((d*x+1)/(2*d*x+b+2))^x/(d*x+1) end proc

 

Array(%id = 18446744074393090278)

 

5

(1)

#
# Define some random variables and get a few samples
#
  X__0:= 1;
  X__1:= RandomVariable(Bernoulli(1/2));
  X__2:= RandomVariable(Normal(0, 1));
  SX__1:= Sample(X__1, N);
  SX__2:= Sample(X__2, N);

X__0 := 1

 

X__1 := _R

 

X__2 := _R0

 

Vector[row](5, {(1) = 1.0, (2) = 1.0, (3) = 0., (4) = 0., (5) = 0.})

 

Vector[row](%id = 18446744074393081726)

(2)

#
# Define a "new" random variable based on previous
# random variables. Observe that every time one asks
# for "samples" from this random variable, one will
# get a "new" list of samples - as expected
#
  X__B:= exp
         ( add
           ( beta[k]*X__||k,
             k = 0 .. 2
           )
         );
  Sample_b:= Sample(X__B, N);
  Sample_b:= Sample(X__B, N);
  Sample_b:= Sample(X__B, N);

X__B := exp(.25+.6*_R-.2*_R0)

 

Vector[row](5, {(1) = .9726008729365625, (2) = 2.150871828664164, (3) = 1.4046239626237194, (4) = 1.4658659638712146, (5) = 2.07645686621087})

 

Vector[row](5, {(1) = .9267272256435518, (2) = .9815610071033524, (3) = 2.234555636247529, (4) = 2.652756505020624, (5) = 2.506187556048575})

 

Vector[row](%id = 18446744074393071734)

(3)

#
# Re-define the new variable in a slightly
# different way (just for convenience really).
# This is not technically a random variable
# because its behaviour is governed by the
# arguments passed to it.
#
# If one passes random variables, then it behaves
# as an rv, and one can generate samples. Each call
# to Sample() will generate a *new* list of samples
# exactly as before.
#
# The 'local k;' construct in the function definition
# is necessary to suppress a (superfluous) warning
# message in Maple 2019. If using Maple 2018 or earlier
# the function definition should be written as
#
# X__B:= (x, y, z)-> exp(add(beta[k]*[x,y,z][k+1], k = 0 .. 2));
#
  X__B:= (x, y, z)-> local k;
                     exp
                     ( add
                       ( beta[k]*[x,y,z][k+1],
                         k = 0 .. 2
                       )
                     );
  Sample_b:= Sample( X__B( X__0, X__1, X__2), N);
  Sample_b:= Sample( X__B( X__0, X__1, X__2), N);
  Sample_b:= Sample( X__B( X__0, X__1, X__2), N);

X__B := proc (x, y, z) options operator, arrow; exp(add(beta[k]*[x, y, z][k+1], k = 0 .. 2)) end proc

 

Vector[row](5, {(1) = 2.138742062366558, (2) = 2.4975703902556288, (3) = 1.1955404063578525, (4) = 1.5578336221095606, (5) = 2.231405041664074})

 

Vector[row](5, {(1) = 2.4545973012034645, (2) = 3.2895866688172197, (3) = 2.6735374323358188, (4) = 2.541300779153873, (5) = 2.768105648812239})

 

Vector[row](%id = 18446744074393063662)

(4)

#
# On the other hand one can now evaluate this new
# function using samples which have been previously
# obtained eg the vectors SX__1 and SX__2 defined
# above. In this case the function itself does *not*
# define a new "random variable". Using the Sample()
# command in this case would rsult in an error because
# one cannot sample something which is not an rv
#
# Obviously repeated "application" will give the same
# answer every time.
#
  seq( X__B(X__0, SX__1[k], SX__2[k]), k=1..N);
  seq( X__B(X__0, SX__1[k], SX__2[k]), k=1..N);

HFloat(2.2461753368121835), HFloat(2.386372343185135), HFloat(1.305439309043902), HFloat(1.3279333023039865), HFloat(1.2430844367475704)

 

HFloat(2.2461753368121835), HFloat(2.386372343185135), HFloat(1.305439309043902), HFloat(1.3279333023039865), HFloat(1.2430844367475704)

(5)

interface(rtablesize=10):

 

Download ranVars.mw

@JAMET 

what "improvement" do you want?

@wlferguson19 

my own fibonacci generator, I would probably go one of the methods shown in the attached

#
# Recursive Fibonacci generator
#
  myFib:= proc(n::integer)
               option remember;
               if   n=1
               then return 1
               elif n=0
               then return 0
               else return myFib(n-1)+myFib(n-2):
               fi:
          end proc:

  seq(myFib(j), j=0..20);

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765

(1)

#
# Explicit Fibonacci generator
#
  f:=n->floor(((1 + sqrt(5))/2)^n/sqrt(5) + 1/2):

  seq( f(j), j=0..20);

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765

(2)

 

 

Download fibon.mw

 

First 70 71 72 73 74 75 76 Last Page 72 of 207