text:
# test constraints
# MPI Declarations
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

# opts = PdePumiInterface.get_defaults()
# # 2D mesh
# opts["order"] = 1
# opts["dimensions"] = 2
# opts["use_DG"] = true
# opts["operator_type"] = "SBPOmega"
# opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
# opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
# opts["numBC"] = 2
#
# # For 2DAirfoil
# opts["BC1"] = [8,11,14,17]
# opts["BC1_name"] = "FarField"
# opts["BC2"] = [5]
# opts["BC2_name"] = "Airfoil"
#
# opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
# opts["jac_type"] = 2
# opts["jac_method"] = 2
# opts["run_type"] = 5
#
# # Create PumiMesh and SBP objects
# sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = PDESolver.createMeshAndOperator(opts, 1)

opts = Dict{ASCIIString, Any}()
opts["order"] = 1
opts["dimensions"] = 2
opts["use_DG"] = true
opts["Tsbp"] = Float64
opts["Tmsh"] = Complex128
opts["operator_type"] = "SBPOmega"
opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
opts["numBC"] = 2

# For 2DAirfoil
opts["BC1"] = [8,11,14,17]
opts["BC1_name"] = "FarField"
opts["BC2"] = [5]
opts["BC2_name"] = "Airfoil"

opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
opts["jac_type"] = 2
opts["jac_method"] = 2
opts["run_type"] = 5

# Create PumiMesh and SBP objects
dofpernode = 1
ref_verts = [-1. 1 -1; -1 -1 1]
Tsbp = opts["Tsbp"]
Tmsh = opts["Tmsh"]
sbp = getTriSBPOmega(degree=opts["order"], Tsbp=opts["Tsbp"])
sbpface = TriFace{Tsbp}(opts["order"], sbp.cub, ref_verts.')
topo = 0
shape_type = 2
mesh = PumiMeshDG2(Tmsh, sbp, opts, sbpface, dofpernode=dofpernode,
                     shape_type=shape_type)

#mesh = PumiMeshDG2{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"],
#                         sbp, opts, sbpface; dofpernode=dofpernode,
#                         coloring_distance=opts["coloring_distance"],
#                         shape_type=shape_type)

#geom_faces = opts["BC2"]
bc_nums = [2]
# Free Form deformation parameters
ndim = 3
order = [2,4,2]  # Order of B-splines in the 3 directions
nControlPts = [2,4,2]

# code here for creating rectilinear map
# Create Mapping object
map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh,
                            bc_nums=bc_nums)

# Create knot vector
calcKnot(map)

# Create Bounding box
offset = [0., 0., 0.5] # No offset in the X & Y direction
ffd_box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

# Control points
controlPoint(map, ffd_box)

# Populate map.xi
calcParametricMappingNonlinear(map, ffd_box, mesh, bc_nums)

@testset "--- Checking Linear-Plane Constraints ---" begin

  di = 2

  # test function for counting number of constraints
  numCnstr = numLinearPlaneConstraints(map, di)
  @test ( numCnstr )== 1*3*4

  # test function for counting variables in each constraint
  cntEq = zeros(Int, numCnstr)
  fncidx = 1
  fncidx = countVarsLinearPlaneConstraints!(map, di, cntEq, fncidx)
  @test ( fncidx )== numCnstr+1
  for k = 1:numCnstr
    @test ( cntEq[k] )== 4
  end

  # test function for filling sparse constraint Jacobian
  fncidx = 1
  ptr = 1
  iLfun = zeros(Int, sum(cntEq))
  jLvar = zeros(Int, sum(cntEq))
  LinG = zeros(Tmsh, sum(cntEq))
  Flow = rand(numCnstr) + 0im
  Fupp = rand(numCnstr) + 0im
  fncidx, ptr = setLinearPlaneConstraints!(map, di, iLfun, jLvar, LinG, Flow,
                                           Fupp, fncidx, ptr)
  @test ( fncidx )== numCnstr+1
  @test ( ptr )== sum(cntEq)+1
  @test isapprox( Flow, zeros(numCnstr), 1e-15)
  @test isapprox( Fupp, zeros(numCnstr), 1e-15)
  ptr = 1
  for k = 1:numCnstr
    for i = 1:cntEq[k]
      @test ( iLfun[ptr] )== k
      ptr += 1
    end
  end
  # !!!!!!!!!! should also include checks on jLvar and LinG, eventually
end

@testset "--- Checking Linear-Corner Constraints ---" begin

  di = 2

  # test function for counting number of constraints
  numCnstr = numLinearCornerConstraints(map, di)
  @test ( numCnstr )== 8

  # test function for counting variables in each constraint
  cntEq = zeros(Int, numCnstr)
  fncidx = 1
  fncidx = countVarsLinearCornerConstraints!(map, di, cntEq, fncidx)
  @test ( fncidx )== numCnstr+1
  for k = 1:numCnstr
    @test ( cntEq[k] )== 4
  end

  # test function for filling sparse constraint Jacobian
  fncidx = 1
  ptr = 1
  iLfun = zeros(Int, sum(cntEq))
  jLvar = zeros(Int, sum(cntEq))
  LinG = zeros(Tmsh, sum(cntEq))
  Flow = rand(numCnstr) + 0im
  Fupp = rand(numCnstr) + 0im
  fncidx, ptr = setLinearCornerConstraints!(map, di, iLfun, jLvar, LinG, Flow,
                                            Fupp, fncidx, ptr)
  @test ( fncidx )== numCnstr+1
  @test ( ptr )== sum(cntEq)+1
  @test isapprox( Flow, zeros(numCnstr), 1e-15)
  @test isapprox( Fupp, zeros(numCnstr), 1e-15)
  ptr = 1
  for k = 1:numCnstr
    for i = 1:cntEq[k]
      @test ( iLfun[ptr] )== k
      ptr += 1
    end
  end
  # !!!!!!!!!! should also include checks on jLvar and LinG, eventually
end

@testset "--- Checking Linear-Stretch Constraints ---" begin

  di = 2

  # test function for counting number of constraints
  numCnstr = numLinearStretchConstraints(map)
  @test ( numCnstr )== (2*4*2 - 1)

  # test function for counting variables in each constraint
  cntEq = zeros(Int, numCnstr)
  fncidx = 1
  fncidx = countVarsLinearStretchConstraints!(map, cntEq, fncidx)
  @test ( fncidx )== numCnstr+1
  for k = 1:numCnstr
    @test ( cntEq[k] )== 2
  end

  # test function for filling sparse constraint Jacobian
  fncidx = 1
  ptr = 1
  iLfun = zeros(Int, sum(cntEq))
  jLvar = zeros(Int, sum(cntEq))
  LinG = zeros(Tmsh, sum(cntEq))
  Flow = rand(numCnstr) + 0im
  Fupp = rand(numCnstr) + 0im
  fncidx, ptr = setLinearStretchConstraints!(map, di, iLfun, jLvar, LinG, Flow,
                                             Fupp, fncidx, ptr)

  @test ( fncidx )== numCnstr+1
  @test ( ptr )== sum(cntEq)+1
  @test isapprox( Flow, zeros(numCnstr), 1e-15)
  @test isapprox( Fupp, zeros(numCnstr), 1e-15)
  ptr = 1
  for k = 1:numCnstr
    for i = 1:cntEq[k]
      @test ( iLfun[ptr] )== k
      ptr += 1
    end
  end
  # !!!!!!!!!! should also include checks on jLvar and LinG, eventually
end

@testset "--- Checking Linear-Root Constraints ---" begin

  di = 2

  # test function for counting number of constraints
  numCnstr = numLinearRootConstraints(map)
  @test ( numCnstr )== 5

  # test function for counting variables in each constraint
  cntEq = zeros(Int, numCnstr)
  fncidx = 1
  fncidx = countVarsLinearRootConstraints!(map, cntEq, fncidx)
  @test ( fncidx )== numCnstr+1
  for k = 1:3
    @test ( cntEq[k] )== 1
  end
  @test ( cntEq[4] )== 2
  @test ( cntEq[5] )== 2

  # test function for filling sparse constraint Jacobian
  fncidx = 1
  ptr = 1
  iLfun = zeros(Int, sum(cntEq))
  jLvar = zeros(Int, sum(cntEq))
  LinG = zeros(Tmsh, sum(cntEq))
  Flow = rand(numCnstr) + 0im
  Fupp = rand(numCnstr) + 0im
  fncidx, ptr = setLinearRootConstraints!(map, di, iLfun, jLvar, LinG, Flow,
                                          Fupp, fncidx, ptr)
  @test ( fncidx )== numCnstr+1
  @test ( ptr )== sum(cntEq)+1
  ptr = 1
  for k = 1:numCnstr
    for i = 1:cntEq[k]
      @test ( iLfun[ptr] )== k
      ptr += 1
    end
  end
  for k = 1:3
    @test ( LinG[k] )== 1.0
  end
  @test ( LinG[4] )== 1.0
  @test ( LinG[5] )== -1.0
  @test ( LinG[6] )== 1.0
  @test ( LinG[7] )== -1.0
end


text:
using Base.Iterators, IterTools

struct Vocab
    v2i::Dict{String,Int16}
    i2v::Vector{String}
    unk::Int16
    eos::Int16
    tokenizer
end

function Vocab(file::String; tokenizer=split, vocabsize=Inf, mincount=1, unk="<unk>", eos="<s>")
    # set unk and eos tokens frequency to inf because
    # we don't want them to be removed from the vocab set
    cdict = Dict(eos => Inf, unk=>Inf)

    # create vocab set and count occurrences
    for l in eachline(file)
        tokens = tokenizer(l)
        map(v -> cdict[v] = get!(cdict, v, 0) + 1, tokens)
    end

    # select vocabs with frequency higher than mincount
    # sort by frequency and delete if vocabsize is determined
    fsorted = sort([ (v, c) for (v, c) in cdict if c >= mincount ], by = x -> x[2], rev = true)

    vocabsize == Inf || length(fsorted) < vocabsize || (fsorted = fsorted[1:vocabsize])

    i2v = [ eos; unk; [ x[1] for x in fsorted[3:end] ] ]
    v2i = Dict( v => i for (i, v) in enumerate(i2v))

    return Vocab(v2i, i2v, v2i[unk], v2i[eos], tokenizer)
end

struct TextReader
    file::String
    vocab::Vocab
end

function Base.iterate(r::TextReader, s=nothing)
    s === nothing && (s = open(r.file))
    eof(s) && return close(s)
    return [[ get(r.vocab.v2i, v, r.vocab.unk) for v in r.vocab.tokenizer(readline(s))] ; [ r.vocab.eos, ] ] , s
end

Base.IteratorSize(::Type{TextReader}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{TextReader}) = Base.HasEltype()
Base.eltype(::Type{TextReader}) = Vector{Int16}

mutable struct TextData
    src::TextReader             # reader for text data
    batchsize::Int              # desired batch size
    batchmajor::Bool            # batch dims (B,T) if batchmajor=false (default) or (T,B) if true.
    dataarray::Array{Int16, 2} # batchified data array for language modeling tasks
    maxsize::Int                # max length of text to read
    bptt::Int                   # how many steps to back propagate through time
    randomize::Bool             # randomize BPTT
end

function TextData(src::TextReader; batchsize = 128, maxlength = typemax(Int),
                batchmajor = false, maxsize = typemax(Int), bptt = 1024, randomize = true)
    datavector = Vector{Int16}()
    nr = lb = length(datavector)

    for l in src
        for i in l
            nr += 1
            if nr > lb
                lb = nr * 2
                resize!(datavector, lb)
            end

            datavector[nr] = i
            if nr > maxsize
                datavector = resize!(datavector, maxsize)
                @goto escape_label
            end
        end
    end
    @label escape_label
    datavector = resize!(datavector, nr)
    N = length(datavector) ÷ batchsize
    batchified_dataarray = reshape(datavector[1:N * batchsize], N, batchsize)'
    TextData(src, batchsize, batchmajor, batchified_dataarray, maxsize, bptt, randomize)
end

