Fitfun now simultaneously fits multiple datasets

Extensions to fitfun for simultaneously fitting multiple datasets

Fitfun - version 6.x, 7.x, 8.1 R. Ghosh, ILL, March 2013

Summary of intended use of mfitfun

This version of fitfun is compatible with earlier versions. In its present pre-compiled form up to 800 parameters can be set and stepped. Up to 20 sets of data (total maximum length 21000 x,y,error) can be treated. The principal objective was to retain the essential simplicity of use of fitfun with a greatly increased range of problems where it is useful to analyse simultaneously several sets of data to find an optimum set of parameters. The reading-in and calculation procedures have to be modified to use several datasets. New COMMON blocks are shown in the examples to allow the programmer to identify the internal values appropriate to each treated dataset. It has been necessary to extend the commands a little to deal with the new functionality.

Recent Changes

March 2013 version 8.1
The data arrays have been expanded to 21000 Other changes are noted here.

September 1999
Version 6.1 includes a larger total number of parameters. Now up to 40 parameters for all 20 sets of data can be accommodated.

mfitfun - introduction

This is an extended version of fitfun to cope with refining parameters some of which are pertinent to a number of sets of data, and parameters which are specific to each set. For neutron scattering experiments the primary use is in analysing several sets of data where isotope substitution, e.g., H20/D20 solvents, have been used for contrast variation in small-angle scattering or reflection measurements. Here the structures to be measured are unchanged, but the contrasts of the components are modified relative to the solvents in use. Alternative uses can be fitting Q dependent quasi-elastic scattering measured at several different angles.

On first entry only one set of data is expected. It is therefore still possible to use the present library without change for use with existing programs, though some of the COMMON blocks used internally in fitfun have now been rearranged and increased in size. An initial set of realistic parameters can be set using a new routine forfit.

The read-in routine supplied by the programmer is expected to provide one or more sequential sets of data; each set can have a short title.

The calculation has information available to change models depending in the set of data in use. For each set of data the local variables appropriate to the set are selected by calling a new routine ftgetp.

When data are displayed several colours are used to distinguish the different data sets, and the residuals from the fitting procedure.

The display of variables resulting from the V command now distinguishes the variables which are set member dependent by showing the number followed by S. These variables may be all be given the same value or defined selectively (or all given the same value, and then one or two modified later). Thus the V command has an extended syntax.

A simple example follows:

is1 7% mtst
 mtst               version  2.0
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>r
 MAXIMUM NUMBER OF DATA AND SETS ARE:  2100   98


 SET   1  GIVE INPUT FILENAME (OR BLANK) : r.spc
 TITLE:     LACU2GE2 1.5 K set#1                                

   41 DATA HAVE BEEN READ IN, TOTAL NOW   41


 SET   2  GIVE INPUT FILENAME (OR BLANK) : s.spc
 TITLE:     LACU2GE2 1.5 K set#2                              

   41 DATA HAVE BEEN READ IN, TOTAL NOW   82


 SET   3  GIVE INPUT FILENAME (OR BLANK) : t.spc
 TITLE:     LACU2GE2X 1.5 K set#3                               

   41 DATA HAVE BEEN READ IN, TOTAL NOW  123


 SET   4  GIVE INPUT FILENAME (OR BLANK) : 
 123 POINTS IN   3 SETS OF DATA HAVE BEEN READ IN
 THERE ARE  123 DATA POINTS IN CURRENT FITTING RANGE
 IN  3 SETS.
 Total number of variables :    9
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>d
v
 Fitfun 6.0  TITLE:    LACU2GE2X 1.5 K                               
 FITTING Y : COUNTS                 VERSUS X : CHANNEL             INPUT SET  1
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )     STEP     % DEVIATION 

    1 S CENTRE       311.2    (   311.2    )  0.1000E-01      0.00
    2   SIGMA        1.825    (   1.825    )   1.000          0.00
    3   FLAT BAK     10.00    (   10.00    )  0.0000E+00      0.00
    4   SLOPE       0.0000E+00(  0.0000E+00)  0.0000E+00      0.00
    5 S HEIGHT       5000.    (   5000.    )   10.00          0.00
  100     MAXIMUM STEP    100.0       300  SUBSTEP   0.01000
  200     ACCURACY      0.1000E-01
 THERE ARE 123 POINTS IN THE CURRENT RANGE SET BY "ONLY" COMMAND.
 CURRENT LIMITS (IF RANGE IS NON-ZERO) ARE : 
  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>f

 mtst>f 50

 FITTING.....

 .....ENDED
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>f
 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>a  
 Fitfun 6.0  TITLE:    LACU2GE2X 1.5 K                               
 FITTING Y : COUNTS                 VERSUS X : CHANNEL             INPUT SET  1
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )     STEP     % DEVIATION 

    1 S CENTRE       313.7    (   311.2    )  0.1000E-01      0.00
    2   SIGMA        1.736    (   1.825    )   1.000          0.44
    3   FLAT BAK     10.00    (   10.00    )  0.0000E+00      0.00
    4   SLOPE       0.0000E+00(  0.0000E+00)  0.0000E+00      0.00
    5 S HEIGHT       8384.    (   5000.    )   10.00          0.75
 FITTING Y : COUNTS                 VERSUS X : CHANNEL             INPUT SET  1
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )     STEP     % DEVIATION 

    1 S CENTRE       314.2    (   314.2    )  0.1000E-01      0.00
    2   SIGMA        1.736    (   1.825    )   1.000          0.44
    3   FLAT BAK     10.00    (   10.00    )  0.0000E+00      0.00
    4   SLOPE       0.0000E+00(  0.0000E+00)  0.0000E+00      0.00
    5 S HEIGHT       8384.    (   5000.    )  -105.0          0.75
 FITTING Y : COUNTS                 VERSUS X : CHANNEL             INPUT SET  1
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )     STEP     % DEVIATION 

    1 S CENTRE       324.1    (   324.0    )  0.1000E-01      0.00
    2   SIGMA        1.736    (   1.825    )   1.000          0.44
    3   FLAT BAK     10.00    (   10.00    )  0.0000E+00      0.00
    4   SLOPE       0.0000E+00(  0.0000E+00)  0.0000E+00      0.00
    5 S HEIGHT       7202.    (   6507.    )   10.00          0.74
  100     MAXIMUM STEP    100.0       300  SUBSTEP   0.01000
  200     ACCURACY      0.1000E-01
 THERE ARE 123 POINTS IN THE CURRENT RANGE SET BY "ONLY" COMMAND.
 CURRENT LIMITS (IF RANGE IS NON-ZERO) ARE : 
  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00
 mtst>f

 TYPE HELP OR OPTION: H,R,D,X,Y,V,A,F,O,P,L,W,C,E,J,S,T,Z
 mtst>e
 DO YOU WISH TO SAVE THE CURRENT PARAMETER VALUES ? (Y)
