1 /* \file joint-chains.c
2
3 Calculates Geyer's reverse logistic regression
4
5 like(param) = Sum_g^G( p(g|param)/ denom)
6
7 denom = Sum_j^chains (n_j P(g|param0_j) / L(param_j)
8
9 \- pick param
10 \- calc chainparamlike
11 \- solve paramlike iteratively
12 \- maximize paramlike
13
14
15 */
16
17 /* Joint chain estimator
18 Peter Beerli January 2000
19
20 Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
21 Copyright 2003-2004 Peter Beerli, Tallahassee FL
22
23 Permission is hereby granted, free of charge, to any person obtaining a copy
24 of this software and associated documentation files (the "Software"), to deal
25 in the Software without restriction, including without limitation the rights
26 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
27 of the Software, and to permit persons to whom the Software is furnished to do
28 so, subject to the following conditions:
29
30 The above copyright notice and this permission notice shall be included in all copies
31 or substantial portions of the Software.
32
33 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
34 INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
35 PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
36 HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
37 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
38 OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
39
40
41 $Id: joint-chains.c 2067 2012-07-27 20:59:32Z beerli $
42 */
43
44 #include "migration.h"
45 #include "broyden.h"
46 #include "combroyden.h"
47 #include "sighandler.h"
48
49 #define EPSILON5 0.00001
50
51 ///
52 /// Precalculate the prob(G|param0) for a replicate and a locus
53 void
create_multiapg0(MYREAL * apg0,nr_fmt * nr,long rep,long locus)54 create_multiapg0 (MYREAL *apg0, nr_fmt * nr, long rep, long locus)
55 {
56 long g, r;
57 MYREAL *tmp;
58 MYREAL tmpmax; //overflow protector
59 MYREAL reps = LOG ((MYREAL) (nr->repstop - nr->repstart));
60
61 tmp = (MYREAL *) mycalloc (nr->repstop, sizeof (MYREAL));
62 // for each tree summary calculate the Sum over replicates
63 for (g = 0; g < nr->world->atl[rep][locus].T; g++)
64 {
65 tmpmax = -MYREAL_MAX;
66 for (r = nr->repstart; r < nr->repstop; r++)
67 {
68 tmp[r] = -reps +
69 probG (nr->world->atl[r][locus].param0,
70 nr->world->atl[r][locus].lparam0,
71 &nr->world->atl[rep][locus].tl[g], nr, locus)
72 - (nr->world->chainlikes[locus][r]);
73
74 if (tmp[r] > tmpmax)
75 tmpmax = tmp[r];
76 }
77 apg0[g] = 0.0;
78 for (r = nr->repstart; r < nr->repstop; r++)
79 apg0[g] += EXP (tmp[r] - tmpmax);
80 apg0[g] = LOG (apg0[g]) + tmpmax;
81 }
82 myfree(tmp);
83 }
84
85 ///
86 /// Solve the Geyer (1994) equation system
87 /// \f[
88 /// L(\Theta_0,M_0) = \sum_{chains} \frac{P(G|\Theta_0,M_0)}{\sum_{chains}\frac{P(G|\Theta_0,M_0)}{L(\Theta_0,M_0)}}
89 /// \f]
90 void
interpolate_like(nr_fmt * nr,long locus)91 interpolate_like (nr_fmt * nr, long locus)
92 {
93 boolean alldone;
94 boolean diffdone=FALSE;
95 long j, z, r;
96 long repdiff = nr->repstop - nr->repstart;
97 MYREAL *newlike, *lparam;
98 MYREAL *diff, *oldiff;
99 MYREAL delta = 0.;
100 MYREAL *oldlike;
101 MYREAL *chainlike = nr->world->chainlikes[locus];
102 //#ifdef INTEGRATEDLIKE
103 //return;
104 //#endif
105 oldlike = (MYREAL *) mycalloc (repdiff, sizeof (MYREAL));
106 memcpy (oldlike, chainlike, sizeof (MYREAL) * repdiff);
107
108 newlike = (MYREAL *) mycalloc (repdiff, sizeof (MYREAL));
109 diff = (MYREAL *) mycalloc (repdiff, sizeof (MYREAL));
110 oldiff = (MYREAL *) mycalloc (repdiff, sizeof (MYREAL));
111 lparam = (MYREAL *) mycalloc (nr->numpop2, sizeof (MYREAL));
112 /* following material is reverse logistic regression
113 a la Geyer 1994 */
114 z = 0;
115 alldone = FALSE;
116 // minimize the difference between oldlike and newlike, using oldiff
117 memset (oldiff, 0, sizeof (MYREAL) * repdiff);
118 while (!alldone && z++ < 10000)
119 {
120 alldone = TRUE;
121
122 for (r = nr->repstart; r < nr->repstop; r++)
123 create_multiapg0 (nr->apg0[r][locus], nr, r, locus);
124
125 for (j = nr->repstart; j < nr->repstop; j++)
126 {
127 newlike[j] = calc_locus_like (nr, nr->atl[j][locus].param0,
128 nr->atl[j][locus].lparam0, locus);
129 diff[j] = fabs (newlike[j] - oldlike[j]);
130 if (delta < diff[j])
131 delta = diff[j];
132 if (delta > EPSILON)
133 alldone = FALSE;
134 }
135
136 if (nr->world->options->progress)
137 {
138 if (z % 100 == 0 && nr->world->options->verbose)
139 printf (" Iteration%6li biggest difference = %f\n", z, delta);
140 }
141
142 memcpy (oldlike, newlike, sizeof (MYREAL) * repdiff);
143 diffdone = TRUE;
144 for (j = nr->repstart; j < nr->repstop; j++)
145 {
146 if (fabs (oldiff[j] - diff[j]) > EPSILON)
147 {
148 diffdone = FALSE;
149 break;
150 }
151 }
152 if (diffdone)
153 break;
154 else
155 memcpy (oldiff, diff, sizeof (MYREAL) * repdiff);
156 // reset delta to the first log(L) difference, this will enter again in the loop
157 // above
158 delta = diff[0];
159 }
160 if ((nr->world->options->progress || diffdone) && delta > EPSILON)
161 {
162 printf (" Reweighting operation converged to\n");
163 printf (" constant multiplier %f in %li cycles\n", delta, z);
164 }
165 myfree(newlike);
166 myfree(oldlike);
167 myfree(diff);
168 myfree(oldiff);
169 myfree(lparam);
170 }
171