text:
module testmod

export hi, func1

mutable struct hi
a::Int64
end

function func1(x::hi)
    return x.a^2
end

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
PajaritoExtras examples and script utilities.
"""
module Examples

using LinearAlgebra
# delete later, affects qr. see https://github.com/JuliaLang/julia/pull/40623
if VERSION < v"1.7.0-DEV.1188"
    const ColumnNorm = Val{true}
end
using SparseArrays

import Random
# this is a workaround for randn's lack of support for BigFloat
Random.randn(R::Type{BigFloat}, dims::Vararg{Int, N} where {N}) = R.(randn(dims...))
function Random.randn(R::Type{Complex{BigFloat}}, dims::Vararg{Int, N} where {N})
    return R.(randn(ComplexF64, dims...))
end

using Test
using Printf
import DataFrames
import CSV
import DataStructures: OrderedDict

import Hypatia
import Hypatia.PolyUtils
import Hypatia.Cones
import Hypatia.Solvers

import MOIPajarito
import PajaritoExtras

import JuMP
const MOI = JuMP.MOI
const MOIU = MOI.Utilities
const JuMPScalar = JuMP.AbstractJuMPScalar

abstract type ExampleInstance end

include("benchmark_utils.jl")
include("JuMP_utils.jl")

const OPT_OR_LIMIT = [MOI.OPTIMAL, MOI.TIME_LIMIT, MOI.ITERATION_LIMIT]

# TODO maybe move to Hypatia
include("WSOS_utils.jl")
include("specnuc_utils.jl")
include(joinpath(pkgdir(Hypatia), "examples", "spectral_functions_JuMP.jl"))
include("SOS2_utils.jl")

include(joinpath(@__DIR__, "..", "test", "model_utils.jl"))

# load all examples
function load_examples(examples::Vector{String})
    for ex in examples
        include(joinpath(@__DIR__, ex, "model.jl"))
    end
end

end

ext: jl


text:
include("../common.jl")

version = v"1.5.1"

name, version, sources, script, platforms, products, dependencies = configure(v"1.5.1")

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

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
# Run this if you do not already have Keras setup on your machine
# Sets up PyCall to use a miniconda distribution
Pkg.add("Conda")
ENV["PYTHON"] = ""
Pkg.build("PyCall")

using Conda
push!(Conda.CHANNELS, "https://conda.anaconda.org/jaikumarm")
Conda.add("keras")


text:
"""
Least-squares spectral estimation toolbox.

For help, see README.md at https://github.com/baggepinnen/LPVSpectral.jl and

