fortran from octave

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

fortran from octave

aurorax

Dear all,

I would like to use a function that computes the multivariate normal probability a in the paper "Numerical Computation of Multivariate Normal Probabilities" from Alan Genz.  
This is already implemented in R, however, I want to use the Fortran function that is called from R in Octave and I can't. The closest approach was to try

 
mkoctfile sadmvnt.f

but the sadmvn function does not work in Octave. Can someone help me?

 
Simple example
 
Inputs en R (it works):
 
library("mnormt")
sadmvn(lower=c(-1.57, -0.9295, -0.0708), 
       upper=c(0.63,1.2705,2.1292), 
       mkean = c(0,0,0), 
       varcov = matrix(c(2.7746,-0.0757,0.0358,-0.0757,3.5240,1.3156,0.0358,1.3156,3.7585),3,3) ) 
 
Required Outputs: 0.08349516 , 1.81793e-07
 
 
R package function:
 
sadmvn <- function(lower, upper, mean, varcov,  
                   maxpts=2000*d, abseps=1e-6, releps=0)
{
  if(any(lower > upper)) stop("lower>upper integration limits")
  if(any(lower == upper)) return(0)
  d <- as.integer(if(is.matrix(varcov)) ncol(varcov) else 1)
  varcov <- matrix(varcov, d, d)
  sd  <- sqrt(diag(varcov))
  rho <- cov2cor(varcov)
  lower <- as.double((lower-mean)/sd)
  upper <- as.double((upper-mean)/sd)
  if(d == 1) return(pnorm(upper)-pnorm(lower))
  infin <- rep(2,d)
  infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)
  infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)
  infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)
  infin <- as.integer(infin)
  if(any(infin == -1)) {
    if(all(infin == -1)) return(1)
    k <- which(infin != -1)
    d <- length(k)
    lower <- lower[k]
    upper <- upper[k]
    if(d == 1) return(pnorm(upper) - pnorm(lower))
    rho <- rho[k, k]
    infin <- infin[k]
    if(d == 2) return(biv.nt.prob(0, lower, upper, rep(0,2), rho))
  }
  lower <- replace(lower, lower == -Inf, 0)
  upper <- replace(upper, upper == Inf, 0)
  correl <- as.double(rho[upper.tri(rho, diag=FALSE)])
  maxpts <- as.integer(maxpts)
  abseps <- as.double(abseps)
  releps <- as.double(releps)
  error  <- as.double(0)
  value  <- as.double(0)
  inform <- as.integer(0)
  result <- .Fortran("sadmvn", d, lower, upper, infin, correl, maxpts,
                     abseps, releps, error, value, inform, PACKAGE="mnormt")
  prob <- result[[10]]
  attr(prob,"error")  <- result[[9]]
  attr(prob,"status") <- switch(1+result[[11]], 
                                "normal completion", "accuracy non achieved", "oversize")
  return(prob)
}
 
 
Fortran inputs:
 
d = 3
lower = -0.5658473 -0.2637628 -0.0188373
upper = 0.2270598 0.3605278 0.5665026
infin = 2 2 2
correl = -0.02420905  0.01108602  0.36149196
maxpts = 600
abseps = 1e-06
releps = 0
error = 0
value = 0
inform = 0

I appreciate very much any suggestion.

Best,



---------------

Aurora González Vidal                      Ph.D. student in Data Analytics for Energy Efficiency                    
T. +34 868 88 7866                             Department of Information and Communication Engineering
[hidden email]                  Faculty of Computer Sciences
sae.saiblogs.inf.um.es                        University of Murcia



