least squares minimisation fortran 77

Question

Trying to get a one-parameter least squares minimisation working in fortran77. Here's the code; it compiles and seems to work except....it gets caught in an infinite loop between values of h1= 1.8E-2 and 3.5E-2. Having a look now but, odds are, I'm not going to have much luck sussing the issue on my own. All help welcome!

      PROGRAM assignment

      ! A program designed to fit experiemental data, using the method
      ! of least squares to minimise the associated chi-squared and
      ! obtain the four control parameters A,B,h1 and h2.
      !*****************************************************************

      IMPLICIT NONE
      INTEGER i
      DOUBLE PRECISION t(17),Ct(17),eCt(17)
      DOUBLE PRECISION h1loop1,h1loop2,deltah,Cs
      DOUBLE PRECISION chisqa,chisqb,dchisq

      OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
      DO i=1,17
          READ(21,*)t(i),Ct(i),eCt(i)
      END DO
      CLOSE(21)

      !Read in data.txt as three one dimensional arrays.
      !*****************************************************************
      !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
      !DO i=1,17
      !   WRITE(21,*)t(i),Ct(i),eCt(i)
      !END DO
      !CLOSE(21)
      !
      !Just to check input file is being read correctly.
      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1loop2= 0.001
      h1loop1= 0.0 !Use initial value of 0 to calculate start-point chisq
      DO 10
           chisqa= 0.0
          DO 20 i= 1, 17
              Cs= exp(-h1loop1*t(i))
              chisqa= chisqa + ((Ct(i) - Cs)/eCt(i))**2
 20       END DO

           chisqb= 0.0
          DO 30 i= 1, 17
              h1loop2= h1loop2 + deltah
              Cs= exp(-h1loop2*t(i))
              chisqb= chisqb + ((Ct(i) - Cs)/eCt(i))**2
 30       END DO

         !Print the two calculated chisq values to screen.
          WRITE(6,*) 'Chi-squared a=',chisqa,'for Lamda1=',h1loop1
          WRITE(6,*) 'Chi-squared b=',chisqb,'for Lamda1=',h1loop2

          dchisq= chisqa - chisqb
          IF (dchisq.GT.0.0) THEN
              h1loop1= h1loop2
          ELSE
              deltah= deltah - ((deltah*2)/100)
          END IF

          IF (chisqb.LE.6618.681) EXIT
 10   END DO

      WRITE(6,*) 'Chi-squared is', chisqb,' for Lamda1 = ', h1loop2
      END PROGRAM assignment

EDIT: Having looked at it again I've decided I have no clue what's screwing it up. Should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!

I don't really understand your question in the edit, your code is missing and end if for the dchisq.gt.0 if at the end of the loop, but what do you refer to with h1 and which neighbor do you mean? — haraldkl, Jan 26 '14 at 21:41
Sorry I missed this. You've probably worked it out by now but, for clarity- h1 is the variable being found via the minimisation, for equation Cs= A*exp(-h1*t) + B*exp(h2*t). I've split it into h1a and h1b here so I can use it in two DO loops. As for neighbour, I just mean a neighbouring value for h1 in the minimisation. Like- if h1 minimises to a point where h1b has a larger chi-squared than h1a, but chi-squared is still large, the program would need to shift back and to a different value of h1b so it doesn't get stuck. Reading that back it doesn't seem clear, but not sure how else to say it. — user3228501, Jan 26 '14 at 22:06
Hopefully the new code upload is easier to understand than I am! — user3228501, Jan 26 '14 at 22:09
You unconditionally set your deltah to 0.001 in your 20 loop. Which essentially renders your conditional setting to a negative value after the second loop kinda useless (its only active for the first i). — haraldkl, Jan 27 '14 at 07:04
Also, setting h1b always to 0.001 in the loops looks strange as well. Is there any reason, why you want to have the code in F77 fixed format? — haraldkl, Jan 27 '14 at 07:29
No good reason other than I'm still figuring out how to get this to work! I've got a newly updated code that'll probably be more to your liking- it compiles fine, however it's getting stuck between two values of chi-squared on an endless loop. Getting sick of this problem now- stayed up all night trying to fix it without much luck. I'll update the first post with the newest code. — user3228501, Jan 27 '14 at 11:18
have look at my updated answer, I changed your code to free format and put some proposed changes which remedies some things that looked weird to me. — haraldkl, Jan 27 '14 at 11:39

