Skip to content

Commit 69340fe

Browse files
authored
Merge pull request #2 from JuliaAstro/ml/rewrite
rewrite base model interface
2 parents dcb743f + 50208f0 commit 69340fe

File tree

17 files changed

+285
-287
lines changed

17 files changed

+285
-287
lines changed

.github/workflows/ci.yml

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ jobs:
1313
fail-fast: false
1414
matrix:
1515
version:
16-
- '1.3'
16+
- '1.5'
1717
- '1'
1818
- 'nightly'
1919
os:
@@ -44,6 +44,17 @@ jobs:
4444
- uses: codecov/codecov-action@v1
4545
with:
4646
file: lcov.info
47+
aqua:
48+
name: Aqua tests
49+
runs-on: ubuntu-latest
50+
steps:
51+
- uses: actions/checkout@v2
52+
- uses: julia-actions/setup-julia@v1
53+
with:
54+
version: '1'
55+
- uses: julia-actions/julia-buildpkg@v1
56+
- run: julia -e 'using Pkg; Pkg.add("Aqua")'
57+
- run: julia --project=@. -e 'using Aqua, PSFModels; Aqua.test_all(PSFModels, ambiguities=false)'
4758
docs:
4859
name: Documentation
4960
runs-on: ubuntu-latest

.github/workflows/preview-cleanup.yml

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
name: Doc Preview Cleanup
2+
3+
on:
4+
pull_request:
5+
types: [closed]
6+
7+
jobs:
8+
doc-preview-cleanup:
9+
runs-on: ubuntu-latest
10+
steps:
11+
- name: Checkout gh-pages branch
12+
uses: actions/checkout@v2
13+
with:
14+
ref: gh-pages
15+
16+
- name: Delete preview and history
17+
run: |
18+
git config user.name "Documenter.jl"
19+
git config user.email "[email protected]"
20+
git rm -rf "previews/PR$PRNUM"
21+
git commit -m "delete preview"
22+
git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree})
23+
env:
24+
PRNUM: ${{ github.event.number }}
25+
26+
- name: Push changes
27+
run: |
28+
git push --force origin gh-pages-new:gh-pages

Project.toml

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,16 @@ version = "0.2.0"
66
[deps]
77
CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
88
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
9+
KeywordCalls = "4d827475-d3e4-43d6-abe3-9688362ede9f"
910
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
1011
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
1112
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1213

1314
[compat]
1415
CoordinateTransformations = "0.6"
1516
Distances = "0.10"
17+
KeywordCalls = "0.2"
1618
RecipesBase = "1"
1719
SpecialFunctions = "0.10, 1"
1820
StaticArrays = "0.12, 1"
19-
julia = "1.3"
21+
julia = "1.5"

README.md

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -81,10 +81,6 @@ m([1.2, 0.4])
8181
It is important to recognize the difference in the order of the dimensions between indexing and calling. Indexing is reverse of the cartesian order, which is the natural way of indexing a multi-dimensional array. If you try calling a model as a function with an index, an error will be thrown.
8282

8383
```julia
84-
# scalar multiplication or division will create a ScaledPSFModel
85-
20 * m # or `m * 20`
86-
m / 20
87-
8884
# evaluate `m` over its indices forming an array
8985
collect(m)
9086

@@ -97,6 +93,5 @@ m .* arr
9793
inds = map(intersect, axes(arr), axes(m))
9894
arr_stamp = @view arr[inds...]
9995
m_stamp = @view m[inds...]
100-
amp = 1.24
10196
resid = sum(abs2, arr_stamp .- amp .* m_stamp) # chi-square loss
10297
```

docs/src/api.md

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,6 @@ using Plots
1010

1111
```@docs
1212
PSFModels.PSFModel
13-
PSFModels.ScaledPSFModel
1413
```
1514

1615
## Gaussian
@@ -21,7 +20,7 @@ PSFModels.Normal
2120
```
2221

2322
```@example plots
24-
model = Gaussian(10)
23+
model = Gaussian(fwhm=10)
2524
plot(model; title="Gaussian(fwhm=10)")
2625
```
2726

@@ -32,7 +31,7 @@ PSFModels.AiryDisk
3231
```
3332

3433
```@example plots
35-
model = AiryDisk(10)
34+
model = AiryDisk(fwhm=10)
3635
plot(model; title="AiryDisk(fwhm=10)")
3736
```
3837

