text:
"""
    Solid(; shape = nothing,
            solidMaterial = nothing,
            massProperties = nothing,
            collision = false,
            contactMaterial = "",
            collisionSmoothingRadius = 0.0,
            visualMaterial = VisualMaterial())

Generate a [Solid](@ref) with physical behavior of a rigid body with mass, visualization and collision properties.
`Solid` is used as `feature` of an [Object3D](@ref).

# Arguments
- `shape`: Defines the shape of the solid. Only [Shapes](@ref) from [---- For Solid and Visual ----](@ref) can be used.
    If undefined the solid is not considered in visualization or [Collision Handling](@ref) and `massProperties` must
    be explicitly defined.

- `solidMaterial`: Defines the material of the solid used for mass properties computation and (if `contactMaterial` is
    undefined) contact force calculation. A pre-defined [Solid material](@ref) from palettes/solidMaterials.json (e.g.
    `"Steel"`) or a user-defined [Solid material](@ref) (e.g. `SolidMaterial(density=7850)`) can be used.

- `massProperties`: Defines the mass properties of the solid. If `massProperties` is undefined, the mass properties are
    computed from `shape` and `solidMaterial`. It is also possible to define only the mass of the solid and compute the
    center of mass and inertia tensor from `shape` (e.g. `MassPropertiesFromShapeAndMass(mass=4.5)`) or explicitly
    define all mass properties (e.g. `MassProperties(mass=4.5, centerOfMass=[1.0,2.0,3.0], Ixx=4.0, Iyy=5.0, Izz=6.0,
    Ixy=4.5, Ixz=4.6, Iyz=5.5)`).

- `collision`: Defines if the solid is considered in [Collision Handling](@ref). `collision=true` requires definition
    of `shape`, mass properties (defined by `solidMaterial` or `massProperties`) and contact material (defined by
    `solidMaterial` or `contactMaterial`).

- `contactMaterial`: Defines the contact material for force computation with [Contact Force Law](@ref) (e.g. `"Steel"`).
    If undefined `solidMaterial` is used as contact material. Only contact material combinations defined in
    palettes/contactPairMaterials.json can be used.

- `collisionSmoothingRadius`: Defines a collision smoothing radius for surface edges.

- `visualMaterial`: Defines the material of the solid used for visualization. A pre-defined [Visual material](@ref)
    from palettes/visualMaterials.json (e.g. `"RedTransparent"`) or a user-defined [Visual material](@ref) (e.g.
    `VisualMaterial(color="DeepSkyBlue4", transparency=0.75)`) can be used.

# Example
```julia
import Modia3D

Solid(; shape = Sphere(diameter=0.5),
        solidMaterial = "DryWood",
        collision = true,
        contactMaterial = "Steel",
        collisionSmoothingRadius = 0.001,
        visualMaterial = VisualMaterial(color="DarkViolet"))
```
"""
struct Solid <: Modia3D.AbstractObject3DFeature
    shape::Union{Modia3D.AbstractGeometry,Nothing}
    solidMaterial::Union{SolidMaterial,Nothing}
    massProperties::Union{MassProperties,Nothing}
    collision::Bool
    contactMaterial::Union{String,Modia3D.AbstractContactMaterial,Nothing}
    collisionSmoothingRadius::Union{Float64,Nothing}
    visualMaterial::Union{Shapes.VisualMaterial,Nothing}

    function Solid(;
        shape::Union{Modia3D.AbstractGeometry,Nothing} = nothing,
        solidMaterial::Union{Modia3D.AbstractMassPropertiesInterface,AbstractString,SolidMaterial,Nothing} = nothing,
        massProperties::Union{Modia3D.AbstractMassProperties, Number, AbstractString, SolidMaterial, Nothing} = nothing,
        collision::Bool = false,
        contactMaterial::AbstractString = "",
        collisionSmoothingRadius::Union{Float64, Nothing}=0.0,
        visualMaterial::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial(),
        visualMaterialConvexDecomposition::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial())

        if collision && isnothing(shape)
            error("For collision/gripping simulations, a shape must be defined.")
        end

        if collision
            # it is a solid material
            if typeof(solidMaterial) == String && contactMaterial == ""
                contactMaterial = solidMaterial
            elseif typeof(solidMaterial) != String && contactMaterial == ""
                error("For collision handling an Object3Ds solidMaterial must be defined as a string or a contactMaterial must be set.")
            end
        else # no collision, therefore contactMaterial is ""
            contactMaterial = ""
        end

        if typeof(solidMaterial) == String
            solidMaterial = solidMaterialPalette[solidMaterial]
        end

        if typeof(visualMaterial) == String
            visualMaterial = Shapes.visualMaterialPalette[visualMaterial]
        end

        if typeof(shape) == FileMesh
            (shape.centroid, shape.longestEdge, shape.objPoints, shape.facesIndizes) = getObjInfos(shape.filename, shape.scaleFactor)
            (shape.volume, shape.centroidAlgo, shape.inertia) = computeMassProperties(shape.objPoints, shape.facesIndizes; bodyCoords=false)
        end

        massProperties = createMassProperties(massProperties, shape, solidMaterial)
        new(shape, solidMaterial, massProperties, collision, contactMaterial, setCollisionSmoothingRadius(shape, collisionSmoothingRadius), visualMaterial)
    end
end

function JSON.show_json(io::JSON.StructuralContext, s::JSON.CommonSerialization, solid::Solid)
   JSON.begin_object(io)
      JSON.show_pair(io, s, "shape"            , solid.shape)
      JSON.show_pair(io, s, "massProperties" , solid.massProperties)
      JSON.show_pair(io, s, "material"       , solid.visualMaterial)
      JSON.show_pair(io, s, "contactMaterial", solid.contactMaterial)
   JSON.end_object(io)
end


text:
export mmul

"""
    mmul(w, x; dims=1)

Generalizes matrix multiplication to possibly more than 2 dims with reshapes: `w(A...,B...)
* x(B...,C...) => y(A...,C...)` where `size(w)=(A...,B...)` etc. and `dims=length(B)`. The last
`dims` dimensions of `w` has to match the first `dims` dimensions of `x`.
"""
function mmul(w, x; dims=1)
    if ndims(w) > 2 || ndims(x) > 2 || dims > 1
        @assert ndims(w) > dims "ndims(w) must be > dims"
        # ndims(x) may be <= dims in which case assume trailing dims assumed 1
        ntuple(i->size(x,i), dims) === size(w)[end-dims+1:end] || throw(DimensionMismatch("mmul: w=$(size(w)) x=$(size(x)) dims=$dims"))
        w2size = (prod(size(w)[1:end-dims]), prod(size(w)[end-dims+1:end]))
        w2 = (size(w) === w2size ? w : reshape(w, w2size))
        x2 = (size(x,1) === size(w2,2) && ndims(x) <= 2 ? x : reshape(x, size(w2,2), :))
        y2 = w2 * x2
        ysize1 = (dims === ndims(w) ? (1,) : size(w)[1:end-dims])
        ysize2 = (dims >= ndims(x) ? () : size(x)[dims+1:end])
        ysize = (ysize1..., ysize2...)
        y = (size(y2) === ysize ? y2 : reshape(y2, ysize))
    else # handle common matmul case without reshapes
        @assert dims === 1 "dims must be >= 1"
        size(x,1) === size(w,2) || throw(DimensionMismatch("mmul: w=$(size(w)) x=$(size(x)) dims=$dims"))
        y = w * x
    end
    return y
end

##
# In RNNs where input has shape (X,B,T) we want the result to be (Y,B,T). In CNNs where
# input has shape (X1,X2,C,B) we want the output to be (Y,B). So in general we want to be
# able to specify how many initial dims of x should be replaced. Does it ever get replaced
# with more than one dim? In my transformer implementation it does.

ext: jl


text:
struct Pose3 <: IncrementalInference.InferenceVariable
  dims::Int
  labels::Vector{String}
  Pose3() = new(6, String["POSE";])
end
struct Point3 <: IncrementalInference.InferenceVariable
  dims::Int
  Point3() = new(3)
end


text:
using GZip

CHUNK_LENGTH = 1024 * 1024 * 128

function write_to_file(accession::String, data::Vector{UInt8})
    GZip.open("output/pfam_split/$(accession).shm.gz", "w") do f
        num_chunks  = Int64(ceil(length(data) / CHUNK_LENGTH))
        for n = 1:num_chunks
            begin_ = ((n - 1) * CHUNK_LENGTH) + 1
            if n === num_chunks
                write(f, data[begin_:end])
            else
                write(f, data[begin_:(n * CHUNK_LENGTH)])
            end
        end
    end
end

GZip.open("input/Pfam-A.rp55.gz") do f
    cur_data::IOBuffer = IOBuffer()
    cur_accession::Union{String, Nothing} = nothing

    for line in eachline(f, keep=true)
        if startswith(line, "# STOCKHOLM 1.0")
            data = take!(cur_data)
            if cur_accession !== nothing
                write_to_file(cur_accession, data)
            end
            cur_accession = nothing
        end
        if startswith(line, "#=GF AC")
            cur_accession = split(line)[3]
        end
        write(cur_data, line)
    end

    if cur_accession !== nothing
        write_to_file(cur_accession, take!(cur_data))
    end
end

max_stars_repo_name: ibivu/protein-glue

max_issues_repo_name: ibivu/protein-glue

max_forks_repo_name: ibivu/protein-glue


text:
using Documenter, Literate, Recombinase

@info "makedocs"

for filename in ["tutorial.jl", "digging_deeper.jl"]
    Literate.markdown(
        joinpath(@__DIR__, "src", filename),
        joinpath(@__DIR__, "src", "generated")
    )
end

makedocs(
   format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
   sitename = "Recombinase.jl",
   pages = [
        "index.md",
        "generated/tutorial.md",
        "interactive_interface.md",
        "generated/digging_deeper.md",
   ]
)

@info "deploydocs"
deploydocs(
    repo = "github.com/piever/Recombinase.jl.git",
)


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function DataFile(TSTART,TSTOP,Ts)
# ----------------------------------------------------------------------------------- #
# DataFile.jl was generated using the Kwatee code generation system.
# DataFile: Stores model parameters as key - value pairs in a Julia Dict()
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 02-26-2016 14:31:22
#
# Input arguments:
# TSTART  - Time start
# TSTOP  - Time stop
# Ts - Time step
#
# Return arguments:
# data_dictionary  - Data dictionary instance (holds model parameters)
# ----------------------------------------------------------------------------------- #

# Initial conditions (default, no initiator)
initial_condition_array = zeros(19)
initial_condition_array[1] = 0.0      # 1 zymosan
initial_condition_array[2] = 1.9      # 2 C4
initial_condition_array[3] = 0.322    # 3 C2
initial_condition_array[8] = 7.57     # 4 C3
initial_condition_array[14] = 0.195   # 5 C5
initial_condition_array[17] = 2.23    # 6 Factor H
initial_condition_array[18] = 0.417   # 7 C4BP

# Load the parameters from disk -
# parameter_array = readdlm("./best_solution.txt")

# JV: initialize a placeholder parameter array.
# The actual parameter ensemble is contained in the data subdir
# pc_array*.dat = number_of_parameters x number_of_trials (each col is a parameter soln)
# rank_array*.dat = number_of_trials x 1  (holds the Pareto rank for each col in pc_array*.dat)
parameter_array = zeros(28,1)

# ---------------------------- DO NOT EDIT BELOW THIS LINE -------------------------- #
data_dictionary = Dict();
data_dictionary["PARAMETER_ARRAY"] = parameter_array;
data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
# ----------------------------------------------------------------------------------- #
return data_dictionary;
end


text:
function get_closest_planar_contact(parent_pose, parent_bbox, child_pose, child_bbox)
    possible_face_pairs = [
        (face1, face2) for face1 = 1:length(S.BOX_SURFACE_IDS)
        for face2 = 1:length(S.BOX_SURFACE_IDS)
    ]
    corresponding_planar_contacts = [
        S.closestApproximatingContact(
            (parent_pose * S.getContactPlane(parent_bbox, S.BOX_SURFACE_IDS[face1])) \
            (child_pose * S.getContactPlane(child_bbox, S.BOX_SURFACE_IDS[face2])),
        ) for (face1, face2) in possible_face_pairs
    ]
    error = [pose_distance(pc.slack, PC.IDENTITY_POSE; r=5.0) for pc in corresponding_planar_contacts]
    best_contact_index = argmin(error)
    return possible_face_pairs[best_contact_index], corresponding_planar_contacts[best_contact_index]
end

export get_closest_planar_contact


text:
mutable struct LockedPriorityQueue{K,V}
    # TODO: find a better queue
    items::Vector{Pair{K,V}}
    nitems::Int
    @atomic minkey::K
    @atomic locked::Bool
end

struct MultiQueue{K,V}
    queues::Vector{LockedPriorityQueue{K,V}}
end

function LockedPriorityQueue{K,V}(buffersize::Integer) where {K,V}
    items = Vector{Pair{K,V}}(undef, buffersize)
    return LockedPriorityQueue{K,V}(items, 0, typemax(K), false)
end

function MultiQueue{K,V}() where {K,V}
    nqueues = 2 * Threads.nthreads()
    queues = [LockedPriorityQueue{K,V}(128) for _ in 1:nqueues]
    return MultiQueue{K,V}(queues)
end

function Base.trylock(pq::LockedPriorityQueue)
    if @atomic :monotonic pq.locked
        return false
    end
    _, ok = @atomicreplace pq.locked false => true
    return ok
end

function Base.unlock(pq::LockedPriorityQueue)
    @_assert pq.locked
    @atomic pq.locked = false
    return
end

function _insert!(xs, n::Integer, i::Integer, x)
    @_assert n < length(xs)
    for j in n:-1:i
        @inbounds xs[j+1] = xs[j]
    end
    xs[i] = x
end

function trypush!(pq::LockedPriorityQueue{K,V}, x::Pair{K,V}) where {K,V}
    @_assert pq.locked
    @_assert first(x) < typemax(K)
    if length(pq.items) == pq.nitems
        # Main.@tlc lpq_full
        return false
    end
    i = searchsortedlast(view(pq.items, 1:pq.nitems), x; by = first, rev = true) + 1
    _insert!(pq.items, pq.nitems, i, x)
    pq.nitems += 1
    if i == 1
        @atomic :monotonic pq.minkey = first(x)
    end
    # Main.@tlc lpq_inserted
    return true
end

function maybepopmin!(pq::LockedPriorityQueue{K,V}) where {K,V}
    pq.minkey < typemax(K) || return nothing
    @_assert pq.locked
    @_assert pq.nitems > 0
    x = pq.items[pq.nitems]
    n = pq.nitems -= 1
    k = n == 0 ? typemax(K) : first(pq.items[n])
    @atomic :monotonic pq.minkey = k
    return Some{V}(last(x))
end

function trypush!(mq::MultiQueue{K,V}, x::Pair{K,V}) where {K,V}
    (; queues) = mq
    nfails = 0
    while true
        # TODO: cheaper RNG?
        a = rand(queues)
        if trylock(a)
            try
                trypush!(a, x) && return true
            finally
                unlock(a)
            end
        end
        nfails += 1
        if nfails > length(queues)  # TODO: what's the good threshold?
            # contention or full; not a good sign anyway
            return false
        end
        GC.safepoint()
    end
end

"""
    maybepopmin!(mq::MultiQueue{K,V}) -> Some(v::V) or nothing

Return approximately minimum element from the multi-queue `mq`; return `nothing`
if empty.
"""
function maybepopmin!(mq::MultiQueue{K,V}) where {K,V}
    (; queues) = mq
    while true
        # TODO: cheaper RNG?
        a = rand(queues)
        b = rand(queues)
        if a.minkey > b.minkey
            a = b
        end
        if a.minkey == typemax(K)
            found = false
            for outer a in queues
                if a.minkey < typemax(K)
                    found = true
                    break
                end
            end
            found || return nothing
        end
        if trylock(a)
            try
                y = maybepopmin!(a)
                y === nothing || return y
            finally
                unlock(a)
            end
        end
        GC.safepoint()
    end
end


text:
using PyPlot, PyCall
include("make_SQ_model.jl")
include("utilities.jl")

##### display plots to gui, or just save
## if they plot to gui, they don't have the size specified in teh code
pygui(true)

experiment_dict = Dict("n_vary"=>Dict(:inpts=>[:exit_distance],:epocs=>1000,:ex_locs=>[1.,5.],:cmp=>["ok","bad"],:legend_loc=>"upper right"),
                       "sense_vary"=>Dict(:inpts=>[:sensor_rwd],:epocs=>1000,:ex_locs=>[-35.,-150.],:cmp=>["ok","bad"],:legend_loc=>"lower right"),
                       "transition_vary"=>Dict(:inpts=>[:tprob],:epocs=>1000,:ex_locs=>[0.25,0.75],:cmp=>["ok","bad"],:legend_loc=>"lower right"),
                       "transition_e_vary"=>Dict(:inpts=>[:tprob,:e_mcts],:epocs=>1000,:ex_locs=>[],:cmp=>["ok","bad"],:legend_loc=>"lower right")
                      )
#  experiment_dict = Dict("transition_vary"=>Dict(:inpts=>[:tprob],:epocs=>1000,:ex_locs=>[0.25,0.75],:cmp=>["bad","ok"],:legend_loc=>"lower right"))
#  experiment_dict = Dict("n_vary"=>Dict(:inpts=>[:exit_distance],:epocs=>1000,:ex_locs=>[1.,5.],:cmp=>["bad"],:legend_loc=>"upper right"))
experiment_dict = Dict("transition_e_vary"=>Dict(:inpts=>[:tprob,:e_mcts],:epocs=>1000,:ex_locs=>[1.,5.],:cmp=>["ok","bad"],:legend_loc=>"upper right"))
#  experiment_dict = Dict("sense_vary"=>Dict(:inpts=>[:sensor_rwd],:epocs=>1000,:ex_locs=>[-35.,-150.],:cmp=>["bad"],:legend_loc=>"lower right"))

