Marsaglia’s MWC (multiply with carry) Random Number Generator for Object Pascal

Originally Posted on December 14, 2019 by hsauro 

In the past, I’ve tended to use the Mersenne Twister for generating pseudo-random numbers. This is a widely used generator because it has a very long period of 2^{19937} – 1. The one issue is that it’s not the fastest generator and this is particularly an issue when drawing random numbers for implementing the Gillespie direct method. A much faster variant is the Multiply-with-carry (MWC) random number generator. This uses much simpler arithmetic which means it’s faster. The algorithm can be described in three lines:

  

  

  

The symbol >> means shift the bits to the right by the given number and << is similar but moves the bits to the left. In object pascal, these operations are represented by shr and shl respectively. The generator needs two values to be initialized, m_z and m_w I've recently implemented the Gillespie method in three languages, C, Go and Object Pascal. Object Pascal has good support for the Mersenne Twister but I had to create my own MWC generator. The code below shows the object pascal implementation:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

// ------------------------------------------------------------ // George Marsaglia's MWC (multiply with carry) generator. // This code translated from C# to Object Pascal and obtained from // https://www.codeproject.com/Articles/25172/Simple-Random-Number-Generation   unit uRandMWC;   interface   Uses System.Types;   // Call this before doing anything procedure initMWC;   // If you want to change the seed call this procedure setSeed (seed : cardinal);   // Call this to get a random number between 0 and 2^32 function  getUint() : cardinal;   // Call this to get a uniformly distributed random number between 0 and 1.0 function getUdouble : double;   implementation   // In Delphi 32 and 64 bit, cardinal is an unsigned 32-bit integer   var   m_z, m_w : cardinal;   procedure setSeed (seed : cardinal); begin   m_z := seed;   m_w := 678934; end;   procedure initMWC; begin   m_z := 467567;   m_w := 125681; end;   function getUint() : cardinal; begin   m_z := 36969 * (m_z and 65535) + (m_z shr 16);   m_w := 18000 * (m_w and 65535) + (m_w shr 16);   result := (m_z shl 16) + m_w; end;   function getUdouble : double; var u : dword; begin   u := getUint();   // The magic number below is 1/(2^32 + 2).   // The result is strictly between 0 and 1.   result := double (u + 1.0) * 2.328306435454494e-10; end;   end.

What about some tests? There are random generator test systems such as TestU01 and dieharder but they are written to run on Linux and I didn't have time to set something up on windows plus the usage docs are not that great. Instead I did three simple sanity tests. The simplest is to check that the mean value for a sample of random numbers is 0.5. Another test is to make sure that, say out of 10,000 random numbers, 10% are between 0 and 0.1, etc. Ideally there should be statistical tests done to check that these values are within certain tolerances. The third test is a visual test where I generate 32768 random numbers and map them on to a bitmap. Each random number is 32 bits long and each bit represents a pixel in the bitmap. 32768 numbers means I generate a bitmap 1024 by 1024 pixels. Thus 32 numbers occupy each row of the bitmap. I wrote a simple Delphi console application that calls the random number generator and creates the bitmap. The main part of this code is shown below. For a comparison I also implemented the NOT to be used randu generator and also the Delphi's built in generator. The approach I take isn't very sophisticated and I am sure it could be made shorter, but for each random number I generate, I convert it into a string of 1's and 0's. I then work my way through the string setting pixels in a bitmap. I do this 32 times for each row in the bitmap.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

program RandomGen;   {$APPTYPE CONSOLE} {$R *.res}   uses   System.SysUtils,   System.Types,   Vcl.Graphics,   uRandMWC in 'uRandMWC.pas';   var r1 : double;       ri : cardinal;     i, j, k : integer; bmp : TBitmap;     sb : string;       iCol : integer;     state : cardinal;  count : integer;   // Get a single bit as position Bit from cardinal aValue function getBit(const aValue: cardinal; const Bit: Byte): Boolean; begin   Result := (aValue and (1 shl Bit)) <> 0; end;   // Converts a cardinal (32-bits) into a string of ones and zeros function getBinaryString (aValue : cardinal) : string; var i : integer; begin   result := '';   for i := 0 to 32 - 1 do       if getBit (aValue, i) then          result := result + '1'       else         result := result + '0' end;   function randu : cardinal; begin    state := (65539 * state) mod 2147483648;    result := state; end;   begin   state := 25; // starting value for radu   writeln ('wait.....');   initMWC();     bmp := TBitmap.Create;   bmp.PixelFormat := pf1bit;   bmp.Width := 1024;   bmp.Height := 1024;   count := 0;   for i := 0 to bmp.Height - 1 do       begin       iCol := 0;       for j := 0 to bmp.Width div 32 - 1 do           begin           // Choose a random number generator           //ri := randu;           ri := getUint();           //ri := random (4294967296-1);           sb := getBinaryString (ri);           for k := 0 to 32 - 1 do               begin               if sb[k+1] = '0' then                  begin                  bmp.Canvas.Pixels[iCol,i] := clWhite;                  end               else                  begin                  bmp.Canvas.Pixels[iCol,i] := clBlack;                  end;               inc (iCol);               end;           end;       inc (count, bmp.Width div 32);       end;     bmp.SaveToFile('image.bmp');   writeln ('count = ', count);   writeln ('Press any key to coninue');   readln; end.

