text:
# beefblup
# Julia package for performing single-variate BLUP to find beef cattle
# breeding values
# (C) 2021 Thomas A. Christensen II
# Licensed under BSD-3-Clause License
# cSpell:includeRegExp #.*
# cSpell:includeRegExp ("""|''')[^\1]*\1

module BeefBLUP

# Import the required packages
using CSV
using DataFrames
using LinearAlgebra
using Dates
using Gtk

# Main entry-level function - acts just like the script
function beefblup()

    # Ask for an input spreadsheet
    path = open_dialog_native(
    "Select a beefblup worksheet",
    GtkNullContainer(),
    ("*.csv", GtkFileFilter("*.csv", name="beefblup worksheet"))
    )

    # Ask for an output text filename
    savepath = save_dialog_native(
    "Save your beefblup results",
    GtkNullContainer(),
    (GtkFileFilter("*.txt", name="Results file"),
    "*.txt")
    )

    # Ask for heritability
    print("What is the heritability for this trait?> ")
    h2 = parse(Float64, readline(stdin))

    beefblup(path, savepath, h2)

end

function beefblup(datafile::String, h2::Float64)
    # Assume the data is named the same as the file without the trailing extension
    dataname = join(split(datafile, ".")[1:end - 1])

    # Create a new results name
    resultsfile = string(dataname, "_results.txt")

    # Pass this info on to the worker
    beefblup(datafile, resultsfile, h2)
end

# Main worker function, can perform all the work if given all the user input
function beefblup(path::String, savepath::String, h2::Float64)

    # Import data from a suitable spreadsheet
    data = DataFrame(CSV.File(path))

    # Make sure the data is in the proper format
    renamecolstospec!(data)

    # Sort the array by date
    sort!(data, :birthdate)

    # Define fields to hold id values for animals and their parents
    numanimals = length(data.id)

    # Calculate the relationship matrix
    A = additiverelationshipmatrix(data.id, data.dam, data.sire)

    # Extract all of the fixed effects
    fixedeffectdata = data[:,5:end-1]

    (X, fixedeffects) = fixedeffectmatrix(fixedeffectdata)

    # Extract the observed data
    Y = convert(Array{Float64}, data[:,end])

    # The random effects matrix
    Z = Matrix{Int}(I, numanimals, numanimals)

    # Remove items where there is no data
    nullobs = findall(isnothing, Y)
    Z[nullobs, nullobs] .= 0

    # Calculate heritability
    λ = (1 - h2) / h2

    # Use the mixed-model equations
    MME = [X' * X X' * Z; Z' * X (Z' * Z) + (inv(A) .* λ)]
    MMY = [X' * Y; Z' * Y]
    solutions = MME \ MMY

    # Find the accuracies
    diaginv = diag(inv(MME))
    reliability = ones(Float64, length(diaginv)) .- diaginv .* λ

    # Find how many traits we found BLUE for
    numgroups = length(reliability) - numanimals

    # Split the BLUP and BLUE results
    β = solutions[1:numgroups]
    μ = solutions[numgroups+1:end]
    μ_reliability = reliability[numgroups+1:end]

    # Extract the names of the traits
    traitname = names(data)[end]

    # Start printing results to output
    fileID = open(savepath, "w")
    write(fileID, "beefblup Results Report\n")
    write(fileID, "Produced using beefblup (")
    write(fileID, "https://github.com/millironx/beefblup")
    write(fileID, ")\n\n")
    write(fileID, "Input:\t")
    write(fileID, path)
    write(fileID, "\nAnalysis performed:\t")
    write(fileID, string(Dates.today()))
    write(fileID, "\nTrait examined:\t")
    write(fileID, traitname)
    write(fileID, "\n\n")

    # Print base population stats
    write(fileID, "Base Population:\n")
    for i in 1:length(fixedeffects)
        effect = fixedeffects[i]
        write(fileID, "\t")
        write(fileID, string(effect.name))
        write(fileID, ":\t")
        write(fileID, string(effect.basetrait))
        write(fileID, "\n")
    end

    write(fileID, "\tMean ")
    write(fileID, traitname)
    write(fileID, ":\t")
    write(fileID, string(solutions[1]))
    write(fileID, "\n\n")

    # Contemporary group adjustments
    counter = 2
    write(fileID, "Contemporary Group Effects:\n")
    for i in 1:length(fixedeffects)
        effect = fixedeffects[i]
        write(fileID, "\t")
        write(fileID, effect.name)
        write(fileID, "\tEffect\n")
        for j in 1:length(effect.alltraits)
            types = effect.alltraits[j]
            write(fileID, "\t")
            write(fileID, string(types))
            write(fileID, "\t")
            write(fileID, string(β[counter]))
            write(fileID, "\n")

            counter = counter + 1
        end
        write(fileID, "\n")
    end
    write(fileID, "\n")

    # Expected breeding values
    write(fileID, "Expected Breeding Values:\n")
    write(fileID, "\tID\tEBV\tReliability\n")
    for i in 1:numanimals
        write(fileID, "\t")
        write(fileID, string(data.id[i]))
        write(fileID, "\t")
        write(fileID, string(μ[i]))
        write(fileID, "\t")
        write(fileID, string(μ_reliability[i]))
        write(fileID, "\n")
    end

    write(fileID, "\n - END REPORT -")
    close(fileID)

end

"""
    fixedeffectmatrix(fixedeffectdata::DataFrame)

Creates contemporary groupings and the fixed-effect incidence matrix based on the fixed
effects listed in `fixedeffectdata`.

Returns a tuple `(X::Matrix{Int}, fixedeffects::Array{FixedEffect})` in which `X` is the
actual matrix, and `fixedeffects` is the contemporary groupings.
"""
function fixedeffectmatrix(fixedeffectdata::DataFrame)
    # Declare an empty return matrix
    fixedeffects = FixedEffect[]

    # Add each trait to the array
    for i in 1:size(fixedeffectdata)[2]
        name = names(fixedeffectdata)[i]
        traits = eachcol(fixedeffectdata)[i]

        if length(unique(traits)) > 1
            push!(fixedeffects, FixedEffect(name, traits))
        else
            @warn string("column '", name, "' does not have any unique animals and will be dropped from analysis")
            pname = propertynames(fixedeffectdata)[i]
            DataFrames.select!(fixedeffectdata, Not(pname))
        end
    end

    X = ones(Int64, (size(fixedeffectdata)[1], 1))

    for i in 1:length(fixedeffects)
        trait = fixedeffects[i]
        for phenotype in trait.alltraits
            X = cat(X, Int64.(fixedeffectdata[:,i] .== phenotype), dims=2)
        end
    end

    return X, fixedeffects
end

"""
    additiverelationshipmatrix(id, dam, sire)

Returns the additive numerator relationship matrix based on the pedigree provided in `dam`
and `sire` for animals in `id`.

"""
function additiverelationshipmatrix(id::AbstractVector, damid::AbstractVector, sireid::AbstractVector)
    # Sanity-check for valid pedigree
    if !(length(id) == length(damid) && length(damid) == length(sireid))
        throw(ArgumentError("id, dam, and sire must be of the same length"))
    end

    # Convert to positions
    dam = indexin(damid, id)
    sire = indexin(sireid, id)

    # Calculate loop iterations
    numanimals = length(dam)

    # Create an empty matrix for the additive relationship matrix
    A = zeros(numanimals, numanimals)

    # Create the additive relationship matrix by the FORTRAN method presented by
    # Henderson
    for i in 1:numanimals
        if !isnothing(dam[i]) && !isnothing(sire[i])
            for j in 1:(i - 1)
                A[j,i] = 0.5 * (A[j,sire[i]] + A[j,dam[i]])
                A[i,j] = A[j,i]
            end
            A[i,i] = 1 + 0.5 * A[sire[i], dam[i]]
        elseif !isnothing(dam[i]) && isnothing(sire[i])
        for j in 1:(i - 1)
            A[j,i] = 0.5 * A[j,dam[i]]
            A[i,j] = A[j,i]
        end
        A[i,i] = 1
    elseif isnothing(dam[i]) && !isnothing(sire[i])
        for j in 1:(i - 1)
            A[j,i] = 0.5 * A[j,sire[i]]
            A[i,j] = A[j,i]
        end
        A[i,i] = 1
    else
        for j in 1:(i - 1)
            A[j,i] = 0
            A[i,j] = 0
        end
        A[i,i] = 1
        end
    end

    return A
end

"""
    renamecolstospec(::DataFrame)

Renames the first four columns of the beefblup data sheet so that they can be referred to by
name instead of by column index, regardless of user input.
"""
function renamecolstospec!(df::DataFrame)
    # Pull out the fixed-effect and observation name
    othernames = propertynames(df)[5:end]

    # Put specification column names and user-defined names together
    allnames = cat([:id, :birthdate, :dam, :sire], othernames, dims=1)

    # Rename in the DataFrame
    rename!(df, allnames, makeunique=true)
    return df
end

struct FixedEffect
    name::String
    basetrait::Any
    alltraits::AbstractArray{Any}
end

function FixedEffect(name::String, incidences)
    basetrait = last(unique(incidences))
    types = unique(incidences)[1:end-1]
    return FixedEffect(name, basetrait, types)
end

end

max_stars_repo_name: MillironX/beefblup

max_issues_repo_name: MillironX/beefblup

max_forks_repo_name: MillironX/beefblup


text:
# ---
# title: 473. Matchsticks to Square
# id: problem473
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search
# link: <https://leetcode.com/problems/matchsticks-to-square/description/>
# hidden: true
# ---
#
# Remember the story of Little Match Girl? By now, you know exactly what
# matchsticks the little match girl has, please find out a way you can make one
# square by using up all those matchsticks. You should not break any stick, but
# you can link them up, and each matchstick must be used **exactly** one time.
#
# Your input will be several matchsticks the girl has, represented with their
# stick length. Your output will either be true or false, to represent whether
# you could make one square using all the matchsticks the little match girl has.
#
# **Example 1:**
#
#
#
#     Input: [1,1,2,2,2]
#     Output: true
#
#     Explanation: You can form a square with length 2, one side of the square came two sticks with length 1.
#
#
# **Example 2:**
#
#
#
#     Input: [3,3,3,3,4]
#     Output: false
#
#     Explanation: You cannot find a way to form a square with all the matchsticks.
#
#
# **Note:**
#
#   1. The length sum of the given matchsticks is in the range of `0` to `10^9`.
#   2. The length of the given matchstick array will not exceed `15`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
paths_file = "$(Pkg.dir("AST"))/src/paths.jl"

if !isfile(paths_file)
    atsp,grasp = try
        ENV["ATSP"], ENV["GRASP"]
    catch
        error("Did not find ATSP and/or GRASP environment variables. Please define them and rebuild AST.")
    end

    isfile("$atsp/hf") || error("Could not find hf exe in $atsp")
    isfile("$grasp/rmcdhf") || error("Could not find rmcdhf exe in $grasp")

    println("Found ATSP2K at $atsp and GRASP2K at $grasp")

    open(paths_file, "w") do file
        write(file, """atsp = \"$atsp\"
    grasp = \"$grasp\"
    """)
    end
end


text:
# Copyright (c) 2021
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using LinearAlgebra;

"""
    SL0(A::Matrix{Float64}, b::Vector{Float64}; maxiter=100, epsilon=0.01, sigma_decrease_factor=0.85)

Find the solution to Ax=b using Iterative Recursive Least Squares.

### Input

- `A`       -- Matrix: Ax=b
- `b`       -- Vector: Ax=b
- `sigma_decrease_factor` -- (optional) number of optmization iterations
- `epsilon` -- (optional) threshold to stop optimizing
- `maxiter` -- (optional) max number of iterations

### Output

Solution to Ax=b (Vector{Float64})

### Algorithm

Smoothed L0 (http://ee.sharif.edu/~SLzero/)

"""
function SL0(A::Matrix{Float64},
             b::Vector{Float64};
             sigma_decrease_factor=.85,
             maxiter=150,
             epsilon=.001)

    # set up the locals
    local mu_0, L, A_pinv, s, x;

    # assign constants
    mu_0 = 2;          # The  value  of  mu_0  scales  the sequence of mu
    L = 3;             # number  of  iterations of the internal (steepest ascent) loop
    A_pinv = pinv(A);  # pseudo-inverse of matrix A defined by A_pinv=A'*inv(A*A')

    # initialize the solution
    s = A_pinv*b;
    sigma = 2*maximum(abs.(s));

    for j = 1:maxiter

        for i = 1:L
            delta = s.*exp.(-abs.(s).^2/sigma^2);
            s -= mu_0*delta;
            s -= A_pinv*(A*s - b);
        end
        sigma *= sigma_decrease_factor;
    end
    x = s;
    i = abs.(x) .< epsilon;
    x[i] = zeros(sum(i));
    return x
end

# A = randn(10, 50); b = randn(10);
# x = SL0(A, b);
# println(x)


text:
using Pkg
pkg"activate .."
push!(LOAD_PATH, "../src/")
using Documenter, Hamburg

makedocs(sitename="Hamburg Documentation")


text:
@timedtestset "addition" begin
    @timedtestset "Correlated GPs" begin
        rng, N, N′, D, gpc = MersenneTwister(123456), 5, 6, 2, GPC()
        X, X′ = ColVecs(randn(rng, D, N)), ColVecs(randn(rng, D, N′))
        f1 = wrap(GP(1, SEKernel()), gpc)
        f2 = wrap(GP(2, SEKernel()), gpc)
        f3 = f1 + f2
        f4 = f1 + f3
        f5 = f3 + f4

        for (n, (fp, fa, fb)) in enumerate([(f3, f1, f2), (f4, f1, f3), (f5, f3, f4)])
            Σp = cov(fa(X)) + cov(fb(X)) + cov(fa(X), fb(X)) + cov(fb(X), fa(X))
            ΣpXX′ = cov(fa(X), fa(X′)) + cov(fb(X), fb(X′)) + cov(fa(X), fb(X′)) +
                cov(fa(X), fb(X′))
            @test mean(fp(X)) ≈ mean(fa(X)) + mean(fb(X))
            @test cov(fp(X)) ≈ Σp
            @test cov(fp(X), fp(X′)) ≈ ΣpXX′
            @test cov(fp(X′), fp(X)) ≈ transpose(ΣpXX′)
            @test cov(fp(X), fa(X′)) ≈ cov(fa(X), fa(X′)) + cov(fb(X), fa(X′))
            @test cov(fp(X′), fa(X)) ≈ cov(fa(X′), fa(X)) + cov(fb(X′), fa(X))
            @test cov(fa(X), fp(X′)) ≈ cov(fb(X), fa(X′)) + cov(fa(X), fa(X′))
            @test cov(fa(X′), fp(X)) ≈ cov(fb(X′), fa(X)) + cov(fa(X′), fa(X))
        end

        @timedtestset "Consistency Tests" begin
            P, Q = 4, 3
            x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
            abstractgp_interface_tests(f3, f1, x0, x1, x2, x3)
            abstractgp_interface_tests(f2 - f1, f1, x0, x1, x2, x3)
        end
    end
    @timedtestset "Verify mean / kernel numerically" begin
        rng, N, D = MersenneTwister(123456), 5, 6
        X = ColVecs(randn(rng, D, N))
        c, f = randn(rng), wrap(GP(5, SEKernel()), GPC())

        @test mean((f + c)(X)) == mean(f(X)) .+ c
        @test mean((f + c)(X)) == c .+ mean(f(X))
        @test cov((f + c)(X)) == cov(f(X))
        @test cov((c + f)(X)) == cov(f(X))

        @test mean((f - c)(X)) == mean(f(X)) .- c
        @test mean((c - f)(X)) == c .- mean(f(X))
        @test cov((f - c)(X)) == cov(f(X))
        @test cov((c - f)(X)) == cov(f(X))

        x = randn(rng, N + D)
        @test mean((f + sin)(x)) == mean(f(x)) + map(sin, x)
        @test mean((sin + f)(x)) == map(sin, x) + mean(f(x))
        @test cov((f + sin)(x)) == cov(f(x))
        @test cov((sin + f)(x)) == cov(f(x))

        @timedtestset "Consistency Tests" begin
            P, Q = 5, 3
            x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
            abstractgp_interface_tests(c + f, f, x0, x1, x2, x3)
        end
    end
    @timedtestset "Standardised Tests (independent sum)" begin
        standard_1D_tests(
            MersenneTwister(123456),
            Dict(:l1=>0.5, :l2=>2.3),
            θ->begin
                gpc = GPC()
                f1 = θ[:l1] * wrap(GP(sin, SEKernel()), gpc)
                f2 = θ[:l2] * wrap(GP(cos, SEKernel()), gpc)
                f3 = f1 + f2
                return f3, f3
            end,
            collect(range(-1.5, 1.5; length=5)),
            collect(range(-1.0, 0.5; length=3)),
        )
    end
    @timedtestset "Standardised Tests (correlated sum)" begin
        standard_1D_tests(
            MersenneTwister(123456),
            Dict(:l1=>0.5, :l2=>2.3),
            θ->begin
                gpc = GPC()
                f1 = θ[:l1] * wrap(GP(sin, SEKernel()), gpc)
                f2 = θ[:l2] * wrap(GP(cos, SEKernel()), gpc)
                f3 = f1 + f2
                f4 = f1 + f3
                f5 = f3 + f4
                return f5, f5
            end,
            collect(range(-1.5, 1.5; length=5)),
            collect(range(-1.0, 0.5; length=3)),
        )
    end
end


text:
module NEAT

import Base.split

include("config.jl")

# track count of various global variables & holds refernce to the config
type Global
    speciesCnt::Int64
    chromosomeCnt::Int64
    nodeCnt::Int64
    innov_number::Int64
    innovations::Dict{(Int64,Int64),Int64}
    cf::Config
    function Global(cf::Config)
        new(0,0,0,0,Dict{(Int64,Int64),Int64}(),cf) # global dictionary
    end
end

include("genome.jl")
include("chromosome.jl")
include("species.jl")
include("population.jl")
include("../networks/nn.jl")

# export Input, Output, Hidden

end # module

max_stars_repo_name: Andy-P/NEAT

max_issues_repo_name: Andy-P/NEAT

max_forks_repo_name: Andy-P/NEAT


text:
module example

include("propwing.jl")
import CCBlade
using PyPlot
close("all")

# ------- common parameters (standard metric units) ------
rho = 1.225  # air density
Vinf = 10.0  # freestream speed
J = 0.5  # prop advance ratio
xtowing = 18.0  # distance from prop to wing (for wake contraction) - assuming it is the same for all props
etaprop_center = [40.0/94, 70.0/94]  # spanwise location of prop (center)
cw = 1.0  # 1 if clockwise, -1 if counter-clockwise (when viewed from upstream)
gamma = 3*pi/180  # flight path angle
alpha_wing_relative_to_fuselage = 2*pi/180 # angle of attack of wing relative to fuselage reference line

alpha = gamma + alpha_wing_relative_to_fuselage  # angle of attack

# --------- propeller --------
# geometry
Rhub = 3*.5
Rtip = 3*3.0

r = 3*[0.526, 0.628, 0.729, 0.831, 0.9132, 0.9586, 1.0332,
     1.1128, 1.1925, 1.2722, 1.3519, 1.4316, 1.5114, 1.5911,
     1.6708, 1.7505, 1.8302, 1.9099, 1.9896, 2.0693, 2.1490, 2.2287,
     2.3084, 2.3881, 2.4678, 2.5475, 2.6273, 2.7070, 2.7867, 2.8661, 2.9410]
chord_prop = 3*[0.6270, 0.6255, 0.6231, 0.6199, 0.6165, 0.6125, 0.6054, 0.5973, 0.5887,
          0.5794, 0.5695, 0.5590, 0.5479, 0.5362, 0.5240, 0.5111, 0.4977,
          0.4836, 0.4689, 0.4537, 0.4379, 0.4214, 0.4044, 0.3867, 0.3685,
          0.3497, 0.3303, 0.3103, 0.2897, 0.2618, 0.1920]

twist_prop = pi/180.0*[40.2273, 38.7657, 37.3913, 36.0981, 34.8803, 33.5899, 31.6400,
                   29.7730, 28.0952, 26.5833, 25.2155, 23.9736, 22.8421, 21.8075,
                   20.8586, 19.9855, 19.1800, 18.4347, 17.7434, 17.1005, 16.5013,
                   15.9417, 15.4179, 14.9266, 14.4650, 14.0306, 13.6210, 13.2343,
                   12.8685, 12.5233, 12.2138]
B = 3  # number of blades

aftype = CCBlade.af_from_aerodynfile("NACA64_A17.dat")
n = length(r)
af = Array{CCBlade.AirfoilData}(n)
for i = 1:n
    af[i] = aftype
end

Dprop = Rtip*2  # propeller diameter
Omega = Vinf/(J*Dprop)*2*pi  # rotation rate

uwake, vwake, T, Q, reff = propwing.propanalysis(Rhub, Rtip, r, chord_prop, twist_prop, B, af, Vinf, Omega, rho, xtowing)

# --------- wing geometry -------
sectionspan = [29, 65]
chord_wing = [43, 26, 11]
twist_wing = [0, 0, 0]*pi/180
tc = [0.13, 0.12, 0.11]
sweep = [25, 30]*pi/180
dihedral = [0, 0]
N = [29, 65]

# L, Di, CP, cl, cllocal, Vinfeff = propwing.winganalysis(sectionspan, chord_wing, twist_wing, tc, sweep, dihedral, N, alpha, rho, Vinf, reff, uwake, vwake, etaprop_center, cw)
# println(L)
# println(Di)
#
# figure()
# plot(CP.y, cl)
# plot(CP.y, cllocal)
#
# # ------ viscous drag -----
# # TODO
# Dp = 0.0
# D = Di + Dp
#
# # ------- add prop thrust to lift --------
# L += T*sin(gamma)  # add prop thrust to lift
# thrust_ratio = T*cos(gamma)/D  # must be larger than 1

# -------- diameter checking --------
nd = 100
dratio = linspace(0.02, 0.8, nd)
CLvec = zeros(nd)
CDivec = zeros(nd)

etaprop_center = [40.5828/94]  # spanwise location of prop (center)

for i = 1:nd
    rprop = reff/reff[end]*dratio[i]/2*sum(sectionspan)
    L, Di, CP, cl, cllocal, Vinfeff, alphaeff = propwing.winganalysis(sectionspan, chord_wing, twist_wing, tc, sweep, dihedral, N, alpha, rho, Vinf, rprop, uwake, vwake, etaprop_center, cw)
    q = 0.5*rho*Vinf^2
    S = 2*sum(CP.chord.*CP.ds)
    CLvec[i] = L/(q*S)
    CDivec[i] = Di/(q*S)
end

# println(CLvec)
# println(CDivec)

CLbad = [0.382905, 0.386424, 0.386774, 0.386081, 0.389872, 0.390167, 0.392373, 0.39475, 0.394911, 0.397192, 0.397608, 0.397071, 0.400073, 0.400052, 0.399881, 0.402262, 0.40314, 0.405263, 0.406519, 0.406577, 0.409258, 0.410152, 0.409989, 0.41318, 0.414047, 0.415491, 0.417441, 0.417875, 0.420247, 0.421134, 0.42098, 0.42363, 0.424618, 0.425658, 0.428337, 0.429073, 0.430862, 0.432741, 0.433258, 0.436002, 0.437109, 0.437732, 0.44078, 0.441917, 0.443226, 0.445345, 0.445991, 0.448085, 0.449637, 0.450339, 0.453482, 0.454775, 0.456139, 0.458766, 0.459777, 0.461669, 0.463841, 0.464765, 0.467362, 0.46903, 0.470428, 0.473236, 0.474219, 0.475866, 0.478265, 0.479343, 0.481874, 0.483983, 0.485173, 0.488129, 0.489759, 0.491443, 0.494142, 0.495331, 0.497667, 0.500125, 0.501265, 0.504214, 0.506247, 0.507659, 0.51031, 0.511974, 0.514132, 0.516882, 0.51815, 0.520937, 0.523536, 0.525125, 0.52798, 0.530019, 0.53199, 0.534845, 0.536503, 0.539055, 0.541797, 0.543271, 0.546195, 0.548757, 0.550623, 0.553543]
CDibad = [0.00588463, 0.00587591, 0.00587341, 0.00590396, 0.00588026, 0.00587782, 0.00588183, 0.00585698, 0.00585025, 0.0058566, 0.00585454, 0.005875, 0.00586265, 0.00586892, 0.00589654, 0.00588845, 0.00588124, 0.00588864, 0.00588322, 0.00589091, 0.00588995, 0.00588725, 0.00590671, 0.00589429, 0.00588971, 0.00589854, 0.00589212, 0.00589244, 0.00589664, 0.00589763, 0.00591354, 0.00591262, 0.00591137, 0.00592136, 0.00591176, 0.00591071, 0.00591893, 0.00591281, 0.00591794, 0.00591655, 0.00591542, 0.00592514, 0.00591515, 0.0059112, 0.00592028, 0.00591486, 0.00591883, 0.00592354, 0.0059209, 0.00592647, 0.00591972, 0.00591554, 0.00591953, 0.00591197, 0.00591032, 0.00591382, 0.0059059, 0.00590737, 0.00590384, 0.00589778, 0.00589761, 0.00588931, 0.00588846, 0.00589156, 0.00588417, 0.00588266, 0.00587735, 0.00586829, 0.00586781, 0.00585689, 0.00584926, 0.00584727, 0.00583654, 0.00583216, 0.00582569, 0.00581349, 0.00580967, 0.00579741, 0.00578652, 0.00578435, 0.00577362, 0.0057639, 0.00575583, 0.00574088, 0.00573432, 0.00572027, 0.00570297, 0.00569357, 0.00567768, 0.00566346, 0.00565299, 0.0056353, 0.00562337, 0.00560788, 0.00558971, 0.00557922, 0.00555981, 0.00554011, 0.00552628, 0.00550529]

