Calculation of derived coords points and bounds is always lazy. #2604

pp-mo · 2017-06-14T13:34:06Z

Addresses #2586

corinnebosley

Very recently made these changes to exactly the same code section. Tests added at that time should cover this PR as well.

Provided those recently-added tests pass, I am happy with this change.

pp-mo · 2017-06-14T14:03:09Z

Very recently made these changes ...

Yes, I also removed the line nd_values = np.array(nd_values) from the _nd_bounds method.
Which is extremely important in this context !

pp-mo · 2017-06-15T10:24:34Z

Errors have come in because some factory calculations are not viable on dask arrays.

~~Please wait for me to fix !*~~

pp-mo · 2017-06-22T11:44:56Z

I think this is finally sorted.
There are many commits, and much to+fro getting this done and fixing outstanding problems it uncovered (!)
As a summary:

fixed basic factory _nd_points + _nd_bounds methods to always yield lazy arrays, even if core data is real : This is the key bit that was required to ensure deferred calculation of derived coords even with real source data, which caused some noticeable slownesses. : LazyArrays used to do this for us, now use Dask, always, for same effect
rewrote some _derive methods to ensure all those calcs are now "dask capable"
added integration tests for OceanSigmaZFactory, covering various usage not handled by the existing unit testing : data with a time dimension; additional non-derived dimensions; awkward transpositions.
removed _shape and _dtype methods : no longer needed since Dask replaced LazyArray
fixed several usages of array.reshape(list(array.shape).append(1)) which was plain wrong (!!)
fixed behaviour testing for nbounds mismatch in iris.tests.test_hybrid, previously wrong due to the above bug
accepted a few more NaN-->masked changes in CMLs

pp-mo · 2017-06-22T11:46:36Z

STATUS:
It's good for review now 😉

Currently waiting a fix for testing problems since we provided numpy 1.13.
Will rebase when that is done.

I will also get on + post some info on the performance improvement that motivated this ...

bjlittle · 2017-06-22T11:49:50Z

@pp-mo See #2616 😉

bjlittle · 2017-06-22T12:10:37Z

lib/iris/aux_factory.py

-                    shape[i] = size
-        return shape
-
-    def _dtype(self, arrays_by_key, **other_args):


@pp-mo 😱 Oh my gosh ... this isn't used anywhere! Who knows when it was relevant and useful ... good spot!

Both _shape and _dtype were needed when making LazyArray-s, so you could tell it what the shape and type the result would be. Biggus or dask does all that for you, of course.

bjlittle · 2017-06-22T12:12:15Z

lib/iris/aux_factory.py

            nd_values_by_key[key] = nd_values
        return nd_values_by_key

-    def _shape(self, nd_values_by_key):


@pp-mo This was only used in OceanSigmaZFactory.make_coord(), so you must have refactored that away ...

Effectively the same calc now exists inside OceanSigmaZFactory._derive(), but now it is done by making an array of all the dependency shapes + taking max size over each dim.

bjlittle · 2017-06-22T13:27:30Z

lib/iris/aux_factory.py

        transpose_order = [pair[0] for pair in sorted_pairs] + [len(dims)]
-        bounds = coord.core_bounds()
+        bounds = coord.lazy_bounds()
        if dims:


@pp-mo There is no equivalent check here as per line 221 to check for a nop transpose, given an increasing dims order ... is it relevant to add that here also? Since we're in this space ...

bjlittle · 2017-06-22T13:44:25Z

lib/iris/aux_factory.py

-                    orography = orography_pts.reshape(
-                        orography_pts_shape.append(1))
+                    bds_shape = list(orography_pts.shape) + [1]
+                    orography = orography_pts.reshape(bds_shape)


@pp-mo Not your doing, but ... the comment on line 408 is no longer relevant. We don't have closures anymore now that we've purged the _LazyArray implementation ... could you please remove all such closure related comments, as they are no longer correct, thanks!

Personally I never liked all these called-only-once functions at all.
So given your prompt I've just removed them all -- hope that suits!

bjlittle · 2017-06-22T14:21:21Z

lib/iris/aux_factory.py

-        nsigma_slice[index] = slice(0, int(nd_points_by_key['nsigma']))
-        nsigma_slice = tuple(nsigma_slice)

+        nsigma, = nd_points_by_key['nsigma']


