#!/usr/bin/python

# Procedure for identifying significant enrichments or depletions in cross-strain transposon library donation data
# The input file must contain four columns, in order:
#   -the genomic location of each probe/data point
#   -the score at each data point (any roughly noramlly distributed statistic will do)
#   -the strand of each data point ("+" or "-")
#   -the homology score at each data point (see the glint distribution for example calculations); should be in [0.0,1.0]
#  The input file should contain all possible data points, with 'nan' present for the score at points that should be excluded from processing

# The peaks will be identified as described at https://tavazoielab.c2b2.columbia.edu/petefred/glint/. Briefly:
#   -the scores of probes will be normalized based on their homology scores to render the distributions at each homology level roughly equivalent
#   -the data set will be normalized to a *background distribution* with a median of 0 and standard deviation of 1
#   -with outlying data masked, an approximate empricial autocorrelation will be calculated to use in significance calling
#   -the running median score over a sliding window centered on each probe will be calculated
#   -significance calls are made at the probe level by calculating p-values based on simulations from a background distribution and then applying a specified FDR
#   -peaks are called by merging all neighboring or nearby significant probes, and their extent reported

import sys
import scipy as sp
import scipy.stats as st
import scipy.optimize as so
import numpy as np
import rpy2.robjects as robjects
import copy
import ConfigParser
import pylab

############## UTILITY FUNCTIONS ####################

def dmod(n1,n2,mod):
  # Return the shortest distance between two sites on a circular genome

  smaller = min(n1%mod,n2%mod)
  bigger = max(n1%mod,n2%mod)

  direct=bigger-smaller
  wrap=smaller + mod - bigger

  return min(direct,wrap)



############# STEPS IN THE ANALYSIS METHOD ###############
def read_tlpeak_data(filename, params):
  # Read the raw score data in a specified input file and do some preprocessing
  

  # read and store each column
  locs = np.loadtxt(filename, dtype='int', usecols=(0,))
  vals = np.loadtxt(filename, usecols=(1,))
  strands = np.loadtxt(filename, usecols=(2,), dtype='string')
  homologies = np.loadtxt(filename, usecols=(3,))

  # do some basic validation
  for strand in strands:
    if (strand != '+' and strand != '-'):
      raise(ValueError("Strand can only be '+' or '-'"))

  strand_real = [int(strand == '+') for strand in strands]

  for loc in locs:
    if (loc < 1 or loc > params['GENOME_LENGTH']):
      raise(ValueError("Genomic locations should be 1-indexed values between 1 and genome_length, inclusive"))

  #for h in homologies:
  #  if (h<0.0 or h>1.0):
  #    raise(ValueError("Homologies must be contained in [0.0,1.0]"))

  # mask out regions of data that are either unspecific or extreme enough to be excluded from correlation calculations
  tail_quantile = params['EXTREME_TAIL_QUANTILE_RAW']
  nan_masks = np.invert(np.isnan(vals))
  extr_val_low = st.scoreatpercentile(vals[nan_masks], tail_quantile)
  extr_val_hi = st.scoreatpercentile(vals[nan_masks], 100-tail_quantile)
  val_masks = np.invert(np.isnan(vals) | np.greater(vals, extr_val_hi) | np.less(vals, extr_val_low))

  # calculate another mask that flags all probes too close to a significant gap
  gap_masks = (homologies > params['MINDIST_FROM_GAP'])

  return {'locs' : locs, 'scores' : vals, 'strands' : strand_real, 'homologies' : homologies, 'val_mask' : val_masks, 'nan_mask' : nan_masks, 'gap_mask' : gap_masks}

def normalize_raw_scores(inpdat, params):
  # Take a dictionary containing raw scores for tlpeak calling and return a version with the scores normalized
  # we do the normalization only on the basis of the 'background' portion of the distribution

  # first center the distribution, and then normalize it 
  val_mask = inpdat['val_mask']
  nan_mask = inpdat['nan_mask']
  new_scores = inpdat['scores'] - np.median(inpdat['scores'][val_mask])
  
  # find quantiles in the background distribution corresponding to
  #  quartiles in the overall distribution
  firstq = st.scoreatpercentile(new_scores[nan_mask], 25)
  thirdq = st.scoreatpercentile(new_scores[nan_mask], 75)
  iqr=thirdq - firstq
  new_scores = (new_scores * (1.34898 / iqr))

  retdat = inpdat.copy()
  retdat['scores'] = new_scores
  return(retdat)


