import numpy as np
import pandas as pd
from pandas import Series,DataFrame
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import csv
from scipy.interpolate import interp1d
import locale
import japanize_matplotlib
from sklearn.linear_model import LinearRegression
import numpy.polynomial.polynomial as P
encoding = locale.getpreferredencoding()
plt.rcParams['figure.subplot.bottom'] = 0.25
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
ax2 = ax1.twinx()
ax3 = ax1.twinx()
covid = pd.read_csv('nhk_news_covid19_domestic_daily_data.csv')
df = DataFrame(covid)
start_date = "2021/9/22"
df_s=df.set_index("日付")[start_date:]
df_r = df_s.reset_index()
df0 = pd.DataFrame(df_r)
ax1.bar(df0["日付"],df0["国内の感染者数_1日ごとの発表数"], label="国内の感染者数_1日ごとの発表数")
ax1.xaxis.set_major_locator(mdates.DayLocator(interval=3))
ax1.tick_params(axis='x', colors ="black", rotation = 90)
ax1.tick_params(axis='y', colors ="navy")
#---------------------------------------
encoding = locale.getpreferredencoding()
df_2 = pd.read_csv('effective_reproduction_number.csv', encoding="ANSI")
n0 = len(df_2)
l_date = df_2['日付'].values[n0-1]
#-----------------------------------
df_3 = DataFrame(df_2)
df_s=df_3.set_index("日付")[start_date:]
df_r = df_s.reset_index()
# num = df_rの行数(rows)
num = len(df_r)
df_4 = pd.DataFrame(df_r)
ml = pd.DataFrame({'no': range(num)})
df_new = pd.concat([df_4, ml], axis=1)
df_new.to_csv('df_new.csv', encoding = encoding)
#特定の列の抽出
df_a=df_new.loc[:,["no", "実効再生産数"]]
df_a.to_csv('df_a.csv', encoding = encoding)
#--------線形再帰--------------------------
lr = LinearRegression()
x = df_a[['no']]
y = df_a[['実効再生産数']]
coef1 = lr.fit(x, y)
co = lr.fit(x, y)
w = lr.coef_
w0 = lr.intercept_
print('lr =', w, w0)
#---------指数回帰曲線------------------
fit2 = np.polyfit(df_a['no'], np.log(df_a['実効再生産数']), 2, w = df_a['実効再生産数'])
print('fit2 =', fit2)
fit22 = np.polyfit(df_a['no'], np.log(df_a['実効再生産数']), 2, w = df_a['実効再生産数']*df_a['実効再生産数'])
print('fit22 =', fit22)
#fit3= np.polyfit(df_a['no'], np.log(df_a['実効再生産数']), 1, w = df_a['実効再生産数']*df_a['実効再生産数'])
#print('fit3 =', fit3)
#-----------------------------------
interp = interp1d(df_a['no'], df_a['実効再生産数'],kind="cubic")
new_t = np.arange(0, num-1, 0.2) #range(num) num-1
#-------------------------
file = pd.DataFrame({'no': new_t, '実効再生産数':interp(new_t)})
#print(file)
ax2.plot(new_t, interp(new_t), label="実効再生産数", color='red')
ax2.tick_params(colors ="red")
ax3.tick_params(colors ="red")
ax2.plot(df_new['no'], df_new['実効再生産数'], 'o', color='orangered')
ymin = 0.5
ymax = 1.7
ax2.set_ylim(ymin, ymax)
ymin, ymax =ax2.set_ylim()
plt.ylim(ymin, ymax)
plt.xlim(0, num+25)
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1+handler2,label1+label2,borderaxespad=0)
ax1.grid(axis='x')
#---------------回帰分析------------------
print('\n')
print('回帰分析 \n')
print(df_2['日付'][n0-1])
print('回帰直線')
print('y =', w0[0],'+', w[0][0],'x1')
print('指数回帰(w=y)')
print('y =','exp(',fit2[2],'+',fit2[1],'*x1+',fit2[0],'*x1*x1)')
print('指数回帰(w=y^2)')
print('y =','exp(',fit22[2],'+',fit22[1],'*x1+',fit22[0],'*x1*x1)')
#---------------------------------
x1 = np.linspace(0,130,130);
plt.plot(x1, w[0]*x1+w0[0], color='red', lw=0.5, label='lienar')
plt.plot(x1, np.exp(fit2[2]+fit2[1]*x1+fit2[0]*x1*x1),linestyle='dashed', color='blue', \
lw=1 , label='exp(2次) + w=y, '+df_2['日付'][n0-1])
plt.plot(x1, np.exp(-4.34656175e-01 + 2.62067149e-03*x1 + 4.83403391e-05*x1*x1), \
linestyle='dashed', color='cyan', lw=1 ,label='exp(2次) + w=y, 2021/12/25')
plt.plot(x1, np.exp(fit22[2]+fit22[1]*x1+fit22[0]*x1*x1), color='black',lw=1 , \
label='exp(2次) + w=y^2, '+df_2['日付'][n0-1])
plt.plot(x1, np.exp(-4.10929621e-01 + 2.05032382e-03*x1 + 5.35490199e-05*x1*x1), \
color='gray', lw=1 , label='exp(2次) + w=y^2, 2021/12/25')
plt.plot(x1, 0*x1+1, color='purple')
plt.legend()
plt.grid()
plt.show()