Fortran77 woes- single parameter least squares minimisation stuck

In summary: Parameters*************************** Aa= 0 Ba= 0 h1a= 0 h2a= 0 OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW') DO i=1,17
  • #1
NathanM
14
0
Hey guys! Having issues with a code I'm writing. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

Here's my code, bless me with your insight!

Code: 
      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,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
      DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
      DOUBLE PRECISION deltah,Cs(17)

        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.
      !*****************************************************************
      !****************************Parameters***************************
      
      Aa= 0
      Ba= 0
      h1a= 0
      h2a= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17
                 h1a= h1a + deltah
                   Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
                     Ctdiff= Ct - Cs
                       c= Ctdiff**2/eCt**2
                         chisqa= chisqa + c(i)
                 deltah= 0.001
 20           END DO
      !Use initial h1 value of 0 to calculate start-point chisq.

      h1b= 0.001
      chisqb= 0
              DO 30 i= 1,17
                   h1b= h1b + deltah
                     Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
                       Ctdiff(i)= Ct(i) - Cs(i)
                         c= Ctdiff(i)**2/eCt(i)**2
                           chisqb= chisqb + c(i)
 30           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)/100

              END IF
           IF (chisqb.LE.6000) EXIT
 10        END DO
           WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
      !*****************************************************************
           
      END PROGRAM assignment
Thanks for the help.
Nathan

EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!EDIT: I've read this is the wrong place for this. Sorry, have opened a thread in the correct place, please remove this one.
 
Last edited:
Technology news on Phys.org
  • #2
Code: 
     1	PROGRAM assignment
     2	
     3	! A program designed to fit experiemental data, using the method
     4	! of least squares to minimise the associated chi-squared and
     5	! obtain the four control parameters A,B,h1 and h2.
     6	!*****************************************************************
     7	
     8	IMPLICIT NONE
     9	INTEGER i,t0
    10	DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
    11	DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
    12	DOUBLE PRECISION deltah,Cs(17)
    13	
    14	OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
    15	DO i=1,17
    16	READ(21,*)t(i),Ct(i),eCt(i)
    17	END DO
    18	CLOSE(21)
    19	
    20	!Read in data.txt as three one dimensional arrays.
    21	!*****************************************************************
    22	!OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
    23	!DO i=1,17
    24	! WRITE(21,*)t(i),Ct(i),eCt(i)
    25	!END DO
    26	!CLOSE(21)
    27	!
    28	!Just to check input file is being read correctly.
    29	!*****************************************************************
    30	!****************************Parameters***************************
    31	
    32	Aa= 0
    33	Ba= 0
    34	h1a= 0
    35	h2a= 0
    36	
    37	!*****************************************************************
    38	!**********************Minimising Lamda1 (h1)*********************
    39	deltah= 0
    40	DO 10
    41	    chisqa= 0
    42	    DO 20 i= 1,17
    43	        h1a= h1a + deltah
    44	        Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
    45	        Ctdiff= Ct - Cs
    46	        c= Ctdiff**2/eCt**2
    47	        chisqa= chisqa + c(i)
    48	        deltah= 0.001
    49	    20 END DO
    50	    
    51	    !Use initial h1 value of 0 to calculate start-point chisq.
    52	    h1b= 0.001
    53	    chisqb= 0
    54	    DO 30 i= 1,17
    55	        h1b= h1b + deltah
    56	        Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
    57	        Ctdiff(i)= Ct(i) - Cs(i)
    58	        c= Ctdiff(i)**2/eCt(i)**2
    59	        chisqb= chisqb + c(i)
    60	    30 END DO
    61	
    62	    !First-step h1 used to find competing chisq for comparison.
    63	    WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
    64	    WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
    65	    !Prints the two calculated chisq values to screen.
    66	    dchisq= chisqa - chisqb
    67	    IF (dchisq.GT.0) THEN
    68	        h1a= h1b
    69	    ELSE IF (dchisq.LE.0) THEN
    70	        deltah= (-deltah*2)/100
    71	    END IF
    72	    IF (chisqb.LE.6000) EXIT
    73	10 END DO
    74	WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
    75	!*****************************************************************
    76	
    77	END PROGRAM assignment
    78	
    79	Read more: https://www.physicsforums.com

