MINGA

Version 2.05

User’s Manual

Third edition

by

Tomasz Plewa

Warsaw University Observatory Al. Ujazdowskie 4 00-478 Warsaw Poland

August 1993

The following document describes the use and the structure of the MINGA program. The MINGA has been created during years 1988–91, and serves in global analysis of the light curves and/or radial velocities of the double star systems using the Wilson–Devinney method of synthesis. However, the MINGA is a general optimizer and may be used to solve any other optimization problem. The computer code is available on request from the author.

Special thanks to: Ryszard Szczerba, Marco Lolli, and Janusz Kaluzny.

Comments on the program, versions with private modifications, and publications referenced to the MINGA program, please send to the author.

History and Philosophy

Primarily, MINGA has been written as the two separate packages: the SDC for global or preliminary minimization, and the MINW as the high quality final (local) minimizer. These two programs were used for solution of HT Vir [10]. The next step was the GRID/MINW package, where as the global minimizer serves the Levenberg–Marquardt algorithm. This version has been successfully applied to the analysis of BC Gru [11], and CT Tau [12]. All analyses concern the contact binaries. Recently, we used MINGA to analyze the UV extinction curves [17]. So that, MINGA seems to be a universal and simple in use.

Our motivation was to create the user–friendly program, which will offer the various ways to achieve the global solution quickly. This is the main advantage of MINGA over the other programs. Two of them, namely, the DC [19] and the CRS [3] codes, are slow, however, the CRS method is able to find the global minimum.

Besides pure (classical) local (or global) methods, as the DC and the CRS, there is still another approach: the hybrid (local and global) schemes. The MINGA belongs to them. Recently, the CRS algorithm, kindly provided by Dr. Guido Russo, has been included into the MINGA. Clustering analysis allowed to merge the MINGA minimization procedures with the Price code of F. Barone, L. Milano, and G. Russo [2]. Performed tests shows, that the MINGA is much faster during global optimization.

1 Components

The MINGA program has been written and developed on the base of three programs:

• MINUIT is a package of procedures used for minimization of a function of many variables. An original source is the CERNLIB library; for description see [7].

• LCE is a program generating the light curves (LC) and/or the radial velocities (RV) of the double star systems (DSS), known as the Wilson–Devinney code (WD). Description of the newest version of this code has been published by R.E. Wilson [18]. The LCE is a basic MINGA function.

• PRICE is a global minimization program developed by F. Barone, L. Milano, and G. Russo which utilizes the algorithm of Price [14]; for description and working examples see [2].

Description presented below is a synthetic compilation of manuals mentioned above, with changes and extensions appropriate for adoption of programs for common use.

1.1 MINGA Version of MINUIT

Solution of the problem of fitting LC/RV DSS data is defined by

min ; p⃗ χ^{2} = (O⃗(p⃗) − C⃗)^{T} V_{obs}(O⃗(p⃗) − C⃗)

where O⃗ –the observed quantities (LC/RV); p⃗ –vector of parameters; V_{obs} –covariance matrix of observations O⃗ (LC/RV). The types and the dimensions of each quantities are

quantity type  dimension
O⃗  vector N_{obs}
C⃗  vector N_{obs}
p⃗  vector N_{par}
V_{obs}  matrix N_{obs}×N_{obs}

where N_{obs} is the number of observations, and N_{par} is the number of parameters. If an information about the correlations between observations is not available, then the matrix V_{obs} is diagonal,

