diff --git a/Project.toml b/Project.toml index 615244c..cb232cf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ComradeBase" uuid = "6d8c423b-a35f-4ef1-850c-862fe21f82c4" authors = ["Paul Tiede and contributors"] -version = "0.9.3" +version = "0.9.4" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/src/ComradeBase.jl b/src/ComradeBase.jl index d18e503..b855ff7 100644 --- a/src/ComradeBase.jl +++ b/src/ComradeBase.jl @@ -19,9 +19,11 @@ export visibility, flux, fieldofview, imagepixels, pixelsizes, IntensityMap, named_dims + include("interface.jl") include("executors.jl") include("domain.jl") +include("multidomain.jl") include("rectigrid.jl") include("unstructured/domain.jl") include("unstructured/map.jl") @@ -32,7 +34,9 @@ const FluxMap2{T, N, E} = Union{ UnstructuredMap{T, <:AbstractVector, E}, } + include("visibilities.jl") +include("modifiers.jl") # include("rrules.jl") @setup_workload begin diff --git a/src/domain.jl b/src/domain.jl index c6e3ac4..5f530bf 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -177,6 +177,12 @@ function MinimalHeader(source, ra, dec, mjd, freq) return MinimalHeader(source, raT, decT, mjdT, freqT) end +function DimensionalData.val(m::AbstractHeader) + n = propertynames(m) + pm = Base.Fix1(getproperty, m) + return NamedTuple{n}(map(pm, n)) +end + """ NoHeader diff --git a/src/fourierdual.jl b/src/fourierdual.jl new file mode 100644 index 0000000..89862c7 --- /dev/null +++ b/src/fourierdual.jl @@ -0,0 +1,60 @@ +""" + $(TYPEDEF) + +This defines an abstract cache that can be used to +hold or precompute some computations. +""" +abstract type AbstractFourierDualDomain <: AbstractDomain end + +abstract type FourierTransform end + + +forward_plan(g::AbstractFourierDualDomain) = getfield(g, :plan_forward) +reverse_plan(g::AbstractFourierDualDomain) = getfield(g, :plan_reverse) +imgdomain(g::AbstractFourierDualDomain) = getfield(g, :imgdomain) +visdomain(g::AbstractFourierDualDomain) = getfield(g, :visdomain) +algorithm(g::AbstractFourierDualDomain) = getfield(g, :algorithm) + +EnzymeRules.inactive(::typeof(forward_plan), args...) = nothing +EnzymeRules.inactive(::typeof(reverse_plan), args...) = nothing + +abstract type AbstractPlan end +getplan(p::AbstractPlan) = getfield(p, :plan) +getphases(p::AbstractPlan) = getfield(p, :phases) +EnzymeRules.inactive(::typeof(getplan), args...) = nothing +EnzymeRules.inactive(::typeof(getphases), args...) = nothing + +function create_plans(algorithm, imgdomain, visdomain) + plan_forward = create_forward_plan(algorithm, imgdomain, visdomain) + plan_reverse = inverse_plan(plan_forward) + return plan_forward, plan_reverse +end + +function create_vismap(arr::AbstractArray, g::AbstractFourierDualDomain) + return ComradeBase.create_map(arr, visdomain(g)) +end + +function create_imgmap(arr::AbstractArray, g::AbstractFourierDualDomain) + return ComradeBase.create_map(arr, imgdomain(g)) +end + +function visibilitymap_analytic(m::AbstractModel, grid::AbstractFourierDualDomain) + return visibilitymap_analytic(m, visdomain(grid)) +end + +function visibilitymap_numeric(m::AbstractModel, grid::AbstractFourierDualDomain) + img = intensitymap_analytic(m, imgdomain(grid)) + vis = applyft(forward_plan(grid), img) + return vis +end + +function intensitymap_analytic(m::AbstractModel, grid::AbstractFourierDualDomain) + return intensitymap_analytic(m, imgdomain(grid)) +end + +function intensitymap_numeric(m::AbstractModel, grid::AbstractFourierDualDomain) + # This is because I want to make a grid that is the same size as the image + # so we revert to the standard method and not what ever was cached + img = intensitymap_numeric(m, imgdomain(grid)) + return img +end diff --git a/src/modifiers.jl b/src/modifiers.jl new file mode 100644 index 0000000..5ad7887 --- /dev/null +++ b/src/modifiers.jl @@ -0,0 +1,553 @@ +export stretched, shifted, rotated, renormed, modify, Stretch, Renormalize, Shift, Rotate, + unmodified, basemodel + +""" + $(TYPEDEF) + +General type for a model modifier. These transform any model +using simple Fourier transform properties. To modify a model +you can use the [`ModifiedModel`](@ref) constructor or the [`modify`](@ref) +function. + +```julia-repl +julia> visanalytic(stretched(Disk(), 2.0, 2.0)) == visanalytic(Disk()) +true +``` + + + +To implement a model transform you need to specify the following methods: +- [`transform_uv`](@ref) +- [`transform_image`](@ref) +- [`scale_uv`](@ref) +- [`scale_image`](@ref) +- [`radialextent`](@ref) +See thee docstrings of those methods for guidance on implementation details. + +Additionally these methods assume the modifiers are of the form + +I(x,y) -> fᵢ(x,y)I(gᵢ(x,y)) +V(u,v) -> fᵥ(u,v)V(gᵥ(u,v)) + +where `g` are the transform_image/uv functions and `f` are the scale_image/uv +function. + +""" +abstract type ModelModifier{T} end + +""" + scale_image(model::AbstractModifier, x, y) + +Returns a number of how to to scale the image intensity at `x` `y` for an modified `model` +""" +function scale_image end + +""" + transform_image(model::AbstractModifier, x, y) + +Returns a transformed `x` and `y` according to the `model` modifier +""" +function transform_image end + +""" + scale_image(model::AbstractModifier, u, u) + +Returns a number on how to scale the image visibility at `u` `v` for an modified `model` +""" +function scale_uv end + +""" + transform_uv(model::AbstractModifier, u, u) + +Returns a transformed `u` and `v` according to the `model` modifier +""" +function transform_uv end + +unitscale(T, ::NotPolarized) = one(T) +unitscale(T, ::IsPolarized) = I + +""" + $(TYPEDEF) + +Container type for models that have been transformed in some way. +For a list of potential modifiers or transforms see `subtypes(ModelModifiers)`. + +# Fields +$(FIELDS) +""" +struct ModifiedModel{M, MT <: Tuple} <: AbstractModel + """base model""" + model::M + """model transforms""" + transform::MT +end + +function Base.show(io::IO, mi::ModifiedModel) + #io = IOContext(io, :compact=>true) + #s = summary(mi) + println(io, "ModifiedModel") + println(io, " base model: ", mi.model) + println(io, " Modifiers:") + for i in eachindex(mi.transform)[1:(end - 1)] + println(io, " $i. ", summary(mi.transform[i])) + end + return print(io, " $(length(mi.transform)). ", summary(mi.transform[end])) +end + +""" + unmodified(model::ModifiedModel) + +Returns the un-modified model + +### Example +```julia-repl +julia> m = stretched(rotated(Gaussian(), π/4), 2.0, 1.0) +julia> umodified(m) == Gaussian() +true +``` +""" +unmodified(model::ModifiedModel) = model.model +unmodified(model::AbstractModel) = model + +""" + basemodel(model::ModifiedModel) + +Returns the ModifiedModel with the last transformation stripped. + +# Example +```julia-repl +julia> basemodel(stretched(Disk(), 1.0, 2.0)) == Disk() +true +``` +""" +basemodel(model::ModifiedModel) = ModifiedModel(model.model, Base.front(model.transform)) +basemodel(model::ModifiedModel{M, Tuple{}}) where {M} = model + +flux(m::ModifiedModel) = flux(m.model) * mapreduce(flux, Base.:*, m.transform) +flux(::ModelModifier{T}) where {T} = one(T) + +radialextent(m::ModifiedModel) = radialextent_modified(radialextent(m.model), m.transform) + +@inline visanalytic(::Type{ModifiedModel{M, T}}) where {M, T} = visanalytic(M) +@inline imanalytic(::Type{ModifiedModel{M, T}}) where {M, T} = imanalytic(M) +@inline ispolarized(::Type{ModifiedModel{M, T}}) where {M, T} = ispolarized(M) + +@inline function ModifiedModel(m, t::ModelModifier) + return ModifiedModel(m, (t,)) +end + +@inline function ModifiedModel(m::ModifiedModel, t::ModelModifier) + model = m.model + t0 = m.transform + return ModifiedModel(model, (t0..., t)) +end + +""" + Base.:∘(m::AbstractModel, t::ModelModifier) + +Modifiers act as transformations on the inputs plus a jacobian term to ensure the +transform is an isometry. Therefore we provide the syntactic sugar of using the +`∘` operator to apply a modifier to a model. This is equivalent to calling `modify(m, t)`. + +## Example +```julia-repl +julia> m = Gaussian() +julia> t = Shift(1.0, 2.0) +julia> m ∘ t == modify(m, t) +true +``` + +Note that because this is composition the order is reversed. So `m ∘ t1 ∘ t2` is equivalent to +`modify(m, t2, t1)`. + +```julia-repl +julia> m ∘ Shift(1.0, 2.0) ∘ Rotate(π/4) == modify(m, Rotate(π/4), Shift(1.0, 2.0)) +true +``` + +""" +@inline Base.:∘(m::AbstractModel, t::ModelModifier) = modify(m, t) +@inline function Base.:∘(m::ModifiedModel, t::ModelModifier) + trf = m.transform + return ModifiedModel(m.model, (t, trf...)) +end + + +@inline doesnot_uv_modify(t::Tuple) = doesnot_uv_modify(Base.front(t)) * + doesnot_uv_modify(last(t)) +@inline doesnot_uv_modify(::Tuple{}) = true + +@inline function modify_uv(model, t::Tuple, p, scale) + pt = transform_uv(model, last(t), p) + scalet = scale_uv(model, last(t), p) + return modify_uv(model, Base.front(t), pt, scalet * scale) +end +@inline modify_uv(model, ::Tuple{}, p, scale) = p, scale + +@inline function modify_image(model, t::Tuple, p, scale) + pt = transform_image(model, last(t), p) + scalet = scale_image(model, last(t), p) + return modify_image(model, Base.front(t), pt, scalet * scale) +end +@inline modify_image(model, ::Tuple{}, p, scale) = p, scale + + +@inline radialextent_modified(r::Real, t::Tuple) = radialextent_modified( + radialextent_modified( + r, + last(t) + ), + Base.front(t) +) +@inline radialextent_modified(r::Real, ::Tuple{}) = r + +""" + modify(m::AbstractModel, transforms...) + +Modify a given `model` using the set of `transforms`. This is the most general +function that allows you to apply a sequence of model transformation for example + +```julia-repl +modify(Gaussian(), Stretch(2.0, 1.0), Rotate(π/4), Shift(1.0, 2.0), Renorm(2.0)) +``` +will create a asymmetric Gaussian with position angle `π/4` shifted to the position +(1.0, 2.0) with a flux of 2 Jy. This is similar to Flux's chain. +""" +function modify(m::AbstractModel, transforms...) + return ModifiedModel(m, transforms) +end + +@inline function visibility_point(m::M, p) where {M <: ModifiedModel} + mbase = m.model + transform = m.transform + ispol = ispolarized(M) + pt, scale = modify_uv(ispol, transform, p, unitscale(Complex{eltype(p.U)}, ispol)) + return scale * visibility_point(mbase, pt) +end + +@inline function intensity_point(m::M, p) where {M <: ModifiedModel} + mbase = m.model + transform = m.transform + ispol = ispolarized(M) + pt, scale = modify_image(ispol, transform, p, unitscale(eltype(p.X), ispol)) + return scale * intensity_point(mbase, pt) +end + + +function visibilitymap_analytic(m::ModifiedModel, p::AbstractSingleDomain) + vis = allocate_vismap(m, p) + visibilitymap_analytic!(vis, m) + return vis +end + + +""" + $(TYPEDEF) + +Shifts the model by `Δx` units in the x-direction and `Δy` units +in the y-direction. + +# Example +```julia-repl +julia> modify(Gaussian(), Shift(2.0, 1.0)) == shifted(Gaussian(), 2.0, 1.0) +true +``` +""" +struct Shift{T, T1, T2} <: ModelModifier{T} + Δx::T1 + Δy::T2 +end + + +function Shift(a::Number, b::Number) + T = promote_type(typeof(a), typeof(b)) + return Shift{T, typeof(a), typeof(b)}(a, b) +end + +function Shift(a::DomainParams{P}, b) where {P} + T = promote_type(P, typeof(b)) + return Shift{T, typeof(a), typeof(b)}(a, b) +end + +function Shift(a, b::DomainParams{P}) where {P} + T = promote_type(P, typeof(a)) + return Shift{T, typeof(a), typeof(b)}(a, b) +end + +function Shift(a::DomainParams{P1}, b::DomainParams{P2}) where {P1, P2} + T = promote_type(P1, P2) + return Shift{T, typeof(a), typeof(b)}(a, b) +end + +""" + $(SIGNATURES) + +Shifts the model `m` in the image domain by an amount `Δx,Δy` +in the x and y directions respectively. +""" +shifted(model, Δx, Δy) = ModifiedModel(model, Shift(Δx, Δy)) + +doesnot_uv_modify(::Shift) = true + +# This is a simple overload to simplify the type system +@inline radialextent_modified(r::Real, t::Shift) = r + max(abs(t.Δx), abs(t.Δy)) + +@inline function transform_image(model, transform::Shift, p) + @unpack_params Δx, Δy = transform(p) + X = p.X - Δx + Y = p.Y - Δy + return update_spat(p, X, Y) +end + +@inline transform_uv(model, ::Shift, p) = p + +@inline scale_image(m, ::Shift, p) = unitscale(p.X, m) +@inline @fastmath function scale_uv(m, transform::Shift, p) + @unpack_params Δx, Δy = transform(p) + (; U, V) = p + T = typeof(Δx) + return cispi( + 2 * + (U * Δx + V * Δy) + ) * + unitscale(T, m) +end + +""" + $(TYPEDEF) + +Renormalizes the flux of the model to the new value `scale*flux(model)`. +We have also overloaded the Base.:* operator as syntactic sugar +although I may get rid of this. + + +# Example +```julia-repl +julia> modify(Gaussian(), Renormalize(2.0)) == 2.0*Gaussian() +true +``` +""" +struct Renormalize{T} <: ModelModifier{T} + scale::T +end + +""" + $(SIGNATURES) + +Renormalizes the model `m` to have total flux `f*flux(m)`. +This can also be done directly by calling `Base.:*` i.e., + +```julia-repl +julia> renormed(m, f) == f*M +true +``` +""" +renormed(model::M, f) where {M <: AbstractModel} = ModifiedModel(model, Renormalize(f)) +@inline doesnot_uv_modify(::Renormalize) = true + +const ModNum = Union{Number, DomainParams} + +Base.:*(model::AbstractModel, f::ModNum) = renormed(model, f) +Base.:*(f::ModNum, model::AbstractModel) = renormed(model, f) +Base.:/(f::ModNum, model::AbstractModel) = renormed(model, inv(f)) +Base.:/(model::AbstractModel, f::ModNum) = renormed(model, inv(f)) +# Dispatch on RenormalizedModel so that I just make a new RenormalizedModel with a different f +# This will make it easier on the compiler. +# Base.:*(model::ModifiedModel, f::Number) = renormed(model.model, model.scale*f) +# Overload the unary negation operator to be the same model with negative flux +Base.:-(model::AbstractModel) = renormed(model, -1) +flux(t::Renormalize) = t.scale + +@inline transform_image(m, ::Renormalize, p) = p +@inline transform_uv(m, ::Renormalize, p) = p + +@inline function scale_image(m, transform::Renormalize, p) + @unpack_params scale = transform(p) + return scale * unitscale(typeof(scale), m) +end + +@inline function scale_uv(m, transform::Renormalize, p) + @unpack_params scale = transform(p) + return scale * unitscale(typeof(scale), m) +end + +@inline radialextent_modified(r::Real, ::Renormalize) = r + +""" + Stretch(α, β) + Stretch(r) + +Stretched the model in the x and y directions, i.e. the new intensity is + Iₛ(x,y) = 1/(αβ) I(x/α, y/β), +where were renormalize the intensity to preserve the models flux. + +# Example +```julia-repl +julia> modify(Gaussian(), Stretch(2.0)) == stretched(Gaussian(), 2.0, 1.0) +true +``` + +If only a single argument is given it assumes the same stretch is applied in both direction. + +```julia-repl +julia> Stretch(2.0) == Stretch(2.0, 2.0) +true +``` +""" +struct Stretch{T, T1, T2} <: ModelModifier{T} + α::T1 + β::T2 +end + +function Stretch(a::Number, b::Number) + T = promote_type(typeof(a), typeof(b)) + return Stretch{T, typeof(a), typeof(b)}(a, b) +end + +function Stretch(a::DomainParams{P}, b) where {P} + T = promote_type(P, typeof(b)) + return Stretch{T, typeof(a), typeof(b)}(a, b) +end + +function Stretch(a, b::DomainParams{P}) where {P} + T = promote_type(P, typeof(a)) + return Stretch{T, typeof(a), typeof(b)}(a, b) +end + +function Stretch(a::DomainParams{P1}, b::DomainParams{P2}) where {P1, P2} + T = promote_type(P1, P2) + return Stretch{T, typeof(a), typeof(b)}(a, b) +end + +Stretch(r) = Stretch(r, r) + +""" + $(SIGNATURES) + +Stretches the model `m` according to the formula + Iₛ(x,y) = 1/(αβ) I(x/α, y/β), +where were renormalize the intensity to preserve the models flux. +""" +stretched(model, α, β) = ModifiedModel(model, Stretch(α, β)) +stretched(model, α) = stretched(model, α, α) + +@inline doesnot_uv_modify(::Stretch) = false +@inline function transform_image(m, transform::Stretch, p) + (; X, Y) = p + @unpack_params α, β = transform(p) + # @show p + pt = update_spat(p, X / α, Y / β) + # @show pt + return pt +end + +@inline function transform_uv(m, transform::Stretch, p) + (; U, V) = p + @unpack_params α, β = transform(p) + return update_spat(p, U * α, V * β) +end + +@inline function scale_image(m, transform::Stretch, p) + @unpack_params α, β = transform(p) + T = typeof(α) + return inv(α * β) * unitscale(T, m) +end + +@inline scale_uv(m, tr::Stretch, p) = unitscale(typeof(getparam(tr, :α, p)), m) + +@inline radialextent_modified(r::Real, t::Stretch) = r * max(t.α, t.β) + +""" + Rotate(ξ) + +Type for the rotated model. This is more fine grained constrol of +rotated model. + +# Example +```julia-repl +julia> modify(Gaussian(), Rotate(2.0)) == rotated(Gaussian(), 2.0) +true +``` +""" +struct Rotate{T} <: ModelModifier{T} + s::T + c::T + function Rotate(ξ::F) where {F <: Real} + s, c = sincos(ξ) + return new{F}(s, c) + end + function Rotate(ξ::DomainParams) + return new{typeof(ξ)}(ξ, ξ) + end +end + +function getparam( + m::Rotate{T}, + s::Symbol, + p + ) where {T <: DomainParams} + m = getproperty(m, s) + mr = Rotate(build_param(m, p)) + return getproperty(mr, s) +end + +""" + $(SIGNATURES) + +Rotates the model by an amount `ξ` in radians in the clockwise direction. +""" +rotated(model, ξ) = ModifiedModel(model, Rotate(ξ)) + +""" + $(SIGNATURES) + +Returns the rotation angle of the rotated `model` +""" +posangle(model::Rotate) = atan(model.s, model.c) + +@inline doesnot_uv_modify(::Rotate) = false + +@inline function transform_image(m, transform::Rotate, p) + @unpack_params s, c = transform(p) + (; X, Y) = p + Xr = c * X - s * Y + Yr = s * X + c * Y + pt = update_spat(p, Xr, Yr) + return pt +end + +@inline function transform_uv(m, transform::Rotate, p) + @unpack_params s, c = transform(p) + (; U, V) = p + Ur = c * U - s * V + Vr = s * U + c * V + return update_spat(p, Ur, Vr) +end + +@inline scale_image(::NotPolarized, model::Rotate, p) = one(typeof(getparam(model, :s, p))) + +@inline function spinor2_rotate(c, s) + u = oneunit(c) + z = zero(s) + c2 = c^2 - s^2 + s2 = 2 * c * s + return SMatrix{4, 4}( + u, z, z, z, + z, c2, s2, z, + z, -s2, c2, z, + z, z, z, u + ) +end + +@inline function scale_image(::IsPolarized, model::Rotate, p) + @unpack_params s, c = model(p) + return spinor2_rotate(c, s) +end + +@inline scale_uv(::NotPolarized, model::Rotate, p) = one(typeof(getparam(model, :s, p))) + +@inline function scale_uv(::IsPolarized, model::Rotate, p) + @unpack_params s, c = model(p) + return spinor2_rotate(c, s) +end +@inline radialextent_modified(r::Real, ::Rotate) = r diff --git a/src/multidomain.jl b/src/multidomain.jl index fc4f435..1c00bcd 100644 --- a/src/multidomain.jl +++ b/src/multidomain.jl @@ -11,21 +11,33 @@ The interface is simple and to extend this with your own time and frequency mode most users will just need to define ```julia -struct MyDomainParam <: DomainParams end -function build_params(param::MyDomainParam, p) +struct MyDomainParam{T} <: DomainParams{T} end +function build_param(param::MyDomainParam{Float64}, p) ... end where `p` is the point where the model will be evaluated at. For an example see the [`TaylorSpectralModel`](@ref). +To evaluate the parameter family at a point `p` in the frequency and time +domain use `build_param(param, p)` or just `param(p)`. + For a model parameterized with a `<:DomainParams` the a use should access the parameters with [`getparam`](@ref) or the `@unpack_params` macro. ``` """ -abstract type DomainParams end +abstract type DomainParams{T} end + +abstract type FrequencyParams{T} <: DomainParams{T} end +abstract type TimeParams{T} <: DomainParams{T} end -abstract type FrequencyParams <: DomainParams end -abstract type TimeParams <: DomainParams end +""" + paramtype(::Type{<:DomainParams}) + +Computes the base parameter type of the DomainParams. If `!<:DomainParams` then it just returns +the type. +""" +@inline paramtype(::Type{<:DomainParams{T}}) where {T} = paramtype(T) +@inline paramtype(T::Type{<:Any}) = T """ getparam(m, s::Symbol, p) @@ -35,6 +47,10 @@ This is similar to getproperty, but allows for the parameter to be a function of domain. Essentially is `m.s <: DomainParams` then `m.s` is evaluated at the parameter `p`. If `m.s` is not a subtype of `DomainParams` then `m.s` is returned. +!!! warn + Developers should not typically overload this function and instead + target [`build_param`](@ref). + !!! warn This feature is experimental and is not considered part of the public stable API. @@ -47,6 +63,66 @@ end return getparam(m, s, p) end +""" + build_param(param::DomainParams, p) + +Constucts the parameters for the `param` model at the point `p` +in the (X/U, Y/V, Ti, Fr) domain. This is a required function for +any `<: DomainParams` and must return a number for the specific +parameter at the point `p`. +""" @inline function build_param(param::Any, p) return param end + +function build_param(param::NTuple, p) + return map(x -> build_param(x, p), param) +end + +function build_param(param::AbstractArray, p) + return map(x -> build_param(x, p), param) +end + +function (m::DomainParams)(p) + return build_param(m, p) +end + +""" + @unpack_params a,b,c,... = m(p) + +Extracts the parameters `a,b,c,...` from the model `m` evaluated at the domain `p`. +This is a macro that essentially lowers to +```julia +a = getparam(m, :a, p) +b = getparam(m, :b, p) +... +``` +For any model that may depend on a `DomainParams` type this macro should be used to +extract the parameters. + +!!! warn + This feature is experimental and is not considered part of the public stable API. + +""" +macro unpack_params(args) + args.head != :(=) && + throw(ArgumentError("Expression needs to be of the form a, b, = c(p)")) + items, suitcase = args.args + items = isa(items, Symbol) ? [items] : items.args + hasproperty(suitcase, :head) || + throw(ArgumentError("RHS of expression must be of form m(p)")) + suitcase.head != :call && throw(ArgumentError("RHS of expression must be of form m(p)")) + m, p = suitcase.args[1], suitcase.args[2] + paraminstance = gensym() + kp = [ + :($key = getparam($paraminstance, Val{$(Expr(:quote, key))}(), $p)) + for key in items + ] + kpblock = Expr(:block, kp...) + expr = quote + local $paraminstance = $m + $kpblock + $paraminstance + end + return esc(expr) +end diff --git a/test/runtests.jl b/test/runtests.jl index c8252ef..89325bd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using OhMyThreads using Enzyme using KernelAbstractions using Polyester +using JET using FiniteDifferences # using ChainRulesCore diff --git a/test/visibilities.jl b/test/visibilities.jl index 0ccba35..db14c4f 100644 --- a/test/visibilities.jl +++ b/test/visibilities.jl @@ -101,3 +101,91 @@ end logclosure_amplitudemap(m, g, g, g, g) @test angle.(bispectrummap(m, g, g, g)) ≈ closure_phasemap(m, g, g, g) end + +@testset "Modifiers" begin + u = randn(100) * 0.5 + v = randn(100) * 0.5 + t = sort(rand(100) * 0.5) + f = fill(230.0e9, 100) + guv = UnstructuredDomain((U = u, V = v, Ti = t, Fr = f)) + gim = imagepixels(10.0, 10.0, 64, 64) + ma = GaussTest() + + @testset "Shifted" begin + mas = shifted(ma, 0.1, 0.1) + gims = imagepixels(10.0, 10.0, 64, 64, -0.1, -0.1) + @test ComradeBase.intensity_point(ma, (X = 0.5, Y = 0.5)) ≈ + ComradeBase.intensity_point(mas, (X = 0.6, Y = 0.6)) + + @test ComradeBase.visibility_point(ma, (U = 0.1, V = 0.1)) ≈ + ComradeBase.visibility_point(mas, (U = 0.1, V = 0.1)) * exp(-2π * 1im * (0.1 * 0.1 + 0.1 * 0.1)) + + @test baseimage(intensitymap(ma, gim)) ≈ baseimage(intensitymap(mas, gims)) + end + + @testset "Renormed" begin + m1 = 3.0 * ma + m2 = ma * 3.0 + m2inv = ma / (1 / 3) + p = (U = 4.0, V = 0.0) + @test visibility(m1, p) == visibility(m2, p) + @test ComradeBase.intensity_point(m1, (X = 0.5, Y = 0.5)) ≈ + ComradeBase.intensity_point(m2, (X = 0.5, Y = 0.5)) + @test ComradeBase.intensity_point(m1, (X = 0.5, Y = 0.5)) ≈ + ComradeBase.intensity_point(m2inv, (X = 0.5, Y = 0.5)) + + @test intensitymap(m1, gim) ≈ intensitymap(m2, gim) + @test visibilitymap(m1, guv) ≈ visibilitymap(m2, guv) + end + + @testset "Stretched" begin + mas = stretched(ma, 5.0, 4.0) + gims = imagepixels(10.0 * 5.0, 10.0 * 4.0, 64, 64) + @test ComradeBase.intensity_point(mas, (X = 0.5, Y = 0.5)) ≈ + ComradeBase.intensity_point(ma, (X = 0.5 / 5, Y = 0.5 / 4)) / 20 + + @test ComradeBase.visibility_point(mas, (U = 0.1, V = 0.1)) ≈ + ComradeBase.visibility_point(ma, (U = 0.1 * 5, V = 0.1 * 4)) + + @test baseimage(intensitymap(ma, gim)) ≈ baseimage(intensitymap(mas, gims)) + @test ComradeBase.radialextent(mas) ≈ ComradeBase.radialextent(ma) * 5.0 + @test flux(mas) ≈ flux(ma) + + end + + @testset "Rotated" begin + ma = stretched(GaussTest(), 2.0, 1.0) + mar = rotated(ma, π / 3) + grs = imagepixels(10.0, 10.0, 64, 64; posang = π / 3) + + @test ComradeBase.intensity_point(mar, (X = 0.5, Y = 0.5)) ≈ + ComradeBase.intensity_point( + ma, + ( + X = 0.5 * cos(π / 3) - 0.5 * sin(π / 3), + Y = 0.5 * sin(π / 3) + 0.5 * cos(π / 3), + ) + ) + + @test ComradeBase.visibility_point(mar, (U = 0.1, V = 0.1)) ≈ + ComradeBase.visibility_point( + ma, ( + U = 0.1 * cos(π / 3) - 0.1 * sin(π / 3), + V = 0.1 * sin(π / 3) + 0.1 * cos(π / 3), + ) + ) + + @test baseimage(intensitymap(ma, gim)) ≈ baseimage(intensitymap(mar, grs)) + + end + + @testset "Compose" begin + m = GaussTest() + m1 = 5 * (m ∘ Rotate(π / 3) ∘ Shift(0.1, 0.1) ∘ Stretch(5.0, 4.0)) + m2 = 5 * modify(m, Stretch(5.0, 4.0), Shift(0.1, 0.1), Rotate(π / 3)) + + @test m1 == m2 + @test_opt 5 * (m ∘ Rotate(π / 3) ∘ Shift(0.1, 0.1) ∘ Stretch(5.0, 4.0)) + @test_opt 5 * modify(m, Stretch(5.0, 4.0), Shift(0.1, 0.1), Rotate(π / 3)) + end +end