###最小二乘法試驗###
import numpy as np
import pymysql
from scipy.optimize import leastsq
from scipy import integrate
###繪圖，看擬合效果###
import matplotlib.pyplot as plt
from sympy import *
 
 
path='E:\PythonProgram\oildata.txt'
 
lieh0 =[]   #初始第一列，油管高度
liev1 =[]   #初始第二列，油槍記錄的體積
 
h_median =[]  # 高度相鄰中位數(shù)
h_subtract =[]   #相鄰高度作差
v_subtract =[]   #相鄰體積作差
select_h_subtr =[]   #篩選的高度作差 ########
select_v_subtr =[]   #篩選的體積作差
 
VdtH=[]      #篩選的V 和 t 的 倒數(shù)。
 
def loadData(Spath,lie0,lie1):
 with open(Spath,'r') as f0:
   for i in f0:
    tmp=i.split()
    lie0.append(float(tmp[0]))
    lie1.append(float(tmp[2]))
 print ("source lie0",len(lie0))
 
 
def connectMySQL():
 db = pymysql.connect(host='10.**.**.**', user='root', passwd='***', db="zabbix", charset="utf8") # 校罐
 cur = db.cursor()
 
 try:
  # 查詢
  cur.execute("SELECT * FROM station_snapshot limit 10 ")
  for row in cur.fetchall():
   # print(row)
   id = row[0]
   snapshot_id = row[1]
   DateTime = row[13]
   attr1V = row[5]
   attr2H = row[6]
   print("id=%d ,snapshot_id=%s,DateTime=%s,attr1V =%s, attr2H=%s ",
     (id, snapshot_id, DateTime, attr1V, attr2H))
 except:
  print("Error: unable to fecth data of station_stock")
 
 try:
  cur.execute("SELECT * FROM can_stock limit 5");
  for row in cur.fetchall():
   # print(row)
   stockid = row[0]
   stationid = row[1]
   DateTime = row[4]
   Volume = row[5]
   Height = row[8]
   print("stockid=%d ,stationid=%s,DateTime=%s,Volume =%f, Height=%f ",
     (stockid, stationid, DateTime, Volume, Height))
 except:
  print("Error: unable to fecth data of can_snapshot")
 
 cur.close()
 db.close()
 
 
def formatData(h_med,h_subtr,v_subtr):
 lh0 = lieh0[:]
 del lh0[0]
 print("lh0 size(): ",len(lh0))
 
 lh1 =lieh0[:]
 del lh1[len(lh1)-1]
 
 print("lh1 size() : ",len(lh1))
 
 lv0 =liev1[:]
 del lv0[0]
 #print (liev1)
 print ("Souce liev1 size() : ",len(liev1))
 print ("lv1 size() :",len(lv0))
 """
 lv1 =liev1[:]
 del lv1[len(lv1)-1]
 print("lv1 size(): ",len(lv1))
 """
 h_med[:] =[(x+y)/2 for x,y in zip(lh0,lh1)]  ###采樣點(Xi,Yi)###
 print("h_med size（） ： ", len(h_med))
 
 h_subtr[:] = [(y-x) for x,y in zip(lh0,lh1)]
 print("h_subtr size（） ： ", len(h_subtr))
 # v_subtr[:] = [(y-x) for x,y in zip(lv0,lv1)]
 v_subtr[:] = lv0
 print("v_subtr size（） ： ", len(v_subtr))
 
 
def removeBadPoint(h_med,h_sub,v_sub):
 for val in h_sub:
  position=h_sub.index(val)
  if 0.01 > val > -0.01:
   del h_sub[position]
   del h_med[position]
   del v_sub[position]
 v_dt_h_ay = [(y/x) for x, y in zip(h_sub, v_sub)]
 return v_dt_h_ay
 
 
 
def selectRightPoint(h_med,h_subtr,v_dt_h_ay):
 for val in v_dt_h_ay:
  pos = v_dt_h_ay.index(val)
  if val > 20 :
   del v_dt_h_ay[pos]
   del h_med[pos]
   del h_subtr[pos]
 for val in v_dt_h_ay:
  ptr = v_dt_h_ay.index(val)
  if val < 14:
   del v_dt_h_ay[ptr]
   del h_med[ptr]
   del h_subtr[ptr]
 
 
