diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..d3b5f83 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,2 @@ +style = "yas" +format_markdown = true \ No newline at end of file diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml new file mode 100644 index 0000000..c3c9fe7 --- /dev/null +++ b/.github/workflows/Format.yml @@ -0,0 +1,12 @@ +name: Format suggestions +on: + pull_request: + # this argument is not required if you don't use the `suggestion-label` input + types: [ opened, reopened, synchronize, labeled, unlabeled ] +jobs: + code-style: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/julia-format@v3 + with: + version: '1' # Set `version` to '1.0.54' if you need to use JuliaFormatter.jl v1.0.54 (default: '1') \ No newline at end of file diff --git a/README.md b/README.md index f765d97..05031c3 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,4 @@ [![Coverage](https://codecov.io/gh/ptiede/ComradeBase.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ptiede/ComradeBase.jl) [![Build Status](https://github.com/ptiede/ComradeBase.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/ptiede/ComradeBase.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) -[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) diff --git a/docs/make.jl b/docs/make.jl index 12aa7ab..b92e739 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -4,21 +4,16 @@ using Documenter DocMeta.setdocmeta!(ComradeBase, :DocTestSetup, :(using ComradeBase); recursive=true) makedocs(; - modules=[ComradeBase], - authors="Paul Tiede and contributors", - repo="https://github.com/ptiede/ComradeBase.jl/blob/{commit}{path}#{line}", - sitename="ComradeBase.jl", - format=Documenter.HTML(; - prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://ptiede.github.io/ComradeBase.jl", - assets=String[], - ), - pages=[ - "Home" => "index.md", - ], -) + modules=[ComradeBase], + authors="Paul Tiede and contributors", + repo="https://github.com/ptiede/ComradeBase.jl/blob/{commit}{path}#{line}", + sitename="ComradeBase.jl", + format=Documenter.HTML(; + prettyurls=get(ENV, "CI", "false") == "true", + canonical="https://ptiede.github.io/ComradeBase.jl", + assets=String[],), + pages=["Home" => "index.md"],) deploydocs(; - repo="github.com/ptiede/ComradeBase.jl", - devbranch="main", -) + repo="github.com/ptiede/ComradeBase.jl", + devbranch="main",) diff --git a/ext/ComradeBaseEnzymeExt.jl b/ext/ComradeBaseEnzymeExt.jl index dcf919a..9600ae3 100644 --- a/ext/ComradeBaseEnzymeExt.jl +++ b/ext/ComradeBaseEnzymeExt.jl @@ -5,18 +5,24 @@ using Enzyme: @parallel const EnzymeThreads = ComradeBase.ThreadsEx{:Enzyme} -function ComradeBase.intensitymap_analytic!(img::IntensityMap{T, N, D, <:ComradeBase.AbstractRectiGrid{D, <:EnzymeThreads}}, s::ComradeBase.AbstractModel) where {T, N, D} +function ComradeBase.intensitymap_analytic!(img::IntensityMap{T,N,D, + <:ComradeBase.AbstractRectiGrid{D, + <:EnzymeThreads}}, + s::ComradeBase.AbstractModel) where {T,N,D} dx, dy = ComradeBase.pixelsizes(img) g = ComradeBase.domainpoints(img) f = Base.Fix1(ComradeBase.intensity_point, s) pimg = parent(img) @parallel for I in CartesianIndices(pimg) - pimg[I] = f(g[I])*dx*dy + pimg[I] = f(g[I]) * dx * dy end return nothing end -function ComradeBase.intensitymap_analytic!(img::UnstructuredMap{T, N, <:ComradeBase.UnstructuredDomain{D, <:EnzymeThreads}}, s::ComradeBase.AbstractModel) where {T, N, D} +function ComradeBase.intensitymap_analytic!(img::UnstructuredMap{T,N, + <:ComradeBase.UnstructuredDomain{D, + <:EnzymeThreads}}, + s::ComradeBase.AbstractModel) where {T,N,D} g = ComradeBase.domainpoints(img) f = Base.Fix1(ComradeBase.intensity_point, s) pimg = parent(img) @@ -26,7 +32,10 @@ function ComradeBase.intensitymap_analytic!(img::UnstructuredMap{T, N, <:Comrade return nothing end -function ComradeBase.visibilitymap_analytic!(vis::ComradeBase.FluxMap2{T,N,<:ComradeBase.AbstractSingleDomain{<:Any, <:EnzymeThreads}}, s::ComradeBase.AbstractModel) where {T,N} +function ComradeBase.visibilitymap_analytic!(vis::ComradeBase.FluxMap2{T,N, + <:ComradeBase.AbstractSingleDomain{<:Any, + <:EnzymeThreads}}, + s::ComradeBase.AbstractModel) where {T,N} dims = axisdims(vis) g = domainpoints(dims) f = Base.Fix1(ComradeBase.visibility_point, s) diff --git a/ext/ComradeBaseOhMyThreadsExt.jl b/ext/ComradeBaseOhMyThreadsExt.jl index cfade0f..9af8526 100644 --- a/ext/ComradeBaseOhMyThreadsExt.jl +++ b/ext/ComradeBaseOhMyThreadsExt.jl @@ -3,7 +3,10 @@ module ComradeBaseOhMyThreadsExt using ComradeBase using OhMyThreads -function ComradeBase.intensitymap_analytic!(img::IntensityMap{T,N,D,<:ComradeBase.AbstractRectiGrid{D, <:OhMyThreads.Scheduler}}, s::ComradeBase.AbstractModel) where {T,N,D} +function ComradeBase.intensitymap_analytic!(img::IntensityMap{T,N,D, + <:ComradeBase.AbstractRectiGrid{D, + <:OhMyThreads.Scheduler}}, + s::ComradeBase.AbstractModel) where {T,N,D} dims = axisdims(img) dx = step(dims.X) dy = step(dims.Y) @@ -11,30 +14,35 @@ function ComradeBase.intensitymap_analytic!(img::IntensityMap{T,N,D,<:ComradeBas pimg = parent(img) f = Base.Fix1(ComradeBase.intensity_point, s) tforeach(CartesianIndices(pimg); scheduler=executor(dims)) do I - pimg[I] = f(g[I])*dx*dy + return pimg[I] = f(g[I]) * dx * dy end return nothing end - -function ComradeBase.intensitymap_analytic!(img::UnstructuredMap{T,<:AbstractVector,<:UnstructuredDomain{D, <:OhMyThreads.Scheduler}}, s::ComradeBase.AbstractModel) where {T,D} +function ComradeBase.intensitymap_analytic!(img::UnstructuredMap{T,<:AbstractVector, + <:UnstructuredDomain{D, + <:OhMyThreads.Scheduler}}, + s::ComradeBase.AbstractModel) where {T,D} dims = axisdims(img) g = domainpoints(dims) f = Base.Fix1(ComradeBase.intensity_point, s) pimg = parent(img) tforeach(CartesianIndices(pimg); scheduler=executor(dims)) do I - pimg[I] = f(g[I]) + return pimg[I] = f(g[I]) end return nothing end -function ComradeBase.visibilitymap_analytic!(vis::ComradeBase.FluxMap2{T,N,<:ComradeBase.AbstractSingleDomain{<:Any, <: OhMyThreads.Scheduler}}, s::ComradeBase.AbstractModel) where {T,N} +function ComradeBase.visibilitymap_analytic!(vis::ComradeBase.FluxMap2{T,N, + <:ComradeBase.AbstractSingleDomain{<:Any, + <:OhMyThreads.Scheduler}}, + s::ComradeBase.AbstractModel) where {T,N} dims = axisdims(vis) g = domainpoints(dims) f = Base.Fix1(ComradeBase.visibility_point, s) pvis = parent(vis) tforeach(CartesianIndices(pvis); scheduler=executor(dims)) do I - pvis[I] = f(g[I]) + return pvis[I] = f(g[I]) end return nothing end diff --git a/scratch/moments.jl b/scratch/moments.jl index cf0ed74..baaeb25 100644 --- a/scratch/moments.jl +++ b/scratch/moments.jl @@ -4,13 +4,11 @@ Computes the flux of a intensity map """ function flux(im::IntensityMap{T,N}) where {T,N} - return sum(im, dims=(:X, :Y)) + return sum(im; dims=(:X, :Y)) end flux(im::IntensityMap{T,2}) where {T} = sum(im) - - """ centroid(im::IntensityMap) @@ -18,22 +16,21 @@ Computes the image centroid aka the center of light of the image. """ function centroid(im::IntensityMap{T,N}) where {T<:Number,N} (X, Y) = named_axiskeys(im) - return mapslices(x->centroid(IntensityMap(x, (;X, Y))), im; dims=(:X, :Y)) + return mapslices(x -> centroid(IntensityMap(x, (; X, Y))), im; dims=(:X, :Y)) end function centroid(im::IntensityMap{T,2})::Tuple{T,T} where {T<:Number} x0 = y0 = zero(T) f = flux(im) - @inbounds for (i,x) in pairs(axiskeys(im,:X)), (j,y) in pairs(axiskeys(im,:Y)) - x0 += x.*im[X=i, Y=j] - y0 += y.*im[X=i, Y=j] + @inbounds for (i, x) in pairs(axiskeys(im, :X)), (j, y) in pairs(axiskeys(im, :Y)) + x0 += x .* im[X=i, Y=j] + y0 += y .* im[X=i, Y=j] end - return x0/f, y0/f + return x0 / f, y0 / f end centroid(im::IntensityMap{<:StokesParams}) = centroid(stokes(im, :I)) - """ second_moment(im::IntensityMap; center=true) @@ -43,28 +40,31 @@ second moment, which is specified by the `center` argument. """ function second_moment(im::IntensityMap{T,N}; center=true) where {T<:Number,N} (X, Y) = named_axiskeys(im) - return mapslices(x->second_moment(IntensityMap(x, (;X, Y)); center), im; dims=(:X, :Y)) + return mapslices(x -> second_moment(IntensityMap(x, (; X, Y)); center), im; + dims=(:X, :Y)) end -second_moment(im::IntensityMap{<:StokesParams}; center=true) = second_moment(stokes(im, :I); center=true) +function second_moment(im::IntensityMap{<:StokesParams}; center=true) + return second_moment(stokes(im, :I); center=true) +end function second_moment(im::IntensityMap{T,2}; center=true) where {T<:Number} xx = zero(T) xy = zero(T) yy = zero(T) f = flux(im) - (;X,Y) = named_axiskeys(im) - for (i,y) in pairs(Y), (j,x) in pairs(X) - xx += x^2*im[j,i] - yy += y^2*im[j,i] - xy += x*y*im[j,i] + (; X, Y) = named_axiskeys(im) + for (i, y) in pairs(Y), (j, x) in pairs(X) + xx += x^2 * im[j, i] + yy += y^2 * im[j, i] + xy += x * y * im[j, i] end if center x0, y0 = centroid(im) xx -= x0^2 yy -= y0^2 - xy -= x0*y0 + xy -= x0 * y0 end return @SMatrix [xx/f xy/f; xy/f yy/f] diff --git a/scratch/stokes_image_old.jl b/scratch/stokes_image_old.jl index 529ca7a..9ece7bf 100644 --- a/scratch/stokes_image_old.jl +++ b/scratch/stokes_image_old.jl @@ -1,4 +1,4 @@ -const StokesIntensityMap{T,N} = KeyedArray{StokesParams{T}, N, <:DataArr} where {T, N} +const StokesIntensityMap{T,N} = KeyedArray{StokesParams{T},N,<:DataArr} where {T,N} baseimage(s::KeyedArray) = parent(parent(s)) @@ -8,9 +8,11 @@ baseimage(s::KeyedArray) = parent(parent(s)) Create an array of full stokes parameters. The image is stored as a `StructArray` of @ref[StokesParams]. """ -function StokesIntensityMap(I::IntensityMap{T}, Q::IntensityMap{T}, U::IntensityMap{T}, V::IntensityMap{T}) where {T} - @assert check_grid(I,Q,U,V) "Intensity grids are not the same across the 4 stokes parameters" - polimg = StructArray{StokesParams{T}}(I=baseimage(I), Q=baseimage(Q), U=baseimage(U), V=baseimage(V)) +function StokesIntensityMap(I::IntensityMap{T}, Q::IntensityMap{T}, U::IntensityMap{T}, + V::IntensityMap{T}) where {T} + @assert check_grid(I, Q, U, V) "Intensity grids are not the same across the 4 stokes parameters" + polimg = StructArray{StokesParams{T}}(; I=baseimage(I), Q=baseimage(Q), U=baseimage(U), + V=baseimage(V)) return IntensityMap(polimg, named_axiskeys(I)) end @@ -20,15 +22,17 @@ Helper function that converts a model from something that compute polarized imag to just a single stokes parameter. This is useful if you just want to fit a single stokes parameter. """ -struct SingleStokes{M, S} <: ComradeBase.AbstractModel +struct SingleStokes{M,S} <: ComradeBase.AbstractModel model::M end -SingleStokes(m::M, param::Symbol) where {M} = SingleStokes{M, param}(m) +SingleStokes(m::M, param::Symbol) where {M} = SingleStokes{M,param}(m) visanalytic(::Type{SingleStokes{M,S}}) where {M,S} = visanalytic(M) -imanalytic(::Type{SingleStokes{M,S}}) where {M,S} = imanalytic(M) -@inline intensity_point(s::SingleStokes{M,S}, x,y) where {M,S} = getproperty(intensity_point(s.model, x,y), S) - +imanalytic(::Type{SingleStokes{M,S}}) where {M,S} = imanalytic(M) +@inline intensity_point(s::SingleStokes{M,S}, x, y) where {M,S} = getproperty(intensity_point(s.model, + x, + y), + S) # # not type piracy because I own StokesParams # function Base.getproperty(s::StokesIntensityMap{T,N,Na}, v::Symbol) where {T,N,Na<:DataArr} @@ -41,28 +45,23 @@ imanalytic(::Type{SingleStokes{M,S}}) where {M,S} = imanalytic(M) # end # end -function StokesIntensityMap( - I::AbstractArray{T,N}, Q::AbstractArray{T,N}, - U::AbstractArray{T,N}, V::AbstractArray{T,N}, - dims::NamedTuple{Na,<:NTuple{N,Any}}) where {T, N, Na} - - polimg = StructArray{StokesParams{T}}((;I,Q,U,V)) +function StokesIntensityMap(I::AbstractArray{T,N}, Q::AbstractArray{T,N}, + U::AbstractArray{T,N}, V::AbstractArray{T,N}, + dims::NamedTuple{Na,<:NTuple{N,Any}}) where {T,N,Na} + polimg = StructArray{StokesParams{T}}((; I, Q, U, V)) return IntensityMap(polimg, dims) end -function StokesIntensityMap( - I::AbstractArray{T,N}, Q::AbstractArray{T,N}, - U::AbstractArray{T,N}, V::AbstractArray{T,N}, - args...) where {T, N} - - polimg = StructArray{StokesParams{T}}((;I,Q,U,V)) +function StokesIntensityMap(I::AbstractArray{T,N}, Q::AbstractArray{T,N}, + U::AbstractArray{T,N}, V::AbstractArray{T,N}, + args...) where {T,N} + polimg = StructArray{StokesParams{T}}((; I, Q, U, V)) return IntensityMap(polimg, args...) end - # simple check to ensure that the four grids are equal across stokes parameters -function check_grid(I,Q,U,V) - named_axiskeys(I) == named_axiskeys(Q) == named_axiskeys(U) == named_axiskeys(V) +function check_grid(I, Q, U, V) + return named_axiskeys(I) == named_axiskeys(Q) == named_axiskeys(U) == named_axiskeys(V) end function stokes(pimg::StokesIntensityMap, v::Symbol) @@ -79,12 +78,13 @@ function _summary(io, x::StokesIntensityMap{T,N,Na}) where {T,N,Na} println(io, ndims(x), "-dimensional") println(io, "StokesIntensityMap{$T, $N, $Na}") for d in 1:ndims(x) - print(io, d==1 ? "↓" : d==2 ? "→" : d==3 ? "◪" : "▨", " ") + print(io, d == 1 ? "↓" : d == 2 ? "→" : d == 3 ? "◪" : "▨", " ") c = AxisKeys.colour(x, d) - AxisKeys.hasnames(x) && printstyled(io, AxisKeys.dimnames(x,d), " ∈ ", color=c) - printstyled(io, length(axiskeys(x,d)), "-element ", AxisKeys.shorttype(axiskeys(x,d)), "\n", color=c) + AxisKeys.hasnames(x) && printstyled(io, AxisKeys.dimnames(x, d), " ∈ "; color=c) + printstyled(io, length(axiskeys(x, d)), "-element ", + AxisKeys.shorttype(axiskeys(x, d)), "\n"; color=c) end - println(io, "Polarizations ", propertynames(parent(x).data)) + return println(io, "Polarizations ", propertynames(parent(x).data)) end Base.show(io::IO, img::StokesIntensityMap) = summary(io, img) diff --git a/src/ComradeBase.jl b/src/ComradeBase.jl index 9eb8a81..a2fc5d8 100644 --- a/src/ComradeBase.jl +++ b/src/ComradeBase.jl @@ -11,20 +11,20 @@ using Reexport @reexport using PolarizedTypes using PrecompileTools -export visibility, - intensitymap, intensitymap!, - visibilitymap, visibilitymap!, - StokesParams, CoherencyMatrix, - flux, fieldofview, imagepixels, pixelsizes, IntensityMap, - named_dims +export visibility, + intensitymap, intensitymap!, + visibilitymap, visibilitymap!, + StokesParams, CoherencyMatrix, + flux, fieldofview, imagepixels, pixelsizes, IntensityMap, + named_dims include("interface.jl") include("domain.jl") include("unstructured_map.jl") include("images/images.jl") -const FluxMap2{T, N, E} = Union{IntensityMap{T,N,<:Any,E}, UnstructuredMap{T,<:AbstractVector,E}} - +const FluxMap2{T,N,E} = Union{IntensityMap{T,N,<:Any,E}, + UnstructuredMap{T,<:AbstractVector,E}} include("visibilities.jl") # include("rrules.jl") @@ -39,20 +39,8 @@ include("visibilities.jl") g = RectiGrid(p) gs = domainpoints(p) imgI = IntensityMap(rand(10, 10), g) - imgI.^2 - end -end - -if !isdefined(Base, :get_extension) - using Requires -end - -@static if !isdefined(Base, :get_extension) - function __init__() - @require OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" include(joinpath(@__DIR__, "..", "ext", "ComradeBaseOhMyThreadsExt.jl")) + imgI .^ 2 end end - - end diff --git a/src/domain.jl b/src/domain.jl index 58b9fc1..3af3c04 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -3,9 +3,8 @@ export RectiGrid, UnstructuredDomain, domainpoints, named_dims, dims, header, axisdims, executor, Serial, ThreadsEx - abstract type AbstractDomain end -abstract type AbstractSingleDomain{D, E} <: AbstractDomain end +abstract type AbstractSingleDomain{D,E} <: AbstractDomain end """ create_map(array, g::AbstractSingleDomain) @@ -45,7 +44,7 @@ Allocate the default map specialized by the grid `g` function allocate_vismap end function allocate_vismap(m::AbstractModel, g::AbstractDomain) - allocate_vismap(ispolarized(typeof(m)), m, g) + return allocate_vismap(ispolarized(typeof(m)), m, g) end function allocate_vismap(::IsPolarized, m::AbstractModel, g::AbstractDomain) @@ -59,7 +58,7 @@ function allocate_vismap(::NotPolarized, m::AbstractModel, g::AbstractDomain) end function allocate_imgmap(m::AbstractModel, g::AbstractDomain) - allocate_imgmap(ispolarized(typeof(m)), m, g) + return allocate_imgmap(ispolarized(typeof(m)), m, g) end function allocate_imgmap(::IsPolarized, m::AbstractModel, g::AbstractDomain) @@ -72,10 +71,6 @@ function allocate_imgmap(::NotPolarized, m::AbstractModel, g::AbstractDomain) return allocate_map(M, g) end - - - - """ domainpoints(g::AbstractSingleDomain) @@ -89,7 +84,6 @@ function domainpoints end # ChainRulesCore.@non_differentiable domainpoints(d::AbstractSingleDomain) # EnzymeRules.inactive(::typeof(domainpoints), args...) = nothing - """ executor(g::AbstractSingleDomain) @@ -99,7 +93,6 @@ executor(g::AbstractSingleDomain) = getfield(g, :executor) # ChainRulesCore.@non_differentiable executor(::AbstractSingleDomain) EnzymeRules.inactive(::typeof(executor), args...) = nothing - """ dims(g::AbstractSingleDomain) @@ -109,7 +102,6 @@ DD.dims(g::AbstractSingleDomain) = getfield(g, :dims) # ChainRulesCore.@non_differentiable DD.dims(::AbstractSingleDomain) EnzymeRules.inactive(::typeof(DD.dims), x::AbstractSingleDomain) = nothing - """ named_dims(g::AbstractSingleDomain) @@ -119,7 +111,6 @@ named_dims(g::AbstractSingleDomain) = NamedTuple{keys(g)}(dims(g)) # ChainRulesCore.@non_differentiable named_dims(::AbstractSingleDomain) EnzymeRules.inactive(::typeof(named_dims), args...) = nothing - """ header(g::AbstractSingleDomain) @@ -128,7 +119,9 @@ Returns the headerinformation of the dimensions `g` header(g::AbstractSingleDomain) = getfield(g, :header) # ChainRulesCore.@non_differentiable header(::AbstractSingleDomain) EnzymeRules.inactive(::typeof(header), args...) = nothing -Base.keys(g::AbstractSingleDomain) = throw(MethodError(Base.keys, "You must implement `Base.keys($(typeof(g)))`")) +function Base.keys(g::AbstractSingleDomain) + throw(MethodError(Base.keys, "You must implement `Base.keys($(typeof(g)))`")) +end """ Serial() @@ -156,18 +149,16 @@ else const schedulers = (:(:dynamic), :(:static)) end - # We index the dimensions not the grid itself Base.getindex(d::AbstractSingleDomain, i::Int) = getindex(dims(d), i) - Base.ndims(d::AbstractSingleDomain) = length(dims(d)) Base.size(d::AbstractSingleDomain) = map(length, dims(d)) Base.length(d::AbstractSingleDomain) = prod(size(d)) Base.firstindex(d::AbstractSingleDomain) = 1 Base.lastindex(d::AbstractSingleDomain) = length(d) Base.axes(d::AbstractSingleDomain) = axes(dims(d)) -Base.iterate(d::AbstractSingleDomain, i::Int = 1) = iterate(dims(d), i) +Base.iterate(d::AbstractSingleDomain, i::Int=1) = iterate(dims(d), i) # Base.front(d::AbstractSingleDomain) = Base.front(dims(d)) # We return the eltype of the dimensions. Should we change this? Base.eltype(d::AbstractSingleDomain) = eltype(basedim(first(dims(d)))) @@ -177,7 +168,6 @@ Base.eltype(d::AbstractSingleDomain) = eltype(basedim(first(dims(d)))) # Base.map(f, args, d::AbstractSingleDomain) = map(f, args, dims(d)) # Base.map(f, d::AbstractSingleDomain, args) = rebuild(typeof(d), map(f, dims(d), args)) - abstract type AbstractHeader end """ @@ -223,19 +213,19 @@ end """ struct NoHeader <: AbstractHeader end - -abstract type AbstractRectiGrid{D, E} <: AbstractSingleDomain{D, E} end +abstract type AbstractRectiGrid{D,E} <: AbstractSingleDomain{D,E} end create_map(array, g::AbstractRectiGrid) = IntensityMap(array, g) -allocate_map(M::Type{<:AbstractArray}, g::AbstractRectiGrid) = IntensityMap(similar(M, size(g)), g) +function allocate_map(M::Type{<:AbstractArray}, g::AbstractRectiGrid) + return IntensityMap(similar(M, size(g)), g) +end function fieldofview(dims::AbstractRectiGrid) - (;X,Y) = dims + (; X, Y) = dims dx = step(X) dy = step(Y) - (X=abs(last(X) - first(X))+dx, Y=abs(last(Y)-first(Y))+dy) + return (X=abs(last(X) - first(X)) + dx, Y=abs(last(Y) - first(Y)) + dy) end - """ pixelsizes(img::IntensityMap) pixelsizes(img::AbstractRectiGrid) @@ -248,28 +238,25 @@ function pixelsizes(keys::AbstractRectiGrid) return (X=step(x), Y=step(y)) end - -struct RectiGrid{D, E, Hd<:AbstractHeader} <: AbstractRectiGrid{D, E} +struct RectiGrid{D,E,Hd<:AbstractHeader} <: AbstractRectiGrid{D,E} dims::D executor::E header::Hd - @inline function RectiGrid(dims::Tuple; executor=Serial(), header::AbstractHeader=NoHeader()) + @inline function RectiGrid(dims::Tuple; executor=Serial(), + header::AbstractHeader=NoHeader()) df = _format_dims(dims) - return new{typeof(df), typeof(executor), typeof(header)}(df, executor, header) + return new{typeof(df),typeof(executor),typeof(header)}(df, executor, header) end - end EnzymeRules.inactive_type(::Type{<:RectiGrid}) = true - -function domainpoints(d::RectiGrid{D, Hd}) where {D, Hd} +function domainpoints(d::RectiGrid{D,Hd}) where {D,Hd} g = map(basedim, dims(d)) N = keys(d) return StructArray(NamedTuple{N}(_build_slices(g, size(d)))) end - function _format_dims(dg::Tuple) return DD.format(dg, map(eachindex, dg)) end @@ -278,34 +265,32 @@ Base.keys(g::RectiGrid) = map(name, dims(g)) @inline RectiGrid(g::RectiGrid) = g - -function Base.show(io::IO, mime::MIME"text/plain", x::RectiGrid{D, E}) where {D, E} +function Base.show(io::IO, mime::MIME"text/plain", x::RectiGrid{D,E}) where {D,E} println(io, "RectiGrid(") println(io, "executor: $(executor(x))") println(io, "Dimensions: ") show(io, mime, dims(x)) - print(io, "\n)") + return print(io, "\n)") end - Base.propertynames(d::RectiGrid) = keys(d) Base.getproperty(g::RectiGrid, p::Symbol) = basedim(dims(g)[findfirst(==(p), keys(g))]) - - # This is needed to prevent doubling up on the dimension -@inline function RectiGrid(dims::NamedTuple{Na, T}; executor=Serial(), header::AbstractHeader=NoHeader()) where {Na, N, T<:NTuple{N, DD.Dimension}} +@inline function RectiGrid(dims::NamedTuple{Na,T}; executor=Serial(), + header::AbstractHeader=NoHeader()) where {Na,N, + T<:NTuple{N, + DD.Dimension}} return RectiGrid(values(dims); executor, header) end @noinline function _make_dims(ks, vs) ds = DD.name2dim(ks) - return map(ds, vs) do d,v - DD.rebuild(d, v) + return map(ds, vs) do d, v + return DD.rebuild(d, v) end end - """ RectiGrid(dims::NamedTuple{Na}; executor=Serial(), header=ComradeBase.NoHeader()) RectiGrid(dims::NTuple{N, <:DimensionalData.Dimension}; executor=Serial(), header=ComradeBase.NoHeader()) @@ -338,25 +323,24 @@ dims = RectiGrid((X(-5.0:0.1:5.0), Y(-4.0:0.1:4.0), Ti([1.0, 1.5, 1.75]), Fr([23 dims = RectiGrid((X = -5.0:0.1:5.0, Y = -4.0:0.1:4.0, Ti = [1.0, 1.5, 1.75], Fr = [230, 345]); executor=ThreadsEx())) ``` """ -@inline function RectiGrid(nt::NamedTuple; executor=Serial(), header::AbstractHeader=ComradeBase.NoHeader()) +@inline function RectiGrid(nt::NamedTuple; executor=Serial(), + header::AbstractHeader=ComradeBase.NoHeader()) dims = _make_dims(keys(nt), values(nt)) return RectiGrid(dims; executor, header) end -function DD.rebuild(::Type{<:RectiGrid}, g, executor=Serial(), header=ComradeBase.NoHeader()) - RectiGrid(g; executor, header) +function DD.rebuild(::Type{<:RectiGrid}, g, executor=Serial(), + header=ComradeBase.NoHeader()) + return RectiGrid(g; executor, header) end - # Define some helpful names for ease typing -const DataNames = Union{<:NamedTuple{(:X, :Y, :T, :F)}, <:NamedTuple{(:X, :Y, :F, :T)}, - <:NamedTuple{(:X, :Y, :T)}, <:NamedTuple{(:X, :Y, :F)}, - <:NamedTuple{(:X,:Y)}} - - +const DataNames = Union{<:NamedTuple{(:X, :Y, :T, :F)},<:NamedTuple{(:X, :Y, :F, :T)}, + <:NamedTuple{(:X, :Y, :T)},<:NamedTuple{(:X, :Y, :F)}, + <:NamedTuple{(:X, :Y)}} # TODO make this play nice with dimensional data -struct UnstructuredDomain{D,E, H<:AbstractHeader} <: AbstractSingleDomain{D,E} +struct UnstructuredDomain{D,E,H<:AbstractHeader} <: AbstractSingleDomain{D,E} dims::D executor::E header::H @@ -364,7 +348,6 @@ end EnzymeRules.inactive_type(::Type{<:UnstructuredDomain}) = true - """ UnstructuredDomain(dims::NamedTuple; executor=Serial(), header=ComradeBase.NoHeader) @@ -389,8 +372,9 @@ Base.lastindex(d::UnstructuredDomain) = lastindex(dims(d)) # Base.front(d::UnstructuredDomain) = UnstructuredDomain(Base.front(StructArrays.components(dims(d))), executor=executor(d), header=header(d)) # Base.eltype(d::UnstructuredDomain) = Base.eltype(dims(d)) -function DD.rebuild(::Type{<:UnstructuredDomain}, g, executor=Serial(), header=ComradeBase.NoHeader()) - UnstructuredDomain(g, executor, header) +function DD.rebuild(::Type{<:UnstructuredDomain}, g, executor=Serial(), + header=ComradeBase.NoHeader()) + return UnstructuredDomain(g, executor, header) end Base.propertynames(g::UnstructuredDomain) = propertynames(domainpoints(g)) @@ -419,7 +403,7 @@ end function Base.summary(io::IO, g::UnstructuredDomain) n = propertynames(domainpoints(g)) printstyled(io, "│ "; color=:light_black) - print(io, "UnstructuredDomain with dims: $n") + return print(io, "UnstructuredDomain with dims: $n") end function Base.show(io::IO, mime::MIME"text/plain", x::UnstructuredDomain) @@ -427,9 +411,10 @@ function Base.show(io::IO, mime::MIME"text/plain", x::UnstructuredDomain) println(io, "executor: $(executor(x))") println(io, "Dimensions: ") show(io, mime, dims(x)) - print(io, "\n)") + return print(io, "\n)") end - create_map(array, g::UnstructuredDomain) = UnstructuredMap(array, g) -allocate_map(M::Type{<:A}, g::UnstructuredDomain) where {A<:AbstractArray} = UnstructuredMap(similar(M, size(g)), g) +function allocate_map(M::Type{<:A}, g::UnstructuredDomain) where {A<:AbstractArray} + return UnstructuredMap(similar(M, size(g)), g) +end diff --git a/src/images/dim_image.jl b/src/images/dim_image.jl index deb0e57..58dddc7 100644 --- a/src/images/dim_image.jl +++ b/src/images/dim_image.jl @@ -3,10 +3,9 @@ const DD = DimensionalData using DimensionalData: AbstractDimArray, NoName, NoMetadata, format, DimTuple, Dimension, XDim, YDim, ZDim, X, Y, Ti - DD.@dim Fr ZDim "frequency" -DD.@dim U XDim "U" -DD.@dim V YDim "V" +DD.@dim U XDim "U" +DD.@dim V YDim "V" export IntensityMap, Fr, X, Y, Ti, U, V @@ -42,25 +41,25 @@ julia> img3 = IntensityMap(data, 10.0, 10.0; header=NoHeader()) Broadcasting, map, and reductions should all just obey the `DimensionalData` interface. """ -struct IntensityMap{T,N,D,G<:AbstractRectiGrid{D},A<:AbstractArray{T,N},R<:Tuple,Na} <: AbstractDimArray{T,N,D,A} +struct IntensityMap{T,N,D,G<:AbstractRectiGrid{D},A<:AbstractArray{T,N},R<:Tuple,Na} <: + AbstractDimArray{T,N,D,A} data::A grid::G refdims::R name::Na - function IntensityMap( - data::A, grid::G, refdims::R, name::Na - ) where {A<:AbstractArray{T,N}, G<:AbstractRectiGrid{D}, R<:Tuple, Na} where {T,N,D} - - new{T,N,D,G,A,R,Na}(data, grid, refdims, name) + function IntensityMap(data::A, grid::G, refdims::R, + name::Na) where {A<:AbstractArray{T,N},G<:AbstractRectiGrid{D}, + R<:Tuple,Na} where {T,N,D} + return new{T,N,D,G,A,R,Na}(data, grid, refdims, name) end end -DD.dims(img::IntensityMap) = dims(getfield(img, :grid)) +DD.dims(img::IntensityMap) = dims(getfield(img, :grid)) DD.refdims(img::IntensityMap) = getfield(img, :refdims) -DD.data(img::IntensityMap) = getfield(img, :data) -DD.name(img::IntensityMap) = getfield(img, :name) -DD.metadata(img::IntensityMap)= header(axisdims(img)) -executor(img::IntensityMap) = executor(axisdims(img)) +DD.data(img::IntensityMap) = getfield(img, :data) +DD.name(img::IntensityMap) = getfield(img, :name) +DD.metadata(img::IntensityMap) = header(axisdims(img)) +executor(img::IntensityMap) = executor(axisdims(img)) # For the `IntensityMap` nothing is AD-able except the data so # let's tell Enzyme this @@ -71,7 +70,8 @@ EnzymeRules.inactive(::typeof(DD.metadata), ::IntensityMap) = nothing EnzymeRules.inactive(::typeof(executor), ::IntensityMap) = nothing @inline function stokes(pimg::IntensityMap{<:StokesParams}, v::Symbol) - IntensityMap(stokes(baseimage(pimg), v), axisdims(pimg), refdims(pimg), name(pimg)) + return IntensityMap(stokes(baseimage(pimg), v), axisdims(pimg), refdims(pimg), + name(pimg)) end function Base.propertynames(img::IntensityMap) @@ -86,8 +86,8 @@ function Base.getproperty(img::IntensityMap, p::Symbol) end end -const SpatialDims = Tuple{<:DD.Dimensions.X, <:DD.Dimensions.Y} -const SpatialIntensityMap{T, A, G} = IntensityMap{T,2,<:SpatialDims, A, G} where {T,A,G} +const SpatialDims = Tuple{<:DD.Dimensions.X,<:DD.Dimensions.Y} +const SpatialIntensityMap{T,A,G} = IntensityMap{T,2,<:SpatialDims,A,G} where {T,A,G} """ IntensityMap(data::AbstractArray, g::AbstractRectiGrid; refdims=(), name=Symbol("")) @@ -95,7 +95,8 @@ const SpatialIntensityMap{T, A, G} = IntensityMap{T,2,<:SpatialDims, A, G} where Creates a IntensityMap with the pixel fluxes `data` on the grid `g`. Optionally, you can specify a set of reference dimensions `refdims` as a tuple and a name for array `name`. """ -function IntensityMap(data::AbstractArray, g::AbstractRectiGrid; refdims=(), name=Symbol("")) +function IntensityMap(data::AbstractArray, g::AbstractRectiGrid; refdims=(), + name=Symbol("")) return IntensityMap(data, g, (), Symbol("")) end @@ -105,18 +106,17 @@ end Creates a IntensityMap with the pixel fluxes `data` and a spatial grid with field of view (`fovx`, `fovy`) and center pixel offset (`x0`, `y0`) and header `header`. """ -function IntensityMap(data::AbstractArray{T}, fovx::Real, fovy::Real, x0::Real=0, y0::Real=0; header=NoHeader()) where {T} +function IntensityMap(data::AbstractArray{T}, fovx::Real, fovy::Real, x0::Real=0, + y0::Real=0; header=NoHeader()) where {T} grid = imagepixels(fovx, fovy, size(data)..., T(x0), T(y0); header) return IntensityMap(data, grid) end - function IntensityMap(data::IntensityMap, g::AbstractRectiGrid) @assert g == axisdims(data) "Dimensions do not agree" return data end - """ axisdims(img::IntensityMap) axisdims(img::IntensityMap, p::Symbol) @@ -129,7 +129,6 @@ axisdims(img::IntensityMap, p::Symbol) = getproperty(axisdims(img), p) EnzymeRules.inactive(::typeof(axisdims), args...) = nothing named_dims(img::IntensityMap) = named_dims(axisdims(img)) - """ header(img::IntensityMap) @@ -139,17 +138,14 @@ header(img::IntensityMap) = header(axisdims(img)) DD._noname(::IntensityMap) = Symbol("") - Base.parent(img::IntensityMap) = DD.data(img) baseimage(x::IntensityMap) = baseimage(parent(x)) -@inline function DD.rebuild( - img::IntensityMap, data, dims::Tuple = dims(img), - refdims = refdims(img), - n = name(img), - metadata = metadata(img), - ) +@inline function DD.rebuild(img::IntensityMap, data, dims::Tuple=dims(img), + refdims=refdims(img), + n=name(img), + metadata=metadata(img)) # @info (typeof(img)) # TODO find why Name is changing type # n2 = n == Symbol("") ? NoName : n @@ -159,37 +155,35 @@ baseimage(x::IntensityMap) = baseimage(parent(x)) return IntensityMap(data, grid, refdims, n) end -@inline function DD.rebuild( - img::IntensityMap; data=DD.data(img), dims::Tuple = dims(img), - refdims = refdims(img), - name = name(img), - metadata = metadata(img), - ) - rebuild(img, data, dims, refdims, name, metadata) +@inline function DD.rebuild(img::IntensityMap; data=DD.data(img), dims::Tuple=dims(img), + refdims=refdims(img), + name=name(img), + metadata=metadata(img),) + return rebuild(img, data, dims, refdims, name, metadata) end function intensitymap_analytic!(img::IntensityMap, s::AbstractModel) dx, dy = pixelsizes(img) g = domainpoints(img) - img .= intensity_point.(Ref(s), g).*dx.*dy + img .= intensity_point.(Ref(s), g) .* dx .* dy return nothing end -function intensitymap_analytic!( - img::IntensityMap{T,N,D,<:ComradeBase.AbstractRectiGrid{D, <:ThreadsEx{S}}}, - s::AbstractModel) where {T,N,D,S} +function intensitymap_analytic!(img::IntensityMap{T,N,D, + <:ComradeBase.AbstractRectiGrid{D, + <:ThreadsEx{S}}}, + s::AbstractModel) where {T,N,D,S} g = domainpoints(img) _threads_intensitymap!(img, s, g, Val(S)) return nothing end - for s in schedulers @eval begin function _threads_intensitymap!(img::IntensityMap, s::AbstractModel, g, ::Val{$s}) dx, dy = pixelsizes(img) Threads.@threads $s for I in CartesianIndices(g) - img[I] = intensity_point(s, g[I])*dx*dy + img[I] = intensity_point(s, g[I]) * dx * dy end end end diff --git a/src/images/images.jl b/src/images/images.jl index 9d18ff1..c338abe 100644 --- a/src/images/images.jl +++ b/src/images/images.jl @@ -1,7 +1,7 @@ export IntensityMap, SpatialIntensityMap, - DataArr, SpatialDataArr, DataNames, - named_axisdims, imagepixels, pixelsizes, domainpoints, - phasecenter, baseimage, stokes + DataArr, SpatialDataArr, DataNames, + named_axisdims, imagepixels, pixelsizes, domainpoints, + phasecenter, baseimage, stokes include("dim_image.jl") @@ -9,7 +9,6 @@ export flux, centroid, second_moment, named_axisdims, axisdims, imagepixels, pixelsizes, domainpoints, phasecenter include("methods.jl") - """ intensitymap(model::AbstractModel, dims::AbstractSingleDomain) @@ -17,12 +16,13 @@ Computes the intensity map or _image_ of the `model`. This returns an `Intensity is a `IntensityMap` with `dims` an [`AbstractSingleDomain`](@ref) as dimensions. """ @inline function intensitymap(s::M, - dims::AbstractDomain - ) where {M<:AbstractModel} + dims::AbstractDomain) where {M<:AbstractModel} return create_imgmap(intensitymap(imanalytic(M), s, dims), dims) end -@inline intensitymap(::IsAnalytic, m::AbstractModel, dims::AbstractDomain) = intensitymap_analytic(m, dims) -@inline intensitymap(::NotAnalytic, m::AbstractModel, dims::AbstractDomain) = intensitymap_numeric(m, dims) +@inline intensitymap(::IsAnalytic, m::AbstractModel, dims::AbstractDomain) = intensitymap_analytic(m, + dims) +@inline intensitymap(::NotAnalytic, m::AbstractModel, dims::AbstractDomain) = intensitymap_numeric(m, + dims) function intensitymap_analytic(s::AbstractModel, dims::AbstractDomain) img = allocate_imgmap(s, dims) @@ -30,8 +30,6 @@ function intensitymap_analytic(s::AbstractModel, dims::AbstractDomain) return img end - - """ intensitymap!(img::AbstractIntensityMap, mode;, executor = SequentialEx()) @@ -44,5 +42,5 @@ done. By default we use the `SequentialEx` which uses a single-core to construct @inline function intensitymap!(img, s::M) where {M} return intensitymap!(imanalytic(M), img, s) end -@inline intensitymap!(::IsAnalytic, img, m::AbstractModel) = intensitymap_analytic!(img, m) +@inline intensitymap!(::IsAnalytic, img, m::AbstractModel) = intensitymap_analytic!(img, m) @inline intensitymap!(::NotAnalytic, img, m::AbstractModel) = intensitymap_numeric!(img, m) diff --git a/src/images/methods.jl b/src/images/methods.jl index aedc704..770b43d 100644 --- a/src/images/methods.jl +++ b/src/images/methods.jl @@ -7,19 +7,20 @@ This is useful for broadcasting a model across an abritrary grid. """ domainpoints(img::IntensityMap) = domainpoints(axisdims(img)) - -struct LazySlice{T, N, A<:AbstractVector{T}} <: AbstractArray{T, N} +struct LazySlice{T,N,A<:AbstractVector{T}} <: AbstractArray{T,N} slice::A dir::Int dims::Dims{N} - @inline function LazySlice(slice::AbstractVector{T}, dim::Int, dims::Dims{N}) where {T, N} + @inline function LazySlice(slice::AbstractVector{T}, dim::Int, + dims::Dims{N}) where {T,N} @assert 1 ≤ dim ≤ N "Slice dimension is not valid. Must be ≤ $N and ≥ 1 and got $dim" - return new{T, N, typeof(slice)}(slice, dim, dims) + return new{T,N,typeof(slice)}(slice, dim, dims) end end @inline Base.size(A::LazySlice) = A.dims -Base.@propagate_inbounds @inline function Base.getindex(A::LazySlice{T, N}, I::Vararg{Int, N}) where {T, N} +Base.@propagate_inbounds @inline function Base.getindex(A::LazySlice{T,N}, + I::Vararg{Int,N}) where {T,N} i = I[A.dir] @boundscheck checkbounds(A.slice, i) return A.slice[i] @@ -37,9 +38,6 @@ end @inline basedim(x::DD.LookupArrays.LookupArray) = basedim(parent(x)) @inline basedim(x) = x - - - """ phasecenter(img::IntensityMap) @@ -47,15 +45,13 @@ Computes the phase center of an intensity map. Note this is the pixels that is i the middle of the image. """ function phasecenter(dims::AbstractRectiGrid) - (;X, Y) = dims - x0 = -(last(X) + first(X))/2 - y0 = -(last(Y) + first(Y))/2 + (; X, Y) = dims + x0 = -(last(X) + first(X)) / 2 + y0 = -(last(Y) + first(Y)) / 2 return (X=x0, Y=y0) end phasecenter(img::IntensityMap) = phasecenter(axisdims(img)) - - # ChainRulesCore.@non_differentiable pixelsizes(img::IntensityMap) """ @@ -78,19 +74,19 @@ to the edge of the last pixel. The `x0`, `y0` offsets shift the image origin ove - `executor=Serial()`: The executor to use for the grid, default is serial execution - `header=NoHeader()`: The header to use for the grid """ -function imagepixels(fovx::Real, fovy::Real, nx::Integer, ny::Integer, x0::Real = 0, y0::Real = 0; executor=Serial(), header=NoHeader()) - @assert (nx > 0)&&(ny > 0) "Number of pixels must be positive" +function imagepixels(fovx::Real, fovy::Real, nx::Integer, ny::Integer, x0::Real=0, + y0::Real=0; executor=Serial(), header=NoHeader()) + @assert (nx > 0) && (ny > 0) "Number of pixels must be positive" - psizex=fovx/nx - psizey=fovy/ny + psizex = fovx / nx + psizey = fovy / ny - xitr = X(LinRange(-fovx/2 + psizex/2 - x0, fovx/2 - psizex/2 - x0, nx)) - yitr = Y(LinRange(-fovy/2 + psizey/2 - y0, fovy/2 - psizey/2 - y0, ny)) + xitr = X(LinRange(-fovx / 2 + psizex / 2 - x0, fovx / 2 - psizex / 2 - x0, nx)) + yitr = Y(LinRange(-fovy / 2 + psizey / 2 - y0, fovy / 2 - psizey / 2 - y0, ny)) grid = RectiGrid((xitr, yitr); executor, header) return grid end - """ fieldofview(img::IntensityMap) fieldofview(img::IntensityMap) @@ -103,24 +99,17 @@ end pixelsizes(img::IntensityMap) = pixelsizes(axisdims(img)) - - """ flux(im::IntensityMap) Computes the flux of a intensity map """ function flux(im::IntensityMap{T,N}) where {T,N} - return sum(im, dims=(:X, :Y)) + return sum(im; dims=(:X, :Y)) end flux(im::SpatialIntensityMap) = sum(parent(im)) - - - - - """ centroid(im::AbstractIntensityMap) @@ -129,19 +118,19 @@ Computes the image centroid aka the center of light of the image. For polarized maps we return the centroid for Stokes I only. """ function centroid(im::IntensityMap{<:Number}) - (;X, Y) = named_dims(im) - return mapslices(x->centroid(IntensityMap(x, RectiGrid((;X, Y)))), im; dims=(:X, :Y)) + (; X, Y) = named_dims(im) + return mapslices(x -> centroid(IntensityMap(x, RectiGrid((; X, Y)))), im; dims=(:X, :Y)) end centroid(im::IntensityMap{<:StokesParams}) = centroid(stokes(im, :I)) function centroid(im::IntensityMap{T,2})::Tuple{T,T} where {T<:Real} f = flux(im) cent = sum(pairs(DimPoints(im)); init=SVector(zero(f), zero(f))) do (I, (x, y)) - x0 = x.*im[I] - y0 = y.*im[I] + x0 = x .* im[I] + y0 = y .* im[I] return SVector(x0, y0) end - return cent[1]./f, cent[2]./f + return cent[1] ./ f, cent[2] ./ f end # function ChainRulesCore.rrule(::typeof(centroid), img::IntensityMap{T,2}) where {T<:Real} @@ -157,7 +146,6 @@ end # return out, _centroid_pullback # end - """ second_moment(im::AbstractIntensityMap; center=true) @@ -168,13 +156,15 @@ second moment, which is specified by the `center` argument. For polarized maps we return the second moment for Stokes I only. """ function second_moment(im::IntensityMap{T,N}; center=true) where {T<:Real,N} - (;X, Y) = named_dims(im) - return mapslices(x->second_moment(IntensityMap(x, RectiGrid((;X, Y))); center), im; dims=(:X, :Y)) + (; X, Y) = named_dims(im) + return mapslices(x -> second_moment(IntensityMap(x, RectiGrid((; X, Y))); center), im; + dims=(:X, :Y)) end # Only return the second moment for Stokes I -second_moment(im::IntensityMap{<:StokesParams}; center=true) = second_moment(stokes(im, :I); center) - +function second_moment(im::IntensityMap{<:StokesParams}; center=true) + return second_moment(stokes(im, :I); center) +end """ second_moment(im::IntensityMap; center=true) @@ -188,21 +178,21 @@ function second_moment(im::IntensityMap{T,2}; center=true) where {T<:Real} xy = zero(T) yy = zero(T) f = flux(im) - for (I, (x,y)) in pairs(DimPoints(im)) - xx += x.^2*im[I] - yy += y.^2*im[I] - xy += x.*y.*im[I] + for (I, (x, y)) in pairs(DimPoints(im)) + xx += x .^ 2 * im[I] + yy += y .^ 2 * im[I] + xy += x .* y .* im[I] end if center x0, y0 = centroid(im) - xx = xx./f - x0.^2 - yy = yy./f - y0.^2 - xy = xy./f - x0.*y0 + xx = xx ./ f - x0 .^ 2 + yy = yy ./ f - y0 .^ 2 + xy = xy ./ f - x0 .* y0 else - xx = xx./f - yy = yy./f - xy = xy./f + xx = xx ./ f + yy = yy ./ f + xy = xy ./ f end return @SMatrix [xx xy; xy yy] diff --git a/src/interface.jl b/src/interface.jl index c56d3ff..1069a01 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -39,10 +39,10 @@ export stokes, IsPolarized, NotPolarized struct IsPolarized end struct NotPolarized end -Base.@constprop :aggressive Base.:*(::IsPolarized, ::IsPolarized) = IsPolarized() -Base.@constprop :aggressive Base.:*(::IsPolarized, ::NotPolarized) = IsPolarized() -Base.@constprop :aggressive Base.:*(::NotPolarized, ::IsPolarized) = IsPolarized() -Base.@constprop :aggressive Base.:*(::NotPolarized, ::NotPolarized) = NotPolarized() +Base.@constprop :aggressive Base.:*(::IsPolarized, ::IsPolarized) = IsPolarized() +Base.@constprop :aggressive Base.:*(::IsPolarized, ::NotPolarized) = IsPolarized() +Base.@constprop :aggressive Base.:*(::NotPolarized, ::IsPolarized) = IsPolarized() +Base.@constprop :aggressive Base.:*(::NotPolarized, ::NotPolarized) = NotPolarized() """ ispolarized(::Type) @@ -51,7 +51,6 @@ Trait function that defines whether a model is polarized or not. """ ispolarized(::Type{<:AbstractModel}) = NotPolarized() - """ $(TYPEDEF) @@ -71,7 +70,6 @@ Extract the specific stokes component `p` from the polarized model `m` """ stokes(m::AbstractPolarizedModel, v::Symbol) = getfield(m, v) - """ DensityAnalytic @@ -123,18 +121,15 @@ If `IsAnalytic()` then it will try to call `intensity_point` to calculate the in """ @inline imanalytic(::Type{<:AbstractModel}) = IsAnalytic() - - #= This is internal function definitions for how to compose whether a model is analytic. We need this for composite models. =# -Base.@constprop :aggressive Base.:*(::IsAnalytic, ::IsAnalytic) = IsAnalytic() -Base.@constprop :aggressive Base.:*(::IsAnalytic, ::NotAnalytic) = NotAnalytic() -Base.@constprop :aggressive Base.:*(::NotAnalytic, ::IsAnalytic) = NotAnalytic() -Base.@constprop :aggressive Base.:*(::NotAnalytic, ::NotAnalytic) = NotAnalytic() - +Base.@constprop :aggressive Base.:*(::IsAnalytic, ::IsAnalytic) = IsAnalytic() +Base.@constprop :aggressive Base.:*(::IsAnalytic, ::NotAnalytic) = NotAnalytic() +Base.@constprop :aggressive Base.:*(::NotAnalytic, ::IsAnalytic) = NotAnalytic() +Base.@constprop :aggressive Base.:*(::NotAnalytic, ::NotAnalytic) = NotAnalytic() # Traits are not differentiable # ChainRulesCore.@non_differentiable visanalytic(::Type) @@ -159,7 +154,6 @@ space and invert it. """ function intensity_point end - """ intensitymap!(buffer::AbstractDimArray, model::AbstractModel) @@ -167,7 +161,6 @@ Computes the intensity map of `model` by modifying the `buffer` """ function intensitymap! end - """ intensitymap(model::AbstractModel, p::AbstractSingleDomain) @@ -208,12 +201,6 @@ Updates the `img` using the model `m` by broadcasting over the analytic [`inten """ function intensitymap_analytic! end - - - - - - """ radialextent(model::AbstractModel) @@ -229,7 +216,6 @@ Computes the complex visibilities at the locations p. """ function visibilitymap end - """ visibilitymap!(vis::AbstractArray, model::AbstractModel, p) @@ -237,7 +223,6 @@ Computes the complex visibilities `vis` in place at the locations p """ function visibilitymap! end - """ _visibilitymap(model::AbstractModel, p) @@ -290,7 +275,6 @@ Computes the visibilties of a `model` in-place, using using the analytic visibil """ function visibilitymap_analytic! end - """ stokes(m::AbstractArray{<:StokesParams}, p::Symbol) diff --git a/src/unstructured_map.jl b/src/unstructured_map.jl index 8cc2937..96199ec 100644 --- a/src/unstructured_map.jl +++ b/src/unstructured_map.jl @@ -1,6 +1,5 @@ export UnstructuredMap - """ UnstructuredMap(data::AbstractVector, dims::UnstructuredDomain) @@ -12,7 +11,7 @@ For instance the locations of the visibilities can be accessed with `axisdims`, as well as the usual `getproperty` and `propertynames` functions. Like with `IntensityMap` during execution the `executor` is used to determine the execution context. """ -struct UnstructuredMap{T, A<:AbstractVector{T}, G<:UnstructuredDomain} <: AbstractVector{T} +struct UnstructuredMap{T,A<:AbstractVector{T},G<:UnstructuredDomain} <: AbstractVector{T} data::A dims::G end @@ -23,23 +22,26 @@ executor(x::UnstructuredMap) = executor(axisdims(x)) baseimage(x::UnstructuredMap) = baseimage(parent(x)) Base.parent(x::UnstructuredMap) = getfield(x, :data) Base.length(x::UnstructuredMap) = length(parent(x)) -Base.IndexStyle(::Type{<:UnstructuredMap{T, A}}) where {T, A} = IndexStyle(A) +Base.IndexStyle(::Type{<:UnstructuredMap{T,A}}) where {T,A} = IndexStyle(A) Base.@propagate_inbounds Base.getindex(a::UnstructuredMap, i::Int) = getindex(parent(a), i) -Base.@propagate_inbounds Base.setindex!(a::UnstructuredMap, v, i::Int) = setindex!(parent(a), v, i) +Base.@propagate_inbounds Base.setindex!(a::UnstructuredMap, v, i::Int) = setindex!(parent(a), + v, i) -UnstructuredMap(data::UnstructuredMap, dims::UnstructuredDomain) = UnstructuredMap(parent(data), dims) +function UnstructuredMap(data::UnstructuredMap, dims::UnstructuredDomain) + return UnstructuredMap(parent(data), dims) +end function stokes(x::UnstructuredMap, p::Symbol) return UnstructuredMap(stokes(parent(x), p), axisdims(x)) end - function Base.similar(m::UnstructuredMap, ::Type{S}) where {S} return UnstructuredMap(similar(parent(m), S), axisdims(m)) end Base.BroadcastStyle(::Type{<:UnstructuredMap}) = Broadcast.ArrayStyle{UnstructuredMap}() -function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{UnstructuredMap}}, ::Type{ElType}) where {ElType} +function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{UnstructuredMap}}, + ::Type{ElType}) where {ElType} # Scan inputs for the time and sites sarr = find_ustr(bc) return UnstructuredMap(similar(parent(sarr), ElType), axisdims(sarr)) @@ -61,26 +63,26 @@ function Base.view(x::UnstructuredMap, I) dims = axisdims(x) g = domainpoints(dims) newdims = UnstructuredDomain(@view(g[I]), executor(dims), header(dims)) - UnstructuredMap(view(parent(x), I), newdims) + return UnstructuredMap(view(parent(x), I), newdims) end Base.@propagate_inbounds function Base.getindex(x::UnstructuredMap, I) dims = axisdims(x) g = domainpoints(dims) newdims = UnstructuredDomain((g[I]), executor(dims), header(dims)) - UnstructuredMap(parent(x)[I], newdims) + return UnstructuredMap(parent(x)[I], newdims) end -function intensitymap_analytic!(img::UnstructuredMap{T, <:Any, <:AbstractSingleDomain}, s::AbstractModel) where {T} +function intensitymap_analytic!(img::UnstructuredMap{T,<:Any,<:AbstractSingleDomain}, + s::AbstractModel) where {T} g = domainpoints(img) img .= intensity_point.(Ref(s), g) return nothing end - -function intensitymap_analytic!( - img::UnstructuredMap{T,<:Any,<:UnstructuredDomain{D, <:ThreadsEx{S}}}, - s::AbstractModel) where {T,D,S} +function intensitymap_analytic!(img::UnstructuredMap{T,<:Any, + <:UnstructuredDomain{D,<:ThreadsEx{S}}}, + s::AbstractModel) where {T,D,S} g = domainpoints(img) _threads_intensitymap!(img, s, g, Val(S)) return nothing @@ -88,7 +90,8 @@ end for s in schedulers @eval begin - function _threads_intensitymap!(img::UnstructuredMap, s::AbstractModel, g, ::Val{$s}) + function _threads_intensitymap!(img::UnstructuredMap, s::AbstractModel, g, + ::Val{$s}) Threads.@threads $s for I in CartesianIndices(g) img[I] = intensity_point(s, g[I]) end diff --git a/src/visibilities.jl b/src/visibilities.jl index c339282..b8b8920 100644 --- a/src/visibilities.jl +++ b/src/visibilities.jl @@ -1,9 +1,8 @@ export visibilitymap, visibilitymap!, - logclosure_amplitude, logclosure_amplitudemap, - amplitude, amplitudemap, - closure_phase, closure_phasemap, - bispectrum, bispectrummap - + logclosure_amplitude, logclosure_amplitudemap, + amplitude, amplitudemap, + closure_phase, closure_phasemap, + bispectrum, bispectrummap """ visibilitymap(m, p) @@ -14,10 +13,8 @@ are expected to have the properties `U`, `V`, and sometimes `T` and `F`. @inline function visibilitymap(m::M, p::AbstractDomain) where {M<:AbstractModel} return create_vismap(_visibilitymap(visanalytic(M), m, p), p) end -@inline _visibilitymap(::IsAnalytic, m::AbstractModel, p) = visibilitymap_analytic(m, p) -@inline _visibilitymap(::NotAnalytic, m::AbstractModel, p) = visibilitymap_numeric(m, p) - - +@inline _visibilitymap(::IsAnalytic, m::AbstractModel, p) = visibilitymap_analytic(m, p) +@inline _visibilitymap(::NotAnalytic, m::AbstractModel, p) = visibilitymap_numeric(m, p) """ visibilitymap!(vis, m, p) @@ -28,8 +25,10 @@ are expected to have the properties `U`, `V`, and sometimes `T` and `F`. @inline function visibilitymap!(vis, m::M) where {M<:AbstractModel} return _visibilitymap!(visanalytic(M), vis, m) end -@inline _visibilitymap!(::IsAnalytic , vis, m::AbstractModel) = visibilitymap_analytic!(vis, m) -@inline _visibilitymap!(::NotAnalytic, vis, m::AbstractModel) = visibilitymap_numeric!(vis, m) +@inline _visibilitymap!(::IsAnalytic, vis, m::AbstractModel) = visibilitymap_analytic!(vis, + m) +@inline _visibilitymap!(::NotAnalytic, vis, m::AbstractModel) = visibilitymap_numeric!(vis, + m) function visibilitymap_analytic(m::AbstractModel, p::AbstractSingleDomain) vis = allocate_vismap(m, p) @@ -45,14 +44,16 @@ function visibilitymap_analytic!(vis, m::AbstractModel) return nothing end -function visibilitymap_analytic!(vis::FluxMap2{T, <:Any, <:ComradeBase.AbstractSingleDomain{<:Any, <:ThreadsEx{S}}}, m::AbstractModel) where {T,S} +function visibilitymap_analytic!(vis::FluxMap2{T,<:Any, + <:ComradeBase.AbstractSingleDomain{<:Any, + <:ThreadsEx{S}}}, + m::AbstractModel) where {T,S} d = axisdims(vis) g = domainpoints(d) _threads_visibilitymap!(parent(vis), m, g, Val(S)) return nothing end - for s in schedulers @eval begin function _threads_visibilitymap!(vis, s::AbstractModel, g, ::Val{$s}) @@ -64,9 +65,6 @@ for s in schedulers end end - - - """ visibility(mimg, p) @@ -92,9 +90,6 @@ end return visibility_point(mimg, p) end - - - """ amplitude(model, p) @@ -119,7 +114,7 @@ If you want to compute the bispectrum over a number of triangles consider using the `bispectrummap` function. """ @inline function bispectrum(model, p1, p2, p3) - return visibility(model, p1)*visibility(model, p2)*visibility(model, p3) + return visibility(model, p1) * visibility(model, p2) * visibility(model, p3) end """ @@ -154,13 +149,9 @@ consider using the `logclosure_amplitudemap` function. a3 = amplitude(model, p3) a4 = amplitude(model, p4) - return log(a1*a2/(a3*a4)) + return log(a1 * a2 / (a3 * a4)) end - - - - """ amplitudemap(m::AbstractModel, p) @@ -169,24 +160,22 @@ The coordinates `p` are expected to have the properties `U`, `V`, and sometimes `Ti` and `Fr`. """ function amplitudemap(m, p) - create_map(_amplitudemap(m, p), p) + return create_map(_amplitudemap(m, p), p) end - function _amplitudemap(m::S, p) where {S} - _amplitudemap(visanalytic(S), m, p) + return _amplitudemap(visanalytic(S), m, p) end function _amplitudemap(::IsAnalytic, m, p::AbstractSingleDomain) g = domainpoints(p) - abs.(visibility_point.(Ref(m), g)) + return abs.(visibility_point.(Ref(m), g)) end function _amplitudemap(::NotAnalytic, m, p::AbstractSingleDomain) - abs.(visibilitymap_numeric(m, p)) + return abs.(visibilitymap_numeric(m, p)) end - """ bispectrummap(m, p1, p2, p3) @@ -194,29 +183,25 @@ Computes the closure phases of the model `m` at the triangles p1, p2, p3, where `pi` are coordinates. """ function bispectrummap(m, - p1::T, - p2::T, - p3::T, - ) where {T<:AbstractSingleDomain} - - _bispectrummap(m, p1, p2, p3) + p1::T, + p2::T, + p3::T) where {T<:AbstractSingleDomain} + return _bispectrummap(m, p1, p2, p3) end # internal method used for trait dispatch function _bispectrummap(m::M, - p1, - p2, - p3 - ) where {M} - _bispectrummap(visanalytic(M), m, p1, p2, p3) + p1, + p2, + p3) where {M} + return _bispectrummap(visanalytic(M), m, p1, p2, p3) end # internal method used for trait dispatch for analytic visibilities function _bispectrummap(::IsAnalytic, m, - p1::T, - p2::T, - p3::T, - ) where {T<:AbstractSingleDomain} + p1::T, + p2::T, + p3::T) where {T<:AbstractSingleDomain} g1 = domainpoints(p1) g2 = domainpoints(p2) g3 = domainpoints(p3) @@ -225,12 +210,11 @@ end # internal method used for trait dispatch for non-analytic visibilities function _bispectrummap(::NotAnalytic, m, - p1::T,p2::T,p3::T - ) where {T<:AbstractSingleDomain} + p1::T, p2::T, p3::T) where {T<:AbstractSingleDomain} vis1 = visibilitymap(m, p1) vis2 = visibilitymap(m, p2) vis3 = visibilitymap(m, p3) - return @. vis1*vis2*vis3 + return @. vis1 * vis2 * vis3 end """ @@ -244,22 +228,22 @@ Computes the closure phases of the model `m` at the triangles p1, p2, p3, where `pi` are coordinates. """ @inline function closure_phasemap(m::AbstractModel, - p1::T,p2::T,p3::T - ) where {T<:UnstructuredDomain} - create_map(_closure_phasemap(m, p1, p2, p3), UnstructuredDomain((;p1=domainpoints(p1), p2=domainpoints(p2), p3=domainpoints(p3)))) + p1::T, p2::T, p3::T) where {T<:UnstructuredDomain} + return create_map(_closure_phasemap(m, p1, p2, p3), + UnstructuredDomain((; p1=domainpoints(p1), p2=domainpoints(p2), + p3=domainpoints(p3)))) end # internal method used for trait dispatch @inline function _closure_phasemap(m::M, p1, p2, p3) where {M<:AbstractModel} - _closure_phasemap(visanalytic(M), m, p1, p2, p3) + return _closure_phasemap(visanalytic(M), m, p1, p2, p3) end # internal method used for trait dispatch for analytic visibilities @inline function _closure_phasemap(::IsAnalytic, m, - p1::UnstructuredDomain, - p2::UnstructuredDomain, - p3::UnstructuredDomain - ) + p1::UnstructuredDomain, + p2::UnstructuredDomain, + p3::UnstructuredDomain) g1 = domainpoints(p1) g2 = domainpoints(p2) g3 = domainpoints(p3) @@ -267,7 +251,7 @@ end end # internal method used for trait dispatch for non-analytic visibilities -function _closure_phasemap(::NotAnalytic, m, p1,p2, p3) +function _closure_phasemap(::NotAnalytic, m, p1, p2, p3) return angle.(bispectrummap(m, p1, p2, p3)) end @@ -283,24 +267,22 @@ Computes the log closure amplitudemap of the model `m` at the quadrangles p1, p2, p3, p4. """ function logclosure_amplitudemap(m::AbstractModel, - p1::T, - p2::T, - p3::T, - p4::T - ) where {T<:AbstractSingleDomain} - glc = UnstructuredDomain((;p1 = domainpoints(p1), p2 = domainpoints(p2), p3 = domainpoints(p3), p4 = domainpoints(p4))) - create_map(_logclosure_amplitudemap(m, p1, p2, p3, p4), glc) + p1::T, + p2::T, + p3::T, + p4::T) where {T<:AbstractSingleDomain} + glc = UnstructuredDomain((; p1=domainpoints(p1), p2=domainpoints(p2), + p3=domainpoints(p3), p4=domainpoints(p4))) + return create_map(_logclosure_amplitudemap(m, p1, p2, p3, p4), glc) end - # internal method used for trait dispatch @inline function _logclosure_amplitudemap(m::M, - p1, - p2, - p3, - p4 - ) where {M<:AbstractModel} - _logclosure_amplitudemap(visanalytic(M), m, p1, p2, p3, p4) + p1, + p2, + p3, + p4) where {M<:AbstractModel} + return _logclosure_amplitudemap(visanalytic(M), m, p1, p2, p3, p4) end # internal method used for trait dispatch for analytic visibilities @@ -318,5 +300,5 @@ end amp2 = _amplitudemap(m, p2) amp3 = _amplitudemap(m, p3) amp4 = _amplitudemap(m, p4) - return @. log(amp1*amp2*inv(amp3*amp4)) + return @. log(amp1 * amp2 * inv(amp3 * amp4)) end diff --git a/test/executors.jl b/test/executors.jl index 1e5a997..298e4cc 100644 --- a/test/executors.jl +++ b/test/executors.jl @@ -19,10 +19,9 @@ function testexvis(img, m, ex) @test img ≈ img2 end - @testset "executors" begin - u = 0.1*randn(60) - v = 0.1*randn(60) + u = 0.1 * randn(60) + v = 0.1 * randn(60) ti = collect(Float64, 1:60) fr = fill(230e9, 60) m = GaussTest() @@ -30,7 +29,7 @@ end @test ThreadsEx() === ThreadsEx(:dynamic) @testset "RectiGrid" begin - pim = (;X=range(-10.0, 10.0, length=64), Y=range(-10.0, 10.0, length=64)) + pim = (; X=range(-10.0, 10.0; length=64), Y=range(-10.0, 10.0; length=64)) gim = RectiGrid(pim) img = intensitymap(m, gim) img0 = copy(img) @@ -43,7 +42,7 @@ end testeximg(img, m, StaticScheduler()) testeximg(img, m, SerialScheduler()) - puv = (U=range(-2.0, 2.0, length=128), V=range(-2.0, 2.0, length=64)) + puv = (U=range(-2.0, 2.0; length=128), V=range(-2.0, 2.0; length=64)) vis = visibilitymap(m, RectiGrid(puv)) vis0 = copy(vis) visibilitymap!(vis, m) @@ -58,15 +57,14 @@ end end @testset "UnstructuredDomain" begin - pim = (;X=randn(64), Y=randn(64)) - puv = (;U=u, V=v, Ti=ti, Fr=fr) + pim = (; X=randn(64), Y=randn(64)) + puv = (; U=u, V=v, Ti=ti, Fr=fr) gim = UnstructuredDomain(pim) img = intensitymap(m, gim) img0 = copy(img) intensitymap!(img, m) @test img ≈ img0 - testeximg(img, m, ThreadsEx()) testeximg(img, m, ThreadsEx(:static)) testeximg(img, m, DynamicScheduler()) @@ -84,13 +82,11 @@ end testexvis(vis, m, StaticScheduler()) testexvis(vis, m, SerialScheduler()) end - end - @testset "executors NotAnalytic" begin - u = 0.1*randn(60) - v = 0.1*randn(60) + u = 0.1 * randn(60) + v = 0.1 * randn(60) ti = collect(Float64, 1:60) fr = fill(230e9, 60) m = GaussTestNA() @@ -98,7 +94,7 @@ end @test ThreadsEx() === ThreadsEx(:dynamic) @testset "RectiGrid" begin - pim = (;X=range(-10.0, 10.0, length=64), Y=range(-10.0, 10.0, length=64)) + pim = (; X=range(-10.0, 10.0; length=64), Y=range(-10.0, 10.0; length=64)) gim = RectiGrid(pim) img = intensitymap(m, gim) @@ -108,7 +104,7 @@ end @test img ≈ intensitymap(m, RectiGrid(pim; executor=StaticScheduler())) @test img ≈ intensitymap(m, RectiGrid(pim; executor=SerialScheduler())) - puv = (U=range(-2.0, 2.0, length=128), V=range(-2.0, 2.0, length=64)) + puv = (U=range(-2.0, 2.0; length=128), V=range(-2.0, 2.0; length=64)) vis = visibilitymap(m, RectiGrid(puv)) @test size(vis) == size(RectiGrid(puv)) @test vis ≈ visibilitymap(m, RectiGrid(puv; executor=ThreadsEx())) @@ -119,8 +115,8 @@ end end @testset "UnstructuredDomain" begin - pim = (;X=randn(64), Y=randn(64)) - puv = (;U=u, V=v, Ti=ti, Fr=fr) + pim = (; X=randn(64), Y=randn(64)) + puv = (; U=u, V=v, Ti=ti, Fr=fr) gim = UnstructuredDomain(pim) img = intensitymap(m, gim) @@ -138,20 +134,19 @@ end @test vis ≈ visibilitymap(m, UnstructuredDomain(puv; executor=StaticScheduler())) @test vis ≈ visibilitymap(m, UnstructuredDomain(puv; executor=SerialScheduler())) end - end @testset "EnzymeExecutors" begin - u = 0.1*randn(60) - v = 0.1*randn(60) + u = 0.1 * randn(60) + v = 0.1 * randn(60) ti = collect(Float64, 1:60) fr = fill(230e9, 60) m = GaussTest() - pim = (;X=range(-10.0, 10.0, length=64), Y=range(-10.0, 10.0, length=64)) + pim = (; X=range(-10.0, 10.0; length=64), Y=range(-10.0, 10.0; length=64)) gim = RectiGrid(pim) - guv = UnstructuredDomain((;U=u, V=v, Ti=ti, Fr=fr)) - guvm = RectiGrid((;U=u, V=v)) + guv = UnstructuredDomain((; U=u, V=v, Ti=ti, Fr=fr)) + guvm = RectiGrid((; U=u, V=v)) img = intensitymap(m, gim) vis = visibilitymap(m, guv) diff --git a/test/images.jl b/test/images.jl index 9275918..81b1015 100644 --- a/test/images.jl +++ b/test/images.jl @@ -1,4 +1,4 @@ -function test_grid_interface(grid::ComradeBase.AbstractSingleDomain{D, E}) where {D,E} +function test_grid_interface(grid::ComradeBase.AbstractSingleDomain{D,E}) where {D,E} @test typeof(executor(grid)) == E arr = zeros(size(grid)) @inferred ComradeBase.create_map(arr, grid) @@ -20,19 +20,19 @@ function test_grid_interface(grid::ComradeBase.AbstractSingleDomain{D, E}) where axes(grid) show(grid) show(IOBuffer(), MIME"text/plain"(), grid) - summary(grid) + return summary(grid) end @testset "AbstractSingleDomain" begin ex = Serial() - prect = (;X=range(-10.0, 10.0, length=128), - Y=range(-10.0, 10.0, length=128), - Fr = [230.0, 345.0], - Ti = sort(rand(24))) - pustr = (;X=range(-10.0, 10.0, length=128), - Y=range(-10.0, 10.0, length=128), - Fr = fill(230e9, 128), - Ti = sort(rand(128))) + prect = (; X=range(-10.0, 10.0; length=128), + Y=range(-10.0, 10.0; length=128), + Fr=[230.0, 345.0], + Ti=sort(rand(24))) + pustr = (; X=range(-10.0, 10.0; length=128), + Y=range(-10.0, 10.0; length=128), + Fr=fill(230e9, 128), + Ti=sort(rand(128))) grect = RectiGrid(prect) pc = phasecenter(grect) @test pc.X ≈ 0.0 @@ -45,23 +45,21 @@ end head = ComradeBase.MinimalHeader("M87", 90.0, 45, 21312, 230e9) g = RectiGrid(prect; header=head) @test header(g) == head - end @testset "IntensityMap" begin - x = X(range(-10.0, 10.0, length=64)) - y = Y(range(-10.0, 10.0, length=64)) + x = X(range(-10.0, 10.0; length=64)) + y = Y(range(-10.0, 10.0; length=64)) t = Ti([0.0, 0.5, 0.8]) f = Fr([86e9, 230e9, 345e9]) - gsp = RectiGrid((x,y)) + gsp = RectiGrid((x, y)) g1 = RectiGrid((x, y, f, t)) g2 = RectiGrid((x, y, t, f)) imp = rand(64, 64, 3, 3) - - img1 = IntensityMap(imp[:,:,1,1], gsp) + img1 = IntensityMap(imp[:, :, 1, 1], gsp) img2 = IntensityMap(imp, g1) img3 = IntensityMap(imp, g2) phasecenter(img2) @@ -74,21 +72,16 @@ end @test_throws ArgumentError img1.Fr - @testset "Slicing" begin @test img1[X=1:1, Y=1:10] isa IntensityMap - @test img1[X=5:10, Y=1:end-10] isa IntensityMap + @test img1[X=5:10, Y=1:(end - 10)] isa IntensityMap @test img2[X=1, Y=1] isa IntensityMap - - - - @test img1[X=1, Y=1] ≈ imp[1,1,1,1] - @test img1[X=1, Y=1:10] ≈ imp[1,1:10,1,1] - @test img1[X=5:10, Y=1:end-10] ≈ imp[5:10,1:end-10,1,1] - @test img1[Y=1, X=1] ≈ imp[1,1,1,1] - @test img2[X=1, Y=1] ≈ imp[1,1,:,:] - + @test img1[X=1, Y=1] ≈ imp[1, 1, 1, 1] + @test img1[X=1, Y=1:10] ≈ imp[1, 1:10, 1, 1] + @test img1[X=5:10, Y=1:(end - 10)] ≈ imp[5:10, 1:(end - 10), 1, 1] + @test img1[Y=1, X=1] ≈ imp[1, 1, 1, 1] + @test img2[X=1, Y=1] ≈ imp[1, 1, :, :] subimg1 = img1[X=5:10, Y=1:20] nk = named_dims(subimg1) @@ -104,10 +97,10 @@ end end @testset "broadcast and map" begin - @test img1.^2 isa typeof(img1) + @test img1 .^ 2 isa typeof(img1) @test cos.(img1) isa typeof(img1) @test img1 .+ img1 isa typeof(img1) - @test cos.(img2[Fr=1,Ti=1]) isa IntensityMap + @test cos.(img2[Fr=1, Ti=1]) isa IntensityMap end @testset "polarized" begin @@ -116,17 +109,16 @@ end imgU = rand(64, 64, 3, 3) imgV = rand(64, 64, 3, 3) - imgP = StructArray{StokesParams}(I=imgI, Q=imgQ, U=imgU, V=imgV) - img1 = IntensityMap(imgP[:,:,1,1], RectiGrid((;X=x,Y=y))) + imgP = StructArray{StokesParams}(; I=imgI, Q=imgQ, U=imgU, V=imgV) + img1 = IntensityMap(imgP[:, :, 1, 1], RectiGrid((; X=x, Y=y))) img2 = IntensityMap(imgP, RectiGrid((x, y, t, f))) - @test flux(img1) ≈ flux(img2)[1,1,1,1] + @test flux(img1) ≈ flux(img2)[1, 1, 1, 1] @test centroid(img1) == centroid(stokes(img1, :I)) @test second_moment(img1) == second_moment(stokes(img1, :I)) end end - function FiniteDifferences.to_vec(k::IntensityMap) v, b = to_vec(DD.data(k)) back(x) = DD.rebuild(k, b(x)) @@ -140,7 +132,6 @@ function FiniteDifferences.to_vec(k::UnstructuredMap) return v, back end - # @testset "ProjectTo" begin # data = rand(32, 32) @@ -149,7 +140,6 @@ end # # test_rrule(centroid, img) - # # pr = ProjectTo(img) # # @test pr(data) == img # # @test pr(NoTangent()) == NoTangent() @@ -172,7 +162,6 @@ end # test_rrule(UnstructuredMap, data, g⊢NoTangent()) # end - # @testset "rrule baseimage" begin # data = rand(32, 24) # g = imagepixels(5.0, 10.0, 32, 24) @@ -182,16 +171,16 @@ end # end @testset "UnstructuredMap" begin - pustr = (;X=range(-10.0, 10.0, length=128), - Y=range(-10.0, 10.0, length=128), - Fr = fill(230e9, 128), - Ti = sort(rand(128))) + pustr = (; X=range(-10.0, 10.0; length=128), + Y=range(-10.0, 10.0; length=128), + Fr=fill(230e9, 128), + Ti=sort(rand(128))) g = UnstructuredDomain(pustr) img = UnstructuredMap(rand(128), g) - @test typeof(img.^2) == typeof(img) - @test img[[1,4,6]] isa UnstructuredMap - @test view(img, [1,4,6]) isa UnstructuredMap + @test typeof(img .^ 2) == typeof(img) + @test img[[1, 4, 6]] isa UnstructuredMap + @test view(img, [1, 4, 6]) isa UnstructuredMap @test img[[6]][1] == img[6] @test @view(img[[6]])[1] == img[6] diff --git a/test/interface.jl b/test/interface.jl index e22edcc..a422b24 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -1,16 +1,15 @@ using ComradeBase: IsAnalytic, NotAnalytic, IsPolarized, NotPolarized @testset "interface" begin - @test IsAnalytic()*NotAnalytic() == NotAnalytic() - @test IsAnalytic()*IsAnalytic() == IsAnalytic() - @test NotAnalytic()*NotAnalytic() == NotAnalytic() - @test NotAnalytic()*IsAnalytic() == NotAnalytic() + @test IsAnalytic() * NotAnalytic() == NotAnalytic() + @test IsAnalytic() * IsAnalytic() == IsAnalytic() + @test NotAnalytic() * NotAnalytic() == NotAnalytic() + @test NotAnalytic() * IsAnalytic() == NotAnalytic() @test ComradeBase.ispolarized(ComradeBase.AbstractPolarizedModel) == IsPolarized() @test ComradeBase.ispolarized(ComradeBase.AbstractModel) == NotPolarized() - @test IsPolarized()*NotPolarized() == IsPolarized() - @test IsPolarized()*IsPolarized() == IsPolarized() - @test NotPolarized()*NotPolarized() == NotPolarized() - @test NotPolarized()*IsPolarized() == IsPolarized() - + @test IsPolarized() * NotPolarized() == IsPolarized() + @test IsPolarized() * IsPolarized() == IsPolarized() + @test NotPolarized() * NotPolarized() == NotPolarized() + @test NotPolarized() * IsPolarized() == IsPolarized() end diff --git a/test/multidomain.jl b/test/multidomain.jl index 39543aa..54c11eb 100644 --- a/test/multidomain.jl +++ b/test/multidomain.jl @@ -1,7 +1,7 @@ function domain4d(N, Nt, Nf) - U_vals = range(-10e9, 10e9, length=N) + U_vals = range(-10e9, 10e9; length=N) U_vals = U_vals' .* ones(N) - V_vals = range(-10e9, 10e9, length=N) + V_vals = range(-10e9, 10e9; length=N) V_vals = V_vals' .* ones(N) U_vals = U_vals' @@ -9,26 +9,29 @@ function domain4d(N, Nt, Nf) U_final = vec(U_vals) V_final = vec(V_vals) - ti = 10 * rand(Nt) |> sort - fr = 1e11 * rand(Nf) |> sort + ti = sort(10 * rand(Nt)) + fr = sort(1e11 * rand(Nf)) # Repeat U and V to match Ti dimensions - U_repeated = repeat(U_final, outer=(length(ti))) - V_repeated = repeat(V_final, outer=(length(ti))) - Ti_repeated = repeat(ti, inner=(Int(length(U_final)))) + U_repeated = repeat(U_final; outer=(length(ti))) + V_repeated = repeat(V_final; outer=(length(ti))) + Ti_repeated = repeat(ti; inner=(Int(length(U_final)))) # Repeat U and V and Ti to match Fr dimensions - U_repeated = repeat(U_repeated, outer=(length(fr))) - V_repeated = repeat(V_repeated, outer=(length(fr))) - Ti_repeated = repeat(Ti_repeated, outer=(length(fr))) - Fr_repeated = repeat(fr, inner=(Int(length(U_repeated)/length(fr)))) - visdomain = UnstructuredDomain((;U=U_repeated, V=V_repeated, Ti=Ti_repeated, Fr=Fr_repeated)) + U_repeated = repeat(U_repeated; outer=(length(fr))) + V_repeated = repeat(V_repeated; outer=(length(fr))) + Ti_repeated = repeat(Ti_repeated; outer=(length(fr))) + Fr_repeated = repeat(fr; inner=(Int(length(U_repeated) / length(fr)))) + visdomain = UnstructuredDomain((; U=U_repeated, V=V_repeated, Ti=Ti_repeated, + Fr=Fr_repeated)) C1 = true for ti_point in ti for fr_point in fr f = visdomain[Ti=ti_point, Fr=fr_point] - g = UnstructuredDomain((;U=U_final, V=V_final, Ti=vcat(fill(ti_point, length(U_final))), Fr=vcat(fill(fr_point, length(U_final))))) + g = UnstructuredDomain((; U=U_final, V=V_final, + Ti=vcat(fill(ti_point, length(U_final))), + Fr=vcat(fill(fr_point, length(U_final))))) C1 = C1 && (domainpoints(f) == domainpoints(g)) end end @@ -38,7 +41,9 @@ function domain4d(N, Nt, Nf) for ti_point in ti for fr_point in fr f = visdomain[Fr=fr_point, Ti=ti_point] - g = UnstructuredDomain((;U=U_final, V=V_final, Ti=vcat(fill(ti_point, length(U_final))), Fr=vcat(fill(fr_point, length(U_final))))) + g = UnstructuredDomain((; U=U_final, V=V_final, + Ti=vcat(fill(ti_point, length(U_final))), + Fr=vcat(fill(fr_point, length(U_final))))) C2 = C2 && (domainpoints(f) == domainpoints(g)) end end @@ -47,9 +52,9 @@ function domain4d(N, Nt, Nf) end function domain3df(N, Nf) - U_vals = range(-10e9, 10e9, length=N) + U_vals = range(-10e9, 10e9; length=N) U_vals = U_vals' .* ones(N) - V_vals = range(-10e9, 10e9, length=N) + V_vals = range(-10e9, 10e9; length=N) V_vals = V_vals' .* ones(N) U_vals = U_vals' @@ -57,19 +62,20 @@ function domain3df(N, Nf) U_final = vec(U_vals) V_final = vec(V_vals) - fr = 1e11 * rand(Nf) |> sort + fr = sort(1e11 * rand(Nf)) # Repeat U and V to match Fr dimensions - U_repeated = repeat(U_final, outer=(length(fr))) - V_repeated = repeat(V_final, outer=(length(fr))) - Fr_repeated = repeat(fr, inner=(Int(length(U_final)))) + U_repeated = repeat(U_final; outer=(length(fr))) + V_repeated = repeat(V_final; outer=(length(fr))) + Fr_repeated = repeat(fr; inner=(Int(length(U_final)))) - visdomain = UnstructuredDomain((;U=U_repeated, V=V_repeated, Fr=Fr_repeated)) + visdomain = UnstructuredDomain((; U=U_repeated, V=V_repeated, Fr=Fr_repeated)) C = true for fr_point in fr f = visdomain[Fr=fr_point] - g = UnstructuredDomain((;U=U_final, V=V_final, Fr=vcat(fill(fr_point, length(U_final))))) + g = UnstructuredDomain((; U=U_final, V=V_final, + Fr=vcat(fill(fr_point, length(U_final))))) C = C && (domainpoints(f) == domainpoints(g)) end @@ -77,9 +83,9 @@ function domain3df(N, Nf) end function domain3dt(N, Nt) - U_vals = range(-10e9, 10e9, length=N) + U_vals = range(-10e9, 10e9; length=N) U_vals = U_vals' .* ones(N) - V_vals = range(-10e9, 10e9, length=N) + V_vals = range(-10e9, 10e9; length=N) V_vals = V_vals' .* ones(N) U_vals = U_vals' @@ -87,32 +93,32 @@ function domain3dt(N, Nt) U_final = vec(U_vals) V_final = vec(V_vals) - ti = 10 * rand(Nt) |> sort + ti = sort(10 * rand(Nt)) # Repeat U and V to match Ti dimensions - U_repeated = repeat(U_final, outer=(length(ti))) - V_repeated = repeat(V_final, outer=(length(ti))) - Ti_repeated = repeat(ti, inner=(Int(length(U_final)))) + U_repeated = repeat(U_final; outer=(length(ti))) + V_repeated = repeat(V_final; outer=(length(ti))) + Ti_repeated = repeat(ti; inner=(Int(length(U_final)))) - visdomain = UnstructuredDomain((;U=U_repeated, V=V_repeated, Ti=Ti_repeated)) + visdomain = UnstructuredDomain((; U=U_repeated, V=V_repeated, Ti=Ti_repeated)) C = true for ti_point in ti f = visdomain[Ti=ti_point] - g = UnstructuredDomain((;U=U_final, V=V_final, Ti=vcat(fill(ti_point, length(U_final))))) + g = UnstructuredDomain((; U=U_final, V=V_final, + Ti=vcat(fill(ti_point, length(U_final))))) C = C && (domainpoints(f) == domainpoints(g)) end return C end - @testset "Test getindex for visdomain" begin - C1, C2 = domain4d(64,10,4) + C1, C2 = domain4d(64, 10, 4) @test C1 @test C2 - C3 = domain3dt(64,10) + C3 = domain3dt(64, 10) @test C3 - C4 = domain3df(64,4) + C4 = domain3df(64, 4) @test C4 -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index d17c9a3..1ce0b53 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,6 @@ using FiniteDifferences # using ChainRulesTestUtils import DimensionalData as DD - @testset "ComradeBase.jl" begin include(joinpath(@__DIR__, "interface.jl")) include(joinpath(@__DIR__, "images.jl")) diff --git a/test/visibilities.jl b/test/visibilities.jl index 98013b8..6f9f380 100644 --- a/test/visibilities.jl +++ b/test/visibilities.jl @@ -1,5 +1,4 @@ - struct GaussTest{T} <: ComradeBase.AbstractModel end GaussTest() = GaussTest{Float64}() @@ -8,18 +7,18 @@ ComradeBase.imanalytic(::Type{<:GaussTest}) = ComradeBase.IsAnalytic() ComradeBase.ispolarized(::Type{<:GaussTest}) = ComradeBase.NotPolarized() function ComradeBase.intensity_point(::GaussTest, p) - (;X, Y) = p - return exp(-(X^2+Y^2)/2)/2π + (; X, Y) = p + return exp(-(X^2 + Y^2) / 2) / 2π end function ComradeBase.visibility_point(::GaussTest, p) u = p.U v = p.V - return complex(exp(-2π^2*(u^2 + v^2))) + return complex(exp(-2π^2 * (u^2 + v^2))) end ComradeBase.flux(::GaussTest{T}) where {T} = one(T) -ComradeBase.radialextent(::GaussTest{T}) where {T} = 5*one(T) +ComradeBase.radialextent(::GaussTest{T}) where {T} = 5 * one(T) struct GaussTestNA{T} <: ComradeBase.AbstractModel end GaussTestNA() = GaussTestNA{Float64}() @@ -29,18 +28,19 @@ ComradeBase.imanalytic(::Type{<:GaussTestNA}) = ComradeBase.NotAnalytic() ComradeBase.ispolarized(::Type{<:GaussTestNA}) = ComradeBase.NotPolarized() function ComradeBase.intensity_point(::GaussTestNA, p) - (;X, Y) = p - return exp(-(X^2+Y^2)/2)/2π + (; X, Y) = p + return exp(-(X^2 + Y^2) / 2) / 2π end function ComradeBase.visibility_point(::GaussTestNA, p) u = p.U v = p.V - return complex(exp(-2π^2*(u^2 + v^2))) + return complex(exp(-2π^2 * (u^2 + v^2))) end # Fake it to for testing -function ComradeBase.intensitymap_numeric(m::GaussTestNA, p::ComradeBase.AbstractSingleDomain) +function ComradeBase.intensitymap_numeric(m::GaussTestNA, + p::ComradeBase.AbstractSingleDomain) return ComradeBase.intensitymap_analytic(m, p) end @@ -48,8 +48,8 @@ function ComradeBase.intensitymap_numeric!(img, m::GaussTestNA) return ComradeBase.intensitymap_analytic!(img, m) end - -function ComradeBase.visibilitymap_numeric(m::GaussTestNA, p::ComradeBase.AbstractSingleDomain) +function ComradeBase.visibilitymap_numeric(m::GaussTestNA, + p::ComradeBase.AbstractSingleDomain) return ComradeBase.visibilitymap_analytic(m, p) end @@ -57,19 +57,16 @@ function ComradeBase.visibilitymap_numeric!(vis, m::GaussTestNA) return ComradeBase.visibilitymap_analytic!(vis, m) end - - ComradeBase.flux(::GaussTestNA{T}) where {T} = one(T) -ComradeBase.radialextent(::GaussTestNA{T}) where {T} = 5*one(T) - +ComradeBase.radialextent(::GaussTestNA{T}) where {T} = 5 * one(T) @testset "visibilities" begin - u = 0.1*randn(60) - v = 0.1*randn(60) + u = 0.1 * randn(60) + v = 0.1 * randn(60) ti = collect(Float64, 1:60) fr = fill(230e9, 60) m = GaussTest() - p = (;U=u, V=v, Ti=ti, Fr=fr) + p = (; U=u, V=v, Ti=ti, Fr=fr) g = UnstructuredDomain(p) @test visibilitymap(m, g) ≈ ComradeBase.visibilitymap_analytic(m, g) @test amplitudemap(m, g) ≈ abs.(ComradeBase.visibilitymap_analytic(m, g)) @@ -79,12 +76,12 @@ ComradeBase.radialextent(::GaussTestNA{T}) where {T} = 5*one(T) end @testset "visibilities not analytic" begin - u = 0.1*randn(60) - v = 0.1*randn(60) + u = 0.1 * randn(60) + v = 0.1 * randn(60) ti = collect(Float64, 1:60) fr = fill(230e9, 60) m = GaussTestNA() - p = (;U=u, V=v, Ti=ti, Fr=fr) + p = (; U=u, V=v, Ti=ti, Fr=fr) g = UnstructuredDomain(p) @test visibilitymap(m, g) ≈ ComradeBase.visibilitymap_analytic(m, g) @test amplitudemap(m, g) ≈ abs.(ComradeBase.visibilitymap_analytic(m, g))