text:
@userplot ImPlot
@recipe function f(h::ImPlot)
    if length(h.args) != 1 || !(typeof(h.args[1]) <: AbstractArray)
        error("Image plots require an arugment that is a subtype of AbstractArray.  Got: $(typeof(h.args))")
    end
    data = only(h.args)
    if !(typeof(data) <: AstroImage)
        data = AstroImage(only(h.args))
    end
    T = eltype(data)
    if ndims(data) != 2
        error("Image passed to `implot` must be two-dimensional.  Got ndims(img)=$(ndims(data))")
    end

    wcsn = get(plotattributes, :wcsn, 1)
    # Show WCS coordinates if wcsticks is true or unspecified, and has at least one WCS axis present.
    showwcsticks = get(plotattributes, :wcsticks, true) &&  !all(==(""), wcs(data, wcsn).ctype)
    showwcstitle = get(plotattributes, :wcstitle, true) &&  length(refdims(data)) > 0 && !all(==(""), wcs(data, wcsn).ctype)

    minx = first(parent(dims(data,1)))
    maxx = last(parent(dims(data,1)))
    miny = first(parent(dims(data,2)))
    maxy = last(parent(dims(data,2)))
    extent = (minx-0.5, maxx+0.5, miny-0.5, maxy+0.5)
    if haskey(plotattributes, :xlims)
        extent = (plotattributes[:xlims]..., extent[3:4]...)
    end
    if haskey(plotattributes, :ylims)
        extent = (extent[1:2]..., plotattributes[:ylims]...)
    end
    if showwcsticks
        wcsg = WCSGrid(data, Float64.(extent), wcsn)
        gridspec = wcsgridspec(wcsg)
    end

    # Use package defaults if not user provided.
    clims   --> _default_clims[]
    stretch --> _default_stretch[]
    cmap    --> _default_cmap[]

    bias = get(plotattributes, :bias, 0.5)
    contrast = get(plotattributes, :contrast, 1)
    platescale = get(plotattributes, :platescale, 1)

    grid := false
    # In most cases, a grid framestyle is a nicer looking default for images
    # but the user can override.
    framestyle --> :box

    if T <: Colorant
        imgv = data
    else
        clims   = plotattributes[:clims]
        stretch = plotattributes[:stretch]
        cmap    = plotattributes[:cmap]
        if T <: Complex
            img = abs.(data)
            img["UNIT"] = "magnitude"
        else
            img = data
        end
        imgv = imview(img; clims, stretch, cmap, contrast, bias)
    end

    # Reduce large images using the same heuristic as Images.jl
    maxpixels = get(plotattributes, :maxpixels, 10^6)
    _length1(A::AbstractArray) = length(eachindex(A))
    _length1(A) = length(A)
    while _length1(imgv) > maxpixels
        imgv = restrict(imgv)
    end

    # We have to do a lot of flipping to keep the orientation corect
    yflip := false
    xflip := false

    # Disable equal aspect ratios if the scales are totally different
    displayed_data_ratio = (extent[2]-extent[1])/(extent[4]-extent[3])
    if displayed_data_ratio >= 7
        aspect_ratio --> :none
    end

    # we have a wcs flag (from the image by default) so that users can skip over
    # plotting in physical coordinates. This is especially important
    # if the WCS headers are mallformed in some way.
    showgrid = get(plotattributes, :xgrid, true) && get(plotattributes, :ygrid, true)
    # Display a title giving our position along unplotted dimensions
    if length(refdims(imgv)) > 0
        if showwcstitle
            refdimslabel = join(map(refdims(imgv)) do d
                # match dimension with the wcs axis number
                i = wcsax(imgv, d)
                ct = wcs(imgv, wcsn).ctype[i]
                label = ctype_label(ct, wcs(imgv, wcsn).radesys)
                if label == "NONE"
                    label = name(d)
                end
                value = pix_to_world(imgv, [1,1]; wcsn, all=true, parent=true)[i]
                unit = wcs(imgv, wcsn).cunit[i]
                if ct == "STOKES"
                    return _stokes_name(_stokes_symbol(value))
                else
                    return @sprintf("%s = %.5g %s", label, value, unit)
                end
            end, ", ")
        else
            refdimslabel = join(map(d->"$(name(d))= $(d[1])", refdims(imgv)), ", ")
        end
        title --> refdimslabel
    end

    # To ensure the physical axis tick labels are correct the axes must be
    # tight to the image
    xl = (first(dims(imgv,1))-0.5)*platescale, (last(dims(imgv,1))+0.5)*platescale
    yl = (first(dims(imgv,2))-0.5)*platescale, (last(dims(imgv,2))+0.5)*platescale
    ylims --> yl
    xlims --> xl

    subplot_i = 0
    # Actual image series (RGB pixels by this point)
    @series begin
        subplot_i += 1
        subplot := subplot_i
        colorbar := false
        aspect_ratio --> 1

        # Note: if the axes are on unusual sides (e.g. y-axis at right, x-axis at top)
        # then these coordinates are not correct. They are only correct exactly
        # along the axis.
        # In astropy, the ticks are actually tilted to reflect this, though in general
        # the transformation from pixel to coordinates can be non-linear and curved.

        if showwcsticks
            xticks --> (gridspec.tickpos1x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,1)), gridspec.tickpos1w))
            xguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,1))], wcs(imgv, wcsn).radesys)

            yticks --> (gridspec.tickpos2x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,2)), gridspec.tickpos2w))
            yguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,2))], wcs(imgv, wcsn).radesys)
        end

        ax1 = collect(parent(dims(imgv,1))) .* platescale
        ax2 = collect(parent(dims(imgv,2))) .* platescale
        # Views of images are not currently supported by plotly() so we have to collect them.
        # ax1, ax2, view(parent(imgv), reverse(axes(imgv,1)),:)
        ax1, ax2, parent(imgv)[reverse(axes(imgv,1)),:]
    end

    # If wcs=true (default) and grid=true (not default), overplot a WCS
    # grid.
    if showgrid && showwcsticks

        # Plot the WCSGrid as a second series (actually just lines)
        @series begin
            subplot := 1
            # Use a default grid color that shows up across more
            # color maps
            if !haskey(plotattributes, :xforeground_color_grid) && !haskey(plotattributes, :yforeground_color_grid)
                gridcolor --> :lightgray
            end

            wcsg, gridspec
        end
    end

    # Disable the colorbar.
    # Plots.jl does not give us sufficient control to make sure the range and ticks
    # are correct after applying a non-linear stretch.
    # We attempt to make our own colorbar using a second plot.
    showcolorbar = !(T <: Colorant) && get(plotattributes, :colorbar, true) != :none
    if T <: Complex
        layout := @layout [
            imgmag{0.5h}
          imgangle{0.5h}
        ]
    end
    if showcolorbar
        if T <: Complex
            layout := @layout [
                imgmag{0.95w, 0.5h}       colorbar{0.5h}
              imgangle{0.95w, 0.5h}  colorbarangle{0.5h}
            ]
        else
            layout := @layout [
                img{0.95w} colorbar
            ]
        end
        colorbar_title = get(plotattributes, :colorbar_title, "")
        if !haskey(plotattributes, :colorbar_title)
            if haskey(header(img), "UNIT")
                colorbar_title = string(img[:UNIT])
            elseif haskey(header(img), "BUNIT")
                colorbar_title = string(img[:BUNIT])
            end
        end

        subplot_i += 1
        @series begin
            subplot := subplot_i
            aspect_ratio := :none
            colorbar := false
            cbimg, cbticks = imview_colorbar(img; clims, stretch, cmap, contrast, bias)
            xticks := []
            ymirror := true
            yticks := cbticks
            yguide := colorbar_title
            xguide := ""
            xlims := Tuple(axes(cbimg, 2))
            ylims := Tuple(axes(cbimg, 2))
            title := ""
            # Views of images are not currently supported by plotly so we have to collect them
            # view(cbimg, reverse(axes(cbimg,1)),:)
            cbimg[reverse(axes(cbimg,1)),:]
        end
    end

    # TODO: refactor to reduce duplication
    if T <: Complex
        img = angle.(data)
        img["UNIT"] = "angle (rad)"
        imgv = imview(img, clims=(-1pi, 1pi),stretch=identity, cmap=:cyclic_mygbm_30_95_c78_n256_s25)
        @series begin
            subplot_i += 1
            subplot := subplot_i
            colorbar := false
            title := ""
            aspect_ratio --> 1

            # Note: if the axes are on unusual sides (e.g. y-axis at right, x-axis at top)
            # then these coordinates are not correct. They are only correct exactly
            # along the axis.
            # In astropy, the ticks are actually tilted to reflect this, though in general
            # the transformation from pixel to coordinates can be non-linear and curved.

            if showwcsticks
                xticks --> (gridspec.tickpos1x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,1)), gridspec.tickpos1w))
                xguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,1))], wcs(imgv, wcsn).radesys)

                yticks --> (gridspec.tickpos2x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,2)), gridspec.tickpos2w))
                yguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,2))], wcs(imgv, wcsn).radesys)
            end

            ax1 = collect(parent(dims(imgv,1))) .* platescale
            ax2 = collect(parent(dims(imgv,2))) .* platescale
            # Views of images are not currently supported by plotly() so we have to collect them.
            # ax1, ax2, view(parent(imgv), reverse(axes(imgv,1)),:)
            ax1, ax2, parent(imgv)[reverse(axes(imgv,1)),:]
        end

        if showcolorbar
            colorbar_title = get(plotattributes, :colorbar_title, "")
            if !haskey(plotattributes, :colorbar_title) && haskey(header(img), "UNIT")
                colorbar_title = string(img[:UNIT])
            end

            @series begin
                subplot_i += 1
                subplot := subplot_i
                aspect_ratio := :none
                colorbar := false
                cbimg, _ = imview_colorbar(img; stretch=identity, clims=(-pi, pi), cmap=:cyclic_mygbm_30_95_c78_n256_s25)
                xticks := []
                ymirror := true
                ax = axes(cbimg,1)
                yticks := ([first(ax), mean(ax), last(ax)], ["-π", "0", "π"])
                yguide := colorbar_title
                xguide := ""
                xlims := Tuple(axes(cbimg, 2))
                ylims := Tuple(axes(cbimg, 2))
                title := ""
                view(cbimg, reverse(axes(cbimg,1)),:)
            end
        end

    end

    return
end

"""
    implot(
        img::AbstractArray;
        clims=Percent(99.5),
        stretch=identity,
        cmap=:magma,
        bias=0.5,
        contrast=1,
        wcsticks=true,
        grid=true,
        platescale=1
    )

Create a read only view of an array or AstroImageMat mapping its data values
to an array of Colors. Equivalent to:

    implot(
        imview(
            img::AbstractArray;
            clims=Percent(99.5),
            stretch=identity,
            cmap=:magma,
            bias=0.5,
            contrast=1,
        ),
        wcsn=1,
        wcsticks=true,
        wcstitle=true,
        grid=true,
        platescale=1
    )

### Image Rendering
See `imview` for how data is mapped to RGBA pixel values.

### WCS & Image Coordinates
If provided with an AstroImage that has WCS headers set, the tick marks and plot grid
are calculated using WCS.jl. By default, use the first WCS coordinate system.
The underlying pixel coordinates are those returned by `dims(img)` multiplied by `platescale`.
This allows you to overplot lines, regions, etc. using pixel coordinates.
If you wish to compute the pixel coordinate of a point in world coordinates, see `world_to_pix`.

* `wcsn` (default `1`) select which WCS transform in the headers to use for ticks & grid
* `wcsticks` (default `true` if WCS headers present) display ticks and labels, and title using world coordinates
* `wcstitle` (default `true` if WCS headers present and `length(refdims(img))>0`). When slicing a cube, display the location along unseen axes in world coordinates instead of pixel coordinates.
* `grid` (default `true`) show a grid over the plot. Uses WCS coordinates if `wcsticks` is true, otherwise pixel coordinates multiplied by `platescale`.
* `platescale` (default `1`). Scales the underlying pixel coordinates to ease overplotting, etc. If `wcsticks` is false, the displayed pixel coordinates are also scaled.

### Defaults
The default values of `clims`, `stretch`, and `cmap` are `extrema`, `identity`, and `nothing`
respectively.
You may alter these defaults using `AstroImages.set_clims!`,  `AstroImages.set_stretch!`, and
`AstroImages.set_cmap!`.
"""
implot