for expr in keys(experiment_dict)
    net_type = expr
    inpts = experiment_dict[expr][:inpts]
    epocs = experiment_dict[expr][:epocs]
    sq_example_locations = experiment_dict[expr][:ex_locs]
    cmp_list = experiment_dict[expr][:cmp]
    legend_loc = experiment_dict[expr][:legend_loc]
    println("Plotting: $(experiment_dict[expr])")
    for compare in cmp_list
        println("running '$compare' data")
        train_fname = "logs_old/$(net_type)_reference_solver_training.csv"
        test_fname = "logs_old/$(net_type)_$(compare)_solver.csv"

        inputs = Dict()
        for i in inpts
            inputs[i] = "ML.Continuous"
        end
        outputs = Dict(:X3_1=>"ML.Continuous",:X3_2=>"ML.Continuous")

        log_fname = "$(net_type)_$(make_label_from_keys(inputs))"
        log_loc = "nn_logs_old/"
        #  log_path = string(log_loc,"/",log_fname,"_net1-$(epocs)

        #  println("##########")
        #  println("$log_fname, $(readdir(log_loc))")
        #  println("##########")
        if !(any([contains(x,string(log_fname,"_")) for x in readdir(log_loc)]))
            println("No nn file exists, making one now")
            #  make_sq_model(net_type,inpts)
            make_sq_model(net_type,inpts,num_epoc=epocs)
        end

        param_files = searchdir(log_loc,log_fname,".params")

        num_epocs = parse(split(match(r"-\d+",param_files[1]).match,"-")[2])

        SQmodel = load_network(string(log_loc,log_fname),num_epocs,string(log_loc,log_fname,"_SQmodel.jld"))

        # get test and make predictions
        test_input, test_output, test_table, input_sch, output_sch = return_data(test_fname, inputs=inputs, outputs=outputs)

        data_mat = ML.featuremat(merge(input_sch,output_sch),test_table)

        info("restoring limits")
        limits = restore_eng_units(SQmodel.range,SQmodel.output_sch)
        info("getting predictions")
        _notused, pred_outputs = SQ_predict(SQmodel,test_input,use_eng_units=true)

        info("restoring test data")
        tst_in_eng = restore_eng_units(test_input,input_sch)
        tst_out_eng_ary = restore_eng_units(test_output,output_sch)

        # make figures
        if length(inputs) == 1
            fig,ax_ary = PyPlot.subplots(1,1,sharex=false)
            fig[:set_size_inches](8.0,6.0)
            fontsize = 20
            PyPlot.grid()

            i1 = collect(keys(inputs))[1]

            idx1 = nearest_to_x(tst_in_eng[i1],sq_example_locations[1])
            idx2 = nearest_to_x(tst_in_eng[i1],sq_example_locations[2])

            scatter_with_conf_bnds(ax_ary,tst_in_eng,tst_out_eng_ary,i1,:X3_1,:X3_2,:red,subsample=[idx1 idx2],label="candidate($(compare))",bar=true)
            #  scatter_with_conf_bnds(ax_ary,tst_in_eng,tst_out_eng_ary,i1,:X3_1,:X3_2,:red,subsample=collect(1:length(tst_in_eng[i1])),label="candidate",bar=true)
            #  ax_ary[:scatter](tst_in_eng[i1],tst_out_eng_ary[:X3_1],color=:red,alpha=0.2)

            ax_ary[:set_xlabel](string(i1),fontsize=fontsize)
            ax_ary[:set_ylabel](string("Reward"),fontsize=fontsize)
            ax_ary[:axhline](limits[:X3_1][2])
            ax_ary[:axhline](limits[:X3_1][1])

            add_sq_annotation(ax_ary,tst_in_eng,tst_out_eng_ary,pred_outputs,idx1,i1,:X3_1,:X3_2,SQmodel,fontsize=fontsize)
            add_sq_annotation(ax_ary,tst_in_eng,tst_out_eng_ary,pred_outputs,idx2,i1,:X3_1,:X3_2,SQmodel,fontsize=fontsize)

            ax_ary[:text](minimum(tst_in_eng[i1]),limits[:X3_1][2],L"r_H",fontsize=fontsize,va="bottom")
            ax_ary[:text](minimum(tst_in_eng[i1]),limits[:X3_1][1],L"r_L",fontsize=fontsize,va="top")

            #  ax_ary[:scatter](tst_in_eng[i1],pred_outputs[:X3_1],color=:blue)
            scatter_with_conf_bnds(ax_ary,tst_in_eng,pred_outputs,i1,:X3_1,:X3_2,:blue,subsample=collect(1:8:length(tst_in_eng[i1])),label="trusted",bar=false)
            ax_ary[:legend](loc=legend_loc,fontsize=fontsize)

            PyPlot.savefig(string("figs/",log_fname,"_",compare,".pdf"),dpi=300,transparent=true)
        elseif length(inputs) > 1
            # import stuff for plotting
            @pyimport mpl_toolkits.axes_grid1 as mpl_kit
            al = mpl_kit.make_axes_locatable

            function trust_surf(x,y,key)
                _notused, pred_outputs = SQ_predict(SQmodel,[x y]',use_eng_units=true)
                return pred_outputs[key][1]
            end

            i_rng = [-1.65;1.71]
            grid_ary = collect(linspace(i_rng[1],i_rng[2],50))
            ts = [trust_surf(x,y,:X3_1) for x in grid_ary, y in grid_ary]'
            tss = [trust_surf(x,y,:X3_2) for x in grid_ary, y in grid_ary]'

            display(ts)
            display(tss)

            # Solver Quality Surface
            fig = plt[:figure](figsize=(6.0,6.0))
            grid = fig[:subplots](2,2)
            div = []
            push!(div,al(grid[1]))
            push!(div,al(grid[2]))
            cax = []
            push!(cax,div[1][:append_axes]("right",size="5%",pad=0.05))
            push!(cax,div[2][:append_axes]("right",size="5%",pad=0.05))

            #  grid = ImageGrid(fig, 111,          # as in plt.subplot(111)
                                 #  nrows_ncols=(3,2),
                                 #  direction="row",
                                 #  aspect = true,
                                 #  axes_pad=0.65,
                                 #  share_all=true,
                                 #  cbar_location="right",
                                 #  cbar_mode="each",
                                 #  cbar_size="3%",
                                 #  cbar_pad=0.10,
                                 #  )
            fontsize = 10

            # make correlation plots
            i1 = inpts[1]
            i2 = inpts[2]

            #  corrplot(data_mat)
            poi1 = [0.2;200.] # point of interest, where X3 will be calculated
            poi2 = [0.7;800.] # point of interest, where X3 will be calculated
            poi1_norm = [(poi1[1]-mean(input_sch[i1]))/std(input_sch[i1]);(poi1[2]-mean(input_sch[i2]))/std(input_sch[i2])]
            poi2_norm = [(poi2[1]-mean(input_sch[i1]))/std(input_sch[i1]);(poi2[2]-mean(input_sch[i2]))/std(input_sch[i2])]

            x_rng = [i_rng[1]*std(input_sch[i1])+mean(input_sch[i1]);i_rng[2]*std(input_sch[i1])+mean(input_sch[i1])]
            y_rng = [i_rng[1]*std(input_sch[i2])+mean(input_sch[i2]);i_rng[2]*std(input_sch[i2])+mean(input_sch[i2])]
            subsample_num = 3

            aspct = diff(x_rng)[1]./diff(y_rng)[1]
            #  error()

            sq_cmap = :viridis
            im1 = grid[1][:imshow](ts,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[x_rng;y_rng],aspect="auto")
            #  im1 = grid[1][:imshow](ts,origin="lower",cmap=sq_cmap,interpolation="bilinear",aspect="auto")
            grid[1][:scatter](poi1[1],poi1[2],color=:black,marker=L"$A$",color=:gray)
            grid[1][:scatter](poi2[1],poi2[2],color=:black,marker=L"$B$",color=:gray)
            grid[1][:set_ylabel]("$(inpts[2])")
            grid[1][:set_title]("Mean Rwd: Trusted Solver")
            fig[:colorbar](im1,cax=cax[1])

            sqs_cmap = :magma
            im2 = grid[2][:imshow](tss,origin="lower",cmap=sqs_cmap,interpolation="bilinear",extent=[x_rng;y_rng],aspect="auto")
            #  im2 = grid[2][:imshow](tss,origin="lower",cmap=sqs_cmap,interpolation="bilinear",aspect="auto")
            grid[2][:scatter](poi1[1],poi1[2],color=:black,s=50,marker=L"$A$",color=:gray)
            grid[2][:scatter](poi2[1],poi2[2],color=:black,s=50,marker=L"$B$",color=:gray)
            grid[2][:set_ylabel]("$(inpts[2])")
            grid[2][:set_xlabel]("$(inpts[1])")
            grid[2][:set_title]("Std: Trusted Solver")
            fig[:colorbar](im2,cax=cax[2])
            #  grid[2][:set_aspect](aspct)

            add_sq_scatter3d_annotation(grid[1],grid[3],test_input,tst_out_eng_ary,i1,i2,poi1_norm,SQmodel,marker=L"$A$",s=50,fontsize=12)
            add_sq_scatter3d_annotation(grid[2],grid[4],test_input,tst_out_eng_ary,i1,i2,poi2_norm,SQmodel,marker=L"$B$",s=50,fontsize=12)

            grid[3][:legend]()

            PyPlot.tight_layout()
            PyPlot.subplots_adjust(hspace=0.4,wspace=0.8)

            PyPlot.savefig(string("figs/",log_fname,"_",compare,".pdf"),dpi=300,transparent=true)
        else
         error("can't support more than 2 inputs yet")
        end
        PyPlot.close_figs()
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
# Extract variables alone a line across multiple frames.
#
# Usage:
#   julia -t 4 demo_lineextracting_mt.jl
# or
#   JULIA_NUM_THREADS=4 julia demo_lineextracting_mt.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob
using JLD2: jldsave

files = glob("bulk*.vlsv", ".")
nfile = length(files)
# Define end points of the line in Earth radius
x1, x2 = 6.0, 20.0
point1 = [x1, 0, 0] .* Vlasiator.RE
point2 = [x2, 0, 0] .* Vlasiator.RE
# Extract cell info along the line
cellids, distances, coords =
   load(files[1]) do meta
      getcellinline(meta, point1, point2)
   end
# WARNING: this may not be exact due to round-off errors in output time stamps!
t = let tstart = load(files[1]).time, tend = load(files[end]).time
   range(round(tstart,digits=1), round(tend, digits=1), length=nfile)
end

n = zeros(Float32, length(cellids), nfile)
v = zeros(Float32, 3, length(cellids), nfile)
p = zeros(Float32, 6, length(cellids), nfile)
b = similar(v)
e = similar(v)

Threads.@threads for i in eachindex(files)
   println("i = $i/$nfile, file = $(files[i])")
   local meta = load(files[i])

   n[:,i] = readvariable(meta, "proton/vg_rho", cellids)
   v[:,:,i] = readvariable(meta, "proton/vg_v", cellids)
   p[1:3,:,i] = readvariable(meta, "proton/vg_ptensor_diagonal", cellids)
   p[4:6,:,i] = readvariable(meta, "proton/vg_ptensor_offdiagonal", cellids)
   b[:,:,i] = readvariable(meta, "vg_b_vol", cellids)
   e[:,:,i] = readvariable(meta, "vg_e_vol", cellids)
end

# Save into binary file
jldsave("vars_sun_earth_line.jld2"; x1, x2, t, n, v, p, b, e)

println("Line extraction finished!")


text:
@testset "$TEST $G" begin

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
z = dfs_tree(g5, 1)

@test ne(z) == 3
@test nv(z) == 4
@test !has_edge(z, 1, 3)

@test topological_sort_by_dfs(g5) == [1, 2, 3, 4]

# TODO add back
# @test !has_cycles(g5)
# @test is_tree(g5)

# g = DG(3)
# add_edge!(g,1,2); add_edge!(g,2,3); add_edge!(g,3,1)
# @test has_cycles(g)
# @test !is_tree(g)

# g = G(4)
# add_edge!(g,1,2); add_edge!(g,2,3); add_edge!(g,2,4);
# @test !has_cycles(g)
# @test is_tree(g)

# add_edge!(g,1,3);
# @test has_cycles(g)
# @test !is_tree(g)

end # testset


text:
@testset "eva.jl" begin
    @testset "computeparamfunction(covariates)" begin
        # function returned for empty covariates should not modify θ
        θ = [1.0]
        fsame = Extremes.computeparamfunction(Variable[])
        θtransformed = fsame(θ)

        @test θ == θtransformed

        # function returned for not empty covariates should modify θ
        θ = [1.0, 1.0]
        ev = Variable("t", collect(1:100))
        ftrans = Extremes.computeparamfunction([ev])
        θtransformed = ftrans(θ)

        @test θ != θtransformed

    end

    x₁ = Variable("x₁", [1])
    x₂ = Variable("x₂", [2])
    x₃ = Variable("x₂", [3])

    μ = 1 + x₂.value[] + x₃.value[]
    ϕ = -.5 + x₁.value[]
    ξ = .1

    y = [6]

    pd = GeneralizedExtremeValue(μ, exp(ϕ), ξ)

    model = BlockMaxima(Variable("y", y), locationcov=[x₂; x₃], logscalecov = [x₁])

    @testset "loglike(model, θ)" begin

        # Test BlockMaxima with known values
        pd = GeneralizedExtremeValue(μ, exp(ϕ), ξ)
        model = BlockMaxima(Variable("y", y), locationcov=[x₂; x₃], logscalecov = [x₁])
        @test Extremes.loglike(model, [1; 1; 1; -.5; 1; .1]) ≈ logpdf(pd,y[])

        # Test ThresholdExceedance with known values
        pd = GeneralizedPareto(exp(ϕ), ξ)
        model = ThresholdExceedance(Variable("y", y), logscalecov = [x₁])
        @test Extremes.loglike(model, [-.5; 1; .1]) ≈ logpdf(pd,y[])

    end

    @testset "quantile(model, θ, p)" begin
        #  p outside of [0, 1] throws
        @test_throws AssertionError quantile(model, [-.5; 1; .1], -1)

        # Test BlockMaxima with known values
        pd = GeneralizedExtremeValue(μ, exp(ϕ), ξ)
        model = BlockMaxima(Variable("y", y), locationcov=[x₂; x₃], logscalecov = [x₁])
        @test Extremes.quantile(model, [1; 1; 1; -.5; 1; .1],.99)[] ≈ quantile(pd,.99)

        # Test ThresholdExceedance with known values
        pd = GeneralizedPareto(exp(ϕ), ξ)
        model = ThresholdExceedance(Variable("y", y), logscalecov = [x₁])
        @test Extremes.quantile(model, [-.5; 1; .1], .99)[] ≈ quantile(pd,.99)

    end

    @testset "validatestationarity(model)" begin
        n = 1000
        y = collect(1:n)
        ev = [Variable("x", rand(n))]
        smodel = nothing

        # non-stationary BlockMaxima model
        model = Extremes.BlockMaxima(Variable("y", y), locationcov = ev, logscalecov = ev, shapecov = ev)
        @test_logs (:warn, "covariates cannot be included in the model when estimating the
            paramters by the probability weighted moment parameter estimation.
            The estimates for the stationary model is returned.") smodel = Extremes.validatestationarity(model)

        @test length(smodel.location.covariate) == 0
        @test length(smodel.logscale.covariate) == 0
        @test length(smodel.shape.covariate) == 0

        # non-stationary ThresholdExceedance model
        model = Extremes.ThresholdExceedance(Variable("y", y), logscalecov = ev, shapecov = ev)
        @test_logs (:warn, "covariates cannot be included in the model when estimating the
            paramters by the probability weighted moment parameter estimation.
            The estimates for the stationary model is returned.") smodel = Extremes.validatestationarity(model)

        @test length(smodel.logscale.covariate) == 0
        @test length(smodel.shape.covariate) == 0

        # stationary BlockMaxima model
        model = Extremes.BlockMaxima(Variable("y", y))
        @test_logs smodel = Extremes.validatestationarity(model)

        @test length(smodel.location.covariate) == 0
        @test length(smodel.logscale.covariate) == 0
        @test length(smodel.shape.covariate) == 0

        # stationary ThresholdExceedance model
        model = Extremes.ThresholdExceedance(Variable("y", y))
        @test_logs smodel = Extremes.validatestationarity(model)

        @test length(smodel.logscale.covariate) == 0
        @test length(smodel.shape.covariate) == 0

    end

    @testset "Base.show(io, obj)" begin
        # print BlockMaxima does not throw
        model = Extremes.BlockMaxima(Variable("y", collect(1:100)))
        buffer = IOBuffer()
        @test_logs Base.show(buffer, model)

        # print ThresholdExceedance does not throw
        model = Extremes.ThresholdExceedance(Variable("y", collect(1:100)))
        buffer = IOBuffer()
        @test_logs Base.show(buffer, model)

    end

    @testset "showparamfun(name, param)" begin
        f(β) = identity(β)
        # stationary param
        param = Extremes.paramfun(Variable[], f)

        @test Extremes.showparamfun("μ", param) == "μ ~ 1"

        # non-stationary param
        n = 100
        param = Extremes.paramfun([Variable("t", collect(1:n)), Variable("x", rand(n))], f)

        @test Extremes.showparamfun("μ", param) == "μ ~ 1 + t + x"

    end

    @testset "validatelength(length, explvariables)" begin
        n = 100
        evn = Variable("tn", collect(1:n))
        evnx = Variable("xn", rand(n))
        evn10 = Variable("tn10", collect(1:n+10))
        # empty explvariables shold not throw
        @test_logs Extremes.validatelength(n, Variable[])

        # explvariables length all == length should not throw
        @test_logs Extremes.validatelength(n, [evn, evnx])

        # first explanatory variable length != length throws
        @test_throws AssertionError Extremes.validatelength(n, [evn10, evn])

        # second explanatory variable length != length throws
        @test_throws AssertionError Extremes.validatelength(n, [evn, evn10])

    end

    include(joinpath("eva", "blockmaxima_test.jl"))
    include(joinpath("eva", "thresholdexceedance_test.jl"))

end


text:
module SamplingUtils

md(m, σ) = Dict("m"=>m, "σ"=>σ)

p_rng=Dict("x"=>50
    ,"ηA"=>0
    ,"ηF"=>0
    ,"muA"=>0
    ,"muF"=>0
    ,"s2"=> md(1.0, 1.5)
    ,"kR5"=>1.0
    ,"kr2"=>1.0
    ,"a15"=>md(1.0,1.5)
    ,"a25"=>md(0.6,1.5)
    ,"a35"=>md(0.3,1.5)
    ,"a5"=>md(2.0,0.5)
    ,"b15"=>md(0.5,1.5)
    ,"b25"=>md(1.0,1.5)
    ,"b35"=>md(2.0,1.5)
    ,"b5"=>md(1.0,0.5)
    ,"a18"=>md(1.0,1.5)
    ,"a28"=>md(0.8,1.5)
    ,"a38"=>md(0.6,1.5)
    ,"a48"=>md(0.3,1.5)
    ,"a8"=>md(2.0,0.5)
    ,"b18"=>md(0.5,1.5)
    ,"b28"=>md(1.0,1.5)
    ,"b38"=>md(2.0,1.5)
    ,"b48"=>md(3.0,1.5)
    ,"b8"=>md(1.0,0.5)
    ,"l05"=>md(2.0, 1.0)
    ,"l15"=>0.010
    ,"l25"=>0.010
    ,"l35"=>0.010
    ,"kP5"=>1.0
    ,"K"=>1000.0
    ,"maxA"=>2.3
    ,"maxF"=>2.1
    ,"s9"=> md(1.0, 1.5)
    ,"kR8"=>1.0
    ,"kR8"=>1.0
    ,"kr9"=>1.0
    ,"l08"=>md(2.0, 1.0)
    ,"l18"=>0.010 * 1
    ,"l28"=>0.010 * 1
    ,"l38"=>0.010 * 1
    ,"l48"=>0.010 * 1
    ,"kP8"=>1.0 * 1
    ,"gA"=>1.0
    ,"gF"=>1.0
    ,"gc"=>10.0
    ,"s05"=>md(0.01,1.5)
    ,"s5"=> md(1.0*10, 1.5)
    ,"s08"=>md(0.01,1.5)
    ,"s8"=> md(1.0*10, 1.5)
    ,"K5A"=>md(1.0,0.5)
    ,"n5A"=>[6.0, 6.00]
    ,"K58"=>md(0.1, 1.5)
    ,"n58"=>[6.0, 6.00]
    ,"K85"=>md(0.1, 1.5)
    ,"n85"=>[6.0, 6.00]
    ,"K8F"=>md(1.0,0.5)
    ,"n8F"=>[6.0, 6.00]
    ,"K9F"=>0.99 * 6.00
    ,"n9F"=>[6.0, 6.00]
    ,"K2A"=>md(1.0, 0.5)
    ,"n2A"=>[6.0, 6.00]
    ,"K25"=>md(0.1, 1.5)
    ,"n25"=>[6.0, 6.00]
    ,"K98"=>md(0.1, 1.5)
    ,"n98"=>[6.0, 6.00]
    ,"y"=>1)

end


text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays, Test

if VERSION >= v"1.2.0-DEV.229"
    sumsz(img) = Base.dims2string(size(img)) * ' '
else
    sumsz(img) = ""
end

const rrstr = VERSION >= v"1.6.0-DEV.1083" ? "reshape, " : ""
rrdim(n) = VERSION >= v"1.6.0-DEV.1083" ? n-1 : n

# N0f8 is shown as either N0f8 or Normed{UInt8, 8}
# RGB is shown as ColorTypes.RGB or RGB
function typestring(::Type{T}) where T
    buf = IOBuffer()
    show(buf, T)
    String(take!(buf))
end
N0f8_str = typestring(N0f8)
N0f16_str = typestring(N0f16)
RGB_str = typestring(RGB)

@testset "show" begin
    rgb32 = rand(RGB{Float32}, 3, 5)
    v = view(rgb32, 2:3, :)
    @test summary(v) == "2×5 view(::Array{RGB{Float32},2}, 2:3, :) with eltype $(RGB_str){Float32}"
    a = channelview(rgb32)
    @test summary(a) == (VERSION >= v"1.6.0-DEV.1083" ? "3×3×5 reinterpret(reshape, Float32, ::Array{RGB{Float32},2}) with eltype Float32" :
                                                        "3×3×5 reinterpret(Float32, ::Array{RGB{Float32},3})")
    num64 = rand(3,5)
    b = colorview(RGB, num64)
    str = summary(b)
    @test occursin("5-element", str) && occursin("reinterpret", str) && occursin("reshape", str) && occursin("RGB{Float64}", str) &&
          occursin("::$(typeof(num64))", str) && occursin("with eltype", str)
    rgb8 = rand(RGB{N0f8}, 3, 5)
    c = rawview(channelview(rgb8))
    str = summary(c)
    @test occursin("3×3×5 rawview(reinterpret($(rrstr)", str) && occursin("N0f8", str) &&
          occursin("::Array{RGB{N0f8},$(rrdim(3))}", str) && occursin("with eltype UInt8", str)
    @test summary(rgb8) == "3×5 Array{RGB{N0f8},2} with eltype $(RGB_str){$(N0f8_str)}"
    rand8 = rand(UInt8, 3, 5)
    d = normedview(PermutedDimsArray(rand8, (2,1)))
    @test summary(d) == "5×3 normedview(N0f8, PermutedDimsArray(::$(typeof(rand8)), (2, 1))) with eltype $(N0f8_str)"
    e = PermutedDimsArray(normedview(rand8), (2,1))
    str = summary(e)
    @test occursin("5×3 PermutedDimsArray(reinterpret", str) && occursin("N0f8", str) &&
          occursin("::$(typeof(rand8))), (2, 1)", str) && occursin("with eltype $(N0f8_str)", str)
    rand16 = rand(UInt16, 3, 5)
    f = PermutedDimsArray(normedview(N0f16, rand16), (2,1))
    str = summary(f)
    @test occursin("5×3 PermutedDimsArray(reinterpret", str) && occursin("N0f16", str) &&
          occursin("::$(typeof(rand16))), (2, 1)", str) && occursin("with eltype $(N0f16_str)", str)
    g = channelview(rgb8)
    etstr = VERSION >= v"1.6.0-DEV.1083" ? " with eltype" : ""
    str = summary(g)
    @test occursin("3×3×5 reinterpret($(rrstr)", str) && occursin("N0f8", str) &&
          occursin("::Array{RGB{N0f8},$(rrdim(3))}", str)
    VERSION >= v"1.6.0-DEV.1083" && @test occursin("with eltype", str)
    h = OffsetArray(rgb8, -1:1, -2:2)
    @test summary(h) == "$(sumsz(h))OffsetArray(::Array{RGB{N0f8},2}, -1:1, -2:2) with eltype $(RGB_str){$(N0f8_str)} with indices -1:1×-2:2"
    i = channelview(h)
    str = summary(i)
    @test occursin("$(sumsz(i))reinterpret($(rrstr)", str) && occursin("N0f8", str) &&
          occursin("OffsetArray(::Array{RGB{N0f8},$(rrdim(3))}", str) && occursin("-1:1, -2:2", str) && occursin("with indices", str)
    c = channelview(rand(RGB{N0f8}, 2))
    o = OffsetArray(c, -1:1, 0:1)
    str = summary(o)
    @test occursin("$(sumsz(o))OffsetArray(reinterpret($(rrstr)", str) && occursin("N0f8,", str) &&
          occursin("::Array{RGB{N0f8},$(rrdim(2))}", str) && occursin("-1:1, 0:1", str) && occursin("with eltype $(N0f8_str)", str) &&
          occursin("with indices", str)
    # Issue #45
    a = collect(tuple())
    @test summary(a) == "0-element $(typeof(a))"
    b = view(a, :)
    @test summary(b) == "0-element view(::$(typeof(a)), :) with eltype Union{}"
end

nothing


text:
export readgw, readedgelist, writegw, writeedgelist

"""
Return line after skipping if pred(line) is true
Default pred: after skipping empty lines and commented lines
"""
function readskipping(fd::IO,
                      pred=line -> length(line)==0 || line[1]=='#' || line[1]=='%')
    while !eof(fd)
        line = strip(readline(fd))
        if pred(line); continue; end
        return line
    end
    return ""
end

"""
    readgw(fd::IO)
    readgw(file::AbstractString) -> SparseMatrixCSC, node list

Reads LEDA format file describing a network. Outputs an undirected network.
An example of a LEDA file is in the examples/ directory.
"""
function readgw(fd::IO)
    line = readskipping(fd)
    strip(line)=="LEDA.GRAPH" || error("Error in line: $line")
    for i = 1:3; readskipping(fd); end
    nverts = parse(Int,readskipping(fd))
    vertices = Array{String}(nverts)
    for i = 1:nverts
        line = readskipping(fd)
        vertices[i] = match(r"\|{(.*?)}\|",line).captures[1]
    end
    nedges = parse(Int,readskipping(fd))
    I = Vector{Int}(nedges)
    J = Vector{Int}(nedges)
    for i = 1:nedges
        line = readskipping(fd)
        caps = match(r"(\d+) (\d+) 0 \|{.*}\|",line).captures
        n1 = parse(Int,caps[1])
        n2 = parse(Int,caps[2])
        n1==n2 && continue
        I[i] = n1
        J[i] = n2
    end
    G = sparse(vcat(I,J), vcat(J,I), 1, nverts, nverts, max)
    return Network(G,vertices)
end
readgw(file::AbstractString, args...) = open(fd -> readgw(fd,args...), file, "r")

"""
    readedgelist(fd::IO; header=false)
    readedgelist(file::AbstractString; header=false) -> SparseMatrixCSC, node list

Read list of edges and output undirected network
"""
function readedgelist(fd::IO; header=false)
    nodes = Set{String}()
    edges = Vector{Tuple{String,String}}()
    header && readline(fd)
    iter = 1
    while !eof(fd)
        line = readskipping(fd)
        isempty(line) && continue
        vals = split(strip(line))
        length(vals)!=2 && error("Error reading file at line $iter. Formatting error?")
        push!(edges,(vals[1],vals[2]))
        push!(nodes,vals[1])
        push!(nodes,vals[2])
        iter += 1
    end
    nodes = collect(nodes)
    nodesnum = indexmap(nodes)
    I = Vector{Int}()
    J = Vector{Int}()
    for vals in edges
        n1 = nodesnum[vals[1]]
        n2 = nodesnum[vals[2]]
        n1==n2 && continue
        push!(I, n1)
        push!(J, n2)
    end
    nverts = length(nodes)
    G = sparse(vcat(I,J), vcat(J,I), 1, nverts, nverts, max)
    Network(G,nodes)
end
readedgelist(file::AbstractString; args...) =
    open(x -> readedgelist(x;args...), file, "r")

"""
    writeedgelist(fd::IO, st::Network; prefix="",suffix="")
    writeedgelist(file::AbstractString, st::Network; prefix="",suffix="")

Write network to file as list of edges.

- `prefix`,`suffix` : Prefix and suffix to each line.
"""
function writeedgelist(fd::IO, st::Network; prefix="",suffix="")
    G = st.G
    nodes = st.nodes
    m,n = size(G)
    for col = 1:n
        for j = nzrange(G,col)
            row = G.rowval[j]
            if row < col
                println(fd, prefix, nodes[row], " ", nodes[col], suffix)
            end
        end
    end
end
writeedgelist(file::AbstractString, args...) =
    open(fd -> writeedgelist(fd, args...), file, "w")

"""
    writegw(fd::IO, st::Network)
    writegw(file::AbstractString, st::Network)

Write undirected network to file as LEDA format.
"""
function writegw(fd::IO, st::Network)
    G = st.G
    gnodes = st.nodes
    print(fd,"LEDA.GRAPH\nvoid\nvoid\n-2\n")
    m,n = size(G)
    println(fd,n)
    for node in gnodes
        println(fd,"|{",node,"}|")
    end
    println(fd,div(nnz(G),2))
    for col = 1:n
        for j = nzrange(G,col)
            row = G.rowval[j]
            if row < col
                println(fd,row," ",col," 0 |{}|")
            end
        end
    end
end
writegw(filename::AbstractString, args...) =
    open(fd -> writegw(fd,args...), filename,"w")


text:
using Flux, Base.Test

@testset "Flux" begin

include("utils.jl")
include("tracker.jl")
include("layers/normalisation.jl")

end


text:
module TestB
using ThinModules

@thinmod module TestA
   include("./c.jl")
   include("./b.jl")

end
end


text:
module SignalScope

include("scope.jl")
include("random.jl")
include("unet.jl")
include("transform.jl")

end # module


text:
module TextureManager

using ModernGL
using Images
using Colors
using FileIO
using MatrixMath

using JLGEngine.GraphicsManager
using ..GLLists
using ..GLDebugControl
using ..GLExtendedFunctions

""" TODO """
type GLTextureBlock
	anchor::Array{GLuint,1}

	GLTextureBlock() = new(GLuint[1])
end

""" TODO """
type GLTexture
	anchor::Array{GLuint,1}
	id::GLuint
	typ::GLenum

	GLTexture() = new(GLuint[1],0,GL_TEXTURE_2D)
end

""" TODO """
function create()
	this = GLTexture()
	glGenTextures(1, this.anchor)
	this.id = this.anchor[1]
	this
end

""" TODO """
bind(this::GLTexture) = glBindTexture(this.typ, this.id)

""" TODO """
function load(this::GLTexture, file::AbstractString)
	img = Images.load(file)
	sz = size(img)
	width = sz[1] #1
	height = sz[2] #1
	a = channelview(img)
	a = reinterpret.(vec(a))
	#a = [0xFF, 0xFF, 0xFF, 0xFF]
	internalformat = GL_RGBA8
	format = GL_RGBA
	typ = GL_UNSIGNED_BYTE

	glActiveTexture(GL_TEXTURE0)
	glBindTexture(this.typ, this.id)
	glTexParameteri(this.typ, GL_TEXTURE_MIN_FILTER, GL_NEAREST) #GL_LINEAR
	glTexParameteri(this.typ, GL_TEXTURE_MAG_FILTER, GL_NEAREST) #GL_LINEAR
	glTexParameteri(this.typ, GL_TEXTURE_WRAP_S, GL_REPEAT) #GL_CLAMP_TO_EDGE)
	glTexParameteri(this.typ, GL_TEXTURE_WRAP_T, GL_REPEAT) #GL_CLAMP_TO_EDGE)
	glTexImage2D(this.typ, 0, internalformat, width, height, 0, format, typ, pointer(a))
	#glGenerateMipmap(this.typ)
end

end #TextureManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
module Tools

# Exporting after each definition

"""
    print_fields(x, l_max=500, padding=15, padfunc=lpad)

Print all fieldnames of `x`, followed their value.
"""
function print_fields(x, l_max=500, padding=15, padfunc=lpad)
    fields = tofn(x)
    for field in fields
        second_part = "$(getfield(x, field))"
        second_part = length(second_part) < l_max ? second_part : "More than $l_max letters."
        println(padfunc("$field = ", padding) *  second_part)
    end
end
export print_fields

"""
    batch(v::AbstractVector, n_batches::Int, shuffle_input=false; check_even=true, return_indices=false)

Create a vector of `n_batches` vectors, containing all of `v`.
"""
function batch(v::AbstractVector, n_batches::Int, shuffle_input=false; check_even=true, return_indices=false)
    @assert length(v) ≥ n_batches "Trying to make $n_batches batches from $(length(v)) elements. This would result in empty arrays of type `Any`, which is likely to cause problems."
    shuffle_input  &&  (v = shuffle(v))
    check_even  &&  length(v) % n_batches != 0   &&   @warn "Number of elements not divisible by number of batches. Batches will be uneven. Set `check_even` to false to silence this warning"
    divs, rems = divrem(length(v), n_batches)
    batchlengths = fill(divs, n_batches)
    batchlengths[end-rems+1:end] .+= 1

    cumsums = pushfirst!(cumsum(batchlengths), 0)
    ranges = [cumsums[i]+1:cumsums[i+1] for i in 1:n_batches]
    if return_indices
        return (ranges, getindex.([v], ranges))
    else
        return getindex.([v], ranges)
    end
end
export batch

"""
    ⊕(args...) = 1/sum(x->1/x, args)

Compute the reciprocal of sum of reciprocal.
Can be used as infix operator.
"""
⊕(args...) = 1/sum(x->1/x, args)
export ⊕

"""
    tofn(x) = x |> typeof |> fieldnames
"""
tofn(x) = x |> typeof |> fieldnames
export tofn

"""
    moving_avg(v::AbstractVector, n)

Iterate over `v`, computing the average of the current and
next n-1 (n total) elements, returning it as a vector.
"""
function moving_avg(v::AbstractVector, n)
    output = Vector{Float64}(undef, length(v)-(n-1))
    for i in eachindex(output)
        output[i] = sum(v[i:i+(n-1)])/(n)
    end
    return output
end
export moving_avg

end


text:
"""
    abstract type AbstractIsotropicDiffusivity <: AbstractTurbulenceClosure end

Abstract type for closures with *isotropic* diffusivities.
"""
abstract type AbstractIsotropicDiffusivity{TD} <: AbstractTurbulenceClosure{TD} end

"""
    viscosity(closure, diffusivities)

Returns the isotropic viscosity associated with `closure`.
"""
function viscosity end

"""
    diffusivity(closure, diffusivities, tracer_indedx) where tracer_index

Returns the isotropic diffusivity associated with `closure` and tracer index `c_idx`.
"""
function diffusivity end

#####
##### Viscosity-stress products
#####

"""
    ν_σᶜᶜᶜ(i, j, k, grid, ν, σᶜᶜᶜ, u, v, w)

Multiply the array `ν` located at `ᶜᶜᶜ` by a function

    `σᶜᶜᶜ(i, j, k, grid, u, v, w)`

at index `i, j, k` and location `ᶜᶜᶜ`.
"""
@inline ν_σᶜᶜᶜ(i, j, k, grid, clock, ν, σᶜᶜᶜ, args...) = νᶜᶜᶜ(i, j, k, grid, clock, ν) * σᶜᶜᶜ(i, j, k, grid, args...)
@inline ν_σᶠᶠᶜ(i, j, k, grid, clock, ν, σᶠᶠᶜ, args...) = νᶠᶠᶜ(i, j, k, grid, clock, ν) * σᶠᶠᶜ(i, j, k, grid, args...)
@inline ν_σᶠᶜᶠ(i, j, k, grid, clock, ν, σᶠᶜᶠ, args...) = νᶠᶜᶠ(i, j, k, grid, clock, ν) * σᶠᶜᶠ(i, j, k, grid, args...)
@inline ν_σᶜᶠᶠ(i, j, k, grid, clock, ν, σᶜᶠᶠ, args...) = νᶜᶠᶠ(i, j, k, grid, clock, ν) * σᶜᶠᶠ(i, j, k, grid, args...)

#####
##### Stress divergences
#####

const AID = AbstractIsotropicDiffusivity

@inline viscous_flux_ux(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), Σ₁₁, U.u, U.v, U.w)
@inline viscous_flux_uy(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), Σ₁₂, U.u, U.v, U.w)
@inline viscous_flux_uz(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), Σ₁₃, U.u, U.v, U.w)

@inline viscous_flux_vx(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶠᶠᶜ(i, j, k, grid, clock, viscosity(closure, args...), Σ₂₁, U.u, U.v, U.w)
@inline viscous_flux_vy(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), Σ₂₂, U.u, U.v, U.w)
@inline viscous_flux_vz(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), Σ₂₃, U.u, U.v, U.w)

@inline viscous_flux_wx(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), Σ₃₁, U.u, U.v, U.w)
@inline viscous_flux_wy(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), Σ₃₂, U.u, U.v, U.w)
@inline viscous_flux_wz(i, j, k, grid, closure::AID, clock, U, args...) = - 2 * ν_σᶜᶜᶜ(i, j, k, grid, clock, viscosity(closure, args...), Σ₃₃, U.u, U.v, U.w)

#####
##### Diffusive fluxes
#####

@inline diffusive_flux_x(i, j, k, grid, closure::AID, c, c_idx, clock, args...) = diffusive_flux_x(i, j, k, grid, clock, diffusivity(closure, c_idx, args...), c)
@inline diffusive_flux_y(i, j, k, grid, closure::AID, c, c_idx, clock, args...) = diffusive_flux_y(i, j, k, grid, clock, diffusivity(closure, c_idx, args...), c)
@inline diffusive_flux_z(i, j, k, grid, closure::AID, c, c_idx, clock, args...) = diffusive_flux_z(i, j, k, grid, clock, diffusivity(closure, c_idx, args...), c)

