APHCM_E

"第二问1"

import pandas as pd

#a)根据所附资料或您收集的资料，建立一个预测核武器数量的数学模型，并预测未来100年拥有核武器的国家;
#b)预测未来100年核武器数量的变化趋势，2123年核武器总数，以及每个国家的核武器数量。

proliferation=pd.read_excel('2022_APHCM_E_Data.xlsx',sheet_name='proliferation')
proliferation

a=[]
count=0
for i in proliferation['Year']:
    count=+1
    if count%10==0:
        a.append(1)
    elif count==1:
        a.append(1)
    else:
        a.append(0)

a.reverse()
proliferation['huizong']=a

temp=proliferation[proliferation['huizong']==1]
temp.to_excel('./Q2/proliferation.xlsx',index=None)
temp


stockpiles=pd.read_excel('2022_APHCM_E_Data.xlsx',sheet_name='stockpiles')
stockpiles

temp1=stockpiles.groupby('Year').sum()
temp1.reset_index(inplace=True,drop=False)
temp1

a=[]
count=0
for i in temp1['Year']:
    count=+1
    if count%10==0:
        a.append(1)
    elif count==1:
        a.append(1)
    else:
        a.append(0)
        
a.reverse()
temp1['huizong']=a

temp1.to_excel('./Q2/stockpiles1.xlsx',index=None)
temp1[temp1['huizong']==1].to_excel('./Q2/stockpiles10.xlsx',index=None)
temp1                

import matplotlib.pyplot as plt
import numpy as np

#用来正常显示中文标签
plt.rcParams['font.sana-serif']=['SimHei']

#用来正常显示负号
plt.rcParams['axes.unicode_minus']=False

#魔法函数：在notebook襄入图表
%matplotlib inline

#定义一个画板，尺寸为（8，5），dpi指定分辨率
plt.figure(figsize=(20,10),dpi=90)

plt.plot(temp1['Year'],temp1['Stockpile'])

plt.title("stockpiles-year diagram")
plt.xlabel("Year")
plt.ylabel('Stockpile')

plt.savefig('./Q2/stockpiles-year diagram.jpg')
plt.show()

import pandas as pd
import warnings
from sklearn.metrics import r2_score
warnings.filterwarnings('ignore')
from sklearn.metrics import mean_squared_error,mean_absolute_error
def mean_absolute_percentage_error(y_true, y_pred):
    
    return np.mean(np.abs((y_true - y_pred) / y_true))*100
def score(y_true, y_pred):
    #MARE
    print("MARE :")
    print(mean_absolute_percentage_error(y_true, y_pred))
    #RMSE
    print("RMSE :")
    print(np.sqrt(mean_squared_error(y_true, y_pred)))
    #MAE
    print("MAE :")
    print(mean_absolute_error(y_true, y_pred))
    #R2
    print("R2 :")
    print(np.abs(score(y_true, y_pred)))
    

import sys
!{sys.executable} -m  pip install --upgrade pip
!{sys.executable} -m pip install tensorflow  -i https://pypi.tuna.tsinghua.edu.cn/simple
!{sys.executable} -m pip install keras

dataset=temp1['Stockpile']
dataset

import matplotlib.pyplot as plt
import numpy as np
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from keras.models import Sequential, load_model

#将整型变为float
dataset = dataset.astype('float32')

from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler(feature_range=(0, 1))
dataset = scaler.fit_transform(dataset.values.reshape(-1, 1))

def create_dataset(dataset,look_back):
#这里的ook_back与timestep相同
    dataX, dataY = [], []
    for i in range(len(dataset)-look_back):
        a = dataset[i:(i+look_back)]
        dataX.append(a)
        dataY.append(dataset[i + look_back ])
    return numpy.array(dataX), numpy.array(dataY)

#训练墩据太少Look_back并不能过大
look_back = 1
trainX,trainY = create_dataset(dataset,look_back)
trainX = numpy.reshape(trainX, (trainX.shape[0], trainX.shape[1], 1))

#create and fit the LSTM network
model = Sequential()
model.add(LSTM(4, input_shape=( None,1)))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
model.fit(trainX, trainY, epochs=100, batch_size=1, verbose=2)
#model.save(os.path.join("DATA","Test" + ".h5"))
#make predictions