[Fredrik Bagge Carlson, Anders Robertsson, Rolf Johansson: "Linear Parameter-Varying Spectral Decomposition". In: 2017 American Control Conference 2017.]
(http://lup.lub.lu.se/record/ac32368e-e199-44ff-b76a-36668ac7d595)
Available at: http://lup.lub.lu.se/record/ac32368e-e199-44ff-b76a-36668ac7d595

This module provide the functions
```
ls_spectral
tls_spectral
ls_windowpsd
ls_windowcsd
ls_cohere
ls_spectral_lpv
ls_sparse_spectral_lpv
ls_windowpsd_lpv
basis_activation_func
```
and re-exports the following from DSP.jl
```
export periodogram, welch_pgram, Windows
```
Periodogram types and SpectralExt type can be plotted using `plot(x::SpectralExt)`
"""
module LPVSpectral
using LinearAlgebra, Statistics, Printf
using DSP, FFTW, StatsBase
using RecipesBase

"""
LPV Spectral estimation result type.

See `ls_spectral_lpv` for additional help.

An object of this type can be plotted if `Plots.jl` is installed. Use regular Plots-syntax, with the additional attributes
```
normalization= :none / :sum / :max
normdim = :freq / :v # Only applies if normalization= :sum or :max
dims = 2 or 3 (default = 2)
```

Fields:
```
Y::AbstractVector
X::AbstractVector
V::AbstractVector
w
Nv
λ
coulomb::Bool
normalize::Bool
x                   # The estimated parameters
Σ                   # Covariance of the estimated parameters
```
"""
struct SpectralExt
    Y::AbstractVector
    X::AbstractVector
    V::AbstractVector
    w
    Nv
    λ
    coulomb::Bool
    normalize::Bool
    x
    Σ
end

include("utilities.jl")
include("windows.jl")
include("mel.jl")
include("plotting.jl")
include("lsfft.jl")
include("autocov.jl")
include("lasso.jl")

# Functions
export ls_spectral,
tls_spectral,
ls_sparse_spectral,
ls_windowpsd,
ls_windowcsd,
ls_cohere,
ls_spectral_lpv,
ls_sparse_spectral_lpv,
ls_windowpsd_lpv,
basis_activation_func,
SpectralExt,
psd,
detrend,
detrend!,
autocov,
autocor

export melspectrogram, mfcc, mel

# Re-export
export Windows, Windows2, Windows3, mapwindows

# ComplexNormal
export ComplexNormal

export cn_V2ΓC,cn_V2ΓC,cn_Vxx,cn_Vyy,cn_Vxy,cn_Vyx,cn_fVxx,cn_fVyy,cn_fVxy,
cn_fVyx,cn_Vs,cn_V,cn_fV,Σ,pdf,affine_transform, rand

using Requires
function __init__()
    @require ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537" include("lasso.jl")
end

end # module


text:
################################################################################
#                  MakieGallery.jl documentation build script                  #
################################################################################
##############################
#      Generic imports       #
##############################

using Documenter, Markdown, Pkg, Random, FileIO, JSON, HTTP, GitHub
using D3TypeTrees, D3Trees

##############################
#      Specific imports      #
##############################

using MakieGallery, AbstractPlotting

import AbstractPlotting: to_string

using MakieGallery: eval_examples, generate_thumbnail, master_url,
                    print_table, download_reference,
                    @cell, @block, @substep

################################################################################
#                              Utility functions                               #
################################################################################

################################################################################
#                                    Setup                                     #
################################################################################

MakieGallery.current_ref_version[] = "master"

empty!(MakieGallery.plotting_backends)
append!(MakieGallery.plotting_backends, ["Makie"])

cd(@__DIR__)
database = MakieGallery.load_database()

pathroot  = normpath(@__DIR__, "..")
docspath  = joinpath(pathroot, "docs")
srcpath   = joinpath(docspath, "src")
buildpath = joinpath(docspath, "build")
genpath   = joinpath(srcpath, "generated")
mediapath = download_reference()

mkpath(genpath)

################################################################################
#                          Syntax highlighting theme                           #
################################################################################

@info("Writing highlighting stylesheet")
open(joinpath(srcpath, "assets", "syntaxtheme.css"), "w") do io
    MakieGallery.Highlights.stylesheet(io, MIME("text/css"), MakieGallery.DEFAULT_HIGHLIGHTER[])
end

################################################################################
#                      Automatic Markdown page generation                      #
################################################################################

########################################
#     Plotting functions overview      #
########################################

@info("Generating functions overview")
path = joinpath(srcpath, "functions-overview.md")
srcdocpath = joinpath(srcpath, "src-functions.md")

plotting_functions = (
    AbstractPlotting.atomic_functions..., contour, arrows,
    barplot, poly, band, slider, vbox, hbox
)

open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plotting functions overview")
    src = read(srcdocpath, String)
    println(io, src, "\n")
    for func in plotting_functions
        fname = to_string(func)
        println(io, "## `$fname`\n")
        println(io, "```@docs")
        println(io, "$fname")
        println(io, "```\n")
        # add previews of all tags related to function
        for example in database
            fname in example.tags || continue
            base_path = joinpath(mediapath, string(example.unique_name))
            thumb = master_url(mediapath, joinpath(base_path, "media", "thumb.jpg"))
            code = master_url(mediapath, joinpath(base_path, "index.html"))
            src_lines = example.file_range
            println(io, """[![library lines $src_lines]($thumb)]($code)""")
        end
        println(io, "\n")
    end
end

cd(docspath)

########################################
#       Plot attributes overview       #
########################################

# automatically generate an overview of the plot attributes (keyword arguments), using a source md file
@info("Generating attributes page")
include("../src/plot_attr_desc.jl")
path = joinpath(srcpath, "plot-attributes.md")
srcdocpath = joinpath(srcpath, "src-plot-attributes.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot attributes")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "## List of attributes")
    print_table(io, plot_attr_desc)
end

########################################
#       Axis attributes overview       #
########################################

# automatically generate an overview of the axis attributes, using a source md file
@info("Generating axis page")
path = joinpath(srcpath, "axis.md")
srcdocpath = joinpath(srcpath, "src-axis.md")
include("../src/Axis2D_attr_desc.jl")
include("../src/Axis3D_attr_desc.jl")

open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    src = read(srcdocpath, String)
    println(io, src)
    print(io)
    # Axis2D section
    println(io, "## `Axis2D`")
    println(io, "### `Axis2D` attributes groups")
    print_table(io, Axis2D_attr_desc)
    print(io)
    for (k, v) in Axis2D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end
    # Axis3D section
    println(io, "## `Axis3D`")
    println(io, "### `Axis3D` attributes groups")
    print_table(io, Axis3D_attr_desc)
    print(io)
    for (k, v) in Axis3D_attr_groups
        println(io, "#### `:$k`\n")
        print_table(io, v)
        println(io)
    end

    # Examples on the bottom of the page
    println(io, """

    ### Examples

    @example_database("Unicode Marker")

    @example_database("Axis + Surface")

    @example_database("Axis theming")
    """)
end

########################################
#     Function signatures overview     #
########################################

# automatically generate an overview of the function signatures, using a source md file
@info("Generating signatures page")
path = joinpath(srcpath, "signatures.md")
srcdocpath = joinpath(srcpath, "src-signatures.md")
open(path, "w") do io
    !ispath(srcdocpath) && error("source document doesn't exist!")
    println(io, "# Plot function signatures")
    src = read(srcdocpath, String)
    println(io, src)
    print(io, "\n")
    println(io, "```@docs")
    println(io, "convert_arguments")
    println(io, "```\n")

    println(io, "See [Plot attributes](@ref) for the available plot attributes.")
end

########################################
#          Colormap reference          #
########################################

@info "Generating colormap reference"

MakieGallery.generate_colorschemes_markdown(; GENDIR = genpath)

########################################
#              Type trees              #
########################################

@info "Generating type trees"
open(joinpath(srcpath, "typetrees.md"), "w") do f

    println(f, "# Type Trees")
    println(f)

    for typ in [AbstractPlotting.Transformable, AbstractPlotting.AbstractCamera]
        println(f, "## $(Symbol(typ))")
        println(f)
        println(f, "```@raw html")
        show(f, MIME("text/html"), TypeTree(typ))
        println(f, "```")
        println(f)
    end

end
################################################################################
#                 Building HTML documentation with Documenter                  #
################################################################################

@info("Running `makedocs` with Documenter.")

makedocs(
    modules = [AbstractPlotting],
    doctest = false, clean = true,
    format = Documenter.HTML(
        prettyurls = false,
        assets = [
            "assets/favicon.ico",
            "assets/syntaxtheme.css"
        ],
    ),
    sitename = "Makie.jl",
    expandfirst = [
        "basic-tutorials.md",
        "statsmakie.md",
        "animation.md",
        "interaction.md",
        "functions-overview.md",
        "scenes.md",
        "signatures.md",
        "plot-attributes.md",
        "generated/colors.md",
        "theming.md",
        "cameras.md",
        "backends.md",
        "axis.md",
        "recipes.md",
        "output.md"
    ],
    pages = Any[
        "Home" => "index.md",
        "Basics" => [
            "basic-tutorials.md",
            "statsmakie.md",
            "animation.md",
            "interaction.md",
            "functions-overview.md",
            "animation.md",
        ],
        "Documentation" => [
            "scenes.md",
            "axis.md",
            "convenience.md",
            "signatures.md",
            "plot-attributes.md",
            "generated/colors.md",
            "theming.md",
            "cameras.md",
            "recipes.md",
            "output.md",
            "backends.md",
            "troubleshooting.md"
        ],
        "Developer Documentation" => [
            "why-makie.md",
            "devdocs.md",
            "gallery.md",
            "AbstractPlotting Reference" => "abstractplotting_api.md",
            "Type Trees" => "typetrees.md"
        ],
    ]
)

################################################################################
#                           Deploying documentation                            #
################################################################################

using Documenter, Base64

ENV["DOCUMENTER_DEBUG"] = "true"

# do this only if local, otherwise let Documenter handle it
if !haskey(ENV, "CI") && !haskey(ENV, "GITHUB_TOKEN")
    @info "Manually deploying - setting environment variables!"
    # Workaround for when deploying locally and silly Windows truncating the env variable
    # on the CI these should be set!
    ENV["CI"] = "no"
    ENV["TRAVIS"] = :lolno
    ENV["TRAVIS_BRANCH"] = "master"
    ENV["TRAVIS_PULL_REQUEST"] = "false"
    ENV["TRAVIS_REPO_SLUG"] = "github.com/JuliaPlots/MakieGallery.jl"
    ENV["TRAVIS_TAG"] = ""
    ENV["TRAVIS_OS_NAME"] = ""
    ENV["TRAVIS_JULIA_VERSION"] = ""
    # ENV["PATH"] = string(ENV["PATH"], Sys.iswindows() ? ";" : ":", Conda.SCRIPTDIR)
    ENV["DOCUMENTER_KEY"] = get(ENV, "DOCUMENTER_KEY", readchomp(joinpath(homedir(), "documenter.key")))
end

############################################
# Set up for pushing preview docs from PRs #
############################################

cd(@__DIR__)

Documenter.authentication_method(::Documenter.GitHubActions) = Documenter.SSH

Base.invokelatest(deploydocs,
    repo = "github.com/JuliaPlots/MakieGallery.jl",
    push_preview = true
)


text:
function EclipseDetector(arg0::PVCoordinatesProvider, arg1::jdouble, arg2::OneAxisEllipsoid)
    return EclipseDetector((PVCoordinatesProvider, jdouble, OneAxisEllipsoid), arg0, arg1, arg2)
end

function g(obj::EclipseDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_occulted(obj::EclipseDetector)
    return jcall(obj, "getOcculted", PVCoordinatesProvider, ())
end

function get_occulted_radius(obj::EclipseDetector)
    return jcall(obj, "getOccultedRadius", jdouble, ())
end

function get_occulting(obj::EclipseDetector)
    return jcall(obj, "getOcculting", OneAxisEllipsoid, ())
end

function get_total_eclipse(obj::EclipseDetector)
    return jcall(obj, "getTotalEclipse", jboolean, ())
end

function with_penumbra(obj::EclipseDetector)
    return jcall(obj, "withPenumbra", EclipseDetector, ())
end

function with_umbra(obj::EclipseDetector)
    return jcall(obj, "withUmbra", EclipseDetector, ())
end


text:
module IntegralHistograms

using StatsBase, Images

include("algorithm.jl")

export
	# main types and functions
	integralimage,
	integralhistogram

end # module


text:
chain_swap = Chain()

push_gate!(chain_swap, Val(:SWAP), 1, 3)
push_gate!(chain_swap, Val(:SWAP), 2, 4)
push_gate!(chain_swap, Val(:SWAP), 2, 3)
push_gate!(chain_swap, Val(:SWAP), 1, 4)
push_gate!(chain_swap, Val(:SWAP), 3, 4)
push_gate!(chain_swap, Val(:SWAP), 4, 3)
push_gate!(chain_swap, Val(:SWAP), 2, 3)

ZXCalculus.simplify_swap!(chain_swap)

@test length(chain_swap) == 9


text:
"""
    AbstractSpacing

Spacing discretization scheme supertype
"""
abstract type AbstractSpacing end

"""
    Uniform()

Uniform discretization scheme.
"""
struct Uniform <: AbstractSpacing end

"""
    Sine()

Sine-spaced discretization scheme.  Using sine-spacing on the right half of a
wing effectively results in cosine spacing once symmetry is applied.
"""
struct Sine <: AbstractSpacing end

"""
    Cosine()

Cosine-spaced discretization scheme.  This is typically one of the most accurate
spacing schemes for spanwise spacing.
"""
struct Cosine <: AbstractSpacing end

"""
   linearinterp(eta, rstart, rend)

Linearly interpolate between `rstart` and `rend` where eta is the fraction
between 0 (rstart) and 1 (rend)
"""
linearinterp(eta, rstart, rend) = (1-eta)*rstart + eta*rend

"""
    spanwise_spacing(n, spacing::AbstractSpacing)

Distribute `n` panel endpoints and `n-1` panel midpoints on the interval between
0 and 1 according to the discretization strategy in `spacing`.
"""
spanwise_spacing

# uniform
function spanwise_spacing(n, ::Uniform)

    eta = range(0, 1.0, length=n)
    eta_mid = linearinterp(0.5, eta[1:end-1], eta[2:end])

    return eta, eta_mid
end

# sine
function spanwise_spacing(n, ::Sine)

    theta = range(0, pi/2, length=n)
    eta = sin.(theta)

    # note that control points are also placed with cosine spacing as this improves accuracy tremendously
    theta_mid = linearinterp(0.5, theta[1:end-1], theta[2:end])
    eta_mid = sin.(theta_mid)

    return eta, eta_mid
end

# cosine
function spanwise_spacing(n, ::Cosine)
    theta = range(0, pi, length=n)
    eta = (1.0 .- cos.(theta))/2.0

    # note that control points are also placed with cosine spacing as this improves accuracy tremendously
    theta_mid = linearinterp(0.5, theta[1:end-1], theta[2:end])
    eta_mid = (1.0 .- cos.(theta_mid))/2.0

    return eta, eta_mid
end

"""
    chordwise_spacing(n, spacing::AbstractSpacing)

Distribute `n` panel edge, `n-1` vortex, and `n-1` control point chordwise
locations on the interval between 0 and 1 according to the discretization
strategy in `spacing`.
"""
chordwise_spacing

# uniform
function chordwise_spacing(n, ::Uniform)

    eta = range(0, 1.0, length=n)

    eta_qtr = linearinterp(0.25, eta[1:end-1], eta[2:end])
    eta_thrqtr = linearinterp(0.75, eta[1:end-1], eta[2:end])

    return eta, eta_qtr, eta_thrqtr
end

# sine
function chordwise_spacing(n, ::Sine)

    theta = range(0, pi/2, length=n)
    eta = sin.(theta)

    eta_qtr = linearinterp(0.25, eta[1:end-1], eta[2:end])
    eta_thrqtr = linearinterp(0.75, eta[1:end-1], eta[2:end])

    return eta, eta_qtr, etq_thrqtr
end

# cosine
function chordwise_spacing(n, ::Cosine)

    theta = range(0, pi, length=n)
    eta = (1.0 .- cos.(theta))/2.0

    eta_qtr = linearinterp(0.25, eta[1:end-1], eta[2:end])
    eta_thrqtr = linearinterp(0.75, eta[1:end-1], eta[2:end])

    return eta, eta_qtr, eta_thrqtr
end

"""
    interpolate_grid(xyz, eta, interp; xdir=0, ydir=1)

Interpolates the grid `xyz` along direction `dir`

# Arguments
 - `xyz`: Grid of size (3, ni, nj)
 - `eta`: New (normalized) coordinates in direction `dir` (0 <= eta <= 1)
 - `interp`: Interpolation method of form `ypt = f(x,y,xpt)`
 - `xdir`: Independent variable direction, defaults to arc length
 - `ydir`: Dependent variable direction `xyz` (`i=1`, `j=2`)
"""
function interpolate_grid(xyz, eta, interp; xdir=0, ydir=1)

    ydim = ydir + 1

    ni = size(xyz, ydim)

    xyz = mapslices(xyz; dims=[1,ydim]) do xyz_i

        # get distance between each grid location
        if xdir == 0
            ds = [norm(xyz_i[:,k] - xyz_i[:,max(1,k-1)]) for k = 1:ni]
        else
            ds = [xyz_i[xdir,k] - xyz_i[xdir,max(1,k-1)] for k = 1:ni]
        end

        # create interpolation vector
        t = cumsum(ds)

        # normalize interpolation vector
        t /= t[end]

        # interpolate x, y, and z
        x = interp(t, xyz_i[1,:], eta)
        y = interp(t, xyz_i[2,:], eta)
        z = interp(t, xyz_i[3,:], eta)

        vcat(x',y',z')
    end

    return xyz
end

"""
    grid_to_surface_panels(xyz; mirror = false, fcore = (c, Δs) -> 1e-3)

Construct a set of panels with associated vortex rings given a potentially curved
lifting surface defined by a grid with dimensions (3, i, j) where `i` corresponds
to the chordwise direction (ordered from leading edge to trailing edge) and `j`
corresponds to the spanwise direction (ordered from left to right).  The leading
edge of each ring vortex will be placed at the 1/4 chord and the control point
will be placed at the 3/4 chord of each panel.

Return a grid with dimensions (3, i, j) containing the panel corners and a matrix
with dimensions (i, j) containing the generated panels.

# Keyword Arguments
- `mirror`:  mirror the geometry across the X-Z plane? defaults to `false`.
- `fcore`: function for setting the finite core size based on the chord length
       (in the x-direction) and/or the panel width (in the y/z directions).
       Defaults to `(c, Δs) -> 1e-3`
"""
function grid_to_surface_panels(xyz;
    mirror = false,
    fcore = (c, Δs) -> 1e-3)

    TF = eltype(xyz)

    nc = size(xyz, 2)-1 # number of chordwise panels
    ns = size(xyz, 3)-1 # number of spanwise panels
    N = (1+mirror)*nc*ns # total number of panels

    # check which side we're working with
    right_side = sum(xyz[2,:,:]) > 0

    # initialize output
    xyz_panels = copy(xyz)
    panels = Matrix{SurfacePanel{TF}}(undef, nc, (1+mirror)*ns)

    # populate each panel
    for j = 1:ns

        # get leading edge panel corners
        r1n = SVector(xyz[1,1,j], xyz[2,1,j], xyz[3,1,j]) # top left
        r2n = SVector(xyz[1,1,j+1], xyz[2,1,j+1], xyz[3,1,j+1]) # top right
        r3n = SVector(xyz[1,2,j], xyz[2,2,j], xyz[3,2,j]) # bottom left
        r4n = SVector(xyz[1,2,j+1], xyz[2,2,j+1], xyz[3,2,j+1]) # bottom right

        # also get chord length for setting finite core size
        cl = norm(xyz[:,end,j] - xyz[:,1,j])
        cr = norm(xyz[:,end,j+1] - xyz[:,1,j+1])
        c = (cl + cr)/2

        for i = 1:nc-1

            # xyz corners of current panel
            r1 = r1n # top left of panel
            r2 = r2n # top right of panel
            r3 = r3n # bottom left of panel
            r4 = r4n # bottom right of panel

            # xyz corners of next panel
            r1n = r3 # top left
            r2n = r4 # top right
            r3n = SVector(xyz[1,i+2,j], xyz[2,i+2,j], xyz[3,i+2,j]) # bottom left
            r4n = SVector(xyz[1,i+2,j+1], xyz[2,i+2,j+1], xyz[3,i+2,j+1]) # bottom right

            # top left corner of ring vortex
            rtl = linearinterp(0.25, r1, r3)

            # top right corner of ring vortex
            rtr = linearinterp(0.25, r2, r4)

            # center of top bound vortex
            rtc = (rtl + rtr)/2

            # bottom left corner of ring vortex
            rbl = linearinterp(0.25, r1n, r3n)

            # bottom right corner of ring vortex
            rbr = linearinterp(0.25, r2n, r4n)

            # center of bottom bound vortex
            rbc = (rbl + rbr)/2

            # control point
            rtop = linearinterp(0.5, r1, r2)
            rbot = linearinterp(0.5, r3, r4)
            rcp = linearinterp(0.75, rtop, rbot)

            # surface normal
            ncp = cross(rcp - rtr, rcp - rtl)
            ncp /= norm(ncp)

            # set finite core size
            Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
            core_size = fcore(c, Δs)

            # get chord length of current panel
            chord = norm((r1 + r2)/2 - (r3 + r4)/2)

            ip = i
            jp = mirror*right_side*ns + j
            panels[ip, jp] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp,
                ncp, core_size, chord)
        end

        # xyz corners of current panel
        r1 = r1n # top left
        r2 = r2n # top right
        r3 = r3n # bottom left
        r4 = r4n # bottom right

        # top left corner of ring vortex
        rtl = linearinterp(0.25, r1, r3)

        # top right corner of ring vortex
        rtr = linearinterp(0.25, r2, r4)

        # center of top bound vortex
        rtc = (rtl + rtr)/2

        # bottom left corner of ring vortex
        rbl = r3

        # bottom right corner of ring vortex
        rbr = r4

        # center of bottom bound vortex
        rbc = (rbl + rbr)/2

        # control point
        rtop = linearinterp(0.5, r1, r2)
        rbot = linearinterp(0.5, r3, r4)
        rcp = linearinterp(0.75, rtop, rbot)

        # surface normal
        ncp = cross(rcp - rtr, rcp - rtl)
        ncp /= norm(ncp)

        # set finite core size
        Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
        core_size = fcore(c, Δs)

        # get chord length of current panel
        chord = norm((r1 + r2)/2 - (r3 + r4)/2)

        ip = nc
        jp = mirror*right_side*ns + j
        panels[ip,jp] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp,
            core_size, chord)
    end

    # other side
    if mirror

        # first reflect xyz
        if right_side
            xyz_l = reverse(xyz_panels, dims=3)
            xyz_l[2,:,:] .= -xyz_l[2,:,:]
            xyz_r = xyz_panels
        else
            xyz_l = xyz_panels
            xyz_r = reverse(xyz_panels, dims=3)
            xyz_r[2,:,:] .= -xyz_r[2,:,:]
        end
        xyz_panels = cat(xyz_l, xyz_r[:,:,2:end], dims=3)

        # then populate remaining panels
        for j = 1:ns
            for i = 1:nc
                ip = i
                jp1 = right_side*ns + j
                jp2 = 2*ns - right_side*ns - j + 1
                panels[ip,jp2] = reflect(panels[ip,jp1])
            end
        end
    end

    return xyz_panels, panels
end

"""
    grid_to_surface_panels(xyz, ns, nc;
        mirror = false,
        fcore = (c, Δs) -> 1e-3,
        spacing_s = Cosine(),
        spacing_c = Uniform(),
        interp_s = (x, y, xpt) -> LinearInterpolation(x, y)(xpt),
        interp_c = (x, y, xpt) -> LinearInterpolation(x, y)(xpt))

Discretize a potentially curved lifting surface defined by a grid with dimensions
(3, i, j) where `i` corresponds to the chordwise direction (ordered from leading
edge to trailing edge) and `j` corresponds to the spanwise direction (ordered from
left to right) into `ns` spanwise and `nc` chordwise panels with associated vortex
rings according to the spanwise discretization scheme `spacing_s` and chordwise
discretization scheme `spacing_c`.  The bound vortex will be placed at the
1/4 chord and the control point will be placed at the 3/4 chord of each panel.

Return a grid with dimensions (3, i, j) containing the interpolated panel corners
and a matrix with dimensions (i, j) containing the generated panels.

# Arguments
 - `xyz`: grid of dimensions (3, i, j) where where `i` corresponds to the
    chordwise direction and `j` corresponds to the spanwise direction.
 - `ns`: number of spanwise panels
 - `nc`: number of chordwise panels

# Keyword Arguments
 - `mirror`:  mirror the geometry across the X-Z plane? defaults to `false`.
 - `fcore`: function for setting the finite core size based on the chord length
        (in the x-direction) and/or the panel width (in the y/z directions).
        Defaults to `(c, Δs) -> 1e-3`
 - `spacing_s`: spanwise discretization scheme, defaults to `Cosine()`
 - `spacing_c`: chordwise discretization scheme, defaults to `Uniform()`
 - `interp_s`: spanwise interpolation function, defaults to linear interpolation
 - `interp_c`: chordwise interpolation function, defaults to linear interpolation
"""
function grid_to_surface_panels(xyz, ns, nc;
    mirror = false,
    fcore = (c, Δs) -> 1e-3,
    spacing_s = Cosine(),
    spacing_c = Uniform(),
    interp_s = (x, y, xpt) -> LinearInterpolation(x, y)(xpt),
    interp_c = (x, y, xpt) -> LinearInterpolation(x, y)(xpt))

    TF = eltype(xyz)

    # get spanwise and chordwise spacing
    etas, etabar = spanwise_spacing(ns+1, spacing_s)
    etac, eta_qtr, eta_thrqtr = chordwise_spacing(nc+1, spacing_c)

    # add leading and trailing edge to grid
    eta_qtr = vcat(0.0, eta_qtr, 1.0)

    # interpolate grid first along `i` then along `j`
    xyz_edge = interpolate_grid(xyz, etac, interp_c, ydir=1)
    xyz_panels = interpolate_grid(xyz_edge, etas, interp_s, ydir=2)

    xyz_bound = interpolate_grid(xyz, eta_qtr, interp_c, ydir=1)
    xyz_corner = interpolate_grid(xyz_bound, etas, interp_s, ydir=2)
    xyz_center = interpolate_grid(xyz_bound, etabar, interp_s, ydir=2)

    xyz_cp = interpolate_grid(xyz, eta_thrqtr, interp_c, ydir=1)
    xyz_cp = interpolate_grid(xyz_cp, etabar, interp_s, ydir=2)

    # check which side we're working with
    right_side = sum(xyz[2,:,:]) > 0

    # initialize output
    panels = Matrix{SurfacePanel{TF}}(undef, nc, (1+mirror)*ns)

    # populate each panel
    for j = 1:ns

        # get chord for setting finite core size

        # also get chord length for setting finite core size
        c = norm(xyz_center[:,end,j] - xyz_center[:,1,j])

        for i = 1:nc
            # note that `i==1` corresponds to the leading edge for `xyz_corner` and `xyz_center`
            rtl = SVector(xyz_corner[1,i+1,j], xyz_corner[2,i+1,j], xyz_corner[3,i+1,j])
            rtc = SVector(xyz_center[1,i+1,j], xyz_center[2,i+1,j], xyz_center[3,i+1,j])
            rtr = SVector(xyz_corner[1,i+1,j+1], xyz_corner[2,i+1,j+1], xyz_corner[3,i+1,j+1])
            rbl = SVector(xyz_corner[1,i+2,j], xyz_corner[2,i+2,j], xyz_corner[3,i+2,j])
            rbc = SVector(xyz_center[1,i+2,j], xyz_center[2,i+2,j], xyz_center[3,i+2,j])
            rbr = SVector(xyz_corner[1,i+2,j+1], xyz_corner[2,i+2,j+1], xyz_corner[3,i+2,j+1])
            rcp = SVector(xyz_cp[1,i,j], xyz_cp[2,i,j], xyz_cp[3,i,j])
            ncp = cross(rcp - rtr, rcp - rtl)
            ncp /= norm(ncp)

            # set finite core size
            Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
            core_size = fcore(c, Δs)

            # set chord length of current panel
            r1 = xyz_panels[:, i, j] # top left
            r2 = xyz_panels[:, i, j+1] # top right
            r3 = xyz_panels[:, i+1, j] # bottom left
            r4 = xyz_panels[:, i+1, j+1] # bottom right
            chord = norm((r1 + r2)/2 - (r3 + r4)/2)

            ip = i
            jp = mirror*right_side*ns + j
            panels[ip, jp] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp,
                ncp, core_size, chord)
        end
    end

    # other side
    if mirror
        # first reflect grid
        if right_side
            xyz_l = reverse(xyz_panels, dims=3)
            xyz_l[2,:,:] .= -xyz_l[2,:,:]
            xyz_r = xyz_panels
        else
            xyz_l = xyz_panels
            xyz_r = reverse(xyz_panels, dims=3)
            xyz_r[2,:,:] .= -xyz_r[2,:,:]
        end
        xyz_panels = cat(xyz_l, xyz_r[:,:,2:end], dims=3)

        # then populate remaining panels
        for j = 1:ns
            for i = 1:nc
                ip = i
                jp1 = right_side*ns + j
                jp2 = 2*ns - right_side*ns - j + 1
                panels[ip,jp2] = reflect(panels[ip,jp1])
            end
        end
    end

    return xyz_panels, panels
end

"""
    wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
        fc = fill(x -> 0, length(xle)),
        mirror = false,
        fcore = (c, Δs) -> 1e-3,
        spacing_s = Cosine(),
        spacing_c = Uniform(),
        interp_s = (x, y, xpt) -> LinearInterpolation(x, y)(xpt))

