source: tags/initial/GDE/MOLPHY/optimtpm.c

Last change on this file was 2, checked in by oldcode, 25 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.1 KB
Line 
1/*
2 * optimtpm.c   Adachi, J.   1994.11.16
3 * Copyright (C) 1992-1994 J. Adachi & M. Hasegawa; All rights reserved.
4 */
5
6double
7abfunc(ab)
8double ab;
9{
10        Node *rp;
11
12        AlphaBeta = ab;
13        tranprobmat();
14        fmlength(Ctree, Distanmat, Maxspc);
15/*      lslength(Ctree, Distanvec, Maxspc); */
16        initpartlkl(Ctree);
17        rp = (Node *)mlikelihood(Ctree);
18        return Ctree->lklhd;
19} /* abfunc */
20
21
22double
23ayrfunc(ayr)
24double ayr;
25{
26        Node *rp;
27
28        AlphaYR = ayr;
29        tranprobmat();
30        fmlength(Ctree, Distanmat, Maxspc);
31/*      lslength(Ctree, Distanvec, Maxspc); */
32        initpartlkl(Ctree);
33        rp = (Node *)mlikelihood(Ctree);
34        return Ctree->lklhd;
35} /* ayrfunc */
36
37
38double
39b12func(b12)
40double b12;
41{
42        Node *rp;
43
44        Beta12 = b12;
45        tranprobmat();
46        fmlength(Ctree, Distanmat, Maxspc);
47/*      lslength(Ctree, Distanvec, Maxspc); */
48        initpartlkl(Ctree);
49        rp = (Node *)mlikelihood(Ctree);
50        return Ctree->lklhd;
51} /* b12func */
52
53
54double
55optimab(x1, x2, x3, f1, f2, f3, func, eps)
56double x1, x2, x3, f1, f2, f3;
57double (*func)();
58double eps;
59{
60        int i;
61        double xa, fa, dl, dr, gold;
62
63        gold = 0.38197;
64        for (i = 0; (x3 - x1) > eps; i++) {
65                if (Info_optn)
66                        printf("%-3d%8.3f%8.3f%8.3f%13.3f%13.3f%13.3f\n",
67                                i,x1,x2,x3,f1,f2,f3);
68                dl = x2 - x1;
69                dr = x3 - x2;
70                if (dl > dr) {
71                        xa = x2 - gold * dl;
72                } else {
73                        xa = x2 + gold * dr;
74                }
75                fa = (*func)(xa);
76                if (dl > dr) {
77                        if (fa > f2) {
78                                x3 = x2; f3 = f2; x2 = xa; f2 = fa;
79                        } else {
80                                x1 = xa; f1 = fa;
81                        }
82                } else {
83                        if (fa > f2) {
84                                x1 = x2; f1 = f2; x2 = xa; f2 = fa;
85                        } else {
86                                x3 = xa; f3 = fa;
87                        }
88                }
89                if (i > 30) break;
90        }
91        if (f3 > f2) {
92                x2 = x3; f2 = f3;
93        } else if (f1 > f2) {
94                x2 = x1; f2 = f1;
95        }
96/*      printf("  opt. para: %.3f  Max. ln likelihood: %.3f\n", x2, f2); */
97        return x2;
98} /* optimab */
99
100
101void
102optimtpm()
103{
104        int i, loop;
105        char buf[64];
106        double x1, x2, x3, f1, f2, f3, ab, ayr, b12, eps1, eps2, eps;
107        double abold, ayrold, b12old, dab, dayr, db12;
108        boolean ayr_mode, b12_mode;
109
110#if 0
111        for (i = 5; i < 500; i += 5) {
112                ab = 0.1 * i;
113                f1 = abfunc(ab);
114                printf("%5.1f %15.3f\n", ab, f1);
115        }
116        exit(1);
117#endif
118
119        Topting = TRUE;
120        eps = 0.005;
121/*      printf("Alpha/Beta: %.3f  AlphaY/AlphaR: %.3f\n", AlphaBeta, AlphaYR); */
122        ayr_mode = FALSE;
123        b12_mode = FALSE;
124        if (Toptim_optn == 1) {
125                eps1 = eps;
126        } else {
127                if (Toptim_optn == 3) {
128                        ayr_mode = TRUE;
129                } else if (Toptim_optn == 5) {
130                        b12_mode = TRUE;
131                } else {
132                        ayr_mode = TRUE;
133                        b12_mode = TRUE;
134                }
135                eps1 = 4.0;
136        }
137        x1 = 10.0;
138        x3 = 15.0;
139        f1 = abfunc(x1);
140        f3 = abfunc(x3);
141        if (f1 > f3) {
142                x2 = x1; f2 = f1;
143                x1 = MINAB;
144                f1 = abfunc(x1);
145        } else {
146                x2 = x3; f2 = f3;
147                x3 = MAXAB;
148                f3 = abfunc(x3);
149        }
150        AlphaBeta = ab = optimab(x1, x2, x3, f1, f2, f3, abfunc, eps1);
151
152        if (ayr_mode) {
153                if (Info_optn) putchar('\n');
154                eps2 = 0.1;
155                x1 = 0.5;
156                x3 = 1.0;
157                f1 = ayrfunc(x1);
158                f3 = ayrfunc(x3);
159                if (f1 > f3) {
160                        x2 = x1; f2 = f1;
161                        x1 = MINAYR;
162                        f1 = ayrfunc(x1);
163                } else {
164                        x2 = x3; f2 = f3;
165                        x3 = MAXAYR;
166                        f3 = ayrfunc(x3);
167                }
168                AlphaYR = ayr = optimab(x1, x2, x3, f1, f2, f3, ayrfunc, eps2);
169
170                dab = 9.0;
171                dayr = 0.2;
172                if (Info_optn) putchar('\n');
173                if (Info_optn) printf("A/B: %.3f  Ay/Ar: %.3f", ab,ayr);
174                if (Info_optn) printf("  %.4f %.4f", dab,dayr);
175                if (Info_optn) putchar('\n');
176                for (loop = 0; dab > eps || dayr > eps; loop++) {
177                        if (loop < 1) {
178                                eps1 = 2.0;  eps2 = 0.05;
179                        } else if (loop < 2) {
180                                eps1 = 0.5;  eps2 = eps;
181                        } else if (loop < 3) {
182                                eps1 = 0.05; eps2 = eps;
183                        } else {
184                                eps1 = eps; eps2 = eps;
185                        }
186       
187                        if (dab < eps) dab = eps;
188                        x1 = ab;
189                        (ab + 1.0) < MAXAB ? (x3 = ab + 1.0) : (x3 = MAXAB);
190                        f1 = abfunc(x1);
191                        f3 = abfunc(x3);
192                        if (f1 > f3) {
193                                x2 = x1; f2 = f1;
194                                (x2 - dab) > MINAB ? (x1 = x2 - dab) : (x1 = MINAB);
195                                f1 = abfunc(x1);
196                        } else {
197                                x2 = x3; f2 = f3;
198                                (x2 + dab) < MAXAB ? (x3 = x2 + dab) : (x3 = MAXAB);
199                                f3 = abfunc(x3);
200                        }
201                        abold = ab;
202                        AlphaBeta = ab = optimab(x1, x2, x3, f1, f2, f3, abfunc, eps1);
203                        dab = fabs(ab - abold);
204       
205                        if (dayr < eps) dayr = eps;
206                        (ayr - 0.02) > MINAYR ? (x1 = ayr - 0.02) : (x1 = MINAYR);
207                        x3 = ayr;
208                        f1 = ayrfunc(x1);
209                        f3 = ayrfunc(x3);
210                        if (f1 > f3) {
211                                x2 = x1; f2 = f1;
212                                (x2 - dayr) > MINAYR ? (x1 = x2 - dayr) : (x1 = MINAYR);
213                                f1 = ayrfunc(x1);
214                        } else {
215                                x2 = x3; f2 = f3;
216                                (x2 + dayr) < MAXAYR ? (x3 = x2 + dayr) : (x3 = MAXAYR);
217                                f3 = ayrfunc(x3);
218                        }
219                        ayrold = ayr;
220                        AlphaYR = ayr = optimab(x1, x2, x3, f1, f2, f3, ayrfunc, eps2);
221                        dayr = fabs(ayr - ayrold);
222                        if (Info_optn) putchar('\n');
223                        if (Info_optn) printf("A/B: %.3f  Ay/Ar: %.3f", ab,ayr);
224                        if (Info_optn) printf("  %.4f %.4f", dab,dayr);
225                        if (Info_optn) putchar('\n');
226                }
227        }
228
229        if (b12_mode) {
230                if (Info_optn) putchar('\n');
231                eps2 = 0.1;
232                x1 = 0.5;
233                x3 = 1.0;
234                f1 = b12func(x1);
235                f3 = b12func(x3);
236                if (f1 > f3) {
237                        x2 = x1; f2 = f1;
238                        x1 = MINB12;
239                        f1 = b12func(x1);
240                } else {
241                        x2 = x3; f2 = f3;
242                        x3 = MAXB12;
243                        f3 = b12func(x3);
244                }
245                Beta12 = b12 = optimab(x1, x2, x3, f1, f2, f3, b12func, eps2);
246
247                dab = 9.0;
248                db12 = 0.2;
249                if (Info_optn) putchar('\n');
250                if (Info_optn) printf("A/B: %.3f  B1/B2: %.3f", ab,b12);
251                if (Info_optn) printf("  %.4f %.4f", dab,db12);
252                if (Info_optn) putchar('\n');
253                for (loop = 0; dab > eps || db12 > eps; loop++) {
254                        if (loop < 1) {
255                                eps1 = 2.0;  eps2 = 0.05;
256                        } else if (loop < 2) {
257                                eps1 = 0.5;  eps2 = eps;
258                        } else if (loop < 3) {
259                                eps1 = 0.05; eps2 = eps;
260                        } else {
261                                eps1 = eps; eps2 = eps;
262                        }
263       
264                        if (dab < eps) dab = eps;
265                        x1 = ab;
266                        (ab + 1.0) < MAXAB ? (x3 = ab + 1.0) : (x3 = MAXAB);
267                        f1 = abfunc(x1);
268                        f3 = abfunc(x3);
269                        if (f1 > f3) {
270                                x2 = x1; f2 = f1;
271                                (x2 - dab) > MINAB ? (x1 = x2 - dab) : (x1 = MINAB);
272                                f1 = abfunc(x1);
273                        } else {
274                                x2 = x3; f2 = f3;
275                                (x2 + dab) < MAXAB ? (x3 = x2 + dab) : (x3 = MAXAB);
276                                f3 = abfunc(x3);
277                        }
278                        abold = ab;
279                        AlphaBeta = ab = optimab(x1, x2, x3, f1, f2, f3, abfunc, eps1);
280                        dab = fabs(ab - abold);
281       
282                        if (db12 < eps) db12 = eps;
283                        (b12 - 0.02) > MINB12 ? (x1 = b12 - 0.02) : (x1 = MINB12);
284                        x3 = b12;
285                        f1 = b12func(x1);
286                        f3 = b12func(x3);
287                        if (f1 > f3) {
288                                x2 = x1; f2 = f1;
289                                (x2 - db12) > MINB12 ? (x1 = x2 - db12) : (x1 = MINB12);
290                                f1 = b12func(x1);
291                        } else {
292                                x2 = x3; f2 = f3;
293                                (x2 + db12) < MAXB12 ? (x3 = x2 + db12) : (x3 = MAXB12);
294                                f3 = b12func(x3);
295                        }
296                        b12old = b12;
297                        Beta12 = b12 = optimab(x1, x2, x3, f1, f2, f3, b12func, eps2);
298                        db12 = fabs(b12 - b12old);
299                        if (Info_optn) putchar('\n');
300                        if (Info_optn) printf("A/B: %.3f  B1/B2: %.3f", ab,b12);
301                        if (Info_optn) printf("  %.4f %.4f", dab,db12);
302                        if (Info_optn) putchar('\n');
303                }
304        }
305
306        if (!Ctacit_optn) {
307                printf("Alpha/Beta: %.3f", AlphaBeta);
308                if (AlphaYR != 1.0) printf("  AlphaY/AlphaR: %.3f", AlphaYR);
309                if (Beta12  != 1.0) printf("  Beta1/Beta2: %.3f", Beta12);
310                putchar('\n');
311        }
312
313#if 0
314        *Modelname = '\0';
315        sprintf(buf, "A/B:%.2f", AlphaBeta);
316        strcat(Modelname, buf);
317        if (AlphaYR != 1.0) {
318                sprintf(buf, " Ay/Ar:%.2f", AlphaYR);
319                strcat(Modelname, buf);
320        }
321        if (Beta12 != 1.0) {
322                sprintf(buf, " B1/B2:%.2f", Beta12);
323                strcat(Modelname, buf);
324        }
325        if (Frequ_optn) strcat(Modelname, " F");
326/*      printf("\"%s\"\n", Modelname); */
327#endif
328
329        Topting = FALSE;
330} /* optimtpm */
Note: See TracBrowser for help on using the repository browser.