"模型验证"
#model = load_model(os.path.join("DATA","Test" + ".h5"))
trainPredict = model.predict(trainX)

#反归一化
trainPredict_ = scaler.inverse_transform(trainPredict)
trainY_= scaler.inverse_transform(trainY)

score(trainPredict_,trainY_)

plt.plot(temp1['Year'].values[1:].trainY_, label='observed data')
plt.plot(temp1['Year'].values[1:].trainPredict_, label='LSTM')
plt.title("stockpiles-year fitting diagram")
plt.xlabel("Year")
plt.ylabel('Stockpile')
plt.legend()
plt.savefig('./Q2/stockpiles-year fitting diagram.jpg')
plt.show()

x_input=trainY[-1]
predict_forword_number=101
predict_list=[]
predict_list.append(x_input)
while len(predict_list) < predict_forword_number:
    x_input = predict_list[-1].reshape((-1, 1, 1))
    yhat = model.predict(x_input,verbose=0)
    #预测新值
    predict_list.append(yhat)
    #取出
    
bb=scaler.inverse_transform(np.array([ i.reshape(-1,1)[:,0].tolist() for i in predict_list]))[1:]
bb    

pd.DataFrame(bb).to_excel('lstapre.xlsx',index=None)

import sys
#安装xgboost库
!{sys.executable} -m pip install xgboost -i https://pypi.tuna.tsinghua.edu.cn/simple
#安装lightgbm库
!{sys.executable} -m pip install lightgbm -i https://pypi.tuna.tsinghua.edu.cn/simple



temp1['rolling5']=temp1['stockpile'].rolling(5).mean()

a=[]
count=0
for i in temp1['Year']:
    count=+1
    if count%5==0:
        a.append(1)
    elif count==1:
        a.append(1)
    else:
        a.append(0)
        
a.reverse()
temp1['aa']=a
temp1



tempp=temp1[temp1['aa']==1]
tempp.reset_index(inplace=True,drop=True)
tempp.dropna(inplace=True)
tempp

tempp.to_excel('./Q2/rolling5.xlsx',index=None) #colling5

dataset=pd.concat([tempp[['Year,rolling5']],tempp['rolling5'].shift(1)],axis=1) #colling5
dataset.dropna(inplace=True,axis=0)
dataset.columns=['Year','X','y']
dataset[y]=dataset['y'].astype('int')
dataset

dataset.to_excel('./Q2/dataset_rolling5.xlsx',index=None)

import xgboost as xgb
import lightgbm as lgb
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
import warnings
warnings.filterwarnings('ignore')

train_X=dataset.iloc[:,1].valus.resreshape(-1, 1)
train_y=dataset.iloc[:,2].valus.resreshape(-1, 1)
text_X=dataset.iloc[::,1].valus.resreshape(-1, 1)
text_y=dataset.iloc[::,2].valus.resreshape(-1, 1)

#高维数据模型业界一般选择线性回归，速度快，准确率高
#线性回归
model_lr = lr.LinearRegression()
model_lr.fit(train_X,train_y)
#print('线性回归')

print(score(model_lr.predict(text_X.reshape(-1, 1)),text_y))

plt.plot(tempp['Year'].values[1:].text_y, label='observed data')
plt.plot(tempp['Year'].values[1:].model_lr.predict(text_X.reshape(-1, 1)), label='LinearRegression')
plt.title("stockpiles-year fitting diagram")
plt.xlabel("Year")
plt.ylabel('Stockpile')
plt.legend()
plt.savefig('./Q2/LinearRegression stockpiles-year-rolling5 fitting diagram.jpg')
plt.show()

aaaa=[]
for i in range(20):
    if i ==0:
        tt=model_lr.predict(text_X[-1].reshape(-1, 1))
    else:
        tt=model_lr.predict(np.array([aaaa[-1]]).reshape(-1, 1))
    aaaa.append(tt[0])
aaaa

#lightgbm回归
model_lgb = lgb.LGBMRegressor()
model_lgb.fit(train_X,train_y)
print('lgbm回归')

print(score(model_lgb.predict(text_X.reshape(-1, 1)),text_y))

