#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 17 12:14:15 2021

@author: rachel
"""

"""
Created on Mon Aug  2 12:52:52 2021

@author: rachel
"""
import numpy as np
from numpy.linalg import inv

#-------------------------------- test data ----------------------------------#
#'''
b=5
c=3
d=2
p_=-0.08
q_=-0.001
r_=-0.5
x= np.arange(0,100,0.1)

y = np.add(np.add(b*np.exp(np.multiply(p_,x)),c*np.exp(np.multiply(q_,x))),d*np.exp(np.multiply(r_,x)))
#y = y+y*np.random.normal(0,0.02,len(x))
plt.plot(x,y)
#'''
#-----------------------------------------------------------------------------#
#def fitSumOfExponentials(x,y):
S = np.zeros(len(x))
SS = np.zeros(len(x))
SSS = np.zeros(len(x))
S[0] = 0
SS[0] = 0
SSS[0] = 0
for k in np.arange(1,len(x)):
    S[k] = S[k-1] + 0.5*(y[k]+y[k-1])*(x[k]-x[k-1])
    SS[k] = SS[k-1] + 0.5*(S[k]+S[k-1])*(x[k]-x[k-1])
    SSS[k] = SSS[k-1] + 0.5*(SS[k]+SS[k-1])*(x[k]-x[k-1])
    
    
M = np.zeros((7,7))
M[0,0] = np.sum(np.square(SSS))
M[0,1] = np.sum(np.multiply(SSS,SS))
M[1,0] = np.sum(np.multiply(SSS,SS))
M[0,2] = np.sum(np.multiply(SSS,S))
M[2,0] = np.sum(np.multiply(SSS,S))
M[0,3] = np.sum(np.multiply(SSS,np.power(x,3)))
M[3,0] = np.sum(np.multiply(SSS,np.power(x,3)))
M[0,4] = np.sum(np.multiply(SSS,np.power(x,2)))
M[4,0] = np.sum(np.multiply(SSS,np.power(x,2)))
M[0,5] = np.sum(np.multiply(SSS,x))
M[5,0] = np.sum(np.multiply(SSS,x))
M[0,6] = np.sum(SSS)
M[6,0] = np.sum(SSS)

M[1,1] = np.sum(np.square(SS))
M[1,2] = np.sum(np.multiply(SS,S))
M[2,1] = np.sum(np.multiply(SS,S))
M[1,3] = np.sum(np.multiply(SS,np.power(x,3)))
M[3,1] = np.sum(np.multiply(SS,np.power(x,3)))
M[1,4] = np.sum(np.multiply(SS,np.power(x,2)))
M[4,1] = np.sum(np.multiply(SS,np.power(x,2)))
M[1,5] = np.sum(np.multiply(SS,x))
M[5,1] = np.sum(np.multiply(SS,x))
M[1,6] = np.sum(SS)
M[6,1] = np.sum(SS)

M[2,2] = np.sum(np.square(S))
M[2,3] = np.sum(np.multiply(S,np.power(x,3)))
M[3,2] = np.sum(np.multiply(S,np.power(x,3)))
M[2,4] = np.sum(np.multiply(S,np.power(x,2)))
M[4,2] = np.sum(np.multiply(S,np.power(x,2)))
M[2,5] = np.sum(np.multiply(S,x))
M[5,2] = np.sum(np.multiply(S,x))
M[2,6] = np.sum(S)
M[6,2] = np.sum(S)

M[3,3] = np.sum(np.power(x,6))
M[3,4] = np.sum(np.power(x,5))
M[4,3] = np.sum(np.power(x,5))
M[3,5] = np.sum(np.power(x,4))
M[5,3] = np.sum(np.power(x,4))
M[3,6] = np.sum(np.power(x,3))
M[6,3] = np.sum(np.power(x,3))

M[4,4] = np.sum(np.power(x,4))
M[4,5] = np.sum(np.power(x,3))
M[5,4] = np.sum(np.power(x,3))
M[4,6] = np.sum(np.power(x,2))
M[6,4] = np.sum(np.power(x,2))

M[5,5] = np.sum(np.power(x,2))
M[5,6] = np.sum(x)
M[6,5] = np.sum(x)
M[6,6] = len(x)

Y = np.zeros((7))

Y[0] = np.sum(np.multiply(SSS,y))
Y[1] = np.sum(np.multiply(SS,y)) 
Y[2] = np.sum(np.multiply(S,y)) 
Y[3] = np.sum(np.multiply(np.power(x,3),y))
Y[4] = np.sum(np.multiply(np.power(x,2),y))
Y[5] = np.sum(np.multiply(x,y))
Y[6] = np.sum(y)


Minv = inv(M)
CS = np.inner(Minv,Y)

A = CS[0]
B = CS[1]
C = CS[2] 
D = CS[3]
E = CS[4]
F = CS[5]
G = CS[6]

# croots = np.roots([1,C,B,A])
# p = np.real(croots[0])
# q = np.real(croots[1])
# r = np.real(croots[2])

p,q,r = rootFinder(-2, CS = CS[0:3])

N = np.zeros((4,4))
N[0,0] = len(x)
N[0,1] = np.sum(np.exp(np.multiply(p,x)))
N[1,0] = np.sum(np.exp(np.multiply(p,x)))
N[0,2] = np.sum(np.exp(np.multiply(q,x)))
N[2,0] = np.sum(np.exp(np.multiply(q,x)))
N[0,3] = np.sum(np.exp(np.multiply(r,x)))
N[3,0] = np.sum(np.exp(np.multiply(r,x)))
N[1,1] = np.sum(np.exp(np.multiply(2*p,x)))
N[1,2] = np.sum(np.exp(np.multiply(p+q,x)))
N[2,1] = np.sum(np.exp(np.multiply(p+q,x)))
N[1,3] = np.sum(np.exp(np.multiply(p+r,x)))
N[3,1] = np.sum(np.exp(np.multiply(p+r,x)))
N[2,2] = np.sum(np.exp(np.multiply(2*q,x)))
N[2,3] = np.sum(np.exp(np.multiply(q+r,x)))
N[3,2] = np.sum(np.exp(np.multiply(q+r,x)))
N[3,3] = np.sum(np.exp(np.multiply(2*r,x)))

Yn = np.zeros((4))

Yn[0] = np.sum(y)
Yn[1] = np.sum(np.multiply(np.exp(np.multiply(p,x)),y)) 
Yn[2] = np.sum(np.multiply(np.exp(np.multiply(q,x)),y))
Yn[3] = np.sum(np.multiply(np.exp(np.multiply(r,x)),y))


Ninv = inv(N)
CSn = np.inner(Ninv,Yn)
lambdas = np.asarray([p,q,r])
print(CSn, p,q,r)

CSn = np.delete(CSn, np.where(min(CSn)))
np.argsort(CSn)
lambdas = lambdas[np.argsort(CSn)]
CSn = CSn[np.argsort(CSn)]
yth =  np.add(np.add(CSn[0]*np.exp(np.multiply(lambdas[0],x)),CSn[1]*np.exp(np.multiply(lambdas[0],x))),CSn[2]*np.exp(np.multiply(lambdas[2],x)))
plt.plot(x,y)
plt.plot(x,yth)
#    return p_result, q_result, dS[0], dS[1]

data = np.genfromtxt ('/home/rachel/Documents/non-homogenous-process/data_long.csv', delimiter=",")
x = data[:,0]
y = 1-data[:,1]
a,b,c,d = fitSumOfExponentials(x, y)

y_exp = np.add(c*np.exp(np.multiply(a,x)), d*np.exp(np.multiply(b,x)))
plt.plot(y, label = "data")
plt.plot(y_exp, label = "calculated")
plt.legend()

#------------------------------ variance -------------------------------------#
l1 = lambdas[0]
l2 = lambdas[1]
l3 = lambdas[2]
A1 = CSn[0]
A2 = CSn[1]
A3 = CSn[2]
L1=l1+l2+l3
L2=l1*l2+l1*l3+l2*l3
L3=l1*l2*l3
S1=A1*l1+A2*l2+A3*l3
S2=A1*l1*l1+A2*l2*l2+A3*l3*l3
S3=A1*l1*l1*l1+A2*l2*l2*l2+A3*l3*l3*l3

K3min = -S1
K1pmin = 1/2 * ( -L1+S2/S1 + np.sqrt(np.power((S1*L1-S2),2)-4*L3*S1)/S1 )
K2pmin = 1/2 * ( -L1+S2/S1 - np.sqrt(np.power((S1*L1-S2),2)-4*L3*S1)/S1 )
K1mmin = 1/2 * (S1-S2/S1 - (np.power(-S1,2)*L1+S1*S2+S1*L2-L3+np.power(S2,2)/S1-S3)/np.sqrt(np.power((S1*L1-S2),2)-4*L3*S1))
K2mmin = 1/2 * (S1-S2/S1 + (np.power(-S1,2)*L1+S1*S2+S1*L2-L3+np.power(S2,2)/S1-S3)/np.sqrt(np.power((S1*L1-S2),2)-4*L3*S1))