Discretize a wing into `ns` spanwise and `nc` chordwise panels with associated
vortex rings according to the spanwise discretization scheme `spacing_s` and
chordwise discretization scheme `spacing_c`.

Return a grid with dimensions (3, i, j) containing the panel corners and a matrix
with dimensions (i, j) containing the generated panels.

# Arguments
 - `xle`: leading edge x-coordinate of each airfoil section
 - `yle`: leading edge y-coordinate of each airfoil section
 - `zle`: leading edge z-coordinate of each airfoil section
 - `chord`: chord length of each airfoil section
 - `theta`: twist of each airfoil section
 - `phi`: dihedral angle of each airfoil section, defined by a right hand rotation about the x-axis
 - `ns`: number of spanwise panels
 - `nc`: number of chordwise panels
 - `fc`: (optional) camber line function y=f(x) of each airfoil section
 - `mirror`:  mirror the geometry across the X-Z plane?, defaults to `false`
 - `fcore`: function for setting the finite core size based on the chord length
        (in the x-direction) and/or the panel width (in the y/z directions).
        Defaults to `(c, Δs) -> 1e-3`
 - `spacing_s`: spanwise discretization scheme, defaults to `Cosine()`
 - `spacing_c`: chordwise discretization scheme, defaults to `Uniform()`
 - `interp_s`: interpolation function between spanwise stations, defaults to linear interpolation
"""
function wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    fc = fill(x->0, length(xle)),
    mirror = false,
    fcore = (c, Δs) -> 1e-3,
    spacing_s = Cosine(),
    spacing_c = Uniform(),
    interp_s = (x, y, xpt) -> LinearInterpolation(x, y)(xpt))

    TF = promote_type(eltype(xle), eltype(yle), eltype(zle), eltype(chord), eltype(theta), eltype(phi))

    N = (1+mirror)*nc*ns

    # get spanwise and chordwise spacing
    etas, etabar = spanwise_spacing(ns+1, spacing_s)
    etac, eta_qtr, eta_thrqtr = chordwise_spacing(nc+1, spacing_c)

    # add leading and trailing edge to grid
    eta_qtr = vcat(0.0, eta_qtr, 1.0)

    # check input dimensions
    n = length(xle)
    @assert n == length(yle) == length(zle) == length(chord) == length(theta) == length(phi)

    # set bound vortex and control point chordwise locations
    xyz_edge = Array{TF, 3}(undef, 3, nc+1, n)
    xyz_bound = Array{TF, 3}(undef, 3, nc+2, n)
    xyz_cp = Array{TF, 3}(undef, 3, nc, n)
    for j = 1:n

        # rotation matrix for twist
        st, ct = sincos(theta[j])
        Rt = @SMatrix [ct 0 st; 0 1 0; -st 0 ct]

        # rotation matrix for dihedral
        sp, cp = sincos(phi[j])
        Rp = @SMatrix [1 0 0; 0 cp -sp; 0 sp cp]

        # location of leading edge
        rle = SVector(xle[j], yle[j], zle[j])

        # panel edge chordwise locations
        for i = 1:nc+1
            # bound vortex location
            xc = etac[i]
            zc = fc[j](xc)

            # location on airfoil
            r = SVector(xc, 0, zc)

            # scale by chord length
            r = chord[j] * r

            # apply twist
            r = Rt * r

            # apply dihedral
            r = Rp * r

            # add leading edge offset
            r = r + rle

            # store final location
            xyz_edge[:,i,j] = r
        end

        # bound vortex chordwise locations (and leading and trailing edges)
        for i = 1:nc+2
            # bound vortex location
            xc = eta_qtr[i]
            zc = fc[j](xc)

            # location on airfoil
            r = SVector(xc, 0, zc)

            # scale by chord length
            r = chord[j] * r

            # apply twist
            r = Rt * r

            # apply dihedral
            r = Rp * r

            # add leading edge offset
            r = r + rle

            # store final location
            xyz_bound[:,i,j] = r
        end

        # control point chordwise locations
        for i = 1:nc
            # bound vortex location
            xc = eta_thrqtr[i]
            zc = fc[j](xc)

            # location on airfoil
            r = SVector(xc, 0, zc)

            # scale by chord length
            r = chord[j] * r

            # apply twist
            r = Rt * r

            # apply dihedral
            r = Rp * r

            # add leading edge offset
            r = r + rle

            # store final location
            xyz_cp[:,i,j] = r
        end
    end

    # interpolate the grids to match the specified spanwise spacing
    xyz_panels = interpolate_grid(xyz_edge, etas, interp_s; ydir=2)
    xyz_corner = interpolate_grid(xyz_bound, etas, interp_s; ydir=2)
    xyz_center = interpolate_grid(xyz_bound, etabar, interp_s; ydir=2)
    xyz_cp = interpolate_grid(xyz_cp, etabar, interp_s; ydir=2)

    # initialize output
    panels = Matrix{SurfacePanel{TF}}(undef,  nc, (1+mirror)*ns)

    # check which side we're working with
    right_side = sum(yle) > 0

    # populate each panel
    for j = 1:ns

        # get chord for setting finite core size
        c = norm(xyz_center[:,end,j] - xyz_center[:,1,j])

        for i = 1:nc
            # note that `i==1` corresponds to the leading edge for `xyz_corner` and `xyz_center`
            rtl = SVector(xyz_corner[1,i+1,j], xyz_corner[2,i+1,j], xyz_corner[3,i+1,j])
            rtc = SVector(xyz_center[1,i+1,j], xyz_center[2,i+1,j], xyz_center[3,i+1,j])
            rtr = SVector(xyz_corner[1,i+1,j+1], xyz_corner[2,i+1,j+1], xyz_corner[3,i+1,j+1])
            rbl = SVector(xyz_corner[1,i+2,j], xyz_corner[2,i+2,j], xyz_corner[3,i+2,j])
            rbc = SVector(xyz_center[1,i+2,j], xyz_center[2,i+2,j], xyz_center[3,i+2,j])
            rbr = SVector(xyz_corner[1,i+2,j+1], xyz_corner[2,i+2,j+1], xyz_corner[3,i+2,j+1])
            rcp = SVector(xyz_cp[1,i,j], xyz_cp[2,i,j], xyz_cp[3,i,j])
            ncp = cross(rcp - rtr, rcp - rtl)
            ncp /= norm(ncp)

            # set finite core size
            Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
            core_size = fcore(c, Δs)

            # set chord length of current panel
            r1 = xyz_panels[:,i,j] # top left
            r2 = xyz_panels[:,i,j+1] # top right
            r3 = xyz_panels[:,i+1, j] # bottom left
            r4 = xyz_panels[:,i+1, j+1] # bottom right
            chord = norm((r1 + r2)/2 - (r3 + r4)/2)

            ip = i
            jp = mirror*right_side*ns + j
            panels[ip, jp] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp,
                ncp, core_size, chord)
        end
    end

    # other side
    if mirror
        # first reflect grid
        if right_side
            xyz_l = reverse(xyz_panels, dims=3)
            xyz_l[2,:,:] .= -xyz_l[2,:,:]
            xyz_r = xyz_panels
        else
            xyz_l = xyz_panels
            xyz_r = reverse(xyz_panels, dims=3)
            xyz_r[2,:,:] .= -xyz_r[2,:,:]
        end
        xyz_panels = cat(xyz_l, xyz_r[:,:,2:end], dims=3)

        # then populate remaining panels
        for j = 1:ns
            for i = 1:nc
                ip = i
                jp1 = right_side*ns + j
                jp2 = 2*ns - right_side*ns - j + 1
                panels[ip,jp2] = reflect(panels[ip,jp1])
            end
        end
    end

    return xyz_panels, panels
end

"""
    update_surface_panels!(surface, grid; fcore = (c, Δs) -> 1e-3)