@pp-mo Subtle unpacking using nsigma, ... but [nsigma] = nd_points_by_key['nsigma'] is an alternative, less subtle pattern ... you're choice, it's a minor point.

Good spot, I do prefer.

bjlittle · 2017-06-22T14:30:33Z

lib/iris/aux_factory.py

                              nd_points_by_key['zlev'],
-                              points_shape,
-                              nsigma_slice)
+                              nsigma,


@pp-mo You could just pass in nd_points_by_key['nsigma'] and do away with the local nsigma ... why did you unpack it? Just for convenience? I can only see it being used on line 866 below ...

Because I wanted it to be seen as an "extra" argument, i.e. a normal Python value, and not another dependency array.

bjlittle · 2017-06-22T14:40:51Z

lib/iris/aux_factory.py

+                zlev, nsigma, coord_dims_func):
+        # Calculate the index of the 'z' dimension in the inputs.
+        # Get the cube dimension...
+        i_levels_cubedim, = coord_dims_func(self.dependencies['zlev'])


@pp-mo Again [i_levels_cubedim] rather than i_levels_cubedim, ... your choice.

pp-mo · 2017-06-22T14:55:26Z

Here's the performance demo.

First I generated some hybrid-height data regridded from iris.tests.stock.realistic_4d to 5* higher horizontal resolutions, getting a cube of dims (time:6, level=70, y=500, x=500).
(Not as easy as it sounded ! See : https://gist.github.com/pp-mo/f393cc2be06d0b457e5dd114a5d8b6e3 )

Then run :

import datetime
import sys

import iris


if __name__ == '__main__':
    # Get the data.
    file_path = 'big_hybrid_height_cube.nc'
    cube = iris.load_cube(file_path)

    print 'cube shape = ', cube.shape

    dependency_names = ('atmosphere_hybrid_height_coordinate',
                        'sigma',
                        'surface_altitude')

    def test():
        for coord_name in dependency_names:
            coord = cube.coord(coord_name)
            print '  {} coord is "{}"'.format(
                coord.name(),
                'lazy' if coord.has_lazy_points() else 'real')

        # Time the altitude calculation.
        start_time = datetime.datetime.now()
        alts = cube.coord('altitude')
        stop_time = datetime.datetime.now()
        time_taken = (stop_time - start_time).total_seconds()
        print '  TIMED: altitude coordinate fetch took {} secs.'.format(time_taken)

    print 'With lazy dependency coords ...'
    test()

    # Realise the aux data.
    for coord in cube.aux_coords:
        _ = coord.points
        _ = coord.bounds

    print 'With realised coords ...'
    test()

Results:

$ git checkout master
$ python hybrid_height_timing.py 
cube shape =  (6, 70, 500, 500)
With lazy dependency coords ...
  atmosphere_hybrid_height_coordinate coord is "lazy"
  sigma coord is "lazy"
  surface_altitude coord is "lazy"
  TIMED: altitude coordinate fetch took 0.010326 secs.
With realised coords ...
  atmosphere_hybrid_height_coordinate coord is "real"
  sigma coord is "real"
  surface_altitude coord is "real"
  TIMED: altitude coordinate fetch took 0.752618 secs.

$ git checkout derived_coords_lazycalc
$ python hybrid_height_timing.py 
cube shape =  (6, 70, 500, 500)
With lazy dependency coords ...
  atmosphere_hybrid_height_coordinate coord is "lazy"
  sigma coord is "lazy"
  surface_altitude coord is "lazy"
  TIMED: altitude coordinate fetch took 0.010303 secs.
With realised coords ...
  atmosphere_hybrid_height_coordinate coord is "real"
  sigma coord is "real"
  surface_altitude coord is "real"
  TIMED: altitude coordinate fetch took 0.014665 secs.

itpp@eld238: /home/h05/itpp/git/iris/iris_main
$

CONTEXT NOTE:
The problem is with getting a derived coord (e.g. 'altitude') when all underlying dependencies have real data ...
Old code (current master) works with "coord.core_xxx()", so it calculates when creating the derived coord. New code (here) uses "coord.lazy_xxx()", so calculation is always deferred.

So in this case, the delay normally only happens after it has fetched the derived coords data.

When all dependencies are realised, you get a delay every time you fetch cube.coord('altitude'), e.g. if you print the cube.

bjlittle · 2017-06-22T15:03:34Z

lib/iris/aux_factory.py

+        derived_cubedims = self.derived_dims(coord_dims_func)
+        i_levels_dim = i_levels_cubedim - sum(
+            i_dim not in derived_cubedims
+            for i_dim in range(i_levels_cubedim))


@pp-mo Would you buy into renaming i_levels_cubedim to zlev_dim, and i_levels_dim to zlev_index or zlev_offset.

There's a lot going on here, and using i_levels_... feels (to me) like its one level if indirection that could be avoided ...

That naming was an attempt to distinguish between the dependency arguments, which are all arrays of the same dimensionality, and the "extra" args which as ordinary Python values.
Hence the "i_" -- it means 'integer', not indirection.

The problem is, we have two contexts for 'dimension' or 'index' : The original cube, and the dependency arguments. That's why I need to calculate 'i_levels_dim' from 'i_levels_cubedim'.

I'll think on ...

Isn't .index your friend here?

Isn't .index your friend here?

Yes, thanks.
I've now returned to that, the way it was done in the original code.

bjlittle · 2017-06-22T15:39:05Z

lib/iris/aux_factory.py

+            [el.shape
+             for el in (sigma, eta, depth, depth_c, zlev)
+             if el.ndim])
+        result_shape = list(np.max(allshapes, axis=0))


