Final Project - Minimization

CodeNotes.txt

@@ -0,0 +1,7 @@
+Notes for fitting code:
+
+The source code is broken up into two files.  
+
+The folder Fitting Code contains the minimization routine and all of the subroutines therein.  The text files contain necessary information for the code to run.  fittingparm.txt gives the parameters that you wanted fitted - zeros indicate the parameter is not fit.  initialparm.txt is where you input your starting potentials.  xexp.txt has the experimental data, and limits.txt contains the parameter bounds.  
+
+Parameter Space Code contains the the code to construct the chi square contour plots.  makefile_grid will compile chigrid which outputs a file with a grid of parameters.  makefile_chi compiles an executable to calculate the chi square value at each point on this grid.  This grid is then run through parse.pl to create individual files for each 2d contour plot.  Contour plots are created with contourplots.py.

Fitting Code/Gradiant.f90

@@ -0,0 +1,117 @@
+   program gradmin
+   implicit none
+   real*8 parm(12),initparm(12),chi,chinew,step,mean(12),sigma(12),w(12)
+   integer nlines,io,i,nfit,j,k
+   real*8, allocatable :: data(:,:)
+   integer, allocatable :: fitparm(:)
+   real*8, allocatable :: grad(:)
+
+   ! define the step size
+   step = 0.0001d0
+
+   ! read in inital parameters
+   nfit = 0
+   open(unit=1,file="fittingparm.txt")
+   do i=1,12
+      read(1,*) parm(i)
+      if (parm(i)==0) then
+         nfit = nfit + 1
+      end if 
+   enddo 
+   close(1)
+   !print *, "nfit", nfit
+
+   allocate (fitparm(nfit))
+
+   ! read in the data points
+   open(unit=2,file="xexp.txt")
+   nlines = 0
+   do 
+      read(2,*,iostat=io)
+      if (io/=0) exit
+      nlines = nlines + 1
+   enddo 
+
+   rewind(2)
+
+   allocate (data(3,nlines))
+
+   do i=1,nlines
+      read(2,*) data(1,i),data(2,i),data(3,i)
+   enddo
+
+   close(2)
+
+   ! read in the initial parameters
+   j = 1
+   open(unit=3,file="initialparm.txt")
+   do i=1,12
+      read(3,*) initparm(i)
+      if (parm(i)==0) then
+         parm(i) = initparm(i)
+      else 
+         fitparm(j) = i
+	 !print *,fitparm(j),j
+	 j = j + 1
+      end if 
+      !print *, parm(i)
+   enddo 
+
+   ! read in the parameter bounds
+   open(unit=8,file="limits.txt")
+   do i=1,12
+      read(8,*) mean(i),sigma(i),w(i)
+   enddo 
+
+   ! count the number of iterations
+   k=0
+   do 
+
+   ! replace values in the input file
+   call inputfile(parm)
+
+   ! run fresco
+   call system("fresco < newfit.in > newfit.out")
+
+   ! calculate chi square
+   call chicalc(data,nlines,chi)
+
+   ! find the gradiant of chi square function
+   allocate (grad(nfit))
+   call calcgrad(parm,fitparm,nfit,data,nlines,grad,mean,sigma,w)
+
+   ! move the points - gradient descent
+   do i=1,12
+   do j=1,nfit
+      if (fitparm(j)==i) then
+         parm(i) = parm(i) - step*grad(j)
+      end if
+   enddo 
+   enddo 
+
+   deallocate (grad)
+
+   ! remake input file and run fresco
+   call inputfile(parm)
+   call system("fresco < newfit.in > newfit.out")
+
+   ! calculate chi square
+   call chicalc(data,nlines,chinew)
+   print *, k, chinew
+
+   ! need a condition for ending - if |newchi - chi| < some value
+   if (abs(chinew - chi) < 0.1d0) then
+      exit
+   end if 
+
+   k = k + 1
+   ! if a minimium isn't found, don't let the code run infinitely
+   if (k>500) exit
+   enddo
+
+   ! print final values
+   do i=1,12
+      print *, "#",parm(i)
+   enddo 
+
+   end program gradmin

Fitting Code/calcgrad.f90

