Currently, Comrade uses a fixed way of defining Likelihood based on the data products. Moreover, all the model parameters are inside the model visibilities while the sigma is fixed based on the errors of the data products. If likelihood function was made modular, we can do variability modeling and multi-frequency modeling.

Currently, our map from gains to visibilities follows a time station order.
To explain this, note that our gain parameters are represented as a flat vector of gains. The order is such that we index over the stations first and then time stamps. This makes sense if you look at each time stamp independently, as our current gain priors typically enforce. However, this order is not ideal when dealing with time-correlated gains.

For time-correlated gains, storing the vectors in a station time ordering would be nice. This means we should first place group the gains for each station and then concatenate the times together. This will simplify certain gain orderings, and I can remove the need for relative gain segmentation, which is slow and not ideal.

The roadmap

• Define a gain ordering abstract type
• Move the current scheme to TimeStation()
• Implement the new scheme as StationTime()
• Remove the ScanSeg{true}() since it is now redundant
• Implement a time correlated gain prior to test this out.

Dev docs for example: https://ptiede.github.io/Comrade.jl/dev/examples/imaging_vis/

In chain, stats = sample(rng, post, AHMC(;metric, autodiff=AD.ZygoteBackend()), 400; nadapts=300, init_params=x0) , AHMC(;metric, autodiff=AD.ZygoteBackend()) does not work and should be changed to AHMC(;metric, autodiff=Val(:Zygote))

Comrade.jl/src/calibration/gains.jl

Line 70 to 81:

struct GainCache{D1<:DesignMatrix, D2<:DesignMatrix, T, S}
    m1::D1
    m2::D2
    times::T
    stations::S
end

function GainCache(st::ScanTable)
    gtime, gstat = gain_stations(st)
    m1, m2 = gain_design(st)
    return GainCache(m1, m2, gtime, gstat)
end

Could this be optimized by reusing the input visibility matrix and threading the calculation? Maybe something like this:

function phasecenter(vis, uu, vv, x0, y0, dx, dy, executor=SequentialEx())
    lenu = length(uu)
    lenv = length(vv)

    @floop executor for i in 1:(lenu*lenv)
        ix = Int(i%lenv+1)
        iy = Int(floor((i-1)/lenu))+1
        vis[ix, iy] = conj(vis[(ix, iy)...])*cispi(2*(uu[ix]*x0 + vv[iy]*y0))
    end
    return vis
end

and then allow the user to swap out the SequentialEx() with a ThreadedEx() if they want.

• Rework data format and data layout to make it easier to handle multifrequency and time-dependent models
• Change the instrument model format so you only specify the cache once and allow for easier tracking of different instrument terms
• rename visibilities to visilibilitymap and embrace the grid system.

When saving fits images, flux appears to be bleeding between adjacent edges (e.g., the left and top edges).

NAXIS1 and NAXIS2 in saving fits seems to be swapped in https://github.com/ptiede/Comrade.jl/blob/32ecaf3aff907d35e21c26e66618ac8d91014486/src/observations/io.jl#LL68-L73

So even if FOV is correct in x and y, number of pixels are most likely exchanged for x and y.

Hi @ptiede , I was going through Comrade's tutorial and found that the latest Comrade with Julia 1.9.3 / 1.9.4 does not export either of "ResponseCache", "JonesModel", "jonescache", "station_tuple", and "jonesmap" defined in src/calibration/jones.jl. I was not sure if the current code does have a grammatical error but fixing the line 2 of jones.jl in the following way seems to fix the issue.

export JonesCache, TrackSeg, ScanSeg, FixedSeg, IntegSeg, jonesG, jonesD, jonesT
export ResponseCache, JonesModel, jonescache, station_tuple, jonesmap

FYI, everything exported in the first line is loaded correctly in the current version. Please check if this fix will help to mitigate the issue as these are critical for polarization modeling.

The current sample method can save checkpoints with scans of chains and tuned HMC sampler configurations in disk, but doesn't support resuming from a checkpoint. The application of comrade HMC samplers will get wider if it has an option to resume a previous run from a checkpoint, as this allows a long run often chucked in multiple MCMC runs.

