c
c    MINE710 - Reading Course - May 2000
c    DGM - Discrete Gaussian Model for change of support (global distribution)
c
c    This program allows the calculation of Hermite polynomials and coefficients,
c    change of support coefficient, and global distribution for block support.
c
c    Program written by: Julian Ortiz C.
c
      program DGM

      parameter(PI=3.1415926535,EXP=2.718281828,MAXDAT=11000,MAXP=5000)

      real y(MAXDAT),z(MAXDAT),zv(MAXDAT),H(MAXP+1,MAXDAT),
	+	 zc(MAXDAT),g(MAXDAT),fi(MAXP),a,b,r1,r2,rm,L
	character datafl*30,pointfl*30,dbgfl*30,outfl*30
      integer p,i

c
c    Open the file with the transformation table (e.g. nscore.trn) and
c    define some parameters required.
c
c    Varblock = Variance at block support, calculated using gammabar.
c
c    Epsilon = difference between the variance obtained using the Hermite 
c			 approximation and the block variance caluclated using gammabar.
c
c    np = number of Hermite polynomials to use
c
c
c
c
      write(*,*) 'DGM - Version 1.0 - May 2000'
	write(*,*) 'Discrete Gaussian Model for Change of Support'
	write (*,*)
      write(*,*) 'File with transformation table (e.g. nscore.trn): '
	read(*,*) datafl
	write (*,*)
	write(*,*) 'Variance at block support (calculated with gammabar):'
	read(*,*) varblock
	write (*,*)
	write(*,*) 'Acceptable error for the block variance calculated usi
	+ng Hermite polynomials (e.g. e=0.0001): '
	read(*,*) epsilon
	write (*,*)
	write(*,*) 'Number of Hermite polynomials to use (e.g. np=100): '
	read(*,*) np
	write (*,*)
	write(*,*) 'Output File to check the anamorphosis at point support
	+ (e.g. point.out): '
 	read(*,*) pointfl
	write (*,*)
	write(*,*) 'Debug File with Hermite Coefficients, fi(p) (e.g. dgm.
	+dbg): '
 	read(*,*) dbgfl
	write (*,*)
	write(*,*) 'DGM Output file with Z values at block support (e.g. d
	+gm.out): '
 	read(*,*) outfl
	write (*,*)

c    
c   Read datafl and store the variable and its transform in z and y 
c
	open(1,file=datafl) 
	open(2,file=pointfl)
      open(3,file=dbgfl)
	open(4,file=outfl)
	a=0.0
	b=1.0
      ndat=0
      do i=1,MAXDAT
        read(1,*,end=99) z(i),y(i)
        ndat=ndat+1
      end do
99    continue
c
c   Initialize fi(p)
c
      do p=1,np+1
        fi(p)=0
      end do
c
c   Calculate the Hermite Polynomials and store them in H(p,i), where i=y(i)
c
      do i=1,ndat
        H(1,i)=1
        H(2,i)=-y(i)
      end do
      do p=1,np-1
        do i=1,ndat
            H(p+2,i)=-y(i)*H(p+1,i)/(sqrt(real(p+1)))
	+		-sqrt(real(p))/sqrt(real(p)+1)*H(p,i)
        end do
      end do
c
c   Calculate the probability for the transformed values and store them in g(i)
c   
      do i=1,ndat
            g(i)=(1/sqrt(2*PI))*(EXP**(-(y(i)**2)/2))
      end do
c
c   Calculate the Hermite Coefficient depending on the order np
c
      do i=1,ndat
		fi(1)=fi(1)+z(i)/ndat
	end do
	do p=2,np+1
      do i=2,ndat
            fi(p)=fi(p)+(z(i-1)-z(i))*(1/sqrt(real(p-1)))*H(p-1,i)*g(i)
      end do
      end do
c
c    Check the anamorphosis
c
	do i=1,ndat
		do p=1,np+1
			zc(i)=zc(i)+fi(p)*H(p,i)
		end do
	end do
c
c   Write on pointfl the calculated z values (zc) and the original z (z)
c   to check the point support anamorphosis
c
	write(2,*)'DGM Output File to check anamorphosis at point support'
	write(2,*) '2'
	write(2,*) 'Approximated Value'
	write(2,*) 'Original Value'
      do i=1,ndat
     		write(2,*) zc(i),z(i)
	end do
	close(2)
c
c   Calculate the variance and write on dbgfl the values of fi(p) and Var
c
	write(3,*) 'Debugging file for DGM'
	write(3,*)
      write(3,*) ' Hermite coefficients:'
      var=0
      do p=1,np+1
		var=var+fi(p)**2
            write(3,*) 'fi(',p-1,')= ',fi(p)
      end do
	write(3,*)
      write(3,*) 'Point Variance (using Hermite polynomials)= '
	+,var-fi(1)**2
	write(3,*)
c
c   Calculate the change of support coefficient such that it minimizes the difference
c   between the sum of fi(p)^2*r^(2*p) and the variance for block support obtained 
c   using gammabar
c	
97   	nit=nit+1
	L=b-a
      rm=(a+b)/2
	r1=a+(L/4)
	r2=b-(L/4)
      varr1=0
	varr2=0
	varrm=0
	do p=2,np+1
		varr1=varr1+(fi(p)**2)*(r1**(2*(p-1)))  
		varr2=varr2+(fi(p)**2)*(r2**(2*(p-1)))  
		varrm=varrm+(fi(p)**2)*(rm**(2*(p-1)))  
	end do
	difr1=abs(varr1-varblock)
	difr2=abs(varr2-varblock)
	difrm=abs(varrm-varblock)
	if(difrm.lt.epsilon) then
		goto 98
	else
	if(difr1.lt.difrm) then
		b=rm
		goto 97
	else
	if(difr1.ge.difrm) then
		if(difr2.lt.difrm) then
			a=rm
			goto 97
		else
		if(difr2.ge.difrm) then
			a=r1
			b=r2
			goto 97
		else 
			write(*,*) 'ERROR'
		end if
		end if
	end if
	end if
	end if
98	continue
c
c    Write on dbgfl the coefficient and the number of iterations required
c
      write(3,*) ' Change of Support Coefficient: r = ',rm
	write(3,*)
	write(3,*) ' Number of Iterations: nit = ',nit
	write(3,*)
c
c    Calculate the anamorphosis for block support
c	  
	do i=1,ndat
		do p=1,np+1
			zv(i)=zv(i)+fi(p)*H(p,i)*rm**(p-1)
		end do
	end do
c
c    Write on outfl the z values and y values for block support
c
	write(4,*) 'DGM Output File'
	write(4,*) '2'
	write(4,*) 'Z values at block support'
	write(4,*) 'Y values (gaussian)' 
      do i=1,ndat
		write(4,*) zv(i),y(i)
	end do
c
c   Terminate Program
c	
      write(*,*) 'DGM - Program terminated'
	end         