text:
function single(iter)
    i = iterate(iter)
    if i === nothing
        throw(ArgumentError("Collection must have exactly one element (input is empty)"))
    end
    (out, state) = i
    if iterate(iter, state) === nothing
        return out
    else
        throw(ArgumentError("Collection must have exactly one element (input contains more than one element)"))
    end
end

function single(::Tuple{})
    throw(ArgumentError("Collection must have exactly one element (input is empty)"))
end
function single(t::Tuple{Any})
    return t[1]
end
function single(::NTuple{N,Any}) where N
    throw(ArgumentError("Collection must have exactly one element (input has $N elements)"))
end


text:
# Pre- and post- processing of behavioral and simulated data.
# file IO
using DelimitedFiles

function log_print!(buffer::IOBuffer, line::String)
    print(line)
    print(buffer, line)
end

function log_print_ln!(buffer::IOBuffer, line::String)
    println(line)
    println(buffer, line)
end

@enum SARSPE_COLS c_state=1 c_action=2 c_reward=3 c_nextState=4 c_puzzle=5 c_episode=6 c_subject=7 c_surprise_trial=8 c_RT_prctile=9 c_action_value=10 c_graph_version=11

function filter_by_person(data::Array{Int,2}, person_id::Int)
    idx = data[:,Int(c_subject)] .== person_id
    data[idx,:]
end

function filter_by_person(data::Array{Int,2}, person_ids::Set)
    filtered_data = Array{Int, 2}(undef, 0, size(data)[2])
    for id in person_ids
        subj_data = filter_by_person(data, id)
        filtered_data = vcat(filtered_data, subj_data)
    end
    filtered_data
end

function filter_by_max_episode(data, max_episode_id)
    idx = data[:,Int(c_episode)].<=max_episode_id
    data[idx,:]
end

function filter_by_episode(data, episode_id)
    idx = data[:,Int(c_episode)].==episode_id
    data[idx,:]
end

function filter_by_max_action(data, max_nr_actions_per_subj)
    allSubj = sort(unique(data[:, Int(c_subject)]))
    (m,n) = size(data)
    allData = zeros(Int,(0,n))
    allSubj = sort(unique(data[:, Int(c_subject)]))
    for subj in allSubj
      data_subj = data[data[:,Int(c_subject)].==subj,:]
      (m,n) = size(data_subj)
      maxIdx = min(m, max_nr_actions_per_subj)
      data_subj = data_subj[1:maxIdx,:]
      allData = vcat(allData, data_subj)
    end
    allData
end

function load_SARSPEI(data_path; subject_id::Int = 0, max_episodes_per_subj::Int = typemax(Int), max_nr_actions_per_subj::Int = typemax(Int), sub_id_offset = 0)
    # sub_id_offset is used to concat multiple files where subject ID-n would reappaear
    data = collect(transpose(readdlm(data_path, ',', Int)))
    if sub_id_offset>0
        println("Warning: subject IDs offset by $sub_id_offset")
        data[:,Int(RL_Fit.c_subject)] += sub_id_offset
    end

    if subject_id != 0
        data=filter_by_person(data,subject_id)
    end
    data = filter_by_max_action(data, max_nr_actions_per_subj)
    data = filter_by_max_episode(data, max_episodes_per_subj)
    reward_rows = data[:,Int(c_reward)] .!=0 # use unit reward
    data[reward_rows,Int(c_reward)] .= 1
    stats = get_data_stats(data)
    data, stats
end

struct Data_stats
    nr_subjects::Int
    all_subject_ids::Array{Int,1}
    nr_puzzles::Int
    nr_episodes::Int
    nr_actions::Int
    max_state_id::Int
    max_action_id::Int
    nr_rewards::Int
end

function to_string(data_stats::Data_stats)
    #pretty print for console
    txt = string("nr_subjects: $(data_stats.nr_subjects)\n",
    "all_subject_ids: $(data_stats.all_subject_ids)\n",
    "nr_puzzles: $(data_stats.nr_puzzles)\n",
    "nr_episodes: $(data_stats.nr_episodes)\n",
    "nr_actions: $(data_stats.nr_actions)\n",
    "max_state_id: $(data_stats.max_state_id)\n",
    "max_action_id: $(data_stats.max_action_id)")
end

function to_short_string(data_stats::Data_stats)
    "S:$(data_stats.nr_subjects)|P:$(data_stats.nr_puzzles)|E:$(data_stats.nr_episodes)|A:$(data_stats.nr_actions)"
end

function get_data_stats(data::Array{Int,2})
    allSubj = sort(unique(data[:, Int(c_subject)]))
    nr_subjects = length(allSubj)
    nr_actions = size(data)[1]
    #count all episodes per subject
    nr_puzzles = 0
    nr_episodes = 0
    for s = allSubj
        idx = data[:,Int(c_subject)].== s
        nr_puzzles += length(unique(data[idx,Int(c_puzzle)]))
        nr_episodes += length(unique(data[idx,Int(c_episode)]))
    end
    #nr of different states and actions
    diff_states = unique( vcat( data[:,Int(c_state)], data[:,Int(c_nextState)] ))

    max_state_id = maximum(diff_states)
    max_action_id = maximum(data[:,Int(c_action)])
    nr_rewards = sum(data[:,Int(c_reward)] .!= 0)
    stats = Data_stats(nr_subjects, allSubj, nr_puzzles, nr_episodes, nr_actions, max_state_id, max_action_id, nr_rewards)
end

struct Step_by_step_signals
    subject_id::Int
    RPE::Array{Float64,1}
    SPE::Array{Float64,1}
    action_selection_prob::Array{Float64,1}
    LL::Float64
    data_stats::Data_stats
    algorithm_name::String
    param_vector::Array{Float64,1}
    param_string::String
    max_bias_action_index::Int
    results_dict::Dict{String,Any}  # optionally add more algo-specific results
end

function to_string(signals::Step_by_step_signals; round_nr_digits=3)
    txt = string("subject_id: $(signals.subject_id)\n",
      "algorithm_name: $(signals.algorithm_name)\n",
      "param_string: $(signals.param_string)\n",
      "param_vector: $(signals.param_vector)\n",
      "LL: $(round(signals.LL,digits=round_nr_digits))\n",
      "RPE\n$(round.(signals.RPE,digits=round_nr_digits))\n",
      "SPE\n$(round.(signals.SPE,digits=round_nr_digits))\n",
      "action_selection_prob\n$(round.(signals.action_selection_prob,digits=round_nr_digits))\n",
      "results_dict\n$(signals.results_dict)")
end

ext: jl


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediapackage_vod

using AWS.Compat
using AWS.UUIDs
"""
    CreateAsset()

Creates a new MediaPackage VOD Asset resource.

# Required Parameters
- `id`: The unique identifier for the Asset.
- `packagingGroupId`: The ID of the PackagingGroup for the Asset.
- `sourceArn`: ARN of the source object in S3.
- `sourceRoleArn`: The IAM role ARN used to access the source S3 bucket.

# Optional Parameters
- `resourceId`: The resource ID to include in SPEKE key requests.
- `tags`:
"""

create_asset(id, packagingGroupId, sourceArn, sourceRoleArn; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/assets", Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId, "sourceArn"=>sourceArn, "sourceRoleArn"=>sourceRoleArn); aws_config=aws_config)
create_asset(id, packagingGroupId, sourceArn, sourceRoleArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/assets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId, "sourceArn"=>sourceArn, "sourceRoleArn"=>sourceRoleArn), args)); aws_config=aws_config)

"""
    CreatePackagingConfiguration()

Creates a new MediaPackage VOD PackagingConfiguration resource.

# Required Parameters
- `id`: The ID of the PackagingConfiguration.
- `packagingGroupId`: The ID of a PackagingGroup.

# Optional Parameters
- `cmafPackage`:
- `dashPackage`:
- `hlsPackage`:
- `mssPackage`:
- `tags`:
"""

create_packaging_configuration(id, packagingGroupId; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_configurations", Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId); aws_config=aws_config)
create_packaging_configuration(id, packagingGroupId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_configurations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id, "packagingGroupId"=>packagingGroupId), args)); aws_config=aws_config)

"""
    CreatePackagingGroup()

Creates a new MediaPackage VOD PackagingGroup resource.

# Required Parameters
- `id`: The ID of the PackagingGroup.

# Optional Parameters
- `authorization`:
- `tags`:
"""

