ptiede / comrade.jl Goto Github PK
View Code? Open in Web Editor NEWLicense: MIT License
License: MIT License
When saving fits images, flux appears to be bleeding between adjacent edges (e.g., the left and top edges).
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
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.
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.
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.
The type of the "model" argument in the "modelimage" function differs between the docs and the source code.
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.
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
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.
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.
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?
Here is a list of things I need to do before releasing 0.7
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!
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.
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
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
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))
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')
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.
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'
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
#180 got into trouble because the docs aren't being generated correctly for old releases. This needs to be fixed.
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, 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.
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?
To install eht-imaging see the [ehtim](https://github.com/achael/ehtim) github repo.
dead link, should point to achael/eht-imaging
Comrade.jl/docs/src/examples/data.jl
Line 17 in 79b66f8
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.
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:
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.
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!
Comrade.jl/src/models/modelimage/fft_alg.jl
Lines 150 to 154 in 252a073
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.
julia>] activate --temp
(jl_z15mWM) pkg> add ComradeGalactic
[...]
julia> using ComradeGalactic
ERROR: LoadError: ArgumentError: Package ComradeGalactic does not have LinearAlgebra in its dependencies:
[...]
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:
It should be
RadioLikelihood(skymodel, dataproducts::EHTObservation...; skymeta=nothing)
instead of
RadioLikelihood(skymodel, obs, dataproducts::EHTObservation...; skymeta=nothing)
There are a number of type instabilities when creating the NFFT cache. These need to be fixed.
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.