haraldkl · Answer 1 · 2014-01-27T07:33:23.820

do i=1

is not starting a loop, for a loop you need to specify an upper bound as well:

do i=1,ub

that's why you get the error message about the doi not having a type, in fixed format spaces are insignificant...

Edit: If you want to have an infinite loop, just skip the "i=" declaration completely. You can use an exit statement to leave the loop, when a certain criterion has been reached:

do
  if (min_reached) EXIT
end do

Edit2: I don't know why you stick to F77 fixed format. Here is your program in free format, with some fixes of places, which looked weird, without digging too much into the details:

PROGRAM assignment

  ! A program designed to fit experiemental data, using the method
  ! of least squares to minimise the associated chi-squared and
  ! obtain the four control parameters A,B,h1 and h2.
  !*****************************************************************

  IMPLICIT NONE

  integer, parameter :: rk = selected_real_kind(15)
  integer, parameter :: nd = 17

  integer :: i,t0

  real(kind=rk) :: t(nd),t2(nd),Ct(nd),eCt(nd),Ctdiff(nd),c(nd)
  real(kind=rk) :: Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
  real(kind=rk) :: deltah,Cs(nd)

  OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
  DO i=1,nd
     READ(21,*) t(i),Ct(i),eCt(i)
  END DO
  CLOSE(21)

  !Read in data.txt as three one dimensional arrays.
  !*****************************************************************
  !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
  !DO i=1,17
  !   WRITE(21,*)t(i),Ct(i),eCt(i)
  !END DO
  !CLOSE(21)
  !
  !Just to check input file is being read correctly.
  !*****************************************************************
  !****************************Parameters***************************

  Aa= 0
  Ba= 0
  h1a= 0
  h2a= 0

  !**********************Minimising Lamda1 (h1)*********************
  deltah= 0.001_rk
  h1b= deltah

  minloop: DO

    chisqa= 0
    DO i= 1,nd
       Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
       Ctdiff(i)= Ct(i) - Cs(i)
       c(i)= Ctdiff(i)**2/eCt(i)**2
       chisqa= chisqa + c(i)
       h1a= h1a + deltah
    END DO

    ! Use initial h1 value of 0 to calculate start-point chisq.

    chisqb= 0
    DO i= 1,nd
       h1b= h1b + deltah
       Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
       Ctdiff(i)= Ct(i) - Cs(i)
       c(i)= Ctdiff(i)**2/eCt(i)**2
       chisqb= chisqb + c(i)
    END DO

    ! First-step h1 used to find competing chisq for comparison.
    WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
    WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
    ! Prints the two calculated chisq values to screen.

    dchisq= chisqa - chisqb
    IF (dchisq.GT.0) THEN
       h1a= h1b
    ELSE IF (dchisq.LE.0) THEN
       deltah= (-deltah*2)/10
    END IF

    IF (chisqb.LE.6000) EXIT minloop

  END DO minloop

  WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b

END PROGRAM assignment

I gotcha. Any idea how to get it to iterate until it reaches a minimum? Like- assuming it's unknown how many iterations would be needed. — user3228501, Jan 26 '14 at 20:17
That's wonderful! Thankyou. Just when I needed it too, about ready to test the first loop. Wish me luck. Again, thanks :) — user3228501, Jan 26 '14 at 21:49
Infinite loop! I'll edit the main question, then going to have a good look at it and see if I can figure out why. — user3228501, Jan 26 '14 at 21:57

least squares minimisation fortran 77

1 Answers1