rc("figure", figsize=(3.5, 2.6))
rc("font", size=10.0)  #, family="CMU Serif")
rc("lines", linewidth=1.5)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
rc("figure.subplot", left=0.18, bottom=0.16, top=0.97, right=0.95)
rc("axes", color_cycle=["348ABD", "A60628", "009E73"])

figure()
plot(dratio, CLbad)
plot(dratio, CLvec, linewidth=1)
xlabel("diameter/semispan")
ylabel(L"C_L")
xlim([0, 0.3])
ylim([0.38, 0.44])
legend(["original", "smooth"], loc="upper left")
# savefig("../MDAO_paper/figures/smoothCL.pdf", transparent=true)

figure()
plot(dratio, CDibad*10000)
plot(dratio, CDivec*10000, linewidth=1)
xlabel("diameter/semispan")
ylabel(L"$C_{Di}$ (counts)")
xlim([0, 0.8])
legend(["original", "smooth"])
# savefig("../MDAO_paper/figures/smoothCDi.pdf", transparent=true)

end


text:
struct CliqueGraph
    membership_mat::SparseMatrixCSC{Bool, Int}  # nodes x num_cliques
    edge_mat::Matrix{Float64}                   # num_cliques x num_cliques
    # TODO: change to bitset?
    active_cliques::Set{Int}
end


