update docstrings and examples

Author: mileslucas

docs/src/api.md

@@ -20,7 +20,7 @@ PSFModels.Normal
 ```
 
 ```@example plots
-model = Gaussian(10)
+model = Gaussian(fwhm=10)
 plot(model; title="Gaussian(fwhm=10)")
 ```
 
@@ -31,7 +31,7 @@ PSFModels.AiryDisk
 ```
 
 ```@example plots
-model = AiryDisk(10)
+model = AiryDisk(fwhm=10)
 plot(model; title="AiryDisk(fwhm=10)")
 ```
 
@@ -42,6 +42,6 @@ PSFModels.Moffat
 ```
 
 ```@example plots
-model = Moffat(10)
+model = Moffat(fwhm=10)
 plot(model; title="Moffat(fwhm=10)")
 ```

docs/src/examples.md

@@ -32,7 +32,7 @@ function model(X::AbstractVector{T}) where T
     position = @view X[1:2] # x, y position
     fwhm     = @view X[3:4] # fwhm_x, fwhm_y
     amp      =       X[5]   # amplitude
-    return amp * Gaussian{T}(position, fwhm)
+    return Gaussian(T, pos=position, fwhm=fwhm, amp=amp)
 end
 
 # objective function
@@ -41,13 +41,13 @@ function loss(X::AbstractVector{T}, target) where T
     all(>(0), X) || return T(Inf)
     # get generative model
     m = model(X)
-    # l2-distance loss (χ² loss) (LossFunctions.jl)
+    # l2-distance loss (χ² loss) (LossFunctions.jl) without cutting out
     stamp = @view m[axes(target)...]
     return value(L2DistLoss(), target, stamp, AggMode.Sum())
 end
 
 # params are [x, y, fwhm_x, fwhm_y, amp]
-test_params = Float32[20, 20, 5, 5, 1]
+test_params = eltype(psf)[20, 20, 5, 5, 1]
 loss(test_params, psf)
 ```
 
@@ -79,7 +79,7 @@ synth_psf = model(best_fit_params)
 
 plot(
     imshow(psf, title="Data"),
-    plot(synth_psf, axes(psf); title="Model"),
+    plot(synth_psf, reverse(axes(psf)); title="Model"),
     cbar=false,
     ticks=false,
     layout=2,

docs/src/index.md

@@ -28,12 +28,12 @@ To import the library
 julia> using PSFModels
 ```
 
-None of the models are exported to avoid namespace clashes, but it can be verbose. You can either import names directly
+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
 
 ```julia
 julia> using PSFModels: Gaussian
 
-julia> model = Gaussian(8)
+julia> model = Gaussian(fwhm=8)
 ```
 
 or you can create an alias for `PSFModels`
@@ -43,9 +43,9 @@ or you can create an alias for `PSFModels`
 using PSFModels
 const M = PSFModels
 # julia version 1.6 or above
-using PSFModels as M
+import PSFModels as M
 
-model = M.Gaussian(10)
+model = M.Gaussian(fwhm=10)
 ```
 
 ```@docs

src/PSFModels.jl

@@ -6,20 +6,21 @@ Statistical models for constructing point-spread functions (PSFs). These models
 ## Models
 
 The following models are currently implemented
-* [`PSFModels.Gaussian`](@ref)
+* [`PSFModels.Gaussian`](@ref)/[`PSFModels.Normal`](@ref)
 * [`PSFModels.AiryDisk`](@ref)
 * [`PSFModels.Moffat`](@ref)
 
 ## Parameters
 
-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).
+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).
 
 ## Usage
 
-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.
+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.
 
 ```jldoctest model
-julia> m = PSFModels.Gaussian(5); # fwhm of 5 pixels, centered at (0, 0)
+julia> m = PSFModels.Gaussian(fwhm=3); # centered at (0, 0)
+
 
 julia> m[0, 0] # [y, x] for indexing
 1.0
@@ -32,20 +33,21 @@ julia> m(0, 0) # (x, y) for evaluating
 
     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)`).
 
-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`
+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`
 
 ```jldoctest model
-julia> size(m)
-(17, 17)
+julia> size(m) # default 'stamp' size is fwhm * 3
+(11, 11)
 
 julia> axes(m)
-(-8:8, -8:8)
+(-5:5, -5:5)
 ```
 
 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
 
 ```jldoctest model
 julia> stamp = collect(m);
+
 ```
 
 these axes are merely a convenience for bounding the model, since they accept any real number as input.
@@ -55,21 +57,22 @@ julia> m[100, 10000] # index-like inputs [y, x]
 0.0
 
 julia> m(2.4, 1.7) # valid for any real (x, y)
-0.38315499005194587
+0.0696156536973086
 ```
 
 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
 
 ```jldoctest
-julia> big_mat = ones(101, 101);
+julia> big_mat = ones(1001, 1001);
+
+julia> model = PSFModels.Gaussian(x=51, y=51, fwhm=2);
 
-julia> model = PSFModels.Gaussian(51, 51, 2); # center of big_mat, fwhm=2
 
 julia> ax = map(intersect, axes(big_mat), axes(model))
 (48:54, 48:54)
 
 julia> cutout = @view big_mat[ax...]
-7×7 view(::Array{Float64,2}, 48:54, 48:54) with eltype Float64:
+7×7 view(::Matrix{Float64}, 48:54, 48:54) with eltype Float64:
  1.0  1.0  1.0  1.0  1.0  1.0  1.0
  1.0  1.0  1.0  1.0  1.0  1.0  1.0
  1.0  1.0  1.0  1.0  1.0  1.0  1.0
@@ -80,18 +83,20 @@ julia> cutout = @view big_mat[ax...]
 
 julia> stamp = @view model[ax...];
 
+
 julia> photsum = sum(cutout .* stamp)
 4.5322418212890625
 ```
-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).
+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).
 
 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).
 
 ```jldoctest
-julia> mbig = PSFModels.Gaussian{BigFloat}(12);
+julia> mbig = PSFModels.Gaussian(BigFloat, fwhm=12);
+
 
 julia> sum(mbig)
-163.07467408408593790971336380361822460116627553361468017101287841796875
+163.07467408408593885562554918859656805096847165259532630443572998046875
 ```
 
 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.
@@ -100,16 +105,9 @@ finally, we provide plotting recipes from [RecipesBase.jl](https://github.com/Ju
 using Plots
 model = PSFModels.Gaussian(8)
 plot(model)              # default axes
-plot(model, 1:5, 1:5)    # custom axes
+plot(model, :, 1:5)    # custom axes (y, x)
 plot(model, axes(other)) # use axes from other array
 ```
-
-!!! tip "Tip: Automatic Differentation"
-    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
-    ```julia
-    # example generative model for position and scalar fwhm
-    model(X::AbstractVector{T}) where {T} = PSFModels.Gaussian{T}(X...)
-    ```
 """
 module PSFModels