https://github.com/ptiede/Comrade.jl/blob/57c70c913ea7ee98e1d438094a9fbf6aafb05745/src/models/modelimage/modelimage.jl#LL69-L80

The type of the "model" argument in the "modelimage" function differs between the docs and the source code.

For fitting different bands at the same time, we need to define MultiRadioLikelihood and use skymodel and vlbimodel for getting residuals and intensity map. Currently skymodel and vlbimodel accept only RadioLikelihood type and not MultiRadioLikelihood .

Hey @ptiede, happy to be reviewing this JOSS submission. I'll be using this issue for working discussions and problems I find!

I'll start with doc typos I've found; this is what I've caught while reading through-

Comrade attempts this imaging uncertainty

this sentence is missing a verb, I think?

Loading Data into Comrade

To install eht-imaging see the [ehtim](https://github.com/achael/ehtim) github repo.

dead link, should point to achael/eht-imaging

During this part of the first tutorial

vis = extract_vis(obs) #complex visibilites

I got some warnings

WARNING: both HypercubeTransform and MeasureTheory export "inverse"; uses of it in module Comrade must be qualified

that seem like they can be addressed.

Making Image of a Black Hole

Very minor comment: for this tutorial, I think it would be clearer if you put the package imports where they are used- for one, if the user isn't interested in, say, the AHMC part of the tutorial, then they won't need to install/load the package. Second, it makes it super obvious that "I need this to enable these things", like you don't need Distributions.jl until you create the priors, and you don't need GalacticOptim.jl and ComradeGalactic.jl until you find the MAP estimate.

(Tangential note: There's something flagrantly un-Julian about having to make 5 unique packages for 5 unique optimization algorithms when the language is build on composability. I don't have a solution, and this is fairly commonplace, but it's obviously a lot of extra work for you as a maintainer and is annoying at minimum as a user)

Problems:

• #142
• ComradeNested and ComradeAdaptMCMC can't be used at the same time due to disparate AbstractMCMC.jl compat (ironically I hold the power to fix this- will look at it soon!)

julia> prob = OptimizationProblem(f, randn(ndim), nothing, lb=fill(-5.0, ndim), ub=fill(5.0, ndim))
ERROR: MethodError: no method matching OptimizationProblem[...]

Check out # from EHTC VI 2019. for some ideas on how to make this better!

not clear what "# from EHTC VI 2019" is supposed to mean, and extraneous period.

Modeling with non-analytic Fourier transforms

m = ExtendedRing(10.0, 8.0)

seems to be invalid syntax, should be

m = ExtendedRing((10.0, 8.0))

julia> plot(m, xlims=(-40.0, 40.0), ylims=(-40.0, 40.0), uvscale=identity)
ERROR: MethodError: no method matching +(::Tuple{Float64, Float64}, ::Int64)

julia> mimage = modelimage(m, image; alg=Comrade.FFT())
ERROR: UndefVarError: FFT not defined

seems like this part of the tutorial is out of date.

Taking a break for now, will keep adding as I continue reviewing the code!

Currently for multi-frequency we just just make separate likelihoods for each frequency however we could instead bake this directly into our model. This is probably required for more complete instrument modeling and better multifrequency models.

Note: It's been a few weeks since I last ran into this issue, so it's possible that there was a patch pushed somewhere upstream that fixed it, though I don't see anything in recent Comrade commits that might have affected it. I'm leaving this here for posterity in case it's helpful to someone, or in case there's a regression. Feel free to close if this is no longer an issue. Update: it happened again today, when I grabbed and ran a fresh copy of the code in this repository: https://github.com/saopeter/challenge1_comrade

I have previously encountered the following error when trying to run https://github.com/ptiede/Comrade.jl/blob/main/examples/imaging.jl directly from a terminal:

** (process:95014): WARNING **: 14:01:35.389: Failed to load shared library 'libpango-1.0.so.0' referenced by the typelib: /usr/lib/libpango-1.0.so.0: undefined symbol: hb_ot_layout_get_horizontal_baseline_tag_for_script
ERROR: PyError (PyImport_ImportModule) <class 'AssertionError'>

This seems to be related to the issue mentioned here: JuliaPlots/Plots.jl#4282

A possible workaround is to first start a Julia REPL and manually type in the following commands:

using Comrade
load_ehtim()
include("imaging.jl")

then the imaging.jl script will run as expected.

There are a number of type instabilities when creating the NFFT cache. These need to be fixed.

Currently we use TypedTables.jl to store our chains and other stuff. While this is a simple table serializing it and making sure it is stable between versions has been annoying. InferenceObjects.jl solves a lot of this. Plus, we can then use Arviz.jl to do all the MCMC post-processing stuff, which is great.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

The ascube and asflat posterior transformations are currently not gpu friendly. These transformations could potentially be made gpu friendly and non allocating. Here's a series of helpful places to look:

• TransformationVariables.jl
  TransformationVariables is the underlying dependency used to handle the posterior transformations.
• ComponentVectors.jl
  The current representation of the data is in terms of named tuples. Representations in terms of componentvectors are likely to be more GPU friendly

A change in commit cdbcc6b causes intensitymap to fail when called on models that are defined with caches. The change is on around line 86 of of src/models/modelimage/nfft_alg.jl where,

@inline function create_cache(alg::ObservedNUFT{<:NFFTAlg}, plan, phases, img)
    timg = IntensityMap(transpose(img.im), img.fovx, img.fovy, img.pulse)
    return NUFTCache(alg, plan, phases, img.pulse, timg)

was changed to

@inline function create_cache(alg::ObservedNUFT{<:NFFTAlg}, plan, phases, img)
    #timg = #IntensityMap(transpose(img.im), img.fovx, img.fovy, img.pulse)
    return NUFTCache(alg, plan, phases, img.pulse, img.im')

In https://ptiede.github.io/Comrade.jl/dev/benchmarks/ : "In fact, Comrade took some design queue's from Themis." -> design cues

In https://ptiede.github.io/Comrade.jl/dev/libs/dynesty/ : "ComradDynesty" -> ComradeDynesty

Here is a list of things I need to do before releasing 0.7

• Add polarization tests (#229 )
• Add jones matrix and RIME tests
• Add missing data checks for coherency matrices (partially handled in VLBILikelihoods 0.1.2 but needs testing)
• Fix docs to use new formalism (#188)
• Add polarization example to the docs (#188)
• Fix all examples to use new formalism (#188)
• Fix all library packages to match new 0.7 formalism (#254 )
• Improve test coverage to a minimum of 85%
• Add conventions page
• Remove unnecessary dependencies (#253)
• Find slowdown (#261)
• Add relative or transitional gain model (#264)

Deserialization of Pigeons.jl checkpoints with NFFT plans causes a segmentation fault upon visibility evaluations. This can be fixed by dispatching off the Serialization.deserialize function and reconstructing the FFT plan.

using Serialization
using Comrade

function Serialization.deserialize(s::AbstractSerializer, type::Type{<:Comrade.NUFTCache{<:Comrade.ObservedNUFT,P,M,PI,I}}) where {P,M,PI,I}

    alg = Serialization.deserialize(s)
    _ = Serialization.deserialize(s)
    phase = Serialization.deserialize(s)
    pulse = Serialization.deserialize(s)
    img = Serialization.deserialize(s)

    newplan = Comrade.plan_nuft(alg, img)
    return Comrade.NUFTCache(alg, newplan, phase, pulse, img)
end

Currently we assume that all d-terms and gains have identical segmentations across all scans. This is usually reasonable, however sometimes the data is such that a single station has unique properties, i.e. time variable d-terms. We should probably partially refactor the Jones matrix stuff and cache stuff to be more station dependent. Eventually everything will get compiled down to the same sparse design matrix formalism before inference but this should allow for us to more easily switch construct complicated models.

Hi @ptiede, thanks for your work on Comrade.jl - I enjoyed testing out the examples and the documentation is very thorough. I have a few minor comments that I hope can be useful:

• I think the paper/docs could benefit from a brief explanation of the likelihood used to connect the prior model to the EHT data. How is this done? What assumptions are made? How does uncertainty on the data/model side come into this? Maybe it is a standard procedure in the field, in which case a reference can be added.

• In a general setting, the user should check the diagnostics of whichever tool they use for optimization/posterior sampling. While this is package-dependent and not part of your code, I think adding some extra steps to the tutorial would make sense to encourage such a workflow. E.g. for AHMC in the black hole image tutorial, the effective sample size and rhat are standard diagnostics.

• I ran the black hole image example several times and once ran into the following error when calling

  plot(sqrt.(meanimg), title="Mean Image")
  ERROR: LoadError: DomainError with -5.97700199247224e-20:
  sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

• I wanted to check out the performance claims in the docs (Benchmarks) but I didn't manage to get the code running. The code snippet in the docs fails with

  ERROR: LoadError: MethodError: no method matching transform(::TransformVariables.TransformTuple{NamedTuple{(:rad, :wid, :a, :b, :f, :sig, :asy, :pa, :x, :y), NTuple{10, TransformVariables.ScaledShiftedLogistic{Float64}}}}, ::NamedTuple{(:f, :rad, :wid, :a, :b, :sig, :asy, :pa, :x, :y), NTuple{10, Float64}})

  I also found some recently updated code in examples/benchmark.jl in the GitHub, but I didn't manage to get that running either:

  ERROR: LoadError: ArgumentError: keys(transformations) == keys(y) must hold. Got
  keys(transformations) => (:rad, :wid, :a, :b, :f, :sig, :asy, :pa, :x, :y)
  keys(y) => (:f, :rad, :wid, :a, :b, :sig, :asy, :pa, :x, :y)

• Does the approach used in this paper with this software also merit a mention in the "similar packages" section of the paper?

• I notice that you have automated tests set up but the CI is currently failing and the coverage is currently 0%, although in usually around 75%. Is this a temporary bug with the workflows/uploading?

I need to fix Comrade._mring_vis to work with Tuples. Currently, this can fail to accumulate properly.
This means I probably need to make the pullback actually non-allocating.

#180 got into trouble because the docs aren't being generated correctly for old releases. This needs to be fixed.

Enzyme is great and from some testing it looks like the most recent versions will work with most if not all the geometric models we currently use. Now because enzyme doesn't fully have its own rule system I can't (yet) fully switch over but I can add special rrules for geometric models to use Enzyme. This should happen soon. Especially for geometric polarized fits where we probably want to have faster gradients.

When running the example in https://ptiede.github.io/Comrade.jl/dev/examples/data/#Loading-Data-into-Comrade, specifically

obseht = Pyehtim.load_uvfits_and_array(
                joinpath(dirname(pathof(Comrade)), "..", "examples", "PolarizedExamples/polarized_gaussian_all_corruptions.uvfits"),
                joinpath(dirname(pathof(Comrade)), "..", "examples", "PolarizedExamples/array.txt"),
                polrep="circ"
                        )

gives an error message

Loading uvfits:  /home/kjwiik/.julia/packages/Comrade/G7U3R/src/../examples/PolarizedExamples/polarized_gaussian_all_corruptions.uvfits
ERROR: Python: FileNotFoundError: [Errno 2] No such file or directory: '/home/kjwiik/.julia/packages/Comrade/G7U3R/src/../examples/PolarizedExamples/polarized_gaussian_all_corruptions.uvfits'

https://github.com/ptiede/Comrade.jl/blob/ab8d15cbd210537503a210adedb416bffe9a2ccd/src/bayes/bayes.jl#LL16-L37

DensityInterface points to https://github.com/JuliaMath/DensityInferface.jl instead of https://github.com/JuliaMath/DensityInterface.jl

Wishlist

Alright so Comrade in its current forms works (kind of) but there are a number of things
I need to do and check. I'll split this into Urgent Moderate Moonshot

Urgent

• Add gains
• Add docstings to all functions
• Create documentation and some examples
• Figure out how to do NFFT caching for numerical visibility functions (note this may require a huge refactor due to poor initial design choices)
• Unit test all functions
• Interface with eht-imaging (I need to be able to read in uvfits files)

Moderate

• I really need to work out some performance and autodiff functionality
• Create movie models intrinsically? I probably want some helper functions
• Polarization and RIME primitives (this will take some thought)

Moonshot

• eht specific inference packages and interfacing with all the great ray-tracing tools out there
• Create a autodiffable ray-tracing and radiative transfer module compatible with GPU's although this should really be a separate package (I would like to have someone more knowledgeable help me with this)
• MPI/Dagger/parallel stuff/tools for people to use. Right now I don't have this.

It should be

RadioLikelihood(skymodel, dataproducts::EHTObservation...; skymeta=nothing)

instead of

RadioLikelihood(skymodel, obs, dataproducts::EHTObservation...; skymeta=nothing)

I was trying out the code given here : https://ptiede.github.io/Comrade.jl/stable/examples/imaging_closures/
But the line obs = scan_average(obs).add_fractional_noise(0.01).flag_uvdist(uv_min=0.1e9) gives me the following error.

ERROR: Python: TypeError: must be real number, not complex
Python stacktrace:
 [1] pandas._libs.lib.maybe_convert_objects
   @ pandas/_libs/lib.pyx:2466
 [2] agg_series
   @ pandas.core.groupby.ops ~/.local/lib/python3.10/site-packages/pandas/core/groupby/ops.py:996
 [3] _python_agg_general
   @ pandas.core.groupby.generic ~/.local/lib/python3.10/site-packages/pandas/core/groupby/generic.py:288
 [4] aggregate
   @ pandas.core.groupby.generic ~/.local/lib/python3.10/site-packages/pandas/core/groupby/generic.py:266
 [5] <dictcomp>
   @ pandas.core.apply ~/.local/lib/python3.10/site-packages/pandas/core/apply.py:421
 [6] agg_dict_like
   @ pandas.core.apply ~/.local/lib/python3.10/site-packages/pandas/core/apply.py:420
 [7] agg
   @ pandas.core.apply ~/.local/lib/python3.10/site-packages/pandas/core/apply.py:163
 [8] aggregate
   @ pandas.core.groupby.generic ~/.local/lib/python3.10/site-packages/pandas/core/groupby/generic.py:1269
 [9] coh_avg_vis
   @ ehtim.statistics.dataframes ~/Code/Julia/Projects/cosmology/.CondaPkg/env/lib/python3.10/site-packages/ehtim/statistics/dataframes.py:198
 [10] avg_coherent
   @ ehtim.obsdata ~/Code/Julia/Projects/cosmology/.CondaPkg/env/lib/python3.10/site-packages/ehtim/obsdata.py:1143
Stacktrace:
 [1] pythrow()
   @ PythonCall ~/.julia/packages/PythonCall/qTEA1/src/err.jl:94
 [2] errcheck
   @ ~/.julia/packages/PythonCall/qTEA1/src/err.jl:10 [inlined]
 [3] pycallargs
   @ ~/.julia/packages/PythonCall/qTEA1/src/abstract/object.jl:211 [inlined]
 [4] pycall(f::Py, args::Float64; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:scan_avg,), Tuple{Bool}}})
   @ PythonCall ~/.julia/packages/PythonCall/qTEA1/src/abstract/object.jl:222
 [5] pycall
   @ ~/.julia/packages/PythonCall/qTEA1/src/abstract/object.jl:218 [inlined]
 [6] #_#11
   @ ~/.julia/packages/PythonCall/qTEA1/src/Py.jl:341 [inlined]
 [7] Py
   @ ~/.julia/packages/PythonCall/qTEA1/src/Py.jl:341 [inlined]
 [8] scan_average(obs::Py)
   @ Pyehtim ~/.julia/packages/Pyehtim/YyZzM/src/Pyehtim.jl:63
 [9] top-level scope
   @ REPL[15]:1

Information on the packages:
[deps]
Comrade = "99d987ce-9a1e-4df8-bc0b-1ea019aa547b"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Pyehtim = "3d61700d-6e5b-419a-8e22-9c066cf00468"
and using julia 1.9.2