def apply_homology_normalization(inpdat, params):
  # Control for the correlation between homology and magnitude of signal
  # To do this, we do a loess estimation of abs(score) vs homology
  #  throughout the genome, and divide the score of each probe by the predicted value
  #  based on its homology

  locs = inpdat['locs']
  scores = inpdat['scores']
  simvals = inpdat['homologies']
  goodflags = inpdat['nan_mask']

  goodlocs = locs[goodflags]
  goodscores = scores[goodflags]
  goodsims = simvals[goodflags]

  # now do a loess curve for everything else
  robjects.globalenv['simvals'] = robjects.FloatVector(goodsims)
  robjects.globalenv['magscores'] = robjects.FloatVector(np.abs(goodscores))
  #robjects.r('png(\"look.png\")')
  #robjects.r('plot(simvals, magscores)')
  robjects.r('h.loess <- loess(\"mags ~ sims\", data=data.frame(mags=magscores, sims=simvals), span=%f)' % params['HOMOLOGY_LOESS_ALPHA'])
  #robjects.r('lines( 1:max(simvals), predict(h.loess, data.frame(sims = 1:max(simvals))), col="red")')
  #robjects.r('dev.off()')
  predictions = robjects.r('predict(h.loess)')

  newscores = np.array(scores, copy=True)
  newscores[goodflags] = goodscores / predictions

  newdat = {'locs' : locs, 'scores' : newscores, 'homologies' : simvals, 'strands' : inpdat['strands'], 'val_mask' : inpdat['val_mask'], 'nan_mask' : inpdat['nan_mask'], 'gap_mask' : inpdat['gap_mask']}

  return(newdat)

def get_ind_windows_probewise(indat, params):
  # for each probe, give the indices of all probes within 
  # 2*params.RAVG_WINDOW_HALFWIDTH_PROBES probes of that probe
  # return the average value at each probe in indat
  # this function is designed to give an appropriate set of probes
  # to consider in calculating of moving average-based statistics later
  # note that we generate this set even for nan probes, although only non-nan
  # probes can contribute -- nan probes are removed from the list of good indices
  # in this version of the function we calculate at the probe, rather than bp or gnee,
  # level

  goodflags = indat['nan_mask']
  locs = indat['locs']

  #goodlocs = locs[goodflags]
  width = params['RAVG_WINDOW_HALFWIDTH_PROBES']

  ret_probeinds = []


  for center_i in range(len(locs)):
    start_i = center_i - 1
    end_i = center_i + 1
    if (end_i >= len(goodflags)):
      end_i = 0
    curr_window = []

    if (goodflags[center_i]):
      curr_window = [center_i]

    goodcount_r = 0
    if (goodflags[start_i]):
      goodcount_r += 1
      curr_window.insert(0,start_i)

    while goodcount_r < width:
      start_i -= 1
      if (goodflags[start_i]):
        goodcount_r += 1
        curr_window.insert(0,start_i)

    goodcount_l = 0
    if (goodflags[end_i]):
      goodcount_l += 1
      curr_window.append(end_i)

    while goodcount_l < width:
      end_i += 1
      if end_i >= len(goodflags):
        end_i = 0
      if (goodflags[end_i]):
        goodcount_l += 1
        curr_window.insert(0,end_i)

    ret_probeinds.append(copy.copy(curr_window))

  return(ret_probeinds)



def calc_moving_median(normdat, window_inds, params):
  # given the set of indices contributing to the moving average at each probe,
  # return the median centered on that point

  scores = normdat['scores']

  newscores = np.zeros(len(window_inds))

  # this can be done much more efficiently, but there's no need for the added complexity
  for i in range(len(window_inds)):
    newscores[i] = np.median( scores[window_inds[i]] )

  return(newscores)