@pp-mo This only works if all the elements have the same number of shape dimensions ... which must be the case, right? ... given the output from _remap and _remap_with_bounds, which aligns the dimensionality of everything or injects 0-d scalars for missing coordinates.

Just looking for re-assurance (and convince me otherwise) but should we ensure that allshapes has equal length elements before doing the np.max ?

Yes, I believe it is guaranteed by _nd_points and _nd_bounds that the results have all-same dimensions, given they are called with the same 'ndim'.

I will try to amend some docstrings (_nd_xxx and _remap_xxx) to make this clearer : it's obvious that _remap_xxx should have docstrings reallly ...

bjlittle · 2017-06-22T15:51:22Z

lib/iris/aux_factory.py

-
+        nsigma_levs = eta + sigma * (da.minimum(depth_c, depth) + eta)
+        # Expand to full shape, as it may sometimes have lower dimensionality.
+        ones_full_result = np.ones(result_shape, dtype=np.int16)


@pp-mo Can't we use nsigma_levs.dtype here instead of np.int16 ...

Also, shouldn't this be da.ones(...) ? To keep all things lazy ...

shouldn't this be da.ones(...) ? To keep all things lazy

I thought I'd only use dask where needed, and I thought it wasn't here because "nsigma_levs" is always lazy, so result is "dask * numpy". However, you just reminded me that this could be creating a large real array, so I'll change !!

Can't we use nsigma_levs.dtype here instead of np.int16

Yes, we can, probably better...

bjlittle · 2017-06-22T15:59:44Z

lib/iris/aux_factory.py

+        # Expand to full shape, as it may sometimes have lower dimensionality.
+        ones_full_result = np.ones(result_shape, dtype=np.int16)
+        ones_nsigma_result = ones_full_result[z_slices_nsigma]
+        result_nsigma_levs = nsigma_levs * ones_nsigma_result


@pp-mo This is getting a tad abstract ... but why is there a need to do nsigma_levs * ones_nsigma_result ?

Is it not sufficient just to do:

nsigma_levs = eta + sigma * (da.minimum(depth_c, depth) + eta) zlev = zlev * da.ones(result_shape, dtype=nsigma_levs.dtype) result = da.concatenate([nsigma_levs, zlev[z_slices_rest]], axis=i_levels_dim)

Is it not sufficient just to do

No, because of the way it has to work with possibly-missing dependencies.

From the CF equation:

k <= nsigma:: z(n,k,j,i) = eta(n,j,i) + sigma(k)*(min(depth_c,depth(j,i))+eta(n,j,i)) k > nsigma:: z(n,k,j,i) = zlev(k)

From _check_dependencies, we always have zlev but maybe only one of sigma or eta.
Thus we always have a 'k' (vertical) dimension in the dependency dims (but not always any 'n' (time) dimension).

If sigma is missing, then the "main" nsigma_levs calculation yields just 'eta(n, 1, j, i)'
-- the 'k' dimension is a 1 because it isn't present in the original eta array.
For the concatenation, this must get "replicated" up to (n, nsigma, j, i).
If not, in numpy you get a different (wrong) result shape
-- and in dask that causes an actual error.
The original code had the assignment result[nsigma_slice] = , which does the right thing "automatically", by broadcasting.

