Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Climate time series #256

Closed
nicolasmerino41 opened this issue Aug 28, 2024 · 30 comments
Closed

Climate time series #256

nicolasmerino41 opened this issue Aug 28, 2024 · 30 comments

Comments

@nicolasmerino41
Copy link
Contributor

nicolasmerino41 commented Aug 28, 2024

Hey there,

Would you mind pointing me out what would be the way of using a climate series in a simulation? I mean, let's say I have a RasterSeries of a climate variable across time. And I want to run a simulation that runs through this climate variation. What'd be the way of implementing this (if feasible)?

Thxs :)
Nico

@rafaqz
Copy link
Member

rafaqz commented Aug 28, 2024

Just pass it as aux=(; climate_aux_name=climatedata, ...), and make sure that it has a Ti dimension with values over the time you simulation runs for tspan. Like they should both be DateTime, or both Int years etc.

Then it should "just work". If you index with get(rule.param, I) and param is Aux{:climate_aux_name}() it will get the correct time at position I.

It does the time lookup each timstep so it's instant within the rule for each cell.

If your time steps don't match the climate data exactly it helps to use sampling=Intervals(Start()) etc in the DD lookup of the To dimensions, so it will find the correct interval with Contains rather than needing the exact time with At. See the DD docs for this stuff.

@nicolasmerino41
Copy link
Contributor Author

Thxs! I'll try it out and see how it goes :)

@rafaqz
Copy link
Member

rafaqz commented Aug 28, 2024

Great. If you ever feel like adding these and other things you learn to the docs somewhere it would be very much appreciated

@nicolasmerino41
Copy link
Contributor Author

Sure! I'll be happy to do so when I figure out how to do it :) I'll add info on MakieOutput and Reflect() for now, and anything else later on.

@rafaqz
Copy link
Member

rafaqz commented Aug 31, 2024

Did you get this working?

@nicolasmerino41
Copy link
Contributor Author

nicolasmerino41 commented Aug 31, 2024 via email

@nicolasmerino41
Copy link
Contributor Author

Ok, so I've been trying to implement this but I ran into another issue that I had been ignoring so far. Since my grids hold complex objects, I've been using Matrix to hold those values, and all the other layers were also Matrix, but since now I'm trying to run a RastersSeries, I need at least one layer to be a RasterSeries, let's say it's the climate. Which leads to the error:
ERROR: ArgumentError: init grid sizes do not match.
I assume there's no way to run a series of matrices, but I believe Rasters cannot hold custom objects either, like:

struct MyStructs256{T <: AbstractFloat} <: FieldVector{2, T}
    a::SVector{256, T}
    b::T

    # Custom constructor for automatic sum calculation
    function MyStructs256(a::SVector{256, T}) where T <: AbstractFloat
        new{T}(a, sum(a))
    end

    # Explicit constructor allowing manual setting of both `a` and `b`
    MyStructs256(a::SVector{256, T}, b::T) where T <: AbstractFloat = new{T}(a, b)
end

Any idea on how to solve this? :)

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Sorry, not sure I get it, I need a lot more context.

But a Raster can hold literally any AbstractArray.

Also don't understand the "run a RasterSeries" part... A RasterSeries is basically just an array of Rasters with some nice helpers. It wont actually work as anything that can "run" in DynamicGrids, I assume you mean as "aux" data? You want to Rasters.combine(series) to get a single Raster (nested arrays of arrays are not so good on GPU so single multidimensional arrays are preferred - but the real reality is I just haven't written DynsmicGrids to get rasters from a series in aux even though it makes sense - it just gets a slice view from a single multidimensional Raster/DimArray. This could be fixed fairly easily).

But again, I just need more code - the only way to be clear in a conversation like this is to post the code

(Oh are you trying to use a RasterSeries as init ? I really dont understand whats happening - init it the template for the array thats used in the simulation - it cant be a nested array. Are you not yet using aux data? that's the way to get static data into a simulation...)

@nicolasmerino41
Copy link
Contributor Author

Ok! I now see what you mean. For some reason I thought Raster could not hold this type of value. I just brought all my grids back to being a Raster.

However, still don't understand how I should write the instructions. Here's the code, given that combined_raster = Rasters.combine(raster_series, Ti)[bounds...].
But I don't know if combined_raster should go only to aux or also to the init; I get error both ways. If it's not in the init I get ERROR: type NamedTuple has no field combined_raster; and if it's in the init I get ERROR: ArgumentError: init grid sizes do not match which I get comes from the time dimension.

Full code:

combined_raster = Rasters.combine(raster_series, Ti)[bounds...]
aux = (; combined_raster=combined_raster)
tspan = Date(2023, 1):Month(1):Date(2023, 8)

init_for_timeseries = (
    state = raster_with_abundances,
    state_richness = raster_richness,
    combined_raster = combined_raster
)

function int_Gr_for_timeseries(state::MyStructs256, self_regulation::AbstractFloat, combined_raster::AbstractFloat)
    return MyStructs256(SVector{256, Float64}(state.a .* self_regulation .* combined_raster))
end

timeseries_rule = Cell{Tuple{:state, :combined_raster}, :state}() do data, (state, combined_raster), I
    merged_state = state + int_Gr_for_timeseries(state, self_regulation, combined_raster)
    return MyStructs256(SVector{256, Float64}(merged_state.a))
end

array_output = ResultOutput(
    init_for_timeseries; tspan = tspan,
    aux = aux,
    mask = raster_sum
)
@time s = sim!(array_output, Ruleset(timeseries_rule))

Sorry my question was not very clear :(

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Ok so seems I don't make these distinctions clear enough 😅

init data is the first frame of the grid in the simulation - its the single time slices of the variables you want to know what happens to during the simulation. They are things that change with time based on your rules.

aux data is variables you already know, either static data that is always the same (e.g. elevation data) or as a time-series (e.g. climate data). It can also be a non-spatial timeseries, e.g. some aggregate economic data in a vector over time. Whatever it is, it is not affected by the simulation. If they are some kind of AbstractDimArray we try to synchronise the time dimension to the simulation time using Contains(current_time) for Intervals and At(current_time) for Points sampling. So in get(data, Aux{:auxname}(), I) we index into the grid location I for the matching time slice, or just take the global value at that time if its just the time dimension.

I'm assuming your combined_raster is a time-series of a variable you already know values for over your time-series.

So to rewrite:

combined_raster = Rasters.combine(raster_series, Ti)[bounds...]
aux = (; combined_raster=combined_raster)
tspan = Date(2023, 1):Month(1):Date(2023, 8)

init_for_timeseries = (
    state = raster_with_abundances,
    state_richness = raster_richness,
)

function int_Gr_for_timeseries(state::MyStructs256, self_regulation::AbstractFloat, combined_raster::AbstractFloat)
    return MyStructs256(SVector{256, Float64}(state.a .* self_regulation .* combined_raster))
end

timeseries_rule = Cell{:state, :state}() do data, state, I
    combined_raster = get(data, Aux(:combined_raster), I) # Important bit here!!
    merged_state = state + int_Gr_for_timeseries(state, self_regulation, combined_raster)
    return MyStructs256(SVector{256, Float64}(merged_state.a))
end

array_output = ResultOutput(
    init_for_timeseries; tspan = tspan,
    aux = aux,
    mask = raster_sum
)
@time s = sim!(array_output, Ruleset(timeseries_rule))

And ultimately to make something reusable you make this a field of your rule struct, and use this get call instead:

combined_raster = get(data, rule.combined_raster, I)

And users can pass Aux(:someaux) or 1.0 or even Grid(:somegrid) for the combined_raster field and use whatever data they want with your rule. get will just handle whatever that is.

@nicolasmerino41
Copy link
Contributor Author

nicolasmerino41 commented Sep 3, 2024

I see! I'll need some time to understand everything you just said but I get the overall functionality :) However, I just ran your code and still throws the ERROR: type NamedTuple has no field combined_raster :(

Now is referring to not finding it in aux, right?

Stacktrace:

getindex(t::@NamedTuple{}, i::Symbol) at [namedtuple.jl](vscode-file://vscode-app/c:/Users/MM-1/AppData/Local/Programs/Microsoft%20VS%20Code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

aux(nt::@NamedTuple{}, ::Aux{:combined_raster}) at [parametersources.jl](vscode-file://vscode-app/c:/Users/MM-1/AppData/Local/Programs/Microsoft%20VS%20Code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

aux at [extent.jl](vscode-file://vscode-app/c:/Users/MM-1/AppData/Local/Programs/Microsoft%20VS%20Code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

...

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

I can't actually test it because its not a complete MWE ;)

But maybe just a typeo somewhere or some variable not updated.

If you make a full MWE I will always test my answers - but I don't have time to generate fake data to test it

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Yeah in this error you have an empty namedtuple for aux:

aux(nt::@NamedTuple{}, ::Aux{:combined_raster})

I'm not sure why

@nicolasmerino41
Copy link
Contributor Author

Me neither jaja Sure, I can make a MWE, it looked like a straight error

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Ohhhh I get it sorry. Yeah this is kind of annoying. You need to pass in the aux as a variable.

On GPUs we need to thin the data a lot or the kernel size overflows and wont compile. So we only pass in the things that we need for each specific rule, and we remove all grids and aux data that isn't used and send the simplest possible thing to the GPU.

By putting the aux in as a variable we can find it, but we don't actually scan the code for manually defined Aux or Grid objects. Probably we should do that too its just a pain.

But this should fix it:

timeseries_rule = let combined_raster_aux=Aux{:combined_raster}()
    Cell{:state, :state}() do data, state, I
        combined_raster = get(data, combined_raster_aux, I) # Important bit here!!
        merged_state = state + int_Gr_for_timeseries(state, self_regulation, combined_raster)
        return MyStructs256(SVector{256, Float64}(merged_state.a))
    end
end

This is kind of awful actually, I didn't think about this so much when I wrote that thinning. Thats why this is in dev still...

(there is some magic happening here in that any variable in the let block will be a field in the anonymous function struct made by the do block. And we can find that struct field with Flatten.jl. passing variables like this also lets MakieOutput make sliders for Param variables using the same trick. But it seems like complet magic that its even different)

@nicolasmerino41
Copy link
Contributor Author

Jajaja, ok, if I was lost before, I definitely am now. But I see what you mean. Now the error switched to ERROR: type NamedTuple has no field state

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Yeah, need that MWE or there are just typos and little bugs

@nicolasmerino41
Copy link
Contributor Author

nicolasmerino41 commented Sep 3, 2024

By: Probably we should do that too its just a pain, you mean we should fix it?
MWE on the way

@nicolasmerino41
Copy link
Contributor Author

using StaticArrays, Serialization, DynamicGrids, Dispersal, Rasters, Dates

struct MyStructs256{T <: AbstractFloat} <: FieldVector{2, T}
    a::SVector{256, T}
    b::T

    # Custom constructor for automatic sum calculation
    function MyStructs256(a::SVector{256, T}) where T <: AbstractFloat
        new{T}(a, sum(a))
    end

    # Explicit constructor allowing manual setting of both `a` and `b`
    MyStructs256(a::SVector{256, T}, b::T) where T <: AbstractFloat = new{T}(a, b)
end

# Define zero and oneunit for MyStructs256
Base.zero(::Type{MyStructs256{T}}) where {T <: AbstractFloat} = MyStructs256(SVector{256, T}(fill(zero(T), 256)), zero(T))
Base.zero(x::MyStructs256{T}) where {T <: AbstractFloat} = MyStructs256(SVector{256, T}(fill(zero(T), 256)), zero(T))
Base.oneunit(::Type{MyStructs256{T}}) where {T <: AbstractFloat} = MyStructs256(fill(oneunit(T), 256), oneunit(T))

# Comparison based on 'b' field
Base.isless(x::MyStructs256, y::MyStructs256) = isless(x.b, y.b)
Base.isless(x::MyStructs256, y::AbstractFloat) = isless(x.b, y)

# Element-wise arithmetic operations ensuring 'b' is recalculated correctly
Base.:+(x::MyStructs256, y::MyStructs256) = MyStructs256(x.a .+ y.a, sum(x.a .+ y.a))
Base.:-(x::MyStructs256, y::MyStructs256) = MyStructs256(x.a .- y.a, sum(x.a .- y.a))
Base.:*(x::MyStructs256, scalar::Real) = MyStructs256(x.a .* scalar, sum(x.a .* scalar))
Base.:/(x::MyStructs256, scalar::Real) = MyStructs256(x.a ./ scalar, sum(x.a ./ scalar))
Base.:-(x::MyStructs256, scalar::Real) = MyStructs256(x.a .- scalar, x.b - scalar * 256)
Base.:+(x::MyStructs256, scalar::Real) = MyStructs256(x.a .+ scalar, x.b + scalar * 256)

# Define what a NaN is for MyStructs256
Base.isnan(x::MyStructs256) = isnan(x.b) || any(isnan, x.a)

# Adding a method in the sum function for MyStructs256
function Base.sum(structs::MyStructs256...)
    # Sum the 'a' vectors
    summed_a = sum([s.a for s in structs])

    # Sum the 'b' values
    summed_b = sum([s.b for s in structs])

    # Create a new MyStructs256 instance with the summed results
    return MyStructs256(summed_a, summed_b)
end

# Adding a method to maximum
# Define maximum for MyStructs256
function Base.maximum(a::MyStructs256, b::MyStructs256)
    return MyStructs256(max.(a.a, b.a))
end

# Define maximum for MyStructs256 with a scalar
function Base.maximum(a::MyStructs256, b::AbstractFloat)
    return MyStructs256(max.(a.a, b))
end

# Define maximum for a scalar with MyStructs256
function Base.maximum(a::AbstractFloat, b::MyStructs256)
    return MyStructs256(max.(a, b.a))
end

# Define maximum for MyStructs256
function Base.maximum(a::MyStructs256)
    return maximum(a.a)
end

# Define maximum for a matrix of MyStructs256
function Base.maximum(a::Matrix{MyStructs256{AbstractFloat}})
    # Extract all `b` values from each MyStructs256 element in the matrix and find the maximum
    return maximum(map(x -> x.b, a))
end

########### HERE STARTS THE RELEVANT PART ##########
combined_raster = deserialize("combined_raster.jls")
raster_with_abundances = deserialize("raster_with_abundances.jls")
raster_sum = deserialize("Objects\\raster_sum.jls")
aux = (; combined_raster=combined_raster)
tspan = Date(2023, 1):Month(1):Date(2023, 8)

state = raster_with_abundances

self_regulation = 0.01
function int_Gr_for_timeseries(state::MyStructs256, self_regulation::AbstractFloat, combined_raster::AbstractFloat)
    return MyStructs256(SVector{256, Float64}(state.a .* self_regulation .* combined_raster))
end

timeseries_rule = let combined_raster_aux=Aux{:combined_raster}()
    Cell{:state, :state}() do data, state, I
        combined_raster = get(data, combined_raster_aux, I) # Important bit here!!
        merged_state = state + int_Gr_for_timeseries(state, self_regulation, combined_raster)
        return MyStructs256(SVector{256, Float64}(merged_state.a))
    end
end

array_output = ResultOutput(
    state; tspan = tspan,
    aux = aux,
    mask = raster_sum
)
@time s = sim!(array_output, Ruleset(timeseries_rule))

@nicolasmerino41
Copy link
Contributor Author

Can I send you the objects by email? github doesn't support .jls

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

I mean do I just separate out the GPU code to thin Aux like that and not bother for CPU so we don't have this error...

Or do I analise the function code to find if there are any Aux or Grid object specified in the code and make sure to include them. This is just a pretty ugly thing to do, but maybe we need it.

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Can I send you the objects by email? github doesn't support .jls

Just make the data with rand inline in the code ;)

(Always best everyone else reading this can see it too)

@nicolasmerino41
Copy link
Contributor Author

nicolasmerino41 commented Sep 3, 2024

Ok, here it is. It works fine and gives the exact same error:

######## TRYING TIME SERIES #########
using StaticArrays, Serialization, DynamicGrids, Dispersal, Rasters, Dates

struct MyStructs256{T <: AbstractFloat} <: FieldVector{2, T}
    a::SVector{256, T}
    b::T

    # Custom constructor for automatic sum calculation
    function MyStructs256(a::SVector{256, T}) where T <: AbstractFloat
        new{T}(a, sum(a))
    end

    # Explicit constructor allowing manual setting of both `a` and `b`
    MyStructs256(a::SVector{256, T}, b::T) where T <: AbstractFloat = new{T}(a, b)
end

# Define zero and oneunit for MyStructs256
Base.zero(::Type{MyStructs256{T}}) where {T <: AbstractFloat} = MyStructs256(SVector{256, T}(fill(zero(T), 256)), zero(T))
Base.zero(x::MyStructs256{T}) where {T <: AbstractFloat} = MyStructs256(SVector{256, T}(fill(zero(T), 256)), zero(T))
Base.oneunit(::Type{MyStructs256{T}}) where {T <: AbstractFloat} = MyStructs256(fill(oneunit(T), 256), oneunit(T))

# Comparison based on 'b' field
Base.isless(x::MyStructs256, y::MyStructs256) = isless(x.b, y.b)
Base.isless(x::MyStructs256, y::AbstractFloat) = isless(x.b, y)

# Element-wise arithmetic operations ensuring 'b' is recalculated correctly
Base.:+(x::MyStructs256, y::MyStructs256) = MyStructs256(x.a .+ y.a, sum(x.a .+ y.a))
Base.:-(x::MyStructs256, y::MyStructs256) = MyStructs256(x.a .- y.a, sum(x.a .- y.a))
Base.:*(x::MyStructs256, scalar::Real) = MyStructs256(x.a .* scalar, sum(x.a .* scalar))
Base.:/(x::MyStructs256, scalar::Real) = MyStructs256(x.a ./ scalar, sum(x.a ./ scalar))
Base.:-(x::MyStructs256, scalar::Real) = MyStructs256(x.a .- scalar, x.b - scalar * 256)
Base.:+(x::MyStructs256, scalar::Real) = MyStructs256(x.a .+ scalar, x.b + scalar * 256)

# Define what a NaN is for MyStructs256
Base.isnan(x::MyStructs256) = isnan(x.b) || any(isnan, x.a)

# Adding a method in the sum function for MyStructs256
function Base.sum(structs::MyStructs256...)
    # Sum the 'a' vectors
    summed_a = sum([s.a for s in structs])

    # Sum the 'b' values
    summed_b = sum([s.b for s in structs])

    # Create a new MyStructs256 instance with the summed results
    return MyStructs256(summed_a, summed_b)
end

# Adding a method to maximum
# Define maximum for MyStructs256
function Base.maximum(a::MyStructs256, b::MyStructs256)
    return MyStructs256(max.(a.a, b.a))
end

# Define maximum for MyStructs256 with a scalar
function Base.maximum(a::MyStructs256, b::AbstractFloat)
    return MyStructs256(max.(a.a, b))
end

# Define maximum for a scalar with MyStructs256
function Base.maximum(a::AbstractFloat, b::MyStructs256)
    return MyStructs256(max.(a, b.a))
end

# Define maximum for MyStructs256
function Base.maximum(a::MyStructs256)
    return maximum(a.a)
end

# Define maximum for a matrix of MyStructs256
function Base.maximum(a::Matrix{MyStructs256{AbstractFloat}})
    # Extract all `b` values from each MyStructs256 element in the matrix and find the maximum
    return maximum(map(x -> x.b, a))
end

########### CREATING DATA ##########
######### raster_with_abundances ########
# Define the size of the raster
size_x, size_y = 125, 76
# Define the x and y dimensions directly
x_dim = X(1:size_x)
y_dim = Y(1:size_y)
# A 125x76 grid of random MyStructs256
raster_matrix = reshape([MyStructs256(SVector{256, Float64}(rand(256))) for _ in 1:(size_x * size_y)], size_x, size_y)
raster_with_abundances = Raster(raster_matrix, dims=(x_dim, y_dim))
######### raster_sum ########
# Create a 125x76 grid of Bool values with first and last rows as false, others true
sum_matrix = [i == 1 || i == size_x ? false : true for i in 1:size_x, j in 1:size_y]
# Create the Bool raster with the specified dimensions
raster_sum = Raster(sum_matrix, dims=(x_dim, y_dim))
######## combined_raster ########
time_dim = Ti([Date(2023, i, 1) for i in 1:8])  # 8 time points
# Create 8 individual rasters with random Int16 values
raster_layers = [
    Raster(rand(Int16, size_x, size_y), dims=(x_dim, y_dim)) for _ in 1:8
]
# Set the first and last rows to -32768 to represent missing values
for raster in raster_layers
    raster[1, :] .= -32768  # First row
    raster[end, :] .= -32768  # Last row
end
# Combine the rasters into a 3D raster
combined_raster = Raster(cat(raster_layers..., dims=3), dims=(x_dim, y_dim, time_dim))

######### HERE STARTS THE RELEVANT PART #############
aux = (; combined_raster=combined_raster)
tspan = Date(2023, 1):Month(1):Date(2023, 8)

init_for_timeseries = (
    state = raster_with_abundances
)

self_regulation = 0.01
function int_Gr_for_timeseries(state::MyStructs256, self_regulation::AbstractFloat, combined_raster::AbstractFloat)
    return MyStructs256(SVector{256, Float64}(state.a .* self_regulation .* combined_raster))
end

timeseries_rule = let combined_raster_aux=Aux{:combined_raster}()
    Cell{:state, :state}() do data, state, I
        combined_raster = get(data, combined_raster_aux, I) # Important bit here!!
        merged_state = state + int_Gr_for_timeseries(state, self_regulation, combined_raster)
        return MyStructs256(SVector{256, Float64}(merged_state.a))
    end
end

array_output = ResultOutput(
    init_for_timeseries; tspan = tspan,
    aux = aux,
    mask = raster_sum
)
@time s = sim!(array_output, Ruleset(timeseries_rule))

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

Ok thats an easy one! By deleting the second argument this is no longer a namedtuple, its just an array:

init_for_timeseries = (
    state = raster_with_abundances
)

And of course you can just pass a single array as init for the unnamed basic case of rule.

Personally I always use ; at the start of a NamedTuple definition so this doesn't happen when i delete the second + arguments:

init_for_timeseries = (; # this makes all the difference 
    state = raster_with_abundances
)

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

And finally for the simulation to run you need this:

function int_Gr_for_timeseries(state::MyStructs256, self_regulation::AbstractFloat, combined_raster)
    return MyStructs256(SVector{256, Float64}(state.a .* self_regulation .* combined_raster))
end

Because combined_raster is not longer AbstractFloat.

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

And maybe you notice I made your code block colorful... it helps to use ```julia as the start fence so github knows to use julia colors. Much easier to read that way

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

We could have better errors like "Youre rule $ruletype wants a grid called :state in init but there are no named grids" or
"Youre rule wants a grid called :state in init but there are only :a, :b, and :c grids"

And the same for aux. Currently they are just the Base julia NamedTuple errors and not very helpful to understand the problem

@nicolasmerino41
Copy link
Contributor Author

Ok! Finally! It worked. Thank you very much Raf, really appreciate your patience :)

@nicolasmerino41
Copy link
Contributor Author

Once I understand the structure behind all you told me, I'm happy to PR #260 if you don't have time. Or perhaps it's very straight forward for you and it'll only take you a sec :)

I'm diving into the GPU part of DG these days so it'll be useful for me to understand it anyway.

@rafaqz
Copy link
Member

rafaqz commented Sep 3, 2024

So in the RuleData constructor we switch from the whole SimData object to just the object needed for a specific rule.

function RuleData(d::AbstractSimData, rule::Rule)
aux_in_rule = Flatten.flatten(rule, Aux)
rule_aux_keys = map(aux_in_rule) do aux
_unwrap(aux)
end
if isnothing(aux(d))
if length(rule_aux_keys) > 0
throw(ArgumentError("no aux data available: expected aux with keys $rule_aux_keys"))
end
return RuleData(d)
end
rule_aux = aux(d)[rule_aux_keys]
# Update the extent to have less aux
rule_extent = ConstructionBase.setproperties(extent(d), (; aux=rule_aux))
# Use the thinned rule extent in RuleData
return RuleData(d; extent=rule_extent)
end

We could conditionally just use all the aux fields based on thin. But we need to pass the thin parameter to a field in SimData as Val{true}() or Val{false} or it wont be type-stable.

Then, making Grid thinning also depend on thin is a bunch more work in maprules. Probably I need to do this lol.

It needs cleaning up anyway its kind of a mess.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants