1 /******************************************************************************
2
3 #### ##### ##### ###### ##### #### ####
4 # # # # # # # # # # # #
5 # # # # # # ##### # # #### #
6 # # ##### # # # ##### # ### #
7 # # # # # # # # # # # ### # #
8 #### # # ##### ###### # # #### ### ####
9
10 ******************************************************************************/
11 /* This file is part of MAPMAKER 3.0b, Copyright 1987-1992, Whitehead Institute
12 for Biomedical Research. All rights reserved. See READ.ME for license. */
13
14 #define INC_LIB
15 #define INC_SHELL
16 #define INC_MISC
17 #include "mapm.h"
18
19 #define set_excluded(a,b,c) \
20 three_pt_excluded[three_pt_index[a]][three_pt_index[b]][three_pt_index[c]] \
21 =EXCLUDED /* OBSOLETE? */
22
23 void left_exclusion(), right_exclusion(), middle_exclusion();
24 bool is_excluded(), is_computed();
25
26 void add_to_order();
27 int compare_markers_to_try();
28 void randomize_markers_to_try();
29 void rank_markers_to_try();
30 void add_order();
31 bool extend_order_by_1();
32 int **add_orders();
33 bool middles_excluded();
34
35 char ***three_pt_excluded=NULL;
36 int three_pt_size=0;
37 int *three_pt_index=NULL;
38
39 #define EXCLUDED 0
40 #define ALLOWED 1
41 #define UNCOMPUTED 2
42
43 void alloc_3pt_matrix(); /* args: int num_loci; sets a global, will free
44 and realloc if it's too small, use to prealloc for efficiency */
45 void free_3pt_matrix();
46
47
48
49 /**************** Two-point, Haplotypes ****************/
50
get_linkage_group(locus,num_loci,linkage_group,group_size,lodbound,thetabound)51 void get_linkage_group(locus,num_loci,linkage_group,group_size,lodbound,
52 thetabound)
53 int *locus, *num_loci, *linkage_group, *group_size;
54 real lodbound,thetabound;
55 {
56 int i, j, linked, *unlinked_list, unlinked_size, oldsize, a, b;
57 real lod, theta, thetam, thetaf;
58
59 /* start linkage_group with first marker in the list of loci */
60 linkage_group[0]=locus[0];
61 *group_size= 1;
62 unlinked_list= locus;
63 *num_loci-= 1;
64 unlinked_size= *num_loci;
65 locus= &unlinked_list[1];
66
67 do {
68 oldsize= unlinked_size;
69 unlinked_size= 0;
70 for (i=0; i<oldsize; i++) { /* for all loci a NOT in the group */
71 a= locus[i];
72 for (j=0, linked=FALSE; j<*group_size; j++) { /* for all in group*/
73 b= linkage_group[j];
74 if (!sex_specific) {
75 get_two_pt(a,b,&lod,&theta);
76 if (lod>=lodbound && theta<=thetabound)
77 { linked=TRUE; break; }
78 } else {
79 get_sex_two_pt(a,b,&lod,&thetam,&thetaf);
80 if (lod>lodbound &&
81 (thetam<=thetabound || thetaf<=thetabound))
82 { linked=TRUE; break; }
83 }
84 }
85 if (linked) {
86 linkage_group[*group_size]= a;
87 ++*group_size;
88 } else {
89 /* Effectively, write back into the locus array, relying on the
90 fact that unlinked_size<=*num_loci */
91 unlinked_list[unlinked_size]= a;
92 ++unlinked_size;
93 }
94 } /* for i */
95 locus= unlinked_list;
96 *num_loci= unlinked_size;
97 } while (unlinked_size!=oldsize && unlinked_size!=0);
98 }
99
100
find_haplo_group(locus,num_loci,haplo_group,num_haplo,old_obs,new_obs)101 void find_haplo_group(locus,num_loci,haplo_group,num_haplo,old_obs,new_obs)
102 int *locus, *num_loci, *haplo_group, *num_haplo;
103 int *old_obs, *new_obs;
104 {
105 int i, num_inf, num_dom, best_i, foo;
106 int best_score, old_score, score;
107
108 /* get the most informative marker */
109 best_score= 0;
110 for (i=0; i<*num_loci; i++) {
111 f2_genotype(locus[i],FALSE,old_obs);
112 num_inf= f2_count_infs(&num_dom,&foo,old_obs);
113 score= 2*num_inf - num_dom;
114 if (score>best_score) { best_score=score; best_i=i; }
115 }
116
117 /* delete it, and make it the haplo group */
118 haplo_group[0]=locus[best_i]; *num_haplo=1;
119 f2_genotype(locus[best_i],FALSE,old_obs);
120 locus[best_i]= locus[--*num_loci];
121 if (*num_loci==0) return;
122
123 do { /* find one to add to the group */
124 old_score= best_score;
125 best_score= 0;
126 for (i=0; i<*num_loci; i++)
127 if (merge_genotypes(locus[i],old_obs,new_obs)) { /* if no recs */
128 num_inf= f2_count_infs(&num_dom,&foo,new_obs);
129 score= 2*num_inf - num_dom;
130 if (score>best_score) { best_score=score; best_i=i; }
131 }
132 if (best_score !=0 ) { /* if we found any with no recs */
133 merge_genotypes(locus[best_i],old_obs,old_obs); /*updates old_obs*/
134 haplo_group[*num_haplo]=locus[best_i]; ++*num_haplo; /* add it */
135 locus[best_i]= locus[--*num_loci]; /* delete it */
136 }
137 } while (*num_loci>0 && best_score!=0); /* got one and more to try */
138 }
139
140
141
142 /**************** Three Point ****************/
143
144
alloc_3pt_matrix(num_group)145 void alloc_3pt_matrix(num_group)
146 int num_group;
147 {
148 if (three_pt_excluded!=NULL && num_group<=three_pt_size) return;
149
150 run {
151 if (three_pt_excluded!=NULL) free_3pt_matrix();
152 three_pt_size=num_group;
153 three_pt_excluded= alloc_char_3d_matrix(num_group,num_group,num_group);
154 array(three_pt_index,raw.num_markers,int);
155 } when_aborting {
156 free_3pt_matrix();
157 relay_messages;
158 }
159 }
160
161
162
free_3pt_matrix()163 void free_3pt_matrix()
164 {
165 free_char_3d_matrix(three_pt_excluded,three_pt_size,three_pt_size);
166 unarray(three_pt_index,int);
167 three_pt_excluded=NULL; three_pt_size=0;
168 }
169
170
171
setup_3pt_data(locus,num_loci,threshold)172 void setup_3pt_data(locus,num_loci,threshold)
173 int *locus, num_loci;
174 real threshold;
175 {
176 int i, j, k;
177 real d1, d2, d3;
178
179 alloc_3pt_matrix(num_loci);
180
181 /* three_pt_index[] is lookup by raw locus# to three_pt_exclusion index */
182 for (i=0; i<raw.num_markers; i++) three_pt_index[i]= -1;
183 for (i=0; i<num_loci; i++) three_pt_index[locus[i]]=i;
184
185 for (i=0; i<num_loci; i++)
186 for (j=0; j<num_loci; j++)
187 for (k=0; k<num_loci; k++)
188 three_pt_excluded[i][j][k]= UNCOMPUTED;
189
190 for (i=0; i<num_loci-2; i++) {
191 for (j=i+1; j<num_loci-1; j++) {
192 for (k=j+1; k<num_loci; k++) {
193
194 if (!restore_triple(locus[i],locus[j],locus[k],&d1,&d2,&d3))
195 continue;
196
197 /* maybe later compute on the fly ADD keep_user_amused
198 compute_triple(locus[i],locus[j],locus[k],&d1,&d2,&d3); */
199
200 if (d1<threshold) {
201 three_pt_excluded[i][j][k]=EXCLUDED;
202 three_pt_excluded[k][j][i]=EXCLUDED;
203 } else {
204 three_pt_excluded[i][j][k]=ALLOWED;
205 three_pt_excluded[k][j][i]=ALLOWED;
206 }
207
208 if (d2<threshold) {
209 three_pt_excluded[i][k][j]=EXCLUDED;
210 three_pt_excluded[j][k][i]=EXCLUDED;
211 } else {
212 three_pt_excluded[i][k][j]=ALLOWED;
213 three_pt_excluded[j][k][i]=ALLOWED;
214 }
215
216 if (d3<threshold) {
217 three_pt_excluded[j][i][k]=EXCLUDED;
218 three_pt_excluded[k][i][j]=EXCLUDED;
219 } else {
220 three_pt_excluded[j][i][k]=ALLOWED;
221 three_pt_excluded[k][i][j]=ALLOWED;
222 }
223 }
224 }
225 }
226 }
227
free_3pt_data()228 void free_3pt_data() { free_3pt_matrix(); }
229
230
start_3pt(locus,num_loci,threshold,order,num_ordered)231 bool start_3pt(locus,num_loci,threshold,order,num_ordered) /* improve this! */
232 int *locus, *num_loci; /* starters are deleted */
233 real threshold;
234 int *order, *num_ordered; /* the starting order */
235 {
236 int i, j, k, num_left, starter[3], temp;
237 real d1, d2, d3, best;
238
239 starter[0]= starter[1]= starter[2]= NO_LOCUS;
240 best= 0.0;
241
242 for (i=0; i<*num_loci-2; i++) {
243 for (j=i+1; j<*num_loci-1; j++) {
244 for (k=j+1; k<*num_loci; k++) {
245 if (!restore_triple(locus[i],locus[j],locus[k],&d1,&d2,&d3))
246 continue;
247
248 if (d2<best && d3<best) {
249 best=min(d2,d3); if (best>threshold) { best=0.0; continue;}
250 starter[0]=i; starter[1]=j; starter[2]=k;
251 } else if (d1<best && d3<best) {
252 best=min(d1,d3); if (best>threshold) { best=0.0; continue;}
253 starter[0]=i; starter[1]=k; starter[2]=j;
254 } else if (d1<best && d2<best) {
255 best=min(d1,d2); if (best>threshold) { best=0.0; continue;}
256 starter[0]=j; starter[1]=i; starter[2]=k;
257 }
258 }
259 }
260 }
261
262 if (best==0.0) return(FALSE);
263
264 /* else, setup order list and delete starters from locus list */
265 order[0]=locus[starter[0]]; locus[starter[0]]=NO_LOCUS;
266 order[1]=locus[starter[1]]; locus[starter[1]]=NO_LOCUS;
267 order[2]=locus[starter[2]]; locus[starter[2]]=NO_LOCUS;
268 *num_ordered= 3;
269
270 num_left=0;
271 for (i=0; i<*num_loci; i++) {
272 temp=locus[i];
273 if (temp!=NO_LOCUS) locus[num_left++]=temp;
274 }
275 *num_loci= num_left;
276 return(TRUE);
277 }
278
279
280
is_excluded(a,b,c)281 bool is_excluded(a,b,c)
282 int a, b, c;
283 { return(three_pt_excluded[three_pt_index[a]][three_pt_index[b]]
284 [three_pt_index[c]]==EXCLUDED); }
285
286
is_computed(a,b,c)287 bool is_computed(a,b,c)
288 int a, b, c;
289 { return(three_pt_excluded[three_pt_index[a]][three_pt_index[b]]
290 [three_pt_index[c]]!=UNCOMPUTED); }
291
292
left_exclusion(excluded,locus,num_loci,newmarker)293 void left_exclusion(excluded,locus,num_loci,newmarker)
294 bool *excluded;
295 int *locus, num_loci, newmarker;
296 {
297 int middle, end, i;
298
299 for (middle=num_loci-2; middle>=0; middle--) {
300 for (end=num_loci-1; end>middle; end--) {
301 if (is_excluded(newmarker,locus[middle],locus[end])) {
302 for (i=0; i<=middle; i++)
303 excluded[i]=TRUE;
304 return;
305 }
306 }
307 }
308 }
309
310
right_exclusion(excluded,locus,num_loci,newmarker)311 void right_exclusion(excluded,locus,num_loci,newmarker)
312 bool *excluded;
313 int *locus, num_loci, newmarker;
314 {
315 int middle,end,i;
316
317 for (middle=1; middle<=num_loci-1; middle++) {
318 for (end=0; end<middle; end++) {
319 if (is_excluded(locus[end],locus[middle],newmarker)) {
320 for (i=num_loci; i>middle; i--)
321 excluded[i]=TRUE;
322 return;
323 }
324 }
325 }
326 }
327
328
middle_exclusion(excluded,locus,num_loci,leftmost,rightmost,init_window_size,newmarker)329 void middle_exclusion(excluded,locus,num_loci,leftmost,rightmost,
330 init_window_size,newmarker)
331 bool *excluded;
332 int *locus, num_loci;
333 int leftmost, rightmost; /* presumably, the left and rightmost intervals? */
334 int init_window_size, newmarker;
335 {
336 int first, last, start, end, left, width, i, j;
337
338 for (width=init_window_size; width>=1; width--) {
339 first=max(1,leftmost-width+1);
340 last= min(rightmost,num_loci-width);
341 for (i=first; i<=last; i++) {
342 start=i;
343 end=i+width-1;
344 if (is_excluded(locus[start-1],newmarker,locus[end])) {
345 for (j=start; j<=end; j++)
346 excluded[j]=TRUE;
347 if (start>leftmost && width>1)
348 middle_exclusion(excluded,locus,num_loci,leftmost,start-1,
349 width-1,newmarker);
350 if (end==rightmost) return;
351 leftmost=end+1;
352 }
353 }
354 }
355 }
356
357
358 /* orig frame: a - b - c - d - e - f - g
359 orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6 num_placed=7
360 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 */
361
three_pt_exclusions(order,num_placed,newmarker,excluded)362 int three_pt_exclusions(order,num_placed,newmarker,excluded)
363 int *order, num_placed, newmarker;
364 bool *excluded; /* side-effected */
365 {
366 int i, j, left, right, num_places;
367 for (i=0; i<num_placed+1; i++) excluded[i]=FALSE;
368
369 left_exclusion(excluded,order,num_placed,newmarker);
370 right_exclusion(excluded,order,num_placed,newmarker);
371
372 left=1; right=num_placed;
373 j=1; while (j<num_placed && excluded[j]) j++;
374 if (j!=num_placed)
375 middle_exclusion(excluded,order,num_placed,left,right,right-left,
376 newmarker);
377 for (num_places=0, i=0; i<num_placed+1; i++)
378 if (!excluded[i]) ++num_places;
379 return(num_places);
380 }
381
382
three_pt_verify(locus,num_loci,window_size)383 bool three_pt_verify(locus,num_loci,window_size)
384 int *locus, num_loci, window_size;
385 {
386 int i, j, k, left, right;
387 if (num_loci<3 || window_size<3) return(TRUE); /* crash? */
388
389 for (i=0; i<num_loci-2; i++)
390 for (k=i+2; k<num_loci && k<i+window_size; k++)
391 for (j=i+1; j<k; j++)
392 if (is_excluded(locus[i],locus[j],locus[k])) return(FALSE);
393
394 return(TRUE);
395 }
396
397
398 /**************** N-Pt Maker Placer ****************/
399
place_locus(map,locus,start,finish,excluded,result,best_pos,best_map,temp)400 void place_locus(map,locus,start,finish,excluded,result,best_pos,best_map,temp)
401 MAP *map; /* map to place locus into, should be ctm-ed! */
402 int locus;
403 int start, finish; /* The leftmost and rightmost frame markers (i) to use */
404 bool *excluded; /* [interval#], one must be FALSE, NOT set */
405 PLACE **result; /* [interval#]->foo side-effected like=NO_LIKE if untried */
406 int *best_pos; /* side effected */
407 MAP *best_map, *temp; /* max_loci must be >=finish-start+2 */
408 {
409 int i, j, k, num_allowed, last, *zero=NULL;
410 real best_like, lod, theta;
411 real theta1, theta2, dist1, dist2, distframe, distscale;
412
413 /* Generate N-point maps from left to right markers with locus
414 inserted into each of the allowed positions, even one */
415
416 for (j=0, num_allowed=0; j<=map->num_loci; j++) {
417 if (!excluded[j]) num_allowed++;
418 result[j]->like=NO_LIKE; result[j]->dist=NO_THETA;
419 result[j]->net_error= 0.0; result[j]->worst_error= 0.0;
420 result[j]->zero=FALSE;
421 }
422
423 if (sex_specific || num_allowed==0 || best_map->max_loci<finish-start+2 ||
424 temp->max_loci<finish-start+2) send(CRASH);
425
426 last=finish+1;
427 best_like=VERY_UNLIKELY;
428 for (i=start; i<=last; i++) /* i=the interval # to place in */
429 if (!excluded[i]) {
430 clean_map(temp); k=0;
431 if (i==0) temp->locus[k++]=locus;
432 for (j=start; j<=finish; j++) { /* locus #s */
433 temp->locus[k++]=map->locus[j];
434 if (j+1==i) temp->locus[k++]=locus;
435 }
436 temp->num_loci= k;
437 init_rec_fracs(temp);
438 converge_to_map(temp);
439 result[i]->like=temp->log_like;
440
441 if (temp->log_like>best_like) { /* keep the best */
442 best_like=temp->log_like; *best_pos=i;
443 mapcpy(best_map,temp,FALSE);
444 }
445
446 /* normalize placement dists by any expansion of this interval
447 if i=first or last interval in the FRAMEWORK, then
448 dist= the correct rec frac, otherwise the placement dist =
449 DL*(Dframe/(DL+DR)) */
450 if (i==start) {
451 result[i]->dist= theta= temp->rec_frac[0][0];
452 if (theta<ZERO_DIST) result[i]->zero=TRUE;
453 } else if (i==last) {
454 result[i]->dist= theta= temp->rec_frac[i-start-1][MALE];
455 if (theta<ZERO_DIST) result[i]->zero=TRUE;
456 } else {
457 theta1= temp->rec_frac[i-start-1][MALE];
458 theta2= temp->rec_frac[i-start][MALE];
459 if (theta1<ZERO_DIST || theta2<ZERO_DIST) result[i]->zero=TRUE;
460 dist1= (*mapfunction->rec_to_dist)(theta1);
461 dist2= (*mapfunction->rec_to_dist)(theta2);
462 distframe= (*mapfunction->rec_to_dist)(map->rec_frac[i-1][0]);
463 distscale= dist1*(distframe/(dist1+dist2));
464 result[i]->dist= (*mapfunction->dist_to_rec)(distscale);
465 }
466
467 /* determine the error lods for any particular placement.
468 Here is a stupid algorithm which looks only at this locus's
469 lods. we realy need to look at flankers too! */
470 result[i]->net_error=0.0; result[i]->worst_error=0.0;
471 if (temp->allow_errors && (i>start || i<last) &&
472 temp->error_rate[i-start]>0.0) {
473 for (j=0; j<raw.data.f2.num_indivs; j++) { /* must be F2 */
474 lod= temp->error_lod[i-start][j];
475 if (lod>result[i]->worst_error) result[i]->worst_error=lod;
476 if (lod>=error_lod_thresh) result[i]->net_error+=lod;
477 }
478 }
479 } /* next interval i */
480 if (best_like==VERY_UNLIKELY) send(CRASH);
481
482 /* normalize the likes */
483 for (i=start; i<=last; i++) if (!excluded[i]) {
484 if (i==*best_pos) result[i]->like=0.0;
485 else result[i]->like-=best_like;
486 }
487
488 /* if all placements with nearly 0 like are zero dist from a framework
489 locus, remove all but the best (thus place_this will place the marker
490 uniquely) */
491 /* here we have a simple (and likely not complete enough) test...
492 are the likelihoods in consecutive intervals really zero,
493 and does each placement have SOME zero dist. */
494 if (result[*best_pos]->zero) {
495 for (i=*best_pos+1; !excluded[i] && i<=last &&
496 result[i]->like>ZERO_PLACE && result[i]->zero; i++)
497 result[i]->like=ZERO_LIKE;
498 for (i=*best_pos-1; !excluded[i] && i>=start &&
499 result[i]->like>ZERO_PLACE && result[i]->zero; i--)
500 result[i]->like=ZERO_LIKE;
501 }
502 }
503
504
npt_exclusions(place,num_loci,like_thresh,one_err_thresh,net_err_thresh,excluded,zero_place,off_ends,error_lod,single_error,next_best,best_unallowed)505 int npt_exclusions(place,num_loci,like_thresh,one_err_thresh,net_err_thresh,
506 excluded,zero_place,off_ends,error_lod,single_error,
507 next_best,best_unallowed)
508 /* return #ints ok */
509 PLACE **place; /* results from above */
510 int num_loci; /* in order */
511 real like_thresh, one_err_thresh, net_err_thresh;
512 bool *excluded; /* contents side-effected */
513 bool *zero_place; /* side-effected, T if all acceptable places are zeros */
514 int *off_ends; /* side-effected, >0 if all acceptable places are end(s) */
515 real *error_lod; /* side-effected, is greatest for all acceptable places */
516 bool *single_error; /* side-effected, T if error_place && err due to summing */
517 real *next_best; /* side-effected, best acceptable but non_zero like (<0) */
518 real *best_unallowed; /* side-effected, best like >threshold */
519 {
520 int i, best_i, num, zero, one, ends;
521 real second, unallowed, this_like, err, x;
522 if (like_thresh>=0.0) send(CRASH);
523
524 num=0; zero=TRUE; err=NO_ERRORS; one=FALSE, ends=2; best_i= -1;
525 second=NO_LIKE; unallowed=NO_LIKE; /* may stay this way */
526
527 for (i=0; i<=num_loci; i++) {
528 this_like=place[i]->like;
529 if (this_like!=NO_LIKE && this_like!=ZERO_LIKE &&
530 this_like>like_thresh) { /* it's acceptable */
531 num++;
532 if (excluded!=NULL) excluded[i]=FALSE;
533 if (this_like==0.0 && best_i==-1) best_i=i;
534 else if (this_like>second) second=this_like;
535 if ((x=place[i]->worst_error)>=one_err_thresh) { err=x; one=TRUE; }
536 else if ((x=place[i]->net_error)>=net_err_thresh) err=x;
537
538 } else { /* it's unacceptable */
539 if (excluded!=NULL) excluded[i]=TRUE;
540 if (i==0 || i==num_loci) ends--;
541 if (this_like!=NO_LIKE) {
542 if (this_like!=ZERO_LIKE) zero=FALSE;
543 if (this_like>unallowed) unallowed=this_like;
544 }
545 }
546 }
547 if (best_i==-1) send(CRASH);
548 if (place[best_i]->worst_error>=one_err_thresh) one=TRUE; /* override */
549
550 if (zero_place!=NULL) *zero_place= (num==1 && zero);
551 if (error_lod!=NULL) *error_lod=err;
552 if (single_error!=NULL) *single_error=one;
553 if (off_ends!=NULL) *off_ends=(ends>0 ? TRUE:FALSE);
554 if (next_best!=NULL) *next_best=second;
555 if (best_unallowed!=NULL) *best_unallowed= unallowed;
556 return(num);
557 }
558
559
560 /**************** Automatic Sub-Order Finder and Stuff ****************/
561
562
informative_subset(locus,num_loci,min_infs,min_theta,codom_only,haplos,subset,num)563 void informative_subset(locus,num_loci,min_infs,min_theta,codom_only,haplos,
564 subset,num)
565 int *locus, num_loci;
566 int min_infs;
567 real min_theta;
568 bool codom_only, haplos;
569 int *subset, *num;
570 {
571 int n_infs, n_dom_obs, n_het_obs, type, i, j, a, b, delete, deleted;
572 int a_infs, b_infs;
573 real lod, theta;
574 *num=0;
575
576 /* filter for informativeness */
577 for (i=0; i<num_loci; i++) {
578 f2_genotype(locus[i],haplos,observations);
579 n_infs=f2_count_infs(&n_dom_obs,&n_het_obs,observations);
580 type=raw.data.f2.cross_type;
581
582 if (n_infs<min_infs) continue;
583 if (type==F3_SELF || type==F2_INTERCROSS) {
584 if (codom_only && n_dom_obs!=0) continue;
585 }
586 subset[(*num)++]=locus[i];
587 }
588 if (*num<2) return;
589
590 /* find zeros - a bit kludgey */
591 do {
592 deleted=FALSE;
593 for (i=0; i<*num-1; i++) {
594 for (j=i+1; j<*num; j++) {
595 a=subset[i]; b=subset[j];
596 get_two_pt(a,b,&lod,&theta);
597 if (theta<min_theta) {
598 /* delete least informative marker, to bad we don't save
599 the infomativeness values */
600 f2_genotype(a,haplos,observations);
601 a_infs=f2_count_infs(&n_dom_obs,&n_het_obs,observations);
602 f2_genotype(b,haplos,observations);
603 b_infs=f2_count_infs(&n_dom_obs,&n_het_obs,observations);
604 /* do something smarter with mix/dr markers */
605 if (a_infs>b_infs) delete=(coin_flip() ? i:j);
606 else if (a_infs>b_infs) delete=j;
607 else delete=i;
608 if (delete<*num-1) subset[delete]=subset[*num-1];
609 *num-=1; deleted=TRUE; break;
610 }
611 }
612 if (deleted) break;
613 }
614 } while (deleted);
615 }
616
617
618
619 #define irand(n) (int)(randnum()*(real)n)
620 #define SEED_FAILED \
621 "Failed to find a starting order at log-likelihood threshold %.2lf.\n"
622 #define SEARCH_INF \
623 "Searching for a unique starting order containing %d of %d informative loci...\n"
624 #define SEARCH_ALL \
625 "Searching for a starting order containing %d of all %d loci...\n"
626 #define MAX_SAFETY 100
627
find_seed_order(is_subset,locus,num_loci,size,max_tries,thresh,map,temp_map,temp)628 bool find_seed_order(is_subset,locus,num_loci,size,max_tries,thresh,map,
629 temp_map,temp)
630 bool is_subset;
631 int *locus, num_loci, size, max_tries;
632 real thresh;
633 MAP *map, *temp_map;
634 bool **temp; /* [num_loci][num_loci] */
635 {
636 int try, pick, match, safety, i, j, k;
637 /* size must be >=num_loci; if size==num_loci, there is only one choice.
638 if size==num_loci-1, there are only num_loci choices, else there are
639 many orders. if size>num_loci it's an error */
640
641 if (num_loci<3) send(CRASH);
642 if (size>num_loci) size=num_loci;
643 try=0;
644 temp_map->num_loci=size; map->num_loci=0;
645
646 if (is_subset) sf(ps,SEARCH_INF,size,num_loci);
647 else sf(ps,SEARCH_ALL,size,num_loci);
648 pr();
649
650 if (size==num_loci) {
651 for (i=0; i<num_loci; i++) temp_map->locus[i]=locus[i];
652 keep_user_amused("subset",try+1,0);
653 if (good_seed(map,temp_map,thresh)) return(TRUE);
654 sf(ps,SEED_FAILED,-thresh); pr();
655 return(FALSE);
656
657 } else if (size==num_loci-1) {
658 for (j=0; j<num_loci; j++) { /* locus to punt */
659 k=0;
660 for (i=0; i<num_loci; i++)
661 if (i!=j) temp_map->locus[k++]=locus[i];
662 keep_user_amused("subset",++try,0);
663 if (good_seed(map,temp_map,thresh)) return(TRUE);
664 if (try>max_tries) break;
665 }
666 sf(ps,SEED_FAILED,-thresh); pr();
667 return(FALSE);
668
669 } else /* size << num_loci */ {
670 do { /* try a subset */
671 safety=0;
672 do { /* find an untried subset */
673 /* first use temp as a flag to indicate loci in order */
674 for (i=0; i<num_loci; i++) temp[try][i]=FALSE;
675 for (i=0; i<size; i++) {
676 do pick=irand(num_loci); while (temp[try][pick]);
677 temp_map->locus[i]=locus[pick]; temp[try][pick]=TRUE;
678 }
679 /* have we tried this subset before? */
680 sort_loci(temp_map->locus,size);
681 for (j=0, match=FALSE; j<try; j++) {
682 for (match=TRUE, k=0; k<size; k++)
683 if (temp_map->locus[k]!=temp[j][k]) { match=FALSE;break;}
684 if (match) break;
685 }
686 /* safety is a very lame way to see when we should give up,
687 because we have tried all possible subsets */
688 safety++;
689 } while (match && safety<=MAX_SAFETY);
690 if (match) break;
691 for (k=0; k<size; k++) temp[try][k]=temp_map->locus[k];
692 keep_user_amused("subset",++try,0);
693 /* for (i=0; i<size; i++)
694 { sf(ps,"%d ",temp_map->locus[i]+1); pr(); } nl(); */
695 if (good_seed(map,temp_map,thresh)) return(TRUE);
696 } while(try<=max_tries);
697 sf(ps,SEED_FAILED,-thresh); pr();
698 return(FALSE);
699 }
700 }
701
702
good_seed(map,temp_map,thresh)703 bool good_seed(map,temp_map,thresh)
704 MAP *map, *temp_map;
705 real thresh;
706 {
707 real best2=VERY_UNLIKELY, best=VERY_UNLIKELY;
708
709 make_compare_seq(temp_map->locus,temp_map->num_loci,0,temp_map->num_loci);
710 for_all_orders(seq,temp_map) {
711 if (use_three_pt &&
712 !three_pt_verify(temp_map->locus,temp_map->num_loci,
713 three_pt_window)) continue;
714 init_rec_fracs(temp_map);
715 converge_to_map(temp_map);
716 if (temp_map->log_like>best)
717 { best2=best; best=temp_map->log_like; mapcpy(map,temp_map,TRUE); }
718 else if (temp_map->log_like>best2) { best2=temp_map->log_like; }
719 }
720 if (best==VERY_UNLIKELY) return(FALSE);
721 else if (best2==VERY_UNLIKELY || best2-best<thresh) {
722 sf(ps,"Got one at log-likelihood %.2lf\n",best-best2); pr();
723 return(TRUE);
724 } else return(FALSE);
725 }
726
727
728
729 /**************** Automatic Order Maker ****************/
730
731 #define BY_NPT 1
732 #define BY_3PT 0
733 #define BY_NPT_ENDS (-1)
734 #define BY_NPT_ERRORS (-2)
735 #define BY_3PT_ENDS (-3)
736 #define NO_PLACES 9999
737
738 #define ALL_EXCLUDED \
739 "Three-point analysis excludes %d %smarker%s from all intervals\n"
740
741 #define sf_locnum(locus,paren) \
742 sf(ps,"%s%d%s%s",(paren ? "(":""),(locus)+1, \
743 (use_haplotypes && haplotyped(locus) ? "+":""),(paren ? ")":""))
744
extend_order(placed,unplaced,num_unplaced,npt_thresh,print_anyway)745 void extend_order(placed,unplaced,num_unplaced,npt_thresh,print_anyway)
746 MAP *placed; /* starting order, must be total long */
747 PLACEME **unplaced;
748 int *num_unplaced;
749 real npt_thresh; /* both are side-effected, and may have NO_LOCUS */
750 bool print_anyway;
751 {
752 int i, j, k, num, total;
753 bool placed_any, contradiction;
754 PLACE **placements=NULL;
755 MAP *temp_map=NULL;
756
757 run {
758 total=placed->num_loci + *num_unplaced;
759 /* if (*num_unplaced==0) return; let it fall through */
760 if (placed->num_loci<2) send(CRASH);
761
762 temp_map=allocate_map(total);
763 parray(placements,total,PLACE);
764 for (i=0; i<*num_unplaced; i++) {
765 unplaced[i]->best_pos=(-1);
766 unplaced[i]->num_places=0;
767 unplaced[i]->best_map->num_loci=0;
768 }
769
770 print("Start: ");
771 for (j=0; j<placed->num_loci; j++) {
772 sf_locnum(placed->locus[j],FALSE); pr();
773 if (j!=placed->num_loci-1) print(" ");
774 } nl();
775
776 placed_any=FALSE;
777 while (*num_unplaced>0)
778 if (!extend_order_by_1(placed,unplaced,num_unplaced,npt_thresh,
779 &contradiction,placements,temp_map)) {
780 sf(ps,"No unique placements for %d %smarker%s\n",
781 *num_unplaced,(placed_any ? "remaining ":""),
782 maybe_s(*num_unplaced)); pr();
783 if (contradiction==TRUE) {
784 sf(ps,"Three-point analysis allows no %s.\n",
785 (placed_any ? "further ordering":"orders at all")); pr();
786 } else if (contradiction==MAYBE) {
787 print("Maximum number of loci in one map reached.\n");
788 }
789 break;
790 } else placed_any=TRUE;
791 if (*num_unplaced==0)
792 { sf(ps,"Uniquely ordered all %d markers\n\n",total); pr(); }
793
794 if (print_anyway || placed_any) {
795 init_not_fixed(placed);
796 converge_to_map(placed);
797 nl();
798 print_long_map(placed,"Map:");
799 if (*num_unplaced>0) {
800 nl();
801 print("Markers placed relative to above map:\n");
802 new_print_placements(placed,unplaced,*num_unplaced);
803 }
804 }
805 } on_exit {
806 unparray(placements,total,PLACE);
807 free_map(temp_map);
808 relay_messages;
809 }
810 }
811
812
extend_order_by_1(placed,unplaced,num_unplaced,npt_thresh,contradiction,placements,temp_map)813 bool extend_order_by_1(placed,unplaced,num_unplaced,npt_thresh,contradiction,
814 placements,temp_map)
815 MAP *placed;
816 PLACEME **unplaced;
817 int *num_unplaced;
818 real npt_thresh;
819 bool *contradiction; /* iff return FALSE, is TRUE if 3pt excludes all */
820 PLACE **placements; /* these are all temps */
821 MAP *temp_map;
822 {
823 int places, best_places, best_pos, best_i, i, j;
824 int left, right, pos, how;
825 real next_best, best_unallowed, error_lod, e1, e2; /* e1, e2=threshs */
826 bool zero_place, off_ends, single_error;
827
828 if (placed->num_loci>=MAX_MAP_LOCI)
829 { *contradiction=MAYBE; return(FALSE); }
830 randomize_markers_to_try(unplaced,*num_unplaced);
831 for (i=0; i<*num_unplaced; i++) unplaced[i]->best_map->num_loci=0;
832
833 /* Try to place remaining loci, into order, in their order */
834 if (use_three_pt) {
835 best_places=NO_PLACES; best_i=-1;
836 for (i=0; i<*num_unplaced; i++) { /* sorted */
837 unplaced[i]->num_places= places=
838 three_pt_exclusions(placed->locus,placed->num_loci,
839 unplaced[i]->locus,unplaced[i]->excluded);
840 unplaced[i]->best_pos=(-1); /* NO_POS? */
841 if (places>0 && places<best_places) {
842 best_places=places; best_i=i;
843 if (places==1 && !middles_excluded(unplaced[i]->excluded,
844 placed->num_loci)) {
845 /* It doesn't get any better than this! */
846 add_to_order(BY_3PT,placed->locus,&placed->num_loci,
847 unplaced,best_i,num_unplaced);
848 return(TRUE);
849 }
850 }
851 } /* for i: loop over num_to_place trying 3pt */
852 if (best_places==NO_PLACES) { *contradiction=TRUE; return(FALSE); }
853
854 /**** We get to here if all allowed positions are off-end, or if no
855 unique placements were found ****/
856 #ifdef DONT_DO_THIS
857 for (i=0; i<*num_unplaced; i++) if (unplaced[i]->num_places==1) {
858 /* a troublesome however unique placement */
859 add_to_order(BY_3PT_ENDS,placed->locus,&placed->num_loci,unplaced,
860 i,num_unplaced);
861 return(TRUE);
862 }
863 #endif
864 /* otherwise fall-through to Npt tests below */
865
866 } else { /* !use_three-pt */
867 for (i=0; i<*num_unplaced; i++) {
868 /* unplaced[i]->num_places= placed->num_loci+1; */
869 for (j=0; j<placed->num_loci; j++) unplaced[i]->excluded[j]=FALSE;
870 }
871 }
872
873 /**** Have No UNIQUE Placements, so try N-pt analysis ****/
874 for (i=0; i<placed->num_loci-1; i++)
875 placed->rec_frac[i][MALE]=placed->rec_frac[i][FEMALE]=NOT_FIXED;
876 init_rec_fracs(placed); converge_to_map(placed);
877 rank_markers_to_try(unplaced,*num_unplaced);
878
879 best_places=NO_PLACES; best_i=-1;
880 for (i=0; i<*num_unplaced; i++) { /* sorted */
881 keep_user_amused("marker",i+1,*num_unplaced);
882 if (unplaced[i]->num_places==0) continue;
883 find_window(placed->locus,placed->num_loci,
884 unplaced[i]->locus,unplaced[i]->excluded,
885 npt_window,&left,&right);
886 place_locus(placed,unplaced[i]->locus,left,right,
887 unplaced[i]->excluded,placements,&unplaced[i]->best_pos,
888 unplaced[i]->best_map,temp_map);
889 unplaced[i]->num_places= places=
890 npt_exclusions(placements,placed->num_loci,npt_thresh,
891 error_single_thresh,error_net_thresh, /* GLOBALS */
892 unplaced[i]->excluded,&zero_place,&off_ends,
893 &error_lod,&single_error,&next_best,&best_unallowed);
894 unplaced[i]->off_end= off_ends;
895 unplaced[i]->errors= (error_lod!=NO_ERRORS);
896
897 /* preferentially add loci which place well */
898 if (places<best_places && !unplaced[i]->off_end &&
899 !unplaced[i]->errors) {
900 best_places=places; best_i=i;
901 if (places==1) { /* Greedy: It doesn't get any better than this! */
902 add_to_order(i+1,placed->locus,&placed->num_loci,unplaced,
903 best_i,num_unplaced);
904 return(TRUE);
905 }
906 }
907 } /* for i: loop over num_to_place trying 3pt */
908
909 /**** We get to here if all placements are troublesome, or if no unique
910 placements were found ****/
911 for (i=0; i<*num_unplaced; i++) if (unplaced[i]->num_places==1) {
912 /* a troublesome however unique placement */
913 how= (unplaced[i]->off_end ? BY_NPT_ENDS:BY_NPT_ERRORS);
914 add_to_order(how,placed->locus,&placed->num_loci,unplaced,i,
915 num_unplaced);
916 return(TRUE);
917 }
918
919 /**** Else we have no unique placements, we fail ****/
920 *contradiction=FALSE; return(FALSE);
921 }
922
923
randomize_markers_to_try(unplaced,num_unplaced)924 void randomize_markers_to_try(unplaced,num_unplaced)
925 PLACEME **unplaced;
926 int num_unplaced;
927 {
928 int i, n_dom, n_het, num;
929
930 for (i=0; i<num_unplaced; i++) {
931 unplaced[i]->num_places=0;
932 f2_genotype(unplaced[i]->locus,use_haplotypes,observations);
933 num=f2_count_infs(&n_dom,&n_het,observations);
934 unplaced[i]->priority=
935 imaxf(num-2*n_dom,0) +
936 irand(raw.data.f2.num_indivs/10);
937 }
938 qsort((QSORT_DATA_PTR_TO*)unplaced,(QSORT_LENGTH)num_unplaced,
939 sizeof(PLACEME*),compare_markers_to_try);
940 }
941
942
rank_markers_to_try(unplaced,num_unplaced)943 void rank_markers_to_try(unplaced,num_unplaced)
944 PLACEME **unplaced;
945 int num_unplaced;
946 {
947 qsort((QSORT_DATA_PTR_TO*)unplaced,(QSORT_LENGTH)num_unplaced,
948 sizeof(PLACEME*),compare_markers_to_try);
949 }
950
951
952 int compare_markers_to_try(a,b)
953 QSORT_COMPARE_PTR_TO(PLACEME*) *a;
954 QSORT_COMPARE_PTR_TO(PLACEME*) *b;
955 {
956 PLACEME *x, *y;
957
958 x=*a; y=*b;
959 /* num_places is three_pt criteria */
960 if (x->num_places < y->num_places) return(-1);
961 else if (x->num_places > y->num_places) return(1);
962 else if (x->priority > y->priority) return(-1);
963 else if (x->priority < y->priority) return(1);
964 else return(0);
965 }
966
967
add_to_order(how,order,num_ordered,unplaced,index,num_unplaced)968 void add_to_order(how,order,num_ordered,unplaced,index,num_unplaced)
969 int how, *order, *num_ordered;
970 PLACEME **unplaced;
971 int index, *num_unplaced;
972 {
973 int i, pos, last, new_marker;
974 PLACEME *temp;
975 new_marker=unplaced[index]->locus;
976
977 for (pos=0; pos<=*num_ordered; pos++)
978 if (!unplaced[index]->excluded[pos]) break;
979 if (pos==*num_ordered+1) send(CRASH);
980
981 new_marker=unplaced[index]->locus;
982 if (pos<*num_ordered)
983 for (i=*num_ordered-1; i>=pos; i--) order[i+1]=order[i];
984 order[pos]=new_marker;
985 *num_ordered+=1;
986
987 last=*num_unplaced-1;
988 if (*num_unplaced==1) *num_unplaced=0;
989 else if (index==last) *num_unplaced-=1;
990 else {
991 temp=unplaced[index];
992 unplaced[index]=unplaced[last];
993 unplaced[last]=temp;
994 *num_unplaced-=1;
995 }
996
997 #ifndef DEBUGGING
998 if (how>0) { sf(ps,"Npt-%d:\t ",how); pr(); }
999 #else
1000 if (how>0) { sf(ps,"Npt:\t "); pr(); }
1001 #endif
1002 else if (how==BY_3PT) print("3pt: ");
1003 else if (how==BY_NPT_ENDS) print("Npt-End: ");
1004 else if (how==BY_NPT_ERRORS) print("Npt-Err: ");
1005 else if (how==BY_3PT_ENDS) print("3pt-End: ");
1006 else send(CRASH);
1007 to_column(9);
1008
1009 for (i=0; i<*num_ordered; i++) { /* print new order */
1010 sf_locnum(order[i],(order[i]==new_marker)); pr();
1011 if (i!=*num_ordered-1) print(" ");
1012 } nl();
1013 }
1014
1015
1016 /* to keep this all straight:
1017 orig frame: a - b - c - d - e - f - g
1018 orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6 num_placed=7
1019 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7
1020 new frame: b - c - d - e - f
1021 new frame index: 0 - 1 - 2 - 3 - 4 */
1022
find_window(order,num_placed,new_locus,excluded,min_window,start,end)1023 void find_window(order,num_placed,new_locus,excluded,min_window,start,end)
1024 int *order; /*NOTUSED*/
1025 int num_placed, new_locus, *excluded;
1026 int min_window; /* #of loci to left and right, should be odd, e.g. 5 */
1027 int *start, *end; /* side-effected with interval indecies */
1028 {
1029 int j, step;
1030 /* simple implementation for now - this should later use informativeness */
1031
1032 step= max((min_window-1)/2,1);
1033 for (j=0; j<=num_placed && excluded[j]; j++) {}
1034 *start= max(j-step,0);
1035 for (j=num_placed; j>0 && excluded[j]; j--) {}
1036 *end= min(j+step-1,num_placed-1); /* Is this right? */
1037 }
1038
1039
count_positions(excluded,num_order,first_place)1040 int count_positions(excluded,num_order,first_place)
1041 int *excluded, num_order, *first_place;
1042 {
1043 int j, num_places;
1044
1045 for (j=0, num_places=0; j<=num_order; j++) if (!excluded[j]) {
1046 if (num_places==0 && first_place!=NULL) *first_place=j;
1047 num_places++;
1048 }
1049 return(num_places);
1050 }
1051
1052
middles_excluded(excluded,num_order)1053 bool middles_excluded(excluded,num_order)
1054 int *excluded, num_order;
1055 {
1056 int j;
1057
1058 for (j=1; j<num_order; j++) if (!excluded[j]) return(FALSE);
1059 return(TRUE);
1060 }
1061
1062
1063 #ifdef OBSOLETE /************************************************************/
make_npt_maps(order,start,finish,newmarker,threshold,excluded,map,like)1064 void make_npt_maps(order,start,finish,newmarker,threshold,excluded,map,like)
1065 int *order, start, finish, newmarker;
1066 real threshold; /* threshold is negative */
1067 bool *excluded;
1068 MAP *map;
1069 real *like;
1070 {
1071 int i, j, k, best_i, maybe_zero, first;
1072 real best_like;
1073
1074 best_like= VERY_UNLIKELY;
1075 for (i=start; i<=finish+1; i++) /* the interval #s */
1076 if (!excluded[i]) {
1077 clean_map(map); k=0;
1078 if (i==0) map->locus[k++]=newmarker;
1079 for (j=start; j<=finish; j++) { /* locus #s */
1080 map->locus[k++]=order[j];
1081 if (j+1==i) map->locus[k++]=newmarker;
1082 }
1083 map->num_loci= k;
1084
1085 init_rec_fracs(map);
1086 converge_to_map(map);
1087 like[i]=map->log_like;
1088 if (like[i]>best_like) { best_like=like[i]; best_i=i; }
1089 }
1090
1091 for (i=start; i<=finish+1; i++) if (!excluded[i]) {
1092 if (like[i] < best_like+threshold) excluded[i]=TRUE;
1093 like[i]-= best_like; /* check this! */
1094 }
1095
1096 /* if all likes are zero just pick one - NEED ALSO TEST FOR DIST */
1097 maybe_zero=TRUE;
1098 for (i=start; i<=finish+1; i++) if (!excluded[i]) {
1099 if (like[i]<=ZERO_PLACE) maybe_zero=FALSE;
1100 }
1101 first=TRUE;
1102 if (maybe_zero) for (i=start; i<=finish+1; i++) if (!excluded[i]) {
1103 if (first) first=FALSE; else excluded[i]=TRUE;
1104 }
1105 }
1106 #endif /* OBSOLETE */
1107
1108
1109 #ifdef EXPERIMENTAL /*******************************************************/
1110
1111 /* to keep this all straight:
1112 orig frame: a - b - c - d - e - f - g
1113 orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6
1114 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7
1115 new frame: y - c - d - e - z
1116 new frame index: 0 - 1 - 2 - 3 - 4
1117 place_in: 0 - 1 - 2 - 3 */
1118
make_pair_maps(order,start,finish,locus,i1,i2,threshold,excluded,map,temp_order)1119 void make_pair_maps(order,start,finish,locus,i1,i2,threshold,excluded,
1120 map,temp_order)
1121 int *order, start, finish;
1122 int *locus, i1, i2; /* the locus array and indecies of the two to place */
1123 real threshold; /* threshold is negative */
1124 bool **excluded;
1125 MAP *map;
1126 int **temp_order; /* used as a temp 1,2=frame, 3..., n... w/places */
1127 {
1128 int j, k, n, n_in_order, n_start, n_orders;
1129 real best;
1130
1131 temp_order[1][0]=NO_LOCUS; /* add two dummy loci */
1132 for (k=1, j=start; j<=finish; j++) temp_order[1][k++]=order[j];
1133 new_order[1][k]=NO_LOCUS; n_start=k+1;
1134
1135 for (k=0, j=start; j<=finish+1; j++) place_in[k++]=!excluded[i1][j];
1136 for (k=0, j=start; j<=finish; j++) temp_order[2][k++]=order[j];
1137 n_in_order=k;
1138
1139 add_orders(temp_order[1],n_start,locus[i1],place_in,&temp_order[2],1,
1140 n_in_order,&temp_order[3],&n_orders);
1141
1142 ++n_in_order;
1143 for (k=0, j=start; j<=finish+1; j++) place_in[k++]=!excluded[i2][j];
1144 n= 3+n_orders;
1145
1146 add_orders(temp_order[1],n_start,locus[i2],place_in,&temp_order[3],
1147 n_orders,n_in_order,&temp_order[n],&n_orders);
1148
1149 for (i=0; i<n_orders; i++) {
1150 sf(ps," %d->\t",i+1); pr();
1151 for (j=0; j<n_in_order+1; j++)
1152 { sf(ps,"%d ",temp_order[n+i][j]); pr(); }
1153 nl();
1154 }
1155
1156 n_in_order-=1; /* ignore dummies on end */
1157 best= VERY_UNLIKELY;
1158 for (i=0; i<n_orders; i++) {
1159 clean_map(map);
1160 for (j=0; j<n_in_order; j++)
1161 map->locus[j]= temp_order[n+i][j+1]; /* watch dummies on end */
1162 map->num_loci= n_in_order;
1163 init_rec_fracs(map);
1164 converge_to_map(map);
1165 }
1166
1167
1168
1169
1170 void add_orders(frame,n_frame,new_marker,place_in,old_order,n_orders,n_loci,
1171 new_order,num_orders)
1172 int *frame; /* frame must have LEFT_DUMMY@0 and RIGHT_DUMMY@n-1 */
1173 int n_frame; /* thus n_frame is real frame len+2 */
1174 int new_marker;
1175 bool *place_in; /* [n]=>T/F, n is the interval right of frame[n], 0..n_f-1 */
1176 int **old_order; /* [order#][position#] */
1177 int n_orders; /* #old orders */
1178 int n_loci; /* in each old order */
1179 int **new_order; /* [order#][position#] */
1180 int *new_orders; /* ptr to int - side-effected */
1181 {
1182 int n, i, j, k, left, right;
1183
1184 k=0;
1185 for (n=0; n<n_orders; n++) {
1186 for (i=0; i<n_frame-1; i++) if (place_in[i]) { /* never off end */
1187 left=frame[i]; right=frame[i+1];
1188 /* j= position of left marker in old_order[n] */
1189 for (j=0; old_order[n][j]!=left; j++) {}
1190
1191 do { /* until j hits right marker in old_order[n] */
1192 add_order(old_order[n],n_loci,j,new_marker,new_order[k++]);
1193 } while (old_order[n][++j]!=right);
1194 }
1195 }
1196 *new_orders= k;
1197 }
1198
1199
1200 void add_order(old_order,n_loci,pos,new_marker,new_order)
1201 int *old_order, n_loci, pos, new_marker, *new_order;
1202 {
1203 int i, d;
1204
1205 for (i=0, d=0; i<n_loci; i++) {
1206 new_order[i+d]= old_order[i];
1207 if (i==pos) { new_order[i+1]= new_marker; d=1; }
1208 }
1209 }
1210
1211
1212 /* needs excluded, a little tweaking */
1213 int **count_orders(frame,n_frame,new_marker,place_in,old_order,n_orders,n_loci,
1214 num_orders)
1215 int *frame; /* frame must have LEFT_DUMMY@0 and RIGHT_DUMMY@n-1 */
1216 int n_frame; /* thus n_frame is real frame len+2 */
1217 int new_marker;
1218 bool *place_in; /* [n]=>T/F, n is the interval right of frame[n], 0..n_f-1 */
1219 int **old_order; /* [order#][position#] */
1220 int n_orders; /* #old orders */
1221 int n_loci; /* in each old order */
1222 int **new_order; /* [order#][position#] */
1223 int *new_orders; /* ptr to int - side-effected */
1224 {
1225 int n, i, j, k, left, right;
1226
1227 k=0;
1228 for (n=0; n<n_orders; n++) {
1229 for (i=0; i<n_frame-1; i++) if (place_in[i]) { /* never off end */
1230 left=frame[i]; right=frame[i+1];
1231 /* j= position of left in old_order[n] */
1232 for (j=0; old_order[n][j]!=left; j++) {}
1233
1234 do { /* until we hit right in old_order[n] */
1235 add_order(old_order[n],n_loci,j,new_marker,new_order[k++]);
1236 } while (old_order[n][++j]!=right);
1237 }
1238 }
1239 *new_orders= k;
1240 }
1241 #endif /* EXPERIMENTAL */
1242
1243