To extract into actual directory: tar -xzf aplcon0.tgz Compile test program: g77 taplcon0.F condfit.F condutil.F or: gfortran taplcon0.F condfit.F condutil.F and then run a.out APLCON manual =============== 1. Introduction --------------- APLCON is a fit program based on "constraint equations". Variables of two types are distinguished: o measured variables, with a measured value, stored by the user in the variable array X(.), and an uncertainty, expressed by a variance and optionally covariances to other parameters. These elements of the covariance matrix of the measured variables have to be stored by the user in the matrix-array V(.) before the fit; o unmeasured variables, or fit parameters, which are not directly measured; through constraint equations (with sufficient rank) they become defined. By convention they are treated like the other measured) variables and the user has to store a reasonable approximate value in the variable-array X(.). By convention their variance in the matrix-array V(.) has to be zero. The variables (measured and unmeasured) are in array X(NVAR) for all NVAR variables. The elements of the covariance matrix of the measured variables are in the matrix-array V(.). The user has to provide statements of the form F(i) = function of X(.) ... for the equality constraints to calculate the values F(i) for the current values of variables X(.). These statements are evaluated repeatedly in a loop during the fit. Convergence is reached for sufficiently small values of the equality constraint values F(i). All arrays X(.) = variables (measured and unmeasured), V(.) = covarinace matrix elements, and F(.) = values of constraint equations are in double precision. A high accuracy, with all calculations in double precision, is necessary for the calculation of the equality constraint values F(i); this accuracy is required, because the derivative calculation during the fit determines first derivatives of the constraint equations w.r.t. to the variables from finite differences. The matrix elements are stored in symmetric storage mode, according to the schema j= 1 2 3 +--------... 1| 1 2| 2 3 index = IJSYM(I,J) i= 3| 4 5 6 4| 7 ... For example the element with index pair (I,J) = (3,2) is stored in the array element with index 5. The utility function IJSYM(I,J) can be used to obtaind the index: IJSYM(3,2) and IJSYM(2,3) will return value 5. The total size of the matrix array for a NVAR-by-NVAR symmetric matrix can be calculated by IJSYM(NVAR,NVAR), which is (NVAR*NVAR+NVAR)/2. In all calculation with the symmetric matrix the symmetry property is used to reduce the amount of computation. At convergence the fitted variables values for measured and unmeasured variables are stored in the variable-array X(.), and the corresponding covariance matrix of all variables is stored in the matrix-array V(.). 2. Initialization ----------------- The dimension parameters are required to initialize the program. NVAR = number of variables NEQS = number of constraint equations CALL APLCON(NVAR,NEQS) Optionally the user can set names (up to 16 chars) for the variables (for up to 128 variables): to set the variable name for X(index) CALL APNAME(index, text for name) ... 3. Declaration and definition of user arrays -------------------------------------------- There are three user arrays: X(NVAR) = array of all (measured and unmeasured) variables V(NSYM) = array for symmetric NVAR-by-NVAR matrix with dimension NSYM F(NEQS) = array for current values of constraint equations The three arrays have to be declared as double precision arrays of sufficient dimension. For a measured variable with index the element X(index) has to be set to the measured value, and the corresponding elements of the matrix array V(...) have be be set to the covariabe matrix elements. For unmeasured variables with index the element X(index) has to be set to an approximate value (used as start value in the fit), and the corresponding elements of the matrix array V(...) have to be set to zero. Example for initialization, and declaration and definition of arrays: PARAMETER (NVAR=3, NEQS=1) DOUBLE PRECISION X(NVAR),V((NAVR*NVAR+NVAR)/2),F(NEQS) DO I=1,IJSYM(NVAR,NVAR) V(I)=0.0D0 END DO X(1)= ... (measured value) X(2)= ... (measured value) X(3)= ... (unmeasured, approximate value) V(IJSYM(1,1))= ... V(IJSYM(1,2))= ... V(IJSYM(2,2))= ... CALL APLCON(NVAR,NEQS) 4. Fit loop with evaluation of constraints The fit requires the repeated evaluation of all constraint equations F(1) ... F(NEQS). After the calculation of the values of all equations with the current values of the variables X(1) ... X(NVAR) the subr. APLOOP has to be called with the arrays X, V and F. If the returned value of the flag IREP is < 0, the calculation has to be repeated. 10 F(1)= ... (function of X(1) ... X(NVAR)) ... CALL APLOOP(X,V,X, IREP) IF(IREP.LT.0) GOTO 10 Within the loop all necessary derivatives of the constraint equations are determined numerically; when all derivatives are ready, a step in parameter space is determined by matrix operations. This sequence is repeated until convergence. After convergence (with value IREP= ) the array X(.) contains the final fitted variable values and the array V(.) the corresponding covariance matrix elements. IREP < 0 continue IREP = 0 convergence IREP = 1 no convergence: IREP = 2 no convergence: too many iterations IREP = 3 no convergence: unphysical parameter values IREP = 4 no convergence: < 0 degrees of freedom IREP = 5 no convergence: insufficient storage 5. Fixed variables and variable transformations Variables, i.e. elements of the array X(.), can be treated as fixed by the call CALL APSTEP(I,0.0D0) ! treat V(I) as fixed By default variables are treated in the least squares sense with variances and covariances given by the elements of the matrix array VX(.). Variables with variance element (diagonal element V(i,i) in the matrix array VX(.)) equal to zero are treated as unmeasured variables. Variables can also be treated as o Poisson-distributed variables: counts (integers) follow the Poisson distribution; the variance is equal to the expected count and this is taken into account by a dynamic change if the variance during the iterative solution. o Log-normal distributed variables: variables with uncertainties given as a relative number (e.g. percentage) follow a log-normal distribution. The original variable is transformed to the log, and the log, which should follow the normal distrinbution, is treated in the fit. CALL APOISS(I) ! Poisson distributed variable CALL APLOGN(I) ! Lognormal distributed variable CALL APSQRT(I) ! sqrt transformed variable 6. Optional calls Between the APLCON-call and the loop the user may change certain default parameters by calls. CALL APRINT(LUNP,IPR) ! set print option CALL APDEPS(EPSIL) ! constraint accuracy CALL APITER(ITEMAX) ! max number of iterations CALL APDERF(DERFAC) ! factor for stepsize (measured) CALL APDERU(DERFAC) ! factor for stepsize (unmeasured) CALL APDLOW(DERFAC) ! factor for lower stepsize limit CALL APSTEP(IA,STEP) ! step size for numerical diff. [ CALL APLIMT(IA,XLOW,XHIG) ! range of variable X(IA) ] LUNP = 6 default print unit JDEBUG = 0 no printout = 1 minimum printout = 2 = 3 = 4 EPSIL = 1.0D-6 default value ITEMAX = 100 default value DERFAC = 1.0D-3 default value for stepsize (measured) = 1.0D-5 default value stepsize (unmeasured) = 1.0D-2 default value for lower stepsize limit STEP(i) = 1.0D-3 * SIGMA(i) for measured X(i) = 1.0D-5 * MAX(1, |X(i)|) initially for unmeasured X(i) =< 1.0D-2 * |X(i)| 7. Printout The vector X(.) of variables and the covariance matrix VX(.) are printed (on unit 6) by CALL CFPRV(6,X,VX,NX) 8. Variable reduction Variables from a vector X(.) and a covariance matrix VX(.) can be selected or reordered, and the selected or reordered variables can be stored in a vector Y(.) and a covariance matrix VY(.) CALL SIMSEL(X,VX,NY,LIST,Y,VY) where LIST(.) - integer array of indices, to be selected. For example if X(.) variables 1,3 and 7 have to be selected, the code would be LIST(1)=1 LIST(2)=3 LIST(3)=7 NY=3 CALL SIMSEL(X,VX,NY,LIST,Y,VY) CALL CFPRV(6,Y,VY,NY) 9. Normalization factor A normalization factor 1 +- epsilon is extracted from a vector X(.) and a covariance matrix VX(.) by SIMTR, and added as a N+1-st variable. CALL SIMTRN(X,VX,NX) ! extract normalization factor CALL CFPRV(6,X,VX,NX+1) ! print result (on unit 6) 10. Propagation of uncertainties The package can also be used for the propagation of uncertainties. It is assumed that variables X(.) with covariance matrix VX(.) exist, and a transformation to new variables Y(.) is done. The covariance matrix VY(.) of the transformed variables Y(.) = Y(X) can be calculated using the first derivatives of Y(.) w.r.t. X(.). This calculation is done with the following code: 10 Y(1)= ... (function of X(1) ... X(NVAR)) ... CALL APROPA(X,VX,Y,VY, IREP) IF(IREP.LT.0) GOTO 10 Appendix A: List of subprograms SUBROUTINE APLCON(NVAR,MCST) ! dimension parameters SUBROUTINE APLOOP(X,VX,F, IRET) ! steering routine for loop SUBROUTINE APROPA(X,VX,Y,VY, IRET) ! error propagation Y = Y(X) SUBROUTINE IPLMAT(X,VX,A,ST) ! check derivative matrix SUBROUTINE IPLDER(X,F, A,ST,XL,FC,HH, JRET) ! derivative calculation SUBROUTINE IPLCON(X,VX,F, AUX,IRET) ! internal steering routine SUBROUTINE JPLCON(X,VX,F, AUX,XS,DX,XP,R,W,DIAG,IRET) ! internal Appendix B: Print flag The default value of the print flag is IPR=5 IPR=0 =1 =2 =3 =4 =5 =6 =7 _________________________________ Program title x x x x x Program parameters x x x x Stop conditions x x x x x x x Error conditions x x x x x x Program exit x x x x x Iterations x x x x Final vector + err. x x x x Correlation matrix x x x Derivative matrix x x Debug x