[关闭]
@semprelibera 2017-06-05T09:20:48.000000Z 字数 30121 阅读 785

附录A 相关代码

Thesis

一、 抓取数据

在PEER钢筋混凝土柱试验数据库中抓取试验信息

1.1 抓取并存储矩形截面柱子的数据

文件名为：scra.py

  1. import re
  2. import requests
  3. import ssl
  4. from pandas import DataFrame
  5. # 由于目标网页用了https协议，需要证书验证防止第三方入侵，这里直接关闭验证，否则会读取失败
  6. ssl._create_default_https_context = ssl._create_unverified_context
  7. from requests.packages.urllib3.exceptions import InsecureRequestWarning
  8. requests.packages.urllib3.disable_warnings(InsecureRequestWarning)
  11. ############################
  12. # Step1, 定义函数getlinkfile#
  13. # 获取链接页的全部超级链接    #
  14. ###########################
  15. def getlinkfile(url):
  16.     # 这是用来保存结果的list
  17.     res = []
  18.     # 读取网页源代码，存储在r中
  19.     r = requests.get(url,verify = False)
  20.     # 检验读取是否成功
  21.     try:
  22.         r.raise_for_status()
  23.     except Exception as exc:
  24.         print('There was a problem: %s' % (exc))
  25.     # 把网页的内容以纯文本的形式存储到变量data中，data的类型为string
  26.     data = r.text
  27.     # 定义正则表达式，用来匹配网页上的超级链接
  28.     hrefRegEx = re.compile(r"(?<=href=\").+?(?=\")|(?<=href=\').+?(?=\')")
  29.     # 匹配所有的超级链接，结果存储在link_list中
  30.     link_list = hrefRegEx.findall(data)
  31.     #mylinkfile = open('linkfile.txt','w')
  32.     # 筛选掉我们不需要的超级链接
  33.     for url in link_list:
  34.         if 'xml&id' not in url:
  35.             res.append('https://nisee.berkeley.edu' + url)
  36.     # 返回结果
  37.     return res
  40. ######################################
  41. # Step2，定义函数getthedata，用来匹配每 #
  42. # 一个矩形截面柱信息页上我们需要信息      #
  43. #####################################
  44. def getthedata(url):
  45.     # 读取网页源代码
  46.     r = requests.get(url,verify = False)
  47.     # 检验读取是否成功
  48.     try:
  49.         r.raise_for_status()
  50.     except Exception as exc:
  51.         print('There was a problem: %s' % (exc))
  52.     # 获取网页源代码，以纯文本形式保存在data中
  53.     data = r.text
  54.     # 用正则表达式分别匹配我们需要的矩形截面柱子的构件信息
  55.     NameReExpr = re.compile("(Name:).*?\">(.*?)<\/td>")
  56.     TypeReExpr = re.compile("(Type:).*?\">(.*?)<\/td>")
  57.     f_cReExpr = re.compile("(Concrete Strength:).*?\">(.*?)\(")
  58.     Tr_steelgradReExpr = re.compile("Transverse Steel:.*?\">(Grade):(.*?)<br>")
  59.     Tr_yStreReExpr = re.compile("Transverse Steel:.*?>(Yield Stress:)(.*?)\(MPa\)")
  60.     Tr_uStreReExpr = re.compile("Transverse Steel:.*?(Strength:)(.*?)\(MPa\)")
  61.     Lg_steelgradReExpr = re.compile("Longitudinal Steel.*?\">(Grade):(.*?)<br>")
  62.     Lg_yStreReExpr = re.compile("Longitudinal Steel.*?(Yield Stress:)(.*?)\(MPa\)")
  63.     Lg_uStreReExpr = re.compile("Longitudinal Steel.*?(Strength:)(.*?)\(MPa\)")
  64.     Geo_WidthReExpr = re.compile("(Width:)(.*?)\(mm\)")
  65.     Geo_DepthReExpr = re.compile("(Depth:)(.*?)\(mm\)")
  66.     Geo_LmeasureExpre = re.compile("(L-Measured:)(.*?)\(mm\)")
  67.     Axi_loadReExpre = re.compile("(Axial Load:).*>(.*?)\(kN\)")
  68.     Lg_DiaReExpre = re.compile("(Diameter:)<.*?\">(.*?)\(mm\)")
  69.     Lg_NumbExpre = re.compile("(Number of Bars:)<.*?\">(.*?)<\/td>")
  70.     Lg_ClearoverExpre = re.compile("(Clear cover:)(.*?)\(mm\)")
  71.     Tr_ShearLagsExpre = re.compile("(Number of Shear Legs:).*?\">(\d)")
  72.     Tr_DiaReExpre = re.compile("(Bar Diameter:)(.*?)\(mm\)")
  73.     Tr_BarSpaceExpre = re.compile("(Hoop Spacing:)(.*?)\(mm\)")
  76.     # 进行匹配，相应的结果以tuple的形式返回
  77.     name = NameReExpr.search(data)
  78.     typeof = TypeReExpr.search(data)
  79.     f_c = f_cReExpr.search(data)
  80.     Tr_steelgrad = Tr_steelgradReExpr.search(data)
  81.     Tr_yStre = Tr_yStreReExpr.search(data)
  82.     Tr_uStre = Tr_uStreReExpr.search(data)
  83.     Lg_steelgrad = Lg_steelgradReExpr.search(data)
  84.     Lg_yStre = Lg_yStreReExpr.search(data)
  85.     Lg_uStre = Lg_uStreReExpr.search(data)
  86.     Width = Geo_WidthReExpr.search(data)
  87.     Depth = Geo_DepthReExpr.search(data)
  88.     Length = Geo_LmeasureExpre.search(data)
  89.     Axi_load = Axi_loadReExpre.search(data)
  90.     Lg_Dia = Lg_DiaReExpre.search(data)
  91.     Lg_Numb = Lg_NumbExpre.search(data)
  92.     Lg_Clear =Lg_ClearoverExpre.search(data)
  93.     Tr_ShearLags = Tr_ShearLagsExpre.search(data)
  94.     Tr_Dia = Tr_DiaReExpre.search(data)
  95.     Tr_BarSpace = Tr_BarSpaceExpre.search(data)
  96.     # 返回结果
  97.     return (name,typeof,f_c,Tr_steelgrad,Tr_yStre,Tr_uStre,Lg_steelgrad,
  98.             Lg_yStre,Lg_uStre,Width,Depth,Length,Axi_load,Lg_Dia,Lg_Numb,Lg_Clear,
  99.             Tr_ShearLags,Tr_Dia,Tr_BarSpace)
  103. #####################################
  104. #Step3、运用之前两个函数进行数据的抓取工作#
  105. #####################################
  107. # 定义字典，mydata，{feature1:[item1,item2,...],feature2:[item1,item2,...],...}，
  108. # 用来存储二维数据表，它可以直接转化为DataFrame
  109. mydata = {}
  110. # 用来计数的变量，由于前254条数据是矩形截面的，这里的代码写的不是很规范，矩形和圆形
  111. # 截面的两种情况，其实可以用条件控制语句实现。
  112. count = 1
  113. # 调用获取超级链接的函数，获取超链接页面中我们需要的全部超级链接，结果保存为list
  114. sourceurl = 'https://nisee.berkeley.edu/spd/servlet/search'
  115. myurllist = getlinkfile(sourceurl)
  116. # 这个list存放读取过的矩形截面链接
  117. myrectlist = [] 
  118. # 这些是mydata字典里面的keys
  119. mykeys = ('Test_Name','Type_of_Column','Concrete_Strength',
  120.               'Tr_steelgrad','Tr_ySre','Tr_uStre','Lg_steelgrad',
  121.               'Lg_yStre','Lg_uStre','Width','Depth','Length','Axi_load',
  122.               'Lg_Dia','Lg_Num','Lg_Clearover','Tr_ShearLags','Tr_Dia',
  123.               'Tr_BarSpace')
  124. # 初始化mydata字典
  125. for mykey in mykeys:
  126.     mydata[mykey] = []
  127. # 用for循环读取所有的矩形截面构件信息
  128. for url in myurllist:
  129.     print('-----------------------------------------------')
  130.     print('Try getting ' + str(count) + ' of all the data')
  132.     # 调用getthedata函数获取信息
  133.     tableattrib = getthedata(url)
  135.     print('Successfully logging in the site!')
  136.     # 这里的if语句作用不大
  137.     if tableattrib[1].group(2) == 'Rectangular':
  138.         # 记录矩形截面柱的链接
  139.         myrectlist.append(url)
  140.         # 读取所有信息
  141.         for (mykey,el) in zip(mykeys,tableattrib):
  142.             mydata[mykey].append(el.group(2))
  145.         count += 1
  146.     # 读到254个的时候就进入圆形截面柱了，跳出循环（用条件控制可以解决此问题）
  147.     if count == 254:
  148.         break
  151. # 把结果存储为DataFrame格式
  152. mydataFrame = DataFrame(mydata)
  153. # 写入csv文件
  154. mydataFrame.to_csv("peerdatabase_rect.csv")
  155. # 用集合差运算获得圆形截面柱的超级链接
  156. spirurllist = list(set(myurllist)-set(myrectlist))
  157. # 把超级链接分为两部分存储，1、矩形截面柱。2、圆形截面柱
  158. rectfile = open("linkfile_rect.txt",'w')
  159. spirfile = open("linkfile_spir.txt",'w')
  160. # 存储矩形截面柱链接
  161. for url1 in myrectlist:
  162.     rectfile.write('https://nisee.berkeley.edu' + url1 + '\n')
  163. # 存储圆形截面柱链接
  164. for url2 in spirurllist:
  165.     spirfile.write('https://nisee.berkeley.edu' + url2 + '\n')
  166. rectfile.close()
  167. spirfile.close()

