<div dir="ltr"><div>Hi,<br><div><br>Looks like I was a little confused. It appears that the nan* versions of functions in numpy just substitute the NaNs in a copy of the original array and so are just convenience methods. I was imagining that they were optimized and handling the NaNs at a lower level. It looks like the "bottleneck" package tries to do this for nansum, nanprod, etc, but I don't know if it is able to take advantage of SSE or not.<br><br></div><div>Anyways, maybe if I elaborate someone can offer suggestions.<br></div><div><br></div><div>My application is: given a large N-dimensional array, and a selected (N-1)-dimensional slice of it, compute the Pearson correlation of that slice versus all other (N-1)-dimensional slices of the array, omitting NaN's from the calculation.<br><br></div><div>Ignoring NaN's for a moment, here is the slow and obvious way to do it:<br></div><div><br></div> from scipy.stats import pearsonr<br> import numpy as np<br> def corrs(data,index,dim):<br> seed = data.take(index,dim)<br> res = np.zeros(data.shape[dim])<br> for i in range(0,data.shape[dim]):<br> res[i] = pearsonr(seed, data.take(i,dim))[0]<br> return res<br><br></div>Doing all the math by hand and using einsum, there is an extremely fast (though fairly cryptic) way of doing this:<br><div><br> def corrs2(data, index, axis):<br> seed = data.take([index], axis=axis)<br><br> sdims = range(0,seed.ndim)<br> ddims = range(0,data.ndim)<br> sample_axes = np.array([i for i in ddims if i != axis])<br><br> seed_mean = np.einsum(seed, sdims, []) / seed.size<br> data_mean = np.einsum(data, ddims, [axis]) / seed.size<br> data_mean.shape = tuple(data.shape[i] if i == axis else 1 for i in ddims) # restore dims after einsum<br><br> seed_dev = np.einsum(seed-seed_mean, sdims, sdims)<br> numerator = np.einsum(seed_dev, ddims, data, ddims, [axis])<br> numerator -= np.einsum(seed_dev, ddims, data_mean, ddims, [axis])<br><br> denominator = np.einsum(data, ddims, data, ddims, [axis])<br> denominator += -2.0*np.einsum(data, ddims, data_mean, ddims, [axis])<br> denominator += np.sum(data_mean**2, axis=sample_axes) * seed.size<br> denominator *= np.einsum(seed_dev**2, sdims, [axis])<br> denominator = np.sqrt(denominator)<br><br> return np.clip(numerator / denominator, -1.0, 1.0)<br><br></div><div>It also doesn't need to make a copy of the array. Re-introducing the requirement to handle NaNs, though, I couldn't find any option besides making a mask array and introducing that explicitly into the calculations. That's why I was imagining an optimized "naneinsum."<br><br></div><div>Are there any existing ways of doing sums and products on numpy arrays that have a fast way of handling NaNs? Or is a mask array the best thing to hope for? I think that for this problem that I could transpose and reshape the N-dimensional array down to 2-dimensions without making an array copy, which might make it easier to interface with some optimizing package that doesn't support multidimensional arrays fully.<br><br></div><div>I am fairly new to making numpy/python fast, and coming from MATLAB am very impressed, though there are a bewlidering number of options when it comes to trying to optimize.<br></div><div><br></div><div>Thanks,<br></div><div>Chris<br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Apr 19, 2016 at 11:12 AM, Chris Barber <span dir="ltr"><<a href="mailto:barberchris01@gmail.com" target="_blank">barberchris01@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Is there any interest in a nan-ignoring version of einsum a la nansum, nanprod, etc? Any idea how difficult it would be to implement?<span class="HOEnZb"><font color="#888888"><br><br></font></span></div><span class="HOEnZb"><font color="#888888">- Chris<br></font></span></div>
</blockquote></div><br></div>