def consolidate_peaks_byprobe(peaklocs_raw, indat, scores, params):
  # consolidate all flagged probes and return peaks passing specified cutoffs
  # any probes that are significant, have the same sign, and are
  #  separated by  less than params.SIG_COMBINE_WINDOWSIZE_PROBES probes
  # are merged into a single peak
  # then, any peaks which cover more than params.MIN_PEAKWIDTH_PROBES probes
  #  are reported (this is, as the name suggests, in base pairs, *not* probes)
  # for each real peak we return a tuple containing the locations and scores of
  #  all probes in it
  # peaks must also have all probes with the same sign -- encounter of a *significant* 
  #  probe of the opposite sign immediately stops peak extension

  locs = indat['locs']

  maxgap = params['SIG_COMBINE_WINDOWSIZE_PROBES']+1

  peak_locs = list(np.arange(len(peaklocs_raw))[peaklocs_raw > 0])

  merged_peaks = []

  if len(peak_locs) == 0:
    return []

  this_peak_locs = [ peak_locs.pop(0) ]

  # handle the special case of only one significant probe
  if len(peak_locs) == 0:
    if (1>params['MIN_PEAKWIDTH_BP']):
      merged_peaks.append((locs[this_peak_locs], [scores[this_peak_locs[0]]]))

    return merged_peaks


  # check if we can wrap the first peak into the back half of the genome
  while( (len(peak_locs) > 0) and  dmod( peak_locs[-1],this_peak_locs[0],len(locs)) < maxgap) and (np.sign(scores[peak_locs[-1]]) == np.sign(scores[this_peak_locs[0]]) ):
    this_peak_locs.insert(0,peak_locs.pop(-1))

  # now see how far this peak extends
  while ( len(peak_locs) > 0 ) and ( abs(this_peak_locs[-1] - peak_locs[0] ) < maxgap) and ( np.sign(scores[peak_locs[0]]) == np.sign(scores[this_peak_locs[-1]])):
    this_peak_locs.append( peak_locs.pop(0) )

  if len(this_peak_locs) >= params['MIN_PEAKWIDTH_PROBES']:
    this_peak_scores = scores[this_peak_locs]
    merged_peaks.append( ( locs[this_peak_locs], this_peak_scores)  )

  # now go forward through the genome and find any other peaks
  while len(peak_locs) > 0:
    this_peak_locs = [ peak_locs.pop(0) ]
    while ( (len(peak_locs) > 0) ) and (abs(this_peak_locs[-1] - peak_locs[0]) < maxgap ) and ( np.sign(scores[peak_locs[0]]) == np.sign(scores[this_peak_locs[-1]])):
      this_peak_locs.append( peak_locs.pop(0))

    #if ( dmod(locs[this_peak_locs[0]], locs[this_peak_locs[-1]], params['GENOME_LENGTH']) > params['MIN_PEAKWIDTH_BP']):
    if len(this_peak_locs) >= params['MIN_PEAKWIDTH_PROBES']:
      this_peak_scores = scores[this_peak_locs]
      merged_peaks.append( ( locs[this_peak_locs], this_peak_scores)  )


  return merged_peaks