text:
function ConfidenceInterval(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return ConfidenceInterval((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_confidence_level(obj::ConfidenceInterval)
    return jcall(obj, "getConfidenceLevel", jdouble, ())
end

function get_lower_bound(obj::ConfidenceInterval)
    return jcall(obj, "getLowerBound", jdouble, ())
end

function get_upper_bound(obj::ConfidenceInterval)
    return jcall(obj, "getUpperBound", jdouble, ())
end

function to_string(obj::ConfidenceInterval)
    return jcall(obj, "toString", JString, ())
end


text:
const ykow = [1.957e-1, 1.947e-1, 1.735e-1, 1.600e-1, 8.44e-2, 6.27e-2,
              4.56e-2,  3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2]

const ukow = [4.0e0, 2.0e0, 1.0e0, 5.0e-1, 2.5e-1, 1.67e-1,
              1.25e-1, 1.0e-1, 8.33e-2, 7.14e-2, 6.25e-2]

const kowalik_osborne = let res_init=zeros(11), jac_init=zeros(11,4), x_init=[.25, .39, .415, .39]

    function res(x, r)
        for i = 1:11
            r[i] = ykow[i] - x[1]*(ukow[i] ^ 2 + x[2] * ukow[i]) / (ukow[i] ^ 2 + x[3]*ukow[i] + x[4])
        end
        return r
    end

    function jac(x, j)
        for i = 1:11
            t = -1.e0 / (ukow[i]^2 + x[3] * ukow[i] + x[4])
            j[i, 1] = t * (ukow[i]^2 + x[2] * ukow[i])
            j[i, 2] = x[1] * ukow[i]*t
            t = t * j[i, 1] * x[1]
            j[i, 3] = ukow[i] * t
            j[i, 4] = t
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("kowalik_osborne", res, jac,
                                                            res_init, jac_init, x_init;
                                                            scale=scale, verbose=verbose,
                                                            print_steps=print_steps)
end


text:
using Dates
using Statistics

function nse(Qobserved::Array{Float64, 1}, Qmodelled::Array{Float64, 1})
    # input is array of modelled and observed data
    #QobservedAverage = sum(Qobserved) / length(Qobserved) #float
    QobservedAverage = ones(length(Qobserved)) * mean(Qobserved)
    Nominator = sum((Qmodelled - Qobserved).^2)
    Denominator = sum((Qobserved - QobservedAverage).^2)
    NashSutcliffe = 1 - (Nominator / Denominator)
    return NashSutcliffe::Float64
end

function lognse(Qobserved::Array{Float64, 1}, Qmodelled::Array{Float64, 1})
    #QobservedAverage = sum(Qobserved) / length(Qobserved) #float
    QobservedAverage = ones(length(Qobserved)) * mean(Qobserved) # average as array
    Nominator = sum((log.(Qobserved)-log.(Qmodelled)).^2)
    Denominator = sum((log.(Qobserved) - log.(QobservedAverage)).^2)
    NashSutcliffelog = 1 - (Nominator / Denominator)
    return NashSutcliffelog::Float64
end

function volumetricefficiency(Qobserved::Array{Float64, 1}, Qmodelled::Array{Float64, 1})
    Nominator = sum(abs.(Qmodelled - Qobserved))
    Denominator = sum(Qobserved)
    Ve = 1 - Nominator / Denominator
    return Ve::Float64
end

function flowdurationcurve(Q::Array{Float64, 1})
    # input as array, the discharge has to be sorted from largest to smallest and assigned value
    # check function with data of former exercises
    SortedQ = sort(Q, rev = true)
    Rank = collect(1 : length(Q))
    #WRONG IT should be current rank/ (total number of events (Q) +1)
    Exceedanceprobability = Rank ./ (length(SortedQ) .+ 1)
    # exceedence probability should not be higher than 1???
    return SortedQ::Array{Float64, 1}, Exceedanceprobability::Array{Float64, 1}
end

# function autocorrelation(Q::Array{Float64, 1}, Timelag::Int64)
#     # the function shifts the normal discharge by a timelag of x days
#     # function from Euser et al 2013
#     Qshort = Q[1 : end - Timelag]
#     Sum_Nominator = 0
#     Sum_Denominator = 0
#     for (i, Current_Discharge) in enumerate(Qshort)
#         Nominator = (Current_Discharge - mean(Qshort)) * (Q[i + Timelag] - mean(Qshort))
#         Denominator = (Current_Discharge - mean(Qshort))^2
#         Sum_Nominator += Nominator
#         Sum_Denominator += Denominator
#     end
#     AC = Sum_Nominator / Sum_Denominator
#     return AC::Float64
# end

function autocorrelation(Q::Array{Float64, 1}, Timelag::Int64)
    Qshifted = Q[1 + Timelag: end]
    Q = Q[1 : end - Timelag]
    Q_average = ones(length(Q)) * mean(Q)
    Nominator = sum((Q -  Q_average) .* (Qshifted - Q_average))
    Denominator = sum((Q - Q_average).^2)
    AC = Nominator / Denominator
    return AC::Float64
end

# function autocorrelationcurve(Q::Array{Float64, 1}, Timeshift::Int64)
#     # this function calculates an autocorrelation curve by shifting the discharge by a day until reaching the timelag
#     # than the autocorrelation value of each time lag is calculated and is plotted against days of the lag
#     # the first entry to evaluate is the entry after the timelag, so for example 31 if timelage=30
#     Q_normal = Q[Timeshift + 1: end]
#     AC = Float64[]
#     for i in 0 : Timeshift
#         Qshifted = Q[Timeshift + 1 - i : end - i]
#         @assert length(Qshifted) == length(Q_normal)
#         Correlation = cor(Q_normal, Qshifted)
#         push!(AC, Correlation)
#     end
#     @assert length(AC) == Timeshift + 1
#     Lags = collect(0: Timeshift)
#     return AC::Array{Float64, 1}, Lags::Array{Int64, 1}
#
#
# end

function autocorrelationcurve(Q::Array{Float64, 1}, Timelag::Int64)
    Q_normal = Q[1: end - Timelag]
    Q_average = ones(length(Q_normal)) * mean(Q_normal)
    AC = Float64[]
    for i in 0 : Timelag
        Qshifted = Q[1 + i : end - Timelag + i]
        @assert length(Qshifted) == length(Q_normal)
        Nominator = sum((Q_normal -  Q_average) .* (Qshifted - Q_average))
        Denominator = sum((Q_normal - Q_average).^2)
        Correlation = Nominator / Denominator
        push!(AC, Correlation)
    end
    @assert length(AC) == Timelag + 1
    Lags = collect(0: Timelag)
    return AC::Array{Float64, 1}, Lags::Array{Int64, 1}
end

function monthlyrunoff(Area, Precipitation::Array{Float64, 1}, Discharge::Array{Float64, 1}, Timeseries::Array{Date,1})
    # function calculates the monthly runoff coefficient of each month in the timeseries
    # discharge is given in m3/s and precipitation in mm/d
    # convert Discharge to mm/d
    #Discharge = Discharge / Area * (1000 * 3600 * 24)
    sum_Precipitation = 0
    sum_Discharge = 0
    monthly_Runoff = Float64[]
    Month = Date[]
    for (i, current_date) in enumerate(Timeseries)
        # the precipitation and discharges are summed up
        sum_Precipitation += Precipitation[i]
        sum_Discharge += Discharge[i]
        #current_date = Date(current_date, dateformat"y,m,d")
        day = Dates.day(current_date)
        # if the last day of the month is reached
        if day == Dates.daysinmonth(current_date)
            if sum_Precipitation != 0
                monthly_Runoff_Current = sum_Discharge / sum_Precipitation
            else
                monthly_Runoff_Current = sum_Discharge / 0.01
            end
            push!(monthly_Runoff, monthly_Runoff_Current)
            sum_Precipitation = 0
            sum_Discharge = 0
            #Month_Current = Dates.format(current_date, "yyyy-mm")
            push!(Month, current_date)
        end
    end
    return monthly_Runoff::Array{Float64,1}, Month::Array{Date,1}
end

function averagemonthlyrunoff(Area::Float64, Precipitation::Array{Float64, 1}, Discharge::Array{Float64, 1}, Timeseries::Array{Date,1})
    # calculates the average runoff coefficient of a certain month over the timeperiod
    # all runoff coefficients of a certain months are summed up
    # assert that
    monthly_Runoff, Month = monthlyrunoff(Area, Precipitation, Discharge, Timeseries)

    @assert length(monthly_Runoff) == length(Month)
    sum_January = 0
    sum_February = 0
    sum_March = 0
    sum_April = 0
    sum_Mai = 0
    sum_June = 0
    sum_July = 0
    sum_August = 0
    sum_Sept = 0
    sum_Okt = 0
    sum_Nov = 0
    sum_Dec = 0
    for (i, current_date) in enumerate(Month)
        #current_date = Date(current_date, dateformat"yyyy-mm")
        current_month = Dates.month(current_date)
        if current_month == 1
            sum_January += monthly_Runoff[i]
        elseif current_month == 2
            sum_February += monthly_Runoff[i]
        elseif current_month == 3
            sum_March += monthly_Runoff[i]
        elseif current_month == 4
            sum_April += monthly_Runoff[i]
        elseif current_month == 5
            sum_Mai += monthly_Runoff[i]
        elseif current_month == 6
            sum_June += monthly_Runoff[i]
        elseif current_month == 7
            sum_July += monthly_Runoff[i]
        elseif current_month == 8
            sum_August += monthly_Runoff[i]
        elseif current_month == 9
            sum_Sept += monthly_Runoff[i]
        elseif current_month == 10
            sum_Okt += monthly_Runoff[i]
        elseif current_month == 11
            sum_Nov += monthly_Runoff[i]
        else
            sum_Dec += monthly_Runoff[i]
        end
    end
    sum_monthly_Runoff = [sum_January, sum_February, sum_March, sum_April, sum_Mai, sum_June, sum_July, sum_August, sum_Sept, sum_Okt, sum_Nov, sum_Dec]
    FirstYear = Dates.year(Month[1])
    LastYear = Dates.year(Month[end])
    Number_Years = LastYear - FirstYear
    @assert Number_Years * 12 == length(monthly_Runoff)
    average_monthly_Runoff = sum_monthly_Runoff / Number_Years
    return average_monthly_Runoff::Vector{Float64}
end

function snowcover(Modeled_Area::Array{Float64, 1}, Observed_Area::Array{Float64, 1})
    # observed data can be -1 which is the error value
    index = findall(x-> x>= 0, Observed_Area)
        Difference = (ones(length(index)) - abs.(Modeled_Area[index] - Observed_Area[index]))
    Mean_Difference = mean(Difference)
    return Mean_Difference::Float64
end

export nse
export lognse
export volumetricefficiency
export flowdurationcurve

function objectivefunctions(Modelled_Discharge::Array{Float64, 1}, Snow_Cover::Float64, Observed_Discharge::Array{Float64, 1}, observed_FDC::Array{Float64, 1}, observed_AC_1day::Float64, observed_AC_90day::Array{Float64, 1}, observed_monthly_runoff::Array{Float64, 1}, Area::Float64, Precipitation::Array{Float64, 1}, Timerseries::Array{Date, 1})
    # calculate the nse, lognse, ve
    NSE = nse(Observed_Discharge, Modelled_Discharge)
    NSElog = lognse(Observed_Discharge, Modelled_Discharge)
    VE = volumetricefficiency(Observed_Discharge, Modelled_Discharge)
    # calculate the flow duration durves
    modelled_FDC = flowdurationcurve(log.(Modelled_Discharge))
    NSE_FDC = nse(observed_FDC, modelled_FDC[1])
    #calculate the autocorrelation curves
    modelled_AC_1day = autocorrelation(Modelled_Discharge, 1)
    modelled_AC_90day = autocorrelationcurve(Modelled_Discharge, 90)
    Reative_Error_AC_1day = 1.0 - abs.(observed_AC_1day - modelled_AC_1day)/ observed_AC_1day
    NSE_AC_90day = nse(observed_AC_90day, modelled_AC_90day[1])
    #calculate the monthly runoff
    #area of whole catchment, precipitation whole catchment
    #timeseries as dates
    Timeseries_Runoff = Timerseries
    modelled_monthly_runoff = monthlyrunoff(Area, Precipitation, Modelled_Discharge, Timeseries_Runoff)
    # calculate the NSE of the monthly runoffs
    NSE_monthly_runoff = nse(observed_monthly_runoff, modelled_monthly_runoff[1])
    #Relative_Error_Runoff = ones(length(observed_average_runoff)) - abs.(observed_average_runoff - modelled_average_runoff) ./ observed_average_runoff
    #Relative_Error_Runoff = mean(Relative_Error_Runoff)
    #snow cover was already calculated for each elevation zone
    # function for snow cover should be maximized
    if NSE > 0 && NSElog > 0 && NSE_FDC > 0 && NSE_AC_90day > 0 && NSE_monthly_runoff > 0
        #store values
        ObjFunctions = [NSE, NSElog, VE, NSE_FDC, Reative_Error_AC_1day, NSE_AC_90day, NSE_monthly_runoff, Snow_Cover]
        Sum = 0
        for Obj in ObjFunctions
            Sum+= (1 - Obj)^2
        end
        Euclidean_Distance = (Sum / length(ObjFunctions))^0.5
    else
        Euclidean_Distance = -9999.0
        ObjFunctions = [-9999.0]
    end

    # volumetric efficiency, NSE should be maximized

    # Euclidean_Distance = ((1-NSE)^2 + (1 - NSElog)^2 + (1 - VE)^2 + (1 - NSE_FDC)^2 + (1 - Reative_Error_AC_1day)^2 + (1 - NSE_AC_90day)^2) + (1 - Reative_Error_Runoff)^2 + (1 - Snow_Cover)^2
    # Euclidean_Distance = Euclidean_Distance / 8

    return Euclidean_Distance::Float64, ObjFunctions::Array{Float64, 1}
end

function objectivefunctions_delete_days(Modelled_Discharge_AC::Array{Float64, 1}, Modelled_Discharge::Array{Float64, 1}, Snow_Cover::Float64, Observed_Discharge::Array{Float64, 1}, observed_FDC::Array{Float64, 1}, observed_AC_1day::Float64, observed_AC_90day::Array{Float64, 1}, observed_monthly_runoff::Array{Float64, 1}, Area::Float64, Precipitation::Array{Float64, 1}, Timerseries::Array{Date, 1})
    # calculate the nse, lognse, ve
    NSE = nse(Observed_Discharge, Modelled_Discharge)
    NSElog = lognse(Observed_Discharge, Modelled_Discharge)
    VE = volumetricefficiency(Observed_Discharge, Modelled_Discharge)
    # calculate the flow duration durves
    modelled_FDC = flowdurationcurve(log.(Modelled_Discharge))
    NSE_FDC = nse(observed_FDC, modelled_FDC[1])
    #calculate the autocorrelation curves
    modelled_AC_1day = autocorrelation(Modelled_Discharge_AC, 1)
    modelled_AC_90day = autocorrelationcurve(Modelled_Discharge_AC, 90)
    Reative_Error_AC_1day = 1.0 - abs.(observed_AC_1day - modelled_AC_1day)/ observed_AC_1day
    NSE_AC_90day = nse(observed_AC_90day, modelled_AC_90day[1])
    #calculate the monthly runoff
    #area of whole catchment, precipitation whole catchment
    #timeseries as dates
    Timeseries_Runoff = Timerseries
    modelled_monthly_runoff = monthlyrunoff(Area, Precipitation, Modelled_Discharge, Timeseries_Runoff)
    # calculate the NSE of the monthly runoffs
    NSE_monthly_runoff = nse(observed_monthly_runoff, modelled_monthly_runoff[1])
    #Relative_Error_Runoff = ones(length(observed_average_runoff)) - abs.(observed_average_runoff - modelled_average_runoff) ./ observed_average_runoff
    #Relative_Error_Runoff = mean(Relative_Error_Runoff)
    #snow cover was already calculated for each elevation zone
    # function for snow cover should be maximized
    if NSE > 0 && NSElog > 0 && NSE_FDC > 0 && NSE_AC_90day > 0 && NSE_monthly_runoff > 0
        #store values
        ObjFunctions = [NSE, NSElog, VE, NSE_FDC, Reative_Error_AC_1day, NSE_AC_90day, NSE_monthly_runoff, Snow_Cover]
        Sum = 0
        for Obj in ObjFunctions
            Sum+= (1 - Obj)^2
        end
        Euclidean_Distance = (Sum / length(ObjFunctions))^0.5
    else
        Euclidean_Distance = -9999.0
        ObjFunctions = [-9999.0]
    end

    # volumetric efficiency, NSE should be maximized

    # Euclidean_Distance = ((1-NSE)^2 + (1 - NSElog)^2 + (1 - VE)^2 + (1 - NSE_FDC)^2 + (1 - Reative_Error_AC_1day)^2 + (1 - NSE_AC_90day)^2) + (1 - Reative_Error_Runoff)^2 + (1 - Snow_Cover)^2
    # Euclidean_Distance = Euclidean_Distance / 8

    return Euclidean_Distance::Float64, ObjFunctions::Array{Float64, 1}
end

function objectivefunctions_projections(Modelled_Discharge::Array{Float64, 1}, Snow_Cover::Float64, Observed_Discharge::Array{Float64, 1}, observed_FDC::Array{Float64, 1}, observed_AC_1day::Float64, observed_AC_90day::Array{Float64, 1}, observed_monthly_runoff::Array{Float64, 1}, Area::Float64, Precipitation::Array{Float64, 1}, Timerseries::Array{Date, 1})
    # calculate the nse, lognse, ve
    NSE = nse(Observed_Discharge, Modelled_Discharge)
    NSElog = lognse(Observed_Discharge, Modelled_Discharge)
    VE = volumetricefficiency(Observed_Discharge, Modelled_Discharge)
    # calculate the flow duration durves
    modelled_FDC = flowdurationcurve(log.(Modelled_Discharge))
    NSE_FDC = nse(observed_FDC, modelled_FDC[1])
    #calculate the autocorrelation curves
    modelled_AC_1day = autocorrelation(Modelled_Discharge, 1)
    modelled_AC_90day = autocorrelationcurve(Modelled_Discharge, 90)
    Reative_Error_AC_1day = 1.0 - abs.(observed_AC_1day - modelled_AC_1day)/ observed_AC_1day
    NSE_AC_90day = nse(observed_AC_90day, modelled_AC_90day[1])
    #calculate the monthly runoff
    #area of whole catchment, precipitation whole catchment
    #timeseries as dates
    Timeseries_Runoff = Timerseries
    modelled_monthly_runoff = monthlyrunoff(Area, Precipitation, Modelled_Discharge, Timeseries_Runoff)
    #print("runoff", size(modelled_monthly_runoff[1]), "\n")
    @assert modelled_monthly_runoff[1] >= zeros(length(modelled_monthly_runoff[1]))
    # calculate the NSE of the monthly runoffs
    #print(observed_monthly_runoff)
    #print("mean1", mean(observed_monthly_runoff),"\n")
    NSE_monthly_runoff = nse(observed_monthly_runoff, modelled_monthly_runoff[1])
    # function for snow cover should be maximized
    ObjFunctions = [NSE, NSElog, VE, NSE_FDC, Reative_Error_AC_1day, NSE_AC_90day, NSE_monthly_runoff, Snow_Cover]
    Sum = 0
    for Obj in ObjFunctions
        Sum+= (1 - Obj)^2
    end
    Euclidean_Distance = (Sum / length(ObjFunctions))^0.5
    return Euclidean_Distance::Float64, ObjFunctions::Array{Float64, 1}
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
"""
	function σTemp(ω0::Real,temp::Number)

Returns the standard deviation due to Doppler broadening of frequency ω0 at temp
"""
function σTemp(ω0::Real,temp::Number)
	kbOverMhC2 = 9.178e-14;
	return sqrt(kbOverMhC2*temp)*ω0
end

export σTemp


text:
#=
urlify:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-27
=#

module URLify
    # parameters
    const max_len = 1000

    # methods
    function replace_space(str)
        # remove leading and trailing spaces
        str = strip(str, [' '])

        # check current lenght and count spaces
        len_str = length(str)
        spc_count = count(i->(i==' '), str)

        # calculate new length
        new_len = len_str + spc_count * 2
        if new_len > max_len
            println("New lenght is longer than maximum")
            return -1
        end

        # start index filling char from end
        index = new_len

        str_sub = split(str, "")

        # fill string array
        for i in (len_str-2):(new_len-3)
            push!(str_sub, "0")
        end

        # fill rest of string from end
        for i in (len_str):-1:1
            if str_sub[i] == " "
                str_sub[index] = "0"
                str_sub[index - 1] = "2"
                str_sub[index - 2] = "%"
                index -= 3
            else
                str_sub[index] = str_sub[i]
                index -= 1
            end
            println("$(i), $(str_sub[i]), $(str_sub)")
        end

        return join(str_sub)
    end

    function main()
        println(replace_space(" Mr John Smith "))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .URLify
    URLify.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
# # Zombie Outbreak
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/agents/zombies.mp4?raw=true" type="video/mp4">
# </video>
# ```
#
# This model showcases an ABM running on a map, using [`OpenStreetMapSpace`](@ref).
#
# ## Constructing the end of days
using Agents

# We'll simulate a zombie outbreak in a city. To do so, we start with an agent which
# satisfies the OSMSpace conditions of having a `pos`ition of type
# `Tuple{Int,Int,Float64}`. For simplicity though we shall build this with the [`@agent`](@ref)
# macro.

@agent Zombie OSMAgent begin
    infected::Bool
end

# To be explicit, this macro builds the following type:
# ```julia
# mutable struct Zombie <: AbstractAgent
#     id::Int
#     pos::Tuple{Int,Int,Float64}
#     infected::Bool
# end
# ```
# where a tuple `(i, j, x)::Tuple{Int,Int,Float64}` means a position
# on the road between nodes `i, j` of the map, having progressed `x` distance along the road.

# The model constructor we build consists of a map, and 100 agents scattered randomly
# around it. They have their own agenda and need to travel to some new destination.
# Unfortunately one of the population has turned and will begin infecting anyone who
# comes close.

function initialise(; map_path = OSM.test_map())
    model = ABM(Zombie, OpenStreetMapSpace(map_path))

    for id in 1:100
        start = random_position(model) # At an intersection
        human = Zombie(id, start, false)
        add_agent_pos!(human, model)
        OSM.random_route!(human, model; limit = 50) # try 50 times to find a random route
    end
    ## We'll add patient zero at a specific (latitude, longitude)
    start = OSM.road((51.5328328, 9.9351811), model)
    finish = OSM.intersection((51.530876112711745, 9.945125635913511), model)
    zombie = add_agent!(start, model, true)
    plan_route!(zombie, finish, model)
    ## This function call creates & adds an agent, see `add_agent!`
    return model
end

# In our model, zombies are seemingly oblivious to their state, since they keep going about their
# business, but start eating people along the way. Perhaps they can finally express their distaste
# for city commuting.

function agent_step!(agent, model)
    ## Each agent will progress slightly along their route
    move_along_route!(agent, model, 0.005)

    if is_stationary(agent, model) && rand(model.rng) < 0.1
        ## When stationary, give the agent a 10% chance of going somewhere else
        OSM.random_route!(agent, model; limit = 50)
        ## Start on new route
        move_along_route!(agent, model, 0.0001)
    end

    if agent.infected
        ## Agents will be infected if they get too close to a zombie.
        map(i -> model[i].infected = true, nearby_ids(agent, model, 0.0003))
    end
end

# ## Visualising the fall of humanity
#
# Plotting this space in a seamless manner is a work in progress. For now we
# use [OSMMakie.jl](https://github.com/fbanning/OSMMakie.jl) and
# a custom routine.

# using OSMMakie
using CairoMakie
ac(agent) = agent.infected ? :green : :black
as(agent) = agent.infected ? 6 : 5

model = initialise()

# fig, ax, plot = osmplot(model.space.map)
ids = model.scheduler(model)
colors = Node([ac(model[i]) for i in ids])
sizes = Node([as(model[i]) for i in ids])
marker = :circle
pos = Node(Point2f[OSM.latlon(model[i].pos, model) for i in ids])
fig, _= scatter(pos; color = colors, size = sizes, marker)

record(fig, "outbreak.mp4", 1:200; framerate = 15) do i
    Agents.step!(model, agent_step!, 1)
    ids = model.scheduler(model)
    colors[] = [ac(model[i]) for i in ids]
    sizes[] = [as(model[i]) for i in ids]
    pos[] = Point2f[OSM.latlon(model[i].pos, model) for i in ids]
end
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../outbreak.mp4" type="video/mp4">
# </video>
# ```


text:
ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python")

using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using ArgParse
using FileIO
using JLD2
using RiskSensitiveSAC
import RobotOS.Time
using GR
GR.inline("png") # https://github.com/JuliaPlots/Plots.jl/issues/1723

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--file_path"
            help = "absolute path to the .jld2 file"
            required = true
            range_tester = isfile
        "--process_all"
            help = "make gifs for all .jld2 files in the same directory"
            action = :store_true
        "--filter_str"
            help = "if process_all is true, apply only to files with fiter_str"
            default = nothing
        "--dtplot"
            help = "plot delta t [s]"
            arg_type = Float64
            required = true
        "--fps"
            help = "frame per second"
            arg_type = Int64
            required = true
        "--xmin"
            help = "minimum x value [m]. default = 0.0"
            arg_type = Float64
            default = 0.0;
        "--xmax"
            help = "maximum x value [m]. default = 20.0"
            arg_type = Float64
            default = 20.0;
        "--ymin"
            help = "minimum y value [m]. default = 0.0"
            arg_type = Float64
            default = 0.0;
        "--ymax"
            help = "maximum y value [m]. default = 20.0"
            arg_type = Float64
            default = 20.0;
        "--height"
            help = "height of the figure. default = 400"
            arg_type = Int64
            default = 400
        "--width"
            help = "width of the figure. default = 600"
            arg_type = Int64
            default = 600
        "--legend"
            help = "legend position. default = topright"
            arg_type = String
            default = "topright"
        "--legendfontsize"
            help = "legend font size. default = 7"
            arg_type = Int64
            default = 7
        "--markersize"
            help = "markersize. default = 8.0"
            arg_type = Float64
            default = 8.0
        "--show_prediction"
            help = "show predictions."
            action = :store_true
        "--show_nominal_trajectory"
            help = "show nominal ego trajectory."
            action = :store_true
        "--show_past_trajectory"
            help = "show past ego trajectory."
            action = :store_true
        "--dummy_pos"
            help = "coordinate to omit from prediction plots."
            arg_type = Float64

    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    file_path = parsed_args["file_path"];
    filter_str = parsed_args["filter_str"];
    if parsed_args["process_all"]
        file_dir = splitdir(file_path)[1];
        files = readdir(file_dir)
        if isnothing(filter_str)
            valid_files = filter(file -> occursin(".jld2", file), files);
        else
            filter_str = parsed_args["filter_str"];
            valid_files = filter((file -> occursin(".jld2", file) && occursin(filter_str, file)), files);
        end
        ii = 1;
        for file in valid_files
            println("Progress: $(ii)/$(length(valid_files))")
            valid_file_path = joinpath(file_dir, file);
            file_content = load(valid_file_path);
            result = file_content["result"];
            gif_file_path = replace(valid_file_path, ".jld2" => ".gif");
            make_gif(result,
                    dtplot=parsed_args["dtplot"],
                    fps=parsed_args["fps"],
                    figsize=(parsed_args["width"], parsed_args["height"]),
                    legendfontsize=parsed_args["legendfontsize"],
                    legend=Symbol(parsed_args["legend"]),
                    xlim=(parsed_args["xmin"], parsed_args["xmax"]),
                    ylim=(parsed_args["ymin"], parsed_args["ymax"]),
                    markersize=parsed_args["markersize"],
                    filename=gif_file_path,
                    show_prediction=parsed_args["show_prediction"],
                    show_nominal_trajectory=parsed_args["show_nominal_trajectory"],
                    show_past_ego_trajectory=parsed_args["show_past_trajectory"],
                    dummy_pos=parsed_args["dummy_pos"]);
            ii += 1;
        end
    else
        file_content = load(file_path);
        result = file_content["result"];
        gif_file_path = replace(file_path, ".jld2" => ".gif");
        make_gif(result,
                dtplot=parsed_args["dtplot"],
                fps=parsed_args["fps"],
                figsize=(parsed_args["width"], parsed_args["height"]),
                legendfontsize=parsed_args["legendfontsize"],
                legend=Symbol(parsed_args["legend"]),
                xlim=(parsed_args["xmin"], parsed_args["xmax"]),
                ylim=(parsed_args["ymin"], parsed_args["ymax"]),
                markersize=parsed_args["markersize"],
                filename=gif_file_path,
                show_prediction=parsed_args["show_prediction"],
                show_nominal_trajectory=parsed_args["show_nominal_trajectory"],
                show_past_ego_trajectory=parsed_args["show_past_trajectory"],
                dummy_pos=parsed_args["dummy_pos"]);
    end
end

main()


text:
using Base.Test

print_with_color(:blue, "Running tests:\n")

tests = (
    "modes.jl",
    "step_size.jl",
    "derivative.jl",
    "second_derivative.jl",
    "gradient.jl",
    "hessian.jl",
)

@testset "All FiniteDiff tests" begin
    for t in tests
        try
            include(t)
            print_with_color(:green, @sprintf("* %s\n", t))
        catch
            print_with_color(:red, @sprintf("* %s\n", t))
        end
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct LCLFilter <: Filter
        lf::Float64
        rf::Float64
        cf::Float64
        lg::Float64
        rg::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a LCL filter outside the converter

# Arguments
- `lf::Float64`: Series inductance in p.u. of converter filter, validation range: (0, nothing)
- `rf::Float64`: Series resistance in p.u. of converter filter, validation range: (0, nothing)
- `cf::Float64`: Shunt capacitance in p.u. of converter filter, validation range: (0, nothing)
- `lg::Float64`: Series inductance in p.u. of converter filter to the grid, validation range: (0, nothing)
- `rg::Float64`: Series resistance in p.u. of converter filter to the grid, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct LCLFilter <: Filter
    "Series inductance in p.u. of converter filter"
    lf::Float64
    "Series resistance in p.u. of converter filter"
    rf::Float64
    "Shunt capacitance in p.u. of converter filter"
    cf::Float64
    "Series inductance in p.u. of converter filter to the grid"
    lg::Float64
    "Series resistance in p.u. of converter filter to the grid"
    rg::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function LCLFilter(lf, rf, cf, lg, rg, ext=Dict{String, Any}(), )
    LCLFilter(lf, rf, cf, lg, rg, ext, [:id_c, :iq_c, :vd_cap, :vq_cap, :id_o, :iq_o], 6, InfrastructureSystemsInternal(), )
end

function LCLFilter(; lf, rf, cf, lg, rg, ext=Dict{String, Any}(), )
    LCLFilter(lf, rf, cf, lg, rg, ext, )
end

# Constructor for demo purposes; non-functional.
function LCLFilter(::Nothing)
    LCLFilter(;
        lf=0,
        rf=0,
        cf=0,
        lg=0,
        rg=0,
        ext=Dict{String, Any}(),
    )
end

"""Get LCLFilter lf."""
get_lf(value::LCLFilter) = value.lf
"""Get LCLFilter rf."""
get_rf(value::LCLFilter) = value.rf
"""Get LCLFilter cf."""
get_cf(value::LCLFilter) = value.cf
"""Get LCLFilter lg."""
get_lg(value::LCLFilter) = value.lg
"""Get LCLFilter rg."""
get_rg(value::LCLFilter) = value.rg
"""Get LCLFilter ext."""
get_ext(value::LCLFilter) = value.ext
"""Get LCLFilter states."""
get_states(value::LCLFilter) = value.states
"""Get LCLFilter n_states."""
get_n_states(value::LCLFilter) = value.n_states
"""Get LCLFilter internal."""
get_internal(value::LCLFilter) = value.internal


text:
#using Pkg
# pkg"add Oceananigans CairoMakie JLD2"
ENV["GKSwstype"] = "100"
# pushfirst!(LOAD_PATH, @__DIR__)
# pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "..")) # add Oceananigans

using Printf
using Statistics
using JLD2

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
using Oceananigans.Grids: xnode, ynode, znode

using Random
Random.seed!(1234)

filename = "zonally_averaged_channel_withGM"

# Architecture
architecture = CPU()

# Domain
const Lx = 1000kilometers # zonal domain length [m]
const Ly = 2000kilometers # meridional domain length [m]
const Lz = 2kilometers    # depth [m]

# number of grid points
Ny = 256
Nz = 80

save_fields_interval = 7days
stop_time = 5years
stop_time = 70 days
Δt₀ = 5minutes

# stretched grid

# we implement here a linearly streched grid in which the top grid cell has Δzₜₒₚ
# and every other cell is bigger by a factor σ, e.g.,
# Δzₜₒₚ, Δzₜₒₚ * σ, Δzₜₒₚ * σ², ..., Δzₜₒₚ * σᴺᶻ⁻¹,
# so that the sum of all cell heights is Lz

# Given Lz and stretching factor σ > 1 the top cell height is Δzₜₒₚ = Lz * (σ - 1) / σ^(Nz - 1)

# σ = 1.1 # linear stretching factor
# Δz_center_linear(k) = Lz * (σ - 1) * σ^(Nz - k) / (σ^Nz - 1) # k=1 is the bottom-most cell, k=Nz is the top cell
# linearly_spaced_faces(k) = k==1 ? -Lz : - Lz + sum(Δz_center_linear.(1:k-1))

# We choose a regular grid though because of numerical issues that yet need to be resolved
grid = RectilinearGrid(architecture,
                       topology = (Flat, Bounded, Bounded),
                       size = (Ny, Nz),
                       halo = (3, 3),
                       y = (0, Ly),
                       z = (-Lz, 0))

# The vertical spacing versus depth for the prescribed grid
#=
plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
     axis=(xlabel = "Vertical spacing (m)",
           ylabel = "Depth (m)"))
=#

@info "Built a grid: $grid."

#####
##### Boundary conditions
#####

α  = 2e-4     # [K⁻¹] thermal expansion coefficient
g  = 9.8061   # [m s⁻²] gravitational constant
cᵖ = 3994.0   # [J K⁻¹] heat capacity
ρ  = 1024.0   # [kg m⁻³] reference density

parameters = (Ly = Ly,
              Lz = Lz,
              Qᵇ = 10 / (ρ * cᵖ) * α * g,          # buoyancy flux magnitude [m² s⁻³]
              y_shutoff = 5/6 * Ly,                # shutoff location for buoyancy flux [m]
              τ = 0.2/ρ,                           # surface kinematic wind stress [m² s⁻²]
              μ = 1 / 30days,                      # bottom drag damping time-scale [s⁻¹]
              ΔB = 8 * α * g,                      # surface vertical buoyancy gradient [s⁻²]
              H = Lz,                              # domain depth [m]
              h = 1000.0,                          # exponential decay scale of stable stratification [m]
              y_sponge = 19/20 * Ly,               # southern boundary of sponge layer [m]
              λt = 7.0days                         # relaxation time scale [s]
)

@inline function buoyancy_flux(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return ifelse(y < p.y_shutoff, p.Qᵇ * cos(3π * y / p.Ly), 0.0)
end

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, discrete_form=true, parameters=parameters)

@inline function u_stress(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return - p.τ * sin(π * y / p.Ly)
end

u_stress_bc = FluxBoundaryCondition(u_stress, discrete_form=true, parameters=parameters)

@inline u_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.u[i, j, 1]
@inline v_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.v[i, j, 1]

u_drag_bc = FluxBoundaryCondition(u_drag, discrete_form=true, parameters=parameters)
v_drag_bc = FluxBoundaryCondition(v_drag, discrete_form=true, parameters=parameters)

b_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc)

u_bcs = FieldBoundaryConditions(top = u_stress_bc, bottom = u_drag_bc)
v_bcs = FieldBoundaryConditions(bottom = v_drag_bc)

#####
##### Coriolis
#####

const f = -1e-4     # [s⁻¹]
const β =  1e-11    # [m⁻¹ s⁻¹]
coriolis = BetaPlane(f₀ = f, β = β)

#####
##### Forcing and initial condition
#####

@inline initial_buoyancy(z, p) = p.ΔB * (exp(z / p.h) - exp(-p.Lz / p.h)) / (1 - exp(-p.Lz / p.h))
@inline mask(y, p) = max(0.0, y - p.y_sponge) / (Ly - p.y_sponge)

@inline function buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)
    timescale = p.λt
    y = ynode(Center(), j, grid)
    z = znode(Center(), k, grid)
    target_b = initial_buoyancy(z, p)
    b = @inbounds model_fields.b[i, j, k]
    return - 1 / timescale  * mask(y, p) * (b - target_b)
end

Fb = Forcing(buoyancy_relaxation, discrete_form = true, parameters = parameters)

# Turbulence closures

κh = 0.5e-5 # [m²/s] horizontal diffusivity
νh = 30.0   # [m²/s] horizontal viscocity
κz = 0.5e-5 # [m²/s] vertical diffusivity
νz = 3e-4   # [m²/s] vertical viscocity

horizontal_diffusivity = AnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz)

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0)

