text:
using CoordinateTransformations
using Dierckx
using HCIToolbox
using ImageTransformations: center
using LinearAlgebra: dot
using Photometry
using ProgressLogging
using StaticKernels

"""
    contrast_curve(alg, cube, angles, psf;
                   fwhm, sigma=5, nbranch=1, theta=0, inner_rad=1,
                   starphot=Metrics.estimate_starphot(cube, fwhm),
                   fc_rad_sep=3, snr=100, k=2, smooth=true,
                   subsample=true, kwargs...)

Calculate the throughput-calibrated contrast. This first processes the algorithmic [`throughput`](@ref) by injecting instances of `psf` into `cube`. These are processed through `alg` and the ratio of the recovered flux to the injected flux is calculated. These companions are injected in resolution elements across the frame, which can be changed via the various keyword arguments.

The throughput can only be calculated for discrete resolution elements, but we typically want a much smoother curve. To accomplish this, we measure the noise (the standard deviation of all resolution elements in an annulus at a given radius) for every pixel in increasing radii. We then interpolate the throughput to this grid and return the subsampled curves.

# Returned Fields
* `distance` - The radial distance (in pixels) for each measurement
* `contrast` - The Gaussian sensitivity
* `contrast_corr` - The Student-t sensitivity
* `noise` - The noise measured for each distance
* `throughput` - The throughput measured for each distance.

# Keyword Arguments
* `sigma` - The confidence limit in terms of Gaussian standard deviations
* `starphot` - The flux of the star. By default calculates the flux in the central core.
**Injection Options** (See also [`throughput`](@ref))
* `nbranch` - number of azimuthal branches to use
* `theta` - position angle of initial branch
* `inner_rad` - position of innermost planet in FWHM
* `fc_rad_sep` - the separation between planets in FWHM for a single reduction
* `snr` - the target signal to noise ratio of the injected planets
**Subsampling Options** (See also [`Metrics.subsample_contrast`](@ref))
* `subsample` - If true, subsamples the throughput measurements to increase density of curve
* `k` - The order of the BSpline used for subsampling the throughput
* `smooth` - If true, will smooth the subsampled noise measurements with a 2nd order Savitzky-Golay filter

Any additional `kwargs` will be passed to `alg` when it is called.

!!! tip
    If you prefer a tabular format, simply pipe the output of this function into any type supporting the Tables.jl interface, e.g.
    ```
    contrast_curve(alg, cube, angles, psf; fwhm=fwhm) |> DataFrame
    ```
"""
function contrast_curve(alg, cube, angles, psf, args...;
        fwhm, sigma=5, nbranch=1, theta=0, inner_rad=1,
        starphot=Metrics.estimate_starphot(cube, fwhm),
        fc_rad_sep=3, snr=100, k=2, smooth=true,
        subsample=true, kwargs...)
    starphot > 0 || @warn "starphot should be positive but was $starphot"
    # measure the noise and throughput in consecutive resolution elements
    # across azimuthal branches
    @info "Calculating Throughput"
    reduced_empty = alg(cube, angles; kwargs...)

    through, meta = throughput(alg, cube, angles, psf, args...;
                               fwhm=fwhm, inner_rad=inner_rad, fc_rad_sep=fc_rad_sep, theta=theta,
                               nbranch=nbranch, snr=snr, reduced_empty=reduced_empty, kwargs...)

    through_mean = mean(through, dims=2) |> vec

    if subsample
        return subsample_contrast(reduced_empty, meta.distance, through_mean;
                                  fwhm=fwhm, starphot=starphot, sigma=sigma, inner_rad=inner_rad,
                                  theta=theta, smooth=smooth, k=k)
    end

    # calculate common terms once
    unit_contrast = @. meta.noise / (through_mean * starphot)
    # gaussian contrast (invalid values become NaN)
    contrast = calculate_contrast.(sigma, unit_contrast)

    # get correction for small-sample statistics
    sigma_corr = correction_factor.(meta.distance, fwhm, sigma)
    # student-t contrast (invalid values become NaN)
    contrast_corr = calculate_contrast.(sigma_corr, unit_contrast)

    return (distance=meta.distance,
            throughput=through_mean,
            contrast=contrast,
            contrast_corr=contrast_corr,
            noise=meta.noise)
end

function correction_factor(radius, fwhm, sigma)
    n_res_els = 2 * π * radius ÷ fwhm
    ss_corr = sqrt(1 + 1 / (n_res_els - 1))
    return quantile(TDist(n_res_els), cdf(Normal(), sigma)) * ss_corr
end

"""
    Metrics.subsample_contrast(empty_frame, distance, throughput;
                               fwhm, starphot, sigma=5, inner_rad=1,
                               theta=0, smooth=true, k=2)

Helper function to subsample and smooth a contrast curve.

Contrast curves, by definition, are calculated with discrete resolution elements. This can cause contrast curves to have very few points instead of appearing as a continuously measured statistic across the pixels. We alleviate this by sub-sampling the throughput (via BSpline interpolation) across each pixel (instead of each resolution element).

The noise can be found efficiently enough, so rather than interpolate we measure the noise in annuli of width `fwhm` increasing in distance by 1 pixel. We measure this noise in `empty_frame`, which should be a 2D reduced ADI frame.

The noise measurements can be noisy, so a 2nd order Savitzky-Golay filter can be applied via `smooth`. This fits a quadratic polynomial over a window of `fwhm/2` points together to reduce high-frequency jitter.

# Examples

Here is an example which calculates the exact contrast curve in addition to a subsampled curve without re-calculating the throughput.

```julia
cube, angles, psf = # load data

alg = PCA(10)
cc = contrast_curve(alg, cube, angles, psf; fwhm=8.4, subsample=false)
reduced_empty = alg(cube, angles)
sp = Metrics.estimate_starphot(cube, 8.4)
cc_sub = Metrics.subsample_contrast(reduced_empty, cc.distance, cc.throughput;
                                    fwhm=8.4, starphot=sp)
```
"""
function subsample_contrast(empty_frame, distance, throughput;
                            fwhm, starphot, sigma=5, inner_rad=1,
                            theta=0, smooth=true, k=2)
    # measure the noise with high sub-sampling-
    # at every pixel instead of every resolution element
    cy, cx = center(empty_frame)
    radii_subsample = first(distance):1:last(distance) + 1
    through_subsample = Spline1D(distance, throughput, k=k)(radii_subsample)
    noise_subsample = @. annulus_noise((empty_frame,), fwhm, cy, cx, radii_subsample, theta)

    if smooth
        window_size = min(length(noise_subsample) - 2, round(Int, 2 * fwhm))
        iseven(window_size) && (window_size += 1)
        width = window_size ÷ 2
        coeffs = savgol_coeffs(window_size, 2) |> Tuple
        smooth_kernel = Kernel{(-width:width,)}(w -> dot(coeffs, Tuple(w)))
        noise_smoothed = map(smooth_kernel, extend(noise_subsample, StaticKernels.ExtensionReplicate()))
    else
        noise_smoothed = noise_subsample
    end

    # calculate common terms once
    unit_contrast = @. noise_smoothed / (through_subsample * starphot)
    # gaussian contrast (invalid values become NaN)
    contrast = calculate_contrast.(sigma, unit_contrast)

    # get correction for small-sample statistics
    sigma_corr = correction_factor.(radii_subsample, fwhm, sigma)
    # student-t contrast (invalid values become NaN)
    contrast_corr = calculate_contrast.(sigma_corr, unit_contrast)

    return (distance=radii_subsample,
            throughput=through_subsample,
            contrast=contrast,
            contrast_corr=contrast_corr,
            noise=noise_smoothed)
end

# simple function to inline the contrast calculation _and_ clipping
@inline function calculate_contrast(k, unit_contrast)
    contrast = k * unit_contrast
    return ifelse(0 ≤ contrast ≤ 1, contrast, oftype(contrast, NaN))
end

function savgol_coeffs(window_size, order=2)
    pos = window_size ÷ 2
    # form vandermonde matrix
    x = -pos:window_size - pos - 1
    powers = 0:order
    A = x' .^ powers
    # solve least squares equation Ac = y
    y = zeros(order + 1)
    y[1] = 1
    return A \ y
end

"""
    Metrics.estimate_starphot(cube, fwhm)
    Metrics.estimate_starphot(frame, fwhm)

Simple utility to estimate the stellar photometry by placing a circular aperture with `fwhm` diameter in the center of the `frame`. If a cube is provided, first the median frame will be found.
"""
function estimate_starphot(frame::AbstractMatrix, fwhm)
    ap = CircularAperture(reverse(center(frame)), fwhm/2)
    return photometry(ap, frame).aperture_sum
end

estimate_starphot(cube::AbstractArray{T,3}, fwhm) where {T} = estimate_starphot(collapse(cube, method=median), fwhm)

"""
    throughput(alg, cube, angles, psf;
               fwhm, nbranch=1, theta=0, inner_rad=1,
               fc_rad_sep=3, snr=100, kwargs...)

Calculate the throughput of `alg` by injecting fake companions into `cube` and measuring the relative photometry of each companion in the reduced frame. The photometry is measured using a circular aperture with a diameter matching the `fwhm`. Any additional `kwargs` will be passed to `alg` when it is called.

# Keyword Arguments
* `nbranch` - number of azimuthal branches to use
* `theta` - position angle of initial branch
* `inner_rad` - position of innermost planet in FWHM
* `fc_rad_sep` - the separation between planets in FWHM for a single reduction
* `snr` - the target signal to noise ratio of the injected planets
* `reduced_empty` - the collapsed residual frame for estimating the noise. Will process using `alg` if not provided.
"""
function throughput(alg, cube::AbstractArray{T,3}, angles, psf_model;
                    fwhm, nbranch=1, theta=0, inner_rad=1, fc_rad_sep=3,
                    snr=100, reduced_empty = nothing, kwargs...) where T
    maxfcsep = size(cube, 2) ÷ (2 * fwhm) - 1
    # too large separation between companions in the radial patterns
    3 ≤ fc_rad_sep ≤ maxfcsep || error("`fc_rad_sep` should lie ∈[3, $(maxfcsep)], got $fc_rad_sep")

    # compute noise in concentric annuli on the empty frame
    if isnothing(reduced_empty)
        reduced_empty = alg(cube, angles; kwargs...)
    end

    cy, cx = center(reduced_empty)

    n_annuli = floor(Int, (cy - fwhm) / fwhm) - 1
    radii = _fix_range(fwhm .* (inner_rad:n_annuli), cube)

    δy, δx = sincosd(theta)
    noise = @. annulus_noise((reduced_empty,), fwhm, cy, cx, radii, theta)

    angle_per_branch = 360 / nbranch
    output = similar(cube, length(radii), nbranch)

    fake_comps_full = zero(reduced_empty)
    @progress name="branch" for branch in 1:nbranch
        θ = theta + angle_per_branch * (branch - 1)
        @progress name="pattern" for init_rad in 1:fc_rad_sep
            slice = init_rad:fc_rad_sep:lastindex(radii)
            fake_comps = zero(reduced_empty)

            cube_fake_comps = copy(cube)

            apertures = map(slice) do ann
                r = radii[ann]
                δy, δx = sincosd(θ)
                x = r * δx + cx
                y = r * δy + cy

                A = snr * noise[ann]

                inject!(fake_comps, psf_model, Polar(r, deg2rad(θ)); A=A)
                fake_comps_full .+= fake_comps
                inject!(cube_fake_comps, angles, psf_model, Polar(r, deg2rad(θ)); A=A)

                return CircularAperture(x, y, fwhm / 2)
            end
            # make sure geometries are maintained
            reduced = alg(cube_fake_comps, angles; kwargs...)

            injected_flux = photometry(apertures, fake_comps).aperture_sum
            recovered_flux = photometry(apertures, reduced .- reduced_empty).aperture_sum
            @. output[slice, branch] = max(zero(T), recovered_flux / injected_flux)
        end
    end

    return output, (distance=radii, fake_comps=fake_comps_full, noise=noise)