#####
##### Support for VerticallyImplicitTimeDiscretization
#####

const VITD = VerticallyImplicitTimeDiscretization
const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded}

@inline z_viscosity(closure::AID, args...) = viscosity(closure, args...)
@inline z_diffusivity(closure::AID, c_idx, args...) = diffusivity(closure, c_idx, args...)

@inline ivd_viscous_flux_uz(i, j, k, grid, closure, clock, U, args...) = - ν_σᶠᶜᶠ(i, j, k, grid, clock, viscosity(closure, args...), ∂xᶠᶜᵃ, U.w)
@inline ivd_viscous_flux_vz(i, j, k, grid, closure, clock, U, args...) = - ν_σᶜᶠᶠ(i, j, k, grid, clock, viscosity(closure, args...), ∂yᶜᶠᵃ, U.v)

# General functions (eg for vertically periodic)
@inline viscous_flux_uz(i, j, k, grid,  ::VITD, closure::AID, args...) = ivd_viscous_flux_uz(i, j, k, grid, closure, args...)
@inline viscous_flux_vz(i, j, k, grid,  ::VITD, closure::AID, args...) = ivd_viscous_flux_vz(i, j, k, grid, closure, args...)
@inline viscous_flux_wz(i, j, k, grid,  ::VITD, closure::AID, args...) = zero(eltype(grid))
@inline diffusive_flux_z(i, j, k, grid, ::VITD, closure::AID, clock, args...) = zero(eltype(grid))

# Vertically bounded grids
#
# For vertically bounded grids, we calculate _explicit_ fluxes on the boundaries,
# and elide the implicit vertical flux component everywhere else. This is consistent
# with the formulation of the tridiagonal solver, which requires explicit treatment
# of boundary contributions (eg boundary contributions must be moved to the right
# hand side of the tridiagonal system).
@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AID, args...)
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  ivd_viscous_flux_uz(i, j, k, grid, closure, args...))
end

@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid, ::VITD, closure::AID, args...)
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  ivd_viscous_flux_vz(i, j, k, grid, closure, args...))
end

@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::AID, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end

@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::AID, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1,
                  diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end


text:
module QuAlgorithmZoo

using LinearAlgebra
using Yao, BitBasis
using YaoExtensions

include("Adam.jl")
include("PhaseEstimation.jl")
include("hamiltonian_solvers.jl")
include("HadamardTest.jl")
include("number_theory.jl")

@deprecate random_diff_circuit variational_circuit

end # module


text:
# Test each op and gradients, along with some combinations
# (TODO) Needs a much cleaner design long term
using cg
using Base.Test
using Base

show_test(test::AbstractString) = println("==============\n Running tests $test\n==============")
i = () -> cg.placeholder([1]) # Shape doesn't matter yet
c = cg.constant

function genarg(shape::AbstractArray, filltype::Symbol ; offset = 0, rand_offset = -0.5)
    if filltype == :range
        # val[j] = i + j
        return reshape(collect(Float64, 1:(*(shape...))), shape...) + offset
    elseif filltype == :rand
        return rand(shape...) + rand_offset
    elseif filltype == :ones
        return ones(shape...)
    else
        @assert false && "invalid filltype"
    end
end

function get_inputs(n::cg.Node)
    G = cg.get_graph([n])
    inputs = filter(x -> isa(x.op, cg.Placeholder), G.nodes)
    collect(inputs)
end

function test_output(out::cg.Node, truth::Function, args::Tuple{cg.Node, Array{Float64, 2}}...)
    # Float64/2 should be TensorValue once type system updated
    session = cg.Session(out)
    for (node, value) in args
        session.values[node] = value
    end
    graphresult = cg.interpret(session, out)
    truthresult = truth([x[2] for x in args]...)
    @test_approx_eq_eps(graphresult, truthresult, 0.001)
end

# Test each operator
function test_gradients(out::cg.Node, shape=[5] ; filltype::Symbol=:range, iters=1, debug=false)
    inputs = get_inputs(out)
    gradients = cg.grad(out, inputs)

    session = cg.Session(out)
    for iter = 1:iters
        for (i,input) = enumerate(inputs)
            session.values[input] = genarg(shape, filltype, offset = i)
        end

        for (input,grad) = zip(inputs, gradients)
            numeric = cg.numeric_grad(session, out, input, 0.0000001)
            symbolic = cg.interpret(session, grad, debug=debug)
            if debug
                @show session.values[input]
                @show session.values[out]
                @show numeric
                @show symbolic
            end
            @test_approx_eq_eps(symbolic, numeric, 0.01)
            for i = 1:length(symbolic)
                @test_approx_eq_eps(symbolic[i], numeric[i], 0.0001)
            end
        end
    end
end

i = () -> cg.placeholder([1]) # Shape doesn't matter yet
function test_scalar_gradients()
    show_test("scalar")
    @show test_gradients(sum(i()))

    # Scalar ops
    @show test_gradients(sum(i() + i()))
    @show test_gradients(sum(i() - i()))
    @show test_gradients(sum(i() * i()))
    @show test_gradients(sum(i() / i()))
    @show test_gradients(sum(i() ^ cg.constant(3.0)))
    @show test_gradients(sum(cg.constant(3.0) ^ i()))

    @show test_gradients(sum(-i()))
    @show test_gradients(sum(sign(i())))
    @show test_gradients(sum(sign(i())), filltype=:rand, iters=100)
    @show test_gradients(sum(exp(i())))
    @show test_gradients(sum(log(i())))
    @show test_gradients(sum(sin(i())))
    @show test_gradients(sum(cos(i())))
    @show test_gradients(sum(abs(i())))

    @show test_gradients(sum(max(i(), i())))
    @show test_gradients(sum(max(i(), i())), filltype=:rand, iters=100)
    @show test_gradients(sum(min(i(), i())))
    @show test_gradients(sum(min(i(), i())), filltype=:rand, iters=100)

    @show test_gradients(sum(cg.sigmoid(i())))
end

function test_shape_gradients()
    show_test("shape")
    # Shape related ops
    shape = [10,15]
    ops = [cg.maximum, cg.sum]
    for op in ops
        @show op
        for graph in [sum(op(i())),
                      sum(op(i(), cg.constant(1))), sum(op(i(), cg.constant(2)))]
            for filltype = [:ones, :rand, :range]
                @show (op, graph, filltype)
                test_gradients(graph, shape, filltype=filltype, iters=100)
            end
        end
    end
end

function test_get_and_set_gradients()
    c = cg.constant
    show_test("get_and_set")
    @show test_gradients(cg.getindex(i(), c(4)))
    @show test_gradients(cg.getindex(cg.setindex(i(), c(3), c(1000)), c(3)))
    @show test_gradients(cg.getindex(cg.setindex(i(), c(3), c(1000)), c(2)))
end

function test_other_gradients()
    show_test("other gradients")
    a = cg.constant(rand(4,5))
    b = i()
    # Mat mul
    @show test_gradients(cg.sum(cg.matmul(cg.t(b), b)), [4, 3])
    @show test_gradients(cg.sum(cg.matmul(b, cg.t(b))), [4, 3])
    @show test_gradients(cg.sum(cg.matmul(cg.t(a), b)), [4, 3])
    @show test_gradients(cg.sum(cg.matmul(cg.t(b), a)), [4, 3])
end

function test_broadcast()
    show_test("Broadcast")
    i1 = i()
    subScalar = c(2.0) - i()
    divScalar =  i() / c(2.0)
    subVec = cg.broadcastop(cg.Sub(), c([1.0, 2.0, 3.0]), i())
    divVec = cg.broadcastop(cg.Div(), i(), c([1.0, 2.0, 3.0]))
    #@show test_output(add, (x) -> [10.0, 10.0] .+ x, (i1, [1.0]'))
    for test in [subScalar, divScalar, subVec, divVec]
        for filltype in [:ones, :range, :rand]
            @show test
            @show filltype
            @show test_gradients(cg.sum(test), [3,5], filltype=filltype)
        end
    end
end

function test_nn()
    show_test("nn")

    for filltype in [:ones, :rand]
        si1 = i()
        softmax = cg.softmax(si1)
        ssi1 = i()
        softmax_stable = cg.softmax_stable(ssi1)
        ci1 = i()
        ci2 = i()
        crossentropy = cg.crossentropy(ci1, ci2)
        sci1 = i()
        sci2 = i()
        softmax_crossentropy = cg.crossentropy(sci1, cg.softmax(sci2))

        softmax_t(x) = exp(x) ./ sum(exp(x), 1)
        crossentropy_t(p, q) = -sum(p .* log(q))
        softmax_crossentropy_t(p, q) = crossentropy_t(p, softmax_t(q))
        shape = [5,10]
        a1 = genarg(shape, :rand, rand_offset=0)
        a2 = genarg(shape, :rand, rand_offset=0)

        @show test_output(softmax, softmax_t, (si1, a1))
        @show test_output(softmax_stable, softmax_t, (ssi1, a1))
        @show test_output(crossentropy, crossentropy_t, (ci1, a1), (ci2, a2))
        @show test_output(softmax_crossentropy, softmax_crossentropy_t, (sci1, a1), (sci2, a2))

        @show test_gradients(cg.sum(broadcast("/", si1, sum(si1, cg.constant(1)))))
        @show test_output(print(cg.sum(broadcast("/", si1, sum(si1)))), x -> sum(x / sum(x)), (si1, a1))

        @show test_gradients(cg.sum(broadcast("/", si1, sum(si1))))
        @show test_gradients(cg.getindex(si1, cg.constant(1)))
        @show test_gradients(sum(si1))

        scal1 = i()
        scal2 = i()
        @show test_gradients(scal1 / (scal1 + scal2), [1])
        @show test_gradients(cg.getindex(broadcast("/", si1, sum(si1)), cg.constant(1)), [2])

        @show test_gradients(cg.sum(softmax))
        @show test_gradients(cg.sum(softmax_stable))
        @show test_gradients(cg.sum(crossentropy))
        @show test_gradients(cg.sum(softmax_crossentropy))
    end
end

function test_sgd_basics()
    show_test("sgd")
    # Test basic optimization - minimize sum((b - a)^2)
    target = rand(10)
    a = cg.constant(target, "a")
    b = cg.variable(cg.constant(zeros(10)), "b") # Parameter to optimize
    c = b - a
    d = c * c
    e = sum(d)
    values = Dict{cg.Node, cg.TensorValue}()
    optimizer = cg.sgd_optimizer(e, [b], cg.constant(0.001, "step_size"))
    #cg.render(cg.get_graph([b]), "graph.png")
    session = cg.Session(optimizer)
    for i = 1:10000
        cg.interpret(session, optimizer)
        if session.values[e][1] <= 0.001
            break
        end
    end
    @test_approx_eq_eps target session.values[b] 0.1
end

# Test that sum gradients work properly
function test_sum()
    show_test("sum")
    @show test_gradients(cg.sum(i()), [3, 5])
    a = cg.placeholder([3, 5])
    b = cg.sum(a, cg.constant(1))
    c = cg.sum(b)
    @show test_gradients(c, [3, 5])
end

test_other_gradients()
test_sum()
test_broadcast()
test_nn()
test_scalar_gradients()
test_shape_gradients()
test_get_and_set_gradients()
test_sgd_basics()


text:
abstract type Animal end
speak(a::Animal) = "$(a.name) says $(sound(a))"

struct Horse <: Animal
  name
end
sound(h::Horse) = "neigh"

struct Cow <: Animal
  name
end
sound(c::Cow) = "moooo"

struct Sheep <: Animal
  name
end
sound(s::Sheep) = "baaaa"

s = Horse("CJ")
@assert speak(s) == "CJ says neigh"
c = Cow("Bessie")
@assert speak(c) == "Bessie says moooo"
@assert speak(Sheep("Little Lamb")) == "Little Lamb says baaaa"

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
type Token
    form::Int
    cat::Int
    head::Int
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
function SimJulia.run( mrs::MRS, showInfo::Bool=false; saveConfig::Bool=true,
    seed::Integer=-1, sysEnt::Bool=false )::Nothing
    if !verifySimulation!( mrs )
        error( "Simulation configuration is inconsistent, cannot run." )
        return
    elseif showInfo
        @info "Simulation configuration is consistent, initialising processes."
    end  # if !verifySimulation!( mrs )

    # Set the random seed of the seed generator.
    mrs.showInfo = showInfo
    mrs.seedRNG = seed < 0 ? MersenneTwister() : MersenneTwister( seed )

    # Initialise the database.
    initialiseMRSresultsDatabase!( mrs )

    # Create folder for tmp databases.
    tmpFolder = joinpath( dirname( mrs.resultsDB.file ), "tmps" )

    if !ispath( tmpFolder )
        mkpath( tmpFolder )
    end  # if !ispath( tmpFolder )

    nThreads = min( mrs.nRuns, mrs.maxThreads, Threads.nthreads() )
    hasSnapshot = mrs.mpSim.dbSize > 0

    if hasSnapshot
        # copySnapshots( mrs, tmpFolder )
        copySnapshots( mrs, tmpFolder, nThreads )
        GC.gc()
    end  # if hasSnapshot

    # Run the simulations.
    tStart = now()
    mrs.nComplete = 0
    currentSim = Channel{Int}( 1 )
    # resultsFree = Channel{Bool}( 1 )
    put!( currentSim, 1 )
    mrs.mpSim.sim = Simulation()
    # put!( resultsFree, true )

    Threads.@threads for threadnum in 1:nThreads
        # runsim( mrs, currentSim, resultsFree, tmpFolder )
        runsim( mrs, threadnum, hasSnapshot, currentSim, tmpFolder, nThreads )
    end  # Threads.@threads for threadnum in ...

    tElapsed = ( now() - tStart ).value / 1000.0
    mrs.showInfo && @info string( "Running ", mrs.nRuns, " replications on ",
        nThreads, " threads took ", tElapsed, " seconds." )

    # Clear tmp databases.
    copyResults( mrs, tmpFolder )
    run( mrs.mpSim )
    GC.gc();  # SQLite keeps files locked until their pointers have been
              #   scrubbed by the GC.
    rm( tmpFolder, recursive=true, force=true )
    return
end  #  run( mrs, showInfo, saveConfig, seed, sysEnt )

include( joinpath( simPrivPath, "multirunsimulation.jl" ) )


text:
#
# include("helpFuncs3D.jl"); n = [10,11,12]; loc1 = 3; loc2 = 4; loc3 = 5; (loc1a,loc2a,loc3a) = cs2loc3D(loc2cs3D(loc1,loc2,loc3,n),n)
#

export loc2cs3D,cs2loc3D

function loc2cs3D(loc::Array{Int64,1},n::Array{Int64,1})
@inbounds cs = loc[1] + (loc[2]-1)*n[1] + (loc[3]-1)*n[1]*n[2];
return cs;
end

function loc2cs3D(loc1::Int64,loc2::Int64,loc3::Int64,n::Array{Int64,1})
@inbounds cs = loc1 + (loc2-1)*n[1] + (loc3-1)*n[1]*n[2];
return cs;
end

function cs2loc3D(loc::Array{Int64,1},cs_loc::Int64,n::Array{Int64,1})
@inbounds loc[1] = mod(cs_loc-1,n[1]) + 1;
@inbounds loc[2] = div(mod(cs_loc-1,n[1]*n[2]),n[1]) + 1;
@inbounds loc[3] = div(cs_loc-1,n[1]*n[2])+1;
return;
end

function cs2loc3D(cs_loc::Int64,n::Array{Int64,1})
@inbounds loc1 = mod(cs_loc-1,n[1])+1;
@inbounds loc2 = div(mod(cs_loc-1,n[1]*n[2]),n[1]) + 1;
@inbounds loc3 = div(cs_loc-1,n[1]*n[2])+1;
return (loc1,loc2,loc3);
end

function get3DNeighbors(loc::Array{Int64,1},n::Array{Int64,1},ans::Array{Int64,2},relevant::Array{Bool,1})
# This function fills the neighbours array ans with neighbours the "C" style, so it doesn't allocate.
@inbounds relevant[1]= loc[1]>1;
@inbounds relevant[2]= loc[1]<n[1];
@inbounds relevant[3]= loc[2]>1;
@inbounds relevant[4]= loc[2]<n[2];
@inbounds relevant[5]= loc[3]>1;
@inbounds relevant[6]= loc[3]<n[3];
######################################################################################
@inbounds ans[1,1] = loc[1]-1;
@inbounds ans[1,2] = loc[2];
@inbounds ans[1,3] = loc[3];
@inbounds ans[2,1] = loc[1]+1;
@inbounds ans[2,2] = loc[2];
@inbounds ans[2,3] = loc[3];
@inbounds ans[3,1] = loc[1];
@inbounds ans[3,2] = loc[2]-1;
@inbounds ans[3,3] = loc[3];
@inbounds ans[4,1] = loc[1];
@inbounds ans[4,2] = loc[2]+1;
@inbounds ans[4,3] = loc[3];
@inbounds ans[5,1] = loc[1];
@inbounds ans[5,2] = loc[2];
@inbounds ans[5,3] = loc[3]-1;
@inbounds ans[6,1] = loc[1];
@inbounds ans[6,2] = loc[2];
@inbounds ans[6,3] = loc[3]+1;
return;
end

function setPaddedArr3D(arr::Array{Bool,3},loc1::Int64,loc2::Int64,loc3::Int64,val::Bool)
# this function assumes that the boundaries are OK and the value is placed in a legal location.
@inbounds arr[loc1+2,loc2+2,loc3+2] = val
end

function getPaddedArr3D(arr::Array{Bool,3},loc1::Int64,loc2::Int64,loc3::Int64)
# this function assumes that the boundaries are OK and the value is placed in a legal location.
@inbounds return arr[loc1+2,loc2+2,loc3+2]
end


text:
# Example of rotation of acquisition and excitation
using Koma
using JLD2
# Scanner
sys = Scanner()
sys.ADC_Δt = 4e-6
# SEQ design
B1 = 6e-6; durRF = π/(2π*γ*B1)
EX = PulseDesigner.RF_hard(B1, durRF, sys; G=[0,0,0])
EPI,_,_,_ = PulseDesigner.EPI_base(40/100, 101, sys)
TE = 25e-3
DELAY = Delay(TE-dur(EPI)/2-dur(EX))
#Rotation
Rz = Koma.rotz(-π/2) #counter-clockwise rotation
Ry = Koma.roty(-π/2) #counter-clockwise rotation
Rx = Koma.rotx(π/2)  #counter-clockwise rotation