Updates the surface panels in `surface` to correspond to the grid coordinates in
`grid`.
"""
function update_surface_panels!(surface, grid; fcore = (c, Δs) -> 1e-3)

    TF = eltype(eltype(surface))

    nc, ns = size(surface) # number of chordwise and spanwise panels

    # populate each panel
    for j = 1:ns

        # get leading edge panel corners
        r1n = SVector(grid[1,1,j], grid[2,1,j], grid[3,1,j]) # top left
        r2n = SVector(grid[1,1,j+1], grid[2,1,j+1], grid[3,1,j+1]) # top right
        r3n = SVector(grid[1,2,j], grid[2,2,j], grid[3,2,j]) # bottom left
        r4n = SVector(grid[1,2,j+1], grid[2,2,j+1], grid[3,2,j+1]) # bottom right

        # also get chord length for setting finite core size
        cl = norm(grid[:,end,j] - grid[:,1,j])
        cr = norm(grid[:,end,j+1] - grid[:,1,j+1])
        c = (cl + cr)/2

        for i = 1:nc-1

            # grid corners of current panel
            r1 = r1n # top left of panel
            r2 = r2n # top right of panel
            r3 = r3n # bottom left of panel
            r4 = r4n # bottom right of panel

            # grid corners of next panel
            r1n = r3 # top left
            r2n = r4 # top right
            r3n = SVector(grid[1,i+2,j], grid[2,i+2,j], grid[3,i+2,j]) # bottom left
            r4n = SVector(grid[1,i+2,j+1], grid[2,i+2,j+1], grid[3,i+2,j+1]) # bottom right

            # top left corner of ring vortex
            rtl = linearinterp(0.25, r1, r3)

            # top right corner of ring vortex
            rtr = linearinterp(0.25, r2, r4)

            # center of top bound vortex
            rtc = (rtl + rtr)/2

            # bottom left corner of ring vortex
            rbl = linearinterp(0.25, r1n, r3n)

            # bottom right corner of ring vortex
            rbr = linearinterp(0.25, r2n, r4n)

            # center of bottom bound vortex
            rbc = (rbl + rbr)/2

            # control point
            rtop = linearinterp(0.5, r1, r2)
            rbot = linearinterp(0.5, r3, r4)
            rcp = linearinterp(0.75, rtop, rbot)

            # surface normal
            ncp = cross(rcp - rtr, rcp - rtl)
            ncp /= norm(ncp)

            # set finite core size
            Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
            core_size = fcore(c, Δs)

            # get chord length of current panel
            chord = norm((r1 + r2)/2 - (r3 + r4)/2)

            surface[i,j] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp,
                ncp, core_size, chord)
        end

        # grid corners of current panel
        r1 = r1n # top left
        r2 = r2n # top right
        r3 = r3n # bottom left
        r4 = r4n # bottom right

        # top left corner of ring vortex
        rtl = linearinterp(0.25, r1, r3)

        # top right corner of ring vortex
        rtr = linearinterp(0.25, r2, r4)

        # center of top bound vortex
        rtc = (rtl + rtr)/2

        # bottom left corner of ring vortex
        rbl = r3

        # bottom right corner of ring vortex
        rbr = r4

        # center of bottom bound vortex
        rbc = (rbl + rbr)/2

        # control point
        rtop = linearinterp(0.5, r1, r2)
        rbot = linearinterp(0.5, r3, r4)
        rcp = linearinterp(0.75, rtop, rbot)

        # surface normal
        ncp = cross(rcp - rtr, rcp - rtl)
        ncp /= norm(ncp)

        # set finite core size
        Δs = sqrt((rtr[2]-rtl[2])^2 + (rtr[3]-rtl[3])^2)
        core_size = fcore(c, Δs)

        # get chord length of current panel
        chord = norm((r1 + r2)/2 - (r3 + r4)/2)

        surface[end, j] = SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp,
            core_size, chord)
    end

    return surface
end

"""
    lifting_line_geometry(grids, xc=0.25)
Construct a lifting line representation of the surfaces in `grids` at the
normalized chordwise location `xc`.  Return the lifting line coordinates
and chord lengths.
# Arguments
 - `grids`: Vector with length equal to the number of surfaces.  Each element of
    the vector is a grid with shape (3, nc, ns) which defines the discretization
    of a surface into panels. `nc` is the number of chordwise panels and `ns` is
    the number of spanwise panels.
 - `xc`: Normalized chordwise location of the lifting line from the leading edge.
    Defaults to the quarter chord
# Return Arguments:
 - `r`: Vector with length equal to the number of surfaces, with each element
    being a matrix with size (3, ns+1) which contains the x, y, and z coordinates
    of the resulting lifting line coordinates
 - `c`: Vector with length equal to the number of surfaces, with each element
    being a vector of length `ns+1` which contains the chord lengths at each
    lifting line coordinate.
"""
function lifting_line_geometry(grids, xc=0.25)
    TF = eltype(eltype(grids))
    nsurf = length(grids)
    r = Vector{Matrix{TF}}(undef, nsurf)
    c = Vector{Vector{TF}}(undef, nsurf)
    for isurf = 1:nsurf
        ns = size(grids[isurf], 3) - 1
        r[isurf] = Matrix{TF}(undef, 3, ns+1)
        c[isurf] = Vector{TF}(undef, ns+1)
    end
    return lifting_line_geometry!(r, c, grids, xc)
end

"""
    lifting_line_geometry!(r, c, grids, xc=0.25)
In-place version of [`lifting_line_geometry`](@ref)
"""
function lifting_line_geometry!(r, c, grids, xc=0.25)
    nsurf = length(grids)
    # iterate through each lifting surface
    for isurf = 1:nsurf
        # extract current grid
        grid = grids[isurf]
        # dimensions of this grid
        nc = size(grid, 2)
        ns = size(grid, 3)
        # loop through each spanwise section
        for j = 1:ns
            # get leading and trailing edges
            le = SVector(grid[1,1,j], grid[2,1,j], grid[3,1,j])
            te = SVector(grid[1,end,j], grid[2,end,j], grid[3,end,j])
            # get quarter-chord
            r[isurf][:,j] = linearinterp(xc, le, te)
            # get chord length
            c[isurf][j] = norm(le - te)
        end
    end
    return r, c
end

"""
    translate(grid, r)

Return a copy of the grid points in `grid` translated the distance specified by vector `r`
"""
translate(grid, r) = grid .+ r

"""
    translate!(grid, r)

Translate the grid points in `grid` the distance specified by vector `r`
"""
function translate!(grid, r)
    grid .+= r
    return grid
end

"""
    rotate(grid, R, r = [0,0,0])

Return a copy of the grid points in `grid` rotated about point `r` using the rotation matrix `R`
"""
rotate(grid, args...) = rotate!(copy(grid), args...)

"""
    rotate!(grid, R, r = [0,0,0])

Rotate the grid points in `grid` about point `r` using the rotation matrix `R`
"""
function rotate!(grid, R, r = (@SVector zeros(3)))
    nc, ns = size(grid, 2), size(grid, 3)

    grid .-= r

    for i = 1:nc, j = 1:ns
        grid[:,i,j] = R*grid[:,i,j]
    end

    grid .+= r

    return grid
end

"""
    trailing_edge_points(surface[s])

Return points on the trailing edge of each surface.
"""
trailing_edge_points(surfaces) = [bottom_left.(surfaces[end, :])..., bottom_right(surfaces[end, end])]

"""
    repeated_trailing_edge_points(surface[s])

Generates a dictionary of the form `Dict((isurf, i) => [(jsurf1, j1), (jsurf2, j2)...]` which
defines repeated trailing edge points.  Trailing edge point `i` on surface
`isurf` is repeated on surface `jsurf1` at point `j1`, `jsurf2` at point `j2`,
and so forth.
"""
repeated_trailing_edge_points

repeated_trailing_edge_points(surface::AbstractMatrix) = repeated_trailing_edge_points([surface])

function repeated_trailing_edge_points(surfaces::AbstractVector{<:AbstractMatrix})

    nsurf = length(surfaces)
    points = trailing_edge_points.(surfaces)
    repeated_points = Dict{NTuple{2, Int}, Vector{NTuple{2, Int}}}()

    # loop through all points
    for isurf = 1:nsurf, i = 1:length(points[isurf])
        # start with empty vector
        repeats = NTuple{2, Int}[]
        # add all repeated points to the vector
        for jsurf = 1:nsurf, j = 1:length(points[jsurf])
            if (isurf, i) != (jsurf, j) && isapprox(points[isurf][i], points[jsurf][j])
                push!(repeats, (jsurf, j))
            end
        end
        # only store vector if nonempty
        if !isempty(repeats)
            repeated_points[(isurf, i)] = repeats
        end
    end

    return repeated_points
end

"""
    flipy(r)

Flip sign of y-component of vector `r` (used for symmetry)
"""
@inline flipy(r) = SVector{3}(r[1], -r[2], r[3])

"""
    on_symmetry_plane(args...; tol=eps())

Test whether points in `args` are on the symmetry plane (y = 0)
"""
on_symmetry_plane

on_symmetry_plane(r; tol=eps()) = isapprox(r[2], 0.0, atol=tol)

@inline function on_symmetry_plane(r, args...; tol=eps())
    return on_symmetry_plane(r; tol=tol) && on_symmetry_plane(args...; tol=tol)
end

"""
    not_on_symmetry_plane(args...; tol=eps())

Test whether and of the points in `args` are not on the symmetry plane (y = 0)
"""
not_on_symmetry_plane(args...; tol=eps()) = !on_symmetry_plane(args...; tol=tol)


text:
module TestMultiReturningFunction

using Modia
@usingModiaPlot

include("$(Modia.path)/models/Blocks.jl")
include("$(Modia.path)/models/Electric.jl")
include("$(Modia.path)/models/Rotational.jl")

function twoCoupledInertias(J1,J2,derw1,derw2,tau0,tau1,tau2)::Array{Float64,1}
    r1 = J1*derw1 - tau0 + tau1
    r2 = J2*derw2 - tau1 + tau2
    return [r1,r2]
end

ThreeCoupledInertias = Model(
    J1 = 1.1,
    J2 = 1.2,
    J3 = 1.3,
    phi1 = Var(init=0.0),
    phi2 = Var(init=0.0),
    phi3 = Var(init=0.0),
    w1   = Var(init=0.0),
    w2   = Var(init=0.0),
    w3   = Var(init=0.0),
    tau1 = Var(start=0.0),
    tau2 = Var(start=0.0),

    # Same model with components from the Rotational.jl library
    inertia1 = Inertia | Map(J=:J1),
    inertia2 = Inertia | Map(J=:J2),
    inertia3 = Inertia | Map(J=:J3),
    connect = :[
        (inertia1.flange_b, inertia2.flange_a)
        (inertia2.flange_b, inertia3.flange_a)
    ],

    equations = :[
        tau0 = sin(time/u"s")
        w1 = der(phi1)
        w2 = der(phi2)
        w3 = der(phi3)
        phi2 = phi1
        phi3 = phi2
        (0,0) = twoCoupledInertias(J1,J2,der(w1),der(w2),tau0,tau1,tau2)
        J3*der(w3) = tau2
        inertia1.flange_a.tau = tau0
        inertia3.flange_b.tau = 0.0
    ]
)

threeCoupledInertias = @instantiateModel(ThreeCoupledInertias, unitless=true, log=false, logDetails=false, logCode=true, logStateSelection=false)

simulate!(threeCoupledInertias, stopTime = 2.0, log=true,
          requiredFinalStates=[0.3029735305821086, 0.3933746028781303, 0.39337460287813036, 0.3029735305821086])

plot(threeCoupledInertias, [("phi1", "phi2", "phi3", "inertia1.phi"),
                            ("w1","w2","w3", "inertia1.w"),
                            ("tau0"),
                            ("tau1", "inertia2.flange_a.tau"),
                            ("tau2", "inertia3.flange_a.tau")])

end


text:
instr = readlines("input.txt")
reg = parse.(unique(first.(split.(instr))))
cmd = Dict("inc" => "+=", "dec" => "-=")
instr = replace.(instr, r"(inc|dec)", x -> cmd[x])
for r in reg
    @eval $r = 0
end

high = 0
for i in instr
    cmd, cond = map(parse, split(i, "if "))
    @eval $cond ? $cmd : nothing
    high = max(high, maximum(eval.(reg)))
end
println(maximum(eval.(reg)))
println(high)


text:
module LifeValStructures

using DataFrames

export qxTables, wxTables, vtTables, expTables, chTables, fTables

export Policy, LifeAssurance, LifeAssuranceUL, TermAssurance, Annuity, EndowmentAssurance, EndowmentAssuranceUL

const MAX_AGE = 108

abstract type typeBasis end

struct qxTables <: typeBasis
    #id::string
    runno::Int64
    spcode::Int64
    datfram::DataFrame
    agemax_m::Int64
    agemax_f::Int64
    function qxTables(runno, spcode, datfram)
        #
        agemax_m=117
        for n in eachrow(datfram)
            if  n.QXM == 1
                agemax_m = n.AGE - 1
                break;
            end
        end
        #
        agemax_f=109
        for n in eachrow(datfram)
            if  n.QXF == 1
                agemax_f = n.AGE - 1
                break;
            end
        end
        return(new(runno, spcode, datfram, agemax_m, agemax_f))
    end
end

struct wxTables <: typeBasis
    #id::string
    runno::Int64
    spcode::Int64
    datfram::DataFrame
    function wxTables(runno, spcode, datfram)
        return(new(runno, spcode, datfram))
    end
end

struct vtTables <: typeBasis
    #id::string
    runno::Int64
    spcode::Int64
    datfram::DataFrame
    function vtTables(runno, spcode, datfram)
        return(new(runno, spcode, datfram))
    end
end

struct expTables <: typeBasis
    #id::string
    runno::Int64
    spcode::Int64
    datfram::DataFrame  #this should be a vector!
    function expTables(runno, spcode, datfram)
        return(new(runno, spcode, datfram))
    end
end

struct chTables <: typeBasis
    #id::string
    runno::Int64
    spcode::Int64
    datfram::DataFrame  #this should be a vector!
    function chTables(runno, spcode, datfram)
        return(new(runno, spcode, datfram))
    end
end

struct fTables <: typeBasis
    #id::string
    runno::Int64
    fundid::Int64
    datfram::DataFrame #this should be a vector!
    function fTables(runno, fundid, datfram)
        return(new(runno, fundid, datfram))
    end
end

# create Policy type
abstract type Policy end

# Life Assurance
# ---------------------------------------------------------------------------- #
struct LifeAssurance <: Policy
    id::Int64
    gender::Bool
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    units_1::Float64
    units_99::Float64
    alloc_1::Float64
    alloc_99::Float64
    projMax::Int64

    function LifeAssurance(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99)
        # check for negative premium
        if premium < 0
            error("NegativePremium Error")
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = 12(MAX_AGE-age)
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99, projMax))
    end
end

struct LifeAssuranceUL <: Policy
    id::Int64
    gender::Bool
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    units_1::Float64
    units_99::Float64
    alloc_1::Float64
    alloc_99::Float64
    projMax::Int64

    function LifeAssuranceUL(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99)
        # check for negative premium
        if premium < 0
            error("NegativePremium Error")
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = 12(MAX_AGE-age)
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99, projMax))
    end
end

struct Annuity <: Policy
    id::Int64
    gender::Bool
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    units_1::Float64
    units_99::Float64
    alloc_1::Float64
    alloc_99::Float64
    projMax::Int64

    function Annuity(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99)
        # check for negative premium
        if premium != 0
            #error("Premium Error")
            #just set it to 0
            premium = 0
        end

        if prem_esc != 0
            #error("Premium Error")
            #just set it to 0
            prem_esc = 0
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = 12(MAX_AGE-age)
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99, projMax))
    end
end

struct TermAssurance <: Policy
    id::Int64
    gender::Bool
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    units_1::Float64
    units_99::Float64
    alloc_1::Float64
    alloc_99::Float64
    projMax::Int64

    function TermAssurance(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99)
        # check for negative premium
        if premium < 0
            error("NegativePremium Error")
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = term-termIF
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99, projMax))
    end
end

struct EndowmentAssurance <: Policy
    id::Int64
    gender::Bool
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    units_1::Float64
    units_99::Float64
    alloc_1::Float64
    alloc_99::Float64
    projMax::Int64

    function EndowmentAssurance(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99)
        # check for negative premium
        if premium < 0
            error("NegativePremium Error")
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = term-termIF
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, alloc_1, alloc_99, projMax))
    end
end

struct EndowmentAssuranceUL <: Policy
    id::Int64
    gender::Int64
    age::Int64
    term::Int64
    termIF::Int64
    premium::Float64
    SumAssured::Float64
    spc::Int64
    prem_esc::Float64
    sa_esc::Float64
    esc_m::Int64
    #units::Vector
    #allocPrem::Vector
    units_1::Float64
    units_99::Float64
    allocPrem_1::Float64
    allocPrem_99::Float64
    projMax::Int64

    function EndowmentAssuranceUL(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, allocPrem_1, allocPrem_99)
        # check for negative premium
        if premium < 0
            error("NegativePremium Error")
        end

        if age>MAX_AGE
            error("MaxAge Error")
        end

        global MAX_AGE
        projMax = term-termIF
        return(new(id, gender, age, term, termIF, premium, SumAssured, spc, prem_esc, sa_esc, esc_m, units_1, units_99, allocPrem_1, allocPrem_99, projMax))
    end
end

end # module


text:
using CSTParser
using Test

import CSTParser: parse, remlineinfo!, span, flisp_parse, typof, kindof, valof

include("parser.jl")
include("interface.jl")
CSTParser.check_base()


text:
export ⊼, atomize

"""
    a ⊼ b

Not AND operator (NAND for short). The NAND operator is a functionally complete boolean operator, meaning all other boolean
operations can be expressed with only NAND operations.
"""
⊼ = LogicalOperation((x, y) -> !(x && y), :⊼, 2, true, true)

"""
    atomize(expr::AbstractExpression)

Converts an expression into its directly converted NAND-only form.
"""
# Reduce an expression to its atomic form using the NAND operation
function atomize(expr::LogicalExpression)
    parts = atomize.(arguments(expr))
    expr_op = operation(expr)

    if expr_op == ¬
        return parts[1] ⊼ parts[1]
    elseif expr_op == ∧
        return atomize(¬(parts[1] ⊼ parts[2]))
    elseif expr_op == ∨
        return atomize(¬parts[1] ⊼ ¬parts[2])
    elseif expr_op == →
        return atomize(¬parts[1] ∨ parts[2])
    elseif expr_op == ⟷
        return atomize((parts[1] ∧ parts[2]) ∨ (¬parts[1] ∧ ¬parts[2]))
    elseif expr_op == ⊼
        return atomize(parts[1]) ⊼ atomize(parts[2])
    end

    throw(ErrorException("Expression contains unsupported atomic operator $(expr_op)"))
end
atomize(s::LogicalSymbol) = s


text:
function gen_dist(scl=5, nvars=2, nmix=nv)
    Product([MixtureModel([Normal(scl*rand(), scl*rand()) for _ in 1:nmix]) for _ in 1:nvars])
end


text:
# Helpful utility functions

# Set the value of a field of a pointer
# Equivalent to s->name = value
function av_setfield{T}(s::Ptr{T}, name::Symbol, value)
    field = findfirst(fieldnames(T), name)
    byteoffset = fieldoffset(T, field)
    S = T.types[field]

    p = convert(Ptr{S}, s+byteoffset)
    a = unsafe_wrap(Array, p,1)
    a[1] = convert(S, value)
end
export av_setfield
function av_pointer_to_field{T}(s::Ptr{T}, name::Symbol)
    field = findfirst(fieldnames(T), name)
    byteoffset = fieldoffset(T, field)
    return s + byteoffset
end
export av_pointer_to_field
av_pointer_to_field(s::Array, name::Symbol) = av_pointer_to_field(pointer(s), name)

function open_stdout_stderr(cmd::Cmd)
    out = Base.PipeEndpoint()
    err = Base.PipeEndpoint()
    cmd_out = Base.PipeEndpoint()
    cmd_err = Base.PipeEndpoint()
    Base.link_pipe(out, true, cmd_out, false)
    Base.link_pipe(err, true, cmd_err, false)

    r = spawn(ignorestatus(cmd), (DevNull, cmd_out, cmd_err))

    Base.close_pipe_sync(cmd_out)
    Base.close_pipe_sync(cmd_err)

    # NOTE: these are not necessary on v0.4 (although they don't seem
    #       to hurt). Remove when we drop support for v0.3.
    Base.start_reading(out)
    Base.start_reading(err)
    return (out, err, r)
end

function readall_stdout_stderr(cmd::Cmd)
    (out, err, proc) = open_stdout_stderr(cmd)
    return (readstring(out), readstring(err))
end
export readall_stdout_stderr

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module TestJuliaR
using Test

include("test_rmachine.jl")

end

max_stars_repo_name: ppalmes/AMLPipelineR

max_issues_repo_name: ppalmes/AMLPipelineR

max_forks_repo_name: ppalmes/AMLPipelineR


text:
############################################################################################
# Import External Packages
using Test
using Random: Random, AbstractRNG, seed!

#using Soss
using LinearAlgebra
using Distributions, Bijectors, DistributionsAD
using ForwardDiff, ReverseDiff, Zygote

############################################################################################
# Import Baytes Packages
using ModelWrappers
using ModelWrappers:
    ModelWrappers,
    FlattenDefault,
    FlattenContinuous,
    FlattenAll,
    UnflattenStrict,
    UnflattenFlexible,
    _get_constraint,
    _checkparams,
    _checkfinite,
    _checkprior,
    _allparam,
    _anyparam,
    _checksampleable,
    construct_flatten,
    construct_transform,
    flatten,
    flattenAD,
    unflatten,
    unflattenAD,
    constrain,
    _get_val,
    _get_constraint,
    log_density,
    _log_density,
    log_density_and_gradient,
    _log_density_and_gradient,
    tag,
    _check,
    flatten_Symmetric,
    Symmetric_from_flatten,
    flatten_Simplex,
    Simplex_from_flatten!,
    Simplex_from_flatten

############################################################################################
# Include Files
include("TestHelper.jl")

############################################################################################
# Run Tests
@testset "All tests" begin

    include("test-flatten/flatten.jl")

    include("test-core.jl")
    include("test-flatten.jl")

    include("test-models.jl")
    include("test-tagged.jl")
    include("test-objective.jl")
    include("test-differentiation.jl")

end


text:
@testset "$TEST $G" begin

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)

y = closeness_centrality(g5; normalize=false)
z = closeness_centrality(g5)

@test y == [0.75, 0.6666666666666666, 1.0, 0.0]
@test z == [0.75, 0.4444444444444444, 0.3333333333333333, 0.0]

g = G(5)
add_edge!(g,1,2)
z = closeness_centrality(g)
@test z[1] == z[2] == 0.25
@test z[3] == z[4] == z[5] == 0.0

end # testset


text:
module SearchLightPostgreSQL

import LibPQ, DataFrames, Logging
import SearchLight

#
# Setup
#

const DEFAULT_PORT = 5432

const COLUMN_NAME_FIELD_NAME = :column_name

function SearchLight.column_field_name()
  COLUMN_NAME_FIELD_NAME
end

const DatabaseHandle = LibPQ.Connection
const ResultHandle   = LibPQ.Result

const TYPE_MAPPINGS = Dict{Symbol,Symbol}( # Julia / Postgres
  :char       => :CHARACTER,
  :string     => :VARCHAR,
  :text       => :TEXT,
  :integer    => :INTEGER,
  :int        => :INTEGER,
  :bigint     => :BIGINT,
  :serial     => :SERIAL,
  :bigserial  => :BIGSERIAL,
  :int4range  => :INT4RANGE,
  :int8range  => :INT8RANGE,
  :float      => :FLOAT,
  :decimal    => :DECIMAL,
  :datetime   => :INTEGER,
  :timestamp  => :TIMESTAMP,
  :timestamptz => :TIMESTAMPTZ,
  :tstzrange  => :TSTZRANGE,
  :time       => :TIME,
  :date       => :DATE,
  :binary     => :BLOB,
  :boolean    => :BOOLEAN,
  :bool       => :BOOLEAN
)

const CONNECTIONS = DatabaseHandle[]

#
# Connection
#

"""
    connect(conn_data::Dict)::DatabaseHandle

