Reaction Diffusion in python

Hi, was looking @laurent_delrieu { https://www.grasshopper3d.com/forum/topics/2d-reaction-diffusion-to-3d?id=2985220%3ATopic%3A1341676&page=1#comments } work on reaction diffusion on C#.

I was wondering if this can be done on python or not?

any lead would be appreciated, Thank you.

For sure it could be replicated. There is no reason it can’t. Is your problem in c# or python part ?

Hey Laurent, I’m much more comfortable in python, but not able to replicate your code in python. Its quite intense coding done in C#.

thanks!

Hello
It is far from intense coding, I speak of the C#. For sure there are not a lot of comments but you only need to make some standard math with arrays and lists. I am not sure it is hard to duplicate in Python.
It is just 100 lines of code. It is not a lot and you are lucky Python is 0 based for index like C#. If someone is fluent in Python (I don’t need to be)

  private void RunScript(int _n, List<double> _A, List<double> _B, double _Da, List<double> _Db, double _f, double _k, double _dt, int _n_step, int _n_jump, ref object Aout, ref object Bout, ref object Astack, ref object Bstack)
  {

    //Number of points/vertex in the mesh
    int n_points;
    n_points = _n * _n;

    double[] A = new double[n_points];
    double[] B = new double[n_points];
    double[] Aprime = new double[n_points];
    double[] Bprime = new double[n_points];
    List<double> A_stacked = new List<double>();
    List<double> B_stacked = new List<double>();

    A = _A.ToArray();
    B = _B.ToArray();

    double dxA;
    double dxB;
    double AB2;
    int i_point;

    for (int i_p = 0; i_p < n_points; i_p++)
    {
      A_stacked.Add(1.0);
      B_stacked.Add(1.0);
    }

    for (int i = 0; i < _n_step; i++)
    {
      for (int j = 0; j < _n_jump; j++)
      {
        for (int i_x = 0; i_x < _n; i_x++)
        {
          for (int i_y = 0; i_y < _n; i_y++)
          {
            //Laplacian calculation using convolution Matrix
            dxA = 0.0;
            dxB = 0.0;

            dxA = dxA + 0.05 * A[modulo(i_x - 1 + (i_y - 1) * _n, n_points)];
            dxB = dxB + 0.05 * B[modulo(i_x - 1 + (i_y - 1) * _n, n_points)];
            dxA = dxA + 0.20 * A[modulo(i_x - 0 + (i_y - 1) * _n, n_points)];
            dxB = dxB + 0.20 * B[modulo(i_x - 0 + (i_y - 1) * _n, n_points)];
            dxA = dxA + 0.05 * A[modulo(i_x + 1 + (i_y - 1) * _n, n_points)];
            dxB = dxB + 0.05 * B[modulo(i_x + 1 + (i_y - 1) * _n, n_points)];

            dxA = dxA + 0.2 * A[modulo(i_x - 1 + (i_y ) *_n, n_points)];
            dxB = dxB + 0.2 * B[modulo(i_x - 1 + (i_y ) *_n, n_points)];
            dxA = dxA - 1.0 * A[modulo(i_x - 0 + (i_y ) *_n, n_points)];
            dxB = dxB - 1.0 * B[modulo(i_x - 0 + (i_y ) *_n, n_points)];
            dxA = dxA + 0.2 * A[modulo(i_x + 1 + (i_y ) *_n, n_points)];
            dxB = dxB + 0.2 * B[modulo(i_x + 1 + (i_y ) *_n, n_points)];

            dxA = dxA + 0.05 * A[modulo(i_x - 1 + (i_y + 1) * _n, n_points)];
            dxB = dxB + 0.05 * B[modulo(i_x - 1 + (i_y + 1) * _n, n_points)];
            dxA = dxA + 0.20 * A[modulo(i_x - 0 + (i_y + 1) * _n, n_points)];
            dxB = dxB + 0.20 * B[modulo(i_x - 0 + (i_y + 1) * _n, n_points)];
            dxA = dxA + 0.05 * A[modulo(i_x + 1 + (i_y + 1) * _n, n_points)];
            dxB = dxB + 0.05 * B[modulo(i_x + 1 + (i_y + 1) * _n, n_points)];
            //Reaction-diffusion equation

            i_point = modulo(i_x + i_y * _n, n_points);
            AB2 = A[i_point] * B[i_point] * B[i_point];
            Aprime[i_point] = A[i_point] + (_Da * dxA - AB2 + _f * (1.0 - A[i_point])) * _dt;
            Bprime[i_point] = B[i_point] + (_Db[i_point] * dxB + AB2 - (_k + _f) * B[i_point]) * _dt;
          }
        }
        for (int i_p = 0; i_p < n_points; i_p++)
        {
          A[i_p] = Math.Min(Math.Max(Aprime[i_p], 0.0), 1.0);
          B[i_p] = Math.Min(Math.Max(Bprime[i_p], 0.0), 1.0);
        }
        if (j == 0)
        {
          for (int i_p = 0; i_p < n_points; i_p++)
          {
            A_stacked.Add(A[i_p]);
            B_stacked.Add(B[i_p]);
          }
        }
      }
    }


    for (int i_p = 0; i_p < n_points; i_p++)
    {
      A_stacked.Add(1.0);
      B_stacked.Add(1.0);
    }
    Aout = A;
    Bout = B;
    Astack = A_stacked;
    Bstack = B_stacked;


  }

  // <Custom additional code> 

  public  int modulo(int n, int size) {
    int retour = n;
    if (n < 0)
    {retour = n + size;}
    if (n >= size)
    {retour = n % size;}
    return retour;
  }
1 Like

Hi,

Here’s a rough Python translation of Laurent’s C# code:

def modulo(n, size):
    if n < 0:
        return n + size
    elif n >= size:
        return n % size
    else:
        return n
    

n_points = _n * _n

A = list(_A)
B = list(_B)
Aprime = [0.0 for _ in range(n_points)]
Bprime = [0.0 for _ in range(n_points)]
A_stacked = []
B_stacked = []

for i_p in range(n_points):
    A_stacked.append(1.0)
    B_stacked.append(1.0)
    
for i in xrange(_n_step):
    for j in xrange(_n_jump):
        for i_x in xrange(_n):
            for i_y in xrange(_n):
                
                # Lacplacian calculation using convolution matrix
                dxA = 0.0
                dxB = 0.0
            
                dxA = dxA + 0.05 * A[modulo(i_x - 1 + (i_y - 1) * _n, n_points)]
                dxB = dxB + 0.05 * B[modulo(i_x - 1 + (i_y - 1) * _n, n_points)]
                dxA = dxA + 0.20 * A[modulo(i_x - 0 + (i_y - 1) * _n, n_points)]
                dxB = dxB + 0.20 * B[modulo(i_x - 0 + (i_y - 1) * _n, n_points)]
                dxA = dxA + 0.05 * A[modulo(i_x + 1 + (i_y - 1) * _n, n_points)]
                dxB = dxB + 0.05 * B[modulo(i_x + 1 + (i_y - 1) * _n, n_points)]

                dxA = dxA + 0.2 * A[modulo(i_x - 1 + (i_y ) * _n, n_points)]
                dxB = dxB + 0.2 * B[modulo(i_x - 1 + (i_y ) * _n, n_points)]
                dxA = dxA - 1.0 * A[modulo(i_x - 0 + (i_y ) * _n, n_points)]
                dxB = dxB - 1.0 * B[modulo(i_x - 0 + (i_y ) * _n, n_points)]
                dxA = dxA + 0.2 * A[modulo(i_x + 1 + (i_y ) * _n, n_points)]
                dxB = dxB + 0.2 * B[modulo(i_x + 1 + (i_y ) * _n, n_points)]

                dxA = dxA + 0.05 * A[modulo(i_x - 1 + (i_y + 1) * _n, n_points)]
                dxB = dxB + 0.05 * B[modulo(i_x - 1 + (i_y + 1) * _n, n_points)]
                dxA = dxA + 0.20 * A[modulo(i_x - 0 + (i_y + 1) * _n, n_points)]
                dxB = dxB + 0.20 * B[modulo(i_x - 0 + (i_y + 1) * _n, n_points)]
                dxA = dxA + 0.05 * A[modulo(i_x + 1 + (i_y + 1) * _n, n_points)]
                dxB = dxB + 0.05 * B[modulo(i_x + 1 + (i_y + 1) * _n, n_points)]
             
                # Reaction-diffusion equation
                i_point = modulo(i_x + i_y * _n, n_points)
                AB2 = A[i_point] * B[i_point] * B[i_point]
                Aprime[i_point] = A[i_point] + (_Da * dxA - AB2 + _f * (1.0 - A[i_point])) * _dt
                Bprime[i_point] = B[i_point] + (_Db[i_point] * dxB + AB2 - (_k + _f) * B[i_point]) * _dt
    
        for i_p in xrange(n_points):
            A[i_p] = min(max(Aprime[i_p], 0.0), 1.0);
            B[i_p] = min(max(Bprime[i_p], 0.0), 1.0);
        
        if j == 0:
            for i_p in xrange(n_points):
                A_stacked.append(A[i_p])
                B_stacked.append(B[i_p])

for i_p in xrange(n_points):
    A_stacked.append(1.0)
    B_stacked.append(1.0)

Aout = A
Bout = B
Astack = A_stacked
Bstack = B_stacked

It was pretty straightforward to port, - even-though what happens is pretty cryptic to me - and it works, however the performance is inferior to the C# version.
This is due to C# being far more forgiving, meaning that you don’t have to write super optimised code (no offence Laurent) in order for things to run rather smoothly. You see, C# is a semi-compiled language, which makes it by nature faster than Python.
Furthermore, Python is performance-wise pretty “badly” implemented in Grasshopper, and Rhino probably isn’t meant to deal with large amounts of information and calculations in the first place.
I came to this conclusion, because in comparison to other CG apps (i.e. Maya, Houdini) that I also write Python scripts for, the performance of Rhino pales. This doesn’t mean that Rhino is an inferior application, quite the contrary, it’s probably just not what it is meant to do. And its price-performance-ratio is much better than these other applications!
Grasshopper just helped opening a window into the dark, bottomless pit that is parametric, avant-garde design. And in the innovation and novelty driven culture of today, it mostly can’t keep up performance-wise with what it has participated in kickstarting in the first place.

Of course you can always write superior, optimised (Python) code, and things will run more smoothly (buy a faster computer, or learn C#, or even better C++*). Doing that is up to you now! :wink:

Thanks for sharing the script, Laurent!

*C++ is fully compiled, which means it’s even faster than C# and it isn’t the baby of Micro$oft :black_flag:

2 Likes

Are they running CPython rather than IronPython? That could explain the difference. CPython is implemented in C with hundreds of contributors optimising the fundamental portions of it…

There are a lot of modulo function calls in that script. That would be an obvious first candidate for optimisation. If n is always >= -size then you should be able to delete the modulo function and use the % operator instead. In any event you don’t need the third else statement

1 Like

If you need speed see this link. I didn’t test it because at the moment I don’t need more speed.

1 Like

Here is the script without all those modulo calls - it is still slow but not so slow:


n_points = _n * _n

A = list(_A)
B = list(_B)
Aprime = [0.0 for _ in range(n_points)]
Bprime = [0.0 for _ in range(n_points)]
A_stacked = []
B_stacked = []

for i_p in range(n_points):
    A_stacked.append(1.0)
    B_stacked.append(1.0)
    
for i in xrange(_n_step):
    for j in xrange(_n_jump):
        for i_x in xrange(_n):
            for i_y in xrange(_n):
                
                # Lacplacian calculation using convolution matrix
                dxA = 0.0
                dxB = 0.0
            
                dxA = dxA + 0.05 * A[((i_x - 1 + (i_y - 1) * _n) % n_points)]
                dxB = dxB + 0.05 * B[((i_x - 1 + (i_y - 1) * _n) % n_points)]
                dxA = dxA + 0.20 * A[((i_x - 0 + (i_y - 1) * _n) % n_points)]
                dxB = dxB + 0.20 * B[((i_x - 0 + (i_y - 1) * _n) % n_points)]
                dxA = dxA + 0.05 * A[((i_x + 1 + (i_y - 1) * _n) % n_points)]
                dxB = dxB + 0.05 * B[((i_x + 1 + (i_y - 1) * _n) % n_points)]

                dxA = dxA + 0.2 * A[((i_x - 1 + (i_y ) * _n) % n_points)]
                dxB = dxB + 0.2 * B[((i_x - 1 + (i_y ) * _n) % n_points)]
                dxA = dxA - 1.0 * A[((i_x - 0 + (i_y ) * _n) % n_points)]
                dxB = dxB - 1.0 * B[((i_x - 0 + (i_y ) * _n) % n_points)]
                dxA = dxA + 0.2 * A[((i_x + 1 + (i_y ) * _n) % n_points)]
                dxB = dxB + 0.2 * B[((i_x + 1 + (i_y ) * _n) % n_points)]

                dxA = dxA + 0.05 * A[((i_x - 1 + (i_y + 1) * _n) % n_points)]
                dxB = dxB + 0.05 * B[((i_x - 1 + (i_y + 1) * _n) % n_points)]
                dxA = dxA + 0.20 * A[((i_x - 0 + (i_y + 1) * _n) % n_points)]
                dxB = dxB + 0.20 * B[((i_x - 0 + (i_y + 1) * _n) % n_points)]
                dxA = dxA + 0.05 * A[((i_x + 1 + (i_y + 1) * _n) % n_points)]
                dxB = dxB + 0.05 * B[((i_x + 1 + (i_y + 1) * _n) % n_points)]
             
                # Reaction-diffusion equation
                i_point = (i_x + i_y * _n)  % n_points
                AB2 = A[i_point] * B[i_point] * B[i_point]
                Aprime[i_point] = A[i_point] + (_Da * dxA - AB2 + _f * (1.0 - A[i_point])) * _dt
                Bprime[i_point] = B[i_point] + (_Db[i_point] * dxB + AB2 - (_k + _f) * B[i_point]) * _dt
    
        for i_p in xrange(n_points):
            A[i_p] = min(max(Aprime[i_p], 0.0), 1.0);
            B[i_p] = min(max(Bprime[i_p], 0.0), 1.0);
        
        if j == 0:
            for i_p in xrange(n_points):
                A_stacked.append(A[i_p])
                B_stacked.append(B[i_p])

for i_p in xrange(n_points):
    A_stacked.append(1.0)
    B_stacked.append(1.0)

Aout = A
Bout = B
Astack = A_stacked
Bstack = B_stacked
4 Likes

One ‘out of the box thinking’ technique i use when my IDE (RunRev - grand child of HyperCard) is too slow is to write a wrapper for a function. That is to offset the hard work to a well compiled exe or native command line. Might not apply here, but there’s always a way.

(until i learn rhinoscript, this is runrev/hypertalk script as i write it pretty much - should be human interpretable into any language)

function RD inputdata
// create input
 put "\\mysrv\share" into inputdatafilepath
 // save inputs for the app - might be easier than passing arguments (powershell for example is a horrible for this)
SavetoFile(inputdata,inputdatafilepath)
// run the external compiled code out of rhino
  _Run "RDonSteroids.exe" & space & filepath
 // get results
 put "\\mysrv\share\myoutputdata.txt" into outfilepath
 put ReadFile(outfilepath) into ouputdata 
 return myoutputdata
end RD

Then you can compile in pure C if you want. :slight_smile:
I learned about Phyton’s performance, thanks.

Hi Laurent. I state that your works are always fantastic. I was taking a look at your script. But one thing is not clear to me. What is the purpose of the n_jump for loop. I thought just n_step loop was enough. If you could make this clear to me, I would be grateful.
thanks

Hello
I thing this was just used to make 3d reaction diffusion.
https://www.grasshopper3d.com/m/discussion?id=2985220%3ATopic%3A1341676
So you are making n steps slice of reaction diffusion
Each slice separated by n step steps.
So no use for classical use…

OK! I got it . So to make the shape grow you created a layer every tot iteration. A last question, What is this chunk of code for?

        for i_p in xrange(n_points):
            A[i_p] = min(max(Aprime[i_p], 0.0), 1.0);
            B[i_p] = min(max(Bprime[i_p], 0.0), 1.0);

Many thanks!

Just to limit the output between 0 and 1

Thanks!