V_{obs} = \{ \begin{matrix}w_{i} & if i = j \\ 0 & otherwise \end{matrix} i, j = 1, N_{obs}

where w_{i}–the weight, formally equal to

w_{i} = \frac{1}{σ_{i}^{2}}

and σ_{i} is the dispersion of the i–th observation.

Changes in Relation to the Original Version of MINUIT

• High–precision inversion of matrices is performed by the procedures adopted from the Numerical Recipes library (LUDCMP, LUBKSB, MPROVE; for description see [13]).

• The internal transformation of variables by trigonometric functions has been removed.

• The GRID procedure serves in global minimization, i.e., a grid of solutions may be done.

Changes in Relation to the Previous Version of MINGA

• The alternative method of global optimization is the CRSshell procedure. It consists of the CRS method and the clustering analysis performed by the SHELL procedure. The selected clusters may be automatically analyzed using ANALYSE procedure which utilizes the MINGA local minimization methods.

Function and Compiler Dependencies

• Dimensioning and the function dependent parts of the MINGA are defined using the INCLUDE files (see Appendix A for more information). The function dependent

     INCLUDE files are

• LCE.INC if the function is LCE,

• NOT LCE.INC for any other function.

Carefully inspect MINGA.FOR and M1A.FOR files to be sure that lines marked by c!! or ? are consistent with the problem to be solved.

• The compiler dependent files are MINGA.FOR and L.FOR. Choose the appropriate compiler metacommands as they are marked in the headers of files.

Minimization Procedures

• SEEK performs the Monte Carlo minimization over the point p⃗_{0}. Points are chosen randomly, uniformly or Gaussian, over hypercube p⃗_{0}±ε(p⃗_{0}), where ε(p⃗_{0}) are current errors of the parameters at point p⃗_{0}.

• >SIMPLX performs a nongradient simplex minimization using algorithm due to Nelder and Mead [9]. We commend the use of the SIMPLX to you, when minimization is far from a minimum (i.e., during an initial phase of minimization).

• >MIGRAD is a variable–metric method (VMM). The unified Fletcher algorithm, switching method [6], is utilized. The MIGRAD is a very fast and effective method, not so far from the minimum. At the greater distance from minimum (where SIMPLX should be used), MIGRAD may have troubles with convergence. The use of this routine we recommend during the final minimization.

• IMPROV is a procedure of nonlocal minimization, which utilizes an information of a shape of a function in the local surroundings of p⃗_{0}, contained in the covariance matrix of parameters V(p⃗_{0}).

• ANNEAL is a procedure of global minimization by simulated annealing method [4]. We modified an original version of algorithm [8] to run with the automatic temperature adjustment.

• >CRSshell is a procedure of global minimization. For description of algorithm see an original paper of Price [14]. The Wilson–Price Code Handbook [2] gives an excellent introduction to the CRS method and general discussion of global minimization problems.

Cluster Analysis

We have introduced the SHELL procedure, which does the cluster analysis. The present version of the SHELL algorithm is rather simple, however, I have found that results (i.e. the locations of clusters) are in good agreement with that, what one can see using two dimensional projections of multidimensional sets of points. The most important thing for the cluster analysis is a proper definition of the metric. The unit value of particular variable should not be too large as the local minimization is \underline{not global}!.

Analysis of the Errors and Correlations

• MINOS is a basic procedure serves estimation of the errors of the parameters ε(p⃗_{0}) for any type of function. However, the MINOS is a very time consuming routine, and we do not recommend it, especially in the case of LCE.

• >HESSE is a procedure, which estimates the parabolic errors ε(p⃗_{0}) on the base of the matrix V(p⃗_{0}),

V (p⃗_{0}) = H^{−1}(p⃗_{0})

H(p⃗_{0}) = \frac{∂^{2}χ^{2}}{∂p⃗^{2}} p⃗_{0}

\underline{Note}

• Symbol > indicates procedures recommended solving the light curves fitting problem.

• If you plan to use MINGA for other function than LCE, then prepare function CHILC(P), where CHILC is the function value, P is a vector of parameters. Be sure that CHILC does \underline{not} uses the MINGA I/O units (see Section 2).

• For the more detailed description of all procedures see [7].

1.2 LCE

The following changes were introduced in relation to the original version of the WD code:

• All the previous versions of the WD code contain an error in procedures, which calculates an index of a vector MMSAVE for component 1. These procedures call the element MMSAVE(0), what causes an error array reference outside limits .

• In the case of an eccentric orbit, it is necessary to calculate a distribution of luminosity over a surface of both components at each photometric phase. On the other hand, in the case of circular orbit (constancy of geometry), such necessity rise only once (at the first photometric phase). LCE works in both cases in the same way, what causes the waste of a computer time. I did simple modifications in a logical structure of procedure BB2.

• Luminosities l_{1} and l_{2} are defined in interval [0,1] (in previous version of LCE [0,4π]).

• An effective wavelength λ is expressed in Ångströms (previously in microns).

• The mean effective temperatures are expressed in Kelvins (previously in 10^{4} K).

• The phases of the LC and RV are defined as

LC [0,1]
RV ]1,2] component 1
RV ]2,3] component 2

• In the case of contact configuration (modes 1 and 3, Ω_{1} = Ω_{2} , Ω–the nondimensional surface potential), it is convenient to use the fill–out parameter

f = \frac{Ω_{inner} − Ω}{Ω_{inner} − Ω_{outer}}

and

