Python embedding via cpp

src/initialization/MOM_state_initialization.F90

@@ -97,6 +97,7 @@ module MOM_state_initialization
 use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd
 use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field
 use MOM_oda_incupd, only: calc_oda_increments, output_oda_incupd_inc
+use MOM_python_embedding, only: runpyF_1_3d
 
 implicit none ; private
 
@@ -643,7 +644,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                     PF, oda_incupd_CSp, restart_CS, Time)
   endif
 
-
+  if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Testing Python embedding ...")
+  Call runpyF_1_3d('test.py', 'py_test1', 'tv%T',  tv%T)
+  if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Stop here to debug!!"); stop
 end subroutine MOM_initialize_state
 
 !> Reads the layer thicknesses or interface heights from a file.

src/initialization/runpyC.C

@@ -0,0 +1,154 @@
+//Inspired from https://github.com/wangsl/python-embedding
+//#include "runpyC.h"
+#include <stdio.h>
+#include <iostream>
+#include <cassert>
+
+#include <Python.h>
+#include <numpy/ndarrayobject.h>
+
+// https://docs.python.org/3/extending/embedding.html
+// https://github.com/dusty-nv/jetson-utils/issues/44
+// export PYTHONPATH=.:$PYTHONPATH
+
+#define FORT(x) x##_
+
+static PyObject *py_module = 0;
+static PyObject *py_function = 0;
+
+inline int init_numpy()
+{
+  if(!PyArray_API) import_array();
+  return PyArray_API ? 1 : 0;
+}
+
+static void runpyC_1_1d(const char *py_script_, const char *py_function_,
+		             double *x, const int &n1)
+{
+  if(!Py_IsInitialized()) Py_Initialize();
+
+  assert(Py_IsInitialized());
+  assert(init_numpy());
+  printf("runpyC_1_1d: n1,%d \n", n1);  
+  printf("From runpyC_1_1d: Trying to call %s in %s \n",py_function_,py_script_);
+  if(!py_module) {
+    PyObject *py_script = PyUnicode_DecodeFSDefault(py_script_);
+    assert(py_script);
+    py_module = PyImport_Import(py_script);
+//Interesting: The following assertion fails if 
+//             the py_script has a python bug even outside the function that is called!
+//             the py_script has import torch
+    assert(py_module);
+    Py_DECREF(py_script);
+  }
+  if(!py_function) {
+    assert(py_module);
+    py_function = PyObject_GetAttrString(py_module, py_function_);
+    assert(py_function);// This assertation fails if py_function is not found in py_script
+    assert(PyCallable_Check(py_function));
+  }
+  //Simle tests
+  //assert(PyObject_CallFunctionObjArgs(py_function, NULL)); //This worked. Printed from the function.
+
+  const npy_intp dim_x [] = { n1 };
+  PyObject *x_py = PyArray_SimpleNewFromData(1, dim_x, NPY_DOUBLE, x);
+  PyObject *pValue;
+  assert(pValue = PyObject_CallFunctionObjArgs(py_function, x_py, NULL));
+  if (pValue != NULL) {
+       printf("Result of call: %f\n", PyFloat_AsDouble(pValue));
+       Py_DECREF(pValue);
+  }
+  Py_DECREF(x_py); x_py = 0;
+
+  // PyRun_SimpleString("import sys; sys.stdout.flush()");
+
+  std::cout.flush();
+}
+
+static void runpyC_1_3d(const char *py_script_, const char *py_function_,
+			const char *varname_,
+			double *x, const int &n1, const int &n2, const int &n3 )
+{
+  if(!Py_IsInitialized()) Py_Initialize();
+
+  assert(Py_IsInitialized());
+  assert(init_numpy());
+  printf("runpyC_1_3d: varname %s , dims %d,%d,%d \n", varname_, n1,n2,n3);  
+  printf("From runpyC_1_3d: Trying to call %s in %s.py \n",py_function_,py_script_);
+  if(!py_module) {
+    PyObject *py_script = PyUnicode_DecodeFSDefault(py_script_);
+    assert(py_script);
+    py_module = PyImport_Import(py_script);
+//Interesting: The following assertion fails if 
+//             the py_script has a python bug even outside the function that is called!
+//             the py_script has import torch
+    assert(py_module && " This assertation fails if the py_script has a python bug..");
+    Py_DECREF(py_script);
+  }
+  if(!py_function) {
+    assert(py_module);
+    py_function = PyObject_GetAttrString(py_module, py_function_);
+    assert(py_function && " This assertation fails if the function is not found in the script."); //This assertation fails if py_function is not found in py_script
+    assert(PyCallable_Check(py_function));
+  }
+  //Simle tests
+  //assert(PyObject_CallFunctionObjArgs(py_function, NULL)); //This worked. Printed from the function.
+
+  const npy_intp dim_x [] = { n1,n2,n3 };
+  PyObject *x_py = PyArray_SimpleNewFromData(3, dim_x, NPY_DOUBLE, x);
+  PyObject *pValue;
+  assert(pValue = PyObject_CallFunctionObjArgs(py_function, x_py, NULL));
+  if (pValue != NULL) {
+       printf("Result of call: %f\n", PyFloat_AsDouble(pValue));
+       Py_DECREF(pValue);
+  }
+  Py_DECREF(x_py); x_py = 0;
+
+  // PyRun_SimpleString("import sys; sys.stdout.flush()");
+
+  std::cout.flush();
+}
+
+static void py_finalize()
+{
+  if(py_function) { Py_DECREF(py_function); py_function = 0; }
+  if(py_module) { Py_DECREF(py_module); py_module = 0; }
+  if(Py_IsInitialized()) assert(!Py_FinalizeEx());
+  std::cout.flush();
+}
+
+// Fortran interface: NOTE: the name of function should be al lowercase!
+extern "C" void FORT(runpyc_1array3d) (const char *py_script, const int &len_py_script,
+		  const char *py_function, const int &len_py_function,
+		  const char *varname, const int &len_varname,
+		  double *x, const int &n1, const int &n2, const int &n3 )
+
+{
+  char *py_script_ = new char [len_py_script+1];
+  assert(py_script);
+  memcpy(py_script_, py_script, len_py_script*sizeof(char));
+  py_script_[len_py_script] = '\0';
+
+  char *py_function_ = new char [len_py_function+1];
+  assert(py_function);
+  memcpy(py_function_, py_function, len_py_function*sizeof(char));
+  py_function_[len_py_function] = '\0';
+
+  char *varname_ = new char [len_varname+1];
+  assert(varname);
+  memcpy(varname_, varname, len_varname*sizeof(char));
+  varname_[len_varname] = '\0';
+
+  runpyC_1_3d(py_script_, py_function_, varname_, x, n1,n2,n3);
+
+  if(py_script_) { delete [] py_script_; py_script_ = 0; }
+  if(py_function_) { delete [] py_function_; py_function_ = 0; }
+  if(varname_) { delete [] varname_; varname_ = 0; }
+}
+
+// Fortran interface: PyFinalize 
+extern "C" void FORT(pyfinalize)()
+{
+  py_finalize();
+}
+