print('\n...........')
plt.plot(temp1['Year'].values[:-1].text_y, label='observed data')
plt.plot(temp1['Year'].values[:-1].model_lgb.predict(text_X.reshape(-1, 1)), label='lgbm')
plt.title("stockpiles-year fitting diagram")
plt.xlabel("Year")
plt.ylabel('Stockpile')
plt.legend()
plt.savefig('./Q2/lgbm stockpiles-year fitting diagram.jpg')
plt.show()

aaaa=[]
for i in range(10):
    if i ==0:
        tt=model_lgb.predict(text_X[-1].reshape(-1, 1))
    else:
        tt=model_lgb.predict(np.array([aaaa[-1]]).reshape(-1, 1))
    aaaa.append(tt[0])
aaaa        

# xgboost回归
model_xgb = xgb.XGBRegressor()

model_xgb.fit(train_X,train_y)
print('xgboost回归')
print(score(model_xgb.predict(text_X.reshape(-1, 1)),text_y))

print('\n...........')
plt.plot(temp1['Year'].values[:-1].text_y, label='observed data')
plt.plot(temp1['Year'].values[:-1].model_xgb.predict(text_X.reshape(-1, 1)), label='XGBRegressor')
plt.title("stockpiles-year fitting diagram")
plt.xlabel("Year")
plt.ylabel('Stockpile')
plt.legend()
plt.savefig('./Q2/XGBRegressor stockpiles-year fitting diagram.jpg')
plt.show()

aaaa=[]
for i in range(10):
    if i ==0:
        tt=model_xgb.predict(text_X[-1].reshape(-1, 1))
    else:
        tt=model_xgb.predict(np.array([aaaa[-1]]).reshape(-1, 1))
    aaaa.append(tt[0])
aaaa

因为爆炸产生的冲击波从中心点向外传播，爆炸的能量越大，
在相同时间内冲击波传播得越远、蘑菇云的半径就越大.

测量从爆炸开始的不同时刻 t所对应的蘑菇云半径 r(t)，
如下表所示:

t/ms	r(t)/m	t/ms	r(t)/m	t/ms	r(t)/m	t/ms	r(t)/m	t/ms	r(t)/m
0.10	11.1	0.80	34.2	1.5	44.4	3.53	61.1	15.0	106.5
0.24	19.9	0.94	36.3	1.65	46.0	3.80	62.9	25.0	130.0
0.38	25.4	1.08	38.9	1.79	46.9	4.07	64.3	34.0	145.0
0.52	28.8	1.22	41.0	1.93	48.7	4.34	65.6	53.0	175.0
0.66	31.9	1.36	42.8	3.26	59.0	4.61	67.3	62.0	185.0

%检验程序
data=[
 0.10 11.1 0.80 34.2 1.5 44.4 3.53 61.1 15.0 106.5
 0.24 19.9 0.94 36.3 1.65 46.0 3.80 62.9 25.0 130.0
 0.38 25.4 1.08 38.9 1.79 46.9 4.07 64.3 34.0 145.0
 0.52 28.8 1.22 41.0 1.93 48.7 4.34 65.6 53.0 175.0
 0.66 31.9 1.36 42.8 3.26 59.0 4.61 67.3 62.0 185.0
 ];
t=[];r=[];
for i=1:5 
    t=[t;data(:,2*i-1)];
    r=[r;data(:,2*i)]; 
end
r1=log(r);t1=log(t); 
s=polyfit(t1,r1,1);
b=s(1);a=s(2);
tt=0:0.1:70;
rr=exp(a)*tt.^b; 
plot(tt,rr,t,r,'+')

%使用非线性模型线性化技术和 polyfit函数进行数据拟合实验。最终确定爆炸能量e
%x=t.^(2/5);y=r;
%polyfit(x,y,1)

%用曲线拟合  
%z=polyfit(x,y,2) 
%f=polyval(z,x); 
%plot(t,r)
%plot(x,f)

%检验程序
x=log10(t*1e-3); 
y=5/2*log10(r)-x; 
plot(x,y,'+') 
xlabel('log10(t)'); 
ylabel('5/2*log10(r)-log10(t)'); 
c=mean(y) 
hold on; 
plot(x,c,'.-'); 
hold off; 
rou0=1.25; 
e=rou0*10^(2*c) 
kiloton=e/4.184e12