syntax.r

merged_fit <- alply(1:n_samp, 1, function(s) {
	h_seg_dat_new$h.capseg.d[[merge_ixs[1]]] <- unlist(h_seg_dat_new$h.capseg.d[merge_ixs])
	h_seg_dat_new$h.capseg.raw[[merge_ixs[1]]] <- unlist(h_seg_dat_new$h.capseg.raw[merge_ixs])
	h_seg_dat_new$h.capseg.annot[[merge_ixs[1]]]$pos <- unlist(lapply(h_seg_dat_new$h.capseg.annot[merge_ixs], "[[", "pos"))
	h_seg_dat_new$gh.wes.allele.d[[merge_ixs[1]]] <- do.call(cbind, h_seg_dat_new$gh.wes.allele.d[merge_ixs])
	h_seg_dat_new$gh.wes.allele.annot[[merge_ixs[1]]]$pos <- unlist(lapply(h_seg_dat_new$gh.wes.allele.annot[merge_ixs], "[[", "pos"))
	if (n_merge_hets >= 2) h_seg_dat_new$seg_phase_pr[[merge_ixs[1]]] <- 0.5
	if (n_merge_hets >= 2) h_seg_dat_new$consensus_phase[[merge_ixs[1]]] <- init_merged_consensus_phase(log_het_phase_prob_list)
	capture_em_fit_new$log.het.phase.prob[[merge_ixs[1]]] <- data_frame("x1" = rep(NA, n_merge_hets))
	capture_em_fit_new$het.phase.prob[[merge_ixs[1]]] <- capture_em_fit_new$log.het.phase.prob[[merge_ixs[1]]]
})

seg_ev <- adply(1:n_samp, 1, function(s) {
	tcr_pr_d <- calc_tcr_pr_d(tcr = h_seg_dat_new$h.capseg.d[[merge_ixs[1]]],
	                          rcov = h_seg_dat_new$h.capseg.raw[[merge_ixs[1]]],
	                          gamma = fit[[s]]$capture.em.fit$Theta$gamma)

	tau_opt <- fit_tau(tau_init, fit[[s]]$capture.em.fit$Theta$gamma,
	                     h_seg_dat_new$h.capseg.d[[merge_ixs[1]]],
	                     h_seg_dat_new$h.capseg.raw[[merge_ixs[1]]])
	tau_ll <- CalcCaptureSegTauLogLik(tau_opt, h_seg_dat_new$h.capseg.d[[merge_ixs[1]]],
	                                  h_seg_dat_new$h.capseg.raw[[merge_ixs[1]]],
	                                  fit[[s]]$capture.em.fit$Theta$gamma)

	snp_ll <- laply(merged_fit_out, function(l) sum(l[["capture.em.fit"]][["SNP_LL"]][merge_ixs[1]]))

	# something wrong starting here
	cons_ll <- laply(merged_fit_out, function(l) sum(l[["capture.em.fit"]][["cons.loglik"]][merge_ixs[1]], na.rm = T))
	f_hat <- aaply(1:n_samp, 1, function(r) return(merged_fit_out[[r]]$capture.em.fit$wes.f[merge_ixs[1], "H1.f.hat"]))
	phases <- extract_phases(merged_fit_out, merge_ixs[1]) %>% mutate(ixs = merge_ixs_s)
	f_pr_d <- aaply(1:n_samp, 1, function(s) {
	  return(calc_f_pr_d(merged_fit_out[[s]]$as.res$h.seg.dat, merge_ixs[1], merged_fit_out[[s]]$capture.em.fit$Theta))
})