text:
#################################################################
######################## BELIEF #################################
#################################################################
# The abstract Belief type implements initialization (initial_belief and create_belief)
# and update (belief) methods for POMDP beliefs.
# For discrete problems, it can be usually be represented as a vector.
# For tools supportng belief updates see POMDPToolbox.jl

abstract Belief <: AbstractDistribution
abstract BeliefUpdater

# returns an example initial belief for the pomdp
@pomdp_func initial_belief(pomdp::POMDP, belief::Belief = create_belief(pomdp))

# allocates and returns an empty problem-native belief structure
@pomdp_func create_belief(pomdp::POMDP)

# creates a belief object of the type used by `updater` (for preallocation purposes)
@pomdp_func create_belief(updater::BeliefUpdater)

# updates the belief given the old belief (belief_old), the action and the observation
@pomdp_func update(updater::BeliefUpdater, belief_old::Belief, action::Action, obs::Observation, belief_new::Belief=create_belief(updater))

# returns a belief that can be updated using `updater` that has a similar distribution to `b` (this conversion may be lossy)
@pomdp_func convert_belief(updater::BeliefUpdater, belief::Belief, new_belief::Belief=create_belief(updater)) = belief

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module KernSmooth

    export locpoly, dpill

    const libkernsmooth = joinpath(Pkg.dir("KernSmooth"), "deps", "libkernsmooth.so")

    include("utils.jl")
    include("locpoly.jl")
    include("dpill.jl")

end

ext: jl


text:
# ---
# title: 853. Car Fleet
# id: problem853
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort
# link: <https://leetcode.com/problems/car-fleet/description/>
# hidden: true
# ---
#
# `N` cars are going to the same destination along a one lane road.  The
# destination is `target` miles away.
#
# Each car `i` has a constant speed `speed[i]` (in miles per hour), and initial
# position `position[i]` miles towards the target along the road.
#
# A car can never pass another car ahead of it, but it can catch up to it, and
# drive bumper to bumper at the same speed.
#
# The distance between these two cars is ignored - they are assumed to have the
# same position.
#
# A _car fleet_ is some non-empty set of cars driving at the same position and
# same speed.  Note that a single car is also a car fleet.
#
# If a car catches up to a car fleet right at the destination point, it will
# still be considered as one car fleet.
#
#
# How many car fleets will arrive at the destination?
#
#
#
# **Example 1:**
#
#
#
#     Input: target = 12, position = [10,8,0,5,3], speed = [2,4,1,1,3]
#     Output: 3
#     **Explanation** :
#     The cars starting at 10 and 8 become a fleet, meeting each other at 12.
#     The car starting at 0 doesn't catch up to any other car, so it is a fleet by itself.
#     The cars starting at 5 and 3 become a fleet, meeting each other at 6.
#     Note that no other cars meet these fleets before the destination, so the answer is 3.
#
#
#
# **Note:**
#
#   1. `0 <= N <= 10 ^ 4`
#   2. `0 < target <= 10 ^ 6`
#   3. `0 < speed[i] <= 10 ^ 6`
#   4. `0 <= position[i] < target`
#   5. All initial positions are different.
#
#
## @lc code=start
using LeetCode

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


text:
# # Classification of MNIST dataset using ConvNet

# In this tutorial, we build a convolutional neural network (ConvNet or CNN) known as [LeNet5](https://en.wikipedia.org/wiki/LeNet)
# to classify [MNIST](http://yann.lecun.com/exdb/mnist/) handwritten digits.

# LeNet5 is one of the earliest CNNs. It was originally used for recognizing handwritten characters. At a high level LeNet (LeNet-5) consists of two parts:

# * A convolutional encoder consisting of two convolutional layers.
# * A dense block consisting of three fully-connected layers.

# The basic units in each convolutional block are a convolutional layer, a sigmoid activation function,
# and a subsequent average pooling operation. Each convolutional layer uses a 5×5 kernel and a sigmoid activation function.
# These layers map spatially arranged inputs to a number of two-dimensional feature maps, typically increasing the number of channels.
# The first convolutional layer has 6 output channels, while the second has 16.
# Each 2×2 pooling operation (stride 2) reduces dimensionality by a factor of 4 via spatial downsampling.
# The convolutional block emits an output with shape given by (width, height, number of channels,  batch size).

# ![LeNet-5](../conv_mnist/docs/LeNet-5.png)

# Source: https://d2l.ai/chapter_convolutional-neural-networks/lenet.html

# >**Note:** The original architecture of Lenet5 used the sigmoind activation function. However, this is a a modernized version since it uses the RELU activation function instead.

# If you need more information about how CNNs work and related technical concepts, check out the following resources:

# * [Gradient-Based Learning Applied to Document Recognition](http://yann.lecun.com/exdb/publis/pdf/lecun-01a.pdf) . This is LeNet5 original paper by Yann LeCunn and others.
# * [Convolutional Neural Networks for Visual Recognition](https://cs231n.github.io/convolutional-networks/).
# * [Neural Networks in Flux.jl with Huda Nassar (working with the MNIST dataset)](https://youtu.be/Oxi0Pfmskus).
# * [Dive into Deep Learning", 2020](https://d2l.ai/chapter_convolutional-neural-networks/lenet.html).

# This example demonstrates Flux’s Convolution and pooling layers, the usage of TensorBoardLogger,
# how to write out the saved model to the file `mnist_conv.bson`,
# and also combines various packages from the Julia ecosystem with Flux.

# To run this example, we need the following packages:

using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using Flux: onehotbatch, onecold, flatten
using Flux.Losses: logitcrossentropy
using Statistics, Random
using Logging: with_logger
using TensorBoardLogger: TBLogger, tb_overwrite, set_step!, set_step_increment!
using ProgressMeter: @showprogress
import MLDatasets
import BSON
using CUDA

# We set default values for the arguments for the function `train`:

Base.@kwdef mutable struct Args
    η = 3e-4             ## learning rate
    λ = 0                ## L2 regularizer param, implemented as weight decay
    batchsize = 128      ## batch size
    epochs = 10          ## number of epochs
    seed = 0             ## set seed > 0 for reproducibility
    use_cuda = true      ## if true use cuda (if available)
    infotime = 1 	     ## report every `infotime` epochs
    checktime = 5        ## Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    tblogger = true      ## log training with tensorboard
    savepath = "runs/"   ## results path
end

# ## Data

# We create the function `get_data` to load the MNIST train and test data from [MLDatasets](https://github.com/JuliaML/MLDatasets.jl) and reshape them so that they are in the shape that Flux expects.

