#!/usr/bin/python

"""Fits an exponential curve to ld data generated by the pipeline. FILE NAME INPUT REQUIRED"""

___author__ = "Emma Fox (e.fox16@imperial.ac.uk)"
__version__ = "0.0.1"


##IMPORTS
import sys
from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit, fit_report
import numpy as np
import matplotlib.pyplot as plt
import pylab

##DATA
ldFile = open(sys.argv[1])

#SNP1,SNP2=np.loadtxt(ldFile, usecols=(0,1), dtype='string', unpack=True)
BPDist,r2Pear,D,DPrime,r2GLS=np.loadtxt(ldFile, usecols=(2,3,4,5,6), unpack=True)

x=BPDist
data=r2Pear	

##FUNCTIONS
#Function
def EXP(paramsEXP, x, data):
	"""Modeling the reverse exponential decay"""
	init = paramsEXP['init']
	lam = paramsEXP['lam']
	model = init * np.exp(-lam * x)
	return model - data

#Create parameter set
paramsEXP = Parameters()
paramsEXP.add('init', value=1, max=1)
paramsEXP.add('lam', value=2)

#Fit with least squares
resultEXP = minimize(EXP, paramsEXP, args=(x,), kws={'data':data})

#Model line
final = data + resultEXP.residual

#Reports
print(fit_report(resultEXP))

#Plot Results
#~ pylab.plot(x, data, 'k+')
#~ pylab.plot(x, final, 'r')
#~ pylab.show()
