Flip location along accumulating dimension for `Accumulating` scans

[glwagner]

Addresses #4354

Things to think about:

• does this apply to all accumulations?
• should we throw an error when you try to accumulate along an already-reduced dimension (eg the location is Nothing)
• can we test this? maybe a simple unit test is enough but it would also be nice to know if this is "correct" somehow
• other concerns I didn't think of

cc @tomchor

[tomchor]

Addresses #4354

Things to think about:

* does this apply to all accumulations?

I think so! Can you point me to the accumulations we currently support? Should be hard to consider them one by one.

* should we throw an error when you try to accumulate along an already-reduced dimension (eg the location is `Nothing`)

Not sure. Maybe at least a warning?

* can we test this? maybe a simple unit test is enough but it would also be nice to know if this is "correct" somehow

I'd say a unit test at least. It also shouldn't be hard to cook up a quantitative example with regular and stretched grids. Something along the lines of

julia> cumsum(a);

julia> cumsum(a) == [i for i in 1:4]
true

?

I can help you out with this if you want.

* other concerns I didn't think of

Possibly testing subsequent accumulations over different dimensions? i.e. accumulating over dim=1, then 2, then 3, should be the same as accumulating over 1, 2, 3.

[glwagner]

I think so! Can you point me to the accumulations we currently support?

Please take a look at the code modified by this PR, it's all in that file.

I can help you out with this if you want.

I'd greatly appreciate help!

Possibly testing subsequent accumulations over different dimensions? i.e. accumulating over dim=1, then 2, then 3, should be the same as accumulating over 1, 2, 3.

Maybe this can be easily tested with an array full of 1s.

[glwagner]

Tests fail because this PR actually changes the meaning of the accumulations. Now we are tasked with defining the first element of the accumulation. For a cumsum this is 0, but other accumulating reductions (which may be invoked rarely, who knows) require other choices, which I believe are equivalent to the neutral element / mask that is used.

[tomchor]

Tests fail because this PR actually changes the meaning of the accumulations. Now we are tasked with defining the first element of the accumulation. For a cumsum this is 0, but other accumulating reductions (which may be invoked rarely, who knows) require other choices, which I believe are equivalent to the neutral element / mask that is used.

To be clear, you're talking about adding more of these methods?:

neutral_element(::typeof(Base.add_sum), T) = convert(T, 0)

[tomchor]

I'm a bit confused about some of the behavior though. Taking one of the tests as an example, you define T as a 2x2x2 CenterField. in the first y and z indices we have:

infil> interior(T, :, 1, 1)
2-element view(::Array{Float64, 3}, 3:4, 3, 3) with eltype Float64:
 1.5
 2.5

And then you define an x-accumulation of T called Tcx, which in the same indices returns:

2-element view(::Array{Float64, 3}, 3:4, 3, 3) with eltype Float64:
 0.0
 2.5

This doesn't really make sense to me. I expected the two elements to be [1.5, 4], which is actually what the current test is testing for. Are you envisioning that in this case the neutral element should be the first element of the array? @glwagner can you clarify?

[glwagner]

I'm a bit confused about some of the behavior though. Taking one of the tests as an example, you define T as a 2x2x2 CenterField. in the first y and z indices we have:

infil> interior(T, :, 1, 1)
2-element view(::Array{Float64, 3}, 3:4, 3, 3) with eltype Float64:
 1.5
 2.5

And then you define an x-accumulation of T called Tcx, which in the same indices returns:

2-element view(::Array{Float64, 3}, 3:4, 3, 3) with eltype Float64:
 0.0
 2.5

This doesn't really make sense to me. I expected the two elements to be [1.5, 4], which is actually what the current test is testing for. Are you envisioning that in this case the neutral element should be the first element of the array? @glwagner can you clarify?

Let's think about this. I think it depends on the location of the parent field.

If we integrate a Center field and we store the result at interfaces, then the first value represents the integral computed over no points, so its 0.

Therefore I expect

T = [1.5, 2.5]
cumint_T = [0, 1.5, 4.0]

This seems right to me and also has the property that the derivative returns you the original array.

If we integrate a Face field, I think that index 1 stores the first value, not the neutral value. For example

cum_int_T = [0, 1.5, 4] # at Face
double_int_T = [1.5, 5.5] # at Center

That's my thought. Does this make sense?

[tomchor]

That makes perfect sense, and it was my thought as well. However, the current behavior is that when cumulatively integrating over a periodic dimension, faces and center grids have the same size. I believe we also need to change this behavior then for this to make sense, no?

[glwagner]

That makes perfect sense, and it was my thought as well. However, the current behavior is that when cumulatively integrating over a periodic dimension, faces and center grids have the same size. I believe we also need to change this behavior then for this to make sense, no?

Are you sure? Changing the location should change the size.

[tomchor]

That makes perfect sense, and it was my thought as well. However, the current behavior is that when cumulatively integrating over a periodic dimension, faces and center grids have the same size. I believe we also need to change this behavior then for this to make sense, no?

Are you sure? Changing the location should change the size.

Not on periodic dimensions:

julia> grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1));

julia> T = CenterField(grid)
4×4×4 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 10×10×10 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, -2:7) with eltype Float64 with indices -2:7×-2:7×-2:7
    └── max=0.0, min=0.0, mean=0.0

julia> Field(CumulativeIntegral(T, dims=1))
4×4×4 Field{Face, Center, Center} on RectilinearGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 4)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── operand: CumulativeIntegral of BinaryOperation at (Center, Center, Center) over dims 1
├── status: time=0.0
└── data: 10×10×10 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, -2:7) with eltype Float64 with indices -2:7×-2:7×-2:7
    └── max=0.0, min=0.0, mean=0.0

It does work for nonperiodic dimensions though:

julia> Field(CumulativeIntegral(T, dims=3))
4×4×5 Field{Center, Center, Face} on RectilinearGrid on CPU
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 5)
├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── operand: CumulativeIntegral of BinaryOperation at (Center, Center, Center) over dims 3
├── status: time=0.0
└── data: 10×10×11 OffsetArray(::Array{Float64, 3}, -2:7, -2:7, -2:8) with eltype Float64 with indices -2:7×-2:7×-2:8
    └── max=0.0, min=0.0, mean=0.0

[glwagner]

Center and Face have the same size on Periodic dimensions.

Also, I am confused whether this makes sense in a Periodic direction. As an example, f(x) = 1 is Periodic but its cumulative integral is not.

[tomchor]

Center and Face have the same size on Periodic dimensions.

Also, I am confused whether this makes sense in a Periodic direction. As an example, f(x) = 1 is Periodic but its cumulative integral is not.

Yes, I guess I haven't been expressing myself very clearly, but that has been precisely my point in these last few comments. The way we treat periodic directions doesn't work well with integrals.

[glwagner]

Center and Face have the same size on Periodic dimensions. Also, I am confused whether this makes sense in a Periodic direction. As an example, f(x) = 1 is Periodic but its cumulative integral is not.

Yes, I guess I haven't been expressing myself very clearly, but that has been precisely my point in these last few comments. The way we treat periodic directions doesn't work well with integrals.

Okay, but you are implying we are doing something wrong. My last question asks: is CumulativeIntegral invalid in Periodic directions? For example, we can throw an error for this case. I'm not sure any accumulation is valid.

[tomchor]

Center and Face have the same size on Periodic dimensions. Also, I am confused whether this makes sense in a Periodic direction. As an example, f(x) = 1 is Periodic but its cumulative integral is not.

Yes, I guess I haven't been expressing myself very clearly, but that has been precisely my point in these last few comments. The way we treat periodic directions doesn't work well with integrals.

Okay, but you are implying we are doing something wrong. My last question asks: is CumulativeIntegral invalid in Periodic directions? For example, we can throw an error for this case. I'm not sure any accumulation is valid.

afaik there's no problem in defining a CumulativeIntegral for a periodic dimension. At least in principle. Like you exemplified yourself, $f(x)=1$ is periodic, but it's integral isn't. If the math works out, I don't see why we can't make the numerics work out too.

I think the trouble for us is that Fields are a bit stiff when it comes to sizes. For example, given a grid with a particular topology, the array sizes for all Fields are predetermined given their location. However, for cumulative integrals, regardless of the topology or cell location, the CumulativeIntegral of a Field should have a size that's larger (by one) in the integrated dimension. I think that's the only way to ensure consistency, and to ensure that the derivative of a CumulativeField equals the original Field. I think you illustrated that nicely in https://github.com/CliMA/Oceananigans.jl/pull/4358#issuecomment-2817115385.

So maybe we need to create a different abstraction to deal with cumulative integrals that is sized differently? Not sure...

Hopefully things are a bit clearer now.

[glwagner]

The problem is that integral of f(x) = 1 is not periodic so it cannot be represented with a periodic Field.

[tomchor]

The problem is that integral of f(x) = 1 is not periodic so it cannot be represented with a periodic Field.

I'm aware. What I'm advocating for is that integrals of periodic Fields shouldn't be periodic Fields. In other words, an integral of a Field shouldn't inherit its original topology.

[glwagner]

The problem is that integral of f(x) = 1 is not periodic so it cannot be represented with a periodic Field.

I'm aware. What I'm advocating for is that integrals of periodic Fields shouldn't be periodic Fields. In other words, an integral of a Field shouldn't inherit its original topology.

What topology would you use? How would you fill the halos of this new topology? What system of rules will be used to change between "related topologies" upon differentiation and integration? For example, the derivative of a Bounded field is not Periodic. Are you proposing to create a NonPeriodicIntegrated topology or something like that? What is a concrete example of some valid mathematical / physical object that this could be used for?

[tomchor]

I don't know how to answer all those questions. I'm not that fast with software design. My only point is that it's definitely possible numerically.

That said, given your questions, it's probably be beyond the scope of this PR. Should we just throw an error when integrating over Periodic dimensions for now?

[glwagner]

I don't know how to answer all those questions. I'm not that fast with software design. My only point is that it's definitely possible numerically.

That said, given your questions, it's probably be beyond the scope of this PR. Should we just throw an error when integrating over Periodic dimensions for now?

Yes. I am also skeptical that one should pursue an abstraction for "integrals of Periodic functions that are not themselves Periodic".

[jagoosw]

Sorry for jumping into this, I might have missed something.

If we integrate a Face field, I think that index 1 stores the first value, not the neutral value. For example

cum_int_T = [0, 1.5, 4] # at Face
double_int_T = [1.5, 5.5] # at Center

That's my thought. Does this make sense?

I think that this isn't consistent with the first/center integral version and doesn't return the origional array when differentiated, and that to be correct the neutral value should be in the first center position. Where solid lines are faces and dashed lines are centers

[glwagner]

I think you're right, so we need

cum_int_T = [0, 1.5, 4] # at Face
double_int_T = [0, 1.5] # at Center

right? And to satisfy the first boundary value, we'd in principle need 5.5 on the first halo point