function get_data(args)
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    xtrain = reshape(xtrain, 28, 28, 1, :)
    xtest = reshape(xtest, 28, 28, 1, :)

    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    train_loader = DataLoader((xtrain, ytrain), batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader((xtest, ytest),  batchsize=args.batchsize)

    return train_loader, test_loader
end

# The function `get_data` performs the following tasks:

# * **Loads MNIST dataset:** Loads the train and test set tensors. The shape of the train data is `28x28x60000` and the test data is `28x28x10000`.
# * **Reshapes the train and test data:**  Notice that we reshape the data so that we can pass it as arguments for the input layer of the model.
# * **One-hot encodes the train and test labels:** Creates a batch of one-hot vectors so we can pass the labels of the data as arguments for the loss function. For this example, we use the [logitcrossentropy](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.logitcrossentropy) function and it expects data to be one-hot encoded.
# * **Creates mini-batches of data:** Creates two DataLoader objects (train and test) that handle data mini-batches of size `128 ` (as defined above). We create these two objects so that we can pass the entire data set through the loss function at once when training our model. Also, it shuffles the data points during each iteration (`shuffle=true`).

# ## Model

# We create the LeNet5 "constructor". It uses Flux's built-in [Convolutional and pooling layers](https://fluxml.ai/Flux.jl/stable/models/layers/#Convolution-and-Pooling-Layers):

function LeNet5(; imgsize=(28,28,1), nclasses=10)
    out_conv_size = (imgsize[1]÷4 - 3, imgsize[2]÷4 - 3, 16)

    return Chain(
            Conv((5, 5), imgsize[end]=>6, relu),
            MaxPool((2, 2)),
            Conv((5, 5), 6=>16, relu),
            MaxPool((2, 2)),
            flatten,
            Dense(prod(out_conv_size), 120, relu),
            Dense(120, 84, relu),
            Dense(84, nclasses)
          )
end

# **Note:** The model can be adapted to any image size and any number of output classes.

function get_data(args)
    xtrain, ytrain = MLDatasets.MNIST(:train)[:]
    xtest, ytest = MLDatasets.MNIST(:test)[:]

# ## Loss function

# We use the function [logitcrossentropy](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.logitcrossentropy) to compute the difference between
# the predicted and actual values (loss).

loss(ŷ, y) = logitcrossentropy(ŷ, y)

# Also, we create the function `eval_loss_accuracy` to output the loss and the accuracy during training:

function eval_loss_accuracy(loader, model, device)
    l = 0f0
    acc = 0
    ntot = 0
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        l += loss(ŷ, y) * size(x)[end]
        acc += sum(onecold(ŷ |> cpu) .== onecold(y |> cpu))
        ntot += size(x)[end]
    end
    return (loss = l/ntot |> round4, acc = acc/ntot*100 |> round4)
end

# ## Utility functions
# We need a couple of functions to obtain the total number of the model's parameters. Also, we create a function to round numbers to four digits.

num_params(model) = sum(length, Flux.params(model))
round4(x) = round(x, digits=4)

# ## Train the model

# Finally, we define the function `train` that calls the functions defined above to train the model.

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.use_cuda && CUDA.functional()

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader = get_data(args)
    @info "Dataset MNIST: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = LeNet5() |> device
    @info "LeNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0 ## add weight decay, equivalent to L2 regularization
        opt = Optimiser(WeightDecay(args.λ), opt)
    end

    ## LOGGING UTILITIES
    if args.tblogger
        tblogger = TBLogger(args.savepath, tb_overwrite)
        set_step_increment!(tblogger, 0) ## 0 auto increment since we manually set_step!
        @info "TensorBoard logging at \"$(args.savepath)\""
    end

    function report(epoch)
        train = eval_loss_accuracy(train_loader, model, device)
        test = eval_loss_accuracy(test_loader, model, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
        if args.tblogger
            set_step!(tblogger, epoch)
            with_logger(tblogger) do
                @info "train" loss=train.loss  acc=train.acc
                @info "test"  loss=test.loss   acc=test.acc
            end
        end
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        @showprogress for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                    ŷ = model(x)
                    loss(ŷ, y)
                end

            Flux.Optimise.update!(opt, ps, gs)
        end

        ## Printing and logging
        epoch % args.infotime == 0 && report(epoch)
        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model.bson")
            let model = cpu(model) ## return model to cpu before serialization
                BSON.@save modelpath model epoch
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
end

# The function `train` performs the following tasks:

# * Checks whether there is a GPU available and uses it for training the model. Otherwise, it uses the CPU.
# * Loads the MNIST data using the function `get_data`.
# * Creates the model and uses the [ADAM optimiser](https://fluxml.ai/Flux.jl/stable/training/optimisers/#Flux.Optimise.ADAM) with weight decay.
# * Loads the [TensorBoardLogger.jl](https://github.com/JuliaLogging/TensorBoardLogger.jl) for logging data to Tensorboard.
# * Creates the function `report` for computing the loss and accuracy during the training loop. It outputs these values to the TensorBoardLogger.
# * Runs the training loop using [Flux’s training routine](https://fluxml.ai/Flux.jl/stable/training/training/#Training). For each epoch (step), it executes the following:
#   * Computes the model’s predictions.
#   * Computes the loss.
#   * Updates the model’s parameters.
#   * Saves the model `model.bson` every `checktime` epochs (defined as argument above.)

# ## Run the example

# We call the  function `train`:

if abspath(PROGRAM_FILE) == @__FILE__
    train()
end

max_stars_repo_name: garibarba/model-zoo

max_issues_repo_name: garibarba/model-zoo

max_forks_repo_name: garibarba/model-zoo


text:
export
    UAVLevelSetEvaluator,
    evaluate_parameters!

const CONFLICT_DISTANCE = 250

type UAVLevelSetEvaluator <: LevelSetEvaluator
    env::Environment
    num_encounters::Int64
    function UAVLevelSetEvaluator(num_encounters::Int64 = 500)
        env = build_environment()
        new(env, num_encounters)
    end
end

function update_parameters!(lse::UAVLevelSetEvaluator, psi_std::Float64,
        theta_std::Float64)
    lse.env.uavs[1].sensor_psi_std = psi_std
    lse.env.uavs[1].sensor_theta_std = theta_std
    return lse
end

function set_seeds!(lse::UAVLevelSetEvaluator, scenario_seed::Int64,
        simulation_seed::Int64)
    srand(lse.env.rng, scenario_seed)
    srand(lse.env.dynamics.rng, scenario_seed)
    srand(lse.env.uavs[1].rng, simulation_seed)
    srand(lse.env.uavs[2].rng, simulation_seed)
    return lse
end

function analyze(encounters::Array{Encounter})
    conflict_count = 0
    for e in encounters
        for horizontal_distance in e.min_horizontal_dists
            if horizontal_distance < CONFLICT_DISTANCE
                conflict_count += 1
                break
            end
        end
    end
    return conflict_count / length(encounters)
end

"""
Description:
    - Run simulations to evaluate the probability of an nmac for given parameters

Args:
    - lse: uavs level set evaluator
    - params: array of parameter values at which to evaluate the target value
        params[1] = standard deviation of relative heading sensor (psi_std)
            - min: 0. rad
            - max: pi / 2 rad
        params[1] = standard deviation of relative angle sensor (theta_std)
            - min: 0. rad
            - max: pi / 2 rad
    - scenario_seed: random seed from which to generate the scenario
    - simulation_seed: random seed to use for simulating the scenario

Returns:
    - value of the target (nmac probability)
        - min: 0.
        - max: 1.
"""
function evaluate_parameters!(lse::UAVLevelSetEvaluator, params::Array{Float64},
        scenario_seed::Int = 1, simulation_seed::Int = 1)
    update_parameters!(lse, params[1], params[2])
    set_seeds!(lse, scenario_seed, simulation_seed)
    encounters = simulate_encounters(lse.env, lse.num_encounters)
    p_nmac = analyze(encounters)
    return p_nmac
end


text:
struct BetaArray <: Distribution{Array{Float64}} end
const beta_array = BetaArray()

function Gen.logpdf(::BetaArray, x::Array{Float64}, a, b)
    pdf((val, a_i, b_i)) = logpdf(beta,val,a_i,b_i)
    logprobs = pdf.(zip(x,a,b))
    sum(logprobs)
end

function Gen.random(::BetaArray, a, b)
    [random(beta,a[i],b[i]) for i in 1:length(a)]
end

function Gen.logpdf_grad(::BetaArray, x::Array{Float64}, a, b)
    pdf_grad((val, a_i, b_i)) = logpdf_grad(beta,val,a_i,b_i)
    grads = pdf_grad.(zip(x,a,b))
    map(i->i[1], grads),map(i->i[2], grads),map(i->i[3], grads)
end

export beta_array


text:
using DoubleFloats, GenericLinearAlgebra, LinearAlgebra
using Test

using Random
phi = Base.MathConstants.golden

include("specialvalues.jl")
include("construct.jl")
include("compare.jl")
include("promote.jl")

include("arithmetic.jl")
include("functions.jl")
include("matmul.jl")
include("complex.jl")

include("concrete_accuracy.jl")
include("function_accuracy.jl")

include("coverage.jl")

include("string.jl")
include("parse.jl")


text:
"""

# link

Compute the link function

### Method
```julia
link(xrange, chain, vars, xbar)
```

### Required arguments
```julia
* `xrange::Turing model`  : Range over which link values are computed
* `chain::Float64`             : Chain samples used
* `vars::Float64`               : Variables in chain used
* `xbar::Float64`               : Mean value of observed predictor
```

### Return values
```julia
* `result`                       : Vector of link values
```

"""
function link(xrange, chain, vars, xbar)
  res = [chain.value[:, vars[1], :] + chain.value[:, vars[2], :] * (x - xbar) for x in xrange]
end


text:
#=
MNIST sample
Taken from model zoo of Julia and modified.
https://github.com/FluxML/model-zoo/blob/master/vision/mnist/mlp.jl
=#

using Printf

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, @epochs
using Base.Iterators: partition
using BSON: @load, @save
using CuArrays
using Random
using Images:channelview

function prepare_dataset(;train=true)
    train_or_test = ifelse(train,:train,:test)
    imgs = MNIST.images(train_or_test)
    labels = MNIST.labels(train_or_test)
    Y = onehotbatch(labels,0:9)
    return imgs, Y
end

function define_model(;hidden::Int)
    mlp = Chain(Conv((3,3), 1=>hidden, relu, pad=1, stride=1),
                MaxPool((2,2)),
                x -> reshape(x, :, size(x, 4)),
                Dense(14*14*8,10),
                softmax)
    return mlp
end

function split_dataset_random(X, Y)
    divide_ratio=0.9
    shuffled_indices = shuffle(1:size(Y)[2])
    divide_idx = round(Int,0.9*length(shuffled_indices))
    train_indices = shuffled_indices[1:divide_idx]
    val_indices = shuffled_indices[divide_idx:end]
    train_X = X[train_indices]
    train_Y = Y[:,train_indices]
    val_X = X[val_indices]
    val_Y = Y[:,val_indices]
    return train_X, train_Y, val_X, val_Y
end

function train()
    println("Start to train")
    epochs = 3
    X, Y = prepare_dataset(train=true)
    train_img, train_Y, val_img, val_Y = split_dataset_random(X, Y)
    model = define_model(hidden=8) |> gpu
    loss(x,y)= crossentropy(model(x),y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    batchsize = 64
    train_X = [cat(permutedims(float.(train_img[batch]))...,dims=4) for batch in partition(1:length(train_img),batchsize)]
    val_X   = [cat(permutedims(float.(val_img[batch]))...,dims=4) for batch in partition(1:length(val_img),batchsize)]
    train_dataset = gpu.([(train_X[i] ,train_Y[:,batch]) for (i, batch) in enumerate(partition(1:size(train_Y)[2],batchsize))])
    val_dataset = gpu.([(val_X[i] ,val_Y[:,batch]) for (i,batch) in enumerate(partition(1:size(val_Y)[2],batchsize))])
    callback_count = 0
    eval_callback = function callback()
        callback_count += 1
        if callback_count == length(train_dataset)
            println("action for each epoch")
            total_loss = 0
            total_acc = 0
            for (vx, vy) in val_dataset
                total_loss += loss(vx, vy)
                total_acc += accuracy(vx, vy)
            end
            total_loss /= length(val_dataset)
            total_acc /= length(val_dataset)
            @show total_loss, total_acc
            callback_count = 0
            pretrained = model |> cpu
            @save "pretrained.bson" pretrained
            callback_count = 0
        end
        if callback_count % 50 == 0
            progress = callback_count / length(train_dataset)
           @printf("%.3f\n", progress)
        end
    end
    optimizer = ADAM(params(model))

    @epochs epochs Flux.train!(loss, train_dataset, optimizer, cb = eval_callback)

    pretrained = model |> cpu
    weights = Tracker.data.(params(pretrained))
    @save "pretrained.bson" pretrained
    @save "weights.bson" weights
    println("Finished to train")
end

function predict()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "pretrained.bson" pretrained
    model = pretrained |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("prepare dataset")
    test_img, Y = prepare_dataset(train=false)
    batchsize=64
    test_X = [cat(permutedims(float.(test_img[batch]))...,dims=4) for batch in partition(1:length(test_img),batchsize)]
    test_label = gpu.([Y[:,batch] for (i,batch) in enumerate(partition(1:size(Y)[2],batchsize))])
    total_acc=0
    cnt = 0
    for (x,y) in zip(test_X,test_label)
        total_acc += accuracy(x|> gpu, y|>gpu)
        cnt += 1
    end
    println("acc = $(total_acc/cnt)")
    println("Done")
end

function predict2()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "weights.bson" weights
    model = define_model(hidden=8)
    Flux.loadparams!(model, weights)
    model = model |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    test_img, Y = prepare_dataset(train=false)
    batchsize=64
    test_X = [cat(permutedims(float.(test_img[batch]))...,dims=4) for batch in partition(1:length(test_img),batchsize)]
    test_label = gpu.([Y[:,batch] for (i,batch) in enumerate(partition(1:size(Y)[2],batchsize))])
    total_acc=0
    cnt = 0
    for (x,y) in zip(test_X,test_label)
        total_acc += accuracy(x|> gpu, y|>gpu)
        cnt += 1
    end
    println("acc = $(total_acc/cnt)")
    println("Done")
end

function main()
    train()
    predict()
    predict2()
end

#main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module WebFuncs

export Mapping, expose!, run, default_port

import HTTP
import JSON
import Unmarshal
import UUIDs
using Sockets: localhost

const default_port = 3000

struct Lambda
    func::Function
    Input::DataType
end

const Mapping = Dict{UUIDs.UUID,Lambda}

function expose!(m::Mapping, func::Function, input_type::DataType=Dict{AbstractString,Any})
    key = UUIDs.uuid4()
    m[key] = Lambda(func, input_type)
    key
end

function expose!(m::Mapping, funcs::Vector{<:Function}, input_types::Vector{DataType})
    func_keys = [UUIDs.uuid4() for _ in 1:length(funcs)]
    for (k, f, t) in zip(func_keys, funcs, input_types)
        m[k] = Lambda(f, t)
    end
    func_keys
end

function expose!(m::Mapping, funcs::Vector{<:Function}, input_type::DataType=Dict{AbstractString,Any})
    type_vec = [input_type for _ in 1:length(funcs)]
    expose!(m, funcs, type_vec)
end

function parse_input(data::Vector{UInt8},DT::DataType=Dict{AbstractString,Any})
    parsed_dict = JSON.Parser.parse(join([Char(v) for v in data]))
	if DT <: Dict
		parsed_dict
	else
		Unmarshal.unmarshal(DT, parsed_dict)
	end
end

function run(m::Mapping, port::Int = default_port)
    # dispatches the request with parsed req body to corresponding Lambda
    r = HTTP.Router()
    for (k,f) in m
        HTTP.register!(r,string("/",k),_func_handler(f))
    end
    HTTP.register!(r, string(localhost), req::HTTP.Request -> HTTP.Response(200, "Hello WebFuncs"))
    server = HTTP.Server(r)
    HTTP.serve(server, localhost, port, verbose = true)
end

function _func_handler(lambda::Lambda)
    function(req::HTTP.Request)
        input = parse_input(req.body, lambda.Input)
        result = Dict(["result" => lambda.func(input)])
        HTTP.Response(200, JSON.json(result))
    end
end

end


text:
using AES, Test

k = [0x2b, 0x7e, 0x15, 0x16, 0x28, 0xae, 0xd2, 0xa6, 0xab, 0xf7, 0x15, 0x88, 0x09, 0xcf, 0x4f, 0x3c]
key = AES128Key(k)

plaintext = "The quick brown fox jumps over the lazy dog."
plaintext8 = "The quick brown fox jumps over the 😴 dog."
iv = [0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f]

cipher = AESCipher(;key_length=128, mode=AES.CTR, key=key)
ct = encrypt(plaintext, cipher; iv=iv)
pt = decrypt(ct, cipher)

cipher = AESCipher(;key_length=128, mode=AES.CTR, key=key)
ct8 = encrypt(plaintext8, cipher; iv=iv)
pt8 = decrypt(ct8, cipher)

ciphertext = [0x04, 0x96, 0x02, 0xec, 0xe8, 0x18, 0x5b, 0xd5,
              0xb1, 0x29, 0x55, 0x9b, 0xf4, 0xd8, 0x82, 0x40,
              0xc9, 0x4e, 0xf6, 0x6e, 0xed, 0x96, 0xfe, 0xda,
              0xa9, 0x11, 0x24, 0x92, 0x1e, 0x85, 0x83, 0x2b,
              0xfe, 0x8f, 0xaf, 0xe5, 0x64, 0xf9, 0xe9, 0x62,
              0xdd, 0xa1, 0xe7, 0xb6]

@testset "ASCII" begin
    @test pt == transcode(UInt8,plaintext)
    @test ct.data == ciphertext
end

@testset "UTF-8" begin
    @test pt8 == transcode(UInt8,plaintext8)
end


text:
# functions to deal with MATLAB arrays

const libmx = open_matlab_library("libmx")
mxfunc(s::Symbol) = Libdl.dlsym(libmx, s)

type MxArray
    ptr::Ptr{Void}
    MxArray(p::Ptr{Void}) = new(p)
end

# delete & duplicate

function delete(mx::MxArray)
    if !(mx.ptr == C_NULL)
        ccall(mxfunc(:mxDestroyArray), Void, (Ptr{Void},), mx.ptr)
    end
    mx.ptr = C_NULL
end

function copy(mx::MxArray)
    pm::Ptr{Void} = ccall(mxfunc(:mxDuplicateArray), Ptr{Void}, (Ptr{Void},), mx.ptr)
    MxArray(pm)
end

# functions to create mxArray from Julia values/arrays

MxRealNum = Union{Float64,Float32,Int32,UInt32,Int64,UInt64,Int16,UInt16,Int8,UInt8,Bool}
MxComplexNum = Union{Complex64, Complex128}
MxNum = Union{MxRealNum, MxComplexNum}

###########################################################
#
#  MATLAB types
#
###########################################################

typealias mwSize UInt
typealias mwIndex Int
typealias mxChar UInt16
typealias mxClassID Cint
typealias mxComplexity Cint

const mxUNKNOWN_CLASS  = convert(mxClassID, 0)
const mxCELL_CLASS     = convert(mxClassID, 1)
const mxSTRUCT_CLASS   = convert(mxClassID, 2)
const mxLOGICAL_CLASS  = convert(mxClassID, 3)
const mxCHAR_CLASS     = convert(mxClassID, 4)
const mxVOID_CLASS     = convert(mxClassID, 5)
const mxDOUBLE_CLASS   = convert(mxClassID, 6)
const mxSINGLE_CLASS   = convert(mxClassID, 7)
const mxINT8_CLASS     = convert(mxClassID, 8)
const mxUINT8_CLASS    = convert(mxClassID, 9)
const mxINT16_CLASS    = convert(mxClassID, 10)
const mxUINT16_CLASS   = convert(mxClassID, 11)
const mxINT32_CLASS    = convert(mxClassID, 12)
const mxUINT32_CLASS   = convert(mxClassID, 13)
const mxINT64_CLASS    = convert(mxClassID, 14)
const mxUINT64_CLASS   = convert(mxClassID, 15)
const mxFUNCTION_CLASS = convert(mxClassID, 16)
const mxOPAQUE_CLASS   = convert(mxClassID, 17)
const mxOBJECT_CLASS   = convert(mxClassID, 18)

const mxREAL    = convert(mxComplexity, 0)
const mxCOMPLEX = convert(mxComplexity, 1)

mxclassid(::Type{Bool})    = mxLOGICAL_CLASS::Cint
mxclassid(::Union{Type{Float64}, Type{Complex128}}) = mxDOUBLE_CLASS::Cint
mxclassid(::Union{Type{Float32}, Type{Complex64}}) = mxSINGLE_CLASS::Cint
mxclassid(::Type{Int8})    = mxINT8_CLASS::Cint
mxclassid(::Type{UInt8})   = mxUINT8_CLASS::Cint
mxclassid(::Type{Int16})   = mxINT16_CLASS::Cint
mxclassid(::Type{UInt16})  = mxUINT16_CLASS::Cint
mxclassid(::Type{Int32})   = mxINT32_CLASS::Cint
mxclassid(::Type{UInt32})  = mxUINT32_CLASS::Cint
mxclassid(::Type{Int64})   = mxINT64_CLASS::Cint
mxclassid(::Type{UInt64})  = mxUINT64_CLASS::Cint

mxcomplexflag{T<:MxRealNum}(::Type{T})    = mxREAL
mxcomplexflag{T<:MxComplexNum}(::Type{T}) = mxCOMPLEX

const classid_type_map = Dict{mxClassID,Type}(
    mxLOGICAL_CLASS => Bool,
    mxCHAR_CLASS    => mxChar,
    mxDOUBLE_CLASS  => Float64,
    mxSINGLE_CLASS  => Float32,
    mxINT8_CLASS    => Int8,
    mxUINT8_CLASS   => UInt8,
    mxINT16_CLASS   => Int16,
    mxUINT16_CLASS  => UInt16,
    mxINT32_CLASS   => Int32,
    mxUINT32_CLASS  => UInt32,
    mxINT64_CLASS   => Int64,
    mxUINT64_CLASS  => UInt64
)

function mxclassid_to_type(cid::mxClassID)
    ty = get(classid_type_map::Dict{mxClassID, Type}, cid, nothing)
    if ty == nothing
        throw(ArgumentError("The input class id is not a primitive type id."))
    end
    ty
end

###########################################################
#
#  Functions to access mxArray
#
#  Part of the functions (e.g. mxGetNumberOfDimensions)
#  are actually a macro replacement of an internal
#  function name as (xxxx_730)
#
###########################################################

# pre-cached some useful functions

const _mx_free = mxfunc(:mxFree)

const _mx_get_classid = mxfunc(:mxGetClassID)
const _mx_get_m = mxfunc(:mxGetM)
const _mx_get_n = mxfunc(:mxGetN)
const _mx_get_nelems = mxfunc(:mxGetNumberOfElements)
const _mx_get_ndims  = mxfunc(:mxGetNumberOfDimensions_730)
const _mx_get_elemsize = mxfunc(:mxGetElementSize)
const _mx_get_data = mxfunc(:mxGetData)
const _mx_get_dims = mxfunc(:mxGetDimensions_730)
const _mx_get_nfields = mxfunc(:mxGetNumberOfFields)
const _mx_get_pr = mxfunc(:mxGetPr)
const _mx_get_pi = mxfunc(:mxGetPi)
const _mx_get_ir = mxfunc(:mxGetIr_730)
const _mx_get_jc = mxfunc(:mxGetJc_730)

const _mx_is_double = mxfunc(:mxIsDouble)
const _mx_is_single = mxfunc(:mxIsSingle)
const _mx_is_int64  = mxfunc(:mxIsInt64)
const _mx_is_uint64 = mxfunc(:mxIsUint64)
const _mx_is_int32  = mxfunc(:mxIsInt32)
const _mx_is_uint32 = mxfunc(:mxIsUint32)
const _mx_is_int16  = mxfunc(:mxIsInt16)
const _mx_is_uint16 = mxfunc(:mxIsUint16)
const _mx_is_int8   = mxfunc(:mxIsInt8)
const _mx_is_uint8  = mxfunc(:mxIsUint8)
const _mx_is_char   = mxfunc(:mxIsChar)

const _mx_is_numeric = mxfunc(:mxIsNumeric)
const _mx_is_logical = mxfunc(:mxIsLogical)
const _mx_is_complex = mxfunc(:mxIsComplex)
const _mx_is_sparse  = mxfunc(:mxIsSparse)
const _mx_is_empty   = mxfunc(:mxIsEmpty)
const _mx_is_struct  = mxfunc(:mxIsStruct)
const _mx_is_cell    = mxfunc(:mxIsCell)

# getting simple attributes

macro mx_get(name, fun, ret, cnv)
    :($(name)(mx::MxArray) = convert($(cnv), ccall($(fun)::Ptr{Void}, $(ret), (Ptr{Void},), mx.ptr)))
end
@mx_get(classid, _mx_get_classid, mxClassID, mxClassID)
@mx_get(nrows, _mx_get_m, UInt, Int)
@mx_get(ncols, _mx_get_n, UInt, Int)
@mx_get(nelems, _mx_get_nelems, UInt, Int)
@mx_get(ndims, _mx_get_ndims, mwSize, Int)
@mx_get(elsize, _mx_get_elemsize, UInt, Int)
@mx_get(nfields, _mx_get_nfields, Cint, Int)

eltype(mx::MxArray)  = mxclassid_to_type(classid(mx))
macro mx_get_ptr(name, fun)
    :($(name)(mx::MxArray) = convert(Ptr{eltype(mx)}, ccall($(fun)::Ptr{Void}, Ptr{Void}, (Ptr{Void},), mx.ptr)))
end
@mx_get_ptr(data_ptr, _mx_get_data)
@mx_get_ptr(real_ptr, _mx_get_pr)
@mx_get_ptr(imag_ptr, _mx_get_pi)

# validation functions

macro mx_pred(name, fun)
    :($(name)(mx::MxArray) = ccall($(fun)::Ptr{Void}, Bool, (Ptr{Void},), mx.ptr))
end
@mx_pred(is_double, _mx_is_double)
@mx_pred(is_single, _mx_is_single)
@mx_pred(is_int64, _mx_is_int64)
@mx_pred(is_uint64, _mx_is_uint64)
@mx_pred(is_int32, _mx_is_int32)
@mx_pred(is_uint32, _mx_is_uint32)
@mx_pred(is_int16, _mx_is_int16)
@mx_pred(is_uint16, _mx_is_uint16)
@mx_pred(is_int8, _mx_is_int8)
@mx_pred(is_uint8, _mx_is_uint8)
@mx_pred(is_numeric, _mx_is_numeric)
@mx_pred(is_logical, _mx_is_logical)
@mx_pred(is_complex, _mx_is_complex)
@mx_pred(is_sparse, _mx_is_sparse)
@mx_pred(is_struct, _mx_is_struct)
@mx_pred(is_cell, _mx_is_cell)
@mx_pred(is_char, _mx_is_char)
@mx_pred(is_empty, _mx_is_empty)

# size function

function size(mx::MxArray)
    nd = ndims(mx)
    pdims = ccall(_mx_get_dims, Ptr{mwSize}, (Ptr{Void},), mx.ptr)
    _dims = unsafe_wrap(Array, pdims, (nd,))
    dims = Array(Int, nd)
    for i = 1 : nd
        dims[i] = convert(Int, _dims[i])
    end
    tuple(dims...)
end

function size(mx::MxArray, d::Integer)
    nd = ndims(mx)
    if d <= 0
        throw(ArgumentError("The dimension must be a positive integer."))
    end

    if nd == 2
        d == 1 ? nrows(mx) :
        d == 2 ? ncols(mx) : 1
    else
        pdims = ccall(_mx_get_dims, Ptr{mwSize}, (Ptr{Void},), mx.ptr)
        _dims = unsafe_wrap(Array, pdims, (nd,))
        d <= nd ? convert(Int, _dims[d]) : 1
    end
end

###########################################################
#
#  functions to create & delete MATLAB arrays
#
###########################################################

# pre-cached functions

const _mx_create_numeric_mat = mxfunc(:mxCreateNumericMatrix_730)
const _mx_create_numeric_arr = mxfunc(:mxCreateNumericArray_730)

const _mx_create_double_scalar = mxfunc(:mxCreateDoubleScalar)
const _mx_create_logical_scalar = mxfunc(:mxCreateLogicalScalar)

const _mx_create_sparse = mxfunc(:mxCreateSparse_730)
const _mx_create_sparse_logical = mxfunc(:mxCreateSparseLogicalMatrix_730)

const _mx_create_char_array = mxfunc(:mxCreateCharArray_730)

const _mx_create_cell_array = mxfunc(:mxCreateCellArray_730)

const _mx_create_struct_matrix = mxfunc(:mxCreateStructMatrix_730)
const _mx_create_struct_array = mxfunc(:mxCreateStructArray_730)

const _mx_get_cell = mxfunc(:mxGetCell_730)
const _mx_set_cell = mxfunc(:mxSetCell_730)

const _mx_get_field = mxfunc(:mxGetField_730)
const _mx_set_field = mxfunc(:mxSetField_730)
const _mx_get_field_bynum = mxfunc(:mxGetFieldByNumber_730)
const _mx_get_fieldname = mxfunc(:mxGetFieldNameByNumber)

# create zero arrays

mxempty() = MxArray(Float64, 0, 0)

function _dims_to_mwSize(dims::Tuple{Vararg{Int}})
    ndim = length(dims)
    _dims = Array(mwSize, ndim)
    for i = 1 : ndim
        _dims[i] = convert(mwSize, dims[i])
    end
    _dims
end

function MxArray{T<:MxNum}(ty::Type{T}, dims::Tuple{Vararg{Int}})
    pm = ccall(_mx_create_numeric_arr, Ptr{Void},
        (mwSize, Ptr{mwSize}, mxClassID, mxComplexity),
        length(dims), _dims_to_mwSize(dims), mxclassid(ty), mxcomplexflag(ty))

    MxArray(pm)
end
MxArray{T<:MxNum}(ty::Type{T}, dims::Int...) = MxArray(ty, dims)

# create scalars

function MxArray(x::Float64)
    pm = ccall(_mx_create_double_scalar, Ptr{Void}, (Cdouble,), x)
    MxArray(pm)
end

function MxArray(x::Bool)
    pm = ccall(_mx_create_logical_scalar, Ptr{Void}, (Bool,), x)
    MxArray(pm)
end

function MxArray{T<:MxRealNum}(x::T)
    pm = ccall(_mx_create_numeric_mat, Ptr{Void},
        (mwSize, mwSize, mxClassID, mxComplexity),
        1, 1, mxclassid(T), mxcomplexflag(T))

    pdat = ccall(_mx_get_data, Ptr{T}, (Ptr{Void},), pm)

    unsafe_wrap(Array, pdat, (1,))[1] = x
    MxArray(pm)
end
MxArray{T<:MxComplexNum}(x::T) = MxArray([x])

# conversion from Julia variables to MATLAB
# Note: the conversion is deep-copy, as there is no way to let
# mxArray use Julia array's memory

function MxArray{T<:MxRealNum}(a::Array{T})
    mx = MxArray(T, size(a))
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
        data_ptr(mx), a, length(a) * sizeof(T))
    mx
end

function MxArray{T<:MxComplexNum}(a::Array{T})
    mx = MxArray(T, size(a))
    na = length(a)
    rdat = unsafe_wrap(Array, real_ptr(mx), na)
    idat = unsafe_wrap(Array, imag_ptr(mx), na)
    for i = 1:na
        rdat[i] = real(a[i])
        idat[i] = imag(a[i])
    end
    mx
end

MxArray(a::BitArray) = MxArray(convert(Array{Bool}, a))
MxArray(a::Range) = MxArray([a;])

# sparse matrix

function mxsparse(ty::Type{Float64}, m::Integer, n::Integer, nzmax::Integer)
    pm = ccall(_mx_create_sparse, Ptr{Void},
        (mwSize, mwSize, mwSize, mxComplexity), m, n, nzmax, mxREAL)
    MxArray(pm)
end

function mxsparse(ty::Type{Bool}, m::Integer, n::Integer, nzmax::Integer)
    pm = ccall(_mx_create_sparse_logical, Ptr{Void},
        (mwSize, mwSize, mwSize), m, n, nzmax)
    MxArray(pm)
end

function _copy_sparse_mat{V,I}(a::SparseMatrixCSC{V,I},
    ir_p::Ptr{mwIndex}, jc_p::Ptr{mwIndex}, pr_p::Ptr{V})

    colptr::Vector{I} = a.colptr
    rinds::Vector{I} = a.rowval
    v::Vector{V} = a.nzval
    n::Int = a.n
    nnz::Int = length(v)

    # Note: ir and jc contain zero-based indices

    ir = unsafe_wrap(Array, ir_p, (nnz,))
    for i = 1 : nnz
        ir[i] = rinds[i] - 1
    end

    jc = unsafe_wrap(Array, jc_p, (n+1,))
    for i = 1 : n+1
        jc[i] = colptr[i] - 1
    end

    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt), pr_p, v, nnz * sizeof(V))