def merge_peaks_bydist(merged_peaks_in, params):
  # do one final merge of nearby peaks, once we've done all of the necessary filtering
  # this time we merge based on genomic distances, and take any nearby peaks of the
  # same sign
  # we also discard any prospective peaks shorter than MIN_PEAKWIDTH_BP

  merged_peaks = copy.deepcopy(merged_peaks_in)

  mergedist = params['SIG_COMBINE_WINDOWSIZE']
  minwidth = params['MIN_PEAKWIDTH_BP']

  new_peaks = []

  if len(merged_peaks) == 0:
    return merged_peaks

  if len(merged_peaks) == 1:
    thispeak = merged_peaks[0]
    if dmod(thispeak[0][0],thispeak[0][1],params['GENOME_LENGTH']) >= minwidth:
      return merged_peaks
    else:
      return []

  def calc_min_peakdist(peak1,peak2,genome_length):
    # return the minimum distance between probes in two peaks
    # we assume probes within each peak are sorted, and that the end of
    #  peak 1 is closest to the start of peak 2

    mydist = dmod(peak1[0][-1],peak2[0][0],genome_length)
    return mydist

  def merge_two_peaks(peak1,peak2):
    # merge peak1 and peak2, assuming peak1 comes first

    return ( np.concatenate( (peak1[0], peak2[0])), np.concatenate( (peak1[1], peak2[1] )))


  new_peaks.append(merged_peaks.pop(0))
  while (len(merged_peaks) > 0) and (calc_min_peakdist(merged_peaks[-1],new_peaks[0],params['GENOME_LENGTH']) < mergedist) and (np.sign(np.mean(merged_peaks[-1][1])) == np.sign(np.mean(new_peaks[0][1]))) :
    add_peak = merged_peaks.pop(-1)
    new_peaks[0] = merge_two_peaks(add_peak,new_peaks[0])

  while (len(merged_peaks) > 0) and (calc_min_peakdist(new_peaks[-1],merged_peaks[0],params['GENOME_LENGTH']) < mergedist) and (np.sign(np.mean(merged_peaks[0][1])) == np.sign(np.mean(new_peaks[-1][1]))) :
    add_peak = merged_peaks.pop(0)
    new_peaks[-1] = merge_two_peaks(new_peaks[-1],add_peak)

  while (len(merged_peaks) > 0):
    new_peaks.append(merged_peaks.pop(0))
    while (len(merged_peaks) > 0) and (calc_min_peakdist(new_peaks[-1],merged_peaks[0],params['GENOME_LENGTH']) < mergedist) and (np.sign(np.mean(merged_peaks[0][1])) == np.sign(np.mean(new_peaks[-1][1]))) :
      add_peak = merged_peaks.pop(0)
      new_peaks[-1] = merge_two_peaks(new_peaks[-1],add_peak)


  # filter out peaks that aren't long enough

  # check the first peak specially
  firstpeak = new_peaks.pop(0)
  new_peaks_2 = []
  if firstpeak[0][0] > firstpeak[0][-1]:
    firstpos = firstpeak[0][0]-params['GENOME_LENGTH']
    lastpos = firstpeak[0][-1]
    if (lastpos-firstpos) >= minwidth:
      new_peaks_2.append(firstpeak)
  else:
    new_peaks.insert(0,firstpeak)
  new_peaks_2 += [peak for peak in new_peaks if ( dmod(np.min(peak[0]), np.max(peak[0]), params['GENOME_LENGTH']) >= minwidth) ] 
  return new_peaks_2

  
def print_peaklocs_formatted(peaklocs_filt, outpeakfile, params):
  # print the set of peaks returned by consolidate_peaks

  ostr=open(outpeakfile, 'w')

  ostr.write("%15s   %15s %15s\n" % ("Start location", "End location", "Mean score"))

  #print peaklocs_filt

  for locs, scores in peaklocs_filt:
    startloc = locs[0]
    endloc = locs[-1]
    scores=np.array(scores)
    meanscore = np.mean(scores[np.isfinite(scores)])
    ostr.write("%15i - %15i % 15f\n" % (startloc, endloc, meanscore))


def find_continuous_chunks(datflags,maxgap=0):
  # given a set of true/false flags for whether each given position is good,
  # return a list of lists, where each element is a set of positions with no
  # 'false' values interspersed
  # continuous sets of 'false' values less than or equal to maxgap are ignored

  this_start = 0
  this_end = 0

  if np.all(datflags):
    return [ range(len(datflags)) ]

  good_windows = []

  this_gaplen = 0

  def update_gaplen(curr_gaplen,newflag):
    # update the current gap length based on the value of the new flag
    if newflag:
      curr_gaplen = 0
    else:
      curr_gaplen += 1

    return curr_gaplen

  if (datflags[0] or datflags[-1] or datflags[1]):
    if not datflags[0]:
      this_gaplen = 1

    while this_gaplen <= maxgap:
      this_start -= 1
      if (this_start+len(datflags)) <= 0:
        break
      this_gaplen = update_gaplen(this_gaplen,datflags[this_start])

    this_start += this_gaplen

    this_gaplen = 0
    this_end=-1

    while this_gaplen <= maxgap:
      this_end += 1
      if this_end >= ( len(datflags) + this_start ):
        break
      this_gaplen = update_gaplen(this_gaplen,datflags[this_end])

    this_end -= (this_gaplen - 1)

    good_windows.append(range(this_start,this_end))

  this_start = this_end

  while (this_start < len(datflags)):

    while (this_start < len(datflags)) and (not(datflags[this_start])):
      this_start += 1

    this_end = this_start
    this_gaplen=0
    while (this_end < len(datflags)) and (this_gaplen <= maxgap):
      this_end += 1
      if this_end == len(datflags):
        break

      this_gaplen = update_gaplen(this_gaplen, datflags[this_end])

    if this_end == len(datflags):
      break

    good_windows.append(range(this_start,this_end-this_gaplen+1))
    this_start = this_end

  return good_windows