sadmvnt.f (61K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: fortran from octave

mmuetzel
Am 22. Oktober 2019 um 15:18 Uhr schrieb "AURORA GONZALEZ VIDAL":
> Dear all,
>
> I would like to use a function that computes the multivariate normal probability a in the paper "Numerical Computation of Multivariate Normal Probabilities" from Alan Genz.  
> This is already implemented in R, however, I want to use the Fortran function that is called from R in Octave and I can't. The closest approach was to try
>

> mkoctfile sadmvnt.f
> but the sadmvn function does not work in Octave. Can someone help me?

Please see the Appendix A of the Octave manual. In particular section "A.1.9 Calling External Code from Oct-Files" which discusses how to interface to Fortran code from Octave.

HTH
Markus



Reply | Threaded
Open this post in threaded view
|

Re: fortran from octave

aurorax

Hello Markus,

thank you for your answer. In such appendix, they indicate to change stop with xstop. Appart from that, there is no other indication given that I only have a .f file, do you know if I need to create a C wrapper as it is written in the appendix?

If I use "mkoctfile sadmvnt.f", it creates two files: mkoctfile sadmvnt.o and mkoctfile sadmvnt.oct, but they are almost empty as it can ben seen in the attached screeshots. I also cannot use the function sadmvn as intended.

Thank you,
Aurora

Markus Mützel <[hidden email]> escribió:

Am 22. Oktober 2019 um 15:18 Uhr schrieb "AURORA GONZALEZ VIDAL":

Dear all,

I would like to use a function that computes the multivariate normal probability a in the paper "Numerical Computation of Multivariate Normal Probabilities" from Alan Genz.  
This is already implemented in R, however, I want to use the Fortran function that is called from R in Octave and I can't. The closest approach was to try

 
mkoctfile sadmvnt.f
but the sadmvn function does not work in Octave. Can someone help me?

Please see the Appendix A of the Octave manual. In particular section "A.1.9 Calling External Code from Oct-Files" which discusses how to interface to Fortran code from Octave.

HTHMarkus






---------------

Aurora González Vidal                      Ph.D. student in Data Analytics for Energy Efficiency                    
T. +34 868 88 7866                             Department of Information and Communication Engineering
[hidden email]                  Faculty of Computer Sciences
sae.saiblogs.inf.um.es                        University of Murcia



Captura de pantalla de 2019-10-25 18-09-49.png (23K) Download Attachment
Captura de pantalla de 2019-10-25 18-08-54.png (29K) Download Attachment
sadmvnt.f (61K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: fortran from octave

aurorax
In reply to this post by mmuetzel

Hello Markus and all

thank you for your answer. In such appendix, they indicate to change stop with xstop. Appart from that, there is no other indication given that I only have a .f file, do you know if I need to create a C wrapper as it is written in the appendix?

If I use "mkoctfile sadmvnt.f", it creates two files: mkoctfile sadmvnt.o and mkoctfile sadmvnt.oct, but they are almost empty as it can ben seen in the attached screeshots. I also cannot use the function sadmvn as intended.

Thank you,
Aurora


Markus Mützel <[hidden email]> escribió:

Am 22. Oktober 2019 um 15:18 Uhr schrieb "AURORA GONZALEZ VIDAL":

Dear all,

I would like to use a function that computes the multivariate normal probability a in the paper "Numerical Computation of Multivariate Normal Probabilities" from Alan Genz.  
This is already implemented in R, however, I want to use the Fortran function that is called from R in Octave and I can't. The closest approach was to try

 
mkoctfile sadmvnt.f
but the sadmvn function does not work in Octave. Can someone help me?

Please see the Appendix A of the Octave manual. In particular section "A.1.9 Calling External Code from Oct-Files" which discusses how to interface to Fortran code from Octave.

HTHMarkus






---------------

Aurora González Vidal                      Ph.D. student in Data Analytics for Energy Efficiency                    
T. +34 868 88 7866                             Department of Information and Communication Engineering
[hidden email]                  Faculty of Computer Sciences
sae.saiblogs.inf.um.es                        University of Murcia



Captura de pantalla de 2019-10-25 18-09-49.png (23K) Download Attachment
Captura de pantalla de 2019-10-25 18-08-54.png (29K) Download Attachment
sadmvnt.f (61K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: fortran from octave

mmuetzel
Am 25. Oktober 2019 um 18:14 Uhr schrieb "AURORA GONZALEZ VIDAL":
> Hello Markus and all
>
> thank you for your answer. In such appendix, they indicate to change stop with xstop. Appart from that, there is no other indication given that I only have a .f file, do you know if I need to create a C wrapper as it is written in the appendix?

I have no background in Fortran. But if I were in your shoes, I'd try to follow the example in the manual and call your Fortran function from a C++ wrapper.
If I understand correctly, you *can* use XSTOPX instead of STOP in your Fortran code. But that is optional (although preferred). Just use the F77_XFCN or F77_FCN macros accordingly.

Markus


Reply | Threaded
Open this post in threaded view
|

Re: fortran from octave

Mike Miller-4
In reply to this post by aurorax
On Fri, Oct 25, 2019 at 18:14:32 +0200, AURORA GONZALEZ VIDAL wrote:
> thank you for your answer. In such appendix, they indicate to change stop
> with xstop. Appart from that, there is no other indication given that I only
> have a .f file, do you know if I need to create a C wrapper as it is written
> in the appendix?

Yes, Octave does not have any capability to call a bare Fortran function
directly, you will need to create a C++ wrapper to pass arguments and
return values between Octave and Fortran.

For an example of such a C++ wrapper to call a pretty complicated
Fortran function, take a look at

  http://hg.code.sf.net/p/octave/control/file/ddb727dc0c90/src/sl_ab01od.cc

--
mike



signature.asc (849 bytes) Download Attachment