n
 %PGPLOT,  Closing file mtst163.ps

Main routine mtst.f
      PROGRAM MTST
C***** BASIC EXAMPLE ROUTINE FOR FITFUN
C----- multiple data sets
C HP    setenv LPATH /lib:/usr/lib:/usr/local/lib:/usr/lib/X11R5
C HP    f77 -o mtst mtst.o mread.o mgauss.o libmfitfun.a -lpgplot -lX11 -lU77
C IRIX  f77 -o mtst mtst.o mread.o mgauss.o libmfitfun.a libpgplot.a -lX11
      EXTERNAL MREAD,MGAUSS
C
C***** MREAD IS A USER SUPPLIED ROUTINE TO READ IN SETS OF SPECTRA EACH CALL
C
C***** MGAUSS IS A USER'S ROUTINE TO CALCULATE THE REQUIRED FUNCTION
C      FOR ALL SETs
C***** INPSET SHOWS WHICH PARAMETERS ARE SHARED BY ALL SETS
      COMMON/TITLES/NAMES(40),TX,TY
      COMMON/SETSET/INPSET(40)
      COMMON/TITLEP/NPARAS
      COMMON/WORK/W(30660)
C***** THIS IS SUFFICIENT FOR 1280 DATA AND 10 PARAMETERS !
      COMMON/VERSION/VERP
C***** VERP IS TYPED BY THE PROGRAM INDICATING THE CURRENTLY RUNNING VERSION
C
      CHARACTER*8 NAMES
      CHARACTER*4 PNAM
      CHARACTER*20 TX,TY
C
      DATA NAMES/'CENTRE  ','SIGMA','FLAT BAK','SLOPE',
     1'HEIGHT',35*' '/
      DATA INPSET/1,0,0,0,1,35*0/
c***** only the centre position and height paras are set dependent
c     shared parameters have 0, set specific 1
      DATA PNAM/'mtst'/
      DATA TX/'CHANNEL'/
      DATA TY/'COUNTS'/
      DATA NPARAS/5/
C
      VERP=2.0
      CALL FITFUN(PNAM,MREAD,MGAUSS)
C
      END

Read-in routine MREAD The main function of MREAD is to provide a contiguous sequence of input data into XOBS, YOBS, and YERR at the same time filling INDSET with an index showing which set of data is where.

      SUBROUTINE MREAD(NPNTS,XOBS,YOBS,YERR,NTEXT)
