home *** CD-ROM | disk | FTP | other *** search
/ PC-Blue - MS DOS Public Domain Library / PC-Blue MS-DOS Public Domain Library - NYACC.iso / vol030 / nonlin.bas < prev    next >
Encoding:
BASIC Source File  |  1987-01-11  |  30.3 KB  |  525 lines
  1. 1000 '*************************************************************************
  2. 1010 '* PROGRAM NONLIN  version 3.0  5/1/82                                   *
  3. 1020 '*                                                                       *
  4. 1030 '* By:       Dave Whitman                                                *
  5. 1040 '*                                                                       *
  6. 1050 '*           Box 383 Baker Lab                                           *
  7. 1060 '*           Cornell U.                                                  *
  8. 1070 '*           Ithaca, NY  14853                                           *
  9. 1080 '*           (607) 256-6479                                              *
  10. 1090 '*                                                                       *
  11. 1100 '* Performs non-linear least squares analysis to determine parameters    *
  12. 1110 '* to fit any function to observed data.                                 *
  13. 1120 '*                                                                       *
  14. 1130 '* Inspired by a FORTRAN program by C.F. Wilcox, which was in turn       *
  15. 1140 '* based on "Rigorous Least Squares Adjustment"; Wentworth, W. E.        *
  16. 1150 '* J . Chem Ed. 42, 96 (1965).                                           *
  17. 1160 '*                                                                       *
  18. 1170 '*************************************************************************
  19. 1180 '* The program requires two input files to work:                         *
  20. 1190 '*                                                                       *
  21. 1200 '*    The first file contains BASIC code which is specific to the        *
  22. 1210 '*    paticular function being fit.  This code will be automatically     *
  23. 1220 '*    merged into the program at run time.  For the merge to work        *
  24. 1230 '*    properly, this file MUST BE SAVED IN ASCII FORMAT.  Otherwise, a   *
  25. 1240 '*    "Bad File Mode" error occurs.                                      *
  26. 1250 '*    To save in ascii mode, use the "a" option in your save command:    *
  27. 1260 '*    Example:  SAVE"function",a                                         *
  28. 1270 '*                                                                       *
  29. 1280 '*    When nonlin prompts: FUNCTION?, input the name of this file.       *
  30. 1290 '*                                                                       *
  31. 1300 '*    Code required:                                                     *
  32. 1310 '*       A line 10050 which sets the variables nvar [# of independant    *
  33. 1320 '*       variables] and np [# of parameters]. Example:                   *
  34. 1330 '*                                                                       *
  35. 1340 '*   10050 nvar = 1: np = 3                                              *
  36. 1350 '*                                                                       *
  37. 1360 '*       A subroutine which fills the array fcalc(nobs) with values of   *
  38. 1370 '*       the function, given observed values of the variables in array   *
  39. 1380 '*       vobs(nobs,nvar) and current guesses at the parameters in p(np). *
  40. 1390 '*       The subroutine should start at line 11000 and end at or before  *
  41. 1400 '*       line 11999. Example: (fitting data to a parabola)               *
  42. 1410 '*                                                                       *
  43. 1420 '*   11000 FOR I = 1 to NOBS                                             *
  44. 1430 '*   11010    fcalc(i) = p(1) * vobs(i,1)^2 + p(2) * vobs(i,1) + p(3)    *
  45. 1440 '*   11020    NEXT I                                                     *
  46. 1450 '*   11030 RETURN                                                        *
  47. 1460 '*                                                                       *
  48. 1470 '*                                                                       *
  49. 1480 '*    Optional code:                                                     *
  50. 1490 '*       Before stopping, nonlin makes a subroutine call to line 20000.  *
  51. 1500 '*       Normally, this subroutine consists of a stub, containing only   *
  52. 1510 '*       a return statement. If any final processing using the least     *
  53. 1520 '*       squares parameter set in p(np) or the variance/covariance       *
  54. 1530 '*       matrix in b(np,np) is desired, supply an appropriate            *
  55. 1540 '*       subroutine starting at line 20000. The only restriction on the  *
  56. 1550 '*       length of this subroutine is the available memory.              *
  57. 1560 '*                                                                       *
  58. 1570 '* The second file contains the observed data to be fit.                 *
  59. 1580 '* When nonlin prompts: INPUT FILE?, enter the name of this file.        *
  60. 1590 '*                                                                       *
  61. 1600 '*    Required data:                                                     *
  62. 1610 '*       nobs: the number of observations being fit                      *
  63. 1620 '*       numit: the # of iterations of the fitting process to be         *
  64. 1630 '*              performed. [ 5-10 is generally sufficient ]              *
  65. 1640 '*       iuserwt: a flag. If iuserwt = 1, nonlin expects all observed    *
  66. 1650 '*                function and variable values to be followed by a       *
  67. 1660 '*                weighting factor. If iuserwt = 0, nonlin automatically *
  68. 1670 '*                sets all initial weights to 1.                         *
  69. 1680 '*      internalwt: a flag. If internalwt = 1, nonlin estimates          *
  70. 1690 '*                weighting factors for each function value based on     *
  71. 1700 '*                the slope of the function.  If internalwt = 0, no      *
  72. 1710 '*                estimate.  Use with care, internal weighting often     *
  73. 1720 '*                causes divergance. Start with internalwt = 0.          *
  74. 1730 '*                                                                       *
  75. 1740 '*      fract: The fraction of the calculated change to apply to each    *
  76. 1750 '*             of the parameters.  Use to restrict changes when function *
  77. 1760 '*             trys to diverge. Normally equal to 1.                     *
  78. 1770 '*                                                                       *
  79. 1780 '*      Repeat the following lines for each observation:                 *
  80. 1790 '*           fobs(i): observed function value                            *
  81. 1800 '*           if iuserwt = 1 include obswt(i): function weight, point i   *
  82. 1810 '*           vobs(i,1): observed value, variable 1.                      *
  83. 1820 '*           if iuserwt = 1 include varwt(i,1): variable weight          *
  84. 1830 '*           vobs(i,2): observed value, variable 2.                      *
  85. 1840 '*           if iuserwt = 1 include varwt(i,2): variable weight          *
  86. 1850 '*           ...vobs(i,nvar): observed value, variable nvar              *
  87. 1860 '*              if iuserwt = 1 include varwt(i,nvar): variable weight    *
  88. 1870 '*                                                                       *
  89. 1880 '*      Repeat the following lines for each parameter:                   *
  90. 1890 '*           pname$(i): the name of the parameter (length <= 8)          *
  91. 1900 '*           p(i): initial guess at the parameter's value                *
  92. 1910 '*                                                                       *
  93. 1911 '*************************************************************************
  94. 1912 '* Note to IPCO users:                                                   *
  95. 1913 '* Included on this disk are the following files for testing and         *
  96. 1914 '* demonstration purposes:                                               *
  97. 1915 '*     DATA       set of test data to be fit to the following functions: *
  98. 1916 '*     FUNC1.BAS  gaussian function                                      *
  99. 1917 '*     FUNC2.BAS  lorentzian function                                    *
  100. 1920 '*************************************************************************
  101. 1921 '*                                                                       *
  102. 1922 '* Note: ocasionally the program diverges, i.e. the fit of the calculated*
  103. 1923 '* function gets worse each iteration, rather than better. If this       *
  104. 1924 '* occurs, it indicates one of two things:                               *
  105. 1925 '*    1. Your initial guesses for the parameters were too far off.       *
  106. 1926 '*       Make a better guess, and try again.                             *
  107. 1927 '*    2. The function you are using really can't adequately describe     *
  108. 1928 '*       your data.                                                      *
  109. 1929 '*************************************************************************
  110. 1930 '* DECLARATIONS:                                                         *
  111. 1940 '* NOBS: INT                                      NUMBER OF OBSERVATIONS *
  112. 1950 '* NVAR: INT                                NUMBER OF VARIABLES OBSERVED *
  113. 1960 '* FOBS: ARRAY(1..NOBS) OF DOUBLE               OBSERVED FUNCTION VALUES *
  114. 1970 '* FCALC:ARRAY(1..NOBS) OF DBL                CALCULATED FUNCTION VALUES *
  115. 1980 '* FTEMP:ARRAY(1..NOBS) OF DBL            HOLDS A SET OF FUNCTION VALUES *
  116. 1990 '*                                         DURING DERIVATIVE CALCULATION *
  117. 2000 '* VOBS: ARRAY(1..NOBS,1..NVAR) OF DBL          OBSERVED VARIABLE VALUES *
  118. 2010 '* NP:   INT                                        NUMBER OF PARAMETERS *
  119. 2020 '* P:    ARRAY(1..NP) OF DBL                            PARAMETER VALUES *
  120. 2030 '* PNAME$: ARRAY(1..NP) OF STRING                        PARAMETER NAMES *
  121. 2040 '* DFDP:ARRAY(1..NOBS,1..NP) OF DBL   PARTIALS OF FUNC W.R.T. PARAMETERS *
  122. 2050 '* DFDV:ARRAY(1..NOBS,1..NVAR) OF DBL     "     "   "     "   VARIABLES  *
  123. 2060 '* OBSWT:ARRAY(1..NOBS) OF DBL             OBSERVATION WEIGHTING FACTORS *
  124. 2070 '* VARWT:ARRAY(1..NOBS,1..NVAR) OF DBL        VARIABLE WEIGHTING FACTORS *
  125. 2080 '* DLAMBDA:ARRAY(1..NOBS) OF DBL                  LAGRANGIAN MULTIPLIERS *
  126. 2090 '* B:   ARRAY(1..NP,1..NP) OF DBL        COEFFICIENTS IN MATRIX EQUATION *
  127. 2100 '* RHS: ARRAY(1..NP) OF DBL           RIGHT HAND SIDE OF MATRIX EQUATION *
  128. 2110 '*                     NOTE: AFTER SOLUTION OF EQN, RHS HOLDS CHANGES TO *
  129. 2120 '*                               PARAMETERS, AND B HOLDS ITS OWN INVERSE *
  130. 2130 '* FRACT: DBL       FRACTION OF CALCULATED CHANGE TO APPLY TO PARAMETERS *
  131. 2140 '* DEVSQ: DBL       SUM OF WEIGHTED SQUARED DEVIATIONS OF CALC. FUNCTION *
  132. 2150 '* DEVSQ1:DBL                                DEVSQ VALUE, LAST ITERATION *
  133. 2160 '* DEVSQ2:DBL                           DEVSQ VALUE, TWO ITERATIONS BACK *
  134. 2170 '* IUSERWT: INT            IF IUSERWT=1, USER SUPPLIES WEIGHTING FACTORS *
  135. 2180 '*                         IF IUSERWT=0, NONLIN ASSUMES ALL WEIGHTS = 1  *
  136. 2190 '* INTERNALWT: INT             IF INTERNALWT=1 NONLIN CALCULATES WEIGHTS *
  137. 2200 '*                                 IF INTERNALWT=0 NO WEIGHT CALCULATION *
  138. 2210 '*************************************************************************
  139. 10000 OPTION BASE 1
  140. 10010 DEFINT I-N
  141. 10020 DEFDBL A-H,O-Z
  142. 10025 ON ERROR GOTO 65000
  143. 10030 CLS: LOCATE 10,30: INPUT "FUNCTION?  ",FUNCTION$: COMMON FUNCTION$
  144. 10040 CHAIN MERGE FUNCTION$,10050   'bring in function specific lines
  145. 10045 ON ERROR GOTO 0
  146. 10050 NVAR = 1: NP = 3
  147. 10060 GOSUB 18000   'initialization routine
  148. 10070 '
  149. 10080 FOR IT = 1 TO NUMIT
  150. 10090    'print progress report on screen
  151. 10100        GOSUB 12000   'subroutine iteration report
  152. 10110    'Test for non-convergance, exit if so
  153. 10120        GOSUB 13000      'subroutine converge
  154. 10130        IF NONCONVERGENCE = 1                                                              THEN PRINT"nonconvergence termination" : GOTO 10300
  155. 10140    'Calculate lagrangian multipliers
  156. 10150        GOSUB 13500    'subroutine lambda
  157. 10160    'If internal weighting desired, calculate new obswts
  158. 10170        IF INTERNALWT = 1                                                                  THEN GOSUB 14500
  159. 10180    'Set up matrix equation to get parameter changes
  160. 10190        GOSUB 15000  'SUBROUTINE SETUP
  161. 10200    'Solve equation for parameter changes
  162. 10210        GOSUB 16000  'subroutine solve
  163. 10220    'Apply changes
  164. 10230        GOSUB 17000  'subroutine deltap
  165. 10240    NEXT IT
  166. 10250 '
  167. 10260 'print final report
  168. 10270    GOSUB 19000      'subroutine report
  169. 10280 'Do any final processing (user supplied)
  170. 10290    GOSUB 20000
  171. 10300 END
  172. 11000 '********************************************************************
  173. 11010 '* SUBROUTINE LORENTZIAN                                            *
  174. 11020 '*                                                                  *
  175. 11030 '* Fills the array fcalc(nobs) with values along a lorentzian of    *
  176. 11040 '* the form:                                                        *
  177. 11050 '*                                                                  *
  178. 11060 '* lor(pos) = H / [(1/W)^2 * (pos - P)^2 + 1]                       *
  179. 11070 '*                                                                  *
  180. 11080 '* where H, W, and P are parameters to be fit.                      *
  181. 11090 '* Assignments: H = p(1)  W = p(2) P = p(3)    pos = vobs(i,1)      *
  182. 11100 '********************************************************************
  183. 11110 '
  184. 11130 WSQ2 = 1# / (P(2) * P(2))
  185. 11160 FOR I = 1 TO NOBS
  186. 11170    FCALC(I) =  P(1) / (WSQ2 * (VOBS(I,1) - P(3))^2 + 1#)
  187. 11230    NEXT I
  188. 11240 RETURN
  189. 12000 '***********************************************************************
  190. 12010 '* SUBROUTINE ITERTATION REPORT                                        *
  191. 12020 '*                                                                     *
  192. 12030 '* Prints out current parameters, function values, and deviation       *
  193. 12040 '***********************************************************************
  194. 12050 GOSUB 12200                   'print parameters
  195. 12060 GOSUB 11000                   'get function values in fcalc
  196. 12070 GOSUB 12500                   'print function values and deviation
  197. 12080 RETURN
  198. 12200 '*********************************************************************
  199. 12210 '* SUBROUTINE PARAMREPORT( IT, NP, P)                                *
  200. 12220 '* 2/7/82 by Dave Whitman                                            *
  201. 12230 '* Prints out current parameter values                               *
  202. 12240 '*********************************************************************
  203. 12250 CLS : LOCATE 4,4 : BEEP  'Wake up,Dave!
  204. 12260 PRINT "Parameters, Iteration";IT : PRINT
  205. 12270 LOCATE 7,2
  206. 12280 COLOR 1 : PRINT "  Name  │  Value   │ Change "; : COLOR 7 : PRINT
  207. 12290 FOR I = 1 TO NP
  208. 12300    LOCATE I+7,2
  209. 12310    PRINT USING  "\      \\\####.##\ \+##.##"; PNAME$(I);"│ "; P(I);                             "  │ ";-1 * RHS(I)
  210. 12320    NEXT I
  211. 12330 PRINT
  212. 12340 RETURN
  213. 12500 '***********************************************************************
  214. 12510 '* SUBROUTINE FUNCTIONREPORT                                           *
  215. 12520 '* Prints obs. and calc. function values, and deviation between them   *
  216. 12530 '***********************************************************************
  217. 12540 IROW = 1 : ICOL = 40 : IOFFSET = 20 : IROOM = 40 : NUMROWS = 20
  218. 12550 LOCATE IROW,ICOL : COLOR 1
  219. 12560 PRINT "  obs. │ calc. "; : COLOR 7 : PRINT
  220. 12570 IF NOBS >= NUMROWS                                                                 THEN LOCATE IROW,ICOL+IOFFSET : COLOR 1:                                             PRINT "  obs. │ calc. "; : COLOR 7 : PRINT
  221. 12580  DEVSQ = 0
  222. 12590  FOR I = 1 TO NOBS
  223. 12600     IF I > IROOM THEN 12630
  224. 12610     IF I <= NUMROWS THEN LOCATE (IROW + I),ICOL                                                     ELSE LOCATE (IROW + I MOD NUMROWS),(ICOL + IOFFSET)
  225. 12620     PRINT USING "###.##\ \###.##"; FOBS(I);" │ ";FCALC(I);
  226. 12630     DEV = FOBS(I) - FCALC(I)
  227. 12640     DEVSQ = DEVSQ + DEV * DEV * OBSWT(I)
  228. 12650     NEXT I
  229. 12660 LOCATE 20,5 : PRINT USING "\         \#####.##";"Σ error² = ";DEVSQ
  230. 12670 IF IT = 1 THEN 12690  'no change in first iteration
  231. 12680 LOCATE 21,5: PRINT USING "\         \#####.##";"Change =   ";DEVSQ-DEVSQ1;
  232. 12690  RETURN
  233. 13000 '*********************************************************************
  234. 13010 '* SUBROUTINE CONVERGE ( ERRSQ, DEVSQ, DEVSQ1, DEVSQ2, NONCONVERGE ) *
  235. 13020 '*                                                                   *
  236. 13030 '* Compares squared deviation of calculated function from observed   *
  237. 13040 '* function with that obtained in the last 2 iterations. If the      *
  238. 13050 '* deviation got worse two iterations in a row, set nonconverge flag.*
  239. 13060 '*********************************************************************
  240. 13070  IF (DEVSQ > DEVSQ1 AND DEVSQ1 > DEVSQ2)                                            THEN NONCONVERGE = 1
  241. 13080  DEVSQ2 = DEVSQ1
  242. 13090  DEVSQ1 = DEVSQ
  243. 13100  RETURN
  244. 13500 '********************************************************************
  245. 13510 '* SUBROUTINE LAMBDA ( DLAMBDA, FCALC, FOBS, OBSWT)                 *
  246. 13520 '*                                                                  *
  247. 13530 '* Calculates lagrangian multipliers for setting up matrix equation *
  248. 13540 '********************************************************************
  249. 13550 FOR I = 1 TO NOBS
  250. 13560    DLAMBDA(I) = (FCALC(I) - FOBS(I)) * OBSWT(I)
  251. 13570    NEXT I
  252. 13580 RETURN
  253. 14000 '*************************************************************
  254. 14010 '* SUBROUTINE VSLOPE( V,DFDV,NVAR )                          *
  255. 14020 '*                                                           *
  256. 14030 '* Calculates the partials of the function w/ r.t. each      *
  257. 14040 '* of the variables at each of the observed points, and      *
  258. 14050 '* stores them in dfdv.                                      *
  259. 14060 '*************************************************************
  260. 14070 GOSUB 11000   'call function( p, v, nobs, nvar, np, fcalc)
  261. 14080 FOR IW = 1 TO NOBS
  262. 14090    FTEMP(IW) = FCALC(IW)
  263. 14100    NEXT IW
  264. 14110 FOR IW = 1 TO NVAR
  265. 14120    FOR JW = 1 TO NOBS
  266. 14130       IFLAG(JW) = 0
  267. 14140       IF ABS(VOBS(JW,IW)) < 1D-20                                                           THEN VOBS(JW,IW) = .0005# : IFLAG(JW) = 1                                       ELSE VOBS(JW,IW) = VOBS(JW,IW) * 1.0005#
  268. 14150 PRINT "modified var:" ;VOBS(JW,IW)
  269. 14160       NEXT JW
  270. 14170    GOSUB 11000    'call function(fcalc)
  271. 14180    FOR JW = 1 TO NOBS
  272. 14190       IF IFLAG(JW) = 1                                                                   THEN DFDV(JW,IW) = (FCALC(JW) - FTEMP(JW)) / .0005# :                                VOBS(JW,IW) = VOBS(JW,IW) - .0005#
  273. 14200       IF IFLAG(JW) <> 1                                                                  THEN DFDV(JW,IW) = (FCALC(JW)-FTEMP(JW)) / (.0005# * VOBS(JW,IW))                    : VOBS(JW,IW) = VOBS(JW,IW) / 1.0005#
  274. 14210       PRINT"dfdv(";JW;IW;")="; DFDV(JW,IW)
  275. 14220       NEXT JW
  276. 14230    NEXT IW
  277. 14240 RETURN
  278. 14500 '********************************************************************
  279. 14510 '* subroutine weigh[ p, nobs, nvar, v, dfdv, varwt, obswt ]         *
  280. 14520 '*                                                                  *
  281. 14530 '* calculates new weights for function values,                      *
  282. 14540 '* using the follwing formula:                                      *
  283. 14550 '* obswt(i) = 1/ Σ [(dfdv)^2 * (1/varwt(v))]                        *
  284. 14560 '* note: obswt(i) = 1/L(i) in Wentworth article                     *
  285. 14570 '********************************************************************
  286. 14580 IF IT = 1 THEN RETURN  'skip weighting on first iteration
  287. 14590 GOSUB 14000 'subroutine vslope
  288. 14600 FOR IW = 1 TO NOBS
  289. 14610    SUM = 0#
  290. 14620    FOR JW = 1 TO NVAR
  291. 14630       SUM = SUM + DFDV(IW,JW)*DFDV(IW,JW)/VARWT(IW,JW)
  292. 14640       NEXT JW
  293. 14650    OBSWT(IW) = 1# / SUM
  294. 14660    NEXT IW
  295. 14670 PRINT "new function weights:"
  296. 14680 FOR IW = 1 TO NOBS
  297. 14690    PRINT OBSWT(IW)
  298. 14700    NEXT IW
  299. 14710 RETURN
  300. 15000 '********************************************************
  301. 15010 '* SUBROUTINE SETUP( B,RHS,dfdp)                        *
  302. 15020 '* Sets up matrix equation to get changes to parameters *
  303. 15030 '********************************************************
  304. 15040 '
  305. 15050 'Get partials of function w.r.t. parameters
  306. 15060     GOSUB 17500     'subroutine pslope
  307. 15070 'Now set up matrices
  308. 15080     FOR I = 1 TO NP
  309. 15090        'Set up right hand side element
  310. 15100        RHS(I) = 0#
  311. 15110        FOR J = 1 TO NOBS
  312. 15120           RHS(I) = RHS(I) + DFDP(J,I) * DLAMBDA(J)
  313. 15130           NEXT J
  314. 15140        'Set up left hand side elements
  315. 15150        FOR J = 1 TO NP
  316. 15160           B(I,J) = 0#
  317. 15170           FOR K = 1 TO NOBS
  318. 15180              B(I,J) = B(I,J) + DFDP(K,I) * DFDP(K,J) * OBSWT(K)
  319. 15190              NEXT K
  320. 15200           NEXT J
  321. 15210        NEXT I
  322. 15220 RETURN
  323. 16000 '*****************************************************
  324. 16010 '* subroutine solve[b#(np,np), rhs(np), np]          *
  325. 16020 '* 1/31/82 by Dave Whitman                           *
  326. 16030 '* solves matrix equations of the form b# x = rhs#   *
  327. 16040 '* inverts b# in place,multiplies rhs# by inverse    *
  328. 16050 '* uses Gauss-Jordan matrix inversion                *
  329. 16060 '* for good results b# and rhs# must be dbl precision*
  330. 16070 '* ref: J.M. McCormick "Numerical Methods in FORTRAN"*
  331. 16080 '*****************************************************
  332. 16090 DETERM = 1#
  333. 16100 FOR I = 1 TO NP
  334. 16110    INDEX(I,3) = 0
  335. 16120    NEXT I
  336. 16130 FOR I = 1 TO NP  'MAIN LOOP
  337. 16140    'search for pivot element
  338. 16150    MAX# = 0#
  339. 16160    FOR J = 1 TO NP
  340. 16170       IF INDEX(J,3) = 1 THEN 16260
  341. 16180       FOR K = 1 TO NP
  342. 16190          IF INDEX(K,3) > 1 THEN 16700
  343. 16200          IF INDEX(K,3) = 1 THEN 16250
  344. 16210          IF MAX# > ABS(B(J,K)) THEN 16250
  345. 16220          IROW = J
  346. 16230          ICOL = K
  347. 16240          MAX# = ABS(B(J,K))
  348. 16250          NEXT K
  349. 16260       NEXT J
  350. 16270   INDEX(ICOL,3) = INDEX(ICOL,3) + 1
  351. 16280   INDEX(I,1) = IROW
  352. 16290   INDEX(I,2) = ICOL
  353. 16300   'interchange rows to put pivot on diagonal
  354. 16310   IF IROW = ICOL THEN 16380  'ALREADY THERE
  355. 16320   DETERM = -1# * DETERM
  356. 16330   FOR J = 1 TO NP
  357. 16340      SWAP B(IROW,J),B(ICOL,J)
  358. 16350      NEXT J
  359. 16360   SWAP RHS(IROW),RHS(ICOL)
  360. 16370   'divide pivot row by pivot element
  361. 16380   PIVOT = B(ICOL,ICOL)
  362. 16390   DETERM = DETERM * PIVOT
  363. 16400   B(ICOL,ICOL) = 1#
  364. 16410   FOR J = 1 TO NP
  365. 16420      B(ICOL,J) = B(ICOL,J)/PIVOT
  366. 16430      NEXT J
  367. 16440   RHS(ICOL) = RHS(ICOL)/PIVOT
  368. 16450   ' reduce non-pivot rows
  369. 16460   FOR J = 1 TO NP
  370. 16470      IF J = ICOL THEN 16540
  371. 16480      T = B(J,ICOL)
  372. 16490      B(J,ICOL) = 0#
  373. 16500      FOR K = 1 TO NP
  374. 16510         B(J,K) = B(J,K) - B(ICOL,K)*T
  375. 16520         NEXT K
  376. 16530      RHS(J) = RHS(J) - RHS(ICOL)*T
  377. 16540      NEXT J
  378. 16550   NEXT I
  379. 16560 'interchange columns
  380. 16570 FOR I = NP TO 1 STEP -1
  381. 16580    IF INDEX(I,1) = INDEX(I,2) THEN 16640
  382. 16590    IROW = INDEX(I,1)
  383. 16600    ICOL = INDEX(I,2)
  384. 16610    FOR J = 1 TO NP
  385. 16620       SWAP B(J,IROW), B(J,ICOL)
  386. 16630       NEXT J
  387. 16640    NEXT I
  388. 16650 'test for singularity
  389. 16660 FOR I = 1 TO NP
  390. 16670    IF INDEX(I,3) <> 1 THEN 16700
  391. 16680    NEXT I
  392. 16690 RETURN
  393. 16700 PRINT"singular matrix error ":RETURN
  394. 17000 '**********************************************************
  395. 17010 '* SUBROUTINE DELTAP ( P, RHS, NP )                       *
  396. 17020 '*                                                        *
  397. 17030 '* Modifies parameters according to changes in rhs        *
  398. 17040 '**********************************************************
  399. 17050 FOR I = 1 TO NP
  400. 17060    P(I) = P(I) - RHS(I) * FRACT
  401. 17070    NEXT I
  402. 17080 RETURN
  403. 17500 '************************************************************************
  404. 17510 '* subroutine pslope[ p, vobs, nobs, np, nvar, dfdp(nobs,np) ]          *
  405. 17520 '* 2/1/82   by Dave Whitman                                             *
  406. 17530 '* calculates partial of the function with                              *
  407. 17540 '* respect to the each of the parameters at                             *
  408. 17550 '* each of the observations, and stores them in dfdp.                   *
  409. 17560 '************************************************************************
  410. 17570 GOSUB 11000                   'call function( fcalc, p, v, nobs, nvar, np)
  411. 17580 FOR IS = 1 TO NOBS
  412. 17590    FTEMP(IS) = FCALC(IS)
  413. 17600    NEXT IS
  414. 17610 FOR IS = 1 TO NP
  415. 17620    TP = P(IS)
  416. 17630    IF TP < 1D-20                                                                      THEN P(IS) = TP + .0005#                                                        ELSE P(IS) = TP * 1.0005#
  417. 17640    GOSUB 11000          'call function( fcalc )
  418. 17650    FOR JS = 1 TO NOBS
  419. 17660       IF TP < 1D-20                                                                      THEN DFDP(JS,IS) = (FCALC(JS) - FTEMP(JS)) / .0005#                             ELSE DFDP(JS,IS) = (FCALC(JS) - FTEMP(JS)) / (.0005# * TP)
  420. 17670       NEXT JS
  421. 17680    P(IS) = TP
  422. 17690    NEXT IS
  423. 17700 RETURN
  424. 18000 '**********************************************************************
  425. 18010 '* SUBROUTINE initialize                                              *
  426. 18020 '*                                                                    *
  427. 18030 '* Prompts for name of input file, then reads problem in              *
  428. 18040 '* from input file.                                                   *
  429. 18050 '**********************************************************************
  430. 18060 '
  431. 18070 KEY OFF : CLS
  432. 18080 LOCATE 13,20 : INPUT "Name of input file?  ",IFILE$
  433. 18090 OPEN IFILE$ FOR INPUT AS #1
  434. 18100 INPUT#1,NOBS
  435. 18110 DIM FOBS(NOBS),FCALC(NOBS),FTEMP(NOBS),OBSWT(NOBS)
  436. 18120 DIM VOBS(NOBS,NVAR),V(NOBS,NVAR),VARWT(NOBS,NVAR),DFDV(NOBS,NVAR)
  437. 18130 DIM P(NP),PNAME$(NP), DFDP(NOBS,NP)
  438. 18140 DIM IFLAG(NOBS),DLAMBDA(NOBS)
  439. 18150 INPUT#1,NUMIT
  440. 18160 INPUT#1,IUSERWT,INTERNALWT
  441. 18170 INPUT#1, FRACT    'fraction of calculated param. change to apply
  442. 18180 FOR I = 1 TO NOBS
  443. 18190    INPUT#1,FOBS(I)
  444. 18200    IF IUSERWT = 1                                                                     THEN INPUT#1,OBSWT(I)                                                           ELSE OBSWT(I) = 1#
  445. 18210    FOR J = 1 TO NVAR
  446. 18220       INPUT#1,VOBS(I,J)
  447. 18230       IF IUSERWT = 1                                                                     THEN INPUT#1,VARWT(I,J)                                                         ELSE VARWT(I,J) = 1#
  448. 18240       NEXT J
  449. 18250    NEXT I
  450. 18260 FOR I = 1 TO NP
  451. 18270    INPUT#1,PNAME$(I), P(I)
  452. 18280    NEXT I
  453. 18290 DEVSQ1 = 1D+20 : DEVSQ2 = 1D+20 'dummy devsqs for non-converge test
  454. 18300 TIME$ = "00:00:00"
  455. 18310 RETURN
  456. 19000 '**********************************************************************
  457. 19010 '* SUBROUTINE REPORT                                                  *
  458. 19020 '*                                                                    *
  459. 19030 '* Prints final report giving observed and calculated values of       *
  460. 19040 '* function, and standard deviation                                   *
  461. 19050 '* Note: assumes NEC 8023 printer. Modify to suit other printers.     *
  462. 19060 '* On NEC, esc X turns on underlining, esc Y turns it off.            *
  463. 19070 '* The following is a partial list of IBM screen charactors, followed *
  464. 19080 '* by the charactor the NEC will print: û = │ ; ┐ = Σ; ╠ = ².         *
  465. 19090 '* Thus "┐ error╠ =" prints as "Σ error² =".                          *
  466. 19100 '**********************************************************************
  467. 19110 GOSUB 11000  'subroutine function
  468. 19130 LPRINT "FINAL REPORT"
  469. 19140 LPRINT "FUNCTION: "; FUNCTION$
  470. 19150 LPRINT "DATA FILE:";IFILE$
  471. 19160 LPRINT : LPRINT"   Function Values"
  472. 19170 LPRINT CHR$(27);"X" "Observed û Calculated"; : LPRINT CHR$(27);"Y"
  473. 19180 DEVSQ = 0#
  474. 19190 FOR I = 1 TO NOBS
  475. 19200     LPRINT USING"####.##\   \####.##";FOBS(I);"  û  ";FCALC(I)
  476. 19210     DEVSQ = DEVSQ + (FOBS(I) - FCALC(I))^2 * OBSWT(I)
  477. 19220     NEXT I
  478. 19230 LPRINT
  479. 19240 LPRINT USING "\        \####.##\                          \+#.##";                  "┐ error╠ = "; DEVSQ; "   Change, last iteration = "; DEVSQ - DEVSQ1
  480. 19250 GOSUB 19500 'subroutine covariance
  481. 19260 LPRINT: LPRINT"FINAL PARAMETERS"
  482. 19270 LPRINT CHR$(27);"X" " Name    û  Value  û Std.Dev."; : LPRINT CHR$(27);"Y"
  483. 19280 FOR I = 1 TO NP
  484. 19290   LPRINT USING  "\       \\\####.##\ \##.###"; PNAME$(I);"û "; P(I);" û ";           SQR(B(I,I))
  485. 19300   NEXT I
  486. 19310 LPRINT
  487. 19320 KEY ON: RETURN
  488. 19500 '**********************************************************************
  489. 19510 '* SUBROUTINE COVARIANCE                                              *
  490. 19520 '*                                                                    *
  491. 19530 '* Calculates estimate of unit variance, then calculates              *
  492. 19540 '* variance/covariance matrix based on inverted B matrix and the      *
  493. 19550 '* variance estimate.                                                 *
  494. 19560 '**********************************************************************
  495. 19570 'estimate unit variance
  496. 19580  IF NOBS <= NP                                                                       THEN VAR = 0#  'should never trust parameters if nobs < np anyways!             ELSE VAR = DEVSQ / (NOBS - NP)
  497. 19590 'convert B to covariance matrix
  498. 19600 FOR I = 1 TO NP
  499. 19610   FOR J = 1 TO NP
  500. 19620      B(I,J) = B(I,J) * VAR
  501. 19630      NEXT J
  502. 19640   NEXT I
  503. 19650 RETURN
  504. 20000 '**********************************************************************
  505. 20010 '* subroutine finalproc                                               *
  506. 20020 '*                                                                    *
  507. 20030 '* Before nonlin stops, it makes a call to a subroutine at line 20000 *
  508. 20040 '* The user may supply a subroutine (in the file with the function    *
  509. 20050 '* subroutine) to do any final calculations using either the final    *
  510. 20060 '* parameter set and/or the variance-covariance matrix.               *
  511. 20070 '**********************************************************************
  512. 20080 RETURN
  513. 65000 ' Trap error of function file not in ascii mode
  514. 65010 IF ERR <> 54 THEN 65090
  515. 65020    CLS: BEEP : LOCATE 5,28
  516. 65030    PRINT "Bad File Mode Error:"
  517. 65040    LOCATE 7,21: PRINT "Function file must be saved in ASCII mode"
  518. 65050    LOCATE 8,15
  519. 65060    PRINT "Read lines 1200-1260 of this program for clarification."
  520. 65070    LOCATE 15,1: LIST 1200-1260
  521. 65080    LOCATE 23,1: STOP
  522. 65090 RESUME
  523. 1200-1260 of this program for clarification."
  524. 65070    LOCATE 15,1: LIST 1200-1260
  525. 650