diff --git a/claudedocs/majority_handling_analysis.md b/claudedocs/majority_handling_analysis.md new file mode 100644 index 0000000..8ea2a9a --- /dev/null +++ b/claudedocs/majority_handling_analysis.md @@ -0,0 +1,119 @@ +# CDF Majority Handling Analysis + +## Summary + +This document explains how CDF files handle row-major vs column-major data layout and how different implementations (cdflib, CDFpp) handle this difference. + +## Background: Row-Major vs Column-Major + +### Data Layout +- **Row-Major** (C-style): Elements of each row are stored contiguously in memory + - Used by: C, C++, Python (NumPy default), Row + - For array `A[i,j,k]`, index `k` varies fastest + +- **Column-Major** (Fortran-style): Elements of each column are stored contiguously + - Used by: Fortran, MATLAB, Julia, R + - For array `A[i,j,k]`, index `i` varies fastest + +### CDF Specification +CDF files can be created with either: +- **Row-Major** (flag bit 0 = 1): Row-majority storage +- **Column-Major** (flag bit 0 = 0): Column-majority storage + +The majority is stored in the CDR (CDF Descriptor Record) flags field: +- Bit 0: Majority (1 = row-major, 0 = column-major) + +## Python cdflib Implementation + +### Dimension Reversal for Column-Major +File: `ref/cdflib/cdflib/cdfread.py:1655-1656` +```python +if self._majority == "Column_major": + dimensions = list(reversed(dimensions)) +``` + +### Array Transposition After Loading +File: `ref/cdflib/cdflib/cdfread.py:1749-1754` +```python +if self._majority == "Column_major": + if dimensions is not None: + axes = [0] + list(range(len(dimensions), 0, -1)) + else: + axes = None + ret = np.transpose(ret, axes=axes) +``` + +### Logic +1. **Before loading**: Reverse dimension sizes if column-major +2. **After loading**: Transpose the array to match Python's row-major layout +3. **Record dimension**: First dimension (axis 0) is always the record count and is NOT reversed + +Example for column-major CDF with dims `[3, 4, 5]` and 10 records: +1. Read dims as `[3, 4, 5]` +2. Reverse to `[5, 4, 3]` for loading +3. Shape becomes `(10, 5, 4, 3)` after adding records +4. Transpose with `axes=[0, 3, 2, 1]` → final shape `(10, 3, 4, 5)` + +## C++ CDFpp Implementation + +### Majority Swap Function +File: [ref/CDFpp/include/cdfpp/cdf-io/majority-swap.hpp](../ref/CDFpp/include/cdfpp/cdf-io/majority-swap.hpp) + +The implementation uses an in-place element reordering strategy: + +```cpp +template +void swap(data_t& data, const shape_t& shape) +{ + // Only swap if dimensions > 2 (excluding record dimension) + if ((dimensions > 2 && !is_string) or (is_string and dimensions > 3)) + { + // Generate access pattern for index mapping + const auto access_patern = _private::generate_access_pattern(record_shape); + + // For each record, reorder elements according to access pattern + for (auto record = 0UL; record < records_count; record++) + { + for (const auto& swap_pair : access_patern) + { + temporary_record[swap_pair.src] = data[offset + swap_pair.dest]; + } + std::memcpy(data.data() + offset, temporary_record.data(), bytes_per_record); + offset += elements_per_record; + } + } +} +``` + +### When Swap is Applied +File: [ref/CDFpp/include/cdfpp/variable.hpp:93-96](../ref/CDFpp/include/cdfpp/variable.hpp#L93-L96) + +```cpp +Variable(const std::string& name, std::size_t number, data_t&& data, shape_t&& shape, + cdf_majority majority = cdf_majority::row, ...) +{ + if (this->majority() == cdf_majority::column) + { + majority::swap(_data(), p_shape); // Swap for column-major files + } +} +``` + +### Logic +1. **Load data as-is** from the file +2. **If column-major**: Apply in-place element reordering +3. **Result**: Data layout matches the C++ row-major convention + +### Index Calculation +The key is the `flat_index` calculation: + +- **Row-major index**: `i*shape[1]*shape[2] + j*shape[2] + k` +- **Column-major index**: `k*shape[1]*shape[0] + j*shape[0] + i` + +The swap function creates a mapping between these two indexing schemes. + +## References + +- CDF Internal Format Specification: https://cdf.gsfc.nasa.gov/ +- Python cdflib: [ref/cdflib/cdflib/cdfread.py](../ref/cdflib/cdflib/cdfread.py) +- CDFpp: [ref/CDFpp/include/cdfpp/cdf-io/majority-swap.hpp](../ref/CDFpp/include/cdfpp/cdf-io/majority-swap.hpp) diff --git a/src/dataset.jl b/src/dataset.jl index ad85ae7..e4aa08b 100644 --- a/src/dataset.jl +++ b/src/dataset.jl @@ -39,14 +39,15 @@ function CDFDataset(filename) end is_compressed(magic_numbers::UInt32) = magic_numbers != 0x0000FFFF +majority(cdf::CDFDataset) = majority(cdf.cdr) # Convenience accessors for the dataset with lazy loading -function Base.getproperty(cdf::CDFDataset{CT}, name::Symbol) where {CT} +@inline function Base.getproperty(cdf::CDFDataset{CT, FST}, name::Symbol) where {CT, FST} name in fieldnames(CDFDataset) && return getfield(cdf, name) if name === :version return version(cdf.cdr) elseif name === :majority - return Majority(cdf.cdr) + return majority(cdf) elseif name === :compression return CT elseif name === :adr diff --git a/src/enums.jl b/src/enums.jl index 6281e43..70ed7cf 100644 --- a/src/enums.jl +++ b/src/enums.jl @@ -1,6 +1,6 @@ -@enum Majority begin - Row = 0 - Column = 1 +@enum Majority::Bool begin + Row = false + Column = true end @enum RecordType::Int8 begin diff --git a/src/loading/variable.jl b/src/loading/variable.jl index daa7e94..2d0da9a 100644 --- a/src/loading/variable.jl +++ b/src/loading/variable.jl @@ -18,6 +18,24 @@ else end end +""" + majority_swap!(data, dims_without_record) + +Convert row-major data layout to column-major (Julia's native layout) in-place. +For row-major CDF files, this reverses the dimension order (except record dimension). +""" +function majority_swap!(data::AbstractArray{T, N}, dims_without_record) where {T, N} + N <= 2 && return data + perm = ntuple(i -> N - i, N - 1) + reversed_dims = reverse(dims_without_record) + temp = similar(data, dims_without_record) + for slc in eachslice(data; dims = N) + permutedims!(temp, reshape(slc, reversed_dims), perm) + copyto!(slc, temp) + end + return data +end + function variable(cdf::CDFDataset, name) vdr = find_vdr(cdf, name) @@ -59,6 +77,7 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r end_idx = findfirst(entry -> entry.first <= last_rec <= entry.last, entries) record_size = prod(dims_without_record) + is_row_major = majority(var) == Row # If the variable is not compressed and the other dimension ranges are the same as the variable range # we can directly read the data into dest if is_no_compression && is_full_record @@ -76,8 +95,9 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r doffs += N_elems end @assert doffs == length(dest) + 1 + is_row_major && majority_swap!(dest, dims_without_record) else - Base.@inbounds Threads.@threads for i in eachindex(start_idx:end_idx) + Base.@inbounds Threads.@threads for i in start_idx:end_idx entry = entries[i] if is_full_record && entry.first >= first_rec && entry.last <= last_rec # full entry @@ -86,6 +106,7 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r total_elems = record_size * length(entry) decompressor = take!(decompressors()) load_cvvr_data!(dest_view, 1, buffer, entry.offset, total_elems, RecordSizeType, compression; decompressor) + is_row_major && majority_swap!(dest_view, dims_without_record) put!(decompressors(), decompressor) else # partial entry @@ -105,6 +126,7 @@ function DiskArrays.readblock!(var::CDFVariable{T, N}, dest::AbstractArray{T}, r # chunk_data = _load_entry_chunk(var, entry, RecordSizeType, buffer, var.compression; decompressor) chunk_array = reshape(chunk, dims_without_record..., :) + is_row_major && majority_swap!(chunk_array, dims_without_record) src_view = view(chunk_array, other_ranges..., local_range) dest_view .= src_view end diff --git a/src/records/cdr.jl b/src/records/cdr.jl index 6d4c854..01db15a 100644 --- a/src/records/cdr.jl +++ b/src/records/cdr.jl @@ -18,7 +18,7 @@ struct CDR{FST} <: Record end version(cdr::CDR; verbose = true) = verbose ? (cdr.version, cdr.release, cdr.increment) : cdr.version -Majority(cdr::CDR) = (cdr.flags & 0x01) != 0 ? Majority(0) : Majority(1) # Row=0, Column=1 +majority(cdr::CDR) = (cdr.flags & 0x01) != 0 ? Row : Column # Row=0, Column=1 is_cdf_v3(cdr::CDR) = cdr.version == 3 """ diff --git a/src/records/records.jl b/src/records/records.jl index 031641b..5ab6b22 100644 --- a/src/records/records.jl +++ b/src/records/records.jl @@ -74,7 +74,6 @@ Decode the CDR flags field into individual boolean flags. function decode_cdr_flags(flags) flags = UInt32(flags) return ( - row_majority = (flags & 0x01) != 0, single_file_format = (flags & 0x02) != 0, checksum_used = (flags & 0x04) != 0, md5_checksum = (flags & 0x08) != 0, @@ -92,7 +91,7 @@ function Base.show(io::IO, cdr::CDR) println(io, " Version: $(cdr.version).$(cdr.release).$(cdr.increment)") println(io, " Encoding: $(cdr.encoding)") println(io, " Flags: 0x$(string(cdr.flags, base = 16, pad = 8))") - println(io, " - Row Majority: $(flag_info.row_majority)") + println(io, " - Majority: $(majority(cdr))") println(io, " - Single File Format: $(flag_info.single_file_format)") println(io, " - Checksum Used: $(flag_info.checksum_used)") println(io, " - MD5 Checksum: $(flag_info.md5_checksum)") diff --git a/src/variable.jl b/src/variable.jl index bd2757b..e4dfedf 100644 --- a/src/variable.jl +++ b/src/variable.jl @@ -18,6 +18,8 @@ end Base.size(var::CDFVariable) = var.dims +@inline majority(var::CDFVariable) = majority(var.parentdataset) + function dst_src_ranges(first, last, entry) overlap_first = max(first, entry.first) overlap_last = min(last, entry.last) diff --git a/test/comprehensive_test.jl b/test/comprehensive_test.jl index f09c4d0..9323938 100644 --- a/test/comprehensive_test.jl +++ b/test/comprehensive_test.jl @@ -64,13 +64,13 @@ const EXPECTED_VARIABLES = Dict( "var3d_counter" => ( shape = (3, 5, 10), # Shape adjusted for Julia column-major data_type = "CDF_DOUBLE", - values = reshape(0.0:(3 * 5 * 10 - 1), 3, 5, 10), + values = permutedims(reshape(0.0:(3 * 5 * 10 - 1), 5, 3, 10), (2, 1, 3)), attributes = Dict("attr1" => "attr1_value", "attr2" => "attr2_value"), ), "var5d_counter" => ( shape = (5, 4, 3, 2, 6), # Shape adjusted for Julia column-major data_type = "CDF_DOUBLE", - values = reshape(0.0:(5 * 4 * 3 * 2 * 6 - 1), 5, 4, 3, 2, 6), + values = permutedims(reshape(0.0:(5 * 4 * 3 * 2 * 6 - 1), 2, 3, 4, 5, 6), (4, 3, 2, 1, 5)), attributes = Dict{String, Any}(), ), "var3d" => ( @@ -113,12 +113,12 @@ const EXPECTED_VARIABLES = Dict( shape = (2, 2, 1), data_type = "CDF_CHAR", attributes = Dict{String, Any}(), - values = ["value[00]"; "value[01]";; "value[10]"; "value[11]";;;], + values = ["value[00]"; "value[10]";; "value[01]"; "value[11]";;;], ), "var4d_string" => ( shape = (3, 2, 2, 1), data_type = "CDF_CHAR", - values = ["value[000]" "value[011]"; "value[001]" "value[100]"; "value[010]" "value[101]";;; "value[110]" "value[201]"; "value[111]" "value[210]"; "value[200]" "value[211]";;;;], + values = ["value[000]"; "value[100]"; "value[200]" ;; "value[010]"; "value[110]"; "value[210]";;; "value[001]"; "value[101]"; "value[201]";; "value[011]"; "value[111]"; "value[211]";;;;], attributes = Dict{String, Any}(), ) ) @@ -162,6 +162,14 @@ end end end +@testset "Variable Slicing" begin + @test ds["var3d_counter"][1, 1, :] == 0:15:135 + @test ds["var3d_counter"][:, :, 2]' == reshape(15:29, 5, 3) + + @test ds["var3d_counter"][2, :, 2] == 20:24 + @test ds["var3d_counter"][:, 1, 2] == 15:5:25 +end + @testset "Variable Varying" begin @test is_record_varying(ds["var"]) @test is_record_varying(ds["var_recvary_string"]) diff --git a/test/perf_test.jl b/test/perf_test.jl index c1ea314..7a20880 100644 --- a/test/perf_test.jl +++ b/test/perf_test.jl @@ -33,11 +33,11 @@ b = [b0, b1, b12, b2, b30, b3, b4, b5] # ┌ Info: Benchmarks # │ b = # │ 8-element Vector{Chairmarks.Sample}: -# │ 539.550 ns (6 allocs: 528 bytes) -# │ 2.083 μs (20 allocs: 29.172 KiB) -# │ 2.000 μs (24 allocs: 29.328 KiB) -# │ 83.896 μs (3777 allocs: 162.469 KiB) -# │ 385.400 ns (7 allocs: 528 bytes) -# │ 9.586 ms (574 allocs: 31.655 MiB) -# │ 467.275 μs (100 allocs: 1.344 MiB) -# └ 20.855 μs (250 allocs: 12.484 KiB) \ No newline at end of file +# │ 537.500 ns (6 allocs: 544 bytes) +# │ 1.425 μs (6 allocs: 28.219 KiB) +# │ 1.317 μs (6 allocs: 28.219 KiB) +# │ 95.062 μs (3777 allocs: 163.078 KiB) +# │ 387.500 ns (7 allocs: 560 bytes) +# │ 10.175 ms (567 allocs: 31.647 MiB, 5.34% gc time) +# │ 462.467 μs (93 allocs: 1.337 MiB) +# └ 23.209 μs (250 allocs: 12.656 KiB) \ No newline at end of file