The results of running this code is below. First randu to show you how bad it is, note the patterns and lines in the bitmap:

The next figure shows the results from Delphi’s builtin random number generator, which isn’t so great either:

Finally, the MWC algorithm, which you can see from the image the distribution of pixels appears to show no pattern:

A simple bistable system

Originally Posted on November 27, 2019 by hsauro

Here is a simple bistable system I sometimes use in class. I like it because it’s easy to explain plus the addition of one extra reaction turns it into a relaxation oscillator.

import tellurium as te
import random
import pylab as plt
import numpy as np

# ---- BISTABLE SYSTEM -----

r = te.loada ('''
     $Xo -> S1;  0.5 + Vmax*S1^n/(15 + S1^n);
      S1 -> ;    k1*S1;

     Xo = 1; S1 = 1;
     n = 4; Vmax = 10; k1 = 2;
''')

plt.xlabel ('Time')
plt.ylabel ('Concentration')
for S1 in np.arange (1.45, 1.55, 0.01):
    r.S1 = S1 
    m = r.simulate (0, 5, 100)
    plt.plot (m['time'], m['[S1]'])
   
r.steadyStateSolver.allow_presimulation = False
for i in range (10):
    r.S1 = random.random()*5
    r.steadyState()
    print ('S1 = ', r.S1, 'eigenvalue = ', r.getFullEigenValues()[0])

If you run the above code you’ll get:

S1 = 5.145227 eigenvalue = -1.84055
S1 = 5.145227 eigenvalue = -1.84055
S1 = 0.251329 eigenvalue = -1.95768
S1 = 5.145227 eigenvalue = -1.84055
S1 = 1.492258 eigenvalue = 3.004970
S1 = 5.145227 eigenvalue = -1.84055
S1 = 1.492258 eigenvalue = 3.004970
S1 = 5.145227 eigenvalue = -1.84055
S1 = 5.145227 eigenvalue = -1.84055
S1 = 0.251329 eigenvalue = -1.957684

The above output is what I call the poor-mans way of finding multiple steady states. I scan across the inital conditions computing the steady state in each case and the corresponding eigenvalue (You can also just assign random values to the initial value for S1, which is better if you don’t know what starting conditions to use). As you can see it found three steady states, two stable (neg eignvalue) and one unstable (1.492, pos eigenvalue).

The plot above shows time course trajectories for about 10 different starting points for S1. The two steady states are clearly evident.

Out of interest here is the addition of the extra reaction at the front to turn it into a relaxation oscillator. You have to make a slight modification to the second reaction’s rate law so that it becomes a function of its reactant (which it wasn’t before)

      $Xo -> So;  k0*Xo;
      So -> S1;   k1*So + Vmax*So*S1^n/(15 + S1^n);
      S1 -> ;     k2*S1;

      Xo = 1; S1 = 1; k0 = 0.07; k1 = 0.04
      n = 4; Vmax = 6; k2 = 0.1;

Simulate Regular Doses using Antimony/SBML in the Tellurium Python Tool

Originally Posted on November 21, 2019 by hsauro

Someone recently asked how to create a model where a molecular species received a regular dose at a fixed interval of 12 hours.

They wanted something where a species A is being degraded but at regular intevals is increased by a fixed amount. See the plot below.

Lucian Smith and I came up with two possible solutions which are shown below. The first solution from Lucian is to set up a variable called trigger that we define using the differential equation trigger’ = 1 (note the apostrophe). The ode ensures that trigger increases linearly during the simulation. We then set up an event that triggers when trigger is greater than 12 hours at which point trigger is set back to zero and we add the dose of value 3 to species A. trigger then starts to integrate again from zero and the cycle repeats.

import tellurium as te

r = te.loada("""
      A ->; k1*A;
      k1 = 0.1; A = 10
      
      trigger = 0
      trigger' = 1

      at (trigger>12): A = A+3, trigger=0
""")

m = r.simulate (0, 100, 100, ['time', 'A'])
r.plot()