First question...are you compiling with Fortran 77 or Fortran 90? I think 90.

I say this because you have typos in the program that may be doing something you did not intend. Recall that Fortran 90 operates on entire arrays when instructed to do so, i.e., no index is passed.

Using numbered linsting above:

line 44 calculates i-th index of Cs and only this index
line 45 performs an array operation!...is this what you want?
line 46 performs an array operation.
line 47 takes a single index from c

what is the point of performing arrays operations if ultimately you seem to just need a single index?

deltah is set to 0.001 at the end of do-loop 20

line 52 sets h1b to 0.001 and
line 53 sets chisqb to 0

this means that do-loop 30 does exactly the same thing on every iteration!

in do-loop 30, you seem to correctly use parenthesis and the i index in your operations, yet, in
line 58 you assign a single scalar to the entire 'c' array (since you did not specify an index)
line 59 will pick up the intended value and assign to chisqb

line 70 may flip the sign of deltah and divide it by 50, but this new value will be used only once at the beginning of do-loop 20; is this what you want? a short step back? just asking
 
  • #3
Hi! Thanks for the reply :) I'm compiling with a program called Force2.0, but it's definitely meant to be fortran 77. With skills like mine, wouldn't be surprised if I've screwed up a few times!

I'll bullet point responses to keep things clear(ish).
(44)- Riiiight. So, having i=1,17 there isn't running through all three equations and then summing? Crap.
(45,46,47)- Yeah I think so- this whole do block is supposed to be calculating chi-squared for values of t(i)=1,17 (ie all seventeen values of t read in from the data file), then adding them all together to find the actual chi-squared for that particular value of lamda 1 (h1 as I've entered it).

NB: I can't believe I didn't notice that the middle two equations are missing their array assignments. Dumbass. All four should be using arrays.

- Yeah I wasn't sure of the placement of deltah at the time of posting, I've moved it outside of the main loop since then. I'll post my newest code at the end of this reply. Likewise with h1b- that's been moved out the way too. I've left chisqb there as it makes sense in my head that chisq should be zeroed each time.

- Another delightful cockup on my part, thanks for pointing it out. These missing arrays must be from when I imported them from my original attempt (tried using functions and subroutines, did not end well), but still can't believe I (nor my friends who had a quick look) didn't notice them.

- hmm no. The idea was it's supposed to run through the whole loop (only operates in the second inner loop in the new code), knock the step back and downsize, but only for that iteration (unless dchisq is negative again). I've edited it to take two steps back now and decrease step size further (though I may change the size of the decrease once the minimisation loop is working).

Here's the newest code (I've added back in the missing array identifiers).

Code: 
      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,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
      DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
      DOUBLE PRECISION deltah,Cs(17)

        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.
      !*****************************************************************
      !****************************Parameters***************************
      
      Aa= 0
      Ba= 0
      h1a= 0
      h2a= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1b= 0.001
      h1a= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17
                 h1a= h1a
                 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)
 20           END DO
      !Use initial h1 value of 0 to calculate start-point chisq.

      chisqb= 0
              DO 30 i= 1,17
                 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)
 30           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-((deltah*2)/100)
              END IF
              
           IF (chisqb.LE.6618.681) EXIT
 10        END DO
           WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
      !*****************************************************************
           
      END PROGRAM assignment
Thankyou for the help :)

EDIT: I imagine I'm going to have to use a function for the chi-squared calculation aren't I?
 
Last edited:
  • #4
NathanM said:
Hey guys! Having issues with a code I'm writing. It doesn't fit the typical format and I may have already accidentally posted it in the wrong forum, so here's a copy of what was there with the newest version of the code. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

