#!/usr/bin/env python
# encoding: utf-8
"""
fdr.py

Created by Nicholas Tatonetti on 2011-01-22.
Copyright (c) 2011 Stanford University. All rights reserved.
"""

import os
import csv
import sys
import numpy
import MySQLdb

if __name__ == '__main__':
    
    db = MySQLdb.connect(host="localhost", port=3307, user="root", passwd="enter_your_password",db="project_aers")
    c = db.cursor()
    
    query = """
    select stitch_id, umls_id
    from gold_drug_ae
    """
    c.execute(query)
    
    gold = dict()
    for sid, uid in c.fetchall():
        if not sid in gold:
            gold[sid] = set()
        gold[sid].add(uid)
    
    query = """
    select stitch_id1, stitch_id2
    from corr_drug_drug
    where estimate > 0.05
    """
    c.execute(query)
    
    corr = dict()
    for sid1, sid2 in c.fetchall():
        if not sid1 in corr:
            corr[sid1] = set()
        if not sid2 in corr:
            corr[sid2] = set()
        
        corr[sid1].add(sid2)
        corr[sid2].add(sid1)
    
    
    query = """
    select stitch_id, umls_id
    from pred_drug_events_%s
    where prr > 2
    """
    
    results = list()
    
    for suffix in ['null','b','e5','e']:
        c.execute(query % suffix)
        
        false = 0
        true = 0
        for sid, uid in c.fetchall():
            found = False
            for corr_sid in corr.get(sid, set()):
                if uid in gold.get(corr_sid, set()):
                    found = True
                    break
            
            if found:
                false += 1
            else:
                true += 1
        
        fdr = false/float(false + true)
        
        # meh, going back to fdr_original.py
    
    
    
    
    
    
    
    query = """
    select stitch_id, umls_id, bg_num/drug_num as ratio, -log10(pvalue) as pvalue
    from pred_drug_events_b
    join umls_list using (umls_id)
    where drug_mean > bg_mean
    and bg_num > drug_num
    """
    c.execute(query)
    drug_data = c.fetchall()
    
    bin_keys = ['0-50', '50-100', '100-200', '200-400', '400-600', '600-100000000']
    
    bins = dict()
    for key in bin_keys:
        bins[key] = set()
    
    for sid, uid, ratio, p in drug_data:
        bin = [k for k in bins.keys() if int(k.split('-')[0]) < ratio <= int(k.split('-')[1])]
        if len(bin) != 1:
            print >> sys.stderr, "Houston, we have a problem with: %f" % ratio
        else:
            bins[bin[0]].add(sid)
    
    odds = dict()
    explainable = dict()
    p_cutoff = 15
    
    for sid, uid, ratio, p in drug_data:
        
        if not sid in odds:
            odds[sid] = [0,0,0,0,None]
        
        if not sid in explainable:
            other_uid_sets = [gold.get(other_drug,set()) for other_drug in corr.get(sid,set())]
            if len(other_uid_sets) > 0:
                explainable[sid] = reduce(operator.or_,other_uid_sets)
            else:
                explainable[sid] = set()
        
        if p > p_cutoff:
            # Predicted TRUE
            if uid in explainable[sid]:
                # EXPLAINED
                odds[sid][0] += 1
            else:
                # NOT EXPLAINED
                odds[sid][2] += 1
        else:
            # Not Predicted
            if uid in explainable[sid]:
                # EXPLINED
                odds[sid][1] += 1
            else:
                # NOT EXPLAINED
                odds[sid][3] += 1
    
    for sid, (a, b, c, d, r) in odds.items():
        if not c == 0 and not d == 0 and not b == 0:
            odds[sid][4] = (a/float(c))/(b/float(d))
    
    for i,key in enumerate(bin_keys):
        odds_ratios = [odds[sid][4] for sid in bins[key] if not odds[sid][4] is None]
        num_predicted = [odds[sid][0] + odds[sid][2] for sid in bins[key] if not odds[sid][4] is None]
        num_not_predicted = [odds[sid][1] + odds[sid][3] for sid in bins[key] if not odds[sid][4] is None]
        print key, numpy.mean(odds_ratios), numpy.std(odds_ratios), numpy.mean(num_predicted)/numpy.mean(num_not_predicted)