text:
#!/usr/bin/env julia

import Celeste: AccuracyBenchmark
import Celeste: ArgumentParse
import Celeste: SDSSIO

const OUTPUT_DIRECTORY = joinpath(splitdir(Base.source_path())[1], "output")

parser = ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "run",
    help="SDSS run #",
    arg_type=Int,
)
ArgumentParse.add_argument(
    parser,
    "camcol",
    help="SDSS camcol #",
    arg_type=Int,
)
ArgumentParse.add_argument(
    parser,
    "field",
    help="SDSS field #",
    arg_type=Int,
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

srand(12345)

rcf = SDSSIO.RunCamcolField(
    parsed_args["run"],
    parsed_args["camcol"],
    parsed_args["field"],
)
strategy = SDSSIO.PlainFITSStrategy(AccuracyBenchmark.SDSS_DATA_DIR)
images = SDSSIO.load_field_images(strategy, [rcf])
@assert length(images) == 5

if !isdir(OUTPUT_DIRECTORY)
    mkdir(OUTPUT_DIRECTORY)
end
filename = joinpath(
    OUTPUT_DIRECTORY,
    @sprintf("sdss_%s_%s_%s_images.fits", rcf.run, rcf.camcol, rcf.field),
)
AccuracyBenchmark.save_images_to_fits(filename, images)
AccuracyBenchmark.append_hash_to_file(filename)


text:
#!/usr/bin/env julia

# The main purpose of this file is to read the config file and copy all of the relevant static
# files into this directory, write out the proper wavelength files, and make sure we can continue
# with the run properly.

using ArgParse

s = ArgParseSettings(description="Initialize a new project directory with the appropriate files. Can also be used to update RADMC-3D input files after making changes to config.yaml")
@add_arg_table s begin
    "--version"
    help = "Print a the version number and exit."
    action = :store_true
    "-M"
    help = "Print out the list of files this program generates."
    action = :store_true
    "--new-project"
    help = "Copy a stock configuration file to this directory. Can be 'standard', 'truncated', 'cavity'"
    default = "no"
    "--fit-every"
    help = "Print out an exclude array that should be manually copied to config.yaml so that only every n-th channel is fit in venus.jl"
    arg_type = Int
    default = 0
    "--prior"
    help = "Copy a default prior.jl template to this current directory, for user modification. CAREFUL! See the docs about this functionality before invoking."
    action = :store_true
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    "--vstart", "-s"
    help = "Starting velocity in km/s."
    arg_type = Float64
    "--vend", "-e"
    help = "Ending velocity in km/s."
    arg_type = Float64
    "--nvel", "-n"
    help = "Number of velocity channels."
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)

assets_dir = Pkg.dir("DiskJockey") * "/assets/"

# The user is going to start modeling a new disk, so copy in the new configuration file.
if parsed_args["new-project"] != "no"
    model = parsed_args["new-project"]

    cp(assets_dir * "config.$(model).yaml", pwd() * "/config.yaml")
    cp(assets_dir * "InitializeWalkers.$(model).ipynb", pwd() * "/InitializeWalkers.ipynb")
    cp(assets_dir * "Makefile", pwd() * "/Makefile")

    println("Copied default config.yaml, InitializeWalkers.ipynb, and Makefile for the $model model to current working directory.")
    println("Exiting")
    quit()
end

if parsed_args["prior"]
    cp(assets_dir * "prior.jl", pwd() * "/prior.jl")
    println("prior.jl copied to current directory. If you didn't mean to do this, delete this file now. Otherwise, edit this file with your favorite text editor to enforce the prior you would like to see for this specific disk. Please see the docs for more information.")
    println("Exiting")
    quit()
end

using DiskJockey.constants

# This is just for new users to test that the package is successfully installed.
if parsed_args["version"]
    println("Your DiskJockey scripts are successfully linked.")
    println("You are running DiskJockey $DISKJOCKEY_VERSION")
    println("Exiting")
    quit()
end

import YAML
config = YAML.load(open(parsed_args["config"]))

if parsed_args["M"]

    str = "amr_grid.inp camera_wavelength_micron.inp gas_temperature.inp gas_velocity.inp lines.inp microturbulence.inp radmc3d.inp wavelength_micron.inp"

    species = config["species"]
    transition = config["transition"]

    if species == "12CO"
        str *= " molecule_co.inp numberdens_co.inp"
    elseif species ==  "13CO"
        str *= " molecule_13co.inp numberdens_13co.inp"
    elseif species == "C18O"
        str *= " molecule_c18o.inp numberdens_c18o.inp"
    end

    println(str)
    quit()
end

using DiskJockey.model
using DiskJockey.visibilities
using HDF5

fit_every = parsed_args["fit-every"]
if fit_every != 0
    dvarr = DataVis(config["data_file"])
    nchan = length(dvarr)
    println("Dataset contains a total of $nchan channels.")

    to_fit = Int[i for i=1:fit_every:nchan]

    # which channels of the dset to exclude
    exclude = filter(x->(!in(x, to_fit)), Int[i for i=1:nchan])

    println("To fit only every $(fit_every)-nd/rd/th channel, place the following line in your `config.yaml file.`")
    println("exclude : $exclude")
    println("Exiting")
    quit()
end

# Create an output directory that will be useful for the results.
out_base = config["out_base"]
if !ispath(out_base)
    println("Creating ", out_base)
    mkdir(out_base)
end

# vstart = parsed_args["vstart"]
# vend = parsed_args["vend"]
# nvel = parsed_args["nvel"]
#
# # If we have specified the velocities, ignore the data files and Doppler shift
# # and create an evenly spaced array of velocities
# vels = linspace(vstart, vend, nvel) # [km/s]
#
# lam0 = lam0s[species*transition]
#
# # convert velocities to wavelengths
# shift_lams = lam0 * (vels/c_kms + 1)

# Do everything common to both models
pars = convert_dict(config["parameters"], config["model"])
grd = config["grid"]
grid = Grid(grd)

# Are we using dust or gas?
# gas
if config["gas"]
    run(`cp $(assets_dir)radmc3d.inp.gas radmc3d.inp`)

    # Need this as a dummy file, apparently.
    run(`cp $(assets_dir)wavelength_micron.inp .`)

    # Copy appropriate molecule files.

    # Which species?
    species = config["species"]
    transition = config["transition"]
    if species == "12CO"
        run(`cp $(assets_dir)molecule_co.inp .`)
        run(`cp $(assets_dir)lines_co.inp lines.inp`)
    elseif species ==  "13CO"
        run(`cp $(assets_dir)molecule_13co.inp .`)
        run(`cp $(assets_dir)lines_13co.inp lines.inp`)
    elseif species == "C18O"
        run(`cp $(assets_dir)molecule_c18o.inp .`)
        run(`cp $(assets_dir)lines_c18o.inp lines.inp`)
    end

    println("Copied over RADMC-3D gas input files.")

    vel = pars.vel # [km/s]
    npix = config["npix"] # number of pixels

    vstart = parsed_args["vstart"]
    vend = parsed_args["vend"]
    nvel = parsed_args["nvel"]

    if vstart != nothing
        # If we have specified the velocities, ignore the data files and Doppler shift
        # and create an evenly spaced array of velocities
        vels = linspace(vstart, vend, nvel) # [km/s]

        lam0 = lam0s[species*transition]

        # convert velocities to wavelengths
        shift_lams = lam0 * (vels/c_kms + 1)
    else
        # read the wavelengths for all data channels
        fid = h5open(config["data_file"], "r")
        freqs = read(fid["freqs"]) # [Hz]
        # Convert from Hz to wavelengths in μm
        lams = cc ./freqs * 1e4 # [μm]
        close(fid)

        # Doppler shift the dataset wavelength according to the velocity in the parameter file
        beta = vel/c_kms # relativistic Doppler formula
        shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    # Doppler shift the dataset wavelength according to the velocity in the parameter file
    # beta = vel/c_kms # relativistic Doppler formula
    # shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]

    # Write everything to the current working directory
    write_grid("", grid)
    write_model(pars, "", grid, species)
    write_lambda(shift_lams, "")
    println("Wrote gas model input files for RADMC-3D.")

# dust
elseif config["dust"]

    # Copy the appropriate files for dust synthesis
    run(`cp $(assets_dir)radmc3d.inp.dust radmc3d.inp`)

    # Write the dust distribution
    write_grid("", grid)
    write_dust(pars, "", grid)

    println("Wrote dust grid and density file for RADMC-3D.")

    # ff = assets_dir * "radmc3d.inp.dust"
    # run(`cp $ff radmc3d.inp`)
    # ff = assets_dir * "stars.inp"
    # run(`cp $ff .`)
    # ff = assets_dir * "wavelength_micron.inp"
    # run(`cp $ff .`)
    # ff = assets_dir * "dustkappa_silicate.inp"
    # run(`cp $ff .`)
    # ff = assets_dir * "dustopac.inp"
    # run(`cp $ff .`)
    #
    # println("Copied over RADMC-3D dust input files.")
end

ext: jl

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
PastaQ.state(optimizer, model::LPDO) =
  optimizerstate(optimizer, getparameters(model))

PastaQ.state(optimizer, model::MPS) =
  PastaQ.state(optimizer, LPDO(model))

PastaQ.state(optimizer, model::MPO) =
  PastaQ.state(optimizer, LPDO(model))

"""
    update!(model, grads, optimizer)

Update a tensor network model
"""
function PastaQ.update!(model, grads, optimizer)
  L = deepcopy(model)
  opt, st = optimizer
  θ, ∇ = getparameters(L; ∇ = copy(grads))
  st, θ′ = Optimisers.update(opt, st, θ, ∇)

  setparameters!(L, θ′)
  model.X[:] = L.X
  return model
end

function getparameters(L::LPDO; ∇ = nothing)
  X = L.X
  N = length(X)
  θ = Complex{Float64}[]
  λ = Complex{Float64}[]
  for j in 1:N
    append!(θ, vec(ITensors.array(X[j])))
    !isnothing(∇) && append!(λ, vec(ITensors.array(∇[j])))
  end
  !isnothing(∇) && return θ, λ
  return θ
end

function setparameters!(L::LPDO, θ::Vector)
  X = L.X
  Y = copy(X)
  N = length(X)

  cnt = 1
  for j in 1:N
    d = dims(X[j])
    tot_d = prod(d)
    ϑ = θ[cnt:cnt+tot_d-1]
    Y[j] = ITensors.itensor(reshape(ϑ, d), inds(X[j])...)
    cnt += tot_d
  end
  X[:] = Y
end


text:
struct CellLocation{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

struct CellSize{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

struct FaceLocation{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

struct MeshStructure
  dimension::Real
  dims::Array{Int,1}
  cellsize::CellSize
  cellcenters::CellLocation
  facecenters::FaceLocation
  corner::Array{Int,1}
  edge::Array{Int,1}
end

struct CellValue # {T<:Real}
  domain::MeshStructure
  value::Union{Array{<:Real}, DenseArray{Bool}} #Union{Array{T}, BitArray{}}
end

struct CellVector # {T<:Real}
  domain::MeshStructure
  xvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
  yvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
  zvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
end

struct FaceValue # {T<:Real}
  domain::MeshStructure
  xvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
  yvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
  zvalue::Union{Array{<:Real}, DenseArray{Bool}} # Array{T}
end

mutable struct BorderValue{T<:Real}
  a::Array{T}
  b::Array{T}
  c::Array{T}
  periodic::Bool
end

struct BoundaryCondition
  domain::MeshStructure
  left::BorderValue
  right::BorderValue
  bottom::BorderValue
  top::BorderValue
  back::BorderValue
  front::BorderValue
end


text:
function clean_whitespaces_and_newlines(contents)
    # Remove lines with only whitespaces:
    contents = join([line for line in split(contents, "\n") if (length(line)≠count(" ", line)) || isempty(line)], "\n")

    # Trim some \n's
    contents = replace(contents, "\n\n\n\n"=>"\n\n")
    contents = replace(contents, "\n\n\n\n"=>"\n\n")
    contents = replace(contents, "\n\n\n"=>"\n\n")
    if startswith(contents, "\n\n")
        chop(contents, head=1, tail=0)
    end
    return contents
end
comment_flags(::Union{Python,Cython}) = ("#",)
doc_string_flags(::Union{Python,Cython}) = ("\"\"\"",)
cond(ts,i,f1,f2) = ts.token_stream[f1][i]≠1 && ts.token_stream[f2][i]≠1

function remove_comments(file, ::Union{Python,Cython})
    println("   -----------------------")
    println("   Removing comments for file $(basename(file))")
    contents = open(f->read(f, String), file)

    println("   Single-line comments...")
    contents = join([rstrip(split(line, "#")[1]) for line in split(contents, "\n")], "\n")
    # Remove lines with only whitespaces:
    println("   Cleaning white spaces...")
    contents = clean_whitespaces_and_newlines(contents)

    contents = join([line for line in split(contents, "\n")], "\n")

    flag_set = FlagSet([
    FlagPair{GreedyType}(
    StartFlag("'''", ["\n", " ", ""], ["\n", " ", ""]),
    StopFlag("'''", ["\n", " ", ""], ["\n", " ", ""])
    ),
    FlagPair{GreedyType}(
    StartFlag("\"\"\"", ["\n", " ", ""], ["\n", " ", ""]),
    StopFlag("\"\"\"", ["\n", " ", ""], ["\n", " ", ""])
    )])
    println("   Creating token stream...")
    ts = TokenStream(contents, flag_set)
    println("   Finished creating token stream...")

    fs = ("'''-'''", "\"\"\"-\"\"\"")
    i_code = [i for i in eachindex(ts.code) if cond(ts, i, fs[1],fs[2])]
    contents = join([SubString(ts.code, i, i) for i in i_code])

    println("   Cleaning white spaces...")
    contents = clean_whitespaces_and_newlines(contents)

    open(file, "w") do io
        print(io, contents)
    end
end


text:
# Chapter 2.6 コンストラクタ
using CSV
using DataFrames
using Plots
using Measures

struct HousePriceData
    df
end

# コンストラクタ
function HousePriceData(path::AbstractString)
    df = CSV.read(path, DataFrame)
    HousePriceData(df)
end

# 建物の等級(MSSubClass)を限定した散布図(横軸GrLibArea, 縦軸SalePrice)をプロットするメソッド
# titles: タイトル(グラフ数のリスト)
# levels: 建物の東急(グラフ数のリスト)
function plot_scatter(data::HousePriceData, titles, levels)
    xrange = extrema(data.df[!, :GrLivArea])
    yrange = extrema(data.df[!, :SalePrice])

    plts = []
    for (title, lv) in zip(titles, levels)
        df_plot = data.df[data.df.MSSubClass .== lv, :]
        fig = scatter(df_plot[!, :GrLivArea],
                      df_plot[!, :SalePrice],
                      title = title,
                      xlims = xrange,
                      ylims = yrange,
                      label = "")
        xlabel!("GrLivArea")
        ylabel!("SalePrice")
        push!(plts, fig)
    end
    plot(plts..., layout = (:, 2), size = [600, 600], margin = 5mm)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

# some preliminary code for constructing hierarchical elements

abstract Hierarchical <: Element

bin(n, k) = prod([(n + 1 - i)/i for i=1:k])

"""
Return Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]

Parameters
----------
n :: Int
    order of polynomial

Returns
-------
function
    Legendgre polynomial of order n in interval ξ ∈ [-1, 1]
"""
function get_legendre_polynomial(n::Int)
    P(xi) = 2^n*sum([xi.^k*bin(n, k)*bin(1/2*(n+k-1), n) for k=0:n])
    P
end

"""
Return derivative of Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]
"""
function get_legendre_polynomial_derivative(n::Int)
    dP(xi) = 2^n*sum([k*xi.^(k-1)*bin(n, k)*bin(1/2*(n+k-1), n) for k=1:n])
    dP
end

"""
Return Legendgre polynomial of order n to inverval ξ ∈ [1, 1].

Parameters
----------
n :: Int
    order of polynomial

Returns
-------
function
    Legendgre polynomial of order n in interval ξ ∈ [-1, 1]

Notes
-----
Uses Bonnet's recursion formula. See
https://en.wikipedia.org/wiki/Legendre_polynomials
"""
function get_legendre_polynomial_recursive(n)
    if n == 0
        P(xi) = 1
    elseif n == 1
        P(xi) = xi
    else
        Pm1 = get_legendre_polynomial(n-1)
        Pm2 = get_legendre_polynomial(n-2)
        P(xi) = 1/n*((2*n-1)*xi*Pm1(xi) - (n-1)*Pm2(xi))
    end
    return P
end

"""
Return derivative of Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]
"""
function get_legendre_polynomial_derivative_recursive(n)
    if n == 0
        P(xi) = 0
    elseif n == 1
        P(xi) = 1
    else
        Pm1 = get_legendre_polynomial_derivative(n-1)
        Pm2 = get_legendre_polynomial_derivative(n-2)
        P(xi) = 1/(n-1)*( (2*(n-1)+1)*xi.*Pm1(xi) - (n+1-1)*Pm2(xi))
    end
    return P
end

"""
Return hierarchical shape function of order N
"""
function get_hierarchial_basis(n)
    if n == 1
        N(xi) = 1/2*(1 - xi)
    elseif n == 2
        N(xi) = 1/2*(1 + xi)
    else
        j = n-1
        Pj = get_legendre_polynomial(j)
        Pjm2 = get_legendre_polynomial(j-2)
        N(xi) = 1/sqrt(2*(2*j-1))*(Pj(xi) - Pjm2(xi))
    end
    return N
end

"""
Return derivative of hierarchical shape function of order N
"""
function get_hierarchial_basis_derivative(n)
    if n == 1
        dN(xi) = -1/2
    elseif n == 2
        dN(xi) = 1/2
    else
        j = n-1
        Pj = get_legendre_polynomial_derivative(j)
        Pjm2 = get_legendre_polynomial_derivative(j-2)
        dN(xi) = 1/sqrt(2*(2*j-1))*(Pj(xi) - Pjm2(xi))
    end
    return dN
end

"""
Set degree of hierarchical element
"""
function set_degree(el::Hierarchical, degree)
    el.degree = degree
end

"""
Get degree of hierarchical element
"""
function get_degree(el::Hierarchical)
    el.degree
end

"""
Hierarchical 1d segment element.
"""
type PSeg <: Hierarchical
    connectivity :: Array{Int, 1}
    fields :: Dict{Any, Any}
    degree :: Int
end
PSeg(connectivity) = PSeg(connectivity, Dict{Any,Any}(), 1)
get_number_of_basis_functions(el::Type{PSeg}) = 2
get_number_of_basis_functions(el::PSeg) = 2 + el.degree - 1
get_element_dimension(el::Type{PSeg}) = 1
function get_basis(el::PSeg, xi)
    m = get_number_of_basis_functions(el)
    out = zeros(m)
    for n=1:m
        N = get_hierarchial_basis(n)
        out[n] = N(xi[1])
    end
    return out
end
function get_dbasisdxi(el::PSeg, xi)
    m = get_number_of_basis_functions(el)
    out = zeros(m)
    for n=1:m
        dN = get_hierarchial_basis_derivative(n)
        out[n] = dN(xi[1])
    end
    return out
end


text:
module GEOTRACES

using OceanographyCruises
using NCDatasets
using Unitful
using Dates
using Match
using Measurements
using DataFrames

const metadatakeys=("lat", "lon", "depth", "cruise", "station", "date") # default

GEOTRACES_IDP21_DiscreteSamples_path() = get(ENV, "GEOTRACES_IDP21_PATH", joinpath(homedir(), "Data/GEOTRACES/GEOTRACES_IDP2021_v1/seawater/netcdf/GEOTRACES_IDP2021_Seawater_Discrete_Sample_Data_v1/GEOTRACES_IDP2021_Seawater_Discrete_Sample_Data_v1.nc"))

include("helper_functions.jl")

"""
    cruisetrack(cruise_name)

Construct GEOTRACES `CruiseTrack` from cruise name.
"""
function cruisetrack(ds::Dataset, cruise)
    ikeep = findall(list_of_cruises(ds) .== cruise)
    ikeep = CartesianIndices((1, length(ds[varname("lat")])))[ikeep]
    metadatakeys=("lat","lon","cruise","station","date")
    GEOTRACESmetadatakeys = varname.(metadatakeys)
    metadata = [metadatakeyvaluepair(ds[k], ikeep) for k in GEOTRACESmetadatakeys]
    m = (; metadata...)
    any(m.cruise .≠ cruise) && error("Not the right cruise!")
    stations = [Station(date=d,lat=y,lon=x,name=string(s)) for (d,y,x,s) in zip(m.date, m.lat, m.lon, m.station)]
    return CruiseTrack(name=cruise, stations=stations)
end

"""
    list_of_cruises()

List of GEOTRACES cruises.
"""
list_of_cruises(ds::Dataset) = reduce.(string, filter.(!=('\0'), eachcol(ds["metavar1"].var[:,:])))

"""
    list_of_stations()

List of GEOTRACES stations.
"""
list_of_stations(ds::Dataset) = string.(reduce.(string, filter.(!=('\0'), eachcol(ds["metavar2"].var[:,:]))))

"""
transect(tracer::String, cruise::String)

The `Transect` of observations of tracer `tracer` along cruise `cruise`.
"""
function transect(ds::Dataset, tracer::String; cruise::String, QCmax=1)
    ts = transects(ds, tracer, QCmax=QCmax)
    i = findall(ts.cruises .== cruise)
    length(i) == 1 ? i = i[1] : error("Multiple identical cruises is an issue!")
    return ts.transects[i]
end

"""
transects(tracer::String, cruise::String)

The `Transects` of observations of tracer `tracer`.
"""
function transects(obs::DataFrame, tracer::String)
    ts = Transect{eltype(getproperty(obs, tracer))}[]
    df_bycruise = groupby(obs, :cruise)
    for (cruise, cruise_df) in pairs(df_bycruise)
        df_byprofile = groupby(cruise_df, [:lat, :lon, :station, :date])
        profiles = [DepthProfile(station = Station(name=string(k.station), lat=k.lat, lon=k.lon, date=k.date),
                                 depths=ustrip.(v.depth),
                                 values=v[:,tracer]) for (k,v) in pairs(df_byprofile)]
        push!(ts, Transect(tracer=tracer, cruise=cruise.cruise, profiles=profiles))
    end
    return Transects(tracer=tracer, cruises=[k.cruise for k in keys(df_bycruise)], transects=ts)
end
function transects(ds::Dataset, tracer::String; QCmax=1)
    obs = observations(ds, tracer; metadatakeys, QCmax)
    return transects(obs, tracer)
end

#function cruisetrack(ds::Dataset, cruise::String)
#    Istations = findall(list_of_cruises(ds) .== cruise)
#    stations_names = list_of_stations(ds)
#    stations = Station[]
#
#    for istation in Istations
#        station = Station(name=stations_names[istation], lat=ds["latitude"][istation], lon=ds["longitude"][istation])
#        push!(stations, station)
#    end
#
#    return CruiseTrack(name=cruise, stations=stations)
#end

"""
    observations(tracer1, tracer2, tracer3, ...)

Returns the GEOTRACES observations of the given tracers.

### Example usage

```
x, y, z, ... = observations(tracer1, tracer2, tracer3, ...)
x = observations(tracer1)
```

### Quality control flags

By default, only observations with quality control flag of `1` are kept.
If you want more data, you must define keyword-argument `QCmax`.

```
x, y, z, ... = observations(tracer1, tracer2, tracer3, ...; QCmax = 2)
x = observations(tracer1; QCmax = 3)
```
"""
function observations(ds::Dataset, tracers::String...; metadatakeys=metadatakeys, QCmax=1)
    dfs = (observations(ds, tracer; metadatakeys, QCmax) for tracer in tracers)
    return innerjoin(dfs..., on=collect(metadatakeys))
end
function observations(ds::Dataset, tracer::String; metadatakeys=metadatakeys, QCmax=1)
    f = unitfunction(ds[varname(tracer)].attrib["units"])
    var = ds[varname(tracer)]
    v = var[:]
    qcv = ds[qcvarname(tracer)].var[:] .|> Char
    ikeep = findall(parse.(Int, qcv) .≤ QCmax)
    ikeep = CartesianIndices(size(qcv))[ikeep]
    GEOTRACESmetadatakeys = varname.(metadatakeys)
    metadata = (metadatakeyvaluepair(ds[k], ikeep) for k in GEOTRACESmetadatakeys)
    df = DataFrame(metadata..., Symbol(tracer)=>f.(float.(view(v, ikeep))))
    #mdf = attach_metadata(df, (GEOTRACESvarname=name(var),))
    return df
end

"""
    qualitycontrols(tracer1, tracer2, tracer3, ...)

Returns the GEOTRACES quality control flag of the given tracers.

GEOTRACES's `Char`s are converted to `Int` by this function.

From GEOTRACES:
```
1 = good quality
2 = not evaluated, not available or unknown quality
3 = questionable/suspect quality
4 = bad quality
9 = missing data
```

Note that there should be no `9` because `missing` data is skipped
by the `observations` function.
"""
function qualitycontrols(ds::Dataset, tracers::String...)
    vs = ((ds[varname(tracer)][:] for tracer in tracers)...,)
    ikeep = findall(i -> !any(ismissing.(getindex.(vs, i))), eachindex(vs[1]))
    qcvs = ((ds[qcvarname(tracer)].var[:] .|> Char for tracer in tracers)...,)
    return ((parse.(Int, view(qcv, ikeep)) for qcv in qcvs)...,)
end
function qualitycontrols(ds::Dataset, tracer::String)
    v = ds[varname(tracer)][:]
    ikeep = findall(!ismissing, v)
    qcv = ds[qcvarname(tracer)].var[:] .|> Char
    return parse.(Int, view(qcv, ikeep))
end

"""
    metadata(tracer1, tracer2, tracer3, ...)

Returns the GEOTRACES metadata for given tracers.
"""
function metadata(ds::Dataset, tracers::String...; metadatakeys=("lat", "lon", "depth"), QCmax=1)
    qcvs = ((ds[qcvarname(tracer)].var[:] .|> Char for tracer in tracers)...,)
    ikeep = findall(i -> all((parse.(Int, getindex.(qcvs, i))) .≤ QCmax), eachindex(qcvs[1]))
    ikeep = CartesianIndices(size(qcvs[1]))[ikeep]
    GEOTRACESmetadatakeys = varname.(metadatakeys)
    metadata = [metadatakeyvaluepair(ds[k], ikeep) for k in GEOTRACESmetadatakeys]
    namedmetadata = (; metadata...)
    return namedmetadata
end

# TODO treatment of uncertainty
# Not all variable have a "standard deviation" in GEOTRACES
# But I could create one from
# - the number of significant digits of each tracer per cruise?
# - the STD along each cruise?
# - the total STD?
function standarddeviations(ds::Dataset, tracer::String)
    _, s = observations(ds, varname(tracer), stdvarname(tracer))
    return s
end

function standarddeviations(ds::Dataset, tracers::String...)
    x_and_s = observations(ds, varname.(tracers)..., stdvarname.(tracers)...)
    return x_and_s[length(tracers)+1:end]
end

function observations_with_std(ds::Dataset, tracer::String)
    x, s = observations(ds, varname(tracer), stdvarname(tracer))
    u = unit(x[1])
    return @. (ustrip(x) ± ustrip(s)) * u
end
function observations_with_std(ds::Dataset, tracers::String...)
    x_and_s = observations(ds, varname.(tracers)..., stdvarname.(tracers)...)
    n = length(tracers)
    us = ((unit(x[1]) for x in x_and_s[1:n])...,)
    return (((ustrip.(x_and_s[i]) .± ustrip.(x_and_s[i+n])) * us[i] for i in 1:n)...,)
end

# open and close ds if not provided
for f in [:observations, :metadata, :transect, :transects, :list_of_cruises,
         :list_of_stations, :variable, :qualitycontrols, :cruisetrack,
         :standarddeviations, :matchingvariables, :observations_with_std]
    @eval begin
        $f(args...; kwargs...) =
        Dataset(GEOTRACES_IDP21_DiscreteSamples_path(), "r") do ds
             $f(ds, args...; kwargs...)
        end
    end
end

end # module


text:
import Base: +, -, *, reshape, size
import Base.Broadcast: broadcasted, Broadcasted, BroadcastStyle

"""
    Zeros()
    Zeros(size...)
    Zeros(Type, size...)

Acts as a stand-in for an array of zeros that can be
used during training which is ignored by the optimisers.

Useful to turn bias off for a forward pass of a layer.

## Examples

```julia
julia> Flux.Zeros(3,3)
3×3 Flux.Zeros{Bool,2}:
 false  false  false
 false  false  false
 false  false  false

julia> Flux.Zeros(Float32, 3,3)
3×3 Flux.Zeros{Float32,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> rand(3,3) .+ Flux.Zeros()
3×3 Array{Float64,2}:
 0.198739  0.490459  0.785386
 0.779074  0.39986   0.66383
 0.854981  0.447292  0.314497

julia> bias_less_conv = Conv((2,2), 1=>3, bias = Flux.Zeros())
Conv((2, 2), 1=>3)
```
"""
struct Zeros{T,N} <: AbstractArray{T,N}
  size::Tuple
end

Zeros(::Type{T}, sz...) where T = Zeros{T,length(sz)}(sz)
Zeros(sz::Integer...) = Zeros(Bool, sz...)

Base.size(xs::Zeros) = xs.size
Base.axes(xs::Zeros) = Base.OneTo.(size(xs))

Base.IndexStyle(::Type{<:Zeros}) = IndexLinear()

Base.getindex(xs::Zeros{T,N}, I::Int) where {T,N} = zero(T)
Base.getindex(xs::Zeros{T,N}, inds::Union{Base.OneTo, Base.UnitRange}) where {T,N} =
              Zeros(T, length(inds))

Base.collect(xs::Zeros{T,N}) where {T,N} = fill(zero(T), size(xs))

@adjoint reshape(xs::Zeros{T}, dims...) where T =
                reshape(xs, dims...), _ -> nothing

# Define basic ops
for f in (:+, :-)
  @eval @inline function $f(a::Union{AbstractArray{<:Number}, Zeros}, b::Zeros)
    @assert size(a) == size(b) throw(DimensionMismatch("dimensions must match"))
    a
  end
end

+(a::Zeros, b::AbstractArray) = b + a
-(a::Zeros, b::AbstractArray) = -b + a

Base.copy(xs::Zeros{T,N}) where {T,N} = xs

# Define broadcasting behaviour
for op in (:+, :-)
  @eval function broadcasted(::typeof($op), a::AbstractArray, b::Zeros)
    bs = Broadcast.broadcast_shape(size(a), size(b))
    size(a) == bs && return a
    sz = similar(a, bs)
    sz .= a
  end
end

broadcasted(::typeof(+), a::Zeros, b::AbstractArray) = broadcasted(+, b, a)
broadcasted(::typeof(-), a::Zeros, b::AbstractArray) = broadcasted(+, -b, a)

function broadcasted(::typeof(*), a::AbstractArray, b::Zeros)
  Zeros(Broadcast.broadcast_shape(size(a), size(b))...)
end

broadcasted(::typeof(*), a::Zeros, b::AbstractArray) = broadcasted(*, b, a)

for op in (:+, :-, :*)
  @eval broadcasted(::typeof($op), a::Zeros, b::Zeros) = Zeros(Broadcast.broadcast_shape(size(a), size(b))...)
end

# Some opportunities to avoid scalar indexing, intermediaries
# Since it replicates a little of what we expect Base to do,
# it should be possible to remove in the future, but for now,
# these help with performance.
broadcasted(::typeof(+), a::AbstractArray, b::Zeros{T,0}) where T = a
broadcasted(::typeof(+), a::Zeros{T,0}, b::AbstractArray) where T = b
broadcasted(::typeof(-), a::AbstractArray, b::Zeros{T,0}) where T = a
broadcasted(::typeof(-), a::Zeros{T,0}, b::AbstractArray) where T = -b
broadcasted(::typeof(*), a::AbstractArray, b::Zeros{T,0}) where T = zero(a)
broadcasted(::typeof(*), a::Zeros{T,0}, b::AbstractArray) where T = zero(b)
broadcasted(::typeof(/), a::Zeros{T,0}, b::AbstractArray) where T = zero(b)


text:
function lanczos_kernel{NumType <: Number}(x::NumType, a::Float64)
    abs(x) < a ? sinc(x) * sinc(x / a): zero(NumType)
end

function sinc_with_derivatives{NumType <: Number}(x::NumType)
    x_pi = pi * x
    sinc_x = sinc(x)
    sinc_x_d = (cos(x_pi) - sinc_x) / x
    sinc_x_h = -pi * (pi * sinc_x + 2 * sinc_x_d / x_pi)
    return sinc_x, sinc_x_d, sinc_x_h
end

# A function without checking a.  Factored out for testing with ForwardDiff..
function lanczos_kernel_with_derivatives_nocheck{NumType <: Number}(
    x::NumType, a::Float64)
    sinc_x, sinc_x_d, sinc_x_h = sinc_with_derivatives(x)
    sinc_xa, sinc_xa_d, sinc_xa_h = sinc_with_derivatives(x / a)