C***** PLACES SETS OF DATA IN SEQUENCE ALONG XOBS,YOBS etc,
C      MAPPING IS IN INDSET
C     NPNTS  INTEGER     out    Total number of data points read in
C     XOBS   REAL        out    Array of X-values
C     YOBS   REAL        out    Array of Y-values
C     YERR   REAL        out    Array of error in Y-values
C     NTEXT  TEXT        out    General title
C The following may be ignored if only one set of data is programmed
C     INDSET INTEGER     out    Array showing to which set each x,y
C                               belongs
C     SNAMES TEXT        out    Short title for each set of data
C     
      PARAMETER (MAXDAT=21000)
      PARAMETER (MAXSET=20)
      COMMON/SETSIN/INDSET(MAXDAT),INSET,LIMSET
      COMMON/INFILE/SNAMES(MAXSET)
      REAL X(180),Y(180),EY(180)
      DIMENSION XOBS(1),YOBS(1),YERR(1)
      CHARACTER*50 NTEXT
      CHARACTER*20 SNAMES
      CHARACTER*20 FNAME
      NTOT=0
      INSET=0
      II=0
      LL=0
      WRITE(6,1) MAXDAT,LIMSET
1     FORMAT(' MAXIMUM NUMBER OF DATA AND SETS ARE: ',2I5)
      DO 50 L=1,LIMSET       
      N=0