f = \{ \begin{matrix}0 & Ω = Ω_{inner} , \\ 1 & Ω = Ω_{outer} , \\ ]0, 1[ & Ω_{inner} > Ω > Ω_{outer} . \end{matrix}

In the case of marginal contact (f ≈ 0), one can use f < 0, e.g., f = −0.01. However, such solution needs a precise, subsequent qualification of configuration (C, SD1, SD2, D, double contact).

• Recent version of LCE offers an extended spot option, which permits a definition of regions (spots), which may be bright (hot), and/or dark (cool). The shape of spot may be circular or \underline{quadratic}. However, in the latter case, the limiting bounds of spot are of constant co–latitudes and longitudes. Definition of co–latitudes and longitudes is exactly the same as in the original WD code. The contrast of temperature of a spot, in relation to surrounding photosphere, is expressed in Kelvins,

dT \{ \begin{matrix}< 0 & cool (dark) spot, \\ = 0 & no spot, \\ > 0 & hot (bright) spot. \end{matrix}

If several spots overlaps each other, then dT is an arithmetic sum of contributions of each spot. The temperature contrast of the circular spot may be limb darkened,

dT (s) = dT (1 − X_{s} + X_{s} cos \frac{s}{R})

X_{s} is the limb darkening coefficient of dT, and s is the running distance from the spot center.

\underline{Note} It is recomended to calculate nonlinear (MINOS-like) errors using GRID option. A and V_{γ} parameters are independent of components’ luminosities, and only weakly depends on secondary’s temperature.

2 I/O Files

The structure of external I/O files used by MINGA is presented below. Each file may be of a type I (input), a type O (output) or a type I/O (input/output). Additional designations are

UD – a file name is defined by user, FIXED – a file name must not be changed.

• *.CMD (I/O=1, I, UD)

The input file containing the parameters of the model (LCE), and the commands controlling a flow of minimization process (see Appendix B for description of available commands). \underline{Note} The *.CMD files presented on distribution diskette on subdirectory \MINGA\EXAMP are examples of use of the MINGA commands.

• *.SOL (I/O=2, O, UD)

The output file containing a history of the minimization process and the results. This file is opened/closed every time before/after the use, for protection of results against damage (e.g., power fail). If a file of an identical name already exists, then results of current run will be added at the end of file. The frequency of output to this file is controlled by command PRINTOUT <density>.

• *.PUN (I/O=3, I/O, UD)

This file is used by procedure HESSE, and after use of commands PUNCH and

 COVARIANCE .

• *.* (I/O=4, O, UD)

This file contains the synthetic values of LC/RV and residuals corresponding to the observational data. This option is under control of the parameter ONE defined by the file MIN.IOS; see below option ONE=1. If a file of identical name already exists, then it will be overwritten.

• CON (I/O=5, I)

• CON (I/O=6, O)

• LCE.ERR (I/O=7, O, FIXED) This file contains the informations about internal errors, which occurred during operation of LCE.

• LCE.TRC (I/O=8, O, FIXED) This file contains a short information about minimization process: the function value, and the values of parameters. The frequency of output to this file is controlled by command TRACE <freq>.

• *.* (I/O=8, O, UD)

This unit is used by CRSshell procedure to save the search array. The file name is controlled by command ARRAY FILE <file name> and the format of file by

command ARRAY FORM <format>.

• MIN.IOS (I/O=10, I, FIXED) This is the basic MINGA file, which contains the parameters controlling the mode of operation of the MINGA, use of a monitor screen and sound, and defining the names of the files *.CMD, *.SOL, and *.PUN. This file name must not be changed in no circumstances.

• The case ONE=1 *.NAM (I/O=4, I, UD) This file contains pairs of the input and output file names (see below) in turn. The number of pairs is unlimited. *.* (I/O=4, O, UD) This file is a counterpart of the file *.CMD, but only the LCE part is used. *.* (I/O=8, O, UD) This is an output file, which contains the LC/RV synthetic data corresponding to the input data, as it has been given by the previous file.

• MIN.RST (I/O=11, I/O, FIXED) It is possible to restart MINGA automatically, with the results adopted from the previous run as input data, using command RESTART . Try to see how it works.

• *.* (I/O=12, I/O, UD)

This file is used by CRSshell routine as output and ANALYSE routine as input and output. The output from the ANALYSE routine depends on the mode of ANALYSE analysis. The file name is defined using command SEED FILE <file name>.

3 Codes of the LCE Parameters

An integration of MINUIT and LCE requires a proper transfer of data between both programs. Particularly, an objective function being minimized is defined by

χ^{2} = \frac{1}{ν} ∑_{i=1}^{N_{lc}} ∑_{j=1}^{N_{i}}(O_{ij} − C_{ij}(p⃗))^{2}w_{ij}

where N_{lc} –the total number of the LC/RV data; N_{i} –the number of observations for i–th LC or RV; ν –the number of degrees of freedom,

ν = ∑_{i=1}^{N_{lc}} N_{i} − N_{1} − N_{2}

where N_{1} is the number of parameters, which are independent on wavelength, and N_{2} is the number of parameters, which depends on the wavelength.

3.1 Parameters Independent on Wavelength

Parameters Which Do Not Describes Spots (0<code<200)

code symbol  code symbol
1 mode  2 iPB
3 if_{AT M1}  4 if_{AT M2}
5 n_{1}  6 n_{2}
7 P  8 ϑ_{e}
9 V_{unit}  10 ph_{N}
11 mode_{w}  12 if_{RAD}
13 e  14 λ_{per}
15 A  16 F_{1}
17 F_{2}  18 V_{γ}
19 ph_{SHIF T}  20 i
21 q  22 Ω_{1} or fill–out
23 Ω_{2}  24 T_{av1}
25 T_{av2}  26 g_{1}
27 g_{2}  28 A_{1}
29 A_{2}

Description of Spots (code<0)

code = −(10 000 komp + 1000 i + k)

komp 1,2
k  the number of spot
i  see below
i  shape
1 2 3 4 5 6
ϑ_{1} φ_{1} ϑ_{2} φ_{2} dT 0 quadratic
ϑ_{0} φ_{0} R X_{s} dT 1 circular

3.2 Parameters Dependent on Wavelength (code>200)

code = 200 + 10 (i − 1) + k i = 1, N_{lc}

k  1  2 3 4 5 6  7  8  9  10
symbol λ [Å] l_{1} l_{2} l_{3} X_{1} X_{2} ZERO  FACTOR N_{obs} WEIGHT

If the WEIGHT is less than 0, then the LCE interpolates synthetic intensities by means of the SPLINE method [13]. Such possibility is useful when you plan to use the true photometric observations, instead of the normal points. The approximate accuracy of the interpolation (the mean and maximum errors) is

ε_{mean} ÷ ε_{max}  configuration
n detached contact
20  0.0004 ÷ 0.004 0.0003 ÷ 0.002
30  0.0002 ÷ 0.003 0.0001 ÷ 0.001

4 Convergence Criteria

An estimation of the convergence is accomplished by use of the following quantities:

• EDM

If function could be represented by a multiparametric hyperparabola, then EDM is a distance to minimum (typically EDM is of order 10^{−5}),

EDM = \frac{1}{2}G⃗^{T}VG⃗

where G⃗ –the gradient; V –the covariance matrix.

• UP

is defined by command ERROR DEF <UP>, and is dependent on the number of parameters N_{par}, and the confidence level P. Typical value of UP suggested for the minimization, is 1.0. Before the analysis of the errors, e.g., before the use of the HESSE command, you should specify UP as below.

UP                    P[%]

N_{par}  1  50  70  90  95  99
1 1.6(-4) 0.45 1.07 2.71 3.84 6.64
2 2.0(-2) 1.39 2.41 4.61 5.99 9.21
3  0.12 2.37 3.67 6.25 7.81 11.35
4  0.30 3.36 4.88 7.78 9.49 13.28
5  0.55 4.35 6.06 9.24 11.07 15.09
6  0.87 5.35 7.23 10.64 12.59 16.81
7  1.24 6.35 8.38 12.02 14.07 18.48
8  1.65 7.34 9.52 13.36 15.51 20.09
9  2.09 8.34 10.66 14.68 16.92 21.67
10  2.56 9.34 11.78 15.99 18.31 23.21
11  3.05 10.34 12.90 17.28 19.68 24.73
12  3.57 11.34 14.01 18.55 21.03 26.22
13  4.11 12.34 15.12 19.81 22.36 27.69
14  4.66 13.34 16.22 21.06 23.69 29.14
15  5.23 14.34 17.32 22.31 25.00 30.58
σ level 0.013 0.67 1.03 1.65 1.96 2.58

• VTEST

defines a relative stability of the covariance matrix of parameters between consecutive minimization steps of the MIGRAD procedure. In other words, VTEST says us something about how fast the function shape changes between the following iterations (typically VTEST is about 0.01 ÷ 0.04),

VTEST = \sqrt{\frac{1}{N_{par}} ∑_{j=1}^{N_{par}} (\frac{V_{jj}^{next} − V_{jj}^{prev}}{V_{jj}^{next} + V_{jj}^{prev}} )^{2}}

V_{jj}^{next} is the j–th diagonal element of a matrix obtained in given iteration, and V_{jj}^{prev} is a similar element obtained in the previous step.

The user may specify his own convergence criteria using the command TOLERANCE <relχ^{2}> <relp⃗> <relG⃗> <EDM > <VTEST >, where the additional quantities are defined by

relχ^{2} = |\frac{χ^{2}_{next} − χ^{2}_{prev}}{max(χ^{2}_{prev}, χ^{2}_{typ}, ε)}|

relp⃗ = ∥\frac{p⃗_{next} − p⃗_{prev}}{typ p⃗} ∥_{∞}

relG⃗ = ∥\frac{typ p⃗ · G⃗}{max(χ^{2}_{prev}, χ^{2}_{typ}, ε)}∥_{∞}

where typ p⃗ = max(∥p⃗∥, ∥p⃗_{typ}∥, ε); p⃗_{typ} –the typical value of parameter; χ^{2}_{typ} –the typical value of function; ε –the machine epsilon. The use of such criteria is recommended in [5].

The MINGA states, that convergence has been achieved, if

                     if (relp⃗.and.EDM ) then
                         if (relχ^{2} .and.2VTEST ) STOP
                         if (relG⃗.and.VTEST ) STOP
                     end if

The MINGA may perform a local minimization as well as a global one (in sense of a search of minimum of χ^{2}(i,q); see, for example, [15] and [10]). Two sets of the convergence criteria are recommended,

Minimization

local global
relp⃗ 10^{−4} 10^{−2}
relχ^{2} 10^{−4} 10^{−3}
relG⃗ 1  10
VTEST 10^{−2} 0.2
EDM  10^{−4} 10^{−2}

5 Derivatives

At present, there are no sufficient examples showing an influence of the algorithm used for approximation of derivatives on the form of solution, or on the convergence of minimization process. Wilson and Biermann [20] found no big differences between the final solutions. However, because of the flatness of a region surroundings minimum and its considerable size [15, 10], one can state, that such situation is rather exceptional. On the other hand, the algorithm of variable–metric minimization of the MIGRAD procedure needs derivatives of high quality. Numerical experiments show, in the case of asymmetric derivatives, that the covariance matrix of parameters V (p⃗_{0}) often losses its positive definiteness. In connection with this, MINGA provides you with a choice of three different algorithms of numerical approximation of derivatives,

type algorithm  order
1 asymmetric O(h)
2 symmetric O(h^{2})
3 3–point rule {O(h^{3}) \atop see [1]}

The choice of a step size, used during calculation of derivatives, is based on a relative change of the function value. Permissible range of a function change is defined by command DF WANT <epsmin> <epsmax>. The step size h is modified, if

                if (relχ^{2}max > epsmax) h= \frac{1}{2}h
                if (relχ^{2}min < epsmin) h=2h
                if (relχ^{2}max > 25epsmax) h= \sqrt{\frac{epsmax}{relχ^{2}max}} h
                if (relχ^{2}max < 0.04epsmin) h= \sqrt{\frac{epsmin}{relχ^{2}max}} h

and

relχ^{2}max = \frac{max(χ^{2}_{+},χ^{2}_{−},χ^{2}_{3})}{χ^{2}_{0}}

relχ^{2}min = \frac{min(χ^{2}_{+},χ^{2}_{−},χ^{2}_{3})}{χ^{2}_{0}}

6 Speed, Accuracy and Compatibility

After optimization of some LCE routines (ATM, LIGHT, BB2), LCE is about 15 ÷ 30% faster, than the original WD code. At present, the use of spots has an negligible influence on speed of calculations (see procedures LUM, LIGHT and SPOT). The speed of calculations is determined mainly by the type of orbit (eccentric, circular), the number of parameters (speed decreases approximately as 1.5 N_{par}), and the grid fineness (speed decreases approximately as n^{2}, see Fig.1),

T ≈ 105.0 − 16.31 n + 1.207 n^{2} (−0.051 n^{3} + 0.001 n^{4}).

The accuracy of LCE depends critically on the grid fineness (see Fig. 2). In the case of an annular eclipse and marginal contact (f ≈ 0), numerical simulations^{1} shows the errors of total intensity (in units of 10^{−4}) and timing, as follows MINGA V2.05 12

1Versus n=60 on 10MHz PC/XT with 10MHz 8087, MS Fortran V4.10 compiler.

Figure 1: Timing of the LCE (SUBROUTINE BB2) versus grid fineness. Figure 2: Accuracy of the LCE versus grid fineness.

\begin{matrix}n = n_{1} = n_{2} & ε_{mean} & ε_{max} & faster× \\ 10 & 16.0 & 104.0 & 22.6 \\ 15 & 6.3 & 37.0 & 11.9 \\ 20 & 3.4 & 21.0 & 7.3 \\ 25 & 2.0 & 12.0 & 5.0 \\ 30 & 1.3 & 7.4 & 3.6 \\ 35 & 0.9 & 6.2 & 2.7 \\ 40 & 0.6 & 3.2 & 2.1 \\ 45 & 0.5 & 2.9 & 1.7 \\ 50 & 0.4 & 2.9 & 1.4 \\ 55 & 0.3 & 1.8 & 1.2 \end{matrix}

The MINGA is compatible with the Microsoft Fortran V5.00 and the SVS Fortran V2.8 compilers.

7 Dimensioning

All the matrices are dimensioned using the PARAMETER statements. The basic dimensions are defined using INCLUDE files as follows

• LCE

ng the grid fineness (FINENESS.INC),

nspomx the number of spots on each component (SPOTS.INC),

ndata the number of observations in various wavelengths (OBS.INC),

nlcmax the maximum number of various wavelengths (OBS.INC).

• MINGA

nn the maximum number of variables (NVARS.INC),

nsearch the maximum dimension of CRSshell search array (SEARCH.INC),

nmax the maximum dimension of the matrix which could be inverted. It is important during the error analysis, when HESSE calculates hessian matrix. Dimension of hessian is equal to the number of simultaneously variable parameters. If it is greater than nmax then the covariance matrix is lost (MATINV.INC).

\underline{Note} In majority of cases user is informed about an exceed of permissible dimension.

8 The MINGA Minimization Commands

The format of commands required by MINGA (file *.CMD) is

                            <name arguments>
                              (A10, 7F10.0)
                               (A10, 1X, A)

The second format is used only if the argument of command is of type CHARACTER (note the extra space after the command name). The name of command may be typed by capital or small letters, and the values of parameters are given as real numbers (even, if they have an integer value). The Appendix B gives a description of commands in alphabetical order in format

• command name and its arguments

• purpose

• description of the arguments Special designations

• Commands designated by the symbol > have been added to the original version of MINUIT for needs of MINGA,

• Commands designated by the symbol ◇ have been added to the previous version of MINGA.

Appendix A The MINGA V2.01 INCLUDE Files

Function Independent Files

MATINV.INC
c Definition of the maximum dimension of matrix to be inverted
        parameter (nmax= 10)
NVARS.INC
c Definition of the maximum number of parameters
        parameter (nn= 50)

SEARCH.INC
c Definition of the CRSshell parameters
        parameter (nsearch=300,nclmax=5,ndneighb=nsearch/10)
USER.INC
c Define Your name, phone number or anything else. It may be
c important for You or the other users too.
c Destroy '?' character.
        character*20 username
        parameter(username=?'                    ')
VERSION.INC
c Definition of the version
        character vnumber*7,vdate*10
        parameter (vnumber='2.05   ',vdate='1993/08/16')

Function Dependent Files

• NOT LCE.INC (list of routines only)

c+
c
c       File NOT_LCE.INC (1990/12/20)
c
c       The following routines are include for function other
c       than LCE
c
c       ERRLC
c       PACK1
c       NOTDTSP
c       UNIT
c-

• LCE files (only if the function is LCE)

LCE.INC
c 1991/02/21
        real*4 function unit(int)
c Definition of the LCE metric.
c                                Wavelength independent parameters
c                                Spot parameters
c                                Wavelength dependent parameters
FINENESS.INC
c Definition of the LCE grid fineness.
c   nmax      60   40   35   33  30  25  20
c   ng       3011 1345 1034 919 762 530 342
        parameter (nmax=35,ng=1035)
OBS.INC
c Definition of the maximum number of LC/RV data.
        parameter (nlcmax= 2,ndata= 450)
SPOTS.INC
c Definition of the maximum number of spots on each component
        parameter (nspomx= 20)

Appendix B The MINGA V2.01 Commands

• ◇ % <comment>

Marks the internal description of the minimization flow.

• > ADD PARAM <num> <nam> <val> <err> <low> <up> <typ> <code>

Adds the LCE parameter to a list of variables (note, that at the last place you must specify the <code> of the LCE parameter).

• ◇ ANALYSE <maxcall> <eps> <VTEST>

Local minimization using SEEK,SIMPLX, and MIGRAD routines. Before execution of the ANALYSE command parameters of the routines mentioned above should be specified. ANALYSE routine may work in two different modes defined by the input file. The name of the input file is specified using command SEED FILE . For description of parameters see command MINIMIZE .

• > ANNEAL <maxcall> <χ^{2}_{stop}> <T_{factor}>

<maxcall> see MIGRAD

<χ^{2}_{stop}> function value to be achieved, <T_{factor}> an initial temperature of system (see [4]).

• ◇ ARRAY FILE <file name>

Defines file name of the file which is used to save the CRS search array.

• ◇ ARRAY FORM <format>

The legal <format> is ASCII (default) or BINARY. In the case of binary format file status is overwrite (see command CRS ). Conversion of the file format may be done with help of the ARR CONV program available on distribution diskette.

• CALL FCN <iflag>

<iflag> is equal to

1 for the first call,

2 calculation of function value and gradient,

3 a) for the last call, see also END ,

b) if on return from FCN iflag is greater than 6, it means, that the values of parameters have been internally redefined by FCN,