bjlittle · 2017-06-22T16:03:15Z

@pp-mo Okay, I'm done for now ... over to you!

pp-mo · 2017-06-23T16:30:32Z

@bjlittle thanks for all the attention !
I'm hoping the latest commit addresses all your earlier points.
This will presumably need at least respinning, possibly rebasing, to get tests passing following #2616

bjlittle

@pp-mo Awesome effort! Thanks for sticking with it!

I'm 👍 for the PR, just rebase once #2616 is merged, and then we can squash and merge!

…working.

pp-mo · 2017-06-26T09:36:54Z

@bjlittle just rebase once #2616 is merged

Done that, let's see if it passes ... 🤞

… warnings tests.

pp-mo · 2017-06-27T14:36:30Z

"Python3 deprecation undermining warnings tests"

That was a bit nasty.
Test code like this was failing in Python 3.5 ...

msg = 'sample text'
with self.assertRaisesRegexp(UserWarning, msg):
    warnings.simple_filter('error')
    call()

... because assertRaisesRegexp() is deprecated in favour of assertRaisesRegex() (with no final "p").
So it was getting the deprecation warning instead of the intended warning.

pp-mo · 2017-06-27T15:14:45Z

At last it passes 🥂
-after many, many, many re-spins for Travis failures ☹️

@pelson are you happy with this, in absence of @bjlittle ?
I really don't want to have to touch this again!

pelson · 2017-06-28T09:44:48Z

@pp-mo - Given the spin-up cost, I honestly can't justify the effort of me merging this today. I completely appreciate your desire for it not to go stale, and am hopeful that @bjlittle will be able to merge tomorrow when he returns. If not, I'll clear a few hours and try to get the ball rolling on it myself in the next few days. Work for you?

corinnebosley · 2017-06-28T09:50:27Z

@pp-mo I can review this. Just give me a bit of time to check what you have done since I was last in.

corinnebosley

Aside from wanting to know what the slices tuple whatever does, I am happy with this.

corinnebosley · 2017-06-28T10:56:07Z

lib/iris/aux_factory.py

+        # Make a slice tuple to index the remaining z-levels.
+        z_slices_rest = [slice(None)] * ndims
+        z_slices_rest[z_dim] = slice(int(nsigma), None)
+        z_slices_rest = tuple(z_slices_rest)


Could you please explain to me what this bit (L.814 to L.820) does? I don't understand.

The z_slices_nsigma thing is a tuple of keys to extract the first 'nsigma' z-levels,
i.e. data[z_slices_nsigma] would be something line data[:, ..., :, :nsigma, :, ..., :]
- recall that here, 'nsigma' is just a number.

This is exactly what in the original code is called nsigma_slice.
I changed over to calculating that within the _derive call, instead of passing in in from the make_coord routine, as I thought it was much clearer to have derivation + use of it in the same place.

Meanwhile, the subsequent bit z_slices_rest is to select the "remaining" z levels, i.e. those beyond the first nsigma :
that is, if data[z_slices_nsigma] does something like data[:, ..., :, 0:nsigma, :, ..., :],
then data[z_slices_rest] is like data[:, ..., :, nsigma:-1, :, ..., :],
In the new calc, we need to extract over the 'remaining' levels to get the concatenate right.

Is that any clearer ?

corinnebosley · 2017-06-28T11:55:35Z

@pp-mo What do you want me to do about this? I am happy to merge, but I can leave it if you would like someone more qualified to check it over first.

pp-mo · 2017-06-28T12:53:34Z

Hi @corinnebosley sorry for slow response -- I've been at a meeting most of this morning.

I am happy to merge, but I can leave it if you would like someone more qualified to check it over first.

I think @bjlittle provisionally approved this anyway according #2604 (review)

That should cover everything except the latest commit, which fixed the warning/error testing bug.

corinnebosley · 2017-06-28T13:18:43Z

@pp-mo I thought that might have been the case, but I checked over everything anyway and couldn't spot any issues. Bearing that in mind, and given that the tests are passing, I will merge this in one hour (at 15:18) unless I get any objections in the meantime.

corinnebosley · 2017-06-28T14:18:59Z