@@ -0,0 +1,99 @@
+   real function gauss(m,s,x)
+   real*8 m,s,x,const,pi
+
+   pi = 4.d0*atan(1.d0)
+
+   !const = 1.d0/(sqrt(2.d0*pi*s**2))
+   const = 1.d0
+
+   gauss = const*exp(-(x-m)**2/(2.d0*s**2))
+
+   return
+   end
+
+   subroutine calcgrad(parm,fparm,nf,data,nl,grad,mean,sigma)
+   use constants
+   use F90library
+   implicit none
+   integer, intent(in) :: nf,nl,fparm(nf)
+   real*8, intent(inout) :: parm(12),grad(nf)
+   real*8, intent(in) :: data(nl),mean(12),sigma(12)
+   integer i,j
+   real*8 jac(nf,nf),chip,chim,step,delta,gaup,gaum,gauss
+   real*8, allocatable :: d(:)
+   real*8, allocatable :: e(:)
+
+   !do i=1,nf
+   !   print*, fparm(i)
+   !enddo 
+
+   ! for each parameter that's being fit, calculate chi^2 at \pm step
+   ! weight the chi^2 by a gaussian to keep the parameters physical
+   step = 0.1d0
+   do i=1,nf
+      ! origianal parm + step, calculate chi square value
+      parm(fparm(i)) = step + parm(fparm(i))
+      call inputfile(parm)
+      call system("fresco < newfit.in > newfit.out")
+      call chicalc(data,nl,chip)
+      gaup = gauss(mean(fparm(i)),sigma(fparm(i)),parm(fparm(i)))
+      !print *, gaup
+
+      ! original parm - step, calculate chi square value
+      parm(fparm(i)) = parm(fparm(i)) - 2*step
+      call inputfile(parm)
+      call system("fresco < newfit.in > newfit.out")
+      call chicalc(data,nl,chim)
+      gaum = gauss(mean(fparm(i)),sigma(fparm(i)),parm(fparm(i)))
+      !print *, gaum
+
+      ! put parameter back at it's original value and move to next
+      parm(fparm(i)) = parm(fparm(i)) + step
+      !print *, parm(fparm(i)),chip,chim
+
+      ! calculate derivative
+      grad(i) = (chip*gaup - chim*gaum)/(2*step)
+      !print *, grad(i),gaup,gaum,chip,chim
+   enddo 
+
+   ! form the approximate jacobian
+   !do i=1,nf
+   !   do j=i,nf
+   !      jac(i,j) = grad(i)*grad(j)
+   !      jac(j,i) = jac(i,j)
+	 !print *, jac(i,j),i,j
+   !   enddo 
+   !enddo 
+
+   ! diagonalize the jacobian and return it in grad
+   ! uses a diagonalization routine recieved from M. Hjorth-Jensen
+   ! define step size
+   !delta = 0.1
+   !allocate (e(nf))
+   !allocate (d(nf))
+   !call tred2(jac,nf,d,e)
+   !call tqli(d,e,nf,jac)
+   !do i=1,nf
+   !   grad(i) = d(i)
+   !   print *, grad(i)
+      !grad(i) = grad(i) + delta
+   !   grad(i) = sqrt(grad(i))
+   !enddo 
+
+   !do i=1,nf
+   !   print *, (jac(i,j),j=1,nf)
+   !enddo 
+
+   ! matrix inversion doesn't work - talk to Jenni about method
+   ! take the inverse of jac after going through tqli
+   !call matinv(jac,nf,1.d0)
+
+   ! reconstruct new parameters
+   !parm = matmul(jac,grad)
+   !do i=1,nf
+   !  print *, parm(i)
+   !enddo 
+   !deallocate(e)
+   !deallocate(d)
+
+   end subroutine calcgrad

Fitting Code/chicalc.f90

@@ -0,0 +1,34 @@
+   subroutine chicalc(data,n,chi)
+   implicit none
+   integer, intent(in) :: n
+   real*8, intent(in) :: data(3,n)
+   real*8, intent(inout) :: chi
+   real*8 res(n),angle,xsth
+   integer i,j
+
+   ! read in calculation from fresco
+   open(unit=8,file="fort.16")
+   ! first 13 lines are plotting options
+   do i=1,13
+      read(8,*)
+   enddo 
+   ! next 181 lines contain data
+   j=1
+   do i=0,180
+      read(8,*) angle, xsth
+      ! if the angles between theory and experiment/data are the same
+      ! calculate the residual
+      if (idnint(data(1,j))==i) then
+         res(j) = (xsth-data(2,j))/data(3,j)
+	 j = j + 1
+      end if 
+   enddo 
+   close(8)
+
+   ! calculate chisquare value
+   chi = 0
+   do i=1,n
+      chi = chi + res(i)**2
+   enddo 
+
+   end subroutine chicalc