C     ROB SALGADO
C
C
C
C     CALCULATE THE VOLUME IN A UNIT n-BALL
C
C     why does this method work?
C     The idea came to me as follows:
C         One can map a family of concentric squares within the unit square
C         into a family of concentric circles within the unit disk, and vice
C         versa.
C         Given a random point (xr,yr) in the unit square, map it to a point on
C         a circle of radius max(xr,yr).  Although for the purpose of this
C         calculation, it is not necessary to specify exactly where on the
C         circle, I chose to map it to a point at obtained by "projection"
C         from the line y=yr (assuming yr>xr) to the circle.
C         So, now I see that a thin shell of points around this circle makes
C         a contribution to the area of the disk as a function of the
C         radius of this determined circle.  Thus, I should weight the
C         points randomly picked in the unit square by a function of the
C         radius.  I looked for an inverse map...
C         Although I have not explicitly derived it analytically for all n,
C         when I apply it to the problem of calculating the enclosed volume
C         it seems to work reasonably well.
C
C     P.S.
C         The program seems to bomb when I use fixed matrices of size
C         greater than (16384,~20)...  I'm not sure why this is so.
C         But, the program is simple enough in structure to either
C         use several matrices or use "two rows or so" of random numbers
C         to make up the components of my position vectors.
C         Since the program seems to run relatively fast, I would like to
C         try to calculate the volume of, say, a 100-ball.
C
C
 




      program mppsphere
 
      integer, parameter :: N=16384, maxdim=20
      real, parameter :: Two31=2147483648.0
      integer IR(N*maxdim),  SR(N*maxdim)
      real v(N,maxdim)
 
      integer dim
      real, parameter :: pi=3.1415926535
 
CMPF  mpl MPF_RANDOM, MPF_SRANDOM
      integer mptm                     ! timer

      integer samp,r
      real logv,truevol
 
      real, array(N) :: t, l, hit, oldhit
      real, array(1000) :: monterad,oldmonterad
      real mean, oldmean, sumdevsq, oldsumdevsq
 
 
      open (unit=3, file='d-sphere')
      write(3,*) 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
      print*, 'Importance d-Sphere Monte-Carlo [ROB SALGADO]'
 
C-----SEED THE RANDOM NUMBER GENERATOR
      SR= [1:maxdim*N]
      call MPF_SRANDOM(SR)
 
      print*, 'Dimension dim=?'
      read*, dim
C-----------------------------COMPUTE TRUE D-VOLUME
C     Use logarithms because partial factors for high dimensions
C     overflows the variables.
      if (2*int(dim/2).eq.dim) goto 8
C-------DIM IS ODD
      logv=(dim+1.0)/2.0*log(2.0)+(dim-1.0)/2.0*log(pi)
      goto 9
C-------DIM IS EVEN
 8    logv=(dim/2.0)*log(2.0*pi)
C-------COMPUTE PARITY-FACTORIAL n(n-2)(n-4)...
 9    do 5 i=dim,1,-2
           logv=logv-log(i*1.0)
 5    continue
      truevol=exp(logv)
 
C-------------------------------------------------|
      print*, 'Number of samples per run?'
      read*, samp
 10   print*, 'Number of runs?'
      read*, r
      if (r.gt.1000) goto 10
 
C------------------------------------------------------(r) RUNS
      write(3,*) '               Trial     Radius'
D     print*,    '               Trial     Radius'
      do 200 j=1,r
         hit=0
         oldhit=0
C---------------------------------(N*samp) SAMPLES

      call mpTimerStart() 
      do 100 i=1,samp
         l=0.0
         t=0.0
 
      call MPF_RANDOM(IR)
      v=reshape( SHAPE=[N,maxdim], SOURCE=IR )/Two31
      v(:,(dim+1):maxdim)=0.0
                                                     ! shape into dim-vectors
      t=maxval(v,2)                                  ! get maximum coordinate
      l=sum(v*v,2)                                   ! get square-norm
 
C------weight the randompoint in the unit d-cube
C      by (t/r)**dim, where t=max(x1,x2,x3,..,xdim) for d-dimensions
C                     where r=sqrt(l)=sqrt(x1^2+x2^2+...+xdim^2)
      hit=hit+(sqrt(t*t/l))**dim
 
C------use the same random numbers to Monte-Carlo the naive way
C
      where ( l < 1.0) oldhit=oldhit+1
 
 100  continue
C------------------------------------------|
      monterad(j)=2**dim*sum(hit)/(N*samp)
      oldmonterad(j)=2**dim*sum(oldhit)/(N*samp)
      write(3,*) j, monterad(j),oldmonterad(j)
D     print*, j, monterad(j),oldmonterad(j)
      mptm=mpTimerElapsed()
D     print*, 'time (sec): ', mptm/1000.0
 200  continue
C-----------------------------------------------------------|
 
C------------------------------STATISTICS
      mean=sum(monterad)/r
      oldmean=sum(oldmonterad)/r
      write(3,*) '----'
      write(3,*) r,' Runs, each with',N*samp,' samples'
      write(3,*) 'Spherical volume in ',dim,'-dimensions'
      sumdevsq=sum((mean-monterad)*(mean-monterad))
      oldsumdevsq=sum((oldmean-oldmonterad)*(oldmean-oldmonterad))
 300  continue
 
 
      write(3,*) oldmean,'=NAIVE      ',sqrt(oldsumdevsq)/n,
     *           ' std dev (mean)'
      write(3,*) mean,'=IMPORTANCE ',sqrt(sumdevsq)/n,' std dev (mean)'
      write(3,*) truevol,'=TRUE volume'
      write(3,*) '% Naive Error     =', 100.0*(oldmean-truevol)/truevol
      write(3,*) '% Importance Error=', 100.0*(mean-truevol)/truevol
 
      end