def plot_acf(vals,maxlag):
  # use R to plot and display an autocorrelation function
  robjects.globalenv['vals'] = robjects.FloatVector(vals)
  robjects.r('acf(vals, lag.max=%i)' % (maxlag,))

def calc_sample_acf(values,maxdist):
  # calculate the acf for a given subset of values
  retval = np.zeros(maxdist)
  retlen = min(len(values)-1,maxdist)
  robjects.globalenv['vals'] = robjects.FloatVector(values)
  myacf = np.ravel(np.array(robjects.r('acf(vals,lag.max=%i,plot=F)$acf' % (maxdist,))))
  retval[:retlen] = myacf[1:]
  return retval

def estimate_acf(values, usewindows, maxdist):
  # given a set of index windows in values, calculate the weighted average
  #  of the acf over those windows
  # for each lag, we weight each contribution by the number of pairs in
  #  that window contributing to that lag

  all_acfs = np.zeros((len(usewindows),maxdist))

  for i,window in enumerate(usewindows):
    all_acfs[i,:] = calc_sample_acf(values[window],maxdist)

  # now average everything together
  allweights = np.zeros(maxdist)
  full_acf = np.zeros(maxdist)

  for i in range(len(usewindows)):
    for j in range(maxdist):
      this_weight = max(0,len(usewindows[i])-j-1)
      full_acf[j] += this_weight * all_acfs[i,j]
      allweights[j] += this_weight

  return full_acf / allweights

def estimate_bg_parameters(scoredict,params):
  # given normalized scores, estimate all needed quantities for the background distribution
  # we estimate the mean/sd based on the median and
  #  mad of the actual distribution
  # Then estimate the acf based on the empirical
  #  acf of all sufficiently long windows
  # return the mean, sd, and empirical acf of the background distribution

  scores  = scoredict['scores']
  goodscores = scores[scoredict['nan_mask']]

  # first estimate the center of the background distribution
  bg_mean = np.median(goodscores)
  print 'bg_mean is %f' % bg_mean

  # now estimate the spread of the background distribution
  bg_sigma = estimate_bg_sigma_pooled(goodscores,bg_mean)
  print 'bg_sigma is %f' % bg_sigma

  # finally, estimate the covariance structure
  # we identify all continuous windows that have complete, non-outlier data (except
  #  with isolated gaps allowed) and use them to get the empirical acf
  #  following the cmarrt method (but with a weighted average taken over these windows)
  use_flags = scoredict['val_mask'] & scoredict['gap_mask']
  cont_windows = find_continuous_chunks(use_flags,maxgap=1)
  use_windows = [w for w in cont_windows if len(w) > params['ACF_WINDOWSIZE']]

  # do imputation for any nan values
  robjects.r('library(zoo)')
  robjects.globalenv['scores'] = robjects.FloatVector(scores)
  robjects.r('scores.clean = na.approx(scores)')
  scores_clean = np.array(robjects.globalenv['scores.clean'])


  # calculate the acf and smooth it
  print "Calculating empirical ACF based on %i windows" % len(use_windows)
  emp_acf = estimate_acf(scores_clean,use_windows,params['ACF_WINDOWSIZE'])
  #loess_weights = np.zeros(len(emp_acf)) ;# weight based on the number of contributing cases

  #for wind in use_windows:
  #  windlen = len(wind)
  #  for i in range(len(loess_weights)):
  #    loess_weights[i] += windlen - i

  #spline_weights = np.sqrt(loess_weights) / bg_sigma

  #robjects.globalenv['emp.acf'] = robjects.FloatVector(emp_acf)
  #robjects.globalenv['loess.weights'] = robjects.FloatVector(loess_weights)
  #robjects.r('myloess = loess(emp.acf ~ seq(length(emp.acf)), weights=loess.weights,span=0.1)')
  #acf_smooth = np.array(robjects.r('predict(myloess)'))

  #myspl = si.splrep(range(len(emp_acf)), emp_acf, spline_weights)
  #acf_smooth = si.splev(range(len(emp_acf)),myspl)

  # now return all the needed values
  return (bg_mean, bg_sigma,emp_acf)

