text:
using DocStringExtensions
"""
Computes the fluxes and changes in storages in a hillslope hydrological response unit of the model taking into account snow redistribution.

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function hillslopeHRU_with_snow_redistribution(hill::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    @assert hill.Total_Interception_Evaporation == 0
    @assert hill.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(hill.Area_Elevations) <= 1.000001 || sum(hill.Area_Elevations) == 0
    @assert hill.Area_HRU >= 0 and <= 1
    @assert hill.Nr_Elevationbands >= 1
    # define Arrays for Snow and Interception storage
    Snow = zeros(hill.Nr_Elevationbands)
    Interception = zeros(hill.Nr_Elevationbands)
    Snow_Cover = zeros(hill.Nr_Elevationbands)
    Total_Melt_Glacier = 0.0
    Snow_Old = copy(storages.Snow)
    # calculate the model components that are height dependend
    #do so in reversed order to account for snow redistribution
    for i in reverse(1 : hill.Nr_Elevationbands)
        # snow component
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(hill.Area_Glacier[i], hill.Precipitation[i], hill.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(hill.Potential_Evaporation_Mean, hill.Precipitation[i], hill.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        hill.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * hill.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        hill.Total_Interception_Evaporation::Float64 += (Interception_Evaporation * hill.Area_Elevations[i])
        Total_Melt_Glacier::Float64 += (Melt_Glacier * hill.Area_Elevations[i])
        # get the extend of snow cover, if there is snow in the snow bucket the snow cover is 1
        if Snow[i] > 1
            Snow_Cover[i] = 1
        else
            Snow_Cover[i] = hill.Area_Glacier[i]
        end
        # snow redistribution
        # if elevation band 2700m or higher
        Current_Elevation = (hill.Elevation_Count[i]-1)* 200 + hill.Catchment_Elevation[1]
        if Current_Elevation >= hill.Snow_Redistribution[2] && Snow[i] > hill.Snow_Redistribution[1]
            # gives the index of the elevation count of 2100 m elevation band
            index_2100 = findall(x->x == (2100 - hill.Catchment_Elevation[1]) / 200 , hill.Elevation_Count)[1]
            # snow to be redistributed according to areal extend of elevation zone
            Snow_redistributed = (Snow[i] - hill.Snow_Redistribution[1]) * hill.Area_Elevations[i]
            # snow that has to be added to lower elevation zones in mm
            Snow_redistributed_mm = Snow_redistributed / sum(hill.Area_Elevations[index_2100: i-1])
            # add the snow to the lower elevation zones
            storages.Snow[index_2100: i-1] = storages.Snow[index_2100: i-1] .+ Snow_redistributed_mm
            # substracte the Snow from the higher elevation zone
            Snow[i] = hill.Snow_Redistribution[1]
            #println("redistribution ", Snow_redistributed)
            #println("sum area below ", sum(hill.Area_Elevations[index_2100: i-1]))
        end

    end
    #soil storage component
    Overlandflow::Float64, Preferentialflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = soilstorage(hill.Total_Effective_Precipitation, hill.Total_Interception_Evaporation, hill.Potential_Evaporation_Mean, storages.Soil, parameters.beta, parameters.Ce, parameters.Ratio_Pref, parameters.Soilstoragecapacity)
    GWflow::Float64 = Preferentialflow
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)

    #calculate total outflow
    Flows = Fast_Discharge + GWflow + hill.Total_Interception_Evaporation + Soil_Evaporation
    @assert hill.Potential_Evaporation_Mean - Soil_Evaporation - hill.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges and evaporation fluxed according to areal fraction
    Fast_Discharge = Fast_Discharge * hill.Area_HRU
    GWflow = GWflow * hill.Area_HRU
    Total_Interception_Evaporation = hill.Total_Interception_Evaporation * hill.Area_HRU
    Soil_Evaporation = Soil_Evaporation * hill.Area_HRU
    # returning all fluxes (evporative, discharge)
    hill_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    #returning all storages
    hill_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)
    # total interception should be zero again for next run
    hill.Total_Interception_Evaporation = 0
    hill.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert hill_out.Fast_Discharge >= 0
    @assert hill_out.GWflow >= 0
    @assert hill_out.Soil_Evaporation >= 0
    @assert hill_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert hill_storages.Fast >= 0
    @assert hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Snow >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Soil >= 0
    @assert ones(hill.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert parameters.Soilstoragecapacity - hill_storages.Soil >= 0

    #assertion water balance
    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    # get total storage of interception and snow over all elevation zones
    for i in 1: hill.Nr_Elevationbands
        Precipitation += hill.Precipitation[i] * hill.Area_Elevations[i]
        Interception_Storage_New += hill_storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_New += hill_storages.Snow[i] * hill.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * hill.Area_Elevations[i]
        #Snow_Storage_Old += storages.Snow[i] * hill.Area_Elevations[i]
        Snow_Storage_Old += Snow_Old[i] * hill.Area_Elevations[i]
    end
    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (hill_storages.Fast - storages.Fast) + (hill_storages.Soil - storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)
    # println("Area ELevations ", hill.Area_Elevations)
    # println("difference snow storage ", Snow_Storage_New -Snow_Storage_Old)
    # println("difference snow storage mm ", hill_storages.Snow - Snow_Old)
    # println("old ", Snow_Old)
    # println("new ", hill_storages.Snow)
    #assert that water balance closes
    #println("check", Precipitation + Total_Melt_Glacier - (Flows + All_Storages))
    @assert -0.00000001 <= Precipitation + Total_Melt_Glacier - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation + Total_Melt_Glacier)  * hill.Area_HRU - (Flows_Area + All_Storages * hill.Area_HRU) <= 0.00000001
    return hill_out::Outflows, hill_storages::Storages, Precipitation, All_Storages

end

"""
Computes the fluxes and changes in storages in a hillslope hydrological response unit of the model.

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function hillslopeHRU(hill::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    @assert hill.Total_Interception_Evaporation == 0
    @assert hill.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(hill.Area_Elevations) <= 1.000001 || sum(hill.Area_Elevations) == 0
    @assert hill.Area_HRU >= 0 and <= 1
    @assert hill.Nr_Elevationbands >= 1
    # define Arrays for Snow and Interception storage
    Snow = zeros(hill.Nr_Elevationbands)
    Interception = zeros(hill.Nr_Elevationbands)
    Snow_Cover = zeros(hill.Nr_Elevationbands)
    Total_Melt_Glacier = 0.0
    # calculate the model components that are height dependend
    for i in 1 : hill.Nr_Elevationbands
        # snow component
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(hill.Area_Glacier[i], hill.Precipitation[i], hill.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(hill.Potential_Evaporation_Mean, hill.Precipitation[i], hill.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        hill.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * hill.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        hill.Total_Interception_Evaporation::Float64 += (Interception_Evaporation * hill.Area_Elevations[i])
        Total_Melt_Glacier::Float64 += (Melt_Glacier * hill.Area_Elevations[i])
        # get the extend of snow cover, if there is snow in the snow bucket the snow cover is 1
        if Snow[i] > 1
            Snow_Cover[i] = 1
        else
            Snow_Cover[i] = hill.Area_Glacier[i]
        end
    end
    #soil storage component
    Overlandflow::Float64, Preferentialflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = soilstorage(hill.Total_Effective_Precipitation, hill.Total_Interception_Evaporation, hill.Potential_Evaporation_Mean, storages.Soil, parameters.beta, parameters.Ce, parameters.Ratio_Pref, parameters.Soilstoragecapacity)
    GWflow::Float64 = Preferentialflow
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)

    #calculate total outflow
    Flows = Fast_Discharge + GWflow + hill.Total_Interception_Evaporation + Soil_Evaporation
    @assert hill.Potential_Evaporation_Mean - Soil_Evaporation - hill.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges and evaporation fluxed according to areal fraction
    Fast_Discharge = Fast_Discharge * hill.Area_HRU
    GWflow = GWflow * hill.Area_HRU
    Total_Interception_Evaporation = hill.Total_Interception_Evaporation * hill.Area_HRU
    Soil_Evaporation = Soil_Evaporation * hill.Area_HRU
    # returning all fluxes (evporative, discharge)
    hill_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    #returning all storages
    hill_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)
    # total interception should be zero again for next run
    hill.Total_Interception_Evaporation = 0
    hill.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert hill_out.Fast_Discharge >= 0
    @assert hill_out.GWflow >= 0
    @assert hill_out.Soil_Evaporation >= 0
    @assert hill_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert hill_storages.Fast >= 0
    @assert hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Snow >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Soil >= 0
    @assert ones(hill.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert parameters.Soilstoragecapacity - hill_storages.Soil >= 0

    #assertion water balance
    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    # get total storage of interception and snow over all elevation zones
    for i in 1: hill.Nr_Elevationbands
        Precipitation += hill.Precipitation[i] * hill.Area_Elevations[i]
        Interception_Storage_New += hill_storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_New += hill_storages.Snow[i] * hill.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_Old += storages.Snow[i] * hill.Area_Elevations[i]
    end
    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (hill_storages.Fast - storages.Fast) + (hill_storages.Soil - storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)

    #assert that water balance closes
    @assert -0.00000001 <= Precipitation + Total_Melt_Glacier - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation + Total_Melt_Glacier)  * hill.Area_HRU - (Flows_Area + All_Storages * hill.Area_HRU) <= 0.00000001
    return hill_out::Outflows, hill_storages::Storages, Precipitation, All_Storages

end

"""
Computes the fluxes and changes in storages in a riparian hydrological response unit of the model taking into account snow redistribution.

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function riparianHRU_with_snow_redistribution(rip::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    # riparian HRU has no preferential flow

    @assert rip.Total_Interception_Evaporation == 0
    @assert rip.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(rip.Area_Elevations) <= 1.000001 || sum(rip.Area_Elevations) == 0
    @assert rip.Area_HRU >= 0 and <= 1
    @assert rip.Nr_Elevationbands >= 1
    Snow = zeros(rip.Nr_Elevationbands)
    Interception = zeros(rip.Nr_Elevationbands)
    Snow_Cover = zeros(rip.Nr_Elevationbands)
    Snow_Old = copy(storages.Snow)
    # println("snow first ", Snow_Old)
    # println("temperature ", rip.Temp_Elevation)
    # println("theshold temp ",  parameters.Temp_Thresh)
    for i in reverse(1 : rip.Nr_Elevationbands)
        # snow component
        @assert rip.Area_Glacier[i] == 0.0
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(rip.Area_Glacier[i], rip.Precipitation[i], rip.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        @assert Melt_Glacier == 0.0
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(rip.Potential_Evaporation_Mean, rip.Precipitation[i], rip.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        rip.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * rip.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        rip.Total_Interception_Evaporation::Float64 += Interception_Evaporation * rip.Area_Elevations[i]
        #get snow cover extent
        if Snow[i] > 1
            Snow_Cover[i] = 1
        end
        # snow redistribution
        # if elevation band 2700m or higher
        Current_Elevation = (rip.Elevation_Count[i]-1)* 200 + rip.Catchment_Elevation[1]
        if Current_Elevation >= rip.Snow_Redistribution[2] && Snow[i] > rip.Snow_Redistribution[1]
            # gives the index of the elevation count of 2100 m elevation band
            index_2100 = findall(x->x == (2100 - rip.Catchment_Elevation[1]) / 200 , rip.Elevation_Count)[1]
            # snow to be redistributed according to areal extend of elevation zone
            Snow_redistributed = (Snow[i] - rip.Snow_Redistribution[1]) * rip.Area_Elevations[i]
            # snow that has to be added to lower elevation zones in mm
            Snow_redistributed_mm = Snow_redistributed / sum(rip.Area_Elevations[index_2100: i-1])
            # add the snow to the lower elevation zones
            storages.Snow[index_2100: i-1] = storages.Snow[index_2100: i-1] .+ Snow_redistributed_mm
            # substracte the Snow from the higher elevation zone
            Snow[i] = rip.Snow_Redistribution[1]
            # println("redistribution ", Snow_redistributed)
            # println("sum area below ", sum(rip.Area_Elevations[index_2100: i-1]))
        end
    end

    #soil storage component
    Overlandflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = ripariansoilstorage(rip.Total_Effective_Precipitation, rip.Total_Interception_Evaporation, rip.Potential_Evaporation_Mean, rip.Riparian_Discharge / rip.Area_HRU, storages.Soil, parameters.beta, parameters.Ce, parameters.Drainagecapacity, parameters.Soilstoragecapacity)
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)
    GWflow = 0
    #calculate total outflow
    Flows = Fast_Discharge + GWflow + rip.Total_Interception_Evaporation + Soil_Evaporation
    @assert rip.Potential_Evaporation_Mean - Soil_Evaporation - rip.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges according to areal fraction
    Fast_Discharge = Fast_Discharge * rip.Area_HRU
    Total_Interception_Evaporation = rip.Total_Interception_Evaporation * rip.Area_HRU
    Soil_Evaporation = Soil_Evaporation * rip.Area_HRU
    # return water flows, evaporation, and states of the storage components
    rip_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    rip_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)

    rip.Total_Interception_Evaporation = 0
    rip.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert rip_out.Fast_Discharge >= 0
    @assert rip_out.GWflow >= 0
    @assert rip_out.Soil_Evaporation >= 0
    @assert rip_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert rip_storages.Fast >= 0
    @assert rip_storages.Interception >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Snow >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Soil >= 0
    @assert rip_storages.Interception <= ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity
    #@assert ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - rip_storages.Interception >= -10^(-10) * ones(rip.Nr_Elevationbands)
    @assert rip_storages.Soil <= parameters.Soilstoragecapacity
    #@assert parameters.Soilstoragecapacity - rip_storages.Soil >= -10^(-10)

    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    for i in 1: rip.Nr_Elevationbands
        Precipitation += rip.Precipitation[i] * rip.Area_Elevations[i]
        Interception_Storage_New += rip_storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_New += rip_storages.Snow[i] * rip.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_Old += Snow_Old[i] * rip.Area_Elevations[i]
    end

    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (rip_storages.Fast - storages.Fast) + (rip_storages.Soil -storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)
    # println("Area ELevations ", rip.Area_Elevations)
    # println("difference snow storage ", Snow_Storage_New -Snow_Storage_Old)
    # println("difference snow storage mm ", rip_storages.Snow - Snow_Old)
    # println("old ", Snow_Old)
    # println("new ", rip_storages.Snow)
    # println("Precipitation ", Precipitation)
    # println("check ", Precipitation + rip.Riparian_Discharge / rip.Area_HRU - (Flows + All_Storages))
    @assert -0.00000001 <= Precipitation + rip.Riparian_Discharge / rip.Area_HRU - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation * rip.Area_HRU + rip.Riparian_Discharge) - (Flows_Area + All_Storages * rip.Area_HRU) <= 0.00000001

    return rip_out, rip_storages, Precipitation, All_Storages
end

"""
Computes the fluxes and changes in storages in a riparian hydrological response unit of the model.

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function riparianHRU(rip::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    # riparian HRU has no preferential flow

    @assert rip.Total_Interception_Evaporation == 0
    @assert rip.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(rip.Area_Elevations) <= 1.000001 || sum(rip.Area_Elevations) == 0
    @assert rip.Area_HRU >= 0 and <= 1
    @assert rip.Nr_Elevationbands >= 1
    Snow = zeros(rip.Nr_Elevationbands)
    Interception = zeros(rip.Nr_Elevationbands)
    Snow_Cover = zeros(rip.Nr_Elevationbands)
    for i in 1 : rip.Nr_Elevationbands
        # snow component
        @assert rip.Area_Glacier[i] == 0.0
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(rip.Area_Glacier[i], rip.Precipitation[i], rip.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        @assert Melt_Glacier == 0.0
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(rip.Potential_Evaporation_Mean, rip.Precipitation[i], rip.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        rip.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * rip.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        rip.Total_Interception_Evaporation::Float64 += Interception_Evaporation * rip.Area_Elevations[i]
        #get snow cover extent
        if Snow[i] > 1
            Snow_Cover[i] = 1
        end
    end

    #soil storage component
    Overlandflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = ripariansoilstorage(rip.Total_Effective_Precipitation, rip.Total_Interception_Evaporation, rip.Potential_Evaporation_Mean, rip.Riparian_Discharge / rip.Area_HRU, storages.Soil, parameters.beta, parameters.Ce, parameters.Drainagecapacity, parameters.Soilstoragecapacity)
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)
    GWflow = 0
    #calculate total outflow
    Flows = Fast_Discharge + GWflow + rip.Total_Interception_Evaporation + Soil_Evaporation
    @assert rip.Potential_Evaporation_Mean - Soil_Evaporation - rip.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges according to areal fraction
    Fast_Discharge = Fast_Discharge * rip.Area_HRU
    Total_Interception_Evaporation = rip.Total_Interception_Evaporation * rip.Area_HRU
    Soil_Evaporation = Soil_Evaporation * rip.Area_HRU
    # return water flows, evaporation, and states of the storage components
    rip_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    rip_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)

    rip.Total_Interception_Evaporation = 0
    rip.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert rip_out.Fast_Discharge >= 0
    @assert rip_out.GWflow >= 0
    @assert rip_out.Soil_Evaporation >= 0
    @assert rip_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert rip_storages.Fast >= 0
    @assert rip_storages.Interception >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Snow >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Soil >= 0
    @assert rip_storages.Interception <= ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity
    #@assert ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - rip_storages.Interception >= -10^(-10) * ones(rip.Nr_Elevationbands)
    @assert rip_storages.Soil <= parameters.Soilstoragecapacity
    #@assert parameters.Soilstoragecapacity - rip_storages.Soil >= -10^(-10)

    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    for i in 1: rip.Nr_Elevationbands
        Precipitation += rip.Precipitation[i] * rip.Area_Elevations[i]
        Interception_Storage_New += rip_storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_New += rip_storages.Snow[i] * rip.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_Old += storages.Snow[i] * rip.Area_Elevations[i]
    end

    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (rip_storages.Fast - storages.Fast) + (rip_storages.Soil -storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)
    @assert -0.00000001 <= Precipitation + rip.Riparian_Discharge / rip.Area_HRU - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation * rip.Area_HRU + rip.Riparian_Discharge) - (Flows_Area + All_Storages * rip.Area_HRU) <= 0.00000001

    return rip_out, rip_storages, Precipitation, All_Storages
end

"""
Computes the fluxes and changes in storages in a hillslope hydrological response unit of the model.

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function hillslopeHRU_future(hill::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    @assert hill.Total_Interception_Evaporation == 0
    @assert hill.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(hill.Area_Elevations) <= 1.000001 || sum(hill.Area_Elevations) == 0
    @assert hill.Area_HRU >= 0 and <= 1
    @assert hill.Nr_Elevationbands >= 1
    # define Arrays for Snow and Interception storage
    Snow = zeros(hill.Nr_Elevationbands)
    Interception = zeros(hill.Nr_Elevationbands)
    Snow_Cover = zeros(hill.Nr_Elevationbands)
    Total_Melt = 0.0
    Total_Melt_Glacier = 0.0
    # calculate the model components that are height dependend
    for i in 1 : hill.Nr_Elevationbands
        # snow component
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(hill.Area_Glacier[i], hill.Precipitation[i], hill.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(hill.Potential_Evaporation_Mean, hill.Precipitation[i], hill.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        hill.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * hill.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        hill.Total_Interception_Evaporation::Float64 += (Interception_Evaporation * hill.Area_Elevations[i])
        Total_Melt_Glacier::Float64 += (Melt_Glacier * hill.Area_Elevations[i])
        Total_Melt += Melt * hill.Area_Elevations[i]
        # get the extend of snow cover, if there is snow in the snow bucket the snow cover is 1
        if Snow[i] > 1
            Snow_Cover[i] = 1
        else
            Snow_Cover[i] = hill.Area_Glacier[i]
        end
    end
    #soil storage component
    Overlandflow::Float64, Preferentialflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = soilstorage(hill.Total_Effective_Precipitation, hill.Total_Interception_Evaporation, hill.Potential_Evaporation_Mean, storages.Soil, parameters.beta, parameters.Ce, parameters.Ratio_Pref, parameters.Soilstoragecapacity)
    GWflow::Float64 = Preferentialflow
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)

    #calculate total outflow
    Flows = Fast_Discharge + GWflow + hill.Total_Interception_Evaporation + Soil_Evaporation
    @assert hill.Potential_Evaporation_Mean - Soil_Evaporation - hill.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges and evaporation fluxed according to areal fraction
    Total_Melt = Total_Melt * hill.Area_HRU
    Fast_Discharge = Fast_Discharge * hill.Area_HRU
    GWflow = GWflow * hill.Area_HRU
    Total_Interception_Evaporation = hill.Total_Interception_Evaporation * hill.Area_HRU
    Soil_Evaporation = Soil_Evaporation * hill.Area_HRU
    # returning all fluxes (evporative, discharge)
    hill_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    #returning all storages
    hill_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)
    # total interception should be zero again for next run
    hill.Total_Interception_Evaporation = 0
    hill.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert hill_out.Fast_Discharge >= 0
    @assert hill_out.GWflow >= 0
    @assert hill_out.Soil_Evaporation >= 0
    @assert hill_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert hill_storages.Fast >= 0
    @assert hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Snow >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Soil >= 0
    @assert ones(hill.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert parameters.Soilstoragecapacity - hill_storages.Soil >= 0

    #assertion water balance
    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    # get total storage of interception and snow over all elevation zones
    for i in 1: hill.Nr_Elevationbands
        Precipitation += hill.Precipitation[i] * hill.Area_Elevations[i]
        Interception_Storage_New += hill_storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_New += hill_storages.Snow[i] * hill.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_Old += storages.Snow[i] * hill.Area_Elevations[i]
    end
    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (hill_storages.Fast - storages.Fast) + (hill_storages.Soil - storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)

    #assert that water balance closes
    @assert -0.00000001 <= Precipitation + Total_Melt_Glacier - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation + Total_Melt_Glacier)  * hill.Area_HRU - (Flows_Area + All_Storages * hill.Area_HRU) <= 0.00000001
    return hill_out::Outflows, hill_storages::Storages, Precipitation, All_Storages, Total_Melt

end

function hillslopeHRU_future_snow_redistribution(hill::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    @assert hill.Total_Interception_Evaporation == 0
    @assert hill.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(hill.Area_Elevations) <= 1.000001 || sum(hill.Area_Elevations) == 0
    @assert hill.Area_HRU >= 0 and <= 1
    @assert hill.Nr_Elevationbands >= 1
    # define Arrays for Snow and Interception storage
    Snow = zeros(hill.Nr_Elevationbands)
    Interception = zeros(hill.Nr_Elevationbands)
    Snow_Cover = zeros(hill.Nr_Elevationbands)
    Total_Melt = 0.0
    Total_Melt_Glacier = 0.0
    Snow_Old = copy(storages.Snow)
    # calculate the model components that are height dependend
    for i in reverse(1 : hill.Nr_Elevationbands)
        # snow component
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(hill.Area_Glacier[i], hill.Precipitation[i], hill.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(hill.Potential_Evaporation_Mean, hill.Precipitation[i], hill.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        hill.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * hill.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        hill.Total_Interception_Evaporation::Float64 += (Interception_Evaporation * hill.Area_Elevations[i])
        Total_Melt_Glacier::Float64 += (Melt_Glacier * hill.Area_Elevations[i])
        Total_Melt += Melt * hill.Area_Elevations[i]
        # get the extend of snow cover, if there is snow in the snow bucket the snow cover is 1
        if Snow[i] > 1
            Snow_Cover[i] = 1
        else
            Snow_Cover[i] = hill.Area_Glacier[i]
        end
        Current_Elevation = (hill.Elevation_Count[i]-1)* 200 + hill.Catchment_Elevation[1]
        if Current_Elevation >= hill.Snow_Redistribution[2] && Snow[i] > hill.Snow_Redistribution[1]
            # gives the index of the elevation count of 2100 m elevation band
            index_2100 = findall(x->x == (2100 - hill.Catchment_Elevation[1]) / 200 , hill.Elevation_Count)[1]
            # snow to be redistributed according to areal extend of elevation zone
            Snow_redistributed = (Snow[i] - hill.Snow_Redistribution[1]) * hill.Area_Elevations[i]
            # snow that has to be added to lower elevation zones in mm
            Snow_redistributed_mm = Snow_redistributed / sum(hill.Area_Elevations[index_2100: i-1])
            # add the snow to the lower elevation zones
            storages.Snow[index_2100: i-1] = storages.Snow[index_2100: i-1] .+ Snow_redistributed_mm
            # substracte the Snow from the higher elevation zone
            Snow[i] = hill.Snow_Redistribution[1]
            #println("redistribution ", Snow_redistributed)
            #println("sum area below ", sum(hill.Area_Elevations[index_2100: i-1]))
        end
    end
    #soil storage component
    Overlandflow::Float64, Preferentialflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = soilstorage(hill.Total_Effective_Precipitation, hill.Total_Interception_Evaporation, hill.Potential_Evaporation_Mean, storages.Soil, parameters.beta, parameters.Ce, parameters.Ratio_Pref, parameters.Soilstoragecapacity)
    GWflow::Float64 = Preferentialflow
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)

    #calculate total outflow
    Flows = Fast_Discharge + GWflow + hill.Total_Interception_Evaporation + Soil_Evaporation
    @assert hill.Potential_Evaporation_Mean - Soil_Evaporation - hill.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges and evaporation fluxed according to areal fraction
    Total_Melt = Total_Melt * hill.Area_HRU
    Fast_Discharge = Fast_Discharge * hill.Area_HRU
    GWflow = GWflow * hill.Area_HRU
    Total_Interception_Evaporation = hill.Total_Interception_Evaporation * hill.Area_HRU
    Soil_Evaporation = Soil_Evaporation * hill.Area_HRU
    # returning all fluxes (evporative, discharge)
    hill_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    #returning all storages
    hill_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)
    # total interception should be zero again for next run
    hill.Total_Interception_Evaporation = 0
    hill.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert hill_out.Fast_Discharge >= 0
    @assert hill_out.GWflow >= 0
    @assert hill_out.Soil_Evaporation >= 0
    @assert hill_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert hill_storages.Fast >= 0
    @assert hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Snow >= zeros(hill.Nr_Elevationbands)
    @assert hill_storages.Soil >= 0
    @assert ones(hill.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - hill_storages.Interception >= zeros(hill.Nr_Elevationbands)
    @assert parameters.Soilstoragecapacity - hill_storages.Soil >= 0

    #assertion water balance
    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    # get total storage of interception and snow over all elevation zones
    for i in 1: hill.Nr_Elevationbands
        Precipitation += hill.Precipitation[i] * hill.Area_Elevations[i]
        Interception_Storage_New += hill_storages.Interception[i] * hill.Area_Elevations[i]
        Snow_Storage_New += hill_storages.Snow[i] * hill.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * hill.Area_Elevations[i]
        #Snow_Storage_Old += storages.Snow[i] * hill.Area_Elevations[i]
        Snow_Storage_Old += Snow_Old[i] * hill.Area_Elevations[i]
    end
    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (hill_storages.Fast - storages.Fast) + (hill_storages.Soil - storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)

    #assert that water balance closes
    @assert -0.00000001 <= Precipitation + Total_Melt_Glacier - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation + Total_Melt_Glacier)  * hill.Area_HRU - (Flows_Area + All_Storages * hill.Area_HRU) <= 0.00000001
    return hill_out::Outflows, hill_storages::Storages, Precipitation, All_Storages, Total_Melt

end
"""
Computes the fluxes and changes in storages in a riparian hydrological response unit of the model taking into account snow redistribution

$(SIGNATURES)

The function returns the flxues leaving the HRU and the new amounts of water stored in each model component. Also it returns the Precipitation the sum of water stored in the HRU during the timestep (new_Storage-old_Storage)
Function needs inputs to be in area independent units (e.g. mm)
"""
function riparianHRU_future(rip::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    # riparian HRU has no preferential flow

    @assert rip.Total_Interception_Evaporation == 0
    @assert rip.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(rip.Area_Elevations) <= 1.000001 || sum(rip.Area_Elevations) == 0
    @assert rip.Area_HRU >= 0 and <= 1
    @assert rip.Nr_Elevationbands >= 1
    Snow = zeros(rip.Nr_Elevationbands)
    Interception = zeros(rip.Nr_Elevationbands)
    Snow_Cover = zeros(rip.Nr_Elevationbands)
    Total_Melt = 0.0
    for i in 1 : rip.Nr_Elevationbands
        # snow component
        @assert rip.Area_Glacier[i] == 0.0
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(rip.Area_Glacier[i], rip.Precipitation[i], rip.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        @assert Melt_Glacier == 0.0
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(rip.Potential_Evaporation_Mean, rip.Precipitation[i], rip.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        rip.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * rip.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        rip.Total_Interception_Evaporation::Float64 += Interception_Evaporation * rip.Area_Elevations[i]
        Total_Melt += Melt * rip.Area_Elevations[i]
        #get snow cover extent
        if Snow[i] > 1
            Snow_Cover[i] = 1
        end
    end

    #soil storage component
    Overlandflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = ripariansoilstorage(rip.Total_Effective_Precipitation, rip.Total_Interception_Evaporation, rip.Potential_Evaporation_Mean, rip.Riparian_Discharge / rip.Area_HRU, storages.Soil, parameters.beta, parameters.Ce, parameters.Drainagecapacity, parameters.Soilstoragecapacity)
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)
    GWflow = 0
    #calculate total outflow
    Flows = Fast_Discharge + GWflow + rip.Total_Interception_Evaporation + Soil_Evaporation
    @assert rip.Potential_Evaporation_Mean - Soil_Evaporation - rip.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges according to areal fraction
    Total_Melt = Total_Melt * rip.Area_HRU
    Fast_Discharge = Fast_Discharge * rip.Area_HRU
    Total_Interception_Evaporation = rip.Total_Interception_Evaporation * rip.Area_HRU
    Soil_Evaporation = Soil_Evaporation * rip.Area_HRU
    # return water flows, evaporation, and states of the storage components
    rip_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    rip_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)

    rip.Total_Interception_Evaporation = 0
    rip.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert rip_out.Fast_Discharge >= 0
    @assert rip_out.GWflow >= 0
    @assert rip_out.Soil_Evaporation >= 0
    @assert rip_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert rip_storages.Fast >= 0
    @assert rip_storages.Interception >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Snow >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Soil >= 0
    @assert rip_storages.Interception <= ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity
    #@assert ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - rip_storages.Interception >= -10^(-10) * ones(rip.Nr_Elevationbands)
    @assert rip_storages.Soil <= parameters.Soilstoragecapacity
    #@assert parameters.Soilstoragecapacity - rip_storages.Soil >= -10^(-10)

    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    for i in 1: rip.Nr_Elevationbands
        Precipitation += rip.Precipitation[i] * rip.Area_Elevations[i]
        Interception_Storage_New += rip_storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_New += rip_storages.Snow[i] * rip.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_Old += storages.Snow[i] * rip.Area_Elevations[i]
    end

    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (rip_storages.Fast - storages.Fast) + (rip_storages.Soil -storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)
    @assert -0.00000001 <= Precipitation + rip.Riparian_Discharge / rip.Area_HRU - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation * rip.Area_HRU + rip.Riparian_Discharge) - (Flows_Area + All_Storages * rip.Area_HRU) <= 0.00000001

    return rip_out, rip_storages, Precipitation, All_Storages, Total_Melt::Float64
end

function riparianHRU_future_snow_redistribution(rip::HRU_Input, storages::Storages, parameters::Parameters)
    # Are_Elevations gives the areal percentage of each elevation band. The sum has to be 1
    # Area_elevations, Precipitation, Temp_elevation, Snowstorage, Interceptionstorage has to be array of length Nr_Elevationbands
    # riparian HRU has no preferential flow

    @assert rip.Total_Interception_Evaporation == 0
    @assert rip.Total_Effective_Precipitation == 0
    @assert 0.9999999 <= sum(rip.Area_Elevations) <= 1.000001 || sum(rip.Area_Elevations) == 0
    @assert rip.Area_HRU >= 0 and <= 1
    @assert rip.Nr_Elevationbands >= 1
    Snow = zeros(rip.Nr_Elevationbands)
    Interception = zeros(rip.Nr_Elevationbands)
    Snow_Cover = zeros(rip.Nr_Elevationbands)
    Total_Melt = 0.0
    Snow_Old = copy(storages.Snow)
    for i in reverse(1 : rip.Nr_Elevationbands)
        # snow component
        @assert rip.Area_Glacier[i] == 0.0
        Melt_Glacier::Float64, Melt::Float64, Snow[i]::Float64 = snow(rip.Area_Glacier[i], rip.Precipitation[i], rip.Temp_Elevation[i], storages.Snow[i], parameters.Meltfactor, parameters.Mm, parameters.Temp_Thresh)
        @assert Melt_Glacier == 0.0
        #interception component
        Effective_Precipitation::Float64, Interception_Evaporation::Float64, Interception[i]::Float64 = interception(rip.Potential_Evaporation_Mean, rip.Precipitation[i], rip.Temp_Elevation[i], storages.Interception[i], parameters.Interceptionstoragecapacity, parameters.Temp_Thresh)
        # the melt, effective precipitation and evaporation can be summed up over all elevations according to the areal extent
        rip.Total_Effective_Precipitation::Float64 += (Effective_Precipitation + Melt) * rip.Area_Elevations[i]
        #global Total_Melt += (Melt * Area_Elevations[i])
        rip.Total_Interception_Evaporation::Float64 += Interception_Evaporation * rip.Area_Elevations[i]
        Total_Melt += Melt * rip.Area_Elevations[i]
        #get snow cover extent
        if Snow[i] > 1
            Snow_Cover[i] = 1
        end
        Current_Elevation = (rip.Elevation_Count[i]-1)* 200 + rip.Catchment_Elevation[1]
        if Current_Elevation >= rip.Snow_Redistribution[2] && Snow[i] > rip.Snow_Redistribution[1]
            # gives the index of the elevation count of 2100 m elevation band
            index_2100 = findall(x->x == (2100 - rip.Catchment_Elevation[1]) / 200 , rip.Elevation_Count)[1]
            # snow to be redistributed according to areal extend of elevation zone
            Snow_redistributed = (Snow[i] - rip.Snow_Redistribution[1]) * rip.Area_Elevations[i]
            # snow that has to be added to lower elevation zones in mm
            Snow_redistributed_mm = Snow_redistributed / sum(rip.Area_Elevations[index_2100: i-1])
            # add the snow to the lower elevation zones
            storages.Snow[index_2100: i-1] = storages.Snow[index_2100: i-1] .+ Snow_redistributed_mm
            # substracte the Snow from the higher elevation zone
            Snow[i] = rip.Snow_Redistribution[1]
            # println("redistribution ", Snow_redistributed)
            # println("sum area below ", sum(rip.Area_Elevations[index_2100: i-1]))
        end
    end

    #soil storage component
    Overlandflow::Float64, Soil_Evaporation::Float64, Soilstorage::Float64 = ripariansoilstorage(rip.Total_Effective_Precipitation, rip.Total_Interception_Evaporation, rip.Potential_Evaporation_Mean, rip.Riparian_Discharge / rip.Area_HRU, storages.Soil, parameters.beta, parameters.Ce, parameters.Drainagecapacity, parameters.Soilstoragecapacity)
    #fast storage
    Fast_Discharge::Float64, Faststorage::Float64 = faststorage(Overlandflow, storages.Fast, parameters.Kf)
    GWflow = 0
    #calculate total outflow
    Flows = Fast_Discharge + GWflow + rip.Total_Interception_Evaporation + Soil_Evaporation
    @assert rip.Potential_Evaporation_Mean - Soil_Evaporation - rip.Total_Interception_Evaporation >= - 5 * eps(Float64)
    # change discharges according to areal fraction
    Total_Melt = Total_Melt * rip.Area_HRU
    Fast_Discharge = Fast_Discharge * rip.Area_HRU
    Total_Interception_Evaporation = rip.Total_Interception_Evaporation * rip.Area_HRU
    Soil_Evaporation = Soil_Evaporation * rip.Area_HRU
    # return water flows, evaporation, and states of the storage components
    rip_out = Outflows(Fast_Discharge, GWflow, Soil_Evaporation, Total_Interception_Evaporation)
    rip_storages = Storages(Faststorage, Interception, Snow, Snow_Cover, Soilstorage)

    rip.Total_Interception_Evaporation = 0
    rip.Total_Effective_Precipitation = 0

    #assertions for the outflows
    @assert rip_out.Fast_Discharge >= 0
    @assert rip_out.GWflow >= 0
    @assert rip_out.Soil_Evaporation >= 0
    @assert rip_out.Interception_Evaporation >= 0
    #assertions for the storages
    @assert rip_storages.Fast >= 0
    @assert rip_storages.Interception >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Snow >= zeros(rip.Nr_Elevationbands)
    @assert rip_storages.Soil >= 0
    @assert rip_storages.Interception <= ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity
    #@assert ones(rip.Nr_Elevationbands) * parameters.Interceptionstoragecapacity - rip_storages.Interception >= -10^(-10) * ones(rip.Nr_Elevationbands)
    @assert rip_storages.Soil <= parameters.Soilstoragecapacity
    #@assert parameters.Soilstoragecapacity - rip_storages.Soil >= -10^(-10)

    Precipitation = 0
    Interception_Storage_New = 0
    Snow_Storage_New = 0
    Interception_Storage_Old = 0
    Snow_Storage_Old = 0
    for i in 1: rip.Nr_Elevationbands
        Precipitation += rip.Precipitation[i] * rip.Area_Elevations[i]
        Interception_Storage_New += rip_storages.Interception[i] * rip.Area_Elevations[i]
        Snow_Storage_New += rip_storages.Snow[i] * rip.Area_Elevations[i]
        Interception_Storage_Old += storages.Interception[i] * rip.Area_Elevations[i]
        #Snow_Storage_Old += storages.Snow[i] * rip.Area_Elevations[i]
        Snow_Storage_Old += Snow_Old[i] * rip.Area_Elevations[i]
    end

    Flows_Area = Fast_Discharge + GWflow + Total_Interception_Evaporation + Soil_Evaporation
    All_Storages = (rip_storages.Fast - storages.Fast) + (rip_storages.Soil -storages.Soil) + (Snow_Storage_New -Snow_Storage_Old) + (Interception_Storage_New - Interception_Storage_Old)
    @assert -0.00000001 <= Precipitation + rip.Riparian_Discharge / rip.Area_HRU - (Flows + All_Storages) <= 0.00000001
    @assert -0.00000001 <= (Precipitation * rip.Area_HRU + rip.Riparian_Discharge) - (Flows_Area + All_Storages * rip.Area_HRU) <= 0.00000001

    return rip_out, rip_storages, Precipitation, All_Storages, Total_Melt::Float64
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
# Test MKL plans
using Pkg
Pkg.activate(".")
using FourierGPE, Test

#transform tests
# using Test

## test Parseval's theorem for wavefunctions
function parsevaltest(L,N)
X,K,dX,dK = makearrays(L,N)
T = makeT(X,K)

ψ = randn(N...) + im*randn(N...)

n1 = sum(abs2.(ψ))*prod(dX)
ϕ = T.Txk*ψ

n2 = sum(abs2.(ϕ))*prod(dK)
return n1,n2
end
##

L = (88.)
N = (128)
n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

L = (30.,20.)
N = (30,34)
n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

## Fails:
L = (30.,20.,10)
N = (30,34,24)
X,K,dX,dK = makearrays(L,N)
# T = makeT(X,K)
FFTW.set_num_threads(Sys.CPU_THREADS)
flags=FFTW.MEASURE
dim = length(X)
N = length.(X)
DX,DK = dfftall(X,K)
dμx = prod(DX)
dμk = prod(DK)
ψtest = ones(N...) |> complex

trans = (plan_fft,plan_fft!,plan_ifft,plan_ifft!)
meas = (dμx,dμx,dμk,dμk)
# flags = FFTW.MEASURE
args = ((ψtest,),(ψtest,),(ψtest,),(ψtest,))

# flags != FFTW.MEASURE && @info "Planning FFTs ..."
Txk,Txk!,Tkx,Tkx! = definetransforms(trans,args,meas,flags)

#Fails:
Mxk,Mxk!,Mkx,Mkx! = makeTMixed(X,K,flags=flags)

# @info "...Plans created."
# return Transforms{dim,j}(Txk,Txk!,Tkx,Tkx!,Mxk,Mxk!,Mkx,Mkx!,crandnpartition(dim,j))

# n1,n2 = parsevaltest(L,N)
# @test n1 ≈ n2

##
L = (30.,20.,10,50.)
N = (30,34,24,22)
n1,n2 = parsevaltest(L,N)
@test n1 ≈ n2

ext: jl


text:
#=
# A set of time steppers
=#

mutable struct Stepper
    type;       # The constant for the stepper type
    Nsteps;     # number of steps
    dt;         # step size
    cfl;        # CFL number
    stages;     # how many stages
    a;          # for RK steppers
    b;          # for RK steppers
    c;          # for RK steppers

    Stepper(t, c) = new(t, 0, 0, c, 0, [], [], []);
end

function init_stepper(x, stepper)
    if stepper.type == EULER_EXPLICIT
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.25;
        end
        stepper.dt = stepper.cfl*dxmin*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == EULER_IMPLICIT
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = (prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == CRANK_NICHOLSON
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.25;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 1;

        return stepper;

    elseif stepper.type == LSRK4
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 5;
        stepper.a = [0; -0.41789047449985195; -1.192151694642677; -1.6977846924715279; -1.5141834442571558];
        stepper.b = [0.14965902199922912; 0.37921031299962726;  0.8229550293869817; 0.6994504559491221; 0.15305724796815198];
        stepper.c = [0.0; 0.14965902199922912; 0.37040095736420475; 0.6222557631344432; 0.9582821306746903];

        return stepper;

    elseif stepper.type == RK4
        dxmin = abs(x[1,2]-x[1,1]); # TODO this only works for similar, square elements
        if stepper.cfl == 0
            stepper.cfl = 0.1;
        end
        stepper.dt = stepper.cfl*dxmin;
        stepper.Nsteps = ceil(prob.end_time/stepper.dt);
        stepper.dt = prob.end_time/stepper.Nsteps;
        stepper.stages = 4;
        stepper.a = [   0   0   0   0;
                        0.5 0   0   0;
                        0   0.5 0   0;
                        0   0   1   0];
        stepper.b = [1/6; 1/3; 1/3; 1/6];
        stepper.c = [0.0; 0.5; 0.5; 1];

        return stepper;
    end
end

function reformat_for_stepper(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs);
        newrhs = copy(rhs);
        for vi=1:length(rhs)
            if length(lhs[1][vi]) > 0 # this dof has a time derivative term
                (newlhs[vi], newrhs[vi]) = reformat_for_stepper((lhs[1][vi], lhs[2][vi]), rhs[vi], stepper);
            else # no time derivative for this dof
                newlhs[vi] = lhs[2][vi];
                newrhs[vi] = rhs[vi];
            end
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_IMPLICIT # lhs1 + dt*lhs2 = dt*rhs + lhs1
            for i=1:length(lhs[2])
                lhs[2][i] = lhs[2][i]*dt; # dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs[2]); # lhs1 + dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs + lhs1

        elseif stepper == EULER_EXPLICIT # lhs1 = dt*rhs - dt*lhs2 + lhs1
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]*dt; # -dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# dt*rhs - dt*lhs2
            append!(newrhs, lhs[1]);# dt*rhs - dt*lhs2 + lhs1

        elseif stepper == CRANK_NICHOLSON # lhs1 + 0.5*dt*lhs2 = dt*rhs - 0.5*dt*lhs2 + lhs1
            lhs2l = copy(lhs[2]);
            lhs2r = copy(lhs[2]);
            for i=1:length(lhs[2])
                lhs2l[i] = Basic(0.5)*lhs2l[i]*dt; # 0.5*dt*lhs2
                lhs2r[i] = Basic(-0.5)*lhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs2l); # lhs1 + 0.5*dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs - 0.5*dt*lhs2 + lhs1
            append!(newrhs, lhs2r);

        elseif stepper == LSRK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            # for i=1:length(rhs[1])
            #     rhs[1][i] = rhs[1][i]; # rhs
            # end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2

        elseif stepper == RK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            # for i=1:length(rhs[1])
            #     rhs[1][i] = rhs[1][i]; # rhs
            # end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2

        end
    end

    return (newlhs, newrhs);
end

function reformat_for_stepper(lhs, rhs, face_lhs, face_rhs,stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];
    newfacelhs = [];
    newfacerhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs);
        newrhs = copy(rhs);
        newfacelhs = copy(face_rhs);
        newfacerhs = copy(face_rhs);
        for vi=1:length(rhs)
            if length(lhs[1][vi]) > 0 # this dof has a time derivative term
                (newlhs[vi], newrhs[vi], newfacelhs[vi], newfacerhs[vi]) = reformat_for_stepper((lhs[1][vi], lhs[2][vi]), rhs[vi], face_lhs[vi], face_rhs[vi], stepper);
            else # no time derivative for this dof
                newlhs[vi] = lhs[2][vi];
                newrhs[vi] = rhs[vi];
                newfacelhs[vi] = face_lhs[vi];
                newfacerhs[vi] = face_rhs[vi];
            end
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_IMPLICIT # lhs1 + dt*lhs2 = dt*rhs + lhs1
            for i=1:length(lhs[2])
                lhs[2][i] = lhs[2][i]*dt; # dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end
            for i=1:length(face_rhs)
                face_rhs[i] = face_rhs[i]*dt; # dt*facelhs
            end
            for i=1:length(face_lhs)
                face_lhs[i] = face_lhs[i]*dt; # dt*facerhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs[2]); # lhs1 + dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs + lhs1
            newfacerhs = copy(face_rhs);# dt*facerhs
            newfacelhs = copy(face_lhs);# dt*facelhs

        elseif stepper == EULER_EXPLICIT # lhs1 = dt*rhs - dt*lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]*dt; # -dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end
            for i=1:length(face_rhs)
                face_rhs[i] = face_rhs[i]*dt; # dt*facerhs
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]*dt; # -dt*facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# dt*rhs - dt*lhs2 + lhs1
            #append!(newrhs, lhs[1]);
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# dt*facerhs - dt*facelhs
            newfacelhs = [];

        elseif stepper == CRANK_NICHOLSON # lhs1 + 0.5*dt*lhs2 = dt*rhs - 0.5*dt*lhs2 + lhs1
            lhs2l = copy(lhs[2]);
            lhs2r = copy(lhs[2]);
            facelhs2l = copy(face_lhs);
            facelhs2r = copy(face_lhs);
            for i=1:length(lhs[2])
                lhs2l[i] = Basic(0.5)*lhs2l[i]*dt; # 0.5*dt*lhs2
                lhs2r[i] = Basic(-0.5)*lhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(face_lhs)
                facelhs2l[i] = Basic(0.5)*facelhs2l[i]*dt; # 0.5*dt*lhs2
                facelhs2r[i] = Basic(-0.5)*facelhs2r[i]*dt; # -0.5*dt*lhs2
            end
            for i=1:length(rhs)
                rhs[i] = rhs[i]*dt; # dt*rhs
            end
            for i=1:length(face_rhs)
                face_rhs[i] = face_rhs[i]*dt; # dt*facelhs
            end

            newlhs = copy(lhs[1]);
            append!(newlhs, lhs2l); # lhs1 + 0.5*dt*lhs2
            newrhs = copy(rhs);
            append!(newrhs, lhs[1]);# dt*rhs - 0.5*dt*lhs2 + lhs1
            append!(newrhs, lhs2r);
            newfacelhs = facelhs2l;
            newfacerhs = face_rhs[1];
            append!(newfacerhs, facelhs2r);

        elseif stepper == LSRK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]; # -facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# facerhs - facelhs
            newfacelhs = [];

        elseif stepper == RK4 # (lhs1) : rhs - lhs2
            for i=1:length(lhs[2])
                lhs[2][i] = -lhs[2][i]; # -lhs2
            end
            for i=1:length(face_lhs)
                face_lhs[i] = -face_lhs[i]; # -facelhs
            end

            newlhs = copy(lhs[1]);# lhs1
            newrhs = copy(rhs);
            append!(newrhs, lhs[2]);# rhs - lhs2
            newfacerhs = copy(face_rhs);
            append!(newfacerhs, face_lhs);# facerhs - facelhs
            newfacelhs = [];

        end
    end

    return (newlhs, newrhs, newfacelhs, newfacerhs);
end

# Special version for FV. Assumes a Dt(u) term that is not explicitly included.
function reformat_for_stepper_fv(flhs, frhs, slhs, srhs, stepper)
    (newflhs, newfrhs) = reformat_for_stepper_fv_flux(flhs, frhs, stepper);
    (newslhs, newsrhs) = reformat_for_stepper_fv_source(slhs, srhs, stepper);

    return (newflhs, newfrhs, newslhs, newsrhs);
end

# Special version for FV. Flux term only
function reformat_for_stepper_fv_flux(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs); # copy just to set up the arrays right
        newrhs = copy(rhs);
        for vi=1:length(rhs);
            (newlhs[vi], newrhs[vi]) = reformat_for_stepper_fv_flux(lhs[vi], rhs[vi], stepper);
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_EXPLICIT || stepper == RK4 || stepper == LSRK4
            # du/dt = -lhs + rhs
            for i=1:length(lhs)
                lhs[i] = -lhs[i];
            end
            # for i=1:length(rhs)
            #     rhs[i] = -rhs[i]; # don't change sign for rhs, because already done by parser
            # end

            newlhs = []; # put everything in rhs for explicit steppers
            newrhs = copy(rhs);
            append!(newrhs, lhs);

        elseif stepper == EULER_IMPLICIT

        elseif stepper == CRANK_NICHOLSON

        end
    end

    return (newlhs, newrhs);
end

# Special version for FV. Assumes a Dt(u) term that is not explicitly included.
function reformat_for_stepper_fv_source(lhs, rhs, stepper)
    # rebuild the expressions depending on type of time stepper
    dt = symbols("dt");
    newlhs = [];
    newrhs = [];

    if typeof(rhs) <: Array && length(rhs) > 0 && typeof(rhs[1]) <: Array # recursively work on subarrays
        newlhs = copy(rhs); # copy just to set up the arrays right
        newrhs = copy(rhs);
        for vi=1:length(rhs);
            (newlhs[vi], newrhs[vi]) = reformat_for_stepper_fv_flux(lhs[vi], rhs[vi], stepper);
        end
    else
        # reformat depending on stepper type
        if stepper == EULER_EXPLICIT || stepper == RK4 || stepper == LSRK4
            # du/dt = -lhs + rhs
            for i=1:length(lhs)
                lhs[i] = -lhs[i];
            end
            for i=1:length(rhs)
                rhs[i] = -rhs[i]; # need to change sign here?
            end

            newlhs = []; # put everything in rhs for explicit steppers
            newrhs = copy(rhs);
            append!(newrhs, lhs);

        elseif stepper == EULER_IMPLICIT

        elseif stepper == CRANK_NICHOLSON

        end
    end

    return (newlhs, newrhs);
end


text:
struct FlattenIterator{Tname, Tvalue}
    name::Tname
    value::Tvalue
end

Base.length(iter::FlattenIterator) = _length(iter.value)
_length(a::AbstractArray) = sum(_length, a)
_length(a::AbstractArray{<:Number}) = length(a)
_length(::Number) = 1

Base.eltype(iter::FlattenIterator{String}) = Tuple{String, _eltype(typeof(iter.value))}
_eltype(::Type{TA}) where {TA <: AbstractArray} = _eltype(eltype(TA))
_eltype(::Type{T}) where {T <: Number} = T

@inline function Base.iterate(iter::FlattenIterator{String, <:Number}, i = 1)
    i === 1 && return (iter.name, iter.value), 2
    return nothing
end
@inline function Base.iterate(
    iter::FlattenIterator{String, <:AbstractArray{<:Number}},
    ind = (1,),
)
    i = ind[1]
    i > length(iter.value) && return nothing
    name = getname(iter, i)
    return (name, iter.value[i]), (i+1,)
end
@inline function Base.iterate(
    iter::FlattenIterator{String, T},
    ind = startind(T),
) where {T <: AbstractArray}
    i = ind[1]
    i > length(iter.value) && return nothing
    name = getname(iter, i)
    local out
    while i <= length(iter.value)
        v = iter.value[i]
        out = iterate(FlattenIterator(name, v), Base.tail(ind))
        out !== nothing && break
        i += 1
    end
    if out === nothing
        return nothing
    else
        return out[1], (i, out[2]...)
    end
end

@inline startind(::Type{<:AbstractArray{T}}) where {T} = (1, startind(T)...)
@inline startind(::Type{<:Number}) = ()
@inline startind(::Type{<:Any}) = throw("Type not supported.")
@inline function getname(iter::FlattenIterator, i::Int)
    name = string(ind2sub(size(iter.value), i))
    name = replace(name, "(" => "[");
    name = replace(name, ",)" => "]");
    name = replace(name, ")" => "]");
    name = iter.name * name
    return name
end
@inline ind2sub(v, i) = Tuple(CartesianIndices(v)[i])


text:
# Method specialisation singletons.
struct GeoPlot end
struct GeoZPlot end
# We only look at arrays with X, Y, Z dims here.
# Otherwise they fall back to DimensionalData.jl recipes
@recipe function f(A::AbstractGeoArray)
    ddplot(A) = DimArray(A; dims=_maybe_mapped(dims(A)))
    max_res = get(plotattributes, :max_res, 1000)
    if !(get(plotattributes, :seriestype, :none) in (:none, :heatmap, :contourf))
        DD.DimensionalPlot(), ddplot(A)
    elseif all(hasdim(A, (SpatialDim, SpatialDim)))
        # Heatmap or multiple heatmaps. Use GD recipes.
        A = _prepare(_subsample(A, max_res))
        GeoPlot(), A
    elseif hasdim(A, ZDim) && ndims(A) == 1
        # Z dim plot, but for spatial data we want Z on the Y axis
        GeoZPlot(), _prepare(A)
    else
        DD.DimensionalPlot(), ddplot(A)
    end
end
# Plot a sinlge 2d map
@recipe function f(::GeoPlot, A::GeoArray{T,2,<:Tuple{D1,D2}}) where {T,D1<:SpatialDim,D2<:SpatialDim}
    # If colorbar is close to symmetric (< 25% difference) use a symmetric
    # colormap and set symmetric limits so zero shows up as a neutral color.

    yguide, xguide = label(dims(A))

    y, x = map(_prepare, dims(A))

    rdt = DD.refdims_title(A; issingle=true)
    :title --> (rdt === "" ? _maybename(A) : _maybename(A) * "\n" * rdt)
    :xguide --> xguide
    :yguide --> yguide
    :grid --> true
    :gridalpha --> 0.2
    # :guidefontsize --> 10
    # :titlefontsize --> 10
    # :tickfontsize --> 6
    # :colorbar_title --> name(A)
    :linewidth --> 1
    :colorbar_titlefontsize --> 9
    :colorbar_tickfontcolor --> RGB(0.3)
    :tickfontcolor --> RGB(0.3)
    :tickcolor --> RGB(0.3)
    :tick_direction --> :out
    :framestyle --> :box
    :foreground_color_axis --> RGB(0.3)
    :foreground_color_border --> RGB(0.3)
    :seriescolor --> :curl

    if all(d -> lookup(d) isa Mapped, (x, y))
        :xlims --> mappedbounds(x)
        :ylims --> mappedbounds(y)
        :aspect_ratio --> :equal
    else
        :xlims --> bounds(A, x)
        :ylims --> bounds(A, y)
        bnds = bounds(A, (D1, D2))
        s1, s2 = map(((l, u),) -> (u - l), bnds) ./ (size(A, D1), size(A, D2))
        square_pixels = s2 / s1
        :aspect_ratio --> square_pixels
    end

    if get(plotattributes, :seriestype, :none) == :contourf
        clims = extrema(skipmissing(A))
        :levels --> range(clims[1], clims[2], length=20)
        index(x), index(y), clamp.(A, clims[1], clims[2])
    else
        :seriestype --> :heatmap
        index(x), index(y), parent(A)
    end
end

# Plot a vertical 1d line
@recipe function f(::GeoZPlot, A::GeoArray{T,1}) where T
    z_dim = dims(A, ZDim)
    yguide = label(z_dim)
    xguide = label(A)
    rdt = DD.refdims_title(A; issingle=true)
    :title --> (rdt == "" ? _maybename(A) : _maybename(A) * "\n" * rdt)
    :xguide --> xguide
    :yguide --> yguide
    :label --> ""
    z = map(_prepare, dims(A))
    parent(A), index(z)
end

# Plot 3d arrays as multiple tiled plots
@recipe function f(::GeoPlot, A::GeoArray{T,3,<:Tuple{<:SpatialDim,<:SpatialDim,D}}) where {T,D}
    nplots = size(A, 3)
    if nplots > 1
        ncols, nrows = _balance_grid(nplots)
        :layout --> (ncols, nrows)
        # link --> :both
        # :colorbar := false
        titles = string.(index(A, D))
        for r in 1:nrows, c in 1:ncols
            i = (r + (c - 1) * nrows)
            @series begin
                :titlefontsize := 9
                :tickfontsize := 7
                subplot := i
                if c != ncols || r != 1
                    :xformatter := _ -> ""
                    :yformatter := _ -> ""
                    :xguide := ""
                    :yguide := ""
                end
                if i <= nplots
                    title := titles[i]
                    GeoPlot(), A[:, :, i]
                else
                    :framestyle := :none
                    :legend := :none
                    []
                end
            end
        end
    else
        GeoPlot(), A[:, :, 1]
    end
end

function _balance_grid(nplots)
    ncols = (nplots - 1) ÷ ceil(Int, sqrt(nplots)) + 1
    nrows = (nplots - 1) ÷ ncols + 1
    return ncols, nrows
end

# We only look at arrays with X, Y, Z dims here.
# Otherwise they fall back to DimensionalData.jl recipes
@recipe function f(st::AbstractGeoStack)
    nplots = length(keys(st))
    if nplots > 1
        ncols, nrows = _balance_grid(nplots)
        :layout --> (ncols, nrows)
        # colorbar := false
        max_res = get(plotattributes, :max_res, 1000/(max(nrows, ncols)))

        l = DD.layers(st)
        for r in 1:nrows, c in 1:ncols
            i = (r + (c - 1) * nrows)
            @series begin
                :titlefontsize := 9
                :tickfontsize := 7
                subplot := i
                if c != ncols || r != 1
                    :xformatter := _ -> ""
                    :yformatter := _ -> ""
                    :xguide := ""
                    :yguide := ""
                end
                if i <= nplots
                    A = l[i]
                    title := string(keys(st)[i])
                    if length(dims(A, (XDim, YDim))) > 0
                        # Get a view of the first slice of the X/Y dimension
                        ods = otherdims(A, (X, Y))
                        if length(ods) > 0
                            od1s = map(d -> DD.basetypeof(d)(firstindex(d)), ods)
                            A = view(A, od1s...)
                        end
                        GeoPlot(), _prepare(_subsample(A, max_res))
                    else
                        framestyle := :none
                        legend := :none
                        []
                    end
                else
                    framestyle := :none
                    legend := :none
                    []
                end
            end
        end
    else
        first(st)
    end
end

# Plots.jl heatmaps pixels are centered.
# So we should center the index, and use the projected value.
_prepare(d::Dimension) = d |> _maybe_shift |> _maybe_mapped
# Convert arrays to a consistent missing value and Forward array order
function _prepare(A::AbstractGeoArray)
    reorder(A, DD.ForwardOrdered) |>
    a -> permutedims(a, DD.commondims(>:, (ZDim, YDim, XDim, TimeDim, Dimension), dims(A))) |>
    a -> replace_missing(a, missing)
end

function _subsample(A, max_res)
    ssdims = dims(A, (XDim, YDim, ZDim))[1:2]
    # Aggregate based on the number of pixels
    s1, s2 = size(A, ssdims[1]), size(A, ssdims[2])
    ag = floor(Int, max(s1, s2) / max_res) + 1
    if ag == 1
        return read(A)
    else
        return A[DD.basetypeof(ssdims[1])(1:ag:s1), DD.basetypeof(ssdims[2])(1:ag:s2)]
    end
end

_maybename(A) = _maybename(name(A))
_maybename(n::Name{N}) where N = _maybename(N)
_maybename(n::NoName) = ""
_maybename(n::Symbol) = string(n)

_maybe_replace_missing(A::AbstractArray{<:AbstractFloat}) = replace_missing(A, eltype(A)(NaN))
_maybe_replace_missing(A) = A

_maybe_shift(d) = _maybe_shift(sampling(d), d)
_maybe_shift(::Intervals, d) = DD.maybeshiftlocus(Center(), d)
_maybe_shift(sampling, d) = d

_maybe_mapped(dims::Tuple) = map(_maybe_mapped, dims)
_maybe_mapped(dim::Dimension) = _maybe_mapped(lookup(dim), dim)
_maybe_mapped(lookup::LookupArray, dim::Dimension) = dim
_maybe_mapped(lookup::Projected, dim::Dimension) = _maybe_mapped(mappedcrs(lookup), dim)
_maybe_mapped(::Nothing, dim::Dimension) = dim
_maybe_mapped(::GeoFormat, dim::Dimension) = convertlookup(Mapped, dim)

# We don't show the Band label for a single-band raster,
# it's just not interesting information.
function DD.refdims_title(refdim::Band; issingle=false)
    if issingle
        ""
    else
        string(name(refdim), ": ", DD.refdims_title(lookup(refdim), refdim))
    end
end


text:
## Testing: Compare it GLMNet on toy problem
@testset "Compare with GLMNet" begin
    N = 100
    srand(1)
    y = collect(1:N) + randn(N)*10
    X = [1:N (1:N)+randn(N)*5 (1:N)+randn(N)*10 (1:N)+randn(N)*20]
    X = X .- mean(X, 1)
    y = y .- mean(y)
    norm_x = sqrt(N)*std(X, 1, corrected=false)
    std_y = std(y, corrected=false)
    yn = y./std_y
    Xn = X./norm_x

    # Test for a single value
    λ = 0.001
    path = GLMNet.glmnet(Xn, yn; lambda=[λ], intercept=false, standardize=false,
              tol=1e-20)
    β, _ = NarmaxLasso.pathwise_coordinate_optimization(Xn, yn, λ, tol=1e-20)
    @test vec(Matrix{Float64}(path.betas)) ≈ β atol=1e-5

    # Test for multiple value
    λ_list = collect(100:-1:1)*0.01
    path = GLMNet.glmnet(Xn, yn; lambda=λ_list, intercept=false,
              standardize=false, tol=1e-20)
    β, _ = NarmaxLasso.pathwise_coordinate_optimization(Xn, yn, λ_list, tol=1e-20)
    @test Matrix{Float64}(path.betas) ≈ β atol=1e-5
end

## Generate nonlinear System Inputs
function nonlinear_function(y, u, v)
    return ((0.8 - 0.5*exp(-y[1]^2))*y[1] - (0.3 + 0.9*exp(-y[1]^2))*y[2] + u[1]
            + 0.2*u[2] + 0.1*u[1]*u[2] + 0.1*v[1] + 0.3*v[2])
end

## Testing: Compare it GLMNet on dynamic example
@testset "Compare with GLMNet on Dynamic Data" begin
    # Generate input
    u = NarmaxLasso.generate_random_input(1000, 5);
    y = NarmaxLasso.simulate(u, nonlinear_function, [1, 2], [1, 2], [1, 2];
                             σv=0.1, σw=0.5);
    # Generate Regressor matrix
    X = NarmaxLasso.lag([u y], [collect(1:30),
                    collect(1:30)])
    N, p = size(X)
    # Normalize
    X = X .- mean(X, 1)
    y = y .- mean(y)
    norm_x = sqrt(N)*std(X, 1, corrected=false)
    std_y = std(y, corrected=false)
    y = y./std_y
    X = X./norm_x
    # Test single point
    λ = 1e-5
    path = GLMNet.glmnet(X, y, lambda=[λ], intercept=false,
                        standardize=false, tol=1e-20)
    β, λ= NarmaxLasso.pathwise_coordinate_optimization(X, y, λ, tol=1e-15)
    @test vec(Matrix{Float64}(path.betas)) ≈ β atol=1e-5
    # Test sequence of points
    λ = collect(100:-1:1)*0.00001
    path = GLMNet.glmnet(X, y, lambda=λ, intercept=false,
                         standardize=false, tol=1e-20)
    β, λ= NarmaxLasso.pathwise_coordinate_optimization(X, y, λ, tol=1e-15)
    @test Matrix{Float64}(path.betas) ≈ β atol=1e-5
end


text:
#Graphics functions
show(io::IO, w::AFWindow) = print(io, "ArrayFire Window")
show(io::IO, w::vcpp"af::Window") = print(io, "ArrayFire Window")

export AFWindow, window, setTitle, setPlot, setImage, setSurfacePlot

window() = AFWindow(af_window())
window(str::AbstractString) = AFWindow(af_window(str))
window(width::Integer, height::Integer, str::AbstractString) = AFWindow(af_window(width, height, str))

setTitle(w::AFWindow, str::AbstractString) = af_setTitle(w, str)
setImage(w::AFWindow, img::AFAbstractArray, str = "") = af_setImage(w, img/255.0, str)
setPlot(w::AFWindow, x::AFAbstractArray, y::AFAbstractArray, title = "") = af_setPlot(w, x, y, title)
setSurfacePlot(w::AFWindow, x::AFAbstractArray) = af_setSurface(w, x)
show(w::AFWindow) = af_show(w)


text:
#data type heirachy tree

#all tradeables have amt and qual
abstract tradeable

#function +(a::tradeable,b::tradeable)
#  return tradeable(a.amt+b.amt,(a.amt*a.qual+b.amt*b.qual)/(a.amt+b.amt))
#end

#all buildings have x::int,y::int,dx::int,dy::int,
abstract building
abstract retail <: building
abstract researchable <: building
abstract occupyable <: building
abstract headquarters <: building
abstract residential <: occupyable
abstract offices <: occupyable
abstract industrial <: building

max_stars_repo_name: JKSano/netgameproj

max_issues_repo_name: JKSano/netgameproj

max_forks_repo_name: JKSano/netgameproj


text:
@testset "@ctypedef" begin
	@eval @ctypedef EmptyStructTypedef @cstruct {
	} __packed__
	@test sizeof(EmptyStructTypedef) == 0

	@eval @ctypedef CcharStructTypedef @cstruct {
		c::Cchar
	} __packed__
	@test sizeof(CcharStructTypedef) == sizeof(Cchar)
	@test :c in propertynames(CcharStructTypedef)

	@eval @ctypedef CintUnionTypedef @cunion {
		i::Cint
	} __packed__
	@test sizeof(CintUnionTypedef) == sizeof(Cint)
	@test :i in propertynames(CintUnionTypedef)

	@eval @ctypedef CenumTypedef @cenum {
		CENUM_TYPEDEF_1,
		CENUM_TYPEDEF_2 = 1 << 8,
		CENUM_TYPEDEF_3 = CENUM_TYPEDEF_2,
	} __packed__
	@test sizeof(CenumTypedef) == 2
	@test CENUM_TYPEDEF_3 != CENUM_TYPEDEF_1
	@test CENUM_TYPEDEF_3 == CENUM_TYPEDEF_2

	@eval @ctypedef CuintArrayTypedef Cuint[32]
	@test sizeof(CuintArrayTypedef) == sizeof(Cuint)*32

	@eval @ctypedef CintTypedef Cint
	@test sizeof(CintTypedef) == sizeof(Cint)
end


text:
using Query, DataFrames, Feather, NamedTuples

testfile = joinpath(Pkg.dir("Feather"),"test", "data", "airquality.feather")

results = @from i in Feather.Source(testfile) begin
    @select i
    @collect DataFrame
end

println(results)


text:
@testset "SubTabular" begin
    @testset "Views of Series" begin
        s = Series([1,2,3])

        @test @inferred(view(s, 2:3))::SubSeries == Series(2=>2,3=>3)
        s2 = view(s, 2:3)
        @test @inferred(s2[2]) === 2
        @test @inferred(s2[2:3]) == s2
        @test @inferred(s2[:]) == s2

        @test @inferred(view(s, :))::SubSeries == s
        s3 = view(s, :)
        @test @inferred(s3[2]) === 2
        @test @inferred(s3[2:3]) == s2
        @test @inferred(s3[:]) == s

        # Views of SubSeries
        @test @inferred(view(view(s, :), :)) == s
        @test @inferred(view(view(s, 2:3), :)) == Series(2=>2,3=>3)
        @test @inferred(view(view(s, :), 2:3)) == Series(2=>2,3=>3)
        @test @inferred(view(view(s, 2:3), 2:3)) == Series(2=>2,3=>3)
    end

    @testset "Views of Table" begin
        t = Table([1 2 3; 4 5 6; 7 8 9; 10 11 12])

        @test @inferred(view(t, 2:4, 3))::SubSeries == Series(2=>6,3=>9,4=>12)
        s2 = view(t, 2:4, 3)
        @test @inferred(s2[2]) === 6
        @test @inferred(s2[2:4]) == s2
        @test @inferred(s2[:]) == s2

        @test @inferred(view(t, 4, 2:3))::SubSeries == Series(2=>11,3=>12)
        s3 = view(t, 4, 2:3)
        @test @inferred(s3[2]) === 11
        @test @inferred(s3[2:3]) == s3
        @test @inferred(s3[:]) == s3

        @test @inferred(view(t, :, 3))::SubSeries == Series(1=>3,2=>6,3=>9,4=>12)
        s4 = view(t, :, 3)
        @test @inferred(s4[2]) === 6
        @test @inferred(s4[2:4]) == s2
        @test @inferred(s4[:]) == s4

        @test @inferred(view(t, 4, :))::SubSeries == Series(1=>10,2=>11,3=>12)
        s5 = view(t, 4, :)
        @test @inferred(s5[2]) === 11
        @test @inferred(s5[2:3]) == s3
        @test @inferred(s5[:]) == s5

        @test @inferred(view(t, 2:4, :))::SubTable == Table([(2=>4, 3=>7, 4=>10), (2=>5, 3=>8, 4=>11), (2=>6, 3=>9, 4=>12)])
        t2 = view(t, 2:4, :)
        @test @inferred(t2[4,3]) === 12
        @test @inferred(t2[:,3]) == s2
        @test @inferred(t2[2:4,3]) == s2
        @test @inferred(t2[4,:]) == t[4,:]
        @test @inferred(t2[4,2:3]) == s3
        @test @inferred(t2[:,1:3]) == t2
        @test @inferred(t2[2:4,1:3]) == t2
        @test @inferred(t2[2:4,:]) == t2
        @test @inferred(t2[:,:]) == t2

        @test @inferred(view(t, :, 2:3))::SubTable == Table(2=>[2,5,8,11], 3=>[3,6,9,12])
        @test @inferred(view(t, 2:4, 2:3))::SubTable == Table(2=>(2=>5, 3=>8, 4=>11), 3=>(2=>6, 3=>9, 4=>12))
        @test @inferred(view(t, :, :))::SubTable == t

        # views of SubTabulars
        @test @inferred(view(view(t, 2:4, 3), 2:4))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, 3), 2:4))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, 2:4, 3), :))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, 3), :))::SubSeries == Series(1=>3,2=>6,3=>9,4=>12)

        @test @inferred(view(view(t, 4, 2:3), 2:3))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, 4, :), 2:3))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, 4, 2:3), :))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, 4, :), :))::SubSeries == Series(1=>10,2=>11,3=>12)

        @test @inferred(view(view(t, 2:4, :), 2:4, 3))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, 2:3), 4, 2:3))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, 2:4, :), 2:4, :))::SubTable == Table([(2=>4, 3=>7, 4=>10), (2=>5, 3=>8, 4=>11), (2=>6, 3=>9, 4=>12)])
        @test @inferred(view(view(t, :, 2:3), :, 2:3))::SubTable == Table(2=>[2,5,8,11], 3=>[3,6,9,12])
        @test @inferred(view(view(t, 2:4, 2:3), 2:4, 2:3))::SubTable == Table(2=>(2=>5, 3=>8, 4=>11), 3=>(2=>6, 3=>9, 4=>12))

        @test @inferred(view(view(t, 2:4, :), :, 3))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, 2:3), 4, :))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, 2:4, :), :, :))::SubTable == Table([(2=>4, 3=>7, 4=>10), (2=>5, 3=>8, 4=>11), (2=>6, 3=>9, 4=>12)])
        @test @inferred(view(view(t, :, 2:3), :, :))::SubTable == Table(2=>[2,5,8,11], 3=>[3,6,9,12])
        @test @inferred(view(view(t, 2:4, 2:3), :, :))::SubTable == Table(2=>(2=>5, 3=>8, 4=>11), 3=>(2=>6, 3=>9, 4=>12))

        @test @inferred(view(view(t, :, :), 2:4, 3))::SubSeries == Series(2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, :), 4, 2:3))::SubSeries == Series(2=>11,3=>12)
        @test @inferred(view(view(t, :, :), :, 3))::SubSeries == Series(1=>3,2=>6,3=>9,4=>12)
        @test @inferred(view(view(t, :, :), 4, :))::SubSeries == Series(1=>10,2=>11,3=>12)
        @test @inferred(view(view(t, :, :), 2:4, :))::SubTable == Table([(2=>4, 3=>7, 4=>10), (2=>5, 3=>8, 4=>11), (2=>6, 3=>9, 4=>12)])
        @test @inferred(view(view(t, :, :), :, 2:3))::SubTable == Table(2=>[2,5,8,11], 3=>[3,6,9,12])
        @test @inferred(view(view(t, :, :), 2:4, 2:3))::SubTable == Table(2=>(2=>5, 3=>8, 4=>11), 3=>(2=>6, 3=>9, 4=>12))
        @test @inferred(view(view(t, :, :), :, :))::SubTable == t
    end
end