end

function MxArray{V<:Union{Float64,Bool},I}(a::SparseMatrixCSC{V,I})
    m::Int = a.m
    n::Int = a.n
    nnz = length(a.nzval)
    @assert nnz == a.colptr[n+1]-1

    mx = mxsparse(V, m, n, nnz)

    ir_p = ccall(_mx_get_ir, Ptr{mwIndex}, (Ptr{Void},), mx.ptr)
    jc_p = ccall(_mx_get_jc, Ptr{mwIndex}, (Ptr{Void},), mx.ptr)
    pr_p = ccall(_mx_get_pr, Ptr{V}, (Ptr{Void},), mx.ptr)

    _copy_sparse_mat(a, ir_p, jc_p, pr_p)
    mx
end

# char arrays and string

function MxArray(s::String)
    wchars = transcode(UInt16, s)
    len = length(wchars)
    pm = ccall(_mx_create_char_array, Ptr{Void}, (mwSize, Ptr{mwSize}), 2, [1, len])
    mx = MxArray(pm)
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
        data_ptr(mx), wchars, len * sizeof(mxChar))
    mx
end

# cell arrays

function mxcellarray(dims::Tuple{Vararg{Int}})
    pm = ccall(_mx_create_cell_array, Ptr{Void}, (mwSize, Ptr{mwSize}),
        length(dims), _dims_to_mwSize(dims))
    MxArray(pm)