@@ -43,6 +42,6 @@ PSFModels.Moffat
4342
```
4443

4544
```@example plots
46-
model = Moffat(10)
45+
model = Moffat(fwhm=10)
4746
plot(model; title="Moffat(fwhm=10)")
4847
```

docs/src/examples.md

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -29,10 +29,11 @@ using LossFunctions
2929
3030
# generative model
3131
function model(X::AbstractVector{T}) where T
32-
position = @view X[1:2] # x, y position
33-
fwhm = @view X[3:4] # fwhm_x, fwhm_y
34-
amp = X[5] # amplitude
35-
return amp * Gaussian{T}(position, fwhm)
32+
x = X[1] # position
33+
y = X[2]
34+
fwhm = @view X[3:4] # fwhm_x, fwhm_y
35+
amp = X[5] # amplitude
36+
return Gaussian(T; x, y, fwhm, amp)
3637
end
3738
3839
# objective function
@@ -41,13 +42,13 @@ function loss(X::AbstractVector{T}, target) where T
4142
all(>(0), X) || return T(Inf)
4243
# get generative model
4344
m = model(X)
44-
# l2-distance loss (χ² loss) (LossFunctions.jl)
45+
# l2-distance loss (χ² loss) (LossFunctions.jl) without cutting out
4546
stamp = @view m[axes(target)...]
4647
return value(L2DistLoss(), target, stamp, AggMode.Sum())
4748
end
4849
4950
# params are [x, y, fwhm_x, fwhm_y, amp]
50-
test_params = Float32[20, 20, 5, 5, 1]
51+
test_params = eltype(psf)[20, 20, 5, 5, 1]
5152
loss(test_params, psf)
5253
```
5354