    return sinc_x * sinc_xa,
           sinc_x_d * sinc_xa + sinc_x * sinc_xa_d / a,
           sinc_x_h * sinc_xa + 2 * sinc_x_d * sinc_xa_d / a +
              sinc_x * sinc_xa_h / (a ^ 2)
end

function lanczos_kernel_with_derivatives{NumType <: Number}(x::NumType, a::Float64)
    if abs(x) > a
        return 0, 0, 0
    end
    return lanczos_kernel_with_derivatives_nocheck(x, a)
end

function bspline_kernel_with_derivatives{NumType <: Number}(x::NumType)
    abs_x = abs(x)
    sign_x = sign(x)
    if abs_x > 2
        return 0, 0, 0
    elseif abs_x > 1
        return (-1 * abs_x^3 + 6 * abs_x^2 - 12 * abs_x + 8) / 6,
               (-3 * abs_x^2 + 12 * abs_x  - 12) * sign_x / 6,
               (-6 * abs_x + 12) / 6
    else
        return (3 * abs_x^3 - 6 * abs_x^2 + 4) / 6,
               (9 * abs_x^2 - 12 * abs_x) * sign_x / 6,
               (18 * abs_x  - 12) / 6
    end
end

function cubic_kernel_with_derivatives{NumType <: Number}(x::NumType, a::Float64)
    abs_x = abs(x)
    sign_x = sign(x)
    if abs_x > 2
        return 0, 0, 0
    elseif abs_x > 1
        return a * abs_x^3      -5 * a * abs_x^2 +  8 * a * abs_x - 4 * a,
               (3 * a * abs_x^2  -10 * a * abs_x +   8 * a) * sign_x,
               6 * a * abs_x    -10 * a
    else
        return (a + 2) * abs_x^3 -    (a + 3) * abs_x^2 + 1,
               (3 * (a + 2) * abs_x^2  -2 * (a + 3) * abs_x) * sign_x,
               6 * (a + 2) * abs_x    -2 * (a + 3)
    end
    return x
end

# TODO: replace the gradient and hessian flags with whether or not the
# sensitive floats in the image have gradients and hessians.
# Interpolate the PSF to the pixel values.
function interpolate!{NumType <: Number}(
        kernel_with_derivs, # A 1-d kernel function
        kernel_width::Int,
        image::Matrix{SensitiveFloat{NumType}},
        psf_image::Matrix{Float64},
        object_loc::Vector{NumType},
        wcs_jacobian::Matrix{Float64},
        calculate_gradient::Bool,
        calculate_hessian::Bool)
    @assert image[1].local_P == length(StarPosParams)

    h_psf_width = (size(psf_image, 1) + 1) / 2.0
    w_psf_width = (size(psf_image, 2) + 1) / 2.0

    param_ids = getids(StarPosParams)

    # These are sensitive floats representing derviatives of the kernel.
    kernel = SensitiveFloat{NumType}(length(StarPosParams), 1,
                                     calculate_gradient, calculate_hessian)
    kernel_h = SensitiveFloat{NumType}(length(StarPosParams), 1,
                                       calculate_gradient, calculate_hessian)

    # Pre-compute terms for transforming derivatives to world coordinates.
    if calculate_gradient
        k_h_grad = wcs_jacobian' * Float64[1, 0]
        k_w_grad = wcs_jacobian' * Float64[0, 1]

        if calculate_hessian
            k_h_hess = wcs_jacobian' * Float64[1 0; 0 0] * wcs_jacobian
            k_w_hess = wcs_jacobian' * Float64[0 0; 0 1] * wcs_jacobian
        end
    end

    # h, w are pixel coordinates.
    for h = 1:size(image, 1), w = 1:size(image, 2)

        # h_psf, w_psf are in psf coordinates.
        # The PSF is centered at object_loc + psf_width.
        h_psf = h - object_loc[1] + h_psf_width
        w_psf = w - object_loc[2] + w_psf_width

        # Centers of indices of the psf matrix, i.e., integer psf coordinates.
        h_ind0, w_ind0 = Int(floor(h_psf)), Int(floor(w_psf))
        h_lower = max(h_ind0 - kernel_width + 1, 1)
        h_upper = min(h_ind0 + kernel_width, size(psf_image, 1))
        for h_ind = (h_lower:h_upper)
            lh_v, lh_d, lh_h = kernel_with_derivs(h_psf - h_ind)
            if lh_v != 0
                clear!(kernel_h)
                kernel_h.v[] = lh_v

