diff --git a/Project.toml b/Project.toml index 113eaab..a8deb8d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,14 +1,12 @@ name = "QuanticsGrids" uuid = "634c7f73-3e90-4749-a1bd-001b8efc642d" authors = ["Ritter.Marc , Hiroshi Shinaoka and contributors"] -version = "0.4.0" +version = "0.5.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] -StaticArrays = "1" julia = "1.6" [extras] diff --git a/src/QuanticsGrids.jl b/src/QuanticsGrids.jl index 972e09b..913d786 100644 --- a/src/QuanticsGrids.jl +++ b/src/QuanticsGrids.jl @@ -1,13 +1,13 @@ module QuanticsGrids -using StaticArrays - export DiscretizedGrid, InherentDiscreteGrid export quantics_to_grididx, quantics_to_origcoord export grididx_to_quantics, grididx_to_origcoord export origcoord_to_quantics, origcoord_to_grididx -include("quantics.jl") +abstract type Grid{D} end + include("grid.jl") +include("grid_discretized.jl") end diff --git a/src/grid.jl b/src/grid.jl index afe3d1e..a89adbb 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -1,446 +1,406 @@ -# Grid for d-dimensional space -abstract type Grid{d} end +struct InherentDiscreteGrid{D} <: Grid{D} + Rs::NTuple{D,Int} + origin::NTuple{D,Int} + step::NTuple{D,Int} + variablenames::NTuple{D,Symbol} + base::Int + indextable::Vector{Vector{Tuple{Symbol,Int}}} + + # Lookup table: lookup_table[variablename_index][bitnumber] -> (site_index, position_in_site) + lookup_table::NTuple{D,Vector{Tuple{Int,Int}}} + maxgrididx::NTuple{D,Int} + + function InherentDiscreteGrid{D}( + Rs::NTuple{D,Int}, + origin::NTuple{D,Int}, + step::NTuple{D,Int}, + variablenames::NTuple{D,Symbol}, + base::Int, + indextable::Vector{Vector{Tuple{Symbol,Int}}} + ) where {D} + if !(D isa Int) + throw(ArgumentError(lazy"Got dimension $D, which is not an Int.")) + end + if base <= 1 + throw(ArgumentError(lazy"Got base = $base. base must be at least 2.")) + end + for d in 1:D + if !(step[d] >= 1) + throw(ArgumentError(lazy"Got step[$d] = $(step[d]). step must be at least 1.")) + end + if !rangecheck_R(Rs[d]; base) + throw(ArgumentError(lazy"Got Rs[$d] = $(Rs[d]) with base = $base. For all gridindices from 1 to base^R to fit into an Int, we must have base ^ R <= typemax(Int)")) + end + end + if !allunique(variablenames) + throw(ArgumentError(lazy"Got variablenames = $variablenames. variablenames must be unique.")) + end + + lookup_table = _build_lookup_table(Rs, indextable, variablenames) + maxgrididx = map(R -> base^R, Rs) + + return new{D}(Rs, origin, step, variablenames, base, indextable, lookup_table, maxgrididx) + end +end +# ============================================================================ +# Helper/utility functions +# ============================================================================ _convert_to_scalar_if_possible(x) = x _convert_to_scalar_if_possible(x::NTuple{1,T}) where {T} = first(x) -_to_tuple(::Val{d}, x::NTuple{d,T}) where {d,T} = x +_to_tuple(::Val{d}, x::NTuple{d}) where {d} = x _to_tuple(::Val{d}, x) where {d} = ntuple(i -> x, d) -function digitmax(d, base, unfoldingscheme::Symbol) - if unfoldingscheme === :fused - return base^d - else - return base +default_step(::Val{D}) where D = ntuple(Returns(1), D) + +default_origin(::Val{D}) where D = ntuple(Returns(1), D) + +function _build_lookup_table(Rs::NTuple{D,Int}, indextable::Vector{Vector{Tuple{Symbol,Int}}}, variablenames::NTuple{D,Symbol}) where D + lookup_table = ntuple(D) do d + if Rs[d] < 0 + throw(ArgumentError(lazy"Got Rs[$d] = $(Rs[d]). Rs must be non-negative.")) + end + Vector{Tuple{Int,Int}}(undef, Rs[d]) + end + + index_visited = [fill(false, Rs[d]) for d in 1:D] + + for (site_idx, quanticsindices) in pairs(indextable) + for (pos_in_site, qindex) in pairs(quanticsindices) + variablename, bitnumber = qindex + var_idx = findfirst(==(variablename), variablenames) + if isnothing(var_idx) + throw(ArgumentError(lazy"Index table contains unknown index $qindex. Valid variablenames are $variablenames.")) + elseif bitnumber > Rs[var_idx] + throw(ArgumentError(lazy"Index table contains quantics index $bitnumber of variable $variablename, but it must be smaller than or equal to the number of quantics indices for that variable, which is $(Rs[var_idx]).")) + elseif index_visited[var_idx][bitnumber] + throw(ArgumentError(lazy"Index table contains quantics index $bitnumber of variable $variablename more than once.")) + end + + lookup_table[var_idx][bitnumber] = (site_idx, pos_in_site) + index_visited[var_idx][bitnumber] = true + end end + + for (var_idx, visited) in enumerate(index_visited) + bitnumber = findfirst(==(false), visited) + if !isnothing(bitnumber) + throw(ArgumentError(lazy"Index table contains no site for quantics index $bitnumber of variable $(variablenames[var_idx]).")) + end + end + + return lookup_table end -function quanticslength(R, d, unfoldingscheme::Symbol) - if unfoldingscheme === :fused - return R - else - return R * d +function rangecheck_R(R::Int; base::Int=2)::Bool + # For all gridindices from 1 to base^R to fit into an Int64, + # we must have base ^ R <= typemax(Int) + result = 1 + for _ in 1:R + result <= typemax(Int) ÷ base || return false + result *= base + end + return true +end + +function _build_indextable(variablenames::NTuple{D,Symbol}, Rs::NTuple{D,Int}, unfoldingscheme::Symbol) where D + if !(unfoldingscheme in (:interleaved, :fused)) + throw(ArgumentError(lazy"Got unfoldingscheme = $unfoldingscheme. Supported are :interleaved and :fused.")) end + indextable = Vector{Tuple{Symbol,Int}}[] + + for bitnumber in 1:maximum(Rs; init=0) + if unfoldingscheme === :interleaved + _add_interleaved_indices!(indextable, variablenames, Rs, bitnumber) + elseif unfoldingscheme === :fused + _add_fused_indices!(indextable, variablenames, Rs, bitnumber) + end + end + + return indextable end -function localdimensions(g::Grid{d}) where {d} - return fill( - digitmax(d, g.base, g.unfoldingscheme), - quanticslength(g.R, d, g.unfoldingscheme), - ) +function _add_interleaved_indices!(indextable, variablenames::NTuple{D,Symbol}, Rs::NTuple{D,Int}, bitnumber) where D + for d in 1:D + bitnumber ∈ 1:Rs[d] || continue + qindex = (variablenames[d], bitnumber) + push!(indextable, [qindex]) + end end -function _rangecheck_R(R; base = 2)::Int - base * (BigInt(base)^R - 1) ÷ (base - 1) <= typemax(Int) || - error("R too large for base $base") +function _add_fused_indices!(indextable, variablenames::NTuple{D,Symbol}, Rs::NTuple{D,Int}, bitnumber) where D + indices_bitnumber = Tuple{Symbol,Int}[] + # Add dimensions in reverse order to match DiscretizedGrid convention + # where the first dimension varies fastest in fused quantics + for d in D:-1:1 + bitnumber ∈ 1:Rs[d] || continue + qindex = (variablenames[d], bitnumber) + push!(indices_bitnumber, qindex) + end + if !isempty(indices_bitnumber) + push!(indextable, indices_bitnumber) + end end -#=== -Conversion between grid index, original coordinate, quantics +function _quantics_to_grididx_general(g::InherentDiscreteGrid{D}, quantics) where D + base = g.base -* grid index => original coordinate (generic) -* grid index => quantics (generic) + return ntuple(D) do d + grididx = 1 + R_d = g.Rs[d] -* original coordinate => grid index (implemented for each grid type) -* original coordinate => quantics (generic) + for bitnumber in 1:R_d + site_idx, pos_in_site = g.lookup_table[d][bitnumber] + quantics_val = quantics[site_idx] + site_len = length(g.indextable[site_idx]) -* quantics => grididx (XXXX) -* quantics => origcoord (generic) -===# + temp = quantics_val - 1 + for _ in 1:(site_len-pos_in_site) + temp = div(temp, base) + end + digit = temp % base -""" -grid index => original coordinate -""" -function grididx_to_origcoord(g::Grid{d}, index::NTuple{d,Int}) where {d} - all(1 .<= index .<= g.base^g.R) || error("1 <= $index <= g.base^g.R") - return _convert_to_scalar_if_possible((index .- 1) .* grid_step(g) .+ grid_min(g)) + grididx += digit * base^(R_d - bitnumber) + end + grididx + end end -""" -grid index => original coordinate -""" -grididx_to_origcoord(g::Grid{1}, index::Int) = grididx_to_origcoord(g, (index,)) - -""" -grid index => quantics -""" -function grididx_to_quantics(g::Grid{d}, grididx::NTuple{d,Int}) where {d} - if g.unfoldingscheme === :fused - return index_to_quantics_fused(grididx, numdigits = g.R, base = g.base) - else - return fused_to_interleaved( - index_to_quantics_fused(grididx, numdigits = g.R, base = g.base), - d, - base = g.base, - ) +function _quantics_to_grididx_base2(g::InherentDiscreteGrid{D}, quantics) where D + return ntuple(D) do d + grididx = 0 + R_d = g.Rs[d] + + for bitnumber in 1:R_d + site_idx, pos_in_site = g.lookup_table[d][bitnumber] + bit_position = length(g.indextable[site_idx]) - pos_in_site + digit = ((quantics[site_idx] - 1) >> bit_position) & 1 + grididx |= digit << (R_d - bitnumber) + end + grididx + 1 end end -""" -grid index => quantics -""" -grididx_to_quantics(g::Grid{1}, grididx::Int) = grididx_to_quantics(g, (grididx,)) +function _grididx_to_quantics_general!(result::Vector{Int}, g::InherentDiscreteGrid{D}, grididx::NTuple{D,Int}) where D + base = g.base + @inbounds for d in 1:D + zero_based_idx = grididx[d] - 1 + R_d = g.Rs[d] -""" -original coordinate => quantics -""" -function origcoord_to_quantics(g::Grid{d}, coordinate) where {d} - return grididx_to_quantics(g, origcoord_to_grididx(g, _to_tuple(Val(d), coordinate))) -end + for bitnumber in 1:R_d + site_idx, pos_in_site = g.lookup_table[d][bitnumber] + site_length = length(g.indextable[site_idx]) + bit_position = R_d - bitnumber + digit = (zero_based_idx ÷ (base^bit_position)) % base -""" -quantics => grid index -""" -function quantics_to_grididx(g::Grid{d}, bitlist) where {d} - (all(1 .<= bitlist) && all(bitlist .<= digitmax(d, g.base, g.unfoldingscheme))) || - error("Range error for bitlist: $bitlist") - quanticslength(g.R, d, g.unfoldingscheme) == length(bitlist) || - error("Length error for bitlist: $bitlist") - - return _convert_to_scalar_if_possible( - quantics_to_index( - bitlist; - base = g.base, - dims = Val(d), - unfoldingscheme = g.unfoldingscheme, - ), - ) + power = site_length - pos_in_site + result[site_idx] += digit * (base^power) + end + end end +function _grididx_to_quantics_base2!(result::Vector{Int}, g::InherentDiscreteGrid{D}, grididx::NTuple{D,Int}) where D + @inbounds for d in 1:D + zero_based_idx = grididx[d] - 1 + R_d = g.Rs[d] -""" -quantics => original coordinate system -""" -function quantics_to_origcoord(g::Grid{d}, bitlist) where {d} - return _convert_to_scalar_if_possible( - grididx_to_origcoord(g, quantics_to_grididx(g, bitlist)), - ) -end + for bitnumber in 1:R_d + site_idx, pos_in_site = g.lookup_table[d][bitnumber] + site_length = length(g.indextable[site_idx]) + bit_position = R_d - bitnumber + digit = (zero_based_idx >> bit_position) & 1 -""" -Make a wrapper function that takes a bitlist as input -""" -function quanticsfunction(::Type{T}, g::Grid{d}, f::Function)::Function where {T,d} - function _f(bitlist)::T - return f(quantics_to_origcoord(g, bitlist)...) + power = site_length - pos_in_site + result[site_idx] += digit << power + end end - return _f end - -@doc raw""" -The InherentDiscreteGrid struct represents a grid for inherently discrete data. -The grid contains values at specific, -equally spaced points, but these values do not represent discretized versions -of continuous data. Instead, they represent individual data points that are -inherently discrete. -The linear size of the mesh is ``base^R``, where ``base`` defaults to 2. -""" -struct InherentDiscreteGrid{d} <: Grid{d} - R::Int - origin::NTuple{d,Int} - base::Int - unfoldingscheme::Symbol - step::NTuple{d,Int} - - function InherentDiscreteGrid{d}( - R::Int, - origin::Union{NTuple{d,Int},Int}; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - step::Union{NTuple{d,Int},Int} = 1, - ) where {d} - _rangecheck_R(R; base = base) - unfoldingscheme in (:fused, :interleaved) || - error("Invalid unfolding scheme: $unfoldingscheme") - origin_ = origin isa Int ? ntuple(i -> origin, d) : origin - step_ = step isa Int ? ntuple(i -> step, d) : step - new(R, origin_, base, unfoldingscheme, step_) - end - +# ============================================================================ +# Constructors +# ============================================================================ + +function InherentDiscreteGrid{D}( + Rs, + origin=default_origin(Val(D)); + unfoldingscheme=:fused, + step=default_step(Val(D)), + base=2 +) where D + Rs = _to_tuple(Val(D), Rs) + origin = _to_tuple(Val(D), origin) + step = _to_tuple(Val(D), step) + variablenames = ntuple(Symbol, D) + indextable = _build_indextable(variablenames, Rs, unfoldingscheme) + InherentDiscreteGrid{D}(Rs, origin, step, variablenames, base, indextable) end function InherentDiscreteGrid( - R::Int, - origin::Int; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - step::Int = 1, -) - return InherentDiscreteGrid{1}( - R, - origin; - base = base, - unfoldingscheme = unfoldingscheme, - step = step, - ) + variablenames::NTuple{D,Symbol}, + indextable::Vector{Vector{Tuple{Symbol,Int}}}; + origin=default_origin(Val(D)), + step=default_step(Val(D)), + base=2 +) where D + Rs = Tuple(map(variablenames) do variablename + count(index -> first(index) == variablename, Iterators.flatten(indextable)) + end) + + return InherentDiscreteGrid{D}(Rs, origin, step, variablenames, base, indextable) end function InherentDiscreteGrid( - R::Int, - origin::NTuple{N,Int}; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - step::Union{Int,NTuple{N,Int}} = 1, -) where {N} - return InherentDiscreteGrid{N}( - R, - origin; - base = base, - unfoldingscheme = unfoldingscheme, - step = step, - ) + variablenames::NTuple{D,Symbol}, + Rs::NTuple{D,Int}; + origin=default_origin(Val(D)), + step=default_step(Val(D)), + base=2, + unfoldingscheme=:fused +) where {D} + origin = _to_tuple(Val(D), origin) + step = _to_tuple(Val(D), step) + indextable = _build_indextable(variablenames, Rs, unfoldingscheme) + return InherentDiscreteGrid{D}(Rs, origin, step, variablenames, base, indextable) end -""" - grid_min(g::InherentDiscreteGrid) +function InherentDiscreteGrid(Rs::NTuple{D,Int}; variablenames=ntuple(Symbol, D), kwargs...) where {D} + return InherentDiscreteGrid(variablenames, Rs; kwargs...) +end -Returns the grid point with minimal coordinate values. -This is equivalent to [`grid_origin`](@ref). -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_min(grid::InherentDiscreteGrid) = _convert_to_scalar_if_possible(grid.origin) +function InherentDiscreteGrid(R::Int, origin; kwargs...) + return InherentDiscreteGrid{length(origin)}(R, origin; kwargs...) +end +# ============================================================================ +# Basic property accessor functions +# ============================================================================ -""" - grid_max(g::InherentDiscreteGrid) +Base.ndims(::InherentDiscreteGrid{D}) where D = D -Returns the grid point with maximal coordinate values. -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_max(grid::InherentDiscreteGrid) = _convert_to_scalar_if_possible( - grid.origin .+ grid_step(grid) .* (grid.base^grid.R .- 1), -) +Base.length(g::InherentDiscreteGrid) = length(g.indextable) -""" - grid_step(g::InherentDiscreteGrid) +grid_Rs(g::InherentDiscreteGrid) = g.Rs -Returns the distance between adjacent grid points in each dimension. -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_step(grid::InherentDiscreteGrid{d}) where {d} = - _convert_to_scalar_if_possible(grid.step) +grid_indextable(g::InherentDiscreteGrid) = g.indextable +grid_base(g::InherentDiscreteGrid) = g.base +grid_variablenames(g::InherentDiscreteGrid) = g.variablenames -""" - grid_origin(g::InherentDiscreteGrid) +grid_step(g::InherentDiscreteGrid) = _convert_to_scalar_if_possible(g.step) -Returns the origin of the grid, as passed to the constructor during grid creation. -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" grid_origin(g::InherentDiscreteGrid) = _convert_to_scalar_if_possible(g.origin) - -""" -Create a grid for inherently discrete data with origin at 1 -""" -function InherentDiscreteGrid{d}( - R::Int; - base::Integer = 2, - step::Union{NTuple{d,Int},Int} = 1, - unfoldingscheme::Symbol = :fused, -) where {d} - InherentDiscreteGrid{d}( - R, - 1; - base = base, - unfoldingscheme = unfoldingscheme, - step = step, - ) +function sitedim(g::InherentDiscreteGrid, site::Int)::Int + if !(site ∈ eachindex(g.indextable)) + throw(DomainError(site, lazy"Site index out of bounds [1, $(length(g.indextable))].")) + end + return g.base^length(g.indextable[site]) end +# ============================================================================ +# Grid coordinate functions +# ============================================================================ -""" -Create a grid for inherently discrete data with origin at an arbitrary point -""" -#function InherentDiscreteGrid{d}( -#R::Int, origin::Union{NTuple{d,Int},Int}; -#base=2, unfoldingscheme::UnfoldingSchemes.UnfoldingScheme=UnfoldingSchemes.fused, -#) where {d} -#InherentDiscreteGrid{d}( -#R, origin; base=base, unfoldingscheme=unfoldingscheme, step=step -#) -#end - - -""" -Convert a coordinate in the original coordinate system to the corresponding grid index -""" -function origcoord_to_grididx( - g::InherentDiscreteGrid, - coordinate::Union{Int,NTuple{N,Int}}, -) where {N} - return _convert_to_scalar_if_possible( - div.(coordinate .- grid_min(g), grid_step(g)) .+ 1, - ) -end +grid_min(g::InherentDiscreteGrid) = _convert_to_scalar_if_possible(g.origin) +grid_max(g::InherentDiscreteGrid) = _convert_to_scalar_if_possible( + grid_origin(g) .+ grid_step(g) .* (g.base .^ g.Rs .- 1), +) -@doc raw""" -The DiscretizedGrid struct represents a grid for discretized continuous data. -This is used for data that is originally continuous, -but has been discretized for computational purposes. -The grid contains values at specific, equally spaced points, which represent discrete -approximations of the original continuous data. -""" -struct DiscretizedGrid{d} <: Grid{d} - R::Int - lower_bound::NTuple{d,Float64} - upper_bound::NTuple{d,Float64} - base::Int - unfoldingscheme::Symbol - includeendpoint::Bool - - function DiscretizedGrid{d}( - R::Int, - lower_bound, - upper_bound; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - includeendpoint::Bool = false, - ) where {d} - _rangecheck_R(R; base = base) - unfoldingscheme in (:fused, :interleaved) || - error("Invalid unfolding scheme: $unfoldingscheme") - return new( - R, - _to_tuple(Val(d), lower_bound), - _to_tuple(Val(d), upper_bound), - base, - unfoldingscheme, - includeendpoint, - ) - end -end +# ============================================================================ +# Core conversion functions +# ============================================================================ +function quantics_to_grididx(g::InherentDiscreteGrid{D}, quantics::AbstractVector{Int}) where D + # TODO: add switch to turn off input validation + if !(length(quantics) == length(g)) + throw(ArgumentError(lazy"Quantics vector must have length $(length(g.indextable)), got $(length(quantics)).")) + end -""" -Create a discretized grid for a 1D space -""" -function DiscretizedGrid( - R::Int, - grid_min::T, - grid_max::T; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - includeendpoint::Bool = false, -) where {T<:Real} - return DiscretizedGrid{1}( - R, - (grid_min,), - (grid_max,); - base = base, - unfoldingscheme = unfoldingscheme, - includeendpoint = includeendpoint, - ) -end + for site in eachindex(quantics) + if !(1 <= quantics[site] <= sitedim(g, site)) + throw(DomainError(quantics[site], lazy"Quantics value for site $site out of range 1:$(sitedim(g, site)).")) + end + end + result = if g.base == 2 + _quantics_to_grididx_base2(g, quantics) + else + _quantics_to_grididx_general(g, quantics) + end -""" -Create a discretized grid for a d-dimensional space -""" -function DiscretizedGrid( - R::Int, - grid_min::NTuple{d,T}, - grid_max::NTuple{d,T}; - base::Integer = 2, - unfoldingscheme::Symbol = :fused, - includeendpoint::Bool = false, -) where {d,T<:Real} - return DiscretizedGrid{d}( - R, - grid_min, - grid_max; - base = base, - unfoldingscheme = unfoldingscheme, - includeendpoint = includeendpoint, - ) + return _convert_to_scalar_if_possible(result) end +function grididx_to_quantics(g::InherentDiscreteGrid{D}, grididx::Int) where D + grididx_tuple = _to_tuple(Val(D), grididx) + return grididx_to_quantics(g, grididx_tuple) +end +function grididx_to_quantics(g::InherentDiscreteGrid{D}, grididx_tuple::NTuple{D,Int}) where D + # TODO: add switch to turn off input validation + for d in 1:D + if !(1 <= grididx_tuple[d] <= g.maxgrididx[d]) + throw(DomainError(grididx_tuple[d], lazy"Grid index out of bounds [1, $(g.maxgrididx[d])].")) + end + end + result = ones(Int, length(g.indextable)) + if g.base == 2 + _grididx_to_quantics_base2!(result, g, grididx_tuple) + else + _grididx_to_quantics_general!(result, g, grididx_tuple) + end + return result +end -""" - grid_min(g::DiscretizedGrid) - - -Returns the grid point with minimal coordinate values. -This is equivalent to [`lower_bound`](@ref). -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_min(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.lower_bound) - +function grididx_to_origcoord(g::InherentDiscreteGrid{D}, grididx) where D + grididx = _to_tuple(Val(D), grididx) + for d in 1:D + if !(grididx[d] ∈ 1:g.maxgrididx[d]) + throw(DomainError(grididx[d], lazy"Grid index out of bounds [1, $(g.maxgrididx[d])].")) + end + end -""" - grid_max(g::DiscretizedGrid) + res = grid_origin(g) .+ (grididx .- 1) .* grid_step(g) -Returns the grid point with maximal coordinate values. - - If `includeendpoint=false` during creation of the grid, this value is dependent on grid resolution. - - If `includeendpoint=true` during creation of the grid, this function is equivalent to [`upper_bound`](@ref). -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_max(g::DiscretizedGrid) = - g.includeendpoint ? _convert_to_scalar_if_possible(g.upper_bound) : - _convert_to_scalar_if_possible(g.upper_bound .- grid_step(g)) + return _convert_to_scalar_if_possible(res) +end -""" - grid_step(g::DiscretizedGrid) +function origcoord_to_grididx(g::InherentDiscreteGrid{D}, coordinate) where {D} + coord_tuple = _to_tuple(Val(D), coordinate) + bounds_lower = grid_min(g) + bounds_upper = grid_max(g) + # TODO: think about the correct bounds to use here + for d in 1:D + if !(bounds_lower[d] <= coord_tuple[d] <= bounds_upper[d]) + throw(DomainError(coord_tuple[d], lazy"Coordinate out of bounds [$(bounds_lower[d]), $(bounds_upper[d])].")) + end + end -Returns the distance between adjacent grid points in each dimension. -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -grid_step(g::DiscretizedGrid{d}) where {d} = _convert_to_scalar_if_possible( - g.includeendpoint ? (g.upper_bound .- g.lower_bound) ./ (g.base^g.R .- 1) : - (g.upper_bound .- g.lower_bound) ./ (g.base^g.R), -) + discrete_idx = div.(coord_tuple .- grid_min(g), grid_step(g)) .+ 1 + discrete_idx = clamp.(discrete_idx, 1, g.base .^ g.Rs) -""" - upper_bound(g::DiscretizedGrid) - -Returns the upper bound of the grid coordinates, as passed to the constructor during grid creation. - - If `includeendpoint=false` during grid creation, this function returns a point that is one grid spacing beyond the last grid point (which can be obtained through [`grid_max`](@ref)). - - If `includeendpoint=true` during grid creation, this function returns the point with maximal coordinate values. This is equivalent to [`grid_max`](@ref). -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -upper_bound(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.upper_bound) - -""" - lower_bound(g::DiscretizedGrid) - -Returns the grid point with minimal coordinate values, as passed to the constructor during grid creation. -This is equivalent to [`grid_min`](@ref). -The return value is scalar for 1D grids, and a `Tuple` otherwise. -""" -lower_bound(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.lower_bound) - -function DiscretizedGrid{d}(R::Int; base = 2, unfoldingscheme::Symbol = :fused) where {d} - return DiscretizedGrid{d}( - R, - ntuple(i -> 0.0, d), - ntuple(i -> 1.0, d); - base = base, - unfoldingscheme = unfoldingscheme, - ) + return _convert_to_scalar_if_possible(discrete_idx) end +function origcoord_to_quantics(g::InherentDiscreteGrid, coordinate) + grididx_to_quantics(g, origcoord_to_grididx(g, coordinate)) +end -""" -Convert a coordinate in the original coordinate system to the corresponding grid index -""" -function origcoord_to_grididx(g::DiscretizedGrid, coordinate::NTuple{N,Float64}) where {N} - all(lower_bound(g) .<= coordinate .<= upper_bound(g)) || error( - "Bound Error: $(coordinate), lower_bound=$(lower_bound(g)), upper_bound=$(upper_bound(g))", - ) - clip(x) = max.(1, min.(x, g.base^g.R)) - return _convert_to_scalar_if_possible( - ((coordinate .- grid_min(g)) ./ grid_step(g) .+ 1) .|> round .|> Int, - ) +function quantics_to_origcoord(g::InherentDiscreteGrid, quantics) + grididx_to_origcoord(g, quantics_to_grididx(g, quantics)) end -function origcoord_to_grididx(g::DiscretizedGrid{1}, coordinate::Float64) - origcoord_to_grididx(g, (coordinate,))[1] +# ============================================================================ +# Other utility functions +# ============================================================================ + +function localdimensions(g::InherentDiscreteGrid)::Vector{Int} + return g.base .^ length.(g.indextable) end diff --git a/src/grid_discretized.jl b/src/grid_discretized.jl new file mode 100644 index 0000000..c104d81 --- /dev/null +++ b/src/grid_discretized.jl @@ -0,0 +1,476 @@ +""" + DiscretizedGrid{D} + +A discretized grid structure for D-dimensional grids with variable resolution, +supporting efficient conversion between quantics, grid indices, and original coordinates. +A `DiscretizedGrid` instance is intended to undergird a quantics tensor train +with a specific index structure, as defined in the `indextable` field. +For example, say indextable is `[[(:a, 1), (:b, 2)], [(:a, 2)], [(:b, 1), (:a, 3)]]`, +then the corresponding tensor train has 3 tensor cores: + + a_1 b_2 a_2 b_1 a_3 + | | | | | + ┌──┴───┴──┐ ┌────┴────┐ ┌──┴───┴──┐ + │ │ │ │ │ │ + │ │────│ │────│ │ + │ │ │ │ │ │ + └─────────┘ └─────────┘ └─────────┘ + +This object may be constructed with +```julia-repl +julia> grid = DiscretizedGrid((:a, :b), [[(:a, 1), (:b, 2)], [(:a, 2)], [(:b, 1), (:a, 3)]]) +DiscretizedGrid{2} with 8×4 = 32 grid points +├─ Variables: (a, b) +├─ Resolutions: (a: 3, b: 2) +├─ Domain: unit square [0, 1)² +├─ Grid spacing: (Δa = 0.125, Δb = 0.25) +└─ Tensor train: 3 sites (dimensions: 4-2-4) +``` +and represents a 2^3 x 2^2 discretization of the unit square in the 2D plane (the x mark grid points): + + 1.0 ┌───────────────────────────────┐ + │ │ + │ │ + 0.75 x x x x x x x x │ + │ │ + │ │ + b 0.5 x x x x x x x x │ + │ │ + │ │ + 0.25 x x x x x x x x │ + │ │ + │ │ + 0.0 x───x───x───x───x───x───x───x───┘ + 0.0 0.25 0.5 0.75 1.0 + + a + +If something other than a unit square is desired, `lower_bound` and `upper_bound` +can be specified. Also, bases different than the default base 2 can be used, +e.g. `base=3` for a ternary grid. + +In addition to the plain methods, there is a convenience layer for conversion +from the original coordinates +```julia-repl +julia> origcoord_to_grididx(grid; a=0.5, b=0.25) +(5, 2) + +julia> origcoord_to_quantics(grid; a=0.5, b=0.25) +3-element Vector{Int64}: + 4 + 1 + 1 +``` +and also from grid indices +```julia-repl +julia> grididx_to_origcoord(grid; a=5, b=2) +(0.5, 0.25) + +julia> grididx_to_quantics(grid; a=5, b=2) +3-element Vector{Int64}: + 4 + 1 + 1 +``` + +For a simpler grid, we can just supply the resolution in each dimension: +```julia-repl +julia> boring_grid = DiscretizedGrid((3, 9)) +DiscretizedGrid{2} with 8×512 = 4096 grid points +├─ Resolutions: (1: 3, 2: 9) +├─ Domain: unit square [0, 1)² +├─ Grid spacing: (Δ1 = 0.125, Δ2 = 0.001953125) +└─ Tensor train: 12 sites (uniform dimension 2) +``` +In this case, variable names are automatically generated as `1`, `2`, etc. +""" +struct DiscretizedGrid{D} <: Grid{D} + discretegrid::InherentDiscreteGrid{D} + lower_bound::NTuple{D,Float64} + upper_bound::NTuple{D,Float64} + + function DiscretizedGrid{D}( + Rs, lower_bound, upper_bound, variablenames, base, indextable, includeendpoint + ) where {D} + lower_bound = _to_tuple(Val(D), lower_bound) + upper_bound = _to_tuple(Val(D), upper_bound) + for d in 1:D + if !(lower_bound[d] < upper_bound[d]) + throw(ArgumentError(lazy"Got (lower_bound[$d], upper_bound[$d]) = $((lower_bound[d], upper_bound[d])). Each lower bound needs to be strictly less than the corresponding upper bound.")) + end + end + + discretegrid = InherentDiscreteGrid(variablenames, indextable; base) + + upper_bound = _adjust_upper_bounds( + upper_bound, lower_bound, includeendpoint, base, Rs, Val(D) + ) + + return new{D}(discretegrid, lower_bound, upper_bound) + end +end + +# ============================================================================ +# Helper/utility functions +# ============================================================================ + +function _adjust_upper_bounds(upper_bound, lower_bound, includeendpoint, base, Rs, ::Val{D}) where D + includeendpoint = _to_tuple(Val(D), includeendpoint) + for d in 1:D + if iszero(Rs[d]) && includeendpoint[d] + throw(ArgumentError(lazy"Got Rs[$d] = 0 and includeendpoint[$d] = true. This is not allowed.")) + end + end + + return ntuple(D) do d + if includeendpoint[d] + upper_bound[d] + (upper_bound[d] - lower_bound[d]) / (base^Rs[d] - 1) + else + upper_bound[d] + end + end +end + +default_lower_bound(::Val{D}) where D = _to_tuple(Val(D), 0.0) + +default_upper_bound(::Val{D}) where D = _to_tuple(Val(D), 1.0) + +function _handle_kwargs_input(g::DiscretizedGrid{D}; kwargs...) where {D} + provided_keys = keys(kwargs) + expected_keys = grid_variablenames(g) + if !(Set(provided_keys) == Set(expected_keys)) + throw(ArgumentError(lazy"Expected keyword arguments $(expected_keys), got $(tuple(provided_keys...))")) + end + + if !all(v -> v isa Real, values(kwargs)) + throw(ArgumentError(lazy"Got kwargs = $kwargs. All keyword argument values must be Real numbers")) + end + + return ntuple(D) do d + variablename = grid_variablenames(g)[d] + kwargs[variablename] + end +end + +# ============================================================================ +# Constructors +# ============================================================================ + +function DiscretizedGrid( + variablenames::NTuple{D,Symbol}, + Rs::NTuple{D,Int}; + lower_bound=default_lower_bound(Val(D)), + upper_bound=default_upper_bound(Val(D)), + base::Int=2, + unfoldingscheme::Symbol=:fused, + includeendpoint=false +) where {D} + indextable = _build_indextable(variablenames, Rs, unfoldingscheme) + + return DiscretizedGrid{D}(Rs, lower_bound, upper_bound, variablenames, base, indextable, includeendpoint) +end + +function DiscretizedGrid(Rs::NTuple{D,Int}; variablenames=ntuple(Symbol, D), kwargs...) where {D} + return DiscretizedGrid(variablenames, Rs; kwargs...) +end + +function DiscretizedGrid{D}( + R, + lower_bound=default_lower_bound(Val(D)), + upper_bound=default_upper_bound(Val(D)); + kwargs... +) where {D} + return DiscretizedGrid(_to_tuple(Val(D), R); lower_bound, upper_bound, kwargs...) +end + +function DiscretizedGrid( + R, + lower_bound::NTuple{D,Real}, + upper_bound::NTuple{D,Real}; + kwargs... +) where {D} + return DiscretizedGrid(_to_tuple(Val(D), R); lower_bound, upper_bound, kwargs...) +end + +function DiscretizedGrid( + variablenames::NTuple{D,Symbol}, + indextable::Vector{Vector{Tuple{Symbol,Int}}}; + lower_bound=default_lower_bound(Val(D)), + upper_bound=default_upper_bound(Val(D)), + base::Int=2, + includeendpoint=false +) where D + Rs = Tuple(map(variablenames) do variablename + count(index -> first(index) == variablename, Iterators.flatten(indextable)) + end) + + return DiscretizedGrid{D}(Rs, lower_bound, upper_bound, variablenames, base, indextable, includeendpoint) +end + +function DiscretizedGrid(R::Int, lower_bound::Real, upper_bound::Real; kwargs...) + return DiscretizedGrid{1}(R, lower_bound, upper_bound; kwargs...) +end + +# ============================================================================ +# Basic property accessor functions +# ============================================================================ + +Base.ndims(::DiscretizedGrid{D}) where D = D + +Base.length(g::DiscretizedGrid) = length(grid_indextable(g)) + +grid_Rs(g::DiscretizedGrid{D}) where D = grid_Rs(g.discretegrid) + +grid_indextable(g::DiscretizedGrid{D}) where D = grid_indextable(g.discretegrid) + +grid_base(g::DiscretizedGrid{D}) where D = grid_base(g.discretegrid) + +grid_variablenames(g::DiscretizedGrid{D}) where D = grid_variablenames(g.discretegrid) + +upper_bound(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.upper_bound) + +lower_bound(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.lower_bound) + +grid_origin(g::DiscretizedGrid) = lower_bound(g) + +function sitedim(g::DiscretizedGrid, site::Int)::Int + if !(site ∈ eachindex(grid_indextable(g))) + throw(DomainError(site, lazy"Site index out of bounds [1, $(length(grid_indextable(g)))]")) + end + return grid_base(g)^length(grid_indextable(g)[site]) +end + +# ============================================================================ +# Grid coordinate functions +# ============================================================================ + +grid_min(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.lower_bound) + +grid_max(g::DiscretizedGrid) = _convert_to_scalar_if_possible(g.upper_bound .- grid_step(g)) + +grid_step(g::DiscretizedGrid) = _convert_to_scalar_if_possible( + (upper_bound(g) .- lower_bound(g)) ./ (grid_base(g) .^ grid_Rs(g)), +) + +function grid_origcoords(g::DiscretizedGrid, d::Int) + if !(1 ≤ d ≤ ndims(g)) + throw(DomainError(d, lazy"Dimension $d out of bounds [1, $(ndims(g))].")) + end + start = grid_min(g)[d] + stop = grid_max(g)[d] + length = grid_base(g)^grid_Rs(g)[d] + return range(start, stop, length) +end + +function grid_origcoords(g::DiscretizedGrid, variablename::Symbol) + d = findfirst(==(variablename), grid_variablenames(g)) + isnothing(d) && throw(ArgumentError(lazy"Variable name :$variablename not found in grid. Available variables: $(grid_variablenames(g))")) + return grid_origcoords(g, d) +end + +# ============================================================================ +# Core conversion functions +# ============================================================================ + +function quantics_to_grididx(g::DiscretizedGrid, quantics::AbstractVector{Int}) + return quantics_to_grididx(g.discretegrid, quantics) +end + +function grididx_to_quantics(g::DiscretizedGrid, grididx) + return grididx_to_quantics(g.discretegrid, grididx) +end + +function grididx_to_origcoord(g::DiscretizedGrid{D}, index) where {D} + index_tuple = _to_tuple(Val(D), index) + for d in 1:D + if !(1 <= index_tuple[d] <= grid_base(g)^grid_Rs(g)[d]) + throw(DomainError(index_tuple[d], lazy"Grid index out of bounds [1, $(grid_base(g) ^ grid_Rs(g)[d])].")) + end + end + + res = ntuple(D) do d + step_d = (g.upper_bound[d] - g.lower_bound[d]) / (grid_base(g)^grid_Rs(g)[d]) + g.lower_bound[d] + (index_tuple[d] - 1) * step_d + end + + return _convert_to_scalar_if_possible(res) +end + +function origcoord_to_grididx(g::DiscretizedGrid{D}, coordinate) where {D} + coord_tuple = _to_tuple(Val(D), coordinate) + + bounds_lower = lower_bound(g) + bounds_upper = upper_bound(g) + for d in 1:D + if !(bounds_lower[d] <= coord_tuple[d] <= bounds_upper[d]) + throw(DomainError(coord_tuple[d], lazy"Coordinate out of bounds [$(bounds_lower[d]), $(bounds_upper[d])] (dimension $d).")) + end + end + + steps = grid_step(g) + indices = ntuple(D) do d + target = coord_tuple[d] + step_d = steps[d] + + continuous_idx = (target - bounds_lower[d]) / step_d + 1 + + discrete_idx = round(Int, continuous_idx) + clamp(discrete_idx, 1, grid_base(g)^grid_Rs(g)[d]) + end + + return _convert_to_scalar_if_possible(indices) +end + +function origcoord_to_quantics(g::DiscretizedGrid{D}, coordinate) where {D} + coord_tuple = _to_tuple(Val(D), coordinate) + grid_idx = origcoord_to_grididx(g, coord_tuple) + return grididx_to_quantics(g, grid_idx) +end + +function quantics_to_origcoord(g::DiscretizedGrid{D}, quantics::AbstractVector{Int}) where {D} + grid_idx = quantics_to_grididx(g, quantics) + return grididx_to_origcoord(g, grid_idx) +end + +# ============================================================================ +# Keyword argument convenience functions +# ============================================================================ + +function origcoord_to_grididx(g::DiscretizedGrid; kwargs...) + coordinate = _handle_kwargs_input(g; kwargs...) + return origcoord_to_grididx(g, coordinate) +end + +function origcoord_to_quantics(g::DiscretizedGrid; kwargs...) + coordinate = _handle_kwargs_input(g; kwargs...) + return origcoord_to_quantics(g, coordinate) +end + +function grididx_to_origcoord(g::DiscretizedGrid; kwargs...) + index = _handle_kwargs_input(g; kwargs...) + return grididx_to_origcoord(g, index) +end + +function grididx_to_quantics(g::DiscretizedGrid; kwargs...) + index = _handle_kwargs_input(g; kwargs...) + return grididx_to_quantics(g, index) +end + +# ============================================================================ +# Other utility functions +# ============================================================================ + +function localdimensions(g::DiscretizedGrid)::Vector{Int} + return grid_base(g) .^ length.(grid_indextable(g)) +end + +function quanticsfunction(::Type{T}, g::DiscretizedGrid, f::F)::Function where {T,F<:Function} + function wrapped_function(quantics)::T + coords = quantics_to_origcoord(g, quantics) + if coords isa Tuple + return f(coords...) + else + return f(coords) + end + end + return wrapped_function +end + +# ============================================================================ +# Display/show methods +# ============================================================================ + +function Base.show(io::IO, ::MIME"text/plain", g::DiscretizedGrid{D}) where D + print(io, "DiscretizedGrid{$D}") + + # Grid resolution and total points + total_points = prod(big(grid_base(g)) .^ grid_Rs(g)) + if D <= 1 + print(io, " with $total_points grid point" * (isone(total_points) ? "" : "s")) + else + print(io, " with $(join(big(grid_base(g)) .^ grid_Rs(g), " × ")) = $total_points grid points") + end + + # Variable names (if meaningful) + if any(name -> !startswith(string(name), r"^\d+$"), grid_variablenames(g)) + var_str = join(grid_variablenames(g), ", ") + print(io, "\n├─ Variables: ($var_str)") + end + + # Resolution per dimension + if D == 1 + print(io, "\n├─ Resolution: $(grid_Rs(g)[1]) bits") + else + res_str = join(["$(grid_variablenames(g)[i]): $(grid_Rs(g)[i])" for i in 1:D], ", ") + print(io, "\n├─ Resolutions: ($res_str)") + end + + # Bounds (only show if not default unit interval/square/cube) + default_lower = default_lower_bound(Val(D)) + default_upper = default_upper_bound(Val(D)) + if lower_bound(g) != default_lower || any(abs.(upper_bound(g) .- default_upper) .> 1e-10) + if D == 1 + print(io, "\n├─ Domain: [$(lower_bound(g)[1]), $(upper_bound(g)[1]))") + else + bounds_str = join(["[$(lower_bound(g)[i]), $(upper_bound(g)[i]))" for i in 1:D], " × ") + print(io, "\n├─ Domain: $bounds_str") + end + + # Grid spacing + step_vals = grid_step(g) + if D == 1 + print(io, "\n├─ Grid spacing: $(step_vals)") + else + step_str = join(["Δ$(grid_variablenames(g)[i]) = $(step_vals[i])" for i in 1:D], ", ") + print(io, "\n├─ Grid spacing: ($step_str)") + end + else + # For unit domain, show appropriate unit description + unit_domain_str = if D == 1 + "unit interval [0, 1)" + elseif D == 2 + "unit square [0, 1)²" + elseif D == 3 + "unit cube [0, 1)³" + else + "unit hypercube [0, 1)^$D" + end + print(io, "\n├─ Domain: $unit_domain_str") + + # Grid spacing + step_vals = grid_step(g) + if D == 1 + print(io, "\n├─ Grid spacing: $(step_vals)") + else + step_str = join(["Δ$(grid_variablenames(g)[i]) = $(step_vals[i])" for i in 1:D], ", ") + print(io, "\n├─ Grid spacing: ($step_str)") + end + end + + # Base (only show if not binary) + if grid_base(g) != 2 + print(io, "\n├─ Base: $(grid_base(g))") + end + + # Tensor structure summary + num_sites = length(grid_indextable(g)) + sitedims = Int[sitedim(g, site) for site in 1:num_sites] + + print(io, "\n└─ Tensor train: $num_sites sites") + if !isempty(sitedims) + if allequal(sitedims) + print(io, " (uniform dimension $(sitedims[1]))") + else + print(io, " (dimensions: $(join(sitedims, "-")))") + end + end +end + +function Base.show(io::IO, g::DiscretizedGrid{D}) where D + total_points = prod(grid_base(g) .^ grid_Rs(g)) + if D == 1 + print(io, "DiscretizedGrid{$D}($(grid_base(g)^grid_Rs(g)[1]) points)") + else + print(io, "DiscretizedGrid{$D}($(join(grid_base(g) .^ grid_Rs(g), "×")) points)") + end +end diff --git a/src/quantics.jl b/src/quantics.jl deleted file mode 100644 index fcafd1c..0000000 --- a/src/quantics.jl +++ /dev/null @@ -1,268 +0,0 @@ -function _checkdigits(base::Integer, digitlist) - maximum(digitlist) <= base || error("maximum(digitlist) <= base") - minimum(digitlist) >= 0 || error("minimum(digitlist) >= 0") -end - -""" - fuse_dimensions(digitlists...) - -Fuse d digitlists that represent a quantics index into a digitlist where each bit -has dimension base^d. This fuses legs for different dimensions that have equal length -scale (see QTCI paper). - -Inverse of [`unfuse_dimensions`](@ref). -""" -function fuse_dimensions(digitlists...; base::Integer = 2) - _checkdigits.(base, digitlists) - result = ones(Int, length(digitlists[1])) - return fuse_dimensions!(result, digitlists...; base = base) -end - - - -function fuse_dimensions!(fused::AbstractArray{<:Integer}, digitlists...; base::Integer = 2) - p = 1 - fused .= 1 - for d in eachindex(digitlists) - @. fused += (digitlists[d] - 1) * p - p *= base - end - return fused -end - - -""" - unfuse_dimensions([base=Val(2)], digitlist, d) - -Unfuse up a fused digitlist with bits of dimension base^d into d digitlists where each bit has dimension `base`. -Inverse of [`fuse_dimensions`](@ref). -""" -function unfuse_dimensions(digitlist, d; base::Integer = 2) - result = [zeros(Int, length(digitlist)) for _ = 1:d] - return unfuse_dimensions!(result, digitlist; base = base) -end - -function unfuse_dimensions!(digitlists, digitlist; base::Integer = 2) - ndim = length(digitlists) - R = length(digitlist) - for i = 1:ndim - for j = 1:R - digitlists[i][j] = - _digit_at_index(digitlist[j], ndim - i + 1; numdigits = ndim, base = base) - end - end - return digitlists -end - - -""" - interleave_dimensions(digitlists...) - -Interleaves the indices of all digitlists into one long digitlist. Use this for -quantics representation of multidimensional objects without fusing indices. -Inverse of [`deinterleave_dimensions`](@ref). -""" -function interleave_dimensions(digitlists...)::Vector{Int} - results = Vector{Int}(undef, length(digitlists[1]) * length(digitlists)) - return interleave_dimensions!(results, digitlists...) -end - -function interleave_dimensions!( - interleaved_digitlist::AbstractArray{<:Integer}, - digitlists..., -) - idx = 1 - for i in eachindex(digitlists[1]) - for d in eachindex(digitlists) - interleaved_digitlist[idx] = digitlists[d][i] - idx += 1 - end - end - return interleaved_digitlist -end - -""" - deinterleave_dimensions(digitlist, d) - -Reverses the interleaving of bits, i.e. yields digitlists for each dimension from -a long interleaved digitlist. Inverse of [`interleave_dimensions`](@ref). -""" -function deinterleave_dimensions(digitlist, d) - return [digitlist[i:d:end] for i = 1:d] -end - - -function deinterleave_dimensions!( - deinterleaved_digitlists::AbstractArray{<:AbstractArray{I}}, - digitlist, -) where {I<:Integer} - d = length(deinterleaved_digitlists) - for i = 1:d - @. deinterleaved_digitlists[i] = digitlist[i:d:end] - end - return deinterleaved_digitlists -end - -""" -Convert a fused quantics representation to an unfused quantics representation -""" -function fused_to_interleaved( - digitlist::AbstractVector{T}, - d; - base = 2, -)::Vector{T} where {T<:Integer} - return interleave_dimensions(unfuse_dimensions(digitlist, d; base = base)...) -end - -""" -Convert an unfused quantics representation to an fused quantics representation -""" -function interleaved_to_fused( - digitlist::AbstractVector{T}; - base::Integer = 2, - d::Int = 1, -)::Vector{T} where {T<:Integer} - fused_digitlist = Vector{T}(undef, length(digitlist) ÷ d) - fuse_dimensions!(fused_digitlist, deinterleave_dimensions(digitlist, d)...) - return fused_digitlist -end - -""" - function quantics_to_index_fused( - digitlist::AbstractVector{<:Integer}; - base::Val{B}=Val(2), dims::Val{d}=Val(1) - )::NTuple{d,Int} where {B, d} - -Convert a d-dimensional index from fused quantics representation to d Integers. - -* `base` base for quantics (default: 2) -* `d` number of dimensions -* `digitlist` base-b representation - -See also [`quantics_to_index`](https://tensors4fields.gitlab.io/quanticstci.jl/dev/#QuanticsTCI.quantics_to_index). -""" -function quantics_to_index_fused( - digitlist::AbstractVector{<:Integer}; - base::Integer = 2, - dims::Val{d} = Val(1), -)::NTuple{d,Int} where {d} - R = length(digitlist) - result = ones(MVector{d,Int}) - - maximum(digitlist) <= base^d || error("maximum(digitlist) <= base^d") - minimum(digitlist) >= 0 || error("minimum(digitlist) >= 0") - - for n = 1:R # from the least to most significant digit - scale = base^(n - 1) # length scale - tmp = digitlist[R-n+1] - 1 - for i = 1:d # in the order of 1st dim, 2nd dim, ... - div_, rem_ = divrem(tmp, base) - result[i] += rem_ * scale - tmp = div_ - end - end - - return tuple(result...) -end - -""" - function quantics_to_index( - digitlist::AbstractVector{<:Integer}; - base::Val{B}=Val(2), dims::Val{d}=Val(1), - unfoldingscheme::Symbol=:fused - )::NTuple{d,Int} where {B, d} - -Convert a d-dimensional index from fused quantics representation to d Integers. - -* `base` base for quantics (default: 2) -* `d` number of dimensions -* `digitlist` base-b representation -* `unfoldingscheme` Unfolding scheme (fused or interleaved) -""" -function quantics_to_index( - digitlist::AbstractVector{<:Integer}; - base::Integer = 2, - dims::Val{d} = Val(1), - unfoldingscheme::Symbol = :fused, -)::NTuple{d,Int} where {d} - if unfoldingscheme === :fused - return quantics_to_index_fused(digitlist, base = base, dims = dims) - else - return quantics_to_index_fused( - interleaved_to_fused(digitlist; base = base, d = d), - base = base, - dims = dims, - ) - end -end - - -""" - function index_to_quantics!(digitlist, index::Integer; base::Integer=2) - -* `digitlist` base-b representation (1d vector) -* `base` base for quantics (default: 2) -""" -function index_to_quantics!(digitlist, index::Integer; base::Integer = 2) - numdigits = length(digitlist) - for i = 1:numdigits - digitlist[i] = mod(index - 1, base^(numdigits - i + 1)) ÷ base^(numdigits - i) + 1 - end - return digitlist -end - -""" - index_to_quantics(index::Integer; numdigits=8, base::Integer=2) - -Does the same as [`index_to_quantics!`](@ref) but returns a new vector. -""" -function index_to_quantics(index::Integer; numdigits = 8, base::Integer = 2) - digitlist = Vector{Int}(undef, numdigits) - return index_to_quantics!(digitlist, index; base = base) -end - -""" - index_to_quantics_fused!(digitlist::AbstractVector{<:Integer}, index::NTuple{D,<:Integer}; base::Integer=2) where {D} - -Does the opposite of [`quantics_to_index_fused`](@ref) -""" -function index_to_quantics_fused!( - digitlist::AbstractVector{<:Integer}, - index::NTuple{D,<:Integer}; - base::Integer = 2, -) where {D} - R = length(digitlist) - ndims = length(index) - digitlist .= 1 - for dim = 1:ndims - for i = 1:R # from the left to right - digitlist[i] += - (base^(dim - 1)) * - (_digit_at_index(index[dim], i; numdigits = R, base = base) - 1) - end - end - return digitlist -end - -function index_to_quantics_fused( - index::NTuple{D,<:Integer}; - numdigits::Integer = 8, - base::Integer = 2, -) where {D} - digitlist = Vector{Int}(undef, numdigits) - return index_to_quantics_fused!(digitlist, index; base = base) -end - - -""" -`index`, `position` and the result are one-based. - -`base`: base -`index`: The integer to look at. -`position=1`: Specify the position of the digit to look at. 1 is the most significant (most left) digit. -`numdigits=8`: Specify the number of digits in the number `index`. -""" -function _digit_at_index(index, position; numdigits = 8, base::Integer = 2) - p_ = numdigits - position + 1 - return mod((index - 1), base^p_) ÷ base^(p_ - 1) + 1 -end diff --git a/test/code_quality_test.jl b/test/code_quality_test.jl index 5d2a488..7eaa138 100644 --- a/test/code_quality_test.jl +++ b/test/code_quality_test.jl @@ -1,23 +1,13 @@ @testitem "Code quality (Aqua.jl)" begin using Aqua - import QuanticsGrids - - @testset "Code quality (Aqua.jl)" begin - Aqua.test_all(QuanticsGrids; unbound_args = false, deps_compat = false) - end - + Aqua.test_all(QuanticsGrids; deps_compat=false) end @testitem "Code linting (JET.jl)" begin - using Test using JET - import QuanticsGrids - if VERSION >= v"1.10" - @testset "Code linting (JET.jl)" begin - JET.test_package(QuanticsGrids; target_defined_modules = true) - end + JET.test_package(QuanticsGrids; target_defined_modules=true) end end diff --git a/test/comprehensive_tests.jl b/test/comprehensive_tests.jl new file mode 100644 index 0000000..d67ea7f --- /dev/null +++ b/test/comprehensive_tests.jl @@ -0,0 +1,404 @@ +@testitem "DiscretizedGrid comprehensive functionality test" begin + variablenames = (:x, :y, :z, :w) + complex_tupletable = [ + [(:x, 1), (:y, 1), (:z, 1)], + [(:w, 1)], + [(:z, 2), (:w, 2), (:x, 3), (:y, 3)], + [(:x, 2), (:y, 2)], + [(:z, 3)], + [(:w, 3), (:x, 4)], + [(:y, 4), (:z, 4), (:w, 4)] + ] + + lower_bound = (-π, 1e-12, 1e8, -1e6) + upper_bound = (2 * π, 1e-10, 1e8 + 1000, 1e6) + base = 3 + + grid = QuanticsGrids.DiscretizedGrid(variablenames, complex_tupletable; + lower_bound=lower_bound, upper_bound=upper_bound, base=base, includeendpoint=true) + + @test ndims(grid) == 4 + @test length(grid) == 7 + @test QuanticsGrids.grid_base(grid) == 3 + @test QuanticsGrids.grid_variablenames(grid) == variablenames + + expected_Rs = (4, 4, 4, 4) + @test QuanticsGrids.grid_Rs(grid) == expected_Rs + + expected_localdims = [27, 3, 81, 9, 3, 9, 27] + @test QuanticsGrids.localdimensions(grid) == expected_localdims + + for (i, expected_dim) in enumerate(expected_localdims) + @test QuanticsGrids.sitedim(grid, i) == expected_dim + end + + lb = QuanticsGrids.lower_bound(grid) + ub = QuanticsGrids.upper_bound(grid) + gmin = QuanticsGrids.grid_min(grid) + gmax = QuanticsGrids.grid_max(grid) + gstep = QuanticsGrids.grid_step(grid) + gorigin = QuanticsGrids.grid_origin(grid) + + @test lb == lower_bound + @test gmin == lower_bound + @test gorigin == lower_bound + + @test ub != upper_bound + @test all(ub .> upper_bound) + + @test all(gmax .≈ upper_bound) + + expected_steps = (ub .- lb) ./ (base .^ QuanticsGrids.grid_Rs(grid)) + @test all(gstep .≈ expected_steps) + + corner_indices = [ + (1, 1, 1, 1), + (base^4, base^4, base^4, base^4), + (1, base^4, 1, base^4), + (base^4, 1, base^4, 1), + (base^2, base^2, base^2, base^2), + ] + + for grididx in corner_indices + quantics = grididx_to_quantics(grid, grididx) + recovered_grididx = quantics_to_grididx(grid, quantics) + @test recovered_grididx == grididx + + origcoord = grididx_to_origcoord(grid, grididx) + recovered_grididx2 = QuanticsGrids.origcoord_to_grididx(grid, origcoord) + @test recovered_grididx2 == grididx + + origcoord2 = quantics_to_origcoord(grid, quantics) + recovered_quantics = origcoord_to_quantics(grid, origcoord2) + @test recovered_quantics == quantics + + @test all(origcoord .≈ origcoord2) + + @test all(lb .<= origcoord .<= ub) + end + + for _ in 1:50 + random_quantics = [rand(1:QuanticsGrids.sitedim(grid, site)) for site in 1:length(grid)] + + grididx = quantics_to_grididx(grid, random_quantics) + origcoord = grididx_to_origcoord(grid, grididx) + back_to_grididx = QuanticsGrids.origcoord_to_grididx(grid, origcoord) + back_to_quantics = grididx_to_quantics(grid, back_to_grididx) + back_to_origcoord = quantics_to_origcoord(grid, back_to_quantics) + + @test back_to_grididx == grididx + @test back_to_quantics == random_quantics + @test all(back_to_origcoord .≈ origcoord) + + @test all(1 .<= grididx .<= (base .^ QuanticsGrids.grid_Rs(grid))) + @test all(lb .<= origcoord .<= ub) + @test length(back_to_quantics) == length(grid) + end + + boundary_coords = [ + lb, + ub, + upper_bound, + gmax, + ] + + for coord in boundary_coords + if all(lb .<= coord .<= ub) + grididx = QuanticsGrids.origcoord_to_grididx(grid, coord) + recovered_coord = grididx_to_origcoord(grid, grididx) + re_recovered_grididx = QuanticsGrids.origcoord_to_grididx(grid, recovered_coord) + @test re_recovered_grididx == grididx + end + end + + epsilon = 1e-10 + for i in 1:4 + coord_inside = ntuple(j -> j == i ? lb[j] + epsilon : (lb[j] + ub[j]) / 2, 4) + if all(lb .<= coord_inside .<= ub) + grididx = QuanticsGrids.origcoord_to_grididx(grid, coord_inside) + @test all(1 .<= grididx .<= (base .^ QuanticsGrids.grid_Rs(grid))) + end + end + + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (lb[1] - 1, lb[2], lb[3], lb[4])) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (ub[1] + 1, ub[2], ub[3], ub[4])) + + @test_throws DomainError grididx_to_origcoord(grid, (0, 1, 1, 1)) + @test_throws DomainError grididx_to_origcoord(grid, (base^4 + 1, 1, 1, 1)) + + invalid_quantics = [28, 1, 1, 1, 1, 1, 1] + @test_throws DomainError quantics_to_grididx(grid, invalid_quantics) + + precision_test_coords = [ + (lb[1] + gstep[1] * 0.5, lb[2], lb[3], lb[4]), + (lb[1], lb[2] + gstep[2] * 0.999999, lb[3], lb[4]), + (ub[1] - gstep[1] * 1e-15, ub[2], ub[3], ub[4]), + ] + + for coord in precision_test_coords + if all(lb .<= coord .<= ub) + grididx = QuanticsGrids.origcoord_to_grididx(grid, coord) + recovered_coord = grididx_to_origcoord(grid, grididx) + re_recovered_grididx = QuanticsGrids.origcoord_to_grididx(grid, recovered_coord) + @test re_recovered_grididx == grididx + end + end + + test_function(x, y, z, w) = sin(x) * cos(y) * exp(-z / 1e8) * (w^2) + quantics_func = QuanticsGrids.quanticsfunction(Float64, grid, test_function) + + for _ in 1:10 + random_quantics = [rand(1:QuanticsGrids.sitedim(grid, site)) for site in 1:length(grid)] + + result1 = quantics_func(random_quantics) + + origcoord = quantics_to_origcoord(grid, random_quantics) + result2 = test_function(origcoord...) + + @test result1 ≈ result2 + end + + simple_grid = DiscretizedGrid((3, 4); base=5, lower_bound=(0.0, -1.0), upper_bound=(1.0, 1.0), unfoldingscheme=:interleaved) + @test QuanticsGrids.grid_Rs(simple_grid) == (3, 4) + @test QuanticsGrids.grid_base(simple_grid) == 5 + @test QuanticsGrids.localdimensions(simple_grid) == [5, 5, 5, 5, 5, 5, 5] + + grid_1d = DiscretizedGrid{1}(8; base=2, lower_bound=(-5.0,), upper_bound=(10.0,)) + @test ndims(grid_1d) == 1 + @test QuanticsGrids.grid_Rs(grid_1d) == (8,) + + for i in [1, 100, 256] + coord = grididx_to_origcoord(grid_1d, i) + back_idx = QuanticsGrids.origcoord_to_grididx(grid_1d, coord) + @test back_idx == i + end + + Rs_test = (3, 3) + grid_fused = DiscretizedGrid(Rs_test; unfoldingscheme=:fused, base=2) + grid_interleaved = DiscretizedGrid(Rs_test; unfoldingscheme=:interleaved, base=2) + + @test length(grid_fused) == 3 + @test length(grid_interleaved) == 6 + + @test QuanticsGrids.grid_Rs(grid_fused) == QuanticsGrids.grid_Rs(grid_interleaved) + + test_grididx = (4, 6) + coord_fused = grididx_to_origcoord(grid_fused, test_grididx) + coord_interleaved = grididx_to_origcoord(grid_interleaved, test_grididx) + @test all(coord_fused .≈ coord_interleaved) + + large_grid = DiscretizedGrid{2}(20; lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0)) + + max_idx = 2^20 + extreme_indices = [1, max_idx ÷ 2, max_idx - 1, max_idx] + for idx in extreme_indices + grididx = (idx, idx) + if all(1 .<= grididx .<= (2, 2) .^ QuanticsGrids.grid_Rs(large_grid)) + quantics = grididx_to_quantics(large_grid, grididx) + recovered_grididx = quantics_to_grididx(large_grid, quantics) + @test recovered_grididx == grididx + end + end + + Rs_comp = (4, 4) + grid_no_endpoint = DiscretizedGrid(Rs_comp; includeendpoint=false) + grid_with_endpoint = DiscretizedGrid(Rs_comp; includeendpoint=true) + + ub_no = QuanticsGrids.upper_bound(grid_no_endpoint) + ub_with = QuanticsGrids.upper_bound(grid_with_endpoint) + @test all(ub_with .> ub_no) + + gmax_no = QuanticsGrids.grid_max(grid_no_endpoint) + gmax_with = QuanticsGrids.grid_max(grid_with_endpoint) + @test all(gmax_with .> gmax_no) + + step_no = QuanticsGrids.grid_step(grid_no_endpoint) + step_with = QuanticsGrids.grid_step(grid_with_endpoint) + @test all(step_no .> 0) && all(step_with .> 0) + + @test isa(grid, DiscretizedGrid{4}) + @test all(isfinite.(QuanticsGrids.lower_bound(grid))) + @test all(isfinite.(QuanticsGrids.upper_bound(grid))) + @test all(QuanticsGrids.grid_step(grid) .> 0) + + @test length(grid.discretegrid.lookup_table) == 4 + @test all(length(grid.discretegrid.lookup_table[d]) == QuanticsGrids.grid_Rs(grid)[d] for d in 1:4) + + base2_grid = DiscretizedGrid((5, 3); base=2, lower_bound=(0.0, -1.0), upper_bound=(3.0, 2.0)) + + for _ in 1:20 + quantics_2 = [rand(1:QuanticsGrids.sitedim(base2_grid, site)) for site in 1:length(base2_grid)] + grididx_2 = quantics_to_grididx(base2_grid, quantics_2) + back_quantics_2 = grididx_to_quantics(base2_grid, grididx_2) + @test back_quantics_2 == quantics_2 + + origcoord_2 = grididx_to_origcoord(base2_grid, grididx_2) + back_grididx_2 = QuanticsGrids.origcoord_to_grididx(base2_grid, origcoord_2) + @test back_grididx_2 == grididx_2 + end + + simple_tupletable_base2 = [ + [(:x, 1)], + [(:x, 2)], + [(:y, 1)], + [(:y, 2)] + ] + simple_grid_base2 = DiscretizedGrid((:x, :y), simple_tupletable_base2; base=2) + + test_quantics_base2 = [[1, 1, 1, 1], [2, 1, 1, 1], [1, 2, 1, 1], [2, 2, 2, 2]] + for q in test_quantics_base2 + grididx = quantics_to_grididx(simple_grid_base2, q) + back_q = grididx_to_quantics(simple_grid_base2, grididx) + @test back_q == q + end + + test_grid_2d = DiscretizedGrid((:x, :y), [[(:x, 1)], [(:y, 1)], [(:x, 2)], [(:y, 2)]]; + lower_bound=(-2.0, -3.0), upper_bound=(5.0, 7.0)) + + grididx_kw = QuanticsGrids.origcoord_to_grididx(test_grid_2d; x=1.0, y=2.0) + grididx_pos = QuanticsGrids.origcoord_to_grididx(test_grid_2d, (1.0, 2.0)) + @test grididx_kw == grididx_pos + + quantics_kw = origcoord_to_quantics(test_grid_2d; x=1.0, y=2.0) + quantics_pos = origcoord_to_quantics(test_grid_2d, (1.0, 2.0)) + @test quantics_kw == quantics_pos + + origcoord_kw = grididx_to_origcoord(test_grid_2d; x=2, y=3) + origcoord_pos = grididx_to_origcoord(test_grid_2d, (2, 3)) + @test all(origcoord_kw .≈ origcoord_pos) + + quantics_from_grididx_kw = grididx_to_quantics(test_grid_2d; x=2, y=3) + quantics_from_grididx_pos = grididx_to_quantics(test_grid_2d, (2, 3)) + @test quantics_from_grididx_kw == quantics_from_grididx_pos + + test_grid_3d = DiscretizedGrid((:a, :b, :c), [[(:a, 1)], [(:b, 1)], [(:c, 1)]]; + lower_bound=(-1.0, 0.0, 10.0), upper_bound=(1.0, 2.0, 20.0)) + + grididx_3d_kw = QuanticsGrids.origcoord_to_grididx(test_grid_3d; a=0.5, b=1.0, c=15.0) + grididx_3d_pos = QuanticsGrids.origcoord_to_grididx(test_grid_3d, (0.5, 1.0, 15.0)) + @test grididx_3d_kw == grididx_3d_pos + + @test_throws ArgumentError QuanticsGrids.origcoord_to_grididx(test_grid_2d; x=1.0) + @test_throws ArgumentError QuanticsGrids.origcoord_to_grididx(test_grid_2d; y=2.0) + @test_throws ArgumentError QuanticsGrids.origcoord_to_grididx(test_grid_2d; x=1.0, y=2.0, z=3.0) + @test_throws ArgumentError QuanticsGrids.origcoord_to_grididx(test_grid_2d; a=1.0, b=2.0) + + @test_throws ArgumentError grididx_to_origcoord(test_grid_2d; x=1) + @test_throws ArgumentError grididx_to_origcoord(test_grid_2d; y=2) + @test_throws ArgumentError grididx_to_origcoord(test_grid_2d; x=1, y=2, z=3) + @test_throws ArgumentError grididx_to_origcoord(test_grid_2d; a=1, b=2) + + for _ in 1:10 + coord_vals = (rand() * 6 - 2, rand() * 9 - 3) + grididx_tuple = QuanticsGrids.origcoord_to_grididx(test_grid_2d, coord_vals) + local grididx_kw = QuanticsGrids.origcoord_to_grididx(test_grid_2d; x=coord_vals[1], y=coord_vals[2]) + @test grididx_tuple == grididx_kw + + quantics_tuple = origcoord_to_quantics(test_grid_2d, coord_vals) + local quantics_kw = origcoord_to_quantics(test_grid_2d; x=coord_vals[1], y=coord_vals[2]) + @test quantics_tuple == quantics_kw + + origcoord_from_grid_tuple = grididx_to_origcoord(test_grid_2d, grididx_tuple) + origcoord_from_grid_kw = grididx_to_origcoord(test_grid_2d; x=grididx_tuple[1], y=grididx_tuple[2]) + @test all(origcoord_from_grid_tuple .≈ origcoord_from_grid_kw) + + quantics_from_grid_tuple = grididx_to_quantics(test_grid_2d, grididx_tuple) + quantics_from_grid_kw = grididx_to_quantics(test_grid_2d; x=grididx_tuple[1], y=grididx_tuple[2]) + @test quantics_from_grid_tuple == quantics_from_grid_kw + end + + # Test for new constructor with variablenames and Rs parameters + @testset "DiscretizedGrid constructor with variablenames and Rs" begin + variablenames = (:a, :b, :c) + Rs = (3, 4, 2) + + # Test with default parameters + grid1 = DiscretizedGrid(variablenames, Rs) + @test QuanticsGrids.grid_variablenames(grid1) == variablenames + @test QuanticsGrids.grid_Rs(grid1) == Rs + @test QuanticsGrids.grid_base(grid1) == 2 # default base + @test ndims(grid1) == 3 + + # Test with custom parameters + lower_bound = (-1.0, 0.0, -5.0) + upper_bound = (1.0, 2.0, 5.0) + base = 3 + + grid2 = DiscretizedGrid(variablenames, Rs; + lower_bound=lower_bound, + upper_bound=upper_bound, + base=base, + unfoldingscheme=:fused, + includeendpoint=true) + + @test QuanticsGrids.grid_variablenames(grid2) == variablenames + @test QuanticsGrids.grid_Rs(grid2) == Rs + @test QuanticsGrids.grid_base(grid2) == base + @test QuanticsGrids.lower_bound(grid2) == lower_bound + + # Test basic functionality works + test_grididx = (2, 3, 1) + origcoord = grididx_to_origcoord(grid2, test_grididx) + back_to_grididx = QuanticsGrids.origcoord_to_grididx(grid2, origcoord) + @test back_to_grididx == test_grididx + + # Test includeendpoint as tuple of bools + variablenames_tuple = (:x, :y) + Rs_tuple = (3, 4) + includeendpoint_tuple = (true, false) # include endpoint for x, not for y + + grid3 = DiscretizedGrid(variablenames_tuple, Rs_tuple; + includeendpoint=includeendpoint_tuple) + + @test QuanticsGrids.grid_variablenames(grid3) == variablenames_tuple + @test QuanticsGrids.grid_Rs(grid3) == Rs_tuple + + # Test that upper bounds are adjusted correctly for each dimension + lb3 = QuanticsGrids.lower_bound(grid3) + ub3 = QuanticsGrids.upper_bound(grid3) + + # For x (includeendpoint=true), upper bound should be adjusted + # For y (includeendpoint=false), upper bound should be unchanged (default 1.0) + @test lb3[1] ≈ 0.0 # default lower bound for x + @test lb3[2] ≈ 0.0 # default lower bound for y + @test ub3[2] ≈ 1.0 # default upper bound for y (unchanged due to includeendpoint=false) + + # Test functionality with tuple includeendpoint + test_grididx3 = (2, 2) + origcoord3 = grididx_to_origcoord(grid3, test_grididx3) + back_to_grididx3 = QuanticsGrids.origcoord_to_grididx(grid3, origcoord3) + @test back_to_grididx3 == test_grididx3 + end +end + +@testitem "integration test 1" begin + # This didn't work before + R = 5 + ndims = 0 + N = 2^R + fermi_frequencies_min = float(-(N - 1)) + fermi_frequencies_max = float(+(N - 1) + 2) # bc includeendpoint=false + bose_frequencies_min = float(-N) + bose_frequencies_max = float(+(N - 2) + 2) # bc includeendpoint=false + + momenta_min = ntuple(_ -> 0.0, ndims) + momenta_max = ntuple(_ -> 2.0pi, ndims) + grid_min = (fermi_frequencies_min, fermi_frequencies_min, bose_frequencies_min, momenta_min..., momenta_min..., momenta_min...) + grid_max = (fermi_frequencies_max, fermi_frequencies_max, bose_frequencies_max, momenta_max..., momenta_max..., momenta_max...) + k_syms = [Symbol("k$(dim)") for dim in 1:ndims] + k´_syms = [Symbol("k'$(dim)") for dim in 1:ndims] + q_syms = [Symbol("q$(dim)") for dim in 1:ndims] + Rs = (R, R, R, ntuple(Returns(R), ndims)..., ntuple(Returns(R), ndims)..., ntuple(Returns(R), ndims)...) + variablenames = (:v, :v´, :w, k_syms..., k´_syms..., q_syms...) + grid = DiscretizedGrid(variablenames, Rs; + lower_bound=grid_min, + upper_bound=grid_max, + unfoldingscheme=:interleaved, + includeendpoint=(true, true, true, ntuple(Returns(false), ndims)..., ntuple(Returns(false), ndims)..., ntuple(Returns(false), ndims)...)) + @test origcoord_to_quantics(grid, 0.0) == [1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] +end + +@testitem "Additional Constructors" begin + g = DiscretizedGrid{2}((3, 4), (-2, 3), (4, 5); unfoldingscheme=:interleaved) + g = DiscretizedGrid((3, 4), (-2, 3), (4, 5); unfoldingscheme=:interleaved) +end diff --git a/test/discretizedgrid_misc_tests.jl b/test/discretizedgrid_misc_tests.jl new file mode 100644 index 0000000..328bd1c --- /dev/null +++ b/test/discretizedgrid_misc_tests.jl @@ -0,0 +1,76 @@ +@testitem "DiscretizedGrid constructors error paths" begin + # upper_bound < lower_bound + @test_throws ArgumentError DiscretizedGrid(10, 0.1, 0.01) + + # upper_bound < lower_bound (multiple) dimensions + @test_throws ArgumentError DiscretizedGrid((10, 4), (0.1, 0.2), (0.01, 0.02)) + + # upper_bound == lower_bound + @test_throws ArgumentError DiscretizedGrid(10, 0.1, 0.1) + @test_throws ArgumentError DiscretizedGrid((10, 4), (0.1, 0.2), (0.1, 0.02)) + + # incompatible indextable/variablenames + @test_throws ArgumentError DiscretizedGrid((:x, :y), [[(:x, 1)], [(:z, 1)]]) + + # R = 0 and includeendpoint=true + @test_throws ArgumentError DiscretizedGrid((0, 5); includeendpoint=true) +end + +@testitem "DiscretizedGrid kwarg handling: _handle_kwargs_input" begin + g = DiscretizedGrid((7, 8, 9); variablenames=(:x, :y, :z)) + quants_kwarg = origcoord_to_quantics(g; x=0.5, z=0.75, y=0.25) + quants_plain = origcoord_to_quantics(g, (0.5, 0.25, 0.75)) + @test quants_kwarg == quants_plain + @test_throws ArgumentError origcoord_to_quantics(g; x=0.5, z=0.75) # missing y + @test_throws ArgumentError origcoord_to_quantics(g; x=0.5, z=0.75, y=0.25, extra=1) # extra keyword argument + @test_throws ArgumentError origcoord_to_quantics(g; x=:abc, y=0.3, z=0.1) + # @test_throws ArgumentError false +end + +@testitem "DiscretizedGrid sitedim" begin + g = DiscretizedGrid((7, 8, 9); variablenames=(:x, :y, :z)) + expected_sitedims = [8, 8, 8, 8, 8, 8, 8, 4, 2] + siteinds = eachindex(QuanticsGrids.grid_indextable(g)) + for i in siteinds + @test QuanticsGrids.sitedim(g, i) == expected_sitedims[i] + end + @test_throws DomainError QuanticsGrids.sitedim(g, last(siteinds) + 1) +end + +@testitem "DiscretizedGrid constructors consistency" begin + Rs = (2, 3, 4) + variablenames = (:x, :y, :z) + indextable = [[(:z, 1), (:y, 1), (:x, 1)], [(:z, 2), (:y, 2), (:x, 2)], [(:z, 3), (:y, 3)], [(:z, 4)]] + g = DiscretizedGrid(Rs; variablenames) + g´ = DiscretizedGrid(variablenames, Rs) + g´´ = DiscretizedGrid(variablenames, indextable) + @test g.discretegrid.lookup_table == g´.discretegrid.lookup_table == g´´.discretegrid.lookup_table +end + +@testitem "DiscretizedGrid 0-dimensional show method" begin + g = DiscretizedGrid(()) + @test try + sprint(show, g) + true + catch e + false + end + + g = DiscretizedGrid((0,)) + @test try + sprint(show, g) + true + catch e + false + end +end + +@testitem "DiscretizedGrid show method" begin + g = DiscretizedGrid((2, 3, 4)) + @test try + sprint(show, g) + true + catch e + false + end +end diff --git a/test/grid_tests.jl b/test/grid_tests.jl index 8911d16..671e472 100644 --- a/test/grid_tests.jl +++ b/test/grid_tests.jl @@ -1,9 +1,7 @@ - @testitem "grid.jl" begin using Test import QuanticsGrids - @testset "quanticsfunction" begin R = 8 grid = QuanticsGrids.DiscretizedGrid{1}(R, 0.0, 1.0) @@ -23,14 +21,14 @@ @testset "2D grid (large R)" for base in [2] R = 62 d = 2 - grid = QuanticsGrids.DiscretizedGrid{d}(R, 0.0, 1.0; base = base) + grid = QuanticsGrids.DiscretizedGrid{d}(R, 0.0, 1.0; base=base) @test QuanticsGrids.grididx_to_quantics(grid, ntuple(i -> base^R, d)) == fill(base^d, R) end @testset "1D grid (too large R)" begin R = 64 - @test_throws ErrorException QuanticsGrids.DiscretizedGrid{1}(R, 0.0, 1.0) + @test_throws ArgumentError QuanticsGrids.DiscretizedGrid{1}(R, 0.0, 1.0) end @testset "grid representation conversion" for R in [10] @@ -83,7 +81,7 @@ step in [(1, 1, 1), (1, 1, 2)], origin in [(1, 1, 1), (1, 1, 2)] - m = QuanticsGrids.InherentDiscreteGrid{3}(5, origin; unfoldingscheme, step = step) + m = QuanticsGrids.InherentDiscreteGrid{3}(5, origin; unfoldingscheme, step=step) @test QuanticsGrids.grid_min(m) == origin @test QuanticsGrids.grid_step(m) == step @@ -144,12 +142,12 @@ a, b; unfoldingscheme, - includeendpoint = true, + includeendpoint=true, ) @test QuanticsGrids.localdimensions(g) == fill(2, R) @test QuanticsGrids.lower_bound(g) == 0.0 - @test QuanticsGrids.upper_bound(g) == 1.0 + @test QuanticsGrids.upper_bound(g) ≈ 1.0 + dx @test QuanticsGrids.grid_min(g) == a @test QuanticsGrids.grid_max(g) == b @@ -185,16 +183,15 @@ for (c, ref) in zip(cs, refs) @inferred(QuanticsGrids.origcoord_to_grididx(g, c)) - @show QuanticsGrids.origcoord_to_grididx(g, c), ref @test all(QuanticsGrids.origcoord_to_grididx(g, c) .== ref) end - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (0.0, 0.0)) - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (0.0, 1.1)) - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (1.1, 0.0)) - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (3.0, 1.1)) - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (1.1, 3.0)) - @test_throws "Bound Error:" QuanticsGrids.origcoord_to_grididx(g, (3.0, 3.0)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (0.0, 0.0)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (0.0, 1.1)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (1.1, 0.0)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (3.0, 1.1)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (1.1, 3.0)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, (3.0, 3.0)) end @testset "2D (includeendpoint)" for unfoldingscheme in [:interleaved, :fused] @@ -208,10 +205,9 @@ a, b; unfoldingscheme, - includeendpoint = true, + includeendpoint=true, ) - @test g.includeendpoint - @test QuanticsGrids.grid_step(g) == dx + @test all(QuanticsGrids.grid_step(g) .≈ dx) if unfoldingscheme === :interleaved @test QuanticsGrids.localdimensions(g) == fill(2, d * R) @@ -255,7 +251,7 @@ @test g3 isa QuanticsGrids.InherentDiscreteGrid{1} @test QuanticsGrids.grid_origin(g3) == 1 - g4 = QuanticsGrids.InherentDiscreteGrid(4, (1, 2, 3); step = (1, 2, 1)) + g4 = QuanticsGrids.InherentDiscreteGrid(4, (1, 2, 3); step=(1, 2, 1)) @test g4 isa QuanticsGrids.InherentDiscreteGrid{3} @test QuanticsGrids.grid_origin(g4) == (1, 2, 3) @test QuanticsGrids.grid_step(g4) == (1, 2, 1) diff --git a/test/inherentdiscretegrid_tests.jl b/test/inherentdiscretegrid_tests.jl new file mode 100644 index 0000000..f7d1c54 --- /dev/null +++ b/test/inherentdiscretegrid_tests.jl @@ -0,0 +1,1023 @@ +@testitem "InherentDiscreteGrid" begin + for unfoldingscheme in [:interleaved, :fused], + step in [(1, 1, 1), (1, 1, 2)], + origin in [(1, 1, 1), (1, 1, 2)] + + m = InherentDiscreteGrid{3}(5, origin; unfoldingscheme, step=step) + @test QuanticsGrids.grid_min(m) == origin + @test QuanticsGrids.grid_step(m) == step + + if unfoldingscheme === :interleaved + @test QuanticsGrids.localdimensions(m) == fill(2, 3 * 5) + else + @test QuanticsGrids.localdimensions(m) == fill(2^3, 5) + end + + for idx in [(1, 1, 1), (1, 1, 2), (1, 25, 1), (14, 1, 1), (25, 25, 25)] + c = QuanticsGrids.grididx_to_origcoord(m, idx) + @test QuanticsGrids.origcoord_to_grididx(m, c) == idx + + q = QuanticsGrids.grididx_to_quantics(m, idx) + if unfoldingscheme === :fused + @test length(q) == 5 + else + @test length(q) == 3 * 5 + end + @test all((1 .<= q) .&& (q .<= 2^3)) + @test QuanticsGrids.quantics_to_origcoord(m, q) == c + end + end +end + +@testitem "InherentDiscreteGrid constructors" begin + # Test basic 1D constructor + R = 4 + grid_1d = InherentDiscreteGrid{1}(R, 1; base=2, step=1, unfoldingscheme=:fused) + @test grid_1d.base == 2 + @test grid_1d.origin == (1,) + @test grid_1d.step == (1,) + + # Test multi-dimensional constructor with tuple origin and step + origin_3d = (5, 10, 15) + step_3d = (2, 3, 1) + grid_3d = InherentDiscreteGrid{3}(R, origin_3d; base=3, step=step_3d, unfoldingscheme=:interleaved) + @test grid_3d.base == 3 + @test grid_3d.origin == origin_3d + @test grid_3d.step == step_3d + + # Test default parameters + grid_default = InherentDiscreteGrid{2}(R, (0, 0)) + @test grid_default.base == 2 + @test grid_default.step == (1, 1) +end + +@testitem "InherentDiscreteGrid convenience constructors" begin + # Test single integer origin constructor + R = 3 + grid_single = InherentDiscreteGrid{2}(R, 5; step=2) + @test grid_single.origin == (5, 5) + @test grid_single.step == (2, 2) + + # Test mixed single/tuple parameters + grid_mixed = InherentDiscreteGrid{3}(R, (1, 2, 3); step=4) + @test grid_mixed.origin == (1, 2, 3) + @test grid_mixed.step == (4, 4, 4) +end + +@testitem "InherentDiscreteGrid basic properties" begin + R = 5 + base = 3 + origin = (10, 20) + step = (2, 5) + + grid = InherentDiscreteGrid{2}(R, origin; base=base, step=step) + + # Test basic properties + @test ndims(grid) == 2 + @test length(grid) == R # Number of quantics sites for fused + @test grid.base == base + + # Test grid boundaries + @test QuanticsGrids.grid_min(grid) == origin + @test QuanticsGrids.grid_max(grid) == origin .+ step .* (base^R - 1) + @test QuanticsGrids.grid_step(grid) == step + @test QuanticsGrids.grid_origin(grid) == origin + + # Test with interleaved scheme + grid_interleaved = InherentDiscreteGrid{2}(R, origin; base=base, step=step, unfoldingscheme=:interleaved) + @test length(grid_interleaved) == 2 * R # Number of quantics sites for interleaved +end + +@testitem "InherentDiscreteGrid coordinate conversions" begin + R = 4 + base = 2 + origin = (5, 10, 15) + step = (2, 3, 1) + + grid = InherentDiscreteGrid{3}(R, origin; base=base, step=step) + + # Test grididx_to_origcoord and origcoord_to_grididx roundtrip + test_grididx = [(1, 1, 1), (1, 1, 2), (1, 8, 1), (4, 1, 1), (16, 16, 16)] + for grididx in test_grididx + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end + + # Test origcoord calculation formula + grididx = (3, 5, 7) + expected_origcoord = origin .+ step .* (grididx .- 1) + @test QuanticsGrids.grididx_to_origcoord(grid, grididx) == expected_origcoord + + # Test boundary cases + min_grididx = (1, 1, 1) + max_grididx = (base^R, base^R, base^R) + @test QuanticsGrids.grididx_to_origcoord(grid, min_grididx) == origin + @test QuanticsGrids.grididx_to_origcoord(grid, max_grididx) == origin .+ step .* (base^R - 1) +end + +@testitem "InherentDiscreteGrid quantics conversions" begin + # Test fused quantics + R = 3 + base = 2 + grid_fused = InherentDiscreteGrid{2}(R, (1, 1); base=base, unfoldingscheme=:fused) + + # Test grididx_to_quantics and quantics_to_grididx roundtrip + test_grididx = [(1, 1), (1, 2), (2, 1), (4, 8), (8, 8)] + for grididx in test_grididx + quantics = QuanticsGrids.grididx_to_quantics(grid_fused, grididx) + @test QuanticsGrids.quantics_to_grididx(grid_fused, quantics) == grididx + @test length(quantics) == R # Fused should have R quantics + @test all(1 .<= quantics .<= base^2) # Each quantics should be in valid range + end + + # Test interleaved quantics + grid_interleaved = InherentDiscreteGrid{2}(R, (1, 1); base=base, unfoldingscheme=:interleaved) + + for grididx in test_grididx + quantics = QuanticsGrids.grididx_to_quantics(grid_interleaved, grididx) + @test QuanticsGrids.quantics_to_grididx(grid_interleaved, quantics) == grididx + @test length(quantics) == 2 * R # Interleaved should have 2*R quantics + @test all(1 .<= quantics .<= base) # Each quantics should be in valid range + end + + quantics_wronglength = [1, 2, 3] + @test_throws ArgumentError QuanticsGrids.quantics_to_grididx(grid_interleaved, quantics_wronglength) +end + +@testitem "InherentDiscreteGrid quantics_to_origcoord and origcoord_to_quantics" begin + R = 4 + origin = (3, 7) + step = (2, 3) + grid = InherentDiscreteGrid{2}(R, origin; step=step) + + # Test quantics_to_origcoord and origcoord_to_quantics roundtrip + for _ in 1:20 + # Generate random valid quantics + quantics = rand(1:2^2, R) # Each quantics element should be valid for the base and dimension + + origcoord = QuanticsGrids.quantics_to_origcoord(grid, quantics) + @test QuanticsGrids.origcoord_to_quantics(grid, origcoord) == quantics + + # Verify the conversion chain: quantics -> grididx -> origcoord + grididx = QuanticsGrids.quantics_to_grididx(grid, quantics) + expected_origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test origcoord == expected_origcoord + end +end + +@testitem "InherentDiscreteGrid localdimensions" begin + R = 4 + base = 3 + + # Test fused scheme + grid_fused_1d = InherentDiscreteGrid{1}(R, 1; base=base, unfoldingscheme=:fused) + @test QuanticsGrids.localdimensions(grid_fused_1d) == fill(base, R) + + grid_fused_3d = InherentDiscreteGrid{3}(R, (1, 1, 1); base=base, unfoldingscheme=:fused) + @test QuanticsGrids.localdimensions(grid_fused_3d) == fill(base^3, R) + + # Test interleaved scheme + grid_interleaved_1d = InherentDiscreteGrid{1}(R, 1; base=base, unfoldingscheme=:interleaved) + @test QuanticsGrids.localdimensions(grid_interleaved_1d) == fill(base, R) + + grid_interleaved_3d = InherentDiscreteGrid{3}(R, (1, 1, 1); base=base, unfoldingscheme=:interleaved) + @test QuanticsGrids.localdimensions(grid_interleaved_3d) == fill(base, 3 * R) + + # Test different bases + for test_base in [2, 4, 5] + grid = InherentDiscreteGrid{2}(R, (1, 1); base=test_base, unfoldingscheme=:fused) + @test QuanticsGrids.localdimensions(grid) == fill(test_base^2, R) + end +end + +@testitem "InherentDiscreteGrid different bases" begin + R = 3 + origin = (1, 1) + + for base in [2, 3, 4, 5] + grid = InherentDiscreteGrid{2}(R, origin; base=base) + + # Test that grid can handle different bases + max_grididx = (base^R, base^R) + min_grididx = (1, 1) + + # Test coordinate conversion works with different bases + max_origcoord = QuanticsGrids.grididx_to_origcoord(grid, max_grididx) + min_origcoord = QuanticsGrids.grididx_to_origcoord(grid, min_grididx) + + @test QuanticsGrids.origcoord_to_grididx(grid, max_origcoord) == max_grididx + @test QuanticsGrids.origcoord_to_grididx(grid, min_origcoord) == min_grididx + + # Test quantics conversion with different bases + quantics = QuanticsGrids.grididx_to_quantics(grid, (2, 3)) + @test all(1 .<= quantics .<= base^2) # For fused scheme + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == (2, 3) + end +end + +@testitem "InherentDiscreteGrid scalar vs tuple return values" begin + R = 3 + + # Test 1D grid returns scalars + grid_1d = InherentDiscreteGrid{1}(R, 5; step=2) + grididx_1d = 3 + origcoord_1d = QuanticsGrids.grididx_to_origcoord(grid_1d, grididx_1d) + @test origcoord_1d isa Int + @test QuanticsGrids.origcoord_to_grididx(grid_1d, origcoord_1d) == grididx_1d + + @test QuanticsGrids.grid_min(grid_1d) isa Int + @test QuanticsGrids.grid_max(grid_1d) isa Int + @test QuanticsGrids.grid_step(grid_1d) isa Int + @test QuanticsGrids.grid_origin(grid_1d) isa Int + + # Test multi-dimensional grid returns tuples + grid_2d = InherentDiscreteGrid{2}(R, (3, 7); step=(1, 2)) + grididx_2d = (2, 4) + origcoord_2d = QuanticsGrids.grididx_to_origcoord(grid_2d, grididx_2d) + @test origcoord_2d isa Tuple{Int,Int} + @test QuanticsGrids.origcoord_to_grididx(grid_2d, origcoord_2d) == grididx_2d + + @test QuanticsGrids.grid_min(grid_2d) isa Tuple{Int,Int} + @test QuanticsGrids.grid_max(grid_2d) isa Tuple{Int,Int} + @test QuanticsGrids.grid_step(grid_2d) isa Tuple{Int,Int} + @test QuanticsGrids.grid_origin(grid_2d) isa Tuple{Int,Int} +end + +@testitem "InherentDiscreteGrid comprehensive conversion test" begin + # Test all possible conversion combinations like in the old grid tests + R = 4 + base = 3 + origin = (2, 5) + step = (3, 2) + + grid = InherentDiscreteGrid{2}(R, origin; base=base, step=step) + + reprs = [:grididx, :quantics, :origcoord] + + transforms = Dict( + (:grididx, :quantics) => QuanticsGrids.grididx_to_quantics, + (:grididx, :origcoord) => QuanticsGrids.grididx_to_origcoord, + (:quantics, :grididx) => QuanticsGrids.quantics_to_grididx, + (:quantics, :origcoord) => QuanticsGrids.quantics_to_origcoord, + (:origcoord, :grididx) => QuanticsGrids.origcoord_to_grididx, + (:origcoord, :quantics) => QuanticsGrids.origcoord_to_quantics, + ) + + initial_grididx = (2, 3) + data = Dict{Symbol,Any}() + data[:grididx] = initial_grididx + + # Build up all representations + while true + newdata = false + for src in reprs, dst in reprs + if src == dst + continue + end + if !(src ∈ keys(data)) + continue + end + + if dst ∈ keys(data) + @test transforms[(src, dst)](grid, data[src]) == data[dst] + else + data[dst] = transforms[(src, dst)](grid, data[src]) + newdata = true + end + end + if newdata == false + break + end + end + + # Verify all conversions are consistent + @test length(data) == 3 + @test haskey(data, :grididx) + @test haskey(data, :quantics) + @test haskey(data, :origcoord) +end + +@testitem "InherentDiscreteGrid edge cases and boundary conditions" begin + R = 3 + base = 2 + + # Test minimum grid size (R=1) + grid_min = InherentDiscreteGrid{1}(1, 10; base=base, step=5) + @test QuanticsGrids.grididx_to_origcoord(grid_min, 1) == 10 + @test QuanticsGrids.grididx_to_origcoord(grid_min, base) == 10 + 5 * (base - 1) + + # Test large step sizes + grid_large_step = InherentDiscreteGrid{2}(R, (0, 0); step=(100, 200)) + grididx = (2, 3) + expected_coord = (100, 400) # (0 + 100*(2-1), 0 + 200*(3-1)) + @test QuanticsGrids.grididx_to_origcoord(grid_large_step, grididx) == expected_coord + + # Test negative origins + grid_neg = InherentDiscreteGrid{2}(R, (-10, -20); step=(2, 3)) + min_coord = QuanticsGrids.grididx_to_origcoord(grid_neg, (1, 1)) + @test min_coord == (-10, -20) + + # Test single dimension + grid_1d = InherentDiscreteGrid{1}(R, 5; step=3) + for i in 1:base^R + coord = QuanticsGrids.grididx_to_origcoord(grid_1d, i) + @test QuanticsGrids.origcoord_to_grididx(grid_1d, coord) == i + end +end + +@testitem "InherentDiscreteGrid consistency with unfolding schemes" begin + R = 4 + origin = (1, 1, 1) + step = (1, 2, 3) + base = 2 + + grid_fused = InherentDiscreteGrid{3}(R, origin; step=step, base=base, unfoldingscheme=:fused) + grid_interleaved = InherentDiscreteGrid{3}(R, origin; step=step, base=base, unfoldingscheme=:interleaved) + + # Test that coordinate conversions are consistent between schemes + test_grididx = [(1, 1, 1), (2, 3, 4), (8, 16, 16)] + for grididx in test_grididx + coord_fused = QuanticsGrids.grididx_to_origcoord(grid_fused, grididx) + coord_interleaved = QuanticsGrids.grididx_to_origcoord(grid_interleaved, grididx) + @test coord_fused == coord_interleaved + + @test QuanticsGrids.origcoord_to_grididx(grid_fused, coord_fused) == grididx + @test QuanticsGrids.origcoord_to_grididx(grid_interleaved, coord_interleaved) == grididx + end + + # Test that quantics differ but origcoord results are the same + grididx = (3, 5, 7) + quantics_fused = QuanticsGrids.grididx_to_quantics(grid_fused, grididx) + quantics_interleaved = QuanticsGrids.grididx_to_quantics(grid_interleaved, grididx) + + @test length(quantics_fused) == R + @test length(quantics_interleaved) == 3 * R + @test quantics_fused != quantics_interleaved # Should be different + + # But they should convert to same coordinates + coord_from_fused = QuanticsGrids.quantics_to_origcoord(grid_fused, quantics_fused) + coord_from_interleaved = QuanticsGrids.quantics_to_origcoord(grid_interleaved, quantics_interleaved) + @test coord_from_fused == coord_from_interleaved +end + +@testitem "InherentDiscreteGrid high-dimensional grids" begin + R = 2 # Keep R small for high dimensions to avoid memory issues + base = 2 + D = 5 # 5-dimensional grid + + origin = ntuple(i -> i, D) # (1, 2, 3, 4, 5) + step = ntuple(i -> i, D) # (1, 2, 3, 4, 5) + + grid = InherentDiscreteGrid{D}(R, origin; step=step, base=base) + + @test ndims(grid) == D + @test QuanticsGrids.grid_min(grid) == origin + @test QuanticsGrids.grid_step(grid) == step + + # Test coordinate conversion + grididx = ntuple(i -> 2, D) # All 2s + expected_origcoord = origin .+ step .* (grididx .- 1) + @test QuanticsGrids.grididx_to_origcoord(grid, grididx) == expected_origcoord + @test QuanticsGrids.origcoord_to_grididx(grid, expected_origcoord) == grididx + + # Test quantics + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + @test length(quantics) == R # Fused scheme +end + +@testitem "InherentDiscreteGrid performance and stress tests" begin + R = 5 + base = 2 + grid = InherentDiscreteGrid{3}(R, (1, 1, 1); base=base) + + # Test random conversions for consistency + for _ in 1:100 + # Generate random valid grid indices + grididx = ntuple(3) do d + rand(1:base^R) + end + + # Test roundtrip conversions + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + + # Test cross-conversions + origcoord_from_quantics = QuanticsGrids.quantics_to_origcoord(grid, quantics) + @test origcoord_from_quantics == origcoord + + quantics_from_origcoord = QuanticsGrids.origcoord_to_quantics(grid, origcoord) + @test quantics_from_origcoord == quantics + end +end + +@testitem "InherentDiscreteGrid grid with zero step (degenerate case)" begin + # This might be an edge case that should either work or throw a meaningful error + R = 3 + origin = (5, 10) + + @test_throws ArgumentError InherentDiscreteGrid{2}(R, origin; step=(0, 1)) +end + +@testitem "InherentDiscreteGrid with custom indextable - basic" begin + # Test with a simple custom indextable + R = 3 + variablenames = (:x, :y) + indextable = [ + [(:x, 1)], # site 1: x_1 + [(:x, 2)], # site 2: x_2 + [(:y, 1)], # site 3: y_1 + [(:x, 3), (:y, 2)] # site 4: x_3, y_2 + ] + origin = (1, 5) + step = (2, 3) + base = 2 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 2) # x has 3 quantics, y has 2 quantics + @test grid.origin == origin + @test grid.step == step + @test grid.base == base + @test length(grid) == 4 # 4 sites in tensor train + + # Test local dimensions: [2^1, 2^1, 2^1, 2^2] = [2, 2, 2, 4] + @test QuanticsGrids.localdimensions(grid) == [2, 2, 2, 4] + + # Test coordinate conversions work + grididx = (3, 2) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + # Test quantics conversion + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test length(quantics) == 4 + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + @test QuanticsGrids.quantics_to_origcoord(grid, quantics) == origcoord +end + +@testitem "InherentDiscreteGrid with complex indextable - mixed sites" begin + # Test with a more complex indextable with mixed site sizes + variablenames = (:a, :b, :c) + indextable = [ + [(:a, 1), (:b, 1)], # site 1: a_1, b_1 (2 indices) + [(:c, 1)], # site 2: c_1 (1 index) + [(:a, 2), (:b, 2), (:c, 2)], # site 3: a_2, b_2, c_2 (3 indices) + [(:a, 3)], # site 4: a_3 (1 index) + [(:b, 3), (:c, 3)] # site 5: b_3, c_3 (2 indices) + ] + origin = (10, 20, 30) + step = (1, 2, 5) + base = 3 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 3, 3) # Each variable has 3 quantics + @test grid.origin == origin + @test grid.step == step + @test grid.base == base + @test length(grid) == 5 # 5 sites in tensor train + + # Test local dimensions: [3^2, 3^1, 3^3, 3^1, 3^2] = [9, 3, 27, 3, 9] + expected_localdims = [9, 3, 27, 3, 9] + @test QuanticsGrids.localdimensions(grid) == expected_localdims + + # Test that all quantics sites have correct dimensions + for (i, expected_dim) in enumerate(expected_localdims) + @test QuanticsGrids.sitedim(grid, i) == expected_dim + end + + @test_throws DomainError QuanticsGrids.sitedim(grid, length(expected_localdims) + 1) + + # Test coordinate conversions for various grid indices + test_grididx = [(1, 1, 1), (2, 3, 1), (9, 27, 27), (5, 10, 15)] + for grididx in test_grididx + if all(grididx .<= base .^ grid.Rs) # Only test valid indices + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test length(quantics) == 5 + @test all(1 .<= quantics[i] .<= expected_localdims[i] for i in 1:5) + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + end + end +end + +@testitem "InherentDiscreteGrid with asymmetric indextable" begin + # Test with highly asymmetric resolution per variable + variablenames = (:x, :y, :z) + indextable = [ + [(:x, 1), (:x, 2), (:x, 3), (:x, 4)], # site 1: all x quantics (4 indices) + [(:y, 1)], # site 2: y_1 (1 index) + [(:y, 2)], # site 3: y_2 (1 index) + [(:z, 1), (:z, 2)] # site 4: z_1, z_2 (2 indices) + ] + origin = (0, 100, -50) + step = (10, 5, 2) + base = 2 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (4, 2, 2) # x has 4 quantics, y and z have 2 each + @test length(grid) == 4 # 4 sites + + # Test local dimensions: [2^4, 2^1, 2^1, 2^2] = [16, 2, 2, 4] + @test QuanticsGrids.localdimensions(grid) == [16, 2, 2, 4] + + # Test that different variables have different resolutions + max_grididx_x = base^4 # 16 + max_grididx_y = base^2 # 4 + max_grididx_z = base^2 # 4 + + # Test boundary cases + grididx = (max_grididx_x, max_grididx_y, max_grididx_z) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + expected_origcoord = origin .+ step .* (grididx .- 1) + @test origcoord == expected_origcoord + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx +end + +@testitem "InherentDiscreteGrid with single-site indextable" begin + # Test edge case: all quantics in a single site + variablenames = (:x, :y) + indextable = [ + [(:x, 1), (:x, 2), (:x, 3), (:y, 1), (:y, 2)] # All quantics in one site + ] + origin = (5, 10) + step = (1, 3) + base = 2 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 2) # x has 3 quantics, y has 2 quantics + @test length(grid) == 1 # Only 1 site + + # Test local dimensions: [2^5] = [32] (5 total quantics in one site) + @test QuanticsGrids.localdimensions(grid) == [32] + + # Test that conversions still work correctly + grididx = (4, 3) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test length(quantics) == 1 + @test 1 <= quantics[1] <= 32 + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx +end + +@testitem "InherentDiscreteGrid with maximum fragmentation" begin + # Test extreme case: each quantics in its own site + variablenames = (:a, :b) + indextable = [ + [(:a, 1)], # site 1: a_1 + [(:a, 2)], # site 2: a_2 + [(:a, 3)], # site 3: a_3 + [(:b, 1)], # site 4: b_1 + [(:b, 2)] # site 5: b_2 + ] + origin = (1, 1) + step = (1, 1) + base = 3 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 2) # a has 3 quantics, b has 2 quantics + @test length(grid) == 5 # 5 sites + + # Test local dimensions: all sites have single quantics, so [3, 3, 3, 3, 3] + @test QuanticsGrids.localdimensions(grid) == fill(3, 5) + + # Test conversions + for grididx in [(1, 1), (3, 3), (9, 9), (27, 9)] + if all(grididx .<= base .^ grid.Rs) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test length(quantics) == 5 + @test all(1 .<= quantics .<= 3) + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + end + end +end + +@testitem "InherentDiscreteGrid complex indextable with base != 2" begin + # Test complex indextable with non-binary base + variablenames = (:u, :v, :w) + indextable = [ + [(:u, 1), (:v, 1)], # site 1: u_1, v_1 + [(:w, 1), (:w, 2)], # site 2: w_1, w_2 + [(:u, 2)], # site 3: u_2 + [(:v, 2), (:u, 3), (:w, 3)] # site 4: v_2, u_3, w_3 + ] + origin = (-5, 0, 10) + step = (3, 7, 2) + base = 5 # Base 5 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 2, 3) # u has 3, v has 2, w has 3 quantics + @test grid.base == 5 + @test length(grid) == 4 + + # Test local dimensions: [5^2, 5^2, 5^1, 5^3] = [25, 25, 5, 125] + @test QuanticsGrids.localdimensions(grid) == [25, 25, 5, 125] + + # Test with base-5 specific values + grididx = (1, 1, 1) # Minimum + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test origcoord == origin + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + # Test maximum valid indices + max_grididx = (base^3, base^2, base^3) # (125, 25, 125) + max_origcoord = QuanticsGrids.grididx_to_origcoord(grid, max_grididx) + expected_max = origin .+ step .* (max_grididx .- 1) + @test max_origcoord == expected_max + @test QuanticsGrids.origcoord_to_grididx(grid, max_origcoord) == max_grididx + + # Test quantics conversions + test_grididx = (10, 5, 20) + quantics = QuanticsGrids.grididx_to_quantics(grid, test_grididx) + @test length(quantics) == 4 + @test all(1 .<= quantics[1] .<= 25) # site 1 + @test all(1 .<= quantics[2] .<= 25) # site 2 + @test all(1 .<= quantics[3] .<= 5) # site 3 + @test all(1 .<= quantics[4] .<= 125) # site 4 + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == test_grididx +end + +@testitem "InherentDiscreteGrid comprehensive indextable test" begin + # Test a realistic complex case similar to what might be used in practice + variablenames = (:k_x, :k_y, :ω, :t) + indextable = [ + [(:k_x, 1), (:k_y, 1)], # site 1: spatial momenta first bits + [(:ω, 1), (:t, 1)], # site 2: frequency and time first bits + [(:k_x, 2)], # site 3: k_x second bit + [(:k_y, 2), (:ω, 2)], # site 4: k_y and ω second bits + [(:t, 2), (:t, 3)], # site 5: time higher resolution + [(:k_x, 3), (:k_y, 3), (:ω, 3)] # site 6: spatial and frequency high bits + ] + origin = (0, 0, -10, 0) + step = (1, 1, 2, 5) + base = 2 + + grid = InherentDiscreteGrid(variablenames, indextable; + origin=origin, step=step, base=base) + + @test grid.variablenames == variablenames + @test grid.Rs == (3, 3, 3, 3) # All variables have 3 quantics + @test length(grid) == 6 # 6 sites + + # Test local dimensions: [2^2, 2^2, 2^1, 2^2, 2^2, 2^3] = [4, 4, 2, 4, 4, 8] + expected_localdims = [4, 4, 2, 4, 4, 8] + @test QuanticsGrids.localdimensions(grid) == expected_localdims + + # Test random conversions for consistency + for _ in 1:50 + # Generate random valid grid indices + grididx = ntuple(4) do d + rand(1:base^3) # Each dimension has 3 quantics, so max index is 2^3 = 8 + end + + # Test all conversion roundtrips + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + quantics = QuanticsGrids.grididx_to_quantics(grid, grididx) + @test length(quantics) == 6 + @test all(1 .<= quantics[i] .<= expected_localdims[i] for i in 1:6) + @test QuanticsGrids.quantics_to_grididx(grid, quantics) == grididx + + # Cross conversions + @test QuanticsGrids.quantics_to_origcoord(grid, quantics) == origcoord + @test QuanticsGrids.origcoord_to_quantics(grid, origcoord) == quantics + end +end + +@testitem "InherentDiscreteGrid indextable validation" begin + # Test that constructor properly validates indextable consistency + variablenames = (:x, :y) + + # Valid indextable + valid_indextable = [ + [(:x, 1), (:y, 1)], + [(:x, 2)], + [(:y, 2)] + ] + + grid = InherentDiscreteGrid(variablenames, valid_indextable; + origin=(1, 1), step=(1, 1)) + @test grid.Rs == (2, 2) # x and y each have 2 quantics + + # Test that unknown variables are detected + invalid_indextable_unknown = [ + [(:x, 1), (:y, 1)], + [(:z, 1)] # z is not in variablenames + ] + + @test_throws ArgumentError InherentDiscreteGrid(variablenames, invalid_indextable_unknown; + origin=(1, 1), step=(1, 1)) + + invalid_indextable_missing_index = [ + [(:x, 1), (:y, 1)], + [(:x, 3), (:y, 2)] + ] + + @test_throws ArgumentError InherentDiscreteGrid( + variablenames, invalid_indextable_missing_index) + + invalid_indextable_repeated_index = [ + [(:x, 1), (:y, 1)], + [(:x, 1), (:y, 2)] + ] + + @test_throws ArgumentError InherentDiscreteGrid( + variablenames, invalid_indextable_repeated_index) +end + +@testitem "InherentDiscreteGrid constructor with variablenames and Rs" begin + # Test basic constructor with explicit variable names and resolutions + variablenames = (:k_x, :k_y, :ω) + Rs = (4, 3, 5) + origin = (1, 10, -5) + step = (2, 1, 3) + + # Test with default parameters + grid1 = InherentDiscreteGrid(variablenames, Rs; origin=origin, step=step) + @test grid1.variablenames == variablenames + @test grid1.Rs == Rs + @test grid1.base == 2 # default base + @test grid1.origin == origin + @test grid1.step == step + @test ndims(grid1) == 3 + + # Test with custom parameters + base = 3 + unfoldingscheme = :interleaved + + grid2 = InherentDiscreteGrid(variablenames, Rs; + origin=origin, + step=step, + base=base, + unfoldingscheme=unfoldingscheme) + + @test grid2.variablenames == variablenames + @test grid2.Rs == Rs + @test grid2.base == base + @test grid2.origin == origin + @test grid2.step == step + + # Test that functionality works + test_grididx = (2, 3, 4) + origcoord = QuanticsGrids.grididx_to_origcoord(grid2, test_grididx) + back_to_grididx = QuanticsGrids.origcoord_to_grididx(grid2, origcoord) + @test back_to_grididx == test_grididx + + # Test quantics conversion + quantics = QuanticsGrids.grididx_to_quantics(grid2, test_grididx) + @test QuanticsGrids.quantics_to_grididx(grid2, quantics) == test_grididx +end + +@testitem "InherentDiscreteGrid constructor with Rs tuple only" begin + # Test constructor that auto-generates variable names + Rs = (3, 4, 2) + origin = (5, 0, -10) + step = (1, 2, 5) + base = 3 + + grid = InherentDiscreteGrid(Rs; origin=origin, step=step, base=base) + + @test grid.Rs == Rs + @test grid.origin == origin + @test grid.step == step + @test grid.base == base + @test ndims(grid) == 3 + + # Should auto-generate variable names as symbols + expected_names = (Symbol(1), Symbol(2), Symbol(3)) + @test grid.variablenames == expected_names + + # Test that grid works correctly + grididx = (2, 3, 1) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx +end + +@testitem "InherentDiscreteGrid 1D constructor with R::Int" begin + # Test 1D constructor similar to DiscretizedGrid(R, lower, upper) + R = 5 + origin = 10 + step = 3 + base = 2 + + grid = InherentDiscreteGrid(R, origin; step=step, base=base) + + @test grid.Rs == (R,) + @test grid.origin == (origin,) + @test grid.step == (step,) + @test grid.base == base + @test ndims(grid) == 1 + + # Should auto-generate 1D variable name + @test grid.variablenames == (Symbol(1),) + + # Test conversions work for 1D case (should return scalars) + grididx = 3 + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test origcoord isa Int + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + # Test expected coordinate calculation + expected_coord = origin + step * (grididx - 1) + @test origcoord == expected_coord +end + +@testitem "InherentDiscreteGrid parametric constructor InherentDiscreteGrid{D}" begin + # Test parametric type constructor like DiscretizedGrid{D} + D = 4 + R = 3 + origin = ntuple(i -> i * 5, D) # (5, 10, 15, 20) + step = ntuple(i -> i, D) # (1, 2, 3, 4) + base = 3 + + grid = InherentDiscreteGrid{D}(R, origin; step=step, base=base, unfoldingscheme=:interleaved) + + @test ndims(grid) == D + @test grid.Rs == ntuple(_ -> R, D) # All dimensions have same R + @test grid.origin == origin + @test grid.step == step + @test grid.base == base + + # Should auto-generate variable names + expected_names = ntuple(Symbol, D) + @test grid.variablenames == expected_names + + # Test that grid works for high dimensions + grididx = ntuple(i -> 2, D) + origcoord = QuanticsGrids.grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx +end + +@testitem "InherentDiscreteGrid default parameter behavior" begin + # Test various combinations of default parameters + Rs = (3, 2) + + # Test with minimal parameters (all defaults) + grid_minimal = InherentDiscreteGrid(Rs) + @test grid_minimal.Rs == Rs + @test grid_minimal.base == 2 # default + @test grid_minimal.origin == (1, 1) # default origin should be (1, 1, ...) + @test grid_minimal.step == (1, 1) # default step should be (1, 1, ...) + + # Test with partial parameters + custom_origin = (5, -2) + grid_partial = InherentDiscreteGrid(Rs; origin=custom_origin) + @test grid_partial.origin == custom_origin + @test grid_partial.step == (1, 1) # should still be default + @test grid_partial.base == 2 # should still be default + + # Test step can be specified as single value or tuple + grid_step_single = InherentDiscreteGrid(Rs; step=3) + @test grid_step_single.step == (3, 3) # should broadcast to all dimensions + + grid_step_tuple = InherentDiscreteGrid(Rs; step=(2, 5)) + @test grid_step_tuple.step == (2, 5) + + # Test origin can be specified as single value or tuple + grid_origin_single = InherentDiscreteGrid(Rs; origin=7) + @test grid_origin_single.origin == (7, 7) # should broadcast + + grid_origin_tuple = InherentDiscreteGrid(Rs; origin=(-1, 3)) + @test grid_origin_tuple.origin == (-1, 3) +end + +@testitem "InherentDiscreteGrid constructor error handling" begin + # Test that constructors properly validate input parameters + + # Test invalid base + @test_throws ArgumentError InherentDiscreteGrid((3, 2); base=1) # base must be > 1 + @test_throws ArgumentError InherentDiscreteGrid((3, 2); base=0) + + # Test mismatched dimensions + @test_throws MethodError InherentDiscreteGrid((3, 2); origin=(1, 1, 1)) # wrong origin length + @test_throws MethodError InherentDiscreteGrid((3, 2); step=(1, 1, 1)) # wrong step length + + # Test negative R values + @test_throws ArgumentError InherentDiscreteGrid((-1, 2)) # negative R + @test_throws ArgumentError InherentDiscreteGrid((3, -2)) + + # Test too large R value + R = 62 + @test ( + try + InherentDiscreteGrid((R,)) + true + catch + false + end + ) + R = 63 + @test_throws ArgumentError InherentDiscreteGrid((R,)) + + # Test repeated variable names + @test_throws ArgumentError InherentDiscreteGrid((3, 3); variablenames=(:x, :x)) + + # Test invalid unfoldingscheme + @test_throws ArgumentError InherentDiscreteGrid((3, 2); unfoldingscheme=:invalid) +end + +@testitem "InherentDiscreteGrid constructor compatibility with existing patterns" begin + # Test that constructors work with patterns used in existing codebase + + # Pattern 1: Like InherentDiscreteGrid{d}(R, origin; kwargs...) + R = 4 + origin_2d = (1, 5) + step_2d = (2, 3) + + grid_pattern1 = InherentDiscreteGrid{2}(R, origin_2d; step=step_2d, base=3) + @test ndims(grid_pattern1) == 2 + @test grid_pattern1.Rs == (R, R) + @test grid_pattern1.origin == origin_2d + @test grid_pattern1.step == step_2d + + # Pattern 2: Like DiscretizedGrid((R1, R2, ...); kwargs...) + Rs_pattern2 = (3, 5, 2) + grid_pattern2 = InherentDiscreteGrid(Rs_pattern2; origin=(0, 0, 0), unfoldingscheme=:interleaved) + @test grid_pattern2.Rs == Rs_pattern2 + + # Pattern 3: Single R for 1D + R_1d = 6 + origin_1d = 10 + grid_pattern3 = InherentDiscreteGrid(R_1d, origin_1d; step=2, base=4) + @test ndims(grid_pattern3) == 1 + @test grid_pattern3.Rs == (R_1d,) + @test grid_pattern3.origin == (origin_1d,) + @test grid_pattern3.base == 4 + + # Test that all patterns produce working grids + for grid in [grid_pattern1, grid_pattern2, grid_pattern3] + # Test basic coordinate conversion works + if ndims(grid) == 1 + test_idx = 2 + elseif ndims(grid) == 2 + test_idx = (2, 3) + else + test_idx = (2, 3, 1) + end + + origcoord = QuanticsGrids.grididx_to_origcoord(grid, test_idx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == test_idx + end +end + +@testitem "InherentDiscreteGrid constructor with mixed parameter types" begin + # Test constructors handle mixed integer/tuple parameters correctly + + # Single R with multi-dimensional origin/step + R = 3 + origin_multi = (5, 10, 15) + step_multi = (1, 2, 3) + + grid_mixed = InherentDiscreteGrid{3}(R, origin_multi; step=step_multi) + @test grid_mixed.Rs == (R, R, R) # R should be replicated + @test grid_mixed.origin == origin_multi + @test grid_mixed.step == step_multi + + # Multi-R with single origin/step values + Rs_multi = (2, 4, 3) + origin_single = 7 + step_single = 2 + + grid_broadcast = InherentDiscreteGrid(Rs_multi; origin=origin_single, step=step_single) + @test grid_broadcast.Rs == Rs_multi + @test grid_broadcast.origin == (origin_single, origin_single, origin_single) + @test grid_broadcast.step == (step_single, step_single, step_single) + + # Test functionality works correctly + grididx = (1, 2, 3) + expected_origcoord = (origin_single, origin_single + step_single * 1, origin_single + step_single * 2) + actual_origcoord = QuanticsGrids.grididx_to_origcoord(grid_broadcast, grididx) + @test actual_origcoord == expected_origcoord +end + +@testitem "InherentDiscreteGrid 0-dimensional" begin + g = InherentDiscreteGrid(()) +end diff --git a/test/origcoord_floating_point_tests.jl b/test/origcoord_floating_point_tests.jl new file mode 100644 index 0000000..d3c03b2 --- /dev/null +++ b/test/origcoord_floating_point_tests.jl @@ -0,0 +1,372 @@ +@testitem "floating point precision - coordinates exactly at boundaries" begin + using QuanticsGrids + # Test behavior when coordinates are exactly at grid boundaries + R = 10 + grid = DiscretizedGrid{1}(R) + + # Test exact lower bound + @test QuanticsGrids.origcoord_to_grididx(grid, 0.0) == 1 + + # Test exact upper bound (should be at grid_max, not upper_bound) + grid_max_val = QuanticsGrids.grid_max(grid) + @test QuanticsGrids.origcoord_to_grididx(grid, grid_max_val) == 2^R + + # Test coordinates that are exactly representable as grid points + step = QuanticsGrids.grid_step(grid) + for i in 1:min(100, 2^R) + exact_coord = (i - 1) * step + @test QuanticsGrids.origcoord_to_grididx(grid, exact_coord) == i + end +end + +@testitem "floating point precision - very small step sizes" begin + using QuanticsGrids + # Test with large R values that create very small step sizes + R = 50 # Creates step size of 1/2^50 ≈ 8.9e-16 + grid = DiscretizedGrid{1}(R) + + step = QuanticsGrids.grid_step(grid) + @test step ≈ 1.0 / 2^R + + # Test coordinates that should map to specific grid indices + # Use Float64 arithmetic carefully to avoid precision loss + test_indices = [1, 2, 3, 2^10, 2^20, 2^R ÷ 2, 2^R - 1, 2^R] + + for idx in test_indices + # Calculate coordinate using exact arithmetic when possible + coord = (idx - 1) * step + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + @test recovered_idx == idx + end +end + +@testitem "floating point precision - catastrophic cancellation" begin + using QuanticsGrids + # Test cases where subtraction coordinate - lower_bound might lose precision + + # Case 1: Very close lower_bound and coordinate values + lower = 1e15 + upper = 1e15 + 1.0 + grid = DiscretizedGrid{1}(10; lower_bound=(lower,), upper_bound=(upper,)) + + # Test coordinates very close to lower bound + eps_values = [1e-15, 1e-14, 1e-13, 1e-12, 1e-10] + for eps in eps_values + coord = lower + eps + if coord <= QuanticsGrids.upper_bound(grid) + idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + @test idx >= 1 && idx <= 2^10 + + # Test round-trip + recovered_coord = grididx_to_origcoord(grid, idx) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, recovered_coord) + @test recovered_idx == idx + end + end + + # Case 2: Very large coordinate values where precision is lost + lower = 1e12 + upper = 1e12 + 1e6 + grid2 = DiscretizedGrid{1}(15; lower_bound=(lower,), upper_bound=(upper,)) + + step = QuanticsGrids.grid_step(grid2) + for i in [1, 100, 1000, 2^15 ÷ 2, 2^15] + coord = lower + (i - 1) * step + if coord <= QuanticsGrids.upper_bound(grid2) + idx = QuanticsGrids.origcoord_to_grididx(grid2, coord) + @test abs(idx - i) <= 1 + end + end +end + +@testitem "floating point precision - round-trip consistency stress test" begin + using QuanticsGrids + # Test round-trip consistency: grididx -> origcoord -> grididx + test_grids = [ + # Standard precision range + DiscretizedGrid{1}(20), + # Very small range + DiscretizedGrid{1}(15; lower_bound=(1e-10,), upper_bound=(1e-9,)), + # Very large range + DiscretizedGrid{1}(12; lower_bound=(1e10,), upper_bound=(1e11,)), + # Negative range + DiscretizedGrid{1}(18; lower_bound=(-1e5,), upper_bound=(-1e4,)), + # Asymmetric around zero + DiscretizedGrid{1}(16; lower_bound=(-1e-6,), upper_bound=(1e-5,)), + ] + + for grid in test_grids + max_idx = 2^QuanticsGrids.grid_Rs(grid)[1] + # Test a representative sample of indices including edge cases + test_indices = unique([1, 2, 3, max_idx ÷ 4, max_idx ÷ 2, 3 * max_idx ÷ 4, max_idx - 2, max_idx - 1, max_idx]) + + for idx in test_indices + # Forward: grididx -> origcoord -> grididx + coord = grididx_to_origcoord(grid, idx) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + @test recovered_idx == idx + end + end +end + +@testitem "floating point precision - coordinates near but not exactly at boundaries" begin + using QuanticsGrids + R = 15 + grid = DiscretizedGrid{1}(R) + + # Test coordinates that are very close to boundaries but not exactly at them + eps_values = [eps(Float64), 2 * eps(Float64), 1e-15, 1e-14, 1e-13, 1e-12] + + lower = QuanticsGrids.lower_bound(grid) + upper = QuanticsGrids.upper_bound(grid) + + for eps_val in eps_values + # Just above lower bound + coord_above_lower = lower + eps_val + if coord_above_lower <= upper + idx = QuanticsGrids.origcoord_to_grididx(grid, coord_above_lower) + @test idx >= 1 + end + + # Just below upper bound + coord_below_upper = upper - eps_val + if coord_below_upper >= lower + idx = QuanticsGrids.origcoord_to_grididx(grid, coord_below_upper) + @test idx <= 2^R + end + end +end + +@testitem "floating point precision - multidimensional edge cases" begin + using QuanticsGrids + # Test 2D grids where floating point errors might compound + R = 12 + grid = DiscretizedGrid((R, R); + lower_bound=(0.0, -1.0), + upper_bound=(1.0, 1.0)) + + step = QuanticsGrids.grid_step(grid) + max_indices = (2^R, 2^R) + + # Test corner cases + corner_indices = [(1, 1), (1, max_indices[2]), (max_indices[1], 1), max_indices] + + for corner_idx in corner_indices + coord = grididx_to_origcoord(grid, corner_idx) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + @test recovered_idx == corner_idx + end + + # Test coordinates computed with potential precision issues + for _ in 1:50 + # Generate random indices + idx = (rand(1:2^R), rand(1:2^R)) + + # Convert to coordinates and back + coord = grididx_to_origcoord(grid, idx) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + + @test recovered_idx == idx + end +end + +@testitem "floating point precision - extreme coordinate ranges" begin + using QuanticsGrids + # Test grids with extreme coordinate ranges that stress floating point arithmetic + + test_cases = [ + # Very small positive range + (1e-16, 1e-15, 10), + # Very large range + (1e15, 1e16, 8), + # Range crossing zero with very small values + (-1e-15, 1e-15, 12), + # Asymmetric range with vastly different magnitudes + (-1e10, 1e-8, 10), + # Range where step size is close to machine epsilon + (0.0, 1e-14, 8), + ] + + for (lower, upper, R) in test_cases + grid = DiscretizedGrid{1}(R; lower_bound=(lower,), upper_bound=(upper,)) + + # Test boundary indices + boundary_indices = [1, 2, 2^R - 1, 2^R] + + for idx in boundary_indices + coord = grididx_to_origcoord(grid, idx) + + # Check that coordinate is within expected bounds + @test QuanticsGrids.lower_bound(grid) <= coord <= QuanticsGrids.upper_bound(grid) + + # Test round-trip + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + @test recovered_idx == idx + end + end +end + +@testitem "floating point precision - step size edge cases" begin + using QuanticsGrids + # Test cases where step size computation might be problematic + + # Case 1: Step size that's not exactly representable in Float64 + grid1 = DiscretizedGrid{1}(10; lower_bound=(0.0,), upper_bound=(1.0 / 3.0,)) + step1 = QuanticsGrids.grid_step(grid1) + + # Test that step size calculations are consistent + @test step1 ≈ (1.0 / 3.0) / 2^10 + + for i in 1:2^10 + coord = grididx_to_origcoord(grid1, i) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid1, coord) + @test recovered_idx == i + end + + # Case 2: Step size involving irrational numbers + sqrt2 = sqrt(2.0) + grid2 = DiscretizedGrid{1}(8; lower_bound=(0.0,), upper_bound=(sqrt2,)) + + for i in [1, 2^4, 2^8] + coord = grididx_to_origcoord(grid2, i) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid2, coord) + @test recovered_idx == i + end + + # Case 3: Very small step size that approaches machine epsilon + grid3 = DiscretizedGrid{1}(4; lower_bound=(0.0,), upper_bound=(eps(Float64) * 100,)) + + for i in 1:2^4 + coord = grididx_to_origcoord(grid3, i) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid3, coord) + @test recovered_idx == i + end +end + +@testitem "floating point precision - boundary coordinate generation" begin + using QuanticsGrids + # Test coordinates that are generated by arithmetic operations and might not be exact + R = 20 + grid = DiscretizedGrid{1}(R) + + step = QuanticsGrids.grid_step(grid) + + # Test coordinates generated by accumulative addition (potential error accumulation) + coord = 0.0 + for i in 1:min(1000, 2^R) + expected_idx = i + actual_idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + + # Allow for small errors due to accumulation + @test abs(actual_idx - expected_idx) <= 1 + + global coord += step + if coord > QuanticsGrids.upper_bound(grid) + break + end + end +end + +@testitem "floating point precision - clamp behavior verification" begin + using QuanticsGrids + # Test that the clamping behavior works correctly for edge cases + R = 8 + grid = DiscretizedGrid{1}(R) + + # Test coordinates that would round to indices outside valid range + step = QuanticsGrids.grid_step(grid) + + # Coordinate slightly before first grid point (should clamp to 1) + coord_before = -step / 2 + if coord_before >= QuanticsGrids.lower_bound(grid) + idx = QuanticsGrids.origcoord_to_grididx(grid, coord_before) + @test idx == 1 + end + + # Coordinate slightly after last grid point (should clamp to max index) + last_coord = grididx_to_origcoord(grid, 2^R) + coord_after = last_coord + step / 2 + if coord_after <= QuanticsGrids.upper_bound(grid) + idx = QuanticsGrids.origcoord_to_grididx(grid, coord_after) + @test idx == 2^R + end + + # Test that coordinates exactly at grid_max map to the last index + grid_max_coord = QuanticsGrids.grid_max(grid) + idx_at_max = QuanticsGrids.origcoord_to_grididx(grid, grid_max_coord) + @test idx_at_max == 2^R +end + +@testitem "floating point precision - includeendpoint behavior edge cases" begin + using QuanticsGrids + # Test floating point behavior with includeendpoint=true vs false + R = 10 + + grid_without = DiscretizedGrid{1}(R; includeendpoint=false) + grid_with = DiscretizedGrid{1}(R; includeendpoint=true) + + # Test that the upper bound is handled differently + upper_bound = 1.0 + + # Without includeendpoint: coordinate at upper_bound should be valid but map to last index + if upper_bound <= QuanticsGrids.upper_bound(grid_without) + idx_without = QuanticsGrids.origcoord_to_grididx(grid_without, upper_bound) + @test idx_without == 2^R + end + + # With includeendpoint: coordinate at upper_bound should map to last index + idx_with = QuanticsGrids.origcoord_to_grididx(grid_with, upper_bound) + @test idx_with == 2^R + + # Test round-trip consistency for both cases + for grid in [grid_without, grid_with] + step = QuanticsGrids.grid_step(grid) + + # Test coordinates near the upper boundary + test_coords = [ + QuanticsGrids.grid_max(grid), + QuanticsGrids.grid_max(grid) - step / 10, + QuanticsGrids.grid_max(grid) + step / 10, # This might be out of bounds + ] + + for coord in test_coords + if QuanticsGrids.lower_bound(grid) <= coord <= QuanticsGrids.upper_bound(grid) + idx = QuanticsGrids.origcoord_to_grididx(grid, coord) + recovered_coord = grididx_to_origcoord(grid, idx) + recovered_idx = QuanticsGrids.origcoord_to_grididx(grid, recovered_coord) + @test recovered_idx == idx + end + end + end +end + +@testitem "floating point precision - mixed precision operations" begin + using QuanticsGrids + # Test scenarios where different precision operations might interact + + # Create grid with parameters that don't divide evenly + grid = DiscretizedGrid{1}(12; lower_bound=(0.1,), upper_bound=(0.9,)) + + step = QuanticsGrids.grid_step(grid) + + # Test coordinates calculated using different methods + for i in [1, 100, 1000, 2^12] + # Method 1: Direct calculation + coord1 = 0.1 + (i - 1) * step + + # Method 2: Via grid function + coord2 = grididx_to_origcoord(grid, i) + + # Both should give same results (within floating point precision) + @test coord1 ≈ coord2 + + # Both should map back to same index + if coord1 <= QuanticsGrids.upper_bound(grid) + idx1 = QuanticsGrids.origcoord_to_grididx(grid, coord1) + @test idx1 == i + end + + idx2 = QuanticsGrids.origcoord_to_grididx(grid, coord2) + @test idx2 == i + end +end diff --git a/test/origcoord_tests.jl b/test/origcoord_tests.jl new file mode 100644 index 0000000..ed1d5d2 --- /dev/null +++ b/test/origcoord_tests.jl @@ -0,0 +1,439 @@ +@testitem "origcoord functions - basic 1D grid" begin + grid = DiscretizedGrid((4,); lower_bound=(0.0,), upper_bound=(1.0,)) + + # Test grididx <-> origcoord conversion + @test grididx_to_origcoord(grid, (1,)) == 0.0 + @test grididx_to_origcoord(grid, (2^4,)) ≈ 1.0 - 1.0 / 2^4 + + # Test origcoord <-> grididx conversion + @test QuanticsGrids.origcoord_to_grididx(grid, 0.0) == 1 + @test QuanticsGrids.origcoord_to_grididx(grid, 0.5) == 2^3 + 1 # 0b1000 + 1 + + # Test round-trip conversions + for grididx in [(1,), (5,), (16,)] + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == only(grididx) + end + + # Test quantics <-> origcoord conversions (should work through existing functions) + quantics = [1, 1, 1, 2] # should correspond to grididx = 2 + expected_origcoord = 1.0 / 2^4 + @test quantics_to_origcoord(grid, quantics) ≈ expected_origcoord + @test origcoord_to_quantics(grid, expected_origcoord) == quantics +end + +@testitem "origcoord functions - 2D grid" begin + grid = DiscretizedGrid((3, 4); lower_bound=(0.0, 1.0), upper_bound=(2.0, 3.0)) + + # Test boundary values + @test grididx_to_origcoord(grid, (1, 1)) == (0.0, 1.0) + @test all(isapprox.(grididx_to_origcoord(grid, (2^3, 2^4)), (2.0 - 2.0 / 2^3, 3.0 - 2.0 / 2^4))) + + # Test middle values + mid_grididx = (2^2 + 1, 2^3 + 1) # Middle of each dimension + expected_origcoord = (0.0 + 2.0 * 2^2 / 2^3, 1.0 + 2.0 * 2^3 / 2^4) + @test all(isapprox.(grididx_to_origcoord(grid, mid_grididx), expected_origcoord)) + + # Test round-trip conversions + for _ in 1:20 + grididx = (rand(1:2^3), rand(1:2^4)) + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end +end + +@testitem "origcoord functions - different base" begin + base = 3 + grid = DiscretizedGrid((2, 3); base=base, lower_bound=(-1.0, 0.0), upper_bound=(1.0, 6.0)) + + # Test boundary values + @test grididx_to_origcoord(grid, (1, 1)) == (-1.0, 0.0) + @test all(isapprox.(grididx_to_origcoord(grid, (base^2, base^3)), (1.0 - 2.0 / base^2, 6.0 - 6.0 / base^3))) + + # Test round-trip conversions + for _ in 1:20 + grididx = (rand(1:base^2), rand(1:base^3)) + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end + + # Test quantics <-> origcoord conversions + for _ in 1:20 + quantics = rand(1:base, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + origcoord = grididx_to_origcoord(grid, grididx) + + @test all(isapprox.(quantics_to_origcoord(grid, quantics), origcoord)) + @test origcoord_to_quantics(grid, origcoord) == quantics + end +end + +@testitem "origcoord functions - boundary checking" begin + grid = DiscretizedGrid((2, 2); lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0)) + + # Test coordinates within bounds + @test QuanticsGrids.origcoord_to_grididx(grid, (0.0, 0.0)) == (1, 1) + @test QuanticsGrids.origcoord_to_grididx(grid, (0.5, 0.5)) == (3, 3) + @test QuanticsGrids.origcoord_to_grididx(grid, (0.999, 0.999)) == (4, 4) + + # Test coordinates outside bounds should throw errors + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (-0.1, 0.5)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (0.5, -0.1)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (1.1, 0.5)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (0.5, 1.1)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (-0.1, -0.1)) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, (1.1, 1.1)) +end + +@testitem "origcoord functions - fused indices grid" begin + # Test with a more complex grid that has fused indices + grid = DiscretizedGrid( + (:x, :y), + [[(:x, 2), (:y, 1)], [(:x, 1)], [(:y, 2)]]; + lower_bound=(0.0, -1.0), + upper_bound=(4.0, 1.0) + ) + + # Test some conversions + for _ in 1:50 + quantics = [rand(1:4), rand(1:2), rand(1:2)] # Based on sitedim values + grididx = quantics_to_grididx(grid, quantics) + origcoord = grididx_to_origcoord(grid, grididx) + + # Test round-trip conversion + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + @test all(isapprox.(quantics_to_origcoord(grid, quantics), origcoord)) + @test origcoord_to_quantics(grid, origcoord) == quantics + end +end + +@testitem "origcoord functions - stress test" begin + # Test high-dimensional grids + grid = DiscretizedGrid((2, 3, 2, 4); base=2, + lower_bound=(0.0, -2.0, 1.0, -5.0), + upper_bound=(1.0, 2.0, 3.0, 5.0)) + + for _ in 1:30 + # Test random grid indices + grididx = (rand(1:2^2), rand(1:2^3), rand(1:2^2), rand(1:2^4)) + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + + # Test random quantics + quantics = rand(1:2, length(grid)) + grididx2 = quantics_to_grididx(grid, quantics) + origcoord2 = grididx_to_origcoord(grid, grididx2) + + @test all(isapprox.(quantics_to_origcoord(grid, quantics), origcoord2)) + @test origcoord_to_quantics(grid, origcoord2) == quantics + end +end + +@testitem "edge cases - numerical precision and floating point" begin + # Test with very small grid steps that might cause floating point issues + grid = DiscretizedGrid((30, 25); + lower_bound=(1e-15, -1e-15), + upper_bound=(1e-14, 1e-14)) + + # Test conversions near boundaries where floating point precision matters + for _ in 1:20 + grididx = (rand(1:2^30), rand(1:2^25)) + origcoord = grididx_to_origcoord(grid, grididx) + recovered_grididx = QuanticsGrids.origcoord_to_grididx(grid, origcoord) + @test recovered_grididx == grididx + end + + # Test coordinates very close to boundaries + eps_coord_x = 1e-15 + 1e-20 # Just above lower bound + eps_coord_y = 1e-14 - 1e-20 # Just below upper bound + @test QuanticsGrids.origcoord_to_grididx(grid, (eps_coord_x, eps_coord_y)) isa Tuple +end + +@testitem "edge cases - maximum R values and large numbers" begin + # Test with largest R values that still fit in computation + # R=60 gives 2^60 ≈ 1.15e18, near but below Int64 max (9.22e18) + large_R = 58 # Be conservative to avoid overflow + grid = DiscretizedGrid((large_R,)) + + # Test extreme indices + min_quantics = ones(Int, large_R) + max_quantics = fill(2, large_R) + + @test quantics_to_grididx(grid, min_quantics) == 1 + @test quantics_to_grididx(grid, max_quantics) == 2^large_R + @test grididx_to_quantics(grid, (1,)) == min_quantics + @test grididx_to_quantics(grid, (2^large_R,)) == max_quantics + + # Test some random middle values + for _ in 1:10 + quantics = rand(1:2, large_R) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + end +end + +@testitem "edge cases - degenerate grids with R=0" begin + # Grid with all R=0 (should only have one grid point each) + grid = DiscretizedGrid((0, 0, 0)) + + # Only valid grid index should be (1, 1, 1) + @test grididx_to_quantics(grid, (1, 1, 1)) == Int[] + @test quantics_to_grididx(grid, Int[]) == (1, 1, 1) + + # Test origcoord functions + @test grididx_to_origcoord(grid, (1, 1, 1)) == (0.0, 0.0, 0.0) + @test QuanticsGrids.origcoord_to_grididx(grid, (0.0, 0.0, 0.0)) == (1, 1, 1) +end + +@testitem "edge cases - mixed R values including zeros" begin + # Grid with mix of zero and non-zero R values + grid = DiscretizedGrid((0, 5, 0, 3, 0); unfoldingscheme=:interleaved) + + # Test conversion with valid quantics length + quantics = [1, 2, 1, 2, 1, 2, 1, 2] # Only bits for dimensions with R>0 + grididx = quantics_to_grididx(grid, quantics) + @test grididx[1] == 1 # R=0 dimensions should always be 1 + @test grididx[3] == 1 # R=0 dimensions should always be 1 + @test grididx[5] == 1 # R=0 dimensions should always be 1 + @test grididx_to_quantics(grid, grididx) == quantics + + # Test boundary conditions + for _ in 1:20 + valid_quantics = rand(1:2, length(grid)) + grididx_local = quantics_to_grididx(grid, valid_quantics) + @test grididx_local[1] == 1 + @test grididx_local[3] == 1 + @test grididx_local[5] == 1 + @test 1 <= grididx_local[2] <= 2^5 + @test 1 <= grididx_local[4] <= 2^3 + end +end + +@testitem "edge cases - extreme base values" begin + # Test with maximum reasonable base + max_base = 16 + grid = DiscretizedGrid((3, 4); base=max_base) + + for _ in 1:20 + quantics = rand(1:max_base, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + @test all(1 .<= grididx .<= (max_base .^ QuanticsGrids.grid_Rs(grid))) + end + + # Test with base=2 (edge case, minimum valid base) + grid2 = DiscretizedGrid((10, 8); base=2) + for _ in 1:20 + quantics = rand(1:2, length(grid2)) + grididx = quantics_to_grididx(grid2, quantics) + @test grididx_to_quantics(grid2, grididx) == quantics + end +end + +@testitem "edge cases - extreme coordinate ranges" begin + # Test with very large coordinate ranges + grid = DiscretizedGrid((10, 8); + lower_bound=(-1e10, -1e15), + upper_bound=(1e10, 1e15)) + + for _ in 1:20 + grididx = (rand(1:2^10), rand(1:2^8)) + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end + + # Test with very small coordinate ranges + grid2 = DiscretizedGrid((5, 6); + lower_bound=(1e-12, -1e-12), + upper_bound=(1e-12 + 1e-15, -1e-12 + 1e-15)) + + for _ in 1:10 + grididx = (rand(1:2^5), rand(1:2^6)) + origcoord = grididx_to_origcoord(grid2, grididx) + recovered = QuanticsGrids.origcoord_to_grididx(grid2, origcoord) + @test recovered == grididx + end +end + +@testitem "edge cases - asymmetric ranges with negative bounds" begin + # Test grids with very asymmetric and negative bounds + grid = DiscretizedGrid((4, 6, 3); + lower_bound=(-1000.0, -0.001, -1e6), + upper_bound=(-999.0, 0.001, 1e6)) + + for _ in 1:30 + quantics = rand(1:2, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + origcoord = grididx_to_origcoord(grid, grididx) + + # Verify coordinates are within bounds + @test all(QuanticsGrids.grid_min(grid) .<= origcoord .<= QuanticsGrids.grid_max(grid)) + + # Test round-trip + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + @test all(isapprox.(quantics_to_origcoord(grid, quantics), origcoord)) + end +end + +@testitem "stress test - maximum complexity fused indices" begin + # Create maximally complex fused index structure + # All variables distributed unevenly across sites + grid = DiscretizedGrid( + (:a, :b, :c, :d, :e, :f, :g, :h), + [ + [(:a, 8), (:b, 7), (:c, 6), (:d, 5)], # 4 vars, site dim = base^4 + [(:e, 4), (:f, 3)], # 2 vars, site dim = base^2 + [(:g, 2)], # 1 var + [(:h, 1)], # 1 var + [(:a, 7), (:c, 5), (:e, 3), (:g, 1)], # 4 vars, complex pattern + [(:b, 6), (:d, 4), (:f, 2)], # 3 vars + [(:a, 6), (:h, 2)], # 2 vars from distant indices + [(:b, 5), (:c, 4), (:d, 3)], # 3 vars + [(:e, 2), (:f, 1)], # 2 vars + [(:a, 5), (:g, 3)], # 2 vars, distant indices + [(:b, 4), (:c, 3), (:d, 2), (:h, 3)], # 4 vars + [(:a, 4), (:e, 1)], # 2 vars + [(:b, 3), (:f, 4)], # 2 vars + [(:a, 3), (:c, 2)], # 2 vars + [(:b, 2), (:d, 1)], # 2 vars + [(:a, 2)], # 1 var + [(:b, 1)], # 1 var + [(:a, 1)], # 1 var + [(:c, 1)] # 1 var + ]; + base=2 + ) + + # Stress test with many conversions + for _ in 1:100 + # Generate valid quantics based on site dimensions + quantics = [ + rand(1:16), # site 1: 2^4 + rand(1:4), # site 2: 2^2 + rand(1:2), # site 3: 2^1 + rand(1:2), # site 4: 2^1 + rand(1:16), # site 5: 2^4 + rand(1:8), # site 6: 2^3 + rand(1:4), # site 7: 2^2 + rand(1:8), # site 8: 2^3 + rand(1:4), # site 9: 2^2 + rand(1:4), # site 10: 2^2 + rand(1:16), # site 11: 2^4 + rand(1:4), # site 12: 2^2 + rand(1:4), # site 13: 2^2 + rand(1:4), # site 14: 2^2 + rand(1:4), # site 15: 2^2 + rand(1:2), # site 16: 2^1 + rand(1:2), # site 17: 2^1 + rand(1:2), # site 18: 2^1 + rand(1:2) # site 19: 2^1 + ] + + grididx = quantics_to_grididx(grid, quantics) + recovered = grididx_to_quantics(grid, grididx) + @test recovered == quantics + + # Test origcoord functions too + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end +end + +@testitem "boundary stress test - coordinates exactly on boundaries" begin + grid = DiscretizedGrid((8, 6); + lower_bound=(-2.0, 3.0), + upper_bound=(5.0, 8.0)) + + # Test coordinates exactly at lower bounds + @test QuanticsGrids.origcoord_to_grididx(grid, (-2.0, 3.0)) == (1, 1) + + # Test coordinates exactly at effective upper bounds (grid_max) + grid_max_coord = (-2.0 + 7.0 * (2^8 - 1) / 2^8, 3.0 + 5.0 * (2^6 - 1) / 2^6) + max_grididx = (2^8, 2^6) + @test QuanticsGrids.origcoord_to_grididx(grid, grid_max_coord) == max_grididx + + # Test coordinates just inside boundaries + eps = 1e-12 + inside_lower = (-2.0 + eps, 3.0 + eps) + inside_upper = (grid_max_coord[1] - eps, grid_max_coord[2] - eps) + + @test QuanticsGrids.origcoord_to_grididx(grid, inside_lower) isa Tuple + @test QuanticsGrids.origcoord_to_grididx(grid, inside_upper) isa Tuple + + # Test coordinates just outside boundaries should error + outside_lower = (-2.0 - eps, 3.0 - eps) + outside_upper = (5.0 + eps, 8.0 + eps) + + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, outside_lower) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(grid, outside_upper) +end + +@testitem "performance stress test - rapid conversions" begin + # Test performance with moderately sized grids and many conversions + grid = DiscretizedGrid((12, 10, 8, 6); base=3) + + # Pre-generate test data to avoid timing allocation + n_tests = 1000 + test_quantics = [rand(1:3, length(grid)) for _ in 1:n_tests] + test_grididx = [(rand(1:3^12), rand(1:3^10), rand(1:3^8), rand(1:3^6)) for _ in 1:n_tests] + + # Test quantics -> grididx -> quantics round trips + for i in 1:n_tests + quantics = test_quantics[i] + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + end + + # Test grididx -> quantics -> grididx round trips + for i in 1:n_tests + grididx = test_grididx[i] + quantics = grididx_to_quantics(grid, grididx) + @test quantics_to_grididx(grid, quantics) == grididx + end + + # Test origcoord conversions + for i in 1:min(n_tests, 200) # Fewer origcoord tests since they're slower + grididx = test_grididx[i] + origcoord = grididx_to_origcoord(grid, grididx) + @test QuanticsGrids.origcoord_to_grididx(grid, origcoord) == grididx + end +end + +@testitem "R=0 dimension behavior - 2D grid with proper semantics" begin + # Test a 2D grid where one dimension has R=0 and the other has R=3 + grid = DiscretizedGrid((0, 3); lower_bound=(-1.0, 2.0), upper_bound=(1.0, 6.0)) + + # R=0 dimension should only accept grid index 1 + @test_nowarn grididx_to_quantics(grid, (1, 1)) + @test_throws DomainError grididx_to_quantics(grid, (2, 1)) + + # All coordinates in R=0 dimension should map to grid index 1 + @test origcoord_to_grididx(grid, (-0.5, 2.0))[1] == 1 + @test origcoord_to_grididx(grid, (0.5, 2.0))[1] == 1 + + # Grid index 1 in R=0 dimension + coord = grididx_to_origcoord(grid, (1, 1)) + @test coord[1] ≈ -1. + + # Quantics should only contain indices for R>0 dimensions + @test length(grididx_to_quantics(grid, (1, 1))) == 3 # Only for R=3 dimension + + # Test grid bounds functions + @test QuanticsGrids.grid_min(grid)[1] ≈ -1.0 # R=0 dimension: midpoint + @test QuanticsGrids.grid_max(grid)[1] ≈ -1.0 # R=0 dimension: midpoint + @test QuanticsGrids.grid_min(grid)[2] ≈ 2.0 # R=3 dimension: lower bound + + # Test grid_origcoords for R=0 dimension + coords_r0 = QuanticsGrids.grid_origcoords(grid, 1) + @test length(coords_r0) == 1 # Only one point for R=0 + @test coords_r0[1] ≈ -1.0 # Should be the midpoint + + coords_r3 = QuanticsGrids.grid_origcoords(grid, 2) + @test length(coords_r3) == 8 # 2^3 points for R=3 + + # Test round-trip consistency + quantics = grididx_to_quantics(grid, (1, 5)) + @test quantics_to_grididx(grid, quantics) == (1, 5) + + @test_throws ArgumentError DiscretizedGrid((3, 0); includeendpoint=true) +end diff --git a/test/quantics_test.jl b/test/quantics_test.jl deleted file mode 100644 index f09f82f..0000000 --- a/test/quantics_test.jl +++ /dev/null @@ -1,103 +0,0 @@ -@testitem "quantics.jl" begin - import QuanticsGrids: fuse_dimensions, unfuse_dimensions - import QuanticsGrids: quantics_to_index_fused, index_to_quantics - import QuanticsGrids: interleave_dimensions, deinterleave_dimensions - - @testset "quantics representation" begin - @testset "fuse_dimensions" begin - x = [1, 1] - y = [2, 2] - fused = fuse_dimensions(x, y) - @test unfuse_dimensions(fused, 2)[1] == x - @test unfuse_dimensions(fused, 2)[2] == y - end - - @testset "quantics to index, 1d" begin - B = 2 - @test quantics_to_index_fused([1, 1, 1, 1]; base = B) == (1,) - @test quantics_to_index_fused([1, 1, 1, 2]; base = B) == (2,) - @test quantics_to_index_fused([1, 1, 2, 1]; base = B) == (3,) - @test quantics_to_index_fused([1, 1, 2, 2]; base = B) == (4,) - @test quantics_to_index_fused([2, 1, 1, 1]; base = B) == (9,) - end - - @testset "quantics to index, 1d (general power base)" for B in [2, 4] - d = 1 - R = 4 # number of digits - - index_reconst = Int[] - for index = 1:B^R - digitlist_ = QuanticsGrids.index_to_quantics(index; numdigits = R, base = B) - push!( - index_reconst, - only(quantics_to_index_fused(digitlist_; base = B, dims = Val(d))), - ) - end - - @test collect(1:B^R) == index_reconst - end - - @testset "quantics to index, 2d, base=3" begin - base = 3 - dim = 2 - - # X_i = Fused quantics index at i (1 <= i <= R) - # x_i = quantics index for the first variable at i (1 <= i <= R) - # y_i = quantics index for the second variable at i (1 <= i <= R) - # - # X_i = (x_i-1) + (base) * (y_i-1) + 1 (column major) - @test quantics_to_index_fused([1, 1]; base = base, dims = Val(dim)) == (1, 1) - @test quantics_to_index_fused([1, 2]; base = base, dims = Val(dim)) == (2, 1) - @test quantics_to_index_fused([1, 3]; base = base, dims = Val(dim)) == (3, 1) - @test quantics_to_index_fused([1, 4]; base = base, dims = Val(dim)) == (1, 2) - @test quantics_to_index_fused([1, 5]; base = base, dims = Val(dim)) == (2, 2) - @test quantics_to_index_fused([1, 6]; base = base, dims = Val(dim)) == (3, 2) - @test quantics_to_index_fused([1, 7]; base = base, dims = Val(dim)) == (1, 3) - @test quantics_to_index_fused([1, 8]; base = base, dims = Val(dim)) == (2, 3) - @test quantics_to_index_fused([1, 9]; base = base, dims = Val(dim)) == (3, 3) - @test quantics_to_index_fused([2, 1]; base = base, dims = Val(dim)) == (4, 1) - end - - @testset "quantics back-and-forth, 2d" begin - base = 2 - dim = 2 - R = 2 - - for j = 1:base^R, i = 1:base^R - index = (i, j) - - digitlist1 = Vector{Int}(undef, R) - QuanticsGrids.index_to_quantics_fused!(digitlist1, index; base = base) - digitlist2 = - QuanticsGrids.index_to_quantics_fused(index; numdigits = R, base = base) - @test digitlist1 == digitlist2 - - index_reconst = QuanticsGrids.quantics_to_index_fused( - digitlist1; - base = base, - dims = Val(dim), - ) - @test index == index_reconst - end - end - - @testset "interleave dimensions" begin - @test [1, 1, 1, 1] == interleave_dimensions([1, 1, 1, 1]) - @test [1, 1, 1, 1, 1, 1, 1, 1] == - interleave_dimensions([1, 1, 1, 1], [1, 1, 1, 1]) - @test [1, 2, 1, 3, 1, 4, 1, 5] == - interleave_dimensions([1, 1, 1, 1], [2, 3, 4, 5]) - @test [1, 2, 11, 1, 3, 12, 1, 4, 13, 1, 5, 14] == - interleave_dimensions([1, 1, 1, 1], [2, 3, 4, 5], [11, 12, 13, 14]) - end - - @testset "deinterleave dimensions" begin - @test deinterleave_dimensions([1, 1, 1, 1], 1) == [[1, 1, 1, 1]] - @test deinterleave_dimensions([1, 1, 1, 1], 2) == [[1, 1], [1, 1]] - @test deinterleave_dimensions([1, 2, 1, 3, 1, 4, 1, 5], 2) == - [[1, 1, 1, 1], [2, 3, 4, 5]] - @test deinterleave_dimensions([1, 2, 11, 1, 3, 12, 1, 4, 13, 1, 5, 14], 3) == - [[1, 1, 1, 1], [2, 3, 4, 5], [11, 12, 13, 14]] - end - end -end diff --git a/test/quantics_tests.jl b/test/quantics_tests.jl new file mode 100644 index 0000000..8e4b79c --- /dev/null +++ b/test/quantics_tests.jl @@ -0,0 +1,303 @@ +@testitem "constructor, square grid" begin + grid = DiscretizedGrid{2}(10; unfoldingscheme=:interleaved) + @test QuanticsGrids.grid_Rs(grid) == (10, 10) + @test only(QuanticsGrids.grid_indextable(grid)[1])[1] == Symbol(1) + @test only(QuanticsGrids.grid_indextable(grid)[1])[2] == 1 + @test only(QuanticsGrids.grid_indextable(grid)[2])[1] == Symbol(2) + @test only(QuanticsGrids.grid_indextable(grid)[2])[2] == 1 +end + +@testitem "constructor, rectangular grid" begin + grid = DiscretizedGrid((3, 5); unfoldingscheme=:interleaved) + @test only(QuanticsGrids.grid_indextable(grid)[6])[1] == Symbol(2) + @test only(QuanticsGrids.grid_indextable(grid)[6])[2] == 3 + @test only(QuanticsGrids.grid_indextable(grid)[7])[1] == Symbol(2) + @test only(QuanticsGrids.grid_indextable(grid)[7])[2] == 4 + @test only(QuanticsGrids.grid_indextable(grid)[8])[1] == Symbol(2) + @test only(QuanticsGrids.grid_indextable(grid)[8])[2] == 5 +end + +@testitem "quantics_to_grididx, rectangular grid" begin + grid = DiscretizedGrid((3, 5); unfoldingscheme=:interleaved) + @test quantics_to_grididx(grid, [1, 2, 1, 2, 1, 2, 1, 2]) == (1, 30) +end + +@testitem "grididx_to_quantics, rectangular grid" begin + grid = DiscretizedGrid((3, 5); unfoldingscheme=:interleaved) + @test grididx_to_quantics(grid, (1, 30)) == [1, 2, 1, 2, 1, 2, 1, 2] +end + +@testitem "quantics_to_grididx ∘ grididx_to_quantics == identity" begin + grid = DiscretizedGrid((5, 3, 17); base=13) + for _ in 1:100 + grididx = ntuple(d -> rand(1:QuanticsGrids.grid_Rs(grid)[d]), ndims(grid)) + @test quantics_to_grididx(grid, grididx_to_quantics(grid, grididx)) == grididx + end +end + +@testitem "grididx_to_quantics ∘ quantics_to_grididx == identity" begin + grid = DiscretizedGrid((48, 31, 62)) + for _ in 1:100 + quantics = rand(1:2, length(grid)) + @test grididx_to_quantics(grid, quantics_to_grididx(grid, quantics)) == quantics + end +end + +@testitem "grididx_to_quantics ∘ quantics_to_grididx == identity, base != 2" begin + base = 7 + grid = DiscretizedGrid((22, 9, 14); base) + for _ in 1:100 + quantics = rand(1:base, length(grid)) + @test grididx_to_quantics(grid, quantics_to_grididx(grid, quantics)) == quantics + end +end + +@testitem "ctor from indextable" begin + grid = DiscretizedGrid((:a, :b, :c), [[(:a, 4)], [(:a, 3)], [(:a, 2)], [(:a, 1)], [(:b, 1)], [(:b, 2)], [(:b, 3)], [(:c, 1)], [(:c, 2)], [(:c, 3)]]) + @test QuanticsGrids.grid_Rs(grid) == (4, 3, 3) + @test QuanticsGrids.lower_bound(grid) == (0.0, 0.0, 0.0) + @test QuanticsGrids.upper_bound(grid) == (1.0, 1.0, 1.0) + @test grid.discretegrid.variablenames == (:a, :b, :c) +end + +@testitem "ctor from indextable, quantics <-> grididx" begin + grid = DiscretizedGrid((:a, :b, :c), [[(:a, 4)], [(:a, 3)], [(:a, 2)], [(:a, 1)], [(:b, 1)], [(:b, 2)], [(:b, 3)], [(:c, 1)], [(:c, 2)], [(:c, 3)]]) + @test quantics_to_grididx(grid, [1, 2, 1, 2, 1, 2, 1, 2, 1, 2]) == (11, 3, 6) + @test grididx_to_quantics(grid, (11, 3, 6)) == [1, 2, 1, 2, 1, 2, 1, 2, 1, 2] +end + +@testitem "ctor from indextable, quantics <-> grididx, fused indices" begin + grid = DiscretizedGrid((:a, :b, :c, :d), [[(:a, 4)], [(:a, 3)], [(:a, 2)], [(:a, 1)], [(:b, 1), (:d, 1)], [(:b, 2)], [(:b, 3)], [(:c, 1), (:d, 2)], [(:c, 2)], [(:c, 3)]]) + @test quantics_to_grididx(grid, [1, 2, 1, 2, 2, 2, 1, 4, 1, 2]) == (11, 3, 6, 4) + @test grididx_to_quantics(grid, (11, 3, 6, 4)) == [1, 2, 1, 2, 2, 2, 1, 4, 1, 2] +end + +@testitem "challenging tests - extreme edge cases" begin + # Test with minimum valid grididx (all 1s) + grid = DiscretizedGrid((10, 5, 8); unfoldingscheme=:interleaved) + min_grididx = (1, 1, 1) + quantics = grididx_to_quantics(grid, min_grididx) + @test all(q -> q == 1, quantics) + @test quantics_to_grididx(grid, quantics) == min_grididx + + # Test with maximum valid grididx + max_grididx = (2^10, 2^5, 2^8) + quantics = grididx_to_quantics(grid, max_grididx) + @test all(q -> q == 2, quantics) + @test quantics_to_grididx(grid, quantics) == max_grididx +end + +@testitem "challenging tests - mixed bases" begin + # Test with base 3 + grid = DiscretizedGrid((4, 6, 3); base=3) + for _ in 1:50 + quantics = rand(1:3, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + @test all(1 .<= grididx .<= (3 .^ QuanticsGrids.grid_Rs(grid))) + end + + # Test with base 5 + grid = DiscretizedGrid((3, 4); base=5, unfoldingscheme=:interleaved) + for _ in 1:50 + grididx = (rand(1:5^3), rand(1:5^4)) + quantics = grididx_to_quantics(grid, grididx) + @test quantics_to_grididx(grid, quantics) == grididx + @test all(1 .<= quantics .<= 5) + end +end + +@testitem "challenging tests - complex fused indices" begin + # Multiple variables fused in single sites + grid = DiscretizedGrid( + (:x, :y, :z), + [ + [(:x, 3), (:y, 2), (:z, 1)], # 3 variables in one site + [(:x, 2)], # single variable + [(:y, 1), (:z, 2)], # 2 variables fused + [(:x, 1)], # single variable + [(:z, 3)] # single variable + ] + ) + + # Test specific known values + @test quantics_to_grididx(grid, [8, 1, 4, 2, 1]) == (6, 4, 7) + @test grididx_to_quantics(grid, (6, 4, 7)) == [8, 1, 4, 2, 1] + + # Test random values + for _ in 1:100 + quantics = [ + rand(1:8), # site 1: base^3 = 8 possibilities + rand(1:2), # site 2: base = 2 possibilities + rand(1:4), # site 3: base^2 = 4 possibilities + rand(1:2), # site 4: base = 2 possibilities + rand(1:2) # site 5: base = 2 possibilities + ] + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + end +end + +@testitem "challenging tests - asymmetric grids" begin + # Very asymmetric grid dimensions + grid = DiscretizedGrid((20, 0, 15, 3)) + + for _ in 1:50 + # Test boundary conditions + quantics = rand(1:2, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + + # Verify grid index bounds + @test 1 <= grididx[1] <= 2^20 + @test grididx[2] == 1 # R=0 means only one possible value + @test 1 <= grididx[3] <= 2^15 + @test 1 <= grididx[4] <= 2^3 + end +end + +@testitem "challenging tests - single dimension edge cases" begin + # 1D grid with large R + grid = DiscretizedGrid((25,)) + + # Test extremes + min_quantics = ones(Int, 25) + max_quantics = fill(2, 25) + + @test quantics_to_grididx(grid, min_quantics) == 1 + @test quantics_to_grididx(grid, max_quantics) == 2^25 + @test grididx_to_quantics(grid, (1,)) == min_quantics + @test grididx_to_quantics(grid, (2^25,)) == max_quantics + + # Test middle values + for _ in 1:20 + quantics = rand(1:2, 25) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + end +end + +@testitem "challenging tests - high dimensional grids" begin + # 8D grid with moderate R values + grid = DiscretizedGrid(ntuple(i -> 4 + (i % 3), 8); base=3) + + for _ in 1:30 + quantics = rand(1:3, length(grid)) + grididx = quantics_to_grididx(grid, quantics) + @test grididx_to_quantics(grid, grididx) == quantics + + # Verify all grid indices are within bounds + for d in 1:8 + @test 1 <= grididx[d] <= 3^QuanticsGrids.grid_Rs(grid)[d] + end + end +end + +@testitem "challenging tests - stress test with complex patterns" begin + # Grid with intentionally complex index table structure + grid = DiscretizedGrid( + (:a, :b, :c, :d, :e), + [ + [(:e, 1)], # site 1 + [(:a, 5), (:c, 4)], # site 2: 2 vars fused + [(:b, 3)], # site 3 + [(:a, 4), (:b, 2), (:d, 3)], # site 4: 3 vars fused + [(:c, 3), (:e, 2)], # site 5: 2 vars fused + [(:a, 3)], # site 6 + [(:b, 1), (:d, 2)], # site 7: 2 vars fused + [(:a, 2), (:c, 2), (:d, 1), (:e, 3)], # site 8: 4 vars fused + [(:a, 1)], # site 9 + [(:c, 1)] # site 10 + ]; + base=3 + ) + + # Stress test with many random conversions + for _ in 1:200 + # Generate valid quantics vector + quantics = [ + rand(1:3), # site 1 + rand(1:9), # site 2: 3^2 = 9 + rand(1:3), # site 3 + rand(1:27), # site 4: 3^3 = 27 + rand(1:9), # site 5: 3^2 = 9 + rand(1:3), # site 6 + rand(1:9), # site 7: 3^2 = 9 + rand(1:81), # site 8: 3^4 = 81 + rand(1:3), # site 9 + rand(1:3) # site 10 + ] + + grididx = quantics_to_grididx(grid, quantics) + recovered = grididx_to_quantics(grid, grididx) + @test recovered == quantics + + # Verify grid indices are reasonable + @test all(1 .<= grididx .<= (3 .^ QuanticsGrids.grid_Rs(grid))) + end +end + +@testitem "fused ordering compatibility with DiscretizedGrid" begin + using QuanticsGrids: DiscretizedGrid + + # Test 2D case - this was the main issue discovered during benchmarking + # where coordinates were being swapped between implementations + R = 3 # Use same R for both dimensions to create equivalent grids + lower = (0.0, 0.0) + upper = (1.0, 1.0) + + # Create equivalent grids + old_grid = DiscretizedGrid{2}(R, lower, upper; unfoldingscheme=:fused) + new_grid = DiscretizedGrid((R, R); lower_bound=lower, upper_bound=upper, unfoldingscheme=:fused) + + # Test specific quantics vectors that revealed the ordering issue + test_cases = [ + [1, 2, 1], # This was giving different coordinates before the fix + [2, 1, 2], + [1, 1, 1], # Edge case: all 1s + [2, 2, 2], # Edge case: all 2s + [2, 1, 1], + [1, 2, 2] + ] + + for quantics in test_cases + # Convert to coordinates using both implementations + old_coord = quantics_to_origcoord(old_grid, quantics) + new_coord = quantics_to_origcoord(new_grid, quantics) + + # Coordinates should be identical (within floating point precision) + @test all(old_coord .≈ new_coord) + + # Also verify the round-trip consistency within each implementation + old_grididx = quantics_to_grididx(old_grid, quantics) + new_grididx = quantics_to_grididx(new_grid, quantics) + @test old_grididx == new_grididx + + # Verify quantics -> grididx -> origcoord is consistent + old_coord_via_grididx = grididx_to_origcoord(old_grid, old_grididx) + new_coord_via_grididx = grididx_to_origcoord(new_grid, new_grididx) + @test all(old_coord_via_grididx .≈ new_coord_via_grididx) + end + + # Test 3D case to ensure the fix works for higher dimensions + R_3d = 2 # Use same R for all dimensions + lower_3d = (0.0, 0.0, 0.0) + upper_3d = (1.0, 1.0, 1.0) + + old_grid_3d = DiscretizedGrid{3}(R_3d, lower_3d, upper_3d; unfoldingscheme=:fused) + new_grid_3d = DiscretizedGrid((R_3d, R_3d, R_3d); lower_bound=lower_3d, upper_bound=upper_3d, unfoldingscheme=:fused) + + # Test several 3D quantics vectors + test_cases_3d = [ + [1, 1], # Length R=2 for fused scheme + [2, 2], + [1, 2], + [2, 1] + ] + + for quantics in test_cases_3d + old_coord = quantics_to_origcoord(old_grid_3d, quantics) + new_coord = quantics_to_origcoord(new_grid_3d, quantics) + @test all(old_coord .≈ new_coord) + end +end diff --git a/test/utilities_tests.jl b/test/utilities_tests.jl new file mode 100644 index 0000000..583263d --- /dev/null +++ b/test/utilities_tests.jl @@ -0,0 +1,523 @@ +@testitem "localdimensions" begin + R = 4 + # Basic test with default base=2 + grid = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,)) + @test QuanticsGrids.localdimensions(grid) == fill(2, R) + + grid_2d = DiscretizedGrid((R, R); lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0), unfoldingscheme=:interleaved) + @test QuanticsGrids.localdimensions(grid_2d) == fill(2, 2R) + + # Test with different base values + grid_base3 = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), base=3) + @test QuanticsGrids.localdimensions(grid_base3) == fill(3, R) + + grid_base4 = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), base=4) + @test QuanticsGrids.localdimensions(grid_base4) == fill(4, R) + + # Test with custom indextable - sites with multiple quantics indices + variablenames = (:x, :y) + # Create indextable where some sites have multiple quantics indices + indextable = [ + [(:x, 1), (:y, 1)], # site 1: 2 indices + [(:x, 2)], # site 2: 1 index + [(:y, 2)], # site 3: 1 index + [(:x, 3), (:y, 3), (:x, 4)] # site 4: 3 indices + ] + + # Test with base=2: Quan ticsGrids.localdimensions should be [2^2, 2^1, 2^1, 2^3] = [4, 2, 2, 8] + grid_custom_base2 = DiscretizedGrid(variablenames, indextable; + lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0), base=2) + @test QuanticsGrids.localdimensions(grid_custom_base2) == [4, 2, 2, 8] + + # Test with base=3: Quan ticsGrids.localdimensions should be [3^2, 3^1, 3^1, 3^3] = [9, 3, 3, 27] + grid_custom_base3 = DiscretizedGrid(variablenames, indextable; + lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0), base=3) + @test QuanticsGrids.localdimensions(grid_custom_base3) == [9, 3, 3, 27] + + # Test edge case: single site with multiple indices + single_site_indextable = [ + [(:x, 1), (:x, 2), (:x, 3)] + ] + grid_single_site = DiscretizedGrid((:x,), single_site_indextable; + lower_bound=(0.0,), upper_bound=(1.0,), base=2) + @test QuanticsGrids.localdimensions(grid_single_site) == [8] # 2^3 = 8 + + # Test with mixed site sizes in higher dimension + complex_indextable = [ + [(:x, 1)], # site 1: 1 index -> base^1 + [(:y, 1), (:z, 1)], # site 2: 2 indices -> base^2 + [(:x, 2), (:y, 2), (:z, 2)], # site 3: 3 indices -> base^3 + [(:x, 3), (:y, 3)] # site 4: 2 indices -> base^2 + ] + grid_complex = DiscretizedGrid((:x, :y, :z), complex_indextable; + lower_bound=(0.0, 0.0, 0.0), upper_bound=(1.0, 1.0, 1.0), base=2) + @test QuanticsGrids.localdimensions(grid_complex) == [2, 4, 8, 4] # [2^1, 2^2, 2^3, 2^2] +end + +@testitem "quanticsfunction" begin + R = 4 + grid = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,)) + + # Test that quanticsfunction should work with DiscretizedGrid + fx(x) = exp(-x) + fq = QuanticsGrids.quanticsfunction(Float64, grid, fx) + + # Test basic functionality - should convert quantics to coordinates and apply function + @test fq(fill(1, R)) == fx(0) # quantics [1,1,1,1] should map to x=0.0 + @test fq(fill(2, R)) == fx(1 - 1 / 2^R) # quantics [2,2,2,2] should map to near x=1.0 + + # Test 2D case + grid_2d = DiscretizedGrid((R, R); lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0)) + f2d(x, y) = x + y + fq_2d = QuanticsGrids.quanticsfunction(Float64, grid_2d, f2d) + + # Test 2D functionality + quantics_origin = fill(1, length(grid_2d)) # Should map to (0.0, 0.0) + @test fq_2d(quantics_origin) == f2d(0.0, 0.0) + + # Test with custom indextable + variablenames = (:x, :y) + indextable = [ + [(:x, 1)], + [(:x, 2)], + [(:y, 1)], + [(:y, 2)] + ] + grid_custom = DiscretizedGrid(variablenames, indextable; + lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0)) + fq_custom = QuanticsGrids.quanticsfunction(Float64, grid_custom, f2d) + @test fq_custom([1, 1, 1, 1]) == f2d(0.0, 0.0) + + # CHALLENGING TEST: Complex multi-variable function with mixed site sizes and base=3 + # Create a complex indextable with different numbers of quantics per site + complex_variablenames = (:x, :y, :z) + complex_indextable = [ + [(:x, 1)], + [(:y, 1), (:z, 1)], + [(:x, 2), (:y, 2)], + [(:z, 2), (:x, 3), (:y, 3)] + ] + + # Create grid with base=3, non-unit domain + grid_complex = DiscretizedGrid(complex_variablenames, complex_indextable; + lower_bound=(-1.0, 2.0, 0.5), upper_bound=(3.0, 5.0, 2.0), base=3) + + # Define a challenging multi-variable function with cross-terms and oscillations + f_complex(x, y, z) = sin(π * x) * cos(2π * y) * exp(-z) + x^2 * y - z^3 + x * y * z + fq_complex = QuanticsGrids.quanticsfunction(Float64, grid_complex, f_complex) + + # Test corner cases: minimum coordinates (all quantics = 1) + @test fq_complex([1, 1, 1, 1]) ≈ f_complex(-1.0, 2.0, 0.5) + + # Test maximum coordinates (all quantics at max value for each site's local dimension) + @test fq_complex([3, 9, 9, 27]) ≈ f_complex(QuanticsGrids.grid_max(grid_complex)...) + + # Test intermediate values with specific quantics combinations + # Choose quantics that should map to approximately middle values + quantics_mid = [2, 5, 5, 14] # Should map to roughly middle of each dimension + coord_mid = quantics_to_origcoord(grid_complex, quantics_mid) + @test fq_complex(quantics_mid) ≈ f_complex(coord_mid...) +end + +@testitem "unfoldingscheme" begin + R = 4 + + # DiscretizedGrid should support unfoldingscheme parameter and property + grid_fused = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), unfoldingscheme=:fused) + grid_interleaved = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), unfoldingscheme=:interleaved) + + # Test 2D grids with different unfolding schemes + grid_2d_fused = DiscretizedGrid((R, R); lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0), unfoldingscheme=:fused) + grid_2d_interleaved = DiscretizedGrid((R, R); lower_bound=(0.0, 0.0), upper_bound=(1.0, 1.0), unfoldingscheme=:interleaved) + + # Test that unfoldingscheme affects quantics behavior + # For 2D grids: fused should have length R, interleaved should have length 2*R + grididx = (2, 2) + quantics_fused = grididx_to_quantics(grid_2d_fused, grididx) + quantics_interleaved = grididx_to_quantics(grid_2d_interleaved, grididx) + + @test length(quantics_fused) == R # fused: length R + @test length(quantics_interleaved) == 2 * R # interleaved: length 2*R for 2D + + # Test that both produce same coordinate when converted back + coord_fused = quantics_to_origcoord(grid_2d_fused, quantics_fused) + coord_interleaved = quantics_to_origcoord(grid_2d_interleaved, quantics_interleaved) + @test coord_fused == coord_interleaved + + # Test default unfoldingscheme (should be :fused) + grid_default = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,)) + + # Test localdimensions with different unfolding schemes + @test QuanticsGrids.localdimensions(grid_2d_fused) == fill(4, R) # 2^2 = 4 for fused + @test QuanticsGrids.localdimensions(grid_2d_interleaved) == fill(2, 2 * R) # 2 for interleaved +end + +@testitem "includeendpoint" begin + R = 4 + + # DiscretizedGrid should support includeendpoint parameter and property + grid_without = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), includeendpoint=false) + grid_with = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,), includeendpoint=true) + + # Test that includeendpoint affects grid behavior + @test QuanticsGrids.grid_max(grid_without) ≈ 1.0 - (1.0 / 2^R) + @test QuanticsGrids.grid_max(grid_with) == 1.0 + + # Test grid_step calculation differences + step_without = QuanticsGrids.grid_step(grid_without) + step_with = QuanticsGrids.grid_step(grid_with) + + @test step_without ≈ 1.0 / 2^R + @test step_with ≈ 1.0 / (2^R - 1) + + # Test that includeendpoint=true allows accessing the upper bound + @test grididx_to_origcoord(grid_with, 2^R) == 1.0 + @test grididx_to_origcoord(grid_without, 2^R) ≈ 1.0 - (1.0 / 2^R) + + # Test 2D case + grid_2d_without = DiscretizedGrid((R, R); + lower_bound=(0.0, 0.0), + upper_bound=(2.0, 3.0), + includeendpoint=false) + grid_2d_with = DiscretizedGrid((R, R); + lower_bound=(0.0, 0.0), + upper_bound=(2.0, 3.0), + includeendpoint=true) + + # Test grid_max behavior in 2D + max_without = QuanticsGrids.grid_max(grid_2d_without) + max_with = QuanticsGrids.grid_max(grid_2d_with) + + @test max_without == (2.0 - 2.0 / 2^R, 3.0 - 3.0 / 2^R) + @test max_with == (2.0, 3.0) + + # Test that the upper bound point is accessible when includeendpoint=true + @test grididx_to_origcoord(grid_2d_with, (2^R, 2^R)) == (2.0, 3.0) + + # CHALLENGING TEST: Complex indextable with mixed site sizes and includeendpoint behavior + # This tests how includeendpoint interacts with custom indextables having different local dimensions + complex_variablenames = (:x, :y, :z) + complex_indextable = [ + [(:x, 1)], + [(:y, 1), (:z, 1)], + [(:x, 2), (:y, 2), (:z, 2)], + [(:x, 3), (:y, 3)] + ] + + # Test with base=3 and asymmetric domain bounds + lower_bound = (-2.5, 1.2, 0.0) + upper_bound = (4.7, 8.1, 3.3) + + grid_complex_without = DiscretizedGrid(complex_variablenames, complex_indextable; + lower_bound, upper_bound, base=3, includeendpoint=false) + grid_complex_with = DiscretizedGrid(complex_variablenames, complex_indextable; + lower_bound, upper_bound, base=3, includeendpoint=true) + + # Test that local dimensions are same regardless of includeendpoint + local_dims_without = QuanticsGrids.localdimensions(grid_complex_without) + local_dims_with = QuanticsGrids.localdimensions(grid_complex_with) + expected_local_dims = [3, 9, 27, 9] # [3^1, 3^2, 3^3, 3^2] + @test local_dims_without == expected_local_dims + @test local_dims_with == expected_local_dims + + # Test grid_max behavior with complex grids + max_without = QuanticsGrids.grid_max(grid_complex_without) + max_with = QuanticsGrids.grid_max(grid_complex_with) + + # Without includeendpoint: should be slightly less than upper_bound bound + expected_max_without = ( + upper_bound[1] - (upper_bound[1] - lower_bound[1]) / 3^3, # x dimension uses 3 quantics + upper_bound[2] - (upper_bound[2] - lower_bound[2]) / 3^3, # y dimension uses 3 quantics + upper_bound[3] - (upper_bound[3] - lower_bound[3]) / 3^2 # z dimension uses 2 quantics + ) + @test all(max_without .≈ expected_max_without) + @test max_with == upper_bound + + # Test maximum quantics combinations for each site's local dimension + max_quantics_without = [3, 9, 27, 9] # maximum for each site + max_quantics_with = [3, 9, 27, 9] # same maximum values + + coord_max_without = quantics_to_origcoord(grid_complex_without, max_quantics_without) + coord_max_with = quantics_to_origcoord(grid_complex_with, max_quantics_with) + + # With includeendpoint=true, max quantics should reach exact upper_bound bound + @test all(coord_max_with .≈ upper_bound) + + # With includeendpoint=false, max quantics should be slightly less than upper_bound bound + @test all(coord_max_without .< upper_bound) + @test all(coord_max_without .≈ expected_max_without) + + # Test quanticsfunction behavior with boundary values + # Define a function that's sensitive to boundary effects + boundary_sensitive_func(x, y, z) = (x - upper_bound[1])^4 + (y - upper_bound[2])^4 + (z - upper_bound[3])^4 + + fq_without = QuanticsGrids.quanticsfunction(Float64, grid_complex_without, boundary_sensitive_func) + fq_with = QuanticsGrids.quanticsfunction(Float64, grid_complex_with, boundary_sensitive_func) + + # Test that includeendpoint=true allows exact evaluation at upper_bound bound + @test fq_with(max_quantics_with) ≈ boundary_sensitive_func(upper_bound...) + + # Test that includeendpoint=false gives a non-zero value (since it can't reach exact upper_bound bound) + result_without = fq_without(max_quantics_without) + @test result_without > 0 # Should be positive since we can't reach the upper_bound bound exactly + @test result_without ≈ boundary_sensitive_func(coord_max_without...) + + # Test step size consistency across different sites + # For complex indextables, step sizes should vary by dimension based on quantics count + step_without = QuanticsGrids.grid_step(grid_complex_without) + step_with = QuanticsGrids.grid_step(grid_complex_with) + + # Expected step sizes for each dimension + expected_step_without = ( + (upper_bound[1] - lower_bound[1]) / 3^3, # x: uses 3 quantics per site + (upper_bound[2] - lower_bound[2]) / 3^3, # y: uses 3 quantics per site + (upper_bound[3] - lower_bound[3]) / 3^2 # z: uses 2 quantics per site + ) + expected_step_with = ( + (upper_bound[1] - lower_bound[1]) / (3^3 - 1), # x: adjusted for includeendpoint + (upper_bound[2] - lower_bound[2]) / (3^3 - 1), # y: adjusted for includeendpoint + (upper_bound[3] - lower_bound[3]) / (3^2 - 1) # z: adjusted for includeendpoint + ) + + @test all(step_without .≈ expected_step_without) + @test all(step_with .≈ expected_step_with) + + # Test that grid_origin is unaffected by includeendpoint (should always be lower_bound) + @test QuanticsGrids.grid_origin(grid_complex_without) == lower_bound + @test QuanticsGrids.grid_origin(grid_complex_with) == lower_bound + + # Test quantics-to-coordinate conversion consistency at boundaries + # Minimum quantics (all 1s) should always map to lower_bound regardless of includeendpoint + min_quantics = [1, 1, 1, 1] + @test all(quantics_to_origcoord(grid_complex_without, min_quantics) .≈ lower_bound) + @test all(quantics_to_origcoord(grid_complex_with, min_quantics) .≈ lower_bound) + + # Test intermediate quantics values behave consistently + mid_quantics = [2, 5, 14, 5] # roughly middle values for each site + coord_mid_without = quantics_to_origcoord(grid_complex_without, mid_quantics) + coord_mid_with = quantics_to_origcoord(grid_complex_with, mid_quantics) + + # Coordinates should be different due to different step sizes + @test coord_mid_without != coord_mid_with + + # But both should be within the domain bounds + @test all(lower_bound .<= coord_mid_without .<= upper_bound) + @test all(lower_bound .<= coord_mid_with .<= upper_bound) +end + +@testitem "boundary error handling" begin + R = 5 + g = DiscretizedGrid{1}(R) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, -0.1) + @test_throws DomainError QuanticsGrids.origcoord_to_grididx(g, 1.1) +end + +@testitem "large grids" begin + R = 62 + g = DiscretizedGrid{1}(R) + @test grididx_to_quantics(g, 2^R) == fill(2, R) +end + +@testitem "grid representation conversion _NEW_" begin + R = 10 + reprs = [:grididx, :quantics, :origcoord] + + varnames = (:x, :y) + tupletable = [ + [(:x, 3), (:y, 2)], + [(:y, 1), (:x, 1)], + [(:x, 2)], + ] + + testset = [ + (DiscretizedGrid{1}(R, -3.2, 4.8), 2), + (DiscretizedGrid{2}(R, (-2.3, 1.2), (9.0, 3.5)), (2, 3)), + (DiscretizedGrid(varnames, tupletable), (7, 2)), + ] + + for (grid, ini_grididx) in testset + + data = Dict{Symbol,Any}() + transforms = Dict( + (:grididx, :quantics) => grididx_to_quantics, + (:grididx, :origcoord) => grididx_to_origcoord, + (:quantics, :grididx) => quantics_to_grididx, + (:quantics, :origcoord) => quantics_to_origcoord, + (:origcoord, :grididx) => QuanticsGrids.origcoord_to_grididx, + (:origcoord, :quantics) => origcoord_to_quantics, + ) + + data[:grididx] = ini_grididx + while true + newdata = false + for src in reprs, dst in reprs + if src == dst + continue + end + if !(src ∈ keys(data)) + continue + end + + if dst ∈ keys(data) + @test transforms[(src, dst)](grid, data[src]) == data[dst] + else + data[dst] = @inferred transforms[(src, dst)](grid, data[src]) + newdata = true + end + end + if newdata == false + break + end + end + end +end + +@testitem "step size and origin" begin + R = 4 + g = DiscretizedGrid{1}(R) + @test QuanticsGrids.grid_step(g) == 1.0 / 2^R + @test QuanticsGrids.grid_origin(g) == 0.0 +end + +@testitem "grid_origcoords" begin + # Test 1D grid with default parameters + R = 4 + grid_1d = DiscretizedGrid{1}(R; lower_bound=(0.0,), upper_bound=(1.0,)) + coords_1d = QuanticsGrids.grid_origcoords(grid_1d, 1) + + # Should return a range with correct length and bounds + @test length(coords_1d) == 2^R + @test coords_1d[1] ≈ QuanticsGrids.grid_min(grid_1d) + @test coords_1d[end] ≈ QuanticsGrids.grid_max(grid_1d) + + # Test that coordinates match grididx_to_origcoord for all indices + for i in 1:2^R + @test coords_1d[i] ≈ grididx_to_origcoord(grid_1d, i) + end + + # Test 2D grid + grid_2d = DiscretizedGrid((3, 4); lower_bound=(-1.0, 2.0), upper_bound=(5.0, 8.0)) + + # Test first dimension + coords_x = QuanticsGrids.grid_origcoords(grid_2d, 1) + @test length(coords_x) == 2^3 + @test coords_x[1] ≈ QuanticsGrids.grid_min(grid_2d)[1] + @test coords_x[end] ≈ QuanticsGrids.grid_max(grid_2d)[1] + + # Test second dimension + coords_y = QuanticsGrids.grid_origcoords(grid_2d, 2) + @test length(coords_y) == 2^4 + @test coords_y[1] ≈ QuanticsGrids.grid_min(grid_2d)[2] + @test coords_y[end] ≈ QuanticsGrids.grid_max(grid_2d)[2] + + # Verify coordinates match grididx_to_origcoord + for i in 1:2^3 + @test coords_x[i] ≈ grididx_to_origcoord(grid_2d, (i, 1))[1] + end + for j in 1:2^4 + @test coords_y[j] ≈ grididx_to_origcoord(grid_2d, (1, j))[2] + end + + # Test with different base + grid_base3 = DiscretizedGrid{1}(3; lower_bound=(0.0,), upper_bound=(2.7,), base=3) + coords_base3 = QuanticsGrids.grid_origcoords(grid_base3, 1) + @test length(coords_base3) == 3^3 + @test coords_base3[1] ≈ QuanticsGrids.grid_min(grid_base3) + @test coords_base3[end] ≈ QuanticsGrids.grid_max(grid_base3) + + # Test with includeendpoint=true + grid_endpoint = DiscretizedGrid{1}(3; lower_bound=(0.0,), upper_bound=(1.0,), includeendpoint=true) + coords_endpoint = QuanticsGrids.grid_origcoords(grid_endpoint, 1) + @test length(coords_endpoint) == 2^3 + @test coords_endpoint[1] ≈ 0.0 + @test coords_endpoint[end] ≈ 1.0 # Should reach exact upper bound + + # Test with custom indextable + variablenames = (:x, :y, :z) + indextable = [ + [(:x, 1), (:y, 1)], # site 1: 2 indices + [(:z, 1)], # site 2: 1 index + [(:x, 2), (:y, 2), (:z, 2)], # site 3: 3 indices + [(:x, 3), (:y, 3)] # site 4: 2 indices + ] + grid_custom = DiscretizedGrid(variablenames, indextable; + lower_bound=(-2.0, 1.0, 0.5), upper_bound=(3.0, 4.0, 2.0), base=2) + + # Test each dimension + coords_x_custom = QuanticsGrids.grid_origcoords(grid_custom, :x) # x dimension (3 quantics) + coords_y_custom = QuanticsGrids.grid_origcoords(grid_custom, :y) # y dimension (3 quantics) + coords_z_custom = QuanticsGrids.grid_origcoords(grid_custom, :z) # z dimension (2 quantics) + + @test length(coords_x_custom) == 2^3 # x has 3 quantics + @test length(coords_y_custom) == 2^3 # y has 3 quantics + @test length(coords_z_custom) == 2^2 # z has 2 quantics + + # Test bounds + @test coords_x_custom[1] ≈ QuanticsGrids.grid_min(grid_custom)[1] + @test coords_y_custom[1] ≈ QuanticsGrids.grid_min(grid_custom)[2] + @test coords_z_custom[1] ≈ QuanticsGrids.grid_min(grid_custom)[3] + @test coords_x_custom[end] ≈ QuanticsGrids.grid_max(grid_custom)[1] + @test coords_y_custom[end] ≈ QuanticsGrids.grid_max(grid_custom)[2] + @test coords_z_custom[end] ≈ QuanticsGrids.grid_max(grid_custom)[3] + + # Test error handling - dimension out of bounds + @test_throws DomainError QuanticsGrids.grid_origcoords(grid_1d, 0) + @test_throws DomainError QuanticsGrids.grid_origcoords(grid_1d, 2) + @test_throws DomainError QuanticsGrids.grid_origcoords(grid_2d, 3) + @test_throws DomainError QuanticsGrids.grid_origcoords(grid_custom, 4) + + # Test that returned range is iterable and has correct type + coords = QuanticsGrids.grid_origcoords(grid_1d, 1) + @test coords isa AbstractRange + @test eltype(coords) <: AbstractFloat + + # Test consistency with grid spacing + grid_spacing = QuanticsGrids.grid_step(grid_2d) + coords_x_spacing = coords_x[2] - coords_x[1] + coords_y_spacing = coords_y[2] - coords_y[1] + @test coords_x_spacing ≈ grid_spacing[1] + @test coords_y_spacing ≈ grid_spacing[2] + + # CHALLENGING TEST: Complex multi-dimensional grid with asymmetric resolutions + complex_Rs = (5, 3, 7, 2) # Different resolutions per dimension + grid_complex = DiscretizedGrid(complex_Rs; + lower_bound=(-10.0, 0.1, 50.0, -2.5), + upper_bound=(15.0, 3.9, 100.0, 7.8), + base=3) + + # Test each dimension has correct properties + for d in 1:4 + coords_d = QuanticsGrids.grid_origcoords(grid_complex, d) + expected_length = 3^complex_Rs[d] + @test length(coords_d) == expected_length + @test coords_d[1] ≈ QuanticsGrids.grid_min(grid_complex)[d] + @test coords_d[end] ≈ QuanticsGrids.grid_max(grid_complex)[d] + + # Test uniform spacing within each dimension + spacings = diff(collect(coords_d)) + @test all(spacing -> spacing ≈ spacings[1], spacings) + @test spacings[1] ≈ QuanticsGrids.grid_step(grid_complex)[d] + + # Test consistency with grididx_to_origcoord + test_indices = [1, expected_length ÷ 2, expected_length] + for idx in test_indices + grid_idx = ntuple(i -> i == d ? idx : 1, 4) + expected_coord = grididx_to_origcoord(grid_complex, grid_idx)[d] + @test coords_d[idx] ≈ expected_coord + end + end + + # Test with very small grid (edge case) + grid_tiny = DiscretizedGrid{1}(1; lower_bound=(0.0,), upper_bound=(1.0,)) + coords_tiny = QuanticsGrids.grid_origcoords(grid_tiny, 1) + @test length(coords_tiny) == 2 + @test coords_tiny[1] ≈ 0.0 + @test coords_tiny[2] ≈ 0.5 # For R=1, max coord should be 1 - 1/2 = 0.5 +end + +@testitem "DiscretizedGrid(R::Int, lower_bound, upper_bound) constructor" begin + grid = DiscretizedGrid(3, -2.0, 3.0) + @test QuanticsGrids.grid_Rs(grid) == (3,) + @test QuanticsGrids.lower_bound(grid) == -2.0 + @test QuanticsGrids.upper_bound(grid) == 3.0 +end