#!/usr/bin/python

# preprocess some data for tlpeakfind by looking through a location/score file and adding strand and homology information
# we mark as nan any data point at a probe on the array that is not in the score file
# also, we fill in the homology column with the distance to the nearest significant gap
# Needed arguments: 
#  -table with locations and scores
#  -output file name
#  -configuration file to be used for tlpeakfind
#  -length of the donor strain genome

import sys
import numpy
from tlpeakfind import dmod
import ConfigParser
import bisect


locscorefile = sys.argv[1]
outfile=sys.argv[2]
strandfile = sys.argv[3]
GENOME_LENGTH = int(sys.argv[4])

locs = list(numpy.loadtxt(locscorefile, dtype='int', usecols=(0,), skiprows=1))
scores = list(numpy.loadtxt(locscorefile, dtype='float', usecols=(1,), skiprows=1))
gapids = list(numpy.loadtxt(locscorefile, dtype='int', usecols=(2,), skiprows=1))

probe_str = open(strandfile)

thisline = probe_str.readline()
thisline = probe_str.readline()

probe_locs = []
probe_strands = []


while (thisline != ""):
  probedir, probeloc = thisline.split()
  probe_locs.append(int(probeloc))
  if probedir == 'Counterclockwise':
    probe_strands.append('-')
  else:
    probe_strands.append('+')

  thisline = probe_str.readline()

if (locs != sorted(locs)):
  raise("Must have locs sorted")

if (probe_locs != sorted(probe_locs)):
  raise("Must have probe_locs sorted")


def find_index(mylist, target, startindex):
  # look through the list for target, stopping if it hits values greater than target
  target_ind = bisect.bisect_left(mylist,target,startindex)

  if mylist[target_ind] != target:
    raise(ValueError("Item not present in list"))

  return target_ind

# find the locations of all gap endpoints
# for each gap, we find the first and last probe inside of that gap
#  and use them as the endpoints

gap_endpoints = []
i=0

while i < len(locs):
  if gapids[i] > 1:
    before_gap = locs[i-1]

    while gapids[i] > 1:
      i += 1

    after_gap = locs[i]

    before_in_probes = bisect.bisect_right(probe_locs, before_gap)
    after_in_probes = bisect.bisect_left(probe_locs, after_gap)

    gap_endpoints.append(probe_locs[before_in_probes])
    gap_endpoints.append(probe_locs[after_in_probes-1])
  else:
    i += 1


# now combine all the data
ostr=open(outfile, 'w')
h_i = 0
p_i = 0
s_i = 0
print "\n"
for p_i in range(len(probe_locs)):
  if p_i%1000 == 0:
    sys.stdout.write("\rWorking on probe %i" % p_i)
    sys.stdout.flush()

  this_loc = probe_locs[p_i]
  this_dir = probe_strands[p_i]

  try:
    s_i = find_index(locs, this_loc, s_i)
    this_score = scores[s_i]
  except:
    this_score = numpy.nan
    gapdist=0
    ostr.write("%i\t%f\t%s\t%f\n" % (this_loc, this_score, this_dir, gapdist))
    continue

  if (gapids[s_i] > 1):
    gapdist = 0
  else:
    gapdists = [dmod(i,this_loc,GENOME_LENGTH) for i in gap_endpoints]
    gapdist = numpy.min(gapdists)

  ostr.write("%i\t%f\t%s\t%f\n" % (this_loc, this_score, this_dir, gapdist))

print ""
ostr.close()