catke = CATKEVerticalDiffusivity()

gerdes_koberle_willebrand_tapering = Oceananigans.TurbulenceClosures.FluxTapering(1e-2)

gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = 1000,
                                                                slope_limiter = gerdes_koberle_willebrand_tapering)
#####
##### Model building
#####

@info "Building a model..."

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = ImplicitFreeSurface(),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = coriolis,
                                    #closure = (horizontal_diffusivity, convective_adjustment, gent_mcwilliams_diffusivity),
                                    closure = (catke, horizontal_diffusivity, gent_mcwilliams_diffusivity),
                                    tracers = (:b, :e, :c),
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    forcing = (; b=Fb))

@info "Built $model."

#####
##### Initial conditions
#####

# resting initial condition
ε(σ) = σ * randn()
bᵢ(x, y, z) = parameters.ΔB * ( exp(z/parameters.h) - exp(-Lz/parameters.h) ) / (1 - exp(-Lz/parameters.h)) + ε(1e-8)
uᵢ(x, y, z) = ε(1e-8)
vᵢ(x, y, z) = ε(1e-8)
wᵢ(x, y, z) = ε(1e-8)

Δy = 100kilometers
Δz = 100
Δc = 2Δy
cᵢ(x, y, z) = exp(-y^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2)

set!(model, b=bᵢ, u=uᵢ, v=vᵢ, w=wᵢ, c=cᵢ)