The second solution is less elegant but increases the trigger value when an event occurs. In the code below when time exceeds the current trigger value we add the dose then increase trigger by 12 hours so that can it can be triggered again.

m = r.simulate (0, 100, 100, ['time', 'A'])
r.plot()

r = te.loada("""
      A ->; k1*A;
      k1 = 0.1; A = 10
      
      trigger = 12
      at (time > trigger): A = A+3, trigger = trigger + 12
""")

m = r.simulate (0, 100, 100, ['time', 'A'])
r.plot()

Generating log ranges in Python

Originally Posted on September 27, 2019 by hsauro

This is a very short blog on how to generate logarithmic ranges in python:

import numpy as np

# Brute force
x = [0.1, 1, 10, 100, 1000]

# Classic explicit loop
x = 0.1
for i in range (5):
    print (x)
    x = x*10
    
# List comprehension, more efficient
x = [10**i for i in range(-1,3)]
print (x)

# Lazy-evaluation (saves memory for very large lists)
for x in (10**i for i in range(-1,3)):
    print (x)
    
# Using the built-in function in numpy
x = np.logspace (-2,2,num=5)
print (x)

One of my best friends passed away today

Originally Posted on September 28, 2019 by hsauro

It is with enormous and deep sorrow that one of my best friends in the world passed away today as a result of Hemangiosarcoma. This is where the spleen becomes cancerous which eventually results in fatal hemorrhaging. There is virtually no treatment for this. She would likely pass away on her own accord within 24 hours. We didn’t know whether she was in pain or not, she rarely complained about anything, but we suspected she was. She looked absolutely miserable. We decided to make the incredibly difficult decision to put her to sleep instead of letting her suffer for another 24 hours. An incredibly sad day for all the family.

We got her 12.5 years ago as a puppy, full of life. We named her Rhody after the blossoming rhododendrons in our garden. She was a black Labrador, the sweetest dog you could imagine. She would wait for me outside to come home from work, and if she was in the house she’d listen for the sound of the bus that brought me home. When I entered the house she was almost always the first to greet me. She had an uncanny skill of measuring time. Her afternoon meal was at 4pm and somehow she knew when that time came. She would come to me or my wife looking at us in anticipation for the late lunch, often exactly at 4pm. The same thing would happen at 8pm for her evening snack.

After coming back from a walk with my wife she would rush (and I mean rush) downstairs to the basement to say hello to me. Like all Labradors, she was obsessed with fetching balls and sticks and loved getting into water. She also knew how to relax and in the evening would lie on her back, legs up. I believe we gave her a happy life, she certainly gave much joy and love in return. She would be ecstatic when we went on trips to cabins we’d rent. We always picked a place where there was a river because jumping into rivers was by far her absolute favorite pastime. It’s difficult to be sure, but she seemed to say thank you for the things we gave her. This was especially the case when I filled her dish with food. She would first walk back from her dish, circle me with a wagging tail while looking at me and return to her food. She also had a thing where she would smile at us, probably mimicking what we did. When my wife and kids were away on a trip and I was on my own I would let her upstairs to sleep in with me on her own bed. She would be so excited to do that. Some say Labradors aren’t the sharpest of dogs, but I can tell you she was a genius when it came to food. She had a wide range of words she knew, these included: rhody, ball, stick, get, fetch, sleep, bed, stay, wait, sit, lie, paw (i.e give me your paw), go pee, let’s go (meaning run), leash, walk, good girl, bad girl, dish, treat, squirrel, cat, rat, water, toy, car (she loved riding in the car), post (which actually meant let’s check out the fig tree on the road by our neighbor's house next to the postbox), and what’s this (i.e check this out). She had a thing about figs. We also have a fig tree in the garden and as the figs developed she would inspect them almost daily to see if they were ripe. Once they were ripe, they would be gone. At Christmas, we always had a wrapped present for her (a bit chewy bone).

This morning at 7.15am (28 Sept), we heard a loud falling sound. Rhody had collapsed, her legs had given away. We rushed downstairs, she was lying flat out at the bottom of the stairs where she usually waited for us to come downstairs in the morning. She couldn’t move and was clearly in a lot of distress. We took her to the vet, her ailment was diagnosed and we faced the abyss of having to lose her. It was an incredibly traumatic moment. You do not want to go through what we went through. It was unbearable. She was still conscious, her faculties were intact but her body was failing her and there was nothing that could be done about it. Imagine having to make the decision to let her go. For my wife and myself, this was the hardest and saddest decision we have ever made.

Rhody, we will miss you so deeply. May you rest in peace, our best friend and most loyal companion. Farewell.