src/initialization/runpyF.F90

@@ -0,0 +1,26 @@
+!Inspired by https://github.com/wangsl/python-embedding
+!
+module MOM_python_embedding
+implicit none
+public runPyF_1_3d
+contains
+
+subroutine runpyF_1_3d(pyScript, pyFunction, varname, x)
+  character(len=*),                   intent(in)    :: pyScript, pyFunction,varname
+  real, dimension(:,:,:),             intent(inout) :: x 
+  integer :: n1,n2,n3, len_pyScript, len_pyFunction , len_varname   
+
+  len_pyScript = Len_Trim(pyScript)
+  if(pyScript(len_pyScript-2:len_pyScript) .eq. '.py') len_pyScript = len_pyScript-3   
+  len_pyFunction = Len_Trim(pyFunction)   
+  len_varname  = Len_Trim(varname)   
+
+  n1 = size(x,1)
+  n2 = size(x,2)       
+  n3 = size(x,3)       
+  call runpyc_1array3d(pyScript, len_pyScript, pyFunction, len_pyFunction,&
+                       varname, len_varname,  x, n1,n2,n3)
+
+end subroutine runpyF_1_3d
+
+end module MOM_python_embedding

src/initialization/test.py

@@ -0,0 +1,73 @@
+#!/bin/env python
+#Inspired from https://github.com/wangsl/python-embedding
+import sys
+import numpy
+import matplotlib.pyplot as plt
+#A test function with one argument
+def py_test1(x) :
+    print(' **** From py_test1 ****')
+    print("py_test1: Shape of the input array x : ", x.shape)
+    print("py_test1: x[4,4,0] x[4,4,-1] : ", x[4,4,0],x[4,4,-1])
+
+    plt.plot(x[4,4,:]);plt.show();
+
+    a=numpy.max(x)
+    print("py_test1: Returning the numpy.max(x): ", a)
+
+    sys.stdout.flush()
+    return a
+
+#A test function with two arguments
+def py_test2(x,y) :
+    print(' **** From py_test2 ****',x,y)
+    sys.stdout.flush()
+    a=sum(x)
+    b=sum(y)
+    return a,b
+
+#A test function with no arguments
+def py_test0():
+    print(' **** From py_test0 ****')
+    sys.stdout.flush()
+    return 
+
+from math import sin
+#import torch
+#import numpy  #This gives all kinds of library errors at runtime
+
+def my_test_torch(x, y) :
+    print(' From Python test: {}'.format(x.size))
+    for i in range(x.size) :
+        x[i] += 1.0
+
+    a = torch.from_numpy(x)
+    print(a)
+
+    b = torch.from_numpy(x.astype(numpy.float32))
+    print(b)
+
+    if torch.cuda.is_available() :
+        b_dev = b.cuda()
+        print(b_dev)
+
+    for i in range(y.size) :
+        y[i] = 2*numpy.float64(b[i]) + 0.1
+
+    print(y)
+
+    sys.stdout.flush()
+
+
+if __name__ == '__main__' :
+    import numpy as np
+
+    x = np.arange(10, dtype=numpy.float64)
+    print(x)
+
+    y = np.arange(10, dtype=numpy.float64)
+    print(y)
+    my_test(x, y)
+
+    print(x)
+
+