function changebatchsize!(d::TextData, newbatchsize::Int)
    (B, N) = size(d.dataarray)
    datavector = reshape(d.dataarray', N * B)
    N = length(datavector) ÷ newbatchsize
    d.dataarray = reshape(datavector[1:N * newbatchsize], N, newbatchsize)'
    d.batchsize = newbatchsize
    return newbatchsize
end

Base.IteratorSize(::Type{TextData}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{TextData}) = Base.HasEltype()
Base.eltype(::Type{TextData}) = Tuple{Array{Int16,2},Array{Int16,2}}

function Base.iterate(d::TextData, state=nothing)
    if state === nothing
        state = (1, d.bptt + 1)
    elseif state == 0
        return nothing
    end

    while true
        x = d.dataarray[:, state[1]:state[2]]
        (d.batchmajor) && (x = x')
        batch = (x[:, 1:end-1], x[:, 2:end])

        if state[2] >= size(d.dataarray, 2)
            return batch, 0
        elseif d.randomize
            seq_len = (rand([1,1,1,1,1,1,1,1,1,0]) == 1) ? d.bptt : (d.bptt ÷ 2) # get half of bptt 10% of the time.
            next_bptt = Int(gaussian(1; mean=seq_len, std=5)[1] ÷ 1)
            starting_idx = state[2]
            ending_idx = state[2] + next_bptt
            state = (starting_idx, ending_idx)
        else
            state = state .+ d.bptt
        end

        if state[2] >= size(d.dataarray, 2)
            state = (state[1], size(d.dataarray, 2))
        end

        return batch, state
    end
end

nothing

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# 1D
@recipe function f(
    prior::NamedTupleDist,
    vsel::Union{Integer, Symbol, Expr};
    intervals = default_credibilities,
    bins = 200,
    nsamples = 10^6,
    colors = default_colors,
    interval_labels = [],
    closed = :left
)
    orientation = get(plotattributes, :orientation, :vertical)
    (orientation != :vertical) ? swap=true : swap = false
    plotattributes[:orientation] = :vertical # without: auto-scaling of axes not correct

    idx = asindex(prior, vsel)
    if length(idx) > 1
        throw(ArgumentError("Symbol :$vsel refers to a multivariate parameter. Use :($vsel[i]) instead."))
    end

    marg = get_marginal_dist(
        prior,
        idx,
        bins = bins,
        nsamples = nsamples,
        closed = closed
    ).result

    xlabel = if isa(vsel, Symbol) || isa(vsel, Expr)
        "$vsel"
    else
        getstring(prior, idx)
    end
    ylabel = "p("*xlabel*")"

    if swap
        xlabel, ylabel = ylabel, xlabel
    end

    xguide := get(plotattributes, :xguide, xlabel)
    yguide := get(plotattributes, :yguide, ylabel)

   @series begin
        seriestype --> :stephist
        linecolor --> :dimgray
        label --> "prior"

        intervals --> intervals
        bins --> bins
        normalize --> true
        colors --> colors
        interval_labels --> interval_labels

        marg, idx
    end

end

# 2D plots
@recipe function f(
    prior::NamedTupleDist,
    vsel::Union{NTuple{2,Integer}, NTuple{2,Union{Symbol, Expr, Integer}}};
    nsamples=10^6,
    intervals = default_credibilities,
    bins = 200,
    nsamples = 10^6,
    colors = default_colors,
    interval_labels = [],
    closed = :left,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict()
)

    xidx = asindex(prior, vsel[1])
    yidx = asindex(prior, vsel[2])

    if length(xidx) > 1
        throw(ArgumentError("Symbol :$(vsel[1]) refers to a multivariate parameter. Use :($(vsel[1])[i]) instead."))
    elseif length(yidx) > 1
        throw(ArgumentError("Symbol :$(vsel[2]) refers to a multivariate parameter. Use :($(vsel[2])[i]) instead."))
    end

    marg = get_marginal_dist(
        prior,
        (xidx, yidx),
        bins = bins,
        closed = closed
    ).result

    xlabel, ylabel = if isa(vsel, Symbol) || isa(vsel, Expr)
        "$(vsel[1])", "$(vsel[2])"
    else
        getstring(prior, xidx), getstring(prior, yidx)
    end

    xguide := get(plotattributes, :xguide, xlabel)
    yguide := get(plotattributes, :yguide, ylabel)

    @series begin
        seriestype --> :smallest_intervals_contour
        label --> "prior"
        intervals --> intervals
        colors --> colors
        diagonal -->  diagonal
        upper --> upper
        right --> right

        marg, (xidx, yidx)
    end
end


text:
jobs=("Academic librarian",
"Accommodation manager",
"Accountant, chartered",
"Accountant, chartered certified",
"Accountant, chartered management",
"Accountant, chartered public finance",
"Accounting technician",
"Actor",
"Actuary",
"Acupuncturist",
"Administrator",
"Administrator, arts",
"Administrator, charities/voluntary organisations",
"Administrator, Civil Service",
"Administrator, education",
"Administrator, local government",
"Administrator, sports",
"Adult guidance worker",
"Adult nurse",
"Advertising account executive",
"Advertising account planner",
"Advertising art director",
"Advertising copywriter",
"Advice worker",
"Aeronautical engineer",
"Agricultural consultant",
"Agricultural engineer",
"Aid worker",
"Air broker",
"Air cabin crew",
"Air traffic controller",
"Airline pilot",
"Ambulance person",
"Amenity horticulturist",
"Analytical chemist",
"Animal nutritionist",
"Animal technologist",
"Animator",
"Applications developer",
"Arboriculturist",
"Archaeologist",
"Architect",
"Architectural technologist",
"Archivist",
"Armed forces logistics/support/administrative officer",
"Armed forces operational officer",
"Armed forces technical officer",
"Armed forces training and education officer",
"Art gallery manager",
"Art therapist",
"Artist",
"Arts administrator",
"Arts development officer",
"Astronomer",
"Audiological scientist",
"Automotive engineer",
"Banker",
"Barrister",
"Barrister's clerk",
"Best boy",
"Biochemist, clinical",
"Biomedical engineer",
"Biomedical scientist",
"Bonds trader",
"Bookseller",
"Brewing technologist",
"Broadcast engineer",
"Broadcast journalist",
"Broadcast presenter",
"Building control surveyor",
"Building services engineer",
"Building surveyor",
"Buyer, industrial",
"Buyer, retail",
"Cabin crew",
"Call centre manager",
"Camera operator",
"Careers adviser",
"Careers information officer",
"Cartographer",
"Catering manager",
"Ceramics designer",
"Charity fundraiser",
"Charity officer",
"Chartered accountant",
"Chartered certified accountant",
"Chartered legal executive (England and Wales)",
"Chartered loss adjuster",
"Chartered management accountant",
"Chartered public finance accountant",
"Chemical engineer",
"Chemist, analytical",
"Child psychotherapist",
"Chiropodist",
"Chiropractor",
"Civil engineer, consulting",
"Civil engineer, contracting",
"Civil Service administrator",
"Civil Service fast streamer",
"Claims inspector/assessor",
"Clinical biochemist",
"Clinical cytogeneticist",
"Clinical embryologist",
"Clinical molecular geneticist",
"Clinical psychologist",
"Clinical research associate",
"Clinical scientist, histocompatibility and immunogenetics",
"Clothing/textile technologist",
"Colour technologist",
"Commercial art gallery manager",
"Commercial horticulturist",
"Commercial/residential surveyor",
"Commissioning editor",
"Communications engineer",
"Community arts worker",
"Community development worker",
"Community education officer",
"Community pharmacist",
"Company secretary",
"Computer games developer",
"Conference centre manager",
"Conservation officer, historic buildings",
"Conservation officer, nature",
"Conservator, furniture",
"Conservator, museum/gallery",
"Consulting civil engineer",
"Contracting civil engineer",
"Control and instrumentation engineer",
"Copy",
"Copywriter, advertising",
"Corporate investment banker",
"Corporate treasurer",
"Counselling psychologist",
"Counsellor",
"Curator",
"Customer service manager",
"Cytogeneticist",
"Dance movement psychotherapist",
"Dancer",
"Data processing manager",
"Database administrator",
"Dealer",
"Dentist",
"Designer, blown glass/stained glass",
"Designer, ceramics/pottery",
"Designer, exhibition/display",
"Designer, fashion/clothing",
"Designer, furniture",
"Designer, graphic",
"Designer, industrial/product",
"Designer, interior/spatial",
"Designer, jewellery",
"Designer, multimedia",
"Designer, television/film set",
"Designer, textile",
"Development worker, community",
"Development worker, international aid",
"Diagnostic radiographer",
"Dietitian",
"Diplomatic Services operational officer",
"Dispensing optician",
"Doctor, general practice",
"Doctor, hospital",
"Dramatherapist",
"Drilling engineer",
"Early years teacher",
"Ecologist",
"Economist",
"Editor, commissioning",
"Editor, film/video",
"Editor, magazine features",
"Editorial assistant",
"Education administrator",
"Education officer, community",
"Education officer, environmental",
"Education officer, museum",
"Educational psychologist",
"Electrical engineer",
"Electronics engineer",
"Embryologist, clinical",
"Emergency planning/management officer",
"Energy engineer",
"Energy manager",
"Engineer, aeronautical",
"Engineer, agricultural",
"Engineer, automotive",
"Engineer, biomedical",
"Engineer, broadcasting (operations)",
"Engineer, building services",
"Engineer, chemical",
"Engineer, civil (consulting)",
"Engineer, civil (contracting)",
"Engineer, communications",
"Engineer, control and instrumentation",
"Engineer, drilling",
"Engineer, electrical",
"Engineer, electronics",
"Engineer, energy",
"Engineer, land",
"Engineer, maintenance",
"Engineer, maintenance (IT)",
"Engineer, manufacturing",
"Engineer, manufacturing systems",
"Engineer, materials",
"Engineer, mining",
"Engineer, petroleum",
"Engineer, production",
"Engineer, site",
"Engineer, structural",
"Engineer, technical sales",
"Engineer, water",
"Engineering geologist",
"English as a foreign language teacher",
"English as a second language teacher",
"Environmental consultant",
"Environmental education officer",
"Environmental health practitioner",
"Environmental manager",
"Equality and diversity officer",
"Equities trader",
"Ergonomist",
"Estate agent",
"Estate manager/land agent",
"Event organiser",
"Exercise physiologist",
"Exhibition designer",
"Exhibitions officer, museum/gallery",
"Facilities manager",
"Farm manager",
"Fashion designer",
"Fast food restaurant manager",
"Field seismologist",
"Field trials officer",
"Film/video editor",
"Financial adviser",
"Financial controller",
"Financial manager",
"Financial planner",
"Financial risk analyst",
"Financial trader",
"Fine artist",
"Firefighter",
"Fish farm manager",
"Fisheries officer",
"Fitness centre manager",
"Food technologist",
"Forensic psychologist",
"Forensic scientist",
"Forest/woodland manager",
"Freight forwarder",
"Furniture conservator/restorer",
"Furniture designer",
"Further education lecturer",
"Futures trader",
"Gaffer",
"Games developer",
"Garment/textile technologist",
"General practice doctor",
"Geneticist, molecular",
"Geochemist",
"Geographical information systems officer",
"Geologist, engineering",
"Geologist, wellsite",
"Geophysical data processor",
"Geophysicist/field seismologist",
"Geoscientist",
"Glass blower/designer",
"Government social research officer",
"Graphic designer",
"Haematologist",
"Health and safety adviser",
"Health and safety inspector",
"Health physicist",
"Health promotion specialist",
"Health service manager",
"Health visitor",
"Herbalist",
"Heritage manager",
"Higher education careers adviser",
"Higher education lecturer",
"Historic buildings inspector/conservation officer",
"Holiday representative",
"Homeopath",
"Horticultural consultant",
"Horticultural therapist",
"Horticulturist, amenity",
"Horticulturist, commercial",
"Hospital doctor",
"Hospital pharmacist",
"Hotel manager",
"Housing manager/officer",
"Human resources officer",
"Hydrogeologist",
"Hydrographic surveyor",
"Hydrologist",
"Illustrator",
"Immigration officer",
"Immunologist",
"Industrial buyer",
"Industrial/product designer",
"Information officer",
"Information systems manager",
"Insurance account manager",
"Insurance broker",
"Insurance claims handler",
"Insurance risk surveyor",
"Insurance underwriter",
"Intelligence analyst",
"Interior and spatial designer",
"International aid/development worker",
"Interpreter",
"Investment analyst",
"Investment banker, corporate",
"Investment banker, operational",
"IT consultant",
"IT sales professional",
"IT technical support officer",
"IT trainer",
"Jewellery designer",
"Journalist, broadcasting",
"Journalist, magazine",
"Journalist, newspaper",
"Land",
"Land/geomatics surveyor",
"Landscape architect",
"Lawyer",
"Lawyer",
"Learning disability nurse",
"Learning mentor",
"Lecturer, further education",
"Lecturer, higher education",
"Legal executive",
"Legal secretary",
"Leisure centre manager",
"Lexicographer",
"Librarian, academic",
"Librarian, public",
"Licensed conveyancer",
"Lighting technician, broadcasting/film/video",
"Lobbyist",
"Local government officer",
"Location manager",
"Logistics and distribution manager",
"Loss adjuster, chartered",
"Magazine features editor",
"Magazine journalist",
"Maintenance engineer",
"Make",
"Management consultant",
"Manufacturing engineer",
"Manufacturing systems engineer",
"Marine scientist",
"Market researcher",
"Marketing executive",
"Materials engineer",
"Mechanical engineer",
"Media buyer",
"Media planner",
"Medical illustrator",
"Medical laboratory scientific officer",
"Medical physicist",
"Medical sales representative",
"Medical secretary",
"Medical technical officer",
"Mental health nurse",
"Merchandiser, retail",
"Merchant navy officer",
"Metallurgist",
"Meteorologist",
"Microbiologist",
"Midwife",
"Minerals surveyor",
"Mining engineer",
"Mudlogger",
"Multimedia programmer",
"Multimedia specialist",
"Museum education officer",
"Museum/gallery conservator",
"Museum/gallery curator",
"Museum/gallery exhibitions officer",
"Music therapist",
"Music tutor",
"Musician",
"Nature conservation officer",
"Naval architect",
"Network engineer",
"Newspaper journalist",
"Nurse, adult",
"Nurse, children's",
"Nurse, learning disability",
"Nurse, mental health",
"Nutritional therapist",
"Occupational hygienist",
"Occupational psychologist",
"Occupational therapist",
"Oceanographer",
"Office manager",
"Oncologist",
"Operational investment banker",
"Operational researcher",
"Operations geologist",
"Ophthalmologist",
"Optician, dispensing",
"Optometrist",
"Orthoptist",
"Osteopath",
"Outdoor activities/education manager",
"Paediatric nurse",
"Paramedic",
"Passenger transport manager",
"Patent attorney",
"Patent examiner",
"Pathologist",
"Pension scheme manager",
"Pensions consultant",
"Personal assistant",
"Personnel officer",
"Petroleum engineer",
"Pharmacist, community",
"Pharmacist, hospital",
"Pharmacologist",
"Photographer",
"Physicist, medical",
"Physiological scientist",
"Physiotherapist",
"Phytotherapist",
"Pilot, airline",
"Planning and development surveyor",
"Plant breeder/geneticist",
"Podiatrist",
"Police officer",
"Politician's assistant",
"Presenter, broadcasting",
"Press photographer",
"Press sub",
"Primary school teacher",
"Print production planner",
"Printmaker",
"Prison officer",
"Private music teacher",
"Probation officer",
"Producer, radio",
"Producer, television/film/video",
"Product designer",
"Product/process development scientist",
"Production assistant, radio",
"Production assistant, television",
"Production designer, theatre/television/film",
"Production engineer",
"Production manager",
"Programme researcher, broadcasting/film/video",
"Programmer, applications",
"Programmer, multimedia",
"Programmer, systems",
"Proofreader",
"Psychiatric nurse",
"Psychiatrist",
"Psychologist, clinical",
"Psychologist, counselling",
"Psychologist, educational",
"Psychologist, forensic",
"Psychologist, occupational",
"Psychologist, prison and probation services",
"Psychologist, sport and exercise",
"Psychotherapist",
"Psychotherapist, child",
"Psychotherapist, dance movement",
"Public affairs consultant",
"Public house manager",
"Public librarian",
"Public relations account executive",
"Public relations officer",
"Publishing copy",
"Publishing rights manager",
"Purchasing manager",
"Quality manager",
"Quantity surveyor",
"Quarry manager",
"Race relations officer",
"Radiation protection practitioner",
"Radio broadcast assistant",
"Radio producer",
"Radiographer, diagnostic",
"Radiographer, therapeutic",
"Ranger/warden",
"Records manager",
"Recruitment consultant",
"Recycling officer",
"Regulatory affairs officer",
"Research officer, government",
"Research officer, political party",
"Research officer, trade union",
"Research scientist (life sciences)",
"Research scientist (maths)",
"Research scientist (medical)",
"Research scientist (physical sciences)",
"Restaurant manager",
"Restaurant manager, fast food",
"Retail banker",
"Retail buyer",
"Retail manager",
"Retail merchandiser",
"Risk analyst",
"Risk manager",
"Runner, broadcasting/film/video",
"Rural practice surveyor",
"Sales executive",
"Sales professional, IT",
"Sales promotion account executive",
"Science writer",
"Scientific laboratory technician",
"Scientist, audiological",
"Scientist, biomedical",
"Scientist, clinical (histocompatibility and immunogenetics)",
"Scientist, forensic",
"Scientist, marine",
"Scientist, physiological",
"Scientist, product/process development",
"Scientist, research (life sciences)",
"Scientist, research (maths)",
"Scientist, research (medical)",
"Scientist, research (physical sciences)",
"Scientist, water quality",
"Secondary school teacher",
"Secretary/administrator",
"Secretary, company",
"Seismic interpreter",
"Senior tax professional/tax inspector",
"Set designer",
"Ship broker",
"Site engineer",
"Social research officer, government",
"Social researcher",
"Social worker",
"Software engineer",
"Soil scientist",
"Solicitor",
"Solicitor, Scotland",
"Sound technician, broadcasting/film/video",
"Special educational needs teacher",
"Special effects artist",
"Speech and language therapist",
"Sport and exercise psychologist",
"Sports administrator",
"Sports coach",
"Sports development officer",
"Sports therapist",
"Stage manager",
"Statistician",
"Structural engineer",
"Sub",
"Surgeon",
"Surveyor, building",
"Surveyor, building control",
"Surveyor, commercial/residential",
"Surveyor, hydrographic",
"Surveyor, insurance",
"Surveyor, land/geomatics",
"Surveyor, minerals",
"Surveyor, mining",
"Surveyor, planning and development",
"Surveyor, quantity",
"Surveyor, rural practice",
"Systems analyst",
"Systems developer",
"Tax adviser",
"Tax inspector",
"Teacher, adult education",
"Teacher, early years/pre",
"Teacher, English as a foreign language",
"Teacher, music",
"Teacher, primary school",
"Teacher, secondary school",
"Teacher, special educational needs",
"Teaching laboratory technician",
"Technical author",
"Technical brewer",
"Technical sales engineer",
"TEFL teacher",
"Telecommunications researcher",
"Television camera operator",
"Television floor manager",
"Television production assistant",
"Television/film/video producer",
"Textile designer",
"Theatre director",
"Theatre manager",
"Theatre stage manager",
"Theme park manager",
"Therapeutic radiographer",
"Therapist, art",
"Therapist, drama",
"Therapist, horticultural",
"Therapist, music",
"Therapist, nutritional",
"Therapist, occupational",
"Therapist, speech and language",
"Therapist, sports",
"Tour manager",
"Tourism officer",
"Tourist information centre manager",
"Town planner",
"Toxicologist",
"Trade mark attorney",
"Trade union research officer",
"Trading standards officer",
"Training and development officer",
"Translator",
"Transport planner",
"Travel agency manager",
"Tree surgeon",
"Veterinary surgeon",
"Video editor",
"Visual merchandiser",
"Volunteer coordinator",
"Warden/ranger",
"Warehouse manager",
"Waste management officer",
"Water engineer",
"Water quality scientist",
"Web designer",
"Wellsite geologist",
"Writer",
"Youth worker")

job() = jobs[rand(1:length(jobs))]


text:
@testset "alias sampling" begin
    probs = [0.05, 0.1, 0.2, 0.3, 0.2, 0.1, 0.05]
    J, q = GeometricFlux.alias_setup(probs)
    samples = [GeometricFlux.alias_sample(J, q) for _ in 1:1000]
    print
    @test length(J) == length(q) == length(probs)
    @test max(q...) <= 1.1
    @test min(q...) >= 0.0
    @test max(J...) <= length(probs)
    @test min(J...) >= 0.0
    @test length(samples) == 1000
end


text:
include("centrality.jl")

using Random
using StatsBase

function collect_data(dataset::String, k::Int64, exact::Bool)
    T = read_data_unweighted(dataset, k, exact)
    Tcc, lcc = largest_component(T)
    ntrials = 100
    z_evecs = zeros(Float64, Tcc.dimension, ntrials)
    convergences = zeros(Bool, ntrials)
    Random.seed!(1234)
    for trial = 1:ntrials
        x_init = rand(Tcc.dimension)
        x_init /= sum(x_init)
        evec, conv = Z_evec_dynsys(Tcc, 1e-5, 500, x_init)
        println("trial $trial converged = $conv")
        z_evecs[:, trial] = evec
        convergences[trial] = conv
    end

    save("results/$dataset-$k-random.jld2",
         Dict("z_evecs" => z_evecs,
              "conv"    => convergences))
end

function compute_similarity(dataset::String, k::Int64)
    data = load("results/$dataset-$k-random.jld2")
    z_evecs = data["z_evecs"]
    ntrials = size(z_evecs, 2)
    max_diff = 0.0
    min_cor  = 1.0
    for i = 1:ntrials
        zi = z_evecs[:, i]
        for j = (i + 1):ntrials
            zj = z_evecs[:, j]
            @show i, j, norm(zi - zj, 1)
            max_diff = max(max_diff, norm(zi - zj, 1))
            min_cor  = min(min_cor, corspearman(zi, zj))
        end
    end
    return max_diff, min_cor
end

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
import Base.show
export project_start

const initFile = "init.jl"

"""
## User interface
"""
function project_start(pName  ::  String)
    proj = get_project(pName);
    project_start(proj);
    return proj
end

"""
## Internals
"""
mutable struct Project
    name :: String
    progDir  ::  String
    # additional directories to put on LOAD_PATH. Empty by default.
    dirsOnPath :: Vector{String}
end

function show(proj :: Project)
    println("Project [$(proj.name)]");
    println("\tin: $(proj.progDir)")
end

## Default directories on path
function shared_dirs()
#    outV = ConfigLH.shared_dir();
	return Vector{String}();
end

function project_dir(pYear :: Integer)
    return joinpath(homedir(),  "Documents",  "projects",  "p$pYear")
end

## Development directory, where code should reside for unregistered packages
function dev_dir()
	return ENV["JULIA_PKG_DEVDIR"]
end

"""
Start a project

Keep in mind that this is executed from within a module. This limits changes
to the environment.
"""
function project_start(proj :: Project)
    println("Starting project $(proj.name)");
    println("Dir:  $(proj.progDir)")
    @assert isdir(proj.progDir)  "Not found: $(proj.progDir)"
#    push!(LOAD_PATH, proj.progDir);
    add_to_path!(proj.dirsOnPath);
    cd(proj.progDir);
    initPath = joinpath(proj.progDir, initFile);
    if isfile(initPath)
        include(initPath);
    end
    return nothing
end

function add_to_path!(newDirV :: Vector{String})
	if !isempty(newDirV)
		for newDir in newDirV
			if !in(newDir, LOAD_PATH)
				push!(LOAD_PATH, newDir);
			end
		end
	end
	return nothing
end

"""
get_project

Project Definitions
Newest at the top
"""
function get_project(pName :: AbstractString)
    if pName == "test"
        progDir = shared_dir();
        proj = Project(pName, progDir, shared_dirs())
	elseif pName in ["ConfigLH", "CommonLH", "ModelParams", "EconLH"]
    	progDir = joinpath(dev_dir(), pName)
    elseif pName == "sampleModel"
        # Code to try out modeling code (setting parameters etc)
        progDir = joinpath(ConfigLH.julia_dir(), "SampleModel")
        proj = Project(pName, progDir, shared_dirs())
    elseif pName == "collegeQuality"
        # Model with college quality (Oksana, Tatyana)
        progDir = joinpath(project_dir(2019), "college_stratification", "CollegeStrat");
        proj = Project(pName, progDir, shared_dirs())
    else
        error("Invalid project $pName")
    end
    return proj
end

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


text:
mutable struct Simulation{M,TS,O,C}
    model::M
    timestepping_params::TS
    output_params::O
    clock::C
end

"""
    Simulation(;
            model = nothing,
            timestepping_params = nothing,
            output_params = OutputParams(),
            pickup_output_update_flag = false)

Construct a WAVI.jl Simulation object.

Keyword arguments
=================

    - `model`: (required) an instance of a `Model`` object
    - `timestepping_params`: (required) an instance of a `TimesteppingParams` object, which stores information relating to timestepping
    - `output_params`: an instance of an `OutputParams` object, which stores information relating to outputting of solutions
    - `pickup_output_update_flag`: a flag which specifies whether to update the output_params upon picking up.
"""
function Simulation(;
                    model = nothing,
                    timestepping_params = nothing,
                    output_params = OutputParams(),
                    pickup_output_update_flag = false)

    (timestepping_params !== nothing) || throw(ArgumentError("You must specify a timestepping parameters"))

    #compute number of timesteps per output (should be robust for Inf output frequency)
    output_params = set_n_iter_out!(output_params,timestepping_params.dt, timestepping_params.n_iter_total)

    #initialize the clock
    clock = Clock(n_iter = 0, time = 0.0)

    #set the timestep in model parameters (fudge to allow model to see the timestep in velocity solve)
    model = set_dt_in_model!(model, timestepping_params.dt)

    #build the simulation
    simulation = Simulation(model, timestepping_params, output_params, clock)

    #pickup?
    pickup!(simulation, pickup_output_update_flag)

    return simulation

end

include("run_simulation.jl")

function set_dt_in_model!(model, dt)
    model = @set model.params.dt = dt
    return model
end

function set_n_iter_out!(output_params, dt,n_iter_total)
    output_params.output_freq == Inf ? n_iter_out = (n_iter_total + 1) : n_iter_out = round(Int,output_params.output_freq/dt)
    output_params = @set output_params.n_iter_out = n_iter_out
    return output_params
end

function pickup!(simulation, pickup_output_update_flag)
    #@unpack model, timestepping_params, clock = simulation

    if simulation.timestepping_params.niter0 > 0
        n_iter_string =  lpad(simulation.timestepping_params.niter0, 10, "0"); #filename as a string with 10 digits
        println("detected niter0 > 0 (niter0 = $(simulation.timestepping_params.niter0)). Looking for pickup...")
        try
            sim_load = load(string("PChkpt_",n_iter_string, ".jld2"), "simulation")
            println("Pickup successful")

            simulation.model = sim_load.model
            simulation.clock = sim_load.clock
            simulation.output_params = sim_load.output_params #pointer based outputting system means we have to use same output parameters after pickup
        catch
            Throw(error("Pickup error, terminating run"))
        end
    end
    return simulation
end


text:
using BisectPy
using Test

@testset "BisectPy.jl" begin
    include("bisect_left.jl")
    include("bisect_right.jl")
end


text:
@testset "Migrations tests" begin

  @safetestset "Create table" begin

  end;

end;


text:
# for gate in [:H, :X, :Y, :Z, :T, :S, :Tdag, :Sdag]
#     @eval begin
#         export $gate
#         @primitive $gate
#     end
# end

# """
#     H

# The Hadamard gate.

# # Definition

# ```math
# \\frac{1}{\\sqrt{2}} \\begin{pmatrix}
# 1 & 1\\
# 1 & -1
# \\end{pmatrix}
# ```
# """
# H

# for gate in [:X, :Y, :Z]
#     str = """
#     $gate

# The Pauli $gate gate.
# """
#     @eval @doc $str $gate
# end

# """
#     T

# The T gate.
# """
# T

# """
#     Tdag

# The dagger of T gate.
# """
# Tdag

# """
#     S

# The S gate.
# """
# S

# """
#     Sdag

# The dagger of S gate.
# """
# Sdag

# export shift, phase, Rx, Ry, Rz, rot

# """
#     shift(θ::Real)

# Phase shift gate.

# # Definition

# ```math
# \\begin{pmatrix}
# 1 & 0\\
# 0 & e^(im θ)
# \\end{pmatrix}
# ```
# """
# @primitive shift(θ::Real) = Diagonal([1.0, exp(im * θ)])

# """
#     phase(theta)

# Global phase gate.

# # Definition
# ```math
# exp(iθ) \\mathbf{I}
# ```
# """
# @primitive phase(θ::T) where {T<:Real} = Diagonal(Complex{T}[exp(im * θ), exp(im * θ)])
# @primitive Rx(θ::Real) = [cos(θ / 2) -im*sin(θ / 2); -im*sin(θ / 2) cos(θ / 2)]
# @primitive Ry(θ::Real) = [cos(θ / 2) -sin(θ / 2); sin(θ / 2) cos(θ / 2)]
# @primitive Rz(θ::Real) = Diagonal([-im * sin(θ / 2) + cos(θ / 2), im * sin(θ / 2) + cos(θ / 2)])

# for axis in [:X, :Y, :Z]
#     gate = Symbol(:R, lowercase(string(axis)))
#     str = """
#         $gate(theta::Real)

#     Return a rotation gate on $axis axis.
#     """

#     @eval @doc $str $gate
# end

# # TODO: specialize this on simulators
# """
#     rot(axis, θ::T, m::Int=size(axis, 1)) where {T <: Real}

# General rotation gate, `axis` is the rotation axis, `θ` is the rotation angle. `m` is the size of rotation space, default
# is the size of rotation axis.
# """
# @primitive function rot(axis, θ::T, m::Int = size(axis, 1)) where {T<:Real}
#     return Diagonal(fill(cos(θ / 2), m)) - im * sin(θ / 2) * axis
# end

# # # TODO: this should be an instruction
# # function time_evolution_stub(circ::Circuit{:time_evolution}, register::ArrayReg, locs::Locations)
# #     H, dt = circ.free
# #     expv!(st, -im*dt, H)
# #     return register
# # end

# # function time_evolution_stub(circ::Circuit{:time_evolution}, register::ArrayReg, locs::Locations, ctrl_locs::Locations)
# # end

# # function time_evolution_krylov_stub(circ::Circuit{:time_evolution}, register::ArrayReg{1}, lcs::Locations)
# #     Ks, H, dt = circ.free
# #     st = statevec(register)
# #     arnoldi!(Ks, H, st)
# #     expv!(st, -im*dt, Ks)
# #     return register
# # end

# # function time_evolution_krylov_stub(circ::Circuit{:time_evolution}, register::ArrayReg, lcs::Locations, ctrl_locs::Locations)
# # end

# # function (::PrimitiveCircuit{:time_evolution})(H, dt)
# #     Circuit{:time_evolution}(time_evolution_stub, (H, dt))
# # end

# # function (::PrimitiveCircuit{:time_evolution})(Ks, H, dt)
# #     Circuit{:time_evolution}(time_evolution_krylov_stub, (Ks, H, dt))
# # end

# # const time_evolution = PrimitiveCircuit{:time_evolution}()


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
@cuda_text """
template<typename T>
__device__ void slowAtomicAdd(T* addr, T val) {}

template<>
__device__ void slowAtomicAdd<float> (float *addr, float val) {
  unsigned int *a = (unsigned int *) addr;
  for (;;) {
    unsigned int old = *a;
    unsigned int swapped = atomicCAS(a, old,
      __float_as_int(__int_as_float(old) + val));
    if (old == swapped) break;
  }
}

template<>
__device__ void slowAtomicAdd<double>(double *addr, double val) {
  unsigned long long *a = (unsigned long long *) addr;
  for (;;) {
    unsigned long long old = *a;
    unsigned long long swapped = atomicCAS(a, old,
      __double_as_longlong(__longlong_as_double(old) + val));
    if (old == swapped) break;
  }
}
"""

### reshape

Base.reshape(a::ANode, i::Int...) = broadcast(a, i)

function Base.reshape(a::ANode, i::Tuple{Vararg{Int}})
  @assert length(a) == prod(i)
  ANode(:reshape, (a,), i)
end

do_forward!(n::ANode{:reshape}, out, in) = copy!(vec(out), vec(in))

gradient_node(n::ANode{:reshape}, wrt::ANode, b::ANode) = reshape(b, size(wrt))

### broadcast

function Base.broadcast(f, a::ANode, b::ANode)
  shape = zeros(Int, max(ndims(a), ndims(b)))
  for i in 1:length(shape)
    if ndims(b) < i || size(b, i) == 1
      shape[i] = size(a, i)
    elseif ndims(a) < i || size(a, i) == 1 || size(b, i) == size(a, i)
      shape[i] = size(b, i)
    else
     error("Incompatible dimensions in broadcast $(size(a)) and $(size(b))")
    end
  end

  return f(broadcast_node(a, (shape...)), broadcast_node(b, (shape...)))
end

broadcast_node(n::ANode, sz::ANode) = broadcast_node(n, size(sz))

function broadcast_node(n::ANode, sz::Tuple{Vararg{Int}})
  if size(n) == sz
    return n
  elseif length(n) == 1
    ANode(:fill, (n,), sz)
  elseif (ndims(n) == length(sz)-1 && all(i->size(n,i) == sz[i], 1:ndims(n))) ||
          (ndims(n) == length(sz) && all(i->size(n,i) == sz[i], 1:ndims(n)-1) && size(n, ndims(n)) == 1)
    ANode(:repeat, (n,), sz)
  elseif (ndims(n) == length(sz)-1 && all(i->size(n,i) == sz[i+1], 1:ndims(n))) ||
          (ndims(n) == length(sz) && all(i->size(n,i) == sz[i], 2:ndims(n)) && size(n, 1) == 1)
    ANode(:repeat_inner, (n,), sz)
  else
    i = ANode(broadcast((x,y)->x, reshape(1:length(n), size(n)), Array(Int, sz)); c=true)
    ANode(:getindex, (n,i), sz)
    error()
  end
end

### fill

_fill{T}(n::Integer, a::T, b::T) = my_fill(Csize_t(n), a, b)

function _fill{T}(n::Csize_t, a::Array{T}, b::Array{T})
  @assert length(a) == n && length(b) == 1
  fill!(a, b[1])
end

@cuda_gsl _fill{T<:Real}(n::Csize_t, a::&T, b::&T) "a[i] = b[0]"

do_forward!(n::ANode{:fill}, out, in) = _fill(Csize_t(length(out)), out, in)

gradient_node(n::ANode{:fill}, wrt::ANode, b::ANode) = sum(b)

### repeat

function _repeat{T}(n::Csize_t, m::Csize_t, a::Array{T}, b::Array{T})
  @assert length(a) == n && length(b) == m && size(a, ndims(a))*m == n
  k = 0
  for i in 1:size(a, ndims(a))
    a[k + (1:length(b))] = b
    k += length(b)
  end
end

@cuda_gsl _repeat{T<:Real}(n::Csize_t, m::Csize_t, a::&T, b::&T) "a[i] = b[i%m]"

do_forward!(n::ANode{:repeat}, out, in) = _repeat(Csize_t(length(out)), Csize_t(length(in)), out, in)

gradient_node(n::ANode{:repeat}, wrt::ANode, b::ANode) = sum(b, ndims(b))

### repeat_inner

function _repeat_inner{T}(n::Csize_t, m::Csize_t, a::Array{T}, b::Array{T})
  @assert length(a) == n && length(b) == m && size(a, 1)*m == n
  nr = n÷m
  for i in 1:nr
    a[i + (0:nr:nr*length(b)-1)] = b
  end
end

@cuda_gsl _repeat_inner{T<:Real}(n::Csize_t, m::Csize_t, a::&T, b::&T) "a[i] = b[i*m/n]"

do_forward!(n::ANode{:repeat_inner}, out, in) = _repeat_inner(Csize_t(length(out)), Csize_t(length(in)), out, in)

gradient_node(n::ANode{:repeat_inner}, wrt::ANode, b::ANode) = sum(b, 1)

### getindex

function _getindex{T,N}(n::Csize_t, a::Array{T,N}, b::Array{T}, c::Array{T,N})
  @assert length(a) == length(c) == n
  for i in 1:n a[i] = b[round(Int, c[i])] end
end

@cuda_gsl _getindex{T<:Real}(n::Csize_t, a::&T, b::&T, c::&T) "a[i] = b[((int) c[i])-1]"

function bprop_getindex{T,N}(n::Csize_t, a::Array{T}, b::Array{T,N}, c::Array{T,N})
  @assert length(b) == length(c) == n
  for i in 1:n a[round(Int, c[i])] += b[i] end
end

@cuda_gsl bprop_getindex{T<:Real}(n::Csize_t, a::&T, b::&T, c::&T) "slowAtomicAdd(&a[((int) c[i])-1], b[i])"
# // a[c[i]-1] += b[i]; // XXX no atomicAdd for doubles -- how to do this right?"

function Base.getindex{I<:Integer}(n::ANode, i::AbstractArray{I})
  @assert all(1 .<= i .<= length(n))
  ANode(:getindex, (n, ANode(i, c=true)), size(i))
end

do_forward!(n::ANode{:getindex}, out, in1, in2) = _getindex(Csize_t(length(out)), out, in1, in2)

gradient_node(n::ANode{:getindex}, wrt::ANode, b::ANode) = ANode(:getindex_bp, (b,n.input[2]), size(n.input[1]))
function do_forward!(n::ANode{:getindex_bp}, out, in1, in2)
  fill!(out, zero(eltype(out)))
  bprop_getindex(Csize_t(length(in1)), out, in1, in2)
end

### getindex (UnitRange)

function _geturindex{T,N}(n::Csize_t, a::Array{T,N}, b::Array{T}, start::Cint, stop::Cint)
  @assert n == length(a) == stop-start+1
  for i in 1:n a[i] = b[start+i-1] end
end

@cuda_gsl _geturindex{T<:Real}(n::Csize_t, a::&T, b::&T, start::Cint, stop::Cint) "a[i] = b[i+start]"

function bprop_geturindex{T,N}(n::Csize_t, a::Array{T}, b::Array{T,N}, start::Cint, stop::Cint)
  @assert length(b) == n
  for i in 1:n a[i+start-1] += b[i] end
end

@cuda_gsl bprop_geturindex{T<:Real}(n::Csize_t, a::&T, b::&T, start::Cint, stop::Cint) "slowAtomicAdd(&a[start+i], b[i])"
# // a[c[i]-1] += b[i]; // XXX no atomicAdd for doubles -- how to do this right?"

function Base.getindex(n::ANode, i::UnitRange)
  @assert all(1 .<= i .<= length(n))
  ANode(:geturindex, (n,), size(i), i)
end

do_forward!(n::ANode{:geturindex}, out, in1) = _geturindex(Csize_t(length(out)), out, in1, Cint(first(n.arg)), Cint(last(n.arg)))

gradient_node(n::ANode{:geturindex}, wrt::ANode, b::ANode) = ANode(:geturindex_bp, (b,), size(n.input[1]), n.arg)
function do_forward!(n::ANode{:geturindex_bp}, out, in1)
  fill!(out, zero(eltype(out)))
  bprop_geturindex(Csize_t(length(in1)), out, in1, Cint(first(n.arg)), Cint(last(n.arg)))
end


text:
using Wavelets

wf = waveletfilter("db4")
wl = waveletls("cdf9/7")
N = 1024;
x0 = rand(N,N);
L = nscales(N)  # int(log2(N)-2)
tn = 10


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 032cb34e-13b0-11eb-0c89-f3ee9ef0175f
using LibSerialPort, COBS, PlutoUI, AbstractPlotting, WGLMakie

# ╔═╡ ee342e94-13b0-11eb-2b33-efdac65a578f
import Pkg; Pkg.add(["COBS", "PlutoUI", "AbstractPlotting", "WGLMakie"])

# ╔═╡ 5673eb22-13b4-11eb-1c43-cb693f7a44b6
Pkg.add(name="LibSerialPort", rev="master")

# ╔═╡ 7fad639a-13b1-11eb-1ce9-4157899b768f
port = only(get_port_list())

# ╔═╡ 7fada75e-13b1-11eb-074b-db27807fc606
baudrate = 115200

# ╔═╡ 7fb45376-13b1-11eb-2986-75dbf701a3d3
sp = LibSerialPort.open(port, baudrate)

# ╔═╡ 7ae388ee-13b1-11eb-001b-0bcef4960180
md"""
$(@bind a PlutoUI.Slider(0:255))
"""

# ╔═╡ f582d05a-13b1-11eb-27b3-ad4302a17c6e
encode(sp, a)

# ╔═╡ e88b926a-13b1-11eb-20df-89d133b41a6a
function toint(msg)
    result = zero(UInt32)
    for c in msg
        result <<= 8
        result += c
    end
    result
end

# ╔═╡ 6db02932-13b3-11eb-3338-eb5a54e7d1a5
signal = Node(0.0)

# ╔═╡ 13c79b70-13b4-11eb-1611-15b327e9577e
@async while true
    t = toint(decode(sp))
    signal[] = t == 0 ? 0.0 : 6e6/4t
    sleep(1/30)
end

# ╔═╡ 3fdf6176-13b3-11eb-3d9d-11dee7a78063
begin
	h = timeseries(signal, history = 30, resolution = (500, 200))
	ylims!(h, 0, 1000)
	h
end

# ╔═╡ Cell order:
# ╠═ee342e94-13b0-11eb-2b33-efdac65a578f
# ╠═5673eb22-13b4-11eb-1c43-cb693f7a44b6
# ╠═032cb34e-13b0-11eb-0c89-f3ee9ef0175f
# ╠═7fad639a-13b1-11eb-1ce9-4157899b768f
# ╠═7fada75e-13b1-11eb-074b-db27807fc606
# ╠═7fb45376-13b1-11eb-2986-75dbf701a3d3
# ╠═7ae388ee-13b1-11eb-001b-0bcef4960180
# ╠═f582d05a-13b1-11eb-27b3-ad4302a17c6e
# ╠═e88b926a-13b1-11eb-20df-89d133b41a6a
# ╠═6db02932-13b3-11eb-3338-eb5a54e7d1a5
# ╠═13c79b70-13b4-11eb-1611-15b327e9577e
# ╠═3fdf6176-13b3-11eb-3d9d-11dee7a78063


text:
function IOU(boxA, boxB)
	xA = max(boxA[1], boxB[1])
	yA = max(boxA[2], boxB[2])
	xB = min(boxA[3], boxB[3])
	yB = min(boxA[4], boxB[4])
	interArea = max(0, xB - xA + 1) * max(0, yB - yA + 1)
	boxAArea = (boxA[3] - boxA[1] + 1) * (boxA[4] - boxA[2] + 1)
	boxBArea = (boxB[3] - boxB[1] + 1) * (boxB[4] - boxB[2] + 1)
	iou = interArea / float(boxAArea + boxBArea - interArea)
	return iou
end

const cv = OpenCV
net = cv.dnn.DetectionModel(joinpath(ENV["OPENCV_TEST_DATA_PATH"], "dnn", "opencv_face_detector.pbtxt"),joinpath(ENV["OPENCV_TEST_DATA_PATH"], "dnn", "opencv_face_detector_uint8.pb"))
size0 = Int32(300)

cv.dnn.setPreferableTarget(net, cv.dnn.DNN_TARGET_CPU)
cv.dnn.setInputMean(net, (104, 177, 123))
cv.dnn.setInputScale(net, 1.)
cv.dnn.setInputSize(net, size0, size0)

img = OpenCV.imread(joinpath(test_dir, "cascadeandhog", "images", "mona-lisa.png"))

classIds, confidences, boxes = cv.dnn.detect(net, img, confThreshold=Float32(0.5))

box = (boxes[1].x, boxes[1].y, boxes[1].x+boxes[1].width, boxes[1].y+boxes[1].height)
expected_rect = (185,101,129+185,169+101)

@test IOU(box, expected_rect) > 0.8

print("dnn test passed\n")


text:
using DataStructures: OrderedDict
import DataFrames
using StatsBase: StatsBase, mean
import Random
import BlackBoxOptim

metricfunc(metric::Symbol) = begin
    if metric == :rmse
        (E, O) -> √mean((E .- O).^2)
    elseif metric == :nrmse
        (E, O) -> √mean((E .- O).^2) / mean(O)
    elseif metric == :rmspe
        (E, O) -> √mean(((E .- O) ./ O).^2)
    elseif metric == :mae
        (E, O) -> mean(abs.(E .- O))
    elseif metric == :mape
        (E, O) -> mean(abs.((E .- O) ./ O))
    # Nash-Sutcliffe model efficiency coefficient (NSE)
    elseif metric == :ef
        (E, O) -> 1 - sum((E .- O).^2) / sum((O .- mean(O)).^2)
    # Willmott's refined index of agreement (d_r)
    elseif metric == :dr
        (E, O) -> let a = sum(abs.(E .- O)),
            b = 2sum(abs.(O .- mean(O)))
            a <= b ? 1 - a/b : b/a - 1
        end
    else
        error("unsupported metric: $metric")
    end
end
metricfunc(metric) = metric

#HACK: handle index columns with non-identical, but compatible units
# https://github.com/JuliaData/DataFrames.jl/issues/2486
normalize!(dfs::DataFrame...; on) = begin
    for i in on
        cols = getindex.(dfs, !, i)
        elts = eltype.(cols)
        t = promote_type(elts...)
        for (d, c, e) in zip(dfs, cols, elts)
            e != t && setindex!(d, convert.(t, c), !, i)
        end
    end
end

"""
    calibrate(S, obs; <keyword arguments>) -> Config | OrderedDict

Obtain a set of parameters for the given system `S` that simulates provided observation `obs` closely as possible. A multitude of simulations are conducted with a differing combination of parameter sets specified by the range of possible values and the optimum is selected based on a choice of evaluation metric. Internally, differential evolution algorithm from BlackboxOptim.jl is used.

# Arguments
- `S::Type{<:System}`: type of system to be calibrated.
- `obs::DataFrame`: observatioan data to be used for calibration.

# Keyword Arguments
## Configuration
- `config=()`: a single base configuration for the system (can't be used with `configs`).
- `configs=[]`: multiple base configurations for the system (can't be used with `config`).

## Layout
- `index=nothing`: variables to construct index columns of the output; default falls back to `context.clock.time`.
- `target`: variables to construct non-index columns of the output.

## Calibration
- `parameters`: parameters with a range of boundary values to be calibrated within.
- `metric=nothing`: evaluation metric (`:rmse`, `:nrmse`, `:mae`, `:mape`, `:ef`, `:dr`); default is RMSE.

## Multi-objective
- `weight=nothing`: weights for calibrating multiple targets; default assumes equal weights.
- `pareto=false`: returns a dictionary containing Pareto frontier instead of a single solution satisfying multiple targets.

## Advanced
- `optim=()`: extra options for `BlackBoxOptim.bboptimize`.

Remaining keyword arguments are passed down to `simulate` with regard to running system `S`.

See also: [`simulate`](@ref), [`evaluate`](@ref), [`@config`](@ref)

# Examples
```julia-repl
julia> @system S(Controller) begin
           a => 0 ~ preserve(parameter)
           b(a) ~ accumulate
       end;

julia> obs = DataFrame(time=10u"hr", b=200);

julia> p = calibrate(S, obs; target=:b, parameters=:S => :a => (0, 100), stop=10)
...
Config for 1 system:
  S
    a = 20.0
```
"""
calibrate(S::Type{<:System}, obs::DataFrame; config=(), configs=[], kwargs...) = begin
    if isempty(configs)
        calibrate(S, obs, [config]; kwargs...)
    elseif isempty(config)
        calibrate(S, obs, configs; kwargs...)
    else
        @error "redundant configurations" config configs
    end
end
calibrate(S::Type{<:System}, obs::DataFrame, configs::Vector; index=nothing, target, parameters, metric=nothing, weight=nothing, pareto=false, optim=(), kwargs...) = begin
    #HACK: use copy due to normalize!
    obs = copy(obs)
    P = configure(parameters)
    K = parameterkeys(P)
    I = parseindex(index, S) |> keys |> collect
    T = parsetarget(target, S) |> keys |> collect
    n = length(T)
    multi = n > 1
    isnothing(metric) && (metric = :rmse)
    metric = metricfunc(metric)
    IC = [t for t in zip(getproperty.(Ref(obs), I)...)]
    IV = parseindex(index, S) |> values |> Tuple
    snap(s) = getproperty.(s, IV) .|> value in IC
    NT = DataFrames.make_unique([propertynames(obs)..., T...], makeunique=true)
    T1 = NT[end-n+1:end]
    residual(c) = begin
        est = simulate(S; config=c, index, target, snap, verbose=false, kwargs...)
        isempty(est) && return repeat([Inf], n)
        normalize!(est, obs, on=I)
        df = DataFrames.innerjoin(est, obs, on=I, makeunique=true)
        r = [metric(df[!, e], df[!, o]) for (e, o) in zip(T, T1)]
    end
    config(X) = parameterzip(K, X)
    cost(X) = begin
        c = config(X)
        l = length(configs)
        R = Vector(undef, l)
        Threads.@threads for i in 1:l
            R[i] = residual(configure(configs[i], c))
        end
        A = eachrow(hcat(R...)) .|> Iterators.flatten .|> collect |> deunitfy
        e = eachrow(hcat(A...)) |> sum
        # BlackBoxOptim expects Float64, not even Int
        e = Float64.(e)
        multi ? Tuple(e) : e[1]
    end
    #FIXME: input parameters units are ignored without conversion
    range = map(p -> Float64.(Tuple(deunitfy(p))), parametervalues(P))
    method = if multi
        agg = isnothing(weight) ? mean : let w = StatsBase.weights(weight); f -> mean(f, w) end
        (Method=:borg_moea, FitnessScheme=BlackBoxOptim.ParetoFitnessScheme{n}(aggregator=agg))
    else
        ()
    end
    optim_default = (;
        MaxSteps=5000,
        TraceInterval=10,
        RandomizeRngSeed=false,
    )
    #HACK: always initialize random seed first on our end regardless of RandomizeRngSeed option
    # https://github.com/robertfeldt/BlackBoxOptim.jl/issues/158
    Random.seed!(0)
    r = BlackBoxOptim.bboptimize(cost;
        SearchRange=range,
        method...,
        optim_default...,
        optim...
    )
    if multi && pareto
        pf = BlackBoxOptim.pareto_frontier(r)
        OrderedDict(BlackBoxOptim.fitness.(pf) .=> config.(BlackBoxOptim.params.(pf)))
    else
        config(BlackBoxOptim.best_candidate(r))
    end
end

export calibrate


text:
# ---
# title: 1002. Find Common Characters
# id: problem1002
# author: Indigo
# date: 2021-01-30
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/find-common-characters/description/>
# hidden: true
# ---
#
# Given an array `A` of strings made only from lowercase letters, return a list
# of all characters that show up in all strings within the list **(including
# duplicates)**.  For example, if a character occurs 3 times in all strings but
# not 4 times, you need to include that character three times in the final
# answer.
#
# You may return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["bella","label","roller"]
#     Output: ["e","l","l"]
#
#
# **Example 2:**
#
#
#
#     Input: ["cool","lock","cook"]
#     Output: ["c","o"]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 100`
#   2. `1 <= A[i].length <= 100`
#   3. `A[i][j]` is a lowercase letter
#
#
## @lc code=start
using LeetCode

function common_chars(A::Vector{String})
    cnt = counter(A[1])
    for i in 2:length(A)
        tmp_cnt = counter(A[i])
        for k in union(keys(tmp_cnt), keys(cnt))
            cnt[k] = min(cnt[k], tmp_cnt[k])
        end
    end
    vcat([fill(k, v) for (k, v) in cnt if v > 0]...)
end
## @lc code=end


text:
# # Jordan-Wigner transform

# Load the interface to electronic structure packages.
# As with Qiskit and OpenFermion driver libraries, `ElectronicStructure` does
# a bit of calculation, as well.
using ElectronicStructure

# Load `QuantumOps`, the implementation of Pauli and Fermionc operators
using QuantumOps

using ElectronicStructure: MolecularData
using ElectronicStructurePySCF: PySCF

using ZChop

# ## Specify geometry of molecules

# Define geometries for three molecules: specification of atoms and their positions.
# These are molecular hydrogen, `LiH`, and water.
geoms = (
    Geometry(Atom(:H, (0., 0., 0.)), Atom(:H, (0., 0., 0.7414))),

    Geometry(Atom(:Li, (0., 0., 0.)), Atom(:H, (0., 0., 1.4))),

    Geometry(Atom(:O, (0., 0., 0.)), Atom(:H, (0.757, 0.586, 0.)),
             Atom(:H, (-0.757, 0.586, 0.)))
);

# Choose one of the geometries, 1, 2, or 3.
geom = geoms[1]

# Choose a orbital basis set.
basis = "sto-3g";
#basis = "631g"

# We have chosen the crudest basis set and the smallest molecule. Otherwise
# the size of the data structures is too large to display in a demonstration.

# ## Compute interaction integrals for the chosen molecule and basis set

# Construct the specification of the electronic structure problem.
mol_spec = MolecularSpec(geometry=geom, basis=basis)

# Do calculations using `PySCF` and populate a `MolecularData` object with results.
# `mol_pyscf` holds a constant, a rank-two tensor, and a rank-four tensor.
mol_pyscf = MolecularData(PySCF, mol_spec);

# ## Include spin orbitals and change the representation

# Create an interaction operator from one- and two-body integrals and a constant.
# This does just a bit of manipulation of `mol_data`; converting space orbitals
# into space-and-spin orbitals. There are options for choosing chemists' or physicists'
# index ordering and block- or inteleaved-spin orbital ordering.
# We take the default here,
# which gives the same as the operator by the same name in OpenFermion.
# The data is still in the form of rank two and four tensors, but the size of each
# dimension is doubled.
iop = InteractionOperator(mol_pyscf);

zchop!(iop)

# Now we convert `iop` to a more-sparse format; a `FermiSum` (alias for `OpSum{FermiOp}`).
# This is sparse in the sense that only non-zero entries in the tensors in `iop` are represented.
# However, it is not as sparse as it could be, in that identity operators on modes are represented
# explicitly.
fermi_op = FermiSum(iop)

# ## Jordan-Wigner transform

using QuantumOps: jordan_wigner

# Compute the Jordan-Wigner transform
qubit_op = jordan_wigner(fermi_op)

nothing;


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3f7b920e-6c8c-11eb-1bc0-e55f9b4e855f
using PlutoUI

# ╔═╡ 5c755c90-6c8d-11eb-118d-694de45f21c4
using Plots

# ╔═╡ 68d0b3ec-6c8d-11eb-11ed-c7698eadc0c8
md"""
# Dubin's Car Model
"""

# ╔═╡ 75a9037e-6c8d-11eb-0656-0d10f277de0d
md"""
Parameters
"""

# ╔═╡ 7ca2c956-6c8d-11eb-28e0-c58a170fd931
begin
	L = 5 # meters
	v = 5 # m/s
	nothing
end

# ╔═╡ e4ca4704-6c8d-11eb-349f-799dd064b7a0
md"""
Model
"""

# ╔═╡ f6ba63fe-6c8d-11eb-1d4e-67fb69c8ffe4
function next_state(s, ϕ; dt = 0.05)
	# s = [x, θ] = [cte, he]
	# x is in meters and θ is in degrees
	# ϕ is steering angle input (deg)
	x, θ = s

	# Dynamics model
	ẋ = v * sind(θ)
	θ̇ = (v / L) * tand(ϕ)

	x′ = x + ẋ * dt
	θ′ = θ + θ̇ * dt

	return [x′, θ′]
end

# ╔═╡ 0de5c448-6c8f-11eb-2875-e1a10740052b
function next_state(s, y, ϕ; dt = 0.05)
	# for plotting purposes, we might also want the y value
	# s = [x, θ] = [cte, he]
	# x is in meters and θ is in degrees
	# ϕ is steering angle input (deg)
	x, θ = s

	# Dynamics model
	ẋ = v * sind(θ)
	ẏ = v * cosd(θ)
	θ̇ = (v / L) * tand(ϕ)

	x′ = x + ẋ * dt
	y′ = y + ẏ * dt
	θ′ = θ + rad2deg(θ̇) * dt

	return [x′, θ′], y′
end

# ╔═╡ 8e0ccc94-6c8f-11eb-242c-3de6730d1c80
md"""
Simulation
"""

# ╔═╡ 985f921e-6c8f-11eb-334c-8f05ba384d18
function sim_dubins(s₀, y₀, get_steering; num_steps = 500, dt = 0.05, ctrl_every = 20)
	xs = zeros(num_steps + 1)
	θs = zeros(num_steps + 1)
	ys = zeros(num_steps + 1)

	xs[1] = s₀[1]
	θs[1] = s₀[2]
	ys[1] = y₀

	ϕ = clamp(get_steering(s₀), -60, 60)

	for i = 2:num_steps + 1
		s = [xs[i - 1], θs[i - 1]]
		if mod(i - 1, ctrl_every) == 0
			ϕ = clamp(get_steering(s), -60, 60)
		end
		s′, y′ = next_state(s, ys[i - 1], ϕ, dt = dt)
		xs[i] = s′[1]
		θs[i] = s′[2]
		ys[i] = y′
	end

	return xs, θs, ys
end

# ╔═╡ 98ef41dc-6c91-11eb-3a9b-0f44554216f6
md"""
# Plotting
"""

# ╔═╡ a1f56856-6c91-11eb-25ad-6f24fe6a9dbf
function plot_sim(xs, ys)
	p = plot(ys, xs, xlabel = "Downtrack Position", ylabel = "Crosstrack Error", label = "", ylims = (-10, 10))
	plot!(p, ys, zeros(length(xs)), linestyle = :dash, label = "")
	return p
end

# ╔═╡ ad7957b0-6c90-11eb-351c-7d09357eb0c2
md"""
# Control
"""

# ╔═╡ d3d80a00-6c90-11eb-1e5f-21af2406beb8
md"""
k1: $(@bind k1 NumberField(-3:0.01:0, default=0))
k2: $(@bind k2 NumberField(-3:0.01:0, default=0))
"""

# ╔═╡ 7e609550-6c91-11eb-37e9-bb7f6a4b235e
function get_control(s)
	return k1 * s[1] + k2 * s[2]
end

# ╔═╡ e124b84a-6c91-11eb-23fb-6bbe1d53606e
md"""
Initial CTE: $(@bind init_cte NumberField(-10:10, default=0))
"""

# ╔═╡ 93b4b832-6c91-11eb-1dbd-7b7d01523a53
xs, θs, ys = sim_dubins([init_cte, 0.0], 0.0, get_control, ctrl_every = 20,
	num_steps = 1000);

# ╔═╡ 19245a0e-6c92-11eb-32f3-4bc2293c6e3c
plot_sim(xs, ys)

# ╔═╡ Cell order:
# ╠═3f7b920e-6c8c-11eb-1bc0-e55f9b4e855f
# ╠═5c755c90-6c8d-11eb-118d-694de45f21c4
# ╟─68d0b3ec-6c8d-11eb-11ed-c7698eadc0c8
# ╟─75a9037e-6c8d-11eb-0656-0d10f277de0d
# ╠═7ca2c956-6c8d-11eb-28e0-c58a170fd931
# ╟─e4ca4704-6c8d-11eb-349f-799dd064b7a0
# ╠═f6ba63fe-6c8d-11eb-1d4e-67fb69c8ffe4
# ╠═0de5c448-6c8f-11eb-2875-e1a10740052b
# ╟─8e0ccc94-6c8f-11eb-242c-3de6730d1c80
# ╠═985f921e-6c8f-11eb-334c-8f05ba384d18
# ╟─98ef41dc-6c91-11eb-3a9b-0f44554216f6
# ╠═a1f56856-6c91-11eb-25ad-6f24fe6a9dbf
# ╟─ad7957b0-6c90-11eb-351c-7d09357eb0c2
# ╟─d3d80a00-6c90-11eb-1e5f-21af2406beb8
# ╠═7e609550-6c91-11eb-37e9-bb7f6a4b235e
# ╟─e124b84a-6c91-11eb-23fb-6bbe1d53606e
# ╠═93b4b832-6c91-11eb-1dbd-7b7d01523a53
# ╠═19245a0e-6c92-11eb-32f3-4bc2293c6e3c

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
function my_retry(f::Function)
    delays = ExponentialBackOff(; n=10, max_delay=30.0)
    return retry(f; delays=delays)()
end


text:
# @testset "Extras" begin
    ##########
    ## countmap
    ##########

    d = @data [missing,3,3]
    w = weights([1.1,2.2,3.3])
    # cm = Dict{Union{Int,Missing}, Int}([(missing, 1), (3, 2)])
    # cmw = Dict{Union{Int,Missing}, Real}([(missing, 1.1), (3, 5.5)])
    cm = Dict{Union{Missing,Int}, Int}([(missing, 1), (3, 2)])
    cmw = Dict{Union{Missing,Int}, Real}([(missing, 1.1), (3, 5.5)])
    @test isequal(countmap(d), cm)
    @test isequal(countmap(d, w), cmw)

    ##########
    ## cut
    ##########

    @test isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
    @test isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
    @test isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
    @test isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
    @test isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
    @test isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
    @test isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

    ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
    bins = [18, 25, 35, 60, 100]
    cats = cut(ages, bins)
    pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
                           "(25,35]", "(18,25]", "(18,25]",
                           "(35,60]", "(25,35]", "(60,100]",
                           "(35,60]", "(35,60]", "(25,35]"])
    @test isequal(cats, pdv)

    ##########
    ## repeat
    ##########

    @test isequal(repeat(@data [3.0, 2.0, missing]; inner = 2, outer = 1),
                  @data [3.0, 3.0, 2.0, 2.0, missing, missing])
    @test isequal(repeat(@pdata ["a", "b", missing]; inner = 2, outer = 1),
                  @pdata ["a", "a", "b", "b", missing, missing])
    @test isequal(repeat(@data [1 2; 3 missing]; inner = [1, 2], outer = [2, 1]),
                  @data [1 1 2 2; 3 3 missing missing; 1 1 2 2; 3 3 missing missing])
    @test isequal(repeat(@pdata [:a :b missing]; inner = [2,1], outer = [1,3]),
                  @pdata [:a :b missing :a :b missing :a :b missing;
                          :a :b missing :a :b missing :a :b missing])
