/****************************************************************************** #### ##### ##### ###### ##### #### #### # # # # # # # # # # # # # # # # # # ##### # # #### # # # ##### # # # ##### # ### # # # # # # # # # # # # ### # # #### # # ##### ###### # # #### ### #### ******************************************************************************/ /* This file is part of MAPMAKER 3.0b, Copyright 1987-1992, Whitehead Institute for Biomedical Research. All rights reserved. See READ.ME for license. */ #define INC_LIB #define INC_SHELL #define INC_MISC #include "mapm.h" #define set_excluded(a,b,c) \ three_pt_excluded[three_pt_index[a]][three_pt_index[b]][three_pt_index[c]] \ =EXCLUDED /* OBSOLETE? */ void left_exclusion(), right_exclusion(), middle_exclusion(); bool is_excluded(), is_computed(); void add_to_order(); int compare_markers_to_try(); void randomize_markers_to_try(); void rank_markers_to_try(); void add_order(); bool extend_order_by_1(); int **add_orders(); bool middles_excluded(); char ***three_pt_excluded=NULL; int three_pt_size=0; int *three_pt_index=NULL; #define EXCLUDED 0 #define ALLOWED 1 #define UNCOMPUTED 2 void alloc_3pt_matrix(); /* args: int num_loci; sets a global, will free and realloc if it's too small, use to prealloc for efficiency */ void free_3pt_matrix(); /**************** Two-point, Haplotypes ****************/ void get_linkage_group(locus,num_loci,linkage_group,group_size,lodbound, thetabound) int *locus, *num_loci, *linkage_group, *group_size; real lodbound,thetabound; { int i, j, linked, *unlinked_list, unlinked_size, oldsize, a, b; real lod, theta, thetam, thetaf; /* start linkage_group with first marker in the list of loci */ linkage_group[0]=locus[0]; *group_size= 1; unlinked_list= locus; *num_loci-= 1; unlinked_size= *num_loci; locus= &unlinked_list[1]; do { oldsize= unlinked_size; unlinked_size= 0; for (i=0; i=lodbound && theta<=thetabound) { linked=TRUE; break; } } else { get_sex_two_pt(a,b,&lod,&thetam,&thetaf); if (lod>lodbound && (thetam<=thetabound || thetaf<=thetabound)) { linked=TRUE; break; } } } if (linked) { linkage_group[*group_size]= a; ++*group_size; } else { /* Effectively, write back into the locus array, relying on the fact that unlinked_size<=*num_loci */ unlinked_list[unlinked_size]= a; ++unlinked_size; } } /* for i */ locus= unlinked_list; *num_loci= unlinked_size; } while (unlinked_size!=oldsize && unlinked_size!=0); } void find_haplo_group(locus,num_loci,haplo_group,num_haplo,old_obs,new_obs) int *locus, *num_loci, *haplo_group, *num_haplo; int *old_obs, *new_obs; { int i, num_inf, num_dom, best_i, foo; int best_score, old_score, score; /* get the most informative marker */ best_score= 0; for (i=0; i<*num_loci; i++) { f2_genotype(locus[i],FALSE,old_obs); num_inf= f2_count_infs(&num_dom,&foo,old_obs); score= 2*num_inf - num_dom; if (score>best_score) { best_score=score; best_i=i; } } /* delete it, and make it the haplo group */ haplo_group[0]=locus[best_i]; *num_haplo=1; f2_genotype(locus[best_i],FALSE,old_obs); locus[best_i]= locus[--*num_loci]; if (*num_loci==0) return; do { /* find one to add to the group */ old_score= best_score; best_score= 0; for (i=0; i<*num_loci; i++) if (merge_genotypes(locus[i],old_obs,new_obs)) { /* if no recs */ num_inf= f2_count_infs(&num_dom,&foo,new_obs); score= 2*num_inf - num_dom; if (score>best_score) { best_score=score; best_i=i; } } if (best_score !=0 ) { /* if we found any with no recs */ merge_genotypes(locus[best_i],old_obs,old_obs); /*updates old_obs*/ haplo_group[*num_haplo]=locus[best_i]; ++*num_haplo; /* add it */ locus[best_i]= locus[--*num_loci]; /* delete it */ } } while (*num_loci>0 && best_score!=0); /* got one and more to try */ } /**************** Three Point ****************/ void alloc_3pt_matrix(num_group) int num_group; { if (three_pt_excluded!=NULL && num_group<=three_pt_size) return; run { if (three_pt_excluded!=NULL) free_3pt_matrix(); three_pt_size=num_group; three_pt_excluded= alloc_char_3d_matrix(num_group,num_group,num_group); array(three_pt_index,raw.num_markers,int); } when_aborting { free_3pt_matrix(); relay_messages; } } void free_3pt_matrix() { free_char_3d_matrix(three_pt_excluded,three_pt_size,three_pt_size); unarray(three_pt_index,int); three_pt_excluded=NULL; three_pt_size=0; } void setup_3pt_data(locus,num_loci,threshold) int *locus, num_loci; real threshold; { int i, j, k; real d1, d2, d3; alloc_3pt_matrix(num_loci); /* three_pt_index[] is lookup by raw locus# to three_pt_exclusion index */ for (i=0; ithreshold) { best=0.0; continue;} starter[0]=i; starter[1]=j; starter[2]=k; } else if (d1threshold) { best=0.0; continue;} starter[0]=i; starter[1]=k; starter[2]=j; } else if (d1threshold) { best=0.0; continue;} starter[0]=j; starter[1]=i; starter[2]=k; } } } } if (best==0.0) return(FALSE); /* else, setup order list and delete starters from locus list */ order[0]=locus[starter[0]]; locus[starter[0]]=NO_LOCUS; order[1]=locus[starter[1]]; locus[starter[1]]=NO_LOCUS; order[2]=locus[starter[2]]; locus[starter[2]]=NO_LOCUS; *num_ordered= 3; num_left=0; for (i=0; i<*num_loci; i++) { temp=locus[i]; if (temp!=NO_LOCUS) locus[num_left++]=temp; } *num_loci= num_left; return(TRUE); } bool is_excluded(a,b,c) int a, b, c; { return(three_pt_excluded[three_pt_index[a]][three_pt_index[b]] [three_pt_index[c]]==EXCLUDED); } bool is_computed(a,b,c) int a, b, c; { return(three_pt_excluded[three_pt_index[a]][three_pt_index[b]] [three_pt_index[c]]!=UNCOMPUTED); } void left_exclusion(excluded,locus,num_loci,newmarker) bool *excluded; int *locus, num_loci, newmarker; { int middle, end, i; for (middle=num_loci-2; middle>=0; middle--) { for (end=num_loci-1; end>middle; end--) { if (is_excluded(newmarker,locus[middle],locus[end])) { for (i=0; i<=middle; i++) excluded[i]=TRUE; return; } } } } void right_exclusion(excluded,locus,num_loci,newmarker) bool *excluded; int *locus, num_loci, newmarker; { int middle,end,i; for (middle=1; middle<=num_loci-1; middle++) { for (end=0; endmiddle; i--) excluded[i]=TRUE; return; } } } } void middle_exclusion(excluded,locus,num_loci,leftmost,rightmost, init_window_size,newmarker) bool *excluded; int *locus, num_loci; int leftmost, rightmost; /* presumably, the left and rightmost intervals? */ int init_window_size, newmarker; { int first, last, start, end, left, width, i, j; for (width=init_window_size; width>=1; width--) { first=max(1,leftmost-width+1); last= min(rightmost,num_loci-width); for (i=first; i<=last; i++) { start=i; end=i+width-1; if (is_excluded(locus[start-1],newmarker,locus[end])) { for (j=start; j<=end; j++) excluded[j]=TRUE; if (start>leftmost && width>1) middle_exclusion(excluded,locus,num_loci,leftmost,start-1, width-1,newmarker); if (end==rightmost) return; leftmost=end+1; } } } } /* orig frame: a - b - c - d - e - f - g orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6 num_placed=7 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 */ int three_pt_exclusions(order,num_placed,newmarker,excluded) int *order, num_placed, newmarker; bool *excluded; /* side-effected */ { int i, j, left, right, num_places; for (i=0; ifoo side-effected like=NO_LIKE if untried */ int *best_pos; /* side effected */ MAP *best_map, *temp; /* max_loci must be >=finish-start+2 */ { int i, j, k, num_allowed, last, *zero=NULL; real best_like, lod, theta; real theta1, theta2, dist1, dist2, distframe, distscale; /* Generate N-point maps from left to right markers with locus inserted into each of the allowed positions, even one */ for (j=0, num_allowed=0; j<=map->num_loci; j++) { if (!excluded[j]) num_allowed++; result[j]->like=NO_LIKE; result[j]->dist=NO_THETA; result[j]->net_error= 0.0; result[j]->worst_error= 0.0; result[j]->zero=FALSE; } if (sex_specific || num_allowed==0 || best_map->max_locimax_locilocus[k++]=locus; for (j=start; j<=finish; j++) { /* locus #s */ temp->locus[k++]=map->locus[j]; if (j+1==i) temp->locus[k++]=locus; } temp->num_loci= k; init_rec_fracs(temp); converge_to_map(temp); result[i]->like=temp->log_like; if (temp->log_like>best_like) { /* keep the best */ best_like=temp->log_like; *best_pos=i; mapcpy(best_map,temp,FALSE); } /* normalize placement dists by any expansion of this interval if i=first or last interval in the FRAMEWORK, then dist= the correct rec frac, otherwise the placement dist = DL*(Dframe/(DL+DR)) */ if (i==start) { result[i]->dist= theta= temp->rec_frac[0][0]; if (thetazero=TRUE; } else if (i==last) { result[i]->dist= theta= temp->rec_frac[i-start-1][MALE]; if (thetazero=TRUE; } else { theta1= temp->rec_frac[i-start-1][MALE]; theta2= temp->rec_frac[i-start][MALE]; if (theta1zero=TRUE; dist1= (*mapfunction->rec_to_dist)(theta1); dist2= (*mapfunction->rec_to_dist)(theta2); distframe= (*mapfunction->rec_to_dist)(map->rec_frac[i-1][0]); distscale= dist1*(distframe/(dist1+dist2)); result[i]->dist= (*mapfunction->dist_to_rec)(distscale); } /* determine the error lods for any particular placement. Here is a stupid algorithm which looks only at this locus's lods. we realy need to look at flankers too! */ result[i]->net_error=0.0; result[i]->worst_error=0.0; if (temp->allow_errors && (i>start || ierror_rate[i-start]>0.0) { for (j=0; jerror_lod[i-start][j]; if (lod>result[i]->worst_error) result[i]->worst_error=lod; if (lod>=error_lod_thresh) result[i]->net_error+=lod; } } } /* next interval i */ if (best_like==VERY_UNLIKELY) send(CRASH); /* normalize the likes */ for (i=start; i<=last; i++) if (!excluded[i]) { if (i==*best_pos) result[i]->like=0.0; else result[i]->like-=best_like; } /* if all placements with nearly 0 like are zero dist from a framework locus, remove all but the best (thus place_this will place the marker uniquely) */ /* here we have a simple (and likely not complete enough) test... are the likelihoods in consecutive intervals really zero, and does each placement have SOME zero dist. */ if (result[*best_pos]->zero) { for (i=*best_pos+1; !excluded[i] && i<=last && result[i]->like>ZERO_PLACE && result[i]->zero; i++) result[i]->like=ZERO_LIKE; for (i=*best_pos-1; !excluded[i] && i>=start && result[i]->like>ZERO_PLACE && result[i]->zero; i--) result[i]->like=ZERO_LIKE; } } int 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) /* return #ints ok */ PLACE **place; /* results from above */ int num_loci; /* in order */ real like_thresh, one_err_thresh, net_err_thresh; bool *excluded; /* contents side-effected */ bool *zero_place; /* side-effected, T if all acceptable places are zeros */ int *off_ends; /* side-effected, >0 if all acceptable places are end(s) */ real *error_lod; /* side-effected, is greatest for all acceptable places */ bool *single_error; /* side-effected, T if error_place && err due to summing */ real *next_best; /* side-effected, best acceptable but non_zero like (<0) */ real *best_unallowed; /* side-effected, best like >threshold */ { int i, best_i, num, zero, one, ends; real second, unallowed, this_like, err, x; if (like_thresh>=0.0) send(CRASH); num=0; zero=TRUE; err=NO_ERRORS; one=FALSE, ends=2; best_i= -1; second=NO_LIKE; unallowed=NO_LIKE; /* may stay this way */ for (i=0; i<=num_loci; i++) { this_like=place[i]->like; if (this_like!=NO_LIKE && this_like!=ZERO_LIKE && this_like>like_thresh) { /* it's acceptable */ num++; if (excluded!=NULL) excluded[i]=FALSE; if (this_like==0.0 && best_i==-1) best_i=i; else if (this_like>second) second=this_like; if ((x=place[i]->worst_error)>=one_err_thresh) { err=x; one=TRUE; } else if ((x=place[i]->net_error)>=net_err_thresh) err=x; } else { /* it's unacceptable */ if (excluded!=NULL) excluded[i]=TRUE; if (i==0 || i==num_loci) ends--; if (this_like!=NO_LIKE) { if (this_like!=ZERO_LIKE) zero=FALSE; if (this_like>unallowed) unallowed=this_like; } } } if (best_i==-1) send(CRASH); if (place[best_i]->worst_error>=one_err_thresh) one=TRUE; /* override */ if (zero_place!=NULL) *zero_place= (num==1 && zero); if (error_lod!=NULL) *error_lod=err; if (single_error!=NULL) *single_error=one; if (off_ends!=NULL) *off_ends=(ends>0 ? TRUE:FALSE); if (next_best!=NULL) *next_best=second; if (best_unallowed!=NULL) *best_unallowed= unallowed; return(num); } /**************** Automatic Sub-Order Finder and Stuff ****************/ void informative_subset(locus,num_loci,min_infs,min_theta,codom_only,haplos, subset,num) int *locus, num_loci; int min_infs; real min_theta; bool codom_only, haplos; int *subset, *num; { int n_infs, n_dom_obs, n_het_obs, type, i, j, a, b, delete, deleted; int a_infs, b_infs; real lod, theta; *num=0; /* filter for informativeness */ for (i=0; ib_infs) delete=(coin_flip() ? i:j); else if (a_infs>b_infs) delete=j; else delete=i; if (delete<*num-1) subset[delete]=subset[*num-1]; *num-=1; deleted=TRUE; break; } } if (deleted) break; } } while (deleted); } #define irand(n) (int)(randnum()*(real)n) #define SEED_FAILED \ "Failed to find a starting order at log-likelihood threshold %.2lf.\n" #define SEARCH_INF \ "Searching for a unique starting order containing %d of %d informative loci...\n" #define SEARCH_ALL \ "Searching for a starting order containing %d of all %d loci...\n" #define MAX_SAFETY 100 bool find_seed_order(is_subset,locus,num_loci,size,max_tries,thresh,map, temp_map,temp) bool is_subset; int *locus, num_loci, size, max_tries; real thresh; MAP *map, *temp_map; bool **temp; /* [num_loci][num_loci] */ { int try, pick, match, safety, i, j, k; /* size must be >=num_loci; if size==num_loci, there is only one choice. if size==num_loci-1, there are only num_loci choices, else there are many orders. if size>num_loci it's an error */ if (num_loci<3) send(CRASH); if (size>num_loci) size=num_loci; try=0; temp_map->num_loci=size; map->num_loci=0; if (is_subset) sf(ps,SEARCH_INF,size,num_loci); else sf(ps,SEARCH_ALL,size,num_loci); pr(); if (size==num_loci) { for (i=0; ilocus[i]=locus[i]; keep_user_amused("subset",try+1,0); if (good_seed(map,temp_map,thresh)) return(TRUE); sf(ps,SEED_FAILED,-thresh); pr(); return(FALSE); } else if (size==num_loci-1) { for (j=0; jlocus[k++]=locus[i]; keep_user_amused("subset",++try,0); if (good_seed(map,temp_map,thresh)) return(TRUE); if (try>max_tries) break; } sf(ps,SEED_FAILED,-thresh); pr(); return(FALSE); } else /* size << num_loci */ { do { /* try a subset */ safety=0; do { /* find an untried subset */ /* first use temp as a flag to indicate loci in order */ for (i=0; ilocus[i]=locus[pick]; temp[try][pick]=TRUE; } /* have we tried this subset before? */ sort_loci(temp_map->locus,size); for (j=0, match=FALSE; jlocus[k]!=temp[j][k]) { match=FALSE;break;} if (match) break; } /* safety is a very lame way to see when we should give up, because we have tried all possible subsets */ safety++; } while (match && safety<=MAX_SAFETY); if (match) break; for (k=0; klocus[k]; keep_user_amused("subset",++try,0); /* for (i=0; ilocus[i]+1); pr(); } nl(); */ if (good_seed(map,temp_map,thresh)) return(TRUE); } while(try<=max_tries); sf(ps,SEED_FAILED,-thresh); pr(); return(FALSE); } } bool good_seed(map,temp_map,thresh) MAP *map, *temp_map; real thresh; { real best2=VERY_UNLIKELY, best=VERY_UNLIKELY; make_compare_seq(temp_map->locus,temp_map->num_loci,0,temp_map->num_loci); for_all_orders(seq,temp_map) { if (use_three_pt && !three_pt_verify(temp_map->locus,temp_map->num_loci, three_pt_window)) continue; init_rec_fracs(temp_map); converge_to_map(temp_map); if (temp_map->log_like>best) { best2=best; best=temp_map->log_like; mapcpy(map,temp_map,TRUE); } else if (temp_map->log_like>best2) { best2=temp_map->log_like; } } if (best==VERY_UNLIKELY) return(FALSE); else if (best2==VERY_UNLIKELY || best2-bestnum_loci + *num_unplaced; /* if (*num_unplaced==0) return; let it fall through */ if (placed->num_loci<2) send(CRASH); temp_map=allocate_map(total); parray(placements,total,PLACE); for (i=0; i<*num_unplaced; i++) { unplaced[i]->best_pos=(-1); unplaced[i]->num_places=0; unplaced[i]->best_map->num_loci=0; } print("Start: "); for (j=0; jnum_loci; j++) { sf_locnum(placed->locus[j],FALSE); pr(); if (j!=placed->num_loci-1) print(" "); } nl(); placed_any=FALSE; while (*num_unplaced>0) if (!extend_order_by_1(placed,unplaced,num_unplaced,npt_thresh, &contradiction,placements,temp_map)) { sf(ps,"No unique placements for %d %smarker%s\n", *num_unplaced,(placed_any ? "remaining ":""), maybe_s(*num_unplaced)); pr(); if (contradiction==TRUE) { sf(ps,"Three-point analysis allows no %s.\n", (placed_any ? "further ordering":"orders at all")); pr(); } else if (contradiction==MAYBE) { print("Maximum number of loci in one map reached.\n"); } break; } else placed_any=TRUE; if (*num_unplaced==0) { sf(ps,"Uniquely ordered all %d markers\n\n",total); pr(); } if (print_anyway || placed_any) { init_not_fixed(placed); converge_to_map(placed); nl(); print_long_map(placed,"Map:"); if (*num_unplaced>0) { nl(); print("Markers placed relative to above map:\n"); new_print_placements(placed,unplaced,*num_unplaced); } } } on_exit { unparray(placements,total,PLACE); free_map(temp_map); relay_messages; } } bool extend_order_by_1(placed,unplaced,num_unplaced,npt_thresh,contradiction, placements,temp_map) MAP *placed; PLACEME **unplaced; int *num_unplaced; real npt_thresh; bool *contradiction; /* iff return FALSE, is TRUE if 3pt excludes all */ PLACE **placements; /* these are all temps */ MAP *temp_map; { int places, best_places, best_pos, best_i, i, j; int left, right, pos, how; real next_best, best_unallowed, error_lod, e1, e2; /* e1, e2=threshs */ bool zero_place, off_ends, single_error; if (placed->num_loci>=MAX_MAP_LOCI) { *contradiction=MAYBE; return(FALSE); } randomize_markers_to_try(unplaced,*num_unplaced); for (i=0; i<*num_unplaced; i++) unplaced[i]->best_map->num_loci=0; /* Try to place remaining loci, into order, in their order */ if (use_three_pt) { best_places=NO_PLACES; best_i=-1; for (i=0; i<*num_unplaced; i++) { /* sorted */ unplaced[i]->num_places= places= three_pt_exclusions(placed->locus,placed->num_loci, unplaced[i]->locus,unplaced[i]->excluded); unplaced[i]->best_pos=(-1); /* NO_POS? */ if (places>0 && placesexcluded, placed->num_loci)) { /* It doesn't get any better than this! */ add_to_order(BY_3PT,placed->locus,&placed->num_loci, unplaced,best_i,num_unplaced); return(TRUE); } } } /* for i: loop over num_to_place trying 3pt */ if (best_places==NO_PLACES) { *contradiction=TRUE; return(FALSE); } /**** We get to here if all allowed positions are off-end, or if no unique placements were found ****/ #ifdef DONT_DO_THIS for (i=0; i<*num_unplaced; i++) if (unplaced[i]->num_places==1) { /* a troublesome however unique placement */ add_to_order(BY_3PT_ENDS,placed->locus,&placed->num_loci,unplaced, i,num_unplaced); return(TRUE); } #endif /* otherwise fall-through to Npt tests below */ } else { /* !use_three-pt */ for (i=0; i<*num_unplaced; i++) { /* unplaced[i]->num_places= placed->num_loci+1; */ for (j=0; jnum_loci; j++) unplaced[i]->excluded[j]=FALSE; } } /**** Have No UNIQUE Placements, so try N-pt analysis ****/ for (i=0; inum_loci-1; i++) placed->rec_frac[i][MALE]=placed->rec_frac[i][FEMALE]=NOT_FIXED; init_rec_fracs(placed); converge_to_map(placed); rank_markers_to_try(unplaced,*num_unplaced); best_places=NO_PLACES; best_i=-1; for (i=0; i<*num_unplaced; i++) { /* sorted */ keep_user_amused("marker",i+1,*num_unplaced); if (unplaced[i]->num_places==0) continue; find_window(placed->locus,placed->num_loci, unplaced[i]->locus,unplaced[i]->excluded, npt_window,&left,&right); place_locus(placed,unplaced[i]->locus,left,right, unplaced[i]->excluded,placements,&unplaced[i]->best_pos, unplaced[i]->best_map,temp_map); unplaced[i]->num_places= places= npt_exclusions(placements,placed->num_loci,npt_thresh, error_single_thresh,error_net_thresh, /* GLOBALS */ unplaced[i]->excluded,&zero_place,&off_ends, &error_lod,&single_error,&next_best,&best_unallowed); unplaced[i]->off_end= off_ends; unplaced[i]->errors= (error_lod!=NO_ERRORS); /* preferentially add loci which place well */ if (placesoff_end && !unplaced[i]->errors) { best_places=places; best_i=i; if (places==1) { /* Greedy: It doesn't get any better than this! */ add_to_order(i+1,placed->locus,&placed->num_loci,unplaced, best_i,num_unplaced); return(TRUE); } } } /* for i: loop over num_to_place trying 3pt */ /**** We get to here if all placements are troublesome, or if no unique placements were found ****/ for (i=0; i<*num_unplaced; i++) if (unplaced[i]->num_places==1) { /* a troublesome however unique placement */ how= (unplaced[i]->off_end ? BY_NPT_ENDS:BY_NPT_ERRORS); add_to_order(how,placed->locus,&placed->num_loci,unplaced,i, num_unplaced); return(TRUE); } /**** Else we have no unique placements, we fail ****/ *contradiction=FALSE; return(FALSE); } void randomize_markers_to_try(unplaced,num_unplaced) PLACEME **unplaced; int num_unplaced; { int i, n_dom, n_het, num; for (i=0; inum_places=0; f2_genotype(unplaced[i]->locus,use_haplotypes,observations); num=f2_count_infs(&n_dom,&n_het,observations); unplaced[i]->priority= imaxf(num-2*n_dom,0) + irand(raw.data.f2.num_indivs/10); } qsort((QSORT_DATA_PTR_TO*)unplaced,(QSORT_LENGTH)num_unplaced, sizeof(PLACEME*),compare_markers_to_try); } void rank_markers_to_try(unplaced,num_unplaced) PLACEME **unplaced; int num_unplaced; { qsort((QSORT_DATA_PTR_TO*)unplaced,(QSORT_LENGTH)num_unplaced, sizeof(PLACEME*),compare_markers_to_try); } int compare_markers_to_try(a,b) QSORT_COMPARE_PTR_TO(PLACEME*) *a; QSORT_COMPARE_PTR_TO(PLACEME*) *b; { PLACEME *x, *y; x=*a; y=*b; /* num_places is three_pt criteria */ if (x->num_places < y->num_places) return(-1); else if (x->num_places > y->num_places) return(1); else if (x->priority > y->priority) return(-1); else if (x->priority < y->priority) return(1); else return(0); } void add_to_order(how,order,num_ordered,unplaced,index,num_unplaced) int how, *order, *num_ordered; PLACEME **unplaced; int index, *num_unplaced; { int i, pos, last, new_marker; PLACEME *temp; new_marker=unplaced[index]->locus; for (pos=0; pos<=*num_ordered; pos++) if (!unplaced[index]->excluded[pos]) break; if (pos==*num_ordered+1) send(CRASH); new_marker=unplaced[index]->locus; if (pos<*num_ordered) for (i=*num_ordered-1; i>=pos; i--) order[i+1]=order[i]; order[pos]=new_marker; *num_ordered+=1; last=*num_unplaced-1; if (*num_unplaced==1) *num_unplaced=0; else if (index==last) *num_unplaced-=1; else { temp=unplaced[index]; unplaced[index]=unplaced[last]; unplaced[last]=temp; *num_unplaced-=1; } #ifndef DEBUGGING if (how>0) { sf(ps,"Npt-%d:\t ",how); pr(); } #else if (how>0) { sf(ps,"Npt:\t "); pr(); } #endif else if (how==BY_3PT) print("3pt: "); else if (how==BY_NPT_ENDS) print("Npt-End: "); else if (how==BY_NPT_ERRORS) print("Npt-Err: "); else if (how==BY_3PT_ENDS) print("3pt-End: "); else send(CRASH); to_column(9); for (i=0; i<*num_ordered; i++) { /* print new order */ sf_locnum(order[i],(order[i]==new_marker)); pr(); if (i!=*num_ordered-1) print(" "); } nl(); } /* to keep this all straight: orig frame: a - b - c - d - e - f - g orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6 num_placed=7 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 new frame: b - c - d - e - f new frame index: 0 - 1 - 2 - 3 - 4 */ void find_window(order,num_placed,new_locus,excluded,min_window,start,end) int *order; /*NOTUSED*/ int num_placed, new_locus, *excluded; int min_window; /* #of loci to left and right, should be odd, e.g. 5 */ int *start, *end; /* side-effected with interval indecies */ { int j, step; /* simple implementation for now - this should later use informativeness */ step= max((min_window-1)/2,1); for (j=0; j<=num_placed && excluded[j]; j++) {} *start= max(j-step,0); for (j=num_placed; j>0 && excluded[j]; j--) {} *end= min(j+step-1,num_placed-1); /* Is this right? */ } int count_positions(excluded,num_order,first_place) int *excluded, num_order, *first_place; { int j, num_places; for (j=0, num_places=0; j<=num_order; j++) if (!excluded[j]) { if (num_places==0 && first_place!=NULL) *first_place=j; num_places++; } return(num_places); } bool middles_excluded(excluded,num_order) int *excluded, num_order; { int j; for (j=1; jlocus[k++]=newmarker; for (j=start; j<=finish; j++) { /* locus #s */ map->locus[k++]=order[j]; if (j+1==i) map->locus[k++]=newmarker; } map->num_loci= k; init_rec_fracs(map); converge_to_map(map); like[i]=map->log_like; if (like[i]>best_like) { best_like=like[i]; best_i=i; } } for (i=start; i<=finish+1; i++) if (!excluded[i]) { if (like[i] < best_like+threshold) excluded[i]=TRUE; like[i]-= best_like; /* check this! */ } /* if all likes are zero just pick one - NEED ALSO TEST FOR DIST */ maybe_zero=TRUE; for (i=start; i<=finish+1; i++) if (!excluded[i]) { if (like[i]<=ZERO_PLACE) maybe_zero=FALSE; } first=TRUE; if (maybe_zero) for (i=start; i<=finish+1; i++) if (!excluded[i]) { if (first) first=FALSE; else excluded[i]=TRUE; } } #endif /* OBSOLETE */ #ifdef EXPERIMENTAL /*******************************************************/ /* to keep this all straight: orig frame: a - b - c - d - e - f - g orig frame index: 0 - 1 - 2 - 3 - 4 - 5 - 6 exclusion index: 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 new frame: y - c - d - e - z new frame index: 0 - 1 - 2 - 3 - 4 place_in: 0 - 1 - 2 - 3 */ void make_pair_maps(order,start,finish,locus,i1,i2,threshold,excluded, map,temp_order) int *order, start, finish; int *locus, i1, i2; /* the locus array and indecies of the two to place */ real threshold; /* threshold is negative */ bool **excluded; MAP *map; int **temp_order; /* used as a temp 1,2=frame, 3..., n... w/places */ { int j, k, n, n_in_order, n_start, n_orders; real best; temp_order[1][0]=NO_LOCUS; /* add two dummy loci */ for (k=1, j=start; j<=finish; j++) temp_order[1][k++]=order[j]; new_order[1][k]=NO_LOCUS; n_start=k+1; for (k=0, j=start; j<=finish+1; j++) place_in[k++]=!excluded[i1][j]; for (k=0, j=start; j<=finish; j++) temp_order[2][k++]=order[j]; n_in_order=k; add_orders(temp_order[1],n_start,locus[i1],place_in,&temp_order[2],1, n_in_order,&temp_order[3],&n_orders); ++n_in_order; for (k=0, j=start; j<=finish+1; j++) place_in[k++]=!excluded[i2][j]; n= 3+n_orders; add_orders(temp_order[1],n_start,locus[i2],place_in,&temp_order[3], n_orders,n_in_order,&temp_order[n],&n_orders); for (i=0; i\t",i+1); pr(); for (j=0; jlocus[j]= temp_order[n+i][j+1]; /* watch dummies on end */ map->num_loci= n_in_order; init_rec_fracs(map); converge_to_map(map); } void add_orders(frame,n_frame,new_marker,place_in,old_order,n_orders,n_loci, new_order,num_orders) int *frame; /* frame must have LEFT_DUMMY@0 and RIGHT_DUMMY@n-1 */ int n_frame; /* thus n_frame is real frame len+2 */ int new_marker; bool *place_in; /* [n]=>T/F, n is the interval right of frame[n], 0..n_f-1 */ int **old_order; /* [order#][position#] */ int n_orders; /* #old orders */ int n_loci; /* in each old order */ int **new_order; /* [order#][position#] */ int *new_orders; /* ptr to int - side-effected */ { int n, i, j, k, left, right; k=0; for (n=0; nT/F, n is the interval right of frame[n], 0..n_f-1 */ int **old_order; /* [order#][position#] */ int n_orders; /* #old orders */ int n_loci; /* in each old order */ int **new_order; /* [order#][position#] */ int *new_orders; /* ptr to int - side-effected */ { int n, i, j, k, left, right; k=0; for (n=0; n