                if calculate_gradient
                    # This is -1 * wcs_jacobian' * [lh_d, 0]
                    # and -1 * wcs_jacobian' * [lh_h 0; 0 0] * wcs_jacobian
                    kernel_h.d[param_ids.u] = -1 * k_h_grad * lh_d
                    if calculate_hessian
                        kernel_h.h[param_ids.u, param_ids.u] = k_h_hess * lh_h;
                    end
                end
                w_lower = max(w_ind0 - kernel_width + 1, 1)
                w_upper = min(w_ind0 + kernel_width, size(psf_image, 2))
                for w_ind = (w_lower:w_upper)
                    lw_v, lw_d, lw_h = kernel_with_derivs(w_psf - w_ind)
                    if lw_v != 0
                        clear!(kernel)
                        kernel.v[] = lw_v
                        # This is -1 * wcs_jacobian' * [0, lw_d]
                        # and -1 * wcs_jacobian' * [0 0; 0 lw_h] * wcs_jacobian
                        if calculate_gradient
                            kernel.d[param_ids.u] = -1 * k_w_grad * lw_d;
                            if calculate_hessian
                                kernel.h[param_ids.u, param_ids.u] = k_w_hess * lw_h;
                            end
                            multiply_sfs!(kernel, kernel_h)
                            add_scaled_sfs!(
                                image[h, w], kernel, psf_image[h_ind, w_ind])
                        else
                            image[h, w].v[] +=
                                kernel.v[] * kernel_h.v[] * psf_image[h_ind, w_ind]
                        end
                    end
                end
            end
        end
    end
end


text:
include("AdaptiveFilter.jl")

"""
This file contains the code for all the plots in the result section
of the thesis in the order in which they appear.
(`latencies` and `alpha_lambda` are not part of the thesis.)
For running, uncomment the desired method at the bottom.
"""

function pretest_sin()
    global training_cycles = 5000
    global task = "sin"
    global show = "out"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6

    ns = 2:2:16
    function mses_n(algo_code)
        global algo = algo_code

        mses = []
        for n_ in ns
            global n = n_
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_n(:af)
    la_mses = mses_n(:lazar)
    n_plot = plot(ns, [af_mses, la_mses],
                  xlabel="n", ylabel="mse", xticks=ns,
                  labels=["AF", "Lazar"])
    scatter!([10, 10], [af_mses[5], la_mses[5]], lab="n=10")
    global n = 10 # reset

    rates = 2:2:16
    function mses_rate(algo_code)
        global show = "out" # dummy, only save end
        global algo = algo_code

        mses = []
        for rate_ in rates
            global rate = rate_
            global κ = (b+offset)/rate/θ   # capacitance (integrator weight)
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_rate(:af)
    la_mses = mses_rate(:lazar)
    rate_plot = plot(rates, [af_mses, la_mses],
                     xlabel="rate", ylabel="mse", xticks=rates,
                     labels=["AF", "Lazar"])
    scatter!([10, 10], [af_mses[5], la_mses[5]], lab="rate=10")
    # reset
    global rate = 10
    global κ = (b+offset)/rate/θ   # capacitance (integrator weight)

    latencies = 0:.3:3
    function mses_lat(algo_code)
        global show = "out" # dummy, only save end
        global algo = algo_code

        mses = []
        for latency_ in latencies
            global latency = - latency_ # mind the minus!
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_lat(:af)
    la_mses = mses_lat(:lazar)
    lat_plot = plot(latencies, [af_mses, la_mses],
                    xlabel="latency", ylabel="mse", xticks=0:.6:3,
    labels=["AF", "Lazar"])
    scatter!([.6, .6], [af_mses[3], la_mses[3]], lab="latency=0.6")
    global latency = -.6 # reset

    alphas = .25:.25:3.25
    function mses_algo(algo_code)
        global algo = algo_code

        mses = []
        for alpha_ in alphas
            global α = alpha_
            global λ = alpha_
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_algo(:af)
    la_mses = mses_algo(:lazar)
    alpha_plot = plot(alphas, [af_mses, la_mses],
                      xlabel="alpha, lambda", ylabel="mse", xticks=.5:.5:3,
                      labels=["AF", "Lazar"])
    scatter!([2, 2], [af_mses[8], la_mses[8]], lab="alpha=lambda=2")
    # reset
    global α = 2
    global λ = 2

    plot(n_plot, rate_plot, lat_plot, alpha_plot,
         title=["A" "B" "C" "D"], titleloc=:left,
         layout=(2,2), size=(500, 500))
    savefig("../../Thesis/figures/pretest-alpha$α-latency$latency.pdf")
end

function latencies()
    global training_cycles = 5000
    global task = "sin"
    global show = "out"
    global α = 2
    global λ = 2
    global η = 1e-4
    latencies = 0:.3:3.
    function mses_lat(algo_code)
        global algo = algo_code

        mses = []
        for latency_ in latencies
            global latency = - latency_ # mind the minus!
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_lat(:af)
    la_mses = mses_lat(:lazar)
    n_plot = plot(latencies, [af_mses, la_mses],
                  xlabel="latency", ylabel="mse", xticks=0:.6:3,
                  labels=["AF", "Lazar"])
    scatter!([.6, .6], [af_mses[3], la_mses[3]], lab="latency=.6")
    savefig("../../Thesis/figures/latencies.pdf")
end

function alpha_lambda()
    global training_cycles = 5000
    global task = "sin"
    global show = "out"
    global latency = -.6
    global η = 1e-4
    alphas = .25:.25:3.5
    function mses_algo(algo_code)
        global algo = algo_code

        mses = []
        for alpha_ in alphas
            global α = alpha_
            global λ = alpha_
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol      = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_algo(:af)
    la_mses = mses_algo(:lazar)
    alpha_plot = plot(alphas, [af_mses, la_mses],
                      xlabel="alpha, lambda", ylabel="mse", xticks=alphas,
                      labels=["AF", "Lazar"],
                      size=(300, 200))
    savefig("../../Thesis/figures/alpha-lambda-latency$latency.pdf")
end

function training_error_sin()
    global training_cycles = 10000
    global task = "sin"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6
    function plot_algo(algo_code, algo_name, point_names)
        global algo = algo_code
        global show = "mse_evol" # to save all periods
        io, tp, u0, tsave = init_training(task, training_cycles)
        sol    = run(io, u0, tsave)

        p1, p2 = (2500, 10000)
        pn1, pn2 = point_names
        mses, sols = running_mse(sol, io, [p1, p2])

        al_evol = plot(mses, title="$algo_name", lab="MSE",
                       titleloc=:left, leg=false,
                       ylims=(0,.8),
                       ylabel="mse", xlabel="periods")
        scatter!([p1, p2], [mses[p1], mses[p2]])
                 #series_annotations = [pn1, pn2])
        annotate!([(p1, mses[p1], text(pn1, :bottom)),
                   (p2, mses[p2], text(pn2, :bottom))])
        al_out1 = plot([t->out(sols[1], t), io.z], sols[1].t,
                       xticks=sols[1].t[1]:tp:sols[1].t[end],
                       title="$pn1 (mse: $(mses[p1]))",
                       titleloc=:left, leg=false)
        al_out2 = plot([t->out(sols[2], t), io.z], sols[2].t,
                       xticks=sols[2].t[1]:tp:sols[2].t[end],
                       title="$pn2 (mse: $(mses[p2]))",
                       titleloc=:left, leg=false,  xlabel="time")

        l = @layout [ a{.5w} [b{.5h}
                              c{.5h}] ]
        plot(al_evol, al_out1, al_out2, layout=l)
    end
    af_plot = plot_algo(:af, "AF", ("A", "B"))
    la_plot = plot_algo(:lazar, "Lazar", ("C", "D"))
    plot(af_plot, la_plot, layout=(2,1), size=(500, 500))
    savefig("../../Thesis/figures/sin-training-cycles20000.pdf")
end

function weights_taps_sin()
    global training_cycles = 10000
    global task = "sin"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6
    function plot_algo(algo_code, algo_name)
        global show = "w_evol"
        global algo = algo_code
        io, tp, u0, tsave = init_training(task, training_cycles)
        sol    = run(io, u0, tsave)
        w_plot = evol_plot("w_evol", sol, io)
        # sol can be used to plot y_end as well
        y_plot = end_plot("y", sol, io, tp)
        # stack w and y horizontally
        return (w_plot, y_plot)
    end
    # stack af and lazar vertically
    af_w, af_y = plot_algo(:af, "AF")
    la_w, la_y = plot_algo(:lazar, "Lazar")
    plot(af_w, af_y, la_w, la_y, layout=(2,2), leg=false)
    savefig("../../Thesis/figures/sin-weights-taps.pdf")
end

function imp_sin()
    global training_cycles = 10000
    global task = "sin"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6

    global algo = :af
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    af_imp = impulse_response(sol)

    global algo = :lazar
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    la_imp = impulse_response(sol)

    af_impi = plot([y_i(af_imp, i) for i in 1:n],
                   af_imp.t, leg=false,
                   labels=["y$i" for i in 1:n],
                   title="AF, individual IR",
                   titleloc=:left, legend=:bottomright)
    la_impi = plot([y_i(la_imp, i) for i in 1:n],
                   la_imp.t,
                   labels=["y$i" for i in 1:n],
                   title="Lazar, individual IR",
                   titleloc=:left, legend=:bottomright)
    both = plot([t->out(af_imp, t), t->out(la_imp, t)],
                af_imp.t,
                label=["AF", "Lazar"], xlabel="time",
                title="Total IR",
                titleloc=:left, legend=:topright)

    plot(af_impi, la_impi, both, layout=(3,1), size=(400,600))

    savefig("../../Thesis/figures/sin-imp.pdf")
end

function task_noise()
    global training_cycles = 10000
    global task = "noise"
    global show = "out"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6
    nstds = 0:.3:1.8
    function mses_algo(algo_code)
        global algo = algo_code

        mses = []
        for nstd in nstds
            global noise_std = nstd
            io, tp, u0, tsave = init_training(task, training_cycles)
            sol    = run(io, u0, tsave)
            _, error = test_error(sol, io, tp)
            push!(mses, error)
        end
        return mses
    end
    af_mses = mses_algo(:af)
    la_mses = mses_algo(:lazar)
    main = plot(nstds, [af_mses, la_mses],
                title="Noise task", titleloc=:left,
                xlabel="noise std", ylabel="mse",
                labels=["AF", "Lazar"], leg=:topleft)
    # again for subplots
    global noise_std = 1.2
    global algo = :lazar
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    la_sol, error = test_error(sol, io, tp)
    global algo = :af
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    af_sol, error = test_error(sol, io, tp)

    sub1 = plot([io.x, io.z], la_sol.t, labels=["input", "target"],
                title="for std = 1.2", titleloc=:left)
    sub2 = plot([t->out(la_sol, t), io.z], la_sol.t,
                title="Lazar", titleloc=:left, leg=false)
    sub3 = plot([t->out(af_sol, t), io.z], af_sol.t,
                title="AF", titleloc=:left, leg=false)
    l = @layout [ a{.4w} [b{.33h}
                          c{.33h}
                          d{.33h}] ]
    plot(main, sub1, sub2, sub3, layout=l)
    savefig("../../Thesis/figures/noise-std-sub.pdf")
end

function task_mix3()
    global training_cycles = 25000
    global test_cycles = 1 # because 1 cycle is already very long
    global task = "mix3"
    global show = "out"
    global α = 3 # mix3 specific
    global λ = 3
    global η = 1e-3
    global latency = -.6
    io, _ = init_task("1of3") # NOTE important for io in plot
    function out_algo(algo_code)
        global algo = algo_code
        global mode = "train"
        io, tp, u0, tsave = init_training(task, training_cycles)
        sol    = run(io, u0, tsave)
        sol2, error = test_error(sol, io, tp)
        return sol2, error
    end
    sol_af, mse_af = out_algo(:af)
    sol_la, mse_la = out_algo(:lazar)
    plot([t->io.z(t), t->out(sol_af, t), t->out(sol_la, t)], sol_af.t,
         labels=["target", "AF (MSE: $mse_af)", "Lazar (MSE: $mse_la)"],
         title="mix3", xlabel="time")
    savefig("../../Thesis/figures/mix3-task.pdf")
    global test_cycles = 3 # reset
end

function task_1of3()
    global training_cycles = 25000
    global task = "1of3"
    global show = "out"
    global α = 3
    global λ = 3
    global η = 1e-3
    global latency = -1 # tp/4 = 4/4 = 1
    io, _ = init_task("1of3") # NOTE important for io in plot
    function out_algo(algo_code)
        global algo = algo_code
        io, tp, u0, tsave = init_training(task, training_cycles)
        sol         = run(io, u0, tsave)
        sol2, error = test_error(sol, io, tp)
        return sol2, error
    end
    sol_af, mse_af = out_algo(:af)
    sol_la, mse_la = out_algo(:lazar)
    plot([io.x, io.z, t->out(sol_af, t), t->out(sol_la, t)], sol_af.t,
         labels=["input", "target",
                 "AF (MSE: $mse_af)", "Lazar (MSE: $mse_la)"],
         title="1of3", xlabel="time")
    savefig("../../Thesis/figures/1of3-task.pdf")
end

function imp_1of3()
    global training_cycles = 25000
    global task = "1of3"
    global show = "out"
    global α = 3
    global λ = 3
    global latency = -1 # tp/4 = 4/4 = 1
    global η = 1e-3