seq = EX + DELAY + EPI*(Ry*Rz)
@save "./EPI_example.seq" seq=seq
plot_seq(seq)


text:
#f(p,x)=KL(p,softmax(x))

#TODO: implement KL(q,p)

function FKLsoftmax(p::Array{Float64,2},x::Array{Float64,2})
    logZ=log(sum(exp(x),1))
    return ([(sum(p.*(log(p)-x))+sum(logZ))/length(x)],logZ)
end

function FKLsoftmax_inplace(value,auxvalue,p::Array{Float64,2},x::Array{Float64,2})
    logZ=log(sum(exp(x),1))
    fill!(value,(sum(p.*(log(p)-x))+sum(logZ))/length(x))
    copy!(auxvalue,logZ)
end

function DKLsoftmax(derivativeIDX,f_c,faux_c,grad_c,grad_n,p::Array{Float64,2},x::Array{Float64,2})
    if derivativeIDX==1
        axpy!(1.0/length(x),grad_c.*(1+log(p)-x+repmat(faux_c,size(x,1),1)),grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,grad_c.*(softmax(x)-p)./length(x),grad_n)
    end
end

if PROC=="GPU"

    function FKLsoftmax(p::CudaArray,x::CudaArray)
        tmp=CudaArray(Float64,size(x))
        exp!(x,tmp)
        onr=CudaArray(Float64,(1,size(x,1))); fill!(onr,1.0)
        logZ=onr*tmp
        log!(logZ,logZ)
        log!(p,tmp)
        axpy!(-1.0,x,tmp)
        out=CudaArray(Float64,(1,1)); fill!(out,0.0)
        axpy!(1./length(x),tmp,out)
        axpy!(1./length(x),sum(logZ),out)
        free(tmp)
        return (out,logZ)
        #return ([(sum(p.*(log(p)-x))+sum(log(Z)))/length(x)],aux)
    end

    function FKLsoftmax_inplace(value,auxvalue,p::CudaArray,x::CudaArray)
        tmp=CudaArray(Float64,size(x))
        exp!(x,tmp)
        onr=CudaArray(Float64,(1,size(x,1))); fill!(onr,1.0);
        logZ=CudaArray(Float64,(1,size(x,2)))
        log!(onr*tmp,logZ)
        log!(p,tmp)
        axpy!(-1.0,x,tmp)
        vmult!(1.0,p,tmp,tmp)
        fill!(value,0.0)
        axpy!(1./length(x),sum(tmp),value)
        axpy!(1./length(x),sum(logZ),value)
        copy!(auxvalue,logZ)
        free(onr); free(logZ); free(tmp)
        #return ([(sum(p.*(log(p)-x))+sum(log(Z)))/length(x)],aux)
    end

    function DKLsoftmax(derivativeIDX,f_c,faux_c,grad_c,grad_n,p::CudaArray,x::CudaArray)
        if derivativeIDX==1
            #axpy!(1.0/length(x),grad_c.*(1+log(p)-x+repmat(faux_c,size(x,1),1)),grad_n)
            tmp=CudaArray(Float64,size(x)); fill!(tmp,1.0)
            logp=CudaArray(Float64,size(p))
            log!(p,logp)
            axpy!(1.0,logp,tmp)
            axpy!(-1.0,x,tmp)
            onc=CudaArray(Float64,size(x,1),1); fill!(onc,1.0)
            axpy!(1.0,onc*faux_c,tmp)
            alphaaxpy!(1./length(x),grad_c,tmp,grad_n)
            free(tmp); free(onc); free(logp)
        elseif derivativeIDX==2
            #axpy!(1.0,grad_c.*(softmax(x)-p)./length(x),grad_n)
            tmp=CudaArray(Float64,size(x))
            softmax!(x,tmp)
            axpy!(-1.0,p,tmp)
            alphaaxpy!(1.0/length(x),grad_c,tmp,grad_n)
            free(tmp)
        end
    end

end

Derivative[FKLsoftmax]=DKLsoftmax
Inplace[FKLsoftmax]=FKLsoftmax_inplace

KLsoftmax(p::ADnode,x::ADnode)=ADnode(FKLsoftmax,[p x])

export KLsoftmax

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function show(io::IO, h::HalfPlane{T}) where T
    compact = get(io, :compact, false)
    !compact && print(io, "HalfPlane{$T}(")
    print(io, "$(h.a)x + $(h.b)y + $(h.c) ≤ 0")
    !compact && print(io, ")")
end

function show(io::IO, inter::Intersection{T}) where T
    padding = get(io, :padding, 0) + 2
    print(io, "Intersection{$T}\n")
    for (i, h) in enumerate(inter.content)
        print(io, " "^padding)
        print(Base.IOContext(io, :compact=>true, :padding=>padding), h)
        i != length(inter.content) && print(io, "\n")
    end
end

function show(io::IO, reunion::Reunion{T}) where T
    padding = get(io, :padding, 0) + 2
    print(io, "Reunion{$T}\n")
    for (i, h) in enumerate(reunion.content)
        print(io, " "^padding)
        print(Base.IOContext(io, :compact=>true, :padding=>padding), h)
        i != length(reunion.content) && print(io, "\n")
    end
end

function show(io::IO, p::ConvexPolygon{T}) where T
    padding = get(io, :padding, 0) + 2
    print(io, "ConvexPolygon{$T} with $(nb_vertices(p)) vertices\n")
    for (i, c) in enumerate(vertices(p))
        print(io, " "^padding)
        print(Base.IOContext(io, :compact=>true, :padding=>padding), c)
        i != nb_vertices(p) && print(io, "\n")
    end
end


text:
using IntArrays

# build an internal structure of the wavelet matrix
function build(::Type{B}, data::AbstractVector{T}, w::Int, destructive::Bool) where {T<:Unsigned,B}
    ivec = IntVector{w}(data)
    if destructive
        # free memory
        empty!(data)
    end
    bits = Vector{B}(undef,w)
    _build!(B, ivec, bits)
    return tuple(bits...)
end

function _build!(::Type{B}, ivec, bits) where B
    w = length(bits)
    len = length(ivec)
    # allocate working space
    bv = BitVector(undef, len)
    ivec′ = similar(ivec)
    for d in 1:w
        # scan d-th bit
        nzeros = 0
        for i in 1:len
            if (ivec[i] >> (w - d)) & 1 == 1
                bv[i] = true  # right
            else
                bv[i] = false # left
                nzeros += 1
            end
        end
        # stably sort integers by bv
        r = nzeros
        l = 0
        for i in 1:len
            if bv[i]
                ivec′[r+=1] = ivec[i]
            else
                ivec′[l+=1] = ivec[i]
            end
        end
        # store the bit vector and go next
        bits[d] = bv
        copy!(ivec, ivec′)
    end
end

# starting points for the rank operation can be precomputed
function locate_sp(a, bits, nzeros)
    w = length(bits)
    sp = 0
    for d in 1:w
        bit = (a >> (w - d)) & 1 == 1
        sp = rank(bit, bits[d], sp)
        if bit
            sp += nzeros[d]
        end
    end
    return sp
end


text:
module Import_Test

using FactCheck, JNeuron

facts() do

    myimport=input(string(dirname(Base.source_path()),"/../examples/data/cell2.asc"));
    @fact length(myimport.sections) --> 198

    blank_neuron=instantiate(myimport);
    @fact length(blank_neuron.secs) --> 198

    blank_neuron2=instantiate(string(dirname(Base.source_path()),"/../examples/data/cell2.asc"));
    @fact length(blank_neuron2.secs) --> 198

    a=zeros(Int64,4)
    for i=1:length(blank_neuron.secs)
        a[blank_neuron.secs[i].mtype]+=1
    end

    @fact a --> [1, 37, 63, 97]

    @fact blank_neuron.secs[198].child --> [1, 38, 43, 62, 63, 80, 87, 96, 101]

end

end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
"""
    function _build_mccormick_relaxation!(m, x, y, z)

McCormick relaxation of binlinear term
```
z >= JuMP.lower_bound(x)*y + JuMP.lower_bound(y)*x - JuMP.lower_bound(x)*JuMP.lower_bound(y)
z >= JuMP.upper_bound(x)*y + JuMP.upper_bound(y)*x - JuMP.upper_bound(x)*JuMP.upper_bound(y)
z <= JuMP.lower_bound(x)*y + JuMP.upper_bound(y)*x - JuMP.lower_bound(x)*JuMP.upper_bound(y)
z <= JuMP.upper_bound(x)*y + JuMP.lower_bound(y)*x - JuMP.upper_bound(x)*JuMP.lower_bound(y)
```
"""
function _build_mccormick_relaxation!(
    m::JuMP.Model,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    z::JuMP.VariableRef,
)::FormulationInfo
    x_lb, x_ub = _variable_domain(x)
    y_lb, y_ub = _variable_domain(y)

    formulation_info = FormulationInfo()

    formulation_info.constraints[:lb_1] =
        @constraint(m, z >= x_lb * y + y_lb * x - x_lb * y_lb)
    formulation_info.constraints[:lb_2] =
        @constraint(m, z >= x_ub * y + y_ub * x - x_ub * y_ub)
    formulation_info.constraints[:ub_1] =
        @constraint(m, z <= x_lb * y + y_ub * x - x_lb * y_ub)
    formulation_info.constraints[:ub_2] =
        @constraint(m, z <= x_ub * y + y_lb * x - x_ub * y_lb)

    return formulation_info
end

"""
    _build_bilinear_relaxation!(m, x, y, z, x_partition, y_partition, pre_base_name)

Build incremental formulation for ``z = xy`` given partition data.
"""
function _build_bilinear_milp_relaxation!(
    m::JuMP.Model,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    z::JuMP.VariableRef,
    x_partition::Vector{<:Real},
    y_partition::Vector{<:Real},
    pre_base_name::AbstractString
)::FormulationInfo

    origin_vs, non_origin_vs = _collect_bilinear_vertices(x_partition, y_partition)
    formulation_info = FormulationInfo()
    num_vars = max(length(x_partition), length(y_partition)) - 1

    # add variables
    delta_1 =
        formulation_info.variables[:delta_1] =
            @variable(m, [1:num_vars],
                lower_bound = 0.0, upper_bound = 1.0,
                base_name = pre_base_name * "delta_1")
    delta_2 =
        formulation_info.variables[:delta_2] =
            @variable(m, [1:num_vars],
                lower_bound = 0.0, upper_bound = 1.0,
                base_name = pre_base_name * "delta_2")
    delta_3 =
        formulation_info.variables[:delta_3] =
            @variable(m, [1:num_vars],
                lower_bound = 0.0, upper_bound = 1.0,
                base_name = pre_base_name * "delta_3")
    z_bin = formulation_info.variables[:z_bin] =
        @variable(m, [1:num_vars], binary = true, base_name = pre_base_name * "z")

    # add x constraints
    formulation_info.constraints[:x] = @constraint(
        m,
        x ==
        origin_vs[1][1] + sum(
            delta_1[i] * (non_origin_vs[i][1] - origin_vs[i][1]) +
            delta_2[i] * (non_origin_vs[i+1][1] - origin_vs[i][1]) +
            delta_3[i] * (origin_vs[i+1][1] - origin_vs[i][1]) for i = 1:num_vars
        )
    )

    # add y constraints
    formulation_info.constraints[:y] = @constraint(
        m,
        y ==
        origin_vs[1][2] + sum(
            delta_1[i] * (non_origin_vs[i][2] - origin_vs[i][2]) +
            delta_2[i] * (non_origin_vs[i+1][2] - origin_vs[i][2]) +
            delta_3[i] * (origin_vs[i+1][2] - origin_vs[i][2]) for i = 1:num_vars
        )
    )

    # add z constraints
    formulation_info.constraints[:z_bin] = @constraint(
        m,
        z ==
        origin_vs[1][3] + sum(
            delta_1[i] * (non_origin_vs[i][3] - origin_vs[i][3]) +
            delta_2[i] * (non_origin_vs[i+1][3] - origin_vs[i][3]) +
            delta_3[i] * (origin_vs[i+1][3] - origin_vs[i][3]) for i = 1:num_vars
        )
    )

    # add first delta constraint
    formulation_info.constraints[:first_delta] =
        @constraint(m, delta_1[1] + delta_2[1] + delta_3[1] <= 1)

    # add linking constraints between delta_1, delta_2 and z
    formulation_info.constraints[:below_z] =
        @constraint(m, [i = 2:num_vars], delta_1[i] + delta_2[i] + delta_3[i] <= z_bin[i-1])
    formulation_info.constraints[:above_z] =
        @constraint(m, [i = 2:num_vars], z_bin[i-1] <= delta_3[i-1])

    return formulation_info