create_packaging_group(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_groups", Dict{String, Any}("id"=>id); aws_config=aws_config)
create_packaging_group(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/packaging_groups", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("id"=>id), args)); aws_config=aws_config)

"""
    DeleteAsset()

Deletes an existing MediaPackage VOD Asset resource.

# Required Parameters
- `id`: The ID of the MediaPackage VOD Asset resource to delete.

"""

delete_asset(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/assets/$(id)"; aws_config=aws_config)
delete_asset(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/assets/$(id)", args; aws_config=aws_config)

"""
    DeletePackagingConfiguration()

Deletes a MediaPackage VOD PackagingConfiguration resource.

# Required Parameters
- `id`: The ID of the MediaPackage VOD PackagingConfiguration resource to delete.

"""

delete_packaging_configuration(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_configurations/$(id)"; aws_config=aws_config)
delete_packaging_configuration(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_configurations/$(id)", args; aws_config=aws_config)

"""
    DeletePackagingGroup()

Deletes a MediaPackage VOD PackagingGroup resource.

# Required Parameters
- `id`: The ID of the MediaPackage VOD PackagingGroup resource to delete.

"""

delete_packaging_group(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_groups/$(id)"; aws_config=aws_config)
delete_packaging_group(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/packaging_groups/$(id)", args; aws_config=aws_config)

"""
    DescribeAsset()

Returns a description of a MediaPackage VOD Asset resource.

# Required Parameters
- `id`: The ID of an MediaPackage VOD Asset resource.

"""

describe_asset(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets/$(id)"; aws_config=aws_config)
describe_asset(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets/$(id)", args; aws_config=aws_config)

"""
    DescribePackagingConfiguration()

Returns a description of a MediaPackage VOD PackagingConfiguration resource.

# Required Parameters
- `id`: The ID of a MediaPackage VOD PackagingConfiguration resource.

"""

describe_packaging_configuration(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations/$(id)"; aws_config=aws_config)
describe_packaging_configuration(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations/$(id)", args; aws_config=aws_config)

"""
    DescribePackagingGroup()

Returns a description of a MediaPackage VOD PackagingGroup resource.

# Required Parameters
- `id`: The ID of a MediaPackage VOD PackagingGroup resource.

"""

describe_packaging_group(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups/$(id)"; aws_config=aws_config)
describe_packaging_group(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups/$(id)", args; aws_config=aws_config)

"""
    ListAssets()

Returns a collection of MediaPackage VOD Asset resources.

# Optional Parameters
- `maxResults`: Upper bound on number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
- `packagingGroupId`: Returns Assets associated with the specified PackagingGroup.
"""

list_assets(; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets"; aws_config=aws_config)
list_assets(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/assets", args; aws_config=aws_config)

"""
    ListPackagingConfigurations()

Returns a collection of MediaPackage VOD PackagingConfiguration resources.

# Optional Parameters
- `maxResults`: Upper bound on number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
- `packagingGroupId`: Returns MediaPackage VOD PackagingConfigurations associated with the specified PackagingGroup.
"""

list_packaging_configurations(; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations"; aws_config=aws_config)
list_packaging_configurations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_configurations", args; aws_config=aws_config)

"""
    ListPackagingGroups()

Returns a collection of MediaPackage VOD PackagingGroup resources.

# Optional Parameters
- `maxResults`: Upper bound on number of records to return.
- `nextToken`: A token used to resume pagination from the end of a previous request.
"""

list_packaging_groups(; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups"; aws_config=aws_config)
list_packaging_groups(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/packaging_groups", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of the tags assigned to the specified resource.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from the response to any request to the resource.

"""

list_tags_for_resource(resource_arn; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("GET", "/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    TagResource()

Adds tags to the specified resource. You can specify one or more tags to add.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from the response to any request to the resource.
- `tags`: A collection of tags associated with a resource

"""

tag_resource(resource_arn, tags; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes tags from the specified resource. You can specify one or more tags to remove.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) for the resource. You can get this from the response to any request to the resource.
- `tagKeys`: A comma-separated list of the tag keys to remove from the resource.

"""

untag_resource(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdatePackagingGroup()

Updates a specific packaging group. You can't change the id attribute or any other system-generated attributes.

# Required Parameters
- `id`: The ID of a MediaPackage VOD PackagingGroup resource.

# Optional Parameters
- `authorization`:
"""

update_packaging_group(id; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)"; aws_config=aws_config)
update_packaging_group(id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mediapackage_vod("PUT", "/packaging_groups/$(id)", args; aws_config=aws_config)


text:
mutable struct NodeTreeProblem
    pd::PopulationData
    bt::BinaryTree
    outgroupnode::Int
    nlevels::Int
    model::JuMP.Model
    assign::JuMP.JuMPArray{JuMP.Variable}
    assign2::JuMP.JuMPDict{JuMP.Variable}
    weight::JuMP.JuMPArray{JuMP.Variable}
    f3formula::JuMP.JuMPDict{JuMP.Variable}
    f3err::JuMP.JuMPDict{JuMP.Variable}
end

function NodeTreeProblem(
    pd::PopulationData,
    bt::BinaryTree;
    binaryencoding::Bool = false,
    nlevels::Int = 1,
    nmixtures::Int = pd.npop,
    solver = Gurobi.GurobiSolver())

    npop = pd.npop
    edges = bt.edges
    leaves = getleaves(bt)
    outgroupnode = leaves[1]
    othernodes = leaves[2:end]
    levels = 1:nlevels

    tree = JuMP.Model(solver=solver)
    if binaryencoding
        @assert nlevels == 1
        JuMP.@variable(tree, assign[1:npop,getleaves(bt),levels] >= 0)
        binaryencodingconstraints(pd, bt, tree, assign)
    else
        JuMP.@variable(tree, assign[1:npop,getleaves(bt),levels], Bin)
    end
    JuMP.@variable(tree, assign2[a=1:npop,b=a:npop,othernodes,othernodes,levels,levels])
    JuMP.@variable(tree, weight[edges] >= 0)
    JuMP.@variable(tree, f3formula[a=1:npop,b=a:npop,u=othernodes,v=othernodes,levels,levels] >= 0)
    JuMP.@variable(tree, f3err[a=1:npop,b=a:npop])

    validtreeconstraints(pd, bt, tree, assign, outgroupnode, nlevels, nmixtures)
    logicalconstraints(pd, bt, tree, assign, assign2, outgroupnode, nlevels)
    errorconstraints(pd, bt, tree, assign2, weight, f3formula, f3err, outgroupnode, nlevels)

    keytoind = Dict()
    ind = 0
    for key in keys(f3err)
        key[1] == pd.outgroup && continue
        key[2] == pd.outgroup && continue
        ind += 1
        keytoind[key] = ind
    end
    mat = zeros(ind,ind)
    for (a1,b1) in keys(keytoind), (a2,b2) in keys(keytoind)
        mat[keytoind[a1,b1],keytoind[a2,b2]] = pd.cov[a1,b1,a2,b2]
    end
    matinv = inv(mat)

    JuMP.@objective(tree, Min,
        sum(matinv[keytoind[a1,b1],keytoind[a2,b2]]*f3err[a1,b1]*f3err[a2,b2]
            for (a1,b1) in keys(keytoind), (a2,b2) in keys(keytoind)))

    NodeTreeProblem(pd, bt, outgroupnode, nlevels, tree, assign, assign2, weight, f3formula, f3err)
end

function logicalconstraints(
    pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    assign::JuMP.JuMPArray{JuMP.Variable},
    assign2::JuMP.JuMPDict{JuMP.Variable},
    outgroupnode::Int,
    nlevels::Int)

    npop = pd.npop
    othernodes = getleaves(bt)[2:end]
    levels = 1:nlevels

    JuMP.@constraint(tree, [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        assign2[a,b,u,v,l,m] <= assign[a,u,l])
    JuMP.@constraint(tree, [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        assign2[a,b,u,v,l,m] <= assign[b,v,m])
    JuMP.@constraint(tree, [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        assign2[a,b,u,v,l,m] >= assign[a,u,l] + assign[b,v,m] - 1)

end

function errorconstraints(
    pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    assign2::JuMP.JuMPDict{JuMP.Variable},
    weight::JuMP.JuMPArray{JuMP.Variable},
    f3formula::JuMP.JuMPDict{JuMP.Variable},
    f3err::JuMP.JuMPDict{JuMP.Variable},
    outgroupnode::Int,
    nlevels::Int)

    npop = pd.npop
    othernodes = getleaves(bt)[2:end]
    bigm = maximum(pd.f3)*2
    levels = 1:nlevels

    JuMP.@expression(tree,
        f3pathsum[u=othernodes,v=othernodes],
        sum(weight[edg] for edg in
            intersect(bt.pathedges[u,outgroupnode],
                      bt.pathedges[v,outgroupnode])))
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        f3formula[a,b,u,v,l,m] <= bigm*assign2[a,b,u,v,l,m])
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        f3formula[a,b,u,v,l,m] >= f3pathsum[u,v] + bigm*assign2[a,b,u,v,l,m] - bigm)
    JuMP.@constraint(tree,
        [a=1:npop,b=a:npop,u=othernodes,v=othernodes,l=levels,m=levels],
        f3formula[a,b,u,v,l,m] <= f3pathsum[u,v])
    JuMP.@constraint(tree, [a=1:npop,b=a:npop],
        f3err[a,b] == pd.f3[a,b] -
                      sum(f3formula[a,b,u,v,l,m]
                          for u=othernodes,v=othernodes,l=levels,m=levels)/nlevels^2)

end


text:
using GLMakie.Makie: Record
using GLMakie.GLFW
using GLMakie.ModernGL
using GLMakie.ShaderAbstractions
using GLMakie.ShaderAbstractions: Sampler
using GLMakie.StaticArrays
using GLMakie.GeometryBasics
using ReferenceTests.RNG

# A test case for wide lines and mitering at joints
@cell "Miter Joints for line rendering" begin
    scene = Scene()

    r = 4
    sep = 4*r
    scatter!(scene, (sep+2*r)*[-1,-1,1,1], (sep+2*r)*[-1,1,-1,1])

    for i=-1:1
        for j=-1:1
            angle = pi/2 + pi/4*i
            x = r*[-cos(angle/2),0,-cos(angle/2)]
            y = r*[-sin(angle/2),0,sin(angle/2)]

            linewidth = 40 * 2.0^j
            lines!(scene, x .+ sep*i, y .+ sep*j, color=RGBAf0(0,0,0,0.5), linewidth=linewidth)
            lines!(scene, x .+ sep*i, y .+ sep*j, color=:red)
        end
    end
    scene
end

@cell "Sampler type" begin
    # Directly access texture parameters:
    x = Sampler(fill(to_color(:yellow), 100, 100), minfilter=:nearest)
    scene = image(x, show_axis=false)
    # indexing will go straight to the GPU, while only transfering the changes
    st = Stepper(scene)
    x[1:10, 1:50] .= to_color(:red)
    Makie.step!(st)
    x[1:10, end] .= to_color(:green)
    Makie.step!(st)
    x[end, end] = to_color(:blue)
    Makie.step!(st)
    st
end
# Test for resizing of TextureBuffer
@cell "Dynamically adjusting number of particles in a meshscatter" begin

    pos = Node(RNG.rand(Point3f0, 2))
    rot = Node(RNG.rand(Vec3f0, 2))
    color = Node(RNG.rand(RGBf0, 2))
    size = Node(0.1*RNG.rand(2))

    makenew = Node(1)
    on(makenew) do i
        pos[] = RNG.rand(Point3f0, i)
        rot[] = RNG.rand(Vec3f0, i)
        color[] = RNG.rand(RGBf0, i)
        size[] = 0.1*RNG.rand(i)
    end

    scene = meshscatter(pos,
        rotations=rot,
        color=color,
        markersize=size,
        limits=FRect3D(Point3(0), Point3(1))
    )

    Record(scene, [10, 5, 100, 60, 177]) do i
        makenew[] = i
    end
end

@cell "Explicit frame rendering" begin
    set_window_config!(renderloop=(screen) -> nothing)
    function update_loop(m, buff, screen)
        for i = 1:20
            GLFW.PollEvents()
            buff .= RNG.rand.(Point3f0) .* 20f0
            m[1] = buff
            GLMakie.render_frame(screen)
            GLFW.SwapBuffers(GLMakie.to_native(screen))
            glFinish()
        end
    end
    fig, ax, meshplot = meshscatter(RNG.rand(Point3f0, 10^4) .* 20f0)
    screen = Makie.backend_display(GLMakie.GLBackend(), fig.scene)
    buff = RNG.rand(Point3f0, 10^4) .* 20f0;
    update_loop(meshplot, buff, screen)
    set_window_config!(renderloop=GLMakie.renderloop)
    fig
end


text:
# Discrete Measure
#
# Represents a generic discrete distribution with finitely many atoms

# Returns a weighted sample with one row for each distinct row in X and an
# associated weight equaling the sum of all weights.  Rows are sorted in
# ascending order if X is univariate.
function compresssample(X::AbstractArray{T}, q::AbstractArray{Float64}) where T
    n = size(X,1);
    p = size(X,2);
    # Add each row of X with its weight to an accumulator
    a = DataStructures.Accumulator{Array{T,2},Float64}()
    for i in 1:n
        push!(a,X[i:i,:],q[i])
    end
    # Collect distinct points
    nnew = length(a);
    Xnew = zeros(Float64,nnew,p);
    qnew = zeros(Float64,nnew,1);
    for (i,x) in enumerate(keys(a))
        Xnew[i,:] = x;
        qnew[i] = a[x];
    end
    # Sort points if univariate
    if p == 1
        order = sortperm(Xnew[:,1]);
        Xnew = Xnew[order,:];
        qnew = qnew[order,:];
    end
    (Xnew, qnew)
end

mutable struct SteinDiscrete
    # n x p matrix of support points
    support::Array{Float64,2}
    # n x 1 vector of weights associated with support points
    weights::Array{Float64,2}
    # Internal constructor collapses X to its unique sorted elements and
    # adjusts q accordingly; sorts X if univariate.
    SteinDiscrete(support, weights) = (
        (support_,weights_) = compresssample(support,weights);
        return new(support_, weights_);
    )
end

# Constructor with uniform weights q_i = 1/n
SteinDiscrete(X) = SteinDiscrete(
    X,
    [1/size(X,1) for i=1:size(X,1)]
);


text:
"""
  `EVENT_STATE`

Enum with values:

- `idle=0`
- `triggered=1`
- `processed=2`
"""
@enum EVENT_STATE idle=0 triggered=1 processed=2

"""
  `Event`

An event is a state machine with three states:

- `idle`
- `triggered`
- `processed`

An event is initially not triggered. Events get trigerred after they are scheduled for processing.

An event has a list of callbacks and a value. A callback can be any function. Once an event gets processed, all callbacks will be invoked. Callbacks can do further processing with the value it has produced.

Failed events, i.e. events having as value an `Exception`, are never silently ignored and will raise this exception upon being processed.

**Fields:**

- `cid :: UInt`
- `callbacks :: PriorityQueue{Function, UInt}`
- `state :: EVENT_STATE`
- `value :: Any`

**Constructor:**

`Event()`
"""
type Event
  cid :: UInt
  callbacks :: PriorityQueue{Function, UInt}
  state :: EVENT_STATE
  value :: Any
  function Event()
    ev = new()
    ev.cid = zero(UInt)
    ev.callbacks = PriorityQueue(Function, UInt)
    ev.state = idle
    ev.value = nothing
    return ev
  end
end

type Process
  task :: Task
  target :: Event
  resume :: Function
  ev :: Event
  function Process(task::Task, target::Event)
    proc = new()
    proc.task = task
    proc.target = target
    proc.resume = append_callback(proc.target, execute, proc)
    proc.ev = Event()
    return proc
  end
end

"""
  `EventKey`

Key for the event heap.

**Fields**:

- `time :: Float64`
- `priority :: Bool`
- `id :: UInt`

**Constructor**:

`EventKey(time :: Float64, priority :: Bool, id :: UInt`)

Only used internally.
"""
immutable EventKey
  time :: TimeType
  priority :: Bool
  id :: UInt
end

"""
  `isless(a::EventKey, b::EventKey) :: Bool`

Compairs two `EventKey`. The criteria in order of importance are:

- time of processing
- priority when time of processing is equal
- scheduling id, i.e. the event that was first scheduled is first processed when time of processing and priority are identical

Only used internally.
"""
function isless(a::EventKey, b::EventKey) :: Bool
  (a.time < b.time) || (a.time == b.time && a.priority > b.priority) || (a.time == b.time && a.priority == b.priority && a.id < b.id)
end

"""
  `Simulation{T<:TimeType}`

Execution environment for a simulation. The passing of time is implemented by stepping from event to event.

**Fields**:

- `time :: T`
- `heap :: PriorityQueue{Event, EventKey}`
- `sid :: UInt`
- `active_proc :: Nullable{Process}`
- `granularity` :: Type

**Constructor**:

`Simulation{T<:TimeType}(initial_time::T)`
`Simulation(initial_time::Number=0)`

An initial_time for the simulation can be specified. By default, it starts at 0.
"""
type Simulation{T<:TimeType}
  time :: T
  heap :: PriorityQueue{Event, EventKey}
  sid :: UInt
  active_proc :: Nullable{Process}
  granularity :: Type
  function Simulation(initial_time::T)
    sim = new()
    sim.time = initial_time
    sim.heap = PriorityQueue(Event, EventKey)
    sim.sid = zero(UInt)
    sim.active_proc = Nullable{Process}()
    sim.granularity = typeof(initial_time.instant.periods)
    return sim
  end
end

function Simulation{T<:TimeType}(initial_time::T)
  Simulation{T}(initial_time)
end

function Simulation(initial_time::Number=0)
  Simulation(SimulationTime(initial_time))
end


text:
using Base.Test

include("atbash-cipher.jl")

@testset "encoding from English to atbash" begin
    @testset "encode yes" begin
        @test encode("yes") == "bvh"
    end

    @testset "encode no" begin
        @test encode("no") == "ml"
    end

    @testset "encode OMG" begin
        @test encode("OMG") == "lnt"
    end

    @testset "encode spaces" begin
        @test encode("O M G") == "lnt"
    end

    @testset "encode mindblowingly" begin
        @test encode("mindblowingly") == "nrmwy oldrm tob"
    end

    @testset "encode numbers" begin
        @test encode("Testing,1 2 3, testing.") == "gvhgr mt123 gvhgr mt"
    end

    @testset "encode deep thought" begin
        @test encode("Truth is fiction.") == "gifgs rhurx grlm"
    end

    @testset "encode all the letters" begin
        @test encode("The quick brown fox jumps over the lazy dog.") == "gsvjf rxpyi ldmul cqfnk hlevi gsvoz abwlt"
    end
end

@testset "decoding from atbash to English" begin
    @testset "decode exercism" begin
        @test decode("vcvix rhn") == "exercism"
    end

    @testset "decode a sentence" begin
        @test decode("zmlyh gzxov rhlug vmzhg vkkrm thglm v") == "anobstacleisoftenasteppingstone"
    end

    @testset "decode numbers" begin
        @test decode("gvhgr mt123 gvhgr mt") == "testing123testing"
    end

    @testset "decode all the letters" begin
        @test decode("gsvjf rxpyi ldmul cqfnk hlevi gsvoz abwlt") == "thequickbrownfoxjumpsoverthelazydog"
    end
end


text:
export ft_setrts

function ft_setrts(ft_handle::Culong)
  ft_status = ccall((:FT_SetRts, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
using Documenter, PhysicalConstant
using FileIO

makedocs(
    modules = [PhysicalConstant],
    clean = false,
    format = Documenter.HTML(),
    sitename = "PhysicalConstant",
    authors = "La Guer.",
    pages    = Any[
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Usage"        => "usage.md",
        "Examples"     => "examples.md",
        "api"          => "api.md",
    ]
)

deploydocs(
    julia = "nightly",
    repo = "github.com/LaGuer/PhysicalConstant.jl.git",
    target = "build",
    deps = Deps.pip("Tornado>=4.0.0,<5.0.0", "mkdocs==0.17.5", "mkdocs-material==2.9.4", "python-markdown-math"),
    make = nothing,
)


text:
mutable struct WorldState
    max_days::Int
    current_state::Array
end

include("Rules.jl")

function WorldState(;
    max_days::Int = 100,
    initial_state::Array = default_state(),
)
    world_size = size(initial_state)
    if length(world_size) != 2
        error("The world is not 2 dimensional")
    elseif world_size[1] < 3 || world_size[2] < 3
        error("The world is too small (less than 3x3)")
    elseif state_check(initial_state) == false
        error("The cells do not have a binary state (0 or 1)")
    end
    WorldState(max_days, initial_state)
end

function state_check(initial_state)
    for i in initial_state
        if i ∉ [0, 1]
            return false
        end
    end
    return true
end

function default_state()
    A = fill(0, 128, 128)
    A[61:63, 61:63] = [
        1 1 0
        0 1 1
        0 1 0
    ]
    return A
end

Base.length(world::WorldState) = world.max_days

function Base.iterate(world::WorldState, state = 0)
    if state > world.max_days
        return nothing
    end
    world.current_state = next(world.current_state)
    world, state + 1
end

max_stars_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_issues_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_forks_repo_name: VasanthManiVasi/Conways-Game-of-Life


text:
using ArgParse
using TMLEEpistasis

function parse_commandline()
    s = ArgParseSettings(
        description = "This program computes estimates of Genetic Epistatis from the UK-Biobank using TMLE."*
                      "Here is a list of the arguments that should be provided, you can also have a look at the "*
                      "test/data and test/config folders to see some examples.",
        commands_are_required = false,
        version = "0.2",
        add_version = true)

    @add_arg_table s begin
        "phenotypes"
            help = "A file (.csv format). The first row contains the column names with `eid` the sample ID"*
                   " and the rest of the columns are phenotypes of interest."
            required = true
        "confounders"
            help = "A file (.csv format) containing the confounding variables values and the sample ids associated"*
                   " with the participants. The first line of the file should contain the columns names and the sample ids "*
                   " column name should be: `SAMPLE_ID`."
            required = true
        "queries"
            help = "A file (.toml format) see: config/sample_query.toml for more information"
            required = true
        "estimator"
            help = "A file (.toml format) describing the tmle estimator to use, see config/sample_estimator.toml"*
                   " for a basic example."
            required = true
        "output"
            help = "A path where the results will be saved (.csv format)"
            required = true
        "--phenotypes-list", "-p"
            help = "A file, one line for each phenotype, containing a restrictions of the phenotypes "*
                   "to consider for the analysis."
            required = false
            arg_type = String
        "--verbosity", "-v"
            help = "Verbosity level"
            arg_type = Int
            default = 1
        "epistasis"
            help = "Run the Epistasis effect size estimation procedure"
            action = :command
        "crossval"
            help = "Run the Cross validation procedure"
            action = :command
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

if parsed_args["%COMMAND%"] == "crossval"
    UKBBVariantRun(parsed_args, run_fn=TMLEEpistasis.PhenotypeCrossValidation)
elseif parsed_args["%COMMAND%"] == "epistasis"
    UKBBVariantRun(parsed_args, run_fn=TMLEEpistasis.PhenotypeTMLEEpistasis)
end


text:
#An implementation of the Enigma I machine in Julia

#This function ensures an index is between 1 and 26, and if not loops it around
function loopIndex(i)
    if(i > 26)
        i -= 26
    end
    if(i < 1)
        i += 26
    end
    return i
end

#This function ensures a character is in the range A - Z, and if not loops it around
function loopChar(c)
    if(c > 'Z')
        c -= 26
    end
    if(c < 'A')
        c += 26
    end
    return c
end

#This method increments a list of characters according to a given int
function incrementMappings(s, i)
    out = ""
    for c in s
        c = loopChar(c + i)
        out = string(out, c)
    end
    return out
end

# This function gets the correct rotor and step given a Key
function getRotor(i)
    #= In this model, only rotors I, II and III are available
    These mappings are the same as the original machines. =#
    rotor1 = "EKMFLGDQVZNTOWYHXUSPAIBRCJ"
    rotor1step = 'Q'
    rotor2 = "AJDKSIRUXBLHWTMCQGZNPYFVOE"
    rotor2step = 'E'
    rotor3 = "BDFHJLCPRTXVZNYEIWGAKMUSQO"
    rotor3step = 'V'
    if i == '1'
        return rotor1, rotor1step
    elseif i == '2'
        return rotor2, rotor2step
    elseif i == '3'
        return rotor3, rotor3step
    end
end

# This function gets the correct reflector. It defaults to UKW - B
function getReflector(i)
    reflectorA = "EJMZALYXVBWFCRQUONTSPIKHGD"
    reflectorB = "YRUHQSLDPXNGOKMIEBFZCWVJAT"
    reflectorC = "FVPJIAOYEDRZXWGCTKUQSBNMHL"
    if i == "A"
        return reflectorA
    elseif i == "C"
        return reflectorC
    else
        return reflectorB
    end
end

#The main function
function enigma()

    # Get plaintext
    println("Please enter the text: ")
    text = readline()
    text = uppercase(text)

    # Select the orders of rotors used
    println("Please enter the rotor order in the form XXX")
    rotorOrder = readline()
    if length(rotorOrder) == 0
        firstRotor, firstRotorStep = getRotor('3')
        secondRotor, secondRotorStep = getRotor('2')
        thirdRotor,thirdRotorStep = getRotor('1')
    else
        firstRotor, firstRotorStep = getRotor(rotorOrder[1])
        secondRotor, secondRotorStep = getRotor(rotorOrder[2])
        thirdRotor,thirdRotorStep = getRotor(rotorOrder[3])
    end

    # The user can set the rotor starting positions
    println("Enter the starting rotor positions in the form XXX, or press enter to use the default AAA")
    startPos = readline()
    if(length(startPos) == 0)
        increment1 = 0
        increment2 = 0
        increment3 = 0
    else
        increment1 = startPos[3] - 'A'
        increment2 = startPos[2] - 'A'
        increment3 = startPos[1] - 'A'
    end

    #= The standard reflector used in the war was reflector UKW-B, so I use it as the default
    The user can also choose UKW-A or UKW-C if they wish =#
    println("Please specify the reflector (A for UKW-A, B for UKW-B or C for UKW-C) or press enter to use the default UKW-B")
    reflectorChoice = readline()
    if length(reflectorChoice) == 0
        reflector = getReflector('B')
    else
        reflector = getReflector(reflectorChoice)
    end

    #=To start, I will also assume fixed ring settings
    TODO: implement ring settings =#
    println("DON'T USE Enter the ring positions in the form XXX, or press enter to use the default AAA")
    ringSettings = readline()
    if(length(ringSettings) == 0)
        ring1 = 0
        ring2 = 0
        ring3 = 0
    else
        ring1 = ringSettings[3] - 'A'
        #increment1 -= ring1
        ring2 = ringSettings[2] - 'A'
        #increment2 -= ring2
        ring3 = ringSettings[1] - 'A'
        #increment3 -= ring3
    end

    # The mappings need to be all incremented according to the ring settings
    firstRotor = incrementMappings(firstRotor, ring1)
    secondRotor = incrementMappings(secondRotor, ring2)
    thirdRotor = incrementMappings(thirdRotor, ring3)

    #Plug board settings (Steckerbrett) directly map given letters at the very beginning and very end of operation
    println("Please Insert plug board settings in the form AB CD EF etc. or press enter to use no plug board")
    plugBoardInput = readline()
    plugBoardInput = uppercase(plugBoardInput)
    plugBoard = split(plugBoardInput)

    #=Now it is time to run the cipher
    The input first goes through the plug board (if used)
    It then goes R1 -> R2 -> R3 -> REF -> R3' -> R2' -> R1'
    Where R(x) denotes the rotor at a given position, ' denotes an inverted rotor and REF denotes a reflector
    The output finally goes though the plug board again (if used) =#
    for char in text

        #Step rotors
        increment1 += 1

        # Plug board
        for s in plugBoard
            if s[1] == char
                char = s[2]
            end
        end

        # First rotor
        index = char - 'A' + increment1 + 1
        index = loopIndex(index)
        char1 = loopChar(firstRotor[index] - increment1)

        # Second rotor
        index = char1 - 'A' + increment2 + 1
        index = loopIndex(index)
        char2 = loopChar(secondRotor[index] - increment2)

        # Third rotor
        index = char2 - 'A' + increment3 + 1
        index = loopIndex(index)
        char3 = loopChar(thirdRotor[index] - increment3)

        # Reflector
        index = char3 - 'A' + 1
        index = loopIndex(index)
        refChar = reflector[index]

        # Third rotor reflected (takes the index from the rotor, and works backwards to the normal alphabet)
        #The increments much be taken into account in reverse
        val = refChar + increment3
        # It must be insured any character used is alphabetical
        val = loopChar(val)
        index = findfirst(isequal(val), thirdRotor)  - increment3
        index = loopIndex(index)
        char4 = loopChar('A' + index - 1)

        # Second rotor reflected (takes the index from the rotor, and works backwards to the normal alphabet)
        val = char4 + increment2
        val = loopChar(val)
        index = findfirst(isequal(val), secondRotor)  - increment2
        index = loopIndex(index)
        char5 = loopChar('A' + index - 1)

        # First rotor reflected (takes the index from the rotor, and works backwards to the normal alphabet)
        val = char5 + increment1
        val = loopChar(val)
        index = findfirst(isequal(val), firstRotor)  - increment1
        index = loopIndex(index)
        char6 = loopChar('A' + index - 1)

        # Plug board
        for s in plugBoard
            if s[2] == char6
                char6 = s[1]
            end
        end

        print(char6)
    end
end
enigma()


text:
"""
    unwrap_nested_exception(err::Exception) -> Bool

Extracts the "core" exception from a nested exception."
"""
unwrap_nested_exception(err::CapturedException) =
    unwrap_nested_exception(err.ex)
unwrap_nested_exception(err::RemoteException) =
    unwrap_nested_exception(err.captured)
unwrap_nested_exception(err) = err

"Fills the result for all registered futures of `node`."
function fill_registered_futures!(state, node, failed)
    if haskey(state.futures, node)
        # Notify any listening thunks
        for future in state.futures[node]
            put!(future, state.cache[node]; error=failed)
        end
        delete!(state.futures, node)
    end
end

"Cleans up any inputs that aren't needed any longer, and returns a `Set{Chunk}`
of all chunks that can now be evicted from workers."
function cleanup_inputs!(state, node)
    to_evict = Set{Chunk}()
    for inp in unwrap_weak_checked.(node.inputs)
        if !istask(inp) && !(inp isa Chunk)
            continue
        end
        if inp in keys(state.waiting_data)
            w = state.waiting_data[inp]
            if node in w
                pop!(w, node)
            end
            if isempty(w)
                if istask(inp) && haskey(state.cache, inp)
                    _node = state.cache[inp]
                    if _node isa Chunk
                        push!(to_evict, _node)
                    end
                elseif inp isa Chunk
                    push!(to_evict, inp)
                end
                delete!(state.waiting_data, inp)
            end
        end
    end
    to_evict
end

"Schedules any dependents that may be ready to execute."
function schedule_dependents!(state, node, failed)
    for dep in sort!(collect(get(()->Set{Thunk}(), state.waiting_data, node)), by=state.node_order)
        dep_isready = false
        if haskey(state.waiting, dep)
            set = state.waiting[dep]
            node in set && pop!(set, node)
            dep_isready = isempty(set)
            if dep_isready
                delete!(state.waiting, dep)
            end
        else
            dep_isready = true
        end
        if dep_isready
            if !failed
                push!(state.ready, dep)
            end
        end
    end
end

"""
Prepares the scheduler to schedule `thunk`, including scheduling `thunk` if
its inputs are satisfied.
"""
function reschedule_inputs!(state, thunk, seen=Set{Thunk}())
    thunk in seen && return
    push!(seen, thunk)
    if haskey(state.cache, thunk) || (thunk in state.ready) || (thunk in state.running)
        return
    end
    w = get!(()->Set{Thunk}(), state.waiting, thunk)
    for input in thunk.inputs
        input = unwrap_weak_checked(input)
        if istask(input) || (input isa Chunk)
            push!(get!(()->Set{Thunk}(), state.waiting_data, input), thunk)
        end
        istask(input) || continue
        if get(state.errored, input, false)
            set_failed!(state, input, thunk)
        end
        haskey(state.cache, input) && continue
        push!(w, input)
        if !((input in state.running) || (input in state.ready))
            reschedule_inputs!(state, input, seen)
        end
    end
    if isempty(w)
        # Inputs are ready
        delete!(state.waiting, thunk)
        if !get(state.errored, thunk, false)
            push!(state.ready, thunk)
        end
    end
end

"Marks `thunk` and all dependent thunks as failed."
function set_failed!(state, origin, thunk=origin)
    filter!(x->x!==thunk, state.ready)
    state.cache[thunk] = ThunkFailedException(thunk, origin, state.cache[origin])
    state.errored[thunk] = true
    fill_registered_futures!(state, thunk, true)
    if haskey(state.waiting_data, thunk)
        for dep in state.waiting_data[thunk]
            haskey(state.waiting, dep) &&
                delete!(state.waiting, dep)
            haskey(state.errored, dep) &&
                continue
            set_failed!(state, origin, dep)
        end
        delete!(state.waiting_data, thunk)
    end
    if haskey(state.waiting, thunk)
        delete!(state.waiting, thunk)
    end
end

"Internal utility, useful for debugging scheduler state."
function print_sch_status(state, thunk; kwargs...)
    iob = IOBuffer()
    print_sch_status(iob, state, thunk; kwargs...)
    seek(iob, 0)
    write(stderr, iob)
end
function print_sch_status(io::IO, state, thunk; offset=0, limit=5, max_inputs=3)
    function status_string(node)
        status = ""
        if get(state.errored, node, false)
            status *= "E"
        end
        if node in state.ready
            status *= "r"
        elseif node in state.running
            status *= "R"
        elseif haskey(state.cache, node)
            status *= "C"
        else
            status *= "?"
        end
        status
    end
    if offset == 0
        println(io, "Ready ($(length(state.ready))): $(join(map(t->t.id, state.ready), ','))")
        println(io, "Running: ($(length(state.running))): $(join(map(t->t.id, collect(state.running)), ','))")
        print(io, "($(status_string(thunk))) ")
    end
    println(io, "$(thunk.id): $(thunk.f)")
    for (idx,input) in enumerate(thunk.inputs)
        if input isa WeakThunk
            input = unwrap_weak(input)
            if input === nothing
                println(io, repeat(' ', offset+2), "(???)")
                continue
            end
        end
        input isa Thunk || continue
        if idx > max_inputs
            println(io, repeat(' ', offset+2), "…")
            break
        end
        status = status_string(input)
        if haskey(state.waiting, thunk) && input in state.waiting[thunk]
            status *= "W"
        end
        if haskey(state.waiting_data, input) && thunk in state.waiting_data[input]
            status *= "w"
        end
        if haskey(state.futures, input)
            status *= "f($(length(state.futures[input])))"
        end
        print(io, repeat(' ', offset+2), "($status) ")
        if limit > 0
            print_sch_status(io, state, input; offset=offset+2, limit=limit-1)
        else
            println(io, "…")
        end
    end
end

function fetch_report(task)
    try
        fetch(task)
    catch err
        @static if VERSION >= v"1.1"
            @static if VERSION < v"1.7-rc1"
                stk = Base.catch_stack(task)
            else
                stk = Base.current_exceptions(task)
            end
            err, frames = stk[1]
            rethrow(CapturedException(err, frames))
        else
            rethrow(task.result)
        end
    end
end

fn_type(x::Chunk) = x.chunktype
fn_type(x) = typeof(x)
function signature(task::Thunk, state)
    inputs = map(x->istask(x) ? state.cache[x] : x, unwrap_weak_checked.(task.inputs))
    Tuple{fn_type(task.f), map(x->x isa Chunk ? x.chunktype : typeof(x), inputs)...}
end

function report_catch_error(err, desc=nothing)
    iob = IOContext(IOBuffer(), :color=>true)
    if desc !== nothing
        println(iob, desc)
    end
    Base.showerror(iob, err)
    Base.show_backtrace(iob, catch_backtrace())
    println(iob)
    seek(iob.io, 0)
    write(stderr, iob)
end


text:
using Test
include( "../src/predictEnsemble.jl" )
println( "running predictEnsemble unit tests." );

predictions = Vector{Matrix{Float64}}();
pred_1      = ones( 100, 3 );
pred_2      = 2.0 * ones( 100, 3 );
pred_3      = 3.0 * ones( 100, 3 );
push!( predictions, pred_1 );
push!( predictions, pred_2 );
push!( predictions, pred_3 );
weights     = zeros( 3 );
fill!( weights, 1/3 );
yEns  = predictEnsemble( predictions, weights );
@test all( yEns .≈ 2.0 )


text:
@testset "conversions.jl" begin
    @testset "Makie" begin
        @testset "frequencies" begin
            probabilities = rand(100)
            frequencies = rand(100)

            # without consistency bounds
            points = convert_arguments(Reliability, probabilities, frequencies)
            @test length(points) == 1
            @test points[1] isa Vector{Point4f0}
            @test map(first, points[1]) ≈ probabilities
            @test map(x -> x[2], points[1]) ≈ frequencies
            @test all(isnan(x[3]) for x in points[1])
            @test all(isnan(x[4]) for x in points[1])

            # with consistency bounds
            low = max.(0, probabilities .- rand.())
            high = min.(1, probabilities .+ rand.())
            points = convert_arguments(Reliability, probabilities, frequencies, low, high)
            @test length(points) == 1
            @test points[1] isa Vector{Point4f0}
            @test map(first, points[1]) ≈ probabilities
            @test map(x -> x[2], points[1]) ≈ frequencies
            @test map(x -> x[3], points[1]) ≈ low
            @test map(x -> x[4], points[1]) ≈ high

            low_high = map(tuple, low, high)
            points = convert_arguments(Reliability, probabilities, frequencies, low_high)
            @test length(points) == 1
            @test points[1] isa Vector{Point4f0}
            @test map(first, points[1]) ≈ probabilities
            @test map(x -> x[2], points[1]) ≈ frequencies
            @test map(x -> x[3], points[1]) ≈ low
            @test map(x -> x[4], points[1]) ≈ high
        end

        @testset "outcomes" begin
            probabilities = rand(100)
            outcomes = rand(Bool, 100)

            for binning in (EqualSize(), EqualMass()),
                consistencybars in (ConsistencyBars(), nothing)

                points = convert_arguments(
                    Reliability,
                    probabilities,
                    outcomes;
                    binning=binning,
                    consistencybars=consistencybars,
                )
                @test length(points) == 1
                @test points[1] isa Vector{Point4f0}
                @test length(points[1]) == binning.n

                meanprobs, meanfreqs, _ = ReliabilityDiagrams.means_and_bars(
                    probabilities, outcomes; binning=binning, consistencybars=nothing
                )
                @test map(first, points[1]) ≈ meanprobs
                @test map(x -> x[2], points[1]) ≈ meanfreqs
                if consistencybars === nothing
                    @test all(isnan(x[3]) for x in points[1])
                    @test all(isnan(x[4]) for x in points[1])
                else
                    @test all(x[3] <= x[4] for x in points[1])
                end
            end
        end
    end

    @testset "Plots" begin
        @testset "frequencies" begin
            probabilities = rand(100)
            frequencies = rand(100)

            # without consistency bounds
            plt = ReliabilityPlot((probabilities, frequencies))
            @test plt isa ReliabilityPlot
            @test plt.x == probabilities
            @test plt.y == frequencies
            @test plt.low_high === nothing

            # with consistency bounds
            low = max.(0, probabilities .- rand.())
            high = min.(1, probabilities .+ rand.())
            low_high = map(tuple, low, high)

            plt = ReliabilityPlot((probabilities, frequencies, low, high))
            @test plt isa ReliabilityPlot
            @test plt.x == probabilities
            @test plt.y == frequencies
            @test plt.low_high == low_high

            plt = ReliabilityPlot((probabilities, frequencies, low_high))
            @test plt isa ReliabilityPlot
            @test plt.x == probabilities
            @test plt.y == frequencies
            @test plt.low_high == low_high
        end

        @testset "outcomes" begin
            probabilities = rand(100)
            outcomes = rand(Bool, 100)
            plt = ReliabilityPlot((probabilities, outcomes))
            @test plt isa ReliabilityPlot
            @test plt.x == probabilities
            @test plt.y == outcomes
            @test plt.low_high === nothing
        end
    end
end


text:
using Sim6502

cs = Chipset()

program = @assemble begin
    LDA, 0x01
    STA, Π(0x0200)
    LDA, 0x05
    STA, Π(0x0201)
    LDA, 0x08
    STA, Π(0x0202)
    BRK
end


text:
#########
## map ##
#########

# Single input
@generated function map{T}(f, a1::StaticArray{T})
    newtype = :(similar_type($a1, promote_op(f, T)))
    exprs = [:(f(a1[$j])) for j = 1:length(a1)]
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# Two inputs
@generated function map{T1,T2}(f, a1::StaticArray{T1}, a2::StaticArray{T2})
    if size(a1) != size(a2)
        error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))")
    end

    newtype = :(similar_type($a1, promote_op(f, T1, T2)))
    exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a1)]
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# TODO these assume linear fast...
@generated function map{T1,T2}(f, a1::StaticArray{T1}, a2::AbstractArray{T2})
    newtype = :(similar_type($a1, promote_op(f, T1, T2)))
    exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a1)]
    return quote
        $(Expr(:meta, :inline))

        if size(a1) != size(a2)
            error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))")
        end

        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function map{T1,T2}(f, a1::AbstractArray{T1}, a2::StaticArray{T2})
    newtype = :(similar_type($a2, promote_op(f, T1, T2)))
    exprs = [:(f(a1[$j], a2[$j])) for j = 1:length(a2)]
    return quote
        $(Expr(:meta, :inline))

        @boundscheck if size(a1) != size(a2)
            error("Dimensions must match. Got sizes $(size(a1)) and $(size(a2))")
        end

        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

# TODO General case involving arbitrary many inputs?

############
## reduce ##
############
@generated function reduce(op, a1::StaticArray)
    if length(a1) == 1
        return :(a1[1])
    else
        expr = :(op(a1[1], a1[2]))
        for j = 3:length(a1)
            expr = :(op($expr, a1[$j]))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

@generated function reduce(op, v0, a1::StaticArray)
    if length(a1) == 0
        return :(v0)
    else
        expr = :(op(v0, a1[1]))
        for j = 2:length(a1)
            expr = :(op($expr, a1[$j]))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

###############
## mapreduce ##
###############

# Single array
@generated function mapreduce(f, op, a1::StaticArray)
    if length(a1) == 1
        return :(f(a1[1]))
    else
        expr = :(op(f(a1[1]), f(a1[2])))
        for j = 3:length(a1)
            expr = :(op($expr, f(a1[$j])))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

@generated function mapreduce(f, op, v0, a1::StaticArray)
    if length(a1) == 0
        return :(v0)
    else
        expr = :(op(v0, f(a1[1])))
        for j = 2:length(a1)
            expr = :(op($expr, f(a1[$j])))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

# Two arrays (e.g. dot has f(a,b) = a' * b, op = +)
@generated function mapreduce(f, op, a1::StaticArray, a2::StaticArray)
    if size(a1) != size(a2)
        error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))")
    end

    if length(a1) == 1
        return :(f(a1[1], a2[1]))
    else
        expr = :(op(f(a1[1], a2[1]), f(a1[2], a2[2])))
        for j = 3:length(a1)
            expr = :(op($expr, f(a1[$j], a2[$j])))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

@generated function mapreduce(f, op, v0, a1::StaticArray, a2::StaticArray)
    if size(a1) != size(a2)
        error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))")
    end

    if length(a1) == 0
        return :(v0)
    else
        expr = :(op(v0, f(a1[1], a2[1])))
        for j = 2:length(a1)
            expr = :(op($expr, f(a1[$j], a2[$j])))
        end
        return quote
            $(Expr(:meta, :inline))
            $expr
        end
    end
end

# TODO General case involving arbitrary many inputs?

###############
## broadcast ##
###############
# Single input version
@inline broadcast(f, a::StaticArray) = map(f, a)

# Two input versions
@generated function broadcast(f, a1::StaticArray, a2::StaticArray)
    if size(a1) == size(a2)
        return quote
            $(Expr(:meta, :inline))
            map(f, a1, a2)
        end
    else
        s1 = size(a1)
        s2 = size(a2)
        ndims = max(length(s1), length(s2))

        s = Vector{Int}(ndims)
        expands1 = Vector{Bool}(ndims)
        expands2 = Vector{Bool}(ndims)
        for i = 1:ndims
            if length(s1) < i
                s[i] = s2[i]
                expands1[i] = false
                expands2[i] = s2[i] > 1
            elseif length(s2) < i
                s[i] = s1[i]
                expands1[i] = s1[i] > 1
                expands2[i] = false
            else
                s[i] = max(s1[i], s1[i])
                @assert s1[i] == 1 || s1[i] == s[i]
                @assert s2[i] == 1 || s2[i] == s[i]
                expands1[i] = s1[i] > 1
                expands2[i] = s2[i] > 1
            end
        end
        s = (s...)
        L = prod(s)

        if s == s1
            newtype = :( similar_type($a1, promote_op(f, $(eltype(a1)), $(eltype(a2)))) )
        else
            newtype = :( similar_type($a1, promote_op(f, $(eltype(a1)), $(eltype(a2))), $s) )
        end

        exprs = Vector{Expr}(L)

        i = 1
        ind = ones(Int, ndims)
        while i <= L
            ind1 = [expands1[j] ? ind[j] : 1 for j = 1:length(s1)]
            ind2 = [expands2[j] ? ind[j] : 1 for j = 1:length(s2)]

            exprs[i] = Expr(:call, :f, Expr(:ref, :a1, ind1...), Expr(:ref, :a2, ind2...))

            i += 1

            ind[1] += 1
            j = 1
            while j < length(s)
                if ind[j] > s[j]
                    ind[j] = 1
                    ind[j+1] += 1
                else
                    break
                end

                j += 1
            end
        end

        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
        end
    end
end

@inline broadcast(f, a::StaticArray, n::Number) = map(x -> f(x, n), a)
@inline broadcast(f, n::Number, a::StaticArray) = map(x -> f(n, x), a)

# Other two-input versions with AbstractArray

##########
## map! ##
##########

# Single input
@generated function map!{F}(f::F, out::StaticArray, a1::StaticArray)
    exprs = [:(out[$j] = f(a1[$j])) for j = 1:length(a1)]
    return quote
        $(Expr(:meta, :inline))
        @inbounds $(Expr(:block, exprs...))
    end
end

# Two inputs
@generated function map!{F}(f::F, out::StaticArray, a1::StaticArray, a2::StaticArray)
    if size(a1) != size(a2)
        error("Dimensions must match. Got sizes $(size(a)) and $(size(a2))")
    end

    if size(a1) != size(a2)
        error("Dimensions must match. Got sizes $(size(out)) and $(size(a1))")
    end

    exprs = [:(out[$j] = f(a1[$j], a2[$j])) for j = 1:length(a1)]
    return quote
        #$(Expr(:meta, :inline))
        @inbounds $(Expr(:block, exprs...))
    end
end

################
## broadcast! ##
################

@inline broadcast!{F}(f::F, out::StaticArray, a::StaticArray) = map!(f, out, a)
@inline broadcast!(f::typeof(identity), out::StaticArray, a::StaticArray) = map!(f, out, a)

# Two input versions
@generated function broadcast!{F}(f::F, out::StaticArray, a1::StaticArray, a2::StaticArray)
    if size(a1) == size(a2) && size(out) == size(a1)
        return quote
            $(Expr(:meta, :inline))
            @inbounds map!(f, out, a1, a2)
        end
    else
        s1 = size(a1)
        s2 = size(a2)
        ndims = max(length(s1), length(s2))

        s = Vector{Int}(ndims)
        expands1 = Vector{Bool}(ndims)
        expands2 = Vector{Bool}(ndims)
        for i = 1:ndims
            if length(s1) < i
                s[i] = s2[i]
                expands1[i] = false
                expands2[i] = s2[i] > 1
            elseif length(s2) < i
                s[i] = s1[i]
                expands1[i] = s1[i] > 1
                expands2[i] = false
            else
                s[i] = max(s1[i], s1[i])
                @assert s1[i] == 1 || s1[i] == s[i]
                @assert s2[i] == 1 || s2[i] == s[i]
                expands1[i] = s1[i] > 1
                expands2[i] = s2[i] > 1
            end
        end
        s = (s...)
        L = prod(s)

        if s != size(out)
            error("Dimension mismatch")
        end

        exprs = Vector{Expr}(L)

        i = 1
        ind = ones(Int, ndims)
        while i <= L
            ind1 = [expands1[j] ? ind[j] : 1 for j = 1:length(s1)]
            ind2 = [expands2[j] ? ind[j] : 1 for j = 1:length(s2)]
            index1 = sub2ind(s1, ind1...)
            index2 = sub2ind(s2, ind2...)

            exprs[i] = :(out[$i] = $(Expr(:call, :f, Expr(:ref, :a1, index1), Expr(:ref, :a2, index2))))

            i += 1

            ind[1] += 1
            j = 1
            while j < length(s)
                if ind[j] > s[j]
                    ind[j] = 1
                    ind[j+1] += 1
                else
                    break
                end

                j += 1
            end
        end

        return quote
            $(Expr(:meta, :inline))
            @inbounds $(Expr(:block, exprs...))
        end
    end
end


text:
##############
# Owner Type #
##############

@ghdef mutable struct Owner
    typ::Union{String, Nothing}
    email::Union{String, Nothing}
    name::Union{String, Nothing}
    login::Union{String, Nothing}
    bio::Union{String, Nothing}
    company::Union{String, Nothing}
    location::Union{String, Nothing}
    avatar_url::Union{HTTP.URI, Nothing}
    gravatar_id::Union{String, Nothing}
    id::Union{Int, Nothing}
    public_repos::Union{Int, Nothing}
    owned_private_repos::Union{Int, Nothing}
    total_private_repos::Union{Int, Nothing}
    public_gists::Union{Int, Nothing}
    private_gists::Union{Int, Nothing}
    followers::Union{Int, Nothing}
    following::Union{Int, Nothing}
    collaborators::Union{Int, Nothing}
    blog::Union{HTTP.URI, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    date::Union{Dates.DateTime, Nothing}
    hireable::Union{Bool, Nothing}
    site_admin::Union{Bool, Nothing}
end

Owner(login::AbstractString, isorg = false) = Owner(Dict("login" => login, "type" => isorg ? "Organization" : "User"))

namefield(owner::Owner) = owner.login

typprefix(isorg) = isorg ? "orgs" : "users"

#############
# Owner API #
#############

isorg(owner::Owner) = something(owner.typ, "") == "Organization"

@api_default function whoami(api::GitHubAPI; options...)
    result = gh_get_json(api, "/user"; options...)
    return Owner(result)
end

@api_default owner(api::GitHubAPI, owner_obj::Owner; options...) = owner(api, name(owner_obj), isorg(owner_obj); options...)

@api_default function owner(api::GitHubAPI, owner_obj, isorg = false; options...)
    result = gh_get_json(api, "/$(typprefix(isorg))/$(name(owner_obj))"; options...)
    return Owner(result)
end

@api_default function users(api::GitHubAPI; options...)
    results, page_data = gh_get_paged_json(api, "/users"; options...)
    return map(Owner, results), page_data
end

@api_default function check_membership(api::GitHubAPI, org, user; public_only = false, options...)
    scope = public_only ? "public_members" : "members"
    resp = gh_get(api, "/orgs/$(name(org))/$scope/$(name(user))"; handle_error = false, allowredirects = false, options...)
    if resp.status == 204
        return true
    elseif resp.status == 404
        return false
    elseif resp.status == 302
        # For convenience, still check public membership. Otherwise, we don't know, so error
        @assert !public_only
        is_public_member = check_membership(org, user; public_only = true, options...)
        is_public_member && return true
        error("Enquiring about an Organization to which you do not have access.\n"*
              "Set `public_only=true` or provide authentication.")
    else
        handle_response_error(resp)
    end
end

@api_default function orgs(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/orgs"; options...)
    return map(Owner, results), page_data
end

@api_default function followers(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/followers"; options...)
    return map(Owner, results), page_data
end

@api_default function following(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/following"; options...)
    return map(Owner, results), page_data
end

@api_default function pubkeys(api::GitHubAPI, owner; options...)
    Base.depwarn("`pubkeys` is deprecated in favor of `sshkeys`, " *
        "which return a vector of keys, instead of a Dict from key-id to key.", :pubkeys)
    results, page_data = sshkeys(api, owner; options...)
    output = Dict{Int,String}([(key["id"], key["key"]) for key in results])
    return output, page_data
end

@api_default function sshkeys(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/keys"; options...)
    output = convert(Vector{Dict{String,Any}}, results)
    return output, page_data
end

@api_default function gpgkeys(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/users/$(name(owner))/gpg_keys"; options...)
    output = convert(Vector{Dict{String,Any}}, results)
    return output, page_data
end

repos(api::GitHubAPI, owner::Owner; options...) = repos(api, name(owner), isorg(owner); options...)

@api_default function repos(api::GitHubAPI, owner, isorg = false; options...)
    results, page_data = gh_get_paged_json(api, "/$(typprefix(isorg))/$(name(owner))/repos"; options...)
    return map(Repo, results), page_data
end

@api_default function teams(api::GitHubAPI, owner; options...)
    results, page_data = gh_get_paged_json(api, "/orgs/$(name(owner))/teams"; options...)
    return map(Team, results), page_data
end


text:
include("bimodal_cauchy_distribution.jl")

const ADE_DefaultOptions = chain(DE_DefaultOptions, ParamsDict(
    # Distributions we will use to generate new F and CR values.
    :fdistr => BimodalCauchy(0.65, 0.1, 1.0, 0.1, clampBelow0 = false),
    :crdistr => BimodalCauchy(0.1, 0.1, 0.95, 0.1, clampBelow0 = false),
    :SearchSpace => RectSearchSpace(1)
))

"""
Specific data and functions for adaptation
An Adaptive DE typically changes parameters of the search dynamically. This is
typically done in the `tell!()` function when we know if the trial vector
was better than the target vector.
"""
mutable struct AdaptiveDiffEvoParameters
    # Distributions we will use to generate new F and CR values.
    # In the literature Cauchy distributions have been used for sampling the
    # `f` and `cr` constants
    # FIXME allow any distribution?
    fdistr::BimodalCauchy
    crdistr::BimodalCauchy

    fs::Vector{Float64}   # One f value per individual in population
    crs::Vector{Float64}  # One cr value per individual in population

    function AdaptiveDiffEvoParameters(fdistr::BimodalCauchy, crdistr::BimodalCauchy)
        @assert !fdistr.clampBelow0 && !crdistr.clampBelow0 # population probs should not be degenerated
        new(fdistr, crdistr,
        Vector{Float64}(), Vector{Float64}()) # start with empty arrays because the population size unknown
    end
end

AdaptiveDiffEvoParameters(options::Parameters) =
    AdaptiveDiffEvoParameters(options[:fdistr], options[:crdistr])

function crossover_parameters(params::AdaptiveDiffEvoParameters, pop::Population, target_index)
    # initialize the f & cr array
    if isempty(params.fs)
        params.fs = rand!(params.fdistr, Vector{Float64}(undef, popsize(pop)))
    end
    if isempty(params.crs)
        params.crs = rand!(params.crdistr, Vector{Float64}(undef, popsize(pop)))
    end
    return (params.crs[target_index], params.fs[target_index])
end

function adjust!(params::AdaptiveDiffEvoParameters, index, is_improved::Bool)
    if !is_improved
        # The trial vector for this target was not better so we change the f and cr constants.
        @inbounds params.fs[index] = rand(params.fdistr)
        @inbounds params.crs[index] = rand(params.crdistr)
    end
end

"""
An Adaptive DE crossover operator changes `cr` and `f` parameters of the search dynamically.
"""
mutable struct AdaptiveDiffEvoRandBin{N} <: DiffEvoCrossoverOperator{N,1}
    params::AdaptiveDiffEvoParameters

    AdaptiveDiffEvoRandBin{N}(params::AdaptiveDiffEvoParameters) where N = new{N}(params)

    AdaptiveDiffEvoRandBin{N}(options::Parameters) where N =
        new{N}(AdaptiveDiffEvoParameters(options))
end

crossover_parameters(xover::AdaptiveDiffEvoRandBin, pop::Population, target_index) =
    crossover_parameters(xover.params, pop, target_index)

adjust!(xover::AdaptiveDiffEvoRandBin, op_index::Int, candi_index::Int,
        new_fitness::F, old_fitness::F, is_improved::Bool) where F =
    adjust!(xover.params, candi_index, is_improved)

const AdaptiveDiffEvoRandBin1 = AdaptiveDiffEvoRandBin{3}
const AdaptiveDiffEvoRandBin2 = AdaptiveDiffEvoRandBin{5}

function adaptive_diffevo(problem::OptimizationProblem,
                options::Parameters, name::String,
                select::IndividualsSelector = SimpleSelector(),
                crossover::DiffEvoCrossoverOperator =
                AdaptiveDiffEvoRandBin1(chain(ADE_DefaultOptions, options)))
  opts = chain(ADE_DefaultOptions, options)
  pop = population(problem, opts)
  DiffEvoOpt(name, pop, select, crossover,
             RandomBound(search_space(problem)))
end

adaptive_de_rand_1_bin(problem::OptimizationProblem,
                       options::Parameters = EMPTY_PARAMS,
                       name = "AdaptiveDE/rand/1/bin") =
    adaptive_diffevo(problem, options, name)

adaptive_de_rand_1_bin_radiuslimited(problem::OptimizationProblem,
                                     options::Parameters = EMPTY_PARAMS,
                                     name = "AdaptiveDE/rand/1/bin/radiuslimited") =
    adaptive_diffevo(problem, options, name,
                     RadiusLimitedSelector(chain(ADE_DefaultOptions, options)[:SamplerRadius]))


text:
"""
    append_gs!(model::JuMP.Model, x, ng::Int)

Append inequality constraints (auto-generated)

# Arguments
    `model::JuMP.Model`: optimization model
    `x::JuMP.variables`: optimization model variables
    `ng::Int`: number of inequality constraints
"""
function append_gs!(model::JuMP.Model, x, ng::Int)
    for ig = 1:ng
        if ig == 1
            @NLconstraint(model, -Inf <= g1(x...) <= 0)
        elseif ig == 2
            @NLconstraint(model, -Inf <= g2(x...) <= 0)
        elseif ig == 3
            @NLconstraint(model, -Inf <= g3(x...) <= 0)
        elseif ig == 4
            @NLconstraint(model, -Inf <= g4(x...) <= 0)
        elseif ig == 5
            @NLconstraint(model, -Inf <= g5(x...) <= 0)
        elseif ig == 6
            @NLconstraint(model, -Inf <= g6(x...) <= 0)
        elseif ig == 7
            @NLconstraint(model, -Inf <= g7(x...) <= 0)
        elseif ig == 8
            @NLconstraint(model, -Inf <= g8(x...) <= 0)
        elseif ig == 9
            @NLconstraint(model, -Inf <= g9(x...) <= 0)
        elseif ig == 10
            @NLconstraint(model, -Inf <= g10(x...) <= 0)
        elseif ig == 11
            @NLconstraint(model, -Inf <= g11(x...) <= 0)
        elseif ig == 12
            @NLconstraint(model, -Inf <= g12(x...) <= 0)
        elseif ig == 13
            @NLconstraint(model, -Inf <= g13(x...) <= 0)
        elseif ig == 14
            @NLconstraint(model, -Inf <= g14(x...) <= 0)
        elseif ig == 15
            @NLconstraint(model, -Inf <= g15(x...) <= 0)
        elseif ig == 16
            @NLconstraint(model, -Inf <= g16(x...) <= 0)
        elseif ig == 17
            @NLconstraint(model, -Inf <= g17(x...) <= 0)
        elseif ig == 18
            @NLconstraint(model, -Inf <= g18(x...) <= 0)
        elseif ig == 19
            @NLconstraint(model, -Inf <= g19(x...) <= 0)
        elseif ig == 20
            @NLconstraint(model, -Inf <= g20(x...) <= 0)
        end
    end
end

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
# Automatically generated using Clang.jl

const P4EST_MAXLEVEL = 30
const P4EST_QMAXLEVEL = 29

# Skipping MacroDefinition: P4EST_ROOT_LEN ( ( p4est_qcoord_t ) 1 << P4EST_MAXLEVEL )
# Skipping MacroDefinition: P4EST_QUADRANT_LEN ( l ) ( ( p4est_qcoord_t ) 1 << ( P4EST_MAXLEVEL - ( l ) ) )
# Skipping MacroDefinition: P4EST_LAST_OFFSET ( l ) ( P4EST_ROOT_LEN - P4EST_QUADRANT_LEN ( l ) )
# Skipping MacroDefinition: P4EST_QUADRANT_INIT ( q ) ( ( void ) memset ( ( q ) , - 1 , sizeof ( p4est_quadrant_t ) ) )

struct p4est_quadrant
    x::p4est_qcoord_t
    y::p4est_qcoord_t
    level::Int8
    pad8::Int8
    pad16::Int16
    p::p4est_quadrant_data
end

const p4est_quadrant_t = p4est_quadrant

struct p4est_tree
    quadrants::sc_array_t
    first_desc::p4est_quadrant_t
    last_desc::p4est_quadrant_t
    quadrants_offset::p4est_locidx_t
    quadrants_per_level::NTuple{31, p4est_locidx_t}
    maxlevel::Int8
end

const p4est_tree_t = p4est_tree

struct p4est_inspect
    use_balance_ranges::Cint
    use_balance_ranges_notify::Cint
    use_balance_verify::Cint
    balance_max_ranges::Cint
    balance_A_count_in::Csize_t
    balance_A_count_out::Csize_t
    balance_comm_sent::Csize_t
    balance_comm_nzpeers::Csize_t
    balance_B_count_in::Csize_t
    balance_B_count_out::Csize_t
    balance_zero_sends::Cint
    balance_zero_receives::Cint
    balance_A::Cdouble
    balance_comm::Cdouble
    balance_B::Cdouble
    balance_ranges::Cdouble
    balance_notify::Cdouble
    balance_notify_allgather::Cdouble
    use_B::Cint
end

const p4est_inspect_t = p4est_inspect

struct p4est_connectivity
    num_vertices::p4est_topidx_t
    num_trees::p4est_topidx_t
    num_corners::p4est_topidx_t
    vertices::Ptr{Cdouble}
    tree_to_vertex::Ptr{p4est_topidx_t}
    tree_attr_bytes::Csize_t
    tree_to_attr::Cstring
    tree_to_tree::Ptr{p4est_topidx_t}
    tree_to_face::Ptr{Int8}
    tree_to_corner::Ptr{p4est_topidx_t}
    ctt_offset::Ptr{p4est_topidx_t}
    corner_to_tree::Ptr{p4est_topidx_t}
    corner_to_corner::Ptr{Int8}
end

const p4est_connectivity_t = p4est_connectivity

struct p4est
    mpicomm::MPI_Comm
    mpisize::Cint
    mpirank::Cint
    mpicomm_owned::Cint
    data_size::Csize_t
    user_pointer::Ptr{Cvoid}
    revision::Clong
    first_local_tree::p4est_topidx_t
    last_local_tree::p4est_topidx_t
    local_num_quadrants::p4est_locidx_t
    global_num_quadrants::p4est_gloidx_t
    global_first_quadrant::Ptr{p4est_gloidx_t}
    global_first_position::Ptr{p4est_quadrant_t}
    connectivity::Ptr{p4est_connectivity_t}
    trees::Ptr{sc_array_t}
    user_data_pool::Ptr{sc_mempool_t}
    quadrant_pool::Ptr{sc_mempool_t}
    inspect::Ptr{p4est_inspect_t}
end

const p4est_t = p4est
const p4est_init_t = Ptr{Cvoid}
const p4est_refine_t = Ptr{Cvoid}
const p4est_coarsen_t = Ptr{Cvoid}
const p4est_weight_t = Ptr{Cvoid}
const P4EST_DIM = 2
const P4EST_FACES = 2P4EST_DIM
const P4EST_CHILDREN = 4
const P4EST_HALF = P4EST_CHILDREN / 2
const P4EST_INSUL = 9
const P4EST_FTRANSFORM = 9
const P4EST_STRING = "p4est"
const P4EST_ONDISK_FORMAT = 0x02000009

@cenum p4est_connect_type_t::UInt32 begin
    P4EST_CONNECT_FACE = 21
    P4EST_CONNECT_CORNER = 22
    P4EST_CONNECT_FULL = 22
end

@cenum p4est_connectivity_encode_t::UInt32 begin
    P4EST_CONN_ENCODE_NONE = 0
    P4EST_CONN_ENCODE_LAST = 1
end

struct p4est_corner_transform_t
    ntree::p4est_topidx_t
    ncorner::Int8
end

struct p4est_corner_info_t
    icorner::p4est_topidx_t
    corner_transforms::sc_array_t
end

const p4est_replace_t = Ptr{Cvoid}
const p4est_geometry_X_t = Ptr{Cvoid}
const p4est_geometry_destroy_t = Ptr{Cvoid}

struct p4est_geometry
    name::Cstring
    user::Ptr{Cvoid}
    X::p4est_geometry_X_t
    destroy::p4est_geometry_destroy_t
end

const p4est_geometry_t = p4est_geometry

struct p4est_ghost_t
    mpisize::Cint
    num_trees::p4est_topidx_t
    btype::p4est_connect_type_t
    ghosts::sc_array_t
    tree_offsets::Ptr{p4est_locidx_t}
    proc_offsets::Ptr{p4est_locidx_t}
    mirrors::sc_array_t
    mirror_tree_offsets::Ptr{p4est_locidx_t}
    mirror_proc_mirrors::Ptr{p4est_locidx_t}
    mirror_proc_offsets::Ptr{p4est_locidx_t}
    mirror_proc_fronts::Ptr{p4est_locidx_t}
    mirror_proc_front_offsets::Ptr{p4est_locidx_t}
end

struct p4est_ghost_exchange
    is_custom::Cint
    is_levels::Cint
    p4est::Ptr{p4est_t}
    ghost::Ptr{p4est_ghost_t}
    minlevel::Cint
    maxlevel::Cint
    data_size::Csize_t
    ghost_data::Ptr{Cvoid}
    qactive::Ptr{Cint}
    qbuffer::Ptr{Cint}
    requests::sc_array_t
    sbuffers::sc_array_t
    rrequests::sc_array_t
    rbuffers::sc_array_t
end

const p4est_ghost_exchange_t = p4est_ghost_exchange

struct p4est_iter_volume_info
    p4est::Ptr{p4est_t}
    ghost_layer::Ptr{p4est_ghost_t}
    quad::Ptr{p4est_quadrant_t}
    quadid::p4est_locidx_t
    treeid::p4est_topidx_t
end

const p4est_iter_volume_info_t = p4est_iter_volume_info
const p4est_iter_volume_t = Ptr{Cvoid}

struct p4est_iter_face_side
    treeid::p4est_topidx_t
    face::Int8
    is_hanging::Int8
    is::p4est_iter_face_side_data
end

const p4est_iter_face_side_t = p4est_iter_face_side

struct p4est_iter_face_info
    p4est::Ptr{p4est_t}
    ghost_layer::Ptr{p4est_ghost_t}
    orientation::Int8
    tree_boundary::Int8
    sides::sc_array_t
end

const p4est_iter_face_info_t = p4est_iter_face_info
const p4est_iter_face_t = Ptr{Cvoid}

struct p4est_iter_corner_side
    treeid::p4est_topidx_t
    corner::Int8
    is_ghost::Int8
    quad::Ptr{p4est_quadrant_t}
    quadid::p4est_locidx_t
    faces::NTuple{2, Int8}
end

const p4est_iter_corner_side_t = p4est_iter_corner_side

struct p4est_iter_corner_info
    p4est::Ptr{p4est_t}
    ghost_layer::Ptr{p4est_ghost_t}
    tree_boundary::Int8
    sides::sc_array_t
end

const p4est_iter_corner_info_t = p4est_iter_corner_info
const p4est_iter_corner_t = Ptr{Cvoid}
const p4est_lnodes_code_t = Int8

struct p4est_lnodes
    mpicomm::MPI_Comm
    num_local_nodes::p4est_locidx_t
    owned_count::p4est_locidx_t
    global_offset::p4est_gloidx_t
    nonlocal_nodes::Ptr{p4est_gloidx_t}
    sharers::Ptr{sc_array_t}
    global_owned_count::Ptr{p4est_locidx_t}
    degree::Cint
    vnodes::Cint
    num_local_elements::p4est_locidx_t
    face_code::Ptr{p4est_lnodes_code_t}
    element_nodes::Ptr{p4est_locidx_t}
end

const p4est_lnodes_t = p4est_lnodes

struct p4est_lnodes_rank
    rank::Cint
    shared_nodes::sc_array_t
    shared_mine_offset::p4est_locidx_t
    shared_mine_count::p4est_locidx_t
    owned_offset::p4est_locidx_t
    owned_count::p4est_locidx_t
end

const p4est_lnodes_rank_t = p4est_lnodes_rank

struct p4est_lnodes_buffer
    requests::Ptr{sc_array_t}
    send_buffers::Ptr{sc_array_t}
    recv_buffers::Ptr{sc_array_t}
end

const p4est_lnodes_buffer_t = p4est_lnodes_buffer

struct p4est_mesh_t
    local_num_quadrants::p4est_locidx_t
    ghost_num_quadrants::p4est_locidx_t
    quad_to_tree::Ptr{p4est_topidx_t}
    ghost_to_proc::Ptr{Cint}
    quad_to_quad::Ptr{p4est_locidx_t}
    quad_to_face::Ptr{Int8}
    quad_to_half::Ptr{sc_array_t}
    quad_level::Ptr{sc_array_t}
    local_num_corners::p4est_locidx_t
    quad_to_corner::Ptr{p4est_locidx_t}
    corner_offset::Ptr{sc_array_t}
    corner_quad::Ptr{sc_array_t}
    corner_corner::Ptr{sc_array_t}
end

struct p4est_mesh_face_neighbor_t
    p4est::Ptr{p4est_t}
    ghost::Ptr{p4est_ghost_t}
    mesh::Ptr{p4est_mesh_t}
    which_tree::p4est_topidx_t
    quadrant_id::p4est_locidx_t
    quadrant_code::p4est_locidx_t
    face::Cint
    subface::Cint
    current_qtq::p4est_locidx_t
end

const p4est_search_local_t = Ptr{Cvoid}
const p4est_search_query_t = p4est_search_local_t
const p4est_search_partition_t = Ptr{Cvoid}
const p4est_search_all_t = Ptr{Cvoid}
const p4est_vtk_context = Cvoid
const p4est_vtk_context_t = p4est_vtk_context


text:
module CUDA

using Libdl
using Base.Threads
import LinearAlgebra.BLAS: scal!, axpy!, gemv, gemv!, gemm, gemm!

if Sys.iswindows()
    const libcuda = Libdl.find_library("nvcuda")
else
    const libcuda = Libdl.find_library("libcuda")
end
const AVAILABLE = !isempty(libcuda)
AVAILABLE || @warn "CUDA cannot be found."
const API_VERSION = Ref{Int}()

function checkstatus(status)
    if status != 0
        ref = Ref{Cstring}()
        ccall((:cuGetErrorString,libcuda), Cint, (Cint,Ptr{Cstring}), status, ref)
        throw(unsafe_string(ref[]))
    end
end

if AVAILABLE
    status = ccall((:cuInit,libcuda), Cint, (Cint,), 0)
    checkstatus(status)

    ref = Ref{Cint}()
    status = ccall((:cuDriverGetVersion,libcuda), Cint, (Ptr{Cint},), ref)
    checkstatus(status)

    API_VERSION[] = Int(ref[])
    @info "CUDA API $(API_VERSION[])"
end

include("define.jl")

macro apicall(f, args...)
    f = get(DEFINE, f.value, f.value)
    quote
        status = ccall(($(QuoteNode(f)),libcuda), Cint, $(map(esc,args)...))
        checkstatus(status)
    end
end

macro unsafe_apicall(f, args...)
    f = get(DEFINE, f.value, f.value)
    quote
        ccall(($(QuoteNode(f)),libcuda), Cint, $(map(esc,args)...))
    end
end

const Cptr = Ptr{Cvoid}
export cstring
cstring(::Type{Int32}) = "int"
cstring(::Type{Float32}) = "float"

include("driver/device.jl")
include("driver/context.jl")
include("driver/stream.jl")
include("driver/memory.jl")
include("driver/module.jl")
include("driver/function.jl")

if AVAILABLE
    const CONTEXTS = Array{CuContext}(undef, ndevices())
    # This must be loaded before kernel.jl and kernels.jl
    include("nvml/NVML.jl")
    include("nvrtc/NVRTC.jl")
    using .NVML
end

include("pointer.jl")
include("array.jl")
include("subarray.jl")
include("kernel.jl")
include("arraymath.jl")
include("broadcast.jl")
include("cat.jl")
include("reduce.jl")
include("devicearray.jl")

include("allocators/atomic_malloc.jl")
include("allocators/cuda_malloc.jl")
include("allocators/mempool_malloc.jl")

const ALLOCATOR = Ref{Any}(CUDAMalloc())

if AVAILABLE
    include("nccl/NCCL.jl")
    include("cublas/CUBLAS.jl")
    include("cudnn/CUDNN.jl")

    using .CUBLAS, .CUDNN
    export CUBLAS, CUDNN
end

end


text:
using LatBo.LB: velocity, momentum, density, equilibrium, D2Q9, D3Q19, LocalQuantities,
        LocalQuantities!

facts("Thermodynamic quantities and functions") do
    for lattice_name in [:D2Q9, :D3Q19]
        context("for $(lattice_name)") do
            lattice = eval(lattice_name)
            ncomponents = size(lattice.celerities, 2)

            context("rho is a geometric series") do
                series(N) = N * (N+1)/2
                @fact density(collect(1:10)) --> series(10)
                @fact density(collect(10:15)) --> series(15) - series(9)
            end

            context("homogeneous populations sum to zero momentum") do
                fᵢ = ones(Float64, ncomponents)
                @fact momentum(fᵢ, lattice.celerities) .== 0 --> all
                fᵢ[1] = 3 # This is the zero momentum component
                @fact momentum(fᵢ, lattice.celerities) .== 0 --> all
                fᵢ[2] = 2 # Now for a negative test
                @fact momentum(fᵢ, lattice.celerities) .!= 0 --> any
            end

            context("velocity from momentum and density ") do
                fᵢ = convert(Array{Float64}, 1.0 + rand(ncomponents))
                actual = velocity(fᵢ, lattice.celerities, 1.)
                half = velocity(fᵢ, lattice.celerities, 0.5)
                μ = momentum(fᵢ, lattice.celerities)
                @fact actual --> roughly(half * 0.5)
                @fact actual --> roughly(μ)
            end

            context("equilibrium function") do
                wᵢ = Float64[1, 2, 3]
                ē = Int64[1 0 1; 0 1 1]
                μ = Float64[1, 2]
                ρ = 1.1

                @fact equilibrium(ρ, zeros(2), ē, wᵢ) --> roughly(wᵢ * ρ)
                @fact equilibrium(ρ, μ, zeros(Int64, 2, 3), wᵢ) --> roughly(wᵢ * (ρ - 1.5  * 5./ρ))

                fᵉ = wᵢ .* (ρ - 1.5dot(μ, μ)/ρ + 4.5(ē.'μ).^2/ρ + 3ē.'μ)
                @fact equilibrium(ρ, μ, ē, wᵢ) --> roughly(fᵉ)
            end

            context("local quantities aggregator") do
                fᵢ = convert(Array{Float64}, 1.0 + rand(Float64, ncomponents))
                μ = momentum(fᵢ, lattice.celerities)
                ν = velocity(fᵢ, lattice.celerities)
                ρ = density(fᵢ)
                feq = equilibrium(lattice, ρ, μ)

                quants = LocalQuantities(fᵢ, lattice)
                @fact quants.density --> roughly(ρ)
                @fact quants.momentum --> roughly(μ)
                @fact quants.velocity --> roughly(ν)
                @fact quants.feq --> roughly(feq)
            end

            context("pre-allocated local aggregator") do
                fᵢ = convert(Array{Float64}, 1.0 + rand(Float64, ncomponents))
                quants = LocalQuantities(fᵢ, lattice)
                pre = LocalQuantities(lattice)
                @fact size(pre.momentum) --> size(quants.momentum)
                @fact size(pre.velocity) --> size(quants.velocity)
                @fact size(pre.feq) --> size(quants.feq)

                LocalQuantities!(pre, fᵢ, lattice)
                @fact pre.density --> roughly(quants.density)
                @fact pre.momentum --> roughly(quants.momentum)
                @fact pre.velocity --> roughly(quants.velocity)
                @fact pre.feq --> roughly(quants.feq)
            end
        end
    end
end


text:
"""
Imbibition cell, IMPES
Buckley Leverett equation
dependent variables: pressure and water saturation
Prepared for educational purposes by ** AAE **
The only input variable is 0<=SF0<=1,
is the degree of wettability
"""
function forced_imb_impes(SF0, SW0)
Nx = 100 # number of cells in x direction
Ny = 40 # number of cells in y direction
W = 0.02 # [m] length of the domain in x direction
H = 0.07 # [m] length of the domain in y direction
A_core=π*W^2/4
V_core=H*A_core
m = createMesh1D(Nx, H)
# m = createMeshCylindrical2D(Nx, Ny, W, H) # creates a 2D mesh
## define the physical parametrs
# all the Corey-type relperm parametrs are defined for the oil-wet and water-wet
# cases
k0 = 0.001e-12 # [m^2] average reservoir permeability
phi0 = 0.35 # average porosity
mu_oil = 2e-3 # [Pa.s] oil viscosity
mu_water = 1e-3 # [Pa.s] water viscosity
krw0_ww = 0.3
krw0_ow = 1.0
kro0_ww = 0.6
kro0_ow = 0.76
nw_ww = 2.4
nw_ow= 2.4
no_ww = 2.0
no_ow= 2.0
sor_ww=0.15
sor_ow=0.2
sor_mw=0.08
swc_ww=0.09
swc_ow=0.09
teta_ow=deg2rad(180)
teta_ww=deg2rad(0)
teta_mw=deg2rad(90)
teta_prime=teta_mw-π/2 # see my report
a=(sor_mw-0.5*(sor_ww-sor_ow)*cos(teta_mw)-0.5*(sor_ww+sor_ow))/(cos(teta_mw)^2-1)
b=0.5*(sor_ww-sor_ow)
c=0.5*(sor_ww+sor_ow)-a

gama_ow=0.03 # N/m
labda=2.4 # for Ekofisk chalk
r_ave=sqrt(k0/phi0) #  meter average pore diameter
pce0= 5e2 #2*gama_ow*cos(teta_ww)/r_ave/10 # Pa capillary entry pressure
# SF0=1.0
SF=createFaceVariable(m, SF0) # 1 is water wet, 0 is oil wet
krw0=krw0_ww*SF+krw0_ow*(1-SF)
kro0=kro0_ww*SF+kro0_ow*(1-SF)
# sor=sor_ww*SF+sor_ow*(1-SF)
swc=swc_ww*SF+swc_ow*(1-SF)
no= no_ww*SF+no_ow*(1-SF)
nw= nw_ww*SF+nw_ow*(1-SF)
teta0=teta_ww*SF0+teta_ow*(1-SF0)
# sor0=(sor_ow-sor_mw)*cos(teta0-teta_prime)^2+sor_mw
sor0=a*cos(teta0)^2+b*cos(teta0)+c
teta=createFaceVariable(m, teta0)
sor=createFaceVariable(m, sor0)
pce=createFaceVariable(m, pce0)

p0 = 100e5 # [bar] pressure
pin = 150e5 # [bar] injection pressure at the left boundary
u_in= 1.0/(24*3600) # [m/s] equal to 1 m/day
sw0 = swc_ww+0.01
# sw0(10:end-10, 10:end-10)=swc+0.2
# sw0 = swc+0.1 # initial water saturation
sw_in = 1

eps1=1e-6
clx=0.05
cly=0.05
V_dp=0.01 # Dykstra-Parsons coef.
perm_val= k0 #field2d(Nx,Ny,k0,V_dp,clx,cly)

k=createCellVariable(m, perm_val)
phi=createCellVariable(m, phi0)

lw = geometricMean(k)/mu_water
lo = geometricMean(k)/mu_oil

## Define the boundaries: all fixed Sw=1, fixed pressure everywhere(?)
BCp = createBC(m) # Neumann BC for pressure
BCs = createBC(m) # Neumann BC for saturation
# left boundary pressure gradient
# BCp.left.a[:]=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)) BCp.left.b[:]=0 BCp.left.c[:]=-u_in
# change the right boandary to constant pressure (Dirichlet)
BCp.left.a[:]=-k0/mu_water
BCp.left.b[:]=0
BCp.left.c[:]=u_in
BCp.right.a[:]=0.0
BCp.right.b[:]=1.0
BCp.right.c[:]=p0
# BCp.top.a[:]=0.0
# BCp.top.b[:]=1.0
# BCp.top.c[:]=p0
# BCp.bottom.a[:]=0.0
# BCp.bottom.b[:]=1.0
# BCp.bottom.c[:]=p0
# change the left boundary to constant saturation (Dirichlet)
# BCs.left.a[:]=0 BCs.left.b[:]=1 BCs.left.c[:]=1.0-sor
# calculate the water saturation at the boundaries
f1(sw)=pc_imb(sw, pce0, swc_ww*SF0+swc_ow*(1-SF0), sor_ww*SF0+sor_ow*(1-SF0),
  teta=teta_ww*SF0+teta_ow*(1-SF0))
sw_b=fzero(f1, 0.5)
println(sw_b)
BCs.right.a[:]=0.0
BCs.right.b[:]=1.0
BCs.right.c[:]=sw_b
BCs.left.a[:]=0.0
BCs.left.b[:]=1.0
BCs.left.c[:]=sw_in

# BCs.top.a[:]=0.0
# BCs.top.b[:]=1.0
# BCs.top.c[:]=sw_b
# BCs.bottom.a[:]=0.0
# BCs.bottom.b[:]=1.0
# BCs.bottom.c[:]=sw_b
## define the time step and solver properties
# dt = 1000 # [s] time step
# dt=(W/Nx)/u_in/20 # [s]
pv_inj=5
t_end = pv_inj*V_core/(u_in*A_core) # [s] final time
dt=H/Nx/u_in
eps_p = 1e-7 # pressure accuracy
eps_sw = 1e-7 # saturation accuracy
## define the variables
sw_old = createCellVariable(m, SW0, BCs)
p_old = createCellVariable(m, p0, BCp)
sw = copyCell(sw_old)
oil_init=domainInt(1-sw_old)
p = copyCell(p_old)
uw = -gradientTerm(p_old) # an estimation of the water velocity
## start the main loop
# generate intial pressure profile (necessary to initialize the fully
# implicit solver)
rec_fact=zeros(1)
t_day=zeros(1)
t = 0.0
dt0=dt
dsw_alwd= 0.01
dp_alwd= 100.0 # Pa
labdaw=createFaceVariable(m,1.0) # initialize
# prog_1=ProgressThresh(t_end)
while (t<t_end)
  print("progress is $((t/t_end*100)) %", "\u1b[1G")
  # for i=1:5
    error_p = 1e5
    error_sw = 1e5
    # Implicit loop
#     while ((error_p>eps_p) || (error_sw>eps_sw))
    while(true)
        # calculate parameters
        pgrad = gradientTerm(p)
#         pcgrad=gradientTerm(pc(sw))
        sw_face = upwindMean(sw, -labdaw.*pgrad) # average value of water saturation
        sw_grad=gradientTerm(sw)
        sw_ave=arithmeticMean(sw)
        pcgrad=faceEval(dpc_imb, sw_ave, pce, swc, sor, teta).*sw_grad
        # solve for pressure at known Sw
        labdao = lo.*faceEval(kro, sw_face, kro0, sor, swc, no)
        labdaw = lw.*faceEval(krw, sw_face, krw0, sor, swc, nw)
        labda = labdao+labdaw
        # compute [Jacobian] matrices
        Mdiffp1 = diffusionTerm(-labda)
        RHSpc1=divergenceTerm(labdao.*pcgrad)
        Mbcp, RHSbcp = boundaryConditionTerm(BCp)
        RHS1 = RHSpc1+RHSbcp # with capillary
        p_new=solvePDE(m, Mdiffp1+Mbcp, RHS1)

        # solve for Sw
        pgrad = gradientTerm(p_new)
        uw=-labdaw.*pgrad
        Mbcsw, RHSbcsw = boundaryConditionTerm(BCs)
        RHS_sw=-divergenceTerm(uw)
        sw_new=solveExplicitPDE(sw_old, dt, RHS_sw, BCs, phi)

        error_p = maximum(abs((internalCells(p_new)-internalCells(p))./internalCells(p_new)))
        error_sw = maximum(abs(internalCells(sw_new)-internalCells(sw)))
        dt_new=dt*min(dp_alwd/error_p, dsw_alwd/error_sw)
        # print("sw_error = $error_sw \n")
        # print("new time step = $dt \n")
        # print(internalCells(sw_new))
        # print(internalCells(sw_old))
        # sleep(1.0)
        # assign new values of p and sw
        if error_sw>dsw_alwd
            dt=dt*(dsw_alwd/error_sw)
            # print("new time step = $dt \n")
        else
            t=t+dt
            p = copyCell(p_new)
            sw = copyCell(sw_new)
            p_old = copyCell(p)
            sw_old = copyCell(sw)
            dt=min(dt*(dsw_alwd/error_sw), 15*dt, t_end-t-eps())
            break
        end
    end

    rec_fact=push!(rec_fact, (oil_init-domainInt(1-sw))/oil_init)
    t_day=push!(t_day, t)
    # print(t)
    # GR.plot(sw.value[2:end-1])
    # GR.imshow(sw.value[2:end-1,2:end-1])
    # visualizeCells(1-sw)
    # GR.plot(t_day/3600/24, rec_fact)
    # xlabel('time [day]')
    # ylabel('recovery factor')
    # title([num2str(t/3600/24) ' day']) drawnow
end
t_day, rec_fact, sw
end # imb_impes

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
x = 2
_test_completion_232_


text:
# conversion methods for CategoricalArrays

function rcopy(::Type{CategoricalArray}, s::Ptr{IntSxp})
    isFactor(s) || error("s is not an R factor")
    refs = UInt32[isNA(x) ? zero(UInt32) : UInt32(x) for x in s]
    levels = rcopy(Array{String}, getattrib(s,Const.LevelsSymbol))
    pool = CategoricalPool(levels, isOrdered(s))
    if anyna(s)
        CategoricalArray{Union{String, Missing}, 1}(refs, pool)
    else
        CategoricalArray{String, 1}(refs, pool)
    end
end

## CategoricalArray to sexp conversion.

function sexp(::Type{RClass{:factor}}, v::CategoricalArray)
    rv = protect(sexp(RClass{:integer}, v.refs))
    try
        for (i,ref) = enumerate(v.refs)
            if ref == 0
                rv[i] = naeltype(IntSxp)
            end
        end
        # due to a bug of CategoricalArrays, we use index(v.pool) instead of index(v)
        setattrib!(rv, Const.LevelsSymbol, CategoricalArrays.index(v.pool))
        setattrib!(rv, Const.ClassSymbol, "factor")
        if CategoricalArrays.isordered(v)
            rv = rcall(:ordered, rv, CategoricalArrays.levels(v))
        end
    finally
        unprotect(1)
    end
    rv
end