Search code examples
pythonfortranf2py

Fortran subroutine returning None in python c.f. working in R


Im trying to implement some probability work that utilises the bivariate normal cdf. In R I can use the pbivnorm package (line 90 specifically), e.g. (with some made up numbers)

library(pbivnorm)
.Fortran("PBIVNORM", as.double(0), c(0,0), as.double(-0.1), as.double(-0.2), as.integer(c(0,0)), as.double(0), as.integer(1), PACKAGE="pbivnorm")

[[1]]
[1] 0.193613

[[2]]
[1] 0 0

[[3]]
[1] -0.1

[[4]]
[1] -0.2

[[5]]
[1] 0 0

[[6]]
[1] 0

[[7]]
[1] 1

I've taken the underlying Fortran code from the package here and compiled it using F2py (named the compiled .so file fortran) so I can e.g. do:

import fortran
fortran.pbivnorm(float(0), [float(0), float(0)], float(-0.1), float(-0.2), [int(0), int(0)], float(0),int(1))

however this returns None. I know that the subroutine technically is changing the values supplied to it so I also supplied various named constants but they remain unchanged when queried.

Also feel like maybe could be to do with fussiness over input types?

Any help much appreciated!


Solution

  • You are right with an assumption that there is a mismatch with the types. Default python types (lists) are not ok for calling this function. Also the variable PROB is used as the output for this function, so you can not just use a literal for it. By looking at the original fortran function signature, your arguments and code should look like this:

    import fortran
    import numpy as np
    
    
    PROB = np.array([0], order="F", dtype=np.double)
    LOWER = np.array([0, 0], order="F", dtype=np.double)
    UPPERA = np.array([-0.1], order="F", dtype=np.double)
    UPPERB = np.array([-0.2], order="F", dtype=np.double)
    INFIN = np.array([0, 0], order="F", dtype=np.int32)
    CORREL = np.array([0], order="F", dtype=np.double)
    LENGTH = np.int32(1)
    
    fortran.pbivnorm(PROB, LOWER, UPPERA, UPPERB, INFIN, CORREL, LENGTH)
    
    print(PROB, LOWER, UPPERA, UPPERB, INFIN, CORREL, LENGTH)