struct WCSGrid
    img::AstroImage
    extent::NTuple{4,Float64}
    wcsn::Int
end

"""
    wcsticks(img, axnum)

Generate nice tick labels for an AstroImageMat along axis `axnum`
Returns a vector of pixel positions and a vector of strings.

Example:
plot(img, xticks=wcsticks(WCSGrid(img), 1), yticks=wcsticks(WCSGrid(img), 2))
"""
function wcsticks(wcsg::WCSGrid, axnum, gs = wcsgridspec(wcsg))
    tickposx = axnum == 1 ? gs.tickpos1x : gs.tickpos2x
    tickposw = axnum == 1 ? gs.tickpos1w : gs.tickpos2w
    return tickposx, wcslabels(
        wcs(wcsg.img, wcsg.wcsn),
        axnum,
        tickposw
    )
end

# Function to generate nice string coordinate labels given a WCSTransform, axis number,
# and a vector of tick positions in world coordinates.
# This is used for labelling ticks and for annotating grid lines.
function wcslabels(w::WCSTransform, axnum, tickposw)

    if length(tickposw) == 0
        return String[]
    end

    # Select a unit converter (e.g. 12.12 -> (a,b,c,d)) and list of units
    if w.cunit[axnum] == "deg"
        if startswith(uppercase(w.ctype[axnum]), "RA")
            converter = deg2hms
            units = hms_units
        else
            converter = deg2dmsmμ
            units = dmsmμ_units
        end
    else
        converter = x->(x,)
        units = ("",)
    end

    # Format inital ticklabel
    ticklabels = fill("", length(tickposw))
    # We only include the part of the label that has changed since the last time.
    # Split up coordinates into e.g. sexagesimal
    parts = map(tickposw) do w
        vals = converter(w)
        return vals
    end

    # Start with something impossible of the same size:
    last_coord = Inf .* converter(first(tickposw))
    zero_coords_i = maximum(map(parts) do vals
        changing_coord_i = findfirst(vals .!= last_coord)
        if isnothing(changing_coord_i)
            changing_coord_i = 1
        end
        last_coord = vals
        return changing_coord_i
    end)

    # Loop through using only the relevant part of the label
    # Start with something impossible of the same size:
    last_coord = Inf .* converter(first(tickposw))
    for (i,vals) in enumerate(parts)
        changing_coord_i = findfirst(vals .!= last_coord)
        if isnothing(changing_coord_i)
            changing_coord_i = 1
        end
        # Don't display just e.g. 00" when we should display 50'00"
        if changing_coord_i > 1 && vals[changing_coord_i] == 0
            changing_coord_i = changing_coord_i -1
        end
        val_unit_zip = zip(vals[changing_coord_i:zero_coords_i],units[changing_coord_i:zero_coords_i])
        ticklabels[i] = mapreduce(*, enumerate(val_unit_zip)) do (coord_i,(val,unit))
            # If the last coordinate we print if the last coordinate we have available,
            # display it with decimal places
            if coord_i + changing_coord_i - 1== length(vals)
                str = @sprintf("%.2f", val)
            else
                str = @sprintf("%02d", val)
            end
            if length(str) > 0
                return str * unit
            else
                return str
            end
        end
        last_coord = vals
    end

    return ticklabels
end

# Extended form of deg2dms that further returns mas, microas.
function deg2dmsmμ(deg)
    d,m,s = deg2dms(deg)
    s_f = floor(s)
    mas = (s - s_f)*1e3
    mas_f = floor(mas)
    μas = (mas - mas_f)*1e3
    return (d,m,s_f,mas_f,μas)
end
const dmsmμ_units = [
    "°",
    "'",
    "\"",
    "mas",
    "μas",
]
const hms_units = [
    "ʰ",
    "ᵐ",
    "ˢ",
]

function ctype_label(ctype,radesys)
    if length(ctype) == 0
        return radesys
    elseif startswith(ctype, "RA")
        return "Right Ascension ($(radesys))"
    elseif startswith(ctype, "GLON")
        return "Galactic Longitude"
    elseif startswith(ctype, "TLON")
        return "ITRS"
    elseif startswith(ctype, "DEC")
        return "Declination ($(radesys))"
    elseif startswith(ctype, "GLAT")
        return "Galactic Latitude"
    # elseif startswith(ctype, "TLAT")
    elseif ctype == "STOKES"
        return "Polarization"
    else
        return ctype
    end
end

"""
    WCSGrid(img::AstroImageMat, ax=(1,2), coords=(first(axes(img,ax[1])),first(axes(img,ax[2]))))

Given an AstroImageMat, return information necessary to plot WCS gridlines in physical
coordinates against the image's pixel coordinates.
This function has to work on both plotted axes at once to handle rotation and general
curvature of the WCS grid projected on the image coordinates.

"""
function WCSGrid(img::AstroImageMat, wcsn=1)
    minx = first(dims(img,2))
    maxx = last(dims(img,2))
    miny = first(dims(img,1))
    maxy = last(dims(img,1))
    extent = (minx-0.5, maxx+0.5, miny-0.5, maxy+0.5)
    @show extent
    return WCSGrid(img, extent, wcsn)
end

# Recipe for a WCSGrid with lines, optional ticks (on by default),
# and optional grid labels (off by defaut).
# The AstroImageMat plotrecipe uses this recipe for grid lines if `grid=true`.
@recipe function f(wcsg::WCSGrid, gridspec=wcsgridspec(wcsg))
    label --> ""
    xs, ys = wcsgridlines(gridspec)

    if haskey(plotattributes, :foreground_color_grid)
        color --> plotattributes[:foreground_color_grid]
    elseif haskey(plotattributes, :foreground_color)
        color --> plotattributes[:foreground_color]
    else
        color --> :black
    end
    if haskey(plotattributes, :foreground_color_text)
        textcolor = plotattributes[:foreground_color_text]
    else
        textcolor = plotattributes[:color]
    end
    annotate = haskey(plotattributes, :gridlabels) && plotattributes[:gridlabels]

    xguide --> ctype_label(wcs(wcsg.img, wcsg.wcsn).ctype[wcsax(wcsg.img, dims(wcsg.img,1))], wcs(wcsg.img, wcsg.wcsn).radesys)
    yguide --> ctype_label(wcs(wcsg.img, wcsg.wcsn).ctype[wcsax(wcsg.img, dims(wcsg.img,2))], wcs(wcsg.img, wcsg.wcsn).radesys)

    xlims --> wcsg.extent[1], wcsg.extent[2]
    ylims --> wcsg.extent[3], wcsg.extent[4]

    grid := false
    tickdirection := :none

    xticks --> wcsticks(wcsg, 1, gridspec)
    yticks --> wcsticks(wcsg, 2, gridspec)

    @series xs, ys

    # We can optionally annotate the grid with their coordinates.
    # These come after the grid lines so they appear overtop.
    if annotate
        @series begin
            # TODO: why is this reverse necessary?
            rotations = reverse(rad2deg.(gridspec.annotations1θ))
            ticklabels = wcslabels(wcs(wcsg.img), 1, gridspec.annotations1w)
            seriestype := :line
            linewidth := 0
            # TODO: we need to use requires to load in Plots for the necessary text control. Future versions of RecipesBase might fix this.
            series_annotations := [
                Main.Plots.text(" $l", :right, :bottom, textcolor, 8, rotation=(-95 <= r <= 95) ? r : r+180)
                for (l, r) in zip(ticklabels, rotations)
            ]
            gridspec.annotations1x, gridspec.annotations1y
        end
        @series begin
            rotations  = rad2deg.(gridspec.annotations2θ)
            ticklabels = wcslabels(wcs(wcsg.img), 2, gridspec.annotations2w)
            seriestype := :line
            linewidth := 0
            series_annotations := [
                Main.Plots.text(" $l", :right, :bottom, textcolor, 8, rotation=(-95 <= r <= 95) ? r : r+180)
                for (l, r) in zip(ticklabels, rotations)
            ]
            gridspec.annotations2x, gridspec.annotations2y
        end

    end

    return
end

# Helper: true if all elements in vector are equal to each other.
allequal(itr) = all(==(first(itr)), itr)

# This function is responsible for actually laying out grid lines for a WCSGrid,
# ensuring they don't exceed the plot bounds, finding where they intersect the axes,
# and picking tick locations at the appropriate intersections with the left and
# bottom axes.
function wcsgridspec(wsg::WCSGrid)
    # Most of the complexity of this function is making sure everything
    # generalizes to N different, possiby skewed axes, where a change in
    # the opposite coordinate or even an unplotted coordinate affects
    # the grid.

    # x and y denote pixel coordinates (along `ax`), u and v are world coordinates roughly along same.
    minx, maxx, miny, maxy = wsg.extent

    # Find the extent of this slice in world coordinates
    posxy = [
        minx minx maxx maxx
        miny maxy miny maxy
    ]
    posuv = pix_to_world(wsg.img, posxy; wsg.wcsn, parent=true)
    (minu, maxu), (minv, maxv) = extrema(posuv, dims=2)

    # In general, grid can be curved when plotted back against the image,
    # so we will need to sample multiple points along the grid.
    # TODO: find a good heuristic for this based on the curvature.
    N_points = 50
    urange = range(minu, maxu, length=N_points)
    vrange = range(minv, maxv, length=N_points)

    # Find nice grid spacings using PlotUtils.optimize_ticks
    # These heuristics can probably be improved
    # TODO: this does not handle coordinates that wrap around
    Q=[(1.0,1.0), (3.0, 0.8), (2.0, 0.7), (5.0, 0.5)]
    k_min = 3
    k_ideal = 5
    k_max = 10

    tickpos2x = Float64[]
    tickpos2w = Float64[]
    gridlinesxy2 = NTuple{2,Vector{Float64}}[]
    # Not all grid lines will intersect the x & y axes nicely.
    # If we don't get enough valid tick marks (at least 2) loop again
    # requesting more locations up to three times.
    local tickposv
    j = 5
    while length(tickpos2x) < 2 && j > 0
        k_min += 2
        k_ideal += 2
        k_max += 2
        j -= 1

        tickposv = optimize_ticks(6minv, 6maxv; Q, k_min, k_ideal, k_max)[1]./6

        empty!(tickpos2x)
        empty!(tickpos2w)
        empty!(gridlinesxy2)
        for tickv in tickposv
            # Make sure we handle unplotted slices correctly.
            griduv = repeat(posuv[:,1], 1, N_points)
            griduv[1,:] .= urange
            griduv[2,:] .= tickv
            posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)

            # Now that we have the grid in pixel coordinates,
            # if we find out where the grid intersects the axes we can put
            # the labels in the correct spot

            # We can use these masks to determine where, and in what direction
            # the gridlines leave the plot extent
            in_horz_ax = minx .<=  posxy[1,:] .<= maxx
            in_vert_ax = miny .<=  posxy[2,:] .<= maxy
            in_axes = in_horz_ax .& in_vert_ax
            if count(in_axes) < 2
                continue
            elseif all(in_axes)
                point_entered = [
                    posxy[1,begin]
                    posxy[2,begin]
                ]
                point_exitted = [
                    posxy[1,end]
                    posxy[2,end]
                ]
            elseif allequal(posxy[1,findfirst(in_axes):findlast(in_axes)])
                point_entered = [
                    posxy[1,max(begin,findfirst(in_axes)-1)]
                    # posxy[2,max(begin,findfirst(in_axes)-1)]
                    miny
                ]
                point_exitted = [
                    posxy[1,min(end,findlast(in_axes)+1)]
                    # posxy[2,min(end,findlast(in_axes)+1)]
                    maxy
                ]
            # Vertical grid lines
            elseif allequal(posxy[2,findfirst(in_axes):findlast(in_axes)])
                point_entered = [
                    minx #posxy[1,max(begin,findfirst(in_axes)-1)]
                    posxy[2,max(begin,findfirst(in_axes)-1)]
                ]
                point_exitted = [
                    maxx #posxy[1,min(end,findlast(in_axes)+1)]
                    posxy[2,min(end,findlast(in_axes)+1)]
                ]
            else
                # Use the masks to pick an x,y point inside the axes and an
                # x,y point outside the axes.
                i = findfirst(in_axes)
                x1 = posxy[1,i]
                y1 = posxy[2,i]
                x2 = posxy[1,i+1]
                y2 = posxy[2,i+1]
                if x2-x1 ≈ 0
                    @warn "undef slope"
                end

                # Fit a line where we cross the axis
                m1 = (y2-y1)/(x2-x1)
                b1 = y1-m1*x1
                # If the line enters via the vertical axes...
                if findfirst(in_vert_ax) <= findfirst(in_horz_ax)
                    # Then we simply evaluate it at that axis
                    x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx
                    x = clamp(x,minx,maxx)
                    y = m1*x+b1
                else
                    # We must find where it enters the plot from
                    # bottom or top
                    x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b1)/m1 : (miny-b1)/m1
                    x = clamp(x,minx,maxx)
                    y = m1*x+b1
                end

                # From here, do a linear fit to find the intersection with the axis.
                point_entered = [
                    x
                    y
                ]

                # Use the masks to pick an x,y point inside the axes and an
                # x,y point outside the axes.
                i = findlast(in_axes)
                x1 = posxy[1,i-1]
                y1 = posxy[2,i-1]
                x2 = posxy[1,i]
                y2 = posxy[2,i]
                if x2-x1 ≈ 0
                    @warn "undef slope"
                end

                # Fit a line where we cross the axis
                m2 = (y2-y1)/(x2-x1)
                b2 = y2-m2*x2
                if findlast(in_vert_ax) > findlast(in_horz_ax)
                    # Then we simply evaluate it at that axis
                    x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx
                    x = clamp(x,minx,maxx)
                    y = m2*x+b2
                else
                    # We must find where it enters the plot from
                    # bottom or top
                    x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b2)/m2 : (miny-b2)/m2
                    x = clamp(x,minx,maxx)
                    y = m2*x+b2
                end

                # From here, do a linear fit to find the intersection with the axis.
                point_exitted = [
                    x
                    y
                ]
            end

            if point_entered[1] == minx
                push!(tickpos2x, point_entered[2])
                push!(tickpos2w, tickv)
            end
            if point_exitted[1] == minx
                push!(tickpos2x, point_exitted[2])
                push!(tickpos2w, tickv)
            end

            posxy_neat = [point_entered  posxy[[1,2],in_axes] point_exitted]
            # posxy_neat = posxy
            # TODO: do unplotted other axes also need a fit?

            gridlinexy = (
                posxy_neat[1,:],
                posxy_neat[2,:]
            )
            push!(gridlinesxy2, gridlinexy)
        end
    end

    # Then do the opposite coordinate
    k_min = 3
    k_ideal = 5
    k_max = 10
    tickpos1x = Float64[]
    tickpos1w = Float64[]
    gridlinesxy1 = NTuple{2,Vector{Float64}}[]
    # Not all grid lines will intersect the x & y axes nicely.
    # If we don't get enough valid tick marks (at least 2) loop again
    # requesting more locations up to three times.
    local tickposu
    j = 5
    while length(tickpos1x) < 2 && j > 0
        k_min += 2
        k_ideal += 2
        k_max += 2
        j -= 1

        tickposu = optimize_ticks(6minu, 6maxu; Q, k_min, k_ideal, k_max)[1]./6

        empty!(tickpos1x)
        empty!(tickpos1w)
        empty!(gridlinesxy1)
        for ticku in tickposu
            # Make sure we handle unplotted slices correctly.
            griduv = repeat(posuv[:,1], 1, N_points)
            griduv[1,:] .= ticku
            griduv[2,:] .= vrange
            posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)

            # Now that we have the grid in pixel coordinates,
            # if we find out where the grid intersects the axes we can put
            # the labels in the correct spot

            # We can use these masks to determine where, and in what direction
            # the gridlines leave the plot extent
            in_horz_ax = minx .<=  posxy[1,:] .<= maxx
            in_vert_ax = miny .<=  posxy[2,:] .<= maxy
            in_axes = in_horz_ax .& in_vert_ax

            if count(in_axes) < 2
                continue
            elseif all(in_axes)
                point_entered = [
                    posxy[1,begin]
                    posxy[2,begin]
                ]
                point_exitted = [
                    posxy[1,end]
                    posxy[2,end]
                ]
            # Horizontal grid lines
            elseif allequal(posxy[1,findfirst(in_axes):findlast(in_axes)])
                point_entered = [
                    posxy[1,findfirst(in_axes)]
                    miny
                ]
                point_exitted = [
                    posxy[1,findlast(in_axes)]
                    maxy
                ]
                # push!(tickpos1x, posxy[1,findfirst(in_axes)])
                # push!(tickpos1w, ticku)
            # Vertical grid lines
            elseif allequal(posxy[2,findfirst(in_axes):findlast(in_axes)])
                point_entered = [
                    minx
                    posxy[2,findfirst(in_axes)]
                ]
                point_exitted = [
                    maxx
                    posxy[2,findfirst(in_axes)]
                ]
            else
                # Use the masks to pick an x,y point inside the axes and an
                # x,y point outside the axes.
                i = findfirst(in_axes)
                x1 = posxy[1,i]
                y1 = posxy[2,i]
                x2 = posxy[1,i+1]
                y2 = posxy[2,i+1]
                if x2-x1 ≈ 0
                    @warn "undef slope"
                end

                # Fit a line where we cross the axis
                m1 = (y2-y1)/(x2-x1)
                b1 = y1-m1*x1
                # If the line enters via the vertical axes...
                if findfirst(in_vert_ax) < findfirst(in_horz_ax)
                    # Then we simply evaluate it at that axis
                    x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx
                    x = clamp(x,minx,maxx)
                    y = m1*x+b1
                else
                    # We must find where it enters the plot from
                    # bottom or top
                    x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b1)/m1 : (miny-b1)/m1
                    x = clamp(x,minx,maxx)
                    y = m1*x+b1
                end

                # From here, do a linear fit to find the intersection with the axis.
                point_entered = [
                    x
                    y
                ]

                # Use the masks to pick an x,y point inside the axes and an
                # x,y point outside the axes.
                i = findlast(in_axes)
                x1 = posxy[1,i-1]
                y1 = posxy[2,i-1]
                x2 = posxy[1,i]
                y2 = posxy[2,i]
                if x2-x1 ≈ 0
                    @warn "undef slope"
                end

                # Fit a line where we cross the axis
                m2 = (y2-y1)/(x2-x1)
                b2 = y2-m2*x2
                if findlast(in_vert_ax) > findlast(in_horz_ax)
                    # Then we simply evaluate it at that axis
                    x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx
                    x = clamp(x,minx,maxx)
                    y = m2*x+b2
                else
                    # We must find where it enters the plot from
                    # bottom or top
                    x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b2)/m2 : (miny-b2)/m2
                    x = clamp(x,minx,maxx)
                    y = m2*x+b2
                end

                # From here, do a linear fit to find the intersection with the axis.
                point_exitted = [
                    x
                    y
                ]
            end

            posxy_neat = [point_entered  posxy[[1,2],in_axes] point_exitted]
            # TODO: do unplotted other axes also need a fit?

            if point_entered[2] == miny
                push!(tickpos1x, point_entered[1])
                push!(tickpos1w, ticku)
            end
            if point_exitted[2] == miny
                push!(tickpos1x, point_exitted[1])
                push!(tickpos1w, ticku)
            end

            gridlinexy = (
                posxy_neat[1,:],
                posxy_neat[2,:]
            )
            push!(gridlinesxy1, gridlinexy)
        end
    end

    # Grid annotations are simpler:
    annotations1w = Float64[]
    annotations1x = Float64[]
    annotations1y = Float64[]
    annotations1θ = Float64[]
    for ticku in tickposu
        # Make sure we handle unplotted slices correctly.
        griduv = posuv[:,1]
        griduv[1] = ticku
        griduv[2] = mean(vrange)
        posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)
        if !(minx < posxy[1] < maxx) ||
            !(miny < posxy[2] < maxy)
            continue
        end
        push!(annotations1w, ticku)
        push!(annotations1x, posxy[1])
        push!(annotations1y, posxy[2])

        # Now find slope (TODO: stepsize)
        # griduv[ax[2]] -= 1
        griduv[2] += 0.1step(vrange)
        posxy2 = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)
        θ = atan(
            posxy2[2] - posxy[2],
            posxy2[1] - posxy[1],
        )
        push!(annotations1θ, θ)
    end
    annotations2w = Float64[]
    annotations2x = Float64[]
    annotations2y = Float64[]
    annotations2θ = Float64[]
    for tickv in tickposv
        # Make sure we handle unplotted slices correctly.
        griduv = posuv[:,1]
        griduv[1] = mean(urange)
        griduv[2] = tickv
        posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)
        if !(minx < posxy[1] < maxx) ||
            !(miny < posxy[2] < maxy)
            continue
        end
        push!(annotations2w, tickv)
        push!(annotations2x, posxy[1])
        push!(annotations2y, posxy[2])

        griduv[1] += 0.1step(urange)
        posxy2 = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true)
        θ = atan(
            posxy2[2] - posxy[2],
            posxy2[1] - posxy[1],
        )
        push!(annotations2θ, θ)
    end

    return (;
        gridlinesxy1,
        gridlinesxy2,
        tickpos1x,
        tickpos1w,
        tickpos2x,
        tickpos2w,

        annotations1w,
        annotations1x,
        annotations1y,
        annotations1θ,

        annotations2w,
        annotations2x,
        annotations2y,
        annotations2θ,
    )
end

# From a WCSGrid, return just the grid lines as a single pair of x & y coordinates
# suitable for plotting.
function wcsgridlines(wcsg::WCSGrid)
    return wcsgridlines(wcsgridspec(wcsg))
end
function wcsgridlines(gridspec::NamedTuple)
    # Unroll grid lines into a single series separated by NaNs
    xs1 = mapreduce(vcat, gridspec.gridlinesxy1, init=Float64[]) do gridline
        return vcat(gridline[1], NaN)
    end
    ys1 = mapreduce(vcat, gridspec.gridlinesxy1, init=Float64[]) do gridline
        return vcat(gridline[2], NaN)
    end
    xs2 = mapreduce(vcat, gridspec.gridlinesxy2, init=Float64[]) do gridline
        return vcat(gridline[1], NaN)
    end
    ys2 = mapreduce(vcat, gridspec.gridlinesxy2, init=Float64[]) do gridline
        return vcat(gridline[2], NaN)
    end

    xs = vcat(xs1, NaN, xs2)
    ys = vcat(ys1, NaN, ys2)
    return xs, ys
end

@userplot PolQuiver
@recipe function f(h::PolQuiver)
    cube = only(h.args)
    bins = get(plotattributes, :bins, 4)
    ticklen = get(plotattributes, :ticklen, nothing)
    minpol = get(plotattributes, :minpol, 0.1)

    i = cube[Pol=At(:I)]
    q = cube[Pol=At(:Q)]
    u = cube[Pol=At(:U)]
    polinten = @. sqrt(q^2 + u^2)
    linpolfrac = polinten ./ i

    binratio=1/bins
    xs = imresize([x for x in dims(cube,1), y in dims(cube,2)], ratio=binratio)
    ys = imresize([y for x in dims(cube,1), y in dims(cube,2)], ratio=binratio)
    qx = imresize(q, ratio=binratio)
    qy = imresize(u, ratio=binratio)
    qlinpolfrac = imresize(linpolfrac, ratio=binratio)
    qpolintenr = imresize(polinten, ratio=binratio)

    # We want the longest ticks to be around 1 bin long by default.
    qmaxlen = quantile(filter(isfinite,qpolintenr), 0.98)
    if isnothing(ticklen)
        a = bins / qmaxlen
    else
        a = ticklen / qmaxlen
    end
    # Only show arrows where the data is finite, and more than a couple pixels
    # long.
    mask = (isfinite.(qpolintenr)) .& (qpolintenr .>= minpol.*qmaxlen)
    pointstmp = map(xs[mask],ys[mask],qx[mask],qy[mask]) do x,y,qxi,qyi
        return ([x, x+a*qxi, NaN], [y, y+a*qyi, NaN])
    end
    xs = reduce(vcat, getindex.(pointstmp, 1))
    ys = reduce(vcat, getindex.(pointstmp, 2))

    colors = qlinpolfrac[mask]
    if !isnothing(colors)
        line_z := repeat(colors, inner=3)
    end

    label --> ""
    color --> :turbo
    framestyle --> :box
    aspect_ratio --> 1
    linewidth --> 1.5
    colorbar --> true
    colorbar_title --> "Linear polarization fraction"

    xl = first(dims(i,2)), last(dims(i,2))
    yl = first(dims(i,1)), last(dims(i,1))
    ylims --> yl
    xlims --> xl

    @series begin
        xs, ys
    end
end

"""
    polquiver(polqube::AstroImage)

Given a data cube (of at least 2 spatial dimensions, plus a polarization axis),
plot a vector field of polarization data.
The tick length represents the polarization intensity, sqrt(q^2 + u^2),
and the color represents the linear polarization fraction, sqrt(q^2 + u^2) / i.

There are several ways you can adjust the appearance of the plot using keyword arguments:
* `bins` (default = 1) By how much should we bin down the polarization data before drawing the ticks? This reduced clutter from higher resolution datasets. Can be fractional.
* `ticklen` (default = bins) How long the 98th percentile arrow should be. By default, 1 bin long. Make this larger to draw longer arrows.
* `color` (default = :turbo) What colorscheme should be used for linear polarization fraction.
* `minpol` (default = 0.2) Hides arrows that are shorter than `minpol` times the 98th percentile arrow to make a cleaner image. Set to 0 to display all data.

Use `implot` and `polquiver!` to overplot polarization data over an image.
"""
polquiver


text:
# COV_EXCL_START
using .CUDA, Cassette

cuda_is_loaded = true

#! format: off
const cudafuns = (
    :cos, :cospi, :sin, :sinpi, :tan,
    :acos, :asin, :atan,
    :cosh, :sinh, :tanh,
    :acosh, :asinh, :atanh,
    :log, :log10, :log1p, :log2,
    :exp, :exp2, :exp10, :expm1, :ldexp,
    :abs,
    :sqrt, :cbrt,
    :ceil, :floor,
)
#! format: on

Cassette.@context CeedCudaContext

@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.kwfunc), f)
    return Core.kwfunc(f)
end
@inline function Cassette.overdub(::CeedCudaContext, ::typeof(Core.apply_type), args...)
    return Core.apply_type(args...)
end
@inline function Cassette.overdub(
    ::CeedCudaContext,
    ::typeof(StaticArrays.Size),
    x::Type{<:AbstractArray{<:Any,N}},
) where {N}
    return StaticArrays.Size(x)
end

for f in cudafuns
    @eval @inline function Cassette.overdub(
        ::CeedCudaContext,
        ::typeof(Base.$f),
        x::Union{Float32,Float64},
    )
        return CUDA.$f(x)
    end
end

function setarray!(v::CeedVector, mtype::MemType, cmode::CopyMode, arr::CuArray)
    ptr = Ptr{CeedScalar}(UInt64(pointer(arr)))
    C.CeedVectorSetArray(v[], mtype, cmode, ptr)
    if cmode == USE_POINTER
        v.arr = arr
    end
end

struct FieldsCuda
    inputs::NTuple{16,Int}
    outputs::NTuple{16,Int}
end

function generate_kernel(qf_name, kf, dims_in, dims_out)
    ninputs = length(dims_in)
    noutputs = length(dims_out)

    input_sz = prod.(dims_in)
    output_sz = prod.(dims_out)

    f_ins = [Symbol("rqi$i") for i = 1:ninputs]
    f_outs = [Symbol("rqo$i") for i = 1:noutputs]

    args = Vector{Union{Symbol,Expr}}(undef, ninputs + noutputs)
    def_ins = Vector{Expr}(undef, ninputs)
    f_ins_j = Vector{Union{Symbol,Expr}}(undef, ninputs)
    for i = 1:ninputs
        if length(dims_in[i]) == 0
            def_ins[i] = :(local $(f_ins[i]))
            f_ins_j[i] = f_ins[i]
            args[i] = f_ins[i]
        else
            def_ins[i] =
                :($(f_ins[i]) = LibCEED.MArray{Tuple{$(dims_in[i]...)},Float64}(undef))
            f_ins_j[i] = :($(f_ins[i])[j])
            args[i] = :(LibCEED.SArray{Tuple{$(dims_in[i]...)},Float64}($(f_ins[i])))
        end
    end
    for i = 1:noutputs
        args[ninputs+i] = f_outs[i]
    end

    def_outs = [
        :($(f_outs[i]) = LibCEED.MArray{Tuple{$(dims_out[i]...)},Float64}(undef))
        for i = 1:noutputs
    ]

    device_ptr_type = Core.LLVMPtr{CeedScalar,LibCEED.AS.Global}

    read_quads_in = [
        :(
            for j = 1:$(input_sz[i])
                $(f_ins_j[i]) = unsafe_load(
                    reinterpret($device_ptr_type, fields.inputs[$i]),
                    q + (j - 1)*Q,
                    a,
                )
            end
        ) for i = 1:ninputs
    ]

    write_quads_out = [
        :(
            for j = 1:$(output_sz[i])
                unsafe_store!(
                    reinterpret($device_ptr_type, fields.outputs[$i]),
                    $(f_outs[i])[j],
                    q + (j - 1)*Q,
                    a,
                )
            end
        ) for i = 1:noutputs
    ]

    qf = gensym(qf_name)
    quote
        function $qf(ctx_ptr, Q, fields)
            gd = LibCEED.gridDim()
            bi = LibCEED.blockIdx()
            bd = LibCEED.blockDim()
            ti = LibCEED.threadIdx()

            inc = bd.x*gd.x

            $(def_ins...)
            $(def_outs...)

            # Alignment for data read/write
            a = Val($(Base.datatype_alignment(CeedScalar)))

            # Cassette context for replacing intrinsics with CUDA versions
            ctx = LibCEED.CeedCudaContext()

            for q = (ti.x+(bi.x-1)*bd.x):inc:Q
                $(read_quads_in...)
                LibCEED.Cassette.overdub(ctx, $kf, ctx_ptr, $(args...))
                $(write_quads_out...)
            end
            return
        end
    end
end

function mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out)
    k_fn = Core.eval(def_module, generate_kernel(qf_name, kf, dims_in, dims_out))
    tt = Tuple{Ptr{Nothing},Int32,FieldsCuda}
    host_k = cufunction(k_fn, tt; maxregs=64)
    return host_k.fun.handle
end
# COV_EXCL_STOP

ext: jl

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
function Glonass(arg0::jdouble, arg1::AbsoluteDate, arg2::AbsoluteDate, arg3::ExtendedPVCoordinatesProvider, arg4::Frame)
    return Glonass((jdouble, AbsoluteDate, AbsoluteDate, ExtendedPVCoordinatesProvider, Frame), arg0, arg1, arg2, arg3, arg4)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Raylib"
version = v"4.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/raysan5/raylib.git",
              "0851960397f02a477d80eda2239f90fae14dec64"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/raylib/src/

atomic_patch -p1 ../../patches/add-missing-header.patch
atomic_patch -p1 ../../patches/make-install-everywhere.patch
atomic_patch -p2 ../../patches/make-ldflags-windows.patch

export CFLAGS="-D_POSIX_C_SOURCE=200112L"
if [[ "${target}" == *-freebsd* ]]; then
    # Allow definition of `u_char`, `u_short`, `u_int`, and `u_long` in sys/types.h
    CFLAGS="${CFLAGS} -D__BSD_VISIBLE"
fi
CFLAGS="${CFLAGS} -DSUPPORT_EVENTS_AUTOMATION -DSUPPORT_FILEFORMAT_BMP -DSUPPORT_FILEFORMAT_JPG"

FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    # raylib.rc.data is broken for x64, remove it from compilation
    sed -i 's+$(RAYLIB_RES_FILE)+ +g' Makefile
    # we need to specify the OS in the flags to make for Windows
    FLAGS+=(OS=Windows_NT)
fi
make -j${nproc} USE_EXTERNAL_GLFW=TRUE PLATFORM=PLATFORM_DESKTOP RAYLIB_LIBTYPE=SHARED RAYLIB_MODULE_RAYGUI=TRUE RAYLIB_MODULE_PHYSAC=TRUE "${FLAGS[@]}"
make install RAYLIB_LIBTYPE=SHARED DESTDIR="${prefix}" RAYLIB_INSTALL_PATH="${libdir}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=p->arch(p)=="armv6l")

# The products that we will ensure are always built
products = [
    LibraryProduct(["libraylib","raylib"], :libraylib)
]

x11_platforms = filter(p ->Sys.islinux(p) || Sys.isfreebsd(p), platforms)
# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="alsa_jll", uuid="45378030-f8ea-5b20-a7c7-1a9d95efb90e"); platforms=filter(Sys.islinux, platforms))
    Dependency(PackageSpec(name="Mesa_jll", uuid="78dcde23-ec64-5e07-a917-6fe22bbc0f45"); platforms=filter(Sys.iswindows, platforms))
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"); platforms=x11_platforms)
    Dependency(PackageSpec(name="Xorg_libXrandr_jll", uuid="ec84b674-ba8e-5d96-8ba1-2a689ba10484"); platforms=x11_platforms)
    Dependency(PackageSpec(name="Xorg_libXi_jll", uuid="a51aa0fd-4e3c-5386-b890-e753decda492"); platforms=x11_platforms)
    Dependency(PackageSpec(name="Xorg_libXcursor_jll", uuid="935fb764-8cf2-53bf-bb30-45bb1f8bf724"); platforms=x11_platforms)
    Dependency(PackageSpec(name="Xorg_libXinerama_jll", uuid="d1454406-59df-5ea1-beac-c340f2130bc3"); platforms=x11_platforms)
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
    Dependency(PackageSpec(name="GLFW_jll", uuid="0656b61e-2033-5cc2-a64a-77c0f6c09b89"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")


text:
export GenNeighbourhood
export GenNeighbourhoodIesimo

using Distributions

"""
    GenNeighbourhood(w::Vector{Real} ,µ::Real, var ::Real)
    Devuelve vecino del vector `w` de acorde a una normal
    de media `µ` y varianza `var`.
"""
function GenNeighbourhood(w::Vector{<:Real} ,µ::Real, var ::Real)::Vector{<:Real}
    distribution = Normal(µ, var)
    z=rand(distribution,size(w))
    w_new = map(
        x -> min(1, max(0,x)),
        w-z
    )
    return w_new
end

"""
GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, index::Int)::Vector{<:Real}
    Devuelve vecino del vector `w` de acorde a una normal
    de media `µ` y varianza `var` donde solo se ha cambia la componente iéxima.
"""
function GenNeighbourhoodIesimo(w::Vector{<:Real} ,µ::Real, var ::Real, index::Int)::Vector{<:Real}
    distribution = Normal(µ, var)
    z=rand(distribution)
    w_new = w
    w_new[index] = min(1, max(0,w_new[index]-z))
    return w_new
end

max_stars_repo_name: BlancaCC/Problema-aprendizaje-pesos

max_issues_repo_name: BlancaCC/Problema-aprendizaje-pesos

max_forks_repo_name: BlancaCC/Problema-aprendizaje-pesos


text:
module Milann

using Flux
import Flux.Tracker: data, @grad, track

using Statistics

# This implemens a MIL version where the bag instances are stored in a continuous tensor and bags are delimited by ranges.

export RangeMIL, segmax, segmean, segmax_naive, segmean_naive, segmaxmean

struct RangeMIL
    premodel
    aggregation
    postmodel
end
@Flux.treelike RangeMIL

function (m::RangeMIL)(X::AbstractArray, B::AbstractArray)
    (nfeatures, ninsts), nbags = size(X), length(B)
    # X: nfeatures x ninsts, B: nbags
    Ypre = m.premodel(X)
    # Ypre: npremodeloutput x ninstances
    Yagg = m.aggregation(Ypre, B)
    # Yagg: npremodeloutput x nbags
    m.postmodel(Yagg)
end

"""
segmax(X, segments)

Compute segmented maximum for instances `X` and bags indexed by `segments`.
"""
function segmax(X, segments)
    Y = similar(X, size(X, 1), length(segments))

    # naive approach, fast on CPU
    #     for (i, seg) in enumerate(segments)
    #         Y[:, i:i] .= mean(view(X, :, seg), dims=2)
    #     end

    # sliced approach much better for GPU
    Y .= 0
    mlen = maximum(length.(segments))
    for i in 1:mlen
        Yindices = [j for (j, e) in enumerate(segments) if length(e) >= i]
        Xindices = map(q -> first(q) + i - 1, filter(e -> length(e) >= i, segments))
        Y[:, Yindices] = max.(view(Y, :, Yindices), view(X, :, Xindices))
    end
    Y
end

"""
segmax_idx(X, segments)

Compute segmented maximum for instances `X` and bags indexed by `segments`. Same as segmax, but returns also the indices of maxima in input (used for grad computation).
"""
function segmax_idx(X, segments)
    Y = similar(X, size(X, 1), length(segments))
    Y2 = similar(X, size(Y)...) # twice the memmory!
    idxs = similar(X, size(Y)...)
    Y .= -Inf
    Y2 .= NaN
    idxs .= 0
    mlen = maximum(length.(segments))
    for i in 1:mlen
        Yindices = [j for (j, e) in enumerate(segments) if length(e) >= i]
        Xindices = map(q -> first(q) + i - 1, filter(e -> length(e) >= i, segments))
        Y2[:, Yindices] = max.(view(Y, :, Yindices), view(X, :, Xindices))
        idxs[:, Yindices] += view(Y, :, Yindices)  .!= view(Y2, :, Yindices) # has been max increased since last iteration?
        tmp = Y; Y = Y2; Y2 = tmp # just copy references
    end
    # idxs now holds indices of maxima relative to bags, make them relative to inputs
    offsets = first.(segments) .- 1 # first indices in X minus 1
    Y, Int.(cpu(idxs) .+ reshape(offsets, 1, length(offsets)))
end

segmax(X::TrackedArray, segments) = track(segmax, X, segments)

function dsegmax(Δ, segments, idxs)
    grads = similar(Δ, size(Δ, 1), last(last(segments)))
    grads .= 0
    for (i, row) in enumerate(eachrow(idxs))
        # @show size(row), minimum(row), size(Δ)
        grads[i, row] .= Δ[i, :]
    end
    grads, nothing
end

@grad function segmax(X, segments)
    Y, idxs = segmax_idx(data(X), segments)
    Y, Δ -> dsegmax(Δ, segments, idxs)
end

"""
segmean(X, segments)

Compute segmented mean for instances `X` and bags indexed by `segments`.
"""
function segmean(X, segments)
    Y = similar(X, size(X, 1), length(segments))

    # naive approach, fast on CPU
    #     for (i, seg) in enumerate(segments)
    #         Y[:, i:i] .= mean(view(X, :, seg), dims=2)
    #     end

    # sliced approach much better for GPU
    Y .= 0
    mlen = maximum(length.(segments))
    for i in 1:mlen
        Yindices = [j for (j, e) in enumerate(segments) if length(e) >= i]
        Xindices = map(q -> first(q) + i - 1, filter(e -> length(e) >= i, segments))
        Y[:, Yindices] += view(X, :, Xindices)
    end
    Y ./= reshape(gpu(length.(segments)), 1, length(segments))
    Y
end

segmean(X::TrackedArray, segments) = track(segmean, X, segments)

function dsegmean(Δ, segments)
    grads = similar(Δ, size(Δ, 1), last(last(segments)))

    # naive approach, fast on CPU
    #     for (i, seg) in enumerate(segments)
    #         grads[:, seg] .= view(Δ, :,i) / length(seg)
    #     end

    # sliced approach much better for GPU
    mlen = maximum(length.(segments))
    for i in 1:mlen
        Yindices = [j for (j, e) in enumerate(segments) if length(e) >= i]
        Xindices = map(q -> first(q) + i - 1, filter(e -> length(e) >= i, segments))
        grads[:, Xindices] .= view(Δ, :, Yindices)
    end
    for i in 2:mlen # normalization
        Xindices = vcat((collect(seg) for seg in segments if length(seg) == i)...)
        if length(Xindices) > 0
            grads[:, Xindices] ./= i
        end
    end

    grads, nothing
end

@grad segmean(X, segments) = segmean(data(X), segments), Δ -> dsegmean(Δ, segments)

"""
segmax_naive(X, segments)

Segmented maximum: naive reference implementation using default AD.
"""
function segmax_naive(X, segments)
    hcat((maximum(X[:,s], dims=ndims(X)) for s in segments)...)
end

"""
    segmean_naive(X, segments)

Segmented mean: naive reference implementation using default AD.
"""
function segmean_naive(X, segments)
    hcat((mean(X[:,s], dims=ndims(X)) for s in segments)...)
end

"""
    segmaxmean(X, segments)

Combine `segmax` and `segmean`.
"""
segmaxmean(X, segments) = vcat(segmax(X, segments), segmean(X, segments))

end # module


text:
module ITKFFI

include("image.jl")
include("filter.jl")

@wraptype Float2D
@wraptype Float3D
@wraptype Double2D
@wraptype Double3D
@wraptype Short2D
@wraptype Short3D

end

max_stars_repo_name: orobix/itkffi

max_issues_repo_name: orobix/itkffi

max_forks_repo_name: orobix/itkffi


text:
module SWScale
  include(joinpath(Pkg.dir("VideoIO"),"src","init.jl"))
  w(f) = joinpath(swscale_dir, f)

  using AVUtil
  include(w("LIBSWSCALE.jl"))
end


text:
for i in 1:2:9
    println(i)
end #> 1 3 5 7 9

typeof(1:1000) #> UnitRange{Int64}

a = split("A,B,C,D",",")
typeof(a) #> Array{SubString{String},1}
show(a) #> SubString{String}["A","B","C","D"]

arr = [100, 25, 37]
arra = Any[100, 25, "ABC"] #> element type Any
arr[1] #> 100
arr[end] #> 37
# arr[6]
#> ERROR: BoundsError: attempt to access 3-element Array{Int64,1} at index [6]
println(eltype(arr)) #> Int64
println(length(arr)) #> 3
println(ndims(arr)) #> 1
println(size(arr,1)) #> 3
println(size(arr)) #> (3,)

arr2 = Array{Int64}(undef, 5)
show(arr2) #> [1, 2, 3, 4, 5]
sizehint!(arr2, 10^5)

arr3 = Float64[] #> 0-element Array{Float64,1}
push!(arr3, 1.0) #> 1-element Array{Float64,1}

arr4 = collect(1:7) #> 7-element Array{Int64,1}
show(arr4) #> [1, 2, 3, 4, 5, 6, 7]

# for in loop and changing the array:
da = [1,2,3,4,5]
for n in da
    n *= 2
end
da #> 5-element Array{Int64,1}: 1 2 3 4 5
for i in 1:length(da)
    da[i] *= 2
end
da #> 5-element Array{Int64,1}: 2 4 6 8 10

arr4 = [1,2,3,4,5,6,7] #> 7-element Array{Int64,1}: [1,2,3,4,5,6,7]
join(arr4, ", ") #> "1, 2, 3, 4, 5, 6, 7"
arr4[1:3] # => [1, 2, 3]
arr4[4:end] # => [4, 5, 6, 7]

arr = [1,2,3,4,5]
arr[2:4] = [8,9,10]
println(arr) #> 1 8 9 10 5

zeros(5)
ones(4)
ones(3, 2)
eqa = range(0, step=10, length=5) #> 0:10:40
println()
show(eqa) #> 0:10:40
fill!(arr, 42) #> [42, 42, 42, 42, 42]
println()
println(Array{Any}(undef, 4)) #> Any[#undef,#undef,#undef,#undef]

v1 = rand(Int32, 5)
println(v1) #> Int32[1735276173,972339632,1303377282,1493859467,-788555652]

b = collect(1:7)
c = [100,200,300]
append!(b,c) # Now b is [1, 2, 3, 4, 5, 6, 7, 100, 200, 300]
pop!(b) #> 300, b is now [1, 2, 3, 4, 5, 6, 7, 100, 200]
push!(b, 42) # b is now [1, 2, 3, 4, 5, 6, 7, 100, 200, 42]
popfirst!(b) #> 1, b is now [2, 3, 4, 5, 6, 7, 100, 200, 42]
pushfirst!(b, 42) # b is now [1, 2, 3, 4, 5, 6, 7, 100, 200, 42]
splice!(b,8) #> 100, b is now [42, 2, 3, 4, 5, 6, 7, 200, 42]
in(42, b) #> true
in(43, b) #> false
println()
println("sorting:")
sort(b) #> [2,3,4,5,6,7,42,42,200], b is not changed
println(b) #> [42,2,3,4,5,6,7,200,42]
sort!(b) #> b is now changed to [2,3,4,5,6,7,42,42,200]
println(b) #> [2,3,4,5,6,7,42,42,200]
println()
arr = [1, 5, 3]
# looping
for e in arr
    print("$e ")
end # prints 1 5 3

arr = [1, 2, 3]
arr .+ 2 #> [3, 7, 5]
arr * 2 #> [2, 10, 6]
a1 = [1, 2, 3]
a2 = [4, 5, 6]
a1 .* a2 #> [4, 10, 18]

using LinearAlgebra
LinearAlgebra.dot(a1, a2) #> 32
sum(a1 .* a2)

repeat([1, 2, 3], inner = [2]) #> [1,1,2,2,3,3]

a = [1,2,4,6]
a1 = a
show(a1)
a[4] = 0
show(a) #> [1,2,4,0]
show(a1) #> [1,2,4,0]
b = copy(a)
b = deepcopy(a)

a = [1,2,3]
function change_array(arr)
	arr[2] = 25
end
change_array(a)
println(a) #>[ 1, 25, 3]

# the splat operator:
arr = ['a', 'b', 'c']
show(join(arr)) #> "abc"
show(string(arr)) #> "['a', 'b', 'c']"
show(string(arr...)) #> "abc"


text:
module ImageUtilities

export save_image_grid
export image_grid_figure
export image_grid_tensor
export image_grid
export color_image
export color_image_tensor

using FluxGAN: GAN

using Flux: cpu
using CairoMakie
using Images
using FileIO
using Dates
using EllipsisNotation
using Random

###########################
# image utility functions #
###########################

# save image grid figure from model
#
function save_image_grid(model::GAN, output_dir::String;
                         layout=(5, 5),
                         date=false,
                         file_info=[],
                         img_res=150,
                         kws...)
    if file_info == []
        info = ""
    else
        info = mapreduce(tag -> string(tag) * "_", *, file_info)
    end
    file = info * "grid_$(layout[1])_$(layout[2])" *
           (date ? "_" * string(today()) : "")
    imgs = image_grid_tensor(model, *(layout...); kws...)
    fig = image_grid_figure(imgs, layout; img_res=img_res)
    save(output_dir * "/" * file * ".png", fig)
end

# generate and return color images tensor from model
#
function image_grid_tensor(model::GAN, n::Int;
                           uncenter=true,
                           seed=false)
    @assert model.hparams.img_size ≠ undef "Image size not defined!"
    if seed Random.seed!(69) end
    imgs = cpu(model.G)(randn(Float32, model.hparams.latent_dim, n))
    if uncenter imgs = @. (imgs + 1f0) / 2f0 end
    imgs = reshape(imgs, model.hparams.img_size..., :)
    color_image_tensor(imgs)
end

# return image grid makie figure from image grid tensor
#
function image_grid_figure(imgs::Array, layout::Tuple{Int,Int};
                           img_res=150,
                           hflip=true)
    n = size(imgs)[end]
    rows, cols = layout
    @assert rows * cols == n "Number of images ≠ grid size!"

    imgs = [imgs[.., i] for i in 1:n]

    if hflip imgs = map(img -> reverse(img, dims=2), imgs) end

    fig = Figure(resolution=(cols * img_res, rows * img_res))
    for i = 1:rows, j = 1:cols
        ax, = image(fig[i,j], imgs[i + rows * (j - 1)])
        hidedecorations!(ax)
        hidexdecorations!(ax, ticks=false)
        hideydecorations!(ax, ticks=false)
    end
    fig
end

# return image grid from GAN and layout
#
function image_grid(model::GAN, layout::Tuple{Int,Int}; kws...)
    imgs = image_grid_tensor(model, *(layout...); kws...)
    image_grid(imgs, layout)
end

# convert images tensor into an image grid
#
function image_grid(imgs::Array{Color, 3}, layout::Tuple{Int,Int})
    h, w, n = size(imgs)
    rows, cols = layout
    @assert rows * cols == n "Number of images ≠ grid size!"
    grid = Array{Color}(undef, h * rows, w * cols)
    for i = 1:rows, j = 1:cols
        img = imgs[:, :, i + rows*(j - 1)]
        grid[((i - 1)*h + 1):i*h, ((j - 1)*w + 1):j*w] = img'
    end
    grid
end

# return color RGB image from image array
#
function color_image(img::Array)
    @assert length(size(img)) ∈ [2, 3] "Image tensor dimension ∉ [2, 3]!"
    if length(size(img)) == 3
        img = permutedims(img, (3, 1, 2))
        img = Matrix(colorview(RGB, img))
    else
        img = Matrix(colorview(Gray, img))
    end
    img
end

# convert tensor of image data into color image tensor
#
function color_image_tensor(imgs::Array)
    color_imgs = Array{Color}(undef, size(imgs)[1:2]..., size(imgs)[end])
    for i = 1:size(imgs)[end]
        color_imgs[:, :, i] = color_image(imgs[.., i])
    end
    color_imgs
end

end

max_stars_repo_name: aarontrowbridge/GAN

max_issues_repo_name: aarontrowbridge/GAN

max_forks_repo_name: aarontrowbridge/GAN


text:
# # Magnetic liquid droplet in magnetic field

# The behaviour of droplets under the action of magnetic fields of different configurations is an important issue in many domains, including microfluidics (Seeman et al. 2012), mechanics of tissues (Douezan et al. 2011; Frasca et al. 2014), studies of dynamic self-assembly (Timonen et al. 2013) and many others. After the successful synthesis of magnetic fluids in the late sixties by (Rosensweig 1985) an exciting story about droplets of magnetic fluid began. At first, the elongation of the droplets in an external field was observed and explained by Arhipenko et al. (1978). Later on, different droplet configurations were experimentally observed in the high-frequency rotating field, which was subject to our numerical study in [^1].

# Here I present the essential ingredients to reproduce the results of the paper with the new tools SurfaceTopology, LaplaceBIE and ElTopo (optional).

using LinearAlgebra
using GeometryTypes
using SurfaceTopology
using LaplaceBIE
using AbstractPlotting, GLMakie
## using ElTopo

# The calculation requires normal vector and vertex areas methods. Latter one gives the area of 1/3 of the vertex ring, so the sum is the area of the droplet.

function normals(vertices,topology)
    n = Point{3,Float64}[]
    for v in 1:length(vertices)
        s = Point(0,0,0)
        for (v1,v2) in EdgeRing(v,topology)
            s += cross(vertices[v1],vertices[v2])
        end
        normal = s ./ norm(s)
        push!(n,normal)
    end
    return n
end

function vertexareas(points,topology)
    vareas = zeros(Float64,length(points))
    for face in Faces(topology)
        v1,v2,v3 = face
        area = norm(cross(points[v3]-points[v1],points[v2]-points[v1])) /2
        vareas[v1] += area/3
        vareas[v2] += area/3
        vareas[v3] += area/3
    end
    return vareas
end

# These are tools to keep track if the calculation is sensible. Since we are working in an incompressible fluid, the volume needs to be constant, and the energy of the droplet decreases until equilibrium is reached.

function surfacevolume(points,topology)
    normal0 = [0,0,1]
    s = 0
    for face in Faces(topology)
        y1 = points[face[1]]
        y2 = points[face[2]]
        y3 = points[face[3]]

        normaly = cross(y2-y1,y3-y1)
        normaly /= norm(normaly)

        area = norm(cross(y2-y1,y3-y1))/2
        areaproj = dot(normaly,normal0)*area
        volume = dot(y1 + y2 + y3,normal0)/3*areaproj

        s += volume
    end
    return s
end

function energy(points,normals,faces,psi,mup,gammap,H0)
    vareas = vertexareas(points,faces)
    Area = sum(vareas)

    s = 0
    for xkey in 1:length(points)
        s += psi[xkey]*dot(H0,normals[xkey]) * vareas[xkey]
    end

    Es = gammap * Area
    Em = 1/8/pi * (1 - mup) * s

    return Es+Em
end

# For calculating the equilibrium of the droplet, we use a curvatureless algorithm for a vicious droplet developed by Zinchenko 1997. He found a way to calculate velocity generated by a surface tension without explicitly calculating the curvature, which makes it easier to implement and from some tests also more stable. The following method accepts surface defined by points, normals and faces (or topology) and surface force as well as surface tension γ and viscosity η (which only affects the scaling of time).

function stokesvelocity(points,normals,faces,forcen,etaP,gammap)
    vareas = vertexareas(points,faces)
    velocityn = zeros(Float64,length(points))

    for xkey in 1:length(points)

        x = points[xkey]
        nx = normals[xkey]
        fx = forcen[xkey]

        s = 0
        for ykey in 1:length(points)
            if ykey==xkey
                continue
            end

            y = points[ykey]
            ny = normals[ykey]
            fy = forcen[ykey]

            ### Need to check a missing 2
            s += vareas[ykey]*1 ./8/pi/etaP* dot(y-x,nx+ny)/norm(y-x)^3*(1-3*dot(y-x,nx)*dot(y-x,ny)/norm(y-x)^2) * gammap
            s += vareas[ykey]*1 ./8/pi/etaP* ( dot(nx,ny)/norm(x-y) + dot(nx,x -y)*dot(ny,x-y)/norm(x-y)^3 )*(fy - fx)
        end
        velocityn[xkey] = s
    end
    return velocityn
end

# Now having methods defined, we can include a sphere and proceed with calculation.

include("sphere.jl")
msh = unitsphere(2)
vertices, faces = msh.vertices, msh.faces

## Now let's do something fun. Visualize the process with Makie in real time.
x = Node(msh)
y = lift(x->x,x)

scene = Scene(show_axis=false)

wireframe!(scene,y,linewidth = 3f0)
mesh!(scene,y, color = :white, shading = false)

display(scene)

## Initial parameters
H0 = [4.,0.,0.]
etap = 1.
gammap = 1.
μ = 10.

t = 0.
Δt = 0.1
N = 100
volume0 = surfacevolume(vertices,faces)

record(scene, "mdrop.gif", 1:N) do i # for i in 1:N
    n = normals(vertices,faces)

    psi = surfacepotential(vertices,n,faces,μ,H0)
    P∇ψ = tangentderivatives(vertices,n,faces,psi)
    Hn = normalderivatives(vertices,n,faces,P∇ψ,μ,H0)

    E = energy(vertices,n,faces,psi,μ,gammap,H0)
    rV = surfacevolume(vertices,faces)/volume0
    @show E,rV

    Ht = [norm(j) for j in P∇ψ]

    # The force generated (M⋅∇)H, which includes a jump of magnetization at the surface and force coming from the whole bulk.
    tensorn = μ*(μ-1)/8/pi * Hn.^2 + (μ-1)/8/pi * Ht.^2

    vn = stokesvelocity(vertices,n,faces,tensorn,etap,gammap)

    vertices .+= n .* vn * Δt

    msh = HomogenousMesh(vertices,faces)

    ### ElTopo stabilization
    ## par = SurfTrack(allow_vertex_movement=true)
    ## msh = stabilize(msh,par)

    push!(x,msh)
    AbstractPlotting.force_update!()

    global vertices, faces = msh.vertices, msh.faces
    global t += Δt
end

# ![](mdrop.gif)

# [^1]: Erdmanis, J. & Kitenbergs, G. & Perzynski, R. & Cebers, A. (2017) Magnetic micro-droplet in rotating field: numerical simulation and comparison with experiment


text:
# This script will usually fail as Stan gives up when the integration fails.

using StanSample, StatsPlots

ProjDir = @__DIR__

ben_model = "
functions {
  real[] dz_dt(real t,       // time
               real[] z,     // system state {prey, predator}
               real[] theta, // parameters
               real[] x_r,   // unused data
               int[] x_i) {
    real u = z[1];
    real v = z[2];

    real alpha = theta[1];
    real beta = theta[2];
    real gamma = theta[3];
    real delta = theta[4];

    real du_dt = (alpha - beta * v) * u;
    real dv_dt = (-gamma + delta * u) * v;

    return { du_dt, dv_dt };
  }
}

data {
  int<lower = 0> N;
  real t0;
  real ts[N];
  real y_init[2];
  real y[N, 2];
}

transformed data {
  real beta = 1.0;
  real gamma = 3.0;
  real delta = 1.0;
}

parameters {
  real alpha;
  real z_init[2];
  real<lower = 0> sigma[2];
}

transformed parameters {
  real theta[4] = { alpha, beta, gamma, delta };
  real z[N, 2] = integrate_ode_rk45(dz_dt, z_init, t0, ts, theta,
                       rep_array(0.0, 0), rep_array(0, 0),
                       1e-5, 1e-3, 500);
}

model {
  alpha ~ normal(1.5, 0.01);

  sigma ~ inv_gamma(3, 3);
  z_init ~ normal(1, 0.01);

  for (k in 1:2) {
    y_init[k] ~ normal(z_init[k], sigma[k]);
    y[ , k] ~ normal(z[, k], sigma[k]);
  }
}
";

sm = SampleModel("Ben", ben_model, tmpdir="$(@__DIR__)/tmp")

obs = reshape([
    2.75487369287842, 6.77861583902452, 0.977380864533721,
    1.87996617432547,6.10484970865812, 1.38241825210688, 1.32745986436496,
    4.35128045090906, 3.28082272129074, 1.02627964892546,
    0.252582832985904, 2.01848763447693, 1.91654508055532,
    0.33796473710305, 0.621805831293324, 3.47392930294277, 0.512736737236118,
    0.304559407399392, 4.56901766791213, 0.91966250277948], (10,2))

ben_data = Dict(
  :N => 10,
  :t0 => 0,
  :ts => 1:10,
  :y => obs,
  :y_init => [1.01789244983237, 0.994539126829031]
)

t = ben_data[:ts]

plot(xlab="t", ylab="prey and pred")
plot!(t, obs[:,1], lab="prey")
plot!(t, obs[:, 2], lab="pred")
savefig("$(ProjDir)/observations.png")

@time rc = stan_sample(sm; data=ben_data)

if success(rc)
  p = read_samples(sm, output_format=:particles)
  display(p)
end


text:
using Flux

struct NeuralNetworkQ{Tm, To, Tp}  <: AbstractQApproximator{Any, Int}
    model::Tm
    opt::To
    ps::Tp
    function NeuralNetworkQ(model::Tm, opt::To) where {Tm, To}
        m = gpu(model)
        ps = params(m)
        new{typeof(m), To, typeof(ps)}(m, opt, ps)
    end
end

(Q::NeuralNetworkQ)(s, ::Val{:dist}) = Q.model(s)
(Q::NeuralNetworkQ)(s) = Q(s, Val(:dist))
(Q::NeuralNetworkQ)(s, ::Val{:argmax}) = reshape(map(i -> i[1], findmax(Q(s).data, dims=1)[2]), :)
(Q::NeuralNetworkQ)(s, ::Val{:max}) = dropdims(maximum(Q(s), dims=1), dims=1)
(Q::NeuralNetworkQ)(s, a::Int) = Q(s)[a]

function (Q::NeuralNetworkQ)(s, a::AbstractArray{Int, 1})
    dist = Q(s)
    inds = CartesianIndex.(a, axes(dist, 2))
    dist[inds]
end

function update!(Q::NeuralNetworkQ, loss)
    Flux.back!(loss)
    Flux.Optimise.update!(Q.opt, Q.ps)
end


text:
module AnnDatas

using HDF5
using SparseArrays
using DataFrames

export AnnData

struct AnnData
    X::SparseMatrixCSC
    obsm::Union{Nothing, Dict{String, Any}}
    obsp::Union{Nothing, Dict{String, Any}}
    uns::Union{Nothing, Dict{String, Any}}
    obs::Union{Nothing, DataFrame}
    var::Union{Nothing, DataFrame}
end

"""
Read a CSR matrix in a SparseMatrixCSC
"""
function read_csr_matrix(g::HDF5.Group)
    attr = attributes(g)
    @assert read(attr["encoding-type"]) == "csr_matrix"
    m, n = read(attr["shape"])

    V           = read(g["data"])
    csr_indices = read(g["indices"]) .+ 1
    csr_indptr  = read(g["indptr"]) .+ 1
    nnz = length(V)

    # easiest way to do this is to go CSR -> COO -> CSC
    I = Vector{Int32}(undef, nnz)
    J = Vector{Int32}(undef, nnz)

    for i in 1:m
        for k in csr_indptr[i]:csr_indptr[i+1]-1
            I[k] = i
            J[k] = csr_indices[k]
        end
    end

    return sparse(I, J, V, m, n)
end

"""
Read a serialized data frame into a DataFrame
"""
function read_dataframe(input::HDF5.File, path::String)
    if !haskey(input, path)
        return nothing
    end
    g = input[path]

    columns = Dict{String, Any}()
    attr = attributes(g)
    @assert read(attr["encoding-type"]) == "dataframe"
    columnorder = read(attr["column-order"])

    for key in keys(g)
        columns[key] = read(g[key])
    end
    df = DataFrame(Dict(key => columns[key] for key in columnorder))
    df[!,"_index"] = read(g["_index"])

    return df
end

"""
General purpose function to read a groups into a dictionary tree structure.
"""
function read_group(input::HDF5.File, path::String)
    if !haskey(input, path)
        return nothing
    end
    g = input[path]

    data = Dict{String, Any}()
    for key in keys(g)
        dataset = g[key]

        if isa(dataset, HDF5.Group)
            attr = attributes(dataset)
            if length(attr) == 0
                # data[key] = read_group(dataset)
                data[key] = read(dataset)
            elseif haskey(attr, "encoding-type") && read(attr["encoding-type"]) == "csr_matrix"
                data[key] = read_csr_matrix(dataset)
            else
                data[key] = read(dataset)
            end
            # TODO: special cases for other types of data
        else
            data[key] = read(dataset)
        end
    end

    return data
end

"""
Read an AnnData struct from the given h5ad filename.
"""
function Base.read(filename::AbstractString, ::Type{AnnData})
    input = h5open(filename)

    if isa(input["X"], HDF5.Group)
        X = read_csr_matrix(input["X"])
    else
        X = read(input["X"])
    end
    obsm = read_group(input, "obsm")
    obsp = read_group(input, "obsp")
    uns  = read_group(input, "uns")
    obs  = read_dataframe(input, "obs")
    var  = read_dataframe(input, "var")

    close(input)

    return AnnData(X, obsm, obsp, uns, obs, var)
end

end # module


text:
module BIFMHelperNodeTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "BIFMHelperNode" begin

    @testset "Creation" begin
        node = make_node(BIFMHelper)

        @test functionalform(node)          === BIFMHelper
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :in)
        @test factorisation(node)           === ((1, 2), )

        @test length(methods(functional_dependencies, (Any, FactorNode{Type{BIFMHelper}}, Int))) === 1

    end

    @testset "Average energy" begin

        node = make_node(BIFMHelper)

        @test score(AverageEnergy(),
                    BIFMHelper,
                    Val{(:out, :in)},
                    (Marginal(MvNormalMeanCovariance([1,1], [2 0; 0 3]), false, false), Marginal(MvNormalMeanCovariance([1,1], [2 0; 0 3]), false, false)),
                    nothing) ≈ entropy(MvNormalMeanCovariance([1,1], [2 0; 0 3]))

        @test score(AverageEnergy(),
                    BIFMHelper,
                    Val{(:out, :in)},
                    (Marginal(MvNormalMeanCovariance([1,2], [2 0; 0 1]), false, false), Marginal(MvNormalMeanPrecision([1,2], [0.5 0; 0 1]), false, false)),
                    nothing) ≈ entropy(MvNormalMeanCovariance([1,2], [2 0; 0 1]))

    end

end

end


text:
__precompile__()

module MicroLogging

# ----- Core API to go in Base -----
export
    ## Logger types
    AbstractLogger, LogLevel, NullLogger,
    # Public logger API:
    #   handle_message, shouldlog, min_enabled_level, catch_exceptions
    # (not exported, as they're not generally called by users)
    #
    ## Log creation
    @debug, @info, @warn, @error, @logmsg,
    ## Logger installation and control
    # TODO: Should some of these go into stdlib ?
    with_logger, current_logger, global_logger, disable_logging

# ----- API to go in StdLib package ? -----
export
    SimpleLogger, # Possibly needed in Base?
    # TODO: configure_logging needs a big rethink (see, eg, python's logger
    # config system)
    configure_logging

# ----- MicroLogging stuff, for now -----
export
    InteractiveLogger

# core.jl includes the code which will hopefully go into Base in 0.7
const core_in_base = isdefined(Base, :CoreLogging)

if core_in_base
    import Logging:
        @debug, @info, @warn, @error, @logmsg,
        AbstractLogger,
        LogLevel, BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
        with_logger, current_logger, global_logger, disable_logging,
        handle_message, shouldlog, min_enabled_level, catch_exceptions,
        SimpleLogger
else
    include("core.jl")
end

include("loggers.jl")

include("config.jl")

if !core_in_base
    include("test.jl")
end

function __init__()
    global_logger(InteractiveLogger(STDERR))
end

end


text:
using Stheno, Test, Random, LinearAlgebra, Statistics

@testset "Stheno" begin

    @testset "util" begin
        include("util/covariance_matrices.jl")
        include("util/woodbury.jl")
        include("util/block_arrays.jl")
        include("util/abstract_data_set.jl")
    end

    @testset "mean_and_kernel" begin
        include("test_util.jl")
        include("mean_and_kernel/mean.jl")
        include("mean_and_kernel/kernel.jl")
        include("mean_and_kernel/compose.jl")
        include("mean_and_kernel/finite.jl")
        include("mean_and_kernel/conditional.jl")
        include("mean_and_kernel/block.jl")
        # include("mean_and_kernel/transform.jl")
        include("mean_and_kernel/input_transform.jl")
        include("mean_and_kernel/zero.jl")
        include("mean_and_kernel/degenerate.jl")
        include("mean_and_kernel/derivative.jl")
        include("mean_and_kernel/conversion.jl")
    end

    @testset "gp" begin
        include("gp/abstract_gp.jl")
        include("gp/gp.jl")
        include("gp/block_gp.jl")
    end

    @testset "linops" begin
        include("linops/indexing.jl")
        include("linops/addition.jl")
        include("linops/product.jl")
        # include("linops/integrate.jl")
        include("linops/project.jl")
        include("linops/conditioning.jl")
    end

    @testset "integration" begin
        include("util/toeplitz_integration.jl")
    end
end


text:
"""
Get the curl of a vector f w.r.t Ds
"""
function curl(vec, Ds)
    [Ds[2](vec[3]) - Ds[3](vec[2]), Ds[3](vec[1]) - Ds[1](vec[3]), Ds[1](vec[2]) - Ds[2](vec[1])]
end

"""
Get the divergence of a function f w.r.t Ds with the option of multiplying each part of the sum by a v
"""
function divergence(Ds, f, v=ones(length(Ds)))
    if f isa AbstractArray
        sum([v[i] * Ds[i](f[i]) for i in eachindex(Ds)])
    else
        sum([v[i] * Ds[i](f) for i in eachindex(Ds)])
    end
end

"""
Print the loss of the loss function
"""

function print_loss(timeCounter, startTime, times, losses)
    deltaT_s = time_ns() #Start a clock when the callback begins
    timeCounter = timeCounter + time_ns() - deltaT_s
    ctime = time_ns() - startTime - timeCounter #This variable is the time to use for the time benchmark plot
    append!(times, ctime/10^9) #Conversion nanosec to seconds

    function cb(p, l)
        append!(losses, l)
        println("Current loss is: $l")
        return false
    end
end

"""
Simulates a collisionless plasma given:

plasma – Struct of type CollisionlessPlasma
lb – lower bound
ub – upper bound
time_lb – lower time boundary
time_ub – upper time boundary
GPU – whether this should be executed on the GPU or not
inner_layers – how many layers in the neural network
strategy – what NeuralPDE training strategy should be used
"""
function solve(plasma::CollisionlessPlasma;
               lb=0.0, ub=1.0, time_lb=lb, time_ub=ub,
               GPU=true, inner_layers=16, strategy=StochasticTraining(100),
               E_bcs=Neumann,
               f_bcs=(a=1, g=0), maxiters=[30,100], normalized=true)
    if lb > ub
        error("lower bound must be larger than upper bound")
    end

    if GPU && strategy == QuadratureTraining()
        error("QuadratureTraining does not have GPU support. Use another strategy (e.g. StochasticTraining(200)) instead")
    end

    losses = []
    times = []
    timeCounter = 0.0
    startTime = time_ns()

    # constants
    dim = 3
    geometry = plasma.geometry.f # this might change with a geometry refactor
    dis = plasma.distributions
    species = [d.species for d in dis]
    consts = Constants()
    μ_0, ϵ_0 = consts.μ_0, consts.ϵ_0

    # get qs, ms, Ps from species
    qs, ms = Float64[], Float64[]
    for s in species
        push!(qs,s.q)
        push!(ms,s.m)
    end
    Ps = [d.P for d in dis]

    if normalized
        qs .= 1
        ms .= 1
        ϵ_0 = 1
        μ_0 = 0
    end

    # variables
    fs = Symbolics.variables(:f, eachindex(species); T=SymbolicUtils.FnType{Tuple,Real})
    Es = Symbolics.variables(:E, 1:dim; T=SymbolicUtils.FnType{Tuple,Real})
    Bs = Symbolics.variables(:B, 1:dim; T=SymbolicUtils.FnType{Tuple,Real})

    # parameters
    @parameters t
    xs,vs = Symbolics.variables(:x, 1:dim), Symbolics.variables(:v, 1:dim)

    # integrals
    _I = Integral(tuple(vs...) in DomainSets.ProductDomain(ClosedInterval(-Inf ,Inf), ClosedInterval(-Inf ,Inf), ClosedInterval(-Inf ,Inf)))

    # differentials
    Dxs = Differential.(xs)
    Dvs = Differential.(vs)
    Dt = Differential(t)

    # domains
    xs_int = xs .∈ Interval(lb, ub)
    vs_int = vs .∈ Interval(lb, ub)
    t_int = t ∈ Interval(time_lb, time_ub)

    domains = [t_int;xs_int;vs_int]

    # helpers
    _Es = [E(t,xs...) for E in Es]
    _Bs = [B(t,xs...) for B in Bs]
    _fs = [f(t,xs...,vs...) for f in fs]

    # divergences
    div_vs = [divergence(Dxs, _f, vs) for _f in _fs]
    div_B = divergence(Dxs, _Bs)
    div_E = divergence(Dxs, _Es)
    Fs = [qs[i]/ms[i] * (_Es + cross(vs,_Bs)) for i in eachindex(qs)]
    divv_Fs = [divergence(Dvs, _fs[i], Fs[i]) for i in eachindex(_fs)]

    # charge and current densities
    ρ = sum([qs[i] * _I(_fs[i]) for i in eachindex(qs)])
    J = sum([qs[i] * _I(_fs[i]) * vs[j] for i in 1:length(_fs), j in 1:length(vs)])

    # system of equations
    vlasov_eqs = Dt.(_fs) .~ .- div_vs .- divv_Fs
    curl_E_eqs = curl(_Es, Dxs) .~ Dt.(_Bs)
    curl_B_eqs = ϵ_0*μ_0 * Dt.(_Es) .- curl(_Bs, Dxs) .~ - μ_0.*J
    div_E_eq = div_E ~ ρ/ϵ_0
    div_B_eq = div_B ~ 0
    eqs = [vlasov_eqs; curl_E_eqs; curl_B_eqs; div_E_eq; div_B_eq]

    # boundary conditions
    div_E0 = sum([Dxs[i](Es[i](time_lb, xs...)) for i in eachindex(Dxs)])
    div_B0 = sum([Dxs[i](Bs[i](time_lb, xs...)) for i in eachindex(Dxs)])

    vlasov_ics = [fs[i](time_lb,xs...,vs...) ~ Ps[i](xs,vs) * geometry(xs) for i in eachindex(fs)]
    div_B_ic = div_B0 ~ 0
    div_E_ic =  div_E0 ~ sum([qs[i] * _I(fs[i](time_lb,xs...,vs...)) for i in eachindex(qs)])/ϵ_0 * geometry(xs)
    E_bcs = [E_bcs(t, xs, Dxs, Es, lb); E_bcs(t, xs, Dxs, Es, ub)]
    f_bcs = [Reflective(t, xs, vs, fs, lb, f_bcs.a, f_bcs.g); Reflective(t, xs, vs, fs, ub, f_bcs.a, f_bcs.g)]
    bcs = [vlasov_ics;div_B_ic; div_E_ic; E_bcs; f_bcs]

    # get variables
    vars_arg = [_fs...; _Es...; _Bs...]
    vars = [fs, Bs, Es]

    dict_vars = Dict()
    for var in vars
        push!(dict_vars, var => [v for v in var])
    end

    # set up PDE System
    @named pde_system = PDESystem(eqs, bcs, domains, [t,xs...,vs...], vars_arg)

    # set up problem
    il = inner_layers
    ps_chains = [FastChain(FastDense(length(domains), il, Flux.σ), FastDense(il,il,Flux.σ), FastDense(il, 1)) for _ in 1:length(_fs)]
    xs_chains = [FastChain(FastDense(length([t, xs...]), il, Flux.σ), FastDense(il,il,Flux.σ), FastDense(il, 1)) for _ in 1:length([_Es; _Bs])]
    chain = [ps_chains;xs_chains]
    initθ = GPU ? map(c -> CuArray(Float64.(c)), DiffEqFlux.initial_params.(chain)) : map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
    discretization = NeuralPDE.PhysicsInformedNN(chain, strategy, init_params=initθ)
    prob = SciMLBase.discretize(pde_system, discretization)

    # solve
    opt = Optim.BFGS()
    res = GalacticOptim.solve(prob, opt, cb = print_loss(timeCounter, startTime, times, losses), maxiters=maxiters[1])
    prob = remake(prob, u0=res.minimizer)
    res = GalacticOptim.solve(prob, ADAM(0.01), cb = print_loss(timeCounter, startTime, times, losses), maxiters=maxiters[2])
    phi = discretization.phi

    return PlasmaSolution(plasma, vars, dict_vars, phi, res, initθ, domains, losses, times)
end

"""
Simulates an electrostatic plasma given:

plasma – Struct of type ElectrostaticPlasma
dim – in how many D and V dimensions should the simulation happen
lb – lower bound
ub – upper bound
time_lb – lower time boundary
time_ub – upper time boundary
GPU – whether this should be executed on the GPU or not
inner_layers – how many layers in the neural network
strategy – what NeuralPDE training strategy should be used
"""
function solve(plasma::ElectrostaticPlasma;
    lb=0.0, ub=1.0, time_lb=lb, time_ub=ub,
    dim=3, GPU=true, inner_layers=16, strategy=StochasticTraining(80*dim),
    E_bcs=Neumann,f_bcs=(a=1, g=0), maxiters=[30,100], normalized=true)

    if lb > ub
        error("lower bound must be larger than upper bound")
    end

    if GPU && strategy == QuadratureTraining()
        error("QuadratureTraining does not have GPU support. Use another strategy (e.g. StochasticTraining(200)) instead")
    end

    losses = []
    times = []
    timeCounter = 0.0
    startTime = time_ns()

    # constants
    geometry = plasma.geometry.f # this might change with a geometry refactor
    dis = plasma.distributions
    species = [d.species for d in dis]
    consts = Constants()
    ϵ_0 = consts.ϵ_0

    # get qs, ms, Ps from species
    qs, ms = Float64[], Float64[]
    for s in species
        push!(qs,s.q)
        push!(ms,s.m)
    end
    Ps = [d.P for d in dis]

    if normalized
        qs .= 1
        ms .= 1
        ϵ_0 = 1
    end

    # variables
    fs = Symbolics.variables(:f, eachindex(species); T=SymbolicUtils.FnType{Tuple,Real})
    Es = Symbolics.variables(:E, 1:dim; T=SymbolicUtils.FnType{Tuple,Real})

    # parameters
    @parameters t
    xs,vs = Symbolics.variables(:x, 1:dim), Symbolics.variables(:v, 1:dim)

    # integrals
    _I = if length(vs) > 1
        intervals = [ClosedInterval(-Inf ,Inf) for _ in 1:length(vs)]
        _I = Integral(tuple(vs...) in DomainSets.ProductDomain(intervals...))
    else
        _I = Integral(first(vs) in DomainSets.ClosedInterval(-Inf ,Inf))
    end

    # differentials
    Dxs = Differential.(xs)
    Dvs = Differential.(vs)
    Dt = Differential(t)

    # domains
    xs_int = xs .∈ Interval(lb, ub)
    vs_int = vs .∈ Interval(lb, ub)
    t_int = t ∈ Interval(time_lb, time_ub)

    domains = [t_int;xs_int;vs_int]

    # helpers
    _Es = [E(t,xs...) for E in Es]
    _fs = [f(t,xs...,vs...) for f in fs]

    # divergences
    div_vs = [divergence(Dxs, _f, vs) for _f in _fs]
    div_E = divergence(Dxs, _Es)
    Fs = [qs[i]/ms[i] * _Es for i in eachindex(qs)]
    divv_Fs = [divergence(Dvs, _fs[i], Fs[i]) for i in eachindex(_fs)]

    # charge density
    ρ = sum([qs[i] * _I(_fs[i]) for i in eachindex(qs)])

    # equations
    vlasov_eqs = Dt.(_fs) .~ .- div_vs .- divv_Fs
    div_E_eq = div_E ~ ρ/ϵ_0
    eqs = [vlasov_eqs; div_E_eq]

    # boundary and initial conditions
    div_E0 = sum([Dxs[i](Es[i](time_lb, xs...)) for i in eachindex(Dxs)])

    vlasov_ics = [fs[i](time_lb,xs...,vs...) ~ Ps[i](xs,vs) * geometry(xs) for i in eachindex(fs)]
    div_E_ic = div_E0 ~ sum([qs[i] * _I(fs[i](time_lb,xs...,vs...)) for i in eachindex(qs)])/ϵ_0 * geometry(xs)

    f_bcs = [Reflective(t, xs, vs, fs, lb, f_bcs.a, f_bcs.g); Reflective(t, xs, vs, fs, ub, f_bcs.a, f_bcs.g)]
    E_bcs = [E_bcs(t, xs, Dxs, Es, lb); E_bcs(t, xs, Dxs, Es, ub)]
    bcs = [vlasov_ics; div_E_ic; E_bcs; f_bcs]

    # set up variables
    vars_arg = [_fs; _Es]
    vars = [fs, Es]

    dict_vars = Dict()
    for var in vars
        push!(dict_vars, var => [v for v in var])
    end

    # set up and return PDE System
    @named pde_system = PDESystem(eqs, bcs, domains, [t,xs...,vs...], vars_arg)

    # set up problem
    il = inner_layers
    ps_chains = [FastChain(FastDense(length(domains), il, Flux.σ), FastDense(il,il,Flux.σ), FastDense(il, 1)) for _ in 1:length(_fs)]
    xs_chains = [FastChain(FastDense(length([t, xs...]), il, Flux.σ), FastDense(il,il,Flux.σ), FastDense(il, 1)) for _ in 1:length(_Es)]
    chain = [ps_chains;xs_chains]
    initθ = GPU ? map(c -> CuArray(Float64.(c)), DiffEqFlux.initial_params.(chain)) : map(c -> Float64.(c), DiffEqFlux.initial_params.(chain))
    discretization = NeuralPDE.PhysicsInformedNN(chain, strategy, init_params=initθ)
    prob = SciMLBase.discretize(pde_system, discretization)

    # solve
    opt = Optim.BFGS()
    res = GalacticOptim.solve(prob, opt, cb = print_loss(timeCounter, startTime, times, losses), maxiters=maxiters[1])
    prob = remake(prob, u0=res.minimizer)
    res = GalacticOptim.solve(prob, ADAM(0.01), cb = print_loss(timeCounter, startTime, times, losses), maxiters=maxiters[2])
    phi = discretization.phi

    return PlasmaSolution(plasma, vars, dict_vars, phi, res, initθ, domains, losses, times)
end


text:
struct BoundingBox
    ramin::Float64
    ramax::Float64
    decmin::Float64
    decmax::Float64

    function BoundingBox(ramin::Float64, ramax::Float64,
                         decmin::Float64, decmax::Float64)
        @assert ramax > ramin "ramax must be greater than ramin"
        @assert decmax > decmin "decmax must be greater than decmin"
        new(ramin, ramax, decmin, decmax)
    end
end

function BoundingBox(ramin::String, ramax::String,
                     decmin::String, decmax::String)
    BoundingBox(parse(Float64, ramin),
                parse(Float64, ramax),
                parse(Float64, decmin),
                parse(Float64, decmax))
end

"""
    SurveyDataSet

Abstract type representing a collection of imaging data and associated
metadata from a survey. Concrete subtypes are for specific
surveys. They provide methods for querying which images are available
and loading Celeste Image objects from disk. In short, they provide an
interface between survey data organized on disk and in-memory Celeste
objects.
"""
abstract type SurveyDataSet end

load_images(::SurveyDataSet, box::BoundingBox) =
    error("load_images() not defined for this SurveyDataSet type")