Here's my code, bless me with your insight!
Code: 
      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,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
      DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
      DOUBLE PRECISION deltah,Cs(17)

        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.
      !*****************************************************************
      !****************************Parameters***************************
      
      Aa= 0
      Ba= 0
      h1a= 0
      h2a= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1b= 0.001
      h1a= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17   !  *** Is part of the line below supposed to be commented out?
                   Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
                     Ctdiff= Ct - Cs  ! *** This is wrong. All three variables are arrays.
                       c= Ctdiff**2/eCt**2
                         chisqa= chisqa + c(i)
 20           END DO
      !Use initial h1 value of 0 to calculate start-point chisq.

      chisqb= 0
              DO 30 i= 1,17
                   h1b= h1b + deltah   ! *** In the line below, are you commenting out part of the line?
                     Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
                       Ctdiff(i)= Ct(i) - Cs(i)
                         c= Ctdiff(i)**2/eCt(i)**2
                           chisqb= chisqb + c(i)
 30           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-((deltah*2)/100)

              END IF
           IF (chisqb.LE.6618.681) EXIT
 10        END DO
           WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
      !*****************************************************************
           
      END PROGRAM assignment

Thanks for the help.
Nathan

EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!

I added some comments where I think you're going wrong. My comments look like this: ! ***
In one comment, you have three assignment statements where you are using array variables, but without indexes. This is wrong. Everywhere you have an array variable in a calculation, it should have an index. E.g., it should be Ct(i), not just Ct.

In my other comments, you have a line of code that might be wrong. It looks like you commented it out, but I'm not sure what you did.

Your variable names are terrible! For example, Aa,Ab,Ba,Bb,h1a,h1b. Can't you come up with some names that are self-explanatory? Also, why does your code meander all across the page? This makes it harder to see which statements are inside IF blocks and DO loops and such.
 
  • #5
Hi! I'll reply using bullet points for ease of reading. As always, thanks for taking the time to read through my mess :)
Comments first!

- RE: commented out variables- yeah I've commented them out whilst I work on the single parameter minimisation. They'll come back in one by one once it's working.
- You caught my cockup too. No clue how I managed to go the last few hours without noticing. Guess I just didn't expect there to be issues with those blocks and ignored them- foolish.As for variable names, well, they made sense to me at the time but if you'd prefer something clearer I can alter them. RE: meandering code, I assume you mean the equation blocks? They weren't like that before, got whinged at by some here when they had a look so changed it, but I'll change it back.
If you mean the code in general, I thought it was okay. Will alter it if you have any specific preferences.

EDIT: Sorry about the thread thing.
 
Last edited:
  • #6
When you post code, it's good practice here at PF to put [code] ... [/code] tags around it. I've done this for you in your code.
NathanM said:
Hi! Thanks for the reply :) I'm compiling with a program called Force2.0, but it's definitely meant to be fortran 77. With skills like mine, wouldn't be surprised if I've screwed up a few times!

I'll bullet point responses to keep things clear(ish).
(44)- Riiiight. So, having i=1,17 there isn't running through all three equations and then summing? Crap.
(45,46,47)- Yeah I think so- this whole do block is supposed to be calculating chi-squared for values of t(i)=1,17 (ie all seventeen values of t read in from the data file), then adding them all together to find the actual chi-squared for that particular value of lamda 1 (h1 as I've entered it).

NB: I can't believe I didn't notice that the middle two equations are missing their array assignments. Dumbass. All four should be using arrays.

- Yeah I wasn't sure of the placement of deltah at the time of posting, I've moved it outside of the main loop since then. I'll post my newest code at the end of this reply. Likewise with h1b- that's been moved out the way too. I've left chisqb there as it makes sense in my head that chisq should be zeroed each time.

- Another delightful cockup on my part, thanks for pointing it out. These missing arrays must be from when I imported them from my original attempt (tried using functions and subroutines, did not end well), but still can't believe I (nor my friends who had a quick look) didn't notice them.

- hmm no. The idea was it's supposed to run through the whole loop (only operates in the second inner loop in the new code), knock the step back and downsize, but only for that iteration (unless dchisq is negative again). I've edited it to take two steps back now and decrease step size further (though I may change the size of the decrease once the minimisation loop is working).

Here's the newest code (I've added back in the missing array identifiers).
Code: 
           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,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
      DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
      DOUBLE PRECISION deltah,Cs(17)

        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.
      !*****************************************************************
      !****************************Parameters***************************
      
      Aa= 0
      Ba= 0
      h1a= 0
      h2a= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1b= 0.001
      h1a= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17
                 h1a= h1a
                 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)
 20           END DO
      !Use initial h1 value of 0 to calculate start-point chisq.

      chisqb= 0
              DO 30 i= 1,17
                 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)
 30           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-((deltah*2)/100)
              END IF
              
           IF (chisqb.LE.6618.681) EXIT
 10        END DO
           WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
      !*****************************************************************
           
      END PROGRAM assignment
