text:
"""
    CharKey(value::Char)

A key representing a printable character.

# Fields

value::Char - The character value of the key
"""
struct CharKey
    """The character value of the key"""
    value::Char
end

"""
    SpecialKey

A key that is not representing a printable character.
"""
@enum SpecialKey begin
    arrowDown
    arrowLeft
    arrowRight
    arrowUp
    backspace
    del
    endKey
    escape
    f1
    f2
    f3
    f4
    f5
    f6
    f7
    f8
    f9
    f10
    f11
    f12
    home
    ins
    pgDn
    pgUp
    returnKey
end

"""
    UnknownKey(value)

A key that is not yet supported.

# Fields

value - The sequence of bytes resulting from pressing the key
"""
struct UnknownKey
    """The sequence of bytes resulting from pressing the key"""
    value::Union{
        UInt8,
        Tuple{UInt8, UInt8},
        Tuple{UInt8, UInt8, UInt8},
        Tuple{UInt8, UInt8, UInt8, UInt8},
        Tuple{UInt8, UInt8, UInt8, UInt8, UInt8},
        Tuple{UInt8, UInt8, UInt8, UInt8, UInt8, UInt8},
        Tuple{UInt8, UInt8, UInt8, UInt8, UInt8, UInt8, UInt8}
    }
end

"""
Modifier keys
"""
@enum ModKey::UInt8 begin
    alt   = 0b001
    ctrl  = 0b010
    shift = 0b100
end

"""
    ModifiedKey

A key with modifiers. Conveniently constructed with `|`, e.g. `ctrl |
CharKey('a')` for Ctrl+A.  Note that not all combinations are useable. For
example, Ctrl+J and Ctrl+M are both identical to [RETURN].

Combining with `|` is nearly always more readable, e.g.
`alt | ctrl | CharKey('q')` instead of `ModKey(0x03, CharKey('q')`

# Fields

mod_bitmap::UInt8 - Modifier key bitmap; From LSB: Alt, Ctrl, Shift
key::Union{CharKey, SpecialKey} - The base key
"""
struct ModifiedKey
    """Modifier key bitmap; From LSB: Alt, Ctrl, Shift"""
    mod_bitmap::UInt8
    """The base key"""
    key::Union{CharKey, SpecialKey}
end

(|)(x::ModKey, y::ModKey) = UInt8(x) | UInt8(y)
(|)(x::ModKey, y::Union{CharKey, SpecialKey, ModifiedKey}) = UInt8(x) | y
(|)(x::UInt8, y::Union{CharKey, SpecialKey}) = ModifiedKey(x, y)
(|)(x::UInt8, y::ModifiedKey) = ModifiedKey(x | y.mod_bitmap, y.key)


text:
using jInvVis
using PyPlot
using LinearAlgebra
import jInv.Utils.clear!

include("setupInversionExample.jl");

close("all")