#####
##### Simulation building
#####

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# add timestep wizard callback
wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=20minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# add progress callback
wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

#####
##### Diagnostics
#####

u, v, w = model.velocities
b, c = model.tracers.b, model.tracers.c

dependencies = (gent_mcwilliams_diffusivity,
                b,
                Val(1),
                model.clock,
                model.diffusivity_fields,
                model.tracers,
                model.buoyancy,
                model.velocities)

using Oceananigans.TurbulenceClosures: diffusive_flux_y, diffusive_flux_z, ∇_dot_qᶜ

vb_op  = KernelFunctionOperation{Center, Face, Center}(diffusive_flux_y, grid, architecture=architecture, computed_dependencies=dependencies)
wb_op  = KernelFunctionOperation{Center, Center, Face}(diffusive_flux_z, grid, architecture=architecture, computed_dependencies=dependencies)
∇_q_op = KernelFunctionOperation{Center, Center, Center}(∇_dot_qᶜ, grid, architecture=architecture, computed_dependencies=dependencies)

vb = Field(vb_op)
wb = Field(wb_op)
∇_q = Field(∇_q_op)

outputs = merge(fields(model), (; vb, wb, ∇_q))

# #####
# ##### Build checkpointer and output writer
# #####

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(5years),
                                                        prefix = filename,
                                                        force = true)

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                      schedule = TimeInterval(save_fields_interval),
                                                      prefix = filename,
                                                      field_slicer = nothing,
                                                      verbose = false,
                                                      force = true)

@info "Running the simulation..."

run!(simulation, pickup=false)

#####
##### Visualization
#####

using CairoMakie

filepath = filename * ".jld2"

zonal_file = jldopen(filepath)

grid = zonal_file["serialized/grid"]

xu, yu, zu = nodes((Face, Center, Center), grid)
xv, yv, zv = nodes((Center,Face,  Center), grid)
xw, yw, zw = nodes((Center, Center, Face), grid)
xc, yc, zc = nodes((Center, Center, Center), grid)

u_timeseries = FieldTimeSeries(filepath, "u", grid=grid)
@show umax = maximum(abs, u_timeseries[:, :, :, :])

v_timeseries = FieldTimeSeries(filepath, "v", grid=grid)
@show umax = maximum(abs, v_timeseries[:, :, :, :])

w_timeseries = FieldTimeSeries(filepath, "w", grid=grid)
@show umax = maximum(abs, w_timeseries[:, :, :, :])

b_timeseries = FieldTimeSeries(filepath, "b", grid=grid)
@show b_timeseries

vb_timeseries = FieldTimeSeries(filepath, "vb", grid=grid)
@show vb_timeseries

wb_timeseries = FieldTimeSeries(filepath, "wb", grid=grid)
@show wb_timeseries

wb_timeseries = FieldTimeSeries(filepath, "wb", grid=grid)
@show wb_timeseries

∇_q_timeseries = FieldTimeSeries(filepath, "∇_q", grid=grid)
@show ∇_q_timeseries

_, _, _, nt = size(b_timeseries)

y_limits = (0, 2000)
z_limits = (-200, 0)

which_iteration = nt-1

fig, ax, hm = heatmap(yc * 1e-3, zc, interior(∇_q_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-1e-7, 1e-7),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "∇⋅q", textsize=30)
xlims!(y_limits[1], y_limits[2])
ylims!(z_limits[1], z_limits[2])
display(fig)

fig, ax, hm = heatmap(yv * 1e-3, zv, interior(vb_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-2e-5, 2e-5),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "v'b'", textsize=30)
xlims!(y_limits[1], y_limits[2])
ylims!(z_limits[1], z_limits[2])
display(fig)

fig, ax, hm = heatmap(yw * 1e-3, zw, interior(wb_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-1e-6, 1e-6),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "w'b'", textsize=30)
xlims!(y_limits[1], y_limits[2])
ylims!(z_limits[1], z_limits[2])
display(fig)

b = b_timeseries[which_iteration]

