#!/usr/bin/python

# preprocessing step for tlpeakfinder: for all large gaps in a given input sequence table,
#  we merge all of the scores from that region into the probe closest to its center
# this way there is not over-representation of a bunch of signals from gap regions that don't exist in the recipient strain

import table
import sys
import numpy

alignment_tab_file = sys.argv[1]
scorefile = sys.argv[2]
outscorefile = sys.argv[3]

alignment_tab=table.read_table(alignment_tab_file, sep=';', header=True, guesstypes=True)
score_tab = table.read_table(scorefile, sep="\t", header=True, guesstypes=True)

all_island_ids = numpy.array(alignment_tab.get_column("gap_id"))
all_islands = list(set(all_island_ids))
all_islands.sort()

colnames = score_tab.ordered_keys

if len(colnames) != 2:
  raise(ValueError("Input table should have two columns: probe locations and probe-level scores"))

score_locs = numpy.array(score_tab.get_column(colnames[0]))
score_vals = numpy.array(score_tab.get_column(colnames[1]))

keepflags = numpy.zeros_like(score_locs) ;# 1/true if we should keep this probe

keepflags[(all_island_ids < 2)] = 1

for island_id in all_islands:
  if (island_id < 2):
    continue

  this_island_probe_flags = (all_island_ids == island_id)
  this_island_locs = score_locs[this_island_probe_flags]
  this_island_scores = score_vals[this_island_probe_flags]
  this_island_inds = numpy.argwhere(this_island_probe_flags)

  island_new_score = numpy.median(this_island_scores)
  numlocs = len(this_island_locs)
  mid_ind = this_island_inds[numlocs/2]
  keepflags[mid_ind] = 1
  score_vals[mid_ind] = island_new_score

keepflags_bool = keepflags > 0
good_locs = score_locs[keepflags_bool]
good_scores = score_vals[keepflags_bool]
good_gaps = all_island_ids[keepflags_bool]

ostr=open(outscorefile, 'w')
ostr.write('%s\t%s\tGAPID\n' % (colnames[0], colnames[1]))

for i in range(len(good_locs)):
  ostr.write('%i\t%f\t%i\n' % (good_locs[i], good_scores[i], good_gaps[i]))

ostr.close()