Thankyou for the help :)

EDIT: I imagine I'm going to have to use a function for the chi-squared calculation aren't I?
Probably be a good idea.
NathanM said:
Hi! I'll reply using bullet points for ease of reading. As always, thanks for taking the time to read through my mess :)
Comments first!

- RE: commented out variables- yeah I've commented them out whilst I work on the single parameter minimisation. They'll come back in one by one once it's working.
- You caught my cockup too. No clue how I managed to go the last few hours without noticing. Guess I just didn't expect there to be issues with those blocks and ignored them- foolish.
Is your code producing the right results now?
NathanM said:
As for variable names, well, they made sense to me at the time but if you'd prefer something clearer I can alter them.
Don't do it for me - do it for yourself. When you review this code in a month or two, having meaningless variable names will make it harder for you to understand the code, even though you wrote it yourself.
NathanM said:
RE: meandering code, I assume you mean the equation blocks?
We call them assignments statements.
NathanM said:
They weren't like that before, got whinged at by some here when they had a look so changed it, but I'll change it back.
If you mean the code in general, I thought it was okay. Will alter it if you have any specific preferences.

EDIT: Sorry about the thread thing.

Minor nit, but this Greek letter - λ - is written "lambda".
 
  • #7
Mark44 said:
When you post code, it's good practice here at PF to put [code] ... [/code] tags around it. I've done this for you in your code.

Done!

Probably be a good idea.

So as the assignment statements stand right now, are they not running through the calculation for every array labelled parameter using the i=1,17 statement? If no, it'll be better to move all of the assignment statements to a function and just execute that 17 times, yes? I'm not great with functions; going to throw one together quick and post it in an edit here, help with any needed corrections would be appreciated.

Is your code producing the right results now?

No, still stuck at two lamda values.

Don't do it for me - do it for yourself. When you review this code in a month or two, having meaningless variable names will make it harder for you to understand the code, even though you wrote it yourself.

Minor nit, but this Greek letter - λ - is written "lambda".

This is true. What might be a great idea at the time might not be later. Not particularly creative, but I've just renamed them after the do loop they're used in.
I know re: lamda, but I didn't want "lamda1, lamda2" lengthening lines out. I usually use h in its place in fortran, since the actual character isn't available.
 
  • #8
I've written this function to replace the assignment statements in the two interior loops, leaving the summation statements in.

Code: 
              DO 20 i= 1,17
                 h1a= h1a
                 chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),Cs(i),eCt(i))
 20           END DO

      DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
      IMPLICIT NONE
      DOUBLE PRECISION Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt

                 Cs= exp(-h1loop1*t)!*Aa !+ Ba*exp(-h2loop1*t(i))
                 Ctdiff= Ct - Cs
                 c= Ctdiff**2/eCt**2

      RETURN
      END

Will this do the trick?

EDIT: It doesn't like the array indexes in the function call-in in the main program body. Not quite sure how to properly call it in...do I just assign the indexes inside the function itself (swear somebody told me you couldn't use them in functions) and write the call-in as c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)?
 
Last edited:
  • #9