Connects to the database and returns a handle.
"""
function SearchLight.connect(conn_data::Dict = SearchLight.config.db_config_settings) :: DatabaseHandle
  dns = String[]

  haskey(conn_data, "host")     && push!(dns, string("host=", conn_data["host"]))
  haskey(conn_data, "hostaddr") && push!(dns, string("hostaddr=", conn_data["hostaddr"]))
  haskey(conn_data, "port")     && push!(dns, string("port=", conn_data["port"]))
  haskey(conn_data, "database") && push!(dns, string("dbname=", conn_data["database"]))
  haskey(conn_data, "username") && push!(dns, string("user=", conn_data["username"]))
  haskey(conn_data, "password") && push!(dns, string("password=", conn_data["password"]))
  haskey(conn_data, "passfile") && push!(dns, string("passfile=", conn_data["passfile"]))
  haskey(conn_data, "connect_timeout") && push!(dns, string("connect_timeout=", conn_data["connect_timeout"]))
  haskey(conn_data, "client_encoding") && push!(dns, string("client_encoding=", conn_data["client_encoding"]))

  push!(CONNECTIONS, LibPQ.Connection(join(dns, " ")))[end]
end

"""
    disconnect(conn::DatabaseHandle)::Nothing

Disconnects from database.
"""
function SearchLight.disconnect(conn::DatabaseHandle = SearchLight.connection()) :: Nothing
  LibPQ.close(conn)
end

function SearchLight.connection()
  isempty(CONNECTIONS) && throw(SearchLight.Exceptions.NotConnectedException())

  CONNECTIONS[end]
end

#
# Data sanitization
#

"""
    escape_column_name(c::String, conn::DatabaseHandle)::String

Escapes the column name.

# Examples
```julia
julia>
```
"""
function SearchLight.escape_column_name(c::String, conn::DatabaseHandle = SearchLight.connection()) :: String
  join(["""\"$(replace(cx, "\""=>"'"))\"""" for cx in split(c, '.')], '.')
end

"""
    escape_value{T}(v::T, conn::DatabaseHandle)::T

Escapes the value `v` using native features provided by the database backend if available.

# Examples
```julia
julia>
```
"""
function SearchLight.escape_value(v::T, _ = nothing)::T where {T} # TODO: deprecate in next major version
  isa(v, Number) ? v : "E'$(replace(string(v), "'"=>"\\'"))'"
end

#
# Query execution
#

function SearchLight.query(sql::String, conn::DatabaseHandle = SearchLight.connection()) :: DataFrames.DataFrame
  @info sql
  result = LibPQ.execute(conn, sql)

  if LibPQ.error_message(result) != ""
    throw(SearchLight.Exceptions.DatabaseAdapterException("$(string(LibPQ)) error: $(LibPQ.errstring(result)) [$(LibPQ.errcode(result))]"))
  end

  result |> DataFrames.DataFrame
end

function SearchLight.to_find_sql(m::Type{T}, q::SearchLight.SQLQuery, joins::Union{Nothing,Vector{SearchLight.SQLJoin}} = nothing)::String where {T<:SearchLight.AbstractModel}
  sql::String = ( string("$(SearchLight.to_select_part(m, q.columns, joins)) $(SearchLight.to_from_part(m)) $(SearchLight.to_join_part(m, joins)) $(SearchLight.to_where_part(q.where)) ",
                      "$(SearchLight.to_group_part(q.group)) $(SearchLight.to_having_part(q.having)) $(SearchLight.to_order_part(m, q.order)) ",
                      "$(SearchLight.to_limit_part(q.limit)) $(SearchLight.to_offset_part(q.offset))")) |> strip
  replace(sql, r"\s+"=>" ")
end

function SearchLight.to_store_sql(m::T; conflict_strategy = :error)::String where {T<:SearchLight.AbstractModel}
  uf = SearchLight.persistable_fields(typeof(m))

  sql = if ! SearchLight.ispersisted(m) || (SearchLight.ispersisted(m) && conflict_strategy == :update)
    pos = findfirst(x -> x == SearchLight.primary_key_name(m), uf)
    pos > 0 && splice!(uf, pos)

    fields = SearchLight.SQLColumn(uf)
    vals = join( map(x -> string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))), uf), ", ")

    "INSERT INTO $(SearchLight.table(typeof(m))) ( $fields ) VALUES ( $vals )" *
        if ( conflict_strategy == :error ) ""
        elseif ( conflict_strategy == :ignore ) " ON CONFLICT DO NOTHING"
        elseif ( conflict_strategy == :update &&
          getfield(m, Symbol(SearchLight.primary_key_name(m))).value !== nothing )
            " ON CONFLICT ($(SearchLight.primary_key_name(m))) DO UPDATE SET $(SearchLight.update_query_part(m))"
        else ""
        end
  else
    "UPDATE $(SearchLight.table(typeof(m))) SET $(SearchLight.update_query_part(m))"
  end

  return string(sql, " RETURNING $(SearchLight.primary_key_name(m))")
end

function SearchLight.delete_all(m::Type{T}; truncate::Bool = true, reset_sequence::Bool = true, cascade::Bool = false)::Nothing where {T<:SearchLight.AbstractModel}
  if truncate
    sql = "TRUNCATE $(SearchLight.table(m))"
    reset_sequence ? sql *= " RESTART IDENTITY" : ""
    cascade ? sql *= " CASCADE" : ""
  else
    sql = "DELETE FROM $(SearchLight.table(m))"
  end

  SearchLight.query(sql)

  nothing
end

function SearchLight.delete(m::T)::T where {T<:SearchLight.AbstractModel}
  SearchLight.ispersisted(m) || throw(SearchLight.Exceptions.NotPersistedException(m))

  "DELETE FROM $(SearchLight.table(typeof(m))) WHERE $(SearchLight.primary_key_name(typeof(m))) = '$(m.id.value)'" |> SearchLight.query

  m.id = SearchLight.DbId()

  m
end

function Base.count(m::Type{T}, q::SearchLight.SQLQuery = SearchLight.SQLQuery())::Int where {T<:SearchLight.AbstractModel}
  count_column = SearchLight.SQLColumn("COUNT(*) AS __cid", raw = true)
  q = SearchLight.clone(q, :columns, push!(q.columns, count_column))

  SearchLight.DataFrame(m, q)[1, Symbol("__cid")]
end

function SearchLight.update_query_part(m::T)::String where {T<:SearchLight.AbstractModel}
  update_values = join(map(x -> "$(string(SearchLight.SQLColumn(x))) = $(string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))) )", SearchLight.persistable_fields(typeof(m))), ", ")

  " $update_values WHERE $(SearchLight.table(typeof(m))).$(SearchLight.primary_key_name(typeof(m))) = '$(m.id.value)'"
end

function SearchLight.column_data_to_column_name(column::SearchLight.SQLColumn, column_data::Dict{Symbol,Any}) :: String
  "$(SearchLight.to_fully_qualified(column_data[:column_name], column_data[:table_name])) AS $(isempty(column_data[:alias]) ? SearchLight.to_sql_column_name(column_data[:column_name], column_data[:table_name]) : column_data[:alias] )"
end

function SearchLight.to_from_part(m::Type{T})::String where {T<:SearchLight.AbstractModel}
  "FROM " * SearchLight.escape_column_name(SearchLight.table(m), SearchLight.connection())
end

function SearchLight.to_where_part(w::Vector{SearchLight.SQLWhereEntity})::String
  where = isempty(w) ?
          "" :
          string("WHERE ",
                (string(first(w).condition) == "AND" ? "TRUE " : "FALSE "),
                join(map(wx -> string(wx), w), " "))

  replace(where, r"WHERE TRUE AND "i => "WHERE ")
end

function SearchLight.to_order_part(m::Type{T}, o::Vector{SearchLight.SQLOrder})::String where {T<:SearchLight.AbstractModel}
  isempty(o) ?
    "" :
    string("ORDER BY ",
            join(map(x -> string((! SearchLight.is_fully_qualified(x.column.value) ?
                                    SearchLight.to_fully_qualified(m, x.column) :
                                    x.column.value), " ", x.direction),
                      o), ", "))
end

function SearchLight.to_group_part(g::Vector{SearchLight.SQLColumn}) :: String
  isempty(g) ?
    "" :
    string(" GROUP BY ", join(map(x -> string(x), g), ", "))
end

function SearchLight.to_limit_part(l::SearchLight.SQLLimit) :: String
  l.value != "ALL" ? string("LIMIT ", string(l)) : ""
end

function SearchLight.to_offset_part(o::Int) :: String
  o != 0 ? string("OFFSET ", string(o)) : ""
end

function SearchLight.to_having_part(h::Vector{SearchLight.SQLWhereEntity}) :: String
  having =  isempty(h) ?
            "" :
            string("HAVING ",
                  (string(first(h).condition) == "AND" ? "TRUE " : "FALSE "),
                  join(map(w -> string(w), h), " "))

  replace(having, r"HAVING TRUE AND "i => "HAVING ")
end

function SearchLight.to_join_part(m::Type{T}, joins::Union{Nothing,Vector{SearchLight.SQLJoin}} = nothing)::String where {T<:SearchLight.AbstractModel}
  joins === nothing && return ""

  join(map(x -> string(x), joins), " ")
end

function Base.rand(m::Type{T}; limit = 1)::Vector{T} where {T<:SearchLight.AbstractModel}
  SearchLight.find(m, SearchLight.SQLQuery(limit = SearchLight.SQLLimit(limit), order = [SearchLight.SQLOrder("random()", raw = true)]))
end

#### MIGRATIONS ####

"""
    create_migrations_table(table_name::String)::Nothing

