#!/usr/bin/env python3
'''Rank correlation processing for a csv tabulation of counts by segment 
   First column is for whole crawl, then 100 columns for segs 0-99
   Each row is counts for some property, e.g. mime-detected or tld

   For example, assuming all.tsv has the whole-crawl warc-only counts
   and s...tsv have the segment counts, all with counts in column 1,

   tr -d ',' <all.tsv |head -100 | while read n m; do printf "%s%s\n" $n $(for i in {0..99}; do printf ",%s" $({ grep -w "w    $m\$" s${i}.tsv || echo NaN ;} | cut -f 1 ) ; done ) ; done > all_100.csv

   will produce such a file with
     * 100 rows, one for each of the top 100 counts
     * 101 columns, 0 for all and 1--100 for segs 0--99

   Usage: python3 -i spearman.py name
     where name.csv has the input
'''
     
import numpy as np
from numpy import loadtxt
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
import pylab

import sys

def qqa():
  # q-q plot for the whole crawl
  sm.qqplot(all, line='s')
  plt.gca().set_title('Rank correlation per segment wrt whole crawl (warc results only)')
  plt.show()

def qqs():
  # q-q plots for the best and worst (by variance) segments
  global xv, xworst, xbest
  xv=[d.variance for d in xd]
  xworst=xv.index(max(xv))
  xbest=xv.index(min(xv))
  print(xbest,xworst)
  sm.qqplot(x[xbest], line='s')
  plt.gca().set_title('Best segment (least variance): %s'%xbest)
  plt.show()
  sm.qqplot(x[xworst], line='s')
  plt.gca().set_title('Worst segment (most variance): %s'%xworst)
  plt.show()

def plot_x(block=True):
  plt.plot([xd[i].mean for i in range(100)],'bx',label='Mean of rank correlation of each segment x all other segments')
  plt.plot([0,99],[xm,xm],'b',label='Mean of segment x segment means')
  plt.plot(all,'rx',label='Rank correlation of segment x whole crawl')
  plt.plot([0,99],[all_m,all_m],'r',label='Mean of segment x whole crawl')
  plt.axis([0,99,0.8,1.0])
  plt.legend(loc='best')
  plt.grid(True)
  plt.show(block=block)

def hist():
  sdd=[(i,xm-(i*xsd)) for i in range(-2,3)]
  fig,hax=plt.subplots() # Thanks to https://stackoverflow.com/a/7769497
  sdax=hax.twiny()
  hax.hist([xd[i].mean for i in range(100)],color='lightblue')
  hax.set_title('Mean of rank correlation of each segment x all other segments')
  for s,v in sdd:
       sdax.plot([v,v],[0,18],'b')
  sdax.set_xlim(hax.get_xlim())
  sdax.set_ylim(hax.get_ylim())
  sdax.set_xticks([v for s,v in sdd])
  sdax.set_xticklabels([str(s) for s,v in sdd])
  plt.show()

def first_diff(ranks):
  # first disagreement with baseline == {1,2,...}
  for i in range(len(ranks)):
    if ranks[i]!=i+1.0:
      return i
  return i+1

def ranks():
  # Combine segment measures:
  #  segID,rank corr. wrt all,inverse variance, mean cross rank corr.,first disagreement
  # convert to ranks, smallest value == highest rank
  all_ranked=stats.rankdata(-all,method='average') # invert since
                                                   #  large corr is good
  x_variance_ranked=stats.rankdata([xd[i].variance for i in range(100)])
                                                  # small corr variance is good
  x_mean_ranked=stats.rankdata([-(xd[i].mean) for i in range(100)])
                                                   # invert since
                                                   #  large mean corr is good
  fd_ranked=stats.rankdata([-first_diff(x_ranks[i]) for i in range(100)])
                                                   # invert since
                                                   #  large first diff is good
  return np.array([[i,
                    all_ranked[i],
                    x_variance_ranked[i],
                    x_mean_ranked[i],
                    fd_ranked[i]] for i in range(100)])

counts=loadtxt(sys.argv[1]+".csv",delimiter=',')
# "If axis=0 (default), then each column represents a variable, with
#        observations in the rows"
# So each column is a sequence of counts, for whole crawl in column 0
#   and for segments 0--99 in columns 1--100
corr=stats.spearmanr(counts,nan_policy='omit').correlation

all=corr[0][1:]
all_s=stats.describe(all)
all_m=all_s.mean

x=np.array([np.concatenate((corr[i][1:i],
                            corr[i][i+1:])) for i in range(1,101)])
# The above, although transposed, works because the correlation matrix
#  is symmetric
xd=[stats.describe(x[i]) for i in range(100)]
xs=stats.describe(np.array([xd[i].mean for i in range(100)]))
xm=xs.mean
xsd=np.sqrt(xs.variance)

x_ranks=[stats.rankdata(-counts[:,i],method='average') for i in range(1,101)]

### I need to review rows, e.g. counts[0] is an array of 101 counts
###   for the most common label in the complete crawl,
###   from the complete crawl and all the segments
### versus columns, e.g. counts[:,0] is an array of 100 decreasing counts
###   for all the labels in the complete crawl