1.2 抓取圆形截面柱子数据

文件名为：scra_spir.py

  1. import re
  2. import ssl
  3. import requests
  4. #import pandas as pd
  5. from pandas import DataFrame
  6. ssl._create_default_https_context = ssl._create_unverified_context
  7. from requests.packages.urllib3.exceptions import InsecureRequestWarning
  8. requests.packages.urllib3.disable_warnings(InsecureRequestWarning)
  10. ###################
  11. #Step1, 定义相关函数#
  12. ###################
  13. # 定义函数getthedataofspir，用来抓取圆形截面柱的信息
  14. def getthedataofspir(url):
  15.     # 流程与抓取矩形截面柱的过程类似
  16.     r = requests.get(url,verify = False)
  17.     # 检验是否读取成功
  18.     try:
  19.         r.raise_for_status()
  20.     except Exception as exc:
  21.         print('There was a problem: %s' % (exc))
  22.     data = r.text
  23.     # 用正则表达式匹配信息
  24.     NameReExpr = re.compile("(Name:).*?\">(.*?)<\/td>")
  25.     TypeReExpr = re.compile("(Type:).*?\">(.*?)<\/td>")
  26.     f_cReExpr = re.compile("(Concrete Strength:).*?\">(.*?)\(")
  27.     Tr_yStreReExpr = re.compile("Transverse Steel:.*?>(Yield Stress:)(.*?)\(MPa\)")
  28.     Lg_yStreReExpr = re.compile("Longitudinal Steel.*?.*?\">(Yield Stress:)(.*?)\(")
  29.     Lg_uStreReExpr = re.compile("Longitudinal Steel.*?(Strength:)(.*?)\(MPa\)")
  30.     Geo_DiameterExpr = re.compile("(Diameter:).*?\">(.*?)\(mm\)")
  31.     Geo_LmeasureExpre = re.compile("(L-Measured:)(.*?)\(mm\)")
  32.     Axi_loadReExpre = re.compile("(Axial Load:).*>(.*?)\(kN\)")
  33.     Lg_DiaReExpre = re.compile("(Diameter:)<.*?\">(.*?)\(mm\)")
  34.     Lg_NumbExpre = re.compile("(Number of Bars:)<.*?\">(.*?)<\/td>")
  35.     Lg_ClearoverExpre = re.compile("(Cover to Center of Hoop Bar).*?\">(.*?)\(mm\)")
  36.     Tr_DiaReExpre = re.compile("(Diameter Spiral).*?\">(.*?)\(mm\)")
  37.     Tr_BarSpaceExpre = re.compile("(Hoop Spacing, Sv).*?\">(.*?)\(mm\)")
  40.     # 搜索需要的信息
  41.     name = NameReExpr.search(data)
  42.     typeof = TypeReExpr.search(data)
  43.     f_c = f_cReExpr.search(data)
  44.     Tr_yStre = Tr_yStreReExpr.search(data)
  45.     Lg_yStre = Lg_yStreReExpr.search(data)
  46.     Lg_uStre = Lg_uStreReExpr.search(data)
  47.     Colum_Dia = Geo_DiameterExpr.search(data)
  48.     Length = Geo_LmeasureExpre.search(data)
  49.     Axi_load = Axi_loadReExpre.search(data)
  50.     Lg_Dia = Lg_DiaReExpre.search(data)
  51.     Lg_Numb = Lg_NumbExpre.search(data)
  52.     Lg_Clear =Lg_ClearoverExpre.search(data)
  53.     Tr_Dia = Tr_DiaReExpre.search(data)
  54.     Tr_BarSpace = Tr_BarSpaceExpre.search(data)
  55.     # 返回信息
  56.     return (name,typeof,f_c,Tr_yStre,
  57.             Lg_yStre,Lg_uStre,Colum_Dia,Length,Axi_load,Lg_Dia,Lg_Numb,Lg_Clear,
  58.             Tr_Dia, Tr_BarSpace)
  60. #################
  61. #Step2, 数据抓取 #
  62. #################
  64. # 定义字典mydata，{feature1:[item1,item2,...],feature2:[item1,item2,...],...},
  65. # 用来存放二维数据表，可以直接转化为DataFrame
  66. mydata = {}
  67. count = 1
  68. # 定义mydata里的keys
  69. mykeys = ('Test_Name','Type_of_Column','Concrete_Strength'
  70.           ,'Tr_ySre',
  71.               'Lg_yStre','Lg_uStre','Colum-Dia','Length','Axi_load',
  72.               'Lg_Dia','Lg_Num','Lg_Clearover','Tr_Dia',
  73.               'Tr_BarSpace')  
  74. # 初始化mydata字典
  75. for mykey in mykeys:
  76.     mydata[mykey] = []
  77. # 读取存放圆形截面柱超级链接的文件
  78. mysprifile = open("linkfile_spir.txt",'r')
  79. while True:
  80.     try:
  81.         print('-----------------------------------------------')
  82.         print('Try getting ' + str(count) + ' of all the data')
  83.         url = mysprifile.readline().strip()
  85.         tableattrib = getthedataofspir(url)
  86.         print('Successfully logging in the site!')
  87.         for (mykey,el) in zip(mykeys,tableattrib):
  88.             mydata[mykey].append(el.group(2))
  89.         count += 1
  90.     except:
  91.         break
  94. mysprifile.close()
  95. # 转化为DataFrame
  96. mydataFrame = DataFrame(mydata)
  97. # 把数据存储到本地硬盘
  98. mydataFrame.to_csv("peerdatabase_spir_new.csv")

1.3 后期处理

  • Pan/Li的数据用，在线OCR软件读取即可
  • 对于PEER数据的文件，手工在csv文件中加入一列，作为承载力
  • 最后把三个文件改名为"01Result_Rect.csv"，"02Result_Spir.csv"，和"03dataofLiBing.csv"

二、 数据清洗

主要目的是把三个文件合并为一个数据表，anaofData
文件名为：datadigging.py

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Fri Apr 14 16:03:11 2017
  6. @author: Bruce
  7. """
  8. import pandas as pd
  9. import Concrete_Strength_ana
  10. import numpy as np
  11. import pylab
  13. # 定义函数cleanofNumber,把带逗号','的字符串数字，转化为浮点数
  14. def cleanofNumber(Acolumn):
  15.     temp = []
  16.     for el in Acolumn:
  17.         temp.append(np.float(el.replace(',','')))
  18.     return temp
  20. # Reading data from raw csv files
  21. # 在原始csv文件中读取数据
  22. # Original data consists of three DataFrames with different data structure
  23. # 原始数据由三个不同数据结构的DataFrame构成
  24. data_rect = pd.read_csv("01Result_Rect.csv",index_col='ind')
  25. data_spir = pd.read_csv("02Result_Spir.csv",index_col= 'ind')
  26. data_libing = pd.read_csv("03dataofLiBing.csv",index_col= 'ind')
  28. #############################
  29. # Step1，处理Pan/Li文献中的数据#
  30. #############################
  32. # 丢掉不需要的列
  33. ana_libing = data_libing.drop(['Reference','Configuration','ms'],axis=1)
  34. # 把以百分号为单位的配筋率转化为无量纲的配筋率
  35. ana_libing['Ratio_Lg'] = ana_libing['Ratio_Lg'] / 100
  36. ana_libing['Ratio_Tr'] = ana_libing['Ratio_Tr'] / 100
  38. #######################
  39. # Step2，处理矩形截面数据#
  40. #######################
  42. # 计算配筋率
  43. data_rect['Ratio_Lg'] = (0.25 * np.pi * data_rect['Lg_Dia']**2) * data_rect['Lg_Num'] / (
  44.         data_rect['Depth'] * data_rect['Width'])
  45. data_rect['Ratio_Tr'] = (0.25 * np.pi * data_rect['Tr_Dia']**2) * data_rect['Tr_ShearLags'] / (
  46.         data_rect['Tr_BarSpace'] * data_rect['Width'])
  47. # 丢掉不需要的列
  48. ana_rect = data_rect.drop(['Lg_Clearover','Lg_Dia','Lg_Num','Lg_steelgrad','Lg_uStre',
  49.                            'Tr_BarSpace','Tr_Dia','Tr_ShearLags',
  50.                            'Tr_steelgrad','Tr_uStre','Type_of_Column'],axis = 1)
  52. #########################
  53. # Step3，把数据中带有字符串#
  54. # 的列全部转化为纯数字     #
  55. ########################
  57. ana_rect['Axi_load'] = cleanofNumber(ana_rect['Axi_load'])
  58. ana_rect['Length'] = cleanofNumber(ana_rect['Length'])
  59. ana_rect['Tr_ySre']=cleanofNumber(ana_rect['Tr_ySre'])
  60. data_spir['Axi_load'] = cleanofNumber(data_spir['Axi_load'])
  61. data_spir['Length'] = cleanofNumber(data_spir['Length'])
  62. data_spir['Colum-Dia'] = cleanofNumber(data_spir['Colum-Dia'])
  63. data_spir['Tr_ySre']=cleanofNumber(data_spir['Tr_ySre'])
  65. #########################
  66. # Step3，处理圆形截面的数据#
  67. ########################
  69. # 计算配筋率
  70. data_spir['Ratio_Lg'] = (data_spir['Lg_Dia']**2) * data_spir['Lg_Num'] / (
  71.         data_spir['Colum-Dia']**2)
  72. data_spir['Ratio_Tr'] = (data_spir['Tr_Dia']**2) * np.pi * 0.5/ (
  73.         (data_spir['Colum-Dia'] - data_spir['Lg_Clearover'] * 2) * data_spir['Tr_BarSpace'])
  74. # 计算轴压比
  75. data_spir['Ratio_AxisLoad'] = data_spir['Axi_load'] * 10**3 / (
  76.         data_spir['Concrete_Strength'] * data_spir['Colum-Dia']**2 * np.pi * 0.25 )
  77. data_spir['Ratio_Aspect'] = data_spir['Length'] / data_spir['Colum-Dia']
  78. # 计算界面剂
  79. data_spir['Area'] = data_spir['Colum-Dia']**2 * np.pi *0.25
  80. # 丢掉不需要的列
  81. ana_spir = data_spir.drop(['Lg_Clearover','Lg_Dia','Lg_Num','Lg_uStre',
  82.                            'Tr_BarSpace','Tr_Dia'
  83.                           ,'Type_of_Column'],axis = 1)
  85. #######################
  86. # Step4，三个数据表的合并#
  87. #######################
  89. # 先合并矩形截面数据和Pan/Li的数据，因为这些都是矩形截面数据
  90. anaOfData = pd.concat([ana_rect,ana_libing])  # just include rectangular section typle
  91. # 计算轴压比
  92. anaOfData['Ratio_AxisLoad'] = anaOfData['Axi_load'] * 10**3 / (
  93.         anaOfData['Concrete_Strength'] * anaOfData['Depth'] * anaOfData['Width'])
  94. # 计算高深比
  95. anaOfData['Ratio_Aspect'] = anaOfData['Length'] / anaOfData['Depth']
  96. # 计算截面积
  97. anaOfData['Area'] = anaOfData['Depth'] * anaOfData['Width']
  98. # 把圆形截面数据也合并进来
  99. anaOfData = pd.concat([anaOfData,ana_spir])
  100. # 筛选出高深比小于4的短柱
  101. anaOfData = anaOfData[anaOfData['Ratio_Aspect'] <= 4 ]   # give me all the short columns
  103. #####################################
  104. # Step5，计算几列数据处理时候需要用的性质#
  105. ####################################
  107. # 平均等效剪应力的计算
  108. anaOfData['v_max'] = anaOfData['V_max'] * 10**3 / anaOfData['Area']
  109. # 归一化剪切承载力
  110. anaOfData['v_max/f_c'] = anaOfData['v_max'] / anaOfData['Concrete_Strength']
  111. # 丢掉无效数据
  112. anaOfData = anaOfData[anaOfData['v_max/f_c'] != 0 ]
  113. # 考虑钢筋屈服强度，对钢筋“配筋率”进行重新计算
  114. anaOfData['Ratio_Tr'] = anaOfData['Ratio_Tr'] * anaOfData['Tr_ySre']/anaOfData['Concrete_Strength']
  115. anaOfData['Ratio_Lg'] = anaOfData['Ratio_Lg'] * anaOfData['Lg_yStre']/anaOfData['Concrete_Strength']
  116. # 丢掉不需要的列
  117. anaOfData = anaOfData.drop(['Colum-Dia','Depth','Width'],axis=1)
  118. # 丢掉有na数据的列
  119. anaOfData = anaOfData.dropna()
  120. # 筛选出轴向受压力的数据，得到最终的数据表anaofData
  121. anaOfData = anaOfData[anaOfData['Ratio_AxisLoad']>=0]
  123. #################################
  124. # Step6，debug代码。总览一下数据分布#
  125. #################################
  126. if __name__ == "__main__":
  128.     # 不加分组，直接绘制散点图
  130.     pylab.figure()
  131.     pylab.plot(anaOfData['Concrete_Strength'],anaOfData['v_max'],'m.',label='experiment data')
  132.     pylab.title('Effect of Concrete Strength')
  133.     pylab.xlabel("$f_c^{\'}$/Mpa")
  134.     pylab.ylabel('$v_{max}$/Mpa')
  135.     pylab.legend()
  136.     pylab.show()
  140.     z1,p1,x1 = Concrete_Strength_ana.linearfit(1/anaOfData['Ratio_Aspect'],anaOfData['v_max/f_c'],0.25,1)
  141.     r2 = Concrete_Strength_ana.getR2(1/anaOfData['Ratio_Aspect'],anaOfData['v_max/f_c'],p1)
  142.     r2=round(r2,2)
  143.     pylab.figure()
  144.     pylab.plot(1/anaOfData['Ratio_Aspect'],anaOfData['v_max/f_c'],'m.',label='experiment data')
  145.     pylab.plot(x1,p1(x1),'k-',label='linear fit with $R^2$= ' +str(r2))
  146.     pylab.title('Effect of Aspect Ratio')
  147.     pylab.xlabel(r"$1/\phi_{asp}$")
  148.     pylab.ylabel('$v$')
  149.     pylab.legend()
  150.     pylab.show()
  153.     anaOfaxis = anaOfData[anaOfData['Ratio_AxisLoad']>=0]
  154.     pylab.figure()
  155.     pylab.plot(anaOfaxis['Ratio_AxisLoad'],anaOfaxis['v_max/f_c'],'b.',label='experiment data')
  156.     pylab.title('Effect of ratio of axis load')
  157.     pylab.xlabel('$\psi_{axis}$')
  158.     pylab.ylabel('$v$')
  159.     pylab.legend()
  160.     pylab.show()
  164.     z2,p2,x2 = Concrete_Strength_ana.linearfit(anaOfData['Ratio_Lg'],anaOfData['v_max/f_c'],0,0.8)
  165.     r2 = Concrete_Strength_ana.getR2(anaOfData['Ratio_Lg'],anaOfData['v_max/f_c'],p2)
  166.     r2 = round(r2,2)
  167.     pylab.figure()
  168.     pylab.plot(anaOfData['Ratio_Lg'],anaOfData['v_max/f_c'],'b.',label = 'experiment data')
  169.     pylab.plot(x2,p2(x2),'k-',label='linear fit with $R^2$= ' + str(r2))
  170.     pylab.title(r'Effect of $\rho_l$')
  171.     pylab.xlabel(r'$\rho_l$')
  172.     pylab.ylabel(r'$v$')
  173.     pylab.legend()
  174.     pylab.show()
  175.     baddata = anaOfData[(anaOfData['Ratio_Lg']>0.7) & (anaOfData['v_max/f_c']<0.075)]
  179.     z3,p3,x3 = Concrete_Strength_ana.linearfit(anaOfData['Ratio_Tr'],anaOfData['v_max/f_c'],0,0.20)
  180.     r2 = Concrete_Strength_ana.getR2(anaOfData['Ratio_Lg'],anaOfData['v_max/f_c'],p3)
  181.     r2 = round(r2,2)
  182.     pylab.figure()
  183.     pylab.plot(anaOfData['Ratio_Tr'],anaOfData['v_max/f_c'],'b.',label = 'experiment data')
  184. #    pylab.plot(x3,p3(x3),'k-',label='linear fit with $R^2$= ' + str(r2))
  185.     pylab.title(r'Effect of $\rho_v$')
  186.     pylab.xlabel(r'$\rho_v$')
  187.     pylab.ylabel(r'$v$')
  188.     pylab.legend()
  189.     pylab.show()
  190.     baddata = anaOfData[anaOfData['v_max/f_c']>0.175]

三、 数据处理

3.1 混凝土强度的影响

文件名为：Concrete_Strength_ana.py

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Sat Apr 22 20:07:44 2017
  5. Mission: Find the influence of concrete strength
  6. @author: Bruce
  7. """
  8. import datadigging
  9. import pylab
  10. import pandas as pd
  11. import numpy as np
  12. ########################
  13. #Step1, 定义需要用到的函数#
  14. ########################
  16. # 根据数据表和列名，绘制数据分布柱状图
  17. def drawHist(data,column):
  18.     """
  19.     To draw a histogram of a column in data.
  20.     data is a DataFrame, and column is an attribute in it.
  21.     """
  22.     pylab.figure()
  23.     pylab.hist(data[column])
  24.     pylab.title('This is the distribution of ' + column)
  25.     pylab.xlabel(column)
  26.     pylab.ylabel("number")
  27.     pylab.show()
  29. # 根据数据表和两个列名绘制二维分布图
  30. def drawHist2D(data,column1,column2):
  31.     """
  32.     To draw a 2D histogram of two columns in data.
  33.     data is a DataFrame, and column1, column2 are two attributes in it.
  34.     """
  36.     pylab.figure()
  37.     pylab.plot(data[column1],data[column2],'b.')
  38.     pylab.title('This is the distribution of ' + column1 + ' and ' + column2)
  39.     pylab.xlabel(column1)
  40.     pylab.ylabel(column2)
  41.     pylab.legend()
  42.     pylab.show()
  44. # 计算模型的回归判定系数
  45. def getR2(x,y,ypre):
  46.     """
  47.     x is an array-like list
  48.     y is an array-like list
  49.     ypre is the prediction model
  50.     --------
  51.     Return
  52.     R^2 of the prediciton model
  53.     --------
  54.     """
  55.     x=np.array(x)
  56.     y=np.array(y)
  57.     preValue = [ypre(el) for el in x]
  58.     preValue = np.array(preValue)
  59.     res = 1 - np.sum((y-preValue)**2) / np.sum((y-np.average(y))**2)
  60.     return res
  62. # 对数据做线性回归
  63. def linearfit(x,y,minimum,maximum):
  64.     """
  65.     Return polynomial, a fucntion, and xrange
  66.     """
  67.     z1 = np.polyfit(x,y,1)
  68.     p1 = np.poly1d(z1)
  69.     xp1 = np.linspace(minimum,maximum,100)
  70.     return z1,p1,xp1
  73. if __name__ == "__main__":
  74.     ############
  75.     # 绘图的过程#
  76.     ###########
  77.     # 此过程可以预想定义函数之后在进行绘图，
  78.     # 可以减少代码的重复量
  80.     # 先根据高深比分类
  81.     anaOfAxisLoad_3 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect'] == 3]
  82.     anaOfAxisLoad_2 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect'] == 2]
  83.     anaOfAxisLoad_4 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect'] == 4]
  84.     anaOfAxisLoad_1 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect']<1.2]
  85.     anaOfAxisLoad_1_5 = datadigging.anaOfData[(datadigging.anaOfData['Ratio_Aspect']<1.7) &
  86.                                               (datadigging.anaOfData['Ratio_Aspect']>1.2)]
  87.     anaOfAxisLoad_2_5 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect']==2.5]
  91.     # 绘制各个分组的关系图
  92.     print('Concrete Strength(aspectratio=3 and axisload = 0.08-0.13)')
  93.     #anaOfAxis_more3 = anaOfAxisLoad_3[
  94.             #        (anaOfAxisLoad_3['Ratio_AxisLoad']>0.08) * (anaOfAxisLoad_3['Ratio_AxisLoad']<0.13)]
  96.     pylab.figure()
  97.     pylab.plot(anaOfAxisLoad_3['Concrete_Strength'],anaOfAxisLoad_3['v_max'],'k^')
  98.     pylab.title('Concrete Strength with aspectratio 3')
  99.     pylab.xlabel('Concrete Strength/MPa')
  100.     pylab.ylabel('$v_{max}/Mpa$')
  101.     pylab.show()
  102.     print('-----------------------------------------------------------------------')
  104.     print('Concrete Strength(aspectratio=2 and axisload = 0.08-0.13)')
  105.     #anaOfAxis_more2 = anaOfAxisLoad_2[
  106.             #        (anaOfAxisLoad_2['Ratio_AxisLoad']>0.08) * (anaOfAxisLoad_3['Ratio_AxisLoad']<0.13)]
  108.     pylab.figure()
  109.     pylab.plot(anaOfAxisLoad_2['Concrete_Strength'],anaOfAxisLoad_2['v_max'],'g^')
  110.     pylab.title('Concrete Strength with aspectratio 2')
  111.     pylab.xlabel('Concrete Strength/MPa')
  112.     pylab.ylabel('$v_{max}/Mpa$')
  113.     pylab.show()
  114.     print('-----------------------------------------------------------------------')
  117.     pylab.figure()
  118.     pylab.plot(anaOfAxisLoad_1['Concrete_Strength'],anaOfAxisLoad_1['v_max'],'g^')
  119.     pylab.title('Concrete Strength with aspectratio 1')
  120.     pylab.xlabel('Concrete Strength/MPa')
  121.     pylab.ylabel('$v_{max}/Mpa$')
  122.     pylab.show()
  123.     print('-----------------------------------------------------------------------')
  125.     anaOfAxis_more2_1 = anaOfAxisLoad_2[
  126.             (anaOfAxisLoad_2['Ratio_AxisLoad']>0) & (anaOfAxisLoad_2['Ratio_AxisLoad']<0.2)]
  128.     z1,p1,xp1 = linearfit(anaOfAxis_more2_1['Concrete_Strength'],anaOfAxis_more2_1['v_max'],0,90)
  129.     pylab.figure()
  130.     pylab.plot(anaOfAxis_more2_1['Concrete_Strength'],anaOfAxis_more2_1['v_max'],'g^'
  131.                ,xp1,p1(xp1),'k-')
  132.     pylab.title('Concrete Strength with aspectratio 2 and axisload = 0-0.2')
  133.     pylab.xlabel('Concrete Strength/MPa')
  134.     pylab.ylabel('$v_{max}/Mpa$')
  135.     pylab.show()
  136.     print('The slope k= ' + str(z1[0]))
  138.     drawHist2D(anaOfAxis_more2_1,'Ratio_Tr','Ratio_Lg')
  142.     print('-----------------------------------------------------------------------')
  145.     anaOfAxis_more2_2 = anaOfAxisLoad_2[
  146.             (anaOfAxisLoad_2['Ratio_AxisLoad']>0.2) & (anaOfAxisLoad_2['Ratio_AxisLoad']<0.5)]
  148.     z2 = np.polyfit(anaOfAxis_more2_2['Concrete_Strength'],anaOfAxis_more2_2['v_max'],1)
  149.     p2 = np.poly1d(z2)
  150.     xp2 = np.linspace(0,125,125)
  153.     pylab.figure()
  154.     pylab.plot(anaOfAxis_more2_2['Concrete_Strength'],anaOfAxis_more2_2['v_max'],'m.',
  155.                label = 'experiment data')
  156.     r2 = getR2(anaOfAxis_more2_2['Concrete_Strength'],anaOfAxis_more2_2['v_max'],p2)
  157.     r2 = round(r2,2)
  158.     pylab.plot(xp2,p2(xp2),'k-',label = 'linear fit with $R^2$ = ' + str(r2))
  159.     pylab.title('Effect of $f_c^{\'}$ with $\phi_{asp}$ = 2 and $\psi_{axis}$ = 0.2-0.5')
  160.     pylab.xlabel('$f_c^{\'}$/MPa')
  161.     pylab.ylabel('$v_{max}$/MPa')
  162.     pylab.legend()
  163.     pylab.show()
  164.     print('The slope k= ' + str(z2[0]))
  165.     print('-----------------------------------------------------------------------')
  168.     pylab.figure()
  169.     pylab.plot(anaOfAxis_more2_2['Concrete_Strength'],anaOfAxis_more2_2['v_max/f_c'],'m.',
  170.                label = 'experiment data')
  172.     pylab.title('Effect of $f_c^{\'}$ with $\phi_{asp}$ = 2 and $\psi_{axis}$ = 0.2-0.5')
  173.     pylab.xlabel('$f_c^{\'}$/MPa')
  174.     pylab.ylabel('$v_{max}$/MPa')
  175.     pylab.ylim(0,0.200)
  176.     pylab.legend()
  177.     pylab.show()
  178.     print('The slope k= ' + str(z2[0]))
  179.     print('-----------------------------------------------------------------------')