def calc_median_bg_dist(bg_mean,bg_sigma,emp_acf,params):
  # calculate a background distribution for the running median
  # we do this by simulating from a distribution parameterized by our 
  #  earlier estimates
  # return the simulations and the resulting kde

  # recall that arima.sim does a burn-in automatically,
  #  so we don't need to do it manually
  sim_length = 2*params['RAVG_WINDOW_HALFWIDTH_PROBES'] + 1

  # first simulate a bunch of values from the null distribution
  #  for this we use an ar(n) model, with n the length of the window
  robjects.globalenv['myacf'] = robjects.FloatVector(emp_acf)
  robjects.globalenv['bg.mean'] = float(bg_mean)
  robjects.globalenv['bg.sigma'] = float(bg_sigma)
  robjects.r('myacf = c(1.0,myacf)')
  robjects.r('coeffs = acf2AR(myacf)')
  #robjects.r('print(coeffs)')
  robjects.r('coeffs.use = coeffs[%i,]' % (params['ACF_WINDOWSIZE']))
  #robjects.r('print(coeffs.use)')
  robjects.r('sim.medians = vector(mode="numeric",length=%i)' % params['N_BG_SIM'])
  print "  Simulating from null distribution..."
  for i in range(params['N_BG_SIM']):
    if i % 1000 == 0:
      sys.stdout.write("\r  Working on %i of %i" % (i,params['N_BG_SIM']))
      sys.stdout.flush()
    robjects.r('sim.vals = arima.sim( list(ar=coeffs.use),%i )' % sim_length)
    robjects.r('sim.vals = (sim.vals * bg.sigma) + bg.mean')
    robjects.r('sim.medians[%i] = median(sim.vals)' % (i+1))

  robjects.r('print(coeffs.use)')
  print ""
  print "Finished with simulations"
  mymedians = np.array(robjects.globalenv['sim.medians'])
  mykde = st.gaussian_kde(mymedians)
  return(mykde,mymedians)



def call_ma_peaks_fdr_fromsims(ma_scores, bg_scores, params):
  # given simulated draws from a background distribution, calculate p values for the ma scores, and do significance calling
  # return an array with 1 in all locations where the ma statistic is significant,
  #  and zero elsewhere
  # to do this, we calculate p values by comparison to the background distribution
  # distribution for the ma statistics, and then apply the fdr specified
  # in params.PEAK_FDR

  bg_sorted = np.sort(bg_scores)
  scores_sort_order = np.argsort(ma_scores)
  revsort_order = np.argsort(scores_sort_order)
  scores_sorted = ma_scores[scores_sort_order]

  print "calculating p values"

  ## older metho if I don't use a kde
  score_ind = 0
  score_ranks = []

  for bg_ind in range(len(bg_sorted)):
    bgval = bg_sorted[bg_ind]
  
    if score_ind >= len(scores_sorted):
      break
  
    while scores_sorted[score_ind] < bgval:
      #print "Assigning score %f less than bgval %f at bgind %i" % (scores_sorted[score_ind],bgval,bg_ind)
      score_ranks.append(bg_ind)
      score_ind += 1
      if score_ind >= len(scores_sorted):
        break
  
    bg_ind += 1
  
  #print "second pass"
  while score_ind < len(scores_sorted):
    score_ranks.append(len(bg_sorted))
    score_ind += 1

  score_ranks = np.array(score_ranks)
  pvals = 2 * np.fmin( score_ranks / float(len(bg_sorted)), 1-(score_ranks / float(len(bg_sorted))))

  pvals=np.array(pvals)
  goodflags = np.zeros_like(pvals)
  pval_sort_order = np.argsort(pvals)
  pvals_sorted = pvals[pval_sort_order]
  fdr = params['PEAK_FDR']
  N = float(len(pvals))

  best_i = -1
  
  for i in range(len(pvals)):
    if pvals_sorted[i] <= (fdr * (i+1) / N):
      best_i = i
 
  if best_i == -1:
    print "Warning: found no significant probes"
    return np.zeros_like(pvals)

  threshold_p = pvals_sorted[best_i]
  print 'p value threshold is %f' % threshold_p

  good_flags_sortedbyval = (pvals <= threshold_p)
  good_flags = good_flags_sortedbyval[revsort_order]


#  score_center = np.median(ma_scores)
#  if np.sum(good_scores < score_center) > 0:
#    lower_thresh = np.max(good_scores[good_scores < score_center])
#  else:
#    lower_thresh=min(ma_scores)-1
#  if np.sum(good_scores > score_center) > 0:
#    upper_thresh = np.max(good_scores[good_scores > score_center])
#  else:
#    upper_thresh=max(ma_scores)+1
 
#  good_scores = ma_scores[good_flags]
#
#  if (np.sum(good_scores > midval) > 0):
#    thresh_hi = min(good_scores[good_scores > midval])
#  else:
#    thresh_hi = np.inf
#
#  if (np.sum(good_scores < midval) > 0):
#    thresh_lo = max(good_scores[good_scores < midval])
#  else:
#    thresh_lo = np.inf


  #print 'value thresholds are %f/%f' % (min(good_scores[good_scores > midval]),max(good_scores[good_scores < midval]))
  #print 'background dist values at this quantile are %f/%f' % (st.scoreatpercentile(bg_sorted, threshold_p/2), st.scoreatpercentile(bg_sorted,1-(threshold_p/2)))

  return good_flags


def estimate_bg_sigma_pooled(scores,meanval):
  # estimate the sd of the background distribution
  #  we do this based on the MAD from a given mean
  #  unlike estimate_bg_sigma, here we pool allof the values instead of just taking one tail

  myscores = scores[np.isfinite(scores)]
  mad_all = np.median(np.abs(myscores - meanval))

  baseline_q = sp.stats.norm.ppf(0.75)

  my_sigma = mad_all / baseline_q

  return my_sigma
  

def estimate_bg_sigma(scores,meanval):
  # estimate the sd of the background distribution
  # we do this based on the median absolute deviation,which should be robust to extreme tails
  # we calculate the MAD for the top and bottom half of the distribution, and take the
  # smaller implied scale parameter
  # (we infer the scale for a normal distirbution that would give that MAD)

  myscores = scores[np.isfinite(scores)]

  top_vals = myscores[myscores>meanval]
  bottom_vals=myscores[myscores<meanval]

  mad_top = np.median(np.abs(top_vals - meanval))
  mad_bottom = np.median(np.abs(bottom_vals - meanval))

  mymad = min(mad_top,mad_bottom)

  baseline_q = sp.stats.norm.ppf(0.75)

  my_sigma = mymad / baseline_q

  return my_sigma
  

def init_params():
  # supply some defaults for all parameters that we need

  myparams = {}

  myparams['EXTREME_TAIL_QUANTILE_RAW'] = 2
  myparams['RAVG_WINDOW_HALFWIDTH_PROBES'] = 10
  myparams['HOMOLOGY_LOESS_ALPHA'] = 0.3
  myparams['PEAK_FDR']=0.01
  myparams['MIN_PEAKWIDTH_BP']=5000
  myparams['MIN_PEAKWIDTH_PROBES']=10
  myparams['SIG_COMBINE_WINDOWSIZE']=5000
  myparams['SIG_COMBINE_WINDOWSIZE_PROBES'] = 1
  myparams['N_BG_SIM']=500000
  myparams['MINDIST_FROM_GAP'] = 55
  myparams['ACF_WINDOWSIZE'] = 41

  return(myparams)


def read_params(conffile, myparams):
  # read parameters from a config file

  int_items = ['RAVG_WINDOW_HALFWIDTH_PROBES', 'MIN_PEAKWIDTH_BP', 'MIN_PEAKWIDTH_PROBES', 'SIG_COMBINE_WINDOWSIZE', 'SIG_COMBINE_WINDOWSIZE_PROBES', 'N_BG_SIM','MINDIST_FROM_GAP','ACF_WINDOWSIZE']
  float_items = ['EXTREME_TAIL_QUANTILE_RAW', 'HOMOLOGY_LOESS_ALPHA', 'PEAK_FDR']

  config=ConfigParser.RawConfigParser()
  config.optionxform = str
  config.read(conffile)

  for name,val in config.items('tlpeakfinder'):

    if int_items.count(name) > 0:
      myparams[name] = int(val)
    elif float_items.count(name) > 0:
      myparams[name] = float(val)
    else:
      raise(ValueError("Unrecognized tlpeakfinder parameter: %s" % name))

  return myparams