def writeFile(h_mp, v_dt_h):
 s='\n'.join(str(num)[1:-1] for num in h_mp)
 v='\n'.join(str(vdt)[1:-1] for vdt in v_dt_h)
 open(r'h_2.txt','w').write(s)
 open(r'v_dt_h.txt','w').write(v)
 print("write h_median: ",len(h_mp))
 # print("V_dt also is (y-x) : ",v_dt_h,end="\n")
 print("Write V_dt_h : ",len(v_dt_h))
# file=open('data.txt','w')
# file.write(str(h_mp));
# file.close
 
 
def integralCalculate(coeff,xspace):
 vCalcute =[]
 x = Symbol('x')
 a, b, c, d = coeff[0]
 y = a * x ** 3 + b * x ** 2 + c * x + d
 i=0
 while (i< len(xspace)-1) :
  m = integrate(y, (x, xspace[i], xspace[i+1]))
  vCalcute.append(abs(m))
  i=i+1
 
 print("求導(dǎo)結(jié)果：",vCalcute)
 print("求導(dǎo)長度 len(VCalcute): ",len(vCalcute))
 return vCalcute
 
 ###需要擬合的函數(shù)func及誤差error###
 
def func(p,x):
 a,b,c,d=p
 return a*x**3+b*x**2+c*x+d #指明待擬合的函數(shù)的形狀，設(shè)定的擬合函數(shù)。
 
def error(p,x,y):
 return func(p,x)-y #x、y都是列表，故返回值也是個列表
 
def leastSquareFitting(h_mp,v_dt_hl):
 p0=[1,2,6,10]  #a,b,c 的假設(shè)初始值，隨著迭代會越來越小
 #print(error(p0,h_mp,v_dt_h,"cishu")) #目標(biāo)是讓error 不斷減小
 #s="Test the number of iteration" #試驗最小二乘法函數(shù)leastsq得調(diào)用幾次error函數(shù)才能找到使得均方誤差之和最小的a~c
 Para=leastsq(error,p0,args=(h_mp,v_dt_hl)) #把error函數(shù)中除了p以外的參數(shù)打包到args中
 a,b,c,d=Para[0]   #leastsq 返回的第一個值是a,b,c 的求解結(jié)果，leastsq返回類型相當(dāng)于c++ 中的tuple
 print(" a=",a," b=",b," c=",c," d=",d)
 plt.figure(figsize=(8,6))
 plt.scatter(h_mp,v_dt_hl,color="red",label="Sample Point",linewidth=3) #畫樣本點
 x=np.linspace(200,2200,1000)
 y=a*x**3+b*x**2+c*x+d
 
 integralCalculate(Para,h_subtract)
 plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #畫擬合曲線
 # plt.plot(h_mp, v_dt_hl,color="blue", label='Origin Line',linewidth=1) #畫連接線
 plt.legend()
 plt.show()
 
def freeParameterFitting(h_mp,v_dt_hl):
 z1 = np.polyfit(h_mp, v_dt_hl, 6) # 第一個擬合，自由度為6
  # 生成多項式對象
 p1 = np.poly1d(z1)
 print("Z1:")
 print(z1)
 print("P1:")
 print(p1)
 print("\n")
 x = np.linspace(400, 1700, 1000)
 plt.plot(h_mp, v_dt_hl, color="blue", label='Origin Line', linewidth=1) # 畫連接線
 plt.plot(x, p1(x), 'gv--', color="black", label='Poly Fitting Line(deg=6)', linewidth=1)
 plt.legend()
 plt.show()
 
def main():
 loadData(path, lieh0, liev1)
 connectMySQL() # 讀取oildata數(shù)據(jù)庫
 
 formatData(h_median, h_subtract, v_subtract)
 
 # 去除被除數(shù)為0對應(yīng)的點，并得到v 和 h 求導(dǎo) 值的列表
 VdtH[:] = removeBadPoint(h_median, h_subtract, v_subtract)
 print("h_median1：", len(h_median))
 
 print("VdtH1 : ", len(VdtH))
 
 # 賽選數(shù)據(jù)，去除異常點
 selectRightPoint(h_median, h_subtract, VdtH)
 print("h_median2：", len(h_median))
 print("h_subtract: ", len(h_subtract))
 print("VdtH2 : ", len(VdtH))
 h_mp = np.array(h_median)
 v_dt_h = np.array(VdtH)
 
 writeFile(h_mp, v_dt_h)
 # 最小二乘法作圖
 leastSquareFitting(h_mp, v_dt_h)
 # 多項式自由參數(shù)法作圖
 freeParameterFitting(h_mp, v_dt_h)
 
if __name__ == '__main__':
 main()