end
mxcellarray(dims::Int...) = mxcellarray(dims)

function get_cell(mx::MxArray, i::Integer)
    pm = ccall(_mx_get_cell, Ptr{Void}, (Ptr{Void}, mwIndex), mx.ptr, i-1)
    MxArray(pm)
end

function set_cell(mx::MxArray, i::Integer, v::MxArray)
    ccall(_mx_set_cell, Void, (Ptr{Void}, mwIndex, Ptr{Void}),
        mx.ptr, i - 1, v.ptr)
end

function mxcellarray(a::Array)
    pm = mxcellarray(size(a))
    for i = 1 : length(a)
        set_cell(pm, i, MxArray(a[i]))
    end
    pm
end

MxArray(a::Array) = mxcellarray(a)

# struct arrays

function _fieldname_array(fieldnames::String...)
    n = length(fieldnames)
    a = Array(Ptr{UInt8}, n)
    for i = 1 : n
        a[i] = Base.unsafe_convert(Ptr{UInt8}, fieldnames[i])
    end
    a
end

function mxstruct(fns::Vector{String})
    a = _fieldname_array(fns...)
    pm = ccall(_mx_create_struct_matrix, Ptr{Void},
        (mwSize, mwSize, Cint, Ptr{Ptr{UInt8}}),
        1, 1, length(a), a)
    MxArray(pm)
end

function mxstruct(fn1::String, fnr::String...)
    a = _fieldname_array(fn1, fnr...)
    pm = ccall(_mx_create_struct_matrix, Ptr{Void},
        (mwSize, mwSize, Cint, Ptr{Ptr{UInt8}}),
        1, 1, length(a), a)
    MxArray(pm)
end

function set_field(mx::MxArray, i::Integer, f::String, v::MxArray)
    ccall(_mx_set_field, Void,
        (Ptr{Void}, mwIndex, Ptr{UInt8}, Ptr{Void}),
        mx.ptr, i-1, f, v.ptr)
end

set_field(mx::MxArray, f::String, v::MxArray) = set_field(mx, 1, f, v)

function get_field(mx::MxArray, i::Integer, f::String)
    pm = ccall(_mx_get_field, Ptr{Void}, (Ptr{Void}, mwIndex, Ptr{UInt8}),
        mx.ptr, i-1, f)
    if pm == C_NULL
        throw(ArgumentError("Failed to get field."))
    end
    MxArray(pm)
end

get_field(mx::MxArray, f::String) = get_field(mx, 1, f)

function get_field(mx::MxArray, i::Integer, fn::Integer)
    pm = ccall(_mx_get_field_bynum, Ptr{Void}, (Ptr{Void}, mwIndex, Cint),
        mx.ptr, i-1, fn-1)
    if pm == C_NULL
        throw(ArgumentError("Failed to get field."))
    end
    MxArray(pm)
end

get_field(mx::MxArray, fn::Integer) = get_field(mx, 1, fn)

function get_fieldname(mx::MxArray, i::Integer)
    p = ccall(_mx_get_fieldname, Ptr{UInt8}, (Ptr{Void}, Cint),
        mx.ptr, i-1)
    unsafe_string(p)
end

if VERSION >= v"0.4.0-dev+980"
    typealias Pairs Union{Pair,NTuple{2}}
else
    typealias Pairs NTuple{2}
end

function mxstruct(pairs::Pairs...)
    nf = length(pairs)
    fieldnames = Array(String, nf)
    for i = 1 : nf
        fn = pairs[i][1]
        fieldnames[i] = string(fn)
    end
    mx = mxstruct(fieldnames)
    for i = 1 : nf
        set_field(mx, fieldnames[i], MxArray(pairs[i][2]))
    end
    mx
end

function mxstruct{T}(d::T)
    names = fieldnames(T)
    names_str = map(string, names)
    mx = mxstruct(names_str...)
    for i = 1:length(names)
        set_field(mx, names_str[i], MxArray(getfield(d, names[i])))
    end
    mx
end

function mxstructarray{T}(d::Array{T})
    names = fieldnames(T)
    names_str = map(string, names)
    a = _fieldname_array(names_str...)

    pm = ccall(_mx_create_struct_array, Ptr{Void}, (mwSize, Ptr{mwSize}, Cint,
        Ptr{Ptr{UInt8}}), ndims(d), _dims_to_mwSize(size(d)), length(a), a)
    mx = MxArray(pm)

    for i = 1:length(d), j = 1:length(names)
        set_field(mx, i, names_str[j],
            MxArray(getfield(d[i], names[j])))
    end
    mx
end

mxstruct(d::Associative) = mxstruct(collect(d)...)
MxArray(d) = mxstruct(d)

###########################################################
#
#  convert from MATLAB to Julia
#
###########################################################

# use deep-copy from MATLAB variable to Julia array
# in practice, MATLAB variable often has shorter life-cycle

function _jarrayx(fun::String, mx::MxArray, siz::Tuple)
    if is_numeric(mx) || is_logical(mx)
        @assert !is_sparse(mx)
        T = eltype(mx)
        if is_complex(mx)
            rdat = unsafe_wrap(Array, real_ptr(mx), siz)
            idat = unsafe_wrap(Array, imag_ptr(mx), siz)
            a = complex(rdat, idat)
        else
            a = Array(T, siz)
            if !isempty(a)
                ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt),
                    a, data_ptr(mx), sizeof(T) * length(a))
            end
        end
        a
        #unsafe_wrap(Array, data_ptr(mx), siz)
    elseif is_cell(mx)
        a = Array(Any, siz)
        for i = 1 : length(a)
            a[i] = Any(get_cell(mx, i))
        end
        a
    else
        throw(ArgumentError("$(fun) only applies to numeric, logical or cell arrays."))
    end
end

jarray(mx::MxArray) = _jarrayx("jarray", mx, size(mx))
jvector(mx::MxArray) = _jarrayx("jvector", mx, (nelems(mx),))

function jmatrix(mx::MxArray)
    if ndims(mx) != 2
        throw(ArgumentError("jmatrix only applies to MATLAB arrays with ndims == 2."))
    end
    _jarrayx("jmatrix", mx, (nrows(mx), ncols(mx)))
end

function jscalar(mx::MxArray)
    if !(nelems(mx) == 1 && (is_logical(mx) || is_numeric(mx)))
        throw(ArgumentError("jscalar only applies to numeric or logical arrays with exactly one element."))
    end
    @assert !is_sparse(mx)
    if is_complex(mx)
        unsafe_wrap(Array, real_ptr(mx), (1,), false)[1] + im*unsafe_wrap(Array, imag_ptr(mx), (1,))[1]
    else
        unsafe_wrap(Array, data_ptr(mx), (1,))[1]
    end
end

function _jsparse{T<:MxRealNum}(ty::Type{T}, mx::MxArray)
    m = nrows(mx)
    n = ncols(mx)
    ir_ptr = ccall(_mx_get_ir, Ptr{mwIndex}, (Ptr{Void},), mx.ptr)
    jc_ptr = ccall(_mx_get_jc, Ptr{mwIndex}, (Ptr{Void},), mx.ptr)
    pr_ptr = ccall(_mx_get_pr, Ptr{T}, (Ptr{Void},), mx.ptr)

    jc_a::Vector{mwIndex} = unsafe_wrap(Array, jc_ptr, (n+1,))
    nnz = jc_a[n+1]

    ir = Array(Int, nnz)
    jc = Array(Int, n+1)

    ir_x = unsafe_wrap(Array, ir_ptr, (nnz,))
    for i = 1 : nnz
        ir[i] = ir_x[i] + 1
    end

    jc_x = unsafe_wrap(Array, jc_ptr, (n+1,))
    for i = 1 : n+1
        jc[i] = jc_x[i] + 1
    end

    pr::Vector{T} = copy(unsafe_wrap(Array, pr_ptr, (nnz,)))
    SparseMatrixCSC(m, n, jc, ir, pr)
end

function jsparse(mx::MxArray)
    if !is_sparse(mx)
        throw(ArgumentError("jsparse only applies to sparse matrices."))
    end
    _jsparse(eltype(mx), mx)
end

function String(mx::MxArray)
    if !(classid(mx) == mxCHAR_CLASS && ((ndims(mx) == 2 && nrows(mx) == 1) || is_empty(mx)))
        throw(ArgumentError("String only applies to char row vectors."))
    end
    transcode(String, unsafe_wrap(Array, Ptr{mxChar}(data_ptr(mx)), ncols(mx), false))
end

function Dict(mx::MxArray)
    if !(is_struct(mx) && nelems(mx) == 1)
        throw(ArgumentError("jdict only applies to a single struct."))
    end
    nf = nfields(mx)
    fnames = Array(String, nf)
    fvals = Array(Any, nf)
    for i = 1 : nf
        fnames[i] = get_fieldname(mx, i)
        pv::Ptr{Void} = ccall(_mx_get_field_bynum,
            Ptr{Void}, (Ptr{Void}, mwIndex, Cint),
            mx.ptr, 0, i-1)
        fx = MxArray(pv)
        fvals[i] = Any(fx)
    end
    Dict(zip(fnames, fvals))
end

function Function(mx::MxArray)
    if classid(mx) != mxFUNCTION_CLASS || nelems(mx) != 1
        throw(ArgumentError("Function only applies to a single function handle."))
    end
    (args...) -> mx(args...)
end

function jvalue(mx::MxArray)
    if is_numeric(mx) || is_logical(mx)
        if !is_sparse(mx)
            nelems(mx) == 1 ? jscalar(mx) :
            ndims(mx) == 2 ? (ncols(mx) == 1 ? jvector(mx) : jmatrix(mx)) :
            jarray(mx)
        else
            jsparse(mx)
        end
    elseif is_char(mx) && (nrows(mx) == 1 || is_empty(mx))
        String(mx)
    elseif is_cell(mx)
        ndims(mx) == 2 ? (ncols(mx) == 1 ? jvector(mx) : jmatrix(mx)) :
        jarray(mx)
    elseif is_struct(mx) && nelems(mx) == 1
        Dict(mx)
    elseif classid(mx) == mxFUNCTION_CLASS && nelems(mx) == 1
        Function(mx)
    else
        mx
    end
end

# deep conversion from MATLAB variable to Julia array

# convert(::Type{Array}, mx::MxArray)  = jarray(mx)
# convert(::Type{Vector}, mx::MxArray) = jvector(mx)
# convert(::Type{Matrix}, mx::MxArray) = jmatrix(mx)
# convert(::Type{Number}, mx::MxArray) = jscalar(mx)::Number
# convert(::Type{String}, mx::MxArray) = String(mx)::String
# convert(::Type{Dict}, mx::MxArray) = jdict(mx)
# convert(::Type{SparseMatrixCSC}, mx::MxArray) = jsparse(mx)


text:
import Base: +, *, zero
using Distances: sqeuclidean, SqEuclidean, Euclidean
using Base.Broadcast: broadcast_shape

abstract type Kernel end

#
# Base Kernels
#

