text:
"""
    execute_sts_weekly_api_call(user_model::PSUserModel, stock_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing, AbstractLogger} = nothing) -> (Union{PSResult{T}, Nothing} where T<:Any)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_weekly_api_call(user_model::PSUserModel, stock_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing, AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # use alpha_vantage_api to access data change function to monthly in this case
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_WEEKLY&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger -
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        if (data_type == :json)
            # process json data -
            data_series_key = "Weekly Time Series"
            return process_raw_json_api_data_sts(api_call_raw_data, data_series_key)

        elseif (data_type == :csv)
            # process csv data -
            return process_raw_csv_api_data(api_call_raw_data)
        else
            # formulate the error message -
            error_message = "$(data_type) is not supported. Supported values are {:json,:csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSError(error)
    end
end

"""
    execute_sts_adjusted_weekly_api_call(user_model::PSUserModel, stock_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing, AbstractLogger} = nothing) -> (Union{PSResult{T}, Nothing} where T<:Any)

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_adjusted_weekly_api_call(user_model::PSUserModel, stock_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing, AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # use alpha_vantage_api to get data
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_WEEKLY_ADJUSTED&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger -
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        # call different handlers depending upon type -
        if (data_type == :json)
            # process json -
            data_series_key = "Weekly Adjusted Time Series"
            return process_raw_json_api_data_sts_adjusted(api_call_raw_data, data_series_key)
        elseif (data_type == :csv)
            return process_raw_csv_api_data(api_call_raw_data)
        else
            # formulate the error message -
            error_message = "$(data_type) is not supported by AlphaVantage. Supported values are {:json,:csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end


text:
check_bitcode(buf1, buf2) = buf1.data == buf2.data

function check_result(ctx, ast, asm)

    mod1 = LLVM.module_from_assembly(ctx, asm)
    mod2 = LLVM.module_from_ast(ctx, ast)

    # check that we can correctly serialize both to assembly / bitcode
    gold = LLVM.module_to_assembly(ctx, mod1)
    test = LLVM.module_to_assembly(ctx, mod2)

    Base.print_with_color(:yellow, gold)
    Base.print_with_color(gold == test ? :green : :red, test)
    Base.print_with_color(:blue, "-"^80 * "\n")

    @fact LLVM.module_to_assembly(ctx, mod1) => LLVM.module_to_assembly(ctx, mod2)
    @fact check_bitcode(LLVM.module_to_bitcode(ctx, mod1),
                        LLVM.module_to_bitcode(ctx, mod2)) => true

    # check that we can lift from assembly to an
    # Ast.Module and perform equality checks on the ast
    @fact LLVM.module_to_ast(ctx, mod1) => ast

    # check that we can produce a valid assembly string from the ast
    @fact LLVM.module_to_assembly(ctx, mod2) => asm

    LLVM.FFI.dispose_module(mod1)
    LLVM.FFI.dispose_module(mod2)
end

test_float_roundtrip{T}(::Type{T}, ctx, n) = begin
    for i = 1:n
        v  = rand(T)
        a1 = LLVM.Ast.ConstFloat(v)

        es = LLVM.EncodeState(ctx)
        c  = LLVM.encode_llvm(es, a1)

        ds = LLVM.DecodeState(ctx)
        a2 = LLVM.decode_llvm(ds, c)

        @assert a1 == a2
    end
end

ext: jl


text:
"""
Forward Euler integration scheme
"""
mutable struct Euler <: AbstractIntegrator
  system::System
  model::AbstractPotential
  function Euler(system::System,model::AbstractPotential)
    new(system,model)
  end
end

function runtraj!(nsteps::Int64,dt::Float64,int::Euler)
  for i in 1:nsteps
    step!(i,dt,int)
  end
end

function runtraj!(nsteps::Int64,dt::Float64,int::Euler,cb::AbstractCallback)
  force!(int.system,int.model)
  for i in 1:nsteps
    if i%cb.every==0
      callback(cb,int.system,int.model,i,dt)
    end
    step!(i,dt,int)
  end
end

function step!(ind::Int64,dt::Float64,int::Euler)
  # forward euler integrator
  # drift
  force!(int.system,int.model)
  addforce!(dt,int.system,int.model)
  int.system.t += dt
end

"""
Stochastic Euler integration scheme
"""
mutable struct StochasticEuler <: AbstractIntegrator
  system::AbstractThermostattedSystem
  model::AbstractPotential
  function StochasticEuler(system::AbstractThermostattedSystem,model::AbstractPotential)
    new(system,model)
  end
end

function runtraj!(nsteps::Int64,dt::Float64,int::StochasticEuler)
  for i in 1:nsteps
    step!(i,dt,int)
  end
end

function step!(ind::Int64,dt::Float64,int::StochasticEuler)
  # forward euler-murayama integrator
  # compute force for drift
  force!(int.system,int.model)
  # generate noise for diffusion
  gennoise!(int.system.dim,int.system.thermostat)
  # update with drift
  addforce!(dt,int.system,int.model)
  # update with diffusion
  int.system.x .+= sqrt(dt)*int.system.thermostat.scale .* int.system.thermostat.rands
  # update time
  int.system.t += dt
end

function runtraj!(nsteps::Int64,dt::Float64,int::StochasticEuler,cb::AbstractCallback)
  for i in 1:nsteps
    step!(i,dt,int,cb)
  end
end

function step!(ind::Int64,dt::Float64,int::StochasticEuler,cb::AbstractCallback)
  # forward euler-murayama integrator
  # compute force for drift
  force!(int.system,int.model)
  # generate noise for diffusion
  gennoise!(int.system.dim,int.system.thermostat)
  # callback
  if ind%cb.every==0
    callback(cb,int.system,int.model,ind,dt)
  end
  # update with drift
  addforce!(dt,int.system,int.model)
  # update with diffusion
  applythermostat!(dt,int.system,int.system.thermostat)
  #int.system.x .+= (sqrt(dt) .* int.system.thermostat.scale .* int.system.thermostat.rands)
  # update time
  int.system.t += dt
end


text: # TODO solve SOS problem to integrate over basic semialgebraic set https://arxiv.org/pdf/0807.2505.pdf


text:
function set_lower_bound(info::Union{VariableInfo, ConstraintInfo}, value::Union{Number, Nothing})
    if value !== nothing
        info.lower_bound !== nothing && throw(ArgumentError("Lower bound has already been set"))
        info.lower_bound = value
    end
    return
end

function set_upper_bound(info::Union{VariableInfo, ConstraintInfo}, value::Union{Number, Nothing})
    if value !== nothing
        info.upper_bound !== nothing && throw(ArgumentError("Upper bound has already been set"))
        info.upper_bound = value
    end
    return
end

function is_empty(model::BaronModel)
    isempty(model.variable_info) && isempty(model.constraint_info)
end

function set_unique_names!(infos, default_base_name::AbstractString)
    names = Set(String[])
    default_name_counter = Ref(1)
    for info in infos
        if info.name === nothing
            base_name = default_base_name
            name_counter = default_name_counter
        elseif info.name ∉ names
            push!(names, info.name)
            continue
        else
            base_name = info.name
            name_counter = Ref(1)
        end
        while true
            name = string(base_name, name_counter[])
            if name ∉ names
                info.name = name
                push!(names, info.name)
                break
            else
                name_counter[] += 1
            end
        end
    end
end

verify_support(c) = c

function verify_support(c::Symbol)
    if c !== :NaN # blocks NaN and +/-Inf
        return c
    end
    error("Got NaN in a constraint or objective.")
end

function verify_support(c::Real)
    if isfinite(c) # blocks NaN and +/-Inf
        return c
    end
    error("Expected number but got $c")
end

function verify_support(c::Expr)
    if c.head == :comparison
        map(verify_support, c.args)
        return c
    end
    if c.head == :call
        if c.args[1] in (:+, :-, :*, :/, :exp, :log)
            return c
        elseif c.args[1] in (:<=, :>=, :(==))
            map(verify_support, c.args[2:end])
            return c
        elseif c.args[1] == :^
            @assert isa(c.args[2], Real) || isa(c.args[3], Real)
            return c
        else # TODO: do automatic transformation for x^y, |x|
            error("Unsupported expression $c")
        end
    end
    return c
end

function print_var_definitions(m::BaronModel, fp, header, condition)
    idx = filter(condition, 1 : length(m.variable_info))
    if !isempty(idx)
        println(fp, header, join([m.variable_info[i].name for i in idx], ", "), ";")
    end
end

wrap_with_parens(x::String) = string("(", x, ")")

to_str(x) = string(x)

struct UnrecognizedExpressionException <: Exception
    exprtype::String
    expr
end
function Base.showerror(io::IO, err::UnrecognizedExpressionException)
    print(io, "UnrecognizedExpressionException: ")
    print(io, "unrecognized $(err.exprtype) expression: $(err.expr)")
end

function to_str(c::Expr)
    if c.head == :comparison
        if length(c.args) == 3
            return join([to_str(c.args[1]), c.args[2], c.args[3]], " ")
        elseif length(c.args) == 5
            return join([c.args[1], c.args[2], to_str(c.args[3]),
                         c.args[4], c.args[5]], " ")
        else
            throw(UnrecognizedExpressionException("comparison", c))
        end
    elseif c.head == :call
        if c.args[1] in (:<=,:>=,:(==))
            if length(c.args) == 3
                return join([to_str(c.args[2]), to_str(c.args[1]), to_str(c.args[3])], " ")
            elseif length(c.args) == 5
                return join([to_str(c.args[1]), to_str(c.args[2]), to_str(c.args[3]), to_str(c.args[4]), to_str(c.args[5])], " ")
            else
                throw(UnrecognizedExpressionException("comparison", c))
            end
        elseif c.args[1] in (:+,:-,:*,:/)
            if all(d->isa(d, Real), c.args[2:end]) # handle unary case
                return wrap_with_parens(string(eval(c)))
            elseif c.args[1] == :- && length(c.args) == 2
                return wrap_with_parens(string("(-$(to_str(c.args[2])))"))
            else
                return wrap_with_parens(string(join([to_str(d) for d in c.args[2:end]], string(c.args[1]))))
            end
        elseif c.args[1] == :^
            if length(c.args) != 3
                throw(UnrecognizedExpressionException("function call", c))
            end
            if c.args[3] isa Real
                return wrap_with_parens(string(to_str(c.args[2]), c.args[1], c.args[3]))
            else
                # BARON does not support x^y natively for x,y variables. Instead
                # we transform to the equivalent expression exp(y * log(x)).
                return to_str(:( exp( $(c.args[3]) * log($(c.args[2])) ) ))
            end
        elseif c.args[1] in (:exp,:log)
            if length(c.args) != 2
                throw(UnrecognizedExpressionException("function call", c))
            end
            return wrap_with_parens(string(c.args[1], wrap_with_parens(to_str(c.args[2]))))
        elseif c.args[1] == :abs
            if length(c.args) != 2
                throw(UnrecognizedExpressionException("function call", c))
            end
            # BARON does not support abs(x) natively for variable x. Instead
            # we transform to the equivalent expression sqrt(x^2).
            return to_str(:( ( $(c.args[2])^2.0 )^(0.5) ))
        else
            throw(UnrecognizedExpressionException("function call", c))
        end
    elseif c.head == :ref
        if c.args[1] == :x
            idx = c.args[2]
            @assert isa(idx, Int)
            # TODO decide is use use defined names
            # might be messy becaus a use can call his variable "sin"
            return "x$idx"
        else
            throw(UnrecognizedExpressionException("reference", c))
        end
    end
end

function write_bar_file(m::BaronModel)
    open(m.problem_file_name, "w") do fp
        # First: process any options
        println(fp, "OPTIONS{")
        for (opt,setting) in m.options
            if isa(setting, AbstractString) # wrap it in quotes
                println(fp, "$opt: \"$setting\";")
            else
                println(fp, "$opt: $setting;")
            end
        end
        println(fp, "}")
        println(fp)

        # Ensure that all variables and constraints have a name
        set_unique_names!(m.variable_info, "x")
        set_unique_names!(m.constraint_info, "e")

        # Next, define variables
        print_var_definitions(m, fp, "BINARY_VARIABLES ",   v->(m.variable_info[v].category == :Bin))
        print_var_definitions(m, fp, "INTEGER_VARIABLES ",  v->(m.variable_info[v].category == :Int))
        print_var_definitions(m, fp, "POSITIVE_VARIABLES ", v->(m.variable_info[v].category == :Cont && m.variable_info[v].lower_bound == 0))
        print_var_definitions(m, fp, "VARIABLE ",           v->(m.variable_info[v].category == :Cont && m.variable_info[v].lower_bound != 0))
        println(fp)

        # Print variable bounds
        # TODO:
        # usage of variable_info.name must be revisited
        # now user names are disabled, but if they are enable
        # then bounds will use these names but
        # expressions will just use default names: "x$(variable_inde.value)"
        if any(info -> info.lower_bound !== nothing, m.variable_info)
            println(fp, "LOWER_BOUNDS{")
            for variable_info in m.variable_info
                l = variable_info.lower_bound
                if l !== nothing
                    println(fp, "$(variable_info.name): $l;")
                end
            end
            println(fp, "}")
            println(fp)
        end
        if any(info -> info.upper_bound !== nothing, m.variable_info)
            println(fp, "UPPER_BOUNDS{")
            for variable_info in m.variable_info
                u = variable_info.upper_bound
                if u !== nothing
                    println(fp, "$(variable_info.name): $u;")
                end
            end
            println(fp, "}")
            println(fp)
        end

        # Now let's declare the equations
        if !isempty(m.constraint_info)
            println(fp, "EQUATIONS ", join([constr.name for constr in m.constraint_info], ", "), ";")
            for c in m.constraint_info
                print(fp, c.name, ": ")
                str = to_str(c.expression)
                # print(fp, str)
                if c.lower_bound == c.upper_bound
                    print(fp, str, " == ", c.upper_bound)
                else
                    if c.lower_bound !== nothing && c.upper_bound !== nothing
                        print(fp, c.lower_bound, " <= ", str, " <= ", c.upper_bound)
                    elseif c.lower_bound !== nothing
                        print(fp, str, " >= ", c.lower_bound)
                    elseif c.upper_bound !== nothing
                        print(fp, str, " <= ", c.upper_bound)
                    else
                        error("Unexpected case writing constraint $c.expression.")
                    end
                end
                println(fp, ";")
            end
            println(fp)
        end

        # Now let's do the objective
        print(fp, "OBJ: ")
        if m.objective_sense == :Feasibility || m.objective_expr === nothing
            print(fp, "minimize 0")
        else
            if m.objective_sense == :Min
                print(fp, "minimize ")
            elseif m.objective_sense == :Max
                print(fp, "maximize ")
            else
                error("Objective sense not recognized.")
            end
            val = to_str(m.objective_expr)
            # might get here as: m.objective_expr = :(())
            if val === nothing
                print(fp, "0")
            else
                print(fp, val)
            end
        end
        println(fp, ";")
        println(fp)

        if any(v -> v.start !== nothing, m.variable_info)
            println(fp, "STARTING_POINT{")
            for var in m.variable_info
                if var.start !== nothing
                    println(fp, "$(var.name): $(var.start)")
                end
            end
            println(fp, "}")
        end
    end
end

function read_results(m::BaronModel)
    m.solution_info = SolutionStatus()

    # First, read the time file to get the solution status
    nodeopt = 0
    open(m.times_file_name, "r") do fp
        spl = split(readchomp(fp))
        if m.objective_sense == :Min
            m.solution_info.dual_bound = parse(Float64, spl[6])
            m.solution_info.objective_value = parse(Float64, spl[7])
        else
            m.solution_info.dual_bound = parse(Float64, spl[7])
            m.solution_info.objective_value = parse(Float64, spl[6])
        end
        m.solution_info.solver_status = BaronSolverStatus(parse(Int, spl[8]))
        m.solution_info.model_status = BaronModelStatus(parse(Int, spl[9]))
        nodeopt = parse(Int, spl[12])
        m.solution_info.wall_time = parse(Float64, spl[end])
    end

    # Next, we read the results file to get the solution
    if nodeopt != -3 # parse as long as there exist a solution
        m.solution_info.feasible_point = fill(NaN, length(m.variable_info))
        open(m.result_file_name, "r") do fp
            while true
                startswith(readline(fp), "The best solution found") && break
                eof(fp) && error("Reached end of results file without finding expected optimal primal solution")
            end
            readline(fp)
            readline(fp)

            while true
                line = chomp(readline(fp))
                parts = split(line)
                isempty(parts) && break
                mt = match(r"\d+", parts[1])
                mt == nothing && error("Cannot find appropriate variable index from $(parts[1])")
                v_idx = parse(Int, mt.match)
                v_val = parse(Float64, parts[3])
                m.solution_info.feasible_point[v_idx] = v_val
            end
        end
    end
    return
end


text:
"""
Script in which I'll store examples
"""

using GPARatScale
using Stheno

using Plots
using Suppressor

function plot_optimized_example()
    SAMPLES = 1  # samples we take for the GP representation

    x, y_obs, x_true, y_true = generate_small_dataset()
    println("Creating f1 GP")
    f1_gp_post = create_optim_gp_post(x, y_obs[1], kernel_structure=Matern52();
    i_log_l=3, debug=true)
    f2_gp_post = create_optim_gp_post(x, y_obs[2], kernel_structure=Matern52())
    f3_gp_post = create_optim_gp_post(x, y_obs[3], kernel_structure=EQ())

    println("Creating f1 GPAR")
    f1_gpar_post = create_optim_gpar_post(
        x,
        y_obs[1];
        time_kernel = Matern52(),
        i_log_time_l=-3,
        i_log_time_var=0.2,
        i_log_noise_sigma=-10,
        multi_input = false,
        debug = true,
    )
    f2_gpar_post = create_optim_gpar_post(
        [x, y_obs[1]],
        y_obs[2];
        time_kernel = Matern52(),
        out_kernel = EQ(),
        multi_input = true,
        debug = false,
    )
    f3_gpar_post = create_optim_gpar_post(
        [x, y_obs[1], y_obs[2]],
        y_obs[3];
        time_kernel = Matern52(),
        out_kernel = Matern52(),
        multi_input = true,
        debug = false,
    )

    # Plotting. Supress warnings caused by stheno
    @suppress begin
        gr();
        overall_plot = plot(layout = (3, 1), legend = false);
        # Plot data
        plot!(overall_plot, x_true, y_true, color = :orange, label = "True")
        scatter!(
            overall_plot,
            x,
            y_obs,
            color = :black,
            label = "Observations",
            markersize = 3,
            markeralpha = 0.8,
        )

        # Plot IGP results
        plot!(
            overall_plot[1],
            f1_gp_post(x_true),
            samples = SAMPLES,
            color = :green,
            fillalpha = 0.5,
            linealpha = 1,
        )
        plot!(
            overall_plot[2],
            f2_gp_post(x_true),
            samples = SAMPLES,
            color = :green,
            fillalpha = 0.3,
            linealpha = 1,
        )
        plot!(
            overall_plot[3],
            f3_gp_post(x_true),
            samples = SAMPLES,
            color = :green,
            fillalpha = 0.3,
            linealpha = 1,
        )

        # Plot GPAR
        # First plot is easy
        y1_mean_vals = mean(f1_gpar_post(x_true))
        plot!(overall_plot[1], x_true, y1_mean_vals, color = :blue, linealpha = 1)
        # Use x and mean of y1 to feed into f2_gpar
        x_y1_true = to_ColVecs([x_true, y1_mean_vals])
        y2_mean_vals = mean(f2_gpar_post(x_y1_true))
        plot!(overall_plot[2], x_true, y2_mean_vals, color = :blue, linealpha = 1)

        # Plot f3
        # Again, collect x, y1, and y2 to feed as input location to f3
        x_y1_y2_true = to_ColVecs([x_true, y1_mean_vals, y2_mean_vals])
        y3_mean_vals = mean(f3_gpar_post(x_y1_y2_true))
        plot!(overall_plot[3], x_true, y3_mean_vals, color = :blue, linealpha = 1)

        display(overall_plot)
    end
end

function plot_temporal_gp_inference_example()
    x, y_obs, x_true, y_true = generate_big_dataset()

    lgssm1, f1_out = get_sde_predictions(
        x,
        y_obs[1],
        x_true;
        kernel_structure = Matern52(),
        debug = true,
    )
    lgssm2, f2_out = get_sde_predictions(x, y_obs[2], x_true; debug=true)
    lgssm3, f3_out = get_sde_predictions(x, y_obs[3], x_true; debug=true)
    # Plotting
    gr();
    overall_plot = plot(layout = (3, 1), legend = false);
    # Plot data
    # # Code used to plot the the magnitude of the observation noise
    # scatter!(overall_plot[1], x, y_obs[1], color=:black, alpha=0.1)
    # scatter!(overall_plot[2], x, y_obs[2], color=:black, alpha=0.1)
    # scatter!(overall_plot[3], x, y_obs[3], color=:black, alpha=0.1)

    plot!(overall_plot, x_true, y_true, color = :orange, label = "True")
    function plot_gp_result(plot_ref, out; ylimits=nothing, standard_devs=3)
        posterior_mean = [f.m[1] for f in out]
        std = [f.P[1] for f in out]
        # Plot mean
        plot!(plot_ref, x_true, posterior_mean, color = :green, linealpha = 1)
        # Plot error bars
        plot!(plot_ref, x_true, [posterior_mean posterior_mean];
        linewidth=0.0,
        linecolor=:green,
        fillrange=[posterior_mean .- standard_devs .* std, posterior_mean .+ standard_devs * std],
        fillalpha=0.3,
        fillcolor=:green,
        ylims = ylimits
        );
    end

    # Plot posterior mean of the IGP results
    plot_gp_result(overall_plot[1], f1_out; ylimits=(-5, 3), standard_devs=5)
    plot_gp_result(overall_plot[2], f2_out; standard_devs=5)
    plot_gp_result(overall_plot[3], f3_out; ylimits=(0, 55), standard_devs=5)

    display(overall_plot)

    return lgssm1, lgssm2, lgssm3, f1_out
end

plot_optimized_example()
# plot_temporal_gp_inference_example()

max_stars_repo_name: willtebbutt/GPAR-at-scale

max_issues_repo_name: willtebbutt/GPAR-at-scale

max_forks_repo_name: willtebbutt/GPAR-at-scale


text:
using ITensors
using Test
using Compat

@testset "IndexSet" begin
  idim = 2
  jdim = 3
  kdim = 4
  ldim = 5
  i = Index(idim,"i")
  j = Index(jdim,"j")
  k = Index(kdim,"k")
  l = Index(ldim,"l")
  @testset "show" begin
    indset = IndexSet(i,j,k)
    @test length(sprint(show, indset)) > 1
  end
  @testset "Basic constructors" begin
    I = IndexSet(i,j,k)
    @test IndexSet(I) === I
    @test l ∈ IndexSet(I..., l)
    @test l ∈ IndexSet(l, I...)
    @test length(IndexSet{2}(i,j)) == 2
  end
  @testset "length of IndexSet and friends" begin
    @test length(typeof(IndexSet{2}(i,j))) == 2
    @test length(IndexSet(i,j)) == 2
    @test length(typeof(IndexSet(i,j))) == 2
  end
  @testset "Convert to Index" begin
    @test Index(IndexSet(i)) === i
    @test_throws ErrorException Index(IndexSet(i, j))
  end
  @testset "Index dimensions" begin
    I = IndexSet(i,j,k)
    @test dim(I) == idim*jdim*kdim
    @test dims(I) == (idim,jdim,kdim)
    @test dim(I,1) == idim
    @test dim(I,2) == jdim
    @test dim(I,3) == kdim

    @test maxdim(I) == max(idim,jdim,kdim)
  end

  @testset "Set operations" begin
    I1 = IndexSet(i, j, k)
    I2 = IndexSet(k, l)
    I3 = IndexSet(j, l)
    @test hassameinds(I1, (k, j, i))
    @test firstsetdiff(I1, I2, I3) == i
    @test isnothing(firstsetdiff(I1, IndexSet(k, j, i)))
    @test setdiff(I1, I2) == [i, j]
    @test setdiff(Order(2), I1, I2) == IndexSet(i, j)
    @test hassameinds(setdiff(I1, I2), IndexSet(i, j))
    @test hassameinds(setdiff(Order(2), I1, I2), IndexSet(i, j))
    @test hassameinds(setdiff(I1, I2), (j, i))
    @test hassameinds(setdiff(Order(2), I1, I2), (j, i))
    @test I1 ∩ I2 == [k]
    @test hassameinds(I1 ∩ I2, IndexSet(k))
    @test firstintersect(I1, I2) == k
    @test isnothing(firstintersect(I1, IndexSet(l)))
    @test intersect(I1, IndexSet(j, l)) == [j]
    @test intersect(Order(1), I1, IndexSet(j, l)) == IndexSet(j)
    @test hassameinds(intersect(I1, IndexSet(j, l)), IndexSet(j))
    @test hassameinds(intersect(Order(1), I1, IndexSet(j, l)), IndexSet(j))
    @test firstintersect(I1, IndexSet(j, l)) == j
    @test intersect(I1, IndexSet(j, k)) == [j, k]
    @test intersect(Order(2), I1, IndexSet(j, k)) == IndexSet(j, k)
    @test hassameinds(intersect(I1, (j, k)), IndexSet(j, k))
    @test hassameinds(intersect(I1, (j, k, l)), (j, k))
    @test filter(I1, "i") == IndexSet(i)
    @test filter(I1; tags = "i") == IndexSet(i)
    @test filter(I1; inds = j) == IndexSet(j)
    @test filter(I1; tags = "i", inds = j) == IndexSet()
    @test filter(I1; plev = 1, inds = j) == IndexSet()
    @test filter(I1; plev = 0, inds = k) == IndexSet(k)
    @test filter(I1; plev = 0) == IndexSet(i, j, k)
    @test filter(I1; inds = l) == IndexSet()
    @test hassameinds(filter(I1, "i"), IndexSet(i))
    @test getfirst(I1, "j") == j
    @test isnothing(getfirst(I1, "l"))
    @test findfirst(I1, i) == 1
    @test findfirst(I1, j) == 2
    @test findfirst(I1, k) == 3
    @test isnothing(findfirst(I1, Index(2)))
  end

  @testset "Set operations with Order" begin
    i,j,k,l = Index.(2, ("i", "j", "k", "l"))

    Iij = IndexSet(i, j)
    Ijl = IndexSet(j, l)
    Ikl = IndexSet(k, l)
    Iijk = IndexSet(i, j, k)

    #
    # setdiff
    # intersect
    # symdiff
    # union
    # filter
    #

    #
    # setdiff
    #

    @test setdiff(Iijk, Ikl) == [i, j]
    @test setdiff(Order(2), Iijk, Ikl) == IndexSet(i, j)

    @test setdiff(Iij, Iijk) == Index{Int}[]
    @test setdiff(Order(0), Iij, Iijk) == IndexSet()

    @test setdiff(Iijk, Ikl; tags = "i") == [i]
    @test setdiff(Order(1), Iijk, Ikl; tags = "i") == IndexSet(i)

    @test setdiff(Iijk, Ikl; tags = not("i")) == [j]
    @test setdiff(Order(1), Iijk, Ikl; tags = not("i")) == IndexSet(j)

    @test setdiff(Iijk, Ijl, Ikl) == [i]
    @test setdiff(Order(1), Iijk, Ijl, Ikl) == IndexSet(i)

    #
    # intersect
    #

    @test intersect(Iijk, Ikl) == [k]
    @test intersect(Order(1), Iijk, Ikl) == IndexSet(k)

    @test intersect(Iijk, Iij) == [i, j]
    @test intersect(Order(2), Iijk, Iij) == IndexSet(i, j)

    @test intersect(Iijk, Iij; tags = "i") == [i]
    @test intersect(Order(1), Iijk, Iij; tags = "i") == IndexSet(i)

    #
    # symdiff
    #

    @test symdiff(Iijk, Ikl) == [i, j, l]
    #@test symdiff(Order(3), Iijk, Ikl) == IndexSet(i, j, l)

    @test symdiff(Iijk, Iij) == [k]
    #@test symdiff(Order(3), Iijk, Iij) == IndexSet(i, j, k)

    #@test symdiff(Iijk, Iij; tags = "i") == [i]
    #@test symdiff(Order(1), Iijk, Iij; tags = "i") == IndexSet(i)

    #
    # union
    #

    @test union(Iijk, Ikl) == [i, j, k, l]
    #@test union(Order(4), Iijk, Ikl) == IndexSet(i, j, k, l)

    @test union(Iijk, Iij) == [i, j, k]
    #@test union(Order(3), Iijk, Iij) == IndexSet(i, j, k)

    #@test union(Iijk, Iij; tags = "i") == [i]
    #@test union(Order(1), Iijk, Iij; tags = "i") == IndexSet(i)
  end

  @testset "intersect index ordering" begin
    I = IndexSet(i,k,j)
    J = IndexSet(j,l,i)
    # Test that intersect respects the ordering
    # of the indices in the first IndexSet
    @test hassameinds(intersect(I,J),IndexSet(i,j))
    @test hassameinds(intersect(J,I),IndexSet(j,i))
  end
  @testset "adjoint" begin
    I = IndexSet(i,k,j)
    @test adjoint(I) == IndexSet(i', k', j')
  end
  @testset "mapprime" begin
    I = IndexSet(i',k'',j)
    @test mapprime(I,1,5) == IndexSet(i^5,k'',j)
    @test mapprime(I,2,0) == IndexSet(i',k,j)

    J = IndexSet(i,j,k')
    @test mapprime(J,0,2) == IndexSet(i'',j'',k')

    J = mapprime(J,1,5)
    @test J == IndexSet(i,j,k^5)
  end
  @testset "strides" begin
    I = IndexSet(i, j)
    @test NDTensors.strides(I) == (1, idim)
    @test NDTensors.stride(I, 1) == 1
    @test NDTensors.stride(I, 2) == idim
  end
  @testset "setprime" begin
    I = IndexSet(i, j)
    J = setprime(I, 2, i)
    @test i'' ∈ J
  end
  @testset "prime" begin
    I = IndexSet(i, j)
    J = prime(I, j)
    @test i ∈ J
    @test j' ∈ J
    J = prime(I; inds = j)
    @test i ∈ J
    @test j' ∈ J
    J = prime(I; inds = not(j))
    @test i' ∈ J
    @test j ∈ J
  end
  @testset "noprime" begin
    I = IndexSet(i'', j')
    J = noprime(I)
    @test i ∈ J
    @test j ∈ J
  end
  @testset "swapprime" begin
    I = IndexSet(i,j)
    @test swapprime(I,0,1) == IndexSet(i',j')
    @test swapprime(I,0,4) == IndexSet(i^4,j^4)
    I = IndexSet(i,j'')
    @test swapprime(I,2,0) == IndexSet(i'',j)
    I = IndexSet(i,j'',k,l)
    @test swapprime(I,2,0) == IndexSet(i'',j,k'',l'')
    I = IndexSet(i,k'',j'')
    @test swapprime(I,2,1) == IndexSet(i,k',j')
    # In-place version:
    I = IndexSet(i,k'',j''')
    I = swapprime(I,2,0)
    @test I == IndexSet(i'',k,j''')
    # With tags specified:
    I = IndexSet(i,k,j)
    @test swapprime(I,0,1,"i") == IndexSet(i',k,j)
    @test swapprime(I,0,1,"j") == IndexSet(i,k,j')

    I = IndexSet(i,i',j)
    @test swapprime(I,0,1,"i") == IndexSet(i',i,j)
    @test swapprime(I,0,1,"j") == IndexSet(i,i',j')
  end
  @testset "swaptags" begin
    i1 = Index(2,"Site,A")
    i2 = Index(2,"Site,B")
    is = IndexSet(i1,i2)
    sis = swaptags(is,"Site","Link")
    for j in sis
      @test !hastags(j,"Site")
      @test hastags(j,"Link")
    end
  end
  @testset "broadcasting" begin
    I = IndexSet(i, j)
    J = prime.(I)
    @test J isa IndexSet
    @test i' ∈ J
    @test j' ∈ J
  end
end

nothing


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    nrows = 2
    ncols = 2

    maingl = GridLayout(scene, nrows, ncols, alignmode = Outside(30))

    axes = [
        maingl[i, j] = LAxis(
            scene,
            aspect=DataAspect(),  # make axis the same aspect as the data (images will not be stretched)
            xautolimitmargin=(0, 0), # don't create margins so that the image aligns with the axis spines
            yautolimitmargin=(0, 0)) # same in y
        for i in 1:nrows, j in 1:ncols]

    for i in 1:nrows, j in 1:ncols
        img = rand(rand(300:700), rand(300:700))
        image!(axes[i, j], img)
    end

end


text:
doc"A single input/output pair"
type DataPoint
  x::VecF
  y::VecF
end

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

# goal: some simple wrapers around datapoint that perform transformations
# we want to: select subsets of values, and transform them (including nop)

doc"""
`Transformation`s wrap an index within a data vector, and descibe how to
calculate Tⱼ(x) by transforming the iᵗʰ index.
"""
abstract Transformation

immutable IdentityTransform <: Transformation idx::Int end
transform(t::IdentityTransform, x::AVec) = x[t.idx]
function transform!(t::IdentityTransform, x_transformed::AVec, x::AVec, i::Int)
  x_transformed[i] = x[t.idx]
end

immutable AbsTransform <: Transformation idx::Int end
transform(t::AbsTransform, x::AVec) = abs(x[t.idx])
function transform!(t::AbsTransform, x_transformed::AVec, x::AVec, i::Int)
  x_transformed[i] = abs(x[t.idx])
end

immutable LogPlus1Transform <: Transformation idx::Int end
transform(t::LogPlus1Transform, x::AVec) = (xi = x[t.idx]; xi > 0.0 ? log(xi+1) : 0.0)
function transform!(t::LogPlus1Transform, x_transformed::AVec, x::AVec, i::Int)
  xi = x[t.idx]
  x_transformed[i] = xi > 0.0 ? log(xi + 1.0) : 0.0
end

immutable SquareTransform <: Transformation idx::Int end
transform(t::SquareTransform, x::AVec) = x[t.idx]^2
function transform!(t::SquareTransform, x_transformed::AVec, x::AVec, i::Int)
  x_transformed[i] = x[t.idx]^2
end

immutable CubeTransform <: Transformation idx::Int end
transform(t::CubeTransform, x::AVec) = x[t.idx]^3
function transform!(t::CubeTransform, x_transformed::AVec, x::AVec, i::Int)
  x_transformed[i] = x[t.idx]^3
end

immutable SignSquareTransform <: Transformation idx::Int end
transform(t::SignSquareTransform, x::AVec) = (xi = x[t.idx]; sign(xi) * xi^2)
function transform!(t::SignSquareTransform, x_transformed::AVec, x::AVec, i::Int)
  xi = x[t.idx]
  x_transformed[i] = sign(xi) * xi^2
end

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

abstract Transformer

immutable IdentityTransformer <: Transformer end
transform(transformer::IdentityTransformer, x) = x
transform!(transformer::IdentityTransformer, x_transformed::AVec, x::AVec) = copy!(x_transformed, x)

immutable VectorTransformer <: Transformer
  transformations::Vector{Transformation}
end
transform{T}(transformer::VectorTransformer, x::AVec{T}) = T[transform(t,x) for t in transformer.transformations]

function transform!(transformer::VectorTransformer, x_transformed::AVec, x::AVec)
  for i in 1:length(x_transformed)
    transform!(transformer.transformations[i], x_transformed, x, i)
  end
end

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

doc"A list of input/output data points"
type DataPoints <: AbstractVector{DataPoint}
  data::Vector{DataPoint}
end

function DataPoints(x::AMatF, y::AMatF, indices::AVecI = 1:nrows(x))
  DataPoints([DataPoint(vec(x[i,:]), vec(y[i,:])) for i in indices])
end

function DataPoints(x::AMat, y::AVec, indices::AVecI = 1:nrows(x))
  DataPoints([DataPoint(float(vec(x[i,:])), float(vec(y[i,:]))) for i in indices])
end

function unzip(dps::DataPoints)
  n = length(dps)
  mx = length(first(dps).x)
  my = length(first(dps).y)
  x = zeros(n, mx)
  y = zeros(n, my)
  for (i,dp) in enumerate(dps)
    x[i,:] = dp.x
    y[i,:] = dp.y
  end
  x, y
end

Base.getindex(dps::DataPoints, i::Int) = dps.data[i]
Base.getindex(dps::DataPoints, a::AVecI) = DataPoints(dps.data[a])
Base.setindex!(dps::DataPoints, dp::DataPoint, i::Int) = (dps.data[i] = dp)
Base.push!(dps::DataPoints, dp::DataPoint) = push!(dps.data, dp)
Base.append!(dps::DataPoints, dps2::DataPoints) = append!(dps.data, dps2.data)
Base.length(dps::DataPoints) = length(dps.data)
Base.size(dps::DataPoints) = size(dps.data)

function splitDataPoints(dps::DataPoints, pct::Real)
  r1, r2 = splitRange(length(dps), pct)
  DataPoints[r1], DataPoints[r2]
end

StatsBase.sample(dps::DataPoints) = dps[sample(1:length(dps))]
StatsBase.sample(dps::DataPoints, n::Int) = dps[sample(1:length(dps), n)]

Base.shuffle(dps::DataPoints) = DataPoints(shuffle(dps.data))
Base.shuffle!(dps::DataPoints) = shuffle!(dps.data)

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

doc"""
Generic approach to sampling data.  Allows for many different approaches in a unified framework.
DataSamplers should implement the StatsBase.sample method which returns a DataPoint object,
and a DataPoints constructor which returns a DataPoints object with unique DataPoints.
"""
abstract DataSampler

# return a DataPoints object with n samples from the sampler
StatsBase.sample(sampler::DataSampler, n::Int) = DataPoints(DataPoint[sample(sampler) for i in 1:n])

doc"Sample the whole data set"
immutable SimpleSampler <: DataSampler
  data::DataPoints
end
Base.print(io::IO, s::SimpleSampler) = print(io, "SimpleSampler{n=$(length(s.data))}")
Base.show(io::IO, s::SimpleSampler) = print(io,s)

StatsBase.sample(sampler::SimpleSampler) = sample(sampler.data)
DataPoints(sampler::SimpleSampler) = sampler.data

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

doc"Sample from a subset of the data (defined by the range)"
immutable SubsetSampler <: DataSampler
  data::DataPoints
  range::AVecI
end
Base.print(io::IO, s::SubsetSampler) = print(io, "SubsetSampler{n=$(length(s.range))}")
Base.show(io::IO, s::SubsetSampler) = print(io,s)

# sample from the range
StatsBase.sample(sampler::SubsetSampler) = sampler.data[sample(sampler.range)]

DataPoints(sampler::SubsetSampler) = sampler.data[sampler.range]

doc"Create two `SubsetSampler`s, one for the first `pct` of the data and the other for the remaining `1-pct`."
function splitDataSamplers(data::DataPoints, pct::Real)
  r1, r2 = splitRange(length(data), pct)
  SubsetSampler(data, r1), SubsetSampler(data, r2)
end

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

"separate a dataset into one DataPoint for each y value"
type StratifiedSampler <: DataSampler
  data::DataPoints
  idxlist::Vector{VecI}
  n::Int  # current bucket
end
Base.print(io::IO, s::StratifiedSampler) = print(io, "StratifiedSampler{n=$(length(s.data)), numbuckets=$(length(s.idxlist)), nextbucket=$(s.n)}")
Base.show(io::IO, s::StratifiedSampler) = print(io,s)

function StratifiedSampler(dps::DataPoints)
  idxlist = VecI[]
  for (i,dp) in enumerate(dps)

    matched = false
    for indices in idxlist
      if dp.y == dps[first(indices)].y
        push!(indices, i)
        matched = true
        break
      end
    end

    if !matched
      push!(idxlist, [i])
    end
  end

  StratifiedSampler(dps, idxlist, 1)
end

# sample from each y value equally
function StatsBase.sample(sampler::StratifiedSampler)

  # get the right bucket
  indices = sampler.idxlist[sampler.n]

  # update the bucket number
  sampler.n += 1
  if sampler.n > length(sampler.idxlist)
    sampler.n = 1 # reset
  end

  # do the sample
  sampler.data[sample(indices)]
end

DataPoints(sampler::StratifiedSampler) = sampler.data

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

type CrossValidationIterator
  dps::DataPoints
  shuffledIndices::VecI
  validrngs::Vector{UnitRange}
end

function makeCVPair(cviter::CrossValidationIterator, i::Int)
  vrng = cviter.validrngs[i]
  trng = setdiff(1:length(cviter.dps), vrng)
  map(rng -> SubsetSampler(cviter.dps, cviter.shuffledIndices[rng]), (trng, vrng))
end

Base.start(cviter::CrossValidationIterator) = 1
Base.done(cviter::CrossValidationIterator, state::Int) = state > length(cviter.validrngs)
Base.next(cviter::CrossValidationIterator, state::Int) = (makeCVPair(cviter,state), state+1)

# returns an iterator returning (trainsampler, validationsampler) tuples for each fold
function crossValidationSets(dps::DataPoints, numFolds::Int)
  n = length(dps)
  shuffledIndices = shuffle(collect(1:n))
  a = round(Int, collect(linspace(0,n,numFolds+1)))
  CrossValidationIterator(dps, shuffledIndices, [a[i]+1:a[i+1] for i in 1:numFolds])
end

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


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const inputText = "Man is distinguished, not only by his reason, but by this singular passion from other animals, which is a lust of the mind, that by a perseverance of delight in the continued and indefatigable generation of knowledge, exceeds the short vehemence of any carnal pleasure."
const encodedMaxLine76 =
"""TWFuIGlzIGRpc3Rpbmd1aXNoZWQsIG5vdCBvbmx5IGJ5IGhpcyByZWFzb24sIGJ1dCBieSB0aGlz
IHNpbmd1bGFyIHBhc3Npb24gZnJvbSBvdGhlciBhbmltYWxzLCB3aGljaCBpcyBhIGx1c3Qgb2Yg
dGhlIG1pbmQsIHRoYXQgYnkgYSBwZXJzZXZlcmFuY2Ugb2YgZGVsaWdodCBpbiB0aGUgY29udGlu
dWVkIGFuZCBpbmRlZmF0aWdhYmxlIGdlbmVyYXRpb24gb2Yga25vd2xlZGdlLCBleGNlZWRzIHRo
ZSBzaG9ydCB2ZWhlbWVuY2Ugb2YgYW55IGNhcm5hbCBwbGVhc3VyZS4="""

# Encode and decode
fname = tempname()
open(fname, "w") do f
    opipe = Base64EncodePipe(f)
    write(opipe,inputText)
    @test close(opipe) === nothing
end

open(fname, "r") do f
    ipipe = Base64DecodePipe(f)
    @test read(ipipe, String) == inputText
    @test close(ipipe) === nothing
end
rm(fname)

# Encode to string and decode
@test String(base64decode(base64encode(inputText))) == inputText

# Decode with max line chars = 76 and padding
ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76))
@test read(ipipe, String) == inputText

# Decode with max line chars = 76 and no padding
ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76[1:end-1]))
@test read(ipipe, String) == inputText

# Decode with two padding characters ("==")
ipipe = Base64DecodePipe(IOBuffer(string(encodedMaxLine76[1:end-2],"==")))
@test read(ipipe, String) == inputText[1:end-1]

# Test incorrect format
ipipe = Base64DecodePipe(IOBuffer(encodedMaxLine76[1:end-3]))
@test_throws ArgumentError read(ipipe, String)

# issue #21314
@test base64decode(chomp("test")) == base64decode("test")

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
using Printf

println("\n")
println("------------------------------------------")
println("--  重複しないように                     --")
println("--  データセットを用意すること           --")
println("------------------------------------------")
println("\n")

function main()
    # 座標点数 ( 格子数 = nr-1 )
    nr    = 200
    nwing = 200

    #=
    fout     = "xy_NACA0012"         # output file name
    NACAname = "NACA0012"            # NACA0012 type
    r    = 2.5                       # O型格子の半径
    cx   = 0.5                       # 翼データはx=0~1のため、その間の0.5をO型の中心とする。
    min_dy = 1.0e-3                  # 半径方向の最小格子幅, m
    dnum = 4                         # NACA0012 type = 4, NACA64A010 type = 6
    =#

    fout     = "xy_hayabusa"         # output file name
    NACAname = "hayabusa_shape"      # NACA0012 type
    r      = 5.0                     # O型格子の半径
    cx     = -0.101                  # O型格子の中心
    min_dy = 1.0e-2                  # 半径方向の最小格子幅, m
    dnum   = 4                       # NACA0012 type = 4, NACA64A010 type = 6

    #=
    fout     = "xy_NACA64A010"       # output file name
    NACAname = "NACA64A010"          # NACA64A010 type
    r    = 2.5                       # O型格子の半径
    cx   = 0.5                       # 翼データはx=0~1のため、その間の0.5をO型の中心とする。
    dnum = 6                         # NACA0012 type = 4, NACA64A010 type = 6
    =#

    #=
    fout     = "xy_NACA4412"       # output file name
    NACAname = "NACA4412"          # NACA64A010 type
    r    = 2.5                       # O型格子の半径
    cx   = 0.5                       # 翼データはx=0~1のため、その間の0.5をO型の中心とする。
    dnum = 6                         # NACA0012 type = 4, NACA64A010 type = 6
    =#

    println("\n")
    println("------------------------------------------")
    println("--  格子の切れ目は下記のデータを使用する --")
    println("--  4d : 下面の最後のデータ             --")
    println("--  6d : 一番最初のデータ               --")
    println("------------------------------------------")
    println("\n")

    if dnum == 4
        surface_len, edgex, edgey, aroundx, aroundy = read_NACA_4d(NACAname)
    elseif dnum == 6
        surface_len, edgex, edgey, aroundx, aroundy = read_NACA_6d(NACAname)
    else
        println("\n")
        println("--  データ読み込みの方法を決めてください --")
        println("--  datanum = 4 or 6                   --")
        println("\n")
        throw(UndefVarError(:x))
    end

    x1, y1 = cal_area1(edgex, edgey, r, cx, nr)
    x2, y2 = cal_area2(surface_len, aroundx, aroundy,nwing)
    x3, y3 = cal_area1(edgex, edgey, r, cx, nr)
    x4, y4 = cal_area4(edgex, edgey, r, cx, nwing)

    x,y = algebraic_grid_gene(nr,nwing,x1,x2,x3,x4,y1,y2,y3,y4)

    # 中央寄せ
    tx = zeros(nr)
    ty = zeros(nr)
    for i in 1:nwing
        for j in 1:nr
            tx[j] = x[i,j]
            ty[j] = y[i,j]
        end
        tx, ty = draw_in_grid(tx, ty, nr, min_dy)
        for j in 1:nr
            x[i,j] = tx[j]
            y[i,j] = ty[j]
        end
    end

    open(fout,"w") do f
        a1 = @sprintf("%8.8e", nwing)
        a2 = @sprintf("%8.8e", nr)
        write(f, a1*" "*a2*"\n")
        for j in 1:nr
            for i in 1:nwing
                ai = string(i)
                aj = string(j)
                a1 = @sprintf("%8.8e", x[i,j])
                a2 = @sprintf("%8.8e", y[i,j])
                write(f, ai*" "*aj*" "*a1*" "*a2*"\n")
            end
        end
    end

    open("External shape","w") do f
        for j in 1:nwing
            a1 = @sprintf("%8.8e", x2[j])
            a2 = @sprintf("%8.8e", y2[j])
            write(f, a1*" "*a2*"\n")
        end
    end
end

function cal_area1(edgex, edgey, r, cx, nr)
    x1 = zeros(nr)
    y1 = zeros(nr)

    dx  = (r - edgex + cx) / (nr-1)

    for i in 1:nr
        x1[i] = edgex + (i-1)*dx
        y1[i] = edgey
    end

    return x1,y1
end

function cal_area2(surface_len, aroundx, aroundy, nwin)
    x2 = zeros(nwin)
    y2 = zeros(nwin)

    dx  = surface_len/(nwin-1)
    naround = length(aroundx)

    pointlen  = 1e-50
    aroundlen = 0.0
    l = 0

    x2[1] = aroundx[1]
    y2[1] = aroundy[1]
    for i in 2:nwin
        pointlen  += dx
        while (pointlen - aroundlen)/(aroundlen) >= 1e-5
            l += 1
            aroundlen += ((aroundx[l]-aroundx[l+1])^2 + (aroundy[l]-aroundy[l+1])^2)^0.5
        end

        p1x = aroundx[l]
        p1y = aroundy[l]
        p2x = aroundx[l+1]
        p2y = aroundy[l+1]
        dis_2top = aroundlen - pointlen
        dis_2to1 = ((p1x-p2x)^2 + (p1y-p2y)^2)^0.5

        x2[i] = dis_2top/dis_2to1*p1x + (dis_2to1-dis_2top)/dis_2to1*p2x
        y2[i] = dis_2top/dis_2to1*p1y + (dis_2to1-dis_2top)/dis_2to1*p2y
    end

    return x2,y2
end

function cal_area4(edgex, edgey, r, cx, nwin)
    x4 = zeros(nwin)
    y4 = zeros(nwin)

    len = 2*pi
    dt  = len/(nwin-1)

    for i in 1:nwin
        t     = len - (i-1)*dt
        x4[i] = r * cos(t) + cx
        y4[i] = r * sin(t)
    end
    return x4,y4
end

#= ----------------------------------------------
 -- draw_in_grid(x, y, n, min_d)
 中央に格子を寄せる
 拡大率alphaとして、隣の格子との比が1:alphaになるように計算する
 半径方向の長さl, 壁面第一格子幅dx, 拡大率alphaの関係は次の等比数列で表すことができる。

 l = dx + dx*alpha + + dx*alpha^2 + ... + + dx*alpha^n

 この関係から、ニュートン法によりalphaを求め、格子を作り直す
 # ----------------------------------------------
=#

function draw_in_grid(x, y, n, min_d)

    l = 0.0         # 直線の長さ
    for i in 2:n
        l += ((x[i] - x[i-1])^2 + (y[i] - y[i-1])^2)^0.5
    end

    if l/(n-1) < min_d
        println("------------------------\n")
        println(" error \n ")
        println(" please change min dy \n ")
        println(" dy = ")
        println(l/(nr-1))
        println(" \n ")
        println(" min dy = ")
        println(min_d)
        println("------------------------")
    end

    alpha = newton_method(n-2, l, min_d)
    #=
    println(" dy = ")
    println(l/(nr-1))
    println(" \n ")
    println(" min dy = ")
    println(min_dy)
    println(" \n ")
    println(" alpha = ")
    println(alpha)
    =#

    pointlen  = 1e-50   # 経路に沿ったその点までの距離
    aroundlen = 1e-50   # pointlenより1点分長くなるようにとった点までの距離
    ite = 0             # 点のiteration
    newx = copy(x)
    newy = copy(y)
    for i in 2:n
        dx = min_d * alpha^(i-2)
        pointlen += dx
        while (pointlen - aroundlen)/(aroundlen) >= 1e-5
            ite += 1
            aroundlen += ((x[ite]-x[ite+1])^2 + (y[ite]-y[ite+1])^2)^0.5
        end

        p1x = x[ite]
        p1y = y[ite]
        p2x = x[ite+1]
        p2y = y[ite+1]
        dis_2top = aroundlen - pointlen
        dis_2to1 = ((p1x-p2x)^2 + (p1y-p2y)^2)^0.5

        newx[i] = dis_2top/dis_2to1*p1x + (dis_2to1-dis_2top)/dis_2to1*p2x
        newy[i] = dis_2top/dis_2to1*p1y + (dis_2to1-dis_2top)/dis_2to1*p2y
    end

    return newx, newy
end

function newton_method(n,l,dx)
	# ニュートン法
	alpha = 1.5
	for i in 1:10000
		old = alpha

		alpha = alpha - f(alpha,n,dx,l)/f_dash(alpha,n,dx,l)

		#if abs(alpha-old) / alpha < 10^(-4)
			#break
		#end
	end
	return alpha
end

function f(alpha,n,dx,l)
    f = 1-alpha^(n+1) - (1-alpha) *l/dx
    return f
end

function f_dash(alpha,n,dx,l)
    f = -(n+1)*alpha^(n) + l/dx
    return f
end

# ----------------------------------------------
# -- algebraic_grid_gene(nx,ny,x1,x2,x3,x4,y1,y2,y3,y4)
# -- 幾何学的計算により、境界適合格子を作成
# ----------------------------------------------
function algebraic_grid_gene(nr,nwing,x1,x2,x3,x4,y1,y2,y3,y4)
    x = zeros(nwing,nr)
    y = zeros(nwing,nr)

    for j in 2:nwing-1
        for i in 2:nr-1
            xi = (i-1)/(nr-1)
            eta  = (j-1)/(nwing-1)
            x[j,i] = (1-eta)*x1[i] + eta*x3[i] + (1-xi)*x2[j] + xi*x4[j] -
                    (xi*eta*x3[nr] + xi*(1-eta)*x1[nr] + eta*(1-xi)*x3[1] + (1-xi)*(1-eta)*x1[1])
            y[j,i] = (1-eta)*y1[i] + eta*y3[i] + (1-xi)*y2[j] + xi*y4[j] -
                    (xi*eta*y3[nr] + xi*(1-eta)*y1[nr] + eta*(1-xi)*y3[1] + (1-xi)*(1-eta)*y1[1])
        end
    end

    # Substitution of bd
    for i in 1:nwing
        x[i,1]  = x2[i]
        y[i,1]  = y2[i]
        x[i,nr] = x4[i]
        y[i,nr] = y4[i]
    end
    for j in 1:nr
        x[1,j]  = x1[j]
        y[1,j]  = y1[j]
        x[nwing,j] = x3[j]
        y[nwing,j] = y3[j]
    end

    return x,y
end

# ----------------------------------------------
# -- read_NACA_4d(NACAname)
# -- データ読み込み、(NACA0012のデータセットの読み方)
# ----------------------------------------------
function read_NACA_4d(NACAname)

    fff=[]
    open(NACAname, "r") do f                  # 全行格納
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)     # 改行分割(\n)

    for i in 1:length(fff)
        fff[i] = replace(fff[i]," \r" => "")  # 改行削除(\r)
    end

    # 上面、下面データ数の読み込み
    temp   = split(fff[2]," ")
    filter!(e->e≠"",temp)

    nt = Int(parse(Float64,temp[1]))
    nb = Int(parse(Float64,temp[2]))

    # 上面、下面ポイントの読み込み
    top = zeros(nt, 2)
    bot = zeros(nb, 2)

    skip = 2
    for i in 1:nt
        temp = split(fff[i+skip]," ")
        filter!(e->e≠"",temp)

        top[i,1] = parse(Float64,temp[1])
        top[i,2] = parse(Float64,temp[2])
    end

    skip = 2 + nt
    for i in 1:nb
        temp = split(fff[i+skip]," ")
        filter!(e->e≠"",temp)

        bot[i,1] = parse(Float64,temp[1])
        bot[i,2] = parse(Float64,temp[2])
    end

    # 表面データの結合
    aroundx = zeros(nt+nb+1)
    aroundy = zeros(nt+nb+1)
    for i in 1:nb
        aroundx[i] = bot[nb+1-i,1]
        aroundy[i] = bot[nb+1-i,2]
    end
    for i in 1:nt
        aroundx[i+nb] = top[i,1]
        aroundy[i+nb] = top[i,2]
    end
    aroundx[nt+nb+1] = bot[nb,1]
    aroundy[nt+nb+1] = bot[nb,2]

    # 表面線長さ、エッジ座標の登録
    surface_len = 0.0
    edgex = 0.0
    edgey = 0.0

    n = nt+nb
    temp = 0.0
    for i in 1:n
        temp += ((aroundx[i]-aroundx[i+1])^2 + (aroundy[i]-aroundy[i+1])^2)^0.5
    end
    surface_len = temp

    edgex = aroundx[1]
    edgey = aroundy[1]

    return surface_len, edgex, edgey, aroundx, aroundy
end

# ----------------------------------------------
# -- read_NACA_6d(NACAname)
# -- データ読み込み、(NACA64A010のデータセットの読み方)
# ----------------------------------------------
function read_NACA_6d(NACAname)

    fff=[]
    open(NACAname, "r") do f                  # 全行格納
        fff = read(f,String)
    end
    fff = split(fff,"\n",keepempty=false)     # 改行分割(\n)

    for i in 1:length(fff)
        fff[i] = replace(fff[i]," \r" => "")  # 改行削除(\r)
    end

    # 表面データの読み込み
    skip = 1
    n    = length(fff)-skip
    aroundx = zeros(n+1)
    aroundy = zeros(n+1)
    taroundx = zeros(n+1)
    taroundy = zeros(n+1)
    for i in 1:n
        temp = split(fff[i+skip]," ")
        filter!(e->e≠"",temp)

        taroundx[i] = parse(Float64,temp[1])
        taroundy[i] = parse(Float64,temp[2])
    end
    taroundx[n+1] = taroundx[1]
    taroundy[n+1] = taroundy[1]

    for i in 1:n+1
        aroundx[i] = taroundx[n+1-(i-1)]
        aroundy[i] = taroundy[n+1-(i-1)]
    end

    # 表面線長さ、エッジ座標の登録
    surface_len = 0.0
    edgex = 0.0
    edgey = 0.0

    temp = 0.0
    for i in 1:n
        temp += ((aroundx[i]-aroundx[i+1])^2 + (aroundy[i]-aroundy[i+1])^2)^0.5
    end
    surface_len = temp

    edgex = aroundx[1]
    edgey = aroundy[1]

    return surface_len, edgex, edgey, aroundx, aroundy
end

# ---------------------------------
main()


text:
"""
$(TYPEDSIGNATURES)

Generate `NonlinearSystem`, dependent variables, and parameters from an `NonlinearProblem`.
"""
function modelingtoolkitize(prob::NonlinearProblem; kwargs...)
    p = prob.p
    has_p = !(p isa Union{DiffEqBase.NullParameters, Nothing})

    _vars = reshape([variable(:x, i) for i in eachindex(prob.u0)], size(prob.u0))

    vars = prob.u0 isa Number ? _vars : ArrayInterfaceCore.restructure(prob.u0, _vars)
    params = if has_p
        _params = define_params(p)
        p isa Number ? _params[1] :
        (p isa Tuple || p isa NamedTuple ? _params :
         ArrayInterfaceCore.restructure(p, _params))
    else
        []
    end

    if DiffEqBase.isinplace(prob)
        rhs = ArrayInterfaceCore.restructure(prob.u0, similar(vars, Num))
        prob.f(rhs, vars, params)
    else
        rhs = prob.f(vars, params)
    end
    out_def = prob.f(prob.u0, prob.p)
    eqs = vcat([0.0 ~ rhs[i] for i in 1:length(out_def)]...)

    sts = vec(collect(vars))

    params = if params isa Number || (params isa Array && ndims(params) == 0)
        [params[1]]
    else
        vec(collect(params))
    end
    default_u0 = Dict(sts .=> vec(collect(prob.u0)))
    default_p = has_p ? Dict(params .=> vec(collect(prob.p))) : Dict()

    de = NonlinearSystem(eqs, sts, params,
                         defaults = merge(default_u0, default_p);
                         name = gensym(:MTKizedNonlinProb),
                         kwargs...)

    de
end


text:
function read_treebank_grammar(number_sentences::Int, Score=LogProb)
    @assert number_sentences in (1, 5, 50, 500, 972)
    grammarfile = joinpath(dirname(@__FILE__), "WSJ_data", "WSJ.$number_sentences.grammar.txt")
    Grammar(readstring(grammarfile), ["START"], Score)
end

function compute_treebank_scores(number_sentences::Int, grammar=read_treebank_grammar(number_sentences))
    @assert number_sentences in (1, 5, 50, 500, 972)
    sentences = map(split, readlines(joinpath(dirname(@__FILE__), "WSJ_data", "WSJ.$number_sentences.sentences.txt")))
    scores = zeros(GeneralizedChartparsing.score_type(grammar), length(sentences))
    for (i,s) in enumerate(sentences)
        scores[i] = score(run_chartparser(s, grammar))
    end
    scores
end

function treebank_goldstandard_scores(number_sentences::Int)
    lines = readlines(joinpath(dirname(@__FILE__), "WSJ_data", "WSJ.$number_sentences.scores.txt"))
    map(parse, lines)
end


text:
function FieldStopOnIncreasing()
    return FieldStopOnIncreasing(())
end

function event_occurred(obj::FieldStopOnIncreasing, arg0::FieldSpacecraftState, arg1::FieldEventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (FieldSpacecraftState, FieldEventDetector, jboolean), arg0, arg1, arg2)
end


text:
using Test

using GradientBoost.Util
using GradientBoost.ML
using DelimitedFiles
using Statistics

using GLM: fit, predict, LinearModel

# Experiment on GBLearner
#
# gbl_func is a function that returns an instantiated GBLearner.
# score_func is a function that takes predictions, actual and returns as score.
# baseline_func is a function that takes labels.
function experiment(gbl_func, score_func, baseline_func,
  num_experiments, instances, labels)

  scores = zeros(Float64, num_experiments)
  for i = 1:num_experiments
    # Obtain training and test set
    (train_ind, test_ind) = holdout(size(instances, 1), 0.2)
    train_instances = instances[train_ind, :]
    test_instances = instances[test_ind, :]
    train_labels = labels[train_ind]
    test_labels = labels[test_ind]

    # Train
    gbl = gbl_func()
    fit!(gbl, train_instances, train_labels)

    # Test
    predictions = predict!(gbl, test_instances)
    score = score_func(predictions, test_labels)
    scores[i] = score
  end

  # Sanity check, score should be less than baseline.
  baseline = baseline_func(labels)
  @test mean(scores) <= baseline

  scores
end

# Error rate
function err_rate(predictions, actual)
  1.0 - mean(predictions .== actual)
end
function baseline_err_rate(labels)
  prop_ones = sum(labels) / length(labels)
  baseline = 1.0 - max(prop_ones, 1 - prop_ones)
end

# Mean squared error
function mse(predictions, actual)
  mean((actual .- predictions) .^ 2.0)
end
function baseline_mse(labels)
  label_mean = mean(labels)
  baseline = mse(label_mean, labels)
end

# Mean absolute deviation
function mad(predictions, actual)
  mean(abs.(actual .- predictions))
end
function baseline_mad(labels)
  label_median = median(labels)
  baseline = mad(label_median, labels)
end

num_experiments = 10

function ML.learner_fit(lf::LossFunction,
  learner::Type{LinearModel}, instances, labels)

  model = fit(learner, instances, labels)
end
function ML.learner_predict(lf::LossFunction,
  learner::Type{LinearModel}, model, instances)

  predict(model, instances)
end

@testset "System tests" begin
  @testset "iris dataset is handled by GBDT" begin
    # Get data
    iris = readdlm(joinpath(dirname(@__FILE__), "iris.csv"),',')
    instances = iris[:, 1:(end-1)]
    labels = iris[:, end]

    # Convert data to required format.
    instances = convert(Matrix{Float64}, instances)
    labels = [species == "setosa" ? 1.0 : 0.0 for species in labels]

    # Train and test multiple times (GBDT)
    function gbl_func()
      gbdt = GBDT(;
        loss_function=BinomialDeviance(),
        sampling_rate=0.6,
        learning_rate=0.1,
        num_iterations=100
      )
      gbl = GBLearner(gbdt, :class)
    end
    experiment(
      gbl_func, err_rate, baseline_err_rate, num_experiments, instances, labels
    )
  end

  @testset "mtcars dataset is handled" begin
    # Get data
    mtcars = readdlm(joinpath(dirname(@__FILE__), "mtcars.csv"),',')
    instances = mtcars[:, 2:end]
    labels = mtcars[:, 1]

    # Convert data to required format.
    instances = convert(Matrix{Float64}, instances)
    labels = convert(Vector{Float64}, labels)

    # Train and test multiple times (MSE)
    gbl_mse_funcs = Function[]
    function mse_gbdt_func()
      gbdt = GBDT(;
        loss_function=LeastSquares(),
        sampling_rate=0.6,
        learning_rate=0.1,
        num_iterations=100
      )
      gbl = GBLearner(gbdt, :regression)
    end
    push!(gbl_mse_funcs, mse_gbdt_func)
    function mse_gbl_func()
      gbl = GBBL(
        LinearModel;
        loss_function=LeastSquares(),
        sampling_rate=0.8,
        learning_rate=0.1,
        num_iterations=100,
      )
      gbl = GBLearner(gbl, :regression)
    end
    push!(gbl_mse_funcs, mse_gbl_func)
    for i = 1:length(gbl_mse_funcs)
      experiment(
        gbl_mse_funcs[i], mse, baseline_mse, num_experiments, instances, labels
      )
    end

    # Train and test multiple times (MAD)
    gbl_mad_funcs = Function[]
    function mad_gbdt_func()
      gbdt = GBDT(;
        loss_function=LeastAbsoluteDeviation(),
        sampling_rate=0.6,
        learning_rate=0.1,
        num_iterations=100
      )
      gbl = GBLearner(gbdt, :regression)
    end
    push!(gbl_mad_funcs, mad_gbdt_func)
    function mad_gbl_func()
      gbl = GBBL(
        LinearModel;
        loss_function=LeastAbsoluteDeviation(),
        sampling_rate=0.8,
        learning_rate=0.1,
        num_iterations=100,
      )
      gbl = GBLearner(gbl, :regression)
    end
    push!(gbl_mad_funcs, mad_gbl_func)
    for i = 1:length(gbl_mad_funcs)
      experiment(
        gbl_mad_funcs[i], mad, baseline_mad, num_experiments, instances, labels
      )
    end
  end
end


text:
using Parameters
using GeometricIntegrators.Equations

ϑ₁(t, q, v) = g₁₁(t,q) * v[1] + A₁(t, q)
ϑ₂(t, q, v) = g₂₂(t,q) * v[2] + A₂(t, q)
ϑ₃(t, q, v) = g₃₃(t,q) * v[3] + A₃(t, q)

function ϑ(t, q, v, θ)
    θ[1] = ϑ₁(t,q,v)
    θ[2] = ϑ₂(t,q,v)
    θ[3] = ϑ₃(t,q,v)
    nothing
end

v¹(t, q, p) = g¹¹(t, q) * (p[1] - A₁(t, q))
v²(t, q, p) = g²²(t, q) * (p[2] - A₂(t, q))
v³(t, q, p) = g³³(t, q) * (p[3] - A₃(t, q))

# ϕ₀(x) = E₀*sin(2π*x)
# ϕ(t,q) = ϕ₀(q[3])

# E₁(t,q) = zero(eltype(q))
# E₂(t,q) = zero(eltype(q))
# E₃(t,q) = - 2π*E₀*cos(2π*q[3])

function hamiltonian(t, q, p, params)
    @unpack μ = params
    0.5 * (g₁₁(t,q) * v¹(t,q,p)^2 + g₂₂(t,q) * v²(t,q,p)^2 + g₃₃(t,q) * v³(t,q,p)^2) + μ*B(t,q) + φ(t,q)
end

function initial_conditions(x₀, v₀)
    u₀ = v₀' * b(0, x₀)
    vpar = u₀ .* b⃗(0, x₀)
    vper = v₀ .- vpar
    μ = vper' * vper / 2 / B(0, x₀)

    (x₀, vpar, (μ = μ,))
end

function pauli_particle_3d_pᵢ(qᵢ, vᵢ, tᵢ=0)
    pᵢ = zero(qᵢ)

    if ndims(qᵢ) == 1
        ϑ(tᵢ, qᵢ, vᵢ, pᵢ)
    else
        for i in 1:size(qᵢ,2)
            @views ϑ(tᵢ, qᵢ[:,i], vᵢ[:,i], pᵢ[:,i])
        end
    end
    pᵢ
end

function pauli_particle_3d_pode_v(t, q, p, v, params)
    v[1] = v¹(t,q,p)
    v[2] = v²(t,q,p)
    v[3] = v³(t,q,p)
    nothing
end

function pauli_particle_3d_pode_f(t, q, p, f, params)
    @unpack μ = params
    f[1] = dA₁dx₁(t,q) * v¹(t,q,p) + dA₂dx₁(t,q) * v²(t,q,p) + dA₃dx₁(t,q) * v³(t,q,p) + E₁(t,q) - μ * dBdx₁(t,q) +
           (dg₁₁dx₁(t,q) * v¹(t,q,p)^2 + dg₂₂dx₁(t,q) * v²(t,q,p)^2 + dg₃₃dx₁(t,q) * v³(t,q,p)^2) / 2
    f[2] = dA₁dx₂(t,q) * v¹(t,q,p) + dA₂dx₂(t,q) * v²(t,q,p) + dA₃dx₂(t,q) * v³(t,q,p) + E₂(t,q) - μ * dBdx₂(t,q) +
           (dg₁₁dx₂(t,q) * v¹(t,q,p)^2 + dg₂₂dx₂(t,q) * v²(t,q,p)^2 + dg₃₃dx₂(t,q) * v³(t,q,p)^2) / 2
    f[3] = dA₁dx₃(t,q) * v¹(t,q,p) + dA₂dx₃(t,q) * v²(t,q,p) + dA₃dx₃(t,q) * v³(t,q,p) + E₃(t,q) - μ * dBdx₃(t,q) +
           (dg₁₁dx₃(t,q) * v¹(t,q,p)^2 + dg₂₂dx₃(t,q) * v²(t,q,p)^2 + dg₃₃dx₃(t,q) * v³(t,q,p)^2) / 2
    nothing
end

pauli_particle_3d_iode_ϑ(t, q, v, θ, params) = ϑ(t, q, v, θ)

function pauli_particle_3d_iode_f(t, q, v, f, params)
    @unpack μ = params
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] + E₁(t,q) - μ * dBdx₁(t,q) +
           (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] + E₂(t,q) - μ * dBdx₂(t,q) +
           (dg₁₁dx₂(t,q) * v[1]^2 + dg₂₂dx₂(t,q) * v[2]^2 + dg₃₃dx₂(t,q) * v[3]^2) / 2
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] + E₃(t,q) - μ * dBdx₃(t,q) +
           (dg₁₁dx₃(t,q) * v[1]^2 + dg₂₂dx₃(t,q) * v[2]^2 + dg₃₃dx₃(t,q) * v[3]^2) / 2
    nothing
end

function pauli_particle_3d_iode_g(t, q, v, f, params)
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] +
           (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] +
           (dg₁₁dx₂(t,q) * v[1]^2 + dg₂₂dx₂(t,q) * v[2]^2 + dg₃₃dx₂(t,q) * v[3]^2) / 2
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] +
           (dg₁₁dx₃(t,q) * v[1]^2 + dg₂₂dx₃(t,q) * v[2]^2 + dg₃₃dx₃(t,q) * v[3]^2) / 2
    nothing
end

function pauli_particle_3d_pode(q₀, v₀, params)
    PODE(pauli_particle_3d_pode_v, pauli_particle_3d_pode_f, q₀, pauli_particle_3d_pᵢ(q₀, v₀);
         parameters=params, h=hamiltonian)
end

pauli_particle_3d_pode(qᵢ=qᵢ, vᵢ=vᵢ) = pauli_particle_3d_pode(initial_conditions(qᵢ, vᵢ)...)

function pauli_particle_3d_iode(q₀::AbstractVector, v₀::AbstractVector, params::NamedTuple)
    IODE(pauli_particle_3d_iode_ϑ, pauli_particle_3d_iode_f, pauli_particle_3d_iode_g,
            q₀, pauli_particle_3d_pᵢ(q₀, v₀);
            parameters=params, h=hamiltonian, v = (t, q, v, params) -> nothing)
end

pauli_particle_3d_iode(q₀::AbstractVector, v₀::AbstractVector, μ::Real) = pauli_particle_3d_iode(q₀, v₀, (μ=μ,))

pauli_particle_3d_iode(qᵢ=qᵢ, vᵢ=vᵢ) = pauli_particle_3d_iode(initial_conditions(qᵢ, vᵢ)...)


text:
using IPNets
using Test
import Sockets: IPAddr, IPv4, IPv6, @ip_str

@testset "IPNets" begin
    ip41 = IPv4("1.2.3.4")
    ip42 = IPv4("5.6.7.8")

    n1 = IPv4Net("1.2.3.0/24")
    n2 = IPv4Net("1.2.3.4/24")
    n3 = IPv4Net("1.2.3.4/26")
    n4 = IPv4Net("5.6.7.0/24")
    n5 = IPv4Net("1.2.3.4/30")

    ip61 = IPv6("2001:1::4")
    ip62 = IPv6("2001:2::8")

    o1 = IPv6Net("2001:1::/64")
    o2 = IPv6Net("2001:1::4/64")
    o3 = IPv6Net("2001:1::4/68")
    o4 = IPv6Net("2001:2::8/64")
    o5 = IPv6Net("2001:1::4/126")

    # ipv4
    @testset "IPv4" begin
        @test_throws ErrorException IPv4Net(ip41, 33)
        @test IPv4Net("1.2.3.4", "255.255.255.0") == n1
        @test IPv4Net("1.2.3.4", 24) == n1
        @test IPv4Net(16909060,24) == n1
        @test IPv4Net((16909060,24)) == n1
        @test IPv4Net("1.2.3.4") == IPv4Net("1.2.3.4/32")
        @test n1 == n2
        @test isless(n1,n3) == false
        @test isless(n1,n4) == true
        @test n1[5] == ip41
        @test isless(ip41, ip42) == true
        @test in(ip42, n4) == true
        @test_throws ErrorException ip42 in o4
        @test contains(n4, ip42) == true
        @test issubset(n3, n2) == true
        @test issubset(n1, n2) == true
        @test IPNets._contiguousbitcount(240,UInt8) == 0x04
        @test IPNets._contiguousbitcount(252,UInt8) == 0x06
        @test_throws ErrorException IPNets._contiguousbitcount(241,UInt8)
        @test endswith(string(n5),"(\"1.2.3.4/30\")") == true
        @test endswith(sprint(print,n5),"(\"1.2.3.4/30\")") == true
        @test endswith(sprint(show,n5),"(\"1.2.3.4/30\")") == true
        @test endswith(string(display,n5),"(\"1.2.3.4/30\")") == true
        @test size(n5) == (4,)
        @test [x for x in n5] == [ip"1.2.3.4", ip"1.2.3.5", ip"1.2.3.6", ip"1.2.3.7"]
        @test lastindex(n5) == 4
        @test minimum(n5) == ip"1.2.3.4"
        @test maximum(n5) == ip"1.2.3.7"
        @test extrema(n5) == (ip"1.2.3.4",ip"1.2.3.7")
        @test getindex(n5,1:2) == [ip"1.2.3.4", ip"1.2.3.5"]
        # @test getindex(n5,(1,)) == ip41
        @test_throws BoundsError getindex(n5, 10)
        @test IPNets.width(IPv4) == 32
        @test_throws BoundsError IPNets._mask2bits(IPv4, UInt64(33))
        @test netmask(n1) == ip"255.255.255.0"
        @test eltype(n1) == IPv4
        @test n5[1] == ip41
        @test getindex(n5,1) == ip41
        @test IPNet("1.2.3.0/24") == n1
    end

    # IPv6
    @testset "IPv6" begin
        @test_throws ErrorException IPv6Net(ip61, 129)
        @test IPv6Net("2001:1::", 64) == o1
        @test IPv6Net(0x20010001000000000000000000000004,64) == o1
        @test IPv6Net((0x20010001000000000000000000000004,64)) == o1
        @test o1 == o2

        @test IPv6Net("2001:1::1") == IPv6Net("2001:1::1/128")
        @test isless(o1,o3) == false
        @test isless(o1,o4) == true
        @test o1[5] == ip61
        @test isless(ip61, ip62) == true
        @test in(ip62, o4) == true
        @test_throws ErrorException ip62 in n4
        @test contains(o4, ip62) == true
        @test endswith(string(o5),"(\"2001:1::4/126\")") == true
        @test endswith(sprint(print,o5),"(\"2001:1::4/126\")") == true
        @test endswith(sprint(show,o5),"(\"2001:1::4/126\")") == true
        @test endswith(string(display,o5),"(\"2001:1::4/126\")") == true
        @test size(o5) == (4,)
        @test [x for x in o5] == [ip"2001:1::4",ip"2001:1::5",ip"2001:1::6",ip"2001:1::7"]
        @test lastindex(o5) == 4
        @test minimum(o5) == ip"2001:1::4"
        @test maximum(o5) == ip"2001:1::7"
        @test extrema(o5) == (ip"2001:1::4",ip"2001:1::7")
        @test getindex(o5,1:2) == [ip"2001:1::4", ip"2001:1::5"]
        # @test getindex(o5,(1,)) == ip61
        @test_throws BoundsError getindex(o5, 10)
        @test IPNets.width(IPv6) == 128
        @test_throws BoundsError IPNets._mask2bits(IPv6, UInt64(129))
        @test eltype(o1) == IPv6
        @test o5[1] == ip61
        @test getindex(o5,1) == ip61
        @test IPNet("2001:1::/64") == o1
    end
end


text:
### BasicMCRange

# BasicMCRange is used for sampling a single parameter via serial Monte Carlo
# It is the most elementary and typical Markov chain Monte Carlo (MCMC) range
# It contains nagivational info (burnin, thinning, number of steps)

immutable BasicMCRange{T<:Integer} <: MCRange{T}
  burnin::T
  thinning::T
  nsteps::T
  postrange::StepRange{T, T}
  npoststeps::T

  function BasicMCRange(postrange::StepRange{T, T})
    burnin = first(postrange)-1
    thinning = postrange.step
    nsteps = last(postrange)

    @assert burnin >= 0 "Number of burn-in iterations should be non-negative"
    @assert thinning >= 1 "Thinning should be >= 1"
    @assert nsteps > burnin "Total number of MCMC iterations should be greater than number of burn-in iterations"

    npoststeps = length(postrange)

    new(burnin, thinning, nsteps, postrange, npoststeps)
  end
end

BasicMCRange{T<:Integer}(postrange::StepRange{T, T}) = BasicMCRange{T}(postrange)

BasicMCRange{T<:Integer}(postrange::UnitRange{T}) = BasicMCRange{T}(first(postrange):1:last(postrange))

BasicMCRange{T<:Integer}(; burnin::T=0, thinning::T=1, nsteps::T=100) = BasicMCRange{T}((burnin+1):thinning:nsteps)

Base.show(io::IO, r::BasicMCRange) =
  print(io, "BasicMCRange: number of steps = $(r.nsteps), burnin = $(r.burnin), thinning = $(r.thinning)")

Base.writemime(io::IO, ::MIME"text/plain", r::BasicMCRange) = show(io, r)


text:
# Model hyperparameters

"""
RAPIDS Docs for Linear Regression:
    https://docs.rapids.ai/api/cuml/stable/api.html#linear-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = LinearRegression()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct LinearRegression <: MMI.Deterministic
    handle = nothing
    algorithm::String = "eig"::(_ in ("svd", "eig", "qr", "svd-qr", "svd-jacobi"))
    fit_intercept::Bool = true
    normalize::Bool = false
    verbose::Bool = false
end

"""
RAPIDS Docs for Ridge Regression:
    https://docs.rapids.ai/api/cuml/stable/api.html#ridge-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = Ridge()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct Ridge <: MMI.Deterministic
    handle = nothing
    alpha::Float64 = 1.0::(_ > 0)
    solver::String = "eig"::(_ in ("svd", "eig", "cd"))
    fit_intercept::Bool = true
    normalize::Bool = false
    verbose::Bool = false
end

"""
RAPIDS Docs for Lasso Regression:
    https://docs.rapids.ai/api/cuml/stable/api.html#lasso-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = Lasso()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct Lasso <: MMI.Deterministic
    handle = nothing
    alpha::Float64 = 1.0::(_ > 0)
    fit_intercept::Bool = true
    normalize::Bool = false
    max_iter::Int = 1000::(_ > 0)
    tol::Float64 = 1e-3::(_ > 0)
    selection::String = "cyclic"::(_ in ("cyclic", "random"))
    verbose::Bool = false
end

"""
RAPIDS Docs for ElasticNet:
    https://docs.rapids.ai/api/cuml/stable/api.html#elasticnet-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = ElasticNet()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct ElasticNet <: MMI.Deterministic
    handle = nothing
    alpha::Float64 = 1.0::(_ > 0)
    l1_ratio::Float64 = 0.5::(_ > 0)
    fit_intercept::Bool = true
    normalize::Bool = false
    max_iter::Int = 1000::(_ > 0)
    tol::Float64 = 1e-3::(_ > 0)
    selection::String = "cyclic"::(_ in ("cyclic", "random"))
    verbose::Bool = false
end

"""
RAPIDS Docs for MBSGDRegressor:
    https://docs.rapids.ai/api/cuml/stable/api.html#mini-batch-sgd-regressor
Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = MBSGDRegressor()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct MBSGDRegressor <: MMI.Deterministic
    handle = nothing
    loss::String = "squared_loss"::(_ in ("squared_loss", ))
    penalty::String = "none"::(_ in ("none", "l1", "l2", "elasticnet"))
    alpha::Float64 = 0.0001::(_ >= 0)
    fit_intercept::Bool = true
    l1_ratio::Float64 = 0.15::(_ >= 0 && _ <= 1)
    batch_size::Int = 32::(_ > 0)
    epochs::Int = 1000::(_ > 0)
    tol::Float64 = 1e-3::(_ > 0)
    shuffle::Bool = true
    eta0::Float64 = 1e-3::(_ > 0)
    power_t::Float64 = 0.5::(_ > 0)
    learning_rate::String = "constant"::(_ in ("adaptive", "constant", "invscaling", "optimal"))
    n_iter_no_change::Int = 5::(_ > 0)
    verbose::Bool = false
end

"""
RAPIDS Docs for RandomForestRegressor:
    https://docs.rapids.ai/api/cuml/stable/api.html#random-forest
Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = RandomForestRegressor()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct RandomForestRegressor <: MMI.Deterministic
    handle = nothing
    n_estimators::Int = 100::(_ > 0)
    split_criterion = 2
    bootstrap::Bool = true
    max_samples::Float64 = 1.0::(_ > 0)
    max_depth::Int = 16::(_ > 0)
    max_leaves::Int = -1
    max_features = "auto"
    n_bins::Int = 128::(_ > 0)
    n_streams::Int = 4::(_ > 0)
    min_samples_leaf::Int = 1::(_ > 0)
    min_samples_split::Int = 2::(_ > 1)
    min_impurity_decrease::Float64 = 0.0::(_ >= 0)
    accuracy_metric::String = "r2"::(_ in ("median_ae", "mean_ae", "mse", "r2"))
    max_batch_size::Int = 4096::(_ > 0)
    random_state = nothing
    verbose::Bool = false
end

"""
RAPIDS Docs for Coordinate Descent:
    https://docs.rapids.ai/api/cuml/stable/api.html#coordinate-descent

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = CD()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct CD <: MMI.Deterministic
    handle = nothing
    loss::String = "squared_loss"::(_ in ("squared_loss", ))
    alpha::Float64 = 1.0::(_ > 0)
    l1_ratio::Float64 = 0.5::(_ > 0)
    fit_intercept::Bool = true
    max_iter::Int = 1000::(_ > 0)
    tol::Float64 = 1e-3::(_ > 0)
    shuffle::Bool = true
    verbose::Bool = false
end

"""
RAPIDS Docs for SVR:
    https://docs.rapids.ai/api/cuml/stable/api.html#support-vector-machines

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = SVR()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct SVR <: MMI.Deterministic
    handle = nothing
    C::Float64 = 1.0::(_ >= 0)
    kernel::String = "rbf"::(_ in ("linear", "poly", "rbf", "sigmoid"))
    degree::Int = 3::(_ > 0)
    gamma = "scale"
    coef0::Float64 = 0.0
    tol::Float64 = 1e-3::(_ > 0)
    epsilon::Float64 = 0.1::(_ >= 0)
    cache_size::Float64 = 1024.0::(_ > 0)
    max_iter::Int = -1
    nochange_steps::Int = 1000::(_ >= 0)
    verbose::Bool = false
end

"""
RAPIDS Docs for LinearSVR:
    https://docs.rapids.ai/api/cuml/stable/api.html#support-vector-machines

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = SVLinearSVR()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct LinearSVR <: MMI.Deterministic
    handle = nothing
    penalty::String = "l2"::(_ in ("l1", "l2"))
    loss = "epsilon_insensitive"::(_ in ("epsilon_insensitive", "squared_epsilon_insensitive"))
    fit_intercept::Bool = true
    penalized_intercept::Bool = true
    max_iter::Int = 1000::(_ > 0)
    linesearch_max_iter::Int = 100::(_ > 0)
    lbfgs_memory::Int = 5::(_ > 0)
    C::Float64 = 1.0::(_ >= 0)
    grad_tol::Float64 = 1e-4::(_ > 0)
    change_tol::Float64 = 1e-5::(_ > 0)
    tol = nothing
    epsilon::Float64 = 0.0::(_ >= 0)
    verbose::Bool = false
end

"""
RAPIDS Docs for Nearest Neighbors Regression:
    https://docs.rapids.ai/api/cuml/stable/api.html#nearest-neighbors-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = rand(100)

model = KNeighborsRegressor()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct KNeighborsRegressor <: MMI.Deterministic
    handle = nothing
    algorithm::String = "brute"::(_ in ("brute",))
    metric::String = "euclidean"
    weights::String = "uniform"::(_ in ("uniform",))
    verbose::Bool = false
end

# Multiple dispatch for initializing models
model_init(mlj_model::LinearRegression) = cuml.LinearRegression(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::Ridge) = cuml.Ridge(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::Lasso) = cuml.Lasso(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::ElasticNet) = cuml.ElasticNet(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::MBSGDRegressor) = cuml.MBSGDRegressor(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::RandomForestRegressor) = cuml.RandomForestRegressor(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::CD) = cuml.CD(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::SVR) = cuml.svm.SVR(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::LinearSVR) = cuml.svm.LinearSVR(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::KNeighborsRegressor) = cuml.KNeighborsRegressor(; mlj_to_kwargs(mlj_model)...)

const CUML_REGRESSION = Union{LinearRegression,
                                Ridge,
                                Lasso,
                                ElasticNet,
                                MBSGDRegressor,
                                RandomForestRegressor,
                                CD,
                                SVR,
                                LinearSVR,
                                KNeighborsRegressor}

# add metadata
MMI.load_path(::Type{<:LinearRegression}) = "$PKG.LinearRegression"
MMI.load_path(::Type{<:Ridge}) = "$PKG.Ridge"
MMI.load_path(::Type{<:Lasso}) = "$PKG.Lasso"
MMI.load_path(::Type{<:ElasticNet}) = "$PKG.ElasticNet"
MMI.load_path(::Type{<:MBSGDRegressor}) = "$PKG.MBSGDRegressor"
MMI.load_path(::Type{<:RandomForestRegressor}) = "$PKG.RandomForestRegressor"
MMI.load_path(::Type{<:CD}) = "$PKG.CD"
MMI.load_path(::Type{<:SVR}) = "$PKG.SVR"
MMI.load_path(::Type{<:LinearSVR}) = "$PKG.LinearSVR"
MMI.load_path(::Type{<:KNeighborsRegressor}) = "$PKG.KNeighborsRegressor"

MMI.input_scitype(::Type{<:CUML_REGRESSION}) = Union{AbstractMatrix, Table(Continuous)}
MMI.target_scitype(::Type{<:CUML_REGRESSION}) = Union{AbstractVector, Table(Continuous)} #should we use count even though we input floats?

MMI.docstring(::Type{<:LinearRegression}) = "cuML's LinearRegression: https://docs.rapids.ai/api/cuml/stable/api.html#linear-regression"
MMI.docstring(::Type{<:Ridge}) = "cuML's Ridge: https://docs.rapids.ai/api/cuml/stable/api.html#ridge-regression"
MMI.docstring(::Type{<:Lasso}) = "cuML's Lasso: https://docs.rapids.ai/api/cuml/stable/api.html#lasso-regression"
MMI.docstring(::Type{<:ElasticNet}) = "cuML's ElasticNet: https://docs.rapids.ai/api/cuml/stable/api.html#elasticnet-regression"
MMI.docstring(::Type{<:MBSGDRegressor}) = "cuML's MBSGDRegressor: https://docs.rapids.ai/api/cuml/stable/api.html#mini-batch-sgd-regressor"
MMI.docstring(::Type{<:AbstractMatrix}) = "cuML's AbstractMatrix: https://docs.rapids.ai/api/cuml/stable/api.html#random-forest"
MMI.docstring(::Type{<:CD}) = "cuML's CD: https://docs.rapids.ai/api/cuml/stable/api.htmll#coordinate-descent"
MMI.docstring(::Type{<:SVR}) = "cuML's SVR: https://docs.rapids.ai/api/cuml/stable/api.html#support-vector-machines"
MMI.docstring(::Type{<:LinearSVR}) = "cuML's LinearSVR: https://docs.rapids.ai/api/cuml/stable/api.html#support-vector-machines"
MMI.docstring(::Type{<:KNeighborsRegressor}) = "cuML's KNeighborsRegressor: https://docs.rapids.ai/api/cuml/stable/api.html#nearest-neighbors-regression"

# fit methods
function MMI.fit(mlj_model::CUML_REGRESSION, verbosity, X, y, w=nothing)
    # fit the model
    model = model_init(mlj_model)
    model.fit(prepare_input(X), prepare_input(y))
    fitresult = model

    # save result
    cache = nothing
    report = ()
    return (fitresult, cache, report)
end

# predict methods
function MMI.predict(mlj_model::CUML_REGRESSION, fitresult, Xnew)
    model  = fitresult
    py_preds = model.predict(prepare_input(Xnew))
    preds = pyconvert(Array, py_preds)

    return preds
end

# Regression metadata
MMI.metadata_pkg.((LinearRegression,
                    Ridge,
                    Lasso,
                    ElasticNet,
                    MBSGDRegressor,
                    RandomForestRegressor,
                    CD,
                    SVR,
                    LinearSVR,
                    KNeighborsRegressor),
    name = "cuML Regression Methods",
    uuid = "2764e59e-7dd7-4b2d-a28d-ce06411bac13", # see your Project.toml
    url  = "https://github.com/tylerjthomas9/RAPIDS.jl",  # URL to your package repo
    julia = false,          # is it written entirely in Julia?
    license = "MIT",        # your package license
    is_wrapper = true,      # does it wrap around some other package?
)


text:
function load_longbase(longbase_filename, m::Model=FRBUS_VAR.model)
    models_path = joinpath(realpath(dirname(@__FILE__)), "models")
    longbase_lines = readlines(joinpath(models_path, longbase_filename))
    function fixdate(x)
        for subs in Dict("-01-01"=>"Q1", "-04-01"=>"Q2", "-07-01"=>"Q3", "-10-01"=>"Q4")
            x = replace(x, subs)
        end
        return x
    end
    start = eval(Meta.parse(fixdate(longbase_lines[2][1:10])))
    data = SimData(start, m.varshks, zeros(Float64, length(longbase_lines)-1, length(m.varshks)))
    names = tuple(Symbol.(split(longbase_lines[1], ",")[2:end])...)
    valid_names = Set(m.variables)
    for line = longbase_lines[2:end]
        date = eval(Meta.parse(fixdate(line[1:10])))
        values = split(line, ",")[2:end]
        @assert length(names) == length(values)
        for (n,v) in zip(names, values)
            if v == "NA" || n ∉ valid_names
                continue
            end
            data[date, n] = parse(Float64, v)
        end
    end
    data
end

nothing


text:
using JuMP
using Clp

LARGE_VALUE = 100000
DEMANDS = [400 300 420 380]
SUPPLY = [500, 600, 200, 300]
COSTS = [100 103 106 109;
        LARGE_VALUE 140 143 146;
        LARGE_VALUE LARGE_VALUE 120 123;
        LARGE_VALUE LARGE_VALUE LARGE_VALUE 150]

m = Model(Clp.Optimizer)
@variable(m, x[1:4, 1:4] ≥ 0)
@objective(m, Min, sum(COSTS .* x))
@constraint(m, sum(x, dims=2) .≤ SUPPLY)
@constraint(m, sum(x, dims=1) .== DEMANDS)
print(m)
optimize!(m)
status = termination_status(m)
println("Solution status: ", status)
println("Objective value: ", objective_value(m))
println("Values: ", value.(x))


text:
#!/usr/bin/env julia

module NRI

import ..ReconstructionEvaluations
const RE = ReconstructionEvaluations;

"""

   nri( om::SparseMatrixCSC; correct=true )

Implementation of the Neural Reconstruction Index. Incorporates
the correction for false positive values if "correct" is true.
"""
function nri( om::SparseMatrixCSC; correct=true )

  TP, segTP = compute_TPs( om )
  FP, segFP = compute_FPs( om; correct=correct )
  FN, segFN = compute_FNs( om )

  NRI = nri( TP, FP, FN )
  segNRI = nri( segTP, segFP, segFN )

  segNRIw = nri_weight( segTP, segFP, segFN )

  NRI, segNRI, segNRIw
end

function nri( om::Array; correct=true )
  nri( sparse(om); correct=correct )
end

function nri( TP, FP, FN )
  2*TP ./ (2*TP + FP + FN)
end

function nri( TP::SparseVector, FP::SparseVector, FN::SparseVector )

  inds = find(TP); #if a segment has no TP's, segNRI=0
  res = sparsevec(Int[],Float64[],length(TP))

  for i in inds res[i] = (2*TP[i]) / (2*TP[i] + FP[i] + FN[i]) end

  res
end

function nri( table1::AbstractArray, table2::AbstractArray; correct=true )
  count_table, A_to_inds, B_to_inds = RE.build_count_table(table1,table2)

  NRI, segNRI, segNRIw = nri( count_table; correct=correct )

  res_segNRI  = sparsevec( Int[], Float64[], maximum(keys(A_to_inds))+1 )
  res_segNRIw = sparsevec( Int[], Float64[], maximum(keys(A_to_inds))+1 )
  for (k,v) in A_to_inds
      res_segNRI[k] = segNRI[v]
      res_segNRIw[k] = segNRIw[v]
  end

  NRI, res_segNRI, res_segNRIw
end

function nri( fname1::AbstractString, fname2::AbstractString )

  @assert isfile(fname1) && isfile(fname2)

  table1, table2 = RE.load_edges(fname1), RE.load_edges(fname2)

  nri( table1, table2 )
end

"""

    compute_TPs( om, offset=true )

Counts the number of true positive paths within each
segment in a network specified by a sparse overlap matrix,
and summarizes their total.

Offset indicates whether segment ids are offset by 1 to include
synapse "insertions" and deletions. This feature is not fully implemented
yet.
"""
function compute_TPs( om, offset=true )

  rs, cs = findn(om); vs = nonzeros(om);

  TPs = choosetwo(vs);
  segTPs = sparsevec(Int[],Int[],size(om,1));

  oT = one(eltype(rs));
  for (r,c,tp) in zip(rs,cs,TPs)

    if r == oT continue end
    if c == oT continue end

    segTPs[r] += tp
  end

  sum(segTPs), segTPs
end

"""

    compute_FPs( om, offset=true; correct=true )

Counts the number of false positive paths within each
segment in a network specified by a sparse overlap matrix,
and summarizes their total.

offset indicates whether segment ids are offset by 1 to include
synapse "insertions" and deletions. This feature is not fully implemented
yet.

correct indicates whether to fully penalize false positives (unlike within
the first draft of MATLAB nri code). The false setting is meant to be
used as a comparison
"""
function compute_FPs( om, offset=true; correct=true )

  rs, cs = findn(om); vs = nonzeros(om);

  sum_over_cols = sparsevec(Int[],Int[],size(om,2))
  for (c,v) in zip(cs, vs) sum_over_cols[c] += v end

  #needs to be a float to account for attributing
  # mergers to gt segments
  segFPs = sparsevec(Int[],Float64[],size(om,1))

  oT = one(eltype(rs))
  for (r,c,v) in zip(rs,cs,vs)

    if c == oT continue end

    other_count = sum_over_cols[c] - v
    segFPs[r] += (v*other_count) / 2

    #fully accounting for synapse insertions (if using the correct vers)
    if r == oT && correct segFPs[r] += choosetwo(v) end
  end

  sum(segFPs), segFPs
end

"""

    compute_FNs( om, offset=true )

Counts the number of false negative paths within each
segment in a network specified by a sparse overlap matrix,
and summarizes their total.

offset indicates whether segment ids are offset by 1 to include
synapse "insertions" and deletions. This feature is not fully implemented
yet.

correct indicates whether to fully penalize false positives (unlike within
the first draft of MATLAB nri code). The false setting is meant to be
used as a comparison
"""
function compute_FNs( om, offset=true )

  rs, cs = findn(om); vs = nonzeros(om);

  sum_over_rows = sparsevec(Int[],Int[],size(om,1))
  for (r,v) in zip(rs, vs) sum_over_rows[r] += v end

  segFNs = sparsevec(Int[],Float64[],size(om,1))

  oT = one(eltype(rs))
  for (r,c,v) in zip(rs,cs,vs)

    if r == oT continue end

    other_count = sum_over_rows[r] - v
    segFNs[r] += (v*other_count) / 2

    #accounting for synapse deletions
    if c == oT segFNs[r] += choosetwo(v) end
  end

  sum(segFNs), segFNs
end

"""

    nri_weight( segTP, segFP, segFN )

Finds the weight for each segment s.t. the full network
NRI is a weighted sum of each individual NRI values. This value
is equal to the sum of "positives" (TP+FP) and "negatives" (TP+FN)
for each segment.
"""
function nri_weight( segTP, segFP, segFN )

  total_weight = 2*sum(segTP) + sum(segFP) + sum(segFN)

  is = Set(find(segTP)); union!(is, find(segFP)); union!(is, find(segFN));

  res = sparsevec(Int[],Float64[],length(segTP))

  for i in is res[i] = (2*segTP[i] + segFP[i] + segFN[i]) / total_weight end

  res
end

"""

    nri_by_class( seg_nri, seg_nriw, class_map )
"""
function nri_by_class( seg_nri, seg_nriw, class_map )

  is = find(seg_nriw)

  class_nri = Dict( k => 0. for k in unique(values(class_map)) )
  class_w   = Dict( k => 0. for k in unique(values(class_map)) )

  for i in is
    if !haskey(class_map,i) continue end

    class_w[class_map[i]] += seg_nriw[i]
  end

  for i in is
    if !haskey(class_map,i) continue end
    class = class_map[i]

    class_nri[class] += (seg_nri[i] * seg_nriw[i]) / class_w[class]
  end

  class_nri, class_w
end

#Utility fns
choosetwo(x) = x * ((x-1)/2)
choosetwo!(x) = for i in eachindex(x) x[i] *= (x[i]-1)/2 end
choosetwo(x::Array) = [ v*((v-1)/2) for v in x ]

end #module end


text:
type Conv2d <: Layer{Float32}
    name::Symbol
    W::Array{Float32, 4} # filter_h × filter_w × in_c × out_c
    b::Vector{Float32} # out_c
    strides::Tuple{Int, Int} # width, height

    parent::Array{Float32, 3} # h × w × in_c
    child::Array{Float32, 3}  # out_h × out_w × out_c
end
get_name(a::Conv2d) = a.name
output(a::Conv2d) = a.child
function add_node!(net::ForwardNet, ::Type{Conv2d},
    name::Symbol,
    parent_index::Int,
    strides::Tuple{Int, Int},
    filter_h::Int,
    filter_w::Int,
    out_c::Int)

    parent = output(net[parent_index])::Array{Float32, 3}
    h, w, in_c = size(parent)

    W = Array{Float32}(filter_h, filter_w, in_c, out_c)
    b = Array{Float32}(out_c)

    # FOR VALID: (padding is always zero)
    # out_w = ceil(float(w - filter_w + 1) / float(strides[2]))
    # out_h = ceil(float(h - filter_h + 1) / float(strides[1]))

    # FOR SAME:
    out_w  =ceil(Int, w / strides[2])
    out_h = ceil(Int, h / strides[1])
    child = Array{Float32}(out_h, out_w, out_c)

    node = Conv2d(name, W, b, strides, parent, child)
    add_node!(net, node, parent_index)
end
function forward!(a::Conv2d)

    h_in = size(a.parent, 1)
    h_out = size(a.child, 1)
    w_in = size(a.parent, 2)
    w_out = size(a.child, 2)
    hh, ww, f_in, f_out = size(a.W)

    s₁, s₂ = a.strides
    pad_h = (s₁ * (h_out - 1) + hh - h_in) / 2
    pad_w = (s₂ * (w_out - 1) + ww - w_in) / 2

    pad_h_lo = floor(Int, pad_h)
    pad_h_hi = ceil(Int, pad_h)
    pad_w_lo = floor(Int, pad_w)
    pad_w_hi = ceil(Int, pad_w)

    # println("parent_size: ", size(a.parent))
    # println("child_size:  ", size(a.child))

    # println("padding: ", pad_h_lo, "  ", pad_h_hi)
    # println("         ", pad_w_lo, "  ", pad_w_hi)

    i_out = 1
    for i_out in 1 : size(a.child, 1)
        i_in_lo = (i_out-1)*s₁ - pad_h_lo + 1
        i_in_hi = i_in_lo + hh - 1

        h_in_lo = 1
        if i_in_lo < 1
            h_in_lo += 1-i_in_lo
            i_in_lo = 1
        end
        h_in_hi = hh
        if i_in_hi > h_in
            h_in_hi -= i_in_hi - h_in
            i_in_hi = h_in
        end

        for j_out in 1 : size(a.child, 2)

            j_in_lo = (j_out-1)*s₂ - pad_w_lo + 1
            j_in_hi = j_in_lo + ww - 1

            # println("j_in_lo: ", j_in_lo)
            # println("j_in_hi: ", j_in_hi)

            w_in_lo = 1
            if j_in_lo < 1
                w_in_lo += 1-j_in_lo
                j_in_lo = 1
            end
            w_in_hi = ww
            if j_in_hi > w_in
                w_in_hi -= j_in_hi - w_in
                j_in_hi = w_in
            end

            # println("j_in_lo: ", j_in_lo)
            # println("j_in_hi: ", j_in_hi)
            # println("w_in_lo: ", w_in_lo)
            # println("w_in_hi: ", w_in_hi)

            for k_out in 1 : f_out

                # println("parent: ", i_in_lo:i_in_hi, ", ", j_in_lo:j_in_hi, ", ", 1:f_in)
                # println("W:      ", h_in_lo:h_in_hi, ", ", w_in_lo,w_in_hi, ", ", 1:f_in, ", ", k_out)

                a.child[i_out,j_out,k_out] = a.b[k_out]

                for (i, i_in) in enumerate(i_in_lo:i_in_hi)
                    h₁ = h_in_lo + i - 1
                    for (j, j_in) in enumerate(j_in_lo:j_in_hi)
                        w₁ = w_in_lo + j - 1
                        for k_in in 1 : f_in
                            W_contrib = a.W[h₁,w₁,k_in,k_out]
                            p_contrib = a.parent[i_in,j_in,k_in]
                            a.child[i_out,j_out,k_out] += W_contrib*p_contrib
                        end
                    end
                end
            end
        end
    end

    a
end
function restore!(a::Conv2d, filename_W::String, filename_b::String)

    filter_h, filter_w, in_c, out_c = size(a.W)

    vec = open(read_binary_vec, filename_W)
    a.W[:] = convert_to_column_major_array(vec, (filter_h, filter_w, in_c, out_c))

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end


text:
module SubModule
    import XUnit

    mutable struct MyTestMetrics <: XUnit.TestMetrics
        default::XUnit.DefaultTestMetrics
    end

    function MyTestMetrics()
        return MyTestMetrics(XUnit.DefaultTestMetrics())
    end

    function XUnit.run_and_gather_test_metrics(fn::Function, ts::XUnit.AbstractTestSet, m::MyTestMetrics; run::Bool=true)
        return XUnit.run_and_gather_test_metrics(fn, ts, m.default; run=run)
    end

    function XUnit.combine_test_metrics(parent::MyTestMetrics, sub::MyTestMetrics)
        return XUnit.combine_test_metrics(parent.default, sub.default)
    end

    function XUnit.save_test_metrics(ts, m::MyTestMetrics)
        return XUnit.save_test_metrics(ts, m.default)
    end
end

module MacroDefMod

import ..SubModule
import XUnit

macro my_testcase(name, code)
    tests_block_location = XUnit.get_block_source(code)
    tests_block = if tests_block_location !== nothing
        Expr(:block, tests_block_location, esc(code))
    else
        quote
            begin
                esc(code)
            end
        end
    end
    quote
        XUnit.@testcase metrics = SubModule.MyTestMetrics $(name) $(tests_block)
    end
end

end