end

_fix_range(radii, cube::AbstractArray{T,3}) where {T} = radii
_fix_range(radii, cube::AnnulusView) = filter(r -> cube.rmin ≤ r ≤ cube.rmax, radii)
function _fix_range(radii, cube::MultiAnnulusView)
    rmin, rmax = extrema(cube.radii) .- cube.width / 2
    filter(r -> rmin ≤ r ≤ rmax, radii)
end

"""
    throughput(alg, cube, angles, psf, position;
               fwhm, snr=100, reduced_empty=nothing,
               verbose=true, kwargs...)

Calculate the throughput of `alg` by injecting `psf` into `cube` at the given `position` and measuring the relative photometry of the companion in the reduced frame. The photometry is measured using a circular aperture with a diameter matching the `fwhm`. Any additional `kwargs` will be passed to `alg` when it is called.

If `position` is a tuple or a vector, it will be parsed as the cartesian coordinate `(x, y)`. If `position` is a `CoordinateTransformations.Polar` it will be parsed as polar coordinates from the center of the cube. Note the `Polar` type expects angles in radians.

# Keyword Arguments
* `fwhm` - the full width at half-maximum
* `snr` - the target signal to noise ratio of the injected planets
* `reduced_empty` - the collapsed residual frame for estimating the noise. Will process using `alg` if not provided.
* `verbose` - show informative messages during process
"""
function throughput(alg, cube::AbstractArray{T,3}, angles, psf_model, position;
    fwhm, snr=100, reduced_empty=nothing, verbose=true, kwargs...) where T

    if reduced_empty === nothing
        verbose && @info "Calculating empty reduced frame"
        reduced_empty = alg(cube, angles; kwargs...)
    end

    cent = center(cube)[2:3]
    r, θ = _coordinates(position, reverse(cent))
    noise = annulus_noise(reduced_empty, fwhm, cent..., r, θ)
    A = snr * noise

    verbose && @info "Injecting companion at r=$r θ=$θ with A=$A"
    fake_comp = inject!(zero(reduced_empty), psf_model, Polar(r, deg2rad(θ)); A=A)
    fake_comp_cube = inject(cube, angles, psf_model, Polar(r, deg2rad(θ)); A=A)

    verbose && @info "Calculating reduced frame with fake companion injected"
    reduced = alg(fake_comp_cube, angles; kwargs...)

    x, y = convert(SVector, Polar(r, deg2rad(θ))) + cent
    ap = CircularAperture(x, y, fwhm/2)

    verbose && @info "Calculating throughput"
    injected_flux = photometry(ap, fake_comp).aperture_sum
    recovered_flux = photometry(ap, reduced .- reduced_empty).aperture_sum
    throughput = max(zero(T), recovered_flux / injected_flux)

    return throughput
end

# get r, θ from cartesian coords
_coordinates(pos::Tuple, center) = _coordinates(SVector(pos), center)
_coordinates(pos::AbstractVector, center) = _coordinates(convert(Polar, pos .- center), center)
_coordinates(pos::Polar, center) = (pos.r, rad2deg(pos.θ))

function annulus_noise(frame, fwhm, cy, cx, r, θ=0)
    δy, δx = sincosd(θ)
    x = r * δx + cx
    y = r * δy + cy
    Metrics.noise(frame, (x, y), fwhm)
end


text:
module OpenSpiel

include("$(@__DIR__)/../deps/deps.jl")

using CxxWrap
import CxxWrap:argument_overloads

@wrapmodule(LIB_OPEN_SPIEL)

include("patch.jl")

# export all
for n in names(@__MODULE__(); all=true)
    if Base.isidentifier(n) &&
        !startswith(String(n), "_") &&
        n ∉ (Symbol(@__MODULE__()), :eval, :include)
        @eval export $n
    end
end

function __init__()
    @initcxx
end

end

ext: jl


text:
# Test complete Golub-Kahan-Lanczos factorization
@testset "Complete Golub-Kahan-Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T, (n,n))
            v = A*rand(T, (n,)) # ensure v is in column space of A
            iter = GKLIterator(wrapop(A), wrapvec2(v), orth)
            verbosity = 1
            fact = @constinferred initialize(iter; verbosity = verbosity)
            while length(fact) < n
                @constinferred expand!(iter, fact; verbosity = verbosity)
                verbosity = 0
            end

            U = hcat(unwrapvec2.(basis(fact, :U))...)
            V = hcat(unwrapvec.(basis(fact, :V))...)
            B = rayleighquotient(fact)
            @test normres(fact) < 10*n*eps(real(T))
            @test U'*U ≈ I
            @test V'*V ≈ I
            @test A*V ≈ U*B
            @test A'*U ≈ V*B'

            @constinferred initialize!(iter, deepcopy(fact); verbosity = 1)
            states = collect(Iterators.take(iter, n)) # collect tests size and eltype?
            @test rayleighquotient(last(states)) ≈ B
        end
    end
end

# Test incomplete Golub-Kahan-Lanczos factorization
@testset "Incomplete Golub-Kahan-Lanczos factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64, Complex{Int})
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            if T == Complex{Int}
                A = rand(-100:100, (N, N)) + im * rand(-100:100, (N, N))
                v = rand(-100:100, (N,))
            else
                A = rand(T,(N,N))
                v = rand(T,(N,))
            end
            iter = @constinferred GKLIterator(wrapop(A), wrapvec2(v), orth)
            krylovdim = 3*n
            fact = @constinferred initialize(iter)
            while normres(fact) > eps(float(real(T))) && length(fact) < krylovdim
                @constinferred expand!(iter, fact)

                Ũ, Ṽ, B, r̃, β, e = fact
                U = hcat(unwrapvec2.(Ũ)...)
                V = hcat(unwrapvec.(Ṽ)...)
                r = unwrapvec2(r̃)
                @test U'*U ≈ I
                @test V'*V ≈ I
                @test norm(r) ≈ β
                @test A*V ≈ U*B + r*e'
                @test A'*U ≈ V*B'
            end

            fact = @constinferred shrink!(fact, div(n,2))
            U = hcat(unwrapvec2.(@constinferred basis(fact, :U))...)
            V = hcat(unwrapvec.(@constinferred basis(fact, :V))...)
            B = @constinferred rayleighquotient(fact)
            r = unwrapvec2(@constinferred residual(fact))
            β = @constinferred normres(fact)
            e = @constinferred rayleighextension(fact)
            @test U'*U ≈ I
            @test V'*V ≈ I
            @test norm(r) ≈ β
            @test A*V ≈ U*B + r*e'
            @test A'*U ≈ V*B'
        end
    end
end


text:
export merge_sort

function merge_sort(array::Vector)
    if length(array) > 1
        divided = divide(array)
        merged1 = merge_sort(divided[1])
        merged2 = merge_sort(divided[2])
        merged = merge(merged1, merged2)

        return merged
    end

    return array
end

function divide(array::Vector)
    split = trunc(Int64, length(array)/2)
    arr1 = array[1:split]
    arr2 = array[split+1:end]

    return (arr1, arr2)
end

function merge(arr1::Vector, arr2::Vector)
    merged = []

    while !(isempty(arr1) || isempty(arr2))
        if arr1[1] <= arr2[1]
            push!(merged, popfirst!(arr1))
        else
            push!(merged, popfirst!(arr2))
        end
    end

    if !isempty(arr1)
        append!(merged, arr1)
    end
    if !isempty(arr2)
        append!(merged, arr2)
    end

    return merged
end


text:
function draw_legend!(specification)
	@unpack P, group_pairs, style_pairs, group_dict, style_dict = specification

	leg, cb = draw_legend!(P, group_pairs, style_pairs, group_dict, style_dict)
end

function draw_legend!(P, groups, styles, group_dict, style_dict0)
	groups_ = collect(keys(group_dict))

	legends_grp = length(group_dict) > 0 ? legends_discrete(P, groups, group_dict) : nothing

	style_dict = deepcopy(style_dict0)
	# special case colorbar
	col_extr = pop!(style_dict, :color, nothing)
	if !isnothing(col_extr)
		cb = (limits = col_extr, title = var_lab(styles[:color]))
		#, width = 30, height = Relative(0.7)
	else
		cb = nothing
	end

	styles_ = collect(keys(style_dict))

	legends_stl = map(styles_) do k

		extr = style_dict[k]

		legend_continuous(P, k, extr, var_lab(styles[k]))
	end

	# Combine
	legends = legends_grp
	if isnothing(legends)
		legends = legends_stl
	else
		if !isnothing(legends_stl)
			legends = [legends; legends_stl]
		end
	end

	if length(legends) > 0
		legends = DataFrame(legends)
		leg = (legends.elements, legends.labels, legends.title)
	else
		leg = nothing
	end

	(; leg, cb)
end

function legends_discrete(P, group_pairs, group_dict)
	group_pairs = clean_groups(group_pairs)

	inv_pairs = invert_pairs(pairs(group_pairs))

	legends = map(inv_pairs) do (var, attrs)
		subdict = filter(kv -> kv[1] in attrs, group_dict)
		nts = invert_dict_of_pairs(subdict)

		legend_discrete(P, var, nts)
	end |> StructArray
end

function legend_discrete(P, var, level_kws_vec)
	legend = map(level_kws_vec) do (level, kws)
		element = legend_element(P; kws...)
		(; level, element)
	end |> StructArray

	(elements = legend.element, labels = legend.level, title = var_lab(var))
end

function legend_continuous(P, attribute, extrema, title, n_ticks=4)

	ticks = MakieLayout.locateticks(extrema..., n_ticks)
	elements = [legend_element(P; attribute => s) for s in ticks]
	labels = string.(ticks)

	(; elements, labels, title)
end

from_default_theme(attr) = AbstractPlotting.current_default_theme()[attr]

line_element(; color     = :black, #default_theme(:color),
			   linestyle = nothing,
			   linewidth = 1.5,
			   kwargs...) =
		LineElement(; color, linestyle, linewidth, kwargs...)

marker_element(; color       = from_default_theme(:color),
		        marker       = from_default_theme(:marker),
			    strokecolor  = :black,
			    markerpoints = [Point2f0(0.5, 0.5)],
		        kwargs...) =
		MarkerElement( ; color, marker, strokecolor, markerpoints, kwargs...)

poly_element(; color = from_default_theme(:color),
	 		   strokecolor = :transparent,
			   kwargs...) =
		PolyElement(; color, strokecolor, kwargs...)

legend_element(::Type{Scatter}; kwargs...) = marker_element(; kwargs...)

legend_element(::Type{Lines}; kwargs...) = line_element(; kwargs...)

legend_element(::Type{BarPlot}; linewidth = 0, strokecolor=:green, kwargs...) = poly_element(; linewidth, kwargs...)

function create_entrygroups(contents::AbstractArray,
    labels::AbstractArray{String},
    title::Optional{String} = nothing)

    if length(contents) != length(labels)
        error("Number of elements not equal: $(length(contents)) content elements and $(length(labels)) labels.")
    end

    entries = [LegendEntry(label, content) for (content, label) in zip(contents, labels)]
    entrygroups = Vector{EntryGroup}([(title, entries)])
end

function create_entrygroups(contentgroups::AbstractArray{<:AbstractArray},
    labelgroups::AbstractArray{<:AbstractArray},
    titles::AbstractArray{<:Optional{String}})

    if !(length(titles) == length(contentgroups) == length(labelgroups))
    error("Number of elements not equal: $(length(titles)) titles,     $(length(contentgroups)) content groups and $(length(labelgroups)) label     groups.")
    end

    entries = [[LegendEntry(l, pg) for (l, pg) in zip(labelgroup, contentgroup)]
        for (labelgroup, contentgroup) in zip(labelgroups, contentgroups)]

    entrygroups = Vector{EntryGroup}([(t, en) for (t, en) in zip(titles, entries)])