def plot_emp_acf(acf,outfile):
  # plot the empirical acf to a figure file

  pylab.figure()
  acf_vals = [1] + list(acf)
  pylab.plot(acf_vals)

  pylab.axhline(0.0, color='g',ls='--')

  pylab.ylabel('Empirical ACF')
  pylab.xlabel('# Probes')
  pylab.savefig(outfile)


########## MAIN FUNCTION (for use from the command line) #############

if (__name__ == "__main__"):
  import sys
  import argparse

  parser=argparse.ArgumentParser(description = "Call peaks from GLINT data", epilog='infile must contain data positions, values, strandedness, and homology information. It is most easily generated with combine_probeinfo_with_gapdist.py')
  parser.add_argument('-c', '--conffile', default=None)
  parser.add_argument('-o', '--outprefix', default='glint_output')
  parser.add_argument('-g', '--genomelength', default=None)
  parser.add_argument('-n', '--no_homology', action='store_true', default=False)

  parser.add_argument('infile')

  args = parser.parse_args()

  infile=args.infile
  outprefix=args.outprefix
  conffile = args.conffile
  genomelength = int(args.genomelength)

  use_homology = not (args.no_homology)

  # set up some needed parameters
  outpeakfile_raw=outprefix + "_tlpeakfinder_peaks_unmerged.txt"
  outpeakfile=outprefix + "_tlpeakfinder_peaks.txt"
  outscorefile=outprefix + "_tlpeakfinder_normed_scores.txt"
  outmafile=outprefix + "_tlpeakfinder_normed_ma_scores.txt"
  outacfplot=outprefix + "_acfplot.pdf"
  #outnhscorefile = outprefix + "_tlpeakfinder_normed_nohcorr_scores.txt"

  params = init_params()
  params['GENOME_LENGTH'] = genomelength
  if conffile is not None:
    read_params(conffile, params)
  else:
    print "Warning! Using default parameters (including for genome size)"
  
  print "Running tlpeakfind with parameters:"
  for key,val in params.iteritems():
    print "\t%s: %s" % (key, val)

  # now run the actual analysis

  print " ---- Normalizing scores ---- "
  rawdat = read_tlpeak_data(infile, params)
  normdat = normalize_raw_scores(rawdat, params)
  if use_homology:
    hndat  = apply_homology_normalization(normdat, params) 
  else:
    hndat = copy.deepcopy(normdat)

  window_inds = get_ind_windows_probewise(hndat, params)

  #  ostr=open('out_inds','w')
  #  for i in range(len(window_inds)):
  #    inds = window_inds[i]
  #    inds.sort()
  #    ostr.write("%i %i %s %s\n" % (i,hndat['locs'][i],inds,[hndat['locs'][j] for j in inds]))
  #
  #  ostr.close()

  np.savetxt(outscorefile, np.vstack((hndat['locs'], hndat['scores'])).transpose(), fmt='%i %.8f', delimiter='\t')
  print " ---- Calculating moving average statistic ---- "
  ma_scores = calc_moving_median(hndat, window_inds, params)
  print " ---- Calculating background distribution ----"
  bg_mean,bg_sigma,emp_acf = estimate_bg_parameters(hndat,params)
  plot_emp_acf(emp_acf, outacfplot)
  bg_dist,bg_sims = calc_median_bg_dist(bg_mean,bg_sigma,emp_acf,params)

  print " ---- Calling peaks ---- "
  peaklocs_raw = call_ma_peaks_fdr_fromsims(ma_scores, bg_sims, params)
  peaklocs_filt = consolidate_peaks_byprobe(peaklocs_raw, hndat, ma_scores, params)
  print_peaklocs_formatted(peaklocs_filt, outpeakfile_raw, params)
  peaklocs_merged = merge_peaks_bydist(peaklocs_filt,params)
  print_peaklocs_formatted(peaklocs_merged, outpeakfile, params)
  np.savetxt(outmafile, np.vstack((hndat['locs'], ma_scores)).transpose(), fmt='%i %.8f', delimiter='\t')