Runs a SQL DB query that creates the table `table_name` with the structure needed to be used as the DB migrations table.
The table should contain one column, `version`, unique, as a string of maximum 30 chars long.
"""
function SearchLight.Migration.create_migrations_table(table_name::String = SearchLight.config.db_migrations_table_name) :: Nothing
  SearchLight.query("CREATE TABLE $table_name (version varchar(30))")

  @info "Created table $table_name"

  nothing
end

function SearchLight.Migration.create_table(f::Function, name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: Nothing
  SearchLight.query(create_table_sql(f, string(name), options))

  nothing
end

function create_table_sql(f::Function, name::String, options::String = "") :: String
  "CREATE TABLE $name (" * join(f()::Vector{String}, ", ") * ") $options" |> strip
end

function SearchLight.Migration.column(name::Union{String,Symbol}, column_type::Union{String,Symbol}, options::Any = ""; default::Any = nothing, limit::Union{Int,Nothing,String} = nothing, not_null::Bool = false) :: String
  "$name $(TYPE_MAPPINGS[column_type] |> string) " *
    (isa(limit, Int) ? "($limit)" : "") *
    (default === nothing ? "" : " DEFAULT $default ") *
    (not_null ? " NOT NULL " : "") *
    string(options)
end

function SearchLight.Migration.column_id(name::Union{String,Symbol} = "id", options::Union{String,Symbol} = ""; constraint::Union{String,Symbol} = "", nextval::Union{String,Symbol} = "") :: String
  "$name SERIAL $constraint PRIMARY KEY $nextval $options"
end

function SearchLight.Migration.add_index(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}; name::Union{String,Symbol} = "", unique::Bool = false, order::Union{String,Symbol} = :none) :: Nothing
  name = isempty(name) ? SearchLight.index_name(table_name, column_name) : name
  SearchLight.query("CREATE $(unique ? "UNIQUE" : "") INDEX $(name) ON $table_name ($column_name)")

  nothing
end

function SearchLight.Migration.add_column(table_name::Union{String,Symbol}, name::Union{String,Symbol}, column_type::Union{String,Symbol}; default::Union{String,Symbol,Nothing} = nothing, limit::Union{Int,Nothing} = nothing, not_null::Bool = false) :: Nothing
  SearchLight.query("ALTER TABLE $table_name ADD $(SearchLight.Migration.column(name, column_type, default = default, limit = limit, not_null = not_null))")

  nothing
end

function SearchLight.Migration.drop_table(name::Union{String,Symbol}) :: Nothing
  SearchLight.query("DROP TABLE $name")

  nothing
end

function SearchLight.Migration.remove_column(table_name::Union{String,Symbol}, name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: Nothing
  SearchLight.query("ALTER TABLE $table_name DROP COLUMN $name $options")

  nothing
end

function SearchLight.Migration.remove_index(name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: Nothing
  SearchLight.query("DROP INDEX $name $options")

  nothing
end

function SearchLight.Migration.create_sequence(name::Union{String,Symbol}) :: Nothing
  SearchLight.query("CREATE SEQUENCE $name")

  nothing
end

function SearchLight.Migration.create_sequence(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: Nothing
  SearchLight.Migration.create_sequence(sequence_name(table_name, column_name))
end

function sequence_name(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: String
  string(table_name) * "__" * "seq_" * string(column_name)
end

function SearchLight.Migration.remove_sequence(name::Union{String,Symbol}, options::Union{String,Symbol}) :: Nothing
  SearchLight.query("DROP SEQUENCE $name $options")

  nothing
end

function SearchLight.Migration.remove_sequence(table_column_name::Tuple{Union{String,Symbol},Union{String,Symbol}}, options::String = "") :: Nothing
  SearchLight.Migration.remove_sequence(sequence_name(string(table_column_name[1]), string(table_column_name[2])), options)
end

function SearchLight.Migration.constraint(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: String
  string("CONSTRAINT ", SearchLight.index_name(table_name, column_name))
end

function SearchLight.Migration.nextval(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}) :: String
  "NEXTVAL('$(sequence_name(table_name, column_name) )')"
end

function SearchLight.Migration.column_id_sequence(table_name::Union{String,Symbol}, column_name::Union{String,Symbol})
  SearchLight.query("ALTER SEQUENCE $(sequence_name(table_name, column_name)) OWNED BY $table_name.$column_name")
end

#### TRANSACTIONS ####

function SearchLight.Transactions.begin_transaction() :: Nothing
  SearchLight.query("BEGIN")

  nothing
end

function SearchLight.Transactions.commit_transaction() :: Nothing
  SearchLight.query("COMMIT")

  nothing
end

function SearchLight.Transactions.rollback_transaction() :: Nothing
  SearchLight.query("ROLLBACK")

  nothing
end

#### GENERATOR ####

function SearchLight.Generator.FileTemplates.adapter_default_config(; database = SearchLight.config.app_env,
                                                                      host = "127.0.0.1",
                                                                      port = 5432,
                                                                      username = "postgres",
                                                                      password = "",
                                                                      env = SearchLight.config.app_env,
                                                                      env_val = """ENV["GENIE_ENV"]""") :: String
  """
  env: $env_val

  $env:
    adapter:  PostgreSQL
    host:     $host
    port:     $port
    database: $database
    username: $username
    password: $password
  """
end

end


text:
using DashBootstrapComponents

breadcrumb = dbc_breadcrumb(
    items = [
        Dict("label" => "Docs", "href" => "/docs", "external_link" => true),
        Dict(
            "label" => "Components",
            "href" => "/docs/components",
            "external_link" => true,
        ),
        Dict("label" => "Breadcrumb", "active" => true),
    ],
)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
ZDDs.jl
===

Implementation of ZDD (Zero-supressed Binary Decision Diagram) proposed by Shin-ichi Minato.
ZDD is a data structure which represents combination sets (family).

The operations below are supported:

- empty
- base
- subset0
- subset1
- change
- intersect
- union
- setdiff

example:
```julia
julia> includes("ZDDs.jl")
julia> using ZDDs
julia> zdd1 = tozdd(((1), (1,2)))
julia> zdd2 = tozdd(((2), (1,2), (2,3,4)))

julia> tofamily(union(zdd1, zdd2))
Set([Set([1]),Set([2,1]),Set([2]),Set([4,2,3])])

julia> tofamily(setdiff(zdd1, zdd2))
Set([Set([1])])

julia> tofamily(setdiff(zdd2, zdd1))
Set([Set([2]),Set([4,2,3])])

julia> tofamily(intersect(zdd1, zdd2))
Set([Set([2,1])])
```

reference:
Minato, S, 2001, Zero-suppressed BDDs and their applications,
International Journal on Software Tools for Technology

"""

module ZDDs

import Base.union
import Base.isempty
import Base.first
import Base.intersect
import Base.setdiff
import Base.length

export
	ZDD,
	tozdd,
	setdiff,
	intersect,
	union,
	Family,
	tofamily,
	print_size,
	subset0,
	subset1,
	change

"""
Node type
===
The type representing nodes, of which ZDDs are comprised.
Node type has 3 fields:

- top::Int
- child0::Node
- child1::Node
"""
type Node
	# label value
	top::Int
	# reference to the root node of the partial graph
	# representing combination sets which don't include *top*s.
	child0::Node
	# reference to the root node of the partial graph
	# representing combination sets which include *top*s.
	child1::Node

	function Node()
		new_node = new()
		new_node.child0 = new_node
		new_node.child1 = new_node
	end
	"""
	This constructor makes an instance with circular references to itself.
	Intented to used as constructor of "terminal" nodes.
	"""
	function Node(top::Int)
		new_node = new()
		new_node.top = top
		new_node.child0 = new_node
		new_node.child1 = new_node
	end

	"""
	This constructor makes an instance with the top value and
	the references to child nodes as given args.
	"""
	function Node(top::Int, child0::Node, child1::Node)
		new_node = new()
		new_node.top = top
		new_node.child0 = child0
		new_node.child1 = child1
		return new_node
	end
end

# Define terminal nodes.
# It is another considerable (and maybe safer) choice to define *Node* as an abstract type
# and to separate terminal nodes and normal nodes into different derived types.
# I thought that defining *Node* as a concrete type is better in view of performance,
# but I've not yet compared the performances of each versions.

# The top values of the terminal nodes are not referenced.

# The true terminal node.
const true_terminal = Node(typemin(Int))

# The false terminal node.
const false_terminal = Node(typemax(Int))

"""
ZDD type
"""
type ZDD
	root::Node
end

"""
returns an empty ZDD (which represents an empty combination set)
"""
function emptyzdd()
	return ZDD(false_terminal)
end

"""
returns a base ZDD (which represents a combination set with only an empty combination)
"""
function basezdd()
	return ZDD(true_terminal)
end

"""
main recursive function of tofamily
"""
function tofamily_sub_zdd(node::Node)
	if isempty(node)
		return Set{Set{Int}}()
	elseif isbase(node)
		return Set([Set{Int}([])])
	end
	child0 = tofamily_sub_zdd(node.child0)
	child1 = tofamily_sub_zdd(node.child1)
	for set in child1
		push!(set, node.top)
	end
	return union(child0, child1)
end

"""
Convert the combination set represented by *ZDD* into *Set{Set{Int}}*.
"""
function tofamily(zdd::ZDD)
	tofamily_sub_zdd(zdd.root)
end

"""
Array with the position currently referenced.
It is used in tozdd_sub! function.
"""
type ArrayWithPos
	parent::AbstractArray
	pos::Int
end

"""Returns true if *pos* of *arr* is over the length of *arr*."""
isempty(arr::ArrayWithPos) = ! (arr.pos <= length(arr.parent))

"""Get the element at *pos* from *arr*."""
first(arr::ArrayWithPos) = arr.parent[arr.pos]

"""Increment *pos* of *arr*."""
function inc!(arr::ArrayWithPos)
	arr.pos += 1
end

"""Decrement *pos* of *arr*."""
function dec!(arr::ArrayWithPos)
	arr.pos -= 1
end

"""
main recursive function of tozdd
"""
function tozdd_sub!(family)
	if isempty(family)
		return false_terminal
	elseif length(family) == 1
		if isempty(first(family))
			return true_terminal
		end
	end

	# Search the maximum value in *family*.
	# *s*s in *family* are already sorted, so checking the first values is enough.
	top = typemin(Int)
	for s in family
		if ! isempty(s)
			top = max(top, first(s))
		end
	end

	# The ver. commented out below is slower than the ver. upon (due to memory allocation within *filter* function)
	# top = mapreduce(first, max, filter(s -> !isempty(s), family))

	# Sort *family* in-place by the condition whether *s* includes *top*.
	# The reason not to use *sort* function is
	#   1) to sort *family* and
	#   2) to get the index of the first combination which includes *top* in sorted family
	# by one-path.
	f = 1
	e = length(family)
	while f < e
		sf = family[f]
		se = family[e]
		if isempty(sf) || first(sf) != top
			f += 1
		else
			if isempty(se) || first(se) != top
				family[f] = family[e]
				family[e] = sf
			else
				e -= 1
			end
		end
	end
	# At the end of the loop, each element in *family[1:f-1]* does not include *top*,
	# and that in *family[f:length(family)]* includes *top*.

	# Create *SubArray*s.
	f0 = sub(family, 1:f-1)
	f1 = sub(family, f:length(family))

	# Increment every *pos* of *s*s in *f1*.
	# After the incrementation, the values in *arr* with index larger than *pos*
	# are less than *top*.
	for s in f1
		inc!(s)
	end
	# Construct sub-ZDDs representing *f0* and *f1*.
	child0 = tozdd_sub!(f0)
	child1 = tozdd_sub!(f1)

	# After sub-ZDDs are constructed, get back *s*s in *f1* to their previous states.
	for s in f1
		dec!(s)
	end

	# Return the node with the maximum value in family as top and
	# the constructed sub-ZDDs as children.
	return getnode(top, child0, child1)
end

"""
Construct ZDD from family.
The family has to be represented as an *Iterable* of *Iterable*s of *Int*s.
"""
function tozdd(family)

	# Sort every set in *family* at first.
	# Because the max values of each set are repeatedly calculated,
	# it's more efficient to sort them at first.
	function preprocess_set(s)
		return ArrayWithPos(sort!(unique(s), rev=true), 1)
	end

	# *tozdd_sub!* assumes *family* does not contain duplication.
	preprocessed_family = unique(map(preprocess_set, family))
	root = tozdd_sub!(preprocessed_family)
	return ZDD(root)
end

"""
Returns true if zdd is an empty set.
"""
isempty(zdd::ZDD) = isempty(zdd.root)

isempty(p::Node) = is(p, false_terminal)

"""
Returns true if zdd is a base set.
"""
isbase(zdd::ZDD) = isbase(zdd.root)