end

# Helpers
function invert_pairs(pairs)
	pairs = pairs |> collect

	attrs = first.(pairs)
	vars = last.(pairs)

	map(unique(vars)) do var
		inds = findall(==(var), vars)
		var => attrs[inds]
	end
end

function invert_dict_of_pairs(dict)
	attrs = collect(keys(dict))

	df = mapreduce(vcat, attrs) do attr
		map(dict[attr]) do (level, characteristics)
			(; level, attr, characteristics)
		end
	end |> DataFrame

	#@chain df begin
	combine(
		groupby(df, :level),
		[:attr, :characteristics] => vec_of_pairs => :kws
		) |>
	x -> NamedTuple.(eachrow(x))

end

vec_of_pairs(attrs, chars) = [[attr => char for (attr, char) in zip(attrs, chars)], ]

@testset "invert_pairs" begin
	nt = (color = :variable, linestyle = :variable)
	prs = pairs(nt)

	inverted = [:variable => [:color, :linestyle]]

	@test invert_pairs(prs) == inverted
end

@testset "invert_dict_of_pairs" begin
	attr_level_char = Dict(
		:color     => ["nominal" => :orange, "real" => :blue],
		:linestyle => ["nominal" => nothing, "real" => :dash]
		)

	level_attr_char = [
		(level = "nominal", kws = [:color     => :orange,
								   :linestyle => nothing])
		(level = "real",    kws = [:color     => :blue,
								   :linestyle => :dash])
		]

	@test invert_dict_of_pairs(attr_level_char) == level_attr_char
end

function clean_groups(group_pairs)
	for key in [:group, :stack, :dodge]
		if haskey(group_pairs, key)
			group_pairs = delete(group_pairs, key)
		end
	end
	group_pairs
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
##
## Visualization of the Mandelbrot set in Julia
##
## Tomáš Kalvoda, FIT CTU, 2015
##

using Images, Colors, FixedPointNumbers
using ProgressMeter

function iter_to_color(i, maxiter)
    clr = convert(RGB, HSV(360.*i/maxiter,1,1))
    return [clr.r, clr.g, clr.b]
end

function getcolor(z::Complex{Float64}, maxiter, radius)
    i = 1; u = 0 + 0im;
    while i <= maxiter
        u = u^2 + z
        if abs(u) > radius
            return iter_to_color(i,maxiter)
        end
        i += 1
    end
    return [1., 1., 1.]
end

function draw_image(;remin=-3.0, remax=3.0, immin=-3.0, immax=3.0, width=400, height=400, maxiter=30, radius=20.)
    redelta = (remax-remin)/(width-1)
    imdelta = (immax-immin)/(height-1)

    img = colorim(Array(Float64, (width,height,3)))
    img["spatialorder"] = ["x", "y"]

    p = Progress(width, 1)

    for j = 0:width-1
        for k = 0:height-1
            img[j+1, k+1, :] = getcolor(remin + j*redelta + 1im*(immin + k*imdelta), maxiter, radius)
        end
        next!(p)
    end

    return fliplr(img)
end

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
Random.seed!(0)

#sig = [1 2; 2 2]
#sig = [1 1; 2 3]
sig = [2 3]
#sig = [2 2]
#sig = [3 2; 2 3]

S = random_S0Graph(sig)
Sthin = complement(forget_S0(complement(S)))

w = random_bounded(S.n)

@time opt1 = dsw(Sthin, w, eps=solver_eps, verbose=1)[1]
@time opt2 = dsw(S, Ψ(S, w), eps=solver_eps)[1]
@test opt1 ≈ opt2*S.n  atol=tol

max_stars_repo_name: dstahlke/NoncommutativeGraphs

max_issues_repo_name: dstahlke/NoncommutativeGraphs

max_forks_repo_name: dstahlke/NoncommutativeGraphs


text:
using Statistics

"""
    linear_regression(X::Matrix, Y::Matrix, alpha::Float64, iterations::Int64)

    Linear Regression Algorithm using Gradient Descent
    to find the best fit line for the data.

    Arguments :
        - X : Training data set of shape [m x n]
        - Y : Labels corresponding to the training data set of shape [m x 1]
        - alpha : Learning rate
        - iterations : Number of iterations to train
    Returns :
        - W : Weight vector
        - b : Bias
        - J_history : Cost Function History

    gradient_descent(W::Matrix, b, X::Matrix, Y::Matrix, alpha::Float64, num_iters::Int64)

    Performs gradient descent to learn W and b

    Arguments
        - `W`: matrix of weights
        - `b`: bias vector
        - `X`: matrix of features
        - `Y`: vector of labels
        - `alpha`: learning rate ( generally ranges from 0.01 to 0.00001 )
        - `num_iters`: number of iterations to run gradient descent

    Returns:
        - `W`: matrix of weights
        - `b`: bias vector
        - `J_history`: vector of cost values after each iteration

    # Contributors:- [Navaneeth Sharma](https://github.com/Navaneeth-Sharma)
"""

function gradient_descent(
    W::Matrix,
    b,
    X::Matrix,
    Y::Matrix,
    alpha::Float64,
    num_iters::Int64,
)
    m = length(Y)
    J_history = []
    for i in 1:num_iters
        W = W - alpha .* (1 / m) * (X * (X'W .+ b - Y'))
        b = b - alpha .* (1 / m) * sum((X'W .+ b - Y'))
        push!(J_history, 1 / 2m * sum((X'W .+ b - Y') .^ 2))
    end
    return W, b, J_history
end

function linear_regression(
    X::Matrix,
    Y::Matrix,
    alpha::Float64,
    num_iters::Int64,
)
    W = zeros(1, size(X, 1))
    b = 0
    W, b, J_history = gradient_descent(W, b, X, Y, alpha, num_iters)
    return W, b, J_history
end

# Function to predict the output of a given data vector
predict(X::Matrix, W::Matrix, b::Float64) = X'W .+ b;

# Example usage
X = [1, 3, 2, 5, 7, 8, 8, 9, 10, 12]
Y = [2, 4, 3, 6, 8, 9, 9, 10, 11, 13]

X = reshape(X, 1, length(X));
Y = reshape(Y, 1, length(Y));

W, b, J_history = linear_regression(X, Y, 0.01, 1000);

output = predict(X, W, b)
println("The predicted is: $output")

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using Makie
using StaticArrays
using Printf

room_count = [0]  # global
v = Vector{Int64}()  # global

struct NODE
    dir_pair::String
    x::Union{Nothing,Float64}
    y::Union{Nothing,Float64}
    function NODE(; x=nothing, y=nothing, dir_pair::String)
        return new(dir_pair, x, y)
    end
end
is_x(n::NODE) = !isnothing(n.x)
is_y(n::NODE) = !isnothing(n.y)

###

mutable struct Attached_Type
    room
    dir::String
    Δ_pos::Float64
    door_pos_in_child::Float64
    function Attached_Type(room, dir::String, Δ_pos::Float64, door_pos_in_child::Float64)
        return new(room, dir, Δ_pos, door_pos_in_child)
    end
end

###

mutable struct Room
    id::Int64
    x::Float64
    y::Float64
    refine_method
    name::String
    attached::Vector{Attached_Type}
    nodes::Vector{NODE}
    scene_element
    is_active::Bool
    global_bounds::NTuple{4,Float64}
    is_pos_x::Bool
    is_pos_y::Bool
    function Room(x, y, refine_method, name::String)
        attached = Vector{Attached_Type}()
        nodes = Vector{NODE}()
        room_count[1] += 1
        return new(room_count[1], x, y, refine_method, name, attached, nodes, nothing, false)
    end
end

get_Δx(r::Room) = r.global_bounds[2] -  r.global_bounds[1]
get_Δy(r::Room) = r.global_bounds[4] -  r.global_bounds[3]
get_x0(r::Room) = r.global_bounds[1]
get_x1(r::Room) = r.global_bounds[2]
get_y0(r::Room) = r.global_bounds[3]
get_y1(r::Room) = r.global_bounds[4]

get_nx(r::Room) = get_nx(r.refine_method)
get_ny(r::Room) = get_ny(r.refine_method)

function search_for_x_constraint(r::Room)
    nodes = r.nodes
    for node_k = nodes
        is_x(node_k) && (return node_k.x)
    end
    return nothing
end
function search_for_y_constraint(r::Room)
    nodes = r.nodes
    for node_k = nodes
        is_y(node_k) && (return node_k.y)
    end
    return nothing
end

is_con_y(r::Room) = !isnothing(search_for_y_constraint(r))
is_con_x(r::Room) = !isnothing(search_for_x_constraint(r))

###

abstract type RefineType end

mutable struct RefineHoriz <: RefineType
    n::Int64
    function RefineHoriz()
        return new(1)
    end
end
mutable struct RefineVert <: RefineType
    n::Int64
    function RefineVert()
        return new(1)
    end
end
mutable struct RefineRectUR <: RefineType
    nx::Int64
    ny::Int64
    function RefineRectUR()
        return new(1, 1)
    end
end
mutable struct RefineRectDR <: RefineType
    nx::Int64
    ny::Int64
    function RefineRectDR()
        return new(1, 1)
    end
end
mutable struct RefineTri <: RefineType end

refine_x!(r::Union{RefineRectUR,RefineRectDR}) = (r.nx += 1; return nothing)
refine_y!(r::Union{RefineRectUR,RefineRectDR}) = (r.ny += 1; return nothing)
refine_x!(r::Union{RefineHoriz,RefineVert})    = (r.n += 1;  return nothing)
refine_y!(r::Union{RefineHoriz,RefineVert})    = (r.n += 1;  return nothing)

get_nx(r::Union{RefineRectUR,RefineRectDR}) = r.nx
get_ny(r::Union{RefineRectUR,RefineRectDR}) = r.ny

get_nx(r::RefineHoriz) = r.n
get_ny(r::RefineHoriz) = 1

get_nx(r::RefineVert) = 1
get_ny(r::RefineVert) = r.n

###

function add_attached!(r::Room, r_new::Room, dir::String, offset::Float64, opening::Float64)
    a = Attached_Type(r_new, dir, offset, opening)
    push!(r.attached, a)
end

plot_room!(r::Room) = plot_room!(r, 0.0, 0.0)

function plot_room!(r::Room, x0::Float64, y0::Float64)
    x1 = r.x + x0
    y1 = r.y + y0
    x = [x0, x1, x1, x0, x0]
    y = [y0, y0, y1, y1, y0]
    plot!(x, y, color="black")
    r.global_bounds = (x0, x1, y0, y1)
    for attached_k = r.attached
        x0_k, y0_k = calc_offset(r, attached_k, x0, y0)
        plot_room!(attached_k.room, x0_k, y0_k)
    end
end

function calc_offset(r::Room, a::Attached_Type, x0::Float64, y0::Float64)
    if a.dir == "r"
        x0 += r.x
        y0 += a.Δ_pos
        door_pos = y0 + a.door_pos_in_child
        push!(r.nodes,      NODE(y=door_pos, dir_pair="r"))
        push!(a.room.nodes, NODE(y=door_pos, dir_pair="l"))
    elseif a.dir == "u"
        y0 += r.y
        x0 += a.Δ_pos
        door_pos = x0 + a.door_pos_in_child
        push!(r.nodes,      NODE(x=door_pos, dir_pair="u"))
        push!(a.room.nodes, NODE(x=door_pos, dir_pair="d"))
    elseif a.dir == "d"
        y0 -= a.room.y
        x0 += a.Δ_pos
        door_pos = x0 + a.door_pos_in_child
        push!(r.nodes,      NODE(x=door_pos, dir_pair="d"))
        push!(a.room.nodes, NODE(x=door_pos, dir_pair="u"))
    elseif a.dir == "l"
        x0 -= a.room.x
        y0 += a.Δ_pos
        door_pos = y0 + a.door_pos_in_child
        push!(r.nodes,      NODE(y=door_pos, dir_pair="l"))
        push!(a.room.nodes, NODE(y=door_pos, dir_pair="r"))
    else
        error("not implemented error")
    end

    return x0, y0