"""
    ZeroKernel <: Kernel

A rank 0 `Kernel` that always returns zero.
"""
struct ZeroKernel{T<:Real} <: Kernel end
ZeroKernel() = ZeroKernel{Float64}()
zero(::Kernel) = ZeroKernel()

# Binary methods.
ew(k::ZeroKernel{T}, x::AV, x′::AV) where {T} = zeros(T, broadcast_shape(size(x), size(x′)))
pw(k::ZeroKernel{T}, x::AV, x′::AV) where {T} = zeros(T, length(x), length(x′))

# Unary methods.
ew(k::ZeroKernel{T}, x::AV) where {T} = zeros(T, length(x))
pw(k::ZeroKernel{T}, x::AV) where {T} = zeros(T, length(x), length(x))

"""
    OneKernel{T<:Real} <: Kernel

A rank 1 constant `Kernel`. Useful for consistency when creating composite Kernels,
but (almost certainly) shouldn't be used as a base `Kernel`.
"""
struct OneKernel{T<:Real} <: Kernel end
OneKernel() = OneKernel{Float64}()

# Binary methods.
ew(k::OneKernel{T}, x::AV, x′::AV) where {T} = ones(T, broadcast_shape(size(x), size(x′)))
pw(k::OneKernel{T}, x::AV, x′::AV) where {T} = ones(T, length(x), length(x′))

# Unary methods.
ew(k::OneKernel{T}, x::AV) where {T} = ones(T, length(x))
pw(k::OneKernel{T}, x::AV) where {T} = ones(T, length(x), length(x))

"""
    ConstKernel{T} <: Kernel

A rank 1 kernel that returns the same value `c` everywhere.
"""
struct ConstKernel{T} <: Kernel
    c::T
end

# Binary methods.
ew(k::ConstKernel, x::AV, x′::AV) = fill(k.c, broadcast_shape(size(x), size(x′))...)
pw(k::ConstKernel, x::AV, x′::AV) = fill(k.c, length(x), length(x′))

# Unary methods.
ew(k::ConstKernel, x::AV) = fill(k.c, length(x))
pw(k::ConstKernel, x::AV) = pw(k, x, x)

"""
    EQ <: Kernel

The standardised Exponentiated Quadratic kernel (no free parameters).
"""
struct EQ <: Kernel end

# Binary methods.
ew(::EQ, x::AV, x′::AV) = exp.(.-ew(SqEuclidean(), x, x′) ./ 2)
pw(::EQ, x::AV, x′::AV) = exp.(.-pw(SqEuclidean(), x, x′) ./ 2)

# Unary methods.
ew(::EQ, x::AV) = exp.(.-ew(SqEuclidean(), x) ./ 2)
pw(::EQ, x::AV) = exp.(.-pw(SqEuclidean(), x) ./ 2)

"""
    PerEQ

The usual periodic kernel derived by mapping the input domain onto the unit circle.
"""
struct PerEQ <: Kernel end

