diff --git a/docs/src/api.md b/docs/src/api.md index 4e9b90de..f5547af7 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -25,7 +25,7 @@ Comrade.caltable(::Comrade.JonesCache, ::AbstractVector) Comrade.caltable(::Comrade.EHTObservation, ::AbstractVector) Comrade.DesignMatrix Comrade.JonesCache -Comrade.TransformCache +Comrade.ResponseCache Comrade.JonesModel Comrade.VLBIModel Comrade.CalPrior diff --git a/examples/imaging_pol.jl b/examples/imaging_pol.jl index 05ef25bb..635db488 100644 --- a/examples/imaging_pol.jl +++ b/examples/imaging_pol.jl @@ -199,7 +199,7 @@ skymeta = (;K, cache, grid) # # First, we will define our deterministic transform cache. Note that this dataset has need # been pre-corrected for feed rotation, so we need to add those into the `tcache`. -tcache = TransformCache(dvis; add_fr=true, ehtim_fr_convention=false) +tcache = ResponseCache(dvis; add_fr=true, ehtim_fr_convention=false) #- # Next we define our cache that maps quantities e.g., gain products, that change from scan-to-scan. scancache = jonescache(dvis, ScanSeg()) diff --git a/ext/ComradeMakieExt.jl b/ext/ComradeMakieExt.jl index feb467ff..ae2b971d 100644 --- a/ext/ComradeMakieExt.jl +++ b/ext/ComradeMakieExt.jl @@ -7,7 +7,4 @@ else using ..Makie end - - - end diff --git a/ext/ComradeVIDAExt.jl b/ext/ComradeVIDAExt.jl new file mode 100644 index 00000000..1a313b44 --- /dev/null +++ b/ext/ComradeVIDAExt.jl @@ -0,0 +1,14 @@ +module ComradeVIDAExt + +using Comrade +if isdefined(Base, :get_extension) + using VIDA +else + using ..VIDA +end + +function Comrade.template_align(div::VIDA.AbstractDivergence, m, lower, upper) + +end + +end diff --git a/lib/ComradeAHMC/src/ComradeAHMC.jl b/lib/ComradeAHMC/src/ComradeAHMC.jl index 3cf916e8..1516e68d 100644 --- a/lib/ComradeAHMC/src/ComradeAHMC.jl +++ b/lib/ComradeAHMC/src/ComradeAHMC.jl @@ -443,7 +443,13 @@ end function load_table(out::String, indices::Union{Base.Colon, UnitRange, StepRange}=Base.Colon(); table="samples") @assert isdir(abspath(out)) "$out is not a directory. This isn't where the HMC samples are stored" @assert isfile(joinpath(abspath(out), "parameters.jld2")) "parameters.jld2 " - return load_table(load(joinpath(abspath(out), "parameters.jld2"), "params"), indices; table) + p = load(joinpath(abspath(out), "parameters.jld2"), "params") + if p.filename != abspath(out) + @warn "filename stored in params does not equal what was passed\n"* + "we will load the path passed\n $(out)." + p = DiskOutput(out, p.nfiles, p.stride, p.nsamples) + end + return load_table(p, indices; table) end diff --git a/src/align.jl b/src/align.jl new file mode 100644 index 00000000..2b2cea94 --- /dev/null +++ b/src/align.jl @@ -0,0 +1,3 @@ +export template_align + +function template_align end diff --git a/src/calibration/calibration.jl b/src/calibration/calibration.jl index dbf2392a..ebc6fdb6 100644 --- a/src/calibration/calibration.jl +++ b/src/calibration/calibration.jl @@ -13,7 +13,7 @@ end """ JonesModel(jones::JonesPairs, refbasis = CirBasis()) - JonesModel(jones::JonesPairs, tcache::TransformCache) + JonesModel(jones::JonesPairs, tcache::ResponseCache) Constructs the intrument corruption model using pairs of jones matrices `jones` and a reference basis diff --git a/src/calibration/jones.jl b/src/calibration/jones.jl index e8cecd5d..6882e012 100644 --- a/src/calibration/jones.jl +++ b/src/calibration/jones.jl @@ -1,5 +1,5 @@ export JonesCache, TrackSeg, ScanSeg, FixedSeg, IntegSeg, jonesG, jonesD, jonesT, - TransformCache, JonesModel, jonescache, station_tuple, jonesmap + ResponseCache, JonesModel, jonescache, station_tuple, jonesmap """ $(TYPEDEF) @@ -879,7 +879,7 @@ on sky reference basis. # Fields $(FIELDS) """ -struct TransformCache{M, B<:PolBasis} <: AbstractJonesCache +struct ResponseCache{M, B<:PolBasis} <: AbstractJonesCache """ Transform matrices for the first stations """ @@ -895,7 +895,7 @@ struct TransformCache{M, B<:PolBasis} <: AbstractJonesCache end """ - TransformCache(obs::EHTObservation; add_fr=true, ehtim_fr_convention=false, ref::PolBasis=CirBasis()) + ResponseCache(obs::EHTObservation; add_fr=true, ehtim_fr_convention=false, ref::PolBasis=CirBasis()) Constructs the cache that holds the transformation from the **chosen** on-sky reference basis to the basis that the telescope measures the electric fields given an observation `obs`. @@ -904,36 +904,47 @@ are included in the transformation cache with the `add_fr` toggle. The user can using the keyword argument `ref` which is the (R,L) circular basis by default. # Notes -We use the following definition for our feed rotations +We use the following definition for our feed rotations in circular basis ``` exp(-iθ) 0 0 exp(iθ) ``` +and in the linear feed basis + +``` + cos(θ) sin(θ) + -sin(θ) cos(θ) +``` + # Warning eht-imaging can sometimes pre-rotate the coherency matrices. As a result the field rotation can sometimes be applied twice. To compensate for this we have added a `ehtim_fr_convention` which will fix this. """ -function TransformCache(obs::EHTObservation; add_fr=true, ehtim_fr_convention=false, ref::PolBasis=CirBasis()) +function ResponseCache(obs::EHTObservation; add_fr=true, ehtim_fr_convention=false, ref::PolBasis=CirBasis()) T1 = StructArray(map(x -> basis_transform(ref, x[1]), obs.data.polbasis)) T2 = StructArray(map(x -> basis_transform(ref, x[2]), obs.data.polbasis)) + # Our feed rotation matrix from extract_FRs always returns the rotation according to + # a circular basis. We then transform to a linear basis before constructing the basis + Tcirc1 = StructArray(map(x -> basis_transform(CirBasis(), x[1]), obs.data.polbasis)) + Tcirc2 = StructArray(map(x -> basis_transform(CirBasis(), x[2]), obs.data.polbasis)) if add_fr field_rotations = extract_FRs(obs; ehtim_fr_convention) - T1 .= field_rotations.m1.*T1 - T2 .= field_rotations.m2.*T2 + @. T1 .= Tcirc1*field_rotations.m1*adjoint(Tcirc1)*T1 + @. T2 .= Tcirc2*field_rotations.m2*adjoint(Tcirc2)*T2 end - return TransformCache{typeof(T1), typeof(ref)}(T1, T2, ref) + return ResponseCache{typeof(T1), typeof(ref)}(T1, T2, ref) end """ - jonesT(tcache::TransformCache) + jonesT(tcache::ResponseCache) Returns a `JonesPair` of matrices that transform from the model coherency matrices basis to the on-sky coherency basis, this includes the feed rotation and choice of polarization feeds. """ -jonesT(tcache::TransformCache) = JonesPairs(tcache.T1, tcache.T2) -JonesModel(jones::J, tcache::TransformCache) where {J} = JonesModel(jones, tcache.refbasis) +jonesT(tcache::ResponseCache) = JonesPairs(tcache.T1, tcache.T2) +JonesModel(jones::J, tcache::ResponseCache) where {J} = JonesModel(jones, tcache.refbasis) """ @@ -1036,6 +1047,7 @@ function extract_FRs(obs::EHTObservation; ehtim_fr_convention=false) pars = tarr.fr_parallactic offs = tarr.fr_offset + # get station names bls = config.data.baseline ant1 = first.(bls) diff --git a/test/Core/bayes.jl b/test/Core/bayes.jl index e7a97663..b326a369 100644 --- a/test/Core/bayes.jl +++ b/test/Core/bayes.jl @@ -97,7 +97,7 @@ end @testset "RadioLikelihood" begin _,vis, amp, lcamp, cphase, coh = load_data() - tcache = TransformCache(coh) + tcache = ResponseCache(coh) lklhd_amp = RadioLikelihood(test_model, amp) lklhd_cp = RadioLikelihood(test_model, cphase) lklhd_lc = RadioLikelihood(test_model, lcamp) diff --git a/test/Core/gains.jl b/test/Core/gains.jl index 39c59e77..03e08c70 100644 --- a/test/Core/gains.jl +++ b/test/Core/gains.jl @@ -216,7 +216,7 @@ end scancache = jonescache(dcoh, ScanSeg()) scancache2 = jonescache(dcoh, segs) - tcache = TransformCache(dcoh; add_fr=true) + tcache = ResponseCache(dcoh; add_fr=true) @@ -282,7 +282,7 @@ end scancache = jonescache(dcoh, ScanSeg()) phasecache= jonescache(dcoh, ScanSeg(); autoref=SingleReference(:AA, 1.0+0.0im)) trackcache= jonescache(dcoh, TrackSeg()) - tcache = TransformCache(dcoh; add_fr=true) + tcache = ResponseCache(dcoh; add_fr=true) dga = CalPrior(station_tuple(dcoh, LogNormal(0.0, 0.1)), scancache) dgp = CalPrior(station_tuple(dcoh, Uniform(0.0, 2π)), phasecache) diff --git a/test/Core/observation.jl b/test/Core/observation.jl index 05bff8c0..d0d32d7a 100644 --- a/test/Core/observation.jl +++ b/test/Core/observation.jl @@ -89,7 +89,7 @@ using Pyehtim plot(m, amp) plot(m, lcamp) plot(m,cphase) - jt = jonesT(TransformCache(dcoh)) + jt = jonesT(ResponseCache(dcoh)) mp = Comrade.VLBIModel(JonesModel(jt), PolarizedModel(m, ZeroModel(), ZeroModel(), 0.1*Gaussian())) plot(mp, dcoh) #residual(m, vis)