end


text:
using ROOTFramework, Cxx

bindings = TTreeBindings()
idx = bindings[:idx] = Ref(zero(Int32))
s = bindings[:s] = Ref(zero(Float64))
a = bindings[:a] = Float64[]

A = Float64[]

open(TChainInput, bindings, "data", "out.root") do input
    @info "Input has $(length(input)) entries."

    n = 0
    @time for _ in input
        n += 1
        @assert idx[] == n
        append!(A, a)
    end
end

@info length(A)

#=
using Plots
stephist(A)
=#


text:
"""
    SigmaRidgeRegressor(; decomposition, groups, σ, center, scale, init_model)

A MLJ model that fits σ-Ridge Regression with `groups` and parameter `σ`.
`center` and `scale` (default `true` for both) control whether the response and
features should be centered and scaled first (make sure that `center=true` if the
model is supposed to have an intercept!). `decomposition` can be one of `:default`,
`:cholesky` or `:woodbury` and determines how the linear system is solved. `init_model`
is the initial model used to define the Method of Moments map from `σ` to `λ`; it defaults
to leave-one-out ridge without groups.
"""
Base.@kwdef mutable struct SigmaRidgeRegressor{G,T,M} <:
                           SigmaRidgeRegression.AbstractGroupRidgeRegressor
    decomposition::Symbol = :default
    groups::G
    σ::T = 1.0
    center::Bool = true
    scale::Bool = true
    init_model::M = LooRidgeRegressor(;ridge=SingleGroupRidgeRegressor(
        decomposition = decomposition,
        groups = groups,
        center = center,
        scale = scale,)
    )
end

_main_hyperparameter(::SigmaRidgeRegressor) = :σ

function _default_hyperparameter_maximum(model::SigmaRidgeRegressor, fitted_machine)
    sqrt(σ_squared_max(fitted_machine.cache.mom))
end

_default_param_min_ratio(::SigmaRidgeRegressor, fitted_machine) = 1e-3
_default_scale(::SigmaRidgeRegressor, fitted_machine)  = :linear

function MMI.fit(m::SigmaRidgeRegressor, verb::Int, X, y)
    @unpack init_model, decomposition, center, scale, groups = m
    init_machine = MLJ.machine(init_model, X, y)
    fit!(init_machine; verbosity = verb)
    mom = MomentTunerSetup(init_machine.cache)

    σ = m.σ
    λs = SigmaRidgeRegression.get_λs(mom, abs2(σ))

    multiridge = MultiGroupRidgeRegressor(;
        groups=groups,
        λs=λs,
        decomposition = decomposition,
        center = center,
        scale = scale,
    )
    multiridge_machine = MLJ.machine(multiridge, X, y)
    fit!(multiridge_machine; verbosity = verb)

    workspace = multiridge_machine.cache
    cache = (workspace = workspace, mom = mom, multiridge_machine = multiridge_machine)

    βs = StatsBase.coef(workspace)
    x_transform = multiridge_machine.fitresult.x_transform
    y_transform = multiridge_machine.fitresult.y_transform
    fitresult = (coef = βs, x_transform = x_transform, y_transform = y_transform)

    # return
    return fitresult, cache, NamedTuple{}()
end

function MMI.update(
    model::SigmaRidgeRegressor,
    verbosity::Int,
    old_fitresult,
    old_cache,
    X,
    y,
)
    @unpack init_model, decomposition, center, scale, groups = model

    workspace = old_cache.workspace
    multiridge_machine = old_cache.multiridge_machine
    mom = old_cache.mom

    σ = model.σ
    λs = SigmaRidgeRegression.get_λs(mom, abs2(σ))
    multiridge = MultiGroupRidgeRegressor(;
        groups=groups,
        λs=λs,
        decomposition = decomposition,
        center = center,
        scale = scale,
    )
    multiridge_machine.model = multiridge
    fit!(multiridge_machine; verbosity = verbosity)

    cache = (workspace = workspace, mom = mom, multiridge_machine = multiridge_machine)

    βs = StatsBase.coef(workspace)
    x_transform = multiridge_machine.fitresult.x_transform
    y_transform = multiridge_machine.fitresult.y_transform
    fitresult = (coef = βs, x_transform = x_transform, y_transform = y_transform)
    return fitresult, cache, NamedTuple{}()
end

"""
    LooSigmaRidgeRegressor(; kwargs...)

Convenience constructors and type-alias for  `LooRidgeRegressor{<:SigmaRidgeRegressor}`.
Equivalent to `LooRidgeRegressor(;ridge= SigmaRidgeRegressor())` with `kwargs`
passed to both `LooRidgeRegressor` and `SigmaRidgeRegressor`.
"""
const LooSigmaRidgeRegressor =  LooRidgeRegressor{<:SigmaRidgeRegressor}

function LooSigmaRidgeRegressor(;
             tuning = DefaultTuning(),
             rng = Random.GLOBAL_RNG, kwargs...)
             sigma_ridge = SigmaRidgeRegressor(;kwargs...)
    LooRidgeRegressor(;ridge=deepcopy(sigma_ridge), tuning=tuning, rng=rng)
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Turing
using Dates
using StatsBase
using Distributions
using CovidSurvey
# using Gadfly
using StatsPlots
using LaTeXStrings
using NamedTupleTools
using StatsFuns
using Random
using TransformVariables, LogDensityProblems
using Optim
using ReverseDiff
using BSON
using TimeSeriesDecompositions
plotlyjs()

# using DataFramesMeta
# plotlyjs()
const dateform = "yyyymmdd"
include( projectdir("scripts/survey/rawsurvey_utils.jl") )

## ============================================================================
region = 4
rawcontacts = readcontacts( "/home/and/data/covidsurvey/rawcontacts.csv", select_region=region )
let
    gp = groupby(rawcontacts[:,[:date]], :date, sort=true)
    cb = combine(gp, nrow)
    @show minimum(cb.nrow)
    histogram(cb.nrow, bins=20)
end
## =============================================================================
# remove outliers
# pgfplotsx()
# scalefontsizes(1/1.5)
threshold = (family=50, colleagues=100, friends=100, strangers=1000)
contacts = filteroutliers(rawcontacts, threshold) #NOTE National: p = 99.9
p = contact_histograms(contacts; bins=0:40, fc=:gray)
# saveto = projectdir("figures","survey_histogram.png")
# savefig(p, saveto)
## =============================================================================
# plot daily aggregated time series
# plotlyjs()
p = plot_confidence_timeseries(contacts);

## =============================================================================
# zero contacts
zerocontacts = countzeros(contacts)
plot_timeseries(zerocontacts)

## =============================================================================
# inference with zero inflated model

surveyproblem = SurveyInferenceProblem(
    ZINegBinomialProblem,
    ctype,
    contacts
)

results = optimize(surveyproblem, LBFGS(),
    Optim.Options(show_trace=true, iterations=2000);
)
contactrates = DataFrame(results)

let
    fname = projectdir("scripts/survey/","survey_inference_results_method=$(results.problemtype)_region=$region.bson")
    save(fname, θ = results)
end
plot_confidence_timeseries(results)
## save and plots results
using FillArrays
BSON.@load projectdir("scripts/survey/survey_inference_results_method=ZINegBinomialProblem_region=$region.bson") results
# save(projectdir("data","inferred_contact_rates.csv"), contactsrates)
# plot_generated_quantities(results; compare_old=true)
pgfplotsx()
p = plot_survey_histograms(results, day=nothing)
savefig(p, projectdir("figures/histogram_survey_data_vs_fake_data.png"))
## load results
contactrates = load(projectdir("data","inferred_contact_rates.csv")) |> DataFrame
## =============================================================================
# decomposition
import TimeSeriesDecompositions as TSD

slowtrends = TSD.decompose(contactrates; fdir=projectdir("figures/"))
save(projectdir("data","smoothed_contact_rates_region=$region.csv"), slowtrends)
# bson("contacts_decomposition=$(results.problemtype)_region=$region", Dict(:slowtrends => slowtrends))
## load results
slowtrends = load(projectdir("data","smoothed_contact_rates_region=$region.csv"))
# BSON.load("survey_inference_results_method=ZINegBinomialProblem.bson")
##
p = plot_comparison_with_SSI(contactrates);
firefox(p; fname="compare_SSI_with_contactrates.html");

## =============================================================================

olddata = CSV.File("/home/and/data/covidsurvey/smoothed_contacts_0425_2m.csv") |> DataFrame
rename!(olddata, :date=>:dates)
p = plot_comparison(olddata, slowtrends);
firefox(p; fname="compare_contactrates_with_olddata.html");

olddata = CSV.File("/home/and/data/covidsurvey/contacts_0425_2m.csv") |> DataFrame
rename!(olddata, :date=>:dates)
p = plot_comparison(olddata, contactrates);
firefox(p; fname="compare_contactrates_with_olddata.html");

##

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# ~/work/julia/Iban/docs$ julia  make.jl

push!(LOAD_PATH,"../src/")

using IbanGen
using Documenter
using Memoize

makedocs(
    sitename = "IbanGen",
    modules = [IbanGen],
    authors = "David Soroko"
)

# see: https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
deploydocs(
    repo = "github.com/sorokod/Iban.jl.git"
)


text:
include("slide.jl")

include("slide_left.jl")
include("slide_right.jl")


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

"""
    initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)

Compressible Euler setup for a double Mach reflection problem.
Involves strong shock interactions as well as steady / unsteady flow structures.
Also exercises special boundary conditions along the bottom of the domain that is a mixture of
Dirichlet and slip wall.
See Section IV c on the paper below for details.

- Paul Woodward and Phillip Colella (1984)
  The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
  [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
"""
@inline function initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)

  if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3)
    phi = pi / 6
    sin_phi, cos_phi = sincos(phi)

    rho =  8
    v1  =  8.25 * cos_phi
    v2  = -8.25 * sin_phi
    p   =  116.5
  else
    rho = 1.4
    v1  = 0
    v2  = 0
    p   = 1
  end

  prim = SVector(rho, v1, v2, p)
  return prim2cons(prim, equations)
end

initial_condition = initial_condition_double_mach_reflection

boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_double_mach_reflection)

# Supersonic outflow boundary condition. Solution is taken entirely from the internal state.
# See `examples/p4est_2d_dgsem/elixir_euler_forward_step_amr.jl` for complete documentation.
@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
                                            surface_flux_function, equations::CompressibleEulerEquations2D)
  # NOTE: Only for the supersonic outflow is this strategy valid
  # Calculate the boundary flux entirely from the internal solution state
  return flux(u_inner, normal_direction, equations)
end

# Special mixed boundary condition type for the :Bottom of the domain.
# It is Dirichlet when x < 1/6 and a slip wall when x >= 1/6
@inline function boundary_condition_mixed_dirichlet_wall(u_inner, normal_direction::AbstractVector,
                                                         x, t, surface_flux_function,
                                                         equations::CompressibleEulerEquations2D)
  if x[1] < 1 / 6
    # From the BoundaryConditionDirichlet
    # get the external value of the solution
    u_boundary = initial_condition_double_mach_reflection(x, t, equations)
    # Calculate boundary flux
    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
  else # x[1] >= 1 / 6
    # Use the free slip wall BC otherwise
    flux = boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, equations)
  end

  return flux
end

boundary_conditions = Dict( :Bottom => boundary_condition_mixed_dirichlet_wall,
                            :Top    => boundary_condition_inflow,
                            :Right  => boundary_condition_outflow,
                            :Left   => boundary_condition_inflow   )

volume_flux = flux_ranocha
surface_flux = flux_lax_friedrichs

polydeg = 4
basis = LobattoLegendreBasis(polydeg)
shock_indicator = IndicatorHennemannGassner(equations, basis,
                                            alpha_max=0.5,
                                            alpha_min=0.001,
                                            alpha_smooth=true,
                                            variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

# Get the unstructured quad mesh from a file (downloads the file if not available locally)
default_mesh_file = joinpath(@__DIR__, "abaqus_double_mach.inp")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/a0806ef0d03cf5ea221af523167b6e32/raw/61ed0eb017eb432d996ed119a52fb041fe363e8c/abaqus_double_mach.inp",
                                      default_mesh_file)
mesh_file = default_mesh_file

mesh = P4estMesh{2}(mesh_file)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.2)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorLöhner(semi, variable=Trixi.density)

amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=0,
                                      med_level=3, med_threshold=0.05,
                                      max_level=6, max_threshold=0.1)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=1,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback)

# positivity limiter necessary for this example with strong shocks
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
                                                     variables=(Trixi.density, pressure))

###############################################################################
# run the simulation
sol = solve(ode, SSPRK43(stage_limiter!),
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary