Make ImageHDU parametric

[giordano]

The parameters are the same as the data it represents: the type of the elements and the number of dimensions of the array. As a result, read(::ImageHDU) is type-stable.

Ref #119.

[kbarbary]

Following up on my comment in #119.

When reading an image or table from a file, a typical procedure is:

f = FITS('myfile.fits', 'r')   # f is a `FITS`: type stable.
hdu = f[1]  # hdu is a Union{ImageHDU, TableHDU, ASCIITableHDU}: not type stable.
data = read(hdu)  # data is Any : not type stable

If I understand correctly, this change would expand the set of types that hdu could be (on the second line), but it doesn't narrow the set of types that data could be, since the type of hdu is not known at compile time. That is, on the third line data will still probably be inferred to be Any (or maybe Array at best).

If one really wants the type of data to be known to the compiler at compile time, currently you can do this:

f = FITS('myfile.fits', 'r')
hdu = f[1]
data = read(hdu)::Array{Float32, 2}

and after this PR, you could alternatively do this:

f = FITS('myfile.fits', 'r')
hdu = f[1]::ImageHDU{Float32, 2}
data = read(hdu)

which doesn't seem like an improvement.

Typically, I think neither of these is really necessary for performance, as long as the next thing you do with data is wrapped in a function call:

perform_my_analysis(data)

without the type assertions, the perform_my_analysis call will be dynamic dispatch, but everything within that function can still be type inferred and fast.

In general, the compiler has no way of inferring what type of data lives in an arbitrary file on disk, so we'll always have either some amount of dynamic dispatch (but usually a negligible amount for performance), or require the programmer to inform the compiler about the data in the file via type assertions.

All that said, I could be missing something about the greater goal of this PR or about the type inference going on here.

[giordano]

When reading an image or table from a file, a typical procedure is:

f = FITS('myfile.fits', 'r')   # f is a `FITS`: type stable.
hdu = f[1]  # hdu is a Union{ImageHDU, TableHDU, ASCIITableHDU}: not type stable.
data = read(hdu)  # data is Any : not type stable

If I understand correctly, this change would expand the set of types that hdu could be (on the second line), but it doesn't narrow the set of types that data could be, since the type of hdu is not known at compile time. That is, on the third line data will still probably be inferred to be Any (or maybe Array at best).

No, hdu = f[1] is still type-unstable, but once we have it read(hdu) is type-stable.

julia> f = FITS("delta.fits")
File: delta.fits
Mode: "r" (read-only)
HDUs: Num  Name  Type   
      1          Image  
      2    RA    Image  
      3    DEC   Image  

julia> @code_warntype f[1]
Body::HDU
1 ── %1  = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %2  = (Base.getfield)(%1, :ptr)::Ptr{Nothing}
│    %3  = (Core.bitcast)(Core.UInt, %2)::UInt64
│    %4  = (%3 === 0x0000000000000000)::Bool
└───       goto #3 if not %4
2 ──       invoke FITSIO.Libcfitsio.error("attempt to access closed FITS file"::String)
└───       $(Expr(:unreachable))
3 ┄─       goto #4
4 ── %9  = (Base.getfield)(f, :hdus)::Dict{Int64,HDU}
│    %10 = invoke Base.ht_keyindex(%9::Dict{Int64,HDU}, _3::Int64)::Int64
│    %11 = (Base.sle_int)(0, %10)::Bool
└───       goto #9 if not %11
5 ── %13 = (Base.getfield)(f, :hdus)::Dict{Int64,HDU}
│    %14 = invoke Base.ht_keyindex(%13::Dict{Int64,HDU}, _3::Int64)::Int64
│    %15 = (Base.slt_int)(%14, 0)::Bool
└───       goto #7 if not %15
6 ── %17 = %new(Base.KeyError, i)::KeyError
│          (Base.throw)(%17)
└───       $(Expr(:unreachable))
7 ┄─ %20 = (Base.getfield)(%13, :vals)::Array{HDU,1}
│    %21 = (Base.arrayref)(false, %20, %14)::HDU
└───       goto #8
8 ──       return %21
9 ── %24 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %25 = (Base.getfield)(%24, :ptr)::Ptr{Nothing}
│    %26 = (Core.bitcast)(Core.UInt, %25)::UInt64
│    %27 = (%26 === 0x0000000000000000)::Bool
└───       goto #11 if not %27
10 ─       invoke FITSIO.Libcfitsio.error("attempt to access closed FITS file"::String)
└───       $(Expr(:unreachable))
11 ┄       goto #12
12 ─ %32 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %33 = invoke FITSIO.fits_get_num_hdus(%32::FITSIO.Libcfitsio.FITSFile)::Int32
│    %34 = (Core.sext_int)(Core.Int64, %33)::Int64
└───       goto #13
13 ─ %36 = (Base.slt_int)(%34, i)::Bool
└───       goto #15 if not %36
14 ─       invoke FITSIO.error("index out of bounds"::String)
└───       $(Expr(:unreachable))
15 ┄ %40 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %41 = invoke FITSIO.fits_movabs_hdu(%40::FITSIO.Libcfitsio.FITSFile, _3::Int64)::Symbol
│    %42 = (%41 === :image_hdu)::Bool
└───       goto #17 if not %42
16 ─ %44 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %45 = invoke FITSIO.ImageHDU(%44::FITSIO.Libcfitsio.FITSFile, _3::Int64)::ImageHDU{_1,_2} where _2 where _1
└───       goto #23
17 ─ %47 = (%41 === :binary_table)::Bool
└───       goto #19 if not %47
18 ─ %49 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %50 = %new(FITSIO.TableHDU, %49, i)::TableHDU
└───       goto #22
19 ─ %52 = (%41 === :ascii_table)::Bool
└───       goto #21 if not %52
20 ─ %54 = (Base.getfield)(f, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│    %55 = %new(FITSIO.ASCIITableHDU, %54, i)::ASCIITableHDU
└───       goto #22
21 ─       invoke FITSIO.error("bad HDU type"::String)
└───       $(Expr(:unreachable))
22 ┄ %59 = φ (#18 => %50, #20 => %55)::Union{ASCIITableHDU, TableHDU}
23 ┄ %60 = φ (#16 => %45, #22 => %59)::Union{ASCIITableHDU, TableHDU, ImageHDU{_1,_2} where _2 where _1}
│    %61 = (Base.getfield)(f, :hdus)::Dict{Int64,HDU}
│    %62 = Base.setindex!::Core.Compiler.Const(setindex!, false)
│    %63 = (isa)(%60, ASCIITableHDU)::Bool
└───       goto #25 if not %63
24 ─ %65 = π (%60, ASCIITableHDU)
│          invoke %62(%61::Dict{Int64,HDU}, %65::ASCIITableHDU, _3::Int64)
└───       goto #28
25 ─ %68 = (isa)(%60, TableHDU)::Bool
└───       goto #27 if not %68
26 ─ %70 = π (%60, TableHDU)
│          invoke %62(%61::Dict{Int64,HDU}, %70::TableHDU, _3::Int64)
└───       goto #28
27 ─       (Base.setindex!)(%61, %60, i)
└───       goto #28
28 ┄ %75 = (Base.getfield)(f, :hdus)::Dict{Int64,HDU}
│    %76 = invoke Base.ht_keyindex(%75::Dict{Int64,HDU}, _3::Int64)::Int64
│    %77 = (Base.slt_int)(%76, 0)::Bool
└───       goto #30 if not %77
29 ─ %79 = %new(Base.KeyError, i)::KeyError
│          (Base.throw)(%79)
└───       $(Expr(:unreachable))
30 ┄ %82 = (Base.getfield)(%75, :vals)::Array{HDU,1}
│    %83 = (Base.arrayref)(false, %82, %76)::HDU
└───       goto #31
31 ─       return %83

julia> hdu = f[1]
File: delta.fits
HDU: 1
Type: Image
Datatype: Float64
Datasize: (128, 128, 32)

julia> @code_warntype read(hdu)
Body::Array{Float64,3}
1 ─ %1  = (Base.getfield)(hdu, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│   %2  = (Base.getfield)(%1, :ptr)::Ptr{Nothing}
│   %3  = (Core.bitcast)(Core.UInt, %2)::UInt64
│   %4  = (%3 === 0x0000000000000000)::Bool
└──       goto #3 if not %4
2 ─       invoke FITSIO.Libcfitsio.error("attempt to access closed FITS file"::String)
└──       $(Expr(:unreachable))
3 ┄       goto #4
4 ─ %9  = (Base.getfield)(hdu, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│   %10 = (Base.getfield)(hdu, :ext)::Int64
│         invoke FITSIO.fits_movabs_hdu(%9::FITSIO.Libcfitsio.FITSFile, %10::Int64)
│   %12 = (Base.getfield)(hdu, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│   %13 = invoke FITSIO.fits_get_img_size(%12::FITSIO.Libcfitsio.FITSFile)::Array{Int64,1}
│   %14 = (Core._apply)(Array{Float64,3}, (array initializer with undefined values,), %13)::Array{Float64,3}
│   %15 = (Base.getfield)(hdu, :fitsfile)::FITSIO.Libcfitsio.FITSFile
│         (Base.arraysize)(%14, 1)
│         (Base.arraysize)(%14, 2)
│         (Base.arraysize)(%14, 3)
│   %19 = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Int64,1}, svec(Any, Int64), :(:ccall), 2, Array{Int64,1}, 3, 3))::Array{Int64,1}
│   %20 = invoke Base.fill!(%19::Array{Int64,1}, 1::Int64)::Array{Int64,1}
│   %21 = (Base.arraylen)(%14)::Int64
│         invoke FITSIO.Libcfitsio.fits_read_pix(%15::FITSIO.Libcfitsio.FITSFile, %20::Array{Int64,1}, %21::Int64, %14::Array{Float64,3})
└──       return %14

Following your example of perform_my_analysis, you can have a method of that function working on hdu instead of data, and it would be type-stable. Yes, this is not a big deal, and yes, there is no significant performance improvement in general.

[kbarbary]

Yes, @code_warntype read(hdu) returns Body::Array{Float64,3}. But in this expression, the type of hdu is known. What happens if you put the three lines in a function, like

f(filename) = read(FITS(filename)[1])

?

[giordano]

Yes, of course that's not going to be type-stable anyway.

[kbarbary]

Yeah, I guess my point is just to ask what is the benefit of making read(hdu) type stable, if the type of hdu is not going to be inferable anyway.

[giordano]

Ok, I guess this can be closed.