"""
Discover adverse events significantly associated with drugs.
"""

import R
import os
import csv
import sys
import MySQLdb
import tempfile
import operator
from namedmatrix import NamedMatrix

db = MySQLdb.connect(host="127.0.0.1", port=3307, user="root", passwd="enter_your_password",db="project_aers")
c = db.cursor()

# Define the drug we are interested in.
# cid = 'CID000005090' # Vioxx
cid = sys.argv[1]

print >> sys.stderr, "%s: Building the drug data matrix and file." % cid
# Populate the drug_data
query = """
select report_id, stitch_id, drug_name, umls_id, adverse_event
from singlet_drug_reports
where stitch_id = '%s';
""" % cid
c.execute(query)

report_data = [row for row in c.fetchall()]
report_ids = sorted(set([r[0] for r in report_data]))
umls_ids = sorted(set([r[3] for r in report_data]))

# drug_data = NamedMatrix(None, report_ids, umls_ids)
drug_data = dict()
for row in report_data:
    report_id = str(row[0])
    umls_id = str(row[3])
    if not drug_data.has_key(report_id):
        drug_data[report_id] = set()
    drug_data[report_id].add(umls_id)

drug_data_file = tempfile.mkstemp(suffix=".csv")[-1]
drug_file = open(drug_data_file,'w')
writer = csv.writer(drug_file)
writer.writerow(['RowID'] + umls_ids)
for report_id in map(str, report_ids):
    row = [1 if umls_id in drug_data[report_id] else 0 for umls_id in umls_ids]
    writer.writerow([report_id] + row)

drug_file.close()

print >> sys.stderr, "%s: Building the background data matrix and file." % cid
# Get the background data.
query = """
select report_id, indication, umls_id, adverse_event
from singlet_ind_reports join
(
	select indication
	from corr_drug_ind
	where stitch_id = '%s'
	and corrected < 0.05
) a
using (indication);
""" % cid
c.execute(query)

drug_report_ids = report_ids
report_data = [row for row in c.fetchall()]
report_ids = sorted(set([r[0] for r in report_data]) - set(drug_report_ids) )

bg_data = dict()
for row in report_data:
    report_id = str(row[0])
    umls_id = str(row[2])
    
    if not bg_data.has_key(report_id):
        bg_data[report_id] = set()
    bg_data[report_id].add(umls_id)
    

bg_data_file = tempfile.mkstemp(suffix=".csv")[-1]
bg_file = open(bg_data_file,'w')
writer = csv.writer(bg_file)
writer.writerow(['RowID'] + umls_ids)
for report_id in map(str,report_ids):
    row = [1 if umls_id in bg_data[report_id] else 0 for umls_id in umls_ids]
    writer.writerow([report_id] + row)

bg_file.close()

print >> sys.stderr, "%s: Running the statistical analysis." % cid

result_file = "/Volumes/UsersHD/Users/nick/Stanford/AltmanLab/aers/part1/%s.csv" % cid
r_string = """
bg_data <- read.csv('%s')
rownames(bg_data) <- bg_data[,1]
bg_data <- bg_data[,-1]

drug_data <- read.csv('%s')
rownames(drug_data) <- drug_data[,1]
drug_data <- drug_data[,-1]

results <- matrix(rep(c(NA,NA,NA),ncol(drug_data)), ncol=3)
colnames(results) <- c("pvalue","drug_mean","bg_mean")
rownames(results) <- colnames(drug_data)

nmeans <- 1000
ssize1 <- round(max(10,0.01*nrow(bg_data)))
ssize2 <- round(max(10,0.01*nrow(drug_data)))
write(paste("%s:","nmeans=",nmeans,"sizes:",ssize1,ssize2, "lengths:",nrow(bg_data),nrow(drug_data)), stderr())

for (i in 1:nrow(results)) {
	#write(paste("Working on ",i,"of ",nrow(results)), stderr())
	vals1 <- bg_data[,i]
	vals2 <- drug_data[,i]
	
	v1means <- rep(0, nmeans)
	v1means <- unlist(lapply(v1means,function(x){mean(sample(vals1,ssize1))}))
	
	v2means <- rep(0, nmeans)
	v2means <- unlist(lapply(v2means,function(x){mean(sample(vals2,ssize2))}))
	
	results[i,1] <- t.test(v1means, v2means)$p.value
	results[i,2] <- mean(v2means)
	results[i,3] <- mean(v1means)
}

write.csv(results,file="%s",quote=F)
""" % (bg_data_file, drug_data_file, cid, result_file)

R.runRString(r_string)

print >> sys.stderr, "%s: Completed run, result saved to %s" % (cid, result_file)

os.system("rm %s" % bg_data_file)
os.system("rm %s" % drug_data_file)