There are several things wrong with what you're doing. The DO loop is presumably part of your main program. It would be in some code that looks like this:
Code: 
PROGRAM Whatever

! Declarations
   ...
   DO 20 i= 1,17
       h1a= h1a
       chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),Cs(i),eCt(i))
20 END DO
   ...
END PROGRAM Whatever

In the code above you have a completely useless assignment statement: h1a = h1a. Here you are assigning the value of h1a to itself. That value is already there, so there's no point in doing this.

In your main program you need a declaration for your c function. BTW, you really should have a better name for it than 'c'. For a function, a reasonable name is a noun that describes the value that is being returned.

Below the main program are any subprograms, including your function, which is one type of subprogram (subroutine being the other).

The main program has the statement that is the call to the function. (In 40+ years I've never heard anybody refer to it as a 'call-in'.) Your function call looks OK - it has the right number of parameters, and they appear to be the correct types. After you have your function declaration in place, I think things should work properly.

It's lamBda - there's a 'B' in the name.
 
  • #10
I wasn't sure if h1a would need stating before the loop or not, now I do.
c works for me. It's an intermediate value that has no use outside of its function, and it's only definition is as the sum of c(i) array- Chi-squared.
I like to be different. That's what it does so works for me to call it that, I'm a physicist first not a programmer.

Here's the code as it now stands.

Code: 
      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,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17)
      DOUBLE PRECISION Aloop1,Aloop2,Bloop1,Bloop2,h1loop1,h1loop2
      DOUBLE PRECISION deltah,Cs(17),h2loop1,h2loop2
      DOUBLE PRECISION chisqa,chisqb,dchisq,c(17,17,17,17,17,17,17)

        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.
      !*****************************************************************
      !****************************Parameters***************************

      Aloop1= 0
      Bloop2= 0
      h1loop1= 0
      h2loop2= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1loop2= 0.001
      h1loop1= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17
                 chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
     >           Cs(i),eCt(i))
 20           END DO
        WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1loop1
      !Use initial h1 value of 0 to calculate start-point chisq.
      
      chisqb= 0
              DO 30 i= 1,17
                 h1loop2= h1loop2 + deltah
                 chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
     >           Cs(i),eCt(i))
 30           END DO
      !First-step h1 used to find competing chisq for comparison.
        WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1loop1
        WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1loop2
      !Prints the two calculated chisq values to screen.
                 dchisq= chisqa - chisqb
              IF (dchisq.GT.0) THEN
                 h1loop1= h1loop2
              ELSE IF (dchisq.LE.0) THEN
                 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

      !*****************************************************************
      !***************************FUNCTIONS*****************************
      DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
      IMPLICIT NONE
      DOUBLE PRECISION Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt

                 Cs= exp(-h1loop1*t)!*Aa !+ Ba*exp(-h2loop1*t(i))
                 Ctdiff= Ct - Cs
                 c= Ctdiff**2/eCt**2

      RETURN
      END FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)

Just took a stab at defining c at the start of the program, wasn't entirely sure but there you go. This presents the following errors on compiling.test.f:46.36:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:46.42:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:46.50:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
C:\Users\Nathan\Desktop\Documents\University of Birmingham MSc Physics and Technology of Nuclear Reactors\Lectures\fortran\test.f:46.55:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:46.65:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:46.71:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:47.23:

> Cs(i),eCt(i))
1
Warning: Extension: REAL array index at (1)
test.f:55.36:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:55.42:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:55.50:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:55.55:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:55.65:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:55.71:

chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:56.23:

> Cs(i),eCt(i))
1
Warning: Extension: REAL array index at (1)
test.f:1.72:

PROGRAM assignment
1
test.f:78.72:

DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
2
Error: Two main PROGRAMs at (1) and (2)
test.f:13: error: size of variable 'c' is too large

Surprised by the last one. I thought fortran 77 could handle seven dimensions to an array. Though perhaps I've just not entered something correctly.
 
  • #11
O.k., stay away from functions, for now. But please, just because you are getting help here, don't stop studying Fortran...you really need to for this assignment and as a physicist!