end

function do_start_smart(v_con::Nothing, r::Room, n::Int64, is_x::Bool)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    v_space = (V1 - V0) / n
    v_start = V0 + v_space / 2
    return v_start, v_space, true
end

function find_con_side(r::Room, v_con::Float64, is_x::Bool)
    v_con = ifelse(is_x, search_for_x_constraint, search_for_y_constraint)(r)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    return abs(v_con - V0) > abs(v_con - V1)
end

function do_start_smart(v_con::Float64, r::Room, n::Int64, is_x::Bool)
    (n == 1) && (return v_con, 0.0, true)
    V0 = ifelse(is_x, get_x0, get_y0)(r)
    V1 = ifelse(is_x, get_x1, get_y1)(r)
    is_pos = abs(v_con - V0) < abs(v_con - V1)
    v0 = v_con
    v1 = V1 + (V0 - v0)
    v_space = (v1 - v0) / (n - 1)
    if !is_pos
        v_space *= -1.0
        v0 = v1
    end
    return v0, v_space, is_pos
end

function do_start_smart(r::Room, is_x::Bool)
    n = ifelse(is_x, get_nx, get_ny)(r)
    v_con = ifelse(is_x, search_for_x_constraint, search_for_y_constraint)(r)
    v_start, v_space, is_pos = do_start_smart(v_con, r, n, is_x)
    return v_start, v_space, n, is_pos
end

function plot_node(r::Room, ::RefineTri)
    n1 = r.nodes[1]
    n2 = r.nodes[2]

    is_y(n1) || error("not implemented")
    is_y(n2) || error("not implemented")
    (length(r.nodes) == 2) || error("not implemented")

    X0 = r.global_bounds[1]
    X1 = r.global_bounds[2]
    ΔX = X1 - X0

    x0 = X0 + ΔX * 0.25
    x1 = X0 + ΔX * 0.75

    y0 = n1.y
    y1 = n2.y
    y1_below = y1 - 2*(y1 - y0)

    x = [x0, x1, x1]
    y = [y0, y1, y1_below]

    return x, y, true, true
end

function plot_node(r::Room, ::Union{RefineRectDR,RefineVert,RefineHoriz,RefineRectUR})
    start_x, space_x, nx, is_pos_x = do_start_smart(r, true)
    start_y, space_y, ny, is_pos_y = do_start_smart(r, false)
    x = Vector{Float64}()
    y = Vector{Float64}()
    for ix = 0:(nx-1)
        for iy = 0:(ny-1)
            push!(x, start_x + space_x * ix)
            push!(y, start_y + space_y * iy)
        end
    end
    return x, y, is_pos_x, is_pos_y
end

function plot_node!(r::Room)
    x, y, is_pos_x, is_pos_y = plot_node(r, r.refine_method)
    r.is_pos_x = is_pos_x
    r.is_pos_y = is_pos_y

    # prevent remove points from previous iterations if necessary
    index_in_v = findall(v .== r.id)
    if length(index_in_v) != 0
        index = only(index_in_v)
        delete!(scene, scene[end - index + 1])
        deleteat!(v, index)
    end
    scatter!(x, y)
    pushfirst!(v, r.id)

    for attached_k = r.attached # recurse
        is_pos_x, is_pos_y = plot_node!(attached_k.room)
    end

    return is_pos_x, is_pos_y
end

###

struct Graph
    points::Vector{SVector{2,Float64}}
    edges::Vector{SVector{2,Int64}}
    function Graph()
        points = Vector{SVector{2,Float64}}()
        edges = Vector{SVector{2,Int64}}()
        return new(points, edges)
    end
end

function add_graph_room!(g::Graph, r::Room, ::RefineTri, ee_parent)
    x, y = plot_node(r, RefineTri())

    len_orig = length(g.points)

    for k = 1:3
        push!(g.points, SVector{2,Float64}(x[k], y[k]))
    end

    for k = 1:3
        push!(g.edges, SVector{2,Int64}(k, rem(k, 3) + 1) .+ len_orig)
    end

    return ExtremeInd(1, 1, 3, 2, len_orig)
end

function add_graph_room!(g::Graph, r::Room, ::Union{RefineRectDR,RefineVert,RefineHoriz,RefineRectUR}, ee_parent)
    ind_2_node(ind_y::Int64, ind_x::Int64) = ind_y + (ind_x - 1) * ny

    start_x, space_x, nx, is_reg_x = do_start_smart(r, true)
    start_y, space_y, ny, is_reg_y = do_start_smart(r, false)
    point_matrix = Matrix{SVector{2,Float64}}(undef, ny, nx)
    edge_vec = Vector{SVector{2,Int64}}()
    for ix = 1:nx
        for iy = 1:ny
            x = start_x + space_x * (ix - 1)
            y = start_y + space_y * (iy - 1)
            point_matrix[iy, ix] = SVector{2,Float64}(x, y)
            index = ind_2_node(iy, ix)
            if ix != nx
                index_right = ind_2_node(iy, ix + 1)
                push!(edge_vec, SVector{2,Int64}(index, index_right))
            end
            if iy != ny
                index_down = ind_2_node(iy + 1, ix)
                push!(edge_vec, SVector{2,Int64}(index, index_down))
            end
        end
    end

    len_orig = length(g.points)

    for point_k = point_matrix
        # point_k = SVector{2,Float64}(point_k[1] + randn() * 0.15, point_k[2] + randn() * 0.15)
        push!(g.points, point_k)
    end

    for edge_k = edge_vec
        push!(g.edges, edge_k .+ len_orig)
    end

    return ExtremeInd(ind_2_node(ny, 1), ind_2_node(1, 1), ind_2_node(ny, nx), ind_2_node(1, nx), len_orig)
end

struct ExtremeInd
    i_ul::Int64
    i_dl::Int64
    i_ur::Int64
    i_dr::Int64
    function ExtremeInd(i_ul, i_dl, i_ur, i_dr, len_orig)
        return new(i_ul + len_orig, i_dl + len_orig, i_ur + len_orig, i_dr + len_orig)
    end
end

function add_graph!(g::Graph, r::Room, a=nothing)
    ee_parent = add_graph_room!(g, r, r.refine_method, a)

    for attached_k = r.attached
        ee_child = add_graph!(g, attached_k.room, attached_k)

        if attached_k.dir == "r"
            if r.is_pos_y
                i_parent = ee_parent.i_dr
            else
                i_parent = ee_parent.i_ur
            end
            i_child = ee_child.i_ul
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "u"
            v_con = search_for_x_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, true), ee_child.i_dr, ee_child.i_dl)
            i_parent = ifelse(find_con_side(r, v_con, true), ee_parent.i_ur, ee_parent.i_ul)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "d"
            v_con = search_for_x_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, true), ee_child.i_ur, ee_child.i_ul)
            i_parent = ifelse(find_con_side(r, v_con, true), ee_parent.i_dr, ee_parent.i_dl)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        elseif attached_k.dir == "l"
            v_con = search_for_y_constraint(attached_k.room)
            i_child = ifelse(find_con_side(attached_k.room, v_con, false), ee_child.i_ur, ee_child.i_dr)
            i_parent = ifelse(find_con_side(r, v_con, false), ee_parent.i_ul, ee_parent.i_dl)
            push!(g.edges, SVector{2,Int64}(i_parent, i_child))
        end

    end

    return ee_parent
end

function print_graph(g::Graph)
    println("===================")
    println()

    for point_k = g.points
        s1 = @sprintf("%0.5f", point_k[1])
        s2 = @sprintf("%0.5f", point_k[2])
        println(s1, ", ", s2)
    end
    println()
    for edge_k = g.edges
        println(edge_k[1] - 1, ", ", edge_k[2] - 1)
    end
    println()
end

function plot_graph(g::Graph)

    # # prevent remove points from previous iterations if necessary
    # index_in_v = findall(v .== -9999)
    # if length(index_in_v) != 0
    #     index = only(index_in_v)
    #     delete!(scene, scene[end - index + 1])
    #     deleteat!(v, index)
    # end

    for edge_k = g.edges
        p1 = g.points[edge_k[1]]
        p2 = g.points[edge_k[2]]
        plot!([p1[1], p2[1]], [p1[2], p2[2]])
    end

    # pushfirst!(v, -9999)
end

###

function create_layout(is_debug::Bool=false)
    width = 1.0

    r1 = Room(2.0, 6.0, RefineRectUR(), "room_1")

    is_debug && (r1.refine_method.nx = 3)
    is_debug && (r1.refine_method.ny = 5)

    r2 = Room(2.0, 2.0, RefineTri(), "room_2")
    add_attached!(r1, r2, "r", 4.5, 1.0)

    r3 = Room(2.0, width, RefineHoriz(), "room_3")
    add_attached!(r2, r3, "r", 0.0, width / 2)

    is_debug && (r3.refine_method.n = 4)

    r_last_hall = r3

    for k = 1:4
        r4 = Room(2.0, width, RefineHoriz(), "room_4")
        add_attached!(r_last_hall, r4, "r", 0.0, width / 2)

        r5 = Room(1.7, 1.7, RefineRectDR(), "room_5")
        add_attached!(r4, r5, "u", 0.1, 1.5)

        r_last_hall = r4

        is_debug && (r_last_hall.refine_method.n = 2)
        is_debug && (r5.refine_method.nx = 2)
        is_debug && (r5.refine_method.ny = 2)
    end

    for k = 1:3
        r6 = Room(width, 2.4, RefineVert(), "room_6")
        Δx = ifelse(k == 1, 1.0, 0.0)

        add_attached!(r_last_hall, r6, "d", Δx, width/2)

        is_debug && (r6.refine_method.n = 2)

        r7 = Room(1.3, width, RefineHoriz(), "room_7")
        add_attached!(r6, r7, "l", 0.3, width / 2)

        is_debug && (r7.refine_method.n = 3)

        r_last_hall = r6
    end

    w8 = 3.0
    r8 = Room(w8, 5.0, RefineRectUR(), "room_8")
    add_attached!(r_last_hall, r8, "d", - w8 + width, w8 - width / 2)

    is_debug && (r8.refine_method.nx = 3)
    is_debug && (r8.refine_method.ny = 5)

    return r1
end

scene = Scene()
r1 = create_layout()
plot_room!(r1)
plot_node!(r1)

### Event-based stuff

function print_active_room(r::Room)
    if r.is_active
        println(r.name, " is active!")
    end

    for attached_k = r.attached
        print_active_room(attached_k.room)
    end
end

function set_all_room_inactive(r::Room)
    r.is_active = false
    for attached_k = r.attached
        set_all_room_inactive(attached_k.room)
    end
end

function set_active_room(r::Room, x::Float32, y::Float32)
    if (r.global_bounds[1] < x < r.global_bounds[2])
        if (r.global_bounds[3] < y < r.global_bounds[4])
            r.is_active = true
        end
    end
    for attached_k = r.attached
        set_active_room(attached_k.room, x, y)
    end
end

function refine_active!(r::Room, fun_name::Function)
    r.is_active && fun_name(r.refine_method)
    for attached_k = r.attached
        refine_active!(attached_k.room, fun_name)
    end
end

display(scene)

xlims!(-1, 14.0)
ylims!(-6, 8.0)