@pp-mo Boom!

pp-mo · 2017-06-28T15:13:45Z

@corinnebosley thanks !

…ools#2604) * Integration tests for OceanSigmaZFactory -- lazy cases not currently working. * Refactor osz (not yet lazy). * Add test with extra cube dims. * Fixed calculation; all working except lazy integration tests. * Enable all-lazy operation; all tests working. * Fix nasty misuses of list.append. * Adjust testing for fixes to aux_factory code. * CML changes: missing altitude points NaN --> mask. * Clarify need for some integration testcases. * Review changes. * Reroute assertRaisesRegexp to prevent Python3 deprecation undermining warnings tests.

pp-mo mentioned this pull request Jun 14, 2017

Fix speed problem with area_weights #2586

Closed

corinnebosley approved these changes Jun 14, 2017

View reviewed changes

corinnebosley mentioned this pull request Jun 21, 2017

Use gallery as basis for performance testing (2 man days' effort) #2573

Closed

pp-mo mentioned this pull request Jun 21, 2017

cube.transpose is an in-place operation #2615

Closed

pp-mo force-pushed the derived_coords_lazycalc branch from 718fad1 to b84ff92 Compare June 22, 2017 11:15

bjlittle self-assigned this Jun 22, 2017

bjlittle reviewed Jun 22, 2017

View reviewed changes

bjlittle added this to the v2.0 milestone Jun 23, 2017

bjlittle approved these changes Jun 23, 2017

View reviewed changes

pp-mo added 3 commits June 26, 2017 10:35

Integration tests for OceanSigmaZFactory -- lazy cases not currently …

0847a4c

…working.

Refactor osz (not yet lazy).

431e41e

Add test with extra cube dims.

90b2573

pp-mo added 7 commits June 26, 2017 10:35

Fixed calculation; all working except lazy integration tests.

d434705

Enable all-lazy operation; all tests working.

e1b021c

Fix nasty misuses of list.append.

507a4df

Adjust testing for fixes to aux_factory code.

f962c1b

CML changes: missing altitude points NaN --> mask.

884faa0

Clarify need for some integration testcases.

48cf11f

Review changes.

8dab554

pp-mo force-pushed the derived_coords_lazycalc branch from 01bd13d to 8dab554 Compare June 26, 2017 09:35

Reroute assertRaisesRegexp to prevent Python3 deprecation undermining…

3a4a536

… warnings tests.

pp-mo force-pushed the derived_coords_lazycalc branch from 507a70f to 3a4a536 Compare June 27, 2017 13:59

corinnebosley approved these changes Jun 28, 2017

View reviewed changes

corinnebosley merged commit 192f26f into SciTools:master Jun 28, 2017

pp-mo deleted the derived_coords_lazycalc branch March 18, 2022 15:41

Calculation of derived coords points and bounds is always lazy. #2604

Calculation of derived coords points and bounds is always lazy. #2604

Uh oh!

Conversation

pp-mo commented Jun 14, 2017

Uh oh!

corinnebosley left a comment

Choose a reason for hiding this comment

Uh oh!

pp-mo commented Jun 14, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pp-mo commented Jun 15, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pp-mo commented Jun 22, 2017

Uh oh!

pp-mo commented Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

bjlittle commented Jun 22, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bjlittle Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bjlittle Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bjlittle Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo commented Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jun 23, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pp-mo Jun 23, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bjlittle Jun 22, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

pp-mo commented Jun 14, 2017 •

edited

Loading

pp-mo commented Jun 15, 2017 •

edited

Loading

pp-mo commented Jun 22, 2017 •

edited

Loading

bjlittle Jun 22, 2017 •

edited

Loading

bjlittle Jun 22, 2017 •

edited

Loading

bjlittle Jun 22, 2017 •

edited

Loading

pp-mo commented Jun 22, 2017 •

edited

Loading

pp-mo Jun 22, 2017 •

edited

Loading

pp-mo Jun 23, 2017 •

edited

Loading

pp-mo Jun 23, 2017 •

edited

Loading

bjlittle Jun 22, 2017 •

edited

Loading

bjlittle Jun 22, 2017 •

edited

Loading

pp-mo Jun 23, 2017 •

edited

Loading

pp-mo commented Jun 23, 2017 •

edited

Loading

pp-mo commented Jun 27, 2017 •

edited

Loading