# Binary methods.
ew(k::PerEQ, x::AV{<:Real}, x′::AV{<:Real}) = exp.(.-2 .* sin.(π .* abs.(x .- x′)).^2)
pw(k::PerEQ, x::AV{<:Real}, x′::AV{<:Real}) = exp.(.-2 .* sin.(π .* abs.(x .- x′')).^2)

# Unary methods.
ew(::PerEQ, x::AV{<:Real}) = ones(eltype(x), length(x))
pw(k::PerEQ, x::AV{<:Real}) = pw(k, x, x)

"""
    Exp <: Kernel

The standardised Exponential kernel.
"""
struct Exp <: Kernel end

# Binary methods
ew(k::Exp, x::AV, x′::AV) = exp.(.-ew(Euclidean(), x, x′))
pw(k::Exp, x::AV, x′::AV) = exp.(.-pw(Euclidean(), x, x′))

# Unary methods
ew(::Exp, x::AV) = exp.(.-ew(Euclidean(), x))
pw(::Exp, x::AV) = exp.(.-pw(Euclidean(), x))

"""
    Matern12

Equivalent to the Exponential kernel.
"""
const Matern12 = Exp

"""
    Matern32 <: Kernel

The Matern kernel with ν = 3 / 2
"""
struct Matern32 <: Kernel end

function _matern32(d)
    d = sqrt(3) * d
    return (1 + d) * exp(-d)
end

# Binary methods
ew(k::Matern32, x::AV, x′::AV) = _matern32.(ew(Euclidean(), x, x′))
pw(k::Matern32, x::AV, x′::AV) = _matern32.(pw(Euclidean(), x, x′))

# Unary methods
ew(k::Matern32, x::AV) = _matern32.(ew(Euclidean(), x))
pw(k::Matern32, x::AV) = _matern32.(pw(Euclidean(), x))

"""
    Matern52 <: Kernel

The Matern kernel with ν = 5 / 2
"""
struct Matern52 <: Kernel end

function _matern52(d)
    d = sqrt(5) * d
    return (1 + d + d^2 / 3) * exp(-d)
end

# Binary methods
ew(k::Matern52, x::AV, x′::AV) = _matern52.(ew(Euclidean(), x, x′))
pw(k::Matern52, x::AV, x′::AV) = _matern52.(pw(Euclidean(), x, x′))

# Unary methods
ew(k::Matern52, x::AV) = _matern52.(ew(Euclidean(), x))
pw(k::Matern52, x::AV) = _matern52.(pw(Euclidean(), x))

"""
    RQ <: Kernel

The standardised Rational Quadratic, with kurtosis `α`.
"""
struct RQ{Tα<:Real} <: Kernel
    α::Tα
end

_rq(d, α) = (1 + d / (2α))^(-α)

# Binary methods.
ew(k::RQ, x::AV, x′::AV) = _rq.(ew(SqEuclidean(), x, x′), k.α)
pw(k::RQ, x::AV, x′::AV) = _rq.(pw(SqEuclidean(), x, x′), k.α)

# Unary methods.
ew(k::RQ, x::AV) = _rq.(ew(SqEuclidean(), x), k.α)
pw(k::RQ, x::AV) = _rq.(pw(SqEuclidean(), x), k.α)

"""
    Linear{T<:Real} <: Kernel

The standardised linear kernel / dot-product kernel.
"""
struct Linear <: Kernel end

# Binary methods
ew(k::Linear, x::AV{<:Real}, x′::AV{<:Real}) = x .* x′
pw(k::Linear, x::AV{<:Real}, x′::AV{<:Real}) = x .* x′'
ew(k::Linear, x::ColVecs, x′::ColVecs) = reshape(sum(x.X .* x′.X; dims=1), :)
pw(k::Linear, x::ColVecs, x′::ColVecs) = x.X' * x′.X

# Unary methods
ew(k::Linear, x::AV{<:Real}) = x.^2
pw(k::Linear, x::AV{<:Real}) = x .* x'
ew(k::Linear, x::ColVecs) = reshape(sum(abs2.(x.X); dims=1), :)
pw(k::Linear, x::ColVecs) = x.X' * x.X

"""
    Poly{Tσ<:Real} <: Kernel

Inhomogeneous Polynomial kernel. `Poly(p, σ²)` creates a `Poly{p}` with variance σ²,
defined as
```julia
k(xl, xr) = (dot(xl, xr) + σ²)^p
```
"""
struct Poly{p, Tσ²<:Real} <: Kernel
    σ²::Tσ²
end
Poly(p::Int, σ²::Real) = Poly{p, typeof(σ²)}(σ²)

_poly(k, σ², p) = (σ² + k)^p
Zygote.@adjoint function _poly(k, σ², p)
    y = _poly(k, σ², p)
    return y, function(Δ)
        d = Δ * p * y / (σ² + k)
        return (d, d, nothing)
    end
end

# Binary methods
ew(k::Poly{p}, x::AV, x′::AV) where {p} = _poly.(ew(Linear(), x, x′), k.σ², p)
pw(k::Poly{p}, x::AV, x′::AV) where {p} = _poly.(pw(Linear(), x, x′), k.σ², p)

# Unary methods
ew(k::Poly{p}, x::AV) where {p} = _poly.(ew(Linear(), x), k.σ², p)
pw(k::Poly{p}, x::AV) where {p} = _poly.(pw(Linear(), x), k.σ², p)

"""
    GammaExp

The γ-Exponential kernel, 0 < γ ⩽ 2, is given by `k(xl, xr) = exp(-||xl - xr||^γ)`.
"""
struct GammaExp{Tγ<:Real} <: Kernel
    γ::Tγ
end

# Binary methods
ew(k::GammaExp, x::AV, x′::AV) = exp.(.-ew(Euclidean(), x, x′).^k.γ)
pw(k::GammaExp, x::AV, x′::AV) = exp.(.-pw(Euclidean(), x, x′).^k.γ)

# Unary methods
ew(k::GammaExp, x::AV) = exp.(.-ew(Euclidean(), x).^k.γ)
pw(k::GammaExp, x::AV) = exp.(.-pw(Euclidean(), x).^k.γ)

"""
    Wiener <: Kernel

The standardised stationary Wiener-process kernel.
"""
struct Wiener <: Kernel end

_wiener(x::Real, x′::Real) = min(x, x′)

# Binary methods
ew(k::Wiener, x::AV{<:Real}, x′::AV{<:Real}) = _wiener.(x, x′)
pw(k::Wiener, x::AV{<:Real}, x′::AV{<:Real}) = _wiener.(x, x′')

# Unary methods
ew(k::Wiener, x::AV{<:Real}) = x
pw(k::Wiener, x::AV{<:Real}) = pw(k, x, x)

"""
    WienerVelocity <: Kernel

The standardised WienerVelocity kernel.
"""
struct WienerVelocity <: Kernel end

_wiener_vel(x::Real, x′::Real) = min(x, x′)^3 / 3 + abs(x - x′) * min(x, x′)^2 / 2

# Binary methods
ew(k::WienerVelocity, x::AV{<:Real}, x′::AV{<:Real}) = _wiener_vel.(x, x′)
pw(k::WienerVelocity, x::AV{<:Real}, x′::AV{<:Real}) = _wiener_vel.(x, x′')

# Unary methods
ew(k::WienerVelocity, x::AV{<:Real}) = ew(k, x, x)
pw(k::WienerVelocity, x::AV{<:Real}) = pw(k, x, x)

"""
    Noise{T<:Real} <: Kernel

The standardised aleatoric white-noise kernel. Isn't really a kernel, but never mind...
"""
struct Noise{T<:Real} <: Kernel end
Noise() = Noise{Int}()

# Binary methods.
ew(k::Noise{T}, x::AV, x′::AV) where {T} = zeros(T, broadcast_shape(size(x), size(x′))...)
pw(k::Noise{T}, x::AV, x′::AV) where {T} = zeros(T, length(x), length(x′))

# Unary methods.
ew(k::Noise{T}, x::AV) where {T} = ones(T, length(x))
pw(k::Noise{T}, x::AV) where {T} = diagm(0=>ones(T, length(x)))

#
# Composite Kernels
#

"""
    Sum{Tkl<:Kernel, Tkr<:Kernel} <: Kernel

Represents the sum of two kernels `kl` and `kr` s.t. `k(x, x′) = kl(x, x′) + kr(x, x′)`.
"""
struct Sum{Tkl<:Kernel, Tkr<:Kernel} <: Kernel
    kl::Tkl
    kr::Tkr
end

+(kl::Kernel, kr::Kernel) = Sum(kl, kr)

# Binary methods
ew(k::Sum, x::AV, x′::AV) = ew(k.kl, x, x′) + ew(k.kr, x, x′)
pw(k::Sum, x::AV, x′::AV) = pw(k.kl, x, x′) + pw(k.kr, x, x′)

# Unary methods
ew(k::Sum, x::AV) = ew(k.kl, x) + ew(k.kr, x)
pw(k::Sum, x::AV) = pw(k.kl, x) + pw(k.kr, x)

"""
    Product{Tkl<:Kernel, Tkr<:Kernel} <: Kernel

Represents the product of two kernels `kl` and `kr` s.t. `k(x, x′) = kl(x, x′) kr(x, x′)`.
"""
struct Product{Tkl<:Kernel, Tkr<:Kernel} <: Kernel
    kl::Tkl
    kr::Tkr
end

*(kl::Kernel, kr::Kernel) = Product(kl, kr)

# Binary methods
ew(k::Product, x::AV, x′::AV) = ew(k.kl, x, x′) .* ew(k.kr, x, x′)
pw(k::Product, x::AV, x′::AV) = pw(k.kl, x, x′) .* pw(k.kr, x, x′)

# Unary methods
ew(k::Product, x::AV) = ew(k.kl, x) .* ew(k.kr, x)
pw(k::Product, x::AV) = pw(k.kl, x) .* pw(k.kr, x)

"""
    Scaled{Tσ²<:Real, Tk<:Kernel} <: Kernel

Scale the variance of `Kernel` `k` by `σ²` s.t. `(σ² * k)(x, x′) = σ² * k(x, x′)`.
"""
struct Scaled{Tσ²<:Real, Tk<:Kernel} <: Kernel
    σ²::Tσ²
    k::Tk
end

*(σ²::Real, k::Kernel) = Scaled(σ², k)
*(k::Kernel, σ²) = σ² * k

# Binary methods.
ew(k::Scaled, x::AV, x′::AV) = k.σ² .* ew(k.k, x, x′)
pw(k::Scaled, x::AV, x′::AV) = k.σ² .* pw(k.k, x, x′)

# Unary methods.
ew(k::Scaled, x::AV) = k.σ² .* ew(k.k, x)
pw(k::Scaled, x::AV) = k.σ² .* pw(k.k, x)

"""
    Stretched{Tk<:Kernel} <: Kernel

Apply a length scale to a kernel. Specifically, `k(x, x′) = k(a * x, a * x′)`.
"""
struct Stretched{Ta<:Union{Real, AV{<:Real}, AM{<:Real}}, Tk<:Kernel} <: Kernel
    a::Ta
    k::Tk
end

stretch(k::Kernel, a::Union{Real, AV{<:Real}, AM{<:Real}}) = Stretched(a, k)

# Binary methods (scalar `a`, scalar-valued input)
ew(k::Stretched{<:Real}, x::AV{<:Real}, x′::AV{<:Real}) = ew(k.k, k.a .* x, k.a .* x′)
pw(k::Stretched{<:Real}, x::AV{<:Real}, x′::AV{<:Real}) = pw(k.k, k.a .* x, k.a .* x′)

# Unary methods (scalar)
ew(k::Stretched{<:Real}, x::AV{<:Real}) = ew(k.k, k.a .* x)
pw(k::Stretched{<:Real}, x::AV{<:Real}) = pw(k.k, k.a .* x)

# Binary methods (scalar and vector `a`, vector-valued input)
function ew(k::Stretched{<:Union{Real, AV{<:Real}}}, x::ColVecs, x′::ColVecs)
    return ew(k.k, ColVecs(k.a .* x.X), ColVecs(k.a .* x′.X))
end
function pw(k::Stretched{<:Union{Real, AV{<:Real}}}, x::ColVecs, x′::ColVecs)
    return pw(k.k, ColVecs(k.a .* x.X), ColVecs(k.a .* x′.X))
end

# Unary methods (scalar and vector `a`, vector-valued input)
ew(k::Stretched{<:Union{Real, AV{<:Real}}}, x::ColVecs) = ew(k.k, ColVecs(k.a .* x.X))
pw(k::Stretched{<:Union{Real, AV{<:Real}}}, x::ColVecs) = pw(k.k, ColVecs(k.a .* x.X))

# Binary methods (matrix `a`, vector-valued input)
function ew(k::Stretched{<:AM{<:Real}}, x::ColVecs, x′::ColVecs)
    return ew(k.k, ColVecs(k.a * x.X), ColVecs(k.a * x′.X))
end
function pw(k::Stretched{<:AM{<:Real}}, x::ColVecs, x′::ColVecs)
    return pw(k.k, ColVecs(k.a * x.X), ColVecs(k.a * x′.X))
end

# Unary methods (scalar and vector `a`, vector-valued input)
ew(k::Stretched{<:AM{<:Real}}, x::ColVecs) = ew(k.k, ColVecs(k.a * x.X))
pw(k::Stretched{<:AM{<:Real}}, x::ColVecs) = pw(k.k, ColVecs(k.a * x.X))

# Create convenience versions of each of the kernels that accept a length scale.
for (k, K) in (
    (:eq, :EQ),
    (:exponential, :Exp),
    (:matern12, :Matern12),
    (:matern32, :Matern32),
    (:matern52, :Matern52),
    (:linear, :Linear),
    (:wiener, :Wiener),
    (:wiener_velocity, :WienerVelocity),
)
    @eval $k() = $K()
    @eval $k(a::Union{Real, AV{<:Real}, AM{<:Real}}) = stretch($k(), a)
    @eval export $k
end

rq(α) = RQ(α)
rq(α, l) = stretch(rq(α), l)
export rq

γ_exponential(γ::Real) = GammaExp(γ)
γ_exponential(γ::Real, l::Union{Real, AV{<:Real}, AM{<:Real}}) = stretch(GammaExp(γ), l)
export γ_exponential

poly(p::Int, σ²::Real) = Poly(p, σ²)
poly(p::Int, σ²::Real, l::Union{Real, AV{<:Real}, AM{<:Real}}) = stretch(Poly(p, σ²), l)
export poly


text:
abstract type AbstractBlockMatrix{T} <: AbstractMatrix{T} end

function nblocks end
function block end

function Base.size(bm::AbstractBlockMatrix)
    n = mapreduce(blk -> LinearAlgebra.checksquare(block(bm, blk)),
                  +, 1:nblocks(bm), init=0)
    return (n, n)
end
function Base.getindex(bm::AbstractBlockMatrix, i::Integer, j::Integer)
    (i < 0 || j < 0) && throw(BoundsError(i, j))
    for k in 1:nblocks(bm)
        blk = block(bm, k)
        n = size(blk, 1)
        if i <= n && j <= n
            return blk[i, j]
        elseif i <= n || j <= n
            return 0
        else
            i -= n
            j -= n
        end
    end
    i, j = (i, j) .+ size(bm)
    throw(BoundsError(i, j))
end
Base.getindex(A::AbstractBlockMatrix, I::Tuple) = getindex(A, I...)


text:
using OrdinaryDiffEq, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear,
       linear, f_2dlinear

@testset "save_idxs Tests" begin
  # scalar, not in-place
  prob = prob_ode_linear
  sol = solve(prob, DP5(); save_idxs=1)
  sol(0.5) # test interpolation

  # vector, not in-place
  prob2 = ODEProblem(linear, [0.5], (0., 1.), 1.01)
  sol2 = solve(prob2, DP5(); save_idxs=1)
  @test sol.t ≈ sol2.t && sol.u ≈ sol2.u
  sol2(0.5)

  sol2b = solve(prob2, DP5(); save_idxs=[1])
  @test sol.t ≈ sol2b.t && sol.u ≈ [u[1] for u in sol2b.u]
  sol2b(0.5)

  # vector, in-place
  prob3 = ODEProblem(f_2dlinear, [0.5], (0., 1.), 1.01)
  sol3 = solve(prob3, DP5(); save_idxs=1)
  @test sol.t ≈ sol3.t && sol.u ≈ sol3.u
  sol3(0.5)

  sol3b = solve(prob3, DP5(); save_idxs=[1])
  @test sol.t ≈ sol3b.t && sol2b.u ≈ sol3b.u
  sol3b(0.5)
end


text:
# We are trying to identify the source of the kink
# after ~900 days in the time integration.

include("../../src/ttv.jl")
include("/Users/ericagol/Computer/Julia/regress.jl")

using PyPlot

# This routine takes derivative of transit times with respect
# to the initial orbital elements.
n = 3
t0 = 7257.93115525
h  = 0.07
tmax = 2000.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')

ntt = zeros(Int64,n)

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
tt  = zeros(n,maximum(ntt))
tt1 = zeros(n,maximum(ntt))
tt_save = zeros(5,n,maximum(ntt))
tt_save_big = zeros(5,n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
dq = ttv_elements!(n,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
tt_save[1,:,:]=tt1

# Create BigFloat versions of the variables:
elements_big = convert(Array{BigFloat,2},elements)
hbig = big(h)
t0big = big(t0)
tmaxbig = big(tmax)
tt1big = big.(tt1)
rstarbig = big(rstar)
dqbig = ttv_elements!(n,t0big,hbig,tmaxbig,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
tt_save_big[1,:,:]=convert(Array{Float64,2},tt1big)

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dq = ttv_elements!(n,t0,h/2.,tmax,elements,tt1,count,0.0,0,0,rstar)
tt_save[2,:,:]=tt1
dqbig = ttv_elements!(n,t0big,hbig,tmaxbig/2,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
tt_save_big[2,:,:]=convert(Array{Float64,2},tt1big)
dq = ttv_elements!(n,t0,h/4.,tmax,elements,tt1,count,0.0,0,0,rstar)
tt_save[3,:,:]=tt1
dqbig = ttv_elements!(n,t0big,hbig,tmaxbig/4,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
tt_save_big[3,:,:]=convert(Array{Float64,2},tt1big)
dq = ttv_elements!(n,t0,h/8.,tmax,elements,tt1,count,0.0,0,0,rstar)
tt_save[4,:,:]=tt1
dqbig = ttv_elements!(n,t0big,hbig/8,tmaxbig,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
tt_save_big[4,:,:]=convert(Array{Float64,2},tt1big)

dq = ttv_elements!(n,t0,h/16.,tmax,elements,tt1,count,0.0,0,0,rstar)
tt_save[5,:,:]=tt1
# Compute the h/16 case in BigFloat precision:
dqbig = ttv_elements!(n,t0big,hbig/16,tmaxbig,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
tt_save_big[5,:,:]=convert(Array{Float64,2},tt1big)

# Make a plot of transit time errors versus stepsize:
ntrans = sum(count)
clf()
sigt = zeros(n-1,5)
tab = 0
#h_list = [h,h/2.,h/4.,h/8.]
h_list = [h,h/2.,h/4.,h/8,h/16]
hlabel = ["h-big(h)","h/2-big(h/2)","h/4-big(h/4)","h/8-big(h/8)","h/16-big(h/16)"]
ch = ["black","red","green","blue"]
for i=2:n
  for j=1:5
    tti1 = tt_save[j,i,1:count[i]]
    tti16 = tt_save_big[j,i,1:count[i]]
    diff = tti1-tti16
    sigt[i-1,j]=std(diff)
    if i == n
      plot(tti1,diff,linestyle="dashed",c=ch[j])
    else
      plot(tti1,diff,label=hlabel[j],c=ch[j])
    end
  end
end
xlabel("Transit time")
ylabel("Difference in transit time")
legend(loc="lower left")

#read(STDIN,Char)
#
## Make a plot of timing errors versus stepsize:
#ntrans = sum(count)
#clf()
#sigt = zeros(n-1,4)
#tab = 0
##h_list = [h,h/2.,h/4.,h/8.]
#h_list = [h,h/2.,h/4.,h/8]
#hlabel = ["h-h/16","h/2-h/16","h/4-h/16","h/8-h/16"]
#ch = ["black","red","green","blue"]
#for i=2:n
#  fn = zeros(Float64,2,count[i])
#  sig = ones(count[i])
#  tti16 = tt_save[5,i,1:count[i]]
#  for j=1:count[i]
#    fn[1,j] = 1.0
#    fn[2,j] = round(Int64,(tti16[j]-elements[i,3])/elements[i,2])
#  end
#  for j=1:4
#    tti1 = tt_save[j,i,1:count[i]]
#    coeff,cov = regress(fn,tti1-tti16,sig)
#    diff = tti1-tti16-coeff[1]-coeff[2]*fn[2,:]
#    sigt[i-1,j]=std(diff)
#    if i == n
#      plot(tti1,diff/h_list[j]^4,linestyle="dashed",c=ch[j])
#    else
#      plot(tti1,diff/h_list[j]^4,label=hlabel[j],c=ch[j])
#    end
#  end
#end
#legend(loc="lower left")
#
#read(STDIN,Char)
#clf()
#
#
## Make a plot of some TTVs:
#loglog(h_list,sigt[1,:]*24.*3600.,".",markersize=15,label="Inner planet")
#loglog(h_list,sigt[1,1]*24.*3600.*(h_list/h[1]).^4,label=L"$\propto h^4$")
#loglog(h_list,sigt[2,:]*24.*3600.,".",markersize=15,label="Outer planet")
#loglog(h_list,sigt[2,1]*24.*3600.*(h_list/h[1]).^4,label=L"$\propto h^4$")
#legend(loc = "upper left")
#ylabel("RMS timing error [sec]")
#xlabel("Step size [day]")
#
#PyPlot.savefig("timing_error_vs_h.pdf",bbox_inches="tight")

max_forks_repo_name: langfzac/nbg-papers


text:
@everywhere begin
    struct SimpleScenario <: StochasticPrograms.AbstractScenarioData
        π::Float64
        d::Vector{Float64}
        q::Vector{Float64}
    end

    function StochasticPrograms.expected(sds::Vector{SimpleScenario})
        sd = SimpleScenario(1,sum([s.π*s.d for s in sds]),sum([s.π*s.q for s in sds]))
    end
end

s1 = SimpleScenario(0.4,[500.0,100],[-24.0,-28])
s2 = SimpleScenario(0.6,[300.0,300],[-28.0,-32])

sds = [s1,s2]

sp = StochasticProgram(sds)

@first_stage sp = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁+x₂ <= 120)
end

@second_stage sp = begin
    @decision x₁ x₂
    s = scenario
    @variable(model, 0 <= y₁ <= s.d[1])
    @variable(model, 0 <= y₂ <= s.d[2])
    @objective(model, Min, s.q[1]*y₁ + s.q[2]*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end

push!(problems,(sp,"Simple"))


text:
function uri2filepath(uri::URI)
    if uri.scheme != "file"
        return nothing
    end

    path = uri.path
    host = uri.authority

    if host!==nothing && host != "" && length(path) > 1
        # unc path: file://shares/c$/far/boo
        value = "//$host$path"
    elseif length(path) >= 3 &&
            path[1] == '/' &&
            isascii(path[2]) && isletter(path[2]) &&
            path[3] == ':'
        # windows drive letter: file:///c:/far/boo
        value = lowercase(path[2]) * path[3:end]
    else
        # other path
        value = path
    end

    if Sys.iswindows()
        value = replace(value, '/' => '\\')
    end

    return value
end

function filepath2uri(path::String)
    isabspath(path) || error("Relative path `$path` is not valid.")

    path = normpath(path)

    if Sys.iswindows()
        path = replace(path, "\\" => "/")
    end

    authority = ""

    if startswith(path, "//")
        # UNC path //foo/bar/foobar
        idx = findnext("/", path, 3)
        if idx===nothing
            authority = path[3:end]
            path = "/"
        else
            authority = path[3:idx.start-1]
            path = path[idx.start:end]
        end
    elseif length(path)>=2 && isascii(path[1]) && isletter(path[1]) && path[2]==':'
        path = string('/', lowercase(path[1]), SubString(path, 2))
    end

    return URI(scheme="file", authority=authority, path=path)
end


text:
using Documenter, BioTools

makedocs(
    format = :html,
    modules = [BioTools.BLAST],
    sitename = "BioTools.jl",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "BLAST" => "blast.md"
    ],
    authors = "The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/BioTools.jl.git",
    deps = nothing,
    make = nothing
)


text:
@everywhere using SMM
@everywhere using Distributions
@everywhere using DataStructures

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# Tests that are too long to be sent to TRAVIS
#---------------------------------------------
# 6d problem
#--------------------------------------------------------------
@testset "Testing smmoptimize with 6d and same magnitude" begin

    # Rermark:
    # It is important NOT to use srand()
    # within the function simulate_empirical_moments!
    # Otherwise, BlackBoxOptim does not find the solution
    #----------------------------------------------------
    # The difference of magniture make it more difficult to find the minimum
    tol6dMean = 0.5

    @everywhere function functionTest6d(x)

        d = MvNormal(x, eye(6))

        output = OrderedDict{String,Float64}()

        draws = rand(d, 1000000)
        for i in 1:6
            output["mean$(i)"] = mean(draws[i,:])
        end

        return output
    end

    t = SMMProblem(options = SMMOptions(maxFuncEvals=1000,saveSteps = 1000, saveName = "parallel6D"))

    #------------------------------------------------------
    dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
    dictEmpiricalMoments["mean1"] = [ -1.0; -1.0]
    dictEmpiricalMoments["mean2"] = [ 1.0; 1.0]
    dictEmpiricalMoments["mean3"] = [ 0.5; 0.5]
    dictEmpiricalMoments["mean4"] = [ -0.5; -0.5]
    dictEmpiricalMoments["mean5"] = [ 0.7; 0.7]
    dictEmpiricalMoments["mean6"] = [-0.7; -0.7]
    set_empirical_moments!(t, dictEmpiricalMoments)

    dictPriors = OrderedDict{String,Array{Float64,1}}()
    dictPriors["mu1"] = [0.2,-3., 3.]
    dictPriors["mu2"] = [-0.2,-2.,2.]
    dictPriors["mu3"] = [-0.3,-2.,2.]
    dictPriors["mu4"] = [-0.4,-2., 2.]
    dictPriors["mu5"] = [0.3,-2., 2.]
    dictPriors["mu6"] = [0.4,-2., 2.]

    set_priors!(t, dictPriors)

    # A. Set the function: parameter -> simulated moments
    set_simulate_empirical_moments!(t, functionTest6d)

    # B. Construct the objective function, using the function: parameter -> simulated moments
    # and moments' weights:
    construct_objective_function!(t)

    # C. Run the optimization
    # This function first modifies t.bbSetup
    # and then modifies t.bbResults

    # [TODO] loading and restarting does
    # not seem to work when using workers in parallel
    # Figure out why
    listBestFitness, listBestCandidates = smmoptimize!(t, verbose = true)

    @test best_candidate(t.bbResults)[1] ≈ -1.0 atol = tol6dMean
    @test best_candidate(t.bbResults)[2] ≈ 1.0 atol = tol6dMean
    @test best_candidate(t.bbResults)[3] ≈ 0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[4] ≈ -0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[5] ≈ 0.7 atol = tol6dMean
    @test best_candidate(t.bbResults)[6] ≈ -0.7 atol = tol6dMean

end

# Ambitious problem, with 18 parameters
# different magnitudes
#--------------------------------------
# source: https://github.com/floswald/MomentOpt.jl/blob/master/src/mopt/Examples.jl
#----------------------------------------------------------------------------------
@testset "Testing smmoptimize with 18d and different magnitude" begin

    # Rermark:
    # It is important NOT to use srand()
    # within the function simulate_empirical_moments!
    # Otherwise, BlackBoxOptim does not find the solution
    #----------------------------------------------------
    # The difference of magniture make it more difficult to find the minimum
    tol18dMean = 0.5

    function functionTest18d(x)

        d = MvNormal(x, eye(18))

        output = OrderedDict{String,Float64}()

        draws = rand(d, 1000000)
        for i in 1:18
            output["mean$(i)"] = mean(draws[i,:])
        end

        return output
    end

    t = SMMProblem(options = SMMOptions(maxFuncEvals=10000,saveSteps = 1000))

    #------------------------------------------------------
    dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
    dictEmpiricalMoments["mean1"] = [ -1.0; -1.0]
    dictEmpiricalMoments["mean2"] = [ 0.01; 0.01]
    dictEmpiricalMoments["mean3"] = [ 0.5; 0.5]
    dictEmpiricalMoments["mean4"] = [ -0.5; -0.5]
    dictEmpiricalMoments["mean5"] = [ 0.008; 0.008]
    dictEmpiricalMoments["mean6"] = [-0.7; -0.7]

    dictEmpiricalMoments["mean7"] = [-1.0; -1.0]
    dictEmpiricalMoments["mean8"] = [0.01; 0.01]
    dictEmpiricalMoments["mean9"] = [ 0.5; 0.5]
    dictEmpiricalMoments["mean10"] = [ -0.5; -0.5]
    dictEmpiricalMoments["mean11"] = [ 0.008; 0.008]
    dictEmpiricalMoments["mean12"] = [-0.7; -0.7]

    dictEmpiricalMoments["mean13"] = [-1.0; -1.0]
    dictEmpiricalMoments["mean14"] =  [0.01; 0.01]
    dictEmpiricalMoments["mean15"] = [ 0.5; 0.5]
    dictEmpiricalMoments["mean16"] = [ -0.5; -0.5]
    dictEmpiricalMoments["mean17"] = [ 0.008; 0.008]
    dictEmpiricalMoments["mean18"] = [-0.7; -0.7]

    set_empirical_moments!(t, dictEmpiricalMoments)

    pb = OrderedDict{String,Array{Float64,1}}()

    pb["p1"] = [0.2, 0., 1.]
    pb["p2"] = [0.002, 0., 0.02]
    pb["p3"] = [-0.3,-2., 2.]
    pb["p4"] = [-0.4,-2., 2.]
    pb["p5"] = [0.008,0.,0.02]
    pb["p6"] = [0.4,-2., 2.]

    pb["p7"] = [0.2,0. ,1.]
    pb["p8"] = [0.002,0., 0.02]
    pb["p9"] = [-0.3,-2., 2.]
    pb["p10"] = [-0.4,-2., 2.]
    pb["p11"] = [0.008,0,0.02]
    pb["p12"] = [0.4,-2., 2.]

    pb["p13"] = [0.2, 0., 1.]
    pb["p14"] = [0.002,0,0.02]
    pb["p15"] = [-0.3,-2., 2.]
    pb["p16"] = [-0.4,-2., 2.]
    pb["p17"] = [0.008,0., 0.02]
    pb["p18"] = [0.4,-2., 2.]

    set_priors!(t, pb)

    # A. Set the function: parameter -> simulated moments
    set_simulate_empirical_moments!(t, functionTest18d)

    # B. Construct the objective function, using the function: parameter -> simulated moments
    # and moments' weights:
    construct_objective_function!(t)

    # C. Run the optimization
    # This function first modifies t.bbSetup
    # and then modifies t.bbResults
    smmoptimize!(t, verbose = true)

    @test best_candidate(t.bbResults)[1] ≈ -1.0 atol = tol6dMean
    @test best_candidate(t.bbResults)[2] ≈ 0.01 atol = tol6dMean
    @test best_candidate(t.bbResults)[3] ≈ 0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[4] ≈ -0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[5] ≈ 0.008 atol = tol6dMean
    @test best_candidate(t.bbResults)[6] ≈ -0.7 atol = tol6dMean

    @test best_candidate(t.bbResults)[7] ≈ -1.0 atol = tol6dMean
    @test best_candidate(t.bbResults)[8] ≈ 0.01 atol = tol6dMean
    @test best_candidate(t.bbResults)[9] ≈ 0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[10] ≈ -0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[11] ≈ 0.008 atol = tol6dMean
    @test best_candidate(t.bbResults)[12] ≈ -0.7 atol = tol6dMean

    @test best_candidate(t.bbResults)[13] ≈ -1.0 atol = tol6dMean
    @test best_candidate(t.bbResults)[14] ≈ 0.01 atol = tol6dMean
    @test best_candidate(t.bbResults)[15] ≈ 0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[16] ≈ -0.5 atol = tol6dMean
    @test best_candidate(t.bbResults)[17] ≈ 0.008 atol = tol6dMean
    @test best_candidate(t.bbResults)[18] ≈ -0.7 atol = tol6dMean

end


text:
using SparseDiffTools
using DiffEqDiffTools: finite_difference_jacobian, finite_difference_jacobian!
using SparsityDetection

using LinearAlgebra, SparseArrays, Test

fcalls = 0

function f(x)
  global fcalls += 1
  dx = similar(x)
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  dx
end

function f(dx,x)
  global fcalls += 1
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      (j-i==-1 || j-i==1) && (A[i,j]=1)
      j-i==0 && (A[i,j]=-2)
  end
  A
end

output = ones(30); input = ones(30)
sparsity_pattern = sparsity!(f,output,input)
true_jac = Float64.(sparse(sparsity_pattern))
colors = matrix_colors(true_jac)
@test colors == repeat(1:3,10)

#Jacobian computed without coloring vector
fcalls = 0
J = finite_difference_jacobian(f, rand(30))
@test J ≈ second_derivative_stencil(30)
@test fcalls == 31

#Jacobian computed with coloring vectors
fcalls = 0
_J = similar(true_jac)
finite_difference_jacobian!(_J, f, rand(30), colorvec = colors)
@test fcalls == 4
@test _J ≈ J

fcalls = 0
_J = similar(true_jac)
_denseJ = collect(_J)
finite_difference_jacobian!(_denseJ, f, rand(30), colorvec = colors, sparsity=_J)
@test fcalls == 4
@test _denseJ ≈ J


text:
"""
Modia module with utility functions.

* Developer: Hilding Elmqvist, Mogram AB
* First version: July 2016
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

"""
module Utilities

using ..Instantiation: GetField, Der, InitVariable, Extends, Instance, Variable, Variability, constant, parameter, discrete, continuous, Connect, vars_of
using ..Execution: split_variables
using Base.Meta: quot, isexpr
using Unitful
using ..Instantiation

using ..ModiaLogging

export @show_io, printSymbolList, showModel, showInstance, checkSizes, printJSON

using Base: show_unquoted

# Version of @show for any stream
macro show_io(io, exs...)
    blk = Expr(:block)
    for ex in exs
        push!(blk.args, :(loglnModia($io, $(sprint(show_unquoted, ex) * " = "),
                                  repr(begin value = $(esc(ex)) end))))
    end
    if !isempty(exs); push!(blk.args, :value); end
    return blk
end

const maxSymbols = 1000

function printSymbolList(label, symbols, numbering=false, vertical=false, A=[])
    logModia(label, ": ")

    if vertical
        loglnModia()
    end

    for i in 1:min(maxSymbols, length(symbols))
        if i > 1
            if vertical
                loglnModia()
            else
                logModia(", ")
            end
        end

        if numbering
            if vertical
                logModia(lpad(i, 5, " "), ": ")
            else
                logModia(i, ": ")
            end
        end

        if vertical
            logModia(rpad(prettyfy(symbols[i]), 30, " "))
        else
            logModia(prettyfy(symbols[i]))
        end

        if A != []
            if A[i] != 0
                logModia("A[$i] = $(A[i])")
            end
        end
    end

    if length(symbols) > maxSymbols
        logModia(", ...")
    end

    loglnModia()
end

function showModel(mod, indent="")
    loglnModia("showModel:")
    #=
    loglnModia("------")
    @show mod
    loglnModia("------")
    =#
    newIndent = string(indent, "  ")
    loglnModia("Instance(model ", mod.name)

    for i in 1:length(mod.initializers)
        if isa(mod.initializers[i], InitVariable)
            logModia(indent, mod.initializers[i].name, " = ")
            loglnModia(mod.initializers[i].init)
        else
            loglnModia(indent, mod.initializers[i])
            if isa(mod.initializers[i], Extends)
                # println("Extends($(mod.initializers[i].fdef.args[2].args[2].args[2].args[1]))")
                # Check that we have an Expr
                if isa(mod.initializers[i].fdef.args[2].args[2].args[2], Expr)
                    partial = mod.initializers[i].fdef.args[2].args[2].args[2].args[1] == :PARTIAL
                    if partial
                        # mod.partial = true
                    end
                else
                    error("Missing parentheses in extends clause: @extends $(mod.initializers[i].fdef.args[2].args[2].args[2])")
                end
            end
        end
    end
    loglnModia(indent, "end)\n")
end

function showVariable(v)
    if isa(v, Variable)
        logModia("Variable(")
        first = true

        if v.variability != continuous
            if !first; logModia(", ") end
            logModia("variability = ", v.variability)
            first = false
        end

        if v.typ != Any
            if !first; logModia(", ") end
            logModia("T = ", v.typ)
            first = false
        end

        if v.value != nothing
            if !first; logModia(", ") end
            logModia("value = ", v.value)
            first = false
        end

        if v.size != nothing
            if !first; logModia(", ") end
            logModia("size = ", v.size)
            first = false
        end

        if v.start != nothing
            if !first; logModia(", ") end
            logModia("start = ", v.start)
            first = false
        end

        if v.unit != NoUnits
            if !first; logModia(", ") end
            logModia("unit = ", v.unit)
            first = false
        end

        if v.min != nothing
            if !first; logModia(", ") end
            logModia("min = ", v.min)
            first = false
        end

        if v.max != nothing
            if !first; logModia(", ") end
            logModia("max = ", v.max)
            first = false
        end

        if v.flow != nothing && v.flow
            if !first; logModia(", ") end
            logModia("flow = ", v.flow)
            first = false
        end

        if v.state != true
            if !first; logModia(", ") end
            logModia("state = ", v.state)
            first = false
        end

        if v.info != ""
            if !first; logModia(", ") end
            logModia("info = \"", v.info, "\"")
            first = false
        end

        loglnModia(")")
    else
        loglnModia(v)
    end
end

function showInstance(inst, indent="")
    # loglnModia("showInstance:")
    #=
    loglnModia("------")
    @show inst
    loglnModia("------")
    =#
    newIndent = string(indent, "  ")
    loglnModia("@model ", inst.model_name, " begin")
    for key in keys(inst.variables)
        if isa(key, Instance)
            logModia(indent, "  ", key, " = ")
            showInstance(inst.variables[key], newIndent)
        else
            @static if VERSION < v"0.7.0-DEV.2005"
                keyname = replace(string(key), ".", "_")
            else
                keyname = replace(string(key), "." => "_")
            end
            logModia(indent, "  ", keyname, " = ")
            showVariable(inst.variables[key])
        end
    end

    if length(inst.equations) > 0
        loglnModia(indent, "@equations begin")
    end

    for e in inst.equations
        loglnModia(indent, "  ", prettyPrint(e))
    end

    loglnModia(indent, "  end")
    loglnModia(indent, "end")
end

function checkSizes(VSizes, ESizes)
    # Check that system matrix is square
    ok = true
    if length(VSizes) != length(ESizes)
        ModiaLogging.increaseLogCategory(:UnbalancedModel)
        # For TestArray5
        #    error("The number of unknowns ($(length(VSizes))) is not equal to the number of equations ($(length(ESizes)))")
    end

    loglnModia()
    # Should check that arrays are not empty
    scalarV = sum(length(zeros(v)) for v in VSizes)
    scalarE = sum(length(zeros(e)) for e in ESizes)

    if scalarV != scalarE
        # error("Scalarized system matrix is not square: $scalarE x $scalarV")
        ModiaLogging.closeLogModia()
        error("The number of scalarized unknowns (= $scalarV) is not equal to the number of scalarized equations (= $scalarE).\n",
              "If option `simulate(<model>, ...; logTranslation=true)` is set, inspect <user>/ModiaResults/<model>.txt for more info.")
        ok = false
    else
        # loglnModia("Scalarized system matrix is square: $scalarE x $scalarV")
        loglnModia("The number of scalarized unknowns (= $scalarV) is equal to the number of scalarized equations (= $scalarE).")
    end

    if false && (sort(VSizes) != sort(ESizes) )  ### Testing
        @show sort(VSizes) sort(ESizes)
        # Show the number of variables (and equations) of different sizes
        differentVSizes = unique(sort(VSizes))
        VSizeCount = [length([s for s in VSizes if s == d ]) for d in differentVSizes]
        differentESizes = unique(sort(ESizes))
        ESizeCount = [length([s for s in ESizes if s == d ]) for d in differentESizes]
        @show differentVSizes VSizeCount differentESizes ESizeCount
        ok = false
        error("The sizes of variables and equations don't match.")
    end
end

# ---------------------------------------------------------------

# Exprimental code for outputting connection structure to be used by Kieler automatic layout algorithm.

function isPort(inst)
    for (n, v) in inst.variables
        if isa(v, Variable) && v.flow
            return true
        end
    end
    return false
end

const size = 30
const sizePort = 3

function printJSON(file, inst::Instance, fullName, name, parent="", indent="", level=1)
    loglnModia(file, indent, "{")
    loglnModia(file, indent, "\"id\": \"$(fullName)\",")
    loglnModia(file, indent, "\"class\": \"$(inst.model_name)\",")

    params, unknowns = split_variables(vars_of(inst))
    if length(params) > 0
        loglnModia(file, indent, "\"parameters\": {")
        indent1 = "  " * indent
        i = 0
        for (n, p) in params
            i += 1
            logModia(file, indent1, "\"", n, "\"", " : ", p)
            if i < length(params)
                loglnModia(file, ",")
            else
                loglnModia(file)
            end
        end
        loglnModia(file, indent, "},")
    end

    loglnModia(file, indent, "\"labels\": [{\"text\": \"$name\"}],")

    if isPort(inst)
        orientation = if length(search(string(parent), "G")) > 0; "\"NORTH\"" elseif length(search(string(name), "p")) > 0 || length(search(string(name), "in")) > 0; "\"WEST\"" else "\"EAST\"" end
        loglnModia(file, indent, "\"properties\": {\"de.cau.cs.kieler.portSide\": $orientation},")
        loglnModia(file, indent, "\"width\": $sizePort,")
        loglnModia(file, indent, "\"height\": $sizePort")
    else
        loglnModia(file, indent, "\"properties\": {\"de.cau.cs.kieler.portConstraints\": \"FIXED_SIDE\"},")
        loglnModia(file, indent, "\"width\": $size,")
        loglnModia(file, indent, "\"height\": $size")
    end

    first = true
    for (n, v) in inst.variables
        if isa(v, Instance) && isPort(v) && level > 1 # length(search(string(n), "n")) == 0
            if first
                loglnModia(file, indent, ", \"ports\": [")
            else
                loglnModia(file, indent, ", ")
            end
            first = false
            printJSON(file, v, string(name) * string(n), n, name, "  " * indent, level + 1)
        end
    end

    if !first
        loglnModia(file, indent, "]")
    end

    first = true
    for (n, v) in inst.variables
        if isa(v, Instance) && (!isPort(v) || level == 1) #  || (isa(v, Instance) && isPort(v) && length(search(string(n), "n")) > 0)
            if first
                loglnModia(file, indent, ", \"children\": [")
            else
                loglnModia(file, indent, ", ")
            end
            first = false
            printJSON(file, v, n, n, name, "  " * indent, level + 1)
        end
    end

    if !first
        loglnModia(file, indent, "]")
    end

    first = true
    id = 0
    for e in inst.equations
        if isa(e, Connect)
            if first
                loglnModia(file, indent, ", \"edges\": [")
            else
                loglnModia(file, indent, ", ")
            end

            first = false
            id += 1
            loglnModia(file, indent, "  {\"id\": \"id$id\",")

            if isa(e.a.base, GetField)
                loglnModia(file, indent, "   \"source\": \"$(e.a.base.name)\",")
                loglnModia(file, indent, "   \"sourcePort\": \"$(string(e.a.base.name) * string(e.a.name))\",")
            else
                loglnModia(file, indent, "   \"source\": \"$(e.a.name)\",")
            end

            if isa(e.b.base, GetField)
                loglnModia(file, indent, "   \"target\": \"$(e.b.base.name)\",")
                loglnModia(file, indent, "   \"targetPort\": \"$(string(e.b.base.name) * string(e.b.name))\"}")
            else
                loglnModia(file, indent, "   \"target\": \"$(e.b.name)\"}")
            end
        end
    end

    if !first
        loglnModia(file, indent, "]")
    end

    loglnModia(file, indent, "}")
end

end


text:
import SymbolicUtils: symtype, term, hasmetadata
struct MTKParameterCtx end

isparameter(x::Num) = isparameter(value(x))
isparameter(x::Symbolic) = getmetadata(x, MTKParameterCtx, false)
isparameter(x) = false

"""
    toparam(s::Sym)

Maps the variable to a paramter.
"""
toparam(s::Symbolic) = setmetadata(s, MTKParameterCtx, true)
toparam(s::Num) = Num(toparam(value(s)))

"""
    tovar(s::Sym)

Maps the variable to a state.
"""
tovar(s::Symbolic) = setmetadata(s, MTKParameterCtx, false)
tovar(s::Num) = Num(tovar(value(s)))

"""
$(SIGNATURES)

Define one or more known variables.
"""
macro parameters(xs...)
    Symbolics._parse_vars(:parameters,
                          Real,
                          xs,
                          x -> x isa Array ? toparam.(x) : toparam(x)
                         ) |> esc
end


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, true),),
    ((2, false),),
    ]
insts["fast"] = [
    ((20, true),),
    ((20, false),),
    ((30, true),),
    ((30, false),),
    ]
insts["various"] = vcat(insts["fast"], [
    ((40, true), nothing, relaxed_tols),
    ((40, false), nothing, relaxed_tols),
    ])
return (StabilityNumber, insts)


text:
#Splits  output of BLAST parsing by readblastXML from BioTools to groups matching one query  (queryname - key of a output dictionary)
function split_by_query(blastres::Array{BLASTResult,1})
    out = Dict{String,Array{BLASTResult,1}}()
    for r in blastres
        if !(r.queryname in keys(out))
            out[r.queryname] = Array{BLASTResult,1}()
        end
        push!(out[r.queryname], r)
    end
    return(out)
end

#check the length of a new interval that are not covered by the one in the array
function unique_length_fraction(matches::Array{Interval{:closed,:closed,Int64},1}, interval::Interval{:closed,:closed,Int64})
    difs = Array{Interval{:closed,:closed,Int64},1}()
    unique_length = interval.right - interval.left + 1
    length_all = deepcopy(unique_length)
    for i in matches
        common = intersect(interval, i)
        if !isempty(common)
            unique_length -= common.right - common.left + 1
        end
    end
    return(unique_length/length_all)
end

#add new interval to interval sets
function add_new_interval(matches::Array{Interval{:closed,:closed,Int64},1}, interval::Interval{:closed,:closed,Int64})
    intersecting = Array{Interval{:closed,:closed,Int64},1}([interval])
    has_intersecting = false

    for i in matches
        if !isempty(intersect(interval, i))
            push!(intersecting, i)
            has_intersecting = true
        end
    end
    #get union of intersecring intervals
    if has_intersecting
        expanded = deepcopy(interval)
        for i in intersecting
            expanded = union(expanded,i)
        end
    end
    out = Array{Interval{:closed,:closed,Int64},1}()
    #collect output
    if has_intersecting
        for i in matches
            if isempty(intersect(i, expanded))
                push!(out, i)
            end
        end
        push!(out, expanded)
    else
        expanded = deepcopy(matches)
        push!(expanded, interval)
        out = expanded
    end
    return(sort!(out, by = x-> x.left))
end

function covered_length(intervals::Array{Interval{:closed,:closed,Int64},1})
    out = 0
    for i in intervals
        out += i.right - i.left +1
    end
    return(out)
end

#collects set of hits that extends the coverage. Only hits that overlaps witha certain length limit with higher scoring hits are considered
function get_coverage_maximising_hits(queryid::String, data::Dict{String,Array{BLASTResult,1}}, max_overlap_frac = 0.5)
    hits = data[queryid]
    hits = sort(hits, by = x -> (x.queryto - x.queryfrom +1), rev = true )
    coverge_expanding_hits = Array{BLASTResult,1}([hits[1]])
    coverge_expanding_hits_intervals = [Interval(hits[1].queryfrom,hits[1].queryto)]
    current_length = covered_length([Interval(hits[1].queryfrom, hits[1].queryto)])
    for h in hits[2:end]
        query_interval = Interval(h.queryfrom,h.queryto)
        intervals_with = add_new_interval(coverge_expanding_hits_intervals, query_interval)
        uniqfrac = unique_length_fraction(coverge_expanding_hits_intervals, query_interval)
        new_length = covered_length(intervals_with)
        if new_length > current_length && uniqfrac < max_overlap_frac
            current_length = new_length
            push!(coverge_expanding_hits, h)
            coverge_expanding_hits_intervals = intervals_with
        end
    end
    return(coverge_expanding_hits)
end