# g = Graph()
# add_graph!(g, r1)
# plot_graph(g)
#
#
# print_graph(g)

if true
    dir = lift(scene.events.keyboardbuttons) do but
        is_lr = ispressed(but, Keyboard.left) || ispressed(but, Keyboard.right)
        is_ud = ispressed(but, Keyboard.up) || ispressed(but, Keyboard.down)

        if is_lr
            println("in keyboard event loop x")
            refine_active!(r1, refine_x!)
            plot_node!(r1)
        end

        if is_ud
            println("in keyboard event loop y")
            refine_active!(r1, refine_y!)
            plot_node!(r1)
        end

        if ispressed(but, Keyboard.space)
            g = Graph()
            add_graph!(g, r1)
            plot_graph(g)
            print_graph(g)
        end
        return 0
    end

    pos = lift(scene.events.mouseposition) do mpos
        # println("HIT")

        x0, x1 = mouseposition(scene)

        set_all_room_inactive(r1)
        set_active_room(r1, x0, x1)
        return 0
    end
end

println()
println()
println("DONE!!!")

# plot!([0.0, 1.0], [0.0, 3.0], linewidth=5)
# pushfirst!(v, 1)
# plot!([0.0, 1.0], [0.0, 3.0] .+ 0.5, linewidth=4)
# pushfirst!(v, 2)
# plot!([0.0, 1.0], [0.0, 3.0] .+ 1.0, linewidth=8)
# pushfirst!(v, 3)
#
# i2 = only(findall(v .== 2))
# delete!(scene, scene[end-i2+1])
# deleteat!(v, i2)
# plot!([0.0, 1.0], [0.0, 3.0] .+ 0.5, color="blue", linewidth=1)
# pushfirst!(v, 2)
#
#
# i2 = only(findall(v .== 3))
# delete!(scene, scene[end-i2+1])
# deleteat!(v, i2)
# plot!([0.0, 1.0], [0.0, 3.0] .+ 1.0, color="red", linewidth=1)
# pushfirst!(v, 3)

# delete!(scene, scene[end])

# plot_room!(r2)
# plot_room!(r3)
#
# r = [r1, r2, r3]

# delete!(scene, scene[end])

# r_in = brute_search_room(r1, x0, x1)
#
# if !isnothing(r_in)
#     println("in_room ", r_in.id)
# end

# println("lift(scene.events.keyboardbuttons): ", lift(scene.events.keyboardbuttons))

# println("scene.events.keyboardbuttons: ", scene.events.keyboardbuttons)

# ispressed(but, Keyboard.space) && println("space")  # return 0

# plot_room!(r1)

# sleep(0.5)

# do stuff
# b1 = button("Test Button")
# lift(b1[end].clicks) do clicks
#     println("Button was clicked!")
#     #Your function goes here
# end

#


text:
include("DGFVModel_kernels.jl")
struct DGFVModel{BL, G, FVR, NFND, NFD, GNF, AS, DS, HDS, D, DD, MD, GF, TF} <:
       SpaceDiscretization
    balance_law::BL
    grid::G
    fv_reconstruction::FVR
    numerical_flux_first_order::NFND
    numerical_flux_second_order::NFD
    numerical_flux_gradient::GNF
    state_auxiliary::AS
    state_gradient_flux::DS
    states_higher_order::HDS
    direction::D
    diffusion_direction::DD
    modeldata::MD
    gradient_filter::GF
    tendency_filter::TF
    check_for_crashes::Bool
end

function DGFVModel(
    balance_law,
    grid,
    fv_reconstruction,
    numerical_flux_first_order,
    numerical_flux_second_order,
    numerical_flux_gradient;
    fill_nan = false,
    check_for_crashes = false,
    state_auxiliary = create_state(
        balance_law,
        grid,
        Auxiliary(),
        fill_nan = fill_nan,
    ),
    state_gradient_flux = create_state(balance_law, grid, GradientFlux()),
    states_higher_order = (
        create_state(balance_law, grid, GradientLaplacian()),
        create_state(balance_law, grid, Hyperdiffusive()),
    ),
    direction = EveryDirection(),
    diffusion_direction = direction,
    modeldata = nothing,
    gradient_filter = nothing,
    tendency_filter = nothing,
)
    # Make sure we are FVM in the vertical
    @assert polynomialorders(grid)[end] == 0
    @assert isstacked(grid.topology)
    state_auxiliary =
        init_state(state_auxiliary, balance_law, grid, direction, Auxiliary())
    DGFVModel(
        balance_law,
        grid,
        fv_reconstruction,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        state_auxiliary,
        state_gradient_flux,
        states_higher_order,
        direction,
        diffusion_direction,
        modeldata,
        gradient_filter,
        tendency_filter,
        check_for_crashes,
    )
end

"""
    (dgfvm::DGFVModel)(tendency, state_prognostic, _, t, α, β)

Uses spectral element discontinuous Galerkin in the horizontal and finite volume
in the vertical to compute the tendency.

Computes the tendency terms compatible with `IncrementODEProblem`

    tendency .= α .* dQdt(state_prognostic, p, t) .+ β .* tendency

The 4-argument form will just compute

    tendency .= dQdt(state_prognostic, p, t)
"""
function (dgfvm::DGFVModel)(tendency, state_prognostic, _, t, α, β)
    device = array_device(state_prognostic)

    FT = eltype(state_prognostic)
    num_state_prognostic = number_states(dgfvm.balance_law, Prognostic())
    num_state_gradient_flux = number_states(dgfvm.balance_law, GradientFlux())
    @assert 0 == number_states(dgfvm.balance_law, Hyperdiffusive())
    num_state_tendency = size(tendency, 2)

    if num_state_prognostic < num_state_tendency && β != 1
        # If we don't operate on the full state, then we need to scale here instead of volume_tendency!
        tendency .*= β
        β = β != 0 # if β==0 then we can avoid the memory load in volume_tendency!
    end

    communicate =
        !(
            isstacked(dgfvm.grid.topology) &&
            typeof(dgfvm.direction) <: VerticalDirection
        )

    update_auxiliary_state!(
        dgfvm,
        dgfvm.balance_law,
        state_prognostic,
        t,
        dgfvm.grid.topology.realelems,
    )

    exchange_state_gradient_flux = NoneEvent()
    exchange_state_prognostic = NoneEvent()

    comp_stream = Event(device)

    if communicate
        exchange_state_prognostic = MPIStateArrays.begin_ghost_exchange!(
            state_prognostic;
            dependencies = comp_stream,
        )
    end

    if num_state_gradient_flux > 0
        ########################
        # Gradient Computation #
        ########################

        comp_stream = launch_volume_gradients!(
            dgfvm,
            state_prognostic,
            t;
            dependencies = comp_stream,
        )

        comp_stream = launch_interface_gradients!(
            dgfvm,
            state_prognostic,
            t;
            surface = :interior,
            dependencies = comp_stream,
        )

        if communicate
            exchange_state_prognostic = MPIStateArrays.end_ghost_exchange!(
                state_prognostic;
                dependencies = exchange_state_prognostic,
                check_for_crashes = dgfvm.check_for_crashes,
            )

            # Update_aux may start asynchronous work on the compute device and
            # we synchronize those here through a device event.
            checked_wait(
                device,
                exchange_state_prognostic,
                nothing,
                dgfvm.check_for_crashes,
            )
            update_auxiliary_state!(
                dgfvm,
                dgfvm.balance_law,
                state_prognostic,
                t,
                dgfvm.grid.topology.ghostelems,
            )
            exchange_state_prognostic = Event(device)
        end

        comp_stream = launch_interface_gradients!(
            dgfvm,
            state_prognostic,
            t;
            surface = :exterior,
            dependencies = (comp_stream, exchange_state_prognostic),
        )

        if dgfvm.gradient_filter !== nothing
            comp_stream = Filters.apply_async!(
                dgfvm.state_gradient_flux,
                1:num_state_gradient_flux,
                dgfvm.grid,
                dgfvm.gradient_filter;
                dependencies = comp_stream,
            )
        end

        if communicate
            if num_state_gradient_flux > 0
                exchange_state_gradient_flux =
                    MPIStateArrays.begin_ghost_exchange!(
                        dgfvm.state_gradient_flux,
                        dependencies = comp_stream,
                    )
            end
        end

        if num_state_gradient_flux > 0
            # Update_aux_diffusive may start asynchronous work on the compute device
            # and we synchronize those here through a device event.
            checked_wait(device, comp_stream, nothing, dgfvm.check_for_crashes)
            update_auxiliary_state_gradient!(
                dgfvm,
                dgfvm.balance_law,
                state_prognostic,
                t,
                dgfvm.grid.topology.realelems,
            )
            comp_stream = Event(device)
        end
    end

    ###################
    # RHS Computation #
    ###################
    comp_stream = launch_volume_tendency!(
        dgfvm,
        tendency,
        state_prognostic,
        t,
        α,
        β;
        dependencies = comp_stream,
    )

    comp_stream = launch_interface_tendency!(
        dgfvm,
        tendency,
        state_prognostic,
        t,
        α,
        β;
        surface = :interior,
        dependencies = comp_stream,
    )

    if communicate
        if num_state_gradient_flux > 0
            exchange_state_gradient_flux = MPIStateArrays.end_ghost_exchange!(
                dgfvm.state_gradient_flux;
                dependencies = exchange_state_gradient_flux,
                check_for_crashes = dgfvm.check_for_crashes,
            )

            # Update_aux_diffusive may start asynchronous work on the
            # compute device and we synchronize those here through a device
            # event.
            checked_wait(
                device,
                exchange_state_gradient_flux,
                nothing,
                dgfvm.check_for_crashes,
            )
            update_auxiliary_state_gradient!(
                dgfvm,
                dgfvm.balance_law,
                state_prognostic,
                t,
                dgfvm.grid.topology.ghostelems,
            )
            exchange_state_gradient_flux = Event(device)
        else
            exchange_state_prognostic = MPIStateArrays.end_ghost_exchange!(
                state_prognostic;
                dependencies = exchange_state_prognostic,
                check_for_crashes = dgfvm.check_for_crashes,
            )

            # Update_aux may start asynchronous work on the compute device and
            # we synchronize those here through a device event.
            checked_wait(
                device,
                exchange_state_prognostic,
                nothing,
                dgfvm.check_for_crashes,
            )
            update_auxiliary_state!(
                dgfvm,
                dgfvm.balance_law,
                state_prognostic,
                t,
                dgfvm.grid.topology.ghostelems,
            )
            exchange_state_prognostic = Event(device)
        end
    end

    comp_stream = launch_interface_tendency!(
        dgfvm,
        tendency,
        state_prognostic,
        t,
        α,
        β;
        surface = :exterior,
        dependencies = (
            comp_stream,
            exchange_state_prognostic,
            exchange_state_gradient_flux,
            # XXX: This is disabled until FVM with hyperdiffusion for DG is implemented: exchange_Qhypervisc_grad,
        ),
    )

    if dgfvm.tendency_filter !== nothing
        comp_stream = Filters.apply_async!(
            tendency,
            1:num_state_tendency,
            dgfvm.grid,
            dgfvm.tendency_filter;
            dependencies = comp_stream,
        )
    end

    # The synchronization here through a device event prevents CuArray based and
    # other default stream kernels from launching before the work scheduled in
    # this function is finished.
    checked_wait(device, comp_stream, nothing, dgfvm.check_for_crashes)
end

function fvm_balance!(
    balance_func,
    m::BalanceLaw,
    state_auxiliary::MPIStateArray,
    grid,
)
    device = array_device(state_auxiliary)

    dim = dimensionality(grid)
    N = polynomialorders(grid)
    Nq = N .+ 1
    Nqj = dim == 2 ? 1 : Nq[2]

    topology = grid.topology
    elems = topology.elems
    nelem = length(elems)
    nvertelem = topology.stacksize
    horzelems = fld1(first(elems), nvertelem):fld1(last(elems), nvertelem)

    event = Event(device)
    event = kernel_fvm_balance!(device, (Nq[1], Nqj))(
        balance_func,
        m,
        Val(nvertelem),
        state_auxiliary.data,
        grid.vgeo,
        horzelems;
        ndrange = (length(horzelems) * Nq[1], Nqj),
        dependencies = (event,),
    )
    wait(device, event)
end


text:
module LogSystem

    using Formatting
    using MPI

    export writeLog

    function writeLog(args...; force :: Bool = false)

        comm = MPI.COMM_WORLD
        rank = MPI.Comm_rank(comm)

        if force || rank == 0
            println(format(args...))
        end
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Plots, CSV
using GPARatScale, Stheno

using Suppressor  # for supressing plot warnings

PLOTTING_SAMPLES = 1

train_data = CSV.read(joinpath(@__DIR__, "datasets\\eeg\\eeg_train.csv"))
test_data = CSV.read(joinpath(@__DIR__, "datasets\\eeg\\eeg_test.csv"))

# Split training data into columns
train_time = train_data[:, 1]
train_f3 = train_data[:, 2]
train_f4 = train_data[:, 3]
train_f5 = train_data[:, 4]
train_f6 = train_data[:, 5]
train_fz = train_data[:, 6]
train_f1 = train_data[:, 7]
train_f2 = train_data[:, 8]
# Split test data into columns
test_time = test_data[:, 1]
test_fz = test_data[:, 2]
test_f1 = test_data[:, 3]
test_f2 = test_data[:, 4]

# Range for which we have data for fz, f1, and f2
data_range = 1:156

# Train IGPs
fz_gp_post = create_optim_gp_post(
    train_time[data_range],
    [train_fz[data_range]...],
    kernel_structure = Matern52();
    # i_log_l=-2.0, i_log_process_var=2.0, i_log_noise_sigma=-2.0,
    debug = true,
)
f1_gp_post = create_optim_gp_post(
    train_time[data_range],
    [train_f1[data_range]...],
    kernel_structure = Matern52();
    # i_log_l=-2.0, i_log_process_var=2.0, i_log_noise_sigma=-2.0,
    debug = true,
)
f2_gp_post = create_optim_gp_post(
    train_time[data_range],
    [train_f2[data_range]...],
    kernel_structure = Matern12();
    i_log_l=-2.0, i_log_process_var=1.0, i_log_noise_sigma=-3.0,
    debug = true,
)

# Train GPAR models
fz_gpar_post = create_optim_gpar_post(
    [
        train_time[data_range], train_f3[data_range], train_f4[data_range],
        train_f5[data_range], train_f6[data_range]
    ],
    [train_fz[data_range]...];
    time_kernel = Matern52(),
    out_kernel = Matern52(),
    i_log_time_l=-3.0, i_log_time_var=1.0, i_log_out_l=6.0,
    i_log_out_var=4.0, i_log_noise_sigma=-2.0,
    multi_input=true,
    debug = true,
)
f1_gpar_post = create_optim_gpar_post(
    [
        train_time[data_range], train_f3[data_range], train_f4[data_range],
        train_f5[data_range], train_f6[data_range], train_fz[data_range]
    ],
    [train_f1[data_range]...];
    time_kernel = Matern52(),
    out_kernel = Matern52(),
    # i_log_time_l=-3.0, i_log_time_var=1.0, i_log_out_l=6.0,
    # i_log_out_var=4.0, i_log_noise_sigma=-2.0,
    multi_input=true,
    debug = true,
)
f2_gpar_post = create_optim_gpar_post(
    [
        train_time[data_range], train_f3[data_range], train_f4[data_range],
        train_f5[data_range], train_f6[data_range], train_fz[data_range],
        train_f1[data_range]
    ],
    [train_f2[data_range]...];
    time_kernel = Matern52(),
    out_kernel = Matern52(),
    # i_log_time_l=-3.0, i_log_time_var=1.0, i_log_out_l=6.0,
    # i_log_out_var=4.0, i_log_noise_sigma=-2.0,
    multi_input=true,
    debug = true,
)

# @suppress begin
# PLOTTING
gr();
overall_plot = plot(layout = (3, 1), legend = false);
# DATA PLOTS
function plot_data(plot_ref, train_y, test_y; xlimit=[0.35, 1])
    scatter!(
        plot_ref,
        train_time[data_range],
        train_y[data_range],
        color = :orange,
        markersize = 2.3,
        markeralpha = 0.8,
        markershape = :circle,
        xlims = xlimit,
        label = "FZ train",
    )
    scatter!(
        plot_ref,
        test_time,
        test_y,
        color = :orange,
        markersize = 1.8,
        markeralpha = 0.8,
        markershape = :square,
        label = "FZ test",
    )
end
# Plot fz
plot_data(overall_plot[1], train_fz, test_fz)
# Plot f1
plot_data(overall_plot[2], train_f1, test_f1)
# Plot f2
plot_data(overall_plot[3], train_f2, test_f2)

# PLOT IGP
plot!(
    overall_plot[1],
    fz_gp_post(train_time),
    samples = PLOTTING_SAMPLES,
    color = :green,
    fillalpha = 0.5,
    linealpha = 1,
)
plot!(
    overall_plot[2],
    f1_gp_post(train_time),
    samples = PLOTTING_SAMPLES,
    color = :green,
    fillalpha = 0.3,
    linealpha = 1,
)
plot!(
    overall_plot[3],
    f2_gp_post(train_time),
    samples = PLOTTING_SAMPLES,
    color = :green,
    fillalpha = 0.3,
    linealpha = 1,
)

# PLOT GPAR
function plot_gpar(plot_ref, means, stds; ylimits=nothing, standard_devs=3)
    # Helper functions for plotting mean of GPAR and standard deviation
    # Plot mean
    plot!(plot_ref, train_time, means, color = :blue, linealpha = 1)
    plot!(plot_ref, train_time, [means means];
    linewidth=0.0,
    linecolor=:blue,
    fillrange=[means .- standard_devs .* stds, means .+ standard_devs * stds],
    fillalpha=0.3,
    fillcolor=:blue,
    );
end
# FZ
fz_input = to_ColVecs(
    [
        train_time, train_f3,  train_f4, train_f5, train_f6
    ])
fz_marginals = marginals(fz_gpar_post(fz_input))

fz_means = [f.μ for f in fz_marginals]
fz_stds = [f.σ for f in fz_marginals]
plot_gpar(overall_plot[1], fz_means, fz_stds)
# F1
f1_input = to_ColVecs(
    [
        train_time, train_f3,  train_f4, train_f5, train_f6, fz_means
    ])
f1_marginals = marginals(f1_gpar_post(f1_input))

f1_means = [f.μ for f in f1_marginals]
f1_stds = [f.σ for f in f1_marginals]
plot_gpar(overall_plot[2], f1_means, f1_stds)
# F2
f2_input = to_ColVecs(
    [
        train_time, train_f3,  train_f4, train_f5, train_f6, fz_means, f1_means
    ])
f2_marginals = marginals(f2_gpar_post(f2_input))

f2_means = [f.μ for f in f2_marginals]
f2_stds = [f.σ for f in f2_marginals]
plot_gpar(overall_plot[3], f2_means, f2_stds)

display(overall_plot)
# end

max_stars_repo_name: willtebbutt/GPAR-at-scale

max_issues_repo_name: willtebbutt/GPAR-at-scale

max_forks_repo_name: willtebbutt/GPAR-at-scale


text:
using DelimitedFiles

input = joinpath(@__DIR__, "input")
raw = readdlm(input, String)

function decimalize(vec)
    reverse!(vec)
    return sum(2^i * x for (i, x) in zip(0:length(vec) - 1, vec))
end

function rating(diagnostics, mcf)
    for i in 1:length(first(diagnostics))
        mc = mcf(diagnostics)
        filter!(d -> d[i] == mc[i], diagnostics)
        if length(diagnostics) == 1
            return first(diagnostics)
        end
    end
end

diagnostics = map(vec(raw)) do diagnostic
    map(n -> parse(Int, n), collect(diagnostic))
end

sums = sum(diagnostics)

gamma = decimalize(sums .>= (length(diagnostics) / 2))
epsilon = decimalize(sums .< (length(diagnostics) / 2))

oxygen = decimalize(rating(copy(diagnostics), ds -> sum(ds) .>= length(ds) / 2))
scrubber = decimalize(rating(copy(diagnostics), ds -> sum(ds) .< length(ds) / 2))

p1 = gamma * epsilon
p2 = oxygen * scrubber

println("-----------------------------------------------------------------------")
println("binary diagnostic -- part one :: $p1")
println("binary diagnostic -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 738234)
@assert(p2 == 3969126)


text:
# TODO: Document this.
macro limit(src, n)
    Expr(
        :call,
        :Limit,
        esc(src),
        n,
    )
end


text:
using LinearAlgebra
using SkyCoords

ra = 89.014303
dec = 13.924912
ω = 37.59

α = deg2rad(ra)
δ = deg2rad(dec)
d = 1000.0 / ω
println("α=",α,"\tδ=",δ,"\td=",d)

pmra = 372.72
pmdec = -483.69
radial_velocity = 0.37

vSUN = [11.1, 232.24, 7.25]#km/s

tropical_year = 31556925.445#seconds in one tropical year
parsec = 30856775812800#km in one parsec
factor = parsec / tropical_year
println("parsec,year conversion factor = ", factor)

#SI_unit = 2.44507E-17
#println("SI mas/year conversion unit:", SI_unit)

#=
μRA = pmra / cos(δ)
μDEC = pmdec

vRA = 4.74 * d * μRA / 1000.0#km/s
vDEC = 4.74 * d * μDEC / 1000.0#km/s
vR = radial_velocity#km/s

println("vRA:", vRA, "\tvREC:", vDEC, "vR:", vR)
println("adj. vRA:", vRA+vSUN[1], "\tadj. vREC:", vDEC+vSUN[2], "adj. vR:", vR-vSUN[3])
=#

η = deg2rad(58.5986320306)
αGC = deg2rad(266.4051)
δGC = deg2rad(-28.936175)
dGC = 8300
xGC = [1, 0, 0]
z = 27
θ = asin(z / dGC)
println("θ=",θ)

#positions
r_icrs = d * [cos(α) * cos(δ), sin(α) * cos(δ), sin(δ)]

#velocities
#pmra = 4.74 * d * pmra / 1000.0#in line with km/s
#pmdec = 4.74 * d * pmdec / 1000.0#in line with km/s
#v_icrs = d * [cos(α) * cos(δ), sin(α) * cos(δ), sin(δ)]
pmra = pmra / tropical_year / 1000
pmdec = pmdec / tropical_year / 1000
d *= parsec

V_icrs = [radial_velocity*cos(α)*cos(δ)-d*sin(α)*pmra-d*cos(α)*sin(δ)*pmdec, radial_velocity*sin(α)*cos(δ)+d*cos(α)*pmra-d*sin(α)*sin(δ)*pmdec, radial_velocity*sin(δ)+d*cos(δ)*pmdec]
println("V_icrs=",V_icrs)

R₁ = Matrix(undef, 3, 3)
R₁[1,1] = cos(δGC)
R₁[1,2] = 0
R₁[1,3] = sin(δGC)
R₁[2,:] = [0 1 0]
R₁[3,1] = - sin(δGC)
R₁[3,2] = 0
R₁[3,3] = cos(δGC)
println("R₁:", R₁)

R₂ = Matrix(undef, 3, 3)
R₂[1,1] = cos(αGC)
R₂[1,2] = sin(αGC)
R₂[1,3] = 0
R₂[2,1] = - sin(αGC)
R₂[2,2] = cos(αGC)
R₂[2,3] = 0
R₂[3,:] = [0 0 1]
println("R₂:", R₂)

R₃ = Matrix(undef, 3, 3)
R₃[1,:] = [1 0 0]
R₃[2,1] = 0
R₃[2,2] = cos(η)
R₃[2,3] = sin(η)
R₃[3,1] = 0
R₃[3,2] = - sin(η)
R₃[3,3] = cos(η)
println("R₃:", R₃)

R = R₃ * R₁ * R₂
println("R:", R)

R = Matrix(undef, 3, 3)
R[1,1] = cos(αGC) * cos(δGC)
R[1,2] = cos(δGC) * sin(αGC)
R[1,3] = sin(δGC)
R[2,1] = - cos(αGC) * sin(δGC) * sin(η) - sin(αGC) * cos(η)
R[2,2] = - sin(αGC) * sin(δGC) * sin(η) + cos(αGC) * cos(η)
R[2,3] = cos(δGC) * sin(η)
R[3,1] = - cos(αGC) * sin(δGC) * cos(η) + sin(αGC) * sin(η)
R[3,2] = - sin(αGC) * sin(δGC) * cos(η) - cos(αGC) * sin(η)
R[3,3] = cos(δGC) * cos(η)
println("R:", R)

H = Matrix(undef, 3, 3)
H[1,1] = cos(θ)
H[1,2] = 0
H[1,3] = sin(θ)
H[2,1] = 0
H[2,2] = 1
H[2,3] = 0
H[3,1] = - sin(θ)
H[3,2] = 0
H[3,3] = cos(θ)
println("H:", H)
H = [cos(θ) 0 sin(θ); 0 1 0; - sin(θ) 0  cos(θ)]

r = R * r_icrs - dGC * xGC
rGC = H * r

println("distance:", d)
println("r_icrs:", r_icrs)
println("r:", r)
println("rGC:", rGC)

#+++++++++++++ velocities ++++++++++++++++++++++++
c1 = ICRSCoords(α, δ)
c2 = convert(GalCoords, c1)
l = c2.l
b = c2.b
println("c1:", c1, "\tc2:", c2)

#transform angular velocities to the galactic frame
C₁ = sin(δGC) * cos(δ) - cos(δGC) * sin(δ) * cos(α - αGC)
C₂ = cos(δGC) * sin(α - αGC)
cosb = sqrt(C₁^2 + C₂^2)

C = Matrix(undef, 2, 2)
C[1,1] = C₁
C[1,2] = C₂
C[2,1] = - C₂
C[2,2] = C₁

μ = [pmra, pmdec]
μ = C * μ / cosb
println(μ)
μl = μ[1]
μb = μ[2]

Vl = 4.74 * μl * d / 1000.0
Vb = 4.74 * μb * d / 1000.0
Vᵣ = radial_velocity

U = Vᵣ * cos(l) * cos(b) - Vl * sin(l) - Vb * cos(l) * sin(b)
V = Vᵣ * sin(l) * cos(b)+ Vl * cos(l) - Vb * sin(l) * sin(b)
W = Vᵣ * sin(b) + Vb * cos(b)
println("U:", U, "\tV:", V, "\tW:", W)
#println("adj. U:", U+vSUN[1], "\tadj. V:", V+vSUN[2], "\tadj. W:", W+vSUN[3])

VGC = [28.46188771, 157.93916086, 17.65189313]
#X = H * R * V_icrs - VGC
#X = R * V_icrs - inv(H) * VGC
println("VGC1 = ",VGC)

#X=[3701.88, -15812.7, 2195.97]
#VGC = H * R * V_icrs - X
VGC = H * R * V_icrs
println("VGC2 = ", VGC)


text:
"""
Package: Forecast

    function d(x::{AbstractVector,AbstractArray,DataFrame},
               or::Int=1,
               la::Int=1;
               center::Bool=false)

Return Laged differences of a given or for Vector, Array and TimeSeries.

# Arguments
- `x`: Vector or Array of data.
- `or`: Order of the differences; number of recursive iterations on the same vector/array.
- `la`: Lag for the difference.
- `center`: Center the result in the response using Missing values.

# Returns
Laged differences Vector or Array of a given order.

# Examples
```julia-repl
julia> x = [1,2,3,4,5];
julia> d(x)
 d(x)
4-element Vector{Int64}:
 1
 1
 1
 1

julia> d(x,2)
3-element Vector{Int64}:
 0
 0
 0

julia> d(x,1,2)
3-element Vector{Int64}:
 2
 2
 2

julia> x = reshape(collect(1:20),10,2);

julia> d(x,2,2)
6×2 Matrix{Int64}:
 0  0
 0  0
 0  0
 0  0
 0  0
 0  0

julia> d(d(x,1,2),1,2) == d(x,2,2)
```
"""
function d(x::AbstractArray,
           or::Int=1,
           la::Int=1;
           center::Bool=false)

    (la == 0) | (or == 0) && (return x)

    n = size(x,1)
    nv = size(x,2)

    @assert 0 <= la & la <= (n-1) "Lag must be larger or equal to size(x,1)"

    dx = x
    for i in or:-1:1
        dx = (dx .- circshift(dx,la))[1+la:end,:]
    end

    if center
        nl = n-size(dx,1)
        dx = vcat(dx, Array{Any}(missing,nl,size(x,2)))
        dx = circshift(dx,(div(nl,2),0))
    end

    return size(dx,2) == 1 ? dx[:,1] : dx

end

function d(df::DataFrame,
           or::Int=1,
           la::Int=1;
           center::Bool=false)

    df = tots(df)
    dfts = df[:,1:1]
    dfx = df[:,2:end]
    dnames = "d[" * string(or) * "," * string(la) * "]_" .* names(dfx)

    dx = d(Array(dfx), or, la; center=center)

    dfdx = DataFrame(reshape(dx,size(dx,1),size(dx,2)),dnames)
    ddfts = df[end-size(dx,1)+1:end,1:1]

    return hcat(ddfts,dfdx)

end


text:
function compute_β(params, α, n)
    I = compute_int_I(params, α, α, n)

    1 / √(2 * sum(I))
end


text:
# aerodyne_status.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.3.1
# 13.12.2016
# Last Edit: 10.03.2020

"""# AerodyneStatus

Parse the StatusW column in an STC file into this custom type

---

### Example

`statuses = AerodyneStatus.(Dstc[:StatusW])`\n
`valve1 = getfield.(statuses,:Valve1)` # Get all 'Valve1' field values from a vector of statuses, Dstc

---

#### Fields (all contain boolean values)\n
* **AutoBG** - Automatic Background
* **AutoCal** - Automatic Calibration
* **FrequencyLock**
* **BinomialFilter**
* **AltMode** - Alternative Mode
* **GuessLast** - First guess is last guess from previous curve fit
* **PowerNorm** - Power Normalization
* **ContRefLock** - Continuous Reference Lock
* **AutoSpectSave** - Automatic Spectrum Save
* **PressureLock**
* **WriteData**
* **RS232** - RS232 output active/inactive
* **ElectronicBGSub** - Subtract electronic background noise
* **Valve 1** - Active/inactive
* **Valve 2** - Active/inactive
* **Valve 3** - Active/inactive
* **Valve 4** - Active/inactive
* **Valve 5** - Active/inactive
* **Valve 6** - Active/inactive
* **Valve 7** - Active/inactive
* **Valve 8** - Active/inactive
"""
mutable struct AerodyneStatus
	AutoBG::Bool
	AutoCal::Bool
	FrequencyLock::Bool
	BinomialFilter::Bool
	AltMode::Bool
	GuessLast::Bool
	PowerNorm::Bool
	ContRefLock::Bool

	AutoSpectSave::Bool
	PressureLock::Bool
	#b11::Bool
	#b12::Bool
	WriteData::Bool
	RS232::Bool
	ElectronicBGSub::Bool
	#b16::Bool

	Valve1::Bool
	Valve2::Bool
	Valve3::Bool
	Valve4::Bool
	Valve5::Bool
	Valve6::Bool
	Valve7::Bool
	Valve8::Bool

	# Constructor for AerodyneStatus type
	function AerodyneStatus(val::Float64)
		AerodyneStatus(Int(val))
	end # End AerodyneStatus(val::Float64) constructor

	# Constructor for AerodyneStatus type
	function AerodyneStatus(val::Int)
		####################
		##  Parse Inputs  ##
		####################
		# Byte 1
		AutoBG = val & 2^(1-1) > 0
		AutoCal = val & 2^(2-1) > 0
		FrequencyLock = val & 2^(3-1) > 0
		BinomialFilter = val & 2^(4-1) > 0
		AltMode = val & 2^(5-1) > 0
		GuessLast = val & 2^(6-1) > 0
		PowerNorm = val & 2^(7-1) > 0
		ContRefLock = val & 2^(8-1) > 0

		# Byte 2
		AutoSpectSave = val & 2^(9-1) > 0
		PressureLock = val & 2^(10-1) > 0
		#b11 = val & 2^(11-1) > 0
		#b12 = val & 2^(12-1) > 0
		WriteData = val & 2^(13-1) > 0
		RS232 = val & 2^(14-1) > 0
		ElectronicBGSub = val & 2^(15-1) > 0
		#b16 = val & 2^(16-1) > 0

		# Byte 3
		Valve1 = val & 2^(17-1) > 0
		Valve2 = val & 2^(18-1) > 0
		Valve3 = val & 2^(19-1) > 0
		Valve4 = val & 2^(20-1) > 0
		Valve5 = val & 2^(21-1) > 0
		Valve6 = val & 2^(22-1) > 0
		Valve7 = val & 2^(23-1) > 0
		Valve8 = val & 2^(24-1) > 0

		new(AutoBG,AutoCal,FrequencyLock,BinomialFilter,AltMode,GuessLast,PowerNorm,ContRefLock,AutoSpectSave,PressureLock,WriteData,RS232,ElectronicBGSub,Valve1,Valve2,Valve3,Valve4,Valve5,Valve6,Valve7,Valve8)
	end # End of constructor
end # End of type

function Base.show(io::IO, status::AerodyneStatus)
	#println("| AutoBG | AutoCal | Frequency Lock | Binomial Filter | AltMode | Guess Last | Power Norm. | Cont. Ref. Lock | Auto Spec. Save | Press. Lock | Write Data | RS232 | Elec. BG Sub. | Valve 1 | Valve 2 | Valve 3 | Valve 4 | Valve 5 | Valve 6 | Valve 7 | Valve 8 |")
	status.Valve1 ? print("| Valve1") : nothing
	status.Valve2 ? print("| Valve2") : nothing
	status.Valve3 ? print("| Valve3") : nothing
	status.Valve4 ? print("| Valve4") : nothing
	status.Valve5 ? print("| Valve5") : nothing
	status.Valve6 ? print("| Valve6") : nothing
	status.Valve7 ? print("| Valve7") : nothing
	status.Valve8 ? print("| Valve8") : nothing
	status.WriteData ? print("| Write Data") : nothing
	status.RS232 ? print("| RS232 Output") : nothing
	status.GuessLast ? print("| Guess Last Fit") : nothing
	status.ContRefLock ? print("| Continuous Reference Lock") : nothing
	status.AutoSpectSave ? print("| Auto Spectrum Save") : nothing
	status.AutoBG ? print("| Auto Background") : nothing
	status.AutoCal ? print("| Auto Calibration") : nothing
	status.FrequencyLock ? print("| Frequency Lock") : nothing
	status.BinomialFilter ? print("| Binomial Filter") : nothing
	status.AltMode ? print("| AltMode") : nothing
	status.PowerNorm ? print("| Power Normalization") : nothing
	status.PressureLock ? print("| Pressure Lock") : nothing
	status.ElectronicBGSub ? print("| Electronic Background Subtraction") : nothing
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
module tensorboard
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module hparams
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("api_pb.jl")
    include("hparams_util_pb.jl")
    include("plugin_data_pb.jl")
  end
end


text:
# If we take 47, reverse and add, 47 + 74 = 121, which is palindromic.
#
# Not all numbers produce palindromes so quickly. For example,
#
# 349 + 943 = 1292,
# 1292 + 2921 = 4213
# 4213 + 3124 = 7337
#
# That is, 349 took three iterations to arrive at a palindrome.
#
# Although no one has proved it yet, it is thought that some numbers, like 196,
# never produce a palindrome. A number that never forms a palindrome through
# the reverse and add process is called a Lychrel number. Due to the
# theoretical nature of these numbers, and for the purpose of this problem, we
# shall assume that a number is Lychrel until proven otherwise. In addition you
# are given that for every number below ten-thousand, it will either (i) become
# a palindrome in less than fifty iterations, or, (ii) no one, with all the
# computing power that exists, has managed so far to map it to a palindrome. In
# fact, 10677 is the first number to be shown to require over fifty iterations
# before producing a palindrome: 4668731596684224866951378664 (53 iterations,
# 28-digits).
#
# Surprisingly, there are palindromic numbers that are themselves Lychrel
# numbers; the first example is 4994.
#
# How many Lychrel numbers are there below ten-thousand?

using ProjectEulerSolutions

# Cycle through list from 11 to n to calculate if a value and its reverse sum
# is a palindrome or not.  Trying a dynamic programming-type answer would not
# help much since the span of the possible numbers to store is quite large.
function p055solution(n_max::Integer=100)::Integer
    count = 0
    count_thresh = 50
    # We start with Int128, but we might do better with BigInt in some choices
    # of n_max.
    if n_max < 11_000
        n_start = Int128(10) # Single digit numbers are already palindromes
    else
        n_start = BigInt(10)
    end
    for n in n_start:n_max
        index = 0
        n2 = reverse_sum(n)
        while(!ispalindrome_integer(n2))
            n2 = reverse_sum_integer(n2)
            index += 1
            if (index == count_thresh)
                count += 1
                break
            end
        end
    end

    return count
end

p055 = Problems.Problem(p055solution)

Problems.benchmark(p055, 10_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using OffsetArrays
using Chain
using Utils

function parse_file(f)
    collect(eachline(f))
end

function load()
    open(joinpath(@__DIR__, "24.txt"), "r") do f
        parse_file(f)
    end
end

problem = load()

function ex1(problem)
    dim = 128 | 1
    @assert dim % 2 != 0 "Later code expects odd here."

    black = falses(2dim, 2dim)
    black = OffsetMatrix(black, 1-dim:dim, 1-dim:dim)
    for line in problem
        cur = 1
        c = r = 0
        while cur <= length(line)
            if line[cur] == 'w'
                c -= 2
            elseif line[cur] == 'e'
                c += 2
            else
                if line[cur] == 'n'
                    r -= 1
                elseif line[cur] == 's'
                    r += 1
                else
                    error("shouldn't get here (n or s)")
                end
                cur += 1
                if line[cur] == 'w'
                    c -= 1
                elseif line[cur] == 'e'
                    c += 1
                else
                    error("shouldn't get here (w or e)")
                end
            end
            cur += 1
        end
        @assert abs(r) % 2 == abs(c) % 2 "row and columns need to have the same LSB"
        black[r, c] = !black[r, c]
    end

    (count = count(values(black)), field = black)
end

function black_neighbors(field, r, c)
    @assert abs(r) % 2 == abs(c) % 2 "row and columns need to have the same LSB"
    coords = [
        (-1, -1),
        (-1, 1),
        (0, -2),
        (0, 2),
        (1, -1),
        (1, 1),
    ]
    @chain coords begin
        map(tup -> (r + tup[1], c + tup[2]), _)
        map(((rr, cc),) -> field[rr,cc], _)
        count
    end
end

function step(field)
    prev = copy(field)
    axr, axc = axes(prev)
    for r in firstindex(axr)+1:lastindex(axr)-1
        if r % 2 == 0
            range = firstindex(axc)+2:2:lastindex(axc)-2
            @assert 0 in range
        else
            range = firstindex(axc)+3:2:lastindex(axc)-3
            @assert 1 in range
        end
        for c in range
            ns = black_neighbors(prev, r, c)
            if prev[r,c] && (ns == 0 || ns > 2)
                field[r,c] = false
            elseif !prev[r,c] && ns == 2
                field[r,c] = true
            end
        end
    end
end

function ex2(field)
    for i in 1:100
        step(field)
    end
    count(field)
end

answer = ex1(problem)
println("Black side up initially: ", answer.count)
println("Black side up after 100 days: ", ex2(answer.field))

using Test

@testset "Lobby Layout" begin
    @testset "example 1" begin
        input1 = """
        sesenwnenenewseeswwswswwnenewsewsw
        neeenesenwnwwswnenewnwwsewnenwseswesw
        seswneswswsenwwnwse
        nwnwneseeswswnenewneswwnewseswneseene
        swweswneswnenwsewnwneneseenw
        eesenwseswswnenwswnwnwsewwnwsene
        sewnenenenesenwsewnenwwwse
        wenwwweseeeweswwwnwwe
        wsweesenenewnwwnwsenewsenwwsesesenwne
        neeswseenwwswnwswswnw
        nenwswwsewswnenenewsenwsenwnesesenew
        enewnwewneswsewnwswenweswnenwsenwsw
        sweneswneswneneenwnewenewwneswswnese
        swwesenesewenwneswnwwneseswwne
        enesenwswwswneneswsenwnewswseenwsese
        wnwnesenesenenwwnenwsewesewsesesew
        nenewswnwewswnenesenwnesewesw
        eneswnwswnwsenenwnwnwwseeswneewsenese
        neswnwewnwnwseenwseesewsenwsweewe
        wseweeenwnesenwwwswnew
        """

        example1 = parse_file(IOBuffer(input1))

        answer = ex1(example1)
        @test ex1(example1).count == 10
        @test ex2(answer.field) == 2208
    end

    @testset "results" begin
        answer = ex1(problem)
        @test answer.count == 332
        @test ex2(answer.field) == 3900
    end
end

max_stars_repo_name: CmdQ/AoC


text:
function get_eop_history(obj::GTODProvider)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_non_interpolating_provider(obj::GTODProvider)
    return jcall(obj, "getNonInterpolatingProvider", GTODProvider, ())
end

function get_transform(obj::GTODProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::GTODProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
module Fractory

using PyPlot
import PyPlot: draw
import Base: show,(*),size

include("AffineMap.jl")
include("IFS.jl")
include("PixelPicture.jl")
include("Fractal.jl")

end


text:
# Autogenerated wrapper script for NASM_jll for i686-w64-mingw32
export nasm, ndisasm

JLLWrappers.@generate_wrapper_header("NASM")
JLLWrappers.@declare_executable_product(nasm)
JLLWrappers.@declare_executable_product(ndisasm)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        nasm,
        "bin\\nasm.exe",
    )

    JLLWrappers.@init_executable_product(
        ndisasm,
        "bin\\ndisasm.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    TrustRegion{H, M, D, R<:Ref, E} <: MultivariateAlgorithm

`TrustRegion` builds a quadratic model within a radius from the current point to approximate
the objective function and uses the model to choose the next step.

# Fields
- `hessian::H`: the method to compute or approximate the Hessian
- `method::M`: the method to compute the appropriate next step
- `Δ₀::D`: the initial radius
- `Δₘₐₓ::D`: the maximum radius
- `η::E`: the threshold to control when `TrustRegion` will take a step
- `ηₛ::E`: the threshold to control when `TrustRegion` shrinks the current radius
- `ηₑ::E`: the threshold to control when `TrustRegion` expands the current radius
- `σₛ::S`: the factor used to shrink the trust region radius
- `σₑ::S`: the factor used to expand the trust region radius
- `Δ::R`: a `Ref` containing the current radius
"""
struct TrustRegion{H, M, D, E, S, R} <: MultivariateAlgorithm
    hessian::H
    method::M
    Δ₀::D
    Δₘₐₓ::D
    η::E
    ηₛ::E
    ηₑ::E
    σₛ::S
    σₑ::S
    Δ::R

    function TrustRegion(hessian, method, Δ₀::D, Δₘₐₓ::D, η, ηₛ, ηₑ, σₛ, σₑ) where {D}
        0 ≤ η < 0.5 || error("`η` must be in the interval [0, 0.5)")
        0 ≤ ηₛ < 0.5 || error("`ηₛ` must be in the interval [0, 0.5)")
        ηₑ ≥ 0.5 || error("`ηₑ` must be larger than or equal to 0.5")
        0 < σₛ < 1 || error("Shrinking factor `σₛ` must be in the interval (0, 1)")
        σₑ > 1 || error("Expanding factor `σₑ` must be larger than 1")
        Δ = Ref{D}()
        H = typeof(hessian)
        M = typeof(method)
        E = typeof(η)
        S = typeof(σₛ)
        R = typeof(Δ)
        new{H, M, D, E, S, R}(hessian, method, Δ₀, Δₘₐₓ, η, ηₛ, ηₑ, σₛ, σₑ, Δ)
    end
end

"""
    TrustRegion(; kwargs...)

Initiate `TrustRegion` algorithm.

# Keywords
- `hessian=Newton()`: the method to compute or approximate the Hessian
- `method=TwoDimSubspace()`: the method to compute the appropriate next step
- `delta=5.0`: the initial radius
- `deltamax=100.0`: the maximum radius
- `threshold=0.1`: the threshold to control when `TrustRegion` will take a step
- `shrinkthreshold=0.25`: the threshold to control when `TrustRegion` shrinks the current radius
- `expandthreshold=0.75`: the threshold to control when `TrustRegion` expands the current radius
- `shrinkfactor=0.25`: the factor used to shrink the trust region radius
- `expandfactor=2.0`: the factor used to expand the trust region radius
"""
function TrustRegion(;
        hessian=Newton(),
        method=TwoDimSubspace(),
        delta=5.0,
        deltamax=100.0,
        threshold=0.1,
        shrinkthreshold=0.25,
        expandthreshold=0.75,
        shrinkfactor=0.25,
        expandfactor=2.0)
    D = promote_type(typeof(delta / 1), typeof(deltamax / 1))
    E = promote_type(typeof(threshold), typeof(shrinkthreshold), typeof(expandthreshold))
    S = promote_type(typeof(shrinkfactor), typeof(expandfactor))
    TrustRegion(hessian, method, D(delta), D(deltamax), E(threshold), E(shrinkthreshold), E(expandthreshold),
        S(shrinkfactor), S(expandfactor))
end

function Base.show(io::IO, tr::TrustRegion)
    println(io, "Trust Region Algorithm:")
    println(io, "  • Hessian approximation method: ", tr.hessian)
    println(io, "  • Subproblem optimization method: ", tr.method)
    println(io, "  • Initial trust region radius: ", tr.Δ₀)
    println(io, "  • Maximum trust region radius: ", tr.Δₘₐₓ)
    println(io, "  • Step threshold: ", tr.η)
    println(io, "  • Trust region shrinkage threshold: ", tr.ηₛ)
    println(io, "  • Trust region expansion threshold: ", tr.ηₑ)
    println(io, "  • Trust region shrinkage factor: ", tr.σₛ)
    print(io, "  • trust region expansion factor: ", tr.σₑ)
end

order(tr::TrustRegion) = order(tr.hessian)


text:
using Optim, StatsFuns

function quap(loglikf, x0, lower, upper)
  res = optimize(loglikf, lower, upper, x0, Fminbox(GradientDescent()))
  (Optim.minimizer(res), res)
end