@@ -79,7 +80,7 @@ synth_psf = model(best_fit_params)
7980
8081
plot(
8182
imshow(psf, title="Data"),
82-
plot(synth_psf, axes(psf); title="Model"),
83+
plot(synth_psf, reverse(axes(psf)); title="Model"),
8384
cbar=false,
8485
ticks=false,
8586
layout=2,

docs/src/index.md

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,12 +28,12 @@ To import the library
2828
julia> using PSFModels
2929
```
3030

31-
None of the models are exported to avoid namespace clashes, but it can be verbose. You can either import names directly
31+
None of the models are exported to avoid namespace clashes, but it can be verbose to continuously rewrite `PSFModels`. You can either import names directly
3232

3333
```julia
3434
julia> using PSFModels: Gaussian
3535

36-
julia> model = Gaussian(8)
36+
julia> model = Gaussian(fwhm=8)
3737
```
3838

3939
or you can create an alias for `PSFModels`
@@ -43,9 +43,9 @@ or you can create an alias for `PSFModels`
4343
using PSFModels
4444
const M = PSFModels
4545
# julia version 1.6 or above
46-
using PSFModels as M
46+
import PSFModels as M
4747

48-
model = M.Gaussian(10)
48+
model = M.Gaussian(fwhm=10)
4949
```
5050

5151
```@docs

src/PSFModels.jl

Lines changed: 23 additions & 92 deletions
Original file line numberDiff line numberDiff line change
@@ -6,20 +6,21 @@ Statistical models for constructing point-spread functions (PSFs). These models
66
## Models
77
88
The following models are currently implemented
9-
* [`PSFModels.Gaussian`](@ref)
9+
* [`PSFModels.Gaussian`](@ref)/[`PSFModels.Normal`](@ref)
1010
* [`PSFModels.AiryDisk`](@ref)
1111
* [`PSFModels.Moffat`](@ref)
1212
1313
## Parameters
1414
15-
In general, the PSFs have a position, a full-width at half-maximum (FWHM) measure, and an amplitude. The position a 1-based pixel coordinate system, where `(1, 1)` represents the *center* of the bottom left pixel. This matches the indexing style of Julia as well as DS9, IRAF, SourceExtractor, and WCS. If a position is not specified, it is set to `(0, 0)`. The FWHM is a consistent scale parameter for the models. All models support a scalar (isotropic) FWHM and a FWHM for each axis (diagonal).
15+
In general, the PSFs have a position, a full-width at half-maximum (FWHM) measure, and an amplitude. The position follows a 1-based pixel coordinate system, where `(1, 1)` represents the *center* of the bottom left pixel. This matches the indexing style of Julia as well as DS9, IRAF, SourceExtractor, and WCS. If a position is not specified, it is set to `(0, 0)`. The FWHM is a consistent scale parameter for the models. All models support a scalar (isotropic) FWHM and a FWHM for each axis (diagonal).
1616
1717
## Usage
1818
19-
Using the models should feel just like an array. In fact, `PSFModels.PSFModel <: AbstractMatrix`. However, no data is stored and no allocations have to be made. In other words, representing the models as matrices is merely a convenience, since typically astronomical data is stored in dense arrays.
19+
Using the models should feel just like an array. In fact, `PSFModels.PSFModel <: AbstractMatrix`. However, no data is stored and no allocations have to be made. In other words, representing the models as matrices is merely a convenience, since typically astronomical data is stored in dense arrays. Another way of thinking of these is a lazy array that applies a function when indexed, rather than returning data stored in memory.
2020
2121
```jldoctest model
22-
julia> m = PSFModels.Gaussian(5); # fwhm of 5 pixels, centered at (0, 0)
22+
julia> m = PSFModels.Gaussian(fwhm=3); # centered at (0, 0)
23+
2324
2425
julia> m[0, 0] # [y, x] for indexing
2526
1.0
@@ -32,34 +33,21 @@ julia> m(0, 0) # (x, y) for evaluating
3233
3334
It's important to note the difference in the axis ordering between the index-style calls and the function-style calls. The index-style calls are reverse cartesian order (e.g., `(z, y, x)`), while function calls are the typical cartesian order `(x, y, z)`. Regardless, the constructors are always in cartesian order (`(x, y, z)`).
3435
35-
To control the amplitude, the best method is using scalar multiplication or division. These operations create another lazy object ([`ScaledPSFModel`](@ref)) that scales the original model without having to broadcast and potentially allocate.
36+
Because the model is a matrix, it needs to have a size. Each model has a bounding box which can be controlled with the `extent` keyword. By default the extent is set by a scalar factor of the FWHM (e.g., `maxsize * FWHM` pixels), centered around the PSF, and rounded up. We can see how this alters the indices from a typical `Matrix`
3637
3738
```jldoctest model
38-
julia> m_scaled = 20 * m;
39-
40-
julia> m_scaled(0, 0)
41-
20.0
42-
43-
julia> m′ = m_scaled / 20;
44-
45-
julia> m′(0, 0)
46-
1.0
47-
```
48-
49-
Because the model is a matrix, it needs to have a size. In this case, the size is `maxsize * FWHM` pixels, centered around the origin, and rounded up. We can see how this alters the indices from a typical `Matrix`
50-
51-
```jldoctest model
52-
julia> size(m)
53-
(17, 17)
39+
julia> size(m) # default 'stamp' size is fwhm * 3
40+
(11, 11)
5441
5542
julia> axes(m)
56-
(-8:8, -8:8)
43+
(-5:5, -5:5)
5744
```
5845
5946
if we want to collect the model into a dense matrix, regardless of the indexing (e.g. to prepare for cross-correlation), we can simply
6047
6148
```jldoctest model
6249
julia> stamp = collect(m);
50+
6351
```
6452
6553
these axes are merely a convenience for bounding the model, since they accept any real number as input.
@@ -69,21 +57,22 @@ julia> m[100, 10000] # index-like inputs [y, x]
6957
0.0
7058
7159
julia> m(2.4, 1.7) # valid for any real (x, y)
72-
0.38315499005194587
60+
0.0696156536973086
7361
```
7462
7563
By bounding the model, we get a cutout which can be applied to arrays with much larger dimensions without having to iterate over the whole matrix
7664
7765
```jldoctest
78-
julia> big_mat = ones(101, 101);
66+
julia> big_mat = ones(1001, 1001);
67+
68+
julia> model = PSFModels.Gaussian(x=51, y=51, fwhm=2);
7969
80-
julia> model = PSFModels.Gaussian(51, 51, 2); # center of big_mat, fwhm=2
8170
8271
julia> ax = map(intersect, axes(big_mat), axes(model))
8372
(48:54, 48:54)
8473
8574
julia> cutout = @view big_mat[ax...]
86-
7×7 view(::Array{Float64,2}, 48:54, 48:54) with eltype Float64:
75+
7×7 view(::Matrix{Float64}, 48:54, 48:54) with eltype Float64:
8776
1.0 1.0 1.0 1.0 1.0 1.0 1.0
8877
1.0 1.0 1.0 1.0 1.0 1.0 1.0
8978
1.0 1.0 1.0 1.0 1.0 1.0 1.0
@@ -94,18 +83,20 @@ julia> cutout = @view big_mat[ax...]
9483
9584
julia> stamp = @view model[ax...];
9685
86+
9787
julia> photsum = sum(cutout .* stamp)
9888
4.5322418212890625
9989
```
100-
Nice- we only had to reduce ~50 pixels instead of ~10,000 to calculate the aperture sum, all in under a microsecond (on my machine).
90+
Nice- we only had to reduce ~50 pixels instead of ~1,000,000 to calculate the aperture sum, all in under a microsecond (on my machine).
10191
10292
Since the models are lazy, that means the type of the output can be specified, as long as it can be converted to from a real number (so no integer types).
10393
10494
```jldoctest
105-
julia> mbig = PSFModels.Gaussian{BigFloat}(12);
95+
julia> mbig = PSFModels.Gaussian(BigFloat, fwhm=12);
96+
10697
10798
julia> sum(mbig)
108-
163.07467408408593790971336380361822460116627553361468017101287841796875
99+
163.07467408408593885562554918859656805096847165259532630443572998046875
109100
```
110101
111102
finally, we provide plotting recipes from [RecipesBase.jl](https://github.com/JuliaPlots/RecipesBase.jl), which can be seen in use in the [API/Reference](@ref) section.
@@ -114,82 +105,22 @@ finally, we provide plotting recipes from [RecipesBase.jl](https://github.com/Ju
114105
using Plots
115106
model = PSFModels.Gaussian(8)
116107
plot(model) # default axes
117-
plot(model, 1:5, 1:5) # custom axes
108+
plot(model, :, 1:5) # custom axes (y, x)
118109
plot(model, axes(other)) # use axes from other array
119110
```
120-
121-
!!! tip "Tip: Automatic Differentation"
122-
Forward-mode AD libraries tend to use dual numbers, which can cause headaches getting the types correct. We recommend using the *primal vector*'s element type to avoid these headaches
123-
```julia
124-
# example generative model for position and scalar fwhm
125-
model(X::AbstractVector{T}) where {T} = PSFModels.Gaussian{T}(X...)
126-
```
127111
"""
128112
module PSFModels
129113

130114
using CoordinateTransformations
131115
using Distances
116+
using KeywordCalls
132117
using SpecialFunctions
133118
using StaticArrays
134119

135-
"""
136-
PSFModels.PSFModel{T} <: AbstractMatrix{T}
137-
138-
Abstract type for PSF models.
139-
140-
In general, all `PSFModel`s have a set of pre-determined axes (the size is set upon creation) but they are lazy. That is, no memory is allocated and the values are calculated on the fly.
141-
142-
# Interface
143-
The interface to define a model is as follows (for an example model `Model`)
144-
145-
| method | description |
146-
|:-------|:------------|
147-
| `Model()` | constructor(s) |
148-
| `Base.size(m::Model)` | size, necessary for `AbstractArray` interface |
149-
| `Base.axes(m::Model)` | axes, necessary for `AbstractArray` interface |
150-
| `(m::Model)(point::AbstractVector)` | evaluate the model at the point in 2d space (x, y) |
151-
152-
browsing through the implementation of [`PSFModels.Gaussian`](@ref) should give a good idea of how to create a model
153-
"""
154-
abstract type PSFModel{T} <: AbstractMatrix{T} end
155-
156-
# always inbounds
157-
Base.checkbounds(::Type{Bool}, ::PSFModel, idx...) = true
158-
Base.checkbounds(::Type{Bool}, ::PSFModel, idx::CartesianIndex) = true
159-
160-
# in general, parse to static vector
161-
(model::PSFModel)(point...) = model(SVector(point))
162-
(model::PSFModel)(point::Tuple) = model(SVector(point))
163-
# disallow index to avoid confusion
164-
(model::PSFModel)(::CartesianIndex) = error("PSF models should be indexed using `getindex`, (equivalently `[]`)")
165-
166-
# getindex just calls model with reversed indices
167-
Base.getindex(model::PSFModel, idx::Vararg{<:Number,2}) = model(reverse(idx))
168-
169-
# broadcasting hack to slurp other axes (doesn't work for numbers)
170-
Broadcast.combine_axes(kern::PSFModel, other) = axes(other)
171-
Broadcast.combine_axes(other, kern::PSFModel) = axes(other)
172-
173-
indices_from_extent(pos, fwhm, maxsize) = indices_from_extent(pos, maxsize .* fwhm)
174-
175-
function indices_from_extent(pos, extent)
176-
halfextent = @. extent / 2
177-
lower = @. round(Int, pos - halfextent)
178-
upper = @. round(Int, pos + halfextent)
179-
return last(lower):last(upper), first(lower):first(upper)
180-
end
181-
182-
# TODO
183-
# function indices_from_extent(pos, fwhm::AbstractMatrix, maxsize)
184-
# halfextent = maxsize .* fwhm ./ 2
185-
# lower = @. floor(Int, pos - halfextent)
186-
# upper = @. ceil(Int, pos - halfextent)
187-
# end
188-
120+
include("core.jl")
189121
include("gaussian.jl")
190122
include("moffat.jl")
191123
include("airy.jl")
192-
include("scaled.jl")
193124
include("plotting.jl")
194125

195126
end # module PSFModels

0 commit comments

Comments
 (0)