pywafo/wafo/source/rind2007/rind_interface.f

! This is a interface-file for Python
! This file contains a interface to RIND a subroutine
! for computing multivariate normal expectations.
! The file is self contained and should compile without errors on (Fortran90)
! standard Fortran compilers.
!
! The interface was written by
!     Per Andreas Brodtkorb
!     Norwegian Defence Research Establishment
!     P.O. Box 115
!     N-3191 Horten
!     Norway
!     Email: Per.Brodtkorb@ffi.no
!
!
! RIND Computes multivariate normal expectations
!
!  E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix))
!  where
!      "Indicator" = I{ H_lo(i) < X(i) < H_up(i), i=1:N_t+N_d }
!      "Jacobian"  = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is
!      "Jacobian"  = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)|
!      "condition" = Xc=xc(:,ix),  ix=1,...,Nx.
!      X = [Xt; Xd; Xc], a stochastic vector of Multivariate Gaussian
!          variables where Xt,Xd and Xc have the length Nt, Nd and Nc,
!          respectively. (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10)
!
!  CALL: [value,error,terror,inform]=rind(S,m,indI,Blo,Bup,INFIN,xc,
!         Nt,SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT,xCutOff,Nc1c2);
!
!
!    VALUE  = estimated value for the expectation as explained above size 1 x Nx
!    ERROR  = estimated sampling error, with 99% confidence level.   size 1 x Nx
!   TERROR  = estimated truncation error
!   INFORM  = INTEGER, termination status parameter: (not implemented yet)
!            if INFORM = 0, normal completion with ERROR < EPS;
!            if INFORM = 1, completion with ERROR > EPS and MAXPTS
!                           function vaules used; increase MAXPTS to
!                           decrease ERROR;
!            if INFORM = 2, N > 100 or N < 1.
!
!         S = Covariance matrix of X=[Xt;Xd;Xc] size Ntdc x Ntdc (Ntdc=Nt+Nd+Nc)
!         m = the expectation of X=[Xt;Xd;Xc]   size N x 1
!      indI = vector of indices to the different barriers in the
!            indicator function,  length NI, where   NI = Nb+1
!             (NB! restriction  indI(1)=0, indI(NI)=Nt+Nd )
! B_lo,B_up = Lower and upper barriers used to compute the integration
!             limits, Hlo and Hup, respectively. size  Mb x Nb
!    INFIN  = INTEGER, array of integration limits flags:  size 1 x Nb   (in)
!             if INFIN(I) < 0, Ith limits are (-infinity, infinity);
!             if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)];
!             if INFIN(I) = 1, Ith limits are [Hlo(I), infinity);
!             if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)].
!        xc = values to condition on            size Nc x Nx
!        Nt = size of Xt
!      SCIS = Integer defining integration method
!             1 Integrate all by SADAPT for Ndim<9 and by KRBVRC otherwise
!             2 Integrate all by SADAPT by Genz (1992) (Fast)
!             3 Integrate all by KRBVRC by Genz (1993) (Fast)
!             4 Integrate all by KROBOV by Genz (1992) (Fast)
!             5 Integrate all by RCRUDE by Genz (1992)
!   XcScale = REAL to scale the conditinal probability density, i.e.,
!              f_{Xc} = exp(-0.5*Xc*inv(Sxc)*Xc + XcScale)
!    ABSEPS = REAL absolute error tolerance.
!    RELEPS = REAL relative error tolerance.
!    COVEPS = REAL error in cholesky factorization
!    MAXPTS = INTEGER, maximum number of function values allowed. This
!             parameter can be used to limit the time. A sensible
!             strategy is to start with MAXPTS = 1000*N, and then
!             increase MAXPTS if ERROR is too large.
!    MINPTS = INTEGER, minimum number of function values allowed
!    SEED   = INTEGER, seed to the random generator used in the integrations
!    NIT    = INTEGER, maximum number of Xt variables to integrate
!   xCutOff = REAL upper/lower truncation limit of the marginal normal CDF
!    Nc1c2  = INTEGER number of times to use the regression equation to restrict
!             integration area. Nc1c2 = 1,2 is recommended.
!
!
!   If  Mb<Nc+1 then B_lo(Mb+1:Nc+1,:) is assumed to be zero.
!   The relation to the integration limits Hlo and Hup are as follows
!    IF INFIN(j)>=0,
!      IF INFIN(j)~=0,  Hlo(i)=Blo(1,j)+Blo(2:Mb,j).'*xc(1:Mb-1,ix),
!      IF INFIN(j)~=1,  Hup(i)=Bup(1,j)+Bup(2:Mb,j).'*xc(1:Mb-1,ix),
!
!   where i=indI(j-1)+1:indI(j), j=2:NI, ix=1:Nx
!
! This file was successfully compiled for matlab 5.3
! using Compaq Visual Fortran 6.1, and Windows 2000 and windows XP.
! The example here uses Fortran90 source.
! First, you will need to modify your mexopts.bat file.
! To find it, issue the command prefdir(1) from the Matlab command line,
! the directory it answers with will contain your mexopts.bat file.
! Open it for editing. The first section will look like:
!
!rem ********************************************************************
!rem General parameters
!rem ********************************************************************
!set MATLAB=%MATLAB%
!set DF_ROOT=C:\Program Files\Microsoft Visual Studio
!set VCDir=%DF_ROOT%\VC98
!set MSDevDir=%DF_ROOT%\Common\msdev98
!set DFDir=%DF_ROOT%\DF98
!set PATH=%MSDevDir%\bin;%DFDir%\BIN;%VCDir%\BIN;%PATH%
!set INCLUDE=%DFDir%\INCLUDE;%DFDir%\IMSL\INCLUDE;%INCLUDE%
!set LIB=%DFDir%\LIB;%VCDir%\LIB
!
! then you are ready to compile this file at the matlab prompt using the following command:
!
!   mex -O -output mexrind2007 intmodule.f  jacobmod.f rind2007.f mexrind2007.f
!


      subroutine set_constants(method,xcscale,abseps,releps,coveps,
     & maxpts,minpts,nit,xcutoff,Nc1c2, NINT1, xsplit)
      use rindmod, only : setconstants
      use rind71mod, only : setdata
      double precision :: xcscale,abseps,releps,coveps,xcutoff,xsplit
      integer method, maxpts, minpts, nit, Nc1c2, NINT1
Cf2py double precision, optional :: xcscale = 0.0e0
Cf2py double precision, optional :: abseps = 0.01e0
Cf2py double precision, optional :: releps = 0.01e0
Cf2py double precision, optional :: coveps = 1.0e-10
Cf2py double precision, optional :: xcutoff = 5.0e0
Cf2py double precision, optional :: xsplit = 5.0e0

Cf2py integer, optional :: method = 3
Cf2py integer, optional :: minpts = 0
Cf2py integer, optional :: maxpts = 40000
Cf2py integer, optional :: nit = 1000
Cf2py integer, optional :: Nc1c2 = 2
Cf2py integer, optional :: nint1 = 2

! Method>0
      call setconstants(method,xcscale,abseps,releps,coveps,
     &     maxpts,minpts,nit,xcutoff,Nc1c2)
! method==0
      call SETDATA(method,xcscale,abseps,releps,coveps,
     &     nit, xCutOff,NINT1,xsplit)
      return
      end subroutine set_constants
      SUBROUTINE show_constants()
      use rindmod
      print *, 'method=', mMethod
      print *, 'xcscale=', mXcScale
      print *, 'abseps=', mAbsEps
      print *, 'releps=', mRelEps
      print *, 'coveps=', mCovEps
      print *, 'maxpts=', mMaxPts
      print *, 'minpts=', mMinPts
      print *, 'nit=',    mNit
      print *, 'xcutOff=', mXcutOff
      print *, 'Nc1c2=',  mNc1c2
      end subroutine show_constants

      SUBROUTINE rind(VALS,ERR,TERR,Big,Ex,Xc,Nt,INDI,Blo,Bup,
     & INFIN,seed1,Ntdc,Nc,Nx,Ni,Mb,Nb,Nx1)
      USE rindmod
      USE rind71mod, only : rind71
      IMPLICIT NONE
      INTEGER :: Ntd,Nj,K,I
      INTEGER :: seed1
      integer :: Nx,Nx1,Nt, Nc,Ntdc,Ni,Nb,Mb
      DOUBLE PRECISION, dimension(Ntdc,Ntdc) :: BIG
      DOUBLE PRECISION, dimension(Ntdc) :: Ex
      DOUBLE PRECISION, dimension(Nc,Nx1) :: Xc
      DOUBLE PRECISION, dimension(Mb,Nb) :: Blo,Bup
      DOUBLE PRECISION, dimension(Nx) :: VALS, ERR,TERR
      INTEGER, dimension(Ni)  :: IndI
      INTEGER, DIMENSION(Nb) :: INFIN
      INTEGER, ALLOCATABLE  :: seed(:)
      INTEGER               :: seed_size
Cf2py integer, intent(hide), depend(Ex) :: Ntdc = len(Ex)
Cf2py integer, intent(hide), depend(Xc) :: Nc = shape(Xc,0)
Cf2py integer, intent(hide), depend(Xc) :: Nx1 = shape(Xc,1)
Cf2py integer, intent(hide), depend(Xc) :: Nx = max(shape(Xc,1),1)
Cf2py integer, intent(hide), depend(Blo) :: Mb = shape(Blo,0), Nb = shape(Blo,1),
Cf2py integer, intent(hide), depend(Indi) :: Ni = len(Indi)
Cf2py depend(Ntdc)  Big
Cf2py depend(Nb)  INFIN
Cf2py depend(Mb,Nb)  Bup
Cf2py double precision, intent(out), depend(Nx) ::  VALS
Cf2py double precision, intent(out), depend(Nx) ::  ERR
Cf2py double precision, intent(out), depend(Nx) ::  TERR

C     print *, 'Ntdc=', Ntdc,' Nt=',Nt,' Nc=',Nc
C     print *, 'Nx=', Nx, 'Mb=', Mb, ' Nb=', Nb, ' Ni=',Ni
C     Ni = Nb+1
C     Nx = max(Nx1,1)
      if (Ni.EQ.Nb+1) then
      else
         print *, '(ni==nb+1) failed: rind:ni=', Ni, ', nb=',Nb
         return
      endif

      Ntd = Ntdc - Nc;
!    Nd  = Ntd - Nt

      IF (Ntd.EQ.INDI(Ni)) THEN
!     Call the computational subroutine.
        IF (mMethod.gt.0) THEN
          CALL random_seed(SIZE=seed_size)
          ALLOCATE(seed(seed_size))
                               !print *,'rindinterface seed', seed1
          CALL random_seed(GET=seed(1:seed_size)) ! get current state
          seed(1:seed_size)=seed1 ! change seed
          CALL random_seed(PUT=seed(1:seed_size))
          CALL random_seed(GET=seed(1:seed_size)) ! get current state
                                !print *,'rindinterface seed', seed
          DEALLOCATE(seed)
          CALL RINDD(VALS,ERR,TERR,Big,Ex,Xc,Nt,INDI,Blo,Bup,INFIN)
        ELSE
          CALL RIND71(VALS,Big,Ex,Xc,Nt,INDI,Blo,Bup)
          ERR(:) = -1
          TERR(:) = -1
        ENDIF
      ELSE
         print *,'INDI(Ni) must equal Nt+Nd!'
      ENDIF

      RETURN
      END SUBROUTINE rind