10    WRITE(6,11) LL+1
11    FORMAT(//' SET ',I3,'  GIVE INPUT FILENAME (OR BLANK) : ',$)
      READ(5,2,END=100) FNAME
C***** EXIT WITH CONTROL-D (VMS CONTROL-Z) IS POSSIBLE HERE
2     FORMAT(A)
      IF(FNAME.EQ.' ') GO TO 51
      J=0
      OPEN(UNIT=10,FILE=FNAME,ACCESS='SEQUENTIAL',STATUS='OLD',ERR=99)
C***** NOTE UNIT IS NOT 1
C***** SIMPLY FORMATTED FILE
      READ(10,3,ERR=99,END=99) NTEXT
3     FORMAT(A)
C***** THIS CAN SERVE AS A TITLE FOR THE SPECTRUM
      DO 5 J=1,180
      READ(10,6,ERR=99,END=99) X(J),Y(J),EY(J)
      N=J
6     FORMAT(3G)
5     CONTINUE
99    CLOSE(UNIT=10)
      IF(J.LE.5) GO TO 12
      LL=LL+1
      SNAMES(LL)=NTEXT
C***** SAVE TITLE OF THIS SET
      DO 20 M=1,N
      II=II+1
      INDSET(II)=LL
C***** set indicator for set number
      XOBS(II)=X(M)
      YOBS(II)=Y(M)
      YERR(II)=EY(M)
20    CONTINUE
      WRITE(6,9) NTEXT
9     FORMAT(' TITLE: ',A/)
      WRITE(6,14) N,II
14    FORMAT(1X,I4,' DATA HAVE BEEN READ IN, TOTAL NOW',I5)
      GO TO 50
12    WRITE(6,98)
98    FORMAT(' INSUFFICIENT DATA TO CONTINUE')
C***** DO NOT RETURN UNTIL SOME GOOD DATA ARE READ
      GO TO 10
50    CONTINUE
51    INSET=LL
      NPNTS=II
      WRITE(6,52) NPNTS,INSET
52    FORMAT(I4,' POINTS IN ',I3,' SETS OF DATA HAVE BEEN READ IN')
      WRITE(6,60)
60    FORMAT(' GIVE SUMMARY TITLE : ',$)
      READ(5,3) NTEXT
      RETURN
100   STOP 'END OF TERMINAL INPUT'
      END
Calculation routine MGAUSS
      SUBROUTINE MGAUSS(NPAR,PARM,NFIT,XUSE,YUSE,YRUSE,YCALC,F)
C
C  CALC OF F AND YCALC FOR GAUSSIAN
C
C  FOR EACH SET MUST COPY OUT CURRENT SETS FROM ALL PARAMETERS PARM 
C  INTO LOCAL PARAMETER ARRAY P 
      PARAMETER (MAXDAT=21000)
      COMMON/SETSIN/INDSET(MAXDAT),INSET,LIMSET
      DIMENSION PARM(NPAR),P(40)
      DIMENSION YCALC(NFIT),XUSE(NFIT),YUSE(NFIT)
     1,YRUSE(NFIT),F(NFIT)
C
C
C***** A LITTLE OPTIMISATION ATTEMPTED....
      MYSET=1
      CALL FTGETP(MYSET,PARM,P)
      DO 10 I=1,NFIT
      IND=INDSET(I)
      IF(IND.EQ.MYSET) GO TO 5
C***** NOT CURRENT PARAMETERS - GET THESE
      CALL FTGETP(IND,PARM,P)
      MYSET=IND
5     CONTINUE
C
C***** YCALC WILL CONTAIN CALCULATED FUNCTION FOR PLOTTING/LISTING
      YCALC(I)= P(5) *
     1   EXP(-0.5*((XUSE(I)-P(1))/P(2))**2)
     1  +P(3) + P(4)*XUSE(I)
      F(I)=YCALC(I)-YUSE(I)
C***** F WILL BE MINIMISED BY VA05A
 10   CONTINUE
C
C
      RETURN
      END

An alternative calculation routine mgauss using fitfun internal COMMONs

This allows existing calculations to be encapsulated easily, simplifying the enhancement to treat several sets of data.

      subroutine mgauss(NPAR,PARM,NFIT,XUSE,YUSE,YRUSE,YCALC,F)
C
C  CALC OF F AND YCALC FOR GAUSSIAN
C
C  FOR EACH SET MUST COPY OUT CURRENT SETS FROM ALL PARAMETERS PARM 
C  INTO LOCAL PARAMETER ARRAY P 
c 
      PARAMETER (MAXDAT=21000)
      PARAMETER (MAXSET=20)
      PARAMETER (MAXPAR=800)
      COMMON/SETLIM/ISETS(MAXSET),ISETE(MAXSET)    
c***** ISETS is an array of starting addresses for each set XUSE,YUSE etc
c      ISETE is an array of end addresses for each set XUSE, YUSE etc 
      COMMON/SETSIN/INDSET(MAXDAT),INSET,LIMSET
      DIMENSION PARM(NPAR),P(40)
      DIMENSION YCALC(NFIT),XUSE(NFIT),YUSE(NFIT)
     1,YRUSE(NFIT),F(NFIT)
C
      DO 10 NSET=1,INSET
      IST=ISETS(NSET)
      LST=ISETE(NSET)
C***** SET UP CURRENT PARAMETERS
      CALL FTGETP(NSET,PARM,P)
C
C
      DO 20 I=IST,LST
C***** CALCULATE FOR EACH INPUT VALUE IN SET NSET
C
C***** YCALC WILL CONTAIN CALCULATED FUNCTION FOR PLOTTING/LISTING
      YCALC(I)= P(5) *
     1   EXP(-0.5*((XUSE(I)-P(1))/P(2))**2)
     1  +P(3) + P(4)*XUSE(I)
      F(I)=YCALC(I)-YUSE(I)
C***** F WILL BE MINIMISED BY VA05A
 20   CONTINUE
C***** THEN TO NEXT SET...
C
C
 10   CONTINUE
C
C
      RETURN
      END

Extensions to commands compared to previous versions

New command A
Without a numerical value following all parameters for all sets are listed.

If A is followed by a set number then this set will be shown when using the V command.

Extensions

Command V n v s k l
This changes parameter number n to value v with an associated step s (about 1.0) if this is to be fitted. If s is negative then parameter n is tied to preceding parameter -s
To tie to a specific set parameter then use -s=-(100*set+par number). To set parameters only for sets of data from k to l then these should be given; if absent then parameters for all sets are updated. Note that it is easiest to set all to a specific value, then again set the value of a specific parameter rather than, more slowly, setting each in turn.
e.g.

     v 5 1000 10           parameter 5 is set to 1000 for all sets
     v 5 10000 100 5 5

rather than the more tedious operations to achieve the same, shown below

     v 5 1000 10 1 1
     v 5 1000 10 2 2
     v 5 1000 10 3 3
     v 5 1000 10 4 4
     v 5 10000 100 5 5

New Library routines

forfit

This routine can be called before the call to fitfun to set up a reasonable set of starting parameters, creating the file PNAM.ffn.

      SUBROUTINE FORFIT(PNAM,PARIN)
C***** PNAM IS PROGRAM NAME (4 CHARACTERS COMPATIBLE WITH FILENAME
C      WHICH WILL BE PNAM.ffn
C      PARIN IS SET OF TRIAL PARAMETERS WITH THE SEQUENCE AS IN NAMES
C      REAL PARIN(40)      IN
C      CHARACTER*4 PNAM    IN

The values the programmer provides in PARIN are placed in a new file PNAM.ffn if this file does not already exist. The initialisation is performed for all set parameters. The total number of parameters, 200, limits the number of sets which can be treated at present.

ftgetp

Using the internal information within fitfun this routine extracts a set of parameters appropriate for fitting a specific set of data.

      SUBROUTINE FTGETP(NSET,PARM,P)
C***** FOR ANY SET NSET RECONSTRUCTS AN APPROPRIATE SET OF PARAS P
C     NSET INTEGER IN  DATASET SEQENCE NUMBER
C     PARM REAL    IN  PARAMETERS FED TO CALSUB
C     P    REAL    OUT PARAMETERS TAKEN FROM PARM
C                      APPROPRIATE TO SET NSET FOR 
C                      CALCULATION

At present the pre-built routines are located in libraries.,