# -*- coding: utf-8 -*-
"""
Created on Tue Jul 22 07:54:05 2014


@author: charleslelosq
Carnegie Institution for Science


"""
import sys
sys.path.append("/Users/charleslelosq/Documents/RamPy/lib-charles/")

import csv
import numpy as np
import scipy
import matplotlib
import matplotlib.gridspec as gridspec
from pylab import *
from StringIO import StringIO
from scipy import interpolate

# to fit spectra we use the lmfit software of Matt Newville, CARS, university of Chicago, available on the web
from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit, fit_report

from spectratools import *  #Charles' libraries and functions

from Tkinter import *
import tkMessageBox
from tkFileDialog import askopenfilename

#### We define a set of functions that will be used for fitting data
#### unfortunatly, as we use lmfit (which is convenient because it can fix or release 
#### easily the parameters) we are not able to use arrays for parameters... 
#### so it is a little bit long to write all the things, but in a way quite robust also...
#### gaussian and pseudovoigt functions are available in spectratools
#### if you need a voigt, fix the gaussian-to-lorentzian ratio to 1 in the parameter definition before
#### doing the data fit
def residual(pars, x, data=None, eps=None):
    # unpack parameters:
    #  extract .value attribute for each parameter
    a1 = pars['a1'].value
    a2 = pars['a2'].value
    
    f1 = pars['f1'].value
    f2 = pars['f2'].value    
    
    l1 = pars['l1'].value
    l2 = pars['l2'].value
    
    # Gaussian model
    
    peak1 = gaussian(x,a1,f1,l1)
    peak2 = gaussian(x,a2,f2,l2)    
    
    model = peak1 + peak2
    
    if data is None:
        return model, peak1, peak2
    if eps is None:
        return (model - data)
    return (model - data)/eps
    
##### CORE OF THE CALCULATION BELOW

#### CALLING THE DATA NAMES
tkMessageBox.showinfo(
            "Open file",
            "Please open the list of spectra")

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
with open(filename) as inputfile:
    results = list(csv.reader(inputfile)) # we read the data list

#### LOOP FOR BEING ABLE TO TREAT MULTIPLE DATA
#### WARNING: OUTPUT ARE AUTOMATICALLY GENERATED IN A DIRECTORY CALLED "DECONV"
#### (see end) THAT SHOULD BE PRESENT !!!!!!!!!!
for lg in range(len(results)):
    name = str(results[lg]).strip('[]')
    name = name[1:-1] # to remove unwanted ""
    sample = np.genfromtxt(name) # get the sample to deconvolute

    # we set here the lower and higher bonds for the interest region
    lb = 4700 ### MAY NEED TO AJUST THAT
    hb = 6000
    
    interestspectra = sample[np.where((sample[:,0] > lb)&(sample[:,0] < hb))]
    ese0 = interestspectra[:,2]/abs(interestspectra[:,1]) #take ese  as a percentage, we assume that the treatment was made correctly for error determination... if not, please put  sigma = None
    interestspectra[:,1] = interestspectra[:,1]/np.amax(interestspectra[:,1])*100 # normalise spectra to maximum, easier to handle after 
    sigma = abs(ese0*interestspectra[:,1]) #calculate good ese
    #sigma = None # you can activate that if you are not sure about the errors

    xfit = interestspectra[:,0] # region to be fitted
    data = interestspectra[:,1] # region to be fitted

    params = Parameters()
    ####################### FOR MELT:
    ####################### COMMENT IF NOT WANTED
    #               (Name,  Value,  Vary,   Min,  Max,  Expr)
    params.add_many(('a1',   1,   True, 0,      None,  None),
                    ('f1',   5200,  True, 750,    None,  None),
                    ('l1',   1,  True, 0,      None,  None),
                    ('a2',   1,  True, 0,      None,  None),
                    ('f2',   5400,  True, None,   None,  None),
                    ('l2',   1,  True, None,   None,  None))  
                         
    result = minimize(residual_melt, params, args=(xfit, data)) # fit data with leastsq model from scipy
    model = fit_report(params) # the report
    yout, peak1,peak2,= residual_melt(params,xfit) # the different peaks
    
    #### We just calculate the different areas up to 4700 cmm-1 and those of the gaussians
    # Select interest areas for calculating the areas of OH and H2Omol peaks
    intarea45 = sample[np.where((sample[:,0]> 4100) & (sample[:,0]<4700))]
    area4500 = np.trapz(intarea45[:,1],intarea45[:,0])
    esearea4500 = 1/sqrt(area4500) # We assume that RELATIVE errors on areas are globally equal to 1/sqrt(Area)
      
    # now for the gaussians
    # unpack parameters:
    #  extract .value attribute for each parameter
    a1 = pars['a1'].value
    a2 = pars['a2'].value
    
    l1 = pars['l1'].value
    l2 = pars['l2'].value
    
    AireG1 = gaussianarea(a1,l1)
    AireG2 = gaussianarea(a2,l2)
    
    ##### WE DO A NICE FIGURE THAT CAN BE IMPROVED FOR PUBLICATION
    fig = figure()
    plot(sample[:,0],sample[:,1],'k-')
    plot(xfit,yout,'r-')
    plot(xfit,peak1,'b-')
    plot(xfit,peak2,'b-')
    
    xlim(lb,hb)
    ylim(0,np.max(sample[:,1]))
    xlabel("Wavenumber, cm$^{-1}$", fontsize = 18, fontweight = "bold")
    ylabel("Absorption, a. u.", fontsize = 18, fontweight = "bold")
    
    text(4000,np.max(intarea45[:,1])+0.03*np.max(intarea45[:,1]),('Area OH: \n'+'%.1f' % area4500),color='b',fontsize = 16)
    text(4650,a1 + 0.05*a1,('Area pic 1$: \n'+ '%.1f' % AireG1),color='b',fontsize = 16)
    text(5000,a2 + 0.05*a2,('OH/H$_2$O$_{mol}$: \n'+'%.3f' % ratioOH_H2O+'\n+/-'+'%.3f' % eseratioOH_H2O),color='r',fontsize = 16)
    
    ##### output of data, fitted peaks, parameters, and the figure in pdf
    ##### all goes into the ./deconv/ folder
    name.rfind('/')
    nameout = name[name.rfind('/')+1::]
    namesample = nameout[0:nameout.find('.')]
    pathint = str('/deconv/') # the output folder
    ext1 = '_ydec.txt'
    ext2 = '_params.txt'
    ext3 = '.pdf'
    pathout1 = pathbeg+pathint+namesample+ext1
    pathout2 = pathbeg+pathint+namesample+ext2
    pathout3 = pathbeg+pathint+namesample+ext3
    matout = np.vstack((xfit,data,yout,peak1,peak2))
    matout = np.transpose(matout)
    np.savetxt(pathout1,matout) # saving the arrays of spectra
    fd = os.open( pathout2, os.O_RDWR|os.O_CREAT ) # Open a file and create it if it do not exist
    fo = os.fdopen(fd, "w+") # Now get a file object for the above file.
    fo.write(model) # write the parameters in it
    fo.close()
    savefig(pathout3) # save the figure
 