It seems that your comments in the source apply to the code above them...that is not how it is typically done; instead, you first tell the reader what you are about to do, then you do it.

From what I can tell, the only variables that need to be arrays are the ones you read from the file; so, I cleaned that.

Deleted unnecessary stuff to force you to add it later with consistent variable naming.

Changed the program a bit to force you to re-read it and focus on the correctness, once again.

And that is it, I am not going to bother to explain EVERYTHING that is wrong with your latest code (you really went backwards); instead, please take a good look at the cleaned-up code that I post below...it at least compiles; but, as simple as it might be, I have not check the math at all or the looping or anything, I just looked at the Fortran.

Code: 
      !
      ! 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.      
      !
      PROGRAM assignment
      
      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)

      !**********************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
 
  • #12
lol oh don't worry, I'm going to be getting a lot more familiar with fortran77 over the next few years. We're only really getting into it now, but it'll be used in my lab reports later in the semester and my dissertation over the summer. Then, I hear reactor control is still done in fortran77...there'll be no escape!

No probs, I'll move the comments around. I'll compare the new code to my own latest and learn from it, thanks. Updated my code, the iteration is still getting stuck between two values. This problem is ridiculous.
 
Last edited:
  • #13
NathanM said:
Just took a stab at defining c at the start of the program, wasn't entirely sure but there you go. This presents the following errors on compiling.


test.f:46.36:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
1
Warning: Extension: REAL array index at (1)
test.f:46.42:

chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),

You have defined 'c' as a DOUBLE PRECISION array in the main program and then again as a DOUBLE PRECISION function. Which is it? That's why you are getting a bunch of compiler errors.
 
  • #14
I would recommend reducing the number of points from 17 down to 3 and performing a hand calculation...I mean, really, hand calculate values one iteration at a time to really see what your equations are yielding, in particular, when to increment h1loop1 and h1loop2.
 

1. What is Fortran77 and why is it causing issues with single parameter least squares minimisation?

Fortran77 is a programming language commonly used for scientific and engineering applications. It is causing issues with single parameter least squares minimisation because it is an older language that does not have the same capabilities and efficiency as newer languages, making it difficult to implement complex algorithms such as least squares minimisation.

2. Why is single parameter least squares minimisation important in scientific research?

Single parameter least squares minimisation is important in scientific research because it is a method used to fit a line or curve to a set of data points, allowing researchers to analyze and model relationships between variables. It is often used in fields such as statistics, physics, and economics to make predictions and draw conclusions from data.

3. What are some common issues encountered when trying to perform least squares minimisation in Fortran77?

Some common issues encountered when trying to perform least squares minimisation in Fortran77 include difficulties with handling arrays and matrices, limited precision and accuracy, and a lack of built-in functions and libraries for performing advanced mathematical operations.

4. Are there any workarounds or solutions for the Fortran77 woes in single parameter least squares minimisation?

Yes, there are some workarounds and solutions for the Fortran77 woes in single parameter least squares minimisation. One option is to use a newer programming language that has better support for mathematical operations and algorithms, such as Python or MATLAB. Another option is to update the Fortran77 code to use newer features and libraries, or to rewrite it in a more modern version of Fortran.

5. What are some resources for further understanding and troubleshooting Fortran77 woes in single parameter least squares minimisation?

Some resources for further understanding and troubleshooting Fortran77 woes in single parameter least squares minimisation include online forums and communities for Fortran programmers, textbooks and guides on Fortran programming, and seeking guidance from experienced Fortran programmers or scientific researchers who have encountered similar issues.

Similar threads

Comp Sci Fortran77 woes- single parameter least squares minimisation stuck
    • Engineering and Comp Sci Homework Help
Replies
1
Views
1K
Fortran Errors with fortran code, help
    • Programming and Computer Science
Replies
7
Views
3K
Least Squares Estimation for two parameters
    • Engineering and Comp Sci Homework Help
Replies
8
Views
1K

Hot Threads

Recent Insights

Back
Top