function startExample()
	tstart = time_ns();

	function plotInputData()
		figure(20)
		imshow(m');colorbar();

		figure(21);
		imshow(reshape(gamma,tuple(size(m)...))');
	end
	plotInputData();

	Dobs, Wd = solveForwardProblemExtendedSources(m, pForp, omega, nfreq);

	# plot observed data
	for i=1:nfreq
		figure(i);
		imshow(Dobs[i][:,:]);
	end

	function plotModelResult(model)
		close(888);
		figure(888);
		imshow(model);colorbar();
		pause(1.0);
	end
	nx = size(m)[1];
	nz = size(m)[2];

	figure(22);
	imshow(mref'); colorbar();

	mc, Dc, pInv, Iact, mback, pMis = solveInverseProblemExtendedSources(pForp, Dobs, Wd, nfreq, nx, nz, mref,
								Mr, 0.5, 0.035,"ES_FWI.dat", true, plotModelResult);
	# Show results
	fullMc = reshape(Iact*pInv.modelfun(mc)[1] + mback,tuple((pInv.MInv.n)...));
	println("Model error:");
	println(norm(fullMc.-m));
	tend = time_ns();
	println("Runtime:");
	println((tend - tstart)/1.0e9);

	#Plot residuals
	figure(23);
	imshow((abs.(m.-fullMc))'); colorbar();
end

startExample()

ext: jl


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

name = "ERFA"
version = v"1.7.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/liberfa/erfa/releases/download/v$(version)/erfa-$(version).tar.gz", "80311221c96b45151f9a68517acc2b54b548b49d740e53cf409a6bc7502bd86b")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd erfa-*/
if [[ ${target} == i686-w64* ]] || [[ ${target} == x86_64-w64* ]]; then
    sed -i 's/LT_INIT/LT_INIT([win32-dll])/' configure.ac;
    sed -i 's/liberfa_la_LDFLAGS = -version-info \$(VI_ALL)/liberfa_la_LDFLAGS = -no-undefined -version-info $(VI_ALL)/' src/Makefile.am;
fi
autoreconf -fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using Knet:
    Knet,
    AutoGrad,
    dir,
    Data,
    minibatch,
    Param,
    @diff,
    value,
    params,
    grad,
    progress,
    progress!,
    KnetArray
using IterTools
using Base.Iterators
using Knet: Ops20.zeroone
using Printf, Random, Test, Statistics
using Plots;
default(fmt = :png);
include("utils.jl")
include("./deeper-cut/variant.data.jl")

# Define Linear Model
struct Linear
    w::Any
    b::Any
end
Linear(i::Int, o::Int) = Linear(param(o, i), param0(o))
(m::Linear)(x) = m.w * x .+ m.b
(m::Linear)(x, y) = Knet.nll(m(x), y)
(m::Linear)(data::Knet.Data) = mean(m(x, y) for (x, y) in data)

# MLP Layer - f: can be identity, relu etc,  - pdrop: dropout rate
struct MLPLayer
    w::Any
    b::Any
    f::Any
    pdrop::Any
end
MLPLayer(i::Int, o::Int, f = relu; pdrop = 0) = MLPLayer(param(o, i), param0(o), f, pdrop)
(l::MLPLayer)(x) = l.f.(l.w * dropout(x, l.pdrop) .+ l.b)

# Define convolutional layer:
struct Conv2
    w::Any
    b::Any
    f::Any
    is_pool_enabled::Bool
end
Conv2(w1, w2, nx, ny, f = relu; is_pool_enabled = true) =
    Conv2(param(w1, w2, nx, ny), param0(1, 1, ny, 1), f, is_pool_enabled)
function (c::Conv2)(x)
    if c.is_pool_enabled
        pool(c.f.(conv4(c.w, x) .+ c.b))
    else
        c.f.(conv4(c.w, x) .+ c.b)
    end
end

# Deconv
struct Deconv
    w::Any
    b::Any
    stride::Any
    padding::Any
    tag::String
end
Deconv(w1, w2, nx, ny; stride = 1, padding = 0, atype = Knet.atype(), tag = "") =
    Deconv(param(w1, w2, nx, ny; atype = atype), param0(1, 1, nx, 1), stride, padding, tag)

Deconv(w, b; stride = 1, padding = 0, atype = Knet.atype(), tag = "") =
    Deconv(param(w; atype = atype), param(b; atype = atype), stride, padding, tag)

function (dc::Deconv)(x)
    dc_res = deconv4(dc.w, x; stride = dc.stride, padding = dc.padding)

 #=   println("deconv - dc res")
    println(sum(isnan.(Array(dc_res))))
    if ( (sum(isnan.(Array(dc_res))) > 0 ))
          global  buggy_w = dc.w
    global buggy_x = x
    global buggy_stride = dc.stride
    global buggy_padding = dc.padding
    end=#

    res = dc_res .+ dc.b
    return res
end

function get_weights_sum(l::Deconv)
    return sum(abs2, l.w)
end

# Define dense layer:
struct Dense
    w::Any
    b::Any
    f::Any
    p::Any
end
Dense(i::Int, o::Int, f = relu; pdrop = 0) = Dense(param(o, i), param0(o), f, pdrop)
Dense(w, b; f = identity, pdrop = 0) =
    Dense(param(w; atype = Knet.atype()), param(b; atype = Knet.atype()), f, pdrop)
(d::Dense)(x) = d.f.(d.w * mat(dropout(x, d.p)) .+ d.b)

# Define a chain of layers and a loss function:
struct DeeperCutOption
    connect_res3_to_res5::Bool

    function DeeperCutOption(; connect_res3_to_res5 = false)
        return new(connect_res3_to_res5)
    end
end

# Deeper Cut Head
struct DeeperCutHead
    part_detection_head::Deconv
    loc_ref_head::Deconv
    tag::String
    is_loc_ref_enabled::Bool

    function DeeperCutHead(;
        part_detection_head = Deconv(
            4,
            4,
            global_num_joints,
            2048;
            padding = 1,
            stride = 2,
            tag = "part_detect_deconv",
        ),
        loc_ref_head = Deconv(
            4,
            4,
            global_num_joints * 4,
            2048;
            padding = 1,
            stride = 2,
            tag = "loc_ref_deconv",
        ),
        tag = "deeper_cut_head",
        is_loc_ref_enabled = false,
    )
        return new(part_detection_head, loc_ref_head, tag, is_loc_ref_enabled)
    end

    function DeeperCutHead(
        p_d_w,
        p_d_b,
        l_r_w,
        l_r_b;
        tag = "deeper_cut_head",
        is_loc_ref_enabled = false,
    )
        part_detection_head =
            Deconv(p_d_w, p_d_b; padding = 1, stride = 2, tag = "part_detect_deconv")
        loc_ref_head = Deconv(l_r_w, l_r_b; padding = 1, stride = 2, tag = "loc_ref_deconv")
        return new(part_detection_head, loc_ref_head, tag, is_loc_ref_enabled)
    end

end

function (deeper_cut_head::DeeperCutHead)(x)
    # TODO: how shall we check if the x is intact between heads
    part_detection_result = deeper_cut_head.part_detection_head(x)

    if deeper_cut_head.is_loc_ref_enabled == true
        loc_ref_result = deeper_cut_head.loc_ref_head(x)

        channel_dim = (size(part_detection_result) |> length) - 1
        combined_result = cat(part_detection_result, loc_ref_result; dims = channel_dim)
        return combined_result
    else
        return part_detection_result
    end
end

function get_weights_sum(l::DeeperCutHead)
    # sum(abs2, l.w)
    sum1 =  get_weights_sum(l.part_detection_head)
    sum2 =  get_weights_sum(l.loc_ref_head)
    return sum1 + sum2
end

# Chain
struct Chain
    layers::Any
    lambda1::Any
    lambda2::Any
    loss::Any
    deeperCutOption::Any
    Chain(layers...; lambda1 = 0, lambda2 = 0, loss = nll, deeperCutOption = nothing) =
        new(layers, lambda1, lambda2, loss, deeperCutOption)
end

# The prediction and average loss do not change
function (c::Chain)(x)
    connection_from3_to5 = nothing
    connection_from3_to5_loc_ref = nothing

    for l in c.layers
        x = l(x)

        if c.deeperCutOption != nothing && c.deeperCutOption.connect_res3_to_res5
            layer_tag = get_object_tag(l)
            if layer_tag == 3
                connection_from3_to5 = l.conv3_for_deepercut_output(x)
                connection_from3_to5_loc_ref = l.conv3_for_deepercut_output_loc_ref(x)
            end
        end

        if c.deeperCutOption != nothing && c.deeperCutOption.connect_res3_to_res5
            layer_tag = get_object_tag(l)
            if layer_tag == "deeper_cut_head"
                if l.is_loc_ref_enabled == true
                    channel_dim = (size(connection_from3_to5) |> length) - 1
                    combined_result = cat(
                        connection_from3_to5,
                        connection_from3_to5_loc_ref;
                        dims = channel_dim,
                    )
                    # TODO: This introduces problems when the dimensions are odd
                    x = x .+ combined_result
                else
                    x = x .+ connection_from3_to5
                end
            end
        end

    end

    x
end

function (c::Chain)(x, y)
    loss = c.loss(c(x), y)
    if training() # Only apply regularization during training, only to weights, not biases.
        c.lambda1 != 0 && (loss += c.lambda1 * sum(sum(abs, l.w) for l in c.layers))
        c.lambda2 != 0 && (loss += c.lambda2 * sum( get_weights_sum(l) for l in c.layers))
    end
    return loss
end
(c::Chain)(d::Data) = mean(c(x, y) for (x, y) in d)
(c::Chain)(d::VariantData) = mean(c(x, y) for (x, y) in d)

# This became redundant now
struct ResLayerConv
    w::Any
    padding::Any
    stride::Any
end
# Random init
ResLayerConv(w1, w2, nx, ny; padding = 0, stride = 1) =
    ResLayerConv(param(w1, w2, nx, ny), padding, stride)
# Predeterminde weights
ResLayerConv(w; padding = 0, stride = 1) = ResLayerConv(param(w), padding, stride)
(rl0::ResLayerConv)(x) = conv4(rl0.w, x; padding = rl0.padding, stride = rl0.stride)

struct BatchNormLayer
    w::Any
    ms::Any

    function BatchNormLayer(pre_w, pre_ms; freeze = false)
        res_mean = popfirst!(pre_ms)
        # Trick to arrange variance value for new(er) batchnorm

        # TODO: DO NOT FORGET TO ADD IT

        res_variance = popfirst!(pre_ms) .^ 2 # .- 1e-5
        ms = bnmoments(momentum=0.9, mean = res_mean, var = res_variance)

        w1 = pre_w[1]
        w2 = pre_w[2]
        w1 = vec(w1)
        w2 = vec(w2)
        w = vcat(w1, w2)
        if freeze
            return new(w, ms)
        else
            param_w = param(w, atype = Knet.atype())
            return new(param_w, ms)
        end
    end

end

function (batch_norm_layer::BatchNormLayer)(x)
    return batchnorm(x, batch_norm_layer.ms, batch_norm_layer.w; eps = 0, training=false)
end

# ResNet 50 initial layer
struct ResLayerX1_50
    batch_layer::Any
    conv_w::Any
    conv_b::Any
    padding::Any
    stride::Any
    pool_window_size::Any
    pool_stride::Any
    pool_padding::Any
end

function ResLayerX1_50(
    w,
    ms;
    padding = 3,
    stride = 2,
    pool_window_size = 3,
    pool_stride = 2,
    pool_padding = 1,
)
    bnl = BatchNormLayer(w[3:4], ms)
    return ResLayerX1_50(
        bnl,
        param(w[1]; atype = Knet.atype()),
        param(w[2]; atype = Knet.atype()),
        padding,
        stride,
        pool_window_size,
        pool_stride,
        pool_padding,
    )
end

function (rlx1_50::ResLayerX1_50)(x)
    o =
        conv4(rlx1_50.conv_w, x; padding = rlx1_50.padding, stride = rlx1_50.stride) .+
        rlx1_50.conv_b
    o = rlx1_50.batch_layer(o)
    o = relu.(o)
    o = pool(
        o;
        window = rlx1_50.pool_window_size,
        stride = rlx1_50.pool_stride,
        padding = rlx1_50.pool_padding,
    )
    return o
end

function get_weights_sum(l::ResLayerX1_50)
    # sum(abs2, l.w)
    # TODO: BATCHNORMLAYER FIX
     sum1 = 0 # sum(abs2, l.batch_layer.w)
    sum2 = sum(abs2, l.conv_w)
    return sum1 + sum2
end

# X0
struct ResLayerX0
    batch_layer::Any
    conv_w::Any
    padding::Any
    stride::Any
    dilation::Any
end
# Predetermined weights
# TODO: should we try to make bnl params??
function ResLayerX0(w, ms; padding = 0, stride = 1, dilation = 1)
    bnl = BatchNormLayer(w[2:3], ms)
    return ResLayerX0(bnl, param(w[1]; atype = Knet.atype()), padding, stride, dilation)
end

function (rlx0::ResLayerX0)(x)
    # batchnorm_as_function(rlx0.batch_w, conv4(rlx0.conv_w, x; padding=rlx0.padding, stride=rlx0.stride), rlx0.ms)
    o = conv4(
        rlx0.conv_w,
        x;
        padding = rlx0.padding,
        stride = rlx0.stride,
        dilation = rlx0.dilation,
    )
    o = rlx0.batch_layer(o)
    return o
end

function get_weights_sum(l::ResLayerX0)
    # sum(abs2, l.w)
     # TODO: BATCHNORMLAYER FIX
    sum1 = 0 #sum(abs2, l.batch_layer.w)
    sum2 = sum(abs2, l.conv_w)
    return sum1 + sum2
end

# X1
struct ResLayerX1
    x0_layer::Any
    is_initial::Bool
end
ResLayerX1(w, ms; padding = 0, stride = 1, is_initial::Bool = false, dilation = 1) =
    ResLayerX1(
        ResLayerX0(w, ms; padding = padding, stride = stride, dilation = dilation),
        is_initial,
    )
function (rlx1::ResLayerX1)(x)
    relu_res = relu.(rlx1.x0_layer(x))
    if rlx1.is_initial
        return pool(relu_res; window = 3, stride = 2, padding = 1)
    else
        return relu_res
    end
end

function get_weights_sum(l::ResLayerX1)
    # sum(abs2, l.w)
    return get_weights_sum(l.x0_layer)
end

# X2
# TODO: can be constructed like Chain
struct ResLayerX2
    x1_a_layer::Any
    x1_b_layer::Any
    x0_c_layer::Any
end
ResLayerX2(w, ms; pads = [0, 1, 0], strides = [1, 1, 1], dilations = [1, 1, 1]) =
    ResLayerX2(
        ResLayerX1(
            w[1:3],
            ms;
            padding = pads[1],
            stride = strides[1],
            dilation = dilations[1],
        ),
        ResLayerX1(
            w[4:6],
            ms;
            padding = pads[2],
            stride = strides[2],
            dilation = dilations[2],
        ),
        ResLayerX0(
            w[7:9],
            ms;
            padding = pads[3],
            stride = strides[3],
            dilation = dilations[3],
        ),
    )
(rlx2::ResLayerX2)(x) = rlx2.x0_c_layer(rlx2.x1_b_layer((rlx2.x1_a_layer(x))))

function get_weights_sum(l::ResLayerX2)
    # sum(abs2, l.w)
    sum1 = get_weights_sum(l.x1_a_layer)
    sum2 = get_weights_sum(l.x1_b_layer)
    sum3 = get_weights_sum(l.x0_c_layer)
    return sum1 + sum2 + sum3
end

# X3
struct ResLayerX3
    x0_a_layer::ResLayerX0
    x2_b_layer::ResLayerX2
end
ResLayerX3(
    w,
    ms;
    pads = [0, 0, 1, 0],
    strides = [2, 2, 1, 1],
    b_layer_dilations = [1, 1, 1],
) = ResLayerX3(
    ResLayerX0(w[1:3], ms; padding = pads[1], stride = strides[1]),
    ResLayerX2(
        w[4:12],
        ms;
        pads = pads[2:4],
        strides = strides[2:4],
        dilations = b_layer_dilations,
    ),
)
function (rlx3::ResLayerX3)(x)
    res_a = rlx3.x0_a_layer(x)
    res_b = rlx3.x2_b_layer(x)
    return relu.(res_a .+ res_b)
end

function get_weights_sum(l::ResLayerX3)
    # sum(abs2, l.w)
    sum1 = get_weights_sum(l.x0_a_layer)
    sum2 = get_weights_sum(l.x2_b_layer)
    return sum1 + sum2
end

# X4
struct ResLayerX4
    x2_layer::Any
end
ResLayerX4(w, ms; pads = [0, 1, 0], strides = [1, 1, 1], dilations = [1, 1, 1]) =
    ResLayerX4(ResLayerX2(w, ms; pads = pads, strides = strides, dilations = dilations))
(rlx4::ResLayerX4)(x) = relu.(x .+ rlx4.x2_layer(x))

function get_weights_sum(l::ResLayerX4)
    # sum(abs2, l.w)
    sum1 = get_weights_sum(l.x2_layer)
    return sum1
end

# X5
struct ResLayerX5
    x3_layer::ResLayerX3
    x4_layers::Any
    is_next_fc::Bool
    conv3_for_deepercut_output::Any
    conv3_for_deepercut_output_loc_ref::Any
    tag::Any
end

function ResLayerX5(
    w,
    ms;
    strides = [2, 2, 1, 1],
    is_next_fc::Bool = false,
    b_layer_dilations = [1, 1, 1],
    b_layer_pads = [0, 1, 0],
    is_conv3_for_deepercut = false,
    tag = nothing,
)
    x3_layer::ResLayerX3 = ResLayerX3(w[1:12], ms; strides = strides)
    x4_layers = []
    for k = 13:9:length(w)
        layer = ResLayerX4(w[k:k+8], ms; dilations = b_layer_dilations, pads = b_layer_pads)
        push!(x4_layers, layer)
    end

    conv3_for_deepercut_output = nothing
    conv3_for_deepercut_output_loc_ref = nothing
    if is_conv3_for_deepercut
        # 14 layer for part detection
        conv3_for_deepercut_output =
            Conv2(1, 1, 512, global_num_joints, identity; is_pool_enabled = false)
        conv3_for_deepercut_output_loc_ref =
            Conv2(1, 1, 512, global_num_joints * 4, identity; is_pool_enabled = false)
    end

    return ResLayerX5(
        x3_layer,
        x4_layers,
        is_next_fc,
        conv3_for_deepercut_output,
        conv3_for_deepercut_output_loc_ref,
        tag,
    )
end

function (rlx5::ResLayerX5)(x)
    x = rlx5.x3_layer(x)
    for l in rlx5.x4_layers
        x = l(x)
    end

    if rlx5.is_next_fc
        return pool(x; stride = 1, window = 7, mode = 2)
    else
        return x
    end
end

function get_weights_sum(l::ResLayerX5)
    # sum(abs2, l.w)
    sum = get_weights_sum(l.x3_layer)
      for layer in l.x4_layers
       sum += get_weights_sum(layer)
    end
    return sum
end


text:
#!/usr/bin/env julia

using JLD2
using FileIO
using ArgParse
using PyCall
using SymmetricTensors
mpl = pyimport("matplotlib")
mc = pyimport("matplotlib.colors")
mpl.rc("text", usetex=true)
mpl.use("Agg")
using PyPlot

function union_size(list1::Vector, list2::Vector)
  return length(list1 ∩ list2)
end

function features_selection_results(data::Dict, d::Int = 0)
  malf_size = data["malf_size"]
  check_size = malf_size

  labels = ["malf"]
  algorithms = ["MEV", "JKN", "JSBS", "JKFS"]
  data_parts_label = sort([ "bands_$(alg)_$l" for l=labels for alg=algorithms])
  plot_data = Dict()
  for label = data_parts_label
    plot_data[label] = zeros(Int, malf_size+1)
  end

  for m=1:length(data["data"]), label = data_parts_label
    list1 = data["data"]["$m"][label][end-check_size+1-d:end]
    list2 = data["data"]["$m"]["malf"]
    no_common = union_size(list1, list2)
    plot_data[label][no_common+1] += 1
  end
plot_data, data_parts_label
end

function theoretical(sizedata::Int, malfsize::Int, found::Int)
  binomial(malfsize, found)*binomial(sizedata-malfsize, malfsize-found)/binomial(sizedata,malfsize)
end

function los(n, ksize, p, r = 500_000)
  ret = []
  for s in 1:r
    b = []
    x = collect(1:n)
    for a in 1:p
      y = rand(x)
      push!(b, y)
      filter!(i -> i != y, x)
    end
    push!(ret, (count(b .<= ksize)))
  end
  [count(ret .== i) for i in 0:ksize]./r, collect(0:ksize)
end

function plotdata(plot_data, data_parts_label, ν, malf_size, var_number, repeating, δ::Int = 0)
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  cols = ["red", "green", "blue", "brown", "grey"]
  p = ["--d", "--o", "--^", "--<", "--v"]
  cla()
  i = 0
  for data_label = data_parts_label
    i += 1
    data_y = plot_data[data_label]
    data_x = collect(0:(length(data_y)-1))
    println(data_y)
    d = replace(data_label, "bands_"=>"")
    plot(data_x, data_y/repeating, p[i], label=replace(d, "_malf"=>""), color = cols[i], linewidth = 0.8, markersize = 3.)
  end
  data_x = 0:malf_size
  data_y = los(var_number, malf_size, malf_size+δ)[1]
  plot(data_x, data_y, "--x", label="rand choice", color = "black", linewidth = 1., markersize = 3.)
  if true
    ax.legend(fontsize = 4.5, loc = 2, ncol = 2)
  end
  subplots_adjust(left = 0.15, bottom = 0.16)
  xlabel("no. selected features", labelpad = 0.)
  ylabel("selection probability", labelpad = 0.)
  savefig("$(ν)_$(δ)jkfs.pdf")
end

function main(args)
  s = ArgParseSettings("description")
  @add_arg_table! s begin
    "file"
    help = "the file name"
    arg_type = String
    "--delta", "-d"
    help = "no. additiona features"
    default = 0
    arg_type = Int
  end
  parsed_args = parse_args(s)
  data = load(parsed_args["file"])
  d = parsed_args["delta"]
  plot_data, data_parts_label = features_selection_results(data, d)
  plotdata(plot_data, data_parts_label, data["ν"], data["malf_size"], data["variables_no"], data["test_number"], d)
end

main(ARGS)


text:
import OpenNI
using OpenNI

include("simenvi.jl")

type NIMJSimMapper
    mj_sim::MJSimEnv

    ni_context::NIContext
    ni_usergen::NIUserGenerator
    ni_userID::OpenNI.NIUserID

    pointMapper::Function

    function NIMJSimMapper(sim::MJSimEnv, c::NIContext, gen::NIUserGenerator, uid, mapper)
        new(sim, c, gen, uid, mapper)
    end
end

function MapNIJointsToMJSim(mapper::NIMJSimMapper, jointPos)
    mjPos = deepcopy(jointPos[:, 1])
    for i=1:length(jointPos[:, 1])
        mjPos[i, :] = mapper.pointMapper(mapper, jointPos[i, :])
    end
    return mjPos
end


text:
import JSON
using GitHub, GitHub.name, GitHub.GitHubString, GitHub.Branch
using Base.Test

# This file tests various GitHubType constructors. To test for proper Nullable
# handling, most fields have been removed from the JSON samples used below.
# Sample fields were selected in order to cover the full range of type behavior,
# e.g. if the GitHubType has a few Nullable{Dates.DateTime} fields, at least one
# of those fields should be present in the JSON sample.

function test_show(g::GitHub.GitHubType)
    tmpio = IOBuffer()
    show(tmpio, g)

    # basically trivial, but proves that things aren't completely broken
    @test repr(g) == takebuf_string(tmpio)

    tmpio = IOBuffer()
    showcompact(tmpio, g)

    @test "$(typeof(g))($(repr(name(g))))" == takebuf_string(tmpio)
end

#########
# Owner #
#########

owner_json = JSON.parse(
"""
{
  "id": 1,
  "email": null,
  "html_url": "https://github.com/octocat",
  "login": "octocat",
  "updated_at": "2008-01-14T04:33:35Z",
  "hireable": false
}
"""
)

owner_result = Owner(
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(GitHubString(owner_json["login"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{Int}(Int(owner_json["id"])),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(owner_json["html_url"])),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(owner_json["updated_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{Dates.DateTime}(),
    Nullable{Bool}(Bool(owner_json["hireable"])),
    Nullable{Bool}()
)

@test Owner(owner_json) == owner_result
@test name(Owner(owner_json["login"])) == name(owner_result)
@test setindex!(GitHub.github2json(owner_result), nothing, "email") == owner_json

test_show(owner_result)

########
# Repo #
########

repo_json = JSON.parse(
"""
{
  "id": 1296269,
  "owner": {
    "login": "octocat"
  },
  "parent": {
    "name": "test-parent"
  },
  "full_name": "octocat/Hello-World",
  "private": false,
  "url": "https://api.github.com/repos/octocat/Hello-World",
  "language": null,
  "pushed_at": "2011-01-26T19:06:43Z",
  "permissions": {
    "admin": false,
    "push": false,
    "pull": true
  }
}
"""
)

repo_result = Repo(
    Nullable{GitHubString}(),
    Nullable{GitHubString}(GitHubString(repo_json["full_name"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{Owner}(Owner(repo_json["owner"])),
    Nullable{Repo}(Repo(repo_json["parent"])),
    Nullable{Repo}(),
    Nullable{Int}(Int(repo_json["id"])),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(repo_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(repo_json["pushed_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{Dates.DateTime}(),
    Nullable{Bool}(),
    Nullable{Bool}(),
    Nullable{Bool}(),
    Nullable{Bool}(),
    Nullable{Bool}(Bool(repo_json["private"])),
    Nullable{Bool}(),
    Nullable{Dict}(repo_json["permissions"])
)

@test Repo(repo_json) == repo_result
@test name(Repo(repo_json["full_name"])) == name(repo_result)
@test setindex!(GitHub.github2json(repo_result), nothing, "language") == repo_json

test_show(repo_result)

##########
# Commit #
##########

commit_json = JSON.parse(
"""
{
  "url": "https://api.github.com/repos/octocat/Hello-World/commits/6dcb09b5b57875f334f61aebed695e2e4193db5e",
  "sha": "6dcb09b5b57875f334f61aebed695e2e4193db5e",
  "html_url": null,
  "commit": {
    "message": "Fix all the bugs",
    "comment_count": 0
  },
  "author": {
    "login": "octocat"
  },
  "parents": [
    {
      "sha": "6dcb09b5b57875f334f61aebed695e2e4193db5e"
    },
    {
      "sha": "7ed9340c309dd91757664cee6d857d161c14e095"
    }
  ],
  "stats": {
    "total": 108
  },
  "files": [
    {
      "filename": "file1.txt"
    },
    {
      "filename": "file2.txt"
    }
  ]
}
"""
)

commit_result = Commit(
    Nullable{GitHubString}(GitHubString(commit_json["sha"])),
    Nullable{GitHubString}(),
    Nullable{Owner}(Owner(commit_json["author"])),
    Nullable{Owner}(),
    Nullable{Commit}(Commit(commit_json["commit"])),
    Nullable{HttpCommon.URI}(HttpCommon.URI(commit_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{Vector{Commit}}(map(Commit, commit_json["parents"])),
    Nullable{Dict}(commit_json["stats"]),
    Nullable{Vector{Content}}(map(Content, commit_json["files"])),
    Nullable{Int}()
)

@test Commit(commit_json) == commit_result
@test name(Commit(commit_json["sha"])) == name(commit_result)
@test setindex!(GitHub.github2json(commit_result), nothing, "html_url") == commit_json

test_show(commit_result)

##########
# Branch #
##########

branch_json = JSON.parse(
"""
{
  "name": "master",
  "sha": null,
  "protection": {
    "enabled": false,
    "required_status_checks": {
      "enforcement_level": "off",
      "contexts": []
    }
  },
  "commit": {
    "sha": "7fd1a60b01f91b314f59955a4e4d4e80d8edf11d"
  },
  "user": {
    "login": "octocat"
  },
  "repo": {
    "full_name": "octocat/Hello-World"
  }
}
"""
)

branch_result = Branch(
    Nullable{GitHubString}(GitHubString(branch_json["name"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{Commit}(Commit(branch_json["commit"])),
    Nullable{Owner}(Owner(branch_json["user"])),
    Nullable{Repo}(Repo(branch_json["repo"])),
    Nullable{Dict}(),
    Nullable{Dict}(branch_json["protection"])
)

@test Branch(branch_json) == branch_result
@test name(Branch(branch_json["name"])) == name(branch_result)
@test setindex!(GitHub.github2json(branch_result), nothing, "sha") == branch_json

test_show(branch_result)

###########
# Comment #
###########

comment_json = JSON.parse(
"""
{
  "url": "https://api.github.com/repos/octocat/Hello-World/comments/1",
  "id": 1,
  "position": null,
  "body": "Great stuff",
  "user": {
    "login": "octocat"
  },
  "created_at": "2011-04-14T16:00:49Z"
}
"""
)

comment_result = Comment(
    Nullable{GitHubString}(GitHubString(comment_json["body"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{Int}(Int(comment_json["id"])),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(comment_json["created_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(comment_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{Owner}(Owner(comment_json["user"]))
)

@test Comment(comment_json) == comment_result
@test name(Comment(comment_json["id"])) == name(comment_result)
@test setindex!(GitHub.github2json(comment_result), nothing, "position") == comment_json

test_show(comment_result)

###########
# Content #
###########

content_json = JSON.parse(
"""
{
  "type": "file",
  "path": "lib/octokit.rb",
  "size": 625,
  "encoding": null,
  "url": "https://api.github.com/repos/octokit/octokit.rb/contents/lib/octokit.rb"
}
"""
)

content_result = Content(
    Nullable{GitHubString}(GitHubString(content_json["type"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(GitHubString(content_json["path"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(content_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{Int}(content_json["size"])
)

@test Content(content_json) == content_result
@test name(Content(content_json["path"])) == name(content_result)
@test setindex!(GitHub.github2json(content_result), nothing, "encoding") == content_json

test_show(content_result)

##########
# Status #
##########

status_json = JSON.parse(
"""
{
  "created_at": "2012-07-20T01:19:13Z",
  "description": "Build has completed successfully",
  "id": 1,
  "context": null,
  "url": "https://api.github.com/repos/octocat/Hello-World/statuses/1",
  "creator": {
    "login": "octocat"
  },
  "statuses": [
    {
      "id": 366962428
    }
  ],
  "repository": {
    "full_name": "JuliaWeb/GitHub.jl"
  }
}
"""
)

status_result = Status(
    Nullable{Int}(Int(status_json["id"])),
    Nullable{Int}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(GitHubString(status_json["description"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(status_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(status_json["created_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{Owner}(Owner(status_json["creator"])),
    Nullable{Repo}(Repo(status_json["repository"])),
    Nullable{Vector{Status}}(map(Status, status_json["statuses"]))
)

@test Status(status_json) == status_result
@test name(Status(status_json["id"])) == name(status_result)
@test setindex!(GitHub.github2json(status_result), nothing, "context") == status_json

test_show(status_result)

###############
# PullRequest #
###############

pr_json = JSON.parse(
"""
{
  "url": "https://api.github.com/repos/octocat/Hello-World/pulls/1347",
  "number": 1347,
  "body": "Please pull these awesome changes",
  "assignee": {
    "login": "octocat"
  },
  "milestone": {
    "id": 1002604,
    "number": 1,
    "state": "open",
    "title": "v1.0"
  },
  "locked": false,
  "created_at": "2011-01-26T19:01:12Z",
  "head": {
    "ref": "new-topic"
  }
}
"""
)

pr_result = PullRequest(
    Nullable{Branch}(),
    Nullable{Branch}(Branch(pr_json["head"])),
    Nullable{Int}(Int(pr_json["number"])),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{Int}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(GitHubString(pr_json["body"])),
    Nullable{GitHubString}(),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(pr_json["created_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{Dates.DateTime}(),
    Nullable{Dates.DateTime}(),
    Nullable{HttpCommon.URI}(HttpCommon.URI(pr_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{Owner}(Owner(pr_json["assignee"])),
    Nullable{Owner}(),
    Nullable{Owner}(),
    Nullable{Dict}(pr_json["milestone"]),
    Nullable{Dict}(),
    Nullable{Bool}(),
    Nullable{Bool}(),
    Nullable{Bool}(pr_json["locked"])
)

@test PullRequest(pr_json) == pr_result
@test name(PullRequest(pr_json["number"])) == name(pr_result)
@test GitHub.github2json(pr_result) == pr_json

test_show(pr_result)

#########
# Issue #
#########

issue_json = JSON.parse(
"""
{
  "url": "https://api.github.com/repos/octocat/Hello-World/issues/1347",
  "number": 1347,
  "title": "Found a bug",
  "user": {
    "login": "octocat"
  },
  "labels": [
    {
      "url": "https://api.github.com/repos/octocat/Hello-World/labels/bug",
      "name": "bug",
      "color": "f29513"
    }
  ],
  "pull_request": {
    "url": "https://api.github.com/repos/octocat/Hello-World/pulls/1347",
    "html_url": "https://github.com/octocat/Hello-World/pull/1347"
  },
  "locked": false,
  "closed_at": null,
  "created_at": "2011-04-22T13:33:48Z"
}
"""
)

issue_result = Issue(
    Nullable{Int}(),
    Nullable{Int}(Int(issue_json["number"])),
    Nullable{Int}(),
    Nullable{GitHubString}(GitHubString(issue_json["title"])),
    Nullable{GitHubString}(),
    Nullable{GitHubString}(),
    Nullable{Owner}(Owner(issue_json["user"])),
    Nullable{Owner}(),
    Nullable{Owner}(),
    Nullable{Dates.DateTime}(Dates.DateTime(chop(issue_json["created_at"]))),
    Nullable{Dates.DateTime}(),
    Nullable{Dates.DateTime}(),
    Nullable{Vector{Dict}}(Vector{Dict}(issue_json["labels"])),
    Nullable{Dict}(),
    Nullable{PullRequest}(PullRequest(issue_json["pull_request"])),
    Nullable{HttpCommon.URI}(HttpCommon.URI(issue_json["url"])),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{HttpCommon.URI}(),
    Nullable{Bool}(Bool(issue_json["locked"]))
)

@test Issue(issue_json) == issue_result
@test name(Issue(issue_json["number"])) == name(issue_result)
@test setindex!(GitHub.github2json(issue_result), nothing, "closed_at") == issue_json

test_show(issue_result)


text:
################################
# Solving the Replicator Equator equation,
# appearing in Corrolary 2.2 of
#  https://www.biorxiv.org/content/10.1101/623330v1
################################

cd(@__DIR__)
using HighDimPDE
using Random
using Test
using Flux
using Revise
using PyPlot
using UnPack
plotting = true

d = 10
x = fill(0.,d)  # initial point

tspan = (0e0, 1e0)
μ(X,p,t) = 0e0 # advection coefficients
σ(X,p,t) = 1e-1 # diffusion coefficients

ss0 = 1 / 4 / π #std g0

σ_sampling = 1.0
alg = MLP(M=5, K=10, L = 5, mc_sample = NormalSampling(σ_sampling) )

##########################
###### PDE Problem #######
##########################
g(x) = (2*π)^(-d/2) * ss0^(- Float64(d) * 5e-1) * exp.(-5e-1 *sum(x.^2)/ ss0) # initial condition
m(x) = - 5e-1 * sum(x.^2)
f(y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) = max(0.0, v_y) * (m(y) - max(0.0, v_z) * m(z) * (2.0 * π)^(d/2) * σ_sampling^d * exp(0.5 * sum(z.^2) / σ_sampling^2)) # nonlocal nonlinear part of the

# defining the problem
prob = PIDEProblem(g, f, μ, σ, tspan, x = x)

# solving
@time xs,ts,sol = solve(prob, alg, verbose = false, multithreading=true)
println("g(x) = " , g(x), "\nResult = ", sol[end])


text:
using Test

import Random

include("robot-name.jl")

# Random names means a risk of collisions.
const history = Set{String}()

isname(x) = occursin(r"^[A-Z]{2}[0-9]{3}$", x)

@testset "one robot" begin
    global r1 = Robot()
    push!(history, name(r1))

    @testset "name of robot is valid" begin
        @test isname(name(r1))
    end

    @testset "names of robot instance are valid and unique in history" begin
        # @testset sets the same seed every time for reproducibility, but we're
        # undoing that deliberately with Random.seed!() to increase the chance
        # of seeing collisions.
        Random.seed!()
        for i in 1:100
            reset!(r1)
            @test isname(name(r1))
            @test !in(name(r1), history)
            push!(history, name(r1))
        end
    end
end

@testset "two robots" begin
    global r2 = Robot()
    global r3 = Robot()

    @testset "names of robots are valid" begin
        @test isname(name(r2))
        @test isname(name(r3))
    end

    @testset "names of robots are not equal" begin
        @test name(r2) != name(r3)
    end

    @testset "names of both robots are unique in history" begin
        @test !in(name(r2), history)
        @test !in(name(r3), history)
    end

    push!(history, name(r2))
    push!(history, name(r3))
end

@testset "many robots" begin
    Random.seed!()
    robots = Robot[]

    for i in 1:100
        push!(robots, Robot())

        @testset "name of robot is valid and unique in history" begin
            @test isname(name(robots[i]))
            @test !in(name(robots[i]), history)
        end

        push!(history, name(robots[i]))
    end

    @testset "fresh names of reset robots are valid and unique in history" begin
        Random.seed!()
        for r in robots
            reset!(r)
            @test isname(name(r))
            @test !in(name(r), history)
            push!(history, name(r))
        end
    end
end


text:
################################################################################
#
#  Print functions
#
################################################################################
function _print_to_file(list::Vector{FieldsTower}, f::String; separator = '#')
  c = open(f, "w+")
  _print_to_file(list, c, separator = separator)
  close(c)
end

function _print_to_file(list::Vector{FieldsTower}, f::IOStream; separator = '#')
  for i = 1:length(list)
    K = list[i].field
    auts = list[i].generators_of_automorphisms
    embs = list[i].subfields
    p = K.pol
    p1 = fmpz[numerator(coeff(p, i)) for i = 0:degree(p)]
    Qx = parent(p)
    imgs = fmpq_poly[Qx(image_primitive_element(x)) for x in auts]
    imgs_coeff = Vector{fmpq}[fmpq[coeff(x1, i) for i = 0:degree(x1)] for x1 in imgs]
    embs_imgs = Tuple{fmpq_poly, fmpq_poly}[(codomain(x2).pol, parent(codomain(x2).pol)(image_primitive_element(x2))) for x2 in embs]
    embs_imgs_coeff = Vector{Tuple{Vector{fmpq}, Vector{fmpq}}}(undef, length(embs_imgs))
    for i = 1:length(embs_imgs)
      y, z = embs_imgs[i]
      embs_imgs_coeff[i] = ([coeff(y, j) for j = 0:degree(y)], [coeff(z, k) for k = 0:degree(z)])
    end
    Base.write(f, "$(p1)$(separator)$(imgs_coeff)$(separator)$(embs_imgs_coeff)\n")
  end
  return nothing
end

Base.convert(::Type{fmpq}, a::Rational{<: Integer}) = fmpq(a)

function _read_from_file(f::String)
  c = open(f, "r")
  v = _read_from_file(c)
  close(c)
  return v
end

function _read_from_file(f::IOStream)
  Qx, x = PolynomialRing(FlintQQ, "x")
  list = Vector{FieldsTower}(undef, countlines(f))
  seekstart(f)
  QQ = number_field(x-1, "b", cached = false, check = false)[1]
  for (k, l) in enumerate(eachline(f))
    a = _parse_fields_tower_blurb(l)
    #a = eval(Meta.parse(replace(l, '#' => ',')))
    #@assert a == olda
    K = number_field(Qx(a[1]), "a", cached = false, check = false)[1]
    auts = Vector{NfToNfMor}(undef, length(a[2]))
    for i = 1:length(auts)
      img = K(Qx(a[2][i]))
      auts[i] = NfToNfMor(K, K, img)
    end
    embs = Vector{NfToNfMor}(undef, length(a[3]))
    def_pol = Qx(a[3][1][1])
    if def_pol == K.pol
      Kcodomain = K
    else
      Kcodomain = number_field(def_pol, cached = false, check = false)[1]
    end
    embs[1] = NfToNfMor(QQ, Kcodomain, Kcodomain(Qx(a[3][1][2])))
    for i = 2:(length(a[3]))
      Kdomain = Kcodomain
      def_pol = Qx(a[3][i][1])
      if def_pol == K.pol
        Kcodomain = K
      else
        Kcodomain = number_field(def_pol, cached = false, check = false)[1]
      end
      embs[i] = NfToNfMor(Kdomain, Kcodomain, Kcodomain(Qx(a[3][i][2])))
    end
    list[k] = FieldsTower(K, auts, embs)
  end
  return list
end

function _parse_as_fmpq_vector(s)
  res = fmpq[]
  i = 1
  while s[i] != '['
    i += 1
  end
  i += 1
  isrational = false
  j = i
  last_parsed = 0
  while true
    while s[j] == ' '
      j += 1
    end
    if s[j] == ']'
      z = parse(fmpz, s[i:j-1])
      if isrational
        push!(res, last_parsed//z)
        isrational = false
      else
        push!(res, fmpq(z))
      end
      break
    end
    if s[j] == ','
      z = parse(fmpz, s[i:j-1])
      if isrational
        push!(res, last_parsed//z)
        isrational = false
      else
        push!(res, fmpq(z))
      end
      i = j + 1
    end
    if s[j] == '/'
      last_parsed = parse(fmpz, s[i:j-1])
      isrational = true
      j += 2
      i = j
    end
    j += 1
  end
  res
end

function _parse_array_of_array_of_fmpq(s)
  res = Vector{fmpq}[]
  i = 1
  while s[i] != '['
    i += 1
  end
  i += 1
  while s[i] != ']'
    k = _find_end_of_array_not_nested(s, i)
    v = _parse_as_fmpq_vector(s[i:k])
    i = k + 1
    push!(res, v)
  end
  return res
end

function _skip_head(s)
  i = 1
  while s[i] != '['
    i += 1
  end
  return i
end

function _parse_array_of_tuples_of_arrays(s)
  res = Tuple{Vector{fmpq}, Vector{fmpq}}[]
  i = _skip_head(s)
  i += 1
  end_sbraket = _find_closing_sbraket(s, i)
  while i < end_sbraket
    while s[i] != '('
      i += 1
    end
    i += 1
    while s[i] != '['
      i += 1
    end
    i += 1
    k = _find_closing_sbraket(s, i)
    vec1 = _parse_as_fmpq_vector(s[i-1:k-1])
    i = k
    while s[i] != ','
      i += 1
    end
    i += 1
    while s[i] != '['
      i += 1
    end
    i += 1
    k = _find_closing_sbraket(s, i)
    vec2 = _parse_as_fmpq_vector(s[i-1:k-1])
    i = k
    while s[i] != ')'
      i += 1
    end
    i += 1
    push!(res, (vec1, vec2))
    while s[i] != ','
      i += 1
      if i == end_sbraket
        return res
      end
    end
    i += 1
  end
  res
end

function _find_closing_sbraket(s, start = 1)
  @assert s[start - 1] == '['
  i = start
  open = 1
  while open != 0
    if s[i] == '['
      open += 1
    elseif s[i] == ']'
      open -= 1
    end
    i += 1
  end
  return i
end

function _find_closing_braket(s, start = 1)
  i = start
  open = 1
  while open != 0
    if s[i] == '('
      open += 1
    elseif s[i] == ')'
      open -= 1
    end
    i += 1
  end
  return i
end

# return value will point to the closing ]
function _find_end_of_array_not_nested(s, start = 1)
  i = start + 1
  while s[i] != ']'
    i += 1
  end
  return i
end

function _parse_fields_tower_blurb(s)
  s_split = split(s, '#')
  @assert length(s_split) == 3
  res1 = _parse_as_fmpq_vector(String(s_split[1]))
  res2 = _parse_array_of_array_of_fmpq(String(s_split[2]))
  res3 = _parse_array_of_tuples_of_arrays(String(s_split[3]))
  return res1, res2, res3
end


text:
#' ---
#' title: Conic Programming
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' This tutorial is aimed at providing a simplistic introduction to conic programming using JuMP.

#' ## What is a Cone?
#' A subset $C$ of a vector space $V$ is a cone if $\forall x \in C$ and positive scalars $\alpha$,
#' the product $\alpha x \in C$. A cone C is a convex cone if $\alpha x + \beta y \in C$,
#' for any positive scalars $\alpha, \beta$, and any $x, y \in C$.

#' ## Conic Programming
#' Conic programming problems are convex optimization problems in which a convex function is minimized
#' over the intersection of an affine subspace and a convex cone.
#' An example of a conic-form minimization problems, in the primal form is:

#' $$
#' \begin{align}
#' & \min_{x \in \mathbb{R}^n} & a_0^T x + b_0 \\
#' & \;\;\text{s.t.} & A_i x + b_i & \in \mathcal{C}_i & i = 1 \ldots m
#' \end{align}
#' $$

#' The corresponding dual problem is:

#' $$
#' \begin{align}
#' & \max_{y_1, \ldots, y_m} & -\sum_{i=1}^m b_i^T y_i + b_0 \\
#' & \;\;\text{s.t.} & a_0 - \sum_{i=1}^m A_i^T y_i & = 0 \\
#' & & y_i & \in \mathcal{C}_i^* & i = 1 \ldots m
#' \end{align}
#' $$

#' where each $\mathcal{C}_i$ is a closed convex cone and $\mathcal{C}_i^*$ is its dual cone.

#' ## Some of the Types of Cones Supported by JuMP

using JuMP
using ECOS
using LinearAlgebra
using Random

Random.seed!(1234);

#' By this point we have used quite a few different solvers.
#' To find out all the different solvers and their supported problem types, check out the
#' [solver table](http://www.juliaopt.org/JuMP.jl/v0.19.0/installation/#Getting-Solvers-1) in the docs.

#' ### Second-Order Cone
#' The Second-Order Cone (or Lorenz Cone) of dimension $n$ is of the form:

#' $$
#' Q^n = \{ (t,x) \in \mathbb{R}^\mbox{n} : t \ge ||x||_2 \}
#' $$

#' A Second-Order Cone rotated by $\pi/4$ in the $(x_1,x_2)$ plane is called a Rotated Second-Order Cone.
#' It is of the form:

#' $$
#' Q_r^n = \{ (t,u,x) \in \mathbb{R}^\mbox{n} : 2tu \ge ||x||_2^2, t,u \ge 0 \}
#' $$

#' These cones are represented in JuMP using the MOI sets `SecondOrderCone` and `RotatedSecondOrderCone`.

#' #### Example: Euclidean Projection on a Hyperplane
#' For a given point $u_{0}$ and a set $K$, we refer to any point $u \in K$
#' which is closest to $u_{0}$ as a projection of $u_{0}$ on $K$.
#' The projection of a point $u_{0}$ on a hyperplane $K = \{u | p' \cdot u = q\}$ is given by

#' $$
#' \begin{align*}
#' & \min & ||u - u_{0}|| \\
#' & \;\;\text{s.t.} & p' \cdot u = q  \\
#' \end{align*}
#' $$

u0 = rand(10)
p = rand(10)
q = rand();

#' We can model the above problem as the following conic program:

#' $$
#' \begin{align*}
#' & \min & t \\
#' & \;\;\text{s.t.} & p' \cdot u = q \\
#' & & (t, u - u_{0}) \in Q^{n+1}
#' \end{align*}
#' $$

#' On comparing this with the primal form of a conic problem we saw above,

#' $$
#' \begin{align*}
#' & x = (t , u) &\\
#' & a_0 = e_1 &\\
#' & b_0 = 0 &\\
#' & A_1 = (0, p) &\\
#' & b_1 = -q &\\
#' & C_1 = \mathbb{R}_- &\\
#' & A_2 = 1 &\\
#' & b_2 = -(0, u_0) &\\
#' & C_2 = Q^{n+1} &
#' \end{align*}
#' $$

#' Thus, we can obtain the dual problem as:

#' $$
#' \begin{align*}
#' & \max & y_1 + (0, u_0)^T y_2 \\
#' & \;\;\text{s.t.} & e_1 - (0,p)^T y_1 - y_2 = 0 \\
#' & & y_1 \in \mathbb{R}_- \\
#' & & y_2 \in Q^{n+1}
#' \end{align*}
#' $$

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, (u - u0)...] in SecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)

#+

@show objective_value(model);
@show value.(u);

#+

e1 = [1, zeros(10)...]
dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(dual_model, y1 <= 0)
@variable(dual_model, y2[1:11])
@objective(dual_model, Max, q * y1 + dot(vcat(0, u0), y2))
@constraint(dual_model, e1 - [0, p...] .* y1 - y2 .== 0)
@constraint(dual_model, y2 in SecondOrderCone())
optimize!(dual_model)

#+

@show objective_value(dual_model);

#' We can also have an equivalent formulation using a Rotated Second-Order Cone:

#' $$
#' \begin{align*}
#' & \min & t \\
#' & \;\;\text{s.t.} & p' \cdot u = q \\
#' & & (t, 1/2, u - u_{0})\in Q_r^{n+2}
#' \end{align*}
#' $$

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, 0.5, (u - u0)...] in RotatedSecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)

#+

@show value.(u);

#' The difference here is that the objective in the case of the Second-Order Cone is $||u - u_{0}||_2$,
#' while in the case of a Rotated Second-Order Cone is $||u - u_{0}||_2^2$.
#' However, the value of x is the same for both.

#' ### Exponential Cone

#' An Exponential Cone is a set of the form:

#' $$
#' K_{exp} = \{ (x,y,z) \in \mathbb{R}^3 : y \exp (x/y) \le z, y > 0 \}
#' $$

#' It is represented in JuMP using the MOI set `ExponentialCone`.

#' #### Example: Entropy Maximization
#' As the name suggests, the entropy maximization problem consists of maximizing the entropy function,
#' $H(x) = -x\log{x}$ subject to linear inequality constraints.

#' $$
#' \begin{align*}
#' & \max & - \sum_{i=1}^n x_i \log x_i \\
#' & \;\;\text{s.t.} & \mathbf{1}' x = 1 \\
#' & & Ax \leq b
#' \end{align*}
#' $$

#' We can model this problem using an exponential cone by using the following transformation:

#' $$
#' t\leq -x\log{x} \iff t\leq x\log(1/x)  \iff (1, x, t) \in K_{exp}
#' $$

#' Thus, our problem becomes,

#' $$
#' \begin{align*}
#' & \max & 1^Tt \\
#' & \;\;\text{s.t.} & Ax \leq b \\
#' & & 1^T x = 1 \\
#' & & (1, x_i, t_i) \in K_{exp} && \forall i = 1 \ldots n \\
#' \end{align*}
#' $$

n = 15;
m = 10;
A = randn(m, n);
b = rand(m, 1);

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b )
# Cannot use the exponential cone directly in JuMP, hence we use MOI to specify the set.
@constraint(model, con[i = 1:n], [1, x[i], t[i]] in MOI.ExponentialCone())

optimize!(model);

#+

@show objective_value(model);

#' ### Positive Semidefinite Cone
#' The set of Positive Semidefinite Matrices of dimension $n$ form a cone in $\mathbb{R}^n$.
#' We write this set mathematically as

#' $$
#' \mathcal{S}_{+}^n = \{ X \in \mathcal{S}^n \mid z^T X z \geq 0, \: \forall z\in \mathbb{R}^n \}.
#' $$

#' A PSD cone is represented in JuMP using the MOI sets
#' `PositiveSemidefiniteConeTriangle` (for upper triangle of a PSD matrix) and
#' `PositiveSemidefiniteConeSquare` (for a complete PSD matrix).
#' However, it is prefferable to use the `PSDCone` shortcut as illustrated below.

#' #### Example: Largest Eigenvalue of a Symmetric Matrix
#' Suppose $A$ has eigenvalues $\lambda_{1} \geq \lambda_{2} \ldots \geq \lambda_{n}$.
#' Then the matrix $t I-A$ has eigenvalues $t-\lambda_{1}, t-\lambda_{2}, \ldots, t-\lambda_{n}$.
#' Note that $t I-A$ is PSD exactly when all these eigenvalues are non-negative,
#' and this happens for values $t \geq \lambda_{1} .$
#' Thus, we can model the problem of finding the largest eigenvalue of a symmetric matrix as:

#' $$
#' \begin{align*}
#' \lambda_{1} = \max t \\
#' \text { s.t. } t I-A \succeq 0
#' \end{align*}
#' $$

#+ tangle = false

using LinearAlgebra
using SCS

A = [3 2 4;
     2 0 2;
     4 2 3]

model = Model(optimizer_with_attributes(SCS.Optimizer, "verbose" => 0))
@variable(model, t)
@objective(model, Min, t)
@constraint(model, t .* Matrix{Float64}(I, 3, 3) - A in PSDCone())

optimize!(model)

#+ tangle = false

@show objective_value(model);

#' ## Other Cones and Functions
#' For other cones supported by JuMP, check out the
#' [MathOptInterface Manual](http://www.juliaopt.org/MathOptInterface.jl/dev/apimanual/#Standard-form-problem-1).
#' A good resource for learning more about functions which can be modelled using cones is the
#' MOSEK Modeling Cookbook[[1]](#c1).

#' ### References
#'
#' 1. MOSEK Modeling Cookbook — MOSEK Modeling Cookbook 3.1. Available at: https://docs.mosek.com/modeling-cookbook/index.html.


text:
@static if VERSION >= v"1.1"
    using Pkg: instantiate
    instantiate()

    using Coverage: process_folder
    using Coverage.Codecov: submit
    using Coverage.LCOV: readfolder

    submit(
        filter(
            let prefix = joinpath(pwd(), "src", "")
                coverage -> startswith(coverage.filename, prefix)
            end,
            readfolder(".")
        )
    )
end


text:
# Zygote fill has issues with non-numbers
ZygoteRules.@adjoint function fill(x::T, dims...) where {T}
    return ZygoteRules.pullback(x, dims...) do x, dims...
        return reshape([x for i in 1:prod(dims)], dims)
    end
end

## Uniform ##

ZygoteRules.@adjoint function uniformlogpdf(a, b, x)
    diff = b - a
    T = typeof(diff)
    if a <= x <= b && a < b
        l = -log(diff)
        da = 1/diff^2
        return l, Δ -> (da * Δ, -da * Δ, zero(T) * Δ)
    else
        n = T(NaN)
        return n, Δ -> (n, n, n)
    end
end

ZygoteRules.@adjoint function Distributions.Uniform(args...)
    return ZygoteRules.pullback(TuringUniform, args...)
end

## Beta ##

function _betalogpdfgrad(α, β, x)
    di = digamma(α + β)
    dα = log(x) - digamma(α) + di
    dβ = log(1 - x) - digamma(β) + di
    dx = (α - 1)/x + (1 - β)/(1 - x)
    return (dα, dβ, dx)
end
ZygoteRules.@adjoint function betalogpdf(α::Real, β::Real, x::Number)
    return betalogpdf(α, β, x), Δ -> (Δ .* _betalogpdfgrad(α, β, x))
end

## Gamma ##

function _gammalogpdfgrad(k, θ, x)
    dk = -digamma(k) - log(θ) + log(x)
    dθ = -k/θ + x/θ^2
    dx = (k - 1)/x - 1/θ
    return (dk, dθ, dx)
end
ZygoteRules.@adjoint function gammalogpdf(k::Real, θ::Real, x::Number)
    return gammalogpdf(k, θ, x), Δ -> (Δ .* _gammalogpdfgrad(k, θ, x))
end

## Chisq ##

function _chisqlogpdfgrad(k, x)
    hk = k/2
    d = digamma(hk)
    dk = (-log(oftype(hk, 2)) - d + log(x))/2
    dx = (hk - 1)/x - one(hk)/2
    return (dk, dx)
end
ZygoteRules.@adjoint function chisqlogpdf(k::Real, x::Number)
    return chisqlogpdf(k, x), Δ -> (Δ .* _chisqlogpdfgrad(k, x))
end

## FDist ##

function _fdistlogpdfgrad(v1, v2, x)
    temp1 = v1 * x + v2
    temp2 = log(temp1)
    vsum = v1 + v2
    temp3 = vsum / temp1
    temp4 = digamma(vsum / 2)
    dv1 = (log(v1 * x) + 1 - temp2 - x * temp3 - digamma(v1 / 2) + temp4) / 2
    dv2 = (log(v2) + 1 - temp2 - temp3 - digamma(v2 / 2) + temp4) / 2
    dx = v1 / 2 * (1 / x - temp3) - 1 / x
    return (dv1, dv2, dx)
end
ZygoteRules.@adjoint function fdistlogpdf(v1::Real, v2::Real, x::Number)
    return fdistlogpdf(v1, v2, x), Δ -> (Δ .* _fdistlogpdfgrad(v1, v2, x))
end

## TDist ##

function _tdistlogpdfgrad(v, x)
    dv = (digamma((v + 1) / 2) - 1 / v - digamma(v / 2) - log(1 + x^2 / v) + x^2 * (v + 1) / v^2 / (1 + x^2 / v)) / 2
    dx = -x * (v + 1) / (v + x^2)
    return (dv, dx)
end
ZygoteRules.@adjoint function tdistlogpdf(v::Real, x::Number)
    return tdistlogpdf(v, x), Δ -> (Δ .* _tdistlogpdfgrad(v, x))
end

## Binomial ##

ZygoteRules.@adjoint function binomlogpdf(n::Int, p::Real, x::Int)
    return binomlogpdf(n, p, x),
        Δ->(nothing, Δ * (x / p - (n - x) / (1 - p)), nothing)
end

## Poisson ##

ZygoteRules.@adjoint function poislogpdf(v::Real, x::Int)
    return poislogpdf(v, x),
        Δ->(Δ * (x/v - 1), nothing)
end

## PoissonBinomial ##

# Zygote loads ForwardDiff, so this dummy adjoint should never be needed.
# The adjoint that is used for `poissonbinomial_pdf_fft` is defined in `src/zygote_forwarddiff.jl`
# ZygoteRules.@adjoint function poissonbinomial_pdf_fft(x::AbstractArray{T}) where T<:Real
#     error("This needs ForwardDiff. `using ForwardDiff` should fix this error.")
# end


text:
"""
    Sinograms module
"""
module Sinograms

const RealU = Number # Real or Unitful

include("fbp2.jl")
include("fbp2_sino_filter.jl")
include("fbp2_back.jl")
include("fbp2_back_fan.jl")
include("fbp_ramp.jl")
include("fbp2_window.jl")
include("fbp2_sino_weight.jl")

include("fbp-par.jl")

end # module


text:
export raw_data_available

function raw_data_available(ROOT::String; format::String="mseed")
    files = glob("*"*format,ROOT)
    N = length(files)
    NET = Array{String,1}(undef,N)
    STA = Array{String,1}(undef,N)
    LOC = Array{String,1}(undef,N)
    CHAN = Array{String,1}(undef,N)
    STARTTIME = Array{DateTime,1}(undef,N)
    ENDTIME = Array{DateTime,1}(undef,N)

    for ii = 1:N
        path, mseed = splitdir(files[ii])
        chan, LOC[ii], starttime, endtime, _ = string.(split(mseed,'.'))
        path, CHAN[ii] = splitdir(path)
        path, STA[ii] = splitdir(path)
        path, NET[ii] = splitdir(path)
        STARTTIME[ii] = DateTime(starttime,"yyyymmddTHHMMSSZ")
        ENDTIME[ii] = DateTime(endtime,"yyyymmddTHHMMSSZ")
    end

    df = DataFrame(NET = NET, STA = STA, LOC = LOC, CHAN = CHAN,
               STARTTIME = STARTTIME, ENDTIME = ENDTIME, FILE = files)
    return df
end

function mseed2jld2(df::DataFrame, DIR::String)

    # check if base dir is in file
    if isdir(DIR) == false
        mkpath(DIR)
    end

    # group stations by network and station
    for station in groupby(df, [:NET, :STA, :CHAN])
        net = unique(station[:NET])[1]
        sta = unique(station[:STA])[1]
        chan = unique(station[:CHAN])[1]
        NET = joinpath(DIR,net)

        if isdir(NET) == false
            mkpath(NET)
        end

        # create JLD2 file and save mseed
        filename = joinpath(NET, "$net.$sta.jld2")
        file = jldopen(filename, "a+")
        group = JLD2.Group(file, chan)
        println("Adding station $net.$sta.$chan $(Dates.now())")

        for day in eachrow(station)
            st = readmseed(day[:FILE])
            group[Dates.format(day[:STARTTIME],"yyyy.mm.dd")] = st
        end

        close(file)
    end
end

function readjld2(filename::String, chan::String, starttime::DateTime,
                  endtime::DateTime)
    file = jldopen(filename,"r")
    days = keys(file[chan])
    dates = [DateTime(z,"yyyy.mm.dd") for z in days]
    ind = findall(x -> (x >= starttime) & (x <= endtime), dates)
    days = days[ind]
    data = Array{SeisData}(undef,length(days))
    for (ii,day) in enumerate(days)
        data[ii] = file[chan][day]
    end
    close(file)
    return data
end

function date_range(dates::AbstractArray, starttime::DateTime, endtime::DateTime)
    return [d for d in dates if (d >= starttime) & (d <= endtime)]
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export center

"""
    center(;kwargs...)
    center(children::Any;kwargs...)
    center(children_maker::Function;kwargs...)

A Center component.
Centers content vertically and horizontally. For more information, see: https://mantine.dev/core/center/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Content that should be centered vertically and horizontally
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `inline` (Bool; optional): Set to true to use inline-flex instead of flex
- `style` (Dict; optional): Inline style override
"""
function center(; kwargs...)
        available_props = Symbol[:children, :id, :className, :inline, :style]
        wild_props = Symbol[]
        return Component("center", "Center", "dash_mantine_components", available_props, wild_props; kwargs...)
end

center(children::Any; kwargs...) = center(;kwargs..., children = children)
center(children_maker::Function; kwargs...) = center(children_maker(); kwargs...)


text:
# Set - unordered collection of unique elements

# Creating sets

s = Set([11, 23, 34, 45]);
println(s);
println(typeof(s));

# Takes care of duplicates
s1 = Set([23, 23, 45, 45, 13]);
println(s1);

# Any types
s2 = Set([[1, 2], [3, 4]]);
println(s2);

s3 = Set(1:3);

# Sets must contain the same type of object

# Union

println(union(s, s1));

# Adding element with union

union!(s3, 6:8);
println(s3);

# Intersection

println(intersect(s, s1));

# Set difference - elements of first set without elements of 2nd set

diff = setdiff(s, s1);
println(diff);

# Check for subsets

isSubset = issubset(s1, s);
println(isSubset);

# Set Equality

println(issetequal(s1, s2));

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
macro unreachable(message="This line of code cannot be reached. Please file an issue.")
    quote
        error($(esc(message)))
    end
end


text:
include_ddp()

# Model
include_model("pendulum")
function f(model::Pendulum, x, u, w)
	mass = model.mass + w[1]
    @SVector [x[2],
              (u[1] / ((mass * model.lc * model.lc))
                - model.g * sin(x[1]) / model.lc
                - model.b * x[2] / (mass * model.lc * model.lc))]
end

n, m, d = 2, 1, 1
model = Pendulum{Midpoint, FixedTime}(n, m, d, 1.0, 0.1, 0.5, 9.81)

# Time
T = 101
h = 0.025

# Initial conditions, controls, disturbances
x1 = [0.0, 0.0]
xT = [π, 0.0] # goal state
ū = [1.0 * rand(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
Q = [(t < T ? Diagonal(1.0e-1 * ones(model.n))
        : Diagonal(1.0 * ones(model.n))) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]

R = [Diagonal(1.0e-2 * ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Problem
prob = problem_data(model, obj, copy(x̄), copy(ū), w, h, T)

# Solve
@time ddp_solve!(prob,
    max_iter = 100, verbose = true)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Visualize
using Plots
plot(t, hcat(x...)', label = "")
plot(hcat(u..., u[end])', linetype = :steppost)

# Simulate policy
include(joinpath(@__DIR__, "simulate.jl"))

# Model
model_sim = model
x1_sim = copy(x1)
T_sim = 10 * T
w_sim = [[10.0] for t = 1:T-1]

# Time
tf = h * (T - 1)
t = range(0, stop = tf, length = T)
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Policy
K = [K for K in prob.p_data.K]

# Simulate
x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate_linear_feedback(
	model_sim,
	K,
    x̄, ū,
	Q, R,
	T_sim, h,
	x1_sim,
	w_sim)

# Visualize
idx = (1:2)
plot(t, hcat(x̄...)[idx, :]',
    width = 2.0, color = :black, label = "")
plot!(t_sim, hcat(x_ddp...)[idx, :]',
    width = 1.0, color = :magenta, label = "")

plot(t, hcat(u..., u[end])',
	width = 2.0, color = :black, label = "",
	linetype = :steppost)
plot!(t_sim, hcat(u_ddp..., u_ddp[end])',
	width = 1.0, color = :magenta, label = "",
	linetype = :steppost)


text:
function power_balance_equation(cur_reactor::AbstractReactor)
  cur_B_power = cur_reactor.mode_scaling[:B]

  cur_R_power = alpha_R_star(cur_reactor)

  cur_I_power = alpha_I_star(cur_reactor)

  cur_G_T = G_PB(cur_reactor)

  cur_G_T /= K_PB(cur_reactor)

  cur_equation = Equation(
    cur_R_power,
    cur_B_power,
    cur_I_power,
    cur_G_T
  )

  cur_equation
end


text:
using Printf

function middle3(n::Integer)
    l = ndigits(n)
    iseven(l) && error("n must have an odd number of digits")
    l < 3 && error("n must have 3 digits or more")
    mid = (l + 1) ÷ 2
    abs(n) ÷ 10^(mid-2) % 10^3
end

for n = [123, 12345, 1234567, 987654321, 10001, -10001, -123,
         -100, 100, -12345, 1, 2, -1, -10, 2002, -2002, 0]
    @printf("%10d -> %s\n", n, try middle3(n) catch e e.msg end)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module CI_Test

"""
    TerribleGreek()

Map lowercase latin characters to symbols that kinda resemble them, ish...
"""
function TerribleGreek()
    return Dict('a' => 'α',
                'b' => 'β',
                'c' => '≪',
                'd' => 'δ',
                'e' => '∈',
                'f' => 'ɾ',
                'g' => 'ʔ',
                'h' => 'ƕ',
                'i' => 'ι',
                'j' => 'ȷ',
                'k' => 'k',
                'l' => 'ł',
                'm' => 'μ',
                'n' => 'π',
                'o' => 'ø',
                'p' => 'ɤ',
                'q' => '♩',
                'r' => 'ρ',
                's' => '𝖲',
                't' => 'τ',
                'u' => '⟰',
                'v' => '∇',
                'w' => '≅',
                'x' => 'Ξ',
                'y' => '¥',
                'z' => 'ζ',)
end

end # module


text:
using Plots
using ADCME
using PyCall
using DelimitedFiles
include("../simulation.jl")
include("../neural_net_multi.jl")
np = pyimport("numpy")

const SRC_CONST = 86400.0 # turns seconds into days
const GRAV_CONST = 0.0    # gravity constant
const CONV_CONST = 2000

################# --- Run Simulation --- #################

Nz = 1              # 1 cell in z direction
Nx = 101            # 100 cell in x direction   {  assert: typeof( (Nx-1) / (Nx_coloc-1) ) == Int  }
Nt = 2300           # number of timesteps       {  assert: typeof( (Nt) / (Nt_coloc-1) ) == Int  }
Δx = 30.0           # = Δz; cell size in x and z direction [m]
Δt = 0.1            # timespet size in [d]

z = (1:Nz)*Δx|>collect        # array of z coordinates
x = (1:Nx)*Δx|>collect        # array of x coordinates
X, Z = np.meshgrid(x, z)      # mesh grid

ρw = 1000.0     # density if wetting fluid [m^3/d]
ρo = 800.0      # density if produciton fluid [m^3/d]
μw = 1.0        # viscosity if wetting fluid [mPa*s]
μo = 10.7       # viscosity if production fluid [mPa*s]

g = GRAV_CONST  # gravitational acceleration

flow = 0.005 * (1/Δx^2)/10.0 * SRC_CONST # injection rate [m^3/d]

qw = zeros(Nt, Nz, Nx)      # injection rate
qw[:, 1, 1] .= flow
qo = zeros(Nt, Nz, Nx)      # produciton rate
qo[:,1,end] .= -flow

sw0 = zeros(Nz, Nx)         # initial saturation [%]

ϕ = ones(Nz, Nx)            # "real" porosity field [%]
ϕ[1, :] = ([610.0, 615.0, 620.0, 623.0, 628.0, 631.0, 635.0, 637.0, 638.0,
            640.0, 108.0, 111.0, 113.0, 120.0, 123.0, 125.0, 130.0, 139.0,
            145.0, 151.0, 155.0, 162.0, 515.0, 511.0, 506.0, 503.0, 498.0,
            495.0, 490.0, 910.0, 901.0, 892.0, 883.0, 875.0, 867.0, 859.0,
            852.0, 845.0, 839.0, 710.0, 709.0, 707.0, 704.0, 700.0, 695.0,
            689.0, 682.0, 674.0, 665.0, 655.0, 644.0, 633.0, 622.0, 610.0,
            689.0, 682.0, 674.0, 665.0, 655.0, 644.0, 633.0, 622.0, 610.0,
            610.0, 610.0, 610.0, 610.0, 630.0, 640.0, 650.0, 660.0, 310.0,
            310.0, 310.0, 330.0, 320.0, 330.0, 330.0, 340.0, 345.0, 350.0,
            355.0, 360.0, 365.0, 369.0, 373.0, 377.0, 381.0, 384.0, 386.0,
            388.0, 689.0, 388.0, 386.0, 384.0, 381.0, 377.0, 373.0, 369.0,
            365.0, 360.0])./ CONV_CONST
K = ones(Nz, Nx)
K[1, :] = ϕ[1, :] .* CONV_CONST   # "real" permeability field [milidarcy]

println("\nRunning simulation...")

# turn all simulation parameters into tf.Tensors
tf_param_real = tfCtxGen(Nz, Nx, Δx, Nt, Δt, Z, X, ρw, ρo, μw, μo, K, g, ϕ, qw, qo, sw0, true)
# create computational graph for simulation
tf_sw_real, tf_p_real = imseq(tf_param_real)

sess = Session(); init(sess)                # initialise simulation session
sw_real = run(sess, tf_sw_real)             # get simulation result for saturation

sess.close()                                # close simulation session
println("\nSimulation session closed!")

################# --- Neural Network to approximate for simulation data --- #################

Nx_coloc = 101  # number of colocation points on x-axis
Nt_coloc = 51  # number of colocation points on time-axis

interpol_x = Int((Nx-1)/(Nx_coloc-1))
interpol_t = Int(Nt/(Nt_coloc-1))
coloc_x = 1:interpol_x:Nx           # colocation point indices in x
coloc_t = 1:interpol_t:Nt+1         # colocation point indices in time

x_dimless = range(0, 1, length=Nx)|>collect    # dimensionless spatial coordinates
t_dimless = range(0, 1, length=Nt+1)|>collect    # dimensionless time

X_star = repeat(x_dimless[coloc_x], Nt_coloc)                                # 1D array(Nx*(Nt+1)) of spatial mesh coordinates
T_star = collect(Iterators.flatten(repeat(t_dimless[coloc_t]', Nx_coloc)))   # 1D array(Nx*(Nt+1)) of time mesh

sw_real_star = collect(reshape(sw_real[coloc_t, 1, coloc_x]', Nt_coloc * Nx_coloc)) # flatten test data

ϵ = 2.5*10e-3                                   # 1/Pe: dimensionless diffusion coefficient
config = [20, 20, 20, 20, 20, 20, 20, 20, 1]    # neural network configuration

# get tf.Tensor struct of input parameters
tf_param_NN = tfCtxGenNN(X_star, T_star, sw_real_star, μw, μo, ϵ, config, Nx_coloc, Nt_coloc)

# neural network function
sw_nn, sw_t_nn, sw_xx_nn, f_w_x_nn = neural_net(tf_param_NN)

# loss function
loss, λ_1, λ_2 = get_loss(tf_param_NN,
                                sw_nn,
                                sw_t_nn,
                                sw_xx_nn,
                                f_w_x_nn,
                                true, true, # calibrating λ1 and λ2, not calibrating λ3
                                true) # including production data in train data

################# --- training --- #################

# Adam optimizer
println("\nTraining with Adam optimizer...")

optimizer_Adam = tf.compat.v1.train.AdamOptimizer(learning_rate=0.001)
train_op = optimizer_Adam.minimize(loss)

sess = Session(); init(sess)

epochs = 5000
losses = []

for i in 1:epochs
    run(sess, train_op)
    if i % 10 == 0
        println("Epoch: ", i, "\tLoss: ", run(sess, loss))
    end
    push!(losses, run(sess, loss))
end

# Scipy L-BFGS-B optimizer
#=
println("\nTraining with L-BFGS-B optimizer...")
i = 0
function cb()
    global i
    if i % 10 == 0
        println("Epoch: $i")
        #push!(losses, loss)
    end
    i += 1
end

opt = ScipyOptimizerInterface(loss; method="L-BFGS-B", options=Dict("maxiter"=> 3000,
                                                                    "maxfun"=> 3000,
                                                                    "ftol"=>1e-12,
                                                                    "gtol"=>1e-12))
ScipyOptimizerMinimize(sess, opt; loss_callback=cb) =#
println("\nTraining with BFGS optimizer...")
BFGS!(sess, loss, 10000)

name = "./model.mat"
save(sess, name)

println("\nTraining finished!")

################# --- Computing errors --- #################

# get final loss and λs
loss_final, λ1, λ2 = run(sess, [loss, λ_1, λ_2])

# get results
X_star_tf = constant(repeat(x_dimless, Nt+1))
T_star_tf = constant(collect(Iterators.flatten(repeat(t_dimless', Nx))))
sw_real_star = collect(reshape(sw_real[:, 1, :]', (Nt+1) * Nx))

sw_pred_star = run(sess, squeeze(fc(hcat(X_star_tf, T_star_tf),
                                            tf_param_NN.config,
                                            tf_param_NN.Theta)) + 1)
# compute error
error = sum(@. sqrt((sw_pred_star - sw_real_star)^2)) / sum(@. sqrt(sw_real_star^2))

println("\nFinal loss:\t", loss_final)
println("L2 error:\t", error)

sess.close()
println("\nTraining session closed!")

################# --- Plot results --- #################

println("\nPlotting results...")

sw_pred_plot = reshape(sw_pred_star, Nx, Nt+1)
sw_real_plot = reshape(sw_real_star, Nx, Nt+1)

# plot predicted and "real" saturation
for i in 1:Int((Nt)/5):Nt+1
    close("all")
    day = (i-1) * Δt
    println(day)
    Plots.plot(x, sw_real_plot[:, i], label="simulation")
    Plots.plot!(x, sw_pred_plot[:, i], label="NN", line=:dash)
    Plots.ylims!((0.0, 1.0))
    name = "./results/"*"$(day)"*".png"
    Plots.savefig(name)
end

# plot porosity
close("all")
Plots.plot(x, ϕ[1, :], label="Porosity real")
Plots.ylims!((0.0, 1.0))
name = "./results/poro_field.png"
Plots.savefig(name)

# plot permeability
close("all")
Plots.plot(x, K[1, :], label="Permeability real")
Plots.ylims!((0.0, 1000.0))
name = "./results/perm_field.png"
Plots.savefig(name)

# plot losses
close("all")
Plots.plot(1:epochs, losses, label="loss", yaxis=:log)
name = "./results/loss.png"
Plots.savefig(name)

# plot lambda
close("all")
Plots.plot((1:Nx_coloc)*interpol_x*Δt, λ1, label="λ1")
Plots.plot!((1:Nx_coloc)*interpol_x*Δt, λ2, label="λ2")
name = "./results/lamda.png"
Plots.savefig(name)

################# --- Saving results --- #################

open("./results/error.txt", "w") do file
    write(file, "Final loss:\t $loss_final")
    write(file, "\nL2 error:\t $error")
    write(file, "\nLambda 1:\t $λ1")
    write(file, "\nLambda 2:\t $λ2")
    write(file, "\nNumber of colocation points:\t $(Nx_coloc*Nt_coloc)")
    write(file, "\nAdam epochs:\t $(epochs)")
end

open("./results/sw_real.txt", "w") do file
    writedlm(file, sw_real_plot')
end

open("./results/sw_pred.txt", "w") do file
    writedlm(file, sw_pred_plot')
end

max_stars_repo_path: 6. calibration heterogeneous lambda vector/nn_prediction.jl


text:
include("../src/global/ImportFolders.jl")

import Config
import Client
import Connection
import Event
import Helper

function _onAfterResponse(ans::Connection.Answer)
  println("after-response\n", ans)
end
function _onBeforeResponse(data::Connection.Command, ans::Connection.Answer)
  ans.data = Array{Int, 1}(2)
  println("before-response\n", data, ans)
end

con = Client.create(ip"127.0.0.1", Config.val(:conServerPort))
Event.on(con.observer, Connection.EVENT_AFTER_REQUEST, _onAfterResponse)
Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, _onBeforeResponse)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
check_match(a, b) = string(a, base=2, pad=32)[17:32] == string(b, base=2, pad=32)[17:32] ? 1 : 0

function puzzle15(astart::Int=883, bstart::Int=879)
    a = copy(astart)
    b = copy(bstart)

    matches = 0
    for iter in 1:40_000_000
        a = rem(a * 16807, 2147483647)
        b = rem(b * 48271, 2147483647)
        matches += check_match(a, b)
    end
    partone = matches

    a = copy(astart)
    b = copy(bstart)
    matches = 0
    for iter in 1:5_000_000
        while true
            a = rem(a * 16807, 2147483647)
            rem(a, 4) == 0 ? break : nothing
        end
        while true
            b = rem(b * 48271, 2147483647)
            rem(b, 8) == 0 ? break : nothing
        end
        matches += check_match(a, b)
    end
    parttwo = matches

    return [partone, parttwo]
end


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

#############
# A sample model
#############
scen = 1
m = StructuredModel(num_scenarios=scen)
@variable(m, -100<=x[1:4]<=100)
@NLobjective(m, Min,  0.1*x[1]^2 + 0.2*x[2]^2 + 0.3*x[3]^2 + 0.4*x[4]^2 +
    + 1.1*x[1]*x[2] + 2.1*x[1]*x[3] + 2.2*x[1]*x[4] +
    + 2.1*x[2]*x[3] + 2.2*x[2]*x[4]
    + 3.1*x[3]*x[4])

for i in getLocalChildrenIds(m)
    bl = StructuredModel(parent=m, id=i)
    @variable(bl, -200<=y1[1:3]<=200)
    @NLconstraint(bl, -100<= 0.1*x[1] + 0.2*x[2]                          <= 100)
    @NLconstraint(bl, -101<=               1.1*x[2] + 1.2*x[3]            <= 101)
    @NLconstraint(bl, -102<=                          2.1*x[3] + 2.2*x[4] <= 102)

    @NLconstraint(bl, -201<= 10.1*x[1] + 10.2*x[2]                           - 10.3*y1[1] + 10.4*y1[2]               <= 201)
    @NLconstraint(bl, -202<=               11.1*x[2] + 11.2*x[3]             - 11.3*y1[1] + 11.4*y1[2] + 11.5*y1[3]  <= 202)
    @NLconstraint(bl, -203<=                           12.1*x[3] + 12.2*x[4]              + 12.4*y1[2] - 12.5*y1[3]  <= 203)

    @NLobjective(bl, Min,
        5.1*y1[1]*x[1] + 5.2*y1[1]*x[2] + 5.3*y1[1]*x[3] + 5.4*y1[1]*x[4]
    +                    6.2*y1[2]*x[2] +                  6.4*y1[2]*x[4]
    + 7.1*y1[3]*x[1] +                  7.3*y1[3]*x[3]
    + 8.1*y1[1]^2    + 8.2*y1[2]^2 + 8.3*y1[3]^2
    + 9.1*y1[1]*y1[2]+ 9.2*y1[1]*y1[3]
    + 10.1*y1[2]*y1[3])
end

solve(m, solver="PipsNlp")


text:
module WellLogs

using PyPlot
using Colors
using Distributions
using Plots
using Statistics
using DataFrames
using RecipesBase

export KhTrack, ResTrack, SwTrack, LithoTrack, PermTrack, VshPhieTrack, PoroTrack,
        Phiefec, SwFunction, PermWR, RhoPorosity, Vshale, AddColPetro, PetroPhysics,
        OGIP, OOIP, meshgrid, VDLTrack, CBLTrack, CCLTrack, Reg, CalTrack,FormationsTrack,
        VshTrack, PhieTrack

#Track functions
include("KhTrack.jl")
include("ResTrack.jl")
include("SwTrack.jl")
include("LithoTrack.jl")
include("PermTrack.jl")
include("VshTrack.jl")
include("PoroTrack.jl")
include("VdlTrack.jl")
include("CblTrack.jl")
include("CCLTrack.jl")
include("CalTrack.jl")
include("FormationsTrack.jl")

#Utility Functions
include("meshgrid.jl")
include("Reg.jl")

#Formula PetroPhysics
include("Phiefec.jl")
include("SwFunction.jl")
include("PermWR.jl")
include("RhoPorosity.jl")
include("Vshale.jl")

#Formula for include PetroPhysics estimates in well logs
include("AddColPetro.jl")
include("PetroPhysics.jl")
include("PropHist.jl")
include("PickettPlot.jl")

#Formula for volumetrics estimates
include("OGIP.jl")
include("OOIP.jl")

include("TopDistance.jl")

#Plot Recipes
include("IDWIRecipe.jl")
# Dictionaries

end # module


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_charging_cost(path::AbstractString, inputs::Dict, dfCharge::DataFrame, dfPrice::DataFrame, dfPower::DataFrame, setup::Dict)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	#calculating charging cost
 	dfChargingcost = DataFrame(Region = dfGen[!,:region], Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone], Cluster = dfGen[!,:cluster], AnnualSum = Array{Union{Missing,Float32}}(undef, G), )
	# the price is already US$/MWh, and dfPower and dfCharge is already in MW, so no scaling is needed
	i = 1
	dfChargingcost_ = (DataFrame([[names(dfCharge)]; collect.(eachrow(dfCharge))], [:column; Symbol.(axes(dfCharge, 1))])[4:T+3,i+1] .*
	DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
	inputs["omega"])
	if i in inputs["FLEX"]
		dfChargingcost_ = (DataFrame([[names(dfPower)]; collect.(eachrow(dfPower))], [:column; Symbol.(axes(dfPower, 1))])[4:T+3,i+1] .*
		DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
		inputs["omega"])
	end
	for i in 2:G
		if i in inputs["FLEX"]
			dfChargingcost_1 = (DataFrame([[names(dfPower)]; collect.(eachrow(dfPower))], [:column; Symbol.(axes(dfPower, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
			inputs["omega"])
		else
			dfChargingcost_1 = (DataFrame([[names(dfCharge)]; collect.(eachrow(dfCharge))], [:column; Symbol.(axes(dfCharge, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1].*
			inputs["omega"])
		end
		dfChargingcost_ = hcat(dfChargingcost_, dfChargingcost_1)
	end
	dfChargingcost = hcat(dfChargingcost, DataFrame(dfChargingcost_', :auto))
 	for i in 1:G
 		dfChargingcost[!,:AnnualSum][i] = sum(dfChargingcost[i,6:T+5])
 	end
	dfChargingcost_annualonly = dfChargingcost[!,1:5]
	CSV.write(joinpath(path, "ChargingCost.csv"), dfChargingcost_annualonly)
	return dfChargingcost
end
#=function write_charging_cost(path::AbstractString, inputs::Dict, dfCharge::DataFrame, dfPrice::DataFrame, dfPower::DataFrame, setup::Dict)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	#calculating charging cost
 	dfChargingcost = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone], Sum = Array{Union{Missing,Float32}}(undef, G))
	# the price is already US$/MWh, and dfPower and dfCharge is already in MW, so no scaling is needed

	for i in 1:G
		if i in inputs["FLEX"]
			dfChargingcost_1 = (DataFrame([[names(dfPower)]; collect.(eachrow(dfPower))], [:column; Symbol.(axes(dfPower, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1] .*
			inputs["omega"])
		else
			dfChargingcost_1 = (DataFrame([[names(dfCharge)]; collect.(eachrow(dfCharge))], [:column; Symbol.(axes(dfCharge, 1))])[4:T+3,i+1] .*
			DataFrame([[names(dfPrice)]; collect.(eachrow(dfPrice))], [:column; Symbol.(axes(dfPrice, 1))])[2:T+1,dfPower[i,:][:Zone]+1] .*
			inputs["omega"])
		end
		dfChargingcost = hcat(dfChargingcost, convert(DataFrame, dfChargingcost_1'))
	end

 	for i in 1:G
 		dfChargingcost[!,:Sum][i] = sum(dfChargingcost[i,6:T+5])
 	end
	CSV.write(joinpath(path, "ChargingCost.csv"), dfChargingcost)
	return dfChargingcost
end
=#

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
package ln

import (
	"math"
	"math/rand"
)

struct Vector
    x::Float64
    y::Float64
    z::Float64
end

type Vector struct {
	X, Y, Z float64
}

func RandomUnitVector() Vector {
	for {
		x := rand.Float64()*2 - 1
		y := rand.Float64()*2 - 1
		z := rand.Float64()*2 - 1
		if x*x+y*y+z*z > 1 {
			continue
		}
		return Vector{x, y, z}.Normalize()
	}
}

func (a Vector) Length() float64 {
	return math.Sqrt(a.X*a.X + a.Y*a.Y + a.Z*a.Z)
}

func (a Vector) Distance(b Vector) float64 {
	return a.Sub(b).Length()
}

func (a Vector) LengthSquared() float64 {
	return a.X*a.X + a.Y*a.Y + a.Z*a.Z
}

func (a Vector) DistanceSquared(b Vector) float64 {
	return a.Sub(b).LengthSquared()
}

func (a Vector) Dot(b Vector) float64 {
	return a.X*b.X + a.Y*b.Y + a.Z*b.Z
}

func (a Vector) Cross(b Vector) Vector {
	x := a.Y*b.Z - a.Z*b.Y
	y := a.Z*b.X - a.X*b.Z
	z := a.X*b.Y - a.Y*b.X
	return Vector{x, y, z}
}

func (a Vector) Normalize() Vector {
	d := a.Length()
	return Vector{a.X / d, a.Y / d, a.Z / d}
}

func (a Vector) Add(b Vector) Vector {
	return Vector{a.X + b.X, a.Y + b.Y, a.Z + b.Z}
}

func (a Vector) Sub(b Vector) Vector {
	return Vector{a.X - b.X, a.Y - b.Y, a.Z - b.Z}
}

func (a Vector) Mul(b Vector) Vector {
	return Vector{a.X * b.X, a.Y * b.Y, a.Z * b.Z}
}

func (a Vector) Div(b Vector) Vector {
	return Vector{a.X / b.X, a.Y / b.Y, a.Z / b.Z}
}

func (a Vector) AddScalar(b float64) Vector {
	return Vector{a.X + b, a.Y + b, a.Z + b}
}

func (a Vector) SubScalar(b float64) Vector {
	return Vector{a.X - b, a.Y - b, a.Z - b}
}

func (a Vector) MulScalar(b float64) Vector {
	return Vector{a.X * b, a.Y * b, a.Z * b}
}

func (a Vector) DivScalar(b float64) Vector {
	return Vector{a.X / b, a.Y / b, a.Z / b}
}

func (a Vector) Min(b Vector) Vector {
	return Vector{math.Min(a.X, b.X), math.Min(a.Y, b.Y), math.Min(a.Z, b.Z)}
}

func (a Vector) Max(b Vector) Vector {
	return Vector{math.Max(a.X, b.X), math.Max(a.Y, b.Y), math.Max(a.Z, b.Z)}
}

func (a Vector) MinAxis() Vector {
	x, y, z := math.Abs(a.X), math.Abs(a.Y), math.Abs(a.Z)
	switch {
	case x <= y && x <= z:
		return Vector{1, 0, 0}
	case y <= x && y <= z:
		return Vector{0, 1, 0}
	}
	return Vector{0, 0, 1}
}

func (a Vector) MinComponent() float64 {
	return math.Min(math.Min(a.X, a.Y), a.Z)
}

func (p Vector) SegmentDistance(v Vector, w Vector) float64 {
	l2 := v.DistanceSquared(w)
	if l2 == 0 {
		return p.Distance(v)
	}
	t := p.Sub(v).Dot(w.Sub(v)) / l2
	if t < 0 {
		return p.Distance(v)
	}
	if t > 1 {
		return p.Distance(w)
	}
	return v.Add(w.Sub(v).MulScalar(t)).Distance(p)
}


text:
LogisticClassifier_ = SKLM.LogisticRegression
@sk_clf mutable struct LogisticClassifier <: MMI.Probabilistic
    penalty::String            = "l2"::(_ in ("l1", "l2", "elasticnet", "none"))
    dual::Bool                 = false
    tol::Float64               = 1e-4::(_ > 0)
    C::Float64                 = 1.0::(_ > 0)
    fit_intercept::Bool        = true
    intercept_scaling::Float64 = 1.0::(_ > 0)
    class_weight::Any          = nothing
    random_state::Any          = nothing
    solver::String             = "lbfgs"::(_ in ("lbfgs", "newton-cg", "liblinear", "sag", "saga"))
    max_iter::Int              = 100::(_ > 0)
    multi_class::String        = "auto"::(_ in ("ovr", "multinomial", "auto"))
    verbose::Int               = 0
    warm_start::Bool           = false
    n_jobs::Option{Int}        = nothing
    l1_ratio::Option{Float64}  = nothing::(_ === nothing || 0 ≤ _ ≤ 1)
end
MMI.fitted_params(m::LogisticClassifier, (f, _, _)) = (
    classes   = f.classes_,
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model(LogisticClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Logistic regression classifier."
    )

# ============================================================================
LogisticCVClassifier_ = SKLM.LogisticRegressionCV
@sk_clf mutable struct LogisticCVClassifier <: MMI.Probabilistic
    Cs::Union{Int,AbstractVector{Float64}} = 10::((_ isa Int && _ > 0) || all(_ .> 0))
    fit_intercept::Bool        = true
    cv::Any                    = 5
    dual::Bool                 = false
    penalty::String            = "l2"::(_ in ("l1", "l2", "elasticnet", "none"))
    scoring::Any               = nothing
    solver::String             = "lbfgs"::(_ in ("lbfgs", "newton-cg", "liblinear", "sag", "saga"))
    tol::Float64               = 1e-4::(_ > 0)
    max_iter::Int              = 100::(_ > 0)
    class_weight::Any          = nothing
    n_jobs::Option{Int}        = nothing
    verbose::Int               = 0
    refit::Bool                = true
    intercept_scaling::Float64 = 1.0::(_ > 0)
    multi_class::String        = "auto"::(_ in ("ovr", "multinomial", "auto"))
    random_state::Any          = nothing
    l1_ratios::Option{AbstractVector{Float64}}=nothing::(_ === nothing || all(0 .≤ _ .≤ 1))
end
MMI.fitted_params(m::LogisticCVClassifier, (f, _, _)) = (
    classes     = f.classes_,
    coef        = f.coef_,
    intercept   = m.fit_intercept ? f.intercept_ : nothing,
    Cs          = f.Cs_,
    l1_ratios   = ifelse(m.penalty == "elasticnet", f.l1_ratios_, nothing),
    coefs_paths = f.coefs_paths_,
    scores      = f.scores_,
    C           = f.C_,
    l1_ratio    = f.l1_ratio_
    )
metadata_model(LogisticCVClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Logistic regression classifier with internal cross-validation."
    )

# ============================================================================
PassiveAggressiveClassifier_ = SKLM.PassiveAggressiveClassifier
@sk_clf mutable struct PassiveAggressiveClassifier <: MMI.Deterministic
    C::Float64            = 1.0::(_ > 0)
    fit_intercept::Bool   = true
    max_iter::Int         = 100::(_ > 0)
    tol::Float64          = 1e-3::(_ > 0)
    early_stopping::Bool  = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int = 5::(_ > 0)
    shuffle::Bool         = true
    verbose::Int          = 0
    loss::String          = "hinge"::(_ in ("hinge", "squared_hinge"))
    n_jobs::Option{Int}   = nothing
    random_state::Any     = 0
    warm_start::Bool      = false
    class_weight::Any     = nothing
    average::Bool         = false
end
MMI.fitted_params(m::PassiveAggressiveClassifier, (f, _, _)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model(PassiveAggressiveClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Passive aggressive classifier."
    )

# ============================================================================
PerceptronClassifier_ = SKLM.Perceptron
@sk_clf mutable struct PerceptronClassifier <: MMI.Deterministic
    penalty::Option{String} = nothing::(_ === nothing || _ in ("l2", "l1", "elasticnet"))
    alpha::Float64          = 1e-4::(_ > 0)
    fit_intercept::Bool     = true
    max_iter::Int           = 1_000::(_ > 0)
    tol::Option{Float64}    = 1e-3
    shuffle::Bool           = true
    verbose::Int            = 0
    eta0::Float64           = 1.0::(_ > 0)
    n_jobs::Option{Int}     = nothing
    random_state::Any       = 0
    early_stopping::Bool    = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int   = 5::(_ > 0)
    class_weight::Any       = nothing
    warm_start::Bool        = false
end
MMI.fitted_params(m::PerceptronClassifier, (f, _, _)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model(PerceptronClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Perceptron classifier."
    )

# ============================================================================
RidgeClassifier_ = SKLM.RidgeClassifier
@sk_clf mutable struct RidgeClassifier <: MMI.Deterministic
    alpha::Float64        = 1.0
    fit_intercept::Bool   = true
    normalize::Bool       = false
    copy_X::Bool          = true
    max_iter::Option{Int} = nothing::(_ === nothing || _ > 0)
    tol::Float64          = 1e-3::(arg>0)
    class_weight::Any     = nothing
    solver::String        = "auto"::(arg in ("auto","svd","cholesky","lsqr","sparse_cg","sag","saga"))
    random_state::Any     = nothing
end
MMI.fitted_params(m::RidgeClassifier, (f, _, _)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model(RidgeClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Ridge regression classifier."
    )

# ============================================================================
RidgeCVClassifier_ = SKLM.RidgeClassifierCV
@sk_clf mutable struct RidgeCVClassifier <: MMI.Deterministic
    alphas::AbstractArray{Float64} = [0.1,1.0,10.0]::(all(0 .≤ _))
    fit_intercept::Bool   = true
    normalize::Bool       = false
    scoring::Any          = nothing
    cv::Int               = 5
    class_weight::Any     = nothing
    store_cv_values::Bool = false
end
MMI.fitted_params(m::RidgeCVClassifier, (f, _, _)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model(RidgeCVClassifier,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Ridge regression classifier."
    )

# ============================================================================
SGDClassifier_ = SKLM.SGDClassifier
@sk_clf mutable struct SGDClassifier <: MMI.Deterministic
    loss::String          = "hinge"::(_ in ("hinge", "log", "modified_huber", "squared_hinge", "perceptron", "squared_loss", "huber", "epsilon_insensitive", "squared_epsilon_insensitive"))
    penalty::String       = "l2"::(_ in ("l1", "l2", "elasticnet", "none"))
    alpha::Float64        = 1e-4::(_ > 0)
    l1_ratio::Float64     = 0.15::(0 ≤ _ ≤ 1)
    fit_intercept::Bool   = true
    max_iter::Int         = 1_000::(_ > 0)
    tol::Option{Float64}  = 1e-3::(_ === nothing || _ > 0)
    shuffle::Bool         = true
    verbose::Int          = 0
    epsilon::Float64      = 0.1::(_ > 0)
    n_jobs::Option{Int}   = nothing
    random_state::Any     = nothing
    learning_rate::String = "optimal"::(_ in ("constant", "optimal", "invscaling", "adaptive"))
    eta0::Float64         = 0.0::(_ ≥ 0)
    power_t::Float64      = 0.5::(_ > 0)
    early_stopping::Bool  = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int = 5::(_ > 0)
    class_weight::Any     = nothing
    warm_start::Bool      = false
    average::Bool         = false
end
ProbabilisticSGDClassifier_ = SKLM.SGDClassifier
@sk_clf mutable struct ProbabilisticSGDClassifier <: MMI.Probabilistic
    loss::String          = "log"::(_ in ("log", "modified_huber")) # only those -> predict proba
    penalty::String       = "l2"::(_ in ("l1", "l2", "elasticnet", "none"))
    alpha::Float64        = 1e-4::(_ > 0)
    l1_ratio::Float64     = 0.15::(0 ≤ _ ≤ 1)
    fit_intercept::Bool   = true
    max_iter::Int         = 1_000::(_ > 0)
    tol::Option{Float64}  = 1e-3::(_ === nothing || _ > 0)
    shuffle::Bool         = true
    verbose::Int          = 0
    epsilon::Float64      = 0.1::(_ > 0)
    n_jobs::Option{Int}   = nothing
    random_state::Any     = nothing
    learning_rate::String = "optimal"::(_ in ("constant", "optimal", "invscaling", "adaptive"))
    eta0::Float64         = 0.0::(_ ≥ 0)
    power_t::Float64      = 0.5::(_ > 0)
    early_stopping::Bool  = false
    validation_fraction::Float64 = 0.1::(0 < _ < 1)
    n_iter_no_change::Int = 5::(_ > 0)
    class_weight::Any     = nothing
    warm_start::Bool      = false
    average::Bool         = false
end
MMI.fitted_params(m::SGDClassifier, (f,_,_)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
# duplication to avoid ambiguity that julia doesn't like
MMI.fitted_params(m::ProbabilisticSGDClassifier, (f,_,_)) = (
    coef      = f.coef_,
    intercept = ifelse(m.fit_intercept, f.intercept_, nothing)
    )
metadata_model.((SGDClassifier,ProbabilisticSGDClassifier),
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    descr   = "Linear classifier with stochastic gradient descent training."
    )


text:
import Base.Random: rand, zeros

export Population, popsize, popfits

@doc """A group of genotypes.
"""
type Population
  genotypes::Vector{Genotype}
end

Population(p::Population) = Population([Genotype(g) for g = p.genotypes])

@doc """rand(::Type{Population}, ls::Landscape, popsize::Int64)

Returns a random population of genotypes from landscape ls
"""
rand(::Type{Population}, ls::Landscape, popsize::Int64) =
    Population([rand(Genotype, ls) for _ = 1:popsize])

@doc """zeros(::Type{Population}, ls::Landscape, popsize::Int64)

Returns a population of zero genotypes from landscape ls
"""
zeros(::Type{Population}, ls::Landscape, popsize::Int64) =
    Population([zeros(Genotype, ls) for _ = 1:popsize])

@doc """popsize(p::Population)

Returns the size of the population.
"""
popsize(p::Population) = length(p.genotypes)

@doc """popfits(p::Population)

Returns a vector of the fitnesses of population p
"""
popfits(p::Population) = Float64[fitness(g) for g = p.genotypes]


text:
"""
    BiPerplex{T <: Real} = Perplexification{Prplex{T}}

A bi-perplex number is a perplexification of the perplex numbers.
"""
const BiPerplex{T <: Real} = Perplexification{Prplex{T}}

function BiPerplex(z::Prplex{T}) where T <: Real
    BiPerplex{T}(z, zero(Prplex{T}))
end

function BiPerplex(a::T, b::T, c::T, d::T) where T <: Real
    BiPerplex{T}(Prplex{T}(a, b), Prplex{T}(c, d))
end

function BiPerplex(a::T, b::T, c::T) where T <: Real
    BiPerplex{T}(Prplex{T}(a, b), Prplex{T}(c, zero(T)))
end

function BiPerplex(a::T, b::T) where T <: Real
    BiPerplex{T}(Prplex{T}(a, b), zero(Prplex{T}))
end

function BiPerplex(a::T) where T <: Real
    BiPerplex{T}(Prplex{T}(a, zero(T)), zero(Prplex{T}))
end

function BiPerplex(a::Real, b::Real)
    BiPerplex(promote(a, b)...)
end

function BiPerplex(a::Real, b::Real, c::Real)
    BiPerplex(promote(a, b, c)...)
end

function BiPerplex(a::Real, b::Real, c::Real, d::Real)
    BiPerplex(promote(a, b, c, d)...)
end

"""
    unreal(z::BiPerplex)

The unreal part of a bi-perplex number is a 3-vector.
"""
function unreal(z::BiPerplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::BiPerplex)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", R: ")
    print(io, z.l.r)
    print(io, ", Q: ")
    print(io, z.r.l)
    print(io, ", QR: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{BiPerplex{T}}) where T <: Real
    BiPerplex{T}(random(Prplex{T}), random(Prplex{T}))
end

function barR(z::BiPerplex{T}) where T <: Real
    BiPerplex{T}(barR(z.l), barR(z.r))
end

function barQ(z::BiPerplex{T}) where T <: Real
    BiPerplex{T}(z.l, -(z.r))
end

function quadrance(z::BiPerplex)
    (z * barQ(z)).l
end

function norm(z::BiPerplex)
    abs2(quadrance(z))
end

function iszerodivisor(z::BiPerplex)
    iszerodivisor(quadrance(z))
end

function inv(z::BiPerplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    barQ(z) * barR(q) / abs2(q)
end


text:
module Random

using Random
using PyCall
using ..Core

const _random = PyNULL()

include("PRNG.jl")
include("distributions.jl")
include("wrappers.jl")

function __init__()
  copy!(_random, pyimport("jax.random"))
end

end