4 calculation of function value only,

5,6 reserved.

• > CHANGE CST <code_{1}> <val_{1}> <code_{2}> <val_{2}> <code_{3}> <val_{3}>

Changes the values of non–wavelength dependent LCE parameters (e.g., grid fineness), <code> the code of parameter (see Section 3), <val> the value of parameter.

• CONTOUR Not implemented.

• COVARIANCE <N_{par}>

Reads elements of covariance matrix dimensioned <N_{par}>×<N_{par}> from the *.CMD file. The elements should be given in consecutive lines in format 7F10.0.

• CRS <maxcall> <memory> <nsave> <append>

Controlled Random Search global optimization [14, 2]. This command should be followed by command SHELL which specifies parameters of the clustering analysis. <maxcall> maximum number of calls (default 1000N_{par}), <memory> maximum dimension of the search array (default 30N_{par}), <nsave> save frequency of the search array (default 100N_{par}), <append> available only if the format of the ARRAY FILE defined by command

 ARRAY FORM is ascii,

0 no,

1 yes (default).

• > DF WANT <relχ^{2}min> <relχ^{2}max>

Definition of the limits of <relχ^{2}> (see Section 5), <relχ^{2}min> the minimum relative change of χ^{2} (default 0.0005), <relχ^{2}max> the maximum relative change of χ^{2} (default 0.005).

