-
Notifications
You must be signed in to change notification settings - Fork 190
New issue
Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? # to your account
Working with a large grid size #2
Comments
Interesting, I hadn't actually anticipated such an issue when I redesigned the code to solve the whole system with a single massive array. But of course, such a large array would indeed take up too much memory. The previous version of the code (v 1.0.1ish) would probably have worked for you, as it actually loops through each grid point and sets up a different Ax=b system to solve each time. The problem of course is that you have to loop through each point in the grid, which in your case is quite large. Regarding the kriging business, you're right about the A matrix, it isn't sparse. Working through subsections of the grid would be reasonable to cut down on the massive-array issue, as long as you of course incorporate all of the data each time. So you'd maybe have to loop through solving a couple of giant matrix equations with array A of shape (nx_sub, ny_sub, n, n). (Also, I think that's the right shape, although I can't recall exactly what I did right off the top of my head right now... The idea is that you have an n x n matrix defined at each of the grid points and each of those matrices is the A matrix in an Ax=b matrix equation at each grid point.) This could possibly be a good solution, as it may be faster than actually having to loop through each point in the (large) grid and still avoid using too much memory. I'll take a look at the code again later this week and see what changes (i.e., provide the option to loop through grid points, enable sub-grid workarounds, or something else) I can make to alleviate this problem. If you need to get to it soon, feeding in sub-grids might be a good way to go. As I said, as long as you feed in all the data to the sub-domains, you should be good. The code shouldn't care if your data are inside or outside of your defined grid. And thanks for your interest in the code! Really glad to see that people are actually using this. |
Actually would've been v 0.2.0 I think, which you might be able to get directly from the pypi website if you want to try to use the loop-through version |
Thanks a lot for the explanations, it makes more sens now! I just submitted to PR with a fix, that reuses the previous version with a loop and speed it up a bit. Here is some result of different optimizations for your code,
Still, even with the version 3, it would still take roughly 2 min to process a case with Edit: added case 3 above. |
I have been thinking about this, and while the above optimisations do speed things up a bit, they do not solve the fundamental problem. The scaling, as the problem size increases just doesn't seem to be right. For instance say we have The logic behind this argument, is that a number of software products are able to do Ordinary Krigging of large grids in a matter of seconds, while with the current implementation in this module, even if it were to be written in a lower level language such as C, that won't be possible, since the bottleneck is the solution of the system Do you think that is is fine, to do define the A matrix only for the nearest |
Kriging with a moving window/neighborhood is certainly a thing that can be done, although there are a few subtleties that would need to be taken into account. I haven't implemented it yet because I personally have never used that approach and it has never seemed to be an urgent need for other people I talked to in developing the code. Also, I have always used kriging for statistically contouring data; as Kitanidis (see the reference in any of the main *.py files) notes, using a moving neighborhood for contouring data is not necessarily a good idea. Of course, N=1000 points is certainly on the upper extreme of a typical kriging problem, so this may be a case in which a moving neighborhood is warranted. For kriging with a moving window to make computational sense, you'd first have to carefully select your variogram model. We'd be assuming that the variogram is spatially stationary (not a bad assumption, and we're already doing that when kriging with all of the data), and we'd have to make sure that the fit of the model variogram to the experimental variogram is better at shorter lags (shorter distances). We'd also want to choose the variogram such that the calculated weights decrease somewhat rapidly with distance and are quite small at the edges of the window. Maintaining somewhat large weights for points on the very edge of the window could possibly lead to weird discontinuities appearing in the final kriged grid (see Kitanidis; if you have a point that's just within the window and a point that's just outside the window, you can probably imagine that a jump in how they're weighted would cause some oddities in the kriging answer). This all means that you'd probably want to use a spherical, gaussian, or exponential variogram model with a relatively short range and a large sill (i.e., the semivariance plateaus to a large value relatively quickly, or equivalently the spatial covariance drops off quickly) or a linear variogram with a relatively large slope (again, semivariance grows quickly/spatial covariance drops off quickly). In real life, the details of the variogram model don't make that big a difference, but I think the choice of variogram model and parameters might be more important in this case. Also, there would be no guarantees that the kriging equations at each grid point would have anything in common when kriging with a moving window. When kriging with all the data, at least the A matrix is the same everywhere, so we can just invert that once. When kriging with a moving window, there's no guarantee that A is the same, or even the same size, at any two grid points. We would then have to set up a separate kriging system at each grid point, and I'm not sure how that would work as a vectorized operation in Python. I also don't know how looping through a large grid and setting up/solving the entirety of a smaller kriging system at each point would compare computationally to solving for everything at once with a single direct call to underlying C/Fortran code. I'm sure you have more thoughts on how to make the actual computations work than I do... Anyways, adding the capability to krige with a moving neighborhood wouldn't necessarily be a bad thing. Ideally, we'd have the option to either use the nearest N points for kriging or use all points within a distance of X from the current grid point (or use the nearest points within a distance of X, up to a total number N). The former would at least keep the size of the kriging system consistent from grid point to grid point, but both features would be desirable. |
Thank you for these very helpful comments. I am not very familiar with the mathematical background for kriging so it is good to know. I completely agree that this should be an optional feature that can be activated so it does not break the current functionality. In terms of the implementation, I was thinking to do the following for the Ordinary Kriging on a moving window.
Let me know whether you think this makes sens for kriging with a moving window. The next step would be to port this version with a loop to Cython for speed. Having looked at various optimization for this, if we want something that is fast, not too difficult to read and doesn't use lots of memory, porting the loop to C with Cython seem to be the optimal solution. It does not change the the installation procedure with |
Hello,
thanks you so much for working on this and making it available on github!
I would like to use this module with a moderately large grid size of the order of nx=1000, ny=1000 and with a number of points n=100. Given that
UniversalKriging
solves the linear system,A X = b
, whereA
is, I believe, an array of shape(nx, ny, n, n)
, this currently requires roughlynx*ny*n**2*8/1e9= 80 GB
of RAM simply to define that matrix in 64bit , and subsequently I'm running out of memory on my laptop.I was wondering if I was possible to optimise this code so it is more memory efficient? For instance, the option would be to port the
execute
method inuk.py
to Cython, so it uses a lower number of temporary numpy arrays for indexing, etc., which would same some memory. However, at the end of the day, it's still necessary to define that sameA
matrix (which is not sparse as far as I can tell) for the linear system. I suppose, it's not possible to do a domain decomposition on the grid (decompose it in regions and do the calculation region by region), right? Is there any other things that could be attempted? I can do some optimisation , I would just need some advice about where to start since I don't know much about Kriging.Thanks,
The text was updated successfully, but these errors were encountered: