diff --git a/.gitignore b/.gitignore index 6786d43..86578cd 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,5 @@ Manifest*.toml ref/ .claude/ data/.* -data/* \ No newline at end of file +data/* +!/data/a_sparse_cdf.cdf diff --git a/data/a_sparse_cdf.cdf b/data/a_sparse_cdf.cdf new file mode 100644 index 0000000..89814fb Binary files /dev/null and b/data/a_sparse_cdf.cdf differ diff --git a/src/loading/variable.jl b/src/loading/variable.jl index 020042b..93b07e6 100644 --- a/src/loading/variable.jl +++ b/src/loading/variable.jl @@ -116,12 +116,11 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r buffer = parent(var.parentdataset) RecordSizeType = recordsize_type(var.parentdataset) entries, vvr_type = read_vvrs(var.vdr) + # vvr record type is the ultimate source of compression + compression = vvr_type == VVR_ ? NoCompression : variable_compression(var.vdr) + sparse_type(var.vdr) == 0 || + return _readblock_sparse!(var, dest, ranges, entries, compression, buffer, RecordSizeType) isempty(entries) && return dest - compression = if !isempty(entries) # # vvr records is the ultimative source - vvr_type == VVR_ ? NoCompression : variable_compression(var.vdr) - else - NoCompression - end record_range = ranges[end] other_ranges = ranges[1:(N - 1)] @@ -195,6 +194,65 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r return dest end +# Sparse records: records absent from the VXR are virtual. Pad sparse (1) fills them +# with the VDR pad value (or the NASA default pad); previous sparse (2) repeats the +# last record of the preceding physical block. +# Mirrors cdflib semantics; spec'd in the CDF IFD (sRecords). +function _readblock_sparse!(var::CDFVariable{T, N}, dest, ranges, entries, compression, buffer, ::Type{RST}) where {T, N, RST} + record_range = ranges[end] + other_ranges = ranges[1:(N - 1)] + dims_without_record = var.dims[1:(N - 1)] + record_size = prod(dims_without_record) + is_row_major = majority(var) == Row + needs_byte_swap = is_big_endian_encoding(var) + use_prev = sparse_type(var.vdr) == 2 + pad = pad_value(var.vdr, T, needs_byte_swap) + + chunk = Vector{T}() + cached = 0 + # load (and permute, for row-major) a physical block once; runs of records from + # the same block (incl. prev-sparse repeats) reuse it + load_block = function (idx) + entry = entries[idx] + if cached != idx + resize!(chunk, record_size * length(entry)) + if compression == NoCompression + load_vvr_data!(chunk, 1, buffer, entry.offset, length(chunk), RST) + else + decompressor = take!(decompressors()) + try + load_cvvr_data!(chunk, 1, buffer, entry.offset, length(chunk), RST, compression; decompressor) + finally + put!(decompressors(), decompressor) + end + end + is_row_major && majority_swap!(reshape(chunk, dims_without_record..., :), dims_without_record) + cached = idx + end + return reshape(chunk, dims_without_record..., :) + end + + blk = 1 + nblocks = length(entries) + for (di, r) in enumerate(record_range) + while blk <= nblocks && entries[blk].last < r + blk += 1 + end + dest_view = _record_view(dest, di) + if blk <= nblocks && entries[blk].first <= r + arr = load_block(blk) + dest_view .= view(arr, other_ranges..., r - entries[blk].first + 1) + elseif use_prev && blk > 1 + arr = load_block(blk - 1) + dest_view .= view(arr, other_ranges..., size(arr, N)) + else + fill!(dest_view, pad) + end + end + needs_byte_swap && _byte_swap!(dest) + return dest +end + function collect_vxr_entries!(entries::Vector{VVREntry}, src, offset, ::Type{FieldSizeT}) where {FieldSizeT} vvr_type = nothing while offset != 0 diff --git a/src/records/vdr.jl b/src/records/vdr.jl index 47b1eae..8a88562 100644 --- a/src/records/vdr.jl +++ b/src/records/vdr.jl @@ -154,4 +154,40 @@ end is_record_varying(vdr) = !is_nrv(vdr) """Whether or not the variable is a non-record variable""" is_nrv(vdr) = (vdr.flags & 0x01) == 0 +has_pad_value(vdr::AbstractVDR) = (vdr.flags & 0x02) != 0 is_compressed(vdr::AbstractVDR) = (vdr.flags & 0x04) != 0 + +# sRecords: 0 = none, 1 = pad sparse, 2 = previous-record sparse +sparse_type(vdr::AbstractVDR) = Int(vdr.s_records) + +# PadValue trails zDimSizes+DimVarys (zVDR) / DimVarys (rVDR), both starting at vdr.pos +_pad_value_pos(vdr::VDR) = vdr.pos + 8 * Int(vdr.num_dims) +_pad_value_pos(vdr::rVDR) = vdr.pos + 4 * Int(vdr.gdr.r_num_dims) + +# Returned in file encoding: virtual records are filled before the final byte swap +# in `readblock!`, which fixes pad and physical data together. +function pad_value(vdr::AbstractVDR, ::Type{T}, needs_byte_swap) where {T} + if has_pad_value(vdr) + buf = vdr.buffer + return GC.@preserve buf unsafe_load(convert(Ptr{T}, pointer(buf, _pad_value_pos(vdr)))) + end + return _file_encode(default_pad(T), needs_byte_swap) +end + +_file_encode(x, needs_byte_swap) = needs_byte_swap ? hton(x) : x +_file_encode(x::StaticString, needs_byte_swap) = x # excluded from byte swap + +# NASA CDF library default pad values (≠ ISTP FILLVAL) +default_pad(::Type{Int8}) = Int8(-127) +default_pad(::Type{Int16}) = Int16(-32767) +default_pad(::Type{Int32}) = Int32(-2147483647) +default_pad(::Type{Int64}) = Int64(-9223372036854775807) +default_pad(::Type{UInt8}) = UInt8(254) +default_pad(::Type{UInt16}) = UInt16(65534) +default_pad(::Type{UInt32}) = UInt32(4294967294) +default_pad(::Type{Float32}) = Float32(-1.0e30) +default_pad(::Type{Float64}) = -1.0e30 +default_pad(::Type{Epoch}) = Epoch(-1.0e30) +default_pad(::Type{Epoch16}) = Epoch16(-1.0e30, -1.0e30) +default_pad(::Type{TT2000}) = TT2000(Int64(-9223372036854775807)) +default_pad(::Type{StaticString{N, UInt8}}) where {N} = StaticString{N, UInt8}(ntuple(_ -> UInt8(' '), Val(N))) diff --git a/test/make_sparse_cdf.py b/test/make_sparse_cdf.py new file mode 100644 index 0000000..9cd85ea --- /dev/null +++ b/test/make_sparse_cdf.py @@ -0,0 +1,45 @@ +# Generate data/a_sparse_cdf.cdf: sparse-record fixtures (pad & prev semantics). +# Run: uv run --with cdflib python test/make_sparse_cdf.py +from pathlib import Path + +import numpy as np +from cdflib import cdfwrite + +out = Path(__file__).parent.parent / "data" / "a_sparse_cdf.cdf" +out.unlink(missing_ok=True) +f = cdfwrite.CDF(out) + +# physical records 0-2, 6-7, 10 (0-based); virtual gaps at 3-5, 8-9 +phys = np.array([0, 1, 2, 6, 7, 10]) +data = np.array([10.0, 11.0, 12.0, 16.0, 17.0, 20.0]) + +base = { + "Data_Type": 45, # CDF_DOUBLE + "Num_Elements": 1, + "Rec_Vary": True, + "Dim_Sizes": [], + "Compress": 0, + "Block_Factor": 3, # force multiple VVR blocks +} +f.write_var({**base, "Variable": "pad_default", "Sparse": "pad_sparse"}, var_data=[phys, data]) +f.write_var( + {**base, "Variable": "pad_explicit", "Sparse": "pad_sparse", "Pad": np.array([-99.0])}, + var_data=[phys, data], +) +f.write_var({**base, "Variable": "prev", "Sparse": "prev_sparse"}, var_data=[phys, data]) + +data2d = np.arange(12.0).reshape(6, 2) +f.write_var( + {**base, "Variable": "pad2d", "Dim_Sizes": [2], "Sparse": "pad_sparse", "Pad": np.array([-99.0])}, + var_data=[phys, data2d], +) +f.close() + +# sanity-check with the cdflib reader +import cdflib + +r = cdflib.CDF(out) +print("pad_default:", r.varget("pad_default")) +print("pad_explicit:", r.varget("pad_explicit")) +print("prev:", r.varget("prev")) +print("pad2d:", r.varget("pad2d").tolist()) diff --git a/test/runtests.jl b/test/runtests.jl index 097ba19..efa6cad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ include("comprehensive_test.jl") include("cdf2_test.jl") include("CommonDataModelExt_test.jl") include("decompress_test.jl") +include("sparse_test.jl") @testset "StaticString" begin include("staticstring.jl") end diff --git a/test/sparse_test.jl b/test/sparse_test.jl new file mode 100644 index 0000000..fcf86a8 --- /dev/null +++ b/test/sparse_test.jl @@ -0,0 +1,32 @@ +# Fixture from test/make_sparse_cdf.py (cdflib writer): physical records 1-3, 7-8, 11 +# (1-based), virtual gaps at 4-6 and 9-10. +@testset "Sparse records" begin + ds = CDFDataset(data_path("a_sparse_cdf.cdf")) + + phys = [10.0, 11.0, 12.0, 16.0, 17.0, 20.0] + + # pad sparse, NASA default pad for CDF_DOUBLE + d = -1.0e30 + @test ds["pad_default"][:] == [10.0, 11.0, 12.0, d, d, d, 16.0, 17.0, d, d, 20.0] + + # pad sparse with explicit VDR pad value + p = ds["pad_explicit"] + @test p[:] == [10.0, 11.0, 12.0, -99.0, -99.0, -99.0, 16.0, 17.0, -99.0, -99.0, 20.0] + @test p[3:5] == [12.0, -99.0, -99.0] # partial range straddling physical/virtual + @test p[4:4] == [-99.0] # purely virtual range + @test p[11:11] == [20.0] + + # previous sparse: repeat last record of preceding block; pad before first block + @test ds["prev"][:] == [10.0, 11.0, 12.0, 12.0, 12.0, 12.0, 16.0, 17.0, 17.0, 17.0, 20.0] + + # 2D: whole virtual record takes the pad value (cdflib's reader only pads the + # first element — spec says the full record) + m = ds["pad2d"][:, :] + @test size(m) == (2, 11) + @test m[:, 1] == [0.0, 1.0] + @test m[:, 4] == [-99.0, -99.0] + @test m[:, 11] == [10.0, 11.0] + @test ds["pad2d"][1:1, 5:7] == [-99.0 -99.0 6.0] + + @test read(ds, "prev", Vector{Float64}) == ds["prev"][:] +end