    global algo = :af
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    af_imp = impulse_response(sol)

    global algo = :lazar
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    la_imp = impulse_response(sol)

    plot([t->out(af_imp, t), t->out(la_imp, t)],
         af_imp.t,
         label=["AF", "Lazar"], xlabel="time",
         title="Total impulse response",
         titleloc=:left, legend=:topright)
    savefig("../../Thesis/figures/imp-1of3.pdf")
end

function imp_mix3()
    global training_cycles = 1000
    global task = "mix3"
    global show = "out"
    global α = 3
    global λ = 3
    global latency = -.6 # tp/4 = 4/4 = 1
    global η = 1e-3

    global algo = :af
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    af_imp = impulse_response(sol)

    global algo = :lazar
    io, tp, u0, tsave = init_training(task, training_cycles)
    sol    = run(io, u0, tsave)
    la_imp = impulse_response(sol)

    plot([t->out(af_imp, t), t->out(la_imp, t)],
         af_imp.t,
         label=["AF", "Lazar"], xlabel="time",
         title="Total impulse response",
         titleloc=:left, legend=:topright)
    savefig("../../Thesis/figures/imp-mix3.pdf")
end

function task_relu()
    global training_cycles = 25000
    global task = "relu"
    global show = "out"
    global α = 2
    global λ = 2
    global η = 1e-4
    global latency = -.6
    #global offset = 1
    io, _ = init_task("relu") # NOTE important for io in plot
    function out_algo(algo_code)
        global algo = algo_code
        io, tp, u0, tsave = init_training(task, training_cycles)
        sol    = run(io, u0, tsave)
        sol2, error = test_error(sol, io, tp)
        return sol2, error
    end
    sol_af, mse_af = out_algo(:af)
    sol_naf, mse_naf = out_algo(:naf)
    sol_la, mse_la = out_algo(:lazar)
    plot([io.z, t->out(sol_af, t), t->r(sol_naf, t), # mind the r (ramp)
          t->out(sol_la, t)], sol_af.t,
         labels=["target", "AF (MSE: $mse_af)", "nAF (MSE: $mse_naf)",
                 "Lazar (MSE: $mse_la)"],
         title="reLU", xlabel="time")
    savefig("../../Thesis/figures/relu-task.pdf")
end

# Pretests
#pretest_sin()
#latencies()
#alpha_lambda()

# Sine recovery
#training_error_sin()
#weights_taps_sin()
#imp_sin()

# Other tasks
#task_noise()
#task_mix3()
#task_1of3()
task_relu()

#imp_1of3()
#imp_mix3()

max_stars_repo_name: toludwig/AdaptiveFilterNeuron

max_issues_repo_name: toludwig/AdaptiveFilterNeuron

max_forks_repo_name: toludwig/AdaptiveFilterNeuron


text:
using Yao, YaoExtensions
using Test, YaoArrayRegister, LinearAlgebra

@testset "test reflect gate" begin
    reg0 = rand_state(3)
    mirror = randn(1<<3)*im; mirror[:] /= norm(mirror)
    rf = YaoExtensions.reflect(mirror)
    reg = copy(reg0)
    apply!(reg, rf)