# end


text:
using ITensors
using OptimKit
using Random
using Zygote

nsites = 20 # Number of sites
nlayers = 3 # Layers of gates in the ansatz
gradtol = 1e-4 # Tolerance for stopping gradient descent

# A layer of the circuit we want to optimize
function layer(nsites, θ⃗)
  RY_layer = [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:nsites]
  CX_layer = [("CX", (n, n + 1)) for n in 1:2:(nsites - 1)]
  return [RY_layer; CX_layer]
end

# The variational circuit we want to optimize
function variational_circuit(nsites, nlayers, θ⃗)
  range = 1:nsites
  circuit = layer(nsites, θ⃗[range])
  for n in 1:(nlayers - 1)
    circuit = [circuit; layer(nsites, θ⃗[range .+ n * nsites])]
  end
  return circuit
end

Random.seed!(1234)

θ⃗ᵗᵃʳᵍᵉᵗ = 2π * rand(nsites * nlayers)
𝒰ᵗᵃʳᵍᵉᵗ = variational_circuit(nsites, nlayers, θ⃗ᵗᵃʳᵍᵉᵗ)

s = siteinds("Qubit", nsites)
Uᵗᵃʳᵍᵉᵗ = ops(𝒰ᵗᵃʳᵍᵉᵗ, s)

ψ0 = MPS(s, "0")

# Create the random target state
ψᵗᵃʳᵍᵉᵗ = apply(Uᵗᵃʳᵍᵉᵗ, ψ0; cutoff=1e-8)

#
# The loss function, a function of the gate parameters
# and implicitly depending on the target state:
#
# loss(θ⃗) = -|⟨θ⃗ᵗᵃʳᵍᵉᵗ|U(θ⃗)|0⟩|² = -|⟨θ⃗ᵗᵃʳᵍᵉᵗ|θ⃗⟩|²
#
function loss(θ⃗)
  nsites = length(ψ0)
  s = siteinds(ψ0)
  𝒰θ⃗ = variational_circuit(nsites, nlayers, θ⃗)
  Uθ⃗ = ops(𝒰θ⃗, s)
  ψθ⃗ = apply(Uθ⃗, ψ0)
  return -abs(inner(ψᵗᵃʳᵍᵉᵗ, ψθ⃗))^2
end

θ⃗₀ = randn!(copy(θ⃗ᵗᵃʳᵍᵉᵗ))

@show loss(θ⃗₀), loss(θ⃗ᵗᵃʳᵍᵉᵗ)

loss_∇loss(x) = (loss(x), convert(Vector, loss'(x)))
algorithm = LBFGS(; gradtol=gradtol, verbosity=2)
θ⃗ₒₚₜ, lossₒₚₜ, ∇lossₒₚₜ, numfg, normgradhistory = optimize(loss_∇loss, θ⃗₀, algorithm)

@show loss(θ⃗ₒₚₜ), loss(θ⃗ᵗᵃʳᵍᵉᵗ)

nothing


text:
@testset "Generate File Template Function Tests" begin

  @test isdefined(Julz, :generate_file_template) == true

end


text:
Base.@kwdef mutable struct Solver <: AMGXObject
    handle::API.AMGX_solver_handle = API.AMGX_solver_handle(C_NULL)
    resources::Union{Resources, Nothing} = nothing
    config::Union{Config, Nothing} = nothing
    mode::Union{Mode, Nothing} = nothing
    bound_matrix::Union{AMGXMatrix, Nothing} = nothing
    function Solver(handle::API.AMGX_solver_handle, resources::Union{Resources, Nothing}, mode::Union{Mode, Nothing},
                    config::Union{Config, Nothing}, bound_matrix::Union{AMGXMatrix, Nothing})
        solver = new(handle, resources, config, mode, bound_matrix)
        finalizer(warn_not_destroyed_on_finalize, solver)
        return solver
    end
end
function dec_refcount_parents(solver::Solver)
    dec_refcount!(solver.resources)
    dec_refcount!(solver.config)
    solver.resources = nothing
    solver.config = nothing
    solver.bound_matrix = nothing
    solver.mode = nothing
    nothing
end
get_api_destroy_call(::Type{Solver}) = API.AMGX_solver_destroy

function create!(solver::Solver, res::Resources, mode::Mode, config::Config)
    solver_handle_ptr = Ref{API.AMGX_solver_handle}()
    amgx_mode = API.AMGX_Mode(Int(mode))
    @checked API.AMGX_solver_create(solver_handle_ptr, res.handle, amgx_mode, config.handle)
    solver.handle = solver_handle_ptr[]
    solver.resources = res
    solver.config = config
    solver.mode = mode
    inc_refcount!(solver.resources)
    inc_refcount!(solver.config)
    return solver
end
Solver(res::Resources, mode::Mode, config::Config) = create!(Solver(), res, mode, config)

function setup!(solver::Solver, matrix::AMGXMatrix)
    @checked API.AMGX_solver_setup(solver.handle, matrix.handle)
    solver.bound_matrix = matrix
    return solver
end

function solve!(sol::AMGXVector, solver::Solver, rhs::AMGXVector; zero_inital_guess::Bool=false)
    if solver.bound_matrix === nothing
        error("no matrix attached to solver")
    elseif solver.bound_matrix.handle == C_NULL
        error("matrix has already been destroyed")
    end
    if zero_inital_guess
        @checked API.AMGX_solver_solve_with_0_initial_guess(solver.handle, rhs.handle, sol.handle)
    else
        @checked API.AMGX_solver_solve(solver.handle, rhs.handle, sol.handle)
    end
    return sol
end

function get_iterations_number(solver::Solver)
    n_ptr = Ref{Cint}()
    @checked API.AMGX_solver_get_iterations_number(solver.handle, n_ptr)
    # for some reason AMGX returns 1 + number of iterations.
    return Int(n_ptr[]) - 1
end

function get_iteration_residual(solver::Solver, iter::Int=get_iterations_number(solver), block_idx::Int=0)
    res_ptr = Ref{Float64}()
    @checked API.AMGX_solver_get_iteration_residual(solver.handle, iter, block_idx, res_ptr)
    return res_ptr[]
end

@enum SolverStatus begin
   SUCCESS = Int(API.AMGX_SOLVE_SUCCESS)
   FAILED = Int(API.AMGX_SOLVE_FAILED)
   DIVERGED = Int(API.AMGX_SOLVE_DIVERGED)
end

function get_status(solver::Solver)
    status_ptr = Ref{API.AMGX_SOLVE_STATUS}()
    @checked API.AMGX_solver_get_status(solver.handle, status_ptr)
    return SolverStatus(Int(status_ptr[]))
end


text:
using DynamicalBilliards, InteractiveDynamics, GLMakie

SV = SVector{2, Float32}
w = 0.4 # stump width
points_right = [
    SV(w, 0),
    SV(w, w),
    SV(2.0, w),
    SV(w, 1.5),
    SV(1.5, 1.5),
    SV(w, 2.25),
    SV(1.0, 2.25),
]
summit = SV(0.0, 3.0)
points_left = reverse!([SV(-p[1], p[2]) for p in points_right])
tree_points = [points_right..., summit, points_left...]

tree_billiard = billiard_vertices(tree_points)

ps = particlebeam(0.0, 0.1, π/4, 1, 0.01, nothing, Float32)

interactive_billiard(tree_billiard, ps; dt = 0.1)


text:
#=
The native lib provides APIs for aggregating public keys and signatures through the
[BLS_ZZZ_ADD_G1 and BLS_ZZZ_ADD_G2](https://github.com/apache/incubator-milagro-crypto-c/blob/2.0.1/include/bls.h.in#L76-L92)
functions. If your application needs it, its perfectly possible to implement the corresponding
BLS_ZZZ_SUB_G1 and BLS_ZZZ_SUB_G2 functions in Julia.
=#

using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
#Pkg.add("https://github.com/blenessy/AMCL.jl")

using MilagroCrypto

function BLS_BLS381_SUB_G1(R1,R2,R)
    P, T = AMCL.ECP_BLS381(undef), AMCL.ECP_BLS381(undef)
    isone(AMCL.ECP_BLS381_fromOctet(P, R1)) || return AMCL.BLS_INVALID_G1
    isone(AMCL.ECP_BLS381_fromOctet(T, R2)) || return AMCL.BLS_INVALID_G1
    AMCL.ECP_BLS381_sub(P, T)
    AMCL.ECP_BLS381_toOctet(R, P, true)
    return AMCL.BLS_OK
end

# get a valid point and serialize it so we have something to test with
P = AMCL.ECP_BLS381(undef)
AMCL.ECP_BLS381_generator(P)
p = octet(AMCL.BFS_BLS381 + 1) # compressed format
AMCL.ECP_BLS381_toOctet(p, P, true)

# p2 = p + p
p2 = octet(AMCL.BFS_BLS381 + 1) # compressed format
AMCL.BLS_BLS381_ADD_G1(p, p, p2) == AMCL.BLS_OK ||
    error("aggregation failed")

# q = p2 - p = p
q = octet(AMCL.BFS_BLS381 + 1) # compressed format
BLS_BLS381_SUB_G1(p2, p, q) == AMCL.BLS_OK ||
   error("failed to subtract")
@assert p == q "BLS_BLS381_SUB_G1 does not work :("

println("Done (nothing bad happened)!")


text:
using RadialBasisFunctionModels
RBF = RadialBasisFunctionModels

sites = [rand(2) for i = 1 : 5]
vals = [rand(2) for i = 1 : 5]

#%%

r = RBFMachine(; features = sites, labels = vals)


text:
using contact_sim

function main()
    q1 = 0.2; u1 = -0.25;
    q2 = -0.4; u2 = 0.2;

    z0 = [q1 u1 q2 u2]
    steps = 4
    fps = 20

    #passive: no controller
    #=
    walker = Walker(1.0, 0.5, 0.02, 1.0, 0.5, 1.0, 0.01, false, 1.0, 100.0, 0.28564, z0[3], z0[4], 0.0)
    z_init = [0.162597833780035  -0.231869638058927  -0.325195667560070   0.03797846807373]
    zpert = z_init
    =#

    #active: use controller
    walker = Walker(1.0, 0.5, 0.02, 1.0, 0.5, 1.0, 0.01, true, 1.0, 100.0, 0.28564, z0[3], z0[4], 0.0)

    #fixed point initial value.
    z_init = [ 0.142819999999995  -0.326813112785275  -0.285640000000000   0.068243824367047]

    #adding slight perturbation
    zpert = z_init + [0. 0.05 -0.1 0.2]

    z, t = simulate_walking(zpert,walker,steps)
    render_walking_trajectory(t,z,walker, steps, fps,true)
end

main()


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2015

using Pkg

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = chomp(readline(open(dir * "/../VERSION")))

start = time_ns()
println( stderr, "Whippet $ver loading and compiling... " )

Pkg.activate(dir * "/..")

using Libz
using ArgParse
using StatsBase
using BioSequences

using Whippet

function parse_cmd()
  s = ArgParseSettings()
  # TODO finish options...
  @add_arg_table s begin
    "--index", "-x"
      help = "Prefix for index 'dir/prefix' (default Whippet/index/graph)"
      arg_type = String
      default  = fixpath( "$(dir)/../../index/graph" )
    "--out", "-o"
      help = "Where should the gzipped output go 'dir/prefix'?"
      arg_type = String
      default  = fixpath( "$(dir)/../../simul" )
    "--gtf"
      help = "Gene anotation file in GTF format"
      arg_type = String
    "--max-complexity", "-k"
      help = "Max complexity KMax allowed per gene (random K for each gene between K1:KMax)?"
      arg_type = Int64
      default  = 8
    "--num-genes", "-n"
      help = "Randomly sample at most -n number of genes from the index to simulate."
      arg_type = Int64
      default  = 5000
  end
  return parse_args(s)
end

function main()

   args = parse_cmd()

   println(stderr, " $( round( (time_ns()-start)/1e9, digits=6 ) ) seconds" )

   println(stderr, "Loading splice graph index... $( args["index"] ).jls")
   @timer const lib = open(deserialize, "$( args["index"] ).jls")

   println(stderr, "Simulating alternative transcripts of $( args["max-complexity"] ) max-complexity..")
   @timer simulate_genes( lib, anno, args["max-complexity"], output=args["out"], gene_num=args["num-genes"] )

end

type SimulTranscript
   seq::Whippet.SGSequence
   nodes::Vector{UInt}
end

immutable SimulGene
   trans::Vector{SimulTranscript}
   gene::String
   complexity::Int
end

istxstart( edge::EdgeType ) = edge == EDGETYPE_SL || edge == EDGETYPE_LS ? true : false
istxstop(  edge::EdgeType ) = edge == EDGETYPE_SR || edge == EDGETYPE_RS ? true : false

function collect_nodes!( st::SimulTranscript, sg::SpliceGraph, r::UnitRange; skip=Vector{Int}(), fromdonor=false, isdependent=false )
   n = r.start
   while n <= r.stop
      if     fromdonor && n in skip
         isdependent = true
      elseif fromdonor && istxstart( sg.edgetype[n] )
         isdependent = true
      elseif fromdonor && istxstop( sg.edgetype[n+1] ) && n < length(sg.nodelen)
         isdependent = true
      elseif isdependent && !( sg.edgetype[n] in (EDGETYPE_RR, EDGETYPE_LR, EDGETYPE_SR) )
      else
         noderange = sg.nodeoffset[n]:(sg.nodeoffset[n]+sg.nodelen[n]-1)
         st.seq *= sg.seq[ noderange ]
         push!( st.nodes, n )
         fromdonor = sg.edgetype[n+1] in (EDGETYPE_LS, EDGETYPE_LR, EDGETYPE_LL) ? true : false
         isdependent = false
      end
      n += 1
   end
end

function simulate_genes( lib, anno, max_comp; output="simul_genes", gene_num=length(lib.graphs) )
   fastaout = open( output * ".fa.gz", "w" )
   nodesout = open( output * ".node.gz", "w" )
   fastastr = ZlibDeflateOutputStream( fastaout )
   nodesstr = ZlibDeflateOutputStream( nodesout )

   for g in sample( 1:length(lib.graphs), min( length(lib.graphs), gene_num ), replace=false )
      comp = min( max_comp, length(lib.graphs[g].nodelen) - 2 )
      sgene = SimulGene( Vector{SimulTranscript}(), lib.names[g], rand(1:comp) )
      simulate_transcripts!( sgene, lib.graphs[g] )
      output_transcripts( fastastr, sgene, lib.graphs[g] )
      output_nodes( nodesstr, lib.names[g], anno.geneset[lib.names[g]].info, lib.graphs[g] )
   end
   close( fastastr )
   close( fastaout )
   close( nodesstr )
   close( nodesout )
end

function simulate_transcripts!( simul::SimulGene, sg::SpliceGraph )
   trans = SimulTranscript( dna"", Vector{UInt}() )
   collect_nodes!( trans, sg, 1:length(sg.nodelen) )
   push!( simul.trans, trans )
   if length(sg.nodelen) > 2
      mod_start = rand( 1:(length(sg.nodelen) - simul.complexity) )
      for i in 1:simul.complexity
         combset = collect( combinations( mod_start:(mod_start+simul.complexity-1), i ) )
         for s in combset
            trans = SimulTranscript( dna"", Vector{UInt}() )
            collect_nodes!( trans, sg, 1:length(sg.nodelen), skip=s )
            push!( simul.trans, trans )
         end
      end
   end
end

function output_transcripts( stream, simul::SimulGene, sg::SpliceGraph )
   prefix = simul.gene * "_K" * string(simul.complexity) * "_"
   used = Set()
   for t in simul.trans
      name = prefix * join( t.nodes, "-" )
      (name in used) && continue
      write( stream, ">$name\n" )
      for i in t.seq
         write( stream, Char(i) )
      end
      write( stream, "\n" )
      push!( used, name )
   end
end

function output_nodes( stream, gname::String, info, sg::SpliceGraph)
   for n in 1:length(sg.nodelen)
      coord = "$(info[1]):$(sg.nodecoord[n])-$(sg.nodecoord[n] + sg.nodelen[n] - 1)"
      write( stream, "$gname\t$n\t$coord\t$(info[2])\n" )
   end
end

main()


text:
"""
    mapmany(f, a...)

Like `map`, but `f(x)` for each `x ∈ a` may return an arbitrary number of values to insert
into the output.

# Example

```jldoctest
julia> mapmany(x -> 1:x, [1,2,3])
6-element Array{Int64,1}:
 1
 1
 2
 1
 2
 3
```
"""
function mapmany(f, a)
    T = eltype(promote_op(f, eltype(a)))
    out = T[]
    for x ∈ a
        append!(out, f(x))
    end
    return out
end

mapmany(f, a, b, c...) = mapmany(x->f(x...), zip(a, b, c...))

"""
    flatten(a)

Takes a collection of collections `a` and returns a collection containing all the elements
of the subcollecitons of `a`. Equivalent to `mapmany(idenity, a)`.

# Example

```jldoctest
julia> flatten([1:1, 1:2, 1:3])
6-element Array{Int64,1}:
 1
 1
 2
 1
 2
 3
```
"""
flatten(x) = mapmany(identity, x)

struct MappedIterator{F, T}
	f::F
	parent::T
end
@inline function iterate(it::MappedIterator)
	x = iterate(it.parent)
	return x === nothing ? nothing : (it.f(x[1]), x[2])
end
@inline function iterate(it::MappedIterator, i)
	x = iterate(it.parent, i)
	return x === nothing ? nothing : (it.f(x[1]), x[2])
end

# Inherit any iteration and indexing properties from parent
Base.IteratorSize(it::MappedIterator) = Base.IteratorSize(it.parent)
Base.length(it::MappedIterator) = length(it.parent)
Base.size(it::MappedIterator) = size(it.parent)
Base.axes(it::MappedIterator) = axes(it.parent)
Base.keys(it::MappedIterator) = keys(it.parent)
@propagate_inbounds Base.getindex(it::MappedIterator, i) = it.parent[i]

struct MappedArray{T, N, F, A <: AbstractArray{<:Any, N}} <: AbstractArray{T, N}
	f::F
	parent::A
end
parent(a::MappedArray) = a.parent
size(a::MappedArray) = size(a.parent)
axes(a::MappedArray) = axes(a.parent)
IndexStyle(a::MappedArray) = IndexStyle(a.parent)
@propagate_inbounds getindex(a::MappedArray{T}, i::Int) where {T} = a.f(a.parent[i])::T
@propagate_inbounds getindex(a::MappedArray{T}, i::Int...) where {T} = a.f(a.parent[i...])::T

"""
    mapview(f, a)

Return a view of `a` where each element is mapped through function `f`. The iteration and
indexing properties of `a` are preserved. Similar to `map(f, a)`, except evaluated lazily.

# Example

```julia
julia> a = [1,2,3];

julia> b = mapview(-, a)
3-element MappedArray{Int64,1,typeof(-),Array{Int64,1}}:
 -1
 -2
 -3

julia> a[1] = 10;

julia> b
3-element MappedArray{Int64,1,typeof(-),Array{Int64,1}}:
 -10
  -2
  -3
```
"""
mapview(f, a) = MappedIterator(f, a)
mapview(f, a::AbstractArray{T, N}) where {T, N} = MappedArray{promote_op(f, T), N, typeof(f), typeof(a)}(f, a)

mapview(::typeof(identity), a) = a
mapview(::typeof(identity), a::AbstractArray) = a


text:
using MPFirst
using Test

@testset "MPFirst.jl" begin
    # Write your tests here.
end


text:
abstract type Sequence end
value(d::Sequence) = d.val
step!(d::Sequence) = d.val

struct ConstantSequence{T} <: Sequence
    val::T
end

mutable struct ExponentialSequence{T} <: Sequence
    val::T
    multiplier::T
end
step!(ED::ExponentialSequence) = ED.val *= ED.multiplier

mutable struct HeavisideSequence{T} <: Sequence
    origVal::T  # Value to start with
    lateVal::T  # Value to switch to
    switchIter::Int
    iter::Int
end
HeavisideSequence(origVal, lateVal, switchIter, iter = 1) = HeavisideSequence(origVal, lateVal, switchIter, iter)
value(HD::HeavisideSequence) = HD.iter < HD.switchIter ? HD.origVal : HD.lateVal
function step!(HD::HeavisideSequence)
    HD.iter += 1
    HD.iter < HD.switchIter ? HD.origVal : HD.lateVal
end

mutable struct ArithmeticSequence{T} <: Sequence
    val::T
    addVal::T
end
step!(AD::ArithmeticSequence) = AD.val += AD.addVal

mutable struct LinearSequence{T} <: Sequence
    startVal::T
    endVal::T
    startIter::Int
    endIter::Int
    currIter::Int
    step::T
    val::T
end
function LinearSequence(startVal::T, endVal::T, endIter::Int; startIter=1, iter = 1) where {T}
    step = (endVal - startVal)/(endIter - startIter)
    val = iter < endIter ? startVal + step*max(iter - startIter, 0) : endVal
    LinearSequence(startVal, endVal, startIter, endIter, 1, step, val)
end

function step!(LD::LinearSequence)
    LD.currIter += 1
    LD.val = LD.currIter < LD.endIter ? LD.startVal + LD.step*max(LD.currIter - LD.startIter, 0) : LD.endVal
end


text:
include("test_2d.jl")
include("test_3d.jl")
include("test_manual.jl")


text:
abstract type AbstractSystem end
struct UnknownSystem <: AbstractSystem end

struct Linux <: AbstractSystem
    id::String
    version_id::String
    os_name::String
    os_version::String
end

struct MacOS <: AbstractSystem
    version_id::String
end

"""
    current_system() -> AbstractSystem

Query the current system information.
"""
function current_system()
    CURRENT_SYS[1] !== UnknownSystem() && return CURRENT_SYS[1]

    if Sys.islinux()
        info = map(readlines("/etc/os-release")) do line
            k, v = split(line, "=")
            Symbol(lowercase(k)) => v
        end
        info = (; info...)
        CURRENT_SYS[1] = Linux(info.id, info.version_id, info.name, info.version)
    elseif Sys.isapple()
        verson_id = on_success(strip, `sw_vers -productVersion`)
        CURRENT_SYS[1] = MacOS(verson_id)
    else
        error("Current system is not supported yet.")
    end
    return CURRENT_SYS[1]
end


text:
module TestUnitsAndUncertainites

using TinyModia
using DifferentialEquations
using ModiaPlot
using Measurements
using Unitful

TestModel = Model(
    T = (0.2 ± 0.02)u"s",
    u0 = 8.5u"kg" ± 0.85u"kg",
    k = (1 ± 0.1)u"kg^-1",
    m = (10 ± 1)u"kg",
    F = (20 ± 2)u"N",

    x = Var(init = (5.5 ± 0.55)u"kg"),
    v = Var(init = (1 ± 0.1)u"m/s"),
    r = Var(init = 1u"m/s"),
    equations = :[
    T*der(x) + x = u0
    y = x
    z = x*k

    F = m * a
    der(v) = a
    der(r) = v ]
)

model = @instantiateModel(TestModel, FloatType = Measurement{Float64})

simulate!(model, Tsit5(), stopTime = 1.0)

plot(model, ["T", "x", "der(x)", "y", "a", "v", "r"])

end


text:
module GalacticNOMAD

using NOMAD, GalacticOptim, GalacticOptim.SciMLBase

export NOMADOpt
struct NOMADOpt end

function __map_optimizer_args!(prob::OptimizationProblem, opt::NOMAD.NomadProblem;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    for j in kwargs
        setproperty!(opt.options, j.first, j.second)
    end

    if !isnothing(maxiters)
        opt.options.max_bb_eval = maxiters
    end

    if !isnothing(maxtime)
        opt.options.max_time = maxtime
    end

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    if !isnothing(abstol)
        @warn "common abstol is currently not used by $(opt)"
    end

    return nothing
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::NOMADOpt;
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    kwargs...)

    local x

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    _loss = function (θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    function bb(x)
        l = _loss(x)
        success = !isnan(l) && !isinf(l)
        count_eval = true
        return (success, count_eval, [l])
    end

    if !isnothing(prob.lcons) | !isnothing(prob.ucons)
        @warn "Linear and nonlinear constraints defined in OptimizationProblem are currently not used by $(opt)"
    end

    bounds = (;)
    if !isnothing(prob.lb)
        bounds = (; bounds..., lower_bound=prob.lb)
    end

    if !isnothing(prob.ub)
        bounds = (; bounds..., upper_bound=prob.ub)
    end

    opt_setup = NOMAD.NomadProblem(length(prob.u0), 1, ["OBJ"], bb; bounds...)

    __map_optimizer_args!(prob, opt_setup, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    opt_res = NOMAD.solve(opt_setup, prob.u0)
    t1 = time()

    SciMLBase.build_solution(prob, opt, opt_res.x_best_feas, first(opt_res.bbo_best_feas); original=opt_res)
end

end


text:
function BrownianBridge(t0,tend,W0,Wend,Z0=nothing,Zend=nothing;kwargs...)
  W = WienerProcess(t0,W0,Z0;kwargs...)
  h = tend-t0
  Wh = Wend-W0
  if Z0 != nothing
    Zh = Zend - Z0
  else
    Zh = nothing
  end
  push!(W.S₁,(h,Wh,Zh))
  W
end

function BrownianBridge!(t0,tend,W0,Wh,Z0=nothing,Zh=nothing;kwargs...)
  W = WienerProcess!(t0,W0,Z0;kwargs...)
  h = tend-t0
  Wh .-= W0
  if Z0 != nothing
    Zh .-= Z0
  else
    Zh = nothing
  end
  push!(W.S₁,(h,Wh,Zh))
  W
end

function GeometricBrownianBridge(μ,σ,t0,tend,W0,Wend,Z0=nothing,Zend=nothing;kwargs...)
  W = GeometricBrownianMotionProcess(μ,σ,t0,W0,Z0;kwargs...)
  h = tend-t0
  Wh = Wend-W0
  if Z0 != nothing
    Zh = Zend - Z0
  else
    Zh = nothing
  end
  push!(W.S₁,(h,Wh,Zh))
  W
end

function GeometricBrownianBridge!(μ,σ,t0,tend,W0,Wh,Z0=nothing,Zh=nothing;kwargs...)
  W = GeometricBrownianMotionProcess!(μ,σ,t0,W0,Z0;kwargs...)
  h = tend-t0
  Wh .-= W0
  if Z0 != nothing
    Zh .-= Z0
  else
    Zh = nothing
  end
  push!(W.S₁,(h,Wh,Zh))
  W
end


text:
## --- Custom objects for holding Chron age data

    # Define a type of struct object to hold data about geochron samples from
    # a region or stratigraphic section of interest
    mutable struct ChronAgeData
        Name::Tuple
        Height::Array{Float64,1}
        Height_sigma::Array{Float64,1}
        Age::Array{Float64,1}
        Age_sigma::Array{Float64,1}
        Age_025CI::Array{Float64,1}
        Age_975CI::Array{Float64,1}
        Age_14C::Array{Float64,1}
        Age_14C_sigma::Array{Float64,1}
        Age_Sidedness::Array{Float64,1}
        Age_DistType::Array{Float64,1}
        Age_Distribution::Array{<:Array,1}
        Params::Array{Float64,2}
        Path::String
        inputSigmaLevel::Int
        Age_Unit::String
        Height_Unit::String
    end

    function NewChronAgeData(nSamples)
        smpl = ChronAgeData(
            ("Sample Names",),
            fill(NaN,nSamples),  # Sample heights
            fill(NaN,nSamples),  # Height_sigma
            fill(NaN,nSamples),  # Sample ages
            fill(NaN,nSamples),  # Sample age uncertainty
            fill(NaN,nSamples),  # Sample age 2.5% CI
            fill(NaN,nSamples),  # Sample age 97.5% CI
            fill(NaN,nSamples),  # Sample 14C ages
            fill(NaN,nSamples),  # Sample 14C uncertainties
            zeros(nSamples), # Sidedness (zeros by default, geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
            zeros(nSamples), # DistType (for Distribution-fitting only: 0=distribution to be fit, 1=single Gaussian)
            Array{Array}(undef,nSamples), # Stationary distribution of eruption age
            fill(NaN,5,nSamples), # Sample age distribution parameters
            "./", # Relative path where we can find .csv data files
            2, # i.e., are the data files 1-sigma or 2-sigma
            "Ma",
            "m",
        )
        return smpl
    end
    export NewChronAgeData

    # Just for backwards compatibility (hopefully not neccessary)
    NewStratAgeData = NewChronAgeData
    export NewStratAgeData

    # A type of object to hold data about hiatuses
    mutable struct HiatusData
        Height::Array{Float64,1}
        Height_sigma::Array{Float64,1}
        Duration::Array{Float64,1}
        Duration_sigma::Array{Float64,1}
        Age_Unit::String
        Height_Unit::String
    end

    function NewHiatusData(nHiatuses)
        hiatus = HiatusData(
            fill(NaN,nHiatuses),  # Height
            fill(NaN,nHiatuses),  # Height_sigma
            fill(NaN,nHiatuses),  # Duration
            fill(NaN,nHiatuses),  # Duration_sigma
            "Ma",
            "m",
        )
        return hiatus
    end
    export NewHiatusData

    # A type of object to specify the configuration of the stratigraphic model
    mutable struct StratAgeModelConfiguration
        resolution::Float64
        burnin::Int
        nsteps::Int
        sieve::Int
        bounding::Float64
    end

    function NewStratAgeModelConfiguration()
        return StratAgeModelConfiguration(NaN, 0, 0, 0, NaN)
    end
    export NewStratAgeModelConfiguration

    struct StratAgeModel
        Height::Array{Float64,1}
        Age::Array{Float64,1}
        Age_sigma::Array{Float64,1}
        Age_Median::Array{Float64,1}
        Age_025CI::Array{Float64,1}
        Age_975CI::Array{Float64,1}
    end