-
Notifications
You must be signed in to change notification settings - Fork 1
/
global_allele_freq.diff
181 lines (169 loc) · 7.17 KB
/
global_allele_freq.diff
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
# this file shows how the global allele frequency calculation was removed from Mendel
# to add it back in need to add back in all the lines that start with (-)
diff --git a/diagnostics.f90 b/diagnostics.f90
index 3d8b16e..b3d8bbe 100644
--- a/diagnostics.f90
+++ b/diagnostics.f90
@@ -1302,7 +1302,7 @@ integer global_mutn_list(MNP,num_tribes)
integer global_list_count(num_tribes)
integer num_falleles(3), num_dalleles(3), lb_limit, dwarn, fwarn
integer par_num_falleles(3), par_num_dalleles(3), mutn_limit
-integer ncount, par_ncount, num_nalleles(3), par_nalleles(3)
+integer ncount, par_ncount, num_nalleles(3), par_num_nalleles(3)
integer nwarn, mutn_thres, mutn, jmax
real*8 dpbin(NB), dpbin_count(NB), dpbin_max, pbin_width
real*8 npbin(NB), npbin_count(NB), npbin_max, x
@@ -1321,11 +1321,7 @@ real bin_fitness(11), bin_center(10)
! mutation is the population size. The polymorphism bin width
! pbin_width value reflects this difference.
-if (is_parallel) then
- pbin_width = 2.*current_pop_size/NB*num_tribes
-else
- pbin_width = 2.*current_pop_size/NB
-end if
+pbin_width = 2.*current_pop_size/NB
if(recombination_model == clonal) pbin_width = pbin_width/2.
@@ -1549,15 +1545,31 @@ endif
! Output global statistics to file caseid.000.plm
if(is_parallel) then
+ call mpi_davg(dpbin, par_dpbin, NB)
+ call mpi_davg(fpbin, par_fpbin, NB)
+ call mpi_davg(npbin, par_npbin, NB)
+
+ call mpi_davg(dpbin_count, par_dpbin_count, NB)
+ call mpi_davg(fpbin_count, par_fpbin_count, NB)
+ call mpi_davg(npbin_count, par_npbin_count, NB)
+
+ call mpi_davg(num_dalleles, par_num_dalleles, 3)
+ call mpi_davg(num_falleles, par_num_falleles, 3)
+ call mpi_davg(num_nalleles, par_num_nalleles, 3)
+
+ dcount = par_num_dalleles(1) + par_num_dalleles(2) + par_dpbin(NB)
+ ncount = par_num_nalleles(1) + par_num_nalleles(2) + par_npbin(NB)
+ fcount = par_num_falleles(1) + par_num_falleles(2) + par_fpbin(NB)
+
if(myid.eq.0 .and. mod(gen,diagnostic_gens)==0 ) then
rewind(21)
write(21,'("# Number of tribes = ",i4)') num_tribes
write(21,'("# generation = ",i8)') gen
write(21,'("# frequency del_normalized fav_normalized", &
- " del_count fav_count neu_normalized neu_count")')
- write(21,'(i11,2f15.11,2f11.0,f15.11,f11.0)') (k, dpbin(k), &
- fpbin(k), npbin(k), dpbin_count(k), fpbin_count(k), &
- npbin_count(k), k=1,NB)
+ " neu_normalized, del_count fav_count neu_count")')
+ write(21,'(i11,2f15.11,2f11.0,f15.11,f11.0)') (k, par_dpbin(k), &
+ par_fpbin(k), par_npbin(k), par_dpbin_count(k), &
+ par_fpbin_count(k), par_npbin_count(k), k=1,NB)
write(21,'("# Allele summary statistics:")')
write(21,'("# Very rare",6x,"Polymorphic",9x,"Fixed", &
11x,"Total")')
@@ -1571,32 +1583,30 @@ if(is_parallel) then
call flush(21)
end if
-else
-
- rewind(11)
- if (mod(gen,plot_allele_gens)==0.and.verbosity>0) then
- write(11,'("# generation = ",i8)') gen
- write(11,'("# frequency del_normalized fav_normalized ", &
- " neu_normalized del_count fav_count neu_count")')
- write(11,'(i11,3f15.11,3f11.0)') (k, dpbin(k), &
- fpbin(k), npbin(k), dpbin_count(k), fpbin_count(k), &
- npbin_count(k), k=1,NB)
- end if
+end if
- ! correct data for total diversity by multiplying by each bin
- rewind(12)
- if (mod(gen,plot_allele_gens)==0.and.verbosity>0) then
- write(12,'("# generation = ",i8)') gen
- write(12,'("# frequency del_normalized fav_normalized ", &
- " neu_normalized del_count fav_count neu_count")')
- write(12,'(i11,3f15.11,3f11.0)') (k, k*dpbin(k), &
- k*fpbin(k), k*npbin(k), k*dpbin_count(k), k*fpbin_count(k), &
- k*npbin_count(k), k=1,NB/2)
- write(12,'(i11,3f15.11,3f11.0)') (k, (NB-k)*dpbin(k), &
- (NB-k)*fpbin(k), (NB-k)*npbin(k), (NB-k)*dpbin_count(k), &
- (NB-k)*fpbin_count(k), (NB-k)*npbin_count(k), k=NB/2,NB)
- end if
+rewind(11)
+if (mod(gen,plot_allele_gens)==0.and.verbosity>0) then
+ write(11,'("# generation = ",i8)') gen
+ write(11,'("# frequency del_normalized fav_normalized ", &
+ " neu_normalized del_count fav_count neu_count")')
+ write(11,'(i11,3f15.11,3f11.0)') (k, dpbin(k), &
+ fpbin(k), npbin(k), dpbin_count(k), fpbin_count(k), &
+ npbin_count(k), k=1,NB)
+end if
+! correct data for total diversity by multiplying by each bin
+rewind(12)
+if (mod(gen,plot_allele_gens)==0.and.verbosity>0) then
+ write(12,'("# generation = ",i8)') gen
+ write(12,'("# frequency del_normalized fav_normalized ", &
+ " neu_normalized del_count fav_count neu_count")')
+ write(12,'(i11,3f15.11,3f11.0)') (k, k*dpbin(k), &
+ k*fpbin(k), k*npbin(k), k*dpbin_count(k), k*fpbin_count(k), &
+ k*npbin_count(k), k=1,NB/2)
+ write(12,'(i11,3f15.11,3f11.0)') (k, (NB-k)*dpbin(k), &
+ (NB-k)*fpbin(k), (NB-k)*npbin(k), (NB-k)*dpbin_count(k), &
+ (NB-k)*fpbin_count(k), (NB-k)*npbin_count(k), k=NB/2,NB)
end if
! Keep a "snapshot" file of latest polymorphism data.
@@ -1611,7 +1621,7 @@ if (verbosity > 0) then
call flush(13)
end if
-if(.not.is_parallel .and. mod(gen,plot_allele_gens)==0.and.verbosity>0) then
+if(mod(gen,plot_allele_gens)==0.and.verbosity>0) then
write(11,"('#',11x,'Allele summary statistics (tracked mutations only):'/ &
'# (Statistics are based on ',i12,' tracked deleterious mutations'/ &
'# ',i12,' tracked favorable mutations'/ &
@@ -1873,51 +1883,6 @@ do i=1,current_pop_size
end do
end do
-! START_MPI
-if(is_parallel) then
-
- ! For parallel cases, gather all the data together
- ! and do a single global analysis on processor 0.
-
- call mpi_gather(mutn_list,MNP,mpi_integer, &
- global_mutn_list,MNP,mpi_integer,0,mycomm,ierr)
- call mpi_gather(mutn_count,MNP,mpi_integer, &
- global_mutn_count,MNP,mpi_integer,0,mycomm,ierr)
- call mpi_gather(list_count,1,mpi_integer, &
- global_list_count,1,mpi_integer,0,mycomm,ierr)
- call mpi_barrier(mycomm,ierr)
-
- if(myid.eq.0) then
-
- ! Rebuild mutn_list and mutn_count arrays from data
- ! aggregated from all the tribes.
-
- mutn_list = 0
- mutn_count = 0
- list_count = 0
-
- do m = 1,num_tribes
- do i = 1, global_list_count(m)
- new_mutn = .true.
- do k = 1, list_count
- if(global_mutn_list(i,m).eq.mutn_list(k)) then
- mutn_count(k) = mutn_count(k) + global_mutn_count(i,m)
- new_mutn = .false.
- end if
- end do
- if(new_mutn) then
- list_count = min(MNP, list_count + 1)
- mutn_list(list_count) = global_mutn_list(i,m)
- mutn_count(list_count) = global_mutn_count(i,m)
- end if
- end do
- end do
-
- end if
-
-end if
-! END_MPI
-
end subroutine count_alleles
subroutine count_string_alleles(xmutn,max_mutn_per_indiv,MNP, &