向量化计算
- 使用np的向量化函数将一次只能返回单个标量的函数，向量化成能接受定制shape数组且可以指定类型的返回。
- 首先观察函数如：
def HSS_one_threshold(obs, pre, threshold=0.1):
'''
HSS - Heidke skill score
Args:
obs (numpy.ndarray): observations
pre (numpy.ndarray): pre
threshold (float) : threshold for rainfall values binaryzation
(rain/no rain)
Returns:
float: HSS value
'''
hits, misses, falsealarms, correctnegatives = prep_clf(obs=obs, pre = pre,
threshold=threshold)
HSS_num = 2 * (hits * correctnegatives - misses * falsealarms)
HSS_den = (misses**2 + falsealarms**2 + 2*hits*correctnegatives +
(misses + falsealarms)*(hits + correctnegatives))
return HSS_num / HSS_den这个HSS函数只能计算单个时次的，如果我想向量化计算不同时刻的hss评分，就只能使用for循环，但显然这不是我的风格。而向量化函数很简单只需要一句：
hss = np.vectorize(HSS_one_threshold,excluded=['threshold'],signature='(m,n),(m,n)->()')np.vectorize这个函数很简单，但定制自己的向量化函数核心却是signature=’(m,n),(m,n)->()’，signatrue标签，这个异常晦涩难懂且很少有人讲清楚这个东西是干什么的，具体请参考文档，我假设你详细阅读完文档后还是不懂，那么接下来请看我的表演，首先signature中有两个核心定义即numpy数组的迭代维和核心计算维，我们假设有一个三维数组（time,m,n）迭代维默认从数组最左边的维开始，核心计算维从右边开始，如果我想通过向量化的hss函数计算time个时间步的hss得分并返回到一个列表中，那么就如我给出的函数所示，计算维为（m,n）那么向量化后的函数接受的参数就是三维数组中的后两维，然后在迭代维time上并行计算，因为numpy底层是用c计算，所以这个签名相当于一个给c接口的参数。而 ->() 表示计算返回的格式（）代表标量当然这个还能根据任务自定义维度，甚至可以跳过中间某个维或者自定义广播形式非常复杂，但普通需要for循环的计算均可以用signture消灭，非常畅快。比如计算某个场与另一个场的相关系数，如果分辨率很小据我观察大多数python计算程序都是用for来遍历这样直接爆炸，但使用signuture自定义可以非常完美的解决这个问题。
利用多核进行多进程计算
显然根据经验io密集型使用多线程，计算密集使用多进程，我曾尝试过多线程开文件套多进程计算但碰到很多问题，但针对本问题使用多进程即可，或者多进程开文件套多进程计算，这需要修改linux默认限制打开文件数。如果我想同时计算不同的如csi，pod的得分就可以使用如下程序：
dbz = ['10dbz','20dbz','30dbz','40dbz','50dbz']
thdlst = [10,20,30,40,50]
for thd_number,key in enumerate(dbz):
"""
多进程实现并行运算
"""
func_lst = [hss,csi,far,pod]
#进程级的字典
share_dict = multiprocessing.Manager().dict()
#进程级的锁
share_lock = multiprocessing.Manager().Lock()
def multiprocess(func,thd):
func = func
print(func._name_)
#print(dir(func))
thd = thd
#全局变量
global unet,rover,tagan
temp_unet = func(unet[:,:,:480],unet[:,:,480:],threshold=thd)
temp_rover = func(rover[:,:,:480],rover[:,:,480:],threshold=thd)
temp_tagan = func(tagan[:,:,:480],tagan[:,:,480:],threshold=thd)
data30 = [temp_unet[4],temp_rover[4],temp_tagan[4]]
data60 = [temp_unet[9],temp_rover[9],temp_tagan[9]]
data30 = pd.Series(data30,['uent','rover','tagan'])
data60 = pd.Series(data60,['uent','rover','tagan'])
#使用当前函数名作为进程字典的key
key = func._name_
share_lock.acquire()
share_dict[key] = (data30,data60)
share_lock.release()
print('the process of %s is finisthed '%(key))
thd = thdlst[thd_number]
process_lst =[]
p = Pool(4)
for func in func_lst :
#p.apply_async(multiprocess, args=(multiprocess,thd))
process = multiprocessing.Process(target=multiprocess,args=(func,thd))
process_lst.append(process)
for pss in process_lst:
pss.start()
for pss in process_lst:
pss.join()首先给一个进程级别的共享列表或者字典用来储存所有返回的得分值，在计算过程中使用相同的数据不同的技巧评分所以将数据在函数中定义为全局变量，不同进程只访问内存不修改，但在写入字典的时候给一个锁防止冲突。这样如果计算指标在几十个以上时数据量在工业的TB或者PB级时，就可以利用服务器的多核把cpu拉满。需要注意的是将函数向量化后新的对象就没有了__name__属性，需要在向量化函数时setattr一个name属性
一些需要注意的小技巧
- 一个进程process或者进程池pool使用完后应该使用close关闭，清除多余内存，而terminte和kill并不会清理内存，容易造成内存泄漏，使用del和gc手动清除大的变量。
- 在计算非对称的数组时可以手动boardcast数组进行向量化的计算，如同学在云检测选点时需要用到这个方法。
使用dask并行计算
具体就是delayed，建图，计算，需要注意的是list返回需要一个节点
代码如下，加速站点数据与格点数据交互运算，当然它也可以对数组分块在compute，与dataframe交互但美这个delayed好用
import xarray as xr
data = xr.open_dataset('/media/workdir/hujh/precipation/reploat/cyclone.nc')
aeradata = data.sel(lat = data.lat[(data.lat>=17.75)&(data.lat<=28.75)],lon=data.lon[(data.lon>=105)&(data.lon<=126.25)])
aeradata = aeradata.rename({"__xarray_dataarray_variable__": "cyc"})
you = aeradata.cyc.sum(['lon','lat'])
qixuan = you[you>0].time
import pandas as pd
station = pd.read_csv('/media/workdir/hujh/precipation/data/neworg_data.csv')
station = station[['id','datetime','pre']]
import numpy as np
station.pre = station.pre.where(station.pre>=50,np.nan)
predata = station
import meteva.base as meb
station = meb.read_station(meb.station_国家站)
station = station[['id','lon','lat']]
keyid = station[(station.lat>=17.75)&(station.lat<=28.75)&(station.lon>=105)&(station.lon<=126.25)].id.values
predata = predata[predata.id.isin(keyid)]
predata = predata.dropna()
predata['datetime'] = pd.to_datetime(predata.datetime)
#两个数据集的时间交集
sttime = predata[['datetime','pre']].copy()
sttime = sttime.set_index('datetime')
stime=sttime[~sttime.index.duplicated(keep='first')]
x = qixuan.time.to_dataframe()
result = pd.concat([x,stime], axis=1, join='inner')
aeradata = aeradata.sel(time =result.index)
predata = predata[predata.datetime.isin(result.index)]
print('len of time--->',len(result.index))
print('process is ok,after to fucker')
def cacualte(df):
global aeradata
global station
presum = 0
presumlst =[]
datetime = df.datetime.iloc[0]
if True:
cycdata = aeradata.sel(time = datetime)
#使用dask并行改造串行
#计算图最终需要一个节点来完成计算
from dask import delayed
for index, row in df.iterrows():
pre = 0
def sumstationpre(row):
print(row.id)
stationid = station[station.id==row.id]
lon,lat = stationid['lon'].values[0],stationid['lat'].values[0]
print(lon,lat)
data = cycdata.cyc
#lat的方向犯了
latvalues = data.lat.values[::-1]
data = data.assign_coords(lat = latvalues)
datas = data.sel(lon=lon,lat = lat,method="nearest")
print(datas)
datas = datas.values
if datas != 0:
pre = row['pre']
else:
pre = 0
return pre
def get_list(lst):
return lst
pre = delayed(sumstationpre)(row)
presumlst.append(pre)
lst = delayed(get_list)(presumlst)
result = lst.compute()
result = np.sum(result)
presum = result
else:
presum = 0
return presum
print('begin groupby')
presum = predata.groupby('datetime').apply(cacualte)
print(presum)
presum.to_csv('qixuanpre.csv')mpipy4
这个我要学习一下，我准备改造这个串行下载数据的代码
#导入模块
import datetime
from metpy.units import units
from siphon.simplewebservice.wyoming import WyomingUpperAir
# 设置下载时段（这里是UTC时刻）
start = datetime.datetime(1979, 1, 1, 0)
end = datetime.datetime(2018, 12, 30, 0)
datelist = []
while start<=end:
datelist.append(start)
start+=datetime.timedelta(hours=12)
# 选择下载站点（以北京为例）
stationlist = ['45004','59758','58968','59316']
# 批量下载
for station in stationlist:
for date in datelist:
try:
df = WyomingUpperAir.request_data(date, station)
df.to_csv(station+'_'+date.strftime('%Y%m%d%H')+'.csv',index=False)
print(f'{date.strftime("%Y%m%d_%H")}下载成功')
except Exception as e:
print(f'{date.strftime("%Y%m%d_%H")}下载失败: {e}')
pass先记录一下笔记
首先它是一个消息接口
可以看作multiprocess里pipe，manger等处理消息的管道等
最主要的是可以和fortran以及c++交互
对于上述代码可以用多进程嵌套多线程改造我试试
mpi实战 我将第一个循环用4个进程干掉，第二个循环用线程池干掉
这个很简单得到当前进程号的序号号取站点就可以了
stationlist = ['45004','59758','58968','59316']
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# 批量下载
station = stationlist[rank]
print('这个进程将执行第{st}站点的下载'.format(st=station))这个不需要规约函数什么的，之前看的一篇帖子太硬了找不到了什么二叉树规约各种降时间复杂度，给我秀的双眼发晕这就是码农强者吗？
看看阻塞和非阻塞
也就是主程序不会被挂起，消息发出去后进一步操作，这个展开实在太复杂了，贴个链接有时间看看：
最终的代码如下，可行思路就是多进程嵌套多线程完成，这里有个思考mpi对比multiprocess似乎更偏向于底层，比如multiprocess里都是高级封装
#导入模块
import datetime
#from metpy.units import units
from siphon.simplewebservice.wyoming import WyomingUpperAir
from concurrent.futures import ThreadPoolExecutor, wait, ALL_COMPLETED, as_completed
from mpi4py import MPI
# 设置下载时段（这里是UTC时刻）
start = datetime.datetime(1979, 1, 1, 0)
end = datetime.datetime(2018, 12, 30, 0)
datelist = []
while start<=end:
datelist.append(start)
start+=datetime.timedelta(hours=12)
# 选择下载站点（以北京为例）
stationlist = ['45004','59758','58968','59316']
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
# 批量下载
station = stationlist[rank]
print('这个进程将执行第{st}站点的下载'.format(st=station))
def download(date,station):
try:
df = WyomingUpperAir.request_data(date, station)
print('开始下载时间',date)
df.to_csv('/media/workdir/hujh/sounding/data/'+station+'/'+station+date.strftime('%Y%m%d%H')+'.csv',index=False)
#print(f'{date.strftime("%Y%m%d_%H")}下载成功')
except Exception as e:
#print(f'{date.strftime("%Y%m%d_%H")}下载失败: {e}')
pass
return 1
executor = ThreadPoolExecutor(max_workers=10)
thread_list = []
for date in datelist:
thrd = executor.submit(download, date,station)
thread_list.append(thrd)
result_list = []
for task in as_completed(thread_list):
result_list.append(task.result())
print('最终下载的文件数量--》',len(result_list))xarray向量化计算
官网文档实在太垃圾，其实就是apply_unfunc这个函数计算个趋势什么的，废话不多说直接贴代码
from scipy import stats
def linear(x):
slope, _, _, p, _ = stats.linregress([i for i in range(len(x))], x)
return (slope,p)
q1_trend = xr.apply_ufunc(
linear,
q1year,
input_core_dims = [['year']],
output_core_dims = [[],[]],
vectorize=True
)这个就是在year这个核心维度上计算，在其他维度上广播，返回的有两个一个是斜率一个是p值，还是比for快多了
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