3.2 高深比的影响

文件名为:Ratio_Aspect_ana.py

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Sun Apr 23 22:53:12 2017
  5. @author: Bruce
  6. """
  8. import datadigging
  9. import Concrete_Strength_ana
  10. import pylab
  11. import pandas as pd
  12. import numpy as np
  13. # 从datadigging文件导入数据
  14. anaofdata = datadigging.anaOfData.copy()
  15. # 定义字典，用于绘图的坐标显示
  16. dicOfKeys = {'Ratio_Aspect': r'$\phi_{asp}$','Ratio_AxisLoad':r'$\psi_{axis}$',
  17.              'Ratio_Lg':r'$\rho_l$','Ratio_Tr':r'$\rho_v$',
  18.              'Concrete_Strength':'$f_c^{\'}$','v_max/f_c':'$v$'}
  20. # 定义函数，绘制实验数据的分布图
  21. def drawHist(data,column,binst=None):
  22.     """
  23.     To draw a histogram of a column in data.
  24.     data is a DataFrame, and column is an attribute in it.
  25.     """
  26.     pylab.figure()
  27.     if binst == None:
  28.         pylab.hist(data[column],color='k')
  29.     else:
  30.         pylab.hist(data[column],bins=binst,color='k')
  32.     pylab.title('The distribution of ' + dicOfKeys[column])
  33.     pylab.xlabel(dicOfKeys[column])
  34.     pylab.ylabel("# of sample")
  35.     pylab.show()
  36. # 定义函数，绘制二维的数据分布图，或者二维散点图
  37. def drawHist2D(data,column1,column2,isHist = True):
  38.     """
  39.     To draw a 2D histogram of two columns in data.
  40.     data is a DataFrame, and column1, column2 are two attributes in it.
  41.     isHist controls whether to draw a histogram or a scatter diagram
  42.     """
  43.     data = data[data[column1]!=np.nan]
  44.     data = data[data[column2]!=np.nan]
  45.     pylab.figure()
  46.     if isHist:
  47.         pylab.hist2d(data[column1],data[column2],bins = 5)
  48.         pylab.colorbar()
  49.     else:
  50.         pylab.plot(data[column1],data[column2],'ko')
  52.     pylab.title('The distribution of ' + column1 + ' and ' + column2)
  53.     pylab.xlabel(column1)
  54.     pylab.ylabel(column2)
  55.     pylab.show()
  57. # 定义函数，用于控制变量，同时绘制图像
  58. def drawRelationship(data,column1,column2,control1,control2,val1_1,val1_2,val2_1,val2_2,
  59.                      style,mytitle,xlabel,ylabel,legend1,xmin,xmax,myvarified,varified = False):
  60.     """
  61.     data is a DataFrame
  62.     To filter the data with the condtion val1_1 < control1 < val1_2 AND
  63.     val2_1 < control2 < val2_2
  64.     The x axis is column1 and y axis is column2
  65.     style is used to set the style of the diagram
  66.     ----------
  67.     Return
  68.     The a DataFrame after filtering, dataofDraw.
  69.     ----------
  70.     """
  71.     dataOfDraw = data[(data[control1]>val1_1) & (data[control1]<val1_2) & 
  72.                       (data[control2]>val2_1) & (data[control2]<val2_2)]
  73.     if column1 == 'Ratio_Aspect':
  74.         z1,p1,x1 = Concrete_Strength_ana.linearfit(1/dataOfDraw[column1],dataOfDraw[column2],xmin,xmax)
  75.         r2 = Concrete_Strength_ana.getR2(1/dataOfDraw[column1],dataOfDraw[column2],p1)
  76.         r2=round(r2,2)
  77.         pylab.figure()
  78.         pylab.plot(1/dataOfDraw[column1],dataOfDraw[column2],style,label = legend1)
  79.         pylab.plot(x1,p1(x1),'k-',label="linear fit with $R^2$= " + str(r2))
  80.     else:
  81.         z1,p1,x1 = Concrete_Strength_ana.linearfit(dataOfDraw[column1],dataOfDraw[column2],xmin,xmax)
  82.         r2 = Concrete_Strength_ana.getR2(dataOfDraw[column1],dataOfDraw[column2],p1)
  83.         r2 = round(r2,2)
  84.         pylab.figure()
  85.         pylab.plot(dataOfDraw[column1],dataOfDraw[column2],style,label = legend1)
  86.         pylab.plot(x1,p1(x1),'k-',label="linear fit with $R^2$= " + str(r2))
  88.     if varified:
  89.         pylab.plot(1/myvarified[column1],myvarified[column2],'rx',label='varified data')
  90.     pylab.title(mytitle)
  91.     pylab.xlabel(xlabel)
  92.     pylab.ylabel(ylabel)
  93.     pylab.xlim(xmin-0.005,xmax+0.005)
  94.     pylab.ylim(0,0.200)
  95.     pylab.legend()
  96.     pylab.show()
  97.     print("k = " + str(z1[0]))
  98.     print("b = " + str(z1[1]))
  99.     return dataOfDraw
  102. if __name__ == "__main__":
  103.     ##################
  104.     # 控制变量，进行绘图#
  105.     ##################
  107.     anaOfAxisLoad_2 = datadigging.anaOfData[datadigging.anaOfData['Ratio_Aspect'] == 2]
  108.     drawHist2D(anaofdata,'Ratio_AxisLoad','Ratio_Aspect',False)
  109.     drawHist2D(anaofdata,'Ratio_Lg','Ratio_Tr')
  110.     drawHist(anaofdata,'Ratio_Lg')
  111.     anaofAspect1 = anaofdata[(anaofdata['Ratio_AxisLoad']<0.9) & (anaofdata['Ratio_AxisLoad']>0.5)]
  114.     myres1 = drawRelationship(anaofdata,'Ratio_Aspect','v_max/f_c',
  115.                               'Ratio_AxisLoad','Ratio_Lg',
  116.                              0.1,0.2,0.30,0.40,'b.',
  117.                              r'Effect of $\phi_{asp}$ with $\psi_{axis}$=0.1-0.2 and $\rho_l$=0.3-0.4',
  118.                              r'1/$\phi_{asp}$','$v$',
  119.                              'experiment data')
  120.     myres2 = drawRelationship(anaofdata,'Ratio_Aspect','v_max/f_c',
  121.                               'Ratio_AxisLoad','Ratio_Lg',
  122.                              0.3,1.0,0.30,0.40,'b.',
  123.                              r'Effect of $\phi_{asp}$ with $\psi_{axis}$=0.3-1.0 and $\rho_l$=0.30-0.40',
  124.                              r'1/$\phi_{asp}$','$v$',
  125.                              'experiment data')
  128.     myres3 = drawRelationship(anaofAspect1,'Ratio_Aspect','v_max/f_c',
  129.                               'Ratio_Tr','Ratio_Lg',
  130.                              0.00,1.0,0.0,1.5,'b.',
  131.                              r'Effect of $\phi_{asp}$ with $\psi_{axis}$=0.5-0.8 and $\rho_v$=0.10-0.20',
  132.                              r'1/$\phi_{asp}$','$v$',
  133.                              'experiment data',xmin=0.22,xmax=1.00)
  135.     drawHist(myres3,'Ratio_AxisLoad')
  138.     myvarified = drawRelationship(myres3,'Ratio_Aspect','v_max/f_c','Ratio_Tr',
  139.                                   'Ratio_AxisLoad',
  140.                              0.06,0.08,0.5,1.0,'b.',
  141.                              r'Effect of $\phi_{asp}$ with $\rho_v$=0.06-0.08 and $\rho_l$=0.3-0.4 and $\psi_{axis}=0-0.2$',
  142.                              r'1/$\phi_{asp}$','$v$',
  143.                              'experiment data',xmin=0.25,xmax=1.00)
  144.     myres4 = drawRelationship(myres3,'Ratio_Aspect','v_max/f_c','Ratio_Tr',
  145.                               'Ratio_AxisLoad',
  146.                              0.06,0.08,0.0,0.3,'b.',
  147.                              r'Effect of $\phi_{asp}$ with $\rho_v$=0.06-0.08 and $\rho_l$=0.3-0.4 and $\psi_{axis}=0-0.2$',
  148.                              r'1/$\phi_{asp}$','$v$',
  149.                              'experiment data',xmin=0.25,xmax=1.00,varified=True)

3.3 轴压比、配筋率的影响

运用先前在Ratio_Aspect_ana.py文件中定义的几个函数drawHist,drawHist2D,drawRelationship，就可以进行数据分组、绘图工作。只需要简单改变输入参数即可。

3.4 机器学习算法

文件名为：machine_learning.py

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. Created on Tue Apr 25 10:22:01 2017
  5. @author: Bruce
  6. """
  8. """
  9. 1.Some basic ideals
  10. a, Hold-out Methon. Divide my data set(D) into two seperate subsets, e.g. 200(S)+100(T)
  11. Conventionly, we use 2/3 - 4/5 to train the LEARNING MACHINE
  12. b, how to apply stratified sampling with respect to my data set? It is continous
  13. distributed.
  14. maybe simply the random.sample(population, k) method to do it(without replacement)?
  15. maybe it is to arbitrary considering the random distribution of originial data set.
  18. c, We need to repeat the random sampling process for many times, and then take
  19. the average value of the output as well as its deviation.
  20. d, consider signal-to-noise ratio(SNR)--feature extraction?
  24. """
  25. #from sklearn import linear_model
  26. from sklearn import neural_network
  27. from sklearn.preprocessing import StandardScaler  
  28. from sklearn.externals import joblib
  29. import datadigging
  30. import numpy as np
  31. import pandas as pd
  32. import pylab
  33. # 从anaOfData中丢去不需要的列，得到机器学习算法需要的数据表data_machine
  34. dicOfKeys = {'Ratio_Aspect': r'$\phi_{asp}$','Ratio_AxisLoad':r'$\psi_{axis}$',
  35.              'Ratio_Lg':r'$\rho_l$','Ratio_Tr':r'$\rho_v$'}
  36. data_machine = datadigging.anaOfData.drop(['Area','Axi_load',
  37.                                            'Concrete_Strength','Length',
  38.                                            'Lg_yStre','Test_Name','Tr_ySre','V_max'
  39.                                            ,'v_max'],axis=1)
  40. data_machine = data_machine.dropna()
  44. # 定义函数，用来绘制一张图，试验数据和模型预测的对比。
  45. def plotOne(totalRes,Model):
  46.     """
  47.     totalRes is the total result of the experiments
  48.     Model is the mathematic model to be evalutated
  49.     """
  50.     myX,myy = getXy(totalRes,'v_max/f_c')
  51.     for col in ["Ratio_Aspect","Ratio_AxisLoad",
  52.                 "Ratio_Tr","Ratio_Lg"]:
  53.         # plot the actual and predicted test group
  54.         pylab.figure()
  55. #        myX = scaler.inverse_transform(myX)
  56.         pylab.plot(myX[col],Model.predict(myX)/myy,'k.',
  57.                    label="experiment data")
  58. #        pylab.plot(myX[col],Model.predict(myX),'b^',
  59. #                   label="model prediction")
  60.         pylab.title("comparison between model prediction and experiment data ")
  61.         pylab.xlabel(dicOfKeys[col])
  62.         pylab.ylabel("v")
  63.         pylab.legend()
  64.         pylab.show()
  66.         # plot the hist for difference between actual and predicted value
  67.     print("This is the difference distribution ")
  68.     pylab.figure()
  69.     pylab.hist(100 * (Model.predict(myX)-myy) / myy,15,color='k')
  70.     pylab.title("$(v_{predicted} - v_{experiment})/v_{expermiment}$   distribution")
  71.     pylab.xlabel("differences/%")
  72.     pylab.ylabel("number")
  73.     pylab.show()
  75. # 另一个绘图函数，临时定义。
  76. def plotOneMy(totalRes,Model,col,lengendlist):
  78.     pylab.figure()
  79.     for mydata,leg in zip(totalRes,lengendlist):
  80.         pylab.plot(mydata[col],Model.predict(mydata),
  81.                    label=leg)
  83.     pylab.title("Effect of " + col)
  84.     pylab.xlabel("$\phi_{asp}$")
  85.     pylab.xlim((0.8,4.2))
  86.     pylab.ylim(0,0.20)
  87.     pylab.ylabel("v_mac/f_c")
  88.     pylab.legend()
  89.     pylab.show()
  91. # 数据采样函数
  92. def sampleOfData(rawdata,num):
  93.     """
  94.     rawdata is the total results of the experiment, a DataFrame
  95.     num is an integer, samping number.
  96.     Break the rawdata into two parts of size, (num) and (total - num).
  97.     --------
  98.     Return a tuple
  99.     --------
  100.     """
  101.     temp = rawdata.reset_index()
  102.     temp = temp.drop('ind',axis = 1)
  103.     myTrain = temp.sample(num)
  104.     myTest = temp.drop(myTrain.index)
  105.     return (myTrain,myTest)
  108. # 把数据表分解为X和y，一个是feature一个是label
  109. def getXy(data,label):
  110.     feature = data.drop(label, axis = 1)
  111.     label = data[label]
  112.     return (feature,label)
  115. # 一次拟合试验
  116. def OneTrial(rawdata,numOfSample,predModel):
  117.     """
  118.     Run one trial
  119.     rawdata is the total result, numOfSample is the sample number
  120.     and preModel is the machine learning model we use, say, MLPResgressor.
  121.     ------
  122.     Return the R2 of the model on testdata, the resultant regression model
  123.     and the corresponding test dataset.
  124.     ------
  125.     """
  126.     mytrain,mytest = sampleOfData(rawdata,numOfSample)
  127.     mytrainX,mytrainy = getXy(mytrain,'v_max/f_c')
  128.     mytestX,mytesty = getXy(mytest,'v_max/f_c')
  129.     model = predModel(activation='logistic',solver='lbfgs',hidden_layer_sizes=(2,))
  130. #    scaler = StandardScaler()
  131. #    scaler.fit(mytrainX)
  132. #    mytrainX = scaler.transform(mytrainX)
  133.     model.fit(mytrainX,mytrainy)
  134. #    mytestX = scaler.transform(mytestX)
  135.     mypre = model.predict(mytestX)
  136.     R2= model.score(mytestX,mytesty)
  138.     return (R2,mytestX,mytesty,mypre,model)
  141. # 一次模拟过程    
  142. def oneSimulation(rawdata,numOfSample,predModel,numOfTrial,res=[],totalRes=[]):
  143.     """
  144.     Run on simulation
  145.     rawdata is the total result, numOfSample is the sample number
  146.     and preModel is the machine learning model we use, say, MLPResgressor.
  147.     numOfTrial is the trial number, res and totalRes are lists to save result
  148.     Result is a tuple generating in OneTrial
  149.     ------
  150.     Return the average R^2
  151.     ------
  152.     """
  153.     print("Running total " + str(numOfTrial) + " trials using " +
  154.           str(predModel.__name__) + " algorithm")
  155.     for i in range(numOfTrial):
  156.         resOfTrial = OneTrial(rawdata,numOfSample,predModel)
  157.         res.append(resOfTrial[0])
  158.         totalRes.append(resOfTrial)
  160.     return sum(res)/len(res)
  165. # Run the simulation and then draw the figure.
  166. def plotFigure():
  167.     """
  168.     --------
  169.     Return the best model in all trials
  170.     --------
  171.     """
  173.     # 定义函数，用来绘图
  174.     def plotOne(totalRes,Model):
  175.         myX,myy = getXy(totalRes,'v_max/f_c')
  178.         for col in ["Ratio_Aspect","Ratio_AxisLoad",
  179.                     "Ratio_Tr","Ratio_Lg"]:
  180.             # plot the actual and predicted test group
  181.             pylab.figure()
  182.             pylab.plot(myX[col],myy,'k.',
  183.                    label="experiment data")
  184.             pylab.plot(myX[col],Model.predict(myX),'b^',
  185.                    label="model prediction")
  186.             pylab.title("Effect of " + col)
  187.             pylab.xlabel(col)
  188.             pylab.ylabel("v_mac/f_c")
  189.             pylab.legend()
  190.             pylab.show()
  192.         # plot the hist for difference between actual and predicted value
  193.         print("This is the difference distribution ")
  194.         pylab.figure()
  195.         pylab.hist(100 * (Model.predict(myX)-myy) / myy,15)
  196.         pylab.title("predicted value - actual value distribution")
  197.         pylab.xlabel("differences/%")
  198.         pylab.ylabel("number")
  199.         pylab.show()
  201.     # 此函数用来返回偏差最小的模型的index
  202.     def smallestIndex(mylist):
  203.         temp = mylist.copy()
  204.         res = -100
  205.         for el in mylist:
  206.             if el > res:
  207.                 res = el
  208.         return temp.index(res)
  210.     # 开始进行回归模拟过程
  211.     res2 = []
  212.     totalRes2 = []
  213.     aveR2_2 = oneSimulation(data_machine,280,neural_network.MLPRegressor,
  214.                             800,res2,totalRes2)
  215.     # 得到最佳模型的index
  216.     plotindex2 = smallestIndex(res2)
  217.     resModel = totalRes2[plotindex2][4]
  220.     # 打印信息
  221.     print("average $R^2$ of neural_network.MLPRegressor is " + 
  222.       str(aveR2_2))
  223.     print("smallest $R^2$ of neural_network.MLPRegressor is "+ 
  224.           str(totalRes2[plotindex2][0]))
  225.     # 绘图
  226.     plotOne(data_machine,resModel)
  227.     print("------------------------------------------------------------" + '\n')
  228.     # 返回结果
  229.     return resModel
  233. ## 多次拟合选最佳的模型，并保存到本地
  234. #myresmodel = plotFigure()
  237. # 读取本地保存的最佳模型
  238. myplotModel = joblib.load('model86_1l2u.model')
  239. # 绘图
  240. plotOne(data_machine,myplotModel)
  243. if __name__ == "__main__":
  244.     # 这一段是用来debug的代码
  245.     mymodel = joblib.load('model85.model')
  246.     def anaofRatioAspect():
  247.         axisload = np.average(data_machine['Ratio_AxisLoad'])
  248.         ratiolg = np.average(data_machine['Ratio_Lg'])
  249.         ratiotr = np.average(data_machine['Ratio_Tr'])
  250.         datalist = []
  251.         for axisl in [axisload-0.2,axisload,axisload+0.2,axisload+0.4,axisload+0.6]:
  252.             data_generate = []
  253.             asp = 1
  254.             while asp <=10:
  255.                 data_generate.append([asp,axisl,ratiolg,ratiotr])
  256.                 asp += 0.01
  257.             data_generate=pd.DataFrame(data_generate)
  258.             data_generate.columns=['Ratio_Aspect','Ratio_AxisLoad','Ratio_Lg','Ratio_Tr']
  259.             datalist.append(data_generate)
  260.         return datalist
  261.     myplot1 = anaofRatioAspect()
  263.     plotOneMy(myplot1,mymodel,'Ratio_Aspect',['$\psi_{axis}$ = 0','$\psi_{axis}$ = 0.2',
  264.     '$\psi_{axis}$ = 0.4','$\psi_{axis}$ = 0.6','$\psi_{axis}$ = 0.8'])
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注