• END <i>

Marks the end of block of data. After execution of command CALL FCN 3, command END forces MINUIT to reading the next block of data. <i> = 0 reads the next block of data and call FCN with iflag=3, ≠ 0 reads the next block of data.

• END RETURN

This command is similar to the command END , but subroutine MINUIT returns to the main program (MINGA).

• ERROR DEF <UP>

See Section 4.

• EXIT

Marks the end of *.CMD file.

• FIX <num_{1} … num_{7}>

Fixes variable (remove from the minimization)–up to seven variables may be fixed simultaneously.

• GRADIENT <choice> <max iter>

The optional check of quality of the derivatives calculated by MINGA and MINUIT,

<choice>

0 check,

1 no checking (default). <max iter> the number of modifications of the step sizes used by function in calculation of gradient (default 3).

• > GRID <N_{par}> <ngrid> <niter> <siter>

• The first line consist of the numbers of <N_{par}> variables as defined by MINUIT part of the file *.CMD,

• The next <ngrid> lines consist of the values of <N_{par}> parameters in turn (grid points).

This command is especially dedicated to make a grid of solutions (see [15]). This is the best way to locate the global minimum of χ^{2} precisely. The best point found is automatically adopted for further minimization. \underline{Note} : If you specify <N_{par}> equal to 0, then GRID works as any other minimization routine. For this purpose omit all lines, which normally appears after this command. <N_{par}> see above, <ngrid> the total number of grid points, <niter> the maximum number of the MIGRAD minimizations in each grid point (default 5), <siter> as above but for the SIMPLX (default -1, i.e. no SIMPLX minimization.

• > G TYPE

<type>

Definition of the algorithm of derivatives approximation,

<type>

1 asymmetric,

2 symmetric,

3 a three–point rule.

• HELP

Prints out the list of available commands.

• HESSE <maxcall> <trials>

Calculation of Hessian (standard error analysis),

<maxcall> see MIGRAD ,

<trials> if HESSE found new minimum during calculation of Hessian, then it would may try to build Hessian no more than <maxcall> times (default 3).

• > H INITIAL <ij> <trials>

Definition of an initial approximation of Hessian, and the minimization algorithm used by MIGRAD ,

<i> type of algorithm (default 0); for references see [16],
0 DFP/BFGS switching method due to Fletcher,
1 DFP  Davidon–Fletcher–Powell update,
2 BFGS  Broyden–Fletcher–Goldfarb–Shanno update,
3 SR1  Symmetric Rank–one upate,
4 exact  full Hessian (no updates),

<j> approximation of Hessian (default 3),
0 H = I,
1 diagH = χ^{2}_{typ},
2 diagH = \begin{matrix}χ^{2} & _{typ} \\ p⃗ & ^{T} & p⃗ \end{matrix},
3 diagH =approximation by finite differences,
4 exact, full Hessian matrix.

<trials> only if <j>=3, see HESSE .

• IMPROVE <maxcall> <attempts>

Searches for a global minimum, if a local one has been already found,

<maxcall> see MIGRAD ,

<attempts> the number of attempts to found new minimum (default N_{par} + 2).

• ◇ INIT RAND

Initializes random number generators giving random number sequence different from that used during the previous run of MINGA (default is do not initialize).

• LIMITS <n_{1}> <l_{1}> <u_{1}> <n_{2}> <l_{2}> <u_{2}>

Redefinition of the limits of variable(s), <n> the number of variable, <l> the lower limit, <u> the upper limit.

• MATOUT

Prints out or calculates the covariance matrix.

• MIGRAD <maxcall> <eps> <VTEST>

The MIGRAD minimization (see also command H INITIAL ) <maxcall> the maximum number of function calls (default 1000), <eps> default <eps>= 0.1<UP>, <VTEST> see Section 4 (default <VTEST>= 0.01<UP>).

• MINIMIZE <maxcall> <eps> <VTEST>

See MIGRAD .
• MINOS <maxcall> i1i2i3i4. . .

The full error analysis (\underline{very} time consuming),

<maxcall> see MIGRAD ,

i1i2i3i4… the numbers of variables to be analyzed, given in format I2 (default is all variables).

• > M ITER <max number of MIGRAD iterations>

Definition of the maximum number of iterations of the MIGRAD procedure (default 10N_{par}), -1 calculate second derivatives and terminate, -2 no MIGRAD minimization.

• ◇ NORMALIZE

Normalizes the function value to unity as it is needed by MINOS routine (default is do not normalize).

• > ONLY <choice>

This command defines the flow of minimization process after the execution of the command MIGRAD ,

<choice>

0 do not use the SIMPLX routine,

1 perform an additional SIMPLX minimization (default).

• PARAMETER

Reads the new values of variables, given in standard format in consecutive lines, from the file *.CMD.

• ◇ PATH

Shows the trace of the ANALYSE minimization (default is do not show).

• PRINTOUT <density>

Definition of the frequency of output to the file *.SOL,

<density>

0 only in minimum,

1 default,

3 of high density (suggested).

• PUNCH

Print out the values of variables and covariance matrix to the file *.PUN.

• > RAND TYPE <type>

Definition of the type of a random numbers generator,

<type>

0 normal–Gaussian (default), 1 uniform.

• RELEASE <number_{1}>…<number_{7}>

Restores the fixed variables of given numbers,

<number>

the number of fixed variable to be restored.

• REM COREL <choice>

To remove correlations between some parameters. <choice> 0 do not remove correlations (default), 1 remove correlations (see routine ZEROV ).

• > RESTART

See description of the I/O files and the definition of file MIN.RST.

• RESTORE <case>

Restores fixed variable,

<case>

0 all the fixed variables, 1 the last fixed variable.

• > SAVE SOL <freq>

Defines the number of function calls after which the output files will be protected against damage (default 1).

• SCAN Not implemented.

• ◇ SEED FILE <file name>

See descriptions of commands ANALYSE and CRS . <file name> has the following structure

• \underline{line #1}: N_{points}, N_{par}, N_{skip}, where N_{points} is the number of starting points, N_{par} the number of parameters, and N_{skip} the number of points to be skipped,

• \underline{line #2}: should be blank,

• \underline{next N_{par} lines}: an initial estimation of errors of parameters. If N_{points}is less than zero then starting points are automatically generated using random number generator, otherwise, there should be N_{points} blocks each starting with one blank line and the next N_{par} lines of the block defining the initial values of parameters.

• SEEK <maxcall>

The Monte Carlo minimization, <maxcall> see MIGRAD (default 100N_{par}).

• > SET PARAM <n_{1}> <val_{1}> <n_{2}> <val_{2}> <n_{3}> <val_{3}>

Redefinition of the values of variables of given numbers; up to three variables may be changed simultaneously.

• ◇ SHELL <nclusters> <neighb> <edmshell> <rtest> <fraction> <both>

Defines the SHELL cluster analysis parameters. Should be executed immediately after execution of the command CRS . <nclusters> the maximum number of clusters which may be adopted for further analysis (default 2), <neighb> the number of neighbors to be searched around each point (default <memory>/10), <edmshell> the cluster analysis will not be performed until EDM (χ^{2}_{max}/χ^{2}_{min} − 1 inside the search array) is reached (default 2 <UP>). The clustering is obtained if within distance of <rtest> around each cluster center there is <fraction> per cent of <memory> points (defaults: <rtest>=2 metric units defined by UNIT routine, <fraction>=95 % of the search array locations), <both> set 0 to analyze only the centers (default), or 1 to analyze the centers and the best points of each cluster.

• > SHOW TRIAL <freq>

Defines the frequency of output to the screen in units of the function calls (default

• ◇ SIGMA <level>

Error level for the MINOS error analysis in units of standard deviation (default 1).

• SIMPLEX <maxcall> <eps>

See MIGRAD .

• > S ITER <max number of SIMPLX iterations>

Similar to M ITER , but for procedure SIMPLX. If the number is less than 0, then SIMPLX only builds up an initial simplex and terminates (default 10N_{par}).

• STANDARD Not implemented.

• > TOLERANCE <relχ^{2}> <relp⃗> <relG⃗> <EDM> <VTEST>

Definition of the convergence criteria (see Section 4). Default values are

relχ^{2} relp⃗  relG⃗ EDM     VTEST
0.0001 0.0001 0.01 0.0001 0.05

• > TRACE <skip>

Definition of the frequency of output to the file *.TRC (default 0, i.e., no output).

• > TYPICAL F <χ^{2}_{typ}>

Definition of the typical function value (function properly normalized is of order unity–default).

A Short List of Commands

ADD PARAM ANALYSE  ANNEAL  ARRAY FILE
ARRAY FORM FIX  CALL FCN  CHANGE CST
CONTOUR  COVARIANCE CRS  DF WANT
END  END RETURN ERROR DEF EXIT
GRADIENT  GRID  G TYPE  HELP
HESSE  H INITIAL  IMPROVE  INIT RAND
LIMITS  MATOUT  MIGRAD  MINIMIZE
MINOS  M ITER  NORMALIZE ONLY
PARAMETER PATH  PRINTOUT PUNCH
RAND TYPE  RELEASE  REM COREL RESTART
RESTORE  SAVE SOL  SCAN  SEED FILE
SEEK  SET PARAM  SHELL  SHOW TRIAL
SIGMA  SIMPLEX  S ITER  STANDARD
TOLERANCE TRACE  TYPICAL F

References

[1] Ash, J.M., Jones, R.L., 1981, Math. Comp., 37, 159.

[2] Barone, F., Milano, L., Russo, G., 1990, Wilson–Price Code Handbook (Universitá di Napoli).

[3] Barone, F., Maceroni, C., Milano, L., Russo, G., 1988, Astron. Astrophys., 197, 347.

[4] C̆erny, V., 1984, Helsinki preprint HU–TFT–84–51.

[5] Dennis, J.E., Schnabel, R.B., 1983, Numerical Methods for Unconstrained Optimization and Nonlinear Equations (Prentice–Hall Inc.).

[6] Fletcher, R., 1970, Comp. J., 13, 317.

[7] James, F., Roos, M., 1971, MINUIT, CERN, Program Library, long write–up D506(D516) (CERN).

[8] Kratka, M., 1989 (private communication).

[9] Nelder, J.A., Mead, R., 1965, Comp. J., 7, 308.

[10] Plewa, T., 1988, Acta Astron., 38, 415.

[11] Plewa, T., and Kałużny, J., 1992, Acta Astron., 42, 103.

[12] Plewa, T., and Włodarczyk, K.J., 1993, Acta Astron. (in press).

[13] Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986, Numerical Recipes (Cambridge University Press).

[14] Price, W.L., 1976, Comp. J., 20, 367.

[15] Rafert, J.B., Markworth, N.L., 1986, Astron. J., 92, 678.

[16] Shanno, D.F., 1980, Math. Comp., 34, 499.

[17] Szczerba, R., Plewa, T., 1990, in Physics and Composition of Interstellar Matter, International Symposium, Bachotek, Poland, eds. J. Krełowski and J. Papaj (Nicolaus Copernicus University).

[18] Wilson, R.E., 1983, SUBROUTINE DC/LC (with eccentric orbit generalization) (University of Florida).

[19] Wilson, R.E., De Luccia, M.R., Johnston, K., Mango, S.A., 1972, Astrophys. J., 177, 191.

[20] Wilson, R.E., Biermann, P., 1976, Astron. Astrophys., 48, 349.

Tomasz Plewa May 1991