isbase(p::Node) = is(p, true_terminal)

let
	global getnode
	global print_size
	unique_table = Dict{Tuple{Int, Node, Node}, Node}()
	"""
	Generates node for given *top* value, *child0* and *child1* subgraph.
	Using this method enables to keep each node unique.
	"""
	function getnode(top::Int, child0::Node, child1::Node)
		if isempty(child1)
			return child0
		else
			key = tuple(top, child0, child1)
			if haskey(unique_table, key)
				return unique_table[key]
			else
				node = Node(top, child0, child1)
				unique_table[key] = node
				return node
			end
		end
	end
	function print_size()
		print(length(unique_table))
		print("\n")
	end
end

"""

"""
function setdiff_sub(p::Node, q::Node)
	if isempty(p) | isempty(q)
		return p
	elseif is(p, q)
		return false_terminal
	elseif isbase(p)
		return setdiff_sub(p, q.child0)
	elseif isbase(q)
		return getnode(p.top, setdiff_sub(p.child0, q), p.child1)
	elseif p.top > q.top
		return getnode(p.top, setdiff_sub(p.child0, q), p.child1)
	elseif p.top < q.top
		return setdiff_sub(p, q.child0)
	else
		# p.top == q.top
		return getnode(
			p.top,
			setdiff_sub(p.child0, q.child0),
			setdiff_sub(p.child1, q.child1))
	end
end

"""
Returns the set of combinations which are in *p* and not in *q*.
"""
function setdiff(p::ZDD, q::ZDD)
	return ZDD(setdiff_sub(p.root, q.root))
end

function subset1_sub(p::Node, v::Int)
	if p.top < v
		return false_terminal
	end
	if p.top == v
		return p.child1
	end
	if p.top > v
		return getnode(
			p.top,
			subset1_sub(p.child0, v),
			subset1_sub(p.child1, v))
	end
end

"""
Select the set of combinations including *v* from *zdd*,
then delete *v* from them and return it.
"""
function subset1(zdd::ZDD, v::Int)
	return ZDD(subset1_sub(zdd.root, v))
end

function subset0_sub(p::Node, v::Int)
	if p.top < v
		return p
	elseif p.top == v
		return p.child0
	elseif p.top > v
		return getnode(
			p.top,
			subset0_sub(p.child0, v),
			subset0_sub(p.child1, v))
	end
end

"""
Select the combinations not including *v* from *zdd*.
"""
function subset0(zdd::ZDD, v::Int)
	return ZDD(subset0_sub(zdd.root, v))
end

function change_sub(p::Node, v::Int)
	if p.top < v
		return getnode(v, false_terminal, p)
	elseif p.top == v
		return getnode(v, p.child1, p.child0)
	elseif p.top > v
		return getnode(
			p.top,
			change_sub(p.child0, v),
			change_sub(p.child1, v))
	end
end

"""
Inverts *v* on each combinations in *zdd*.
"""
function change(zdd::ZDD, v::Int)
	return ZDD(change_sub(zdd.root, v))
end

function union_sub(p::Node, q::Node)
	if isempty(p)
		return q
	elseif isempty(q)
		return p
	elseif is(p, q)
		return p
	elseif isbase(q)
		return getnode(p.top, union_sub(p.child0, q), p.child1)
	elseif isbase(p)
		return getnode(q.top, union_sub(p, q.child0), q.child1)
	elseif p.top > q.top
		return getnode(p.top, union_sub(p.child0, q), p.child1)
	elseif p.top < q.top
		return getnode(q.top, union_sub(p, q.child0), q.child1)
	else
		return getnode(p.top, union_sub(p.child0, q.child0), union_sub(p.child1, q.child1))
	end
end

"""
Returns the union set of *p* and *q*.
"""
function union(p::ZDD, q::ZDD)
	return ZDD(union_sub(p.root, q.root))
end

function intersect_sub(p::Node, q::Node)
	if isempty(p) | isempty(q)
		return false_terminal
	elseif is(p, q)
		return p
	elseif p.top > q.top
		return intersect_sub(p.child0, q)
	elseif p.top < q.top
		return intersect_sub(p, q.child0)
	elseif p.top == q.top
		return getnode(p.top, intersect_sub(p.child0, q.child0), intersect_sub(p.child1, q.child1))
	end
end

"""
Returns the intersection set of *p* and *q*.
"""
function intersect(p::ZDD, q::ZDD)
	return ZDD(intersect_sub(p.root, q.root))
end

function length_sub(p::Node)
	if isempty(p)
		return 0
	elseif isbase(p)
		return 1
	else
		return length_sub(p.child0) + length_sub(p.child1)
	end
end

"""
Returns the number of combinations in *zdd*.
"""
function length(zdd::ZDD)
	return length_sub(zdd.root)
end

end


text:
function before_first_sunday(current_time::DateTime)
  return day(current_time)<Dates.Sunday && day(current_time)>=dayofweek(current_time)
end

function get_utc_time(current_time::DateTime, airport_code::ASCIIString)
  if airport_code=="ADL"
    # NB doesn't work for before 2am
    # First Sunday in October	First Sunday in April
    if month(current_time)>Dates.April && month(current_time)<Dates.October
      hour_offset=9
    elseif month(current_time)<Dates.April || month(current_time)>Dates.October
      hour_offset=10 # have daylight saving
    elseif month(current_time)==Dates.April
      # special case for April
      if before_first_sunday(current_time)
        hour_offset=10
      else
        hour_offset=9
      end
    elseif month(current_time)==Dates.October
      # special case for October
      if before_first_sunday(current_time)
        hour_offset=9
      else
        hour_offset=10
      end
    else
      error("logic error on month")
    end
    min_offset=30
  elseif airport_code=="HND"
    hour_offset=10
    min_offset=0
  else
    error("Invalid airport code")
  end
  return current_time-Dates.Hour(hour_offset)-Dates.Minute(min_offset)
end

function make_sure_data_loaded(current_time::DateTime)
  if (!haskey(got_temp_year, year(current_time)))
    # we don't have this year loaded
    # Open the file, if we cant give up
    # Station list in ftp://ftp.ncdc.noaa.gov/pub/data/inventories/ISH-HISTORY.TXT
    wanted_site="946720"
    wanted_year=string(year(current_time))
    wanted_file="$wanted_site-99999-$wanted_year.gz"
    if !isfile(wanted_file)
      # We need to fetch the file
      download("ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/$wanted_year/$wanted_file", wanted_file)
      if !isfile(wanted_file)
        error("Tried to download $source_file and that didn't work, giving up")
      end
    end
    # File is gzipped, open inside a pipe
    dt=readlines(`gzip -cd $wanted_file`)
    # The file is in a fixed format described at http://www1.ncdc.noaa.gov/pub/data/noaa/isd-lite/isd-lite-format.pdf
    # 2009 01 01 01   204    88 10127   220    87 -9999 -9999 -9999\n
    # yyyy mm dd hh  temp  dewp press widir wispd skyco rain1 rain6
    for thiselem in dt
      if thiselem[14:19]!=" -9999" # Handle the missing value
        push!(adelaide_temp_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[14:19])/10))
      end
      if thiselem[20:25]!=" -9999" # Handle the missing value
        push!(adelaide_dewp_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[20:25])/10))
      end
      if thiselem[38:43]!=" -9999" # Handle the missing value
        push!(adelaide_wind_speed_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[38:43])/10))
      end
    end
    # Read in the data without altering the offset??
    sort!(adelaide_temp_obs)
    sort!(adelaide_dewp_obs)
    sort!(adelaide_wind_speed_obs)
    got_temp_year[year(current_time)]=true
  end
end

function get_adelaide_temp(current_time::DateTime)
  make_sure_data_loaded(get_utc_time(current_time,"ADL"))
  # find the time immediately before the current time and return the adelaide airport observation
  return adelaide_temp_obs[searchsortedlast(adelaide_temp_obs,get_utc_time(current_time,"ADL"))].mtemp
  # For the moment, just return the average max temperature for the month
  # average_mean_max=[28.1,28.1,25.5,22.2,18.6,15.9,14.9,15.9,18.2,21.0,23.9,25.7]
  #  return average_mean_max[month(current_time)]
end

function get_adelaide_dewp(current_time::DateTime)
  make_sure_data_loaded(get_utc_time(current_time,"ADL"))
  # find the time immediately before the current time and return the adelaide airport observation
  return adelaide_dewp_obs[searchsortedlast(adelaide_dewp_obs,get_utc_time(current_time,"ADL"))].mtemp
end

function get_adelaide_wind_speed(current_time::DateTime)
  make_sure_data_loaded(get_utc_time(current_time,"ADL"))
  # find the time immediately before the current time and return the adelaide airport observation
  return adelaide_wind_speed_obs[searchsortedlast(adelaide_wind_speed_obs,get_utc_time(current_time,"ADL"))].mtemp
end

function make_sure_tokyo_loaded(current_time::DateTime)
  if !haskey(got_tokyo_year, year(current_time))
    # we don't have this year loaded
    # Open the file, if we cant give up
    wanted_site="476710" # RJTT Haneda Tokyo International
    wanted_year=string(year(current_time))
    wanted_file="$wanted_site-99999-$wanted_year.gz"
    if !isfile(wanted_file)
      # We need to fetch the file
      download("ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/$wanted_year/$wanted_file", wanted_file)
      if !isfile(wanted_file)
        error("Tried to download $source_file and that didn't work, giving up")
      end
    end
    # File is gzipped, open inside a pipe
    dt=readlines(`gzip -cd $wanted_file`)
    # The file is in a fixed format described at http://www1.ncdc.noaa.gov/pub/data/noaa/isd-lite/isd-lite-format.pdf
    # 2009 01 01 01   204    88 10127   220    87 -9999 -9999 -9999\n
    # yyyy mm dd hh  temp  dewp press widir wispd skyco rain1 rain6
    for thiselem in dt
      if thiselem[14:19]!=" -9999" # Handle the missing value
        push!(tokyo_temp_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[14:19])/10))
      end
      if thiselem[20:25]!=" -9999" # Handle the missing value
        push!(tokyo_dewp_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[20:25])/10))
      end
      if thiselem[38:43]!=" -9999" # Handle the missing value
        push!(tokyo_wind_speed_obs,met_obs(DateTime(int(thiselem[1:4]),int(thiselem[6:7]),int(thiselem[9:10]),int(thiselem[12:13]),0,0,0),int(thiselem[38:43])/10))
      end
    end
    # Read in the data without altering the offset??
    sort!(tokyo_temp_obs)
    sort!(tokyo_dewp_obs)
    sort!(tokyo_wind_speed_obs)
    got_tokyo_year[year(current_time)]=true
  end
end

function get_tokyo_temp(current_time::DateTime)
  make_sure_tokyo_loaded(get_utc_time(current_time,"HND"))
  # find the time immediately before the current time and return the adelaide airport observation
  return tokyo_temp_obs[searchsortedlast(tokyo_temp_obs,get_utc_time(current_time,"HND"))].mtemp
end

function get_tokyo_dewp(current_time::DateTime)
  make_sure_tokyo_loaded(get_utc_time(current_time,"HND"))
  # find the time immediately before the current time and return the adelaide airport observation
  return tokyo_dewp_obs[searchsortedlast(tokyo_dewp_obs,get_utc_time(current_time,"HND"))].mtemp
end

function get_tokyo_wind_speed(current_time::DateTime)
  make_sure_tokyo_loaded(get_utc_time(current_time,"HND"))
  # find the time immediately before the current time and return the adelaide airport observation
  return tokyo_wind_speed_obs[searchsortedlast(tokyo_wind_speed_obs,get_utc_time(current_time,"HND"))].mtemp
end

got_temp_year=Dict{Int,Bool}()
got_tokyo_year=Dict{Int,Bool}()

type met_obs
 mdate::DateTime
 mtemp::Float64
end
Base.isless(x::met_obs,y::met_obs) = isless(x.mdate,y.mdate) # this is how we sort them
Base.isless(x::DateTime,y::met_obs) = isless(x,y.mdate) # make searchsorted and friends work
adelaide_temp_obs=met_obs[]
adelaide_dewp_obs=met_obs[]
adelaide_wind_speed_obs=met_obs[]
tokyo_temp_obs=met_obs[]
tokyo_dewp_obs=met_obs[]
tokyo_wind_speed_obs=met_obs[]

max_stars_repo_name: softgrow/mfsmodel-weather

max_issues_repo_name: softgrow/mfsmodel-weather

max_forks_repo_name: softgrow/mfsmodel-weather