    v0, v1 = vec(reg.state), vec(reg0.state)
    @test mat(rf) * statevec(reg0) ≈ statevec(apply!(copy(reg0), rf))
    @test statevec(rf.psi)' * v0 ≈ statevec(rf.psi)' * v1
    @test v0 - statevec(rf.psi)'*v0* statevec(rf.psi) ≈
        -(v1- statevec(rf.psi)' * v1 * statevec(rf.psi))

    @test mat(ComplexF32, rf) == ComplexF32.(mat(rf))

    @test rf == YaoExtensions.reflect(rf.psi)
    # copy do not occur on register
    @test copy(rf).psi === rf.psi

    @test ishermitian(rf)
    @test isreflexive(rf)
    @test isunitary(rf)
end


text:
function density_biased_sample(data, labels, num_samples=25; from_oversample=false, kde_bw="scott")
    indices = collect(1:length(labels))
    inlier_indices = indices[labels .== :inlier]
    d = data[:, inlier_indices]
    p = multi_kde(d, kde_bw)(d)
    ranking = sortperm(p)
    if from_oversample
        selection = ranking[1:min(end, num_samples * 2)]
        selection = selection[randperm(length(selection))[1:min(end, num_samples)]]
    else
        selection = ranking[1:num_samples]
    end
    return inlier_indices[selection]
end

import NearestNeighbors

function nn_biased_sample(data, labels, num_samples=25; kwargs...)
    indices = collect(1:length(labels))
    inlier_indices = indices[labels .== :inlier]
    d = data[:, inlier_indices]
    kdtree = NearestNeighbors.KDTree(d)
    neighbors = NearestNeighbors.knn(kdtree, d, length(inlier_indices), true)[1]

    selection = Set([rand(1:length(inlier_indices))])
    last = first(selection)
    for _ in 1:num_samples-1
        for i in 2:length(inlier_indices)
            candidate = neighbors[last][i]
            if candidate ∉ selection
                push!(selection, candidate)
                last = candidate
                break
            end
        end
    end
    selection = collect(selection)
    return inlier_indices[selection]
end

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
# low-level C version, use CCVDenseMatrix instead
@runonce immutable CCVDenseMatrixType{T}
    typ::Cint
    sig::UInt64
    refcount::Cint
    rows::Cint
    cols::Cint
    step::Cint
    tag::Float64 # can actually be anything, but mapping to C
                 #  requires this field to be 64 bits
    data::Ptr{T}
end

@runonce type CCVDenseMatrix{T} <: AbstractArray{T,2}
    ptr::Ptr{CCVDenseMatrixType{T}}
    function CCVDenseMatrix(ptr::Ptr{CCVDenseMatrixType{T}})
        mat = new(ptr)
        finalizer(mat, mat -> free(mat))
        mat
    end
end

function CCVDenseMatrix(rows::Int, cols::Int, typ::UInt32)
    T = julia_type(typ)
    ptr = ccall((:ccv_dense_matrix_new, libccv), Ptr{CCVDenseMatrixType{T}},
              (Cint, Cint, Cint, Ptr{Void}, UInt64),
              Cint(rows), Cint(cols), Cint(typ), C_NULL, UInt64(0))
    return CCVDenseMatrix(ptr)
end

function free(mat::CCVDenseMatrix)
    if mat.ptr != C_NULL
        ccall((:ccv_matrix_free, libccv), Void,
              (Ptr{CCVDenseMatrixType},), mat.ptr)
        mat.ptr = C_NULL
    end
end

CCVDenseMatrix(rows::Int, cols::Int) = CCVDenseMatrix(rows, cols, CCV_64F)

function Base.show(io::IO, mat::CCVDenseMatrix)
    struct = unsafe_load(mat.ptr)
    print(io, "CCVDenseMatrix($(struct.rows),$(struct.cols))")
end

# for REPL
Base.writemime(io::IO, ::MIME"text/plain", mat::CCVDenseMatrix) = show(io, mat)

Base.linearindexing(::Type{CCVDenseMatrix}) = Base.LinearFast()

function Base.getindex(mat::CCVDenseMatrix, i::Integer)
    struct = unsafe_load(mat.ptr)
    offset = i
    if i < 1 || i > struct.rows * struct.cols
        throw(BoundsError("Index $(i) is outside of " *
                          "matrix length $(struct.rows * struct.cols)"))
    end
    return unsafe_load(struct.data, offset)
end

function Base.getindex(mat::CCVDenseMatrix, i::Integer, j::Integer)
    struct = unsafe_load(mat.ptr)
    offset = struct.cols * (i-1) + j
    if i < 1 || i > struct.rows || j < 1 || j > struct.cols
        throw(BoundsError("Index $((i, j)) is outside of " *
                          "matrix size $((struct.rows, struct.cols))"))
    end
    return unsafe_load(struct.data, offset)
end

function Base.setindex!{T}(mat::CCVDenseMatrix{T}, v::T, i::Integer, j::Integer)
    struct = unsafe_load(mat.ptr)
    offset = struct.cols * (i-1) + j
    if i < 1 || i > struct.rows || j < 1 || j > struct.cols
        throw(BoundsError("Index $((i, j)) is outside of " *
                          "matrix size $((struct.rows, struct.cols))"))
    end
    return unsafe_store!(struct.data, v, offset)
end

function Base.setindex!{T}(mat::CCVDenseMatrix{T}, v::T, i::Integer)
    struct = unsafe_load(mat.ptr)
    offset = i
    if i < 1 || i > struct.rows * struct.cols
        throw(BoundsError("Index $((i, j)) is outside of " *
                          "matrix size $((struct.rows, struct.cols))"))
    end
    return unsafe_store!(struct.data, v, offset)
end

function Base.size(mat::CCVDenseMatrix)
    struct = unsafe_load(mat.ptr)
    return (Int(struct.rows), Int(struct.cols))
end

function matrix_free(mat::CCVDenseMatrix)
    ccall((:ccv_matrix_free, libccv), Void, (Ptr{CCVDenseMatrix},), pointer(CCVDenseMatrix[mat]))
end

"Low-level function to create CCVDenseMatrix from Julia 1D array"
function from_vector{T}(data::Vector{T}, rows::Int, cols::Int)
    @assert rows*cols == length(data) "rows * cols != length of data!"
    matout = Ref{Ptr{CCVDenseMatrixType{T}}}()
    scanline = sizeof(T)
    typ = ccv_type(T) | CCV_IO_ANY_RAW
    code = ccall((:ccv_read_impl, libccv), Cint,
                 (Ptr{Void}, Ptr{Ptr{CCVDenseMatrixType{T}}}, Cint, Cint, Cint, Cint),
                 data, matout, typ, rows, cols, scanline)
    if code != CCV_IO_FINAL
        throw(ErrorException("ccv_read_impl() exited with error: $code"))
    else
        ptr = matout[]
        return CCVDenseMatrix{T}(ptr)
    end
end

function Base.convert{T}(::Type{CCVDenseMatrix}, arr::Array{T,2})
    tarr = arr'  # convert underlying memory from column-major to row-major layout
    rows, cols = size(arr)
    return from_vector(reshape(tarr, rows*cols), rows, cols)
end

Base.cconvert{T}(::Type{Ptr{CCVDenseMatrixType{T}}}, mat::CCVDenseMatrix{T}) = mat
Base.unsafe_convert{T}(::Type{Ptr{CCVDenseMatrixType{T}}}, mat::CCVDenseMatrix{T}) = mat.ptr


text:
##################################################################
# Source code taken from: https://github.com/jbytecode/MCCGA
# Credits to corresponding author.
##################################################################
include("utils.jl")

mutable struct MCCGA <: AbstractParameters
    N::Int # population size
    maxsamples::Int
    mutation::Float64
    probvector::Vector{Float64}
    use_local_search::Bool
end

"""
    MCCGA(;N, maxsamples)

### Parameters:

- `N` population size. Default is 100.
- `maxsamples` maximum number of samples. Default is 10000.

### Description

MCCGA method implements the Machine-coded genetic algorithms for real valued optimization problems.
The algorithm is based on the concept of a compact genetic algorithm but with a machine-coded representation
using IEEE-754 floating point encoding standard. In the first stage of the algorithm, `maxsamples` number
of samples are generated within the range of function domain. This process is required to obtain a vector
of probabilities for each single bit of the IEEE-754 representation. In classical CGAs, the initial vector
of probabilities is generated using the constant probability of 0.5 whereas in MCCGA, the probability of ith bit having
a value of 1 depends on the function domain. The second step performs a CGA search but with IEEE-754 bits again. Since
CGA does not use a population of solutions but a single vector of probabilities, the parameter `N` does not really mean
number of solutions. Instead, it means the amount of mutation in each iteration, e.g. 1/N.
In the last stage, a local search is performed for fine-tuning. In this implementation, `Optim` package is required.

### References

- Harik, G. R., Lobo, F. G., & Goldberg, D. E. (1999). The compact genetic algorithm. IEEE transactions on evolutionary computation, 3(4), 287-297.
- Satman, M. H. & Akadal, E. (2020). Machine Coded Compact Genetic Algorithms for Real Parameter Optimization Problems . Alphanumeric Journal , 8 (1) , 43-58 . DOI: 10.17093/alphanumeric.576919
- Mehmet Hakan Satman, Emre Akadal, Makine Kodlu Hibrit Kompakt Genetik Algoritmalar Optimizasyon Yöntemi, Patent, TR, 2022/01, 2018-GE-510239

### Example

```jldoctest

julia> import Optim

julia> f, bounds, solutions = Metaheuristics.TestProblems.rastrigin();

julia> result = optimize(f, bounds, MCCGA())

+=========== RESULT ==========+
  iteration: 42833
    minimum: 0
  minimizer: [-5.677669786379456e-17, 2.7942451898022582e-39, -3.60925916059986e-33, -6.609510861086017e-34, 2.998586759675011e-32, -1.8825832500775007e-38, -3.0729484147585938e-31, 1.7675578057632446e-38, 5.127944874215823e-16, 1.9623770480857484e-19]
    f calls: 86404
 total time: 3.2839 s
+============================+
```
"""
function MCCGA(;
        N = 100,
        maxsamples =10_000,
        mutation = 1 / N,
        use_local_search = true,
        information = Information(),
        options = Options()
    )

    parameters = MCCGA(N, maxsamples, mutation, [], use_local_search)

    Algorithm(
        parameters,
        information = information,
        options = options,
    )
end

function initialize!(
        status,
        parameters::MCCGA,
        problem,
        information,
        options,
        args...;
        kargs...
    )

    if options.iterations == 0
        options.iterations = 100_000_000
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = options.iterations * parameters.N + 1
    end

    lower = problem.bounds[1,:]
    upper = problem.bounds[2,:]

    parameters.probvector = initialprobs(lower, upper, maxsamples = parameters.maxsamples)

    # sample a vector to create an initial State
    x = sample(parameters.probvector) |> floats
    initial_sol = create_solution(x, problem)
    return State(initial_sol, [initial_sol for i in 1:parameters.N])

end

function update_state!(
        status,
        parameters::MCCGA,
        problem,
        information,
        options,
        args...;
        kargs...
    )

    probvector = parameters.probvector
    chsize = length(probvector)

    # parents
    ch1 = sample(probvector)
    ch2 = sample(probvector)

    # evaluate cost function
    sol1 = create_solution(floats(ch1), problem)
    sol2 = create_solution(floats(ch2), problem)

    sol_winner = sol1
    winner = ch1
    loser  = ch2

    # check if ch2 (sol2) is better that ch1 (sol1)
    if is_better(sol2, sol1)
        sol_winner = sol2
        winner = ch2
        loser = ch1
    end

    # save in population the winner (informative only)
    # this is not used in the algorithm
    status.population[1 + status.iteration % parameters.N] = sol_winner

    # save best solution found so far
    is_better(sol_winner, status.best_sol) &&  (status.best_sol = sol_winner)

    for i = 1:chsize
        if winner[i] != loser[i]
            if winner[i] == 1
                probvector[i] += parameters.mutation
            else
                probvector[i] -= parameters.mutation
            end
            if probvector[i] > 1
                probvector[i] = 1
            elseif probvector[i] < 0
                probvector[i] = 0
            end
        end
    end
end

function final_stage!(
        status::State{xf_indiv}, # unconstrained case
        parameters::MCCGA,
        problem,
        information,
        options,
        args...;
        kargs...
    )

    if !parameters.use_local_search
        status.final_time = time()
        return
    end

    costfunction(x) = evaluate(x, problem)

    probvector = parameters.probvector
    sampledvector = sample(probvector)

    # initial solution for the local search
    initial_solution = floats(sampledvector)

    options.debug && @info "Running NelderMead..."
    local_result = Optim.optimize(costfunction, initial_solution, Optim.NelderMead())
    # display Nelder-Mead result
    options.debug && display(local_result)
    options.debug && @info "NelderMead done!"

    # save best solution found so far!
    sol = create_child(local_result.minimizer, local_result.minimum)

    status.final_time = time()
    status.f_calls = problem.f_calls

    if is_better(sol, status.best_sol)
        status.best_sol = sol
    end

    return
end

function final_stage!(
        status,
        parameters::MCCGA,
        problem,
        information,
        options,
        args...;
        kargs...
    )

    status.final_time = time()
    # nothing to do for constrained or multi-objective
    @warn "MCCGA has been designed for unconstrained problems (local search not performed)."
end

function stop_criteria!(status, parameters::MCCGA, problem, information, options)
    # check budget limitation
    if status.stop
        return
    end

    mutation = parameters.mutation
    status.stop = all(x -> (x <= mutation) || (x >= 1.0 - mutation), parameters.probvector)

    status.stop && options.debug && @info "MCCGA stopping criteria is met."
end


text:
using FileIO
using Test
using UUIDs

const testdata_filename = "deleteme.klb"
const testdata_dtype = UInt16
const testdata_size = [521, 521, 64]
const testdata_lower_bounds = UInt32[1, 1, 1, 1, 1]
const testdata_upper_bounds = UInt32[64, 64, 32, 1, 1]
const testdata_bounds = cat(testdata_lower_bounds, testdata_upper_bounds; dims=2)

add_format(format"KLB", (), ".klb", [:KLB => UUID("8bb66c0d-974d-412d-8f46-f9b8d1ef37d0")])

save(testdata_filename, rand(testdata_dtype, testdata_size...));
@test isfile(testdata_filename) == true

A = load(testdata_filename)
@test eltype(A) == testdata_dtype
@test prod(size(A)) == prod(testdata_size)
@test sum(A) != 0

A = zeros(testdata_dtype, testdata_size...)
load(testdata_filename, inplace=A, nochecks=true)
@test sum(A) != 0

A = load(testdata_filename, bounds=testdata_bounds)
@test eltype(A) == testdata_dtype
@test prod(size(A)) == prod(testdata_upper_bounds)
@test sum(A) != 0

A = zeros(testdata_dtype, testdata_upper_bounds...)
load(testdata_filename, bounds=testdata_bounds, inplace=A, nochecks=true)
@test sum(A) != 0

using KLB

header = klbheader(testdata_filename)
@test header["datatype"] == testdata_dtype
@test prod(header["imagesize"]) == prod(testdata_size)

A = KLB.loadklb(testdata_filename)
@test eltype(A) == testdata_dtype
@test prod(size(A)) == prod(testdata_size)
@test sum(A) != 0

A = zeros(testdata_dtype, testdata_size...)
KLB.loadklb!(A, testdata_filename)
@test sum(A) != 0

A = KLB.loadklb(testdata_filename, testdata_lower_bounds, testdata_upper_bounds)
@test eltype(A) == testdata_dtype
@test prod(size(A)) == prod(testdata_upper_bounds)
@test sum(A) != 0

A = zeros(testdata_dtype, testdata_upper_bounds...)
KLB.loadklb!(A, testdata_filename, testdata_lower_bounds, testdata_upper_bounds)
@test sum(A) != 0

A = zeros(testdata_dtype, testdata_upper_bounds...)
KLB.loadklb!(A, testdata_filename, (testdata_lower_bounds...,), (testdata_upper_bounds...,))
@test sum(A) != 0

rm(testdata_filename)
@test isfile(testdata_filename) == false


text:
using LearnConvection
using Plots

default_modify_predictor_fn(x, 𝒟, time_index) = x

modify_pred_fns = [
    default_modify_predictor_fn,
    append_tke,
    partial_temp_profile(1:16),
    partial_temp_profile(17:32),
]

f = append_tke

problems = [
            Slack("KPP"; parameters=KPP.Parameters(), modify_predictor_fn=f),
            Slack("TKE"; parameters=TKEMassFlux.TKEParameters(), modify_predictor_fn=f),
            Residual("KPP"; parameters=KPP.Parameters(), modify_predictor_fn=f),
            Residual("TKE"; parameters=TKEMassFlux.TKEParameters(), modify_predictor_fn=f),
            Sequential("TKE"; parameters=TKEMassFlux.TKEParameters(), modify_predictor_fn=f),
            Sequential("KPP"; parameters=KPP.Parameters(), modify_predictor_fn=f),
            Sequential("T"; modify_predictor_fn=f),
            Sequential("dT"; modify_predictor_fn=f),
]

## Interpolation

train = ["general_strat_4_profiles.jld2", "general_strat_32_profiles.jld2"]
validate = ["general_strat_12_profiles.jld2", "general_strat_24_profiles.jld2"]
test = ["general_strat_8_profiles.jld2", "general_strat_16_profiles.jld2", "general_strat_20_profiles.jld2", "general_strat_28_profiles.jld2"]

## Extrapolation

train = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2", "general_strat_12_profiles.jld2", "general_strat_16_profiles.jld2"]
validate = ["general_strat_20_profiles.jld2", "general_strat_24_profiles.jld2"]
test = ["general_strat_28_profiles.jld2", "general_strat_32_profiles.jld2"]

descriptor = "extrapolation_append_tke"
D=32
N=4

for problem in problems

    println("--*--*--*--*--*--$(problem)--*--*--*--*--*--")

    𝒟_train     = LearnConvection.Data.data(train, problem; D=D, N=N);
    𝒟_validate  = LearnConvection.Data.data(validate, problem; D=D, N=N);
    𝒟_test      = LearnConvection.Data.data(test, problem; D=D, N=N);

    train_validate_test(𝒟_train, 𝒟_validate, 𝒟_test, problem; log_γs=-2.0:0.1:2.0, distances=[euclidean_distance], descriptor=descriptor)

end

# problems = [
#             Residual("KPP"; parameters=KPP.Parameters(), modify_predictor_fn=f, convective_adjust=true),
#             Residual("TKE"; parameters=TKEMassFlux.TKEParameters(), modify_predictor_fn=f, convective_adjust=true),
#             Sequential("TKE"; parameters=TKEMassFlux.TKEParameters(), modify_predictor_fn=f, convective_adjust=true),
#             Sequential("KPP"; parameters=KPP.Parameters(), modify_predictor_fn=f, convective_adjust=true),
#             Sequential("T"; modify_predictor_fn=f, convective_adjust=true),
#             Sequential("dT"; modify_predictor_fn=f, convective_adjust=true),
# ]
#
# descriptor = "interpolation_conv_adj"
#
# for problem in problems
#
#     println("--*--*--*--*--*--$(problem)--*--*--*--*--*--")
#
#     𝒟_train     = LearnConvection.Data.data(train, problem; D=D, N=N);
#     𝒟_validate  = LearnConvection.Data.data(validate, problem; D=D, N=N);
#     𝒟_test      = LearnConvection.Data.data(test, problem; D=D, N=N);
#
#     train_validate_test(𝒟_train, 𝒟_validate, 𝒟_test, problem; log_γs=-2.0:0.1:2.0, distances=[euclidean_distance], descriptor="extrapolation_conv_adj")
#
# end

# 𝒢 = LearnConvection.GaussianProcess.model(𝒟_train; kernel=get_kernel(2,0.0,0.0,euclidean_distance))
# predict(𝒢, 𝒟_train; postprocessed=true)
# LearnConvection.GaussianProcess.get_me_true_check(𝒢, 𝒟_validate)
#
# k=5
# distance=euclidean_distance
# get_min_gamma(k, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.4:0.1:0.4)
# get_min_gamma_alpha(5, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.4:0.1:0.4)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
#using Pkg
#Pkg.add("Luxor")
#Pkg.add("Colors")

using Luxor, Colors
Drawing(1000,1000,"turtletwo.png")
origin()
background("midnightblue")

turtle = Turtle()
Pencolor(turtle, "red")
Penwidth(turtle, 4)
n = 5
for i in 1:4000
    global n
    Forward(turtle, n)
    Turn(turtle, 75)
    HueShift(turtle)
    n += 0.85
end
fontsize(20)
#Message(turtle, "finished")
finish()

Drawing(1000,1000,"turtlethree.png")
origin()
background("orange")

turtle = Turtle()
Pencolor(turtle, "blue")
Penwidth(turtle, 0.5)
n = 5
for i in 1:8000
    global n
    Forward(turtle, n)
    Turn(turtle, 95)
    HueShift(turtle)
    n += 0.85
end
fontsize(20)
#Message(turtle, "finished")
finish()

function draw(w,t)
 Drawing(1000,1000,"turtlethree.png")
 origin()
 background("midnightblue")
 turtle = Turtle()
 Pencolor(turtle, "blue")
 Penwidth(turtle, w)
 n=5
 for i in 1:8000
  Forward(turtle, n)
  Turn(turtle, t)
  HueShift(turtle)
  n += 0.85
 end
 fontsize(20)
 #Message(turtle, "finished")
 finish()
end

# Need to use the function now it's defined
draw(2,33.333)

Drawing(1000,1000,"turtletest.png")
origin()
background("midnightblue")

turtle = Turtle()
Pencolor(turtle, "red")
Penwidth(turtle, 4)
n = 5
for i in 1:50
    global n
    Forward(turtle, n)
    HueShift(turtle)
    n += 5
end
fontsize(20)
#Message(turtle, "finished")
finish()


text:
# STFT/ISTFT

function blackman(n::Integer)
    const a0, a1, a2 = 0.42, 0.5, 0.08
    t = 2*pi/(n-1)
    [a0 - a1*cos(t*k) + a2*cos(t*k*2) for k=0:n-1]
end

function hanning(n::Integer)
    [0.5*(1-cos(2*pi*k/(n-1))) for k=0:n-1]
end

# countframes returns the number of frames that will be processed.
function countframes{T<:Number}(x::Vector{T}, framelen::Int, hopsize::Int)
    div(length(x) - framelen, hopsize) + 1
end

# splitframes performs overlapping frame splitting.
function splitframes{T<:Number}(x::Vector{T},
                                framelen::Int=1024,
                                hopsize::Int=framelen/2)
    const N = countframes(x, framelen, hopsize)
    frames = Array(eltype(x), framelen, N)

    for i=1:N
        frames[:,i] = x[(i-1)*hopsize+1:(i-1)*hopsize+framelen]
    end

    return frames
end

# stft performs the Short-Time Fourier Transform (STFT) for real signals.
function stft{T<:Real}(x::Vector{T},
                       framelen::Int=1024,
                       hopsize::Int=div(framelen,2),
                       window=hanning(framelen))
    frames = splitframes(x, framelen, hopsize)

    const freqbins = div(framelen, 2) + 1
    spectrogram = Array(Complex64, freqbins, size(frames,2))
    for i=1:size(frames,2)
        spectrogram[:,i] = rfft(frames[:,i] .* window)
    end

    return spectrogram
end

# istft peforms the Inverse STFT to recover the original signal from STFT
# coefficients.
function istft{T<:Complex}(spectrogram::Matrix{T},
                           framelen::Int=1024,
                           hopsize::Int=div(framelen,2),
                           window=hanning(framelen))
    const numframes = size(spectrogram, 2)

    expectedlen = framelen + (numframes-1)*hopsize
    reconstructed = zeros(expectedlen)
    windowsum = zeros(expectedlen)
    const windowsquare = window .* window

    # Overlapping addition
    for i=1:numframes
        s, e = (i-1)*hopsize+1, (i-1)*hopsize+framelen
        r = irfft(spectrogram[:,i], framelen)
        reconstructed[s:e] += r .* window
        windowsum[s:e] += windowsquare
    end

    # Normalized by window
    for i=1:endof(reconstructed)
        # avoid zero division
        if windowsum[i] > 1.0e-7
            reconstructed[i] /= windowsum[i]
        end
    end

    return reconstructed
end


text:
# Format outputs of probSNFE and solveSNFE

function show(io::IO, p::ProbOutput1D)
    print(io,
        " ├─ Domain:       Ω × [0,T] = [$(p.Ω.x[1]),$(p.Ω.x[end])] × [0,$(p.Ω.t[end])]", '\n',
        " ├─ Spatial step: dx   = $(p.Ω.dx)", '\n',
        " ├─ Time step:    dt   = $(p.Ω.dt)", '\n',
        " ├─ Velocity:     v    = $(p.Ω.v)", '\n',
        " ├─ Delay rings:  umax = $(p.Ω.rings1D)")
end

function show(io::IO, p::ProbOutput2D)
    print(io,
        " ├─ Domain:       Ω × [0,T] = [$(p.Ω.x[1]),$(p.Ω.x[end])]^2 × [0,$(p.Ω.t[end])]", '\n',
        " ├─ Spatial step: dx   = $(p.Ω.dx)", '\n',
        " ├─ Time step:    dt   = $(p.Ω.dt)", '\n',
        " ├─ Velocity:     v    = $(p.Ω.v)", '\n',
        " ├─ Delay rings:  umax = $(p.Ω.rings2D)")
end

function show(io::IO, s::SolveOutDet1D)
    if issubset(s.tsaved,s.t)
        print(io,
              " ├─ Solution saved at all time instants", '\n',
              " ├─ 1D deterministic solution")
    else
        print(io,"Solution not saved at all time instants. Choose appropriate time instants.")
    end
end
function show(io::IO, s::SolveOutDet2D)
    if issubset(s.tsaved,s.t)
        print(io,
        " ├─ Solution saved at all time instants", '\n',
        " ├─ 2D deterministic solution")
    else
        print(io,"Solution not saved at all time instants. Choose appropriate time instants.")
    end
end
function show(io::IO, s::SolveOutSto1D)
    if issubset(s.tsaved,s.t)
        print(io,
        " ├─ Solution saved at all time instants", '\n',
        " ├─ 1D stochastic solution")
    else
        print(io,"Solution not saved at all time instants. Choose appropriate time instants.")
    end
end
function show(io::IO, s::SolveOutSto2D)
    if issubset(s.tsaved,s.t)
        print(io,
        " ├─ Solution saved at all time instants", '\n',
        " ├─ 2D stochastic solution")
    else
        print(io,"Solution not saved at all time instants. Choose appropriate time instants.")
    end
end


text:
# Autogenerated wrapper script for iODBC_debug_jll for x86_64-apple-darwin14
export iodbctestw, libiodbcdSYM, libiodbcinstdSYM, libiodbcinst, libiodbc, iodbctest

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `iodbctestw`
const iodbctestw_splitpath = ["bin", "iodbctestw"]

# This will be filled out by __init__() for all products, as it must be done at runtime
iodbctestw_path = ""

# iodbctestw-specific global declaration
function iodbctestw(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(iodbctestw_path)
    end
end

# Relative path to `libiodbcdSYM`
const libiodbcdSYM_splitpath = ["lib", "libiodbc.2.dylib.dSYM", "Contents", "Resources", "DWARF", "libiodbc.2.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libiodbcdSYM_path = ""

# libiodbcdSYM-specific global declaration
# This will be filled out by __init__()
libiodbcdSYM = ""

# Relative path to `libiodbcinstdSYM`
const libiodbcinstdSYM_splitpath = ["lib", "libiodbcinst.2.dylib.dSYM", "Contents", "Resources", "DWARF", "libiodbcinst.2.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libiodbcinstdSYM_path = ""

# libiodbcinstdSYM-specific global declaration
# This will be filled out by __init__()
libiodbcinstdSYM = ""

# Relative path to `libiodbcinst`
const libiodbcinst_splitpath = ["lib", "libiodbcinst.2.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libiodbcinst_path = ""

# libiodbcinst-specific global declaration
# This will be filled out by __init__()
libiodbcinst_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libiodbcinst = "@rpath/libiodbcinst.2.dylib"

# Relative path to `libiodbc`
const libiodbc_splitpath = ["lib", "libiodbc.2.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libiodbc_path = ""

# libiodbc-specific global declaration
# This will be filled out by __init__()
libiodbc_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libiodbc = "@rpath/libiodbc.2.dylib"

# Relative path to `iodbctest`
const iodbctest_splitpath = ["bin", "iodbctest"]

# This will be filled out by __init__() for all products, as it must be done at runtime
iodbctest_path = ""

# iodbctest-specific global declaration
function iodbctest(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(iodbctest_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"iODBC_debug")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global iodbctestw_path = normpath(joinpath(artifact_dir, iodbctestw_splitpath...))

    push!(PATH_list, dirname(iodbctestw_path))
    global libiodbcdSYM_path = normpath(joinpath(artifact_dir, libiodbcdSYM_splitpath...))

    global libiodbcdSYM = libiodbcdSYM_path
    global libiodbcinstdSYM_path = normpath(joinpath(artifact_dir, libiodbcinstdSYM_splitpath...))

    global libiodbcinstdSYM = libiodbcinstdSYM_path
    global libiodbcinst_path = normpath(joinpath(artifact_dir, libiodbcinst_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libiodbcinst_handle = dlopen(libiodbcinst_path)
    push!(LIBPATH_list, dirname(libiodbcinst_path))

    global libiodbc_path = normpath(joinpath(artifact_dir, libiodbc_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libiodbc_handle = dlopen(libiodbc_path)
    push!(LIBPATH_list, dirname(libiodbc_path))

    global iodbctest_path = normpath(joinpath(artifact_dir, iodbctest_splitpath...))

    push!(PATH_list, dirname(iodbctest_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
module Autoreload

using Compat

export arequire, areload,  aoptions_set, aoptions_get, smart_reload
export @aimport, @ausing

include("constants.jl")
include("smart_types.jl")
include("files.jl")
include("dependencies.jl")

DEBUG = false

info_debug(msg) = DEBUG ? info(msg) : nothing

function toggle_debug()
    global DEBUG
    DEBUG = !DEBUG
end

function aoptions_set(;kwargs...)
    global options
    for (key, value) in kwargs
        option_used = false
        for (i, (opt_key, opt_value)) in enumerate(options)
            if key==opt_key
                options[key] = value
                option_used = true
            end
        end
        if !option_used
            error("Option $key not recognized")
        end
    end
    options
end

aoptions_get() = copy(options)

function remove_file(filename)
    pop!(files, filename)
end

function arequire(filename=""; command= :on, depends_on=UTF8String[])
    if isempty(filename)
        return collect(keys(files))
    end
    original_filename = filename
    filename = find_file(standarize(filename), constants=options[:constants])
    filename == nothing && error("File $original_filename not found")
    if command in [:on, :on_depends]
        was_reloading = false
        if filename in keys(files)
            was_reloading = files[filename].should_reload
            remove_file(filename)
        end
        if command == :on
            should_reload = true
        else
            should_reload = was_reloading
        end
        files[filename] = AFile(should_reload, reload_mtime(filename), UTF8String[])
        parsed_file = parse_file(filename)
        auto_depends = extract_deps(parsed_file)
        auto_depends = [_.name for _ in auto_depends] #todo respect reload flag
        for i in 1:length(auto_depends)
            if !isabspath(auto_depends[i])
                auto_depends[i] = joinpath(dirname(filename), auto_depends[i])
            end
        end
        depends_on = vcat(auto_depends, depends_on)
        for d in depends_on
            d = standarize(d)
            if !haskey(files, d)
                d = find_file(d, base_file = filename)
                if d == nothing
                    error("Dependent file not found")
                end
                arequire(d, command = :on_depends)
            end
            push!(files[filename].deps, d)
        end
        if files[filename].should_reload
            include(find_file(filename))
        end
    elseif command == :off
        if haskey(files, filename)
            remove_file(filename)
        end
    else
        error("Command $command not recognized")
    end
    return
end

function smart_reload(file; kwargs...)
    global suppress_warnings
    info_debug("Initiating smart reload")
    original_file = file
    if !isabspath(file)
        file = find_file(standarize(file))
        if file == nothing
            error("File $(original_file) not found")
        end
    end
    suppress_warnings = true
    cd(dirname(file)) do
        info_debug("parsing $file")
        parsed = parse_file(file; kwargs...)
        info_debug("extracting modules")
        module_paths = extract_modules(parsed)
        info_debug("looping")
        for (module_name, e) in module_paths
            info_debug("reloading $module_name")
            reload_module(module_name, e)
            info_debug("done reloading $module_name")
        end
    end
    suppress_warnings = false
    return
end

function try_reload(file; kwargs...)
    smart_reload(file; kwargs...)
    # try
    #     smart_reload(file; kwargs...)
    # catch err
    #     if options[:verbose_level] == :warn
    #         warn("Could not autoreload $(file):\n$(err)")
    #     end
    # end
end

function is_equal_dict(d1::Dict, d2::Dict)
    for key in keys(d1)
        if d1[key]!=d2[key]
            return false
        end
    end
    return true
end

function areload(command= :force; kwargs...)
    if command == :use_state
        if options[:state] == :off
            return
        end
    elseif command == :off
        options[:state] = :off
        return
    elseif command == :on
        options[:state] = :on
    end
    dependencies = get_dependency_graph()
    file_order = topological_sort(dependencies)
    should_reload = [filename=>false for filename in file_order]
    marked_for_mtime_update = UTF8String[]
    for (i, file) in enumerate(file_order)
        file_time = files[file].mtime
        if reload_mtime(file) > file_time
            should_reload[file] = true
            # files[file].mtime = reload_mtime(file)
            push!(marked_for_mtime_update, file)
        end
    end
    old_reload = copy(should_reload)
    while true
        for (child, parents) in dependencies
            for parent in parents
                if should_reload[parent]
                    should_reload[child] = true
                end
            end
        end
        if is_equal_dict(old_reload, should_reload)
            break
        end
        old_reload = copy(should_reload)
    end

    for file in file_order
        if should_reload[file] && files[file].should_reload
            try_reload(file; constants=options[:constants], kwargs...)
        end
    end

    for file in marked_for_mtime_update
        files[file].mtime = reload_mtime(file)
    end

    return
end

function areload_hook()
    areload(:use_state)
end

if isdefined(Main, :IJulia)
    IJulia = Main.IJulia
    try
        IJulia.push_preexecute_hook(areload_hook)
    catch err
        warn("Could not add IJulia hooks:\n$(err)")
    end
end

include("macros.jl")

end


text:
export updatecontent!
################################################################
# Updaters
################################################################
# Update based on the elm
################################################################
# Nothing remove
"""
    updatecontent!(value, elm::Node|Nothing, node, argAmlType)
    updatecontent!(value, elms::Vector{Node|Nothing}, node, argAmlType)

In the given `node`, update the content of all the given `elm` or `elms` to `value`.

    updatecontent!(value, name::String, node, argAmlType)

In the given `node`, update the content of all the elements with given `name` to `value`.
"""
function updatecontent!(value::Nothing, elm::Node, node::Node, argAmlType::Type)
    unlink!(elm)
end

# Nothing Nothing
function updatecontent!(value::Nothing, elm::Nothing, node::Node, argAmlType::Type)
    # nothing
end

################################################################
# Number and String
function updatecontent!(value::T, elm::Node, node::Node, argAmlType::Type) where {T<:Union{AbstractString, Number}}
    elm.content = value
end

# Number and String Vector
function updatecontent!(values::Vector{T}, elms::Vector{Node}, node::Node, argAmlType::Type) where {T<:Union{AbstractString, Number}}
    for (i, elm) in enumerate(elms)
        elm.content = values[i]
    end
end

################################################################
# General Type
function updatecontent!(value::T, elm::Node, node::Node, argAmlType::Type)  where {T}
    if hasmethod(string, Tuple{T})
        elm.content = string(value)
    else
        unlink!(elm)
        link!(node, value.aml)
    end
end

# General Type Vector
function updatecontent!(values::Vector{T}, elms::Vector{Node}, node::Node, argAmlType::Type)  where {T}
    for (i, elm) in enumerate(elms)
        if isnothing(values[i])
            unlink!(elm)
        else
            if hasmethod(string, Tuple{T})
                elm.content = string(values[i])
            else
                unlink!(elm)
                link!(node, values[i].aml)
            end
        end
    end
end

################################################################
# Update based on the name
################################################################
# Documents
function updatecontent!(value, name::String, doc::Document, argAmlType::Type{<:AbsNode})
    updatecontent!(value, name, root(doc), argAmlType)
end
################################################################
# Nodes
# Single Updater

function updatecontent!(value, name::String, node::Node, argAmlType::Type{<:Union{AbsNormal, AbsText}})
    elm = findfirst(name, node, argAmlType)
    if !isnothing(elm)
        updatecontent!(value, elm, node, argAmlType)
    else
        # addnode! if nothing is found
        addnode!(node, name, value, argAmlType)
    end
end

# Fast attribute updater
function updatecontent!(value, name::String, node::Node, argAmlType::Type{<:AbsAttribute})
    if haskey(node, name)
        node[name] = value
    else
        # addnode! if nothing is found
        addnode!(node, name, value, argAmlType)
    end
end

################################################################
# Vector update

function updatecontent!(values::Vector, name::String, node::Node, argAmlType::Type{<:Union{AbsNormal, AbsText}})
    elms = findall(name, node, argAmlType)
    if !isnothing(elms)
        updatecontent!(values, elms, node, argAmlType)
    else
        # addnode! if nothing is found
        addnode!(node, name, values, argAmlType)
    end
end

# Fast attribute updater
function updatecontent!(values::Vector, name::String, node::Node, argAmlType::Type{<:AbsAttribute})
    if haskey(node, name)
        node[name] .= values
    else
        # addnode! if nothing is found
        addnode!(node, name, values, argAmlType)
    end
end

################################################################
# Dict Updating
function updatecontent!(value::Type{AbstractDict}, name, node, argAmlType)
    throw(MethodError("Dicts are not supported for extraction/updating"))
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis
using ConstrainedControl
using LinearAlgebra

path = "examples/examples_files/sawyer_arm.urdf"
mech = Mechanism(path, floating=false, g = 0.0)

setPosition!(mech,mech.eqconstraints["right_j0"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j1"],[-pi/2])
setPosition!(mech,mech.eqconstraints["right_j2"],[-pi/2])
setPosition!(mech,mech.eqconstraints["right_j3"],[pi/2])

xd=[mech.bodies[i].state.xc for i=1:length(mech.bodies)]
qd=[mech.bodies[i].state.qc for i=1:length(mech.bodies)]

setPosition!(mech,mech.eqconstraints["right_j0"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j1"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j2"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j3"],[0.0])

Q = [diagm(ones(12))*1000.0 for i=1:7]
R = [ones(1,1) for i=1:7]

ConstrainedDynamics.deactivateConstraints!(mech)
lqr = LQR(mech, getid.(mech.bodies), getid.(mech.eqconstraints), Q, R, 20., xd=xd, qd=qd)
ConstrainedDynamics.activateConstraints!(mech)

storage = simulate!(mech,20,lqr,record = true)
visualize(mech,storage)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _eq2hor(ra::T, dec::T, jd::T, lat::T, lon::T, altitude::T,
                 pressure::T, temperature::T, ws::Bool, B1950::Bool,
                 precession::Bool, nutate::Bool, aberration::Bool,
                 refract::Bool) where {T<:AbstractFloat}

    j_now = (jd - J2000) / JULIANYEAR + 2000

    if precession
        if B1950
            ra, dec = precess(ra, dec, 1950, j_now, FK4=true)
        else
            ra, dec = precess(ra, dec, 2000, j_now)
        end
    end
    dra1, ddec1, eps, d_psi = co_nutate(jd, ra, dec)

    if aberration
        dra2, ddec2 = co_aberration(jd, ra, dec, eps)
        ra += dra2 / 3600
        dec += ddec2 / 3600
    end

    if nutate
       ra += dra1
       dec += ddec1
    end
    last = 15 * ct2lst(lon, jd) + d_psi * cos(eps) / 3600
    ha = mod(last - ra, 360)
    alt, az = hadec2altaz(ha, dec, lat, ws=ws)

    if refract
        alt = co_refract(alt, altitude, pressure, temperature, to_observe = true)
    end
    return alt, az, ha
end

eq2hor(ra::Real, dec::Real, jd::Real, lat::Real, lon::Real, altitude::Real=0;
       ws::Bool=false, B1950::Bool=false, precession::Bool=true, nutate::Bool=true,
       aberration::Bool=true, refract::Bool=true, pressure::Real=NaN,
       temperature::Real=NaN) =
           _eq2hor(promote(float(ra), float(dec), float(jd), float(lat), float(lon),
                           float(altitude), float(temperature), float(pressure))..., ws,
                   B1950, precession, nutate, aberration, refract)

eq2hor(ra::Real, dec::Real, jd::Real, obsname::AbstractString; kwargs...) =
    eq2hor(ra, dec, jd, observatories[obsname].latitude, observatories[obsname].longitude,
           observatories[obsname].altitude; kwargs...)

"""
    eq2hor(ra, dec, jd[, obsname; ws=false, B1950=false, precession=true, nutate=true,
           aberration=true, refract=true, pressure=NaN, temperature=NaN]) -> alt, az, ha

    eq2hor(ra, dec, jd, lat, lon[, altitude=0; ws=false, B1950=false,
           precession=true, nutate=true, aberration=true, refract=true,
           pressure=NaN, temperature=NaN]) -> alt, az, ha

### Purpose

Convert celestial  (ra-dec) coords to local horizon coords (alt-az).

### Explanation

This code calculates horizon (alt,az) coordinates from equatorial (ra,dec) coords.
It performs precession, nutation, aberration, and refraction corrections.

### Arguments

This function has two base methods.  With one you can specify the name of the observatory,
if present in `AstroLib.observatories`, with the other one you can provide the coordinates
of the observing site and, optionally, the altitude.

Common mandatory arguments:

* `ra`: right ascension of object, in degrees
* `dec`: declination of object, in degrees
* `jd`: julian date

Other positional arguments:

* `obsname`: set this to a valid observatory name in `AstroLib.observatories`.

or

* `lat`: north geodetic latitude of location, in degrees.
* `lon`: AST longitude of location, in degrees. You can specify west longitude with a
  negative sign.
* `altitude`: the altitude of the observing location, in meters.  It is `0` by default

Optional keyword arguments:

* `ws` (optional boolean keyword): set this to `true` to get the azimuth measured
* `B1950` (optional boolean keyword): set this to `true` if the ra and dec
  are specified in B1950 (FK4 coordinates) instead of J2000 (FK5). This is `false` by
  default
* `precession` (optional boolean keyword): set this to `false` for no precession
  correction, `true` by default
* `nutate` (optional boolean keyword): set this to `false` for no nutation,
  `true` by default
* `aberration` (optional boolean keyword): set this to `false` for no aberration
  correction, `true` by default
* `refract` (optional boolean keyword): set this to `false` for no refraction
  correction, `true` by default
* `pressure` (optional keyword): the pressure at the observing location, in millibars.
  Default value is `NaN`
* `temperature` (optional keyword): the temperature at the observing location, in Kelvins.
  Default value is `NaN`

### Output

* `alt`: altitude of horizon coords, in degrees
* `az`: azimuth angle measured East from North (unless ws is `true`), in degrees
* `ha`: hour angle, in degrees

### Example

```jldoctest
julia> using AstroLib

julia> alt_o, az_o = eq2hor(ten(6,40,58.2)*15, ten(9,53,44), 2460107.25, ten(50,31,36),
                            ten(6,51,18), 369, pressure = 980, temperature=283)
(16.423991509721567, 265.60656932130564, 76.11502253130612)

julia> adstring(az_o, alt_o)
" 17 42 25.6  +16 25 26"
```

### Notes

Code of this function is based on IDL Astronomy User's Library.
"""
eq2hor

# TODO: Make eq2hor type-stable, which it isn't currently because of keyword arguments
# Note that the inner function `_eq2hor` is type stable


text:
struct Dataset{N<:Name, F<:Format, T<:Task}
    shuffle::Bool
    seed::Int64
    format_args
    task_args

    function Dataset(
        N::Type{<:Name};
        shuffle::Bool = false,
        seed::Int64 = 1234,
        kwargs...
    )
        F = format(N)
        T = task(N)
        return new{N, F, T}(shuffle, seed, args(F; kwargs...), args(T; kwargs...))
    end
end

function Base.show(io::IO, data::Dataset{N, F, T}) where {N, F, T}
    args = (shuffle = data.shuffle, seed = data.seed)
    print(io, nameof(N), merge(args, data.format_args, data.task_args))
    return
end

function load(
    dataset::Dataset{N, F, T},
    type::Symbol;
    dopostprocess::Bool = true,
) where {N<:Name, F<:Format, T<:Task}

    seed = dataset.seed
    shuffle = dataset.shuffle
    obsdim = samplesdim(F)
    data = data_shuffle(loadraw(N, type), obsdim, shuffle; seed)

    if dopostprocess
        return postprocess(dataset, data)
    else
        return data
    end
end

function postprocess(d::Dataset{N, F, T}, data) where {N<:Name, F<:Format, T<:Task}
    data = postprocess(N, T, data; d.task_args...)
    data = postprocess(N, F, data; d.format_args...)
    return data
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider