! 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=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