Как получить среднее и стандартное значение нескольких сохраненных массивов numpy с помощью поэлементной операции

У меня есть папка с 1000 сжатыми файлами numpy (npz), представляющими результаты моделирования данных. В каждом файле есть два массива a и b с одинаковым размером, формой и типом данных. В качестве окончательного результата я хочу получить массивы среднего и стандартного отклонения по элементам a, b и c (которые я создаю в приведенном ниже примере), принимая во внимание все моделирование, то есть:

mean_a = np.mean(a1,a2,a3,...a1000)

std_a = np.std(a1,a2,a3...a1000) и др.

Мне удалось получить средние значения, но без прямой поэлементной операции. Больше всего я борюсь с ЗППП. Я пытался добавить все массивы в списки, но у меня проблема с памятью. Есть идеи, как мне действовать? См. Ниже, чего я добился до сих пор. Заранее спасибо!!

import glob
import numpy as np
import os 

simulation_runs = 10
simulation_range = np.arange(simulation_runs)

npFiles = [npFile for npFile in glob.iglob(os.path.join(outDir, "sc0*.npz"))]

a_accum = np.empty([885, 854], dtype=np.float32)
b_accum = np.empty([885, 854], dtype=np.float32)    
c_accum = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_accum  = a + a_accum
    b_accum = b + b_accum   
    c_accum = c + b_accum   

aMean = a_accum/len(simulation_range)
bMean= b_accum/len(simulation_range)
cMean = c_accum/len(simulation_range)

person renan-brso    schedule 03.01.2019    source источник
comment
Поскольку вы используете numpy, исследовали ли вы np.mean() и np.std()?   -  person kcw78    schedule 03.01.2019
comment
Да, я предполагаю, что так и должно быть, но у меня возникла проблема, как собрать все 1000 массивов вместе, чтобы использовать как np.std (), так и np.mean (). Я попытался добавить каждый массив в список, но возникает ошибка памяти. Так что, возможно, есть способ сделать это внутри цикла или даже другой способ. Я новичок в программировании, поэтому не могу представить другой альтернативы = /   -  person renan-brso    schedule 03.01.2019


Ответы (1)


Во-первых, это проще всего, если у вас есть (ssh) доступ к машине с большим объемом памяти. Может быть, ты даже сможешь обойтись без него. 885 * 854 * (1000 симуляций) * (4 байта на float32) = 2,8 ГиБ, поэтому, если вы выполняете a, b и c по отдельности, у вас должно быть достаточно памяти на разумной машине. В этом случае просто поместите их в массив и используйте np.mean и np.std:

a = np.zeros((1000,885,854), dtype=np.float32)
for run, i in enumerate(npFiles):
    a[i]=np.load(run)['scc']
amean = a.mean(axis=0)
astd = a.std(axis=0)

Аналогично для b и c.

В противном случае наиболее элегантным вариантом является сохранение данных в формате, который можно легко загрузить с ленивой загрузки. dask был специально разработан для этого, но может потребоваться некоторое время для изучения (хотя в конечном итоге оно того стоит). Вы также можете сохранить его в файлах netcat и использовать xarray как своего рода интерфейс для dask, может быть, это даже более удобно.

Если вам нужно только среднее значение, std, вы можете сделать это вручную. Формула для std:

std = sqrt(mean(abs(x - x.mean())**2))

Итак, поскольку у вас уже есть средства, процедура будет работать очень похоже на то, что вы уже делали: (непроверено)

import numpy as np
import os 

simulation_runs = 10
simulation_range = np.arange(simulation_runs)

npFiles = [npFile for npFile in glob.iglob(os.path.join(outDir, "sc0*.npz"))]

a_accum = np.empty([885, 854], dtype=np.float32)
b_accum = np.empty([885, 854], dtype=np.float32)    
c_accum = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_accum  = a + a_accum
    b_accum = b + b_accum   
    c_accum = c + b_accum   

aMean = a_accum/len(simulation_range)
bMean= b_accum/len(simulation_range)
cMean = c_accum/len(simulation_range)


a_sumsq = np.empty([885, 854], dtype=np.float32)
b_sumsq = np.empty([885, 854], dtype=np.float32)    
c_sumsq = np.empty([885, 854], dtype=np.float32)    

for run, i in enumerate(npFiles):
    npData = np.load(i)
    a = npData['scc'] 
    b = npData['bcc']
    c = a+b
    a_sumsq += (a-aMean)**2
    b_sumsq += (b-bMean)**2
    c_sumsq += (c-cMean)**2

a_std = np.sqrt(a_sumsq/(len(npFiles)-1)) # The -1 is to get an unbiased estimator
b_std = np.sqrt(b_sumsq/(len(npFiles)-1))
c_std = np.sqrt(c_sumsq/(len(npFiles)-1))
person Ewoud    schedule 03.01.2019
comment
Спасибо! Я пытаюсь заставить это работать, но все еще не понимаю. Когда интерпретатор считывает строку с a_sumsq + = (a-aMean) ** 2 ... он возвращает: '' FloatingPointError: недопустимое значение, обнаруженное в add '' '. Есть идеи, как это решить? - person renan-brso; 08.01.2019
comment
Есть ли в ваших данных NaN? Что, если вы попытаетесь распечатать все переменные непосредственно перед ошибкой, что-то не так? - person Ewoud; 15.01.2019
comment
Значения NaN появляются, когда я начинаю вычислять стандартное значение для каждого прогона. Нет Nan при вычислении среднего значения. Насчет 1-го метода ... Мне пришлось изменить прогоны с 1000 на 10000, поэтому я попытался и обнаружил ошибку памяти, к сожалению - person renan-brso; 15.01.2019
comment
Я не могу помочь вам понять, откуда берутся NaN, не глядя на данные. Но если первый метод сработал, вы можете просто взять сразу 10 фрагментов из 1000 файлов, заполнить их в некоторый массив и вычислить среднеквадратическое значение стандартных показателей ваших файлов. Это будет то же самое, что и фактический std. - person Ewoud; 15.01.2019
comment
Просто обновление: мне удалось избавиться от NaN, заменив numpy.empty на numpy.zeros (не знаю почему) ... Итак, я успешно использовал ваше второе решение! Спасибо еще раз - person renan-brso; 18.01.2019