fig, ax, hm = heatmap(yc * 1e-3, zc, interior(b)[1, :, :],
                    #   colormap=:balance, colorrange=(-1e-7, 1e-7),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "b", textsize=30)
xlims!(y_limits[1], y_limits[2])
ylims!(z_limits[1], z_limits[2])
display(fig)

fig, ax, hm = heatmap(yu * 1e-3, zu, interior(u_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-0.5, 0.5),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "u", textsize=30)
# xlims!(y_limits[1], y_limits[2])
# ylims!(z_limits[1], z_limits[2])
display(fig)

fig, ax, hm = heatmap(yv * 1e-3, zv, interior(v_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-0.1, 0.1),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "v", textsize=30)
# xlims!(y_limits[1], y_limits[2])
# ylims!(z_limits[1], z_limits[2])
display(fig)

fig, ax, hm = heatmap(yw * 1e-3, zw, interior(w_timeseries)[1, :, :, which_iteration],
                      colormap=:balance, colorrange=(-1e-4, 1e-4),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "w", textsize=30)
# xlims!(y_limits[1], y_limits[2])
# ylims!(z_limits[1], z_limits[2])
display(fig)

b_z = Field(Center, Center, Face, grid)

b_z .= ∂z(b)

fig, ax, hm = heatmap(yw * 1e-3, zw, interior(b_z)[1, :, :],
                      colormap=:curl, colorrange=(-5e-4, 5e-4),
                      axis=(xlabel="y", ylabel="z", xticklabelsize=20,  yticklabelsize=20, xlabelsize=25, ylabelsize=25))
Colorbar(fig[1, 2], hm, ticklabelsize=20)
fig[0, :] = Label(fig, "∂b/∂z", textsize=30)
xlims!(y_limits[1], y_limits[2])
ylims!(z_limits[1], z_limits[2])
display(fig)

umax = 1
ulims = (-umax, umax) .* 0.8
ulevels = vcat([-umax], range(ulims[1], ulims[2], length=31), [umax])

ylims = (0, grid.Ly) .* 1e-3
zlims = (-grid.Lz, 0)

anim = @animate for i in 1:length(b_timeseries.times)-1
    b = b_timeseries[i]
    u = u_timeseries[i]

    b_z .= ∂z(b)

    b_yz = interior(b)[1, :, :]
    u_yz = interior(u)[1, :, :]

    b_z_yz = interior(b_z)[1, :, :]

    @show bmax = max(1e-9, maximum(abs, b_yz))

    blims = (-bmax, bmax) .* 0.8

    blevels = vcat([-bmax], range(blims[1], blims[2], length=31), [bmax])

    u_yz_plot = contourf(yu * 1e-3, zu, u_yz',
                         xlabel = "y (km)",
                         ylabel = "z (m)",
                         aspectratio = :equal,
                         linewidth = 0,
                         levels = ulevels,
                         clims = ulims,
                         xlims = ylims,
                         ylims = zlims,
                         color = :balance)

    bz_yz_plot = heatmap(yw * 1e-3, zw, b_z_yz',
                         xlabel = "y (km)",
                         ylabel = "z (m)",
                         aspectratio = :equal,
                         xlims = ylims,
                         ylims = zlims,
                         color = :balance)

    contour!(u_yz_plot,
             yc * 1e-3, zc, b_yz',
             linewidth = 1,
             color = :black,
             levels = blevels,
             legend = :none)
end

mp4(anim, filename*".mp4", fps = 8) # hide


text:
Random.seed!(515)
n, p = 100, 5
X  = randn(n, p)
X1 = hcat(X, ones(100))
θ  = randn(p+1)
y  = X1 * θ .+ 0.1 .* randn(n)

@testset "GPRegressor" begin
    gpr = GaussianProcessRegressor(random_state = 1)
    res, _, _ = fit(gpr, 1, X, y)
    @test res[1].score(X,y) ≈ 1.0
    @test norm(predict(gpr, res, X) .- y) / norm(y) ≤ 1e-10 # overfitting to the max
    fp = fitted_params(gpr, res)
    @test keys(fp) == (:X_train, :y_train, :kernel, :L, :alpha, :log_marginal_likelihood_value)
end

Xc2, yc2 = gen_classif(classes=["A", "B", "C"])

@testset "GPClassif" begin
    m, f = simple_test_classif_prob(GaussianProcessClassifier(), Xc2, yc2)
    fp = fitted_params(m, f)
    @test keys(fp) == (:kernel, :log_marginal_likelihood_value, :classes, :n_classes)
    infos = info_dict(m)
    @test infos[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test infos[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test !isempty(infos[:docstring])
end


text:
using Documenter, MiniOPF

makedocs(;
    modules=[MiniOPF],
    format=:html,
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/lsindoni/MiniOPF.jl/blob/{commit}{path}#L{line}",
    sitename="MiniOPF.jl",
    authors="Lorenzo Sindoni",
    assets=[],
)

deploydocs(;
    repo="github.com/lsindoni/MiniOPF.jl",
    target="build",
    julia="0.6",
    deps=nothing,
    make=nothing,
)


text:
using Latexify
using ParameterizedFunctions
using LaTeXStrings
using Test

f = @ode_def feedback begin
    dx = y/c_1 - x
    dy = x^c_2 - y
    end c_1 c_2

test_array = ["x/y * d" :x ; :( (t_sub_sub - x)^(2*p) ) 3//4 ]

@test latexify(f; starred=true) ==
raw"\begin{align*}
\frac{dx}{dt} =& \frac{y}{c_{1}} - x \\
\frac{dy}{dt} =& x^{c_{2}} - y
\end{align*}
"

#### Test the setting of default keyword arguments.

@test latexify("x * y") ==
raw"$x \cdot y$"

set_default(cdot = false)

@test latexify("x * y") ==
raw"$x y$"

@test get_default() == Dict{Symbol,Any}(:cdot => false)

set_default(cdot = true, transpose = true)

@test get_default() == Dict{Symbol,Any}(:cdot => true,:transpose => true)
@test get_default(:cdot) == true
@test get_default(:cdot, :transpose) == (true, true)
@test get_default([:cdot, :transpose]) == Bool[1, 1]

reset_default()
@test get_default() == Dict{Symbol,Any}()

@test latexify("x * y") ==
raw"$x \cdot y$"


text:
using DataArrays

## Convenience functions
## unlike the charts, these take an IO object to plot

## Plot a function
## GoogleCharts does not like Inf values, but handles NaN gracefully. We replace
## args like [:title=>"My title"]
function plot(io::Union(IO, String, Nothing), f::Function, a::Real, b::Real, args::Dict)
    n = 250
    xs = linspace(a, float(b), n)
    ys = map(f, xs)
    ys[ ys.== Inf] = NaN

    d = DataFrame(x=xs,y=ys)

    chart = line_chart(d, merge(args, [:curveType => "function"]), nothing, nothing)
    ## Trade off (likely from ignorance)
    ## returning `chart` works with Interact
    ## returning the following allows multiple plot call per cell
    #    io == nothing ? redisplay(chart)  : render(io, chart)
    chart
end

plot(f::Function, a::Real, b::Real, args::Dict) = plot(nothing, f, a, b, args)
function plot(f::Function, a::Real, b::Real; kwargs...)
    d = Dict()
    [d[s] = v for (s, v) in kwargs]
    plot(nothing, f, a, b, d)
end

## 1 or more functions at once
function  plot(io::Union(IO, String, Nothing), fs::Vector{Function}, a::Real, b::Real, args::Dict)
    n = 250
    xs = linspace(a, float(b), n)
    d = DataFrame()
    d[:x] =  xs

    for i in 1:length(fs)
        ys = float(map(fs[i], xs))
        ys[ ys.== Inf] = NaN
        d[symbol("f$i")] = ys
    end

    line_chart(d, merge(args, [:curveType => "function"]), nothing, nothing)
end

plot(fs::Vector{Function}, a::Real, b::Real, args::Dict) = plot(nothing, fs, a, b, args)
function plot(fs::Vector{Function}, a::Real, b::Real; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    plot(nothing, fs, a, b, d)
end

## parametric plots -- tuples of functions
function plot(io::Union(IO, String, Nothing), fs::Tuple, a::Real, b::Real, args::Dict)
    u = linspace(a, b, 250)
    x = map(fs[1], u)
    y = map(fs[2], u)
    args[:curveType] = "function"
    plot(io, x, y, args)
end
plot( fs::Tuple, a::Real, b::Real, args::Dict)=plot(nothing, fs, a, b, args)
function plot(io::Union(IO, String, Nothing), fs::Tuple, a::Real, b::Real; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    plot(io, fs, a, b, d)
end
plot( fs::Tuple, a::Real, b::Real; kwargs...)=plot(nothing, fs, a, b;  kwargs...)

## plot x,y
function plot{S <: Real, T <: Real}(io::Union(IO, String, Nothing),
                                    x::Union(DataArray{S, 1}, Range1{S}, Vector{S}),
                                    y::Union(DataArray{T, 1}, Range1{T}, Vector{T}), args::Dict)
    if !(length(x)  == length(y)) error("Lengths don't match") end
    d = DataFrame(x=x, y=y)

    line_chart(d, args, nothing, nothing)
end

function plot{S <: Real, T <: Real}(io::Union(IO, String, Nothing),
                                    x::Union(DataArray{S, 1}, Range1{S}, Vector{S}),
                                    y::Union(DataArray{T, 1}, Range1{T}, Vector{T});
                                    kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    plot(io, x, y, d)
end

VectorLike = Union(DataArray, Range1, Vector)
plot(x::VectorLike, y::VectorLike, args::Dict) = plot(nothing, x, y, args)
function plot(x::VectorLike, y::VectorLike; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    plot(nothing, x, y, d)
end

## Plots with data frames
SymOrExpr = Union(Symbol, Expr)
plot(x::SymOrExpr, y::SymOrExpr, data::DataFrame, args::Dict) = plot(data[string(x)], data[string(y)], args)
function plot(x::SymOrExpr, y::SymOrExpr, data::DataFrame; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    plot(x, y, data, d)
end

## Scatterplots
function scatter(io::Union(IO, String, Nothing),
                 x::VectorLike, y::VectorLike,
                 args::Dict)
    d = DataFrame(x=x, y=y)
    scatter_chart(d, args)
end
scatter(x::VectorLike, y::VectorLike, args::Dict) = scatter(nothing, x, y, args)

function scatter(io::Union(IO, String, Nothing),
                 x::VectorLike, y::VectorLike; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    scatter(x, y, d)
end
scatter(x::VectorLike, y::VectorLike; kwargs...) = scatter(nothing, x, y; kwargs...)

function scatter(io::Union(IO, String, Nothing), x::SymOrExpr, y::SymOrExpr, data::DataFrame, args::Dict)
    scatter(io, data[x], data[y], args)
end

scatter(x::SymOrExpr, y::SymOrExpr, data::DataFrame, args::Dict) = scatter(nothing, data[x], data[y], args)

function scatter(io::Union(IO, String, Nothing),
                 x::SymOrExpr, y::SymOrExpr, data::DataFrame; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    !haskey(d, :hAxis) && (d[:hAxis] = [:title=>string(x)])
    !haskey(d, :vAxis) && (d[:vAxis] = [:title=>string(y)])

    scatter(io,data[x], data[y], d)
end
scatter(x::SymOrExpr, y::SymOrExpr, data::DataFrame; kwargs...) = scatter(nothing, x, y, data; kwargs...)

function NaNWrap(idx::Integer, gp::GroupedDataFrame)
    [[i == idx ? gp[i][:,2] : rep(NaN, size(gp[i])[1]) for i in 1:length(gp)]...]
end
## Assume dataframes are in [x,y, group,...] format
## using RDatasets
## iris = data("datasets", "iris")
## d=iris[:, [2,3,6]]
## gp = groupby(d, "i")
## scatter(gp)
function scatter(io::Union(IO, String, Nothing), gp::GroupedDataFrame, args::Dict)
    n = length(gp)
    d = DataFrame(x = [[gp[i][:,1] for i in 1:n]...])
    for i in 1:n
        nm = symbol("x$i")
        d[nm] = GoogleCharts.NaNWrap(i, gp)
    end

    scatter_chart(d, merge([:hAxis=>[:title=>names(gp[1])[1]], :vAxis=>[:title=>names(gp[1])[2]]], args))
end
scatter(gp::GroupedDataFrame, args::Dict) = scatter(nothing, gp, args)
function scatter(io::Union(IO, String, Nothing), gp::GroupedDataFrame; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    scatter(io, gp, d)
end
scatter(gp::GroupedDataFrame; kwargs...)  = scatter(nothing, gp; kwargs...)

## Boxplots are not right! I don't know how to add points in a combo chart!
function boxplot_stats{T<:Number}(x::Vector{T}; coef=1.5)
    sort!(x)
    fivenum = quantile(x, 0:(.25):1)
    outliers = T[]

    IQR = fivenum[4] - fivenum[2]
    upper =  fivenum[4] + coef * IQR
    lower =  fivenum[2] - coef * IQR

    if fivenum[5] > upper
        fivenum[5] = upper
        outliers = [outliers, x[ x.> upper]]
    end

    if fivenum[1] < lower
        fivenum[1] = lower
        outliers = [outliers, x[ x.< lower]]
    end

    (fivenum, outliers)
end

## Boxplot only shows five number summary, no marking of outliers
## in fact, no marking of the median!
function boxplot(io::Union(IO, String, Nothing), x::Vector, args::Dict)
    # stats = boxplot_stats(x)
    sort!(x)
    stats = quantile(x, 0:1/4:1)
    data = DataFrame(names=["x"],
                     IQR = stats[1], # really min?
                     q3 = stats[4],
                     q1 = stats[2],
                     max= stats[5]
                     )

    candlestick_chart(data, merge(args, [:legend=>nothing]))
end
boxplot(x::Vector, args::Dict) = boxplot(nothing, x, args)

function boxplot(io::Union(IO, String, Nothing), x::Vector; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    boxplot(io, x, d)
end
boxplot(x::Vector; kwargs...) = boxplot(nothing, x; kwargs...)

## Not great, as ordering in d is not guaranteed
function boxplot(io::Union(IO, String, Nothing), d::Dict, args::Dict)
    ## not efficient, but whatever
    nms = String[string(k) for k in keys(d)]
    vals = [quantile(v, u)::Real for (k,v) in d, u in 0:.25:1]

    data = DataFrame(names=nms,
                     IQR = vals[:,1],
                     q3 = vals[:,4],
                     q1 = vals[:,2],
                     max = vals[:,5])

    candlestick_chart(data, merge(args, [:legend=>nothing]))
end
boxplot(d::Dict, args::Dict) = boxplot(nothing, d, args)
function boxplot(io::Union(IO, String, Nothing), D::Dict; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    boxplot(io, D, d)
end
boxplot(d::Dict; kwargs...) = boxplot(nothing, d; kwargs...)

## XXX This is broken. How to easily get names from GroupedDataFrame
function boxplot(io::Union(IO, String, Nothing), gp::GroupedDataFrame, args::Dict)
    n = length(gp)
    nms = (gp |> sum)[:,1] ## hack!
    vals = [quantile(v[:,1], u) for v in gp, u in 0:.25:1]

    data = DataFrame(names=nms,
                     IQR = vals[:,1],
                     q3 = vals[:,4],
                     q1 = vals[:,2],
                     max = vals[:,5])

    candlestick_chart(data, merge(args, [:title=>"boxplot", :legend=>nothing]))
end
boxplot(gp::GroupedDataFrame, args::Dict) = boxplot(gp, args)
function boxplot(io::Union(IO, String, Nothing), gp::GroupedDataFrame; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    boxplot(io, gp, d)
end
boxplot(gp::GroupedDataFrame; kwargs...) = boxplot(nothing, gp; kwargs...)

### histogram
function histogram(io::Union(IO, String, Nothing), x::Vector, args::Dict; n::Integer=0)

    if n == 0
        n = iceil(log2(length(x)) + 1) #  sturges from R
    end
    bins, counts = hist(x, n)

    centers = (bins[1:end-1] .+ bins[2:end]) / 2
    data = DataFrame(x=centers, counts=counts)

    column_chart(data, merge(args, [:legend=>nothing,
                                            :hAxis=>[:maxValue=>maximum(bins), :minValue=>minimum(bins)], :bar=>[:groupWidth=>"99%"]]))

end
histogram(x::Vector, args::Dict; n::Integer=0) = histogram(nothing, io, args, n=n)

function histogram(io::Union(IO, String, Nothing), x::Vector; n::Integer=0, kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    histogram(io, x, d, n=n)
end
histogram(x::Vector; n::Integer=0, kwargs...) = histogram(nothing, x; n=n, kwargs...)

function histogram(io::Union(IO, String, Nothing), x::Vector; kwargs...)
    d = Dict()
    [d[s] = v for (s,v) in kwargs]
    histogram(io, x, d)
end
histogram(x::Vector; kwargs...)  = histogram(nothing, x; kwargs...)

##################################################

## Surface plot is different. XXX We need to be able to pass in arguments below
surface_tpl = """
var tooltipStrings = {{{:tooltips}}};
var {{:chart_id}}_data = {{{:datatable}}}
         var surfacePlot = new greg.ross.visualisation.SurfacePlot(document.getElementById("surfacePlotDiv_{{:chart_id}}"));

         // Don't fill polygons in IE. It's too slow.
         var fillPly = true;

         // Define a colour gradient.
         var colour1 = {red:0, green:0, blue:255};
         var colour2 = {red:0, green:255, blue:255};
         var colour3 = {red:0, green:255, blue:0};
         var colour4 = {red:255, green:255, blue:0};
         var colour5 = {red:255, green:0, blue:0};
         var colours = [colour1, colour2, colour3, colour4, colour5];

         // Axis labels.
         var xAxisHeader = "X";
         var yAxisHeader = "Y";
         var zAxisHeader = "Z";

         var options = {
           xPos: {{:xPos}},
           yPos: {{:yPos}},
           width: {{:width}},
           height: {{:height}},
           colourGradient: colours,
           fillPolygons: fillPly, tooltips: tooltipStrings,
           xTitle: xAxisHeader,yTitle: yAxisHeader, zTitle: zAxisHeader,
           restrictXRotation: false
         };

        surfacePlot.draw({{:chart_id}}_data, options);
"""

type SurfacePlot
    x::Dict
end

function render(io, p::SurfacePlot)
    plt = Mustache.render(surface_tpl, p.x)

    tpl = Mustache.template_from_file(Pkg.dir("GoogleCharts", "tpl", "surface.html"))
    f = tempname() * ".html"
    io = open(f, "w")
    Mustache.render(io, tpl, [:surfaceplot=>plt, :chart_id=>p.x[:chart_id]])
    close(io)
    open_url(f)
end

function Base.display(io::IO, p::SurfacePlot)
    if io === STDOUT
        render(nothing, p)
    else
        writemime(io, "text/html", p)
    end
end
function writemime(io::IO, ::MIME"text/html", p::SurfacePlot)
    plt = Mustache.render(surface_tpl, p.x)
    out = """

"""
    print(io, out)
end

## XXX This needs a way to pass options in ...
function surfaceplot(io::Union(IO, String, Nothing), f::Function, x::Vector, y::Vector;
                     xPos::Real=0,
                     yPos::Real=0,
                     width::Int=600,
                     height::Int=400

)
    chart_id = get_id()
    m = Float64[f(x,y) for x in x, y in y]
    d = convert(DataFrame, m)

    tool_tip(x,y) = "(" * (map(u -> round(u, 2), [x,y,f(x,0)]) |> u -> join(u, ", ")) * ")"
    tooltips = [tool_tip(x,y) for x in x, y in y] |> u -> reshape(u, length(x)*length(y)) |> JSON.json
    chart = SurfacePlot([:datatable => make_data_array(chart_id, d),
                         :tooltips=>tooltips,
                         :chart_id=>chart_id,
                         :xPos=>xPos, :yPos=>yPos,
                         :width=>width, :height=>height
                         ])
    chart
end
surfaceplot(f::Function, x::Vector, y::Vector; kwargs...)   = surfaceplot(nothing, f, x, y; kwargs...)


text:
"""
    SimpleWeightedDiGraph{T, U}

A type representing a directed graph with weights of type `U`.

Note that adding or removing vertices or edges is not particularly performant;
see MetaGraphs.jl for possible alternatives.
"""
mutable struct SimpleWeightedDiGraph{T<:Integer, U<:Real} <: AbstractSimpleWeightedGraph{T, U}
    weights::SparseMatrixCSC{U, T} # indexed by [dst, src]
    function SimpleWeightedDiGraph{T, U}(adjmx::SparseMatrixCSC{U,T}; permute=true) where T <: Integer where U <: Real
        dima,dimb = size(adjmx)
        isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
        permute ? new{T, U}(permutedims(adjmx)) : new{T, U}(adjmx)
    end

end

SimpleWeightedDiGraph{T}(adjmx::SparseMatrixCSC{U, T}; permute=true) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph{T, U}(adjmx; permute=permute)

SimpleWeightedDiGraph(adjmx::SparseMatrixCSC{U, T}; permute=true) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph{T, U}(adjmx; permute=permute)

SimpleWeightedDiGraph(m::AbstractMatrix{U}) where U <: Real =
    SimpleWeightedDiGraph{Int, U}(SparseMatrixCSC{U, Int}(m))
SimpleWeightedDiGraph{T}(m::AbstractMatrix{U}) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph{T, U}(SparseMatrixCSC{U, T}(m))
SimpleWeightedDiGraph{T, U}(m::AbstractMatrix) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph{T, U}(SparseMatrixCSC{U, T}(m))

SimpleWeightedDiGraph(g::SimpleWeightedDiGraph) = SimpleWeightedDiGraph(g.weights, false)
SimpleWeightedDiGraph{T,U}(g::SimpleWeightedDiGraph) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph(SparseMatrixCSC{U, T}(g.weights); permute=false)

ne(g::SimpleWeightedDiGraph) = nnz(g.weights)

function SimpleWeightedDiGraph{T,U}(n::Integer = 0) where T<:Integer where U<:Real
    weights = spzeros(U, T, T(n), T(n))
    return SimpleWeightedDiGraph{T, U}(weights)
end

# Graph()
SimpleWeightedDiGraph() = SimpleWeightedDiGraph{Int, Float64}()

# Graph(6), Graph(0x5)
SimpleWeightedDiGraph(n::T) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(n)

# Graph(UInt8)
SimpleWeightedDiGraph(::Type{T}) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(zero(T))

# Graph(UInt8, Float32)
SimpleWeightedDiGraph(::Type{T}, ::Type{U}) where T<:Integer where U<:Real = SimpleWeightedDiGraph{U, T}(zero(T))

# DiGraph(AbstractSimpleGraph)
function SimpleWeightedDiGraph(g::LightGraphs.SimpleGraphs.AbstractSimpleGraph, ::Type{U}=Float64) where U <: Real
    T = eltype(g)
    return SimpleWeightedDiGraph{T}(adjacency_matrix(g, U))
end

# DiGraph(AbstractSimpleGraph, defaultweight)
function SimpleWeightedDiGraph(g::LightGraphs.SimpleGraphs.AbstractSimpleGraph, x::U) where U <: Real
    T = eltype(g)
    return SimpleWeightedDiGraph{T, U}(x.*adjacency_matrix(g, U))
end

# DiGraph(srcs, dsts, weights)
function SimpleWeightedDiGraph(i::AbstractVector{T}, j::AbstractVector{T}, v::AbstractVector{U}; combine = +) where T<:Integer where U<:Real
    m = max(maximum(j), maximum(i))
    SimpleWeightedDiGraph{T, U}(sparse(j, i, v, m, m, combine); permute=false)
end

edgetype(::SimpleWeightedDiGraph{T, U}) where T<:Integer where U<:Real = SimpleWeightedGraphEdge{T,U}

edges(g::SimpleWeightedDiGraph) = (SimpleWeightedEdge(x[2], x[1], x[3]) for x in zip(findnz(g.weights)...))
weights(g::SimpleWeightedDiGraph) = g.weights'

inneighbors(g::SimpleWeightedDiGraph, v::Integer) = g.weights[v,:].nzind

# add_edge! will overwrite weights.
function add_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    T = eltype(g)
    U = weighttype(g)
    s_, d_, w = Tuple(e)

    if w == zero(U)
        @warn "Note: adding edges with a zero weight to this graph type has no effect." maxlog=1 _id=:swd_add_edge_zero
        return false
    end

    s = T(s_)
    d = T(d_)
    (s in vertices(g) && d in vertices(g)) || return false
    @inbounds g.weights[d, s] = w
    return true
end

function rem_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    has_edge(g, e) || return false
    U = weighttype(g)
    @inbounds g.weights[dst(e), src(e)] = zero(U)
    return true
end

copy(g::SimpleWeightedDiGraph) =  SimpleWeightedDiGraph(copy(g.weights'))

==(g::SimpleWeightedDiGraph, h::SimpleWeightedDiGraph) = g.weights == h.weights

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleWeightedDiGraph}) = true
is_directed(::Type{SimpleWeightedDiGraph{T,U}}) where T where U = true
is_directed(g::SimpleWeightedDiGraph) = true


text:
# """
#     GDALInitGCPs(int,
#                  GDAL_GCP *) -> void
# """
# function initgcps(arg1::Integer,arg2)
#     ccall((:GDALInitGCPs,libgdal),Void,(Cint,Ptr{GDAL_GCP}),arg1,arg2)
# end

# """
#     GDALDeinitGCPs(int,
#                    GDAL_GCP *) -> void
# """
# function deinitgcps(arg1::Integer,arg2)
#     ccall((:GDALDeinitGCPs,libgdal),Void,(Cint,Ptr{GDAL_GCP}),arg1,arg2)
# end

# """
#     GDALDuplicateGCPs(int,
#                       const GDAL_GCP *) -> GDAL_GCP *
# """
# function duplicategcps(arg1::Integer,arg2)
#     ccall((:GDALDuplicateGCPs,libgdal),Ptr{GDAL_GCP},(Cint,Ptr{GDAL_GCP}),arg1,arg2)
# end

# """
# Generate Geotransform from GCPs.

# Given a set of GCPs perform first order fit as a geotransform.

# Due to imprecision in the calculations the fit algorithm will often return
# non-zero rotational coefficients even if given perfectly non-rotated inputs.
# A special case has been implemented for corner corner coordinates given in
# TL, TR, BR, BL order. So when using this to get a geotransform from 4 corner
# coordinates, pass them in this order.

# ### Parameters
# * `nGCP`             the number of GCPs being passed in.
# * `pasGCPs`          the list of GCP structures.
# * `padfGeoTransform` the six double array in which the affine geotransformation
# will be returned.
# * `bApproxOK`        If `FALSE` the function will fail if the geotransform is
# not essentially an exact fit (within 0.25 pixel) for all GCPs.

# ### Returns
# `TRUE` on success or `FALSE` if there aren't enough points to prepare a
# geotransform, the pointers are ill-determined or if `bApproxOK` is `FALSE` and
# the fit is poor.
# """
# gcpstogeotransform(nGCP::Integer,
#                     pasGCPs::Ptr{GDAL.GDAL_GCP},
#                     padfGeoTransform::Ptr{Cdouble},
#                     bApproxOK::Integer) =
#     GDALGCPsToGeoTransform(nGCP, pasGCPs, padfGeoTransform, bApproxOK)::Cint

"""
Invert Geotransform.

This function will invert a standard 3x2 set of GeoTransform coefficients. This
converts the equation from being pixel to geo to being geo to pixel.

### Parameters
* `gt_in`       Input geotransform (six doubles - unaltered).
* `gt_out`      Output geotransform (six doubles - updated).

### Returns
`TRUE` on success or `FALSE` if the equation is uninvertable.
"""
function invgeotransform!(gt_in::Vector{Cdouble}, gt_out::Vector{Cdouble})
    result = Bool(GDAL.invgeotransform(pointer(gt_in), pointer(gt_out)))
    result || error("Geotransform coefficients is uninvertable")
    gt_out
end

invgeotransform(gt_in::Vector{Cdouble}) = invgeotransform!(gt_in, Array(Cdouble, 6))

"""
Apply GeoTransform to x/y coordinate.

Applies the following computation, converting a (pixel,line) coordinate into a
georeferenced `(geo_x,geo_y)` location.
```C
    *pdfGeoX =  padfGeoTransform[0] +
                dfPixel * padfGeoTransform[1] +
                dfLine * padfGeoTransform[2];

    *pdfGeoY =  padfGeoTransform[3] +
                dfPixel * padfGeoTransform[4] +
                dfLine * padfGeoTransform[5];
```
### Parameters
* `padfGeoTransform`  Six coefficient GeoTransform to apply.
* `dfPixel`           input pixel position.
* `dfLine`            input line position.
* `pdfGeoX`           location for `geo_x` (easting/longitude)
* `pdfGeoY`           location for `geo_y` (northing/latitude)

"""
function applygeotransform(geotransform::Vector{Cdouble},
                           pixel::Cdouble,
                           line::Cdouble)
    geo_xy = Array(Cdouble, 2)
    geo_x = pointer(geo_xy);
    geo_y = geo_x + sizeof(Cdouble)
    GDAL.applygeotransform(pointer(geotransform), pixel, line, geo_x, geo_y)
    geo_xy
end

"""
Compose two geotransforms.

The resulting geotransform is the equivelent to `padfGT1` and then `padfGT2`
being applied to a point.

### Parameters
* `gt1`     the first geotransform, six values.
* `gt2`     the second geotransform, six values.
* `gtout`   the output geotransform, six values, may safely be the same
array as `gt1` or `gt2`.

"""
function composegeotransform!(gt1::Vector{Cdouble},
                              gt2::Vector{Cdouble},
                              gtout::Vector{Cdouble})
    GDAL.composegeotransform(pointer(gt1), pointer(gt2), pointer(gtout))
    gtout
end

composegeotransform(gt1::Vector{Cdouble}, gt2::Vector{Cdouble}) =
    composegeotransform!(gt1, gt2, Array(Cdouble, 6))


text:
##########################################################################
# Assignment types
type AssignmentResults <: Results
    results::Dict{ASCIIString, Any}
end
AssignmentResults() =
    AssignmentResults(Dict{ASCIIString, Any}())
Base.getindex(p::AssignmentResults, k::ASCIIString) =
    getindex(p.results, k)
Base.setindex!(p::AssignmentResults, v, k::ASCIIString) =
    setindex!(p.results, v, k)

immutable Assignment
    cell_assignment::Vector{Int} # cell_assignment[k] is the BS serving MS k
    cell_assignment_inverse::Vector{IntSet} # cell_assignment_inverse[i] is the IntSet of MSs that BS i is serving

    cluster_assignment::Vector{IntSet} # cluster_assignment[k] is the IntSet of BSs coordinating to MS k
    cluster_assignment_inverse::Vector{IntSet} # cluster_assignment_inverse[i] is the IntSet of MSs that BS i is coordinating to
end

Assignment() = Assignment(Array(Int, 0), [IntSet()], [IntSet()], [IntSet()])

Assignment(cell_assignment::Vector{Int}, num_BSs::Int) =
    Assignment(cell_assignment, ones(length(cell_assignment), num_BSs))

Assignment(cell_assignment::Matrix) =
    Assignment([ (findin(cell_assignment[k,:], [true]))[1] for k = 1:size(cell_assignment, 1) ],
               ones(length(cell_assignment), size(cell_assignment, 2)))

Assignment(cell_assignment::Vector{Int}, cluster_assignment::Matrix) =
    Assignment(cell_assignment,
               [ IntSet(findin(cell_assignment, i)) for i = 1:size(cluster_assignment, 2) ],
               [ IntSet(findin(cluster_assignment[k,:], [true])) for k = 1:size(cluster_assignment, 1) ],
               [ IntSet(findin(cluster_assignment[:,i], [true])) for i = 1:size(cluster_assignment, 2) ])

Assignment(cell_assignment::Matrix, cluster_assignment::Matrix) =
    Assignment([ findin(cell_assignment[k,:], [true])[1] for k = 1:size(cell_assignment, 1) ],
               [ IntSet(findin(cell_assignment[:,i], [true])) for i = 1:size(cell_assignment, 2) ],
               [ IntSet(findin(cluster_assignment[k,:], [true])) for k = 1:size(cluster_assignment, 1) ],
               [ IntSet(findin(cluster_assignment[:,i], [true])) for i = 1:size(cluster_assignment, 2) ])

# Helper functions to get lists of BSs or MSs
active_BSs(assignment) = unique(assignment.cell_assignment)
num_served_MSs(BS_id, assignment) = length(assignment.cell_assignment_inverse[BS_id])

serving_BS_id(MS_id, assignment) =
    assignment.cell_assignment[MS_id]
served_MS_ids(BS_id, assignment) =
    assignment.cell_assignment_inverse[BS_id]
served_MS_ids_except_me(MS_id, BS_id, assignment) =
    setdiff(assignment.cell_assignment_inverse[BS_id], MS_id)

coordinated_BS_ids(MS_id, assignment) =
    assignment.cluster_assignment[MS_id]
coordinated_MS_ids(BS_id, assignment) =
    assignment.cluster_assignment_inverse[BS_id]

function require_equal_num_MSs_per_cell(assignment)
    BS_ind_max = maximum(assignment.cell_assignment)

    Kc_test = length(served_MS_ids(1, assignment))
    for BS_ind = 2:BS_ind_max
        (length(served_MS_ids(BS_ind, assignment)) == Kc_test) || Lumberjack.error("BSs must all serve equal number of MSs.")
    end
end


text:
module GradientsTestsuite

using ..KernelAbstractions
using ..KernelGradients
using ..Enzyme
using ..Test

include("matmul.jl")

function testsuite(backend, backend_str, backend_mod, AT, DAT)
    @testset "Matmul" begin
        matmul_testsuite(backend, AT)
    end
end

end