`random()` call within custom kernels not working

[VeckoTheGecko]

Parcels version: 2.3.1 OS: Windows 10

I've been trying to specify random behaviour within custom parcels kernels, however including any call to parcels.rng.random or random.random causes errors both when using JITParticle and ScipyParticle for the particle set. From prior parcels code I have seen having kernels with random behaviour does seem to work in an older version, however I'm not sure what I would need to do to get this to work in this version of parcels.

The following is a trimmed example of my code, with the file random_not_working.nc available here (for 7 days after posting).

from parcels import ParticleSet, JITParticle, AdvectionRK4, ErrorCode, FieldSet, ScipyParticle
from parcels.rng import random
from datetime import timedelta
from pathlib import Path
import numpy as np

variables = {"U": "U", "V": "V"}
dimensions = {"lon": "x_mesh", "lat": "y_mesh", "time": "time"}
field_set = FieldSet.from_netcdf(
    "random_not_working.nc",
    variables,
    dimensions,
    mesh="flat",
    interp_method="cgrid_velocity",
    gridindexingtype="mitgcm",
    allow_time_extrapolation=True 
    )
array = np.array([25_000, 25_000])
pset = ParticleSet(fieldset=field_set, pclass=JITParticle, lon=array, lat=array)

def kernel_random(particle, fieldset, time):
    """
    Kernel with some random behaviour
    """
    rand = random()

combined_kernel = AdvectionRK4 + pset.Kernel(kernel_random)
pset.execute(combined_kernel,
            runtime=timedelta(days=24),  # runtime controls the interval of the plots
            dt=timedelta(minutes=5),
            )  # the recovery kernel

Where I get the following error messages for JITParticle

---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\compilation\codecompiler.py:229, in CCompiler._create_compile_process_(self, cmd, src, log)
    228 try:
--> 229     subprocess.check_call(cmd, stdout=logfile, stderr=logfile)
    230 except OSError:

File ~\anaconda3\envs\parcels_hons\lib\subprocess.py:369, in check_call(*popenargs, **kwargs)
    368         cmd = popenargs[0]
--> 369     raise CalledProcessError(retcode, cmd)
    370 return 0

CalledProcessError: Command '['gcc', '-Wall', '-fPIC', '-std=gnu11', '-IC:\\Users\\nickh\\anaconda3\\envs\\parcels_hons\\lib\\site-packages\\parcels\\include', '-I.', '-g', '-O3', '-DDOUBLE_COORD_VARIABLES', '-m64', '-o', 'C:\\Users\\nickh\\AppData\\Local\\Temp\\parcels-tmp\\lib8750aac502f3e2c53dd5ae53b52131ed_0.dll', 'C:\\Users\\nickh\\AppData\\Local\\Temp\\parcels-tmp\\8750aac502f3e2c53dd5ae53b52131ed_0.c', '-shared', '-lm', '-m64']' returned non-zero exit status 1.

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Input In [9], in <cell line: 22>()
     19     rand = random()
     21 combined_kernel = AdvectionRK4 + pset.Kernel(kernel_random)
---> 22 pset.execute(combined_kernel,
     23             runtime=timedelta(days=24),  # runtime controls the interval of the plots
     24             dt=timedelta(minutes=5),
     25             )

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\particleset\baseparticleset.py:340, in BaseParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
    338         self.kernel.remove_lib()
    339         cppargs = ['-DDOUBLE_COORD_VARIABLES'] if self.collection.lonlatdepth_dtype else None
--> 340         self.kernel.compile(compiler=GNUCompiler(cppargs=cppargs, incdirs=[path.join(get_package_dir(), 'include'), "."]))
    341         self.kernel.load_lib()
    343 # Set up the interaction kernel(s) if not set and given.

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\kernel\basekernel.py:253, in BaseKernel.compile(self, compiler)
    251         if self.src_file is not None:
    252             all_files_array.append(self.src_file)
--> 253         compiler.compile(self.src_file, self.lib_file, self.log_file)
    254 if len(all_files_array) > 0:
    255     logger.info("Compiled %s ==> %s" % (self.name, self.lib_file))

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\compilation\codecompiler.py:290, in GNUCompiler_SS.compile(self, src, obj, log)
    287 lib_pathdir = os.path.dirname(obj)
    288 obj = os.path.join(lib_pathdir, lib_pathfile)
--> 290 super(GNUCompiler_SS, self).compile(src, obj, log)

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\compilation\codecompiler.py:268, in CCompiler_SS.compile(self, src, obj, log)
    266 with open(log, 'w') as logfile:
    267     logfile.write("Compiling: %s\n" % " ".join(cc))
--> 268 self._create_compile_process_(cc, src, log)

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\compilation\codecompiler.py:242, in CCompiler._create_compile_process_(self, cmd, src, log)
    235                 with open(log, 'r') as logfile2:
    236                     err = """Error during compilation:
    237 Compilation command: %s
    238 Source/Destination file: %s
    239 Log file: %s
    240 
    241 Log output: %s""" % (" ".join(cmd), src, logfile.name, logfile2.read())
--> 242                 raise RuntimeError(err)
    243         return True

RuntimeError: Error during compilation:
Compilation command: gcc -Wall -fPIC -std=gnu11 -IC:\Users\nickh\anaconda3\envs\parcels_hons\lib\site-packages\parcels\include -I. -g -O3 -DDOUBLE_COORD_VARIABLES -m64 -o C:\Users\nickh\AppData\Local\Temp\parcels-tmp\lib8750aac502f3e2c53dd5ae53b52131ed_0.dll C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c -shared -lm -m64
Source/Destination file: C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c
Log file: C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.log

Log output: C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c:1:0: warning: -fPIC ignored for target (all code is position independent)
 #include "parcels.h"
 ^
C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c: In function 'AdvectionRK4kernel_random':
C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c:127:10: warning: implicit declaration of function 'random' [-Wimplicit-function-declaration]
   rand = random();
          ^
C:\Users\nickh\AppData\Local\Temp\parcels-tmp\8750aac502f3e2c53dd5ae53b52131ed_0.c:83:82: warning: variable 'rand' set but not used [-Wunused-but-set-variable]
   type_coord u1, v1, lon1, lat1, u2, v2, lon2, lat2, u3, v3, lon3, lat3, u4, v4, rand;
                                                                                  ^
C:\Users\nickh\AppData\Local\Temp\cc9K0GYe.o: In function `AdvectionRK4kernel_random':
C:/Users/nickh/AppData/Local/Temp/parcels-tmp/8750aac502f3e2c53dd5ae53b52131ed_0.c:127: undefined reference to `random'
collect2.exe: error: ld returned 1 exit status

and ScipyParticle

---------------------------------------------------------------------------
KernelError                               Traceback (most recent call last)
Input In [12], in <cell line: 22>()
     19     rand = random()
     21 combined_kernel = AdvectionRK4 + pset.Kernel(kernel_random)
---> 22 pset.execute(combined_kernel,
     23             runtime=timedelta(days=24),  # runtime controls the interval of the plots
     24             dt=timedelta(minutes=5),
     25             )

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\particleset\baseparticleset.py:452, in BaseParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, moviedt, recovery, output_file, movie_background_field, verbose_progress, postIterationCallbacks, callbackdt)
    450 # If we don't perform interaction, only execute the normal kernel efficiently.
    451 if self.interaction_kernel is None:
--> 452     self.kernel.execute(self, endtime=next_time, dt=dt, recovery=recovery, output_file=output_file,
    453                         execute_once=execute_once)
    454 # Interaction: interleave the interaction and non-interaction kernel for each time step.
    455 # E.g. Inter -> Normal -> Inter -> Normal if endtime-time == 2*dt
    456 else:
    457     cur_time = time

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\kernel\kernelsoa.py:227, in KernelSOA.execute(self, pset, endtime, dt, recovery, output_file, execute_once)
    225 recovery_kernel = recovery_map[p.state]
    226 p.set_state(StateCode.Success)
--> 227 recovery_kernel(p, self.fieldset, p.time)
    228 if p.isComputed():
    229     p.reset_state()

File ~\anaconda3\envs\parcels_hons\lib\site-packages\parcels\tools\statuscodes.py:121, in recovery_kernel_error(particle, fieldset, time)
    119 """Default error kernel that throws exception"""
    120 msg = "Error: %s" % particle.exception if particle.exception else None
--> 121 raise KernelError(particle, fieldset=fieldset, msg=msg)

KernelError: 0
Particle P[6](lon=25000.000000, lat=25000.000000, depth=0.000000, time=0.000000)
Time: 2022-01-10T18:00:00.000000000,    timestep dt: 300.000000
Error: 'module' object is not callable

[erikvansebille]

Parcels uses its own ParcelsRandom library, see e.g. the Diffusion Kernels at https://github.com/OceanParcels/parcels/blob/dd3ba5b3f99b106b275feb87389a63b4b03f4216/parcels/application_kernels/advectiondiffusion.py#L29-L30

The full list of supported random functions is

• ParcelsRandom.seed()
• ParcelsRandom.random()
• ParcelsRandom.uniform(low, high),
• ParcelsRandom.randint(low, high)
• ParcelsRandom.normalvariate(loc, scale)
• ParcelsRandom.expovariate(lamb)
• ParcelsRandom.vonmisesvariate(mu, kappa)

The code is at https://github.com/OceanParcels/parcels/blob/752efb1f57c06fa9e9cc6e1e4ab871c35ffa3137/parcels/rng.py#L137-L195

[VeckoTheGecko]

This is really weird. Replacing

from parcels.rng import random
...
rand = random()

with

import parcels.rng as ParcelsRandom
...
rand = ParcelsRandom.random()

fixes all the issues. From my understanding, there shouldn't be anything different between the two pieces of code (as random() is imported from parcels.rng). Is there something I'm missing @erikvansebille ? Is this error reproducible? Could this be due to how it compiles to C?

I posted in the Python Discord #help-donut channel about this and I think they're equally stumped 😂

[erikvansebille]

Yes indeed, it's intentional that you have to use the ParcelsRandom methods; this has to do with the code converted to C, which can't shadow the built-in random function. Long technical story...

A warning indeed is a good idea, just like we issue a warning when someone uses numpy in a kernel. https://github.com/OceanParcels/parcels/blob/34eb1cfeb7fc92963ce32f75788c3cdef2466b10/parcels/compilation/codegenerator.py#L318-L321